Pythonを使用したLU分解3

Question

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)) 

です。ありがとう！

Answer 1

matrixMulの方法が正しくありません。これは、2つの単位行列を乗算され

matrixMul([[1, 0], [0, 1]], [[1, 0], [0, 1]]) 

、および[[1, 0], [0, 1]]を返す必要があります：これを試してみてください。私はそれを実行すると、[[1, 0], []]を返します。つまり、luの内部にはA2が1行しかなく、残りは空です。したがって、i == 1とj == 0の場合のインデックスエラー。

この理由は、TBがzipオブジェクトであるためです。それらは一度だけ繰り返すことができ、イテレータを消費します。実際にはTBオブジェクトが必要だとは思っていませんが、Bの要素を反復処理するだけです。

>>> matrixMul([[1, 0], [0, 1]], [[1, 0], [0, 1]]) 
>>> [[1, 0], [0, 1]] 

EDIT：

def matrixMul(A, B): 
    return [[sum(ea*eb for ea,eb in zip(a,b)) for b in B] for a in A] 

これは、所望の出力を返しところで

を、あなたのソリューションを持つ他の問題が残っています。あなたとNumPyのバージョンはまだ一致しません。しかし、ここの解決策はインデックスエラーを修正します。