Boost 数学関係ライブラリの使い方

Last update : Jan. 13th, 2005

この文書は，線形演算ライブラリ boost::numeric::ublas の使い方の一部を簡単に説明したものです．

どうも boost ― uBLAS については日本語の説明書きがないようです．頼みの日本語解説書^[2]も uBLAS はたった 2 ページ．Web をあさっても私の希望にあう解説は見あたりません．仕方がないので英語のオリジナルドキュメントと格闘しました．その結果，なんとか連立1次方程式を解くことと，逆行列を求めることはできるようになったので，私と同じようなお悩みを抱えて Web を巡回中のあなたに私のメモ書きを公開することにしました．

ここで説明する ublas は，Boost C++ Library の一部です．Boost は膨大な数の C++ ライブラリの集合体です．非常にパワフルかつフリーです．あなたもぜひ Boost にハマってください（笑）．日本語の解説書も出ています（こちらを参照）．全てのライブラリを網羅とまではいきませんが，あなたの Boost ライフをハッピーにすること間違いなしです．内容の割にお値段も 2,800 円と手頃なため，是非手元に置いておくことをおすすめします．

メリットは大きく三つあります．

第一に，計算が高速らしいことです．これは expression template^（下注） という C++ の高度なテクニックによるものです^[2]p.353．

第二に，様々なタイプの行列やベクトルに適したクラスが定義されていることです．行列のほとんどの要素がゼロである行列は「疎行列」と呼ばれますが，このような行列の場合は値が 0 でない要素の分だけメモリを確保すれば十分です．また，対称行列やエルミート行列では，( i , j ) 要素の値は ( j , i ) 要素を見ることで簡単に分かりますから，まともに全要素分のメモリを確保するとその約半分は無駄になります．boost::numeric::ublas では，これらの行列に対する適切な内部表現を持ったクラスが提供されています．

第三に，ライセンスが極めて緩いものになっていることです．Boost は誰でも自由に無料で使用できます．Boost を使用して制作したソフトウェアも全く自由に配布でき，GPL のようにソースコードの公開義務もありません．バイナリを配布するときは著作権表示すら必要ないという自由さです．このようにもともと広く使われる素地を持っていることから，将来的にデファクトスタンダードとなることが期待されます．

expression template 技法

偉そうに注を付けていますが，実は私はこのテクニックについて何も知りません．文献[2]によれば，expression template とは，「式の情報ををテンプレートの構造として保持する」ような技法らしいです．

筆者が理解している範囲で具体例を示しますが，あくまで私の理解を文章にしただけなので，間違っている可能性があります．ご了解の上お読みください．

例えば行列 A と B があって，prod(A,B)+Aと書いたとしても，実は A・B + A はすぐには計算されず，この式は「A と B を掛けて A を足した」という情報を示す型のオブジェクトを返します．ここで，別な行列オブジェクト C を用意し，C = prod(A,B)+Aとすると，代入演算子が実行される段階で初めて掛け算と足し算が発生するそうです．それも，「まず掛け算をやって，その結果に足し算をする」のではなく，同時に処理してしまうらしいです．

つまり，実際の計算を本当にその値が必要になるときまで遅らせることができるということです．また，演算のたびに発生する一時オブジェクトが不要になることもあり，高速化が期待できるとされています^{（筆者注）}．これ以上の詳しい説明は私の脳味噌には期待しないでください（笑）．詳しく知りたい方はググるなり本（日本語の本はまだないそうですが）を買うなりして調べてください．

（筆者注） これは行列演算に限っては逆に遅くなる可能性もあるような…

既に Boost がインストール済みであることを前提とします．Boost のインストール方法はこちら^[3]．

この文書の記述は，Boost バージョン 1.32.0 に基づいています．

この説明する全てのクラスは，特記以外 boost::numeric::ublas 名前空間で定義されています．また，ヘッダファイルをインクルードする必要がある場合は，その旨記載します．

あなたは、金星や火星に高くジャンプだろうか？

