C/C++ Minpack

C/C++ Minpack source code on GitHub

What is Minpack?

This is the official description of Minpack, from the original ReadMe file: Minpack includes software for solving nonlinear equations and nonlinear least squares problems. Five algorithmic paths each include a core subroutine and an easy-to-use driver. The algorithms proceed either from an analytic specification of the Jacobian matrix or directly from the problem functions. The paths include facilities for systems of equations with a banded Jacobian matrix, for least squares problems with a large amount of data, and for checking the consistency of the Jacobian matrix with the functions.

The original authors of the FORTRAN version are Jorge More', Burt Garbow, and Ken Hillstrom from Argonne National Laboratory, and the code can be obtained from Netlib.

Minpack is probably the best open-source implementation of the Levenberg-Marquardt algorithm (in fact, it is even better, since it adds to L-M automatic variables scaling). There is another open-source L-M implementation in C/C++, levmar by Manolis Lourakis, but unfortunately is is released under the GPL, which restricts its inclusion in commercial software. Minpack is licensed under a BSD-like license (available in the distribution).

What about CMinpack?

In July 2002 (before levmar), Manolis Lourakis (lourakis at ics forth gr) released a C version of Minpack, called CMinpack, obtained from the FORTRAN version using f2c and some limited manual editing. However, this version had several problems, which came from the FORTRAN version:

  1. All the function prototypes were following the original FORTRAN call conventions, so that all parameters are passed by reference (if a function needs an int parameter, you have to pass a pointer to this int).
  2. There were lots of static variables in the code, thus you could not optimize a function which required calling Minpack to be evaluated (a minimization-of-minimization problem for example): The Minpack code is not reentrant.
  3. If the function to be optimized has to use extra parameters or data (this is the case most of the time), the only way to access them was though global variables, which is very bad, especially if you want to use the same function with different data in different threads: The Minpack code is not MT-Safe.
  4. There was no C/C++ include file.
  5. Examples and tests were missing from the distribution, although there are some FORTRAN examples in the documentation.

Why is C/C++ Minpack better?

I took a dozen of hours to rework all these problems, and came out with a pure C version of Minpack, with has standard (ISO C99) parameters passing, is fully reentrant, multithread-safe, and has a full set of examples and tests:

  1. Input variables are now passed by value, output variables are passed by reference. The keyword "const" is used as much as possible for constant arrays. The return value of each function is now used to get the function status (it was obtained via the IFLAG or INFO parameter in Minpack).
  2. All non-const static variables were removed, and the code was tested after that. Luckily, Minpack didn't use the nastiest feature in FORTRAN: all local variables are static, so that a function can behave differently when you call it several times.
  3. The function to be minimized and all the Minpack functions now take an extra "void*" argument, which can be used to pass any pointer-to-struct or pointer-to-class, and you can put all you extra parameters and data in that struct. Just cast this pointer to the appropriate pointer type in your function, and there they are! There is no need for global variables anymore. Be careful if you access the same object from different threads, though (a solution is to protect this extra data with a mutex).
  4. The Debian project did a C include file for Minpack. It still needed some work (add consts and C++ compatibility), so I did this work, and used the include file for the FORTRAN version as the base for my C/C++ version.
  5. The Debian project also translated all the FORTRAN examples to C. I worked from these to produce examples which also call my C/C++ version of Minpack instead of the FORTRAN version. Also included in the distribution are reference output files produced by the test runs (for comparison).

If you use C/C++ Minpack for a publication, you should cite it as:

