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
018 package org.apache.commons.math3.optimization.general;
019
020 import org.apache.commons.math3.exception.MathIllegalStateException;
021 import org.apache.commons.math3.analysis.UnivariateFunction;
022 import org.apache.commons.math3.analysis.solvers.BrentSolver;
023 import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
024 import org.apache.commons.math3.exception.util.LocalizedFormats;
025 import org.apache.commons.math3.optimization.GoalType;
026 import org.apache.commons.math3.optimization.PointValuePair;
027 import org.apache.commons.math3.optimization.SimpleValueChecker;
028 import org.apache.commons.math3.optimization.ConvergenceChecker;
029 import org.apache.commons.math3.util.FastMath;
030
031 /**
032 * Non-linear conjugate gradient optimizer.
033 * <p>
034 * This class supports both the Fletcher-Reeves and the Polak-Ribière
035 * update formulas for the conjugate search directions. It also supports
036 * optional preconditioning.
037 * </p>
038 *
039 * @version $Id: NonLinearConjugateGradientOptimizer.java 1462503 2013-03-29 15:48:27Z luc $
040 * @deprecated As of 3.1 (to be removed in 4.0).
041 * @since 2.0
042 *
043 */
044 @Deprecated
045 public class NonLinearConjugateGradientOptimizer
046 extends AbstractScalarDifferentiableOptimizer {
047 /** Update formula for the beta parameter. */
048 private final ConjugateGradientFormula updateFormula;
049 /** Preconditioner (may be null). */
050 private final Preconditioner preconditioner;
051 /** solver to use in the line search (may be null). */
052 private final UnivariateSolver solver;
053 /** Initial step used to bracket the optimum in line search. */
054 private double initialStep;
055 /** Current point. */
056 private double[] point;
057
058 /**
059 * Constructor with default {@link SimpleValueChecker checker},
060 * {@link BrentSolver line search solver} and
061 * {@link IdentityPreconditioner preconditioner}.
062 *
063 * @param updateFormula formula to use for updating the β parameter,
064 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
065 * ConjugateGradientFormula#POLAK_RIBIERE}.
066 * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()}
067 */
068 @Deprecated
069 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) {
070 this(updateFormula,
071 new SimpleValueChecker());
072 }
073
074 /**
075 * Constructor with default {@link BrentSolver line search solver} and
076 * {@link IdentityPreconditioner preconditioner}.
077 *
078 * @param updateFormula formula to use for updating the β parameter,
079 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
080 * ConjugateGradientFormula#POLAK_RIBIERE}.
081 * @param checker Convergence checker.
082 */
083 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
084 ConvergenceChecker<PointValuePair> checker) {
085 this(updateFormula,
086 checker,
087 new BrentSolver(),
088 new IdentityPreconditioner());
089 }
090
091
092 /**
093 * Constructor with default {@link IdentityPreconditioner preconditioner}.
094 *
095 * @param updateFormula formula to use for updating the β parameter,
096 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
097 * ConjugateGradientFormula#POLAK_RIBIERE}.
098 * @param checker Convergence checker.
099 * @param lineSearchSolver Solver to use during line search.
100 */
101 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
102 ConvergenceChecker<PointValuePair> checker,
103 final UnivariateSolver lineSearchSolver) {
104 this(updateFormula,
105 checker,
106 lineSearchSolver,
107 new IdentityPreconditioner());
108 }
109
110 /**
111 * @param updateFormula formula to use for updating the β parameter,
112 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
113 * ConjugateGradientFormula#POLAK_RIBIERE}.
114 * @param checker Convergence checker.
115 * @param lineSearchSolver Solver to use during line search.
116 * @param preconditioner Preconditioner.
117 */
118 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
119 ConvergenceChecker<PointValuePair> checker,
120 final UnivariateSolver lineSearchSolver,
121 final Preconditioner preconditioner) {
122 super(checker);
123
124 this.updateFormula = updateFormula;
125 solver = lineSearchSolver;
126 this.preconditioner = preconditioner;
127 initialStep = 1.0;
128 }
129
130 /**
131 * Set the initial step used to bracket the optimum in line search.
132 * <p>
133 * The initial step is a factor with respect to the search direction,
134 * which itself is roughly related to the gradient of the function
135 * </p>
136 * @param initialStep initial step used to bracket the optimum in line search,
137 * if a non-positive value is used, the initial step is reset to its
138 * default value of 1.0
139 */
140 public void setInitialStep(final double initialStep) {
141 if (initialStep <= 0) {
142 this.initialStep = 1.0;
143 } else {
144 this.initialStep = initialStep;
145 }
146 }
147
148 /** {@inheritDoc} */
149 @Override
150 protected PointValuePair doOptimize() {
151 final ConvergenceChecker<PointValuePair> checker = getConvergenceChecker();
152 point = getStartPoint();
153 final GoalType goal = getGoalType();
154 final int n = point.length;
155 double[] r = computeObjectiveGradient(point);
156 if (goal == GoalType.MINIMIZE) {
157 for (int i = 0; i < n; ++i) {
158 r[i] = -r[i];
159 }
160 }
161
162 // Initial search direction.
163 double[] steepestDescent = preconditioner.precondition(point, r);
164 double[] searchDirection = steepestDescent.clone();
165
166 double delta = 0;
167 for (int i = 0; i < n; ++i) {
168 delta += r[i] * searchDirection[i];
169 }
170
171 PointValuePair current = null;
172 int iter = 0;
173 int maxEval = getMaxEvaluations();
174 while (true) {
175 ++iter;
176
177 final double objective = computeObjectiveValue(point);
178 PointValuePair previous = current;
179 current = new PointValuePair(point, objective);
180 if (previous != null && checker.converged(iter, previous, current)) {
181 // We have found an optimum.
182 return current;
183 }
184
185 // Find the optimal step in the search direction.
186 final UnivariateFunction lsf = new LineSearchFunction(searchDirection);
187 final double uB = findUpperBound(lsf, 0, initialStep);
188 // XXX Last parameters is set to a value close to zero in order to
189 // work around the divergence problem in the "testCircleFitting"
190 // unit test (see MATH-439).
191 final double step = solver.solve(maxEval, lsf, 0, uB, 1e-15);
192 maxEval -= solver.getEvaluations(); // Subtract used up evaluations.
193
194 // Validate new point.
195 for (int i = 0; i < point.length; ++i) {
196 point[i] += step * searchDirection[i];
197 }
198
199 r = computeObjectiveGradient(point);
200 if (goal == GoalType.MINIMIZE) {
201 for (int i = 0; i < n; ++i) {
202 r[i] = -r[i];
203 }
204 }
205
206 // Compute beta.
207 final double deltaOld = delta;
208 final double[] newSteepestDescent = preconditioner.precondition(point, r);
209 delta = 0;
210 for (int i = 0; i < n; ++i) {
211 delta += r[i] * newSteepestDescent[i];
212 }
213
214 final double beta;
215 if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) {
216 beta = delta / deltaOld;
217 } else {
218 double deltaMid = 0;
219 for (int i = 0; i < r.length; ++i) {
220 deltaMid += r[i] * steepestDescent[i];
221 }
222 beta = (delta - deltaMid) / deltaOld;
223 }
224 steepestDescent = newSteepestDescent;
225
226 // Compute conjugate search direction.
227 if (iter % n == 0 ||
228 beta < 0) {
229 // Break conjugation: reset search direction.
230 searchDirection = steepestDescent.clone();
231 } else {
232 // Compute new conjugate search direction.
233 for (int i = 0; i < n; ++i) {
234 searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
235 }
236 }
237 }
238 }
239
240 /**
241 * Find the upper bound b ensuring bracketing of a root between a and b.
242 *
243 * @param f function whose root must be bracketed.
244 * @param a lower bound of the interval.
245 * @param h initial step to try.
246 * @return b such that f(a) and f(b) have opposite signs.
247 * @throws MathIllegalStateException if no bracket can be found.
248 */
249 private double findUpperBound(final UnivariateFunction f,
250 final double a, final double h) {
251 final double yA = f.value(a);
252 double yB = yA;
253 for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) {
254 final double b = a + step;
255 yB = f.value(b);
256 if (yA * yB <= 0) {
257 return b;
258 }
259 }
260 throw new MathIllegalStateException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH);
261 }
262
263 /** Default identity preconditioner. */
264 public static class IdentityPreconditioner implements Preconditioner {
265
266 /** {@inheritDoc} */
267 public double[] precondition(double[] variables, double[] r) {
268 return r.clone();
269 }
270 }
271
272 /** Internal class for line search.
273 * <p>
274 * The function represented by this class is the dot product of
275 * the objective function gradient and the search direction. Its
276 * value is zero when the gradient is orthogonal to the search
277 * direction, i.e. when the objective function value is a local
278 * extremum along the search direction.
279 * </p>
280 */
281 private class LineSearchFunction implements UnivariateFunction {
282 /** Search direction. */
283 private final double[] searchDirection;
284
285 /** Simple constructor.
286 * @param searchDirection search direction
287 */
288 public LineSearchFunction(final double[] searchDirection) {
289 this.searchDirection = searchDirection;
290 }
291
292 /** {@inheritDoc} */
293 public double value(double x) {
294 // current point in the search direction
295 final double[] shiftedPoint = point.clone();
296 for (int i = 0; i < shiftedPoint.length; ++i) {
297 shiftedPoint[i] += x * searchDirection[i];
298 }
299
300 // gradient of the objective function
301 final double[] gradient = computeObjectiveGradient(shiftedPoint);
302
303 // dot product with the search direction
304 double dotProduct = 0;
305 for (int i = 0; i < gradient.length; ++i) {
306 dotProduct += gradient[i] * searchDirection[i];
307 }
308
309 return dotProduct;
310 }
311 }
312 }