Borland C++ Compiler 5.6.4 では一部機能が動作しないそうです．何が使えて何が使えないのかは私には分かりません．gcc も 2.95.3 は問題があるそうなので 3.3.3 以降を使いましょう．Visual C++ 6.0, 同 .net, 同 .net 2003 では問題は出ないようです^[1][2]．本文書に掲載のサンプルコードは Visual C++ .net 2003 でコンパイルと実行確認を行っています．

Boost 1.32.0 インストール時の注意

こちらを参照してください．

行列・ベクトル変数の宣言

行列やベクトルの変数を宣言するだけのサンプルコードを示します．基本的に，テンプレート引数に要素の型（float, double, int, std::complex など）を指定し，インスタンス名の後に行数と列数を指定して宣言します．

#include <boost/numeric/ublas/matrix.hpp> // (1) 普通の行列用のヘッダ #include <boost/numeric/ublas/triangular.hpp> // (2) 三角行列用のヘッダ #include <boost/numeric/ublas/symmetric.hpp> // (3) 対称行列用のヘッダ #include <boost/numeric/ublas/hermitian.hpp> // (4) エルミート行列用のヘッダ #include <boost/numeric/ublas/matrix_sparse.hpp> // (5) 疎行列用のヘッダ #include <boost/numeric/ublas/vector.hpp> // (6) ベクトル用のヘッダ #include <boost/numeric/ublas/vector_sparse.hpp> // (7) 疎ベクトル用のヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 using namespace boost::numeric; // boost::numeric 名前空間を使用(※下注) // (※注) // uBLAS の vector クラスは STL の vector クラスと名前が衝突しています． // このため，使用にあたっては名前空間を明示し，ublas::vector<double> ... // としなければなりません．もちろん std::vector を使わないのであれば // 問題はないのですが． int main() { // 行列クラス名<要素の型> インスタンス名(行数,列数) で宣言 // ただし，正方行列であることが自明な場合は行数と列数を // 指定するのではなく，単に次数を指定する． // 後ろのコメントに (番号) を付したのものは対応するヘッダを // インクルードする必要がある． matrix<double> mat_d(3,4); // (1) double を要素とする通常の 3 行 4 列の行列 matrix<float> mat_f(2,3); // (1) float を要素とする通常の 2 行 3 列の行列 triangular_matrix<double> mat_t(3,3); // (2) double を要素とする 3 行 3 列の三角行列 symmetric_matrix<double> mat_sm(3); // (3) double を要素とする 3 次対称行列 hermitian_matrix<double> mat_he(3); // (4) double を要素とする 3 次エルミート行列 sparse_matrix<int> mat_sp(10,1000); // (5) int を要素とする 10 行 1000 列の疎行列 vector<float> vec_f(4); // (6) float を要素とする通常の 4 次ベクトル sparse_vector<double> vec_sp(100); // (7) double を要素とする 100 次疎ベクトル // よく使う行列のタイプは typedef しておくとよい． typedef matrix<double> dmatrix; dmatrix mat_tdef(3,3); // typedef により，matrix<double> mat_tdef(3,3) と同じ意味 return 0; }

上記の例以外に，対角要素に近い要素だけが非零であるような「帯行列」を扱うbanded_matrixクラステンプレートが用意されています．

要素の参照と代入

行列の ( i , j ) 要素を参照するには，mat(i,j)とするか，通常の 2 次元配列と同様に mat[i][j]とします．サンプルコードを以下に示します．i , j はいずれも 0 から始まるインデックスです．

ベクトルの i 番目の要素を参照するには，vec(i)とするか，通常の 1 次元配列と同様に vec[i]とします．i はもちろん 0 から始まるインデックスです．

#include <boost/numeric/ublas/matrix.hpp> // 普通の行列用のヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 #include <iostream> using namespace std; int main() { typedef matrix<double> dmatrix; dmatrix mat1(2,2); // (i,j)も[i][j]も OK!! mat1(0,0) = 1.0; mat1(0,1) = 2.0; mat1[1][0] = 3.0; mat1[1][1] = 4.0; return 0; }

