2015/04/16(木)kv-0.4.20
今回は、主に端点に無限大を持つような区間に関するbug-fixです。区間演算ライブラリは、端点に無限大を持たせることによって半無限区間や全区間を取り扱えるように設計されています。区間演算の実装は、上端下端型と、中心半径型の2種類が考えられますが、本ライブラリでは半無限区間や全区間を扱えるという上端下端型のメリットを重視し、上端下端型を採用しています。
端点に無限大が入ると、無限大-無限大、無限大/無限大、無限大*0などでNaNが発生する危険性があり、これらを注意深く排除する必要があります。今回、乗算において思慮が浅かった部分が発見され、きちんと考えなおした上で修正しました。詳細なアルゴリズムは 端点に無限大を持つ区間演算まとめを見て下さい。無限大を扱わない場合は問題ありませんが、無限大を扱った計算において、0.4.19以前のライブラリの計算結果は精度保証されていなかった可能性があります。0.4.20へのアップデートをお勧めします。
また、sinh, coshに-無限大を入れた場合の不具合も修正しました。
2015/03/23(月)kv-0.4.19
主な変更点は、
- mpfrである程度以上仮数部の長い数を扱ったとき、e, ln2, piの定数の 精度が不十分なため区間幅が小さくならない問題を修正。
- atanhの端点処理をちゃんとやるようにした。
- ddとmpfrにceilとldexpを追加。
区間演算では、テンプレートを使って、interval<double>ならdoubleにふさわしい精度で、interval<dd>ならdouble-doubleにふさわしい精度で数学関数の計算を行う仕組みになっています。この仕組みで困るのが定数の計算で、定数は入力を持たないのでどの精度で計算すべきか判断できません。そこで、(numeric_limits<T>の使い勝手に倣って) constants<T>::pi(), constants<T>::e(), constants<T>::ln2()を用意し、Tの型によって計算精度を変える仕組みを採用しています。で、これらの実装は、例えばpiなら、
template <class T> struct constants< interval<T> > {
static interval<T> pi() {
static const interval<T> tmp(
"3.1415926535897932384626433832795028841971693993751",
"3.1415926535897932384626433832795028841971693993752"
);
return tmp;
}
}
のように文字列で上端下端を与えそれを上下に丸めることによって得ていました。この定数は10進数50桁で円周率の上界下界になっています。doubleとddしか無ければこれだけ精度があれば十分だったのですが、mpfrで50桁以上の精度を指定した場合、(精度保証は出来ているものの)定数の区間幅が広すぎてまともな精度が出ない事態になっていました。ここの仕組みを変えて、きちんとその型なりの定数の計算方法を書くようにしました。2015/03/12(木)kv-0.4.18
0.4.17で、lgammaにdouble以外の型を入れるとコンパイルできないバグを混ぜてしまっていたので修正しました。いつもならもう少し修正が貯まってから公開するのですが、某所でlgammaが必要だったのですぐに公開しました。ついでに、第1種Bessel関数(整数次のみ)の簡単な実装を追加しました。
2015/03/07(土)kv-0.4.17
主な変更点は、
- 特殊関数の説明文を書いた。特にガンマ関数の計算方法の詳細を書いた。
- 内部で例外として使っていたrange_errorをdomain_errorに変更。
- 数学関数を含む常微分方程式が正しく解けないことがあったバグを修正。
2番めの例外についてちょっと補足しておきます。区間演算やaffine arithmeticなどでゼロ除算や負の数の平方根などが出てきたとき、プログラムをただちに停止せず例外をthrowするようになっています。そして、常微分方程式や全解探索などの上位プログラムでは解くべき問題の関数評価をするときにtry-catchでこの例外をcatchし、例外が起きてもすぐに諦めずにステップ幅を小さくしたり区間を分割したりするような実装になっています。この例外として、0.4.16以前はrange_errorという例外を使っていたのですが、どうやらdomain_errorの方がふさわしそうなので、プログラム全体でdomain_errorを使うように変更しました。
普通に使ってるだけならあまり関係のない変更ですが、例えばaffine.hppをコピーして改造して使っているような方がいれば、同じように変更した方がいいでしょう。単に一括置換するだけで大丈夫ななずです。
忙しい1月2月を過ぎて、やりたいこともいろいろ溜まっているのでここから頑張ろう。
2015/01/27(火)INTLAB version 9の重要性
(なぜかページが2つあって微妙に中身が違うので両方載せておこう。)
まったく予想してなかったのですが、このINTLABがversion 9にアップデートして、新たにOctaveに対応したそうです。Octaveは、フリーの有名なmatlabクローンです。と言っても、完璧な互換性を目指してるわけでは無さそう。
INTLABはRump先生の大変な努力で作成された素晴らしいソフトウェアで、この分野のデファクトスタンダードと言えるものです。ところが、最近はきちんと動かないことが多く、周囲の研究者に聞くと、matlabのversionを2012aで留めておかないと、いろいろまずいことがおきるのだそうです。
これは、matlabのせいと言うよりは、matlabの中で使っているBLASの問題です。
BLAS(=Basic Linear Algebra Subprograms)は、LAPACKから呼び出される行列積などの基本的なプログラムの集合体です。近年の計算機では計算速度に比較して相対的にメインメモリが遅くなっており、キャッシュメモリを使ってその遅さをカバーしています。そのため、行列積などのプログラムはキャッシュのヒット率を極限まで高めるきわめて技巧的なプログラミングが要求され、高速なプログラムは誰でも書けるものではなくなってしまいました。実際、素人が書いた単なる三重ループとBLASの持つ行列積とでは速度が10倍以上違うのが普通です。また、CPUの性能を極限まで引き出すにはそのCPUに合わせたプログラムが必要になります。そのため、行列積などの基本的なプログラムはインターフェースだけ揃えておいて中身は自由に差し替えられるように分離するようになり、それがBLASです。
最近のmatlabで使っているBLASはIntel Math Kernel Libraryに含まれているもので、これの挙動が精度保証に使っている手法と合わなくなってしまったため、きちんと精度保証が出来なくなってしまったというわけです。
より詳しく説明します。浮動小数点数(double)を係数に持つ行列A,Bがあったとき、行列積ABを精度保証付きで計算する問題を考えます。三重ループで行列積を計算するプログラムを元に、計算過程を全て区間演算に置き換えれば精度保証付きの結果を得ることが出来ます。ところが、このようなことをすると高速なBLASの恩恵にあずかることが出来ません。そこで、
- 丸めの向きを下向きに変える
- C1=ABをBLASで計算する。
- 丸めの向きを上向きに変える。
- C2=ABをBLASで計算する。
- (丸めの向きを元に戻す)
BLASの中でマルチスレッドあるいはマルチコア計算を行っており、他のスレッドまたはコアに丸めの向きの情報が伝わらないことがあって、この場合は精度保証されていない結果が得られてしまうことになります。また、
BLASの中でStrassenアルゴリズムなどの「減算を含む」計算が使われている場合は、「計算式の中の全ての計算が下向き丸めならば最終的な計算結果が下向き丸めになる」という仮定を満たさなくなってしまい、このケースも精度保証が破れてしまいます。Strassenアルゴリズムはある程度大きい行列でないと効果がないことが知られていますが、計算機の発展に伴い大きな問題が解けるようになってきてそろそろStrassenアルゴリズムが有効に機能する段階に入りつつあると考えられています。
上記のような精度保証の破綻は、自分でBLASを書いたならば当然そういうことがありえるのかどうか判断できますが、closed sourceな他人の書いたBLASの場合はそうはいきません。Intel Math Kernel Libraryはまさにclosed sourceであり、2012年に(何があったか窺い知ることは出来ないが)何らかの変更が行われ精度保証を破綻させてしまったというわけです。
さらに考えてみれば、2012aまでのmatlabでも「観察の結果問題が無いように見える」だけであり、実は精度保証が破綻しているケースがあるかも知れません。使っているアルゴリズムも分からずソースコードも見られないプロダクトを使って精度保証が出来ると考えるほうがどうかしてるとさえ思えます。
さて、ここまでお読みいただければ、INTLABがOctaveに移植されたことの重要性がよく分かると思います。Octaveはopen sourceであり、(インストールの仕方にも拠るでしょうが)使われるBLASもatlasやopenblasなどのopen sourceなものです。中身を見れば精度保証に問題があるかどうか判断できるし、問題があるならば修正することも出来ます。
INTLABはversion 7から有償になってしまったので簡単に試せないのが残念ですが、Octaveでどのくらいちゃんと動くものか、近いうちに試してみたいと思います。
※上に書いたような問題点は、ずっと前から感じていたことですが、代替手段が無かったこともあり、あまり明確に否定してしまうと線形計算の精度保証の分野の研究者が困ると思い、あまり口にしないで来ました。INTLABがOctaveで動くようになったことで、書くなら今しかないと思い、この文章を書きました。