Sample Efficient Toeplitz Covariance Estimation

TL;DR

This paper develops new sample-efficient algorithms for estimating Toeplitz covariance matrices, significantly reducing sample and entry access requirements, especially when the covariance is low-rank.

Contribution

It introduces novel bounds and algorithms for Toeplitz covariance estimation, leveraging structure and low-rank properties to improve efficiency over previous methods.

Findings

Non-asymptotic bounds on sample complexity for Toeplitz covariance estimation.

Algorithms based on sparse rulers outperform generic methods, especially for low-rank matrices.

Sample complexity scales sublinearly with dimension in low-rank cases.

Abstract

We study the sample complexity of estimating the covariance matrix $T$ of a distribution $\mathcal{D}$ over $d$-dimensional vectors, under the assumption that $T$ is Toeplitz. This assumption arises in many signal processing problems, where the covariance between any two measurements only depends on the time or distance between those measurements. We are interested in estimation strategies that may choose to view only a subset of entries in each vector sample $x \sim \mathcal{D}$, which often equates to reducing hardware and communication requirements in applications ranging from wireless signal processing to advanced imaging. Our goal is to minimize both 1) the number of vector samples drawn from $\mathcal{D}$ and 2) the number of entries accessed in each sample. We provide some of the first non-asymptotic bounds on these sample complexity measures that exploit $T$'s Toeplitz structure, and by doing so, significantly improve on results for generic covariance matrices. Our bounds follow from a novel analysis of classical and widely used estimation algorithms (along with some new variants), including methods based on selecting entries from each vector sample according to a so-called sparse ruler. In many cases, we pair our upper bounds with matching or nearly matching lower bounds. In addition to results that hold for any Toeplitz $T$, we further study the important setting when $T$ is close to low-rank, which is often the case in practice. We show that methods based on sparse rulers perform even better in this setting, with sample complexity scaling sublinearly in $d$. Motivated by this finding, we develop a new covariance estimation strategy that further improves on all existing methods in the low-rank case: when $T$ is rank-$k$ or nearly rank-$k$, it achieves sample complexity depending polynomially on $k$ and only logarithmically on $d$.

Tables

Table 1. Table 1 : Overview of our results on estimating a d 𝑑 d -dimensional Toeplitz covariance from samples to spectral norm error ε ​ ‖ T ‖ 2 𝜀 subscript norm 𝑇 2 \varepsilon\|T\|_{2} (Problem  2.1 ). The above results hold for success probability 2 / 3 2 3 2/3 and we let O ~ ​ ( ⋅ ) ~ 𝑂 ⋅ \widetilde{O}(\cdot) and Ω ~ ​ ( ⋅ ) ~ Ω ⋅ \widetilde{\Omega}(\cdot) hide poly ⁡ ( log ⁡ d , log ⁡ k , log ⁡ 1 / ε ) poly 𝑑 𝑘 1 𝜀 \operatorname{poly}(\log d,\log k,\log 1/\varepsilon) factors. See the corresponding theorems for full statements of the results.

Algorithm	ESC	VSC upper bound	TSC lower bound
general $d\times d$ Toeplitz matrix
Full samples	$d$	$\tilde{O}\left(1/{\epsilon }^2\right)$ (Thm. 2.4)	$\Omega \left(d/{\epsilon }^2\right)$ (folklore)
$\Theta \left(\sqrt{d}\right)$-sparse ruler	$\Theta \left(\sqrt{d}\right)$	$\tilde{O}\left(d/{\epsilon }^2\right)$ (Thm. 2.6)	$\Omega \left(d^{3/2}/{\epsilon }^2\right)$ (Thm. 4.5)
$\Theta \left(d^{\alpha }\right)$-sparse ruler, $\alpha \in [\frac{1}{2},1]$	$\Theta \left(d^{\alpha }\right)$	$\tilde{O}\left(d^{2-2\alpha }/{\epsilon }^2\right)$ (Thm. 4.4)	$\Omega \left(d^{3-3\alpha }/{\epsilon }^2\right)$ (Thm. 4.5)
rank-$k$ Toeplitz matrix
Prony’s method	$2k$	$\tilde{O}\left(1/{\epsilon }^2\right)$ (Thm. C.8)	See below
Any non-adaptive method	–	–	$\Omega \left(k/\log k\right)$ (Thm. 6.1)
approximately rank-$k$ Toeplitz matrix
$\Theta \left(\sqrt{d}\right)$-sparse ruler	$\Theta \left(\sqrt{d}\right)$	$\tilde{O}\left(\min (k^2,d)/{\epsilon }^2\right)$ (Thm. 2.8)	$\Omega \left(\sqrt{d}k/{ϵ}^2\right)$ (Thm. 4.5)
Algorithm 3 (sparse Fourier trans. based)	$\tilde{O}\left(k^2\right)$	$\tilde{O}\left(k+1/{\epsilon }^2\right)$ (Thm. 2.10)	See below
Any adaptive method	–	–	$\Omega \left(k\right)$ (Thm. 6.6)

Equations

$$\Vert T-\tilde{T}{\Vert }_2\le \epsilon \Vert T{\Vert }_2,$$

$$F(f_j{)}^{\top }=\left[1,e^{-2\pi i{f}_j},e^{-2\pi i2f_j},\dots ,e^{-2\pi i(d-1)f_j}\right].$$

$$\Vert \mathrm{Toep}(a){\Vert }_2$$

$${\tau }_j(A):=a_j{\left(A^{\ast }A\right)}^{+}{a}_j^{\ast }.$$

$${\tau }_j(A)=\underset{y\in C^d\text{ s.t. }{y}^{T}A=a_j}{\min }\Vert y{\Vert }_2^2.$$

$${\tau }_j(A)=\underset{y\in C^s}{\max }\frac{{\left|{\left(Ay\right)}_j\right|}^2}{\Vert Ay{\Vert }_2^2}.$$

$$d_{\mathrm{KL}}(N(0,M_1),N(0,M_2))=\Theta \left({\Vert I-M_2^{-1/2}{M}_1{M}_2^{-1/2}\Vert }_F^2\right).$$

$$d_{\mathrm{TV}}(N(0,M_1),N(0,M_2))=O\left({\Vert I-M_1^{-1/2}{M}_2{M}_1^{-1/2}\Vert }_F\right).$$

$$\mathrm{Pr}\left[\left|\frac{1}{n}\sum _{\ell =1}^n{x^{(\ell )}}^{T}A{x}^{(\ell )}-\mathrm{tr}(A)\right|>t\right]\le 2\exp \left(-c\cdot n\min \left(\frac{t^2}{\Vert A{\Vert }_F^2},\frac{t}{\Vert A{\Vert }_2}\right)\right).$$

$$R_{\alpha }^{(1)}=\{1,2,\dots ,d^{\alpha }\},\text{mbox}and{R}_{\alpha }^{(2)}=\{d,d-d^{1-\alpha },d-2d^{1-\alpha },\dots ,d-(d^{\alpha }-1)d^{1-\alpha }\}.$$

$$\Delta (R_{\alpha })\le 2d^{2-2\alpha }+d^{1-\alpha }(1+\log (\lceil d^{2\alpha -1}\rceil ))=d^{2-2\alpha }+O(d^{1-\alpha }\cdot \log d).$$

$$\mathrm{Pr}\left[{\Vert T-\tilde{T}\Vert }_2>\epsilon \Vert T{\Vert }_2\right]\le \frac{C{d}^2}{\epsilon }\exp \left(-cn\kappa \right).$$

$$\underset{e}{\mathrm{Pr}}\left[\exists x\in [0,1]:|L_e(x)|\ge \epsilon \Vert T{\Vert }_2\right]\le \frac{C{d}^2}{\epsilon }\exp \left(-cn\kappa \right).$$

$$\underset{e}{\mathrm{Pr}}\left[\exists s\in \{0,\dots ,d-1\}:|e_s|>10\Vert T{\Vert }_2\right]$$

$$\underset{e}{\mathrm{Pr}}\left[\exists x\in N:|L_e(x)|>\frac{\epsilon }{2}\Vert T{\Vert }_2\right]$$

$$1-\frac{4Q{d}^2}{{\epsilon }^2}\exp \left(-cn\kappa \right).$$

$$|L_e^{\prime }(x)|$$

$$\le 4\pi d^2\Vert e{\Vert }_{\infty }\le 40\pi d^2\Vert T{\Vert }_2,$$

$$|L(x)-L(x^{\prime })|$$

$$\le \frac{\epsilon }{2}\Vert T{\Vert }_2+40\pi d^2\Vert T{\Vert }_2\cdot |x-x^{\prime }|\le \epsilon \Vert T{\Vert }_2,$$

$$e_s=a_s-{\tilde{a}}_s=\frac{1}{n\cdot |R_s|}\sum _{\ell =1}^n\sum _{(j,k)\in R_s}\left[T_{j,k}-x_j^{(\ell )}\cdot x_k^{(\ell )}\right].$$

$$E\left[(x_j^{(\ell )}{)}^2\cdot (x_k^{(\ell )}{)}^2\right]=T_{j,j}{T}_{k,k}+2T_{j,k}^2\le 3\Vert T{\Vert }_2^2.$$

$$M_{j,k}=\frac{\cos (2\pi sx)}{|R_s|}.$$

$$L_e(x)=e_0+2\sum _{s=1}^{d-1}{e}_s\cos (2\pi sx)=\mathrm{tr}(T_R-{\tilde{T}}_R,M_R)=\mathrm{tr}(T_R-{\overline{T}}_R,M_R).$$

$$\sum _{(j,k)\in R\times R:|j-k|=s}{\tilde{T}}_{j,k}=\sum _{(j,k)\in R\times R:|j-k|=s}{\overline{T}}_{j,k},$$

$$\mathrm{tr}({\tilde{T}}_R,M_R)$$

$$=\sum _{s=0}^{d-1}\left[\sum _{(j,k)\in R\times R:|j-k|=s}{\tilde{T}}_{j,k}\cdot M_{j,k}\right]$$

$$=\sum _{s=0}^{d-1}\left[\sum _{(j,k)\in R\times R:|j-k|=s}{\overline{T}}_{j,k}\cdot M_{j,k}\right]=\mathrm{tr}({\overline{T}}_R,M).$$

$$\mathrm{Pr}\left[|\mathrm{tr}(T-{\overline{T}}_R,M)|>\frac{\epsilon }{2}\Vert T{\Vert }_2\right]$$

$$\mathrm{tr}(T_R-{\overline{T}}_R,M)$$

$$=\mathrm{tr}(M^{\prime })-\frac{1}{n}\sum _{\ell =1}^n{y^{(\ell )}}^T{M}^{\prime }{y}^{(\ell )},$$

$$\mathrm{Pr}\left[|\mathrm{tr}(T_R-{\overline{T}}_R,M)|>\frac{\epsilon }{2}\Vert T{\Vert }_2\right]$$

$$\le 2\exp \left(-cn\min \left(\frac{{\epsilon }^2\Vert T{\Vert }_2^2}{4\Vert M^{\prime }{\Vert }_F^2},\frac{\epsilon \Vert T{\Vert }_2}{2\Vert M^{\prime }{\Vert }_2}\right)\right)$$

$$\le 2\exp \left(-cn\min \left(\frac{\Vert T{\Vert }_2^2}{\Vert T_R{\Vert }_2^2}\cdot \frac{{\epsilon }^2}{4\Vert M{\Vert }_F^2},\frac{\Vert T{\Vert }_2}{\Vert T_R{\Vert }_2}\cdot \frac{\epsilon }{2\Vert M{\Vert }_2}\right)\right)$$

Full text

We provide some of the first non-asymptotic bounds on these sample complexity measures that exploit $T$’s Toeplitz structure, and by doing so, significantly improve on results for generic covariance matrices. These bounds follow from a novel analysis of classical and widely used estimation algorithms (along with some new variants), including methods based on selecting entries from each vector sample according to a so-called sparse ruler.

In addition to results that hold for any Toeplitz $T$, we further study the important setting when $T$ is close to low-rank, which is often the case in practice. We show that methods based on sparse rulers perform even better in this setting, with sample complexity scaling sublinearly in $d$. Motivated by this finding, we develop a new covariance estimation strategy that further improves on existing methods in the low-rank case: when $T$ is rank-$k$ or nearly rank-$k$, it achieves sample complexity depending polynomially on $k$ and only logarithmically on $d$.

Our results utilize tools from random matrix sketching, leverage score based sampling techniques for continuous time signals, and sparse Fourier transform methods. In many cases, we pair our upper bounds with matching or nearly matching lower bounds.

1 Introduction

Estimating the covariance matrix of a distribution $\mathcal{D}$ over vectors in $\operatorname{\mathbb{C}}^{d}$ given independent samples $x^{(1)},\ldots,x^{(n)}\sim\mathcal{D}$ is a fundamental statistical problem. In signal processing, many applications require this problem to be solved under the assumption that $\mathcal{D}$’s covariance matrix, $T\in\operatorname{\mathbb{R}}^{d\times d}$, is a symmetric Toeplitz matrix. Specifically, $T_{a,b}=T_{c,d}$ whenever $|a-b|=|c-d|$.

Toeplitz structure naturally arises when entries in the random vector correspond to measurements on a spatial or temporal grid, and the covariance between measurements only depends on the distance between them, as in a stationary process. To name just a few applications (see [RATL16] for more), estimation algorithms for Toeplitz covariance matrices are used in:

• •

Direction-of-arrival (DOA) estimation for signals received by antenna arrays [KV96], which allows, for example, cellular networks to identify and target signal transmissions based on the geographic location of devices [DMP*+*14, BL16].

• •

Spectrum sensing for cognitive radio [MLJ09, CTE18], which allows a receiver to estimate which parts of the frequency spectrum are in use at any given time.

• •

Medical and radar imaging processing [SOM89, RSTL88, Fuh91, BVR*+*08, AM12, CE18].

In these applications, the goal is to find an approximation to $T$ using as few samples as possible. In contrast to generic covariance estimation, there is significant interest in algorithms that consider only a subset of the entries in each of the samples $x^{(1)},\ldots,x^{(n)}$. This subset can be chosen in any way, deterministically or randomly. This suggests two separate measures of sample complexity:

1. Vector sample complexity (VSC).

How many vector samples $x^{(1)},\ldots,x^{(n)}\sim\mathcal{D}$ does an algorithm require to estimate the Toeplitz covariance matrix $T$ up to a specified tolerance?

2. Entry sample complexity (ESC).

How many entries, $s$, does the algorithm view in each $x^{(\ell)}$?

There is typically a trade off between VSC and ESC, which makes different algorithms suited to different applications. For example, in direction-of-arrival estimation, fewer vector samples means a shorter acquisition time, while fewer entry samples means that fewer active receiving antennas are needed. An algorithm with high ESC and low VSC might give optimal performance in terms of acquisition speed, while one with low ESC and higher VSC would minimize hardware costs. In many cases, it is also natural to simply consider the combination of VSC and ESC:

3. Total sample complexity (TSC).

How many total entries $n\cdot s$ does the algorithm read across all sampled vectors to estimate the covariance $T$?

In this paper we introduce new estimation algorithms and analyze some common practical approaches, to obtain non-asymptotic bounds (sometimes tight) on these fundamental sample complexity measures. Surprisingly, relatively little was known previously in the non-asymptotic setting, beyond results for general covariance matrices without Toeplitz structure. Intuitively, since Toeplitz matrices are parameterized by just $O(d)$ variables (one for each diagonal) we might expect to learn them with fewer samples.

1.1 Our contributions in brief

We first prove that it is in fact possible to estimate Toeplitz covariance matrices with lower vector sample complexity than general covariance matrices. We focus on the common case when $x^{(1)},\ldots,x^{(n)}$ are drawn from a $d$-dimensional Gaussian distribution with covariance $T$ and consider the goal of recovering an approximation $\tilde{T}$ with spectral norm error $\|T-\tilde{T}\|_{2}\leq\varepsilon\|T\|_{2}$. We show that with entry sample complexity (ESC) $s=d$ (i.e., we view all entries in each vector sample) a fast and simple algorithm222Compute the sample covariance of $x^{(1)},\ldots,x^{(n)}$ and average its diagonals to make it Toeplitz. achieves vector sample complexity (VSC) $O\left(\frac{\log(d/\varepsilon)\cdot\log d}{\varepsilon^{2}}\right)$.

This bound is significantly less than the $\Omega\left({d}/{\varepsilon^{2}}\right)$ required for generic covariance matrices, and serves as a baseline for the rest of our results, which focus on methods with ESC $s<d$. In this setting, we give a number of new upper and lower bounds, which we review in detail, along with our main techniques, in Section 2. We briefly describe our contributions here and summarize the bounds in Table 1. Several of our results apply in the important practical setting (see e.g., [KV96, CEC13, QP17]) where $T$ is both Toeplitz and either exactly or approximately rank-$k$ for some $k\ll d$ .

First, in Section 4, we give a non-asymptotic analysis of a widely used estimation scheme that reads $\Theta(\sqrt{d})$ entries per vector according to a sparse ruler [Mof68, RATL16]. Sparse ruler methods are important because, up to constant factors, they minimize ESC among all methods that read a fixed subset of entries from each vector sample. This is a natural constraint because, in many applications, the subset of entries read is encoded in the signal acquisition hardware. For algorithms that read a fixed subset of entries, it is easy to exhibit Toeplitz covariance matrices that can only be estimated if this subset has size $\Omega(\sqrt{d})$.

While optimal in terms of ESC, at least without additional assumptions on $T$, the VSC of sparse ruler methods was not previously well understood.We prove a VSC bound that is nearly linear in $d$, which we show is tight up to logarithmic factors. Moreover, we introduce a new class of ruler based methods that read anywhere between $\Theta(\sqrt{d})$ and $d$ entries per vector, smoothly interpolating between the ideal ESC of standard sparse rulers and the ideal $O(\log^{2}d)$ VSC of fully dense sampling.

Beyond these results, we theoretically confirm a practical observation: sparse ruler methods perform especially well when $T$ is low-rank or close to low-rank. The VSC of these methods decreases with rank, ultimately allowing for total sample complexity sublinear in $d$. Inspired by this finding, Section 5 is devoted to developing methods specifically designed for estimating low-rank, or nearly low-rank, Toeplitz covariance matrices. We develop algorithms based on combining classical tools from harmonic analysis (e.g., the Carathéodory-Fejér-Pisarenko decomposition of Toeplitz matrices) with methods from randomized numerical linear algebra [Woo14, DM16].

In particular, our work builds on connections between random matrix sketching, leverage score based sampling, and sampling techniques for continuous Fourier signals, which have found a number of recent applications [CKPS16, CP18, CP19, AKM*+*17, AKM*+*19]. Ultimately, we develop an estimation algorithm with total sample complexity depending just logarithmically on $d$ and polynomially on $k$ for any $T$ that is approximately rank-$k$. The method reads $\tilde{O}(k^{2})$ entries from each $x^{(\ell)}$ using a fixed pattern, which is constructed randomly. To the best of our knowledge, this approach is the first general Toeplitz covariance estimation algorithm with ESC $<\sqrt{d}$. It provides a potentially powerful alternative to sparse ruler based methods, so tightening our results (potentially to linear in $k$) is an exciting direction for future work. In Section 6 we demonstrate that a linear dependence on $k$ is unavoidable for total sample complexity, even for adaptive sampling algorithms.

1.2 Comparison to prior work

An in-depth discussion of prior work is included in Section 2. Our results on general full-rank Toeplitz covariance matrices give much tighter bounds on VSC than prior work. While some results address the variance of estimators (including those based on sparse rulers) for each entry in $T$ [AL12], the challenge is that obtaining a bound on $\|T-\tilde{T}\|_{2}$ requires understanding potentially strong correlations between the errors on different matrix entries. A naive analysis leads to a VSC bound of $\Omega(d^{2}/\epsilon^{2})$ for both full samples and sparse ruler samples, which is significantly worse than our respective bounds of $\tilde{O}(1/\epsilon^{2})$ and $\tilde{O}(d/\epsilon^{2})$ proven in Theorems 2.4 and 2.6. Most closely related to our work is that of Qiao and Pal [QP17]. While they consider a different error metric, their analysis yields VSC $O(d/\epsilon^{2})$ for sparse ruler samples to estimate $T$ to relative error in the Frobenius norm. In this case, $O(d/\epsilon^{2})$ samples can be obtained without considering correlations between the errors on different matrix entries. Qiao and Pal also consider the case when $T$ is approximately low-rank, and like us show that an even lower ESC is possible in this setting. They actually give ESC $O(\sqrt{k}),$ which goes beyond any of our bounds, however their analysis requires strong assumptions on the Carathéodory-Fejér-Pisarenko decomposition of $T$.

For low-rank or nearly low-rank covariance matrices, our total sample complexity bounds depend polynomially in the rank $k$ and sublinearly on $d$, in fact just logarithmically in Theorems C.8 and 2.10. When $T$ is not assumed to be Toeplitz, TSC must depend linearly in $d$, even when $T$ is low-rank. This dependence is reflected in prior work, which obtains TSC bounds of at best $O(dk/\epsilon^{2})$, for any ESC between $2$ and $d$ [GRESS16]. Our work critically takes advantage of Toeplitz and low-rank structure simultaneously to surpass such bounds.

2 Discussion of Results

Formally, we study sampling schemes and recovery algorithms for the following problem:

Problem 2.1** (Covariance estimation).**

For a positive semidefinite matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, given query access to the entries of i.i.d. samples drawn from a $d$-dimension normal distribution with covariance $T$, i.e., $x^{(1)},\ldots,x^{(n)}\sim\mathcal{N}(0,T)$, return $\tilde{T}$ satisfying, with probability $\geq 1-\delta$:

[TABLE]

*where $\lVert\cdot\rVert_{2}$ denotes the operator (spectral) norm. The total sample complexity (TSC) * of computing $\tilde{T}$ is the total number of entries that the algorithm must read from $x^{(1)},\ldots,x^{(n)}$, combined.

As discussed, in addition to total sample complexity, we would like to understand the tradeoff between the number of vector samples used by an algorithm and the maximum number of entries viewed in each vector individually (VSC and ESC, respectively). Problem 2.1 assumes a normal distribution for simplicity, but our results also hold when the $x^{(j)}$ are distributed as $x^{(j)}\sim T^{1/2}y^{(j)}$, where $y^{(j)}$ is an isotropic sub-gaussian random variable. We also conjecture that similar bounds hold for more general classes of random variables. Finally, we note that when $T$ is Toeplitz, our algorithms will return $\tilde{T}$ that is also Toeplitz, which is useful in many applications.

2.1 What is known

Standard matrix concentration results can be used to bound the error of estimating a generic covariance matrix $T$ when $\tilde{T}$ is set to the empirical covariance $\frac{1}{n}\sum_{j=1}^{n}x^{(j)}{x^{(j)}}^{T}$ [Ver12]. Computing this estimate requires reading all $d$ entries of each sample, yielding the following standard bound:

Claim 2.2** (General covariance estimation).**

For any positive semidefinite $T\in\operatorname{\mathbb{R}}^{d\times d}$, Problem 2.1 can be solved with vector sample complexity $O\left(\frac{d+\log(1/\delta)}{\varepsilon^{2}}\right)$, entry sample complexity $d$, and total sample complexity $O\left(\frac{d^{2}+d\log(1/\delta)}{\varepsilon^{2}}\right)$. These bounds are achieved by setting $\tilde{T}=\frac{1}{n}\sum_{j=1}^{n}x^{(j)}{x^{(j)}}^{T}$.

