(cache) 精度保証付き数値計算（その1）

コンピュータ上で実数を表現する際には浮動小数点数を使うのですが、浮動小数点数の計算では誤差が発生します。

簡単な例を見てみます。

#include <cstdio>

int main(void) {
    float a = 0.0;
    for(int i=0; i<10000; ++i) a += 0.01;
    printf("%.10f\n", a);
}

という0.01を10000回足すプログラムを実行すると結果は100.0029525757となり、期待される100.000000000に比べて0.003ほどの誤差が発生しています。

浮動小数点数計算での誤差を抑える一番簡単な方法はfloatではなくdoubleなどのより精度の高い型を使って計算精度を上げることですが、どうしても限界はあります。
他にも問題ごとにテクニックは存在しますが、誤差を完全に無くすことはできません。
正確な計算のためには誤差がどの程度含まれているのかを定量的に評価する必要があります。
誤差を扱う技術の一つが今回扱う精度保証付き数値計算です。

浮動小数点数の丸め処理

浮動小数点数を表す際にはIEEE754という規格が使われています。
浮動小数点数は単精度なら32bit, 倍精度なら64bitの範囲内で実数を近似するので、どうしても丸める必要があります。
この規格では5種類の丸めアルゴリズムが規定されています。

最近接丸め(偶数) (Round to nearest, ties to even): 表現可能な最も近い値へ丸めます。表現可能な2つの値の中間の値の場合は、仮数の最下位ビットが0になるほうを採用します。これをデフォルトにすることが推奨されています。
最近接丸め（0から遠いほうへ） (Round to nearest, ties away from zero): 表現可能な最も近い値へ丸めます。表現可能な2つの値の中間の値の場合は、正の値ならより大きいほう、負の値ならより小さいほうの値を採用します。
Round toward 0, 0方向への丸め: 実数rについて、表現可能な値の中で絶対値がr以下でrに最も近い値を採用します
切り上げ (Rounding up, Round toward +∞): 実数rについて、表現可能な値の中でr以上で最も小さい値を採用します
切り捨て (Rounding down, Round toward −∞): 実数rについて、表現可能な値の中でr以下で最も大きい値を採用します

C99以降なら#include <fenv.h>、C++11以降なら#include <cfenv>してfesetround()を呼び出せば丸めの方法を指定できます。

区間演算

誤差を含む値を処理するときに基本となる方法は、値を表現するときに値が含まれる範囲の上限と下限を指定することです。

値が $\underline{x}$ 以上 $\overline{x}$ 以下であるときはを $[\underline{x}, \overline{x}]$ と表すことにします。
なお、簡単のために単に [x] と書いて $[\underline{x}, \overline{x}]$ を表すこともあります。

区間同士の演算を、演算後も真の値が上限〜下限の範囲に収まるように定義することを考えます。これは浮動小数点数の丸めをうまく利用することで実現できます。切り捨てを $\triangledown$ 、切り上げを $\vartriangle$ で表すと、

${\displaystyle \begin{align} [x] + [y] &= [\triangledown(\underline{x}+\underline{y}), \vartriangle(\overline{x}+\overline{y})] \\ [x] - [y] &= [\triangledown(\underline{x}-\overline{y}), \vartriangle(\overline{x}-\underline{y})] \\ [x] \times [y] &= [\text{min}\{\triangledown\underline{x}\underline{y}, \triangledown\underline{x}\overline{y}, \triangledown\overline{x}\underline{y}, \triangledown\overline{x}\overline{y}\}, \text{max}\{\vartriangle\underline{x}\underline{y}, \vartriangle\underline{x}\overline{y}, \vartriangle\overline{x}\underline{y}, \vartriangle\overline{x}\overline{y}\}] \\ [x] \div [y] &= [\text{min}\{\triangledown\underline{x}/\underline{y}, \triangledown\underline{x}/\overline{y}, \triangledown\overline{x}/\underline{y}, \triangledown\overline{x}/\overline{y}\}, \text{max}\{\vartriangle\underline{x}/\underline{y}, \vartriangle\underline{x}/\overline{y}, \vartriangle\overline{x}/\underline{y}, \vartriangle\overline{x}/\overline{y}\}] \\ \end{align} }$

とすることで区間演算についての四則演算が実現できます。

実装

これを実装したpfloatクラス（pはpreciseのp）とその利用例を以下に示します。（-std=c++11などのオプションをつけてC++11以降で動かしてください）

#include <iostream>
#include <iomanip>
#include <algorithm>
#include <cfenv>
#include <cmath>

class pfloat {
    float lb, ub; //lower bound, upper bound
public:
    pfloat(float f): lb(f), ub(f){}
    pfloat(float lb_, float ub_): lb(lb_), ub(ub_){}

    float getlb(void) {
        return lb;
    }
    float getub(void) {
        return ub;
    }
    float getcenter(void) {
        fesetround(FE_UPWARD);
        return (ub-lb)/2.0 + lb;
    }
    float getdelta(void) {
        fesetround(FE_UPWARD);
        return (ub-lb)/2.0;
    }

