liboctaveで主成分分析 -特異値分解を使う-

C++ と liboctave でPCAをするプログラム．データの共分散行列を求めて，これを特異値分解するというやりかた．一般的に，固有値分解のよりも精度がいいらしい (Principal components analysis, Wikipedia, The Free Encyclopedia)．
Octave-Forge princomp の octave ソースを参照した．

#include <octave/config.h>
#include <octave/Matrix.h>
#include <iostream>
#include <fstream>

// 共分散行列を計算
// x は各行がサンプルデータ
inline Matrix get_covariance (Matrix x)
{
  const int n(x.rows());
  x = x - ColumnVector(n,1.0) * x.sum(0).row(0) / static_cast<double>(n);
  return x.transpose() * x / static_cast<double>(n-1);
}

using namespace std;

int main(int argc, char**argv)
{
  Matrix Data(1000,3); // データサイズは固定
  // データの読み込み (各行がサンプルデータ)
  ifstream ifs("sample.dat");
  ifs >> Data;
  ifs.close();
  // 共分散行列の計算
  Matrix Cov (get_covariance(Data));
  cout<< "Cov= "<<endl<< Cov;
  // 特異値分解
  SVD  svd (Cov, SVD::economy);
  cout<<"svd.singular_values()= "<<endl<< svd.singular_values();
  cout<<"svd.right_singular_matrix()= "<<endl<< svd.right_singular_matrix();
  // もとのデータを1次元に圧縮
  ofstream ofs("reduced1.dat");
  Matrix Tr (Cov.rows(),1); // 変換行列
  Tr.insert (svd.right_singular_matrix().column(0), 0,0);
  ofs<< Data * Tr * Tr.transpose();
  ofs.close();
  // もとのデータを2次元に圧縮
  ofs.open("reduced2.dat");
  Tr.resize (Cov.rows(),2);
  Tr.insert (svd.right_singular_matrix().column(0), 0,0);
  Tr.insert (svd.right_singular_matrix().column(1), 0,1);
  ofs<< Data * Tr * Tr.transpose();
  ofs.close();
  return 0;
}