# Square roots of positive definite matrices

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.

edit retag close merge delete

Sort by ยป oldest newest most voted

The method principal_square_root is present in Sage and seems to work (at a quick glance at the source code) the way your function does, but it runs pretty quickly for that matrix B:

sage: B = matrix([[21,33,42],[33,83,65],[42,65,89]])
sage: %time B.principal_square_root()
CPU times: user 35.7 ms, sys: 1.07 ms, total: 36.8 ms
Wall time: 36.2 ms
[2.477638707052645? 1.988175360680636? 3.302796568439680?]
[1.988175360680636? 8.125503327539394? 3.608788496064663?]
[3.302796568439680? 3.608788496064663? 8.066485010100398?]


Main part of the source code (see https://github.com/sagemath/sage/blob...):

d, L = self.eigenmatrix_left()
return L.inverse() * diagonal_matrix([sqrt(a) for a in d.diagonal()]) * L


Edit: In response to your comment, I'm not sure which "exact computation" you mean. Here is what I tried:

sage: A = matrix(QQbar, [[2,-1,0],[-1,2,1],[0,1,2]])
sage: B = matrix([[21,33,42],[33,83,65],[42,65,89]])
sage: C = A.principal_square_root()
sage: D = C.inverse() * B * C.inverse()
sage: %time D.principal_square_root()
CPU times: user 770 ms, sys: 2.14 ms, total: 772 ms
Wall time: 772 ms
[3.029530081?  3.22867050?  2.19389687?]
[3.228670494? 6.748880792? 0.881494325?]
[2.193896864? 0.881494325? 5.074002154?]

more

Thanks! While it works for $B$, it still fails for the matrix in my context. I think nested roots are too problematic somehow. Can you see if you can make this exact computation work? That would be plenty helpful!

Edit: I deeply apologise. The issue actually arises when you swap $A$ and $B$, I got it wrong. I tend to make such mistakes too often. I would prefer adding a screenshot of the code, but my karma seems to be too low for that.

( 2023-09-01 16:36:26 +0200 )edit

I now see a problem. I interrupted the computation since it was taking so long, so I don't know if it ends in an error. I agree with your guess: it's difficult to do computations in QQbar.

( 2023-09-02 17:52:59 +0200 )edit

In that case, do you believe there is any alternative method?

( 2023-09-02 22:47:07 +0200 )edit