If you are looking for a __dense__ Levenberg-Marquardt C/C++ implementation,
please have a look at levmar.

# Introduction

This site provides `sparseLM`, a general-purpose software package for large-scale, arbitrarily
sparse non-linear least squares
that is distributed under the GNU General Public License.
`sparseLM` implements a sparse variant of the
Levenberg-Marquardt optimization algorithm
in ANSI C and includes a Matlab MEX-interface.

# Motivation

Non-linear model fitting
is a standard procedure for inferring the values of parameters from observed data in a manner that
reduces the impact of measurement error. This is achieved by employing non-linear least squares to
minimize the total error pertaining to overdetermined sets of observations. The total error is expressed
by a sum-of-squares cost function (i.e., a $L$_{2} norm), whose minimizer represents
the statistically optimal estimate of the sought parameters under the assumption of Gaussian noise.
Owing to their non-convexity, $L$_{2} cost functions are minimized
with iterative non-linear optimization techniques, of which the Levenberg-Marquardt (LM) algorithm
has become the de facto standard. LM operates by repeatedly linearizing the function to be minimized
in the neighborhood of the current minimizer estimate and computing an improvement to it through the
solution of a linear system defined with the aid of the Jacobian
and known as the *normal equations*.
Considering that each computation of the solution to a dense linear system has complexity
$O(N3)$ in the number of unknown parameters $N$, it is clear that
general purpose LM implementations are computationally very demanding when employed to minimize
functions involving a large number of parameters.

Fortunately, when dealing with large estimation problems arising in practice, the corresponding error typically exhibits lack of interdependence among certain subgroups of the parameters to be estimated. This observation translates to Jacobians for the least squares minimization that are sparse, that is, consist of mostly zero elements. In turn, sparse Jacobians yield normal equation systems with sparse block structure. By avoiding storing and operating on zero elements of the normal matrix during the course of LM, substantial memory and execution time benefits can be gained. The downside is that taking advantage of arbitrary sparseness is problem-specific and requires considerable implementation effort for developing LM variants adapted to particular sparse problems.

# Overview

During the last few years, effective mechanisms for dealing with arbitrary sparseness have emerged in the form of a number of
algorithms and corresponding implementations for the direct solution
of large sparse linear systems of equations. Owing to a number of attractive properties (no need for preconditioners,
fast convergence, computation of exact rather than approximate solutions, mature implementations), sparse linear solvers are
well-suited as general-purpose sparse linear solvers. `sparseLM` builds upon this technology to fulfill the need for
a high quality sparse Levenberg-Marquardt solver designed for general-purpose, arbitrarily sparse non-linear least squares
minimization. The software has been designed with the
twofold objective of exploiting sparseness for maximizing performance while shielding the user from the algorithmic
details associated with direct linear solvers. It accepts sparse Jacobians encoded in either compressed row storage
(CRS) or compressed column
storage (CCS)
(a.k.a. *Harwell-Boeing*) format,
allowing user applications to choose the representation that is most natural to them. It also offers the possibility
of numerically approximating the Jacobian using forward finite differences on data provided by successive invocations
of the function to be minimized. In this case, only the sparsity pattern of the Jacobian should be specified by the
user and not its numerical values. Although not recommended due to performance considerations, `sparseLM` can
even be instructed to attempt the automatic detection of the sparsity pattern corresponding to the Jacobian of the
function to be minimized. The core LM algorithm employed by `sparseLM` is identical to that in the
`levmar` __dense__ LM package which is available here.
As a result, `sparseLM` has a similar "look and feel" with `levmar`.
More details on `sparseLM` can be found in this ECCV10 paper
(bibtex entry).

`sparseLM` includes interfaces to a wide variety of direct sparse codes, the list of which currently consists of
HSL's MA57,
MA77, MA47 & MA27,
PARDISO,
SuperLU,
TAUCS,
UMFPACK,
CHOLMOD,
LDL,
CSparse,
SPOOLES
and
MUMPS.

A few of the above solvers (e.g. PARDISO, MUMPS) use shared-memory parallelism to improve numerical performance,
whereas others (such as MA77 & TAUCS) support out-of-core processing for dealing with very large
matrices that cannot fit in RAM.
CHOLMOD, an open-source sparse Cholesky factorization package, is chosen as `sparseLM`'s default solver.
Header files for CHOLMOD and libraries compiled with MSVC 2005 under WinXP are
here.

To the best of my knowledge, `sparseLM` is the first and currently the only package
of its kind to be released as free software.

# Available Functions

`sparseLM` offers several user-callable functions whose suffix denotes whether the Jacobian is analytic
or approximate and specifies the sparse format it is stored in. More specifically, `sparseLM` includes
the functions below:

`sparselm_dercrs()`: analytic Jacobian supplied in CRS`sparselm_derccs()`: analytic Jacobian supplied in CCS`sparselm_difcrs()`: finite difference approximated Jacobian, zero pattern supplied in CRS`sparselm_difccs()`: finite difference approximated Jacobian, zero pattern supplied in CCS

Both CRS and CCS employ C's native zero-based indexing convention. Notice that using finite differences to approximate the Jacobian results in repetitive evaluations of the function to be fitted, thus negatively influencing performance. To reduce the total number of these evaluations, a finite difference approximation to a Jacobian is estimated using a scheme that exploits sparseness to compute several of the Jacobian's columns with a single function evaluation. All functions basically 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.

Briefly, the arguments of the functions are as follows: pointers to routines evaluating the vector function f and its
Jacobian J, pointers to the initial estimate of the parameter vector p and the measurement vector x, the
dimension of p, the number of initial elements of p that are to be kept fixed to their initial values, the dimension of x,
the number of nonzero elements of J, the number of nonzero elements of the approximate Hessian J^t*J, the maximum number of
iterations, a pointer to an 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, 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 `splmdemo.c` sample program gives working examples of minimizing several nonlinear functions.
The exact prototypes of `sparselm_dercrs()` and `sparselm_difcrs()` are given next, along with a description of
each argument. I and O in these descriptions denote input and output arguments, respectively.
Being similar, the prototypes of `sparselm_derccs()` and `sparselm_difccs()` are omitted.

`sparselm_dercrs()`

```
.:: toggle display ::.
```

`sparselm_difcrs()`

```
.:: toggle display ::.
```

Note that passing `NULL` as the value of the `opts` argument of both functions results in default settings being
used. Consult Madsen et al's lecture notes
for details on the roles of `opts` and `info` arguments.
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.
For CCS Jacobians, functions `sparselm_derccs()` and `sparselm_difccs()` can be used.

In addition to the above, `sparseLM` also offers some utility functions that often prove handy.
These include functions `splm_ccsm2crsm()` & `splm_crsm2ccsm()` which implement conversion
of a sparse matrix between CRS and CCS formats. Furthermore, function `splm_stm2ccsm()` allows a Jacobian
that has been constructed element-by-element in sparse triplet
format to be converted to CCS for use by `sparselm_derccs()`/`sparselm_difccs()`.

# Contact Address

If you find this package useful or have any comments/questions/suggestions, please contact me at

Be warned that although I try to reply to most messages, I might take long to do so.

In case that you use `sparseLM` in your published work, please include a reference to this paper:
[ bibtex entry ].

hits since Sat Dec 4 11:45:23 EET 2010