NonLinearConjugateGradientOptimizer.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.  *      http://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. package org.apache.commons.math4.legacy.optim.nonlinear.scalar.gradient;

  20. import org.apache.commons.math4.legacy.analysis.MultivariateFunction;
  21. import org.apache.commons.math4.legacy.exception.MathInternalError;
  22. import org.apache.commons.math4.legacy.exception.MathUnsupportedOperationException;
  23. import org.apache.commons.math4.legacy.exception.TooManyEvaluationsException;
  24. import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
  25. import org.apache.commons.math4.legacy.optim.ConvergenceChecker;
  26. import org.apache.commons.math4.legacy.optim.OptimizationData;
  27. import org.apache.commons.math4.legacy.optim.PointValuePair;
  28. import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GoalType;
  29. import org.apache.commons.math4.legacy.optim.nonlinear.scalar.GradientMultivariateOptimizer;

  31. /**
  32.  * Non-linear conjugate gradient optimizer.
  33.  * <br>
  34.  * This class supports both the Fletcher-Reeves and the Polak-Ribière
  35.  * update formulas for the conjugate search directions.
  36.  * It also supports optional preconditioning.
  37.  * <br>
  38.  * Line search must be setup via {@link org.apache.commons.math4.legacy.optim.nonlinear.scalar.LineSearchTolerance}.
  39.  * <br>
  40.  * Constraints are not supported: the call to
  41.  * {@link #optimize(OptimizationData[]) optimize} will throw
  42.  * {@link MathUnsupportedOperationException} if bounds are passed to it.
  43.  *
  44.  * @since 2.0
  45.  */
  46. public class NonLinearConjugateGradientOptimizer
  47.     extends GradientMultivariateOptimizer {
  48.     /** Update formula for the beta parameter. */
  49.     private final Formula updateFormula;
  50.     /** Preconditioner (may be null). */
  51.     private final Preconditioner preconditioner;

  53.     /**
  54.      * Available choices of update formulas for the updating the parameter
  55.      * that is used to compute the successive conjugate search directions.
  56.      * For non-linear conjugate gradients, there are
  57.      * two formulas:
  58.      * <ul>
  59.      *   <li>Fletcher-Reeves formula</li>
  60.      *   <li>Polak-Ribière formula</li>
  61.      * </ul>
  62.      *
  63.      * On the one hand, the Fletcher-Reeves formula is guaranteed to converge
  64.      * if the start point is close enough of the optimum whether the
  65.      * Polak-Ribière formula may not converge in rare cases. On the
  66.      * other hand, the Polak-Ribière formula is often faster when it
  67.      * does converge. Polak-Ribière is often used.
  68.      *
  69.      * @since 2.0
  70.      */
  71.     public enum Formula {
  72.         /** Fletcher-Reeves formula. */
  73.         FLETCHER_REEVES,
  74.         /** Polak-Ribière formula. */
  75.         POLAK_RIBIERE
  76.     }

  78.     /**
  79.      * Constructor with default {@link IdentityPreconditioner preconditioner}.
  80.      *
  81.      * @param updateFormula formula to use for updating the &beta; parameter,
  82.      * must be one of {@link Formula#FLETCHER_REEVES} or
  83.      * {@link Formula#POLAK_RIBIERE}.
  84.      * @param checker Convergence checker.
  85.      *
  86.      * @since 3.3
  87.      */
  88.     public NonLinearConjugateGradientOptimizer(final Formula updateFormula,
  89.                                                ConvergenceChecker<PointValuePair> checker) {
  90.         this(updateFormula,
  91.              checker,
  92.              new IdentityPreconditioner());
  93.     }

  95.     /**
  96.      * @param updateFormula formula to use for updating the &beta; parameter,
  97.      * must be one of {@link Formula#FLETCHER_REEVES} or
  98.      * {@link Formula#POLAK_RIBIERE}.
  99.      * @param checker Convergence checker.
  100.      * @param preconditioner Preconditioner.
  101.      *
  102.      * @since 3.3
  103.      */
  104.     public NonLinearConjugateGradientOptimizer(final Formula updateFormula,
  105.                                                ConvergenceChecker<PointValuePair> checker,
  106.                                                final Preconditioner preconditioner) {
  107.         super(checker);

  109.         this.updateFormula = updateFormula;
  110.         this.preconditioner = preconditioner;
  111.     }

  113.     /**
  114.      * {@inheritDoc}
  115.      */
  116.     @Override
  117.     public PointValuePair optimize(OptimizationData... optData)
  118.         throws TooManyEvaluationsException {
  119.         // Set up base class and perform computation.
  120.         return super.optimize(optData);
  121.     }

  123.     /** {@inheritDoc} */
  124.     @Override
  125.     protected PointValuePair doOptimize() {
  126.         final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
  127.         final double[] point = getStartPoint();
  128.         final GoalType goal = getGoalType();
  129.         final MultivariateFunction func = getObjectiveFunction();
  130.         final int n = point.length;
  131.         double[] r = computeObjectiveGradient(point);
  132.         if (goal == GoalType.MINIMIZE) {
  133.             for (int i = 0; i < n; i++) {
  134.                 r[i] = -r[i];
  135.             }
  136.         }

  138.         // Initial search direction.
  139.         double[] steepestDescent = preconditioner.precondition(point, r);
  140.         double[] searchDirection = steepestDescent.clone();

  142.         double delta = 0;
  143.         for (int i = 0; i < n; ++i) {
  144.             delta += r[i] * searchDirection[i];
  145.         }

  147.         createLineSearch();

  149.         PointValuePair current = null;
  150.         while (true) {
  151.             incrementIterationCount();

  153.             final double objective = func.value(point);
  154.             PointValuePair previous = current;
  155.             current = new PointValuePair(point, objective);
  156.             if (previous != null &&
  157.                 checker.converged(getIterations(), previous, current)) {
  158.                 // We have found an optimum.
  159.                 return current;
  160.             }

  162.             final double step = lineSearch(point, searchDirection).getPoint();

  164.             // Validate new point.
  165.             for (int i = 0; i < point.length; ++i) {
  166.                 point[i] += step * searchDirection[i];
  167.             }

  169.             r = computeObjectiveGradient(point);
  170.             if (goal == GoalType.MINIMIZE) {
  171.                 for (int i = 0; i < n; ++i) {
  172.                     r[i] = -r[i];
  173.                 }
  174.             }

  176.             // Compute beta.
  177.             final double deltaOld = delta;
  178.             final double[] newSteepestDescent = preconditioner.precondition(point, r);
  179.             delta = 0;
  180.             for (int i = 0; i < n; ++i) {
  181.                 delta += r[i] * newSteepestDescent[i];
  182.             }

  184.             final double beta;
  185.             switch (updateFormula) {
  186.             case FLETCHER_REEVES:
  187.                 beta = delta / deltaOld;
  188.                 break;
  189.             case POLAK_RIBIERE:
  190.                 double deltaMid = 0;
  191.                 for (int i = 0; i < r.length; ++i) {
  192.                     deltaMid += r[i] * steepestDescent[i];
  193.                 }
  194.                 beta = (delta - deltaMid) / deltaOld;
  195.                 break;
  196.             default:
  197.                 // Should never happen.
  198.                 throw new MathInternalError();
  199.             }
  200.             steepestDescent = newSteepestDescent;

  202.             // Compute conjugate search direction.
  203.             if (getIterations() % n == 0 ||
  204.                 beta < 0) {
  205.                 // Break conjugation: reset search direction.
  206.                 searchDirection = steepestDescent.clone();
  207.             } else {
  208.                 // Compute new conjugate search direction.
  209.                 for (int i = 0; i < n; ++i) {
  210.                     searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
  211.                 }
  212.             }
  213.         }
  214.     }

  216.     /**
  217.      * {@inheritDoc}
  218.      */
  219.     @Override
  220.     protected void parseOptimizationData(OptimizationData... optData) {
  221.         // Allow base class to register its own data.
  222.         super.parseOptimizationData(optData);

  224.         checkParameters();
  225.     }

  227.     /** Default identity preconditioner. */
  228.     public static class IdentityPreconditioner implements Preconditioner {
  229.         /** {@inheritDoc} */
  230.         @Override
  231.         public double[] precondition(double[] variables, double[] r) {
  232.             return r.clone();
  233.         }
  234.     }

  236.     // Class is not used anymore (cf. MATH-1092). However, it might
  237.     // be interesting to create a class similar to "LineSearch", but
  238.     // that will take advantage that the model's gradient is available.
  239. //     /**
  240. //      * Internal class for line search.
  241. //      * <p>
  242. //      * The function represented by this class is the dot product of
  243. //      * the objective function gradient and the search direction. Its
  244. //      * value is zero when the gradient is orthogonal to the search
  245. //      * direction, i.e. when the objective function value is a local
  246. //      * extremum along the search direction.
  247. //      * </p>
  248. //      */
  249. //     private class LineSearchFunction implements UnivariateFunction {
  250. //         /** Current point. */
  251. //         private final double[] currentPoint;
  252. //         /** Search direction. */
  253. //         private final double[] searchDirection;

  255. //         /**
  256. //          * @param point Current point.
  257. //          * @param direction Search direction.
  258. //          */
  259. //         public LineSearchFunction(double[] point,
  260. //                                   double[] direction) {
  261. //             currentPoint = point.clone();
  262. //             searchDirection = direction.clone();
  263. //         }

  265. //         /** {@inheritDoc} */
  266. //         public double value(double x) {
  267. //             // current point in the search direction
  268. //             final double[] shiftedPoint = currentPoint.clone();
  269. //             for (int i = 0; i < shiftedPoint.length; ++i) {
  270. //                 shiftedPoint[i] += x * searchDirection[i];
  271. //             }

  273. //             // gradient of the objective function
  274. //             final double[] gradient = computeObjectiveGradient(shiftedPoint);

  276. //             // dot product with the search direction
  277. //             double dotProduct = 0;
  278. //             for (int i = 0; i < gradient.length; ++i) {
  279. //                 dotProduct += gradient[i] * searchDirection[i];
  280. //             }

  282. //             return dotProduct;
  283. //         }
  284. //     }

  286.     /**
  287.      * @throws MathUnsupportedOperationException if bounds were passed to the
  288.      * {@link #optimize(OptimizationData[]) optimize} method.
  289.      */
  290.     private void checkParameters() {
  291.         if (getLowerBound() != null ||
  292.             getUpperBound() != null) {
  293.             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
  294.         }
  295.     }
  296. }
Created with JaCoCo 0.8.10.202304240956Code Coverage Report for Apache Commons Math 4.0-SNAPSHOT