I'm working on some code to find the inverse of a covariance matrix of arbitrary dimension. Obviously, the matrix will always be symmetric, and I think it's that fact that I need to take advantage of to make this easier, but am failing. Any advice other than brute force elementary operations and checks?


  • $\begingroup$ If this is for homework or self-study, you should add the self-study tag to your post. $\endgroup$ – Patrick Coulombe Feb 17 '14 at 4:03
  • 1
    $\begingroup$ (1) Why do you need the inverse? (2) There are some very common algorithms for matrix inversion. If the matrix is sure to be symmetric positive definite, you could use Cholesky decomposition (it's relatively easy to invert the triangular factor), but there are more stable approaches that are suitable even if it's only positive semi-definite, or nearly so. If you have the original variables from which the covariance is computed, the QR decomposition is generally a better place to start. $\endgroup$ – Glen_b -Reinstate Monica Feb 17 '14 at 4:17
  • 1
    $\begingroup$ Also see here $\endgroup$ – Glen_b -Reinstate Monica Feb 17 '14 at 4:22
  • $\begingroup$ Cholesky decomposition is a way to use the fact that covariance matrix is nonnegative definite and symmetric. Complexity for Cholesky decomposition seems to be smaller than that of other ways to find the inverse matrix. $\endgroup$ – Alexey Zaytsev Feb 17 '14 at 10:30

The symmetric matrix $A$ taken to the power of $\alpha$ can be computed by using spectral decomposition, i. e. $\text A^{\alpha} = \Gamma \Lambda^{\alpha} \Gamma'$. With $\Gamma$ being the matrix with the normalized eigenvectors of $A$ and $\Lambda$ being a diagonal matrix with the eigenvalues of $A$. Hence, the inverse of $A$ is computed by setting $\alpha=-1$. Now, one gets the inverse of the diagonal matrix $\Lambda$ by simply taking the inverse of every element of the diagonal matrix, i. e. $\Lambda^{-1} = diag(1/\lambda_1,\dots,1/\lambda_p) $.

An example for an R-Code looks like this:

A <- cov(X)    #Covariance of X
E <- eigen(A)  #Eigenvectors and -values of A
Lambda <- diag(E$values^-1)    #Diagonal matrix with inverse of eigenvalues
Gamma <- E$vectors  #Eigenvectors
A_Inverse <- Gamma%*%Lambda%*%t(Gamma)  #compute the inverse
round(A_Inverse%*%A,5)  #check

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.