PowellOptimizer.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.direct;

  20. import org.apache.commons.math3.util.FastMath;
  21. import org.apache.commons.math3.util.MathArrays;
  22. import org.apache.commons.math3.analysis.UnivariateFunction;
  23. import org.apache.commons.math3.analysis.MultivariateFunction;
  24. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  25. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  26. import org.apache.commons.math3.optimization.GoalType;
  27. import org.apache.commons.math3.optimization.PointValuePair;
  28. import org.apache.commons.math3.optimization.ConvergenceChecker;
  29. import org.apache.commons.math3.optimization.MultivariateOptimizer;
  30. import org.apache.commons.math3.optimization.univariate.BracketFinder;
  31. import org.apache.commons.math3.optimization.univariate.BrentOptimizer;
  32. import org.apache.commons.math3.optimization.univariate.UnivariatePointValuePair;
  33. import org.apache.commons.math3.optimization.univariate.SimpleUnivariateValueChecker;

  35. /**
  36.  * Powell algorithm.
  37.  * This code is translated and adapted from the Python version of this
  38.  * algorithm (as implemented in module {@code optimize.py} v0.5 of
  39.  * <em>SciPy</em>).
  40.  * <br/>
  41.  * The default stopping criterion is based on the differences of the
  42.  * function value between two successive iterations. It is however possible
  43.  * to define a custom convergence checker that might terminate the algorithm
  44.  * earlier.
  45.  * <br/>
  46.  * The internal line search optimizer is a {@link BrentOptimizer} with a
  47.  * convergence checker set to {@link SimpleUnivariateValueChecker}.
  48.  *
  49.  * @deprecated As of 3.1 (to be removed in 4.0).
  50.  * @since 2.2
  51.  */
  52. @Deprecated
  53. public class PowellOptimizer
  54.     extends BaseAbstractMultivariateOptimizer<MultivariateFunction>
  55.     implements MultivariateOptimizer {
  56.     /**
  57.      * Minimum relative tolerance.
  58.      */
  59.     private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
  60.     /**
  61.      * Relative threshold.
  62.      */
  63.     private final double relativeThreshold;
  64.     /**
  65.      * Absolute threshold.
  66.      */
  67.     private final double absoluteThreshold;
  68.     /**
  69.      * Line search.
  70.      */
  71.     private final LineSearch line;

  73.     /**
  74.      * This constructor allows to specify a user-defined convergence checker,
  75.      * in addition to the parameters that control the default convergence
  76.      * checking procedure.
  77.      * <br/>
  78.      * The internal line search tolerances are set to the square-root of their
  79.      * corresponding value in the multivariate optimizer.
  80.      *
  81.      * @param rel Relative threshold.
  82.      * @param abs Absolute threshold.
  83.      * @param checker Convergence checker.
  84.      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
  85.      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
  86.      */
  87.     public PowellOptimizer(double rel,
  88.                            double abs,
  89.                            ConvergenceChecker<PointValuePair> checker) {
  90.         this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
  91.     }

  93.     /**
  94.      * This constructor allows to specify a user-defined convergence checker,
  95.      * in addition to the parameters that control the default convergence
  96.      * checking procedure and the line search tolerances.
  97.      *
  98.      * @param rel Relative threshold for this optimizer.
  99.      * @param abs Absolute threshold for this optimizer.
  100.      * @param lineRel Relative threshold for the internal line search optimizer.
  101.      * @param lineAbs Absolute threshold for the internal line search optimizer.
  102.      * @param checker Convergence checker.
  103.      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
  104.      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
  105.      */
  106.     public PowellOptimizer(double rel,
  107.                            double abs,
  108.                            double lineRel,
  109.                            double lineAbs,
  110.                            ConvergenceChecker<PointValuePair> checker) {
  111.         super(checker);

  113.         if (rel < MIN_RELATIVE_TOLERANCE) {
  114.             throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
  115.         }
  116.         if (abs <= 0) {
  117.             throw new NotStrictlyPositiveException(abs);
  118.         }
  119.         relativeThreshold = rel;
  120.         absoluteThreshold = abs;

  122.         // Create the line search optimizer.
  123.         line = new LineSearch(lineRel,
  124.                               lineAbs);
  125.     }

  127.     /**
  128.      * The parameters control the default convergence checking procedure.
  129.      * <br/>
  130.      * The internal line search tolerances are set to the square-root of their
  131.      * corresponding value in the multivariate optimizer.
  132.      *
  133.      * @param rel Relative threshold.
  134.      * @param abs Absolute threshold.
  135.      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
  136.      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
  137.      */
  138.     public PowellOptimizer(double rel,
  139.                            double abs) {
  140.         this(rel, abs, null);
  141.     }

  143.     /**
  144.      * Builds an instance with the default convergence checking procedure.
  145.      *
  146.      * @param rel Relative threshold.
  147.      * @param abs Absolute threshold.
  148.      * @param lineRel Relative threshold for the internal line search optimizer.
  149.      * @param lineAbs Absolute threshold for the internal line search optimizer.
  150.      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
  151.      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
  152.      * @since 3.1
  153.      */
  154.     public PowellOptimizer(double rel,
  155.                            double abs,
  156.                            double lineRel,
  157.                            double lineAbs) {
  158.         this(rel, abs, lineRel, lineAbs, null);
  159.     }

  161.     /** {@inheritDoc} */
  162.     @Override
  163.     protected PointValuePair doOptimize() {
  164.         final GoalType goal = getGoalType();
  165.         final double[] guess = getStartPoint();
  166.         final int n = guess.length;

  168.         final double[][] direc = new double[n][n];
  169.         for (int i = 0; i < n; i++) {
  170.             direc[i][i] = 1;
  171.         }

  173.         final ConvergenceChecker<PointValuePair> checker
  174.             = getConvergenceChecker();

  176.         double[] x = guess;
  177.         double fVal = computeObjectiveValue(x);
  178.         double[] x1 = x.clone();
  179.         int iter = 0;
  180.         while (true) {
  181.             ++iter;

  183.             double fX = fVal;
  184.             double fX2 = 0;
  185.             double delta = 0;
  186.             int bigInd = 0;
  187.             double alphaMin = 0;

  189.             for (int i = 0; i < n; i++) {
  190.                 final double[] d = MathArrays.copyOf(direc[i]);

  192.                 fX2 = fVal;

  194.                 final UnivariatePointValuePair optimum = line.search(x, d);
  195.                 fVal = optimum.getValue();
  196.                 alphaMin = optimum.getPoint();
  197.                 final double[][] result = newPointAndDirection(x, d, alphaMin);
  198.                 x = result[0];

  200.                 if ((fX2 - fVal) > delta) {
  201.                     delta = fX2 - fVal;
  202.                     bigInd = i;
  203.                 }
  204.             }

  206.             // Default convergence check.
  207.             boolean stop = 2 * (fX - fVal) <=
  208.                 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
  209.                  absoluteThreshold);

  211.             final PointValuePair previous = new PointValuePair(x1, fX);
  212.             final PointValuePair current = new PointValuePair(x, fVal);
  213.             if (!stop && checker != null) {
  214.                 stop = checker.converged(iter, previous, current);
  215.             }
  216.             if (stop) {
  217.                 if (goal == GoalType.MINIMIZE) {
  218.                     return (fVal < fX) ? current : previous;
  219.                 } else {
  220.                     return (fVal > fX) ? current : previous;
  221.                 }
  222.             }

  224.             final double[] d = new double[n];
  225.             final double[] x2 = new double[n];
  226.             for (int i = 0; i < n; i++) {
  227.                 d[i] = x[i] - x1[i];
  228.                 x2[i] = 2 * x[i] - x1[i];
  229.             }

  231.             x1 = x.clone();
  232.             fX2 = computeObjectiveValue(x2);

  234.             if (fX > fX2) {
  235.                 double t = 2 * (fX + fX2 - 2 * fVal);
  236.                 double temp = fX - fVal - delta;
  237.                 t *= temp * temp;
  238.                 temp = fX - fX2;
  239.                 t -= delta * temp * temp;

  241.                 if (t < 0.0) {
  242.                     final UnivariatePointValuePair optimum = line.search(x, d);
  243.                     fVal = optimum.getValue();
  244.                     alphaMin = optimum.getPoint();
  245.                     final double[][] result = newPointAndDirection(x, d, alphaMin);
  246.                     x = result[0];

  248.                     final int lastInd = n - 1;
  249.                     direc[bigInd] = direc[lastInd];
  250.                     direc[lastInd] = result[1];
  251.                 }
  252.             }
  253.         }
  254.     }

  256.     /**
  257.      * Compute a new point (in the original space) and a new direction
  258.      * vector, resulting from the line search.
  259.      *
  260.      * @param p Point used in the line search.
  261.      * @param d Direction used in the line search.
  262.      * @param optimum Optimum found by the line search.
  263.      * @return a 2-element array containing the new point (at index 0) and
  264.      * the new direction (at index 1).
  265.      */
  266.     private double[][] newPointAndDirection(double[] p,
  267.                                             double[] d,
  268.                                             double optimum) {
  269.         final int n = p.length;
  270.         final double[] nP = new double[n];
  271.         final double[] nD = new double[n];
  272.         for (int i = 0; i < n; i++) {
  273.             nD[i] = d[i] * optimum;
  274.             nP[i] = p[i] + nD[i];
  275.         }

  277.         final double[][] result = new double[2][];
  278.         result[0] = nP;
  279.         result[1] = nD;

  281.         return result;
  282.     }

  284.     /**
  285.      * Class for finding the minimum of the objective function along a given
  286.      * direction.
  287.      */
  288.     private class LineSearch extends BrentOptimizer {
  289.         /**
  290.          * Value that will pass the precondition check for {@link BrentOptimizer}
  291.          * but will not pass the convergence check, so that the custom checker
  292.          * will always decide when to stop the line search.
  293.          */
  294.         private static final double REL_TOL_UNUSED = 1e-15;
  295.         /**
  296.          * Value that will pass the precondition check for {@link BrentOptimizer}
  297.          * but will not pass the convergence check, so that the custom checker
  298.          * will always decide when to stop the line search.
  299.          */
  300.         private static final double ABS_TOL_UNUSED = Double.MIN_VALUE;
  301.         /**
  302.          * Automatic bracketing.
  303.          */
  304.         private final BracketFinder bracket = new BracketFinder();

  306.         /**
  307.          * The "BrentOptimizer" default stopping criterion uses the tolerances
  308.          * to check the domain (point) values, not the function values.
  309.          * We thus create a custom checker to use function values.
  310.          *
  311.          * @param rel Relative threshold.
  312.          * @param abs Absolute threshold.
  313.          */
  314.         LineSearch(double rel,
  315.                    double abs) {
  316.             super(REL_TOL_UNUSED,
  317.                   ABS_TOL_UNUSED,
  318.                   new SimpleUnivariateValueChecker(rel, abs));
  319.         }

  321.         /**
  322.          * Find the minimum of the function {@code f(p + alpha * d)}.
  323.          *
  324.          * @param p Starting point.
  325.          * @param d Search direction.
  326.          * @return the optimum.
  327.          * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
  328.          * if the number of evaluations is exceeded.
  329.          */
  330.         public UnivariatePointValuePair search(final double[] p, final double[] d) {
  331.             final int n = p.length;
  332.             final UnivariateFunction f = new UnivariateFunction() {
  333.                     public double value(double alpha) {
  334.                         final double[] x = new double[n];
  335.                         for (int i = 0; i < n; i++) {
  336.                             x[i] = p[i] + alpha * d[i];
  337.                         }
  338.                         final double obj = PowellOptimizer.this.computeObjectiveValue(x);
  339.                         return obj;
  340.                     }
  341.                 };

  343.             final GoalType goal = PowellOptimizer.this.getGoalType();
  344.             bracket.search(f, goal, 0, 1);
  345.             // Passing "MAX_VALUE" as a dummy value because it is the enclosing
  346.             // class that counts the number of evaluations (and will eventually
  347.             // generate the exception).
  348.             return optimize(Integer.MAX_VALUE, f, goal,
  349.                             bracket.getLo(), bracket.getHi(), bracket.getMid());
  350.         }
  351.     }
  352. }
Created with JaCoCo 0.7.2.201409121644