Class-Based Levenberg-Marquardt Fitter #90

Closed
wants to merge 19 commits into from
@tecki

This is the same pull request as
#88
rebased onto the current master. The current
master already contains qr_multiply, and
this this pull request gets much more readable.

This pull request contains a new Class-Based
Levenberg-Marquardt fitter. It is the original minpack
lmder/lmdif algorithm translated to python and
strongly pythonised. It features:

  • a cool class-based interface
  • a function based interface if you don't like it so cool
  • the possibility to tell the algorithm that it went too far with the parameters while optimizing
  • the possibilty to stop fitting from another thread

The reason we used to have leastsq in FORTRAN
was that: first, it was already there, and 2nd, it is
fast. But this speed comes at a cost: when we
call python from FORTRAN, a lot of magic has to
be done. This magic prevents us, for example, to
properly pass exceptions through the FORTRAN
code.

The speed disadvantage from using Python is
marginal: since all numerical functions (this is
where most of the time is spent) are still in
FORTRAN, we don't loose much. I even re-coded
the last remaining inner loop which has no
simple FORTRAN equivalent in cython so we
also don't loose there. In the end, the Python code
could even be faster, as it uses LAPACK code,
which is typically pretty well optimized thanks to
ATLAS, while MINPACK still does everything by
hand, which might be sub-optimal (especially on
multi-core machines which are used very efficiently
by ATLAS).

Python also has the advantage of being much better
readable, and thus adopting for your needs is much
simpler. This is why I was able to add the "Parameter
ran away!" functionality.

tecki added some commits Oct 1, 2011
@tecki tecki ENH: add Levenberg-Marquardt fitting class
Add a new module lmfit which contains a class-based
Levenberg-Marquardt type fitter, and an additional
function interface to it.
25a5531
@tecki tecki ENH: add cython-compiled C-file and compile it
Added the result of compiling the pyx-file _qrsolv.pyx with
cython to C. Add the necessary lines into setup.py so that
it will get compiled to a module.
2211416
@tecki tecki ENH: introduce new classes and functions
add the new Fit class to the scipy.optimize namespace.
The function leastsq is not added as is is overwritten by
importing minpack later.
93823a6
@tecki tecki ENH: add tests for the new Fit ac9eb67
@tecki tecki TEST: demonstrate some details. DONT PULL THIS
This commit replaces scipy.optimize.leastsq from
minpack.leastsq to lmfit.leastsq. This is a demonstration that
minpack.leastsq and lmfit.leastsq are actually equivalent.
lmfit.leastsq passes the minpack tests, but only by re-writing
ValueErrors into TypeErrors, which I consider ugly.

