boost::numeric::ublas 線形代数ライブラリの使い方
連立方程式を解く・逆行列を求める

DT Specials -> Boost -> boost::numeric::ublas 線形代数ライブラリの使い方

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 円と手頃なため,是非手元に置いておくことをおすすめします.

目次

boost::numeric::ublas を使うメリット

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

第一に,計算が高速らしいことです.これは 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 など)を指定し,インスタンス名の後に行数と列数を指定して宣言します.

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

要素の参照と代入

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

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

ストリームへの出力

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

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

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

サイズの変更

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

単位行列

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

計算

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

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

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

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

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

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

発展

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 次元ベクトル同士以外の外積ってどういう定義なんでしょうね…

連立方程式を解く

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

x+2y+3z=24, 2x+y-3z=-9, 4x-3y+2z=2

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

行列((1,2,3),(2,1,-3),(4,-3,2)) 縦ベクトル(x1,x2,x3) = 縦ベクトル(24,-9,2)

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

出力は次のようになり,解 (x1,x2,x3) = (1,4,5) が得られました.

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

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

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

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

逆行列を求める

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

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

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

更新履歴

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


Copyright 2004 TAKAHASHI Daisuke (130km/h Software)
Mail :