# Computing $(X^TX)^{-1}X^Ty$ in OLS

Let $A\in\mathbb{R}^{n \times n}$ be a dense symmetric positive-definite matrix (the $X^TX$ from here) and $b$ a vector in $\mathbb{R}^n$.

I need to compute $A^{-1}b$.

Two questions:

1. Could you recommend an efficient and numerically stable algorithm for computing $A^{-1}b$ for $n \approx 1000$?

2. Let $\tilde{A_i}$ denote the matrix obtained from $A$ by removing its $i$-th row and $i$-th column. Is there an algorithm that, having been allowed to pre-process $A$ in some way, would enable me to quickly compute $\tilde{A_i}^{-1}\tilde b$ for any $i \in \{1,2,\ldots,n\}$ and any $\tilde b \in \mathbb{R}^{n-1}$?

-

Regarding your second question, here is a way to do this for $n=i$, for simplicity of notation. Let $\alpha = A_{nn}, \, a = (A_{1n},\dots, A_{n-1,n})^T$ and therefore $$A = \begin{pmatrix} \tilde A & a \\ a^T & \alpha \end{pmatrix}$$ Also let $A^{-1}$ be partitioned in the same way, $$A^{-1} = \begin{pmatrix} \tilde C & c \\ c^T & \gamma \end{pmatrix}$$ I understand you want $\tilde A^{-1}$ and have already computed $A^{-1}$. Set $d = \tilde A^{-1}a$ and $\delta = c^Ta$. These quantities already appear in $A^{-1}$, since $$A^{-1} = \begin{pmatrix} \tilde A^{-1} + \frac{1}{\alpha - \delta}dd^T & -\frac{1}{\alpha - \delta} d \\ -\frac{1}{\alpha - \delta} d^T & \frac{1}{\alpha - \delta} \end{pmatrix}$$ as a direct calculation shows. Therefore, $$\tilde A^{-1} = \tilde C - \frac{1}{\alpha - \delta} dd^T = \tilde C - \frac{1}{\gamma}cc^T$$ which is an $O(n^2)$ update of $A^{-1}$.

-
However, as a general rule of thumb, downdates are more prone to numerical instability than updates. –  Ｊ. Ｍ. May 12 '11 at 13:33

To add to @onestop's answer, another efficient way is to use QR decomposition. The added benefit is that QR decomposition can be applied directly to $X$, and not to $X^TX$.

I think the QR decomposition can be made to work for your second question, it is definitely straightforward for column removal.

However the question is why do you need it? There are readily available libraries which perform these operations very efficiently. Do you want to reimplement them? As far as I understand there is a lot of non-trivial fine-tuning of code even if the algorithm is basically the best for the job, so chances are pretty high that you might not get full benefits by implementing algorithm which is supposedly theoretically superior.

Here is the link to the chapter in a book, which discusses modifications needed for solving the second question in case of QR decomposition. Although it states that the methods apply for least squares problems, it can be applied for system of linear equations.

I am pretty sure that this should be a standard problem, so maybe someone will give a more suitable reference.

-
@mpiktas: Thanks for this. Just to clarify, I don't have $X$. I only have $X^TX$ -- see stats.stackexchange.com/questions/9341/… –  NPE Apr 14 '11 at 10:59
@aix, QR decomposition can still be applied, although Cholesky might be more appropriate. –  mpiktas Apr 14 '11 at 11:07
@mpiktas: To address your second point, I am not saying I am going to code up my own SVD/QR/Cholesky/whatever (unless I really have to, that is). –  NPE Apr 14 '11 at 11:08
Apart from Björck, the other neat reference for solving least-squares problems is the (classic) Lawson/Hanson book. –  Ｊ. Ｍ. May 12 '11 at 11:49
@J. M., thanks! This was first book I read on this topic, but I forgot the authors and the exact title. Very good book. It also has Fortran code examples, which was very useful for me too. –  mpiktas May 12 '11 at 12:00
If $A$ is symmetric and positive definite, then we can solve $Ax = b$ [for $x$] by first computing the Cholesky decomposition $A = LL^\mathrm{T}$, then solving $Ly = b$ for $y$, and finally solving $L^\mathrm{T} x = y$ for $x$.