Claim 2.2 applies for generic covariance matrices, which are specified by $\Theta(d^{2})$ parameters instead of the $\Theta(d)$ required for Toeplitz matrices. For other matrices specified by $\Theta(d)$ parameters, it is possible to improve on Claim 2.2. When $T$ is diagonal, the entries of a sample $x^{(j)}\sim\mathcal{N}(0,T)$ are independent Gaussians with variances equal to $T$’s diagonal. So each diagonal entry can be approximated to relative error $1\pm\varepsilon$ with probability $\geq 1-\delta/d$ using $O\left({\log(d/\delta)}/{\varepsilon^{2}}\right)$ samples. Applying a union bound, Problem 2.1 can be solved with $O\left({\log(d/\delta)}/{\varepsilon^{2}}\right)$ full vector samples.

Closer to our setting, it is possible to prove the same bound for any circulant covariance matrix: i.e., a Toeplitz matrix where each row is a cyclic permutation of the first. To see why, we use the fact that any circulant matrix, $T$, can be written as $T=FDF^{*}$ where $F\in\operatorname{\mathbb{C}}^{d\times d}$ is the DFT matrix and $D$ is diagonal. If we transform a sample $x\sim\mathcal{N}(0,T)$ by forming $F^{*}x$, we obtain a random vector with diagonal covariance, allowing us to apply the diagonal result. Overall we have:

Claim 2.3** (Diagonal and circulant covariance estimation).**

For any diagonal or circulant positive semidefinite $T\in\operatorname{\mathbb{R}}^{d\times d}$, Problem 2.1 can be solved with vector sample complexity $O\left(\frac{\log(d/\delta)}{\varepsilon^{2}}\right)$, entry sample complexity $d$, and total sample complexity $O\left(\frac{d\log(d/\delta)}{\varepsilon^{2}}\right)$. These bounds are achieved by setting $\tilde{T}=\mathrm{diag}\left(\frac{1}{n}\sum_{j=1}^{n}x^{(j)}{x^{(j)}}^{T}\right)$ when $T$ is diagonal, or $\tilde{T}=F\mathrm{diag}\left(\frac{1}{n}\sum_{j=1}^{n}F^{*}x^{(j)}{x^{(j)}}^{T}F\right)F^{*}$ when $T$ is circulant. Here $\mathrm{diag}(A)$ returns the matrix $A$ but with all off-diagonal entries set to 0.

See Appendix C for a more formal proof of this known result. Unfortunately, Claim 2.3 does not extend to general Toeplitz covariance matrices. While all Toeplitz matrices have some Fourier structure (as will be discussed later, they can be diagonalized by Fourier matrices with “off-grid” frequencies), they cannot be diagonalized using a single known basis like circulant matrices.

2.2 What is new: full rank matrices

Our first result is that Claim 2.3 can nevertheless be matched for Toeplitz matrices using a very simple algorithm: we compute the empirical covariance matrix $\overline{T}$ of $x^{(1)},\ldots,x^{(n)}$ and then form $\tilde{T}$ by averaging its diagonals to obtain a Toeplitz matrix (see Algorithm 1). Surprisingly, unlike the method for circulant matrices, this algorithm does not explicitly use the Fourier structure of the Toeplitz matrix, but achieves a similar sample complexity, up to logarithmic factors:

Theorem 2.4** (Near linear sample complexity).**

For any positive semidefinite Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, Problem 2.1 can be solved by Algorithm 1 with vector sample complexity $O\left(\frac{\log d\log\left(d/\varepsilon\delta\right)}{\varepsilon^{2}}\right)$, entry sample complexity $d$, and total sample complexity $O\left(\frac{d\log d\log\left(d/\varepsilon\delta\right)}{\varepsilon^{2}}\right)$.

The averaging method used for Theorem 2.4 has been suggested before as a straightforward way to improve a naive sample covariance estimator [CRZ13]. Our theorem shows that this improvement can be very significant. For example, when $T$ is the identity, the sample covariance truly requires ${\Omega}(d)$ vector samples to converge. Averaging reduces this complexity to $O(\log^{2}d)$.

Moreover, averaging gives a natural strategy for solving the Toeplitz problem with entry sample complexity $o(d)$. In particular, dating back to the work of [Mof68], who sought to minimize the number of receivers in directional antenna arrays, signal processing applications have widely applied the averaging method with less than $d$ samples taken from $x\sim\mathcal{N}(0,T)$ according to a so called sparse ruler [RATL16].

The idea is elegant and ingenious: a symmetric Toeplitz matrix $T$ is fully described by a vector $a\in\operatorname{\mathbb{R}}^{d}$ where $a_{i}$ is the value on the $i^{\text{th}}$ diagonal of $T$, i.e., the covariance of samples spaced $i$ steps apart. So, to estimate each entry in $a$, we only need to view a subset of entries $R\subseteq[d]$ from $x\sim\mathcal{N}(0,T)$ such that there are two entries, $u,v$ in $R$ at distance $s$ apart, for every distance $s\in 0,\ldots,d-1$. If we average $x^{(\ell)}_{u}x^{(\ell)}_{v}$ for $\ell\in[n]$, we will eventually converge on an estimate of $a_{s}$. It was noticed by [Mof68] and others [PBNH85] that such a subset $R$ can always be found with $|R|=\Theta(\sqrt{d})$. Formally we define:

Definition 2.1** (Ruler).**

A subset $R\subseteq[d]$ is a ruler if for all $s=0,\ldots,d-1$, there exist $j,k\in R$ so that $s=|j-k|$. We let $R_{s}\vcentcolon=\{(j,k)\in R\times R:|j-k|=s\}$ denote the set of (ordered) pairs in $R\times R$ with distance $s$. We say $R$ is “sparse” if $|R|<d$.

Sparse rulers are also called “sparse linear arrays” in the literature [WZY16] and are closely related to other constructions like Golumb rulers, which are also used in signal processing applications [Bab53]. It is known that any sparse ruler must have $|R|\gtrsim 1.557\sqrt{d}$ [EG48, Lee56], and it is possible to nearly match this bound for all $d$ with a linear time constructible ruler:

Claim 2.5** ($\Theta(\sqrt{d})$ Sparse Ruler).**

For any $d$, there is an explicit ruler $R$ of size $|R|=2\lceil\sqrt{d}\rceil-1$.

Proof.

It suffices to take $R=\{1,\ldots,\lceil\sqrt{d}\rceil\}\cup\{d,d-\lceil\sqrt{d}\rceil,d-2\lceil\sqrt{d}\rceil,\ldots,d-(\sqrt{d}-2)\lceil\sqrt{d}\rceil\}$. ∎

Using samples taken from $x^{(1)},\ldots,x^{(\ell)}$ according to a fixed $\Theta(\sqrt{d})$ sparse ruler and averaging to estimate $a$ (formalized in Algorithm 1) gives a $\tilde{T}$ that clearly converges to $T$ as $n\rightarrow\infty$. However, while commonly used for its potential to save on hardware cost, power, and space in applications, it has been unknown how much the sparse ruler strategy sacrifices in terms of vector sample complexity compared to looking at all entries in $x^{(1)},\ldots,x^{(\ell)}$. We bound this complexity:

Theorem 2.6** (Sparse ruler sample complexity).**

For any positive semidefinite Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, Problem 2.1 can be solved by Algorithm 1 using any $\Theta(\sqrt{d})$ sparse ruler (e.g., the construction from Claim 2.5) with vector sample complexity $O\left(\frac{d\log\left(d/\varepsilon\delta\right)}{\varepsilon^{2}}\right)$, entry sample complexity $\Theta(\sqrt{d})$, and total sample complexity $O\left(\frac{d^{3/2}\log\left(d/\varepsilon\delta\right)}{\varepsilon^{2}}\right)$.

More generally, in Section 4 we introduce a new class of rulers with sparsity $\Theta(d^{\alpha})$ for $\alpha\in[1/2,1]$. These rulers give an easy way to trade off entry sample complexity and vector sample complexity, with bounds interpolating smoothly between Theorem 2.4 and Theorem 2.6.

Note that the overall sample complexity of Theorem 2.6 is worse than that of Theorem 2.4 by roughly a $\sqrt{d}$ factor. This is in fact inherent: we provide a lower bound in Theorem 4.5.

Theorem 2.7** (Informal, see Theorem 4.5).**

For any $\varepsilon>0$ sufficiently small, any algorithm that takes samples from a ruler $R$ with sparsity $|R|=O(\sqrt{d})$ and solves Problem 2.1 with success probability $1-\delta\geq 1/10$ requires vector sample complexity $\Omega(d/\varepsilon^{2})$.

Thus, in terms of total sample complexity, the ruler-based algorithms are worse than taking full samples. However, the “hard” case for the sparse ruler methods appears to be when $T$ is close to an identity matrix – i.e., it is nearly rank $d$. In practice, it is much more common for $T$ to be nearly $k$ rank for some $k\ll d$. For example, this is the case in DOA estimation when an antenna array is detecting the direction of just $k$ different signal sources [KV96, SAL12, CEC13]. Experimentally, the performance of sparse ruler methods suffers less when $T$ is low-rank. In fact, in terms of TSC, they often significantly outperform algorithms that look at all entries in each $x^{(\ell)}$ (see Figure 3).

2.3 What is new: low-rank matrices

We confirm these experimental findings theoretically by proving a tighter bound for the sparse ruler estimation algorithm when $T$ is low-rank (or nearly low-rank). In terms of total sample complexity, our bound surpasses the algorithm that looks at all $d$ entries in each $x^{(\ell)}$. In fact, for $k\lesssim d^{1/4}$ we show that sublinear total sample complexity is possible. This isn’t even possible for diagonal covariance matrices that are close to low-rank.333To see that this is the case, consider diagonal $T$ with just a single entry placed in a random position equal to $1$ and the rest equal to [math]. Then $\Omega(d)$ samples are necessary to identify this entry, even though $T$ is rank $1$. Overall we obtain:

Theorem 2.8** (Sublinear sparse ruler sample complexity).**

For any positive semidefinite Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$ with rank $k$, Problem 2.1 can be solved by Algorithm 1 using the $\Theta(\sqrt{d})$ sparse ruler from Claim 2.5 with vector sample complexity $O\left(\frac{\min(k^{2},d)\log\left({d}/{\varepsilon\delta}\right)}{\varepsilon^{2}}\right)$, entry sample complexity $\Theta(\sqrt{d})$, and total sample complexity $O\left(\frac{\sqrt{d}\min(k^{2},d)\log\left({d}/{\varepsilon\delta}\right)}{\varepsilon^{2}}\right)$.

Importantly, Theorem 2.8 also holds when $T$ is close to rank-$k$, which is the typical case in practice. For example, if $T_{k}$ is $T$’s optimal $k$-rank approximation given by the singular value decomposition, the bound stated here also holds when $\frac{\sqrt{d}}{k^{2}}\cdot\lVert T-T_{k}\rVert_{F}^{2}\leq\lVert T\rVert_{2}^{2}$ or $\frac{d}{k^{2}}\cdot\lVert T-T_{k}\rVert_{2}^{2}\leq\lVert T\rVert_{2}^{2}$. See Section 4 for a formal statement and extensions to rulers with sparsity between $\Theta(\sqrt{d})$ and $d$. By achieving total sample complexity sublinear in $d$, Theorem 2.8 raises an intriguing question.

Can we learn a low-rank Toeplitz matrix with sample complexity independent of $d$?

In particular, the $O(\sqrt{d})$ dependence on the total sample complexity in Theorem 2.8 arises from the sparsity of the ruler. While sparse rulers with less than $\sqrt{d}$ elements do not exist, we might hope to design alternative approaches for selecting entries from $x^{(1)},\ldots,x^{(n)}$ which are specifically tailored to the low-rank setting. Doing so is the focus of our remaining results and requires more explicitly leveraging the Fourier structure in Toeplitz covariance matrices.

Specifically, the powerful Vandermonde decomposition theorem [CF11, Pis73] (see Fact 3.1) implies a connection between low-rank Toeplitz matrices and Fourier sparse functions: $T$ is rank-$k$ if and only if its columns are spanned by $k$ vectors obtained from evaluating $k$ complex sinusoids on the integer grid $0,\ldots,d-1$. I.e., every column of $T$ has a $k$-sparse Fourier representation (with arbitrary frequencies in the range $[0,1]$) and it follows that any sample $x^{(\ell)}\sim\mathcal{N}(0,T)$ does as well.

This observation allows for a resolution of the question above. In particular, Prony’s method [DP95] implies that any sum of $k$ complex sinusoids can be recovered from exactly $2k$ values. Accordingly, when $T$ is exactly rank-$k$, we can select $2k$ arbitrary entries in $x^{(\ell)}$, use Prony’s method to recover the remaining entries, and run our algorithm from Theorem 2.4 to approximate $T$. This approach achieves sample complexity linear in $k$ and logarithmic in $d$. In Section 5 we improve the logarithmic dependence to obtain a result that is fully independent of the dimension:

Theorem 2.9** (Dimension indep. complexity via Prony’s method).**

For any positive semidefinite Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$ with rank $k$, Problem 2.1 can be solved by Algorithm 2 with vector sample complexity $O\left(\frac{\log(k/\delta)}{\varepsilon^{2}}\right)$, entry sample complexity $\Theta(k)$, and total sample complexity $O\left(\frac{k\log(k/\delta)}{\varepsilon^{2}}\right)$.

Unfortunately, Prony’s method is notoriously sensitive to noise [SP95]: it fails to recover a sum of $k$ complex sinusoids from $2k$ values when those values are perturbed by even a very small amount of noise. The effect in our setting is that Theorem 2.9 relies heavily on $T$ being exactly rank-$k$.

To handle the more typical case when $T$ is only close to low-rank, we first prove that in this setting, $T$ is still closely approximated by a $k$-sparse Fourier representation. We then build on recent work on provably interpolating Fourier sparse functions with noise [CKPS16, CP18] to show how to recover this representation. Our proofs leverage several tools from randomized numerical linear algebra, including column subset selection bounds and the projection-cost preserving sketches of [CEM*+*15].

While this approach is somewhat involved (see Section 5 for the full development) a key takeaway is how it collects samples from each $x^{(\ell)}$. In contrast to ruler based methods, we generate a fixed pattern based on randomly sampling a subset of entries in $1,\ldots,d$ (see Figure 4). The randomized strategy concentrates samples near the edges of $1,\ldots,d$, producing a pattern that is intuitively similar to many sparse ruler designs, but much sparser. Overall, we obtain:

Theorem 2.10** (Sparse Fourier trans. sample complexity).**

For any positive semidefinite Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$ with rank $k$, Problem 2.1 can be solved with failure probability $\delta\leq\frac{1}{10}$ by Algorithm 3 with vector sample complexity $O\left(\frac{\log d\log(d/\varepsilon)}{\varepsilon^{2}}+k\right)$ and entry sample complexity $\Theta(k^{2}\log k\log(d/\varepsilon))$. The same bounds hold when $T$ is close to low-rank, with $\left(\lVert T-T_{k}\rVert_{2}+\frac{\mathrm{tr}(T-T_{k})}{k}\right)\cdot\mathrm{tr}(T)=O(\varepsilon^{2}\left\lVert T\right\rVert_{2}^{2})$.

See Section 5 for a full statement of Theorem 2.10. As a corollary, in Theorem 5.8, we also give a bound with sample complexity depending polynomially on the stable rank of $T$, $s=\frac{\mathrm{tr}(T)}{\lVert T\rVert_{2}}$.

We believe that improving Theorem 2.10 is an interesting and potentially impactful direction for future work. For example, it may be possible to improve the dependence on $k$ in the total sample complexity. We do give two incomparable lower bounds which argue that at least a linear dependence on $k$ is unavoidable. The first is against any non-adaptive algorithm:

Theorem 2.11** (Informal, see Theorem 6.1).**

Any algorithm that solves Problem 2.1 for any rank-$k$ Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$ by inspecting the entries of $n$ samples $x^{(1)},\ldots,x^{(n)}$ non-adaptively (i.e. the algorithm chooses, possibly randomly, a fixed set of positions to be read before seeing the samples), has total sample complexity $\Omega(k/\log k)$.

In particular, all of our upper bounds use non-adaptive sampling strategies, and so this theorem shows that they cannot be improved beyond potentially a $\tilde{O}(k^{2})$ factor. We further demonstrate that even adaptive algorithms cannot avoid this lower bound, if a small amount of noise is added:

Theorem 2.12** (Informal, see Theorem 6.6).**

Any algorithm that solves Problem 2.1 for any approximately rank-$k$ Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, by inspecting the entries of $n$ samples in a potentially adaptive manner, has total sample complexity $\Omega(k)$.

For the formal definition of what it means to be “approximately rank-$k$” in this context, see Theorem 6.6. The proof of this is inspired by lower bounds on sparse Fourier transforms and goes via information theory. We highlight one particularly interesting consequence of these bounds. By taking $k=\Omega(d)$ in the above lower bounds, one immediately obtains that no algorithm (even potentially adaptive), can avoid total sample complexity $\Omega(d)$ for Problem 2.1, without any rank constraints. Thus, in the metric of TSC, with no rank constraints, the simple algorithm of Theorem 2.4 is in nearly optimal for Problem 2.1. See Section 6 for a more detailed discussion.

2.4 Related work

The problem of estimating Toeplitz covariance matrices using all $d$ entries of each sample $x^{(\ell)}$ is a classical statistical problem. Methods for finding maximum likelihood or minimax estimators have been studied in a variety of settings [Bur67, KM81, BLW82] and we refer the reader to [BS97] for an overview. These methods seek to extract a Toeplitz matrix from the sample covariance to improve estimation performance. Existing approaches range from very simple (e.g., diagonal averaging as in our Theorem 2.4 [CRZ13]) to more complex (e.g., EM-like methods [MS87]). We are not aware of any work that establishes a generic non-asymptotic bound like Theorem 2.4.

Work on Toeplitz covariance estimation with the goal of minimizing entry sample complexity has focused on sparse ruler based sampling. This approach has been known for decades [Mof68, PBNH85] and is widely applied in signal processing applications: we refer the reader to [RATL16] for an excellent overview. There has been quite a bit of interest in finding rulers of minimal size – asymptotically, $\Theta(\sqrt{d})$ is optimal, but in many applications minimizing the leading constant is important [Lee56, Wic63, RSTL88]. There has been less progress on understanding rates of convergence for sparse ruler methods. [WZY17] provides some initial results on bounding $\|T-\tilde{T}\|_{2}$, but we are not aware of any strong relative error bounds like those of Theorem 2.6 or Theorem 2.8. [QP17] gives the closest analog, proving non-asymptotic vector sample complexity bounds for sparse ruler sampling. They measure error in terms of the $\ell_{2}$ norm of the first column of $T$, which can be translated to a Frobenius norm bound of the form $\lVert T-\tilde{T}\rVert_{F}\leq\epsilon\lVert T\rVert_{F}$, but not directly to a spectral norm bound.

More related to the algorithmic contributions of Section 5, are methods that make the additional assumption that $T$ is low-rank, or close to low-rank. The most well-known algorithms in this setting are the MUSIC [Sch81, Sch86, Pis73] and ESPRIT methods [RK89], although other approaches have been studied [BM86, FB96, ME09]. As discussed, Toeplitz $T$ is low-rank exactly when its columns are Fourier sparse functions. Accordingly, such methods are closely related to sparse recovery algorithms for “off-grid” frequencies. Most theoretical bounds for “off-grid” sparse recovery require strong separation assumptions that no frequencies are too close together [TBSR13, CFG14, TBR15, BCG*+*15]. A recent exception is [CKPS16], which we build on in this work since a lack of separation assumptions is essential to obtaining bounds for general low-rank Toeplitz matrices. Under separation assumptions on the frequencies in $T$’s Vandermonde decomposition, we note that [QP17] gives an algorithm with $O(\sqrt{k})$ ESC and $\operatorname{poly}(d,k,\epsilon)$ VSC when $T$ is close to rank-$k$, using a sparse recovery based approach.

Finally, we note that covariance estimation where $<d$ entries of $x^{(l)}$ are available has been studied without a Toeplitz assumption [GRESS16]. Without Toeplitz and thus Fourier structure, this work is more closely related to the matrix completion problem [CR09] than our setting.

2.5 Open questions

Our work initiates a comprehensive study of non-asymptotic complexity bounds for the recovery of Toeplitz covariance matrices, and leaves a number of open questions.

Efficiency. All of our estimators run in polynomial time except for Algorithm 3, which is used to obtain the bound of Theorem 2.10. An interesting open question is if one can use techniques from e.g., [CKPS16] to obtain polynomial time estimators for this setting.

Tighter bounds. While we present non-asymptotic sample complexity upper and lower bounds for many of the settings considered in this paper, and they are within polynomial factors of each other, it is an interesting open question to exactly characterize the sample complexities of these problems and to fully understand the tradeoffs between entry and vector sample complexity. Moreover, it is unclear if the form of the tail error we obtain for approximately rank-$k$ Toeplitz matrices is optimal, and it would be interesting to understand this better. A particularly interesting question is if the $O(\sqrt{k})$ entry sample complexity bound of [QP17] can be obtained for general low-rank $T$, without the additional assumptions made in that work.

Other settings / measures of error. It would be interesting to consider continuous analogs of this problem, for instance, on continuous or infinite length stationary signals. There has been substantial recent work on this subject [CTE18, AL12, LDTN11], and we believe that similar insights can be brought to bear here. There are also other interesting measures of error, for instance, relative Loewner ordering, or total variation distance the sampled and estimated distributions, which corresponds to recovering $T$ in the so-called Mahalanobis distance. Many of our results can be directly applied to these settings when $T$ is well-conditioned. However, it would be interesting to understand the sample complexities of these problems in more detail.

3 Notation and Preliminaries

Indexing. For a postive integer $z$ let $[z]\vcentcolon=\{1,\ldots,z\}$. For a set $S$, let $S^{z}$ denote the set of subsets of $S$ with $z$ elements. For a vector $v\in\operatorname{\mathbb{C}}^{d}$, let $v_{S}\in\operatorname{\mathbb{R}}^{|S|}$ denote the restriction of $v$ to the coordinates in $S$. For $T\in\operatorname{\mathbb{R}}^{d\times d}$ and set $R\subseteq[d]$ let $T_{R}\in\operatorname{\mathbb{R}}^{|R|\times|R|}$ denote the principle submatrix of $T$ with rows and columns indexed by $R$. Note that even when $T$ is Toeplitz, $T_{R}$ may not be.

Toeplitz matrix operations. For any vector $a=[{a}_{0},\ldots,{a}_{d-1}]\in\operatorname{\mathbb{R}}^{d}$, let $\mathrm{Toep}(a)\in\operatorname{\mathbb{R}}^{d\times d}$ denote the symmetric Toeplitz matrix whose entries are $\mathrm{Toep}(a)_{j,k}=a_{|j-k|}$. For any $M\in\operatorname{\mathbb{R}}^{d\times d}$, let $\mathrm{diag}(M)$ be the matrix given by setting all of $M$’s off diagonal entries to [math]. Let $\mathrm{avg}(M)$ be the symmetric Toeplitz matrix given by averaging the diagonals of $M$. Namely, $\mathrm{avg}(M)=\mathrm{Toep}(a)$ where for $s\in\{0,\ldots d-1\}$, $a_{s}=\frac{1}{\left|j,k\in[d]:|j-k|=s\right|}\sum_{j,k\in[d]:|j-k|=s}M_{j,k}$.

