![]() |
Home | Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Data Structures | File List | Namespace Members | Data Fields | Globals | Related Pages |
00001 /*====================================================================== 00002 00003 This file is part of the elastix software. 00004 00005 Copyright (c) University Medical Center Utrecht. All rights reserved. 00006 See src/CopyrightElastix.txt or http://elastix.isi.uu.nl/legal.php for 00007 details. 00008 00009 This software is distributed WITHOUT ANY WARRANTY; without even 00010 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 00011 PURPOSE. See the above copyright notices for more information. 00012 00013 ======================================================================*/ 00014 00015 #ifndef __itkQuasiNewtonLBFGSOptimizer_h 00016 #define __itkQuasiNewtonLBFGSOptimizer_h 00017 00018 #include "itkScaledSingleValuedNonLinearOptimizer.h" 00019 #include "itkLineSearchOptimizer.h" 00020 #include <vector> 00021 00022 namespace itk 00023 { 00054 class QuasiNewtonLBFGSOptimizer : public ScaledSingleValuedNonLinearOptimizer 00055 { 00056 public: 00057 00058 typedef QuasiNewtonLBFGSOptimizer Self; 00059 typedef ScaledSingleValuedNonLinearOptimizer Superclass; 00060 typedef SmartPointer<Self> Pointer; 00061 typedef SmartPointer<const Self> ConstPointer; 00062 00063 itkNewMacro(Self); 00064 itkTypeMacro(QuasiNewtonLBFGSOptimizer, ScaledSingleValuedNonLinearOptimizer); 00065 00066 typedef Superclass::ParametersType ParametersType; 00067 typedef Superclass::DerivativeType DerivativeType; 00068 typedef Superclass::CostFunctionType CostFunctionType; 00069 typedef Superclass::ScaledCostFunctionType ScaledCostFunctionType; 00070 typedef Superclass::MeasureType MeasureType; 00071 typedef Superclass::ScalesType ScalesType; 00072 00073 typedef Array<double> RhoType; 00074 typedef std::vector<ParametersType> SType; 00075 typedef std::vector<DerivativeType> YType; 00076 typedef Array<double> DiagonalMatrixType; 00077 typedef LineSearchOptimizer LineSearchOptimizerType; 00078 00079 typedef LineSearchOptimizerType::Pointer LineSearchOptimizerPointer; 00080 00081 typedef enum { 00082 MetricError, 00083 LineSearchError, 00084 MaximumNumberOfIterations, 00085 InvalidDiagonalMatrix, 00086 GradientMagnitudeTolerance, 00087 ZeroStep, 00088 Unknown } StopConditionType; 00089 00090 virtual void StartOptimization(void); 00091 virtual void ResumeOptimization(void); 00092 virtual void StopOptimization(void); 00093 00095 itkGetConstMacro(CurrentIteration, unsigned long); 00096 itkGetConstMacro(CurrentValue, MeasureType); 00097 itkGetConstReferenceMacro(CurrentGradient, DerivativeType); 00098 itkGetConstMacro(InLineSearch, bool); 00099 itkGetConstReferenceMacro(StopCondition, StopConditionType); 00100 itkGetConstMacro(CurrentStepLength, double); 00101 00103 itkSetObjectMacro(LineSearchOptimizer, LineSearchOptimizerType); 00104 itkGetObjectMacro(LineSearchOptimizer, LineSearchOptimizerType); 00105 00107 itkGetConstMacro(MaximumNumberOfIterations, unsigned long); 00108 itkSetClampMacro(MaximumNumberOfIterations, unsigned long, 00109 1, NumericTraits<unsigned long>::max()); 00110 00116 itkGetConstMacro(GradientMagnitudeTolerance, double); 00117 itkSetMacro(GradientMagnitudeTolerance, double); 00118 00122 itkSetClampMacro(Memory,unsigned int,0,NumericTraits<unsigned int>::max()); 00123 itkGetConstMacro(Memory,unsigned int); 00124 00125 00126 protected: 00127 QuasiNewtonLBFGSOptimizer(); 00128 virtual ~QuasiNewtonLBFGSOptimizer(){}; 00129 00130 // \todo: should be implemented 00131 void PrintSelf(std::ostream& os, Indent indent) const {}; 00132 00133 DerivativeType m_CurrentGradient; 00134 MeasureType m_CurrentValue; 00135 unsigned long m_CurrentIteration; 00136 StopConditionType m_StopCondition; 00137 bool m_Stop; 00138 double m_CurrentStepLength; 00139 00141 bool m_InLineSearch; 00142 00143 RhoType m_Rho; 00144 SType m_S; 00145 YType m_Y; 00146 00147 unsigned int m_Point; 00148 unsigned int m_PreviousPoint; 00149 unsigned int m_Bound; 00150 00151 itkSetMacro(InLineSearch, bool); 00152 00157 virtual void ComputeDiagonalMatrix(DiagonalMatrixType & diag_H0); 00158 00165 virtual void ComputeSearchDirection( 00166 const DerivativeType & gradient, 00167 ParametersType & searchDir); 00168 00173 virtual void LineSearch( 00174 const ParametersType searchDir, 00175 double & step, 00176 ParametersType & x, 00177 MeasureType & f, 00178 DerivativeType & g ); 00179 00182 virtual void StoreCurrentPoint( 00183 const ParametersType & step, 00184 const DerivativeType & grad_dif); 00185 00190 virtual bool TestConvergence(bool firstLineSearchDone); 00191 00192 private: 00193 QuasiNewtonLBFGSOptimizer(const Self&); //purposely not implemented 00194 void operator=(const Self&); //purposely not implemented 00195 00196 unsigned long m_MaximumNumberOfIterations; 00197 double m_GradientMagnitudeTolerance; 00198 LineSearchOptimizerPointer m_LineSearchOptimizer; 00199 unsigned int m_Memory; 00200 00201 00202 }; // end class QuasiNewtonLBFGSOptimizer 00203 00204 00205 00206 } // end namespace itk 00207 00215 /* LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION */ 00216 /* JORGE NOCEDAL */ 00217 /* *** July 1990 *** */ 00218 00219 /* This subroutine solves the unconstrained minimization problem */ 00220 00221 /* min F(x), x= (x1,x2,...,xN), */ 00222 00223 /* using the limited memory BFGS method. The routine is especially */ 00224 /* effective on problems involving a large number of variables. In */ 00225 /* a typical iteration of this method an approximation Hk to the */ 00226 /* inverse of the Hessian is obtained by applying M BFGS updates to */ 00227 /* a diagonal matrix Hk0, using information from the previous M steps. */ 00228 /* The user specifies the number M, which determines the amount of */ 00229 /* storage required by the routine. The user may also provide the */ 00230 /* diagonal matrices Hk0 if not satisfied with the default choice. */ 00231 /* The algorithm is described in "On the limited memory BFGS method */ 00232 /* for large scale optimization", by D. Liu and J. Nocedal, */ 00233 /* Mathematical Programming B 45 (1989) 503-528. */ 00234 00235 /* The user is required to calculate the function value F and its */ 00236 /* gradient G. In order to allow the user complete control over */ 00237 /* these computations, reverse communication is used. The routine */ 00238 /* must be called repeatedly under the control of the parameter */ 00239 /* IFLAG. */ 00240 00241 /* The steplength is determined at each iteration by means of the */ 00242 /* line search routine MCVSRCH, which is a slight modification of */ 00243 /* the routine CSRCH written by More' and Thuente. */ 00244 00245 /* The calling statement is */ 00246 00247 /* CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG) */ 00248 00249 /* where */ 00250 00251 /* N is an INTEGER variable that must be set by the user to the */ 00252 /* number of variables. It is not altered by the routine. */ 00253 /* Restriction: N>0. */ 00254 00255 /* M is an INTEGER variable that must be set by the user to */ 00256 /* the number of corrections used in the BFGS update. It */ 00257 /* is not altered by the routine. Values of M less than 3 are */ 00258 /* not recommended; large values of M will result in excessive */ 00259 /* computing time. 3<= M <=7 is recommended. Restriction: M>0. */ 00260 00261 /* X is a DOUBLE PRECISION array of length N. On initial entry */ 00262 /* it must be set by the user to the values of the initial */ 00263 /* estimate of the solution vector. On exit with IFLAG=0, it */ 00264 /* contains the values of the variables at the best point */ 00265 /* found (usually a solution). */ 00266 00267 /* F is a DOUBLE PRECISION variable. Before initial entry and on */ 00268 /* a re-entry with IFLAG=1, it must be set by the user to */ 00269 /* contain the value of the function F at the point X. */ 00270 00271 /* G is a DOUBLE PRECISION array of length N. Before initial */ 00272 /* entry and on a re-entry with IFLAG=1, it must be set by */ 00273 /* the user to contain the components of the gradient G at */ 00274 /* the point X. */ 00275 00276 /* DIAGCO is a LOGICAL variable that must be set to .TRUE. if the */ 00277 /* user wishes to provide the diagonal matrix Hk0 at each */ 00278 /* iteration. Otherwise it should be set to .FALSE., in which */ 00279 /* case LBFGS will use a default value described below. If */ 00280 /* DIAGCO is set to .TRUE. the routine will return at each */ 00281 /* iteration of the algorithm with IFLAG=2, and the diagonal */ 00282 /* matrix Hk0 must be provided in the array DIAG. */ 00283 00284 /* DIAG is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE., */ 00285 /* then on initial entry or on re-entry with IFLAG=2, DIAG */ 00286 /* it must be set by the user to contain the values of the */ 00287 /* diagonal matrix Hk0. Restriction: all elements of DIAG */ 00288 /* must be positive. */ 00289 00290 /* IPRINT is an INTEGER array of length two which must be set by the */ 00291 /* user. */ 00292 00293 /* IPRINT(1) specifies the frequency of the output: */ 00294 /* IPRINT(1) < 0 : no output is generated, */ 00295 /* IPRINT(1) = 0 : output only at first and last iteration, */ 00296 /* IPRINT(1) > 0 : output every IPRINT(1) iterations. */ 00297 00298 /* IPRINT(2) specifies the type of output generated: */ 00299 /* IPRINT(2) = 0 : iteration count, number of function */ 00300 /* evaluations, function value, norm of the */ 00301 /* gradient, and steplength, */ 00302 /* IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of */ 00303 /* variables and gradient vector at the */ 00304 /* initial point, */ 00305 /* IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of */ 00306 /* variables, */ 00307 /* IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.*/ 00308 00309 00310 /* EPS is a positive DOUBLE PRECISION variable that must be set by */ 00311 /* the user, and determines the accuracy with which the solution*/ 00312 /* is to be found. The subroutine terminates when */ 00313 00314 /* ||G|| < EPS max(1,||X||), */ 00315 00316 /* where ||.|| denotes the Euclidean norm. */ 00317 00318 /* XTOL is a positive DOUBLE PRECISION variable that must be set by */ 00319 /* the user to an estimate of the machine precision (e.g. */ 00320 /* 10**(-16) on a SUN station 3/60). The line search routine will*/ 00321 /* terminate if the relative width of the interval of uncertainty*/ 00322 /* is less than XTOL. */ 00323 00324 /* W is a DOUBLE PRECISION array of length N(2M+1)+2M used as */ 00325 /* workspace for LBFGS. This array must not be altered by the */ 00326 /* user. */ 00327 00328 /* IFLAG is an INTEGER variable that must be set to 0 on initial entry*/ 00329 /* to the subroutine. A return with IFLAG<0 indicates an error, */ 00330 /* and IFLAG=0 indicates that the routine has terminated without*/ 00331 /* detecting errors. On a return with IFLAG=1, the user must */ 00332 /* evaluate the function F and gradient G. On a return with */ 00333 /* IFLAG=2, the user must provide the diagonal matrix Hk0. */ 00334 00335 /* The following negative values of IFLAG, detecting an error, */ 00336 /* are possible: */ 00337 00338 /* IFLAG=-1 The line search routine MCSRCH failed. The */ 00339 /* parameter INFO provides more detailed information */ 00340 /* (see also the documentation of MCSRCH): */ 00341 00342 /* INFO = 0 IMPROPER INPUT PARAMETERS. */ 00343 00344 /* INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF */ 00345 /* UNCERTAINTY IS AT MOST XTOL. */ 00346 00347 /* INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE */ 00348 /* REQUIRED AT THE PRESENT ITERATION. */ 00349 00350 /* INFO = 4 THE STEP IS TOO SMALL. */ 00351 00352 /* INFO = 5 THE STEP IS TOO LARGE. */ 00353 00354 /* INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.*/ 00355 /* THERE MAY NOT BE A STEP WHICH SATISFIES */ 00356 /* THE SUFFICIENT DECREASE AND CURVATURE */ 00357 /* CONDITIONS. TOLERANCES MAY BE TOO SMALL. */ 00358 00359 /* IFLAG=-2 The i-th diagonal element of the diagonal inverse */ 00360 /* Hessian approximation, given in DIAG, is not */ 00361 /* positive. */ 00362 00363 /* IFLAG=-3 Improper input parameters for LBFGS (N or M are */ 00364 /* not positive). */ 00365 00366 00367 /* ON THE DRIVER: */ 00368 00369 /* The program that calls LBFGS must contain the declaration: */ 00370 00371 /* EXTERNAL LB2 */ 00372 00373 /* LB2 is a BLOCK DATA that defines the default values of several */ 00374 /* parameters described in the COMMON section. */ 00375 00376 /* COMMON: */ 00377 00378 /* The subroutine contains one common area, which the user may wish to */ 00379 /* reference: */ 00380 00381 /* awf added stpawf */ 00382 00383 /* MP is an INTEGER variable with default value 6. It is used as the */ 00384 /* unit number for the printing of the monitoring information */ 00385 /* controlled by IPRINT. */ 00386 00387 /* LP is an INTEGER variable with default value 6. It is used as the */ 00388 /* unit number for the printing of error messages. This printing */ 00389 /* may be suppressed by setting LP to a non-positive value. */ 00390 00391 /* GTOL is a DOUBLE PRECISION variable with default value 0.9, which */ 00392 /* controls the accuracy of the line search routine MCSRCH. If the */ 00393 /* function and gradient evaluations are inexpensive with respect */ 00394 /* to the cost of the iteration (which is sometimes the case when */ 00395 /* solving very large problems) it may be advantageous to set GTOL */ 00396 /* to a small value. A typical small value is 0.1. Restriction: */ 00397 /* GTOL should be greater than 1.D-04. */ 00398 00399 /* STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which */ 00400 /* specify lower and uper bounds for the step in the line search. */ 00401 /* Their default values are 1.D-20 and 1.D+20, respectively. These */ 00402 /* values need not be modified unless the exponents are too large */ 00403 /* for the machine being used, or unless the problem is extremely */ 00404 /* badly scaled (in which case the exponents should be increased). */ 00405 00406 00407 /* MACHINE DEPENDENCIES */ 00408 00409 /* The only variables that are machine-dependent are XTOL, */ 00410 /* STPMIN and STPMAX. */ 00411 00412 00413 /* GENERAL INFORMATION */ 00414 00415 /* Other routines called directly: DAXPY, DDOT, LB1, MCSRCH */ 00416 00417 /* Input/Output : No input; diagnostic messages on unit MP and */ 00418 /* error messages on unit LP. */ 00419 00420 00421 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/ 00422 00423 00424 #endif //#ifndef __itkQuasiNewtonLBFGSOptimizer_h 00425
Generated on 11-05-2012 for elastix by ![]() |
![]() |