最終更新: 2016/8/7 (初版: 2013/12/26)

精度保証付き数値計算の必要性

数値計算

数値計算は、さまざまな方程式を計算機で数値的に解く技術です。我々人類は世界の仕組みを表現する方程式を発見し、それを利用することによって便利な世の中を作り上げてきました。そのような方程式のうちの多くは、 (少なくとも現時点では、恐らく将来も)解を表すきれいな式の形で解くことができません。しかしそれでもその解を知りたいという欲求に答えるため、計算機で数値の形で解く技術 (= 数値計算) を発展させてきました。

これにより、

気象観測の結果を元に明日の天気を予測する。
建物や橋を実際に作る前に十分な強度があるかどうか計算する。
自動車や飛行機などの空気の流れを実際に制作する前に計算する。

のようなことが可能になっています。

数値計算の誤差

このように、もはや現代社会に必要不可欠なものになっている数値計算ですが、さまざまな理由で誤差が混入します。

今の数値計算では、一つの実数を計算機で表現するのに有効数字約16桁(正確には2進数で有効数字53桁)の浮動小数点形式を使っています。浮動小数点形式は、

6.02 × 10²³

のようにベキ乗を使って小数点を1番左の桁と左から2番目の桁の間に移動する形式のことです。実際には10のベキでなく2のベキを使っています。もう少し詳しくいうと、

±X.XXXX × 2^YYY

の±の部分を符号部、 X.XXXXの部分を仮数部、YYYの部分を指数部といいます。計算機は、符号(±)、X.XXXX、YYYをメモリに記憶します。 X.XXXXの部分の長さが53で、log(2⁵³)≈15.95なので、約16桁というわけです。ここの長さが有限なので、計算の結果がこの長さに収まりきらない場合、 (四捨五入のような操作を行って)末尾の桁が削られてしまいます。この誤差を丸め誤差と言います。

また、解が何らかの無限反復の収束先として表されていたり、無限級数の和として表されているような場合、計算機で無限回計算することは出来ないので有限回で止めたものを解をせざるを得ません。これに起因する誤差を打ち切り誤差と言います。

更に、微分方程式など、求める解そのものが「関数」である場合を考えます。関数の形を完全に再現するには無限個の実数が必要であり、計算機は無限個の実数を記憶することはできないので、何らかの手段で問題を有限個の実数を求める問題に置き換えます。これを離散化といいます。このとき、元の問題の解と離散化された問題の解にはどうしてもズレが生じます。この誤差を離散化誤差と言います。

精度保証付き数値計算とは

従来、計算結果に混入した誤差がどのくらいの大きさなのか、大体の推定をすることはあっても、数学的に厳密に評価するというようなことは行われていませんでした。精度保証付き数値計算とは、前述した丸め誤差、打ち切り誤差、離散化誤差などのすべての誤差を数学的に厳密に把握するような数値計算法の総称です。

簡単な例を示してみましょう。2次方程式

a x² + b x + c = 0

で、a=1, b=10¹⁵, c=10¹⁴として、解のうち大きい方

(-b + sqrt(b² - 4ac)) / (2a)

を計算する問題を考えます。また同時に分子を有理化した式

(2c) / (-b - sqrt(b² - 4ac))

でも計算してみます。まずは通常の近似計算。

#include <iostream>
#include <cmath>

int main()
{
    double a, b, c;

    std::cout.precision(17);

    a = 1.;
    b = 1e15;
    c = 1e14;

    std::cout << (-b + sqrt(b * b - 4. * a * c)) / (2. * a) << "\n";
    std::cout << (2 * c) / (-b - sqrt(b * b - 4. * a * c)) << "\n";
}

計算結果は、

-0.125
-0.10000000000000002

のようになりました。数学的には同じになるはずなのに、大きく異なる計算結果が得られました。

次に、区間演算と呼ばれる丸め誤差を把握する技術を使って、この計算に混入した丸め誤差を評価してみます。 kvライブラリはC++で書かれた精度保証付き数値計算のためのライブラリで、その区間演算機能を使ってみます。演算子多重定義を利用して、わずかな変更で区間演算ができます。

#include <iostream>
#include <cmath>
#include <kv/interval.hpp>
#include <kv/rdouble.hpp>

int main()
{
    kv::interval<double> a, b, c;

    std::cout.precision(17);

    a = 1.;
    b = 1e15;
    c = 1e14;

    std::cout << (-b + sqrt(b * b - 4. * a * c)) / (2. * a) << "\n";
    std::cout << (2 * c) / (-b - sqrt(b * b - 4. * a * c)) << "\n";
}

実行結果は、

[-0.1875,-0.0625]
[-0.10000000000000004,-0.099999999999999991]

となりました。これは、この閉区間に真の値が含まれていることを表しています。区間の幅を見れば、前者の計算は大きな誤差を含んだ計算で、後者の計算は誤差が小さいことが分かります。

2次方程式は誤差混入の原理も明らかになっていて大きな誤差を回避する方法もよく知られています。しかし、例えば最新のスーパーコンピュータでは1秒に10¹⁷回もの計算を行なうことができます。多くの場合数値計算では「前の計算結果を使って次の計算を行なう」ので、一度どこかで大きな誤差が混入してしまうとそれ以降の計算は信頼できないものになってしまいます。大量の計算を行なう過程のどこかで致命的な誤差の混入が発生していないことを誰が保証してくれるのでしょうか? これからますます「計算の速度」のみならず「計算の品質」が重要になってくると思っています。

Lorenz方程式

最後に、有名なLorenz方程式の軌道を精度保証付きで計算してみます。

Lorenz方程式は、カオス的な振る舞いを示す常微分方程式として有名です。

x' = 10(y - x)
y' = 28x - y - xz
z' = (-8/3)z + xy

これを、初期値 (x,y,z) =(15, 15, 36) として、

matlab
精度保証付き数値計算

で計算し比較してみます。

matlabの方は、 lorenz.mを作成し、

options=odeset('RelTol', 1e-15)
[T, X] = ode45(@lorenz, [0 28.5], [15 15 36], options);
dlmwrite('lorenz-matlab.txt', horzcat(T, X))

のように最も標準的なode45を用いて最高精度(10^-15)を指定しました。

精度保証付き数値計算の方は、 verify.ccを作成しこれを kvライブラリと共にコンパイルしたものを用いました。
(実際にはプロットさせるためのコードを挿入している。)

計算した解を、時刻の幅6ずつの区間に区切って観察してみます。

t=6まででは、両者の軌道はぴったり重なっていて区別がつきません。精度保証の方は真の解の上限と下限を示す2本の軌道を描いているのですが、目でみることはできません。

t=12まででも同様。両者はぴったり重なっています。

このあたりから様子が違ってきます。t=13-14あたりからわずかにずれが認められます。 t=15あたりではっきりとずれが認められ、これ以降はまるで異なる軌道を描くことになります。

このあたりは全くずれています。

精度保証の方もいつまでも計算できるわけではなく、t=27くらいから軌道の上限と下限が分離し、幅を持った解になっていることが見えるようになります。

世間で信頼されているソフトウェアに最高精度を指定しても、このように何の警告も出さずに真の解とは全く違った軌道を出力することがあります。精度保証付き数値計算では、黙って嘘を付くことは決して無く、計算が無理な場合は区間の幅が広くなることによって警告を発してくれます。

全ての数値計算を精度保証付き数値計算に置き換えるべく、我々は日夜努力しています。