計算機イプシロンのこと

浮動小数点数について等しいかどうかを調べる場合、普通に等値比較をおこなうのではなく、計算機イプシロンepsを使って「|x - y| / |x| < eps」とか「|x - y| / max(|x|, |y|) < eps」としなさい、という話を聞いたり読んだりしたことがある。でも何かおかしい気がする。それ以前に計算機イプシロンが何なのかもよくわからない。なので、そのことについてメモにまとめる。

浮動小数点数の等値判定の仕方

文章が長いので、浮動小数点数の等値判定の仕方を先に書いておく。

適当な定数Cを取って
\frac{|x-y|}{\max(|x|,|y|)} < C
で等値判定をする。または同じことだけど、適当な定数cと計算機イプシロンεを使って
\frac{|x-y|}{\max(|x|,|y|)} < c \epsilon
で等値判定をする。定数Cやcの決め方は状況次第。

倍精度計算でC=10-5〜10-10ぐらいに取っても、桁落ち等でこれより誤差が大きくなってしまうことも珍しくない(倍精度での計算機イプシロンは10-16ぐらい)。そのような誤差を許容するのかそれとも誤差が減るように計算方法を変えるのかに関して、決まった答えがあるわけではない。

具体例

たとえば0.12 - \frac{0.77}{7.0}を正確に計算すると答えは0.01。
でもこの式を浮動小数点数演算で計算したときの0.01との相対誤差は、計算機イプシロンε=2-52の2倍ちょっとになっている。つまりこのような簡単な計算の場合でも、引き算によるちょっとした桁落ちで計算精度が落ちる。だから相対誤差の許容値は計算機イプシロンよりももっと大きくしないと現実的に役に立たない。


単体テストで等値判定をする場合でも、プログラムの明白な誤りを検知することが目的ならCはそれほど小さくなくても充分だろうし、「この計算アルゴリズムならこれくらいの精度が出るべき」というのをテストしたいならそれに合わせてCやcを決める必要がある。
あと0との等値判定は別に扱う必要があって、xが0と等しいかはどうかは別の定数C'を使って
|x| < C'
で判定する。この定数C'の値をいくつにするかも状況次第。

計算機イプシロンの定義

そもそも計算機イプシロンには二種類の異なる定義がある。しかも、その二つがごっちゃに使われている。

計算機イプシロン定義:その1

定義1. 1.0より大きい最小の浮動小数点数と1.0との差

例えばIEEE754倍精度の場合、1.0より大きい最小の数(1.0の次の数)は

(1.0000000000000000000000000000000000000000000000000001)2×20

である(右下の添字「2」は2進数表記していることを表す)。よって計算機イプシロン

eps = (0.0000000000000000000000000000000000000000000000000001)2×20
= 2-52≒ 2.220446049250313e-16

となる。この値になるのは、IEEE754倍精度では仮数部が52ビットで、正規化数の場合は小数点以下の部分が52桁になるから。
この値は、問題にしている浮動小数点数システムが、どのくらい近くにある二つの数を区別して扱っているかを表している。つまり一種の分解能・解像度のようなものと解釈できる。ただし指数の大小に応じて、隣の数との隙間も大きくなったり小さくなったりするので、この値は「相対的な」分解能を表している。

計算機イプシロンの定義:その2

定義2. 「1.0 + eps > 1.0」が真となる最小の浮動小数点数eps

こちらの定義の方が良く使われている気がする。ただし正確には定義2のepsの値は、丸めモードによって違ってくる。「上への丸め」モードの場合なら、eps = 2-1074(最小の正の浮動小数点数)となるし、「下への丸め」モードと「0への丸め」モードの場合、1.0+epsが1.0の次の数以上の大きさにならないと切り捨てにより1.0になるので、epsの値は定義1の値と等しくなる。
おそらく定義2.は標準の丸め(最も近い点への丸め)の場合を考えている。その場合、eps=2-53≒1.1102230246251565e-16となる(定義1の値の半分)。
:正確には2-53ではなく「2-53より大きい最小の数」。なぜなら1.0 + 2-53の結果は、偶数点への丸めにより1.0になるから。ただし、計算してみたら何か挙動が変な感じになる。

