数値計算ライブラリ 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

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

using namespace std;

 int main()
{
  Matrix m(2, 2);
  m(0,0) = 3; m(0,1) = 2;
  m(1,0) = 2; m(1,1) = 0;
  cout << "Original Matrix" << endl <<  m << endl;

  EIG eig(m);
  cout << "Eigen Vectors" << endl <<  eig.eigenvectors() << endl;
  cout << "Eigen Values" << endl <<  eig.eigenvalues() << endl;
  cout << "Recomposed Matrix" << endl << eig.eigenvectors()
                                       *ComplexMatrix(ComplexDiagMatrix(eig.eigenvalues()))
                                       *eig.eigenvectors().inverse() << endl;

  return 0;
}

特異値分解

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

using namespace std;

 int main()
{
  Matrix m(3, 2);
  m(0,0) = 3; m(0,1) = 2;
  m(1,0) = 2; m(1,1) = 0;
  m(2,0) = 4; m(2,1) = 2;
  cout << "Original Matrix" << endl <<  m << endl;

  SVD svd(m);
  cout << "Left Singular Matrix" << endl <<  svd.left_singular_matrix() << endl;
  cout << "Singular Values" << endl <<  svd.singular_values() << endl;
  cout << "Right Singular Matrix" << endl <<  svd.right_singular_matrix() << endl;
  cout << "Recomposed Matrix" << endl << svd.left_singular_matrix()
                                              *Matrix(DiagMatrix(svd.singular_values()))
                                              *svd.right_singular_matrix() << endl;


  return 0;
}
ページの先頭へ

Octave C++で多変量解析

平均と共分散・相関係数

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

#include "octaveio.h"

using namespace std;

 int main(int argc, char **argv)
{
  Matrix X;
  "X.dat" >> X;         // 各列がサンプルデータ
  int dim = X.rows();    // ベクトルの次元数
  int nsample = X.cols(); // サンプルベクトル数
  cout << "X\n" << X << endl;

  // 平均ベクトルの計算
  ColumnVector mean(dim);
  mean = (X.sum(1) / nsample).column(0); // X.sum(1)/nsample だとMatrix型のまま
  cout << "mean\n" << mean << endl;

  // 共分散行列の計算
  Matrix Cov(dim, dim);
  Matrix Xd(dim, nsample); // 平均ベクトルを各列から引いた行列
  for(int i=0; i < nsample; i++){
    Xd.insert(X.column(i) - mean, 0, i);
  }
  Cov = Xd*Xd.transpose() / nsample; // 不偏量を考えるなら nsample - 1 で割る
  cout << "Cov\n" << Cov << endl;

  // 相関係数行列の計算
  Matrix Corr(dim, dim);
  Matrix Xn(dim, nsample); // 各行を正規化した行列
  for(int i=0; i < dim; i++){
    Xn.insert(X.row(i)/sqrt(X.row(i)*X.row(i).transpose()), i, 0);
  }
  Corr = Xn*Xn.transpose();
  cout << "Corr\n" << Corr << endl;

  return 0;
}

主成分分析

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

#include "octaveio.h"

using namespace std;

 int main(int argc, char **argv)
{
  Matrix X;
  "X.dat" >> X;         // 各列がサンプルデータ
  int dim = X.rows();    // ベクトルの次元数
  int nsample = X.cols(); // サンプルベクトル数
  cout << "X\n" << X << endl;

  ColumnVector mean(dim);  // 平均ベクトル
  Matrix Xd(dim, nsample); // 平均ベクトルを各列から引いた行列
  Matrix Cov(dim, dim);    // 共分散行列
  mean = (X.sum(1) / nsample).column(0);
  for(int i=0; i < nsample; i++){
    Xd.insert(X.column(i) - mean, 0, i);
  }
  Cov = Xd*Xd.transpose() / nsample;

  // 主成分分析(固有値分解を用いて)
  EIG eig(Cov);
  cout << "eigen vectors\n" << eig.eigenvectors() << endl; // 主成分
  cout << "eigen values\n" << eig.eigenvalues() << endl;   // 固有値（昇順）

  // 主成分分析(特異値分解を用いて)
  SVD svd(Xd);
  cout << "eigen vectors\n" << svd.left_singular_matrix() << endl;
  ColumnVector eigenvalues(dim);
  for(int i=0; i < dim; i++){ // Xdを直接分解したので共分散行列の固有値のサンプル数倍
    eigenvalues(i) = svd.singular_values()(i,i)*svd.singular_values()(i,i) / nsample;
  }
  cout << "eigen values\n" << eigenvalues << endl; // 固有値（降順）

  return 0;
}

正準相関分析

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

#include "octaveio.h"

using namespace std;

 int main()
{
  // 2組の2次元変数群(x_i, y_i (i=1,...,nsample))を読み込む
  // ただし平均が0であることを仮定する
  Matrix X;
  Matrix Y;
  "X.dat" >> X;
  "Y.dat" >> Y;

  cout << "X\n" <<  X << endl;
  cout << "Y\n" <<  Y << endl;

  int dim = X.rows();
  int nsample = X.cols();

  // それぞれの変数群を白色化するために，主成分分析を行う
  Matrix Ux, Uy;
  DiagMatrix Sx_inv, Sy_inv;
  EIG eig_X(X*X.transpose()/nsample);
  EIG eig_Y(Y*Y.transpose()/nsample);
  Ux = real(eig_X.eigenvectors());
  Uy = real(eig_Y.eigenvectors());
  ColumnVector Sx = real(eig_X.eigenvalues());
  ColumnVector Sy = real(eig_Y.eigenvalues());
  Sx_inv.resize(dim,dim);
  Sy_inv.resize(dim,dim);
  for(int i=0; i < dim; i++){
    Sx_inv(i,i) = 1.0/sqrt(Sx(i));
    Sy_inv(i,i) = 1.0/sqrt(Sy(i));
  }

  // 白色化する行列 (XX'/nsample)^(-1/2), (YY'/nsample)^(-1/2)
  Matrix Wx = Ux*(Matrix)Sx_inv*Ux.transpose();
  Matrix Wy = Uy*(Matrix)Sy_inv*Uy.transpose();

  // 白色化された変量群
  Matrix Xw = Wx * X;
  Matrix Yw = Wy * Y;

  // 白色化された変量群間の相互相関を特異値分解する
  SVD svd_R(Yw*Xw.transpose()/nsample, SVD::economy);
  // 係数行列
  Matrix A = svd_R.right_singular_matrix().transpose() * Wx;
  Matrix B = svd_R.left_singular_matrix().transpose() * Wy;
  cerr << "A\n" << A << endl;
  // 正準相関係数
  DiagMatrix Sr = svd_R.singular_values();
  cerr << "Sr\n" << Sr << endl;

  //--------------------------------------------------
  // 検証
  //--------------------------------------------------
  // 正準変量群
  Matrix Xc = A * X;
  Matrix Yc = B * Y;
  cout << "Xc\n" << Xc << endl;
  cout << "Yc\n" << Yc << endl;
  // 正準変量群間の相互相関は，正準相関係数が対角に並んだ行列
  cout << "Xc*Yc'/nsample\n" << Xc*Yc.transpose()/nsample << endl;

  return 0;
}
ページの先頭へ