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.math3.optimization.general;

  20. import org.apache.commons.math3.exception.MathIllegalStateException;
  21. import org.apache.commons.math3.analysis.UnivariateFunction;
  22. import org.apache.commons.math3.analysis.solvers.BrentSolver;
  23. import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
  24. import org.apache.commons.math3.exception.util.LocalizedFormats;
  25. import org.apache.commons.math3.optimization.GoalType;
  26. import org.apache.commons.math3.optimization.PointValuePair;
  27. import org.apache.commons.math3.optimization.SimpleValueChecker;
  28. import org.apache.commons.math3.optimization.ConvergenceChecker;
  29. import org.apache.commons.math3.util.FastMath;

  31. /**
  32.  * Non-linear conjugate gradient optimizer.
  33.  * <p>
  34.  * This class supports both the Fletcher-Reeves and the Polak-Ribi&egrave;re
  35.  * update formulas for the conjugate search directions. It also supports
  36.  * optional preconditioning.
  37.  * </p>
  38.  *
  39.  * @deprecated As of 3.1 (to be removed in 4.0).
  40.  * @since 2.0
  41.  *
  42.  */
  43. @Deprecated
  44. public class NonLinearConjugateGradientOptimizer
  45.     extends AbstractScalarDifferentiableOptimizer {
  46.     /** Update formula for the beta parameter. */
  47.     private final ConjugateGradientFormula updateFormula;
  48.     /** Preconditioner (may be null). */
  49.     private final Preconditioner preconditioner;
  50.     /** solver to use in the line search (may be null). */
  51.     private final UnivariateSolver solver;
  52.     /** Initial step used to bracket the optimum in line search. */
  53.     private double initialStep;
  54.     /** Current point. */
  55.     private double[] point;

  57.     /**
  58.      * Constructor with default {@link SimpleValueChecker checker},
  59.      * {@link BrentSolver line search solver} and
  60.      * {@link IdentityPreconditioner preconditioner}.
  61.      *
  62.      * @param updateFormula formula to use for updating the &beta; parameter,
  63.      * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
  64.      * ConjugateGradientFormula#POLAK_RIBIERE}.
  65.      * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()}
  66.      */
  67.     @Deprecated
  68.     public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) {
  69.         this(updateFormula,
  70.              new SimpleValueChecker());
  71.     }

  73.     /**
  74.      * Constructor with default {@link BrentSolver line search solver} and
  75.      * {@link IdentityPreconditioner preconditioner}.
  76.      *
  77.      * @param updateFormula formula to use for updating the &beta; parameter,
  78.      * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
  79.      * ConjugateGradientFormula#POLAK_RIBIERE}.
  80.      * @param checker Convergence checker.
  81.      */
  82.     public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
  83.                                                ConvergenceChecker<PointValuePair> checker) {
  84.         this(updateFormula,
  85.              checker,
  86.              new BrentSolver(),
  87.              new IdentityPreconditioner());
  88.     }


  91.     /**
  92.      * Constructor with default {@link IdentityPreconditioner preconditioner}.
  93.      *
  94.      * @param updateFormula formula to use for updating the &beta; parameter,
  95.      * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
  96.      * ConjugateGradientFormula#POLAK_RIBIERE}.
  97.      * @param checker Convergence checker.
  98.      * @param lineSearchSolver Solver to use during line search.
  99.      */
  100.     public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
  101.                                                ConvergenceChecker<PointValuePair> checker,
  102.                                                final UnivariateSolver lineSearchSolver) {
  103.         this(updateFormula,
  104.              checker,
  105.              lineSearchSolver,
  106.              new IdentityPreconditioner());
  107.     }

  109.     /**
  110.      * @param updateFormula formula to use for updating the &beta; parameter,
  111.      * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
  112.      * ConjugateGradientFormula#POLAK_RIBIERE}.
  113.      * @param checker Convergence checker.
  114.      * @param lineSearchSolver Solver to use during line search.
  115.      * @param preconditioner Preconditioner.
  116.      */
  117.     public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
  118.                                                ConvergenceChecker<PointValuePair> checker,
  119.                                                final UnivariateSolver lineSearchSolver,
  120.                                                final Preconditioner preconditioner) {
  121.         super(checker);

  123.         this.updateFormula = updateFormula;
  124.         solver = lineSearchSolver;
  125.         this.preconditioner = preconditioner;
  126.         initialStep = 1.0;
  127.     }

  129.     /**
  130.      * Set the initial step used to bracket the optimum in line search.
  131.      * <p>
  132.      * The initial step is a factor with respect to the search direction,
  133.      * which itself is roughly related to the gradient of the function
  134.      * </p>
  135.      * @param initialStep initial step used to bracket the optimum in line search,
  136.      * if a non-positive value is used, the initial step is reset to its
  137.      * default value of 1.0
  138.      */
  139.     public void setInitialStep(final double initialStep) {
  140.         if (initialStep <= 0) {
  141.             this.initialStep = 1.0;
  142.         } else {
  143.             this.initialStep = initialStep;
  144.         }
  145.     }

  147.     /** {@inheritDoc} */
  148.     @Override
  149.     protected PointValuePair doOptimize() {
  150.         final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
  151.         point = getStartPoint();
  152.         final GoalType goal = getGoalType();
  153.         final int n = point.length;
  154.         double[] r = computeObjectiveGradient(point);
  155.         if (goal == GoalType.MINIMIZE) {
  156.             for (int i = 0; i < n; ++i) {
  157.                 r[i] = -r[i];
  158.             }
  159.         }

  161.         // Initial search direction.
  162.         double[] steepestDescent = preconditioner.precondition(point, r);
  163.         double[] searchDirection = steepestDescent.clone();

  165.         double delta = 0;
  166.         for (int i = 0; i < n; ++i) {
  167.             delta += r[i] * searchDirection[i];
  168.         }

  170.         PointValuePair current = null;
  171.         int iter = 0;
  172.         int maxEval = getMaxEvaluations();
  173.         while (true) {
  174.             ++iter;

  176.             final double objective = computeObjectiveValue(point);
  177.             PointValuePair previous = current;
  178.             current = new PointValuePair(point, objective);
  179.             if (previous != null && checker.converged(iter, previous, current)) {
  180.                 // We have found an optimum.
  181.                 return current;
  182.             }

  184.             // Find the optimal step in the search direction.
  185.             final UnivariateFunction lsf = new LineSearchFunction(searchDirection);
  186.             final double uB = findUpperBound(lsf, 0, initialStep);
  187.             // XXX Last parameters is set to a value close to zero in order to
  188.             // work around the divergence problem in the "testCircleFitting"
  189.             // unit test (see MATH-439).
  190.             final double step = solver.solve(maxEval, lsf, 0, uB, 1e-15);
  191.             maxEval -= solver.getEvaluations(); // Subtract used up evaluations.

  193.             // Validate new point.
  194.             for (int i = 0; i < point.length; ++i) {
  195.                 point[i] += step * searchDirection[i];
  196.             }

  198.             r = computeObjectiveGradient(point);
  199.             if (goal == GoalType.MINIMIZE) {
  200.                 for (int i = 0; i < n; ++i) {
  201.                     r[i] = -r[i];
  202.                 }
  203.             }

  205.             // Compute beta.
  206.             final double deltaOld = delta;
  207.             final double[] newSteepestDescent = preconditioner.precondition(point, r);
  208.             delta = 0;
  209.             for (int i = 0; i < n; ++i) {
  210.                 delta += r[i] * newSteepestDescent[i];
  211.             }

  213.             final double beta;
  214.             if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) {
  215.                 beta = delta / deltaOld;
  216.             } else {
  217.                 double deltaMid = 0;
  218.                 for (int i = 0; i < r.length; ++i) {
  219.                     deltaMid += r[i] * steepestDescent[i];
  220.                 }
  221.                 beta = (delta - deltaMid) / deltaOld;
  222.             }
  223.             steepestDescent = newSteepestDescent;

  225.             // Compute conjugate search direction.
  226.             if (iter % n == 0 ||
  227.                 beta < 0) {
  228.                 // Break conjugation: reset search direction.
  229.                 searchDirection = steepestDescent.clone();
  230.             } else {
  231.                 // Compute new conjugate search direction.
  232.                 for (int i = 0; i < n; ++i) {
  233.                     searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
  234.                 }
  235.             }
  236.         }
  237.     }

  239.     /**
  240.      * Find the upper bound b ensuring bracketing of a root between a and b.
  241.      *
  242.      * @param f function whose root must be bracketed.
  243.      * @param a lower bound of the interval.
  244.      * @param h initial step to try.
  245.      * @return b such that f(a) and f(b) have opposite signs.
  246.      * @throws MathIllegalStateException if no bracket can be found.
  247.      */
  248.     private double findUpperBound(final UnivariateFunction f,
  249.                                   final double a, final double h) {
  250.         final double yA = f.value(a);
  251.         double yB = yA;
  252.         for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) {
  253.             final double b = a + step;
  254.             yB = f.value(b);
  255.             if (yA * yB <= 0) {
  256.                 return b;
  257.             }
  258.         }
  259.         throw new MathIllegalStateException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH);
  260.     }

  262.     /** Default identity preconditioner. */
  263.     public static class IdentityPreconditioner implements Preconditioner {

  265.         /** {@inheritDoc} */
  266.         public double[] precondition(double[] variables, double[] r) {
  267.             return r.clone();
  268.         }
  269.     }

  271.     /** Internal class for line search.
  272.      * <p>
  273.      * The function represented by this class is the dot product of
  274.      * the objective function gradient and the search direction. Its
  275.      * value is zero when the gradient is orthogonal to the search
  276.      * direction, i.e. when the objective function value is a local
  277.      * extremum along the search direction.
  278.      * </p>
  279.      */
  280.     private class LineSearchFunction implements UnivariateFunction {
  281.         /** Search direction. */
  282.         private final double[] searchDirection;

  284.         /** Simple constructor.
  285.          * @param searchDirection search direction
  286.          */
  287.         public LineSearchFunction(final double[] searchDirection) {
  288.             this.searchDirection = searchDirection;
  289.         }

  291.         /** {@inheritDoc} */
  292.         public double value(double x) {
  293.             // current point in the search direction
  294.             final double[] shiftedPoint = point.clone();
  295.             for (int i = 0; i < shiftedPoint.length; ++i) {
  296.                 shiftedPoint[i] += x * searchDirection[i];
  297.             }

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

  302.             // dot product with the search direction
  303.             double dotProduct = 0;
  304.             for (int i = 0; i < gradient.length; ++i) {
  305.                 dotProduct += gradient[i] * searchDirection[i];
  306.             }

  308.             return dotProduct;
  309.         }
  310.     }
  311. }
Created with JaCoCo 0.7.2.201409121644