    friend pfloat operator+(pfloat x, pfloat y) {
        fesetround(FE_DOWNWARD);
        float l = x.lb + y.lb;
        fesetround(FE_UPWARD);
        float u = x.ub + y.ub;
        return pfloat(l, u);
    }

    friend pfloat operator-(pfloat x, pfloat y) {
        fesetround(FE_DOWNWARD);
        float l = x.lb - y.ub;
        fesetround(FE_UPWARD);
        float u = x.ub - y.lb;
        return pfloat(l, u);
    }
    friend pfloat operator*(pfloat x, pfloat y) {
        fesetround(FE_DOWNWARD);
        float l = std::min(std::min(x.lb*y.lb, x.lb*y.ub), std::min(x.ub*y.lb, x.ub*y.ub));
        fesetround(FE_UPWARD);
        float u = std::max(std::max(x.lb*y.lb, x.lb*y.ub), std::max(x.ub*y.lb, x.ub*y.ub));
        return pfloat(l, u);
    }
    friend pfloat operator/(pfloat x, pfloat y) {
        fesetround(FE_DOWNWARD);
        float l = std::min(std::min(x.lb/y.lb, x.lb/y.ub), std::min(x.ub/y.lb, x.ub/y.ub));
        fesetround(FE_UPWARD);
        float u = std::max(std::max(x.lb/y.lb, x.lb/y.ub), std::max(x.ub/y.lb, x.ub/y.ub));
        return pfloat(l, u);
    }
    friend pfloat operator+(pfloat x) {
        return x;
    }
    friend pfloat operator-(pfloat x) {
        return pfloat(-x.ub, -x.lb);
    }
    friend pfloat sqrt(pfloat x) {
        fesetround(FE_DOWNWARD);
        float l = sqrt(x.lb);
        fesetround(FE_UPWARD);
        float u = sqrt(x.ub);
        return pfloat(l, u);
    }

    friend std::ostream& operator<<(std::ostream &s, pfloat x) {
        return s << "["<< x.lb << ", " << x.ub << "]";
    }
};

using namespace std;

int main(void) {
    cout << fixed << setprecision(10);
    pfloat a(0.0);
    for(int i=0; i<10000; ++i) a = a+0.01;
    cout << a << endl;
    cout << a.getcenter() << " " << a.getdelta() << endl;
}

結果は

[99.9714889527, 100.0152282715]
99.9933624268 0.0218696595

となり、先ほどの誤差を含んだ計算について上限下限を定量的に調べることができました。

2次方程式の解を求める

そもそも精度保証付き数値計算について調べようと思ったきっかけは

「２次方程式の解の公式を使って数値解を求めてはいけない」ことを示すコードを書きました https://t.co/vC3DCr5hYq
— TM (@tmaehara) 2016年7月5日

というツイートでした。

これは「数値計算の常識」第2章に書かれている、値が近い数字同士の引き算をすると桁落ちを起こして精度が下がるという例です。

この例は $x^2+(1+\varepsilon)x+\frac{\varepsilon}{2} = 0$ という2次方程式の解を求めています。
解の公式を使うと $x = \frac{-(1+\varepsilon) \pm \sqrt{(1+\varepsilon)^2-4*\frac{\varepsilon}{2}}}{2} = \frac{-(1+\varepsilon) \pm \sqrt{1+\varepsilon^2}}{2}$ となり、足し算の方において近い数同士の引き算が発生するため、解の公式を用いた方の解は精度が低くなります。
（@tmaeharaさんによる精度保証解は-0.0000499975000000062)

これを回避するには、プログラムにあるような分子の有理化や、精度が高い方の解 x_1 を先にもとめて解と係数の関係から $x_2=\frac{c}{ax_1}$ とするようにします。

ではこのプログラムを精度保証付き数値計算で行う例を見てみます。main関数だけ示します。

int main(void) {
    float eps = 1e-4;
    pfloat a=1, b=1+eps, c=eps/2;
    pfloat D = b*b - 4*a*c;
 
	pfloat x1 = (-b + sqrt(D))/(2*a); // 解の公式を直接使うバージョン
	pfloat x2 = (2*c)/(-b - sqrt(D)); // 分子有理化したバージョン
    printf("[%e, %e]\n", x1.getlb(), x1.getub());
    printf("[%e, %e]\n", x2.getlb(), x2.getub());
}

結果は

[-5.000829e-05, -4.994869e-05]
[-4.999750e-05, -4.999749e-05]

となり、確かに後者の方が精度が高いことが分かります。

まとめと次回予告

この投稿では丸めの方向を切り替えるという単純な方法による精度保証付き数値計算の方法を示しました。
この方法は分かりやすいのですが、基本演算のたびに丸めの方向を切り替える必要がありオーバーヘッドが大きくなります。
これを回避するために、精度保証をもっと大きな計算の単位ごとに行って高速にするテクニックがいろいろあるそうなので、次回はそれらを実装したものについて書こうと思います（やる気が残っていれば）。

ちなみにkv - C++ Numerical Verification Libraries by kashiという精度保証付き数値計算のライブラリも存在するので、真面目に精度保証したい人はこういうのを使うといいのではないでしょうか。