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.algebra.AlgebraException;
19 import com.irurueta.algebra.Matrix;
20 import com.irurueta.numerical.EvaluationException;
21 import com.irurueta.numerical.GradientFunctionEvaluatorListener;
22 import com.irurueta.numerical.LockedException;
23 import com.irurueta.numerical.MultiDimensionFunctionEvaluatorListener;
24 import com.irurueta.numerical.NotAvailableException;
25 import com.irurueta.numerical.NotReadyException;
26
27 /**
28 * This class searches for a multi dimension function local minimum.
29 * The local minimum is searched by starting the algorithm at a start point
30 * The implementation of this class is based on Numerical Recipes 3rd ed.
31 * Section 10.9 page 521.
32 */
33 public class QuasiNewtonMultiOptimizer extends MultiOptimizer {
34 /**
35 * Maximum number of iterations.
36 */
37 public static final int ITMAX = 200;
38
39 /**
40 * Machine precision.
41 */
42 public static final double EPS = 1e-12;
43
44 /**
45 * Convergence criterion on x values.
46 */
47 public static final double TOLX = 4.0 * EPS;
48
49 /**
50 * Scaled maximum step length allowed in line searches.
51 */
52 public static final double STPMX = 100.0;
53 public static final double ALF = 1e-4;
54 public static final double TOLX2 = 1e-12;
55
56 /**
57 * Constant defining default tolerance or accuracy to be achieved on the
58 * minimum being estimated by this class.
59 */
60 public static final double DEFAULT_TOLERANCE = 3e-8;
61
62 /**
63 * Minimum allowed tolerance value.
64 */
65 public static final double MIN_TOLERANCE = 0.0;
66
67 /**
68 * n-dimensional start point.
69 */
70 private double[] p;
71
72 /**
73 * The fractional tolerance in the function value such that failure to
74 * decrease by more than this amount on one iteration signals done-ness.
75 */
76 private double tolerance;
77
78 /**
79 * Member contains number of iterations that were needed to estimate a
80 * minimum.
81 */
82 private int iter;
83
84 /**
85 * Listener to obtain gradient values for the multi dimension function being
86 * evaluated.
87 * If the gradient is unknown (e.g. doesn't have a closed expression), the
88 * provided listener could use a GradientEstimator to obtain one.
89 */
90 GradientFunctionEvaluatorListener gradientListener;
91
92 /**
93 * value of the function at the minimum.
94 */
95 double fret;
96
97 /**
98 * Empty constructor.
99 */
100 public QuasiNewtonMultiOptimizer() {
101 super();
102 tolerance = DEFAULT_TOLERANCE;
103 p = null;
104 iter = 0;
105 gradientListener = null;
106 }
107
108 /**
109 * Constructor.
110 *
111 * @param listener Listener to evaluate a multidimensional function.
112 * @param gradientListener Listener to obtain gradient value for the
113 * multidimensional function being evaluated.
114 * @param tolerance Tolerance or accuracy to be expected on estimated local
115 * minimum.
116 * @throws IllegalArgumentException Raised if tolerance is negative.
117 */
118 public QuasiNewtonMultiOptimizer(
119 final MultiDimensionFunctionEvaluatorListener listener,
120 final GradientFunctionEvaluatorListener gradientListener, final double tolerance) {
121 super(listener);
122 internalSetTolerance(tolerance);
123 p = null;
124 iter = 0;
125 this.gradientListener = gradientListener;
126 }
127
128 /**
129 * Constructor.
130 *
131 * @param listener Listener to evaluate a multidimensional function.
132 * @param gradientListener Listener to obtain gradient value for the
133 * multidimensional function being evaluated.
134 * @param startPoint Start point where algorithm will be started. Start point
135 * should be close to the local minimum to be found. Provided array must
136 * have a length equal to the number of dimensions of the function being
137 * evaluated, otherwise and exception will be raised when searching for the
138 * minimum.
139 * @param tolerance Tolerance or accuracy to be expected on estimated local
140 * minimum.
141 * @throws IllegalArgumentException Raised if tolerance is negative.
142 */
143 public QuasiNewtonMultiOptimizer(
144 final MultiDimensionFunctionEvaluatorListener listener,
145 final GradientFunctionEvaluatorListener gradientListener, final double[] startPoint,
146 final double tolerance) {
147 super(listener);
148 internalSetTolerance(tolerance);
149 p = startPoint;
150 iter = 0;
151 this.gradientListener = gradientListener;
152 }
153
154 /**
155 * This function estimates a function minimum.
156 * Implementations of this class will usually search a local minimum close
157 * to a start point.
158 * Given a starting point p[0..n-1], the Broyden-Fletcher-Goldfarb-Sharno
159 * variant of Davidon Fletcher-Powell minimization is performed on a
160 * function whose value and gradient are provided by respective listeners.
161 * The convergence requirement on zeroing the gradient is provided by
162 * tolerance. This method estimates the location of the minimum, sets the
163 * number of iterations that were required and the minimum value of the
164 * function at the minimum.
165 * The internal method lnsrch is called within this method to perform
166 * approximate line minimizations.
167 *
168 * @throws LockedException Raised if this instance is locked, because
169 * estimation is being computed.
170 * @throws NotReadyException Raised if this instance is not ready, because
171 * a listener, a gradient listener and a start point haven't been provided
172 * @throws OptimizationException Raised if the algorithm failed because of
173 * lack of convergence or because function couldn't be evaluated.
174 */
175 @Override
176 public void minimize() throws LockedException, NotReadyException, OptimizationException {
177
178 if (isLocked()) {
179 throw new LockedException();
180 }
181 if (!isReady()) {
182 throw new NotReadyException();
183 }
184
185 locked = true;
186
187 var validResult = false;
188 final var n = p.length;
189
190 try {
191 double den;
192 double fac;
193 double fad;
194 double fae;
195 double fp;
196 final double stpmax;
197 var sum = 0.0;
198 double sumdg;
199 double sumxi;
200 double temp;
201 double test;
202 final var dg = new double[n];
203 final var g = new double[n];
204 final var hdg = new double[n];
205 final var pnew = new double[n];
206 final var xi = new double[n];
207 final var fretArray = new double[1];
208 final var check = new boolean[1];
209
210 final var hessin = new Matrix(n, n);
211 fp = listener.evaluate(p);
212 gradientListener.evaluateGradient(p, g);
213
214 for (var i = 0; i < n; i++) {
215 for (var j = 0; j < n; j++) {
216 hessin.setElementAt(i, j, 0.0);
217 }
218 hessin.setElementAt(i, i, 1.0);
219 xi[i] = -g[i];
220 sum += p[i] * p[i];
221 }
222 stpmax = STPMX * Math.max(Math.sqrt(sum), n);
223 for (var its = 0; its < ITMAX; its++) {
224 iter = its;
225 lnsrch(p, fp, g, xi, pnew, fretArray, stpmax, check);
226 fret = fretArray[0];
227 fp = fret;
228 for (var i = 0; i < n; i++) {
229 xi[i] = pnew[i] - p[i];
230 p[i] = pnew[i];
231 }
232 test = 0.0;
233 for (var i = 0; i < n; i++) {
234 temp = Math.abs(xi[i]) / Math.max(Math.abs(p[i]), 1.0);
235 if (temp > test) {
236 test = temp;
237 }
238 }
239 if (test < TOLX) {
240 // minimum found
241 validResult = true;
242 break;
243 }
244
245 System.arraycopy(g, 0, dg, 0, n);
246
247 gradientListener.evaluateGradient(p, g);
248
249 test = 0.0;
250 den = Math.max(Math.abs(fret), 1.0);
251
252 for (var i = 0; i < n; i++) {
253 temp = Math.abs(g[i]) * Math.max(Math.abs(p[i]), 1.0) / den;
254 if (temp > test) {
255 test = temp;
256 }
257 }
258 if (test < tolerance) {
259 // minimum found
260 validResult = true;
261 break;
262 }
263
264 for (var i = 0; i < n; i++) {
265 dg[i] = g[i] - dg[i];
266 }
267
268 for (var i = 0; i < n; i++) {
269 hdg[i] = 0.0;
270 for (var j = 0; j < n; j++) {
271 hdg[i] += hessin.getElementAt(i, j) * dg[j];
272 }
273 }
274 fac = fae = sumdg = sumxi = 0.0;
275 for (var i = 0; i < n; i++) {
276 fac += dg[i] * xi[i];
277 fae += dg[i] * hdg[i];
278 sumdg += sqr(dg[i]);
279 sumxi += sqr(xi[i]);
280 }
281
282 if (fac > Math.sqrt(EPS * sumdg * sumxi)) {
283 fac = 1.0 / fac;
284 fad = 1.0 / fae;
285 for (var i = 0; i < n; i++) {
286 dg[i] = fac * xi[i] - fad * hdg[i];
287 }
288
289 for (var i = 0; i < n; i++) {
290 for (var j = i; j < n; j++) {
291 hessin.setElementAt(i, j, hessin.getElementAt(i, j) + fac * xi[i] * xi[j]
292 - fad * hdg[i] * hdg[j] + fae * dg[i] * dg[j]);
293 hessin.setElementAt(j, i, hessin.getElementAt(i, j));
294 }
295 }
296 }
297 for (var i = 0; i < n; i++) {
298 xi[i] = 0.0;
299 for (var j = 0; j < n; j++) {
300 xi[i] -= hessin.getElementAt(i, j) * g[j];
301 }
302 }
303
304 if (iterationCompletedListener != null) {
305 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
306 }
307 }
308
309 if (!validResult) {
310 // too many iterations
311 locked = false;
312 throw new OptimizationException();
313 }
314 } catch (final AlgebraException | EvaluationException e) {
315 throw new OptimizationException(e);
316 } finally {
317 locked = false;
318 }
319
320 // set result
321 xmin = p;
322 resultAvailable = true;
323 fmin = fret;
324 }
325
326 /**
327 * Returns boolean indicating whether this instance is ready to start the
328 * estimation of a local minimum.
329 * An instance is ready once a listener, a gradient listener and a start
330 * point are provided.
331 *
332 * @return True if this instance is ready, false otherwise.
333 */
334 @Override
335 public boolean isReady() {
336 return isListenerAvailable() && isGradientListenerAvailable() && isStartPointAvailable();
337 }
338
339 /**
340 * Returns tolerance or accuracy to be expected on estimated local minimum.
341 *
342 * @return Tolerance or accuracy to be expected on estimated local minimum.
343 */
344 public double getTolerance() {
345 return tolerance;
346 }
347
348 /**
349 * Internal method to set tolerance or accuracy to be expected on estimated
350 * local minimum.
351 * This method does not check whether this instance is locked.
352 *
353 * @param tolerance Tolerance or accuracy to be expected on estimated local
354 * minimum.
355 * @throws IllegalArgumentException Raised if provided tolerance is
356 * negative.
357 */
358 private void internalSetTolerance(final double tolerance) {
359 if (tolerance < MIN_TOLERANCE) {
360 throw new IllegalArgumentException();
361 }
362 this.tolerance = tolerance;
363 }
364
365 /**
366 * Sets tolerance or accuracy to be expected on estimated local minimum.
367 *
368 * @param tolerance Tolerance or accuracy to be expected on estimated local
369 * minimum.
370 * @throws LockedException Raised if this instance is locked.
371 * @throws IllegalArgumentException Raised if provided tolerance is
372 * negative.
373 */
374 public void setTolerance(final double tolerance) throws LockedException {
375 if (isLocked()) {
376 throw new LockedException();
377 }
378 internalSetTolerance(tolerance);
379 }
380
381 /**
382 * Returns boolean indicating whether start point has already been provided
383 * and is ready for retrieval.
384 *
385 * @return True if available, false otherwise.
386 */
387 public boolean isStartPointAvailable() {
388 return p != null;
389 }
390
391 /**
392 * Sets start point where algorithm will be started. Start point should be
393 * close to the local minimum to be found.
394 *
395 * @param startPoint Start point where algorithm will be started.
396 * @throws LockedException Raised if this instance is locked.
397 */
398 public void setStartPoint(final double[] startPoint) throws LockedException {
399 if (isLocked()) {
400 throw new LockedException();
401 }
402 p = startPoint;
403 }
404
405 /**
406 * Returns start point where algorithm will be started. Start point should
407 * be close to the local minimum to be found.
408 *
409 * @return Start point where algorithm will be started.
410 * @throws NotAvailableException Raised if start point has not yet been
411 * provided and is not available.
412 */
413 public double[] getStartPoint() throws NotAvailableException {
414 if (!isStartPointAvailable()) {
415 throw new NotAvailableException();
416 }
417 return p;
418 }
419
420 /**
421 * Returns gradient listener in charge of obtaining gradient values for the
422 * function to be evaluated.
423 *
424 * @return Gradient listener.
425 * @throws NotAvailableException Raised if gradient listener has not yet
426 * been provided.
427 */
428 public GradientFunctionEvaluatorListener getGradientListener() throws NotAvailableException {
429 if (!isGradientListenerAvailable()) {
430 throw new NotAvailableException();
431 }
432 return gradientListener;
433 }
434
435 /**
436 * Sets gradient listener in charge of obtaining gradient values for the
437 * function to be evaluated.
438 *
439 * @param gradientListener Gradient listener.
440 * @throws LockedException Raised if this instance is locked.
441 */
442 public void setGradientListener(final GradientFunctionEvaluatorListener gradientListener) throws LockedException {
443 if (isLocked()) {
444 throw new LockedException();
445 }
446 this.gradientListener = gradientListener;
447 }
448
449 /**
450 * Returns boolean indicating whether a gradient listener has already been
451 * provided and is available for retrieval.
452 *
453 * @return True if available, false otherwise.
454 */
455 public boolean isGradientListenerAvailable() {
456 return gradientListener != null;
457 }
458
459 /**
460 * Return number of iterations that were needed to estimate a minimum.
461 *
462 * @return number of iterations that were needed.
463 */
464 public int getIterations() {
465 return iter;
466 }
467
468 /**
469 * Computes the squared value of provided double.
470 *
471 * @param x Value to be squared.
472 * @return Squared value.
473 */
474 private double sqr(final double x) {
475 return x * x;
476 }
477
478 /**
479 * Internal method to search for a minimum along a line.
480 * NOTE: comments of params might be incorrect.
481 *
482 * @param xold previous point.
483 * @param fold previous function evaluation.
484 * @param g gradient.
485 * @param p point.
486 * @param x point.
487 * @param f function.
488 * @param stpmax stores point maximum.
489 * @param check checks line search.
490 * @throws OptimizationException Raised if convergence was not achieved.
491 */
492 private void lnsrch(final double[] xold, final double fold, final double[] g, final double[] p, final double[] x,
493 final double[] f, final double stpmax, final boolean[] check) throws OptimizationException {
494 double a;
495 double alam;
496 var alam2 = 0.0;
497 final double alamin;
498 double b;
499 double disc;
500 var f2 = 0.0;
501 double rhs1;
502 double rhs2;
503 var slope = 0.0;
504 var sum = 0.0;
505 double temp;
506 double test;
507 double tmplam;
508 int i;
509 final var n = xold.length;
510 check[0] = false;
511
512 for (i = 0; i < n; i++) {
513 sum += p[i] * p[i];
514 }
515 sum = Math.sqrt(sum);
516
517 if (sum > stpmax) {
518 for (i = 0; i < n; i++) {
519 p[i] *= stpmax / sum;
520 }
521 }
522 for (i = 0; i < n; i++) {
523 slope += g[i] * p[i];
524 }
525 if (slope >= 0.0) {
526 throw new OptimizationException();
527 }
528 test = 0.0;
529 for (i = 0; i < n; i++) {
530 temp = Math.abs(p[i]) / Math.max(Math.abs(xold[i]), 1.0);
531 if (temp > test) {
532 test = temp;
533 }
534 }
535 alamin = TOLX2 / test;
536 alam = 1.0;
537 for (; ; ) {
538 for (i = 0; i < n; i++) {
539 x[i] = xold[i] + alam * p[i];
540 }
541 try {
542 f[0] = listener.evaluate(x);
543 } catch (final EvaluationException e) {
544 throw new OptimizationException(e);
545 }
546
547 if (alam < alamin) {
548 System.arraycopy(xold, 0, x, 0, n);
549 check[0] = true;
550 return;
551 } else if (f[0] <= fold + ALF * alam * slope) {
552 return;
553 } else {
554 if (alam == 1.0) {
555 tmplam = -slope / (2.0 * (f[0] - fold - slope));
556 } else {
557 rhs1 = f[0] - fold - alam * slope;
558 rhs2 = f2 - fold - alam2 * slope;
559 a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2);
560 b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2);
561 if (a == 0.0) {
562 tmplam = -slope / (2.0 * b);
563 } else {
564 disc = b * b - 3.0 * a * slope;
565 if (disc < 0.0) {
566 tmplam = 0.5 * alam;
567 } else if (b <= 0.0) {
568 tmplam = (-b + Math.sqrt(disc)) / (3.0 * a);
569 } else {
570 tmplam = -slope / (b + Math.sqrt(disc));
571 }
572 }
573 if (tmplam > 0.5 * alam) {
574 tmplam = 0.5 * alam;
575 }
576 }
577 }
578 alam2 = alam;
579 f2 = f[0];
580 alam = Math.max(tmplam, 0.1 * alam);
581 }
582 }
583 }