Linear algebra notation. For a matrix $A$, let $A^{T}$ and $A^{*}$ denote the transpose and Hermitian transpose, respectively. When all of $A$’s entries are real, $A^{*}=A^{T}$. For a vector $y\in\operatorname{\mathbb{C}}^{d}$ let $\|y\|_{2}=\sqrt{y^{*}y}$ denote the $\ell_{2}$ norm. For a matrix $A$ with $d$ columns, let $\|A\|_{2}=\sup_{x\in C^{d}}\|Ax\|_{2}/\|x\|_{2}$ denote spectral norm. Let $\|A\|_{F}$ denote the Frobenius norm.

A Hermitian matrix $A\in\operatorname{\mathbb{C}}^{d\times d}$ is positive semidefinite (PSD) if for all $x\in\operatorname{\mathbb{C}}^{d}$, $x^{*}Ax\geq 0$. Let $\preceq$ denote the Loewner ordering: $A\preceq B$ indicates that $B-A$ is PSD. Let $A^{+}$ denote the Moore-Penrose pseudoinverse of $A$. $A^{+}=V\Sigma^{-1}U^{*}$ where $U\Sigma V^{*}=A$ is the compact singular value decomposition of $A$. When $A$ is PSD, let $A^{1/2}=U\Sigma^{1/2}$ where $\Sigma^{1/2}$ is obtained by taking the entrywise square root of $\Sigma.$ Let $A_{k}$ denote the projection of $A$ onto its top $k$-singular vectors: $A_{k}=U_{k}\Sigma_{k}V_{k}^{*}$ where $U_{k},V_{k}\in\operatorname{\mathbb{C}}^{d\times k}$ contain the first $k$ columns of $U,V$ respectively and $\Sigma_{k}\in\operatorname{\mathbb{R}}^{k\times k}$ is diagonal matrix containing $A$’s largest $k$ singular values. Note that $A_{k}$ is an optimal low-rank approximation in that: $A_{k}=\operatorname*{arg\,min}_{\mathrm{rank}\,k\ M}\lVert A-M\rVert_{F}=\operatorname*{arg\,min}_{\mathrm{rank}\,k\ M}\lVert A-M\rVert_{2}$.

Gaussians, TV distance, KL divergence. We let $x\sim\mathcal{N}(m,T)$ indicates that $x$ is drawn from a Gaussian distribution with mean $m$ and covariance $T$. When $B\in\operatorname{\mathbb{C}}^{d\times d}$, it is an elementary fact that $x\sim\operatorname{\mathcal{N}}(0,BB^{*})$, can be written as $x=BUg$ where $g\sim\operatorname{\mathcal{N}}(0,I)$ and $U\in\operatorname{\mathbb{C}}^{d\times d}$ is a fixed unitary matrix. When $B$ is real, $U$ can be any orthogonal matrix in $\operatorname{\mathbb{R}}^{d\times d}$. For any two distributions $F,G$, let $d_{\mathrm{TV}}(F,G)=\frac{1}{2}\int|dF-dG|$ denote their total variation distance, and let $d_{\mathrm{KL}}(F,G)=\int\log\frac{dF}{dG}dF$ denote the KL divergence [CT12].

Probability notation. We let $\Pr[g]$ denote the probability that event $g$ occurs. Sometimes we add an additional subscript to remind the reader of any random events that effect $g$, and thus its probability. For example, we might write $\Pr_{x}[g]$ if $x$ is a random variable and $g$ depends on $x$. We use $:$ to denote “such that”. For example, $\Pr[\exists x\in\mathcal{X}:x>10]$ denotes the “probability that there exists some element $x$ in the set $\mathcal{X}$ such that $x$ is greater than $10$.”

Miscellaneous. Let $|z|=\sqrt{z^{*}z}$ denote the magnitude of a complex number $z$. Let $i$ denote $\sqrt{-1}$. When stating asymptotic bounds we use $\tilde{O}(f)$ as shorthand for $O(f\log^{c}f)$ for some constant $c$, and similarly we let $\widetilde{\Omega}(f)=\Omega(f/\log^{c}f)$ for some fixed $c$.

3.1 Fourier analytic tools

We repeatedly make use of the Fourier structure of Toeplitz matrices. To that end, we often work with asymmetric Fourier transform matrices whose rows correspond to “on-grid” integer points (e.g., samples in time domain) and whose columns correspond to “off-grid” frequencies.

Definition 3.1** (Fourier matrix).**

For a set $S=\{f_{1},\ldots f_{s}\}\subset\operatorname{\mathbb{C}}$, we define the Fourier matrix $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$ as the matrix with $j^{\text{th}}$ column given by $F(f_{j})$, where:

[TABLE]

While Toeplitz matrices cannot in general be diagonalized by the discrete Fourier transform like circulant matrices, the Carathéodory-Fejér-Pisarenko decomposition (also known as the Vandermonde decomposition) implies that they can be diagonalized by a Fourier matrix with off-grid frequencies [CF11, Pis73]. A general version of the classical result holding for all ranks $r\leq d$ is:

Fact 3.1** (Vandermonde Decomposition, Cor. 1 of [Cyb82]).**

Any PSD Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$ with rank $r\leq d$ can be written as $T=F_{S}DF_{S}^{*}$ where $D\in\operatorname{\mathbb{R}}^{r\times r}$ is positive and diagonal and $F_{S}\in\operatorname{\mathbb{C}}^{d\times r}$ is a Fourier matrix (as in Definition 3.1) with frequencies $S=\{f_{1},\ldots f_{r}\}\subset[0,1]$. Furthermore, $f_{1},\ldots f_{r}$ come in conjugate pairs with equal weights: i.e., for every $j\in[r]$, there is a $j^{\prime}$ (possibly with $j^{\prime}=j$) with $e^{-2\pi if_{j}}=e^{2\pi if_{j^{\prime}}}$ and $D_{j,j}=D_{j^{\prime},j^{\prime}}$.

One useful consequence of the Vandermonde decomposition is that, as $d\rightarrow\infty$, $T$ is diagonalized by the discrete time Fourier transform. We can take advantage of this property by noting that any Toeplitz matrix $\mathrm{Toep}(a)$ can be extended to be arbitrarily large by padding $a$ with [math]’s. The spectral norm of the extension upper bounds that of $\mathrm{Toep}(a)$. Accordingly, we can bound $\|\mathrm{Toep}(a)\|_{2}$ by the supremum of $a$’s DTFT. We refer the reader to [Mec07] for a formal argument, ultimately giving:

Fact 3.2** (See e.g., [Mec07]).**

For any $a=[a_{0},\ldots,a_{d-1}]\in\operatorname{\mathbb{R}}^{d}$,

[TABLE]

3.2 Linear algebra tools

Our methods with sample complexity logarithmic in $d$ for low-rank (or close to low-rank) $T$ do not select entries from each $x\sim\mathcal{N}(0,T)$ according to a fixed pattern like a sparse ruler. Instead, they use a randomly chosen pattern. We view this approach as randomly selecting a subset of rows from the $d\times n$ matrix of samples $[x^{(1)},\ldots,x^{(n)}]$. To find a subset that preserves properties of this matrix, we turn to tools from randomized numerical linear algebra and matrix sketching [Mah11, DM16]. Specifically, we adapt ideas from random selection algorithms based on importance sampling – i.e., when rows are sampled independently at random with non-uniform probabilities.

In particular, work on matrix sketching and graph sparsification has repeatedly used importance sampling probabilities proportional to the statistical leverage scores [SS11, DMM06], which allow for random matrix compressions that preserve significant information [CEM*+*15, MM17]. Formally, the leverage score is a measure defined for any row in a matrix $A\in\operatorname{\mathbb{C}}^{d\times s}$:

Definition 3.2**.**

The leverage score, $\tau_{j}(A)$, of the $j^{\text{th}}$ row $a_{j}\in\operatorname{\mathbb{C}}^{d\times 1}$ in $A\in\operatorname{\mathbb{C}}^{d\times s}$ is defined as:

[TABLE]

Note that that $\sum_{j=1}^{d}\tau_{j}(A)=\mathrm{tr}\left(A\left(A^{*}A\right)^{+}A^{*}\right)=\mathrm{rank}(A)\leq\min(d,s)$.

The leverage score can be seen as measuring the “uniqueness” of any row. This intuition is captured by the following lemmas, which are used extensively in prior work on leverage scores:

Fact 3.3** (Minimization Characterization).**

[TABLE]

Fact 3.4** (Maximization Characterization).**

Let $\left(Ay\right)_{j}$ denote the $j^{\text{th}}$ entry of the vector $Ay$.

[TABLE]

One consequence of Fact 3.3 is that $\tau_{j}(A)\leq 1$ for all $j$ since we can choose $y$ to be the $j^{\text{th}}$ standard basis vector. We defer further discussion and proofs of these bounds to Appendix A, where we also state a standard matrix concentration bound for sampling by leverage scores.

In our applications, we cannot compute leverage scores explicitly, since we use them to sample rows from $[x^{(1)},\ldots,x^{(n)}]$ without ever looking at this whole matrix. To cope with this challenge we take advantage of the Fourier structure of Toeplitz matrices (Claim 3.1), which ensures that each $x^{(\ell)}$ can be written as a Fourier sparse or nearly Fourier sparse function when $T$ is low-rank or close to low-rank. This allows us to extend recent work on nearly tight a priori leverage scores estimates for Fourier sparse matrices [CKPS16, CP18, CP19, AKM*+*19] to our setting. See Appendix A.

3.3 Distances between Gaussians

In our lower bound proofs, we repeatedly use the following bounds on the TV distance and KL divergence between Gaussians.

Fact 3.5**.**

Let $M_{1},M_{2}\in\operatorname{\mathbb{R}}^{d\times d}$ be PSD. Then,

[TABLE]

By Pinsker’s inequality [CT12], this implies that

[TABLE]

3.4 Concentration bounds

We make use of the following well known concentration inequality for quadratic forms:

Claim 3.6** (Hanson-Wright Inequality – see e.g., [RV13]).**

Let $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,I)$ be independent. There exists a universal constant $c>0$ so that for all $A\in\operatorname{\mathbb{R}}^{d\times d}$ and all $t>0$,

[TABLE]

4 Covariance Estimation by Ruler

We begin by establishing Theorems 2.4 and 2.6, which can be analyzed in a unified way. In particular, both bounds are achieved by Algorithm 1 (see below) with the ruler $R$ chosen to either equal $[d]$ for Theorem 2.4, or to be a suitable sparse ruler for Theorem 2.6 (see Definition 2.1). Algorithm 1 approximates $T=\mathrm{Toep}(a)$ by individually approximating each entry ${a}_{s}$ in $a$ by an average over all entries in each $x^{(j)}$ with distance $s$ from each other. Letting ${\tilde{a}}$ denote this approximation, the algorithm returns the Toeplitz matrix $\tilde{T}=\mathrm{Toep}({\tilde{a}})$.

When $R=[d]$, Algorithm 1 reduces to approximating $T$ by the empirical covariance with its diagonals averaged to make it symmetric and Toeplitz. In the other extreme, we can approximate $T$ using a sparse ruler with $|R|=O(\sqrt{d})$ (see Claim 2.5). It is easy to see the following:

Claim 4.1**.**

Algorithm 1 has vector sample complexity $n$, entry sample complexity $|R|$, and total sample complexity $n|R|$. It runs in time $O\left(n|R|^{2}\right)$ and outputs a symmetric Toeplitz matrix $\tilde{T}$.

We prove a generic bound on the accuracy of Algorithm 1 for any $n$ and ruler $R$ (the accuracy will depend on properties of $R$). We then instantiate this bound for specific choices of these parameters to give Theorem 2.4, Theorem 2.6, and also Theorem 2.8 when $T$ is approximately low-rank. We first define a quantity characterizing how extensively a ruler covers distances $0,\ldots,d-1$:

Definition 4.1** (Coverage coefficient, $\Delta(R)$).**

For any ruler $R\subseteq[d]$ let $\Delta(R)\vcentcolon=\sum_{s=0}^{d-1}\frac{1}{|R_{s}|}.$

The coverage coefficient of a ruler is smaller when distances are represented more frequently by the ruler, so intuitively, the sample complexity will scale with it. It is not hard to check that when $R=[d]$, $\Delta(R)=\Theta(\log d)$ and for any $R$ with $|R|=\Theta(\sqrt{d})$, $\Delta(R)=\Theta(d)$. It is also possible to construct rulers which interpolate between $\Theta(\sqrt{d})$ and $d$ sparsity:

Definition 4.2**.**

For $\alpha\in[1/2,1]$, let $R_{\alpha}$ be the ruler given by $R_{\alpha}=R^{(1)}_{\alpha}\cup R^{(2)}_{\alpha}$, where444We assume for simplicity of exposition that $d^{\alpha}$ and $d^{1-\alpha}$ are integers.

[TABLE]

For such rulers we prove the following general bound in Appendix B:

Lemma 4.2**.**

For any $\alpha\in[1/2,1]$, $R_{\alpha}$ has size $|R_{\alpha}|\leq 2d^{\alpha}$, and moreover

[TABLE]

Note that for $\alpha=1$ we recover a bound of $O(\log d)$ and for $\alpha=1/2$ a bound of $O(d)$. With the coverage coefficient defined, we prove our general bound for any ruler $R$:

Theorem 4.3** (Accuracy of Algorithm 1).**

Let $\tilde{T}$ be the output of Algorithm 1 run with ruler $R$. For any $\varepsilon\in(0,1]$, let $\kappa\vcentcolon=\min\left(1,\frac{\varepsilon^{2}\cdot\lVert T\rVert_{2}^{2}}{\Delta(R)\cdot\lVert T_{R}\rVert_{2}^{2}}\right)$. There exist universal constants $C,c>0$ such that

[TABLE]

In particular, if $n={\Theta}\left(\frac{\log\left({d/\varepsilon\delta}\right)}{\kappa}\right)$, then $\left\|T-\tilde{T}\right\|_{2}\leq\varepsilon\lVert T\rVert_{2}$ with probability at least $1-\delta$.

Proof.

As noted in Claim 4.1, $\tilde{T}$ is a symmetric Toeplitz matrix and so $T-\tilde{T}$ is a symmetric Toeplitz matrix. Let $e\in\operatorname{\mathbb{R}}^{d}$ be its associated vector, so that $e$ is a random variable. By Fact 3.2, to prove the theorem it suffices to show that there exist $C,c>0$ with:

[TABLE]

Let $N=\{0,\frac{\varepsilon}{Qd^{2}},\frac{2\varepsilon}{Qd^{2}},\ldots,1\}$, for some universal constant $Q$ sufficiently large ($Q=80\pi$ suffices). We will prove that there exists a universal constant $c>0$ so that the following two bounds hold:

[TABLE]

We first show how (2) and (3) together imply the desired bound. Condition on the event that $|{e}_{s}|\leq 10\|T\|_{2}$ for all $s$ and the event that $|L_{e}(x)|\leq\frac{\varepsilon}{2}\|T\|_{2}$ for all $x\in N$. By (2) and (3) and a union bound, we know that this happens with probability at least

[TABLE]

Observe that for all $x\in[0,1]$, we have

[TABLE]

where in the last line we use that we are conditioning on $|{e}_{s}|\leq 10\lVert T\rVert_{2}$ for all $s$. Now, for any $x\in[0,1]$, let $x^{\prime}\in N$ be so that $|x-x^{\prime}|\leq\varepsilon/(Qd^{2})$. Then, as long as $Q\geq 80\pi$, we have

[TABLE]

which proves the theorem, letting and $C=4Q$.

It thus suffices to prove (2) and (3). We first prove (2). We can write for any $s\in 0,\ldots,d-1$:

[TABLE]

Observe that for each $\ell=1,\ldots,n$, we have that $x^{(\ell)}_{j}\cdot x^{(\ell)}_{k}$ is a subexponential random variable with mean $T_{j,k}$ and second moment which can be written using Isserlis’s theorem as

[TABLE]

The desired bound then follows from standard bounds on subexponential random variables [Wai19].

We now turn our attention to (3). Fix $x\in[0,1]$, and associate to it the Toeplitz matrix $M\in\operatorname{\mathbb{R}}^{d\times d}$, where for all $s\in 0,\ldots,d-1$ and $(j,k)\in[d]\times[d]$ with $|j-k|=s$, we let

[TABLE]

Letting $\overline{T}$ be the empirical covariance matrix $\overline{T}\vcentcolon=\frac{1}{n}\sum_{\ell=1}^{n}x^{(\ell)}{x^{(\ell)}}^{T}$, we can see that

[TABLE]

The last equality follows from the fact that $M$ is symmetric Toeplitz and $\tilde{T}_{R}$ is just obtained by averaging the diagonal entries of $\overline{T}_{R}$. That is, for any $s\in 0,\ldots d-1$,

[TABLE]

and thus

[TABLE]

We will show that for any $\varepsilon\in[0,1]$, there exists $c>0$ such that:

[TABLE]

Notice that if we let $x^{(\ell)}_{R}\in\operatorname{\mathbb{R}}^{|R|}$ denote the $\ell^{th}$ sample restricted to the indices in $R$ we have $x^{(\ell)}_{R}\sim\operatorname{\mathcal{N}}(0,T_{R})$. Therefore, if we let $y^{(\ell)}=T_{R}^{-1/2}x^{(\ell)}_{R}$, the $y^{(\ell)}$ are i.i.d., with $y^{(\ell)}\sim\operatorname{\mathcal{N}}(0,I_{R\times R})$. Moreover, we can rewrite the quantity in (5) as:

[TABLE]

where $M^{\prime}=T_{R}^{1/2}MT_{R}^{1/2}$. Observe that $\|M\|_{F}^{2}=\sum_{s=0}^{d-1}|R_{s}|\cdot\frac{\cos(2\pi sx)^{2}}{|R_{s}|^{2}}\leq\Delta(R)$. Additionally, $\|M^{\prime}\|_{F}\leq\|T_{R}\|_{2}\|M\|_{F}$, and $\|M^{\prime}\|_{2}\leq\|T_{R}\|_{2}\|M\|_{2}$. Thus, applying the Hanson-Wright inequality (Claim 3.6):

[TABLE]

recalling that $\kappa\vcentcolon=\min\left(1,\frac{\varepsilon^{2}\cdot\lVert T\rVert_{2}^{2}}{\Delta(R)\cdot\lVert T_{R}\rVert_{2}^{2}}\right)$. This gives (5). Combining with (4), for any fixed $x\in[0,1]$:

[TABLE]

and so (3) follows by a union bound over all $x\in N$. This completes the proof of Theorem 4.3. ∎

4.1 Applications of the general bound, full rank matrices

We describe a few specific instantiations of Theorem 4.3. When $R=[d]$ is the full ruler, $\Delta(R)=O(\log d)$ by Lemma 4.2 and $\|T_{R}\|_{2}=\|T\|_{2}$. This immediately yields:

Theorem 2.4 (Near linear sample complexity – full theorem).

Let $\tilde{T}$ be the output of Algorithm 1 run with the full ruler $R=[d]$. There exist universal constants $C,c>0$ such that, for any $\varepsilon\in(0,1]$:

[TABLE]

In particular, if $n={\Omega}\left(\frac{\log\left({d}/{\varepsilon\delta}\right)\log d}{\varepsilon^{2}}\right)$, then $\left\|T-\tilde{T}\right\|_{2}\leq\varepsilon\lVert T\rVert_{2}$ with probability at least $1-\delta$.

We can also apply Theorem 4.3 to any ruler with $|R|=\Theta(\sqrt{d})$ (e.g., the ruler in Claim 2.5, or more optimal constructions). In this case, $\Delta(R)=O(d)$ and can bound $\|T_{R}\|_{2}\leq\|T\|_{2}$ to give:

Theorem 2.6 (Sparse ruler sample complexity – full theorem).

Let $\tilde{T}$ be the output of Algorithm 1 run with any $\Theta(\sqrt{d})$-sparse ruler $R$. There exist universal constants $C,c>0$ such that, for any $\varepsilon\in(0,1]$:

[TABLE]

*In particular, for $n={O}\left(\frac{d\log\left({d/\varepsilon\delta}\right)}{\varepsilon^{2}}\right)$, $\left\|T-\tilde{T}\right\|_{2}\leq\varepsilon\lVert T\rVert_{2}$ with probability at least $1-\delta$. *

Finally, we record a bound for rulers $R_{\alpha}$ with any $\alpha\in[1/2,1]$. This result smoothly interpolates between Theorems 2.4 and 2.6, giving a natural way of trading between entry sample complexity and vector sample complexity, which might be valuable in practice. For example, by setting $\alpha=3/4$, we obtain an algorithm with ESC $O(d^{3/4})$ and VSC $\tilde{O}(d^{1/2})$.

Theorem 4.4**.**

Let $\alpha\in[1/2,1]$ be fixed. Let $\tilde{T}$ the output of Algorithm 1 run with the ruler $R_{\alpha}$ of Def. 4.2. There exist universal constants $C,c>0$ such that, for any $\varepsilon\in(0,1]$:

[TABLE]

In particular, for any $\delta>0$, if

[TABLE]

then $\left\|T-\tilde{T}\right\|_{2}\leq\varepsilon\lVert T\rVert_{2}$ with probability at least $1-\delta$.

As discussed in Section 2, while obtaining optimal entry sample complexity for any method that reads a fixed subset of entries of each sample, Theorem 2.6 gives total sample complexity that is worse than Theorem 2.4 by roughly a $\sqrt{d}$ factor. In Appendix B.1 we prove a lower bound demonstrating that this is tradeoff is inherent:

Theorem 4.5**.**

Let $R\subseteq[d]$ be a sparse ruler with $|R|\leq d^{\alpha}$, for some $\alpha\in[1/2,1]$, and let $\varepsilon>0$ be sufficiently small. Let $\mathscr{A}$ be a (possibly randomized) algorithm that takes $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,T)$ for some unknown, positive semidefinite and Toeplitz $T$, and for all $\ell$, looks only at the coordinates of $x^{(\ell)}$ in $R$ before outputting $\tilde{T}$. If, for any $T$, $\|\tilde{T}-T\|_{2}\leq\varepsilon\lVert T\rVert_{2}$ with probability $\geq 1/10$, then we must have $n=\Omega(d^{3-4\alpha}/\varepsilon^{2})$. The algorithm thus requires total sample complexity $\Omega(d^{3-3\alpha}/\varepsilon^{2})$.

Note that when $\alpha=1/2$, Theorem 4.5 implies that we require $\Omega(d/\varepsilon^{2})$ vector samples, which nearly matches the upper bound of Theorem 2.6 and confirms the simulation results in Figure 3. Interestingly in Theorem 4.5 we get some tradeoff for all $\alpha\in[1/2,1]$, but it does not match the upper bound of Theorem 4.4. We leave it as an open question to resolve this gap.

4.2 Applications of the general bound, low-rank matrices

As we saw in Figure 3, sparse ruler estimation can actually outperform full ESC methods in terms of total sample complexity when $T$ is low-rank. We support this observation theoretically, giving a much tighter bound than Theorem 2.6 when $T$ is close to low-rank and $R$ is the sparse ruler of Claim 2.5 or one of the family of rulers of Definition 4.2.

To do so, we critically use that a low-rank matrix cannot concentrate significant mass on more than a few small principal submatrices. A version of this fact, showing that a low-rank matrix cannot concentrate on many diagonal entries is shown in [MMW19]. We use as similar argument.

Lemma 4.6**.**

