[SciPy-user] singular matrix linalg.basic.LinAlgError in optimize.leastsq

[Taro Sato]

I frequently use optimize.leastsq, and there is this annoying tendency
for the routine to throw an exception when the resulting covariant
matrix ends up singular.  The error message follows.

Basically, I need full_output=1 to just obtain a stopping condition in
the output parameter ier, but the way the code prepare the full
output, it also tries to compute covariant matrix.  The problem is
when it fails, it throws an exception which is not handled within
minpack.py:

-----------------------------------------------------------------
Traceback (most recent call last):
  File "./test.py", line 27, in ?
    if __name__ == '__main__': main()
  File "./test.py", line 23, in main
    p = O.leastsq(residual, [10.,.1,.1], args=(x,y),full_output=1)
  File "/usr/lib/python2.3/site-packages/scipy/optimize/minpack.py",
line 271, in leastsq
    cov_x = sl.inv(dot(transpose(R),R))
  File "/usr/lib/python2.3/site-packages/scipy/linalg/basic.py", line
221, in inv
    if info>0: raise LinAlgError, "singular matrix"
scipy.linalg.basic.LinAlgError: singular matrix
----------------------------------------------------------------

I'm not sure when exactly a covariant matrix cannot be computed this
way.  My test code (attached below and used to produce the error
above) seems to indicate when the initial guess parameters are way off
from the best parameters, the above error can result.

The existing behavior of leastsq is very annoying, since to obtain a
stopping condition that can be handled easily (i.e., in ier as an
integer), I must use full_output=1, but when the computation of a
covarant matrix is not well behaving, it simply explodes; there's no
way to get all the other info that are relevant for knowing why the
fitting failed.

Is there any better way to handle the situation (which might include
suggesting a modification to minpack.py)?  Also, it appears to make
more sense (at least to me) for leastsq to return ier when
full_output=0, rather than mesg, as it makes error handling easier...a
string message is human friendly but not really friendly to coders....
:)

By the way, in the mailing list archive, I noticed there was a similar
request about a year ago for adding an exception handling in this
case.

Thank you for your time,
Taro


---CODE BEGINS (test.py)-------------------------------------------

#!/usr/bin/env python
import numpy as N
import scipy.optimize as O



def fgauss(lamb, params):
    lambc, a, b = params
    return 1.+a*N.exp(-0.5*((lamb-lambc)/b)**2)


def residual(p, x, y):
    return y - fgauss(x, p)


def main():
    p0 = [100., 50., 5.]
    x = N.arange(200).astype(float)
    y0 = fgauss(x, p0)
    y = N.zeros(y0.shape).astype(float)
    for i in xrange(y.size): y[i] = N.random.poisson(y0[i])

    p = O.leastsq(residual, [10.,.1,.1], args=(x,y),full_output=1)
    print p[0],p[3:]


if __name__ == '__main__': main()


---CODE ENDS-------------------------------------------