この記事では, pythonによる特異値分解の方法を端的にまとめました.
特異値分解の厳密な理論については, 以下の記事を参考にしてください.
簡単に特異値分解をまとめると次のようになります. （行列の成分は実数であることを仮定します. ）
行列 \(A\) に対して, \[\begin{aligned} A \boldsymbol{v} =\sigma \, \boldsymbol{u}, \quad \boldsymbol{u}^T\, A =\sigma \, \boldsymbol{v}^T \quad （\boldsymbol{v}\neq \boldsymbol{0}, \boldsymbol{u}\neq \boldsymbol{0}） \end{aligned}\]をみたす正の数 \(\sigma\) を特異値といい, \(\boldsymbol{v}, \boldsymbol{u}\) を右特異ベクトル, 左特異ベクトルという.
特異値を大きい方から順に \[\sigma_1\geq \sigma_2\geq \cdots \geq \sigma_r\] と並べる. \(\boldsymbol{v}_i, \boldsymbol{u}_i\) は特異値 \(\sigma_i\) の特異ベクトルで, その集合\(\{\boldsymbol{v}_i\},\{\boldsymbol{u}_i\} \) が正規直交系となるようにとる. このとき, \(U=\begin{pmatrix}\boldsymbol{u}_1 &\cdots & \boldsymbol{u}_r\end{pmatrix}\) ,\(V=\begin{pmatrix}\boldsymbol{v}_1 &\cdots & \boldsymbol{v}_r\end{pmatrix}\) , \(\Sigma = \mathrm{diag} \begin{pmatrix}\sigma_1 &\cdots & \sigma_r\end{pmatrix}\) とおけば, \[\begin{aligned} A= U \Sigma V^T \end{aligned}\]と分解できる. これを特異値分解という.
※ このとき \(U, V\) は直交行列になっている.
python で特異値分解を行うには scipy の中にあるsvd関数を用います.
svd(A, full_matrices=False)
で行列 \(A\) に対して, 3つ組 U, S, VT を出力する.
U, VT は上述した行列 \(U\) と \(V^T\). Sは特異値のリスト \( (\sigma_1, \sigma_2, \cdots )\).
※VTは \(V\) の転置なので注意
例
import numpy as np
from scipy.linalg import svd
A=np.array([[1,2,4,1],[3,-1,0,4],[2,1,-2,1]])
U, S, VT=svd(A,full_matrices=False)
print(U)
print(S)
print(VT)
#出力
#[[ 0.36584702 0.89021417 0.27143082]
# [ 0.87490457 -0.22952896 -0.42644864]
# [ 0.31732939 -0.39349103 0.86282493]]
#[5.59956229 4.6943452 2.14663115]
#[[ 0.64741123 0.03109508 0.14799895 0.74698601]
# [-0.12469359 0.34434329 0.92618641 -0.08976602]
# [ 0.33435403 0.85349326 -0.29810738 -0.26624919]]実際に積をとって A に一致するか確かめてみます.
近似値なので微妙に誤差が生じることがあります.
np.dot(np.dot(U,np.diag(S)),VT)
#出力
#array([[ 1.00000000e+00, 2.00000000e+00, 4.00000000e+00,
# 1.00000000e+00],
# [ 3.00000000e+00, -1.00000000e+00, -6.44670797e-16,
# 4.00000000e+00],
# [ 2.00000000e+00, 1.00000000e+00, -2.00000000e+00,
# 1.00000000e+00]])svd関数以外にも scipy.sparse.linalg にsvds関数というのもあります. 使い方は似ていますが, オプションで, 特異ベクトルを持ってくる個数を指定できたりします.
コメント