精度保証付き数値計算の原理と実装

精度保証付き数値計算の原理と実装
柏木 雅英
[email protected]
http://verifiedby.me/
早稲田大学 基幹理工学部 応用数理学科
第 153 回リスク工学セミナー@筑波大学
(2016 年 11 月 21 日)
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
1 / 72
発表の概要
数値計算の誤差
精度保証付き数値計算について
区間演算とその実装
Krawczyk 法
KV ライブラリの紹介
概要
区間演算
自動微分、4 倍精度演算、mpfr、affine arithmeric
ベキ級数演算 (psa)
数値積分
常微分方程式
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
2 / 72
数値計算
方程式は解けない
人類は現象を方程式で記述するという方法論を発見し、それを利用して
様々な技術を生み出したり未来を予測したりできるようになった。しか
し、ほとんどの方程式はきれいな式の形で解くことはできない。
数値計算
解析的に解くことが不可能な数学上の問題を計算機で数値的に解く手法
のこと。
コンピュータの発展とともに様々な分野で用いられ、現代社会に欠
かせない技術である。
気象予報、建築物の強度解析、自動車の衝突、化学反応、流体などの
シミュレーション。
天体、ロケットなどの軌道計算。
コンピュータグラフィックス。
世界最初のコンピュータと呼ばれる ENIAC は、ミサイルの弾道の数
値計算を目的に開発された。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
3 / 72
計算機の実数の表現 (IEEE754 倍精度浮動小数点数)
実数の表現
6.02 × 1023 のような、「1 に近い数 × ベキ数」の形で記憶する。浮
動小数点形式という。
2 進数で記憶し、64bit を使う。
64bit を
s
e10
···
e9
e0
m51
m50
···
m0
のように分割し、それは
(−1)s (1 +
51
∑
∑10
mi 2i−52 ) × 2(
i
i=0 ei 2 )−1023
i=0
で計算される数値と対応する。
計算精度 (有効数字) は、2−52 ≃ 2.22 × 10−16 ≃ 10−15.65 なのでおよ
そ 16 桁。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
4 / 72
丸め誤差
丸め誤差
有効数字がおよそ 16 桁なので、例えば 1 ÷ 3 を計算すると
1 ÷ 3 = 0.3333333333333333148 . . .
のように少し誤差が入る。これを丸め誤差という。
2 進法を使っているので、1 ÷ 10 = 0.1 でも
1 ÷ 10 = 0.100000000000000005551 . . .
のように誤差が入る。
1 回だけの計算なら大したことないが、誤差が積み重なると想像以
上に大きく誤差が拡大することがある。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
5 / 72
2 次方程式の丸め誤差
2 次方程式
2 次方程式 ax 2 + bx + c = 0 の解の公式は、
√
−b ± b 2 − 4ac
x=
2a
大きな誤差が見られる例
a = 1, b = 1015 , c = 1014
のときの大きい方の解を計算する。
√
−b + b 2 − 4ac
(解の公式)
= −0.125
2a
(解の公式の分子を有理化)
2c
√
= −0.10000000000000002
−b − b 2 − 4ac
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
6 / 72
連立一次方程式の誤差
連立一次方程式
(
64919121 −159018721
41869520.5 −102558961
真の解
(
x
y
)(
x
y
)
(
=
1
0
)
ガウスの消去法で計算した解
)
(
=
205117922
83739041
)
(
x
y
)
(
=
106018308.0071325
43281793.001783125
)
間違った解の残差を計算すると · · ·
(
64919121 −159018721
41869520.5 −102558961
柏木 雅英 (早稲田大学)
)(
106018308.0071325
43281793.001783125
精度保証付き数値計算の原理と実装
)
(
−
1
0
)
(
=
リスク工学セミナー
0
0
)
7 / 72
湾岸戦争におけるパトリオットミサイルの事故
http://www.sydrose.com/case100/298/
http://www.fas.org/spp/starwars/gao/im92026.htm
湾岸戦争中、イラク軍のスカッド・ミサイル (Scud Missile) の迎撃のた
め、アメリカ軍のパトリオット・ミサイル (Patriot Missile) が発射された
が、その内部コンピュータのわずかな誤差からスカッド・ミサイルを捉
えることが出来ず迎撃は失敗した。
発生日時 1991 年 2 月 25 日
発生場所 サウジアラビア、Dharan
死者 28 名、負傷者約 100 名
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
8 / 72
精度保証付き数値計算
精度保証付き数値計算
計算をすると同時に計算結果の数学的に厳密な誤差評価をも計算する数
値計算法
必要な技術
区間演算 (切り捨てと切り上げを併用して計算に混入した誤差を把
握する。関数の像の評価も行う。)
不動点定理 (解の存在する十分条件を区間演算で確認して、方程式
の解の存在と存在範囲を保証する)
自動微分も重要。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
9 / 72
区間演算 (1)
区間演算
精度保証付き数値計算の基本技術。
数値を、計算機で表現可能な浮動小数点数を両端に持つ閉区間
X = [a, b] で表現する。
区間同士の演算は、集合値演算として計算結果として有り得る値を
包含するように行う。
IEEE754 標準でも定義されている「方向付き丸め」を利用する。
「丸め誤差の把握」と「関数の値域の評価」の 2 つの役割がある。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
10 / 72
区間演算 (2)
X = [a, b], Y = [c, d]
·, · はそれぞれ下向き丸め、上向き丸め
加算 X + Y = [a+c, b+d]
減算 X − Y = [a−d, b−c]
乗算 X × Y =
d ≤0
c < 0, d > 0
b≤0
[b×d, a×c]
[a×d, a×c]
a < 0, b > 0 [b×c, a×c] [min(a×d, b×c), max(a×c, b×d)]
a≥0
[b×c, a×d]
[b×c, b×d]
Y <0
Y >0
X < 0 [b/c, a/d] [a/c, b/d]
除算 X /Y =
X ∋ 0 [b/d, a/d] [a/c, b/c]
X > 0 [b/d, a/c] [a/d, b/c]
√
√ √
X = [ a, b]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
c≥0
[a×d, b×c]
[a×d, b×d]
[a×c, b×d]
リスク工学セミナー
11 / 72
区間演算 (3)
10 進数有効数字 3 桁で (1 ÷ 3) × 3 を計算すると · · ·
・普通の数値計算
・区間演算
最初は幅のない区間からスタート。
1 ÷ 3 = 0.333
[1, 1] ÷ [3, 3] = [0.333, 0.334]
0.333 × 3 = 0.999
[0.333, 0.334] × [3, 3] = [0.999, 1.01]
四捨五入による誤差がでる。
柏木 雅英 (早稲田大学)
常に区間内に真の値を含みながら計
算が進む。
精度保証付き数値計算の原理と実装
リスク工学セミナー
12 / 72
2 次方程式の例を区間演算で計算する
2 次方程式
2 次方程式 ax 2 + bx + c = 0 の解の公式は、
√
−b ± b 2 − 4ac
x=
2a
大きな誤差が見られる例
a = 1, b = 1015 , c = 1014
のときの大きい方の解を計算する。
√
−b + b 2 − 4ac
(解の公式)
= [−0.1875, −0.0625]
2a
2c
√
(解の公式の分子を有理化)
=
−b − b 2 − 4ac
[−0.10000000000000004, −0.099999999999999991]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
13 / 72
区間演算 (4)
丸めの向きの変更方法
IEEE754 で丸めの向きを変更可能にすることが要請されているが、
その方法は CPU やコンパイラによって違う。
X86 だと FPU と SSE2 の 2 種類の演算器で丸めの向きの変更方法が
異なるなど複雑。
最近は C99 準拠のコンパイラが増えたので、fenv.h と fesetround
を使えば簡単になった。
#i n c l u d e <i o s t r e a m >
#i n c l u d e <f e n v . h>
i n t main ( )
{
d o u b l e x =1 , y =10 , z ;
std : : cout . p r e c i s i o n (17) ;
f e s e t r o u n d (FE TONEAREST) ;
z = x / y;
s t d : : c o u t << z << s t d : : e n d l ;
f e s e t r o u n d (FE DOWNWARD) ;
z = x / y;
柏木 雅英 (早稲田大学)
s t d : : c o u t << z << s t d : : e n d l ;
f e s e t r o u n d (FE UPWARD) ;
z = x / y;
s t d : : c o u t << z << s t d : : e n d l ;
}
0.10000000000000001
0.099999999999999991
0.10000000000000001
精度保証付き数値計算の原理と実装
リスク工学セミナー
14 / 72
区間演算 (5)
自分で実装するときに気をつけるべき注意点
(入力誤差) プログラム中の定数 (double x = 0.1; の「0.1」など)
はコンパイル時に 10 進数 →2 進数変換が行われ、その丸めの向きは
制御できない。
(出力誤差) 同様に、std::cout << x << std::endl; などの表示
のときの 2 進数 →10 進数変換の丸めの向きも制御できない。
(コンパイラの最適化) 最適化によって演算の順序が変えられたり定
数をコンパイル時に計算されてしまうなどなどして、-O3 などで最
適化を強くかけると丸めの向きの変更が効かないことがある。=⇒
volatile を使うなどして適切に最適化を抑制する。
(数学関数) 数学関数で精度が保証されており丸めの向きの変更も出
来るのは sqrt のみであり、sqrt 以外の関数は全く信用できないの
で、自作する必要がある。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
15 / 72
区間演算 (6)
区間演算の過大評価
区間演算は、確かに真の値を含むものの、想像以上に区間幅が広がることがある。
f (x) = x 2 − 2x,
x ∈ [0.9, 1.1]
区間演算の結果は [−1.39, −0.59] となるが、真の像は [−1, −0.99]。
過大評価を緩和する手段:
平均値形式
f (I ) を直接評価するより、c = mid(I ) として
f (c) + f ′ (I )(I − c)
を計算した方がよい場合が多い。
affine arithmetic
(後述)
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
16 / 72
Krawczyk 法 (Krawczyk(1969). Kahan(1968) が先?)
区間演算を用いて、非線形方程式 f (x) = 0,
証する。
f : Rn → Rn の解の存在を保
Krawczyk 法
I ⊂ Rn は区間ベクトル (候補者集合)、c = mid(I )、R ≃ f ′ (c)−1 、E は単位
行列とし、
K (I ) = c − Rf (c) + (E − Rf ′ (I ))(I − c)
としたとき、K (I ) ⊂ int(I ) ならば I に f (x) = 0 の解が唯一存在する。
証明
g (x) = x − Rf (x) に対して平均値形式と縮小写像原理を適用する。
f ′ (I ) (区間 I におけるヤコビ行列を包含する区間行列) が必要 =⇒ 自動微分
の必要性
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
17 / 72
Krawczyk 法の例
例題
{
x12 + x22 − 1 = 0
の解が区間 I =
x1 − x2 = 0
計算例
(
c=
f ′ (I ) =
0.7
0.7
(
)
(
[0.6, 0.8]
[0.6, 0.8]
)
にあることを示す。
(
)
0.4 0.5
、R =
≃ f ′ (c)−1 とすると、
0.4 −0.5
)
[1.2, 1.6] [1.2, 1.6]
となり、
1
−1
(
′
c − Rf (c) + (E − Rf (I ))(I − c) =
[0.68, 0.736]
[0.68, 0.736]
)
(
⊂
[0.6, 0.8]
[0.6, 0.8]
)
なので、解の一意存在が保証された。
機械的にできるのが嬉しい!
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
18 / 72
Krawczyk 法の応用
近似解 c を元にした解の存在保証
R ≃ f ′ (c)−1 、r = 2∥Rf (c)∥ (Newton 法の修正量の 2 倍) とし、


[−1, 1]


..
I =c +r

.
[−1, 1]
を候補者集合として Krawczyk 法を使うとよい。
区間ベクトル I 内の全解探索
区間ベクトル I での解の存在定理 (Krawczyk 法)
区間ベクトル I での解の非存在定理 (例えば f (I ) ̸∋ 0 なら I に解なし)
の 2 つの定理を両方試し、両方共失敗したならば区間を分割する、という作
業をどちらかが成立するまで再帰的に繰り返す。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
19 / 72
KV ライブラリ
http://verifiedby.me/kv/ で公開中。
作成開始は 2007 年秋頃。公開開始は 2013 年 9 月 18 日。
言語は C++。boost C++ Libraries (http://www.boost.org/) も
必要。
全てヘッダファイルで記述されており、インストールはヘッダファ
イルをどこかに配置するだけ。
計算に使う数値の型を double 以外の型に容易に変更することが出
来る。
オープンソースである。精度保証付き数値計算の結果が「証明」で
あると主張するならば、精度保証付き数値計算のプログラムが公開
されていないという状態はありえない。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
20 / 72
kv - a C++ Library for Verified Numerical Computation
http://verifiedby.me/kv/
kv - a C++ Library for Verified Numerical Computation
http://verifiedby.me/kv/
5. 区間演算(interval)
English | 日本語
最終更新: 2016/11/20
6. 4倍精度演算(dd)
kv - C++による精度保証付き数値計算ライブラリ
7. MPFRラッパー
8. 複素数演算(complex)
柏木 雅英
9. 自動微分(autodif)
1. はじめに
10. Affine Arithmetic (affine)
本ページでは、精度保証付き数値計算を行うためにC++で作成した ライブラリ群を公開している。
11. ベキ級数演算 (psa)
特に非線形計算の精度保証を行うとき、template機能によって 複雑な数値型をすっきり記述でき、なおかつ "zero-overhead principle" で 計算速度が遅くならない
C++は、非常に適していると言える (ほぼ唯一無二であると作者は考えている。)。
12. Krawczyk法による非線形方程式の解の精度保証
2007年秋頃〜2013年春頃の間は、区間演算を行うのにboostに含まれている intervalライブラリを用いて開発していたが、 boost.intervalは残念ながら不完全な部
分が多く ライブラリ本体に手を入れざるを得なかった。 boost全体がアップデートする度にinterval部分にpatchを当てるのも面倒になって きたので、interval部分は全
て自作することにした。 本ライブラリはboost.intervalは使っていないが、線形計算を行うboost.ublasなど、 部分的にまだboostを使っている。
13. 非線形方程式の全解探索
boost.intervalを使っていた頃の古い情報はもうアップデートしないが、 一応 ここに保存しておく。
14. 常微分方程式の初期値問題の精度保証
2. 動作環境
15. 初期値問題solverと射撃法による境界値問題の精度保証
C++とboostが動くことが必要。boostが動かないといけないので、あまり古い コンパイラでは動作しないだろう。
16. 数値積分
区間演算を実現するために丸めモードの変更を行っているので、 CPUとコンパイラには制限がある。 詳細は 区間演算(interval) の項を参照。
17. 特殊関数の精度保証
また、 4倍精度演算 を行うときにIntelのFPUは いろいろ問題を抱えており、一応対策しているが、特に32bitモードでは問題が 発生する可能性があることを注意してお
く。
18. その他の機能
一応、次の環境で動作を確認したことがあるが、主に開発は ubuntu 14.04 64bitで行っており、その他は コンパイルが通るかチェックする程度である。
19. 関数オブジェクトによる問題の記述
ubuntu 16.04 (64bit) + gcc 5.3
ubuntu 16.04 (64bit) + clang 3.8
ubuntu 14.04 (64bit) + gcc 4.8
ubuntu 10.04 (64bit) + gcc 4.4
ubuntu 10.04 (32bit) + gcc 4.4
windows7 (64bit) + Visual Studio 2015
windows7 (64bit) + Visual Studio 2013
windows7 (64bit) + Visual Studio 2008
windows7 (64bit) + cygwin + gcc
Windows7 (64bit) + MSYS2 (64bit) + gcc
Mac OS X snow leopard + gcc
Mac OS X Yosemite + gcc
ubuntu 12.04 on MK802(ARM) + gcc
raspberry pi + raspbian + gcc
Sharp NetWalker PC-Z1 + ubuntu 9.04 + gcc
Intel Edison + gcc
20. その他
20.1 boostとは
20.2 行列計算 (boost.ublas)
20.3 数値型のプログラミング
20.4 おまけ: シンプルな区間演算ライブラリ
21. kvライブラリのwebデモ
kvライブラリをweb上で試せるデモ。
3. ダウンロードとインストール
22. おわりに
ダウンロード: kv-0.4.37.tar.gz (2016年11月20日公開)
(古いversionはこちら)
本ライブラリは、 HIKMOT のC++パートで使われ、その高速性に大きく寄与しています。 また、当研究室の学生の日々の研究にも使われ、 その意見を反映しながら開
発を進めています。
ヘッダファイルのみで動くように作られている。よって、ライブラリをmakeする 等のインストール作業は必要ない。 archiveを展開するとkv, test, exampleの3つの
directoryが作られるが、 本体は kv 以下。kv以下をどこか (current directoryでも/usr/local/includeでも) に置いておくだけで良い。 動作確認は、kv及びboostが
include pathに入った状態でtest以下または example以下の適当な.ccファイルをコンパイル出来ればOK。 例えば、
大勢の人に使われて鍛えられないとライブラリは成長しないので、 なるべく多くの方に使って頂いてご意見を頂ければありがたいです。
本ライブラリの開発には、NTT未来ねっと研究所の柏木啓一郎氏の多大なる協力を 得ています。ここに感謝の意を表します。
単にarchiveを展開した状態
boostは/usr/local/includeにある
(2015年10月5日追記) このソフトウェアはMITライセンスに基づいて公開されています。
更新履歴
ならば、例えばtest以下で
c++ -I.. -I/usr/local/include test-interval.cc
開発にご協力下さった方々
とやってエラーが出なければ問題ない。 (「kv及びboost directoryが置いてあるdirectory」を指定することに注意)
kv - a C++ Library for Verified Numerical Computation /
コンパイルオプションは、-O3等で最適化を最大にし、 -DNDEBUGを付けることを推奨。 どちらも実行速度に大きく影響する。 (NDEBUGマクロの意味は boost.ublas
の項を参照。) また、IntelのCPUで64bit OSの場合は、-DKV_FASTROUND を付けると区間演算が速くなる。
([2015/12/17]に追加) 新しい環境やコンパイルオプションを使うときは、test/test-rounding.ccをコンパイル/実行して、 加減乗除と平方根で丸めの向きがちゃんと
変わっているか簡易チェックするとよい。
kvライブラリが提供する機能は全てkv名前空間の中にあり、他のライブラリと ぶつからないように配慮されている。
4. 構成ファイルの役割一覧
1/2
2016年11⽉20⽇ 21:54
柏木 雅英 (早稲田大学)
2/2
精度保証付き数値計算の原理と実装
2016年11⽉20⽇ 21:54
リスク工学セミナー
21 / 72
KV ライブラリの主な機能
数値型
区間演算 (多数の数学関数
含む)
複素数演算
4 倍精度 (double-double) 演算
affine arithmeric
MPFR ラッパー
ベキ級数演算
自動微分
アプリケーション
Krawczyk 法による非線形方程式の精度保証
非線形方程式の全解探索
常微分方程式の初期値問題
常微分方程式の境界値問題
数値積分
特殊関数
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
22 / 72
なぜ C++を選んだか?
y = (x+1) * (x-2) + log(x);
同じ表記の数式 (プログラム) に対して、
double
interval
自動微分型
内部が interval な自動微分型
ベキ級数型
多倍長数
多倍長数 interval
etc
など、様々な特殊な動作をする「数値型」を「流し込む」ことが
多い。
python, ruby, matlab などの「型の無い」言語を使って記述すると楽
だが、実行時に演算を行う度に内部では型の判定が行われることに
なり、速度が低下する。
C++の template 機能を使えば、型を仮定しない generic な記述を行
いながら、実行時ではなくコンパイル時に型の判定を全て終わらせ
るため、速度が低下しない。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
23 / 72
C++のテンプレート機能 (テンプレート関数)
テンプレート無し
#i n c l u d e <i o s t r e a m >
v o i d swap ( i n t& a , i n t& b ) {
i n t tmp ;
}
tmp = a ;
a = b;
b = tmp ;
テンプレートあり
t e m p l a t e <c l a s s T> v o i d swap (T& a , T& b ) {
T tmp ;
}
v o i d swap ( d o u b l e& a , d o u b l e& b ) {
d o u b l e tmp ;
tmp = a ;
a = b;
b = tmp ;
}
tmp = a ;
a = b;
b = tmp ;
i n t main ( )
{
i n t a =1 , b =2;
}
i n t main ( )
{
i n t a =1 , b =2;
swap ( a , b ) ;
s t d : : c o u t << a << ” ” << b << ”\n ” ;
swap ( a , b ) ;
s t d : : c o u t << a << ” ” << b << ”\n ” ;
d o u b l e x =1. , y =2.;
swap ( x , y ) ;
s t d : : c o u t << x << ” ” << y << ”\n ” ;
d o u b l e x =1 . , y = 2 . ;
}
swap ( x , y ) ;
s t d : : c o u t << x << ” ” << y << ”\n ” ;
柏木 雅英 (早稲田大学)
#i n c l u d e <i o s t r e a m >
精度保証付き数値計算の原理と実装
リスク工学セミナー
24 / 72
C++のテンプレート機能 (テンプレートクラス)
テンプレート無し
#i n c l u d e <i o s t r e a m >
class pair int {
int a , b;
public :
p a i r i n t ( i n t x , i n t y ) : a ( x ) , b ( y ) {}
};
void print () {
s t d : : c o u t << a << ” ” << b << ”\n ” ;
}
t e m p l a t e <c l a s s T> c l a s s
T a, b;
public :
pair {
void print () {
s t d : : c o u t << a << ” ” << b << ”\n ” ;
void print () {
s t d : : c o u t << a << ” ” << b << ”\n ” ;
}
柏木 雅英 (早稲田大学)
テンプレートあり
#i n c l u d e <i o s t r e a m >
p a i r (T x , T y ) : a ( x ) , b ( y ) {}
pa i r do u b l e ( double x , double y ) : a ( x ) , b
( y ) {}
i n t main ( )
{
p a i r i n t p (1 , 2) ;
p . print () ;
class pair double {
double a , b ;
public :
};
pair double q(1. , 2.) ;
q . print () ;
}
};
}
i n t main ( )
{
p a i r <i n t > p ( 1 , 2 ) ;
p . print () ;
p a i r <d o u b l e> q ( 1 . , 2 . ) ;
q . print () ;
}
精度保証付き数値計算の原理と実装
リスク工学セミナー
25 / 72
行列ベクトル計算
boost.ublas
行列ベクトル計算は、boost (http://www.boost.org/) に含まれている
ublas を用いている。
ublas は、行列やベクトルの成分の型がテンプレートになっているので、区
間行列等が自然に扱える。
名前は ublas だが、BLAS の機能を全て持っているという意味で、BLAS 的な
高速性を持つわけではない。
KV ライブラリにおける線形計算
線形計算においては、例えば行列積 C = A × B を
(1) 丸めの向きを下向きに変更してから C = A × B を計算
(2) 丸めの向きを上向きに変更してから C = A × B を計算
のような手順で計算することによって、丸めの向きの変更回数を減らし高速
な BLAS を利用できる。
KV ライブラリでは double 以外の型を自然に利用できることを重視したた
め、現在の version ではこのような技術は全く使われていない
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
26 / 72
区間演算 (interval)
上端下端型の区間演算を行う。
exp, log, sin, cos, tan, sinh, cosh, tanh, asin, acos, atan, asinh, acosh,
atanh, expm1, log1p, abs, pow などの精度保証付きの数学関数を
持つ。
10 進文字列との丸め方向指定付き相互変換を持ち、正しく入出力が
出来る。
上端と下端に用いる数値型はテンプレートになっており、double 以
外の型も使える。例えば double-double 型や MPFR を使える。ただ
し、上向き下向き双方の丸めに対応した加減乗除、平方根、文字列
との相互変換の方法を定義する必要がある。
サポートする環境は、C99 準拠の fesetround が使えること。x86 の場
合のより高速なオプションや、丸めの変更を全く行わないオプショ
ンもある。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
27 / 72
区間演算プログラムの例
#i n c l u d e <kv / i n t e r v a l . hpp> // 区間演算
#i n c l u d e <kv / r d o u b l e . hpp> // d o u b l e の方向付き丸めを定義
i n t main ( )
{
kv : : i n t e r v a l <d o u b l e> s , x ;
std : : cout . p r e c i s i o n (17) ;
s = 0;
f o r ( i n t i =1; i <=1000; i ++) {
x = i;
s += 1/ x ;
}
s t d : : c o u t << s << ”\n ” ;
}
[7.485470860549956 ,7.4854708605508238]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
28 / 72
使い方
解凍
$ ls
kv − 0 . 4 . 3 7 . t a r . gz
$ t a r x f z kv − 0 . 4 . 3 7 . t a r . gz
$ ls
kv −0.4.37/ kv − 0 . 4 . 3 7 . t a r . gz
$ cd kv −0.4.37
$ ls
LICENSE . t x t README. t x t
example
kv
test
必要なのは kv ディレクトリ以下。適当な場所に配置する。
compile & run
$ ls
i n t e r v a l . c c kv /
$ c++ −I . −O3 i n t e r v a l . c c
$ . / a . out
[7.485470860549956 ,7.4854708605508238]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
29 / 72
区間演算プログラム (double-double)
#i n c l u d e <kv / i n t e r v a l . hpp> // 区間演算
#i n c l u d e <kv / dd . hpp> // d o u b l e−d o u b l e
#i n c l u d e <kv / r d d . hpp> // dd の方向付き丸めを定義
i n t main ( )
{
kv : : i n t e r v a l <kv : : dd> s , x ;
std : : cout . p r e c i s i o n (34) ;
s = 0;
f o r ( i n t i =1; i <=1000; i ++) {
x = i;
s += 1/ x ;
}
s t d : : c o u t << s << ”\n ” ;
}
[7.485470860550344912656518204308257 ,7.485470860550344912656518204360964]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
30 / 72
Krawczyk 法による解の精度保証
#i n c l u d e <kv / kraw−a p p r o x . hpp>
namespace ub = b o o s t : : n u m e r i c : : u b l a s ;
s t r u c t Func {
t e m p l a t e <c l a s s T> ub : : v e c t o r <T> o p e r a t o r ( ) ( c o n s t ub : : v e c t o r <T>& x ) {
ub : : v e c t o r <T> y ( 2 ) ;
y (0) = x (0) ∗ x (0) − x (1) − 1 . ;
y (1) = ( x (0) − 2.) ∗ ( x (0) − 2.) − x (1) − 1 . ;
return y ;
}
};
i n t main ( ) {
ub : : v e c t o r <d o u b l e> x ;
ub : : v e c t o r < kv : : i n t e r v a l <d o u b l e> > i x ;
std : : cout . p r e c i s i o n (17) ;
x . r e s i z e (2) ;
x ( 0 ) = 1 . 0 1 ; x ( 1 ) = 0 . 0 1 ; // ( 1 . 0 1 , 0 . 0 1 ) を初期値としてニュートン法を３回行い、
kv : : k r a w c z y k a p p r o x ( Func ( ) , x , i x , 3 , 1 ) ; // 候補者集合を作り、解の存在をチェック
}
newton0 : [ 2 ] ( [ 1 , 1 ] , [ − 9 . 9 9 9 9 9 9 9 9 9 9 8 5 3 6 7 9 e −05 , −9.9999999999853678 e −05])
newton1 : [ 2 ] ( [ 1 , 1 ] , [ 2 . 4 2 8 6 1 2 8 6 6 3 6 7 5 2 9 9 e −17 ,2.42861286636753 e −17])
newton2 : [ 2 ] ( [ 1 , 1 ] , [ − 3 . 1 2 2 5 0 2 2 5 6 7 5 8 2 5 2 8 e −17 , −3.1225022567582527 e −17])
I : [ 2 ] ( [ 0 . 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 1 1 , 1 . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 9 ] , [ − 3 . 9 2 0 4 7 5 0 5 5 7 0 7 5 8 4 1 e −16 ,3.29 59746043559 335 e −16])
K : [ 2 ] ( [ 0 . 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 7 7 , 1 . 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 ] , [ − 3 . 1 2 2 5 0 2 2 5 6 7 5 8 4 1 0 6 e −17 ,1.90 81958235745 036 e −16])
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
31 / 72
全解探索の例
#i n c l u d e <kv / a l l s o l . hpp>
namespace ub = b o o s t : : n u m e r i c : : u b l a s ;
s t r u c t Func { // 解きたい問題を関数オブジェクトの形で記述
t e m p l a t e <c l a s s T> ub : : v e c t o r <T> o p e r a t o r ( ) ( c o n s t ub : : v e c t o r <T>& x ) {
ub : : v e c t o r <T> y ( 2 ) ;
y (0) = x (0) ∗ x (1) − cos ( x (1) ) ;
y (1) = x (0) − x (1) + 1;
return y ;
}
};
i n t main ( )
{
ub : : v e c t o r < kv : : i n t e r v a l <d o u b l e> > x ( 2 ) ;
std : : cout . p r e c i s i o n (17) ;
x ( 0 ) = kv : : i n t e r v a l <d o u b l e >(−1000, 1 0 0 0 ) ;
x ( 1 ) = kv : : i n t e r v a l <d o u b l e >(−1000, 1 0 0 0 ) ;
kv : : a l l s o l ( Func ( ) , x ) ; // 全解探索
}
[2]([ −1.964111328125 , −1.47607421875] ,[ −0.66169175448117435 , −0.47607421875]) ( ex )
[2]([ −1.5500093499272621 , −1.5500093499272609] ,[ −0.55000934992726192 , −0.55000934992726113])
( ex : improved )
[ 2 ] ( [ − 0 . 0 1 1 9 6 2 8 9 0 6 2 5 , 0 . 4 7 6 0 7 4 2 1 8 7 5 ] , [ 0 . 9 8 8 0 3 7 1 0 9 3 7 5 , 1 . 4 7 6 0 7 4 2 1 8 7 5 ] ) ( ex )
[ 2 ] ( [ 0 . 2 5 1 1 5 1 8 3 5 2 2 0 7 6 4 5 , 0 . 2 5 1 1 5 1 8 3 5 2 2 0 7 6 5 0 7 ] , [ 1 . 2 5 1 1 5 1 8 3 5 2 2 0 7 6 4 2 , 1 . 2 5 1 1 5 1 8 3 5 2 2 0 7 6 5 4 ] ) ( ex :
improved )
n e t e s t : 4 9 , e x t e s t : 3 , ne : 2 3 , e x : 2 , g i v e u p : 0
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
32 / 72
自動微分 (autodif)
Bottom-Up 型の自動微分を実装している。一階微分のみ。(一変数関
数であれば、ベキ級数型 (psa) で高階微分を行える。)
#i n c l u d e <kv / a u t o d i f . hpp>
namespace ub = b o o s t : : n u m e r i c : : u b l a s ;
// 関数の定義
t e m p l a t e <c l a s s T> ub : : v e c t o r <T> f u n c ( c o n s t
ub : : v e c t o r <T>& x ) {
ub : : v e c t o r <T> y ( 2 ) ;
y (0) = 2. ∗ x (0) ∗ x (0) ∗ x (1) − 1 . ;
y (1) = x (0) + 0.5 ∗ x (1) ∗ x (1) − 2 . ;
return y ;
}
i n t main ( )
{
ub : : v e c t o r <d o u b l e> v1 , v2 ;
ub : : v e c t o r < kv : : a u t o d i f <d o u b l e> > va1 ,
va2 ;
ub : : m a t r i x<d o u b l e> m;
v1 . r e s i z e ( 2 ) ;
v1 ( 0 ) = 5 . ; v1 ( 1 ) = 6 . ;
// 自動微分型の初期化
va1 = kv : : a u t o d i f <d o u b l e >:: i n i t ( v1 ) ;
// 関数呼び出し
va2 = f u n c ( va1 ) ;
// 自動微分型を分解
kv : : a u t o d i f <d o u b l e >:: s p l i t ( va2 , v2 , m) ;
// f ( 5 , 6 )
s t d : : c o u t << v2 << ”\n ” ;
// J a c o b i a n m a t r i x a t ( 5 , 6 )
s t d : : c o u t << m << ”\n ” ;
}
[2](299 ,21)
[2 ,2]((120 ,50) ,(1 ,6) )
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
33 / 72
double-double(dd)
twosum: 2 つの数の和を 2 つの仮数部の重なりのない数の和に変換
するアルゴリズム。(例: 1234 + 5.432 → 1239 + 0.432)
twoproduct: 2 つの数の積を 2 つの仮数部の重なりのない数の和に変
換するアルゴリズム。(例: 1234 × 5.432 → 6703 + 0.088)
twosum と twoproduct を組み合わせると、2 つの倍精度数を用いた擬
似的な 4 倍精度計算を実現できる。一般的な多倍長演算よりかな
り高速。
単体 (dd.hpp) で使った時は近似計算。
dd.hpp と rdd.hpp(方向付き丸めでの dd 型の演算を定義) を併用
し、interval 型の内部型として dd 型を使うと、端点に dd 型を持つ
4 倍精度区間演算が可能。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
34 / 72
mpfr
高精度浮動小数点計算が行える有名な MPFR ライブラリの簡単な
wrapper。
単体 (mpfr.hpp) で使った時は近似計算。
kv::mpfr<106>のようにパラメータとして仮数部長を指定して使う。
mpfr.hpp と rmpfr.hpp を併用し、interval 型の内部型として
mpfr 型を使うと、端点に mpfr 型を持つ区間演算が可能。但し、
MPFR の機能を使うのは加減乗除と平方根のみであり、せっかく
MPFR が持っている優秀な数学関数群は一切用いられない。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
35 / 72
Affine Arithmetic(affine)
Affine Arithmetic とは
区間演算の過大評価を抑制できる。その代わり計算時間がかかる。
全ての変数について、入力変数またはノイズに関する依存性を保持
するため、区間幅の爆発を防げる。
全ての数値は x0 + x1 ε1 + x2 ε2 + · · · + xn εn のような Affine 形式で表
現される。εi は −1 ≤ εi ≤ 1 を動くダミー変数であり、その係数に
より依存性を表現する。
乗除算や数学関数などの非線形演算が出現する度にダミー変数の数
が増え、計算が遅くなる。
内部型は double だけでなく dd や mpfr も入れられる (interval と同様
の仕様)。
ダミー変数の数を削減する機能を持っている。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
36 / 72
Affine Arithmetic の使用例
QRT(Quispel-Roberts-Thompson) 写像
1 + αxn
xn−1 xnσ
σ = 2、α = 2、x0 = x1 = 1 として、区間演算と Affine Arithmeric でどこ
まで計算できるか試してみた。
三項間漸化式 xn+1 =
#i n c l u d e <kv / i n t e r v a l . hpp>
#i n c l u d e <kv / r d o u b l e . hpp>
i n t main ( )
{
int i ;
kv : : i n t e r v a l <d o u b l e> x , y , z ;
std : : cout . p r e c i s i o n (17) ;
x = 1.;
y = 1.;
f o r ( i =2; i <=10000; i ++) {
z = (1 + 2 ∗ y ) / ( x ∗ y ∗ y ) ;
s t d : : c o u t << i << ” ” << z << ”\n ” ;
x = y;
y = z;
}
}
柏木 雅英 (早稲田大学)
#i n c l u d e <kv / a f f i n e . hpp>
i n t main ( )
{
int i ;
kv : : a f f i n e <d o u b l e> x , y , z ;
std : : cout . p r e c i s i o n (17) ;
x = 1.;
y = 1.;
f o r ( i =2; i <=10000; i ++) {
z = (1 + 2 ∗ y ) / ( x ∗ y ∗ y ) ;
s t d : : c o u t << i << ” ” << t o i n t e r v a l (
z ) << ”\n ” ;
x = y;
y = z;
}
}
精度保証付き数値計算の原理と実装
リスク工学セミナー
37 / 72
区間演算と Affine Arithmeric の計算結果
n
区間演算
Affine Arithmetic
2
3
4
5
6
7
8
[3, 3]
[0.77777777777777767, 0.7777777777777778]
[1.408163265306122, 1.4081632653061232]
[2.4744801512287302, 2.4744801512287369]
[0.68995395922102109, 0.68995395922102732]
[2.020393474742363, 2.0203934747424169]
[1.7898074714307314, 1.7898074714308816]
[3, 3]
[0.77777777777777756, 0.77777777777777824]
[1.4081632653061197, 1.4081632653061247]
[2.4744801512287271, 2.4744801512287414]
[0.6899539592210222, 0.68995395922102621]
[2.0203934747423817, 2.0203934747423987]
[1.7898074714307978, 1.7898074714308153]
.
.
.
31
32
33
34
35
36
37
38
39
40
41
42
43
.
.
.
[0.70098916182277204, 0.70941982097935608]
[1.7816188152293368, 1.8444838202503787]
[1.890688867011997, 2.1073484458445711]
[0.58372124794988644, 0.81879304568504608]
[1.5341327531940911, 4.0942614522583929]
[0.29640395761996329, 6.6882779916662063]
[0.0086967943592607538, 106.66548725824453]
[1.3369859317919986 × 10−5 , 9560542.5436595381]
[1.0257053348148149 × 10−16 , 1.229984619387229 × 1019]
[6.9138205018986439 × 10−46 , 1.7488703313159241 × 1056 ]
[2.6581843623384974 × 10−132 , 7.1339291414655989 × 10162 ]
[0, ∞]
−
.
.
.
[0.70519175616865292, 0.70519175616868424]
[1.8127715215496742, 1.8127715215497711]
[1.9960405520559754, 1.9960405520560838]
[0.69119381312156691, 0.69119381312160023]
[2.4982930525184534, 2.4982930525186054]
[1.3900085495715059, 1.390008549571586]
[0.78309678534845506, 0.78309678534849637]
[3.0105168251706007, 3.0105168251708015]
[0.98924416180811902, 0.98924416180817921]
[1.0109923081577889, 1.0109923081578503]
[2.9887738208443805, 2.9887738208445911]
[0.7726251804826496, 0.77262518048269758]
[1.4265918581977079, 1.4265918581978075]
.
.
.
9999
10000
.
.
.
−
−
.
.
.
[0.76071510659932817, 0.76071510667899534]
[1.4727965248961243, 1.4727965251850226]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
38 / 72
区間幅のグラフ
1e+10
interval
dd
affine
100000
1
1e-05
1e-10
1e-15
1e-20
1e-25
1e-30
10
20
30
40
50
60
70
n と区間幅の関係。dd は擬似 4 倍精度区間演算によるもの。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
39 / 72
区間幅のグラフ (n = 10000 まで)
1e+10
interval
dd
affine
100000
1
1e-05
1e-10
1e-15
1e-20
1e-25
1e-30
0
1000
2000
3000
4000
5000
6000
7000
8000
9000
10000
n と区間幅の関係。(n = 10000 まで)
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
40 / 72
ベキ級数演算 (psa)
ベキ級数演算のためのライブラリ。常微分方程式の初期値問題や数値積分の
精度保証に使う。高階微分の計算にも使える。
高次の項を捨ててしまう Type-I PSA と、高次の項を捨てずに区間係数とし
て残す Type-II PSA がある。
Type-II PSA では、最高次の項の係数のみを (幅の広い) 区間とし、それ以外
の低次の項の係数は点区間 (あるいは丸め誤差程度の幅の狭い区間) とする
点に特徴がある。
PSA の例 (積)
1 + 2t − 3t 2 と 1 − t + t 2 の積
Type-I PSA
Type-II PSA
定義域は決めなくてよい
定義域 = [0, 0.1]
1 + t − 4t 2
1 + t + [−4, −3.5]t 2
(1 + 2t − 3t 2 )(1 − t + t 2 ) = 1 + t − 4t 2 + 5t 3 − 3t 4
= 1 + t + (−4 + 5t − 3t 2 )t 2 ∈ 1 + t + [−4, −3.5]t 2
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
41 / 72
Type-I PSA
ベキ級数型
x0 + x1 t + x2 t 2 + · · · + xn t n
ベキ級数型同士の加減乗除や数学関数を考える。
演算結果の n 次までの項を残し、n + 1 次以降は切り捨てる。
1 変数関数の高階微分を計算する自動微分法とほとんど同じ。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
42 / 72
Type-II PSA
ベキ級数型
x0 + x1 t + x2 t 2 + · · · + xn t n
固定された有限閉区間 D = [t1 , t2 ] 上で定義される。
演算結果は n 次までしか保持しないが、n + 1 次以降の項の影響は n
次の項の係数を区間にすることで吸収する。
係数 x0 , · · · , xn は区間。
ただし多くの場合、x0 , · · · , xn−1 は幅の狭い区間、xn は幅の広
い区間。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
43 / 72
Type-II PSA の演算規則 (1)
x(t) = x0 + x1 t + x2 t 2 + · · · + xn t n
y (t) = y0 + y1 t + y2 t 2 + · · · + yn t n
加減算
x(t) ± y (t) = (x0 ± y0 ) + (x1 ± y1 )t + · · · (xn ± yn )t n
加算の例
x(t) = 1 + 2t − 3t 2
y (t) = 1 − t + t 2
x(t) + y (t) = 2 + t − 2t 2
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
44 / 72
Type-II PSA の演算規則 (2)
乗算
(1) まず、打ち切り無しで乗算を行う。
x(t) × y (t)
=
zk
=
z0 + z1 t + · · · + z2n t 2n
∑
min(k,n)
xi yk−i
i=max(0,k−n)
(2) 2n 次から n 次に減次する。
m 次から n 次への減次
x0 + x1 t + x2 t 2 + · · · + xm t m =⇒ z0 + z1 t + · · · + zn t n
zi
=
zn
=
xi (0 ≤ i ≤ n − 1)
{ m
}
∑ i−n
xi t
|t∈D
i=n
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
45 / 72
Type-II PSA の演算規則 (3)
乗算の例
定義域を D = [0, 0.1] とする。
x(t) = 1 + 2t − 3t 2
y (t) = 1 − t + t 2
x(t) × y (t) = 1 + t − 4t 2 + 5t 3 − 3t 4
= 1 + t + (−4 + 5t − 3t 2 )t 2
{
}
∈ 1 + t + −4 + 5t − 3t 2 | t ∈ [0, 0.1] t 2
= 1 + t + [−4, −3.5]t 2
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
46 / 72
Type-II PSA の演算規則 (4)
sin などの数学関数
その関数を g として、
g (x0 + x1 t + · · · + xn t n )
n−1
∑
1 (i)
= g (x0 ) +
g (x0 )(x1 t + · · · + xn t n )i
i!
i=1
({ n
})
∑
1 (n)
+ g
xi t i | t ∈ D
(x1 t + · · · xn t n )n
n!
i=0
のように g の点 x0 での剰余項付きの Taylor 展開に代入することによっ
て得る。この計算中に現れる加算や乗算は Type-II PSA で行う。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
47 / 72
Type-II PSA の演算規則 (5)
除算
x ÷ y = x × (1/y ) と乗算と逆数関数に分解
不定積分
∫
t
x(t)dt = x0 t +
0
柏木 雅英 (早稲田大学)
x1 2
xn n+1
t + ···
t
2
n+1
精度保証付き数値計算の原理と実装
リスク工学セミナー
48 / 72
精度保証付き数値積分
積分区間内に特異点を持たない数値積分の方法を示す。区間 [xi , xi + ∆t]
における積分
∫
xi +∆t
f (t)dt
xi
を次のように計算する。
(1) n 次のベキ級数
(+0t 2 + · · · 0t n )
x(t) = 0 + t
に対して、
∫
y (t) =
t
f (xi + x(t))dt
0
を [0, ∆t] を定義域とした Type-II PSA で計算する。
(2) 計算結果 y (t) を
y (t) = y1 t + y2 t 2 + · · · yn+1 t n+1
とすると、積分値は y (∆t) で得られる。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
49 / 72
ステップ幅 ∆t の決定
ε0 を 1 ステップで混入する誤差の目標値とする。例えば machine epsilon。
(1) Type-I PSA を用いて Taylor 展開を計算し、その係数を見て適切なステッ
プ幅 ∆t0 を推定する。Type-I PSA で計算された Taylor 展開を
x0 + x1 t + x2 t 2 + · · · + xn−1 t n−1 + xn t n
として、
1
∆t0 =
ε0n
1
1
max(|xn−1 | n−1 , |xn | n )
とする。
(2) ステップ幅 ∆t0 を用いて Type-II PSA を使って精度保証付きで 1 ステッ
プの計算を行う。
(3) ステップ幅 ∆t0 で計算して実際に混入した誤差を ε として、新しいス
テップ幅を
(ε )1
0 n
∆t1 = ∆t0
ε
で推定する。ただし、n は Taylor 展開の次数。
(4) ステップ幅 ∆t1 を用いて Type-II PSA を使って精度保証付きで 1 ステッ
プの計算を行う。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
50 / 72
精度保証付き数値積分 (特異点無し) の例 (1)
∫
10
sin(x)
dx
+ 1 + 2−10
cos(x 2 )
0
使用プログラム
#i n c l u d e <i o s t r e a m >
#i n c l u d e <kv / d e f i n t . hpp>
t y p e d e f kv : : i n t e r v a l <d o u b l e> i t v ;
s t r u c t Func {
t e m p l a t e <c l a s s T> T o p e r a t o r ( ) ( c o n s t T& x ) {
r e t u r n s i n ( x ) / ( c o s ( x∗x ) + 1 . + pow ( 2 . , −10) ) ;
}
};
i n t main ( ) {
std : : cout . p r e c i s i o n (17) ;
s t d : : c o u t << kv : : d e f i n t a u t o s t e p ( Func ( ) , ( i t v ) 0 . , ( i t v ) 1 0 . , 1 0 ) << ”\n ” ;
}
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
51 / 72
精度保証付き数値積分 (特異点無し) の例 (2)
∫
0
10
kv-0.4.23
octave 3.8.1
intlab 9
Mathematica 10.1.0
matlab 2007b
keisan (ロンバーグ)
keisan (Tanh-Sinh)
keisan (ガウス-ルジャンドル)
intde2 by ooura
python + scipy
CASIO fx-5800P
柏木 雅英 (早稲田大学)
sin(x)
dx
+ 1 + 2−10
cos(x 2 )
[38.383526264535703,38.383526264649654]
38.3837105761501
[38.34845927756175, 38.41859325162576]
0.0608979
38.383519835854528
38.324147930794
38.24858948837754677984
116.448156707725851273
32.4641
36.48985372847387
38.38352669
精度保証付き数値計算の原理と実装
リスク工学セミナー
52 / 72
精度保証付き数値積分 (特異点無し) の例 (3)
1
integrand
step size
1000
0.1
500
step size
integrand
0.01
0
0.001
-500
0.0001
-1000
1e-05
0
柏木 雅英 (早稲田大学)
2
4
6
精度保証付き数値計算の原理と実装
8
10
リスク工学セミナー
53 / 72
常微分方程式の初期値問題
一階連立常微分方程式
dx
= f (x, t),
dt
x(t0 ) = x0
x ∈ Rl , t ∈ R
初期値問題の精度保証アルゴリズムは、t0 < t1 < t2 < . . . に対して、
x(ti ) を元に x(ti+1 ) を精度保証付きで計算する方法 (短い区間での
精度保証)
短い区間での精度保証法を利用して、長い区間に渡って (区間幅の膨
らみを抑制しながら) 接続する方法
の 2 つに分けて考えることが出来る。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
54 / 72
短い区間での精度保証 (1)
平行移動 & 両辺を積分で不動点形式に
∫
x(t) = v +
t
f (x(t), t + ti )dt
0
(v = x(ti ),
t ∈ [0, ti+1 − ti ])
解の Taylor 展開の生成
Type-I PSA 型の変数 X0 = v , T = t を用いて、k = 0 とし、
(1) 次数 k の Type-I PSA で以下を計算
∫ t
Xk+1 = v +
f (Xk , T + ti )dt
0
(2) 次数 k = k + 1 とする。
を n 回繰り返すと、Xn として解の n 次の Taylor 展開が得られる。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
55 / 72
短い区間での精度保証 (2)
解の存在保証
Type-II PSA の定義域を D = [0, ti+1 − ti ] と設定し、Type-I PSA の反復で得られ
た n 次の Taylor 近似
Xn = x0 + x1 t + x2 t 2 + · · · + xn t n
と T = t を用いて、
(1) Xn の最終項の係数を膨らませた候補者集合
Y = x0 + x1 t + x2 t 2 + · · · + V t n
を作成する。
∫t
(2) v + 0 f (Y , T + ts )dt を次数 n の Type-II PSA で計算し、n + 1 次から n 次に
減次したものを
Y1 = x0 + x1 t + x2 t 2 + · · · + V1 t n
とする。n − 1 次までの係数は Xn と全く同じになることに注意。
(3) V1 ⊂ V なら Y1 内に解の存在が保証される。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
56 / 72
短い区間での精度保証 (3)
候補者集合の作成
候補者集合の作成は、例えば次の手順で行う。
∫t
(1) v + 0 f (Xn , T + ts )dt を次数 n の Type-II PSA で計算し、n + 1 次か
ら n 次に減次したものを Y0 = x0 + x1 t + · · · + V0 t n とする。
(2) r = ||V0 − xn || とし、
V = xn + 2r ([−1, 1], . . . , [−1, 1])T
とする。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
57 / 72
短い区間での精度保証例
dx
= −x 2
dt
x(0) = 1, t ∈ [0, 0.1]
展開の次数: n = 2、10 進 3 桁演算。
(Type-I PSA による Taylor 展開の生成)
X0
X1
=
1
∫
∫
t
Y0 = 1 − t + [0.9, 1]t 2
t
(−X02 )dt = 1 +
=
1+
=
1−t
∫ t
∫ t
Y0 = 1 − t + [0.8, 1.2]t 2
(−(1 − t)2 )dt
1+
(−X12 )dt = 1 +
0
0
∫ t
(Type-II PSA による精度保証)
1+
(−(1 − 2t))dt
0
X2
2 次に減次して、
=
=
(−1)dt
0
r = ||[0.9, 1] − 1|| = 0.1 なので、
0
=
∫
1 − t + t2
=
(−X22 )dt
∫
2 次に減次して、
0
t
(−(1 − t + t 2 )2 )dt
1+
∫
1 − t + t 2 + [−1.133, −0.786]t 3
t
1+
=
(−Y02 )dt
0
(候補者集合の生成)
∫
t
1+
Y1 = 1 − t + [0.886, 1]t 2
0
t
(−(1 − 2t + [2.8, 3]t ))dt
2
=
1+
=
1 − t + t 2 + [−1, −0.933]t 3
0
柏木 雅英 (早稲田大学)
[0.886, 1] ⊂ [0.8, 1.2] なので、Y1 内に真の解
が存在する。
精度保証付き数値計算の原理と実装
リスク工学セミナー
58 / 72
長い区間に渡って接続する (1)
推進写像
t = ts における値 v = x(ts ) に対して、x(te ) を対応させる写像
ϕts ,te : Rs → Rs ,
ϕts ,te : x(ts ) 7→ x(te )
を推進写像と呼ぶことにする。
初期値に関する変分方程式
x ∗ (t) を v を初期値とした与式の解とすると、
d
y (t) = fx (x ∗ (t), t)y (t),
dt
y (ts ) = I , t ∈ [ts , te ]
y ∈ Rs×s
を解くことによって、推進写像の微分を ϕ′ts ,te (v ) = y (te ) で計算出来る。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
59 / 72
長い区間に渡って接続する (2)
時刻 t0 < t1 < t2 < · · · における解の包含を Xi とする。
平均値形式
Xi+1 = ϕti ,ti+1 (mid(Xi )) + ϕ′ti ,ti+1 (Xi )(Xi − mid(Xi ))
単純にこのまま計算すると、wrapping effect によって区間幅が激し
く増大する。
affine arithmetic を使って接続する
ステップ幅の調整は、数値積分の場合と同じアルゴリズムを使う。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
60 / 72
例題: van der Pol 方程式
van del Pol 方程式
d 2x
dx
− µ(1 − x 2 )
+x =0
2
dt
dt
一階に直す
dx
dt
dy
dt
柏木 雅英 (早稲田大学)
= y
= µ(1 − x 2 )y − x
精度保証付き数値計算の原理と実装
リスク工学セミナー
61 / 72
常微分方程式の初期値問題の例
#i n c l u d e <kv / ode−m a f f i n e . hpp>
namespace ub = b o o s t : : n u m e r i c : : u b l a s ;
t y p e d e f kv : : i n t e r v a l <d o u b l e> i t v ;
c l a s s VDP { // 解きたい問題の右辺を関数オブジェクトで記述
public :
t e m p l a t e <c l a s s T> ub : : v e c t o r <T> o p e r a t o r ( ) ( c o n s t ub : : v e c t o r <T>& x , T t ) {
ub : : v e c t o r <T> y ( 2 ) ;
y (0) = x (1) ;
y ( 1 ) = 1 . ∗ ( 1 . − x ( 0 ) ∗x ( 0 ) ) ∗ x ( 1 ) − x ( 0 ) ;
return y ;
}
};
i n t main ( )
{
ub : : v e c t o r <i t v > x ;
i t v end ;
std : : cout . p r e c i s i o n (17) ;
x . r e s i z e (2) ;
x ( 0 ) = 1 . ; // 初期値
x (1) = 1;
end = 1 0 0 . ; // 終了時刻
o d e l o n g m a f f i n e (VDP( ) , x ,
s t d : : c o u t << x << ”\n ” ;
}
i t v ( 0 . ) , end ) ; // 初期値問題を解く(0−end )
[2]([2.0077904809521123 ,2.0077904809521399] ,[ −0.056051438751190147 , −0.056051438750530216])
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
62 / 72
境界値問題の例 (1)
3
2
1
0
-1
-2
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
2.5
van der Pol 方程式 (µ = 1) の周期解を Poincaré Map に対する Krawczyk
法で精度保証した例。周期は、
T = [6.6632868593231044, 6.6632868593231534]
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
63 / 72
境界値問題の例 (2)
#i n c l u d e <kv / p o i n c a r e m a p . hpp>
#i n c l u d e <kv / kraw−a p p r o x . hpp>
namespace ub = b o o s t : : n u m e r i c : : u b l a s ;
t y p e d e f kv : : i n t e r v a l <d o u b l e> i t v d ;
c l a s s VDP {
public :
t e m p l a t e <c l a s s T> ub : : v e c t o r <T> o p e r a t o r ( ) ( ub : : v e c t o r <T> x , T t ) {
ub : : v e c t o r <T> y ( 2 ) ;
y (0) = x (1) ;
y ( 1 ) = 1 . ∗ ( 1 . − x ( 0 ) ∗x ( 0 ) ) ∗ x ( 1 ) − x ( 0 ) ;
return y ;
}
};
c l a s s VDPPoincareSection {
public :
t e m p l a t e <c l a s s T> T o p e r a t o r ( ) ( ub : : v e c t o r <T> x ) {
T y;
y = x (0) − 0 . ;
return y ;
}
};
i n t main ( )
{
ub : : v e c t o r <d o u b l e> x ;
ub : : v e c t o r <i t v d > i x ;
std : : cout . p r e c i s i o n (17) ;
VDP f ;
VDPPoincareSection g ;
kv : : PoincareMap<VDP, V D P P o i n c a r e S e c t i o n , d o u b l e> h ( f , g , ( i t v d ) 0 . ) ;
x . r e s i z e (3) ; x (0) = 0 . ; x (1) = 1 . ; x (2) = 6.28;
kv : : k r a w c z y k a p p r o x ( h , x , i x , 1 0 , 0 ) ;
s t d : : c o u t << i x << s t d : : e n d l ;
}
[ 3 ] ( [ − 5 . 4 5 8 7 3 4 5 6 8 7 1 0 3 1 5 7 e −30 ,5.458734568710315 e
−30] ,[2.1727136926224956 ,2.1727136926225979] ,[6.6632868593231044 ,6.6632868593231534])
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
64 / 72
ローレンツ方程式
dx
dt
dy
dt
dz
dt
= 10(y − x)
= 28x − y − xz
z
8
= − z + xy
3
45
40
35
30
25
20
15
10
5
-20
-15
(x(0), y (0), z(0)) = (36, 15, 15)
-10
-5
0
x
5
10
15
-10
-15
-20
20-25
-5
0
5
10
15
20
25
y
カオス現象の発見で有名な方程式
matlab という世界で最も信頼されている数値計算ソフトウェアと、
自作の精度保証付き数値計算プログラムで解き比べてみる。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
65 / 72
計算結果をグラフで見ると · · ·
matlab 側は ode45 に最高精度 (相対誤差 10−15 ) を指定、精度保証側
は kv ライブラリ。
横軸が時刻 t 、縦軸は x, y , z の値。
最初はぴったり重なっているように見える。
t = 13 あたりからわずかにずれが見えはじめ、t = 15 あたりで決定
的に運命が分かれ、それ以降は異なる人生を歩んでいるかのよう。
t = 28 付近になると、精度保証付き数値計算プログラムの方のグラ
フに幅が出てくる。実はずっと上限と下限の 2 本のグラフを描画し
ていた。このへんまでが double での限界。
double-double を使えば t = 50 あたりまで行ける。
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
66 / 72
ローレンツ方程式 (t = 6 まで)
50
40
30
20
10
0
-10
-20
-30
0
1
2
verified(x)
verified(y)
柏木 雅英 (早稲田大学)
3
verified(z)
matlab(x)
4
5
6
matlab(y)
matlab(z)
精度保証付き数値計算の原理と実装
リスク工学セミナー
67 / 72
ローレンツ方程式 (t = 12 まで)
50
40
30
20
10
0
-10
-20
-30
6
7
8
verified(x)
verified(y)
柏木 雅英 (早稲田大学)
9
verified(z)
matlab(x)
10
11
12
matlab(y)
matlab(z)
精度保証付き数値計算の原理と実装
リスク工学セミナー
68 / 72
ローレンツ方程式 (t = 18 まで)
50
40
30
20
10
0
-10
-20
-30
12
13
14
verified(x)
verified(y)
柏木 雅英 (早稲田大学)
15
verified(z)
matlab(x)
16
17
18
matlab(y)
matlab(z)
精度保証付き数値計算の原理と実装
リスク工学セミナー
69 / 72
ローレンツ方程式 (t = 24 まで)
50
40
30
20
10
0
-10
-20
-30
18
19
20
verified(x)
verified(y)
柏木 雅英 (早稲田大学)
21
verified(z)
matlab(x)
22
23
24
matlab(y)
matlab(z)
精度保証付き数値計算の原理と実装
リスク工学セミナー
70 / 72
ローレンツ方程式 (t = 28.5 まで)
50
40
30
20
10
0
-10
-20
-30
24
25
verified(x)
verified(y)
柏木 雅英 (早稲田大学)
26
verified(z)
matlab(x)
27
28
matlab(y)
matlab(z)
精度保証付き数値計算の原理と実装
リスク工学セミナー
71 / 72
まとめ (KV ライブラリ)
http://verifiedby.me/kv/ で公開中。
作成開始は 2007 年秋頃。公開開始は 2013 年 9 月 18 日。
言語は C++。boost C++ Libraries も必要。
全てヘッダファイルで記述されており、インストールはヘッダファ
イルをどこかに配置するだけ。オープンソース。
計算に使う数値の型を double 以外の型に容易に変更することが出
来る。
(数値型) 区間演算 (多数の数学関数含む)、4 倍精度 (double-double)
演算、MPFR ラッパー、複素数演算、自動微分、affine arithmeric、
ベキ級数演算
(アプリケーション) Krawczyk 法による非線形方程式の精度保証、非
線形方程式の全解探索、常微分方程式の初期値問題、常微分方程式
の境界値問題、数値積分、特殊関数
皆様のご利用をお待ちしております!
柏木 雅英 (早稲田大学)
精度保証付き数値計算の原理と実装
リスク工学セミナー
72 / 72