線形方程式はもう一度数値解を計算する手間で精度保証できる

Accuracy of Numerical Solution for Linear Equation can be Verified with n^{3}/3 flops

大石 進一 Shin’ichi Oishi

2000 年 8 月 30 日

アブストラクト本稿では精度保証付き数値計算を行うための環境を整備するための web 情報を詳細に記述する。また、線形連立方程式の数値解の精度保証が、ガウスの消去法による数値解を求めるのと同じ手間で可能となることを示す。この結果、線形系については精度保証付き数値計算がほぼ実用段階に到達したことを明らかにする。

abstract In this article, it is shown how to build computational environment to perform numerical simulation with guaranteed accuracy. For the purpose, detailed web information is cited. It is also shown that verification of numerical solution of the simultaneous linear equation can be performed with the computational n^{3}/3 flops which is same as that for calculating a numerical solution of the problem by the Gauss elimination.

1 はじめに

結果の精度が保証された数値計算を精度保証付き数値計算という。区間演算や丸め制御による数値計算によって、精度保証付き数値計算が可能であることが示されている。本学会誌でも講座として精度保証付きシミュレーション [1]∼[4] が連載され、その全体像が把握できるので参照されたい。また、精度保証付き数値計算の専門書 [5],[6],[7] が出版され、深くこの分野が学べるようにもなりつつある。数値計算とは精度保証付き数値計算であるという立場からの数値計算の入門書も出版されている [8]。一方、現在普通に行われている数値計算は結果の精度を保証することは考えられていない。これを近似計算と呼ぶことにしよう。近似計算の結果を基にして、近似解の近くに真の解が存在し、その真の解と近似解の誤差を精度保証付き数値計算で計算することを、数値計算結果の検証 (verification) と呼ぶ。また、近似解を得た後で、その誤差を計算するので、数値計算結果の検証は事後誤差解析 (a posteriori error estimate) の一種となる。一方、あまり規模の大きい計算には向かないが、近似解を求めることなく、すべての計算を精度保証付きで進めることによって、問題を解く方法も考えられている。

精度保証付き数値計算という言葉は、事後誤差解析的に真の解の存在を示し、更に近似解の誤差を計算する場合も、すべての計算を精度保証付きで進めることによって、問題を解く場合も含めて指すことが多い。本稿では、文献 [1]∼[8] を補足して、精度保証付き数値計算を行うための計算機環境を整えるための web 情報や、連立一次方程式の精度保証が、数値解をもう一回求める手間で行えるなどの精度保証付き数値計算の研究の最近の話題につき解説することを目的としている。特に、つぎの点について重点的に解説する。

1. 精度保証付き数値計算のアルゴリズムが革新的に進展し、実用段階に到達した部分が多く現れている。

2. 精度保証付き数値計算を行うための環境を整備するためのツールを web から調達し、構築する方法を示す。

3. 精度保証付き数値計算の諸問題に対する参考文献を web 上から集めるためのサイト情報を提供する。

以上を要するに、現代の研究においては、 web の利用が不可欠になっている。この点を重視した精度保証付き数値計算の解説になることを目指して本稿は執筆されている。

2 数値計算ツール

精度保証付き数値計算はやはり数値計算である。したがって、まず、数値計算のツールとしてどのようなものがあるかを知ることが重要となる。まず、数値計算に関連するサイトを紹介しよう。よくある質問とその解答のサイトとして FAQ Numerical Analysis and Associated Fields Resource Guide がある。この中には数値計算パッケージや数値計算ソフトウエアへのリンクもあり、基本的な数値計算ツールを揃えるのに便利である。また、最適化に関連する数値計算のサイトとしては Bibliography for Optimization with Sensitivity Analysis がある。

2.1 数値線形代数の直接解法ライブラリ

