Package com.imsl.math

Class ConjugateGradient

java.lang.Object
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}\; n\gg 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,\rm{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 & \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:

  • Constructor Details

    • 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 Details