liboctave で数値計算

〜 C++ で octave の行列演算機能を利用する 〜

概要
C++ で行列演算がメインの数値計算をしたいとき(ロボティクスとか機械学習の分野などで),いちいち自分で逆行列とかの基本的な演算を自分で実装してたら埓があかない.何かライブラリを利用したい.ここでは Matlab と互換性が(それなりに)ある octave を C++ から利用する方法, liboctave を解説する. liboctave は linux でも(一応)windows でも使えて,フリーでそこそこ信頼性もあるので,割と使いやすいのではないでしょうか.名前空間無いけどな.
目次

リンク

まずウェブ上で参考になる情報から. octave 関連では,

日本語マニュアルの質はあまりよくない.重要な情報が抜けている場合がある.なるべく本家のマニュアルを参照するようにした方がよい.

liboctave 関連の情報はあまり多くないが,

あたりが適当だろう. ちなみに,ほかの数値計算ライブラリにGSLIntel Math Kernel Library (MKL) などがあって,これらはCからでも利用できる.だが,複雑な行列計算を書くと,コードが繁雑になる.GSLにどのような機能があるかは,次を参照するとよい.

インストール

Debian

以下のパッケージをインストールすればよい.

 sudo apt-get -f install octave3.0 octave3.0-headers octave3.0-htmldoc
 sudo apt-get -f install atlas3-base libatlas3gf-base
 sudo apt-get -f install libfftw3-3 libfftw3-dev

バージョンは適宜変更すること.

The Other Linux

たいていのディストリビューションには octave 関連のバイナリパッケージが用意されている(はず).

Ubuntu 12.04 LTS 32 bit (12.07.06 update)

上記のDebianの項のようにインストールを進めてもよいが,インストール自体はできるものの,Matrixクラスも呼べない程度にしか構成されない(リンク等を編集するとできるんだろうが。。。).

現時点では,OctaveのHPからソースをダウンロードして,ビルドしインストールすることを推める.

下記リンクの「Ubuntu で, コマンドを手入力して,Octave のソースコードのダウンロードとビルドとインストールを行う」の項あたりから,必要なものを実行すると良い.

Windows

以下のリンク先を参照するとよい(と思う.試してない).

Mac

頑張れ.

liboctave の基礎:コンパイルと型と基本演算

liboctave には「行列型 Matrix 」や「列ベクトル型 ColumnVector 」がクラスとして実装されていて, これらの間の足し算や割算が operator として定義されている. だから,行列の複雑な演算がとてもシンプルに書くことができるのは,容易に想像がつくだろう.例えば,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A,B;
  ColumnVector x,y,z;
  // A,B,x,y,z に適当に値を代入
  z = B*(A*x+y);

みたいな感じだ.このように liboctave では,行列やベクトルの演算が,オブジェクト指向プログラミングによってとても直観的に実装できるように作られている.

以下では必要なヘッダとコンパイル方法, liboctave の基本型と基本的な演算について解説する.

ヘッダとコンパイルと実行

liboctave を利用するには,基本的には

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  #include <octave/config.h>
  #include <octave/Matrix.h>

を include すればよい.ひとつ具体例を示す

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  // hello-octave.cpp
  #include <octave/config.h>
  #include <octave/Matrix.h>
  #include <iostream>
  using namespace std;
  int main(int argc, char**argv)
  {
    Matrix A(2,3), B(2,2,1.0);
    ColumnVector x(3), y(2), z(2);
    for(int r(0); r<2; ++r)
      for(int c(0); c<3; ++c)
        A(r,c) = r+c;
    for(int r(0); r<3; ++r) x(r)=(r+1.0)*2.0;
    for(int r(0); r<2; ++r) y(r)=(r+1.0)*r;
    z = B*(A*x+y);
    cout<<"A="<<endl;
    cout<<A;
    cout<<"B="<<endl;
    cout<<B;
    cout<<"x= "<<x.transpose()<<endl;
    cout<<"y="<<y.transpose()<<endl;
    cout<<"z="<<z.transpose()<<endl;
    return 0;
  }

これをコンパイルするには,

 g++ -g -Wall -I/usr/include/octave-`octave-config -v` \
   -L/usr/lib/octave-`octave-config -v` -loctave -lcruft \
   hello-octave.cpp

