scipy/scipy — pull request #2532 by argriffing (2013-06-01 07:09:40+00:00)
Title: ENH: sqrtm improvement

This PR implements yet another matrix function improvement published by Nick Higham and his crew. This time it is for the matrix square root. http://eprints.ma.man.ac.uk/1951/

According to that pdf, the matrix square root function is allegedly useful not only for its own sake, but also as an ingredient in algorithms to compute the matrix logarithm, matrix pth roots, and arbitrary matrix powers.

The following table suggests that using the optional `blocksize` argument added by this PR can increase `sqrtm` speed, especially for large `float64` matrices. When the `blocksize` is set to 1 (its default), the `sqrtm` function uses its current implementation rather than using the new blocked algorithm implemented by this PR.

The PR also includes tests and the informal benchmark that generated the table below, except without the slow `n=1024` cases.

---

Matrix Square Root
==============================================================
shape | block | dtype | time
| size | | (seconds)
--------------------------------------------------------------
(16, 16) | 1 | float64 | 0.00
(16, 16) | 32 | float64 | 0.01
(16, 16) | 1 | complex128 | 0.01
(16, 16) | 32 | complex128 | 0.00
(64, 64) | 1 | float64 | 0.09
(64, 64) | 32 | float64 | 0.05
(64, 64) | 1 | complex128 | 0.09
(64, 64) | 32 | complex128 | 0.04
(256, 256) | 1 | float64 | 2.73
(256, 256) | 32 | float64 | 0.47
(256, 256) | 1 | complex128 | 2.96
(256, 256) | 32 | complex128 | 0.71
(1024, 1024) | 1 | float64 | 152.54
(1024, 1024) | 32 | float64 | 8.47
(1024, 1024) | 1 | complex128 | 158.15
(1024, 1024) | 32 | complex128 | 16.01