最小二乗法におけるLU、QR、SVDの安定性比較
概要
行列 $A \in \mathbb{R}^{m \times n}$について、$m \ge n$であるとしよう。LU分解とQR分解、そして特異値分解で最小二乗法を行うとき、三つの方法はいずれも後方安定的である。ただし、後方安定性を満たす条件が異なる。行列がフルランクであるとは、$A$のランク$r$が$\min (m, n)$と等しいという意味である。
説明
最小二乗法の実質的な実装として行列分解を使う方法はいくつもあるが、一般にはLU分解(コレスキー分解)による方法とQR分解による方法、そしてSVD分解による方法の三つが特に広く使われる。参考までに、LU分解であれコレスキー分解であれガウス消去法であれ、本質的には前進代入と後退代入を通じて線形システムを解く方法であるため、最小二乗法の実装という側面ではいずれも同じものと見なして差し支えない。実際、コレスキー分解のflopsは$mn^{2} + {\frac{ 1 }{ 3 }} n^{3}$でLU分解より少し良いが、現実的には安定性のためにピボッティングpivotingを含むPLU分解を使うことになる。
Julia4やMATLAB5の実装を見ると、事実上$A$において$m \ne n$のときはQR分解、$m = n$のときはLU分解を使うと見てよい。
LU分解
行列の条件数 $\kappa = \kappa(A)$を考えてみると、アルゴリズムの相対誤差の上限は$\kappa$にかかっているが、LU分解であれコレスキー分解であれ、標準方程式$A^{\top} A \mathbf{x} = A^{\top} \mathbf{b}$を立てる時点ですでに$\kappa \left( A^{\top} A \right) = \kappa^{2}$となり、ただでさえ良くない可能性のある条件数がさらに悪化する。
$A = LU$について$\rho = {\frac{ \max \left| U_{ij} \right| }{ \max \left| A_{ij} \right| }}$を$A$の成長因子growth factorというが、せいぜいピボッティングLUで$\rho = O (1)$であり、$A$の相対誤差が$\rho$とマシンイプシロン $\epsilon$について$O \left( \rho \epsilon \right)$であるとき後方安定性を持つ。
一目見ても感じられるだろうが、結局のところ後方安定的ではあるものの、くどくどと必要な条件が多くて事実上不安定的である。これはガウス消去法の基本行演算を行うとき、数を掛けたり足したりしながら生じる誤差がだんだん大きくなる可能性があるためである。
QR分解
- Step 2-4. 次を計算する。 $$ \mathbf{q}_{j} = {{ \mathbf{v}_{j} } \over {r_{jj} }} $$
QR分解も実装する方式がいくつもあるため、必ずしもこのような理由で安定性を失うとは言えないが、少なくともこのようにグラム・シュミット直交化を使う場合には、ベクトルを$r_{jj}$で割るなどの過程で問題が発生し得る。まさに$A$が悪条件であるか否か、言い換えればこのような特異性singularityが問題の中に潜んでいるかどうかにかかっていると言える。
しかしLUに比べるほどではなく、最も汎用的に使われるだけあって安定性を過度に心配する必要はない。参考までに、QR分解でもだいたいの場合ピボッティングが基本的に入る。
SVD
ランク欠損rank-deficiencyがある状況でも安定的なアルゴリズムはSVDだけである。数式的に$A = U \Sigma V^{\top}$において$\Sigma$の逆行列 $\Sigma^{-1}$を求めるときには、$A$の特異値のうち$0$に近い部分を丸ごと除去した縮小特異値分解reduced SVDを試みることができる。これはピボッティングに依存しなければならない他の分解法と違い、十分に小さい閾値を明示的に設定して数値的な問題を根本から封じ込めるのと同じである。基本的なコストが高いため最小二乗法を実装するためによく使われるわけではないが、その分安定性においては決して劣らないのである。
速度とのトレードオフ
三つの方法の速度を比較してみて安定性を比較してみると、次のように非常に明確にトレードオフを確認できる。
| 方法 | 速度 | 安定性 |
|---|---|---|
| LU | 速い | 不安定 |
| QR | 普通 | 安定的 |
| SVD | 遅い | 非常に安定的 |
Trefethen, L. N., & Bau, D. (2022). Numerical linear algebra. Society for Industrial and Applied Mathematics. p142, 153 ↩︎
Trefethen, L. N., & Bau, D. (2022). Numerical linear algebra. Society for Industrial and Applied Mathematics. p141 ↩︎
Trefethen, L. N., & Bau, D. (2022). Numerical linear algebra. Society for Industrial and Applied Mathematics. p143 ↩︎
https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#Base.:\-Tuple{AbstractMatrix,%20AbstractVecOrMat} ↩︎
https://kr.mathworks.com/help/matlab/ref/double.mldivide.html ↩︎
