PCA の練習

練習に http://archive.ics.uci.edu/ml/datasets/Iris のデータをさくっと PCA してみました。

#include <iostream>
#include <string>
#include <utility>
#include <vector>
#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <eigen3/Eigen/Dense>
// #include <redsvd/redsvd.hpp>
using namespace std;
using namespace boost;

int main() {
  vector<pair<string, vector<float>>> data;
  for (string line; getline(cin, line); ) {
    vector<string> row;
    split(row, line, is_any_of(","));
    if (row.size() < 2) continue;

    auto p = make_pair(move(row.back()), vector<float>());
    row.pop_back();
    for (const auto& val : row) p.second.push_back(lexical_cast<float>(val));
    data.push_back(move(p));
  }

  Eigen::MatrixXf data_matrix(data.size(), data[0].second.size());
  for (size_t i = 0; i < data.size(); ++i) {
    for (size_t j = 0; j < data[i].second.size(); ++j) {
      data_matrix(i, j) = data[i].second[j];
    }
  }

  Eigen::JacobiSVD<Eigen::MatrixXf> svd(
      data_matrix, Eigen::ComputeThinU | Eigen::ComputeThinV);
  const auto& pca_dat = svd.matrixU() * svd.singularValues().asDiagonal();
  // const auto& pca_dat = REDSVD::RedPCA(data_matrix, 2).scores();
  for (int i = 0; i < pca_dat.rows(); ++i) {
    cout << data[i].first;
    for (int j = 0; j < pca_dat.cols(); ++j) cout << ',' << pca_dat(i, j);
    cout << endl;
  }
}

g++4.6 -std=c++0x でコンパイルできます。RedSVD を使う方はコメントアウトしてあります (svd, pca_dat の定義を消してコメントを戻せば RedSVD 版になる)。RedSVD の方は上位 2 次元だけ出力しています。Eigen は魔法ですね。チュートリアルも親切。でもリファレンス読むのは難しい。