4

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?

Community
  • 1
user603
  • 21,225
  • 3
  • 71
  • 135

2 Answers2

2

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.

Dikran Marsupial
  • 46,962
  • 5
  • 121
  • 178
2

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$.

Jack Poulson
  • 136
  • 2
  • > would you know of any implementation of QR updates in R? [there are a couple implementations of chol-update). – user603 Dec 28 '11 at 17:59
  • 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. – Jack Poulson Dec 28 '11 at 18:03
  • Jack Poulson:> i moved it to stakeoverflow ( http://stackoverflow.com/questions/8659412/efficient-use-of-choleski-decomposition-in-r ). 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 – user603 Dec 28 '11 at 18:08
  • 1
    Answered there as well. – Jack Poulson Dec 28 '11 at 18:29
  • 1
    +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?). – Dikran Marsupial Jan 16 '12 at 17:16