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

Question

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}$?

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

Answer 2

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.

Answer 3

The standard answer to your first question is Cholesky decomposition. To quote the Wikipedia article:

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

I'm not really clear what you're after in your second question. Doesn't a solution to the first question also provide a solution to the second? Surely it's not difficult computationally to remove a row and column from a matrix, and surely that's all the 'pre-proccessing' required??