4. 平方根法（Cholesky Decomposition）

平方根法（Cholesky Decomposition）是一种专门用于对称正定矩阵的分解方法。它将矩阵分解为一个下三角矩阵及其转置的乘积，即 ( A = LL^T )，其中 ( L ) 是一个下三角矩阵。

下面是使用Python实现Cholesky分解并求解线性方程组的代码。

实现步骤

1. Cholesky分解：将对称正定矩阵 ( A ) 分解为下三角矩阵 ( L ) 和其转置矩阵 ( L^T )。

2. 前向替代（Forward Substitution）：解决方程 ( Ly = b )。

3. 回代（Backward Substitution）：解决方程 ( L^Tx = y )。

代码实现

import numpy as np

def cholesky_decomposition(A):
    """
    Perform Cholesky decomposition of matrix A.

    Parameters:
    A : 2D list or numpy array
        Symmetric positive-definite matrix

    Returns:
    L : numpy array
        Lower triangular matrix such that A = L * L.T
    """
    n = len(A)
    L = np.zeros_like(A)

    for i in range(n):
        for j in range(i + 1):
            sum_k = sum(L[i][k] * L[j][k] for k in range(j))
            if i == j:
                L[i][j] = np.sqrt(A[i][i] - sum_k)
            else:
                L[i][j] = (A[i][j] - sum_k) / L[j][j]

    return L

def forward_substitution(L, b):
    """
    Perform forward substitution to solve Ly = b.

    Parameters:
    L : 2D list or numpy array
        Lower triangular matrix
    b : 1D list or numpy array
        Right-hand side vector

    Returns:
    y : numpy array
        Solution vector
    """
    n = len(b)
    y = np.zeros(n)

    for i in range(n):
        y[i] = (b[i] - sum(L[i][j] * y[j] for j in range(i))) / L[i][i]

    return y

def backward_substitution(L_T, y):
    """
    Perform backward substitution to solve L.T x = y.

    Parameters:
    L_T : 2D list or numpy array
        Transpose of the lower triangular matrix
    y : 1D list or numpy array
        Solution vector from forward substitution

    Returns:
    x : numpy array
        Solution vector
    """
    n = len(y)
    x = np.zeros(n)

    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - sum(L_T[i][j] * x[j] for j in range(i + 1, n))) / L_T[i][i]

    return x

def solve_cholesky_decomposition(A, b):
    """
    Solve the system of linear equations Ax = b using Cholesky Decomposition.

    Parameters:
    A : 2D list or numpy array
        Symmetric positive-definite matrix
    b : 1D list or numpy array
        Right-hand side vector

    Returns:
    x : numpy array
        Solution vector
    """
    L = cholesky_decomposition(A)
    y = forward_substitution(L, b)
    x = backward_substitution(L.T, y)
    return x

# Example usage
A = np.array([[4, 12, -16],
              [12, 37, -43],
              [-16, -43, 98]], dtype=float)
b = np.array([1, 2, 3], dtype=float)

x = solve_cholesky_decomposition(A, b)
print("Solution:", x)

代码说明

1. Cholesky分解（cholesky_decomposition函数）：

   • 对称正定矩阵 ( A ) 被分解为下三角矩阵 ( L ) 和其转置矩阵 ( L^T )。

2. 前向替代（forward_substitution函数）：

   • 使用下三角矩阵 ( L ) 和右端向量 ( b ) 求解中间变量 ( y )。

3. 回代（backward_substitution函数）：

   • 使用下三角矩阵的转置 ( L^T ) 和中间变量 ( y ) 求解最终解 ( x )。

4. 求解线性方程组（solve_cholesky_decomposition函数）：

   • 综合上述步骤，完成Cholesky分解并求解线性方程组。

示例

给定对称正定矩阵 (A) 和右端向量 (b)：

A = [[4, 12, -16],
     [12, 37, -43],
     [-16, -43, 98]]
b = [1, 2, 3]

运行上述代码，求解线性方程组 (Ax = b)，可以得到解向量 x。