ストリームへの出力

<<演算子を使って，標準出力やファイルストリームへ行列やベクトルを出力することができます．ただし，ヘッダファイル のインクルードが必要です．

（注） を行列クラスやベクトルクラスが定義されているヘッダよりも先にインクルードしてはいけません．これをやってしまうとコンパイルエラーになります．VC++ .net 2003 では 「<< のオーバーロードが曖昧です」と怒られました．

47の2連続した整数を見つけるためにどのように

#include <boost/numeric/ublas/matrix.hpp> // 普通の行列用のヘッダ #include <boost/numeric/ublas/io.hpp> // ストリーム入出力用ヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 #include <iostream> using namespace std; int main() { typedef matrix<double> dmatrix; dmatrix mat1(2,2); // (i,j)も[i][j]も OK!! mat1(0,0) = 1.0; mat1(0,1) = 2.0; mat1[1][0] = 3.0; mat1[1][1] = 4.0; cout << "mat1 = " << mat1 << endl; return 0; }

出力結果は以下のようになります．読み方は[行数,列数] ((1行目....),(2行目 ....),....)のようになります．

mat1 = [2,2]((1,2),(3,4))

サイズの変更

resize()メンバ関数を使って，行列やベクトルのサイズを変更することができます．引数はベクトルに対しては変更後の次数，行列に対しては行数を列数を与える場合と次数をただ 1 つ与える場合があります．

単位行列

ublas::identity_matrix<要素の型>(次数)で単位行列への const 参照が得られます．ヘッダ が必要です．

計算

和・差・積・定数倍・定数分の1

和，差はそれぞれ+ 演算子，- 演算子を用いて記述します．定数倍は*演算子，定数分の 1 は/演算子です．

「行列どうしの積」や「ベクトルと行列の積」は prod() 関数を用います．ベクトルどうしの内積はinner_prod()を，外積にはouter_prod()を用います．

サイズの異なる行列どうしの和や差を取ろうとしたり，積が定義されない組み合わせでprod()を呼び出したりすると，例外が投げられます．

では，サンプルコードを示します．

#include <boost/numeric/ublas/matrix.hpp> // 普通の行列用のヘッダ #include <boost/numeric/ublas/vector.hpp> // 普通のベクトル用のヘッダ #include <boost/numeric/ublas/io.hpp> // ストリーム入出力用ヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 #include <iostream> using namespace std; int main() { typedef matrix<double> dmatrix; typedef boost::numeric::ublas::vector<double> dvector; // ↑今回は特に必要ないが，std::vectorを使う場合は // 混同されないように名前空間を明示する． dmatrix mat1(2,2); mat1(0,0) = 1.0; mat1(0,1) = 2.0; mat1[1][0] = 3.0; mat1[1][1] = 4.0; dmatrix mat2(2,2); mat2(0,0) = 3.0; mat2(0,1) = 4.0; mat2[1][0] = 5.0; mat2[1][1] = 7.0; dvector vec(3); vec(0) = 2.0; vec(1) = -3.0; vec(2) = -1.0; // 和 cout << "mat1 + mat2 = " << mat1 + mat2 << endl; // 差 cout << "mat1 - mat2 = " << mat1 - mat2 << endl; // 積 cout << "mat1 * mat2 = " << prod(mat1,mat2) << endl; cout << "mat2 * mat1 = " << prod(mat2,mat1) << endl; cout << "mat1 * vec = " << prod(mat1,vec) << endl; cout << "vec * mat1 = " << prod(vec ,mat1) << endl; // 定数倍 cout << "mat1 * 5 = " << mat1 * 5.0 << endl; // 定数分の1 cout << "mat1 / 5 = " << mat1 / 5.0 << endl; return 0; }

出力結果は以下のようになります．

