筆者は「精度保証付き数値計算」という分野で研究をしている大学院生です．
「数値計算は分かるけど」「精度保証付き数値計算？ナニソレ？」という方がほとんどだと思います．
「精度保証付き数値計算」の研究自体は30年ほど前から盛んに行われていますが，世間に浸透しているとは言えない状況です．
自分の研究分野が世間に知られていないのは何か少し寂しい感じがするので「精度保証付き数値計算」を少しでも広めるべく記事を投稿することにしました．(シリーズ化するかも知れません)
本日は「精度保証付き数値計算」というワードだけでも覚えていただければ幸いです．
今回は"数値計算に潜むとんでもないリスク"に関してカジュアルにお話します．
そして筆者の研究分野である「精度保証付き数値計算」の必要性を知ってもらえればなと思います．
この記事を読み終える頃には計算機を信頼できなくなっているかも知れません(笑)
※不安を煽ることが目的の記事ではないです．
0. 注意
本記事では説明をする際にPython3を使用します．普段はMATLABを愛用していますが「読めねーよ」と言われそうなので，自分の使える言語の中でも可読性の高いものを・・・と思いPython3を使用することにしました．
また，なるべく専門用語は使わずに書いたつもりなのですが，もし「ここがよく分からない」というところがあれば教えて欲しいです．随時加筆修正を加えます．
1. さて，突然ですが・・・
数値計算には誤差が含まれるのをご存知ですか？
「そんなの知っとるわ！」という方もブラウザバックせず少々お付き合いください．
有名なもので
print(0.1+0.1+0.1)
0.30000000000000004
がありますね．これは0.1を2進表現にした際に有限桁で表せないことに起因します．
これはあくまで一例で，数値計算には他にも様々な誤差があります．(「丸め誤差」「打ち切り誤差」「離散化誤差」などでググってみてください．)
数値計算に誤差が含まれることについてはご存知の方が大半だとは思いますが，意外と軽視されていると感じるのでその"とんでもないリスク"に関して紹介します．
2. まずはこれを見てくれ
2次方程式 $x^2+10^{15}x+10^{14}=0$ の解のうち大きい方を求めよ．
こんな問題を出されたらあなたはどうしますか？
当然，手計算するのは大変ですから数値計算に頼ることでしょう．
では，実装の方法についてはどうでしょうか？
実は実装の仕方によっては解が大きく異なることがあります．
解の公式を用いて，
$x=\frac{-b+\sqrt{b^2-4ac}}{2a}$ を計算した場合と，分子を有理化した
$x=\frac{2c}{-b-\sqrt{b^2-4ac}}$ を計算した場合を見てみましょう．
数学的には同じ解が出るはずですが・・・
a=1
b=10**15
c=10**14
print((-b+(b**2-4*a*c)**.5)/(2*a))
print((2*c)/(-b-(b**2-4*a*c)**.5))
-0.125
-0.10000000000000002
※筆者は驚くと猫になってしまいます．
2次方程式の解の計算でさえこのように誤差が生じ，結局どっちが真の解かも分かりません．
普段計算機でしている計算の結果が本当に正しいものなのか，ちょっと心配になってきますね．
実は2次方程式に関しては誤差混入の原理が明らかになっていて誤差を回避する方法も知られています．
数値計算は前の計算結果を使って次の計算を行なうので，どこか一箇所でも誤差が混入してしまうとそれ以降の計算は信頼できないものになります．
大量の計算を行なう過程のどこかで致命的な誤差の混入が発生していたら…と考えると恐ろしくて夜しか眠れません．
さて，ここからが本題です．(唐突)
本当に紹介したかったのは次の例で，筆者がこれを知った時は驚きのあまり研究室で叫びました．
計算機不信に陥るので，覚悟してください(笑)
3. 数値計算に潜むとんでもないリスク
では，いよいよ本題の"数値計算に潜むとんでもないリスク"を
「Rumpの例題[Rump 1988]」を例に紹介します．
[Rumpの例題]
$f(a,b)=333.75b^6+a^2(11a^2b^2-b^6-121b^4-2)+5.5b^8+\frac{a}{2b}$ に $(a,b)=(77617.0,33096.0)$ を代入した値は？
IBMのメインフレームS/370を使用して，単精度，倍精度，拡張精度で数値計算を行うと，それぞれ
- 単精度： $f(a,b)\approx1.172603\dots$
- 倍精度： $f(a,b)\approx1.1726039400531\dots$
- 拡張精度： $f(a,b)\approx1.172603940053178\dots$
となるそうな．
「ふむふむなるほど．まぁじゃあ解は約$1.17$ってことね……」
実はこれは大きな間違いで，実際に頑張って手計算すると真の解は $f(a,b)=-0.82739605\dots$ になります．
符号すら違う………？
1つ1つは小さな誤差でも塵も積もれば山となる…で大変なことになりました．
このRumpの例題を通して我々が学ばなくてはならないことは
- 数値計算にはとんでもないリスクが潜んでいること
- 複数の計算機で似通った解が出力されたからといって，それが真の解とは限らないこと
の2点です．
飛行機を飛ばすのにも数値計算は使われているわけですが，人命がかかっている数値計算で「あ，計算間違えちゃった(笑)」では済まされませんから，これからは「計算の速度」だけでなく「計算の品質」が重要になってくると筆者は考えます．
4. 精度保証付き数値計算の必要性
世間で信頼されているソフトウェアでさえ，上記の例のように何の警告も出さずに間違った解を出力することがあります．
「精度保証付き数値計算」では数学の理論を用いて数値計算に含まれる誤差を厳密に把握し，正確な計算が不可能な場合はなんらかの形で警告を発してくれます．黙って嘘をつくことはありません．
夢は大きく全ての数値計算を精度保証付き数値計算に置き換えることをモチベーションに日々研究をしています．目指せ！ストレスフリーな社会ならぬエラーフリーな社会！
5. 終わりに
冒頭でもお話しましたが「精度保証付き数値計算」の研究はあまり広くは知られていないというのが現状です．
筆者の指導教員は「精度保証の方法を追求するのではなく，精度保証を使って何ができるのかを示していかないといけない」とおっしゃっていましたが，全く同感です．
この分野の研究をするには
- 計算機の知識が必要
- 高度な数学の知識が必要
- 理論だけでなくプログラミングによる実装力が必要
とハードルが高く，(筆者も最初は苦労しました)
それだけで敬遠されがちですが，なるべくそういった難しい話はやめて精度保証付き数値計算を使って何を実現したのかを世に発信していき，みなさんにはその成果を使っていただくだけという状態にしていくことが目標です．
Qiitaでは，上記のような難しい話をしつつも，何か役に立つものが作れた時は紹介していこうかなと考えています．当分は難しい話になりそうです(笑)
今考えているもの
最後に
この記事で「ふ〜ん，精度保証付き数値計算………必要じゃん」と思ってくださった方や精度保証に少しでも興味を持ってくださった方がいれば幸いです．
ここまで読んでいただきありがとうございました！
6. 参考
精度保証付き数値計算(Wikipedia)
精度保証付き数値計算の必要性(柏木先生の投稿)
TakLAB日誌(筆者の所属する研究室のブログ，筆者も記事を書いています)