Bounds Implementation

This section describes the implementation of Parameter bounds. The MINPACK-1 implementation used in scipy.optimize.leastsq for the Levenberg-Marquardt algorithm does not explicitly support bounds on parameters, and expects to be able to fully explore the available range of values for any Parameter. Simply placing hard constraints (that is, resetting the value when it exceeds the desired bounds) prevents the algorithm from determining the partial derivatives, and leads to unstable results.

Instead of placing such hard constraints, bounded parameters are mathematically transformed using the formulation devised (and documented) for MINUIT. This is implemented following (and borrowing heavily from) the leastsqbound from J. J. Helmus. Parameter values are mapped from internally used, freely variable values $P_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}$ to bounded parameters $P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}$. When both min and max bounds are specified, the mapping is:

$$\begin{array}{rcl}{P}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}& =& \arcsin (\frac{2(P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}-\mathrm{m}\mathrm{i}\mathrm{n})}{(\mathrm{m}\mathrm{a}\mathrm{x}-\mathrm{m}\mathrm{i}\mathrm{n})}-1)\\ P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}& =& \mathrm{m}\mathrm{i}\mathrm{n}+(\sin (P_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}})+1)\frac{(\mathrm{m}\mathrm{a}\mathrm{x}-\mathrm{m}\mathrm{i}\mathrm{n})}{2}\end{array}$$

With only an upper limit max supplied, but min left unbounded, the mapping is:

$$\begin{array}{rcl}{P}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}& =& \sqrt{(\mathrm{m}\mathrm{a}\mathrm{x}-P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}+1{)}^2-1}\\ P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}& =& \mathrm{m}\mathrm{a}\mathrm{x}+1-\sqrt{P_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}^2+1}\end{array}$$

With only a lower limit min supplied, but max left unbounded, the mapping is:

$$\begin{array}{rcl}{P}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}& =& \sqrt{(P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}-\mathrm{m}\mathrm{i}\mathrm{n}+1{)}^2-1}\\ P_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}& =& \mathrm{m}\mathrm{i}\mathrm{n}-1+\sqrt{P_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}^2+1}\end{array}$$

With these mappings, the value for the bounded Parameter cannot exceed the specified bounds, though the internally varied value can be freely varied.

It bears repeating that code from leastsqbound was adopted to implement the transformation described above. The challenging part (thanks again to Jonathan J. Helmus!) here is to re-transform the covariance matrix so that the uncertainties can be estimated for bounded Parameters. This is included by using the derivate $d{P}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{n}\mathrm{a}\mathrm{l}}/d{P}_{\mathrm{b}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{d}}$ from the equations above to re-scale the Jacobin matrix before constructing the covariance matrix from it. Tests show that this re-scaling of the covariance matrix works quite well, and that uncertainties estimated for bounded are quite reasonable. Of course, if the best fit value is very close to a boundary, the derivative estimated uncertainty and correlations for that parameter may not be reliable.

The MINUIT documentation recommends caution in using bounds. Setting bounds can certainly increase the number of function evaluations (and so computation time), and in some cases may cause some instabilities, as the range of acceptable parameter values is not fully explored. On the other hand, preliminary tests suggest that using max and min to set clearly outlandish bounds does not greatly affect performance or results.