mat1 + mat2 = [2,2]((4,6),(8,11)) mat1 - mat2 = [2,2]((-2,-2),(-2,-3)) mat1 * mat2 = [2,2]((13,18),(29,40)) mat2 * mat1 = [2,2]((15,22),(26,38)) mat1 * vec = [2](-4,-6) vec * mat1 = [2](-7,-11) mat1 * 5 = [2,2]((5,10),(15,20)) mat1 / 5 = [2,2]((0.2,0.4),(0.6,0.8))

発展

double よりも小さな型を要素とする行列やベクトルに対し，積を double 精度で計算させるprec_prod(), prec_inner_prod()の各関数が用意されています．外積だけ prec_ 版がないのは…なぜ？

prec_ 版が存在する理由は，おおよそ次のようなものであると考えられます．多次元ベクトル同士で内積をとると多数の項の加算が発生するため，特に最後の方の項を加算する段階で情報落ちが発生する危険があります．

これを回避する（おそらく）最善の策は次のようなものです．まず最初に，加算される項を正のグループ(P)と負のグループ(N)に分けます．P，N それぞれの総和を絶対値の小さい順に足していく方法で求めます．最後に P + N をもって最終的な答えとするものです．しかし，この方法は項の分類やソートにかかる計算量が多いため，最高の精度が求められる場合以外にはお勧めできません．

double 精度で計算した場合，1000 項くらいまでの加算なら，情報落ちは大抵問題にならないため，boost::numeric::ublas では double 精度で計算する関数を提供することでこれに代えているのだと思われます．

連星系は何ですか

prec_outer_prod()がないのは，外積では足し合わされる項が少ないからでしょう．
# そういえば 3 次元ベクトル同士以外の外積ってどういう定義なんでしょうね…

連立方程式を解く

以下のような連立方程式を解くことを考えましょう．

行列とベクトルを使って書くとこんな感じです．

サンプルコードを以下に示します．筆者が確認した限りでは，boost::numeric::ublas に直接線形方程式を解く関数は用意されていません．そこで，LU分解，前進消去，後退代入を行うコードを自分で記述します．

#include <boost/numeric/ublas/matrix.hpp> // 普通の行列用のヘッダ #include <boost/numeric/ublas/vector.hpp> // 普通のベクトル用のヘッダ #include <boost/numeric/ublas/triangular.hpp> // 三角行列用のヘッダ．前進消去，後退代入に必要 #include <boost/numeric/ublas/lu.hpp> // LU分解，前進消去，後退代入用のヘッダ #include <boost/numeric/ublas/io.hpp> // ストリーム入出力用ヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 #include <iostream> using namespace std; int main() { typedef boost::numeric::ublas::vector<double> dvector; typedef boost::numeric::ublas::matrix<double> dmatrix; // 連立１次方程式 // x1 + 2*x2 + 3*x3 = 24 // 2*x1 + x2 - 3*x3 = -9 // 4*x1 - 3*x2 + 2*x3 = 2 // を解くサンプル // 右辺の定数ベクトル dvector b(3); b[0] = 24.0; b[1] = -9.0; b[2] = 2.0; // 左辺の係数行列 dmatrix 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] = -3.0; A[2][0] = 4.0; A[2][1] = -3.0; A[2][2] = 2.0; // LU分解の際の枢軸選択によって起こる行交換の情報が入る行列 permutation_matrix<> pm(A.size1()); // 係数行列と定数ベクトルの表示 cout << "A = " << A << endl; cout << "b = " << b << endl; // LU 分解します．ここで A は LU 分解の結果で置き換わります． // この関数通過後の A は，下三角部分が L 行列，上三角部分が U 行列に // なっています．ただし，対角成分は U 行列に属し，L 行列の対角成分は // 1 とみなします． // pm には行交換の情報が入ります．この結果が[2,2,2]だった場合， // まず最初に0行目と2行目が入れ替わり，その状態で次に1行目と2行目が // 入れ替わり，さらにその状態から2行目と2行目が入れ替わったことを // 示します． lu_factorize(A,pm); cout << "行交換情報 : " << pm << endl; // 行交換を考慮した前進消去と後退代入 // この関数によって，b は解ベクトルで置き換わります． // もとの b を取っておきたい場合は事前にコピーしておかなければ // なりません． lu_substitute(A,pm,b); // 解の表示 cout << "Ax = b の解ベクトル x = " << b << endl; return 0; }

