Crout矩阵分解的Python实现

Question

我想写一个适用于3x3矩阵的Crout矩阵分解的Python实现。我只允许使用numpy包。这里是我的尝试至今：

import numpy as np 

def crout(A: np.ndarray): 

    L = np.zeros((3, 3)) 
    U = np.zeros((3, 3)) 

    for k in range(0, 3): 
     U[k, k] = 1 

     for j in range(k, 3): 
      sum0 = sum(L[k, s] * U[s, j] for s in range(1, k-1)) 
      L[k, j] = A[k, j] - sum0 

     for j in range(k, 3): 
      sum1 = sum(L[k, s] * U[s, j] for s in range(1, k-1)) 
      U[k, j] = (A[k, j] - sum1)/L[k, k] 

    print(" L =", '\n', L, '\n', " U =", '\n', U) 
    return L, U 

A = np.array([[60.0, 30.0, 20.0], [30.0, 20.0, 15.0], [20.0, 15.0, 12.0]]) 
crout(A) 

我使用的矩阵A我的函数，它应该产生的尝试：

而是给

所以，显然输出是错误的。我的索引是错误还是存在另一个问题？

Answer 1

看起来你翻转的j和K在第二个for循环

def crout(A): 

    L = np.zeros((3, 3)) 
    U = np.zeros((3, 3)) 

    for k in range(0, 3): 
     U[k, k] = 1 

     for j in range(k, 3): 
      sum0 = sum(L[k, s] * U[s, j] for s in range(1, k-1)) 
      #reversed 
      L[j, k] = A[k, k] - sum0 

     for j in range(k, 3): 
      sum1 = sum(L[k, s] * U[s, j] for s in range(1, k-1)) 
      U[k, j] = (A[k, j] - sum1)/L[k, k] 

    print L 
    print U 
    return L, U