QR 分解 — Householder・Givens・MGS・列ピボット

概要

QR 分解 $A = QR$ は、$m \times n$ 行列 $A$（$m \geq n$）を直交行列 $Q$（$m \times m$、$Q^\top Q = I$）と上三角行列 $R$（$m \times n$）の積に分解する。 直交変換は条件数を変えず、丸め誤差の伝播も最小限であるため、QR は数値計算の中で最も安定な行列分解の一つである。

QR を構築する古典的手段は次の 3 つで、それぞれ得手不得手がある:

• Householder 反射 — 密行列を一括三角化、フロップ数最小
• Givens 回転 — 単一要素をゼロ化、疎・帯・Hessenberg 構造で有利
• Gram-Schmidt 直交化 — 列ベクトル単位で直交基底を逐次構成、Arnoldi 等の Krylov 法で必須

本記事ではこれらの数学的構造、数値安定性、ブロック化による BLAS-3 最適化、 列ピボットによる rank revealing 拡張、最小二乗や固有値計算への応用を順に述べる。

API の使い方は LinAlg API — QR 分解 を参照。

Householder 反射

Householder 反射は単位ベクトル $\mathbf{v}$ に対して

$$H = I - 2 \frac{\mathbf{v} \mathbf{v}^\top}{\mathbf{v}^\top \mathbf{v}}$$

と定義される直交対称行列で、$\mathbf{v}$ に直交する超平面に関する鏡映を表す。 任意のベクトル $\mathbf{x}$ に対し、$\mathbf{v} = \mathbf{x} \pm \|\mathbf{x}\|_2 \mathbf{e}_1$ と選べば $H\mathbf{x} = \mp \|\mathbf{x}\|_2 \mathbf{e}_1$ となり、$\mathbf{x}$ を第 1 軸に一気に揃えられる。 符号は $\mathrm{sign}(x_1)$ を採ることで $\mathbf{v}$ の桁落ちを回避する（Wilkinson の符号則）。

三角化アルゴリズムは次のとおり:

1. 第 1 列 $\mathbf{a}_1$ に対して反射子 $H_1$ を構成し、$H_1 A$ で第 1 列を $\mathbf{e}_1$ に揃える
2. 残りの $(m-1) \times (n-1)$ 部分行列に対して再帰的に同じ操作を $H_2, H_3, \ldots, H_n$ と繰り返す
3. 合成 $H_n \cdots H_2 H_1 A = R$ が上三角、よって $Q = H_1^\top H_2^\top \cdots H_n^\top = H_1 H_2 \cdots H_n$（$H_i$ は対称）

フロップ数は $2mn^2 - \tfrac{2}{3}n^3$ で、密行列の QR 分解として漸近的に最適。 $Q$ は明示的に組み立てずに、反射ベクトル $\mathbf{v}_1, \ldots, \mathbf{v}_n$ を $R$ の下三角部分に packed format で保持し、 必要に応じて $Q\mathbf{x}$ や $Q^\top \mathbf{b}$ を行列-ベクトル積で評価する暗黙的 Q表現が標準である。 これにより、$Q$ を陽に格納すれば $O(m^2)$ 必要なメモリが $O(mn)$ で済む。

ブロック Householder と WY 表現

素朴な Householder QR は $H_k$ をひとつずつ右から行列に適用するため、各適用が Level-2 BLAS（行列-ベクトル積）となる。 WY 表現は連続する $b$ 個の反射子をまとめて単一の行列乗算 (Level-3 BLAS, gemm) に書き換える。

$k$ 番目までの反射子の積を

$$H_1 H_2 \cdots H_k = I - W_k T_k W_k^\top$$

と表す。ここで $W_k \in \mathbb{R}^{m \times k}$ は反射ベクトルを並べた行列、$T_k \in \mathbb{R}^{k \times k}$ は上三角行列で、 漸化式

$$T_k = \begin{pmatrix} T_{k-1} & -\tau_k T_{k-1} W_{k-1}^\top \mathbf{v}_k \\ 0 & \tau_k \end{pmatrix}$$

で更新する（$\tau_k = 2 / (\mathbf{v}_k^\top \mathbf{v}_k)$）。 ブロックサイズ $b$ ごとにパネル分解を行い、残りのトレーリング行列への適用を $(I - WTW^\top) A_{[:, b{:}n]}$ という形で一回の $\mathrm{gemm}$ 二回で済ます。フロップ数は同じだが、メモリアクセス比が $O(1/b)$ に下がり、 L2 キャッシュに収まるブロックでは 5〜10 倍の高速化が得られる。

sangi は MKL/OpenBLAS が利用可能なときは dgeqrf 経由でブロック QR を呼び、 ない場合は内部で同等の WY 表現実装にフォールバックする。

Givens 回転と fast Givens

Givens 回転は 2 つの行（または列）に対して局所的にユニタリ変換を施す:

$$G(i, k, \theta) = \begin{pmatrix} \cdots & & & \\ & c & \cdots & s & \\ & \vdots & \ddots & \vdots & \\ & -s & \cdots & c & \\ & & & \cdots \end{pmatrix}, \quad c = \cos\theta, \, s = \sin\theta$$

