応用数学 ch6メモ（主軸変換とその応用：主成分分析、特異値分解）

$m^2$ 個の画素を持つ任意の画像がある $m^2$ 個の画像 $\boldsymbol{e}_1, \boldsymbol{e}_2, \ldots, \boldsymbol{e}_{m^2}$ の線形結合でただ一通りに表せるとき， $\left\{e_i\right\}$ を基底と呼び，各 $e_i$ を基底画像と呼ぶ。これらが互いに直交するとき，その基底を直交基底と呼び，さらに各画像のノルムがすべて 1 のとき正規直交基底と呼ぶ．

画像の展開¶

一部の基底画像によって

\boldsymbol{f} \approx c_1 \boldsymbol{e}_1+c_2 \boldsymbol{e}_2+\cdots+c_n \boldsymbol{e}_n, \quad n<m^2

(8)

と近似することを展開とよぶ。

近似の尺度として二乗誤差 $\left\|\boldsymbol{f}-\sum_{i=1}^n c_i \boldsymbol{e}_i\right\|^2$ を用いると、各係数 $c_i$ は

c_i=\left(\boldsymbol{f}, \boldsymbol{e}_i\right), \quad i=1, \ldots, n

(9)

証明

この誤差関数を最小化する $c_i$ を求める：

E\left(c_1, \ldots, c_n\right)=\left\|f-\sum_{i=1}^n c_i e_i\right\|^2

(10)

これは内積を用いて次のように展開できる：

\begin{gathered} E=\left(f-\sum_{i=1}^n c_i e_i, f-\sum_{j=1}^n c_j e_j\right) \\ =(f, f)-2 \sum_{i=1}^n c_i\left(f, e_i\right) + \sum_{i=1}^n \sum_{j=1}^n c_i c_j\left(e_i, e_j\right) \end{gathered}

(11)

ここで基底が正規直交なので：

\left(e_i, e_j\right)= \begin{cases}1 & i=j \\ 0 & i \neq j\end{cases}

(12)

よって $\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$ となり、したがって誤差関数は

E=\|f\|^2-2 \sum_{i=1}^n c_i\left(f, e_i\right)+\sum_{i=1}^n c_i^2

(13)

となる。これを各 $c_k$ について偏微分してゼロとおくと

\frac{\partial E}{\partial c_k}=-2\left(f, e_k\right)+2 c_k=0 \quad \Rightarrow \quad c_k=\left(f, e_k\right)

(14)

画像の基底¶

自明な基底¶

$m \times m$ 画像の $(i, j)$ 画素のみが 1 で，残りの画素がすべて 0 である画像を $\boldsymbol{e}_{i j}$ とすると，任意の画像 $\boldsymbol{f}=\left(f_{i j}\right)$ は次のように $\left\{\boldsymbol{e}_{i j}\right\}$ の線形結合で表せる．

\boldsymbol{f}=\sum_{i, j=0}^{m-1} f_{i j} e_{i j}

(15)

この $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]]

import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import hadamard

# アダマール行列を作成（サイズは2のべき乗）
n = 4
H = hadamard(2**n)

# プロット
plt.figure(figsize=(3, 3))
plt.imshow(H, cmap='gray', interpolation='nearest')
plt.title(f"Hadamard Matrix (size {2**n}×{2**n})")
plt.show()

離散コサイン変換の基底¶

次のように定義した基底 $\left\{\boldsymbol{e}_{k l}\right\}, k=0,1, \ldots, m-1 ; l=0,1, \ldots, m-1$ ，は正規直交基底である。

\boldsymbol{e}_{k l}=\left(\frac{2 \alpha_k \alpha_l}{m} \cos \frac{\pi k(2 p+1)}{2 N} \cos \frac{\pi l(2 q+1)}{2 N}\right)

(16)

右辺のカッコの中は画像 $\boldsymbol{e}_{k l}$ の $(p, q)$ 画素の値を表す．ただし，次のように置いた。

