SVD — Golub-Reinsch・Jacobi・分割統治・ランダム化

概要

特異値分解 (SVD) は $m \times n$ 行列 $A$ に対して

$$A = U \Sigma V^\top, \quad U \in \mathbb{R}^{m \times m}, \, V \in \mathbb{R}^{n \times n} \text{（直交）}, \, \Sigma = \mathrm{diag}(\sigma_1 \geq \sigma_2 \geq \cdots \geq 0)$$

として、ランク・条件数・低ランク近似・擬似逆行列など線形代数のほぼ全ての中心的量を統一的に抽出する。 固有値分解と異なり、対称性も正方性も要求せず、任意の長方行列に対して常に存在する。

SVD アルゴリズムは大きく分けて

• Golub-Reinsch — 標準的な 密 SVD、二段階法
• Jacobi SVD — 高相対精度・並列向き、反復回数多
• 分割統治 (D&C) — 大きな $n$ で漸近的に速い
• ランダム化 SVD — 低ランク近似に特化、$k \ll \min(m,n)$ で圧倒的に速い

本記事ではそれぞれの内部構造、Demmel-Kahan zero-shift QR による高精度化、 Eckart-Young 定理を背景とした低ランク近似、PCA / 擬似逆行列 / 最小二乗への応用を解説する。

API の使い方は LinAlg API — SVD および ランダム化 SVD を参照。

Golub-Reinsch 法

1970 年に Golub と Reinsch が提案した標準的な 密 SVD アルゴリズムで、二段構成である。

第 1 段: Householder による上二重対角化

両側から Householder 反射を交互に適用して、$A$ を上二重対角行列 $B$ に変換する:

$$U_1^\top A V_1 = B = \begin{pmatrix} \alpha_1 & \beta_1 & & \\ & \alpha_2 & \beta_2 & \\ & & \ddots & \beta_{n-1} \\ & & & \alpha_n \end{pmatrix}$$

$U_1$ は左から $n$ 個、$V_1$ は右から $n-2$ 個の Householder の積。 フロップ数は $4mn^2 - \tfrac{4}{3}n^3$。 この段階で問題が二重対角に縮約されるため、以降の反復は $n$ 次元の問題になる。

第 2 段: 陰的シフト QR 反復で対角化

$B$ の特異値は $B^\top B$ の固有値の平方根に等しい。 $B^\top B$ は対称三重対角行列なので、対称 QR 反復で対角化できる。 しかし $B^\top B$ を陽に組み立てると条件数が $\kappa(B)^2$ に悪化するため、 実用アルゴリズムは $B$ に対して直接 Givens 回転を適用し、$B^\top B$ の QR 反復と等価な変換を陰的に実行する。 Wilkinson シフトを用いると 1 ステップあたり O(n) のコストで、平均 $O(n)$ 反復で全特異値が確定する。

合計コストは $O(mn^2 + n^3)$ で、密 SVD のベンチマーク。sangi の svd_decomposition はこの方式を採用する。

Demmel-Kahan zero-shift QR

素朴な Wilkinson シフト QR は収束が速いが、二重対角の特異値の相対精度は

$$|\hat{\sigma}_i - \sigma_i| / \sigma_i = O(u \cdot \sigma_1 / \sigma_i)$$

と悪化し、最小特異値 $\sigma_n$ では条件数倍の劣化を受ける。 Demmel と Kahan (1990) はシフト 0 で QR 反復を回すと、 二重対角の数値計算で乗除のみが現れ加減算による桁落ちが起こらないため、 全特異値の相対精度が $O(nu)$ に抑えられることを示した。

実際の実装はハイブリッド戦略で、ほとんどのステップは Wilkinson シフトで速く収束させ、 絶対値が小さくなった最後の数ステップだけ zero-shift に切り替える。 LAPACK の xBDSQR ルーチンがこの戦略を採用、sangi の SVD 内部でも同等の処理を行う。

この高相対精度は、特に小さな特異値が物理的に意味を持つ問題（信号処理の SVD ノイズ除去、 ナノスケール解析のスケール分離、悪条件最小二乗の正則化パラメータ選択）で重要である。

Jacobi SVD

Jacobi 法は対称行列の固有値計算（後述の固有値ページ参照）に古典的に使われるが、 SVD にも自然に拡張できる。Forsythe-Henrici-Hestenes による one-sided Jacobi SVD はとくに高精度・並列向き。

one-sided Jacobi の手続き

