はじめに
機械学習ではPCAやLatent Semantic Analysis（LSA）、隠れマルコフモデルの学習、行列パラメータを持つモデルのトレースノルム正則化で、SVD（特異値分解）が必要になるときがあります。一回だけSVDすれば良い場合でも行列が巨大になるとつらいですし、行列パラメータを持つモデルでパラメータ更新のたびにSVDするときには行列のサイズがそこそこでもつらくなってきます。そこで特異値分解の高速解法について調べてみました。
元ネタはブログ記事[1]とそこで参照されている論文[2]です。特異値分解の高速化のための前処理についてのもので、もう8年も前の記事ですが、当時読んで「すげー！」と感動していたのを今頃必要に駆られて試してみたのでした。
特異値分解の復習
特異値分解
$M \times N$行列$\boldsymbol{A}$に対して$\boldsymbol{u} \in \mathbb{R}^M$、$\boldsymbol{v} \in \mathbb{R}^N$と実数$\sigma$が
\begin{align}
\boldsymbol{A} \boldsymbol{v} = \sigma \boldsymbol{u}
\end{align}
かつ
\begin{align}
\boldsymbol{A}^T \boldsymbol{u} = \sigma \boldsymbol{v}
\end{align}
を満たすとき、$\sigma$を特異値、$\boldsymbol{u}$を左特異ベクトル、$\boldsymbol{v}$を右特異ベクトルと呼びます。$\boldsymbol{u}$と$\boldsymbol{v}$は正規化されたものを考えます。また、異なる特異値に属する左特異ベクトル同士、右特異ベクトル同士は直交するように選びます。
多くても$\min(M,N)$組の特異値・左特異ベクトル・右特異ベクトルの組を持ちます。
$M > N$で$N$個の特異値があるとき
\begin{align}
\boldsymbol{U} &= (\boldsymbol{u}_1, \cdots, \boldsymbol{u}_N) \\
\boldsymbol{V} &= (\boldsymbol{v}_1, \cdots, \boldsymbol{v}_N) \\
\boldsymbol{\sigma} &= (\sigma_1, \cdots, \sigma_N)^T
\end{align}
と置くと、$\boldsymbol{A} \boldsymbol{v} = \sigma \boldsymbol{u}$の分解を並べると
\begin{align}
&\boldsymbol{A} \boldsymbol{V} = \boldsymbol{U} {\rm diag}(\boldsymbol{\sigma}) \\
&\Leftrightarrow
\boldsymbol{A} \boldsymbol{V} \boldsymbol{V}^T
= \boldsymbol{U} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{V}^T \\
&\Leftrightarrow
\boldsymbol{A}
= \boldsymbol{U} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{V}^T \\
\end{align}
となり、これが特異値分解（Singular Value Decomposition：SVD）です。右特異値ベクトル同士は直交することから$\boldsymbol{V} \boldsymbol{V}^T = \boldsymbol{I}$を使っています。
特異値分解の求め方
基本的な考え方として、$\boldsymbol{A} \boldsymbol{A}^T$は
\begin{align}
\boldsymbol{A} \boldsymbol{A}^T
&= \boldsymbol{U} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{V}^T \boldsymbol{V} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{U}^T \\
&= \boldsymbol{U} {\rm diag}(\sigma_1^2, \cdots, \sigma_D^2) \boldsymbol{U}^T \\
\end{align}
なので、$\boldsymbol{U}$と$\boldsymbol{\sigma}$は$\boldsymbol{A} \boldsymbol{A}^T$の固有値分解により求めることができます。あるいは$\boldsymbol{V}$を$\boldsymbol{A}^T \boldsymbol{A}$の固有値分解により求めることもできます。$\boldsymbol{u}_i$と$\boldsymbol{v}_i$は$\boldsymbol{A} \boldsymbol{v} = \sigma \boldsymbol{u}$の関係にあるので、一方が求まれば他方が求まります。
上位k番目までのSVD（k-SVD）を求める際にはLazySVD [3]を用いるのが良いようです。この手法は1-SVDをk回繰り返すものです。
今回は特異値分解そのものでなく、高速化のための前処理について説明します。
特異値分解の高速化のための前処理
概要
[2]の手法を大雑把に言うと、$M \times N$行列（$M>N$とします）$\boldsymbol{A}$のランクが$R \ll N$のとき、サイズ$R \times N$に圧縮した行列をSVDし、この結果を元の行列$\boldsymbol{A}$のSVDに変換する、というものです。
ランク$R$が$N$より小さければ小さいほどこの高速化の恩恵が得られます が、$R \approx N$のときにはメリットはありません。しかしながら、実データや行列パラメータをトレースノルム正則化する場合などでは、$R \ll N$が想定されるので、大幅な高速化が期待できます。
ここで $R$は後段のk-SVDで取り出したい$k$とは異なる ことに注意が必要です。実データではランク$R$は不明ですが、[2]では許容値$\epsilon$で打ち切られるところまでインクリメンタルに基底を計算するアルゴリズムが紹介されています。
[1]では取り出したい$k$を指定して、高速にSVDするライブラリredsvdが公開されています。内部的に$R=k$として、$k$本までの基底で圧縮行列を構成しているようなので、redsvdを使う場合は、本当に取り出したい$k$ではなく、 大きめのランク$R$を指定してredsvdしてから$k$番目まで取り出したほうが良さそうです。このように求めたredsvdの最小特異値（$R$番目の特異値）が、取り出したい$k$番目の特異値に比較して十分小さければ問題ないと思われます。
手法
[2]のおおまかな流れとしては、
- $\boldsymbol{A}$の列に対する$R$本の正規直交基底$\{\boldsymbol{q}_r \}_{r=1}^R$を列に並べた行列$\boldsymbol{Q} = (\boldsymbol{q}_1,\cdots,\boldsymbol{q}_R) \in \mathbb{R}^{M \times R}$を求める。
- $\boldsymbol{A} = \boldsymbol{Q} \boldsymbol{Q}^T \boldsymbol{A}$なので、$\boldsymbol{Q}^T \boldsymbol{A} \in \mathbb{R}^{R \times N}$について特異値分解：$\boldsymbol{Q}^T \boldsymbol{A} = \boldsymbol{\tilde{U}} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{V}$を求める。
- $\boldsymbol{A} = \boldsymbol{Q} \boldsymbol{\tilde{U}} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{V}$なので、$\boldsymbol{U} = \boldsymbol{Q} \boldsymbol{\tilde{U}}$とすれば$\boldsymbol{A}$の特異値分解：$\boldsymbol{A} = \boldsymbol{U} {\rm diag}(\boldsymbol{\sigma}) \boldsymbol{V}$が求まる。
となります。$\boldsymbol{A} = \boldsymbol{Q} \boldsymbol{Q}^T \boldsymbol{A}$について補足すると、$\{\boldsymbol{q}_r \}_{r=1}^R$が$\boldsymbol{A}$の列ベクトルの正規直行基底なことから、$\boldsymbol{A}$の列ベクトルは$\boldsymbol{a}_n = \sum_{r=1}^R (\boldsymbol{q}_r^T \boldsymbol{a}_n) \boldsymbol{q}_r$と書けて、これを行列表現したものです。
前処理としてやることは$\boldsymbol{Q}$を求めるだけなので、問題は正規直交基底$\{\boldsymbol{q}_r\}_{r=1}^R$をどうやって効率的に求めるか、です。
これが[2]の論文のすごいところなのですが、ものすごく簡単です。
乱数ベクトル$\{\boldsymbol{w}_r\}_{r=1}^R$を生成したら、$\boldsymbol{y}_r = \boldsymbol{A} \boldsymbol{w}_r$とすると${\boldsymbol{y}_r}_{r=1}^R$は高確率で線形独立になるので、これを正規直交化すれば$\{\boldsymbol{q}_r\}_{r=1}^R$が得られる、というものです。
ランク$R$が高次元になるほど、ランク$R$と同数の$R$本の乱数ベクトルから基底ベクトルを求めるのでは誤差が生じそうですが、$l$本オーバーサンプリングすると指数的に高確率（$1-6 l^{-l}$（$l \geq 2$））で決定的手法と同等のError boundを達成できるようです。$l$は$5$とか$10$で十分とのことです。
ランク$R$でなく誤差$\parallel \boldsymbol{A} - \boldsymbol{Q} \boldsymbol{Q}^T \boldsymbol{A} \parallel$に対する許容値$\epsilon$を指定して、許容値$\epsilon$で打ち切られるまで、基底をインクリメンタルに生成するアルゴリズムも紹介されています。[1]には、オーバーサンプリングありでのアルゴリズムが記載されているのですが、ここではわかりやすさのためオーバーサンプリングなしで書きます。
\begin{align}
&r = 0 {,\;} \boldsymbol{Q} = () \\
&{\rm while \; not \; broken \; do} \\
&\quad r = r+1 \\
&\quad {\rm Draw \; standard \; Gaussian \; vector: \;} \boldsymbol{w}_r \\
&\quad \boldsymbol{q}_r = \boldsymbol{A} \boldsymbol{w}_r \\
&\quad {\rm for \;} s=1 {\rm \; to \;} r-1 {\rm \; do}\\
&\quad \quad \boldsymbol{q}_r = \boldsymbol{q}_r - (\boldsymbol{q}_s^T \boldsymbol{q}_r) \boldsymbol{q}_s\\
&\quad {\rm end \; for} \\
&\quad {\rm if\;} \parallel \boldsymbol{q}_r \parallel < \epsilon {\rm \; then \; break} \\
&\quad \boldsymbol{q}_r = \frac{\boldsymbol{q}_r}{\parallel \boldsymbol{q}_r \parallel} \\
&\quad \boldsymbol{Q} = (\boldsymbol{Q}, \boldsymbol{q}_r) \\
&{\rm end \; while}
\end{align}
真ん中あたりのfor文はグラム・シュミットの正規直交化法です。
たったこれだけで、許容値$\epsilon$の下でのランク$R$が求まり、$R \ll N$であればSVDの大幅な高速化ができます。
参考文献
[1] DO++ "行列分解ライブラリredsvdを公開しました"
[2] HALKO, Nathan; MARTINSSON, Per-Gunnar; TROPP, Joel A. Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions. 2009.
[3] ALLEN-ZHU, Zeyuan; LI, Yuanzhi. LazySVD: Even faster SVD decomposition yet without agonizing pain. In: Advances in Neural Information Processing Systems. 2016. p. 974-982.