View Javadoc
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   */
17  
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;
23  
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;
29  
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 {
82  
83      /** Key for the <a href="#context">exception context</a>. */
84      public static final String OPERATOR = "operator";
85  
86      /** Key for the <a href="#context">exception context</a>. */
87      public static final String VECTOR = "vector";
88  
89      /**
90       * {@code true} if positive-definiteness of matrix and preconditioner should
91       * be checked.
92       */
93      private boolean check;
94  
95      /** The value of &delta;, for the default stopping criterion. */
96      private final double delta;
97  
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     }
113 
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     }
131 
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     }
142 
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);
162 
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);
172 
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 }