出力は次のようになり，解 (x₁,x₂,x₃) = (1,4,5) が得られました．

A = [3,3]((1,2,3),(2,1,-3),(4,-3,2)) b = [3](24,-9,2) 行交換情報 : [3](2,2,2) Ax = b の解ベクトル x = [3](1,4,5)

係数行列が同じである複数の連立方程式を一度に解く

次は，以下に示すような連立方程式の組を一度に解きます．

この連立方程式の組は，行列を用いて以下のように一つの式にまとめることができます．

では，サンプルコードを示します．

#include <boost/numeric/ublas/matrix.hpp> // 普通の行列用のヘッダ #include <boost/numeric/ublas/triangular.hpp> // 三角行列用のヘッダ．前進消去，後退代入に必要 #include <boost/numeric/ublas/lu.hpp> // LU分解，前進消去，後退代入用のヘッダ #include <boost/numeric/ublas/io.hpp> // ストリーム入出力用ヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 #include <iostream> using namespace std; int main() { typedef boost::numeric::ublas::vector<double> dvector; typedef boost::numeric::ublas::matrix<double> dmatrix; // 係数行列が同じである連立１次方程式の組 // x1 + 2*x2 + 3*x3 = 24 // 2*x1 + x2 - 3*x3 = -9 // 4*x1 - 3*x2 + 2*x3 = 2 // と // y1 + 2*y2 + 3*y3 = 9 // 2*y1 + y2 - 3*y3 = -6 // 4*y1 - 3*y2 + 2*y3 = 17 // を一度に解くサンプル // 右辺の定数行列(=縦ベクトルの組) dmatrix B(3,2); B[0][0] = 24.0; B[1][0] = -9.0; B[2][0] = 2.0; B[0][1] = 9.0; B[1][1] = -6.0; B[2][1] = 17.0; // 左辺の係数行列 dmatrix 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] = -3.0; A[2][0] = 4.0; A[2][1] = -3.0; A[2][2] = 2.0; // LU分解の際の枢軸選択によって起こる行交換の情報が入る行列 permutation_matrix<> pm(A.size1()); // 係数行列と定数ベクトルの表示 cout << "A = " << A << endl; cout << "B = " << B << endl; // LU 分解します． lu_factorize(A,pm); // 行交換を考慮した前進消去と後退代入 lu_substitute(A,pm,B); // 解の表示 cout << "AX = B の解行列 X = " << B << endl; return 0; }

逆行列を求める

連立方程式を解く方法を応用して，逆行列を求めることができます．逆行列を求めることは，次の等式を満たす行列　X を求めることに他なりません．

ではサンプルコードを示します．

#include <boost/numeric/ublas/matrix.hpp> // 普通の行列用のヘッダ #include <boost/numeric/ublas/triangular.hpp> // 三角行列用のヘッダ．前進消去，後退代入に必要 #include <boost/numeric/ublas/lu.hpp> // LU分解，前進消去，後退代入用のヘッダ #include <boost/numeric/ublas/io.hpp> // ストリーム入出力用ヘッダ using namespace boost::numeric::ublas; // boost::numeric::ublas 名前空間を使用 #include <iostream> using namespace std; int main() { typedef boost::numeric::ublas::vector<double> dvector; typedef boost::numeric::ublas::matrix<double> dmatrix; // 逆行列を求めるサンプル // 右辺の定数行列を3次単位行列に dmatrix B = identity_matrix<double>(3); // 左辺の係数行列 dmatrix 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] = -3.0; A[2][0] = 4.0; A[2][1] = -3.0; A[2][2] = 2.0; dmatrix CopyOfA(A); //元の行列のコピーを取っておく // LU分解の際の枢軸選択によって起こる行交換の情報が入る行列 permutation_matrix<> pm(A.size1()); // 係数行列と定数ベクトルの表示 cout << "A = " << A << endl; cout << "B = " << B << endl; // LU 分解します． lu_factorize(A,pm); // 行交換を考慮した前進消去と後退代入 lu_substitute(A,pm,B); // 解の表示 cout.precision(4); //小数点以下4桁まで表示 cout << "A の逆行列 X = " << B << endl; cout << "A * X = " << prod(CopyOfA,B) << endl; return 0; }

