[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
c p^t *(jac^t *jac)*p = r^t *r,
c
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?
-Brendan
-----------------
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