7
$\begingroup$

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!

Improve this question
$\endgroup$
7
  • 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
    Commented Jan 3, 2012 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
    Commented Jan 3, 2012 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 ♦
    Commented Jan 3, 2012 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
    Commented Jan 3, 2012 at 15:45
  • 3
    $\begingroup$ Only slightly tongue-in-cheek: By using chol2inv or qr.solve. $\endgroup$
    – cardinal
    Commented Jan 3, 2012 at 15:52

1 Answer 1

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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