行列を「きれいな形」に分解できたら、計算がどれほど楽になるでしょうか。たとえば連立方程式 $\bm{A}\bm{x} = \bm{b}$ を解くとき、$\bm{A}$ がもし直交行列と上三角行列の積で書けたら、後退代入だけで解が求まります。また、最小二乗法で回帰直線をフィッティングする際にも、正規方程式を直接解くより数値的に安定な方法が手に入ります。
この「直交行列 $\times$ 上三角行列」への分解こそが QR分解（QR decomposition）です。QR分解は線形代数の中でも特に実用性が高く、以下のような場面で中心的な役割を果たします。
- 最小二乗法: 過決定な連立方程式の数値的に安定な求解
- 固有値計算: QRアルゴリズムによる行列の固有値・固有ベクトルの反復計算
- 信号処理: ビームフォーミングや適応フィルタリングにおける直交分解
- 機械学習: 主成分分析（PCA）やグラム行列の正則化
本記事の内容
- QR分解の直感的な理解と数学的定義
- グラム・シュミット直交化法の理論と導出
- 修正グラム・シュミット法による数値安定性の改善
- ハウスホルダー変換によるQR分解
- Pythonによる実装と最小二乗法への応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
QR分解とは — 直感的な理解
まず、QR分解を直感的に理解しましょう。3次元空間に3本の線形独立なベクトルが与えられた場面を想像してください。これらのベクトルは互いに斜めを向いているかもしれませんが、「1本目はそのままの方向を保ち、2本目は1本目と直交するように修正し、3本目は1本目・2本目の両方と直交するように修正する」という手順で、互いに直交する3本のベクトルを作ることができます。
これが グラム・シュミットの直交化法 であり、QR分解の核心です。元の行列 $\bm{A}$ の列ベクトルを直交化して得られる直交行列 $\bm{Q}$ と、直交化の過程で生じる係数をまとめた上三角行列 $\bm{R}$ に分解するのです。
日常的なアナロジーで言えば、QR分解は「傾いた本棚の本を1冊ずつ垂直に立て直す」作業に似ています。最初の本（第1列）を基準にまっすぐ立て、2冊目は1冊目と直角になるように立て、3冊目は1冊目と2冊目の両方と直角になるように立てます。この「立て直し方のレシピ」が行列 $\bm{R}$ に記録されます。
歴史的には、グラム・シュミット直交化法は19世紀のヨルゲン・ペダーセン・グラム（1883年）とエルハルト・シュミット（1907年）にその名を由来します。彼らは関数空間における直交系の構成を研究していましたが、その離散版が行列のQR分解として現代の数値線形代数の基盤となりました。
QR分解の直感を掴んだところで、次に数学的な定義を見ていきましょう。
QR分解の数学的定義
QR分解とは、$m \times n$ 行列 $\bm{A}$（$m \geq n$）を直交行列（またはユニタリ行列）$\bm{Q}$ と上三角行列 $\bm{R}$ の積に分解することです。
$$ \begin{equation} \bm{A} = \bm{Q}\bm{R} \end{equation} $$
ここで各行列のサイズと性質を整理します。
| 行列 | サイズ | 性質 |
|---|---|---|
| $\bm{A}$ | $m \times n$ | 分解対象の行列（$m \geq n$、列フルランクを仮定） |
| $\bm{Q}$ | $m \times m$（フルQR）または $m \times n$（薄いQR） | 直交行列: $\bm{Q}^T\bm{Q} = \bm{I}$ |
| $\bm{R}$ | $m \times n$（フルQR）または $n \times n$（薄いQR） | 上三角行列 |
なぜこの形に分解するのでしょうか。直交行列 $\bm{Q}$ は逆行列が転置と一致する（$\bm{Q}^{-1} = \bm{Q}^T$）ため、逆行列計算が不要になります。また上三角行列 $\bm{R}$ の連立方程式は後退代入で $O(n^2)$ で解けます。LU分解と比較して、QR分解は数値的に安定（ピボット選択が不要）であり、長方形行列にも適用できるという利点があります。
フルQR分解と薄いQR分解
実用上は、$m > n$ の場合に2種類のQR分解を区別します。
フルQR分解: $\bm{A} = \bm{Q}_{\text{full}} \bm{R}_{\text{full}}$ で、$\bm{Q}_{\text{full}}$ は $m \times m$ の直交行列、$\bm{R}_{\text{full}}$ は $m \times n$ の行列です。$\bm{R}_{\text{full}}$ の下側 $m – n$ 行はゼロになります。
薄い（reduced）QR分解: $\bm{A} = \bm{Q}_{\text{thin}} \bm{R}_{\text{thin}}$ で、$\bm{Q}_{\text{thin}}$ は $m \times n$（列が正規直交）、$\bm{R}_{\text{thin}}$ は $n \times n$ の上三角行列です。計算コストとメモリの面から、多くの応用ではこちらが使われます。
QR分解の一意性
$\bm{A}$ が列フルランクのとき、$\bm{R}$ の対角成分を全て正とする条件を付ければ、薄いQR分解は一意に定まります。これは、直交化の過程でベクトルの「向き」を揃える自由度を固定することに対応します。
それでは、QR分解を実際にどう計算するのか、最も基本的なグラム・シュミット法から見ていきましょう。
グラム・シュミット直交化法
アルゴリズムの導出
グラム・シュミット法は、行列 $\bm{A}$ の列ベクトル $\bm{a}_1, \bm{a}_2, \ldots, \bm{a}_n$ から正規直交ベクトル $\bm{q}_1, \bm{q}_2, \ldots, \bm{q}_n$ を逐次的に構成するアルゴリズムです。
核心となるアイデアは射影の除去です。新しいベクトル $\bm{a}_k$ から、既に構成した正規直交ベクトル $\bm{q}_1, \ldots, \bm{q}_{k-1}$ の方向成分を取り除くと、残りは $\bm{q}_1, \ldots, \bm{q}_{k-1}$ と直交する成分だけになります。
ステップ1: 最初のベクトルは単純に正規化します。
$$ \bm{q}_1 = \frac{\bm{a}_1}{\|\bm{a}_1\|} $$
ステップ2: $\bm{a}_2$ から $\bm{q}_1$ 方向の成分を除去し、正規化します。
$$ \tilde{\bm{q}}_2 = \bm{a}_2 – (\bm{q}_1^T \bm{a}_2)\bm{q}_1, \quad \bm{q}_2 = \frac{\tilde{\bm{q}}_2}{\|\tilde{\bm{q}}_2\|} $$
ここで $(\bm{q}_1^T \bm{a}_2)\bm{q}_1$ は $\bm{a}_2$ の $\bm{q}_1$ 方向への射影です。これを引くことで $\bm{q}_1$ と直交する成分だけが残ります。
一般のステップ $k$: $\bm{a}_k$ から $\bm{q}_1, \ldots, \bm{q}_{k-1}$ の全ての方向成分を除去します。
$$ \tilde{\bm{q}}_k = \bm{a}_k – \sum_{j=1}^{k-1} (\bm{q}_j^T \bm{a}_k)\bm{q}_j, \quad \bm{q}_k = \frac{\tilde{\bm{q}}_k}{\|\tilde{\bm{q}}_k\|} $$
QR分解との関係
グラム・シュミット法の各ステップを行列の言葉で書き直すと、QR分解が自然に現れます。
$\bm{a}_k$ の式を変形すると
$$ \bm{a}_k = \sum_{j=1}^{k-1} (\bm{q}_j^T \bm{a}_k)\bm{q}_j + \|\tilde{\bm{q}}_k\|\bm{q}_k $$
ここで $r_{jk} = \bm{q}_j^T \bm{a}_k$（$j < k$）および $r_{kk} = \|\tilde{\bm{q}}_k\|$ と定義すると
$$ \bm{a}_k = \sum_{j=1}^{k} r_{jk} \bm{q}_j $$
これを全ての列についてまとめると
$$ \bm{A} = \bm{Q}\bm{R} $$
となります。$\bm{R}$ は上三角行列であることに注目してください。$r_{jk}$ は $j \leq k$ でのみ定義されるので、$j > k$ の成分はゼロです。これは直交化の過程で、各ベクトルは自分より前のベクトルの情報しか使わないことを反映しています。
具体例: 3×3行列のQR分解
具体的な数値例でグラム・シュミット法を実行してみましょう。
$$ \bm{A} = \begin{pmatrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 0 & 1 & 1 \end{pmatrix} $$
ステップ1: $\bm{a}_1 = (1, 1, 0)^T$ を正規化します。
$$ r_{11} = \|\bm{a}_1\| = \sqrt{1^2 + 1^2 + 0^2} = \sqrt{2} $$
$$ \bm{q}_1 = \frac{1}{\sqrt{2}}(1, 1, 0)^T $$
ステップ2: $\bm{a}_2 = (1, 0, 1)^T$ から $\bm{q}_1$ 成分を除去します。
$\bm{q}_1$ への射影係数を計算すると
$$ r_{12} = \bm{q}_1^T \bm{a}_2 = \frac{1}{\sqrt{2}}(1 \cdot 1 + 1 \cdot 0 + 0 \cdot 1) = \frac{1}{\sqrt{2}} $$
射影成分を引くと
$$ \tilde{\bm{q}}_2 = \bm{a}_2 – r_{12}\bm{q}_1 = \begin{pmatrix} 1 \\ 0 \\ 1 \end{pmatrix} – \frac{1}{\sqrt{2}} \cdot \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix} = \begin{pmatrix} 1/2 \\ -1/2 \\ 1 \end{pmatrix} $$
正規化すると
$$ r_{22} = \|\tilde{\bm{q}}_2\| = \sqrt{1/4 + 1/4 + 1} = \sqrt{3/2} $$
$$ \bm{q}_2 = \frac{1}{\sqrt{3/2}}\begin{pmatrix} 1/2 \\ -1/2 \\ 1 \end{pmatrix} = \frac{1}{\sqrt{6}}\begin{pmatrix} 1 \\ -1 \\ 2 \end{pmatrix} $$
ステップ3: $\bm{a}_3 = (0, 1, 1)^T$ から $\bm{q}_1$ 成分と $\bm{q}_2$ 成分を除去します。
各射影係数を計算すると
$$ r_{13} = \bm{q}_1^T \bm{a}_3 = \frac{1}{\sqrt{2}}(0 + 1 + 0) = \frac{1}{\sqrt{2}} $$
$$ r_{23} = \bm{q}_2^T \bm{a}_3 = \frac{1}{\sqrt{6}}(0 – 1 + 2) = \frac{1}{\sqrt{6}} $$
射影成分を引くと
$$ \tilde{\bm{q}}_3 = \bm{a}_3 – r_{13}\bm{q}_1 – r_{23}\bm{q}_2 = \begin{pmatrix} 0 \\ 1 \\ 1 \end{pmatrix} – \frac{1}{2}\begin{pmatrix} 1 \\ 1 \\ 0 \end{pmatrix} – \frac{1}{6}\begin{pmatrix} 1 \\ -1 \\ 2 \end{pmatrix} = \begin{pmatrix} -2/3 \\ 2/3 \\ 2/3 \end{pmatrix} $$
正規化して
$$ r_{33} = \|\tilde{\bm{q}}_3\| = \sqrt{4/9 + 4/9 + 4/9} = \frac{2\sqrt{3}}{3} $$
$$ \bm{q}_3 = \frac{1}{\sqrt{3}}\begin{pmatrix} -1 \\ 1 \\ 1 \end{pmatrix} $$
最終的に、$\bm{Q}$ と $\bm{R}$ は
$$ \bm{Q} = \begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{6} & -1/\sqrt{3} \\ 1/\sqrt{2} & -1/\sqrt{6} & 1/\sqrt{3} \\ 0 & 2/\sqrt{6} & 1/\sqrt{3} \end{pmatrix}, \quad \bm{R} = \begin{pmatrix} \sqrt{2} & 1/\sqrt{2} & 1/\sqrt{2} \\ 0 & \sqrt{3/2} & 1/\sqrt{6} \\ 0 & 0 & 2\sqrt{3}/3 \end{pmatrix} $$
$\bm{Q}^T\bm{Q} = \bm{I}$ を確認すれば、直交性が保たれていることがわかります。
数値例でQR分解の具体的な計算が理解できたところで、次に古典的グラム・シュミット法の問題点と、それを改善する修正グラム・シュミット法を見ていきましょう。
修正グラム・シュミット法 — 数値安定性の改善
古典的グラム・シュミット法の問題
理論上は美しい古典的グラム・シュミット（CGS）法ですが、浮動小数点演算で実行すると深刻な問題が生じます。丸め誤差の蓄積により、得られる $\bm{Q}$ の直交性が大きく崩れてしまうのです。
直感的に言えば、CGSでは $\bm{a}_k$ から全ての射影を「一度に」引いています。しかし、$\bm{q}_1$ 成分の除去で丸め誤差が生じると、その誤差が $\bm{q}_2$ 成分の除去にも影響し、誤差が累積します。特に $\bm{A}$ の列ベクトルがほぼ線形従属（条件数が大きい）な場合、この問題は深刻です。
条件数 $\kappa(\bm{A})$ が大きい行列に対して、CGSで得られる $\bm{Q}$ の直交性の誤差は
$$ \|\bm{Q}^T\bm{Q} – \bm{I}\| \sim O(\varepsilon_{\text{mach}} \kappa(\bm{A})^2) $$
となることが知られています（$\varepsilon_{\text{mach}} \approx 10^{-16}$ は機械イプシロン）。条件数が $10^8$ 程度の行列では、直交性が完全に失われます。
修正グラム・シュミット（MGS）法
修正グラム・シュミット（MGS）法は、この問題を大幅に改善する巧妙な変更です。アイデアは、射影を「一度に全部」引くのではなく、「1つずつ順番に」引くことです。
CGSでは
$$ \tilde{\bm{q}}_k = \bm{a}_k – \sum_{j=1}^{k-1} (\bm{q}_j^T \bm{a}_k)\bm{q}_j $$
と全ての射影を元のベクトル $\bm{a}_k$ に対して計算しますが、MGSでは
$$ \bm{v}_k^{(0)} = \bm{a}_k $$
$$ \bm{v}_k^{(j)} = \bm{v}_k^{(j-1)} – (\bm{q}_j^T \bm{v}_k^{(j-1)})\bm{q}_j, \quad j = 1, 2, \ldots, k-1 $$
$$ \bm{q}_k = \frac{\bm{v}_k^{(k-1)}}{\|\bm{v}_k^{(k-1)}\|} $$
と各ステップで「最新の」ベクトルに対して射影を計算します。数学的にはCGSと同じ結果を与えますが、数値的には大きな違いがあります。
MGSにおける直交性の誤差は
$$ \|\bm{Q}^T\bm{Q} – \bm{I}\| \sim O(\varepsilon_{\text{mach}} \kappa(\bm{A})) $$
とCGSの $\kappa(\bm{A})^2$ から $\kappa(\bm{A})$ に改善されます。
なぜ改善されるのでしょうか。MGSでは $\bm{q}_j$ 成分を除去した後のベクトル $\bm{v}_k^{(j)}$ に対して次の射影を計算するため、前のステップの丸め誤差で残った $\bm{q}_j$ 方向の成分がある程度再除去されるのです。一方CGSでは、全ての射影を元の $\bm{a}_k$ に対して計算するため、このような自己修正メカニズムが働きません。
MGSでも条件数が非常に大きい行列では不十分な場合があります。そこで、さらに安定な方法として次に紹介するハウスホルダー変換が使われます。
ハウスホルダー変換によるQR分解
ハウスホルダー反射の定義
ハウスホルダー変換（Householder transformation）は、ベクトルを鏡面反射させる直交変換です。鏡に映すように、あるベクトルを特定の方向に「ひっくり返す」操作と考えることができます。
単位ベクトル $\bm{v}$（$\|\bm{v}\| = 1$）に対して、ハウスホルダー行列は次のように定義されます。
$$ \begin{equation} \bm{H} = \bm{I} – 2\bm{v}\bm{v}^T \end{equation} $$
この行列は $\bm{v}$ に直交する超平面に関する鏡面反射を表します。重要な性質として
- 直交性: $\bm{H}^T\bm{H} = \bm{I}$（直交行列）
- 対称性: $\bm{H}^T = \bm{H}$
- 対合性: $\bm{H}^2 = \bm{I}$（2回反射すると元に戻る）
が成り立ちます。対合性は $\bm{H}^2 = (\bm{I} – 2\bm{v}\bm{v}^T)(\bm{I} – 2\bm{v}\bm{v}^T) = \bm{I} – 4\bm{v}\bm{v}^T + 4\bm{v}(\bm{v}^T\bm{v})\bm{v}^T = \bm{I}$ と $\bm{v}^T\bm{v} = 1$ から確認できます。
QR分解への適用
ハウスホルダー変換を使ったQR分解のアイデアは、列ベクトルを順次 $\bm{e}_1 = (1, 0, \ldots, 0)^T$ 方向に「倒す」ことです。
行列 $\bm{A}$ の第1列 $\bm{a}_1$ を $\alpha \bm{e}_1$（$\alpha = \pm\|\bm{a}_1\|$）に写すハウスホルダーベクトルは
$$ \bm{u}_1 = \bm{a}_1 – \alpha \bm{e}_1, \quad \bm{v}_1 = \frac{\bm{u}_1}{\|\bm{u}_1\|} $$
数値安定性のため、$\alpha = -\text{sign}(a_{11})\|\bm{a}_1\|$ とします（キャンセレーションを避けるため）。すると
$$ \bm{H}_1 \bm{A} = \begin{pmatrix} \alpha & * & \cdots & * \\ 0 & & & \\ \vdots & & \bm{A}’ & \\ 0 & & & \end{pmatrix} $$
第1列が $(\alpha, 0, \ldots, 0)^T$ になりました。次に、残りの部分行列 $\bm{A}’$（$(m-1) \times (n-1)$ 行列）の第1列に対して同様のハウスホルダー変換 $\bm{H}_2’$ を適用します。これを $n$ ステップ繰り返すと
$$ \bm{H}_n \cdots \bm{H}_2 \bm{H}_1 \bm{A} = \bm{R} $$
各 $\bm{H}_k$ は直交行列なので、$\bm{Q} = \bm{H}_1 \bm{H}_2 \cdots \bm{H}_n$ も直交行列です。よって
$$ \bm{A} = \bm{Q}\bm{R} $$
ハウスホルダー法の優位性
ハウスホルダー法はグラム・シュミット法と比較して、数値安定性が格段に優れています。得られる $\bm{Q}$ の直交性の誤差は
$$ \|\bm{Q}^T\bm{Q} – \bm{I}\| \sim O(\varepsilon_{\text{mach}}) $$
と条件数に依存しません。これは、各ステップで適用するハウスホルダー行列が「完全な」直交行列であり、丸め誤差を除いて直交性が正確に保たれるためです。
計算量は $O(mn^2 – n^3/3)$ であり、MGSの $O(mn^2)$ と同程度ですが、メモリアクセスパターンがBLAS-3（行列-行列演算）に適しているため、実際の実行速度はMGSを上回ることが多いです。
理論面での考察が済んだので、次にこれらのアルゴリズムをPythonで実装し、数値的な振る舞いを確認しましょう。
Pythonでの実装
グラム・シュミット法の実装
まず、古典的グラム・シュミット法（CGS）と修正グラム・シュミット法（MGS）を実装します。数値安定性の違いを後で比較するため、両方を実装しておきます。
import numpy as np
import matplotlib.pyplot as plt
def qr_classical_gram_schmidt(A):
"""古典的グラム・シュミット法によるQR分解"""
m, n = A.shape
Q = np.zeros((m, n))
R = np.zeros((n, n))
for k in range(n):
# 全ての射影を元のベクトルに対して計算（CGS）
v = A[:, k].copy()
for j in range(k):
R[j, k] = Q[:, j] @ A[:, k]
v -= R[j, k] * Q[:, j]
R[k, k] = np.linalg.norm(v)
Q[:, k] = v / R[k, k]
return Q, R
def qr_modified_gram_schmidt(A):
"""修正グラム・シュミット法によるQR分解"""
m, n = A.shape
Q = np.zeros((m, n))
R = np.zeros((n, n))
V = A.copy().astype(float)
for k in range(n):
R[k, k] = np.linalg.norm(V[:, k])
Q[:, k] = V[:, k] / R[k, k]
# 最新のベクトルに対して射影を計算（MGS）
for j in range(k + 1, n):
R[k, j] = Q[:, k] @ V[:, j]
V[:, j] -= R[k, j] * Q[:, k]
return Q, R
# 具体例: 3x3行列のQR分解
A = np.array([[1, 1, 0],
[1, 0, 1],
[0, 1, 1]], dtype=float)
Q_mgs, R_mgs = qr_modified_gram_schmidt(A)
print("=== 修正グラム・シュミット法 ===")
print("Q =")
print(np.round(Q_mgs, 4))
print("\nR =")
print(np.round(R_mgs, 4))
print("\nQ^T Q (直交性の確認) =")
print(np.round(Q_mgs.T @ Q_mgs, 10))
print("\nQR - A (復元誤差) =")
print(np.round(Q_mgs @ R_mgs - A, 10))
上のコードでは、3行3列の行列に対してQR分解を実行し、$\bm{Q}^T\bm{Q}$ が単位行列になることと、$\bm{Q}\bm{R}$ が元の $\bm{A}$ を正確に復元することを確認しています。修正グラム・シュミット法では射影の計算順序が古典的方法と異なり、各ステップで更新済みのベクトルに対して内積を取ることで、丸め誤差の累積が抑えられています。
ハウスホルダー変換の実装
次にハウスホルダー変換によるQR分解を実装します。
import numpy as np
def qr_householder(A):
"""ハウスホルダー変換によるQR分解"""
m, n = A.shape
R = A.copy().astype(float)
Q = np.eye(m)
for k in range(min(m - 1, n)):
# 部分ベクトル
x = R[k:, k]
# ハウスホルダーベクトルの構成
alpha = -np.sign(x[0]) * np.linalg.norm(x)
u = x.copy()
u[0] -= alpha
v = u / np.linalg.norm(u)
# ハウスホルダー変換を適用
R[k:, k:] -= 2.0 * np.outer(v, v @ R[k:, k:])
Q[:, k:] -= 2.0 * np.outer(Q[:, k:] @ v, v)
return Q, R[:n, :] # 薄いQR分解
# ハウスホルダー法の結果
Q_hh, R_hh = qr_householder(A)
print("=== ハウスホルダー法 ===")
print("Q =")
print(np.round(Q_hh[:, :3], 4))
print("\nR =")
print(np.round(R_hh, 4))
print("\nQ^T Q (直交性の確認) =")
print(np.round(Q_hh[:, :3].T @ Q_hh[:, :3], 10))
ハウスホルダー法では、行列の各列に対してハウスホルダー反射行列を構成し、行列全体に左から適用していきます。出力される $\bm{Q}$ は累積的なハウスホルダー変換の積として得られます。
数値安定性の比較
3つの方法（CGS, MGS, ハウスホルダー）の数値安定性を条件数の異なる行列で比較しましょう。
import numpy as np
import matplotlib.pyplot as plt
def qr_classical_gram_schmidt(A):
m, n = A.shape
Q = np.zeros((m, n))
R = np.zeros((n, n))
for k in range(n):
v = A[:, k].copy()
for j in range(k):
R[j, k] = Q[:, j] @ A[:, k]
v -= R[j, k] * Q[:, j]
R[k, k] = np.linalg.norm(v)
Q[:, k] = v / R[k, k]
return Q, R
def qr_modified_gram_schmidt(A):
m, n = A.shape
Q = np.zeros((m, n))
R = np.zeros((n, n))
V = A.copy().astype(float)
for k in range(n):
R[k, k] = np.linalg.norm(V[:, k])
Q[:, k] = V[:, k] / R[k, k]
for j in range(k + 1, n):
R[k, j] = Q[:, k] @ V[:, j]
V[:, j] -= R[k, j] * Q[:, k]
return Q, R
def qr_householder(A):
m, n = A.shape
R = A.copy().astype(float)
Q = np.eye(m)
for k in range(min(m - 1, n)):
x = R[k:, k]
alpha = -np.sign(x[0]) * np.linalg.norm(x)
u = x.copy()
u[0] -= alpha
v = u / np.linalg.norm(u)
R[k:, k:] -= 2.0 * np.outer(v, v @ R[k:, k:])
Q[:, k:] -= 2.0 * np.outer(Q[:, k:] @ v, v)
return Q[:, :n], R[:n, :]
# 条件数を変えて直交性の誤差を比較
np.random.seed(42)
m, n = 50, 30
cond_numbers = np.logspace(1, 15, 15)
errors_cgs = []
errors_mgs = []
errors_hh = []
errors_np = []
for cond in cond_numbers:
# 指定した条件数の行列を生成
U, _ = np.linalg.qr(np.random.randn(m, m))
V, _ = np.linalg.qr(np.random.randn(n, n))
S = np.zeros((m, n))
singular_values = np.logspace(0, -np.log10(cond), n)
np.fill_diagonal(S, singular_values)
A = U @ S @ V.T
# 各手法でQR分解
Q_cgs, _ = qr_classical_gram_schmidt(A)
Q_mgs, _ = qr_modified_gram_schmidt(A)
Q_hh, _ = qr_householder(A)
Q_np, _ = np.linalg.qr(A)
# 直交性誤差 ||Q^T Q - I||
errors_cgs.append(np.linalg.norm(Q_cgs.T @ Q_cgs - np.eye(n)))
errors_mgs.append(np.linalg.norm(Q_mgs.T @ Q_mgs - np.eye(n)))
errors_hh.append(np.linalg.norm(Q_hh.T @ Q_hh - np.eye(n)))
errors_np.append(np.linalg.norm(Q_np.T @ Q_np - np.eye(n)))
fig, ax = plt.subplots(figsize=(10, 6))
ax.loglog(cond_numbers, errors_cgs, "ro-", linewidth=2, markersize=6, label="Classical GS")
ax.loglog(cond_numbers, errors_mgs, "bs-", linewidth=2, markersize=6, label="Modified GS")
ax.loglog(cond_numbers, errors_hh, "g^-", linewidth=2, markersize=6, label="Householder")
ax.loglog(cond_numbers, errors_np, "kd-", linewidth=2, markersize=6, label="numpy.linalg.qr")
# 理論的な参照線
eps = np.finfo(float).eps
ax.loglog(cond_numbers, eps * cond_numbers**2, "r:", alpha=0.4, label=r"$O(\epsilon \kappa^2)$")
ax.loglog(cond_numbers, eps * cond_numbers, "b:", alpha=0.4, label=r"$O(\epsilon \kappa)$")
ax.axhline(eps, color="gray", linestyle="--", alpha=0.3, label=r"$\epsilon_{\rm mach}$")
ax.set_xlabel("Condition number $\\kappa(A)$", fontsize=13)
ax.set_ylabel("Orthogonality error $\\|Q^TQ - I\\|$", fontsize=13)
ax.set_title("Numerical Stability of QR Decomposition Methods", fontsize=14)
ax.legend(fontsize=10, loc="upper left")
ax.grid(True, alpha=0.3, which="both")
ax.set_ylim(1e-17, 1e5)
plt.tight_layout()
plt.savefig("qr_stability.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、3つのQR分解アルゴリズムの数値安定性の違いが明確に読み取れます。
-
古典的グラム・シュミット法（赤丸）: 直交性誤差が条件数の2乗に比例して増大しています。条件数が $10^8$ を超えると $\|\bm{Q}^T\bm{Q} – \bm{I}\|$ が $O(1)$ 程度になり、直交性がほぼ完全に失われています。赤い点線の $O(\varepsilon\kappa^2)$ のラインに沿っていることが確認できます。
-
修正グラム・シュミット法（青四角）: 条件数の1乗に比例する誤差成長を示しています。CGSと比較して格段に安定ですが、条件数が $10^{15}$ 付近では誤差が $O(1)$ に達します。青い点線の $O(\varepsilon\kappa)$ のラインに対応しています。
-
ハウスホルダー法（緑三角）とNumPy（黒ひし形）: 条件数によらず、直交性誤差が機械イプシロン $\varepsilon_{\text{mach}} \approx 10^{-16}$ 近傍に留まっています。NumPyの
linalg.qrもハウスホルダー法に基づいているため、ほぼ同じ挙動を示しています。
この結果は、実用的なアプリケーションではハウスホルダー法を使うべきであること、そしてグラム・シュミット法を使う場合は必ず修正版（MGS）を用いるべきことを示しています。
最小二乗法への応用
QR分解の最も重要な応用の一つが、最小二乗法です。過決定系 $\bm{A}\bm{x} \approx \bm{b}$（$m > n$）の最小二乗解を求める問題を考えましょう。
正規方程式 $\bm{A}^T\bm{A}\bm{x} = \bm{A}^T\bm{b}$ を直接解くと、条件数が2乗になる（$\kappa(\bm{A}^T\bm{A}) = \kappa(\bm{A})^2$）ため数値的に不安定です。QR分解を使えば、この問題を回避できます。
$\bm{A} = \bm{Q}\bm{R}$ と分解すると
$$ \|\bm{A}\bm{x} – \bm{b}\|^2 = \|\bm{Q}\bm{R}\bm{x} – \bm{b}\|^2 $$
$\bm{Q}$ は直交行列なのでノルムを保存します。つまり $\|\bm{Q}\bm{y}\| = \|\bm{y}\|$ が成り立ちます。両辺に $\bm{Q}^T$ を左からかけると
$$ \|\bm{R}\bm{x} – \bm{Q}^T\bm{b}\|^2 $$
を最小化すればよいことになります。$\bm{R}$ は上三角行列なので、後退代入で解けます。
$$ \bm{R}\bm{x} = \bm{Q}^T\bm{b} $$
import numpy as np
import matplotlib.pyplot as plt
# 最小二乗法のデモ: 多項式フィッティング
np.random.seed(42)
n_data = 50
x_data = np.linspace(0, 2 * np.pi, n_data)
y_true = np.sin(x_data)
y_data = y_true + 0.3 * np.random.randn(n_data)
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 多項式次数を変えたフィッティング
ax = axes[0]
ax.scatter(x_data, y_data, s=30, alpha=0.6, color="gray", label="Data")
ax.plot(x_data, y_true, "k--", linewidth=1.5, alpha=0.5, label="True: sin(x)")
x_plot = np.linspace(0, 2 * np.pi, 200)
colors = ["blue", "green", "red"]
degrees = [3, 5, 10]
for deg, color in zip(degrees, colors):
# ヴァンデルモンド行列（デザイン行列）
A = np.vander(x_data, deg + 1, increasing=True)
# QR分解による最小二乗法
Q, R = np.linalg.qr(A)
coeffs = np.linalg.solve(R, Q.T @ y_data)
y_fit = np.vander(x_plot, deg + 1, increasing=True) @ coeffs
residual = np.linalg.norm(y_data - A @ coeffs)
ax.plot(x_plot, y_fit, color=color, linewidth=2,
label=f"degree {deg} (residual={residual:.3f})")
ax.set_xlabel("x", fontsize=12)
ax.set_ylabel("y", fontsize=12)
ax.set_title("Polynomial Fitting via QR Decomposition", fontsize=13)
ax.legend(fontsize=9, loc="lower left")
ax.grid(True, alpha=0.3)
# (b) 正規方程式 vs QR分解の数値安定性比較
ax = axes[1]
degrees_range = range(1, 20)
errors_ne = [] # 正規方程式
errors_qr = [] # QR分解
for deg in degrees_range:
A = np.vander(x_data, deg + 1, increasing=True)
# 正規方程式: (A^T A) x = A^T b
try:
coeffs_ne = np.linalg.solve(A.T @ A, A.T @ y_data)
err_ne = np.linalg.norm(A @ coeffs_ne - y_data)
except np.linalg.LinAlgError:
err_ne = np.nan
# QR分解
Q, R = np.linalg.qr(A)
coeffs_qr = np.linalg.solve(R, Q.T @ y_data)
err_qr = np.linalg.norm(A @ coeffs_qr - y_data)
errors_ne.append(err_ne)
errors_qr.append(err_qr)
ax.semilogy(list(degrees_range), errors_ne, "ro-", linewidth=2,
markersize=6, label="Normal Equations")
ax.semilogy(list(degrees_range), errors_qr, "b^-", linewidth=2,
markersize=6, label="QR Decomposition")
ax.set_xlabel("Polynomial degree", fontsize=12)
ax.set_ylabel("Residual norm $\\|Ax - b\\|$", fontsize=12)
ax.set_title("Numerical Stability: Normal Eq. vs QR", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which="both")
plt.tight_layout()
plt.savefig("qr_least_squares.png", dpi=150, bbox_inches="tight")
plt.show()
このグラフから、QR分解を用いた最小二乗法の利点が読み取れます。
-
左図（多項式フィッティング）: 3次、5次、10次の多項式を sin(x) のノイズ付きデータにフィッティングしています。次数が上がるほど残差は小さくなりますが、10次では過学習の兆候が見え始めます。QR分解による求解は全ての次数で安定に動作しています。
-
右図（数値安定性の比較）: 多項式の次数を上げていくと、ヴァンデルモンド行列の条件数が急激に増大します。正規方程式（赤丸）は次数15程度から数値誤差が暴走し始めますが、QR分解（青三角）は次数19まで安定した結果を返しています。これは正規方程式が $\kappa(\bm{A})^2$ の条件数で解を求めるのに対し、QR分解は $\kappa(\bm{A})$ の条件数で解を求めるためです。
まとめ
本記事では、QR分解の理論とアルゴリズム、Pythonでの実装について解説しました。
- QR分解は行列を直交行列 $\bm{Q}$ と上三角行列 $\bm{R}$ の積に分解する手法であり、連立方程式の求解、最小二乗法、固有値計算など幅広く応用される
- グラム・シュミット法は射影の除去によって列ベクトルを逐次直交化するアルゴリズムであり、QR分解の直感的な理解に最適である
- 修正グラム・シュミット法（MGS）は射影の計算順序を変更することで、直交性誤差を $O(\varepsilon\kappa^2)$ から $O(\varepsilon\kappa)$ に改善する
- ハウスホルダー変換はベクトルの鏡面反射を利用し、条件数に依存しない $O(\varepsilon)$ の直交性誤差を実現する最も安定なアルゴリズムである
- 最小二乗法への応用では、QR分解は正規方程式よりも数値的に安定であり、高次多項式フィッティングなど条件数が大きくなる場面で特に重要である
次のステップとして、以下の記事も参考にしてください。
- コレスキー分解の理論と実装 — 正定値行列に特化した分解手法
- ムーア・ペンローズ擬逆行列の理論 — QR分解とSVDの関係
- 行列の微分公式を完全まとめ — 行列演算の微分理論