1 最小二乗解
定義 1.1. A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n 、b ∈ R m b\in\R^m b ∈ R m 、x ∈ R n x\in\R^n x ∈ R n に対して、b − A x ∈ R m b-Ax\in\R^m b − A x ∈ R m をA x = b Ax=b A x = b に関するx x x の 残差 (residual ) という。
命題 1.2. A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n 、b ∈ R m b\in\R^m b ∈ R m とする。x ∈ R n x\in\R^n x ∈ R n について、次の三条件は同値である。
x x x はA x = b Ax=b A x = b の最小二乗解である。すなわち、任意のy ∈ R n y\in\R^n y ∈ R n に対して∥ b − A x ∥ 2 ≤ ∥ b − A y ∥ 2 \lVert b-Ax\rVert_2\le\lVert b-Ay\rVert_2 ∥ b − A x ∥ 2 ≤ ∥ b − A y ∥ 2 が成り立つ。
残差b − A x b-Ax b − A x はA A A の列空間{ A y ∣ y ∈ R n } \{Ay\mid y\in\R^n\} { A y ∣ y ∈ R n } のすべての元に直交する。
A T A x = A T b A^{\mathsf T}Ax=A^{\mathsf T}b A T A x = A T b が成り立つ。
さらに次が成り立つ。
A x = b Ax=b A x = b の最小二乗解は存在し、その一つをx 0 x_0 x 0 とすると、最小二乗解の全体はx 0 + ker A = { x 0 + z ∣ z ∈ ker A } x_0+\ker A=\{x_0+z\mid z\in\ker A\} x 0 + ker A = { x 0 + z ∣ z ∈ ker A } である。
A x = b Ax=b A x = b の最小二乗解がただ一つであることと、A A A の列が一次独立であることは同値である。A A A の列が一次独立ならば、A T A A^{\mathsf T}A A T A は正則であり、ただ一つの最小二乗解は( A T A ) − 1 A T b (A^{\mathsf T}A)^{-1}A^{\mathsf T}b ( A T A ) − 1 A T b である。
証明. §D3.17 定理 2.2 により、条件 (a) と条件 (c) は同値であり、(1) と、(2) の前半が成り立つ。A A A の列空間はA A A の列a 1 , … , a n a_1,\dots,a_n a 1 , … , a n で張られるから、条件 (b) はa j T ( b − A x ) = 0 a_j^{\mathsf T}(b-Ax)=0 a j T ( b − A x ) = 0 (1 ≤ j ≤ n 1\le j\le n 1 ≤ j ≤ n )と同値であり、これはA T ( b − A x ) = 0 A^{\mathsf T}(b-Ax)=0 A T ( b − A x ) = 0 、すなわち条件 (c) と同値である。A A A の列が一次独立ならばker A = { 0 } \ker A=\{0\} ker A = { 0 } であり、§D3.18 補題 1.1 によりker ( A T A ) = ker A = { 0 } \ker(A^{\mathsf T}A)=\ker A=\{0\} ker ( A T A ) = ker A = { 0 } であるから、A T A A^{\mathsf T}A A T A は正則である。このとき条件 (c) の解は( A T A ) − 1 A T b (A^{\mathsf T}A)^{-1}A^{\mathsf T}b ( A T A ) − 1 A T b だけである。▨
2 特異値と条件数
補題 2.1. R n \R^n R n とR m \R^m R m に Euclid ノルム∥ ⋅ ∥ 2 \lVert\cdot\rVert_2 ∥ ⋅ ∥ 2 を入れ、M ∈ R m × n M\in\R^{m\times n} M ∈ R m × n の作用素ノルム(§E20.2 定義 1.1 )を∥ M ∥ 2 \lVert M\rVert_2 ∥ M ∥ 2 と書く。p : = min { m , n } p:=\min\{m,n\} p := min { m , n } とし、M M M の特異値(§D3.18 定理 1.2 )をσ 1 ( M ) ≥ ⋯ ≥ σ p ( M ) \sigma_1(M)\ge\dots\ge\sigma_p(M) σ 1 ( M ) ≥ ⋯ ≥ σ p ( M ) と書く。
∥ M ∥ 2 = σ 1 ( M ) = ∥ M T ∥ 2 ≤ ∥ M ∥ F \lVert M\rVert_2=\sigma_1(M)=\lVert M^{\mathsf T}\rVert_2\le\lVert M\rVert_F ∥ M ∥ 2 = σ 1 ( M ) = ∥ M T ∥ 2 ≤ ∥ M ∥ F が成り立つ。
m ≥ n m\ge n m ≥ n でありM M M の列が一次独立であるとする。σ n ( M ) > 0 \sigma_n(M)>0 σ n ( M ) > 0 であり、任意のh ∈ R n h\in\R^n h ∈ R n に対して∥ M h ∥ 2 ≥ σ n ( M ) ∥ h ∥ 2 \lVert Mh\rVert_2\ge\sigma_n(M)\lVert h\rVert_2 ∥ M h ∥ 2 ≥ σ n ( M ) ∥ h ∥ 2 が成り立ち、等号を満たすh ≠ 0 h\ne0 h = 0 が存在する。さらに∥ ( M T M ) − 1 ∥ 2 = σ n ( M ) − 2 \lVert(M^{\mathsf T}M)^{-1}\rVert_2=\sigma_n(M)^{-2} ∥( M T M ) − 1 ∥ 2 = σ n ( M ) − 2 、∥ ( M T M ) − 1 M T ∥ 2 = σ n ( M ) − 1 \lVert(M^{\mathsf T}M)^{-1}M^{\mathsf T}\rVert_2=\sigma_n(M)^{-1} ∥( M T M ) − 1 M T ∥ 2 = σ n ( M ) − 1 である。
m ≥ n m\ge n m ≥ n でありM M M の列が一次独立であるとし、E ∈ R m × n E\in\R^{m\times n} E ∈ R m × n が∥ E ∥ 2 < σ n ( M ) \lVert E\rVert_2<\sigma_n(M) ∥ E ∥ 2 < σ n ( M ) を満たすとする。このときM + E M+E M + E の列は一次独立であり、σ n ( M + E ) ≥ σ n ( M ) − ∥ E ∥ 2 \sigma_n(M+E)\ge\sigma_n(M)-\lVert E\rVert_2 σ n ( M + E ) ≥ σ n ( M ) − ∥ E ∥ 2 が成り立つ。
証明. §D3.18 定理 1.2 により、直交行列U ∈ R m × m U\in\R^{m\times m} U ∈ R m × m 、V ∈ R n × n V\in\R^{n\times n} V ∈ R n × n と、( i , i ) (i,i) ( i , i ) 成分がσ i ( M ) \sigma_i(M) σ i ( M ) (1 ≤ i ≤ p 1\le i\le p 1 ≤ i ≤ p )でその他の成分が0 0 0 のΣ ∈ R m × n \Sigma\in\R^{m\times n} Σ ∈ R m × n が存在してM = U Σ V T M=U\Sigma V^{\mathsf T} M = U Σ V T である。h ∈ R n h\in\R^n h ∈ R n に対してk : = V T h k:=V^{\mathsf T}h k := V T h と置くと、直交行列は Euclid ノルムを保つので∥ k ∥ 2 = ∥ h ∥ 2 \lVert k\rVert_2=\lVert h\rVert_2 ∥ k ∥ 2 = ∥ h ∥ 2 であり、
∥ M h ∥ 2 2 = ∥ Σ k ∥ 2 2 = ∑ i = 1 p σ i ( M ) 2 k i 2 \lVert Mh\rVert_2^2=\lVert\Sigma k\rVert_2^2=\sum_{i=1}^p\sigma_i(M)^2k_i^2 ∥ M h ∥ 2 2 = ∥ Σ k ∥ 2 2 = i = 1 ∑ p σ i ( M ) 2 k i 2 である。
右辺はσ 1 ( M ) 2 ∥ k ∥ 2 2 \sigma_1(M)^2\lVert k\rVert_2^2 σ 1 ( M ) 2 ∥ k ∥ 2 2 以下であり、h h h をV V V の第1 1 1 列とすると等号が成り立つから、∥ M ∥ 2 = σ 1 ( M ) \lVert M\rVert_2=\sigma_1(M) ∥ M ∥ 2 = σ 1 ( M ) である。M T = V Σ T U T M^{\mathsf T}=V\Sigma^{\mathsf T}U^{\mathsf T} M T = V Σ T U T は§D3.18 定理 1.2 の形のM T M^{\mathsf T} M T の分解であるから、§D3.18 命題 2.1 によりM T M^{\mathsf T} M T の特異値はσ 1 ( M ) , … , σ p ( M ) \sigma_1(M),\dots,\sigma_p(M) σ 1 ( M ) , … , σ p ( M ) であり、∥ M T ∥ 2 = σ 1 ( M ) \lVert M^{\mathsf T}\rVert_2=\sigma_1(M) ∥ M T ∥ 2 = σ 1 ( M ) である。M M M の第i i i 行をμ i \mu_i μ i とすると、Cauchy–Schwarz の不等式により∥ M h ∥ 2 2 = ∑ i ( μ i h ) 2 ≤ ∑ i ∥ μ i ∥ 2 2 ∥ h ∥ 2 2 = ∥ M ∥ F 2 ∥ h ∥ 2 2 \lVert Mh\rVert_2^2=\sum_i(\mu_ih)^2\le\sum_i\lVert\mu_i\rVert_2^2\lVert h\rVert_2^2=\lVert M\rVert_F^2\lVert h\rVert_2^2 ∥ M h ∥ 2 2 = ∑ i ( μ i h ) 2 ≤ ∑ i ∥ μ i ∥ 2 2 ∥ h ∥ 2 2 = ∥ M ∥ F 2 ∥ h ∥ 2 2 であり、∥ M ∥ 2 ≤ ∥ M ∥ F \lVert M\rVert_2\le\lVert M\rVert_F ∥ M ∥ 2 ≤ ∥ M ∥ F である。これで(1) は示された。
m ≥ n m\ge n m ≥ n でありM M M の列が一次独立ならば、M M M の階数はn n n であり、§D3.18 定理 1.2 によりσ n ( M ) > 0 \sigma_n(M)>0 σ n ( M ) > 0 である。p = n p=n p = n であるから、上の等式の右辺はσ n ( M ) 2 ∥ k ∥ 2 2 \sigma_n(M)^2\lVert k\rVert_2^2 σ n ( M ) 2 ∥ k ∥ 2 2 以上であり、h h h をV V V の第n n n 列とすると等号が成り立つ。M T M = V diag ( σ 1 ( M ) 2 , … , σ n ( M ) 2 ) V T M^{\mathsf T}M=V\operatorname{diag}(\sigma_1(M)^2,\dots,\sigma_n(M)^2)V^{\mathsf T} M T M = V diag ( σ 1 ( M ) 2 , … , σ n ( M ) 2 ) V T であるから
( M T M ) − 1 = V diag ( σ 1 ( M ) − 2 , … , σ n ( M ) − 2 ) V T , ( M T M ) − 1 M T = V S U T (M^{\mathsf T}M)^{-1}=V\operatorname{diag}(\sigma_1(M)^{-2},\dots,\sigma_n(M)^{-2})V^{\mathsf T},\qquad(M^{\mathsf T}M)^{-1}M^{\mathsf T}=VSU^{\mathsf T} ( M T M ) − 1 = V diag ( σ 1 ( M ) − 2 , … , σ n ( M ) − 2 ) V T , ( M T M ) − 1 M T = V S U T である。ここでS ∈ R n × m S\in\R^{n\times m} S ∈ R n × m は( i , i ) (i,i) ( i , i ) 成分がσ i ( M ) − 1 \sigma_i(M)^{-1} σ i ( M ) − 1 (1 ≤ i ≤ n 1\le i\le n 1 ≤ i ≤ n )でその他の成分が0 0 0 の行列である。直交行列は Euclid ノルムを保つので、二つの行列の作用素ノルムは対角行列diag ( σ i ( M ) − 2 ) \operatorname{diag}(\sigma_i(M)^{-2}) diag ( σ i ( M ) − 2 ) とS S S の作用素ノルムに等しく、上と同じ計算によりそれぞれσ n ( M ) − 2 \sigma_n(M)^{-2} σ n ( M ) − 2 とσ n ( M ) − 1 \sigma_n(M)^{-1} σ n ( M ) − 1 である。これで(2) は示された。
E E E が(3) の仮定を満たすとする。任意のh ∈ R n h\in\R^n h ∈ R n に対して、(2) と§E20.2 補題 1.2 (1) により
∥ ( M + E ) h ∥ 2 ≥ ∥ M h ∥ 2 − ∥ E h ∥ 2 ≥ ( σ n ( M ) − ∥ E ∥ 2 ) ∥ h ∥ 2 \lVert(M+E)h\rVert_2\ge\lVert Mh\rVert_2-\lVert Eh\rVert_2\ge(\sigma_n(M)-\lVert E\rVert_2)\lVert h\rVert_2 ∥( M + E ) h ∥ 2 ≥ ∥ M h ∥ 2 − ∥ E h ∥ 2 ≥ ( σ n ( M ) − ∥ E ∥ 2 ) ∥ h ∥ 2 である。h ≠ 0 h\ne0 h = 0 ならば右辺は正であるから、M + E M+E M + E の列は一次独立である。(2) をM + E M+E M + E に適用し、等号を満たすh ≠ 0 h\ne0 h = 0 をとると、σ n ( M + E ) ∥ h ∥ 2 = ∥ ( M + E ) h ∥ 2 ≥ ( σ n ( M ) − ∥ E ∥ 2 ) ∥ h ∥ 2 \sigma_n(M+E)\lVert h\rVert_2=\lVert(M+E)h\rVert_2\ge(\sigma_n(M)-\lVert E\rVert_2)\lVert h\rVert_2 σ n ( M + E ) ∥ h ∥ 2 = ∥( M + E ) h ∥ 2 ≥ ( σ n ( M ) − ∥ E ∥ 2 ) ∥ h ∥ 2 である。▨
定義 2.2. m ≥ n m\ge n m ≥ n とし、A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n の列が一次独立であるとする。κ 2 ( A ) : = σ 1 ( A ) / σ n ( A ) \kappa_2(A):=\sigma_1(A)/\sigma_n(A) κ 2 ( A ) := σ 1 ( A ) / σ n ( A ) をA A A の 条件数 (condition number of a matrix with full column rank ) という。m = n m=n m = n のとき、§E20.5 命題 4.4 (4) により、この値は Euclid ノルムに関する正則行列の条件数κ 2 ( A ) \kappa_2(A) κ 2 ( A ) に一致する。
命題 2.3. m ≥ n m\ge n m ≥ n とし、A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n の列が一次独立であるとする。A T A A^{\mathsf T}A A T A は実対称正定値行列であり、その固有値はσ 1 ( A ) 2 , … , σ n ( A ) 2 \sigma_1(A)^2,\dots,\sigma_n(A)^2 σ 1 ( A ) 2 , … , σ n ( A ) 2 であって、κ 2 ( A T A ) = κ 2 ( A ) 2 \kappa_2(A^{\mathsf T}A)=\kappa_2(A)^2 κ 2 ( A T A ) = κ 2 ( A ) 2 が成り立つ。
証明. ( A T A ) T = A T A (A^{\mathsf T}A)^{\mathsf T}=A^{\mathsf T}A ( A T A ) T = A T A であり、x ≠ 0 x\ne0 x = 0 ならば列の一次独立性からA x ≠ 0 Ax\ne0 A x = 0 であってx T A T A x = ∥ A x ∥ 2 2 > 0 x^{\mathsf T}A^{\mathsf T}Ax=\lVert Ax\rVert_2^2>0 x T A T A x = ∥ A x ∥ 2 2 > 0 であるから、A T A A^{\mathsf T}A A T A は実対称正定値行列である。§D3.18 命題 2.1 により、σ 1 ( A ) 2 ≥ ⋯ ≥ σ n ( A ) 2 \sigma_1(A)^2\ge\dots\ge\sigma_n(A)^2 σ 1 ( A ) 2 ≥ ⋯ ≥ σ n ( A ) 2 はA T A A^{\mathsf T}A A T A のn n n 個の固有値である。§E20.5 命題 4.4 (5) によりκ 2 ( A T A ) = σ 1 ( A ) 2 / σ n ( A ) 2 = κ 2 ( A ) 2 \kappa_2(A^{\mathsf T}A)=\sigma_1(A)^2/\sigma_n(A)^2=\kappa_2(A)^2 κ 2 ( A T A ) = σ 1 ( A ) 2 / σ n ( A ) 2 = κ 2 ( A ) 2 である。▨
3 Householder 反射と QR 法
定義 3.1. ℓ ∈ N ≥ 1 \ell\in\NN ℓ ∈ N ≥ 1 とし、v ∈ R ℓ v\in\R^\ell v ∈ R ℓ をv ≠ 0 v\ne0 v = 0 を満たすベクトルとする。
H v : = I ℓ − 2 v T v v v T H_v:=I_\ell-\frac{2}{v^{\mathsf T}v}vv^{\mathsf T} H v := I ℓ − v T v 2 v v T をv v v の定める Householder 反射 (Householder reflector ) という。
補題 3.2. ℓ ∈ N ≥ 1 \ell\in\NN ℓ ∈ N ≥ 1 とする。
v ∈ R ℓ ∖ { 0 } v\in\R^\ell\setminus\{0\} v ∈ R ℓ ∖ { 0 } に対して、H v H_v H v は対称な直交行列であってH v 2 = I ℓ H_v^2=I_\ell H v 2 = I ℓ を満たす。H v v = − v H_vv=-v H v v = − v であり、v v v に直交する任意のw ∈ R ℓ w\in\R^\ell w ∈ R ℓ に対してH v w = w H_vw=w H v w = w である。
x ∈ R ℓ ∖ { 0 } x\in\R^\ell\setminus\{0\} x ∈ R ℓ ∖ { 0 } とし、α : = ∥ x ∥ 2 \alpha:=\lVert x\rVert_2 α := ∥ x ∥ 2 、x 1 ≥ 0 x_1\ge0 x 1 ≥ 0 のときs : = 1 s:=1 s := 1 、x 1 < 0 x_1<0 x 1 < 0 のときs : = − 1 s:=-1 s := − 1 と置き、v : = x + s α e 1 v:=x+s\alpha e_1 v := x + s α e 1 とする。このときs v 1 = ∣ x 1 ∣ + α > 0 sv_1=|x_1|+\alpha>0 s v 1 = ∣ x 1 ∣ + α > 0 、v T v = 2 α s v 1 v^{\mathsf T}v=2\alpha sv_1 v T v = 2 α s v 1 であり、任意のy ∈ R ℓ y\in\R^\ell y ∈ R ℓ に対して
H v y = y − v T y α s v 1 v H_vy=y-\frac{v^{\mathsf T}y}{\alpha sv_1}v H v y = y − α s v 1 v T y v
が成り立ち、H v x = − s α e 1 H_vx=-s\alpha e_1 H v x = − s α e 1 である。
(2) の記号でv ′ : = x − s α e 1 v':=x-s\alpha e_1 v ′ := x − s α e 1 と置く。s v 1 ′ = − ( ∑ i = 2 ℓ x i 2 ) / ( ∣ x 1 ∣ + α ) sv'_1=-\bigl(\sum_{i=2}^\ell x_i^2\bigr)/(|x_1|+\alpha) s v 1 ′ = − ( ∑ i = 2 ℓ x i 2 ) / ( ∣ x 1 ∣ + α ) であり、v ′ ≠ 0 v'\ne0 v ′ = 0 ならばH v ′ x = s α e 1 H_{v'}x=s\alpha e_1 H v ′ x = s α e 1 である。
証明. H v T = H v H_v^{\mathsf T}=H_v H v T = H v である。v v T v v T = ( v T v ) v v T vv^{\mathsf T}vv^{\mathsf T}=(v^{\mathsf T}v)vv^{\mathsf T} v v T v v T = ( v T v ) v v T から
H v 2 = I ℓ − 4 v T v v v T + 4 ( v T v ) 2 ( v T v ) v v T = I ℓ H_v^2=I_\ell-\frac{4}{v^{\mathsf T}v}vv^{\mathsf T}+\frac{4}{(v^{\mathsf T}v)^2}(v^{\mathsf T}v)vv^{\mathsf T}=I_\ell H v 2 = I ℓ − v T v 4 v v T + ( v T v ) 2 4 ( v T v ) v v T = I ℓ であり、H v T H v = H v 2 = I ℓ H_v^{\mathsf T}H_v=H_v^2=I_\ell H v T H v = H v 2 = I ℓ である。H v v = v − 2 v = − v H_vv=v-2v=-v H v v = v − 2 v = − v であり、v T w = 0 v^{\mathsf T}w=0 v T w = 0 ならばH v w = w H_vw=w H v w = w である。
(2) の記号で、s x 1 = ∣ x 1 ∣ sx_1=|x_1| s x 1 = ∣ x 1 ∣ であるからs v 1 = ∣ x 1 ∣ + α sv_1=|x_1|+\alpha s v 1 = ∣ x 1 ∣ + α であり、x ≠ 0 x\ne0 x = 0 からα > 0 \alpha>0 α > 0 である。v T v = α 2 + 2 s α x 1 + α 2 = 2 α ( α + ∣ x 1 ∣ ) = 2 α s v 1 v^{\mathsf T}v=\alpha^2+2s\alpha x_1+\alpha^2=2\alpha(\alpha+|x_1|)=2\alpha sv_1 v T v = α 2 + 2 s α x 1 + α 2 = 2 α ( α + ∣ x 1 ∣ ) = 2 α s v 1 であるから、H v H_v H v の定義によりH v y = y − ( v T y / ( α s v 1 ) ) v H_vy=y-(v^{\mathsf T}y/(\alpha sv_1))v H v y = y − ( v T y / ( α s v 1 )) v である。v T x = α 2 + s α x 1 = α s v 1 v^{\mathsf T}x=\alpha^2+s\alpha x_1=\alpha sv_1 v T x = α 2 + s α x 1 = α s v 1 であるからH v x = x − v = − s α e 1 H_vx=x-v=-s\alpha e_1 H v x = x − v = − s α e 1 である。
s v 1 ′ = ∣ x 1 ∣ − α = ( x 1 2 − α 2 ) / ( ∣ x 1 ∣ + α ) = − ( ∑ i ≥ 2 x i 2 ) / ( ∣ x 1 ∣ + α ) sv'_1=|x_1|-\alpha=(x_1^2-\alpha^2)/(|x_1|+\alpha)=-\bigl(\sum_{i\ge2}x_i^2\bigr)/(|x_1|+\alpha) s v 1 ′ = ∣ x 1 ∣ − α = ( x 1 2 − α 2 ) / ( ∣ x 1 ∣ + α ) = − ( ∑ i ≥ 2 x i 2 ) / ( ∣ x 1 ∣ + α ) である。v ′ T x = α 2 − s α x 1 = α ( α − ∣ x 1 ∣ ) v'^{\mathsf T}x=\alpha^2-s\alpha x_1=\alpha(\alpha-|x_1|) v ′ T x = α 2 − s α x 1 = α ( α − ∣ x 1 ∣ ) 、v ′ T v ′ = 2 α 2 − 2 s α x 1 = 2 α ( α − ∣ x 1 ∣ ) v'^{\mathsf T}v'=2\alpha^2-2s\alpha x_1=2\alpha(\alpha-|x_1|) v ′ T v ′ = 2 α 2 − 2 s α x 1 = 2 α ( α − ∣ x 1 ∣ ) であり、v ′ ≠ 0 v'\ne0 v ′ = 0 ならばv ′ T v ′ > 0 v'^{\mathsf T}v'>0 v ′ T v ′ > 0 であるから、H v ′ x = x − v ′ = s α e 1 H_{v'}x=x-v'=s\alpha e_1 H v ′ x = x − v ′ = s α e 1 である。▨
例 3.3. F = F ( 2 , 53 , − 1022 , 1023 ) F=F(2,53,-1022,1023) F = F ( 2 , 53 , − 1022 , 1023 ) 、fl \operatorname{fl} fl を最近接偶数丸めとし、x : = ( 1 , 2 − 30 ) ∈ F 2 x:=(1,2^{-30})\in F^2 x := ( 1 , 2 − 30 ) ∈ F 2 とする。α = 1 + 2 − 60 \alpha=\sqrt{1+2^{-60}} α = 1 + 2 − 60 、s = 1 s=1 s = 1 である。fl ( x 1 x 1 ) = 1 \operatorname{fl}(x_1x_1)=1 fl ( x 1 x 1 ) = 1 、fl ( x 2 x 2 ) = 2 − 60 \operatorname{fl}(x_2x_2)=2^{-60} fl ( x 2 x 2 ) = 2 − 60 であり、1 + 2 − 60 1+2^{-60} 1 + 2 − 60 の最近接点は、1 1 1 との距離2 − 60 2^{-60} 2 − 60 が1 + 2 − 52 1+2^{-52} 1 + 2 − 52 との距離より小さいので1 1 1 だけである。したがって二乗和の計算値は1 1 1 、α \alpha α の計算値はα ^ = 1 \hat\alpha=1 α ^ = 1 である。
補題 3.2 (3) のv ′ v' v ′ の第1 1 1 成分はv 1 ′ = 1 − 1 + 2 − 60 = − 2 − 60 / ( 1 + 1 + 2 − 60 ) ∈ ( − 2 − 61 , − 2 − 62 ) v'_1=1-\sqrt{1+2^{-60}}=-2^{-60}/(1+\sqrt{1+2^{-60}})\in(-2^{-61},-2^{-62}) v 1 ′ = 1 − 1 + 2 − 60 = − 2 − 60 / ( 1 + 1 + 2 − 60 ) ∈ ( − 2 − 61 , − 2 − 62 ) であるが、計算値はfl ( 1 − α ^ ) = 0 \operatorname{fl}(1-\hat\alpha)=0 fl ( 1 − α ^ ) = 0 であり、第1 1 1 成分の相対誤差は1 1 1 である。計算値のベクトルv ^ ′ = ( 0 , 2 − 30 ) \hat v'=(0,2^{-30}) v ^ ′ = ( 0 , 2 − 30 ) の定める反射はH v ^ ′ = diag ( 1 , − 1 ) H_{\hat v'}=\operatorname{diag}(1,-1) H v ^ ′ = diag ( 1 , − 1 ) であり、H v ^ ′ x = ( 1 , − 2 − 30 ) H_{\hat v'}x=(1,-2^{-30}) H v ^ ′ x = ( 1 , − 2 − 30 ) の第2 2 2 成分は0 0 0 でない。
補題 3.2 (2) のv v v では、計算値はv ^ = ( fl ( 1 + α ^ ) , 2 − 30 ) = ( 2 , 2 − 30 ) \hat v=(\operatorname{fl}(1+\hat\alpha),2^{-30})=(2,2^{-30}) v ^ = ( fl ( 1 + α ^ ) , 2 − 30 ) = ( 2 , 2 − 30 ) であり、v 1 = 1 + 1 + 2 − 60 v_1=1+\sqrt{1+2^{-60}} v 1 = 1 + 1 + 2 − 60 に対するv ^ 1 \hat v_1 v ^ 1 の相対誤差は2 − 62 2^{-62} 2 − 62 より小さい。v ^ \hat v v ^ の定める反射について
H v ^ x = ( − 1 − 2 2 62 + 1 , − 2 − 30 2 62 + 1 ) H_{\hat v}x=\Bigl(-1-\frac{2}{2^{62}+1},\ -\frac{2^{-30}}{2^{62}+1}\Bigr) H v ^ x = ( − 1 − 2 62 + 1 2 , − 2 62 + 1 2 − 30 ) であり、第2 2 2 成分の絶対値は2 − 92 2^{-92} 2 − 92 より小さい。s v 1 = ∣ x 1 ∣ + α sv_1=|x_1|+\alpha s v 1 = ∣ x 1 ∣ + α は正の二数の和であり、s v 1 ′ = ∣ x 1 ∣ − α sv'_1=|x_1|-\alpha s v 1 ′ = ∣ x 1 ∣ − α は補題 3.2 (3) により絶対値が∑ i ≥ 2 x i 2 / ( ∣ x 1 ∣ + α ) \sum_{i\ge2}x_i^2/(|x_1|+\alpha) ∑ i ≥ 2 x i 2 / ( ∣ x 1 ∣ + α ) の量を二つの近い正の数の差として与える。
定義 3.4. m ≥ n ≥ 1 m\ge n\ge1 m ≥ n ≥ 1 を整数とし、A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n 、b ∈ R m b\in\R^m b ∈ R m とする。作業行列W = ( w i j ) ∈ R m × ( n + 1 ) W=(w_{ij})\in\R^{m\times(n+1)} W = ( w ij ) ∈ R m × ( n + 1 ) をW : = ( A b ) W:=(A\ b) W := ( A b ) で初期化し、k = 1 , … , n k=1,\dots,n k = 1 , … , n の順に次の段k k k を行う計算式を、A A A と右辺b b b の Householder QR 法 (Householder QR factorization ) という。段k k k ではℓ : = m − k + 1 \ell:=m-k+1 ℓ := m − k + 1 と置き、x : = ( w k k , w k + 1 , k , … , w m k ) ∈ R ℓ x:=(w_{kk},w_{k+1,k},\dots,w_{mk})\in\R^\ell x := ( w k k , w k + 1 , k , … , w mk ) ∈ R ℓ とする。
x = 0 x=0 x = 0 ならば、段k k k では何も行わない。
x ≠ 0 x\ne0 x = 0 ならば、q : = x 1 x 1 q:=x_1x_1 q := x 1 x 1 と置き、i = 2 , … , ℓ i=2,\dots,\ell i = 2 , … , ℓ の順に乗算x i x i x_ix_i x i x i の後に加算を行ってq q q をq + x i x i q+x_ix_i q + x i x i に置き換え、α : = q \alpha:=\sqrt q α := q と置く。x 1 ≥ 0 x_1\ge0 x 1 ≥ 0 のときs : = 1 s:=1 s := 1 、x 1 < 0 x_1<0 x 1 < 0 のときs : = − 1 s:=-1 s := − 1 と置き、v 1 : = x 1 + s α v_1:=x_1+s\alpha v 1 := x 1 + s α 、v i : = x i v_i:=x_i v i := x i (2 ≤ i ≤ ℓ 2\le i\le\ell 2 ≤ i ≤ ℓ )、p : = α ⋅ ( s v 1 ) p:=\alpha\cdot(sv_1) p := α ⋅ ( s v 1 ) と置く。s v 1 sv_1 s v 1 はv 1 v_1 v 1 の符号をs s s に従って反転した値であり、演算に数えない。
x ≠ 0 x\ne0 x = 0 ならば、w k k w_{kk} w k k を− s α -s\alpha − s α に、w i k w_{ik} w ik (k < i ≤ m k<i\le m k < i ≤ m )を0 0 0 に置き換える。
x ≠ 0 x\ne0 x = 0 ならば、j = k + 1 , … , n + 1 j=k+1,\dots,n+1 j = k + 1 , … , n + 1 の順に次を行う。y : = ( w k j , w k + 1 , j , … , w m j ) y:=(w_{kj},w_{k+1,j},\dots,w_{mj}) y := ( w k j , w k + 1 , j , … , w mj ) とし、ω : = v 1 y 1 \omega:=v_1y_1 ω := v 1 y 1 と置き、i = 2 , … , ℓ i=2,\dots,\ell i = 2 , … , ℓ の順に乗算v i y i v_iy_i v i y i の後に加算を行ってω \omega ω をω + v i y i \omega+v_iy_i ω + v i y i に置き換え、τ : = ω / p \tau:=\omega/p τ := ω / p と置く。i = 1 , … , ℓ i=1,\dots,\ell i = 1 , … , ℓ について、乗算τ v i \tau v_i τ v i の後に減算を行ってw k + i − 1 , j w_{k+i-1,j} w k + i − 1 , j をy i − τ v i y_i-\tau v_i y i − τ v i に置き換える。
段n n n の後のW W W の第1 1 1 行から第n n n 行、第1 1 1 列から第n n n 列の部分をR ∈ R n × n R\in\R^{n\times n} R ∈ R n × n 、第n + 1 n+1 n + 1 列をc ∈ R m c\in\R^m c ∈ R m とし、( R , c ) (R,c) ( R , c ) をこの計算式の出力という。
補題 3.5. m ≥ n ≥ 1 m\ge n\ge1 m ≥ n ≥ 1 とし、Q ∈ R m × m Q\in\R^{m\times m} Q ∈ R m × m を直交行列、T ∈ R n × n T\in\R^{n\times n} T ∈ R n × n を対角成分がすべて0 0 0 でない上三角行列、b ∈ R m b\in\R^m b ∈ R m とする。A : = Q ( T 0 ) A:=Q\begin{pmatrix}T\\0\end{pmatrix} A := Q ( T 0 ) と置き、Q T b Q^{\mathsf T}b Q T b の第1 1 1 成分から第n n n 成分をc 1 ∈ R n c_1\in\R^n c 1 ∈ R n 、第n + 1 n+1 n + 1 成分から第m m m 成分をc 2 ∈ R m − n c_2\in\R^{m-n} c 2 ∈ R m − n とする。このときA A A の列は一次独立であり、A x = b Ax=b A x = b のただ一つの最小二乗解x x x はT x = c 1 Tx=c_1 T x = c 1 のただ一つの解であって、∥ b − A x ∥ 2 = ∥ c 2 ∥ 2 \lVert b-Ax\rVert_2=\lVert c_2\rVert_2 ∥ b − A x ∥ 2 = ∥ c 2 ∥ 2 が成り立つ。
証明. Q Q Q の第1 1 1 列から第n n n 列からなる行列をQ 1 Q_1 Q 1 とするとA = Q 1 T A=Q_1T A = Q 1 T 、Q 1 T Q 1 = I n Q_1^{\mathsf T}Q_1=I_n Q 1 T Q 1 = I n である。A h = 0 Ah=0 A h = 0 ならばT h = Q 1 T A h = 0 Th=Q_1^{\mathsf T}Ah=0 T h = Q 1 T A h = 0 であり、T T T は対角成分がすべて0 0 0 でない上三角行列であるから正則であって、h = 0 h=0 h = 0 である。したがってA A A の列は一次独立である。任意のz ∈ R n z\in\R^n z ∈ R n についてQ T ( b − A z ) = ( c 1 − T z c 2 ) Q^{\mathsf T}(b-Az)=\begin{pmatrix}c_1-Tz\\c_2\end{pmatrix} Q T ( b − A z ) = ( c 1 − T z c 2 ) であり、Q Q Q は Euclid ノルムを保つので
∥ b − A z ∥ 2 2 = ∥ c 1 − T z ∥ 2 2 + ∥ c 2 ∥ 2 2 ≥ ∥ c 2 ∥ 2 2 \lVert b-Az\rVert_2^2=\lVert c_1-Tz\rVert_2^2+\lVert c_2\rVert_2^2\ge\lVert c_2\rVert_2^2 ∥ b − A z ∥ 2 2 = ∥ c 1 − T z ∥ 2 2 + ∥ c 2 ∥ 2 2 ≥ ∥ c 2 ∥ 2 2 である。等号はT z = c 1 Tz=c_1 T z = c 1 のときに限って成り立ち、T T T は正則であるからT z = c 1 Tz=c_1 T z = c 1 の解x x x はただ一つである。したがってx x x はA x = b Ax=b A x = b のただ一つの最小二乗解であり、∥ b − A x ∥ 2 = ∥ c 2 ∥ 2 \lVert b-Ax\rVert_2=\lVert c_2\rVert_2 ∥ b − A x ∥ 2 = ∥ c 2 ∥ 2 である。▨
命題 3.6. m ≥ n ≥ 1 m\ge n\ge1 m ≥ n ≥ 1 とし、A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n 、b ∈ R m b\in\R^m b ∈ R m とする。A A A と右辺b b b の Householder QR 法を厳密算術で実行し、その出力を( R , c ) (R,c) ( R , c ) とする。段k k k のx x x が0 0 0 ならばH k : = I m H_k:=I_m H k := I m と置き、0 0 0 でないならば、段k k k のv v v によりH k : = diag ( I k − 1 , H v ) H_k:=\operatorname{diag}(I_{k-1},H_v) H k := diag ( I k − 1 , H v ) と置く。Q : = H 1 H 2 ⋯ H n Q:=H_1H_2\cdots H_n Q := H 1 H 2 ⋯ H n は直交行列であり、R R R は上三角行列であって
A = Q ( R 0 ) , b = Q c A=Q\begin{pmatrix}R\\0\end{pmatrix},\qquad b=Qc A = Q ( R 0 ) , b = Q c が成り立つ。A A A の列が一次独立ならば、R R R の対角成分はすべて0 0 0 でなく、A x = b Ax=b A x = b のただ一つの最小二乗解はR x = ( c 1 , … , c n ) R x=(c_1,\dots,c_n) R x = ( c 1 , … , c n ) のただ一つの解である。
証明. 段k k k の前の作業行列をW ( k − 1 ) W^{(k-1)} W ( k − 1 ) 、後の作業行列をW ( k ) W^{(k)} W ( k ) とする。段k k k は第1 1 1 列から第k − 1 k-1 k − 1 列を変えず、第k k k 列の第k + 1 k+1 k + 1 行から第m m m 行を0 0 0 にする(定義 3.4 (1) の場合は既に0 0 0 である)から、k k k に関する帰納法により、W ( k ) W^{(k)} W ( k ) の第j j j 列(j ≤ k j\le k j ≤ k )の第j + 1 j+1 j + 1 行から第m m m 行は0 0 0 である。特にR R R は上三角行列である。
W ( k ) = H k W ( k − 1 ) W^{(k)}=H_kW^{(k-1)} W ( k ) = H k W ( k − 1 ) が成り立つ。実際、定義 3.4 (1) の場合はH k = I m H_k=I_m H k = I m である。x ≠ 0 x\ne0 x = 0 の場合、H k H_k H k は第1 1 1 行から第k − 1 k-1 k − 1 行を変えない。j < k j<k j < k ならば第j j j 列の第k k k 行から第m m m 行は0 0 0 であるから、補題 3.2 (1) によりH k H_k H k は第j j j 列を変えない。補題 3.2 (2) によりH v x = − s α e 1 H_vx=-s\alpha e_1 H v x = − s α e 1 であり、これは定義 3.4 (3) の書き込みに一致する。j > k j>k j > k ならば、p = α s v 1 p=\alpha sv_1 p = α s v 1 であるから定義 3.4 (4) の値はy − ( v T y / ( α s v 1 ) ) v = H v y y-(v^{\mathsf T}y/(\alpha sv_1))v=H_vy y − ( v T y / ( α s v 1 )) v = H v y である。
したがってW ( n ) = H n ⋯ H 1 ( A b ) W^{(n)}=H_n\cdots H_1(A\ b) W ( n ) = H n ⋯ H 1 ( A b ) である。補題 3.2 (1) により各H k H_k H k は対称な直交行列であるから、Q Q Q は直交行列であってQ T = H n ⋯ H 1 Q^{\mathsf T}=H_n\cdots H_1 Q T = H n ⋯ H 1 である。W ( n ) W^{(n)} W ( n ) の第1 1 1 列から第n n n 列は( R 0 ) \begin{pmatrix}R\\0\end{pmatrix} ( R 0 ) 、第n + 1 n+1 n + 1 列はc c c であるから、A = Q ( R 0 ) A=Q\begin{pmatrix}R\\0\end{pmatrix} A = Q ( R 0 ) 、b = Q c b=Qc b = Q c を得る。A A A の列が一次独立ならば、( R 0 ) = Q T A \begin{pmatrix}R\\0\end{pmatrix}=Q^{\mathsf T}A ( R 0 ) = Q T A の列も一次独立であり、R R R は正則であるから、上三角行列R R R の対角成分はすべて0 0 0 でない。最後の主張は補題 3.5 をT = R T=R T = R に適用したものである。▨
例 3.7. 四点( t i , y i ) = ( 0 , 1 ) , ( 3 , 6 ) , ( 4 , 3 ) , ( 7 , 5 ) (t_i,y_i)=(0,1),(3,6),(4,3),(7,5) ( t i , y i ) = ( 0 , 1 ) , ( 3 , 6 ) , ( 4 , 3 ) , ( 7 , 5 ) に直線y = x 1 + x 2 t y=x_1+x_2t y = x 1 + x 2 t を当てはめる。A A A の第i i i 行を( 1 , t i ) (1,t_i) ( 1 , t i ) 、b : = ( 1 , 6 , 3 , 5 ) b:=(1,6,3,5) b := ( 1 , 6 , 3 , 5 ) とすると、四点は一つの直線上にないのでA x = b Ax=b A x = b は解をもたない。
Householder QR 法を厳密算術で実行する。段1 1 1 ではx = ( 1 , 1 , 1 , 1 ) x=(1,1,1,1) x = ( 1 , 1 , 1 , 1 ) 、α = 2 \alpha=2 α = 2 、s = 1 s=1 s = 1 、v = ( 3 , 1 , 1 , 1 ) v=(3,1,1,1) v = ( 3 , 1 , 1 , 1 ) 、p = 6 p=6 p = 6 であり、第2 2 2 列ではω = 14 \omega=14 ω = 14 、τ = 7 / 3 \tau=7/3 τ = 7/3 、右辺ではω = 17 \omega=17 ω = 17 、τ = 17 / 6 \tau=17/6 τ = 17/6 であるから、段1 1 1 の後の作業行列は
( − 2 − 7 − 15 / 2 0 2 / 3 19 / 6 0 5 / 3 1 / 6 0 14 / 3 13 / 6 ) \begin{pmatrix}-2&-7&-15/2\\0&2/3&19/6\\0&5/3&1/6\\0&14/3&13/6\end{pmatrix} − 2 0 0 0 − 7 2/3 5/3 14/3 − 15/2 19/6 1/6 13/6 である。段2 2 2 ではx = ( 2 / 3 , 5 / 3 , 14 / 3 ) x=(2/3,5/3,14/3) x = ( 2/3 , 5/3 , 14/3 ) 、α = 5 \alpha=5 α = 5 、s = 1 s=1 s = 1 、v = ( 17 / 3 , 5 / 3 , 14 / 3 ) v=(17/3,5/3,14/3) v = ( 17/3 , 5/3 , 14/3 ) 、p = 85 / 3 p=85/3 p = 85/3 であり、右辺ではω = 85 / 3 \omega=85/3 ω = 85/3 、τ = 1 \tau=1 τ = 1 であるから、出力は
R = ( − 2 − 7 0 − 5 ) , c = ( − 15 2 , − 5 2 , − 3 2 , − 5 2 ) R=\begin{pmatrix}-2&-7\\0&-5\end{pmatrix},\qquad c=\Bigl(-\frac{15}2,-\frac52,-\frac32,-\frac52\Bigr) R = ( − 2 0 − 7 − 5 ) , c = ( − 2 15 , − 2 5 , − 2 3 , − 2 5 ) である。R R R の対角成分は負であり、D : = − I 2 D:=-I_2 D := − I 2 と置くと、§D3.17 定義 1.1 の意味の QR 分解の上三角因子はD R = ( 2 7 0 5 ) DR=\begin{pmatrix}2&7\\0&5\end{pmatrix} D R = ( 2 0 7 5 ) である。命題 3.6 により、最小二乗解はR x = ( − 15 / 2 , − 5 / 2 ) Rx=(-15/2,-5/2) R x = ( − 15/2 , − 5/2 ) の解x = ( 2 , 1 / 2 ) x=(2,1/2) x = ( 2 , 1/2 ) である。残差はb − A x = ( − 1 , 5 / 2 , − 1 , − 1 / 2 ) b-Ax=(-1,5/2,-1,-1/2) b − A x = ( − 1 , 5/2 , − 1 , − 1/2 ) であり、∥ b − A x ∥ 2 2 = 17 / 2 = ( 3 / 2 ) 2 + ( 5 / 2 ) 2 \lVert b-Ax\rVert_2^2=17/2=(3/2)^2+(5/2)^2 ∥ b − A x ∥ 2 2 = 17/2 = ( 3/2 ) 2 + ( 5/2 ) 2 はc c c の第3 3 3 、第4 4 4 成分の二乗和に等しい。
4 浮動小数点算術での Householder QR 法
定義 4.1. F F F を浮動小数点数系、fl \operatorname{fl} fl をF F F の最近接丸めとする。x ∈ F x\in F x ∈ F がx ≥ 0 x\ge0 x ≥ 0 かつx ≤ N max \sqrt x\le N_{\max} x ≤ N m a x を満たすとき、fl ( x ) \operatorname{fl}(\sqrt x) fl ( x ) を 正しく丸めた平方根 (correctly rounded square root ) の結果といい、実数x \sqrt x x をその演算の厳密な結果という。加算・減算・乗算・除算と平方根を定められた順序で有限回行って値を定める計算式について、入力がF F F の元であるとき、各演算x ∘ y x\circ y x ∘ y をfl ( x ∘ y ) \operatorname{fl}(x\circ y) fl ( x ∘ y ) に、各平方根x \sqrt x x をfl ( x ) \operatorname{fl}(\sqrt x) fl ( x ) に順に置き換えて計算することを、計算式をF F F とfl \operatorname{fl} fl の浮動小数点算術で実行するという。その実行が次の条件を満たすとき、実行は範囲条件を満たすという。
各乗算・各除算・各平方根の厳密な結果は、0 0 0 であるかF F F の正規範囲にある。
各加算と各減算の厳密な結果の絶対値はN max N_{\max} N m a x 以下である。
平方根を含まない計算式については、この範囲条件は§E20.5 定義 1.1 の範囲条件と一致する。
補題 4.2. 0 ≤ u < 1 0\le u<1 0 ≤ u < 1 を実数とし、実数k ≥ 0 k\ge0 k ≥ 0 に対してB k : = [ ( 1 − u ) k , ( 1 − u ) − k ] B_k:=[(1-u)^k,(1-u)^{-k}] B k := [( 1 − u ) k , ( 1 − u ) − k ] と置く。
∣ δ ∣ ≤ u |\delta|\le u ∣ δ ∣ ≤ u を満たす実数δ \delta δ に対して1 + δ ∈ B 1 1+\delta\in B_1 1 + δ ∈ B 1 である。
実数j , k ≥ 0 j,k\ge0 j , k ≥ 0 とa ∈ B j a\in B_j a ∈ B j 、a ′ ∈ B k a'\in B_k a ′ ∈ B k に対して、a a ′ ∈ B j + k aa'\in B_{j+k} a a ′ ∈ B j + k 、a − 1 ∈ B j a^{-1}\in B_j a − 1 ∈ B j 、a ∈ B j / 2 \sqrt a\in B_{j/2} a ∈ B j /2 である。j ≤ k j\le k j ≤ k ならばB j ⊂ B k B_j\subset B_k B j ⊂ B k である。
a 1 , … , a N ∈ B k a_1,\dots,a_N\in B_k a 1 , … , a N ∈ B k とし、実数λ 1 , … , λ N ≥ 0 \lambda_1,\dots,\lambda_N\ge0 λ 1 , … , λ N ≥ 0 が∑ i λ i > 0 \sum_i\lambda_i>0 ∑ i λ i > 0 を満たすならば、∑ i λ i a i / ∑ i λ i ∈ B k \sum_i\lambda_ia_i/\sum_i\lambda_i\in B_k ∑ i λ i a i / ∑ i λ i ∈ B k である。
整数k ≥ 1 k\ge1 k ≥ 1 がk u < 1 ku<1 k u < 1 を満たすとし、γ k : = k u / ( 1 − k u ) \gamma_k:=ku/(1-ku) γ k := k u / ( 1 − k u ) と置く。任意のa ∈ B k a\in B_k a ∈ B k に対して∣ a − 1 ∣ ≤ γ k |a-1|\le\gamma_k ∣ a − 1∣ ≤ γ k である。
補題 4.3. F F F を浮動小数点数系、u u u をその単位丸め誤差、fl \operatorname{fl} fl をF F F の最近接丸めとし、実数k ≥ 0 k\ge0 k ≥ 0 に対してB k : = [ ( 1 − u ) k , ( 1 − u ) − k ] B_k:=[(1-u)^k,(1-u)^{-k}] B k := [( 1 − u ) k , ( 1 − u ) − k ] と置く。ℓ ∈ N ≥ 1 \ell\in\NN ℓ ∈ N ≥ 1 、x ∈ F ℓ x\in F^\ell x ∈ F ℓ 、x ≠ 0 x\ne0 x = 0 とし、α \alpha α 、s s s 、v v v を補題 3.2 (2) のとおりに置く。x x x から定義 3.4 (2) の計算式をF F F とfl \operatorname{fl} fl の浮動小数点算術で実行すると範囲条件を満たすとし、α \alpha α 、v v v 、p p p の計算値をα ^ \hat\alpha α ^ 、v ^ \hat v v ^ 、p ^ \hat p p ^ とする。このときc 1 ∈ B ℓ / 2 + 1 c_1\in B_{\ell/2+1} c 1 ∈ B ℓ /2 + 1 、c 2 ∈ B ℓ / 2 + 2 c_2\in B_{\ell/2+2} c 2 ∈ B ℓ /2 + 2 と∣ δ p ∣ ≤ u |\delta_p|\le u ∣ δ p ∣ ≤ u を満たす実数δ p \delta_p δ p が存在して
α ^ = c 1 α , v ^ = v + ( c 2 − 1 ) v 1 e 1 , p ^ = c 1 c 2 ( 1 + δ p ) α s v 1 \hat\alpha=c_1\alpha,\qquad\hat v=v+(c_2-1)v_1e_1,\qquad\hat p=c_1c_2(1+\delta_p)\,\alpha sv_1 α ^ = c 1 α , v ^ = v + ( c 2 − 1 ) v 1 e 1 , p ^ = c 1 c 2 ( 1 + δ p ) α s v 1 が成り立つ。特にα ^ > 0 \hat\alpha>0 α ^ > 0 、s v ^ 1 > 0 s\hat v_1>0 s v ^ 1 > 0 、p ^ > 0 \hat p>0 p ^ > 0 である。
証明. 範囲条件により各積x i x i x_ix_i x i x i は0 0 0 であるか正規範囲にあるから、§E20.1 系 3.2 (1) または§E20.1 系 3.2 (2) により、∣ μ i ∣ ≤ u |\mu_i|\le u ∣ μ i ∣ ≤ u を満たすμ i \mu_i μ i が存在してfl ( x i x i ) = x i 2 ( 1 + μ i ) \operatorname{fl}(x_ix_i)=x_i^2(1+\mu_i) fl ( x i x i ) = x i 2 ( 1 + μ i ) である。q q q の計算値q ^ \hat q q ^ はこれらの逐次和であり、§E20.3 系 1.4 (1) により各項の深さはℓ − 1 \ell-1 ℓ − 1 以下であるから、§E20.3 定理 1.3 (1) によりq ^ = ∑ i x i 2 e i \hat q=\sum_ix_i^2e_i q ^ = ∑ i x i 2 e i であって、各e i e_i e i は1 + μ i 1+\mu_i 1 + μ i と高々ℓ − 1 \ell-1 ℓ − 1 個の1 + δ 1+\delta 1 + δ (∣ δ ∣ ≤ u |\delta|\le u ∣ δ ∣ ≤ u )の積である。u ≤ 1 / 2 u\le1/2 u ≤ 1/2 であるから、補題 4.2 (1) と補題 4.2 (2) によりe i ∈ B ℓ e_i\in B_\ell e i ∈ B ℓ であり、α 2 = ∑ i x i 2 > 0 \alpha^2=\sum_ix_i^2>0 α 2 = ∑ i x i 2 > 0 であるから、補題 4.2 (3) によりq ^ / α 2 ∈ B ℓ \hat q/\alpha^2\in B_\ell q ^ / α 2 ∈ B ℓ であり、特にq ^ > 0 \hat q>0 q ^ > 0 である。q ^ > 0 \sqrt{\hat q}>0 q ^ > 0 であり、範囲条件によりq ^ \sqrt{\hat q} q ^ は正規範囲にあるから、§E20.1 定理 2.2 (2) により∣ δ α ∣ ≤ u |\delta_\alpha|\le u ∣ δ α ∣ ≤ u を満たすδ α \delta_\alpha δ α が存在してα ^ = q ^ ( 1 + δ α ) \hat\alpha=\sqrt{\hat q}(1+\delta_\alpha) α ^ = q ^ ( 1 + δ α ) である。c 1 : = α ^ / α = ( q ^ / α 2 ) 1 / 2 ( 1 + δ α ) c_1:=\hat\alpha/\alpha=(\hat q/\alpha^2)^{1/2}(1+\delta_\alpha) c 1 := α ^ / α = ( q ^ / α 2 ) 1/2 ( 1 + δ α ) は、補題 4.2 (2) によりB ℓ / 2 + 1 B_{\ell/2+1} B ℓ /2 + 1 に属する。
補題 3.2 (2) によりs v 1 = ∣ x 1 ∣ + α sv_1=|x_1|+\alpha s v 1 = ∣ x 1 ∣ + α であり、x 1 + s α ^ = s ( ∣ x 1 ∣ + c 1 α ) = c ′ v 1 x_1+s\hat\alpha=s(|x_1|+c_1\alpha)=c'v_1 x 1 + s α ^ = s ( ∣ x 1 ∣ + c 1 α ) = c ′ v 1 である。ここでc ′ : = ( ∣ x 1 ∣ ⋅ 1 + α c 1 ) / ( ∣ x 1 ∣ + α ) c':=(|x_1|\cdot1+\alpha c_1)/(|x_1|+\alpha) c ′ := ( ∣ x 1 ∣ ⋅ 1 + α c 1 ) / ( ∣ x 1 ∣ + α ) は1 1 1 とc 1 c_1 c 1 の正の重みによる平均であるから、補題 4.2 (3) によりB ℓ / 2 + 1 B_{\ell/2+1} B ℓ /2 + 1 に属する。x 1 x_1 x 1 とs α ^ s\hat\alpha s α ^ はF F F の元であり、範囲条件によりその和の絶対値はN max N_{\max} N m a x 以下であるから、§E20.1 系 3.3 により∣ δ v ∣ ≤ u |\delta_v|\le u ∣ δ v ∣ ≤ u を満たすδ v \delta_v δ v が存在してv ^ 1 = c ′ ( 1 + δ v ) v 1 \hat v_1=c'(1+\delta_v)v_1 v ^ 1 = c ′ ( 1 + δ v ) v 1 である。c 2 : = c ′ ( 1 + δ v ) c_2:=c'(1+\delta_v) c 2 := c ′ ( 1 + δ v ) はB ℓ / 2 + 2 B_{\ell/2+2} B ℓ /2 + 2 に属し、i ≥ 2 i\ge2 i ≥ 2 ではv ^ i = x i = v i \hat v_i=x_i=v_i v ^ i = x i = v i であるから、v ^ = v + ( c 2 − 1 ) v 1 e 1 \hat v=v+(c_2-1)v_1e_1 v ^ = v + ( c 2 − 1 ) v 1 e 1 である。
s v ^ 1 = c 2 s v 1 > 0 s\hat v_1=c_2sv_1>0 s v ^ 1 = c 2 s v 1 > 0 であり、積α ^ ⋅ s v ^ 1 = c 1 c 2 α s v 1 \hat\alpha\cdot s\hat v_1=c_1c_2\alpha sv_1 α ^ ⋅ s v ^ 1 = c 1 c 2 α s v 1 は正であって、範囲条件により正規範囲にあるから、§E20.1 系 3.2 (2) により∣ δ p ∣ ≤ u |\delta_p|\le u ∣ δ p ∣ ≤ u を満たすδ p \delta_p δ p が存在してp ^ = c 1 c 2 ( 1 + δ p ) α s v 1 \hat p=c_1c_2(1+\delta_p)\alpha sv_1 p ^ = c 1 c 2 ( 1 + δ p ) α s v 1 である。c 1 , c 2 , 1 + δ p c_1,c_2,1+\delta_p c 1 , c 2 , 1 + δ p は正であるから、α ^ > 0 \hat\alpha>0 α ^ > 0 、p ^ > 0 \hat p>0 p ^ > 0 である。▨
補題 4.4. F F F を浮動小数点数系、u u u をその単位丸め誤差、fl \operatorname{fl} fl をF F F の最近接丸めとし、ℓ ∈ N ≥ 1 \ell\in\NN ℓ ∈ N ≥ 1 が( 2 ℓ + 6 ) u < 1 (2\ell+6)u<1 ( 2 ℓ + 6 ) u < 1 を満たすとする。t ℓ : = γ 2 ℓ + 6 = ( 2 ℓ + 6 ) u / ( 1 − ( 2 ℓ + 6 ) u ) t_\ell:=\gamma_{2\ell+6}=(2\ell+6)u/(1-(2\ell+6)u) t ℓ := γ 2 ℓ + 6 = ( 2 ℓ + 6 ) u / ( 1 − ( 2 ℓ + 6 ) u ) と置く。x ∈ F ℓ x\in F^\ell x ∈ F ℓ 、x ≠ 0 x\ne0 x = 0 、y ∈ F ℓ y\in F^\ell y ∈ F ℓ とし、α \alpha α 、s s s 、v v v を補題 3.2 (2) のとおりに置く。x x x から定義 3.4 (2) の計算式を、続けてy y y から定義 3.4 (4) の計算式(ω \omega ω 、τ \tau τ とy i − τ v i y_i-\tau v_i y i − τ v i の計算)をF F F とfl \operatorname{fl} fl の浮動小数点算術で実行すると範囲条件を満たすとし、α \alpha α の計算値をα ^ \hat\alpha α ^ 、y i − τ v i y_i-\tau v_i y i − τ v i の計算値を並べたベクトルをy ^ ′ \hat y' y ^ ′ とする。
∣ α ^ − α ∣ ≤ t ℓ α |\hat\alpha-\alpha|\le t_\ell\alpha ∣ α ^ − α ∣ ≤ t ℓ α であり、したがって∥ − s α ^ e 1 − H v x ∥ 2 ≤ t ℓ ∥ x ∥ 2 \lVert-s\hat\alpha e_1-H_vx\rVert_2\le t_\ell\lVert x\rVert_2 ∥ − s α ^ e 1 − H v x ∥ 2 ≤ t ℓ ∥ x ∥ 2 である。
∥ y ^ ′ − H v y ∥ 2 ≤ 3 t ℓ ∥ y ∥ 2 \lVert\hat y'-H_vy\rVert_2\le3t_\ell\lVert y\rVert_2 ∥ y ^ ′ − H v y ∥ 2 ≤ 3 t ℓ ∥ y ∥ 2 である。
証明. 実数k ≥ 0 k\ge0 k ≥ 0 に対してB k : = [ ( 1 − u ) k , ( 1 − u ) − k ] B_k:=[(1-u)^k,(1-u)^{-k}] B k := [( 1 − u ) k , ( 1 − u ) − k ] と置き、補題 4.3 のc 1 c_1 c 1 、c 2 c_2 c 2 、δ p \delta_p δ p をとる。c 1 ∈ B ℓ / 2 + 1 ⊂ B 2 ℓ + 6 c_1\in B_{\ell/2+1}\subset B_{2\ell+6} c 1 ∈ B ℓ /2 + 1 ⊂ B 2 ℓ + 6 であるから、補題 4.2 (4) により∣ α ^ − α ∣ = α ∣ c 1 − 1 ∣ ≤ t ℓ α |\hat\alpha-\alpha|=\alpha|c_1-1|\le t_\ell\alpha ∣ α ^ − α ∣ = α ∣ c 1 − 1∣ ≤ t ℓ α である。補題 3.2 (2) によりH v x = − s α e 1 H_vx=-s\alpha e_1 H v x = − s α e 1 であるから、(1) が成り立つ。
β : = 1 / ( α s v 1 ) \beta:=1/(\alpha sv_1) β := 1/ ( α s v 1 ) と置く。補題 3.2 (2) によりv T v = 2 α s v 1 v^{\mathsf T}v=2\alpha sv_1 v T v = 2 α s v 1 であるから、H v = I ℓ − β v v T H_v=I_\ell-\beta vv^{\mathsf T} H v = I ℓ − β v v T 、β ∥ v ∥ 2 2 = 2 \beta\lVert v\rVert_2^2=2 β ∥ v ∥ 2 2 = 2 である。ι 1 : = 1 \iota_1:=1 ι 1 := 1 、ι i : = 0 \iota_i:=0 ι i := 0 (i ≥ 2 i\ge2 i ≥ 2 )と置くと、補題 4.3 によりv ^ i = c 2 ι i v i \hat v_i=c_2^{\iota_i}v_i v ^ i = c 2 ι i v i 、1 / p ^ = β / ( c 1 c 2 ( 1 + δ p ) ) 1/\hat p=\beta/(c_1c_2(1+\delta_p)) 1/ p ^ = β / ( c 1 c 2 ( 1 + δ p )) である。範囲条件と§E20.1 系 3.2 (1) 、§E20.1 系 3.2 (2) によりfl ( v ^ i y i ) = v ^ i y i ( 1 + μ i ) \operatorname{fl}(\hat v_iy_i)=\hat v_iy_i(1+\mu_i) fl ( v ^ i y i ) = v ^ i y i ( 1 + μ i ) (∣ μ i ∣ ≤ u |\mu_i|\le u ∣ μ i ∣ ≤ u )であり、§E20.3 系 1.4 (1) と§E20.3 定理 1.3 (1) により、ω \omega ω の計算値はω ^ = ∑ j v ^ j y j d j \hat\omega=\sum_j\hat v_jy_jd_j ω ^ = ∑ j v ^ j y j d j であって、各d j d_j d j は1 + μ j 1+\mu_j 1 + μ j と高々ℓ − 1 \ell-1 ℓ − 1 個の1 + δ 1+\delta 1 + δ (∣ δ ∣ ≤ u |\delta|\le u ∣ δ ∣ ≤ u )の積である。同じく範囲条件と§E20.1 系 3.2 (1) 、§E20.1 系 3.2 (2) により、∣ δ τ ∣ , ∣ δ i ∣ ≤ u |\delta_\tau|,|\delta_i|\le u ∣ δ τ ∣ , ∣ δ i ∣ ≤ u を満たす実数が存在して、τ \tau τ の計算値はτ ^ = ( ω ^ / p ^ ) ( 1 + δ τ ) \hat\tau=(\hat\omega/\hat p)(1+\delta_\tau) τ ^ = ( ω ^ / p ^ ) ( 1 + δ τ ) 、積τ v i \tau v_i τ v i の計算値はz ^ i = τ ^ v ^ i ( 1 + δ i ) \hat z_i=\hat\tau\hat v_i(1+\delta_i) z ^ i = τ ^ v ^ i ( 1 + δ i ) である。これらを代入すると
z ^ i = β v i ∑ j = 1 ℓ v j y j g i j , g i j : = c 2 ι i + ι j − 1 d j ( 1 + δ τ ) ( 1 + δ i ) c 1 ( 1 + δ p ) \hat z_i=\beta v_i\sum_{j=1}^\ell v_jy_jg_{ij},\qquad g_{ij}:=\frac{c_2^{\iota_i+\iota_j-1}d_j(1+\delta_\tau)(1+\delta_i)}{c_1(1+\delta_p)} z ^ i = β v i j = 1 ∑ ℓ v j y j g ij , g ij := c 1 ( 1 + δ p ) c 2 ι i + ι j − 1 d j ( 1 + δ τ ) ( 1 + δ i ) である。ι i + ι j − 1 ∈ { − 1 , 0 , 1 } \iota_i+\iota_j-1\in\{-1,0,1\} ι i + ι j − 1 ∈ { − 1 , 0 , 1 } であり、補題 4.2 (1) と補題 4.2 (2) によりg i j ∈ B ( ℓ / 2 + 2 ) + ℓ + 2 + ( ℓ / 2 + 1 ) + 1 = B 2 ℓ + 6 g_{ij}\in B_{(\ell/2+2)+\ell+2+(\ell/2+1)+1}=B_{2\ell+6} g ij ∈ B ( ℓ /2 + 2 ) + ℓ + 2 + ( ℓ /2 + 1 ) + 1 = B 2 ℓ + 6 であるから、補題 4.2 (4) により∣ g i j − 1 ∣ ≤ t ℓ |g_{ij}-1|\le t_\ell ∣ g ij − 1∣ ≤ t ℓ である。
z : = β v ( v T y ) z:=\beta v(v^{\mathsf T}y) z := β v ( v T y ) と置くとH v y = y − z H_vy=y-z H v y = y − z である。Cauchy–Schwarz の不等式により
∣ z ^ i − z i ∣ = β ∣ v i ∣ ∣ ∑ j v j y j ( g i j − 1 ) ∣ ≤ t ℓ β ∣ v i ∣ ∑ j ∣ v j ∣ ∣ y j ∣ ≤ t ℓ β ∣ v i ∣ ∥ v ∥ 2 ∥ y ∥ 2 |\hat z_i-z_i|=\beta|v_i|\Bigl|\sum_jv_jy_j(g_{ij}-1)\Bigr|\le t_\ell\beta|v_i|\sum_j|v_j||y_j|\le t_\ell\beta|v_i|\lVert v\rVert_2\lVert y\rVert_2 ∣ z ^ i − z i ∣ = β ∣ v i ∣ j ∑ v j y j ( g ij − 1 ) ≤ t ℓ β ∣ v i ∣ j ∑ ∣ v j ∣∣ y j ∣ ≤ t ℓ β ∣ v i ∣ ∥ v ∥ 2 ∥ y ∥ 2 であり、∥ z ^ − z ∥ 2 ≤ t ℓ β ∥ v ∥ 2 2 ∥ y ∥ 2 = 2 t ℓ ∥ y ∥ 2 \lVert\hat z-z\rVert_2\le t_\ell\beta\lVert v\rVert_2^2\lVert y\rVert_2=2t_\ell\lVert y\rVert_2 ∥ z ^ − z ∥ 2 ≤ t ℓ β ∥ v ∥ 2 2 ∥ y ∥ 2 = 2 t ℓ ∥ y ∥ 2 である。y i y_i y i と− z ^ i -\hat z_i − z ^ i はF F F の元であり、範囲条件により§E20.1 系 3.3 が適用されるので、∣ δ i ′ ∣ ≤ u |\delta'_i|\le u ∣ δ i ′ ∣ ≤ u を満たす実数が存在してy ^ i ′ = ( y i − z ^ i ) ( 1 + δ i ′ ) \hat y'_i=(y_i-\hat z_i)(1+\delta'_i) y ^ i ′ = ( y i − z ^ i ) ( 1 + δ i ′ ) である。したがってy ^ ′ − H v y = ( z − z ^ ) + ( δ i ′ ( y i − z ^ i ) ) i \hat y'-H_vy=(z-\hat z)+(\delta'_i(y_i-\hat z_i))_i y ^ ′ − H v y = ( z − z ^ ) + ( δ i ′ ( y i − z ^ i ) ) i である。補題 3.2 (1) により∥ H v y ∥ 2 = ∥ y ∥ 2 \lVert H_vy\rVert_2=\lVert y\rVert_2 ∥ H v y ∥ 2 = ∥ y ∥ 2 であるから、∥ y − z ^ ∥ 2 ≤ ∥ y − z ∥ 2 + ∥ z − z ^ ∥ 2 ≤ ( 1 + 2 t ℓ ) ∥ y ∥ 2 \lVert y-\hat z\rVert_2\le\lVert y-z\rVert_2+\lVert z-\hat z\rVert_2\le(1+2t_\ell)\lVert y\rVert_2 ∥ y − z ^ ∥ 2 ≤ ∥ y − z ∥ 2 + ∥ z − z ^ ∥ 2 ≤ ( 1 + 2 t ℓ ) ∥ y ∥ 2 であり、
∥ y ^ ′ − H v y ∥ 2 ≤ ( 2 t ℓ + u ( 1 + 2 t ℓ ) ) ∥ y ∥ 2 \lVert\hat y'-H_vy\rVert_2\le\bigl(2t_\ell+u(1+2t_\ell)\bigr)\lVert y\rVert_2 ∥ y ^ ′ − H v y ∥ 2 ≤ ( 2 t ℓ + u ( 1 + 2 t ℓ ) ) ∥ y ∥ 2 である。t ℓ ≥ ( 2 ℓ + 6 ) u ≥ 8 u t_\ell\ge(2\ell+6)u\ge8u t ℓ ≥ ( 2 ℓ + 6 ) u ≥ 8 u であり、( 2 ℓ + 6 ) u < 1 (2\ell+6)u<1 ( 2 ℓ + 6 ) u < 1 からu < 1 / 8 u<1/8 u < 1/8 であるので、t ℓ ( 1 − 2 u ) ≥ 8 u ⋅ 3 / 4 ≥ u t_\ell(1-2u)\ge8u\cdot3/4\ge u t ℓ ( 1 − 2 u ) ≥ 8 u ⋅ 3/4 ≥ u 、すなわちu ( 1 + 2 t ℓ ) ≤ t ℓ u(1+2t_\ell)\le t_\ell u ( 1 + 2 t ℓ ) ≤ t ℓ である。これで(2) は示された。▨
定理 4.5. F F F を浮動小数点数系、u u u をその単位丸め誤差、fl \operatorname{fl} fl をF F F の最近接丸めとし、整数m ≥ n ≥ 1 m\ge n\ge1 m ≥ n ≥ 1 が( 2 m + 6 ) u < 1 (2m+6)u<1 ( 2 m + 6 ) u < 1 を満たすとする。整数j ≥ 0 j\ge0 j ≥ 0 がj u < 1 ju<1 j u < 1 を満たすときγ j : = j u / ( 1 − j u ) \gamma_j:=ju/(1-ju) γ j := j u / ( 1 − j u ) と置き、
t : = γ 2 m + 6 , η m , n : = ( 1 + 3 t ) n − 1 t:=\gamma_{2m+6},\qquad\eta_{m,n}:=(1+3t)^n-1 t := γ 2 m + 6 , η m , n := ( 1 + 3 t ) n − 1 と置く。A ∈ F m × n A\in F^{m\times n} A ∈ F m × n 、b ∈ F m b\in F^m b ∈ F m とし、A A A の第j j j 列をa j a_j a j と書く。A A A と右辺b b b の Householder QR 法をF F F とfl \operatorname{fl} fl の浮動小数点算術で実行すると範囲条件を満たすとし、その出力を( R ^ , c ^ ) (\widehat R,\widehat c) ( R , c ) とする。
R ^ \widehat R R は上三角行列である。各k k k についてI m I_m I m であるかdiag ( I k − 1 , H v ) \operatorname{diag}(I_{k-1},H_v) diag ( I k − 1 , H v ) (H v H_v H v は Householder 反射)の形である直交行列H k H_k H k の積Q : = H 1 ⋯ H n ∈ R m × m Q:=H_1\cdots H_n\in\R^{m\times m} Q := H 1 ⋯ H n ∈ R m × m と、Δ A ∈ R m × n \Delta A\in\R^{m\times n} Δ A ∈ R m × n 、Δ b ∈ R m \Delta b\in\R^m Δ b ∈ R m が存在して
A + Δ A = Q ( R ^ 0 ) , b + Δ b = Q c ^ A+\Delta A=Q\begin{pmatrix}\widehat R\\0\end{pmatrix},\qquad b+\Delta b=Q\widehat c A + Δ A = Q ( R 0 ) , b + Δ b = Q c
が成り立ち、Δ A \Delta A Δ A の第j j j 列Δ a j \Delta a_j Δ a j とΔ b \Delta b Δ b は
∥ Δ a j ∥ 2 ≤ η m , n ∥ a j ∥ 2 ( 1 ≤ j ≤ n ) , ∥ Δ b ∥ 2 ≤ η m , n ∥ b ∥ 2 \lVert\Delta a_j\rVert_2\le\eta_{m,n}\lVert a_j\rVert_2\quad(1\le j\le n),\qquad\lVert\Delta b\rVert_2\le\eta_{m,n}\lVert b\rVert_2 ∥ Δ a j ∥ 2 ≤ η m , n ∥ a j ∥ 2 ( 1 ≤ j ≤ n ) , ∥ Δ b ∥ 2 ≤ η m , n ∥ b ∥ 2
を満たす。特に∥ Δ A ∥ F ≤ η m , n ∥ A ∥ F \lVert\Delta A\rVert_F\le\eta_{m,n}\lVert A\rVert_F ∥ Δ A ∥ F ≤ η m , n ∥ A ∥ F である。
3 n t ≤ 1 3nt\le1 3 n t ≤ 1 ならばη m , n ≤ 6 n t \eta_{m,n}\le6nt η m , n ≤ 6 n t である。
証明. W ( 0 ) : = ( A b ) W^{(0)}:=(A\ b) W ( 0 ) := ( A b ) と置き、1 ≤ k ≤ n 1\le k\le n 1 ≤ k ≤ n に対して段k k k の後の作業行列の計算値をW ( k ) ∈ F m × ( n + 1 ) W^{(k)}\in F^{m\times(n+1)} W ( k ) ∈ F m × ( n + 1 ) とする。W ( k − 1 ) W^{(k-1)} W ( k − 1 ) の第k k k 列の第k k k 行から第m m m 行をx ( k ) ∈ F m − k + 1 x^{(k)}\in F^{m-k+1} x ( k ) ∈ F m − k + 1 とする。x ( k ) = 0 x^{(k)}=0 x ( k ) = 0 ならばH k : = I m H_k:=I_m H k := I m と置き、x ( k ) ≠ 0 x^{(k)}\ne0 x ( k ) = 0 ならば、x ( k ) x^{(k)} x ( k ) に対して補題 3.2 (2) のv v v をv ( k ) v^{(k)} v ( k ) とし、H k : = diag ( I k − 1 , H v ( k ) ) H_k:=\operatorname{diag}(I_{k-1},H_{v^{(k)}}) H k := diag ( I k − 1 , H v ( k ) ) と置く。補題 3.2 (1) によりH k H_k H k は対称な直交行列である。E k : = W ( k ) − H k W ( k − 1 ) E_k:=W^{(k)}-H_kW^{(k-1)} E k := W ( k ) − H k W ( k − 1 ) と置く。
段k k k は第1 1 1 列から第k − 1 k-1 k − 1 列を変えず、第k k k 列の第k + 1 k+1 k + 1 行から第m m m 行を0 0 0 にする(定義 3.4 (1) の場合は既に0 0 0 である)から、k k k に関する帰納法により、W ( k ) W^{(k)} W ( k ) の第j j j 列(j ≤ k j\le k j ≤ k )の第j + 1 j+1 j + 1 行から第m m m 行は0 0 0 である。
主張 4.5.1. 1 ≤ k ≤ n 1\le k\le n 1 ≤ k ≤ n 、1 ≤ j ≤ n + 1 1\le j\le n+1 1 ≤ j ≤ n + 1 ならば、∥ E k e j ∥ 2 ≤ 3 t ∥ W ( k − 1 ) e j ∥ 2 \lVert E_ke_j\rVert_2\le3t\lVert W^{(k-1)}e_j\rVert_2 ∥ E k e j ∥ 2 ≤ 3 t ∥ W ( k − 1 ) e j ∥ 2 である。
証明. x ( k ) = 0 x^{(k)}=0 x ( k ) = 0 ならば段k k k は作業行列を変えず、H k = I m H_k=I_m H k = I m であるからE k = 0 E_k=0 E k = 0 である。x ( k ) ≠ 0 x^{(k)}\ne0 x ( k ) = 0 とし、ℓ : = m − k + 1 \ell:=m-k+1 ℓ := m − k + 1 と置く。ℓ ≤ m \ell\le m ℓ ≤ m であるから( 2 ℓ + 6 ) u < 1 (2\ell+6)u<1 ( 2 ℓ + 6 ) u < 1 であり、γ j \gamma_j γ j はj j j について単調非減少であるから、補題 4.4 のt ℓ t_\ell t ℓ はt t t 以下である。段k k k とH k H_k H k はどちらも第1 1 1 行から第k − 1 k-1 k − 1 行を変えないので、E k e j E_ke_j E k e j の第1 1 1 成分から第k − 1 k-1 k − 1 成分は0 0 0 である。j < k j<k j < k ならば、W ( k − 1 ) e j W^{(k-1)}e_j W ( k − 1 ) e j の第k k k 成分から第m m m 成分は0 0 0 であるから、補題 3.2 (1) によりH k W ( k − 1 ) e j = W ( k − 1 ) e j H_kW^{(k-1)}e_j=W^{(k-1)}e_j H k W ( k − 1 ) e j = W ( k − 1 ) e j であり、段k k k は第j j j 列を変えないのでE k e j = 0 E_ke_j=0 E k e j = 0 である。j = k j=k j = k ならば、E k e j E_ke_j E k e j の第k k k 成分から第m m m 成分は− s α ^ e 1 − H v ( k ) x ( k ) -s\hat\alpha e_1-H_{v^{(k)}}x^{(k)} − s α ^ e 1 − H v ( k ) x ( k ) であり、補題 4.4 (1) によりそのノルムはt ∥ x ( k ) ∥ 2 ≤ t ∥ W ( k − 1 ) e k ∥ 2 t\lVert x^{(k)}\rVert_2\le t\lVert W^{(k-1)}e_k\rVert_2 t ∥ x ( k ) ∥ 2 ≤ t ∥ W ( k − 1 ) e k ∥ 2 以下である。j > k j>k j > k ならば、W ( k − 1 ) e j W^{(k-1)}e_j W ( k − 1 ) e j の第k k k 成分から第m m m 成分をy y y とすると、E k e j E_ke_j E k e j の第k k k 成分から第m m m 成分はy ^ ′ − H v ( k ) y \hat y'-H_{v^{(k)}}y y ^ ′ − H v ( k ) y であり、補題 4.4 (2) によりそのノルムは3 t ∥ y ∥ 2 ≤ 3 t ∥ W ( k − 1 ) e j ∥ 2 3t\lVert y\rVert_2\le3t\lVert W^{(k-1)}e_j\rVert_2 3 t ∥ y ∥ 2 ≤ 3 t ∥ W ( k − 1 ) e j ∥ 2 以下である。▨
H k H_k H k は Euclid ノルムを保つから、主張 4.5.1 により∥ W ( k ) e j ∥ 2 ≤ ∥ H k W ( k − 1 ) e j ∥ 2 + ∥ E k e j ∥ 2 ≤ ( 1 + 3 t ) ∥ W ( k − 1 ) e j ∥ 2 \lVert W^{(k)}e_j\rVert_2\le\lVert H_kW^{(k-1)}e_j\rVert_2+\lVert E_ke_j\rVert_2\le(1+3t)\lVert W^{(k-1)}e_j\rVert_2 ∥ W ( k ) e j ∥ 2 ≤ ∥ H k W ( k − 1 ) e j ∥ 2 + ∥ E k e j ∥ 2 ≤ ( 1 + 3 t ) ∥ W ( k − 1 ) e j ∥ 2 であり、k k k に関する帰納法により∥ W ( k − 1 ) e j ∥ 2 ≤ ( 1 + 3 t ) k − 1 ∥ W ( 0 ) e j ∥ 2 \lVert W^{(k-1)}e_j\rVert_2\le(1+3t)^{k-1}\lVert W^{(0)}e_j\rVert_2 ∥ W ( k − 1 ) e j ∥ 2 ≤ ( 1 + 3 t ) k − 1 ∥ W ( 0 ) e j ∥ 2 である。W ( k ) = H k W ( k − 1 ) + E k W^{(k)}=H_kW^{(k-1)}+E_k W ( k ) = H k W ( k − 1 ) + E k をk = 1 , … , n k=1,\dots,n k = 1 , … , n について順に代入すると
W ( n ) = H n ⋯ H 1 W ( 0 ) + ∑ k = 1 n H n ⋯ H k + 1 E k W^{(n)}=H_n\cdots H_1W^{(0)}+\sum_{k=1}^nH_n\cdots H_{k+1}E_k W ( n ) = H n ⋯ H 1 W ( 0 ) + k = 1 ∑ n H n ⋯ H k + 1 E k である。Q : = H 1 H 2 ⋯ H n Q:=H_1H_2\cdots H_n Q := H 1 H 2 ⋯ H n は直交行列であり、H k 2 = I m H_k^2=I_m H k 2 = I m からQ H n ⋯ H k + 1 = H 1 ⋯ H k QH_n\cdots H_{k+1}=H_1\cdots H_k Q H n ⋯ H k + 1 = H 1 ⋯ H k であるので、Δ W : = ∑ k = 1 n H 1 ⋯ H k E k \Delta W:=\sum_{k=1}^nH_1\cdots H_kE_k Δ W := ∑ k = 1 n H 1 ⋯ H k E k と置くとQ W ( n ) = W ( 0 ) + Δ W QW^{(n)}=W^{(0)}+\Delta W Q W ( n ) = W ( 0 ) + Δ W である。直交行列は Euclid ノルムを保つから、
∥ Δ W e j ∥ 2 ≤ ∑ k = 1 n ∥ E k e j ∥ 2 ≤ ∑ k = 1 n 3 t ( 1 + 3 t ) k − 1 ∥ W ( 0 ) e j ∥ 2 = η m , n ∥ W ( 0 ) e j ∥ 2 \lVert\Delta We_j\rVert_2\le\sum_{k=1}^n\lVert E_ke_j\rVert_2\le\sum_{k=1}^n3t(1+3t)^{k-1}\lVert W^{(0)}e_j\rVert_2=\eta_{m,n}\lVert W^{(0)}e_j\rVert_2 ∥ Δ W e j ∥ 2 ≤ k = 1 ∑ n ∥ E k e j ∥ 2 ≤ k = 1 ∑ n 3 t ( 1 + 3 t ) k − 1 ∥ W ( 0 ) e j ∥ 2 = η m , n ∥ W ( 0 ) e j ∥ 2 である。j ≤ n j\le n j ≤ n ならばW ( n ) W^{(n)} W ( n ) の第j j j 列の第j + 1 j+1 j + 1 行から第m m m 行は0 0 0 であるから、第1 1 1 列から第n n n 列は( R ^ 0 ) \begin{pmatrix}\widehat R\\0\end{pmatrix} ( R 0 ) でR ^ \widehat R R は上三角行列であり、第n + 1 n+1 n + 1 列はc ^ \widehat c c である。Δ W \Delta W Δ W の第1 1 1 列から第n n n 列をΔ A \Delta A Δ A 、第n + 1 n+1 n + 1 列をΔ b \Delta b Δ b とすると、A + Δ A = Q ( R ^ 0 ) A+\Delta A=Q\begin{pmatrix}\widehat R\\0\end{pmatrix} A + Δ A = Q ( R 0 ) 、b + Δ b = Q c ^ b+\Delta b=Q\widehat c b + Δ b = Q c と列ごとの評価が成り立ち、∥ Δ A ∥ F 2 = ∑ j ∥ Δ a j ∥ 2 2 ≤ η m , n 2 ∥ A ∥ F 2 \lVert\Delta A\rVert_F^2=\sum_j\lVert\Delta a_j\rVert_2^2\le\eta_{m,n}^2\lVert A\rVert_F^2 ∥ Δ A ∥ F 2 = ∑ j ∥ Δ a j ∥ 2 2 ≤ η m , n 2 ∥ A ∥ F 2 である。これで(1) は示された。
3 n t ≤ 1 3nt\le1 3 n t ≤ 1 とする。1 + 3 t ≤ e 3 t 1+3t\le e^{3t} 1 + 3 t ≤ e 3 t からη m , n ≤ e 3 n t − 1 \eta_{m,n}\le e^{3nt}-1 η m , n ≤ e 3 n t − 1 である。0 ≤ y ≤ 1 0\le y\le1 0 ≤ y ≤ 1 に対して、e y e^y e y の凸性からe y ≤ ( 1 − y ) + y e e^y\le(1-y)+ye e y ≤ ( 1 − y ) + y e であり、e y − 1 ≤ ( e − 1 ) y ≤ 2 y e^y-1\le(e-1)y\le2y e y − 1 ≤ ( e − 1 ) y ≤ 2 y である。y = 3 n t y=3nt y = 3 n t としてη m , n ≤ 6 n t \eta_{m,n}\le6nt η m , n ≤ 6 n t を得る。▨
定理 4.7. F F F を浮動小数点数系、u u u をその単位丸め誤差、fl \operatorname{fl} fl をF F F の最近接丸めとし、整数m ≥ n ≥ 1 m\ge n\ge1 m ≥ n ≥ 1 が( 2 m + 6 ) u < 1 (2m+6)u<1 ( 2 m + 6 ) u < 1 を満たすとする。γ j \gamma_j γ j 、t t t 、η m , n \eta_{m,n} η m , n を定理 4.5 のとおりに置き、η m , n ′ : = η m , n + γ n ( 1 + η m , n ) \eta'_{m,n}:=\eta_{m,n}+\gamma_n(1+\eta_{m,n}) η m , n ′ := η m , n + γ n ( 1 + η m , n ) と置く。A ∈ F m × n A\in F^{m\times n} A ∈ F m × n の列が一次独立であるとし、b ∈ F m b\in F^m b ∈ F m とする。A A A と右辺b b b の Householder QR 法をF F F とfl \operatorname{fl} fl の浮動小数点算術で実行すると範囲条件を満たすとし、その出力を( R ^ , c ^ ) (\widehat R,\widehat c) ( R , c ) 、c ^ \widehat c c の第1 1 1 成分から第n n n 成分をc ^ 1 \widehat c_1 c 1 とする。η m , n ∥ A ∥ F < σ n ( A ) \eta_{m,n}\lVert A\rVert_F<\sigma_n(A) η m , n ∥ A ∥ F < σ n ( A ) を仮定する。
R ^ \widehat R R は対角成分がすべて0 0 0 でない上三角行列である。
R ^ x = c ^ 1 \widehat Rx=\widehat c_1 R x = c 1 の後退代入をF F F とfl \operatorname{fl} fl の浮動小数点算術で実行すると範囲条件を満たすとし、その計算値をx ^ \hat x x ^ とする。A A A の第j j j 列をa j a_j a j と書く。このとき、Δ A ′ ∈ R m × n \Delta A'\in\R^{m\times n} Δ A ′ ∈ R m × n とΔ b ∈ R m \Delta b\in\R^m Δ b ∈ R m が存在して、Δ A ′ \Delta A' Δ A ′ の第j j j 列Δ a j ′ \Delta a'_j Δ a j ′ が∥ Δ a j ′ ∥ 2 ≤ η m , n ′ ∥ a j ∥ 2 \lVert\Delta a'_j\rVert_2\le\eta'_{m,n}\lVert a_j\rVert_2 ∥ Δ a j ′ ∥ 2 ≤ η m , n ′ ∥ a j ∥ 2 (1 ≤ j ≤ n 1\le j\le n 1 ≤ j ≤ n )を、Δ b \Delta b Δ b が∥ Δ b ∥ 2 ≤ η m , n ∥ b ∥ 2 \lVert\Delta b\rVert_2\le\eta_{m,n}\lVert b\rVert_2 ∥ Δ b ∥ 2 ≤ η m , n ∥ b ∥ 2 を満たし、A + Δ A ′ A+\Delta A' A + Δ A ′ の列は一次独立であって、x ^ \hat x x ^ は( A + Δ A ′ ) x = b + Δ b (A+\Delta A')x=b+\Delta b ( A + Δ A ′ ) x = b + Δ b のただ一つの最小二乗解である。
証明. 定理 4.5 (1) のQ Q Q 、Δ A \Delta A Δ A 、Δ b \Delta b Δ b をとる。補題 2.1 (1) により∥ Δ A ∥ 2 ≤ ∥ Δ A ∥ F ≤ η m , n ∥ A ∥ F < σ n ( A ) \lVert\Delta A\rVert_2\le\lVert\Delta A\rVert_F\le\eta_{m,n}\lVert A\rVert_F<\sigma_n(A) ∥ Δ A ∥ 2 ≤ ∥ Δ A ∥ F ≤ η m , n ∥ A ∥ F < σ n ( A ) であるから、補題 2.1 (3) によりA + Δ A = Q ( R ^ 0 ) A+\Delta A=Q\begin{pmatrix}\widehat R\\0\end{pmatrix} A + Δ A = Q ( R 0 ) の列は一次独立であり、( R ^ 0 ) = Q T ( A + Δ A ) \begin{pmatrix}\widehat R\\0\end{pmatrix}=Q^{\mathsf T}(A+\Delta A) ( R 0 ) = Q T ( A + Δ A ) の列も一次独立である。したがってR ^ \widehat R R は正則であり、上三角行列R ^ \widehat R R の対角成分はすべて0 0 0 でない。これで(1) は示された。
n ≤ m n\le m n ≤ m と( 2 m + 6 ) u < 1 (2m+6)u<1 ( 2 m + 6 ) u < 1 からn u < n / ( 2 m + 6 ) < 1 / 2 nu<n/(2m+6)<1/2 n u < n / ( 2 m + 6 ) < 1/2 であり、γ n < 1 \gamma_n<1 γ n < 1 である。§E20.5 定理 3.8 (1) をT = R ^ ∈ M n ( F ) T=\widehat R\in M_n(F) T = R ∈ M n ( F ) 、c = c ^ 1 ∈ F n c=\widehat c_1\in F^n c = c 1 ∈ F n に適用すると、∣ E ∣ ≤ γ n ∣ R ^ ∣ |E|\le\gamma_n|\widehat R| ∣ E ∣ ≤ γ n ∣ R ∣ を満たすE ∈ R n × n E\in\R^{n\times n} E ∈ R n × n が存在して( R ^ + E ) x ^ = c ^ 1 (\widehat R+E)\hat x=\widehat c_1 ( R + E ) x ^ = c 1 である。∣ E ∣ ≤ γ n ∣ R ^ ∣ |E|\le\gamma_n|\widehat R| ∣ E ∣ ≤ γ n ∣ R ∣ からE E E は上三角行列であり、R ^ + E \widehat R+E R + E の対角成分は∣ r ^ i i + e i i ∣ ≥ ( 1 − γ n ) ∣ r ^ i i ∣ > 0 |\hat r_{ii}+e_{ii}|\ge(1-\gamma_n)|\hat r_{ii}|>0 ∣ r ^ ii + e ii ∣ ≥ ( 1 − γ n ) ∣ r ^ ii ∣ > 0 を満たす。Δ A ′ : = Δ A + Q ( E 0 ) \Delta A':=\Delta A+Q\begin{pmatrix}E\\0\end{pmatrix} Δ A ′ := Δ A + Q ( E 0 ) と置くと、A + Δ A ′ = Q ( R ^ + E 0 ) A+\Delta A'=Q\begin{pmatrix}\widehat R+E\\0\end{pmatrix} A + Δ A ′ = Q ( R + E 0 ) 、b + Δ b = Q c ^ b+\Delta b=Q\widehat c b + Δ b = Q c である。補題 3.5 をT = R ^ + E T=\widehat R+E T = R + E と右辺b + Δ b b+\Delta b b + Δ b に適用すると、Q T ( b + Δ b ) = c ^ Q^{\mathsf T}(b+\Delta b)=\widehat c Q T ( b + Δ b ) = c であるから、A + Δ A ′ A+\Delta A' A + Δ A ′ の列は一次独立であり、ただ一つの最小二乗解は( R ^ + E ) x = c ^ 1 (\widehat R+E)x=\widehat c_1 ( R + E ) x = c 1 の解x ^ \hat x x ^ である。Q ( E 0 ) Q\begin{pmatrix}E\\0\end{pmatrix} Q ( E 0 ) の第j j j 列のノルムはE E E の第j j j 列のノルムに等しく、∣ E ∣ ≤ γ n ∣ R ^ ∣ |E|\le\gamma_n|\widehat R| ∣ E ∣ ≤ γ n ∣ R ∣ からγ n ∥ R ^ e j ∥ 2 = γ n ∥ ( A + Δ A ) e j ∥ 2 ≤ γ n ( 1 + η m , n ) ∥ a j ∥ 2 \gamma_n\lVert\widehat Re_j\rVert_2=\gamma_n\lVert(A+\Delta A)e_j\rVert_2\le\gamma_n(1+\eta_{m,n})\lVert a_j\rVert_2 γ n ∥ R e j ∥ 2 = γ n ∥( A + Δ A ) e j ∥ 2 ≤ γ n ( 1 + η m , n ) ∥ a j ∥ 2 以下である。これと定理 4.5 (1) の∥ Δ a j ∥ 2 ≤ η m , n ∥ a j ∥ 2 \lVert\Delta a_j\rVert_2\le\eta_{m,n}\lVert a_j\rVert_2 ∥ Δ a j ∥ 2 ≤ η m , n ∥ a j ∥ 2 から∥ Δ a j ′ ∥ 2 ≤ η m , n ′ ∥ a j ∥ 2 \lVert\Delta a'_j\rVert_2\le\eta'_{m,n}\lVert a_j\rVert_2 ∥ Δ a j ′ ∥ 2 ≤ η m , n ′ ∥ a j ∥ 2 である。▨
例 4.8. F = F ( 2 , 53 , − 1022 , 1023 ) F=F(2,53,-1022,1023) F = F ( 2 , 53 , − 1022 , 1023 ) 、fl \operatorname{fl} fl を最近接偶数丸め、u = 2 − 53 u=2^{-53} u = 2 − 53 、ε : = 2 − 27 \varepsilon:=2^{-27} ε := 2 − 27 とし、
A : = ( 1 1 1 ε 0 0 0 ε 0 0 0 ε ) , b : = ( 3 ε ε ε ) = A ( 1 1 1 ) A:=\begin{pmatrix}1&1&1\\\varepsilon&0&0\\0&\varepsilon&0\\0&0&\varepsilon\end{pmatrix},\qquad b:=\begin{pmatrix}3\\\varepsilon\\\varepsilon\\\varepsilon\end{pmatrix}=A\begin{pmatrix}1\\1\\1\end{pmatrix} A := 1 ε 0 0 1 0 ε 0 1 0 0 ε , b := 3 ε ε ε = A 1 1 1 とする。最小二乗解はx = ( 1 , 1 , 1 ) x=(1,1,1) x = ( 1 , 1 , 1 ) 、残差は0 0 0 である。J J J を成分がすべて1 1 1 の3 × 3 3\times3 3 × 3 行列とするとA T A = J + ε 2 I 3 A^{\mathsf T}A=J+\varepsilon^2I_3 A T A = J + ε 2 I 3 であり、固有値は3 + ε 2 , ε 2 , ε 2 3+\varepsilon^2,\varepsilon^2,\varepsilon^2 3 + ε 2 , ε 2 , ε 2 であるから、命題 2.3 によりA A A の特異値は3 + ε 2 , ε , ε \sqrt{3+\varepsilon^2},\varepsilon,\varepsilon 3 + ε 2 , ε , ε 、κ 2 ( A ) = 3 + ε 2 / ε ≈ 2.32 × 10 8 \kappa_2(A)=\sqrt{3+\varepsilon^2}/\varepsilon\approx2.32\times10^8 κ 2 ( A ) = 3 + ε 2 / ε ≈ 2.32 × 1 0 8 である。
正規方程式:A T A A^{\mathsf T}A A T A の対角成分の計算値はfl ( 1 + 2 − 54 ) \operatorname{fl}(1+2^{-54}) fl ( 1 + 2 − 54 ) であり、1 + 2 − 54 1+2^{-54} 1 + 2 − 54 の最近接点は、1 1 1 との距離2 − 54 2^{-54} 2 − 54 が1 + 2 − 52 1+2^{-52} 1 + 2 − 52 との距離3 ⋅ 2 − 54 3\cdot2^{-54} 3 ⋅ 2 − 54 より小さいので1 1 1 である。非対角成分は1 1 1 であるから、浮動小数点算術で作ったA T A A^{\mathsf T}A A T A はJ J J であり、A A A の列は一次独立であるがJ J J は階数1 1 1 の特異行列である。A T b A^{\mathsf T}b A T b の各成分の厳密な値は3 + 2 − 54 3+2^{-54} 3 + 2 − 54 であり、その最近接点は3 3 3 であるから、A T b A^{\mathsf T}b A T b の計算値は( 3 , 3 , 3 ) (3,3,3) ( 3 , 3 , 3 ) である。J J J に§E20.5 定理 5.2 (5) の式を適用するとc 11 = 1 c_{11}=1 c 11 = 1 、c 21 = 1 c_{21}=1 c 21 = 1 であり、c 22 c_{22} c 22 の根号の中は1 − 1 = 0 1-1=0 1 − 1 = 0 であって正でない。
Gram–Schmidt 法: 古典 Gram–Schmidt 法(CGS)はj = 1 , 2 , 3 j=1,2,3 j = 1 , 2 , 3 の順に、r i j : = q i T a j r_{ij}:=q_i^{\mathsf T}a_j r ij := q i T a j (i < j i<j i < j )を元の列a j a_j a j との内積として計算し、w : = a j − ∑ i < j r i j q i w:=a_j-\sum_{i<j}r_{ij}q_i w := a j − ∑ i < j r ij q i 、r j j : = ∥ w ∥ 2 r_{jj}:=\lVert w\rVert_2 r j j := ∥ w ∥ 2 、q j : = w / r j j q_j:=w/r_{jj} q j := w / r j j と置く。修正 Gram–Schmidt 法(MGS)はw : = a j w:=a_j w := a j と置き、i = 1 , … , j − 1 i=1,\dots,j-1 i = 1 , … , j − 1 の順にr i j : = q i T w r_{ij}:=q_i^{\mathsf T}w r ij := q i T w 、w : = w − r i j q i w:=w-r_{ij}q_i w := w − r ij q i と更新してからr j j r_{jj} r j j 、q j q_j q j を同じく定める。q 1 , … , q i − 1 q_1,\dots,q_{i-1} q 1 , … , q i − 1 が正規直交ならば MGS のq i T w q_i^{\mathsf T}w q i T w はq i T a j q_i^{\mathsf T}a_j q i T a j に等しいので、厳密算術では二つの計算式は§D3.14 定理 2.1 と同じq j q_j q j を与える。内積とノルムを逐次和と正しく丸めた平方根で計算し、CPython の float と math.sqrt で実行して次の値を観察した。
CGS ではr ^ 11 = 1 \hat r_{11}=1 r ^ 11 = 1 、q ^ 1 = ( 1 , ε , 0 , 0 ) \hat q_1=(1,\varepsilon,0,0) q ^ 1 = ( 1 , ε , 0 , 0 ) 、r ^ 12 = r ^ 13 = 1 \hat r_{12}=\hat r_{13}=1 r ^ 12 = r ^ 13 = 1 、r ^ 23 = 0 \hat r_{23}=0 r ^ 23 = 0 であり、q ^ 2 = ( 0 , − g , g , 0 ) \hat q_2=(0,-g,g,0) q ^ 2 = ( 0 , − g , g , 0 ) 、q ^ 3 = ( 0 , − g , 0 , g ) \hat q_3=(0,-g,0,g) q ^ 3 = ( 0 , − g , 0 , g ) (g = 0.7071067811865475 g=0.7071067811865475 g = 0.7071067811865475 )である。q ^ 2 T q ^ 3 = g 2 ≈ 0.5 \hat q_2^{\mathsf T}\hat q_3=g^2\approx0.5 q ^ 2 T q ^ 3 = g 2 ≈ 0.5 、∥ Q ^ T Q ^ − I ∥ F ≈ 0.707 \lVert\widehat Q^{\mathsf T}\widehat Q-I\rVert_F\approx0.707 ∥ Q T Q − I ∥ F ≈ 0.707 である。Q ^ T b \widehat Q^{\mathsf T}b Q T b の計算値は( 3 , 0 , 0 ) (3,0,0) ( 3 , 0 , 0 ) であり、後退代入の計算値はx ^ = ( 3 , 0 , 0 ) \hat x=(3,0,0) x ^ = ( 3 , 0 , 0 ) 、相対誤差∥ x ^ − x ∥ 2 / ∥ x ∥ 2 = 2 \lVert\hat x-x\rVert_2/\lVert x\rVert_2=\sqrt2 ∥ x ^ − x ∥ 2 / ∥ x ∥ 2 = 2 である。
MGS では∥ Q ^ T Q ^ − I ∥ F ≈ 8.60 × 10 − 9 \lVert\widehat Q^{\mathsf T}\widehat Q-I\rVert_F\approx8.60\times10^{-9} ∥ Q T Q − I ∥ F ≈ 8.60 × 1 0 − 9 である。Q ^ T b \widehat Q^{\mathsf T}b Q T b を内積で計算して解くとx ^ ≈ ( 3 , − 4.53 × 10 − 17 , 9.06 × 10 − 17 ) \hat x\approx(3,-4.53\times10^{-17},9.06\times10^{-17}) x ^ ≈ ( 3 , − 4.53 × 1 0 − 17 , 9.06 × 1 0 − 17 ) 、相対誤差は約1.41 1.41 1.41 である。b b b を第4 4 4 列として同じ更新c i : = q ^ i T w c_i:=\hat q_i^{\mathsf T}w c i := q ^ i T w 、w : = w − c i q ^ i w:=w-c_i\hat q_i w := w − c i q ^ i で処理するとc ≈ ( 3 , 1.58 × 10 − 8 , 9.13 × 10 − 9 ) c\approx(3,1.58\times10^{-8},9.13\times10^{-9}) c ≈ ( 3 , 1.58 × 1 0 − 8 , 9.13 × 1 0 − 9 ) であり、相対誤差は約1.92 × 10 − 16 1.92\times10^{-16} 1.92 × 1 0 − 16 である。
Householder QR 法(右辺b b b )ではR ^ \widehat R R の第1 1 1 行は( − 1 , − 1 , − 1 ) (-1,-1,-1) ( − 1 , − 1 , − 1 ) 、c ^ ≈ ( − 3 , 1.58 × 10 − 8 , 9.13 × 10 − 9 , 0 ) \widehat c\approx(-3,1.58\times10^{-8},9.13\times10^{-9},0) c ≈ ( − 3 , 1.58 × 1 0 − 8 , 9.13 × 1 0 − 9 , 0 ) であり、x ^ \hat x x ^ の相対誤差は約3.20 × 10 − 16 3.20\times10^{-16} 3.20 × 1 0 − 16 、保存した反射から作ったQ ^ \widehat Q Q について∥ Q ^ T Q ^ − I ∥ F ≈ 4.95 × 10 − 16 \lVert\widehat Q^{\mathsf T}\widehat Q-I\rVert_F\approx4.95\times10^{-16} ∥ Q T Q − I ∥ F ≈ 4.95 × 1 0 − 16 である。
x ^ = ( 3 , 0 , 0 ) \hat x=(3,0,0) x ^ = ( 3 , 0 , 0 ) ではx ^ − x = ( 2 , − 1 , − 1 ) \hat x-x=(2,-1,-1) x ^ − x = ( 2 , − 1 , − 1 ) は( 1 , 1 , 1 ) (1,1,1) ( 1 , 1 , 1 ) に直交し、A T A A^{\mathsf T}A A T A の固有値ε 2 \varepsilon^2 ε 2 の固有ベクトルであるから、残差はb − A x ^ = A ( x − x ^ ) = ( 0 , − 2 ε , ε , ε ) b-A\hat x=A(x-\hat x)=(0,-2\varepsilon,\varepsilon,\varepsilon) b − A x ^ = A ( x − x ^ ) = ( 0 , − 2 ε , ε , ε ) 、∥ b − A x ^ ∥ 2 = ε ∥ x ^ − x ∥ 2 = 6 ε ≈ 1.83 × 10 − 8 \lVert b-A\hat x\rVert_2=\varepsilon\lVert\hat x-x\rVert_2=\sqrt6\varepsilon\approx1.83\times10^{-8} ∥ b − A x ^ ∥ 2 = ε ∥ x ^ − x ∥ 2 = 6 ε ≈ 1.83 × 1 0 − 8 、∥ b − A x ^ ∥ 2 / ∥ b ∥ 2 ≈ 6.1 × 10 − 9 \lVert b-A\hat x\rVert_2/\lVert b\rVert_2\approx6.1\times10^{-9} ∥ b − A x ^ ∥ 2 / ∥ b ∥ 2 ≈ 6.1 × 1 0 − 9 である。
5 重み付き最小二乗問題
定義 5.1. A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n 、b ∈ R m b\in\R^m b ∈ R m とし、W ∈ R m × m W\in\R^{m\times m} W ∈ R m × m を実対称正定値行列とする。x ∈ R n x\in\R^n x ∈ R n が任意のy ∈ R n y\in\R^n y ∈ R n に対して
( b − A x ) T W ( b − A x ) ≤ ( b − A y ) T W ( b − A y ) (b-Ax)^{\mathsf T}W(b-Ax)\le(b-Ay)^{\mathsf T}W(b-Ay) ( b − A x ) T W ( b − A x ) ≤ ( b − A y ) T W ( b − A y ) を満たすとき、x x x をA x = b Ax=b A x = b の重みW W W に関する 重み付き最小二乗解 (weighted least squares solution ) という。
命題 5.2. A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n 、b ∈ R m b\in\R^m b ∈ R m とし、W ∈ R m × m W\in\R^{m\times m} W ∈ R m × m を実対称正定値行列とする。
R m \R^m R m に標準内積を入れるとW W W は正作用素であり、その正の平方根W 1 / 2 W^{1/2} W 1/2 (§E3.38 定理 1.1 )は実対称行列であって( W 1 / 2 ) T W 1 / 2 = W (W^{1/2})^{\mathsf T}W^{1/2}=W ( W 1/2 ) T W 1/2 = W を満たす。W W W の Cholesky 因子L L L について( L T ) T L T = W (L^{\mathsf T})^{\mathsf T}L^{\mathsf T}=W ( L T ) T L T = W である。W = diag ( w 1 , … , w m ) W=\operatorname{diag}(w_1,\dots,w_m) W = diag ( w 1 , … , w m ) ならば、w i > 0 w_i>0 w i > 0 でありW 1 / 2 = diag ( w 1 , … , w m ) W^{1/2}=\operatorname{diag}(\sqrt{w_1},\dots,\sqrt{w_m}) W 1/2 = diag ( w 1 , … , w m ) である。
C ∈ R m × m C\in\R^{m\times m} C ∈ R m × m がC T C = W C^{\mathsf T}C=W C T C = W を満たすならば、C C C は正則であり、任意のx ∈ R n x\in\R^n x ∈ R n に対して( b − A x ) T W ( b − A x ) = ∥ C b − C A x ∥ 2 2 (b-Ax)^{\mathsf T}W(b-Ax)=\lVert Cb-CAx\rVert_2^2 ( b − A x ) T W ( b − A x ) = ∥ C b − C A x ∥ 2 2 が成り立つ。
A A A の列が一次独立であるとし、C ∈ R m × m C\in\R^{m\times m} C ∈ R m × m がC T C = W C^{\mathsf T}C=W C T C = W を満たすとする。このときC A CA C A の列は一次独立であり、A x = b Ax=b A x = b の重みW W W に関する重み付き最小二乗解はただ一つ存在して、( C A ) x = C b (CA)x=Cb ( C A ) x = C b のただ一つの最小二乗解( A T W A ) − 1 A T W b (A^{\mathsf T}WA)^{-1}A^{\mathsf T}Wb ( A T W A ) − 1 A T W b に等しい。
証明. W W W は実対称であり、任意のx x x に対してx T W x ≥ 0 x^{\mathsf T}Wx\ge0 x T W x ≥ 0 であるから、標準内積に関する正作用素である。§E3.38 定理 1.1 により正作用素W 1 / 2 W^{1/2} W 1/2 が存在して( W 1 / 2 ) 2 = W (W^{1/2})^2=W ( W 1/2 ) 2 = W であり、正作用素は自己随伴であるからW 1 / 2 W^{1/2} W 1/2 は実対称行列であって( W 1 / 2 ) T W 1 / 2 = W (W^{1/2})^{\mathsf T}W^{1/2}=W ( W 1/2 ) T W 1/2 = W である。§E20.5 定理 5.2 (3) によりW = L L T W=LL^{\mathsf T} W = L L T である。W W W が対角行列ならばw i = e i T W e i > 0 w_i=e_i^{\mathsf T}We_i>0 w i = e i T W e i > 0 であり、diag ( w i ) \operatorname{diag}(\sqrt{w_i}) diag ( w i ) は正作用素であってその二乗はW W W であるから、§E3.38 定理 1.1 の一意性によりW 1 / 2 W^{1/2} W 1/2 に等しい。これで(1) は示された。
C T C = W C^{\mathsf T}C=W C T C = W とする。C x = 0 Cx=0 C x = 0 ならばx T W x = ∥ C x ∥ 2 2 = 0 x^{\mathsf T}Wx=\lVert Cx\rVert_2^2=0 x T W x = ∥ C x ∥ 2 2 = 0 であり、W W W は正定値であるからx = 0 x=0 x = 0 である。したがってC C C は正則である。( b − A x ) T C T C ( b − A x ) = ∥ C ( b − A x ) ∥ 2 2 (b-Ax)^{\mathsf T}C^{\mathsf T}C(b-Ax)=\lVert C(b-Ax)\rVert_2^2 ( b − A x ) T C T C ( b − A x ) = ∥ C ( b − A x ) ∥ 2 2 である。これで(2) は示された。
(3) の仮定の下で、C A h = 0 CAh=0 C A h = 0 ならばC C C の正則性からA h = 0 Ah=0 A h = 0 であり、h = 0 h=0 h = 0 である。(2) により、x x x が重み付き最小二乗解であることと、x x x が( C A ) x = C b (CA)x=Cb ( C A ) x = C b の最小二乗解であることは同値である。命題 1.2 (2) をC A CA C A とC b Cb C b に適用すると、最小二乗解はただ一つであり、( C A ) T C A = A T W A (CA)^{\mathsf T}CA=A^{\mathsf T}WA ( C A ) T C A = A T W A 、( C A ) T C b = A T W b (CA)^{\mathsf T}Cb=A^{\mathsf T}Wb ( C A ) T C b = A T W b から、それは( A T W A ) − 1 A T W b (A^{\mathsf T}WA)^{-1}A^{\mathsf T}Wb ( A T W A ) − 1 A T W b である。▨
例 5.3. 直線y = x 1 + x 2 t y=x_1+x_2t y = x 1 + x 2 t を四つの観測( t i , y i ) = ( 0 , 1 ) , ( 1 , 2 ) , ( 2 , 4 ) , ( 3 , 3 ) (t_i,y_i)=(0,1),(1,2),(2,4),(3,3) ( t i , y i ) = ( 0 , 1 ) , ( 1 , 2 ) , ( 2 , 4 ) , ( 3 , 3 ) に当てはめる。第1 1 1 群(t = 0 , 1 t=0,1 t = 0 , 1 )の測定の標準偏差を0.1 0.1 0.1 、第2 2 2 群(t = 2 , 3 t=2,3 t = 2 , 3 )の標準偏差を1 1 1 とし、重みを標準偏差の二乗の逆数W : = diag ( 100 , 100 , 1 , 1 ) W:=\operatorname{diag}(100,100,1,1) W := diag ( 100 , 100 , 1 , 1 ) とする。A A A の第i i i 行を( 1 , t i ) (1,t_i) ( 1 , t i ) 、b : = ( 1 , 2 , 4 , 3 ) b:=(1,2,4,3) b := ( 1 , 2 , 4 , 3 ) とする。命題 5.2 (1) によりC : = W 1 / 2 = diag ( 10 , 10 , 1 , 1 ) C:=W^{1/2}=\operatorname{diag}(10,10,1,1) C := W 1/2 = diag ( 10 , 10 , 1 , 1 ) であり、
C A = ( 10 0 10 10 1 2 1 3 ) , C b = ( 10 20 4 3 ) CA=\begin{pmatrix}10&0\\10&10\\1&2\\1&3\end{pmatrix},\qquad Cb=\begin{pmatrix}10\\20\\4\\3\end{pmatrix} C A = 10 10 1 1 0 10 2 3 , C b = 10 20 4 3 である。命題 5.2 (3) により、重み付き最小二乗解は( C A ) x = C b (CA)x=Cb ( C A ) x = C b の最小二乗解であり、Householder QR 法を( C A , C b ) (CA,Cb) ( C A , C b ) に適用して計算することができる。その値はx W = ( 11906 / 11801 , 11599 / 11801 ) ≈ ( 1.00890 , 0.98288 ) x_W=(11906/11801,\ 11599/11801)\approx(1.00890,0.98288) x W = ( 11906/11801 , 11599/11801 ) ≈ ( 1.00890 , 0.98288 ) 、残差は
b − A x W = 1 11801 ( − 105 , 97 , 12100 , − 11300 ) ≈ ( − 0.0089 , 0.0082 , 1.0253 , − 0.9575 ) b-Ax_W=\frac{1}{11801}(-105,\ 97,\ 12100,\ -11300)\approx(-0.0089,\ 0.0082,\ 1.0253,\ -0.9575) b − A x W = 11801 1 ( − 105 , 97 , 12100 , − 11300 ) ≈ ( − 0.0089 , 0.0082 , 1.0253 , − 0.9575 ) である。重みなしの最小二乗解は( 13 / 10 , 4 / 5 ) (13/10,4/5) ( 13/10 , 4/5 ) 、残差は( − 0.3 , − 0.1 , 1.1 , − 0.7 ) (-0.3,-0.1,1.1,-0.7) ( − 0.3 , − 0.1 , 1.1 , − 0.7 ) である。C ( b − A x W ) ≈ ( − 0.089 , 0.082 , 1.025 , − 0.958 ) C(b-Ax_W)\approx(-0.089,0.082,1.025,-0.958) C ( b − A x W ) ≈ ( − 0.089 , 0.082 , 1.025 , − 0.958 ) の各成分は残差を各観測の標準偏差で割った値であり、その絶対値は1.03 1.03 1.03 以下である。
6 最小二乗解の感度
定理 6.1. m ≥ n m\ge n m ≥ n とし、A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n の列が一次独立であるとする。σ 1 : = σ 1 ( A ) \sigma_1:=\sigma_1(A) σ 1 := σ 1 ( A ) 、σ n : = σ n ( A ) \sigma_n:=\sigma_n(A) σ n := σ n ( A ) 、G : = A T A G:=A^{\mathsf T}A G := A T A と置き、b ∈ R m b\in\R^m b ∈ R m 、x x x をA x = b Ax=b A x = b の最小二乗解、r : = b − A x r:=b-Ax r := b − A x とする。Δ A ∈ R m × n \Delta A\in\R^{m\times n} Δ A ∈ R m × n とΔ b ∈ R m \Delta b\in\R^m Δ b ∈ R m が∥ Δ A ∥ 2 ≤ σ n / 2 \lVert\Delta A\rVert_2\le\sigma_n/2 ∥ Δ A ∥ 2 ≤ σ n /2 を満たすとし、f : = Δ b − Δ A x f:=\Delta b-\Delta Ax f := Δ b − Δ A x と置く。
A + Δ A A+\Delta A A + Δ A の列は一次独立であり、( A + Δ A ) T ( A + Δ A ) (A+\Delta A)^{\mathsf T}(A+\Delta A) ( A + Δ A ) T ( A + Δ A ) の逆行列は∥ ( ( A + Δ A ) T ( A + Δ A ) ) − 1 ∥ 2 ≤ 4 / σ n 2 \lVert((A+\Delta A)^{\mathsf T}(A+\Delta A))^{-1}\rVert_2\le4/\sigma_n^2 ∥(( A + Δ A ) T ( A + Δ A ) ) − 1 ∥ 2 ≤ 4/ σ n 2 を満たす。
( A + Δ A ) x ′ = b + Δ b (A+\Delta A)x'=b+\Delta b ( A + Δ A ) x ′ = b + Δ b のただ一つの最小二乗解をx ′ x' x ′ とし、Δ x : = x ′ − x \Delta x:=x'-x Δ x := x ′ − x と置く。このとき
Δ x = G − 1 A T f + G − 1 ( Δ A ) T r + ρ \Delta x=G^{-1}A^{\mathsf T}f+G^{-1}(\Delta A)^{\mathsf T}r+\rho Δ x = G − 1 A T f + G − 1 ( Δ A ) T r + ρ
であり、ρ ∈ R n \rho\in\R^n ρ ∈ R n は
∥ ρ ∥ 2 ≤ 4 ∥ Δ A ∥ 2 σ n 2 ( ( 2 σ 1 + ∥ Δ A ∥ 2 ) ( ∥ f ∥ 2 σ n + ∥ Δ A ∥ 2 ∥ r ∥ 2 σ n 2 ) + ∥ f ∥ 2 ) \lVert\rho\rVert_2\le\frac{4\lVert\Delta A\rVert_2}{\sigma_n^2}\Bigl((2\sigma_1+\lVert\Delta A\rVert_2)\Bigl(\frac{\lVert f\rVert_2}{\sigma_n}+\frac{\lVert\Delta A\rVert_2\lVert r\rVert_2}{\sigma_n^2}\Bigr)+\lVert f\rVert_2\Bigr) ∥ ρ ∥ 2 ≤ σ n 2 4 ∥ Δ A ∥ 2 ( ( 2 σ 1 + ∥ Δ A ∥ 2 ) ( σ n ∥ f ∥ 2 + σ n 2 ∥ Δ A ∥ 2 ∥ r ∥ 2 ) + ∥ f ∥ 2 )
を満たす。
一次の項L : = G − 1 A T f + G − 1 ( Δ A ) T r L:=G^{-1}A^{\mathsf T}f+G^{-1}(\Delta A)^{\mathsf T}r L := G − 1 A T f + G − 1 ( Δ A ) T r は∥ L ∥ 2 ≤ ∥ f ∥ 2 / σ n + ∥ Δ A ∥ 2 ∥ r ∥ 2 / σ n 2 \lVert L\rVert_2\le\lVert f\rVert_2/\sigma_n+\lVert\Delta A\rVert_2\lVert r\rVert_2/\sigma_n^2 ∥ L ∥ 2 ≤ ∥ f ∥ 2 / σ n + ∥ Δ A ∥ 2 ∥ r ∥ 2 / σ n 2 を満たす。x ≠ 0 x\ne0 x = 0 ならば、κ : = κ 2 ( A ) \kappa:=\kappa_2(A) κ := κ 2 ( A ) について
∥ L ∥ 2 ∥ x ∥ 2 ≤ κ ( ∥ Δ b ∥ 2 σ 1 ∥ x ∥ 2 + ∥ Δ A ∥ 2 σ 1 ) + κ 2 ∥ Δ A ∥ 2 σ 1 ⋅ ∥ r ∥ 2 σ 1 ∥ x ∥ 2 \frac{\lVert L\rVert_2}{\lVert x\rVert_2}\le\kappa\Bigl(\frac{\lVert\Delta b\rVert_2}{\sigma_1\lVert x\rVert_2}+\frac{\lVert\Delta A\rVert_2}{\sigma_1}\Bigr)+\kappa^2\frac{\lVert\Delta A\rVert_2}{\sigma_1}\cdot\frac{\lVert r\rVert_2}{\sigma_1\lVert x\rVert_2} ∥ x ∥ 2 ∥ L ∥ 2 ≤ κ ( σ 1 ∥ x ∥ 2 ∥ Δ b ∥ 2 + σ 1 ∥ Δ A ∥ 2 ) + κ 2 σ 1 ∥ Δ A ∥ 2 ⋅ σ 1 ∥ x ∥ 2 ∥ r ∥ 2
が成り立つ。
証明. A ′ : = A + Δ A A':=A+\Delta A A ′ := A + Δ A 、G ′ : = A ′ T A ′ G':=A'^{\mathsf T}A' G ′ := A ′ T A ′ と置く。∥ Δ A ∥ 2 ≤ σ n / 2 < σ n \lVert\Delta A\rVert_2\le\sigma_n/2<\sigma_n ∥ Δ A ∥ 2 ≤ σ n /2 < σ n であるから、補題 2.1 (3) によりA ′ A' A ′ の列は一次独立であってσ n ( A ′ ) ≥ σ n / 2 \sigma_n(A')\ge\sigma_n/2 σ n ( A ′ ) ≥ σ n /2 であり、補題 2.1 (2) により∥ G ′ − 1 ∥ 2 = σ n ( A ′ ) − 2 ≤ 4 / σ n 2 \lVert G'^{-1}\rVert_2=\sigma_n(A')^{-2}\le4/\sigma_n^2 ∥ G ′ − 1 ∥ 2 = σ n ( A ′ ) − 2 ≤ 4/ σ n 2 である。これで(1) は示された。
命題 1.2 条件 (c) によりA T r = 0 A^{\mathsf T}r=0 A T r = 0 であり、G ′ x ′ = A ′ T ( b + Δ b ) G'x'=A'^{\mathsf T}(b+\Delta b) G ′ x ′ = A ′ T ( b + Δ b ) である。b + Δ b − A ′ x = r + f b+\Delta b-A'x=r+f b + Δ b − A ′ x = r + f であるから、
G ′ Δ x = A ′ T ( b + Δ b ) − G ′ x = A ′ T ( r + f ) = A T f + ( Δ A ) T r + ( Δ A ) T f G'\Delta x=A'^{\mathsf T}(b+\Delta b)-G'x=A'^{\mathsf T}(r+f)=A^{\mathsf T}f+(\Delta A)^{\mathsf T}r+(\Delta A)^{\mathsf T}f G ′ Δ x = A ′ T ( b + Δ b ) − G ′ x = A ′ T ( r + f ) = A T f + ( Δ A ) T r + ( Δ A ) T f である。g : = A T f + ( Δ A ) T r g:=A^{\mathsf T}f+(\Delta A)^{\mathsf T}r g := A T f + ( Δ A ) T r と置くとL = G − 1 g L=G^{-1}g L = G − 1 g であり、G ′ − 1 − G − 1 = G ′ − 1 ( G − G ′ ) G − 1 G'^{-1}-G^{-1}=G'^{-1}(G-G')G^{-1} G ′ − 1 − G − 1 = G ′ − 1 ( G − G ′ ) G − 1 から
ρ : = Δ x − L = G ′ − 1 ( G − G ′ ) G − 1 g + G ′ − 1 ( Δ A ) T f \rho:=\Delta x-L=G'^{-1}(G-G')G^{-1}g+G'^{-1}(\Delta A)^{\mathsf T}f ρ := Δ x − L = G ′ − 1 ( G − G ′ ) G − 1 g + G ′ − 1 ( Δ A ) T f である。G ′ − G = A T Δ A + ( Δ A ) T A + ( Δ A ) T Δ A G'-G=A^{\mathsf T}\Delta A+(\Delta A)^{\mathsf T}A+(\Delta A)^{\mathsf T}\Delta A G ′ − G = A T Δ A + ( Δ A ) T A + ( Δ A ) T Δ A であり、補題 2.1 (1) により∥ A T ∥ 2 = σ 1 \lVert A^{\mathsf T}\rVert_2=\sigma_1 ∥ A T ∥ 2 = σ 1 、∥ ( Δ A ) T ∥ 2 = ∥ Δ A ∥ 2 \lVert(\Delta A)^{\mathsf T}\rVert_2=\lVert\Delta A\rVert_2 ∥( Δ A ) T ∥ 2 = ∥ Δ A ∥ 2 であるから、∥ G ′ − G ∥ 2 ≤ ( 2 σ 1 + ∥ Δ A ∥ 2 ) ∥ Δ A ∥ 2 \lVert G'-G\rVert_2\le(2\sigma_1+\lVert\Delta A\rVert_2)\lVert\Delta A\rVert_2 ∥ G ′ − G ∥ 2 ≤ ( 2 σ 1 + ∥ Δ A ∥ 2 ) ∥ Δ A ∥ 2 である。補題 2.1 (2) により∥ G − 1 A T ∥ 2 = 1 / σ n \lVert G^{-1}A^{\mathsf T}\rVert_2=1/\sigma_n ∥ G − 1 A T ∥ 2 = 1/ σ n 、∥ G − 1 ∥ 2 = 1 / σ n 2 \lVert G^{-1}\rVert_2=1/\sigma_n^2 ∥ G − 1 ∥ 2 = 1/ σ n 2 であるから、
∥ L ∥ 2 ≤ ∥ f ∥ 2 σ n + ∥ Δ A ∥ 2 ∥ r ∥ 2 σ n 2 \lVert L\rVert_2\le\frac{\lVert f\rVert_2}{\sigma_n}+\frac{\lVert\Delta A\rVert_2\lVert r\rVert_2}{\sigma_n^2} ∥ L ∥ 2 ≤ σ n ∥ f ∥ 2 + σ n 2 ∥ Δ A ∥ 2 ∥ r ∥ 2 である。これらと∥ G ′ − 1 ∥ 2 ≤ 4 / σ n 2 \lVert G'^{-1}\rVert_2\le4/\sigma_n^2 ∥ G ′ − 1 ∥ 2 ≤ 4/ σ n 2 からρ \rho ρ の評価を得る。これで(2) は示された。
∥ L ∥ 2 \lVert L\rVert_2 ∥ L ∥ 2 の評価は上で得た。∥ f ∥ 2 ≤ ∥ Δ b ∥ 2 + ∥ Δ A ∥ 2 ∥ x ∥ 2 \lVert f\rVert_2\le\lVert\Delta b\rVert_2+\lVert\Delta A\rVert_2\lVert x\rVert_2 ∥ f ∥ 2 ≤ ∥ Δ b ∥ 2 + ∥ Δ A ∥ 2 ∥ x ∥ 2 を代入して∥ x ∥ 2 \lVert x\rVert_2 ∥ x ∥ 2 で割り、1 / σ n = κ / σ 1 1/\sigma_n=\kappa/\sigma_1 1/ σ n = κ / σ 1 、1 / σ n 2 = κ 2 / σ 1 2 1/\sigma_n^2=\kappa^2/\sigma_1^2 1/ σ n 2 = κ 2 / σ 1 2 を用いると(3) の第二の不等式を得る。▨
例 6.3. ε : = 10 − 3 \varepsilon:=10^{-3} ε := 1 0 − 3 、a 1 : = ( 1 , 1 , 1 ) a_1:=(1,1,1) a 1 := ( 1 , 1 , 1 ) 、a 2 : = ( 1 − ε , 1 , 1 + ε ) a_2:=(1-\varepsilon,1,1+\varepsilon) a 2 := ( 1 − ε , 1 , 1 + ε ) とし、A : = ( a 1 a 2 ) ∈ R 3 × 2 A:=(a_1\ a_2)\in\R^{3\times2} A := ( a 1 a 2 ) ∈ R 3 × 2 とする。A T A = ( 3 3 3 3 + 2 ε 2 ) A^{\mathsf T}A=\begin{pmatrix}3&3\\3&3+2\varepsilon^2\end{pmatrix} A T A = ( 3 3 3 3 + 2 ε 2 ) であり、σ 1 ( A ) ≈ 2.449 \sigma_1(A)\approx2.449 σ 1 ( A ) ≈ 2.449 、σ 2 ( A ) ≈ 9.9999992 × 10 − 4 \sigma_2(A)\approx9.9999992\times10^{-4} σ 2 ( A ) ≈ 9.9999992 × 1 0 − 4 、κ 2 ( A ) ≈ 2.449 × 10 3 \kappa_2(A)\approx2.449\times10^3 κ 2 ( A ) ≈ 2.449 × 1 0 3 である。
残差が0 0 0 のモデル:b 1 : = A ( 1 , 1 ) = ( 2 − ε , 2 , 2 + ε ) b_1:=A(1,1)=(2-\varepsilon,2,2+\varepsilon) b 1 := A ( 1 , 1 ) = ( 2 − ε , 2 , 2 + ε ) とする。A x = b 1 Ax=b_1 A x = b 1 の最小二乗解はx = ( 1 , 1 ) x=(1,1) x = ( 1 , 1 ) 、残差は0 0 0 である。Δ b : = 10 − 8 ( − 1 , 0 , 1 ) \Delta b:=10^{-8}(-1,0,1) Δ b := 1 0 − 8 ( − 1 , 0 , 1 ) に対して、Δ A = 0 \Delta A=0 Δ A = 0 の場合の定理 6.1 (2) のρ \rho ρ は0 0 0 であり、A T ( − 1 , 0 , 1 ) = ( 0 , 2 ε ) A^{\mathsf T}(-1,0,1)=(0,2\varepsilon) A T ( − 1 , 0 , 1 ) = ( 0 , 2 ε ) からΔ x = ( A T A ) − 1 A T Δ b = ( − 10 − 5 , 10 − 5 ) \Delta x=(A^{\mathsf T}A)^{-1}A^{\mathsf T}\Delta b=(-10^{-5},10^{-5}) Δ x = ( A T A ) − 1 A T Δ b = ( − 1 0 − 5 , 1 0 − 5 ) である。∥ Δ b ∥ 2 / ∥ b 1 ∥ 2 ≈ 4.1 × 10 − 9 \lVert\Delta b\rVert_2/\lVert b_1\rVert_2\approx4.1\times10^{-9} ∥ Δ b ∥ 2 / ∥ b 1 ∥ 2 ≈ 4.1 × 1 0 − 9 の摂動に対して∥ Δ x ∥ 2 / ∥ x ∥ 2 = 10 − 5 \lVert\Delta x\rVert_2/\lVert x\rVert_2=10^{-5} ∥ Δ x ∥ 2 / ∥ x ∥ 2 = 1 0 − 5 である。残差の大きいモデル: 同じb 1 b_1 b 1 に定数y = z y=z y = z を当てはめる。係数行列はa 1 a_1 a 1 であり、κ 2 ( a 1 ) = 1 \kappa_2(a_1)=1 κ 2 ( a 1 ) = 1 、最小二乗解はz = 2 z=2 z = 2 、残差は( − ε , 0 , ε ) (-\varepsilon,0,\varepsilon) ( − ε , 0 , ε ) でノルムは2 ε ≈ 1.41 × 10 − 3 \sqrt2\,\varepsilon\approx1.41\times10^{-3} 2 ε ≈ 1.41 × 1 0 − 3 である。同じΔ b \Delta b Δ b に対してa 1 T Δ b = 0 a_1^{\mathsf T}\Delta b=0 a 1 T Δ b = 0 であるからΔ z = 0 \Delta z=0 Δ z = 0 であり、任意のΔ b \Delta b Δ b に対して∣ Δ z ∣ = ∣ a 1 T Δ b ∣ / 3 ≤ ∥ Δ b ∥ 2 / 3 |\Delta z|=|a_1^{\mathsf T}\Delta b|/3\le\lVert\Delta b\rVert_2/\sqrt3 ∣Δ z ∣ = ∣ a 1 T Δ b ∣/3 ≤ ∥ Δ b ∥ 2 / 3 である。
非零の残差:b 2 : = b 1 + ( 1 , − 2 , 1 ) b_2:=b_1+(1,-2,1) b 2 := b 1 + ( 1 , − 2 , 1 ) とする。( 1 , − 2 , 1 ) (1,-2,1) ( 1 , − 2 , 1 ) はa 1 a_1 a 1 とa 2 a_2 a 2 に直交するから、A x = b 2 Ax=b_2 A x = b 2 の最小二乗解はx = ( 1 , 1 ) x=(1,1) x = ( 1 , 1 ) 、残差はr = ( 1 , − 2 , 1 ) r=(1,-2,1) r = ( 1 , − 2 , 1 ) である。δ : = 10 − 8 \delta:=10^{-8} δ := 1 0 − 8 とし、Δ A \Delta A Δ A の第1 1 1 列を0 0 0 、第2 2 2 列をδ ( 1 , − 2 , 1 ) \delta(1,-2,1) δ ( 1 , − 2 , 1 ) とすると、∥ Δ A ∥ 2 = 6 δ ≤ σ 2 ( A ) / 2 \lVert\Delta A\rVert_2=\sqrt6\,\delta\le\sigma_2(A)/2 ∥ Δ A ∥ 2 = 6 δ ≤ σ 2 ( A ) /2 である。f = − Δ A x = − δ ( 1 , − 2 , 1 ) f=-\Delta Ax=-\delta(1,-2,1) f = − Δ A x = − δ ( 1 , − 2 , 1 ) はA A A の列に直交するからA T f = 0 A^{\mathsf T}f=0 A T f = 0 であり、定理 6.1 (2) の一次の項は( A T A ) − 1 ( Δ A ) T r = ( A T A ) − 1 ( 0 , 6 δ ) = ( 3 δ / ε 2 ) ( − 1 , 1 ) = ( − 0.03 , 0.03 ) (A^{\mathsf T}A)^{-1}(\Delta A)^{\mathsf T}r=(A^{\mathsf T}A)^{-1}(0,6\delta)=(3\delta/\varepsilon^2)(-1,1)=(-0.03,0.03) ( A T A ) − 1 ( Δ A ) T r = ( A T A ) − 1 ( 0 , 6 δ ) = ( 3 δ / ε 2 ) ( − 1 , 1 ) = ( − 0.03 , 0.03 ) である。有理数で計算した厳密な変化は
Δ x = 299999997 10000000003 ( − 1 , 1 ) ≈ ( − 0.029999999691 , 0.029999999691 ) \Delta x=\frac{299999997}{10000000003}(-1,1)\approx(-0.029999999691,\ 0.029999999691) Δ x = 10000000003 299999997 ( − 1 , 1 ) ≈ ( − 0.029999999691 , 0.029999999691 ) である。同じΔ A \Delta A Δ A を残差0 0 0 のb 1 b_1 b 1 に対して加えると、一次の項は0 0 0 であり、厳密な変化はΔ x = 3 10000000003 ( 1 , − 1 ) ≈ ( 3.0 × 10 − 10 , − 3.0 × 10 − 10 ) \Delta x=\frac{3}{10000000003}(1,-1)\approx(3.0\times10^{-10},-3.0\times10^{-10}) Δ x = 10000000003 3 ( 1 , − 1 ) ≈ ( 3.0 × 1 0 − 10 , − 3.0 × 1 0 − 10 ) である。∥ Δ A ∥ 2 / σ 1 ( A ) ≈ 1.0 × 10 − 8 \lVert\Delta A\rVert_2/\sigma_1(A)\approx1.0\times10^{-8} ∥ Δ A ∥ 2 / σ 1 ( A ) ≈ 1.0 × 1 0 − 8 の摂動に対して、b 2 b_2 b 2 では係数の相対変化が0.03 0.03 0.03 であり、b 1 b_1 b 1 では3.0 × 10 − 10 3.0\times10^{-10} 3.0 × 1 0 − 10 である。
注意 6.4. m ≥ n m\ge n m ≥ n とし、A ∈ R m × n A\in\R^{m\times n} A ∈ R m × n の列が一次独立、b ∈ R m b\in\R^m b ∈ R m 、D = diag ( d 1 , … , d n ) D=\operatorname{diag}(d_1,\dots,d_n) D = diag ( d 1 , … , d n ) を正則な対角行列とし、A ~ : = A D \tilde A:=AD A ~ := A D と置く。A ~ z = A ( D z ) \tilde Az=A(Dz) A ~ z = A ( D z ) であるからA ~ \tilde A A ~ の列空間はA A A の列空間に等しく、A ~ \tilde A A ~ の列は一次独立である。命題 1.2 条件 (a) により、z z z がA ~ z = b \tilde Az=b A ~ z = b の最小二乗解であることとx : = D z x:=Dz x := D z がA x = b Ax=b A x = b の最小二乗解であることは同値であり、当てはめた値A ~ z = A x \tilde Az=Ax A ~ z = A x と残差は一致する。係数についてはx j = d j z j x_j=d_jz_j x j = d j z j であるから、成分ごとの相対変化∣ Δ x j ∣ / ∣ x j ∣ = ∣ Δ z j ∣ / ∣ z j ∣ |\Delta x_j|/|x_j|=|\Delta z_j|/|z_j| ∣Δ x j ∣/∣ x j ∣ = ∣Δ z j ∣/∣ z j ∣ (x j ≠ 0 x_j\ne0 x j = 0 )は尺度によらないが、∥ Δ x ∥ 2 / ∥ x ∥ 2 \lVert\Delta x\rVert_2/\lVert x\rVert_2 ∥ Δ x ∥ 2 / ∥ x ∥ 2 と∥ Δ z ∥ 2 / ∥ z ∥ 2 \lVert\Delta z\rVert_2/\lVert z\rVert_2 ∥ Δ z ∥ 2 / ∥ z ∥ 2 (x ≠ 0 x\ne0 x = 0 )は一般に異なる。摂動については、A A A の摂動E E E はA ~ \tilde A A ~ の摂動E ~ : = E D \tilde E:=ED E ~ := E D に対応し、( A + E ) x = b (A+E)x=b ( A + E ) x = b の最小二乗解は( A ~ + E ~ ) z = b (\tilde A+\tilde E)z=b ( A ~ + E ~ ) z = b の最小二乗解のD D D 倍である。E ~ e j = d j E e j \tilde Ee_j=d_jEe_j E ~ e j = d j E e j 、A ~ e j = d j A e j \tilde Ae_j=d_jAe_j A ~ e j = d j A e j であるから、列ごとの相対的な大きさ∥ E e j ∥ 2 / ∥ A e j ∥ 2 \lVert Ee_j\rVert_2/\lVert Ae_j\rVert_2 ∥ E e j ∥ 2 / ∥ A e j ∥ 2 は尺度によらず、∥ E ∥ 2 / ∥ A ∥ 2 \lVert E\rVert_2/\lVert A\rVert_2 ∥ E ∥ 2 / ∥ A ∥ 2 は一般に異なる。
説明変数の単位を変える例として、t = ( 1 , 2 , 3 ) t=(1,2,3) t = ( 1 , 2 , 3 ) (単位は m)と1000 t 1000t 1000 t (単位は mm)を考え、A A A の第i i i 行を( 1 , 1000 t i ) (1,1000t_i) ( 1 , 1000 t i ) 、D : = diag ( 1 , 10 − 3 ) D:=\operatorname{diag}(1,10^{-3}) D := diag ( 1 , 1 0 − 3 ) とすると、A ~ \tilde A A ~ の第i i i 行は( 1 , t i ) (1,t_i) ( 1 , t i ) である。κ 2 ( A ) ≈ 5.72 × 10 3 \kappa_2(A)\approx5.72\times10^3 κ 2 ( A ) ≈ 5.72 × 1 0 3 、κ 2 ( A ~ ) ≈ 6.79 \kappa_2(\tilde A)\approx6.79 κ 2 ( A ~ ) ≈ 6.79 である。∥ A ∥ 2 / ∥ A e 1 ∥ 2 ≈ 2.16 × 10 3 \lVert A\rVert_2/\lVert Ae_1\rVert_2\approx2.16\times10^3 ∥ A ∥ 2 / ∥ A e 1 ∥ 2 ≈ 2.16 × 1 0 3 であるから、∥ E ∥ 2 ≤ ϵ ∥ A ∥ 2 \lVert E\rVert_2\le\epsilon\lVert A\rVert_2 ∥ E ∥ 2 ≤ ϵ ∥ A ∥ 2 を満たす摂動E E E は、第1 1 1 列をそのノルムの2.16 × 10 3 ϵ 2.16\times10^3\epsilon 2.16 × 1 0 3 ϵ 倍まで動かしうる。定理 4.7 (2) のΔ A ′ \Delta A' Δ A ′ は列ごとに∥ Δ a j ′ ∥ 2 ≤ η m , n ′ ∥ a j ∥ 2 \lVert\Delta a'_j\rVert_2\le\eta'_{m,n}\lVert a_j\rVert_2 ∥ Δ a j ′ ∥ 2 ≤ η m , n ′ ∥ a j ∥ 2 を満たすので、Δ A ′ D \Delta A'D Δ A ′ D はA ~ \tilde A A ~ に対して同じ列ごとの評価を満たし、補題 2.1 (1) により∥ Δ A ′ D ∥ 2 ≤ η m , n ′ ∥ A ~ ∥ F \lVert\Delta A'D\rVert_2\le\eta'_{m,n}\lVert\tilde A\rVert_F ∥ Δ A ′ D ∥ 2 ≤ η m , n ′ ∥ A ~ ∥ F である。∥ Δ A ′ D ∥ 2 ≤ σ n ( A ~ ) / 2 \lVert\Delta A'D\rVert_2\le\sigma_n(\tilde A)/2 ∥ Δ A ′ D ∥ 2 ≤ σ n ( A ~ ) /2 ならば、摂動( Δ A ′ D , Δ b ) (\Delta A'D,\Delta b) ( Δ A ′ D , Δ b ) に定理 6.1 (3) をA ~ \tilde A A ~ とκ 2 ( A ~ ) \kappa_2(\tilde A) κ 2 ( A ~ ) で適用してΔ z \Delta z Δ z を評価し、Δ x = D Δ z \Delta x=D\Delta z Δ x = D Δ z により、元の係数x 1 x_1 x 1 (切片)とx 2 = 10 − 3 z 2 x_2=10^{-3}z_2 x 2 = 1 0 − 3 z 2 (mm あたりの傾き)の変化として報告することができる。
7 演習
解答. 0 ≤ u < 1 0\le u<1 0 ≤ u < 1 であるから0 < 1 − u ≤ 1 0<1-u\le1 0 < 1 − u ≤ 1 であり、B k B_k B k の元はすべて正である。
補題 4.2 (1) を示す。1 − u ≤ 1 + δ ≤ 1 + u 1-u\le1+\delta\le1+u 1 − u ≤ 1 + δ ≤ 1 + u であり、( 1 + u ) ( 1 − u ) = 1 − u 2 ≤ 1 (1+u)(1-u)=1-u^2\le1 ( 1 + u ) ( 1 − u ) = 1 − u 2 ≤ 1 から1 + u ≤ ( 1 − u ) − 1 1+u\le(1-u)^{-1} 1 + u ≤ ( 1 − u ) − 1 である。
補題 4.2 (2) を示す。( 1 − u ) j ≤ a ≤ ( 1 − u ) − j (1-u)^j\le a\le(1-u)^{-j} ( 1 − u ) j ≤ a ≤ ( 1 − u ) − j と( 1 − u ) k ≤ a ′ ≤ ( 1 − u ) − k (1-u)^k\le a'\le(1-u)^{-k} ( 1 − u ) k ≤ a ′ ≤ ( 1 − u ) − k の辺々を掛けて( 1 − u ) j + k ≤ a a ′ ≤ ( 1 − u ) − ( j + k ) (1-u)^{j+k}\le aa'\le(1-u)^{-(j+k)} ( 1 − u ) j + k ≤ a a ′ ≤ ( 1 − u ) − ( j + k ) である。正の数の逆数をとると不等号が反転するので( 1 − u ) j ≤ a − 1 ≤ ( 1 − u ) − j (1-u)^j\le a^{-1}\le(1-u)^{-j} ( 1 − u ) j ≤ a − 1 ≤ ( 1 − u ) − j であり、⋅ \sqrt{\cdot} ⋅ は正の数の上で単調増加であるから( 1 − u ) j / 2 ≤ a ≤ ( 1 − u ) − j / 2 (1-u)^{j/2}\le\sqrt a\le(1-u)^{-j/2} ( 1 − u ) j /2 ≤ a ≤ ( 1 − u ) − j /2 である。0 < 1 − u ≤ 1 0<1-u\le1 0 < 1 − u ≤ 1 とj ≤ k j\le k j ≤ k から( 1 − u ) k ≤ ( 1 − u ) j (1-u)^k\le(1-u)^j ( 1 − u ) k ≤ ( 1 − u ) j 、( 1 − u ) − j ≤ ( 1 − u ) − k (1-u)^{-j}\le(1-u)^{-k} ( 1 − u ) − j ≤ ( 1 − u ) − k であり、B j ⊂ B k B_j\subset B_k B j ⊂ B k である。
補題 4.2 (3) を示す。Λ : = ∑ i λ i > 0 \Lambda:=\sum_i\lambda_i>0 Λ := ∑ i λ i > 0 と置くと、λ i ≥ 0 \lambda_i\ge0 λ i ≥ 0 と( 1 − u ) k ≤ a i ≤ ( 1 − u ) − k (1-u)^k\le a_i\le(1-u)^{-k} ( 1 − u ) k ≤ a i ≤ ( 1 − u ) − k からΛ ( 1 − u ) k ≤ ∑ i λ i a i ≤ Λ ( 1 − u ) − k \Lambda(1-u)^k\le\sum_i\lambda_ia_i\le\Lambda(1-u)^{-k} Λ ( 1 − u ) k ≤ ∑ i λ i a i ≤ Λ ( 1 − u ) − k であり、Λ \Lambda Λ で割って主張を得る。
補題 4.2 (4) を示す。§E20.1 補題 4.1 をδ 1 = ⋯ = δ k = − u \delta_1=\dots=\delta_k=-u δ 1 = ⋯ = δ k = − u とε 1 = ⋯ = ε k = 1 \varepsilon_1=\dots=\varepsilon_k=1 ε 1 = ⋯ = ε k = 1 に適用すると、( 1 − u ) k = 1 + θ (1-u)^k=1+\theta ( 1 − u ) k = 1 + θ 、∣ θ ∣ ≤ γ k |\theta|\le\gamma_k ∣ θ ∣ ≤ γ k であるから( 1 − u ) k ≥ 1 − γ k (1-u)^k\ge1-\gamma_k ( 1 − u ) k ≥ 1 − γ k である。ε 1 = ⋯ = ε k = − 1 \varepsilon_1=\dots=\varepsilon_k=-1 ε 1 = ⋯ = ε k = − 1 に適用すると、( 1 − u ) − k = 1 + θ ′ (1-u)^{-k}=1+\theta' ( 1 − u ) − k = 1 + θ ′ 、∣ θ ′ ∣ ≤ γ k |\theta'|\le\gamma_k ∣ θ ′ ∣ ≤ γ k であるから( 1 − u ) − k ≤ 1 + γ k (1-u)^{-k}\le1+\gamma_k ( 1 − u ) − k ≤ 1 + γ k である。したがってa ∈ B k a\in B_k a ∈ B k ならば1 − γ k ≤ a ≤ 1 + γ k 1-\gamma_k\le a\le1+\gamma_k 1 − γ k ≤ a ≤ 1 + γ k である。▨