最終更新: 2013/7/4

boostで精度保証しよう

柏木 雅英
学生の卒論準備のために研究室内向けに作成したページですが、 一般公開することにしました。

1. はじめに

精度保証付き数値計算を行う方法はいろいろあるが、 特に非線形計算でC++を使った本格的なプログラムを書くことを 考えたとき、boostに含まれているintervalクラスとublasクラスを 使うのは一つの有力な方法だと思われる。

ここに、そのboost.intervalとboost.ublasを使って作った 比較的単純な非線形方程式の全解探索のプログラムがある。 これを動かしてみて、boostに慣れてみて欲しい。

2. boostとは

boostとは、C++の巨大なテンプレートクラスライブラリ集である。 C++はオブジェクト指向言語に分類されるが、それ以上にテンプレートによる 記述力が他の言語の追随を許さぬレベルにある点が特色である。

テンプレートとは、簡単に言えば、変数の型を未定にしたまま クラスや関数を定義する機能である。C言語のマクロ機能を発達させたものという 見方も出来る。

C++にはSTLと呼ばれる標準テンプレートクラスライブラリがあるが、 boostはそれに飽き足らぬ人達が標準化を待たずに勝手に開発した テンプレートの集合体である。

boostの和書として、

が出版されている。

3. boostのインストール

boostは巨大なクラスライブラリの集合体だが、 の二種類に分けられる。 boostの全機能を使おうと思えば、リンクするライブラリをmakeする 必要がある。その方法は適当に検索するとか本を見るとかして下さい。

intervalとublasは両方ともヘッダファイルだけで使えるクラスであり、 そのようなクラスしか使わないならば、boostパッケージのヘッダファイルだけを 抜き出してどこかに置いておくだけで問題ない。 boostパッケージを展開して出来る「boost」というディレクトリをどこかに 置いておき、コンパイル時に「-I(その場所)」とすればよい。

例えば/usr/local/includeに置いたなら、

cc -I/usr/local/include hoge.cc

カレントディレクトリに置いたなら

cc -I. hoge.cc

という感じ。

3.1 intervalへのpatch

ublasの方は頻繁に更新されているが、実はintervalの方は ほとんど更新されておらず、やや荒削りな部分が見受けられる。 プログラムを作成しながら、どうしてもここは直して欲しいという 部分がいくつか見付かったので、本体に直接手を入れることにした。

以下にその変更点をpatchの形で示す。

boostディレクトリの見えている場所で、 patch -p0 < boost-1.xx.x.diffとすればパッチを当てられる。
このpatchは、以下の内容を含む。 柏木の開発しているプログラムは、もはやpatchの当たったintervalを 前提に書いている。
手抜きするならここに を置いておくのでこれ使ってもいいよ。

4. 区間演算(interval)

4.1 概要

数値を、[下端, 上端] のような2つの数の組で表される閉区間で表し、 「結果として有り得る数値の範囲を包含するように」演算を行う方法である。 区間の上端と下端を表す数の型はテンプレートになっており、 doubleだけで無く他の型を使うことも出来る。

4.2 使い方

使い方のサンプル (test-interval.cc) を見て欲しい。基本的な使い方はこれで十分と思う。 詳細は Boost Library Documentation の中の Intervalの項を参照のこと。 日本語の解説も参考になるかも。

boost::numeric::とかstd::とか長くて面倒と 思うだろうが、using namespaceを使うのが常套手段だろうか (test-interval2.cc)。 stdならともかくintervalなんて得体の知れないものを 名前空間に引っ張り出すのは気持ち悪いと思うなら、 typedefを使う手もある (test-interval3.cc)。 namespaceの別名を使って短縮する方法もある。 (test-interval4.cc)。

4.3 crlibmと組み合わせる

区間演算を行う上でsinなどの数学関数をどうするかは重大な問題である。 boost::intervalではsin_down, sin_upなどのsinを下向き丸め、上向き丸めで 計算するような関数を「自分で用意すれば」sinの区間演算を行うことが出来る。 そこで、 CRlibm(Correctly Rounded mathematical library) と組み合わせて、boostでの数学関数の区間演算を行うことにする。

crlibmのmakeは苦労するかもだが頑張って下さい。 現状ではたぶんUNIX系OSのみで、windows不可。

boostでcrlibmを使うために、 use-crlibm.hpp を作成した。 test-crlibm.cc のようにuse-crlibm.hppをincludeし、interval型はpolicyを変更して crlibmを使うようにする。 コンパイルするときは 

g++ -I/usr/local/include -L/usr/local/lib test-crlibm.cc -lcrlibm
のように、crlibmのheaderとlibraryのありかをきちんと教えて-lcrlibmを付ける。

