我通過以下鏈接實現了矩陣的 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)我哪里錯了？
1 回答
ibeautiful
TA貢獻1993條經驗 獲得超6個贊
看來您可能在第一個內層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 分解。希望這可以幫助。
添加回答
舉報
0/150
提交
取消