@misc{cminpack,
  title={C/C++ Minpack},
  author={Devernay, Fr{\'e}d{\'e}ric},
  year={2007},
  howpublished = "\url{http://devernay.github.io/cminpack}",
}

Distribution

The distribution contains:

It is distributed under the original Minpack license (see the file CopyrightMINPACK.txt in the distribution).

Download

All releases are listed on the GitHub releases page. Versions 1.3.5 and later link to the GitHub source archive; versions up to 1.3.4 link to the original distribution tarball attached to each release.

1.3.14 (latest) · 1.3.13 · 1.3.12 · 1.3.11 · 1.3.10 · 1.3.9 · 1.3.8 · 1.3.7 · 1.3.6 · 1.3.5 · 1.3.4 · 1.3.3 · 1.3.2 · 1.3.1 · 1.3.0 · 1.2.2 · 1.2.1 · 1.2.0 · 1.1.5 · 1.1.4 · 1.1.3 · 1.1.2 · 1.1.1 · 1.0.4 · 1.0.3 · 1.0.2 · 1.0.1

GitHub repository: https://github.com/devernay/cminpack

Building CMinpack

CMinpack can be built with CMake. By default, CMake will build a single precision library named cminpacks, and a double precision library named cminpack. You can choose to build only one of the single, double, or long double precision variants by setting CMINPACK_PRECISION to one of "s", "d", or "ld" respectively (e.g. cmake -DCMINPACK_PRECISION=d ..).

The source distribution also contains a Makefile which can be used to build cminpack for double precision (the default, make or make double) for single-precision (make float), half-precision (make half), long double precision (make longdouble), and even CUDA (make cuda), or using LAPACK for linear algebra (make lapack, double or single precision only). Unit tests are also provided, but results may depend on the platform (make checkdouble checkfloat checkhalf checklongdouble checklapack).

Using CMinpack

The CMinpack calls have the same name as the FORTRAN functions, in lowercase (e.g. lmder(...)). See the links to the documentation below, or take a look at the simple examples in the examples directory of the distribution. The simple examples are named after the function they call: tlmder.c is the simple example for lmder.

If you want to use the single precision CMinpack, you should define __cminpack_float__ before including cminpack.h. __cminpack_half__ has to be defined for the half-precision version (and the code needs to be compiled with a C++ compiler).

The single-precision versions of the functions are prefixed by "s" (as in "slmder(...)"), and the half-precision are prefixed by "h".

CMinpack defines __cminpack_real__ as the floating point type, and the __cminpack_func__() macro can be used to call CMinpack functions independently of the precision used (as in the examples). However, you shouldn't use these macros in your own code, since your code is probably designed for a specific precision, and you should prefer calling directly lmder(...) or slmder(...).

C++ bindings

include/cminpackcpp.hpp is a header-only C++ wrapper (issue #74). The C API takes the user callback as a plain function pointer plus an opaque void *p for user data, which a capturing lambda, a functor, or a std::function cannot satisfy. The wrapper adds overloads in namespace cminpack that accept any callable whose signature matches the C callback minus the leading void *p; the callable is forwarded through the void *p slot by a small template trampoline, so no global state is needed. It adds no allocation and no virtual dispatch, and the callable only has to outlive the solver call.

The wrapped entry points are hybrd1/hybrd and hybrj1 (nonlinear equations), and lmdif1/lmdif, lmder1/lmder and lmstr1/lmstr (least squares). Every argument after the callable is identical to the C function — you just drop the (fcn, p) pair and pass your callable first:

#include <cminpackcpp.hpp>
#include <vector>

std::vector<double> y = /* measured data */;

// signature = the C callback minus the leading void *p
auto residual = [&](int m, int n, const double *x, double *fvec, int iflag) {
    for (int i = 0; i < m; ++i)
        fvec[i] = model(x, i) - y[i];   // captures y, no global needed
    return 0;
};

int info = cminpack::lmdif1(residual, m, n, x, fvec, tol, iwa, wa, lwa);

A class member function works too, as long as it is wrapped in a callable that supplies the instance — a lambda capturing this (or the object) is the simplest, and std::bind/std::function also work. A bare pointer-to-member (&Class::residual) is not itself callable, so it must be wrapped:

struct LineModel {
    std::vector<double> y;
    int residual(int m, int n, const double *x, double *fvec, int iflag) const {
        for (int i = 0; i < m; ++i) fvec[i] = model(x, i) - y[i];
        return 0;
    }
};

LineModel obj = /* ... */;
auto cb = [&obj](int m, int n, const double *x, double *fvec, int iflag) {
    return obj.residual(m, n, x, fvec, iflag);      // supplies the instance
};
int info = cminpack::lmdif1(cb, m, n, x, fvec, tol, iwa, wa, lwa);
// std::bind is equivalent:
//   auto cb = std::bind(&LineModel::residual, &obj, _1, _2, _3, _4, _5);

The header is precision-agnostic: like cminpack.h it selects the real type via __cminpack_real__, so defining __cminpack_float__ (or the long-double macro) before including it targets that variant — link against the matching cminpack library. A worked example exercising a capturing lambda, a stateful functor and a std::function is in examples/tcppwrap.cpp.

Testing

Two test routes are available:

The cross-check (examples/crosscheck.py, and the CMake crosscheck_* tests) is a regression gate against the original FORTRAN MINPACK: on the intensive driver problems, cminpack must converge on every problem FORTRAN converges on. It compares against committed FORTRAN reference outputs (examples/ref/*.fortran.ref), so no Fortran compiler is needed. Where FORTRAN itself does not converge (the problems pushed from 10x/100x starting points), a different result is accepted — pure-C and f2c, like different compilers, take different (equally valid) iteration paths there, for the reasons in the next section. Iteration-count and last-digit differences never fail the build.

Python 3 is optional. It is not needed to build cminpack, nor for the standard or smoke tests (which use cmpfiles); only the cross-check uses it. When Python 3 is missing, both build systems skip only that cross-check — printing a disclaimer that coverage is incomplete — and run everything else.

Numerical differences from FORTRAN MINPACK

On the difficult test problems (the Moré/Garbow/Hillstrom functions exercised by the intensive driver programs examples/*drv*), cminpack and the original FORTRAN MINPACK can report different iteration, function- and Jacobian-evaluation counts, and last-digit differences in the results. This is expected and harmless: the problems still converge to equally valid solutions. Test messages reporting such differences are normal, not a sign of a broken build.

The differences do not come from the algorithm. dpmpar returns identical machine constants on both sides, so the convergence tolerance is the same, and cminpack's own pure C (src/) is a cleaned-up rewrite of the f2c output (src/f2c/) — so even these two can reach different (equally valid) results on the hardest problems when the compiler contracts FMAs differently in each. The dominant cause is floating-point contraction: gcc and gfortran fuse a*b + c into a fused multiply-add (FMA) at different places, so intermediate values differ by one unit in the last place (ULP). On the ill-conditioned problems a single ULP early in the iteration can flip a trust-region accept/reject decision, and the two runs then follow different paths to (equally valid) results. Compiling both sides with -ffp-contract=off removes most of the divergence.

(The enorm scaling constants rdwarf/rgiant, which cminpack tunes to the IEEE range rather than MINPACK's 1980 values, are a separate and minor effect: they only change the norm of out-of-range vectors in the last bit. See src/enorm.c for the full explanation.)

Because of this, examples/crosscheck.py and the CMake crosscheck_* tests are a regression gate against the original FORTRAN MINPACK: on the intensive driver problems, cminpack must converge on every problem FORTRAN converges on (compared against the committed examples/ref/*.fortran.ref, so no Fortran compiler is needed). Following the original MINPACK drivers, a problem is gated only where FORTRAN converged (exit parameter info < 5); where FORTRAN itself does not converge (info >= 5) a different result is accepted. On a gated problem the build fails if either its own solver reports non-convergence (info >= 5) or FORTRAN drove the residual to ~0 but this build did not (which catches a genuine divergence the exit code alone would miss). Iteration-count and last-digit differences are reported for information only and never fail the build. A handful of deliberately-extreme problems have compiler-dependent, coin-flip convergence (FORTRAN itself converges on them only at some optimization levels); those are listed, with the rationale, in examples/crosscheck_exclude.txt and skipped by the gate. USE_BLAS/USE_LAPACK builds (double or single precision only) take a genuinely different numeric path — a BLAS dnrm2 in enorm, BLAS dot/trsv/swap/rot in lmpar/qrsolv, and above all a LAPACK geqp3/geqrf QR factorization with different column pivoting (a different algorithm, not merely a last-bit FMA effect). They can flip a different, implementation-dependent set of problems, so they get an additional list, examples/crosscheck_exclude_blas.txt, applied only when those options are on.

Documentation

References

LMDER/LMDIF

HYBRJ/HYBRD

Simulating box constraints

Note that box constraints can easily be simulated in C++ Minpack, using a change of variables in the function (that hint was found in the lmfit documentation).

For example, say you want xmin[j] < x[j] < xmax[j], just apply the following change of variable at the beginning of fcn on the variables vector, and also on the computed solution after the optimization was performed:

  for (j = 0; j < 3; ++j) {
    real xmiddle = (xmin[j]+xmax[j])/2.;
    real xwidth = (xmax[j]-xmin[j])/2.;
    real th =  tanh((x[j]-xmiddle)/xwidth);
    x[j] = xmiddle + th * xwidth;
    jacfac[j] = 1. - th * th;
  }

This change of variables preserves the variables scaling, and is almost the identity near the middle of the interval.

Of course, if you use lmder, lmder1, hybrj or hybrj1, the Jacobian must be also consistent with that new function, so the column of the original Jacobian corresponding to x1 must be multiplied by the derivative of the change of variable, i.e jacfac[j].

Similarly, each element of the covariance matrix must be multiplied by jacfac[i]*jacfac[j].

For examples on how to implement this in practice, see the portions of code delimited by "#ifdef BOX_CONSTRAINTS" in the following source files: tlmderc.c, thybrj.c, tchkderc.c.

Equivalence table with other libraries

The following table may be useful if you need to switch to or from another library.

Equivalence table between MINPACK and NAG, NPL, SLATEC, levmar, GSL, SciPy and Eigen. A blank or * cell means that library has no direct equivalent; see the note on Ceres Solver below the table.
MINPACKNAGNPLSLATEClevmarGSLSciPyEigen
lmdifE04FCFLSQNDNDNLS1dlevmar_dif leastsqLevenbergMarquardt (num. diff)
lmdif1E04FYFLSNDN1DNLS1Edlevmar_dif leastsqLevenbergMarquardt (num. diff)
lmderE04GDF/E04GBFLSQFDNDNLS1dlevmar_dergsl_multifit_fdfsolver_lmsderleastsqLevenbergMarquardt
lmder1E04GZFLSFDN2DNLS1Edlevmar_der leastsqLevenbergMarquardt
lmstr**DNLS1  **
hybrdC05NCF*DNSQ gsl_multiroot_fsolver_hybridsfsolveHybridNonLinearSolver (num. diff)
hybrd1C05NBF*DNSQE  fsolveHybridNonLinearSolver (num. diff)
hybrjC05PCF*DNSQ gsl_multiroot_fdfsolver_hybridsjfsolveHybridNonLinearSolver
hybrj1C05PBF*DNSQE  fsolveHybridNonLinearSolver
covarE04YCF*DCOV gsl_multifit_covarleastsq (cov_x)*
chkderC05ZAF    check_grad*

Ceres Solver. Ceres Solver is not listed as a column because it does not map to MINPACK function-by-function: you build a ceres::Problem and call ceres::Solver. Its trust-region minimizer implements the same Levenberg-Marquardt method as lmder/lmdif (with a Dogleg alternative), selected via Solver::Options::trust_region_strategy_type = LEVENBERG_MARQUARDT (default) or DOGLEG (solving docs). Parameter covariance — the analogue of covar — is provided by ceres::Covariance (covariance docs), and the analogue of chkder is the Solver::Options::check_gradients option. Ceres solves non-linear least squares only, so there is no direct equivalent of the hybrd/hybrj non-linear-equation solvers (a system F(x)=0 would be handled by minimizing ||F(x)||2).

Other MINPACK implementations

History

Future work

There is now a very powerful alternative to MINPACK, which is the Ceres Solver. You may want to consider using Ceres for any new project.

The main feature that's missing on cminpack is the possibility to add constraints on variables. Simple boundary constraints should be enough, as implemented in ALGLIB or MPFIT, and they can easily be implemented using the hack above (section "Simulating box constraints").

levmar also has linear constraints, but they should not be necessary since linear constraints can be changed to box constraints by a simple change of variables. If you really need nonlinear constraints, and no reparameterization of variables (which may be able to linearize these constraints), you should consider using NLopt instead of cminpack.

Please file a GitHub issue for any suggestion or request.

Frédéric Devernay

Valid XHTML 1.1! Valid CSS! Level Double-A conformance icon, W3C-WAI Web Content Accessibility Guidelines 1.0