    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
    018    package org.apache.commons.math3.optimization.general;
    019
    020    import org.apache.commons.math3.exception.MathIllegalStateException;
    021    import org.apache.commons.math3.analysis.UnivariateFunction;
    022    import org.apache.commons.math3.analysis.solvers.BrentSolver;
    023    import org.apache.commons.math3.analysis.solvers.UnivariateSolver;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025    import org.apache.commons.math3.optimization.GoalType;
    026    import org.apache.commons.math3.optimization.PointValuePair;
    027    import org.apache.commons.math3.optimization.SimpleValueChecker;
    028    import org.apache.commons.math3.optimization.ConvergenceChecker;
    029    import 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 1422230 2012-12-15 12:11:13Z erans $
    040     * @deprecated As of 3.1 (to be removed in 4.0).
    041     * @since 2.0
    042     *
    043     */
    044    @Deprecated
    045    public 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) {
    181                    if (checker.converged(iter, previous, current)) {
    182                        // We have found an optimum.
    183                        return current;
    184                    }
    185                }
    186
    187                // Find the optimal step in the search direction.
    188                final UnivariateFunction lsf = new LineSearchFunction(searchDirection);
    189                final double uB = findUpperBound(lsf, 0, initialStep);
    190                // XXX Last parameters is set to a value close to zero in order to
    191                // work around the divergence problem in the "testCircleFitting"
    192                // unit test (see MATH-439).
    193                final double step = solver.solve(maxEval, lsf, 0, uB, 1e-15);
    194                maxEval -= solver.getEvaluations(); // Subtract used up evaluations.
    195
    196                // Validate new point.
    197                for (int i = 0; i < point.length; ++i) {
    198                    point[i] += step * searchDirection[i];
    199                }
    200
    201                r = computeObjectiveGradient(point);
    202                if (goal == GoalType.MINIMIZE) {
    203                    for (int i = 0; i < n; ++i) {
    204                        r[i] = -r[i];
    205                    }
    206                }
    207
    208                // Compute beta.
    209                final double deltaOld = delta;
    210                final double[] newSteepestDescent = preconditioner.precondition(point, r);
    211                delta = 0;
    212                for (int i = 0; i < n; ++i) {
    213                    delta += r[i] * newSteepestDescent[i];
    214                }
    215
    216                final double beta;
    217                if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) {
    218                    beta = delta / deltaOld;
    219                } else {
    220                    double deltaMid = 0;
    221                    for (int i = 0; i < r.length; ++i) {
    222                        deltaMid += r[i] * steepestDescent[i];
    223                    }
    224                    beta = (delta - deltaMid) / deltaOld;
    225                }
    226                steepestDescent = newSteepestDescent;
    227
    228                // Compute conjugate search direction.
    229                if (iter % n == 0 ||
    230                    beta < 0) {
    231                    // Break conjugation: reset search direction.
    232                    searchDirection = steepestDescent.clone();
    233                } else {
    234                    // Compute new conjugate search direction.
    235                    for (int i = 0; i < n; ++i) {
    236                        searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
    237                    }
    238                }
    239            }
    240        }
    241
    242        /**
    243         * Find the upper bound b ensuring bracketing of a root between a and b.
    244         *
    245         * @param f function whose root must be bracketed.
    246         * @param a lower bound of the interval.
    247         * @param h initial step to try.
    248         * @return b such that f(a) and f(b) have opposite signs.
    249         * @throws MathIllegalStateException if no bracket can be found.
    250         */
    251        private double findUpperBound(final UnivariateFunction f,
    252                                      final double a, final double h) {
    253            final double yA = f.value(a);
    254            double yB = yA;
    255            for (double step = h; step < Double.MAX_VALUE; step *= FastMath.max(2, yA / yB)) {
    256                final double b = a + step;
    257                yB = f.value(b);
    258                if (yA * yB <= 0) {
    259                    return b;
    260                }
    261            }
    262            throw new MathIllegalStateException(LocalizedFormats.UNABLE_TO_BRACKET_OPTIMUM_IN_LINE_SEARCH);
    263        }
    264
    265        /** Default identity preconditioner. */
    266        public static class IdentityPreconditioner implements Preconditioner {
    267
    268            /** {@inheritDoc} */
    269            public double[] precondition(double[] variables, double[] r) {
    270                return r.clone();
    271            }
    272        }
    273
    274        /** Internal class for line search.
    275         * <p>
    276         * The function represented by this class is the dot product of
    277         * the objective function gradient and the search direction. Its
    278         * value is zero when the gradient is orthogonal to the search
    279         * direction, i.e. when the objective function value is a local
    280         * extremum along the search direction.
    281         * </p>
    282         */
    283        private class LineSearchFunction implements UnivariateFunction {
    284            /** Search direction. */
    285            private final double[] searchDirection;
    286
    287            /** Simple constructor.
    288             * @param searchDirection search direction
    289             */
    290            public LineSearchFunction(final double[] searchDirection) {
    291                this.searchDirection = searchDirection;
    292            }
    293
    294            /** {@inheritDoc} */
    295            public double value(double x) {
    296                // current point in the search direction
    297                final double[] shiftedPoint = point.clone();
    298                for (int i = 0; i < shiftedPoint.length; ++i) {
    299                    shiftedPoint[i] += x * searchDirection[i];
    300                }
    301
    302                // gradient of the objective function
    303                final double[] gradient = computeObjectiveGradient(shiftedPoint);
    304
    305                // dot product with the search direction
    306                double dotProduct = 0;
    307                for (int i = 0; i < gradient.length; ++i) {
    308                    dotProduct += gradient[i] * searchDirection[i];
    309                }
    310
    311                return dotProduct;
    312            }
    313        }
    314    }