This is just a nice demonstration and should not
be pulled into mainline scipy. Maybe in the far future,
when everything is really well tested.
fe4432f
@tecki tecki BUG: don't use C99 calls
The code used fsqrtf and fabsf for float variables.
We're using the double equivalents, hoping that the compiler
will optimize properly...
0884683
@tecki tecki DOC: document _qrsolv.pyx 017e28f
@tecki tecki ENH: add cythoned _qrsolv.c 997f39d
@pv pv commented on an outdated diff Oct 3, 2011
scipy/optimize/lmfit.py
+ This function is supposed to be a drop-in replacement for
+ `scipy.optimize.leastsq`. See the documentation there.
+ """
+ class MyFit(Fit):
+ def func(self, x):
+ return func(x, *args)
+
+ if Dfun is not None:
+ def jacobian(self, x, fvec, ret):
+ if col_deriv:
+ return Dfun(x)
+ else:
+ return Dfun(x).T
+ fit = MyFit()
+ try:
+ ret = fit.fit(x0, *kwargs)
@pv
SciPy member
pv added a note Oct 3, 2011

**kwargs

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on the diff Oct 3, 2011
scipy/optimize/lmfit.py
+ f = self.func(x)
+ if not hasattr(self, "eps"):
+ self.eps = finfo(f.dtype).eps
+ fdjac = self.jacobian(x, f, empty((len(f), len(x)), dtype=f.dtype))
+ r, jpvt = qr(fdjac, pivoting=True, mode="r")
+ condition = abs(r.diagonal()) <= r[0, 0] * eps
+ if condition.any():
+ nsing = condition.nonzero()[0][0]
+ else:
+ nsing = len(jpvt)
+ r = r[:nsing, :nsing]
+ ri = eye(r.shape[0], dtype=r.dtype)
+ ri = solve_triangular(r, ri.T, overwrite_b=True)
+ cov = dot(ri.T, ri)
+ ret = zeros_like(fdjac)
+ j1, j2 = meshgrid(jpvt[:nsing], jpvt[:nsing])
@pv
SciPy member
pv added a note Oct 3, 2011

Better to use np.ix_/np.ogrid instead of meshgrid.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on the diff Oct 3, 2011
scipy/optimize/lmfit.py
+
+ qtb : vector, length `N`
+ contains the first `N` elements of the vector ``Q.T b``.
+
+ Returns
+ -------
+ x : vector, length `N`
+ contains the least squares solution of the system
+ ``A x = b, D x = 0``.
+
+ s : matrix, `N` by `N`
+ contains the aforementioned matrix `S`. """
+ s = empty((r.shape[0], r.shape[1] + 1), dtype=r.dtype)
+ s[:, :-1] = r
+ s[:, -1] = qtb
+ # Eliminate the diagonal matrix D using a givens rotation.
@pv
SciPy member
pv added a note Oct 3, 2011

Coding style nitpick: empty line before comment. Use more empty lines to separate blocks of code that belong logically together.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv
SciPy member
pv commented Oct 3, 2011

Comments:

  • I don't think it is a good idea to go and replace the leastsq function with a new implementation. The new implementation is probably not bug-for-bug backwards compatible, and I believe our test suite is too poor to check for this. A better approach would be to add a new function, for instance named least_squares or lsq_fit or something else -- and deprecate the old function. This also gives more leeway on fixing any flaws with the existing interface.
  • The name Fit for the class is too general: one might also want to do fitting in L1-norm. LMFit or LSQFit would be better.
@rgommers
SciPy member

This seems to overlap with http://newville.github.com/lmfit-py/. I'm under the impression that Matt Newville would like that code to land in scipy at some point, which should be considered before merging this pull request.

@pv
SciPy member
pv commented Oct 3, 2011

The lmfit-py stuff is partly orthogonal to this, as it in the end calls Scipy's leastsq. However, there may be some overlap in the objective function specification.

@rgommers
SciPy member

True. It looked more similar to me after a quick glance than it is.

lmfit.py is probably not the best filename here though.

@rgommers
SciPy member

Could you mark the Cython-generated file as binary in .gitattributes? I got the spinning colorwheel for a minute after trying to look at the diff.

@rgommers rgommers commented on the diff Oct 30, 2011
scipy/optimize/lmfit.py
+ Set to None (the default) if no specific parameter is a problem.
+ """
+ def __init__(self, number=None):
+ Exception.__init__(self)
+ self.number = number
+
+class Fit(object):
+ """ Least-Squares Fitting using Levenberg-Marquardt
+
+ This class contains methods to minimize the sum of the squares
+ of `M` nonlinear functions with `N` parameters by a modification of the
+ Levenberg-Marquardt algorithm. The user must overwrite the method
+ `func` that calculates the functions, and, optionally, the method
+ `jacobian` to calculate its partial derivatives.
+
+ Parameters
@rgommers
SciPy member

These are Returns, not Parameters.

@dlax
SciPy member
dlax added a note Mar 17, 2012

In fact, it seems these are Attributes.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers
SciPy member

It would be good if the parameters of the fit method were named the same as those in the proposed unified optimization interface where possible. See minimize() in #94 and https://github.com/dlaxalde/scipy/wiki/Improvements-to-optimize-package

@rgommers
SciPy member

PR 94 is now merged by the way, so that interface probably won't change anymore.

