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
018package org.apache.commons.math3.optim.nonlinear.scalar.gradient;
019
020import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
021import org.apache.commons.math3.exception.MathInternalError;
022import org.apache.commons.math3.exception.TooManyEvaluationsException;
023import org.apache.commons.math3.exception.MathUnsupportedOperationException;
024import org.apache.commons.math3.exception.util.LocalizedFormats;
025import org.apache.commons.math3.optim.OptimizationData;
026import org.apache.commons.math3.optim.PointValuePair;
027import org.apache.commons.math3.optim.ConvergenceChecker;
028import org.apache.commons.math3.optim.nonlinear.scalar.GoalType;
029import org.apache.commons.math3.optim.nonlinear.scalar.GradientMultivariateOptimizer;
030import org.apache.commons.math3.optim.nonlinear.scalar.LineSearch;
031
032
033/**
034 * Non-linear conjugate gradient optimizer.
035 * <br/>
036 * This class supports both the Fletcher-Reeves and the Polak-Ribi?re
037 * update formulas for the conjugate search directions.
038 * It also supports optional preconditioning.
039 * <br/>
040 * Constraints are not supported: the call to
041 * {@link #optimize(OptimizationData[]) optimize} will throw
042 * {@link MathUnsupportedOperationException} if bounds are passed to it.
043 *
044 * @since 2.0
045 */
046public class NonLinearConjugateGradientOptimizer
047    extends GradientMultivariateOptimizer {
048    /** Update formula for the beta parameter. */
049    private final Formula updateFormula;
050    /** Preconditioner (may be null). */
051    private final Preconditioner preconditioner;
052    /** Line search algorithm. */
053    private final LineSearch line;
054
055    /**
056     * Available choices of update formulas for the updating the parameter
057     * that is used to compute the successive conjugate search directions.
058     * For non-linear conjugate gradients, there are
059     * two formulas:
060     * <ul>
061     *   <li>Fletcher-Reeves formula</li>
062     *   <li>Polak-Ribi?re formula</li>
063     * </ul>
064     *
065     * On the one hand, the Fletcher-Reeves formula is guaranteed to converge
066     * if the start point is close enough of the optimum whether the
067     * Polak-Ribi?re formula may not converge in rare cases. On the
068     * other hand, the Polak-Ribi?re formula is often faster when it
069     * does converge. Polak-Ribi?re is often used.
070     *
071     * @since 2.0
072     */
073    public static enum Formula {
074        /** Fletcher-Reeves formula. */
075        FLETCHER_REEVES,
076        /** Polak-Ribi?re formula. */
077        POLAK_RIBIERE
078    }
079
080    /**
081     * The initial step is a factor with respect to the search direction
082     * (which itself is roughly related to the gradient of the function).
083     * <br/>
084     * It is used to find an interval that brackets the optimum in line
085     * search.
086     *
087     * @since 3.1
088     * @deprecated As of v3.3, this class is not used anymore.
089     * This setting is replaced by the {@code initialBracketingRange}
090     * argument to the new constructors.
091     */
092    @Deprecated
093    public static class BracketingStep implements OptimizationData {
094        /** Initial step. */
095        private final double initialStep;
096
097        /**
098         * @param step Initial step for the bracket search.
099         */
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}