この文章は筑波大学大学院システム情報工学研究科リスク工学専攻の紀要「リスク工学研究」Vol. 16, pp.27-30に掲載された内容をJuliaを用いた計算例を加えて加筆修正したものです。
「数値計算・数値解析」という言葉、日頃よく耳にすると思いますが、「数値計算・数値解析」とは何でしょうか。本稿は、まずこの二つの言葉を明確に区別する事から始めます。今日、筆者が知る限り、「数値解析」は同じ言葉で二つの使い方をされます。一つは、「数値解析」とは数値を用いて代数的操作によって解くことができない数学の問題を解決する手法とする使い方、数値計算・数値実験・数値シミュレーションなどと同義のように使用されます。もう一つは、「数値解析」とは応用数学の一分野で、上記の数学問題を、数値を用いて近似的に解く手法に関する数学的な概念を研究する分野とする使い方。筆者は後者の立場で、数値解析という言葉を使います。すなわち、数値を用いた計算を数値計算、数値計算に関する数学分野を数値解析とします。
さて、区別をはっきりさせたところで質問があります。数値計算の結果はいつも正しいのでしょうか？普段、数値計算を頻繁に利用するけれども、その計算結果を闇雲に信じてはいませんか。あるいは「数値計算には誤差がある」という知識がある人は計算結果は正しいかもしれないし、間違っているときもあるだろうと認識されているかもしれません。本稿の目的は、数値計算の結果が間違える事がある例をいくつか紹介し、数値計算の不正確さ・不確実性を明確にし、数値計算による間違いのリスクを明らかにする事です。
では、古典的な電卓（電子卓上計算機）の計算結果から紹介します。お手元の電卓 に以下の計算をさせてみて下さい。
$$ \begin{array}{rrrl} &100 & 000 & 000\\ -)& 9 & 999 & 999.99 \\ \hline \end{array} $$答えは
$$ \begin{array}{rrrl} &100 & 000 & 000\\ -)& 9 & 999 & 999.99 \\ \hline &\color{red}{90} &\color{red} {000} &\color{red} {000.0}~(!?) \end{array} $$流石にこれはまずい。一目で間違えに気づきます。筆者の手元のiPhoneの電卓は9桁の表示領域があり、これに2桁加えた11桁の固定小数で計算が実行されます。すると
$$ \begin{array}{rrrl} &100 & 000 & 000.\color{red} {00}|\\ -)& 9 & 999 & 999.99|\color{red} {00}\\ \hline &90 & 000 & 000.01 \end{array} $$一見、正しい答えが計算されているように見えますが、表示できる数値が9桁のため
$$ \begin{array}{rrrl} &100 & 000 & 000.00\\ -)& 9 & 999 & 999.99{00}\\ \hline &90 & 000 & 000.0\not{1} \end{array} $$このように最後の「１」は消えてしまいます。数値計算が間違えた時、一番の問題は結果が間違っていても実行した計算機は知らせてくれない事です。この問題は結果を見れば間違いに気づきますが、少し複雑化した次の場合はどうでしょうか。
$3$ 以上の自然数 $n$ について、
$$ x^n + y^n = z^n $$となる $0$ でない自然数 $(x, y, z)$ の組が存在しない。
これは17世紀、数学者のフェルマーによって古代ギリシャの数学者ディオファントスの著作「算術」の余白に「この定理に関して、私は真に驚くべき証明を見つけたが、この余白はそれを書くには狭すぎる」というメモとともに書かれた定理で、長い間数学の未解決問題として残っていた問題です。そして1995年、この定理はイギリスの数学者ワイルズによって完全に証明され （文献1）、その証明の顛末を含めて、大変脚光を浴びました。
筆者はこの定理の証明に関して、これ以上書ける事はないのですが、数値計算を使うと以下の主張ができます。
$(x,y,z)=(139,954,2115)$ において
$$ 139^3 + 954^3 = 2115^3 $$が成り立つ！？
これはフェルマーの最終定理の反例です。根拠として以下を示します。
x = Int32(139)
y = Int32(954)
z = Int32(2115)
if x ^ 3 + y ^ 3 == z ^ 3
println("Counter example of Fermat's theorem")
println("(x, y, z) = ($x, $y, $z)")
end
もちろん（！）筆者の数値計算結果は間違っており、これは符号付き整数の計算機内での表現が引き起こす間違いが原因です。正しくは
$$ (x^3,y^3,z^3)=(2685619,868250664,9460870875) $$となり反例にはなりません。数の表現範囲（精度）が足りなかったのです。
では十分精度があれば、正しい結果が得られるのでしょうか？先の例では32ビット符号付き整数という数値を用いて失敗しました。もっと数の表現範囲が増える64ビット浮動小数点数（binary64と呼ばれる計算機上での小数の標準的表現方法）や128ビットのdouble-double (dd)型（四倍精度には仮数部が足りない擬似四倍精度）、あるいはそれ以上の桁数の精度を使って計算したら…
次の例は Rumpの例題（文献2）として有名な例題です。
$2$ 変数 $a$, $b$ を引数に持つ非線形関数
$$ f(a,b)=(333.75-a^2)b^6+a^2(11a^2b^2-121b^4-2)+5.5b^8+\frac{a}{2b} $$について $a=77617$, $b=33096$ のときの値を求めよ。
この関数、最大8次の多項式があるので、手計算はまずやりたくありません。そこで数値計算の出番ですが、先ほど懲りましたので、今回はいくつかの精度（Float32、Float64、dd型、真の四倍精度に相当する仮数部113桁のBigFloat型）で計算をしてみます。
注意：Juliaでdd型を実装するパッケージとしてここではDoubleFloatsを使いました。setprecision(BigFloat, 106)としても同じ計算ができます。
using DoubleFloats
# f(a, b) = (333.75 - a ^ 2) * b^6 + a ^ 2 * (11a ^ 2 * b ^ 2 - 121b ^ 4 - 2) + 5.5b ^ 8 + a / (2b)
function f(a::T, b::T) where T
return (convert(T,333.75) - a^2) * b^6 + a^2 * (convert(T,11) * a^2 * b^2 - convert(T,121) * b^4 - convert(T,2)) + convert(T,5.5) * b^8 + a / (convert(T,2)*b)
end
a = Float32(77617)
b = Float32(33096)
c = f(a,b)
println("float: $c")
a = Float64(77617)
b = Float64(33096)
c = f(a,b)
println("double: $c")
a = Double64(77617)
b = Double64(33096)
print("dd: ")
showall(f(a,b))
println("")
setprecision(113)
a = BigFloat(77617)
b = BigFloat(33096)
c = f(a,b)
println("qd: $c")
計算結果はどの数値も $1.17260\dots$ であり、おそらくこの数値が正しい値であろうと数値計算の結果から予想できます。しかし、真の値は $f(a,b)=-0.827396059946821\dots$ 符号さえも合っていません。一体何が起こっているのでしょうか。
Rumpの例題は $a=77617$, $b=33096$ のとき $a^2=5.5b^2+1$ という等式が成立する事を巧みに利用して、「桁落ち」という数値計算の誤差が発生するように作られています。従って、単に精度を上げれば計算が正しいとは限らない。さらに恐ろしい事に、数値計算をしている際に突如このような間違いが起らないと言い切れないという事です。
ちなみにRumpの例題は仮数部を122桁に設定した小数を使うと正しく計算できます。binary128の113桁ではちょっと足りないというRump先生の匠の技であることがうかがえます。そして普通、そこまで多数桁の計算をする気にはならない。
setprecision(122)
a = BigFloat(77617)
b = BigFloat(33096)
c = f(a,b)
println("BigFloat: $c")
またRumpの例題の亜種として柏木の例題「$a=10^9$ に対して、$(a^2+1)(a+1)(a-1)-a^4+1/7$ のときの値を求めよ」というもの（文献3）もあります。
function g(a::T) where T
return (a^2 + convert(T,1))*(a + convert(T,1))*(a - convert(T,1))-a^4 + convert(T,1)/7.
end
a = Float32(10^9)
@show g(a)
a = Float64(10^9)
@show g(a)
a = Double64(10^9)
@show g(a)
setprecision(120)
a = BigFloat(10^9)
@show g(a);
ここまでは関数値の評価などの数値演算に起こる間違いの事例を紹介しました。もう少し実用的な数値計算についても考えます。基本的な数値計算として、連立一次方程式の求解を考えましょう。数値計算法としては掃き出し法やGaussの消去法という名前の計算方法を使います。いま、連立一次方程式
$$ \left\{\begin{array}{c} 64919121 x-159018721 y=1 \\ 41869520.5 x-102558961 y=0 \end{array}\right. $$が与えられた（文献4）とします。連立一次方程式の求解も手計算は複雑になるため、できれば数値計算を実行したい。実行速度はとても速く、解は $x=1.45867\times 10^8$, $y=5.95501\times 10^7$ と瞬時に計算されます。しかし、実際の解は
$$ x=205117922,\quad y=83739041. $$またしても数値計算は計算間違いを指摘することができずに、間違った結果を返してしまいました。連立一次方程式の求解は様々な数値解法に現れる最も基本的な数値計算です。さらに大規模な問題に対しては数値計算の結果を確認しづらい面もあります。そのような問題を数値計算する際、解が正しいのか疑い出すと怖くて数値計算なんてできなくなってしまいます。
A = [64919121 -159018721; 41869520.5 -102558961]
b = [1; 0]
x = A\b
一体どのようにしたら数値計算の間違いを正しく認識する事ができるのでしょうか。そのための一つの方法として、数学の知識を使う事が考えられます。電卓の例は自明にしても、フェルマーの最終定理に対してはこの定理が数学的に証明されている事から反例は原則ありえません。Rumpの例題に対しては、「桁落ち」という数値計算に起こる誤差を知っておくと怪しいと考える事ができます。さらに連立一次方程式では「行列の条件数」という数を考えると
$$ \mathrm{cond}(A) = 1.52008\times 10^{16} $$という数値が計算され、これは数値解の相対的な誤差が最大10進16桁ほど混入する可能性を示唆しています。先に述べた64ビット浮動小数点数が、概ね16桁の精度を持つと言われているため、この誤差の大きさは致命的です。すなわち解は1桁も合わない事になります。実際に前節の例では1桁も数値が合いませんでした。
using LinearAlgebra
cond(A)
このように、数学の知識と数値計算を合わせる事で、数値計算に対するリスクを制御しようというのが一つのアプローチ方法です。そしてこの考え方が筆者の研究活動の根源です。 数値計算の名誉挽回のために、これらは人工的に作成された問題である事を強調しておきます。しかし、このような問題が実際の数値計算において起きうる事は数値計算のリスクとして認識して欲しいところです。今日の計算機で使用されている64ビット浮動小数点数の規格の制定に尽力したW. M. Kahanの言葉を借りると、
浮動小数点演算によって得られた結果と真値に大きな差が生じることは非常に稀であり、つねに心配するにはあまりにも稀であるが、だからといって無視できるほど稀なわけではない。
という事になります。非常にオブラートに包まれた感じがしますが、毎日の数値計算に対して杞憂するほど、数値計算は間違えない事も事実です。
最後に、筆者の研究テーマである「精度保証付き数値計算」について紹介して、本稿を擱筆したいと思います。「精度保証付き数値計算」とは数学的に正しい結果を数値計算によって導く手法全般、数値計算結果の品質保証だけでなく、計算機を援用する数学解析手法も「精度保証付き数値計算」といいます。精度保証付き数値計算を実行するためには数値計算に生じる全ての誤差（図1、離散化誤差、打ち切り誤差、丸め誤差等）をすべて把握する必要があります。手法の基本原則として、区間の端点に浮動小数点数を用いる（機械）区間演算という演算規則を使用して、数値の代わりに区間を利用して数値計算を実行します。精度保証付き数値計算は計算を間違えないので、先のような例に遭遇すると区間の幅が著しく増大したり、解の検証が失敗して、エラーとなり警告を発します。これにより数値計算の間違いを見落とす事がありません。
一方で、課題として、本来の数値計算に比べて実行速度が遅い（検算を伴うため、原則、数値計算よりも実行時間を要する）、警告を発して計算が度々失敗するなど、本来の数値計算に比べて「有用性」の面で劣る部分があります。従って、精度保証付き数値計算の有用性を世の中に浸透させ、多くの人に使ってもらうための研究が、今後、本分野の研究者には必須です。また、「精度保証付き数値計算は難しい」、「数学も計算機の知識も必要で敷居が高い」、「細かな職人芸が必要になる」、「プログラミングによる実装力が必要」という理由で研究対象としても敬遠されがちです。
筆者の研究に対する基本姿勢は、精度保証付き数値計算を使って何ができるかを示す事であり、様々な数理モデルに現れる微分方程式（偏微分方程式・遅延微分方程式・常微分方程式）を対象に、精度保証付き数値計算による計算機援用証明を研究主題としています。数理モデルの間違えない数値計算は、すなわち現象の正しい理解へとつながり、モデルの信頼性検証が可能になります。今後も精度保証付き数値計算を使ってできる事を発信し続けます。
本資料はルンプの例題のS.M. Rump先生、柏木雅英先生との雑談が元になっています。他にも精度保証付き数値計算の研究に従事されている先生方には日頃から様々なことをご教示いただいています。ここに感謝の意を表します。また、以下のような文献・Web ページ等を参考にこの文章は書いています。