tecki added some commits Mar 16, 2012
@tecki tecki ENH add minimize function
add a function with a call syntax like optimize.minimize
4302a70
@tecki tecki ENH add more tests
add tests to test if lm fitting works with explicit jacobians
and the new function minimize.
d24d835
@tecki tecki ENH add convenience classes
Add a convenience class to easily turn a python function
into a fittable function, and a decorator that does the same.
2fb5190
@tecki tecki ENH add tests for convenience classes e93c97b
@dlax dlax commented on the diff Mar 17, 2012
scipy/optimize/lmfit.py
+ if error.number is not None:
+ self.scale[error.number] *= 2
+ else:
+ delta = 0.5 * min(delta, 10 * pnorm)
+ par *= 2
+ if self.running:
+ continue
+ else:
+ return x
+ if self.iterations == 0:
+ delta = min(delta, pnorm)
+ testfnorm = norm(testf)
+ actual_reduction = -1
+ if 0.1 * testfnorm < fnorm:
+ actual_reduction = 1 - (testfnorm / fnorm) ** 2
+ temp1 = norm(dot(r, p[ipvt])) / fnorm
@dlax
SciPy member
dlax added a note Mar 17, 2012

if fnorm == 0.0, you get a ZeroDivisionError

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
tecki added some commits Mar 17, 2012
@tecki tecki FIX dont crash if solution is exact
if the solution found by fitting was exact,
a division by zero occured. Fixed that.
6feafa8
@tecki tecki DOC updated documentation ede6ba5
@tecki tecki FIX fixed typo
thanks to pv
b665d16
@tecki tecki Revert "TEST: demonstrate some details. DONT PULL THIS"
This reverts commit fe4432f.

This patch had replaced the minpack leastsq with the lmfit
leastsq. This was intended for testing purposes only, so
that one can proove that indeed both functions work the same
and give the same result. As people started complaining
I wanted to kill good ole minpack leastsq (which I didn't
intend) I hereby revert the commit.
403527f
@tecki tecki FIX typo
increase number of function evals, don't slash it.
5988882
@tecki tecki ENH improve readability of code
make code nicer by introducing more comments and
grouping of lines.
bf8869f
@tecki tecki FIX mark cython generated code binary 1141a6d
@tecki

I was thinking a lot about pv's comment on naming the class: on one hand,
LMFit sounds good, but what about my other classes? LMFunction? lmfitfunction?
that sounds awkward. Also, this is nothing else than prefixing, which is
considered not so nice in python (see PEP 8). Maybe just let the function names
as is, but don't import them into the optimize namespace, but only let them be
used as lmfit.Fit, lmfit.Function etc.?

@pv
SciPy member
pv commented Mar 17, 2012

What would be useful here is considering this work into two independent parts: first, the new implementation of Levenberg-Marquardt, with some new features (including the addition of the InvalidParameters communication mechanism). And second, a new Fit/Function/fitfunction user interface for specifying fitting problems. The second part is pretty similar in purpose to curve_fit.

I think these two parts should not be tightly coupled to each other. This is connected to the class naming: names like Fit or fitfunction are OK for the user interface, but not OK for the algorithm itself. If you look at the lmfit-py package, it allows for anneal or bfgs to be used in the fitting (sure, these are probably often less appropriate choices than LM). The user should be able to swap in another fitting algorithm (say, 1-norm fitting) easily, without having to change the way the problem is specified, which is why I'm arguing for loose coupling.

I don't think adding sub-namespaces to scipy.optimize would be good.

The way out of this, as I see it:

  1. Rename the class Fit that has the implementation to _LM (not available in the public API), and provide a function-based interface leastsq_lm. (Or, split our current leastsq function into an internal function with MINPACK implementation, and allow this implementation as the second option.) This part is "uncontroversial" and can then be just merged in without any further thought.

  2. The InvalidParameter could be considered as a common interface element in all of the algorithms (although I guess not many of them support that).

  3. Move the class-based fitting interfaces to a new file _fitting.py. Provide a class named Fit there, which can directly make use of the functionality in _LM, but can also be extended to provide for other fitting algorithms. How to exactly do this requires some thought. method keyword argument is one option. Or, one could offer one base class per algorithm. The latter option requires choosing names for each of the algorithms provided.

