最小二乗問題：正規方程式・QR・特異値分解と感度

0. この記事の要点

測った点の数が未知数より多く、すべての式を同時には満たせないとき、「いちばんずれが小さい答え」をどう選び、どう正確に計算するかを扱います。たとえば 3 つの測定点に 1 本の直線を引くと、普通は 3 点すべてを通りません。そこで各点とのずれの 2 乗和を最小にする直線を選びます。これが最小二乗法です。この記事では、その答えの性質、計算法ごとの精度の違い、データに雑音があるときの対処を順に見ます。

• $\|A\boldsymbol{x}-\boldsymbol{b}\|_2$ を最小にする $\boldsymbol{x}$ は、残差が $A$ の像に直交する $\boldsymbol{x}$、すなわち正規方程式 $A^{*}A\boldsymbol{x}=A^{*}\boldsymbol{b}$（$A^{*}$ は共役転置）の解です。列が一次独立なら一意です。
• 特異値分解から作る擬似逆 $A^{+}$ について、$A^{+}\boldsymbol{b}$ は最小二乗解のうちノルムが最小のものです。$A^{+}$ は Moore–Penrose の 4 条件（$AXA=A$ など 4 つの等式。第 4 節）で一意に決まります。
• QR 分解 $A=Q_1R_1$ を使えば $\boldsymbol{x}=R_1^{-1}Q_1^{\mathsf{T}}\boldsymbol{b}$ で、Householder QR による解法は後退安定です（後退安定とは、計算結果が「入力をわずかに動かした問題の厳密な答え」になっていること。計算した解は、各列と右辺を相対 $\widetilde\gamma_{mn}$ だけ動かした問題の厳密な最小二乗解）。
• 解の感度は $\varepsilon\bigl(2\kappa/\cos\theta+\kappa^2\tan\theta\bigr)+O(\varepsilon^2)$ です（$\theta$ は $\boldsymbol{b}$ と像のなす角）。残差が大きい問題では $\kappa^2$ が避けられず、残差が小さい問題では $\kappa$ で済みます。
• 正規方程式は $\kappa(A^{*}A)=\kappa(A)^2$ の問題を解くので、残差が小さくても $\kappa^2u$ の誤差を出しえます。数値例では次数 11 の多項式の当てはめで、QR の $8\times10^{-11}$ に対し正規方程式は $0.24$ でした。
• 階数落ちに近い問題では、打ち切り特異値分解（小さい特異値を捨てる方法）と Tikhonov 正則化（捨てずに滑らかに弱める方法。弱める度合いがフィルタ因子 $\sigma_i^2/(\sigma_i^2+\lambda^2)$）で小さい特異値の方向を抑えます。

この記事で使う事前知識

• 直交射影: ベクトルを部分空間の中で最も近い点に落とす操作。落とした先とのずれは部分空間に直交します。→ 内積空間とグラム・シュミット直交化
• 条件数: 入力の小さな変化が答えをどれだけ拡大するかの目安。大きいほど計算しにくい問題です。→ 条件数
• 特異値: 行列が空間を各方向にどれだけ伸ばすかを表す非負の数。特異値分解 $A=U\Sigma V^{*}$ の $\Sigma$ の対角成分です。→ 特異値分解（定理 4.1）［行列の分解と標準形］
• QR 分解: 行列を直交行列と上三角行列の積に分けること。→ QR 分解
• Householder 鏡映: 1 本のベクトルを鏡に映すように反転させる直交変換で、QR 分解の部品として使います。→ Householder 鏡映
• ユニタリ行列: ベクトルの長さを変えない行列（実数なら直交行列）。→ ユニタリ行列

1. 動機

観測点 $t_1,\ldots,t_m$ でのデータ $b_i$ に、未知の係数を $n$ 個もつモデルを当てはめるとき、方程式は $m>n$ 本で一般には解がありません。そこで「残差の 2 乗和を最小にする」係数を探すのが最小二乗法です。Gauss が天体の軌道決定に使って以来、データ解析・逆問題・統計の回帰の基本です。

数値計算の観点では、同じ最小二乗解を求めるのに 3 通りの道があります。正規方程式を Cholesky 分解で解く、QR 分解を使う、特異値分解を使う、の 3 つです。前者ほど速く、後者ほど頑健です。この記事では、どれがどんな場合にどれだけ精度を失うかを、問題の感度（条件数と残差の大きさ）とアルゴリズムの後退安定性に分けて説明します。行列のノルムと条件数 の物差しと、QR 分解 の後退誤差解析をそのまま使います。

2. 準備

$\mathbb{K}=\mathbb{R}$ または $\mathbb{C}$、$A\in\mathbb{K}^{m\times n}$、$\boldsymbol{b}\in\mathbb{K}^m$ とし、ノルムはすべて 2 ノルムで $\|\cdot\|$ と書きます。内積は $\langle\boldsymbol{x},\boldsymbol{y}\rangle=\boldsymbol{y}^{*}\boldsymbol{x}$ です。$A$ の特異値分解 $A=U\Sigma V^{*}$、$\operatorname{rank}A=r$、特異値 $\sigma_1\ge\cdots\ge\sigma_r>0$ は 特異値分解（定理 4.1）［行列の分解と標準形］ のとおりです。$U,V$ の列を $\boldsymbol{u}_i,\boldsymbol{v}_i$ と書きます。浮動小数の解析（第 5 節）では実行列を扱い、記号は QR 分解 と同じです。

列が一次独立（$r=n$）な $A$ について、条件数 $\kappa(A)=\sigma_1/\sigma_n$ と書きます。正方の正則行列ではこれは $\kappa_2(A)$ に一致します（特異行列までの距離（定理 6.4）［行列のノルムと線形方程式の条件数］）。

3. 最小二乗解の特徴づけ

まず「ずれが最小」を幾何で言い換えます。$A\boldsymbol{x}$ が動ける範囲は $A$ の像（列ベクトルが張る部分空間）で、$\boldsymbol{b}$ に最も近いのはそこへ垂線を下ろした足です。下の定理は、この見方と正規方程式が同じことを言っていると示します。

定義 3.1（最小二乗問題）↗

$\|A\boldsymbol{x}-\boldsymbol{b}\|$ を最小にする $\boldsymbol{x}\in\mathbb{K}^n$ を最小二乗解といい、$\boldsymbol{r}=\boldsymbol{b}-A\boldsymbol{x}$ をその残差という。$A^{*}A\boldsymbol{x}=A^{*}\boldsymbol{b}$ を正規方程式という。

定理 3.2（最小二乗解の特徴づけ）✓

1. 最小二乗解は存在する。
2. $\boldsymbol{x}$ について次は同値である: (a) $\boldsymbol{x}$ は最小二乗解、(b) $A\boldsymbol{x}=P\boldsymbol{b}$（$P$ は $A$ の像 $W=\operatorname{im}A$ への直交射影）、(c) 残差 $\boldsymbol{b}-A\boldsymbol{x}$ が $W$ に直交する、(d) $A^{*}A\boldsymbol{x}=A^{*}\boldsymbol{b}$。
3. 最小二乗解が一意 $\iff$ $\operatorname{rank}A=n$ $\iff$ $A^{*}A$ が正則。このとき解は $(A^{*}A)^{-1}A^{*}\boldsymbol{b}$。
4. 最小二乗解の全体は $\boldsymbol{x}_0+\ker A$（$\boldsymbol{x}_0$ は 1 つの最小二乗解）である。