まず、線形計算のためのツールとして、メモリ階層(キャッシュヒット率)を考えた、信頼性のあるパッケージとして 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 Library) によって吸収してある。すなわち、線形代数の基本演算を BLAS の関数にし、これを標準化している。そして、 LAPACK 本体を BLAS を用いて書くことによってポータブル化を図っている。自分の使用する CPU に対して最適化されている BLAS(アセンブラ) を用いることによって、 C で書かれた BLAS を用いた場合の数十倍の高速化が達成されている。線形方程式の直接解法に関しては、最適化 BLAS と LAPACK の組み合わせは、最も優れた組み合わせの一つとなる。 Pentium II 上の Linux 用の最適化 BLAS がある。現在はまだ Pentium III 用の BLAS は開発されていないようである。 C のコンパイラの特性を考えてループ展開などの技法を駆使して、 C で書かれているにもかかわらず、コンパイル時にループ展開の数などのパラメータを最適化することによって最適化 BLAS に匹敵するオブジェクトコードを生成する試みもある。 ATLAS や PHIPACK などである。尚、 PHIPACK は dgemm という行列の乗算プログラムだけが対象にされている。原理的には、どのような CPU と OS であっても、 C コンパイラがあれば、 C でかかれた ATLAS や PHIPACK をコンパイルするだけで最適化 BLAS に匹敵するコードが得られる。線形代数計算の直接解法に関しては最も高速で、信頼性を求めるなら LAPACK と最適化 BLAS を実装するのが有力な解となる。 LAPACK++は C++ version であるが、 LAPACK の一部の関数しか実装されていない。また、分散マシン上の LAPACK として ScaLAPACK が開発されている。数値計算パッケージの紹介としてはつぎのようなサイトもある GAMS List of Packages。

注意

1. 数値計算のコードを載せた本として、しばしば利用される向きのある Numerical Recipes は数値計算の専門家でない著者による書であるが、つぎのような批判 Why not use Numerical Recipes 他がある。これらの批判はこの本に載せてあるコードが最適なものに比べて実行速度が数十倍遅く、また、信頼性も低い場合があるという点に集中している。 Netlib などの他の優良なソフトを利用すべきであろう。

2. 本稿では反復解法は扱わないが、数値線形代数の反復解法のサーベイがある。

2.2 数値計算パッケージ

LAPACK などの数値計算ライブラリは高速、高信頼であるが、 C や FORTRAN から呼び出して使うの、プログラム開発に時間がかかる。そこで、 LAPACK などのライブラリを基礎にして、使い易い GUI をもたせた数値計算用インタプリタパッケージが開発されている。 The MathWorks, Inc. から発売されている Matlab はその代表である。欧米では、研究者から、学生までこのようなパッケージを利用して数値計算をするのが一般的である。

