はじめに
最近ロバスト主成分分析(Robust Principal Component Analysis; RPCA)に興味があり、色々と情報を探していた。 記事もいくつか見つかる。例えば以下の記事など。
上の記事ではロバストPCAを丁寧な説明と共に実装しており、とても勉強になった。 自分でもロバストPCAを実装してみたくなったわけで、記事のタイトルにつながる。
ロバストPCA
密(dense)なデータ行列 \(D\) が与えられたとき、下記の最適化問題を解く。
ここで、\( 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次元空間に射影した結果を示している。
スパース成分Sの「左下部分」を拡大して、低ランク成分Lと再び並べた結果を図2に示す。
低ランク成分Lが対角線上に並び、外れ値以外のスパース成分Sの多くが直交軸上に並ぶのは興味深い。SにはL1ノルムに関する制約がかかっているので、座標成分のうち一方の次元を0につぶす力が働くと考えれば、もっともな結果ではある。
おわりに
cvxpyでさくっと実装して結果を確認できるのはありがたい。今後は他の色々なデータで試してみる予定である。特にデータサイズに対して体感どれくらい「遅い」のかを把握しておきたい。