18

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.

tchakravarty
  • 345
  • 3
  • 8
  • 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}$. – tchakravarty 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 '14 at 07:25
  • 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. – roni May 12 '15 at 03:50

3 Answers3

12

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.

Brian Borchers
  • 18,719
  • 1
  • 39
  • 70
  • 7
    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
6

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
suvrit
  • 161
  • 1
  • 4
  • 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. – user67724 Sep 24 '23 at 16:28
0

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.

power
  • 521
  • 1
  • 4
  • 8