Consider any partition $R_{1}\cup\ldots\cup R_{t}=[d]$. For any $T\in\operatorname{\mathbb{R}}^{d\times d}$, let $T_{R_{j}}$ be the principal sub-matrix corresponding to $R_{j}$. If $T$ is rank-$k$, then for some $S\subseteq[t]$ with $|S|\leq\frac{k}{\varepsilon}$, for all $\ell\in[t]\setminus S$:

[TABLE]

where $T_{(R_{\ell},[d])}$ is the submatrix of $T$ given by selecting the rows in $R_{\ell}$ and all columns.

Proof.

Letting $\tau_{j}(T)$ be $T$’s $j^{\text{th}}$ leverage score (Definition 3.2) we have $\sum_{j=1}^{d}\tau_{j}(T)=\mathrm{rank}(T)=k$. Thus, excluding a set $S\subseteq[t]$ with $|S|\leq\frac{k}{\varepsilon}$ indices, for all $\ell\in[t]\setminus S$, we have $\sum_{j\in R_{\ell}}\tau_{j}(T)\leq\varepsilon$. Applying the maximization characterization of Fact 3.4, for any $k\in R_{\ell}$, letting $e_{k}$ be the $k^{th}$ standard basis vector and $T_{(k,[d])}$ denote the $k^{\text{th}}$ column of $T$:

[TABLE]

Applying this bound to all $k\in R_{\ell}$ we have $\frac{\lVert T_{R_{\ell}}\rVert_{F}^{2}}{\lVert T_{(R_{\ell},[d])}\rVert_{F}^{2}}\leq\varepsilon$, which completes the lemma. ∎

In Appendix B we use Lemma 4.6 to show:

Lemma 4.7**.**

Let $\alpha\in[1/2,1]$. For any $k\leq d$, any PSD Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, and the sparse ruler $R_{\alpha}$ of Def. 4.2,555For $\alpha=1/2$ the bound also applies to the construction of Claim 2.5, which is essentially identical to that of Def. 4.2, but applies when $\sqrt{d}$ is not an integer.

[TABLE]

where $\displaystyle T_{k}=\operatorname*{arg\,min}_{\mathrm{rank}-k\ M}\lVert T-M\rVert_{F}=\operatorname*{arg\,min}_{\mathrm{rank}-k\ M}\lVert T-M\rVert_{2}$. If $T$ is rank-$k$, $\lVert T-T_{k}\rVert_{F}^{2}=\lVert T-T_{k}\rVert_{2}^{2}=0$.

Plugging Lemma 4.7 in Theorem 4.3 we obtain:

Theorem 2.8 (Sublinear sparse ruler sample complexity – full theorem).

Let $\alpha\in[1/2,1]$. Let $\tilde{T}$ be the output of Algorithm 1 run with the sparse ruler $R_{\alpha}$ of Def. 4.2. There exist universal constants $C,c>0$ such that, for any $k\leq d$ and $\varepsilon\in(0,1]$:

[TABLE]

This bound implies, for example, that if if we take $\alpha=1/2$, and $\frac{\sqrt{d}}{k^{2}}\cdot\lVert T-T_{k}\rVert_{F}^{2}\leq\lVert T\rVert_{2}^{2}$ or $\frac{d}{k^{2}}\cdot\lVert T-T_{k}\rVert_{2}^{2}\leq\lVert T\rVert_{2}^{2}$, and $n={O}\left(\frac{k^{2}\cdot\log\left({d}/{\varepsilon\delta}\right)}{\varepsilon^{2}}\right)$, then $\left\|T-\tilde{T}\right\|_{2}\leq\varepsilon\lVert T\rVert_{2}$ with probability at least $1-\delta$. Note that, unlike Theorem 2.6, which applies to all sparse rulers, we only prove Theorem 2.6 for the specific family of sparse rulers $R_{\alpha}$ (it also applies to the essentially identical ruler of Claim 2.5 for $\alpha=1/2$). We conjecture that it should hold more generally.

Proof.

By Lemma 4.7 we have

[TABLE]

By Lemma 4.2, we also have $\Delta(R_{\alpha})\leq 2d^{2-2\alpha}+O(d^{1-\alpha}\cdot\log d)$. Applying Theorem 4.3 with

[TABLE]

and plugging in our bound on $\Delta(R_{\alpha})$ gives the theorem. ∎

5 Covariance Estimation by Sparse Fourier Transform

In this section we describe an estimation algorithm that further improves on sparse ruler methods for low-rank or nearly low-rank $T$. When $k$ is fixed, we almost eliminate dependence on $d$ in our total sample complexity, reducing the $\sqrt{d}$ in Theorem 2.8 to a logarithmic dependence.

Our approach heavily uses the Fourier structure of Toeplitz matrices, in particular that a rank-$k$ Toeplitz matrix can be written in the span of $k$ off-grid frequencies. Using sparse Fourier transform techniques, we can recover (approximations to) these frequencies using very few samples.

5.1 Recovering exactly low-rank matrices via Prony’s method

We first consider the case when $T$ is exactly rank-$k$ for some $k\leq d$. In this case, we give a short argument that Problem 2.1 can be solved with $O(k\log k/\varepsilon^{2})$ total sample complexity, which is completely independent of the ambient dimension $d$. Our main contribution will be to extend this claim to the practical setting where $T$ is not precisely low-rank. However, we discuss the exactly low-rank case first to exposit the main idea.

Consider $x^{(1)},\ldots,x^{(n)}$ drawn independently from $\mathcal{N}(0,T)$ and let $X\in\operatorname{\mathbb{R}}^{d\times n}$ be the matrix with the samples as its columns. Via the Vandermonde decomposition of Lemma 3.1 we can write $T$ as $T=F_{S}DF_{S}^{*}$ where $D\in\operatorname{\mathbb{R}}^{k\times k}$ is positive and diagonal and $F_{S}\in\operatorname{\mathbb{C}}^{d\times k}$ is a Fourier matrix with columns corresponding to some set of $k$ frequencies $S=\{f_{1},\ldots f_{k}\}$. Note that $X$ is distributed as $T^{1/2}G$ where each entry of $G$ is drawn i.i.d. from $\mathcal{N}(0,1)$. Equivalently, there exists some unitary matrix $U\in\operatorname{\mathbb{C}}^{d\times d}$ so that $X$ is distributed as $F_{S}D^{1/2}UG=F_{S}Z$ where we let $Z\in\operatorname{\mathbb{C}}^{k\times n}$ denote $D^{1/2}UG$. Note that each column of $F_{S}Z$ is Fourier $k$-sparse. This structure lets us recover any column of $Z$ from $2k$ entrywise measurements of the corresponding column in $X$, using Prony’s method [DP95]. Moreover, the method lets us simultaneously recover the frequencies $S$, giving us a representation of $X=F_{S}Z$.

For any vector $x\in\operatorname{\mathbb{R}}^{d}$, let $P_{k}(x)\in\operatorname{\mathbb{R}}^{k\times k}$ denote the matrix whose entries are given by $(P_{k}(x))_{i,j}=x_{i+j-1}$, and let $b_{k}(x)\in\operatorname{\mathbb{R}}^{k}$ be given by $(b_{k}(x))_{i}=x_{k+i}$. Observe that $P_{k}(x)$ and $b_{k}(x)$ can be simultaneously formed from $2k$ entrywise observations from $x$. We require the following classical lemma, underlying Prony’s method, which states that a signal with a $k$-sparse Fourier transform can be recovered from $2k$ measurements.

Lemma 5.1** ([DP95]).**

Let $S$ be a set of frequencies with $|S|=k$, and let $x=F_{S}y$ for some unknown $y\in\operatorname{\mathbb{R}}^{k}$. Let $c\in\operatorname{\mathbb{R}}^{k}$ be the solution to the linear equation $P_{k}(x)c=-b_{k}(x)$, and define the polynomial $p(t)=\sum_{s=1}^{k}c_{s}t^{s}$. Let $R=R(x)=\{r_{1},\ldots,r_{k^{\prime}}\}$ be the roots of $p$. Then (1) we have $R\subseteq S$, and (2) we have that $x$ is in the column space of $F_{R}$.

For simplicity, in this section we show how Lemma 5.1, combined with an exact polynomial root finding oracle, yields an algorithm that recovers $T$ using only $O(k\log k/\varepsilon^{2})$ samples and time. In reality, such an oracle does not exist, so in Appendix C.2 we show how to use approximate root finding algorithms to actually instantiate this algorithm. We give two algorithms: the first pays logarithmically in $d$ in terms of sample complexity, and polynomially in $d$ in terms of runtime. The second recovers the $O(k\log k/\varepsilon^{2})$ sample complexity and pays just a $\log\log d$ factor and a $\log\log\kappa$ factor in runtime, where $\kappa$ is the condition number of $T$ (the ratio between its largest and smallest non-zero eigenvalues). It is stronger except in the extreme case when $\kappa=\Omega(2^{2^{d}})$.

Theorem 2.9 (Dimension independent sample complexity – full theorem).

Let $\widehat{T}$ be the output of Algorithm 2 run on samples $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,T)$ for a rank-$k$ Toeplitz covariance $T=F_{S}DF_{S}^{*}$. Suppose there is an oracle $\mathscr{O}$ which takes a degree $k$ polynomial $p$ with roots $S=\{z_{1},\ldots,z_{k}\}$ in the unit complex disc and exactly returns the set $S$. Then there is a universal constant $C>0$ so that for any $\varepsilon\in(0,1)$:

[TABLE]

Moreover, the algorithm requires $O(kn)$ total samples, $O(n)$ queries to $\mathscr{O}$, and $\operatorname{poly}(k,n)$ additional runtime. In particular, for any $\delta>0$, if we let $n=\Theta(\log(k/\delta)/\varepsilon^{2})$, then the algorithm outputs $\hat{T}$ so that $\Pr\left[\lVert T-\widehat{T}\rVert_{2}>\varepsilon\lVert T\rVert_{2}\right]<\delta$, using $\Theta(k\log(k/\delta)/\varepsilon^{2})$ entrywise samples and $\operatorname{poly}(k,1/\varepsilon,\log 1/\delta)$ runtime.

Proof.

The sample and runtime bounds are immediate by inspection of Algorithm 2. Note that the algorithm does not explicitly output $\widehat{T}$, as that would require $d^{2}$ time. Instead, it outputs a set of $k$ frequencies $R$ and a $k\times k$ diagonal matrix $\hat{D}$. As in the pseudocode, we let $\widehat{T}=F_{R}\hat{D}F_{R}^{*}$. Thus the remainder of this section is dedicated to a proof of correctness.

As discussed, for all $j\in[n]$, we can write $x^{(j)}=F_{S}D^{1/2}Ug^{(j)}$, where $g^{(j)}\sim\operatorname{\mathcal{N}}(0,I)$ and $U$ is unitary. Letting $y^{(j)}=D^{1/2}Ug^{(j)}$, we see that with probability $1$, $y^{(j)}$ has full support. Since $F_{S}$ has full column rank (or $T$ would have rank $<k$), we thus see that with probability $y^{(j)}$ will not fall in the span of $F_{R}$ for any $R\subset S$. Therefore, by Lemma 5.1, the oracle $\mathscr{O}$ outputs the full set $\{e^{2\pi if_{1}},\ldots,e^{2\pi if_{k}}\}$, where $S=\{f_{1},\ldots,f_{k}\}$. In particular $R_{1}=R_{j}=S$ for all $j\in[n]$ with probability $1$. As a result, with probability $1$, we have that $\hat{y}^{(j)}=y^{(j)}$ for all $j=1,\ldots,n$. The random variable $Z_{j\ell}\vcentcolon=|y^{(j)}_{\ell}|^{2}=D_{\ell\ell}\left|\left(Ug^{(j)}\right)_{\ell}\right|^{2}$ is a sub-exponential random variable with expecation $D_{\ell\ell}$. Hence, if we let $Z_{\ell}=\frac{1}{n}\sum_{j\in[n]}Z_{j\ell}$, we have

[TABLE]

for some universal constant $C$. In Line 4 of Algorithm 2 we set $\hat{D}=\mathrm{diag}(Z_{1},\ldots,Z_{k})$. By a union bound we thus have

[TABLE]

In other words, if the event of (7) holds, for any vector $w$, $|w^{T}(\hat{D}-D)w|\leq\varepsilon\cdot w^{T}Dw.$ This gives that for any $w$, $|w^{T}F_{S}(\hat{D}-D)F_{S}^{*}w|\leq\varepsilon\cdot w^{T}F_{S}DF_{S}^{*}w$ and thus $\left\lVert F_{S}\hat{D}F_{S}^{*}-T\right\rVert_{2}\leq\varepsilon\left\lVert T\right\rVert_{2}$. Thus by (7) we have as claimed:

[TABLE]

∎

5.2 Approximately low-rank matrices

Unlike in the exact low-rank case, when $T$ is just close to low-rank, the samples $x^{(1)},\ldots,x^{(n)}$ do not exactly have a $k$-sparse Fourier representation, and we cannot recover them exactly using Prony’s method. However, we can prove that these samples are still approximately Fourier $k$-sparse. Thus, we will be able to recover an approximation to the samples and the empirical covariance matrix. We start with a short lemma, which guarantees that a good approximation to the empirical covariance matrix $XX^{T}$ suffices to obtain a good approximation to $T$.

Lemma 5.2**.**

Consider PSD Toeplitz $T\in\operatorname{\mathbb{R}}^{d\times d}$ and let $X\in\operatorname{\mathbb{R}}^{d\times n}$ have columns drawn i.i.d. from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. If $n={\Omega}\left(\frac{\log\left({d}/{\varepsilon\delta}\right)\log d}{\varepsilon^{2}}\right)$ then with probability at least $1-\delta$, for any $B\in\operatorname{\mathbb{R}}^{d\times d}$:

[TABLE]

where $\mathrm{avg}(B)$ is the symmetric Toeplitz matrix given by averaging the diagonals of $B$ (see Sec. 3).

Proof.

We have by triangle inequality:

[TABLE]

where the last inequality follows from the fact that $\mathrm{avg}(\cdot)$ can only decrease Frobenius norm. The lemma follows by applying Theorem 2.4, which shows that, for $n={\Omega}\left(\frac{\log\left({d}/{\varepsilon\delta}\right)\log d}{\varepsilon^{2}}\right)$, with probability at least $1-\delta$, $\lVert T-\mathrm{avg}(XX^{T})\rVert_{2}\leq\varepsilon\lVert T\rVert_{2}$. ∎

5.2.1 Existence of frequency based low-rank approximation

With Lemma 5.2 in place, our goal is to show how to obtain $B$ with small $\lVert B-XX^{T}\rVert_{F}$ using a sublinear number of reads from each sample of $\mathcal{N}(0,T)$ (i.e., by loooking at a sublinear number of rows of $X$). For general $X$, this would be impossible. However, since the columns of $X$ are distributed as $\mathcal{N}(0,T)$, if $T$ is close to low-rank we can argue that there is a $B$ that well approximates $XX^{T}$ and further is spanned by a small number of frequencies. Using sparse Fourier transform techniques similar at a high level to the MUSIC algorithm [Sch81], we can sample-efficiently recover such a $B$, without reading all of $X$.

We again start with the Vandermonde decomposition of Lemma 3.1, writing $T$ as $T=F_{S}DF_{S}^{*}$ where $D\in\operatorname{\mathbb{R}}^{d\times d}$ is nonnegative and diagonal and $F_{S}$ is a Fourier matrix with columns corresponding to some set of frequencies $S=\{f_{1},\ldots f_{d}\}$. If $T$ were circulant, $f_{1},\ldots,f_{d}$ would be ‘on-grid’ orthogonal frequencies and the columns of $F_{S}$ would be eigenvectors of $T$. Thus, the columns corresponding to the largest $k$ entries of $D$ (the top $k$ eigenvalues) would span an optimal rank-$k$ approximation of $T$. While this is not the case for general Toeplitz matrices, we can still apply the well known technique of column subset selection to prove the existence of a set of $O(k)$ frequencies spanning a near optimal low-rank approximation to $T$. Specifically in Appendix C we show:

Lemma 5.3** (Frequency-based low-rank approximation).**

For any PSD Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, rank $k$, and $m\geq ck$ for some fixed constant $c$, there exists $M=\{f_{1},\ldots,f_{m}\}\subset[0,1]$ such that, letting $F_{M}\in\operatorname{\mathbb{C}}^{d\times m}$ be the Fourier matrix with frequencies $M$ (Def. 3.1) and $Z=F_{M}^{+}T^{1/2}$, we have 1): $\lVert F_{M}^{+}\rVert_{2}^{2}\leq\frac{2}{\beta}$ and 2):

[TABLE]

Note that a bound on $\lVert F_{M}^{+}\rVert_{2}$ will eventually be important in arguing that we can discretize $[0,1]$ to search for $M$ with frequencies on some finite grid, without incurring too much error. If we did not require this bound, we could obtain (8) and (9) without the additive error depending on $\beta$.

Lemma 5.3 shows the existence of $m=O(k)$ frequencies that span an $m$-rank approximation to $T^{1/2}$ that is nearly as good as the best rank-$k$ approximation. However, to apply Lemma 5.2 we must identify a subset of frequencies $M$ such that $F_{M}Z$ approximates not $T^{1/2}$ itself but $X$ with columns drawn from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. $X$ is distributed as $T^{1/2}G$ where each entry of $G$ is distributed as $\mathcal{N}(0,1/\sqrt{n})$. Thus $X$ can be viewed as a random ‘sketch’ of $T^{1/2}$, approximating many of its properties [DM16, Woo14]. In particular, we employ a projection-cost-preserving sketch property, following from Theorems 12 and 27 of [CEM*+*15]666[CEM*+*15] considers only real valued matrices, however their proofs can be extended to complex matrices. Alternatively, the bound can be shown by rotating $F_{S}$ to a real matrix. See Appendix C., which in particular implies that if $T^{1/2}$ can be well approximated by projecting onto a set of $m$ columns of $F$ (i.e., $m$ frequencies), then as long as $n\approx O(m)$, so can $X$:

Lemma 5.4** (Projection-cost-preserving sketch).**

Consider PSD $T\in\operatorname{\mathbb{R}}^{d\times d}$ and $X\in\operatorname{\mathbb{R}}^{d\times n}$ with columns drawn i.i.d. from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. For any rank $m$ and $\gamma,\delta\in(0,1]$, if $n\geq\frac{c(m+\log(1/\delta))}{\gamma^{2}}$ for sufficiently large $c$, then with probability $\geq 1-\delta$, for all $F_{M}\in\operatorname{\mathbb{C}}^{d\times m}$:

[TABLE]

Applying Lemma 5.4 along with the frequency subset selection result of Lemma 5.3 we conclude (see Appendix C for full proof):

Lemma 5.5**.**

Consider PSD Toeplitz $T\in\operatorname{\mathbb{R}}^{d\times d}$ and $X\in\operatorname{\mathbb{R}}^{d\times n}$ with columns drawn i.i.d. from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. For any rank $k$, $\varepsilon,\delta\in(0,1]$, $m\geq c_{1}k$, and $n\geq c_{2}\left(m+\log(1/\delta)\right)$ for sufficiently large $c_{1},c_{2}$, with probability $\geq 1-\delta$, there exists $M=\{f_{1},\ldots,f_{m}\}\subset[0,1]$ such that, letting $Z=F_{M}^{+}X$, $\lVert Z\rVert_{2}^{2}\leq\frac{c_{3}\cdot d^{2}\lVert T\rVert_{2}}{\varepsilon^{2}}$ for some fixed $c_{3}$ and:

[TABLE]

5.2.2 Computing a frequency based low-rank approximation

With Lemma 5.5 in place, it remains to show that, with few samples from $X$, we can actually find an $M\subset[0,1]$ and $W=ZZ^{*}$ (approximately) satisfying the bound of the lemma. To do this, we will consider a large but finite number of possible sets obtained by discretizing $[0,1]$ and will brute force search over these sets. The critical step required to perform this search is to compute a nearly optimal $W$ for each $M$ while only looking at a subset of rows in $X$. We can argue that it is possible to do so if we could sample $X$ using the row leverage scores of $F_{M}$ (Def. 3.2), which are well known to be useful in approximately solving regression problems, such as the one required to compute $W$ [DMM06, Woo14]. Since leverage score sampling succeeds with high probability, by a union bound this would allow us to find a near optimal $W$ for every subset $M$ in our brute-force search, and thus identify a near optimal subset.

Unfortunately, the leverage scores of $F_{M}$ depend on $M$ and will be different for different frequency sets. Naively, we will have to sample from a different distribution for each $M$. Over the course of the search, this could require reading all $d$ rows of the $X$ (i.e., all entries in our samples). To deal with this challenge, we extend recent work on a priori leverage score bounds for Fourier sparse functions. As we prove in Appendix A, there is a fixed function that well approximates the leverage scores of any $F_{M}$ with $|M|=m$. See Figure 4 for an illustration of the sampling distribution corresponding to this function. This fact lets us sample a single set of $\tilde{O}(m)$ rows using this distribution and, by a union bound, find a near optimal $W$ for all $M$ in our search.

We start by claiming that discretizing frequencies to a net does not introduce too much additional error. Formally, we prove the following in Appendix C:

Lemma 5.6**.**

Consider PSD Toeplitz $T\in\operatorname{\mathbb{R}}^{d\times d}$ and $X\in\operatorname{\mathbb{R}}^{d\times n}$ with columns drawn i.i.d. from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. For any rank $k$ and $\delta,\varepsilon\in(0,1]$, consider $m\geq c_{1}k$, $n\geq c_{2}\left(m+\log(1/\delta)\right)$, and $N=\{0,\alpha,2\alpha,\ldots 1\}$ for $\alpha=\frac{\varepsilon^{2}}{c_{3}d^{3.5}}$ for sufficiently large constants $c_{1},c_{2},c_{3}$. With probability $\geq 1-\delta$, there exists $M=\{f_{1},\ldots,f_{m}\}\subset N$ with:

[TABLE]

We next show how to sample efficiently find $M\subset N$ and $W\in\operatorname{\mathbb{C}}^{d\times d}$ that nearly minimize $\lVert F_{M}WF_{M}^{*}-XX^{T}\rVert_{F}$ and thus (approximately) satisfy the guarantee of Lemma 5.6.

Lemma 5.7**.**

Consider any $N\subset[0,1]$, rank $m$, and $B\in\operatorname{\mathbb{R}}^{d\times d}$. Let $S_{1}\in\operatorname{\mathbb{R}}^{s_{1}\times d}$ and $S_{2}\in\operatorname{\mathbb{R}}^{s_{2}\times d}$ be independent sampling matrices sampled using the scheme of Claim A.1 with parameters $\varepsilon=\frac{1}{c_{1}}$ and $\delta=\frac{1}{c_{2}N^{m}}$ and the leverage score distribution of Cor. C.2 with frequency set size $2m$ for sufficiently large constants $c_{1},c_{2}$. With probability at least $97/100$, $\max(s_{1},s_{2})\leq cm^{2}\log m\log(|N|)$ for sufficiently large $c$ and letting $\displaystyle\tilde{M},\tilde{W}=\operatorname*{arg\,min}_{W\in\operatorname{\mathbb{C}}^{m\times m},M\in N^{m}}\lVert S_{1}F_{M}WF_{M}^{*}S_{2}^{T}-S_{1}BS_{2}^{T}\rVert_{F}$,

[TABLE]

Proof.

