pythonのmath.acos()関数でValueError: math domain error

pythonのmath.acos()関数でValueError: math domain errorが出た話。 浮動小数点演算で誤差が乗って、定義域の範囲外の数を引数に入れていた。

2次元座標の上の3つの点から、なす角を求めようとした。
より正確に言えば、2つのベクトル \boldsymbol{a, b}のなす角θを求めようとした。  \boldsymbol{a \cdot b} = | \boldsymbol{a}| | \boldsymbol{b}| \cos \thetaなので、
ベクトルa, bの内積と、aの長さとbの長さを求めて、最後に逆三角関数arccosを使うという方法で角度を求めた。 以下のページに詳しく説明されているので、詳細は割愛する。

3点からなる角度を求める 画像処理ソリューション

Pythonのversionは3.7.6。

math domain error 発生状況

import numpy as np
import math

def calc_angle(a, b):
    cosine_val = (a[0] * b[0] + a[1] * b[1]) / (math.sqrt(a[0]**2 + a[1]**2) * math.sqrt(b[0]**2 + b[1]**2))
    ans = math.degrees(math.acos(cosine_val))
    return ans

上記の関数を考える。 math.acosラジアンの角度を返すので、ラジアンから度数法に変換するために、最後にdegreesを適用してから値を返している。

例えば、原点から真横に行くベクトルと、斜め右上にいくベクトルのなす角度は……

a = [1, 0]
b = [1, 1]
calc_angle(a, b)

# 45.00000000000001

45度だ。

2つ目の例として、、原点から真横に行くベクトルと、斜め左上にいくベクトルのなす角度は……

a = [1, 0]
b = [-1,  1]
calc_angle(a, b)

# 135.0

135度だ。

いいですね。

では、2次元座標上に適当にx, yを取り、その中点をmidとする。 midからxに向かうベクトルと、midからyに向かうベクトルのなす角はいくつだろうか。

下の図を見れば分かる通り、これは180度である。

ではやってみよう。

x = np.array([3, 4])
y = np.array([-2, -7])
mid = (x+y)/2
calc_angle(x-mid, y-mid)

# 180.0

期待通り。

x = np.array([0.123, 4.567])
y = np.array([3.141, 2.718])
mid = (x+y)/2
calc_angle(x-mid, y-mid)

# 179.99999914622634

数値誤差は出たけど、180度だ。

x = np.array([3.14159265, 2.71828183])
y = np.array([1.23456789, 9.8765432])
mid = (x+y)/2
print(calc_angle(x-mid, y-mid))
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-27-58af6d8dd291> in <module>
      2 y = np.array([1.23456789, 9.8765432])
      3 mid = (x+y)/2
----> 4 print(calc_angle(x-mid, y-mid))

<ipython-input-2-3bbda96d5884> in calc_angle(a, b)
      1 def calc_angle(a, b):
      2     cosine_val = (a[0] * b[0] + a[1] * b[1]) / (math.sqrt(a[0]**2 + a[1]**2) * math.sqrt(b[0]**2 + b[1]**2))
----> 3     ans = math.degrees(math.acos(cosine_val))
      4     return ans

ValueError: math domain error

ValueError: math domain errorというエラーが出た。なぜだろうか?

原因

arccosの引数を見るとエラーの理由がわかる。関数の途中にprint(cosine_val)と入れてみると、以下のようになる。

-1.0000000000000002

つまり、エラーが発生したしくみは以下の通りである。

  • math.acosの引数に取れるのは-1以上1以下の引数で、それ以外だとDomain Errorを返す
  • 引数は、無限精度の演算であれば-1になるはずだが、浮動小数点の演算で誤差が発生し、-1よりわずかに小さい値になった
  • その結果、domain errorが発生した

再現しようと思っていろんな数で試したけど、xやyの桁数を多くしたほうが良いみたい? xやyとして、0から1までの乱数を適当に代入してやると、高確率で発生する。


domain errorというのは、関数の定義域から外れたエラーのことであろう。
math --- 数学関数 — Python 3.8.3rc1 ドキュメントには ValueError: math domain errorのことは特に記載がない。

ついでにmathモジュールの関数について、いくつか同様にdomain errorを発生させてみた。

sqrt(平方根)。

math.sqrt(-1)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-34-2037b8d41d70> in <module>
----> 1 math.sqrt(-1)

ValueError: math domain error

log(対数関数)。

math.log(-1)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-35-8f964a7914e2> in <module>
----> 1 math.log(-1)

ValueError: math domain error

acosh(逆双曲線関数)。

math.acosh(0.99) # 逆双曲線関数arccoshの定義域は1以上の実数
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-38-c0f54ec4e970> in <module>
----> 1 math.acosh(0.99)

ValueError: math domain error

factorial(階乗関数)。

math.factorial(-5)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-33-a46d876612ec> in <module>
----> 1 math.factorial(-5)

ValueError: factorial() not defined for negative values

いや、なんでこれだけdomain errorじゃないねん!

math.factorial(2.34)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-32-044b0561db2f> in <module>
----> 1 math.factorial(2.34)

ValueError: factorial() only accepts integral values

いや、なんでfactorialだけdomain errorじゃないねん!!


浮動小数点の演算でうまく行かなかった話では、以下も参照:

linus-mk.hatenablog.com

それでは。