Computing Krylov iterates in the time of matrix multiplication
Vincent Neiger, Clément Pernet, Gilles Villard
TL;DR
This work advances the efficient computation of Krylov structures and related canonical forms for a matrix $A$ by replacing costly exponentiation with polynomial-matrix methods. It introduces a deterministic framework that uses the fraction description $S(x)T(x)^{-1}=(I-xA)^{-1}U$ and high-order lifting to extract maximal Krylov bases, achieving $O(n^{\omega}\log\log n)$ complexity when $m=O(n)$ and $O(n^{\omega})$ when $m=O(n/\log^c n)$. The authors connect maximal Krylov indices to diagonal degrees of a Hermite form and leverage kernel bases to compute these indices efficiently, yielding direct consequences for the Frobenius normal form, Kalman decomposition, and accelerated matrix exponentiation under suitable log-constraints. The approach blends kernel bases, Hermite forms, and truncated polynomial-matrix operations to surpass Keller-Gehrig’s prior bound and to enable deterministic improvements for classical linear-algebra tasks that are central in control theory and symbolic computation.
Abstract
Krylov methods rely on iterated matrix-vector products $A^k u_j$ for an $n\times n$ matrix $A$ and vectors $u_1,\ldots,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, Au_1, A^2u_1,\ldots$, 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^ω\log(n))$ comes from repeated squarings of $A$ and logarithmically many Gaussian eliminations. Here $ω>2$ is a feasible exponent for matrix multiplication over the base field. We present an algorithm computing the maximal Krylov basis in $O(n^ω\log\log(n))$ field operations when $m \in O(n)$, and even $O(n^ω)$ 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^ω(\log\log(n))^2)$, and therefore matrix exponentiation~$A^k$ can be performed in the latter complexity if $\log(k) \in O(n^{ω-1-\varepsilon})$ for some fixed $\varepsilon>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.