Let $\displaystyle\hat{M},\hat{W}=\operatorname*{arg\,min}_{W\in\operatorname{\mathbb{C}}^{m\times m},M\in N^{m}}\lVert F_{M}\tilde{W}F_{M}^{*}-B\rVert_{F}.$ We will show that with probability $\geq 97/100$, for all $M\in N^{m}$ and $W\in\operatorname{\mathbb{C}}^{m\times m}$,

[TABLE]

Conditioning on (10) holding, the lemma follows by applying it twice to bound:

[TABLE]

We now prove that (10) holds with probability $\geq 97/100$. By triangle inequality, for any $M,W$:

[TABLE]

Note that by design of the sampling scheme in Claim A.1, both $S_{1}$ and $S_{2}$ preserve the squared Frobenius norm in expectation. I.e., for any matrix $C\in\operatorname{\mathbb{C}}^{d\times p}$, $\operatorname*{\mathbb{E}}[\lVert S_{1}C\rVert_{F}^{2}]=\operatorname*{\mathbb{E}}[\lVert S_{2}C\rVert_{F}^{2}]=\lVert C\rVert_{F}^{2}$. Thus, by Markov’s inequality, with probability $\geq 99/100$,

[TABLE]

In turn, fixing $S_{1}$, with probability $99/100$,

[TABLE]

Overall, with probability $\geq 98/100$,

[TABLE]

Additionally, by Corollary C.2, for any $M\in N^{m}$, $S_{1}$ is sampled using leverage scores upper bounds of $[F_{M},F_{\hat{M}}]$ which sum to $O(m\log m)$. Thus if we apply Claim A.1 with error parameter $\varepsilon=\frac{1}{32}$ and some failure probability $\delta$, with probability $\geq 1-\delta$, $s_{1}\leq cm\log m\cdot\log\left(\frac{m}{\delta}\right)$ for sufficiently large $c$, and for all $y\in\operatorname{\mathbb{C}}^{2m}$, $\lVert S_{1}[F_{M},F_{\hat{M}}]y\rVert_{2}^{2}\in(1\pm 1/32)\lVert[F_{M},F_{\hat{M}}]y\rVert_{2}^{2}$. In particular this yields that for all $W\in\operatorname{\mathbb{C}}^{m\times m}$:

[TABLE]

Similarly, fixing $S_{1}$, with probability $\geq 1-\delta$ over the random choice of $S_{2}$:

[TABLE]

Combining these two bounds, for any $M$, with probability $\geq 1-2\delta$, for all $W\in\operatorname{\mathbb{C}}^{m\times m}$:

[TABLE]

Setting $\delta=\frac{1}{200\cdot|N|^{m}}$ gives that, by a union bound, with probability $\geq 99/100$: (13) holds for all $M\in N^{m}$ and $\max(s_{1},s_{2})\leq cm\log m\log\left(m\cdot|N|^{m}\right)=2cm^{2}\log m\log(|N|)$ for sufficiently large $c$. Plugging (13) and (12) back into (5.2.2) and applying a union bound gives that with probability $\geq 97/100$, (10) holds for all $M\in N^{m}$ and $W\in\operatorname{\mathbb{C}}^{m\times m}$ simultaneously, completing the lemma. ∎

With Lemmas 5.6, and 5.7 in place, we are ready to give our method, detailed in Algorithm 3, and its analysis.

Theorem 2.10 (Sparse Fourier transform sample complexity – full theorem).

Consider PSD Toeplitz $T\in\operatorname{\mathbb{R}}^{d\times d}$. For any rank $k$, Algorithm 3 run on inputs drawn from $\mathcal{N}(0,T)$ with $c_{1},c_{2}$ sufficiently large, rank $m\geq c_{3}k$, $n\geq c_{4}\max\left(\frac{\log\left(\frac{d}{\varepsilon}\right)\log d}{\varepsilon^{2}},k\right)$, and discretization $\alpha=\frac{\varepsilon^{2}}{c_{5}d^{3.5}}$ for sufficiently large constants $c_{3},c_{4},c_{5}$ with probability $\geq 19/20$, (1): returns $\overline{T}$ that satisfies

[TABLE]

and (2): has total sample complexity $O\left(k^{3}\log k\cdot\log\left(\frac{d}{\varepsilon}\right)+\frac{k^{2}\log k\cdot\log d\cdot\log^{2}\left(\frac{d}{\varepsilon}\right)}{\varepsilon^{2}}\right)$.

While this above error guarantee appears somewhat unusual, its can provide fairly strong bounds. See Theorem 5.8 for an example giving error $\epsilon\lVert T\rVert_{2}$ when $T$ has low stable rank $\mathrm{tr}(T)/\left\lVert T\right\rVert_{2}$.

Proof.

We first note that in Line 4, we have (see, e.g. [FT07]):

[TABLE]

If $c_{1},c_{2}$ are sufficiently large, for our setting of $\varepsilon,\delta$, by Lemma 5.7 we thus have with probability $97/100$: $\max(s_{1},s_{2})\leq cm^{2}\log m\log(1/\alpha)$ for some constant $c$ and:

[TABLE]

Further, by Lemma 5.6 for our setting of $m$, $n$, and $\alpha$, with probability $\geq 99/100$ we have:

[TABLE]

Thus by a union bound, with probability $\geq 96/100$,

[TABLE]

Finally, for our setting of $n\geq\frac{c_{4}\log\left(\frac{d}{\varepsilon}\right)\log d}{\varepsilon^{2}}$ by Lemma 5.2 and another union bound, with probability $\geq 95/100$:

[TABLE]

Finally, we note that by our bound on $\max(s_{1},s_{2})$ the entry sample complexity is bounded by

[TABLE]

which completes the theorem. ∎

Finally, we can consider the case when $T$ has low stable rank: $\frac{\mathrm{tr}(T)}{\lVert T\rVert_{2}}\leq s$. In this case $\lVert T-T_{k}\rVert_{2}\leq\frac{\mathrm{tr}(T)}{k}\leq\frac{s\lVert T\rVert_{2}}{k}$. Setting $k=\frac{cs^{2}}{\varepsilon^{2}}$ for large enough $c$ we can apply Theorem 2.10 to give:

Theorem 5.8** (Toeplitz Covariance Estimation via Sparse Fourier Transform – Low Stable Rank).**

Consider PSD $T\in\operatorname{\mathbb{R}}^{d\times d}$ with stable rank $\frac{\mathrm{tr}(T)}{\lVert T\rVert_{2}}\leq s$. Algorithm 3 run on inputs drawn from $\mathcal{N}(0,T)$ with $c_{1},c_{2}$ sufficiently large, rank $m\geq\frac{c_{3}s^{2}}{\varepsilon^{2}}$, $n\geq c_{4}\max\left(\frac{\log\left(\frac{d}{\varepsilon}\right)\log d}{\varepsilon^{2}},m\right)$, and discretization $\alpha=\frac{\varepsilon^{2}}{c_{5}d^{3.5}}$ for sufficiently large constants $c_{3},c_{4},c_{5}$ probability $\geq 19/20$, (1): returns $\overline{T}$ that satisfies

[TABLE]

and (2): has entrywise sample complexity $O\left(\frac{s^{6}}{\varepsilon^{6}}\log^{2}\left(\frac{d}{\varepsilon}\right)+\frac{s^{4}}{\varepsilon^{4}}\log^{4}\left(\frac{d}{\varepsilon}\right)\right)$.

6 Lower Bounds

In this section we give lower bounds for Toeplitz covariance estimation, demonstrating that polynomial dependencies on the rank are unavoidable in the TSC of the problem. This validates our low-rank bounds in Section 5, and also, when we take the rank to be linear in the dimension, demonstrates that a linear dependence on the dimension is unavoidable.

All lower bounds in this section will apply to the case when the matrix is in fact circulant. When $T$ is a symmetric circulant matrix, it is not hard to show that (with the possible exception of the all-ones eigenvector), all of its eigenvectors come in conjugate pairs, and the corresponding eigenvalues for these eigenvectors are the same. More concretely, one can show that for all $j=1,\ldots,\lfloor d/2\rfloor$, the following two vectors are orthonormal eigenvectors of $T$ with the same eigenvalue:

[TABLE]

and the matrix

[TABLE]

is an orthonormal matrix that diagonalizes $T$. For simplicity of exposition in the remainder of this section we will assume that $d$ is odd, however, the ideas immediately transfer over to the case where $d$ is even. For any subset $S\subseteq\{1,\ldots,\lfloor d/2\rfloor\}$, let $D_{S}\in\operatorname{\mathbb{R}}^{d\times d}$ denote the diagonal matrix which has diagonal entries $(D_{S})_{2j,2j}=(D_{S})_{2j+1,2j+1}=1$ for all $j\in S$, and [math] elsewhere, and let $\Gamma_{S}=\Phi D_{S}\Phi$. This matrix is PSD, circulant, and has eigenvalues which are either [math] and $1$, and the nonzero eigenvalues have corresponding eigenvectors exactly $u_{j},v_{j}$ for $j\in S$.

6.1 Non-adaptive lower bound in the noiseless case

In this section we consider the case where the algorithm is non-adaptive, and $T$ is exactly rank-$k$. Note that in this setting, the Prony’s method based algorithm of Theorem 2.9 gives an upper bound of $O(k\log k)$ total samples using a non-adaptive algorithm (with ESC $O(k)$). We show that this TSC bound is tight up to the factor of $\log^{2}k$.

Theorem 6.1**.**

Let $k=\omega(1)$. Any (potentially randomized) algorithm that reads a total of $m$ entries of $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,T)$ non-adaptively where $T$ is an (unknown) rank-$k$ circulant matrix and outputs $\widehat{T}$ so that $\lVert T-\widehat{T}\rVert_{2}\leq\frac{1}{10}\lVert T\rVert_{2}$ with probability $\geq 2/3$ must satisfy $m=\Omega(k/\log k)$.

We observe that the constant $9/10$ is arbitrary here, and can be replaced with any constant larger than $1/2$, by changing the constant in the dependence between $m$ and $k$. Moreover, observe that by letting $k=\Theta(d)$, we get a nearly linear total sample complexity lower bound in the full-rank case, nearly matching Theorem 2.4. It is an interesting question for future work if this can be made tight up to constants (i.e. without log factors).

By Yao’s minimax principle, it suffices to demonstrate a distribution over $\widehat{T}$ so that any deterministic algorithm that succeeds with probability $\geq 9/10$ (over the choice of $\widehat{T}$ and the random samples $x^{(1)},\ldots,x^{(n)}$) requires $m=\Omega(k/\log k)$. Without loss of generality, assume that $n=m$: since the algorithm inspects at most $m$ entries, chosen ahead of time, clearly we may assume it takes at most $m$ full samples from the distribution. If it takes fewer samples, it can simply ignore the remaining samples. Observe that any non-adaptive algorithm for learning $T$ to error $\varepsilon$ is fully characterized by the following:

• •

subsets $S_{1},\ldots,S_{m}\subseteq[d]$, where $S_{i}$ is the set of positions that the algorithm inspects for sample $i$, and

• •

a function $f:\operatorname{\mathbb{R}}^{S_{1}}\times\ldots\operatorname{\mathbb{R}}^{S_{m}}\to\operatorname{\mathbb{R}}^{d\times d}$ so that if $x^{(1)},\ldots,x^{(m)}\sim\operatorname{\mathcal{N}}(0,T)$, with probability $\geq 2/3$

[TABLE]

We will show that the following random distinguishing problem requires $m=\Omega(k/\log k)$:

Problem 6.2**.**

We consider the following two player game with parameter $k^{\prime}$.

• •

First, Alice chooses subsets $S_{1},\ldots,S_{m}\subseteq[d]$, where $\sum_{i=1}^{m}|S_{i}|\leq m$.

• •

Then, Bob chooses two random subsets $R_{0},R_{1}$ of $[(d-1)/2]$ as follows: to form $R_{0}$, he includes every $i\in[(d-1)/2]$ with probability $O(\frac{k^{\prime}\log k^{\prime}}{d})$. Then, to form $R_{1}$, if $R_{0}=\varnothing$, he lets $R_{1}=\varnothing$. Otherwise, he chooses a random $i\in R_{0}$, and removes it.

• •

If $R_{0}=\varnothing$ or $R_{1}=\varnothing$, we say Alice succeeds.

• •

Otherwise, Bob draws $\ell\sim\mathrm{Ber}(1/2)$, and generates $x^{(1)},\ldots,x^{(m)}\sim\operatorname{\mathcal{N}}(0,\Gamma_{\ell})$, where $\Gamma_{\ell}=\Gamma_{R_{\ell}}$, for $\ell=\{0,1\}$, and gives Alice the vectors $(x^{(1)})_{S_{1}},\ldots,(x^{(m)})_{S_{m}}$, as well as $R_{0}$ and $R_{1}$.

• •

Alice can then do any arbitrary deterministic postprocessing on $(x^{(1)})_{S_{1}},\ldots,(x^{(m)})_{S_{m}}$ and $R_{0}$ and $R_{1}$ and outputs $\ell^{\prime}\in\{0,1\}$, and we say she succeeds if $\ell=\ell^{\prime}$. Otherwise, we say she fails.

Theorem 6.1 will be a simple consequence of the following two lemmata:

Lemma 6.3**.**

Let $k=\omega(1)$, and suppose there exists an (potentially randomized) algorithm, that reads $m$ entries non-adaptively from $x^{(1)},\ldots,x^{(m)}\sim\operatorname{\mathcal{N}}(0,T)$ where $T$ is an (unknown) rank-$k$ circulant matrix and outputs $\widehat{T}$ so that $\lVert T-\widehat{T}\rVert_{2}\leq\frac{1}{10}\lVert T\rVert_{2}$ with probability $\geq 2/3$. Then, Alice can succeed at Problem 6.2 where $k^{\prime}=O(k/\log k)$ with probability $\geq 1/2+1/5$ with $m$ queries.

Lemma 6.4**.**

Alice cannot succeed at Problem 6.2 with probability greater than $1/2+1/5$ unless $m=\Omega(k^{\prime})$.

Observe that by combining these two lemmata, we immediately obtain Theorem 6.1. We now prove these two lemmata in turn.

Proof of Lemma 6.3.

Our reduction will be the trivial one. Any non-adaptive learning algorithm immediately gives a routine for the distinguishing problem: simply run the learning algorithm, obtain output $\widehat{T}$, and output $\operatorname*{arg\,min}_{\ell\in\{0,1\}}\lVert\widehat{T}-\Gamma_{\ell}\rVert_{2}$. This algorithm clearly also only pays for $m$ samples. We now show that any algorithm with the guarantees as in Theorem 6.1, when transferred in this way, immediately yields a solution to Problem 6.2.

In Problem 6.2, condition on the event that $|R_{0}|\neq\varnothing,|R_{1}|\neq\varnothing,$ and $|R_{1}|<|R_{0}|\leq k/2$. By basic Chernoff bounds, we know this happens with probability $\geq 99/100$ for $k=\omega(1)$. Conditioned on this event, we know that $\lVert\Gamma_{\ell}\rVert_{2}=1$, $\lVert\Gamma_{0}-\Gamma_{1}\rVert_{2}=1$, and $\mathrm{rank}(\Gamma_{\ell})\leq k$ for $\ell\in\{0,1\}$. Thus, if this event occurs, if the samples are drawn from $\Gamma_{\ell}$, then the guarantees of the supposed algorithm imply that with probability $\geq 2/3$, it outputs $\widehat{T}$ so that $\lVert\widehat{T}-\Gamma_{\ell}\rVert_{2}\leq 1/10$. By triangle inequality, this implies that $\lVert\widehat{T}-\Gamma_{1-\ell}\rVert_{2}\geq 9/10$. Hence in this case, the distinguishing algorithm succeeds so long as the learning algorithm succeeds, so overall the distinguishing algorithm succeeds with probability at least $2/3-99/100\geq 1/2+1/5$. ∎

Proof of Lemma 6.4.

Let $Y$ denote the distribution of the entries of the samples that the algorithm looks at, if we stack them to form a single, $m$-dimensional column vector. It is not hard to see that if the samples are drawn from $\operatorname{\mathcal{N}}(0,\Gamma_{\ell})$, for $\ell\in\{0,1\}$, then $Y\sim\operatorname{\mathcal{N}}(0,\Xi_{\ell})$, where

[TABLE]

We will show that if $m=O(k^{\prime})$, with probability $99/100$ over the choice of $R_{0}$ and $R_{1}$,

[TABLE]

We can then condition on the event that both $R_{0},R_{1}\neq\varnothing$, as in this case the distributions are not well-formed (and Alice succeeds by default). By a Chernoff bound, it is not hard to see that this occurs with overwhelming probability as long as $k=\omega(1)$. This implies that conditioned on this choice of $R_{0}$ and $R_{1}$, with probability $9/10$, one cannot distinguish between $Y\sim\operatorname{\mathcal{N}}(0,\Xi_{0})$ and $Y\sim\operatorname{\mathcal{N}}(0,\Xi_{1})$. Since the algorithm is only allowed to take $Y\sim\operatorname{\mathcal{N}}(0,\Xi_{\ell})$ and then do post-processing on it to distinguish between $R_{0}$ and $R_{1}$, by data processing inequalities (see e.g. [CT12]), this implies that the algorithm cannot succeed with probability more than $1/2+1/10+1/100<1/2+1/5$ over the choice of both $R_{0},R_{1}$, and $Y$.

We proceed to prove the above total variation distance bound. It is well-known (e.g. via Pinsker’s inequality) that for any two PSD matrices $M_{1},M_{2}$, we have $d_{\mathrm{TV}}(\operatorname{\mathcal{N}}(0,M_{1}),\operatorname{\mathcal{N}}(0,M_{2}))\leq O\left(\left\lVert I-M_{1}^{-1/2}M_{2}M_{1}^{-1/2}\right\rVert_{F}\right)$. This can further be upper bounded by

[TABLE]

Instantiating this bound for our case, and using the shared block-diagonal structure of our matrices $\Xi_{0}$ and $\Xi_{1}$, we get that

[TABLE]

We will bound both terms on the RHS of (16) separately. We first bound the spectral norm. Let $\Phi_{S_{i}}$ denote the restriction of $\Phi$ to the columns in $S_{i}$. Observe that $\Phi_{S_{i}}$ has orthonormal columns and hence its row leverage scores (Def. 3.2) are all identically $|S_{i}|/d\leq m/d=O(k^{\prime}/d)$. Then, $(\Gamma_{0})_{S_{i}}$ is formed by subsampling the rows of $\Phi_{S_{i}}$ according to $R_{0}$, and taking the outer product of this matrix with itself. Since each entry is included in $R_{0}$ with probability $\frac{ck^{\prime}\log k^{\prime}}{d}$ for some constant $c$, by Claim A.1 we have that with probability $\geq 1-1/(k^{\prime})^{2}$, $\lVert\frac{d}{ck^{\prime}\log k^{\prime}}\left(\Gamma_{0}\right)_{S_{i}}-\Phi_{S_{i}}\Phi_{S_{i}}^{\top}\rVert_{2}\leq 1/10$. In particular, since the smallest (nonzero) singular value of $\Phi_{S_{i}}$ is at least the smallest singular value of $\Phi$, which is 1, we obtain that for every $i$, we have $(\Gamma_{0})_{S_{i}}\succeq\frac{9ck^{\prime}\log k^{\prime}}{10d}\cdot I$ with probability at least $1-1/(k^{\prime})^{2}$. Thus, by a union bound, if $m\leq k^{\prime}$, we obtain that with probability at least $1-1/k^{\prime}$,

[TABLE]

We now turn our attention to the second term on the RHS of (16). Let $j\in[(d-1)/2]$ be the unique index in which $R_{0},R_{1}$ differ. Then, observe that

[TABLE]

Hence

[TABLE]

Plugging (17) and (18) into (16) yields that for $k^{\prime}$ sufficiently large, with probability $\geq 99/100$, if $m=O(k^{\prime})$,

[TABLE]

which completes the proof. ∎

6.2 Adaptive lower bound in the noisy case

In this section, we consider a somewhat different setting. We allow the algorithm now to be adaptive; that is, the choice of entry it can inspect is no longer specified ahead of time, as it was in the previous section, but rather can now depend on the answers it has seen so far. However, we also ask the algorithm to solve a slightly harder question, where the covariance is still circulant, but is now only approximately low-rank. By using information theoretic techniques inspired by [HIKP12], we demonstrate that in this setting, a linear dependence on the rank is still necessary. Note that by taking the rank to be linear in $d$, this implies that the near linear sample complexity result of Theorem 2.4 is tight for full-rank matrices, up to logarithmic factors. This is in spite of that fact that the algorithm used, and in fact all or algorithms, make non-adaptive queries.

Problem 6.5**.**

Let $\alpha\in[0,1]$ be a parameter. Let $T\in\operatorname{\mathbb{R}}^{d\times d}$ be an (unknown) PSD Toeplitz matrix, and let $k\leq d$ be a known integer. Given samples $X_{1},\ldots,X_{m}\sim\operatorname{\mathcal{N}}(0,T)$, and $m$ adaptively chosen entrywise lookups into $X_{1},\ldots,X_{m}$, output $\widehat{T}$ so that with probability $\geq 1/10$, we have $\left\lVert T-\widehat{T}\right\rVert_{2}\leq\frac{1}{10}\left\lVert T\right\rVert_{2}+\alpha d\left\lVert T-T_{k}\right\rVert_{2}$, where $T_{k}=\operatorname*{arg\,min}_{\mathrm{rank}-k\ M}\lVert T-M\rVert_{2}$.

As in the previous section, the constants in the statement of Problem 6.5 are more or less arbitrary, and the results hold for any (sufficiently small) constants. Our main result in this section is:

Theorem 6.6**.**

Let $\alpha\geq 10/k$. Any algorithm for Problem 6.5 requires

[TABLE]

In particular, we note a few interesting regimes.

• •

When $\alpha$ is a small constant, we get a linear lower bound, on the sample complexity, even for any algorithm achieving a weak tail guarantee of $\lVert T-\widehat{T}\rVert_{2}\leq\frac{1}{10}\lVert T\rVert_{2}+\frac{d}{10}\lVert T-T_{k}\rVert_{2}$. Notice that (up to constants in the tail), this is weaker than guarantee achieved in Theorem 2.10.

• •

When $\alpha=10/k$, this gives a super-linear lower bound of $\widetilde{\Omega}(k\log(k/d))$ for algorithms that achieve a stronger tail guarantee. In analogy to results on sparse Fourier transforms, we conjecture that this bound is in fact tight (up to log log factors).

• •

When $\alpha=O(1/k)$ and additionally $k=\Omega(d)$, this implies that any algorithm for solving Problem 2.1 requires total sample complexity $\Omega(d)$. In particular, this shows that the simple algorithm of Theorem 2.4 based on taking full samples is in fact optimal for in terms of total sample complexity, up to logarithmic factors.

Our lower bound goes through information theory, via the following game between Alice and Bob.

Problem 6.7**.**

Let $\alpha$ be as in Problem 6.5. We consider the following two-player game.

• •

First, Alice selects a subset $S\subseteq[(d-1)/2]$ of size $|S|=k/2$ uniformly at random, then generates samples $x^{(1)},\ldots,x^{(m)}\sim\operatorname{\mathcal{N}}(0,\Gamma_{S}+\frac{1}{\alpha d}I)$.

• •

Then, Bob makes $m$, sequentially adaptive queries $(i_{1},j_{1}),\ldots,(i_{m},j_{m})$. More precisely, for all times $t\leq m$, he sees $x^{(i_{t})}_{j_{t}}$, then decides $(i_{t+1},j_{t+1})$ as a deterministic function of the answers he has seen so far.

