10

I was wondering what the difference between the two methods scipy.optimize.leastsq and scipy.optimize.least_squares is?

When I implement them they yield minimal differences in chi^2:

>>> solution0 = ((p0.fun).reshape(100,100))
>>> # p0.fun are the residuals of my fit function np.ravel'ed as returned by least_squares
>>> print(np.sum(np.square(solution0))) 
0.542899505806

>>> solution1 = np.square((median-solution1))
>>> # solution1 is the solution found by least_sq, it does not yield residuals thus I have to subtract it from the median to get the residuals (my special case)
>>> print(np.sum((solution1)))
0.54402852325

Could anybody expand on that or point out where I can find an alternative documentation, the one from scipy is a bit cryptic.

MartyMacGyver
  • 9,483
  • 11
  • 47
  • 67
Sebastiano1991
  • 867
  • 1
  • 10
  • 26

3 Answers3

5

From the docs for least_squares, it would appear that leastsq is an older wrapper.

See also

leastsq   A legacy wrapper for the MINPACK implementation of the Levenberg-Marquadt algorithm.

So you should just use least_squares. It appears that least_squares has additional functionality. Foremost among them is that the default "method" (i.e. algorithm) used is different:

  • trf : Trust Region Reflective algorithm, particularly suitable for large sparse problems with bounds. Generally robust method.
  • dogbox : dogleg algorithm with rectangular trust regions, typical use case is small problems with bounds. Not recommended for problems with rank-deficient Jacobian.
  • lm : Levenberg-Marquardt algorithm as implemented in MINPACK. Doesn’t handle bounds and sparse Jacobians. Usually the most efficient method for small unconstrained problems.

Default is trf. See Notes for more information.

The old leastsq algorithm was only a wrapper for the lm method, which—as the docs say—is good only for small unconstrained problems.

The difference you see in your results might be due to the difference in the algorithms being employed.

Praveen
  • 6,872
  • 3
  • 43
  • 62
4

The key reason for writing the new Scipy function least_squares is to allow for upper and lower bounds on the variables (also called "box constraints"). This was a highly requested feature.

This apparently simple addition is actually far from trivial and required completely new algorithms, specifically the dogleg (method="dogleg" in least_squares) and the trust-region reflective (method="trf"), which allow for a robust and efficient treatment of box constraints (details on the algorithms are given in the references to the relevant Scipy documentation ).

Also important is the support for large-scale problems and sparse Jacobians.

When bounds on the variables are not needed, and the problem is not very large, the algorithms in the new Scipy function least_squares have little, if any, advantage with respect to the Levenberg-Marquardt MINPACK implementation used in the old leastsq one.

However, the very same MINPACK Fortran code is called both by the old leastsq and by the new least_squares with the option method="lm". For this reason, the old leastsq is now obsoleted and is not recommended for new code.

divenex
  • 15,176
  • 9
  • 55
  • 55
  • Is it possible to provide different bounds on the variables. My problem requires the first half of the variables to be positive and the second half to be in [0,1]. It is hard to make this fix? – linello Feb 08 '19 at 16:47
  • 1
    Of course, every variable has its own bound: `bounds=[(L1, L2,...), (U1, U2,...)]` – divenex Feb 08 '19 at 20:56
2

In least_squares you can give upper and lower boundaries for each variable

There are some more features that leastsq does not provide if you compare the docstrings

dnalow
  • 974
  • 4
  • 14