线性代数行变换代码异常求助：行约简结果与预期不符

1. 核心问题1：整数数据类型导致浮点数运算截断 #

你定义矩阵A和b时指定了dtype=int，但行变换中涉及非整数运算（比如最后一步乘以-1/17），NumPy整数数组会自动截断小数部分，直接导致结果失真。例如最后一步对行3的缩放操作，整数类型会将-1/17的结果强制转为0，这是实际输出最后一行最后一个元素错误的直接原因。

2. 核心问题2：行变换步骤逻辑混乱 #

你的行变换顺序不符合行阶梯形矩阵（REF）的构建规则（从左到右、从上到下依次确定主元，消去主元下方的所有非零元素），存在大量无效或错误操作：

• 交换行1和行3后，未优先处理第二列的主元消元，反而先执行行2与行3的无意义相加，导致主元列下方未被正确消去；
• 重复执行AddRows(A_ref, 2, 3, 1)这类操作，打乱了已有的消元结果；
• 未按顺序对每一列的主元进行标准化（将主元转为1），后续消元失去基准。

修正后的代码及正确行变换步骤 #

首先将数据类型改为float64支持浮点数运算，再按照标准行约简步骤执行：

import numpy as np

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

def MultiplyRow(M, row_num, row_num_multiple):
    M_new = M.copy()
    M_new[row_num] = M_new[row_num] * row_num_multiple
    return M_new

def AddRows(M, row_num_1, row_num_2, row_num_1_multiple):
    M_new = M.copy()   
    M_new[row_num_2] = row_num_1_multiple * M_new[row_num_1] + M_new[row_num_2]
    return M_new

def SwapRows(M, row_num_1, row_num_2):
    M_new = M.copy()
    M_new[[row_num_1, row_num_2]] = M_new[[row_num_2, row_num_1]]
    return M_new

def augmented_to_ref(A, b):    
    A_system = np.hstack((A, b.reshape((4, 1))))
    print("初始增广矩阵：")
    print(A_system)

    # 步骤1：交换行0和行1，使第一列主元为1
    A_ref = SwapRows(A_system, 0, 1)
    
    # 步骤2：消去第一列主元下方的所有元素
    A_ref = AddRows(A_ref, 0, 1, -2)  # 行1 = -2*行0 + 行1
    A_ref = AddRows(A_ref, 0, 2, 1)   # 行2 = 1*行0 + 行2
    A_ref = AddRows(A_ref, 0, 3, -1)  # 行3 = -1*行0 + 行3
    
    # 步骤3：将第二列主元标准化为1
    A_ref = MultiplyRow(A_ref, 1, -1/5)
    
    # 步骤4：消去第二列主元下方的所有元素
    A_ref = AddRows(A_ref, 1, 2, -4)  # 行2 = -4*行1 + 行2
    A_ref = AddRows(A_ref, 1, 3, 1)   # 行3 = 1*行1 + 行3
    
    # 步骤5：将第三列主元标准化为1
    A_ref = MultiplyRow(A_ref, 2, 5/17)
    
    # 步骤6：消去第三列主元下方的元素
    A_ref = AddRows(A_ref, 2, 3, -17/5)
    
    # 步骤7：将第四列主元标准化为1
    A_ref = MultiplyRow(A_ref, 3, -1/17)
    
    return A_ref

A_ref = augmented_to_ref(A, b)
print("\n行阶梯形矩阵（REF）：")
print(np.round(A_ref, decimals=0))

验证结果 #

修正后运行代码，取整后的输出与预期完全一致：

[[ 1.  2. -1. -1.  3.]
 [ 0.  1.  4.  3. 22.]
 [ 0.  0.  1.  3.  7.]
 [ 0.  0.  0.  1.  1.]]

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