• •

After $m$ queries, Bob must output $S^{\prime}$ as a deterministic output of the answers he has seen so far. We say that Bob succeeds if $S^{\prime}=S$, otherwise we say he fails.

As before, the theorem will be a simple consequence of the following two lemmata along with Yao’s minimax principal:

Lemma 6.8**.**

Let $\alpha\geq 10/k$. Suppose there exists an algorithm for Problem 6.5 that succeeds with probability $\geq 1/10$ with $m$ entrywise lookups. Then, there exists a strategy for Bob to succeed at Problem 6.7 with probability $\geq 1/10$ with the same number of samples.

Lemma 6.9**.**

Let $\alpha\geq 10/k$. Any deterministic strategy for Bob for Problem 6.7 that succeeds with probability $\geq 1/10$ requires

[TABLE]

We now prove these two lemmata in turn.

Proof of Lemma 6.8.

The reduction is again straightforward: given an algorithm for Problem 6.5 and an instance of Problem 6.7, we run the algorithm to obtain some matrix $\widehat{T}$, and we output $\operatorname*{arg\,min}_{|S|=k/2}\lVert\widehat{T}-\Gamma_{S}\rVert_{2}$. Let $S^{*}$ be the true subset, and suppose the algorithm succeeds. For all $S$ with $|S|=k/2$, we have $\Gamma_{S}$ is rank-$k/2$, and $\lVert\Gamma_{S}+\frac{1}{\alpha d}I\rVert_{2}=1+\frac{1}{\alpha d}\leq 1+1/10$. Thus, if the algorithm succeeds, we have $\lVert\widehat{T}-(\Gamma_{S^{*}}+\alpha I)\rVert_{2}\leq\frac{1}{5}+\frac{1}{\alpha d}\ll 1/2$ by our choice of $\alpha$. Since $\lVert\Gamma_{S}-\Gamma_{S^{\prime}}\rVert_{2}=1$ for all $S\neq S^{\prime}$, we conclude that $\operatorname*{arg\,min}_{|S|=k/2}\lVert\widehat{T}-\Gamma_{S}\rVert_{2}=S^{*}$, in this case, so the algorithm successfully solves the distinguishing problem.

∎

Before we prove Lemma 6.9, we first need the following moment bound.

Lemma 6.10**.**

Let $M\in\operatorname{\mathbb{R}}^{d\times d}$ have rows $M_{1},\ldots,M_{d}$, and let $s=\max_{j=1,\ldots,d}\left\lVert M_{j}\right\rVert_{2}^{2}$. Let $X_{1},\ldots,X_{m}\sim\operatorname{\mathcal{N}}(0,MM^{\top})$. Then

[TABLE]

Proof.

For all $i\in[m]$ and $j\in[d]$, we have that $X_{i,j}\sim\operatorname{\mathcal{N}}(0,\lVert M_{j}\rVert_{2}^{2})$. Thus, by standard Gaussian concentration, there exists some universal constant $C$ so that for all $t>0$,

[TABLE]

Thus, by a union bound, and since probability is bounded by $1$, we know that for all $t>0$, we have

[TABLE]

Hence

[TABLE]

for some other constant $C^{\prime}$, where (a) follows since by explicit calculation, we have that for all $t\geq 2s\log(md)/C$, we have $md\leq\exp\left(\frac{Ct}{2s}\right)$. ∎

Proof of Lemma 6.9.

We will assume familiarity with the basics of information theory (see e.g. [CT12]). We will use $I(X;Y)$ to denote mutual information, $H$ to denote discrete entropy, and $h$ to denote differential entropy.

Suppose that there exists a protocol for Bob that succeeds with probability $\geq 9/10$. Since there exist $\binom{(d-1)/2}{k}$ possibilities for $S,$ by Fano’s inequality, we have $H(S|S^{\prime})\leq 1+\frac{9}{10}\log\binom{(d-1)/2}{k/2}$. Thus,

[TABLE]

The main work will be to show a corresponding upper bound, namely,

[TABLE]

Combining (19) and (20) and solving for $m$ immediately yields Lemma 6.9. Thus to complete the proof it suffices to verify Equation 20. WLOG we may assume that $(i_{t},j_{t})\neq(i_{t^{\prime}},j_{t^{\prime}})$ for all $t\neq t^{\prime}$, as otherwise that look-up give exactly no additional information. Observe that since $x^{(\ell)}\sim\operatorname{\mathcal{N}}(0,\Gamma_{S}+\frac{1}{\alpha d}I)$, we can write $x^{(\ell)}=y^{(\ell)}+\eta^{(\ell)}$, where $y^{(\ell)}\sim\operatorname{\mathcal{N}}(0,\Gamma_{S})$, and $\eta^{(\ell)}\sim\operatorname{\mathcal{N}}(0,\alpha I)$ are independent. Let $X^{\left(t\right)}=x^{(i_{t})}_{j_{t}}$, $Y^{\left(t\right)}=y^{(i_{t})}_{j_{t}}$ and $\eta^{\left(t\right)}=\eta^{(i_{t})}_{j_{t}}$, so that $X^{\left(t\right)}=Y^{\left(t\right)}+\eta^{\left(t\right)}$. Observe that since all the $(i,j)$ pairs are distinct, we have that each $\eta^{\left(t\right)}$ is a completely independent Gaussian. By the data processing inequality, we have that $I(S;S^{\prime})\leq I(S;X^{\left(1\right)},\ldots,X^{\left(m\right)})$. Then, for all $t=1,\ldots,m$, we have

[TABLE]

where (a) follows since $\eta^{\left(t\right)}$ is independent of all other quantities in that expression, and (b) follows by monotonicity of information.

We now seek to upper bound the variance of the random variable $X^{\left(t\right)}=x^{(i_{t})}_{j_{t}}$. The main difficulty for doing so is that $(i_{t},j_{t})$ is allowed to depend arbitrarily on $S$, and the previous $X^{\left(\ell\right)}$ for $\ell<t$. We will circumvent this by simply noting that the second moment of this random variable is certainly upper bounded by the second moment of the maximum over all possible $(i,j)$ of $\left(x^{(i)}_{j}\right)^{2}$. Observe that for any $S$, we can write $\Gamma_{S}=\Phi D_{S}(\Phi D_{S})^{\top}$, and moreover, the maximum squared $\ell_{2}$ norm of any row in $\Phi D_{S}$ is $k/d$. Hence, we obtain

[TABLE]

where (a) follows since $\eta$ is independent of $S$, and (b) follows from two applications of Lemma 6.10, and our choice of $\alpha$. Thus, by taking an expectation over $S$, we obtain that

[TABLE]

Since Gaussians maximize the entropy of a distribution with fixed variance, we can upper bound (21) by

[TABLE]

Equation (20) immediately follows from this and the chain rule for information. ∎

Conclusion

Our work provides some of the first non-asymptotic bounds on the sample complexity of recovering Toeplitz covariance matrices. We analyze several classical techniques, including those based on obtaining samples according to a sparse ruler, and Prony’s method. Additionally, we are able to improve on these techniques with new algorithms. For estimating full-rank covariance matrices, we introduce a new class of rulers with sparsity between $\sqrt{d}$ and $d$, that let us smoothly trade between the optimal entry sample complexity of methods based on $O(\sqrt{d})$-sparse rulers, and the optimal vector sample complexity of methods based on a fully dense ruler.

For estimating rank-$k$ covariance matrices for $k\ll d$, we introduce a randomly constructed sampling set, which can be used with a new recovery algorithm to achieve entry sample complexity and total sample complexity that scales just logarithmically with $d$. This surpasses the $O(\sqrt{d})$ limitation for methods based on sparse rulers, while offering significantly improved robustness in comparison to Prony’s method: the method works even when $T$ is only close to rank-$k$.

Technically, our work requires a combination of classical tools from harmonic analysis with techniques from theoretical computer science and randomized numerical linear algebra. We hope that it can help initiate increased collaboration and cross-fertilization of ideas between these research areas and the signal processing community.

Acknowledgements

We thank Eric Price for helpful discussion and clarification on [CKPS16]. We also thank Haim Avron for helpful conversations on the techniques used in Section 5.

Appendix A Leverage Score Facts and Properties

We begin by proving the alternative characterizations of leverage scores from Fact 3.3 and Fact 3.4.

Proof of Fact 3.4.

The minimization problem is simply an underconstrained regression problem and accordingly, the least norm solution can be obtained by setting $y=a_{j}(A^{*}A)^{+}A^{*}$. This yields

[TABLE]

proving equivalence to our original Definition 3.2. ∎

Fact 3.3 implies that upper bounds on $\tau_{j}(A)$ can be obtained by exhibiting a linear combination of rows in $A$ that reconstructs $a_{j}$. The squared $\ell_{2}$ norm of the coefficients in this linear combination upper bounds the leverage score. More generally, when $a_{j}$ “aligns well” with other rows in $A$, its leverage score is smaller, indicating that it is less “unique”.

Proof of Fact 3.4.

Our proof follows the one for continuous operators in [AKM*+*17]. For $j\in[d]$, let $m=\max_{y\in\operatorname{\mathbb{C}}^{s}}\frac{\left|\left(Ay\right)_{j}\right|^{2}}{\lVert Ay\rVert_{2}^{2}}$. We want to establish that for row vector $a_{j}\in\operatorname{\mathbb{C}}^{d\times 1}$,

[TABLE]

To do so, we will argue separately that $m\geq\tau_{j}(A)$ and $m\leq\tau_{j}(A)$. To see the former, we simply plug in $y=(A^{*}A)^{+}a_{j}^{*}$ to the maximization problem to obtain:

[TABLE]

To show that $m\leq\tau_{j}(A)$ we first note that we can parameterize our maximization problem:

[TABLE]

Since $(A^{*}A)^{+}$ is PSD and thus has a Hermitian square root $Z=Z^{*}$ with $(A^{*}A)^{+}=Z^{*}Z$, we can write $a_{j}(A^{*}A)^{+}w=(a_{j}Z^{*})(Zw)$ and apply Cauchy-Schwarz inequality. Specifically, for any $w\in\operatorname{\mathbb{C}}^{s}$, we have that:

[TABLE]

So we conclude that $m\leq\tau_{j}(A)$, which completes the proof. ∎

Fact 3.4 is “dual” to Fact 3.3. It implies that the leverage score of a row $a_{j}$ is higher if we can find a vector $y$ whose inner product with $a_{j}$ has large squared magnitude in comparison to its inner product with other rows in the matrix (the sum of which is $\|Ay\|_{2}^{2}$).

A.1 Matrix concentration sampling bound

We restate a matrix concentration result for sampling rows of a matrix by leverage scores, which is by now standard. Note that the result holds when rows are sampled according to any set of probabilities that upper bound the leverage scores. However, the number of rows sampled depends on the sum of probabilities used.

Claim A.1** (Leverage score sampling – e.g., Lemma 4 in the arXiv version of [CLM*+*15]).**

Given a matrix $A\in C^{d\times n}$, let $\tilde{\tau}_{1},\ldots,\tilde{\tau}_{d}\in[0,1]$ be a set of values such that, for all $j\in[d]$,

[TABLE]

There exists fixed constants $c,C$ such that if we:

set $p_{i}=\min\left(1,\tilde{\tau}_{j}c\log(d/\delta)/\varepsilon^{2}\right)$ for all $j\in[d]$, and 2. 2.

construct a matrix $B\in C^{s\times n}$ by selecting each row $a_{j}$ independently with probability $p_{i}$ and adding $\frac{1}{\sqrt{p_{j}}}a_{j}$ as a row to $B$ if it is selected

then with probability $(1-\delta)$:

$(1-\varepsilon)B^{*}B\preceq A^{*}A\preceq(1+\varepsilon)B^{*}B$, and 2. 2.

$s\leq\frac{C\log(d/\delta)}{\varepsilon^{2}}\cdot\sum_{j=1}^{d}\tilde{\tau_{j}}$.

Recall that $\preceq$ denotes the PSD Loewner ordering, which equivalently implies that for all $x\in\operatorname{\mathbb{C}}^{d}$,

[TABLE]

Finally, note that we could equivalently state this claim as sampling a matrix $S\in\operatorname{\mathbb{R}}^{s\times d}$ where $S$ has a row with entry $j$ equal to $1/\sqrt{p_{j}}$, and zeros elsewhere, for every $a_{j}$ which was sampled. Then we have $B=SA$. This is how the claim is used in Algorithm 3.

While Claim A.1 is stated for real value matrices in A.1, it holds for general complex $A$ as well: it follows Corollary 5.2 in [Tro12] which is proven for complex matrices.

Appendix B Additional Proofs: Ruler Based Methods

Lemma 4.2.

For any $\alpha\in[1/2,1]$, letting $R_{\alpha}$ be defined as in Def. 4.2, we have $|R_{\alpha}|\leq 2d^{\alpha}$, and moreover:

[TABLE]

Proof.

The bound on the size of $R_{\alpha}$ is immediate. It suffices to bound the coverage coefficient. We break up the set of distances into two intervals $A_{1}=\{1,d-d^{\alpha}\}$ and $A_{2}=\{d-d^{\alpha},\ldots,d-1\}$. For every $s\in A_{1}$, there exist at least $\lfloor d^{2\alpha-1}\rfloor$ elements $r_{2}\in R^{(2)}_{\alpha}$ so that $s<r_{2}\leq s+d^{\alpha}$. For each such element, there is at least one element $r_{1}\in\operatorname{\mathbb{R}}^{(1)}_{\alpha}$ so that $r_{2}-r_{1}=s$. Hence, for every $s\in A_{1}$, we have that $|(R_{\alpha})_{s}|\geq\lfloor d^{2\alpha-1}\rfloor$. This yields

[TABLE]

We now turn our attention to $A_{2}$. We further subdivide $A_{2}$ into $A_{2}=\bigcup_{j=1}^{\lceil d^{2\alpha-1}\rceil}B_{j}$, where $B_{j}=\{d-(j-1)d^{1-\alpha}-1,\ldots,\max(d-jd^{1-\alpha},d-d^{\alpha})\}$. For each $j\in[1,\ldots,\lceil 2^{2\alpha-1}\rceil]$, there exist $j$ elements $r\in R^{(2)}_{\alpha}$ so that $r>d-(j-1)d^{1-\alpha}-1$. This implies that $|(R_{\alpha})_{s}|\geq j$ for all $s\in B_{j}$. Therefore,

[TABLE]

Combining the bounds on $A_{1}$ and $A_{2}$ yield the desired claim. ∎

Lemma 4.7.

Let $\alpha\in[1/2,1]$. For any $k\leq d$, any PSD Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, and the sparse ruler $R_{\alpha}$ defined in Definition 4.2,

[TABLE]

where $\displaystyle T_{k}=\operatorname*{arg\,min}_{\mathrm{rank}-k\ M}\lVert T-M\rVert_{F}=\operatorname*{arg\,min}_{\mathrm{rank}-k\ M}\lVert T-M\rVert_{2}$. If $T$ is rank-$k$, $\lVert T-T_{k}\rVert_{F}^{2}=\lVert T-T_{k}\rVert_{2}^{2}=0$.

Proof.

Recall $R_{\alpha}^{(1)}\vcentcolon=\{1,\ldots,d^{\alpha}\}$ and $R^{(2)}_{\alpha}=\{d,d-d^{1-\alpha},d-2d^{1-\alpha},\ldots,d-(d^{\alpha}-1)d^{1-\alpha}\}$. We can bound:

[TABLE]

since for any $x=[x_{1},x_{2}]$, letting $y=[x_{1},-x_{2}]$ we have

[TABLE]

Squaring both sides gives (22). So to prove the lemma it suffices to show that:

[TABLE]

We first show the bound for $T_{R_{\alpha}^{(1)}}$. The bound for $T_{R_{\alpha}^{(2)}}$ follows analogously. We first claim:

[TABLE]

There are $d^{1-\alpha}$ disjoint principal submatrices of $T$ that are identical to $T_{R_{\alpha}^{(1)}}$. They correspond to the sets $R_{j}=\{jd^{\alpha}+1,jd^{\alpha}+2,\ldots,(j+1)d^{\alpha}\}$ for $j\in\{0,\ldots,d^{1-\alpha}-1\}$. We can apply Lemma 4.6 to the low-rank approximation $T_{k}$ with partition $R_{0}\cup\ldots\cup R_{d^{1-\alpha}-1}$ and $\varepsilon=\frac{2k}{d^{1-\alpha}}$. Letting $S$ be the set with $|S|\leq\frac{k}{\varepsilon}\leq\frac{d^{1-\alpha}}{2}$ whose existence is guaranteed by the lemma:

[TABLE]

Since each $T_{R_{j}}$ is identical, and since $|\{0,\ldots,d^{1-\alpha}-1\}\setminus S|\geq\frac{k}{\varepsilon}=\frac{d^{1-\alpha}}{2}$ we have:

[TABLE]

where the last bound follows from (25). This gives the bound of (24) for $T_{R_{\alpha}^{(1)}}$. We now show

[TABLE]

which will complete (23). We proceed similarly to the proof of (24). We have that

[TABLE]

where the last bound again follows from (25). Putting this together with (24) yields (23) for $T_{R_{\alpha}^{(1)}}$. An identical argument shows the same bound for $T_{R_{\alpha}^{(2)}}$, completing the lemma. ∎

B.1 A lower bound for estimation with sparse rulers: Proof of Theorem 4.5

In this section we give a lower bound for Toeplitz covariance estimation via sparse ruler based measurements. We show that the bound for rulers with sparsity $\Theta(\sqrt{d})$ given in Theorem 2.6 cannot be improved by more than a logarithmic factor. We also give a more general tradeoff for more dense rulers, however, it appears to be loose. Closing this gap is an interesting open direction.

Recall that any two distributions $F,G$, we let $d_{\mathrm{KL}}(F,G)\vcentcolon=\int\log\frac{dF}{dG}dF$ denote the KL divergence between the two distributions [CT12]. Our main tool for demonstrating this lower bound will be the following classical lemma:

Lemma B.1** (Assouad’s lemma [Ass83]).**

Let $\mathscr{C}=\{D_{z}\}_{z\in\{1,-1\}^{r}}$ be a family of $2^{r}$ probability distributions. Let $\varepsilon>0$, and suppose there exist $\alpha=\alpha(\varepsilon),\beta=\beta(\varepsilon)>0$ satisfying:

• •

For all $x,y\in\{1,-1\}^{r}$, we have $d_{\mathrm{TV}}(D_{x},D_{y})\geq\alpha\left\lVert x-y\right\rVert_{1}$.

• •

For all $x,y\in\{1,-1\}^{r}$ with Hamming distance $1$, we have $d_{\mathrm{KL}}(D_{x},D_{y})\leq\beta$.

Then, no algorithm that takes samples $x^{(1)},\ldots,x^{(n)}$ from an unknown $D_{z}\in\mathscr{C}$ and outputs $\hat{D}$ can satisfy

[TABLE]

with probability $\geq 1/10$, for all $D_{z}\in\mathscr{C}$.

B.1.1 Lower bound construction

We now begin to describe our construction. For any set $S\subseteq[d/2]$, any $\eta\in(0,1/(2d))$, and any $\sigma\in\{1,-1\}^{S}$, let $a_{S,\eta}(\sigma)\in\operatorname{\mathbb{R}}^{d}$ be the vector given by

[TABLE]

Let $T_{S,\eta}(\sigma)=\mathrm{Toep}(a_{S,\eta}(\sigma))\in\operatorname{\mathbb{R}}^{d\times d}$. By our restriction on $\eta$, we know that $\frac{1}{2}I\preceq T_{S,\eta}(\sigma)\preceq 2I$ for all $\sigma$. Therefore all of these matrices are positive semidefinite Toeplitz covariance matrices. Let

[TABLE]

For simplicity of notation, in settings where $S$ and $\eta$ understood, we will often drop them from the subscripts. Their meanings will be clear from context. We first observe:

Lemma B.2**.**

For all $\sigma,\sigma^{\prime}\in\{1,-1\}^{S}$ with $\sigma\neq\sigma^{\prime}$, we have $\left\lVert\mathrm{Toep}(\sigma)-\mathrm{Toep}(\sigma^{\prime})\right\rVert_{2}\geq 2\eta$.

Proof.

Let $S^{\prime}=\{j\in S:\sigma(j)\neq\sigma^{\prime}(j)\}$. Then we have:

[TABLE]

since $|S^{\prime}|\geq 1$, from which we conclude that $\left\lVert\mathrm{Toep}(\sigma)-\mathrm{Toep}(\sigma^{\prime})\right\rVert_{2}\geq 2\eta$, as claimed. ∎

For any ruler $R\subseteq[d]$, we also let

[TABLE]

denote the distribution of samples from any element of $\mathscr{F}_{S,\eta}$ when restricted to the ruler $R$. As above, in settings where $S$ and $\eta$ are understood, we drop them from the subscripts and denote this set simply by $\mathscr{F}_{R}$. Note that as with the original covariance matrices $\mathrm{Toep}(\sigma)$ we also have that $\frac{1}{2}I\preceq(\mathrm{Toep}(\sigma))_{R}\preceq 2I$ for all $R$.

We now show that recovering a covariance in $\mathscr{F}$ to $\varepsilon$ spectral norm error implies we must recover the covariance to good Frobenius norm error. Specifically, for any $M\in\operatorname{\mathbb{R}}^{d\times d}$, let $\gamma(M)\in\{1,-1\}^{S}$ be the vector given by $\gamma(M)_{s}=\mathrm{sign}\left(\sum_{|i-j|=s}M_{i,j}\right)$. That is, for all diagonals $s\in S$, it simply takes the average of the $s$-th diagonal of $M$ and outputs the sign. We then have:

Lemma B.3**.**

Let $\sigma\in\{1,-1\}^{S}$, and let $M\in\operatorname{\mathbb{R}}^{d\times d}$ be so that $\left\lVert M-\mathrm{Toep}(\sigma)\right\rVert_{2}\leq\xi$. Moreover, assume that for all $s\in S$, we have $|R_{s}|\leq B$. Then

[TABLE]

Since all the matrices are well-conditioned, this implies that

[TABLE]

Proof.

For all $s\in S$, let $\alpha_{s}=\frac{1}{2(d-s)}\sum_{|i-j|=s}M_{i,j}$ be the average of $M$’s entries along the $s$-th diagonal. Note that for all $s\in S$, we have that

[TABLE]

Let $S^{\prime}=\{j\in S:\sigma_{s}\neq\gamma(M)_{s}\}$. By the above calculation, we know that

[TABLE]

from which we deduce

[TABLE]

or $|S^{\prime}|\leq\xi^{2}/\eta^{2}$. This implies that

[TABLE]

from which the claim follows. ∎

In particular, Lemma B.3 implies that to prove Theorem 4.5, it suffices to demonstrate that the learning problem on distributions restricted to $R$, $\mathscr{F}_{R}$ is hard. With these tools in hand, we are now ready to prove the full lower bound:

Proof of Theorem 4.5.

Suppose that there was an algorithm violating the theorem. We will say that, on samples from $\operatorname{\mathcal{N}}(0,T)$, the algorithm succeeds if it outputs $\tilde{T}$ so that $\|\tilde{T}-T\|_{2}<\varepsilon\lVert T\rVert_{2}$.

By assumption, we have that $|R\times R|\leq d^{2\alpha}$. This implies that there exist $\geq 3d/4$ elements $s\in[d]$ so that $|R_{s}|\leq 2d^{2\alpha-1}$. In particular, this means there is a set $S\subseteq[d/2]$ of size at least $|S|\geq d/4$ so that for all $s\in S$, we have $|R_{s}|\leq 2d^{2\alpha-1}$.

