Eigenで局所線形埋め込み
C++の行列演算ライブラリのEigenを使うと比較的簡単にアルゴリズムを実装することができます。
次元削減手法の一つの局所線形埋め込み(LLE)を実装してみました。
ずっと以前に書いていたコードを放置していたのですが、ここで紹介します。
ソースと使い方
https://github.com/xiangze/LLE_eigen
本体部分は以下のようになります。
残念ながらまだ疎行列版(RedSVD)には対応していません。
#include <eigen3/Eigen/Dense> #include <eigen3/Eigen/SVD> #include <eigen3/Eigen/LU> #ifdef _SPARSE #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET #include <eigen3/Eigen/Sparse> #include <redsvd/redsvd.hpp> #endif #include "iostream" using namespace Eigen; #define DBL_MAX (100000) #define MINVAL (.00001) //data :input data //d; distance matrix of data //out: output (dimenstion compressed expression of data) //n: neigher number void lle_eigen(MatrixXf & dist,MatrixXf & data,MatrixXf & out, int nnum){ #ifdef _SPARSE SparseMatrix<float> W(dist.rows(), dist.cols()); W.setZero(); #else MatrixXf W= MatrixXf::Zero(dist.rows(), dist.cols()); #endif VectorXf::Index m; VectorXf vd; for (int i = 0;i < data.rows();i++){ // neigherhood matrix MatrixXf nm(nnum,data.cols()); vd=dist.row(i); vd(i) = DBL_MAX;//self=max for (int r = 0;r < nnum;++r){ vd.minCoeff(&m); for (int c = 0;c < data.cols();c++) nm(r, c) = data(m, c) - data(i, c); vd(m)=DBL_MAX; } //local covariance MatrixXf lc = nm*nm.transpose(); if(nnum>data.cols()) lc+=MINVAL*MatrixXf::Identity(lc.rows(),lc.cols()); VectorXf v1=VectorXf::Constant(lc.rows(),1.0); FullPivLU<MatrixXf> lu(lc); VectorXf vw=lu.solve(v1);//lc*vw=v1 double wsum = vw.sum(); vd = dist.row(i); vd(i) = DBL_MAX;//self=max for (int r=0;r<vw.size();r++){ vd.minCoeff(&m); #ifdef _SPARSE W.insert(i, m) = vw(r)/wsum; #else W(i, m) = vw(r)/wsum; #endif vd(m) = DBL_MAX; } } #ifdef _SPARSE SparseMatrix<float>sm = -W; for (int r = 0;r < sm.rows();r++) sm.insert(r,r)= sm(r,r)+1; sm = sm.transpose()*sm; REDSVD::RedSVD redSVD(sm, out.cols()+1); // for (int r = 0;r < out.rows();r++) // for (int c = 0;c < out.cols();c++) // out(r, c) = evecs(r,c +1); #else MatrixXf sm = -W; sm+= MatrixXf::Identity(sm.rows(),sm.cols()); sm = sm.transpose()*sm; SelfAdjointEigenSolver<MatrixXf> eigensolver(sm); MatrixXf evecs=eigensolver.eigenvectors(); for (int r = 0;r < out.rows();r++) for (int c = 0;c < out.cols();c++) out(r, c) = evecs(r,c +1); #endif return; }
処理の手順は
- 重み行列の計算
- 対角化
に分かれます。
まずデータの各点に対しnnumで指定された数だけ近隣点(neighbor)をとってきて、それを用いて各点を表現できるような重みの行列Wを計算します。
式で書くと
を満たすのがWとなります。Wの拘束条件
を用いると
となります。を各要素として持つ行列がlocal covariance行列lcに相当します。
さらに拘束条件に未定乗数をかけて加えたもの
をWijに対して微分した値が0になる場合が極値となるので
となります。λはあとでWijを規格化することを考えると1としてよく、結局
を解けばいいことになります。コードではこれをsolve(ここではLU分解)で計算しています。
各点iに対してWが求まるとこれを滑らかにつなげて1枚の低次元空間にしなければ行けません。各点iの低次元での座標をβiとするとデータ点xiと同様に
βiは近隣点のβjによってWで重みづけられて表現されます。これを各iに関して足したものの総和が最小になる場合
を求めなければ行けません。
の条件で再び未定乗数法を使うと
これをβで微分した
を解く問題に帰着されます。これは行列の固有値問題なので解くことができます。ただし最大固有値は自明な固有ベクトル(すべての係数が同じ)となるのでこれをのぞいた結果が低次元空間での各点の座標となります。
結果
を低次元(2次元)に埋め込みます。巻きの位相を色で表示しています。
埋め込む際に近隣点(neighbor)の数をいくつにするかという任意性があり、それによって埋め込みがうまくできるかどうかが決定してしまいます。
neighbor=10の場合
neighbor=20の場合
neighbor=50の場合
今回のデータでは20の場合が最も分離がうまくいっているようなのですが、あらかじめ最適な値を知ることは難しいのがLLEの難点です。
Reference
- 提唱者Roweis先生のページmatlabコードがあります。
- モーションデータの様々な多様体学習手法での2次元への埋め込みGSL(GNU Scientifc Library)で実装されています。
ここの記述を参考にさせていただきました。使っているライブラリは違うもののコードの見た感じほとんど同じになってしまいました。。。
こちらではIsomap,MDSなど他の手法ではうまくいっているようなのですが、LLEではなぜか結果が一カ所につぶれて見えてしまいます。
- 人体を典型とした物体の動きの分類には注目が集まっております。多様体学習のような次元削減手法を用いた人体の動きの分類、解析とその応用分野については【コラム】動画から歩き方の個性を測る「歩容解析 (Gait Analysis )」:犯罪捜査への応用にむけてで詳しく考察されています。
- scikit-learnでの手書き文字データの埋め込み
- カーネル多変量解析
カーネル多変量解析―非線形データ解析の新しい展開 (シリーズ確率と情報の科学)
- 作者: 赤穂昭太郎
- 出版社/メーカー: 岩波書店
- 発売日: 2008/11/27
- メディア: 単行本
- 購入: 7人 クリック: 180回
- この商品を含むブログ (32件) を見る
他の次元削減手法とともに日本語で詳しく解説されています。