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 1422230 2012-12-15 12:11:13Z erans $
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) {
181 if (checker.converged(iter, previous, current)) {
182 // We have found an optimum.
183 return current;
184 }
185 }
186
187 // Find the optimal step in the search direction.
188 final UnivariateFunction lsf = new LineSearchFunction(searchDirection);
189 final double uB = findUpperBound(lsf, 0, initialStep);
190 // XXX Last parameters is set to a value close to zero in order to
191 // work around the divergence problem in the "testCircleFitting"
192 // unit test (see MATH-439).
193 final double step = solver.solve(maxEval, lsf, 0, uB, 1e-15);
194 maxEval -= solver.getEvaluations(); // Subtract used up evaluations.
195
196 // Validate new point.
197 for (int i = 0; i < point.length; ++i) {
198 point[i] += step * searchDirection[i];
199 }
200
201 r = computeObjectiveGradient(point);
202 if (goal == GoalType.MINIMIZE) {
203 for (int i = 0; i < n; ++i) {
204 r[i] = -r[i];
205 }
206 }
207
208 // Compute beta.
209 final double deltaOld = delta;
210 final double[] newSteepestDescent = preconditioner.precondition(point, r);
211 delta = 0;
212 for (int i = 0; i < n; ++i) {
213 delta += r[i] * newSteepestDescent[i];
214 }
215
216 final double beta;
217 if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) {
218 beta = delta / deltaOld;
219 } else {
220 double deltaMid = 0;
221 for (int i = 0; i < r.length; ++i) {
222 deltaMid += r[i] * steepestDescent[i];
223 }
224 beta = (delta - deltaMid) / deltaOld;
225 }
226 steepestDescent = newSteepestDescent;
227
228 // Compute conjugate search direction.
229 if (iter % n == 0 ||
230 beta < 0) {
231 // Break conjugation: reset search direction.
232 searchDirection = steepestDescent.clone();
233 } else {
234 // Compute new conjugate search direction.
235 for (int i = 0; i < n; ++i) {
236 searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
237 }
238 }
239 }
240 }
241
242 /**
243 * Find the upper bound b ensuring bracketing of a root between a and b.
244 *
245 * @param f function whose root must be bracketed.
246 * @param a lower bound of the interval.
247 * @param h initial step to try.
248 * @return b such that f(a) and f(b) have opposite signs.
249 * @throws MathIllegalStateException if no bracket can be found.
250 */
251 private double findUpperBound(final UnivariateFunction f,
252 final double a, final double h) {
253 final double yA = f.value(a);
254 double yB = yA;
255 for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) {
256 final double b = a + step;
257 yB = f.value(b);
258 if (yA * yB <= 0) {
259 return b;
260 }
261 }
262 throw new MathIllegalStateException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH);
263 }
264
265 /** Default identity preconditioner. */
266 public static class IdentityPreconditioner implements Preconditioner {
267
268 /** {@inheritDoc} */
269 public double[] precondition(double[] variables, double[] r) {
270 return r.clone();
271 }
272 }
273
274 /** Internal class for line search.
275 * <p>
276 * The function represented by this class is the dot product of
277 * the objective function gradient and the search direction. Its
278 * value is zero when the gradient is orthogonal to the search
279 * direction, i.e. when the objective function value is a local
280 * extremum along the search direction.
281 * </p>
282 */
283 private class LineSearchFunction implements UnivariateFunction {
284 /** Search direction. */
285 private final double[] searchDirection;
286
287 /** Simple constructor.
288 * @param searchDirection search direction
289 */
290 public LineSearchFunction(final double[] searchDirection) {
291 this.searchDirection = searchDirection;
292 }
293
294 /** {@inheritDoc} */
295 public double value(double x) {
296 // current point in the search direction
297 final double[] shiftedPoint = point.clone();
298 for (int i = 0; i < shiftedPoint.length; ++i) {
299 shiftedPoint[i] += x * searchDirection[i];
300 }
301
302 // gradient of the objective function
303 final double[] gradient = computeObjectiveGradient(shiftedPoint);
304
305 // dot product with the search direction
306 double dotProduct = 0;
307 for (int i = 0; i < gradient.length; ++i) {
308 dotProduct += gradient[i] * searchDirection[i];
309 }
310
311 return dotProduct;
312 }
313 }
314 }