写真から立体を再現（7）　実際に特異値分解(SVD)してみる

　Numpyの特異値分解は、結果の解釈に注意が必要です。説明します。
　行列Aは、以下のように特異値分解（SVD）されるのでした。
f:id:Inuosann:20200318223038p:plain:w350
過去の記事にあります。
www.omoshiro-suugaku.com
Pythonのモジュール、Numpyの特異値分解では上のようでなく、別の形の結果が得られます。まず、区分行列の計算で以下が成立することを確認してください。ここで r は A のランクです。
f:id:Inuosann:20200318223059p:plain:w400
Numpyではこの形で特異値分解が得られます。＊のところには何か数値が入りますが、区分行列のかけ算の結果、ここに何が入っても結果は変わりません。
　実際に試してみましょう。この例では、各サイズは m = 2, n = 3, r = 2 ということになります。

import sys
import numpy as np
A = np.array([[3, 1, 2], [3, 2, 1]])
U, S, VT = np.linalg.svd(A)
print('A = \n', A)
print('rank =  ', np.linalg.matrix_rank(A))
U, S, VT = np.linalg.svd(A) #特異値分解
print('U = \n', U)
print('S = \n', S)
print('VT = \n', VT)
W = VT[0:2, 0:3] #VTの1行目と2行目を取り出す
B = np.dot(np.dot(U, np.diag(S)), W) #かけ算して元に戻るか、確認
print('B = \n', B)
sys.exit()

実行した結果は次の通りです。
A =
[[3 1 2]
[3 2 1]]
rank = 2
U =
[[-0.70710678 -0.70710678]
[-0.70710678 0.70710678]]
S =
[5.19615242 1. ]
VT =
[[-8.16496581e-01 -4.08248290e-01 -4.08248290e-01]
[ 1.94289029e-16 7.07106781e-01 -7.07106781e-01]
[-5.77350269e-01 5.77350269e-01 5.77350269e-01]]
B =
[[3. 1. 2.]
[3. 2. 1.]]

区分行列の式から、svd() の結果の U はそのまま使えます（n - r = 0 に注意）。Sは特異値（正）のリストが求まるだけなので、diag(S) として対角行列に直します。また VT から1，2行目を取り出せばよいことが分かります(r = 2 に注意)。コードでは1，2行目を取り出して作った行列をWとしています。つまり、求める特異値分解は
A = U * diag(S) * W
ということになります。念のためかけ算して元に戻っていることが分かります。なお、svd() には今回のコードでは使っていないパラメータがあり、それによっては上とは異なる形の結果が得られます。
　行列の例は何度も紹介している以下の本に載っていたものを使わせていただきました。