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 */ 017package org.apache.commons.math3.optim.nonlinear.scalar.noderiv; 018 019import org.apache.commons.math3.util.FastMath; 020import org.apache.commons.math3.util.MathArrays; 021import org.apache.commons.math3.exception.NumberIsTooSmallException; 022import org.apache.commons.math3.exception.NotStrictlyPositiveException; 023import org.apache.commons.math3.exception.MathUnsupportedOperationException; 024import org.apache.commons.math3.exception.util.LocalizedFormats; 025import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; 026import org.apache.commons.math3.optim.PointValuePair; 027import org.apache.commons.math3.optim.ConvergenceChecker; 028import org.apache.commons.math3.optim.nonlinear.scalar.MultivariateOptimizer; 029import org.apache.commons.math3.optim.nonlinear.scalar.LineSearch; 030import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair; 031 032/** 033 * Powell's algorithm. 034 * This code is translated and adapted from the Python version of this 035 * algorithm (as implemented in module {@code optimize.py} v0.5 of 036 * <em>SciPy</em>). 037 * <br/> 038 * The default stopping criterion is based on the differences of the 039 * function value between two successive iterations. It is however possible 040 * to define a custom convergence checker that might terminate the algorithm 041 * earlier. 042 * <br/> 043 * Line search is performed by the {@link LineSearch} class. 044 * <br/> 045 * Constraints are not supported: the call to 046 * {@link #optimize(OptimizationData[]) optimize} will throw 047 * {@link MathUnsupportedOperationException} if bounds are passed to it. 048 * In order to impose simple constraints, the objective function must be 049 * wrapped in an adapter like 050 * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionMappingAdapter 051 * MultivariateFunctionMappingAdapter} or 052 * {@link org.apache.commons.math3.optim.nonlinear.scalar.MultivariateFunctionPenaltyAdapter 053 * MultivariateFunctionPenaltyAdapter}. 054 * 055 * @since 2.2 056 */ 057public class PowellOptimizer 058 extends MultivariateOptimizer { 059 /** 060 * Minimum relative tolerance. 061 */ 062 private static final double MIN_RELATIVE_TOLERANCE = 2 * FastMath.ulp(1d); 063 /** 064 * Relative threshold. 065 */ 066 private final double relativeThreshold; 067 /** 068 * Absolute threshold. 069 */ 070 private final double absoluteThreshold; 071 /** 072 * Line search. 073 */ 074 private final LineSearch line; 075 076 /** 077 * This constructor allows to specify a user-defined convergence checker, 078 * in addition to the parameters that control the default convergence 079 * checking procedure. 080 * <br/> 081 * The internal line search tolerances are set to the square-root of their 082 * corresponding value in the multivariate optimizer. 083 * 084 * @param rel Relative threshold. 085 * @param abs Absolute threshold. 086 * @param checker Convergence checker. 087 * @throws NotStrictlyPositiveException if {@code abs <= 0}. 088 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. 089 */ 090 public PowellOptimizer(double rel, 091 double abs, 092 ConvergenceChecker<PointValuePair> checker) { 093 this(rel, abs, FastMath.sqrt(rel), FastMath.sqrt(abs), checker); 094 } 095 096 /** 097 * This constructor allows to specify a user-defined convergence checker, 098 * in addition to the parameters that control the default convergence 099 * checking procedure and the line search tolerances. 100 * 101 * @param rel Relative threshold for this optimizer. 102 * @param abs Absolute threshold for this optimizer. 103 * @param lineRel Relative threshold for the internal line search optimizer. 104 * @param lineAbs Absolute threshold for the internal line search optimizer. 105 * @param checker Convergence checker. 106 * @throws NotStrictlyPositiveException if {@code abs <= 0}. 107 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. 108 */ 109 public PowellOptimizer(double rel, 110 double abs, 111 double lineRel, 112 double lineAbs, 113 ConvergenceChecker<PointValuePair> checker) { 114 super(checker); 115 116 if (rel < MIN_RELATIVE_TOLERANCE) { 117 throw new NumberIsTooSmallException(rel, MIN_RELATIVE_TOLERANCE, true); 118 } 119 if (abs <= 0) { 120 throw new NotStrictlyPositiveException(abs); 121 } 122 relativeThreshold = rel; 123 absoluteThreshold = abs; 124 125 // Create the line search optimizer. 126 line = new LineSearch(this, 127 lineRel, 128 lineAbs, 129 1d); 130 } 131 132 /** 133 * The parameters control the default convergence checking procedure. 134 * <br/> 135 * The internal line search tolerances are set to the square-root of their 136 * corresponding value in the multivariate optimizer. 137 * 138 * @param rel Relative threshold. 139 * @param abs Absolute threshold. 140 * @throws NotStrictlyPositiveException if {@code abs <= 0}. 141 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. 142 */ 143 public PowellOptimizer(double rel, 144 double abs) { 145 this(rel, abs, null); 146 } 147 148 /** 149 * Builds an instance with the default convergence checking procedure. 150 * 151 * @param rel Relative threshold. 152 * @param abs Absolute threshold. 153 * @param lineRel Relative threshold for the internal line search optimizer. 154 * @param lineAbs Absolute threshold for the internal line search optimizer. 155 * @throws NotStrictlyPositiveException if {@code abs <= 0}. 156 * @throws NumberIsTooSmallException if {@code rel < 2 * Math.ulp(1d)}. 157 */ 158 public PowellOptimizer(double rel, 159 double abs, 160 double lineRel, 161 double lineAbs) { 162 this(rel, abs, lineRel, lineAbs, null); 163 } 164 165 /** {@inheritDoc} */ 166 @Override 167 protected PointValuePair doOptimize() { 168 checkParameters(); 169 170 final GoalType goal = getGoalType(); 171 final double[] guess = getStartPoint(); 172 final int n = guess.length; 173 174 final double[][] direc = new double[n][n]; 175 for (int i = 0; i < n; i++) { 176 direc[i][i] = 1; 177 } 178 179 final ConvergenceChecker<PointValuePair> checker 180 = getConvergenceChecker(); 181 182 double[] x = guess; 183 double fVal = computeObjectiveValue(x); 184 double[] x1 = x.clone(); 185 while (true) { 186 incrementIterationCount(); 187 188 double fX = fVal; 189 double fX2 = 0; 190 double delta = 0; 191 int bigInd = 0; 192 double alphaMin = 0; 193 194 for (int i = 0; i < n; i++) { 195 final double[] d = MathArrays.copyOf(direc[i]); 196 197 fX2 = fVal; 198 199 final UnivariatePointValuePair optimum = line.search(x, d); 200 fVal = optimum.getValue(); 201 alphaMin = optimum.getPoint(); 202 final double[][] result = newPointAndDirection(x, d, alphaMin); 203 x = result[0]; 204 205 if ((fX2 - fVal) > delta) { 206 delta = fX2 - fVal; 207 bigInd = i; 208 } 209 } 210 211 // Default convergence check. 212 boolean stop = 2 * (fX - fVal) <= 213 (relativeThreshold * (FastMath.abs(fX) + FastMath.abs(fVal)) + 214 absoluteThreshold); 215 216 final PointValuePair previous = new PointValuePair(x1, fX); 217 final PointValuePair current = new PointValuePair(x, fVal); 218 if (!stop && checker != null) { // User-defined stopping criteria. 219 stop = checker.converged(getIterations(), previous, current); 220 } 221 if (stop) { 222 if (goal == GoalType.MINIMIZE) { 223 return (fVal < fX) ? current : previous; 224 } else { 225 return (fVal > fX) ? current : previous; 226 } 227 } 228 229 final double[] d = new double[n]; 230 final double[] x2 = new double[n]; 231 for (int i = 0; i < n; i++) { 232 d[i] = x[i] - x1[i]; 233 x2[i] = 2 * x[i] - x1[i]; 234 } 235 236 x1 = x.clone(); 237 fX2 = computeObjectiveValue(x2); 238 239 if (fX > fX2) { 240 double t = 2 * (fX + fX2 - 2 * fVal); 241 double temp = fX - fVal - delta; 242 t *= temp * temp; 243 temp = fX - fX2; 244 t -= delta * temp * temp; 245 246 if (t < 0.0) { 247 final UnivariatePointValuePair optimum = line.search(x, d); 248 fVal = optimum.getValue(); 249 alphaMin = optimum.getPoint(); 250 final double[][] result = newPointAndDirection(x, d, alphaMin); 251 x = result[0]; 252 253 final int lastInd = n - 1; 254 direc[bigInd] = direc[lastInd]; 255 direc[lastInd] = result[1]; 256 } 257 } 258 } 259 } 260 261 /** 262 * Compute a new point (in the original space) and a new direction 263 * vector, resulting from the line search. 264 * 265 * @param p Point used in the line search. 266 * @param d Direction used in the line search. 267 * @param optimum Optimum found by the line search. 268 * @return a 2-element array containing the new point (at index 0) and 269 * the new direction (at index 1). 270 */ 271 private double[][] newPointAndDirection(double[] p, 272 double[] d, 273 double optimum) { 274 final int n = p.length; 275 final double[] nP = new double[n]; 276 final double[] nD = new double[n]; 277 for (int i = 0; i < n; i++) { 278 nD[i] = d[i] * optimum; 279 nP[i] = p[i] + nD[i]; 280 } 281 282 final double[][] result = new double[2][]; 283 result[0] = nP; 284 result[1] = nD; 285 286 return result; 287 } 288 289 /** 290 * @throws MathUnsupportedOperationException if bounds were passed to the 291 * {@link #optimize(OptimizationData[]) optimize} method. 292 */ 293 private void checkParameters() { 294 if (getLowerBound() != null || 295 getUpperBound() != null) { 296 throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT); 297 } 298 } 299}