証明（定理 3.2）

$W=\operatorname{im}A$ は $\mathbb{K}^m$ の有限次元部分空間です。

(a)⇔(b): $\{A\boldsymbol{x}\}$ は $W$ 全体を動くので、最良近似（定理 7.2）［内積空間とグラム・シュミット直交化］ により $\|A\boldsymbol{x}-\boldsymbol{b}\|$ が最小になるのは $A\boldsymbol{x}=P\boldsymbol{b}$ のとき、かつそのときに限ります。$P\boldsymbol{b}\in W$ なので $A\boldsymbol{x}=P\boldsymbol{b}$ となる $\boldsymbol{x}$ があり、1 が従います。

(b)⇔(c): 直交分解定理（定理 7.1）［内積空間とグラム・シュミット直交化］ の分解 $\boldsymbol{b}=\boldsymbol{w}+\boldsymbol{w}'$（$\boldsymbol{w}\in W$、$\boldsymbol{w}'\in W^{\perp}$）は一意で $\boldsymbol{w}=P\boldsymbol{b}$ です。$\boldsymbol{b}=A\boldsymbol{x}+(\boldsymbol{b}-A\boldsymbol{x})$ で $A\boldsymbol{x}\in W$ なので、$\boldsymbol{b}-A\boldsymbol{x}\in W^{\perp}$ と $A\boldsymbol{x}=P\boldsymbol{b}$ は同値です。

(c)⇔(d): $W$ は $A$ の列 $\boldsymbol{a}_j$ で張られるので、$\boldsymbol{y}\perp W$ は全 $j$ で $\boldsymbol{a}_j^{*}\boldsymbol{y}=0$、すなわち $A^{*}\boldsymbol{y}=\boldsymbol{0}$ と同値です。$\boldsymbol{y}=\boldsymbol{b}-A\boldsymbol{x}$ とすれば (d) です。

3: 最小二乗解は (b) の方程式 $A\boldsymbol{x}=P\boldsymbol{b}$ の解の全体で、その差は $\ker A$ を動きます（4 もこれで示された）。よって一意 $\iff\ker A=\{\boldsymbol{0}\}\iff\operatorname{rank}A=n$（次元定理（定理 8.3）［行列と連立一次方程式］）。また $A^{*}A\boldsymbol{x}=\boldsymbol{0}$ なら $\|A\boldsymbol{x}\|^2=\boldsymbol{x}^{*}A^{*}A\boldsymbol{x}=0$ なので $\ker A^{*}A=\ker A$ で、$n$ 次正方行列 $A^{*}A$ が正則 $\iff\ker A^{*}A=\{\boldsymbol{0}\}$ です（正則性の判定（定理 7.1）［行列式とその性質］）。このとき (d) から解は $(A^{*}A)^{-1}A^{*}\boldsymbol{b}$ です。

∎

残差 $\boldsymbol{r}$ と $A\boldsymbol{x}$ は直交するので、Pythagoras の関係 $\|\boldsymbol{b}\|^2=\|A\boldsymbol{x}\|^2+\|\boldsymbol{r}\|^2$ が成り立ちます。そこで $0\le\theta\le\pi/2$ を $\sin\theta=\|\boldsymbol{r}\|/\|\boldsymbol{b}\|$ で定めると $\|A\boldsymbol{x}\|=\|\boldsymbol{b}\|\cos\theta$ で、$\theta$ は $\boldsymbol{b}$ と像 $W$ のなす角です。$\theta$ が小さいほどデータはモデルによく合っています。

像 W を真横から見た図。b から W へ下ろした垂線の足が Ax（= Pb）で、垂線が残差 r。b と W のなす角が θ

4. 擬似逆と最小ノルム解

定義 4.1（擬似逆行列）↗

$A=U\Sigma V^{*}$ を特異値分解とし、$\Sigma^{+}\in\mathbb{R}^{n\times m}$ を $(\Sigma^{+})_{ii}=1/\sigma_i$（$i\le r$）、他の成分 0 の行列とする。$A^{+}=V\Sigma^{+}U^{*}$ を $A$ の擬似逆行列（Moore–Penrose 逆）という。

特異値分解は一意ではありませんが、$A^{+}$ が分解の選び方によらないことは 命題 4.3 で示します。

定理 4.2（擬似逆による最小ノルム解）✓

$\boldsymbol{x}^{+}=A^{+}\boldsymbol{b}$ は最小二乗解であり、最小二乗解の中でノルムが最小のものはこれただ 1 つである。具体的に

$$\boldsymbol{x}^{+}=\sum_{i=1}^r\frac{\boldsymbol{u}_i^{*}\boldsymbol{b}}{\sigma_i}\boldsymbol{v}_i,\qquad\min_{\boldsymbol{x}}\|A\boldsymbol{x}-\boldsymbol{b}\|^2=\sum_{i=r+1}^m|\boldsymbol{u}_i^{*}\boldsymbol{b}|^2.$$

証明（定理 4.2）

$\boldsymbol{y}=V^{*}\boldsymbol{x}$、$\boldsymbol{c}=U^{*}\boldsymbol{b}$ とおきます。ユニタリ行列の特徴づけ（命題 3.2）［スペクトル定理］ の 3 により $U^{*}$ はノルムを保ち、$V$ は $\mathbb{K}^n$ 上の全単射でノルムを保つので

$$\|A\boldsymbol{x}-\boldsymbol{b}\|^2=\|U^{*}(U\Sigma V^{*}\boldsymbol{x}-\boldsymbol{b})\|^2=\|\Sigma\boldsymbol{y}-\boldsymbol{c}\|^2=\sum_{i=1}^r|\sigma_iy_i-c_i|^2+\sum_{i=r+1}^m|c_i|^2,\qquad\|\boldsymbol{x}\|=\|\boldsymbol{y}\|$$

です。第 2 の和は $\boldsymbol{y}$ によらないので、最小は $y_i=c_i/\sigma_i$（$i\le r$）のとき、かつそのときに限り達成され、$y_i$（$i>r$）は任意です。そのうち $\|\boldsymbol{y}\|$ が最小なのは $y_i=0$（$i>r$）のときだけです。このとき $\boldsymbol{y}=\Sigma^{+}\boldsymbol{c}$、$\boldsymbol{x}=V\Sigma^{+}U^{*}\boldsymbol{b}=A^{+}\boldsymbol{b}$ で、成分で書けば主張の式です。

∎

次の命題は、擬似逆を特異値分解の取り方に頼らずに特徴づけるものです。4 つの等式を満たす行列がただ 1 つであることを示すので、証明 (2) では $X$ と $Y$ の 2 つが条件を満たすと仮定し、条件を 1 つずつ使って両方が同じ $XAY$ に等しいことを示します。

