グラム・シュミットの直交化法

ツールには 10⁻⁸ 離れた三つの列と、それらに対する古典的グラム・シュミットの答えが表示されています。アルゴリズムを手計算で追い、30度の誤差を生む原因となった一つの丸め誤差を見つけ出してください。そして、それぞれの計算手法が実用に耐えなくなる限界がどこにあるかを説明してください。

この問題をツールで開く

1. それぞれの列は、ベクトル (1, 0, 0, 0) の異なる位置に一つだけ ε を挿入したものです。AᵀA はすべての成分が1である 3 × 3 行列に ε²I を足したものであり、その固有値は 3 + ε² が一つと ε² が二つです。したがって特異値は √(3 + ε²) と ε になり、条件数は推定値ではなく厳密に求まります。表の先頭行には 1.7321 × 10⁸ と表示されています。

   $$A=\left(\begin{array}{ccc}1& 1& 1\\ & & \\ \epsilon & 0& 0\\ & & \\ 0& \epsilon & 0\\ & & \\ 0& 0& \epsilon \end{array}\right),\epsilon =10^{-8},\kappa (A)=\frac{\sqrt{3+{\epsilon }^2}}{\epsilon }=1.7321\times 10^8$$

2. ‖a₁‖² は 1 + ε² であり、ε² は 10⁻¹⁶ です。倍精度小数が1を保持するとき、その次に保持できる数は 2.2204 × 10⁻¹⁶ だけ先にあるため、その間隔の半分である 1.1102 × 10⁻¹⁶ 未満の数は、1に足された時点で消えてしまいます。10のマイナス16乗はこの値より小さくなります。平方根は1となり、それで割っても何も変わらないため、q₁ は元の列そのものとして出力され、ε² は計算の記録から完全に消去されてしまいます。この問題における他のすべての丸め誤差は末位の1単位分に過ぎませんが、この丸め誤差だけは30度分の重みを持っているのです。

   $$\Vert a_1{\Vert }^2=1+{\epsilon }^2=1+10^{-16},10^{-16}<1.1102\times 10^{-16}⟹\mathrm{fl}(1+{\epsilon }^2)=1,q_1=a_1$$

3. 次は列 2 ですが、ここより前にはベクトルが一つしかないため、二つの手法はここで一致します。射影の係数は q₁ᵀa₂ = 1 となり、それを引くと (0, −ε, ε, 0) が残り、その長さは ε√2 = 1.4142 × 10⁻⁸ になります。倍精度小数の認識する限り列 2 は長さ1として計算に入ってきましたが、そのうち 1.4142 × 10⁻⁸ を除くすべてがすでに q₁ によって覆われていたことになります。

   $$r_{12}=q_1^{\mathsf{T}}{a}_2=1,v_2=a_2-q_1=(0,-\epsilon ,\epsilon ,0{)}^{\mathsf{T}},\Vert v_2\Vert =\epsilon \sqrt{2}=1.4142\times 10^{-8}$$

4. 列 3 を、まずは古典的の手法で見てみましょう。q₁ に対する係数は再び1です。q₂ に対する係数は q₂ᵀa₃ ですが、q₂ は完全に第2、第3成分に収まっているのに対し、a₃ の ε は第4成分にあるため、その内積は 0 となります。これは q₂ が列 3 から取り除くもの の行に示されている通りです。q₁ だけを引くと (0, −ε, 0, ε) が残り、その長さは再び ε√2 となるため、列 3 の残り も同じく 1.4142 × 10⁻⁸ と表示されます。

   $$r_{13}=q_1^{\mathsf{T}}{a}_3=1,r_{23}=q_2^{\mathsf{T}}{a}_3=0,v_3=a_3-q_1=(0,-\epsilon ,0,\epsilon {)}^{\mathsf{T}},q_3=\frac{1}{\sqrt{2}}(0,-1,0,1{)}^{\mathsf{T}}$$

