🕌

numpyで特異値分解

2024/02/04に公開

「Pythonによる数値計算入門」で最小二乗法をnumpy.linalg.solveで解く例があったので、それのsvd版を書いておきます。

import numpy as np

# y = 2x + 3に誤差を加えたデータ
# 数値はPythonによる数値計算入門より

N = 5
X = np.zeros(N)
Y = np.zeros(N)
Z = np.zeros(N)

X[0] = 0.0; Y[0] = 3.00
X[1] = 1.0; Y[1] = 4.51
X[2] = 2.0; Y[2] = 7.54
X[3] = 3.0; Y[3] = 8.43
X[4] = 4.0; Y[4] = 11.29

Z += 1

X, -Y, Z

(array([0., 1., 2., 3., 4.]),
array([ -3. , -4.51, -7.54, -8.43, -11.29]),
array([1., 1., 1., 1., 1.]))

# 行方向に積む
A = np.vstack((X, -Y, Z)).T
A

array([[ 0. , -3. , 1. ],
[ 1. , -4.51, 1. ],
[ 2. , -7.54, 1. ],
[ 3. , -8.43, 1. ],
[ 4. , -11.29, 1. ]])

# 特異値分解
U, S, V = np.linalg.svd(A, full_matrices=True)
U, S, V

(array([[-0.16585796, 0.76990057, -0.05739962, -0.5002555 , -0.3552379 ],
[-0.26295701, 0.41604903, 0.47881684, 0.17075764, 0.70663504],
[-0.44082209, 0.24508048, -0.57009142, 0.64387872, -0.07763639],
[-0.5049771 , -0.18336807, 0.61268906, 0.20099165, -0.54368074],
[-0.67380914, -0.37479064, -0.25893552, -0.5153725 , 0.26991999]]),
array([17.81907654, 1.46097052, 0.25841126]),
array([[-0.30050766, 0.94682632, -0.11495676],
[-0.78239553, -0.1757807 , 0.59745994],
[ 0.54548362, 0.26948295, 0.79361613]]))

a, b, c = V[-1]  # 最小の特異値に対応するベクトル == Vの最後の行が解
a / b, b / b, c / b  # yの係数が1になるように調整

(2.024186046157259, 1.0, 2.944958643197044)

もともとの式が2x + 3なので、妥当そうな値ですね。

X * a / b + c / b  # 検算

array([ 2.94495864, 4.96914469, 6.99333074, 9.01751678, 11.04170283])

与えられたYのデータと見比べても妥当そうです。

colabのリンクも貼っておきます。

参考文献

Pythonによる数値計算入門 (実践 Pythonライブラリー)

参考文献

Discussion