NumPy数组幂运算浮点误差过大的解决办法求助

我有一组24位A/D转换器的响应数据，准备用于多项式回归，但出现严重浮点误差，导致回归系数完全无法使用。

• 数据归一化后范围在-1到1之间，四次项的浮点误差约为0.15
• 因外部限制，无法降低回归项的阶数
• 目标是保证精度，速度和内存都不是问题
• 尝试过使用gmpy2库的mpfr类，但未取得效果，附上相关代码
• 分别用QR分解回归和SVD风格回归验证，得到的系数几乎一致，说明问题不在回归方法本身

相关代码 #

原始数据加载与特征构造 #

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

gmpy2版本数据加载 #

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

回归实现代码 #

def do_qr_regression(X, y):
    # QR decomposition
    Q, R = np.linalg.qr(X)
    beta = np.linalg.inv(R).dot(Q.T).dot(y)

    # predict using coefficients
    yhat = X.dot(beta)

    return beta, yhat


def do_svd_regression(X, y):
    # calculate coefficients
    beta = np.linalg.pinv(X).dot(y)
    # predict using coefficients
    yhat = X.dot(beta)

    return beta, yhat


if __name__ == '__main__':
    data, ref = load_data("path\\to\\file.txt")

    coefs_qr, values = do_qr_regression(data, ref)
    residuals_qr: np.array = values - ref

    coefs_svd, values = do_svd_regression(data, ref)
    residuals_svd: np.array = values - ref

    plt.plot(residuals_qr, 'r-', label='QR')
    plt.plot(residuals_svd, 'g-', label='SVD')

    plt.xlabel('Index')
    plt.ylabel('Residuals')  
    plt.title('Comparison of Residuals')

    plt.legend()
    plt.grid()
    plt.show()

    coefs_qr = np.array(coefs_qr)
    coefs_svd = np.array(coefs_svd)

    # Combine the coefficients into a single array for easier handling
    coefficients = np.column_stack((coefs_qr, coefs_svd))

    fig, ax = plt.subplots()

    # Hide the axes
    ax.axis('tight')
    ax.axis('off')

    table = ax.table(cellText=coefficients, colLabels=['QR', 'SVD'], loc='center')

    plt.show()

1. 修复gmpy2使用方式：全程保留高精度 #

你的gmpy2代码无效的核心原因是中途将mpfr类型转成了float，导致高精度计算的结果被截断回普通双精度。要全程用mpfr处理所有计算，包括后续的回归步骤（numpy不支持mpfr，需要改用支持高精度的线性代数工具）。

修改后的gmpy2特征构造示例（不转float）：

import gmpy2
from gmpy2 import mpfr
def load_data_gmpy_fixed(path):
    data = np.loadtxt(path, delimiter=',', skiprows=15)
    b = data[:, 0]
    P = data[:, 1]
    T = data[:, 2]

    ctx = gmpy2.context()
    ctx.precision = 2000
    gmpy2.set_context(ctx)

    # 归一化并保留mpfr类型
    P = [mpfr(p) / mpfr(2**24) for p in P]
    T = [mpfr(t) / mpfr(2**24) for t in T]

    # 构造特征，全程用mpfr
    col_ones = [mpfr(1)] * len(P)
    T2 = [t**2 for t in T]
    T3 = [t**3 for t in T]
    PT = [p*t for p,t in zip(P,T)]
    PT2 = [p*t**2 for p,t in zip(P,T)]

    # 转成sympy Matrix（方便后续高精度线性代数运算）
    from sympy import Matrix
    A = Matrix([col_ones, T, T2, T3, P, PT, PT2]).T
    b_mat = Matrix([mpfr(val) for val in b])

    return A, b_mat

然后用sympy的QR分解求解回归系数：

def qr_regression_high_precision(X, y):
    Q, R = X.QRdecomposition()
    beta = R.inv() * Q.T * y
    yhat = X * beta
    return beta, yhat

2. 用正交多项式替换普通幂次项，缓解矩阵病态性 #

多项式回归的高次幂项之间高度相关，会导致特征矩阵病态，放大浮点误差。改用正交多项式（如Legendre多项式）构造特征，能让特征矩阵列之间正交，大幅提升数值稳定性，即使使用普通双精度也能显著降低误差。

以T的特征为例，替换为Legendre多项式：

def legendre_features(T):
    # 假设T已归一化到[-1,1]
    L0 = np.ones_like(T)
    L1 = T
    L2 = (3*T**2 - 1)/2
    L3 = (5*T**3 - 3*T)/2
    L4 = (35*T**4 - 30*T**2 + 3)/8
    return L0, L1, L2, L3, L4

然后用这些正交特征替换原来的T^2、T3等项，构造特征矩阵后再做回归，误差会显著降低。

3. 优化回归求解的数值稳定性 #

即使使用普通双精度，也可以通过更稳定的求解方式降低误差：

• 不要直接对R求逆，改用np.linalg.solve(R, Q.T.dot(y))，solve函数比inv更稳定
• 对于SVD回归，numpy的pinv已经是稳定实现，但可以指定rcond参数来过滤过小的奇异值，避免噪声放大

修改后的QR回归函数：

def do_qr_regression_stable(X, y):
    Q, R = np.linalg.qr(X)
    # 用solve替代inv，数值更稳定
    beta = np.linalg.solve(R, Q.T.dot(y))
    yhat = X.dot(beta)
    return beta, yhat

4. 使用任意精度线性代数库 #

如果需要极致精度，可以使用sympy或mpmath这类支持任意精度的库，它们的矩阵运算完全基于高精度浮点数，能彻底避免普通双精度的浮点误差。

比如用mpmath实现回归：

import mpmath as mp
mp.mp.dps = 2000  # 设置精度为2000位小数

def mpmath_regression(X_np, y_np):
    # 转成mpmath矩阵
    X = mp.matrix(X_np.tolist())
    y = mp.matrix(y_np.tolist())
    # 用伪逆求解
    beta = mp.pinv(X) * y
    yhat = X * beta
    # 转回numpy数组（可选）
    beta_np = np.array([float(val) for val in beta])
    yhat_np = np.array([float(val) for val in yhat])
    return beta_np, yhat_np

内容的提问来源于stack exchange，提问作者javery