Python高精度联立方程求解器需求：替代np.linalg.solve方案咨询

需要在Python环境下实现20位小数精度的联立方程求解，当前使用np.linalg.solve(X,Y)和scipy.linalg.solve的精度无法满足要求，即使小规模5元方程组也存在明显误差（如下列代码及运行结果所示），可接受代码运行速度较慢。

用户测试代码 #

import numpy as np

n=5

Y = np.random.rand(n)
X = np.tile(np.array(range(1,n+1)),n)
X = X.reshape((n,n),order='F')

for c in range(n) : 
  X[:,c] = X[:,c]**c

A = np.linalg.solve(X,Y)

predicted_Y = X@A

table = [(y,pred_y,y-pred_y) for y,pred_y in zip(Y,predicted_Y)]

print('y                        predicted_y              difference')
for c1,c2,c3 in table : 
  print(f"%.20f | %.20f | %.20f" % (c1, c2, c3))

测试结果（误差示例） #

y                        predicted_y              difference
0.68935295599312118586 | 0.68935295599312118586 | 0.00000000000000000000
0.72899266151307307027 | 0.72899266151307240413 | 0.00000000000000066613
0.18770646040141103494 | 0.18770646040141256150 | -0.00000000000000152656
0.02144867791874205398 | 0.02144867791873661389 | 0.00000000000000544009
0.54517050144884360297 | 0.54517050144883372198 | 0.00000000000000988098

1. mpmath #

mpmath是专注于高精度数值计算的库，可自定义精度位数，完全覆盖20位小数的需求，内置线性方程组求解接口。

示例代码 #

import mpmath as mp

# 设置精度为50位（远高于20位需求）
mp.mp.dps = 50

n = 5
# 生成高精度矩阵与向量
Y = [mp.mpf(np.random.rand()) for _ in range(n)]
X = mp.matrix(n, n)
for c in range(n):
    for r in range(n):
        X[r, c] = (r+1)**c

# 求解方程组
A = mp.lu_solve(X, Y)

# 计算预测值并输出
predicted_Y = X * A
print('y                        predicted_y              difference')
for y, pred_y in zip(Y, predicted_Y):
    diff = y - pred_y
    print(f"{mp.nstr(y, 20)} | {mp.nstr(pred_y, 20)} | {mp.nstr(diff, 20)}")

2. SymPy #

SymPy是符号计算库，支持精确代数运算，可先求解方程组的精确符号解，再转换为指定精度的浮点数。

示例代码 #

import sympy as sp
import numpy as np

n = 5
# 定义符号变量
A_vars = sp.symbols(f'a0:{n}')

# 构建方程组
Y = np.random.rand(n)
X = np.tile(np.array(range(1, n+1)), n).reshape((n, n), order='F')
for c in range(n):
    X[:, c] = X[:, c]**c

equations = []
for i in range(n):
    expr = sum(X[i][j] * A_vars[j] for j in range(n)) - Y[i]
    equations.append(sp.Eq(expr, 0))

# 求解并转换为20位精度浮点数
solution = sp.solve(equations, A_vars)
A = [sp.N(solution[var], 20) for var in A_vars]
predicted_Y = [sum(X[i][j] * A[j] for j in range(n)) for i in range(n)]

# 输出结果
print('y                        predicted_y              difference')
for y, pred_y in zip(Y, predicted_Y):
    diff = y - pred_y
    print(f"{sp.N(y, 20)} | {sp.N(pred_y, 20)} | {sp.N(diff, 20)}")

3. 标准库decimal模块 #

无需第三方库，通过设置高精度上下文，手动实现高斯消元算法即可完成求解。

示例代码（高斯消元实现） #

from decimal import Decimal, getcontext
import numpy as np

# 设置精度为50位
getcontext().prec = 50

n = 5
# 转换为Decimal类型数据
Y = [Decimal(np.random.rand()) for _ in range(n)]
X = [[Decimal((r+1)**c) for c in range(n)] for r in range(n)]

# 高斯消元求解函数
def gauss_solve(matrix, rhs):
    n = len(rhs)
    aug = [row + [rhs[i]] for i, row in enumerate(matrix)]
    
    for col in range(n):
        # 选取主元
        pivot = max(range(col, n), key=lambda i: abs(aug[i][col]))
        aug[col], aug[pivot] = aug[pivot], aug[col]
        
        # 归一化主行
        factor = aug[col][col]
        for j in range(col, n+1):
            aug[col][j] /= factor
        
        # 消去其他行
        for i in range(n):
            if i != col and aug[i][col] != 0:
                factor = aug[i][col]
                for j in range(col, n+1):
                    aug[i][j] -= factor * aug[col][j]
    
    return [row[-1] for row in aug]

# 求解并计算预测值
A = gauss_solve(X, Y)
predicted_Y = [sum(X[i][j] * A[j] for j in range(n)) for i in range(n)]

# 输出结果
print('y                        predicted_y              difference')
for y, pred_y in zip(Y, predicted_Y):
    diff = y - pred_y
    print(f"{y.normalize():.20f} | {pred_y.normalize():.20f} | {diff.normalize():.20f}")

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