5.3 Block Krylov Iteration

Subspace Iteration can dramatically improve on the performance of the Randomized SVD when $\vec{A}$ has a heavy tail. The key observation is that the singular value tail of $(\vec{A}\vec{A}^\T)^q\vec{A}$ is much lighter than that of $\vec{A}$. This is because the polynomial $x^{2q+1}$ pushes tail singular values of $\vec{A}$ to zero.

From a computational perspective, there isn’t that much special about $x^{2q+1}$. In fact, we can replace this monomial with any polynomial of the same degree containing only odd degree terms, without significantly changing the cost of the algorithm. Thus, we might consider an approximation of the form:

It’s conceivable that, if we choose a good polynomial, this method could be more accurate than the subspace iteration. For instance, Chebyshev polynomials grow substantially faster than monomials.

To see why, note that the approximation $\vec{Q} = \Call{orth}(p(\vec{A}\vec{A}^\T)\vec{A}\vec{\Omega})$ scales the sketch along the $i$-th singular direction by $\sigma_i\, p(\sigma_i^2)$, so we want $|p|$ to be large at the leading singular values and small on the tail. Rescaling so the tail (squared) fills $[0,1]$ and the singular values we want to capture sit just above 1, the ideal polynomial stays small on $[0,1]$ but shoots up just past it — so a faster-growing polynomial gives a better approximation at the same degree.

The figure below compares the monomial $x^d$ with the Chebyshev polynomial $T_d(x)$. Both equal 1 at $x=1$ and stay within $[-1,1]$ on the tail $[0,1]$, but for $x$ just above 1 the Chebyshev polynomial is far larger. The right panel fixes an evaluation point $x$ slightly above 1 and plots each polynomial’s value as the degree grows. Drag the sliders to change the degree $d$ and the evaluation point $x$.

Block Krylov Iteration¶

The difficulty with the approach in (5.7) is that we don’t know how to choose a good polynomial in advance, since we don’t know the singular values of $\vec{A}$.

Observe that for any polynomial $p(x)$ with $\deg(p)<t$,

where

In particular, note that, using the same number of matrix-vector products with $\vec{A}$ and $\vec{A}^\T$ as required to compute $p(\vec{A}\vec{A}^\T)\vec{A}\vec{\Omega}$, we can actually compute an orthonormal basis for the entire block Krylov subspace $\mathcal{K}_t(\vec{A}\vec{A}^\T,\vec{A}\vec{\Omega})$.

This suggests the Randomized Block Krylov Iteration (RBKI) approximation:

Intuitively, by projecting onto a larger space, we can get a better approximation.

Block-size versus depth¶

The orthonormal basis for $\mathcal{K}_t(\vec{A}\vec{A}^\T,\vec{A}\vec{\Omega})$ has up to $tb$ columns, where $b$ is the number of columns in $\vec{\Omega}$, so we may hope to obtain a good rank-$k$ approximation even if $b<k$. However, this raises the practical question of how the block-size should be set. We tend to observe two competing factors:

• Smaller is better: While the worst-case matrix-vector product complexity for block sizes $b=1$ and $b=k$ match, the convergence of Krylov methods on a given instance $\vec{A}$ is highly dependent on $\vec{A}$'s spectrum, so “better-than-worst-case” performance is often observed. Theoretical evidence suggests that, in the absence of small singular value gaps, choosing $b=1$ typically requires fewer matrix-vector products to reach a desired level of accuracy for non-worst-case instances; see Meyer et al., 2024.

• Bigger is better: Due to parallelism, caching, and dedicated hardware accelerators in modern computing systems, it is typically much faster to compute multiple matrix-vector products all at once (i.e., grouped into a matrix-matrix product) instead of one after the other. This was observed in our experiments on the cost of NLA. RBKI groups matrix-vector products into blocks of size $b$, so for a fixed number of total matrix-vector products, we expect the method to run faster if $b$ is larger. The benefits of blocking into matrix-matrix products can be even more pronounced when $\vec{A}$ is too large to store in fast memory, in which case the cost of reading $\vec{A}$ from slow memory is a bottleneck, and using a larger block size $b\gg 1$ can be essentially “free”.

These two effects are in competition with one another, and one observes (e.g. as in the experiment below) that the best block size $b$ is often $1\ll b\ll k$.

Convergence Guarantees¶

Theoretical work answering this question is provided in Chen et al., 2025, which builds off of Meyer et al., 2024. In particular, Chen et al., 2025 asserts that, up to logarithmic factors in $n$ and certain spectral properties of $\vec{A}$, the number of matrix-vector products required to obtain a near-optimal approximation can be bounded independent of the block-size.

This means that users are free to choose the block-size based on their computational environments.

Numerical Experiment¶

We plot the value of $\varepsilon$ so that $\|\vec{A} -\llbracket \widehat{\vec{A}}\rrbracket_k \|_\F \leq (1+\varepsilon) \|\vec{A} - \llbracket \vec{A} \rrbracket_k\|_\F$.