数値計算ライブラリ liboctave を用いたC++プログラミング (Windows / Linux)

Main Page >> Memo/Link >> liboctave

liboctave for windows MSVC or linux g++ 2008

2008年夏辺りのliboctaveの使い方についてのメモ

liboctave以外には,IT++などのライブラリがあり,(これはliboctaveと異なり,そもそもAPIとして設計されたものなので)ドキュメントが豊富です.この情報もメモしておきます.

以下に2001,2002年あたりの古い情報がありましたが,あまりにも古くなってきたので使用例を除いて削除しました.

Octave C++で数値計算例

固有値分解

EIG

00001 #include <iostream>
00002 #include <octave/config.h>
00003 #include <octave/Matrix.h>
00004 
00005 using namespace std;
00006 
00007 int main()
00008 {
00009   Matrix m(2, 2);
00010   m(0,0) = 3; m(0,1) = 2;
00011   m(1,0) = 2; m(1,1) = 0;
00012   cout << "Original Matrix" << endl <<  m << endl;
00013 
00014   EIG eig(m);
00015   cout << "Eigen Vectors" << endl <<  eig.eigenvectors() << endl;
00016   cout << "Eigen Values" << endl <<  eig.eigenvalues() << endl;
00017   cout << "Recomposed Matrix" << endl << eig.eigenvectors()
                                       *ComplexMatrix(ComplexDiagMatrix(eig.eigenvalues()))
                                       *eig.eigenvectors().inverse() << endl;
00018 
00019   return 0;
00020 }

特異値分解

00001 #include <iostream>
00002 #include <octave/config.h>
00003 #include <octave/Matrix.h>
00004 
00005 using namespace std;
00006 
00007 int main()
00008 {
00009   Matrix m(3, 2);
00010   m(0,0) = 3; m(0,1) = 2;
00011   m(1,0) = 2; m(1,1) = 0;
00012   m(2,0) = 4; m(2,1) = 2;
00013   cout << "Original Matrix" << endl <<  m << endl;
00014 
00015   SVD svd(m);
00016   cout << "Left Singular Matrix" << endl <<  svd.left_singular_matrix() << endl;
00017   cout << "Singular Values" << endl <<  svd.singular_values() << endl;
00018   cout << "Right Singular Matrix" << endl <<  svd.right_singular_matrix() << endl;
00019   cout << "Recomposed Matrix" << endl << svd.left_singular_matrix()
                                              *Matrix(DiagMatrix(svd.singular_values()))
                                              *svd.right_singular_matrix() << endl;

00020 
00021   return 0;
00022 }

Octave C++で多変量解析

平均と共分散・相関係数

00001 #include <iostream>
00002 #include <octave/config.h>
00003 #include <octave/Matrix.h>
00004 
00005 #include "octaveio.h"
00006 
00007 using namespace std;
00008 
00009 int main(int argc, char **argv)
00010 {
00011   Matrix X;
00012   "X.dat" >> X;         // 各列がサンプルデータ
00013   int dim = X.rows();    // ベクトルの次元数
00014   int nsample = X.cols(); // サンプルベクトル数
00015   cout << "X\n" << X << endl;
00016 
00017   // 平均ベクトルの計算
00018   ColumnVector mean(dim);
00019   mean = (X.sum(1) / nsample).column(0); // X.sum(1)/nsample だとMatrix型のまま
00020   cout << "mean\n" << mean << endl;
00021 
00022   // 共分散行列の計算
00023   Matrix Cov(dim, dim);
00024   Matrix Xd(dim, nsample); // 平均ベクトルを各列から引いた行列
00025   for(int i=0; i < nsample; i++){
00026     Xd.insert(X.column(i) - mean, 0, i);
00027   }
00028   Cov = Xd*Xd.transpose() / nsample; // 不偏量を考えるなら nsample - 1 で割る
00029   cout << "Cov\n" << Cov << endl;
00030 
00031   // 相関係数行列の計算
00032   Matrix Rate(dim, dim);
00033   Matrix Xn(dim, nsample); // 各行を正規化した行列
00034   for(int i=0; i < dim; i++){
00035     Xn.insert(X.row(i)/sqrt(X.row(i)*X.row(i).transpose()), i, 0);
00036   }
00037   Rate = Xn*Xn.transpose();
00038   cout << "Rate\n" << Rate << endl;
00039 
00040   return 0;
00041 }

主成分分析

