Iterative Solution of Symmetric Positive-Definite Linear Systems Using Conjugate Gradients

Hard

The Conjugate Gradient (CG) method is one of the most elegant and powerful iterative algorithms in numerical linear algebra. It is specifically designed to solve systems of linear equations of the form Ax = b, where A is a symmetric positive-definite (SPD) matrix. Unlike direct methods such as Gaussian elimination or LU decomposition, CG operates iteratively, making it highly efficient for large, sparse systems where direct methods become computationally prohibitive. Mathematical Foundation The CG method is rooted in optimization theory. Solving Ax = b when A is SPD is mathematically equivalent to minimizing the quadratic functional: $$f(x) = \frac{1}{2}x^\top A x - b^\top x$$ The minimum of this function occurs precisely at the solution x* = A⁻¹b. The CG algorithm navigates through the solution space by generating a sequence of A-conjugate search directions, ensuring that each iteration makes progress toward the solution without undoing previous work. Algorithm Overview The algorithm maintains several key vectors at each iteration: x_k: The current approximation to the solution r_k: The residual vector (r_k = b - Ax_k), measuring how far we are from the true solution p_k: The search direction, constructed to be A-conjugate to all previous directions α_k: The step length along the current search direction β_k: The scaling factor for updating the search direction At each iteration, the algorithm: Computes the step length α_k that minimizes f(x) along the current search direction Updates the solution estimate: x_{k+1} = x_k + α_k * p_k Updates the residual: r_{k+1} = r_k - α_k * A * p_k Computes β_k and updates the search direction for the next iteration Convergence Properties For an n×n SPD matrix, the CG method is guaranteed to converge to the exact solution (in exact arithmetic) in at most n iterations. In practice, the method often converges much faster, especially when the matrix has clustered eigenvalues or when a preconditioner is applied. The convergence rate depends on the condition number κ(A) = λ_max / λ_min, where λ_max and λ_min are the largest and smallest eigenvalues of A. Your Task Implement a function that solves the linear system Ax = b using the Conjugate Gradient method. Your implementation should: Accept a symmetric, positive-definite matrix A and a right-hand side vector b Support an optional initial guess for the solution (default to the zero vector) Iterate until either the residual norm falls below a specified tolerance or a maximum iteration count is reached Return the approximate solution vector

Example 1

Input

A = [[4.0, 1.0], [1.0, 3.0]]
b = [1.0, 2.0]
max_iterations = 5
convergence_threshold = 1e-8

Output

[0.09090909, 0.63636364]

Explanation

The 2×2 matrix A is symmetric positive-definite (both eigenvalues are positive). Starting from the zero vector, the CG method computes: **Iteration 1:** • Initial residual r₀ = b - A·x₀ = [1.0, 2.0] (since x₀ = [0, 0]) • Initial search direction p₀ = r₀ = [1.0, 2.0] • Compute A·p₀ = [4·1 + 1·2, 1·1 + 3·2] = [6, 7] • Step length α₀ = (r₀·r₀) / (p₀·A·p₀) = 5 / 20 = 0.25 • Update x₁ = x₀ + α₀·p₀ = [0.25, 0.5] • Update r₁ = r₀ - α₀·A·p₀ = [1, 2] - 0.25·[6, 7] = [-0.5, 0.25] **Iteration 2:** • The algorithm continues, refining the solution further • After 2 iterations (n = 2), the exact solution is reached The final solution x = [0.09090909, 0.63636364] satisfies Ax = b with the given matrix.

Example 2

Input

A = [[4.0, 1.0, 0.0], [1.0, 5.0, 1.0], [0.0, 1.0, 4.0]]
b = [1.0, 2.0, 3.0]
max_iterations = 10
convergence_threshold = 1e-8

Output

[0.19444444, 0.22222222, 0.69444444]

Explanation

This is a tridiagonal SPD matrix, which is common in finite difference discretizations. The matrix structure: A = | 4 1 0 | | 1 5 1 | | 0 1 4 | The off-diagonal elements represent coupling between adjacent variables, while the dominant diagonal ensures positive-definiteness. The CG method converges in at most 3 iterations (the matrix dimension), with the solution vector balancing the right-hand side values according to the system's structure.

Example 3

Input

A = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
b = [5.0, 10.0, 15.0]
max_iterations = 5
convergence_threshold = 1e-8

Output

[5.0, 10.0, 15.0]

Explanation

When A is the identity matrix, the system Ax = b simplifies to Ix = b, meaning x = b. The CG method converges in exactly 1 iteration: • Starting from x₀ = [0, 0, 0], the residual r₀ = b = [5, 10, 15] • Since A = I, the product A·r₀ = r₀, and the step length α₀ = 1 • Update x₁ = x₀ + 1·r₀ = [5, 10, 15] — the exact solution This demonstrates that for well-conditioned systems (condition number = 1 for identity), CG converges immediately.

Constraints

1 ≤ n ≤ 100 (matrix dimension n×n)
1 ≤ max_iterations ≤ 1000
1e-15 ≤ convergence_threshold ≤ 1e-4
The matrix A is guaranteed to be symmetric positive-definite
-10⁶ ≤ A[i][j] ≤ 10⁶
-10⁶ ≤ b[i] ≤ 10⁶
If initial_guess is not provided, use the zero vector

Mathematical Foundation

The CG method is rooted in optimization theory. Solving Ax = b when A is SPD is mathematically equivalent to minimizing the quadratic functional:

$$f(x) = \frac{1}{2}x^\top A x - b^\top x$$

The minimum of this function occurs precisely at the solution x* = A⁻¹b. The CG algorithm navigates through the solution space by generating a sequence of A-conjugate search directions, ensuring that each iteration makes progress toward the solution without undoing previous work.

Algorithm Overview

The algorithm maintains several key vectors at each iteration:

x_k: The current approximation to the solution

r_k: The residual vector (r_k = b - Ax_k), measuring how far we are from the true solution

p_k: The search direction, constructed to be A-conjugate to all previous directions

α_k: The step length along the current search direction

β_k: The scaling factor for updating the search direction

At each iteration, the algorithm:

Computes the step length α_k that minimizes f(x) along the current search direction

Updates the solution estimate: x_{k+1} = x_k + α_k * p_k

Updates the residual: r_{k+1} = r_k - α_k * A * p_k

Computes β_k and updates the search direction for the next iteration

Convergence Properties

For an n×n SPD matrix, the CG method is guaranteed to converge to the exact solution (in exact arithmetic) in at most n iterations. In practice, the method often converges much faster, especially when the matrix has clustered eigenvalues or when a preconditioner is applied. The convergence rate depends on the condition number κ(A) = λ_max / λ_min, where λ_max and λ_min are the largest and smallest eigenvalues of A.

Your Task

Implement a function that solves the linear system Ax = b using the Conjugate Gradient method. Your implementation should:

Accept a symmetric, positive-definite matrix A and a right-hand side vector b

Support an optional initial guess for the solution (default to the zero vector)

Iterate until either the residual norm falls below a specified tolerance or a maximum iteration count is reached

Return the approximate solution vector