はじめに

本稿はこちらの参考書の1.8節、Randomized Singular Value Decompositionの個人的備忘として記録。前回記事までは、特異値分解と固有値分解について扱ってきたが、このランダム化特異値分解は3章のスパース性と圧縮や7章のDMDと関連が強い項目なので、強い気持ちで精読するために独立させた記事とした。

概要

最近の数値計算やデータサイエンスでは巨大データ行列を正確かつ効率的に分解することがひとつコーナーストーンとなっている。多くの場合、対象の行列が内部に持つ低ランク構造を抽出することに焦点を当てているが、近年ではランダムサンプリングに基づく手法を用いることで極めて効率的な行列分解が出来ることが示されている。

この所謂、ランダム数値法と呼ばれる手法は、決定論的手法の$1/n$の計算量で正確な行列分解が可能であり、数値計算分野における線形代数を一変させる可能性を秘めている。さらに、測定値が増え、測定空間の次元が膨大になっても、データの本質的なランクはそれほど上がらないことがよくある。

故に、ランダム化手法による計算コストの節約は、今後数年〜数十年の間にデータの増加とともにより重要となってくる。

ランダム化線形代数

ランダム化線形代数はSingular Value Decompotision; SVDよりも一般的な概念である。ランダム化SVDに加えて、主成分分析、ピボット化LU分解、ピボット化QR分解、およびDMDのためのランダム化アルゴリズムが開発されている。
ここでは、$n>m$のtall-skinnyな行列を扱うが、short-fatな行列にも容易に一般化できる。

Step 0: 階数$r<m$を決める
Step 1: ランダム射影$\boldsymbol{P}$を用いて列空間をサンプリングし、その列が$\boldsymbol{X}$の列空間に近似する行列$\boldsymbol{Q}$を求める。即ち、$\boldsymbol{X}\simeq\boldsymbol{QQ}^*\boldsymbol{X}$
Step 2: 行列$\boldsymbol{X}$を$\boldsymbol{Q}$の部分空間へ射影する。つまり、$\boldsymbol{Y}=\boldsymbol{Q}^*\boldsymbol{X}$
そして、これで得た行列$\boldsymbol{Y}$を分解する。
Step 3: 行列$\boldsymbol{Q}$と$\boldsymbol{Y}$から得たモードを用いて、元の高次元モード$\boldsymbol{U}=\boldsymbol{QU}_{\boldsymbol{Y}}$を再構築

ランダム化特異値分解アルゴリズム

ここ20年ほどで、低ランクSVDを計算するためのランダム化アルゴリズムがいくつか提案されてきた。これらの手法は、構造化サンプリング行列を採用し、高速な行列演算を実現することで改良された。
ここでは、Halko, Martinsson, Troppらの手法を紹介する。

Step 1

行列$\boldsymbol{X}$の列空間をサンプリングするランダム射影$\boldsymbol{P}\in\mathbb{R}^{m\times r}$を構築する。

\begin{align}
\boldsymbol{Z}=\boldsymbol{XP}
\end{align}

階数$r<m$なので、行列$\boldsymbol{Z}$は$\boldsymbol{X}$よりもサイズは小さくなる。

ここで、ランダムな射影行列$\boldsymbol{P}$を用いるので、行列$\boldsymbol{X}$の重要な成分を射影することは難しいと思えるが、実際には高い確率で$\boldsymbol{X}$の列空間をうまく近似することができる。これは、ランダム化技法の持つ特性であり、適切に選ばれたランダム行列を用いれば、元の行列の情報をうまく保持したまま次元削減することができる。

故に、$\boldsymbol{Z}$の低ランクQR分解を行うことで、$\boldsymbol{X}$の直交規定を得ることが出来る。直交行列$\boldsymbol{Q}\in\mathbb{R}^{n\times r}$と上三角行列$\boldsymbol{R}\in\mathbb{R}^{r\times r}$を用いて、

