scipy.optimize.leastsq with bound constraints

Question

I am looking for an optimisation routine within scipy/numpy which could solve a non-linear least-squares type problem (e.g., fitting a parametric function to a large dataset) but including bounds and constraints (e.g. minima and maxima for the parameters to be optimised). At the moment I am using the python version of mpfit (translated from idl...): this is clearly not optimal although it works very well.

An efficient routine in python/scipy/etc could be great to have ! Any input is very welcome here :-)

thanks!

Supportive

Answer 1

scipy.optimize.least_squares in scipy 0.17 (January 2016) handles bounds; use that, not this hack.

---

Bound constraints can easily be made quadratic, and minimized by leastsq along with the rest.
Say you want to minimize a sum of 10 squares Σ f_i(p)^2, so your func(p) is a 10-vector [f0(p) ... f9(p)],
and also want 0 <= p_i <= 1 for 3 parameters.
Consider the "tub function" max( - p, 0, p - 1 ), which is 0 inside 0 .. 1 and positive outside, like a \_____/ tub.
If we give leastsq the 13-long vector

[ f0(p), f1(p), ... f9(p), w*tub(p0), w*tub(p1), w*tub(p2) ]

with w = say 100, it will minimize the sum of squares of the lot: the tubs will constrain 0 <= p <= 1. General lo <= p <= hi is similar.
The following code is just a wrapper that runs leastsq with e.g. such a 13-long vector to minimize.

# leastsq_bounds.py
# see also test_leastsq_bounds.py on gist.github.com/denis-bz

from __future__ import division
import numpy as np
from scipy.optimize import leastsq

__version__ = "2015-01-10 jan  denis"  # orig 2012


#...............................................................................
def leastsq_bounds( func, x0, bounds, boundsweight=10, **kwargs ):
    """ leastsq with bound conatraints lo <= p <= hi
    run leastsq with additional constraints to minimize the sum of squares of
        [func(p) ...]
        + boundsweight * [max( lo_i - p_i, 0, p_i - hi_i ) ...]

    Parameters
    ----------
    func() : a list of function of parameters `p`, [err0 err1 ...]
    bounds : an n x 2 list or array `[[lo_0,hi_0], [lo_1, hi_1] ...]`.
        Use e.g. [0, inf]; do not use NaNs.
        A bound e.g. [2,2] pins that x_j == 2.
    boundsweight : weights the bounds constraints
    kwargs : keyword args passed on to leastsq

    Returns
    -------
    exactly as for leastsq,
http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

    Notes
    -----
    The bounds may not be met if boundsweight is too small;
    check that with e.g. check_bounds( p, bounds ) below.

    To access `x` in `func(p)`, `def func( p, x=xouter )`
    or make it global, or `self.x` in a class.

    There are quite a few methods for box constraints;
    you'll maybe sing a longer song ...
    Comments are welcome, test cases most welcome.

"""
    # Example: test_leastsq_bounds.py

    if bounds is not None  and  boundsweight > 0:
        check_bounds( x0, bounds )
        if "args" in kwargs:  # 8jan 2015
            args = kwargs["args"]
            del kwargs["args"]
        else:
            args = ()
#...............................................................................
        funcbox = lambda p: \
            np.hstack(( func( p, *args ),
                        _inbox( p, bounds, boundsweight ))) 
    else:
        funcbox = func
    return leastsq( funcbox, x0, **kwargs )


def _inbox( X, box, weight=1 ):
    """ -> [tub( Xj, loj, hij ) ... ]
        all 0  <=>  X in box, lo <= X <= hi
    """
    assert len(X) == len(box), \
        "len X %d != len box %d" % (len(X), len(box))
    return weight * np.array([
        np.fmax( lo - x, 0 ) + np.fmax( 0, x - hi )
            for x, (lo,hi) in zip( X, box )])

# def tub( x, lo, hi ):
#     """ \___/  down to lo, 0 lo .. hi, up from hi """
#     return np.fmax( lo - x, 0 ) + np.fmax( 0, x - hi )

#...............................................................................
def check_bounds( X, box ):
    """ print Xj not in box, loj <= Xj <= hij
        return nr not in
    """
    nX, nbox = len(X), len(box)
    assert nX == nbox, \
        "len X %d != len box %d" % (nX, nbox)
    nnotin = 0
    for j, x, (lo,hi) in zip( range(nX), X, box ):
        if not (lo <= x <= hi):
            print "check_bounds: x[%d] %g is not in box %g .. %g" % (j, x, lo, hi)
            nnotin += 1
    return nnotin

Motorboating

