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    }