2012-09-30

C++の行列演算ライブラリのEigenを使うと比較的簡単にアルゴリズムを実装することができます。
次元削減手法の一つの局所線形埋め込み(LLE)を実装してみました。
ずっと以前に書いていたコードを放置していたのですが、ここで紹介します。

ソースと使い方

本体部分は以下のようになります。
残念ながらまだ疎行列版(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を計算します。
式で書くと
$min_{W} ||x^i -\sum_{j \in neighbor of xi} W_{ij}x^j ||^2$
を満たすのがWとなります。Wの拘束条件
$\sum_{j} W_ij =1$
を用いると
$||x^i -\sum_{j \in neighbor of xi} W_{ij}x^j ||^2 =||\sum_{j \in neighbor of xi}W_{ij}(x^i - x^j) ||^2$
$= \sum_{j}\sum_{k} W_{ij}W_{ik}(x^i - x^j)\cdot (x^i - x^k)$
となります。 $(x^i - x^j)\cdot (x^i - x^k)$ を各要素として持つ行列がlocal covariance行列lcに相当します。
さらに拘束条件に未定乗数をかけて加えたもの
$\sum_{j} \sum_{k} W_{ij}W_{ik}(x^i - x^j)\cdot (x^i - x^k) +\lambda(\sum_{j} W_ij - 1)$
をWijに対して微分した値が0になる場合が極値となるので
$\sum_{k \in neighbor of xi} W_{ik}(x^i - x^j)\cdot (x^i - x^k) -\lambda=0$
となります。λはあとでWijを規格化することを考えると1としてよく、結局
$\sum_{k \in neighbor of xi} W_{ik}(x^i - x^j)\cdot (x^i - x^k) =1$
を解けばいいことになります。コードではこれをsolve(ここではLU分解)で計算しています。

各点iに対してWが求まるとこれを滑らかにつなげて１枚の低次元空間にしなければ行けません。各点iの低次元での座標をβiとするとデータ点xiと同様に
βiは近隣点のβjによってWで重みづけられて表現されます。これを各iに関して足したものの総和が最小になる場合
$min_{\beta} \sum_{i} (\beta_i- \sum_{k \in neighbor of xi} W_{ij} \beta_j)^2$
を求めなければ行けません。
$||\beta||^2=1$ の条件で再び未定乗数法を使うと
$||(I-W)\beta||^2-\lambda (||\beta||^2 -1 )$
これをβで微分した
$(I-W)^T(I-W)\beta = \lambda \beta$
を解く問題に帰着されます。これは行列 $(I-W)^T(I-W)$ の固有値問題なので解くことができます。ただし最大固有値は自明な固有ベクトル(すべての係数が同じ)となるのでこれをのぞいた結果が低次元空間での各点の座標となります。

結果

３次元空間内のスイスロール状のデータ

を低次元(2次元)に埋め込みます。巻きの位相を色で表示しています。
埋め込む際に近隣点(neighbor)の数をいくつにするかという任意性があり、それによって埋め込みがうまくできるかどうかが決定してしまいます。
neighbor=10の場合

neighbor=20の場合

neighbor=50の場合