diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java | 231 |
1 files changed, 231 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java b/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java new file mode 100644 index 0000000..61dab67 --- /dev/null +++ b/src/main/java/org/apache/commons/math3/linear/ConjugateGradient.java @@ -0,0 +1,231 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.commons.math3.linear; + +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.exception.MaxCountExceededException; +import org.apache.commons.math3.exception.NullArgumentException; +import org.apache.commons.math3.exception.util.ExceptionContext; +import org.apache.commons.math3.util.IterationManager; + +/** + * This is an implementation of the conjugate gradient method for {@link RealLinearOperator}. It + * follows closely the template by <a href="#BARR1994">Barrett et al. (1994)</a> (figure 2.5). The + * linear system at hand is A · x = b, and the residual is r = b - A · x. + * + * <h3><a id="stopcrit">Default stopping criterion</a></h3> + * + * <p>A default stopping criterion is implemented. The iterations stop when || r || ≤ δ || + * b ||, where b is the right-hand side vector, r the current estimate of the residual, and δ + * a user-specified tolerance. It should be noted that r is the so-called <em>updated</em> residual, + * which might differ from the true residual due to rounding-off errors (see e.g. <a + * href="#STRA2002">Strakos and Tichy, 2002</a>). + * + * <h3>Iteration count</h3> + * + * <p>In the present context, an iteration should be understood as one evaluation of the + * matrix-vector product A · x. The initialization phase therefore counts as one iteration. + * + * <h3><a id="context">Exception context</a></h3> + * + * <p>Besides standard {@link DimensionMismatchException}, this class might throw {@link + * NonPositiveDefiniteOperatorException} if the linear operator or the preconditioner are not + * positive definite. In this case, the {@link ExceptionContext} provides some more information + * + * <ul> + * <li>key {@code "operator"} points to the offending linear operator, say L, + * <li>key {@code "vector"} points to the offending vector, say x, such that x<sup>T</sup> + * · L · x < 0. + * </ul> + * + * <h3>References</h3> + * + * <dl> + * <dt><a id="BARR1994">Barret et al. (1994)</a> + * <dd>R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. M. Donato, J. Dongarra, V. Eijkhout, R. + * Pozo, C. Romine and H. Van der Vorst, <a + * href="http://www.netlib.org/linalg/html_templates/Templates.html"><em> Templates for the + * Solution of Linear Systems: Building Blocks for Iterative Methods</em></a>, SIAM + * <dt><a id="STRA2002">Strakos and Tichy (2002) + * <dt> + * <dd>Z. Strakos and P. Tichy, <a + * href="http://etna.mcs.kent.edu/vol.13.2002/pp56-80.dir/pp56-80.pdf"><em>On error estimation + * in the conjugate gradient method and why it works in finite precision + * computations</em></a>, Electronic Transactions on Numerical Analysis 13: 56-80, 2002 + * </dl> + * + * @since 3.0 + */ +public class ConjugateGradient extends PreconditionedIterativeLinearSolver { + + /** Key for the <a href="#context">exception context</a>. */ + public static final String OPERATOR = "operator"; + + /** Key for the <a href="#context">exception context</a>. */ + public static final String VECTOR = "vector"; + + /** {@code true} if positive-definiteness of matrix and preconditioner should be checked. */ + private boolean check; + + /** The value of δ, for the default stopping criterion. */ + private final double delta; + + /** + * Creates a new instance of this class, with <a href="#stopcrit">default stopping + * criterion</a>. + * + * @param maxIterations the maximum number of iterations + * @param delta the δ parameter for the default stopping criterion + * @param check {@code true} if positive definiteness of both matrix and preconditioner should + * be checked + */ + public ConjugateGradient(final int maxIterations, final double delta, final boolean check) { + super(maxIterations); + this.delta = delta; + this.check = check; + } + + /** + * Creates a new instance of this class, with <a href="#stopcrit">default stopping criterion</a> + * and custom iteration manager. + * + * @param manager the custom iteration manager + * @param delta the δ parameter for the default stopping criterion + * @param check {@code true} if positive definiteness of both matrix and preconditioner should + * be checked + * @throws NullArgumentException if {@code manager} is {@code null} + */ + public ConjugateGradient( + final IterationManager manager, final double delta, final boolean check) + throws NullArgumentException { + super(manager); + this.delta = delta; + this.check = check; + } + + /** + * Returns {@code true} if positive-definiteness should be checked for both matrix and + * preconditioner. + * + * @return {@code true} if the tests are to be performed + */ + public final boolean getCheck() { + return check; + } + + /** + * {@inheritDoc} + * + * @throws NonPositiveDefiniteOperatorException if {@code a} or {@code m} is not positive + * definite + */ + @Override + public RealVector solveInPlace( + final RealLinearOperator a, + final RealLinearOperator m, + final RealVector b, + final RealVector x0) + throws NullArgumentException, + NonPositiveDefiniteOperatorException, + NonSquareOperatorException, + DimensionMismatchException, + MaxCountExceededException { + checkParameters(a, m, b, x0); + final IterationManager manager = getIterationManager(); + // Initialization of default stopping criterion + manager.resetIterationCount(); + final double rmax = delta * b.getNorm(); + final RealVector bro = RealVector.unmodifiableRealVector(b); + + // Initialization phase counts as one iteration. + manager.incrementIterationCount(); + // p and x are constructed as copies of x0, since presumably, the type + // of x is optimized for the calculation of the matrix-vector product + // A.x. + final RealVector x = x0; + final RealVector xro = RealVector.unmodifiableRealVector(x); + final RealVector p = x.copy(); + RealVector q = a.operate(p); + + final RealVector r = b.combine(1, -1, q); + final RealVector rro = RealVector.unmodifiableRealVector(r); + double rnorm = r.getNorm(); + RealVector z; + if (m == null) { + z = r; + } else { + z = null; + } + IterativeLinearSolverEvent evt; + evt = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), xro, bro, rro, rnorm); + manager.fireInitializationEvent(evt); + if (rnorm <= rmax) { + manager.fireTerminationEvent(evt); + return x; + } + double rhoPrev = 0.; + while (true) { + manager.incrementIterationCount(); + evt = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), xro, bro, rro, rnorm); + manager.fireIterationStartedEvent(evt); + if (m != null) { + z = m.operate(r); + } + final double rhoNext = r.dotProduct(z); + if (check && (rhoNext <= 0.)) { + final NonPositiveDefiniteOperatorException e; + e = new NonPositiveDefiniteOperatorException(); + final ExceptionContext context = e.getContext(); + context.setValue(OPERATOR, m); + context.setValue(VECTOR, r); + throw e; + } + if (manager.getIterations() == 2) { + p.setSubVector(0, z); + } else { + p.combineToSelf(rhoNext / rhoPrev, 1., z); + } + q = a.operate(p); + final double pq = p.dotProduct(q); + if (check && (pq <= 0.)) { + final NonPositiveDefiniteOperatorException e; + e = new NonPositiveDefiniteOperatorException(); + final ExceptionContext context = e.getContext(); + context.setValue(OPERATOR, a); + context.setValue(VECTOR, p); + throw e; + } + final double alpha = rhoNext / pq; + x.combineToSelf(1., alpha, p); + r.combineToSelf(1., -alpha, q); + rhoPrev = rhoNext; + rnorm = r.getNorm(); + evt = + new DefaultIterativeLinearSolverEvent( + this, manager.getIterations(), xro, bro, rro, rnorm); + manager.fireIterationPerformedEvent(evt); + if (rnorm <= rmax) { + manager.fireTerminationEvent(evt); + return x; + } + } + } +} |