MATCOM(Matlab to C++ Compiler は Matlab から C++へのコンパイラ、すなわち Matlab コードから MEX ファイルあるいは C++コードを生成するフリーなツールである。生成された C++コードを最適化された C++コンパイラでコンパイルすると Matlab より数段高速なコードとなる。 MAT<LIB> は Matlab コンパチな C++行列クラスライブラリである。同じサイトからダウンロードできる。 W. Kahan による Matlab’s Loss is Nobody’s Gain は Matlab に対する面白いコメントである。

また、つぎのようなフリーなパッケージソフトがある。

1. Octave は MATLAB に最も近い MATLAB クローンである。 Linux 用のソフトであるが、 Windows に実装することもできる。まだ、オブジェクトプログラミングが実装されていない点では MATLAB に譲るが、線型方程式の解や固有値を計算する組み込み関数の計算速度、信頼性では MATLAB に勝るとも劣らないと思われれる。グラフィックスは gnuplot を利用している。また、 ATLAS を使って Octave を高速化することもできる (Tuning up Octave by ATLAS)。 ATLAS による高速化の効果は大変大きく、 Alpha/EV56/533MHz で、 2 つの 1000 × 1000 行列の積の計算に reference BLAS(C で書かれた高速化していない BLAS) と Octave の組み合わせでは 54.7 秒かかっていたものが、 ATLAS を使うと 4.0 秒で計算できるようになったとの報告が書かれている。

2. Scilab は typed list によりオブジェクト指向プログラミングができる環境があり、 MATLAB クローンとしては一番完成している。グラフィックスも独自の優れたものが実装されている。 MAPLE とのリンクもできる。 Linux でも Windows でも実装できる。

3. Tela は大規模数値計算のプロトタイピング言語となることを目指して作られており、グラフィックス、線形代数、 FFT などを装備して偏微分方程式を解くことを念頭にしている。

4. RLaB は行列数学用パッケージ。

5. Euler は実数、複素数、区間などを扱える Matlab クローンな言語。 Linux 版はフリーで、 Windows 版はシェアウエアである。

6. Prophet は Unix 上のライフサイエンス用言語。 7. Yorick は型チェックをしない C ライクな Linux 上の数値計算インタプリタ

言語。

8. SpYorick は Yorick のプラグインで疎行列処理をする。

9. PETSc は MPI を用いた偏微分方程式の並列数値計算用のツール。 unix 上の C や C++から呼び出して使う。

10. Ascend はフリーな強い型付けされたオブジェクト指向モデル記述言語。

11. Algae(Alki が前身) は疎系に対して Matlab より高速であることをねらった Matlab クローン的なインタプリタ数値計算言語。

著者は、 Matlab と Octave と Scilab を利用している。 Matlab 互換の Scilab による線形連立方程式の解法の例を次に示す。

-->A=[1 2;2 1]
A  =
!   1.    2. !
!   2.    1. !
-->b=[1;1]
 b  =
!   1. !
!   1. !
-->x=A\b
 x  =
!   0.3333333 !
!   0.3333333 !
-->A*x-b
 ans  =
!   0. !
!   0. !

このように非常に簡単に (かつ高速に) 線形方程式が解ける。詳細は著者による Matalab, Octave, Scilab の使い方の html ファイルなどを参照されたい (そのテキストも出版される [9])。

Demmel の著書 [10] の Matlab コードがあるのでよいプログラム例として参考にされたい。

3 演算

精度保証付き数値計算の基礎は、コンピュータ上での四則演算である。高速なハードウエアでの演算が実現されている浮動小数点演算とソフト的に演算を実行する多倍長演算に分けて議論しよう。

3.1 浮動小数点演算

現在の浮動小数点演算は IEEE754 規格に基ずくものが標準的である。これは優れた数値計算学者である Kahan らの努力により 1985 年に制定された。その内容については、当事者である W. Kahan の解説を読むのがよいであろう。制定の歴史を語った文書もある (An Interview with the Old Man of Floating-Point)。現在では Intel の Pentium CPU を含むほとんどのマシンが IEEE754 に従っているが、どのような CPU が例外かなどは Interval FAQ の記録 What machines support IEEE754? が参考になる。また、 IEEE754 にしたがってどのように浮動小数点演算が実装されているかの例は、例えば、 Intel の PentiumIII のマニュアルなどを見るとよい。

与えれた CPU が IEEE754 規格を満たしているかどうかをチェックするプログラムが存在する (TestFloat)。 TestFloat1 はソフトウエアでシミュレートされた浮動小数点システム SoftFoat と対象としているハードウエアの実行結果を比較して、違いがあるとバグの候補として出力するようなソフトウエアである。

IEEE754 によって標準化されているのは 2 進の倍精度浮動小数点数である。浮動小数点数を利用した精度保証付き計算には IEEE754 の丸め (以下に示すものの中で上と下への丸め) を利用する。 c を実数 (c ∈ ℝ) とする。

(1) 上向きの丸め (round upward) c 以上の浮動小数点数の中で最も小さい数に丸める。これを △ : ℝ → 𝔽 と表す。プログラム中では up と表すことにする。

(2) 下向きの丸め (round downward) c 以下の浮動小数点数の中で最も大きい数に丸める。これを ▽ : ℝ → 𝔽 と表す。プログラム中では down と表すことにする。

これ以外にも近似計算のデフォールトの丸めである最近点への丸め (round to nearest) とチョッピング (切り捨て)(round toward 0) の丸めがある。 IEEE754 では浮動小数点数演算 (𝔽 上での四則演算) は丸めとの関係によりつぎのように定義されている。 · ∈ {+, −, ×, /} ， ◯ ∈ {△, ▽} のとき

x ⨀ y = ◯(x · y) (任意の x，y ∈ 𝔽 について)  (1)

この式は，左辺の浮動小数点数の四則演算の結果 x ⨀ y は，右辺の数学的に正しい (実数としての) 四則演算の結果 x · y を指定された丸めを行って得られた数 ◯(x · y) に一致するように計算する規格を表している。また，平方根も

(\sqrt{x})_{fp} = ◯(\sqrt{x}) (任意の x ∈ 𝔽 について)  (2)

と，浮動小数点数演算によって計算された平方根 (\sqrt{x})_{fp} は，正確な実数演算で計算された平方根\sqrt{x} を指定された丸めの方向へ丸めた数となる規格である。注意すべきことは，指数関数や三角関数などはこのような規格をみたすようには規定されていないことである。この事に関しては初等関数などの実装の研究者 J.-M. Muller の The Table Maker’s Dilemma:our search for worst cases が参考になる。定理証明システム (HOL Light) を使って、指数関数などのマシンへの実装が IEEE754 にしたがっているかどうかをチェックするような試みもある (Floating Point Verification)。このように、初等関数の値を精度保証付きで求めるためには別途工夫が必要である。

3.2 区間演算

浮動小数点数は実数の近似であるので、一つの浮動小数点数だけでは、誤差を表現することができない。そこで、求めたい実数の下界と上界を表す 2 つの浮動小数点数のペアで誤差まで表現することが考えられた。すなわち、実数を両端を浮動小数点数とする区間で表すという考え方である。このような考え方の芽は数学においては 1900 年代の半ばころには存在していたが、具体的に区間を基礎データ型とする数値計算を唱えたのは須永 (九大) である。これを区間解析という。現代の精度保証付き数値計算はすべて区間解析をベースとすると考えてよい。その意味で精度保証付き数値計算の発祥地は日本である。これは 1960 年代のことであるが、区間解析が実用に近づくのは非常に最近であるので、数値計算技術としては 40 年近くの紆余曲折がある。これは、区間解析の基本演算である区間演算が単純でない性質をもっているからであり、これを基に効率的な精度保証数値計算を実行するのは難しいからであった。すなわち、前章までに紹介した、プロの設計した数値計算ライブリや数値計算言語の高速性と信頼性を保ちながら精度保証付き数値計算を実行する方法が見え出すまでに長い年月がかかったのである。本章では、区間解析の入門から最新の状況までを眺めてみよう。

3.2.1 区間解析の基礎

区間解析の情報を集めたサイトは区間解析の勉強に便利である。しかし、このサイトの区間解析の入門記事には須永の名前がなく Moore が区間解析のアイディアを出したと勘違いしてる。須永の業績は国際的に引用されることも多いが、まだ十分ではない。このサイトの中でインターバル計算用の数値計算ライブラリやパッケージが紹介されているが、お勧めは S.Rump の Matlab のツールボックス INTLAB である。丸めの制御、区間演算、自動スロープ計算 (微係数の近似)、線形方程式の精度保証、精度保証された初等関数と区間への拡張、非線形方程式の精度保証などがユーザフレンドリに実用的な速度で実行できる。数学的に精度の保証された初等関数を実用的にどのように実装するかのもでるとしても C++の区間演算のクラスライブラリとして PROFILE/BIAS がよい。 FORTRAN90 用区間演算ライブラリもある。

3.3 区間演算

さて、区間演算の難しい点について述べよう。以下、区間は IEEE754 に従う倍精度浮動小数点数を両端とする区間とする。区間演算は区間を基本的なデータ型として、その四則演算のことである。 2 つの区間 [a, b], [c, d] に対して例えば区間演算の和は

down; d = a + c; up; u = b + d;

と定義される。減算、乗除算も同様である。区間演算をいちいち浮動小数点演算と区別して書くのは大変なので、通常演算子多重定義によって区間型を定義して利用することが多い。ここで問題となるのは、既存の高速なパッケージでは四則演算の多重定義に対応していないことである。また、丸めの向きの変更の命令 down;,up; などの使用も考えていないことである。よって、現在の所、区間型を定義して演算子多重定義により区間演算を実行すると、つぎのような不利益が考えられる。

1. 高速パッケージを利用できなくなる。このため、自作のルーチンを使うとすると、線形方程式の解法では 20 倍前後は遅いルーチンを使うことになる。

2. 演算子多重定義によるオーバーヘッドが入り込む。

3. 丸め演算によりパイプライン処理などが崩れる可能性がある。

これらの効果により、通常の数値計算に対して、その近似解の精度を検証しようとすると、簡単に 100 倍程度の計算時間がかかるのが常であった。また、浮動小数点演算を区間演算に置き換えただけでは、区間幅が増大してしまい、ゼロを含む区間による除算などが発生して、最後まで計算が進まないことも経験された。これに対処するための様々な技法が開発され単純な区間解析から現代の区間解析に変遷してきているというのが実情である。区間幅の拡大を防ぐために、自動微分の考え方やアファイン形式の考え方などが発生した。これらによって、区間幅の増大はある程度制御できるようになったが、実行速度の向上は望めなかった。これらは非線形系への適用に関する話題であるが、つい最近になって、線形系に対しては、非常に効果のある高速化の手法が提案された。以下、章を改めてその話題を論じる。

4 数値線形代数の精度保証

著者はドイツの研究者 Rump との共同研究で、丸めの制御計算による精度保証法を開発し、線形連立方程式の数値解の精度保証が、数値解を求めるのと同じ計算量で理論的にも、実際計算 (この場合には計算時間が同程度となるという意味) でも実行できることを示した [11]。人間の計算でも、答えの検算は、もう一度やり直すのが普通であるから、数値解を求めるのと同じ計算時間で、近似解と真の解の誤差がシャープに求まるというのは理論的にもこれ以上高速化が望めないのではないかと推察される (Demmel のプレプリントによれば、厳密に条件数を計算する計算量の下限は O(n^{3})flops になることが示される)。丸めの制御計算による精度保証方式とは、区間演算を演算子多重定義することなしに実装することを考えて考案されたもので、線形演算に区間演算を制御する場合に適用できる。

基礎は内積計算である。 n次元ベクトルu = (u_{1}, u_{2}, ⋯, u_{n})^{t} とv = (v_{1}, v_{2}, ⋯, v_{n})^{t} を考える。ただし、 u_{i}, v_{i} は i = 1, 2, ⋯, n に対して倍精度浮動小数点数であるとする。このとき、 u と v の内積 z を計算することを考えよう。よく考えると

down; zd = u + v; up; zu = u + v;

とすると z ∈ [z_{d}, z_{u}] となることがわかる。ただし、以上の計算は Matlab 的に計算したものとする。 Matlab ではベクトル u, v に対して u + v はベクトルの和として演算子多重定義されており、内部では LAPACK(場合により LINPACK) のベクトルの和のライブラリ関数が呼び出されて計算される。もちろん、最高に高速化を望む場合は、 C や FORTLAN でプログラムを書いて LAPACK と速い BLAS の組み合わせで、 LAPACK の和の計算関数を呼び出して計算する。丸めの制御計算方式とは、従来の高速パッケージ関数の関数を呼び出して精度保証計算を進める方式で、従来の高速パッケージでは使用されていなかった丸めの制御命令が、関数呼び出しの前と後に位置し、パーケージの呼び出し (ひいては高速化) の妨げにならないような計算方式のことを指すことにした。以上の内積の精度保証計算が、通常の近似内積計算の 2 倍の計算速度で実行できることはプログラムを見てすぐにわかる。内積計算が高速に精度保証できると、いわゆる数値線形代数のアルゴリズムは高速に精度保証することができる。例えば A を n × n 行列、 b を n ベクトルとすると、 Ax = b の数値解の精度保証は基本的に内積計算の精度保証方式を応用して実行できる。しかし、単純のプログラムを作成すると精度保証に数値解を得る計算時間の 9 倍の手間がかかる。これを、後退誤差解析の最新の成果を利用して、改良すると、精度保証の計算時間を 1/9 倍にできることを著者は Rump との共同研究で示した。すなわち、精度保証に要する計算時間が数値解を計算する手間と同じであることが示された。

以下、少し詳細にこの結果を紹介しよう。連立一次方程式

Ax = b,  (3)

を考える。ただし、 A は n × n 実行列で x と b は n 次元ベクトルとする。 x̃ を (3) の数値解とする。このとき、つぎの定理はよく知られている。

定理 1 R を A の逆行列の近似とする。もし、

∥RA − I∥_{∞} < 1  (4)

が成立するなら、 A^{−1} が存在して、次の評価が成立する:

∥x^{∗} − x̃∥_{∞} ≤ \frac{∥R(Ax̃ − b)∥_{∞}}{1 − ∥RA − I∥_{∞}} ,  (5)

ただし、 x^{∗} は (3) の真の解とする。 □ RA − I と Ax̃ − b の包み込みを計算するプログラムは次のようになる。

down(); \underline{G} = RA − I; \underline{r} = Ac − b; up(); \overline{G} = RA − I; \overline{r} = Ac − b;

このプログラムを実行すると

RA − I ∈ [\underline{G}, \overline{G}], Ax̃ − b ∈ [\underline{r}, \overline{r}].  (6)

となる。これから ∥A^{−1}∥_{∞} と ∥x^{∗} − x̃∥_{∞} の上界を計算するプログラムは次のようになる。

Gmax = max \{|\underline{G}|, |\overline{G}|\}; rmax = max \{|\underline{r}|, |\overline{r}|\}; up();

Rnorm = max_{i=1,2,···,n} ∑_{j=1}^{n} |R|_{ij};

Gnorm = max_{i=1,2,···,n} ∑_{j=1}^{n} Gmax_{ij};

rnorm = max_{i=1,2,···,n} {rmax_{i}};

down();

D = 1 − Gnorm;

if(D > 0){ up(); Ainvnorm = \frac{Rnorm}{D}; s = Ainvnormrnorm;

}

else print(”D ≤ 0”);

このプログラムを実行したとき、 D > 0 であれば、定理 1 より

∥x^{∗} − c∥ ≤ s.  (7)

となることがわかる。 c は数値解である。

このプログラムの計算量を評価してみる。逆行列 R の数値計算 n^{3}flops かかり、 RA − I ∈ [\underline{G}, \overline{G}] という包み込みの計算に 2n^{3}flops かかるので、総合して 3n^{3}flops の計算量がかかることがわかる。これは、ガウス消去による数値解の計算量 n^{3}/3 の 9 倍の手間である。数値解を得る手間の 9 倍で精度保証できることは、従来 100 倍ぐらいはかかっていたことに比べれば、高速化されたことになるが、まだまだ精度保証のコストが高いといえる。そこで、精度保証の高速化を図ることを考えよう。そこで、次のような方針を取ることとする。逆行列 R を計算する代わりに L̃ と Ũ の逆行列 X_{L} と X_{U} を計算する。ただし、 L̃ と Ũ はガウス消去によって A を LU 分解して得られた L と U とする。そして、 R の代わりに X_{U} · X_{L} を用いる。ただし、行列の積は計算しないものとする。 X_{L} の計算量は\frac{1}{6}n^{3} flops で X_{U} の計算量も\frac{1}{6}n^{3} flops となる。以下、これ以外の精度保証の計算には O(n^{2}) の計算量の計算しか現れないことを示そう。

定理 1 から

∥x^{∗} − x̃∥_{∞} ≤ \frac{∥(X_{U} · X_{L})(Ax̃ − b)∥_{∞}}{1 − ∥(X_{U} · X_{L})A − I∥_{∞}} .  (8)

が成り立つ。ここで、 ∥(X_{U} · X_{L})(Ax̃ − b)∥ という項は行列とベクトルの積となるので O(n^{2}) flops で計算できる。高速化の基本的なアイディアは次である。 L, U, X_{L} 及び X_{U} をガウスの消去アルゴリズムで最近点への丸めのモードで計算したとする。このとき、最近得られた Higham による Wilkinson 型の後退誤差解析を使うと

∥X_{U}X_{L}A − I∥_{∞}
= 2γ_{n}∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} + γ_{n}∥ |X_{U}| |Ũ| ∥_{∞}.  (9)

