Python行约简求矩阵行列式与np.linalg.det结果误差问题咨询

你遇到的误差问题来自三个核心原因：

• 整型数组溢出
  你输入的整型矩阵默认使用numpy固定精度整型存储（通常为int32或int64），25阶、元素取值范围在±1000的矩阵，行列式理论值最高可达(1000)^25=1e75，远超过int64的最大值（约9.2e18）。计算过程中整数溢出后会被截断为错误的溢出值，这是你测试25阶整型矩阵时结果差了几十个数量级的核心原因。
  修复方案是计算前先将输入矩阵转换为双精度浮点型：

A = np.matrix.copy(A).astype(np.float64)

• 未使用部分主元消去，数值稳定性差
  你当前的find_non_zero函数仅选择当前列第一个非零元素作为主元，没有选择当前列中对角线及以下绝对值最大的元素作为主元。如果选中的主元绝对值极小，消元时计算的factor = A[row, col] / A[col, col]会非常大，后续减法操作会把微小的舍入误差放大数个数量级。矩阵阶数越高，这种误差的累积效应越明显，就是70×70以上矩阵误差暴涨的主要原因。
  numpy的np.linalg.det底层调用LAPACK的带部分主元的LU分解实现，数值稳定性远优于你当前的朴素高斯消去实现，因此误差更小。

• 浮点运算的累积误差
  行约简过程是O(n³)复杂度的浮点运算，每一步的除法、乘法、减法都会引入微小的舍入误差，矩阵阶数越高，运算次数越多，误差累积量越大。最后对n个对角线元素求乘积的操作，也会进一步放大总误差。

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