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 }