命題 4.3（Moore–Penrose の条件）✓

1. $X=A^{+}$ は次の 4 条件を満たす: $AXA=A$、$XAX=X$、$(AX)^{*}=AX$、$(XA)^{*}=XA$。
2. 4 条件を満たす $X\in\mathbb{K}^{n\times m}$ はただ 1 つである。とくに $A^{+}$ は特異値分解の選び方によらず、$A$ が実行列なら $A^{+}$ も実行列である。
3. $\operatorname{rank}A=n$ なら $A^{+}=(A^{*}A)^{-1}A^{*}$、$A^{+}A=I_n$。
4. 任意の $\boldsymbol{b}$ について $AA^{+}\boldsymbol{b}=P\boldsymbol{b}$（$P$ は $\operatorname{im}A$ への直交射影）。

証明（命題 4.3）

(1) $\Sigma\Sigma^{+}\Sigma=\Sigma$、$\Sigma^{+}\Sigma\Sigma^{+}=\Sigma^{+}$ で、$\Sigma\Sigma^{+}=\operatorname{diag}(I_r,0)\in\mathbb{R}^{m\times m}$、$\Sigma^{+}\Sigma=\operatorname{diag}(I_r,0)\in\mathbb{R}^{n\times n}$ は実対称です。$U^{*}U=I$、$V^{*}V=I$ を使うと $AA^{+}A=U\Sigma\Sigma^{+}\Sigma V^{*}=A$、$A^{+}AA^{+}=A^{+}$、$AA^{+}=U\Sigma\Sigma^{+}U^{*}$、$A^{+}A=V\Sigma^{+}\Sigma V^{*}$ はエルミートです。

(2) $X,Y$ がともに 4 条件を満たすとします。条件を順に使って

$$X=XAX=X(AX)^{*}=XX^{*}A^{*}=XX^{*}(AYA)^{*}=XX^{*}A^{*}Y^{*}A^{*}=X(AX)^{*}(AY)^{*}=XAXAY=XAY,$$$$Y=YAY=(YA)^{*}Y=A^{*}Y^{*}Y=(AXA)^{*}Y^{*}Y=A^{*}X^{*}A^{*}Y^{*}Y=(XA)^{*}(YA)^{*}Y=XAYAY=XAY$$

なので $X=Y$ です。別の特異値分解から作った $A^{+}$ も (1) により 4 条件を満たすので一致します。$A$ が実なら $\overline{A^{+}}$ も（条件の複素共役をとれば）4 条件を満たすので $\overline{A^{+}}=A^{+}$ です。

(3) 定理 3.2 の 3 により $M=A^{*}A$ は正則です。$X=M^{-1}A^{*}$ は $XA=I_n$ を満たし、$AXA=A$、$XAX=X$、$XA=I_n$ はエルミート、$AX=AM^{-1}A^{*}$ は $M^{-1}$ がエルミート（$M$ がエルミートだから）なのでエルミートです。(2) により $X=A^{+}$ です。

(4) 定理 4.2 により $A^{+}\boldsymbol{b}$ は最小二乗解で、定理 3.2 の 2 (b) から $AA^{+}\boldsymbol{b}=P\boldsymbol{b}$ です。

∎

命題 4.4（正規方程式の条件数）✓

$\operatorname{rank}A=n$ なら、$A^{*}A$ の特異値は $\sigma_1^2,\ldots,\sigma_n^2$ で、$\kappa_2(A^{*}A)=\kappa(A)^2$。また $\|A^{+}\|=1/\sigma_n$、$\|(A^{*}A)^{-1}\|=1/\sigma_n^2$。

証明（命題 4.4）

$A^{*}A=V\Sigma^{*}\Sigma V^{*}$ で、$\Sigma^{*}\Sigma=\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2)$ です。これは $A^{*}A$ の特異値分解なので、特異値は $\sigma_i^2$ です。特異行列までの距離（定理 6.4）［行列のノルムと線形方程式の条件数］ の 2 ノルムの部分から $\kappa_2(A^{*}A)=\sigma_1^2/\sigma_n^2=\kappa(A)^2$ です。2 ノルムの公式（命題 4.2）［行列のノルムと線形方程式の条件数］ から $\|A^{*}A\|=\sigma_1^2$ なので、$\|(A^{*}A)^{-1}\|=\kappa_2(A^{*}A)/\|A^{*}A\|=1/\sigma_n^2$ です。$A^{+}=V\Sigma^{+}U^{*}$ では ユニタリ不変性（命題 4.4）［行列のノルムと線形方程式の条件数］ から $\|A^{+}\|=\|\Sigma^{+}\|$ で、$\|\Sigma^{+}\boldsymbol{y}\|^2=\sum_{i\le n}|y_i|^2/\sigma_i^2\le\|\boldsymbol{y}\|^2/\sigma_n^2$（等号は $\boldsymbol{y}=\boldsymbol{e}_n$）なので $\|A^{+}\|=1/\sigma_n$ です。

∎

したがって $\kappa(A)=\|A\|\|A^{+}\|$ で、正方の場合の $\kappa_2$ の自然な拡張になっています。正規方程式を作った時点で条件数は 2 乗になります。$\kappa(A)=10^8$ なら $\kappa(A^{*}A)=10^{16}$ で、倍精度では正規方程式の解に 1 桁も保証がありません。

5. QR 分解による解法と後退安定性

この節では $A\in\mathbb{R}^{m\times n}$、$m\ge n$ とします。

定理 5.1（QR 分解による最小二乗解）✓

$\operatorname{rank}A=n$ とし、$A=QR$ を QR 分解（QR 分解の存在（定理 4.1）［QR 分解］）、$Q=(Q_1\ Q_2)$（$Q_1\in\mathbb{R}^{m\times n}$）、$R=\binom{R_1}{0}$ と分ける。このとき $R_1$ は正則で、最小二乗解は $\boldsymbol{x}=R_1^{-1}Q_1^{\mathsf{T}}\boldsymbol{b}$、最小の残差ノルムは $\|Q_2^{\mathsf{T}}\boldsymbol{b}\|$ である。

証明（定理 5.1）

$Q$ は正則なので $\operatorname{rank}R_1=\operatorname{rank}R=\operatorname{rank}A=n$ で、$R_1$ は正則です。直交行列 $Q^{\mathsf{T}}$ はノルムを保つので

$$\|A\boldsymbol{x}-\boldsymbol{b}\|^2=\|Q^{\mathsf{T}}A\boldsymbol{x}-Q^{\mathsf{T}}\boldsymbol{b}\|^2=\|R_1\boldsymbol{x}-Q_1^{\mathsf{T}}\boldsymbol{b}\|^2+\|Q_2^{\mathsf{T}}\boldsymbol{b}\|^2$$

です。第 2 項は $\boldsymbol{x}$ によらず、第 1 項は $\boldsymbol{x}=R_1^{-1}Q_1^{\mathsf{T}}\boldsymbol{b}$ のときだけ 0 になります。

