go home Home | Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Data Structures | File List | Namespace Members | Data Fields | Globals | Related Pages
itkQuasiNewtonLBFGSOptimizer.h
Go to the documentation of this file.
1/*=========================================================================
2 *
3 * Copyright UMC Utrecht and contributors
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
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
19#ifndef itkQuasiNewtonLBFGSOptimizer_h
20#define itkQuasiNewtonLBFGSOptimizer_h
21
24#include <vector>
25
26namespace itk
27{
59{
60public:
62
65 using Pointer = SmartPointer<Self>;
66 using ConstPointer = SmartPointer<const Self>;
67
68 itkNewMacro(Self);
70
71 using Superclass::ParametersType;
72 using Superclass::DerivativeType;
73 using Superclass::CostFunctionType;
75 using Superclass::MeasureType;
77
78 using RhoType = Array<double>;
79 using SType = std::vector<ParametersType>;
80 using YType = std::vector<DerivativeType>;
81 using DiagonalMatrixType = Array<double>;
83
85
87 {
95 };
96
97 void
99
100 virtual void
102
103 virtual void
105
107 itkGetConstMacro(CurrentIteration, unsigned long);
108 itkGetConstMacro(CurrentValue, MeasureType);
109 itkGetConstReferenceMacro(CurrentGradient, DerivativeType);
110 itkGetConstMacro(InLineSearch, bool);
111 itkGetConstReferenceMacro(StopCondition, StopConditionType);
112 itkGetConstMacro(CurrentStepLength, double);
113
117
119 itkGetConstMacro(MaximumNumberOfIterations, unsigned long);
120 itkSetClampMacro(MaximumNumberOfIterations, unsigned long, 1, NumericTraits<unsigned long>::max());
121
127 itkGetConstMacro(GradientMagnitudeTolerance, double);
128 itkSetMacro(GradientMagnitudeTolerance, double);
129
133 itkSetMacro(Memory, unsigned int);
134 itkGetConstMacro(Memory, unsigned int);
135
136protected:
138 ~QuasiNewtonLBFGSOptimizer() override = default;
139
140 // \todo: should be implemented
141 void
142 PrintSelf(std::ostream & os, Indent indent) const override
143 {}
144
145 DerivativeType m_CurrentGradient;
146 MeasureType m_CurrentValue{ 0.0 };
147 unsigned long m_CurrentIteration{ 0 };
149 bool m_Stop{ false };
150 double m_CurrentStepLength{ 0.0 };
151
153 bool m_InLineSearch{ false };
154
158
159 unsigned int m_Point{ 0 };
160 unsigned int m_PreviousPoint{ 0 };
161 unsigned int m_Bound{ 0 };
162
163 itkSetMacro(InLineSearch, bool);
164
169 virtual void
171
178 virtual void
179 ComputeSearchDirection(const DerivativeType & gradient, ParametersType & searchDir);
180
185 virtual void
186 LineSearch(const ParametersType searchDir, double & step, ParametersType & x, MeasureType & f, DerivativeType & g);
187
190 virtual void
191 StoreCurrentPoint(const ParametersType & step, const DerivativeType & grad_dif);
192
197 virtual bool
198 TestConvergence(bool firstLineSearchDone);
199
200private:
201 unsigned long m_MaximumNumberOfIterations{ 100 };
204 unsigned int m_Memory{ 5 };
205};
206
207} // end namespace itk
208
216/* LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION */
217/* JORGE NOCEDAL */
218/* *** July 1990 *** */
219
220/* This subroutine solves the unconstrained minimization problem */
221
222/* min F(x), x= (x1,x2,...,xN), */
223
224/* using the limited memory BFGS method. The routine is especially */
225/* effective on problems involving a large number of variables. In */
226/* a typical iteration of this method an approximation Hk to the */
227/* inverse of the Hessian is obtained by applying M BFGS updates to */
228/* a diagonal matrix Hk0, using information from the previous M steps. */
229/* The user specifies the number M, which determines the amount of */
230/* storage required by the routine. The user may also provide the */
231/* diagonal matrices Hk0 if not satisfied with the default choice. */
232/* The algorithm is described in "On the limited memory BFGS method */
233/* for large scale optimization", by D. Liu and J. Nocedal, */
234/* Mathematical Programming B 45 (1989) 503-528. */
235
236/* The user is required to calculate the function value F and its */
237/* gradient G. In order to allow the user complete control over */
238/* these computations, reverse communication is used. The routine */
239/* must be called repeatedly under the control of the parameter */
240/* IFLAG. */
241
242/* The steplength is determined at each iteration by means of the */
243/* line search routine MCVSRCH, which is a slight modification of */
244/* the routine CSRCH written by More' and Thuente. */
245
246/* The calling statement is */
247
248/* CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) */
249
250/* where */
251
252/* N is an INTEGER variable that must be set by the user to the */
253/* number of variables. It is not altered by the routine. */
254/* Restriction: N>0. */
255
256/* M is an INTEGER variable that must be set by the user to */
257/* the number of corrections used in the BFGS update. It */
258/* is not altered by the routine. Values of M less than 3 are */
259/* not recommended; large values of M will result in excessive */
260/* computing time. 3<= M <=7 is recommended. Restriction: M>0. */
261
262/* X is a DOUBLE PRECISION array of length N. On initial entry */
263/* it must be set by the user to the values of the initial */
264/* estimate of the solution vector. On exit with IFLAG=0, it */
265/* contains the values of the variables at the best point */
266/* found (usually a solution). */
267
268/* F is a DOUBLE PRECISION variable. Before initial entry and on */
269/* a re-entry with IFLAG=1, it must be set by the user to */
270/* contain the value of the function F at the point X. */
271
272/* G is a DOUBLE PRECISION array of length N. Before initial */
273/* entry and on a re-entry with IFLAG=1, it must be set by */
274/* the user to contain the components of the gradient G at */
275/* the point X. */
276
277/* DIAGCO is a LOGICAL variable that must be set to .TRUE. if the */
278/* user wishes to provide the diagonal matrix Hk0 at each */
279/* iteration. Otherwise it should be set to .FALSE., in which */
280/* case LBFGS will use a default value described below. If */
281/* DIAGCO is set to .TRUE. the routine will return at each */
282/* iteration of the algorithm with IFLAG=2, and the diagonal */
283/* matrix Hk0 must be provided in the array DIAG. */
284
285/* DIAG is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE., */
286/* then on initial entry or on re-entry with IFLAG=2, DIAG */
287/* it must be set by the user to contain the values of the */
288/* diagonal matrix Hk0. Restriction: all elements of DIAG */
289/* must be positive. */
290
291/* IPRINT is an INTEGER array of length two which must be set by the */
292/* user. */
293
294/* IPRINT(1) specifies the frequency of the output: */
295/* IPRINT(1) < 0 : no output is generated, */
296/* IPRINT(1) = 0 : output only at first and last iteration, */
297/* IPRINT(1) > 0 : output every IPRINT(1) iterations. */
298
299/* IPRINT(2) specifies the type of output generated: */
300/* IPRINT(2) = 0 : iteration count, number of function */
301/* evaluations, function value, norm of the */
302/* gradient, and steplength, */
303/* IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of */
304/* variables and gradient vector at the */
305/* initial point, */
306/* IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of */
307/* variables, */
308/* IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.*/
309
310/* EPS is a positive DOUBLE PRECISION variable that must be set by */
311/* the user, and determines the accuracy with which the solution*/
312/* is to be found. The subroutine terminates when */
313
314/* ||G|| < EPS max(1,||X||), */
315
316/* where ||.|| denotes the Euclidean norm. */
317
318/* XTOL is a positive DOUBLE PRECISION variable that must be set by */
319/* the user to an estimate of the machine precision (e.g. */
320/* 10**(-16) on a SUN station 3/60). The line search routine will*/
321/* terminate if the relative width of the interval of uncertainty*/
322/* is less than XTOL. */
323
324/* W is a DOUBLE PRECISION array of length N(2M+1)+2M used as */
325/* workspace for LBFGS. This array must not be altered by the */
326/* user. */
327
328/* IFLAG is an INTEGER variable that must be set to 0 on initial entry*/
329/* to the subroutine. A return with IFLAG<0 indicates an error, */
330/* and IFLAG=0 indicates that the routine has terminated without*/
331/* detecting errors. On a return with IFLAG=1, the user must */
332/* evaluate the function F and gradient G. On a return with */
333/* IFLAG=2, the user must provide the diagonal matrix Hk0. */
334
335/* The following negative values of IFLAG, detecting an error, */
336/* are possible: */
337
338/* IFLAG=-1 The line search routine MCSRCH failed. The */
339/* parameter INFO provides more detailed information */
340/* (see also the documentation of MCSRCH): */
341
342/* INFO = 0 IMPROPER INPUT PARAMETERS. */
343
344/* INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF */
345/* UNCERTAINTY IS AT MOST XTOL. */
346
347/* INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE */
348/* REQUIRED AT THE PRESENT ITERATION. */
349
350/* INFO = 4 THE STEP IS TOO SMALL. */
351
352/* INFO = 5 THE STEP IS TOO LARGE. */
353
354/* INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.*/
355/* THERE MAY NOT BE A STEP WHICH SATISFIES */
356/* THE SUFFICIENT DECREASE AND CURVATURE */
357/* CONDITIONS. TOLERANCES MAY BE TOO SMALL. */
358
359/* IFLAG=-2 The i-th diagonal element of the diagonal inverse */
360/* Hessian approximation, given in DIAG, is not */
361/* positive. */
362
363/* IFLAG=-3 Improper input parameters for LBFGS (N or M are */
364/* not positive). */
365
366/* ON THE DRIVER: */
367
368/* The program that calls LBFGS must contain the declaration: */
369
370/* EXTERNAL LB2 */
371
372/* LB2 is a BLOCK DATA that defines the default values of several */
373/* parameters described in the COMMON section. */
374
375/* COMMON: */
376
377/* The subroutine contains one common area, which the user may wish to */
378/* reference: */
379
380/* awf added stpawf */
381
382/* MP is an INTEGER variable with default value 6. It is used as the */
383/* unit number for the printing of the monitoring information */
384/* controlled by IPRINT. */
385
386/* LP is an INTEGER variable with default value 6. It is used as the */
387/* unit number for the printing of error messages. This printing */
388/* may be suppressed by setting LP to a non-positive value. */
389
390/* GTOL is a DOUBLE PRECISION variable with default value 0.9, which */
391/* controls the accuracy of the line search routine MCSRCH. If the */
392/* function and gradient evaluations are inexpensive with respect */
393/* to the cost of the iteration (which is sometimes the case when */
394/* solving very large problems) it may be advantageous to set GTOL */
395/* to a small value. A typical small value is 0.1. Restriction: */
396/* GTOL should be greater than 1.D-04. */
397
398/* STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which */
399/* specify lower and uper bounds for the step in the line search. */
400/* Their default values are 1.D-20 and 1.D+20, respectively. These */
401/* values need not be modified unless the exponents are too large */
402/* for the machine being used, or unless the problem is extremely */
403/* badly scaled (in which case the exponents should be increased). */
404
405/* MACHINE DEPENDENCIES */
406
407/* The only variables that are machine-dependent are XTOL, */
408/* STPMIN and STPMAX. */
409
410/* GENERAL INFORMATION */
411
412/* Other routines called directly: DAXPY, DDOT, LB1, MCSRCH */
413
414/* Input/Output : No input; diagnostic messages on unit MP and */
415/* error messages on unit LP. */
416
417/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
418
419#endif //#ifndef itkQuasiNewtonLBFGSOptimizer_h
A base class for LineSearch optimizers.
SmartPointer< Self > Pointer
ITK version of the lbfgs algorithm ...
itkGetModifiableObjectMacro(LineSearchOptimizer, LineSearchOptimizerType)
std::vector< DerivativeType > YType
LineSearchOptimizerPointer m_LineSearchOptimizer
ITK_DISALLOW_COPY_AND_MOVE(QuasiNewtonLBFGSOptimizer)
virtual void ComputeSearchDirection(const DerivativeType &gradient, ParametersType &searchDir)
virtual void ComputeDiagonalMatrix(DiagonalMatrixType &diag_H0)
~QuasiNewtonLBFGSOptimizer() override=default
virtual void StoreCurrentPoint(const ParametersType &step, const DerivativeType &grad_dif)
virtual bool TestConvergence(bool firstLineSearchDone)
void StartOptimization() override
void PrintSelf(std::ostream &os, Indent indent) const override
std::vector< ParametersType > SType
virtual void LineSearch(const ParametersType searchDir, double &step, ParametersType &x, MeasureType &f, DerivativeType &g)
A cost function that applies a scaling to another cost function.


Generated on 2023-01-13 for elastix by doxygen 1.9.6 elastix logo