という評価が得られることが示せる。これを使うと ∥X_{U}X_{L}A−I∥_{∞} の評価が O(n^{2}) flops で計算できる。

この公式を導くためには次の Higham による補題が必要となる。尚、 Higham の著書 [12] の追加の情報やその 1000 に及ぶ文献表を切り出した dvi ファイルもこの web サイトにある。

補題 1 (Higham) 最近点への丸めのモードの下で、 A ∈ ℝ^{n×n} に対してガウスの消去法で LU 分解 L̃ と Ũ を計算したとする。このとき、

L̃Ũ = A + ΔA, |ΔA| ≤ γ_{n}|L̃| |Ũ|,  (10)

が成立する。ただし、

γ_{n} = \frac{nu}{1 − nu}.  (11)

ここに、 u は単位の丸めで IEEE 754 の倍精度浮動少雨数点数に対しては 2^{−53} で与えられる。 □

補題 2 T を三角正方行列とする。 T の左逆 X_{T} をガウスの消去法で計算したとする。このとき、

|X_{T} T − I| ≤ γ_{n}|X_{T} | |T |.  (12)

が成立する。 □ この二つの補題を使うと、

∥X_{U}X_{L}A − I∥_{∞}
≤ ∥X_{U}X_{L}(A − L̃Ũ)∥_{∞} + ∥X_{U}X_{L}L̃Ũ − I∥_{∞}
≤ γ_{n}∥ |X_{U}| |X_{L}||L̃| |Ũ| ∥_{∞} + γ_{n}∥ |X_{U}| |Ũ| ∥_{∞}
+γ_{n}∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞}
= 2γ_{n}∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} + γ_{n}∥ |X_{U}| |Ũ| ∥_{∞},  (13)

