#まえがき
二次方程式の解き方は覚えていますか？
a x ^2 + b x + c = 0
の解は
x = {-b \pm \sqrt{b^2 - 4ac}\over 2a}
ですね。
これをPythonで解いてみましょう。
解の公式をそのまま出力するものが以下の関数です。
import numpy
def qeq(a, b, c):
return((-b + numpy.sqrt(b**2 - 4 * a * c))/(2 * a),
(-b - numpy.sqrt(b**2 - 4 * a * c))/(2 * a))
例えば
x^2 + 5x + 6 = 0
を計算すると
print(qeq(1, 5, 6)) # -2.0, -3.0
という答えが出ますね。
#何が問題になるのか
x^2 + 1.000000001 x + 0.000000001 = 0
という二次方程式を考えてみましょう。
x = -0.000000001, -1
の二つが解となります。
同様の関数で計算してみると・・・？
print(qeq(1, 1.000000001, 0.000000001))
# -1.0000000272292198e-09, -1.0
(-1.0000000272292198e-09は
-0.0000000010000000272292198のこと。)
9桁目でずれが生じてしまいました。
有効数字8桁はちょっと少なすぎ。
桁落ちとは？
有効数字という概念があります。
ざっくりいうと小数の取り扱い時に計算の正確さを保証できる桁数のこと。
1.001なら有効数字4桁、
1.000でも有効数字4桁。
1は有効数字1桁、
0.001も有効数字1桁。
ここまでいいでしょうか。
ここからが問題。
1.001 - 1.000
の有効桁数は？
4桁から4桁を引いたんだから、 4桁！って言いたいところですが、
= 0.001
で1桁。これが桁落ちです。
桁落ちは極めて近い数同士での引き算で起こります。例えば、
2.001 - 1.000 = 1.001
なら4桁のまま。
何がダメだった？
今回の場合は
-b + \sqrt{b^2 - 4ac}
ここが問題。
-bは-1.000000001、$\sqrt{b^2 - 4ac}$は0.999999999...
合わせて0.000000002くらい。
二つの値が近すぎたため、9桁の有効数字が失われました。
これで解の有効数字は8桁あったということは、元は有効数字が17桁ほどあったということでしょうか。
#解決策
計算方法を変えてみましょう。
足し算では桁落ちは起こらないので、解の一つはそのまま解の公式で求めることができます。
$(i)b\geq 0の時$
x = {-b -\sqrt{b^2 - 4ac}\over 2a}
$(ii)b< 0の時$
x = {-b +\sqrt{b^2 - 4ac}\over 2a}
もう一つの解を求めるために、以下の公式が役に立ちます。
a x ^2 + b x + c = 0
の解を$\alpha, \beta$とすると
\alpha \beta = {c \over a}
証明は省略。
これでもう一つの解を求めることができますね。
#実装してみよう
import numpy # numpy.sign()という関数を使う。引数の正負を判断する。
def qeq(a, b, c):
alpha = (-b - numpy.sign(b) * numpy.sqrt(b**2 - 4 * a * c))/(2 * a)
beta = (c / a) / alpha
return(alpha, beta)
print(qeq(1, 1.000000001, 0.000000001))
# (-1.0, -1e-09)
桁落ちがなくなったので正確な値が計算できました。
まとめ
値が極めて近い数字同士を引き算すると情報が失われ、精度が落ちる。なんとか工夫して回避しよう。
『機械学習のエッセンス』という本を参考にしました。引き続きこの本についてまとめていきます。
前回：丸め誤差に気をつけよう
https://qiita.com/NNNiNiNNN/items/42107479407d9bb14419
次回：計算過程に気をつけよう
https://qiita.com/NNNiNiNNN/items/2f17f874f8a234bd10ad