Manolis Lourakis
Institute of Computer Science,
Foundation for Research and Technology - Hellas,
Heraklion, Crete, Greece
Last updated May 9, 2008
A sparse variant of the Levenberg-Marquardt algorithm implemented by levmar has been applied to bundle adjustment, a computer vision/photogrammetry problem that typically involves several thousand variables; please have a look at sba for more details.
This page provides GPL native ANSI C implementations of the Levenberg-Marquardt optimization algorithm, usable also from C++, Matlab, Perl and Python, and explains their use. Both unconstrained and constrained (under linear equations and box constraints) Levenberg-Marquardt variants are included. The Levenberg-Marquardt (LM) algorithm is an iterative technique that finds a local minimum of a function that is expressed as the sum of squares of nonlinear functions. It has become a standard technique for nonlinear least-squares problems and can be thought of as a combination of steepest descent and the Gauss-Newton method. When the current solution is far from the correct one, the algorithm behaves like a steepest descent method: slow, but guaranteed to converge. When the current solution is close to the correct solution, it becomes a Gauss-Newton method.
The lmder routine from Minpack, implemented in the early '80s at the Argonne National Lab, is perhaps the most widely used free implementation of the LM algorithm. lmder is written in FORTRAN77 and over the years has proven to be a reliable piece of software. Considering that FORTRAN routines can be called from C/C++, one might wonder about the motivation for writing a version of LM in C. Well, the problem is that when FORTRAN is called from C, the programmer should be aware of (and conform to) several rules regarding name mangling, argument passing, multidimensional array memory layout, linkage conventions, etc, that are unnatural compared to ordinary C rules. A second reason is that this approach takes for granted that a FORTRAN compiler for the target programming environment is available, which might not necessarily be the case. Another reason has to do with failure to understand the inner workings of a FORTRAN implementation: Occasionally, when it is necessary to precisely understand what the FORTRAN code does, certain pieces of it might seem incomprehensible to programmers without any knowledge of FORTRAN. Automatic FORTRAN to C translators (e.g. f2c) do not solve the problem since the produced C code is pretty illegible to “uninitiated” humans. Moreover, documentation describing the mathematics upon which the implementation is based might be obscure or inaccessible. Last but not least, a candidate LM implementation in C should be free and technically sound. For example, the C variant of the LM algorithm presented in the “Numerical Recipes” book (i.e. mrqmin), is not always a viable choice: Besides its being copyrighted, it is reputed to lack robustness.
For the above reasons, I have developed the levmar package which includes C implementations of LM flavors that are also usable with C++. levmar includes double and single precision LM implementations, both with analytic and finite difference approximated Jacobians. It is provided free of charge, under the terms of the GNU General Public License. The mathematical theory behind unconstrained levmar is described in detail in the lecture notes entitled Methods for Non-Linear Least Squares Problems, by K. Madsen, H.B. Nielsen and O. Tingleff, Technical University of Denmark; Matlab implementations of the algorithms presented in the lecture notes are also available. Note however that the formulation of the minimization problem adopted here is slightly different from that described in the lecture notes. I have also written a short note, providing a quick overview of the material in the lecture notes.
To deal with linear equation constraints, levmar employs variable elimination based on QR factorization, as described in ch. 15 of the book Numerical Optimization by Nocedal and Wright. For the box-constrained case, levmar implements the algorithm proposed by C. Kanzow, N. Yamashita and M. Fukushima, Levenberg-Marquardt methods for constrained nonlinear equations with strong local convergence properties, Journal of Computational and Applied Mathematics 172, 2004, pp. 375-397.
levmar provides the following two options regarding the solution of the linear systems formed by the augmented normal equations:
Use of LAPACK or an equivalent vendor library is strongly recommended; if none is installed at your site I suggest getting CLAPACK, the f2c'ed version of LAPACK. MSWin users have the additional option of downloading precompiled LAPACK/BLAS libraries from here. You might also consider combining CLAPACK with ATLAS, (the automatically tuned BLAS implementation) or GotoBLAS (K. Goto's high-performance BLAS). On the other hand, LAPACK's use is not mandatory for compiling the unconstrained routines and the second option makes levmar totally self-contained. The linear equation constrained routines, however, cannot be compiled without LAPACK.
Interfaces for using levmar from high-level programming environments such as Matlab, Perl and Python are also available; please refer to the FAQ for more details.
Change LogApr. 30, 2008: A CMakeLists.txt script supports cross-platform compilation using the CMake project/makefile generator.
Apr. 29, 2008: levmar has been integrated into Astrometry.net.
Mar. 26, 2008: levmar has been integrated into StochFit.
Dec. 17, 2007: levmar can now handle simultaneous box & linear constraints.
Dec. 17, 2007: levmar now has a matlab MEX interface. This might also work with GNU octave but it has not been tested.
Nov. 27, 2006: Alastair Tse has written a Python interface to levmar.
Jun. 30, 2006: John Lapeyre has written a PDL (perl) interface to levmar.
Dec. 08, 2005: levmar has been integrated into SciGraphica, CellML, Stimfit, Autopano.
Feb. 11, 2005: Linear and box constraints support, covariance of a fit estimated, new routines xlevmar_chkjac() for verifying a Jacobian.
Feb. 11, 2005: Attention old users: Compared to versions 1.2 or earlier, levmar v. 2.0 and later contain minor changes to the API, see the change log.
levmar offers several user-callable functions obeying the following naming convention: The first letter (d or s) specifies double or single precision and the suffix (_der or _dif) denotes analytic or approximate Jacobian. If present, the lec, bc and blec components imply linear equation, box and simultaneous box and linear equation constraints, respectively. More specifically, levmar includes the functions below:
Notice that using finite differences to approximate the Jacobian results in repetitive evaluations of the function to be fitted. Aiming to reduce the total number of these evaluations, the xxxxxxx_dif functions implement secant approximations to the Jacobian using Broyden's rank one updates. All functions solve the same problem, i.e. they seek the parameter vector p that best describes (in terms of the L2 norm) the measurements vector x. More precisely, given a vector function f : R^m --> R^n with n>=m, they compute a p such that f(p) ~= x, i.e. the squared norm ||e||^2=||x-f(p)||^2 is minimized. Also, box constraints of the form lb[i]<=p[i]<=ub[i] and linear constrains of the form A*p=b with A kxm and b kx1, can be enforced on the minimizer. More details on the unconstrained LM algorithm implemented by levmar can be found in this note.
Briefly, the arguments of the functions are as follows: pointers to routines evaluating the vector function f and its Jacobian (if applicable), pointers to the initial estimate of the parameter vector p and the measurement vector x, the dimension m of p, the dimension n of x, the maximum number of iterations, a pointer to a 3 element array whose elements specify certain minimization options, a pointer to an array into the elements of which various pieces of information regarding the result of minimization will be returned, a pointer to working memory, a pointer to a covariance matrix and a pointer to application-specific data structures that might be necessary for computing the function to be minimized and its derivatives. More details can be found below. Additionally, the included lmdemo.c sample program gives working examples of minimizing several nonlinear functions. The exact prototypes of dlevmar_der(), dlevmar_dif() are given next, along with a description of each argument. slevmar_der() and slevmar_dif() are identical, except that doubles are changed to floats. I and O in the descriptions below denote input and output arguments respectively. For descriptions of the arguments to the constrained optimization routines, see the comments in the code.
/*
* This function seeks the mx1 parameter vector p that best describes the nx1 measurements vector x.
* All computations are double precision.
* An analytic Jacobian is required. In case the latter is unavailable or expensive to compute,
* use dlevmar_dif() below.
*
* Returns the number of iterations (>=0) if successful, -1 if failed
*
*/
int dlevmar_der(
void (*func)(double *p, double *hx, int m, int n, void *adata), /* functional relation describing measurements.
* A p \in R^m yields a \hat{x} \in R^n
*/
void (*jacf)(double *p, double *j, int m, int n, void *adata), /* function to evaluate the Jacobian \part x / \part p */
double *p, /* I/O: initial parameter estimates. On output contains the estimated solution */
double *x, /* I: measurement vector. NULL implies a zero vector */
int m, /* I: parameter vector dimension (i.e. #unknowns) */
int n, /* I: measurement vector dimension */
int itmax, /* I: maximum number of iterations */
double opts[4], /* I: minim. options [\tau, \epsilon1, \epsilon2, \epsilon3]. Respectively the scale factor for initial \mu,
* stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2. Set to NULL for defaults to be used
*/
double info[LM_INFO_SZ],
/* O: information regarding the minimization. Set to NULL if don't care
* info[0]= ||e||_2 at initial p.
* info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, \mu/max[J^T J]_ii ], all computed at estimated p.
* info[5]= # iterations,
* info[6]=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased \mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
* info[7]= # function evaluations
* info[8]= # Jacobian evaluations
*/
double *work, /* I: pointer to working memory, allocated internally if NULL. If !=NULL, it is assumed to point to
* a memory chunk at least LM_DER_WORKSZ(m, n)*sizeof(double) bytes long
*/
double *covar, /* O: Covariance matrix corresponding to LS solution; Assumed to point to a mxm matrix.
* Set to NULL if not needed.
*/
void *adata) /* I: pointer to possibly needed additional data, passed uninterpreted to func & jacf.
* Set to NULL if not needed
*/
/*
* Similar to dlevmar_der() except that the Jacobian is approximated internally with the aid of finite differences.
* Broyden's rank one updates are used to compute secant approximations to the Jacobian, effectively avoiding to call
* func several times for computing the finite difference approximations.
* If the analytic Jacobian is available, use dlevmar_der() above.
*
* Returns the number of iterations (>=0) if successful, -1 if failed
*
*/
int dlevmar_dif(
void (*func)(double *p, double *hx, int m, int n, void *adata), /* functional relation describing measurements.
* A p \in R^m yields a \hat{x} \in R^n
*/
double *p, /* I/O: initial parameter estimates. On output contains the estimated solution */
double *x, /* I: measurement vector. NULL implies a zero vector */
int m, /* I: parameter vector dimension (i.e. #unknowns) */
int n, /* I: measurement vector dimension */
int itmax, /* I: maximum number of iterations */
double opts[5], /* I: opts[0-4] = minim. options [\tau, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively the
* scale factor for initial \mu, stopping thresholds for ||J^T e||_inf, ||Dp||_2 and ||e||_2 and the
* step used in difference approximation to the Jacobian. If \delta<0, the Jacobian is approximated
* with central differences which are more accurate (but slower!) compared to the forward differences
* employed by default. Set to NULL for defaults to be used.
*/
double info[LM_INFO_SZ],
/* O: information regarding the minimization. Set to NULL if don't care
* info[0]= ||e||_2 at initial p.
* info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, \mu/max[J^T J]_ii ], all computed at estimated p.
* info[5]= # iterations,
* info[6]=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased \mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
* info[7]= # function evaluations
* info[8]= # Jacobian evaluations
*/
double *work, /* I: working memory, allocated internally if NULL. If !=NULL, it is assumed to point to
* a memory chunk at least LM_DIF_WORKSZ(m, n)*sizeof(double) bytes long
*/
double *covar, /* O: Covariance matrix corresponding to LS solution; Assumed to point to a mxm matrix.
* Set to NULL if not needed.
*/
void *adata) /* I: pointer to possibly needed additional data, passed uninterpreted to func.
* Set to NULL if not needed
*/
Note that passing NULL as the value of the last four arguments of both functions results in default values being used. Consult Madsen et al's lecture notes for detailed explanation of the roles of opts and info arguments. Argument work is meant to be used to avoid the overhead associated with repeatedly allocating and freeing memory chunks of the same size, when several minimizations of the same dimension are carried out. Instead of letting the functions allocate the required memory themselves, you can use the macros LM_DER_WORKSZ and LM_DIF_WORKSZ to calculate the appropriate memory size, allocate it and pass it to each invocation. lmdemo.c provides such an example. The last argument (i.e. adata) is intended to help avoid direct use of globals in the routines computing the function to be minimized and its Jacobian. A structure containing pointers to appropriate data structures can be set up and a pointer to it can be passed to the LM function which then passes it uninterpreted to each call of the user-supplied routines.
The source code is distributed in a gzip'ed tar file. It has been tested under Linux with gcc/icc/ lcc, under cygwin with gcc and under Windows with MSVC. Remember to include lm.h before calling any of the routines from your C/C++ programs. Refer to lmdemo.c for working examples of minimizing different nonlinear functions. The core functionality of levmar is provided by source files Axb_core.c, lm_core.c, lmbc_core.c, lmlec_core.c, lmblec_core.c and misc_core.c. With the aid of appropriate definitions, these files are included twice by Axb.c, lm.c, lmbc.c, lmlec.c, lmblec.c and misc.c, respectively, to provide single and double precision LM versions. Thus, when building levmar, the *_core.c files should not be compiled directly.
A mathematical explanation, plus pseudocode, of the unconstrained LM algorithm implemented by levmar can be found in this note.
For historical reasons, older versions of levmar are also available:
|
|
This content is made available under the GNU Free Documentation License.
|