I've received this error when I've tried to implement it (python 2.7): File "[...]/leastsq_bounds.py", line 49, in leastsq_bounds return leastsq( funcbox, x0, **kwargs ) File "[...]/minpack.py", line 369, in leastsq shape, dtype = _check_func('leastsq', 'func', func, x0, args, n) File "[...]/minpack.py", line 20, in _check_func res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) TypeError: <lambda>() takes exactly 1 argument (5 given) – Headsman

@f_ficarola, sorry, args= was buggy; please cut/paste and try it again – Motorboating

Thank you for the quick reply, denis. However, in the meantime, I've found this: scipy.optimize.minimize(residualsModel, x0, args=(arg1, arg2, ...), method='SLSQP', bounds=[(xmin,xmax)]). Is there any difference with your code? – Headsman

Anyway, I now receive a different error: File "leastsq_bounds.py", line 58, in leastsq_bounds return leastsq( funcbox, x0, **kwargs ) File "minpack.py", line 369, in leastsq shape, dtype = _check_func('leastsq', 'func', func, x0, args, n) File "minpack.py", line 20, in _check_func res = atleast_1d(thefunc(*((x0[:numinputs],) + args))) File "leastsq_bounds.py", line 55, in <lambda> _inbox( p, bounds, boundsweight ))) File "leastsq_bounds.py", line 66, in _inbox "len X %d != len box %d" % (len(X), len(box)) AssertionError: len X 1 != len box 2 – Headsman

@f_ficarola, 1) SLSQP does bounds directly (box bounds, == <= too) but minimizes a scalar func(); leastsq minimizes a sum of squares, quite different. Which do you have, how many parameters and variables ? 2) what is np.shape( bounds ) ? it should be nparam x 2 -- run leastsq_bounds.check_bounds( x0, bounds) first – Motorboating

Yes, sorry... You're definitely right. I was passing my bounds in a wrong way (i.e., [xmin, xmax] instead of [[xmin, xmax]]). Thank you again, denis. – Headsman

Your January 2016 edits had already been stated in my answer, which you seem to have missed... – Romaic

Does using the tubs not give problems when considering your algorithm? If it uses gradients to find the minimum it might have trouble finding the right search area since the derivative remains mostly constant. – Silicate

@Jan, it gives scipy.optimize.leastsq a list of functions to minimize all together, yours and tubs; see the example with 13-long vector, also Penalty method. – Motorboating

Answer 2

scipy has several constrained optimization routines in scipy.optimize. The constrained least squares variant is scipy.optimize.fmin_slsqp

Stocky

Thanks! Will test this vs mpfit in the coming days for my problem and will report asap! – Supportive

Just tried slsqp. I may not be using it properly but basically it does not do much good. It does seem to crash when using too low epsilon values. And otherwise does not change anything (or almost) in my input parameters. I'll do some debugging, but looks like it is not that easy to use (so far). Will try further. – Supportive

In fact I just get the following error ==> Positive directional derivative for linesearch (Exit mode 8). This is why I am not getting anywhere.... not very useful. Any hint? – Supportive

Can you get it to work for a simple problem, say fitting y = mx + b + noise? scipy.org/Cookbook/FittingData is a fairly comprehensive tour of scipy.optimize and has examples. – Stocky

Answer 3

The capability of solving nonlinear least-squares problem with bounds, in an optimal way as mpfit does, has long been missing from Scipy.

This much-requested functionality was finally introduced in Scipy 0.17, with the new function scipy.optimize.least_squares.

This new function can use a proper trust region algorithm to deal with bound constraints, and makes optimal use of the sum-of-squares nature of the nonlinear function to optimize.

Notes:

The solution proposed by @denis has the major problem of introducing a discontinuous "tub function". This renders the scipy.optimize.leastsq optimization, designed for smooth functions, very inefficient, and possibly unstable, when the boundary is crossed.

The use of scipy.optimize.minimize with method='SLSQP' (as @f_ficarola suggested) or scipy.optimize.fmin_slsqp (as @matt suggested), have the major problem of not making use of the sum-of-square nature of the function to be minimized. These functions are both designed to minimize scalar functions (true also for fmin_slsqp, notwithstanding the misleading name). These approaches are less efficient and less accurate than a proper one can be.

Romaic

Answer 4

Have a look at: http://lmfit.github.io/lmfit-py/, it should solve your problem.

Claman

Thanks for the tip: one issue is that I would like to be able to have a self-consistent python module including the bounded non-lin least-sq part. This means either that the user will have to install lmfit too or that I include the entire package in my module. I will thus try fmin_slsqp first as this is an already integrated function in scipy. But lmfit seems to do exactly what I would need! – Supportive

Consider that you already rely on SciPy, which is not in the standard library. lmfit is on pypi and should be easy to install for most users. – Ersatz

link is broken! – Pincus