2016/03/25(金)kv-0.4.30とstiffなODE

精度保証
先日、とある人と常微分方程式の精度保証に関する情報交換をしたのですが、そのときにあるstiffな初期値問題が話題になりました。

y=(y1,y2,y3)として、

y1' = -0.04y1 + 104y2y3
y2' = 0.04y1 - 104y2y3 - 3*107y22
y3' = 3*107y22

y1(0)=1, y2(0)=0, y3(0)=0

というものです。Robertson's Problemとして有名なもののようです。係数に非常に大きな数が見えるので、確かにヤバそうです。これをkvのODE solver(ode_maffine)で解かせてみると、確かにステップ幅が非常に小さくなってしまい、全く先に進みません。次数を変えたりいろいろ試してみても、どうにもうまく行きません。stiffな奴は本気で取り組まないといけないなあと考えていたのですが、ふと思いついて以前からkvに含まれていたode_maffine2というsolverで解かせて見ると、なんと楽に解けることが判明しました。ode_maffineに比べると1000倍くらい刻み幅が大きく取れるので、結果として同じ時刻まで計算するのに1000倍くらい速いということになります。

ODE solverは、いわゆるwrappring effectを防ぐために、解軌道の初期値に関する微分を利用しています。ode_maffineは、それを得るのに与えられた微分方程式の初期値に関する変分方程式を同時に精度保証付きで解いています。ode_maffine2はそこを少し手抜きして、推進写像を低次の部分と高次の部分に分け、定次の部分の初期値に関する微分を自動微分で計算し、高次の部分は微分を計算せず直接区間評価する、という手法で計算しています。この手法はLohnerも使っているもので、新しいものではありません。普通のODEで試すとode_maffineとode_maffine2は計算時間に大差はなく、またode_maffine2は推進写像の正確なヤコビ行列を計算しないのでshooting methodに使いづらいこともあり、ode_maffineの方を主に使っていました。

stiffなODEでこのような大きな差が出た理由を推測してみます。n2の大きさの変分方程式を精度保証する必要のあるode_maffineと、nの大きさの与えられた方程式を精度保証しさえすればよいode_maffine2では、元々精度保証可能な刻み幅の限界に大きな差があった。通常のstiffでないODEでは1ステップ毎に混入する誤差をmachine epsilon以下にするための刻み幅の限界の方が精度保証限界より小さいのでそちらに制限されてしまい、差が見えない。しかし、stiffなODEになると、精度保証可能な刻み幅の限界値が小さくなり、その差が顕著に見えるようになった、と推測できます。

問題の方程式をode_maffineとode_maffine2で解かせて見ると、t=0.00390625 (1/256)までの解を精度保証するのに、ode_maffineだとおよそ3分12秒かかったのに対して、ode_maffine2だとわずか0.02秒で終了しました。そこまでのtの分割数は前者が10842なのに対して後者はわずか6でした。また、ode_maffine2だと、t=40までの精度保証をさせても、およそ3分40秒、分割数13281でした。そのときの解は、およそ次の図のような感じでした。

aaa.png


t=0からの出発直後に大きな変化があるので、あえてtを対数で表示しています。

また、このことに関連して、kvライブラリをkv-0.4.30にアップデートしました。本質的にはほとんど変わっていませんが、刻み幅の制御方法を少し変えたのと、ODEのweb demoでode_maffineとode_maffine2を選べるようにしました。example/rober.ccに今回の例題も入れておきました。
コメント（0件）

2016/03/16(水)kv-0.4.29

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

細かいバグフィックスが主で、大きな変更や機能追加はありません。

小野 寛生様のご指摘により、conv-double.hpp, conv-dd.hppで配列の外側にアクセスする可能性のある危険な書き方を直しました。ご指摘ありがとうございました。また、Visual C++では丸めの向きを変更するのに_controlfpを用いていますが、_controlfp_sが推奨されている点もご指摘いただき、そのように直しました。また、ついでにVisual C++ 2013でmpfrを用いたものを除く全てのサンプルファイルがコンパイル可能なことを確認しました。

Visual C++での丸めの向きの変更に関しては、非標準である_controlfpではなくVisual C++ 2013で追加された標準規格のfesetroundを使いたいとずっと思っていました。しかし、最近上向き丸めと下向き丸めが逆になっているというとんでもないバグを抱えていることが分かり、当分これは使えないと判断しました。

fesetroundの丸めモードが逆になる

困ったものです。interval-simpleというkvの区間演算部分だけを切り出した簡単なライブラリを公開していますが、こちらは単純化のためfesetroundを用いた実装をしていたのですが、丸めの向きが逆では致命的過ぎるのでifdefで分けてVisual C++のときは_controlfp_sを使うように直しました。

(2016/4/4に追記)

Visual Studio 2015 update2 が出たので試してみましたが、同じ結果でした。x64では相変わらず上向き丸め下向き丸めが逆になります。
コメント（2件）

2016/01/28(木)googleのAlphaGoが人類を超えた!?

その他
昨日の深夜に発表され、瞬く間に大騒ぎになったgoogleの囲碁AIがプロに勝った件、メモ代わりにまとめておきます。

最初に知ったのは、現在最強の囲碁プログラムの一つ、Crazy Stoneの作者のつぶやきでした。そこにyoutubeへのリンクが。まもなく、NHKの記事が出て、それらによると、googleがネイチャーに論文を発表し、AlphaGoと名付けられた囲碁AIが史上初めて人間のプロ棋士にハンデ無しで勝ったということでした。
  • 中身はモンテカルロ木探索にディープラーニングを組み合わせたもの
  • 中国のプロに5戦全勝
  • 3月に現在人類最強の一人と思われるイ・セドルと対局
ネイチャーの論文は
にありますが、基本的に有料です。

googleのブログの記事はこちら
開発の中心となったのは、最近Googleに買収されたイギリスのAIを研究するDeepMind社のDemis Hassabis氏。DeepMind社の記事はこちら。
この記事中に、無料で読める論文のプレプリントへのリンクがあります。

まだ論文は読んでいませんが、既存の最強の囲碁ソフトに500戦499勝とか、既存の最強の囲碁ソフトに4子おかせて勝率8割とか、恐ろしい噂が聞こえてきます。プロ棋士に勝った棋譜を見た感想もとても強いというものが多く、3月のセドルとの対局が楽しみです。
コメント（0件）

2015/12/23(水)kv-0.4.28

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

今回はライブラリ本体の変更はほとんどなく、highderiv.hpp(高階微分の計算)、test-rounding.cc(丸めモードの変更が正常に動いているか簡易チェック)を追加したくらいです。webのドキュメントは、psaやaffineにいくらか加筆しています。明日、応用数理セミナーの講演でこのライブラリの話をするのでその前に区切りとして全ての変更を放出しておきたかったという理由もあります。

test-rounding.ccの追加は、前回日記で触れたcygwinのsqrt問題があったので、とりあえず軽くチェックしようと思ったためです。いろいろな環境でテストしてみたところ、概ね以下のような感じです。
  • Linux, Macは64bit OSで64bitコンパイルしていれば問題なし
  • x86の32bitモードはFPUが使われてしまうためIEEE754に従っておらず、お勧め出来ない。
  • Windowsでは、cygwin, msys2で最適化をしないと、sqrtに丸めの向きが変わらない問題が発生。VC++の64bitなら問題なし。
sqrtはVC++の32bitモードでも問題が発生することが分かっており、VC++由来のライブラリを使っている環境ではsqrtの方向付き丸めをエミュレーションで何とかすることも考えてみます。
コメント（0件）
OK キャンセル 確認 その他