1 /*
2 * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.optimization;
17
18 import com.irurueta.numerical.EvaluationException;
19 import com.irurueta.numerical.GradientFunctionEvaluatorListener;
20 import com.irurueta.numerical.LockedException;
21 import com.irurueta.numerical.MultiDimensionFunctionEvaluatorListener;
22 import com.irurueta.numerical.NotAvailableException;
23 import com.irurueta.numerical.NotReadyException;
24
25 /**
26 * This class searches for a multi dimension function local minimum.
27 * The local minimum is searched by starting the algorithm at a start point
28 * and a given direction which should be close and point to the local minimum to
29 * be found to achieve the best accuracy with the lowest number of iterations.
30 * NOTE: this algorithm might not have proper convergence in some situations,
31 * but it is ensured to provide faster convergence than other algorithms such
32 * as Brent because gradient information is also provided.
33 * The implementation of this class is based on Numerical Recipes 3rd ed.
34 * Section 10.8 page 515.
35 */
36 public class ConjugateGradientMultiOptimizer extends LineMultiOptimizer {
37 /**
38 * Constant defining default tolerance or accuracy to be achieved on the
39 * minimum being estimated by this class.
40 */
41 public static final double DEFAULT_TOLERANCE = 3e-8;
42
43 /**
44 * Minimum allowed tolerance value.
45 */
46 public static final double MIN_TOLERANCE = 0.0;
47
48 /**
49 * Maximum allowed iterations.
50 */
51 public static final int ITMAX = 200;
52
53 /**
54 * Constant defining a value to be considered as machine precision.
55 */
56 public static final double EPS = 1e-18;
57
58 /**
59 * Convergence criterion for the zero gradient test.
60 */
61 public static final double GTOL = 1e-8;
62
63 /**
64 * Defines whether Polak-Ribiere is used if true, otherwise Fletcher-Reeves
65 * will be used.
66 */
67 public static final boolean DEFAULT_USE_POLAK_RIBIERE = true;
68
69 /**
70 * The fractional tolerance in the function value such that failure to
71 * decrease by more than this amount on one iteration signals done-ness.
72 */
73 private double tolerance;
74
75 /**
76 * Member contains number of iterations that were needed to estimate a
77 * minimum.
78 */
79 private int iter;
80
81 /**
82 * Value of the function at the minimum.
83 */
84 private double fret;
85
86 /**
87 * Listener to obtain gradient values for the multi dimension function being
88 * evaluated.
89 * If the gradient is unknown (e.g. doesn't have a closed expression), the
90 * provided listener could use a GradientEstimator to obtain one.
91 */
92 private GradientFunctionEvaluatorListener gradientListener;
93
94 /**
95 * Boolean indicating whether Polak-Ribiere method is used if true,
96 * otherwise Fletcher-Reeves will be used.
97 */
98 private boolean usePolakRibiere;
99
100 /**
101 * Empty constructor.
102 */
103 public ConjugateGradientMultiOptimizer() {
104 super();
105 tolerance = DEFAULT_TOLERANCE;
106 usePolakRibiere = DEFAULT_USE_POLAK_RIBIERE;
107 gradientListener = null;
108 iter = 0;
109 }
110
111 /**
112 * Constructor.
113 *
114 * @param listener Listener to evaluate a multi-dimension function.
115 * @param gradientListener Listener to obtain gradient value for the
116 * multi-dimension function being evaluated.
117 * @param point Start point where algorithm will be started. Start point
118 * should be close to the local minimum to be found. Provided array must
119 * have a length equal to the number of dimensions of the function being
120 * evaluated, otherwise and exception will be raised when searching for the
121 * minimum.
122 * @param direction Direction to start looking for a minimum. Provided array
123 * must have the same length as the number of dimensions of the function
124 * being evaluated. Provided direction is considered as a vector pointing
125 * to the minimum to be found.
126 * @param tolerance Tolerance or accuracy to be expected on estimated local
127 * minimum.
128 * @param usePolakRibiere True if Polak-Ribiere method is used, otherwise
129 * Fletcher-Reeves will be used.
130 * @throws IllegalArgumentException Raised if provided point and direction
131 * don't have the same length or if provided tolerance is negative.
132 */
133 public ConjugateGradientMultiOptimizer(
134 final MultiDimensionFunctionEvaluatorListener listener,
135 final GradientFunctionEvaluatorListener gradientListener, final double[] point, final double[] direction,
136 final double tolerance, final boolean usePolakRibiere) {
137 super(listener, point, direction);
138 internalSetTolerance(tolerance);
139 this.usePolakRibiere = usePolakRibiere;
140 this.gradientListener = gradientListener;
141 iter = 0;
142 }
143
144 /**
145 * Constructor.
146 *
147 * @param listener Listener to evaluate a multidimensional function.
148 * @param gradientListener Listener to obtain gradient value for the
149 * multidimensional function being evaluated.
150 * @param point Start point where algorithm will be started. Start point
151 * should be close to the local minimum to be found. Provided array must
152 * have a length equal to the number of dimensions of the function being
153 * evaluated, otherwise and exception will be raised when searching for the
154 * minimum.
155 * @param tolerance Tolerance or accuracy to be expected on estimated local
156 * minimum.
157 * @param usePolakRibiere True if Polak-Ribiere method is used, otherwise
158 * Fletcher-Reeves will be used.
159 * @throws IllegalArgumentException Raised if tolerance is negative.
160 */
161 public ConjugateGradientMultiOptimizer(
162 final MultiDimensionFunctionEvaluatorListener listener,
163 final GradientFunctionEvaluatorListener gradientListener, final double[] point, final double tolerance,
164 final boolean usePolakRibiere) {
165 super(listener);
166 internalSetStartPoint(point);
167 internalSetTolerance(tolerance);
168 this.usePolakRibiere = usePolakRibiere;
169 this.gradientListener = gradientListener;
170 iter = 0;
171 }
172
173 /**
174 * This function estimates a function minimum.
175 * Implementations of this class will usually search a local minimum close
176 * to a start point and will start looking into provided start direction.
177 * Minimization of a function f. Input consists of an initial starting point
178 * p. The initial matrix xi, whose columns contain the initial set of
179 * directions, is set to the identity. Returned is the best point found, at
180 * which point fret is the minimum function value and iter is the number of
181 * iterations taken.
182 *
183 * @throws LockedException Raised if this instance is locked, because
184 * estimation is being computed.
185 * @throws NotReadyException Raised if this instance is not ready, because
186 * a listener, a gradient listener and a start point haven't been provided.
187 * @throws OptimizationException Raised if the algorithm failed because of
188 * lack of convergence or because function couldn't be evaluated.
189 */
190 @SuppressWarnings("DuplicatedCode")
191 @Override
192 public void minimize() throws LockedException, NotReadyException, OptimizationException {
193
194 if (isLocked()) {
195 throw new LockedException();
196 }
197 if (!isReady()) {
198 throw new NotReadyException();
199 }
200
201 locked = true;
202
203 final var n = p.length;
204
205 // set vector of directions
206 if (!isDirectionAvailable()) {
207 xi = new double[n];
208 } else {
209 if (xi.length != n) {
210 xi = new double[n];
211 }
212 }
213
214 var validResult = false;
215 try {
216 double gg;
217 double dgg;
218
219 final var g = new double[n];
220 final var h = new double[n];
221
222 var fp = listener.evaluate(p);
223 gradientListener.evaluateGradient(p, xi);
224
225 for (var j = 0; j < n; j++) {
226 g[j] = -xi[j];
227 h[j] = g[j];
228 xi[j] = h[j];
229 }
230 for (var its = 0; its < ITMAX; its++) {
231 iter = its;
232 fret = linmin();
233 if (2.0 * Math.abs(fret - fp) <= tolerance * (Math.abs(fret) + Math.abs(fp) + EPS)) {
234 // minimum found
235 validResult = true;
236
237 if (iterationCompletedListener != null) {
238 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
239 }
240 break;
241 }
242
243 fp = fret;
244
245 gradientListener.evaluateGradient(p, xi);
246
247 var test = 0.0;
248 final var den = Math.max(Math.abs(fp), 1.0);
249 for (var j = 0; j < n; j++) {
250 final var temp = Math.abs(xi[j]) * Math.max(Math.abs(p[j]), 1.0) / den;
251
252 if (temp > test) {
253 test = temp;
254 }
255 }
256 if (test < GTOL) {
257 // minimum found
258 validResult = true;
259
260 if (iterationCompletedListener != null) {
261 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
262 }
263 break;
264 }
265
266 dgg = gg = 0.0;
267 for (var j = 0; j < n; j++) {
268 gg += g[j] * g[j];
269
270 if (isPolakRibiereEnabled()) {
271 // This statement for Polak-Ribiere
272 dgg += (xi[j] + g[j]) * xi[j];
273 } else {
274 // This statement for Fletcher-Reeves
275 dgg += xi[j] * xi[j];
276 }
277 }
278
279 if (gg == 0.0) {
280 // minimum found
281 validResult = true;
282
283 if (iterationCompletedListener != null) {
284 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
285 }
286 break;
287 }
288
289 final var gam = dgg / gg;
290 for (var j = 0; j < n; j++) {
291 g[j] = -xi[j];
292 h[j] = g[j] + gam * h[j];
293 xi[j] = h[j];
294 }
295
296 if (iterationCompletedListener != null) {
297 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
298 }
299 }
300
301 if (!validResult) {
302 // too many iterations
303 locked = false;
304 throw new OptimizationException();
305 }
306 } catch (final EvaluationException e) {
307 throw new OptimizationException(e);
308 } finally {
309 locked = false;
310 }
311
312 // set result
313 xmin = p;
314 resultAvailable = true;
315 fmin = fret;
316 }
317
318 /**
319 * Returns boolean indicating whether this instance is ready to start the
320 * estimation of a local minimum.
321 * An instance is ready once a listener, a gradient listener and a start
322 * point are provided.
323 *
324 * @return True if this instance is ready, false otherwise.
325 */
326 @Override
327 public boolean isReady() {
328 return isListenerAvailable() && isGradientListenerAvailable() && isStartPointAvailable();
329 }
330
331 /**
332 * Returns tolerance or accuracy to be expected on estimated local minimum.
333 *
334 * @return Tolerance or accuracy to be expected on estimated local minimum.
335 */
336 public double getTolerance() {
337 return tolerance;
338 }
339
340 /**
341 * Sets tolerance or accuracy to be expected on estimated local minimum.
342 *
343 * @param tolerance Tolerance or accuracy to be expected on estimated local
344 * minimum.
345 * @throws LockedException Raised if this instance is locked.
346 * @throws IllegalArgumentException Raised if provided tolerance is
347 * negative.
348 */
349 public void setTolerance(final double tolerance) throws LockedException {
350 if (isLocked()) {
351 throw new LockedException();
352 }
353 internalSetTolerance(tolerance);
354 }
355
356 /**
357 * Returns gradient listener in charge of obtaining gradient values for the
358 * function to be evaluated.
359 *
360 * @return Gradient listener.
361 * @throws NotAvailableException Raised if gradient listener has not yet
362 * been provided.
363 */
364 public GradientFunctionEvaluatorListener getGradientListener() throws NotAvailableException {
365 if (!isGradientListenerAvailable()) {
366 throw new NotAvailableException();
367 }
368 return gradientListener;
369 }
370
371 /**
372 * Sets gradient listener in charge of obtaining gradient values for the
373 * function to be evaluated.
374 *
375 * @param gradientListener Gradient listener.
376 * @throws LockedException Raised if this instance is locked.
377 */
378 public void setGradientListener(final GradientFunctionEvaluatorListener gradientListener) throws LockedException {
379 if (isLocked()) {
380 throw new LockedException();
381 }
382
383 this.gradientListener = gradientListener;
384 }
385
386 /**
387 * Returns boolean indicating whether a gradient listener has already been
388 * provided and is available for retrieval.
389 *
390 * @return True if available, false otherwise.
391 */
392 public boolean isGradientListenerAvailable() {
393 return gradientListener != null;
394 }
395
396 /**
397 * Returns boolean indicating whether Polak-Ribiere method is used or
398 * Fletcher-Reeves is used instead.
399 *
400 * @return If true, Polak-Ribiere method is used, otherwise Fletcher-Reeves
401 * is used.
402 */
403 public boolean isPolakRibiereEnabled() {
404 return usePolakRibiere;
405 }
406
407 /**
408 * Sets boolean indicating whether Polak-Ribiere method or Fletcher-Reeves
409 * method is used.
410 * If provided value is true, Polak-Ribiere method will be used, otherwise
411 * Flecther-Reeves method will be used.
412 *
413 * @param useIt Boolean to determine method.
414 * @throws LockedException Raised if this instance is locked.
415 */
416 public void setUsePolakRibiere(final boolean useIt) throws LockedException {
417 if (isLocked()) {
418 throw new LockedException();
419 }
420 usePolakRibiere = useIt;
421 }
422
423 /**
424 * Sets start point where local minimum is searched nearby.
425 *
426 * @param point Start point to search for a local minimum.
427 * @throws LockedException Raised if this instance is locked.
428 */
429 public void setStartPoint(final double[] point) throws LockedException {
430 if (isLocked()) {
431 throw new LockedException();
432 }
433 internalSetStartPoint(point);
434 }
435
436 /**
437 * Return number of iterations that were needed to estimate a minimum.
438 *
439 * @return number of iterations that were needed.
440 */
441 public int getIterations() {
442 return iter;
443 }
444
445 /**
446 * Internal method to set tolerance or accuracy to be expected on estimated
447 * local minimum.
448 * This method does not check whether this instance is locked.
449 *
450 * @param tolerance Tolerance or accuracy to be expected on estimated local
451 * minimum.
452 * @throws IllegalArgumentException Raised if provided tolerance is
453 * negative.
454 */
455 private void internalSetTolerance(final double tolerance) {
456 if (tolerance < MIN_TOLERANCE) {
457 throw new IllegalArgumentException();
458 }
459 this.tolerance = tolerance;
460 }
461
462 /**
463 * Internal method to set start point where local minimum is searched
464 * nearby.
465 * This method does not check whether this instance is locked.
466 *
467 * @param point Start point to search for a local minimum.
468 */
469 private void internalSetStartPoint(final double[] point) {
470 p = point;
471 }
472 }