5. q₂ と q₃ はどちらも、先頭の1を共有する列から q₁ を引くことで作られています。そのため、両者とも第2成分に −1/√2 を持ち、それ以外に共通する成分はありません。それらの内積は半分になります。ツールは 0.5000 および 60.00° と出力しますが、±1/√2 が2進数の分数ではないため、内部に保存されている値は 0.4999999999999999 となり、代数計算が導く半分という値からは末位の1単位だけずれています。これがこの算術計算で起こりうる誤差の大きさです。あの30度の誤差はここから来たのではありません。古典的の手法による直交性の喪失の行には 0.7071 と表示されます。これは √(2 × 0.5²) であり、同じ半分という値が QᵀQ の両隅で一度ずつ数えられた結果です。

   $$q_2^{\mathsf{T}}{q}_3=\frac{1}{2}[(-1)(-1)+(1)(0)]=\frac{1}{2},\arccos \frac{1}{2}=60^{\circ },\sqrt{2\cdot 0.5^2}=0.7071$$

6. 今度は同じ列に対して、修正版を適用してみましょう。q₁ の寄与を引いた後の残りは (0, −ε, 0, ε) となり、この手法では a₃ ではなくこれに対して計算を行います。これと q₂ の内積をとると、ゼロではなく ε/√2 = 7.0711 × 10⁻⁹ が得られます。これを取り除くと残るのは (0, −ε/2, −ε/2, ε) で、その長さは ε√(3/2) = 1.2247 × 10⁻⁸ となり、q₂ᵀq₃ は半分ではなく 1.1102 × 10⁻¹⁶ になります。修正版による直交性の喪失の行には 1.1547 × 10⁻⁸ と表示されています。そこに残っているのは、q₁ が取り除けなかった分だけです。つまり引き算によって q₁ 方向に −ε² が残り、それをわずか ε√2 の長さしかない残りで割ることで、再び ε/√2 に戻ってしまったのです。

   $$r_{23}=q_2^{\mathsf{T}}{v}_3=\frac{\epsilon }{\sqrt{2}}=7.0711\times 10^{-9},v_3\mapsto (0,-\frac{\epsilon }{2},-\frac{\epsilon }{2},\epsilon {)}^{\mathsf{T}},\Vert v_3\Vert =\epsilon \sqrt{\frac{3}{2}}=1.2247\times 10^{-8},q_2^{\mathsf{T}}{q}_3=0$$

7. ツールが提供するすべての列数において、同じ議論が成り立ちます。古典的の手法では、先頭の1を共有する列から q₁ を引くことでそれ以降のすべての q を構成するため、それらの間のすべてのペアの角度が60°になります。その数は (n − 1)(n − 2)/2 組であり、直交性の喪失はそのペア数の半分の平方根に等しくなります。二つの列しかない場合はそのようなペアは存在せず、二つの手法に共通する ε の大きさの傾きだけが残ることになります。

   $$\begin{array}{rrrrr}n& \kappa (A)& \#\{60^{\circ }\}& \Vert Q_C^{\mathsf{T}}{Q}_C-I{\Vert }_F& \Vert Q_M^{\mathsf{T}}{Q}_M-I{\Vert }_F\\ & & & & \\ 2& 1.4142\times 10^8& 0& 1.0000\times 10^{-8}& 1.0000\times 10^{-8}\\ & & & & \\ 3& 1.7321\times 10^8& 1& 0.7071& 1.1547\times 10^{-8}\\ & & & & \\ 4& 2.0000\times 10^8& 3& 1.2247& 1.2247\times 10^{-8}\\ & & & & \\ 5& 2.2361\times 10^8& 6& 1.7321& 1.2649\times 10^{-8}\end{array}$$

8. 古典的の手法の当てはめ直線は κ = 2.35 × 10⁸ で喪失が1に達します。そのため、古典的の曲線はパネルの右端ですでに平坦になっています。一方、修正版の手法では外挿ではなく測定された限界値が得られます。というのも、この曲線は 10⁻⁸ を過ぎると向きを変えるため、そこまで当てはめ直線を伸ばしてしまうと、その下にあるパネルの表示と矛盾してしまうからです。条件数を動かした724点全体において、直交性の喪失が 1.25 × (計算機イプシロン) × κ を超えることはなく、この限界値が1に達するのは κ = 3.60 × 10¹⁵ のときだけです。なお、計算機イプシロンの逆数は 4.50 × 10¹⁵ です。

   $$\begin{array}{rcl}{\log }_{10}\kappa =\frac{16.5755}{1.9801}& =& 8.3710⟹\kappa =2.35\times 10^8\\ & & \\ 1.25\cdot (2.2204\times 10^{-16})\kappa & =& 1⟹\kappa =3.60\times 10^{15}\\ & & \\ (2.2204\times 10^{-16}{)}^{-1}& =& 4.50\times 10^{15}\end{array}$$