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.exception.NumberIsTooSmallException;
22  import org.apache.commons.math3.exception.NotStrictlyPositiveException;
23  import org.apache.commons.math3.exception.MathUnsupportedOperationException;
24  import org.apache.commons.math3.exception.util.LocalizedFormats;
25  import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
26  import org.apache.commons.math3.optim.PointValuePair;
27  import org.apache.commons.math3.optim.ConvergenceChecker;
28  import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer;
29  import org.apache.commons.math3.optim.nonlinear.scalar.LineSearch;
30  import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
31  
32  /**
33   * Powell's algorithm.
34   * This code is translated and adapted from the Python version of this
35   * algorithm (as implemented in module {@code optimize.py} v0.5 of
36   * <em>SciPy</em>).
37   * <br/>
38   * The default stopping criterion is based on the differences of the
39   * function value between two successive iterations. It is however possible
40   * to define a custom convergence checker that might terminate the algorithm
41   * earlier.
42   * <br/>
43   * Line search is performed by the {@link LineSearch} class.
44   * <br/>
45   * Constraints are not supported: the call to
46   * {@link #optimize(OptimizationData[]) optimize} will throw
47   * {@link MathUnsupportedOperationException} if bounds are passed to it.
48   * In order to impose simple constraints, the objective function must be
49   * wrapped in an adapter like
50   * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter
51   * MultivariateFunctionMappingAdapter} or
52   * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter
53   * MultivariateFunctionPenaltyAdapter}.
54   *
55   * @version $Id: PowellOptimizer.java 1579346 2014-03-19 18:43:39Z erans $
56   * @since 2.2
57   */
58  public class PowellOptimizer
59      extends MultivariateOptimizer {
60      /**
61       * Minimum relative tolerance.
62       */
63      private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
64      /**
65       * Relative threshold.
66       */
67      private final double relativeThreshold;
68      /**
69       * Absolute threshold.
70       */
71      private final double absoluteThreshold;
72      /**
73       * Line search.
74       */
75      private final LineSearch line;
76  
77      /**
78       * This constructor allows to specify a user-defined convergence checker,
79       * in addition to the parameters that control the default convergence
80       * checking procedure.
81       * <br/>
82       * The internal line search tolerances are set to the square-root of their
83       * corresponding value in the multivariate optimizer.
84       *
85       * @param rel Relative threshold.
86       * @param abs Absolute threshold.
87       * @param checker Convergence checker.
88       * @throws NotStrictlyPositiveException if {@code abs <= 0}.
89       * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
90       */
91      public PowellOptimizer(double rel,
92                             double abs,
93                             ConvergenceChecker<PointValuePair> checker) {
94          this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
95      }
96  
97      /**
98       * This constructor allows to specify a user-defined convergence checker,
99       * in addition to the parameters that control the default convergence
100      * checking procedure and the line search tolerances.
101      *
102      * @param rel Relative threshold for this optimizer.
103      * @param abs Absolute threshold for this optimizer.
104      * @param lineRel Relative threshold for the internal line search optimizer.
105      * @param lineAbs Absolute threshold for the internal line search optimizer.
106      * @param checker Convergence checker.
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                            double lineRel,
113                            double lineAbs,
114                            ConvergenceChecker<PointValuePair> checker) {
115         super(checker);
116 
117         if (rel < MIN_RELATIVE_TOLERANCE) {
118             throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true);
119         }
120         if (abs <= 0) {
121             throw new NotStrictlyPositiveException(abs);
122         }
123         relativeThreshold = rel;
124         absoluteThreshold = abs;
125 
126         // Create the line search optimizer.
127         line = new LineSearch(this,
128                               lineRel,
129                               lineAbs,
130                               1d);
131     }
132 
133     /**
134      * The parameters control the default convergence checking procedure.
135      * <br/>
136      * The internal line search tolerances are set to the square-root of their
137      * corresponding value in the multivariate optimizer.
138      *
139      * @param rel Relative threshold.
140      * @param abs Absolute threshold.
141      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
142      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
143      */
144     public PowellOptimizer(double rel,
145                            double abs) {
146         this(rel, abs, null);
147     }
148 
149     /**
150      * Builds an instance with the default convergence checking procedure.
151      *
152      * @param rel Relative threshold.
153      * @param abs Absolute threshold.
154      * @param lineRel Relative threshold for the internal line search optimizer.
155      * @param lineAbs Absolute threshold for the internal line search optimizer.
156      * @throws NotStrictlyPositiveException if {@code abs <= 0}.
157      * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
158      */
159     public PowellOptimizer(double rel,
160                            double abs,
161                            double lineRel,
162                            double lineAbs) {
163         this(rel, abs, lineRel, lineAbs, null);
164     }
165 
166     /** {@inheritDoc} */
167     @Override
168     protected PointValuePair doOptimize() {
169         checkParameters();
170 
171         final GoalType goal = getGoalType();
172         final double[] guess = getStartPoint();
173         final int n = guess.length;
174 
175         final double[][] direc = new double[n][n];
176         for (int i = 0; i < n; i++) {
177             direc[i][i] = 1;
178         }
179 
180         final ConvergenceChecker<PointValuePair> checker
181             = getConvergenceChecker();
182 
183         double[] x = guess;
184         double fVal = computeObjectiveValue(x);
185         double[] x1 = x.clone();
186         while (true) {
187             incrementIterationCount();
188 
189             double fX = fVal;
190             double fX2 = 0;
191             double delta = 0;
192             int bigInd = 0;
193             double alphaMin = 0;
194 
195             for (int i = 0; i < n; i++) {
196                 final double[] d = MathArrays.copyOf(direc[i]);
197 
198                 fX2 = fVal;
199 
200                 final UnivariatePointValuePair optimum = line.search(x, d);
201                 fVal = optimum.getValue();
202                 alphaMin = optimum.getPoint();
203                 final double[][] result = newPointAndDirection(x, d, alphaMin);
204                 x = result[0];
205 
206                 if ((fX2 - fVal) > delta) {
207                     delta = fX2 - fVal;
208                     bigInd = i;
209                 }
210             }
211 
212             // Default convergence check.
213             boolean stop = 2 * (fX - fVal) <=
214                 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) +
215                  absoluteThreshold);
216 
217             final PointValuePair previous = new PointValuePair(x1, fX);
218             final PointValuePair current = new PointValuePair(x, fVal);
219             if (!stop && checker != null) { // User-defined stopping criteria.
220                 stop = checker.converged(getIterations(), previous, current);
221             }
222             if (stop) {
223                 if (goal == GoalType.MINIMIZE) {
224                     return (fVal < fX) ? current : previous;
225                 } else {
226                     return (fVal > fX) ? current : previous;
227                 }
228             }
229 
230             final double[] d = new double[n];
231             final double[] x2 = new double[n];
232             for (int i = 0; i < n; i++) {
233                 d[i] = x[i] - x1[i];
234                 x2[i] = 2 * x[i] - x1[i];
235             }
236 
237             x1 = x.clone();
238             fX2 = computeObjectiveValue(x2);
239 
240             if (fX > fX2) {
241                 double t = 2 * (fX + fX2 - 2 * fVal);
242                 double temp = fX - fVal - delta;
243                 t *= temp * temp;
244                 temp = fX - fX2;
245                 t -= delta * temp * temp;
246 
247                 if (t < 0.0) {
248                     final UnivariatePointValuePair optimum = line.search(x, d);
249                     fVal = optimum.getValue();
250                     alphaMin = optimum.getPoint();
251                     final double[][] result = newPointAndDirection(x, d, alphaMin);
252                     x = result[0];
253 
254                     final int lastInd = n - 1;
255                     direc[bigInd] = direc[lastInd];
256                     direc[lastInd] = result[1];
257                 }
258             }
259         }
260     }
261 
262     /**
263      * Compute a new point (in the original space) and a new direction
264      * vector, resulting from the line search.
265      *
266      * @param p Point used in the line search.
267      * @param d Direction used in the line search.
268      * @param optimum Optimum found by the line search.
269      * @return a 2-element array containing the new point (at index 0) and
270      * the new direction (at index 1).
271      */
272     private double[][] newPointAndDirection(double[] p,
273                                             double[] d,
274                                             double optimum) {
275         final int n = p.length;
276         final double[] nP = new double[n];
277         final double[] nD = new double[n];
278         for (int i = 0; i < n; i++) {
279             nD[i] = d[i] * optimum;
280             nP[i] = p[i] + nD[i];
281         }
282 
283         final double[][] result = new double[2][];
284         result[0] = nP;
285         result[1] = nD;
286 
287         return result;
288     }
289 
290     /**
291      * @throws MathUnsupportedOperationException if bounds were passed to the
292      * {@link #optimize(OptimizationData[]) optimize} method.
293      */
294     private void checkParameters() {
295         if (getLowerBound() != null ||
296             getUpperBound() != null) {
297             throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT);
298         }
299     }
300 }