Blocked Schur Algorithms for Computing the Matrix Square Root

Abstract

Abstract. The Schur method for computing a matrix square root reduces the matrix to the Schur triangular form and then computes a square root of the triangular matrix. We show that by using either standard blocking or recursive blocking the computation of the square root of the triangular matrix can be made rich in matrix multiplication. Numerical experiments making appropriate use of level 3 BLAS show significant speedups over the point algorithm, both in the square root phase and in the algorithm as a whole. In parallel implementations, recursive blocking is found to provide better performance than standard blocking when the parallelism comes only from threaded BLAS, but the reverse is true when parallelism is explicitly expressed using OpenMP. The excellent numerical stability of the point algorithm is shown to be preserved by blocking. These results are extended to the real Schur method. Blocking is also shown to be effective for multiplying triangular matrices.

Matrix Functions: A Short Course

“…Recently, Deadman, Higham, and Ralha [18] have developed blocked versions of the recurrence (6.2) that give substantially better performance on modern computers.…”

An Algorithm for the Matrix Lambert $W$ Function

“…The matrix exponentials and logarithms in step 4 are evaluated using the algorithms of Al-Mohy and Higham [1], [2], which exploit triangularity. The matrix square roots in step 4b are evaluated using the Björck-Hammarling recurrence [8], with the efficient blocked implementation of [16].…”

Estimating the condition number of f(A)b

“…-the principal matrix logarithm, available in SciPy as scipy.linalg.logm, computed using inverse scaling and squaring [3]; -the matrix sine, scipy.linalg.sinm, computed using the Schur-Parlett algorithm [8]; -the matrix cosine, scipy.linalg.cosm, also computed using the SchurParlett algorithm; -the matrix square root, scipy.linalg.sqrtm, computed using a blocked version of the Björck-Hammarling algorithm [5], [12]; -the matrix cube root, scipy.linalg.fractional matrix power, computed using the Schur-Padé algorithm [21].…”