求排查：输入为RREF的矩阵时零空间基计算错误问题

问题现象 #

给定的RREF计算函数，在输入已处于行最简形的矩阵时，无法生成正确的零空间基：

• 输入矩阵：

  numpy.array([[1,0,3,0],[0,0,2,1]], dtype=float)
  

• 错误输出：生成3个零空间基向量（实际应为2个），且向量值不符合线性代数规则
• 正确零空间基应为：[0, 1, 0, 0] 和 [1.5, 0, -0.5, 1]
• 输入非RREF矩阵时，输出结果符合预期。

错误原因分析 #

问题出在RREF转换的循环逻辑中，主元行的递增规则错误：
原代码中，无论当前列是否找到主元，每次循环都会同时递增i（行索引）和j（列索引）。这会导致：

1. 当遇到全0列时，i依然递增，直接跳过后续行的主元处理
2. 对于输入已为RREF的矩阵，第二行的主元列（列2）未被识别，l列表（记录主元列）仅包含列0，错误地将列2判定为自由变量列
3. 最终计算出的零空间基维度错误，向量值也不符合线性代数规则。

具体来说，输入矩阵的第二行主元在列2，但循环执行到i=1, j=1时，列1全0，j递增到2，但i同时递增到2，此时i >= m（矩阵行数为2），循环直接退出，完全没处理列2的主元。

修复方案 #

修改循环内的i递增逻辑：仅在找到主元并完成当前行的归一化、消元操作后，才递增i；若当前列全0，仅递增j，保持i不变。

修复后的代码：

import numpy as np

def r_r_e_f(matrix):
    a = matrix.copy()
    m, n = np.shape(a)
    i = 0
    j = 0
    rank = 0
    l = []
    while i <= (m-1) and j <= (n-1):
        p = np.max(np.abs(a[i:m, j]))
        max_p = np.argmax(np.abs(a[i:m, j]))
        k = max_p + i
        if p == 0 or p <= 1e-18:
            # 当前列全0，仅移动列索引
            j += 1
        else:
            rank += 1
            # 交换当前行与主元行
            a[[i, k]] = a[[k, i]]
            # 主元归一化
            a[i, :] = a[i, :] / a[i, j]
            l.append(j)
            # 消去其他行的主元
            for x in range(m):
                if x != i:
                    a[x, :] = a[x, :] - a[x, j] * a[i, :]
            # 仅在处理完主元行后移动行索引
            i += 1
            j += 1
    # 确定自由变量列
    ncols = [col for col in range(n) if col not in l]
    nbasis = []
    for col in ncols:
        basis = np.zeros(n)
        basis[col] = 1
        for row in range(len(l)):
            pcol = l[row]
            basis[pcol] = -a[row, col]
        nbasis.append(basis)
    return a, nbasis

# 测试已为RREF的矩阵
matrix = np.array([[1,0,3,0],[0,0,2,1]], dtype=float)
rref, null_space = r_r_e_f(matrix)
print("原矩阵:")
print(matrix)
print("RREF结果: ")
print(rref)
print("零空间基: ")
if not null_space:
    print("零空间为零向量")
else:
    for bv in null_space:
        print(bv.reshape(-1, 1))

验证结果 #

修复后运行测试代码，输出的零空间基为：

[[0.]
 [1.]
 [0.]
 [0.]]
[[ 1.5]
 [ 0. ]
 [-0.5]
 [ 1. ]]

与预期结果完全一致，且输入非RREF矩阵时的输出依然正确。

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