Hankel matrices are an important class of highly-structured matrices, arising across computational mathematics, engineering, and theoretical computer science. It is well-known that positive semidefinite (PSD) Hankel matrices are always approximately low-rank. In particular, a celebrated result of Beckermann and Townsend shows that, for any PSD Hankel matrix \(H \in \mathbb{R}^{n \times n}\) and any \(\epsilon \gt 0\), letting \(H_k\) be the best rank-\(k\) approximation of \(H\) (obtained via truncated singular value decomposition), \(\|H - H_k\|_F \le \epsilon \|H\|_F\) for \(k = O(\log n \log (1/\epsilon))\). I.e., the optimal low-rank approximation error decays exponentially in the rank-\(k\). As such, PSD Hankel matrices are natural targets for low-rank approximation algorithms. We give the first such algorithm that runs in sublinear time. In particular, we show how to compute, in \(\operatorname{polylog}(n, 1/\epsilon)\) time, a factored representation of a rank-\(O(\log n \log(1/\epsilon))\) Hankel matrix \(\widehat{H}\) matching the error guarantee of Beckermann and Townsend up to constant factors. We further show that our algorithm is robust – given input \(H + E\) where \(E \in \mathbb{R}^{n \times n}\) is an arbitrary non-Hankel noise matrix, we obtain error \(\|H - \widehat{H}\|_F \le O(\|E\|_F) + \epsilon \|H\|_F\). Towards this algorithmic result, our first contribution is a structure-preserving existence result – we show that there exists a rank-\(k\) Hankel approximation to \(H\) matching the error bound of Beckermann and Townsend. Our result can be interpreted as a finite-dimensional analog of the widely applicable AAK theorem, which shows that the optimal low-rank approximation of an infinite Hankel operator is itself Hankel. Armed with our existence result, and leveraging the well-known Vandermonde structure of Hankel matrices, we achieve our sublinear time algorithm using a sampling-based approach that relies on universal ridge leverage score bounds for Vandermonde matrices.

Motivated by the problem of fast processing of attention matrices, we study fast algorithms for computing matrix-vector products for asymmetric Gaussian Kernel matrices $K\in \mathbb{R}^{n\times n}$. $K$'s columns are indexed by a set of $n$ keys $k_1,k_2\ldots, k_n\in \mathbb{R}^d$, rows by a set of $n$ queries $q_1,q_2,\ldots,q_n\in \mathbb{R}^d $, and its $i,j$ entry is $K_{ij} = e^{-\|q_i-k_j\|_2^2/2\sigma^2}$ for some bandwidth parameter $\sigma>0$. Given a vector $x\in \mathbb{R}^n$ and error parameter $\epsilon>0$, our task is to output a $y\in \mathbb{R}^n$ such that $\|Kx-y\|_2\leq \epsilon \|x\|_2$ in time subquadratic in $n$ and linear in $d$. Our algorithms rely on the following modelling assumption about the matrices $K$: the sum of the entries of $K$ scales linearly in $n$, as opposed to worst case quadratic growth. We validate this assumption experimentally, for Gaussian kernel matrices encountered in various settings such as fast attention computation in LLMs. Under this assumption, we obtain the first subquadratic time algorithm for kernel matrix-vector multiplication for unrestricted vectors.

Large language models (LLMs) have achieved impressive success, but their high memory requirements present challenges for long-context token generation. In this paper we study the streaming complexity of attention approximation, a key computational primitive underlying token generation. Our main contribution is BalanceKV, a streaming algorithm for $\epsilon$-approximating attention computations based on geometric process for selecting a balanced collection of Key and Value tokens as per Banaszczyk's vector balancing theory. We complement our algorithm with space lower bounds for streaming attention computation. Besides strong theoretical guarantees, BalanceKV exhibits empirically validated performance improvements over existing methods, both for attention approximation and end-to-end performance on various long context benchmarks.

We present a sublinear time algorithm for computing a near optimal low-rank approximation to any positive semidefinite (PSD) Toeplitz matrix T ∈ ℝd×d, given noisy access to its entries. In particular, given entrywise query access to T + E for an arbitrary noise matrix E ∈ ℝd×d, integer rank k≤d, and error parameter δ > 0, our algorithm runs in time poly(k, log(d/δ)) and outputs (in factored form) a Toeplitz matrix with rank poly(k, log(d/δ)) satisfying, for some fixed constant C,

We present a sublinear query algorithm for outputting a near-optimal low-rank approximation to any positive semidefinite Toeplitz matrix T ∈ ℝd×d. In particular, for any integer rank k ≤ d and ε, δ > 0, our algorithm makes Õ (k2 · log(1/δ) · poly(1/ε)) queries to the entries of T and outputs a rank Õ (k · log(1/δ)/ε) matrix d×d such that ||T – ||F ≤ (1 + ε) · ||T - Tk ||F + δ||Τ||F. Here, || · ||F is the Frobenius norm and Tk is the optimal rank-k approximation to T, given by projection onto its top k eigenvectors. Õ(·) hides polylog(d) factors. Our algorithm is structure-preserving, in that the approximation is also Toeplitz. A key technical contribution is a proof that any positive semidefinite Toeplitz matrix in fact has a near-optimal low-rank approximation which is itself Toeplitz. Surprisingly, this basic existence result was not previously known. Building on this result, along with the well-established off-grid Fourier structure of Toeplitz matrices [Cybenko'82], we show that Toeplitz with near optimal error can be recovered with a small number of random queries via a leverage-score-based off-grid sparse Fourier sampling scheme.