2018/09/10(月)「精度保証付き数値計算の基礎」チュートリアル

精度保証
今年の7月に

精度保証付き数値計算の基礎

という本を出しました。私が執筆したのは、第4章「数学関数の精度保証」と、第7章「常微分方程式の精度保証付き数値解法」です。

で、今回、

「精度保証付き数値計算の基礎」チュートリアル

で話させていただきました。そこで使ったスライドを上げておきます。内容は第4章の数学関数の精度保証の話で、発表時間は20分だったので、簡単な内容です。割と多くの方に来ていただいてとても嬉しかったです。

なお、今回のチュートリアルは6章までで、7章のODEと8章のPDEについては年末に別の機会を準備中です。
コメント（0件）

2018/06/17(日)kv-0.4.44

精度保証
およそ半年ぶりに、kvライブラリを0.4.44にアップデートしました。

変更点は、
  • dd.hpp, rdd-hwround.hpp, rdd-nohwround.hppを修正し、主にddのオーバーフローに 関するバグを修正。
  • affine.hppのepsilon_reduceで、与えられたaffineベクトルが1次元だった場合に正しく 動かないバグを修正。
  • rkf78.hpp追加
  • その他細かい修正
です。

今回は、dd (double double) を使った区間演算に関する重大なバグフィックスを含んでいます。加減乗除の結果がオーバーフローするかしないかギリギリの値の場合の処理に問題があり、確率は低いものの、精度保証が正しく行われない可能性があるバグです。

以下、その詳細について説明します。ddでの区間演算は、例えば加算では真値と同じか真値よりも大きい値を計算するadd_up、真値と同じか真値よりも小さい値を計算するadd_downによって実現されています。例えばadd_downは、次のようなアルゴリズムによって実現されていました。簡単のため、python-likeに書いています。dd数(x1,x2)とdd数(y1,y2)の和を下向き丸めで計算するものです。
def dd_add_down(x1, x2, y1, y2) :
    z1, z2 = twosum(x1, y1)
    if z1 == -float("inf") :
        return z1, 0
    if z1 == float("inf") :
        return sys.float_info.max, sys.float_info.max * 2.**(-54)
    down()
    z2 = z2 + x2 + y2
    near()
    z1, z2 = twosum(z1, z2)
    return z1, z2
無誤差のtwosumはそのまま計算し、誤差の入る可能性のある「z2 + x2 + y2」の部分は下向き丸めで計算することで、下向き丸めの結果を得ています。これに加えて、赤字で示したオーバーフローの処理を行っています。ddライブラリではIEEE754と同様の無限大の取り扱いを出来るように設計しています。これは、端点に無限大を許す区間演算の実装には不可欠です。無限大は(inf,0)または(-inf,0)のような内部表現にしています。また、twosumは
def twosum(a, b) :
    x = a + b
    if (abs(a) > abs(b)):
        tmp = x - a
        y = b - tmp
    else:
        tmp = x - b
        y = a - tmp
    return x, y
のような実装になっていますが、a,bが無限大でなくa+bがオーバーフローした場合、(inf,-inf)を返すことが分かります。add_downのオーバーフロー処理は、
  • 負の方向へオーバフローした場合は、(-inf,0)を返す。
  • 正の方向へオーバーフローした場合は、下向き丸めなので無限大を返すのは誤りで、ddで表現できる正の最大数を返す。
という処理を行っていました。

自分はこれで抜けは無いと安易に考えていたのですが、研究室の修士の学生だった宮崎敏文さんにまずい場合があるのではないかと指摘され、確かに問題があることが分かりました。

問題は3つあります。
  1. 最初の問題は、「add_downで正の方向にオーバーフローするとddでの正の最大数を返すが、x1+y1は確かにddの正の最大数より大きいが、x2とy2が負の数で、全体の合計はddの正の最大数を下回ることがあるのではないか」というものです。宮崎さんの指摘はこの問題でした。
  2. 二つ目の問題は、「最初のtwosumの後は無限大の処理をしているが、二つ目のtwosumの後は何もしていない。最初のtwosumではオーバーフローしないが、二つ目のtwosumでオーバーフローする可能性があるのでは」というものです。
  3. 三つ目の問題は、「add_downの引数の片方がinf、もう片方が非infだったとき、その和はinfになるべきだが、ddの正の最大数が返ってくる」というものです。
これらの問題点を解決するため、次のように改修しました。
def dd_add_down(x1, x2, y1, y2) :
    if abs(x1) == float("inf") or abs(y1) == float("inf") :
        return x1 + y1, 0
    z1, z2 = twosum(x1, y1)
    if z1 == -float("inf") :
        return z1, 0
    if z1 == float("inf") :
        z1 = sys.float_info.max
        z2 = sys.float_info.max * 2.**(-54)
    down()
    z2 = z2 + x2 + y2
    near()
    z1, z2 = twosum(z1, z2)
    if z1 == -float("inf") :
        return z1, 0
    if z1 == float("inf") :
        z1 = sys.float_info.max
        z2 = sys.float_info.max * 2.**(-54)
    return z1, z2
