    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.optim.nonlinear.scalar.noderiv;
    018
    019    import org.apache.commons.math3.util.FastMath;
    020    import org.apache.commons.math3.util.MathArrays;
    021    import org.apache.commons.math3.analysis.UnivariateFunction;
    022    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    023    import org.apache.commons.math3.exception.NotStrictlyPositiveException;
    024    import org.apache.commons.math3.optim.MaxEval;
    025    import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
    026    import org.apache.commons.math3.optim.PointValuePair;
    027    import org.apache.commons.math3.optim.ConvergenceChecker;
    028    import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
    029    import org.apache.commons.math3.optim.univariate.BracketFinder;
    030    import org.apache.commons.math3.optim.univariate.BrentOptimizer;
    031    import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
    032    import org.apache.commons.math3.optim.univariate.SimpleUnivariateValueChecker;
    033    import org.apache.commons.math3.optim.univariate.SearchInterval;
    034    import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
    035
    036    /**
    037     * Powell algorithm.
    038     * This code is translated and adapted from the Python version of this
    039     * algorithm (as implemented in module {@code optimize.py} v0.5 of
    040     * <em>SciPy</em>).
    041     * <br/>
    042     * The default stopping criterion is based on the differences of the
    043     * function value between two successive iterations. It is however possible
    044     * to define a custom convergence checker that might terminate the algorithm
    045     * earlier.
    046     * <br/>
    047     * The internal line search optimizer is a {@link BrentOptimizer} with a
    048     * convergence checker set to {@link SimpleUnivariateValueChecker}.
    049     *
    050     * @version $Id: PowellOptimizer.java 1413594 2012-11-26 13:16:39Z erans $
    051     * @since 2.2
    052     */
    053    public class PowellOptimizer
    054        extends MultivariateOptimizer {
    055        /**
    056         * Minimum relative tolerance.
    057         */
    058        private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
    059        /**
    060         * Relative threshold.
    061         */
    062        private final double relativeThreshold;
    063        /**
    064         * Absolute threshold.
    065         */
    066        private final double absoluteThreshold;
    067        /**
    068         * Line search.
    069         */
    070        private final LineSearch line;
    071
    072        /**
    073         * This constructor allows to specify a user-defined convergence checker,
    074         * in addition to the parameters that control the default convergence
    075         * checking procedure.
    076         * <br/>
    077         * The internal line search tolerances are set to the square-root of their
    078         * corresponding value in the multivariate optimizer.
    079         *
    080         * @param rel Relative threshold.
    081         * @param abs Absolute threshold.
    082         * @param checker Convergence checker.
    083         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    084         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    085         */
    086        public PowellOptimizer(double rel,
    087                               double abs,
    088                               ConvergenceChecker<PointValuePair> checker) {
    089            this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
    090        }
    091
    092        /**
    093         * This constructor allows to specify a user-defined convergence checker,
    094         * in addition to the parameters that control the default convergence
    095         * checking procedure and the line search tolerances.
    096         *
    097         * @param rel Relative threshold for this optimizer.
    098         * @param abs Absolute threshold for this optimizer.
    099         * @param lineRel Relative threshold for the internal line search optimizer.
    100         * @param lineAbs Absolute threshold for the internal line search optimizer.
    101         * @param checker Convergence checker.
    102         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    103         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    104         */
    105        public PowellOptimizer(double rel,
    106                               double abs,
    107                               double lineRel,
    108                               double lineAbs,
    109                               ConvergenceChecker<PointValuePair> checker) {
    110            super(checker);
    111
    112            if (rel < MIN_RELATIVE_TOLERANCE) {
    113                throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
    114            }
    115            if (abs <= 0) {
    116                throw new NotStrictlyPositiveException(abs);
    117            }
    118            relativeThreshold = rel;
    119            absoluteThreshold = abs;
    120
    121            // Create the line search optimizer.
    122            line = new LineSearch(lineRel,
    123                                  lineAbs);
    124        }
    125
    126        /**
    127         * The parameters control the default convergence checking procedure.
    128         * <br/>
    129         * The internal line search tolerances are set to the square-root of their
    130         * corresponding value in the multivariate optimizer.
    131         *
    132         * @param rel Relative threshold.
    133         * @param abs Absolute threshold.
    134         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    135         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    136         */
    137        public PowellOptimizer(double rel,
    138                               double abs) {
    139            this(rel, abs, null);
    140        }
    141
    142        /**
    143         * Builds an instance with the default convergence checking procedure.
    144         *
    145         * @param rel Relative threshold.
    146         * @param abs Absolute threshold.
    147         * @param lineRel Relative threshold for the internal line search optimizer.
    148         * @param lineAbs Absolute threshold for the internal line search optimizer.
    149         * @throws NotStrictlyPositiveException if {@code abs <= 0}.
    150         * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
    151         */
    152        public PowellOptimizer(double rel,
    153                               double abs,
    154                               double lineRel,
    155                               double lineAbs) {
    156            this(rel, abs, lineRel, lineAbs, null);
    157        }
    158
    159        /** {@inheritDoc} */
    160        @Override
    161        protected PointValuePair doOptimize() {
    162            final GoalType goal = getGoalType();
    163            final double[] guess = getStartPoint();
    164            final int n = guess.length;
    165
    166            final double[][] direc = new double[n][n];
    167            for (int i = 0; i < n; i++) {
    168                direc[i][i] = 1;
    169            }
    170
    171            final ConvergenceChecker<PointValuePair> checker
    172                = getConvergenceChecker();
    173
    174            double[] x = guess;
    175            double fVal = computeObjectiveValue(x);
    176            double[] x1 = x.clone();
    177            int iter = 0;
    178            while (true) {
    179                ++iter;
    180
    181                double fX = fVal;
    182                double fX2 = 0;
    183                double delta = 0;
    184                int bigInd = 0;
    185                double alphaMin = 0;
    186
    187                for (int i = 0; i < n; i++) {
    188                    final double[] d = MathArrays.copyOf(direc[i]);
    189
    190                    fX2 = fVal;
    191
    192                    final UnivariatePointValuePair optimum = line.search(x, d);
    193                    fVal = optimum.getValue();
    194                    alphaMin = optimum.getPoint();
    195                    final double[][] result = newPointAndDirection(x, d, alphaMin);
    196                    x = result[0];
    197
    198                    if ((fX2 - fVal) > delta) {
    199                        delta = fX2 - fVal;
    200                        bigInd = i;
    201                    }
    202                }
    203
    204                // Default convergence check.
    205                boolean stop = 2 * (fX - fVal) <=
    206                    (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
    207                     absoluteThreshold);
    208
    209                final PointValuePair previous = new PointValuePair(x1, fX);
    210                final PointValuePair current = new PointValuePair(x, fVal);
    211                if (!stop) { // User-defined stopping criteria.
    212                    if (checker != null) {
    213                        stop = checker.converged(iter, previous, current);
    214                    }
    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                }
    223
    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                }
    230
    231                x1 = x.clone();
    232                fX2 = computeObjectiveValue(x2);
    233
    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;
    240
    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];
    247
    248                        final int lastInd = n - 1;
    249                        direc[bigInd] = direc[lastInd];
    250                        direc[lastInd] = result[1];
    251                    }
    252                }
    253            }
    254        }
    255
    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            }
    276
    277            final double[][] result = new double[2][];
    278            result[0] = nP;
    279            result[1] = nD;
    280
    281            return result;
    282        }
    283
    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();
    305
    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            }
    320
    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                    };
    342
    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(new MaxEval(Integer.MAX_VALUE),
    349                                new UnivariateObjectiveFunction(f),
    350                                goal,
    351                                new SearchInterval(bracket.getLo(),
    352                                                   bracket.getHi(),
    353                                                   bracket.getMid()));
    354            }
    355        }
    356    }