精度保証付き数値計算
数値計算ツールFAQからページの名前を変更しました(2001/1/20)。<トピックス>
[2002/6/12] 石井(九大)さんからLorenzアトラクタの存在検証がW.Tuckerによってなされていることを教えていただいた。これはSmaleの第14問題(Lorenz微分方程式のアトラクタはWilliams, Guckenheimer, とYorkeが提案したgeometric Lorenz attractorであるか?)と関連する。そして、 Tucker (2002) はこれを肯定的に解いた(pdf,180k)。（中尾先生のご指摘、関連のページ)
[2002/1/27] Algae3.6.2になり( 11/1/01) patch Minor installation improvements によって、Redhat Linux7.2など、Rlabに似た数値計算ツールAlgaeがソースから問題なくインストールできるようになった。
[2001/12/29] Kearfottから次のようなMooreの初期の論文がダウンロードできるというメールがきた。
Sun Microsystems has bought the copyrights to two additional early papers by Moore, and I have made them available from the URL: http://interval.louisiana.edu/Moores_early_papers/bibliography.html The citations to these papers in that short bibliography are as follows: [6] R. E. Moore. The automatic analysis and control of error in digital computing based on the use of interval numbers. In L. B. Rall, editor, Error in Digital Computation, Vol. I, chapter 2, pages 61-130. John Wiley and Sons, Inc., New York, 1965. John Wiley and Sons, Inc. has given permission for this paper to be posted indefinitely on the non-password protected website http://interval.louisiana.edu/Moores_early_papers/Moore_in_Rall_V1.pdf for use by the interval research community. [7] R. E. Moore. Automatic local coordinate transformations to reduce the growth of error bounds in interval computation of solutions of ordinary differential equations. volume II of Error in Digital Computation, pages 103-140. John Wiley and Sons, Inc., New York, New York, 1965. John Wiley and Sons, Inc. has given permission for this paper to be posted indefinitely on the non-password protected website http://interval.louisiana.edu/Moores_early_papers/Moore_in_Rall_V2.pdf for use by the interval research community.
数値計算ツールに関する若干の補足情報
OctaveをWindowsにインストールする(大澤 洋祐君)
有理数演算のオブジェクト指向言語
著者の本に関連するページへのリンク
「Linux数値計算ツール」への質問
ドイツでの精度保証付き数値計算ワークショップと3大学訪問記
ドイツ精度保証数値計算の旅 (2000/12/10から2000/12/22)
精度保証付き数値計算の研究者のページ
山本哲朗先生
中尾充宏先生
西田 孝明先生
山本野人さん
渡部 善隆 さん
柏木雅英さん
村重 淳先生
神澤雄智さん
相馬隆郎さん
中谷祐介さん
数値計算インタプリタを作ろう(2001/8/18開設, 精度保証付き数値計算ツールSlabが完成している．公開(2002/8/15))
Slab:精度保証モードをもつMATLAB的な数値計算ツール(日本応用数理学会講演スライド(2002/9/19))
ヒルベルトの第18問題の一つが次の問題である。これが精度保証付き数値計算で解決した。
Thomas C. Hales:``Cannonballs and Honeycombs'', Notices of the AMS, Vol. 47, No.4 (2000/4) pp.440-449 には立方体に球を充填する際に，どのような方法が最も多く球を充填できるかという問題に対する，ケプラーの予想が正しいことの計算機援用証明が示されている。ケプラーは天体物理の学者であると同時に微積分などにも多くの貢献をした優れた数学者であったが，300年以上証明がつけられていなかったケプラー予想の証明が精度保証付き数値計算によって，つい最近解かれた。充填問題は情報理論でも応用のある基本問題である。
B. M. Brown, D. K. R. McCormack, and A. Zettle:``On the existence of an eigenvalue below the essential spectrum'', Proc. Royal Society London, Vol. 455 (1999) pp.2229-2234ではSturm -Liouville問題の未解決問題が精度保証付き数値計算で解決されている。
いよいよ，精度保証付き数値計算の本格的な応用が始まった。
CygwinとOctavecygwinをインストールして，その上にoctaveをインストールする。tarボールを展開して，インストールしたが，make一回ではうまくいかない。しかし，2度目にはうまくmakeが成功し，make installもうまくいく。数回makeしないとだめなのが一般的のようだ。Linux上でのインストールに比べて，非常にコンパイル時間がかかる。これが速いときもあるという情報もある。大体30分はかかる。mkoctfileはうまくいかないようだ。これは卒論の課題にしたが，うまく解決できなかった。丸めの制御が，したがって，mkoctfileではできない。ソースをいじって，丸めの制御の命令を組み込む方がはやい気がするが，まだ実行していない。octaveのcygwin上での実行速度は，Linux上の1.5倍から2倍ほど遅くなるようだ。 (2001/3/29)
[2001/4/17] Scilab2.6が出ている。
[2001/4/17] MLDで面白いページ「Linux Applications for Science & Engineering -- Linux で科学しよう！」が作られている。
精度保証付き数値計算とは何かについては著者による次の資料を参照されたい。1. 精度保証付き数値計算入門 (Scilab, Octave, Matlabを使った精度保証付き数値計算 (旧))
2. 高速精度保証付き数値計算
3. 計算機で微分方程式の厳密解を求める(``科学'(岩波書店)'1996年6月号 pp.437-445)
精度保証付き数値計算を行うためには，まず，数値計算ツールについて熟知する必要があります。
線形計算のためのツールとして、メモリ階層(キャッシュヒット率)を考えた、 信頼性のあるパッケージとしてLAPACK(レイパックと発音する)がある。 これは、Oak Ridge National Laboratory,Knoxville, TennesseeとAT&T Bell Laboratories, Murray Hill, NJに置かれている Netlib サイトからダウンロードできる。日本のミラーサイトとしては phase がある。LAPACKのユーザ マニュアル第3版がある。第2版は日本語訳がある。LAPACKはもともとFORTRANで書かれているライブラリであるが、 CLAPACKはそのC版である。LAPACKはCPUによる違いをBLAS(Basic Linear ALgebraic Subprograms)によって吸収してある。すなわち、 線形代数の基本演算をBLASの関数にし、これを標準化している。そして、LAPACKの関数をBLASを用いて書くことによってポータブル化を実現している。自分の使用するCPUに対して最適化されているBLAS(アセンブラ)を用いることによって、Cで書かれたBLASを 用いた場合の十倍程高速化が達成されている。線形方程式の直接解法に関しては、最適化BLASとLAPACKの組み合わせは、最も 優れた組み合わせの一つとなる。Pentium II上のLinux用の最適化 BLASがある。 現在はまだPentium III用のBLASは開発されていないようである(最近ベータ版が出た)。Cのコンパイラの 特性を考えてループ展開などの技法を駆使して、Cで書かれているにもかかわらず、コンパイル時にループ展開の数などのパラメータを 最適化することによって最適化BLASに匹敵するオブジェクトコードを生成する試みもある。ATLAS である。原理的には、どのようなCPUとOSであっても、Cコンパイラがあれば、 CでかかれたATLASをコンパイルするだけで最適化BLASに匹敵するコードが得られる。線形代数計算の直接解法に関しては 最も高速で、信頼性を求めるならLAPACKと最適化BLASを実装するのが有力な解となる。 LAPACK++はC++ versionであるが、LAPACKの一部の関数しか実装されていない。また、分散コンピュータ上のLAPACKとして ScaLAPACKが開発されている。数値計算パッケージの紹介としてはつぎのようなサイトもあるGAMS List of Packages。
ATLASについてはここに少し詳しい解説をおいた。また詳しい数値計算ツールの解説については次の資料を作成したので参照されたい。
また，関連して dgemmの高速化に対応するプログラミングをするという荻田氏の行列乗算高速化プロジェクト（このプロジェクトで使用されたCとFORTRANのサブプログラムも公開されたようである)も面白い。
精度保証を高速に行うためには，数値計算ツールを高速化する必要がある。以下，各種数値計算ツールについて高速化の観点を中心にして解説する。
(1) Matlab LAPACKなどの数値計算ライブラリは高速、高信頼であるが、CやFORTRANから呼び出して使うので、 プログラム開発に時間がかかる。 そこで、LAPACKなどのライブラリを基礎にして、使い易いGUIをもたせた数値計算用インタプリタパッケージが開発されている。 The MathWorks, Inc.から発売されているMatlabはその代表である。 欧米では、研究者から、学生までこのようなパッケージを利用して数値計算をするのが 一般的である。
Matlab6ではLAPACKへの対応が進み，ATLASによる高速化が果たされている。
(2) Octave Octaveは MATLABに最も近いMATLABクローンである。Linux用のソフトであるが、Windowsに実装することもできる。 まだ、オブジェクトプログラミ ングが実装されていない点ではMATLABに譲るが、 線型方程式の解や固有値を計算する組み込み関数の計算速度、信頼性ではMATLABに勝るとも劣らないと思われれる。 グラフィックスはgnuplotを利用している。また、OctaveはLAPACK対応が完全ではないがATLASを使ってOctaveを高速化できる部分も多い。
実際，Octave2.1.xxはATLASをコンパイルし、生成されたすべての*.aを/usr/local/libに入れたあとで./configure, make, make installするとA*Bやlu(A)などは高速化される。ただし、A\bは高速化されないようである。
(CPU:Intel Celeron 500MHz MEM:192MB OS:Vine Linux2.1，谷口君のレポートより転載，a\bもうまくやれば高速化できる)
ATLAS を適用しない場合 ATLAS を適用した場合 a*b 52.030 5.523 lu(a) 17.431 2.602 qr(a) 66.939 27.024 a\b 16.064 16.236 そこで、lu分解をした後、後退代入と、前進代入をoctaveのプログラミング機能で実現したm-ファイル(lusol.m)を作成した。Ax=bを考える。PentiumIII600MHzで1000行1000列の行列のlu分解は2.5秒程度である。lusolを使ってlusol(A,b)と解くと、4秒程度となる。A\bで解くと15秒程なので、これでも大幅な高速化になる。
(3) Scilab Scilabは typed listによりオブジェクト指向プログラミ ングができる環境があり、MATLABクローンとしては一番完成している。 グラフィックスも独自の優れたものが実装されている。MAPLEとのリンクもできる。LinuxでもWindowsでも実装できる。
(4) RlabはMatlab クローンではないことをうたっているが，Matlabにかなり近い数値計算インタープリタ言語。開発がここ10年ぐらいと比較的最近なため，線形数値パッケージとしては，LAPACKとBLASが使われている。
(2001/6/9) glibc2.2ではソースからのコンパイルがうまくいくことがわかった。また，丸めの方向指定の関数も組み込めた。以下，インストールの記録である。
Red Hatの最新版が非常にきれいにインストールできたので，気分良くRlabのソースを読み始める。Rlabはbisonとflexとgcを使ったインタープリタで，コンパイラの講義で使いたくなるような実にオーソドックスな作りである。
さて，Rlabへ丸めのモードの変更の命令を組み込むことを，まず，試みた。ルートディレクトリのbltin1.cに
#include <fpu_control.h>
static int _RoundDown = _FPU_RC_DOWN|_FPU_DEFAULT;
static int _RoundUp = _FPU_RC_UP|_FPU_DEFAULT;
static int _RoundNear = _FPU_RC_NEAREST|_FPU_DEFAULT;
Ent *
Down (int nargs, Datum args[])
{
Ent *rent;
/* Check n_args */
if (nargs != 0) {
rerror ("down: no arguments allowed");
}
asm volatile("fldcw _RoundDown");
rent = ent_Create ();
ent_data (rent) = mdr_CreateScalar (0.0);
ent_type (rent) = MATRIX_DENSE_REAL;
return (rent);
}
Ent *
Near (int nargs, Datum args[])
{
Ent *rent;
/* Check n_args */
if (nargs != 0) {
rerror ("near: no arguments allowed");
}
asm volatile("fldcw _RoundNear");
rent = ent_Create ();
ent_data (rent) = mdr_CreateScalar (0.0);
ent_type (rent) = MATRIX_DENSE_REAL; return (rent);
}
Ent *
Up (int nargs, Datum args[])
{
Ent *rent = 0; clock_t toctime; double elapsed;
/* Check n_args */
if (nargs != 0) {
rerror ("up: no arguments allowed");
}
asm volatile("fldcw _RoundUp");
rent = ent_Create ();
ent_data (rent) = mdr_CreateScalar (0.0);
ent_type (rent) = MATRIX_DENSE_REAL;
return (rent);
}
という3つの関数を加える。そして，bltin1.hにおいて
extern Ent *Up (int nargs, Datum args[]);
extern Ent *Down (int nargs, Datum args[]);
extern Ent *Near (int nargs, Datum args[]);
という宣言をする。そして，init.cにおいて
{BLTIN, "up", Up},
{BLTIN, "down", Down},
{BLTIN, "near", Near},
という3行を加える。この修正をした後make,make installすれば，up()で上へ，down()で下へ，near()で最近点への丸めのモードの切り換えができるようになる。
Rlabはcode.cに書かれた仮想計算機を動かすインタプリタである。
(5) 数値計算ツール Euler (2001/3/5)
EulerはDr. R. Grothmann が開発したものである。Dr. R. Grothmann のホームページに詳細があるので参照されたい。Windows版とLinux(Unix)版がある。ソースが含まれているのはLinux版である。両方試してみたが，以下では，ソースをいじるので，Linux版についてのみ議論する。ほぼ，Matlabクローンであるが，Matlabとずいぶん機能が違う。Eulerだけの特長としては，完全精度内積を実装していることである。ソースをftpで取得する。 eulersrc.zipがそうである。これを，ホームにおいて解凍する。以下，Pentium III 600MHzマシンにVine Linux2.0という組み合わせでのインストール結果を報告する。
unzip eulersrc.zip
すると，./euler/というディレクトリができるので，./euler/sourse/にcdする。まず，makefileの中のX11をX11R6に変更する(1個所)。次に，メモリサイズの指定を大きくする。そのために，./euler/sourse/x11/の中の，sysdepx11.cppの中で次の変更をする。例えば，
long memsize = 1024*1024l を long memsize = 64*1024*1024l
とした。次に，IEEE754の丸めの制御ができるようにする。そのために，丸めの制御用の組み込み関数を作った。以下，./euler/sourse/の中のファイルを変更していく。extend.hに
viod msetround(header *hd);
という行を一行加える。extend.cppのbuilintyp builtin_list[] = {{"setround", 1, msetround},...}を加えた(橙色の部分を加える)。更に，funcs.hに
void msetround(header *hd);
という一行を加える。そして，funcs.cppに
#include<fpu_control.h>
と
int _RoundUp = _FPU_RC_UP|DEFAULT;
int _RoundDown = _FPU_RC_DOWN|DEFAULT;
int _RoundNear = _FPU_RC_NEAREST|DEFAULT;
を加える。そして，関数の中で，
double rsetround (double x)
{ if (x < -1 || x > 1) {error=1;return 2;}
{ if (x == 1) {asm volatile("fldcw _RoundUp"); return 1;}
{ if (x == -1) {asm volatile("fldcw _RoundDown"); return -1;}
{ if (x == 0) {asm volatile("fldcw _RoundNear"); return 0;}
}
void csetround(double *x, double *xi, double *y, double *yi)
{
}
void msetround(header *hd)
{ spread1(rsetround,csetround,hd);
}
という3つの関数を加える。これで，変更は終了。./euler/sourse/でmakeをする。すると，./euler/progs/にxeulerが実行ファイルとして生成される。このディレクトリにおいて./xeulerでeulerは実行できる。このように変更を加えたxeulerにおいては丸めの変更が実装されている。すなわち，
> setround(1)
で上への丸め，
> setround(-1)
で下への丸め，
> setround(0)
で最近点への丸めとなる。
> A=random(1000,1000);
> b=random(1000,1);
> t=time();c=A\b;time()-t
とすると，60秒ぐらいで解が得られる。高速化したOctaveの20倍(高速化していないOctaveの4倍)ぐらい遅いが，面白い特長があるので，使い道があるであろう。何しろ，以上のように簡単に組み込み関数が作れるのは魅力的である。
精度保証付き数値計算で必要になるIEEE754(その日本語訳)の丸めの制御法はOS，プログラミング言語により異なる。
0. What machines support IEEE 754?
1. Cによる丸めの制御を行う方法を示す。(Cの丸めの制御)
2. Linux上のOctaveの丸めの制御(Octfileの利用)
3. Windows上のScilabの丸め(Dll) (関連する廣瀬君のレポート)
4. 荻田 武史氏(FortranおよびC言語における浮動小数点演算の丸め方向指定方法)
5. Macでの丸め(数値演算ユニットの制御)
6. MaTXでの丸め(古賀先生)
7. 桂田 祐史氏(FreeBSD, Solaris 2.6 (SPARC) などの環境で、 丸めモードの制御をするための命令 Up(), Down(), Near() を実現する round.h )
9. Matlab6における丸めの制御
2回程Matlab6についての丸めの制御法について質問を受けたので，以下にその方法を示す。丸めの制御の命令は
>> system_dependent('setround',mode);
である。ただし，
mode=-Infのときround downへ
mode=+Infのときround upへ
mode=0のときround to zeroへ(チョッピング)
mode='nearest'のときround to nearestへ(default mode)
である。 (2001/3/8)
10. 関連のリンク
Intel上のCでのIEEE754のdoubleの64bitを書き出すプログラム(print the bits of double real number)
IEEE Standard 754 Floating Point Numbers
非線形問題の精度保証付き数値計算で必要になる区間演算についてはオブジェクト指向言語による区間クラスの実装が便利である。
1. C++の区間クラス
2. Sun FORTRAN95コンパイラには区間演算が定義されている。
線形問題の精度保証付き数値計算は著者による丸めの制御方式精度保証法が高速で有効である。
著者によるMATLABの精度保証ライブラリ(m-files)を公開していきます。
九工大の古賀さんによりMATXで丸めのモードの変更をする関数が装備され，著者の方法で行列方程式を精度保証保証する関数が作られた。
Javaで精度保証付き数値計算 Javaによる数値計算関連のリンク
Javaを科学技術計算用に用いようとするときに参考になる サイトとして
Java Numerics (ここには色々なパッケージの情報がある)
などがある。
Javaはプラットフォームインディペンデントという特徴がある。したがって，精度保証プログラムを証明の一部と見なすとき，現状のCPUやOSに依存しないで記述できているという特長がある。Javaによる精度保証はその意味で今後注目を集めるような気がする。
数値計算関係のリンク集へのリンク
Mathtools.net is a technical computing portal for all scientific and engineering needs. The portal is free and contains over 20,000 useful links to technical computing programmers, covering MATLAB, Java, Excel, C/C+, Fortran and others.
Help for Octave Matrix Numeric Tools and also MatLab Matrix Laboratory & Scilab Links
Linux Road Linuxと科学ソフトウェア総合情報サイト 数値計算システム Octave