高精度求解非方阵线性方程组A*b=c的精度优化问题

我正在求解线性方程组A*b = c，其中矩阵A和c的定义如下：

A = np.block([[X_matrix - rho_0, np.zeros((3, 3)), np.zeros((3, 3))],
              [np.zeros((3, 3)), X_matrix - rho_1, np.zeros((3, 3))],
              [np.zeros((3, 3)), np.zeros((3, 3)), X_matrix - rho_2],
              [I, I, I]])

c = np.block([[np.zeros((3, 3))],
              [np.zeros((3, 3))],
              [np.zeros((3, 3))],
              [I]])

由于A是非方阵，我只能使用伪逆求解：

A_inv = np.linalg.pinv(A)
# Solve for b
b = np.matmul(A_inv, c)

但该方法精度不足，执行print(np.linalg.norm(np.matmul(A, b) - c))得到的结果为0.212871998042824。请问是否有更高精度的方法求解b？

我已尝试np.linalg.lstsq，但结果与np.linalg.pinv(A)一致；使用mpmath的A**-1时出现错误：ValueError: only powers of square matrices are defined。

高精度解决方案建议 #

• 切换高精度数据类型与mpmath伪逆：NumPy默认float64精度有限，可尝试用np.float128（环境支持前提下），或者用mpmath的高精度矩阵计算。注意mpmath非方阵伪逆需调用mp.pinv()而非A**-1，示例代码：

import mpmath as mp

# 设置精度，比如50位小数
mp.mp.dps = 50

# 转换为mpmath矩阵
A_mp = mp.matrix(A.tolist())
c_mp = mp.matrix(c.tolist())

# 计算伪逆并求解
A_inv_mp = mp.pinv(A_mp)
b_mp = A_inv_mp * c_mp

# 可选：转回NumPy高精度数组
b = np.array(b_mp.tolist(), dtype=np.float128)

• 正则化缓解病态问题：若A的条件数（np.linalg.cond(A)）很大，说明矩阵病态，会放大计算误差。可加入正则化项求解带约束的最小二乘解：

lambda_ = 1e-6  # 可调的正则化参数
A_T_A = np.matmul(A.T, A)
reg_matrix = lambda_ * np.eye(A.shape[1])
b_reg = np.matmul(np.matmul(np.linalg.inv(A_T_A + reg_matrix), A.T), c)

• 手动拆解方程组减少误差：观察A和c的分块结构，可拆分方程组手动求解：
  1. 前9行对应子方程：(X_matrix - rho_0)*b1 = 0、(X_matrix - rho_1)*b2 = 0、(X_matrix - rho_2)*b3 = 0
  2. 最后3行：b1 + b2 + b3 = I
     若X_matrix - rho_i奇异，先求每个子方程的通解，再代入最后一行确定通解中的参数，这种方式能避免整体矩阵运算的误差累积。

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