Computing standard errors for linear regression problems without calculating inverse

Question

Is there a speedier way to calculate standard errors for linear regression problems, than by inverting $X'X$? Here I assume we have regression:

$$y=X\beta+\varepsilon,$$

where $X$ is $n\times k$ matrix and $y$ is $n\times 1$ vector.

For finding least squares problem solution it is impractical to do anything with $X'X$, you can use QR or SVD decompositions on matrix $X$ directly. Or alternatively you can use gradient methods. But what about standard errors? We really only need the diagonal of $(X'X)^{-1}$ (and naturally LS solution to calculate the estimate of standard error of $\varepsilon$). Are there any specific methods for standard error calculation?

Answer 1

Let's suppose that you solved your least squares problem using the singular value decomposition (SVD) of $X$, given by

$$X = U\Sigma V',$$

where $U$ and $V$ are unitary, and $\Sigma$ is diagonal.

Then

$$X'X = V \Sigma^2 V'.$$

$(X'X)^{-1}$ exists iff $X$ is full rank (or has strictly positive singular values), in which case

$$(X'X)^{-1} = V \Sigma^{-2} V'.$$

(See an answer I gave to a related question on Math.SE.)

If you already have $\Sigma$ and $V$, calculating $(X'X)^{-1}$ requires inverting and squaring a diagonal matrix ($n$ operations for an $n \times n$ matrix), scaling the columns (or rows) of a matrix ($n^2$ operations), and a single matrix multiply (unfortunately $\mathcal{O}(n^{3})$). This method will be well-behaved numerically.

There are fast methods for obtaining the diagonal elements of the inverse of a sparse matrix (see work by Yousef Saad's group and work by Lin Lin, et al). However, in your case, $X'X$ is probably not sparse (even if $X$ is), and even if it were, it would likely be ill-conditioned enough that these fast methods would yield inaccurate results.