00001 #include <iostream>
00002 #include <octave/config.h>
00003 #include <octave/Matrix.h>
00004 
00005 #include "octaveio.h"
00006 
00007 using namespace std;
00008 
00009 int main(int argc, char **argv)
00010 {
00011   Matrix X;
00012   "X.dat" >> X;         // 各列がサンプルデータ
00013   int dim = X.rows();    // ベクトルの次元数
00014   int nsample = X.cols(); // サンプルベクトル数
00015   cout << "X\n" << X << endl;
00016 
00017   ColumnVector mean(dim);  // 平均ベクトル
00018   Matrix Xd(dim, nsample); // 平均ベクトルを各列から引いた行列
00019   Matrix Cov(dim, dim);    // 共分散行列
00020   mean = (X.sum(1) / nsample).column(0);
00021   for(int i=0; i < nsample; i++){
00022     Xd.insert(X.column(i) - mean, 0, i);
00023   }
00024   Cov = Xd*Xd.transpose() / nsample;
00025 
00026   // 主成分分析(固有値分解を用いて)
00027   EIG eig(Cov);
00028   cout << "eigen vectors\n" << eig.eigenvectors() << endl; // 主成分
00029   cout << "eigen values\n" << eig.eigenvalues() << endl;   // 固有値(昇順)
00030 
00031   // 主成分分析(特異値分解を用いて)
00032   SVD svd(Xd);
00033   cout << "eigen vectors\n" << svd.left_singular_matrix() << endl;
00034   ColumnVector eigenvalues(dim);
00035   for(int i=0; i < dim; i++){ // Xdを直接分解したので共分散行列の固有値のサンプル数倍
00036     eigenvalues(i) = svd.singular_values()(i,i)*svd.singular_values()(i,i) / nsample;
00037   }
00038   cout << "eigen values\n" << eigenvalues << endl; // 固有値(降順)
00039 
00040   return 0;
00041 }

正準相関分析

00001 #include <iostream>
00002 #include <octave/config.h>
00003 #include <octave/Matrix.h>
00004 
00005 #include "octaveio.h"
00006 
00007 using namespace std;
00008 
00009 int main()
00010 {
00011   // 2組の2次元変数群(x_i, y_i (i=1,...,nsample))を読み込む
00012   // ただし平均が0であることを仮定する
00013   Matrix X;
00014   Matrix Y;
00015   "X.dat" >> X;
00016   "Y.dat" >> Y;
00017 
00018   cout << "X\n" <<  X << endl;
00019   cout << "Y\n" <<  Y << endl;
00020 
00021   int dim = X.rows();
00022   int nsample = X.cols();
00023 
00024   // それぞれの変数群を白色化するために,主成分分析を行う
00025   Matrix Ux, Uy;
00026   DiagMatrix Sx_inv, Sy_inv;
00027   EIG eig_X(X*X.transpose()/nsample);
00028   EIG eig_Y(Y*Y.transpose()/nsample);
00029   Ux = real(eig_X.eigenvectors());
00030   Uy = real(eig_Y.eigenvectors());
00031   ColumnVector Sx = real(eig_X.eigenvalues());
00032   ColumnVector Sy = real(eig_Y.eigenvalues());
00033   Sx_inv.resize(dim,dim);
00034   Sy_inv.resize(dim,dim);
00035   for(int i=0; i < dim; i++){
00036     Sx_inv(i,i) = 1.0/sqrt(Sx(i));
00037     Sy_inv(i,i) = 1.0/sqrt(Sy(i));
00038   }
00039 
00040   // 白色化する行列 (XX'/nsample)^(-1/2), (YY'/nsample)^(-1/2)
00041   Matrix Wx = Ux*(Matrix)Sx_inv*Ux.transpose();
00042   Matrix Wy = Uy*(Matrix)Sy_inv*Uy.transpose();
00043 
00044   // 白色化された変量群
00045   Matrix Xw = Wx * X;
00046   Matrix Yw = Wy * Y;
00047 
00048   // 白色化された変量群間の相互相関を特異値分解する
00049   SVD svd_R(Yw*Xw.transpose()/nsample, SVD::economy);
00050   // 係数行列
00051   Matrix A = svd_R.right_singular_matrix().transpose() * Wx;
00052   Matrix B = svd_R.left_singular_matrix().transpose() * Wy;
00053   cerr << "A\n" << A << endl;
00054   // 正準相関係数
00055   DiagMatrix Sr = svd_R.singular_values();
00056   cerr << "Sr\n" << Sr << endl;
00057 
00058   //--------------------------------------------------
00059   // 検証
00060   //--------------------------------------------------
00061   // 正準変量群
00062   Matrix Xc = A * X;
00063   Matrix Yc = B * Y;
00064   cout << "Xc\n" << Xc << endl;
00065   cout << "Yc\n" << Yc << endl;
00066   // 正準変量群間の相互相関は,正準相関係数が対角に並んだ行列
00067   cout << "Xc*Yc'/nsample\n" << Xc*Yc.transpose()/nsample << endl;
00068 
00069   return 0;
00070 }