scipy/scipy — pull request #3556 by argriffing (2014-04-16 21:29:58+00:00)
Title: ENH: matrix square root for positive semi-definite matrices

Inspired by https://github.com/scipy/scipy/issues/3549 this PR adds `funm_psd` to evaluate general matrix functions of positive semi-definite matrices, and in particular adds `sqrtm_psd` specialized to the matrix square root. It includes benchmarks and tests. For real positive semi-definite matrices near size N=1000, `sqrtm_psd` is about twice as fast as the less specific `sqrtm` and correctly returns real (as opposed to complex) matrices for real positive semi-definite input matrices.

Interface advice would be welcome, for example maybe it would be better to add more kwargs to existing functions rather than to add new separate functions specialized to positive semi-definite matrices. Another alternative would be to just provide a recipe for these functions in the scipy documentation rather than adding new functionality, because the implementations are thin wrappers around `eigh`.

```
Positive Semi-Definite Matrix Square Root
================================================================
shape | method | dtype | time
| | | (seconds)
----------------------------------------------------------------
(64, 64) | sqrtm | float64 | 0.07
(64, 64) | sqrtm_psd | float64 | 0.02
(64, 64) | sqrtm | complex128 | 0.20
(64, 64) | sqrtm_psd | complex128 | 0.01
(256, 256) | sqrtm | float64 | 1.14
(256, 256) | sqrtm_psd | float64 | 0.55
(256, 256) | sqrtm | complex128 | 1.57
(256, 256) | sqrtm_psd | complex128 | 0.62
(1024, 1024) | sqrtm | float64 | 41.72
(1024, 1024) | sqrtm_psd | float64 | 17.55
(1024, 1024) | sqrtm | complex128 | 71.01
(1024, 1024) | sqrtm_psd | complex128 | 23.55
```