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.
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 doxygen 1.7.6.1 elastix logo