$A$ の列ベクトル $\mathbf{a}_i, \mathbf{a}_j$ のペアに対して、$2 \times 2$ の平面回転 $J_{ij}$ を右から作用させて 新しい 2 列が直交するように選ぶ:

$$A J_{ij} = A', \quad (\mathbf{a}'_i)^\top \mathbf{a}'_j = 0$$

すべての列ペアに対して一巡したものを 1 スイープと呼び、全ペアが既に直交になるまで繰り返す。 収束後、各列のノルム $\|\mathbf{a}'_i\|_2$ が特異値 $\sigma_i$、正規化した列が $U$、累積した回転積が $V^\top$ を与える。

高相対精度

各回転は 2 列のなす角を直接計算するため、行列要素の相対誤差が増幅されない。 Demmel-Veselić (1992) の解析により、全特異値が相対精度 $O(u \cdot \kappa(A_s))$ で得られる （$A_s$ は $A$ の列スケーリング後の行列で、$\kappa(A_s) \leq \kappa(A)$）。 ill-graded な問題（特異値が桁違いに散らばっている）では Golub-Reinsch を凌駕する。

並列性

列ペアを非干渉に分割すれば、1 スイープあたり $\lceil n/2 \rceil$ 個の独立な $2 \times 2$ 問題に分解でき、並列性が高い。 Round-robin スケジューリングで全ペアを巡回する。 反復回数は 5-20 スイープ、$O(n)$ オーダーで密の Golub-Reinsch より定数倍大きい。

分割統治 SVD

大きな二重対角行列 $B$ に対しては、分割統治 (divide-and-conquer; D&C) SVD が漸近的に速い。 LAPACK の xGESDD ルーチンがこの方式を採用する。

$B$ を真ん中で 2 つの部分二重対角行列 $B_1, B_2$ と接続項に分割し、 $B_1, B_2$ の SVD を再帰的に求めて結合する。結合段は secular equation

$$1 + \rho \sum_{i=1}^{n} \frac{\zeta_i^2}{d_i - \lambda} = 0$$

を解く問題に帰着し、各根が結合行列の固有値を与える。secular 根は Gragg-Reinsch の手法で $O(1)$ 反復で確定する。

計算量は典型的に $O(n^2 \log n)$ で、Golub-Reinsch の $O(n^3)$ より良い。 実用上、$n$ が 1000 を超える領域で D&C が Golub-Reinsch を抜き、$n = 10^4$ では 5-10 倍速。 sangi の SVD 関数は MKL/OpenBLAS を経由する経路で LAPACKE_dgesvd (Golub-Reinsch) にディスパッチする。 大行列での LAPACKE_dgesdd (D&C) への切替は今後の課題である。

ランダム化 SVD

$A$ の主成分が上位 $k$ 個の特異値・特異ベクトルだけでよく、$k \ll \min(m, n)$ である場合、 ランダム化 SVD は $O(mn k)$ で近似 SVD を計算する。 Halko-Martinsson-Tropp (2011) が誤差解析を完成させた現代的アルゴリズム。

3 段階アルゴリズム

1. 値域捕捉: ランダムガウス行列 $\Omega \in \mathbb{R}^{n \times (k+p)}$（$p$ は oversampling、典型 5-10）を取り、 $Y = A \Omega$ を計算。$Y$ の列は $A$ の値域空間を高確率で張る。
2. 正規直交化: $Y$ の QR 分解で $Q$（$m \times (k+p)$ の直交列）を得る。$Q$ は $A$ の値域の近似基底。
3. 小型 SVD: $B = Q^\top A \in \mathbb{R}^{(k+p) \times n}$ の標準 SVD $B = \tilde{U} \Sigma V^\top$ を計算し、 $U = Q \tilde{U}$ を取って $A \approx U \Sigma V^\top$ を得る。

誤差バウンド

Halko et al. の解析により、期待誤差は次で抑えられる:

$$\mathbb{E} \|A - Q Q^\top A\|_2 \leq \left(1 + \frac{4 \sqrt{k+p}}{p-1} \sqrt{\min(m,n)}\right) \sigma_{k+1}(A)$$

$\sigma_{k+1}$ が小さい（スペクトル減衰が速い）ほど近似が精確になる。 パワー反復 $Y = (AA^\top)^q A \Omega$ を加えると残差項が $\sigma_{k+1}^{2q+1}$ で減衰し、低ランク近似の品質が劇的に改善する。

使いどころ

$10^5 \times 10^5$ 規模の行列で上位 50-500 個の特異値を取る、画像・推薦行列の低ランク近似、 大規模 PCA、データ圧縮、潜在意味解析（LSA）に広く使われる。 sangi の randomized_svd がこの実装を提供する。

