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

Brendan Simons brendansimons at yahoo.ca
Sun Mar 13 21:49:27 CST 2005

I am scratching my brain trying to find the quickest way to construct 
the covariance matrix corresponding to a least squares fit, and I've 
come across some documentation quirks that are confusing me.

The method scipy.optimize.leastsq returns, in its optional infodict 
parameter, a value called fjac which according to the docstring is:

	the orthogonal matrix, q, produced by the
	QR facotrization of the final approximate
           Jacobi matrix, stored column wise.

But scipy.optimize.leastsq calls the trusty MINPACK fortran routine 
lmder, who's comments say:

c       fjac is an output m by n array. the upper n by n submatrix
c         of fjac contains an upper triangular matrix r with
c         diagonal elements of nonincreasing magnitude such that
c               p^t *(jac^t *jac)*p = r^t *r,
c         where p is a permutation matrix and jac is the final
c         calculated jacobian. column j of p is column ipvt(j)
c         (see below) of the identity matrix. the lower trapezoidal
c         part of fjac contains information generated during
c         the computation of r.

And a recent comment to this list said:

On 10-Feb-05, at 11:23 PM, scipy-user-request at scipy.net wrote:

> R. Padraic Springuel wrote:
>> Is there a way to get the uncertainty in the resultant parameters out
>> of the optimize.leastsq fitting function?
> optimize.leastsq? reveals that it can return a dictionary of optional
> values, among which is the Jacobian matrix fjac at the final step.  The
> diagonal elements of fjac times its transpose can be taken as the
> squares of the errors in the corresponding coefficients.

So fjac is either:
1) the jacobian at the final step,
2) The Q of the QR factorization of the jacobian.
3) or the R of the QR factorization of the jacobian.

Why is this important?  Because the jacobian at the final step is 
needed to compute the covariance matrix, which gives the uncertainty in 
the fitted parameters, which is VERY important.

So, any MINPACK experts out there?  Can you explain to me what fjac is 
in clearer terms?  Is the scipy.optimize.leastsq docstring accurate?

PS, Here's a bit of linear algebra background:

lmder is an implementation of the Levenberg-Marquardt algorithm.    The 
levenberg marquardt method minimizes the sum of the squares of the 
residual functions by iteratively solving its derivative:

( J^t * J ) * x = J^t * f( x ),
for x:
x = ( J^t * J )^-1 * J^t * f( x )

where J is the Jacobian matrix (the derivatives of all the functions 
wrt each parameter x), x is the parameters to solve for, and f(x) is a 
vector of the residuals.   (There's another bit here, where the 
diagonals of ( J^t * J ) are multiplied by a damping factor, but that's 
not important for the present problem)

That first term on the right side of that last equation, ( J^t * J 
)^-1,  is important, because *that's* the covariance matrix of the 
parameters x.

Here's where I start to get confused.  I *think* lmder solves the above 
function with QR factorization, so that the above equation becomes:

x = ( R^t * Q^t * Q * R )^-1* R^t * Q^t * f( x )
   = (R^t * R)^-1* R^t * Q^t * f( x )
   = R^-1 * Q^T * f( x )

where Q, R are the results of the QR factorization of J.

If lmder is indeed doing this (I've looked at the code, but my Fortran 
skills are poor), than I should be able to get the value of R at the 
last step, and calculate my covariance from
(R^t * R)^-1

But at this point my head is hurting too much to know where to look for 
this matrix.   Am I way off base?

More information about the SciPy-user mailing list