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.  *      https://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. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.optim.nonlinear.scalar.noderiv;

  24. import org.hipparchus.exception.LocalizedCoreFormats;
  25. import org.hipparchus.exception.MathIllegalArgumentException;
  26. import org.hipparchus.exception.MathRuntimeException;
  27. import org.hipparchus.optim.ConvergenceChecker;
  28. import org.hipparchus.optim.OptimizationData;
  29. import org.hipparchus.optim.PointValuePair;
  30. import org.hipparchus.optim.nonlinear.scalar.GoalType;
  31. import org.hipparchus.optim.nonlinear.scalar.LineSearch;
  32. import org.hipparchus.optim.nonlinear.scalar.MultivariateOptimizer;
  33. import org.hipparchus.optim.univariate.UnivariatePointValuePair;
  34. import org.hipparchus.util.FastMath;

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

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

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

  119.         if (rel < MIN_RELATIVE_TOLERANCE) {
  120.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
  121.                                                    rel, MIN_RELATIVE_TOLERANCE);
  122.         }
  123.         if (abs <= 0) {
  124.             throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL_BOUND_EXCLUDED,
  125.                                                    abs, 0);
  126.         }
  127.         relativeThreshold = rel;
  128.         absoluteThreshold = abs;

  130.         // Create the line search optimizer.
  131.         line = new LineSearch(this,
  132.                               lineRel,
  133.                               lineAbs,
  134.                               1d);
  135.     }

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

  153.     /**
  154.      * Builds an instance with the default convergence checking procedure.
  155.      *
  156.      * @param rel Relative threshold.
  157.      * @param abs Absolute threshold.
  158.      * @param lineRel Relative threshold for the internal line search optimizer.
  159.      * @param lineAbs Absolute threshold for the internal line search optimizer.
  160.      * @throws MathIllegalArgumentException if {@code abs <= 0}.
  161.      * @throws MathIllegalArgumentException if {@code rel < 2 * FastMath.ulp(1d)}.
  162.      */
  163.     public PowellOptimizer(double rel,
  164.                            double abs,
  165.                            double lineRel,
  166.                            double lineAbs) {
  167.         this(rel, abs, lineRel, lineAbs, null);
  168.     }

  170.     /** {@inheritDoc} */
  171.     @Override
  172.     protected PointValuePair doOptimize() {
  173.         checkParameters();

  175.         final GoalType goal = getGoalType();
  176.         final double[] guess = getStartPoint();
  177.         final int n = guess.length;

  179.         final double[][] direc = new double[n][n];
  180.         for (int i = 0; i < n; i++) {
  181.             direc[i][i] = 1;
  182.         }

  184.         final ConvergenceChecker<PointValuePair> checker
  185.             = getConvergenceChecker();

  187.         double[] x = guess;
  188.         double fVal = computeObjectiveValue(x);
  189.         double[] x1 = x.clone();
  190.         while (true) {
  191.             incrementIterationCount();

  193.             double fX = fVal;
  194.             double delta = 0;
  195.             int bigInd = 0;

  197.             for (int i = 0; i < n; i++) {
  198.                 final double[] d = direc[i].clone();

  200.                 final double fX2 = fVal;

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

  208.                 if ((fX2 - fVal) > delta) {
  209.                     delta = fX2 - fVal;
  210.                     bigInd = i;
  211.                 }
  212.             }

  214.             // Default convergence check.
  215.             boolean stop = 2 * (fX - fVal) <=
  216.                 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
  217.                  absoluteThreshold);

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

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

  239.             x1 = x.clone();
  240.             final double fX2 = computeObjectiveValue(x2);

  242.             if (fX > fX2) {
  243.                 double t = 2 * (fX + fX2 - 2 * fVal);
  244.                 double temp = fX - fVal - delta;
  245.                 t *= temp * temp;
  246.                 temp = fX - fX2;
  247.                 t -= delta * temp * temp;

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

  256.                     final int lastInd = n - 1;
  257.                     direc[bigInd] = direc[lastInd];
  258.                     direc[lastInd] = result[1];
  259.                 }
  260.             }
  261.         }
  262.     }

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

  285.         final double[][] result = new double[2][];
  286.         result[0] = nP;
  287.         result[1] = nD;

  289.         return result;
  290.     }

  292.     /**
  293.      * @throws MathRuntimeException if bounds were passed to the
  294.      * {@link #optimize(OptimizationData...)}  optimize} method.
  295.      */
  296.     private void checkParameters() {
  297.         if (getLowerBound() != null ||
  298.             getUpperBound() != null) {
  299.             throw new MathRuntimeException(LocalizedCoreFormats.CONSTRAINT);
  300.         }
  301.     }
  302. }
Created with JaCoCo 0.8.12.202403310830