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.optimization.general; 19 20 import org.apache.commons.math3.exception.MathIllegalStateException; 21 import org.apache.commons.math3.analysis.UnivariateFunction; 22 import org.apache.commons.math3.analysis.solvers.BrentSolver; 23 import org.apache.commons.math3.analysis.solvers.UnivariateSolver; 24 import org.apache.commons.math3.exception.util.LocalizedFormats; 25 import org.apache.commons.math3.optimization.GoalType; 26 import org.apache.commons.math3.optimization.PointValuePair; 27 import org.apache.commons.math3.optimization.SimpleValueChecker; 28 import org.apache.commons.math3.optimization.ConvergenceChecker; 29 import org.apache.commons.math3.util.FastMath; 30 31 /** 32 * Non-linear conjugate gradient optimizer. 33 * <p> 34 * This class supports both the Fletcher-Reeves and the Polak-Ribière 35 * update formulas for the conjugate search directions. It also supports 36 * optional preconditioning. 37 * </p> 38 * 39 * @deprecated As of 3.1 (to be removed in 4.0). 40 * @since 2.0 41 * 42 */ 43 @Deprecated 44 public class NonLinearConjugateGradientOptimizer 45 extends AbstractScalarDifferentiableOptimizer { 46 /** Update formula for the beta parameter. */ 47 private final ConjugateGradientFormula updateFormula; 48 /** Preconditioner (may be null). */ 49 private final Preconditioner preconditioner; 50 /** solver to use in the line search (may be null). */ 51 private final UnivariateSolver solver; 52 /** Initial step used to bracket the optimum in line search. */ 53 private double initialStep; 54 /** Current point. */ 55 private double[] point; 56 57 /** 58 * Constructor with default {@link SimpleValueChecker checker}, 59 * {@link BrentSolver line search solver} and 60 * {@link IdentityPreconditioner preconditioner}. 61 * 62 * @param updateFormula formula to use for updating the β parameter, 63 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link 64 * ConjugateGradientFormula#POLAK_RIBIERE}. 65 * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()} 66 */ 67 @Deprecated NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula)68 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) { 69 this(updateFormula, 70 new SimpleValueChecker()); 71 } 72 73 /** 74 * Constructor with default {@link BrentSolver line search solver} and 75 * {@link IdentityPreconditioner preconditioner}. 76 * 77 * @param updateFormula formula to use for updating the β parameter, 78 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link 79 * ConjugateGradientFormula#POLAK_RIBIERE}. 80 * @param checker Convergence checker. 81 */ NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula, ConvergenceChecker<PointValuePair> checker)82 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula, 83 ConvergenceChecker<PointValuePair> checker) { 84 this(updateFormula, 85 checker, 86 new BrentSolver(), 87 new IdentityPreconditioner()); 88 } 89 90 91 /** 92 * Constructor with default {@link IdentityPreconditioner preconditioner}. 93 * 94 * @param updateFormula formula to use for updating the β parameter, 95 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link 96 * ConjugateGradientFormula#POLAK_RIBIERE}. 97 * @param checker Convergence checker. 98 * @param lineSearchSolver Solver to use during line search. 99 */ NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula, ConvergenceChecker<PointValuePair> checker, final UnivariateSolver lineSearchSolver)100 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula, 101 ConvergenceChecker<PointValuePair> checker, 102 final UnivariateSolver lineSearchSolver) { 103 this(updateFormula, 104 checker, 105 lineSearchSolver, 106 new IdentityPreconditioner()); 107 } 108 109 /** 110 * @param updateFormula formula to use for updating the β parameter, 111 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link 112 * ConjugateGradientFormula#POLAK_RIBIERE}. 113 * @param checker Convergence checker. 114 * @param lineSearchSolver Solver to use during line search. 115 * @param preconditioner Preconditioner. 116 */ NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula, ConvergenceChecker<PointValuePair> checker, final UnivariateSolver lineSearchSolver, final Preconditioner preconditioner)117 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula, 118 ConvergenceChecker<PointValuePair> checker, 119 final UnivariateSolver lineSearchSolver, 120 final Preconditioner preconditioner) { 121 super(checker); 122 123 this.updateFormula = updateFormula; 124 solver = lineSearchSolver; 125 this.preconditioner = preconditioner; 126 initialStep = 1.0; 127 } 128 129 /** 130 * Set the initial step used to bracket the optimum in line search. 131 * <p> 132 * The initial step is a factor with respect to the search direction, 133 * which itself is roughly related to the gradient of the function 134 * </p> 135 * @param initialStep initial step used to bracket the optimum in line search, 136 * if a non-positive value is used, the initial step is reset to its 137 * default value of 1.0 138 */ setInitialStep(final double initialStep)139 public void setInitialStep(final double initialStep) { 140 if (initialStep <= 0) { 141 this.initialStep = 1.0; 142 } else { 143 this.initialStep = initialStep; 144 } 145 } 146 147 /** {@inheritDoc} */ 148 @Override doOptimize()149 protected PointValuePair doOptimize() { 150 final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker(); 151 point = getStartPoint(); 152 final GoalType goal = getGoalType(); 153 final int n = point.length; 154 double[] r = computeObjectiveGradient(point); 155 if (goal == GoalType.MINIMIZE) { 156 for (int i = 0; i < n; ++i) { 157 r[i] = -r[i]; 158 } 159 } 160 161 // Initial search direction. 162 double[] steepestDescent = preconditioner.precondition(point, r); 163 double[] searchDirection = steepestDescent.clone(); 164 165 double delta = 0; 166 for (int i = 0; i < n; ++i) { 167 delta += r[i] * searchDirection[i]; 168 } 169 170 PointValuePair current = null; 171 int iter = 0; 172 int maxEval = getMaxEvaluations(); 173 while (true) { 174 ++iter; 175 176 final double objective = computeObjectiveValue(point); 177 PointValuePair previous = current; 178 current = new PointValuePair(point, objective); 179 if (previous != null && checker.converged(iter, previous, current)) { 180 // We have found an optimum. 181 return current; 182 } 183 184 // Find the optimal step in the search direction. 185 final UnivariateFunction lsf = new LineSearchFunction(searchDirection); 186 final double uB = findUpperBound(lsf, 0, initialStep); 187 // XXX Last parameters is set to a value close to zero in order to 188 // work around the divergence problem in the "testCircleFitting" 189 // unit test (see MATH-439). 190 final double step = solver.solve(maxEval, lsf, 0, uB, 1e-15); 191 maxEval -= solver.getEvaluations(); // Subtract used up evaluations. 192 193 // Validate new point. 194 for (int i = 0; i < point.length; ++i) { 195 point[i] += step * searchDirection[i]; 196 } 197 198 r = computeObjectiveGradient(point); 199 if (goal == GoalType.MINIMIZE) { 200 for (int i = 0; i < n; ++i) { 201 r[i] = -r[i]; 202 } 203 } 204 205 // Compute beta. 206 final double deltaOld = delta; 207 final double[] newSteepestDescent = preconditioner.precondition(point, r); 208 delta = 0; 209 for (int i = 0; i < n; ++i) { 210 delta += r[i] * newSteepestDescent[i]; 211 } 212 213 final double beta; 214 if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) { 215 beta = delta / deltaOld; 216 } else { 217 double deltaMid = 0; 218 for (int i = 0; i < r.length; ++i) { 219 deltaMid += r[i] * steepestDescent[i]; 220 } 221 beta = (delta - deltaMid) / deltaOld; 222 } 223 steepestDescent = newSteepestDescent; 224 225 // Compute conjugate search direction. 226 if (iter % n == 0 || 227 beta < 0) { 228 // Break conjugation: reset search direction. 229 searchDirection = steepestDescent.clone(); 230 } else { 231 // Compute new conjugate search direction. 232 for (int i = 0; i < n; ++i) { 233 searchDirection[i] = steepestDescent[i] + beta * searchDirection[i]; 234 } 235 } 236 } 237 } 238 239 /** 240 * Find the upper bound b ensuring bracketing of a root between a and b. 241 * 242 * @param f function whose root must be bracketed. 243 * @param a lower bound of the interval. 244 * @param h initial step to try. 245 * @return b such that f(a) and f(b) have opposite signs. 246 * @throws MathIllegalStateException if no bracket can be found. 247 */ findUpperBound(final UnivariateFunction f, final double a, final double h)248 private double findUpperBound(final UnivariateFunction f, 249 final double a, final double h) { 250 final double yA = f.value(a); 251 double yB = yA; 252 for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) { 253 final double b = a + step; 254 yB = f.value(b); 255 if (yA * yB <= 0) { 256 return b; 257 } 258 } 259 throw new MathIllegalStateException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH); 260 } 261 262 /** Default identity preconditioner. */ 263 public static class IdentityPreconditioner implements Preconditioner { 264 265 /** {@inheritDoc} */ precondition(double[] variables, double[] r)266 public double[] precondition(double[] variables, double[] r) { 267 return r.clone(); 268 } 269 } 270 271 /** Internal class for line search. 272 * <p> 273 * The function represented by this class is the dot product of 274 * the objective function gradient and the search direction. Its 275 * value is zero when the gradient is orthogonal to the search 276 * direction, i.e. when the objective function value is a local 277 * extremum along the search direction. 278 * </p> 279 */ 280 private class LineSearchFunction implements UnivariateFunction { 281 /** Search direction. */ 282 private final double[] searchDirection; 283 284 /** Simple constructor. 285 * @param searchDirection search direction 286 */ LineSearchFunction(final double[] searchDirection)287 LineSearchFunction(final double[] searchDirection) { 288 this.searchDirection = searchDirection; 289 } 290 291 /** {@inheritDoc} */ value(double x)292 public double value(double x) { 293 // current point in the search direction 294 final double[] shiftedPoint = point.clone(); 295 for (int i = 0; i < shiftedPoint.length; ++i) { 296 shiftedPoint[i] += x * searchDirection[i]; 297 } 298 299 // gradient of the objective function 300 final double[] gradient = computeObjectiveGradient(shiftedPoint); 301 302 // dot product with the search direction 303 double dotProduct = 0; 304 for (int i = 0; i < gradient.length; ++i) { 305 dotProduct += gradient[i] * searchDirection[i]; 306 } 307 308 return dotProduct; 309 } 310 } 311 } 312