4.4 入出力に注意

boost.intervalは、残念なことに10進数と2進数の相互変換に関して ほとんどサポートが無い。

数値の入力に関して言えば、 test-interval5.ccのように 例えば10進数の0.1を区間演算で扱うには注意が必要である。

数値の出力については、boost/numeric/interval/io.hppを includeすれば区間を表示できるが、これは単にC++標準の方法で(最近点への丸めで) 下端、上端を表示しているだけであり、そのつもりで見る必要がある。 正しくは、下端は下向き丸めで、上端は上向き丸めで10進数に変換すべきで あるが、そうはなっていない。

4.5 丸めの向きの変更

区間演算を行う際の丸めの向きの変更はboost.intervalの内部で きちんと行われるが、自作プログラムの一部で丸めの向きを変更したいことも あるであろう。その場合は、 test-rounding.ccのように 行えばよい。 ただし、丸めの向きの変更はコンパイラの最適化が大敵で、 普通にコンパイルしたら正しく動いているように見えたものが -O3 だとうまく動かないようなことがよく起こる。 force_rounding(与えられた引数をvolatile変数に代入することに よって、最適化による計算式の除去を防ぐ)を使うなどして、 慎重に記述する必要がある。

5. 行列計算(ublas)

5.1 概要

ublasは、行列計算を行うテンプレートライブラリである。 行列やベクトルの各要素はテンプレートになっており、doubleや floatばかりでなくintervalを入れたり、自作の型を入れることも出来る。 自動微分などの特殊な型を使うことが多い非線形計算の精度保証には うってつけのライブラリであると言える。

なお、blasという名前から高速に動作するようなイメージを持つかも 知れないが、blasと同等のインターフェースを持つという意味で ublasと名付けられており、特に高速なわけではない (自分でforループを回して書いたプログラムと同程度)。

5.2 使い方

詳細は Boost Library Documentation の中の ublasの項を参照のこと。 日本語の解説も参考になるかも。 boost::numeric::ublas 線形代数ライブラリの使い方 も大変参考になる。

なお、STLの一部としてstd::vectorが存在するため、 using namespace std;している場合などに boost::numeric::ublas::vectorと衝突してエラーになりやすい (test-ublas2.cc)。 test-ublas3.ccのように ublasのvectorは明示的に指定することにしてusing namespaceに頼らないのが 安全かも知れない。

5.3 マクロNDEBUGについて

ublasではExpression Templateと呼ばれる技法を用い、 余分な一時変数を作成しないことで高速化を図っている。 但し、この機能はコンパイル時に -DNDEBUGを付けて NDEBUGマクロを定義した場合のみ有効になる。 NDEBUGはC++においてrelease版(debug版で無い)を makeする場合に付加するマクロとして一般的であり、実際Visual C++で release版としてbuildすると自動的にこのマクロが付加される。

なお、ublasでは、NDEBUGマクロが定義されていない場合は 配列の添字溢れのチェックが行われるが、NDEBUGを定義すると チェックが省かれて高速化されるようになっている。

6. 全解探索プログラムを動かす

大きめのプログラムとして、非線形方程式の全解探索を行う プログラムのサンプルを示す。高速化のための凝った工夫を取り除いた シンプルなものである。 必要なのはこれらのファイル。 test-allsol-simple.cc をコンパイル、実行すれば、 その中で定義された非線形関数の全解探索が行われる。

test-allsol-simple.cc での解きたい関数の定義の仕方が分かりにくいかも知れない。 解きたい関数ごとにクラスを定義し、そのクラスでは operator()をオーバーロードする。 そのクラスのインスタンスfに対してf(x)と書くと operator()が引数xで呼び出されるので、 そこに関数の定義を書き、インスタンスfを関数とみなす。 これは関数オブジェクトと呼ばれる技法で、関数ポインタを使うよりも 一般的にオーバーヘッドが少ない。

7. 数値型のプログラミング

ここまでは、既存のintervalとublasの使い方を説明してきたが、 精度保証付き数値計算では、自動微分などそれ以外の「数値っぽい振る舞いをする型」 を自作する必要がある。 特に自作型を入れ子にする場合など、注意すべき点も多いので、 その方法を 数値型のプログラミング にまとめた。

8. Affine Arithmetic

Affine Arithmetic にまとめた。

9. おわりに

まだ書いてない。自動微分とかpsaとかについて書く予定。 アプリケーションとして、全解探索とかODEとか数値積分とかも。

ublasで内部型がdoubleの場合にlapack/blasを使う方法、 intervalの内部型としてmpfr型を使い任意精度区間演算を行う方法、 なども書いていきたい。

変更履歴

変更履歴参照。