すなわち、まず引数に無限大が含まれている場合を先に処理し(3.の対策)、次にtwosum後に正方向にオーバーフローして正の最大数に置き換えた場合はreturnせずにそのまま計算を継続するようにし(2.の対策)、最後に2回目のtwosumに対しても同じようにオーバーフロー処理を行なう(2.の対策)ようにしました。

これによって、次のように問題が修正されました。

1.の問題は、例えば(x_1,x_2)=(2^{1023}-2^{970},-(2^{969}-2^{916})), (y_1,y_2)=(2^{1023},-2^{969})という入力で発生します。
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
    std::cout.precision(32);
    kv::interval<kv::dd> x, y, z;

    x.lower().a1 = pow(2., 1023) - pow(2., 970);
    x.lower().a2 = -(pow(2., 969) - pow(2., 916));
    x.upper() = x.lower();
    std::cout << x << "\n";

    y.lower().a1 = pow(2., 1023);
    y.lower().a2 = -pow(2., 969);
    y.upper() = y.lower();
    std::cout << y << "\n";

    z = x + y;
    std::cout << z << "\n";
}
これを改修前のkvで計算すると
[8.9884656743115780417662938029053e+307,8.9884656743115780417662938029054e+307]
[8.9884656743115790396864485702651e+307,8.9884656743115790396864485702652e+307]
[1.797693134862315807937289714053e+308,inf]
となりますが、正しい値は1.797693134862315708145274237317049\cdots \times 10^{307}なので下端が真値より大きくなってしまっています。改修後は、
[8.9884656743115780417662938029053e+307,8.9884656743115780417662938029054e+307]
[8.9884656743115790396864485702651e+307,8.9884656743115790396864485702652e+307]
[1.797693134862315708145274237317e+308,inf]
と正しく計算されていることが分かります。

2.の問題は、例えば(x_1,x_2)=(2^{1023},2^{970}), (y_1,y_2)=(2^{1023}-2^{971},2^{969}-2^{916})という入力で発生します。
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
    std::cout.precision(32);
    kv::interval<kv::dd> x, y, z;

    x.lower().a1 = pow(2., 1023);
    x.lower().a2 = pow(2., 970);
    x.upper() = x.lower();
    std::cout << x << "\n";

    y.lower().a1 = pow(2., 1023) - pow(2., 971);
    y.lower().a2 = pow(2., 969) - pow(2., 916);
    y.upper() = y.lower();
    std::cout << y << "\n";

    z = x + y;
    std::cout << z << "\n";
}
これを実行すると、改修前は
[8.988465674311580536566680721305e+307,8.9884656743115805365666807213051e+307]
[8.9884656743115780417662938029052e+307,8.9884656743115780417662938029053e+307]
[inf,inf]
のように[inf,inf]という不正な区間が発生していましたが、改修後は
[8.988465674311580536566680721305e+307,8.9884656743115805365666807213051e+307]
[8.9884656743115780417662938029052e+307,8.9884656743115780417662938029053e+307]
[1.797693134862315807937289714053e+308,inf]
のように[ddの最大数、inf]という正しい結果を返すようになりました。

3.の問題は、区間演算では発生しないのですが、add_downを明示的に呼ぶと再現できます。
#include <kv/interval.hpp>
#include <kv/dd.hpp>
#include <kv/rdd.hpp>

int main()
{
        std::cout.precision(32);
        kv::dd x, y, z;

        x = std::numeric_limits<kv::dd>::infinity();
        std::cout << x << "\n";
        y = 0;
        std::cout << y << "\n";
        z = kv::rop<kv::dd>::add_down(x, y);
        std::cout << z << "\n";
}
これを実行すると、
inf
0
1.797693134862315807937289714053e+308
のようにinf+0の下向き丸めが正の最大数になってしまいますが、改修でちゃんと
inf
0
inf
と正常になりました。

以前からこのバグには気付いていたのですが、多忙でなかなか改修できず、モヤモヤした気分が続いていました。確率が低いとはいえ、「精度保証が保証になっていない」バグを放置してはならないと考えていますので。しかし、double-double絡みのバグが目立つのは困りものです。針の穴を通すような繊細さが要求されると感じています。
コメント（0件）

2017/11/08(水)kv-0.4.43とAVX-512による区間演算

精度保証
久しぶりにkvライブラリを0.4.43にアップデートしました。

