Using Singular Value Decomposition to Compute Variance Covariance Matrix from linear regression model 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.
 A: 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.)
