平方根法求解对称正定阵系数阵线方程组

### 平方根法（Cholesky分解）求解对称正定矩阵的线性方程组

对于给定的对称正定矩阵 \( A \)，以及向量 \( b \)，目标是找到满足 \( Ax = b \) 的未知向量 \( x \)[^1]。

通过 Cholesky 分解，可以将矩阵 \( A \) 表达为两个下三角矩阵相乘的形式：

\[ A = G G' \]

其中 \( G \) 是一个下三角矩阵。原方程组 \( Ax = b \) 可以转换成两步计算过程:

1. 解决下三角形方程组 \( Gy = b \)
2. 解决上三角形方程组 \( G'x = y \)

下面展示 Python 实现此算法的方式：

import numpy as np

def cholesky_decomposition(A):
    n = len(A)
    L = np.zeros((n, n))
    
    for i in range(n):
        for j in range(i+1):
            sumk = sum(L[i][k] * L[j][k] for k in range(j))
            
            if (i == j):  # Diagonal elements
                L[i][j] = np.sqrt(max(0, A[i][i] - sumk))
            else:
                L[i][j] = (A[i][j] - sumk) / L[j][j]
                
    return L


def solve_cholesky(A, b):
    L = cholesky_decomposition(A)
    LT = L.T
    
    # Solve Ly = b using forward substitution
    y = np.linalg.solve(L, b)
    
    # Solve LTx = y using backward substitution
    x = np.linalg.solve(LT, y)
    
    return x


# Example usage with a symmetric positive definite matrix and vector b
A_example = np.array([[6., 1., 0.], [1., 4., 1.], [0., 1., 14.]])
b_example = np.array([6., 24., 322.])

solution = solve_cholesky(A_example, b_example)
print(solution)

上述代码实现了基于 Cholesky 分解来解决对称正定线性系统的函数 `solve_cholesky` 和用于执行实际分解操作的支持函数 `cholesky_decomposition`[^2]。