resources https en wikipedia orgwikiconjugategradientmethod https en wikipedia orgwikidefinitesymmetricmatrix returns true if input matrix is symmetric positive definite returns false otherwise for a matrix to be spd all eigenvalues must be positive import numpy as np matrix np array 4 12401784 5 01453636 0 63865857 5 01453636 12 33347422 3 40493586 0 63865857 3 40493586 5 78591885 ismatrixspdmatrix true matrix np array 0 34634879 1 96165514 2 18277744 0 74074469 1 19648894 1 34223498 0 7687067 0 06018373 1 16315631 ismatrixspdmatrix false ensure matrix is square if matrix not symmetric exit right away get eigenvalues and eignevectors for a symmetric matrix check sign of all eigenvalues np all returns a value of type np bool returns a symmetric positive definite matrix given a dimension input dimension gives the square matrix dimension output spdmatrix is an diminesion x dimensions symmetric positive definite spd matrix import numpy as np dimension 3 spdmatrix createspdmatrixdimension ismatrixspdspdmatrix true returns solution to the linear system np dotspdmatrix x b input spdmatrix is an nxn symmetric positive definite spd matrix loadvector is an nx1 vector output x is an nx1 vector that is the solution vector import numpy as np spdmatrix np array 8 73256573 5 02034289 2 68709226 5 02034289 3 78188322 0 91980451 2 68709226 0 91980451 1 94746467 b np array 5 80872761 3 23807431 1 95381422 conjugategradientspdmatrix b array0 63114139 0 01561498 0 13979294 ensure proper dimensionality initialize solution guess residual search direction set initial errors in solution guess and residual set iteration counter to threshold number of iterations save this value so we only calculate the matrixvector product once the main algorithm update search direction magnitude update solution guess calculate new residual calculate new krylov subspace scale calculate new a conjuage search direction calculate errors update variables update number of iterations testconjugategradient self running tests create linear system with spd matrix and known solution xtrue numpy solution our implementation ensure both solutions are close to xtrue and therefore one another returns true if input matrix is symmetric positive definite returns false otherwise for a matrix to be spd all eigenvalues must be positive import numpy as np matrix np array 4 12401784 5 01453636 0 63865857 5 01453636 12 33347422 3 40493586 0 63865857 3 40493586 5 78591885 _is_matrix_spd matrix true matrix np array 0 34634879 1 96165514 2 18277744 0 74074469 1 19648894 1 34223498 0 7687067 0 06018373 1 16315631 _is_matrix_spd matrix false ensure matrix is square if matrix not symmetric exit right away get eigenvalues and eignevectors for a symmetric matrix check sign of all eigenvalues np all returns a value of type np bool_ returns a symmetric positive definite matrix given a dimension input dimension gives the square matrix dimension output spd_matrix is an diminesion x dimensions symmetric positive definite spd matrix import numpy as np dimension 3 spd_matrix _create_spd_matrix dimension _is_matrix_spd spd_matrix true returns solution to the linear system np dot spd_matrix x b input spd_matrix is an nxn symmetric positive definite spd matrix load_vector is an nx1 vector output x is an nx1 vector that is the solution vector import numpy as np spd_matrix np array 8 73256573 5 02034289 2 68709226 5 02034289 3 78188322 0 91980451 2 68709226 0 91980451 1 94746467 b np array 5 80872761 3 23807431 1 95381422 conjugate_gradient spd_matrix b array 0 63114139 0 01561498 0 13979294 ensure proper dimensionality initialize solution guess residual search direction set initial errors in solution guess and residual set iteration counter to threshold number of iterations save this value so we only calculate the matrix vector product once the main algorithm update search direction magnitude update solution guess calculate new residual calculate new krylov subspace scale calculate new a conjuage search direction calculate errors update variables update number of iterations test_conjugate_gradient self running tests create linear system with spd matrix and known solution x_true numpy solution our implementation ensure both solutions are close to x_true and therefore one another
from typing import Any

import numpy as np


def _is_matrix_spd(matrix: np.ndarray) -> bool:
    
    
    assert np.shape(matrix)[0] == np.shape(matrix)[1]

    
    if np.allclose(matrix, matrix.T) is False:
        return False

    
    eigen_values, _ = np.linalg.eigh(matrix)

    
    
    return bool(np.all(eigen_values > 0))


def _create_spd_matrix(dimension: int) -> Any:
    
    random_matrix = np.random.randn(dimension, dimension)
    spd_matrix = np.dot(random_matrix, random_matrix.T)
    assert _is_matrix_spd(spd_matrix)
    return spd_matrix


def conjugate_gradient(
    spd_matrix: np.ndarray,
    load_vector: np.ndarray,
    max_iterations: int = 1000,
    tol: float = 1e-8,
) -> Any:
    
    
    assert np.shape(spd_matrix)[0] == np.shape(spd_matrix)[1]
    assert np.shape(load_vector)[0] == np.shape(spd_matrix)[0]
    assert _is_matrix_spd(spd_matrix)

    
    x0 = np.zeros((np.shape(load_vector)[0], 1))
    r0 = np.copy(load_vector)
    p0 = np.copy(r0)

    
    error_residual = 1e9
    error_x_solution = 1e9
    error = 1e9

    
    iterations = 0

    while error > tol:
        
        w = np.dot(spd_matrix, p0)

        

        
        alpha = np.dot(r0.T, r0) / np.dot(p0.T, w)
        
        x = x0 + alpha * p0
        
        r = r0 - alpha * w
        
        beta = np.dot(r.T, r) / np.dot(r0.T, r0)
        
        p = r + beta * p0

        
        error_residual = np.linalg.norm(r - r0)
        error_x_solution = np.linalg.norm(x - x0)
        error = np.maximum(error_residual, error_x_solution)

        
        x0 = np.copy(x)
        r0 = np.copy(r)
        p0 = np.copy(p)

        
        iterations += 1
        if iterations > max_iterations:
            break

    return x


def test_conjugate_gradient() -> None:
    
    
    dimension = 3
    spd_matrix = _create_spd_matrix(dimension)
    x_true = np.random.randn(dimension, 1)
    b = np.dot(spd_matrix, x_true)

    
    x_numpy = np.linalg.solve(spd_matrix, b)

    
    x_conjugate_gradient = conjugate_gradient(spd_matrix, b)

    
    assert np.linalg.norm(x_numpy - x_true) <= 1e-6
    assert np.linalg.norm(x_conjugate_gradient - x_true) <= 1e-6


if __name__ == "__main__":
    import doctest

    doctest.testmod()
    test_conjugate_gradient()