3
$\begingroup$

Given the matrix $X_{n,p}$ from the least squares problem $$ \mathbf{X} \cdot \mathbf{\beta} = z $$

Where the normal equation is:

$$ \mathbf{\hat{\beta}} = \left(\mathbf{X}^T \mathbf{X}\right)^{-1} \mathbf{X}^T z $$

I was very happy when I found the existence of the Woodbury matrix identity unfortunantly I am struggling to use it (don't know if it's possible) for my problem.
$$ {(A+UCV)}^{-1}=A^{-1}-A^{-1}U{(C^{-1}+VA^{-1}U)}^{-1}VA^{-1} $$

The Problem

I want to compute a new $(X^TX)^{-1}$ after removing the first $k$ rows of $X$. I heard maybe it's called the leave-one-out (k-out?) statistics.

I found that for the my case the Woodbury formula is something like: $$ {((X^TX)+UCV)}^{-1}=(X^TX)^{-1}-(X^TX)^{-1}U{(C^{-1}+V(X^TX)^{-1}U)}^{-1}V(X^TX)^{-1} $$ where $+UCV$ should somehow subtract the first $k$ rows.

If someone can give some help or point to some direction or references.

New contributor
eusoubrasileiro is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
3
$\begingroup$

You've basically laid out the key facts, I think you just need a hint on how to fit them all together. Here's a quick-and-dirty overview.

I think it's easier to see how to accomplish your goal if you build up from the Sherman-Morrison formula, which is just a special case of the Woodbury matrix identity. The Sherman-Morrison formula is a rank-1 update, while the Woodbury identity is a rank-$r$ update.

We have a matrix $X_{n \times p}$ with $n$ samples/observations of $p$ variables/features and $X$ is full rank. The product $X^\top X$ can be viewed as a sum of outer products. Denote $x_j$ the $j$th column of $X^\top$ (i.e. the transpose of the $j$th row of $X$). Suppose we leave out one row $k$. We have

$$ \begin{align} X^\top X &= \sum_j x_j x_j^\top \\ &= x_k x_k^\top + \sum_{j\neq k} x_j x_j^\top \\ X^\top X - x_k x_k^\top &= \sum_{j\neq k} x_j x_j^\top. \end{align} $$

Relating this to the Sherman-Morrison formula can be done by inspection. Sherman-Morrison gives us $$ (A + uv^\top)^{-1} = A^{-1} - \frac{A^{-1}uv^\top A^{-1}}{1+v^\top A^{-1} u}, $$

so we just need to make appropriate substitutions:

$$ \begin{align} A &= X^\top X \\ u &= -x_k \\ v^\top &= x_k^\top. \end{align} $$

And of course we can repeat this for $r > 1$ indices and then we are splitting $A=X^\top X$ into the sum of two non-empty sets of outer products, $k\in \mathcal{S}$ and its complement. This leads us to the Woodbury identity, because now we have a rank-$r$ update to $A$. (Naturally, we can't leave out too many rows because then we have non-invertible matrix problems, and the procedure will blow up if the "denominator" is too close to 0, signaling that removing these rows is causing the matrix to become ill-conditioned.)

So the Woodbury identity will use

$$\begin{aligned} C &= I_{r\times r}\\ U &= -X_{k\in\mathcal{S}}^\top \\ V &= X_{k\in\mathcal{S}}. \end{aligned}$$

One caveat here is that we haven't characterized the loss of precision incurred by using floating-point arithmetic. Before implementing this in code, I would recommend studying the numerical conditioning of this procedure.

A colleague observes that eventually, for $r=|\mathcal{S}|$ too large, this becomes more expensive than the original problem. A better alternative is to form a QR factorization. This procedure is faster and more accurate and has its own update capabilities. I believe this is outlined in Golub & van Loan but I don't have my copy handy.

$\endgroup$
  • $\begingroup$ thank you thousands! I really got curious about what you said about QR factorization: "... is faster .. and has its own update capabilities." Would be something like using the Woodbury identiy on $R$ here $ \mathbf{\hat{\beta}} = R^{-1} Q^T z $ ? Golub & van Loan is this one here amazon.com/Computations-Hopkins-Studies-Mathematical-Sciences/…? Thank you millions! $\endgroup$ – eusoubrasileiro 39 mins ago

Your Answer

eusoubrasileiro is a new contributor. Be nice, and check out our Code of Conduct.

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.