[SciPy-dev] code for incremental least squares

josef.pktd@gmai... josef.pktd@gmai...
Wed Feb 17 13:01:37 CST 2010

On Wed, Feb 17, 2010 at 12:35 PM, Nathaniel Smith <njs@pobox.com> wrote:
> On Wed, Feb 17, 2010 at 7:10 AM,  <josef.pktd@gmail.com> wrote:
>> In the description of some estimators and the optimization routines
>> used, I have read that the updating is done on the cholesky or QR
>> decomposition of the inverse matrix ( (x'x)^{-1} or inverse Hessian).
>> Except for simple least squares problem the inverse is often more
>> important than the the original X'X.
>> Are there algorithms available for this, and is there an advantage for
>> either way? From what I have seen in incremental_ls, the updating
>> works on QR of x'x and then the inverse of the QR decomposition is
>> taken. Assuming a well behaved non-singular X'X.

Maybe, I'm getting ahead of myself and mixing oranges and apples.

One of the estimators that I have been looking at recently is a F-GLS
with a general covariance matrix
as in mixed/random effects models or panel data. I was trying to
figure out how to work more efficiently with the matrix decompositions
to calculate the covariance matrix of the parameter estimates.

> I'm not sure what you mean here -- when using QR for OLS, you take the
> QR decomposition of X, not X'X (this is what gives the improved
> numerical stability -- if X is ill-conditioned, then X'X is
> ill-conditioned squared, so you want to avoid it).

I didn't realize QR is on X and not on X'X. I think, I was looking
more at XtXIncrementalLS than at QRIncrementalLS. Because of the block
structure of the heteroscedasticity, you update separately by groups,
if my reading is correct.
I assume now XtXIncrementalLS is your second implementation that
doesn't use QR. I think, I understand now the basic structure of the
QR method for OLS.

>  http://en.wikipedia.org/wiki/Linear_least_squares#Orthogonal_decomposition_methods
> You do need to explicitly calculate inv(X'X) at the end in order to
> get standard error estimates, and it's nice if you can re-use any
> decompositions you have lying around to make this faster (for the QR
> case, inv(X'X) = inv(R'Q'QR) = inv(R'R), which can take advantage of
> R's triangular structure).

I saw this, and it reminded me of the possibility to update the
decomposition of the inverse directly. This part, I really have to
look at later when I work on GLS again.

> But I'm not sure what you're referring to when you say "some
> estimators", or why those routines want to incrementally compute the
> inverse Hessian, so I can't help you there.

most likely apples and oranges again. I was also working on time
series models where the optimization problem is non-linear, but the
underlying structure is linear, eg. arma with maximum likelihood

Since I have only a vague idea what the low level linear algebra does,
I better stop until I find time to figure this out. I just learn it as
the need or the opportunity arises. (Similar to my fight with fft a
while ago.)



> -- Nathaniel
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev

More information about the SciPy-Dev mailing list