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

[Brendan Simons]

>Brendan,
>
>Thanks for your detailed comments.  I would love to 
>see the covariance 
>matrix as a standard optional output.  If you can
help 
>us get there from 
>the code that would be great. 
>
>As I look at the Fortran docstring, I think it's 
>pretty clear that the 
>leastsq Python docstring is wrong (sorry about that 
>since I probably 
>wrote it...)  fjac is basically the R matrix (once
you 
>apply some 
>permutations), thus, I'm pretty sure you can use the 
>upper triangular 
>part of the fjac matrix with the ipvt vector (also in

>infodict) to 
>caculate the covariance matrix as you describe.
>
>So, I would construct a permutation matrix P (using 
>ipvt) like this:
>
>n = len(ipvt)
>p = mat(take(eye(n),ipvt))
>
>and then get R using some thing like
>
>r = mat(MLab.triu(fjac))
>R = r * p.T
>
>Then covariance matrix estimate is
>
>C = (R.T * R).I
>
>
>This is totally untested, but I think you are on the 
>right track to use 
>fjac as the R matrix of the QR factorization.
>
>Thanks for your help,
>
>-Travis

Travis,
I too would like to see the covariance in the leastsq
optional outputs.  Using your suggestion, I whipped up
the following script, which fits a simple quadratic
(chosen because I can calculate the jacobian easily by
hand):

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

"""script to test the fjac output matrix of
scipy.optimize.leastsq
Placed in public domain by Brendan Simons
March 2005
"""

from scipy import *
import MLab
from scipy.optimize import leastsq

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

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

def J(p, x):
    """the jacobian of a quadratic"""
    result = zeros([m,n], typecode='f')
    result[:, 0] = ones(m, typecode='f')
    result[:, 1] = x
    result[:, 2] = x*x
    return result

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

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

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

pMin, infodict, ier, mesg = leastsq(resid, p0, 
    args=(y_meas, x), full_output=True)
fjac = infodict['fjac']
ipvt = infodict['ipvt']

#here's travis' suggestion
perm = mat(take(eye(n),ipvt-1))
r = mat(MLab.triu(fjac[:,:n]))
R = r * perm.T
C = (R.T * R).I

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

print 'pMin'
print pMin
print '\nfjac'
print fjac
print'\nipvt'
print ipvt
print '\nperm'
print perm
print '\nr'
print r
print '\nR'
print R
print '\nR_true'
print R_true
print '\nC'
print C
print '\nC_true'
print C_true
print '\nC_true2'
print C_true2

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

Note that I had to make the following modifications to
your algorithm:

a) The Fortran docstring says r is in the first nxn
submatrix of fjac, so I truncated fjac with [:, :n]
before applying MLab.triu

b) the identity matrix columns returned by ipvt are
1-based, not 0-based (Fortran is 1-based I guess).  So
I subtracted 1 in 
which results in perm = mat(take(eye(n),ipvt-1))

The results from the script are as follows:

#--- begin script output----
pMin
[ 4.13714286  2.93428571  5.00714286]

fjac
[[-31.28897539  -0.03196014  -0.12784056  -0.28764125 
-0.51136222  -0.79900346]
 [ -7.19103119   1.81357942   0.59589038   0.51365984 
 0.17797854  -0.41115309]
 [ -1.75780755   1.30104626  -1.10335453   0.03523648 
-0.29732131  -0.95308465]]

ipvt
[3 2 1]

perm
Matrix([[ 0.,  0.,  1.],
       [ 0.,  1.,  0.],
       [ 1.,  0.,  0.]])

r
Matrix([[-31.28897539,  -0.03196014,  -0.12784056],
       [  0.        ,   1.81357942,   0.59589038],
       [  0.        ,   0.        ,  -1.10335453]])

R
Matrix([[ -0.12784056,  -0.03196014, -31.28897539],
       [  0.59589038,   1.81357942,   0.        ],
       [ -1.10335453,   0.        ,   0.        ]])

R_true
[[ -2.44948983  -6.12372446 -22.45365715]
 [  0.           4.18329954  20.916502  ]
 [  0.           0.           6.11010218]
 [  0.           0.           0.        ]
 [  0.           0.           0.        ]
 [  0.           0.           0.        ]]

C
Matrix([[ 8.21428627e-001, -2.69897976e-001,
-3.08050728e-003],
       [-2.69897976e-001,  3.92718047e-001, 
7.01607645e-004],
       [-3.08050728e-003,  7.01607645e-004, 
1.03332016e-003]])

C_true
[[ 0.8214277  -0.58928472  0.08928555]
 [-0.58928496  0.72678488 -0.13392843]
 [ 0.08928559 -0.13392843  0.02678569]]

C_true2
[[ 0.82142925 -0.58928621  0.08928578]
 [-0.58928627  0.72678584 -0.13392855]
 [ 0.08928579 -0.13392855  0.0267857 ]]


#--- end script output----

So you can see it doesn't work as expected. C does not
equal C_true.  More troublesome is that *no* columns
of fjac resemble R_true, no matter what permutation
matrix is applied.  At this point I'm lost as to what
is what.  Any suggestions?

-Brendan
--
PS, If and when we figure out the answer, I'd like to
submit the change to CVS.  However, never having done
such a thing, I will need some instruction.


______________________________________________________________________ 
Post your free ad now! http://personals.yahoo.ca