Ask Your Question

Square roots of positive definite matrices

asked 2023-09-01 02:00:32 +0200

yeetcode gravatar image

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 flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2023-09-01 03:40:05 +0200

updated 2023-09-01 18:49:34 +0200

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

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?]
edit flag offensive delete link 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.

yeetcode gravatar imageyeetcode ( 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.

John Palmieri gravatar imageJohn Palmieri ( 2023-09-02 17:52:59 +0200 )edit

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

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

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2023-09-01 02:00:32 +0200

Seen: 223 times

Last updated: Sep 01 '23