結果は下の通りです．A × X が確かに単位行列になっています．

A = [3,3]((1,2,3),(2,1,-3),(4,-3,2)) B = [3,3]((1,0,0),(0,1,0),(0,0,1)) A の逆行列 X = [3,3]((0.1014,0.1884,0.1304),(0.2319,0.1449,-0.1304),(0.1449,-0.1594,0.04348)) A * X = [3,3]((1,5.551e-017,0),(-1.11e-016,1,8.327e-017),(0,0,1))

Jan. 13th, 2005

prec_inner_prod, prec_prod 関数に関する記述を追加

Nov. 25th, 2004

日本語が変な箇所を修正しました．内容自体に変更はありません．

Nov. 10th, 2004

初版

[1] boost.org
[1'] Boost 翻訳プロジェクト 公開ページ
[2] 稲葉一浩 : Boost C++ Library プログラミング，秀和システム，2004年
[3] Let's boost
[4] boost user - uBLAS

These are our most popular posts:

Boost 数学関係ライブラリの使い方

2005年1月13日 ... boost::numeric::ublas 線形代数ライブラリの使い方連立方程式を解く・逆行列を求める . DT Specials ... 間違いなしです．内容の割にお値段も2800 円と手頃なため，是非 手元に置いておくことをおすすめします． ... 文献[2]によれば，expression template とは ，「式の情報ををテンプレートの構造として保持する」ような技法らしいです． ... 何が使え て何が使えないのかは私には分かりません．gcc も 2.95.3 は問題があるそうなので 3.3.3 以降を使いましょう． .... 行列とベクトルを使って書くとこんな感じです． read more

Books:FOP:BugReport

オートマトンについては、出力の一部を入力にフィードバックするというより巧妙な ところがある。 ... する。「proc、パターン、→、コマンド」と書く。 ... one often provides a pure function for the base case (trees of one element) and an expression for trees of pairs, ... などと定義を入れないと、以降の補題4.1、補題4.2の意味が定まらない。 ... 型のオーバーロードがある(ため、テストでどの型の値を与えるべきか決定できない)から である。 ..... 他の多くのデータ型」ではそれができるのかできないのか、いずれにも解釈 できる。 read more

SQL 言語

以下に示すのは、リレーショナル代数やリレーショナル論理の一部でない SQLによって 提供されるいくつかの追加機能のリスト です。 ... { UNION | INTERSECT | EXCEPT [ ALL ] } select ] [ ORDER BY expression [ ASC | DESC | USING operator ] [, . ... 例えば、 PART テーブル の 2倍の値段がどれくらいなのか知りたければ、以下のよ うな問い合わせを発行することになります。 .... GROUP BY と集約演算子を使っている 問い合わせが意味をなすためには、 ターゲットリストの中に直接記述できる属性は GROUP BY で ... read more

Q＆A (初級者コース)/11 - RjpWiki

複数のデータセットにフィットするために; 対応分析について; バイオインフォマティクスで の利用について; 時系列データのグラフ化 ... 各ベクトル要素の差の2乗を求めたい; Mac で、sourceをすると、不正なマルチバイト文字列がある; クラスターから主要な語句を ... expression()でsymbol等が表示できなくなりました; ヒートマップの目盛りの値の変更 について; data.frameへのapply適用で .... Rファイルを開いたら、自動でオンにするため の何か良い解決方法などはございますでしょうか？ ...... にその意味を込めたつもりで した． read more