Compute square root of positive semidefinite matrix

[stumpc5]

This is more a question, I do not have checked how to actually do it.

An (Hermitian) positive semidefinite matrix M has a unique positive semidefinite square root T, this is, M = T^2.

I would like to have this square root implemented as I need it for some further computations. This is already possible numerically as in

sage: W = ReflectionGroup(['B',5])
sage: M = W.invariant_form()
sage: M
[ 1 -1  0  0  0]
[-1  2 -1  0  0]
[ 0 -1  2 -1  0]
[ 0  0 -1  2 -1]
[ 0  0  0 -1  2]

sage: D, S = M.eigenmatrix_right()
sage: S

[                   1                    1                    1                    1                    1]
[ 0.9189859472289948?  0.3097214678905701? -0.7153703234534297?  -1.830830026003773?  -2.682507065662362?]
[ 0.7635211184333675? -0.5943511444371404?  -1.203615623775565?  0.5211085581132027?   3.513337091666136?]
[ 0.5462003494572026?  -1.088155921225222?  0.3727855977717922?   1.397877389115792?  -3.228707415119565?]
[ 0.2846296765465703? -0.8308300260037728?   1.309721467890570?  -1.682507065662363?   1.918985947228995?]

sage: T = S*diagonal_matrix([sqrt(x) for x in D.diagonal()])*S^-1
sage: T^2 == M
True
sage: T.is_positive_definite()
True
sage: T.is_hermitian()
True

I wonder whether this guarantees the matrix square root to be positive definite whenever one starts with a positive definite matrix.

In this and a few other examples it actually is, but the documentation does not say anything about it.

CC: @dimpase @nthiery @sagetrac-jmichel @tscrim

Component: linear algebra

Keywords: matrix root, days93.51

Issue created by migration from https://trac.sagemath.org/ticket/23424

[dimpase]

comment:2

The example in the ticket description is actually an exact computation, just displayed in a slightly misleading way. Indeed, one computes in QQbar here:

sage: d1=D[1][1]
sage: d1
0.5857864376269049?
sage: d1.minpoly()
x^2 - 4*x + 2
sage: d1.sqrt().minpoly()
x^4 - 4*x^2 + 2
sage: type(d1)
<class 'sage.rings.qqbar.AlgebraicNumber'>

Thus, what is the question? This way it seems like an implementation is already there.

[stumpc5]

comment:3

Oh, I indeed missed this. Updated description.

[stumpc5]

Description changed:

--- 
+++ 
@@ -5,26 +5,32 @@
 I would like to have this square root implemented as I need it for some further computations. This is already possible numerically as in

-sage: W = RootSystem(['A',3]) -sage: C = W.cartan_matrix() -sage: D, S = C.eigenmatrix_right() -sage: D +sage: W = ReflectionGroup(['B',5]) +sage: M = W.invariant_form() +sage: M +[ 1 -1 0 0 0] +[-1 2 -1 0 0] +[ 0 -1 2 -1 0] +[ 0 0 -1 2 -1] +[ 0 0 0 -1 2]

-[ 2 0 0] -[ 0 0.5857864376269049? 0] -[ 0 0 3.414213562373095?] +sage: D, S = M.eigenmatrix_right() sage: S

-[ 1 1 1] -[ 0 1.414213562373095? -1.414213562373095?] -[ -1 1 1] +[ 1 1 1 1 1] +[ 0.9189859472289948? 0.3097214678905701? -0.7153703234534297? -1.830830026003773? -2.682507065662362?] +[ 0.7635211184333675? -0.5943511444371404? -1.203615623775565? 0.5211085581132027? 3.513337091666136?] +[ 0.5462003494572026? -1.088155921225222? 0.3727855977717922? 1.397877389115792? -3.228707415119565?] +[ 0.2846296765465703? -0.8308300260037728? 1.309721467890570? -1.682507065662363? 1.918985947228995?] + sage: T = Sdiagonal_matrix([sqrt(x) for x in D.diagonal()])S^-1 -sage: T

-[ 1.360388263624736? -0.3826834323650898? -0.05382529874835926?] -[ -0.3826834323650898? 1.306562964876377? -0.3826834323650898?] -[-0.05382529874835926? -0.3826834323650898? 1.360388263624736?] -sage: T^2 == C +sage: T^2 == M +True +sage: T.is_positive_definite() +True +sage: T.is_hermitian() True

-It would be great to (1) have this exact (in some field extension in which the square roots live) and (2) to have the positive semidefinite decomposition.
+I wonder whether this guarantees the matrix square root to be positive definite whenever one starts with a positive definite matrix.
+
+In this and a few other examples it actually is, but the documentation does not say anything about it.

[dimpase]

comment:4

IMHO it is more of a mathematical question whether you get a p.s.d. T.

In the basis of the columns of S, one has the quadratic form, given by M in the standard basis, diagonal, given by D, with nonnegative entries (i.e. M=SDS^-1). The square root T is a power series in M, and thus given by T=SD^1/2S^-1, a p.s.d. matrix with eigenvalues given by D^1/2, thus all positive whenever D>0.

What can still go wrong is that T is not equal to T^*.

However we know that D=D and so M=M=SDS^-1=S^-1DS, implying SSD=DSS. Thus any power of D commutes with SS, thus SSD^1/2=D^1/2SS implying T=T.

Anything I still miss?

[stumpc5]

Changed keywords from matrix root to matrix root, days93.51

[stumpc5]

comment:5

duplicate of #25464