cpcg-intro

Functions in this package serve the purpose of solving for x in Ax = b, where A is a n × n symmetric and positive definite matrix, b is a n × 1 column vector.

To improve scalability of conjugate gradient methods for larger matrices, the C++ Armadillo templated linear algebra library is used for the implementation. The package also provides flexibility to have user-specified preconditioner options to cater for different optimization needs.

This vignette will walk through some simple examples for using main functions in the package.

1. cgsolve(): Conjugate gradient method

The idea of conjugate gradient method is to find a set of mutually conjugate directions for the unconstrained problem arg minxf(x) where f(x) = 0.5yTΣy − yx + z and z is a constant. The problem is equivalent to solving Σx = y.

This function implements an iterative procedure to reduce the number of matrix-vector multiplications. The conjugate gradient method improves memory efficiency and computational complexity, especially when Σ is relatively sparse.

Example: Solve Ax = b where $A = \begin{bmatrix} 4 & 1 \\ 1 & 3 \end{bmatrix}$, $b = \begin{bmatrix} 1 \\ 2 \end{bmatrix}$.

test_A <- matrix(c(4,1,1,3), ncol = 2)
test_b <- matrix(1:2, ncol = 1)

cgsolve(test_A, test_b, 1e-6, 1000)

2. pcgsolve(): Preconditioned conjugate gradient method

When the condition number for Σ is large, the conjugate gradient (CG) method may fail to converge in a reasonable number of iterations. The Preconditioned Conjugate Gradient (PCG) Method applies a precondition matrix C and approaches the problem by solving: C−1Σx = C−1y where the symmetric and positive-definite matrix C approximates Σ and C−1Σ improves the condition number of Σ.

Choices for the preconditioner include: Jacobi preconditioning (Jacobi), symmetric successive over-relaxation (SSOR), and incomplete Cholesky factorization (ICC).
Example revisited: Now we solve the same problem using incomplete Cholesky factorization of A as preconditioner.

test_A <- matrix(c(4,1,1,3), ncol = 2)
test_b <- matrix(1:2, ncol = 1)

pcgsolve(test_A, test_b, "ICC")

Check Github repo and cPCG: Efficient and Customized Preconditioned Conjugate Gradient Method for more information.