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 18 package org.apache.commons.math3.optim.nonlinear.scalar.gradient; 19 20 import org.apache.commons.math3.analysis.solvers.UnivariateSolver; 21 import org.apache.commons.math3.exception.MathInternalError; 22 import org.apache.commons.math3.exception.TooManyEvaluationsException; 23 import org.apache.commons.math3.exception.MathUnsupportedOperationException; 24 import org.apache.commons.math3.exception.util.LocalizedFormats; 25 import org.apache.commons.math3.optim.OptimizationData; 26 import org.apache.commons.math3.optim.PointValuePair; 27 import org.apache.commons.math3.optim.ConvergenceChecker; 28 import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; 29 import org.apache.commons.math3.optim.nonlinear.scalar.GradientMultivariateOptimizer; 30 import org.apache.commons.math3.optim.nonlinear.scalar.LineSearch; 31 32 33 /** 34 * Non-linear conjugate gradient optimizer. 35 * <br/> 36 * This class supports both the Fletcher-Reeves and the Polak-Ribière 37 * update formulas for the conjugate search directions. 38 * It also supports optional preconditioning. 39 * <br/> 40 * Constraints are not supported: the call to 41 * {@link #optimize(OptimizationData[]) optimize} will throw 42 * {@link MathUnsupportedOperationException} if bounds are passed to it. 43 * 44 * @since 2.0 45 */ 46 public class NonLinearConjugateGradientOptimizer 47 extends GradientMultivariateOptimizer { 48 /** Update formula for the beta parameter. */ 49 private final Formula updateFormula; 50 /** Preconditioner (may be null). */ 51 private final Preconditioner preconditioner; 52 /** Line search algorithm. */ 53 private final LineSearch line; 54 55 /** 56 * Available choices of update formulas for the updating the parameter 57 * that is used to compute the successive conjugate search directions. 58 * For non-linear conjugate gradients, there are 59 * two formulas: 60 * <ul> 61 * <li>Fletcher-Reeves formula</li> 62 * <li>Polak-Ribière formula</li> 63 * </ul> 64 * 65 * On the one hand, the Fletcher-Reeves formula is guaranteed to converge 66 * if the start point is close enough of the optimum whether the 67 * Polak-Ribière formula may not converge in rare cases. On the 68 * other hand, the Polak-Ribière formula is often faster when it 69 * does converge. Polak-Ribière is often used. 70 * 71 * @since 2.0 72 */ 73 public static enum Formula { 74 /** Fletcher-Reeves formula. */ 75 FLETCHER_REEVES, 76 /** Polak-Ribière formula. */ 77 POLAK_RIBIERE 78 } 79 80 /** 81 * The initial step is a factor with respect to the search direction 82 * (which itself is roughly related to the gradient of the function). 83 * <br/> 84 * It is used to find an interval that brackets the optimum in line 85 * search. 86 * 87 * @since 3.1 88 * @deprecated As of v3.3, this class is not used anymore. 89 * This setting is replaced by the {@code initialBracketingRange} 90 * argument to the new constructors. 91 */ 92 @Deprecated 93 public static class BracketingStep implements OptimizationData { 94 /** Initial step. */ 95 private final double initialStep; 96 97 /** 98 * @param step Initial step for the bracket search. 99 */ 100 public BracketingStep(double step) { 101 initialStep = step; 102 } 103 104 /** 105 * Gets the initial step. 106 * 107 * @return the initial step. 108 */ 109 public double getBracketingStep() { 110 return initialStep; 111 } 112 } 113 114 /** 115 * Constructor with default tolerances for the line search (1e-8) and 116 * {@link IdentityPreconditioner preconditioner}. 117 * 118 * @param updateFormula formula to use for updating the β parameter, 119 * must be one of {@link Formula#FLETCHER_REEVES} or 120 * {@link Formula#POLAK_RIBIERE}. 121 * @param checker Convergence checker. 122 */ 123 public NonLinearConjugateGradientOptimizer(final Formula updateFormula, 124 ConvergenceChecker<PointValuePair> checker) { 125 this(updateFormula, 126 checker, 127 1e-8, 128 1e-8, 129 1e-8, 130 new IdentityPreconditioner()); 131 } 132 133 /** 134 * Constructor with default {@link IdentityPreconditioner preconditioner}. 135 * 136 * @param updateFormula formula to use for updating the β parameter, 137 * must be one of {@link Formula#FLETCHER_REEVES} or 138 * {@link Formula#POLAK_RIBIERE}. 139 * @param checker Convergence checker. 140 * @param lineSearchSolver Solver to use during line search. 141 * @deprecated as of 3.3. Please use 142 * {@link #NonLinearConjugateGradientOptimizer(Formula,ConvergenceChecker,double,double,double)} instead. 143 */ 144 @Deprecated 145 public NonLinearConjugateGradientOptimizer(final Formula updateFormula, 146 ConvergenceChecker<PointValuePair> checker, 147 final UnivariateSolver lineSearchSolver) { 148 this(updateFormula, 149 checker, 150 lineSearchSolver, 151 new IdentityPreconditioner()); 152 } 153 154 /** 155 * Constructor with default {@link IdentityPreconditioner preconditioner}. 156 * 157 * @param updateFormula formula to use for updating the β parameter, 158 * must be one of {@link Formula#FLETCHER_REEVES} or 159 * {@link Formula#POLAK_RIBIERE}. 160 * @param checker Convergence checker. 161 * @param relativeTolerance Relative threshold for line search. 162 * @param absoluteTolerance Absolute threshold for line search. 163 * @param initialBracketingRange Extent of the initial interval used to 164 * find an interval that brackets the optimum in order to perform the 165 * line search. 166 * 167 * @see LineSearch#LineSearch(MultivariateOptimizer,double,double,double) 168 * @since 3.3 169 */ 170 public NonLinearConjugateGradientOptimizer(final Formula updateFormula, 171 ConvergenceChecker<PointValuePair> checker, 172 double relativeTolerance, 173 double absoluteTolerance, 174 double initialBracketingRange) { 175 this(updateFormula, 176 checker, 177 relativeTolerance, 178 absoluteTolerance, 179 initialBracketingRange, 180 new IdentityPreconditioner()); 181 } 182 183 /** 184 * @param updateFormula formula to use for updating the β parameter, 185 * must be one of {@link Formula#FLETCHER_REEVES} or 186 * {@link Formula#POLAK_RIBIERE}. 187 * @param checker Convergence checker. 188 * @param lineSearchSolver Solver to use during line search. 189 * @param preconditioner Preconditioner. 190 * @deprecated as of 3.3. Please use 191 * {@link #NonLinearConjugateGradientOptimizer(Formula,ConvergenceChecker,double,double,double,Preconditioner)} instead. 192 */ 193 @Deprecated 194 public NonLinearConjugateGradientOptimizer(final Formula updateFormula, 195 ConvergenceChecker<PointValuePair> checker, 196 final UnivariateSolver lineSearchSolver, 197 final Preconditioner preconditioner) { 198 this(updateFormula, 199 checker, 200 lineSearchSolver.getRelativeAccuracy(), 201 lineSearchSolver.getAbsoluteAccuracy(), 202 lineSearchSolver.getAbsoluteAccuracy(), 203 preconditioner); 204 } 205 206 /** 207 * @param updateFormula formula to use for updating the β parameter, 208 * must be one of {@link Formula#FLETCHER_REEVES} or 209 * {@link Formula#POLAK_RIBIERE}. 210 * @param checker Convergence checker. 211 * @param preconditioner Preconditioner. 212 * @param relativeTolerance Relative threshold for line search. 213 * @param absoluteTolerance Absolute threshold for line search. 214 * @param initialBracketingRange Extent of the initial interval used to 215 * find an interval that brackets the optimum in order to perform the 216 * line search. 217 * 218 * @see LineSearch#LineSearch(MultivariateOptimizer,double,double,double) 219 * @since 3.3 220 */ 221 public NonLinearConjugateGradientOptimizer(final Formula updateFormula, 222 ConvergenceChecker<PointValuePair> checker, 223 double relativeTolerance, 224 double absoluteTolerance, 225 double initialBracketingRange, 226 final Preconditioner preconditioner) { 227 super(checker); 228 229 this.updateFormula = updateFormula; 230 this.preconditioner = preconditioner; 231 line = new LineSearch(this, 232 relativeTolerance, 233 absoluteTolerance, 234 initialBracketingRange); 235 } 236 237 /** 238 * {@inheritDoc} 239 */ 240 @Override 241 public PointValuePair optimize(OptimizationData... optData) 242 throws TooManyEvaluationsException { 243 // Set up base class and perform computation. 244 return super.optimize(optData); 245 } 246 247 /** {@inheritDoc} */ 248 @Override 249 protected PointValuePair doOptimize() { 250 final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker(); 251 final double[] point = getStartPoint(); 252 final GoalType goal = getGoalType(); 253 final int n = point.length; 254 double[] r = computeObjectiveGradient(point); 255 if (goal == GoalType.MINIMIZE) { 256 for (int i = 0; i < n; i++) { 257 r[i] = -r[i]; 258 } 259 } 260 261 // Initial search direction. 262 double[] steepestDescent = preconditioner.precondition(point, r); 263 double[] searchDirection = steepestDescent.clone(); 264 265 double delta = 0; 266 for (int i = 0; i < n; ++i) { 267 delta += r[i] * searchDirection[i]; 268 } 269 270 PointValuePair current = null; 271 while (true) { 272 incrementIterationCount(); 273 274 final double objective = computeObjectiveValue(point); 275 PointValuePair previous = current; 276 current = new PointValuePair(point, objective); 277 if (previous != null && checker.converged(getIterations(), previous, current)) { 278 // We have found an optimum. 279 return current; 280 } 281 282 final double step = line.search(point, searchDirection).getPoint(); 283 284 // Validate new point. 285 for (int i = 0; i < point.length; ++i) { 286 point[i] += step * searchDirection[i]; 287 } 288 289 r = computeObjectiveGradient(point); 290 if (goal == GoalType.MINIMIZE) { 291 for (int i = 0; i < n; ++i) { 292 r[i] = -r[i]; 293 } 294 } 295 296 // Compute beta. 297 final double deltaOld = delta; 298 final double[] newSteepestDescent = preconditioner.precondition(point, r); 299 delta = 0; 300 for (int i = 0; i < n; ++i) { 301 delta += r[i] * newSteepestDescent[i]; 302 } 303 304 final double beta; 305 switch (updateFormula) { 306 case FLETCHER_REEVES: 307 beta = delta / deltaOld; 308 break; 309 case POLAK_RIBIERE: 310 double deltaMid = 0; 311 for (int i = 0; i < r.length; ++i) { 312 deltaMid += r[i] * steepestDescent[i]; 313 } 314 beta = (delta - deltaMid) / deltaOld; 315 break; 316 default: 317 // Should never happen. 318 throw new MathInternalError(); 319 } 320 steepestDescent = newSteepestDescent; 321 322 // Compute conjugate search direction. 323 if (getIterations() % n == 0 || 324 beta < 0) { 325 // Break conjugation: reset search direction. 326 searchDirection = steepestDescent.clone(); 327 } else { 328 // Compute new conjugate search direction. 329 for (int i = 0; i < n; ++i) { 330 searchDirection[i] = steepestDescent[i] + beta * searchDirection[i]; 331 } 332 } 333 } 334 } 335 336 /** 337 * {@inheritDoc} 338 */ 339 @Override 340 protected void parseOptimizationData(OptimizationData... optData) { 341 // Allow base class to register its own data. 342 super.parseOptimizationData(optData); 343 344 checkParameters(); 345 } 346 347 /** Default identity preconditioner. */ 348 public static class IdentityPreconditioner implements Preconditioner { 349 /** {@inheritDoc} */ 350 public double[] precondition(double[] variables, double[] r) { 351 return r.clone(); 352 } 353 } 354 355 // Class is not used anymore (cf. MATH-1092). However, it might 356 // be interesting to create a class similar to "LineSearch", but 357 // that will take advantage that the model's gradient is available. 358 // /** 359 // * Internal class for line search. 360 // * <p> 361 // * The function represented by this class is the dot product of 362 // * the objective function gradient and the search direction. Its 363 // * value is zero when the gradient is orthogonal to the search 364 // * direction, i.e. when the objective function value is a local 365 // * extremum along the search direction. 366 // * </p> 367 // */ 368 // private class LineSearchFunction implements UnivariateFunction { 369 // /** Current point. */ 370 // private final double[] currentPoint; 371 // /** Search direction. */ 372 // private final double[] searchDirection; 373 374 // /** 375 // * @param point Current point. 376 // * @param direction Search direction. 377 // */ 378 // public LineSearchFunction(double[] point, 379 // double[] direction) { 380 // currentPoint = point.clone(); 381 // searchDirection = direction.clone(); 382 // } 383 384 // /** {@inheritDoc} */ 385 // public double value(double x) { 386 // // current point in the search direction 387 // final double[] shiftedPoint = currentPoint.clone(); 388 // for (int i = 0; i < shiftedPoint.length; ++i) { 389 // shiftedPoint[i] += x * searchDirection[i]; 390 // } 391 392 // // gradient of the objective function 393 // final double[] gradient = computeObjectiveGradient(shiftedPoint); 394 395 // // dot product with the search direction 396 // double dotProduct = 0; 397 // for (int i = 0; i < gradient.length; ++i) { 398 // dotProduct += gradient[i] * searchDirection[i]; 399 // } 400 401 // return dotProduct; 402 // } 403 // } 404 405 /** 406 * @throws MathUnsupportedOperationException if bounds were passed to the 407 * {@link #optimize(OptimizationData[]) optimize} method. 408 */ 409 private void checkParameters() { 410 if (getLowerBound() != null || 411 getUpperBound() != null) { 412 throw new MathUnsupportedOperationException(LocalizedFormats.CONSTRAINT); 413 } 414 } 415 }