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   * @since 2.0
45   */
46  public class NonLinearConjugateGradientOptimizer
47      extends GradientMultivariateOptimizer {
48      /** Update formula for the beta parameter. */
49      private final Formula updateFormula;
50      /** Preconditioner (may be null). */
51      private final Preconditioner preconditioner;
52      /** Line search algorithm. */
53      private final LineSearch line;
54  
55      /**
56       * Available choices of update formulas for the updating the parameter
57       * that is used to compute the successive conjugate search directions.
58       * For non-linear conjugate gradients, there are
59       * two formulas:
60       * <ul>
61       *   <li>Fletcher-Reeves formula</li>
62       *   <li>Polak-Ribière formula</li>
63       * </ul>
64       *
65       * On the one hand, the Fletcher-Reeves formula is guaranteed to converge
66       * if the start point is close enough of the optimum whether the
67       * Polak-Ribière formula may not converge in rare cases. On the
68       * other hand, the Polak-Ribière formula is often faster when it
69       * does converge. Polak-Ribière is often used.
70       *
71       * @since 2.0
72       */
73      public enum Formula {
74          /** Fletcher-Reeves formula. */
75          FLETCHER_REEVES,
76          /** Polak-Ribière formula. */
77          POLAK_RIBIERE
78      }
79  
80      /**
81       * The initial step is a factor with respect to the search direction
82       * (which itself is roughly related to the gradient of the function).
83       * <br/>
84       * It is used to find an interval that brackets the optimum in line
85       * search.
86       *
87       * @since 3.1
88       * @deprecated As of v3.3, this class is not used anymore.
89       * This setting is replaced by the {@code initialBracketingRange}
90       * argument to the new constructors.
91       */
92      @Deprecated
93      public static class BracketingStep implements OptimizationData {
94          /** Initial step. */
95          private final double initialStep;
96  
97          /**
98           * @param step Initial step for the bracket search.
99           */
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 &beta; 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 &beta; 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 &beta; 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 &beta; 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 &beta; 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 }