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

We exhibit a randomized algorithm which, given a square matrix A∈ C^n×n with ‖ A‖ ≤ 1 and δ> 0 , computes with high probability an invertible V and diagonal D such that ‖ A- VDV^{- 1}‖ ≤ δ using O(T_MM(n) log ²(n/ δ)) arithmetic operations, in finite arithmetic with O(log ⁴(n/ δ) log n) bits of precision. The computed similarity V additionally satisfies ‖ V‖ ‖ V^{- 1}‖ ≤ O(n^2.5/ δ). Here T_MM(n) is the number of arithmetic operations required to multiply two n× n complex matrices numerically stably, known to satisfy T_MM(n) = O(n^ω+η) for every η> 0 where ω 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¹⁰/ δ²) 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 on n) O(n³) 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 and QR factorization up to polylogarithmic factors. The proof rests on two new ingredients. (1) We show that adding a small complex Gaussian perturbation to any matrix splits its pseudospectrum into n small 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.