∎

正規方程式 $R_1^{\mathsf{T}}R_1\boldsymbol{x}=R_1^{\mathsf{T}}Q_1^{\mathsf{T}}\boldsymbol{b}$ を作らずに、三角方程式 $R_1\boldsymbol{x}=Q_1^{\mathsf{T}}\boldsymbol{b}$ を直接解くのが要点です。条件数は $\kappa(R_1)=\kappa(A)$（$R_1=Q_1^{\mathsf{T}}A$ と $A$ は同じ特異値をもつ）で、2 乗になりません。実際の計算では $Q$ を作らず、Householder 鏡映を $\boldsymbol{b}$ に作用させます。

次の系は、Householder QR の後退安定性（QR 分解の章）を最小二乗問題へ持ち込むものです。$u$ は丸め誤差の単位、$\gamma_k$ は $k$ 回の演算で積もる誤差の見積もりで、条件 $86mnu<1$ は、これらの見積もりの分母が正になって式が意味をもつための仮定です。結論は「計算した解は、各列と右辺をごくわずか動かした問題の厳密な最小二乗解になっている」という形です。

系 5.2（Householder QR による最小二乗法の後退安定性）✓

$A\in\mathbb{F}^{m\times n}$、$\boldsymbol{b}\in\mathbb{F}^m$、$m\ge n$ とし、Householder QR で $\widehat R$ を、同じ鏡映を $\boldsymbol{b}$ に作用させて $\widehat{\boldsymbol{c}}$ を計算し、上三角の $\widehat R_1$（$\widehat R$ の最初の $n$ 行）について $\widehat R_1\boldsymbol{x}=\widehat{\boldsymbol{c}}_1$（$\widehat{\boldsymbol{c}}$ の最初の $n$ 成分）を後退代入で解いて $\widehat{\boldsymbol{x}}$ を得たとする。アンダーフロー・オーバーフローが起きず、$\widehat R_1$ の対角成分が 0 でなく、$86mnu<1$ とする。このとき $\Delta A$、$\Delta\boldsymbol{b}$ で

$$\|\Delta\boldsymbol{a}_j\|\le\gamma_{44mn}\|\boldsymbol{a}_j\|\ (j=1,\ldots,n),\qquad\|\Delta\boldsymbol{b}\|\le\gamma_{43mn}\|\boldsymbol{b}\|$$

を満たし、$\widehat{\boldsymbol{x}}$ が $\|(A+\Delta A)\boldsymbol{x}-(\boldsymbol{b}+\Delta\boldsymbol{b})\|$ を最小にするものがある。

証明（系 5.2）

Householder QR の後退安定性（定理 5.3）［QR 分解］ の 1・2 により、ある直交行列 $Q$ について $A+\Delta A_1=Q\widehat R$、$\|\Delta\boldsymbol{a}_{1,j}\|\le\gamma_{43mn}\|\boldsymbol{a}_j\|$、$\widehat{\boldsymbol{c}}=Q^{\mathsf{T}}(\boldsymbol{b}+\Delta\boldsymbol{b})$、$\|\Delta\boldsymbol{b}\|\le\gamma_{43mn}\|\boldsymbol{b}\|$ です（$Q$ は両方で同じ）。三角方程式の後退誤差（補題 3.1）［Gauss 消去と LU 分解］（上三角の後退代入）により $(\widehat R_1+\Delta R)\widehat{\boldsymbol{x}}=\widehat{\boldsymbol{c}}_1$、$|\Delta R|\le\gamma_n|\widehat R_1|$ です。

$\Delta A=\Delta A_1+Q\binom{\Delta R}{0}$ とおくと $A+\Delta A=Q\binom{\widehat R_1+\Delta R}{0}$ です。第 $j$ 列について、$\|Q\binom{\Delta\boldsymbol{r}_j}{\boldsymbol{0}}\|=\|\Delta\boldsymbol{r}_j\|\le\gamma_n\|\widehat{\boldsymbol{r}}_j\|$ で、$\|\widehat{\boldsymbol{r}}_j\|=\|Q\widehat{\boldsymbol{r}}_j\|=\|\boldsymbol{a}_j+\Delta\boldsymbol{a}_{1,j}\|\le(1+\gamma_{43mn})\|\boldsymbol{a}_j\|$ なので

$$\|\Delta\boldsymbol{a}_j\|\le\bigl(\gamma_{43mn}+\gamma_n(1+\gamma_{43mn})\bigr)\|\boldsymbol{a}_j\|\le\gamma_{43mn+n}\|\boldsymbol{a}_j\|\le\gamma_{44mn}\|\boldsymbol{a}_j\|$$

です（γₙ の計算規則（補題 3.2）［浮動小数点数と丸め誤差］）。最後に、任意の $\boldsymbol{x}$ について

$$\|(A+\Delta A)\boldsymbol{x}-(\boldsymbol{b}+\Delta\boldsymbol{b})\|^2=\|(\widehat R_1+\Delta R)\boldsymbol{x}-\widehat{\boldsymbol{c}}_1\|^2+\|\widehat{\boldsymbol{c}}_2\|^2\ge\|\widehat{\boldsymbol{c}}_2\|^2$$

で、$\boldsymbol{x}=\widehat{\boldsymbol{x}}$ で等号です。

∎

次の注意では 2 種類の誤差を使い分けます。後退誤差は「計算結果が、入力をどれだけ動かした問題の答えか」、前進誤差は「計算結果が本当の答えからどれだけ離れているか」です。後退安定でも、問題の感度（条件数）が大きければ前進誤差は大きくなりえます。

注意 5.3

正規方程式の精度. 正規方程式を浮動小数で作ると、$\mathrm{fl}(A^{\mathsf{T}}A)$ の各成分は内積の丸め誤差 $\gamma_m|\boldsymbol{a}_i|^{\mathsf{T}}|\boldsymbol{a}_j|$ を含み、Cholesky 分解は Cholesky 分解の後ろ向き誤差（定理 4.2）［行列計算の認証と最適化の証明書］ の意味で後退安定です。合わせると、計算した解はおおよそ、$A^{\mathsf{T}}A$ に相対 $O(mnu)$ の摂動を加えた（右辺 $A^{\mathsf{T}}\boldsymbol{b}$ の丸めも含む）正規方程式の厳密解とみなせ、線形方程式の摂動定理（定理 6.3）［行列のノルムと線形方程式の条件数］ を $A^{\mathsf{T}}A$ に使えば前進誤差は $\kappa(A^{\mathsf{T}}A)u=\kappa(A)^2u$ の程度になりえます。これは残差が 0 でも起こります。一方 QR による解法の前進誤差は、系 5.2 と次節の 定理 6.2 から $\kappa u+\kappa^2u\tan\theta$ の程度で、残差が小さい問題では $\kappa u$ にとどまります。$A^{\mathsf{T}}A$ に加わる摂動は、小さい $\Delta A$ による $(A+\Delta A)^{\mathsf{T}}(A+\Delta A)-A^{\mathsf{T}}A$ の形に書けるとは限りません。実際、もし後退安定なら残差 0 の問題の前進誤差は $\kappa u$ の程度にとどまるはずなので、$\kappa^2u$ の誤差が出る例（例 6.3）は、正規方程式が最小二乗問題として後退安定でないことを示しています。

