特異値分解と主成分分析

R を使った特異値分解・主成分分析・バイプロットの Python 版です（ただしバイプロットについてはまだ書いていません）。Python のあやめ（iris）でも主成分分析を扱っています。

一般の $m \times n$ 行列 $A$ の特異値分解とは、$A$ を次のような形の三つの行列の積に分解することです。


$A$ = $U$ $S$ $V^T$

ここで $U$ は $m \times m$、$V$ は $n \times n$ の直交行列です。$S$ は上（または左）の正方形部分が対角行列で、残りは全部 0 です。対角要素（特異値）は左上ほど大きい値になるように並べます。

以下では $m > n$ つまり $A$ が縦長の場合だけ考えます。このとき $S$ の下側は全部 0 なので、この部分とそれに対応する $U$ の右側の部分は、なくてもかまいません。そこで、


$A$ = $U$ $S$ $V^T$

のように限定してもかまいません。以下の NumPy や SciPy の関数では full_matrices=False というオプションでこのような縮小版になります。

特異値分解は Python では numpy.linalg.svd か scipy.linalg.svd か sklearn.decomposition.TruncatedSVD で行います。

例としてあやめ（iris）のデータを考えます。

from sklearn.datasets import load_iris

iris = load_iris()
iris.data

array([[5.1, 3.5, 1.4, 0.2],
       [4.9, 3. , 1.4, 0.2],
       [4.7, 3.2, 1.3, 0.2],
       ...,
       [5.9, 3. , 5.1, 1.8]])

iris.data.shape

(150, 4)

4列はそれぞれ sepal length、sepal width、petal length、petal width を表します（単位はいずれも cm）。それぞれの平均を求めます：

iris.data.mean(axis=0)

array([5.84333333, 3.05733333, 3.758     , 1.19933333])

各列からその平均を引いたものを $A$ とします：

A = iris.data - iris.data.mean(axis=0)
A

array([[-7.43333333e-01,  4.42666667e-01, -2.35800000e+00, -9.99333333e-01],
       [-9.43333333e-01, -5.73333333e-02, -2.35800000e+00, -9.99333333e-01],
       [-1.14333333e+00,  1.42666667e-01, -2.45800000e+00, -9.99333333e-01],
       ...,
       [ 5.66666667e-02, -5.73333333e-02,  1.34200000e+00,  6.00666667e-01]])

まずは NumPy の関数を使って特異値分解してみます：

import numpy as np

U, s, Vh = np.linalg.svd(A, full_matrices=False)
U

array([[-1.06937444e-01, -5.31164840e-02,  8.17734010e-03,  1.20053534e-03],
       [-1.08133305e-01,  2.94357038e-02,  6.16531816e-02,  5.25472619e-02],
       [-1.15099407e-01,  2.41054172e-02, -5.24368218e-03,  1.05959887e-02],
       ...,
       [ 5.53860977e-02,  4.70071528e-02, -1.06310369e-01, -8.22694053e-02]])

s

array([25.09996044,  6.01314738,  3.41368064,  1.88452351])

Vh

array([[ 0.36138659, -0.08452251,  0.85667061,  0.3582892 ],
       [-0.65658877, -0.73016143,  0.17337266,  0.07548102],
       [ 0.58202985, -0.59791083, -0.07623608, -0.54583143],
       [ 0.31548719, -0.3197231 , -0.47983899,  0.75365743]])

SciPy でのやりかたもほぼ同じです：

from scipy import linalg
  
U, s, Vh = linalg.svd(A, full_matrices=False)

本当に $USV^T$ が $A$ に（誤差範囲内で）一致するか確認します：

np.allclose(U @ np.diag(s) @ Vh, A)

True

大丈夫なようです。

同じ問題を主成分分析した結果との対応を見てみましょう：

from sklearn.decomposition import PCA

model = PCA(n_components=2)
x = model.fit_transform(iris.data)
x

array([[-2.68412563,  0.31939725],
       [-2.71414169, -0.17700123],
       [-2.88899057, -0.14494943],
       ...,
       [ 1.39018886, -0.28266094]])

これは特異値分解の $US$ の最初の2列に一致します（ただし2列目の符号がすべて逆です）：

U @ np.diag(s)

array([[-2.68412563e+00, -3.19397247e-01,  2.79148276e-02,  2.26243707e-03],
       [-2.71414169e+00,  1.77001225e-01,  2.10464272e-01,  9.90265503e-02],
       [-2.88899057e+00,  1.44949426e-01, -1.79002563e-02,  1.99683897e-02],
       ...,
       [ 1.39018886e+00,  2.82660938e-01, -3.62909648e-01, -1.55038628e-01]])

主成分分析の回転行列は

model.components_

array([[ 0.36138659, -0.08452251,  0.85667061,  0.3582892 ],
       [ 0.65658877,  0.73016143, -0.17337266, -0.07548102]])

ですが、これはさきほど求めた特異値分解の Vh の最初の2行と一致します（ただし2行目の符号がすべて逆です）。このように、$U$ の列と $V^T$ の行の符号には不定性があります。