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