経済形 SVD と打切り SVD

$m \geq n$ の場合、完全 SVD は $U \in \mathbb{R}^{m \times m}$ を持つが、 最後の $m-n$ 列は零特異値に対応し常に零空間の基底にすぎない。 実用上は最初の $n$ 列のみ残した経済形 (thin) SVD で十分:

$$A = U_n \Sigma_n V^\top, \quad U_n \in \mathbb{R}^{m \times n}, \, \Sigma_n \in \mathbb{R}^{n \times n}, \, V \in \mathbb{R}^{n \times n}$$

$m \gg n$ では完全 $U$ が $O(m^2)$ メモリを食うため、経済形が必須。

打切り SVD と Eckart-Young 定理

特異値を最初の $k$ 個だけ残した

$$A_k = \sum_{i=1}^{k} \sigma_i \mathbf{u}_i \mathbf{v}_i^\top$$

を打切り SVD と呼ぶ。Eckart-Young (1936) の定理は、Frobenius ノルムおよびスペクトル（2）ノルムにおいて $A_k$ が $A$ のランク $k$ 最良近似であることを述べる:

$$\min_{\mathrm{rank}(B) \leq k} \|A - B\|_F = \sqrt{\sigma_{k+1}^2 + \cdots + \sigma_r^2}, \quad \min_{\mathrm{rank}(B) \leq k} \|A - B\|_2 = \sigma_{k+1}$$

これが PCA、画像圧縮、推薦システム、潜在意味解析の理論的基盤である。

応用

最小二乗法（rank-deficient）

$A$ が rank-deficient あるいは悪条件の場合、SVD で擬似逆行列を構成する:

$$A^+ = V \Sigma^+ U^\top, \quad \Sigma^+_{ii} = \begin{cases} 1/\sigma_i & \sigma_i > \tau \\ 0 & \sigma_i \leq \tau \end{cases}$$

$\tau$ は数値ランクの閾値（典型的に $\tau = \mathrm{eps} \cdot \sigma_1 \cdot \max(m,n)$）。 $A\mathbf{x} = \mathbf{b}$ の最小ノルム最小二乗解は $\mathbf{x} = A^+ \mathbf{b}$。 sangi の svd_solve は rcond パラメータでこの閾値を制御する。

主成分分析（PCA）

中心化されたデータ行列 $X \in \mathbb{R}^{n \times d}$（$n$ サンプル、$d$ 次元）に対して、 共分散 $C = X^\top X / n$ の固有値分解が PCA に等しい。 SVD $X = U \Sigma V^\top$ を直接計算すれば、$V$ の列が主成分軸、$\sigma_i^2 / n$ が分散寄与となり、 $C$ を明示的に組まずに同じ結果が得られる。 $X$ が悪条件のとき、$C$ では条件数が二乗されて精度を失うため、SVD 経由が標準。

条件数とランク

$\kappa_2(A) = \sigma_1 / \sigma_n$（$\sigma_n > 0$）が 2 ノルム条件数。 SVD のみがこれを直接与える。 数値ランクは「絶対値が $\tau$ 以上の特異値の個数」として定義され、ノイズ混入データの実効ランクを安定に判定できる。

低ランク近似

画像圧縮、推薦行列補完、潜在意味解析、グラフ埋め込み、データ可視化など、 $A_k = \sum_{i=1}^{k} \sigma_i \mathbf{u}_i \mathbf{v}_i^\top$ の形で大規模データを縮約する手法は全て Eckart-Young 定理を背景とする。 大規模問題ではランダム化 SVD と組み合わせるのが定石。

sangi の選択

参考文献

• Golub, G. H. and Reinsch, C. (1970). "Singular value decomposition and least squares solutions". Numer. Math., 14, 403–420.
• Eckart, C. and Young, G. (1936). "The approximation of one matrix by another of lower rank". Psychometrika, 1, 211–218.
• Demmel, J. and Kahan, W. (1990). "Accurate singular values of bidiagonal matrices". SIAM J. Sci. Stat. Comput., 11(5), 873–912.
• Demmel, J. and Veselić, K. (1992). "Jacobi's method is more accurate than QR". SIAM J. Matrix Anal. Appl., 13(4), 1204–1245.
• Halko, N., Martinsson, P.-G., and Tropp, J. A. (2011). "Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions". SIAM Review, 53(2), 217–288.
• Trefethen, L. N. and Bau, D. (1997). Numerical Linear Algebra, SIAM. Lectures 4–5, 31–32.