基于前后向替换的边框形式Cholesky分解Python实现需求

I'll walk through implementing the two required functions using your existing forward_subs and back_subs functions, following the bordered form approach and leveraging substitution for solving triangular systems.

Approach for cholesky_factor (Bordered Form) #

The bordered form builds the upper triangular Cholesky factor ( R ) column by column:

1. Start with the first diagonal element: ( R[0,0] = \sqrt{A[0,0]} )
2. For each subsequent column ( j ):
   • Extract the first ( j ) elements of column ( j ) from ( A )
   • Solve the lower triangular system ( R_{0:j,0:j}^T c = A_{0:j,j} ) using forward substitution (since ( R^T ) is lower triangular)
   • Compute the diagonal element ( R[j,j] = \sqrt{A[j,j] - c^T c} )
   • Assign ( c ) to the first ( j ) positions of column ( j ) in ( R )

Approach for solve_system_cholesky #

To solve ( Ax = b ) where ( A = R^T R ):

1. First solve ( R^T y = b ) using forward substitution (lower triangular system)
2. Then solve ( R x = y ) using back substitution (upper triangular system)

Full Code #

import numpy as np
from copy import copy

# Existing forward substitution function (solves Gx = b for lower triangular G)
def forward_subs(G, b):
    x = []
    for i in range(len(b)):
        x.append(b[i])
        for j in range(i):
            x[i] = x[i] - (G[i][j] * x[j])
        x[i] = x[i] / G[i][i]
    return x

# Existing back substitution function (solves Gx = b for upper triangular G)
def back_subs(G, b):
    n = b.size
    x = np.zeros_like(b)
    if G[n-1, n-1] == 0:
        raise ValueError("Matrix is singular")
    x[n-1] = b[n-1] / G[n-1, n-1]
    for i in range(n-2, -1, -1):
        bb = 0
        for j in range(i+1, n):
            bb += G[i, j] * x[j]
        x[i] = (b[i] - bb) / G[i][i]
    return x

# Bordered Cholesky decomposition function
def cholesky_factor(A):
    n = A.shape[0]
    R = np.zeros_like(A)
    
    for j in range(n):
        if j == 0:
            # Initialize first diagonal element
            R[j, j] = np.sqrt(A[j, j])
        else:
            # Extract the upper part of column j from A
            a = A[:j, j]
            # Solve lower triangular system R^T[:j,:j] * c = a
            c = np.array(forward_subs(R[:j, :j].T, a))
            # Assign c to the upper part of column j in R
            R[:j, j] = c
            # Compute diagonal element ensuring positive definiteness
            R[j, j] = np.sqrt(A[j, j] - np.dot(c, c))
    
    return R

# Solve linear system using Cholesky factor
def solve_system_cholesky(R, b):
    # Step 1: Solve R^T y = b using forward substitution
    y = np.array(forward_subs(R.T, b))
    # Step 2: Solve R x = y using back substitution
    x = back_subs(R, y)
    return x

# Test with the example matrix
if __name__ == "__main__":
    A = np.array([[4, 0], [0, 9]], dtype=np.float64)
    b = np.array([4, 18], dtype=np.float64)
    
    # Compute Cholesky factor
    R = cholesky_factor(A)
    print("Cholesky Factor R:")
    print(R)
    
    # Solve Ax = b
    x = solve_system_cholesky(R, b)
    print("\nSolution x:")
    print(x)
    
    # Verify solution correctness
    print("\nVerification Ax = b:")
    print(np.dot(A, x))

Test Output Explanation #

• The Cholesky factor ( R ) will be ( \begin{bmatrix} 2 & 0 \ 0 & 3 \end{bmatrix} ), which is correct since ( R^T R = A )
• The solution ( x ) for ( b = [4, 18] ) is ( [1, 2] ), as ( 41 + 02 = 4 ) and ( 01 +92=18 )
• The verification step confirms that ( Ax ) equals ( b ), validating the solution

内容的提问来源于stack exchange，提问作者sjoe17