Python特異値分解
行列
の特異値分解（Singular Value Decomposition, SVD）については、岡本（2009）、pp. 47-49、あるいは岡本（2019）、pp. 67-72で説明している。
で与えられる。特異値分解はumpyの関数numpy.linalg.svdにより得ることができる。
リスト１にサンプルスクリプトを示す。スクリプトファイルは圧縮ファイルFilesSVD.zipとして用意した。ダウンロード展開すれば自由に利用できる。
特異値分解は次のコードで行っている。
U, S, Vh = np.linalg.svd(a, full_matrices = False)
行列aの特異値分解が、U、S、Vhとして得られる。パラメータfull_matrices = Falseについては、本ウェブサイトの後半で説明している。
特異値Sは降順で並べられている。リスト１では、０であるかどうかの判定を次のコードで行っている。
ck_idx = S > 1.0e-7
特異値が０でないものを集めて、次のコードで特異値分解が得られていることを確認している。
ck_a = U @ np.diag(S) @ Vh
print('ck_a =\n', ck_a)
import numpy as np
a =
np.array([[1,2,3,4],[8,7,6,5],[1,1,1,1]])
print('a =\n', a)
U, S, Vh = np.linalg.svd(a,
full_matrices = False) # 特異値分解
print('U =\n', U)
print('S =\n', S)
print('Vh =\n', Vh)
ck_idx = S > 1.0e-7
# ０でないと見做せるか？
print('ck_idx =\n', ck_idx)
U = U[:,ck_idx]
# 特異値０に対応するものを除く
print('U =\n', U)
S = S[ck_idx]
# 特異値０に対応するものを除く
print('S =\n', S)
Vh = Vh[ck_idx]
# 特異値０に対応するものを除く
print('Vh =\n', Vh)
ck_a = U @ np.diag(S) @ Vh
# 特異値分解の確認
print('ck_a =\n', ck_a)
リスト１を実行すると、以下の出力を得る。
a =
[[1 2 3 4]
[8 7 6 5]
[1 1 1 1]]
U =
[[ 0.33741061 0.9349363 -0.10976426]
[ 0.93075201 -0.34878718 -0.10976426]
[ 0.14090696 0.06512768 0.98787834]]
S =
[1.41310647e+01 2.88322914e+00
1.74222008e-16]
Vh =
[[ 0.560774 0.51878556 0.47679711 0.43480867]
[-0.62091265 -0.17567455 0.26956356 0.71480166]
[ 0.26126424 0.01045822 -0.80470915 0.53298669]]
ck_idx =
[ True True False]
U =
[[ 0.33741061 0.9349363 ]
[ 0.93075201 -0.34878718]
[ 0.14090696 0.06512768]]
S =
[14.13106471 2.88322914]
Vh =
[[ 0.560774 0.51878556 0.47679711 0.43480867]
[-0.62091265 -0.17567455 0.26956356 0.71480166]]
ck_a =
[[1. 2. 3. 4.]
[8. 7. 6. 5.]
[1. 1. 1. 1.]]
リスト１におけるSVDのパラメータfull_matricesの比較を行うスクリプトをリスト２に用意した。
次のコードは、リスト１と同じくfull_matrices = Falseと設定されている。
U, S, Vh = np.linalg.svd(a, full_matrices = False)
これに対して、比較のために次のコードも用意した。
Uc, Sc, Vhc = np.linalg.svd(a, full_matrices = True)
パラメータがfull_matrices = Falseのとき、式（１）に合わせた形でUとVhが得られるが、パラメータの指定がfull_matrices = Trueのときは正方行列の形（必要ならば拡張して）でUcとVhcが得られる。この正方行列（直交行列）への拡張は、固有分解の場合について岡本（2009）、pp. 105-106で論じられているが、参考になると思う。
import numpy as np
a =
np.array([[1,2,3,4],[8,7,6,5],[1,1,1,1]])
print('a =\n', a)
U, S, Vh = np.linalg.svd(a,
full_matrices = False) # 特異値分解
print('U =\n', U)
print('S =\n', S)
print('Vh =\n', Vh)
Uc, Sc, Vhc = np.linalg.svd(a,
full_matrices = True)
# 特異値分解
print('Uc =\n', Uc)
print('S =\n', Sc)
print('Vhc =\n', Vhc)
ck_Vh = Vhc.T @ Vhc
print('ck_Vh =\n', ck_Vh)
リスト２を実行すると、以下の出力を得る。
a =
[[1 2 3 4]
[8 7 6 5]
[1 1 1 1]]
U =
[[ 0.33741061 0.9349363 -0.10976426]
[ 0.93075201 -0.34878718 -0.10976426]
[ 0.14090696 0.06512768 0.98787834]]
S =
[1.41310647e+01 2.88322914e+00
1.74222008e-16]
Vh =
[[ 0.560774 0.51878556 0.47679711 0.43480867]
[-0.62091265 -0.17567455 0.26956356 0.71480166]
[ 0.26126424 0.01045822 -0.80470915 0.53298669]]
Uc =
[[ 0.33741061 0.9349363 -0.10976426]
[ 0.93075201 -0.34878718 -0.10976426]
[ 0.14090696 0.06512768 0.98787834]]
S =
[1.41310647e+01 2.88322914e+00
1.74222008e-16]
Vhc =
[[ 0.560774 0.51878556 0.47679711 0.43480867]
[-0.62091265 -0.17567455 0.26956356 0.71480166]
[ 0.26126424 0.01045822 -0.80470915 0.53298669]
[-0.48139485 0.83659466 -0.22900478 -0.12619503]]
ck_Vh =
[[ 1.00000000e+00 -1.11022302e-16 -3.05311332e-16
-8.32667268e-17]
[-1.11022302e-16 1.00000000e+00 -5.55111512e-17
-9.71445147e-17]
[-3.05311332e-16 -5.55111512e-17 1.00000000e+00 -1.70002901e-16]
[-8.32667268e-17 -9.71445147e-17
-1.70002901e-16 1.00000000e+00]]
解VhcがVhの直交行列への拡張であることが確認できる・
引用文献
岡本安晴（2009）統計学を学ぶための数学入門［下］．培風館
岡本安晴（2019）いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――．丸善出版