# Accurate way of getting the square root inverse of a positive definite symmetric matrix

What is the most accurate algorithm to get the square root inverse of a positive definite symmetric matrix? I am not looking as much for efficiency, though using quadruple precision computation is out of the question.

By square root inverse I mean, given my matrix $$J$$, its square root inverse $$J^{-1/2}$$ is such that

$$J^{-1/2} J J^{-1/2} = I$$.

Another way to put it is

$$J^{-1/2} J^{-1/2} = J^{-1}$$.

In particular, I'm interested in the positive definite $$J^{-1/2}$$.

• Do you need the square root? Often one can do with a (say Cholesky) factorization. – Jan Jun 2 '19 at 4:56
• I forgot to mention that my matrix is real, positive definite. So, I'd like to get the real, positive definite square root of the inverse of that matrix. – Javier Garcia Jun 2 '19 at 13:52
• Well you can do it by diagonalisation which in my experience is pretty stable for symmetric matrix, but I'm not numerical analyst. But as the previous comment says id you can do what you want with the Choleski factors that will be much quicker – Ian Bush Jun 2 '19 at 14:26
• @JavierGarcia OK if you need the symmetric positive square root the other factorizations are out of the game. Maybe you can rephrase your question accordingly. – Jan Jun 3 '19 at 9:19

As suggested in the comments, the eigendecomposition $$\mathbf J = \mathbf Q \mathbf D \mathbf Q^{T}$$ can be used to generate the matrix $$\mathbf J^{-1/2}$$, just take the eigenvectors from $$\mathbf J$$ (denoted $$\mathbf Q$$) and the inverse square root of the eigenvalues (denoted $$\mathbf D^{-1/2}$$) .. this is a "scalar"/"entrywise" operation because $$\mathbf D$$ is diagonal.

A short snippet of matlab code can demonstrate:

% Make random NxN SPD matrix J
N = 10;
A = rand(N);
J = A'*A;
% Compute eigendecomposition, J = Q*D*Q'
[Q,D] = eig(J);
% Form Z := J^{-1/2}
Z = Q*diag(diag(D).^(-1/2))*Q';
% Verify "square root" property
norm(Z*J*Z - eye(N))