Python特異値分解
numpy.linalg.svd()のラッピング
行列
の特異値分解（Singular Value Decomposition, SVD）は次式
で与えられる（岡本安晴「統計学を学ぶための数学入門［下］」培風館、第4.4節）。この関数を
u, s, vh = numpy.linalg.svd(A)
と実行すれば、uに左特異ベクトル、sに特異値、vhに右特異ベクトルが返される。左特異ベクトルはuに列ベクトルとして設定されているが、これは行ベクトルとした方が１つ１つを取り出すときに便利である。特異値は０も含めてsに設定されるが、０でない固有値の数と、特異値は０でないものだけを返す方が、利用し易い。これらのことを考えて、関数svdの結果を少し変えた関数svd_wを用意した。行列Aに対して次のスクリプト
r, urow, s,
vrow = svd_w(A)
を実行すると、rに０でない特異値の数、urowに左特異ベクトルを行ベクトルとして表したものをsに設定された特異値に対応する順番に並べたもの、sには特異値を大きさの降順に並べたもの、vrowには右特異ベクトルを行ベクトルとしてsに設定された特異値に対応する順番に並べたものが返される。関数svd_wをファイルsvd4me.pyに以下のように宣言した。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝ svd4me.py ＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import numpy as np
def svd_w(A):
"""
Return r, urow, s, vrow
r: Rank of A
urow: Left singular vectors as row vectors
s: Singular
value in descending order
vrow: Right sigular vectors as row vectors
"""
u, s, vh = np.linalg.svd(A,
full_matrices = False)
r =
np.linalg.matrix_rank(np.diag(s))
M, N = A.shape
urow = None
vrow = None
if r == 0:
s =
np.array([0])
urow
= np.full((1, M), 0.0)
vrow
= np.full((1, N), 0.0)
else:
s =
s[:r]
urow
= u.T[:r]
vrow
= vh[:r]
return r, urow, s, vrow
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
関数svd_wを用いたスクリプト例を以下のように用意した。階数２の
型行列Aが作成され、svd_wによって特異値分解が求められている。求めた特異値と特異ベクトルに対して、特異値分解の式（X.1）および式（X.2）に基づいて元の行列Aの再構成が行われ、正しく特異値分解が行われていることを確認している。
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import numpy as np
from svd4me import svd_w
u1 = np.array([[1, 1, 1,
1]]).transpose()
u2 = np.array([[1, 1, -1,
-1]]).transpose()
mu1 = 3.0
mu2 = 2.0
v1 = np.array([[1, 1, -1]])
v2 = np.array([[1, 0, 1]])
A = mu1 * u1 @ v1 + mu2 * u2 @ v2
print('A =\n', A)
r, urow, s, vrow = svd_w(A)
print('r = ', r)
print('urow = \n', urow)
print('s = \n', s)
print('vrow = \n', vrow)
m, n = A.shape
ckA = np.full((m, n), 0.0)
for i in range(r):
ckA += s[i] *
np.array([urow[i]]).transpose() @ \
np.array([vrow[i]])
print('ckA =\n', ckA)
ckA = urow.transpose() @ np.diag(s) @
vrow
print('ckA =\n', ckA)
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
上のスクリプトを実行すると、以下の出力が得られる。
A =
[[ 5. 3. -1.]
[ 5.
3. -1.]
[ 1.
3. -5.]
[ 1.
3. -5.]]
r = 2
urow =
[[-0.5 -0.5 -0.5 -0.5]
[-0.5 -0.5 0.5
0.5]]
s =
[10.39230485 5.65685425]
vrow =
[[-0.57735027 -0.57735027 0.57735027]
[-0.70710678 0.
-0.70710678]]
ckA =
[[ 5. 3. -1.]
[ 5.
3. -1.]
[ 1.
3. -5.]
[ 1.
3. -5.]]
ckA =
[[ 5. 3. -1.]
[ 5.
3. -1.]
[ 1.
3. -5.]
[ 1.
3. -5.]]
本ウェブサイトのPythonスクリプトファイルは、圧縮ファイルFilesSVD.zipにまとめた。
なお、特異値分解は＜岡本安晴「いまさら聞けないPythonでデータ分析――多変量解析、ベイズ統計分析（PyStan、PyMC）――」丸善出版＞でも簡潔に説明している。