Up

Python特異値分解

numpy.linalg.svd()のラッピング

 

行列の特異値分解(Singular Value Decomposition, SVD)は次式

で与えられる(岡本安晴「統計学を学ぶための数学入門[下]」培風館、第4.4節)。この関数を

u, s, vh = numpy.linalg.svd(A)

と実行すれば、uに左特異ベクトル、sに特異値、vhに右特異ベクトルが返される。左特異ベクトルはuに列ベクトルとして設定されているが、これは行ベクトルとした方が1つ1つを取り出すときに便利である。特異値は0も含めてsに設定されるが、0でない固有値の数と、特異値は0でないものだけを返す方が、利用し易い。これらのことを考えて、関数svdの結果を少し変えた関数svd_wを用意した。行列Aに対して次のスクリプト

r, urow, s, vrow = svd_w(A)

を実行すると、rに0でない特異値の数、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を用いたスクリプト例を以下のように用意した。階数2の型行列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――丸善出版>でも簡潔に説明している。

 

Up