# QR 分解：Householder 変換・Givens 回転・Gram–Schmidt の安定性 > QR 分解の存在と一意性を Householder 鏡映で示し、鏡映 1 回の丸め誤差から Householder QR と Givens QR の後退安定性を明示的な定数つきで証明する。古典・修正 Gram–Schmidt が直交性を失う様子を条件数と並べて数値で確かめる。 > https://rikai.mugen-giken.com/mathematics/numerical-linear-algebra/qr-factorization ## 0. この記事の要点 QR 分解は、行列を「ベクトルの長さを変えない部分 $Q$」と「三角形の形をした部分 $R$」の積に分ける方法です。長さを変えない変換だけで計算すれば誤差が膨らみにくいので、最小二乗問題や固有値計算の土台として使われます。この記事では、その作り方（鏡映と回転）、コンピュータの丸め誤差があっても計算結果が「入力をほんの少し動かした問題の厳密な答え」になっていること、そして Gram–Schmidt 法では直交性が崩れうることを順に見ていきます。 - 任意の $m\times n$ 実行列は $A=QR$（$Q$ は直交行列、$R$ は上三角）と書けます。列が一次独立なら、$R$ の対角を正にした薄い QR 分解は一意（薄い QR 分解とは $Q$ を $m\times n$ の長方形のまま使う形）で、その $Q$ の列は Gram–Schmidt の直交化と一致します。 - **Householder 鏡映** $P=I-\beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}$ はベクトルを $\mp\|\boldsymbol{x}\|_2\boldsymbol{e}_1$ に写す対称な直交行列です。符号を $x_1$ と同じ向きに選ぶと $v_1$ の計算で桁落ち（近い数の引き算で有効数字が減ること）が起きません。これを $n$ 回掛ける Householder QR の演算量は $2mn^2-2n^3/3$ です。 - 浮動小数で計算した鏡映は「厳密な鏡映 + 相対 $\gamma_{14p+29}$ の誤差」で、これを積み重ねると、計算した $\widehat R$ について**ある厳密な直交行列 $Q$ で $A+\Delta A=Q\widehat R$、$\|\Delta\boldsymbol{a}_j\|_2\le\gamma_{43mn}\|\boldsymbol{a}_j\|_2$** が成り立ちます（列ごとの後退安定性。計算結果が、入力を少しだけ動かした問題の厳密な答えになっているということです）。Givens 回転による QR も同じ形です。 - 古典 Gram–Schmidt と修正 Gram–Schmidt は厳密な算術では同じ $Q,R$ を与えますが、浮動小数では直交性 $\|I-\widehat Q^{\mathsf{T}}\widehat Q\|$ がそれぞれおよそ $\kappa^2u$、$\kappa u$ まで崩れ、Householder では $u$ の程度にとどまります。後退誤差 $\|A-\widehat Q\widehat R\|/\|A\|$ はどれも $u$ の程度で、違いは直交性に出ます。 ## 1. 動機 [最小二乗問題](/mathematics/numerical-linear-algebra/least-squares#def-least-squares)・固有値計算・[Krylov 部分空間](/mathematics/numerical-linear-algebra/lanczos-and-conjugate-gradient#def-krylov)法は、どれも「列ベクトルを直交化する」操作を中心に持ちます。直交変換は 2 ノルムを保つので、誤差を拡大しません（[行列のノルムと条件数](/mathematics/numerical-linear-algebra/norms-and-conditioning) の ）。ただし、それは**厳密に**直交な行列を掛けた場合の話で、計算機の中の「直交行列」は丸め誤差を含みます。 教科書で最初に習う Gram–Schmidt の直交化は、厳密な算術では申し分のない算法ですが、浮動小数では列がほとんど一次従属のとき直交性が崩れます。これに対して Householder 鏡映や Givens 回転は、「計算した結果が、少しだけ動かしたデータに厳密な直交変換を施した結果になっている」ことを証明できます。この記事の中心は、この後退誤差解析を、丸めの誤差モデルから省略せずに組み立てることです。 ## 2. 準備 この記事では実数の行列を扱います（複素数でも、転置を共役転置に、符号 $\operatorname{sgn}x_1$ を $x_1/|x_1|$ に替えれば同じです）。$\operatorname{sgn}(t)=1$（$t\ge0$）、$-1$（$t<0$）と約束します。直交行列は $Q^{\mathsf{T}}Q=I$ を満たす実正方行列で、 により $QQ^{\mathsf{T}}=I$ かつ $\|Q\boldsymbol{x}\|_2=\|\boldsymbol{x}\|_2$ です。 **浮動小数の約束.** $\mathbb{F}$ は 2 進浮動小数の集合、$u$ は単位丸めで、計算の途中でアンダーフロー・オーバーフローは起きないとします。このとき四則と平方根は により $\mathrm{fl}(x\circ y)=(x\circ y)(1+\delta)$、$|\delta|\le u$ を満たします。$ku<1$ のとき $\gamma_k=ku/(1-ku)$ です（）。内積の丸め誤差は 、$\gamma$ の足し算・掛け算の規則は を使います。以下「$1+\theta$、$|\theta|\le\gamma_k$」の形の因子を $1+\theta_k$ と略記することがあります。 後退誤差解析では定数の正確な値はあまり重要でないので、次の記法が便利です。 $c$ を問題の大きさ（$m,n$）とデータによらない正の整数定数とし、$cku<1$ のとき $\widetilde\gamma_k=\gamma_{ck}=\dfrac{cku}{1-cku}$ とおく。本記事の定理では $c$ の値を明示する（Householder と Givens の QR では $c=43$ で足りる）。 $k\ge1$、$2ku<1$、$|\theta|\le\gamma_k$ とする。 1. $(1+\theta)^{-1}=1+\theta'$ と書くと $|\theta'|\le\gamma_{2k}$。 2. $\sqrt{1+\theta}=1+\theta''$ と書くと $|\theta''|\le\gamma_k$。 $2ku<1$ から $\gamma_k=ku/(1-ku)<1$ なので $1+\theta>0$ です。(1) $\theta'=-\theta/(1+\theta)$ なので $$ |\theta'|\le\frac{\gamma_k}{1-\gamma_k}=\frac{ku/(1-ku)}{(1-2ku)/(1-ku)}=\frac{ku}{1-2ku}\le\frac{2ku}{1-2ku}=\gamma_{2k}. $$ (2) $\sqrt{1+\theta}-1=\theta/(\sqrt{1+\theta}+1)$ で、分母は $1$ 以上なので $|\theta''|\le|\theta|\le\gamma_k$ です。 ## 3. Householder 鏡映と Givens 回転 ここで目指すのは、ベクトル $\boldsymbol{x}$ を「鏡」に映して座標軸 $\boldsymbol{e}_1$ の上に載せること、つまり第 1 成分以外を 0 にすることです。鏡映は長さを変えないので、列ごとにこれを繰り返せば、長さを保ったまま行列の左下を 0 にしていけます。Givens 回転は、鏡の代わりに 2 成分だけを回して 0 を作る道具です。 $\boldsymbol{0}\ne\boldsymbol{v}\in\mathbb{R}^p$ について $P=I-\beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}$、$\beta=2/(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v})$ を **Householder 鏡映**（Householder 変換）という。$\boldsymbol{x}\in\mathbb{R}^p\setminus\{\boldsymbol{0}\}$ に対して $$ \sigma=\|\boldsymbol{x}\|_2,\qquad\boldsymbol{v}=\boldsymbol{x}+\operatorname{sgn}(x_1)\,\sigma\,\boldsymbol{e}_1,\qquad\alpha=-\operatorname{sgn}(x_1)\,\sigma $$ とおいたものを **$\boldsymbol{x}$ の Householder 鏡映**という。 1. 任意の $\boldsymbol{v}\ne\boldsymbol{0}$ について、$P=I-\beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}$ は対称な直交行列で $P^2=I$。$P\boldsymbol{v}=-\boldsymbol{v}$、$\boldsymbol{w}\perp\boldsymbol{v}$ なら $P\boldsymbol{w}=\boldsymbol{w}$。 2. $\boldsymbol{x}\ne\boldsymbol{0}$ の Householder 鏡映では $|v_1|=|x_1|+\sigma\ge\sigma>0$、$\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v}=2\sigma|v_1|$（したがって $\beta=1/(\sigma|v_1|)$）で、$P\boldsymbol{x}=\alpha\boldsymbol{e}_1$。 (1) $P^{\mathsf{T}}=P$ は明らかです。$\beta\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v}=2$ を使うと $$ P^{\mathsf{T}}P=P^2=I-2\beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}+\beta^2\boldsymbol{v}(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v})\boldsymbol{v}^{\mathsf{T}}=I-2\beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}+2\beta\boldsymbol{v}\boldsymbol{v}^{\mathsf{T}}=I $$ で、$P$ は直交行列です。$P\boldsymbol{v}=\boldsymbol{v}-\beta(\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v})\boldsymbol{v}=-\boldsymbol{v}$、$\boldsymbol{v}^{\mathsf{T}}\boldsymbol{w}=0$ なら $P\boldsymbol{w}=\boldsymbol{w}$ です。 (2) $v_1=x_1+\operatorname{sgn}(x_1)\sigma$ で、$x_1$ と $\operatorname{sgn}(x_1)\sigma$ は同符号（または $x_1=0$）なので $|v_1|=|x_1|+\sigma$ です。$\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v}=\sigma^2+2|x_1|\sigma+\sigma^2=2\sigma(\sigma+|x_1|)=2\sigma|v_1|$。また $\boldsymbol{v}^{\mathsf{T}}\boldsymbol{x}=\sigma^2+|x_1|\sigma=\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v}/2$ なので $\beta\boldsymbol{v}^{\mathsf{T}}\boldsymbol{x}=1$、$P\boldsymbol{x}=\boldsymbol{x}-\boldsymbol{v}=-\operatorname{sgn}(x_1)\sigma\boldsymbol{e}_1$ です。
x=(3,4) Px=(−5,0) e₁ の方向 鏡（v に直交）
$P$ は $\boldsymbol{v}$ に直交する超平面に関する鏡映です。$\boldsymbol{x}$ を $\|\boldsymbol{x}\|_2\boldsymbol{e}_1$ と $-\|\boldsymbol{x}\|_2\boldsymbol{e}_1$ のどちらに写すかは自由ですが、$\boldsymbol{v}=\boldsymbol{x}-\operatorname{sgn}(x_1)\sigma\boldsymbol{e}_1$ を選ぶと、$\boldsymbol{x}$ が $\boldsymbol{e}_1$ の方向に近いとき $v_1=x_1-\operatorname{sgn}(x_1)\sigma$ が近い数の引き算になり、桁落ちで $\boldsymbol{v}$ の相対誤差が大きくなります。定義の符号はこれを避ける選び方で、後の の証明で $|x_1|+\sigma$ の相対誤差が $\sigma$ の相対誤差以下になることに効いています。 $1\le i $G=G(i,k;c,s)$ は直交行列で、$G\boldsymbol{y}$ は $\boldsymbol{y}$ の第 $i$・第 $k$ 成分だけを $(cy_i+sy_k,\ -sy_i+cy_k)$ に替える。$c=a/r$、$s=b/r$ のとき、$y_i=a$、$y_k=b$ なら $(G\boldsymbol{y})_i=r$、$(G\boldsymbol{y})_k=0$。 成分の式は行列の積の定義どおりです。$G^{\mathsf{T}}G$ は第 $i$・$k$ 行列以外で単位行列に一致し、その $2\times2$ の部分は $$ \begin{pmatrix}c&-s\\s&c\end{pmatrix}\begin{pmatrix}c&s\\-s&c\end{pmatrix}=\begin{pmatrix}c^2+s^2&0\\0&c^2+s^2\end{pmatrix}=I $$ です。$ca+sb=(a^2+b^2)/r=r$、$-sa+cb=(-ab+ab)/r=0$ です。 Householder 鏡映は 1 回で多くの成分を 0 にし、Givens 回転は 1 成分ずつ、しかも 2 行しか動かさずに 0 にします。密な行列では前者、すでに 0 の多い行列（対角の 1 つ下より下がすべて 0 の Hessenberg 行列や、帯行列）では後者が向いています。 ## 4. QR 分解の存在・一意性と演算量 前節の鏡映を第 1 列、第 2 列、…と順に使って $A$ の左下を 0 にしていくと、上三角の $R$ が残り、掛けた鏡映の積が $Q$ になります。次の定理は、これがどんな行列でもできること、そして列が一次独立なら Gram–Schmidt の直交化と同じ結果になることを述べます。 1. 任意の $A\in\mathbb{R}^{m\times n}$ について、直交行列 $Q\in\mathbb{R}^{m\times m}$ と上三角（$i>j$ で $r_{ij}=0$）の $R\in\mathbb{R}^{m\times n}$ で $A=QR$ となるものがある。$Q$ は高々 $\min(m,n)$ 個の Householder 鏡映の積にとれる。 2. $A$ の列 $\boldsymbol{a}_1,\ldots,\boldsymbol{a}_n$ が一次独立（したがって $m\ge n$）なら、列が正規直交な $Q_1\in\mathbb{R}^{m\times n}$ と対角成分が正の上三角行列 $R_1\in\mathbb{R}^{n\times n}$ で $A=Q_1R_1$ となるもの（**薄い QR 分解**）がただ 1 組存在する。$Q_1$ の第 $k$ 列は $\boldsymbol{a}_1,\ldots,\boldsymbol{a}_n$ を Gram–Schmidt の直交化で正規直交化した第 $k$ のベクトルであり、$r_{kk}=\boldsymbol{q}_k^{\mathsf{T}}\boldsymbol{a}_k$ である。 (1) $A^{(1)}=A$ とし、$k=1,\ldots,\min(m,n)$ について次を行います。$A^{(k)}$ の第 $k$ 列の第 $k$〜$m$ 成分を $\boldsymbol{x}\in\mathbb{R}^{m-k+1}$ とし、$\boldsymbol{x}\ne\boldsymbol{0}$ なら $\widetilde P_k$ を $\boldsymbol{x}$ の Householder 鏡映、$\boldsymbol{x}=\boldsymbol{0}$ なら $\widetilde P_k=I$ として、$P_k=\operatorname{diag}(I_{k-1},\widetilde P_k)$、$A^{(k+1)}=P_kA^{(k)}$ とおきます。 により $P_k$ は対称な直交行列で、$A^{(k+1)}$ の第 $k$ 列の第 $k+1$ 成分以降は 0 です。$P_k$ は第 $1$〜$k-1$ 行を変えず、第 $j0$ です。したがって $\boldsymbol{q}_1,\ldots,\boldsymbol{q}_n$ は の条件 1・2・3 を満たし、同定理の一意性の主張により Gram–Schmidt の出力と一致します。$R_1=Q_1^{\mathsf{T}}A$ も決まるので、薄い QR 分解は一意です。 **Householder QR の算法.** 実際の計算では $Q$ を作らず、鏡映 $P_k$ をベクトル $\boldsymbol{v}_k$ と $\beta_k$ で持ち、残りの列 $B$（第 $k$〜$m$ 行、第 $k+1$〜$n$ 列）に $\widetilde P_kB=B-(\beta_k\boldsymbol{v}_k)(\boldsymbol{v}_k^{\mathsf{T}}B)$ と作用させます。第 $k$ 列そのものは計算せず、対角に $\alpha_k$ を書き、その下に 0 を書きます（ の 2）。 $m\ge n$ とし、浮動小数の加減乗除を 1 回 1 演算と数える。上の算法で鏡映を残りの列に作用させる演算は、$p\times q$ の行列 $B$ に $B-(\beta\boldsymbol{v})(\boldsymbol{v}^{\mathsf{T}}B)$ を $4pq$ 回で計算するとき合計 $$ \sum_{k=1}^{n}4(m-k+1)(n-k)=2mn^2-\frac23n^3-2mn+2n^2-\frac43n $$ 回であり、鏡映のベクトルの計算は合計 $O(mn)$ 回である。したがって全体は $2mn^2-\frac23n^3+O(mn)$ である。 $\boldsymbol{w}^{\mathsf{T}}=\boldsymbol{v}^{\mathsf{T}}B$ は $q$ 個の長さ $p$ の内積で $q(2p-1)$ 回、$\beta\boldsymbol{v}$ は $p$ 回、$B-(\beta\boldsymbol{v})\boldsymbol{w}^{\mathsf{T}}$ は $2pq$ 回で、主要部は $4pq$ です（$p$ と $q$ の 1 次の項は $O(p+q)$ で、$k$ について足しても $O(mn)$）。第 $k$ 段では $p=m-k+1$、$q=n-k$ です。$j=n-k$（$0\le j\le n-1$）とおくと $p=m-n+1+j$ で $$ \sum_{j=0}^{n-1}4(m-n+1+j)j=2(m-n+1)n(n-1)+\frac23n(n-1)(2n-1) $$ を展開すると主張の式です。$m\ge n$ なので $-2mn+2n^2-\frac43n=O(mn)$ です。$\boldsymbol{v}_k$・$\beta_k$ の計算は第 $k$ 段で $O(m-k+1)$ 回、合計 $O(mn)$ です。 $m=n$ なら $\frac43n^3$ で、[LU 分解](/mathematics/numerical-linear-algebra/lu-and-pivoting#def-lu)（$\frac23n^3$）の 2 倍です。$Q$ を明示的に作るには、さらに鏡映を単位行列に作用させる $O(m^2n)$ の手間がかかります。最小二乗問題では $Q^{\mathsf{T}}\boldsymbol{b}$ だけが要るので、$Q$ を作らずに鏡映を $\boldsymbol{b}$ に作用させます。 $H\in\mathbb{R}^{n\times n}$ を上 Hessenberg 行列（$i>j+1$ で $h_{ij}=0$）とします。$k=1,\ldots,n-1$ の順に、第 $k$・$k+1$ 行に Givens 回転 $G_k=G(k,k+1;c_k,s_k)$ を掛けて $(k+1,k)$ 成分を 0 にします。第 $k$ 段の時点で第 $k$ 列の第 $k+2$ 行以降はもとから 0、それより左の列の第 $k$・$k+1$ 行は 0 なので、回転は第 $k$〜$n$ 列の 2 行だけを動かし、1 列あたり乗算 4 回・加算 2 回の 6 演算です。合計 $\sum_{k=1}^{n-1}6(n-k+1)=3n^2+O(n)$ 演算で、$G_{n-1}\cdots G_1H=R$ が上三角になります。Householder 鏡映をそのまま使うと $\frac43n^3$ 演算かかるので、0 の構造を使う Givens 回転が有利です。この計算は QR 法（固有値計算）の 1 反復の中身です。 ## 5. 丸め誤差の解析と後退安定性 ### 5.1. 鏡映 1 回の丸め誤差 $\boldsymbol{x}\in\mathbb{F}^p$、$\boldsymbol{x}\ne\boldsymbol{0}$ に対し、Householder 鏡映を次のように計算します（$\widehat{\ }$ は計算値）。 $$ \widehat s=\mathrm{fl}(\boldsymbol{x}^{\mathsf{T}}\boldsymbol{x}),\quad\widehat\sigma=\mathrm{fl}\bigl(\sqrt{\widehat s}\bigr),\quad\widehat v_1=\mathrm{fl}\bigl(x_1+\operatorname{sgn}(x_1)\widehat\sigma\bigr),\quad\widehat v_i=x_i\ (i\ge2),\quad\widehat\beta=\mathrm{fl}\bigl(1/\mathrm{fl}(\widehat\sigma|\widehat v_1|)\bigr),\quad\widehat\alpha=-\operatorname{sgn}(x_1)\widehat\sigma. $$ 内積 $\widehat s$ は先頭から順に足します。$\boldsymbol{b}\in\mathbb{F}^p$ への作用は $$ \widehat t=\mathrm{fl}(\widehat{\boldsymbol{v}}^{\mathsf{T}}\boldsymbol{b}),\quad\widehat\tau=\mathrm{fl}(\widehat\beta\,\widehat t),\quad\widehat y_i=\mathrm{fl}\bigl(b_i-\mathrm{fl}(\widehat\tau\,\widehat v_i)\bigr) $$ で計算します。比べる相手は、同じ $\boldsymbol{x}$ から**厳密に**作った鏡映 $P$（）です。 次の補題は、計算した鏡映が厳密な鏡映からどれだけずれるかを測ります。証明は、各演算に標準モデルの誤差因子 $(1+\delta)$ を掛けて追いかけ、現れる $\gamma$ の添字を足し合わせていく作業です。添字 $14p+29$ はその合計で、数そのものより「$p$ の 1 次式で済む」ことが大事です。 $p\ge1$、$43pu<1/2$ とし、上の計算でアンダーフロー・オーバーフローが起きないとする。$\boldsymbol{v},\beta,\alpha,P$ を $\boldsymbol{x}$ から厳密に作った量とする。 1. $\widehat v_1=v_1(1+\theta_v)$、$|\theta_v|\le\gamma_{p+2}$、$\widehat\beta=\beta(1+\theta_\beta)$、$|\theta_\beta|\le\gamma_{4p+8}$。 2. 任意の $\boldsymbol{b}\in\mathbb{F}^p$ について $\|\widehat{\boldsymbol{y}}-P\boldsymbol{b}\|_2\le\gamma_{14p+29}\|\boldsymbol{b}\|_2$。 3. $\|\widehat\alpha\boldsymbol{e}_1-P\boldsymbol{x}\|_2\le\gamma_{p+1}\|\boldsymbol{x}\|_2$。 以下、$\gamma$ の和・積の評価には を、$|\delta|\le u\le\gamma_1$ の形の丸めには を断らずに使います。現れる添字はすべて $43p$ 以下なので、規則の前提 $ku<1$ は満たされます。 **$\widehat\sigma$.** により $\widehat s=\sum_ix_i^2(1+\theta^{(i)})$、$|\theta^{(i)}|\le\gamma_p$ です。各項が非負なので $\widehat s=\sigma^2(1+\theta_s)$、$\theta_s=\sum_ix_i^2\theta^{(i)}/\sigma^2$、$|\theta_s|\le\gamma_p$ と書けます。 の 2 で $\sqrt{1+\theta_s}=1+\eta$、$|\eta|\le\gamma_p$ なので、$\widehat\sigma=\sigma(1+\eta)(1+\delta)=\sigma(1+\theta_\sigma)$、$|\theta_\sigma|\le\gamma_{p+1}$ です。 **(3).** の 2 より $P\boldsymbol{x}=\alpha\boldsymbol{e}_1$ で、$|\widehat\alpha-\alpha|=|\widehat\sigma-\sigma|\le\gamma_{p+1}\sigma=\gamma_{p+1}\|\boldsymbol{x}\|_2$ です。 **(1).** $\operatorname{sgn}(x_1)x_1=|x_1|$ なので $\widehat v_1=\operatorname{sgn}(x_1)\bigl(|x_1|+\sigma+\sigma\theta_\sigma\bigr)(1+\delta_1)$ です。$|x_1|+\sigma+\sigma\theta_\sigma=(|x_1|+\sigma)(1+\theta')$、$\theta'=\sigma\theta_\sigma/(|x_1|+\sigma)$ とおくと $|\theta'|\le|\theta_\sigma|\le\gamma_{p+1}$ で（ここで符号の選び方が効く）、$v_1=\operatorname{sgn}(x_1)(|x_1|+\sigma)$ から $\widehat v_1=v_1(1+\theta')(1+\delta_1)=v_1(1+\theta_v)$、$|\theta_v|\le\gamma_{p+2}$ です。次に $\mathrm{fl}(\widehat\sigma|\widehat v_1|)=\widehat\sigma|\widehat v_1|(1+\delta_2)$、$\widehat\beta=(1+\delta_3)/(\widehat\sigma|\widehat v_1|(1+\delta_2))$ で、$\beta=1/(\sigma|v_1|)$（ の 2）を使うと $$ \widehat\beta=\beta\,(1+\delta_3)(1+\delta_2)^{-1}(1+\theta_\sigma)^{-1}(1+\theta_v)^{-1} $$ です。逆数の因子は、$(1+\delta_2)^{-1}=1-\delta_2/(1+\delta_2)$ で $|\delta_2|/(1-u)\le\gamma_1$ なので $1+\theta_1$、残りの 2 つは の 1 でそれぞれ $1+\theta_{2p+2}$、$1+\theta_{2p+4}$ となり、積は $1+\theta_\beta$、$|\theta_\beta|\le\gamma_{1+1+(2p+2)+(2p+4)}=\gamma_{4p+8}$ です。 **(2).** 厳密な値を $T=\boldsymbol{v}^{\mathsf{T}}\boldsymbol{b}$、$\boldsymbol{z}=\beta T\boldsymbol{v}$ とおくと $P\boldsymbol{b}=\boldsymbol{b}-\boldsymbol{z}$ です。 (i) $\widehat{\boldsymbol{v}}$ は $\boldsymbol{v}$ と第 1 成分だけが違い、$|\widehat{\boldsymbol{v}}|\le(1+\gamma_{p+2})|\boldsymbol{v}|$ です。内積の丸め誤差から $\widehat t=\widehat{\boldsymbol{v}}^{\mathsf{T}}\boldsymbol{b}+e$、$|e|\le\gamma_p|\widehat{\boldsymbol{v}}|^{\mathsf{T}}|\boldsymbol{b}|$ で、$\widehat{\boldsymbol{v}}^{\mathsf{T}}\boldsymbol{b}=T+\theta_vv_1b_1$ なので $\widehat t=T+e'$、 $$ |e'|\le\bigl(\gamma_{p+2}+\gamma_p(1+\gamma_{p+2})\bigr)|\boldsymbol{v}|^{\mathsf{T}}|\boldsymbol{b}|\le\gamma_{2p+2}\|\boldsymbol{v}\|_2\|\boldsymbol{b}\|_2 $$ です（最後は Cauchy–Schwarz の不等式）。 (ii) $\widehat\tau=\widehat\beta\,\widehat t(1+\delta_4)=\beta(1+\theta_a)(T+e')$、$|\theta_a|\le\gamma_{4p+9}$ です。$\widehat q_i:=\mathrm{fl}(\widehat\tau\widehat v_i)=\widehat\tau\widehat v_i(1+\delta_{5,i})$ で、$\widehat v_i(1+\delta_{5,i})=v_i+f_i$ と書くと $|f_i|\le\gamma_{p+3}|v_i|$、したがって $\|\boldsymbol{f}\|_2\le\gamma_{p+3}\|\boldsymbol{v}\|_2$ です。よって $\widehat{\boldsymbol{q}}=\beta(1+\theta_a)(T+e')(\boldsymbol{v}+\boldsymbol{f})$ です。 (iii) $(1+\theta_a)(T+e')(\boldsymbol{v}+\boldsymbol{f})-T\boldsymbol{v}$ を展開すると、$\theta_a,e',\boldsymbol{f}$ の少なくとも 1 つを含む 7 項の和です。$|T|\le\|\boldsymbol{v}\|_2\|\boldsymbol{b}\|_2$、$|e'|\le\gamma_{2p+2}\|\boldsymbol{v}\|_2\|\boldsymbol{b}\|_2$、$\|\boldsymbol{f}\|_2\le\gamma_{p+3}\|\boldsymbol{v}\|_2$ を各項に使って足すと $$ \|\widehat{\boldsymbol{q}}-\boldsymbol{z}\|_2\le\beta\|\boldsymbol{v}\|_2^2\|\boldsymbol{b}\|_2\bigl[(1+\gamma_{4p+9})(1+\gamma_{2p+2})(1+\gamma_{p+3})-1\bigr]\le2\gamma_{7p+14}\|\boldsymbol{b}\|_2 $$ です（$\beta\|\boldsymbol{v}\|_2^2=2$）。 (iv) $\widehat y_i=(b_i-\widehat q_i)(1+\delta_{6,i})$ なので $\|\widehat{\boldsymbol{y}}-(\boldsymbol{b}-\widehat{\boldsymbol{q}})\|_2\le u\|\boldsymbol{b}-\widehat{\boldsymbol{q}}\|_2$ です。$\boldsymbol{b}-\widehat{\boldsymbol{q}}=P\boldsymbol{b}-(\widehat{\boldsymbol{q}}-\boldsymbol{z})$ と $\|P\boldsymbol{b}\|_2=\|\boldsymbol{b}\|_2$（ の 1）から $$ \|\widehat{\boldsymbol{y}}-P\boldsymbol{b}\|_2\le u\|\boldsymbol{b}\|_2+(1+u)\|\widehat{\boldsymbol{q}}-\boldsymbol{z}\|_2\le\bigl(u+(1+u)\gamma_{14p+28}\bigr)\|\boldsymbol{b}\|_2\le\gamma_{14p+29}\|\boldsymbol{b}\|_2 $$ です（$2\gamma_{7p+14}\le\gamma_{14p+28}$）。 定数 $14p+29$ は評価を簡単にするための粗いもので、丁寧に数えれば小さくできます。大事なのは、誤差が $\boldsymbol{b}$ の大きさに比例し、$p$ の 1 次でしか増えないことです。 ### 5.2. 直交変換の積み重ね $Q_1,\ldots,Q_r\in\mathbb{R}^{p\times p}$ を直交行列、$\delta\ge0$ とし、$\boldsymbol{y}_0=\boldsymbol{b}$ と $\boldsymbol{y}_1,\ldots,\boldsymbol{y}_r$ が $\|\boldsymbol{y}_k-Q_k\boldsymbol{y}_{k-1}\|_2\le\delta\|\boldsymbol{y}_{k-1}\|_2$ を満たすとする。このとき $$ \|\boldsymbol{y}_r-Q_r\cdots Q_1\boldsymbol{b}\|_2\le\bigl((1+\delta)^r-1\bigr)\|\boldsymbol{b}\|_2 $$ であり、$\delta\le\gamma_k$、$rku<1$ なら右辺は $\gamma_{rk}\|\boldsymbol{b}\|_2$ 以下である。 $\|Q_k\boldsymbol{y}_{k-1}\|_2=\|\boldsymbol{y}_{k-1}\|_2$ なので $\|\boldsymbol{y}_k\|_2\le(1+\delta)\|\boldsymbol{y}_{k-1}\|_2$、帰納的に $\|\boldsymbol{y}_k\|_2\le(1+\delta)^k\|\boldsymbol{b}\|_2$ です。$\boldsymbol{d}_k=\boldsymbol{y}_k-Q_k\cdots Q_1\boldsymbol{b}$ とおくと $\boldsymbol{d}_k=(\boldsymbol{y}_k-Q_k\boldsymbol{y}_{k-1})+Q_k\boldsymbol{d}_{k-1}$ なので $$ \|\boldsymbol{d}_k\|_2\le\delta(1+\delta)^{k-1}\|\boldsymbol{b}\|_2+\|\boldsymbol{d}_{k-1}\|_2,\qquad\|\boldsymbol{d}_r\|_2\le\sum_{k=1}^r\delta(1+\delta)^{k-1}\|\boldsymbol{b}\|_2=\bigl((1+\delta)^r-1\bigr)\|\boldsymbol{b}\|_2 $$ です（$\boldsymbol{d}_0=\boldsymbol{0}$）。後半は、 の 2 を繰り返して $(1+\gamma_k)^r\le1+\gamma_{rk}$ とすればよいです。 ### 5.3. Householder QR と Givens QR の後退安定性 **算法（浮動小数）.** $A\in\mathbb{F}^{m\times n}$、$m\ge n$ とし、$k=1,\ldots,n$ について、現在の行列の第 $k$ 列の第 $k$〜$m$ 成分を $\boldsymbol{x}$ とします。$\boldsymbol{x}=\boldsymbol{0}$ なら何もしません。そうでなければ 5.1 の方法で $\widehat{\boldsymbol{v}},\widehat\beta,\widehat\alpha$ を計算し、第 $k+1$〜$n$ 列の第 $k$〜$m$ 行に鏡映を作用させ、第 $k$ 列の第 $k$ 成分を $\widehat\alpha$、第 $k+1$〜$m$ 成分を $0$ とします。最後の行列を $\widehat R$ とします。各段で使った $\boldsymbol{x}$（計算値）から**厳密に**作った鏡映を $\widetilde P_k$（$\boldsymbol{x}=\boldsymbol{0}$ なら $I$）、$P_k=\operatorname{diag}(I_{k-1},\widetilde P_k)$ とします。 各段の誤差は 5.1 の補題で、段を重ねたときの誤差の増え方は 5.2 の補題で抑えられます。次の定理は、この 2 つを列ごとに組み合わせたものです。 上の算法でアンダーフロー・オーバーフローが起きず、$86mnu<1$ とする。$Q=P_1P_2\cdots P_n$ は直交行列で、 1. $A+\Delta A=Q\widehat R$、$\|\Delta\boldsymbol{a}_j\|_2\le\gamma_{43mn}\|\boldsymbol{a}_j\|_2$（$j=1,\ldots,n$、$\Delta\boldsymbol{a}_j$ は $\Delta A$ の第 $j$ 列）。とくに $\|\Delta A\|_F\le\gamma_{43mn}\|A\|_F$。 2. $\boldsymbol{b}\in\mathbb{F}^m$ に同じ鏡映を順に作用させた計算値 $\widehat{\boldsymbol{c}}$ について、$\widehat{\boldsymbol{c}}=Q^{\mathsf{T}}(\boldsymbol{b}+\Delta\boldsymbol{b})$、$\|\Delta\boldsymbol{b}\|_2\le\gamma_{43mn}\|\boldsymbol{b}\|_2$。 $P_k$ は により対称な直交行列なので、$Q$ は直交行列で $Q^{\mathsf{T}}=P_n\cdots P_1$ です。第 $k$ 段の $\boldsymbol{x}$ の長さは $p=m-k+1\le m$ で、$14p+29\le43m$、$p+1\le43m$ です。また $86mnu<1$ から $43pu\le43mu<1/2$ なので が使えます。 (1) 第 $j$ 列を固定し、第 $k$ 段の後の第 $j$ 列（計算値）を $\boldsymbol{y}_k$、$\boldsymbol{y}_0=\boldsymbol{a}_j$ とします。 - $k の 2 から $\|\boldsymbol{y}_k-P_k\boldsymbol{y}_{k-1}\|_2\le\gamma_{43m}\|\boldsymbol{y}_{k-1}\|_2$。 - $k=j$: 第 $j$ 列が $\boldsymbol{x}$ そのもので、 の 3 から同じ評価（$\boldsymbol{x}=\boldsymbol{0}$ なら誤差 0）。 - $k>j$: 第 $j$ 列は変わらず、第 $k$ 行以降は 0 なので $P_k\boldsymbol{y}_{k-1}=\boldsymbol{y}_{k-1}$、誤差 0。 （$r=n$、$\delta=\gamma_{43m}$、$43mnu<1$）により $\widehat{\boldsymbol{r}}_j:=\boldsymbol{y}_n$ は $\|\widehat{\boldsymbol{r}}_j-Q^{\mathsf{T}}\boldsymbol{a}_j\|_2\le\gamma_{43mn}\|\boldsymbol{a}_j\|_2$ を満たします。$\Delta\boldsymbol{a}_j:=Q\widehat{\boldsymbol{r}}_j-\boldsymbol{a}_j=Q(\widehat{\boldsymbol{r}}_j-Q^{\mathsf{T}}\boldsymbol{a}_j)$ とおけば、$Q$ は 2 ノルムを保つので主張の評価です。Frobenius ノルムの評価は列ごとの評価の 2 乗和です。 (2) 各段で の 2 を使い、(1) と同様に を適用して $\Delta\boldsymbol{b}:=Q\widehat{\boldsymbol{c}}-\boldsymbol{b}$ とおけばよいです。 この定理は「$\widehat R$ は $A$ の厳密な QR 分解の $R$ に近い」とは言っていません。言っているのは「$\widehat R$ は、各列を相対 $\gamma_{43mn}$ だけ動かした行列の、厳密な QR 分解の $R$ である」ことです。前進誤差（計算結果そのものと厳密な答えとのずれ）は QR 分解の条件数に依存します。また $Q$ は計算で現れる行列ではなく、解析のための厳密な直交行列です。計算機の中で $Q$ を作った場合については次の系が成り立ちます。 の仮定のもとで、$I_m$ の各列に鏡映を $P_n,\ldots,P_1$ の順に作用させて計算した $\widehat Q$ は $\widehat Q=Q(I+\Delta)$、$\|\Delta\|_F\le\sqrt m\,\gamma_{43mn}$ を満たし、 $$ \|\widehat Q^{\mathsf{T}}\widehat Q-I\|_2\le2\sqrt m\,\gamma_{43mn}+m\gamma_{43mn}^2. $$ $\boldsymbol{e}_i$ に $P_n,\ldots,P_1$ の順に作用させる計算に、 の 2 と （$Q_k$ を $P_{n+1-k}$ とする）を使うと、第 $i$ 列は $P_1\cdots P_n\boldsymbol{e}_i+\boldsymbol{d}_i=Q(\boldsymbol{e}_i+Q^{\mathsf{T}}\boldsymbol{d}_i)$、$\|\boldsymbol{d}_i\|_2\le\gamma_{43mn}$ です。$\Delta$ の第 $i$ 列を $Q^{\mathsf{T}}\boldsymbol{d}_i$ とすれば $\|\Delta\|_F^2\le m\gamma_{43mn}^2$。$\widehat Q^{\mathsf{T}}\widehat Q-I=\Delta+\Delta^{\mathsf{T}}+\Delta^{\mathsf{T}}\Delta$ で、 の 3 から $\|\Delta\|_2\le\|\Delta\|_F$、劣乗法性から $\|\Delta^{\mathsf{T}}\Delta\|_2\le\|\Delta\|_2^2$ です。 **Givens QR の算法（浮動小数）.** $k=1,\ldots,n$、$i=k+1,\ldots,m$ の順に、現在の $(i,k)$ 成分 $b$ が 0 でなければ、$(k,k)$ 成分 $a$ と合わせて $\widehat r=\mathrm{fl}\bigl(\sqrt{\mathrm{fl}(\mathrm{fl}(a^2)+\mathrm{fl}(b^2))}\bigr)$、$\widehat c=\mathrm{fl}(a/\widehat r)$、$\widehat s=\mathrm{fl}(b/\widehat r)$ を計算し、第 $k+1$〜$n$ 列の第 $k$・$i$ 行を $(x,y)\mapsto(\mathrm{fl}(\mathrm{fl}(\widehat cx)+\mathrm{fl}(\widehat sy)),\ \mathrm{fl}(\mathrm{fl}(\widehat cy)-\mathrm{fl}(\widehat sx)))$ で置き換え、$(k,k)$ 成分を $\widehat r$、$(i,k)$ 成分を $0$ とします。比べる相手は、計算値の $a,b$ から厳密に作った回転 $G=G(k,i;a/r,b/r)$ です。 上の算法でアンダーフロー・オーバーフローが起きず、$36mnu<1$ とする。使った回転（厳密）の積を $Q^{\mathsf{T}}=G_N\cdots G_1$（$N\le mn$）とすると、$A+\Delta A=Q\widehat R$、$\|\Delta\boldsymbol{a}_j\|_2\le\gamma_{18mn}\|\boldsymbol{a}_j\|_2$。 **回転 1 回.** $\mathrm{fl}(a^2)+\mathrm{fl}(b^2)$ の丸めを含めて $\mathrm{fl}(\mathrm{fl}(a^2)+\mathrm{fl}(b^2))=r^2(1+\theta_2)$（$a^2(1+\delta)+b^2(1+\delta')=r^2(1+\eta)$、$|\eta|\le u$ と和の丸め）なので、 の 2 と平方根の丸めで $\widehat r=r(1+\theta_r)$、$|\theta_r|\le\gamma_3$。 の 1 で $\widehat c=c(1+\delta)(1+\theta_r)^{-1}=c(1+\theta_c)$、$|\theta_c|\le\gamma_7$、同様に $\widehat s=s(1+\theta_s)$、$|\theta_s|\le\gamma_7$ です。作用させた成分は $\widehat x'=cx(1+\theta_9)+sy(1+\theta_9')$、$\widehat y'=cy(1+\theta_9'')-sx(1+\theta_9''')$ と書けるので、誤差ベクトルの成分は $\gamma_9(|c||x|+|s||y|)$、$\gamma_9(|s||x|+|c||y|)$ 以下です。Cauchy–Schwarz の不等式で $(|c||x|+|s||y|)^2\le(c^2+s^2)(x^2+y^2)$ なので、誤差の 2 ノルムは $\sqrt2\gamma_9\|(x,y)\|_2\le\gamma_{18}\|(x,y)\|_2$ です。0 にした列では、厳密な回転の結果 $(r,0)$ との差は $|\widehat r-r|\le\gamma_3r\le\gamma_{18}\|(a,b)\|_2$ です。どちらの場合も、2 行以外の成分は変わらず、列全体で $\|\widehat{\boldsymbol{y}}-G\boldsymbol{y}\|_2\le\gamma_{18}\|\boldsymbol{y}\|_2$ です。 **積み重ね.** 第 $k$ 段の回転は第 $k$ 行と第 $i>k$ 行だけを動かし、第 $j により各 $G_\ell$ は直交行列です。（$\delta=\gamma_{18}$、$r=N$、$18Nu\le18mnu<1$）から、 の証明と同じく $\|\widehat{\boldsymbol{r}}_j-Q^{\mathsf{T}}\boldsymbol{a}_j\|_2\le\gamma_{18mn}\|\boldsymbol{a}_j\|_2$ で、$\Delta\boldsymbol{a}_j=Q\widehat{\boldsymbol{r}}_j-\boldsymbol{a}_j$ とおけばよいです。 の前提（$2\cdot3u<1$）や $\gamma_{18}$ の定義の前提は $36mnu<1$ に含まれます。 互いに重ならない行の組への回転は同時に（順序によらず）行えるので、回転を $O(m+n)$ 個の「段」にまとめると、1 列が受ける誤差の積み重ねは段の数で数えられ、評価は $\widetilde\gamma_{m+n}$ にまで改善されます（Higham, §19.6）。上の定理はこの改善を使わない素朴な版です。 ## 6. Gram–Schmidt の直交化と直交性の喪失 列が一次独立な $A=(\boldsymbol{a}_1,\ldots,\boldsymbol{a}_n)\in\mathbb{R}^{m\times n}$ について、$k=1,\ldots,n$ の順に次で $\boldsymbol{q}_k$ と $R$ の第 $k$ 列を作る。 - **古典 Gram–Schmidt（CGS）**: $r_{ik}=\boldsymbol{q}_i^{\mathsf{T}}\boldsymbol{a}_k$（$i CGS は係数をすべて元の $\boldsymbol{a}_k$ との内積で計算し、MGS は「それまでの方向を取り除いた残り」との内積で計算します。 厳密な算術では、CGS と MGS は同じ $\boldsymbol{q}_1,\ldots,\boldsymbol{q}_n$ と $R$ を与え、それは の 2 の薄い QR 分解（対角が正）に一致する。 $k$ についての帰納法で、両者の $\boldsymbol{q}_1,\ldots,\boldsymbol{q}_{k-1}$ が一致して正規直交であると仮定します。MGS で $\boldsymbol{w}_k^{(i)}=\boldsymbol{a}_k-\sum_{j の漸化式そのものなので、$\boldsymbol{w}_k\ne\boldsymbol{0}$ で $\boldsymbol{q}_1,\ldots,\boldsymbol{q}_k$ は正規直交です。作り方から $\boldsymbol{a}_k=\sum_{i0$ で、 の 2 の一意性から主張が従います。 浮動小数では両者は違う結果を与えます。違いが目に見える最小の例は 3 列です（2 列では CGS と MGS は同じ演算を行います。演習を参照）。 次の例で丸められるのは $\varepsilon^2=10^{-16}$ です。倍精度では $1$ に足しても桁に収まらず消えるので、$1+\varepsilon^2$ は $1$ になります。この消えた分が、直交性の喪失のもとになります。 $\varepsilon>0$ を $\mathrm{fl}(1+\varepsilon^2)=1$ となるほど小さく（倍精度で $\varepsilon=10^{-8}$）とり、 $$ A=\begin{pmatrix}1&1&1\\\varepsilon&0&0\\0&\varepsilon&0\\0&0&\varepsilon\end{pmatrix} $$ とします（列は一次独立）。計算では $\|\boldsymbol{a}_1\|_2=\sqrt{1+\varepsilon^2}$ が $1$ に丸められ、$\widehat{\boldsymbol{q}}_1=(1,\varepsilon,0,0)^{\mathsf{T}}$ です。$\boldsymbol{q}_1^{\mathsf{T}}\boldsymbol{a}_2=1$ で $\boldsymbol{w}_2=(0,-\varepsilon,\varepsilon,0)^{\mathsf{T}}$、$\widehat{\boldsymbol{q}}_2=(0,-1,1,0)^{\mathsf{T}}/\sqrt2$。 - CGS: $\boldsymbol{q}_1^{\mathsf{T}}\boldsymbol{a}_3=1$、$\boldsymbol{q}_2^{\mathsf{T}}\boldsymbol{a}_3=0$ なので $\boldsymbol{w}_3=(0,-\varepsilon,0,\varepsilon)^{\mathsf{T}}$、$\widehat{\boldsymbol{q}}_3=(0,-1,0,1)^{\mathsf{T}}/\sqrt2$ で、$\widehat{\boldsymbol{q}}_2^{\mathsf{T}}\widehat{\boldsymbol{q}}_3=1/2$。 - MGS: $\boldsymbol{q}_1$ を除いた $\boldsymbol{w}_3^{(2)}=(0,-\varepsilon,0,\varepsilon)^{\mathsf{T}}$ との内積 $\boldsymbol{q}_2^{\mathsf{T}}\boldsymbol{w}_3^{(2)}=\varepsilon/\sqrt2$ を使うので $\boldsymbol{w}_3=(0,-\varepsilon/2,-\varepsilon/2,\varepsilon)^{\mathsf{T}}$、$\widehat{\boldsymbol{q}}_3=(0,-1,-1,2)^{\mathsf{T}}/\sqrt6$ で、$\widehat{\boldsymbol{q}}_2^{\mathsf{T}}\widehat{\boldsymbol{q}}_3=0$。 どちらも $\widehat{\boldsymbol{q}}_1^{\mathsf{T}}\widehat{\boldsymbol{q}}_2=-\varepsilon/\sqrt2$ 程度の誤差は残ります。CGS では $\boldsymbol{q}_2$ の方向の成分が $\boldsymbol{a}_3$ の段階で 0 と計算され、$\boldsymbol{q}_1$ を引いた後に初めて現れる $\boldsymbol{q}_2$ 方向の成分が取り除かれません。倍精度で実行すると、CGS の $\widehat Q^{\mathsf{T}}\widehat Q$ の $(2,3)$ 成分は $0.5$、MGS では $0.0$（表示 10 桁）でした。この行列は Läuchli の例として知られる形です（Björck, Ch. 2）。 **MGS の直交性（主張のみ、定理環境にしない）.** 列の数を $n$、$\kappa_2(A)u$ が十分小さいとき、MGS の $\widehat Q$ は $\|I-\widehat Q^{\mathsf{T}}\widehat Q\|_2\le c_1(m,n)\,\kappa_2(A)\,u$ を満たします（Björck の教科書 Ch. 2）。CGS では喪失が $\kappa_2(A)^2u$ の程度まで悪化しうること、$\kappa_2(A)^2u$ が十分小さければその程度で抑えられることが知られています（Giraud–Langou–Rozložník 2005 の概説）。丸めの解析が長いのでこの記事では証明せず、以後の証明でも使いません。 **再直交化.** CGS の各段で $\boldsymbol{w}_k$ に対してもう一度同じ射影を行う（CGS2、"twice is enough"）と、$\kappa_2(A)u$ が十分小さい範囲で直交性は $u$ の程度に回復します（同じ概説を参照）。CGS2 は内積をまとめて行列ベクトル積で計算できるので、並列計算機では MGS より速いことが多く、Krylov 部分空間法の Arnoldi 過程でも使われます。 **MGS と Householder の関係.** $A$ の上に $n\times n$ の零行列を積んだ $\binom{O}{A}$ に Householder QR を施すと、浮動小数でも MGS と同じ数値が現れます（Björck–Paige 1992）。この同値性により、MGS の $\widehat R$ と $\widehat Q$ の組は「直交性は崩れても、最小二乗問題を解くには後退安定」に使えることが示されます。 特異値を $1$ から $1/\kappa$ まで等比に並べた $60\times20$ の行列 $A=U\Sigma V^{\mathsf{T}}$（$U,V$ は乱数行列の QR 分解の $Q$）で、CGS・MGS・CGS2・`numpy.linalg.qr`（LAPACK の Householder QR）の直交性 $\|I-\widehat Q^{\mathsf{T}}\widehat Q\|_2$ を比べます。 ```python # Python 3.12.10 / NumPy 2.5.2 ／乱数種 20261006 import numpy as np def cgs(A, passes=1): # 古典 Gram-Schmidt（passes=2 で再直交化） m, n = A.shape; Q = np.zeros((m, n)); R = np.zeros((n, n)) for k in range(n): v = A[:, k].copy() for _ in range(passes): c = Q[:, :k].T @ v; v -= Q[:, :k] @ c; R[:k, k] += c R[k, k] = np.linalg.norm(v); Q[:, k] = v / R[k, k] return Q, R def mgs(A): # 修正 Gram-Schmidt m, n = A.shape; V = A.copy(); Q = np.zeros((m, n)) R = np.zeros((n, n)) for k in range(n): R[k, k] = np.linalg.norm(V[:, k]); Q[:, k] = V[:, k] / R[k, k] R[k, k+1:] = Q[:, k] @ V[:, k+1:] V[:, k+1:] -= np.outer(Q[:, k], R[k, k+1:]) return Q, R rng = np.random.default_rng(20261006) m, n = 60, 20 U, _ = np.linalg.qr(rng.standard_normal((m, n))) V, _ = np.linalg.qr(rng.standard_normal((n, n))) loss = lambda Q: np.linalg.norm(np.eye(n) - Q.T @ Q, 2) print("kappa CGS MGS CGS2 Householder") for e in range(2, 13, 2): A = U @ np.diag(np.logspace(0, -e, n)) @ V.T Qs = [cgs(A)[0], mgs(A)[0], cgs(A, 2)[0], np.linalg.qr(A)[0]] print(f"1e{e:<3d} " + " ".join(f"{loss(Q):.1e}" for Q in Qs)) ``` ```text kappa CGS MGS CGS2 Householder 1e2 6.3e-14 9.3e-15 4.6e-16 6.8e-16 1e4 8.0e-10 3.3e-13 4.6e-16 9.9e-16 1e6 1.3e-05 4.9e-11 5.4e-16 8.7e-16 1e8 2.3e-01 1.3e-09 4.8e-16 1.3e-15 1e10 3.9e+00 3.6e-07 6.0e-16 7.9e-16 1e12 5.4e+00 2.6e-06 7.0e-16 9.9e-16 ``` 続けて、後退誤差 $\|A-\widehat Q\widehat R\|_2/\|A\|_2$ と、mpmath の 50 桁で計算した QR 分解に対する Householder の $\widehat Q$・$\widehat R$ の前進誤差（列の符号をそろえて比べる）を見ます。 ```python # （上のブロックの続き。cgs・mgs・U・V・m・n を使う） import mpmath as mp # mpmath 1.3.0 mp.mp.dps = 50 nrm = lambda X: np.linalg.norm(X, 2) print("kappa CGS MGS CGS2 Househ. Q fwd R fwd") for e in (2, 6, 10): A = U @ np.diag(np.logspace(0, -e, n)) @ V.T QRs = [cgs(A), mgs(A), cgs(A, 2), np.linalg.qr(A)] Q, R = QRs[3] Qm, Rm = mp.qr(mp.matrix(A.tolist())) # 50 桁の QR Qe = np.array(Qm.tolist(), dtype=float)[:, :n] Re = np.array(Rm.tolist(), dtype=float)[:n, :] d = np.sign(np.diag(R)) * np.sign(np.diag(Re)) # 符号をそろえる bw = [nrm(A - Qi @ Ri) / nrm(A) for Qi, Ri in QRs] fq = nrm(Q * d - Qe) fr = nrm(d[:, None] * R - Re) / nrm(Re) print(f"1e{e:<3d} " + " ".join(f"{x:.1e}" for x in bw + [fq, fr])) ``` ```text kappa CGS MGS CGS2 Househ. Q fwd R fwd 1e2 9.1e-17 9.5e-17 9.3e-17 4.2e-16 7.1e-15 2.4e-16 1e6 7.2e-17 7.5e-17 8.0e-17 3.2e-16 1.9e-11 3.0e-16 1e10 7.3e-17 6.3e-17 7.9e-17 4.3e-16 1.0e-07 2.1e-16 ``` **検算.** (a) 既知値: $A$ の 2 ノルムの条件数は構成から厳密に $\kappa$ です。 が予言するとおり、Householder の直交性は $\kappa$ によらず $10^{-15}$ の程度です。(b) 別精度・別手法: 4 通りの直交化を並べ、Householder の $\widehat Q,\widehat R$ は mpmath の 50 桁の QR 分解と比べました。$\widehat R$ の前進誤差は $u$ の程度、$\widehat Q$ の前進誤差は $\kappa u$ の程度（$\kappa u=1.1\times10^{-14},10^{-10},10^{-6}$ に対し $7\times10^{-15},2\times10^{-11},1\times10^{-7}$）で、$\widehat Q$ は「真の $Q$ からは $\kappa u$ ずれるが、それ自体はほぼ直交」です。これは の「近いデータの厳密な分解」という形の主張と合っています。(c) 次数: $\kappa$ を $10^2$ 倍するごとに、CGS の喪失はおよそ $10^4$ 倍（$\kappa^2u$ に比例、$\kappa=10^8$ 以降は $O(1)$ で飽和）、MGS はおよそ $10^2$ 倍（$\kappa u$ に比例）になり、CGS2 と Householder は変わりません。後退誤差はどの方法でも $u$ の程度で、違いは直交性にだけ出ます。 ## 7. 演習 $\boldsymbol{x}=(3,4)^{\mathsf{T}}$ の Householder 鏡映 $\boldsymbol{v},\beta,P$ を求め、$P\boldsymbol{x}$ を計算せよ。また $\boldsymbol{x}=(1,10^{-9})^{\mathsf{T}}$ について、逆の符号 $\boldsymbol{v}'=\boldsymbol{x}-\operatorname{sgn}(x_1)\|\boldsymbol{x}\|_2\boldsymbol{e}_1$ を倍精度で計算すると何が起きるか説明せよ。 $\sigma=5$、$\boldsymbol{v}=(8,4)^{\mathsf{T}}$、$\boldsymbol{v}^{\mathsf{T}}\boldsymbol{v}=80=2\sigma|v_1|$、$\beta=1/40$ で $$ P=I-\frac1{40}\begin{pmatrix}64&32\\32&16\end{pmatrix}=\begin{pmatrix}-0.6&-0.8\\-0.8&0.6\end{pmatrix},\qquad P\boldsymbol{x}=(-5,0)^{\mathsf{T}} $$ で、 の 2 のとおり $\alpha=-5$ です。後半: $\|\boldsymbol{x}\|_2=\sqrt{1+10^{-18}}$ は倍精度で $1$ に丸められるので $v_1'=1-1=0$ となり、$\boldsymbol{v}'=(0,10^{-9})^{\mathsf{T}}$ です。厳密には $v_1'=1-\sqrt{1+10^{-18}}\approx-5\times10^{-19}$ なので、$v_1'$ の相対誤差は 100% です。この $\boldsymbol{v}'$ から作った鏡映は第 2 成分の符号を反転するだけで、$\boldsymbol{x}$ を $\boldsymbol{e}_1$ の方向に写しません（$(1,-10^{-9})$ になる）。定義の符号では $v_1=2$ で桁落ちは起きません。 $A=\begin{pmatrix}1&1\\0&0\end{pmatrix}$ について、$R$ の対角成分が非負の QR 分解 $A=QR$ が 2 通り以上あることを示し、 の 2 で列の一次独立性が必要な理由を説明せよ。 $Q=I$、$R=A$ と、$Q'=\operatorname{diag}(1,-1)$、$R'=Q'^{\mathsf{T}}A=A$ はどちらも $A=QR$ で、$R$ の対角は $(1,0)$ で非負です。$Q\ne Q'$ なので一意ではありません。$r_{22}=0$ のとき $\boldsymbol{q}_2$ は $\boldsymbol{a}_1,\boldsymbol{a}_2$ のどちらにも現れず、$\boldsymbol{q}_1$ に直交する任意の単位ベクトルを選べます。一意性の証明は $r_{kk}>0$（すなわち $\boldsymbol{a}_k\notin\operatorname{span}(\boldsymbol{a}_1,\ldots,\boldsymbol{a}_{k-1})$）を使っていました。 列が 2 本（$n=2$）のとき、CGS と MGS は浮動小数でもまったく同じ演算を行い、同じ結果を与えることを示せ。$n=3$ では最初にどこで違いが生じるか。 $k=1$ はどちらも $\boldsymbol{q}_1=\mathrm{fl}(\boldsymbol{a}_1/\|\boldsymbol{a}_1\|_2)$ です。$k=2$ では CGS が $r_{12}=\mathrm{fl}(\widehat{\boldsymbol{q}}_1^{\mathsf{T}}\boldsymbol{a}_2)$、$\boldsymbol{w}_2=\mathrm{fl}(\boldsymbol{a}_2-r_{12}\widehat{\boldsymbol{q}}_1)$ を計算し、MGS は $\boldsymbol{w}_2^{(1)}=\boldsymbol{a}_2$ に対して同じ式を計算します。したがって演算は同一です。$n=3$ では $k=3$ の $r_{23}$ で、CGS は $\widehat{\boldsymbol{q}}_2^{\mathsf{T}}\boldsymbol{a}_3$、MGS は $\widehat{\boldsymbol{q}}_2^{\mathsf{T}}\widehat{\boldsymbol{w}}_3^{(2)}$ を計算し、ここで初めて違いが出ます。 はまさにこの $r_{23}$ の違い（$0$ と $\varepsilon/\sqrt2$）が直交性を左右する例です。 ## 参考文献 - N. J. Higham, *Accuracy and Stability of Numerical Algorithms*, 2nd ed., SIAM, 2002 — Ch. 19（Householder・Givens QR の後退誤差解析、Gram–Schmidt）。 - G. H. Golub and C. F. Van Loan, *Matrix Computations*, 4th ed., Johns Hopkins University Press, 2013 — §5.1–5.2。 - L. N. Trefethen and D. Bau III, *Numerical Linear Algebra*, SIAM, 1997 — Lectures 7–10, 16。 - Å. Björck, *Numerical Methods for Least Squares Problems*, SIAM, 1996 — Ch. 2（Gram–Schmidt の丸め誤差、Läuchli の例）。 - Å. Björck and C. C. Paige, "Loss and recapture of orthogonality in the modified Gram–Schmidt algorithm", *SIAM J. Matrix Anal. Appl.* 13 (1992), 176–190. - L. Giraud, J. Langou, and M. Rozložník, "The loss of orthogonality in the Gram–Schmidt orthogonalization process", *Comput. Math. Appl.* 50 (2005), 1069–1075.