JMSLTM Numerical Library 6.0

com.imsl.math
Class ConjugateGradient

java.lang.Object
  extended by com.imsl.math.ConjugateGradient
All Implemented Interfaces:
Serializable

public class ConjugateGradient
extends Object
implements Serializable

Solves a real symmetric definite linear system using the conjugate gradient method with optional preconditioning.

Class ConjugateGradient solves the symmetric positive or negative definite linear system Ax = b using the conjugate gradient method with optional preconditioning. This method is described in detail by Golub and Van Loan (1983, Chapter 10), and in Hageman and Young (1981, Chapter 7).

The preconditioning matrix M is a matrix that approximates A, and for which the linear system Mz=r is easy to solve. These two properties are in conflict; balancing them is a topic of current research. If no preconditioning matrix is specified, M is set to the identity, i.e. M=I.

The number of iterations needed depends on the matrix and the error tolerance. As a rough guide,

{rm{itmax}}=sqrt{n};mbox{for}; ngg 1,

where n is the order of matrix A.

See the references for details.

Let M be the preconditioning matrix, let b,p,r,x and z be vectors and let tau be the desired relative error. Then the algorithm used is as follows:

             hspace*{1cm}lambda = -1

             hspace*{1cm}p_0 = x_0

             hspace*{1cm}r_1 = b-Ap_0

             hspace*{1cm}mbox{for } k=1,ldots,mbox{{tt{itmax}}}

                  hspace*{2cm}z_k = M^{-1}r_k

                  hspace*{2cm}mbox{if }k=1 mbox{ then}

                       hspace*{2.5cm}beta_k = 1

                       hspace*{2.5cm}p_k = z_k

                  hspace*{2cm}mbox{else}

                       hspace*{2.5cm}beta_k=(z_k^Tr_k)/(z_{k-1}^Tr_{k-1})

                       hspace*{2.5cm}p_k=z_k+beta_kp_{k-1}

                  hspace*{2cm}mbox{endif}

                  hspace*{2cm}alpha_k=(r_{k}^Tz_{k})/(p_k^TAp_k)

                  hspace*{2cm}x_k=x_{k-1}+alpha_kp_k

                  hspace*{2cm}r_{k+1}=r_k-alpha_kAp_k

                  hspace*{2cm}mbox{if }(|Ap_k|_2 leq tau (1-lambda)|x_k|_2) mbox{ then}

                       hspace*{2.5cm}mbox{recompute }lambda

                       hspace*{2.5cm}mbox{if }(|Ap_k|_2 leq tau(1-lambda)|x_k|_2)mbox{ exit}

                  hspace*{2cm}mbox{endif}

             hspace*{1cm}mbox{endfor}

Here, lambda is an estimate of lambda_{max}(Gamma), the largest eigenvalue of the iteration matrix Gamma=I-M^{-1}A. The stopping criterion is based on the result (Hageman and Young 1981, pp. 148-151)

frac{|x_k-x|_M}{|x|_M} leq left(frac{1}{1-lambda_{max}(Gamma)}right)left(frac{|z_k|_M}{|x_k|_M}right),

where

|x|_M^2=x^TMx,.

It is also known that

lambda_{max}(T_1) leq lambda_{max}(T_2) leq ldots leq lambda_{max}(Gamma) lt 1,

where the T_l are the symmetric, tridiagonal matrices

T_l = left[ begin{array}{ccccc}
                        mu_1 & omega_2 & & & \
                        omega_2 & mu_2 & omega_3 & & \
                        & omega_3 & mu_3 & raisebox{-1ex}{ ddots } & \
                        & & ddots & ddots & omega_l \
                        & & & omega_l & mu_l
                end{array} right]

with mu_1=1-1/alpha_1 and, for k=2,ldots,l,

mu_k=1-beta_k/alpha_{k-1}-1/alpha_k quad mbox{and} quad omega_k=sqrt{beta_k}/alpha_{k-1}.

Usually, the eigenvalue computation is needed for only a few of the iterations.

See Also:
Example without preconditioning, Example with different preconditioners, Serialized Form

Nested Class Summary
static interface ConjugateGradient.Function
          Public interface for the user supplied function to ConjugateGradient.
static class ConjugateGradient.NoConvergenceException
          The conjugate gradient method did not converge within the allowed maximum number of iterations.
static class ConjugateGradient.NotDefiniteAMatrixException
          The input matrix A is indefinite, that is the matrix is not positive or negative definite.
static class ConjugateGradient.NotDefiniteJacobiPreconditionerException
          The Jacobi preconditioner is not strictly positive or negative definite.
static class ConjugateGradient.NotDefinitePreconditionMatrixException
          The Precondition matrix is indefinite.
static interface ConjugateGradient.Preconditioner
          Public interface for the user supplied function to ConjugateGradient used for preconditioning.
