はじめに

情報工学科出身なので、特異値分解は大学2,3年生の頃に履修しているのだが、就職して2年ほど統計・機械学習分野から離れていたため、今回の復習を備忘録として記録。参考書籍はこちら。
本稿では特異値分解の基礎および実際の画像を使った復元を実装してみる。

概要

そもそも特異値分解(Singular Value Decomposition；SVD)とはなんだっけ、というところから理解を始める。
ずばりその目的は、高次元データをシステマチックに低次元へ近似することにある。
実際のデータは特徴量が多く、高次元ベクトルのデータ分析が求められることが常であるため、特異値分解はデータサイエンス領域において非常に重要な役割を担い、さまざまな技術の基礎となっている；動的モード分解、主成分分析、固有直交分解など。
また、次元圧縮以外にも劣決定性/優決定性問題や擬似逆行列などにおいても強力な応用を持つ。

導入

早速、数学的な導入を簡単に行う。
まずデータセット$\boldsymbol{X}\in\mathbb{C}^{n\times m}$を用意する。

\boldsymbol{X}=
\begin{bmatrix}
| & | & & | \\
\boldsymbol{x}_1 & \boldsymbol{x}_2 & \cdots & \boldsymbol{x}_m \\
| & | & & |
\end{bmatrix}

例えばここで、$n>>m$の時はtall-skinnyな行列と呼び、反対に$n<<m$の時short-fatな行列と呼ばれる。
特異値分解は任意の複素行列$\boldsymbol{X}$において、以下のような行列分解が一意に定まることを保証する。

\boldsymbol{X}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^*

ここで行列$\boldsymbol{U}\in\mathbb{C}^{n\times n}$と$\boldsymbol{V}\in\mathbb{C}^{m\times m}$はユニタリ行列であり、正規直交な列ベクトルを持つ。
行列$\boldsymbol{\Sigma}\in\mathbb{C}^{n\times m}$は非負対角矩形行列である。

計算方法

$n\geq m$の時、これはデータセット$\boldsymbol{X}$がサンプル数よりも特徴量数の方が多いことを示唆するが、行列$\boldsymbol{\Sigma}$の階数を考えると、

\boldsymbol{\Sigma}=
\begin{bmatrix}
\hat{\boldsymbol{\Sigma}}\\
\boldsymbol{0}
\end{bmatrix}

となり、2つの行列$\hat{\boldsymbol{\Sigma}}$と$\boldsymbol{0}$に分解できるため、計算量を削減できる。(what we call economy SVD)

\boldsymbol{X}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^*
=
\begin{bmatrix}
\hat{\boldsymbol{U}} & \hat{\boldsymbol{U}}^\perp
\end{bmatrix}
\begin{bmatrix}
\hat{\boldsymbol{\Sigma}}\\
\boldsymbol{0}
\end{bmatrix}
\boldsymbol{V}^*
=\hat{\boldsymbol{U}}\hat{\boldsymbol{\Sigma}}\boldsymbol{V}^*

因みに、Pythonのnumpyパッケージでは以下のようにfull_matrices引数をFalseにすることで、上記の計算手法を取ってくれる。

import numpy as np
X = np.random.rand(10, 7)
U, S, Vh = np.linalg.svd(X, full_matrices=False)

行列近似

前述の通り、SVDの最も使われる用途は低ランク近似である。行列$\boldsymbol{\Sigma}$が非負対角矩形行列であることから、

\boldsymbol{X}=\sum_{k=1}^{m}\sigma_k \boldsymbol{u}_k \boldsymbol{v}_k^*
=\sigma_1 \boldsymbol{u}_1 \boldsymbol{v}_1^*+\sigma_2 \boldsymbol{u}_2 \boldsymbol{v}_2^*+\cdots+\sigma_m \boldsymbol{u}_m \boldsymbol{v}_m^*

ここで、$\sigma_k$は行列$\boldsymbol{\Sigma}$の$k$番目の要素。$\boldsymbol{u}_k$と$\boldsymbol{v}_k$はそれぞれ行列$\boldsymbol{U}$と$\boldsymbol{V}$の$k$番目の列ベクトルを表す。
また、$\boldsymbol{u}_k\boldsymbol{v}_k^* $はテンソル積であることに留意。
さて、ここで$\sigma_k$を降順に並び替え、$\sigma_1\geq\sigma_2\geq\cdots\geq\sigma_k\geq 0$を保証してみると、行列$\sigma_k \boldsymbol{u}_k \boldsymbol{v}_k^*$が持つ情報はそれよりも前の項より小さいことになる。
つまり適当な階数$r$で切れば、行列$\boldsymbol{X}$の近似を得ることができる。

