View Javadoc

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;
18  
19  import org.apache.commons.math3.util.FastMath;
20  import org.apache.commons.math3.util.MathArrays;
21  import org.apache.commons.math3.analysis.UnivariateFunction;
22  import org.apache.commons.math3.exception.NumberIsTooSmallException;
23  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
24  import org.apache.commons.math3.exception.MathUnsupportedOperationException;
25  import org.apache.commons.math3.exception.util.LocalizedFormats;
26  import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
27  import org.apache.commons.math3.optim.MaxEval;
28  import org.apache.commons.math3.optim.PointValuePair;
29  import org.apache.commons.math3.optim.ConvergenceChecker;
30  import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
31  import org.apache.commons.math3.optim.univariate.BracketFinder;
32  import org.apache.commons.math3.optim.univariate.BrentOptimizer;
33  import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
34  import org.apache.commons.math3.optim.univariate.SimpleUnivariateValueChecker;
35  import org.apache.commons.math3.optim.univariate.SearchInterval;
36  import org.apache.commons.math3.optim.univariate.UnivariateObjectiveFunction;
37  
38  /**
39   * Powell's algorithm.
40   * This code is translated and adapted from the Python version of this
41   * algorithm (as implemented in module {@code optimize.py} v0.5 of
42   * <em>SciPy</em>).
43   * <br/>
44   * The default stopping criterion is based on the differences of the
45   * function value between two successive iterations. It is however possible
46   * to define a custom convergence checker that might terminate the algorithm
47   * earlier.
48   * <br/>
49   * The internal line search optimizer is a {@link BrentOptimizer} with a
50   * convergence checker set to {@link SimpleUnivariateValueChecker}.
51   * <br/>
52   * Constraints are not supported: the call to
53   * {@link #optimize(OptimizationData[]) optimize} will throw
54   * {@link MathUnsupportedOperationException} if bounds are passed to it.
55   * In order to impose simple constraints, the objective function must be
56   * wrapped in an adapter like
57   * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
58   * MultivariateFunctionMappingAdapter} or
59   * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
60   * MultivariateFunctionPenaltyAdapter}.
61   *
62   * @version $Id: PowellOptimizer.java 1462503 2013-03-29 15:48:27Z luc $
63   * @since 2.2
64   */
65  public class PowellOptimizer
66      extends MultivariateOptimizer {
67      /**
68       * Minimum relative tolerance.
69       */
70      private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
71      /**
72       * Relative threshold.
73       */
74      private final double relativeThreshold;
75      /**
76       * Absolute threshold.
77       */
78      private final double absoluteThreshold;
79      /**
80       * Line search.
81       */
82      private final LineSearch line;
83  
84      /**
85       * This constructor allows to specify a user-defined convergence checker,
86       * in addition to the parameters that control the default convergence
87       * checking procedure.
88       * <br/>
89       * The internal line search tolerances are set to the square-root of their
90       * corresponding value in the multivariate optimizer.
91       *
92       * @param rel Relative threshold.
93       * @param abs Absolute threshold.
94       * @param checker Convergence checker.
95       * @throws NotStrictlyPositiveException if {@code abs <= 0}.
96       * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
97       */
98      public PowellOptimizer(double rel,
99                             double abs,
100                            ConvergenceChecker<PointValuePair> checker) {
101         this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
102     }
103 
104     /**
105      * This constructor allows to specify a user-defined convergence checker,
106      * in addition to the parameters that control the default convergence
107      * checking procedure and the line search tolerances.
108      *
109      * @param rel Relative threshold for this optimizer.
110      * @param abs Absolute threshold for this optimizer.
111      * @param lineRel Relative threshold for the internal line search optimizer.
112      * @param lineAbs Absolute threshold for the internal line search optimizer.
113      * @param checker Convergence checker.
114      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
115      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
116      */
117     public PowellOptimizer(double rel,
118                            double abs,
119                            double lineRel,
120                            double lineAbs,
121                            ConvergenceChecker<PointValuePair> checker) {
122         super(checker);
123 
124         if (rel < MIN_RELATIVE_TOLERANCE) {
125             throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
126         }
127         if (abs <= 0) {
128             throw new NotStrictlyPositiveException(abs);
129         }
130         relativeThreshold = rel;
131         absoluteThreshold = abs;
132 
133         // Create the line search optimizer.
134         line = new LineSearch(lineRel,
135                               lineAbs);
136     }
137 
138     /**
139      * The parameters control the default convergence checking procedure.
140      * <br/>
141      * The internal line search tolerances are set to the square-root of their
142      * corresponding value in the multivariate optimizer.
143      *
144      * @param rel Relative threshold.
145      * @param abs Absolute threshold.
146      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
147      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
148      */
149     public PowellOptimizer(double rel,
150                            double abs) {
151         this(rel, abs, null);
152     }
153 
154     /**
155      * Builds an instance with the default convergence checking procedure.
156      *
157      * @param rel Relative threshold.
158      * @param abs Absolute threshold.
159      * @param lineRel Relative threshold for the internal line search optimizer.
160      * @param lineAbs Absolute threshold for the internal line search optimizer.
161      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
162      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
163      */
164     public PowellOptimizer(double rel,
165                            double abs,
166                            double lineRel,
167                            double lineAbs) {
168         this(rel, abs, lineRel, lineAbs, null);
169     }
170 
171     /** {@inheritDoc} */
172     @Override
173     protected PointValuePair doOptimize() {
174         checkParameters();
175 
176         final GoalType goal = getGoalType();
177         final double[] guess = getStartPoint();
178         final int n = guess.length;
179 
180         final double[][] direc = new double[n][n];
181         for (int i = 0; i < n; i++) {
182             direc[i][i] = 1;
183         }
184 
185         final ConvergenceChecker<PointValuePair> checker
186             = getConvergenceChecker();
187 
188         double[] x = guess;
189         double fVal = computeObjectiveValue(x);
190         double[] x1 = x.clone();
191         while (true) {
192             incrementIterationCount();
193 
194             double fX = fVal;
195             double fX2 = 0;
196             double delta = 0;
197             int bigInd = 0;
198             double alphaMin = 0;
199 
200             for (int i = 0; i < n; i++) {
201                 final double[] d = MathArrays.copyOf(direc[i]);
202 
203                 fX2 = fVal;
204 
205                 final UnivariatePointValuePair optimum = line.search(x, d);
206                 fVal = optimum.getValue();
207                 alphaMin = optimum.getPoint();
208                 final double[][] result = newPointAndDirection(x, d, alphaMin);
209                 x = result[0];
210 
211                 if ((fX2 - fVal) > delta) {
212                     delta = fX2 - fVal;
213                     bigInd = i;
214                 }
215             }
216 
217             // Default convergence check.
218             boolean stop = 2 * (fX - fVal) <=
219                 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
220                  absoluteThreshold);
221 
222             final PointValuePair previous = new PointValuePair(x1, fX);
223             final PointValuePair current = new PointValuePair(x, fVal);
224             if (!stop && checker != null) { // User-defined stopping criteria.
225                 stop = checker.converged(getIterations(), previous, current);
226             }
227             if (stop) {
228                 if (goal == GoalType.MINIMIZE) {
229                     return (fVal < fX) ? current : previous;
230                 } else {
231                     return (fVal > fX) ? current : previous;
232                 }
233             }
234 
235             final double[] d = new double[n];
236             final double[] x2 = new double[n];
237             for (int i = 0; i < n; i++) {
238                 d[i] = x[i] - x1[i];
239                 x2[i] = 2 * x[i] - x1[i];
240             }
241 
242             x1 = x.clone();
243             fX2 = computeObjectiveValue(x2);
244 
245             if (fX > fX2) {
246                 double t = 2 * (fX + fX2 - 2 * fVal);
247                 double temp = fX - fVal - delta;
248                 t *= temp * temp;
249                 temp = fX - fX2;
250                 t -= delta * temp * temp;
251 
252                 if (t < 0.0) {
253                     final UnivariatePointValuePair optimum = line.search(x, d);
254                     fVal = optimum.getValue();
255                     alphaMin = optimum.getPoint();
256                     final double[][] result = newPointAndDirection(x, d, alphaMin);
257                     x = result[0];
258 
259                     final int lastInd = n - 1;
260                     direc[bigInd] = direc[lastInd];
261                     direc[lastInd] = result[1];
262                 }
263             }
264         }
265     }
266 
267     /**
268      * Compute a new point (in the original space) and a new direction
269      * vector, resulting from the line search.
270      *
271      * @param p Point used in the line search.
272      * @param d Direction used in the line search.
273      * @param optimum Optimum found by the line search.
274      * @return a 2-element array containing the new point (at index 0) and
275      * the new direction (at index 1).
276      */
277     private double[][] newPointAndDirection(double[] p,
278                                             double[] d,
279                                             double optimum) {
280         final int n = p.length;
281         final double[] nP = new double[n];
282         final double[] nD = new double[n];
283         for (int i = 0; i < n; i++) {
284             nD[i] = d[i] * optimum;
285             nP[i] = p[i] + nD[i];
286         }
287 
288         final double[][] result = new double[2][];
289         result[0] = nP;
290         result[1] = nD;
291 
292         return result;
293     }
294 
295     /**
296      * Class for finding the minimum of the objective function along a given
297      * direction.
298      */
299     private class LineSearch extends BrentOptimizer {
300         /**
301          * Value that will pass the precondition check for {@link BrentOptimizer}
302          * but will not pass the convergence check, so that the custom checker
303          * will always decide when to stop the line search.
304          */
305         private static final double REL_TOL_UNUSED = 1e-15;
306         /**
307          * Value that will pass the precondition check for {@link BrentOptimizer}
308          * but will not pass the convergence check, so that the custom checker
309          * will always decide when to stop the line search.
310          */
311         private static final double ABS_TOL_UNUSED = Double.MIN_VALUE;
312         /**
313          * Automatic bracketing.
314          */
315         private final BracketFinder bracket = new BracketFinder();
316 
317         /**
318          * The "BrentOptimizer" default stopping criterion uses the tolerances
319          * to check the domain (point) values, not the function values.
320          * We thus create a custom checker to use function values.
321          *
322          * @param rel Relative threshold.
323          * @param abs Absolute threshold.
324          */
325         LineSearch(double rel,
326                    double abs) {
327             super(REL_TOL_UNUSED,
328                   ABS_TOL_UNUSED,
329                   new SimpleUnivariateValueChecker(rel, abs));
330         }
331 
332         /**
333          * Find the minimum of the function {@code f(p + alpha * d)}.
334          *
335          * @param p Starting point.
336          * @param d Search direction.
337          * @return the optimum.
338          * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
339          * if the number of evaluations is exceeded.
340          */
341         public UnivariatePointValuePair search(final double[] p, final double[] d) {
342             final int n = p.length;
343             final UnivariateFunction f = new UnivariateFunction() {
344                     public double value(double alpha) {
345                         final double[] x = new double[n];
346                         for (int i = 0; i < n; i++) {
347                             x[i] = p[i] + alpha * d[i];
348                         }
349                         final double obj = PowellOptimizer.this.computeObjectiveValue(x);
350                         return obj;
351                     }
352                 };
353 
354             final GoalType goal = PowellOptimizer.this.getGoalType();
355             bracket.search(f, goal, 0, 1);
356             // Passing "MAX_VALUE" as a dummy value because it is the enclosing
357             // class that counts the number of evaluations (and will eventually
358             // generate the exception).
359             return optimize(new MaxEval(Integer.MAX_VALUE),
360                             new UnivariateObjectiveFunction(f),
361                             goal,
362                             new SearchInterval(bracket.getLo(),
363                                                bracket.getHi(),
364                                                bracket.getMid()));
365         }
366     }
367 
368     /**
369      * @throws MathUnsupportedOperationException if bounds were passed to the
370      * {@link #optimize(OptimizationData[]) optimize} method.
371      */
372     private void checkParameters() {
373         if (getLowerBound() != null ||
374             getUpperBound() != null) {
375             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
376         }
377     }
378 }