平方根法求解对称正定方程组

1.对称正定矩阵三角分解或Cholesky分解定理

如果$A$为$n$阶对称正定矩阵，则存在一个实的非奇异下三角矩阵$L$使$A=L{L}^T$，当限定$L$的对角元素为正时，这种分解是唯一的.
下面用直接分解方法来确定计算$L$元素的递推公式.因为
$$A=\left[\begin{array}{cccc}{l}_{11}& & & \\ l_{21}& l_{22}& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \dots & l_{nn}\end{array}\right]\left[\begin{array}{cccc}{l}_{11}& l_{21}& \dots & l_{n1}\\ & l_{22}& \dots & l_{n2}\\ & & \ddots & ⋮\\ & & & l_{nn}\end{array}\right]$$
其中$l_{ii}>0(i=1,2,\dots ,n).$由矩阵乘法及$l_{jk}=0(当j<k时),得$
$$a_{ij}=\sum _{k=1}^n{l}_{ik}{l}_{jk}=\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk}+l_{jj}{l}_{ij},$$
于是得到解对称正定方程组$Ax=b$的平方根计算公式:
对于$j=1,2,\dots ,n$
$(1)l_{jj}=(a_{jj}-\sum _{k=1}^{j-1}{l}_{jk}^2{)}^{\frac{1}{2}};$
$(2)l_{ij}=(a_{ij}-\sum _{k=1}^{j-1}{l}_{ik}{l}_{jk})/l_{jj},i=j+1,\dots ,n.$
求解$Ax=b,$即求解两个三角形方程组:
1.$Ly=b,求y;$ 2.$L^{T}x=y,求x.$
$(3)y_i=(b_i-\sum _{k=1}^{i-1}{l}_{ik}{y}_k)/l_{ii},i=1,2,\dots ,n.$
$(4)x_i=(b_i-\sum _{k=i+1}^n{l}_{ki}{x}_k)/l_{ii},i=n,n-1,\dots ,1.$

2.Python实现平方根法

# 自己原创平方根法
def cholesky_decomposition(sym_pos_define_matrix: np.ndarray, right_hand_side_vector: np.ndarray):
    rows, columns = sym_pos_define_matrix.shape
    # judge if it is a square matrix
    if rows == columns:
        # 判断是否对称
        if (sym_pos_define_matrix.T == sym_pos_define_matrix).all():
            for k in range(rows):  # 判断各阶顺序主子式是否为0,即判定是否正定
                if det(sym_pos_define_matrix[:k + 1, :k + 1]) <= 0:
                    raise Exception("cannot decompose")
            else:
                # 初始化单位下三角矩阵L
                square_ones_matrix = np.ones((rows, columns))
                lower_triangle_matrix = np.tril(square_ones_matrix)
                for j in range(rows):
                    prod = 0
                    for k in range(j):
                        prod += lower_triangle_matrix[j, k] * lower_triangle_matrix[j, k]
                    lower_triangle_matrix[j, j] = np.sqrt(sym_pos_define_matrix[j, j] - prod)

                    for i in range(j + 1, rows):
                        prod = 0
                        for k in range(j):
                            prod += lower_triangle_matrix[i, k] * lower_triangle_matrix[j, k]
                        lower_triangle_matrix[i, j] = (sym_pos_define_matrix[i, j] - prod) / lower_triangle_matrix[j, j]

        else:
            raise Exception("error,the input matrix must be a symmetric-matrix")
    else:
        raise Exception("ERROR:please pass a square matrix.")
    return inv(lower_triangle_matrix.transpose()) @ inv(lower_triangle_matrix) @ right_hand_side_vector

3.平方根解对称正定矩阵方程组

$$\left[\begin{array}{cccccccc}4& 2& -4& 0& 2& 4& 0& 0\\ 2& 2& -1& -2& 1& 3& 2& 0\\ -4& -1& 14& 1& -8& -3& 5& 6\\ 0& -2& 1& 6& -1& -4& -3& 3\\ 2& 1& -8& -1& 22& 4& -10& -3\\ 4& 3& -3& -4& 4& 11& 1& -4\\ 0& 2& 5& -3& -10& 1& 14& 2\\ 0& 0& 6& 3& -3& -4& 2& 19\end{array}\right]\left[\begin{array}{c}x1\\ x2\\ x3\\ x4\\ x5\\ x6\\ x7\\ x8\end{array}\right]=\left[\begin{array}{c}0\\ -6\\ 20\\ 23\\ 9\\ -22\\ -15\\ 45\end{array}\right]$$

4.测试

if __name__ == '__main__':
    symmetric_positive_define_matrix2 = np.array([[4, 2, -4, 0, 2, 4, 0, 0],
                                                  [2, 2, -1, -2, 1, 3, 2, 0],
                                                  [-4, -1, 14, 1, -8, -3, 5, 6],
                                                  [0, -2, 1, 6, -1, -4, -3, 3],
                                                  [2, 1, -8, -1, 22, 4, -10, -3],
                                                  [4, 3, -3, -4, 4, 11, 1, -4],
                                                  [0, 2, 5, -3, -10, 1, 14, 2],
                                                  [0, 0, 6, 3, -3, -4, 2, 19]],dtype=np.float64)
    column_vector = np.array([0, -6, 20, 23, 9, -22, -15, 45],dtype=np.float64).reshape((8, 1))
    print(cholesky_decomposition(symmetric_positive_define_matrix2, column_vector))

5.运行截图