static class ConjugateGradient.SingularPreconditionMatrixException
          The Precondition matrix is singular.
 
Constructor Summary
ConjugateGradient(int n, ConjugateGradient.Function argF)
          Conjugate gradient constructor.
 
Method Summary
 int getIterations()
          Returns the number of iterations needed by the conjugate gradient algorithm.
 double[] getJacobi()
          Returns the Jacobi preconditioning matrix.
 int getMaxIterations()
          Returns the maximum number of iterations allowed.
 double getRelativeError(double errorRelative)
          Returns the relative error used for stopping the algorithm.
 void setJacobi(double[] diagonal)
          Defines a Jacobi preconditioner as the preconditioning matrix, that is, M is the diagonal of A.
 void setMaxIterations(int maxIterations)
          Sets the maximum number of iterations allowed.
 void setRelativeError(double tolerance)
          Sets the relative error used for stopping the algorithm.
 double[] solve(double[] b)
          Solves a real symmetric positive or negative definite system Ax=b using a conjugate gradient method with or without preconditioning.
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Constructor Detail

ConjugateGradient

public ConjugateGradient(int n,
                         ConjugateGradient.Function argF)
Conjugate gradient constructor.

Parameters:
n - an int scalar value defining the order of the matrix.
argF - a Function that defines the user-supplied function which computes z = Ap. If argF implements Preconditioner then right preconditioning is performed using this user supplied function. Otherwise, no preconditioning is performed. Note that argF can be used to act upon the coefficients of matrix A stored in different storage modes.
Method Detail

getIterations

public int getIterations()
Returns the number of iterations needed by the conjugate gradient algorithm.

Returns:
an int value indicating the number of iterations needed.

getJacobi

public double[] getJacobi()
Returns the Jacobi preconditioning matrix.

Returns:
a double vector diagonal containing the diagonal of the Jacobi preconditioner M, that is, diagonal[i]=A_{i,i}, A the input matrix.

getMaxIterations

public int getMaxIterations()
Returns the maximum number of iterations allowed.

Returns:
an int value specifying the maximum number of iterations allowed.

getRelativeError

public double getRelativeError(double errorRelative)
Returns the relative error used for stopping the algorithm.

Returns:
a double containing the relative error.

setJacobi

public void setJacobi(double[] diagonal)
Defines a Jacobi preconditioner as the preconditioning matrix, that is, M is the diagonal of A.

Parameters:
diagonal - a double vector containing the diagonal of A as the Jacobi preconditioner M, that is, diagonal[i]=A_{i,i}.
Throws:
IllegalArgumentException - is thrown if the length of vector diagonal is not equal to the order n of input matrix A.

setMaxIterations

public void setMaxIterations(int maxIterations)
Sets the maximum number of iterations allowed.

Parameters:
maxIterations - an int value specifying the maximum number of iterations allowed. By default, maxIterations = max(1000, sqrt{n}).
Throws:
IllegalArgumentException - is thrown if maxIterations is less than or equal to 0.

setRelativeError

public void setRelativeError(double tolerance)
Sets the relative error used for stopping the algorithm.

Parameters:
tolerance - a double specifying the relative error. By default, tolerance = 1.49e-08, the square root of the precision.
Throws:
IllegalArgumentException - is thrown if tolerance is less than 0.

solve

public double[] solve(double[] b)
               throws ConjugateGradient.SingularPreconditionMatrixException,
                      ConjugateGradient.NotDefinitePreconditionMatrixException,
                      SingularMatrixException,
                      ConjugateGradient.NotDefiniteAMatrixException,
                      ConjugateGradient.NoConvergenceException,
                      ConjugateGradient.NotDefiniteJacobiPreconditionerException
Solves a real symmetric positive or negative definite system Ax=b using a conjugate gradient method with or without preconditioning.

Parameters:
b - a double vector of length n containing the right-hand side.
Returns:
a double vector of length n containing the approximate solution to the linear system.
Throws:
IllegalArgumentException - is thrown if the length of b is not consistent with the order n of A.
ConjugateGradient.SingularPreconditionMatrixException - is thrown if the preconditioning matrix is singular.
ConjugateGradient.NotDefinitePreconditionMatrixException - is thrown if the preconditioning matrix is not definite.
SingularMatrixException - is thrown if input matrix A is singular.
ConjugateGradient.NotDefiniteAMatrixException - is thrown if matrix A is not definite.
ConjugateGradient.NoConvergenceException - is thrown if the algorithm is not convergent within maxIterations iterations.
ConjugateGradient.NotDefiniteJacobiPreconditionerException - is thrown if the Jacobi preconditioner is not definite.

JMSLTM Numerical Library 6.0

Copyright © 1970-2009 Visual Numerics, Inc.
Built September 1 2009.