ConjugateGradient.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.linear;

  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.exception.MathIllegalStateException;
  27. import org.hipparchus.exception.NullArgumentException;
  28. import org.hipparchus.util.IterationManager;

  30. /**
  31.  * <p>
  32.  * This is an implementation of the conjugate gradient method for
  33.  * {@link RealLinearOperator}. It follows closely the template by <a
  34.  * href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The linear system at
  35.  * hand is A &middot; x = b, and the residual is r = b - A &middot; x.
  36.  * </p>
  37.  * <p><strong>Default stopping criterion</strong></p>
  38.  * <p>
  39.  * A default stopping criterion is implemented. The iterations stop when || r ||
  40.  * &le; &delta; || b ||, where b is the right-hand side vector, r the current
  41.  * estimate of the residual, and &delta; a user-specified tolerance. It should
  42.  * be noted that r is the so-called <em>updated</em> residual, which might
  43.  * differ from the true residual due to rounding-off errors (see e.g. <a
  44.  * href="#STRA2002">Strakos and Tichy, 2002</a>).
  45.  * </p>
  46.  * <p><strong>Iteration count</strong></p>
  47.  * <p>
  48.  * In the present context, an iteration should be understood as one evaluation
  49.  * of the matrix-vector product A &middot; x. The initialization phase therefore
  50.  * counts as one iteration.
  51.  * </p>
  52.  * <p><strong><a id="context">Exception context</a></strong></p>
  53.  * <p>
  54.  * Besides standard {@link MathIllegalArgumentException}, this class might throw
  55.  * {@link MathIllegalArgumentException} if the linear operator or
  56.  * the preconditioner are not positive definite.
  57.  * </p>
  58.  * <ul>
  59.  * <li>key {@code "operator"} points to the offending linear operator, say L,</li>
  60.  * <li>key {@code "vector"} points to the offending vector, say x, such that
  61.  * x<sup>T</sup> &middot; L &middot; x &lt; 0.</li>
  62.  * </ul>
  63.  * <p><strong>References</strong></p>
  64.  * <dl>
  65.  * <dt>Barret et al. (1994)</dt>
  66.  * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra,
  67.  * V. Eijkhout, R. Pozo, C. Romine and H. Van der Vorst,
  68.  * <a href="http://www.netlib.org/linalg/html_templates/Templates.html"><em>
  69.  * Templates for the Solution of Linear Systems: Building Blocks for Iterative
  70.  * Methods</em></a>, SIAM</dd>
  71.  * <dt>Strakos and Tichy (2002)</dt>
  72.  * <dd>Z. Strakos and P. Tichy, <a
  73.  * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf">
  74.  * <em>On error estimation in the conjugate gradient method and why it works
  75.  * in finite precision computations</em></a>, Electronic Transactions on
  76.  * Numerical Analysis 13: 56-80, 2002</dd>
  77.  * </dl>
  78.  *
  79.  */
  80. public class ConjugateGradient
  81.     extends PreconditionedIterativeLinearSolver {

  83.     /** Key for the <a href="#context">exception context</a>. */
  84.     public static final String OPERATOR = "operator";

  86.     /** Key for the <a href="#context">exception context</a>. */
  87.     public static final String VECTOR = "vector";

  89.     /**
  90.      * {@code true} if positive-definiteness of matrix and preconditioner should
  91.      * be checked.
  92.      */
  93.     private boolean check;

  95.     /** The value of &delta;, for the default stopping criterion. */
  96.     private final double delta;

  98.     /**
  99.      * Creates a new instance of this class, with <a href="#stopcrit">default
  100.      * stopping criterion</a>.
  101.      *
  102.      * @param maxIterations the maximum number of iterations
  103.      * @param delta the &delta; parameter for the default stopping criterion
  104.      * @param check {@code true} if positive definiteness of both matrix and
  105.      * preconditioner should be checked
  106.      */
  107.     public ConjugateGradient(final int maxIterations, final double delta,
  108.                              final boolean check) {
  109.         super(maxIterations);
  110.         this.delta = delta;
  111.         this.check = check;
  112.     }

  114.     /**
  115.      * Creates a new instance of this class, with <a href="#stopcrit">default
  116.      * stopping criterion</a> and custom iteration manager.
  117.      *
  118.      * @param manager the custom iteration manager
  119.      * @param delta the &delta; parameter for the default stopping criterion
  120.      * @param check {@code true} if positive definiteness of both matrix and
  121.      * preconditioner should be checked
  122.      * @throws NullArgumentException if {@code manager} is {@code null}
  123.      */
  124.     public ConjugateGradient(final IterationManager manager,
  125.                              final double delta, final boolean check)
  126.         throws NullArgumentException {
  127.         super(manager);
  128.         this.delta = delta;
  129.         this.check = check;
  130.     }

  132.     /**
  133.      * Returns {@code true} if positive-definiteness should be checked for both
  134.      * matrix and preconditioner.
  135.      *
  136.      * @return {@code true} if the tests are to be performed
  137.      * @since 1.4
  138.      */
  139.     public final boolean shouldCheck() {
  140.         return check;
  141.     }

  143.     /**
  144.      * {@inheritDoc}
  145.      *
  146.      * @throws MathIllegalArgumentException if {@code a} or {@code m} is
  147.      * not positive definite
  148.      */
  149.     @Override
  150.     public RealVector solveInPlace(final RealLinearOperator a,
  151.                                    final RealLinearOperator m,
  152.                                    final RealVector b,
  153.                                    final RealVector x0)
  154.         throws MathIllegalArgumentException, NullArgumentException,
  155.         MathIllegalStateException {
  156.         checkParameters(a, m, b, x0);
  157.         final IterationManager manager = getIterationManager();
  158.         // Initialization of default stopping criterion
  159.         manager.resetIterationCount();
  160.         final double rmax = delta * b.getNorm();
  161.         final RealVector bro = RealVector.unmodifiableRealVector(b);

  163.         // Initialization phase counts as one iteration.
  164.         manager.incrementIterationCount();
  165.         // p and x are constructed as copies of x0, since presumably, the type
  166.         // of x is optimized for the calculation of the matrix-vector product
  167.         // A.x.
  168.         final RealVector x = x0;
  169.         final RealVector xro = RealVector.unmodifiableRealVector(x);
  170.         final RealVector p = x.copy();
  171.         RealVector q = a.operate(p);

  173.         final RealVector r = b.combine(1, -1, q);
  174.         final RealVector rro = RealVector.unmodifiableRealVector(r);
  175.         double rnorm = r.getNorm();
  176.         RealVector z;
  177.         if (m == null) {
  178.             z = r;
  179.         } else {
  180.             z = null;
  181.         }
  182.         IterativeLinearSolverEvent evt;
  183.         evt = new DefaultIterativeLinearSolverEvent(this,
  184.             manager.getIterations(), xro, bro, rro, rnorm);
  185.         manager.fireInitializationEvent(evt);
  186.         if (rnorm <= rmax) {
  187.             manager.fireTerminationEvent(evt);
  188.             return x;
  189.         }
  190.         double rhoPrev = 0.;
  191.         while (true) {
  192.             manager.incrementIterationCount();
  193.             evt = new DefaultIterativeLinearSolverEvent(this,
  194.                 manager.getIterations(), xro, bro, rro, rnorm);
  195.             manager.fireIterationStartedEvent(evt);
  196.             if (m != null) {
  197.                 z = m.operate(r);
  198.             }
  199.             final double rhoNext = r.dotProduct(z);
  200.             if (check && (rhoNext <= 0.)) {
  201.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_POSITIVE_DEFINITE_OPERATOR);
  202.             }
  203.             if (manager.getIterations() == 2) {
  204.                 p.setSubVector(0, z);
  205.             } else {
  206.                 p.combineToSelf(rhoNext / rhoPrev, 1., z);
  207.             }
  208.             q = a.operate(p);
  209.             final double pq = p.dotProduct(q);
  210.             if (check && (pq <= 0.)) {
  211.                 throw new MathIllegalArgumentException(LocalizedCoreFormats.NON_POSITIVE_DEFINITE_OPERATOR);
  212.             }
  213.             final double alpha = rhoNext / pq;
  214.             x.combineToSelf(1., alpha, p);
  215.             r.combineToSelf(1., -alpha, q);
  216.             rhoPrev = rhoNext;
  217.             rnorm = r.getNorm();
  218.             evt = new DefaultIterativeLinearSolverEvent(this,
  219.                 manager.getIterations(), xro, bro, rro, rnorm);
  220.             manager.fireIterationPerformedEvent(evt);
  221.             if (rnorm <= rmax) {
  222.                 manager.fireTerminationEvent(evt);
  223.                 return x;
  224.             }
  225.         }
  226.     }
  227. }
Created with JaCoCo 0.8.12.202403310830