気象情報解析特論 第7回 特異値分解（理論編）¶
担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
特異値分解(SVD)とは¶
特異値分解は，正方行列$X$について行った固有値分解
$$XE = E\Lambda \iff X = E\Lambda E^T$$の，非正方行列にも適用できる一般化であると考えることができます。
つまり，$M$行$N$列の実行列$X$について，特異値分解(Singular Value Decomposition; SVD)を行うとは，
$$XV = U\Sigma \iff X = U \Sigma V^T$$となるような直交行列$U$（$M$行$M$列）と$V$（$N$行$N$列）と，対角成分が非負であり非対角成分が0である行列$\Sigma$（$M$行$N$列）を求めることです。 このとき，$U$と$V$の列ベクトルをそれぞれ左特異ベクトルと右特異ベクトルといい，$\Sigma$の対角成分は普通左上から降順で並べて特異値と言います。
（そのような分解が存在することの証明らしいです（すみませんが神山は追ってないです） https://risalc.info/src/svd.html ）
（本当は色々と等価な分解方法がありえますが，本資料ではこのように定義します。詳しくは，たとえばhttps://qiita.com/kidaufo/items/0f3da4ca4e19dc0e987e あたりでしょうか）
ちなみに$X$が正方行列の際は，$E=U=V,\Lambda = \Sigma$となり，固有値分解と一致します。
特異ベクトルと特異値の求め方¶
特異ベクトルと特異値を手計算で求める方法のうち，もっともわかりやすいのは，固有値展開を使う方法です。
特異ベクトルの求め方：$U$は，$XX^T$の単位固有ベクトルを固有値の大きい順に左から並べた行列であり，$V$は$X^T X$の単位固有ベクトルを左から並べた行列として求めることができます（注：固有ベクトルの符号はうまくいくように選ぶ必要があります→D問題）。
$X = U \Sigma V^T$ゆえ，$X^T=V \Sigma U^T$です（直交行列の性質$U^{-1} = U^T$, $V^{-1} = V^T$に注意）。よって，
$$X X^T = (U \Sigma V^T)(V \Sigma^T U^T) = U \Sigma \Sigma^T U^T \ \ \ \ (1)$$であるので，確かに$U$は$XX^T$の固有ベクトルを並べた行列です。$XX^T$は実対称行列ですので，固有ベクトルを並べた行列が直交行列となることとも整合します（復習：線形代数学4）。$V$についても，
$$X^T X = (V \Sigma^T U^T)(U \Sigma V^T) = V \Sigma^T \Sigma V^T \ \ \ \ (2)$$ですので，$U$と同様です。
特異値の求め方：特異値は，$XX^T$と$X^TX$に共通する非ゼロ固有値の平方根として求めることができます。
前回の授業で，$X$を転置してから主成分分析を行っても，数学的に等価な結果が得られることがわかりました。$XX^T$と$X^TX$というのは，共分散行列（に定数をかけたもの）そのものですから，$XX^T$と$X^TX$の非ゼロ固有値は等しいということになります。
式(1)と(2)から，$XX^T$と$X^TX$の固有値はそれぞれ対角行列$\Sigma \Sigma^T$と$\Sigma^T \Sigma$の対角成分ですので，確かにこれも大丈夫そうです。
X = np.array([[1, 2, 3], [4, 5, 6]])
print(X)
$X$を特異値分解します。
U, S, VT = np.linalg.svd(X) #特異値分解，Sは特異値が降順に入っているベクトル
V = VT.T #Vの転置
Sigma = np.zeros(X.shape) # Xと同じ大きさのfloat64型行列Sigmaを用意
for ss in range(0, S.shape[0]):
Sigma[ss, ss] = S[ss] #特異値を対角行列に並べる
では，実際に特異値分解ができているかどうかを確かめてみましょう。
X_SVD = U@Sigma@V.T # 行列の積は@で書ける
print(X_SVD)
ちゃんと展開できていますね。
続いて，左特異ベクトルを並べた行列$U$が，$XX^T$の固有ベクトルを並べた行列$E$と等しいことを確かめてみます。
print(U)
## X@X.Tを固有値展開する
[Lam, E] = np.linalg.eig(X@X.T)
# なぜか中身は実数なのに複素数型にキャストされる（未解決...）ので，実数型に直す
Lam = np.real(Lam) # 固有値の入ったベクトル
E = np.real(E) # 固有ベクトルが入った行列
# 固有値の大きい順にソート
index = Lam.argsort()[::-1]
Lam = Lam[index]
E = E[:, index] # PC時系列を並べた行列
# 対角成分に固有値を並べる
L = np.diag(Lam) # 対角成分に大きい順に固有値の入った行列
###必要ならば，固有ベクトルの符号はそれぞれ自由につけ変えて良い
##第1主成分の基底
#e1 = -E[:, 0]
#E[:, 0] = e1
##第2主成分の基底
#e2 = -E[:, 1]
#E[:, 1] = e2
print(E)
$U$と$E$も一致しています！
最後に，特異値が$XX^T$の非ゼロ固有値の平方根と等しいことを確かめてみましょう。
print(Sigma)
print(np.sqrt(L)) # np.sqrtは，要素ごとに平方根を取ってくれる関数
こちらも大丈夫そうです。
気象データ解析の文脈では何を意味するか¶
特異値分解が気象データ解析の文脈で出てくるのは，主に以下の3通りの場合です。
東京と宇都宮の気温を例に，説明してみましょう。
東京の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv
宇都宮の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Utsunomiya_temp.csv
（出典は第3回に記述）
## 下準備
tokyo_temp = np.genfromtxt("Tokyo_temp.csv", # ファイルのパスを書く
delimiter=",", # 区切り文字
usecols=(0, 1, 2) # 読み込みたい列番号
)
utsu_temp = np.genfromtxt("Utsunomiya_temp.csv", # ファイルのパスを書く
delimiter=",", # 区切り文字
usecols=(0, 1, 2) # 読み込みたい列番号
)
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
tokyo = tokyo_temp[:, 2]
utsu = utsu_temp[:, 2]
# 今回は1990年から2019年の30年分のデータを用いる
tokyo = tokyo[(1990 <= y)*(y <= 2019)]
utsu = utsu[(1990 <= y)*(y <= 2019)]
m = m[(1990 <= y)*(y <= 2019)]
y = y[(1990 <= y)*(y <= 2019)]
# 気候値の計算
tokyoc= np.zeros((12))
utsuc= np.zeros((12))
for mm in range(1, 13):
tokyoc[mm-1] = np.nanmean(tokyo[m==mm], 0)
utsuc[mm-1] = np.nanmean(utsu[m==mm], 0)
# 偏差の計算
tokyoa = np.zeros((tokyo.shape))
utsua = np.zeros((utsu.shape))
for yy in range(1990, 2020):
for mm in range(1, 13):
tokyoa[(y==yy)*(m==mm)] = tokyo[(y==yy)*(m==mm)] - tokyoc[mm-1]
utsua[(y==yy)*(m==mm)] = utsu[(y==yy)*(m==mm)] - utsuc[mm-1]
# 今回は温暖化には興味がないのでデトレンド
tokyoa = signal.detrend(tokyoa)
utsua = signal.detrend(utsua)
# データ行列の定義
X = np.array([tokyoa, utsua]) #データ行列
1. 共分散行列$C = XX^T/N$を特異値分解すれば，$U$も$V$もEOFを並べた行列¶
まずもっとも簡単な場合として，共分散行列自体を特異値分解する場合，それは主成分分析そのものです。
なぜなら，共分散行列は正方行列ゆえ，特異値分解と固有値分解は一致するからです。
N = tokyoa.shape[0] #データ行列の大きさ
C = X@X.T/N # 共分散行列の計算
まず，固有値分解しています。
## 固有値分解
[Lam, E] = np.linalg.eig(C)
# 固有値の大きい順にソート
index = Lam.argsort()[::-1]
Lam = Lam[index]
E = E[:, index] # PC時系列を並べた行列
L = np.diag(Lam) # 対角成分に大きい順に固有値の入った行列
print(E)
print(L)
次に，特異値分解してみます。
## 特異値分解
U, S, VT = np.linalg.svd(C) #特異値分解，Sは特異値が降順に入っているベクトル
V = VT.T #Vの転置
Sigma = np.zeros(C.shape) # Xと同じ大きさのfloat64型行列Sigmaを用意
for ss in range(0, C.shape[0]):
Sigma[ss, ss] = S[ss] #特異値を対角行列に並べる
print(U)
print(Sigma)
print(V) # U=Vになっている
確かに，どちらも同じ結果になりました！
2. データ行列Xそのものを特異値分解すると，$U$はEOFを，$V$はPC時系列を並べた行列¶
次に，共分散行列ではなくデータ行列そのものを特異値分解するとどうなるでしょうか。
## 特異値分解
U, S, VT = np.linalg.svd(X) #特異値分解，Sは特異値が降順に入っているベクトル
V = VT.T #Vの転置
Sigma = np.zeros(X.shape) # Xと同じ大きさのfloat64型行列Sigmaを用意
for ss in range(0, X.shape[0]):
Sigma[ss, ss] = S[ss] #特異値を対角行列に並べる
print(U)
まず，ここで出てくる$U$は，先ほどすでに求めた共分散行列$C = XX^T/N$の固有ベクトル（=EOF）を並べた行列と等しくなります($U=E$)。特異値分解の性質として，「左特異ベクトルは$XX^T$の固有ベクトル」でしたのでこれは当然です。
発展的な注：どうせ固有ベクトルは単位行列として求めますので，共分散行列にかかるファクター$1/N$は無視して大丈夫ということになります。
では，$V$は何を意味するかというと，これはPC時系列です。固有値分解で求めたPC時系列は，EOFに元のデータを射影した時系列$Z = E^T X$ですので，
$$Z = E^T X = E^T U\Sigma V^T = E^T E \Sigma V^T = \Sigma V^T$$という関係になります。よって$Z$は，$V$の$i$番目(i = 1, 2, ..., rank(X))の列ベクトルに$i$番目の特異値のスカラー倍をかけた行列であるといえます。
発展的な注：行列の非ゼロ特異値の数は，その行列のランクと一致します。
すなわち，$Z$と$V$は，標準偏差等で規格化して（そして符号も合わせて）しまえば，どちらも同じPC時系列を与えます。
# EOFに元のデータを射影して得られた，正規化PC1時系列
Z = E.T@X
pc1_eig = Z[0, :]
pc1_eig = pc1_eig/np.std(pc1_eig)
# 右特異ベクトルとして得られた，正規化PC1時系列
pc1_svd = V[:, 0]
pc1_svd = pc1_svd/np.std(pc1_svd)
# 固有ベクトルに元のデータを射影して得られたPC1時系列を青で描画
plt.plot(pc1_eig, 'b')
plt.show()
# 右特異ベクトルとして得られたPC1時系列を赤で描画
plt.plot(pc1_svd, 'r')
plt.show()
確かに全く同じ時系列になっていますね！
ちなみに，$C$の固有値と$X$の特異値については $$ZZ^T = E^TXX^TE = E^TCEN = \Lambda N$$ $$ZZ^T = \Sigma V^T V \Sigma^T = \Sigma \Sigma^T$$ ゆえ， $$\Sigma \Sigma^T = \Lambda N$$ すなわち$C$の$i$番目の固有値を$\lambda_i$, そして$X$の$i$番目の特異値を$\sigma_i$と書くことにすると， $$\sigma_i^2 = \lambda_iN$$ の関係が成り立ちます。
3. 二つのデータ行列$X$と$Y$の共分散行列$C_{XY} = XY^T/N$を特異値分解＝最大共分散分析¶
この使い方は，すごく大雑把にいうと主成分分析の2変数バージョンです。日本の大気海洋コミュニティで「SVD解析を行った」という場合は，大抵この最大共分散分析（Maximum Covariance Analysis; MCA）を指します。
主成分分析は，データ行列$X$の分散を最大にするような基底を見つける分析方法でした。
それに対し，最大共分散分析は，（サンプリング次元の長さが等しい）データ行列$X$と$Y$の共分散を最大にするような基底を見つける分析方法です。たとえば，海面水温と気圧のSVD1は，そのSVD1に射影して求まる時系列の共分散が最大となるような空間パターンです。
主成分分析は「1変数の卓越した変動」を探す解析でした。特異値分解は「2変数の関係が深い変動」を探す解析，といってもいいでしょう。
...と言われても全くピンとこないと思うので，これは次回に回します。
課題7（締切：次回の授業開始前まで）¶
この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。
A問題. 適当な行列$X$を自分で定めて，np.linalg.svdを用いて特異値分解を行ってください。また上記と同様に，求めた特異値分解から元の行列を復元できることを確かめてください。
B問題. A問題で求めた$X$について，
・左特異ベクトルを並べた行列$U$が，$XX^T$の固有ベクトルを並べた行列と（固有ベクトルの符号をうまくとれば）等しくなること
・右特異ベクトルを並べた行列$V$が，$X^TX$の固有ベクトルを並べた行列と（固有ベクトルの符号をうまくとれば）等しくなること
・特異値が$XX^T$の非ゼロ固有値の平方根と等しいこと
・特異値が$X^TX$の非ゼロ固有値の平方根と等しいこと
の4つを確かめてください。
ヒント：本来ゼロになるべきところで，数値誤差で微小な負の値になってしまった場合，平方根を取れなくてエラーが出ます。0に置き換えて平方根を取ってみてください。
C問題. EOFに射影して求めたPC2と，右特異ベクトルとして求めたPC2は，正規化して符号を合わせると一致することを確かめてください。
D問題. 次の行列$X$を手計算で特異値分解してください。それが出来たら，手計算またはnumpyで$U\Sigma V^T$の行列計算を行って正しく分解できていることを確かめてください（うまくいかない場合は，一部の固有ベクトルの符号をひっくり返してみましょう）。次に，np.linalg.svdを用いた特異値分解の結果と比べて，固有ベクトルの符号を適切にひっくり返せば，自分の手計算の結果と一致することを確かめてください。
ヒント：$\Sigma$は一つに定まりますが，実は$U$と$V$の組は何通りもあり得ます。たとえば，$\vec{v}$が単位固有ベクトルならば，その逆ベクトル$-\vec{v}$も単位固有ベクトルですので，それだけでも色々な選び方がありうるのがわかります。
提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論7，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。