C++ と liboctave でPCAを適用するプログラムについて．固有値分解を使う方法と特異値分解を使う方法を紹介する．

固有値分解を使う †

データの共分散行列を求めて，その行列を固有値分解するという，テキストに載っているやりかた．

固有値分解された結果をソートしたベクトル std::vector<int> idx を作っているが，これは固有値分解された結果がソートされていないから．

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

// 共分散行列を計算
// 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);
}

// greater 関数オブジェクト．ただし比較に ref(i) を使う
template <typename T>
struct __type_ref_greater
{
  const T &ref;
  __type_ref_greater (const T &r) : ref(r) {};
  bool operator() (int i, int j)
    {return ref(i)>ref(j);};
};

// x を大きい順に並びかえたインデックスを得る
void get_sorted_indices (const ColumnVector &x, std::vector<int> &idx)
{
  idx.resize (x.length());
  int i(0);
  for(std::vector<int>::iterator itr(idx.begin()); itr!=idx.end(); ++itr,++i) *itr=i;
  __type_ref_greater<ColumnVector> ref_greater(x);
  std::sort(idx.begin(), idx.end(), ref_greater);
}

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;
  // 固有値分解
  EIG eig(Cov);
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // 大きい順にソートしたインデックスを得る
  vector<int> idx;
  get_sorted_indices (real(eig.eigenvalues()), idx);
  // もとのデータを1次元に圧縮
  ofstream ofs("reduced1.dat");
  Matrix Tr (eig.eigenvectors().rows(),1); // 変換行列
  Tr.insert (real(eig.eigenvectors().column(idx[0])), 0,0);
  ofs<< Data * Tr * Tr.transpose();
  ofs.close();
  // もとのデータを2次元に圧縮
  ofs.open("reduced2.dat");
  Tr.resize (eig.eigenvectors().rows(),2);
  Tr.insert (real(eig.eigenvectors().column(idx[0])), 0,0);
  Tr.insert (real(eig.eigenvectors().column(idx[1])), 0,1);
  ofs<< Data * Tr * Tr.transpose();
  ofs.close();
  return 0;
}

特異値分解を使う †

データの共分散行列を求めて，これを特異値分解するというやりかた．一般的に，固有値分解のよりも精度がいいらしい (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;
}

サンプルデータ生成 †

PCA (主成分分析) にかけるデータの生成プログラム．

3次元空間上の平面データを作って， x,y,z それぞれにガウスノイズを加えてデータを生成するプログラムです． gnuplot とかで表示するとボリュームがあるように見えます． コンパイルには liboctave 以外に， boost が必要です．

#include <octave/config.h>
#include <octave/Matrix.h>
#include <iostream>
#include <boost/random.hpp>

namespace loco_rabbits
{

template< class T >
inline T _square( const T &val )
{
  return val*val;
}

class TAddGaussianNoise
{
private:
  boost::variate_generator< boost::mt19937, boost::normal_distribution<> >  ndrand;
public:
  TAddGaussianNoise(const double &std_dev=1.0)
    : ndrand (boost::mt19937(static_cast<unsigned long>(std::time(0))),
              boost::normal_distribution<>(0.0, _square(std_dev)))  {};
  void operator() (ColumnVector &x)
    {for(int r(x.dim1()-1); r>=0; --r)  x(r)+=ndrand();};
};

inline double u_rand (const double &max)
{
  return (max)*static_cast<double>(rand()) / static_cast<double>(RAND_MAX);
}

inline double u_rand (const double &min, const double &max)
{
  return u_rand(max - min) + min;
}

};  // end of namespace loco_rabbits

using namespace std;
using namespace loco_rabbits;

int main(int argc, char**argv)
{
  srand((unsigned)time(NULL));
  ColumnVector x(3);
  TAddGaussianNoise  addGaussianNoise(1.0);
  for(int i(0); i<1000; ++i)
  {
    x(0)=u_rand(-5.0,5.0);
    x(1)=x(0)+u_rand(-10.0,10.0);
    x(2)= 0.1*x(0) + 0.2*x(1);
    addGaussianNoise(x);
    cout<<x.transpose()<<endl;
  }

  return 0;
}