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