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.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&egrave;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 &beta; 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
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 &beta; parameter,
78       * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
79       * ConjugateGradientFormula#POLAK_RIBIERE}.
80       * @param checker Convergence checker.
81       */
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 &beta; 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       */
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 &beta; 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      */
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      */
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
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      */
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} */
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          */
287         public LineSearchFunction(final double[] searchDirection) {
288             this.searchDirection = searchDirection;
289         }
290 
291         /** {@inheritDoc} */
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 }