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   * @version $Id: NonLinearConjugateGradientOptimizer.java 1462503 2013-03-29 15:48:27Z luc $
40   * @deprecated As of 3.1 (to be removed in 4.0).
41   * @since 2.0
42   *
43   */
44  @Deprecated
45  public class NonLinearConjugateGradientOptimizer
46      extends AbstractScalarDifferentiableOptimizer {
47      /** Update formula for the beta parameter. */
48      private final ConjugateGradientFormula updateFormula;
49      /** Preconditioner (may be null). */
50      private final Preconditioner preconditioner;
51      /** solver to use in the line search (may be null). */
52      private final UnivariateSolver solver;
53      /** Initial step used to bracket the optimum in line search. */
54      private double initialStep;
55      /** Current point. */
56      private double[] point;
57  
58      /**
59       * Constructor with default {@link SimpleValueChecker checker},
60       * {@link BrentSolver line search solver} and
61       * {@link IdentityPreconditioner preconditioner}.
62       *
63       * @param updateFormula formula to use for updating the &beta; parameter,
64       * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
65       * ConjugateGradientFormula#POLAK_RIBIERE}.
66       * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()}
67       */
68      @Deprecated
69      public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) {
70          this(updateFormula,
71               new SimpleValueChecker());
72      }
73  
74      /**
75       * Constructor with default {@link BrentSolver line search solver} and
76       * {@link IdentityPreconditioner preconditioner}.
77       *
78       * @param updateFormula formula to use for updating the &beta; parameter,
79       * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
80       * ConjugateGradientFormula#POLAK_RIBIERE}.
81       * @param checker Convergence checker.
82       */
83      public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
84                                                 ConvergenceChecker<PointValuePair> checker) {
85          this(updateFormula,
86               checker,
87               new BrentSolver(),
88               new IdentityPreconditioner());
89      }
90  
91  
92      /**
93       * Constructor with default {@link IdentityPreconditioner preconditioner}.
94       *
95       * @param updateFormula formula to use for updating the &beta; parameter,
96       * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
97       * ConjugateGradientFormula#POLAK_RIBIERE}.
98       * @param checker Convergence checker.
99       * @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 &beta; 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 }