Eigenの特異値分解を利用した主成分分析

今日はEigenを利用してデータ分析の手法として良く知られた主成分分析を実装してみます．基本的な部分は以下の記事と同じで，C言語やVBAから利用可能なDLLをビルドしようというものです．

EigenのMapオブジェクトを使ってC言語の配列を操作する

さて，主成分分析は $m$ 次元のデータベクタ ${\bf a}_i$ ，( $i=1,\dots,n$ )について，その共分散行列の固有値問題を解くことで行うことができます．しかしながら，特に $m>>n$ である場合，共分散行列の計算に時間がかかるなどの問題があります．こうした場合は特異値分解を利用した方法が推奨されます．具体的には平均減算を行ったデータベクタを横に並べた以下の $m\times n$ 行列 ${\bf M}$ の特異値問題を解くものです．
${\bf M} = \begin{bmatrix} ({\bf a}_1-{\bf\mu})& ({\bf a}_2-{\bf\mu})& \dots & ({\bf a}_n-{\bf\mu})\end{bmatrix}$
ただし， ${\bf \mu}$ は ${\bf a}_i$ ，( $i=1,\dots,n$ )の平均ベクタです．この ${\bf M}$ の左特異行列が主成分ベクタに，特異値の二乗が主成分に対応します．

Eigenを利用する場合，実装は非常に簡単であり，その骨子は以下のようになります．

配列を受け取ってMapオブジェクトで行列とみなす
rowwise()とmean()で平均ベクタを計算
平均減算した行列をJacobiSVDオブジェクトで特異値分解
左特異行列から主成分ベクタを，特異値の二乗から主成分を得る

Mapを利用することでコピー操作が不要になり，また，組み込み関数を利用するだけでfor文などを利用した反復処理さえ不要になります．以下はデータベクタを列優先で $m\times n$ 行列として与え，平均ベクタと $r$ 本の主成分ベクタ，それに対応する主成分を計算する関数です．

// This code is licensed under the GPL.
// simplepca.cpp
#include "eigen_vba.h"

// Eigen
#include <Eigen/Dense>
using namespace Eigen;

#ifdef _DEBUG
#include <iostream>
using namespace std;
#endif

// 特異値分解を用いた単純な主成分分析
EIGEN_VBA_API int CALL_VBADLL SimplePCA_(int m, int n, int r,
					 double *_a,
					 double *_mu, double *_lambda, double *_v){
  
  // Map オブジェクト (既定通りの列優先)
  Map<MatrixXd> a(_a, m, n);
  Map<VectorXd> mu(_mu, m);
  Map<VectorXd> lambda(_lambda, r);
  Map<MatrixXd> u(_v, m, r);

  // 平均を計算する
  mu = a.rowwise().mean();
 
  // 平均減算した行列を特異値分解する
  // (ComputeThinUを与えて特異値と左特異ベクタだけを計算する)
  JacobiSVD<MatrixXd> SVD((a.colwise() - mu), ComputeThinU);
  
  // 特異値の二乗と左特異ベクタを返す
  lambda = (SVD.singularValues()).block(0, 0, r, 1).array().square();
  u =  (SVD.matrixU()).block(0, 0, m, r);

#ifdef _DEBUG
  cout << "A=[" << a << "];" << endl;
  cout << "mu=[" << mu << "];" << endl;
  cout << "lambda=[" << lambda << "];" << endl;
  cout << "V=[" << u << "];" << endl;
#endif

  return 0;
}

ビルドの方法は前回の記事を参照してください．乱数で構成した $4\times 3$ 行列の第二主成分までの計算結果は以下です．

A=[0.423394  0.104033 0.561979
   0.988998  0.477972 0.18587
   0.0909832 0.281566 0.924881
   0.155299  0.271587 0.481722];

mu=[0.363135
    0.550947
    0.432477
    0.302869];

lambda=[0.749016
        0.128381];

V=[ 0.170522 0.830388
   -0.631153 0.5002
    0.70674  0.245461
    0.270346 0.00231725];

平均ベクタ，主成分，主成分ベクタが正しく計算されていることが分かります．