Opened 4 years ago

Closed 3 years ago

#23424 closed enhancement (duplicate)

Compute square root of positive semidefinite matrix

Reported by: stumpc5 Owned by:
Priority: major Milestone: sage-8.1
Component: linear algebra Keywords: matrix root, days93.51
Cc: dimpase, nthiery, jmichel, tscrim Merged in:
Authors: Reviewers:
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Status badges

Description (last modified by 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 = T2.

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
sage: T.is_positive_definite()
sage: T.is_hermitian()

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.

Change History (5)

comment:1 Changed 4 years ago by stumpc5

  • Cc jmichel added

comment:2 Changed 4 years ago by dimpase

  • Status changed from new to needs_info

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
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.

comment:3 Changed 4 years ago by stumpc5

  • Description modified (diff)

Oh, I indeed missed this. Updated description.

comment:4 Changed 4 years ago by dimpase

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=SD1/2S-1, a p.s.d. matrix with eigenvalues given by D1/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 S*SD=DS*S. Thus any power of D commutes with S*S, thus S*SD1/2=D1/2S*S implying T=T*.

Anything I still miss?

comment:5 Changed 3 years ago by stumpc5

  • Cc tscrim added
  • Keywords days93.51 added
  • Resolution set to duplicate
  • Status changed from needs_info to closed

duplicate of #25464

Note: See TracTickets for help on using tickets.