Take the 2-minute tour ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

A common problem in statistics is computing the square root inverse of a symmetric positive definite matrix. What would be the most efficient way of computing this?

I came across some literature (which I haven't read yet), and some incidental R code here, which I will reproduce here for convenience

# function to compute the inverse square root of a matrix
fnMatSqrtInverse = function(mA) {
  ei = eigen(mA)
  d = ei$values
      d = (d+abs(d))/2
      d2 = 1/sqrt(d)
      d2[d == 0] = 0
      return(ei$vectors %*% diag(d2) %*% t(ei$vectors))

I am not entirely sure I understand the line d = (d+abs(d))/2. Is there a more efficient way of computing the matrix square root inverse? The R eigen function calls LAPACK.

share|improve this question
Provided $d$ is real, $(d+|d|)/2$ is equal to $\max(d,0)$. This effectively gets rid of any negative eigenvalues that the matrix may have. Do you need all the entries of the matrix $A^{-1/2}$, or is it sufficient to be able to multiply $A^{-1/2}$ by an arbitrary vector $x$? –  Daniel Shapero Dec 21 '13 at 22:06
@DanielShapero Thanks for your comment. So if I have a PSD matrix, I don't need that line? My application requires computing quadratic forms such as $A^{-1/2}BA^{-1/2}$. –  fg nu Dec 21 '13 at 22:18
I am not familiar with R, but given line 7 I assume that it has logical indexing like Matlab. If so, I suggest you to rewrite line 5 as d[d<0] = 0, which is more expressive. –  Federico Poloni Apr 14 at 7:25

3 Answers 3

up vote 7 down vote accepted

The code that you've posted uses the eigenvalue decomposition of the symmetric matrix to compute $A^{-1/2}$.

The statement


effectively takes any negative entry in d and sets it to 0, while leaving non-negative entries alone. That is, any negative eigenvalue of $A$ is treated as though it was 0. In theory, the eigenvalues of A should all be non-negative, but in practice it's common to see small negative eigenvalues when you compute the eigenvalues of a supposedly positive definite covariance matrix that is nearly singular.

If you really need the inverse of the symmetric matrix square root of $A$, and $A$ is reasonably small (no bigger than say 1,000 by 1,000), then this is about as good as any method you might use.

In many cases you can instead use a Cholesky factor of the inverse of the covariance matrix (or practically the same, the Cholesky factor of the covariance matrix itself.) Computing the Cholesky factor is typically an order of magnitude faster than computing the eigenvalue decomposition for dense matrices and vastly more efficient (both in computational time and required storage) for large and sparse matrices. Thus using the Cholesky factorization becomes very desirable when $A$ is large and sparse.

share|improve this answer
Brian's answer gives good advice: use the Cholesky factor instead (if you can). There is another optimization you can do on top of that: don't compute your PSD matrix $A$. Often you get $A$ from a computation like $A=B^TB$, with $B$ rectangular. In this case, it is sufficient to compute a QR decomposition of $B$ to obtain the Cholesky factor $R$ of $A$, with much better accuracy. –  Federico Poloni Dec 28 '13 at 16:51

if the eigenvalues of a matrix are negative then how could we compute the square root of that matrix? however if eigenvalues are positive it could be computed easily.

share|improve this answer
The question explicitly refers to symmetric positive definite matrices, for which the eigenvalues are always positive, but it asks for the square root of the inverse of the matrix. Can you expand on how this can "be computed easily"? –  Christian Clason Jul 2 at 7:28

Optimise your code:

Option 1 - Optimise your R code:
a. You can apply() a function to d that will both max(d,0) and d2[d==0]=0 in one loop.
b. Try operating on ei$values directly.

Option 2 - Use C++:
Rewrite the whole function in C++ with RcppArmadillo. You will still be able to call it from R.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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