\alpha_k=\left\{\begin{array}{cc} \frac{1}{\sqrt{2}} & k=0 \\ 1 & k \neq 0 \end{array}\right.

(17)

この基底による展開を，画像の離散コサイン変換と呼ぶ．これは，式（4．79）の離散コサイン変換を $x, y$ の 2 方向に適用したものである。図 6.8 に $8 \times 8$ 画像の場合の基底を示す．

画像の固有空間¶

アダマール変換や離散コサイン変換の基底画像は規則的な模様であり、画像としての意味はない

意味のある画像からの正規直交基底を作ることを考える

多数の顔写真 $\left\{\boldsymbol{f}_\alpha\right\}, \alpha=1, \ldots, N$ があるとする。これらの平均をとれば、平均的な顔の画像

g_0=\frac{1}{N} \sum_{\alpha=1}^N f_\alpha

(18)

が得られる。

また個々の画像と平均顔との差分

f_\alpha^{\prime}=f_\alpha-g_0

(19)

をとれば、個々の画像の平均からの差異を意味する

もし個々の顔の差の代表的な特徴が見つかれば、 $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$ に関して展開して次の表現が得られる

\boldsymbol{f}_\alpha=\boldsymbol{g}_0+c_{\alpha 1} \boldsymbol{g}_1+c_{\alpha 2} \boldsymbol{g}_2+c_{\alpha 3} \boldsymbol{g}_3+\cdots

(20)

各画像 $f_\alpha$ の平均画像からの差 $f_\alpha^{\prime}$ をよく表す画像が $g$ であるとすると，これは $f_\alpha^{\prime}$ が

\boldsymbol{f}_\alpha^{\prime} \approx C_\alpha \boldsymbol{g}

(21)

の形で近似できるという意味である。

近似誤差を二乗誤差 $\| \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点を使うと二乗誤差は

\begin{aligned} J & =\sum_{\alpha=1}^N\left\|\boldsymbol{f}_\alpha^{\prime}-C_\alpha \boldsymbol{g}\right\|^2 =\sum_{\alpha=1}^N\left(\boldsymbol{f}_\alpha^{\prime}-C_\alpha \boldsymbol{g}, \boldsymbol{f}_\alpha^{\prime}-C_\alpha \boldsymbol{g}\right) \\ & =\sum_{\alpha=1}^N\left(\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{f}_\alpha^{\prime}\right)-2 C_\alpha\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)+C_\alpha^2(\boldsymbol{g}, \boldsymbol{g})\right) \\ & =\sum_{\alpha=1}^N\left(\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{f}_\alpha^{\prime}\right)-2\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)^2+\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)^2(\boldsymbol{g}, \boldsymbol{g})\right) \\ & =\sum_{\alpha=1}^N\left(\left\|\boldsymbol{f}_\alpha^{\prime}\right\|^2-\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)^2\right) \end{aligned}

(22)

と整理できる。

$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}}$ とすると次のようになる。

\begin{aligned} \sum_{\alpha=1}^N\left(\boldsymbol{f}_\alpha^{\prime}, \boldsymbol{g}\right)^2 & =\sum_{\alpha=1}^N\left(\tilde{\boldsymbol{f}}_\alpha^{\prime}, \tilde{\boldsymbol{g}}\right)^2=\sum_{\alpha=1}^N \tilde{\boldsymbol{g}}^{\top} \tilde{\boldsymbol{f}}_\alpha^{\prime} \tilde{\boldsymbol{f}}_\alpha^{\prime \top} \tilde{\boldsymbol{g}} \\ & =\left(\tilde{\boldsymbol{g}},\left(\sum_{\alpha=1}^N \tilde{\boldsymbol{f}}_\alpha^{\prime} \tilde{\boldsymbol{f}}_\alpha^{\prime\top}\right) \tilde{\boldsymbol{g}}\right) =(\tilde{\boldsymbol{g}}, \tilde{\boldsymbol{M}} \tilde{\boldsymbol{g}}) \end{aligned}

(23)

ここで

\tilde{\boldsymbol{M}}:=\sum_{\alpha=1}^N \tilde{\boldsymbol{f}}_\alpha^{\prime} \tilde{\boldsymbol{f}}_\alpha^{\prime \top}

(24)

である。行列 $\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{M}=\sum_{\alpha=1}^N \boldsymbol{p}_\alpha \boldsymbol{p}_\alpha^{\top}

(25)

の固有値，固有ベクトルを計算することである．

ベクトル $\boldsymbol{p}_1, \ldots, \boldsymbol{p}_N$ を列とする $n \times N$ 行列を

\boldsymbol{P}=\left(\boldsymbol{p}_1 \quad \boldsymbol{p}_2\quad \cdots \quad \boldsymbol{p}_N\right)

(26)

とすると、 $\boldsymbol{M}=\boldsymbol{P} \boldsymbol{P}^{\top}$ となる（ベクトルの直積）

