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

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

  Matrix A,B;
  ColumnVector x,y,z;
  // A,B,x,y,z に適当に値を代入
  z = B*(A*x+y);

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

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

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

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

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

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

  // 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\` を加えておく必要があるようだ.ちゃんと実行できたら,次のように結果が出力される.

  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 を typedef したものだ.
  • ComplexColumnVector : Complex 型の列ベクトル(縦長のベクトル)
  • ComplexRowVector : Complex 型の行ベクトル(横長のベクトル)
  • ComplexMatrix : Complex 型の行列

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

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

宣言 (初期化)

列・行ベクトルの定義:

  ColumnVector x(int N);
  ColumnVector y(int N, double init);
  RowVector    u(int N);
  RowVector    v(int N, double init);

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

行列の定義:

  Matrix A(int r, int c);
  Matrix A(int r, int c, double init);

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

コピーコンストラクタ:

  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 をコピー
入出力ストリーム演算

出力:

  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;

このとき標準出力に,

  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 の派生クラスに出力できる.

入力:

  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 でもよい.例からわかるように,行列やベクトルのサイズはあらかじめわかっている必要がある.利用頻度は低い.文字列解析が含まれるので,実行速度の低下を招くから.

要素アクセス
  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=)
  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 では,ほぼ 数学的な感覚通り に行列やベクトルの演算が実行できる.

かけ算:

  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;

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

  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 メンバ関数を紹介する(自身を変化させないメンバ関数).
コメントは説明,後ろの括弧内は戻り値の型を表す.

  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;

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

  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 メンバ関数を紹介する(自身を変化させるメンバ関数):

  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);

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

  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 の値が変化していることに注意.