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

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?

• Does X'X happen to be sparse? Jun 27 '17 at 14:48
• No. Is dense and medium size (1000 - 4000 rows) Jun 27 '17 at 15:20

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.
• One more comment before we leave the subject. I believe there is a small amount of additional structure that you could exploit during this process.That triangular backsolve, $R^{-T}e_i$, only operates on the trailing subtriangle of $R$ and the tail of $e_i$. For example, when processing the "next-to-last" $e_i$, the problem is effectively only 2x2. You can exploit this via judiciously "slicing" the arguments you pass into dtrsv(). The dot product, dnrm2(), can be sliced the same way. This should yield 2x reduction in flops. [At least, on those two phases. The cost of dgeqrf() is unchanged] Jun 30 '17 at 12:52