main(){
  double x = 1.1102230246251565e-16; // 2^(-53)
  double y = 1.1102230246251568e-16; // 2^(-53)より大きい最小の数
  printf("(1) %30.25e\n", 1.0 + x);   
  printf("(2) %30.25e\n", 1.0 + 1.1102230246251565e-16); // (1)と(2)は同じはず
  printf("(3) %30.25e\n", 1.0 + y);
  printf("(4) %30.25e\n", 1.0 + 1.1102230246251568e-16); // (3)と(4)は同じはず
}

実行結果

(1) 1.0000000000000000000000000e+00
(2) 1.0000000000000000000000000e+00
(3) 1.0000000000000000000000000e+00
(4) 1.0000000000000002220446049e+00

なぜか? 拡張倍精度のせい? 最適化のせい?

定義1と定義2の混同

よくあるのが、説明の上では定義2を使っているのに、実際の数値は定義1の値を使っている、というもの。例えばRubyのマニュアルには

EPSILON
1.0 + Float::EPSILON != 1.0 となる最小の値
constant Float::EPSILON (Ruby 1.9.3)

と、定義2と同値の説明をしているけど、実際の値は定義1のものを使っている。
そのためEPSILONよりも小さい値でも不等式が成立する。

irb(main):001:0> Float::EPSILON
=> 2.22044604925031e-16
irb(main):002:0> 1.0 + (Float::EPSILON * 0.51) != 1.0
=> true

(定義2の計算機イプシロンの意味)

(かなり脱線するが、手間がかかったので一応残す)
定義2の計算機イプシロンがいったい何を表す数値なのか良くわからなかった。で、その答えになることがDavid Goldberg「What Every Computer Scientist Should Know About Floating-Point Arithmetic」に書いてあった。
この文章での計算機イプシロンの説明は一見すると定義1とも定義2とも違っている。この文章では計算機イプシロンとは、
実数をそれに最も近い浮動小数点に丸めたときに生じる相対誤差の上限
だと説明している。でも実はこの値は定義2の値と等しい。なぜなら

  1. 丸め誤差が最も大きくなるのは、実数の値が浮動小数点で表せる二つの値のちょうど中間の場合である。このとき生じる丸め誤差の大きさは、浮動小数点数の指数がeで浮動小数点数が1.??...??×2eの形のとき、2e-53となる(誤差は1.??...??の次の桁の丸めで発生するので)。
  2. 指数がeとなるxは、1.0×2e≦x<2.0×2eの範囲。
  3. よって丸めによる相対誤差が最も大きくなるのは1.0×2eのときで、2e-53/(1.0×2e)=2-53

となり、定義2の値と等しくなる。
e=0の場合で考える。丸めによる相対誤差が最も大きくなるのは、真の値が1.0と1.0の次の値のちょうど中間のとき。中間の値を求めるには、「1.0 + x」が1.0へ切り下げになる場合と1.0の次の値へ切り上げになる場合の境界を求めればよい。これは1.0 +x > 1.0が真になる値と偽になる値の境界を求めることに等しい。で、このとき求めたxは、丸めの絶対誤差の値に等しく、同時に丸めの相対誤差の値になっている。
ということで、定義2の計算機イプシロンは丸めで生じる相対誤差の上限を表している。

ついでに「What Every Computer Scientist Should Know About Floating-Point Arithmetic」で計算機イプシロンについての説明が出ていた節「Relative Error and Ulps」を引用しておく。何かわかりにくい。

相対誤差とulp