$c$ と $s$ は $\mathbf{x} = (x_i, x_k)^\top$ に対して $c = x_i / r$, $s = x_k / r$, $r = \sqrt{x_i^2 + x_k^2}$ と取り、 $G^\top \mathbf{x} = (r, 0)^\top$ で第 $k$ 成分を消去する。1 回の Givens 回転は 1 要素のみをゼロ化するので、 $m \times n$ の密行列の QR には $m n / 2$ 個の Givens を要し、Householder の約 2 倍のフロップを必要とする。

fast Givens（sqrt 回避）

$\sqrt{x_i^2 + x_k^2}$ は 1 回の平方根を要し、組込みでも数命令を消費する。 fast Givens 変換は対角スケーリング行列を別途保持して、回転行列を有理数行列 $M = \mathrm{diag}(d_1, d_2)^{-1} \tilde{M}$ の形で表現し、 平方根を回避する。代償としてスケール行列 $D$ の追跡が必要で、最終的に $A \cdot D^{1/2}$ で復元する。

Givens が有利な場面

密行列では Householder に劣るが、次の構造を持つ行列では Givens が圧勝する:

• Hessenberg 行列: 下副対角のみ消せばよく、$n-1$ 個の Givens で完了。QR アルゴリズムの内ループはこの形。
• 帯行列: バンド外の零を保ったままバンド内のみで Givens を適用できる。Householder では fill-in が出る。
• 疎行列: 既存の非零パターンに沿った fill-in 最小化が局所操作で可能。
• 更新型 QR: 既存の QR に 1 行を追加・削除する Sherman-Morrison 風の更新が Givens の連続で実行できる。

sangi は固有値計算（QR アルゴリズム、Hessenberg 縮約後）と疎 QR で Givens を、密 QR で Householder を採用する。

古典 vs 修正 Gram-Schmidt

Gram-Schmidt 直交化は、列ベクトルを順に直交化して QR を構成する。 ベクトル単位で並列性が低い反面、列を一本ずつ生成できるため Krylov 部分空間法（Arnoldi・GMRES）では必須となる。

古典 Gram-Schmidt (CGS)

各 $\mathbf{a}_k$ に対し、先行の $\mathbf{q}_1, \ldots, \mathbf{q}_{k-1}$ への射影を一括計算する:

$$r_{jk} = \mathbf{q}_j^\top \mathbf{a}_k, \quad \tilde{\mathbf{q}}_k = \mathbf{a}_k - \sum_{j=1}^{k-1} r_{jk} \mathbf{q}_j, \quad r_{kk} = \|\tilde{\mathbf{q}}_k\|_2, \quad \mathbf{q}_k = \tilde{\mathbf{q}}_k / r_{kk}$$

$r_{jk}$ の計算が並列化しやすい（独立な内積）反面、$\mathbf{q}_j$ 自体に丸め誤差があると射影量の見積もりが狂い、 計算済みの直交性 $\mathbf{q}_i^\top \mathbf{q}_j$ が $O(u \kappa^2(A))$ で崩れる。 条件数が $1/\sqrt{u}$ を超える行列では直交性が完全に失われる。

修正 Gram-Schmidt (MGS)

$\mathbf{a}_k$ から各 $\mathbf{q}_j$ への射影を逐次除去し、その都度更新したベクトルを次の射影に使う:

v = a_k
for j = 1 ... k-1:
    r_jk = q_j^T v
    v = v - r_jk * q_j
r_kk = ||v||, q_k = v / r_kk

数学的には CGS と等価だが、各除去後に「すでに直交化済み」の状態でベクトルを使うため、 直交性誤差が $O(u \kappa(A))$ に抑えられる。並列性は低下するが、安定性が条件数の 1 乗で抑えられる利点が大きい。

sangi の Arnoldi 反復・GMRES は MGS を採用する。 さらに数値的に厳しい場合は再直交化（一度 MGS した後、もう一度 MGS を回す）を行い、誤差を $O(u)$ に下げる。

経済形 QR と完全 QR

$m \geq n$ の場合、$A = QR$ には次の 2 形式がある:

• 完全 QR: $Q \in \mathbb{R}^{m \times m}$（直交行列）、$R \in \mathbb{R}^{m \times n}$（上三角、下 $m-n$ 行は零）
• 経済形 (thin) QR: $Q_1 \in \mathbb{R}^{m \times n}$（列直交）、$R_1 \in \mathbb{R}^{n \times n}$（上三角）。$Q_1$ は完全 $Q$ の最初の $n$ 列に対応。

多くの応用では $Q_1$ さえあれば足りる:

• 最小二乗 $\min \|A\mathbf{x} - \mathbf{b}\|_2$ は $R_1 \mathbf{x} = Q_1^\top \mathbf{b}$ の三角ソルブで解ける
• $A$ の列空間への射影行列は $Q_1 Q_1^\top$（$Q Q^\top = I$ ではなく）
• 低ランク近似や Krylov 基底も $Q_1$ のみで構成

