Fast way to compute the diagonal elements of the inverse of covariance matrix

Question

We have a variance-covariance matrix denoted with $X^TX$, where $X$ is the design matrix. In linear regression we can estimate beta coefficients with normal equations like $\hat{\beta} = (X^TX)^{-1}X^Ty$ and we can also compute the variance of betas with $Var(\hat{\beta}) = \hat{\sigma}^2(X^TX)^{-1}$, where $\hat{\sigma}^2$ is estimated from sample.

I implemented linear regression in my library and I used QR factorization for solving betas and also I used QR factorization for computing $(X^TX)^{-1}$, and then took only square roots of the diagonal elements of it. The last one I computed by solving $(X^TX)A = I$, where $A$ is what I am searching for.

Is there a faster way to compute only the diagonal elements of the inverse of $X^TX$? I know about Cholesky, I did not consider it, since it is also $O(n^3)$ and is less stable numerically. Is there a shortcut only for those diagonal elements?

Answer 1

You can write the (i,i)'th diagonal entry of $(X^TX)^{-1}$ as the product $d=e_i^T\cdot(X^TX)^{-1}\cdot e_i$, where $e_i$ is the i'th euclidean vector (all zeroes, except for a single one at the i'th position).

If you already possess the QR factorization $X=QR$, then $(X^TX)^{-1} = (R^TR)^{-1}$ and $d=e_i^T\cdot(R^TR)^{-1}\cdot e_i = (R^{-T}e_i)^T(R^{-T}e_i)$. If you introduce the vector $y = R^{-T}e_i$, then $d=y^Ty$.

This procedure maps readily onto LAPACK/BLAS: (a) Compute $X=QR$ using dgeqrf(), inplace. Afterward, the desired $R$ will be in triu($X$) (b) For each diagonal $d$, form $e_i$ then backsolve it by $R^{-T}$ from the left using dtrsv(), yielding $y$. Then $d=y^Ty$. This is dnrm2().

This might be exactly what you are doing, but your OP suggested you were perhaps forming the inverse $(X^TX)^{-1}$ explicitly, which I don't think it necessary. The computational cost of this posts approach is asymptotically the same $O(n^3)$, though maybe it saves you a temporary.