Revision history [back]

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/develop/src/sage/matrix/matrix2.pyx#L15004):

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

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/develop/src/sage/matrix/matrix2.pyx#L15004):

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

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?]