6. 最小二乗問題の感度

この節では、$A$ や $\boldsymbol{b}$ を少し動かしたとき解がどれだけ動くかを調べます。解は $A^{+}\boldsymbol{b}$ なので、$A^{+}$ を $A$ の関数として微分すれば動きの大きさが分かります。そのための準備が最初の補題です。

補題 6.1（擬似逆の微分）✓

$A\in\mathbb{K}^{m\times n}$ を $\operatorname{rank}A=n$、$E\in\mathbb{K}^{m\times n}$ とする。ある $\varepsilon_0>0$ があって、$|\varepsilon|<\varepsilon_0$ の実数 $\varepsilon$ で $\operatorname{rank}(A+\varepsilon E)=n$ となり、$X(\varepsilon)=(A+\varepsilon E)^{+}$ の各成分の実部と虚部は $\varepsilon$ について何回でも微分可能で

$$X'(0)=-A^{+}EA^{+}+(A^{*}A)^{-1}E^{*}(I-AA^{+}).$$

証明（補題 6.1）

$M(\varepsilon)=(A+\varepsilon E)^{*}(A+\varepsilon E)=A^{*}A+\varepsilon(E^{*}A+A^{*}E)+\varepsilon^2E^{*}E$ の行列式は $\varepsilon$ の（複素係数の）多項式で、$\varepsilon=0$ で 0 でありません（定理 3.2 の 3）。多項式は連続なので、ある $\varepsilon_0>0$ について $|\varepsilon|<\varepsilon_0$ で $\det M(\varepsilon)\ne0$ です。すると 定理 3.2 の 3 から $\operatorname{rank}(A+\varepsilon E)=n$、命題 4.3 の 3 から $X(\varepsilon)=M(\varepsilon)^{-1}(A+\varepsilon E)^{*}$ です。余因子行列の基本関係式（命題 5.4）［行列式とその性質］ により $M(\varepsilon)^{-1}=\operatorname{adj}M(\varepsilon)/\det M(\varepsilon)$ の成分は、分母が 0 にならない多項式の商なので、$X(\varepsilon)$ の成分は何回でも微分できます。

$M(\varepsilon)X(\varepsilon)=(A+\varepsilon E)^{*}$ を $\varepsilon=0$ で微分すると（積の微分）$M'(0)X(0)+M(0)X'(0)=E^{*}$、すなわち

$$X'(0)=M^{-1}\bigl(E^{*}-(E^{*}A+A^{*}E)A^{+}\bigr)=M^{-1}E^{*}(I-AA^{+})-A^{+}EA^{+}$$

です（$M=A^{*}A$、$M^{-1}A^{*}=A^{+}$）。

∎

次の定理の $\theta$ は第 3 節の図の角です。残差が 0 なら $\theta=0$、データがモデルに合わないほど $\theta$ は $\pi/2$ に近づきます。評価の第 1 項 $2\kappa/\cos\theta$ は線形方程式と同じ型の増幅、第 2 項 $\kappa^2\tan\theta$ は残差が大きいときだけ効く増幅です。

定理 6.2（最小二乗問題の感度）✓

$\operatorname{rank}A=n$、$\boldsymbol{x}=A^{+}\boldsymbol{b}$、$\boldsymbol{r}=\boldsymbol{b}-A\boldsymbol{x}$ とし、$\theta$ を $\sin\theta=\|\boldsymbol{r}\|/\|\boldsymbol{b}\|$（$0\le\theta<\pi/2$ を仮定する、すなわち $A\boldsymbol{x}\ne\boldsymbol{0}$）で定める。$\|E\|\le\|A\|$、$\|\boldsymbol{f}\|\le\|\boldsymbol{b}\|$ とし、$\boldsymbol{x}(\varepsilon)=(A+\varepsilon E)^{+}(\boldsymbol{b}+\varepsilon\boldsymbol{f})$ とおく。このとき $\varepsilon\to0$ で

$$\frac{\|\boldsymbol{x}(\varepsilon)-\boldsymbol{x}\|}{\|\boldsymbol{x}\|}\le|\varepsilon|\left(\frac{2\kappa(A)}{\cos\theta}+\kappa(A)^2\tan\theta\right)+O(\varepsilon^2).$$

証明（定理 6.2）

補題 6.1 により $\boldsymbol{x}(\varepsilon)=X(\varepsilon)(\boldsymbol{b}+\varepsilon\boldsymbol{f})$ の各成分の実部・虚部 $g$ は $|\varepsilon|<\varepsilon_0$ で 2 回連続微分可能です。$0<\varepsilon_1<\varepsilon_0$ を固定すると、$g''$ は $[-\varepsilon_1,\varepsilon_1]$ で連続なので 最大値・最小値の定理（定理 2.3）［平均値の定理とテイラーの定理］ により有界で、Taylor の定理（定理 5.3）［平均値の定理とテイラーの定理］（$n=1$）から $|g(\varepsilon)-g(0)-\varepsilon g'(0)|\le C\varepsilon^2$ です。よって $\boldsymbol{x}(\varepsilon)-\boldsymbol{x}=\varepsilon\boldsymbol{x}'(0)+O(\varepsilon^2)$ で、$\|\boldsymbol{x}'(0)\|$ を評価すれば十分です。

補題 6.1 と $(I-AA^{+})\boldsymbol{b}=\boldsymbol{b}-P\boldsymbol{b}=\boldsymbol{r}$（命題 4.3 の 4、定理 3.2 の 2）から

$$\boldsymbol{x}'(0)=X'(0)\boldsymbol{b}+A^{+}\boldsymbol{f}=-A^{+}E\boldsymbol{x}+(A^{*}A)^{-1}E^{*}\boldsymbol{r}+A^{+}\boldsymbol{f}$$

です。命題 4.4 の $\|A^{+}\|=1/\sigma_n$、$\|(A^{*}A)^{-1}\|=1/\sigma_n^2$ と、$\|E^{*}\|=\|E\|\le\sigma_1$（2 ノルムの公式（命題 4.2）［行列のノルムと線形方程式の条件数］）、$\|\boldsymbol{f}\|\le\|\boldsymbol{b}\|$ を使うと

