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.math.optimization.direct;
19  
20  import org.apache.commons.math.util.FastMath;
21  import org.apache.commons.math.util.MathArrays;
22  import org.apache.commons.math.analysis.UnivariateRealFunction;
23  import org.apache.commons.math.analysis.MultivariateRealFunction;
24  import org.apache.commons.math.exception.NumberIsTooSmallException;
25  import org.apache.commons.math.exception.NotStrictlyPositiveException;
26  import org.apache.commons.math.optimization.GoalType;
27  import org.apache.commons.math.optimization.RealPointValuePair;
28  import org.apache.commons.math.optimization.ConvergenceChecker;
29  import org.apache.commons.math.optimization.MultivariateRealOptimizer;
30  import org.apache.commons.math.optimization.univariate.BracketFinder;
31  import org.apache.commons.math.optimization.univariate.BrentOptimizer;
32  import org.apache.commons.math.optimization.univariate.UnivariateRealPointValuePair;
33  
34  /**
35   * Powell algorithm.
36   * This code is translated and adapted from the Python version of this
37   * algorithm (as implemented in module {@code optimize.py} v0.5 of
38   * <em>SciPy</em>).
39   * <br/>
40   * The default stopping criterion is based on the differences of the
41   * function value between two successive iterations. It is however possible
42   * to define a custom convergence checker that might terminate the algorithm
43   * earlier.
44   *
45   * @version $Id$
46   * @since 2.2
47   */
48  public class PowellOptimizer
49      extends BaseAbstractScalarOptimizer<MultivariateRealFunction>
50      implements MultivariateRealOptimizer {
51      /**
52       * Minimum relative tolerance.
53       */
54      private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
55      /**
56       * Relative threshold.
57       */
58      private final double relativeThreshold;
59      /**
60       * Absolute threshold.
61       */
62      private final double absoluteThreshold;
63      /**
64       * Line search.
65       */
66      private final LineSearch line;
67  
68      /**
69       * This constructor allows to specify a user-defined convergence checker,
70       * in addition to the parameters that control the default convergence
71       * checking procedure and the line search tolerances.
72       *
73       * @param rel Relative threshold.
74       * @param abs Absolute threshold.
75       * @param checker Convergence checker.
76       * @throws NotStrictlyPositiveException if {@code abs <= 0}.
77       * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
78       */
79      public PowellOptimizer(double rel,
80                             double abs,
81                             ConvergenceChecker<RealPointValuePair> checker) {
82          super(checker);
83  
84          if (rel < MIN_RELATIVE_TOLERANCE) {
85              throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
86          }
87          if (abs <= 0) {
88              throw new NotStrictlyPositiveException(abs);
89          }
90          relativeThreshold = rel;
91          absoluteThreshold = abs;
92  
93          // Line search tolerances can be much lower than the tolerances
94          // required for the optimizer itself.
95          final double minTol = 1e-4;
96          final double lsRel = Math.min(FastMath.sqrt(relativeThreshold), minTol);
97          final double lsAbs = Math.min(FastMath.sqrt(absoluteThreshold), minTol);
98          line = new LineSearch(lsRel, lsAbs);
99      }
100 
101     /**
102      * The parameters control the default convergence checking procedure, and
103      * the line search tolerances.
104      *
105      * @param rel Relative threshold.
106      * @param abs Absolute threshold.
107      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
108      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
109      */
110     public PowellOptimizer(double rel,
111                            double abs) {
112         this(rel, abs, null);
113     }
114 
115     /** {@inheritDoc} */
116     @Override
117     protected RealPointValuePair doOptimize() {
118         final GoalType goal = getGoalType();
119         final double[] guess = getStartPoint();
120         final int n = guess.length;
121 
122         final double[][] direc = new double[n][n];
123         for (int i = 0; i < n; i++) {
124             direc[i][i] = 1;
125         }
126 
127         final ConvergenceChecker<RealPointValuePair> checker
128             = getConvergenceChecker();
129 
130         double[] x = guess;
131         double fVal = computeObjectiveValue(x);
132         double[] x1 = x.clone();
133         int iter = 0;
134         while (true) {
135             ++iter;
136 
137             double fX = fVal;
138             double fX2 = 0;
139             double delta = 0;
140             int bigInd = 0;
141             double alphaMin = 0;
142 
143             for (int i = 0; i < n; i++) {
144                 final double[] d = MathArrays.copyOf(direc[i]);
145 
146                 fX2 = fVal;
147 
148                 final UnivariateRealPointValuePair optimum = line.search(x, d);
149                 fVal = optimum.getValue();
150                 alphaMin = optimum.getPoint();
151                 final double[][] result = newPointAndDirection(x, d, alphaMin);
152                 x = result[0];
153 
154                 if ((fX2 - fVal) > delta) {
155                     delta = fX2 - fVal;
156                     bigInd = i;
157                 }
158             }
159 
160             // Default convergence check.
161             boolean stop = 2 * (fX - fVal) <=
162                 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
163                  absoluteThreshold);
164 
165             final RealPointValuePair previous = new RealPointValuePair(x1, fX);
166             final RealPointValuePair current = new RealPointValuePair(x, fVal);
167             if (!stop) { // User-defined stopping criteria.
168                 if (checker != null) {
169                     stop = checker.converged(iter, previous, current);
170                 }
171             }
172             if (stop) {
173                 if (goal == GoalType.MINIMIZE) {
174                     return (fVal < fX) ? current : previous;
175                 } else {
176                     return (fVal > fX) ? current : previous;
177                 }
178             }
179 
180             final double[] d = new double[n];
181             final double[] x2 = new double[n];
182             for (int i = 0; i < n; i++) {
183                 d[i] = x[i] - x1[i];
184                 x2[i] = 2 * x[i] - x1[i];
185             }
186 
187             x1 = x.clone();
188             fX2 = computeObjectiveValue(x2);
189 
190             if (fX > fX2) {
191                 double t = 2 * (fX + fX2 - 2 * fVal);
192                 double temp = fX - fVal - delta;
193                 t *= temp * temp;
194                 temp = fX - fX2;
195                 t -= delta * temp * temp;
196 
197                 if (t < 0.0) {
198                     final UnivariateRealPointValuePair optimum = line.search(x, d);
199                     fVal = optimum.getValue();
200                     alphaMin = optimum.getPoint();
201                     final double[][] result = newPointAndDirection(x, d, alphaMin);
202                     x = result[0];
203 
204                     final int lastInd = n - 1;
205                     direc[bigInd] = direc[lastInd];
206                     direc[lastInd] = result[1];
207                 }
208             }
209         }
210     }
211 
212     /**
213      * Compute a new point (in the original space) and a new direction
214      * vector, resulting from the line search.
215      * The parameters {@code p} and {@code d} will be changed in-place.
216      *
217      * @param p Point used in the line search.
218      * @param d Direction used in the line search.
219      * @param optimum Optimum found by the line search.
220      * @return a 2-element array containing the new point (at index 0) and
221      * the new direction (at index 1).
222      */
223     private double[][] newPointAndDirection(double[] p,
224                                             double[] d,
225                                             double optimum) {
226         final int n = p.length;
227         final double[][] result = new double[2][n];
228         final double[] nP = result[0];
229         final double[] nD = result[1];
230         for (int i = 0; i < n; i++) {
231             nD[i] = d[i] * optimum;
232             nP[i] = p[i] + nD[i];
233         }
234         return result;
235     }
236 
237     /**
238      * Class for finding the minimum of the objective function along a given
239      * direction.
240      */
241     private class LineSearch extends BrentOptimizer {
242         /**
243          * Automatic bracketing.
244          */
245         private final BracketFinder bracket = new BracketFinder();
246 
247         /**
248          * @param rel Relative threshold.
249          * @param abs Absolute threshold.
250          */
251         LineSearch(double rel,
252                    double abs) {
253             super(rel, abs);
254         }
255 
256         /**
257          * Find the minimum of the function {@code f(p + alpha * d)}.
258          *
259          * @param p Starting point.
260          * @param d Search direction.
261          * @return the optimum.
262          * @throws org.apache.commons.math.exception.TooManyEvaluationsException
263          * if the number of evaluations is exceeded.
264          */
265         public UnivariateRealPointValuePair search(final double[] p, final double[] d) {
266             final int n = p.length;
267             final UnivariateRealFunction f = new UnivariateRealFunction() {
268                     public double value(double alpha) {
269                         final double[] x = new double[n];
270                         for (int i = 0; i < n; i++) {
271                             x[i] = p[i] + alpha * d[i];
272                         }
273                         final double obj = PowellOptimizer.this.computeObjectiveValue(x);
274                         return obj;
275                     }
276                 };
277 
278             final GoalType goal = PowellOptimizer.this.getGoalType();
279             bracket.search(f, goal, 0, 1);
280             // Passing "MAX_VALUE" as a dummy value because it is the enclosing
281             // class that counts the number of evaluations (and will eventually
282             // generate the exception).
283             return optimize(Integer.MAX_VALUE, f, goal,
284                             bracket.getLo(), bracket.getHi(), bracket.getMid());
285         }
286     }
287 }