改进平方根法解对称正定方程组

1.改进平方根法

对称阵的三角分解定理

设$A$为$n$阶对称矩阵，且$A$的所有顺序主子式均不为零，则$A$可唯一分解为$$A=LD{L}^T,$$
其中$L$为单位下三角矩阵，$D$为对角矩阵.即
$$A=\left[\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& \dots & 1\end{array}\right]\left[\begin{array}{cccc}{d}_1& & & \\ & d_2& & \\ & & \ddots & \\ & & & d_n\end{array}\right]\left[\begin{array}{cccc}1& l_{21}& \dots & l_{n1}\\ & 1& \dots & l_{n2}\\ & & \ddots & ⋮\\ & & & 1\end{array}\right]$$
由矩阵乘法，并注意$l_{jj}=1,l_{jk}=0(j<k),得$
$$a_{ij}=\sum _{k=1}^n(LD{)}_{ik}(L^T{)}_{kj}=\sum _{k=1}^n{l}_{ik}{d}_k{l}_{jk}=\sum _{k=1}^{j-1}{l}_{ik}{d}_k{l}_{jk}+l_{ij}{d}_j{l}_{jj}$$
于是得到计算$L$的元素及$D$的对角元素公式:
对于$i=1,2,\dots ,n.$
$(1)l_{ij}=(a_{ij}-\sum _{k=1}^{j-1}{l}_{ik}{d}_k{l}_{jk})/d_j,j=1,2,\dots ,i-1$
$(2)d_i=a_{ii}-\sum _{k=1}^{i-1}{l}_{ik}^2{d}_k.(i)$
为了避免重复计算，引进$$t_{ij}=l_{ij}{d}_j,$$由(i)式得到按行计算$L,T$元素的公式:$$d_1=a_{11}.$$
对于$i=2,3,\dots ,n.$
$(1)t_{ij}=a_{ij}-\sum _{k=1}^{j-1}{t}_{ik}{l}_{jk},j=1,2,\dots ,i-1;$
$(2)l_{ij}=t_{ij}/d_j,j=1,2,\dots ,i-1;$
$(3)d_i=a_{ii}-\sum _{k=1}^{i-1}{t}_{ik}{l}_{ik}.$
求解$Ly=b,D{L}^{T}x=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/d_n;\\ x_i=y_i/d_i-\sum _{k=i+1}^n{l}_{ki}{x}_k,i=n-1,n-2,\dots ,1.\end{array}$

2.Python实现改进平方根

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

                return inv(
                    lower_triangle_matrix @ diagonal_matrix @ lower_triangle_matrix.T) @ right_hand_side_vector
                # 自己写的发现第一个解x0=1.,而不是1.11111111，有精度损失，所以采用矩阵运算
                # # Ly=b,此处不指定双精度浮点数据类型，则可能会损失精度，导致结果偏差太大
                # y = np.ones_like(right_hand_side_vector, dtype=np.float64)
                # for i in range(rows):
                #     prod = 0
                #     for k in range(i):
                #         prod += lower_triangle_matrix[i, k] * y[k, 0]
                #     y[i, 0] = right_hand_side_vector[i, 0] - prod
                # # DL^Tx=y，此处不指定双精度浮点数据类型，则可能会损失精度，导致结果偏差太大
                # x = np.ones_like(right_hand_side_vector, dtype=np.float64)
                # x[rows - 1, 0] = y[rows - 1, 0] / diagonal_matrix[rows - 1, rows - 1]
                # for i in range(rows - 2, 0, -1):
                #     prod = 0
                #     for k in range(i + 1, rows):
                #         prod += lower_triangle_matrix[k, i] * x[k, 0]
                #     x[i, 0] = y[i, 0] / diagonal_matrix[i, i] - prod
                # return x

                # for i in range(rows):
                #     for j in range(i):
                #         prod = 0
                #         for k in range(j):
                #             prod += lower_triangle_matrix[i, k] * diagonal_matrix[k, k] * lower_triangle_matrix[j, k]
                #         lower_triangle_matrix[i, j] = (symmetric_positive_define_matrix[i, j] - prod) / diagonal_matrix[
                #             j, j]
                #     prod = 0
                #     for k in range(i):
                #         prod += lower_triangle_matrix[i, k] * diagonal_matrix[k, k] * lower_triangle_matrix[i, k]
                #     diagonal_matrix[i, i] = symmetric_positive_define_matrix[i, i] - prod
        else:
            raise Exception("error,the input matrix must be a symmetric-matrix")
    else:
        raise Exception("ERROR:please pass a square matrix.")

3.用改进平方根解对称正定方程组

$$\left[\begin{array}{ccc}2& -1& 1\\ -1& -2& 3\\ 1& 3& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}4\\ 5\\ 6\end{array}\right].$$

4.测试

if __name__ == '__main__':
    # 改进平方根法测试成功，来源详见李庆扬数值分析第5版P177，e.g.10
    symmetric_positive_define_matrix1 = np.array([[2, -1, 1], [-1, -2, 3], [1, 3, 1]])
    column_vector = np.array([4, 5, 6]).reshape((3, 1))
    print(improved_square_root(symmetric_positive_define_matrix1, column_vector))

5.运行截图