今回の主なアップデートは2つです。一つは、
の記事で書いたライブラリをkvに入れて、誰も使ってなかったと思われるgnuplot.hppを削除しました。

もう一つが重要で、AVX-512を使った高速な区間演算の実装です。AVX-512は、Intelの最新のCPUに搭載されている機能で、現時点では、
  • Xeon Phi Knights Landing
  • Skylake-X (Core i9 7900Xなど)
  • Skylake-SP (Xeon Platinum/Gold/Silver/Bronze など)
の限られたハイエンドCPUにしか搭載されていません。2018年末に予定されているCannonlakeで一般のCPUに搭載される予定になっています。AVX-512はいわゆるSIMD命令で、現在のAVX2 (256bitレジスタ) を更に拡張してレジスタ長を512bitとし、doubleなら8つのデータを同時に処理できるようになります。

ところで、AVX-512にはひっそりと(?)、「命令そのものに丸めモード指定を含んでいるような加減乗除と平方根」が搭載されました。精度保証付き数値計算で非常に重要な区間演算には、上向き丸めや下向き丸めでの加減乗除と平方根が必要です。従来は、演算を行なう前に丸めモードの変更命令を発行しており、これが高速化の妨げになっていました。直接丸めモードを指定できる演算を使えるならば、高速化が期待できます。具体的には、intrinsic命令で
  • _mm_add_round_sd(x, y, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
  • _mm_add_round_sd(x, y, _MM_FROUND_TO_NEG_INF | _MM_FROUND_NO_EXC);
  • _mm_sub_round_sd(x, y, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
  • _mm_sub_round_sd(x, y, _MM_FROUND_TO_NEG_INF | _MM_FROUND_NO_EXC);
  • _mm_mul_round_sd(x, y, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
  • _mm_mul_round_sd(x, y, _MM_FROUND_TO_NEG_INF | _MM_FROUND_NO_EXC);
  • _mm_div_round_sd(x, y, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
  • _mm_div_round_sd(x, y, _MM_FROUND_TO_NEG_INF | _MM_FROUND_NO_EXC);
  • _mm_sqrt_round_sd(x, x, _MM_FROUND_TO_POS_INF | _MM_FROUND_NO_EXC);
  • _mm_sqrt_round_sd(x, x, _MM_FROUND_TO_NEG_INF | _MM_FROUND_NO_EXC);
の命令を使えば、加減乗除と平方根の上向き丸めと下向き丸めを直接実行できます。

新しいコンパイルオプションを設け、-DKV_USE_AVX512と付けると、これらのAVX-512命令を使って区間演算を実行するようにしました。

core i9 7900X + ubuntu 16.04 + gcc 5.4.0 という環境で、examples以下にあるtest-nishi.ccというある非線形回路の全ての解を求める計算量の多いプログラムを使って、ベンチマークをしてみました。すると、
計算精度(端点の型)コンパイルオプション計算時間
double-O310.81 sec
double-O3 -DKV_USE_AVX512 -mavx512f5.55 sec
dd-O354.10 sec
dd-O3 -DKV_USE_AVX512 -mavx512f21.42 sec
のように高速化されました。なお、kvには他にもさまざまなコンパイルオプションがあり、今回改めてそれらをまとめたページ(丸めモードの変え方とコンパイルオプションまとめ)を作成しました。このページの末尾により詳細なベンチマークのデータがあります。

AVX-512は8個のデータを並列で処理する能力がありますが、今回のプログラムでは全く並列化していません。うまく並列処理するようにプログラムを書けば更に高速化できる可能性もあります。まだAVX-512が使えるPCは身近ではありませんが、一年も経てば役立つようになるのではないでしょうか。
コメント（0件）

2017/08/23(水)kv-0.4.42

精度保証
kvライブラリを0.4.42にアップデートしました。

今回は大きな変更はありません。
  • コンパイル時にマクロKV_USE_TPFMAを定義すると、twoproductの計算にfma命令を使うようになる。
  • dka.hpp (Durand Kerner Aberth法) の重解時の安定性が向上。
  • 精度保証付き数値計算とkvライブラリの概要を説明するスライドを掲載
くらいです。概要スライドは結構手間が掛かっていますので、読んで頂ければ幸いです。概要スライドの英語版も9割がた出来ているのですが、掲載は次回以降に。
コメント（0件）

2017/08/08(火)二重振り子の精度保証付き数値計算

精度保証
少し前(5月頃)、二重振り子のシミュレーション動画がtwitterで流行したことがありました。最初のtweetがこれ。
わずかに異なる初期値をもった50本の二重振り子のシミュレーションを動画にしたもので、途中まで完全に重なっているように見える振り子が一気にバラける様子がとても面白い。

次は、それを三重振り子にしたもの。
gmpを用いて高精度計算をしており、ねとらぼで記事にもなりました。

どちらも常微分方程式の計算にはルンゲクッタ法を用いています。わずかな初期値のずれが後に大きな違いをもたらすことがとてもよく分かる動画ですが、一方で、ルンゲクッタ法で計算された値も真値とはわずかにずれており、当然その誤差も同様に後で大きな違いをもたらすことになります。だとすると、果たして意味のある計算になっているのかという疑問が生じます。

そこで、二重振り子の軌道を精度保証付きで計算し、ルンゲクッタ法とどのくらいずれるのか検証してみました。そのtweetがこれです。
このように、4次のルンゲクッタ法(刻み幅2^{-7})と精度保証付き計算とを比較してみると、t=8あたりから先は完全に違う軌道になってしまっていることが分かりました。

以下では、twitterに書けなかった細かい情報を備忘録も兼ねて書いておこうと思います。

使った式は、
\begin{align*} & (m_1+m_2)l_1\ddot{\theta_1} + m_2l_2\ddot{\theta_2}\cos(\theta_1-\theta_2) + m_2l_2\dot{\theta_2}^2\sin(\theta_1-\theta_2) +(m_1+m_2)g \sin \theta_1= 0 \\ & l_1l_2\ddot{\theta_1}\cos(\theta_1-\theta_2) + l_2^2\ddot{\theta_2} -l_1l_2\dot{\theta_1}^2\sin(\theta_1-\theta_2) + g l_2 \sin \theta_2= 0 \end{align*}
のようなものです(Wikipediaの記事と全く同じです)。舞台設定は
bbb.jpg

の通り。この式を\ddot{\theta_1}, \ddot{\theta_2}の2変数連立一次方程式と見て解いて正規形に直し、\dot{\theta_1}=\phi_1, \dot{\theta_2}=\phi_2のように置き直して1階4変数の常微分方程式とします。パラメータは
\begin{align*} m_1 &= m_2 = 1 \\ l_1 &= l_2 = 1 \\ g &= 9.8 \\ \end{align*}
とし、初期値は
\begin{align*} \theta_1(0) &= \frac{3}{4}\pi \\ \theta_2(0) &= \frac{3}{4}\pi \\ \dot{\theta_1}(0) &= 0 \\ \dot{\theta_2}(0) &= 0 \end{align*}
としました。

kvライブラリを用いて、この初期値問題を精度保証付きで計算するプログラムと、4次のルンゲクッタ法で計算するプログラムを書きました。両方のプログラムをzipで上げておきます ( doublependulum.zip )。4次のルンゲクッタ法の刻み幅は固定で2^{-7}としました。精度保証付きの方は全ての時刻の解を多項式の形で連続的に計算しているのですが、ルンゲクッタ法と同じ2^{-7}の刻み幅で密出力させています。この計算結果を動画にしたものが、先のtweetというわけです。

プログラムを走らせて得られた生データを上げておきます ( doublependulum-result.zip )。t=1のときの両プログラムのデータを抜き出すと、
\theta_1\theta_2
kv[ 0.14046540255551162,0.14046540255558421][ -0.71766989300731332,-0.71766989300724726]
rk40.1404652582264235-0.71766977132932186
となっていました。精度保証付きの方は全てのデータを区間で保持しており、この閉区間の中に真の解が存在することが数学的に保証されています。これを見ると、精度保証付き計算ではこの時点で13桁の精度を保っていますが、ルンゲクッタ法はすでに6桁程度しか精度がないことが分かります。例えば最初のtweetの動画では初期値に10^{-12}の摂動を与えていますが、t=1の時点ですでにルンゲクッタ法そのものの誤差の方が支配的になっているとするならば、本当に初期値の摂動の影響を計算していると言えるのでしょうか。

このデータをグラフにしてみました。
theta.png

t=8付近でずれが目に見えるようになり、そのあとは完全に異なる軌道を描いています。

精度保証付き計算もいつまでも高精度を保っていられるわけではなく、t=17付近を見るとわずかにグラフに幅があるのが分かります。この後一気に区間幅が爆発し、t=17.05付近で計算不能に陥ります。この精度保証付き計算プログラムはgmpのような高精度数を使用しているわけではなく、内部で使っているのは普通のdoubleです。試しに内部の数値型をdd(擬似4倍精度数)に変えてみたところ、t=30くらいまでは普通に破綻すること無く計算出来ました。

二重振り子の動画を見るだけなら、精度保証してもしなくても動きの「印象」はほとんど変わらないのですが、カオス系に対して数値シミュレーションで何らかの数学的な結論を得ようと思うならば、誤差が明確に把握できる精度保証付き数値計算は必須だと思います。みなさま、精度保証付き数値計算を使いましょう！
コメント（0件）
OK キャンセル 確認 その他