I incidently(*) discovered that using MatrixPower to find the square root matrix of a relatively small matrix, returns wrong results, even for symetric positive definite matrices.
Trying to understand how MatrixPower proceeds using showstat is not that easy, but something I'm almost sure is that it doesn't ise an Eigenvectors decomposition.
Here is an example where I'm seeking for a matrix M such that C = M^2 and C is symetric definite positive of dimension 40.
(I didn't check if problems already occurs for smaller sizes)
(*) I was very bored by the prohibitive execution time of MatrixPower and looked for an alternative.
The matrix M that MatrixPower generates is far from verifying C = M^2 (using the option 'shape'='symmetric' doesn't improve anything).
The naive alternative based on an Eigenvectors decomposition works perfectly well while being 200 times faster.