「線形代数の半歩先」をPythonで写経 ～ ３章 特異値分解、固有値分解、主成分分析
第３部「ならべた数に応用を」
書籍の著者 大久保 潤 先生
この記事は、書籍「線形代数の半歩先」の 第３部「ならべた数に応用を」に掲載の「データサイエンスと機械学習」に関する Python写経活動 のドキュメンタリーです。
第３部は機械学習でおなじみの「回帰」と「次元削減」を線形代数と接続しています。
機械学習に慣れた方は書籍の内容がすっと腹落ちするのではないでしょうか？
この記事では、第17話、第18話の 特異値分解、固有値分解、主成分分析 の計算に取り組みます！
数学素人なので、どうぞお手柔らかにお願いいたします。
では書籍を開いて線形代数の旅に出発です🚀
はじめに
このブログシリーズは、書籍「線形代数の半歩先 データサイエンス・機械学習に挑む前の30話」（講談社サイエンティフィク、「テキスト」と呼びます）の Python 写経の実践を通じて得た個人的な知見を書きます。
書籍の紹介と引用表記はリンク先の記事に掲載しています。
第３部 ならべた数に応用を
テキストの第17話「行列の特別な分解」と第18話「直交の技術」は、主成分分析からはじまります。
特異値分解や固有値分解といった行列分解へ進んで、行列の特別な側面を学びます。
では、この記事のコードで利用するライブラリのインポートからスタートしましょう。
（注）この記事では可視化はありません！
### インポート
# 数値計算
import numpy as np
# 主成分分析(PCA)
from sklearn.decomposition import PCA１．特異値分解
① 特異値分解の概要
テキスト第17話「行列の特別な分解」に進みます。
ある行列を別の行列やベクトルの積で表現するのが行列の分解の醍醐味です。
特異値分解ではデータ行列 $${X}$$ を $${U, \Sigma, V}$$ の３つの要素に分解します。
$$
X = U \Sigma V^{\top}
$$
$${U}$$
左特異ベクトルと呼ばれる列ベクトルで構成される行列
（直交行列であり、各列・各行が正規直交基底となる）$${\Sigma}$$
特異値 $${\sigma}$$ で構成される対角行列$${V}$$
右特異ベクトルと呼ばれる列ベクトルで構成される行列
（直交行列であり、各列・各行が正規直交基底となる）
各行列の列は特異値の大きな順序で並んでいます。
各行列の行・列の数合わせが複雑です。
ぜひ、テキストの図をご参照ください！
② 特異値分解の実装
テキストの図 17.3「特異値分解で軸を選ぶ」に示された特異値分解の具体例を実装します。
行列の形や数値を見ることで、特異値分解のイメージが広がるかと。
ところでテキストの第１刷・第２刷では、図 17.3 が誤植になってます。
次のサポートページで訂正情報をご確認ください。
では実装しましょう！
データ行列 $${X}$$ を中心化した $${\widetilde{X}}$$（X_tilde）の特異値分解を numpy.linalg の svd() で実行します。
中心化は、説明変数の各値から当該説明変数の平均値を差し引いて、説明変数の平均値をゼロにするデータ変換のことです。
テキストの式（17.5）では $${\widetilde{X} = \widetilde{U} \widetilde{\Sigma} \widetilde{V}^{\top}}$$ と示されています。
### p.148 図17.3 特異値分解で軸を選ぶ
# [1] データセットX
X = np.array([[8, 6, 7, 0],
[7, 8, 8, 0],
[9, 8, 7, 0],
[0, 0, 1, 9],
[0, 0, 0, 8]])
print('[1] データセットX:')
print(X)
# ランク R の設定
R = 4
# [2] 中心化 X_tilde
X_tilde = X - X.mean(axis=0)
print('[2] X(中心化後):')
print(X_tilde)
# [3] 特異値分解
print('-'*10, '特異値分解', '-'*10)
U, S, VT = np.linalg.svd(X_tilde)
# U, Σ, V.Tの結果表示
print('U:')
print(U[:, :R].round(3))
print('Σ:')
print(np.diag(S[:R]).round(3)) # 対角行列に変換
print('V.T')
print(VT[:R, :].round(3))【実行結果】
numpy を使うことであっさり求まりました！
■ コードの補足説明
データの中心化
データ行列 $${X}$$ から $${X}$$ の列ごとの平均値を X.mean(axis=0) で計算して、引き算します。
mean の引数 axis=0 で「列の平均」を計算できます。
X_tilde = X - X.mean(axis=0)特異値分解
numpy.linalg の svd() に 中心化済みのデータ行列 $${\widetilde{X}}$$（X_tilde）を与えて、X_tilde の特異値分解を実行します。
計算結果（戻り値）は $${U}$$、$${\Sigma}$$、$${V^{\top}}$$ の順で取得できます。
$${\Sigma}$$ はテキストの対角行列の形式ではなく、平坦な配列の形式になっています。
$${V}$$ は転置済みの $${V^{\top}}$$ で取得します。
U, S, VT = np.linalg.svd(X_tilde)◆ ◆ ◆
③ 特異値分解と主成分分析の関係
続いて特異値分解の結果と主成分分析の関係を探りましょう。
テキストの 149 ページに次のことが書かれています。
$${\widetilde{V}^{\top}}$$ の行ベクトルの１行目（転置前だと列ベクトルの１列目）が、主成分分析で選ばれる最初の軸の係数を与えます。
$${\widetilde{U}}$$ の１列目は、主成分分析の最初の軸での各データの「座標」と関係します。
まず１点目の $${V}$$ と主成分分析の最初の軸の係数の関係を見てみましょう。
scikit-learn の PCA を用いてデータ行列 $${\widetilde{X}}$$ の主成分分析を実行します。
### sklearnのPCAで係数ベクトル（固有ベクトル）と特異値を取得
# PCAの実行
pca = PCA(n_components=R)
pca.fit(X_tilde)
# 主成分（固有ベクトル）の表示 ★特異値分解のV.Tと一致する
print('主成分.T:') # shape=(主成分の数R, 特徴量の数4)
print(pca.components_.round(3))
# 特異値の表示 ★特異値分解のΣと一致する
print('特異値:')
print(pca.singular_values_.round(3))
# 参考：主成分の分散（固有値）の表示
print('参考：主成分の分散（固有値）:')
print(pca.explained_variance_.round(3))【実行結果】
主成分.T で示したデータ行列 $${\widetilde{X}}$$ の主成分の転置は、特異値分解の $${\widetilde{V}^{\top}}$$ と一致しています！
$${V}$$ は主成分分析の軸の係数であることが分かりました。
■ コードの補足説明（主成分分析）
scikit-learn の PCA を利用します。
pca = PCA(n_components=R)
pca.fit(X_tilde)１行目：PCAの器 pca を作ります。引数は主成分数です。
２行目：主成分分析の実行です。X_tilde の主成分等を算出します。
器 pca に「. xxx」の形式でアクセスして、主成分分析のさまざまな数値を取得します。
主成分の取得：.components_
pca.components_特異値の取得：.singular_values_
pca.singular_values_分散の取得：.explained_variance_
pca.explained_variance_◆ ◆ ◆
続いて$${U}$$ と主成分分析の軸での各データの「座標」の関係を見てみましょう。
テキストの色付きエリアに「特異値を$${U}$$に掛け算したものが、それぞれのデータの座標に対応」するとあります。
掛け算しましょう。
### p.149「座標としての解釈の補足」にあるUの補正
print('特異値をUに掛け算したものがそれぞれのデータの座標に対応する')
print((U[:, :R] * S[:R]).round(3))【実行結果】
この数値群は主成分分析界隈では「主成分得点」と呼ばれる量に相当します。
主成分得点を計算してみましょう。
PCA で計算したものと、主成分$${\times}$$データで計算したもの（ときどき$${z}$$で表現されるもの）を並べます（両者同じです。）
### 上記の「Uの補正」はいわゆる「主成分得点」
# PCAで主成分得点を計算
score = pca.transform(X_tilde)
print('PCAで主成分得点を計算：')
print(score.round(3), '\n')
# sum(主成分@データ）で計算
print('sum(主成分×データ)で計算：')
print((X_tilde @ pca.components_.T).round(3))【実行結果】
さきほどの「掛け算の結果」と一致しました！
行がデータ点、列が主成分（左から順に第１主成分, …, 第４主成分）、値が主成分得点です。
主成分得点は、各データ点が各主成分の軸上のどの位置にあるかを示す指標です。
■ コードの補足説明（主成分分析の続き）
器 pca に対して変換操作「.transform」をすることで、引数で与えた X_tilde の主成分得点を取得できます。
score = pca.transform(X_tilde)２．低ランク近似
① 低ランク近似の概要
さきほどの特異値分解では、ランク＝ $${4}$$ を指定しました。
特異値が４つ並んでいますよね！
データは特異値の大きな順に並んでいて、最初から２つの特異値がとても大きく「最初の２つの影響が大きく残る」だろうとテキストは説明します。
そして、テキストの図 17.4「当たらずとも⋯遠からず？」では最初の２つの特異値に関係する$${U,\Sigma,V}$$ を使って、データ行列 $${X}$$ の「近似」を探ります。
特異値の数を減らして元のデータ行列に近いデータを得られるなら、データ圧縮が可能になります！
② 低ランク近似の実装
では図 17.4 を実装しましょう。
なお、こちらの特異値分解ではデータを中心化していません。
まずはデータ設定と特異値分解から。
### 図17.4 低ランク近似
# データセットX
X = np.array([[8, 6, 7, 0],
[7, 8, 8, 0],
[9, 8, 7, 0],
[0, 0, 1, 9],
[0, 0, 0, 8]])
print('データセットX:')
print(X)
# 特異値分解
U, S, VT = np.linalg.svd(X)
# U, Σ, V.Tの結果表示
print('U:')
print(U.round(2))
print('Σ:')
print(np.diag(S).round(2))
print('V.T')
print(VT.round(2))【実行結果】
図 17.4 の上段の各値と「絶対値」が合っているかと思います。
続いてランク２の低ランク近似です。
図 17.4 の下段に示された行列の積を用いて「再構成」をします。
また再構成後の行列のランクを numpy.linalg の matrix_rank() で求めています。
### 低ランク近似による再構成１ 図17.4に示された行列計算
# ランクの指定
R = 2
# 低ランク近似X_rk2の算出
X_rk2 = (U[:, :R] @ np.diag(S[:R]) @ VT[:R, :])
# 結果の表示
print('ランク:', np.linalg.matrix_rank(X_rk2))
print(X_rk2.round(2))【実行結果】
確かに行列のランクは２になりました。
データ行列 $${X}$$ とかなり近い感じがします！
もう一つの低ランク近似をテキストの式（17.8）で計算したいと思います。
もちろん、上の低ランク近似の結果と一致します！
$$
X = \sum_{r=1}^R \sigma_r \bm u_r \bm v_r^{\top} = \sum_{r=1}^R \sigma_r \ket{\bm u_r} \bra{\bm v_r}
$$
### 低ランク近似による再構成２ 式（17.8）を利用
# ランクの指定
R = 2
# ヘルパー関数の定義
vec = lambda x: np.array(x).reshape(-1, 1)
# 低ランク近似X_rk2の算出 式（17.8）
X_rk2_2 = np.sum(
[S[r] * vec(U[:, r]) @ vec(VT[r, :]).T for r in range(R)], axis=0)
# 結果の表示
print('ランク：', np.linalg.matrix_rank(X_rk2_2))
print(X_rk2_2.round(2))【実行結果】
３．正規性・直交性の確認
テキスト第18話「直交の技術」に進みます。
固有ベクトルが直交していること、正規直交基底になっていることを踏まえて、いろんな議論を繰り広げます。
正規とはベクトルの長さ（ノルム）が１であること、直交とはベクトル間の内積が０であること、です。
ということで「主成分分析」に関連する固有ベクトルのノルムと内積を計算しましょう。
主成分分析の固有値・固有ベクトルは、データ行列の「共分散行列」に対する固有値分解で求めることもできます。
今回は共分散行列を用いて固有値・固有ベクトルを計算しましょう。
前に作成した中心化したデータ行列 $${\widetilde{X}}$$（X_tilde）を利用して、共分散行列を計算し、固有値・固有ベクトルを求めます。
### データの共分散行列を 1/(N-1)･XᵀX で算出して固有値分解する 参考：式（18.20）
# データ行列のデータ点の数の取得
N = X_tilde.shape[0]
# データ行列の共分散行列の算出
cov = 1 / (N - 1) * X_tilde.T @ X_tilde
print('共分散行列:')
print(cov)
# 共分散行列の固有値・固有ベクトルの算出
eigval, eigvec = np.linalg.eig(cov)
sort_index = eigval.argsort()[::-1]
# 結果の表示
print('固有値:')
print(eigval[sort_index].round(3))
print('固有ベクトル:')
print(eigvec[:, sort_index].round(3))【実行結果】
さきほどの PCA で求めた固有値と一致し、主成分の転置した行列と固有ベクトルが一致しています。
では固有ベクトルの長さ・ノルムを計算しましょう。
eigvec の４つの列ベクトルそれぞれが固有ベクトルです。
### p.156 固有ベクトルのノルムは１（正規直交基底なのです）
print('固有ベクトルのノルム:')
print(np.linalg.norm(eigvec, axis=0))【実行結果】
４つの固有ベクトルのノルムは１です！
正規性があります！
続いて固有ベクトル間の内積を計算して、直交性を確認しましょう。
直交するベクトル同士の内積は０になります。
### p.156 固有ベクトル間の内積は０（正規直交基底なのです）
# 固有ベクトルの内積を計算
dots = np.array([
[eigvec[:, i].T @ eigvec[:, j] for j in range(4)] for i in range(4)])
# 結果の表示
print('固有ベクトルの内積:')
print(dots.round(10))【実行結果】
対角成分は自分自身との内積なので１です（ノルムの二乗＝内積＝１）。
その他の成分は異なる固有ベクトル間の内積であり、すべて０です！
直交性を確認できました。
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！