応用数学 ch6メモ（主軸変換とその応用：主成分分析、特異値分解）#
6.1 主成分分析#
主軸変換#
ベクトル\(x\)を単位ベクトル\(u\)の延長上の線\(l\)に射影した長さは、それらの内積
に等しい
原点\(O\)を平均とする\(N\)点\(\{x_{\alpha} \}\)を\(l\)上に射影した長さの平均は常に零である
原点\(O\)を平均とする\(N\)点\(\{x_{\alpha} \}\)を\(l\)上に射影した長さの 二乗の平均 は
となる。このとき、
を分散共分散行列とよぶ。
分散\(S\)を最大化する方向\(u_1\)を 主方向 と呼ぶ。
共分散行列Vは半正値対称行列
\(N\)点\(\{x_{\alpha} \}\)の主方向はVの最大固有値\(\lambda_1\)に対する固有ベクトル\(u_1\)の方向であり、その方向の分散は\(\lambda_1\)である
主軸を座標系にとるとその座標系では共分散行列はどうなるか#
元の座標系の点 \(\boldsymbol{x}_\alpha\) は，新しい座標系では \(\boldsymbol{x}_\alpha^{\prime}=\boldsymbol{U}^{\top} \boldsymbol{x}_\alpha\) と書ける．その共分散行列は次のようになる。
すなわち、共分散行列\(V\)の固有値を対角要素とする対角行列になる。
このように主軸方向に新しい座標系をとることを 主軸変換 という。
6.2 画像の表現#
画像の展開#
\(m\times m\)画素の大きさの画像があるとし、\((i,j)\)画素の濃淡値が\(f_{ij}\)の画像を\(\boldsymbol{f}_{ij} = (f_{ij})\)とする。
画像の内積とノルム#
画像\(\boldsymbol{f}_{ij} = (f_{ij})\)、\(\boldsymbol{g}_{ij} = (g_{ij})\)の内積と画像\(\boldsymbol{f}_{ij} = (f_{ij})\)のノルムを次のように定義する
内積が0となる2つの画像は互いに直交するという
画像の基底#
\(m^2\) 個の画素を持つ任意の画像がある \(m^2\) 個の画像 \(\boldsymbol{e}_1, \boldsymbol{e}_2, \ldots, \boldsymbol{e}_{m^2}\) の線形結合でただ一通りに表せるとき，\(\left\{e_i\right\}\) を基底と呼び，各 \(e_i\) を基底画像と呼ぶ。 これらが互いに直交するとき，その基底を直交基底と呼び，さらに各画像の ノルムがすべて 1 のとき正規直交基底と呼ぶ．
画像の展開#
一部の基底画像によって
と近似することを 展開 とよぶ。
近似の尺度として二乗誤差\(\left\|\boldsymbol{f}-\sum_{i=1}^n c_i \boldsymbol{e}_i\right\|^2\)を用いると、各係数\(c_i\)は
この誤差関数を最小化する \(c_i\) を求める：
これは内積を用いて次のように展開できる：
ここで基底が正規直交なので：
よって\(\sum_{i=1}^n \sum_{j=1}^n c_i c_j\left(e_i, e_j\right)\)の項は \(\sum_{i=1}^n c_i^2\) となり、したがって誤差関数は
となる。これを各 \(c_k\) について偏微分してゼロとおくと
画像の基底#
自明な基底#
\(m \times m\) 画像の \((i, j)\) 画素のみが 1 で，残りの画素がすべて 0 である画像を \(\boldsymbol{e}_{i j}\) とすると，任意の画像 \(\boldsymbol{f}=\left(f_{i j}\right)\) は次のように \(\left\{\boldsymbol{e}_{i j}\right\}\) の線形結合で表せる．
この \(m^2\) 枚の画像 \(\left\{\boldsymbol{e}_{i j}\right\}, i, j=0, \ldots, m-1\) を自明な基底という．しかし，各 \(\boldsymbol{e}_{i j}\) は 1 点のみが光る画像であり，画像としての意味はない．
アダマール変換の基底#
次のような基底を考える．
すべての画素が 1 の画像を考える．
画像を \((m / 2) \times(m / 2)\) 画素の 4 枚の小画像に分割する
その内の 2 枚のすべて の画素を 1 ，残りの 2 枚のすべての画素を -1 とする。ただし，すでに作った画像と同じ，または 1 と -1 を反転すると同じになるものは除く
各々 の \((m / 2) \times(m / 2)\) 小画像をさらに 4 枚の \((m / 4) \times(m / 4)\) 小画像に分割
以下同様にする
このような基底による画像の展開を アダマール変換という
アダマール変換（Hadamard Transform）とは、 「入力ベクトルに対して、非常にシンプルな +1 と -1 のみを使った線形変換」
行列要素がすべて +1 か -1（ゼロはない）
直交行列（行同士・列同士が直交している）
フーリエ変換の簡易版ともいえる（周波数解析に近い）
from scipy.linalg import hadamard
print(hadamard(2))
[[ 1 1]
[ 1 -1]]
離散コサイン変換の基底#
次のように定義した基底 \(\left\{\boldsymbol{e}_{k l}\right\}, k=0,1, \ldots, m-1 ; l=0,1, \ldots, m-1\) ，は正規直交基底である。
右辺のカッコの中は画像 \(\boldsymbol{e}_{k l}\) の \((p, q)\) 画素の値を表す．ただし，次のように置 いた。
この基底による展開を，画像の離散コサイン変換と呼ぶ．これは，式（4．79） の離散コサイン変換を \(x, y\) の 2 方向に適用したものである。図 6.8 に \(8 \times 8\)画像の場合の基底を示す．
画像の固有空間#
アダマール変換や離散コサイン変換の基底画像は規則的な模様であり、画像としての意味はない
意味のある画像からの正規直交基底を作ることを考える
多数の顔写真 \(\left\{\boldsymbol{f}_\alpha\right\}, \alpha=1, \ldots, N\) があるとする。これらの平均をとれば、平均的な顔の画像
が得られる。
また個々の画像と平均顔との差分
をとれば、個々の画像の平均からの差異を意味する
もし個々の顔の差の代表的な特徴が見つかれば、\(f_\alpha^{\prime}\)はその画像の定数倍で近似できる。
例えば顔の縦横比が違いなら、額と顎の部分がプラスで、頬の両側がマイナスの画像を代表的な差画像\(g_1\)を作ることができる
そして\(f_\alpha^{\prime}\)と\(g_1\)の定数倍で最良に近似した画像との差画像をさらにべつの代表的な特徴\(g_2\)の定数倍で近似して・・・を繰り返していく。このとき\(g_1,g_2,g_3,\dots\)が正規直交系であるなら、画像\(f_\alpha^{\prime}\)を\(\boldsymbol{g}_1, \boldsymbol{g}_2, \boldsymbol{g}_3, \ldots\)に関して展開して次の表現が得られる
各画像 \(f_\alpha\) の平均画像からの差 \(f_\alpha^{\prime}\) をよく表す画像が \(g\) であるとすると， これは \(f_\alpha^{\prime}\) が
の形で近似できるという意味である。
近似誤差を二乗誤差 \(\| \boldsymbol{f}_\alpha^{\prime} - C_\alpha \boldsymbol{g}\|^2\) で表すとすると、この解は解析的に \(C_\alpha=\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)\) となる（しかし\(g\)が未知なのでこれを求める必要がある）。
また、\(\boldsymbol{g}\) に適当な定数を掛けて \(\|\boldsymbol{g}\|=1\) であるように定めてあると仮定する。
この2点を使うと二乗誤差は
と整理できる。
\(J\)を最小にするには、\(\sum_{\alpha=1}^N\left(f_\alpha^{\prime}, g\right)^2\)を最大化すればいい。
\(m \times m\) 画像 \(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\) を \(m^2\) 次元列べクトル に並べ換えたものをそれぞれ \(\tilde{\boldsymbol{f}}_\alpha^{\prime}, \tilde{\boldsymbol{g}}\) とすると次のようになる。
ここで
である。行列 \(\tilde{\boldsymbol{M}}\) の固有値を \(\lambda_1, \ldots, \lambda_{m^2}\) ，対応する単位固有ベクトルを \(\tilde{\boldsymbol{g}}_1, \ldots, \tilde{\boldsymbol{g}}_{m^2}\) と する．各ベクトル \(\tilde{\boldsymbol{g}}_\kappa\) の \(m^2\) 個の成分を \(m \times m\) 画素に配置し直したものを \(\boldsymbol{g}_\kappa\) とし， \(\boldsymbol{g}_1, \boldsymbol{g}_2, \ldots, \boldsymbol{g}_{m^2}\) を 固有画像 と呼ぶ。これらは対称行列の固有ベクトル から作られているので，正規直交基底となる．
\(\sum_{\alpha=1}^N\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)^2\)を最大化する単位ベクトル\(\tilde{\boldsymbol{g}}\)は\(\tilde{\boldsymbol{M}}\)の最大固有値に対する固有ベクトルである。したがって，\(\tilde{\boldsymbol{g}}\) は符号を除いて \(\tilde{\boldsymbol{g}}_1\) に一致する。
このような、画像における主成分分析は画像処理の領域では 固有空間法 と呼ばれる。
6.3 特異値分解#
固有画像を作るには\(\tilde{\boldsymbol{M}}\)の固有値と固有ベクトルを計算する必要があるが、\(m^2\times m^2\)と巨大な行列であるため計算量が非常に多くなる。
固有値が0の固有ベクトルは応用上必要ないため、\(0\)でない固有値と固有ベクトルのみ計算できれば良い。そこで効率よく計算する方法が特異値分解。
問題は，\(N\) 個の \(n\left(=m^2\right)\) 次元ベクトル \(\boldsymbol{p}_1, \ldots, \boldsymbol{p}_N\) があるとき，\(n \times n\) 行列
の固有値，固有ベクトルを計算することである．
ベクトル \(\boldsymbol{p}_1, \ldots, \boldsymbol{p}_N\) を列とする \(n \times N\) 行列を
とすると、\(\boldsymbol{M}=\boldsymbol{P} \boldsymbol{P}^{\top}\)となる（ベクトルの直積）
行列\(M\)の固有値・固有ベクトルは、
で定義される行列\(N\)の行列の固有値・固有ベクトルを計算する方法でも求められる
\(M\)と\(N\)のゼロでない固有値は一致する
\(M\)の固有値と固有ベクトルをそれぞれ\(\lambda, u\)とすると
左から\(P^T\)をかけると
\(N\)の固有値と固有ベクトルをそれぞれ\(\lambda, v\)とすると
左から\(P\)をかけると
\(M\)の固有値\(\lambda (> 0)\)に対する単位ベクトルを\(u\)、\(N\)の固有値\(\lambda(> 0)\)に対する単位ベクトルを\(v\)とすると
となる。
例えば\(M\)の固有値\(\lambda\)に対する固有ベクトル\(\boldsymbol{P v}\)のノルムの二乗は
であるため、\(\sqrt{\lambda}\)で割って単位ベクトルにする
行列 \(\boldsymbol{M}\left(=\boldsymbol{P P}^{\top}\right)\) の固有値を大きい順に並べたものを \(\lambda_1, \lambda_2, \ldots, \lambda_n\) とし，対応する \(n\) 次元単位固有ベクトルの正規直交系を \(\boldsymbol{u}_1, \boldsymbol{u}_2, \ldots, \boldsymbol{u}_n\) とする。こ れらを列として並べた \(n \times n\) 行列を
とする．同様に，行列 \(\boldsymbol{N}\left(=\boldsymbol{P}^{\top} \boldsymbol{P}\right)\) の固有値を大きい順に並べたものを \(\lambda_1, \lambda_2, \ldots, \lambda_N\) とし，対応する \(N\) 次元単位固有ベクトルの正規直交系を \(v_1, v_2, \ldots, v_N\) とする．これらを列として並べた \(N \times N\) 行列を
とする．行列 \(\boldsymbol{M}\) および \(\boldsymbol{N}\) の（共通の）ランク（=0でない固有値の個数） を \(r\) とすると，\(i>r\) に対して \(\lambda_i=0\) である。
Note
行列\(P\)は直交行列\(U,V\)によって次のように書ける
\(\boldsymbol{M}\) の単位固有ベクトル \(\boldsymbol{u}_1, \boldsymbol{u}_2, \ldots, \boldsymbol{u}_r\) の符号を上式の右辺の \(\pm\) がすべて + に なるように選べば
右から \(\boldsymbol{V}^{\top}\) を掛けると
これを行列\(P\)の 特異値分解 といい、\(\sigma_i\)を 特異値 という。
例#
import numpy as np
P = np.array([
[1, 2],
[3, 4],
[5, 6],
])
M = P @ P.T
print("M:\n", M)
eigenvalues, eigenvectors = np.linalg.eig(M)
print("固有値:", eigenvalues.round(2))
print("固有ベクトル:\n", eigenvectors.round(2))
M:
[[ 5 11 17]
[11 25 39]
[17 39 61]]
固有値: [90.74 0.26 0. ]
固有ベクトル:
[[-0.23 -0.88 0.41]
[-0.52 -0.24 -0.82]
[-0.82 0.4 0.41]]
N = P.T @ P
print("N:\n", N)
eigenvalues, eigenvectors = np.linalg.eig(N)
print("固有値:", eigenvalues.round(2))
print("固有ベクトル:\n", eigenvectors.round(2))
N:
[[35 44]
[44 56]]
固有値: [ 0.26 90.74]
固有ベクトル:
[[-0.78 -0.62]
[ 0.62 -0.78]]
u1 = P @ eigenvectors[0] / np.sqrt(eigenvalues[0])
u2 = P @ eigenvectors[1] / np.sqrt(eigenvalues[1])
u2
array([-0.09974884, -0.13444826, -0.16914768])
import numpy as np
def svd_from_eig(A):
AtA = A.T @ A
# 固有値分解（V: 右特異ベクトル）
eigenvalues, V = np.linalg.eig(AtA)
# 固有値の実数部分とインデックスを昇順で並び替え
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
V = V[:, idx]
# 特異値の計算（負値を0に修正）
singular_values = np.sqrt(np.clip(eigenvalues.real, 0, None))
# 左特異ベクトルの計算
U = []
for i in range(len(singular_values)):
if singular_values[i] > 1e-10:
u = A @ V[:, i] / singular_values[i]
else:
# ランク欠損時の直交補完（ここではゼロベクトル）
u = np.zeros(A.shape[0])
U.append(u)
U = np.column_stack(U)
# 必要ならUの直交補完を追加（ここでは省略）
# Σ 行列の構築
m, n = A.shape
Sigma = np.zeros((m, n))
for i in range(min(m, n)):
Sigma[i, i] = singular_values[i]
return U, Sigma, V
U, S, V = svd_from_eig(P)
# 検証
print("U:\n", U)
print("Σ:\n", S)
print("V:\n", V)
print("再構成:\n", U @ S @ V.T)
U:
[[-0.2298477 0.88346102]
[-0.52474482 0.24078249]
[-0.81964194 -0.40189603]]
Σ:
[[9.52551809 0. ]
[0. 0.51430058]
[0. 0. ]]
V:
[[-0.61962948 -0.78489445]
[-0.78489445 0.61962948]]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[31], line 45
43 print("Σ:\n", S)
44 print("V:\n", V)
---> 45 print("再構成:\n", U @ S @ V.T)
ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 3 is different from 2)