$m \gg n$ の場面（例: 1 万行 × 100 列の最小二乗）では完全 $Q$ を保持すると $O(m^2)$ のメモリが必要だが、経済形では $O(mn)$ で済む。 sangi の qr_decomposition は完全 QR を返すが、最小二乗ソルバの内部では暗黙的 $Q$ 表現を用い、 $Q$ を陽に組み立てずに $Q^\top \mathbf{b}$ のみを計算してメモリを節約している。

列ピボット QR（rank revealing）

ランクが不完全あるいは数値的に低い行列に対しては、列を並べ替えながら QR を進めると、 $|R_{kk}|$ の階段状の落ち込みからランクを読み取れる。 最も古典的な戦略が Businger-Golub 列ピボット QR である。

Businger-Golub の手続き

第 $k$ 段において、残り部分行列 $A_{[k:m, k:n]}$ の各列の 2 ノルム

$$c_j^{(k)} = \|A_{[k:m, j]}^{(k)}\|_2, \quad j = k, \ldots, n$$

を計算し、最大の $c_p^{(k)}$ を持つ列 $p$ を $k$ 列目と交換してから Householder 反射を適用する。 交換は置換行列 $\Pi$ で記録し、最終的に

$$A \Pi = Q R, \quad |R_{11}| \ge |R_{22}| \ge \cdots \ge |R_{nn}|$$

を得る。数値ランクは $|R_{kk}| < \tau |R_{11}|$ となる最初の $k$ で読み取れる（$\tau$ は機械精度や問題依存のしきい値）。

列ノルムの増分更新

素朴に毎段で全列のノルムを再計算すると $O(mn)$/段になり、QR 全体が $O(mn^2)$ から $O(mn^2)$ に変わらないが定数が 2 倍以上に増える。 Businger-Golub では既存ノルムから差分更新

$$c_j^{(k+1)} = \sqrt{(c_j^{(k)})^2 - r_{kj}^2}$$

を行い、桁落ちが疑われる場合のみ再計算するハイブリッド戦略でコストを抑える。

応用

列ピボット QR は SVD 並みの数値的ランク推定を $O(mn^2)$ で提供し、最小二乗の rank-deficient 解、 部分集合選択（変数選択）、低ランク近似、TLS 近似に広く使われる。 ランダム化 SVD の前段で出現することも多い（後述の SVD 記事参照）。

応用

最小二乗法

$m > n$ の過決定系 $\min \|A\mathbf{x} - \mathbf{b}\|_2$ は、QR 分解 $A = QR$ により

$$\|A\mathbf{x} - \mathbf{b}\|_2^2 = \|R\mathbf{x} - Q^\top \mathbf{b}\|_2^2 = \|R_1 \mathbf{x} - (Q^\top \mathbf{b})_{[1:n]}\|_2^2 + \|(Q^\top \mathbf{b})_{[n+1:m]}\|_2^2$$

と分解できる。第 2 項は $\mathbf{x}$ に依存せず、最小化は $R_1 \mathbf{x} = (Q^\top \mathbf{b})_{[1:n]}$ という三角系を解くだけで済む。 正規方程式 $A^\top A \mathbf{x} = A^\top \mathbf{b}$ より条件数が良く（$\kappa(A^\top A) = \kappa(A)^2$ に対し QR は $\kappa(A)$）、sangi は最小二乗の既定実装に QR を採用する。

Arnoldi 反復

非対称行列に対する Krylov 部分空間 $\mathcal{K}_k(A, \mathbf{v}) = \mathrm{span}(\mathbf{v}, A\mathbf{v}, \ldots, A^{k-1}\mathbf{v})$ の正規直交基底を、 各ステップで $A \mathbf{q}_j$ から MGS で射影成分を取り除いて構築するのが Arnoldi 反復である。 副産物として得られる上 Hessenberg 行列 $H_k$ が、$A$ の Krylov 部分空間への射影を表し、GMRES や Krylov 固有値ソルバの基幹となる。

QR アルゴリズム（固有値計算）

$A_0 = A$ とおいて

$$A_k = Q_k R_k, \quad A_{k+1} = R_k Q_k$$

を反復すると、$A_k$ が（一般には複素 Schur 形の実版である）擬似 Schur 形に収束する。 実用版では Hessenberg 縮約 + Francis のダブルシフトでシフトを与え、Givens 回転で 1 段あたり $O(n^2)$ に落とす。 詳細は 固有値分解の記事を参照。

sangi の選択

参考文献

• Golub, G. H. and Van Loan, C. F. (2013). Matrix Computations, 4th ed., Johns Hopkins University Press. Chapter 5.
• Householder, A. S. (1958). "Unitary triangularization of a nonsymmetric matrix". J. ACM, 5(4), 339–342.
• Businger, P. and Golub, G. H. (1965). "Linear least squares solutions by Householder transformations". Numer. Math., 7, 269–276.
• Schreiber, R. and Van Loan, C. F. (1989). "A storage-efficient WY representation for products of Householder transformations". SIAM J. Sci. Stat. Comput., 10(1), 53–57.
• Björck, Å. (1996). Numerical Methods for Least Squares Problems, SIAM.