5つの行列分解（コレスキー分解、LU分解、QR分解、スペクトル分解、特異値分解）
機械学習などでもよく利用される5つの行列分解について簡単に説明する。 理論的な説明を行い、Pythonによる実装を行う。
以下の素晴らしい記事を参考にした。
ソースコード
github
- jupyter notebook形式のファイルはこちら
google colaboratory
- google colaboratory で実行する場合はこちら
実行環境
OSはmacOSである。LinuxやUnixのコマンドとはオプションが異なることに注意。
!sw_vers
ProductName: macOS
ProductVersion: 13.5.1
BuildVersion: 22G90
!python -V
Python 3.9.17
基本的なライブラリをインポートし watermark を利用してそのバージョンを確認する。 ついでに乱数のseedの設定を行う。
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import random
import scipy
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
seed = 123
random_state = 123
random.seed(seed)
np.random.seed(seed)
from watermark import watermark
print(watermark(python=True, watermark=True, iversions=True, globals_=globals()))
Python implementation: CPython
Python version : 3.9.17
IPython version : 8.17.2
scipy : 1.11.2
numpy : 1.25.2
matplotlib: 3.8.1
Watermark: 2.4.3
コレスキー分解
コレスキー分解は正定値対称行列を三角行列に分解する手法である。 正定値対称行列$ \mathbf{A} $に対し、コレスキー分解は$ \mathbf{A} = \mathbf{L} \mathbf{L}^\top $となる。 ここで$ \mathbf{L} $は下三角行列であり、対角要素は正の値をとる。
例えば$ 3 \times 3 $の正定値対称行列をコレスキー分解すると、下三角行列1つで表される。
計算量は$ O(n^3) $程度であり、ガウス消去に比べて効率が良い場合もある。
コレスキー分解が存在するためには行列が正定値かつ対称である必要がある。正定値対称行列の判定は固有値がすべて正かどうかでわかる。 コレスキー分解を用いることで、回帰問題やカーネル行列の操作が高速化する場合がある。
コレスキー分解により、逆行列を明示的に計算せずに方程式を解けるメリットがある。 数値計算上、行列要素が大きく変動する場合は計算の安定性に注意が必要となる。
以下にPythonでの簡単な例を示す。
import numpy as np
# 正定値対称行列を用意
A = np.array([[4.0, 2.0], [2.0, 3.0]])
# NumPyにはcholesky関数がある
L = np.linalg.cholesky(A)
# 結果を表示
print("A =\n", A.round(2))
print("L =\n", L.round(2))
print("L * L^T =\n", (L @ L.T).round(2))
A =
[[4. 2.]
[2. 3.]]
L =
[[2. 0. ]
[1. 1.41]]
L * L^T =
[[4. 2.]
[2. 3.]]
上記の例では、$ \mathbf{A} = \begin{pmatrix} 4 & 2 \\ 2 & 3 \end{pmatrix} $である。
そのコレスキー分解$ \mathbf{L} $は下三角行列になる。
ここで$ \mathbf{L} \mathbf{L}^\top $を計算すると元の行列に戻ることが確認できる。
LU分解
LU分解は、任意の正方行列を上下三角行列に分解する代表的な手法である。 FEMなどの数値計算においてもよく利用される。
行列$ \mathbf{A} $を、下三角行列$ \mathbf{L} $と上三角行列$ \mathbf{U} $に分解する。したがって$ \mathbf{A} = \mathbf{L} \mathbf{U} $という形で表せる。
ただしピボット選択のため、実際の計算では置換行列$ \mathbf{P} $も介在する。 よって$ \mathbf{P} \mathbf{A} = \mathbf{L} \mathbf{U} $となることが多い。 行列のサイズが$n$のとき、LU分解の計算量は$ O(n^3) $程度とされる。 ガウスの消去法を体系的に再利用した分解方式とも言える。
線形連立方程式$ \mathbf{A}\mathbf{x} = \mathbf{b} $の解法で用いられるケースが多い。
実際、LU分解を行うと前進代入と後退代入で解を求められる。 行列を何度も使う場合、分解を一度するだけで解法の速度が向上する。
以下にPythonによる簡単な例を示す。
import numpy as np
import scipy.linalg as la
# 例示行列
A = np.array([[2.0, 1.0, 1.0], [4.0, 3.0, 3.0], [8.0, 7.0, 9.0]])
# LU分解 (部分ピボットあり)
P, L, U = la.lu(A)
print("A =\n", A)
print("P =\n", P)
print("L =\n", L.round(2))
print("U =\n", U.round(2))
print("P*A =\n", (P @ A).round(2))
print("L*U =\n", (L @ U).round(2))
A =
[[2. 1. 1.]
[4. 3. 3.]
[8. 7. 9.]]
P =
[[0. 1. 0.]
[0. 0. 1.]
[1. 0. 0.]]
L =
[[1. 0. 0. ]
[0.25 1. 0. ]
[0.5 0.67 1. ]]
U =
[[ 8. 7. 9. ]
[ 0. -0.75 -1.25]
[ 0. 0. -0.67]]
P*A =
[[4. 3. 3.]
[8. 7. 9.]
[2. 1. 1.]]
L*U =
[[8. 7. 9.]
[2. 1. 1.]
[4. 3. 3.]]
行列が大きいほどLU分解の効率が重要となり、並列計算やブロック分解技術が研究されている。 大規模数値計算のような高速数値計算システムにもLU分解は組み込まれている。
QR分解
QR分解は、行列を直交行列(正規直交行列)と上三角行列の積に分解する手法である。
$m \times n$行列$ \mathbf{A} $に対し、$ \mathbf{A} = \mathbf{Q} \mathbf{R} $と分解できる。 ここで$ \mathbf{Q} $は列が直交ベクトルで構成された$m \times m$行列(もしくは$m \times n$)である。 $ \mathbf{R} $は$n \times n$の上三角行列となる(列数が$n$の場合)。
主に最小二乗法や回帰分析での基礎となる計算で用いられる。
$ \mathbf{A}\mathbf{x} \approx \mathbf{b} $を解くとき、QR分解によって安定した解を得られる。 勾配降下法などと組み合わせて、数値的に安定な解の導出に使われることも多い。 QR分解のアルゴリズムにはグラム・シュミット法やハウスホルダー変換などがある。 特にハウスホルダー変換を用いる方法は数値的に安定であるとされる。 計算量は$ O(mn^2) $程度になることが多い。
行列$ \mathbf{Q} $が直交行列という性質から、転置行列をかける操作が容易になる。
以下にPythonでのQR分解の例を示す。
import numpy as np
# m×n行列Aを用意
A = np.array([[1.0, 1.0], [1.0, 2.0], [1.0, 3.0]])
# QR分解
Q, R = np.linalg.qr(A)
print("A =\n", A.round(2))
print("Q =\n", Q.round(2))
print("R =\n", R.round(2))
print("Q*R =\n", (Q @ R).round(2))
A =
[[1. 1.]
[1. 2.]
[1. 3.]]
Q =
[[-0.58 0.71]
[-0.58 -0. ]
[-0.58 -0.71]]
R =
[[-1.73 -3.46]
[ 0. -1.41]]
Q*R =
[[1. 1.]
[1. 2.]
[1. 3.]]
ここでは$ m = 3, n = 2 $の行列を例示している。 結果の$ \mathbf{Q} $は$ 3 \times 3 $、$ \mathbf{R} $は$ 3 \times 2 $だが、実質的に上三角部分を持つ。 行列$ \mathbf{Q} $は直交行列なので、$ \mathbf{Q}^\top \mathbf{Q} = \mathbf{I} $が成立する。 QR分解を用いると、回帰分析の正規方程式を安定的に解ける。
スペクトル分解
スペクトル分解(固有値分解)は、正方行列を固有値と固有ベクトルによって対角化する操作である。 特に対称行列の場合は固有値が実数になり、直交行列による対角化が可能となる。
$n \times n$行列$ \mathbf{A} $の固有ベクトルを列ベクトルとして並べた行列を$ \mathbf{V} $とする。 固有値を対角成分にもつ行列を$ \mathbf{\Lambda} $とすると、$ \mathbf{A} \mathbf{V} = \mathbf{V} \mathbf{\Lambda} $となる。 もし$ \mathbf{V} $が正則であれば、$ \mathbf{A} = \mathbf{V} \mathbf{\Lambda} \mathbf{V}^{-1} $が成り立つ。 対称行列なら$ \mathbf{V} $は直交行列となり、$ \mathbf{A} = \mathbf{V} \mathbf{\Lambda} \mathbf{V}^\top $が成り立つ。
スペクトル分解は行列の特性を調べる際に重要であり、固有値の分布で行列の性質を把握できる。 例えば機械学習においては、共分散行列やグラフラプラシアンなどの固有値解析が多用される。
また、計算量が$ O(n^3) $と高いため、大規模行列では効率化手法や近似的アプローチが考えられる。 以下にPythonによる固有値分解の例を示す。
import numpy as np
# 固有値分解可能な正方行列
A = np.array([[4.0, 1.0], [1.0, 3.0]])
# 固有値と固有ベクトル
eigvals, eigvecs = np.linalg.eig(A)
print("A =\n", A)
print("Eigenvalues =\n", eigvals.round(2))
print("Eigenvectors =\n", eigvecs.round(2))
# 再構成
Lambda = np.diag(eigvals)
V = eigvecs
A_reconstructed = V @ Lambda @ np.linalg.inv(V)
print("Reconstructed A =\n", A_reconstructed.round(2))
A =
[[4. 1.]
[1. 3.]]
Eigenvalues =
[4.62 2.38]
Eigenvectors =
[[ 0.85 -0.53]
[ 0.53 0.85]]
Reconstructed A =
[[4. 1.]
[1. 3.]]
ここで$ \mathbf{A} = \begin{pmatrix} 4 & 1 \\ 1 & 3 \end{pmatrix} $である。 固有値として2つの実数を得られ、その固有ベクトルで対角化可能な例となっている。 高次元のスペクトル分解は行列のサイズが巨大になるため、反復法や近似手法が活用されることが多い。 グラフ固有ベクトルを用いたクラスタリングなど、多岐にわたる応用が存在する。
特異値分解(SVD)
特異値分解は、任意の$m \times n$行列$ \mathbf{A} $を3つの行列に分解する有用な手法である。
$ \mathbf{A} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^\top $という形で表す。
ここで$ \mathbf{U} $は$m \times m$の直交行列、$ \mathbf{\Sigma} $は特異値を対角に並べた$m \times n$の対角行列、$ \mathbf{V} $は$n \times n$の直交行列である。
$ \mathbf{\Sigma} $の対角要素である特異値は非負の実数となる。
特異値分解は行列のランクを明らかにし、近似行列を得るのに役立つ。 低ランク近似を行うことで、大規模行列の圧縮やノイズ除去が可能となる。
特に推薦システムの行列補完問題や潜在因子モデルでSVDが鍵となる。 さらに画像処理や自然言語処理でも、データの次元圧縮にSVDが用いられる。
特異値分解の計算量は$ O(\min(mn^2, m^2n)) $程度のオーダーで比較的重い。
SVDの大きな利点として、任意の行列に適用でき、数値的に安定な点が挙げられる。 一方で計算コストとメモリ使用量が増大するので、大規模な場合は近似SVDが利用される。
以下にPythonでの特異値分解例を示す。
import numpy as np
A = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]])
U, S, Vt = np.linalg.svd(A, full_matrices=True)
print("A =\n", A)
print("U =\n", U.round(2))
print("S =\n", S.round(2))
print("V^T =\n", Vt.round(2))
# 特異値を対角行列に拡張
Sigma = np.zeros_like(A, dtype=float)
np.fill_diagonal(Sigma, S)
A_reconstructed = U @ Sigma @ Vt
print("Reconstructed A =\n", A_reconstructed.round(2))
A =
[[1. 2. 3.]
[4. 5. 6.]
[7. 8. 9.]]
U =
[[-0.21 0.89 0.41]
[-0.52 0.25 -0.82]
[-0.83 -0.39 0.41]]
S =
[16.85 1.07 0. ]
V^T =
[[-0.48 -0.57 -0.67]
[-0.78 -0.08 0.63]
[-0.41 0.82 -0.41]]
Reconstructed A =
[[1. 2. 3.]
[4. 5. 6.]
[7. 8. 9.]]
ここでは3×3の小さな行列でSVDを行った。
$ \mathbf{U} $は3×3の直交行列、$ \mathbf{\Sigma} $は3×3で対角要素が特異値となる。 $ \mathbf{V}^\top $は3×3の直交行列で、転置をとった形で得られる。 分解された行列の積を計算し、ちゃんと元の行列に戻ることが確認できる。
結論
この記事では、行列分解の代表的手法であるコレスキー分解、LU分解、QR分解、スペクトル分解、特異値分解を概説した。 行列分解は、機械学習やAIの大規模データ解析において計算効率を向上させる中心的な技術である。
行列が正定値対称であればコレスキー分解が有効であり、LU分解は一般的な正方行列の解法に用いられる。 QR分解は回帰や最小二乗問題に適し、スペクトル分解は固有値解析や対称行列に対して有用である。 最後に特異値分解は任意の行列に適用可能であり、低ランク近似による次元削減が可能となる。
個人的な備忘録を兼ねて、行列分解の基本的な手法をまとめた。