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

Question

我已经通过以下链接实现了矩阵的LU分解的标准方程/算法：(1(和(2(这完美地返回了如下所示的平方矩阵的LU分解。然而，我的问题是，它也给出了一个Divide by Zero warning。

此处编码：

import numpy as np
def 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之前求值；您也没有正确地计算求和项CCD_ 5，并且没有正确地将CCD_。在进行一些其他语法修改后，您可以如下重写您的函数：

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实现不包括旋转，只有当scipy.linalg.lu返回的置换矩阵P是单位矩阵时，这两个比较才成立，即scipy没有执行任何置换，这确实是您的矩阵A的情况。scipy算法中确定的置换矩阵是为了优化生成矩阵的条件个数，以减少舍入误差。最后，你可以简单地验证A = LU，如果因子分解做得正确，情况总是如此：

A = np.random.rand(10,10)
L, U = LUDecomposition(A)
np.allclose(A, np.dot(L, U))
>>>  True

然而，就数值效率和精度而言，我不建议您使用自己的函数来计算LU分解。希望这能有所帮助。