$$\frac{\|\boldsymbol{x}'(0)\|}{\|\boldsymbol{x}\|}\le\frac{\sigma_1}{\sigma_n}+\frac{\sigma_1\|\boldsymbol{r}\|}{\sigma_n^2\|\boldsymbol{x}\|}+\frac{\|\boldsymbol{b}\|}{\sigma_n\|\boldsymbol{x}\|}=\kappa+\kappa^2\frac{\|\boldsymbol{r}\|}{\sigma_1\|\boldsymbol{x}\|}+\kappa\frac{\|\boldsymbol{b}\|}{\sigma_1\|\boldsymbol{x}\|}$$

です。$\sigma_1\|\boldsymbol{x}\|\ge\|A\boldsymbol{x}\|=\|\boldsymbol{b}\|\cos\theta>0$（第 3 節の Pythagoras の関係）なので、右辺は $\kappa+\kappa^2\tan\theta+\kappa/\cos\theta\le2\kappa/\cos\theta+\kappa^2\tan\theta$ 以下です。

∎

この評価の読み方は次のとおりです。

• 残差が 0（$\theta=0$）なら $2\kappa$ で、正方の線形方程式と同じ $\kappa$ の程度です。
• 残差が大きいと $\kappa^2\tan\theta$ が効きます。これは $A$ の摂動で像の向きが $\kappa$ に比例して傾き、それが大きな残差を通して解に $\kappa$ 倍で伝わるためで、どんな算法でも避けられない問題そのものの感度です。
• $\theta$ が $\pi/2$ に近い（データがモデルにほとんど合わない）と $1/\cos\theta$ が効きます。

後退安定な算法（系 5.2、$\varepsilon\approx\widetilde\gamma_{mn}$）の前進誤差は、この評価から $\widetilde\gamma_{mn}(2\kappa/\cos\theta+\kappa^2\tan\theta)$ の程度と見積もれます。

例 6.3（多項式の当てはめ：3 つの解法の精度）

$[0,1]$ の等間隔 100 点で、単項式基底 $1,t,\ldots,t^{n-1}$ の Vandermonde 行列 $A$ に当てはめます。真の係数をすべて 1 とし、残差 0 の右辺 $\boldsymbol{b}=A\boldsymbol{1}$ と、像に直交する成分（乱数ベクトルの射影、mpmath の 50 桁で計算）を加えた残差の大きい右辺を比べます。参照解は浮動小数のデータ $(A,\boldsymbol{b})$ に対する mpmath の 50 桁の qr_solve です。

# Python 3.12.10 / NumPy 2.5.2 / SciPy 1.18.1 / mpmath 1.3.0
# 乱数種 20261006
import numpy as np, mpmath as mp
from scipy.linalg import cho_factor, cho_solve, solve_triangular
mp.mp.dps = 50
rng = np.random.default_rng(20261006)
t = np.linspace(0, 1, 100)
print(" n  kappa    tan    normal   QR       lstsq")
for n in (4, 8, 12):
    A = np.vander(t, n, increasing=True)        # 単項式基底
    Am = mp.matrix(A.tolist())
    z = mp.matrix(rng.standard_normal(100).tolist())
    rp = z - Am * mp.qr_solve(Am, z)[0]          # 像に直交する成分
    for s in (0.0, 1.0):                          # 残差 0 と大きい残差
        b = A @ np.ones(n) + s * np.array(rp.tolist(), dtype=float)[:, 0]
        bm = mp.matrix(b.tolist())
        xr = np.array(mp.qr_solve(Am, bm)[0].tolist(), dtype=float)[:, 0]
        x1 = cho_solve(cho_factor(A.T @ A), A.T @ b)   # 正規方程式
        Q, R = np.linalg.qr(A)
        x2 = solve_triangular(R, Q.T @ b)               # Householder QR
        x3 = np.linalg.lstsq(A, b, rcond=None)[0]       # SVD
        r = b - A @ xr
        tan = np.linalg.norm(r) / np.linalg.norm(b - r)
        err = [np.linalg.norm(x - xr) / np.linalg.norm(xr)
               for x in (x1, x2, x3)]
        print(f"{n:2d}  {np.linalg.cond(A):.1e}  {tan:.1e}  "
              + "  ".join(f"{e:.1e}" for e in err))

 n  kappa    tan    normal   QR       lstsq
 4  1.2e+02  8.0e-17  8.3e-13  5.5e-15  1.3e-15
 4  1.2e+02  4.4e-01  5.7e-13  2.1e-14  2.3e-14
 8  1.2e+05  5.1e-17  1.5e-07  4.1e-12  4.4e-12
 8  1.2e+05  2.7e-01  2.0e-07  4.4e-09  4.4e-09
12  1.2e+08  9.0e-17  2.4e-01  8.3e-11  6.5e-09
12  1.2e+08  2.2e-01  4.1e-01  2.8e-03  2.8e-03

検算. (a) 既知値: 残差 0 の場合、参照解は真の係数（すべて 1）に浮動小数のデータの丸めの影響を加えたものです。(b) 別精度・別手法: 参照解は mpmath の 50 桁の QR、比べた 3 解法は Cholesky（正規方程式）・Householder QR・特異値分解（lstsq）です。(c) 次数: $\kappa$ が $10^3$ 倍になるごとに、正規方程式の誤差は残差によらずおよそ $10^5$〜$10^6$ 倍（$\kappa^2u$ に比例、$n=12$ で $\kappa^2u\approx1.6$ に達して桁が残らない）、QR の誤差は残差 0 では $\kappa u$ 以下（$n=12$ で $\kappa u\approx1.3\times10^{-8}$ に対し $8.3\times10^{-11}$）、残差が大きいと $n=8\to12$ で $6\times10^5$ 倍と $\kappa^2$ の増え方に移ります。QR と lstsq の誤差はどれも 定理 6.2 の評価 $u(2\kappa/\cos\theta+\kappa^2\tan\theta)$ の範囲内です（$n=12$ の残差 0 で lstsq の誤差が QR より大きいのも、$\kappa u$ の範囲内のばらつきです）。

7. 階数落ちに近い問題と正則化

$\sigma_n$ がほとんど 0 だと、$\boldsymbol{x}^{+}=\sum_i(\boldsymbol{u}_i^{*}\boldsymbol{b}/\sigma_i)\boldsymbol{v}_i$ の小さい $\sigma_i$ の項が、$\boldsymbol{b}$ に含まれる雑音を $1/\sigma_i$ 倍に拡大します。そこで、小さい特異値の方向を捨てるか、抑えるかします。

命題 7.1（打ち切り特異値分解の解）✓

$1\le k\le r$ とし、$A_k=\sum_{i\le k}\sigma_i\boldsymbol{u}_i\boldsymbol{v}_i^{*}$ とする。$\boldsymbol{x}_k:=A_k^{+}\boldsymbol{b}$ は $\|A_k\boldsymbol{x}-\boldsymbol{b}\|$ の最小ノルムの最小二乗解で、

$$\boldsymbol{x}_k=\sum_{i=1}^k\frac{\boldsymbol{u}_i^{*}\boldsymbol{b}}{\sigma_i}\boldsymbol{v}_i,\qquad\|\boldsymbol{x}_k\|^2=\sum_{i=1}^k\frac{|\boldsymbol{u}_i^{*}\boldsymbol{b}|^2}{\sigma_i^2},\qquad\|\boldsymbol{b}-A\boldsymbol{x}_k\|^2=\sum_{i=k+1}^m|\boldsymbol{u}_i^{*}\boldsymbol{b}|^2.$$

とくに $\|\boldsymbol{x}_k\|$ は $k$ について非減少、残差は非増加である。

証明（命題 7.1）

$\Sigma_k$ を $\Sigma$ の第 $k+1$ 以降の対角成分を 0 にした行列とすると、$A_k=U\Sigma_kV^{*}$ は $A_k$ の特異値分解です。定理 4.2 を $A_k$ に使えば最初の 2 式が出ます（2 番目は $\boldsymbol{v}_i$ の正規直交性）。$A\boldsymbol{v}_i=\sigma_i\boldsymbol{u}_i$ より $A\boldsymbol{x}_k=\sum_{i\le k}(\boldsymbol{u}_i^{*}\boldsymbol{b})\boldsymbol{u}_i$ で、$U$ の列は正規直交基底（ユニタリ行列の特徴づけ（命題 3.2）［スペクトル定理］ の 4）なので $\boldsymbol{b}=\sum_{i\le m}(\boldsymbol{u}_i^{*}\boldsymbol{b})\boldsymbol{u}_i$、差のノルムが 3 番目の式です。

∎

$A_k$ は Eckart–Young の定理（命題 4.3）［行列の分解と標準形］ により階数 $k$ 以下で $A$ に最も近い行列で、$\|A-A_k\|=\sigma_{k+1}$ です。「$\sigma_{k+1}$ 以下の違いはデータの誤差の範囲で区別できない」と判断できるとき、$A$ を $A_k$ に置き換えるのが打ち切り特異値分解の考え方です。numpy.linalg.lstsq の rcond は、$\sigma_i\le\mathtt{rcond}\cdot\sigma_1$ の特異値を捨てる閾値です。

命題 7.2（Tikhonov 正則化とフィルタ因子）✓

$\lambda>0$ とする。$\|A\boldsymbol{x}-\boldsymbol{b}\|^2+\lambda^2\|\boldsymbol{x}\|^2$ を最小にする $\boldsymbol{x}$ はただ 1 つで、

$$\boldsymbol{x}_\lambda=(A^{*}A+\lambda^2I)^{-1}A^{*}\boldsymbol{b}=\sum_{i=1}^r\frac{\sigma_i^2}{\sigma_i^2+\lambda^2}\cdot\frac{\boldsymbol{u}_i^{*}\boldsymbol{b}}{\sigma_i}\boldsymbol{v}_i.$$

さらに $\|\boldsymbol{x}_\lambda\|\le\|\boldsymbol{x}^{+}\|$ かつ $\|\boldsymbol{x}_\lambda\|\le\|\boldsymbol{b}\|/(2\lambda)$。

証明（命題 7.2）

$K=\binom{A}{\lambda I}\in\mathbb{K}^{(m+n)\times n}$、$\boldsymbol{d}=\binom{\boldsymbol{b}}{\boldsymbol{0}}$ とおくと $\|K\boldsymbol{x}-\boldsymbol{d}\|^2=\|A\boldsymbol{x}-\boldsymbol{b}\|^2+\lambda^2\|\boldsymbol{x}\|^2$ です。$K\boldsymbol{x}=\boldsymbol{0}$ なら $\lambda\boldsymbol{x}=\boldsymbol{0}$ なので $\operatorname{rank}K=n$ で、定理 3.2 の 3 により最小解はただ 1 つ、$(K^{*}K)^{-1}K^{*}\boldsymbol{d}=(A^{*}A+\lambda^2I)^{-1}A^{*}\boldsymbol{b}$ です。特異値分解から $A^{*}A+\lambda^2I=V(\Sigma^{*}\Sigma+\lambda^2I)V^{*}$、$A^{*}\boldsymbol{b}=V\Sigma^{*}U^{*}\boldsymbol{b}=\sum_{i\le r}\sigma_i(\boldsymbol{u}_i^{*}\boldsymbol{b})\boldsymbol{v}_i$ なので、$\boldsymbol{v}_i$ の係数は $\sigma_i(\boldsymbol{u}_i^{*}\boldsymbol{b})/(\sigma_i^2+\lambda^2)$ で、主張の式です。フィルタ因子 $\sigma_i^2/(\sigma_i^2+\lambda^2)$ は $1$ 未満なので $\|\boldsymbol{x}_\lambda\|\le\|\boldsymbol{x}^{+}\|$（定理 4.2 の式と比べる）。また $\sigma_i^2+\lambda^2\ge2\sigma_i\lambda$ から係数の絶対値は $|\boldsymbol{u}_i^{*}\boldsymbol{b}|/(2\lambda)$ 以下で、$\|\boldsymbol{x}_\lambda\|^2\le\sum_i|\boldsymbol{u}_i^{*}\boldsymbol{b}|^2/(4\lambda^2)\le\|\boldsymbol{b}\|^2/(4\lambda^2)$ です。

∎

フィルタ因子は $\sigma_i\gg\lambda$ でほぼ 1、$\sigma_i\ll\lambda$ でほぼ $\sigma_i^2/\lambda^2$ です。打ち切り特異値分解の「0 か 1 か」のフィルタを滑らかにしたもので、$\lambda$ が打ち切りの閾値の役をします。計算では $K$ に QR 分解を使えば、$A^{*}A$ を作らずに済みます。

例 7.3（打ち切りと Tikhonov 正則化の比較）

特異値を $1,10^{-1},\ldots,10^{-7}$ と $10^{-10},10^{-12}$ にした $30\times10$ の行列を作り、真の解を上位 8 個の右特異ベクトルの和（ノルム $\sqrt8\approx2.83$）、右辺に大きさ $10^{-8}$ の雑音を加えます。lstsq の rcond と Tikhonov の $\lambda$ を変え、解のノルム・残差・真の解との誤差、および 命題 7.1・命題 7.2 の式で構成に使った $U,\Sigma,V$ から直接計算した解との相対差（formula）を並べます。

# Python 3.12.10 / NumPy 2.5.2 ／乱数種 20261006
import numpy as np
rng = np.random.default_rng(20261006)
m, n = 30, 10
U, _ = np.linalg.qr(rng.standard_normal((m, n)))
V, _ = np.linalg.qr(rng.standard_normal((n, n)))
s = np.r_[np.logspace(0, -7, 8), 1e-10, 1e-12]   # 特異値
A = U @ np.diag(s) @ V.T
xt = V[:, :8] @ np.ones(8)          # 真の解（上位 8 方向）
b = A @ xt + 1e-8 * rng.standard_normal(m)       # 雑音つき
nr, c = np.linalg.norm, U.T @ b
print("method     param   |x|       |Ax-b|    |x-xt|    formula")
for rc, k in ((None, 10), (1e-11, 9), (1e-9, 8)):
    x = np.linalg.lstsq(A, b, rcond=rc)[0]
    xf = V[:, :k] @ (c[:k] / s[:k])               # 打ち切り SVD の式
    print(f"lstsq   {str(rc):>7}  {nr(x):.2e}  {nr(A@x-b):.2e}"
          f"  {nr(x-xt):.2e}  {nr(x-xf)/nr(xf):.0e}")
for lam in (1e-12, 1e-9, 1e-8, 1e-6):
    K = np.vstack([A, lam * np.eye(n)])          # [A; lam I]
    x = np.linalg.lstsq(K, np.r_[b, np.zeros(n)], rcond=None)[0]
    xf = V @ (s**2 / (s**2 + lam**2) * c / s)    # フィルタ因子の式
    print(f"Tikhonov {lam:.0e}  {nr(x):.2e}  {nr(A@x-b):.2e}"
          f"  {nr(x-xt):.2e}  {nr(x-xf)/nr(xf):.0e}")

method     param   |x|       |Ax-b|    |x-xt|    formula
lstsq      None  3.99e+02  4.89e-08  3.99e+02  1e-03
lstsq     1e-11  2.33e+02  4.89e-08  2.33e+02  4e-07
lstsq     1e-09  2.87e+00  5.42e-08  1.01e-01  2e-10
Tikhonov 1e-12  2.84e+02  4.89e-08  2.83e+02  8e-04
Tikhonov 1e-09  3.68e+00  5.41e-08  2.31e+00  3e-07
Tikhonov 1e-08  2.86e+00  5.42e-08  9.34e-02  3e-09
Tikhonov 1e-06  2.50e+00  5.29e-07  1.11e+00  2e-11

すべての特異値を使う解（rcond=None）は残差が最小ですが、ノルムが 400 に膨らみ、真の解から大きく外れます。$10^{-10},10^{-12}$ の方向を捨てる rcond=1e-9 では誤差が $0.1$（雑音 $10^{-8}$ を $\sigma_8=10^{-7}$ で割った大きさ）まで下がり、残差はほとんど増えません。Tikhonov では $\lambda=10^{-8}$ 付近が最もよく、$\lambda$ を大きくしすぎると（$10^{-6}$）上位の方向まで縮められて残差と誤差が増えます。

検算. (a) 既知値: 特異値分解を構成から知っているので、命題 7.1・命題 7.2 の式の値（formula 列の比較相手）が計算できます。(b) 別手法: lstsq（特異値分解）と、拡大行列 $K$ に対する lstsq の結果が式の値と一致します。(c) 次数: 一致の度合いは $u/\sigma_{\min}$（残した最小の特異値）の程度で、捨てる特異値を増やすと $10^{-3}\to4\times10^{-7}\to2\times10^{-10}$ と改善します。浮動小数で $A$ を作った時点で特異値が絶対 $u\|A\|$ 程度動くので、$\sigma=10^{-12}$ の方向は相対 $10^{-4}$ の不確かさをもつ、という見積もりと合っています。

8. 演習

演習 8.1易

3 点 $(t,y)=(0,1),(1,2),(2,2)$ に直線 $y=c_0+c_1t$ を最小二乗で当てはめ、係数と残差を求めよ。残差が $A$ の像に直交することを確かめよ。

解答

$A=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}$、$\boldsymbol{b}=(1,2,2)^{\mathsf{T}}$ で、$A^{\mathsf{T}}A=\begin{pmatrix}3&3\\3&5\end{pmatrix}$、$A^{\mathsf{T}}\boldsymbol{b}=(5,6)^{\mathsf{T}}$。正規方程式（定理 3.2）を解いて $c_0=7/6$、$c_1=1/2$。残差は $\boldsymbol{r}=(-1/6,1/3,-1/6)^{\mathsf{T}}$ で、$A^{\mathsf{T}}\boldsymbol{r}=(0,0)^{\mathsf{T}}$、すなわち両方の列に直交します。$\|\boldsymbol{r}\|^2=1/6$。

