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.optimization.general;
019
020import org.apache.commons.math3.exception.MathIllegalStateException;
021import org.apache.commons.math3.analysis.UnivariateFunction;
022import org.apache.commons.math3.analysis.solvers.BrentSolver;
023import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
024import org.apache.commons.math3.exception.util.LocalizedFormats;
025import org.apache.commons.math3.optimization.GoalType;
026import org.apache.commons.math3.optimization.PointValuePair;
027import org.apache.commons.math3.optimization.SimpleValueChecker;
028import org.apache.commons.math3.optimization.ConvergenceChecker;
029import org.apache.commons.math3.util.FastMath;
030
031/**
032 * Non-linear conjugate gradient optimizer.
033 * <p>
034 * This class supports both the Fletcher-Reeves and the Polak-Ribi&egrave;re
035 * update formulas for the conjugate search directions. It also supports
036 * optional preconditioning.
037 * </p>
038 *
039 * @version $Id: NonLinearConjugateGradientOptimizer.java 1462503 2013-03-29 15:48:27Z luc $
040 * @deprecated As of 3.1 (to be removed in 4.0).
041 * @since 2.0
042 *
043 */
044@Deprecated
045public class NonLinearConjugateGradientOptimizer
046    extends AbstractScalarDifferentiableOptimizer {
047    /** Update formula for the beta parameter. */
048    private final ConjugateGradientFormula updateFormula;
049    /** Preconditioner (may be null). */
050    private final Preconditioner preconditioner;
051    /** solver to use in the line search (may be null). */
052    private final UnivariateSolver solver;
053    /** Initial step used to bracket the optimum in line search. */
054    private double initialStep;
055    /** Current point. */
056    private double[] point;
057
058    /**
059     * Constructor with default {@link SimpleValueChecker checker},
060     * {@link BrentSolver line search solver} and
061     * {@link IdentityPreconditioner preconditioner}.
062     *
063     * @param updateFormula formula to use for updating the &beta; parameter,
064     * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
065     * ConjugateGradientFormula#POLAK_RIBIERE}.
066     * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()}
067     */
068    @Deprecated
069    public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) {
070        this(updateFormula,
071             new SimpleValueChecker());
072    }
073
074    /**
075     * Constructor with default {@link BrentSolver line search solver} and
076     * {@link IdentityPreconditioner preconditioner}.
077     *
078     * @param updateFormula formula to use for updating the &beta; parameter,
079     * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
080     * ConjugateGradientFormula#POLAK_RIBIERE}.
081     * @param checker Convergence checker.
082     */
083    public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
084                                               ConvergenceChecker<PointValuePair> checker) {
085        this(updateFormula,
086             checker,
087             new BrentSolver(),
088             new IdentityPreconditioner());
089    }
090
091
092    /**
093     * Constructor with default {@link IdentityPreconditioner preconditioner}.
094     *
095     * @param updateFormula formula to use for updating the &beta; parameter,
096     * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
097     * ConjugateGradientFormula#POLAK_RIBIERE}.
098     * @param checker Convergence checker.
099     * @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}