LU分解中除以零警告-Doolittle算法工作

Question

我通过以下链接实现了矩阵的 LU 分解的标准方程/算法：（1）和（2）这完美地返回了如下方矩阵的 LU 分解。但是，我的问题是-它也给出了Divide by Zero warning.代码在这里：import numpy as npdef LUDecomposition (A):    L = np.zeros(np.shape(A),np.float64)    U = np.zeros(np.shape(A),np.float64)    acc = 0    L[0,0]=1    for i in np.arange(len(A)):        for k in range(i,len(A)):            for j in range(0,i):                acc += L[i,j]*U[j,k]            U[i,k] = A[i,k]-acc            for m in range(k+1,len(A)):                if m==k:                    L[m,k]=1                else:                    L[m,k] = (A[m,k]-acc)/U[k,k]            acc=0    return (L,U)A = np.array([[-4, -1, -2],              [-4, 12,  3],              [-4, -2, 18]])L, U = LUDecomposition (A)我哪里错了？

Answer 1

看来您可能在第一个内层for循环方面犯了一些缩进错误：U必须在之前进行评估L；您也没有正确计算求和项acc，也没有正确地将对角项设置L为 1。在进行一些其他语法修改后，您可以按如下方式重写您的函数：

def LUDecomposition(A):

    n = A.shape[0]

    L = np.zeros((n,n), np.float64)

    U = np.zeros((n,n), np.float64)

    for i in range(n):

        # U

        for k in range(i,n):

            s1 = 0  # summation of L(i, j)*U(j, k) 

            for j in range(i):

                s1 += L[i,j]*U[j,k]

            U[i,k] = A[i,k] - s1

        # L

        for k in range(i,n):

            if i==k:

                # diagonal terms of L 

                L[i,i] = 1

            else:

                s2 = 0 # summation of L(k, j)*U(j, i) 

                for j in range(i):

                    s2 += L[k,j]*U[j,i]

                L[k,i] = (A[k,i] - s2)/U[i,i]

    return L, U

与scipy.linalg.lu作为可靠参考A相比，这一次给出了矩阵的正确输出：

import numpy as np

from scipy.linalg import lu

A = np.array([[-4, -1, -2],

              [-4, 12,  3],

              [-4, -2, 18]])

L, U = LUDecomposition(A)

P, L_sp, U_sp = lu(A, permute_l=False)

P

>>> [[1., 0., 0.],

     [0., 1., 0.],

     [0., 0., 1.]])

L

>>> [[ 1.          0.          0.        ]

     [ 1.          1.          0.        ]

     [ 1.         -0.07692308  1.        ]]

np.allclose(L_sp, L))

>>>  True

U

>>> [[-4.         -1.         -2.        ]

     [ 0.         13.          5.        ]

     [ 0.          0.         20.38461538]]

np.allclose(U_sp, U))

>>>  True

注意：与 scipy lapack getrf 算法不同，此 Doolittle 实现不包括旋转，这两个比较只有在P返回的置换矩阵scipy.linalg.lu是单位矩阵时才为真，即scipy 没有执行任何置换，这确实是您的矩阵的情况A. 在 scipy 算法中确定的置换矩阵旨在优化结果矩阵的条件数，以减少舍入误差。最后，您可以简单地验证一下A = LU，如果分解正确，情况总是如此：

A = np.random.rand(10,10)

L, U = LUDecomposition(A)

np.allclose(A, np.dot(L, U))

>>>  True

尽管如此，就数值效率和准确性而言，我不建议您使用自己的函数来计算 LU 分解。希望这可以帮助。