NumPy 高精度运算

Question

我正在使用numpy和pyfits来操作光谱，需要高精度计算（一个值可能高达10^12，需要保留8-10个小数位）。为此，"decimal"数据类型非常适合（float64不够精确），但不幸的是，numpy.interp不支持它：

File ".../modules/manip_fits.py", line 47, in get_shift
pix_shift = np.interp(x, xp, fp)-fp
File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1053, in interp
return compiled_interp(x, xp, fp, left, right)
TypeError: array cannot be safely cast to required type

我使用的代码的简化版本：

fp = np.array(range(new_wave.shape[-1]),dtype=Decimal)
pix_shift = np.empty_like(wave,dtype=Decimal)
      x = wave
  xp = new_wave
 pix_shift = np.interp(x, xp, fp)-fp

其中'wave'和'new_wave'是代表1D频谱的一维numpy数组。该代码需要将我的光谱沿着X轴(即波长)移动。

我最大的问题是在代码的后面，我会通过从所有光谱的总和构建的模板光谱来除以我的光谱，以分析差异，但由于我没有足够的小数位数，所以出现了舍入误差。有什么想法吗？

谢谢！

更新:

测试示例：

import numpy as np
from decimal import *
getcontext().prec = 12

wave = np.array([Decimal(xx*np.pi) for xx in range(0,10)],dtype=np.dtype(Decimal))
new_wave = np.array([Decimal(xx*np.pi+0.5) for xx in range(0,10)],dtype=np.dtype(Decimal))

fp = np.array(range(new_wave.shape[-1]),dtype=Decimal)
pix_shift = np.empty_like(wave,dtype=Decimal)

x = wave
xp = new_wave
pix_shift = np.interp(x, xp, fp)-fp

错误信息为：

Traceback (most recent call last):
  File "untitled.py", line 16, in <module>
    pix_shift = np.interp(x, xp, fp)-fp
  File "/usr/lib/python2.7/dist-packages/numpy/lib/function_base.py", line 1053, in interp
    return compiled_interp(x, xp, fp, left, right)
TypeError: array cannot be safely cast to required type

这是最接近使用适合格式的真实光谱数据的翻译。

更新2： 使用十进制数打印的我的光谱数据的一些典型值：

  18786960689.118938446044921875
  18473926205.282184600830078125
  18325454516.792461395263671875
  18400241010.149127960205078125
2577901751996.03857421875
2571812230557.63330078125
2567431795280.80712890625

我遇到的问题是，在对它们进行操作时，会产生四舍五入误差。例如，我通过对所有光谱进行求和来创建一个模板。然后我使用这个模板来对每个光谱进行归一化。例如：

Spectra = np.array([Spectrum1, Spectrum2, ...])
Template = np.nansum(Spectra, axis= 0)

NormSpectra = Spectra/Template

假设模板是恰当地代表了这颗星，我希望只返回光谱上的噪声。我尝试将每个光谱归一化为其总通量。

(Spectrum1 = Spectrum1/np.nansum(Spectrum1), ...) 

对于模板来说，使用Decimal可能会很好，但是会导致舍入误差更加严重。

对我来说，使用Decimal没问题，但是我需要“移动”我的光谱使得所有的光谱特征/线都能够对齐。

希望这样讲清楚了？

- jorgehumberto

---

“something like 8-10 decimal places on a value which might go as high as 10^12” 可能会对 float32 影响较大，但远未达到 float64 的极限。您可以提供一个示例以展示 float64 不够用的情况吗？您也可以尝试缩小问题的规模。 - jorgeca

Decimal不是numpy的数据类型对象，也不能转换为其一种类型，因此您永远无法将其用作数据类型。NumPy所做的是看到“哦，您提供了一个任意类，最好创建一个Python-Object数组”。这样做还算可以，但是在需要将其强制转换回NumPy数据类型以执行计算时就有问题了。从您的代码中我得到了一个更有用的错误消息：TypeError: Cannot cast array data from dtype('O') to dtype('float64') according to the rule 'safe'。这支持了这个理论。如果您创建的这个对象数组对scipy.interp1d也有效，那么我会感到非常惊讶，因为... - kampu

因为数据类型为 'O' 的数组不被视为数值类型，NumPy 将不会尝试使用它进行数值计算。 - kampu

很抱歉，我不知道如何从您提供的数据中重现舍入误差。顺便问一下，您是否尝试过使用dtype=np.longdouble的额外位数？此外，mpmath可能对您有所帮助（它是一个用纯Python编写的多精度浮点运算库，还可以使用快速编译后端）。 - jorgeca

你的数据有多少有效数字？问题仅在于舍入误差吗？如果你需要少于14个小数位且需要精确求和，你可以使用Karatsuba算法。 - unixmin

Answer 1

在典型用例中，双精度浮点数可以期望得到约15个有效数字。如何确保np.float64是正确的？如果您确定这不足够，可以尝试np.float128（又名np.longdouble）。但是，您的问题似乎比这更深刻：它似乎是一个病态问题（通常涉及大数除以小数）。这不是您想要的结果。增加精度应该在一定程度上解决问题，但您将遇到一些数据需要使用float256 / float512等来避免病态舍入误差。我建议您解释您的问题，而不是解决方案，这样我们就可以希望在每种情况下找到另一种解决方法（XY Problem）。

- user1940040