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 }