演習 8.2標準

$\boldsymbol{a}\in\mathbb{K}^m$、$\boldsymbol{c}\in\mathbb{K}^n$ を 0 でないベクトルとし、$A=\boldsymbol{a}\boldsymbol{c}^{*}$ とする。$A^{+}=\boldsymbol{c}\boldsymbol{a}^{*}/(\|\boldsymbol{a}\|^2\|\boldsymbol{c}\|^2)$ であることを示せ。

解答

$X=\boldsymbol{c}\boldsymbol{a}^{*}/(\|\boldsymbol{a}\|^2\|\boldsymbol{c}\|^2)$ とおくと、$AX=\boldsymbol{a}\boldsymbol{a}^{*}/\|\boldsymbol{a}\|^2$、$XA=\boldsymbol{c}\boldsymbol{c}^{*}/\|\boldsymbol{c}\|^2$ はエルミートで、$AXA=(\boldsymbol{a}\boldsymbol{a}^{*}/\|\boldsymbol{a}\|^2)\boldsymbol{a}\boldsymbol{c}^{*}=A$、$XAX=(\boldsymbol{c}\boldsymbol{c}^{*}/\|\boldsymbol{c}\|^2)X=X$ です。命題 4.3 の 2 の一意性から $X=A^{+}$ です。

