I did not find an implementation for square roots of positive definite matrices, so I decided to write one of my own. I diagonalised it as $A = PDP^{-1}$ by using D, P = A.eigenmatrix_right(), took an entry-wise square root of $D$, and made my function return that. Here is the exact implementation:
``` Code:
def matrix_root(A):
D,P = A.eigenmatrix_right();
for i in range(A.nrows()):
D[i,i] = abs(sqrt(D[i,i]));
return (P*D*P.inverse())
```
However, for some inputs $A$, the code generates a runtime error. In particular, letting $C = \sqrt{A}$ and computing $\sqrt{C^{-1} B C^{-1}}$ for $$B = \begin{pmatrix} 21 & 33 & 42 \\ 33 & 83 & 65 \\ 42 & 65 & 89 \end{pmatrix}$$
makes the code run for far too long. I have lost the error log now, and will need me to run the code till it declares a runtime error to find it again. But if I understand things correctly, it was doing an RREF computation and it was trying to figure out if two algebraic numbers cancel each other out, and it somehow managed to run into an algebraic number whose minimal polynomial was just 1? I have no idea what's going on.
I am assuming it is some mix of decimal precision issues and the general difficulty in computations in QQbar. I am quite confused and am very unsure as to how to go about doing this whole thing.
