[SciPy-user] Re: optimize.leastsq -> fjac and the covariance redux

[Brendan Simons]

>>
> I made the change and now when using full_output for leastsq, the 
> covariance matrix is returned (as a separate output).
> Please test it and report back problems.
>
> -Travis
>

Hmmm, some issues.  I ran the following updated test script:

#----------begin test script--------------

from scipy import *
from scipy.optimize import leastsq

x = array([0., 1., 2., 3., 4., 5.], Float64)
coeffs =  [4., 3., 5., 2.]
yErrs = array([0.1, 0.2, -0.1, 0.05, 0, -.02], Float64)
m = len(x)
n = len(coeffs)

def evalY(p, x):
     """a simple cubic"""
     return x**3 * p[3] + x**2 * p[2] + x * p[1] + p[0]

def J(p, y, x):
     """the jacobian of a cubic (not actually
     a function of y, p since resid is linear in p)
     """
     result = zeros([m,n], Float64)
     result[:, 0] = ones(m, Float64)
     result[:, 1] = x
     result[:, 2] = x**2
     result[:, 3] = x**3
     return result

def resid(p, y, x):
     return y - evalY(p, x)

y_true = evalY(coeffs, x)
y_meas = y_true + yErrs

p0 = [3., 2., 4., 1.]

pMin, C, infodict, ier, mesg = leastsq(resid, p0,
	args=(y_meas, x),  full_output=True)

#check against know algorithms for computing C
JAtPMin = J(pMin, None, x)
Q_true, R_true = linalg.qr(JAtPMin)
C_true = linalg.inv(dot(transpose(JAtPMin), JAtPMin))
C_true2 = linalg.inv(dot(transpose(R_true), R_true))

print 'pMin'
print pMin
print '\nC'
print C
print '\nC_true'
print C_true
print '\nC_true2'
print C_true2

#-------end test script---------

Which is basically the same as before, but with a cubic instead of a 
quadratic.  Here's what I get:

#------ begin script output------

pMin
[ 4.1315873   2.95965608  4.99325397  2.00185185]

C
[[ 3.62323451 -1.22354445  0.21141966 -1.71957585]
  [-1.22354445  0.96031733 -0.04629627  0.43650769]
  [ 0.21141966 -0.04629627  0.01543209 -0.11574069]
  [-1.71957585  0.43650769 -0.11574069  0.89484084]]

C_true
[[ 0.96031746 -1.22354497  0.43650794 -0.0462963 ]
  [-1.22354497  3.62323633 -1.71957672  0.21141975]
  [ 0.43650794 -1.71957672  0.89484127 -0.11574074]
  [-0.0462963   0.21141975 -0.11574074  0.0154321 ]]

C_true2
[[ 0.96031746 -1.22354497  0.43650794 -0.0462963 ]
  [-1.22354497  3.62323633 -1.71957672  0.21141975]
  [ 0.43650794 -1.71957672  0.89484127 -0.11574074]
  [-0.0462963   0.21141975 -0.11574074  0.0154321 ]]

#----------end script output-----------

Here C has all the right outputs, but in the wrong order.  Evidently 
we're applying the
ipvt permutation matrix wrong somehow.

Here's something worse.  If I include the jacobian in the leastsq call 
as follows
pMin, C, infodict, ier, mesg = leastsq(resid, p0, args=(y_meas, x), 
Dfun = J, full_output=True)

the resulting pMin is the same as p0, and C is way off of C_true.  What 
gives?

-Brendan