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}$.
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))
