人和计算机用都高效且简洁的LU分解

矩阵的LU三角分解

算法

1. 初始化矩阵 $L=\left[\begin{array}{lll}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$, $U=\left[\begin{array}{lll}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$, $R=\left[\begin{array}{lll}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]$
2. 更新 $U$, 其中 $U[i]=A[i]-L[i]\cdot U$
3. 更新 $R$, 其中 $R[i]=U[i],R[i][i]=0$
4. 更新 $L$, 其中 $L[:,i]=\frac{A[:,i]-L\cdot R[:,i]}{U[i][i]}$
5. 重复步骤2到步骤4

def LU(A):
    n = A.shape[0]
    L = np.diag(np.ones(n))
    U = np.zeros(A.shape)
    # 引入R来简化更新L的逻辑, A = LU = L(D+R)
    R = np.zeros(A.shape)
    for i in range(n):
        # 更新第i行U
        dU = np.dot(L[i], U)
        U[i] = A[i] - dU
        R[i] = U[i]
        R[i][i] = 0
        # 更新第i列L
        dL = np.dot(L, R[:, i]).reshape(-1, 1)
        L[:, [i]] = (A[:, [i]] - dL) / U[i][i]
    return L, U