4. The Conjugate Gradient Method¶

4.1. The `cg` Module¶

class cg.CG(op, **kwargs)¶

Bases: pykrylov.generic.generic.KrylovMethod

A pure Python implementation of the conjugate gradient (CG) algorithm. The conjugate gradient algorithm may be used to solve symmetric positive definite systems of linear equations, i.e., systems of the form

A x = b

where the operator A is square, symmetric and positive definite. This is equivalent to solving the unconstrained convex quadratic optimization problem

minimize -<b,x> + 1/2 <x, Ax>

in the variable x.

CG performs 1 operator-vector product, 2 dot products and 3 daxpys per iteration.

If a preconditioner is supplied, it needs to solve one preconditioning system per iteration. Our implementation is standard and follows [Kelley] and [Templates].

solve(rhs, **kwargs)¶

Solve a linear system with rhs as right-hand side by the CG method. The vector rhs should be a Numpy array.

Keywords:

check_symmetric:
guess:	Initial guess (Numpy array). Default: 0.
matvec_max:	Max. number of operator-vector produts. Default: 2n.
	Ensure operator is symmetric. Default: False.
check_curvature:
	Ensure operator is positive definite. Default: True.
store_resids:	Store full residual vector history. Default: False.
store_iterates:	Store full iterate history. Default: False.

4.2. Example¶

The following script runs through a list of symmetric positive definite matrices in Matrix Market format, builds the right-hand side so that the exact solution of the system is the vector of ones, and applies the conjugate gradient to each of them. The sparse matrices are handled by way of the Pysparse package.

from pykrylov.cg import CG
from pykrylov.linop import PysparseLinearOperator
from pysparse import spmatrix
from pysparse.pysparseMatrix import PysparseMatrix as sp
import numpy as np
from math import sqrt

matrices = ['1138bus.mtx', 'bcsstk08.mtx', 'bcsstk09.mtx', 'bcsstk10.mtx']
matrices += ['bcsstk11.mtx', 'bcsstk18.mtx', 'bcsstk19.mtx' ]

hdr = '%15s  %5s  %5s  %8s  %8s  %8s' % ('Name', 'Size', 'Mult',
                                         'Resid0', 'Resid', 'Error')
print hdr

for matName in matrices:
    A = sp( matrix=spmatrix.ll_mat_from_mtx(matName) )
    n = A.shape[0]
    e = np.ones(n)
    rhs = A * e
    cg = CG(PysparseLinearOperator(A), matvec_max=2*n)
    cg.solve(rhs)
    err = np.linalg.norm(cg.bestSolution - e)/sqrt(n)
    print '%15s  %5d  %5d  %8.2e  %8.2e  %8.2e' % (matName, n, cg.nMatvec,
                                                   cg.residNorm0,
                                                   cg.residNorm, err)

The script above produces the following output:

        Name   Size   Mult    Resid0     Resid     Error
 1138bus.mtx   1138   1759  1.46e+03  1.44e-03  1.30e-05
bcsstk08.mtx   1074   1255  8.74e+10  7.41e+04  7.54e-02
bcsstk09.mtx   1083    180  3.17e+08  1.84e+02  8.59e-06
bcsstk10.mtx   1086   1753  8.10e+07  7.30e+01  1.04e-03
bcsstk11.mtx   1473   1689  5.43e+09  4.88e+03  5.68e-02
bcsstk18.mtx  11948   9729  2.79e+11  2.49e+05  2.85e-01
bcsstk19.mtx    817    468  1.34e+15  1.26e+09  8.09e-01

Note that the default relative stopping tolerance is 1.0e-6 and is achieved in all cases.

4. The Conjugate Gradient Method¶

4.1. The cg Module¶

4.2. Example¶

4.1. The `cg` Module¶