Pseudospectral Shattering, the Sign Function, and Diagonalization in Nearly Matrix Multiplication Time

Abstract

AbstractWe exhibit a randomized algorithm which, given a square matrix$$A\in \mathbb {C}^{n\times n}$$A∈Cn×nwith$$\Vert A\Vert \le 1$$‖A‖≤1and$$\delta >0$$δ>0, computes with high probability an invertibleVand diagonalDsuch that$$ \Vert A-VDV^{-1}\Vert \le \delta $$‖A-VDV-1‖≤δusing$$O(T_\mathsf {MM}(n)\log ^2(n/\delta ))$$O(TMM(n)log2(n/δ))arithmetic operations, in finite arithmetic with$$O(\log ^4(n/\delta )\log n)$$O(log4(n/δ)logn)bits of precision. The computed similarityVadditionally satisfies$$\Vert V\Vert \Vert V^{-1}\Vert \le O(n^{2.5}/\delta )$$‖V‖‖V-1‖≤O(n2.5/δ). Here$$T_\mathsf {MM}(n)$$TMM(n)is the number of arithmetic operations required to multiply two$$n\times n$$n×ncomplex matrices numerically stably, known to satisfy$$T_\mathsf {MM}(n)=O(n^{\omega +\eta })$$TMM(n)=O(nω+η)for every$$\eta >0$$η>0where$$\omega $$ωis the exponent of matrix multiplication (Demmel et al. in Numer Math 108(1):59–91, 2007). The algorithm is a variant of the spectral bisection algorithm in numerical linear algebra (Beavers Jr. and Denman in Numer Math 21(1-2):143–169, 1974) with a crucial Gaussian perturbation preprocessing step. Our result significantly improves the previously best-known provable running times of$$O(n^{10}/\delta ^2)$$O(n10/δ2)arithmetic operations for diagonalization of general matrices (Armentano et al. in J Eur Math Soc 20(6):1375–1437, 2018) and (with regard to the dependence onn)$$O(n^3)$$O(n3)arithmetic operations for Hermitian matrices (Dekker and Traub in Linear Algebra Appl 4:137–154, 1971). It is the first algorithm to achieve nearly matrix multiplication time for diagonalization in any model of computation (real arithmetic, rational arithmetic, or finite arithmetic), thereby matching the complexity of other dense linear algebra operations such as inversion andQRfactorization up to polylogarithmic factors. The proof rests on two new ingredients. (1) We show that adding a small complex Gaussian perturbation toanymatrix splits its pseudospectrum intonsmall well-separated components. In particular, this implies that the eigenvalues of the perturbed matrix have a large minimum gap, a property of independent interest in random matrix theory. (2) We give a rigorous analysis of Roberts’ Newton iteration method (Roberts in Int J Control 32(4):677–687, 1980) for computing the sign function of a matrix in finite arithmetic, itself an open problem in numerical analysis since at least 1986.