行列 $M$ の固有値・固有ベクトルは、

\boldsymbol{N}=\boldsymbol{P}^{\top} \boldsymbol{P}

(27)

で定義される行列 $N$ の行列の固有値・固有ベクトルを計算する方法でも求められる

$M$ と $N$ のゼロでない固有値は一致する

$M$ の固有値と固有ベクトルをそれぞれ $\lambda, u$ とすると

\underbrace{ PP^T }_{=M} u = \lambda u

(28)

左から $P^T$ をかけると

\underbrace{ P^T P}_{=N} P^T u = \lambda P^T u\\ \iff N (P^T u) = \lambda (P^T u)

(29)

$N$ の固有値と固有ベクトルをそれぞれ $\lambda, v$ とすると

\underbrace{ P^T P }_{=N} v = \lambda v

(30)

左から $P$ をかけると

\underbrace{ P P^T }_{=M} P v = \lambda P v\\ \iff M (P v) = \lambda (P v)

(31)

$M$ の固有値 $\lambda (> 0)$ に対する単位ベクトルを $u$ 、 $N$ の固有値 $\lambda(> 0)$ に対する単位ベクトルを $v$ とすると

\boldsymbol{u}=\frac{ \pm \boldsymbol{P} \boldsymbol{v}}{\sqrt{\lambda}}, \quad \boldsymbol{v}=\frac{ \pm \boldsymbol{P}^{\top} \boldsymbol{u}}{\sqrt{\lambda}}

(32)

となる。

例えば $M$ の固有値 $\lambda$ に対する固有ベクトル $\boldsymbol{P v}$ のノルムの二乗は

\left\|\boldsymbol{P v}\right\|^2=(\boldsymbol{P v}, \boldsymbol{P v})=\left(\boldsymbol{v}, \boldsymbol{P}^{\top} \boldsymbol{P v}\right)=(\boldsymbol{v}, \boldsymbol{N v})=(\boldsymbol{v}, \lambda \boldsymbol{v})=\lambda\|\boldsymbol{v}\|^2=\lambda

(33)

であるため、 $\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{U}=\left(\begin{array}{llll} \boldsymbol{u}_1 & \boldsymbol{u}_2 & \cdots & \boldsymbol{u}_n \end{array}\right)

(34)

とする．同様に，行列 $\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{V}=\left(\begin{array}{llll} \boldsymbol{v}_1 & \boldsymbol{v}_2 & \cdots & \boldsymbol{v}_N \end{array}\right)

(35)

とする．行列 $\boldsymbol{M}$ および $\boldsymbol{N}$ の（共通の）ランク（=0でない固有値の個数）を $r$ とすると， $i>r$ に対して $\lambda_i=0$ である。

\boldsymbol{P} \boldsymbol{V}=\left(\begin{array}{llll} \boldsymbol{P} \boldsymbol{v}_1 & \boldsymbol{P} \boldsymbol{v}_2 & \cdots & \boldsymbol{P} \boldsymbol{v}_N \end{array}\right)=\left(\begin{array}{lllll} \pm \sqrt{\lambda_1} \boldsymbol{u}_1 & \cdots & \pm \sqrt{\lambda_r} \boldsymbol{u}_r & \mathbf{0} & \cdots & \mathbf{0} \end{array}\right)

(37)

$\boldsymbol{M}$ の単位固有ベクトル $\boldsymbol{u}_1, \boldsymbol{u}_2, \ldots, \boldsymbol{u}_r$ の符号を上式の右辺の ± がすべて + になるように選べば

\boldsymbol{P} \boldsymbol{V}=\left(\boldsymbol{u}_1 \cdots \boldsymbol{u}_n\right)\left(\begin{array}{ccc} \sqrt{\lambda_1} & & \\ & \ddots & \\ & & \sqrt{\lambda_r} \end{array}\right)=\boldsymbol{U}\left(\begin{array}{ccc} \sigma_1 & & \\ & \ddots & \\ & & \\ & & \sigma_r \end{array}\right)

(38)

右から $\boldsymbol{V}^{\top}$ を掛けると

\boldsymbol{P}=\boldsymbol{U}\left(\begin{array}{ccc} \sigma_1 & & \\ & \ddots & \\ & & \sigma_r \end{array}\right) \boldsymbol{V}^{\top}

(39)

これを行列 $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)