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.  */
  17. package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;

  19. import org.apache.commons.math3.util.FastMath;
  20. import org.apache.commons.math3.util.MathArrays;
  21. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  22. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  23. import org.apache.commons.math3.exception.MathUnsupportedOperationException;
  24. import org.apache.commons.math3.exception.util.LocalizedFormats;
  25. import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
  26. import org.apache.commons.math3.optim.PointValuePair;
  27. import org.apache.commons.math3.optim.ConvergenceChecker;
  28. import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
  29. import org.apache.commons.math3.optim.nonlinear.scalar.LineSearch;
  30. import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;

  32. /**
  33.  * Powell's algorithm.
  34.  * This code is translated and adapted from the Python version of this
  35.  * algorithm (as implemented in module {@code optimize.py} v0.5 of
  36.  * <em>SciPy</em>).
  37.  * <br/>
  38.  * The default stopping criterion is based on the differences of the
  39.  * function value between two successive iterations. It is however possible
  40.  * to define a custom convergence checker that might terminate the algorithm
  41.  * earlier.
  42.  * <br/>
  43.  * Line search is performed by the {@link LineSearch} class.
  44.  * <br/>
  45.  * Constraints are not supported: the call to
  46.  * {@link #optimize(OptimizationData[]) optimize} will throw
  47.  * {@link MathUnsupportedOperationException} if bounds are passed to it.
  48.  * In order to impose simple constraints, the objective function must be
  49.  * wrapped in an adapter like
  50.  * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
  51.  * MultivariateFunctionMappingAdapter} or
  52.  * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
  53.  * MultivariateFunctionPenaltyAdapter}.
  54.  *
  55.  * @since 2.2
  56.  */
  57. public class PowellOptimizer
  58.     extends MultivariateOptimizer {
  59.     /**
  60.      * Minimum relative tolerance.
  61.      */
  62.     private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
  63.     /**
  64.      * Relative threshold.
  65.      */
  66.     private final double relativeThreshold;
  67.     /**
  68.      * Absolute threshold.
  69.      */
  70.     private final double absoluteThreshold;
  71.     /**
  72.      * Line search.
  73.      */
  74.     private final LineSearch line;

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

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

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

  125.         // Create the line search optimizer.
  126.         line = new LineSearch(this,
  127.                               lineRel,
  128.                               lineAbs,
  129.                               1d);
  130.     }

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

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

  165.     /** {@inheritDoc} */
  166.     @Override
  167.     protected PointValuePair doOptimize() {
  168.         checkParameters();

  170.         final GoalType goal = getGoalType();
  171.         final double[] guess = getStartPoint();
  172.         final int n = guess.length;

  174.         final double[][] direc = new double[n][n];
  175.         for (int i = 0; i < n; i++) {
  176.             direc[i][i] = 1;
  177.         }

  179.         final ConvergenceChecker<PointValuePair> checker
  180.             = getConvergenceChecker();

  182.         double[] x = guess;
  183.         double fVal = computeObjectiveValue(x);
  184.         double[] x1 = x.clone();
  185.         while (true) {
  186.             incrementIterationCount();

  188.             double fX = fVal;
  189.             double fX2 = 0;
  190.             double delta = 0;
  191.             int bigInd = 0;
  192.             double alphaMin = 0;

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

  197.                 fX2 = fVal;

  199.                 final UnivariatePointValuePair optimum = line.search(x, d);
  200.                 fVal = optimum.getValue();
  201.                 alphaMin = optimum.getPoint();
  202.                 final double[][] result = newPointAndDirection(x, d, alphaMin);
  203.                 x = result[0];

  205.                 if ((fX2 - fVal) > delta) {
  206.                     delta = fX2 - fVal;
  207.                     bigInd = i;
  208.                 }
  209.             }

  211.             // Default convergence check.
  212.             boolean stop = 2 * (fX - fVal) <=
  213.                 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
  214.                  absoluteThreshold);

  216.             final PointValuePair previous = new PointValuePair(x1, fX);
  217.             final PointValuePair current = new PointValuePair(x, fVal);
  218.             if (!stop && checker != null) { // User-defined stopping criteria.
  219.                 stop = checker.converged(getIterations(), previous, current);
  220.             }
  221.             if (stop) {
  222.                 if (goal == GoalType.MINIMIZE) {
  223.                     return (fVal < fX) ? current : previous;
  224.                 } else {
  225.                     return (fVal > fX) ? current : previous;
  226.                 }
  227.             }

  229.             final double[] d = new double[n];
  230.             final double[] x2 = new double[n];
  231.             for (int i = 0; i < n; i++) {
  232.                 d[i] = x[i] - x1[i];
  233.                 x2[i] = 2 * x[i] - x1[i];
  234.             }

  236.             x1 = x.clone();
  237.             fX2 = computeObjectiveValue(x2);

  239.             if (fX > fX2) {
  240.                 double t = 2 * (fX + fX2 - 2 * fVal);
  241.                 double temp = fX - fVal - delta;
  242.                 t *= temp * temp;
  243.                 temp = fX - fX2;
  244.                 t -= delta * temp * temp;

  246.                 if (t < 0.0) {
  247.                     final UnivariatePointValuePair optimum = line.search(x, d);
  248.                     fVal = optimum.getValue();
  249.                     alphaMin = optimum.getPoint();
  250.                     final double[][] result = newPointAndDirection(x, d, alphaMin);
  251.                     x = result[0];

  253.                     final int lastInd = n - 1;
  254.                     direc[bigInd] = direc[lastInd];
  255.                     direc[lastInd] = result[1];
  256.                 }
  257.             }
  258.         }
  259.     }

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

  282.         final double[][] result = new double[2][];
  283.         result[0] = nP;
  284.         result[1] = nD;

  286.         return result;
  287.     }

  289.     /**
  290.      * @throws MathUnsupportedOperationException if bounds were passed to the
  291.      * {@link #optimize(OptimizationData[]) optimize} method.
  292.      */
  293.     private void checkParameters() {
  294.         if (getLowerBound() != null ||
  295.             getUpperBound() != null) {
  296.             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
  297.         }
  298.     }
  299. }
Created with JaCoCo 0.7.2.201409121644