  We present a sublinear query algorithm for outputting a near-optimal low-rank
approximation to any positive semidefinite Toeplitz matrix $T \in \mathbb{R}^{d
\times d}$. In particular, for any integer rank $k \leq d$ and $\epsilon,\delta
> 0$, our algorithm makes $\tilde{O} \left (k^2 \cdot \log(1/\delta) \cdot
\text{poly}(1/\epsilon) \right )$ queries to the entries of $T$ and outputs a
rank $\tilde{O} \left (k \cdot \log(1/\delta)/\epsilon\right )$ matrix
$\tilde{T} \in \mathbb{R}^{d \times d}$ such that $\| T - \tilde{T}\|_F \leq
(1+\epsilon) \cdot \|T-T_k\|_F + \delta \|T\|_F$. Here, $\|\cdot\|_F$ is the
Frobenius norm and $T_k$ is the optimal rank-$k$ approximation to $T$, given by
projection onto its top $k$ eigenvectors. $\tilde{O}(\cdot)$ hides
$\text{polylog}(d) $ factors. Our algorithm is \emph{structure-preserving}, in
that the approximation $\tilde{T}$ 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 $\tilde{T}$ 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.