Let $\eta>0$ be a parameter to be fixed later, and consider the family $\mathscr{F}=\mathscr{F}_{S,\eta}$. Since $\mathscr{A}$ succeeds with probability $\geq 1/10$ on any Toeplitz covariance $T$, it clearly succeeds with at least this probability when the Toeplitz matrix is guaranteed to come from $\mathscr{F}$.

Since the matrices in $\mathscr{F}$ have spectral norm at most $2$, by Lemma B.3, this implies that there is an algorithm with the following guarantee: given $n$ samples $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,\mathrm{Toep}(\sigma)_{R})$ for some $\mathrm{Toep}(\sigma)_{R}\in\mathscr{F}_{R}$, the algorithm outputs $\widetilde{M}\in\operatorname{\mathbb{R}}^{R\times R}$ so that with probability $\geq 1/10$, we have $d_{\mathrm{TV}}(\operatorname{\mathcal{N}}(0,\mathrm{Toep}(\sigma)_{R}),\operatorname{\mathcal{N}}(0,\widetilde{M}))\leq cd^{\alpha-1/2}\varepsilon$, for some small constant $c$.

We now show that this contradicts Assouad’s lemma (Lem. B.1) unless the number of samples taken satisfies $n=\Omega(d^{3-4\alpha}/\varepsilon^{2})$. Indeed, to do so it suffices to compute the constants $\alpha,\beta$ in the statement of the lemma. In general, for any $\sigma,\sigma^{\prime}\in V$ so that $\sigma\neq\sigma^{\prime}$, we have

[TABLE]

Since our matrices are well-conditioned, it follows that for all $\sigma,\sigma^{\prime}\in V$, we have

[TABLE]

Combining these bounds, this implies that for all $\sigma,\sigma^{\prime}\in V$, we have

[TABLE]

and if $\sigma,\sigma^{\prime}$ are neighboring, then we have

[TABLE]

Thus, by Assouad’s lemma (Lemma B.1) we know that no algorithm can can take $n$ samples $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,\mathrm{Toep}(\sigma)_{R})$ for some $\mathrm{Toep}(\sigma)_{R}\in\mathscr{F}_{R}$, and output $\widetilde{M}\in\operatorname{\mathbb{R}}^{R\times R}$ so that with probability $\geq 1/10$,

[TABLE]

Setting $\eta=O(\varepsilon d^{\alpha-1})$ gives that no algorithm can achieve

[TABLE]

However, this contraticts our assumption, which gives an algorithm achieving error $cd^{\alpha-1/2}\varepsilon$ with probability $\geq 1/10$, unless $n=\Omega(d^{3-4\alpha}/\varepsilon^{2})$. This completes the proof.

∎

Appendix C Additional Proofs: Fourier Methods

Proof of Claim 2.3.

For diagonal $T$, let $E$ denote the empirical covariance $E=\frac{1}{n}\sum_{j=1}^{n}(x^{(j)})(x^{(j)})^{T}$. The $\ell^{\text{th}}$ diagonal of each sample outerproduct $(x^{(j)})(x^{(j)})^{T}_{\ell,\ell}$ is the square of a Gaussian random variable with mean [math] and variance $T_{\ell,\ell}$ – i.e., a chi-squared random variable with mean $T_{\ell,\ell}$. $E_{\ell,\ell}$ is the average of $n$ chi-squared random variables. So by a standard chi-squared tail bound [Wai19],

[TABLE]

as long as $n\geq c\log(d/\delta)/\varepsilon^{2}$ for a fixed constant $c$. From a union bound, it follows that, for all $\ell$, $E_{\ell,\ell}\in(1\pm\varepsilon)T_{\ell,\ell}$ with probability $(1-\delta)$, and thus,

[TABLE]

The argument for circulant matrices is essentially the same, but requires a transformation. In particular, since any positive semidefinite circulant matrix $T$ can be written as $T=FDF^{*}$ where $F$ is the unitary discrete Fourier transform matrix and $D$ is a positive diagonal matrix, we see that:

[TABLE]

So it is natural to approximate $D$ by the empirical average $A=\mathrm{diag}\left(\frac{1}{n}\sum_{j=1}^{n}F^{*}(x^{(j)})(x^{(j)})^{T}F\right)$. We just need to bound how well $A_{\ell,\ell}$ concentrates around $D_{\ell,\ell}$. To do so, we claim that for any sample $x$, $\left[F^{*}xx^{T}F\right]_{\ell,\ell}$ is still distributed as a chi-squared random variable. In particular,

[TABLE]

So $\left[F^{*}xx^{T}F\right]_{\ell,\ell}$ can be written as $g^{2}$ where $g$ is a weighted sum of entries in $x_{w}$. Since $x_{w}$ is a normal random vector, it follows that $g$ is normal (even though $x_{w}$ has correlated entries). Accordingly, $\left[F^{*}xx^{T}F\right]_{\ell,\ell}$ is a chi-squared random variable and, as we already argued, has mean $T_{\ell,\ell}$.

Again from chi-squared concentration and a union bound, we have that $\|\mathrm{diag}(A)-D\|_{2}\leq\varepsilon\|D\|_{2}$ with probability $(1-\delta)$ as long as $n\geq c\log(d/\delta)/\varepsilon^{2}$. Since $F$ is unitary, this implies that $\|F\mathrm{diag}(A)F^{*}-FDF^{*}\|_{2}\leq\varepsilon\|FDF^{*}\|_{2}$, which proves the claim. ∎

Lemma 5.3 (Frequency-based low-rank approximation).

For any PSD Toeplitz matrix $T\in\operatorname{\mathbb{R}}^{d\times d}$, rank $k$, and $m\geq ck$ for some fixed constant $c$, there exists $M=\{f_{1},\ldots,f_{m}\}\subset[0,1]$ such that, letting $F_{M}\in\operatorname{\mathbb{C}}^{d\times m}$ be the Fourier matrix with frequencies $M$ (Def. 3.1) and $Z=F_{M}^{+}T^{1/2}$, we have 1): $\lVert F_{M}^{+}\rVert_{2}^{2}\leq\frac{2}{\beta}$ and 2):

[TABLE]

Proof.

Let $\overline{T}=T+\frac{\beta\lVert T\rVert_{2}}{d}\cdot I$ and $\overline{T}=F_{S}DF_{S}^{*}$ be the Vandermonde decomposition of $\overline{T}$ of Lemma 3.1 with $S=\{f_{1},\ldots f_{d}\}$. We can bound the entries of $D$ since, letting $(F_{S})_{j}$ denote the $j^{th}$ column of $F_{S}$:

[TABLE]

so $D_{j,j}\leq\lVert\overline{T}\rVert_{2}/d$. Note that $\overline{T}$ is full rank and thus $F_{S}$ is full-rank. For any $M\subseteq S$ we have $\displaystyle\sigma_{\min}(F_{M})=\min_{x:\lVert x\rVert_{2}=1}\lVert F_{M}x\rVert_{2}^{2}\geq\min_{x:\lVert x\rVert_{2}=1}\lVert F_{S}x\rVert_{2}^{2}=\sigma_{\min}(F_{S}).$ Thus,

[TABLE]

We now proceed to prove (30) and (31). Let $\overline{Z}=F_{M}^{+}\overline{T}^{1/2}$. For any $M$:

[TABLE]

Similarly,

[TABLE]

To prove (30) we will show that there exists $M$ with:

[TABLE]

Combined with (C) and (C), (34) gives $\lVert F_{M}Z-T^{1/2}\rVert_{F}\leq\sqrt{\frac{3}{2}\lVert T^{1/2}-T^{1/2}_{k}\rVert_{F}}+3\sqrt{\beta}\lVert T^{1/2}\rVert_{2}$. (30) then follows from the AMGM inequality. Following a similar argument, we note that for any $M$:

[TABLE]

and

[TABLE]

To prove (31) we thus must show that there exists $M$ with

[TABLE]

(31) follows from combining (37) with (C), (35) and (36) and applying AMGM.

We proceed to prove (34) and (37), which in combination complete the lemma. Let $U_{1}=\overline{T}^{-1/2}F_{S}D^{1/2}$. We have $U_{1}U_{1}^{*}=\overline{T}^{-1/2}F_{S}DF_{S}^{*}\overline{T}^{-1/2}=I$ so $U_{1}$ is unitary. Additionally, for any $M$,

[TABLE]

So to show the lemma, it suffices to prove (34) and (37) for lefthand side $F_{S}D^{1/2}-F_{M}(F_{M}^{+}(F_{S}D^{1/2}))$, which is just the original the lefthand side with the rotation $U_{1}$ applied. If $F$ were real, the existence of column submatrix $F_{M}$ with $m\leq ck$ for some fixed $c$ satisfying (34) is shown e.g. in Theorem 1.1 of [GS12]. The existence of $F_{M}$ with $m\leq ck$ satisfying (37) is a corollary of Theorem 27 of [CEM*+*15]. Simply taking the union of these two column subsets yields a subset satisfying both guarantees simultaneously since projection to a larger set can only reduce spectral and Frobenius norm respectively. In our setting $F_{S}$ is complex, however, the proofs of the above results can be seen to still hold. Alternatively, we can prove the lemma directly from the existing results for real matrices by noticing that, by Lemma 3.1, the frequencies of $F_{S}$ come in conjugate pairs with equal corresponding entries in $D$. Thus $F_{S}$ can be rotated via a $2\times 2$ block diagonal matrix to a real matrix $F^{\prime}$ to which these results can be applied.

Specifically, as noted, it suffices to prove (34) and (37) for $F_{S}D^{1/2}-F_{M}(F_{M}^{+}(F_{S}D^{1/2}))$. Since by Lemma 3.1 the frequencies of $F_{S}$ come in conjugate pairs with equal corresponding entries in $D$, there exists unitary $U_{2}$ such that $F_{S}D^{1/2}U_{2}=F_{S}U_{2}D^{1/2}=F^{\prime}D^{1/2}$ where $F^{\prime}$ is real. In particular, if $f_{j}$ is conjugate to itself ($e^{2\pi if_{j}}=e^{-2\pi if_{j}}$) then the $j^{th}$ column of $U_{2}$ is just the $j^{th}$ standard basis vector. For $j\neq j^{\prime}$ that are conjugate to each other, on the principal submatrix corresponding to $j,j^{\prime}$, $U_{2}$ equals:

[TABLE]

Since $F^{\prime}$ is real, prior work [GS12, CEM*+*15] implies that, for some universal constant $c$, there exists a subset $M^{\prime}$ of $m^{\prime}\leq ck$ columns of $F^{\prime}$ (which we denote $F_{M^{\prime}}^{\prime}$) such that $F_{M^{\prime}}^{\prime}({F_{M^{\prime}}^{\prime}}^{+}(F^{\prime}D^{1/2}))-F^{\prime}D^{1/2}$ satisfies (34) and (37). We can then write $F^{\prime}_{M^{\prime}}$ in the span of at most $2ck$ columns of $F_{S}$ using the transformation $U_{2}$. This proves the existence of $F_{M}$ such that $F_{M}(F_{M}^{+}(F_{S}D^{1/2}))-F_{S}D^{1/2}$ satisfies (34) and (37), completing the lemma. ∎

Lemma 5.5.

Consider PSD $T\in\operatorname{\mathbb{R}}^{d\times d}$ and $X\in\operatorname{\mathbb{R}}^{d\times n}$ with columns drawn i.i.d. from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. For any rank $k$, $\varepsilon,\delta\in(0,1]$, $m\geq c_{1}k$, and $n\geq c_{2}\left(m+\log(1/\delta)\right)$ for sufficiently large $c_{1},c_{2}$, with probability $\geq 1-\delta$, there exists $M=\{f_{1},\ldots,f_{m}\}\subset[0,1]$ such that, letting $Z=F_{M}^{+}X$, $\lVert Z\rVert_{2}^{2}\leq\frac{c_{3}\cdot d^{2}\lVert T\rVert_{2}}{\varepsilon^{2}}$ for some fixed $c_{3}$ and:

[TABLE]

Proof.

Applying Lemmas 5.3 and 5.4 (with error parameters $\beta=\frac{\varepsilon^{2}}{d\cdot 12\cdot 32}$ and $\gamma=1/2$ respectively), for sufficiently large $c_{1},c_{2}$ and $m\geq c_{1}k$, $n\geq c_{2}(m+\log(1/\delta))$, with probability $\geq 1-\delta/2$ there exists $M=\{f_{1},\ldots,f_{m}\}\subset[0,1]$ such that, letting $Z=F_{M}^{+}X$,

[TABLE]

Using that for any $A,B$, $\lVert AA^{*}-BB^{*}\rVert_{F}\leq\lVert AB^{*}-BB^{*}\rVert_{F}+\lVert AA^{*}-AB^{*}\rVert_{F}\leq\lVert A-B\rVert_{2}\cdot(\lVert A\rVert_{F}+\lVert B\rVert_{F})$ we can bound:

[TABLE]

We have $\lVert X\rVert_{F}=\lVert T^{1/2}G\rVert_{F}$ where $G$ has columns distributed as $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,I)$. By standard Gaussian concentration, for $n\geq c\log(1/\delta)$ for large enough $c$, $\lVert X\rVert_{F}\leq 2\lVert T^{1/2}\rVert_{F}$ with probability $\geq 1-\delta/2$. Conditioning on this event we also have $\lVert F_{M}Z\rVert_{F}=\lVert F_{M}F_{M}^{+}Z\rVert_{F}\leq\lVert X\rVert_{F}\leq 2\lVert T^{1/2}\rVert_{F}$ since this is just a projection of $X$ onto the column space of $F_{M}$. Plugging back into (C):

[TABLE]

Applying our stable rank bound, we have $\lVert T^{1/2}\rVert_{F}^{2}=\mathrm{tr}(T)\leq s\lVert T\rVert_{2}$. Further, $\lVert T^{1/2}-T_{k}^{1/2}\rVert_{2}^{2}\leq\frac{\lVert T^{1/2}\rVert_{F}^{2}}{k}\leq\frac{s}{k}\lVert T\rVert_{2}$. Finally, we can bound the above by:

[TABLE]

which gives the error bound. It just remains to bound $\lVert Z\rVert_{2}^{2}\leq\lVert F_{M}^{+}\rVert_{2}^{2}\cdot\lVert X\rVert_{F}^{2}\leq\frac{2}{\beta}\cdot 4\lVert T^{1/2}\rVert_{F}^{2}\leq\frac{8d\lVert T\rVert_{2}}{\beta}\leq\frac{(8\cdot 12\cdot 32)d^{2}\lVert T\rVert_{2}}{\varepsilon^{2}}$. The total success probability, union bounding over the event that a suitable $M$ exists and that $\lVert F_{M}Z\rVert_{F}\leq\lVert X\rVert_{F}\leq 2\lVert T^{1/2}\rVert_{F}$ is at least $1-\delta$. ∎

Lemma 5.6.

Consider PSD $T\in\operatorname{\mathbb{R}}^{d\times d}$ and $X\in\operatorname{\mathbb{R}}^{d\times n}$ with columns drawn i.i.d. from $\frac{1}{\sqrt{n}}\cdot\mathcal{N}(0,T)$. For any rank $k$ and $\delta,\varepsilon\in(0,1]$, consider $m\geq c_{1}k$, $n\geq c_{2}\left(m+\log(1/\delta)\right)$, and $N=\{0,\alpha,2\alpha,\ldots 1\}$ for $\alpha=\frac{\varepsilon^{2}}{c_{3}d^{3.5}}$ for sufficiently large constants $c_{1},c_{2},c_{3}$. With probability $\geq 1-\delta$, there exists $M=\{f_{1},\ldots,f_{m}\}\subset N$ with:

[TABLE]

Proof.

If suffices to show that, for the subset $M=\{f_{1},\ldots,f_{m}\}\subset[0,1]$ shown to exist (with probability $\geq 1-\delta$) in Lemma 5.5, letting $\hat{M}=\{\hat{f}_{1},\ldots,\hat{f}_{m}\}\subset N$ be the set obtained by rounding each element of $M$ to the nearest element in $N$, and letting $Z=F_{M}^{+}X$,

[TABLE]

$M$ has $\lVert F_{M}ZZ^{*}F_{M}^{*}-XX^{T}\rVert_{F}\leq 10\sqrt{\lVert T-T_{k}\rVert_{2}\cdot\mathrm{tr}(T)+\frac{\mathrm{tr}(T-T_{k})\cdot\mathrm{tr}(T)}{k}}+\frac{\varepsilon}{2}\lVert T\rVert_{2}$. Thus by triangle inequality, we will have:

[TABLE]

We thus turn to proving (39). We can bound:

[TABLE]

By Lemma 5.5, $\lVert Z\rVert_{2}^{2}\leq\frac{c\cdot d^{2}\lVert T\rVert_{2}}{\varepsilon^{2}}$ for some constant $c$. Also, $\lVert F_{M}Z\rVert_{2}\leq\lVert F_{M}F_{M}^{+}X\rVert_{F}\leq\lVert X\rVert_{F}\leq 2\lVert T\rVert_{F}\leq 2\sqrt{d}\lVert T\rVert_{2}$ with probability $\geq 1-\delta$ as argued in the proof of Lemma 5.5. Thus:

[TABLE]

To prove the lemma it thus suffices to show that $\lVert F_{\hat{M}}-F_{M}\rVert_{2}^{2}\leq\lVert F_{\hat{M}}-F_{M}\rVert_{F}^{2}\leq\frac{\varepsilon^{4}}{32cd^{3}}$. We have

[TABLE]

Since $(k-1)<d$, for $\alpha=\frac{1}{2\pi d}\cdot\frac{\varepsilon^{2}}{6\sqrt{c}d^{2.5}}$, $\left|(F_{\hat{M}}-F_{M})_{kj}\right|^{2}\leq\frac{\varepsilon^{4}}{36cd^{5}}$ and so $\lVert F_{\hat{M}}-F_{M}\rVert_{F}^{2}\leq\frac{\varepsilon^{4}}{36cd^{3}}$, giving the required bound. ∎

C.1 A priori leverage bounds for Fourier matrices

In Section 5, we apply Claim A.1 to the sample matrix $X\in\operatorname{\mathbb{R}}^{d\times n}$, which has columns drawn according to $\mathcal{N}(0,T)$. Since we do not want to examine all $d$ entries in each sample, we cannot afford to explicitly compute $X$’s leverage scores. Instead, we use that $X$ can be written as $F_{S}Y$ (exactly or approximately) for some matrix $Y$, where $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$ is a Fourier matrix with few columns. Accordingly, it suffices to bound the leverage scores of $F_{S}$.

Surprisingly this can be done without any computation: the leverage scores of a Fourier matrix $F_{S}$ with any frequencies $S=\{f_{1},\ldots f_{s}\}$ (see Definition 3.1) can be upper bounded by a simple closed-form formula which does not depend on $S$. I.e. the upper bound is valid no matter how large, small, clustered, or separated the frequencies are. Moreover, this closed form upper bound is nearly tight (up to a logarithmic factor) for any $S$. In particular, in our setting where $F_{s}$ is full rank and $s\leq d$, the true leverage scores satisfy $\sum_{j=1}^{d}\tau_{j}(F_{S})=s$ while the closed form produces leverage score upper bounds which sum to just $O(s\log s)$.

Formally, we prove two results which can be combined to obtain a nearly tight upper bound:

Lemma C.1** (Leverage score upper bounds).**

For any Fourier matrix $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$ with $s\leq d$,

[TABLE]

and, for some fixed constant $c$,

[TABLE]

For continuous Fourier operators, a bound similar to (40) was proven in [CP18] and improved in [AKM*+*19] using a different technique, which we adapt. We prove (41) using techniques from [CKPS16] which establishes a similar bound for continuous operators. These results were very recently improved by an $s$ factor in [CP19], so it might also be possible to improve our bound as well. However, polynomial improvements in (41) will only lead to constant factor improvements in our final application of Lemma C.1.

For discrete (i.e. on-grid) Fourier matrices with $s$ frequencies, the leverage score of every row is simply equal to $\frac{s}{d}$ since all columns are orthogonal and all rows have squared norm $s/d$. This bound gives an approach to establishing sparse recovery results for subsampled discrete Fourier matrices (e.g., via the restricted isometry property).

Beginning with the results of [CKPS16], bounds like Lemma C.1, on the other hand, can be used to establish similar results for off-grid Fourier matrices. In contrast to earlier work on “off-grid” sparse recovery problems (e.g., [TBSR13, CFG14, TBR15, BCG*+*15]), results based on such bounds require no assumptions on the frequencies in $S$, including no “separation” assumption that $|f_{i}-f_{j}|$ is not too small for all $i,j$. We refer the reader to [CKPS16] and [PS15] for a more in-depth discussion of this issue. In the context of our work, Lemma C.1 allows us to avoid a separation assumption on the frequencies in $T$’s Vandermonde decomposition, which in necessary in some related work [CCG15].

Before proving Lemma C.1 we state a simple corollary:

Corollary C.2**.**

For any positive integers $d$ and $s\leq d$, there is an explicit set of values $\tilde{\tau}_{1}^{(s)},\ldots,\tilde{\tau}_{d}^{(s)}\in(0,1]$ such that, for any Fourier matrix $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$,

[TABLE]

Proof.

Letting $c$ be the constant from Lemma C.1, we can satisfy the conditions of Corollary C.2 by choosing:

[TABLE]

By Lemma C.1 and the fact that $\tau_{j}(F_{s})\leq 1$ for all $s$ (by Fact 3.3) we immediately have $\tilde{\tau}_{j}^{(s)}\geq\tau_{j}(F_{s})$. So we are just left to bound $\sum_{j=1}^{d}\tilde{\tau}_{j}^{(s)}$. Let $q=\frac{cs^{6}\log^{3}(s+1)}{d}$ denote the right hand side of (41).

[TABLE]

The last inequality follows from the fact that for any integer $z$, $\log z\leq\sum_{j=1}^{z}\frac{1}{j}\leq 1+\log z$. ∎

All our applications of leverage score bounds go through Corollary C.2: we sample rows from $X$ randomly according to the closed form probabilities $\tilde{\tau}_{1},\ldots,\tilde{\tau}_{d}$. The upper bound on the sum of these probabilities allows us to bound how many samples need to be taken in expectation.

Next, we need to prove the upper bounds of (40) and (41), which are handled separately. Both rely on the following claim, which allow us to restrict our attention to $\tau_{j}(F_{S})$ for $j\leq(d+1)/2$:

Claim C.3** (Leverage score symmetry).**

For any Fourier matrix $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$ with $s\leq d$,

[TABLE]

Proof.

Let $a_{1},\ldots,a_{d}\in\operatorname{\mathbb{C}}^{s}$ denote the rows of $F_{S}$ and let $f_{1},\ldots,f_{s}$ denote the frequencies of its columns. By the minimization characterization of the leverage scores from Fact 3.3, we know that for any $j\in[d]$ there is some $y\in\operatorname{\mathbb{C}}^{d}$ with $\|y\|_{2}^{2}=\tau_{j}(F_{S})$ and $y^{T}F_{S}=a_{j}$. I.e., for all $w\in[s]$,

[TABLE]

To establish (43) refer to Definition 3.1 for the entries of $F_{S}$. Taking the complex conjugate of both sides, we have:

[TABLE]

Finally, multiplying both sides by $e^{-2\pi if_{w}d-1}$ yields:

