特異値分解・LU分解・コレスキー分解

はじめに

正方行列でない行列に対しても固有値のような性質を利用したい場合があります。そのときは特異値分解を行い、特異値を使用します。
他には、行列の計算量を減らすために、行列を分解する場合があります。

参考

• 高校数学の美しい物語『特異値分解の定義，性質，具体例』，https://mathtrain.jp/svd
• 東條遼平，『LU分解』，http://hooktail.org/computer/index.php?LU%CA%AC%B2%F2
• SAM猫，『非負値行列因子分解（NMF）によるレコメンドのちょっとした例』，http://smrmkt.hatenablog.jp/entry/2014/08/23/211555

手順

特異値分解

$A$を$(m,n)$行列とします。$U$は$AA^\dagger$の固有ベクトルを列に並べた$(m,m)$ユニタリ行列です。$V$は$A^\dagger A$の固有ベクトルを列に並べた$(n,n)$ユニタリ行列です。$\Sigma$は特異値を大きさ順に対角に並べた$(m,n)$行列です。

A = U \Sigma V^\dagger

これらはsp.linalg.svdで求めることができます。
戻り値のUは左特異ベクトルを列にもつユニタリ行列、sは特異値、Vhは右特異ベクトルを行に持つユニタリ行列です。なのでVhは上の$V^\dagger$に相当します。

from scipy import linalg

A = np.array([[3,6,1,8],
              [6,2,1,1],
              [7,8,2,2]])

U, s, Vh = sp.linalg.svd(A)
m, n = A.shape

S = sp.linalg.diagsvd(s,m,n)

print(U@S@Vh)

# array([[3., 6., 1., 8.],
#        [6., 2., 1., 1.],
#        [7., 8., 2., 2.]])

LU分解

$A$を正方行列として、連立方程式$y = Ax$を解くことを考えます。$A$を下三角行列$L$と上三角行列$U$に$A = LU$と分解することで、連立方程式を解きやすくします。

数値計算では計算の安定性のために、行列の要素を並び替える置換行列$P$を用いて$A=PLU$とします。

A = np.array([[6, 4, 1],
              [1, 8, -2],
              [3, 2, 0]])
b = np.array([7, 6, 8])

# 正方行列をPLU分解する
(LU,piv) = linalg.lu_factor(A)
P, L, U = linalg.lu(A)

# 解を求める
x = linalg.lu_solve((LU,piv),b)

print('解:',x)
print(P@L@U)

# 解: [ 4. -2. -9.]

コレスキー分解

Aが正定値エルミート行列のときは、$A=LL^\dagger$と下三角行列$L$を用いて分解することで、さらに高速に計算できます。

A = np.array([[6, 4, 1],
              [4, 8, 2],
              [1, 2, 1]])
b = np.array([7, 6, 8])

L = linalg.cholesky(A)

t = linalg.solve(L, b)
x = linalg.solve(L.T.conj(), t)

print(x)
# 解: [ 0.45955989 -0.42470278 16.08144135]

非負値行列因子分解

ある行列AをA≒WHと近似した時に、その近似後の行列W、Hの要素が全部正になるようなもので分解する手法です。

from sklearn.decomposition import NMF

A = np.array([[1,1,1],
              [2,2,2],
              [3,3,3],
              [4,4,4]])

model = NMF(n_components=3, init='random', random_state=0)

W = model.fit_transform(A)
H = model.components_

W # 4行、n_components列の形になる。
# array([[0.08704015, 0.15523454, 0.74630123],
#       [0.88812181, 0.40340471, 0.01531961],
#       [0.48961523, 1.34083547, 0.        ],
#       [1.16745373, 1.07765565, 0.55878315]])

H # n_components行、3列の形になる。
# array([[1.46415955, 1.46434578, 1.465027  ],
#        [1.70285176, 1.70273674, 1.70231602],
#        [0.81509515, 0.81503439, 0.81481215]])

W@H
array([[1.0000886 , 1.00004161, 0.99986973],
       [1.99977738, 1.99989545, 2.00032733],
       [3.00011886, 3.00005582, 2.99982524],
       [3.99988777, 3.99994729, 4.00016501]])