と所望の公式が導出できる。ここで e = (1, 1, ⋯, 1)^{t} ∈ ℝ^{n} とすると明らかに

∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} = ∥(|X_{U}|(|X_{L}|(|L̃|(|Ũ|e))))∥_{∞}  (14)

及び

∥ |X_{U}| |Ũ| ∥_{∞} = ∥(|X_{U}|(|Ũ|e))∥_{∞}.  (15)

が成立する。これは ∥ |X_{U}| |X_{L}| |L̃| |Ũ| ∥_{∞} と ∥ |X_{U}| |Ũ| ∥_{∞} が O(n^{2}) で計算できることを示している。こうして、 ∥x^{∗} − x̃∥_{∞} が\frac{1}{3}n^{3} + O(n^{2}) flops で計算できることがわかった。

例を示そう。上述の方法により，近似解とその精度保証をする C 言語プログラムを作りその計算時間を測定した。実行時間の測定結果をまとめて表 1に示す。また，参考のために， Intel の Pentium 用に最適化された BLAS を用いた LAPACK を利用して作成したプログラムによる，同じパソコン (Intel PentiumIII 600MHz) での計算時間も示す。

固有値や得意の数値計算の精度保証の方式も丸めの制御方式で高速に実行できることも著者は示している [13]。こうして、数値線形代数に現れる問題の精度保証は、丸めの制御方式による精度保証で高速に実行できることが明らかになってきた。実装についても、ユーザの数値計算を普段実行してる環境の中で、精度保証付き数値計算が実行できることように考えられた方式なので、実用性も高いと考えられる。表 1: LU 分解を利用した近似解とその高速精度保証のための計算時間 (sec)