[TABLE]

In other words, if we define $\tilde{y}\in\operatorname{\mathbb{C}}^{d}$ by setting $\tilde{y}_{g}=(y_{d+1-g})^{*}$ for $g\in[d]$ then we have $\tilde{y}^{T}F_{S}=a_{d-j+1}$. Of course $\|\tilde{y}\|_{2}^{2}=\|{y}\|_{2}^{2}$, so we conclude that $\tau_{d+1-j}(F_{S})\leq\tau_{j}(F_{S})$. Since this bound holds for all $j\in[d]$, we also have $\leq\tau_{j}(F_{S})\leq\tau_{d+1-j}(F_{S})$,, which establishes the claim. ∎

Proof of Lemma C.1.

With Claim C.3 we are ready to prove the main result: as long as we establish (40) and (41) for $j\leq(d+1)/2$ then the claim implies that it holds for all $j\in[d]$.

Proof of (40): Again let $a_{1},\ldots,a_{d}\in\operatorname{\mathbb{C}}^{s}$ denote the rows of $F_{S}$ and suppose $S=\{f_{1},\ldots,f_{s}\}$. Let $F_{S}^{(1,j)}\in\operatorname{\mathbb{C}}^{j\times s}$ contain the first $j$ rows of $F_{S}$. As discussed in Definition 3.2, $\sum_{z=1}^{j}\tau_{z}(F_{S}^{(1,j)})=\mathrm{rank}(F_{S}^{(1,j)})\leq s$. So it must be that there exists some $m\in[j]$ such that $\tau_{m}(F_{S}^{(1,j)})\leq\frac{s}{j}$. I.e. there is some integer $m\leq j$ whose leverage score in $F_{S}^{(1,j)}$ is upper bounded by exactly the leverage score upper bound that we want to prove for $j$.

To take advantage of this fact, we use Fact 3.3, which implies that since $a_{m}$’s leverage score is small, it can be written as a linear combination of $a_{1},\ldots,a_{j}$ with small coefficients. We can “shift” this linear combination to write $a_{j}$ as a linear combination of $a_{1+j-m},\ldots,a_{2j-m}$ with the same coefficients. Since $j\leq\frac{d}{2}$ and thus $2j-m\leq d$, this is a valid linear combination of rows in $F_{S}$ to form $a_{j}$, which implies an upper bound on $\tau_{j}(F_{S})$.

Formally, Fact 3.3 implies that there is some $y\in\operatorname{\mathbb{C}}^{j}$ with $\|y\|_{2}^{2}\leq\frac{s}{j}$ and $a_{m}=y^{T}F_{S}^{(1,j)}$. I.e., $\forall\,w\in[s]$,

[TABLE]

Then we claim that for any $w\in[s]$ we also have:

[TABLE]

The second equality follows from (44). Now let $\tilde{y}\in\operatorname{\mathbb{C}}^{d}$ be a vector which is zero everywhere, except that for $z\in[j]$ we set $\tilde{y}_{z+j-m}={y}_{z}$777Note that this operation is valid for $j\leq(d+1)/2$ because then $z+j-m\leq 2j-m\leq d$ (since $m\geq 1)$. From (46) we have that $\tilde{y}^{T}F_{s}=a_{j}$. Furthermore, $\|\tilde{y}\|_{2}^{2}=\|y\|_{2}^{2}$, and thus $\|\tilde{y}\|_{2}^{2}\leq\frac{s}{j}$. By Fact 3.3, we conclude that $\tau_{j}(F_{S})\leq\frac{s}{j}$ for $j\leq(d+1)/2$, as desired.

Proof of (41): While the proof above uses the minimization characterization of the leverage scores, here we will apply the maximization characterization of Fact 3.4. It follows from this characterization that (41) can be proven by showing that, for any Fourier matrix $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$ and any $y$ which can be written as $F_{S}x$ for some $x\in\operatorname{\mathbb{C}}^{s}$,

[TABLE]

In other words, we need to establish that any Fourier sparse function (which can be written as $F_{S}x$ for some set of frequencies $S$ and coefficients $x$) cannot be too concentrated at any point on the integer grid. In particular, the squared magnitude at a point $|y_{j}|^{2}$ cannot exceed the average squared magnitude $\|y\|_{2}^{2}/d$ by more than $\tilde{O}(s^{6})$.

Tools for proving this sort of bound were recently developed in [CKPS16]. We will rely on one particular result from that paper, which we restate below in a less general form than what was originally proven:

Claim C.4** (Corollary of Claim 5.2 in the arXiv version of [CKPS16]).**

For any Fourier matrix $F_{S}\in\operatorname{\mathbb{C}}^{d\times s}$ and any $y$ that can be written $y=F_{S}x$ for some $x\in\operatorname{\mathbb{C}}^{s}$, there is a fixed constant $c$ such that, if $m=cs^{2}\log(s+1)$, then for any $\tau\in 1,\ldots,\lfloor d/2m\rfloor$ and $j\leq(d+1)/2$, there always exist $C_{1},\ldots,C_{m}\in\operatorname{\mathbb{C}}$ such that:

$|C_{z}|\leq 11$* for all $z\in[m]$.* 2. 2.

$y_{j}=\sum_{z=1}^{m}C_{z}y_{\left[j+z\tau\right]}$.

Claim C.4 implies that, if $y$ has an $s$ sparse Fourier transform, then any value $y_{j}$ can be written as a small coefficient sum of entries of $y$ on any uniform grid of size $\tilde{O}(s^{2})$ that starts at $j$. The implication is that $|y_{j}|$ cannot be that much larger than the value of $y$ on this uniform grid, which is a first step to proving (47), which asserts that $|y_{j}|$ cannot be much larger than the average absolute value of all entries in $y$. Formally, from Cauchy-Schwarz inequality we have that, for any $\tau\in 1,\ldots,\lfloor d/2m\rfloor$,

[TABLE]

It follows that:

[TABLE]

Since $m=O(s^{2}\log(s+1))$, (41) follows. ∎

C.2 Prony’s method with inexact root-finding

In this section, we give algorithms that recover rank-$k$ Toeplitz covariance matrices, but when we only have access to an inexact root-finding algorithm. We will use the following standard subroutine for approximately finding roots of complex univariate polynomials as a black-box.

Lemma C.5** ([Pan02]).**

For all integers $k$, all $\beta\geq k\log k$, there exists an algorithm $\mathtt{FindRoots}$ which takes in a degree-$k$ polynomial $p:\mathbb{C}\to\mathbb{C}$ so that if $z_{1},\ldots,z_{k}$ are the roots of $p$, and satisfy $|z_{i}|\leq 1$ for all $i=1,\ldots,k$, then $\mathtt{FindRoots}(q,k,\beta)$ returns $z_{1}^{*},\ldots,z_{k}^{*}$ so that there exists a permutation $\pi:[k]\to[k]$ so that $|z_{\pi(i)}-z_{i}^{*}|\leq 2^{2-\beta/k}$. Moreover, the algorithm runs in time $O(k\log^{2}k\cdot(\log^{2}k+\log\beta))$.

We will also require the following bound on how large coordinates of Fourier sparse functions can get. This can be viewed as a discrete analog of Lemma 5.5 in [CKPS16].

Lemma C.6**.**

Let $z\in\operatorname{\mathbb{R}}^{d}$ be given by $z=F_{R}y$ for some $R$ with $|R|\leq m$. Then, we have

[TABLE]

Proof.

Let $m_{t}=(1+1/m)^{t}m$. We proceed by induction, by showing the stronger claim that

[TABLE]

for all $t$. The desired claim then follows by repeatedly applying this claim. The base case $t=0$ is trivial. For any $t>0$, assume that (48) holds for $[m_{t}]$. Let $j\in[(1+1/m)m_{t}]$. Then we can write $j$ as $j=j_{0}+m\tau$ for some $j_{0},\tau\in m_{t}/m$, so that $j_{0}+(k-1)\tau\in[m_{t}]$. By Lemma 5.3 in [CKPS16], we can write $z_{j}=\sum_{j=0}^{m-1}a_{j}z_{j_{0}+j\tau}$, for some coefficients $a_{j}$ satisfying $|a_{j}|\leq(2m)^{m}$. Therefore, we have

[TABLE]

where the last line follows from (40). Thus overall we have

[TABLE]

as claimed. ∎

For completeness, in Algorithm 4, we give the formal pseudocode for Prony’s method using FindRoots, although it is almost identical to the pseudocode in Algorithm 2. The two differences are that (1) we round the frequencies we get to a grid, and (2) we regress to find the coefficients of the frequencies, rather than exactly solving a system of linear equations.

Before we analyze the algorithm, it will be helpful to establish some additional notation. For any set $S=\{f_{1},\ldots,f_{k}\}$ of $k$ frequencies, let $\overline{S}=\{f_{1}^{\prime},\ldots,f_{k^{\prime}}^{\prime}\}$ formed by rounding the elements of $S$ to the nearest integer multiple of $2^{3-\beta/k}$. Let $\rho:[k]\to[k^{\prime}]$ be the function which takes $j\in[k]$ and maps it to the unique element $\ell\in[k^{\prime}]$ so that $|e^{2\pi if_{j}}-e^{2\pi if^{\prime}_{\ell}}|\leq 2^{3-\beta/k}$. That is, it takes each frequency to the rounded version that the algorithm recovers. For any $y\in\operatorname{\mathbb{R}}^{k}$, let $y^{\prime}\in\operatorname{\mathbb{R}}^{k^{\prime}}$ be given by $(y^{\prime})_{j}=\sum_{\ell\in\rho^{-1}(j)}y_{\ell}$. Then, putting together the guarantees of Lemmata C.5 and 5.1 yields:

Corollary C.7**.**

Let $\beta\geq k\log k$, and let $S$ be set of $k$ frequencies. Let $y\in\operatorname{\mathbb{R}}^{k}$, let $x=F_{S}y$, and let $I=\mathrm{supp}(y)$. Then, on input $x$, Algorithm 4 outputs the set of frequencies $S^{\prime}=\{e^{2\pi if_{\rho(j)}^{\prime}}:j\in I\}$, and $\hat{y}\in\operatorname{\mathbb{R}}^{|S^{\prime}|}$ so that $\lVert x-F_{S^{\prime}}\hat{y}\rVert_{2}\leq 2^{-\beta/k}d^{\left(k^{O(1)}\right)}\left\lVert y\right\rVert_{2}$. Moreover, Algorithm 4 requires only $2k$ entrywise queries to $x$.

Proof.

By Lemma 5.1, and the linear independence of columns of Fourier matrices, we know that for all $j\in I$, there exists $f\in S^{\prime}$ so that $|e^{2\pi if_{j}}-e^{2\pi if}|\leq 2^{3-\beta/k}$. Since the exact roots of the polynomial in Lemma 5.1 is exactly the set $\{e^{2\pi if_{j}}:j\in I\}$, the inexact roots that the algorithm of Lemma C.5 outputs must be rounded to the set $S^{\prime}$. Hence Algorithm 4 will output this set of frequencies.

Letting $F^{\prime}$ denote the restriction of $F_{S^{\prime}}$ to the first $2k$ rows, we have

[TABLE]

Thus, $\hat{y}$ satisfies $\left\lVert x_{[2k]}-F^{\prime}\hat{y}\right\rVert_{2}\leq O\left(2^{-\beta/k}k^{2}\left\lVert y\right\rVert_{2}\right)$. Let $z=x-F_{S^{\prime}}\hat{y}$. The vector $z$ has Fourier support of size at most $2k$, and moreover, the above calculation demonstrates that $\left\lVert z_{[2k]}\right\rVert_{2}\leq O(2^{-\beta/k}k^{2}\left\lVert y\right\rVert_{2})$. Thus, by Lemma C.6, we have that $\left\lVert x-F_{S^{\prime}}\hat{y}\right\rVert_{2}=\left\lVert z\right\rVert_{2}\leq 2^{-\beta/k}d^{\left(k^{O(1)}\right)}\left\lVert y\right\rVert_{2}$, as claimed. ∎

C.2.1 A sample-efficient and $\operatorname{poly}(d)$ time algorithm

Building on Algorithm 4, we first demonstrate how to give a simple algorithm which requires $O(k\log d/\varepsilon^{2})$ samples and runtime $\operatorname{poly}(d,k,1/\varepsilon)$ to learn a rank-$k$ Toeplitz matrix. The algorithm is very straightforward: for every sample, simply use Algorithm 4 to recover it to very good accuracy from $2k$ entries, and then run the full sample algorithm (i.e. Algorithm 1 with the full ruler) on the set of samples we recover. The formal pseudocode is given in Algorithm 5.

We prove:

Theorem C.8**.**

Let $\varepsilon>0$. Let $T\in\operatorname{\mathbb{R}}^{d\times d}$ be a PSD rank-$k$ Toeplitz matrix, let $\beta>k\log k$, and let $\eta=2^{-\beta/k}d^{\left(k^{O(1)}\right)}$. Given $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,T)$, Algorithm 5 outputs $\hat{T}$ so that

[TABLE]

Moreover, the algorithm has entrywise sample complexity $2k$, and runs in time $\operatorname{poly}(n,k,d,\log\beta)$. In particular, for any $\delta>0$, if we let $n=\Theta(\log(d/(\varepsilon\delta))\log d/\varepsilon^{2})$ and $\beta=k^{O(1)}\log d+k\log 1/\varepsilon$, then Algorithm 5 outputs $\widehat{T}$ so that $\Pr\left[\left\lVert\widehat{T}-T\right\rVert_{2}>2\varepsilon\left\lVert T\right\rVert_{2}\right]<\delta$, and the algorithm a total sample complexity of $\Theta(k\log d/\varepsilon^{2})$, and runs in time $\operatorname{poly}(k,d,1/\varepsilon,\log 1/\delta)$.

Proof.

For all $j$, we have $x^{(j)}=F_{S}D^{1/2}Ug^{(j)}$ for some $S$ with $|S|\leq k$, a fixed unitary matrix $U$, and $g^{(j)}\sim\operatorname{\mathcal{N}}(0,I)$. Let $y^{(j)}=D^{1/2}Ug^{(j)}$. By Corollary C.7, we know that $\left\lVert x^{(j)}-\hat{x}^{(j)}\right\rVert_{2}\leq\eta\left\lVert y^{(j)}\right\rVert_{2}$. Hence, we have

[TABLE]

where the last line follows since $F_{S}$ has maximum singular value at most $1$. Since $U$ is unitary, we know that $\left\lVert y^{(j)}\right\rVert_{2}^{2}$ is distributed as a non-uniform $\chi^{2}$-random variable with $d$ degrees of freedom, with coefficients given by the diagonal entries of $D$. As argued in the proof of Lemma 5.3, we know that each diagonal entry of $D$ is bounded by $\left\lVert T\right\rVert_{2}/d$. By standard concentration results for subexponential random variables, we know that for all $t>0$,

[TABLE]

for some universal constant $c$. Moreover, by Theorem 2.4, we know that

[TABLE]

for universal constants $c,C$. Thus, combining these two bounds, and observing that averaging can only decrease difference in Frobenius norm, yields (49), as claimed. ∎

C.3 An $O(\log\log d)$ time algorithm for well-conditioned Toeplitz matrices

In this section, we demonstrate an algorithm with extremely efficient runtime in terms of dimension, when the ratio between the largest and smallest non-zero singular value of the matrix is bounded away from zero. Specifically, for any rank $k$ PSD matrix $M\in\operatorname{\mathbb{R}}^{d\times d}$ with non-zero eigenvalues $\sigma_{1}\geq\ldots\geq\sigma_{k^{\prime}}>0$, where $k^{\prime}\leq d$, let $\kappa(M)=\sigma_{1}/\sigma_{k^{\prime}}$.

Our algorithm for this setting will be almost identical to the algorithm presented in the main text that assumes an exact root finding oracle. For each sample $x$, we will apply Algorithm 4 to obtain $S^{\prime}$ and $y^{\prime}$ so that $x$ is well-approximated by $F_{S^{\prime}}y^{\prime}$. We then simply estimate the diagonals of the $k\times k$ matrix given by the second moment $y$ vectors, and output the set of frequencies as well as the set of diagonal entries estimated. The formal pseudocode is presented in Algorithm 6.

Our guarantee is the following:

Theorem C.9**.**

Let $\varepsilon>0$. Let $T\in\operatorname{\mathbb{R}}^{d\times d}$ be a PSD rank-$k$ Toeplitz matrix. Given $x^{(1)},\ldots,x^{(n)}\sim\operatorname{\mathcal{N}}(0,T)$ and $\kappa(T)$, Algorithm 6 outputs $R$ and $\hat{D}$ so that if we let $\hat{T}=F_{R}\hat{D}F_{R}^{*}$, then

[TABLE]

Moreover, the algorithm has entrywise sample complexity $2k$, and runs in time

[TABLE]

In particular, for any $\delta>0$, if we let $n=\Theta(\log(1/\delta)\log k/\varepsilon^{2})$, then Algorithm 5 outputs $\widehat{T}$ so that $\Pr\left[\left\lVert\widehat{T}-T\right\rVert_{2}>2\varepsilon\left\lVert T\right\rVert_{2}\right]<\delta$, and the algorithm requires $\Theta(k\log(1/\delta)\log d/\varepsilon^{2})$ entrywise samples, and runs in time

[TABLE]

We first prove the following geometric lemma, which states that any Toeplitz matrix with bounded $\kappa$ cannot have frequencies which are too close together:

Lemma C.10**.**

Let $T\in R^{d\times d}$ be a PSD Toeplitz matrix, and let $T=F_{S}DF_{S}$ be the Vandermonde decomposition of $T$. Then

[TABLE]

Proof.

Let $j,j^{\prime}\in S$ achieve the minimum on the LHS, and WLOG assume that $D_{j,j}\geq D_{j^{\prime},j^{\prime}}$. Let $\Delta=|1-e^{2\pi i(f_{j}-f_{j^{\prime}})}|=\frac{1}{2}|e^{2\pi if_{j}}-e^{2\pi if_{j^{\prime}}}|^{2}$. By calculation, observe that

[TABLE]

As argued in Lemma 5.3, the largest eigenvalue of $T$ is at least $d\cdot D_{j,j}$. Moreover, since the vectors $\{F(f_{\ell}):\ell\in S\}$ are linearly independent and span the range of $T$, we know that the smallest non-zero eigenvalue of $T$ is given by $\min_{x\in U}\frac{x^{*}Tx}{\left\lVert x\right\rVert_{2}^{2}}$, where $U=\mathrm{span}\left(F(f_{\ell}):\ell\in S\}\right)$. Now we can write $F(f_{j})=u+v$, where $u\in\mathrm{span}\left(\{F(f_{\ell}):\ell\in S\setminus\{j^{\prime}\}\}\right)$, and $v$ is orthogonal to $F(f_{\ell})$, for all $\ell\neq j^{\prime}$. Moreover, we have that $\left\lVert u\right\rVert_{2}^{2}\geq\left|\left\langle F(f_{j}),F(f_{j^{\prime}})\right\rangle\right|$, and so by the Pythagorean theorem we must have $\left\lVert v\right\rVert_{2}^{2}\leq O(d\Delta)$, and clearly $v\in\mathrm{span}\left(\{F(f_{\ell}):\ell\in S\}\right)$. Therefore, the smallest nonzero eigenvalue of $T$ is at most

[TABLE]

Since $D_{j^{\prime},j^{\prime}}\leq D_{j,j}$, we conclude that the ratio between the largest eigenvalue of $T$ and the smallest nonzero eigenvalue of $T$ is at least $\Omega(1/(d\Delta))$, as claimed. ∎

Proof of Theorem C.9.

Let $T=F_{S}DF_{S}^{*}$ be the Vandermonde decomposition of $T$, so that $S=\{f_{1},\ldots,f_{k^{\prime}}\}$ for some $k^{\prime}\leq k$. For simplicity of exposition assume that $k^{\prime}=k$. By our choice of $\beta$, the guarantees of Lemma C.5, and Lemma C.10, we have that if $R=\{f^{\prime}_{1},\ldots,f^{\prime}_{k^{\prime\prime}}\}$ is the set of roots obtained by Algorithm 6, then $k^{\prime\prime}=|R|=k$, and moreover, for all $j\in R$, there exists a unique $\ell\in S$ so that $|e^{2\pi if^{\prime}_{j}}-e^{2\pi if_{\ell}}|\leq 2^{3-\beta/k}$. In particular, by our choice of $\beta$, we have that for all $j,j^{\prime}\in R$, $|e^{2\pi if^{\prime}_{j}}-e^{2\pi if^{\prime}_{j^{\prime}}}|\geq\Omega\left(\frac{1}{\sqrt{d\kappa(T)}}\right)$.

Let $j^{*}=\operatorname*{arg\,max}_{j\in[k]}D_{j,j}$, and let $\lambda_{1}\geq\ldots\geq\lambda_{k}$ be the nonzero eigenvalues of $T$. As argued before, we have $\lambda_{1}\geq d\cdot D_{j^{*}j^{*}}$. Moreover, since $D_{j^{*}}F_{S}F_{S}^{*}\succeq T$, we conclude that $\lambda_{k}\leq D_{j^{*}j^{*}}\sigma_{k}$, where $\sigma_{k}$ is the $k$th largest singular value of $F_{S}F_{S}^{*}$. Putting these bounds together, we conclude that $\sigma_{k}\geq\lambda_{k}/D_{j^{*}j^{*}}\geq d/\kappa(T)$. By our choice of $\beta$, we can also conclude that the smallest singular value of $F_{R}$ is at least $\Omega\left(\sqrt{\frac{d}{\kappa(T)}}\right)$.

Let $\eta=2^{-\beta/k}d^{\left(k^{O(1)}\right)}$, and let $y^{(j)}$ and $g^{(j)}$ be as before, so that $x^{(j)}=F_{S}y^{(j)}$ and $y^{(j)}=D^{1/2}Ug^{(j)}$, and $g^{(j)}$ is a standard normal Gaussian. By Corollary C.7, we know that for all $j\in[n]$, we must have $\left\lVert x^{(j)}-F_{R}\hat{y}^{(j)}\right\rVert_{2}\leq\eta\left\lVert y\right\rVert_{2}$. By a triangle inequality, this in turn implies that $\left\lVert F_{R}(y^{(j)}-\hat{y}^{(j)})\right\rVert_{2}\leq O(\eta)\cdot\left\lVert y\right\rVert_{2}$. By our bound on the smallest singular value of $F_{R}$, we conclude that $\left\lVert y^{(j)}-\hat{y}^{(j)}\right\rVert_{2}\leq O(\eta\sqrt{\kappa(T)/d})\cdot\left\lVert y\right\rVert_{2}$.

Let $Z_{\ell}=\frac{1}{n}\sum_{j\in[n]}|y^{(j)}_{\ell}|^{2}$. As argued in the proof of Theorem 2.9, we know that for any $\ell\in[k]$, we must have

[TABLE]

for some constant $c$. By standard concentration results for subexponential random variables, we know that for all $t>0$,

[TABLE]

for some universal constant $c$. Combining these two bounds with union bounds, we conclude that

[TABLE]

and thus

[TABLE]

To conclude, we note that $\left\lVert F_{S}F_{S}^{*}\right\rVert_{2}\leq\left\lVert F_{S}F_{S}^{*}\right\rVert_{F}\leq kd$, and conditioned on our good event,

[TABLE]

and hence

[TABLE]

The result then follows by adjusting $\varepsilon$, and from our choice of $\beta$. ∎