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.math.optimization.general;
19
20 import org.apache.commons.math.ConvergenceException;
21 import org.apache.commons.math.FunctionEvaluationException;
22 import org.apache.commons.math.analysis.UnivariateRealFunction;
23 import org.apache.commons.math.analysis.solvers.BrentSolver;
24 import org.apache.commons.math.analysis.solvers.UnivariateRealSolver;
25 import org.apache.commons.math.optimization.GoalType;
26 import org.apache.commons.math.optimization.OptimizationException;
27 import org.apache.commons.math.optimization.DifferentiableMultivariateRealOptimizer;
28 import org.apache.commons.math.optimization.RealPointValuePair;
29 import org.apache.commons.math.optimization.SimpleVectorialValueChecker;
30
31 /**
32 * Non-linear conjugate gradient optimizer.
33 * <p>
34 * This class supports both the Fletcher-Reeves and the Polak-Ribière
35 * update formulas for the conjugate search directions. It also supports
36 * optional preconditioning.
37 * </p>
38 *
39 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
40 * @since 2.0
41 *
42 */
43
44 public class NonLinearConjugateGradientOptimizer
45 extends AbstractScalarDifferentiableOptimizer
46 implements DifferentiableMultivariateRealOptimizer {
47
48 /** Update formula for the beta parameter. */
49 private final ConjugateGradientFormula updateFormula;
50
51 /** Preconditioner (may be null). */
52 private Preconditioner preconditioner;
53
54 /** solver to use in the line search (may be null). */
55 private UnivariateRealSolver solver;
56
57 /** Initial step used to bracket the optimum in line search. */
58 private double initialStep;
59
60 /** Simple constructor with default settings.
61 * <p>The convergence check is set to a {@link SimpleVectorialValueChecker}
62 * and the maximal number of iterations is set to
63 * {@link AbstractScalarDifferentiableOptimizer#DEFAULT_MAX_ITERATIONS}.
64 * @param updateFormula formula to use for updating the β parameter,
65 * must be one of {@link ConjugateGradientFormula#FLETCHER_REEVES} or {@link
66 * ConjugateGradientFormula#POLAK_RIBIERE}
67 */
68 public NonLinearConjugateGradientOptimizer(final ConjugateGradientFormula updateFormula) {
69 this.updateFormula = updateFormula;
70 preconditioner = null;
71 solver = null;
72 initialStep = 1.0;
73 }
74
75 /**
76 * Set the preconditioner.
77 * @param preconditioner preconditioner to use for next optimization,
78 * may be null to remove an already registered preconditioner
79 */
80 public void setPreconditioner(final Preconditioner preconditioner) {
81 this.preconditioner = preconditioner;
82 }
83
84 /**
85 * Set the solver to use during line search.
86 * @param solver solver to use during line search, may be null
87 * to remove an already registered solver and fall back to the
88 * default {@link BrentSolver Brent solver}.
89 */
90 public void setLineSearchSolver(final UnivariateRealSolver solver) {
91 this.solver = solver;
92 }
93
94 /**
95 * Set the initial step used to bracket the optimum in line search.
96 * <p>
97 * The initial step is a factor with respect to the search direction,
98 * which itself is roughly related to the gradient of the function
99 * </p>
100 * @param initialStep initial step used to bracket the optimum in line search,
101 * if a non-positive value is used, the initial step is reset to its
102 * default value of 1.0
103 */
104 public void setInitialStep(final double initialStep) {
105 if (initialStep <= 0) {
106 this.initialStep = 1.0;
107 } else {
108 this.initialStep = initialStep;
109 }
110 }
111
112 /** {@inheritDoc} */
113 @Override
114 protected RealPointValuePair doOptimize()
115 throws FunctionEvaluationException, OptimizationException, IllegalArgumentException {
116 try {
117
118 // initialization
119 if (preconditioner == null) {
120 preconditioner = new IdentityPreconditioner();
121 }
122 if (solver == null) {
123 solver = new BrentSolver();
124 }
125 final int n = point.length;
126 double[] r = computeObjectiveGradient(point);
127 if (goalType == GoalType.MINIMIZE) {
128 for (int i = 0; i < n; ++i) {
129 r[i] = -r[i];
130 }
131 }
132
133 // initial search direction
134 double[] steepestDescent = preconditioner.precondition(point, r);
135 double[] searchDirection = steepestDescent.clone();
136
137 double delta = 0;
138 for (int i = 0; i < n; ++i) {
139 delta += r[i] * searchDirection[i];
140 }
141
142 RealPointValuePair current = null;
143 while (true) {
144
145 final double objective = computeObjectiveValue(point);
146 RealPointValuePair previous = current;
147 current = new RealPointValuePair(point, objective);
148 if (previous != null) {
149 if (checker.converged(getIterations(), previous, current)) {
150 // we have found an optimum
151 return current;
152 }
153 }
154
155 incrementIterationsCounter();
156
157 double dTd = 0;
158 for (final double di : searchDirection) {
159 dTd += di * di;
160 }
161
162 // find the optimal step in the search direction
163 final UnivariateRealFunction lsf = new LineSearchFunction(searchDirection);
164 final double step = solver.solve(lsf, 0, findUpperBound(lsf, 0, initialStep));
165
166 // validate new point
167 for (int i = 0; i < point.length; ++i) {
168 point[i] += step * searchDirection[i];
169 }
170 r = computeObjectiveGradient(point);
171 if (goalType == GoalType.MINIMIZE) {
172 for (int i = 0; i < n; ++i) {
173 r[i] = -r[i];
174 }
175 }
176
177 // compute beta
178 final double deltaOld = delta;
179 final double[] newSteepestDescent = preconditioner.precondition(point, r);
180 delta = 0;
181 for (int i = 0; i < n; ++i) {
182 delta += r[i] * newSteepestDescent[i];
183 }
184
185 final double beta;
186 if (updateFormula == ConjugateGradientFormula.FLETCHER_REEVES) {
187 beta = delta / deltaOld;
188 } else {
189 double deltaMid = 0;
190 for (int i = 0; i < r.length; ++i) {
191 deltaMid += r[i] * steepestDescent[i];
192 }
193 beta = (delta - deltaMid) / deltaOld;
194 }
195 steepestDescent = newSteepestDescent;
196
197 // compute conjugate search direction
198 if ((getIterations() % n == 0) || (beta < 0)) {
199 // break conjugation: reset search direction
200 searchDirection = steepestDescent.clone();
201 } else {
202 // compute new conjugate search direction
203 for (int i = 0; i < n; ++i) {
204 searchDirection[i] = steepestDescent[i] + beta * searchDirection[i];
205 }
206 }
207
208 }
209
210 } catch (ConvergenceException ce) {
211 throw new OptimizationException(ce);
212 }
213 }
214
215 /**
216 * Find the upper bound b ensuring bracketing of a root between a and b
217 * @param f function whose root must be bracketed
218 * @param a lower bound of the interval
219 * @param h initial step to try
220 * @return b such that f(a) and f(b) have opposite signs
221 * @exception FunctionEvaluationException if the function cannot be computed
222 * @exception OptimizationException if no bracket can be found
223 */
224 private double findUpperBound(final UnivariateRealFunction f,
225 final double a, final double h)
226 throws FunctionEvaluationException, OptimizationException {
227 final double yA = f.value(a);
228 double yB = yA;
229 for (double step = h; step < Double.MAX_VALUE; step *= Math.max(2, yA / yB)) {
230 final double b = a + step;
231 yB = f.value(b);
232 if (yA * yB <= 0) {
233 return b;
234 }
235 }
236 throw new OptimizationException("unable to bracket optimum in line search");
237 }
238
239 /** Default identity preconditioner. */
240 private static class IdentityPreconditioner implements Preconditioner {
241
242 /** {@inheritDoc} */
243 public double[] precondition(double[] variables, double[] r) {
244 return r.clone();
245 }
246
247 }
248
249 /** Internal class for line search.
250 * <p>
251 * The function represented by this class is the dot product of
252 * the objective function gradient and the search direction. Its
253 * value is zero when the gradient is orthogonal to the search
254 * direction, i.e. when the objective function value is a local
255 * extremum along the search direction.
256 * </p>
257 */
258 private class LineSearchFunction implements UnivariateRealFunction {
259 /** Search direction. */
260 private final double[] searchDirection;
261
262 /** Simple constructor.
263 * @param searchDirection search direction
264 */
265 public LineSearchFunction(final double[] searchDirection) {
266 this.searchDirection = searchDirection;
267 }
268
269 /** {@inheritDoc} */
270 public double value(double x) throws FunctionEvaluationException {
271
272 // current point in the search direction
273 final double[] shiftedPoint = point.clone();
274 for (int i = 0; i < shiftedPoint.length; ++i) {
275 shiftedPoint[i] += x * searchDirection[i];
276 }
277
278 // gradient of the objective function
279 final double[] gradient = computeObjectiveGradient(shiftedPoint);
280
281 // dot product with the search direction
282 double dotProduct = 0;
283 for (int i = 0; i < gradient.length; ++i) {
284 dotProduct += gradient[i] * searchDirection[i];
285 }
286
287 return dotProduct;
288
289 }
290
291 }
292
293 }