n C MATLAB LAPACK(最適化 BLAS 利用)
100 0  0.04  0.02
200 0.3  0.18  0.05
300 0.7  0.62  0.14
400 2  1.3  0.24
500 4  2.7  0.47
600 8  5  0.7
700 13  8  1.1
800 19  12  1.5
900 24  19  2.1
1000 37  24  2.8

5 今後の方向

精度保証付き数値計算は、数値線形代数を皮切りに実用化され始めた。その実用化の方向が着実に進むためには基本的な問題から解決しなければいけない問題も多い。以下に、そのような問題について議論しよう。

5.1 区間演算の実装

区間演算自身をより高速に実行できる環境を作るアプローチがある。区間演算を FORTRAN などの基本演算として採用するように働きかける努力の結果、 Fortran-90 は区間計算用に設計されたことは既に述べた。更に、プロポーザルの段階であるが、区間版 BLAS(INTBLAS) 開発の提案がある。

5.2 Java による数値計算

Java による数値計算に関しては、 Java をあと少し改善してくれれば、色々な意味で期待したいというのが数値計算に関わる研究者の大方の声であると思われる (The Evolution of Numerical Computing in Java, Why Java for High-Performance Computing is not an oxymoron, Dagstuhl Seminar(2000年8月): High Performance Computing and Java。 Kahan の手厳しいノートがある (How Java’s Floating-Point Hurts Everyone Everywhere)。

このような方向の Java の提案として RealJava がある。

また、現状の Java で数値計算ライブラリも開発されている。数値計算学者の G. W. Stewart の書いた Java の行列計算用パッケージ Jampack(A JAVA PACK-AGE FOR MATRIX COMPUTATIONS) がある。他の行列計算用パッケージとして JAMA : A Java Matrix Package) もあり、 Jampack とは異なる特徴をもっている。 Java BLAS の提案 ( (Draft Proposal for Java BLAS Interface) もされている。 Java を科学技術計算用に用いようとするときに参考になるサイトとして (Java for Scientific Computing)、 Java Numerics などがある。 Java LAPACK (version 0.3 beta) も公開された。

Java で精度保証付き数値計算を実行するための最低限の要求としては、 Java において丸めのモードの変更をできるようにして欲しいことが挙げられる。現在の Java の仕様では、 native method を利用して丸めの制御を行うことになる。この場合、丸めの制御部分は使用しているコンピュータや OS に依存するコードになる。

5.3 大規模な問題

精度保証付き数値計算によって、大規模な問題を解くことも可能になりつつある。本稿では、フルな係数行列をもつ線形連立方程式の直接解法による数値解の精度保証が n^{3}/3flops で可能であること示したが、構造のある係数行列の場合には、より高速な精度保証方式を作ることができる。すなわち、数値解を得る手間と同程度の計算量で精度保証するという指導原理によって、新しい精度保証方式を作り出すことが可能な分野が拡大することが期待される。

精度保証付き数値計算の応用分野も本格化し、流体系の非線形偏微分方程式の解析にも成功を収めつつある。

謝辞

本稿は、当学会に設置された精度保証付きシミュレーション研究専門委員会の活動の一環として執筆された。当専門委員会の活動を暖かくご支援頂いている理事会やご協力いただいている各位に深謝いたします。

参考文献

[1] 柏木雅英: “精度保証付き数値計算 [1]–区間解析”, シミュレーション 18 4, 261/267 (1999)

[2] 大石進一: “精度保証付き数値計算 [2]–線形方程式の高速精度保証とプログラミング技法”, シミュレーション 19 1, 39/45 (2000)

[3] 神澤雄智、相馬隆郎: “精度保証付き数値計算 [3]–常微分方程式の精度保証”, シミュレーション 19 2, 30/35 (2000)

[4] 中尾充宏: “精度保証付き数値計算 [4]–偏微分方程式の精度保証”, シミュレーション 19 3 (2000)

[5] 大石進一: 非線形解析入門、コロナ社 (2000)

[6] 中尾充宏、山本野人: 精度保証付き数値計算、日本評論社 (1998)

[7] 大石進一: 精度保証付き数値計算、コロナ社 (2000)

[8] 大石進一: 数値計算、裳華房 (1999)

[9] 大石進一: “Linux 数値計算ツール”, コロナ社 (2000.9)

[10] J. Demmel: “Applied Numerical Linear Algebra”, SIAM (1997)

[11] Shin’ichi Oishi and Siegfried M. Rump: “Fast verification of solutions of matrix equations”, submitted to Numer. Math..

[12] Nicholas J. Higham: “Accuracy and Stability of Numerical Algorithms”, SIAM (1996).

[13] Shin’ichi OISHI: “Fast Enclosure of Matrix Eigenvalues and Singular Values via Rounding Mode Controlled Computation”, to appear in Linear Algebra and its Applications