\begin{align}
\boldsymbol{Z}=\boldsymbol{QR}
\end{align}

Step 2

低ランク基底$\boldsymbol{Q}$を用いて、$\boldsymbol{X}$をより小さな空間に射影することができる。

\begin{align}
\boldsymbol{Y}&=\boldsymbol{Q}^*\boldsymbol{X}
\end{align}

また、特異値$\sigma_k$が$k>r$で急激に小さくなる場合、$\boldsymbol{X}＝\boldsymbol{QY}$が成立し、より良い結果が得られる。
さて、ここで$\boldsymbol{Y}$に対してSVDを行うと、

\begin{align}
\boldsymbol{Y}=\boldsymbol{U}_{\boldsymbol{Y}}\boldsymbol{\Sigma V}^*
\end{align}

行列$\boldsymbol{Q}$は直交行列なので、$\boldsymbol{X}$と$\boldsymbol{Y}$のSVDした時の特異値行列$\boldsymbol{\Sigma}$とユニタリ行列$\boldsymbol{V}$は保存される。

直交変換の特異値保存についてはこちら

直交行列$\boldsymbol{Q}\in\mathbb{R}^{n\times r}$による直交変換のSVD左辺を展開する。

\begin{align}
\boldsymbol{Q}\left(\boldsymbol{U\Sigma V}^*\right)
&=\left(\boldsymbol{QU}\right)\boldsymbol{\Sigma V}^*
\end{align}

ここで、$\boldsymbol{Q}$は直交行列なので、それとユニタリ行列との積$\boldsymbol{QU}$もまた直交行列となる。即ち、$\boldsymbol{U}_{\boldsymbol{Y}}=\boldsymbol{QU}$であり、特異値行列と右ユニタリ行列は保存される。

Step 3

先ほどの結果から元の左ユニタリ行列$\boldsymbol{U}$は以下のように構成される。

\begin{align}
\boldsymbol{U}_{\boldsymbol{Y}}=\boldsymbol{QU}
\end{align}

以降の節からは、rSVDの特性と課題について述べる。

オーバーサンプリング

実は先述では嘘を吐いていたらしく、第一の問題として多くの行列$\boldsymbol{X}$は、正確な低ランク構造を持っていないことが一般的らしい。これは、指定したランク$r$よりも高いモード$\sigma_{r+1},\sigma_{r+2},\cdots,\sigma_n$に対しても非ゼロ特異値が存在することを意味する。
そのため、スケッチ行列$\boldsymbol{Z}$は厳密には$\boldsymbol{X}$の列空間を完全に網羅することは出来ない。一般的に、ランダムな射影行列$\boldsymbol{P}$の列数を$r$から$r+p$に増やすと、結果が大幅に改善される（$p$は通常、約5～10）。このオーバーサンプリングにより、スケッチ行列$\boldsymbol{Z}$の特異値スペクトル${\sigma_i}_{i=1}^n$の分散が減少し、より精度の高い結果が得られるようになる。

べき乗法

次の問題として、特異値スペクトル$\lbrace\sigma_i\rbrace_{i=1}^n$の減衰が遅く、切り捨てた特異値にデータ$\boldsymbol{X}$の有意な分散が含まれる場合である。こういう場合は$\boldsymbol{X}$を$q$乗反復処理することで、特異値減衰がより速い新しい行列$\boldsymbol{X}^{(q)}$を作ることができる。

\begin{align}
\boldsymbol{X}^{(q)}
&=\left(\boldsymbol{XX}^*\right)^q\boldsymbol{X} \\
&=\lbrace\left(\boldsymbol{U\Sigma V}^*\right)\left(\boldsymbol{U\Sigma V}^*\right)^*\rbrace^q\boldsymbol{X} \\
&=\left(\boldsymbol{U\Sigma}^2\boldsymbol{U}^*\right)^q\boldsymbol{X} \\
&=\left(\boldsymbol{U\Sigma}^{2q}\boldsymbol{U}^*\right)\left(\boldsymbol{U\Sigma V}^*\right) \\
&=\boldsymbol{U\Sigma}^{2q+1}\boldsymbol{V}^*
\end{align}