The motivation here is that the interfaces scipy.optimize provides for specifying problems should be the same, independent of what algorithm is used. Previously, we only had the "raw" function-based one and curve_fit. I believe it is best to think carefully about how the new interfaces should look like, rather than providing a smörgåsbrod of interfaces (a different one for each algorithm).

@dlax
SciPy member
@argriffing

What's going on in this pull request?

@pv
SciPy member
pv commented May 12, 2013

It's stalled: the algorithmic part is OK, the new interfaces proposed controversial. The algorithmic part is a translation of the MINPACK Fortran code currently used for leastsq, why merging this has not been high on the priority list.

However, this could perhaps be extended to Levenberg-Marquardt supporting sparse Jacobians, which is on the wishlist. Another retranslation (though without Cython) is called MPFIT, but it has some added features.

@Tillsten

One possible advantage of this is that the minpack code uses it's own linalg-rountines in some places. Using LAPACK instead could give a speed boast as well. Also i one could make the method to solve the linear system a parameter, e.g. SVD, QR, LU and Cholsky.

@pv pv added the PR label Feb 19, 2014
@pv pv removed the PR label Aug 13, 2014
@sturlamolden

MINPACK (used for leastsq) uses an unoptimized QR factorization derived from LINPACK. It is a major botteneck. I would not be surprised if a plain Python implementation that calls LAPACK is faster. We could also reuse the recent Cython code for updating QR factorizations. Another problem is that the MINPACK Fortran 77 code is not reentrant. The advantage of the MINPACK version is that is well-tried and known to be very numerically stabile.

SVD does not make a lot of sense as the diagonal bias added to the covariance matrix ensures that it is not ill-conditioned. If it is, Levenberg-Marquardt will detect it (residual sum of squares increases due to rounding error) and more bias is added to the diagonal.

@rgommers
SciPy member
@sturlamolden

No, I was not, but that is good to hear :)

@dashesy

Great to hear about the GSoC project! I think many people will be happy. Would be interesting to compare this with lmfit (which works great) in terms of performance and accuracy of fit. Also the api of lmfit is very intuitive, so another point of comparison would be the api design.

@sturlamolden

lmfit is just a wrapper for the optimizers already in SciPy. Wrapping scipy.optimize.leastsq with a different API and writing a new Levenberg-Marquardt algorithm from scratch with NumPy and LAPACK is not in the same league. lmfit cannot do any better than what's currently in SciPy.

@dashesy

@sturlamolden It is not too trivial, working around the bounds requires some transformation before and after the fit, it is not a simple bound check because that would throw off the gradient. Also mathematical expression can be used using ast for constraints (that is non-trivial). Even if it was just a nice api It would be more useful if the new Levenberg-Marquardt algorithm could help lmfit to be simpler (native support for bounds) and faster. Using LAPACK will definitely help because of the matrix decompositions for example. In addition the api of lmfit is nice and mature now, so I thought a comparison (when GSoC project is done) would be nice.

@charris
SciPy member

It would be helpful if all those interested could review and comment on Nikolay's work as it progresses. You can follow some of his plans at https://nmayorov.wordpress.com/

@ev-br
SciPy member

OK, the main part of the GSoC work has been merged, gh-5044.

@nmayorov

Perhaps it's time to close this PR :)? By the way it is the oldest one currently.

@rgommers
SciPy member

Perhaps it's time to close this PR :)? By the way it is the oldest one currently.

@nmayorov I agree. The algorithmic part here is still potentially useful, but by closing this PR that doesn't get lost. So if anyone wants to use that code still, they can still grab it. Hence closing.

@rgommers rgommers closed this Jan 18, 2016
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment