おっぱいそん!

pythonを数値計算に使うときのテクニックとかをまとめていきたい。

複素数型が実数に近かったら、実数型にする

complex型のndarrayが与えられた時に、もし複素数の部分が小さいなら、floatなどの実数の型で近似する方法。
そういう関数、ありそうだと思ってググったらNumpyの関数にnumpy.real_if_closeというのがあった。

numpy.real_if_close(a, tol=100)

この関数はNdarray aの虚部a_imagを見て、それがa_imag≦マシンイプシロン*tol(Tolerance)なら、虚部を無視した実部だけのNdarrayを返す。
a = a_real + i*a_imag → a = a_real

マシンイプシロン
np.finfo(np.float).eps
で調べることができ、Pythonでは普通、2.22044604925e-16になっているみたいです。
tolはデフォルトでは100に設定されているので、a_imag≦2.22044604925e-14であれば虚部を無視することになります(tolはオプションで変更できます)。


このように、np.real_if_closeは虚部の絶対値をみて、小さいかどうかを判断する仕様になっている。
しかし、普通aの絶対値が大きくなると、誤差の絶対値も大きくなるので、虚部/絶対値が小さいかどうかを見るほうが、良い場合が多いと思う。:

print np.finfo(np.float).eps
a_real = np.array([1.0e4])
a_imag =  np.array([3e-12j])
a = a_real + a_imag
print a.dtype
a = np.real_if_close(a, tol=100*np.absolute(a))
print a.dtype

2.22044604925e-16
complex128
float64

Numpyは欲しい関数がたいてい用意されているのですごい!
http://docs.scipy.org/doc/numpy/reference/generated/numpy.real_if_close.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.absolute.html

2015.6.26 追記
オプションtolで1より小さい値を設定すると、real_if_closeは配列を常に実数型にしてしまうようです(本来は複素数型のままにして欲しい)。
そもそも、tolが1より小さいとa_imagがゼロでない値を持っていたら、最小値はマシンイプシロンなので、常に複素数のままになるべきですが、たぶんそういうtolを想定して関数が作られておらず、バグってしまうのだと思います。
上のtol=100*np.absolute(a)という設定はaの絶対値が大きい場合を想定していますが、逆に絶対値が小さい場合は使うときに注意してください。