# 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!

• R is calculating standard error without inverting the matrix. R is doing QR decomposition. See this [link] . See page no 2 and 3 : cran.r-project.org/doc/contrib/Leisch-CreatingPackages.pdf – vinux Jan 3 '12 at 14:21
• Yes, please don't use Cholesky for OLS, especially if the data matrix is of any appreciable size. The accumulation of numerical roundoff error is often practically significant. – cardinal Jan 3 '12 at 15:10
• You might be able to capitalize on special circumstances to speed up the calculation. See, for instance, stats.stackexchange.com/q/6920 – whuber Jan 3 '12 at 15:29
• If $X = Q R$, then $X^T X = R^T Q^T Q R = R^T R$. Recall that the covariance matrix of linear regression is $\sigma^2 (X^T X)^{-1}$. But $R$ is an upper triangular matrix, so $(R^T R)^{-1}$ can be found by one backward and one forward substitution. – cardinal Jan 3 '12 at 15:45
• Only slightly tongue-in-cheek: By using chol2inv or qr.solve. – cardinal Jan 3 '12 at 15:52

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.