001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017
018 package org.apache.commons.math3.optimization.direct;
019
020 import org.apache.commons.math3.util.FastMath;
021 import org.apache.commons.math3.util.MathArrays;
022 import org.apache.commons.math3.analysis.UnivariateFunction;
023 import org.apache.commons.math3.analysis.MultivariateFunction;
024 import org.apache.commons.math3.exception.NumberIsTooSmallException;
025 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
026 import org.apache.commons.math3.optimization.GoalType;
027 import org.apache.commons.math3.optimization.PointValuePair;
028 import org.apache.commons.math3.optimization.ConvergenceChecker;
029 import org.apache.commons.math3.optimization.MultivariateOptimizer;
030 import org.apache.commons.math3.optimization.univariate.BracketFinder;
031 import org.apache.commons.math3.optimization.univariate.BrentOptimizer;
032 import org.apache.commons.math3.optimization.univariate.UnivariatePointValuePair;
033 import org.apache.commons.math3.optimization.univariate.SimpleUnivariateValueChecker;
034
035 /**
036 * Powell algorithm.
037 * This code is translated and adapted from the Python version of this
038 * algorithm (as implemented in module {@code optimize.py} v0.5 of
039 * <em>SciPy</em>).
040 * <br/>
041 * The default stopping criterion is based on the differences of the
042 * function value between two successive iterations. It is however possible
043 * to define a custom convergence checker that might terminate the algorithm
044 * earlier.
045 * <br/>
046 * The internal line search optimizer is a {@link BrentOptimizer} with a
047 * convergence checker set to {@link SimpleUnivariateValueChecker}.
048 *
049 * @version $Id: PowellOptimizer.java 1422313 2012-12-15 18:53:41Z psteitz $
050 * @deprecated As of 3.1 (to be removed in 4.0).
051 * @since 2.2
052 */
053 @Deprecated
054 public class PowellOptimizer
055 extends BaseAbstractMultivariateOptimizer<MultivariateFunction>
056 implements MultivariateOptimizer {
057 /**
058 * Minimum relative tolerance.
059 */
060 private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d);
061 /**
062 * Relative threshold.
063 */
064 private final double relativeThreshold;
065 /**
066 * Absolute threshold.
067 */
068 private final double absoluteThreshold;
069 /**
070 * Line search.
071 */
072 private final LineSearch line;
073
074 /**
075 * This constructor allows to specify a user-defined convergence checker,
076 * in addition to the parameters that control the default convergence
077 * checking procedure.
078 * <br/>
079 * The internal line search tolerances are set to the square-root of their
080 * corresponding value in the multivariate optimizer.
081 *
082 * @param rel Relative threshold.
083 * @param abs Absolute threshold.
084 * @param checker Convergence checker.
085 * @throws NotStrictlyPositiveException if {@code abs <= 0}.
086 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}.
087 */
088 public PowellOptimizer(double rel,
089 double abs,
090 ConvergenceChecker<PointValuePair> checker) {
091 this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker);
092 }
093
094 /**
095 * This constructor allows to specify a user-defined convergence checker,
096 * in addition to the parameters that control the default convergence
097 * checking procedure and the line search tolerances.
098 *
099 * @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) { // User-defined stopping criteria.
215 if (checker != null) {
216 stop = checker.converged(iter, previous, current);
217 }
218 }
219 if (stop) {
220 if (goal == GoalType.MINIMIZE) {
221 return (fVal < fX) ? current : previous;
222 } else {
223 return (fVal > fX) ? current : previous;
224 }
225 }
226
227 final double[] d = new double[n];
228 final double[] x2 = new double[n];
229 for (int i = 0; i < n; i++) {
230 d[i] = x[i] - x1[i];
231 x2[i] = 2 * x[i] - x1[i];
232 }
233
234 x1 = x.clone();
235 fX2 = computeObjectiveValue(x2);
236
237 if (fX > fX2) {
238 double t = 2 * (fX + fX2 - 2 * fVal);
239 double temp = fX - fVal - delta;
240 t *= temp * temp;
241 temp = fX - fX2;
242 t -= delta * temp * temp;
243
244 if (t < 0.0) {
245 final UnivariatePointValuePair optimum = line.search(x, d);
246 fVal = optimum.getValue();
247 alphaMin = optimum.getPoint();
248 final double[][] result = newPointAndDirection(x, d, alphaMin);
249 x = result[0];
250
251 final int lastInd = n - 1;
252 direc[bigInd] = direc[lastInd];
253 direc[lastInd] = result[1];
254 }
255 }
256 }
257 }
258
259 /**
260 * Compute a new point (in the original space) and a new direction
261 * vector, resulting from the line search.
262 *
263 * @param p Point used in the line search.
264 * @param d Direction used in the line search.
265 * @param optimum Optimum found by the line search.
266 * @return a 2-element array containing the new point (at index 0) and
267 * the new direction (at index 1).
268 */
269 private double[][] newPointAndDirection(double[] p,
270 double[] d,
271 double optimum) {
272 final int n = p.length;
273 final double[] nP = new double[n];
274 final double[] nD = new double[n];
275 for (int i = 0; i < n; i++) {
276 nD[i] = d[i] * optimum;
277 nP[i] = p[i] + nD[i];
278 }
279
280 final double[][] result = new double[2][];
281 result[0] = nP;
282 result[1] = nD;
283
284 return result;
285 }
286
287 /**
288 * Class for finding the minimum of the objective function along a given
289 * direction.
290 */
291 private class LineSearch extends BrentOptimizer {
292 /**
293 * Value that will pass the precondition check for {@link BrentOptimizer}
294 * but will not pass the convergence check, so that the custom checker
295 * will always decide when to stop the line search.
296 */
297 private static final double REL_TOL_UNUSED = 1e-15;
298 /**
299 * Value that will pass the precondition check for {@link BrentOptimizer}
300 * but will not pass the convergence check, so that the custom checker
301 * will always decide when to stop the line search.
302 */
303 private static final double ABS_TOL_UNUSED = Double.MIN_VALUE;
304 /**
305 * Automatic bracketing.
306 */
307 private final BracketFinder bracket = new BracketFinder();
308
309 /**
310 * The "BrentOptimizer" default stopping criterion uses the tolerances
311 * to check the domain (point) values, not the function values.
312 * We thus create a custom checker to use function values.
313 *
314 * @param rel Relative threshold.
315 * @param abs Absolute threshold.
316 */
317 LineSearch(double rel,
318 double abs) {
319 super(REL_TOL_UNUSED,
320 ABS_TOL_UNUSED,
321 new SimpleUnivariateValueChecker(rel, abs));
322 }
323
324 /**
325 * Find the minimum of the function {@code f(p + alpha * d)}.
326 *
327 * @param p Starting point.
328 * @param d Search direction.
329 * @return the optimum.
330 * @throws org.apache.commons.math3.exception.TooManyEvaluationsException
331 * if the number of evaluations is exceeded.
332 */
333 public UnivariatePointValuePair search(final double[] p, final double[] d) {
334 final int n = p.length;
335 final UnivariateFunction f = new UnivariateFunction() {
336 public double value(double alpha) {
337 final double[] x = new double[n];
338 for (int i = 0; i < n; i++) {
339 x[i] = p[i] + alpha * d[i];
340 }
341 final double obj = PowellOptimizer.this.computeObjectiveValue(x);
342 return obj;
343 }
344 };
345
346 final GoalType goal = PowellOptimizer.this.getGoalType();
347 bracket.search(f, goal, 0, 1);
348 // Passing "MAX_VALUE" as a dummy value because it is the enclosing
349 // class that counts the number of evaluations (and will eventually
350 // generate the exception).
351 return optimize(Integer.MAX_VALUE, f, goal,
352 bracket.getLo(), bracket.getHi(), bracket.getMid());
353 }
354 }
355 }