高斯消元程序中Numpy数组未正确更新问题排查

我正在编写名为gaussian_elim的Python函数，该函数接收n×n的Numpy数组A和n×1的Numpy数组b，对增广矩阵[A|b]执行高斯消元操作，返回形式为M=[U|c]的n×(n+1)矩阵（其中U为n×n上三角矩阵）。但测试简单2×2矩阵时，发现消元步骤未正确执行：打印显示M[i,:] - m*M[j,:]的计算结果正确，但赋值后M的第i行并未更新为该正确值，实际输出与预期不符。

函数代码 #

import numpy as np

def gaussian_elim(A,b):
    """
    A: n x n numpy array
    b: n x 1 numpy array 

    Applies Gaussian Elimination to the system Ax = b.
    Returns a matrix of the form M = [U|c], where U is upper triangular.
    """
    
    n = len(b) 
    b = b.reshape(-1, 1)        # Ensure b is a column vector of shape (n, 1)
    M = np.hstack((A,b))        #Form the n x (n+1) augmented matrix M := [A|b]
    
    #For each pivot:
    for j in range(n-1):  #j = 0,1,...,n-2
    
        #For each row under the pivot:
        for i in range(j+1,n):               #i = j + 1, j + 2,..., n-1

            if (M[j,j] == 0):
                print("Error! Zero pivot encountered!")
                return 

            #The multiplier for the the i-th row 
            m = M[i,j] / M[j,j]
 
            print("M[i,:] = ", M[i,:])
            print("M[j,:] = ", M[j,:])
            print("m = ", m)
            print("M[i,:] - m*M[j,:] = ", M[i,:] - m*M[j,:])

            #Eliminate entry M[i,j] (the first nonzero entry of the i-th row)
            M[i,:] = M[i,:] - m*M[j,:]
            
            print("M[i,:] = ", M[i,:])   #Make sure that i-th row of M is correct (it's not!)
       
    return M

测试代码 #

A = np.array([[3,-2],[1,5]])
b = np.array([1,1])
gaussian_elim(A,b)

测试输出 #

M[i,:] =  [1 5 1]
M[j,:] =  [ 3 -2  1]
m =  0.3333333333333333
M[i,:] - m*M[j,:] =  [0.         5.66666667 0.66666667]    <-- 这个结果是正确的！
M[i,:] =  [0 5 0]                                          <-- 为什么这行和上面的结果不一样？

array([[ 3, -2,  1],
       [ 0,  5,  0]])

预期输出是array([[ 3, -2, 1],[0. 5.66666667 0.66666667]])，请问为什么第二行没有正确更新？

问题根源 #

你的数组A和b都是整数类型（numpy.int64），通过np.hstack构造的增广矩阵M也继承了整数类型。当你计算M[i,:] - m*M[j,:]时，虽然临时结果是浮点数，但赋值给整数类型的数组行时，Numpy会自动将浮点数截断为整数，导致结果被舍入，出现[0 5 0]这种错误结果。

解决办法 #

有两种简单的修复方式：

1. 创建数组时指定浮点类型：
   在定义A和b时，添加dtype=np.float64参数：

   A = np.array([[3,-2],[1,5]], dtype=np.float64)
   b = np.array([1,1], dtype=np.float64)
   

2. 构造增广矩阵时转换为浮点类型：
   在构造M时，直接转换为浮点数组，避免后续赋值时的类型截断：

   M = np.hstack((A,b)).astype(np.float64)
   

修复后的测试输出 #

应用上述任一方法后，函数会输出预期结果：

M[i,:] =  [1. 5. 1.]
M[j,:] =  [ 3. -2.  1.]
m =  0.3333333333333333
M[i,:] - m*M[j,:] =  [0.         5.66666667 0.66666667]
M[i,:] =  [0.         5.66666667 0.66666667]

array([[ 3.        , -2.        ,  1.        ],
       [ 0.        ,  5.66666667,  0.66666667]])

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