丸め誤差浮動小数点演算につきものなので、この誤差を測る方法を持っているのは重要である。β=10、p=3という浮動小数点フォーマットを考える[引用者注:βは基数、pは仮数部の桁数。ここの例なら10進で仮数部3桁]。なおこのセクションを通じてこのフォーマットを使う。浮動小数点演算の結果が3.12×10-2で、無限の精度で計算したときの答えが0.0314なら、最終桁を単位として2の誤差があるのは明らかである。同様に、実数0.031459が3.14×10-2で表される場合、最終桁を単位として0.159の誤差がある。一般に、zを表すのに浮動小数点数d.d…d×βeを使った場合、最終桁を単位として|d.d…d - (z/βe)|βp-1の誤差になる。「最終桁を単位として(unit in the last place)」の略称としてulpという語を使うことにする。正確な結果に一番近い浮動小数点数を演算結果にしても、0.5ulpの誤差がありえる。浮動小数点数とそれが近似している実数との差を測る別のやり方は相対誤差である。相対誤差は単に、二つの数の差を実数の方で割ったものである。例えば、3.14159を3.14×100で近似したとき相対誤差は 0.00159/3.14159 ≒ 0.0005である。

0.5ulpに対応する相対誤差を計算するため、次のことを見る。ある実数をそれにできる限り近い浮動小数点数d.dd…dd×βeで近似したとすると、誤差は最大で0.00…00β'×βeである。ただしβ'=β/2で、この浮動小数点数はp桁の仮数部を持ち、誤差の仮数部にはp個の0がある。この最大誤差は((β/2)β-p)×βeに等しい。 この誤差は絶対誤差としてはd.dd…dd×βeの形のどの数でも同じだが、この形の数が取る値の範囲はβeからβ×βeなので、相対誤差の範囲は((β/2)β-p)×βeeと((β/2)β-p)×βee+1の間となる。つまり
(2)   \frac{1}{2}\beta^{-p}\leq \frac{1}{2}\mbox{ulp} \leq \frac{\beta}{2}\beta^{-p}
である。つまり0.5ulpに相当する相対誤差は因子β倍の範囲で変化する。この因子の大きさをwobbleと呼ぶ。(2)の上限をε=(β/2)β-pとおけば、ある実数を最も近い浮動小数点に丸めたときの相対誤差はけっしてεを越えないと言うことができる。このεは計算機イプシロンと言われる。
上の例で、相対誤差は0.00159/3.14159 ≒ 0.0005である。こういう小さい数を避けるため普通、相対誤差はεを掛けた形で書かれる。εはこの場合ではε=(β/2)β-p=5(10)-3 = 0.005なので、相対誤差は(0.00159/3.14159)/0.005)ε≒0.1εと表される。
ulpと相対誤差の違いを説明するため、実数x=12.35を考える。これは\tilde{x}=1.24×101で近似される。誤差は0.5ulpで、相対誤差は0.8εである。次に8\tilde{x}を計算する。正確な値は8x=98.8で、計算した値は8\tilde{x}=9.92×101である。誤差は4.0ulpになったが、相対誤差は0.8εのままである。ulpで測った誤差は八倍になったのに、相対誤差は同じである。一般に、基数がβのとき、同一の相対誤差をulpで表すとβ倍の範囲で揺れ(wobble)を起こす。逆に上の式(2)で示したように、同じ0.5ulpの誤差でも相対誤差にするとβ倍の範囲で揺れが起こる。
丸め誤差を測る最も自然なやり方はulpで測ることである。例えば、最も近い浮動小数点数に丸めると、誤差は0.5ulp以下である。しかし、様々な公式における丸め誤差を分析する場合、相対誤差の方が良い測り方になる。そのことの良い例が、このセクションの定理9の分析である。εは、最も近い浮動小数点数への丸めの影響をwobble βの範囲で過大評価しうるので、公式の誤差の評価はβが小さい計算機の方が小さくなる。
丸め誤差の大きさのオーダーだけに興味がある場合は、ulpとεは取り違えて使うことができる。多くてもβ倍の違いしかないので。例えば、ある浮動小数点数n ulpの誤差がある場合、誤差の混入している桁数はlogβnである。もし計算結果の相対誤差がnεなら、
(3)   誤差に汚染された桁数 ≒ logβn

浮動小数点数の等値判定に計算機イプシロンを(どう)使うべきか

ようやく本題。「相対誤差が計算機イプシロン以下である」ということを調べるのは、何か変だと思う理由を書く。まず1.0とそれと数直線上で隣り合った数を比較した場合、どれくらいの相対誤差があるのかを考える。大きい方の隣りには、

