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.NotReadyException;
23
24 /**
25 * This class searches for a multi dimension function local minimum.
26 * The local minimum is searched by starting the algorithm at a start point
27 * and a given direction which should be close and point to the local minimum to
28 * be found to achieve the best accuracy with the lowest number of iterations.
29 * NOTE: this algorithm might not have proper convergence in some situations,
30 * but it is ensured to provide faster convergence than other algorithms such
31 * as Brent because gradient information is also provided.
32 * The implementation of this class is based on Numerical Recipes 3rd ed.
33 * Section 10.8 page 518.
34 */
35 public class DerivativeConjugateGradientMultiOptimizer extends DerivativeLineMultiOptimizer {
36
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 * Boolean indicating whether Polak-Ribiere method is used if true,
88 * otherwise Fletcher-Reeves will be used.
89 */
90 private boolean usePolakRibiere;
91
92 /**
93 * Empty constructor.
94 */
95 public DerivativeConjugateGradientMultiOptimizer() {
96 super();
97 tolerance = DEFAULT_TOLERANCE;
98 usePolakRibiere = DEFAULT_USE_POLAK_RIBIERE;
99 }
100
101 /**
102 * Constructor.
103 *
104 * @param listener Listener to evaluate a multi-dimension function.
105 * @param gradientListener Listener to obtain gradient value for the
106 * multi-dimension function being evaluated.
107 * @param point Start point where algorithm will be started. Start point
108 * should be close to the local minimum to be found. Provided array must
109 * have a length equal to the number of dimensions of the function being
110 * evaluated, otherwise and exception will be raised when searching for the
111 * minimum.
112 * @param direction Direction to start looking for a minimum. Provided array
113 * must have the same length as the number of dimensions of the function
114 * being evaluated. Provided direction is considered as a vector pointing
115 * to the minimum to be found.
116 * @param tolerance Tolerance or accuracy to be expected on estimated local
117 * minimum.
118 * @param usePolakRibiere True if Polak-Ribiere method is used, otherwise
119 * Fletcher-Reeves will be used.
120 * @throws IllegalArgumentException Raised if provided point and direction
121 * don't have the same length or if provided tolerance is negative.
122 */
123 public DerivativeConjugateGradientMultiOptimizer(
124 final MultiDimensionFunctionEvaluatorListener listener,
125 final GradientFunctionEvaluatorListener gradientListener, final double[] point, final double[] direction,
126 final double tolerance, final boolean usePolakRibiere) {
127 super(listener, gradientListener, point, direction);
128 internalSetTolerance(tolerance);
129 this.usePolakRibiere = usePolakRibiere;
130 }
131
132 /**
133 * This function estimates a function minimum.
134 * Implementations of this class will usually search a local minimum close
135 * to a start point and will start looking into provided start direction.
136 * Minimization of a function f. Input consists of an initial starting point
137 * p. The initial matrix xi, whose columns contain the initial set of
138 * directions, is set to the identity. Returned is the best point found, at
139 * which point fret is the minimum function value and iter is the number of
140 * iterations taken.
141 * Note: the starting point p will be modified after execution of this
142 * method.
143 *
144 * @throws LockedException Raised if this instance is locked, because
145 * estimation is being computed.
146 * @throws NotReadyException Raised if this instance is not ready, because
147 * a listener, a gradient listener and a start point haven't been provided.
148 * @throws OptimizationException Raised if the algorithm failed because of
149 * lack of convergence or because function couldn't be evaluated.
150 */
151 @Override
152 @SuppressWarnings("Duplicates")
153 public void minimize() throws LockedException, NotReadyException, OptimizationException {
154 if (isLocked()) {
155 throw new LockedException();
156 }
157 if (!isReady()) {
158 throw new NotReadyException();
159 }
160
161 locked = true;
162
163 final var n = p.length;
164
165 // set vector of directions
166 if (!isDirectionAvailable()) {
167 xi = new double[n];
168 } else {
169 if (xi.length != n) {
170 xi = new double[n];
171 }
172 }
173
174 var validResult = false;
175 try {
176 double gg;
177 double dgg;
178
179 final var g = new double[n];
180 final var h = new double[n];
181
182 var fp = listener.evaluate(p);
183 gradientListener.evaluateGradient(p, xi);
184
185 for (var j = 0; j < n; j++) {
186 g[j] = -xi[j];
187 h[j] = g[j];
188 xi[j] = h[j];
189 }
190 for (var its = 0; its < ITMAX; its++) {
191 iter = its;
192 fret = linmin();
193 if (2.0 * Math.abs(fret - fp) <= tolerance * (Math.abs(fret) + Math.abs(fp) + EPS)) {
194 // minimum found
195 validResult = true;
196
197 if (iterationCompletedListener != null) {
198 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
199 }
200 break;
201 }
202
203 fp = fret;
204
205 gradientListener.evaluateGradient(p, xi);
206
207 var test = 0.0;
208 final var den = Math.max(Math.abs(fp), 1.0);
209 for (var j = 0; j < n; j++) {
210 final var temp = Math.abs(xi[j]) * Math.max(Math.abs(p[j]), 1.0) / den;
211
212 if (temp > test) {
213 test = temp;
214 }
215 }
216 if (test < GTOL) {
217 // minimum found
218 validResult = true;
219
220 if (iterationCompletedListener != null) {
221 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
222 }
223 break;
224 }
225
226 dgg = gg = 0.0;
227 for (var j = 0; j < n; j++) {
228 gg += g[j] * g[j];
229
230 if (isPolakRibiereEnabled()) {
231 // This statement for Polak-Ribiere
232 dgg += xi[j] + g[j] * xi[j];
233 } else {
234 // This statement for Fletcher-Reeves
235 dgg += xi[j] * xi[j];
236 }
237 }
238
239 if (gg == 0.0) {
240 // minimum found
241 validResult = true;
242
243 if (iterationCompletedListener != null) {
244 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
245 }
246 break;
247 }
248
249 final var gam = dgg / gg;
250 for (var j = 0; j < n; j++) {
251 g[j] = -xi[j];
252 h[j] = g[j] + gam * h[j];
253 xi[j] = h[j];
254 }
255
256 if (iterationCompletedListener != null) {
257 iterationCompletedListener.onIterationCompleted(this, its, ITMAX);
258 }
259 }
260
261 if (!validResult) {
262 // too many iterations
263 locked = false;
264 throw new OptimizationException();
265 }
266 } catch (final EvaluationException e) {
267 throw new OptimizationException(e);
268 } finally {
269 locked = false;
270 }
271
272 // set result
273 xmin = p;
274 resultAvailable = true;
275 fmin = fret;
276 }
277
278 /**
279 * Returns boolean indicating whether this instance is ready to start the
280 * estimation of a local minimum.
281 * An instance is ready once a listener, a gradient listener and a start
282 * point are provided.
283 *
284 * @return True if this instance is ready, false otherwise.
285 */
286 @Override
287 public boolean isReady() {
288 return isListenerAvailable() && isGradientListenerAvailable() && isStartPointAvailable();
289 }
290
291 /**
292 * Returns tolerance or accuracy to be expected on estimated local minimum.
293 *
294 * @return Tolerance or accuracy to be expected on estimated local minimum.
295 */
296 public double getTolerance() {
297 return tolerance;
298 }
299
300 /**
301 * Sets tolerance or accuracy to be expected on estimated local minimum.
302 *
303 * @param tolerance Tolerance or accuracy to be expected on estimated local
304 * minimum.
305 * @throws LockedException Raised if this instance is locked.
306 * @throws IllegalArgumentException Raised if provided tolerance is
307 * negative.
308 */
309 public void setTolerance(final double tolerance) throws LockedException {
310 if (isLocked()) {
311 throw new LockedException();
312 }
313 internalSetTolerance(tolerance);
314 }
315
316 /**
317 * Returns boolean indicating whether Polak-Ribiere method is used or
318 * Fletcher-Reeves is used instead.
319 *
320 * @return If true, Polak-Ribiere method is used, otherwise Fletcher-Reeves
321 * is used.
322 */
323 public boolean isPolakRibiereEnabled() {
324 return usePolakRibiere;
325 }
326
327 /**
328 * Sets boolean indicating whether Polak-Ribiere method or Fletcher-Reeves
329 * method is used.
330 * If provided value is true, Polak-Ribiere method will be used, otherwise
331 * Flecther-Reeves method will be used.
332 *
333 * @param useIt Boolean to determine method.
334 * @throws LockedException Raised if this instance is locked.
335 */
336 public void setUsePolakRibiere(final boolean useIt) throws LockedException {
337 if (isLocked()) {
338 throw new LockedException();
339 }
340 this.usePolakRibiere = useIt;
341 }
342
343 /**
344 * Sets start point where local minimum is searched nearby.
345 *
346 * @param point Start point to search for a local minimum.
347 * @throws LockedException Raised if this instance is locked.
348 */
349 public void setStartPoint(final double[] point) throws LockedException {
350 if (isLocked()) {
351 throw new LockedException();
352 }
353 internalSetStartPoint(point);
354 }
355
356 /**
357 * Return number of iterations that were needed to estimate a minimum.
358 *
359 * @return number of iterations that were needed.
360 */
361 public int getIterations() {
362 return iter;
363 }
364
365 /**
366 * Internal method to set tolerance or accuracy to be expected on estimated
367 * local minimum.
368 * This method does not check whether this instance is locked.
369 *
370 * @param tolerance Tolerance or accuracy to be expected on estimated local
371 * minimum.
372 * @throws IllegalArgumentException Raised if provided tolerance is
373 * negative.
374 */
375 private void internalSetTolerance(final double tolerance) {
376 if (tolerance < MIN_TOLERANCE) {
377 throw new IllegalArgumentException();
378 }
379 this.tolerance = tolerance;
380 }
381
382 /**
383 * Internal method to set start point where local minimum is searched
384 * nearby.
385 * This method does not check whether this instance is locked.
386 *
387 * @param point Start point to search for a local minimum.
388 */
389 private void internalSetStartPoint(final double[] point) {
390 p = point;
391 }
392 }