不选主元的矩阵三角分解法

1.矩阵三角分解原理

设$A$为非奇异矩阵，且有分解式$$A=LU,$$
$$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1n}\\ a_{21}& a_{21}& \dots & a_{2n}\\ ⋮& ⋮& & ⋮\\ a_{n1}& a_{n1}& \dots & a_{nn}\end{array}\right]$$
其中$L$为单位下三角矩阵,$U$为上三角矩阵,即
$$A=\left[\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \dots & 1\end{array}\right]\left[\begin{array}{cccc}{u}_{11}& u_{12}& \dots & u_{1n}\\ & u_{22}& \dots & u_{2n}\\ & & \ddots & ⋮\\ & & & u_{nn}\end{array}\right]$$
下面说明$L,U$的元素可以由$n$步直接定出,其中第$r$步定出$U$的第$r$行和$L$的第$r$列.由上述分解式有$$a_{1i}=u_{1i},i=1,2,\dots ,n,$$
得到$U$的第$1$行元素;
$$a_{i1}=l_{i1}{u}_{11},l_{i1}=a_{i1}/u_{11},i=2,3,\dots ,n,$$
得到$L$的第$1$列元素.
设已经定出$U$的第1行到第$r-1$列元素.有上述分解式，利用矩阵乘法(注意当$r<k,l_{rk}=0$),有$$a_{ri}=\sum _{k=1}^n{l}_{rk}{u}_{ki}=\sum _{k=1}^{r-1}{l}_{rk}{u}_{ki}+u_{ri},$$
故$$u_{ri}=a_{ri}-\sum _{k=1}^{r-1}{l}_{rk}{u}_{ki},i=r,r+1,\dots ,n.$$
又由该分解式有$$a_{ir}=\sum _{k=1}^n{l}_{ik}{u}_{kr}=\sum _{k=1}^{r-1}{l}_{ik}{u}_{kr}+l_{ir}{u}_{rr}.$$
总结上述讨论，得到直接三角分解法解$Ax=b$(要求$A$的所有顺序主子式都不为零)的计算公式.

1. $u_{1i}=a_{1i}(i=1,2,\dots ,n),l_{i1}=a_{i1}/u_{11},i=2,3,\dots ,n.$
   计算$U$的第$r$行,$L$的第$r$列元素$(r=2,3,\dots ,n):$
2. $u_{ri}=a_{ri}-\sum _{k=1}^{r-1}{l}_{rk}{u}_{ki},i=r,r+1,\dots ,n;$
3. $l_{ir}=(a_{ir}-\sum _{k=1}^{r-1}{l}_{ik}{u}_{kr})/u_{rr},i=r+1,\dots ,n,\text{且}r\ne n.$
   求解$Ly=b,Ux=y$的计算公式:
