我需要实现一个LU分解,然后将其与numpy中的np.linalg.solve
函数进行比较。
代码中的函数(见下文)运行没有任何问题,但是当我使用它来解决矩阵时,我不断收到错误:
IndexError: list index out of range
在线上:
L[i][j] = (A2[i][j] - s2) / U[j][j]
这是整个代码:
def matrixMul(A, B):
TB = zip(*B)
return [[sum(ea*eb for ea,eb in zip(a,b)) for b in TB] for a in A]
def pivotize(m):
#Creates the pivoting matrix for m.
n = len(m)
ID = [[float(i == j) for i in range(n)] for j in range(n)]
for j in range(n):
row = max(range(j, n), key=lambda i: abs(m[i][j]))
if j != row:
ID[j], ID[row] = ID[row], ID[j]
return ID
def lu(A):
#Decomposes a nxn matrix A by PA=LU and returns L, U and P.
n = len(A)
L = [[0.0] * n for i in range(n)]
U = [[0.0] * n for i in range(n)]
P = pivotize(A)
A2 = matrixMul(P, A)
for j in range(n):
L[j][j] = 1.0
for i in range(j+1):
s1 = sum(U[k][j] * L[i][k] for k in range(i))
U[i][j] = A2[i][j] - s1
for i in range(j, n):
s2 = sum(U[k][j] * L[i][k] for k in range(j))
L[i][j] = (A2[i][j] - s2) / U[j][j]
return (L)
A = np.array([[1,1,3],[5,3,1],[2,3,1]])
b = np.array([2,3,-1])
print('LU factorization: ', lu(A))
A = np.array([[1,1,3],[5,3,1],[2,3,1]])
b = np.array([2,3,-1])
print('Internal solver : ', np.linalg.solve(A,b))
有什么想法吗?谢谢!
您的matrixMul
方法不正确。试试这个:
matrixMul([[1, 0], [0, 1]], [[1, 0], [0, 1]])
这是将两个单位矩阵相乘,应该返回[[1, 0], [0, 1]]
。当我运行它时,它会返回[[1, 0], []]
.这意味着A2
,在您的lu
内只有一行,其余的都是空的。因此,i == 1
和j == 0
时的索引错误。
失败的原因是TB
是一个zip
对象。这些只能迭代一次,消耗迭代器。我实际上认为您根本不需要TB
对象,只需迭代B
元素即可。
def matrixMul(A, B):
return [[sum(ea*eb for ea,eb in zip(a,b)) for b in B] for a in A]
这将返回所需的输出:
>>> matrixMul([[1, 0], [0, 1]], [[1, 0], [0, 1]])
>>> [[1, 0], [0, 1]]
编辑:
顺便说一句,您的解决方案还有其他问题。您的版本和 NumPy 的版本仍然不匹配。但是这里的解决方案将修复您的索引错误。