(1.0000000000000000000000000000000000000000000000000001)2×20

があり、小さい方の隣には、

(1.1111111111111111111111111111111111111111111111111111)2×2-1 =
(0.11111111111111111111111111111111111111111111111111111)2×20

がある。1.0と大きい方との相対誤差は2-52で、小さい方との相対誤差は2-53になる。これは1.0の代わりに1.0×2eで考えても同じ結果になる。さらに、隣り合った数の間での相対誤差はこの二つの値の間に入ることもわかる。
それから、2-52は定義1の計算機イプシロンeps1に等しく、2-53は定義2の計算機イプシロンeps2に等しい。つまり

eps2 ≦ 隣りの数との相対誤差 ≦ eps1

が成り立つ。
そうすると定義2の値を使った場合、すぐ隣りの数と比較しても相対誤差は計算機イプシロン以下にはならない(1.0×2eの形の数の場合は例外)。つまり定義2の値を使った場合、単に等しいかどうかを判定しているのとほとんど変わらないことになる。
一方、定義1の値を使うと、相対誤差が計算機イプシロン以下になるのは、自分自身と比較した場合あるいはすぐ隣りにある数と比較した場合だけになる。よって定義1の値を使った場合、厳密な等値比較よりちょっとだけゆるめた等値比較をしていることになる。


で疑問なのは、そもそも何で浮動小数点数の等値比較にこのやり方を使わないといけないのかということ。特に計算機イプシロンの説明として定義2の方を使っている場合、ほとんど意味のないことを推奨しているように思える。普通の等値比較とほとんど変わらないから。実際にはライブラリで与えられている計算機イプシロンの値が定義1の方だったおかげで、単なる等値比較とは違うことをやっているという結果になっているかもしれない。でもその場合、間違った値による比較方法を使っていることになって、ますますその方法を使う根拠がわからなくなる。
計算機イプシロンとして定義1を与えた場合でもやっぱり、その値を等値比較に用いる理由がわからない(ちょっとゆるめの等値比較が欲しいという理由以外には)。
ただ計算機イプシロンの値が等値比較とまったく無関係だとは思わない(特に定義2の方)。相対誤差が計算機イプシロンのc倍以下なら等しいとする、みたいな比較をおこなうことは理解できる。

浮動小数点数の等値判定はどのようにすべきか

定数cと計算機イプシロンEを使って
\frac{|x-y|}{\max(|x|,|y|)} < cE
によってxとyが等しいかを判定するとして、問題は定数cの選び方。厳密な精度判定が必要なわけではないなら適当に大きめの値を使う(計算によっては大きな桁落ち誤差が割と簡単に発生する。その場合に許容値を大きくするか計算方法を変えるかは状況次第)。
R5RSの6.2.5節、R6RSの11.7.4.3節に次のようにある。

小さな不正確さが結果に影響を与えうるので、結果は信頼できないかもしれない。特に=とzero?で。疑わしければ、数値計算の専門家に相談せよ。

また、The Test Tools: floating-point comparison algorithmsでは、次のように言われている。

浮動小数点比較アルゴリズム

(...)浮動小数点数u、vと許容値eに対して(...)

|u - v| / |u| ≦ e  かつ   |u - v| / |v| ≦ e (1')
|u - v| / |u| ≦ e  または |u - v| / |v| ≦ e (2')
許容値の選択の仕方

領域特有の要求が無い場合、許容値として、比較する値の「相対丸め誤差」の予測上限値の和を選ぶことができる。
(...)最初の二つ[型昇格と数値演算]の相対丸め誤差は「計算機イプシロン」の1/2を越えないことが証明できる。(...)
最後に、1/2*epsilonルールは推移的でないことに注意。言い換えると、二つの数値演算を組み合わせたとき、2*1/2*epsilonをはるかに越える丸め誤差が生じるうる。結局、許容値をどのように選択するかの一般ルールはなく、常識および領域特有/問題特有の知識を使って許容値を決める必要がある。