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!

  • 1
    $\begingroup$ R is calculating standard error without inverting the matrix. R is doing QR decomposition. See this [link][1] . See page no 2 and 3 [1]: cran.r-project.org/doc/contrib/Leisch-CreatingPackages.pdf $\endgroup$ – vinux Jan 3 '12 at 14:21
  • 2
    $\begingroup$ 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. $\endgroup$ – cardinal Jan 3 '12 at 15:10
  • 1
    $\begingroup$ You might be able to capitalize on special circumstances to speed up the calculation. See, for instance, stats.stackexchange.com/q/6920 $\endgroup$ – whuber Jan 3 '12 at 15:29
  • 8
    $\begingroup$ 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. $\endgroup$ – cardinal Jan 3 '12 at 15:45
  • 3
    $\begingroup$ Only slightly tongue-in-cheek: By using chol2inv or qr.solve. $\endgroup$ – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.