NumPy 数组中高浮点错误的解决方法

Question

我有代表 24 位 A/D 转换器响应的数据。我正在将数据设置为进入下面函数中的多项式回归，但我得到了有害的浮点误差，以至于它使得我的回归系数稍后完全无法使用。

数据经过标准化后，范围在 -1 和 1 之间。我得到的浮点误差大约为 0.15，最高可达四度。由于无关的原因，我无法降低回归项的阶数，但是感谢 @Robert Dodier 指出，一般来说，当添加更高的阶数时，这种方法会变得相当不稳定。

有没有办法减少浮点误差？速度和内存都不是这里的主要目标：如果可能的话，准确性需要接近精确。数据集非常小（约 400 行）。

我还尝试使用 gmpy2 库的“mpfr”类，但这似乎没有什么区别。为了完整起见，我也包含了该版本。

在此函数之后，我只需将数据输入 QR 分解回归即可：

Q, R = np.linalg.qr(X)
beta = np.linalg.inv(R).dot(Q.T).dot(y)

处理数据后，我还检查了浮点误差是否可以归因于回归，但是当我运行 SVD 式回归时，我得到的系数与 QR 式几乎相同，所以这不是问题.

以下是这些功能：

import numpy as np
def load_data(path):
    data = np.loadtxt(path, delimiter=',')

    # Features
    A = data[:, 1:3]

    # Normalize the data
    A /= (2 ** 24)

    # Extract columns P and T
    P = A[:, 0]
    T = A[:, 1]

    # Compute new columns based on P and T
    col_ones = np.ones_like(P, dtype=float)
    T2 = T ** 2
    T3 = T ** 3
    PT = P * T
    PT2 = P * (T ** 2)

    # ... continue making features

    # Combine the new columns into a new array
    A = np.column_stack((col_ones, T, T2, T3, P, PT, PT2)) # ... add higher degree features
    b = data[:, 0]  # ref

    return A, b

import gmpy2
from gmpy2 import mpfr
def load_data_gmpy(path):
    data = np.loadtxt(path, delimiter=',', skiprows=15)

    # Features
    b = data[:, 0]  # ref
    P = data[:, 1]
    T = data[:, 2]

    gmpy2.set_context(gmpy2.context())
    with gmpy2.local_context() as ctx:
        ctx.precision = 2000

        # Normalize the data
        P = [mpfr(p) / (2 ** 24) for p in P]
        T = [mpfr(t) / (2 ** 24) for t in T]

        # Compute new columns based on P and T
        col1 = np.ones_like(P, dtype=float)
        T2 =   np.array([float(t ** 2) for t in T])
        T3 =   np.array([float(t ** 3) for t in T])
        PT =   np.array([float(p * t) for p, t in zip(P, T)])
        PT2 =  np.array([float(p * (t ** 2)) for p, t in zip(P, T)])
        # ... continue making features

        # Convert the results back to numpy arrays
        # add higher-degree features here...
        A = np.column_stack((col1, np.array(T, dtype=float), T2, T3, np.array(P, dtype=float), PT, PT2)) 

        return A, b

Answer 1

通常有 2 种方法来防止浮点错误：

1. 使用整数/小数表示（a/b），并且仅在执行完所有加法、减法、乘法等后才在最后一步转换为浮点数。由于您使用的是 24 位 ADC，我怀疑这将是最好的方法：将 24 位数据读入更大的缓冲区（例如 64 位，但如果它们是您需要的有符号值，请小心一些位移以防止弄乱标志）。仅使用整数运行回归，并且仅在尝试向用户显示值时才使用浮点数

2. 如果您必须转换为浮点数并使用浮点数，唯一真正的选择是使用 double 而不是 float （因为您不受内存限制，这应该是一个简单的测试）并在加法、乘法之前对浮点数进行排序您还可以研究卡汉求和算法，这是做同样事情的更有效方法，但您说性能不是问题。我还没有尝试过对 Kahan 进行基准排序，但从之前的经验来看，它们非常相似。

就像其他人提到的那样，发布您的数据（因为您说它很小，大约 400 行应该可以在此处发布原始数据）和具有可重现错误的代码，那么推荐适当的修复或解决方案会容易得多。