实现方阵LU分解算法时NumPy二维数组切片赋值未按预期更新

问题根因 #

你遇到的切片赋值不生效问题由两个原因共同导致：

• 输入数组的类型限制：你定义的A默认是整数（int）类型，乘子A[i+1:,i] / A[i,i]的计算结果是浮点数，直接赋值给int类型的数组切片时，NumPy会自动把浮点数截断为整数，最终赋值结果和预期不符。
• 循环范围错误：外层循环遍历了整个矩阵的所有行索引range(A.shape[0])，当遍历到最后一个行索引时，i+1:的切片是空数组，没有实际运算意义，后续给L矩阵赋值时还会触发索引溢出报错。

修复方案 #

对代码做两处调整即可：

1. 函数入口处先对输入矩阵做拷贝，同时转换为浮点类型，既避免修改原输入数组，也解决类型截断问题
2. 把外层循环范围修改为range(A.shape[0]-1)，只遍历到倒数第二行即可，最后一行不需要执行消元操作

修复后代码 #

import numpy as np

def LUGAUSS(A):
    if A.shape[0] != A.shape[1]:
        return "Invalid Matrix. A must be a square marix."
    # 转换为浮点类型的拷贝，不修改原数组同时避免整数截断
    A = A.astype(np.float64).copy()
    multipliers = dict()
    # 循环到倒数第二行即可
    for i in range(A.shape[0]-1):
        print('i',i)
        if A[i,i] == 0:
            return "Pivot is zero"
        else:
            multipliers[(i+1,i)] = A[i+1:,i] / A[i,i]
            A[i+1:,i] =  multipliers[(i+1,i)]
            A[i+1:,i+1:] = A[i+1:,i+1:] -  A[i+1:,i].reshape(-1,1) * A[i,i+1:].reshape(1,-1)
           
    L = np.eye(A.shape[0])
    for x in range(L.shape[1]-1):
        L[x+1:,x] = multipliers[(x+1,x)]

    U = A.copy()
    for x in range(U.shape[1]):
        U[x+1:,x] = 0 

    return (L,U,multipliers)

# 测试调用
A = np.array([[1, -3, 5, 2], [1, 0, 1, -1], [6, 1, -9, 2], [1, 0, -6, 3]])
L,U,m = LUGAUSS(A)
print("L矩阵：\n", L)
print("U矩阵：\n", U)

运行修复后代码即可得到正确的LU分解结果，切片赋值也会按照预期更新为浮点数值。

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