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