はじめに

この記事では、Random Projectionを用いた行列の低ランク近似を逐次的に処理する方法についてまとめる。
参考文献は以下の通りである。なお、Qiitaにはすでに大変参考になる記事があるので、ご参照いただきたい（前処理による特異値分解の高速化）。

Erichson et al.(2018), Randomized Nonnegative Matrix Factorization
Halko et al.(2011), Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions

基本自分への備忘録なので、適当なこと抜かしてたらすいません。

低ランク近似

機械学習などの文脈で現れるデータ行列は大規模なものが多く、しばしば計算コストやメモリ消費量が問題となる。一方、こうした行列は実際のサイズに比して低ランクである場合が多く、元の行列の低ランク近似を求めることで大幅に計算コストを削減できる。こうした処理は次元削減として知られる。

よく知られた低ランク近似の手法は特異値分解であり、上位の少数の支配的な特異値に対応する部分だけを残すことで元の行列の低ランク近似を求めることができる（Truncated SVD）。しかし、特異値分解の処理はかなり重く、行列のサイズが大きくなると特異値分解を計算することそのものが難しくなる。また、大規模なデータ行列はその全体をメモリに格納できないため、メモリ消費を抑えつつ処理を行う必要がある。

ここでは、乱数を使った低ランク近似の手法について説明する。特に、行列全体をメモリに格納できない場合のアルゴリズムについて述べる。逐次処理によって特異値分解を少ないメモリ量で行うアルゴリズム（Incremental SVD）もあるが、ここでは扱わない。

Random Projection

データ行列を$\mathbf{X} \in \mathbb{R}^{n \times m}$とし、$r \triangleq \operatorname{rank} \mathbf{X} \ll m \ll n$が成り立つとする。すなわち、$\mathbf{X}$は縦に細長い巨大な密行列であり、かつそのランクは列数よりもはるかに小さいとする。

このとき、以下のような近似が成り立つ。

\mathbf{X} \approx \mathbf{Q} \mathbf{Q}^{\top}\mathbf{X}

ここで、$\mathbf{Q} \in \mathbb{R}^{n \times r}$は$\mathbf{X}$の列空間の直交基底を$r$本並べた直交行列であり、上式は$r$次元部分空間への直交射影を表す。

\mathbf{B} = \mathbf{Q}^{\top}\mathbf{X}

とおくと、上式は

\mathbf{X} \approx \mathbf{Q} \mathbf{B}

という近似的な行列分解を計算しているともいえる。$\mathbf{B} \in \mathbb{R}^{r \times m}$は元の行列に対して小規模な行列であり、$\mathbf{B}$に対する各種の計算コストは小さい。したがって、直交行列$\mathbf{Q}$が求められれば次元削減が達成される。

上式は上位$r$個の特異値のみを残すTruncated SVDと本質的には同じであり、実際、特異値分解を$\mathbf{X} = \mathbf{U \Sigma V^{\top}}$とすると、直交行列は$\mathbf{Q} = \mathbf{U}_{[:,:r]}$で求められる。しかし、大規模行列の場合は特異値分解の計算コストが非常に大きくなるため、低ランク近似のために特異値分解を計算することは本末転倒といえる。

そこで、乱数ベクトルを使って直交基底を近似的に求めることを考える。原理はとてもシンプルである。まず、次式によって$\mathbf{X}$の列空間のベクトル$\mathbf{y}$を生成する。

\mathbf{y} = \mathbf{X} \mathbf{w} = \sum_{i=1}^{m} w_{i} \mathbf{x}_{i}

ここで、$\mathbf{x}_{i} \in \mathbb{R}^{n}$は$\mathbf{X}$の$i$列目の列ベクトルであり、上式は$\mathbf{y}$が$\mathbf{x}_{i}$によって張られる空間に含まれることを表す。これを$r$本まとめて行列で書くと、

\mathbf{Y} = \mathbf{X} \mathbf{\Omega}

ここで、$\Omega \in \mathbb{R}^{m \times r}$は乱数行列であり、正規乱数や一様乱数から生成する。

さて、相異なる乱数ベクトルは非常に高い確率で線形独立となるから、上式で生成した行列$\mathbf{Y} \in \mathbb{R}^{n \times r}$の各列は非常に高い確率で$\mathbf{X}$の列空間の基底となることが期待できる。したがって、求める直交行列はQR分解

\mathbf{Y} = \mathbf{Q} \mathbf{R}

で計算することができる。なお、実際の計算では以下のような精度向上のためのテクニックも合わせて使用する。

オーバーサンプリング

実際の計算では$r$はターゲットランクとしてユーザーが指定する。しかし、データ行列のランクは厳密に$r$になるわけではなく、0でない微小な特異値が含まれることによって誤差が生じる。そこで、精度向上のためのために使用する乱数ベクトルの本数を$l = r + p$として$p$本だけ多めに設定する。おおむね$p = 10$程度で十分である。

べき乗反復

通常、低ランク行列の特異値を降順に並べると、その値は急激に小さくなるため、少ないターゲットランクでもよい近似が得られることが期待できる。しかし、特異値の減衰速度がゆるやかな行列の場合、基底ベクトルの近似精度が下がってしまう。そのような場合、微小な特異値を0に近づけるため、次のような前処理を行う。

\mathbf{Y} = \left( \left( \mathbf{X}\mathbf{X}^{\top} \right)^{q} \mathbf{X} \right) \mathbf{\Omega}

