Efficient computation of the matrix square root inverse

Question

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.

Answer 1

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

The statement

d=(d+abs(d))/2

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.

Answer 2

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));
  n=size(M,1);
  if scal
    tm = trace(M);
    M  = M / tm;
  else
    tm = 1;
  end
  nm = norm(M,'fro');

  % to compute inv(sqrt(M)) make change here
  R=chol(M+5*eps*I);

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

    iU=U\I;
    mu=sqrt(sqrt(norm(iU,1)/norm(U,1)*norm(iU,inf)/norm(U,inf)));
    U=0.5*(mu*U+iU'/mu);

   if (err(k) < 1e-12), break; end
  end
  X=sqrt(tm)*R'*U;
  X = 0.5*(X+X');
end

Answer 3

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.