のようにする.生成されたa.out を実行してみよう.ただし linux で, liboctave をシェアドライブラリとしてリンクしている場合は, LD_LIBRARY_PATH に /usr/lib/octave-\`octave-config -v\` を加えておく必要があるようだ.ちゃんと実行できたら,次のように結果が出力される.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
   0 1 2
   1 2 3
  B=
   1 1
   1 1
  x=  2 4 6
  y= 0 2
  z= 46 46

liboctave の基本型と基本演算

liboctave には以下のような型(クラス)が用意されている:

  • ColumnVector : double 型の列ベクトル(縦長のベクトル)
  • RowVector : double 型の行ベクトル(横長のベクトル)
  • Matrix : double 型の行列
  • Complex : 複素数型. std::complex<double> を typedef したものだ.
  • ComplexColumnVector : Complex 型の列ベクトル(縦長のベクトル)
  • ComplexRowVector : Complex 型の行ベクトル(横長のベクトル)
  • ComplexMatrix : Complex 型の行列

これ以外にも多数の型が定義されている.octave3.0-3.0.0: Class List に liboctave で定義されている全クラスが記載されているので,参照するとよい.

複素数版は double 版とほぼ同じ使い方なので,以下では double 版を中心に解説する.

宣言 (初期化)

列・行ベクトルの定義

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(int N);
  ColumnVector y(int N, double init);
  RowVector    u(int N);
  RowVector    v(int N, double init);

それぞれ,サイズ N,初期値 init の列・行ベクトルを生成する.初期値は指定しなくてもよい.

行列の定義

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(int r, int c);
  Matrix A(int r, int c, double init);

r 行 c 列,初期値 init の行列を生成する.初期値は指定しなくてもよい.

コピーコンストラクタ

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(3,2.5); // サイズ3, 各要素2.5の列ベクトル
  RowVector    y(x);  // x をコピー (サイズ3, 各要素2.5の「列」ベクトル)
  Matrix       A(x);  // x をコピー (サイズ(3,1), 各要素2.5の「行列」)
  Matrix       B(y);  // x をコピー (サイズ(1,3), 各要素2.5の「行列」)
  ColumnVector x2(x);  // x をコピー
  RowVector    y2(y);  // y をコピー
  Matrix       A2(A);  // A をコピー

入出力ストリーム演算

出力

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(3,2.5);
  RowVector    y(5,3.2);
  Matrix       A(2,3,0.5);
  cout<<"x="<<endl<<x;
  cout<<"y="<<endl<<y<<endl;
  cout<<"A="<<endl<<A;

このとき標準出力に,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  x=
  2.5
  2.5
  2.5
  y=
   3.2 3.2 3.2 3.2 3.2
  A=
   0.5 0.5 0.5
   0.5 0.5 0.5

と表示される.ファイルストリームなど,すべての std::ostream の派生クラスに出力できる.

入力

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(3);
  RowVector    y(5);
  Matrix       A(2,3);
  stringstream ss1("1 2 3"); ss1 >> x;
  stringstream ss2("5 4 3 2 1"); ss2 >> y;
  stringstream ss3("1 2 3 4 5 6"); ss3 >> A;

stringstream を使うためには, sstream を include する必要がある. stringstream の代わりに, cin でもよい.例からわかるように,行列やベクトルのサイズはあらかじめわかっている必要がある.利用頻度は低い.文字列解析が含まれるので,実行速度の低下を招くから.

要素アクセス

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x;
  RowVector    y;
  Matrix       A;
  x(r) = 3.0;  // 列ベクトル x の r 番目の要素
  y(c) = 3.0;  // 行ベクトル y の c 番目の要素
  A(r,c) = x(r)*y(c);
    // 行列 A の (r,c) 要素に x(r)*y(c) を代入

行列の要素は, A(行,列) の順. 各インデクスは0から始まることに注意 . C/C++の配列と同じである.

代入 (operator=)

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(2,1.0), y(5,3.0);
  RowVector z;
  x=y;
  z=y;
  Matrix A;
  // A=y;  <-- これはコンパイルエラー

サイズが異なっていても代入できる. ColumnVector を RowVector に代入することもできる(逆も可).ただし, Matrix に ColumnVector を代入したり, ColumnVector に Matrix を代入したりはできない.

演算 (operator+,-,*)

liboctave では,ほぼ 数学的な感覚通り に行列やベクトルの演算が実行できる.

かけ算:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(3); x(0)=1.0; x(1)=2.0; x(2)=3.0;
  RowVector    y(3); y(0)=3.0; y(1)=2.0; y(2)=1.0;
  Matrix       A(2,3);
  A(0,0)=1.0; A(0,1)=2.0; A(0,2)=3.0;
  A(1,0)=4.0; A(1,1)=3.0; A(1,2)=2.0;
  double       a(2.0);
  cout<<"x="<<endl<<x;
  cout<<"y="<<endl<<y<<endl;
  cout<<"A="<<endl<<A;
  cout<<"a="<<a<<endl;
  cout<<"A*x="<<endl<< A*x;
  cout<<"y*x= "<< y*x <<endl;
  cout<<"x*y= "<<endl<< x*y;
  cout<<"y*(x*y)= "<<endl<< y*(x*y)<<endl;
  cout<<"a*x="<<endl<<a*x;
  cout<<"a*y="<<endl<<a*y<<endl;
  cout<<"a*A="<<endl<<a*A;

結果(見やすくインデントしてある)

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  x=
    1
    2
    3
  y=
    3 2 1
  A=
    1 2 3
    4 3 2
  a=2
  A*x=
    14
    16
  y*x= 10
  x*y=
    3 2 1
    6 4 2
    9 6 3
  y*(x*y)=
    30 20 10
  a*x=
    2
    4
    6
  a*y=
    6 4 2
  a*A=
    2 4 6
    8 6 4

と表示される. 行ベクトルと列ベクトルの積が内積になっている ことも確認しておこう.

同様に足し算・引き算も直観的に使えるので,省略.

const メンバ関数

比較的よく使う const メンバ関数を紹介する(自身を変化させないメンバ関数). コメントは説明,後ろの括弧内は戻り値の型を表す.

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(4); x(0)=1.0; x(1)=2.0; x(2)=3.0; x(3)=4.0;
  ColumnVector y(2); y(0)=11.0; y(1)=12.0;
  Matrix       A(2,4);
  A(0,0)=1.0; A(0,1)=2.0; A(0,2)=3.0; A(0,3)=4.0;
  A(1,0)=5.0; A(1,1)=6.0; A(1,2)=7.0; A(1,3)=8.0;
  cout<<"x="<<endl<<x;
  cout<<"y="<<endl<<y;
  cout<<"A="<<endl<<A;
  // x のサイズ(int) :
  cout<<"x.length()= "<< x.length()<<endl;
  // x のサイズ(int, 非標準?) :
  cout<<"x.dim1()= "<< x.dim1()<<endl;
  // A の行数(int) :
  cout<<"A.rows()= "<< A.rows()<<endl;
  // A の列数(int) :
  cout<<"A.cols()= "<< A.cols()<<endl;
  // x の転置ベクトル(RowVector) :
  cout<<"x.transpose()="<<endl<< x.transpose()<<endl;
  // A の転置行列(Matrix) :
  cout<<"A.transpose()="<<endl<< A.transpose();
  // x の最小要素(double) :
  cout<<"x.min()= "<< x.min()<<endl;
  // x の最大要素(double) :
  cout<<"x.max()= "<< x.max()<<endl;
  // x の最後に y を付け足す(ColumnVector) :
  cout<<"x.stack(y)= "<<endl<< x.stack(y);
  // x(0)からx(2)までを抜き出す(ColumnVector) :
  cout<<"x.extract(0,2)= "<<endl<< x.extract(0,2);
  // x(1)から2個の要素を抜き出す(ColumnVector) :
  cout<<"x.extract_n(1,2)= "<<endl<< x.extract_n(1,2);
  // A の右端に列ベクトル y や行列 A を付け足す(Matrix) :
  cout<<"A.append(y)= "<<endl<< A.append(y);
  cout<<"A.append(A)= "<<endl<< A.append(A);
  // A の下端に行ベクトル x.transpose() や行列 A を付け足す(Matrix) :
  cout<<"A.stack(x.transpose())= "<<endl<< A.stack(x.transpose());
  cout<<"A.stack(A)= "<<endl<< A.stack(A);
  // A(0,1), A(1,2) で作られる四角形(行列)を抽出する(Matrix) :
  cout<<"A.extract(0,1,1,2)= "<<endl<< A.extract(0,1,1,2);
  // A(0,0) を基点にして, 2行3列分の要素(行列)を抽出する(Matrix) :
  cout<<"A.extract_n(0,0,2,3)= "<<endl<< A.extract_n(0,0,2,3);
  // おまけ :
  cout<<"A.extract_n(0,0,2,1).stack(y)= "<<endl<< A.extract_n(0,0,2,1).stack(y);
  // Aのすべての行の和を計算(Matrix) :
  cout<<"A.sum(0)= "<<endl<< A.sum(0);
  // Aのすべての列の和を計算(Matrix) :
  cout<<"A.sum(1)= "<<endl<< A.sum(1);
  // Aの2列目を抽出
  cout<<"A.column(2)= "<<endl<< A.column(2);
  // Aの1行目を抽出
  cout<<"A.row(1)= "<< A.row(1)<<endl;

結果(見やすくインデントしてある)

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  x=
    1
    2
    3
    4
  y=
    11
    12
  A=
    1 2 3 4
    5 6 7 8
  x.length()= 4
  x.dim1()= 4
  A.rows()= 2
  A.cols()= 4
  x.transpose()=
    1 2 3 4
  A.transpose()=
    1 5
    2 6
    3 7
    4 8
  x.min()= 1
  x.max()= 4
  x.stack(y)=
    1
    2
    3
    4
    11
    12
  x.extract(0,2)=
    1
    2
    3
  x.extract_n(1,2)=
    2
    3
  A.append(y)=
    1 2 3 4 11
    5 6 7 8 12
  A.append(A)=
    1 2 3 4 1 2 3 4
    5 6 7 8 5 6 7 8
  A.stack(x.transpose())=
    1 2 3 4
    5 6 7 8
    1 2 3 4
  A.stack(A)=
    1 2 3 4
    5 6 7 8
    1 2 3 4
    5 6 7 8
  A.extract(0,1,1,2)=
    2 3
    6 7
  A.extract_n(0,0,2,3)=
    1 2 3
    5 6 7
  A.extract_n(0,0,2,1).stack(y)=
    1
    5
    11
    12
  A.sum(0)=
    6 8 10 12
  A.sum(1)=
    10
    26
  A.column(2)=
    3
    7
  A.row(1)=  5 6 7 8

RowVector については省略したが,(たぶん)同じようなことができる.

非 const メンバ関数

非 const メンバ関数を紹介する(自身を変化させるメンバ関数)

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  ColumnVector x(4); x(0)=1.0; x(1)=2.0; x(2)=3.0; x(3)=4.0;
  ColumnVector y(2); y(0)=11.0; y(1)=12.0;
  Matrix       A(2,4);
  A(0,0)=1.0; A(0,1)=2.0; A(0,2)=3.0; A(0,3)=4.0;
  A(1,0)=5.0; A(1,1)=6.0; A(1,2)=7.0; A(1,3)=8.0;
  Matrix       B(2,2);
  B(0,0)=11.0; B(0,1)=12.0;
  B(1,0)=13.0; B(1,1)=14.0;
  cout<<"x="<<endl<<x;
  cout<<"y="<<endl<<y;
  cout<<"A="<<endl<<A;
  cout<<"B="<<endl<<B;
  // x(2) の位置に y を挿入する :
  cout<<"x.insert(y,2)="<<endl<< x.insert(y,2);
  // A(0,1) の位置に列ベクトル y を挿入する :
  cout<<"A.insert(y,0,1)="<<endl<< A.insert(y,0,1);
  // A(0,1) の位置に行ベクトル y.transpose() を挿入する :
  cout<<"A.insert(y.transpose(),0,1)="<<endl<< A.insert(y.transpose(),0,1);
  // A(0,1) の位置に行列 B を挿入する :
  cout<<"A.insert(B,0,1)="<<endl<< A.insert(B,0,1);
  // x(2) から x(3) までを 2.5 で埋める :
  cout<<"x.fill(2.5, 2,3)="<<endl<< x.fill(2.5, 2,3);
  // x を 4.5 で埋める :
  cout<<"x.fill(4.5)="<<endl<< x.fill(4.5);
  // A(0,1), A(1,2) で作られる四角形(行列)を 3.0 で埋める :
  cout<<"A.fill(3.0, 0,1, 1,2)="<<endl<< A.fill(3.0, 0,1, 1,2);
  // A を 5.0 で埋める :
  cout<<"A.fill(5.0)="<<endl<< A.fill(5.0);

結果(見やすくインデントしてある)

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  x=
    1
    2
    3
    4
  y=
    11
    12
  A=
    1 2 3 4
    5 6 7 8
  B=
    11 12
    13 14
  x.insert(y,2)=
    1
    2
    11
    12
  A.insert(y,0,1)=
    1 11 3 4
    5 12 7 8
  A.insert(y.transpose(),0,1)=
    1 11 12 4
    5 12 7 8
  A.insert(B,0,1)=
    1 11 12 4
    5 13 14 8
  x.fill(2.5, 2,3)=
    1
    2
    2.5
    2.5
  x.fill(4.5)=
    4.5
    4.5
    4.5
    4.5
  A.fill(3.0, 0,1, 1,2)=
    1 3 3 4
    5 3 3 8
  A.fill(5.0)=
    5 5 5 5
    5 5 5 5

「非const=自身を変化させる」ので, insert や fill によって x や A の値が変化していることに注意.

逆行列,疑似逆行列

行列(Matrix/ComplexMatrix) A が正則の場合には A.inverse() で逆行列が求められる.非正則の場合には A.pseudo_inverse() で疑似逆行列(ムーア・ペンローズ型の一般逆行列)が求められる. 数学的知識については以下のリンク/書籍を参照.

逆行列

逆行列の例

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(3,3);
  A(0,0)=1.0; A(0,1)=2.0;  A(0,2)=3.0;
  A(1,0)=2.0; A(1,1)=-1.0; A(1,2)=1.0;
  A(2,0)=4.0; A(2,1)=3.0;  A(2,2)=2.0;
  cout<<"A="<<endl<<A;
  cout<<"A.inverse()="<<endl<< A.inverse();
  cout<<"A * A.inverse()="<<endl<< A * A.inverse();

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    1 2 3
    2 -1 1
    4 3 2
  A.inverse()=
    -0.2 0.2 0.2
    -0 -0.4 0.2
    0.4 0.2 -0.2
  A * A.inverse()=
    1 1.11022e-16 -1.11022e-16
    0 1           -1.11022e-16
    0 0           1

逆行列計算時のエラー

逆行列の計算に失敗したときにエラーを出力させるには, int info; A.inverse(info) のようにする.ここで整数 info が 0 なら逆行列の計算に成功,それ以外(-1 になるようだ)だと失敗したと判断できる.例

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(3,3);
  A(0,0)=1.0; A(0,1)=2.0;  A(0,2)=3.0;
  A(1,0)=2.0; A(1,1)=-1.0; A(1,2)=1.0;
  A(2,0)=4.0; A(2,1)=3.0;  A(2,2)=2.0;
  int info;
  cout<<"A="<<endl<<A;
  cout<<"A.inverse(info)="<<endl<< A.inverse(info);
  cout<<"info= "<<info<<endl;
  cout<<"A.fill(0.0, 0,1, 2,1)="<<endl<< A.fill(0.0, 0,1, 2,1);
  cout<<"A.inverse(info)="<<endl<< A.inverse(info);
  cout<<"info= "<<info<<endl;

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    1 2 3
    2 -1 1
    4 3 2
  A.inverse(info)=
    -0.2 0.2 0.2
    -0 -0.4 0.2
    0.4 0.2 -0.2
  info= 0
  A.fill(0.0, 0,1, 2,1)=
    1 0 3
    2 0 1
    4 0 2
  A.inverse(info)=
    1 0 3
    2 0 1
    4 0 2
  info= -1

疑似逆行列

ムーア・ペンローズ型の一般逆行列もサクっと計算できるよ!

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(2,3);
  A(0,0)=1.0; A(0,1)=2.0;  A(0,2)=3.0;
  A(1,0)=2.0; A(1,1)=-1.0; A(1,2)=1.0;
  cout<<"A="<<endl<<A;
  cout<<"A.pseudo_inverse()="<<endl<< A.pseudo_inverse();
  cout<<"A * A.pseudo_inverse()="<<endl<< A * A.pseudo_inverse();

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    1 2 3
    2 -1 1
  A.pseudo_inverse()=
    -1.38778e-17 0.333333
    0.2 -0.266667
    0.2 0.0666667
  A * A.pseudo_inverse()=
    1 0
    0 1

行列式

行列式 (determinant) の計算は,例えば多次元のガウス関数を扱うときなどに必要だ.行列 (Matrix/ComplexMatrix) A に対して A.determinant().value() で行列式が計算できる.行列式計算中に発生したエラーの判定などは,逆行列の場合と同じようにして判定できる

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(4,4);
  stringstream ss (
    "  1.0  2.0  3.0  4.0"
    "  2.0 -1.0  1.0  0.0"
    "  4.0  3.0  2.0  1.0"
    " -4.0 -1.0 -5.0  0.0"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // 行列式を得る(double) :
  cout<<"A.determinant().value()= "<< A.determinant().value()<<endl;
  // 行列式を計算し,エラー情報をinfoに格納 :
  cout<<"A.determinant(info).value()= "<< A.determinant(info).value()<<endl;
  cout<<"info= "<<info<<endl;
  cout<<"A.fill(0.0, 2,0, 2,3)="<<endl<< A.fill(0.0, 2,0, 2,3);
  // 行列式の計算に失敗する場合 :
  cout<<"A.determinant(info).value()= "<< A.determinant(info).value()<<endl;
  cout<<"info= "<<info<<endl;

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    1 2 3 4
    2 -1 1 0
    4 3 2 1
    -4 -1 -5 0
  A.determinant().value()= 120
  A.determinant(info).value()= 120
    info= 0
  A.fill(0.0, 2,0, 2,3)=
    1 2 3 4
    2 -1 1 0
    0 0 0 0
    -4 -1 -5 0
  A.determinant(info).value()= 0
    info= -1

この行列式を返すメンバ関数の宣言は

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
 DET Matrix::determinant (octave_idx_type &info) const;

となっていて,DET クラス を返すことがわかる. DET クラスの主なメンバ関数には以下のものがある

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(4,4);
  stringstream ss (
    "  4.0  2.0  3.0  1.0"
    "  0.0 -1.0  1.0  2.0"
    "  1.0  3.0  2.0  4.0"
    "  0.0 -1.0 -5.0 -4.0"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // DET オブジェクトの生成
  DET det(A.determinant(info));
  // 行列式の値を得る(double) :
  cout<<"det.value()= "<< det.value()<<endl;
  // 係数を得る(int) :
  cout<<"det.coefficient()= "<< det.coefficient()<<endl;
  // 指数を得る(int) :
  cout<<"det.exponent()= "<< det.exponent()<<endl;

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    4 2 3 1
    0 -1 1 2
    1 3 2 4
    0 -1 -5 -4
  det.value()= -120
  det.coefficient()= -1.2
  det.exponent()= 2

行列分解(固有値分解,LU分解,..)

固有値分解などの各種行列分解の方法について解説する.

固有値分解

固有値分解(Eigendecomposition) の用途: 主成分分析など多数. 正方行列 に対してのみ適用可能(非正方の場合は特異値分解).

行列 (Matrix) A に対する固有値分解は,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  EIG  eig (const Matrix& A, bool calc_eigenvectors=true);
  EIG  eig (const ComplexMatrix& A, bool calc_eigenvectors=true);

によって得られる.ここでEIG クラス のオブジェクト eig には固有値分解の結果,すなわち固有値 (eigenvalues) と固有ベクトル (eigenvectors) が格納される. calc_eigenvectors は固有ベクトルを計算するかどうかを表すフラグである(省略してもよい).

EIG オブジェクトは以下の public メンバ関数を持つ

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  // 固有値を返す(ComplexColumnVector) :
  ComplexColumnVector  eigenvalues (void) const;
  // 固有ベクトルを返す(ComplexMatrix) :
  ComplexMatrix  eigenvectors (void) const;

個々の固有値には eig.eigenvalues()(i), i:int のようにしてアクセスできる(iの範囲はAの行数と同じ). eig.eigenvalues() は固有値を縦に並べた列ベクトル (ComplexColumnVector) である.個々の固有ベクトルには eig.eigenvectors().column(i), i:int のようにしてアクセスできる. eig.eigenvectors() は固有ベクトル(列ベクトル)を横に並べた行列 (ComplexMatrix) である.

固有値・固有ベクトルのいずれも,要素の一部または全部が複素数になる場合がある.実部を取り出す必要がある場合は, real 関数を使う. real(eig.eigenvalues()) とすると実数の ColumnVector , real(eig.eigenvalues()(i)) とすると実数 (double), real(eig.eigenvectors()) とすると実数の Matrix , real(eig.eigenvectors().column(i)) とすると実数の ColumnVector がそれぞれ得られる.

注意! octave のマニュアルにあるように, The eigenvalues returned by eig are not ordered. ,つまり eig.eigenvalues() はソートされていない.サンプルデータの共分散行列に EIG を適用すると昇順にソートされているように見えるが,保証はない.

主成分分析などで固有値の大きいものから順に取得したい場合, eig.eigenvalues() をソートすればよいのだが,固有ベクトルのソートも同時に行わなければならない.この対処のひとつとして, eig.eigenvalues() を直接ソートせずに, {std::vector<int> idx, idx[n] はn番目に大きい固有値のインデックス} となるようなインデックスベクトル idx を定義する方法が考えられる.具体的な実装は主成分分析のサンプルプログラムを参照.ただし,主成分分析をする場合は特異値分解を用いることが多いようで(精度がいいらしい),この場合, liboctave では固有値と固有ベクトルはソートされたものが得られるから,特に気にする必要はない.

サンプル

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(4,4);
  stringstream ss (
    "  4.0  2.0  3.0  1.0"
    "  0.0 -1.0  1.0  2.0"
    "  1.0  3.0  2.0  4.0"
    "  0.0 -1.0 -5.0 -4.0"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // EIG オブジェクトの生成 :
  EIG  eig (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // A を対称にする :
  A= 0.5*(A+A.transpose());
  cout<<"A="<<endl<<A;
  // 再度固有値分解 :
  eig= EIG(A,info);
  cout<<"info= "<<info<<endl;
  cout<<"eig.eigenvalues()= "<<endl<< eig.eigenvalues();
  cout<<"eig.eigenvectors()= "<<endl<< eig.eigenvectors();
  // 固有値,固有ベクトルを実数値で抽出する :
  cout<<"real(eig.eigenvalues()(2))= "<< real(eig.eigenvalues()(2))<<endl;
  cout<<"real(eig.eigenvectors().column(2))= "<<endl<< real(eig.eigenvectors().column(2));

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    4 2 3 1
    0 -1 1 2
    1 3 2 4
    0 -1 -5 -4
  info= 0
  eig.eigenvalues()=
    (4.4782,0)
    (-0.524237,3.27929)
    (-0.524237,-3.27929)
    (-2.42973,0)
  eig.eigenvectors()=
    (0.9733,0) (0.0591816,-0.335765) (0.0591816,0.335765) (-0.168139,0)
    (-0.00621851,0) (0.0758315,0.282186) (0.0758315,-0.282186) (0.796196,0)
    (0.197962,0) (0.450522,0.382927) (0.450522,-0.382927) (0.0231463,0)
    (-0.116014,0) (-0.669908,0) (-0.669908,-0) (-0.580745,0)
  A=
    4 1 2 0.5
    1 -1 2 0.5
    2 2 2 -0.5
    0.5 0.5 -0.5 -4
  info= 0
  eig.eigenvalues()=
    (-4.24974,0)
    (-1.80298,0)
    (1.25496,0)
    (5.79776,0)
  eig.eigenvectors()=
    (0.0718392,0) (-0.0357097,0) (0.618216,0) (-0.781904,0)
    (0.2315,0) (0.870549,0) (-0.331435,0) (-0.280539,0)
    (-0.17344,0) (-0.402036,0) (-0.70642,0) (-0.556108,0)
    (-0.95455,0) (0.28149,0) (0.0945012,0) (-0.0258393,0)
  real(eig.eigenvalues()(2))= 1.25496
  real(eig.eigenvectors().column(2))=
    0.618216
    -0.331435
    -0.70642
    0.0945012

eig.eigenvalues() などの結果が (r,i) のように書かれているのは,複素数だからだ. r が実部で, i が虚部である.

特異値分解

特異値分解(Singular value decomposition, SVD) の用途: 主成分分析,疑似逆行列の計算,機械学習,パタン認識,数値解析など.任意の形の行列を分解できる.

SingularValue Decomposition -- LAPACK Users Guide

行列 (Matrix/ComplexMatrix) A に対する特異値分解は

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
SVD  svd (const Matrix& A, SVD::type svd_type=SVD::std);
SVD  svd (const Matrix& A, int& info, SVD::type svd_type=SVD::std);
ComplexSVD  svd (const ComplexMatrix& A, SVD::type svd_type=SVD::std);
ComplexSVD  svd (const ComplexMatrix& A, int& info, SVD::type svd_type=SVD::std);

によって得られる. octave では,行列Aは

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A = USV*

のように分解される.ここでU, V はユニタリ行列, V* は V の共役転置を表す. S は対角成分が A の特異値(対角成分以外はゼロ)である. S は大きい順にソートされている (この情報は octave のマニュアルに基づくものではなく, Singular Value Decomposition -- LAPACK Users' Guide に記述されていることからの判断である).

ここでSVD クラス, ComplexSVD クラス のオブジェクトには特異値分解の結果が代入される. info は計算後にゼロ以外なら何らかのエラーが発生したことを示すフラグである. SVD::type は列挙型で,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  00032 class
  00033 OCTAVE_API
  00034 SVD
  00035 {
  00036 public:
  00038   enum type
  00039     {
  00040       std,
  00041       economy,
  00042       sigma_only
  00043     };

のように定義されている.デフォルト値SVD::std は通常の特異値分解(U,S,Vすべて計算), SVD::sigmaonly は S のみ計算, SVD::economy はエコノミーサイズ(UとVの不要な列や行を消去)で計算する.

SVD オブジェクトは以下の public メンバ関数を持つ

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  // S を返す :
  DiagMatrix   singular_values (void) const;
  // U を返す :
  Matrix   left_singular_matrix (void) const;
  // V を返す :
  Matrix   right_singular_matrix (void) const;

サンプル

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(3,6);
  stringstream ss (
    "  4   0   4   6   5   4"
    "  9   6   2   4   1   0"
    "  1   6   9   0   7   5"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // SVD オブジェクトの生成 :
  SVD  svd (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"svd.singular_values()= "<<endl<< svd.singular_values();
  cout<<"svd.left_singular_matrix()= "<<endl<< svd.left_singular_matrix();
  cout<<"svd.right_singular_matrix()= "<<endl<< svd.right_singular_matrix();
  // USV* を計算(Aに一致するか?) :
  cout<<"svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()= "<<endl
  << svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose();
  // エコノミーモードで特異値分解 :
  svd= SVD(A, info, SVD::economy);
  cout<<"info= "<<info<<endl;
  cout<<"svd.singular_values()= "<<endl<< svd.singular_values();
  cout<<"svd.left_singular_matrix()= "<<endl<< svd.left_singular_matrix();
  cout<<"svd.right_singular_matrix()= "<<endl<< svd.right_singular_matrix();
  // USV* を計算(Aに一致するか?) :
  cout<<"svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()= "<<endl
  << svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose();

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    4 0 4 6 5 4
    9 6 2 4 1 0
    1 6 9 0 7 5
  info= 0
  svd.singular_values()=
    17.6341 0 0 0 0 0
    0 9.51377 0 0 0 0
    0 0 6.12601 0 0 0
  svd.left_singular_matrix()=
    -0.509954 0.0830227 0.856186
    -0.498851 0.782325 -0.372982
    -0.700782 -0.617312 -0.357534
  svd.right_singular_matrix()=
    -0.410016 0.710098 -0.0472778 -0.556602 -0.109191 -0.060699
    -0.408176 0.104068 -0.715489 0.416853 0.235993 0.284919
    -0.529915 -0.384608 -0.0879892 0.0729675 -0.59143 -0.45652
    -0.286668 0.381283 0.595035 0.646778 0.00468758 -0.00584839
    -0.451064 -0.32834 0.229385 -0.224744 0.730871 -0.22679
    -0.314376 -0.289525 0.267234 -0.205602 -0.219999 0.80948
  svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()=
    4 -8.88178e-16 4 6 5 4
    9 6 2 4 1 -1.77636e-15
    1 6 9 2.22045e-16 7 5
  info= 0
  svd.singular_values()=
    17.6341 0 0
    0 9.51377 0
    0 0 6.12601
  svd.left_singular_matrix()=
    -0.509954 0.0830227 0.856186
    -0.498851 0.782325 -0.372982
    -0.700782 -0.617312 -0.357534
  svd.right_singular_matrix()=
    -0.410016 0.710098 -0.0472778
    -0.408176 0.104068 -0.715489
    -0.529915 -0.384608 -0.0879892
    -0.286668 0.381283 0.595035
    -0.451064 -0.32834 0.229385
    -0.314376 -0.289525 0.267234
  svd.left_singular_matrix()*svd.singular_values()*svd.right_singular_matrix().transpose()=
    4 -8.88178e-16 4 6 5 4
    9 6 2 4 1 -1.77636e-15
    1 6 9 2.22045e-16 7 5

コレスキー分解

コレスキー分解(Cholesky decomposition) の用途: サンプルデータの前処理, Ax = b の数値解,モンテカルロ法,カルマンフィルタなど(一部 Cholesky decomposition - Wikipedia, the free encyclopedia より). 正定値エルミート行列 (実数値行列なら, 正定値対称行列 ) にのみ適用可能.

行列 (Matrix/ComplexMatrix) A に対するコレスキー分解は

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  CHOL  chol (const Matrix& A);
  CHOL  chol (const Matrix& A, octave_idx_type& info);
  ComplexCHOL  chol (const ComplexMatrix& A);
  ComplexCHOL  chol (const ComplexMatrix& A, octave_idx_type& info);

によって得られる.ここでCHOL クラス, ComplexCHOL クラス のオブジェクトにはコレスキー分解の結果の行列が代入されている. info は計算後にゼロ以外なら何らかのエラーが発生したことを示すフラグである.

CHOL オブジェクトは以下の public メンバ関数を持つ:

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  // コレスキー分解された行列を返す :
  Matrix  chol_matrix (void) const;

備考

liboctave のコレスキー分解は, A = chol.chol_matrix().transpose() * chol.chol_matrix() となるように計算される.

サンプル

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(4,4);
  stringstream ss (
    "  4.0  2.0  3.0  1.0"
    "  0.0 -1.0  1.0  2.0"
    "  1.0  3.0  2.0  4.0"
    "  0.0 -1.0 -5.0 -4.0"); ss >> A;
  int info;
  cout<<"A="<<endl<<A;
  // CHOL オブジェクトの生成 :
  CHOL  chol (A, info);
  cout<<"info= "<<info<<endl;
  cout<<"chol.chol_matrix()= "<<endl<< chol.chol_matrix();
  cout<<"chol.chol_matrix().transpose()*chol.chol_matrix()= "<<endl
  << chol.chol_matrix().transpose()*chol.chol_matrix();
  // A を正定値対称にする :
  A= A.transpose()*A;
  cout<<"A="<<endl<<A;
  // 再度コレスキー分解 :
  chol= CHOL(A,info);
  cout<<"info= "<<info<<endl;
  cout<<"chol.chol_matrix()= "<<endl<< chol.chol_matrix();
  cout<<"chol.chol_matrix().transpose()*chol.chol_matrix()= "<<endl
  << chol.chol_matrix().transpose()*chol.chol_matrix();

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    4 2 3 1
    0 -1 1 2
    1 3 2 4
    0 -1 -5 -4
  info= -1
  chol.chol_matrix()=
    2 1 1.5 0.5
    0 -2 1 2
    1 3 2 4
    0 -1 -5 -4
  chol.chol_matrix().transpose()*chol.chol_matrix()=
    5 5 5 5
    5 15 10.5 12.5
    5 10.5 32.25 30.75
    5 12.5 30.75 36.25
  A=
    17 11 14 8
    11 15 16 16
    14 16 39 33
    8 16 33 37
  info= 0
  chol.chol_matrix()=
    4.12311 2.66789 3.3955 1.94029
    0 2.80755 2.47232 3.85515
    0 0 4.62149 3.65263
    0 0 0 2.24309
  chol.chol_matrix().transpose()*chol.chol_matrix()=
    17 11 14 8
    11 15 16 16
    14 16 39 33
    8 16 33 37

前半の行列は対称でも正定値でもないため,コレスキー分解に失敗しているが,後半の行列は正定値対称行列なので,コレスキー分解に成功している.

LU分解

LU分解(LU decomposition) の用途: 連立方程式を解くなど.正方行列でなくともよい.

行列 (Matrix/ComplexMatrix) A に対するLU分解は

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  LU  fact (const Matrix& A);
  ComplexLU  fact (const ComplexMatrix& A);

によって得られる. octave では,行列Aは

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  PA = LU

のように分解される.ここでL, U はLU分解された行列, P は置換行列 Permutation matrix) である. P は直行行列なので,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A = P'LU

のようにして A を復元できる( P' は P の転置行列).

余談だが, octave で

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  [l,u]=lu(A)

とした場合には, l=P'L, u=U が代入され,

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  [l,u,p]=lu(A)

とした場合には, l=L, u=U, p=P が代入される.

ここでLU クラス, ComplexLU クラス のオブジェクトにはLU分解の結果が代入される.

LU オブジェクトは以下の public メンバ関数を持つ

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  // L を返す(Matrix) :
  Matrix  L (void) const;
  // U を返す(Matrix) :
  Matrix  U (void) const;
  // P を返す(Matrix) :
  Matrix  P (void) const;

サンプル

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  Matrix A(3,3);
  stringstream ss (
    "  3 1 2"
    "  5 0 0"
    "  4 8 1"); ss >> A;
  cout<<"A="<<endl<<A;
  // LU オブジェクトの生成 :
  LU  fact (A);
  cout<<"fact.L()= "<<endl<< fact.L();
  cout<<"fact.U()= "<<endl<< fact.U();
  cout<<"fact.P()= "<<endl<< fact.P();
  // P'LU を計算(Aに一致するか?) :
  cout<<"fact.P().transpose()*fact.L()*fact.U()= "<<endl
    << fact.P().transpose()*fact.L()*fact.U();
  // 行列Aを変更
  A.resize(3,6);
  ss.clear(); ss.str(
    "  4   0   4   6   5   4"
    "  9   6   2   4   1   0"
    "  1   6   9   0   7   5"); ss >> A;
  cout<<"A="<<endl<<A;
  // 再度LU分解 :
  fact= LU(A);
  cout<<"fact.L()= "<<endl<< fact.L();
  cout<<"fact.U()= "<<endl<< fact.U();
  cout<<"fact.P()= "<<endl<< fact.P();
  // P'LU を計算(Aに一致するか?) :
  cout<<"fact.P().transpose()*fact.L()*fact.U()= "<<endl
    << fact.P().transpose()*fact.L()*fact.U();

結果

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  A=
    3 1 2
    5 0 0
    4 8 1
  fact.L()=
    1 0 0
    0.8 1 0
    0.6 0.125 1
  fact.U()=
    5 0 0
    0 8 1
    0 0 1.875
  fact.P()=
    0 1 0
    0 0 1
    1 0 0
  fact.P().transpose()*fact.L()*fact.U()=
    3 1 2
    5 0 0
    4 8 1
  A=
    4 0 4 6 5 4
    9 6 2 4 1 0
    1 6 9 0 7 5
  fact.L()=
    1 0 0
    0.111111 1 0
    0.444444 -0.5 1
  fact.U()=
    9 6 2 4 1 0
    0 5.33333 8.77778 -0.444444 6.88889 5
    0 0 7.5 4 8 6.5
  fact.P()=
    0 1 0
    0 0 1
    1 0 0
  fact.P().transpose()*fact.L()*fact.U()=
    4 0 4 6 5 4
    9 6 2 4 1 0
    1 6 9 0 7 5

FFT

Matrix::fourier, ComplexMatrix::fourier を呼び出せばおしまい....なのか??

特殊関数

高速化のためのコンパイルオプション

C++ で書いたコードを octave から呼び出す

liboctave+C++ で可能なことを増やすために

octave は Matlab には劣るかもしれないが,かなり多くのことができる.それらのうち, liboctave で実装されているものはもちろん C++ から利用できるのだが,そうでないもの,例えばスクリプトやDLDなどによって提供されているものについても, C++ から利用できる場合がある.そのための力技テクニックを紹介する.

liboctave に関するドキュメントはとても少ない.よって liboctave で可能なことを増やすためには,どうしてもコードを読み,場合によっては改変する必要がある.

情報源

octaveマニュアルを C++ のために読む

octaveは liboctave+C++ で書かれているから,octaveマニュアルに書かれていることは,すべて C++ 上で実現できる. octave で利用できる関数が, C++ ではどのように実装されているかを探す方法を述べておこう.

octaveがデフォルトで使える関数は次の4つのうちいずれかの形式で提供されており,関数の説明箇所でどの形式かが明記されている.

  • [Function File] The function described is defned using Octave commands stored in a text file. See Section 11.7 [Function Files], page 113.
  • [Built-in Function] The function described is written in a language like C++, C, or Fortran, and is part of the compiled Octave binary.
  • [Loadable Function] The function described is written in a language like C++, C, or Fortran. On systems that support dynamic linking of user-supplied functions, it may be automatically linked while Octave is running, but only if it is needed. See Appendix A [Dynami-cally Linked Functions], page 441.
  • [Mapping Function] The function described works element-by-element for matrix and vector arguments. [octaveマニュアル p.21 から抜粋]

この情報をもとにしてC++で使いたい関数を調べ,利用する. octave のソースコードを調べるので, octave のソースをどこかに展開しておく必要がある.以下では octave-src-dir は, octave のソースディレクトリを指すものとする.

[Function File] の利用

基本的に, octave-src-dir/scripts に"関数名.m"というファイル名で置かれている.このスクリプトを liboctave+C++ に翻訳する という力技で,C++から利用できるようになる. octave のスクリプトを C++ のソースに自動翻訳するプログラムも存在するようだが,信頼性は低いようだ.

[Built-in Function] の利用

基本的に`octave-src-dir/src/data.cc` で定義されている.例えば columns というビルトイン関数は, data.cc で

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  01302 DEFUN (columns, args, ,
  01303   "-*- texinfo -*-\n\
  01304 @deftypefn {Built-in Function} {} columns (@var{a})\n\
  01305 Return the number of columns of @var{a}.\n\
  01306 @seealso{size, numel, rows, length, isscalar, isvector, and ismatrix}\n\
  01307 @end deftypefn")
  01308 {
  01309   octave_value retval;
  01310
  01311   if (args.length () == 1)
  01312     retval = args(0).columns ();
  01313   else
  01314     print_usage ();
  01315
  01316   return retval;
  01317 }

のように定義されている.前半6行はoctave に関数を提供するためのマクロ,残りが関数の中身である. octave_value retval は関数の戻り値で, args(0) は関数の最初の引数(行列)である. retval = args(0).columns (); で,行列の列数を retval に代入していることがわかる.

このように, data.cc からビルトイン関数の定義場所を探して,それを C++ 向けに少し書き直せば,比較的容易に C++ からビルトイン関数が利用できるようになる.

[Loadable Function] の利用

これは「C++で記述した関数を odet.value()= ctave から利endl用する」方法によって提供されている関数である.この提供方法はビルトイン関数とほとんど同じで,既に C++ で実装されているため,容易に書き直せるだろう.これらの関数は octave-src-dir/src/DLD-FUNCTIONS に置かれていて,多くの関数が 関数名.cc で定義されている.

[Mapping Function] の利用

多くの関数がoctave-src-dir/src/mappers.cc で定義されているようだ.例えば octave の isspace は mappers.cc で

length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  00103 static int
  00104 xisspace (int c)
  00105 {
  00106   return isspace (c);
  00107 }
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1360.
length() used on @lines (did you mean "scalar(@lines)"?) at /usr/bin/code2html line 1370.
  00511   DEFUN_MAPPER (isspace, xisspace, 0, 0, 0, 0, 0, 0.0, 0.0, 0, 0,
  00512     "-*- texinfo -*-\n\
  00513 @deftypa*x=efn {Mapping Function} {} isspace (@var{s})\n\
  00514 Return 1 for whitespace characters (space, formfeed, newline,\n\
  00515 carriage return, tab, and vertical tab).\n\
  00516 @end deftypefn");

のように定義されている.前半は xisspace というC++の関数を定義しており,後半は octave の isspace を C++ の xisspace にマップしていると考えられる.

見つからない場合

ほとんどの関数のC++による実装が上記の探し方で見つかるはずだが,それでも見つからない場合は,利用したい関数名を octave のソースディレクトリで grep コマンドなどによって調べよう.

リンク

数値計算関係のリンク.

MacWiki - 数値計算:


トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2018-08-30 (木) 07:17:04 (19d)