矩阵的LU分解

矩阵LU分解的python实现

import numpy as np


def lu(A):
    # 方阵的维度
    n = A.shape[0]
    # 求L,U
    U = A
    L = np.eye(n)
    # 利用高斯消元的方法一次性求出L,U
    for i in range(n - 1):
        L[i + 1:, i] = U[i + 1:, i] / U[i, i]
        U[i + 1:, i + 1:] -= L[i + 1:, i].reshape(-1, 1).dot(U[i, i + 1:].reshape(1, -1))
        U[i + 1:, i] = 0
    return L, U


if __name__ == '__main__':
    # example
    A = np.array(
        [[2, 1, 4, 5, 2],
         [1, 2, -3, 0, 1],
         [2, 3, -7, 4, 5],
         [1, 2, 2, 5, 7],
         [9, 3, 5, 1, -3]
         ], dtype=np.float64
    )
    L, U = lu(A)
    print(L)
    print(U)
    print(L.dot(U))

测试实例：

# L
[[ 1.          0.          0.          0.          0.        ]
 [ 0.5         1.          0.          0.          0.        ]
 [ 1.          1.33333333  1.          0.          0.        ]
 [ 0.5         1.         -1.15384615  1.          0.        ]
 [ 4.5        -1.          4.15384615 -4.38        1.        ]]
# U
[[ 2.          1.          4.          5.          2.        ]
 [ 0.          1.5        -5.         -2.5         0.        ]
 [ 0.          0.         -4.33333333  2.33333333  3.        ]
 [ 0.          0.          0.          7.69230769  9.46153846]
 [ 0.          0.          0.          0.         16.98      ]]
# A=LU
[[ 2.  1.  4.  5.  2.]
 [ 1.  2. -3.  0.  1.]
 [ 2.  3. -7.  4.  5.]
 [ 1.  2.  2.  5.  7.]
 [ 9.  3.  5.  1. -3.]]