通常、$q = 2 \sim 4$程度で十分である。見ての通り、上式は余分な行列積が必要となるため、行列の規模が大きくなると計算コストも大きくなる。そのため、特異値の減衰が速いことが事前にわかっている問題であれば、あえて使用する必要はないと思われる。

メモリ削減のための逐次処理

基本的なアルゴリズムは上述した通りだが、$\mathbf{X}$がメモリに格納できない程大規模な行列であった場合、$\mathbf{Y} = \mathbf{X} \mathbf{\Omega}$のような行列積を通常通り計算することはできない。

そこで、$\mathbf{X}$を列ベクトルごとに分割して保存しておき、メモリに格納可能な複数本のベクトルのみを読み出して来て逐次的に処理する必要がある。たとえば、上記の行列積であれば、

\begin{align}
& \mathbf{Y} := \mathbf{0}^{n \times l} \\
& \mathrm{for} \ J := J_{1}, J_{2}, \cdots, J_{b} \\
& \qquad \mathbf{Y} := \mathbf{Y} + \mathbf{X}_{[:,J]} \mathbf{\Omega}_{[J,:]}
\end{align}

のように計算できる。ここで、$J_{i}$は列のインデックスを$b$個に分割したリストである。たとえば、$m = 100,\ b = 10$とすると、$J_{1} = [1,2, \cdots, 10]$、$J_{10} = [91,92, \cdots, 100]$である。$b = 1$のときは列ごとの更新となり、いわゆるrank-1 updateとなる。

$\mathbf{X}_{[:,J]}$のアクセスにはファイルの読み込みが必要であるから、それが増えるほど計算コストもかさむ。たとえば、上述のべき乗反復を行うと$q$回だけ余分にデータアクセスが生じるため、近似精度と計算コストのトレードオフから実施するかどうかを決める必要がある。

サンプル

アルゴリズムのまとめも兼ねて、以下にpythonによるサンプルコードを示す。データ行列の列ブロックに対してスライスでアクセスしているが、この部分は実際にはファイルの読み込みが必要になる。

import numpy as np
from scipy.linalg import qr

def calc_qb(X, r, p=10, q=2, b=10) :
    n, m = X.shape
    l = r + p
    Omega = np.random.randn(m, l)
    col_idx = np.array_split(range(m), b)
    
    # 基底の生成
    Y = np.zeros((n,l))
    for col in col_idx :
        Y += X[:,col] @ Omega[col,:]
    
    # べき乗反復
    for k in range(q) :
        Q, _ = qr(Y, mode='economic')
        Y = np.zeros((n,l))
        for col in col_idx :
            x = X[:,col]
            z = x.T @ Q
            Y += x @ z
    
    # 低ランク近似
    Q, _ = qr(Y, mode='economic')
    B = np.zeros((l,m))
    for col in col_idx :
        x = X[:,col]
        B[:,col] = Q.T @ x
        
    return Q, B

テスト

簡単なテスト計算をやってみる。次のような空間と時間の関数を考える。

$$
f(x,t) = \sum_{k=1}^{5} \frac{1}{k} \cos(k \pi t) \sin(k \pi x)
$$

この関数の値を適当な離散点で評価した結果を行列としてまとめる。このとき、行方向に空間、列方向に時間としてデータを並べると、この行列のランクは$5$になる。明らかに上式は$\sin(k \pi x)$という基底の線形結合になっているからだ。

実際、以下のコードで行列を作成してランクを確認してみると、結果は5になった。

t0 = 0
te = 1.0
L = 1.0
nx = 1000
nt = 100
kmax = 5

x = np.linspace(0,L,nx)
t = np.linspace(t0,te,nt)
X, T = np.meshgrid(x,t)

def f(x,t) :
    res = 0
    for k in range(1,kmax+1) : 
        res += 1.0/k*np.cos(k*np.pi*t)*np.sin(k*np.pi*x)
    return res

D = f(X,T).T

# 行列のサイズとランクの確認
print(D.shape)
print(np.linalg.matrix_rank(D))

# 実行結果
# (1000, 100)
# 5

ちなみにnp.linalg.matrix_rankは特異値の大きさから数値的ランクを判定する関数なので、行列が大きくなると気軽には使えないものだったりする。

最後に低ランク近似を計算してみると、以下のようになった。

r = kmax
Q, B = calc_qb(D,r,p=10,q=2,b=10)
print(np.linalg.norm(D - Q@B))

# 実行結果
# 1.0385441066479628e-13

差のノルムは丸め誤差の影響か毎回異なったが、概ね$10^{-13}$のオーダーであり、十分な近似精度があるといえる。今回の場合は三角関数という直交基底から生成したデータであったためか、オーバーサンプリングやべき乗反復を行わなくても結果の精度に変化はなかった。より複雑な実データであればこうしたパラメータも影響するものと思われる。

まとめ

Random Projectionを用いた低ランク近似手法についてまとめた。冒頭の参考記事の通り、この手法は特異値分解の高速化などに応用されており、大規模データの次元削減の手法として有用である。今後も機械学習の分野をはじめ、幅広い分野で重宝されるものと思われる。

追記

書き忘れていたが、上記を利用して高速に特異値分解を計算する関数がscikit-leranやdaskなどに用意されているので、特異値分解の高速化が目的であれば上記のアルゴリズムを自分で実装する必要はなかったりする。ありがたや、ありがたや。