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


Generated on 04-09-2015 for elastix by doxygen 1.8.9.1 elastix logo