[go: up one dir, main page]
More Web Proxy on the site http://driver.im/

ロバスト主成分分析をcvxpyで実装してみた話

はじめに

最近ロバスト主成分分析(Robust Principal Component Analysis; RPCA)に興味があり、色々と情報を探していた。 記事もいくつか見つかる。例えば以下の記事など。

qiita.com

上の記事ではロバストPCAを丁寧な説明と共に実装しており、とても勉強になった。 自分でもロバストPCAを実装してみたくなったわけで、記事のタイトルにつながる。

ロバストPCA

密(dense)なデータ行列 \(D\) が与えられたとき、下記の最適化問題を解く。


\begin{align*}
\mathrm{min} \; \|L\|_{∗} + \lambda \|S\|_{1} \;\;\; \mathrm{subject}\; \mathrm{to}\;  D = L + S
\end{align*}

ここで、\( L \) は低ランク行列、\( S \) はスパース行列、 \( || \cdot ||_{∗} \) は核ノルム(行列の特異値の和),\( || \cdot ||_{1} \) はL1ノルム(要素ごとの絶対値の和)を表す。\( \lambda \) は正則化の強さを決めるハイパパラメタである。得られた低ランク行列に対して特異値分解をかけてデータを射影し、次元削減が完了する。

実装

凸最適化問題のソルバの一つであるcvxpyで実装した例はこちら。

このデモスクリプトでは \( [0, 1) \times [0, 1) \) の範囲に各次元独立で乱数を発生させ、それに外れ値を加えたものを密行列としている。

実装のコアの部分を抜粋すると、こうなる:

m, n = dense.shape
low_rank = cp.Variable((m, n))
sparse = cp.Variable((m, n))
objective = cp.Minimize(cp.norm(low_rank, "nuc") + lambda_ * cp.norm(sparse, 1))
problem = cp.Problem(objective, [low_rank + sparse == dense])
problem.solve()

要するに最適化対象となるLとSのための変数(Variable)をそれぞれ用意し、目的関数(objective)と 制約条件(D = L + S)をソルバに与えたら、あとはお任せ、である。とても簡単!

実行結果

図1に実行結果を示す。

図1:実行結果

図1上段の右側、外れ値が入っている様子が見て取れる。図1の中段で低ランク成分とスパース成分に分解した結果が示されている。外れ値がスパース成分に分離されていることから、最適化はうまくいったことが分かる。図1の下段は分解後の成分から再構成した結果と1次元空間に射影した結果を示している。

スパース成分Sの「左下部分」を拡大して、低ランク成分Lと再び並べた結果を図2に示す。

図2: 低ランク成分とスパース成分の分解(Sの拡大図を添えて)

低ランク成分Lが対角線上に並び、外れ値以外のスパース成分Sの多くが直交軸上に並ぶのは興味深い。SにはL1ノルムに関する制約がかかっているので、座標成分のうち一方の次元を0につぶす力が働くと考えれば、もっともな結果ではある。

おわりに

cvxpyでさくっと実装して結果を確認できるのはありがたい。今後は他の色々なデータで試してみる予定である。特にデータサイズに対して体感どれくらい「遅い」のかを把握しておきたい。