このべき乗法はランダム化分解の精度を劇的に改善する。
しかし、データ行列$\boldsymbol{X}$を$q$回演算しなければならないため、計算コストは高くなる。

誤差境界の調整

rSVDの最も重要な特性の1つは、誤差境界の調整可能性である。これは特異値スペクトル${\sigma_i}_{i=1}^n$、ランク$r$、オーバーサンプリングパラメータ$p$、および冪乗反復数$q$の関数となる。
決定論的アルゴリズムで達成可能な最良の誤差境界は以下の通り。

\begin{align}
\|\boldsymbol{X}-\boldsymbol{QY}\|_2\geq\sigma_{r+1}\left(\boldsymbol{X}\right)
\end{align}

ここで、$|\cdot|$はスペクトルノルム（最大特異値）を表す。この境界は、階数$r$の最良の近似でも、元の行列$\boldsymbol{X}$の次の切り捨て特異値以上の誤差を持つことを示す。つまり、ランダム化手法では、誤差の期待値を束縛することが可能である。

\begin{align}
\mathbb{E}\left(\|\boldsymbol{X}-\boldsymbol{QY}\|_2\right)
\leq
\left(
    1+
    \sqrt{\frac{r}{p-1}}+
    \frac{e\sqrt{r+p}}{p}\,
    \sqrt{m-r}    
\right)
^
{1/(2p+1)}
\sigma_{k+1}
\left(
    \boldsymbol{X}
\right)
\end{align}

ランダム射影行列の選択

ランダム行列$\boldsymbol{P}$の適切な選択法はいくつか用意がある。

ガウスランダム射影（=$\boldsymbol{P}$の要素はガウスランダム変数）は、有利な数学的特性とスケッチ行列$\boldsymbol{Z}$で抽出される情報の豊富さから、よく用いられる手法である。特に、ガウスランダム行列$\boldsymbol{P}$が、Xの重要な情報を射影するような悪い挙動をする可能性は非常に低い。しかし、ガウシアン投影は、生成、保存、計算にコストがかかる。

一様ランダム行列もよく用いられるが、同様の制限があり、Rademacher行列のようにエントリが等確率で+1または-1になるような代替案も提案されている。

構造化ランダム射影行列は、計算コストを$\mathcal{O}(nm\log{r})$まで削減し、効率的なスケッチ行列$\boldsymbol{Z}$を提供することができる。

さらにもう1つの選択肢はスパース射影行列$\boldsymbol{P}$で、これは保存と計算を改善しますが、スケッチ行列$\boldsymbol{Z}$に含まれる情報が少なくなるというトレードオフが存在する。これは最速の手法であるが、$\boldsymbol{X}$の構造が列部分空間に明確に局在している場合、列サンプリングによって情報が失われる可能性があるため、注意して使用すべきである。

rSVDの例

ここで、高解像度の画像データを用いてランダム化SVDのデモを行う。この実装はあくまでrSVDのアルゴリズム説明用であり、速度、データ転送、精度などの指標は気にしない実装となっている。
1. 元データ、2. 階数$r=400$のSVD、3. rSVD、の3つのデータの比較を行う。まずはライブラリのインポートと元データの読み込みから。

# import libraries
from matplotlib.image import imread
import matplotlib.pyplot as plt
import numpy as np
plt.rcParams['figure.figsize'] = [16, 6]
plt.rcParams.update({'font.size': 18})

# read data & check to display
X = imread('../../data/banta.jpg')
X = np.mean(X, -1)
plt.set_cmap('gray')
plt.imshow(X)
plt.axis('off')
plt.show()

出力結果はこちら。

画像サイズは$4032\times3024$である。RGBの平均を取ってグレースケールとしていることに留意。
続いて、通常のSVDを実施する関数を実装

