Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

はじめに

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

参考

手順

特異値分解

$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]])

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away