C++ の行列ライブラリ Eigen で SVD と LSI
C++ の行列ライブラリ Eigen(アイゲン)で SVD と LSI。
インストール
% tar zvxf eigen-eigen-3.0.3.tar.bz2 % cd eigen-eigen-3.0.3 % mkdir build_dir % cd build_dir % cmake ../ % sudo make install
サンプルを動かす
Eigen: Getting started にあるサンプルを動かしてみる。インクルードパスを追加するのを忘れずに(例：-I /usr/local/include/eigen3)。
#include <Eigen/Core> #include <iostream> using Eigen::MatrixXd; int main() { MatrixXd m(2, 2); m(0, 0) = 3; m(1, 0) = 2.5; m(0, 1) = -1; m(1, 1) = m(1, 0) + m(0, 1); std::cout << m << std::endl; }
動いた。
% ./eigen 3 -1 2.5 1.5
行列の初期化
Comma-initialization を使うと楽。
Matrix<float, 5, 6> c; c << 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1; std::cout << c << std::endl;
SVD
#include <Eigen/SVD>
して
JacobiSVD< Matrix<float, 5, 6> > svd(c, Eigen::ComputeFullU | Eigen::ComputeFullV); std::cout << "singular values" << std::endl << svd.singularValues() << std::endl; std::cout << "matrix U" << std::endl << svd.matrixU() << std::endl; std::cout << "matrix V" << std::endl << svd.matrixV() << std::endl;
できた。コンパイルに時間がかかるようになった。
low-rank approximation
型名が面倒になってきたので typedef する。
typedef Matrix<float, 5, 6> Matrixf56; typedef JacobiSVD<Matrixf56> SVDOfMatrixf56;
C の近似 C_3 を計算する。Singular Values から対角行列作る方法が泥臭い。もっといい方法あるかも。
static Matrixf56 LowRankApprox(const SVDOfMatrixf56::SingularValuesType& sv, const SVDOfMatrixf56::MatrixVType& v, int rank) { Matrixf56 y = Matrixf56::Zero(); for (int i = 0; i < rank && i < sv.size(); i++) { y(i, i) = sv(i); } return y * v.transpose(); }
c_3 = --1.6189 -0.604876 -0.440347 -0.965693 -0.70302 -0.262673 0.456717 0.842566 0.296174 -0.997319 -0.350573 -0.646747 0.356714 -0.954712 0.569498 -0.259686 0.154906 -0.414592 0 0 0 -0 0 0 0 0 0 -0 0 0
できた。
Cosine Similarity
類似度を計算してみよう。
単語文書行列の列（すなわち文書)同士の類似度を計算する。
float Cosine(const Matrixf56& c, int i, int j) { const VectorXf& col_i = c.col(i); const VectorXf& col_j = c.col(j); float a = col_i.transpose() * col_j; float b = col_i.norm() * col_j.norm(); return a / b; }
類似度が近似された C_k と C がどの程度一致するのか比較を出してみた。
rank が 3 と 2 の間に超えられない壁があって k = 3 の近似は比較的良い結果に見える。
c c_4 c_3 c_2 0 1.000000 1.000000 1.000000 1.000000 1 0.408248 0.422235 0.422235 0.781837 2 0.577350 0.630847 0.785422 0.950136 3 0.408248 0.418053 0.418053 0.474432 4 0.577350 0.608732 0.750473 0.740118 5 0.000000 -0.010525 -0.012915 0.110596
コード
#include <Eigen/Core> #include <Eigen/SVD> #include <iostream> using namespace Eigen; typedef Matrix<float, 5, 6> Matrixf56; typedef JacobiSVD<Matrixf56> SVDOfMatrixf56; static Matrixf56 LowRankApprox(const SVDOfMatrixf56::SingularValuesType& sv, const SVDOfMatrixf56::MatrixVType& v, int rank) { Matrixf56 y = Matrixf56::Zero(); for (int i = 0; i < rank && i < sv.size(); i++) { y(i, i) = sv(i); } return y * v.transpose(); } float Cosine(const Matrixf56& c, int i, int j) { const VectorXf& col_i = c.col(i); const VectorXf& col_j = c.col(j); float a = col_i.transpose() * col_j; float b = col_i.norm() * col_j.norm(); return a / b; } int main() { Matrixf56 c; c << 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1; SVDOfMatrixf56 svd(c, Eigen::ComputeFullU | Eigen::ComputeFullV); const SVDOfMatrixf56::MatrixUType& u = svd.matrixU(); const SVDOfMatrixf56::MatrixVType& v = svd.matrixV(); const SVDOfMatrixf56::SingularValuesType& sv = svd.singularValues(); std::cout << "singular values" << std::endl << sv << std::endl; std::cout << "matrix V" << std::endl << v << std::endl; std::cout << "matrix U" << std::endl << u << std::endl; Matrixf56 c_2 = LowRankApprox(sv, v, 2); Matrixf56 c_3 = LowRankApprox(sv, v, 3); Matrixf56 c_4 = LowRankApprox(sv, v, 4); std::cout << "c_3 = " << std::endl << c_3 << std::endl; std::cout << "cosine(0, 1) = " << std::endl << Cosine(c_3, 0, 1) << std::endl; printf(" c c_4 c_3 c_2 \n"); for (int i = 0; i < c_3.row(0).size(); i++) { printf("%d %5f %5f %5f %5f\n", i, Cosine(c, 0, i), Cosine(c_4, 0, i), Cosine(c_3, 0, i), Cosine(c_2, 0, i)); } return 0; }
