kashiの日記

2015/10/15(木)MATLAB 2015b での精度保証

INTLAB version 9の重要性の記事で書いたように、マルチコアを持つPC (いまどきはほとんど全部そうですが) のmatlab上でintlabを使おうとすると、2012b以降は丸めの向きの変更が正しく行われずに精度保証が出来ない状態が続いていました。原因はこのバージョンから内蔵のIntel Math Kernel Libraryが変更されたことで、マルチスレッド計算を行うときに元スレッドの丸めモードを他スレッドに伝えなくなってしまったことによるものです。これを避けるため、

BLASを他のものに差し替える。速度は遅くなってしまうことが多い。
matlabを2012aで留めておく。
マルチコア使用を禁止し、シングルコアで動作させる。当然遅くなる。

のような対処が必要でした。

で、最近出たmatlab 2015bですが、Intel MKLが再び丸め変更が効くようになったという噂があり、3年ぶりにちゃんと精度保証の出来るmatlabになったかと期待されていました。しかし、Rump先生のテストによると結果は残念なものだったそうです。確かにIntel MKLはちゃんと動くようになったそうですが、今度は行列a, bに対する a .* b のような演算で丸め変更が効かなくなってしまったそうです。このような簡単な計算で外部のBLASを呼ぶとは考えにくく、恐らくこれはmatlab内部で完結した処理と考えられます。matlab自身でマルチスレッドでデータを分割して計算するような実装になっていて、そのとき他スレッドに丸めの向きの情報を伝えていないと予想されます。

たまたまうまく行っていた「他人のふんどしで精度保証する」方法論はそろそろ限界なんじゃないでしょうか。自分のコントロールできる範囲内の事柄で精度が保証できている証拠を構成するようにアルゴリズム全体を設計するべきだと思います。

コメント（0件）

2015/10/15(木)kv-0.4.25

精度保証

kvライブラリを0.4.25にアップデートしました。

今回の主な変更は、

内部型の異なる区間の間の相互変換を行うinterval-conv.hpp追加
eig.hppで重複固有値があるとうまく計算できなかった問題を修正
ガウスの超幾何関数を追加

などです。いずれの変更も、ライブラリが実際の問題に使われる中で必要になったりしたのがきっかけとなっており、とてもいい傾向だと思います。

最近はほんの遊びで追加したmpfrラッパーが使われる頻度が高くなってきており、もう少しきちんと設計しなおした方がいいかもなあ、なんて考えています。それも含めて、混合精度演算をスマートに行える全体の構成をきちんと考える必要を感じています。

コメント（0件）

2015/10/11(日)数学関数の精度を調べる

精度保証

ふと思い立って、C言語の(というかlibmの)数学関数の精度を調べてみました。

現代の一般的なプログラミング環境では、浮動小数点演算はIEEE 754 Std.に基づいて行われるのが一般的です。そしてIEEE 754 Std.においては、演算に完全精度 (フォーマットで表現可能な範囲内で最善の精度、0.5ulp) を要求しているのは

加算
減算
乗算
除算
平方根

の5種類の演算のみで、それ以外の演算の精度については何も規定されていません。つまり、どんな値が計算されても文句は言えません。

従って、C言語などから使うexp, sinなどの数学関数の返す値は信用することは出来ません。kvライブラリの区間演算でも数学関数は全て自力で実装しています。区間演算ライブラリの一つである PROFIL/BIAS では数学関数をC言語の組み込み関数で計算した後succ/predで区間を広げて計算結果としていますが、これは誤りであり、精度保証付き数値計算を行うためのツールとしては役に立ちません。

Intelの64bit環境では、数学関数は全てソフトウェアで計算されるのが一般的です。従って、計算精度はlibmの実装に依存することになります。たとえ完全精度を出す義務が無かったとしても、長い間のアルゴリズムの改良により完全精度が出るようになっているかも知れません。そこで、

CPU: Intel Core M 5Y70
OS: Ubuntu 14.04 LTS 64bit

という環境で、数学関数が実際にどのくらいの精度が出ているのか調べてみました。比較対象は、kvライブラリのinterval<dd>の計算結果をinterval<double>に直したもので、その区間の外側の値が計算されたら、誤差が1ulpより大きかったことになります。下に示したような入力範囲で100万個の乱数を発生させて調べました。結果を示します。

