Computing Krylov iterates in the time of matrix multiplication

Krylov methods rely on iterated matrix-vector products Akuj for an n × n matrix A and vectors u1, …, um. The space spanned by all iterates Akuj admits a particular basis — the maximal Krylov basis — which consists of iterates of the first vector u1, Au1, A2u1, …, 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 <Formula format="inline"><TexMath><?TeX $\mathchoice{O\left(n^\omega \log (n)\right)}{O(n^\omega \log (n))}{O(n^\omega \log (n))}{O(n^\omega \log (n))}$?></TexMath><AltText>Math 1</AltText><File name="issac24-45-inline1" type="svg"/></Formula> 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 <Formula format="inline"><TexMath><?TeX $\mathchoice{O\left(n^\omega \operatorname{loglog}(n)\right)}{O(n^\omega \operatorname{loglog}(n))}{O(n^\omega \operatorname{loglog}(n))}{O(n^\omega \operatorname{loglog}(n))}$?></TexMath><AltText>Math 2</AltText><File name="issac24-45-inline2" type="svg"/></Formula> field operations when <Formula format="inline"><TexMath><?TeX $m \in \mathchoice{O\left(n\right)}{O(n)}{O(n)}{O(n)}$?></TexMath><AltText>Math 3</AltText><File name="issac24-45-inline3" type="svg"/></Formula>, and even <Formula format="inline"><TexMath><?TeX $\mathchoice{O\left(n^\omega \right)}{O(n^\omega)}{O(n^\omega)}{O(n^\omega)}$?></TexMath><AltText>Math 4</AltText><File name="issac24-45-inline4" type="svg"/></Formula> as soon as <Formula format="inline"><TexMath><?TeX $m\in \mathchoice{O\left(n/\log (n)^c\right)}{O(n/\log (n)^c)}{O(n/\log (n)^c)}{O(n/\log (n)^c)}$?></TexMath><AltText>Math 5</AltText><File name="issac24-45-inline5" type="svg"/></Formula> 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 <Formula format="inline"><TexMath><?TeX $\mathchoice{O\left(n^\omega (\operatorname{loglog}(n))^2\right)}{O(n^\omega (\operatorname{loglog}(n))^2)}{O(n^\omega (\operatorname{loglog}(n))^2)}{O(n^\omega (\operatorname{loglog}(n))^2)}$?></TexMath><AltText>Math 6</AltText><File name="issac24-45-inline6" type="svg"/></Formula>, and therefore matrix exponentiation Ak can be performed in the latter complexity if <Formula format="inline"><TexMath><?TeX $\log (k) \in \mathchoice{O\left(n^{\omega -1-\varepsilon }\right)}{O(n^{\omega -1-\varepsilon })}{O(n^{\omega -1-\varepsilon })}{O(n^{\omega -1-\varepsilon })}$?></TexMath><AltText>Math 7</AltText><File name="issac24-45-inline7" type="svg"/></Formula> for some fixed ε > 0. A key idea for these improvements is to rely on fast algorithms for m × m polynomial matrices of average degree n/m, involving high-order lifting and minimal kernel bases.