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.

  • $\begingroup$ 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$? $\endgroup$ Commented Dec 21, 2013 at 22:06
  • $\begingroup$ @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}$. $\endgroup$ Commented Dec 21, 2013 at 22:18
  • $\begingroup$ 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. $\endgroup$ Commented Apr 14, 2014 at 7:25
  • $\begingroup$ Is this code correct ? I ran it on a simple example in matlab and found the answer to be wrong. My matrix is positive definite but definitely not symmetric. Please see my answer below : I have transferred the code to matlab. $\endgroup$
    – roni
    Commented May 12, 2015 at 3:50

3 Answers 3


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.

  • 7
    $\begingroup$ 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. $\endgroup$ Commented Dec 28, 2013 at 16:51

In my experience, the polar-Newton method of Higham works much faster (see Chapter 6 of Functions of Matrices by N. Higham). In this short note of mine there are plots that compare this method to first-order methods. Also, citations to several other matrix-square-root approaches are presented, though mostly the polar Newton iteration seems to work the best (and avoids doing eigenvector computations).

% compute the matrix square root; modify to compute inverse root.
function X = PolarIter(M,maxit,scal)
  fprintf('Running Polar Newton Iteration\n');
  skip = floor(maxit/10);
  I = eye(size(M));
  if scal
    tm = trace(M);
    M  = M / tm;
    tm = 1;
  nm = norm(M,'fro');

  % to compute inv(sqrt(M)) make change here

  % computes the polar decomposition of R
  U=R; k=0;
  while (k < maxit)
    % err(k) = norm((R'*U)^2-M,'fro')/nm;
    %if (mod(k,skip)==0)
    %  fprintf('%d: %E\n', k, out.err(k));


   if (err(k) < 1e-12), break; end
  X = 0.5*(X+X');
  • $\begingroup$ What is the advantage of this iterative scheme over Cholesky decomposition? Especially, since this algorithm performs a Cholesky decomposition of a perturbation of M, it cannot be faster than Cholesky. $\endgroup$
    – user67724
    Commented Sep 24, 2023 at 16: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.


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.