演習 8.3難

$0<\delta<1$、$A=\begin{pmatrix}1&0\\0&\delta\\0&0\end{pmatrix}$、$\boldsymbol{b}=(1,0,1)^{\mathsf{T}}$ とし、$A$ の $(3,2)$ 成分だけを $\varepsilon$ に変えた $A+\varepsilon E$（$\|E\|=1=\|A\|$）の最小二乗解を求めよ。定理 6.2 の $\kappa^2\tan\theta$ の項が改善できないことを確かめよ。

解答

$\boldsymbol{x}=(1,0)^{\mathsf{T}}$、$\boldsymbol{r}=(0,0,1)^{\mathsf{T}}$、$\|A\boldsymbol{x}\|=\|\boldsymbol{r}\|=1$ なので $\tan\theta=1$、$\kappa=1/\delta$ です。摂動後の正規方程式は対角で、第 2 成分は $(\delta\cdot0+\varepsilon\cdot1)/(\delta^2+\varepsilon^2)$、第 1 成分は $1$ のままです。よって

$$\frac{\|\boldsymbol{x}(\varepsilon)-\boldsymbol{x}\|}{\|\boldsymbol{x}\|}=\frac{\varepsilon}{\delta^2+\varepsilon^2}=\varepsilon\kappa^2+O(\varepsilon^3)$$

で、$\varepsilon$ の 1 次の係数は $\kappa^2\tan\theta$ に一致します。評価の $2\kappa/\cos\theta$ の部分を除いて、$\kappa^2\tan\theta$ の項は定数倍も含めて改善できません。

参考文献

• N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002 — Ch. 20（最小二乗問題の感度と QR 法の後退安定性）。
• Å. Björck, Numerical Methods for Least Squares Problems, SIAM, 1996 — Ch. 1–2。
• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013 — §5.3–5.5、§6.1。
• L. N. Trefethen and D. Bau III, Numerical Linear Algebra, SIAM, 1997 — Lectures 11, 18, 19。
• P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, 1998 — Ch. 3–5（打ち切り特異値分解、フィルタ因子、Tikhonov 正則化）。
• R. Penrose, “A generalized inverse for matrices”, Proc. Cambridge Philos. Soc. 51 (1955), 406–413.