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 * @deprecated As of 3.1 (to be removed in 4.0).
040 * @since 2.0
041 *
042 */
043@Deprecated
044public class NonLinearConjugateGradientOptimizer
045    extends AbstractScalarDifferentiableOptimizer {
046    /** Update formula for the beta parameter. */
047    private final ConjugateGradientFormula updateFormula;
048    /** Preconditioner (may be null). */
049    private final Preconditioner preconditioner;
050    /** solver to use in the line search (may be null). */
051    private final UnivariateSolver solver;
052    /** Initial step used to bracket the optimum in line search. */
053    private double initialStep;
054    /** Current point. */
055    private double[] point;
056
057    /**
058     * Constructor with default {@link SimpleValueChecker checker},
059     * {@link BrentSolver line search solver} and
060     * {@link IdentityPreconditioner preconditioner}.
061     *
062     * @param updateFormula formula to use for updating the &beta; parameter,
063     * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
064     * ConjugateGradientFormula#POLAK_RIBIERE}.
065     * @deprecated See {@link SimpleValueChecker#SimpleValueChecker()}
066     */
067    @Deprecated
068    public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) {
069        this(updateFormula,
070             new SimpleValueChecker());
071    }
072
073    /**
074     * Constructor with default {@link BrentSolver line search solver} and
075     * {@link IdentityPreconditioner preconditioner}.
076     *
077     * @param updateFormula formula to use for updating the &beta; parameter,
078     * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
079     * ConjugateGradientFormula#POLAK_RIBIERE}.
080     * @param checker Convergence checker.
081     */
082    public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula,
083                                               ConvergenceChecker<PointValuePair> checker) {
084        this(updateFormula,
085             checker,
086             new BrentSolver(),
087             new IdentityPreconditioner());
088    }
089
090
091    /**
092     * Constructor with default {@link IdentityPreconditioner preconditioner}.
093     *
094     * @param updateFormula formula to use for updating the &beta; parameter,
095     * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
096     * ConjugateGradientFormula#POLAK_RIBIERE}.
097     * @param checker Convergence checker.
098     * @param lineSearchSolver Solver to use during line search.
099     */
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}