[SciPy-Dev] matrix_sqrt for singular symmetric square matrices

[Pauli Virtanen]

su, 2018-10-28 kello 08:56 +0100, Gael Varoquaux kirjoitti:
> On Sat, Oct 27, 2018 at 04:16:58PM -0400, josef.pktd at gmail.com wrote:
> >     A matrix square-root would be very useful.  As noted here, the
> > matrix
> >     square-root is typically non-unique:
> >     https://en.wikipedia.org/wiki/Square_root_of_a_matrix
> > Which basis is chosen is irrelevant in most use cases that I know.
> 
> In some case, the following considerations matter:
> 
> SDP matrices (symmetric definite positive matrices) form an algebraic
> structure linked to a group. The square root matrix on this group
> should
> be also SDP, hence it should be also symmetric. 

Right indeed if R=R' the two equations are the same. scipy.linalg.sqrtm
computes the principal matrix square root, which then should be SDP
too.

It can handle singular cases. I would guess the algorithm in the
presence of eigenvalues on the negative real line gives "continuous"
behavior, i.e., you get the principal square root of a matrix with
eigenvalues pushed to either side of the branch cut in some
(uncontrolled) way giving either of the +/- i sqrt(|z|) and the result
probably is close to Hermitian still in the same way as the eigenvalue
decomposition.

So for the question whether there's something to add to scipy here, I'm
not so sure --- computation of the principal sqrtm, Cholesky, LDLt, and
eigenvalue decomposition is there.

	Pauli


> The easiest way to build
> it is using the eigen-value decomposition of the original matrix:
> 
>     M = U' @ S @ U
> 
> using '@' to denote the matrix product, and S = np.diag(s) is the
> diagonal matrix of eigenvalues. The matrix square root is then given
> by:
> 
>    sqrt(M) = U' @ np.diag(np.sqrt(s)) @ U
> 
> Unlike with using a Cholesky decomposition to obtain a square-root
> matrix, the operation defined above is smooth.
> 
> Such considerations are typically important in signal processing, to
> manipulate covariance matrices for whitening and averaging.
> 
> Gaël
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at python.org
> https://mail.python.org/mailman/listinfo/scipy-dev