I have a design matrix of p regressors, n observations, and I am trying to compute the sample variance-covariance matrix of the parameters. I am trying to directly calculate it using svd.

I am using R, when I take svd of the design matrix, I get three components: a matrix $U$ which is $n \times p$, a matrix $D$ which is $1\times 3$ (presumably eigenvalues), and a matrix $V$ which is $3\times 3$. I diagonalized $D$, making it a $3\times 3$ matrix with 0's in the off-diagonals.

Supposedly, the formula for covariance is: $V D^2 V'$, however, the matrix does not match, nor is it even close to R's built in function, vcov. Does anyone have any advice/references? I admit that I am a bit unskilled in this area.


First, recall that under assumptions of multivariate normality of the linear-regression model, we have that $$ \hat{\beta} \sim \mathcal{N}( \beta, \sigma^2 (X^T X)^{-1} ) . $$

Now, if $X = U D V^T$ where the right-hand side is the SVD of X, then we get that $X^T X = V D U^T U D V = V D^2 V^T$. Hence, $$ (X^T X)^{-1} = V D^{-2} V^T . $$

We're still missing the estimate of the variance, which is $$ \hat{\sigma}^2 = \frac{1}{n - p} (y^T y - \hat{\beta}^T X^T y) . $$

Though I haven't checked, hopefully vcov returns $\hat{\sigma}^2 V D^{-2} V^T$.

Note: You wrote $V D^2 V^T$, which is $X^T X$, but we need the inverse for the variance-covariance matrix. Also note that in $R$, to do this computation you need to do

vcov.matrix <- var.est * (v %*% d^(-2) %*% t(v))

observing that for matrix multiplication we use %*% instead of just *. var.est above is the estimate of the variance of the noise.

(Also, I've made the assumptions that $X$ is full-rank and $n \geq p$ throughout. If this is not the case, you'll have to make minor modifications to the above.)


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.