4. $\{\begin{array}{l}{y}_1=b_1,\\ y_i=b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k,i=2,3,\dots ,n;\end{array}$
5. $\{\begin{array}{l}{x}_n=y_n/u_{nn},\\ x_i=(y_i-\sum _{k=i+1}^n{u}_{ik}{x}_k)/u_{ii},i=n-1,n-2,\dots ,1.\end{array}$

2.Python实现不选主元矩阵三角分解

# 自己原创
def triangle_decomposition(coefficient_matrix: np.ndarray):
    """
    实现方程Ax=b系数矩阵A的Doolittle decomposition
    :param coefficient_matrix: 初始系数矩阵A
    :return: 单位下三角矩阵L,上三角矩阵U
    """
    # first step: evaluate the decomposition condition
    rows, columns = coefficient_matrix.shape
    if rows == columns:  # judge if it is a square matrix
        for k in range(rows):  # 判断各阶顺序主子式是否为0
            if det(coefficient_matrix[:k + 1, :k + 1]) == 0:
                raise Exception("cannot decompose ")
        else:
            # directive triangle  decomposition
            # 初始化
            square_ones_matrix = np.ones((rows, columns))
            lower_triangle_matrix = np.tril(square_ones_matrix)
            upper_triangle_matrix = np.triu(square_ones_matrix)
            # 计算第1行的U(1,i),第1列的L(i,1)
            upper_triangle_matrix[0, :] = coefficient_matrix[0, :]
            lower_triangle_matrix[1:, 0] = coefficient_matrix[1:, 0] / upper_triangle_matrix[0, 0]
            # 计算第row行的U(row,i),第row列的L(i,row)
            for row in range(1, rows - 1):

                for i in range(row, rows):
                    # 定义一个累积和变量
                    row_i_prod = 0
                    for k in range(row):
                        row_i_prod += lower_triangle_matrix[row, k] * upper_triangle_matrix[k, i]
                    upper_triangle_matrix[row, i] = coefficient_matrix[row, i] - row_i_prod

                for i in range(row + 1, rows):
                    # 定义一个累积和变量
                    i_row_prod = 0
                    for k in range(row):
                        i_row_prod += lower_triangle_matrix[i, k] * upper_triangle_matrix[k, row]
                    lower_triangle_matrix[i, row] = (coefficient_matrix[i, row] - i_row_prod) / upper_triangle_matrix[
                        row, row]

            # 定义一个累积和变量,单独计算n行n列的U(n,n)
            row_i_prod = 0
            for k in range(rows - 1):
                row_i_prod += lower_triangle_matrix[rows - 1, k] * upper_triangle_matrix[k, rows - 1]
            upper_triangle_matrix[rows - 1, rows - 1] = coefficient_matrix[rows - 1, rows - 1] - row_i_prod

            return lower_triangle_matrix, upper_triangle_matrix
    else:
        raise Exception("ERROR:please pass a square matrix.")
# 自己原创
def directive_triangle_decomposition(coefficient_matrix: np.ndarray, right_hand_side_vector: np.ndarray,
                                     decomposition_function=triangle_decomposition):
    """
    实现方程Ax=b系数矩阵A directive triangle decomposition
    :param decomposition_function: 三角分解函数
    :param coefficient_matrix: 初始系数矩阵A
    :param right_hand_side_vector: 初始常数列向量b
    :return: 解向量x
    """
    l, u = decomposition_function(coefficient_matrix)
    rows, columns = coefficient_matrix.shape
    # Ly=b，此处不指定双精度浮点数据类型，则可能会损失精度，导致结果偏差太大
    y = np.ones_like(right_hand_side_vector, dtype=np.float64)
    y[0, 0] = right_hand_side_vector[0, 0]
    for row in range(1, rows):
        # 定义一个累积和变量
        i_k_prod = 0
        for k in range(row):
            i_k_prod += l[row, k] * y[k, 0]
        y[row, 0] = right_hand_side_vector[row, 0] - i_k_prod

    # Ux=y，此处不指定双精度浮点数据类型，则可能会损失精度，导致结果偏差太大
    x = np.ones_like(right_hand_side_vector, dtype=np.float64)
    x[rows - 1] = y[rows - 1] / u[rows - 1, rows - 1]
    for row in range(rows - 2, 0, -1):
        i_k_prod = 0
        for k in range(row + 1, rows):
            i_k_prod += u[row, k] * x[k, 0]
        x[row, 0] = (y[row, 0] - i_k_prod) / u[row, row]

    return x

3.解方程

用直接三角分解法解$$\left[\begin{array}{ccc}1& 2& 3\\ 2& 5& 2\\ 3& 1& 5\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}14\\ 18\\ 20\end{array}\right]$$

4.测试

if __name__ == '__main__':
# 直接三角分解法-不选主元三角分解法测试成功，来源详见李庆扬数值分析第5版P153，e.g.5
    square_matrix = np.array([[1, 2, 3], [2, 5, 2], [3, 1, 5]])
    column_vector = np.array([14, 18, 20]).reshape((3, 1))
    print(directive_triangle_decomposition(square_matrix, column_vector))

5.运行结果