Computing Krylov iterates in the time of matrix multiplication

Abstract

Krylov methods rely on iterated matrix-vector products $A^k{u}_j$ for an $n\times n$ matrix $A$ and vectors $u_1,\dots ,u_m$. The space spanned by all iterates $A^k{u}_j$ admits a particular basis -- the \emph{maximal Krylov basis} -- which consists of iterates of the first vector $u_1,A{u}_1,A^2{u}_1,\dots$, until reaching linear dependency, then iterating similarly the subsequent vectors until a basis is obtained. Finding minimal polynomials and Frobenius normal forms is closely related to computing maximal Krylov bases. The fastest way to produce these bases was, until this paper, Keller-Gehrig's 1985 algorithm whose complexity bound $O(n^{\ddot{I}}‰\log (n))$ comes from repeated squarings of $A$ and logarithmically many Gaussian eliminations. Here $\ddot{I}‰>2$ is a feasible exponent for matrix multiplication over the base field. We present an algorithm computing the maximal Krylov basis in $O(n^{\ddot{I}}‰\log \log (n))$ field operations when $m\in O(n)$, and even $O(n^{\ddot{I}}‰)$ as soon as $m\in O(n/\log (n{)}^c)$ for some fixed real $c>0$. As a consequence, we show that the Frobenius normal form together with a transformation matrix can be computed deterministically in $O(n^{\ddot{I}}‰(\log \log (n){)}^2)$, and therefore matrix exponentiation~$A^k$ can be performed in the latter complexity if $\log (k)\in O(n^{\ddot{I}‰-1-\epsilon })$ for some fixed $\epsilon >0$. A key idea for these improvements is to rely on fast algorithms for $m\times m$ polynomial matrices of average degree $n/m$, involving high-order lifting and minimal kernel bases.

10 pages, 5 algorithms