\boldsymbol{X}\simeq\tilde{\boldsymbol{X}}=\sum_{k=1}^r\sigma_k \boldsymbol{u}_k\boldsymbol{v}_k^*
=\sigma_1 \boldsymbol{u}_1 \boldsymbol{v}_1^*+\sigma_2 \boldsymbol{u}_2 \boldsymbol{v}_2^*+\cdots+\sigma_r \boldsymbol{u}_r \boldsymbol{v}_r^*
=\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^*

Eckart-Youngの定理

前節において行列の低ランク近似について説明したが、一旦それを置いておいて行列のノルム、つまり2つの行列がどのくらい違っているかを測るための尺度を導入する。
行列$\boldsymbol{A}\in\mathbb{C}^{n\times m}, n\geq m$を考えたとき、その特異値も$\sigma_k$と表現して

\|\boldsymbol{A}\|_F=\sqrt{\sum_i^m\sum_j^n\mid a_{ij}\mid^2}
=\sqrt{\mathrm{tr}(A^*A)}=\sqrt{\sum_{i=1}^m \sigma^2_i}

このノルムはフロベニウスノルムと呼ばれ、これを目的関数とした場合の最小化問題の解が前項で述べたSVD(what we call truncated SVD)である、という内容がEckart-Youngの定理である。

\newcommand{\argmax}{\mathop{\rm arg~max}\limits}
\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
\argmin_{\tilde{\boldsymbol{X}}, \mathrm{s.t. rank}(\tilde{\boldsymbol{X}})=r}\|\boldsymbol{X}-\tilde{\boldsymbol{X}}\|_F
=\tilde{\boldsymbol{U}}\tilde{\boldsymbol{\Sigma}}\tilde{\boldsymbol{V}}^*

解釈としてはこれが絶対誤差のようなものに該当するので、階数$r$のSVD近似率は以下のように表せられる。

\frac{\|\boldsymbol{X}-\tilde{\boldsymbol{X}}\|_F}{\|\boldsymbol{X}\|_F}

画像データの低ランク近似

上記を踏まえて実際に画像データの低ランク近似をPythonを使って行なってみる。
用いる画像はこちら。沖縄県読谷村に位置する星野リゾート内にあるバンタカフェで撮った写真である。

余談であるが、読谷村は村ながら人口約42,000であり、これは全国の村の中で首位である。(2023年7月時点)
まずは、ライブラリのインポートとファイルの読み込みを行う。

### import libraries
import imageio
from matplotlib.image import imread
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from tqdm import tqdm

### read data & check to display
fig = imread('../data/banta.jpg')
plt.imshow(fig)
plt.show()

留意点として、画像はBGRなので$3024\times 4032\times 3$のテンソルデータである。
したがってBGRそれぞれについてSVDを実施する。

### Take the Singular Value Decompose
U, S, Vh = list(), list(), list()
for c in range(3):
    u, s, vh = np.linalg.svd(fig[:,:,c], full_matrices=False)
    s = np.diag(s)
    U.append(u)
    S.append(s)
    Vh.append(vh)

次に階数$r$を1から200まで変化させて低ランク近似を行う。

### Approximate matrix with truncated SVD for various ranks r
image_list = list()
for r in tqdm(range(1, 200)):
    appx = list()
    for c in range(3):
        appx.append(np.dot(np.dot(U[c][:,:r], S[c][:r,:r]), Vh[c][:r,:]))
    appx = np.stack(appx, axis=-1)
    appx = Image.fromarray(np.uint8(appx))
    image_list.append(appx)

出力された画像を束ねて.gifとして保存する。

### create the animation (.gif)
with imageio.get_writer('appx.gif', duration=50) as fd:
    for image in image_list:
        fd.append_data(image)

.gifデータは重たいので、階数$r=1,10,100$の場合だけをこちらに掲載する。
$r=1$、近似率=29.6%

$r=10$、近似率=8.7%

$r=100$、近似率=2.3%

近似率の計算はフロベニスノルムを3階テンソルに拡張したものとして、以下のように実装した。
即席の指標としてなので、コードの計算速度も遅く、そもそも安易にフロベニウスノルムをテンソルに拡張して良いのかも分からない。(詳しい方是非教えて下さい。)

def frobenius_norm(X):
    return np.sqrt(np.sum(np.power(X.flatten().astype(np.int32), 2)))

ご覧の通り、太陽光が反射している机の部分(白飛び)や陰になって黒くなっている部分は情報として価値が低い部分なため、階数を大きくしてもノイズが残っている状態である。
それと同様に$r=1$の画像を見てみると、ほぼ垂直方向に伸びているストローや背景の遮光カーテンや木々なども横方向に伸びており、これも情報が少ないと解釈されていることが見て取れる。

とは言え、元が$3024\times 4032$のデータであったにも関わらず$r=100$で十分に復元出来ていると言っても過言ではない。

まとめ

今回はここまで。
次回は主成分分析への応用について綴る予定。