関数名	下限	上限	結果
exp	-700	700	OK
log	1e-308	1e308	OK
sin	-1e10	1e10	OK
cos	-1e10	1e10	OK
tan	-1e10	1e10	OK
asin	-1	1	OK
acos	-1	1	OK
atan	-10	10	OK
sinh	-700	700	NG
cosh	-700	700	OK
tanh	-700	700	NG
asinh	-1000	1000	OK
acosh	1	1000	OK
atanh	-1	1	NG

このように、sinh, tanh, atanhで精度不足が見つかりました。言うまでも無いことですが、OKとなった他の関数について、このテストで精度不足が見つからなかったからと言って精度に問題が無いことの証明にはなりません。また、このテストでは1ulp以上の誤差を検出しますが、一般に完全精度と言えば0.5ulp以下の精度のことを指します。

sinh, tanh, atanhについて、精度不足となった例を示しておきます。

sinh(0.74251904152333736419677734375)=0.8126541801799147535234624228905886...
と計算されたが、真の値はdoubleの区間
[0.8126541801799148645457648854062426,0.8126541801799149755680673479218968]
の間にある。

tanh(0.48066521994769573211669921875)=0.4467762045341044374602290645270841...
と計算されたが、真の値はdoubleの区間
[0.4467762045341044929713802957849111,0.4467762045341045484825315270427382]
の間にある。

atanh(0.2377787786535918712615966796875)=0.2424184572356452571639806592429522...
と計算されたが、真の値はdoubleの区間
[0.2424184572356452016528294279851252,0.2424184572356452294084050436140388]
の間にある。

もう記録は残っていませんが、10年以上前に調べたときはもっと精度に問題があった記憶があるので、libmの品質は上がっているのかもしれません。案外精度良かったんだな、という印象です。

コメント（0件）

2015/10/05(月)kv-0.4.24

精度保証

kvライブラリを0.4.24にアップデートしました。

今回の主な変更点は、

ライセンスを明記し、MITライセンスとした。
interval.hppのproper_subsetを修正し、proper_subset([0,1],[0,2])が真になるようにした。
beta.hppを追加。
eig.hppを追加(高安亮紀氏の提供による)。

などです。

ライセンスは今まで曖昧にしていましたが、某所から問い合わせがあって、MITライセンスと明記することにしました。大雑把に言えば、著作権表示を残しさえすれば自由にしていいよ、責任は取らないよ、というライセンスと理解しています。「煮るなり焼くなり好きにして」という気分だったのですが、やはり明記されていないと使いにくいということはありそうです。

proper_subset(X,Y)は、XがYの真部分集合のとき真、という関数ですが、X=[0,1], Y=[0,2]のようなとき偽を返していたという問題がありました(X=[0.5,1]なら真、すなわち両側に余裕がないと真になっていなかった)。この関数を使っていたのはKrawczyk法の部分だけですが、今まではこのバグによりわずかに性能が低下していた可能性があります。

eig.hppは、非対称実行列の固有値と固有ベクトルを(近似的に)計算するeig関数と、非対称実行列の固有値を精度保証付きで計算するveig関数を含んでいます。高安亮紀氏に提供していただきました。ありがとうございます。

(2015年10月6日追記)

kv-0.4.15、シンプルな区間演算ライブラリで公開した簡単な区間演算ライブラリを、0.4.24相当にアップデートしました。

interval-simple-0.4.24.tar.gz

特徴を再掲しておきます。

kvライブラリの区間演算の部分を切り出して、簡単な区間演算ライブラリを作ってみました。
テンプレートを排除しました。従って内部型はdoubleのみですが、プログラムが読みやすくなっています。
テンプレートを排除したおかげでboostに依存しません。単体で使えます。
必要なファイルはinterval.hpp一つです。
内部型のカスタマイズが出来ないこと以外は、kvライブラリに含まれるinterval.hppと機能は同等です。数学関数も多数あります。
丸めの変更は単純にfesetroundを使っています。
文字列との方向丸め付き相互変換機能もあります。従って、文字列からその文字列の表す数値を含むように区間を初期化できるし、表示も真の値を含むように外側に丸めるようになっています。
volatileを使った最適化対策はきちんとやってます。

ちょっと区間演算を使うだけならkvライブラリは少し大掛かりなので、気軽に使うならこのくらいでいいかも、と思ってます。また、ライブラリを読んで仕組みを知りたい人には読みやすいと思います。

コメント（0件）