How to calculate standard errors in OLS without inverting the X'X matrix? This may seem a trivial question but I haven't found a satisfactory answer anywhere. I need to compute standard errors in a OLS regression $y = X\beta + u$ in R from scratch. How can I do this without inverting $X'X$ matrix? (The reason is that speed is important.) Suppose that I already have factorized the $X'X$ matric by Cholesky and already have $\hat \beta$ and $\hat \sigma$. Thanks!
 A: I had the same problem as I wanted to use the most efficient solvers available in my econometrics package. I developed an algorithm to solve for both the $\beta$ and the inverse of the linear predictor (normal matrix) for linear least squares (it also applies to WLS, Ridge regression, etc.)
So here is the pseudo-code,
Given $X$ and $y$,
1) use an efficient solver to estimate $\beta$.
2) Store $X^{\top}X + \Gamma \Gamma^{\top}$
3) Create a new $Ax=b$ problem with:


*

*x = vector of values of values that fully characterize $\left(X^{\top}X + \Gamma \Gamma^{\top}\right)^{-1}$. It is symmetric so there are $\left(\frac{k (k + 1)}{2}\right)$

*b = vec($\beta$, vec($I_{k}$))

*A = Linear predictor that establishes the following restrictions:


*

*$A X^{\top} y = \beta$

*$\left(X^{\top}X + \Gamma \Gamma^{\top}\right)^{-1}\left(X^{\top}X + \Gamma \Gamma^{\top}\right) = I$
Lastly, solve the new $Ax=b$ for $x$ using an efficient solver and recreate $A$ with it.
This way, you don't have to calculate the Cholesky decomposition or anything. A good solver could be lsmr which can give more precision and a more efficient implementation.