# normal SVD reconstruction
def SVD(X):
    return np.linalg.svd(X, full_matrices=False)

同様にして、rSVDを実装

# randomized SVD reconstruction
def rSVD(X, r=400, q=1, p=5):
    # Step 1: Sample column space of X with P matrix
    ny = X.shape[1]
    P = np.random.randn(ny, r+p)
    Z = X @ P
    for k in range(q):
        Z = X @ (X.T @ Z)
    Q, R = np.linalg.qr(Z, mode='reduced')

    # Step 2: Compute SVD on projected Y = Q.T @ X
    Y = Q.T @ X
    Uy, S, Vh = np.linalg.svd(Y, full_matrices=0)
    U = Q @ Uy
    return U, S, Vh

これらの関数から得られた元画像行列の分解行列から指定の階数で画像を再構築して、エラー率を計算する関数を書き、実際に出力結果を比較してみる。

# calculate a reconstruct image and an error ratio
def reconstruct(X, U, S, Vh, r=400):
    rX = U[:, :(r+1)] @ np.diag(S[:(r+1)]) @ Vh[:(r+1), :]
    error = np.linalg.norm(X-rX, ord=2) / np.linalg.norm(X, ord=2)
    return rX, error

# normal SVD reconstruction
U, S, Vh = SVD(X)
X_SVD, err_SVD = reconstruct(X, U, S, Vh)

# randomized SVD reconstruction
rU, rS, rVh = rSVD(X)
X_rSVD, err_rSVD = reconstruct(X, rU, rS, rVh)

# Plot
fig, axs = plt.subplots(1, 2, figsize=(18, 15))
axs[0].imshow(X_SVD, aspect='equal')
axs[0].set_title(f'Normal SVD Reconstruction\nError: {err_SVD*100:.4f}%')
axs[0].axis('off')

axs[1].imshow(X_rSVD, aspect='equal')
axs[1].set_title(f'Randomized SVD Reconstruction\nError: {err_rSVD*100:.4f}%')
axs[1].axis('off')

plt.show()

出力結果は以下の通り。

確かに、ランダムサンプリングによる元画像の列空間の再現性が肌感よりも良いことが分かる。

次に、べき乗パラメータ$p$を変化させた時のエラー率の推移を観察してみる。ガウス分布から生成される適当な行列$\boldsymbol{X}$を新しく作成して、その特異値を$1$から$0$まで線形に減少するように書き換えて、再度$X$を構成し、それを用いてべき乗法の効果を調べる。

# illustrate power iterations
X = np.random.randn(1000,100)
U, S, Vh = SVD(X)
S = np.arange(1,0,-0.01)
S = np.arange(1,0,-0.01)
X = U @ np.diag(S) @ Vh

# Define color map
color_list = np.array([[0, 0, 2/3],
                       [0, 0, 1],
                       [0, 1/3, 1],
                       [0, 2/3, 1],
                       [0, 1, 1],
                       [1/3, 1, 2/3],
                       [2/3, 1, 1/3],
                       [1, 1, 0],
                       [1, 2/3, 0],
                       [1, 1/3, 0],
                       [1, 0, 0],
                       [2/3, 0, 0]])

plt.plot(S, 'o-', color='k', linewidth=2, label='SVD')

Y = X
for q in range(1, 6):
    Y = X @ X.T @ Y
    Uq, Sq, VTq = SVD(Y)
    plt.plot(
        Sq, '-o', color=tuple(color_list[2*q+1]), linewidth=2, label='rSVD, q = '+str(q))

plt.legend()
plt.show()

出力結果は以下の通り。

このように、パラメータ$q$が大きくなるにつれて収束スピードが早くなっていることが分かった。ところが、実際の画像データだと、$q=1$でサチったりするので、むやみに大きくする必要はないと考えている。

おわりに

ランダム化抽出の特性が理解できた。これをベースにスパース性圧縮やDMDを理解していく。