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 018package org.apache.commons.math3.optim.nonlinear.scalar.gradient; 019 020import org.apache.commons.math3.analysis.solvers.UnivariateSolver; 021import org.apache.commons.math3.exception.MathInternalError; 022import org.apache.commons.math3.exception.TooManyEvaluationsException; 023import org.apache.commons.math3.exception.MathUnsupportedOperationException; 024import org.apache.commons.math3.exception.util.LocalizedFormats; 025import org.apache.commons.math3.optim.OptimizationData; 026import org.apache.commons.math3.optim.PointValuePair; 027import org.apache.commons.math3.optim.ConvergenceChecker; 028import org.apache.commons.math3.optim.nonlinear.scalar.GoalType; 029import org.apache.commons.math3.optim.nonlinear.scalar.GradientMultivariateOptimizer; 030import org.apache.commons.math3.optim.nonlinear.scalar.LineSearch; 031 032 033/** 034 * Non-linear conjugate gradient optimizer. 035 * <br/> 036 * This class supports both the Fletcher-Reeves and the Polak-Ribière 037 * update formulas for the conjugate search directions. 038 * It also supports optional preconditioning. 039 * <br/> 040 * Constraints are not supported: the call to 041 * {@link #optimize(OptimizationData[]) optimize} will throw 042 * {@link MathUnsupportedOperationException} if bounds are passed to it. 043 * 044 * @since 2.0 045 */ 046public class NonLinearConjugateGradientOptimizer 047 extends GradientMultivariateOptimizer { 048 /** Update formula for the beta parameter. */ 049 private final Formula updateFormula; 050 /** Preconditioner (may be null). */ 051 private final Preconditioner preconditioner; 052 /** Line search algorithm. */ 053 private final LineSearch line; 054 055 /** 056 * Available choices of update formulas for the updating the parameter 057 * that is used to compute the successive conjugate search directions. 058 * For non-linear conjugate gradients, there are 059 * two formulas: 060 * <ul> 061 * <li>Fletcher-Reeves formula</li> 062 * <li>Polak-Ribière formula</li> 063 * </ul> 064 * 065 * On the one hand, the Fletcher-Reeves formula is guaranteed to converge 066 * if the start point is close enough of the optimum whether the 067 * Polak-Ribière formula may not converge in rare cases. On the 068 * other hand, the Polak-Ribière formula is often faster when it 069 * does converge. Polak-Ribière is often used. 070 * 071 * @since 2.0 072 */ 073 public enum Formula { 074 /** Fletcher-Reeves formula. */ 075 FLETCHER_REEVES, 076 /** Polak-Ribière formula. */ 077 POLAK_RIBIERE 078 } 079 080 /** 081 * The initial step is a factor with respect to the search direction 082 * (which itself is roughly related to the gradient of the function). 083 * <br/> 084 * It is used to find an interval that brackets the optimum in line 085 * search. 086 * 087 * @since 3.1 088 * @deprecated As of v3.3, this class is not used anymore. 089 * This setting is replaced by the {@code initialBracketingRange} 090 * argument to the new constructors. 091 */ 092 @Deprecated 093 public static class BracketingStep implements OptimizationData { 094 /** Initial step. */ 095 private final double initialStep; 096 097 /** 098 * @param step Initial step for the bracket search. 099 */ 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}