This is somewhat related to this question.

I have two full rank matrices $A_1,A_2$ each of dimension $p \times p$ and a p vector $y$.

These matrices are closely related in the sense that matrix $A_2$ has $p-1$ row in common with matrix $A_1$.

Now, I'm interested in the vector $\beta_2$, where $\beta_2=(A_2'A_2)^{-1}(A_2'y)$.

My question: Is there is a fast way to get $\beta_2|(\beta_1,A_1,A_2,y)$ in R/matlab?


While I largely agree with what @DikranMarsupial said, it is typically a bad idea to form $A' A$ since it squares the condition number of $A$. It is more stable to recognize that changing the $i$th row of $A$ from $a'$ to $\hat a'$ is equivalent to perform the rank-one update $\hat A = A + e_i (\hat a-a)'$, where $e_i$ is a standard basis vector.

Thus, given a QR decomposition of $A$, the QR decomposition of $\hat A$ can be computed in $\mathcal{O}(n^2)$ work using, for example, MATLAB's qrupdate.

Roughly speaking, this approach should yield twice as many digits of accuracy as an approach which explicitly forms $A' A$.

  • $\begingroup$ > would you know of any implementation of QR updates in R? [there are a couple implementations of chol-update). $\endgroup$
    – user603
    Dec 28 '11 at 17:59
  • $\begingroup$ To be perfectly honest, I've never even used R, but this question was linked from your SciComp question. I was about to answer your question there, but you deleted it. $\endgroup$ Dec 28 '11 at 18:03
  • $\begingroup$ Jack Poulson:> i moved it to stakeoverflow ( stackoverflow.com/questions/8659412/… ). Im not sure it belongs to SciComp because im not sure the problem i have is a theoretical one (it could be due to me not using the right R functions). You are most welcome to answer it there $\endgroup$
    – user603
    Dec 28 '11 at 18:08
  • 1
    $\begingroup$ Answered there as well. $\endgroup$ Dec 28 '11 at 18:29
  • 1
    $\begingroup$ +1 for QR, I am to used to dealing with kernel methods, where you start with a kernel matrix, which is a covariance matrix in an implicit feature space. Sometimes I forget people make linear models as well ;o) P.S. ISTR that Morven Gentleman (?) published a nice paper on an incremental/decremental QR for linear regression, which would be ideal for this (title was something to do with Givens rotations without square roots?). $\endgroup$ Jan 16 '12 at 17:16

Changing a row of $A$ will give a rank-1 update of $A'A$, recomputing the inverse of a matrix following a rank-1 update can be achieved using the Sherman-Woodbury-Morrison formula. However this is numerically unstable, so it is better to perform the computation using a cholesky decomposition instead, which can be updated in a similar, but more numerically stable manner. IIRC MATLAB has some commands built in for this, but I can't remember the names ('lookfor cholesky' ought to find them). See also Matthias Seeger's software for Updating the Cholesky Decomposition for Low Rank Modifications. Matthias is a very bright guy, I suspect this is probably the best approach.

UPDATE - The command is cholupdate:

R1 = cholupdate(R,x) where R = chol(A) is the original Cholesky factorization of A, returns the upper triangular Cholesky factor of A + x*x', where x is a column vector of appropriate length. cholupdate uses only the diagonal and upper triangle of R. The lower triangle of R is ignored

so R1 = chol(A1'*A1); followed by R2=cholupdate(R1,x); ought to do the trick, where x is the differrence between the row of A1 that is replaced to get A2 and the new row in A2. I can't check as I don't have access to MATLAB at home, but it will be something of that nature. You can then get $\beta_2$ by the standard method for regression with Cholesy decomposition.

  • $\begingroup$ I've updated the post, I think it os what you want, but I can't check it until I go bak to work on the 3rd. $\endgroup$ Dec 26 '11 at 14:14
  • $\begingroup$ note update chol(A1'*A1) rather than just chol(A1). $\endgroup$ Dec 26 '11 at 15:09
  • $\begingroup$ I think you will find that this method yields a rank-2 update of A' A. $\endgroup$ Dec 28 '11 at 19:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.