What are recommended ways of doing nonlinear least squares, min $\sum err_i(p)^2$, with box constraints $lo_j <= p_j <= hi_j$ ? It seems to me (fools rush in) that one could make the box constraints quadratic, and minimize $$ \sum_i err_i(p)^2 + C * \sum_j tub( p_j, lo_j, hi_j )^2 $$ where $tub( x, lo, hi )$ is the "tub function" shaped like \___/, $max( lo - x, 0, x - hi )$.
Does this work in theory, work in practice ?
(There seem to be many many theoretical papers on NLS+, but my interest is practical —
real or realistic test cases would help me to choose among methods.)

(Experts, please add tags: "least-squares" ?)

  • 5
    $\begingroup$ Replacing strict constraints with penalty functions is a common technique in numerical optimization. It seems that what you are proposing is a particular form of that replacement. You can read all about similar techniques, e.g., here: stanford.edu/~boyd/cvxbook $\endgroup$ Commented Mar 24, 2012 at 18:01
  • $\begingroup$ You could use a suitable parametrization of $p$ to satisfy the box constraints (for example $p_i=\min(\max(lo_j, p'_j),hi_j)$. With respect to NLS solvers, Levenberg-Marquardt is good enough most of the time, perhaps combined with a global stochastic optimizer like simulated annealing. Some commercial toolboxes also seem to offer trust region methods based on adaptive response surface models, which looks like a reasonable generalization of Levenberg-Marquardt to me. $\endgroup$ Commented Mar 24, 2012 at 23:55

4 Answers 4


Adding squared penalty terms to get rid of constraints is a simple approach giving an accuracy of order 1/penalty factor only. Hence it is not recommended for high accuracy unless you let the penalty go to infinity during the computation. But a high penalty factor makes the Hessian very ill-conditioned, which limits the total accuracy achievable without taking into account the constraints explicitly.

Note that bound constraints are much easier to handle than general constraints, whence they are virtually never converted to penalties.

The solver L-BFGS-B (used with an about 5-dimensional history) usually solves bound constrained problems very reliably and fast in both low high dimensions. Exceptions are misconvergence on problems that can become very flat far off the solutions, where it is easy to get stuck with a descent method.

We made lots of experiments on very diverse functions in many different dimensions, with many different solvers available, as we needed a very robust bound-constrained solver as part of our global optimization software. L-BFGS-B clearly stands out as general purpose method, though of course on sme problems other solvers perform significantly better. Thus I'd recommend L-BFGS-B as a first choice, and would try alternative techniques just in case L-BFGS-B handles your particular class of problems poorly.

  • $\begingroup$ L-BFGS is available in IPOPT, I revised my answer. $\endgroup$
    – Ali
    Commented May 21, 2012 at 11:23

I would simply use the general-purpose NLP solver IPOPT. It is the most robust solver among those I have tried.

Unless you have some very special requirements, there is no reason why you should insist on a problem specific solver that only works for NLS with box-constraints.

A change in requirements (e.g. adding nonlinear constraints) would cause a major headache with a problem specific solver. You will have no such problems if you use the general-purpose IPOPT.

UPDATE: You can try L-BFGS with IPOPT, see under Quasi-Newton in the documentation.

The solution procedure can become faster at the expense of spoiling the remarkable robustness of IPOPT. In my opinion, use the exact derivatives if they are available. I would start messing with approximations (such as L-BFGS) only if I had proven performance issues.

  • $\begingroup$ I don't know how well IPOPT works, but your suggestion reminds me on similar statements by advocates of the downhill simplex method. Because nonlinear least squares is a common problem class, outright rejecting using one of the existing NLS solvers seems a bit suspicious to me. $\endgroup$ Commented Mar 24, 2012 at 23:37
  • $\begingroup$ @ThomasKlimpel Well, denis should give us more details, then we could help him choosing the right solver. :) Or he can check it for himself and find out which one fits his needs best. IPOPT seems to be a good solver to start with. $\endgroup$
    – Ali
    Commented Mar 25, 2012 at 7:59
  • $\begingroup$ @Ali, can you point to some "real or realistic test cases" please ? $\endgroup$
    – denis
    Commented Mar 25, 2012 at 15:56
  • $\begingroup$ @denis I could but I have no intention to do so, it would throw you off the track. The only thing that matters is how IPOPT handles your problem. Unless you have some very special requirements, it should solve it nicely. IPOPT has interfaces to MATLAB, C++, C, Fortran, R, AMPL, CUTEr. Pick one interface and test what happens with your problem :) Testing a problem specific solver would not be easier either. $\endgroup$
    – Ali
    Commented Mar 25, 2012 at 18:16
  • $\begingroup$ @Thomas Klimpel, guess I wasn't clear: I'm not rejecting, not asking about packages, but asking for insights or test cases: why might this trivial method not work well ? $\endgroup$
    – denis
    Commented Mar 26, 2012 at 9:35

(Years later) two solvers that handle box constraints:

  • Scipy least_squares has 3 methods, with extensive doc:

    1. 'trf’: Trust Region Reflective
    2. 'dogbox'
    3. 'lm': a legacy wrapper for MINPACK, without box constraints.
  • ceres
  • 1
    $\begingroup$ The Scipy one explicitly says that the Levenberg-Marquardt algorithm can't handle box constraints. $\endgroup$
    – tholy
    Commented Feb 16, 2018 at 22:26

The R minpack.lm CRAN package provides a Levenberg-Marquardt implementation with box constraints.

In general, Levenberg-Marquardt is much better suited than L-BFGS-B for least-squares problems. It will converge (much) better on challenging problems. It will also be much faster than the general purpose IPOPT, as it is tailored to non-linear least-squares problems.

The R package chooses a very straightforward projection approach to enforce the constraints (see the source code). Depending on the LM implementation you are using, it could be simple to include.

Now, the suggestion in the comments of using a transformation, (for example a sine transformation as in scipy) is also a good, simple alternative to transform your unconstrained LM algorithm into a constrained one. You will also need to include the transformation in the Jacobian if the Jacobian is analytic.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.