Python使用高斯-若尔当消元法求矩阵逆的代码问题求解

问题原因分析 #

• 错误1：行交换逻辑冗余错误
  原代码行交换部分先循环交换b矩阵k到n列的元素，随后又执行了整行交换b[k], b[i] = b[i], b[k]，相当于对b矩阵做了两次交换，导致行顺序错乱；且行交换本质是整行操作，不需要限制列范围。
• 错误2：b矩阵变换的列遍历范围错误
  对a矩阵做缩放、消元时仅遍历k到n列没问题（k左侧元素已经完成化简无需修改），但b矩阵需要同步对所有列执行变换，原代码仅修改k到n列，导致b矩阵前几列未参与计算，结果错误。
• 可选优化：原代码中判断a[i][k] == 0的逻辑不严谨，浮点计算场景下应该用绝对值小于阈值判断，避免精度问题导致逻辑异常。

修正后代码 #

M = [[0, 2, 1], [4, 0, 1], [-1, 2, 0]]
P = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
n = len(P)
def inverse(a, b):
    for k in range(n):
        # 选主元+行交换
        if abs(a[k][k]) < 1.0e-12:
            for i in range(k+1, n):
                if abs(a[i][k]) > abs(a[k][k]):
                    # 整行交换，不需要限制列范围
                    a[k], a[i] = a[i], a[k]
                    b[k], b[i] = b[i], b[k]
                    break
        pivot = a[k][k]
        # 主元行缩放：a仅修改k到n列，b修改所有列
        for j in range(k, n):
            a[k][j] /= pivot
        for j in range(n):
            b[k][j] /= pivot
        # 消去其他行的k列元素
        for i in range(n):
            if i == k or abs(a[i][k]) < 1.0e-12:
                continue
            factor = a[i][k]
            for j in range(k, n):
                a[i][j] -= factor * a[k][j]
            for j in range(n):
                b[i][j] -= factor * b[k][j]
    return a, b

a, b = inverse(M, P)
print("化简后的原矩阵：")
for row in a:
    print(row)
print("逆矩阵：")
for row in b:
    print(row)

输出结果 #

化简后的原矩阵：
[1.0, 0.0, 0.0]
[0.0, 1.0, 0.0]
[0.0, 0.0, 1.0]
逆矩阵：
[-0.3333333333333333, 0.3333333333333333, 0.3333333333333333]
[-0.16666666666666666, 0.16666666666666666, 0.6666666666666666]
[1.3333333333333333, -0.3333333333333333, -1.3333333333333333]

和参考逆矩阵[[-1/3, 1/3, 1/3], [-1/6, 1/6, 2/3], [4/3, -1/3, -4/3]]完全一致。

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