A Universal Sampling Method for Reconstructing Signals with Simple Fourier Transforms

TL;DR

This paper introduces a universal sampling method for reconstructing signals with simple Fourier structures, achieving near-optimal sample complexity and broad applicability across various signal classes.

Contribution

It formalizes the relationship between Fourier spectral constraints and sample complexity, and presents a universal, efficient sampling strategy applicable to multiple signal types.

Findings

Achieves near-optimal sample complexity based on statistical dimension.

Provides a universal sampling strategy effective for all signal classes.

Offers a simple, efficient algorithm matching state-of-the-art for bandlimited and sparse signals.

Abstract

Reconstructing continuous signals from a small number of discrete samples is a fundamental problem across science and engineering. In practice, we are often interested in signals with 'simple' Fourier structure, such as bandlimited, multiband, and Fourier sparse signals. More broadly, any prior knowledge about a signal's Fourier power spectrum can constrain its complexity. Intuitively, signals with more highly constrained Fourier structure require fewer samples to reconstruct. We formalize this intuition by showing that, roughly, a continuous signal from a given class can be approximately reconstructed using a number of samples proportional to the *statistical dimension* of the allowed power spectrum of that class. Further, in nearly all settings, this natural measure tightly characterizes the sample complexity of signal reconstruction. Surprisingly, we also show that, up to logarithmic factors, a universal non-uniform sampling strategy can achieve this optimal complexity for *any class of signals*. We present a simple and efficient algorithm for recovering a signal from the samples taken. For bandlimited and sparse signals, our method matches the state-of-the-art. At the same time, it gives the first computationally and sample efficient solution to a broad range of problems, including multiband signal reconstruction and kriging and Gaussian process regression tasks in one dimension. Our work is based on a novel connection between randomized linear algebra and signal reconstruction with constrained Fourier structure. We extend tools based on statistical leverage score sampling and column-based matrix reconstruction to the approximation of continuous linear operators that arise in signal reconstruction. We believe that these extensions are of independent interest and serve as a foundation for tackling a broad range of continuous time problems using randomized methods.

Tables

Table 1. Table 1 : Statistical dimension upper bounds for common Fourier interpolation problems. Our result (Theorem 1 ) requires O ​ ( s μ , ϵ ⋅ log ⁡ s μ , ϵ ) 𝑂 ⋅ subscript 𝑠 𝜇 italic-ϵ subscript 𝑠 𝜇 italic-ϵ O(s_{\mu,\epsilon}\cdot\log s_{\mu,\epsilon}) samples.

Fourier prior, $\mu$	Statistical dimension, $s_{\mu ,ϵ}$	Proof
$k$-sparse	$k$	Since ${𝒦}_{\mu }$ has rank $k$.
bandlimited to $[-F,F]$	$O\left(FT+\log \left(1/ϵ\right)\right)$	Theorem 48.
multiband, widths $F_1,\dots ,F_s$	$O\left({\sum }_i{F}_{i}T+s\log \left(1/ϵ\right)\right)$	Theorem 53.999Just as Theorem 48 intuitively matches the Nyquist sampling rate, Theorem 53 intuitively matches the Landau rate for asymptotic recovery of multiband functions [Lan67a].
Gaussian, variance $F$	$O\left(FT\sqrt{\log \left(1/ϵ\right)}+\log \left(1/ϵ\right)\right)$	Theorem 54.
Cauchy-Lorentz, scale $F$	$O\left(FT\sqrt{1/ϵ}+\sqrt{1/ϵ}\right)$	Theorem 55.

Equations

$$\hat{y}(\xi )={\int }_{-\infty }^{\infty }y(t)e^{-2\pi it\xi }dt.$$

$$\left[{\mathcal{F}}_{\mu }^{\ast }g\right](t)\stackrel{\mathrm{def}}{=}{\int }_{\mathbb{R}}g(\xi )e^{2\pi i\xi t}d\mu (\xi ).$$

$$\Vert y-\tilde{y}{\Vert }_T^2\le ϵ\Vert x{\Vert }_{\mu }^2+C\Vert n{\Vert }_T^2,$$

$$[{\mathcal{K}}_{\mu }z](t)\stackrel{\mathrm{def}}{=}{\int }_{\xi \in \mathbb{R}}{e}^{2\pi i\xi t}\left[\frac{1}{T}{\int }_{s\in [0,T]}z(s)e^{-2\pi i\xi s}ds\right]d\mu (\xi ).$$

$$s_{\mu ,ϵ}\stackrel{\mathrm{def}}{=}\mathrm{tr}({\mathcal{K}}_{\mu }({\mathcal{K}}_{\mu }+ϵ{\mathcal{I}}_T{)}^{-1}),$$

$$s_{\mu ,ϵ}=\sum _{i=1}^{\infty }\frac{{\lambda }_i\left({\mathcal{K}}_{\mu }\right)}{{\lambda }_i\left({\mathcal{K}}_{\mu }\right)+ϵ}.$$

$$k_{\mu }(t_1,t_2)={\int }_{\xi \in \mathbb{R}}{e}^{-2\pi i(t_1-t_2)\xi }d\mu (\xi ).$$

$$\left[{\mathcal{F}}_{\mu }f\right](\xi )=\frac{1}{T}{\int }_0^{T}f(t)e^{-2\pi it\xi }dt.$$

$$\left[{\mathcal{F}}_{\mu }^{\ast }g\right](t)\stackrel{\mathrm{def}}{=}{\int }_{\mathbb{R}}g(\xi )e^{2\pi i\xi t}d\mu (\xi ).$$

$${\mathcal{K}}_{\mu }={\mathcal{F}}_{\mu }^{\ast }{\mathcal{F}}_{\mu }.$$

$$[{\mathcal{K}}_{\mu }z](t)=\frac{1}{T}{\int }_0^T{k}_{\mu }(s,t)z(s)ds,$$

$${\mathcal{F}}_{\mu }(\xi ,t)=\sqrt{\frac{p(\xi )}{T}}\cdot e^{-2\pi it\xi }.$$

$$\underset{g\in L_2(\mu )}{\min }\Vert {\mathcal{F}}_{\mu }^{\ast }g-(y+n){\Vert }_T^2+ϵ\Vert g{\Vert }_{\mu }^2.$$

$$\Vert {\mathcal{F}}_{\mu }^{\ast }\tilde{g}-(y+n){\Vert }_T^2+ϵ\Vert \tilde{g}{\Vert }_{\mu }^2\le C\cdot \underset{g\in L_2(\mu )}{\min }\left[\Vert {\mathcal{F}}_{\mu }^{\ast }g-(y+n){\Vert }_T^2+ϵ\Vert g{\Vert }_{\mu }^2\right].$$

$$\Vert {\mathcal{F}}_{\mu }^{\ast }\tilde{g}-y{\Vert }_T^2\le 2Cϵ\Vert x{\Vert }_{\mu }^2+2(C+1)\Vert n{\Vert }_T^2.$$

$$\Vert {\mathcal{F}}_{\mu }^{\ast }\tilde{g}-y{\Vert }_T-\Vert n{\Vert }_T$$

$$\Vert {\mathcal{F}}_{\mu }^{\ast }\tilde{g}-y{\Vert }_T$$

$$\Vert {\mathcal{F}}_{\mu }^{\ast }\tilde{g}-y{\Vert }_T^2$$

$${\tau }_{\mu ,ϵ}(t)=\frac{1}{T}\cdot \underset{\{\alpha \in L_2(\mu ):\Vert \alpha {\Vert }_{\mu }>0\}}{\max }\frac{{\left|[{\mathcal{F}}_{\mu }^{\ast }\alpha ](t)\right|}^2}{\Vert {\mathcal{F}}_{\mu }^{\ast }\alpha {\Vert }_T^2+ϵ\Vert \alpha {\Vert }_{\mu }^2}.$$

$${\int }_0^T{\tau }_{\mu ,ϵ}(t)dt=s_{\mu ,ϵ}\stackrel{\mathrm{def}}{=}\mathrm{tr}({\mathcal{K}}_{\mu }({\mathcal{K}}_{\mu }+ϵ{\mathcal{I}}_T{)}^{-1}).$$

$${\tau }_{\mu ,ϵ}(t)=\frac{1}{T}\cdot ⟨{\phi }_t,({\mathcal{G}}_{\mu }+ϵ{\mathcal{I}}_{\mu }{)}^{-1}{\phi }_t{⟩}_{\mu }.$$

$${\tau }_{\mu ,ϵ}(t)=\frac{1}{T}\cdot \underset{\beta \in L_2(T)}{\min }\frac{\Vert {\mathcal{F}}_{\mu }\beta -{\phi }_t{\Vert }_{\mu }^2}{ϵ}+\Vert \beta {\Vert }_T^2.$$

$$\left[\mathbf{F}g\right](\xi )=\sum _{j=1}^s{w}_j\cdot g(j)\cdot e^{-2\pi i\xi t_j}$$

$$\tilde{g}=\underset{g\in L_2(\mu )}{\mathrm{argmin}}\left[\Vert {\mathbf{F}}^{\ast }g-(\mathbf{y}+\mathbf{n}){\Vert }_2^2+ϵ\Vert g{\Vert }_{\mu }^2\right]$$

$$\Vert {\mathcal{F}}_{\mu }^{\ast }\tilde{g}-(y+n){\Vert }_T^2+ϵ\Vert \tilde{g}{\Vert }_{\mu }^2\le 3\underset{g\in L_2(\mu )}{\min }\left[\Vert {\mathcal{F}}_{\mu }^{\ast }g-(y+n){\Vert }_T^2+ϵ\Vert g{\Vert }_{\mu }^2\right].$$

$$s_{\mu ,cϵ}\le \frac{1}{c}{s}_{\mu ,ϵ}.$$

$$\Vert \tilde{y}-y{\Vert }_T^2\le 6ϵ\Vert x{\Vert }_{\mu }^2+8\Vert n{\Vert }_T^2.$$

$$\tilde{g}=\underset{g\in L_2(\mu )}{\mathrm{argmin}}\left[\Vert {\mathbf{F}}^{\ast }g-(\mathbf{y}+\mathbf{n}){\Vert }_2^2+ϵ\Vert g{\Vert }_{\mu }^2\right].$$

$$\Vert \tilde{y}-y{\Vert }_T^2\le 6ϵ\Vert x{\Vert }_{\mu }^2+8\Vert n{\Vert }_T^2.$$

$$\tilde{g}=\mathbf{F}(\mathbf{K}+ϵ\mathbf{I}{)}^{-1}(\mathbf{y}+\mathbf{n})=\mathbf{F}(\mathbf{K}+ϵ\mathbf{I}{)}^{-1}\bar{\mathbf{y}}$$

$$\tilde{y}={\mathcal{F}}_{\mu }^{\ast }\tilde{g}$$

Full text

We formalize this intuition by showing that, roughly, a continuous signal from a given class can be approximately reconstructed using a number of samples proportional to the statistical dimension of the allowed power spectrum of that class. We prove that, in nearly all settings, this natural measure tightly characterizes the sample complexity of signal reconstruction.

Surprisingly, we also show that, up to logarithmic factors, a universal non-uniform sampling strategy can achieve this optimal complexity for any class of signals. We present a simple, efficient, and general algorithm for recovering a signal from the samples taken. For bandlimited and sparse signals, our method matches the state-of-the-art. At the same time, it gives the first computationally and sample efficient solution to a broad range of problems, including multiband signal reconstruction and kriging and Gaussian process regression tasks in one dimension.

Our work is based on a novel connection between randomized linear algebra and the problem of reconstructing signals with constrained Fourier structure. We extend tools based on statistical leverage score sampling and column-based matrix reconstruction to the approximation of continuous linear operators that arise in the signal reconstruction problem. We believe that these extensions are of independent interest and serve as a foundation for tackling a broad range of continuous time problems using randomized methods.

1 Introduction

Consider the following fundamental function fitting problem, pictured in Figure 1. We can access a continuous signal $y(t)$ at any time $t\in[0,T]$. We wish to select a finite set of sample times $t_{1},\ldots,t_{q}$ such that, by observing the signal values $y(t_{1}),\ldots,y(t_{q})$ at those samples, we are able to find a good approximation $\tilde{y}$ to $y$ over the entire range $[0,T]$. We also study the problem in a noisy setting, where for each sample $t_{i}$, we only observe $y(t_{i})+n(t_{i})$ for some fixed noise function $n$.

We seek to understand:

How many samples $q$ are required to approximately reconstruct $y$ and how should we select these samples? 2. 2.

After sampling, how can we find and represent $\tilde{y}$ in a computationally efficient way?

Answering these questions requires assumptions about the underlying signal $y$. In particular, for the information at our samples $t_{1},\ldots,t_{q}$ to be useful in reconstructing $y$ on the entirety of $[0,T]$, the signal must be smooth, structured, or otherwise “simple” in some way.

Across science and engineering, by far one of the most common ways in which structure arises is through various assumptions about $\hat{y}$, the Fourier transform of $y$:

[TABLE]

Our goal is to understand signal reconstruction under natural constraints on the complexity of $\hat{y}$.

1.1 Classical sampling theory and bandlimited signals

Classically, the most standard example of such a constraint is requiring $y$ to be bandlimited, meaning that $\hat{y}$ is only non-zero for frequencies $\xi$ with $|\xi|\leq F$ for some bandlimit $F$. In this case, we recall the famous sampling theory of Nyquist, Shannon, and others [Whi15, Kot33, Nyq28, Sha49]. This theory shows that $y$ can be reconstructed exactly using sinc interpolation (i.e, Whittaker-Shannon interpolation) if $1/2F$ uniformly spaced samples of $y$ are taken per unit of time (the ‘Nyquist rate’).

Unfortunately, this theory is asymptotic: it requires infinite samples over the entire real line to interpolate $y$, even at a single point. When a finite number of samples are taken over an interval $[0,T]$, sinc interpolation is not a good reconstruction method, either in theory or in practice [Xia01].222Approximation bounds can be obtained by truncating the Whittaker-Shannon method; however, they are weak, depending polynomially, rather than logarithmically, on the desired error $\epsilon$ (see Appendix A, Example 25).

This well-known issue was resolved through a seminal line of work by Slepian, Landau, and Pollak [SP61, LP61, LP62], who presented a set of explicit basis functions for interpolating bandlimited functions when a finite number of samples are taken from a finite interval. Their so-called “prolate spheroidal wave functions” can be combined with numerical quadrature methods [XRY01, STR06, KZW*+*17] to obtain sample efficient (and computationally efficient) methods for bandlimited reconstruction. Overall, this work shows that roughly $O(FT+\log(1/\epsilon))$ samples from $[0,T]$ are required to interpolate a signal with bandlimit $F$ to accuracy $\epsilon$ on that same interval.333We formalize our notion of accuracy in Section 2.

1.2 More general Fourier structure

While the aforementioned line of work is beautiful and powerful, in today’s world, we are interested in far more general constraints than bandlimits. For example, there is wide-spread interest in Fourier-sparse signals [Don06], where $\hat{y}$ is only non-zero for a small number of frequencies, and multiband signals, where the Fourier transform is confined to a small number of intervals. Methods for recovering signals in these classes have countless applications in communication, imaging, statistics, and a wide variety of other disciplines [Eld15].

More generally, in statistical signal processing, a prior distribution, specified by some probability measure $\mu$, is often assumed on the frequency content of $y$ [EU06, RvdVU05]. For signals with bandlimit $F$, $\mu$ would be the uniform probability measure on $[-F,F]$. Alternatively, instead of assuming a hard bandlimit, a zero-centered Gaussian prior on $\hat{y}$ can encode knowledge that higher frequencies are less likely to contribute significantly to $y$, although they may still be present. Such a prior naturally suits a Bayesian approach to signal reconstruction [HS93] and, in fact, is essentially equivalent to assuming $y$ is a stationary stochastic process with a certain covariance function (see Section 3 and Appendix G). Under various names, including “Gaussian process regression” and “kriging,” likelihood estimation under a covariance prior is the dominant statistical approach to fitting continuous signals in many scientific disciplines, from geostatistics to economics to medical imaging [Rip05, RW06].

1.3 Our contributions

Despite their clear importance, accurate methods for fitting continuous signals under most common Fourier transform priors are not well understood, even 50 years after the groundbreaking work of Slepian, Landau, and Pollak on the bandlimited problem. The only exception is Fourier sparse signals: the noiseless interpolation problem can be solved using classical methods [dP95, Pis73, BM86], and recent work has resolved the much more difficult noisy case [CKPS16, CP18].

In this paper, we address the problem far more generally. Our contributions are as follows:

We tightly characterize the information theoretic sample complexity of reconstructing $y$ under any Fourier transform prior, specified by probability measure $\mu$. In essentially all settings, we can prove that this complexity scales nearly linearly with a natural statistical dimension parameter associated with $\mu$. See Theorem 1. 2. 2.

We present a method for sampling from $y$ that achieves the aforementioned statistical dimension bound to within a polylogarithmic factor. Our approach is randomized and universal: we prove that it is possible to draw $t_{1},\ldots,t_{q}$ from a fixed non-uniform distribution over $[0,T]$ that is independent of $\mu$, i.e., “spectrum-blind.” In other words, the same sampling scheme works for bandlimited, sparse, or more general priors. See Theorem 2. 3. 3.

We show that $y$ can be recovered from $t_{1},\ldots,t_{q}$ using a simple, efficient, and completely general interpolation method. In particular, we just need to solve a kernel ridge regression problem using $y(t_{1}),\ldots,y(t_{q})$, with an appropriately chosen kernel function for $\mu$. This method runs in $O(q^{3})$ time and is already widely used for signal reconstruction in practice, albeit with suboptimal strategies for choosing $t_{1},\ldots,t_{q}$. See Theorem 3.

Overall, this approach gives the first finite sample, provable approximation bounds for all common Fourier-constrained signal reconstruction problems beyond bandlimited and sparse functions.

Our results are obtained by drawing on a rich set of tools from randomized numerical linear algebra, including sampling methods for approximate regression and deterministic column-based low-rank approximation methods [BSS14, CNW16]. Many of these methods view matrices as sums of rank-$1$ outer products and approximate them by sampling or deterministically selecting a subset of these outer products. We adapt these tools to the approximation of continuous operators, which can be written as the (weak) integral of rank-$1$ operators. For example, our universal time domain sampling distribution is obtained using the notion of statistical leverage [SS11, AM15, DM16], extended to a continuous Fourier transform operator that arises in the signal reconstruction problem. We hope that, by extending many of the fundamental contributions of randomized numerical linear algebra to build a toolkit for ‘randomized operator theory’, our work will offer a starting point for progress on many signal processing problems using randomized methods.

2 Formal statement of results

As suggested, we formally capture Fourier structure through any probability measure $\mu$ over the reals.444Formally, we consider the measure space $(\mathbb{R},\mathcal{B},\mu)$ where $\mathcal{B}$ is the Borel $\sigma$-algebra on $\mathbb{R}$. We often refer to $\mu$ as a “prior”, although we will see that it can be understood beyond the context of Bayesian inference. The simplicity of a set of constraints will be quantified by a natural statistical dimension parameter for $\mu$, defined in Section 2.1.

For signals with bandlimit $F$, $\mu$ is the uniform probability measure on $[-F,F]$. For multiband signals, it is uniform on the union of $k$ intervals, while for Fourier-sparse functions, $\mu$ is uniform on a union of $k$ Dirac measures. More general priors are visualized in Figure 2. Those based on Gaussian or Cauchy-Lorentz distributions are especially common in scientific applications, and we will discuss examples shortly. For now, we begin with our main problem formulation.

Problem 1**.**

Given a known probability measure $\mu$ on $\mathbb{R}$, for any $t\in[0,T]$, define the inverse Fourier transform of a function $g(\xi)$ with respect to $\mu$ as

[TABLE]

Suppose our input $y$ can be written as $y=\mathcal{F}^{*}_{\mu}\,x$ for some frequency domain function $x(\xi)$ and, for any chosen $t$, we can observe $y(t)+n(t)$ for some fixed noise function $n(t)$. Then, for error parameter $\epsilon$, our goal is to recover an approximation $\tilde{y}$ satisfying

[TABLE]

where $\|x\|_{\mu}^{2}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\int_{\mathbb{R}}|x(\xi)|^{2}\,d\mu(\xi)$ is the energy of the function $x$ with respect to $\mu$, while $\|z\|_{T}^{2}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\frac{1}{T}\int_{0}^{T}|z(t)|^{2}dt$, so that $\|y-\tilde{y}\|_{T}^{2}$ is our mean squared error and $\|n\|_{T}^{2}$ is the mean squared noise level. $C\geq 1$ is a fixed positive constant.

Unlike the $\|x\|_{\mu}^{2}$ term in (2), which we can control by adjusting $\epsilon$, we can never hope to recover $y$ to accuracy better than $\|n\|_{T}^{2}$. Accordingly, we consider $\|n\|_{T}^{2}$ to be small and are happy with any solution of Problem 1 that is within a constant factor of optimal – i.e., where $C=O(1)$.

Problem 1 captures signal reconstruction under all standard Fourier transform constraints, including bandlimited, multiband, and sparse signals.555For sparse or multiband signals, Problem 1 assumes frequency or band locations are known a priori. There has been significant work on algorithms that can recover $y$ when these locations are not known [ME09, Moi15, PS15, CKPS16]. Understanding this more complicated problem in the multiband case is an important future direction. The error in (2) naturally scales with the average energy of the signal over the allowed frequencies. For more general priors, $\|x\|_{\mu}^{2}$ will be larger when $y$ contains a significant component of frequencies with low density in $\mu$.666Informally, decreasing $d\mu(\xi)$ by a factor of $c>1$ requires increasing $x(\xi)$ by a factor of $c$ to give the same time domain signal. This increases $x(\xi)^{2}$ by a factor of $c^{2}$ and so increases its contribution to $\|x\|_{\mu}^{2}$ by a factor of $c^{2}/c=c$. For a given number of samples, we would thus incur larger error in (2) in comparison to a signal that uses more “likely” frequencies.

As an alternative to Problem 1, we can formulate signal fitting from a Bayesian perspective. We assume that $n$ is independent random noise, and $y$ is a stationary stochastic process with expected power spectral density $\mu$. This assumption on $y$’s power spectral density is equivalent to assuming that $y$ has covariance function (a.k.a. autocorrelation) $\hat{\mu}(t)$, which is the type of prior used in kriging and Gaussian process regression. While we focus on the formulation of Problem 1 in this work, we give an informal discussion of the Bayesian setup in Appendix G.

Examples and applications

As discussed in Section 1.2, “hard constraint” versions of Problem 1, such as bandlimited, sparse, and multiband signal reconstruction, have many applications in communications, imaging, audio, and other areas of engineering. Generalizations of the multiband problem to non-uniform measures (see Figure 2(d)) are also useful in various communication problems [ME10].

On the other hand, “soft constraint” versions of the problem are widely applied in scientific applications. In medical imaging, images are often denoised by setting $\mu$ to a heavy-tailed Cauchy-Lorentz measure on frequencies [Fud89, LH95, BWvOGD01]. This corresponds to assuming an exponential covariance function for spatial correlation. Exponential covariance and its generalization, Matérn covariance, are also common in the earth and geosciences [Rip89, Rip05], as well as in general image processing [PPV02, RVU06].

A Gaussian prior $\mu$, which corresponds to Gaussian covariance, is also used to model both spatial and temporal correlation in medical imaging [FJT94, WMN*+*96] and is very common in machine learning. Other choices for $\mu$ are practically unlimited. For example, the popular ArcGIS kriging library also supports the following covariance functions: circular, spherical, tetraspherical, pentaspherical, rational quadratic, hole effect, k-bessel, and j-bessel, and stable [ESR18].

2.1 Sample complexity

With Problem 1 defined, our first goal is to characterize the number of samples required to reconstruct $y$, as a function of the accuracy parameter $\epsilon$, the range $T$, and the measure $\mu$. We do so using what we refer to as the Fourier statistical dimension of $\mu$, which corresponds to the standard notion of statistical or ‘effective dimension’ for regularized function fitting problems [HTF02, Zha05].

Definition 2** (Fourier statistical dimension).**

For a probability measure $\mu$ on $\mathbb{R}$ and time length $T$, define the kernel operator $\mathcal{K}_{\mu}:L_{2}(T)\rightarrow L_{2}(T)$777$L_{2}(T)$ denotes the complex-valued square integrable functions with respect to the uniform measure on $[0,T]$. as:

[TABLE]

Note that $\mathcal{K}_{\mu}$ is self-adjoint, positive semidefinite and trace-class.888See Section 3 for a formal explanation of these facts. The Fourier statistical dimension for $\mu$, $T$, and error $\epsilon$ is denoted by $s_{\mu,\epsilon}$ and defined as:

[TABLE]

where $\mathcal{I}_{T}$ is the identity operator on $L_{2}(T)$. Letting $\lambda_{i}(\mathcal{K}_{\mu})$ denote the $i^{th}$ largest eigenvalue of $\mathcal{K}_{\mu}$, we may also write

[TABLE]

Note that $\mathcal{K}_{\mu}$ and $s_{\mu,\epsilon}$ as defined above, and $\mathcal{F}_{\mu}$ as defined in Problem 1 all depend on $T$ and thus could naturally be denoted $\mathcal{F}_{\mu,T}$, $\mathcal{K}_{\mu,T}$, and $s_{\mu,\epsilon,T}$. However, since $T$ is fixed throughout our results, for conciseness we do not use $T$ in our notation for these and related notions.

It is not hard to see that $s_{\mu,\epsilon}$ increases as $\epsilon$ decreases, meaning that we will require more samples to obtain a more accurate solution to Problem 1. The operator $\mathcal{K}_{\mu}$ corresponds to taking the Fourier transform of a time domain input $z(t)$, scaling that transform by $\mu$, and then taking the inverse Fourier transform. Readers familiar with the literature on bandlimited signal reconstruction will recognize $\mathcal{K}_{\mu}$ as the natural generalization of the frequency limiting operator studied in the work of Landau, Slepian, and Pollak on prolate spheroidal wave functions [SP61, LP61, LP62]. In that work, it is established that a quantity nearly identical to $s_{\mu,\epsilon}$ bounds the sample complexity of solving Problem 1 for bandlimited functions.

Our first technical result is that this is actually true for any prior $\mu$.

Theorem 1** (Main result, sample complexity).**

For any probability measure $\mu$, Problem 1 can be solved using $q=O\left(s_{\mu,\epsilon}\cdot\log s_{\mu,\epsilon}\right)$ noisy signal samples $y(t_{1})+n(t_{1}),\ldots,y(t_{q})+n(t_{q})$.

What does Theorem 1 imply for common classes of functions with constrained Fourier transforms? Table 1 includes a list of upper bounds on $s_{\mu,\epsilon}$ for many standard priors.

A complexity of $O(s_{\mu,\epsilon}\cdot\log s_{\mu,\epsilon})$ equates to $\tilde{O}(k)$ samples for $k$-sparse functions and $\tilde{O}(FT+\log 1/\epsilon)$ for bandlimited functions. Up to log factors, these bounds are tight for these well studied problems. In Section 6, we show that Theorem 1 is actually tight for all common Fourier transform priors: $\Omega(s_{\mu,\epsilon})$ time points are required for solving Problem 1 as long as $s_{\mu,\epsilon}$ grows slower than $1/\epsilon^{p}$ for some $p<1$. This property holds for all $\mu$ in Table 1. We conjecture that our lower bound can be extended to hold even without this weak assumption.

To compliment the sample complexity bound of Theorem 1, we introduce a universal method for selecting samples $t_{1},\ldots,t_{q}$ that nearly matches this complexity. Our method selects samples at random, in a way that does not depend on the specific prior $\mu$.

Theorem 2** (Main result, sampling distribution).**

For any sample size $q$, there is a fixed probability density $p_{q}$ over $[0,T]$ such that, if $q$ time points $t_{1},\ldots,t_{q}$ are selected independently at random according to $p_{q}$, and $q\geq c\cdot s_{\mu,\epsilon}\cdot\log^{2}s_{\mu,\epsilon}$ for some fixed constant $c$, then it is possible to solve Problem 1 with probability 99/100 using the noisy signal samples $y(t_{1})+n(t_{1}),\ldots,y(t_{q})+n(t_{q})$.101010In Section 5.4, we formally quantify the tradeoff between success probability and sample complexity.

Theorem 2 is our main technical contribution. By achieving near optimal sample complexity with a universal distribution, it shows that wide range of Fourier constrained interpolation problems considered in the literature are more closely related than previously understood.

Moreover, $p_{q}$ (which is formally defined in Theorem 17) is very simple to describe and sample from. As may be intuitive from results on polynomial interpolation, bandlimited approximation, and other function fitting problems, it is more concentrated towards the endpoints of $[0,T]$, so our sampling scheme selects more time points in these regions. The density is shown in Figure 3.

2.2 Algorithmic complexity

While Theorem 2 immediately yields an approach for selecting samples $t_{1},\ldots,t_{q}$, it is only useful if we can efficiently solve Problem 1 given the noisy measurements $y(t_{1})+n(t_{1}),\ldots,y(t_{q})+n(t_{q})$. We show that this is possible for a broad class of constraint measures. Specifically, we need only assume that we can efficiently compute the positive-definite kernel function111111When $y$ is real valued, it makes sense to consider symmetric $\mu$. In this case, $k_{\mu}$ is also real valued. However, in general it may be complex valued.:

[TABLE]

The above integral can be approximated via numerical quadrature, but for many of the aforementioned applications, it has a closed-form. For example, when $\mu$ is supported on just $k$ frequencies, it is a sum of these frequencies. When $\mu$ is uniform on $[-F,F]$, $k_{\mu}(t_{1},t_{2})=\operatorname{sinc}(2\pi F(t_{1}-t_{2}))$. For multiband signals with $s$ bands, $k_{\mu}(t_{1},t_{2})$ is a sum of $s$ modulated sinc functions. In fact, $k_{\mu}(t_{1},t_{2})$ has a closed-form for all $\mu$ illustrated in Figure 2. Further details are discussed in Appendix F. Assuming a subroutine for computing $k_{\mu}(t_{1},t_{2})$, our main algorithmic result is as follows:

Theorem 3**.**

*(Main result, algorithmic complexity)

There is an algorithm that solves Problem 1 with probability $99/100$ which uses $O\left(s_{\mu,\epsilon}\cdot\log^{2}(s_{\mu,\epsilon})\right)$ time domain samples (sampled according to the distribution given by Theorem 2) and runs in $\tilde{O}(s_{\mu,\epsilon}^{\omega}+s_{\mu,\epsilon}^{2}\cdot Z)$ time, assuming the ability to compute $k_{\mu}(t_{1},t_{2})$ for any $t_{1},t_{2}\in[0,T]$ in $Z$ time.121212For conciseness, we use $\tilde{O}(z)$ to denote $\tilde{O}(z\log^{c}z)$, where $c$ is some fixed constant (usually $\leq 2$). In formal theorem statements we give $c$ explicitly. $\omega<2.373$ is the current exponent of fast matrix multiplication [Wil12]. The algorithm returns a representation of $\tilde{y}(t)$ that can be evaluated in $\tilde{O}(s_{\mu,\epsilon}\cdot Z)$ time for any $t$.*

For bandlimited, Gaussian, or Cauchy-Lorentz priors $\mu$, $Z=O(1)$. For $s$ sparse signals or multiband signals with $s$ blocks, $Z=O(s)$.

We note that, while Theorem 3 holds when $\tilde{O}\left(s_{\mu,\epsilon}\right)$ samples are taken, $s_{\mu,\epsilon}$ may be not be known and thus it may be unclear how to set the sample size. In our full statement of the Theorem in Section 5.4 we make it clear that any upper bound on $s_{\mu,\epsilon}$ suffices to set the sample size. The sample complexity will depend on how tight this upper bound is. In Appendix E we give upper bounds on $s_{\mu,\epsilon}$ for a number of common $\mu$, which can be plugged into Theorem 3.

2.3 Our approach

Theorems 1, 2, and 3 are achieved through a simple and practical algorithmic framework. In Section 4, we show that Problem 1 can be modeled as a least squares regression problem with $\ell_{2}$ regularization. As long as we can compute $k_{\mu}(t_{1},t_{2})$, we can solve this problem using kernel ridge regression, a popular function fitting technique in nonparametric statistics [STC04].

Naively, the kernel regression problem is infinite dimensional: it needs to be solved over the continuous time domain $[0,T]$ to solve our signal reconstruction problem. This is where sampling comes in. We need to discretize the problem and establish that our solution over a fixed set of time samples nearly matches the solution over the continuous interval. To bound the error of discretization, we turn to a tool from randomized numerical linear algebra: statistical leverage score sampling [SS11, DM16]. We show how to randomly discretize Problem 1 by sampling time points with probability proportional to an appropriately defined non-uniform leverage score distribution on $[0,T]$. The required number of samples is $O(s_{\mu,\epsilon}\log s_{\mu,\epsilon})$, which proves Theorem 1.

Unfortunately, the leverage score distribution does not have a closed-form, varies depending on $\epsilon$, $T$, and $\mu$, and likely cannot be sampled from exactly. To prove Theorem 2, we show that for any $\mu$, for large enough $q$, the closed form distribution $p_{q}$ upper bounds the leverage score distribution. This upper bound closely approximates the true leverage score distribution and, therefore, can be used in its place during sampling, losing only a $\log s_{\mu,\epsilon}$ factor in the sample complexity.

The leverage score distribution roughly measures, for each time point $t$, how large $|y(t)|^{2}$ can be compared to $\|y\|_{T}^{2}$ when $y$’s Fourier transform is constrained by $\mu$ (i.e., when $\|x\|_{\mu}^{2}$ as defined in Problem 1 is bounded). To upper bound this measure we turn to another powerful result from the randomized numerical linear algebra literature: every matrix contains a small subset of columns that span a near-optimal low-rank approximation to that matrix [Sar06, BMD09, DR10]. In other words, every matrix admits a near-optimal low-rank approximation with sparse column support. By extending this result to continuous linear operators, we prove that the smoothness of a signal whose Fourier transform has $\|x\|_{\mu}^{2}$ bounded can be bounded by the smoothness of an $O(s_{\mu,\epsilon})$ sparse Fourier function. This lets us apply recent results of [CKPS16, CP18] that bound $|y(t)|^{2}$ in terms of $\|y\|_{T}^{2}$ for any sparse Fourier function $y$ . Intuitively, our result shows that the simplicity of sparse Fourier functions governs the simplicity of any class of Fourier constrained functions.

The above argument yields Theorem 2. Since we can sample from $p_{q}$ in $O(1)$ time, we can efficiently sample the time domain to $O(s_{\mu,\epsilon}\cdot\log^{2}s_{\mu,\epsilon})$ points and then solve Problem 1 by applying kernel ridge regression to these points, which takes $\tilde{O}(s_{\mu,\epsilon}^{\omega}+s_{\mu,\epsilon}^{2}\cdot Z)$ time, assuming the ability to compute $k_{\mu}(\cdot,\cdot)$ in $Z$ time. This yields the algorithmic result of Theorem 3.

2.4 Roadmap

The rest of this paper is devoted to proving Theorems 1, 2, and 3, and is structured as follows:

Section 3

We lay out basic notation that is used throughout the paper.

Section 4

We reduce Problem 1 to a kernel ridge regression problem and explain how to randomly discretize and solve this problem via leverage score sampling, proving Theorem 1.

Section 5

We give an upper bound on the leverage score distribution for general priors, proving Theorems 2 and 3.

Section 6

We prove that, under a mild assumption, the statistical dimension tightly characterizes the sample complexity of solving Problem 1, and thus that our results are nearly optimal.

Section 7

We conclude with a discussion of open questions.

We defer an in depth overview of related work to Appendix A. In Appendix B we give operator theory preliminaries. In Appendix C we prove our extensions of a number of randomized linear algebra primitives to continuous operators. In Appendix D, we bound the statistical dimension for the important case of bandlimited functions. We use this result in Appendix E to prove statistical dimension bounds for multiband, Gaussian, and Cauchy-Lorentz priors (shown in Table 1). In Appendix F, we show how to compute the kernel function $k_{\mu}$ for these common priors. In Appendix G, we discuss a Bayesian approach to signal reconstruction under a Fourier transform prior.

3 Notation

Let $\mu$ be a probability measure on $(\mathbb{R},\mathcal{B})$, where $\mathcal{B}$ is the Borel $\sigma$-algebra on $\mathbb{R}$. Let $L_{2}(\mu)$ denote the space of complex-valued square integrable functions with respect to $\mu$. For $a,b\in L_{2}(\mu)$, let $\langle a,b\rangle_{\mu}$ denote $\int_{\xi\in\mathbb{R}}a(\xi)^{*}b(\xi)\,d\mu(\xi)$ where for any $x\in\mathbb{C}$, $x^{*}$ is its complex conjugate. Let $\|a\|_{\mu}^{2}$ denote $\langle a,a\rangle_{\mu}$. Let $\mathcal{I}_{\mu}$ denote the identity operator on $L_{2}(\mu)$. Note that for any $\mu$, $L_{2}(\mu)$ is a separable Hilbert space and thus has a countably infinite orthonormal basis [HN01].

We overload notation and use $L_{2}(T)$ to denote the space of complex-valued square integrable functions with respect to the uniform probability measure on $[0,T]$. It will be clear from context that $T$ is not a measure. For $a,b\in L_{2}(T)$, let $\langle a,b\rangle_{T}$ denote $\frac{1}{T}\int_{0}^{T}a(t)^{*}b(t)\,dt$ and let $\|a\|_{T}^{2}$ denote $\langle a,a\rangle_{T}$. Let $\mathcal{I}_{T}$ denote the identity operator on $L_{2}(T)$.

Define the Fourier transform operator $\mathcal{F}_{\mu}:L_{2}(T)\rightarrow L_{2}(\mu)$ as:

[TABLE]

The adjoint of $\mathcal{F}_{\mu}$ is the unique operator $\mathcal{F}_{\mu}^{*}:L_{2}(\mu)\rightarrow L_{2}(T)$ such that for all $f\in L_{2}(T),g\in L_{2}(\mu)$ we have $\langle g,\mathcal{F}_{\mu}\,f\rangle_{\mu}=\langle\mathcal{F}_{\mu}^{*}\,g,f\rangle_{T}$. It is not hard to see that $\mathcal{F}_{\mu}^{*}$ is the inverse Fourier transform operator with respect to $\mu$ as defined in Section 2, equation (1):

[TABLE]

Note that the kernel operator $\mathcal{K}_{\mu}:L_{2}(T)\rightarrow L_{2}(T)$ originally defined in (3) is equal to

[TABLE]

$\mathcal{K}_{\mu}$ is self-adjoint, positive semidefinite and trace-class and an integral operator with kernel $k_{\mu}$:

[TABLE]

where $k_{\mu}$ is as defined in (6). The trace of $\mathcal{K}_{\mu}$ is equal to $1$.131313Since the kernel is a Fourier transform of a probability measure, it is Hermitian positive definite (Bochner’s Theorem). Then we can conclude that $\mathcal{K}_{\mu}$ is trace-class from Mercer’s theorem, and calculate $\operatorname{tr}(\mathcal{K}_{\mu})=\frac{1}{T}\int^{T}_{0}k_{\mu}(t,t)dt=1$. We will also make use of the Gram operator: $\cal{G}_{\mu}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\cal{F}_{\mu}\cal{F}^{*}_{\mu}$. ${\cal G}_{\mu}$ is also self-adjoint, positive semidefinite, and trace-class.

Remark: It may be useful for the reader to informally regard $\mathcal{F}_{\mu}$ as an infinite matrix with rows indexed by $\xi\in\mathbb{R}$ and columns indexed by $t\in[0,T]$. Following the definition of $\mathcal{F}_{\mu}$ above, and assuming that $\mu$ has a density $p$, this infinite matrix has entries given by:

[TABLE]

The results we apply on leverage score sampling can all be seen as extending results for finite matrices from the randomized numerical linear algebra literature to this infinite matrix.

4 Function fitting with least squares regression

Least squares regression provides a natural approach to solving the interpolation task of Problem 1. In particular, consider the following regularized minimization problem over functions $g\in L_{2}(\mu)$141414The fact that the minimum is attainable is a simple consequence of the extreme value theorem, since the search space can be restricted to $\|g\|_{\mu}^{2}\leq\|(y+n)\|^{2}_{T}/\epsilon$.:

[TABLE]

The first term encourages us to find a function $g$ whose inverse Fourier transform is close to our measured signal $y+n$. The second term encourages us to find a low energy solution – ultimately, we solve (10) based on only a small number of samples $y(t_{1}),\ldots,y(t_{k})$, and smoother, lower energy solutions will better generalize to the entire interval $[0,T]$. We remark that it is well known that least squares approximations benefit from regularization even in the noiseless case [CDL13].

We first state a straightforward fact: if we minimize (10), even to a coarse approximation, then we are able to solve Problem 1.

Claim 4**.**

Let $y=\mathcal{F}_{\mu}^{*}x$, $n\in L_{2}(T)$ be an arbitrary noise function, and for any $C\geq 1$, let $\tilde{g}\in L_{2}(\mu)$ be a function satisfying:

[TABLE]

Then

[TABLE]

Proof.

Since $y=\mathcal{F}_{\mu}^{*}x$, $\min_{g\in L_{2}(\mu)}\left[\|\mathcal{F}_{\mu}^{*}g-(y+n)\|_{T}^{2}+\epsilon\|g\|_{\mu}^{2}\right]\leq\|n\|_{T}^{2}+\epsilon\|{x}\|_{\mu}^{2}$. Thus, $\|\mathcal{F}_{\mu}^{*}\tilde{g}-(y+n)\|_{T}^{2}\leq C\epsilon\|x\|_{\mu}^{2}+C\|n\|_{T}^{2}$. The claim then follows via triangle inequality:

[TABLE]

∎

Claim 4 shows that approximately solving the regression problem in (10), with regularization parameter $\epsilon$ gives a solution to Problem 1 with parameter $2C\epsilon$ (decreasing the regularization parameter to $\frac{\epsilon}{2C}$ will let us solve with parameter $\epsilon$). But how can we solve the regression problem efficiently? Not only does the problem involve a possibly infinite dimensional parameter vector $g$, but the objective function also involves the continuous time interval $[0,T]$.

4.1 Random discretization via leverage function sampling

The first step is to deal with the latter challenge, i.e., that of a continuous time domain. We show that it is possible to randomly discretize the time domain of (10), thereby reducing our problem to a regression problem on a finite set of times $t_{1},\ldots,t_{q}$. In particular, we can sample time points with probability proportional to the so-called ridge leverage function, a specific non-uniform distribution that has been applied widely in randomized algorithms for regression and other linear algebra problems on discrete matrices [AM15, CLV16, CMM17, MM17, MW17].

While we cannot compute the leverage function explicitly for our problem, an issue highlighted in [Bac17], our main result (Theorem 2) uses a simple, but very accurate, closed form approximation in its place. We start with the definition of the ridge leverage function:

Definition 3** (Ridge leverage function).**

For a probability measure $\mu$ on $\mathbb{R}$, time length $T>0$, and $\epsilon\geq 0$, we define the $\epsilon$-ridge leverage function for $t\in[0,T]$ as151515Formally $L_{2}(T)$ is a space of equivalence classes of functions that differ at a set of points with measure [math]. For notational simplicity, here and throughout we use $\mathcal{F}_{\mu}^{*}\alpha$ to denote the specific representative of the equivalence class $\mathcal{F}_{\mu}^{*}\alpha\in L_{2}(T)$ given by (8). In this way, we can consider the pointwise value $[\mathcal{F}_{\mu}^{*}\alpha](t)$, which we could alternatively express as $\langle\varphi_{t},\alpha\rangle_{\mu}$, for $\varphi_{t}(\xi)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}e^{-2\pi it\xi}$.:

[TABLE]

Intuitively, the ridge leverage function at time $t$ is an upper bound of how much a function can “blow up” at $t$ when its Fourier transform is constrained by $\mu$. The denominator term $\|\mathcal{F}_{\mu}^{*}\alpha\|_{T}^{2}$ is the average squared magnitude of the function $F_{\mu}^{*}\alpha$, while the numerator term, $|[\mathcal{F}_{\mu}^{*}\alpha](t)|^{2}$, is the squared magnitude at $t$. The regularization term $\epsilon\|\alpha\|_{\mu}^{2}$ reflects the fact that, to solve (10), we only need to bound the smoothness for functions with bounded Fourier energy under $\mu$. As observed in [PBV18], the ridge leverage function can be viewed as a type of Christoffel function, studied in the literature on orthogonal polynomials and approximation theory [PBV18, Nev86, Tot00, BE12].

The larger the leverage “score” $\tau_{\mu,\epsilon}(t)$, the higher the probability we will sample time $t$, to ensure that our sample points well reflect any possibly significant components or ‘spikes’ of the function $y$. Ultimately, the integral of the ridge leverage function $\int_{0}^{T}\tau_{\mu,\epsilon}(t)dt$ determines how many samples we require to solve (10) to a given accuracy. Theorem 5 below states the already known fact that the ridge leverage function integrates to the statistical dimension [AKM*+*17], which will ultimately allow us to achieve the $\tilde{O}(s_{\mu,\epsilon})$ sample complexity bound of Theorems 1 and 2. Theorem 5 also gives two alternative characterizations of the leverage function that will prove useful. The theorem is proven in Appendix C, using techniques for finite matrices, adapted to the operator setting.

Theorem 5** (Leverage function properties).**

Let $\tau_{\mu,\epsilon}(t)$ be the ridge leverage function (Definition 3) and define $\varphi_{t}\in L_{2}(\mu)$ by $\varphi_{t}(\xi)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}e^{-2\pi it\xi}$. We have:

• •

The ridge leverage function integrates to the statistical dimension:

[TABLE]

• •

Inner Product characterization:

[TABLE]

• •

Minimization Characterization:

[TABLE]

In Theorem 6, we give our formal statement that the ridge leverage function can be used to randomly sample time domain points to discretize the regression problem in (10) and solve it approximately. While complex in appearance, readers familiar with randomized linear algebra will recognize Theorem 6 as closely analogous to standard approximate regression results for leverage score sampling from finite matrices [CW13]. As discussed, since we are typically unable to sample according to the true ridge leverage function, we give a general result, showing that sampling with any upper bound function with a finite integral suffices.

Theorem 6** (Approximate regression via leverage function sampling).**

Assume that $\epsilon\leq\|\mathcal{K}_{\mu}\|_{\mathrm{op}}$.161616If $\epsilon>\|\mathcal{K}_{\mu}\|_{\mathrm{op}}$ then (10) is solved to a constant approximation factor by the trivial solution $g=0$. Consider a measurable function $\tilde{\tau}_{\mu,\epsilon}(t)$ with $\tilde{\tau}_{\mu,\epsilon}(t)\geq\tau_{\mu,\epsilon}(t)$ for all $t$ and let $\tilde{s}_{\mu,\epsilon}=\int_{0}^{T}\tilde{\tau}_{\mu,\epsilon}(t)dt$. Let $s=c\cdot\tilde{s}_{\mu,\epsilon}\cdot\left(\log\tilde{s}_{\mu,\epsilon}+1/\delta\right)$ for sufficiently large fixed constant $c$ and let $t_{1},\ldots,t_{s}$ be time points selected by drawing each randomly from $[0,T]$ with probability proportional to $\tilde{\tau}_{\mu,\epsilon}(t)$. For $j\in 1,\ldots,s$, let $w_{j}=\sqrt{\frac{1}{sT}\cdot\frac{\tilde{s}_{\mu,\epsilon}}{\tilde{\tau}_{\mu,\epsilon}(t_{j})}}$. Let $\mathbf{F}:\mathbb{C}^{s}\rightarrow L_{2}(\mu)$ be the operator defined by:

[TABLE]

and $\mathbf{y},\mathbf{n}\in\mathbb{R}^{s}$ be the vectors with $\mathbf{y}(j)=w_{j}\cdot y(t_{j})$ and $\mathbf{n}(j)=w_{j}\cdot n(t_{j})$. Let:

[TABLE]

With probability $\geq 1-\delta$:

[TABLE]

A generalized version of this result is proven in Appendix C, which holds even when $\tilde{g}$ is only an approximate minimizer of (15).

Theorem 6 shows that $\tilde{g}$ obtained from solving the discretized regression problem provides an approximate solution to (10) and by Claim 4, $\tilde{y}=\mathcal{F}_{\mu}^{*}\tilde{g}$ solves Problem 1 with parameter $\Theta(\epsilon)$. If we have $\tilde{\tau}_{\mu,\epsilon}(t)=\tau_{\mu,\epsilon}(t)$, Theorem 6 combined with Claim 4 shows that Problem 1 with parameter $\Theta(\epsilon)$ can be solved with sample complexity $O\left(s_{\mu,\epsilon}\cdot\log s_{\mu,\epsilon}\right)$, since by (12), $\int_{0}^{T}\tau_{\mu,\epsilon}(t)dt=s_{\mu,\epsilon}$. Note that, by simply decreasing the regularization parameter in (10) by a constant factor, we can solve Problem 1 with parameter $\epsilon.$ The asymptotic complexity is identical since, by (14), for any $c\leq 1$ and any $t\in[0,T]$, $\tau_{\mu,c\epsilon}(t)\leq\frac{1}{c}\tau_{\mu,\epsilon}(t)$ and so:

[TABLE]

This proves the sample complexity result of Theorem 1. However, since it is not clear that sampling according to $\tau_{\mu,\epsilon}(t)$ can be done efficiently (or at all), it does not yet give an algorithm yielding this complexity.171717We conjecture that the existential sample complexity can in fact be upper bounded by $O(s_{\mu,\epsilon})$ by adapting deterministic sampling methods for finite matrices to the operator setting [CNW16], like we do in Lemma 46. This issue will be addressed in Section 5, where we prove Theorem 2.

We prove Theorem 6 in Appendix C. We show that leverage function sampling satisfies, with good probability, an affine embedding guarantee: that $\|\mathbf{F}^{*}{g}-(\mathbf{y}+\mathbf{n})\|_{2}^{2}+\epsilon\|{g}\|_{\mu}^{2}$ closely approximates $\|\mathcal{F}_{\mu}^{*}{g}-(y+n)\|_{T}^{2}+\epsilon\|{g}\|_{\mu}^{2}$ for all $g\in L_{2}(\mu)$. Thus, a (near) optimal solution to the discretized problem, $\min_{g\in L_{2}(\mu)}\left[\|\mathbf{F}^{*}g-(\mathbf{y}+\mathbf{n})\|_{2}^{2}+\epsilon\|g\|_{\mu}^{2}\right]$, gives a near optimal solution to the original problem, $\min_{g\in L_{2}(\mu)}\left[\|\mathcal{F}_{\mu}^{*}g-(y+n)\|_{T}^{2}+\epsilon\|g\|_{\mu}^{2}\right]$. Our proof of the affine embedding property is analogous to existing proofs for finite dimensional matrices [CW13, ACW17].

4.2 Efficient solution of the discretized problem

Given an upper bound on the ridge leverage function $\tilde{\tau}_{\mu,\epsilon}(t)\geq\tau_{\mu,\epsilon}(t)$, we can apply Theorem 6 to approximately solve the ridge regression problem of (10) and therefore Problem 1 by Claim 4. In Section 5 we show how to obtain such an upper bound for any $\mu$ using a universal distribution.

First, however, we demonstrate how to apply Theorem 6 algorithmically. Specifically, we show how to solve the randomly discretized problem of (15) efficiently. Combined with Theorem 6 and our bound on $\tau_{\mu,\epsilon}(t)$ given in Section 5, this yields a randomized algorithm (Algorithm 1) for Problem 1. The formal analysis of Algorithm 1 is given in Theorem 7.

Theorem 7** (Efficient signal reconstruction given leverage function upper bounds).**

Assume that $\epsilon\leq\|\mathcal{K}_{\mu}\|_{\mathrm{op}}$.181818As discussed for Theorem 6, if $\epsilon>\|\mathcal{K}_{\mu}\|_{\mathrm{op}}$, Problem 1 is trivially solved by $\tilde{y}=0$. Algorithm 1 returns $t_{1},\ldots,t_{s}\in[0,T]$ and $\mathbf{z}\in\mathbb{C}^{s}$ such that $\tilde{y}(t)=\sum_{i=1}^{s}\mathbf{z}(i)\cdot k_{\mu}(t_{i},t)$ (as computed in Algorithm 2) satisfies with probability $\geq 1-\delta$:

[TABLE]

Suppose we can sample $t\in[0,T]$ with probability proportional to $\tilde{\tau}_{\mu,\epsilon}(t)$ in time $W$ and compute the kernel function $k_{\mu}(t_{1},t_{2})=\int_{\xi\in\mathbb{R}}e^{-2\pi i(t_{1}-t_{2})}d\mu(\xi)$ in time $Z$. Algorithm 1 queries $y+n$ at $s$ points and runs in $O\left(s\cdot W+s^{2}\cdot Z+s^{\omega}\right)$ time191919Here $\omega<2.373$ is the exponent of fast matrix multiplication. $s^{\omega}$ is the theoretically fastest runtime required to invert a dense $s\times s$ matrix. We note that the $s^{\omega}$ term may be thought of as $s^{3}$ in practice, and potentially could be accelerated using a variety of techniques for fast (regularized) linear system solvers. where $s=O\left(\tilde{s}_{\mu,\epsilon}\cdot\left(\log\tilde{s}_{\mu,\epsilon}+1/\delta\right)\right)$. Algorithm 2 evaluates $\tilde{y}(t)$ in $O(s\cdot Z)$ time for any $t$.

Proof.

In Step 2 of Algorithm 1, $t_{1},\ldots,t_{s}$ are sampled according to $\tilde{\tau}_{\mu,\epsilon}(t)$, which upper bounds $\tau_{\mu,\epsilon}(t)$. We can thus apply Theorem 6. If the constant $c$ in Step 1 is set large enough, with probability $\geq 1-\delta$, letting $\mathbf{F},\mathbf{y},$ and $\mathbf{n}$ be as defined in that theorem, (16) holds for

[TABLE]

Therefore, letting $\tilde{y}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\mathcal{F}_{\mu}^{*}\tilde{g}$ and applying Claim 4, with probability $\geq 1-\delta$,

[TABLE]

Further, the minimizer $\tilde{g}$ is indeed unique and can be written as (see Lemma 38 in Appendix C):

[TABLE]

where $\mathbf{K}=\mathbf{F}^{*}\mathbf{F}$ is as defined in Step 3 of Algorithm 1 and $\mathbf{\bar{y}}=\mathbf{y}+\mathbf{n}$ is formed in Step 4. If we let $\mathbf{\bar{z}}=(\mathbf{K}+\epsilon\mathbf{I})^{-1}\mathbf{\bar{y}}$ and let $\mathbf{z}$ have $\mathbf{z}(i)=\mathbf{\bar{z}}(i)\cdot w_{i}$ as in Steps 5 and 6, we can see that:

[TABLE]

giving the expression returned in Algorithm 2. Combined with (18), this completes the accuracy bound of the theorem. The runtime and sample complexity bounds follow from observing that:

• •

$s\cdot W$ time is required to sample $t_{1},\ldots,t_{s}$ in Step 2.

• •

$s^{2}\cdot Z$ time is required to form $\mathbf{K}$ in Step 3.

• •

$s$ queries to $y+n$ are required to form $\mathbf{\bar{y}}$ in Step 4.

• •

$O(s^{\omega})$ time is required to compute $\mathbf{\bar{z}}:=(\mathbf{K}+\epsilon\mathbf{I})^{-1}\mathbf{\bar{y}}$ in Step 5. This runtime could potentially be improved with a variety of fast system solvers. We take $s^{\omega}$ as a simple upper bound.

• •

$O(s\cdot Z)$ time is required to compute $k(t_{1},t),\ldots,k(t_{s},t)$ to evaluate $\tilde{y}(t)$ in Algorithm 2.

This completes the proof of Theorem 7. ∎

**Remark: ** As discussed, in Section 5 we will give a ridge leverage function upper bound that can be sampled from in $W=O(1)$ time and closely bounds the true leverage function for any $\mu$, giving $\tilde{s}_{\mu,\epsilon}=O(s_{\mu,\epsilon}\cdot\log s_{\mu,\epsilon})$. Using this upper bound to sample time domain points, our sample complexity $s$ is thus within a $O(\log s_{\mu,\epsilon})$ factor of the best possible using Theorem 6, which we would achieve if sampling using the true ridge leverage function.

In Appendix D we prove a tighter leverage function bound than the one in Section 5 for bandlimited signals, removing the logarithmic factor in this case. It is not hard to see that for general $\mu$ we can also achieve optimal sample complexity by further subsampling $t_{1},\ldots,t_{s}$ using the ridge leverage scores of $\mathbf{K}^{1/2}$. These scores can be computed in $\tilde{O}(s\cdot s_{\mu,\epsilon}^{2})$ time using known techniques for finite kernel matrices [MM17]. Subsampling $O\left(\frac{s_{\mu,\epsilon}\log s_{\mu,\epsilon}}{\delta^{2}}\right)$ time domain points according to these scores lets us approximately solve the discretized problem of (15) to error $(1+\delta)$.

Applying the more general version of Theorem 6 stated in Appendix C, this yields an approximate solution to (10) and thus to Problem 1. For constant $\delta$, we need just $O(s_{\mu,\epsilon}\cdot\log s_{\mu,\epsilon})$ time samples to to solve the subsampled regression problem, matching the best possible sample complexity of Theorem 6. By the lower bound given in Section 6, Theorem 24, this complexity is within a $O(\log s_{\mu,\epsilon})$ factor of optimal in nearly all settings. We conjecture that one can in fact achieve within an $O(1)$ factor of the optimal sample complexity by applying deterministic selection methods to $\mathbf{F}$ [CNW16], similar to the techniques used to prove Lemma 46.

5 A near-optimal spectrum blind sampling distribution

In the previous section, we showed how to solve Problem 1 given the ability to sample time points according to the ridge leverage function $\tau_{\mu,\epsilon}$. In general, this function depends strongly on $T$, $\mu$, and $\epsilon$, and it is not clear if it can be computed or sampled from directly.

Nevertheless, in this section we show that it is possible to efficiently obtain samples from a function that very closely approximates the true leverage function for any constraint measure $\mu$. In particular we describe a set of closed form functions $\tilde{\tau}_{\alpha}(t)$, each parameterized by $\alpha>0$. $\tilde{\tau}_{\alpha}$ upper bounds the leverage function $\tau_{\mu,\epsilon}$ for any $\mu$ and $\epsilon$, as long as the statistical dimension $s_{\mu,\epsilon}\leq O(\alpha)$. Our upper bound satisfies

[TABLE]

which means it can be used in place of the true ridge leverage function to give near optimal sample complexity via Theorem 6 and 7. This result is proven formally in Theorem 17, which as a consequence immediately yields our main technical result, Theorem 2. The majority of this section is devoted towards building tools necessary for proving Theorem 17.

5.1 Uniform leverage bound via Fourier sparsification

We seek a simple closed form function that upper bounds the leverage function $\tau_{\mu,\epsilon}$. Ultimately, we want this upper bound to be very tight, but a natural first question is whether it should exists at all. Is it possible to prove any finite upper bound on $\tau_{\mu,\epsilon}$ without using specific knowledge of $\mu$?

We answer this first question by showing that $\tau_{\mu,\epsilon}$ can be upper bounded by a constant function. Specifically, we show that for $t\in[0,T]$, $\tau_{\mu,\epsilon}(t)\leq C$ for $C=\operatorname{poly}(s_{\mu,\epsilon})$. This upper bound depends on the statistical dimension, but importantly, it does not depend on $\mu$. Formally we show:

Theorem 8** (Uniform leverage function bound).**

For all $t\in[0,T]$ and $\epsilon\leq 1$202020If $\epsilon>1=\operatorname{tr}(\mathcal{K}_{\mu})$, Problem 1 is trivially solved by returning $\tilde{y}=0$.

[TABLE]

While Theorem 8 appears to give a relatively weak bound, proving this statement is a key technical challenge. Ultimately, it is used in Section 5.3 as one of two main ingredients in proving the much tighter leverage function bound that yields Theorem 17 and Theorem 2.

Towards a proof of Theorem 8, we consider the operator $\mathcal{F}_{\mu}$ defined in Section 3. Since $\mathcal{F}_{\mu}$ has statistical dimension $s_{\mu,\epsilon}$, $\mathcal{K}_{\mu}=\mathcal{F}_{\mu}^{*}\mathcal{F}_{\mu}$ can have at most $2s_{\mu,\epsilon}$ eigenvalues $\geq\epsilon$:

[TABLE]

Thus, if we project $\mathcal{F}_{\mu}$ onto the span of $\mathcal{K}_{\mu}$’s top $2s_{\mu,\epsilon}$ eigenfunctions (when $\mu$ is uniform on an interval these are the prolate spherical wave functions of Slepian and Pollak [SP61]) we will approximate $\mathcal{K}_{\mu}$ up to its small eigenvalues. The total mass of these eigenvalues is bounded by:

[TABLE]

Alternatively, instead of projecting onto the span of the eigenfunctions, we can approximate $\mathcal{K}_{\mu}$ nearly optimally by projecting $\mathcal{F}_{\mu}$ onto the span of a subset of $O(s_{\mu,\epsilon})$ of its “rows” – i.e. frequencies in the support of $\mu$. For finite linear operators, is well known that such a subset exists: the problem of finding these subsets has been studied extensively in the literature on randomized low-rank matrix approximation under the name column subset selection [Sar06, BMD09, DR10]. In Appendix C we show that an analogous result extends to the continuous operator $\mathcal{F}_{\mu}$:

Theorem 9** (Frequency subset selection).**

For some $s\leq\lceil 36\cdot s_{\mu,\epsilon}\rceil$ there exists a set of distinct frequencies $\xi_{1},\ldots,\xi_{s}\in\mathbb{C}$ such that, if $\mathbf{C}_{s}:L_{2}(T)\rightarrow\mathbb{C}^{s}$ and $\mathbf{Z}:L_{2}(\mu)\rightarrow\mathbb{C}^{s}$ are defined by:

[TABLE]

then

[TABLE]

Note that, if $\varphi_{t}\in L_{2}(\mu)$ is defined $\varphi_{t}(\xi)=e^{-2\pi it\xi}$ and $\boldsymbol{\phi}_{t}\in\mathbb{C}^{s}$ is defined $\boldsymbol{\phi}_{t}(j)=\varphi_{t}(\xi_{j})$, we have:

[TABLE]

Leverage function bound proof sketch. With Theorem 9 in place, we explain how to use this result to prove Theorem 8, i.e., to establish a universal bound on the leverage function of $\mathcal{F}_{\mu}$. For the sake of exposition, we use the term “row” of an operator $\mathcal{A}:L_{2}(\mu)\rightarrow L_{2}(T)$ to refer to the corresponding operator restricted to some time $t$. We use the term “column” of an operator as the adjoint of a row of $\mathcal{A}^{*}:L_{2}(T)\rightarrow L_{2}(\mu)$, i.e., the adjoint operator restricted to some frequency $\xi$.

By Theorem 9, $\mathbf{C}_{s}^{*}\mathbf{Z}:L_{2}(\mu)\rightarrow L_{2}(T)$ (the projection of $\mathcal{F}_{\mu}^{*}$ onto the range of $\mathbf{C}_{s}$) closely approximates the operator $\mathcal{F}_{\mu}^{*}$ yet has columns spanned by just $O(s_{\mu,\epsilon})$ frequencies: $\xi_{1},\ldots,\xi_{s}$. Thus, for any $\alpha\in L_{2}(\mu)$, $\mathbf{C}_{s}^{*}\mathbf{Z}\alpha\in L_{2}(T)$ is just an $O(s_{\mu,\epsilon})$ sparse Fourier function. Using the maximization characterization of Definition 3, we can thus bound the time domain ridge leverage function of $\mathbf{C}_{s}^{*}\mathbf{Z}$ by appealing to known smoothness bounds for Fourier sparse functions [CP18], even for $\epsilon=0$. When $\epsilon=0$, the ridge leverage function is known as the standard leverage function in the randomized numerical linear algebra literature, and we will refer to them as such.

We can use a similar argument to bound the row norms of the residual operator $[\mathcal{F}_{\mu}^{*}-\mathbf{C}_{s}^{*}\mathbf{Z}]$. The columns of this residual operator are each spanned by $O(s_{\mu,\epsilon})$ frequencies, and so are again sparse Fourier functions whose smoothness we can bound. This smoothness ensures that no row can have norm significantly higher than average.

Finally, we note that the time domain ridge leverage function of $\mathcal{F}_{\mu}$ is approximated to within a constant factor by the sum of the standard row leverage function of $\mathbf{C}_{s}^{*}\mathbf{Z}$ along with row norms of $\mathcal{F}_{\mu}-\mathbf{C}_{s}^{*}\mathbf{Z}$. This gives us a bound on $\mathcal{F}_{\mu}$’s ridge leverage function. We prove this formally below:

Theorem 10** (Ridge leverage function approximation).**

Let $\mathbf{C}_{s}$ and $\mathbf{Z}$ be the operators guaranteed to exist by Theorem 9. Let $\ell(t)$ be the standard leverage function of $t$ in $\mathbf{C}_{s}^{*}\mathbf{Z}$:222222Analogously to how $[\mathcal{F}_{\mu}^{*}\alpha](t)$ is used in Definition 3, while $L_{2}(T)$ is formally a space of equivalence classes of functions, here we use $\mathbf{C}^{*}_{s}\mathbf{Z}\alpha$ to denote the specific representative of the equivalence class $\mathbf{C}^{*}_{s}\mathbf{Z}\alpha\in L_{2}(T)$ given by $[\mathbf{C}^{*}_{s}\mathbf{Z}\alpha](t)=\sum_{j=1}^{s}[\mathbf{Z}\alpha](j)\cdot e^{2\pi i\xi_{j}t}=\langle\boldsymbol{\phi}_{t},\mathbf{Z}\alpha\rangle_{\mathbb{C}^{s}}$. In this way, we can consider the pointwise value $[\mathbf{C}^{*}_{s}\mathbf{Z}\alpha](t)$.

[TABLE]

Let $r(t)$ be the residual:

[TABLE]

where $\varphi_{t}$ and $\boldsymbol{\phi}_{t}$ are as defined in Theorem 9. Then for all $t$:

[TABLE]

Proof.

For any $\alpha\in L_{2}(\mu)$ we can write $[\mathcal{F}_{\mu}^{*}\alpha](t)=\langle\varphi_{t},\alpha\rangle_{\mu}$ and $[\mathbf{C}_{s}^{*}\mathbf{Z}\alpha](t)=\langle\boldsymbol{\phi}_{t},\mathbf{Z}\alpha\rangle_{\mathbb{C}^{s}}=\langle\mathbf{Z}^{*}\boldsymbol{\phi}_{t},\alpha\rangle_{\mu}$. By the maximization characterization of the ridge leverage function in Definition 3,

[TABLE]

where the second to last line follows from observing that due to Cauchy-Schwarz,

[TABLE]

and that, letting $\mathcal{P}_{s}=\mathbf{C}_{s}^{*}(\mathbf{C}_{s}\mathbf{C}_{s}^{*})^{-1}\mathbf{C}_{s}$:

[TABLE]

In the above, the inequality is due to the fact that $\mathcal{P}_{s}$ is an orthogonal projection, so $\mathcal{P}_{s}\preceq{\cal I}_{\mu}$. This completes the proof. ∎

With Theorem 10 in place, we now bound $\bar{\tau}_{\mu,\epsilon}(t)=2\left(\ell(t)+\frac{r(t)}{\epsilon}\right)$, which yields a uniform bound on the true ridge leverage scores.

Lemma 11**.**

Let $\ell(t),r(t)$ be as defined in Theorem 10 and $\bar{\tau}_{\mu,\epsilon}(t)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}2\cdot\left(\ell(t)+\frac{r(t)}{\epsilon}\right)$. For all $t\in[0,T]$:

[TABLE]

Combining Lemma 11 with Theorem 10 yields Theorem 8. We just simplify the constants by noting that for $\epsilon\leq 1$, $s_{\mu,\epsilon}\geq\frac{\operatorname{tr}(\mathcal{K}_{\mu})}{2}=\frac{1}{2}$ and so $36s_{\mu,\epsilon}+2\leq 40s_{\mu,\epsilon}$.

Proof of Lemma 11.

We separately bound the leverage score $\ell(t)$ and residual $r(t)$ components of $\bar{\tau}_{\mu,\epsilon}(t)$ using a similar argument based on the smoothness of sparse Fourier functions for both. Specifically, for both bounds we employ the following smoothness bound of Chen et al.:

Lemma 12** (Follows from Lemma 5.1 of [CKPS16]).**

For any $f(t)=\sum_{j=1}^{k}v_{j}e^{2\pi i\xi_{j}t}$,

[TABLE]

Proof.

This follows from Lemma 5.1 of [CKPS16], which gives the bound without an explicit constant. It is not hard to check that their proof gives the constant of $1540$ stated above. ∎

Bounding the leverage scores $\ell(t)$ of $\mathbf{C}_{s}^{*}\mathbf{Z}$.

For every $\alpha\in L_{2}(\mu)$, $\mathbf{C}_{s}^{*}\mathbf{Z}\alpha$ is an $s=O(s_{\mu,\epsilon})$ sparse Fourier function. Specifically, we have:

[TABLE]

for frequencies $\xi_{1},\ldots,\xi_{s}\in\mathbb{C}$ given by Theorem 9. We can thus directly apply Lemma 12 giving for any $t\in[0,T]$:

[TABLE]

Bounding the residuals $r(t)$.

We first give some intuition. To bound the squared row norms of the residual $\mathcal{F}_{\mu}^{*}-\mathbf{C}_{s}^{*}\mathbf{Z}$ we show that each “column” of this residual is an $s+1=O(s_{\mu,\epsilon})$ sparse Fourier function. Thus, applying Lemma 12, no entry’s squared value can significantly exceed the average squared value in the column. This lets us show that no squared row norm $r(t)$ can significantly exceed the average squared row norm, which is bounded by Theorem 9.

Concretely, define $\vartheta_{\xi}\in L_{2}(T)$ by $\vartheta_{\xi}(t)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}e^{2\pi it\xi}$, and notice that given $g\in L_{2}(T)$ the function $\xi\mapsto\langle\vartheta_{\xi},g\rangle_{T}$ is equal to $\mathcal{F}_{\mu}g$ in the $L_{2}(T)$ sense (i.e., is a member of the equivalence class $\mathcal{F}_{\mu}g$). For $\xi\in\mathbb{R}$, let $\mathbf{z}_{\xi}\in\mathbb{C}^{s}$ be given by $\mathbf{z}_{\xi}(j)=\langle\vartheta_{\xi},\mathbf{C}^{*}_{s}(\mathbf{C}_{s}\mathbf{C}_{s}^{*})^{-1}\mathbf{e}_{j}\rangle^{*}_{T}$ where $\mathbf{e}_{j}$ is the $j^{th}$ standard basis vector in $\mathbb{C}^{s}$. The function $\xi\mapsto\langle z_{\xi},\boldsymbol{\phi}_{t}\rangle=\sum_{j=1}^{s}\mathbf{z}_{\xi}^{*}(j)e^{-2\pi i\xi_{j}t}$ is equal in the $L_{2}(\mu)$ sense to $\mathbf{Z}^{*}\boldsymbol{\phi}_{t}$. Let us define:

[TABLE]

For a fixed $t$, consider the function $\xi\mapsto r_{\xi}(t)$, which we denote by $r_{\cdot}(t)$. We have $r_{\cdot}(t)=\varphi_{t}-\mathbf{Z}^{*}\boldsymbol{\phi}_{t}$, again in the $L_{2}(\mu)$ sense. Thus, we can write

[TABLE]

Further, for a fixed $\xi$, if we consider the function $t\mapsto{r}_{\xi}(t)$, which we denote by ${r}_{\xi}(\cdot)$, we notice that it is a $s+1=O(s_{\mu,\epsilon})$ sparse Fourier function, so applying Lemma 12 we have for any $\xi\in\mathbb{R}$ and $t\in[0,T]$:

[TABLE]

Combining (24) with (5.1) we can thus bound for any $t\in[0,T]$:

[TABLE]

where the last bound again follows from (5.1). By Theorem 9 we have $\frac{1}{T}\int_{w\in[0,T]}\|\varphi_{w}-\mathbf{Z}^{*}\boldsymbol{\phi}_{w}\|_{\mu}^{2}\,dw\leq 4\epsilon\cdot s_{\mu,\epsilon}$. Plugging into (5.1) and using that we can choose $s\leq 36\cdot s_{\mu,\epsilon}+1$, for all $t\in[0,T]$:

[TABLE]

Combining (5.1) and (26) completes the proof of Lemma 11 since $\bar{\tau}_{\mu,\epsilon}(t)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}2\cdot\left(\ell(t)+\frac{r(t)}{\epsilon}\right)$ and thus

[TABLE]

∎

Theorem 8 gives a universal uniform bound on the ridge leverage scores corresponding to measure $\mu$ in terms of $s_{\mu,\epsilon}$. If we directly sample time points according to the uniform distribution over $[0,T]$, this theorem shows that $\operatorname{poly}(s_{\mu,\epsilon})$ samples and $\operatorname{poly}(s_{\mu,\epsilon})$ runtime suffice to apply Theorem 7 and solve Problem 1 with good probability. This is already a surprising result, showing that the simplest sampling scheme, uniform random sampling, can give bounds in terms of the optimal complexity $s_{\mu,\epsilon}$ for any $\mu$. Existing methods with similar complexity, such as those that interpolate bandlimited signals using prolate spheroidal wave functions [XRY01, STR06] require nonuniform sampling. Methods that use uniform sampling, such as truncated Whittaker-Shannon, have sample complexity depending polynomially rather than logarithmically on the desired error $\epsilon$.

5.2 Gap-based leverage score bound

Our final result gives a much tighter bound on the ridge leverage scores than the uniform bound of Theorem 8. The key idea is to show that the bound is loose for $t$ bounded away from the edges of $[0,T]$. Specifically we have:

Theorem 13** (Gap-Based Leverage Score Bound).**

For all $t$,

[TABLE]

Proof.

Consider $t\in[0,T/2]$. We will show that $\tau_{\mu,\epsilon}(t)\leq\frac{s_{\mu,\epsilon}}{t}.$ A symmetric proof will hold for $t\in[T/2,T]$, giving the theorem. We define an auxiliary operator: $\mathcal{F}_{\mu,t}:L_{2}(T)\rightarrow L_{2}(\mu)$ which is given by restricting the integration in $\mathcal{F}_{\mu}$ to $[0,t]$. Specifically, for $f\in L_{2}(T)$ we have:

[TABLE]

We can see that $[\mathcal{F}_{\mu,t}^{*}g](s)=\int_{\mathbb{R}}g(\xi)e^{2\pi is\xi}\,d\mu(\xi)$ for $s\in[0,t]$ and $[\mathcal{F}_{\mu,t}^{*}g](s)=0$ for $s\in(t,T]$. We will use the leverage score of some $s\in[0,t]$ in the restricted operator $\mathcal{F}_{\mu,t}$ to upper bound those of $t$ in $\mathcal{F}_{\mu}$. We start by defining these scores analogously to Definition 3 for $\mathcal{F}_{\mu}$.

Definition 4** (Restricted ridge leverage scores).**

For probability measure $\mu$ on $\mathbb{R}$, time length $T$, $t\in[0,T]$ and $\epsilon\geq 0$, define the $\epsilon$-ridge leverage score of $s\in[0,t]$ in $\mathcal{F}_{\mu,t}$ as:

[TABLE]

We have the following leverage score properties, analogous to those given for $\mathcal{F}_{\mu}$ in Theorem 5:

Theorem 14** (Restricted leverage score properties).**

Let $\tau_{\mu,\epsilon,t}(s)$ be as defined in Definition 4.

• •

The leverage scores integrate to the statistical dimension:

[TABLE]

• •

Inner Product Characterization: Letting $\varphi_{s}\in L_{2}(\mu)$ have $\varphi_{s}(\xi)=e^{-2\pi is\xi}$ for $s\in[0,t]$,

[TABLE]

• •

Minimization Characterization:

[TABLE]

We first show that the restricted leverage scores of Definition 4 are not too large on average.

Claim 15** (Restricted statistical dimension bound).**

[TABLE]

Proof.

Via (28) we have $\int_{0}^{t}\tau_{\mu,\epsilon,t}(s)\,ds=s_{\mu,\epsilon,t}$ which we can write as:

[TABLE]

From Claim 35 we have $\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*}\preceq\mathcal{F}_{\mu}\mathcal{F}_{\mu}^{*}=\mathcal{G}_{\mu}$. Since $(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*}+\epsilon{\cal I}_{\mu})^{-1/2}(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*})(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*}+\epsilon{\cal I}_{\mu})^{-1/2}={\cal I}_{\mu}-\epsilon(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*}+\epsilon{\cal I}_{\mu})^{-1}$ and $(\mathcal{G}_{\mu}+\epsilon{\cal I}_{\mu})^{-1/2}\mathcal{G}_{\mu}(\mathcal{G}_{\mu}+\epsilon{\cal I}_{\mu})^{-1/2}={\cal I}_{\mu}-\epsilon(\mathcal{G}_{\mu}+\epsilon{\cal I}_{\mu})^{-1}$ we have from Claim 27 $(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*}+\epsilon{\cal I}_{\mu})^{-1/2}(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*})(\mathcal{F}_{\mu,t}\mathcal{F}_{\mu,t}^{*}+\epsilon{\cal I}_{\mu})^{-1/2}\preceq(\mathcal{G}_{\mu}+\epsilon{\cal I}_{\mu})^{-1/2}\mathcal{G}_{\mu}(\mathcal{G}_{\mu}+\epsilon{\cal I}_{\mu})^{-1/2}$, and since the trace is monotone for trace-class operators, we have

[TABLE]

which gives the claim. ∎

From Claim 15 we immediately have:

Claim 16**.**

There exists $s^{\star}\in[0,t]$ with $\tau_{\mu,\epsilon,t}(s^{\star})\leq\frac{s_{\mu,\epsilon}}{t}$.

Proof.

Assume for the sake of contradiction that $\tau_{\mu,\epsilon,t}(s)>\frac{s_{\mu,\epsilon}}{t}$ for all $s\in[0,t]$. The by (28),

[TABLE]

This contradicts Claim 15, giving the claim. ∎

We now show that the leverage score of $s^{\star}$ in $\mathcal{F}_{\mu,t}$ upper bounds the leverage score of $t$ in $\mathcal{F}_{\mu}$, completing the proof of Theorem 13. We apply the minimization characterization of Theorem 14, equation (30), showing that by simply shifting an optimal solution for $s^{\star}$ we can show the existence of a good solution for $t$, upper bounding its leverage score by that of $s^{\star}$ and giving $\tau_{\mu,\epsilon}(t)\leq\tau_{\mu,\epsilon,t}(s^{\star})\leq\frac{s_{\mu,\epsilon}}{t}$ by Claim 16.

Formally, by Claim 16 and (30), there is some $\beta^{\star}\in L_{2}(T)$ achieving:

[TABLE]

We can assume without loss of generality that $\beta^{\star}(s)=0$ for $s\notin[0,t]$, since $\mathcal{F}_{\mu,t}\beta^{\star}$ is unchanged if we set $\beta^{\star}(s)=0$ on this range and since doing this cannot increase $\|\beta\|_{T}^{2}$. Now, let $\bar{\beta}\in L_{2}(T)$ be given by $\bar{\beta}(s)=\beta^{\star}(s-(t-s^{\star}))$. That is, $\bar{\beta}$ is just $\beta^{\star}$ shifted from the range $[0,t]$ to the range $[t-s^{\star},2t-s^{\star}]$. Note that since we are assuming $t\leq T/2$, $[t-s^{\star},2t-s^{\star}]\subset[0,T]$. For any $\xi$:

[TABLE]

Now,

[TABLE]

Combined with (33) this gives:

[TABLE]

Finally, noting that $\|\bar{\beta}\|_{T}=\|\beta^{\star}\|_{T}$ and applying the minimization characterization of Theorem 5, the bound in (5.2) along with (32) gives:

[TABLE]

which completes the theorem.

∎

5.3 Nearly tight leverage score bound

Combining Theorems 8 and 13 gives our tight, spectrum blind leverage score bound:

Theorem 17** (Spectrum Blind Leverage Score Bound).**

For any $\alpha,T\geq 0$ let $\tilde{\tau}_{\alpha}(t)$ be given by:

[TABLE]

For any probability measure $\mu$, $T\geq 0$, $0\leq\epsilon\leq 1$ and $t\in[0,T]$, if $\alpha\geq 256\cdot s_{\mu,\epsilon}$:

[TABLE]

A visualization of $\tilde{\tau}_{\alpha}$ is given in Figure 3.

Proof.

The fact that $\tau_{\mu,\epsilon}(t)\leq\tilde{\tau}_{\alpha}(t)$ follows from Theorems 8 and 13:

• •

For $t\in[T/\alpha^{6},T(1-1/\alpha^{6})]$, by Theorem 13 if $\alpha\geq 256\cdot s_{\mu,\epsilon}$ we have

[TABLE]

• •

For $t\in[0,T/\alpha^{6}]\cup[T(1-1/\alpha^{6}),T]$, by Theorem 8 we can bound,

[TABLE]

for $\alpha\geq 256\cdot s_{\mu,\epsilon}$. Note that the second inequality uses that $\log^{3}(40x)\leq 64x$ for any $x$.

The integral of the approximate scores $\tilde{s}_{\alpha}$ is bounded as:

[TABLE]

where the last inequality follows since for $\epsilon\leq 1$, $s_{\mu,\epsilon}\geq 1/2$ and so $\alpha\geq 128$. ∎

5.4 Putting it all together: generic signal reconstruction

Finally, we combine the leverage score bound of Theorem 17 with Theorem 7 to give our main algorithmic result, Theorem 3 (and as a corollary, Theorem 2). We state the full theorem below:

Theorem 3 (Main result, algorithmic complexity).

Consider any measure $\mu$, for which we can compute the kernel function $k_{\mu}(t_{1},t_{2})=\int_{\xi\in\mathbb{R}}e^{-2\pi i(t_{1}-t_{2})}d\mu(\xi)$ for any $t_{1},t_{2}\in[0,T]$ in time $Z$.

Let $\tilde{\tau}_{\alpha}(t)$ be as defined in Theorem 17. For any $\epsilon\leq\|\mathcal{K}_{\mu}\|_{\mathrm{op}}$ and $T>0$, let $\tilde{\tau}_{\mu,\epsilon}(t)=\tilde{\tau}_{\alpha}(t)$ for $\alpha=\beta\cdot s_{\mu,\epsilon}$ with $\beta\geq 256$. Algorithm 1 applied with $\tilde{\tau}_{\mu,\epsilon}(t)$ and failure probability $\delta$ returns $t_{1},\ldots,t_{s}\in[0,T]$ and $\mathbf{z}\in\mathbb{C}^{s}$ such that $\tilde{y}(t)=\sum_{i=1}^{s}\mathbf{z}(i)\cdot k_{\mu}(t_{i},t)$ solves Problem 1 with parameter $6\epsilon$ and probability $\geq 1-\delta$. That is, with probability of at least $1-\delta$:

[TABLE]

The algorithm queries $y+n$ at $s$ points and runs in $O\left(s^{2}\cdot Z+s^{\omega}\right)$ time where

[TABLE]

The output $\tilde{y}(t)$ can be evaluated in $O(s\cdot Z)$ time for any $t$ using Algorithm 2.

Note that if we want to solve Problem 1 with parameter $\epsilon$, it suffices to apply Theorem 3 with parameter $\epsilon^{\prime}=\epsilon/6$. The asymptotic complexity will be identical since, by (17), $s_{\mu,\epsilon/6}\leq 6s_{\mu,\epsilon}$.

Proof.

The theorem follows directly from Theorem 7, along with Theorem 17 which shows that, for $\alpha=\beta\cdot s_{\mu,\epsilon}$ with $\beta\geq c_{1}$ and $\tilde{\tau}_{\mu,\epsilon}(t)=\tilde{\tau}_{\alpha}(t)$ we have:

$\tilde{\tau}_{\mu,\epsilon}(t)\geq\tau_{\mu,\epsilon}(t)$ for all $t\in[0,T]$. 2. 2.

$\tilde{s}_{\mu,\epsilon}=\int_{0}^{T}\tilde{\tau}_{\mu,\epsilon}(t)\,dt=O\left(\beta\cdot s_{\mu,\epsilon}\log(\beta\cdot s_{\mu,\epsilon})\right)$.

The runtime bound follows after noting that we can sample according to $\tau_{\alpha}$ in $W=O(1)$ time using inverse transform sampling since it is straightforward to derive an explicit expression for the CDF and compute the inverse (see (35)).

∎

6 Lower bound

We conclude by showing that the statistical dimension $s_{\mu,\epsilon}$ tightly characterizes the sample complexity of solving Problem 1, under a very mild assumption on $\mu$ that holds for all natural constraints we discuss in this paper. Thus, Theorem 1 is tight up to logarithmic factors.

We first define a quantity, $n_{\mu,\epsilon}$ that gives a natural lower bound on $s_{\mu,\epsilon}$. For any $\mu,\epsilon$, let

[TABLE]

That is, $n_{\mu,\epsilon}$ is the number of eigenvalues of $\mathcal{K}_{\mu}$ that are larger than $\epsilon$. As shown in (19), we always have $n_{\mu,\epsilon}\leq 2s_{\mu,\epsilon}$. We first prove that solving Problem 1 requires $\Omega(n_{\mu,\epsilon})$ samples. We then show that, under a very mild constraint on $\mu$ (which holds for all $\mu$ we consider including sparse, bandlimited, multiband, Gaussian, and Cauchy-Lorentz), $n_{\mu,\epsilon}=\Omega(s_{\mu,\epsilon})$. Thus, $s_{\mu,\epsilon}$ gives a tight bound on the query complexity of solving Problem 1.

Theorem 18** (Lower bound in terms of eigenvalue count).**

Consider a measure $\mu$, an error parameter $\epsilon>0$, and any (possibly randomized) algorithm that solves Problem 1 with probability $\geq 2/3$ for any function $y$ and makes at most $r$ (possibly adaptive) queries on any input. Then $r\geq n_{\mu,72\epsilon}/20$.

Proof.

We describe a distribution on inputs $y$ on which any deterministic algorithm that takes $r=o(n_{\mu,72\epsilon})$ samples on any input fails with probability $>1/3$. The theorem then follows by Yao’s principle.

Notation: Let $v_{1},\ldots,v_{n_{\mu,72\epsilon}}\in L_{2}(\mu)$ be the eigenfunctions of $\mathcal{G}_{\mu}$ corresponding to its top $n_{\mu,72\epsilon}$ eigenvalues. Let $\mathbf{Z}:L_{2}(\mu)\rightarrow\mathbb{C}^{n_{\mu,72\epsilon}}$ be the operator with $v_{i}$ as its $i^{th}$ row – i.e., $[\mathbf{Z}g](i)=\langle v_{i},g\rangle_{\mu}$. Note that $\mathbf{Z}$ has orthonormal rows. Let $\mathbf{D}\in\mathbb{R}^{n_{\mu,72\epsilon}\times n_{\mu,72\epsilon}}$ be a diagonal matrix with $\mathbf{D}_{ii}=\sqrt{\lambda_{i}(\mathcal{K}_{\mu})}$. Let $\mathbf{U}=\mathcal{F}_{\mu}^{*}\mathbf{Z}^{*}\mathbf{D}^{-1}$. We can see that $\mathbf{Z}\mathcal{F}_{\mu}\mathcal{F}_{\mu}^{*}\mathbf{Z}^{*}=\mathbf{Z}\mathcal{G}_{\mu}\mathbf{Z}^{*}=\mathbf{D}^{2}$ and hence, $\mathbf{U}^{*}\mathbf{U}=\mathbf{D}^{-1}\mathbf{Z}\mathcal{F}_{\mu}\mathcal{F}_{\mu}^{*}\mathbf{Z}^{*}\mathbf{D}^{-1}=\mathbf{I}$. While not needed for our proof, we can check that $\mathbf{U}:\mathbb{C}^{n_{\mu,72\epsilon}}\rightarrow L_{2}(T)$ is an operator with columns corresponding to all eigenfunctions of $\mathcal{K}_{\mu}$ with eigenvalue $\geq 72\epsilon$.

Hard Input Distribution: Let $\mathbf{c}\in\mathbb{R}^{n_{\mu,72\epsilon}}$ be a random vector with each entry distributed independently as a Gaussian: $\mathbf{c}(i)\sim\mathcal{N}(0,\frac{1}{n_{\mu,72\epsilon}})$. Let $\mathbf{\bar{c}}=\mathbf{D}^{-1}\mathbf{c}$, $x=\mathbf{Z}^{*}\mathbf{\bar{c}}$, and the random input be $y=\mathcal{F}_{\mu}^{*}x$. That is, $y=\mathcal{F}_{\mu}^{*}\mathbf{Z}^{*}\mathbf{D}^{-1}\mathbf{c}=\mathbf{U}\mathbf{c}$ is a random linear combination of the top eigenfunctions of $\mathcal{K}_{\mu}$. While, formally, $\mathcal{F}_{\mu}^{*}x\in L_{2}(T)$ is an equivalence class of functions, since our input model requires that $y$ admits pointwise evaluation, we will abuse notation, letting $y$ denote the member of this class with $y(t)=\langle\varphi_{t},\mathbf{Z}^{*}\mathbf{D}^{-1}\mathbf{c}\rangle_{\mu}=\langle\mathbf{D}^{-1}\mathbf{Z}\varphi_{t},\mathbf{c}\rangle$, where $\varphi_{t}(\xi)=e^{-2\pi it\xi}$.

We prove that accurately reconstructing $y$ drawn from the hard input distribution yields an accurate reconstruction of the random vector $\mathbf{c}$. Since $\mathbf{c}$ is $n_{\mu,72\epsilon}$ dimensional, this reconstruction requires $\Omega(n_{\mu,72\epsilon})$ samples, giving us a lower bound for accurately reconstructing $y$.

Claim 19**.**

For random $x$ distributed as described above, with probability $\geq 5/6$, $\|x\|_{\mu}^{2}\leq\frac{1}{12\epsilon}$.

Proof.

[TABLE]

We then bound $\|\mathbf{\bar{c}}\|_{2}^{2}\leq\|\mathbf{c}\|_{2}^{2}/\lambda_{n_{\mu,72\epsilon}}(\mathcal{K}_{\mu})\leq\frac{\|\mathbf{c}\|_{2}^{2}}{72\epsilon}$ since $\lambda_{n_{\mu,72\epsilon}}(\mathcal{K}_{\mu})\geq 72\epsilon$ by definition. Finally, note that $\|\mathbf{c}\|_{2}^{2}$ is a Chi-squared random variable, with $\mathbb{E}[\|\mathbf{c}\|_{2}^{2}]=1$. So loosely, by Markov’s inequality, with probability $\geq 5/6$, ${\|\mathbf{c}\|_{2}^{2}}\leq 6$, which gives the claim. ∎

From Claim 19 we have:

Claim 20**.**

Given random input $y=\mathcal{F}_{\mu}^{*}x$ generated as described above, with probability $\geq 5/6$, to solve Problem 1, an algorithm must return a representation of $\tilde{y}$ with $\|y-\tilde{y}\|_{T}^{2}\leq\frac{1}{12}$.

Proof.

Solving Problem 1 requires finding a representation of $\tilde{y}$ with $\|y-\tilde{y}\|_{T}^{2}\leq\epsilon\|x\|_{\mu}^{2}+C\|n\|_{T}^{2}$. By Claim 19 and the fact that for our input $\|n\|_{T}^{2}=0$, with probability $\geq 5/6$ one has that $\epsilon\|x\|_{\mu}^{2}+C\|n\|_{T}^{2}\leq\frac{1}{12}$, yielding the claim. ∎

We next show that finding a $\tilde{y}$ satisfying the condition of Claim 20 is at least as hard as finding an accurate approximation to $\mathbf{c}$.

Claim 21**.**

For $\tilde{y}$ with $\|y-\tilde{y}\|_{T}^{2}\leq\frac{1}{12}$, $\mathbf{\tilde{c}}=\mathbf{U}^{*}\tilde{y}$ satisfies $\|\mathbf{c}-\mathbf{\tilde{c}}\|_{2}^{2}\leq\frac{1}{12}$.

Proof.

Recalling that $y=\mathbf{U}\mathbf{c}$, for $\mathbf{\tilde{c}}=\mathbf{U}^{*}\tilde{y}$ we have:

[TABLE]

Recalling that $\mathbf{U}^{*}\mathbf{U}=\mathbf{I}$ we thus have:

[TABLE]

The second to last inequality follows since $\mathbf{U}^{*}\mathbf{U}=\mathbf{I}$ and $\mathbf{U}\mathbf{U}^{*}$ are finite rank, so are compact and share the same non-zero eigenvalues. Thus, $\mathbf{U}\mathbf{U}^{*}\preceq{\cal I}_{T}$ [HN01, Lemma 8.26]. This completes the claim. ∎

Combining Claims 20 and 21 we have:

Claim 22**.**

If a deterministic algorithm solves Problem 1 with probability $\geq 2/3$ over our random input $y=\mathbf{U}\mathbf{c}$, then with probability $\geq 1/2$, letting $\tilde{y}$ be the output of the algorithm, $\mathbf{\tilde{c}}=\mathbf{U}^{*}\tilde{y}$ satisfies $\|\mathbf{c}-\mathbf{\tilde{c}}\|_{2}^{2}\leq\frac{1}{12}$.

Proof.

If an algorithm solves Problem 1 probability $\geq 2/3$ then by Claim 20, it returns $\tilde{y}$ with $\|y-\tilde{y}\|_{T}^{2}\leq\frac{1}{12}$ with probability $\geq 2/3-1/6=1/2$. Thus, by Claim 21, $\mathbf{\tilde{c}}$ satisfies $\|\mathbf{c}-\mathbf{\tilde{c}}\|_{2}^{2}\leq\frac{1}{12}$ with probability $\geq 1/2$. ∎

Finally, we complete the proof of Theorem 18 by arguing that if $\tilde{y}$ is formed using $o(n_{\mu,72\epsilon})$ queries, then for $\mathbf{\tilde{c}}=\mathbf{U}^{*}\tilde{y}$, $\|\mathbf{c}-\mathbf{\tilde{c}}\|_{2}^{2}>\frac{1}{12}$ with good probability. Thus the bound of Claim 22 cannot hold and so $\tilde{y}$ cannot be a solution to Problem 1 with good probability.

Assume for the sake of contradiction that there is a deterministic algorithm solving Problem 1 with probability $\geq 2/3$ over the random input $\mathbf{U}\mathbf{c}$ that makes $r=\frac{n_{\mu,72\epsilon}}{20}$ queries on any input (note that if there exists an algorithm that makes fewer queries on some inputs, we can always modify it to make exactly $\frac{n_{\mu,72\epsilon}}{20}$ queries and return the same output.)

As discussed, each query to $y$ is a query to $y(t)=\langle\mathbf{D}^{-1}\mathbf{Z}\varphi_{t},\mathbf{c}\rangle$. Consider a deterministic function $Q$, that is given input $\mathbf{V}\in\mathbb{C}^{i\times n_{\mu,72\epsilon}}$ (for any positive integer $i$) and outputs $Q(\mathbf{V})\in\mathbb{C}^{n_{\mu,72\epsilon}\times n_{\mu,72\epsilon}}$ such that $Q(\mathbf{V})$ has orthonormal rows with the first $i$ spanning the $i$ rows of $\mathbf{V}$. For example, $Q$ may run Gram-Schmidt orthogonalization on $\mathbf{V}$ fixing its first $\operatorname{rank}(\mathbf{V})\leq i$ rows and then fill out the remaining $n_{\mu,72\epsilon}-\operatorname{rank}(\mathbf{V})$ rows using some canonical approach. Letting $\mathbf{D}^{-1}\mathbf{Z}\varphi_{t_{1}},\ldots,\mathbf{D}^{-1}\mathbf{Z}\varphi_{t_{r}}$ denote the queries made by our algorithm on random input $\mathbf{c}$, let $\mathbf{Q}^{i}=Q([\mathbf{D}^{-1}\mathbf{Z}\varphi_{t_{1}},\ldots,\mathbf{D}^{-1}\mathbf{Z}\varphi_{t_{i}}]^{*})$. That is $\mathbf{Q}^{i}$ is an orthonormal matrix whose first $i$ rows span our first $i$ queries. Note that since our algorithm is deterministic, $\mathbf{Q}^{i}$ is a deterministic function of the random input $\mathbf{c}$. We have the following claim:

Claim 23**.**

Conditioned on the queries $y(t_{1}),\ldots y(t_{r})$, for $j>r$, each $[\mathbf{Q}^{r}\mathbf{c}](j)$ is distributed independently as $\mathcal{N}\left(0,\frac{1}{n_{\mu,72\epsilon}}\right)$.

Proof.

We prove the claim via induction on the number of queries considered. For the base case set $i=1$. $\mathbf{Q}^{1}$ is a deterministic matrix (since the choice of our first query is made determinstically before seeing any input) and so by the rotational invariance of the Gaussian distribution, the entries of $\mathbf{Q}^{1}\mathbf{c}$ are distributed independently as $\mathcal{N}\left(0,\frac{1}{n_{\mu,72\epsilon}}\right)$ (the same as the entries of $\mathbf{c}$). The first row of $\mathbf{Q}^{1}$ spans our first query, and thus this row is just equal to $\mathbf{D}^{-1}\mathbf{Z}\varphi_{t_{1}}$ scaled to have unit norm. Thus $y(t_{1})=\mathbf{D}^{-1}\mathbf{Z}\varphi_{t_{1}}\mathbf{c}$ is just a fixed scaling of $[\mathbf{Q}^{1}\mathbf{c}](1)$. So conditioning on $y(t_{1})$, we still have [$\mathbf{Q}^{1}\mathbf{c}](j)$ for $j>1$ distributed independently as $\mathcal{N}\left(0,\frac{1}{n_{\mu,72\epsilon}}\right)$.

Now, consider $i>1$. By the inductive assumption, conditioned on $y(t_{1}),\ldots y(t_{i-1})$, for $j\geq i$, $[\mathbf{Q}^{i-1}\mathbf{c}](j)$, are distributed independently as $\mathcal{N}\left(0,\frac{1}{n_{\mu,72\epsilon}}\right)$. We can see that both $\mathbf{Q}^{i-1}$ and $\mathbf{Q}^{i}$ are fixed conditioned on $y(t_{1}),\ldots y(t_{i-1})$ (since the $i^{th}$ query is chosen deterministically, possibly adaptively as a function of the previously seen queries $y(t_{1}),\ldots y(t_{i-1})$). Additionally, since they share their first $i-1$ rows, the remaining $n_{\mu,72\epsilon}-i+1$ rows of $\mathbf{Q}^{i-1}$ and $\mathbf{Q}^{i}$ have the same rowspans. Thus we can write $\mathbf{Q}^{i}=[\mathbf{I};\mathbf{R}]\mathbf{Q}^{i-1}$ where $\mathbf{R}\in\mathbb{C}^{n_{\mu,72\epsilon}-i+1\times n_{\mu,72\epsilon}-i+1}$ is some fixed rotation with $\mathbf{R}^{*}\mathbf{R}=\mathbf{I}$. Thus, by the rotational invariance of the Gaussian, for all $j\geq i$, $[\mathbf{Q}^{i}\mathbf{c}](j)$ are distributed independently as $\mathcal{N}\left(0,\frac{1}{n_{\mu,72\epsilon}}\right)$ (the same as $[\mathbf{Q}^{i-1}\mathbf{c}](j)$). Further conditioning on $y(t_{i})$, which is a deterministic function of $[\mathbf{Q}^{i}\mathbf{c}](i)$ and $y(t_{1})\ldots y(t_{i-1})$, we still have that for $j>i$, $[\mathbf{Q}^{i}\mathbf{c}](j)$ are distributed independently as $\mathcal{N}\left(0,\frac{1}{n_{\mu,72\epsilon}}\right)$. This completes the inductive step and so the claim. ∎

Armed with Claim 23 we can compute:

[TABLE]

where the last line follows since, conditioned on $y(t_{1})...,y(t_{r})$, $\mathbf{Q}^{r}\mathbf{\tilde{c}}$ is fixed and for $i\geq r+1$, $\mathbf{Q}^{r}\mathbf{c}(i)$ are distributed independently as Gaussians centered around [math] (by Claim 23). So the probability of the sum of differences being small is only smaller than if we replaced each $\mathbf{Q}^{r}\mathbf{\tilde{c}}(i)$ by [math].

Now, conditioned on $y(t_{1})...,y(t_{r})$, $\sum_{i=r+1}^{n_{\mu,72\epsilon}}\left|[\mathbf{Q}^{r}\mathbf{c}](i)\right|^{2}$ is a Chi-squared random variable with

[TABLE]

For $r=\frac{n_{\mu,72\epsilon}}{20}$, we thus have $\mathbb{E}\left[\sum_{i=r+1}^{n_{\mu,72\epsilon}}\left|[\mathbf{Q}^{r}\mathbf{c}](i)\right|^{2}\mid y(t_{1})...,y(t_{r})\right]\geq\frac{19}{20}$. We can loosely upper bound the probability in (37), using that for a Chi-squared random variable $X$ with $k$ degrees of freedom, $\Pr[X\leq\delta\mathbb{E}[X]]\leq(\delta e^{1-\delta})^{k/2}\leq(\delta e^{1-\delta})^{1/2}$. So,

[TABLE]

Plugging back into (37) gives:

[TABLE]

However, we have assumed that our algorithm solves Problem 1 with probability $\geq 2/3$, and hence, by Claim 22, $\Pr\left[\|\mathbf{c}-\mathbf{\tilde{c}}\|_{2}^{2}\leq\frac{1}{12}\right]\geq\frac{1}{2}$. This is a contradiction, yielding the theorem.

∎

6.1 Statistical Dimension Lower Bound

We now use Theorem 18 to prove that the statistical dimension tightly characterizes the sample complexity of solving Problem 1 for any constraint measure $\mu$ satisfying a simple condition: we must have $s_{\mu,\epsilon}=O(1/\epsilon^{p})$ for some $p<1$. Note that this assumption holds for all $\mu$ considered in this work (including bandlimited, multiband, sparse, Gaussian, and Cauchy-Lorentz), where $s_{\mu,\epsilon}$ either grows as $\log(1/\epsilon)$ or $1/\sqrt{\epsilon}$. Also note that by (5) we can always bound $s_{\mu,\epsilon}\leq\operatorname{tr}(\mathcal{K}_{\mu})/\epsilon=1/\epsilon$. So this assumption holds whenever we have a nontrivial upper bound on $s_{\mu,\epsilon}$.

Theorem 24** (Statistical Dimension Lower Bound).**

For any probability measure $\mu$, suppose that $s_{\mu,\epsilon}=O(1/\epsilon^{p})$ for some constant $p<1$. Consider any (possibly randomized) algorithm that solves Problem 1 with probability $\geq 2/3$ for any function $y$ and any $\epsilon>0$ and makes at most $r_{\mu,\epsilon}$ (possibly adaptive) queries on any input. Then $r_{\mu,\epsilon}=\Omega(s_{\mu,\epsilon})$.232323Here we follow the Hardy-Littlewood definition [HL14], using $f(\epsilon)=\Omega(g(\epsilon))$ to denote that $\limsup_{x\to\infty}\frac{f(\epsilon)}{g(\epsilon)}>0$. Thus the lower bound shows that, for some fixed constant $c>0$, for every $\epsilon$, there is at least some $\epsilon^{\prime}<\epsilon$ where the number of queries used by any algorithm solving Problem 1 with probability $\geq 2/3$ is at least $c\cdot s_{\mu,\epsilon}$. In other words, the lower bound rules out the possibility that the number of queries is $o(s_{\mu,\epsilon})$.

Proof.

We simply prove that for this class of measures, $n_{\mu,72\epsilon}=\Omega(s_{\mu,\epsilon})$ and then apply Theorem 18. It suffices to show that $n_{\mu,\epsilon}=\Omega(s_{\mu,c\epsilon})$ for any fixed constant $c\geq 1$ since by (17), $s_{\mu,c\epsilon}\geq\frac{s_{\mu,\epsilon}}{c}$. Thus $n_{\mu,\epsilon}=\Omega(s_{\mu,c\epsilon})$ gives that $n_{\mu,72\epsilon}=\Omega(s_{\mu,72c\epsilon})=\Omega(s_{\mu,\epsilon})$, giving the theorem.

Let $c_{p}=2^{\frac{4}{1-p}}>1$. Assume for the sake of contradiction that $n_{\mu,\epsilon}=o(s_{\mu,c_{p}\epsilon})$. By this assumption, there is some fixed $\epsilon_{0}$ such that,

[TABLE]

We can bound:

[TABLE]

and thus by (38) have for any $\epsilon\leq\epsilon_{0}$:

[TABLE]

Now we also have:

[TABLE]

Combined with (39) this gives that for any $\epsilon\leq\epsilon_{0}$:

[TABLE]

By (40) we in turn have that, for every $\epsilon\leq\epsilon_{0}$,

[TABLE]

Using that $\lfloor\log_{c_{p}}\epsilon_{0}/\epsilon\rfloor\geq\log_{c_{p}}\epsilon_{0}/\epsilon-1$ and that $c_{p}=2^{\frac{4}{1-p}}\geq 16$ we can then bound, for all $\epsilon\leq\epsilon_{0}$:

[TABLE]

Note that $\left(\frac{c_{p}}{4}\right)^{\log_{c_{p}}\epsilon_{0}-1}$ is a constant independent of $\epsilon$. Thus, the above contradicts the assumption that $s_{\mu,\epsilon}=O(1/\epsilon^{p})$, giving the theorem. ∎

Remark We remark that a similar technique to Theorem 24 can be used to show that $n_{\mu,\epsilon}=\Omega(s_{\mu,\epsilon}/\epsilon^{p})$ for any $p>0$, without any assumptions on $s_{\mu,\epsilon}$.

7 Conclusion and Open Problems

We view our work as the starting point for further exploring the application of techniques from the randomized numerical linear algebra literature (such as leverage score sampling, column based matrix reconstruction, and random projection) in signal processing. We lay out a number of open directions that we consider interesting below:

• •

The most immediate question is to generalize our results for interpolation over an interval to higher dimensional spaces. Fourier constrained interpolation in two or three dimensions is important in many areas, such as the earth and geosciences [Rip89, Rip05] and image processing [PPV02, RVU06]. Interpolation in even higher dimensions is common in Gaussian process methods in machine learning. We believe that our techniques should extend to higher dimensions in a similar manner to prior related work on kernel approximation [AKM*+*17].

• •

We have considered a simple signal reconstruction problem, where we wish to reconstruct a function over a fixed interval given sample access at points in that interval. There are many interesting variations of this problem. For example, can better bounds be achieved if samples can be taken from the interval $[0,T]$, but we only consider reconstruction error over a subset of this interval? In this setting, can uniform sampling give optimal bounds? How can one formulate a similar reconstruction problem and adapt our techniques to the streaming setting, where we hope to estimate a signal at any given point in time using measurements at past samples (and perhaps must limit memory/computation at any given time)? Can our techniques be extended to the setting where the error is averaged using a non-uniform measure in time domain? This question is especially relevant for applications to machine learning, where we may wish to approximate the signal well on average on input points drawn from some non-uniform distribution. In traditional supervised learning, reconstruction would be performed using points drawn from this same distribution. However, in an active learning setting, we may be allowed to drawn points from some other distribution, such as the leverage score distribution, which yields better error.

• •

In our work we have assumed knowledge of the constraint $\mu$. However, as discussed, in the case of sparse and multiband signal reconstruction, it is important to learn $\mu$ (i.e., the locations of the frequencies or frequency bands) as part of the reconstruction process. Understanding how to do this, perhaps by combining existing techniques [ME09, Moi15, PS15, CKPS16] with our own is an important direction. More generally, in many applications, $\mu$ is derived from the signal itself, by estimating the signal’s autocorrelation, which corresponds to our kernel function $k_{\mu}$. Can our techniques be used to give bounds in this setting?

• •

Can our techniques be extended to learning signals giving constraints on other transforms such as the short-time Fourier transform (the signal’s spectrogram), the wavelet transform, etc.? More generally, can leverage score sampling be used to approximate these transforms and to approximately apply filters or other signal modifications based on them?

• •

What is the connection between our randomized leverage score sampling method and deterministic ‘sampling’ methods such as Chebyshev interpolation for low-degree polynomials, uniform sampling for bandlimited signal reconstruction, and non-uniform “multicoset” sampling schemes considered in the signal processing literature [FB96, VB00, Bre08, ME09]. Can our results be made deterministic, perhaps using deterministic sampling methods for operator approximation like those employed in our proof of Lemma 46?

Acknowledgements

We thank Ron Levie for helpful discussions on weak integrals in Hilbert spaces and Zhao Song for discussions on smoothness bounds for sparse Fourier functions. We also thank Yonina Eldar for helpful discussions and pointers to related work.

Haim Avron’s work is supported in part by Israel Science Foundation (grant no. 1272/17) and United States-Israel Binational Science Foundation (grant no. 2017698). Michael Kapralov is supported in part by ERC Starting Grant 759471.

Appendix A Prior work on Fourier constrained interpolation

As mentioned in Section 1, constrained interpolation problems similar to Problem 1 have been studied for decades in a number of different communities, often with widely varying computational models, assumptions, and goals. We discuss the most relevant prior work here.

Bandlimited functions. The most well studied special case of Problem 1 is when $\mu$ is uniform on an interval $[-F,F]$, which corresponds to reconstructing a bandlimited function from discrete samples. Work on this problem goes back to famous results of Nyquist, Shannon, and others [Whi15, Nyq28, Kot33, Sha49], who showed that it suffices to sample time uniformly with frequency $O(1/F)$. While this rate roughly suggests that $O(FT)$ samples should be required to reconstruct a signal on $[0,T]$, this does not follow directly: common reconstruction methods like Whittaker-Shannon sinc intetrpolation infer $y(t)$ from an infinite sum of past and future samples from $y$. It is possible to perform approximate reconstruction by truncating this sum, however the number of samples required to gives error $\epsilon$ will be large: $\Omega(1/\epsilon)$. See Example 25 at the end of the section.

Progress on the finite time reconstruction question beyond truncated Whittaker-Shannon began with the pioneering work Slepain, Landau, and Pollak, who study the operator $\mathcal{K}_{\mu}$ for uniform, bandlimited measures $\mu$ [SP61, LP61, LP62]. They bound the number of eigenvalues of $\mathcal{K}_{\mu}$ above $\epsilon$, a quantity that is at most a constant factor larger than our $s_{\mu,\epsilon}$. Using this bound, it is possible to argue that Problem 1 can be solved via regression onto at most $O(s_{\mu,\epsilon})$ prolate spheroidal wave functions (PSWFs).

While the prolate spheroidal wave functions cannot be explicitly represented and used directly in a regression algorithm, later work presents practical methods for working with them using quadrature rules and a finite number of time samples [XRY01, STR06, OR14]. For the noiseless version of Problem 1, that work, combined with the statistical dimension bound of Landau and Widom [LW80], yields algorithms that take roughly $\tilde{O}\left(FT+\log(1/\epsilon)\right)$ samples and $\tilde{O}\left(\left(FT+\log(1/\epsilon)\right)^{\omega}\right)$ time, matching our results up to log factors.242424Landau and Widom’s bounds can be used to show that $s_{\mu,\epsilon}=\tilde{\Theta}(FT+\log(1/\epsilon))$, however their result only holds asymptotically as $FT$ goes to infinity. Ours holds for all values of $F,T,\epsilon$ – see Theorem 48.

We note that existing quadrature methods access the function $f(t)$ at a pre-determined set of time domain points. Thus, they are inherently not robust to noise, since the noise function $n(t)$ of Problem 1 can place arbitrarily bad corruptions on $\mathcal{F}_{\mu}x$ at the pre-determined sample points. To the best of our knowledge our work is the first to solve Problem 1 for bandlimited signals in the adversarial noise setting.

Sparse functions. Signal interpolation has also been studied extensively when $y$ is assumed to have a sparse Fourier transform: this is the basic problem of compressed sensing and sparse recovery. While most results in the compressed sensing literature are for discrete functions and address sparsity in the discrete Fourier transform, there has been interest in extending that work to the continuous case [DW12]. Furthermore, there are a number of results on the continuous problem that predate compressed sensing: reader’s may be familiar with Prony’s method [dP95], Pisarenko’s method [Pis73], the matrix pencil method [BM86], or the MUSIC algorithm [Sch86].

While these methods do not provide direct guarantees for Problem 1, recently Chen, Kane, Price, and Song study a formulation of the sparse signal interpolation problem that closely matches our formulation [CKPS16]. Follow-up work in [CP18] achieves a sample complexity of $\tilde{O}(k\log^{2}k)$, exactly matching our bounds. In fact, our proofs for general constraint measures rely directly on two essentially lemmas on the smoothness of Fourier-sparse functions from [CKPS16] and [CP18].

We note that most compressed sensing type results, including those of [CKPS16, CP18], are distinguished from our work in that they can also learn the support of $\mu$ – our methods assume that this support is known a priori. We believe that our methods can be combined with existing techniques for learning the Fourier support and view this is an interesting open direction.

Multiband functions. Due to applications in radio, radar, medical imaging, and many other areas, there has been substantial interest in sample efficient algorithms for reconstructing multiband functions [Eld15]. Landau [Lan67a, Lan67b] was the first to characterize the sample complexity of reconstructing such functions in the sense of the Nyquist-Shannon sampling theorem, showing that to recover a signal with $s$ frequency bands of widths $F_{1},\ldots,F_{s}$, the average sampling rate must be at least $1/\sum_{i}F_{i}$. Unlike bandlimited interpolation, it is not obvious how to construct sampling schemes that achieve this optimal rate, and doing so has been the subject of a rich line of work on non-uniform “multicoset” sampling schemes [FB96, VB00, Bre08, ME09].

As in the bandlimited setting, the rate of the infinite time-horizon problem suggests, but does not imply, that the finite-time problem can be solved with roughly $\sum_{i}F_{i}T$ samples. Via a direct analysis of $\mathcal{K}_{\mu}$, results on prolate spheroidal wave functions can be used to upper bound the Fourier statistical dimension for a multiband support by $O\left(\sum_{i}F_{i}T+s\log(1/\epsilon)\right)$ [LW80], matching our bound in Theorem 53. However, we are unaware of existing work that solves Problem 1 with a number of samples matching this statistical dimension bound. We suspect that, as in the bandlimited case, in the noiseless setting, our methods could be matched via a combination of numerical quadrature and PSWF regression. There has been some initial work in that direction [LH12].

General constraints. Beyond the three standard settings discussed above, there has been an effort to understand the complexity of approximately reconstructing functions with more general Fourier transform constraints. In the discrete setting, model based compressed sensing has proven to be a powerful framework [BCDH10, HIS15]. Similar ideas have been extended to continuous functions [EM09, Eld09]. However, the constraints considered in prior work do not correspond with those captured by Problem 1. We are interested in a more refined understanding of how sample complexity depends on a the complexity of a function’s representation in the Fourier basis. Model based compressed sensing focuses on functions that can be sparsely represented in other bases.

Leverage score sampling. Finally, we note that, beyond widespread applications in randomized numerical linear algebra, there has been prior work studying leverage score sampling schemes for discretizing continuous operators, which is the approach we take to solving Problem 1. See for example [CM17, Bac17, PBV18]. Our main contribution is demonstrating how to actually upper bound the leverage scores for operators of interest, which is the missing ingredient that typically prevents such sampling results from being algorithmic. We think that the tools presented in this paper offer a powerful approach to discretization in general, with significant potential for future research. We use similar methods in our recent work on randomized approximation schemes for Gaussian kernel matrices [AKM*+*17].

Example 25**.**

Truncated Whittaker-Shannon interpolation requires $\Omega(1/\epsilon)$ samples to approximate $y(t)$ with bandlimit $F=1/2$ on $[0,1]$ up to error $\epsilon$ (i.e., to solve Problem 1, outputting $\tilde{y}(t)$ with $\|y-\tilde{y}\|_{[0,1]}^{2}\leq\epsilon\|\hat{y}\|_{\mu}^{2}$ where $\mu$ is the uniform probability measure on $[-1/2,1/2]$ and $\hat{y}$ is the Fourier transform of $y$ with $y=\mathcal{F}_{\mu}^{*}\hat{y}$.)

Proof.

Let $\mathcal{E}$ be the set of even integers in $[\lfloor 1/2\epsilon\rfloor,\lfloor 1/\epsilon\rfloor]$. Note that $|\mathcal{E}|=\Theta(1/\epsilon)$. Let $y$ be a sum of $\Theta(1/\epsilon)$ standard sinc functions centered at the points in $\mathcal{E}$:

[TABLE]

The Fourier transform $\hat{y}_{k}(\xi)$ is the box on $[-1/2,1/2]$ multiplied by $e^{-2\pi ik}$. Thus the Fourier transform $\hat{y}(\xi)=\sum_{k\in\mathcal{E}}\hat{y}_{k}(\xi)$ is supported on $[-1/2,1/2]$ and so the Nyquist rate is $1$ and Whittaker-Shannon interpolation reconstructs $y(t)$ as a sum of sinc functions centered at the integers:

[TABLE]

We can see that this reconstruction is exact since $y(k)=0$ for all integer $k$ except $y(k)=1$ for $k\in\mathcal{E}$. However, if we approximate $y(t)$ on the range $[0,1]$ by truncating the Whittaker-Shannon sum to $\leq\lfloor 1/\epsilon\rfloor$ samples centered at [math], we will not include the terms corresponding to $k\in[\lfloor 1/2\epsilon\rfloor,\lfloor 1/\epsilon\rfloor]\supseteq\mathcal{E}$ and so will have $\tilde{y}(t)=0$ and so $\|y-\tilde{y}\|_{[0,1]}^{2}=\|y\|_{[0,1]}^{2}$. Since $\mathcal{E}$ is the set of even integers in $[\lfloor 1/2\epsilon\rfloor,\lfloor 1/\epsilon\rfloor]$:

[TABLE]

Finally we note that the for $j\neq k$ $\langle\hat{y}_{k},\hat{y}_{j}\rangle_{\mu}=\int_{-1/2}^{1/2}e^{-2\pi i(j-k)\xi}d\xi=0$. Thus

[TABLE]

Combined with (A) this gives:

[TABLE]

which completes the lower bound. ∎

Appendix B Operator theory preliminaries

Throughout the paper, we use the term operator for linear transformation between two Hilbert spaces. In this section we discuss and prove basic results on operators that we use throughout the paper.

B.1 Basic definitions and the Loewner partial ordering

Consider two Hilbert spaces ${\cal H}_{1}$ and ${\cal H}_{2}$ with inner products $\langle\cdot,\cdot\rangle_{\mathcal{H}_{1}}$ and $\langle\cdot,\cdot\rangle_{\mathcal{H}_{2}}$. We denote by ${\mathbb{B}}({\cal H}_{1},{\cal H}_{2})$ the set of bounded operators from ${\cal H}_{1}$ to ${\cal H}_{2}$, and abbreviate ${\mathbb{B}}({\cal H})$ if ${\cal H}_{1}={\cal H}_{2}={\cal H}$. We denote by $\mathbb{B}_{TC}({\cal H})$ and $\mathbb{B}_{HS}({\cal H})$ the set of trace-class and Hilbert-Schmidt operators (respectively) on ${\cal H}$ (i.e. from ${\cal H}$ to ${\cal H}$). Note that $\mathbb{B}_{TC}({\cal H})\subset\mathbb{B}_{HS}({\cal H})\subset{\mathbb{B}}({\cal H})$. Recall that for operators, boundedness is equivalent to continuity. The open mapping theorem states that if ${\cal A}$ is invertible, then ${\cal A}^{-1}$ is bounded. This implies that a compact operator is not invertible. If ${\cal A}\in{\mathbb{B}}({\cal H})$ and ${\cal B}\in\mathbb{B}_{TC}({\cal H})$ then ${\cal A}{\cal B},{\cal B}{\cal A}\in\mathbb{B}_{TC}({\cal H})$ and $\operatorname{tr}({\cal A}{\cal B})=\operatorname{tr}({\cal B}{\cal A})$.

We call self-adjoint ${\cal A}$ positive semidefinite (or simply positive) and write $\mathcal{A}\succeq 0$ if $\langle x,{\cal A}x\rangle_{\mathcal{H}}\geq 0$ for all $x\in{\cal H}$. We write ${\cal A}\succsim 0$ if ${\cal A}$ is positive definite, i.e. $\langle x,{\cal A}x\rangle>0$ for all $x\in{\cal H}$. We denote ${\cal A}\succ 0$ is ${\cal A}$ if ${\cal A}$ is strictly positive, i.e. there exist a $c>0$ such that ${\cal A}\succsim c\cdot{\cal I}_{\mathcal{H}}$ where ${\cal I}_{\mathcal{H}}$ is the identity operator on $\mathcal{H}$. Note that for operators on finite dimensional Hilbert spaces, ${\cal A}\succsim 0$ if and only if ${\cal A}\succ 0$, but this is not always the case for infinite dimensional Hilbert spaces. The notation for ${\cal A}\succeq{\cal B}$, ${\cal A}\succsim{\cal B}$, and ${\cal A}\succ{\cal B}$ follow in the standard way.

If ${\cal A}\succeq 0$ is self-adjoint and bounded, then it possesses a unique self-adjoint bounded square root ${\cal A}^{1/2}\succeq 0$ [Wou66]. Furthermore, if ${\cal A}$ is strictly positive then so is ${\cal A}^{1/2}$. This implies that if ${\cal A}$ is strictly positive and bounded, then ${\cal A}$ is bounded below and that the inverse of of the square root of ${\cal A}$ is ${\cal A}^{-1/2}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}({\cal A}^{-1})^{1/2}$. Lidskii’s theorem states the that trace of a trace-class operator is the sum of its eigenvalues.

Many of the following claims are well known of matrices, and the proofs in most cases, but not all, mirror the matrix case. However, for the operator case we need to be more careful with the conditions due to the aforementioned distinction between $\succsim$ and $\succ$.

Claim 26**.**

Suppose that ${\cal A}$ is a self-adjoint bounded positive semidefinite operator on an Hilbert space ${\cal H}$. For every $\epsilon>0$, the operator ${\cal A}+\epsilon{\cal I}_{\cal H}$ is bounded, strictly positive and invertible, and the inverse is bounded.

Proof.

The operator ${\cal A}+\epsilon{\cal I}_{\cal H}$ is the sum of two bounded operators, and so it is bounded. It is also clearly bounded below, since ${\cal A}+\epsilon{\cal I}_{\cal H}\succeq\epsilon{\cal I}_{\cal H}\succ 0$. A continuous (i.e., bounded) bounded-below operator is invertible, so ${\cal A}+\epsilon{\cal I}_{\cal H}$ is invertible. The inverse is bounded due to the open mapping theorem. ∎

Claim 27**.**

Suppose that $0\prec{\cal A}\preceq{\cal I}_{\cal H}$ for a self-adjoint operator ${\cal A}$. Then, ${\cal A}^{-1}\succeq{\cal I}_{\cal H}$.

Proof.

For every $x\in{\cal H}$ we have $\langle x,{\cal A}x\rangle_{\mathcal{H}}\leq\langle x,x\rangle_{\mathcal{H}}$. Given $y$, let $x={\cal A}^{-1/2}y$. Then $\langle y,y\rangle_{\mathcal{H}}=\langle{\cal A}^{1/2}x,{\cal A}^{1/2}x\rangle_{\mathcal{H}}=\langle x,{\cal A}x\rangle_{\mathcal{H}}\leq\langle x,x\rangle_{\mathcal{H}}=\langle{\cal A}^{-1/2}y,{\cal A}^{-1/2}y\rangle_{\mathcal{H}}=\langle y,{\cal A}^{-1}y\rangle_{\mathcal{H}}$ so ${\cal A}^{-1}\succeq{\cal I}_{\cal H}$. ∎

Claim 28**.**

Suppose that ${\cal A}\in{\mathbb{B}}({\cal H})$ and that ${\cal B}\succeq 0$ is self-adjoint trace-class operator. Then, ${\cal B}^{1/2}{\cal A}{\cal B}^{1/2}$ is trace-class, and $\operatorname{tr}({\cal B}^{1/2}{\cal A}{\cal B}^{1/2})=\operatorname{tr}({\cal A}{\cal B})$.

Proof.

Since ${\cal B}$ is trace-class, ${\cal B}^{1/2}\in\mathbb{B}_{HS}({\cal H})$. This implies that ${\cal A}{\cal B}^{1/2}$ is also Hilbert-Schmidt. Thus, ${\cal B}^{1/2}{\cal A}{\cal B}^{1/2}$ is the product of two Hilbert-Schmidt operators, so it is trace-class. The trace equality follows from the cyclic property of the trace. ∎

Claim 29**.**

Suppose that ${\cal A}\succ 0$ is a self-adjoint bounded operator, and that ${\cal B}\succeq 0$ is self-adjoint trace-class operator, both on a separable Hilbert space ${\cal H}$. Suppose we have $\operatorname{tr}({\cal A}{\cal B})\leq 1$. Then, ${\cal B}\preceq{\cal A}^{-1}$.

Proof.

Due to the cyclicity of the trace $\operatorname{tr}({\cal A}^{1/2}{\cal B}{\cal A}^{1/2})\leq 1$. The operator ${\cal A}^{1/2}{\cal B}{\cal A}^{1/2}$ is positive semidefinite, so due to Lidskii’s theorem it’s largest eigenvalue $\leq 1$. For ${\cal A}^{1/2}{\cal B}{\cal A}^{1/2}$, the largest eigenvalue is equal to the operator norm , so for any $y$,

[TABLE]

Since ${\cal A}^{1/2}$ is invertible, with inverse ${\cal A}^{-1/2}$, the conclusion of the claim follows. ∎

Claim 30**.**

Let ${\cal A},{\cal B}$ be two self-adjoint, bounded, strictly positive operators. If ${\cal A}\preceq{\cal B}$ then ${\cal A}^{-1}\succeq{\cal B}^{-1}$.

Proof.

Since ${\cal B}$ is bounded and strictly positive, then it is invertible and has an invertible square root. For any $y\in{\cal H}$ let $x={\cal B}^{-1/2}y$. We have

[TABLE]

So ${\cal B}^{-1/2}{\cal A}{\cal B}^{-1/2}\preceq{\cal I}_{\cal H}$. Since both ${\cal A}$ and ${\cal B}$ are strictly positive, then ${\cal B}^{-1/2}{\cal A}{\cal B}^{-1/2}$ is also strictly positive. Thus, according to Claim 27, ${\cal B}^{1/2}{\cal A}^{-1}{\cal B}^{1/2}\succeq{\cal I}_{\cal H}$, from which the claim easily follows. ∎

Claim 31**.**

Suppose that ${\cal A}\succ 0$ and ${\cal A}\succeq{\cal B}$. Then for any $0\leq c<1$ we have ${\cal A}-c{\cal B}\succ 0$.

Proof.

Suppose by contradiction that ${\cal A}-c{\cal B}\not\succ 0$. Then for any $\epsilon>0$ there exists a $x$ with unit norm ($\langle x,x\rangle_{\mathcal{H}}=1$) such that $\langle x,({\cal A}-c{\cal B})x\rangle_{\mathcal{H}}\leq\epsilon$. We have $\langle x,{\cal B}x\rangle_{\mathcal{H}}\geq(\langle x,{\cal A}x\rangle_{\mathcal{H}}-\epsilon)/c$, and since $\langle x,{\cal A}x\rangle_{\mathcal{H}}$ is bounded away from zero and $c<1$, for sufficiently small $\epsilon$ we have $\langle x,{\cal B}x\rangle_{\mathcal{H}}>\langle x,{\cal A}x\rangle_{\mathcal{H}}$ so $\langle x,({\cal A}-{\cal B})x\rangle_{\mathcal{H}}<0$ which contradicts the assumption that ${\cal A}\succeq{\cal B}$. ∎

Definition 5**.**

Given $x\in{\cal H}_{1}$ and $y\in{\cal H}_{2}$, we define the operator $x\otimes y:{\cal H}_{2}\to{\cal H}_{1}$ by

[TABLE]

Claim 32**.**

Let ${\cal H}$ be a separable Hilbert space, and assume that ${\cal A}\in\mathbb{B}({\cal H})$ and $v\in{\cal H}$. Then, $\langle v,{\cal A}v\rangle_{\cal H}=\operatorname{tr}({\cal A}(v\otimes v))$. (We remark that ${\cal A}(v\otimes v)$ is trace-class since $v\otimes v$ has finite-rank and ${\cal A}$ is bounded.)

Proof.

Let $e_{1},e_{2},\ldots$ be an orthonormal basis for ${\cal H}$. Write $v=\sum^{\infty}_{i=1}\alpha_{i}e_{i}$. On one hand we have

[TABLE]

On the other we have

[TABLE]

so the two terms are equal. ∎

B.2 Weak integrals of operators

We are going to work with operator-valued random variables. To reason about the expected value, we need a notion of an integral of operator-valued functions. We use a generalization of the concept of weak integrals (also called Pettis integral) of vector-valued functions [Pet38].

Definition 6**.**

Let ${\cal H}_{1},{\cal H}_{2}$ be two separable Hilbert spaces, $G$ a measurable space and $\mu$ a measure on $G$, and consider a mapping $\mathcal{A}:G\to\mathbb{B}({\cal H}_{1},{\cal H}_{2})$. If the mapping $(x,z)\mapsto\int_{G}\langle x,{\cal A}(\xi)z\rangle_{{\cal H}_{2}}\,d\mu(\xi)$ is a bounded sesquilinear map in $x,z$, then we say that ${\cal A}$ is a weakly integrable operator valued function and the weak operator integral is defined to be the unique bounded operator

[TABLE]

such that for all $x$ and $z$ we have

[TABLE]

The existence and uniqueness of such an operator is guaranteed by the Riesz representation theorem for sesquilinear maps [Hel69, Page 92, Theorem 5].252525We remark that [Hel69, Page 92, Theorem 5] is stated and proved only for sesquilinear forms on the same Hilbert space (i.e., ${\cal H}_{1}={\cal H}_{2}$). However, it is easy to verify that the result also holds for sesquilinear forms between two Hilbert spaces.

Claim 33**.**

Suppose that ${\cal A}:G\to\mathbb{B}({\cal H}_{1},{\cal H}_{2})$ is weakly integrable operator valued function, and ${\cal S}\in\mathbb{B}({\cal H}_{2}),{\cal T}\in\mathbb{B}({\cal H}_{1})$. Then $\xi\mapsto{\cal T}{\cal A}(\xi){\cal S}$ is also a weakly integrable operator valued function and

[TABLE]

Proof.

Recall that $(x,z)\mapsto\int_{G}\langle x,{\cal A}(\xi)z\rangle_{{\cal H}_{2}}\,d\mu(\xi)$ is bounded, so there exists a $\gamma$ such that for every $x\in{\cal H}_{2},z\in{\cal H}_{1}$ we have

[TABLE]

Let $x\in{\cal H}_{2},z\in{\cal H}_{1}$. We have

[TABLE]

where we used the fact that both ${\cal S}$ and ${\cal T}$ are bounded. So the mapping $(x,z)\mapsto\int_{G}\langle x,{\cal T}{\cal A}(\xi){\cal S}z\rangle_{{\cal H}_{2}}d\mu(\xi)$ is bounded and $\xi\mapsto{\cal T}{\cal A}(\xi){\cal S}$ is weakly integrable.

We now show that the value of the integral is ${\cal T}\left(\int_{G}{\cal A}(\xi)d\mu(\xi)\right){\cal S}$. Again, for any $x\in{\cal H}_{2},z\in{\cal H}_{1}$:

[TABLE]

By definition of $\int_{G}{\cal A}(\xi)d\mu(\xi)$ we have

[TABLE]

so indeed $\int_{G}{\cal T}{\cal A}(\xi){\cal S}d\mu(\xi)={\cal T}\left(\int_{G}{\cal A}(\xi)d\mu(\xi)\right){\cal S}$. ∎

Claim 34**.**

Let $\rho,\mu$ be two, possibly different, probability measures, on $\mathbb{R}$, and let ${\cal A}\in{\mathbb{B}}(L_{2}(\rho))$ be self-adjoint and positive semi-definite, and let ${\cal B}\in\mathbb{B}_{TC}(L_{2}(\rho))$. Assume that there exists an orthonormal basis for $L_{2}(\rho)$ consisting of eigenvectors of ${\cal A}$. Given a mapping $\eta\in\mathbb{R}\mapsto v_{\eta}\in L_{2}(\rho)$ such that ${\cal B}=\int_{\mathbb{R}}(v_{\eta}\otimes v_{\eta})d\mu(\eta)$ we have:

[TABLE]

Proof.

Let $e_{1},e_{2},\ldots$ be an orthonormal basis for $L_{2}(\rho)$ consisting of eigenvectors of ${\cal A}$. Using Claim 32, we have

[TABLE]

where the exchange of the integral and infinite sum in the third equality is justified by Tonelli’s Theorem. In order to apply Tonelli’s theorem we need to show that $\langle e_{i},{\cal A}(v_{\eta}\otimes v_{\eta})e_{i}\rangle_{\rho}\geq 0$ for every $i$ and $\eta$. This is indeed the case since $\langle e_{i},{\cal A}(v_{\eta}\otimes v_{\eta})e_{i}\rangle_{\rho}=\langle{\cal A}e_{i},(v_{\eta}\otimes v_{\eta})e_{i}\rangle_{\rho}=\lambda_{i}\langle e_{i},(v_{\eta}\otimes v_{\eta})e_{i}\rangle_{\rho}\geq 0$ where $\lambda_{i}$ is the eigenvalue corresponding to $e_{i}$. Note that since ${\cal A}$ is self-adjoint and positive semi-definite, $\lambda_{i}$ is real and non-negative. We also used the immediate fact that $v_{\eta}\otimes v_{\eta}$ is positive semi-definite. ∎

Remark: One way to guarantee that there is an orthonormal basis of eigenvectors of ${\cal A}$ is to require ${\cal A}$ to be compact. However, it is quite possible for ${\cal A}$ not to be compact, and still have an orthonormal basis of eigenvectors. In fact, we primarily apply Claim 34 to operators of the form $({\cal C}+\epsilon{\cal I})^{-1}$ where ${\cal C}$ is compact, and such operators have an orthonormal basis of eigenvectors (since they share eigenvectors with ${\cal C}$).

We say that a weakly integrable ${\cal A}(\cdot)$ is self-adjoint if ${\cal A}(\xi)$ is self-adjoint for all $\xi$. It is easy to verify that if ${\cal A}(\cdot)$ is self-adjoint, then $\int_{G}A(\xi)d\mu(\xi)$ is self-adjoint as well.

Claim 35**.**

Suppose that ${\cal A},{\cal B}:G\to\mathbb{B}({\cal H})$ are two self-adjoint weakly integrable operator valued functions. If, with respect to a measure $\mu$ on $G$, ${\cal A}(\xi)\preceq{\cal B}(\xi)$ almost everywhere, then $\int_{G}{\cal A}(\xi)d\mu(\xi)\preceq\int_{G}{\cal B}(\xi)d\mu(\xi)$.

Proof.

For every $x\in{\cal H}$,

[TABLE]

so indeed $\int_{G}{\cal A}(\xi)d\mu(\xi)\preceq\int_{G}{\cal B}(\xi)d\mu(\xi)$. ∎

Claim 36**.**

Suppose that ${\cal B}:G\to\mathbb{B}({\cal H})$ is a self-adjoint weakly integrable operator valued function. Consider another self-adjoint operator valued function ${\cal A}:G\to\mathbb{B}({\cal H})$. If for every $\xi\in G$ we have $0\preceq{\cal A}(\xi)\preceq{\cal B}(\xi)$, then ${\cal A}$ is weakly integrable and $\int_{G}{\cal A}(\xi)d\mu(\xi)\preceq\int_{G}{\cal B}(\xi)d\mu(\xi)$.

Proof.

We need to prove only that ${\cal A}$ is weakly integrable, since the integral bound follows from Claim 35. A sesquilinear form is bounded if and only if the associated quadratic form is bounded [Hel69, Page 92, Theorem 3], so we need to show that the integral of the quadratic form associated with ${\cal A}$ is bounded. Since ${\cal A}(\xi)$ is always positive semidefinite, for any $x$

[TABLE]

and since the integral of the quadratic form associated with ${\cal B}$ is bounded (since ${\cal B}$ is weakly integrable) we conclude that integral quadratic form associated with ${\cal A}$ is bounded, so indeed ${\cal A}$ is weakly integrable. ∎

B.3 Concentration of random operators

Let ${\cal A}:G\to\mathbb{B}({\cal H})$ be a weakly integrable operator valued function. If the underlying measure $\mu$ is a probability measure, then we shall call ${\cal A}$ a random operator, and write

[TABLE]

Certain matrix concentration results can be generalized to the case that ${\cal A}$ is a random operator which takes only self-adjoint Hilbert-Schmidt values. The underlying reason is that Hilbert-Schmidt operators can be well-approximated using finite rank operators. The basic technique is outlined in [Min17, Section 3.2]. We use this technique to prove the following lemma.

Lemma 37**.**

Suppose that ${\cal H}$ is a separable Hilbert space, and that ${\cal B}$ is a fixed self-adjoint Hilbert-Schmidt operator on ${\cal H}$. Let ${\cal R}$ be a self-adjoint Hilber-Schmidt random operator that satisfies

[TABLE]

Let ${\cal M}$ be another self-adjoint trace-class operator such that $\mathbb{E}({\cal R}^{2})\preceq{\cal M}$. Form the operator sampling estimator

[TABLE]

where each ${\cal R}_{k}$ is an independent copy of ${\cal R}$. Then, for all $t>\sqrt{\|{\cal M}\|_{\mathrm{op}}/n}+2L/3n$,

[TABLE]

Proof.

Let $e_{1},e_{2},\dots$ be the eigenvectors of ${\cal M}$, ordered according to the magnitude of the corresponding eigenvalue, and let ${\cal P}_{j}$ be the orthogonal projector on the span of $e_{1},e_{2},\dots,e_{j}$. Consider the finite-rank operators ${\cal R}^{(j)}={\cal P}_{j}{\cal R}{\cal P}_{j}$, ${\cal R}^{(j)}_{k}={\cal P}_{j}{\cal R}_{k}{\cal P}_{j}$, $\bar{{\cal R}}^{(j)}_{n}={\cal P}_{j}\bar{{\cal R}}_{n}{\cal P}_{j}$, ${\cal B}^{(j)}={\cal P}_{j}{\cal B}{\cal P}_{j}$ and ${\cal M}^{(j)}={\cal P}_{j}{\cal M}{\cal P}_{j}$. We will apply on these operator sequences the matrix version of the current lemma [AKM*+*17]262626The lemma in [AKM*+*17] is stated as a bound on $\Pr\left(\|\bar{{\cal R}}_{n}-{\cal B}\|_{\mathrm{op}}\geq t\right)$, while for operators strict inequality is necessary. It is easy to verify that the matrix version of the Lemma continues to hold for $\Pr\left(\|\bar{{\cal R}}_{n}-{\cal B}\|_{\mathrm{op}}>t\right)$.

Due to linearity of weak operator integrals we have $\mathbb{E}({\cal R}^{(j)})={\cal P}_{j}{\cal B}^{(j)}{\cal P}_{j}$. We can bound the operator norm of ${\cal R}^{(j)}$: $\|{\cal R}^{(j)}\|_{\mathrm{op}}\leq\|{\cal P}_{j}{\cal R}{\cal P}_{j}\|_{\mathrm{op}}\leq\|{\cal P}_{j}\|_{\mathrm{op}}^{2}\|{\cal R}\|_{\mathrm{op}}\leq L$ since the operator norm of a projection operator is $1$. Using the fact that ${\cal P}_{j}\preceq{\cal I}_{\cal H}$ and so ${\cal R}{\cal P}_{j}{\cal R}\preceq{\cal R}^{2}$ we have

[TABLE]

Now applying the aforementioned matrix version of the current lemma272727Technically, the aforementioned concentration result is for matrices, while here we deal with abstract operators on finite dimensional subspaces. We can address this issue by using the corresponding transformation matrices, but we find that to be tedious details. we find that

[TABLE]

Due to the way we constructed ${\cal P}_{j}$, and ${\cal M}$ being trace-class, we have $\operatorname{tr}({\cal M}^{(j)})\to\operatorname{tr}({\cal M})$ as $j\to\infty$. Furthermore, since ${\cal M}$ is trace-class, ${\cal P}_{j}{\cal M}\to{\cal M}$ uniformly [HN01, Theorem 9.21], and so ${\cal M}^{(j)}\to{\cal M}$ while implies that $\|{\cal M}^{(j)}\|_{\mathrm{op}}\to\|{\cal M}\|_{\mathrm{op}}$. Thus, the entire right side of (43) converges to the right side of (42), so

[TABLE]

Let $G$ and $\mu$ denote the underlying probability space and probability measure. Let $f_{j}$ now denote the indicator function for the event $\|\bar{{\cal R}}^{(j)}_{n}-{\cal B}^{(j)}\|_{\mathrm{op}}>t$, and $f$ the indicator for the event $\|\bar{{\cal R}}_{n}-{{\cal B}}\|_{\mathrm{op}}>t$. Again, due to the fact that $\bar{{\cal R}}_{n}-{{\cal B}}$ is Hilbert-Schmidt we have $\bar{{\cal R}}^{(j)}_{n}-{\cal B}^{(j)}\to\bar{{\cal R}}_{n}-{{\cal B}}$, while implies that that for any $\xi\in G$, $f(\xi)=\liminf_{j\to\infty}f_{j}(\xi)$. Now due to Fatou’s lemma:

[TABLE]

∎

Appendix C Properties of the ridge leverage scores

C.1 Basic facts about leverage scores

In this section we prove Theorem 5, giving fundamental properties of the ridge leverage scores that we use both in bounding these scores and proving that ridge leverage score sampling can be used to solve the regularized regression problem of (10) (and hence Problem 1 by Claim 4).

Theorem 5 (Leverage Function Properties).

Letting $\tau_{\mu,\epsilon}(t)$ be the ridge leverage function of Definition 3, that is

[TABLE]

and let $\varphi_{t}\in L_{2}(\mu)$ be defined by $\varphi_{t}(\xi)=e^{-2\pi it\xi}$, we have the following basic properties:

• •

The leverage scores integrate to the statistical dimension:

[TABLE]

• •

Inner Product characterization:

[TABLE]

• •

Minimization Characterization:

[TABLE]

Proof.

Recall, that given $t\in[0,T]$, we defined $\varphi_{t}(\xi)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}e^{-2\pi it\xi}$ ($\varphi_{t}\in L_{2}(\mu)$). It is easy to verify that:

[TABLE]

To prove the equality between Equations (44), (46), and (47), we first show that the right hand side of (46) is equal to the right hand side of (47) and then show that the right hand side of (46) is equal to the right hand side of (44).

First, we need an auxiliary lemma regarding the solution of regularized least squares problems. If $\cal A$ is matrix with full column rank or a one-to-one linear operator between finite-dimensional Hilbert spaces, and $b$ some vector, then $F(x)=\|{\cal A}x-b\|^{2}$ has a unique minimizer. In infinite dimension spaces, this remains true if only the co-domain of $\cal A$ is infinite dimensional. However, if both the domain and co-domain are infinite dimensional there might not be a minimizer even if the $\cal A$ is bounded: the range of $\cal A$ might not be closed, so it is possible that $\|{\cal A}x-b\|>0$ for every $x$, but also that there exists a series $\{x_{n}\}$ such that $\|\mathcal{A}x_{n}-b\|\to 0$ as $n\to\infty$. However, once we introduce a ridge term (i.e., minimize $F(x)=\|\mathcal{A}x-b\|^{2}+\lambda\|x\|^{2}$ for some $\lambda>0$) there is always a unique minimizer (as long as $\mathcal{A}$ is bounded), due to the extreme value theorem (since we can bound the search domain). Furthermore, we can write an analytic expression for the minimizer in an analogous way to the finite dimensional case, as the following lemma shows.

Lemma 38** (Regularized Least Squares Minimizer).**

Let $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$ be two Hilbert spaces, and $\mathcal{A}:\mathcal{H}_{1}\to\mathcal{H}_{2}$ be a bounded linear operator. Let $b\in\mathcal{H}_{2}$ and $\lambda>0$. The function

[TABLE]

has a unique minimizer which is $x^{\star}=\mathcal{A}^{*}(\mathcal{A}\mathcal{A}^{*}+\lambda\mathcal{I}_{\mathcal{H}_{2}})^{-1}b$.

Proof.

Consider the Hilbert space $\mathcal{H}_{1}\times\mathcal{H}_{2}$ equipped with the inner product

[TABLE]

Define the operator $\mathcal{T}:\mathcal{H}_{1}\to\mathcal{H}_{1}\times\mathcal{H}_{2},\quad\mathcal{T}(x)=(\sqrt{\lambda}x,\mathcal{A}x)$. Let $y=(0,b)\in\mathcal{H}_{1}\times\mathcal{H}_{2}$. We have $F(x)=\|\mathcal{T}x-y\|^{2}_{\mathcal{H}_{1}\times\mathcal{H}_{2}}$. Thus, we need to show that there is a unique point $\tilde{y}\in\operatorname{range}(\mathcal{T})$ that minimizes $\|\tilde{y}-y\|^{2}_{\mathcal{H}_{1}\times\mathcal{H}_{2}}$ and that $\tilde{y}=\mathcal{T}x^{\star}$ for $x^{\star}=\mathcal{A}^{*}(\mathcal{A}\mathcal{A}^{*}+\lambda\mathcal{I}_{\mathcal{H}_{2}})^{-1}b$.

The operator $\mathcal{T}$ is a bounded linear operator, so it is continuous. We also have that for every $x\in\mathcal{H}_{1}$, $\|\mathcal{T}x\|^{2}_{\mathcal{H}_{1}\times\mathcal{H}_{2}}\geq\lambda\|x\|^{2}_{\mathcal{H}_{1}}$ where $\lambda>0$, so $\mathcal{T}$ is bounded from below. So $\mathcal{T}$ has a closed range [AA02, Theorem 2.5]. Thus, there is a unique $\tilde{y}\in\operatorname{range}(\mathcal{T})$ that minimizes $\|\tilde{y}-y\|^{2}_{\mathcal{H}_{1}\times\mathcal{H}_{2}}$, and that $\tilde{y}$ is the unique element of $\operatorname{range}(\mathcal{T})$ with the property $y-\tilde{y}\perp\operatorname{range}(\mathcal{T})$ [HN01, Theorem 6.13]. So it suffices to show that for every $x\in\mathcal{H}_{1}$ we have $y-\mathcal{T}x^{\star}\perp\mathcal{T}x$. We compute:

[TABLE]

So indeed, for every $x\in\mathcal{H}_{1}$ we have $y-\mathcal{T}x^{\star}\perp\mathcal{T}x$ and $x^{\star}$ is the unique minimizer. ∎

Using Lemma 38 we now proceed with the proof of Theorem 5.

Corollary 39**.**

Let

[TABLE]

Then,

[TABLE]

Claim 40**.**

We have

[TABLE]

so the right hand side of (46) is equal to the right hand side of (47).

Proof.

We compute:

[TABLE]

and

[TABLE]

Summing the last equalities completes the proof. ∎

Claim 41**.**

We have

[TABLE]

so the right hand side of (46) is equal to the right hand side of (44).

Proof.

We can reformulate the previous claim as :

[TABLE]

Let the optimal solution be $\beta^{\star}$ and $u^{\star}$. We have $\varphi_{t}=\mathcal{F}_{\mu}\beta^{\star}+\sqrt{\epsilon}u^{\star}$, hence for any $0\neq\alpha\in L_{2}(\mu)$:

[TABLE]

where the last inequality follows from Cauchy-Schwarz inequality. Using Cauchy-Schwarz again:

[TABLE]

So for every $0\neq\alpha\in L_{2}(\mu)$:

[TABLE]

We conclude by showing that the maximum value is attained. Let $\alpha^{\star}=(\mathcal{G}_{\mu}+\epsilon\mathcal{I}_{\epsilon})^{-1}\varphi_{t}$. We have

[TABLE]

and finally,

[TABLE]

∎

We now turn to showing that the leverage function integrates to the statistical dimension.

Claim 42**.**

[TABLE]

Proof.

It follows from Eq. (48) and Claim 34 that $\int_{0}^{T}\tau_{\mu,\epsilon}(t)=\operatorname{tr}((\mathcal{G}_{\mu}+\epsilon\mathcal{I}_{\mu})^{-1}\mathcal{G}_{\mu})$. The claim follows by noting that ${\cal K}_{\mu}$ and ${\cal G}_{\mu}$ have the same eigenvalues (both operators are compact self adjoint operators, so their spectrum consists of only eigenvalues, and it is easy to verify that if $x$ is an eigenvector $\mathcal{K}_{\mu}$ then $\mathcal{F}_{\mu}x$ is eigenvector of $\mathcal{G}_{\mu}$).

∎

We thus have completed the proof of Theorem 5. ∎

C.2 Operator Approximation via Leverage Score Sampling

Analogs of the following concentration result are well known for matrices. Accordingly, the proof is an adaptation of standard proofs for finite matrix approximation by leverage score sampling, where matrix concentration results are replaced with operator concentration results. A similar strategy was employed in [Bac17].

Lemma 43**.**

Consider the conditions of Theorem 6, and denote $\widehat{\cal G}_{\mu}=\mathbf{F}\mathbf{F}^{*}$. Let $\Delta\leq 1/2$ and $\epsilon\leq\|\mathcal{G}_{\mu}\|_{\mathrm{op}}$. If $s\geq\frac{8}{3}\Delta^{-2}\tilde{s}_{\mu,\epsilon}\ln(16\tilde{s}_{\mu,\epsilon}^{2}/\delta)$, then

[TABLE]

with probability of at least $1-\delta$.

Proof.

The condition (49) is equivalent to

[TABLE]

By composing with $({\cal G}_{\mu}+\epsilon\mathcal{I}_{\mu})^{-1/2}$ on the left and right, we find that the condition is equivalent to

[TABLE]

Noticing that

[TABLE]

and that

[TABLE]

we understand that $\widehat{{\cal G}}_{\mu}=\sum^{s}_{j=1}w^{2}_{j}(\varphi_{t_{j}}\otimes\varphi_{t_{j}})$. Let

[TABLE]

Note that ${\cal X}_{j}$ is self-adjoint and Hilbert-Schmidt (since it has finite rank). We have

[TABLE]

Since the time samples are drawn randomly, ${\cal X}_{1},\dots,{\cal X}_{s}$ are random operators. We also have, using Claim 33,

[TABLE]

Write $w(t)=\sqrt{\frac{\tilde{s}_{\mu,\epsilon}}{T\cdot\tilde{\tau}_{\mu,\epsilon}(t)}}$. For every $x,z\in L_{2}(\mu)$,

[TABLE]

which shows that

[TABLE]

so,

[TABLE]

Next, we bound the operator norm of ${\cal X}_{j}$. The random operator only takes values that are both positive semidefinite and rank one, so the operator norm of ${\cal X}_{j}$ is equal to the trace of the operator. Thus, we have

[TABLE]

where the last line follows since $\tilde{\tau}_{\mu,\epsilon}(t_{j})\geq\tau_{\mu,\epsilon}(t_{j})$ by assumption.

The final ingredient for applying Lemma 37 is to bound ${\cal X}^{2}_{j}$. Again, using the fact that $\tilde{\tau}_{\mu,\epsilon}(t_{j})\geq\tau_{\mu,\epsilon}(t_{j})$ we have:

[TABLE]

So, using Claim 36,

[TABLE]

Noticing that $\operatorname{tr}({\cal M})=\tilde{s}_{\mu,\epsilon}\cdot s_{\mu,\epsilon}$ and $\|{\cal M}\|_{\mathrm{op}}=\frac{\|\mathcal{G}_{\mu}\|_{\mathrm{op}}}{\|\mathcal{G}_{\mu}\|_{\mathrm{op}}+\epsilon}\geq 1/2$ by our assumption that $\epsilon\leq\|\mathcal{G}_{\mu}\|_{\mathrm{op}}$, and Lemma 37 we have:

[TABLE]

∎

C.3 Approximate Discretization via Leverage Score Sampling

With the operator approximation bound of Lemma 43 in place, we can prove Theorem 6, which shows that we can approximately solve the regression problem of (10) (and by Claim 4 solve Problem 1) by sampling time domain points via over-approximations to their ridge leverage scores.

Theorem 6 (Approximate Regression via Leverage Score Sampling).

Assume $\epsilon\leq\|\mathcal{K}_{\mu}\|_{\mathrm{op}}$ and consider a measurable $\tilde{\tau}_{\mu,\epsilon}(t)$ with $\tilde{\tau}_{\mu,\epsilon}(t)\geq\tau_{\mu,\epsilon}(t)$ for all $t$ and let $\tilde{s}_{\mu,\epsilon}=\int_{0}^{T}\tilde{\tau}_{\mu,\epsilon}(t)dt$. Let $s=c\cdot\left(\tilde{s}_{\mu,\epsilon}\cdot[\log(\tilde{s}_{\mu,\epsilon})+1/\delta]\right)$ for sufficiently large fixed constant $c$ and let $t_{1},\ldots,t_{s}$ be time points selected by drawing each randomly from $[0,T]$ with probability proportional to $\tilde{\tau}_{\mu,\epsilon}(t)$. For $j\in 1,\ldots,s$ let $w_{j}=\sqrt{\frac{1}{sT}\cdot\frac{\tilde{s}_{\mu,\epsilon}}{\tilde{\tau}_{\mu,\epsilon}(t_{j})}}$. Let $\mathbf{F}:\mathbb{C}^{s}\rightarrow L_{2}(\mu)$ be the operator defined by:

[TABLE]

and $\mathbf{y},\mathbf{n}\in\mathbb{R}^{s}$ be the vector with $\mathbf{y}(j)=w_{j}\cdot y(t_{j})$ and $\mathbf{n}(j)=w_{j}\cdot n(t_{j})$. With probability $\geq 1-\delta$:

• •

For any $\beta\geq 0$, if $\tilde{g}\in L_{2}(\mu)$ satisfies 282828We can see that the adjoint $\mathbf{F}^{*}:L_{2}(\mu)\rightarrow\mathbb{C}^{s}$ is given by $\left[\mathbf{F}^{*}\,g\right](j)=w_{j}\cdot\int_{\mathbb{R}}g(\xi)e^{2\pi i\xi t_{j}}\,d\mu(\xi)$.

[TABLE]

then

[TABLE]

So $\tilde{g}$ provides an approximate solution to (10) and by Claim 4, $\tilde{y}=\mathcal{F}_{\mu}^{*}\tilde{g}$ solves Problem 1 with parameter $\Theta(\epsilon)$.

Proof.

Throughout the proof we will let $\bar{y}=y+n$ and $\mathbf{\bar{y}}=\mathbf{y}+\mathbf{n}$. Let

[TABLE]

By Lemma 38, $g^{\star}=\mathcal{F}_{\mu}(\mathcal{K}_{\mu}+\lambda\mathcal{I}_{T})^{-1}\bar{y}$. Denote the optimal error as $b^{\star}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\mathcal{F}_{\mu}^{*}g^{\star}-\bar{y}$ and the optimal cost as $B^{\star}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\|\mathcal{F}_{\mu}^{*}g^{\star}-\bar{y}\|_{T}^{2}+\epsilon\|g^{\star}\|_{\mu}^{2}$.

Reduction to Affine Embedding

We prove that, for all $g\in L_{2}(\mu)$, ridge leverage score sampling lets us approximate the value of the objective function of (10) when evaluated at $g$. In the randomized linear algebra literature, this is known as an affine embedding guarantee. Specifically, we show that, with probability $\geq 1-\delta$, for all $g\in L_{2}(\mu)$,

[TABLE]

where $\alpha$ is some fixed value independent of $g$ (but which depends on $\mathbf{F}$ and $\mathbf{\bar{y}}$) with $|\alpha|\leq\frac{1}{\delta}\cdot B^{\star}$.

It is not hard to see that (55) gives the theorem. For any $\tilde{g}\in L_{2}(\mu)$ satisfying:

[TABLE]

we can apply (55) to give the main claim of the theorem:

[TABLE]

Thus, we focus our attention to proving that the affine embedding guarantee of (55) holds with probability $\geq 1-\delta$.

Expression of Error in Terms of $g-g^{\star}$

We begin by showing how, for any $g\in L_{2}(\mu)$, the cost $\|\mathcal{F}_{\mu}^{*}g-\bar{y}\|_{T}^{2}+\epsilon\|g\|_{\mu}^{2}$ can be written as a function of the deviation from the optimum: $g-g^{\star}$.

Claim 44** (Expression for Excess Cost).**

For any $g\in L_{2}(\mu)$:

[TABLE]

recalling that $B^{\star}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\|\mathcal{F}_{\mu}^{*}g^{\star}-\bar{y}\|_{T}^{2}+\epsilon\|g^{\star}\|_{\mu}^{2}$ is the optimum cost of the ridge regression problem.

Proof.

Following Lemma 38 we define $\mathcal{T}:L_{2}(\mu)\to L_{2}(\mu)\times L_{2}(T),\quad\mathcal{T}g=(\sqrt{\epsilon}g,\mathcal{F}_{\mu}^{*}g)$. For any $g$, $\|\mathcal{F}_{\mu}^{*}g-\bar{y}\|_{T}^{2}+\epsilon\|g\|_{\mu}^{2}=\|\mathcal{T}g-(0,\bar{y})\|_{L_{2}(\mu)\times L_{2}(T)}^{2}$. Again, as in Lemma 38 we know $g^{\star}$ is the unique minimizer of this function with the property that $(0,y)-\mathcal{T}g^{\star}\perp\operatorname{range}(\mathcal{T})$ [HN01, Theorem 6.13]. We can thus decompose:

[TABLE]

which gives the claim. ∎

Bounding The Sampling Error

We now show that Claim 44 holds approximately, even after sampling. This almost immediately yields the affine embedding bound of (55).

Let $\tilde{B}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\|\mathbf{F}^{*}g^{\star}-\mathbf{\bar{y}}\|_{2}^{2}+\epsilon\|g^{\star}\|_{\mu}^{2}$ be the error of the optimal solution in our randomly discretized regression problem. We can write the discretized objective function value for any $g\in L_{2}(\mu)$ as:

[TABLE]

Let $\mathcal{\bar{F}}_{\mu}:L_{2}(T)\times L_{2}(\mu)\rightarrow L_{2}(\mu)$ be the operator $\mathcal{\bar{F}}_{\mu}(f,g)=\mathcal{F}_{\mu}f+\sqrt{\epsilon}\cdot g$. We can see that $\mathcal{\bar{F}}_{\mu}^{*}:L_{2}(\mu)\rightarrow L_{2}(T)\times L_{2}(\mu)$ is given by $\mathcal{\bar{F}}_{\mu}^{*}g=(\mathcal{F}_{\mu}^{*}g,\sqrt{\epsilon}\cdot g)$. Further, we see that $\mathcal{\bar{F}}_{\mu}\mathcal{\bar{F}}_{\mu}^{*}=\mathcal{G}_{\mu}+\epsilon\mathcal{I}_{\mu}$. We can write:

[TABLE]

where $\mathcal{\bar{P}}_{\mu}=\mathcal{\bar{F}}_{\mu}^{*}(\mathcal{G}_{\mu}+\epsilon\mathcal{I}_{\mu})^{-1}\mathcal{\bar{F}}_{\mu}$. Note that $\mathcal{\bar{P}}_{\mu}$ is self adjoint. Correspondingly, let $\mathbf{\bar{F}}:\mathbb{C}^{s}\times L_{2}(\mu)\rightarrow L_{2}(\mu)$ be given by $\mathbf{\bar{F}}(f,g)=\mathbf{F}f+\sqrt{\epsilon}\cdot g$. We have $\mathbf{\bar{F}}^{*}g=(\mathbf{F}^{*}g,\sqrt{\epsilon}\cdot g)$. We can also write $\mathbf{\bar{P}}=\mathbf{\bar{F}}^{*}(\mathcal{G}_{\mu}+\epsilon\mathcal{I}_{\mu})^{-1}\mathcal{\bar{F}}_{\mu}$, and observe that $\mathbf{\bar{F}}^{*}=\mathbf{\bar{P}}\mathcal{\bar{F}}_{\mu}^{*}$.

With this notation in place we can rewrite the last term of (57) as:

[TABLE]

Using the fact that $\Re(z)\leq|z|$ for all $z\in\mathbb{C}$, and applying Cauchy-Schwarz to (C.3) and plugging back into (57) we have:

[TABLE]

We now bound $\|\mathbf{\bar{P}}^{*}(\mathbf{F}^{*}g^{\star}-\mathbf{\bar{y}},\sqrt{\epsilon}g^{\star})\|_{L_{2}(T)\times L_{2}(\mu)}$. If we had not sampled, this would equal:

[TABLE]

since $g^{\star}$ is the optimum of $\|\mathcal{\bar{F}}_{\mu}^{*}g-(\bar{y},0)\|_{L_{2}(T)\times L_{2}(\mu)}$ and thus $\mathcal{\bar{F}}_{\mu}^{*}g^{\star}-(\bar{y},0)$ is orthogonal to $\operatorname{range}(\mathcal{\bar{F}}_{\mu}^{*})$. We will show that after sampling, while the norm no longer equals [math], it is still small. The bound we give is analogous to standard approximate matrix multiplication results for finite dimensional matrices. Specifically, our proof follows that of Lemma 4 in [DKM06].

Claim 45** (Approximate Operator Application).**

With probability $\geq 1-\delta$:

[TABLE]

Proof.

For conciseness let $\mathcal{H}$ denote the space $L_{2}(T)\times L_{2}(\mu)$. Let $\varphi_{t}\in L_{2}(\mu)$ be given by $\varphi_{t}(\xi)=e^{-2\pi it\xi}$. Let $b^{\star}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\mathcal{F}_{\mu}^{*}g^{\star}-\bar{y}$ and $\mathbf{b}^{\star}\in\mathbb{C}^{s}$ be given by $\mathbf{b}^{\star}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\mathbf{F}^{*}g^{\star}-\mathbf{\bar{y}}.$ We can see that $\mathbf{b}^{\star}(j)=w_{j}\cdot[\langle\varphi_{t_{j}},g^{\star}\rangle_{\mu}-\bar{y}(t_{j})]$. We have:

[TABLE]

where the last equality follows since by (52), for any $\mathbf{x}\in\mathbb{C}^{s}$, $\mathbf{F}\mathbf{x}=\sum_{j=1}^{s}\varphi_{t_{j}}\cdot w_{j}\cdot\mathbf{x}(j)$. To simplify (61) we first bound, for any $g\in L_{2}(\mu)$, $\mathbb{E}\left[\langle g,\varphi_{t_{j}}\cdot w_{j}\cdot\mathbf{b}^{\star}(j)\rangle_{\mu}\right]$, recalling that $\mathbf{b}^{\star}(j)=w_{j}\cdot[\langle\varphi_{t_{j}},g^{\star}\rangle_{\mu}-\bar{y}(t_{j})]$. Let $p(t)=\frac{\tilde{\tau}_{\mu,\epsilon}(t)}{\tilde{s}_{\mu,\epsilon}}$ be the density with which we sample our time points $t_{1},\ldots,t_{s}$ and $w(t)=\sqrt{\frac{1}{sT\cdot p(t)}}$ be the reweighting factor we apply if we sample time $t$ (so $w_{j}=w(t_{j})$).

First we argue that we can apply Fubini’s theorem to switch the order of the double integration in $\mathbb{E}\left[\langle g,\varphi_{t_{j}}\cdot w_{j}\cdot\mathbf{b}^{\star}(j)\rangle_{\mu}\right]$ (over random instantiations of $\varphi_{t_{j}}\cdot w_{j}\cdot\mathbf{b}^{\star}(j)$ and within the inner product). Letting for $z\in L_{2}(\mu)$, $|z|\in L_{2}(\mu)$ be given by $|z|(\eta)=|z(\eta)|$ we have:

[TABLE]

which, noting that $\|\varphi_{t_{j}}\|_{\mu}=1$ gives:

[TABLE]

where the last line follows since $g\in L_{2}(\mu)$ so $\|g\|_{\mu}<\infty$ and since $\frac{1}{T}\int_{0}^{T}|\langle\varphi_{t},g^{\star}\rangle_{\mu}-y(t)|\,dt\leq\frac{1}{T}\int_{0}^{T}\left(|\langle\varphi_{t},g^{\star}\rangle_{\mu}-y(t)|^{2}+1\right)\,dt=\|\mathcal{F}_{\mu}^{*}g^{\star}-y\|_{T}^{2}+T\leq\|y\|_{T}^{2}<\infty$. Since we have established that $\mathbb{E}\left[\langle|g|,|\varphi_{t_{j}}\cdot w_{j}\cdot\mathbf{b}^{\star}(j)|\rangle_{\mu}\right]$ is finite we can apply Fubini’s theorem to compute:

[TABLE]

This in turn gives that

[TABLE]

and so for any $g\in L_{2}(\mu)$:

[TABLE]

Further, since $t_{1},\ldots,t_{s}$ are independent, the above gives that for $j\neq k$:

[TABLE]

We can apply (C.3) and (64) to expand out (61), giving:

[TABLE]

where the second to last line follows from (C.3) which gives

[TABLE]

Given the bound of (65) we can now expand out, using the fact that time $t$ is sampled with probability proportional to $\tilde{\tau}_{\mu,\epsilon}(t)$:

[TABLE]

Since $s=\Omega\left(\frac{\tilde{s}_{\mu,\epsilon}}{\delta}\right)$ we thus have via Markov’s inequality, with probability $\geq 1-\delta$,

[TABLE]

which completes the claim. Note that $64$ is an arbitrarily chosen constant, which can be made as small as we want by increasing the sample size $s$ by a constant factor. ∎

Plugging Claim 45 back into (59) gives:

[TABLE]

Applying the operator approximation bound of Lemma 43 with error $\Delta=1/4$ then gives:

[TABLE]

Finally, applying Claim 44 gives:

[TABLE]

Note that $\mathbb{E}[\tilde{B}]=B^{\star}$. So writing $\alpha=\tilde{B}-B^{\star}$ we have $|\alpha|\leq\frac{1}{\delta}\cdot B^{\star}$ with probability $1-\delta$. This completes the theorem. ∎

C.4 Frequency Subset Selection

We now prove the frequency subset selection guarantee Theorem 9 used in Section 5.1 to bound the leverage scores for general constraints $\mu$, by showing that $\mathcal{F}_{\mu}^{*}$ can be well approximated by an operator whose columns are spanned by just $O(s_{\mu,\epsilon})$ frequencies.

Theorem 9 (Frequency Subset Selection).

For some $s\leq\lceil 36\cdot s_{\mu,\epsilon}\rceil$ there exists a set of distinct frequencies $\xi_{1},\ldots,\xi_{s}\in\mathbb{C}$ such that, letting $\mathbf{C}_{s}:L_{2}(T)\rightarrow\mathbb{C}^{s}$ be defined by:

[TABLE]

and $\mathbf{Z}=(\mathbf{C}_{s}\mathbf{C}_{s}^{*})^{-1}\mathbf{C}_{s}\mathcal{F}_{\mu}^{*}$, for $\varphi_{t}\in L_{2}(\mu),\boldsymbol{\phi}_{t}\in\mathbb{C}^{s}$ with $\varphi_{t}(\xi)=e^{-2\pi it\xi}$ and $\boldsymbol{\phi}_{t}(j)=\varphi_{t}(\xi_{j})$:

[TABLE]

Our proof relies on the following spectral error bound for weighted frequency subset selection:

Lemma 46** (Frequency Subset Selection – Direct Spectral Approximation).**

For some $s\leq\lceil 36\cdot s_{\mu,\epsilon}\rceil$ there exists a set of distinct frequencies $\xi_{1},\ldots,\xi_{s}\in\mathbb{R}$ and positive weights $w_{1},\ldots,w_{s}\in\mathbb{R}$ such that letting $\mathbf{\bar{C}}_{s}:L_{2}(T)\rightarrow\mathbb{C}^{s}$ be given by:

[TABLE]

and letting $\widehat{\cal K}_{\mu}=\mathbf{\bar{C}}_{s}^{*}\mathbf{\bar{C}}_{s}$, we have

[TABLE]

Proof.

We prove a more general statement, in which we are given $0<\Delta<1$ and we select $s=\lceil 9s_{\mu,\epsilon}/\Delta^{2}\rceil$ frequencies $\xi_{1},\ldots,\xi_{s}\in\mathbb{R}$ and weights $w_{1},\ldots,w_{s}\in\mathbb{R}$ such that

[TABLE]

The claim follows by setting $\Delta=1/2$. We can assume that $\xi_{1},\ldots,\xi_{s}$ are distinct, since if $\xi_{i},\xi_{j}$ are equal, we can simply remove $\xi_{j}$ and update $w_{i}\leftarrow\sqrt{w^{2}_{i}+w^{2}_{j}}$, leaving $\widehat{\cal K}_{\mu}$ unchanged and only decreasing $s$.

The last condition is equivalent to

[TABLE]

Multiplying with $({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}$ on the left and right, we find that the condition is equivalent to:

[TABLE]

To shorten notation, we write ${\cal Z}=({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}{\cal K}_{\mu}({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}$ and $\widehat{\cal Z}=({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}\widehat{\cal K}_{\mu}({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}$. Given $\xi\in\mathbb{R}$, we define $\vartheta_{\xi}(t)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}e^{2\pi it\xi}$ ($\vartheta_{\xi}\in L_{2}(T)$). It is easy to verify that

[TABLE]

and

[TABLE]

Further define $\bar{\vartheta}_{\xi}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}\vartheta_{\xi}$. Since $({\cal K}_{\mu}+\epsilon{\cal I}_{T})^{-1/2}$ is self-adjoint and bounded, we have

[TABLE]

and

[TABLE]

We prove the existence of $\xi_{1},\dots,\xi_{s}$ and $w_{1},\dots,w_{s}$ using the deterministic selection process known as “BSS” [BSS14].292929We remark that unlike the process described in [BSS14], our existence proof does not trivially translate to an algorithm, since it involves a search over an infinite domain. Nevertheless, for our needs, existence suffices. In particular, we use a process that in essence is the same as the one described in [CNW16, Theorem 5 (arxiv version)]. Indeed, since $\|{\cal Z}\|_{\mathrm{op}}\leq 1$ and $\operatorname{tr}({\cal Z})=s_{\mu,\epsilon}$ the aforementioned results would suffice if we were dealing with matrices instead of operators. The rest of the proof extends these results to the operator case. Let

[TABLE]

and for $j=0,1,\dots,s$,

[TABLE]

The process we shall describe iteratively selects $\xi_{1},\xi_{2},\dots$ and unscaled weights $\tilde{w}_{1},\tilde{w}_{2},\ldots$ such that if we define $\widehat{{\cal Z}}^{(j)}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\sum^{j}_{i=1}\tilde{w}_{i}(\bar{\vartheta}_{\xi_{i}}\otimes\bar{\vartheta}_{\xi_{i}})$ the invariant

[TABLE]

is held. Let us write $s=\lceil 9s_{\mu,\epsilon}/\Delta^{2}\rceil$, so $s=Cs_{\mu,\epsilon}/\Delta^{2}$ for $C\geq 9$. If indeed we are able to select the frequencies and weights for $s$ steps such that this invariant holds, we shall have

[TABLE]

where we used the fact that ${\cal Z}\preceq{\cal I}_{T}$. Since $C\geq 9$ we have

[TABLE]

so by defining $w_{i}=\sqrt{\frac{3\Delta}{Cs_{\mu,\epsilon}}\tilde{w}_{i}}$ for $i=1,\dots,s$ we shall then have $\widehat{{\cal Z}}=\frac{3\Delta}{Cs_{\mu,\epsilon}}\widehat{{\cal Z}}^{(s)}$ thereby establishing the desired bound.

Thus, it suffices to show that we can select frequencies and weights iteratively so that (68) is maintained. In fact, the iterative selection process will maintain two additional invariants:

[TABLE]

All the invariants hold for $j=0$. Eq. (68) trivially holds for $j=0$. As for the integral,

[TABLE]

and similarly for the second invariant. In the above, the last equality is due to Claim 34.

Suppose by induction that the invariants for $j$. We prove that it is possible to pick a frequency $\xi$ and weight $w>0$ such that if we set $\xi_{j+1}=\xi$ and $\tilde{w}_{j+1}=w$ then the invariants will hold for $j+1$.

Fix $j$. For $t\geq 0$, let us denote

[TABLE]

where $M_{u}$ is defined for any $t$ (since the inverted operator is strictly positive and bounded, so invertible), and $M_{l}$ is defined for $t<1$. We can define $M_{l}$ for $t<1$ since $\widehat{{\cal Z}}^{(j)}-{\cal X}^{(j)}_{l}-t{\cal Z}\succ 0$ for $t<1$ as we now show. Due to Claim 34:

[TABLE]

Since $\widehat{{\cal Z}}^{(j)}-{\cal X}^{(j)}_{l}\succ 0$ (induction assumption), $(\widehat{{\cal Z}}^{(j)}-{\cal X}^{(j)}_{l})^{-1}$ is bounded so according to Claim 29, ${\cal Z}\preceq\widehat{{\cal Z}}^{(j)}-{\cal X}^{(j)}_{l}$, and then Claim 31 implies that $\widehat{{\cal Z}}^{(j)}-{\cal X}^{(j)}_{l}-t{\cal Z}\succ 0$.

Consider some fixed $\xi$. We first claim that for $w<1/\langle\bar{\vartheta}_{\xi},M_{u}(\delta_{u})\bar{\vartheta}_{\xi}\rangle_{T}$ we have $M_{u}(\delta_{u})^{-1}-w(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})\succ 0$. Obviously, the last statement holds for $w=0$, and due to continuity of $w\mapsto\langle x,(M_{u}(\delta_{u})^{-1}-w(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi}))x\rangle_{T}$ with respect to $w$, it will also hold for some interval around [math]. Let $w^{\star}$ be the maximal value such that for all $w\in[0,w^{\star})$ we have $M_{u}(\delta_{u})^{-1}-w(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})\succ 0$. Our goal is to show that $w^{\star}\geq 1/\langle\bar{\vartheta}_{\xi},M_{u}(\delta_{u})\bar{\vartheta}_{\xi}\rangle_{T}$. Assume by contradiction that $w^{\star}<1/\langle\bar{\vartheta}_{\xi},M_{u}(\delta_{u})\bar{\vartheta}_{\xi}\rangle_{T}$. For every $w\in[0,w^{\star})$, the operator $M_{u}(\delta_{u})^{-1}-w(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})$ is invertible, and we can apply a operator pseudo-inversion lemma due to Deng [Den11, Theorem 2.1] to find that

[TABLE]

Since we assumed $w^{\star}<1/\langle\bar{\vartheta}_{\xi},M_{u}(\delta_{u})\bar{\vartheta}_{\xi}\rangle_{T}$, clearly, there exists a $K$ such that for all $w\in[0,w^{\star})$ we have:

[TABLE]

Note that $M_{u}(\delta_{u})^{-1}-w^{\star}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})$ is not strictly positive for otherwise due to continuity we could have extended the interval, so there exists a $x$ with norm $1$ such that $\langle x,(M_{u}(\delta_{u})^{-1}-w^{\star}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi}))x\rangle<1/2K$. Let $w_{1},w_{2},\dots$ be a sequence which converges to $w^{\star}$, and let $y_{i}=(M_{u}(\delta_{u})^{-1}-w_{i}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})^{1/2})x$. We now have $\langle y_{i},y_{i}\rangle_{T}=\langle x,(M_{u}(\delta_{u})^{-1}-w_{i}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})x\rangle_{T}\to\langle x,(M_{u}(\delta_{u})^{-1}-w^{\star}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi}))x\rangle_{T}<1/2K$ as $i\to\infty$. However $\langle y_{i},(M_{u}(\delta_{u})^{-1}-w_{i}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi}))^{-1}y_{i}\rangle_{T}=\langle x_{i},x_{i}\rangle_{T}=1$ which contradicts the bound on $(M_{u}(\delta_{u})^{-1}-w_{i}(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi}))^{-1}$.

Thus, if we picked $\xi$ and $w<1/\langle\bar{\vartheta}_{\xi},M_{u}(\delta_{u})\bar{\vartheta}_{\xi}\rangle_{T}$ for the step, we shall have $\widehat{{\cal Z}}^{(j+1)}-{\cal X}^{(j)}_{l}=M_{u}(\delta_{u})^{-1}-w(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})\succ 0$ as required, and the upper invariant will translate to

[TABLE]

which is equivalent to

[TABLE]

The induction hypothesis is

[TABLE]

so the upper invariant is held if

[TABLE]

Consider any $\eta\in\mathbb{R}$, and let $f_{\eta}(y)\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\langle\bar{\vartheta}_{\eta},M_{u}(y)\bar{\vartheta}_{\eta}\rangle_{T}$. Using the operator inversion formula, we have for any $t_{2}\geq t_{1}$:

[TABLE]

From this equation we see that

[TABLE]

Furthermore, since for $t_{2}>t_{1}$ we have ${\cal I}_{T}+t_{2}{\cal Z}^{1/2}M_{u}(t_{1}){\cal Z}^{1/2}\succeq{\cal I}_{T}+t_{1}{\cal Z}^{1/2}M_{u}(t_{1}){\cal Z}^{1/2}$ and both operators are strictly positive and bounded, then $({\cal I}_{T}+t_{1}{\cal Z}^{1/2}M_{u}(t_{1}){\cal Z}^{1/2})^{-1}\preceq({\cal I}_{T}+t_{2}{\cal Z}^{1/2}M_{u}(t_{1}){\cal Z}^{1/2})^{-1}$, and we can easily verify that $f_{\eta}$ is convex. Thus,

[TABLE]

After integrating on both sides, we have the bound

[TABLE]

Using this bound in (69) and rearranging, we find that for any $\xi$, the upper invariant is held if we select $w$ such that

[TABLE]

Note that if this is held, we also have $w<1/\langle\bar{\vartheta}_{\xi},M_{u}(\delta_{u})\bar{\vartheta}_{\xi}\rangle_{T}$, as previously required.

We now consider the lower invariants. If we picked $\xi$ and $w>0$ for the step, then $\widehat{{\cal Z}}^{(j+1)}-{\cal X}_{l}^{(j)}=M_{l}(\delta_{l})^{-1}+w(\bar{\vartheta}_{\xi}\otimes\bar{\vartheta}_{\xi})\succeq M_{l}(\delta_{l})^{-1}\succ 0$ as long $\delta_{l}<1$ which holds for our choice of $\delta_{l}$. So the left part of (68) will hold regardless of how we choose $\xi$ and $w>0$. As for the lower trace bound, it translates to

[TABLE]

Applying another variant of operator pseudo-inversion lemma [Oga88, Theorem 2], we find that the last condition is equivalent to

[TABLE]

The induction hypothesis is

[TABLE]

so the lower invariant is held if

[TABLE]

Similarly to before, by using the convexity of each integrand, we can bound

[TABLE]

Using this bound in (71) and rearranging, we find that for any $\xi$, the lower invariant is held if we select $w$ such that

[TABLE]

Thus, we need to show that there exists a $\xi$ and $w$ such that both (70) and (72) hold. However, for a given $\xi$, such a $w$ will surely exist if

[TABLE]

Thus, it it suffices to show that there exists a $\xi$ for which the last inequality holds. To show that such a $\xi$ exists, we will show that the inequality holds for the integral of both sides with respect to $\mu$ measure. This will guarantee the existence of such a $\xi$ since the Lebesgue integral is strictly positive for non-negative functions. We compute:

[TABLE]

Similarly,

[TABLE]

${\cal Z}$ is self-adjoint and positive definite, so the operator pseudo-inversion lemma [Oga88, Theorem 2] implies that $M_{u}(\delta_{u})\preceq M_{u}(0)$, so by the induction hypothesis

[TABLE]

We now consider the lower invariant. We already showed that ${\cal Z}\preceq\widehat{{\cal Z}}^{(j)}-{\cal X}^{(j)}_{l}$, so as long as $\delta_{l}\leq 1/2$ we will have:

[TABLE]

where we used Claim 30. So there will be a gap in the value of the integrals (as desired), if

[TABLE]

which is the case for our selection of $\delta_{l}$ and $\delta_{u}$. ∎

From Lemma 46 we can prove a stronger spectral error bound for the projection onto the range of $\mathbf{\bar{C}}_{s}$.

Lemma 47** (Frequency Subset Selection – Projection Based Spectral Approximation).**

For some $s\leq\lceil 36\cdot s_{\mu,\epsilon}\rceil$ there exists a set of distinct frequencies $\xi_{1},\ldots,\xi_{s}\in\mathbb{C}$ such that letting $\mathbf{C}_{s}:L_{2}(T)\rightarrow\mathbb{C}^{s}$ and $\mathbf{Z}:L_{2}(\mu)\rightarrow\mathbb{C}^{s}$ be defined as in Theorem 9 and $\widehat{\cal G}_{\mu}=\mathbf{Z}^{*}\mathbf{C}_{s}\mathbf{C}_{s}^{*}\mathbf{Z}$,

[TABLE]

Proof.

Let $\xi_{1},\ldots,\xi_{s}\in\mathbb{C}$ and $w_{1},\ldots,w_{s}\in\mathbb{R}$ be the frequencies and weights shown to exist in Lemma 46 and let $\mathbf{\bar{C}}_{s}$ be as defined in that lemma (note that $\mathbf{\bar{C}}_{s}$ is identical to $\mathbf{C}_{s}$ except with its rows weighted by $w_{1},\ldots,w_{s}$.) First note that for any $g\in L_{2}(\mu)$,

[TABLE]

where the inequality follows from observing that $\mathbf{C}_{s}^{*}(\mathbf{C}_{s}\mathbf{C}_{s}^{*})^{-1}\mathbf{C}_{s}$ is an orthogonal projection. Thus $\widehat{\cal G}_{\mu}\preceq{\cal G}_{\mu}$. It remains to show that ${\cal G}_{\mu}\preceq\widehat{\cal G}_{\mu}+\epsilon\mathcal{I}_{\mu}$. Let $\mathcal{\bar{P}}=\mathcal{I}_{T}-\mathbf{C}_{s}^{*}(\mathbf{C}_{s}\mathbf{C}^{*})^{-1}\mathbf{C}_{s}$ be the projection to the orthogonal complement of $\mathbf{C}_{s}^{*}$’s range and let $\widehat{\cal K}_{\mu}=\mathbf{\bar{C}}_{s}^{*}\mathbf{\bar{C}}_{s}$ be as defined in Lemma 46. Rearranging the guarantee of Lemma 46 gives

[TABLE]

which immediately gives

[TABLE]

Note that $\mathbf{\bar{C}}_{s}\mathcal{\bar{P}}=0$ (since $\mathcal{\bar{P}}$ is an orthogonal projection onto $\ker(\mathbf{C}_{s})=\ker(\mathbf{\bar{C}}_{s})$) and so $\mathcal{\bar{P}}\widehat{\cal K}_{\mu}\mathcal{\bar{P}}=0$, giving:

[TABLE]

Note that $\mathcal{\bar{P}}{\cal K}_{\mu}\mathcal{\bar{P}}=\mathcal{\bar{P}}\mathcal{F}_{\mu}^{*}\mathcal{F}_{\mu}\mathcal{\bar{P}}$ and

[TABLE]

Thus by (74) we also have $\mathcal{G}_{\mu}-\widehat{\mathcal{G}_{\mu}}\preceq\epsilon\mathcal{I}_{\mu}$ (since the norm of an operator and its adjoint are the same so $\mathcal{\bar{P}}{\cal K}_{\mu}\mathcal{\bar{P}}\preceq\epsilon\mathcal{I}_{T}\implies\mathcal{F}_{\mu}\mathcal{\bar{P}}\mathcal{F}_{\mu}^{*}\preceq\epsilon\mathcal{I}_{\mu}$), which completes the lemma. ∎

Finally, from Lemma 47 we can prove the frequency subset selection guarantee of Theorem 9.

Proof of Theorem 9.

We consider the same set of frequencies $\xi_{1},\ldots,\xi_{s}$ shown to exist in Lemma 47 and the corresponding operators $\mathbf{C}_{s}$, $\mathbf{Z}$. We show that these frequencies satisfy the guarantee of Theorem 9. First, we note that

[TABLE]

(In the above, we abuse notation and use $\boldsymbol{\phi}_{t}$ to denote both the vector defined in the Theorem statement, and the operator $x\in\mathbb{C}\mapsto x\boldsymbol{\phi}_{t}$). From Claim 34:

[TABLE]

We have,

[TABLE]

From Claim 33:

[TABLE]

Next, consider $\frac{1}{T}\int^{T}_{0}\boldsymbol{\phi}_{t}\otimes\varphi_{t}\,dt$. For any $\alpha$,

[TABLE]

where the integral on the left is a weak vector integral. Since for every $g\in L_{2}(T)$,

[TABLE]

and for every $\alpha\in L_{2}(\mu)$, $[\mathcal{F}_{\mu}^{*}\alpha](t)=\langle\varphi_{t},\alpha\rangle_{\mu}$, we have $\frac{1}{T}\int^{T}_{0}\boldsymbol{\phi}_{t}\otimes\varphi_{t}\,dt=\mathbf{C}_{s}\mathcal{F}_{\mu}^{*}$, so

[TABLE]

Combining the previous observations, we find that

[TABLE]

Let $v_{1},\ldots,v_{2s_{\mu,\epsilon}}\in L_{2}(\mu)$ be the eigenfunctions of $\mathcal{G}_{\mu}$ corresponding to its top $2s_{\mu,\epsilon}$ eigenvalues. Define $\mathbf{X}:L_{2}(\mu)\rightarrow\mathbb{C}^{2s_{\mu,\epsilon}}$ as: for $g\in L_{2}(\mu)$, $[\mathbf{X}g](j)=\langle v_{j},g\rangle_{\mu}$. Note that

[TABLE]

since $\mathbf{C}_{s}\mathbf{Z}\mathbf{Z}^{*}\mathbf{C}_{s}\succeq\mathbf{C}_{s}\mathbf{Z}\mathbf{X}^{*}\mathbf{X}\mathbf{Z}^{*}\mathbf{C}_{s}$ ($\mathbf{X}^{*}\mathbf{X}$ is a projection, so $\mathbf{X}^{*}\mathbf{X}\preceq{\cal I}_{\mu}$). So we can bound:

[TABLE]

Let $i_{\epsilon}$ be the smallest $i$ with $\lambda_{i}(\mathcal{G}_{\mu})\leq\epsilon$. We have:

[TABLE]

Thus we can bound $\operatorname{tr}(\mathcal{G}_{\mu}-\mathbf{X}^{*}\mathbf{X}{\mathcal{G}_{\mu}}\mathbf{X}^{*}\mathbf{X})$ as:

[TABLE]

where the last bound follows from the fact that $s_{\mu,\epsilon}\geq\sum_{i=i_{\epsilon}+1}^{\infty}\frac{\lambda_{i}(\mathcal{G}_{\mu})}{\lambda_{i}(\mathcal{G}_{\mu})+\epsilon}\geq\sum_{i=i_{\epsilon}+1}^{\infty}\frac{\lambda_{i}(\mathcal{G}_{\mu})}{2\epsilon}$.

We can also bound $\operatorname{tr}(\mathbf{X}^{*}\mathbf{X}(\mathcal{G}_{\mu}-\widehat{\mathcal{G}_{\mu}})\mathbf{X}^{*}\mathbf{X})$ using Lemma 47. Since $\mathcal{G}_{\mu}\leq\widehat{\mathcal{G}_{\mu}}+\epsilon\mathcal{I}_{\mu}$ we have:

[TABLE]

Plugging (76) and (77) back into (C.4) we have:

[TABLE]

which completes the theorem. ∎

Appendix D Tight Statistical Dimension Bound for Bandlimited Functions

In Section 5 we demonstrate, perhaps surprisingly, that a simple function $\tilde{\tau}_{\mu,\epsilon}(t)$ (defined in Theorem 17) exists for any $\mu$ that upper bounds $\tau_{\mu,\epsilon}(t)$ and has $\tilde{s}_{\mu,\epsilon}=\tilde{O}(s_{\mu,\epsilon})$. Combined with Theorem 7 this yields our main algorithmic result Theorem 3, which shows that we can achieve $O\left(s_{\mu,\epsilon}\log^{2}(s_{\mu,\epsilon}\right)$ sample complexity with just $\tilde{O}\left(s_{\mu,\epsilon}^{\omega}\right)$ runtime.

Instantiating Theorem 3 using the approximate ridge leverage function of Theorem 17 requires an upper bound on $s_{\mu,\epsilon}$. In this section we show how to bound $s_{\mu,\epsilon}$ when $\mu$ is uniform measure on some interval – i.e., when our interpolation problem is over bandlimited functions. In Section E we leverage this result to bound $s_{\mu,\epsilon}$ for a number of other important priors, including for multiband, Gaussian, and Cauchy-Lorentz.

Beyond letting us upper bound $s_{\mu,\epsilon}$ to apply Theorem 3, our proof for bandlimited functions is constructive, giving a simple upper bound for $\tau_{\mu,\epsilon}(t)$ for any $t$. This upper bound can be plugged directly into Algorithm 1 and Theorem 7 to give a tightening of Theorem 3 by a logarithmic factor in the bandlimited case. Like our general result, the proof is based on the definition of leverage scores given in (11). This definition makes it clear that, to upper bound $\tau_{\mu,\epsilon}(t)$, it suffices to show that a function with Fourier support controlled by $\mu$ cannot “spike” too extremely at time $t$.

For bandlimited functions, we obtain a smoothness bound by introducing and applying a Bernstein type smoothness bound for low-degree polynomials and relying on the fact that any bandlimited function is well approximated by a low-degree polynomial. This approach mirrors the general proof in Section 5, which uses a more sophisticated smoothness bound for sparse Fourier functions.

Our result for bandlimited function is as follows:

Theorem 48**.**

Let $\mu$ be the uniform measure on $[-F,F]$. Let $q=\lceil 16\pi eFT+2\log(1/\epsilon)+11\rceil$. For all $t\in[0,T]$, let the approximate ridge leverage function $\tilde{\tau}_{\mu,\epsilon}$ equal:

[TABLE]

For any $\epsilon\leq 1,F,T$, $\tilde{\tau}_{\mu,\epsilon}(t)$ satisfies:

$\tilde{\tau}_{\mu,\epsilon}(t)\geq\tau_{\mu,\epsilon}(t)$. 2. 2.

$\int_{0}^{T}\tilde{\tau}_{\mu,\epsilon}(t)dt\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\tilde{s}_{\mu,\epsilon}=O\left(FT+\log(1/\epsilon)\right).$**

Thus we have $s_{\mu,\epsilon}\leq\tilde{s}_{\mu,\epsilon}=O\left(FT+\log(1/\epsilon)\right)$.

Combined with Theorem 7, Theorem 48 immediately gives:

Corollary 49**.**

Let $\mu$ be the uniform measure on $[-F,F]$. Using $\tilde{\tau}_{\mu,\epsilon}$ as defined in Theorem 48, Algorithm 1 returns $t_{1},\ldots,t_{s}\in[0,T]$ and $\mathbf{z}\in\mathbb{C}^{s}$ such that $\tilde{y}(t)=\sum_{i=1}^{s}\mathbf{z}(i)\cdot k_{\mu}(t_{i},t)$ satisfies with probability $\geq 1-\delta$:

[TABLE]

The algorithm queries $y+n$ at $s$ points and runs in $O(s^{\omega})$ time where $s=O([FT+\log(1/\epsilon)]\cdot[\log(FT+\log(1/\epsilon))+1/\delta])$. The output $\tilde{y}(t)$ can be evaluated using Algorithm 2 in $O(s)$ time.

Proof.

The corollary follows immediately from Theorem 7 after noting that

• •

$Z=O(1)$ since, as shown in Appendix F, $k_{\mu}(t_{1},t_{2})=\frac{\sin(2\pi F(t_{1}-t_{2}))}{2\pi F(t_{1}-t_{2})}$ and so can be computed in $O(1)$ arithmetic operations.

• •

$W=O(1)$ since to sample points proportional to $\tilde{\tau}_{\mu,\epsilon}(t)$, we must sample a mixture of the uniform distribution and the distribution with density proportional to $\frac{1}{\sqrt{\min(t,T-t)/T}}$. It suffices to show that we can sample from the later in $O(1)$ time, and in fact that we can sample $t\in[0,1/2]$ with probability proportional to $\frac{1}{\sqrt{t}}$ in $O(1)$ time, since we can then symmetrize and scale this distribution. We can accomplish this with inverse transform sampling. Our density is $p(t)=\frac{1}{2\sqrt{2t}}$ and so its cumulative distribution function is $C(t)=\sqrt{t/2}$. Thus we can sample $z$ uniformly in $[0,1]$ and return $C^{-1}(z)=2z^{2}$, which will be a sample from the desired distribution. This can be done with $O(1)$ arithmetic computations.

∎

D.1 Smoothness bounds for polynomials

As mentioned our main techniques tool is a Bernstein type smoothness bounds for low-degree polynomials. In general, low-degree polynomials are smoother than high-degree polynomials, and thus cannot spike as sharply. There are a number of ways to formalize this statement. The well known Markov brother’s inequality and Bernstein inequality bound the maximum derivative of a polynomial by a function of the polynomial’s degree and it’s maximum value on an interval.

To bound leverage scores, we are interested in a slightly different metric of smoothness. In particular, we need to bound the maximum squared value of a polynomial by its average squared value on $[0,T]$. We can use standard properties of the Legendre polynomials to prove: s

Claim 50**.**

For any degree $d$ polynomial $p(\cdot)$ with complex coefficients and $t\in[0,T]$, let $r=\frac{\min\left(t,T-t\right)}{T}$. Then:

[TABLE]

This bound is tighter for points near the center of the interval $[0,T]$ and goes to infinity near the edges. Using the Markov brother’s inequality, it’s possible to obtain a fixed up bound of $O(d^{2})$, which is tighter for small values of $r$. However, this won’t be necessary for our purposes. We note that, when $t=T/2$, the upper bound on $p(t)^{2}$ improves to $O(d)$ times the average squared value of $p$, quadratically better than an $O(d^{2})$ bound. This improvement is nearly optimal: the upper bound of Claim 50 is matched up to a logarithmic factor by an appropriately scaled and shifted Chebyshev polynomial of the first kind applied to $[T/2-t]^{2}$ (see e.g. [FMMS16] for a construction).

Proof of Claim 50.

The claim follows from properties of the standard orthogonal Legendre polynomials, which are denoted $P_{0},P_{1},\ldots$ and defined via the recurrence relation:

[TABLE]

The Legendre polynomials are orthogonal over the interval $[-1,1]$ with respect to the constant weight function. In particular, they satisfy

[TABLE]

where $\delta_{m,n}$ is the Kronecker delta function. Additionally, for $x\in[-1,1]$, $|P_{j}(x)|\leq 1$ for all $j$.

Using these facts we can show that for any degree $d$ polynomial $p(\cdot)$, interval $[a,b]$, and $x\in[a,b]$:

[TABLE]

where $r=\frac{\min\left(|a-x|,|b-x|\right)}{(b-a)}$. Setting $a=0$ and $b=T$ gives the claim.

We begin by noting that, without loss of generality, we can assume that $a=-1$ and $b=1$. In particular, shift and stretch $p(x)$ by defining $g(x)=p\left(\frac{2(x-a)}{b-a}-1\right)$. $g$ has degree $d$ and the maximum of $|g(x)|^{2}$ for $x\in[-1,1]$ is the same as the maximum of $|p(x)|^{2}$ for $x\in[a,b]$. Additionally, $\frac{\int_{-1}^{1}|g(t)|^{2}\,dt}{2}=\frac{\int_{a}^{b}|p(t)|^{2}\,dt}{(b-a)}$. Accordingly, to prove the claim it suffices to prove that, for any degree $d$ polynomial $g$,

[TABLE]

Our proof depends on a Bernstein type inequality for Legendre polynomials, which can be found in [Lor83]. Specifically, for all $j=0,1,2,\ldots$ and any $x\in[-1,1]$ it holds that:

[TABLE]

Writing $g$ in the Legendre basis:

[TABLE]

we have from (80) that

[TABLE]

and thus

[TABLE]

The last equality step follows from (78). Finally, let $q=\min\left(|-1-x|,|1-x|\right)$ and note that

[TABLE]

As defined, $r=q/2$ Plugging into (D.1) we have a final bound of

[TABLE]

which establishes (79) and thus the claim. ∎

D.2 Smoothness bounds for bandlimited functions

With Claim 50 in place, we are now ready to prove our main result for bandlimited functions.

Proof of Theorem 48.

Following Definition 3, our goal is to choose $\tilde{\tau}_{\mu,\epsilon}$ to satisfy:

[TABLE]

for any $\alpha$. Let $z=\mathcal{F}_{\mu}\alpha$. Expanding $e^{-2i\pi\xi t}$ using its Maclaurin series and letting $d$ be some degree parameter that we will fix later, we write $z$ as the sum of two functions, $a$ and $b$:

[TABLE]

Note that $a$ is a degree $d$ polynomial with complex coefficients. So by Claim 50,

[TABLE]

Turning our attention to $b$, we see that:

[TABLE]

The second to last step uses Cauchy-Schwarz inequality. Finally using that for all $j$, $j!\geq(j/e)^{j}$, for any $d\geq 4\pi eFT$:

[TABLE]

So, if we take $d=\lceil 4\pi eFT+\log(1/\epsilon)/2+1\rceil$, it follows from (85) and (86) that

[TABLE]

It follows that $\|b\|_{T}\leq\frac{\sqrt{\epsilon}}{2}\|\alpha\|_{\mu}$. Using the decomposition of (D.2) and the fact that for any real nonnegative $c,d$, $c^{2}+d^{2}\leq(c+d)^{2}$, and for any complex $e,f$, $|e+f|^{2}\leq 2|e|^{2}+2|f|^{2}$:

[TABLE]

It follows from (84) that:

[TABLE]

In Theorem 48 we set $q=\lceil 16\pi eFT+2\log(1/\epsilon)+11\rceil$. We have $q\geq 4\cdot\lceil 4\pi eFT+\log(1/\epsilon)/2+2\rceil=4(d+1)$ since, for any $x$, $\lceil 4x+3\rceil\geq 4\lceil x\rceil$. Recalling that $z=\mathcal{F}_{\mu}\alpha$, it follows $\tilde{\tau}_{\mu,\epsilon}$ defined in that theorem satisfies (82) for any $\alpha$. It remains to bound the total measure of our approximate ridge leverage function, $\tilde{s}_{\mu,\epsilon}$. To do so, note that:

[TABLE]

We can compute:

[TABLE]

This bound establishes the theorem. ∎

Appendix E Statistical dimension for common Fourier constraints

In this section we leverage Theorem 48 to give upper bounds on the statistical dimensions of a number common priors $\mu$ used for Fourier constrained interpolation, including multiband, Gaussian, and Cauchy-Lorentz priors. We start by giving two simple lemmas that we use to translate our bound for bandlimited functions to these more general priors.

Lemma 51** (Statistical Dimension of Sum of Measures).**

Let $\mu_{1},\mu_{2},\cdots\mu_{s}$ be finite measures on $\mathbb{R}$. Let $\mu$ be a probability measure defined by $\mu=\mu_{1}+\mu_{2}+\cdots+\mu_{s}$.

[TABLE]

Proof.

We can see from Definition 2 that for $\mu=\mu_{1}+\ldots+\mu_{s}$ the kernel operator $\mathcal{K}_{\mu}$ satisfies $\mathcal{K}_{\mu}=\sum_{i=1}^{s}\mathcal{K}_{\mu_{i}}$. We can thus bound:

[TABLE]

The second to last inequality follows since $0\preceq\mathcal{K}_{\mu_{i}}\preceq\mathcal{K}_{\mu}$, so $0\prec\mathcal{K}_{\mu_{i}}+\epsilon\mathcal{I}_{T}\preceq\mathcal{K}_{\mu}+\epsilon\mathcal{I}_{T}$ and $(\mathcal{K}_{\mu}+\epsilon\mathcal{I}_{T})^{-1}\preceq(\mathcal{K}_{\mu_{i}}+\epsilon\mathcal{I}_{T})^{-1}$ by Claim 30. Letting $e_{1},e_{2}$ be an orthonormal basis for $L_{2}(T)$, we thus have:

[TABLE]

This completes the lemma. ∎

Lemma 52** (Statistical Dimension of Scaled Measures).**

Let $\mu$ be a measure on $\mathbb{R}$. For any parameter $\gamma>0$, we have:

[TABLE]

Proof.

From Definition 2, we can see that $\mathcal{K}_{(\mu/\gamma)}=\frac{1}{\gamma}\mathcal{K}_{\mu}$ and thus has eigenvalues equal to $\lambda_{1}(\mathcal{K}_{\mu})/\gamma,\lambda_{2}(\mathcal{K}_{\mu})/\gamma,\ldots$ We can thus compute:

[TABLE]

∎

We now use Lemmas 51 and 52 to prove our statistical dimension bounds. We first start with multiband Fourier constraints, showing that the statistical dimension is roughly proportional to the total length of all the frequency bands times the time domain window size, intuitively matching the Landau rate for asymptotic recovery of multiband functions [Lan67a].

Theorem 53** (Multiband Statistical Dimension).**

Consider a set of $s$ disjoint frequency bands, $I_{1},I_{2},\cdots,I_{s}$, and suppose that the length of the band $I_{i}$ is denoted by $F_{i}$. Let $\mu$ be the measure which induces a uniform probability density on $I_{1}\cup I_{2}\cup\cdots\cup I_{s}$. We have:

[TABLE]

Proof.

For every $i$, let $\mu_{i}$ be the measure defined by $\mu_{i}(A)=\mu(A\cap I_{i})$. Note that we have $\mu=\sum_{i}\mu_{i}$ and so can invoke Lemma 51, giving:

[TABLE]

If $\mu_{i}$ gave a uniform probability measure on frequency band $I_{i}$ (i.e., if we had $\mu_{i}(\mathbb{R})=1$), we could use the result of Theorem 48 to bound $s_{\mu_{i},\epsilon}=O(F_{i}T+\log(1/\epsilon))$. This is not the case, but we can instead let $\gamma_{i}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\mu_{i}(\mathbb{R})\leq 1$. By Lemma 52,

[TABLE]

Now $\mu_{i}/\gamma_{i}$ is a uniform probability measure on $I_{i}$, so we can invoke Theorem 48 giving:

[TABLE]

Plugging this bound in (87) and using that $\gamma_{i}\leq 1$ we obtain:

[TABLE]

completing the theorem.∎

We next bound the statistical dimension of Gaussian measure.

Theorem 54** (Gaussian Statistical Dimension).**

Let $\mu$ induce the Gaussian probability distribution with standard deviation $F$ defined by $d\mu(\xi)=\frac{1}{\sqrt{2\pi F^{2}}}e^{-\xi^{2}/2F^{2}}d\xi$. We have:

[TABLE]

Proof.

Let $I_{h}$ be the interval defined by $I_{h}=\{\xi\in\mathbb{R}:|\xi|\leq F\sqrt{\log(1/\epsilon)}\}$. We decompose $\mu$ into two measures $\mu_{h}$ and $\mu_{t}$ as follows:

[TABLE]

We can see that $\mu=\mu_{h}+\mu_{t}$ and so by Lemma 51, $s_{\mu,\epsilon}\leq s_{\mu_{t},\epsilon}+s_{\mu_{h},\epsilon}.$ For $\mu_{t}$ we have:

[TABLE]

where the last bound follows from a Chernoff bound, giving $1-\operatorname{erf}(x)\leq 2e^{-x^{2}}$ [Wai18]. This lets us crudely bound:

[TABLE]

where the first ineguality is because $\|(\mathcal{K}_{\mu_{t}}+\epsilon\mathcal{I}_{T})^{-1}\|_{\mathrm{op}}\leq 1/\epsilon$.

We next bound the statistical dimension of $\mu_{h}$. Let $\tilde{\mu}_{h}$ be a uniform measure on $I_{h}$, with $d\mu(\xi)=\frac{1}{\sqrt{2\pi F^{2}}}d\xi$ for all $\xi\in I_{h}$. Note that $d\tilde{\mu}_{h}(\xi)\geq d{\mu}_{h}(\xi)$ for all $\xi\in I_{h}$ which gives that $K_{\mu_{h}}\preceq K_{{\tilde{\mu}}_{h}}$ and so $s_{\mu_{h},\epsilon}\leq s_{\tilde{\mu}_{h},\epsilon}.$

Let $\gamma\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\tilde{\mu}_{h}(\mathbb{R})=\sqrt{\frac{2\log(1/\epsilon)}{\pi}}$. By Lemma 52, $s_{\tilde{\mu},\epsilon}=s_{(\tilde{\mu}/\gamma),(\epsilon/\gamma)}.$ Since $\tilde{\mu}/\gamma$ is a uniform probability measure on $I_{h}$, we can invoke Theorem 48 to give:

[TABLE]

where the last equality follows from the fact that $\gamma=O(\sqrt{\log(1/\epsilon)})$. Combining (88) and (89) and applying Lemma 51 we have:

[TABLE]

which completes the theorem. ∎

Finally, we bound the statistical dimension of the Cauchy-Lorentz measure.

Theorem 55**.**

Let $\mu$ induce the Cauchy-Lorentz probability distribution with scale parameter $F$ defined by $d\mu(\xi)=\frac{1}{\pi F\left[1+\left(\frac{\xi}{F}\right)^{2}\right]}d\xi$. We have:

[TABLE]

Proof.

As in the proof of Theorem 54 we define $I_{h}$ to be the interval $I_{h}=\{\xi\in\mathbb{R}:|\xi|\leq F/\sqrt{\epsilon}\}$. We decompose $\mu$ into two measures $\mu_{h}$ and $\mu_{t}$ as follows:

[TABLE]

We have $\mu=\mu_{h}+\mu_{t}$ by Lemma 51, $s_{\mu,\epsilon}\leq s_{\mu_{t},\epsilon}+s_{\mu_{h},\epsilon}.$ For $\mu_{t}$ we have:

[TABLE]

As in (88) we can thus bound:

[TABLE]

We next bound the statistical dimension of $\mu_{h}$. Let $\tilde{\mu}_{h}$ be a uniform measure on $I_{h}$ with $d\mu(\xi)=\frac{1}{\pi F}$ for all $\xi\in I_{h}$. As in the proof of Theorem 54, $d\tilde{\mu}_{h}(\xi)\geq d{\mu}_{h}(\xi)$ for all $\xi\in I_{h}$ which gives that $K_{\mu_{h}}\preceq K_{{\tilde{\mu}}_{h}}$ and so $s_{\mu_{h},\epsilon}<s_{\tilde{\mu}_{h},\epsilon}.$

Let $\gamma\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}\tilde{\mu}_{h}(\mathbb{R})=\frac{2}{\pi\sqrt{\epsilon}}$. By Lemma 52, $s_{\tilde{\mu},\epsilon}=s_{(\tilde{\mu}/\gamma),(\epsilon/\gamma)}.$ Since $\tilde{\mu}/\gamma$ is a uniform probability measure on $I_{h}$, we can invoke Theorem 48 to give:

[TABLE]

where the last equality follows from the fact that $\gamma=O(1/\sqrt{\epsilon})$. Combining (90) and (91) and applying Lemma 51 we have:

[TABLE]

which completes the theorem. ∎

Appendix F Kernel computation for common Fourier constraints

Algorithm 1 and the corresponding Theorem 3 assumes the ability to compute the kernel function $k_{\mu}(t_{1},t_{2})=\int_{\xi\in\mathbb{R}}e^{-2\pi i(t_{1}-t_{2})}d\mu(\xi)$. In this section we give close forms for the kernel functions of popular measures $\mu$, including all those whose statistical dimension we bound in Appendix E.

Bandlimited Fourier Constraint:

When $\mu$ is the uniform measure on frequencies in $[-F,F]$,

[TABLE]

So, $k_{\mu}$ is the sinc kernel.

Multiband Fourier Constraint:

Consider a set of $s$ disjoint frequency bands, $I_{1},I_{2},\cdots,I_{s}$, where $I_{j}=[c_{j}-F_{j},c_{j}+F_{j}]$. Let $\mu$ be the uniform measure on $I_{1}\cup I_{2}\cup\cdots\cup I_{s}$. Then we have:

[TABLE]

Gaussian Fourier Constraint:

When $\mu$ induces the Gaussian probability distribution with standard deviation $F$ defined by $d\mu(\xi)=\frac{1}{\sqrt{2\pi F^{2}}}e^{-\xi^{2}/2F^{2}}d\xi$, then

[TABLE]

So, $k_{\mu}$ is the Gaussian kernel.

Cauchy-Lorentz Fourier Constraint:

When $\mu$ induces the Cauchy-Lorentz probability density with scale parameter $F$ defined by $d\mu(\xi)=\frac{1}{\pi F\left[1+\left(\frac{\xi}{F}\right)^{2}\right]}\cdot d\xi$, we have:

[TABLE]

In the machine learning literature, $k_{\mu}$ is known as the Laplacian kernel.

Appendix G Signal Reconstruction as Bayesian Estimation

In this section, we show how, as an alternative to Problem 1, we can formulate signal fitting as a Bayesian estimation problem, where the signal $y$ is a stationary stochastic process and the measure $\mu$ (which we assume to be symmetric about [math] throughout this section so that $k_{\mu}(t_{1},t_{2})$ is real valued) corresponds to a prior on $y$’s power spectral density. This form of prior is commonly used in statistical signal processing, kriging, and machine learning applications [HS93, Rip05, RW06]: We first define a stationary Gaussian process:

Definition 7** (Stationary Gaussian Process [RW06]).**

A stochastic process $y:\mathbb{R}\rightarrow\mathbb{R}$ is a Gaussian process if for any finite collection $t_{1},\ldots t_{s}\in\mathbb{R}$, $y(t_{1}),\ldots,y(t_{s})$ is distributed as a multivariate Gaussian. $y$ is a stationary Gaussian process if the mean $\mathbb{E}[y(t)]$ is independent of $t$ and the autocorrelation $\mathbb{E}[y(t_{1})\cdot y(t_{2})]$ depends only on $t_{1}-t_{2}$.

We now define the specific Gaussian process prior we consider:

Definition 8** (Gaussian Process Prior).**

Consider a symmetric probability density function $p_{\mu}:\mathbb{R}\rightarrow\mathbb{R}^{+}$ and the associated measure $\mu$ corresponding to $p_{\mu}$. We say that a stochastic process $y:\mathbb{R}\rightarrow\mathbb{R}$ is distributed according to $\mathcal{D}_{\mu}$ if $y$ is distributed as a stationary Gaussian process (Definition 7) with with mean $\mathbb{E}[y(t)]=0$ and autocorrelation function $\mathbb{E}[y(t_{1})\cdot y(t_{2})]=k_{\mu}(t_{1},t_{2})$ for any $t_{1},t_{2}$, where $k_{\mu}$ is defined in (6).

As discussed, the prior of Definition 8 amounts to a prior on the power spectral density of $y$, with the expected power spectral density given by $p_{\mu}$. Formally:

Claim 56** (Equivalent Power Spectral Density Prior).**

Consider $y$ distributed as in Definition 8. Suppose that $p_{\mu}$ is bounded. For every $T>0$, let $\hat{y}_{T}:\mathbb{R}\rightarrow\mathbb{R}$ be the truncated Fourier transform of $y$, a.k.a. the amplitude spectral density of $y$:

[TABLE]

For every $T$, $\hat{y}_{T}$ is a Gaussian process with $\mathbb{E}[\hat{y}(t)]=0$. Also as $T$ goes to infinity, the covariance of $\hat{y}_{T}$ converges to a diagonal covariance given by $p_{\mu}$. That is, for any $\xi_{1},\ldots,\xi_{s}\in\mathbb{R}$, $\lim\limits_{T\rightarrow\infty}[\hat{y}_{T}(\xi_{1}),\ldots,\hat{y}_{T}(\xi_{s})]\sim\mathcal{N}(0,\mathbf{P})$ where $\mathbf{P}$ is a diagonal matrix with $\mathbf{P}_{i,i}=p_{\mu}(\xi_{i})$.

Proof.

$\hat{y}$ is a Gaussian process since it is a linear transformation of a Gaussian process, $y$ [Ras04]. We first check that for every $T$, the mean of this random process is zero at every point $\xi$.

[TABLE]

where the application of Fubini’s theorem in second line above is valid because $k_{\mu}(0)=1$ and hence for every $t\in\mathbb{R}$, $\mathbb{E}[|y(t)|]<\infty$. Now in order to show that the covariance of $\hat{y}_{T}$ converges to being diagonal we check the covariance of $\hat{y}_{T}$ at two arbitrary points $\xi_{1}\neq\xi_{2}$, as $T$ goes to infinity,

[TABLE]

Now note that for every fixed $T$,

[TABLE]

Therefore,

[TABLE]

where the last equality above used Plancherel theorem. Now we switch the order of two integrals,

[TABLE]

Let $\tilde{\mu}$ be the measure which induces the probability density $p_{\mu}(\cdot+\xi_{2})$, hence it satisfies $d{\tilde{\mu}}(\xi)=p_{\mu}(\xi+\xi_{2})d{\xi}$. Now if we take the limit of covariance as $T\rightarrow\infty$ we get that,

[TABLE]

For ease of notation we call the integrand in above $f_{T}(\xi)=T\operatorname{sinc}{(T(\xi-\xi_{1}+\xi_{2}))}{\operatorname{sinc}{(T\xi)}}$. Remember, the assumption is that $\xi_{1}\neq\xi_{2}$. The sequence $\{f_{T}(\xi)\}$ converges pointwise to zero for all $\xi\in\mathbb{R}\setminus\{\xi_{1},\xi_{2}\}$. On points $\xi_{1},\xi_{2}$ it is also bounded by $\frac{1}{|\xi_{2}-\xi_{1}|}$. Therefore, the sequence $\{f_{T}(\xi)\}$ converges pointwise to zero $\tilde{\mu}$-almost everywhere. Also the seguence $\{f_{T}\}$ is $\tilde{\mu}$-almost dominated by an integrable function g in the sense that for all $T\geq 1$,

[TABLE]

$g$ exists since $|f_{T}(\xi)|\leq T\cdot\frac{2}{T|\xi|+1}\cdot\frac{2}{T|\xi-\xi_{1}+\xi_{2}|+1}\mathbin{\stackrel{{\scriptstyle\rm def}}{{=}}}h_{T}(\xi)$ for every $\xi\in\mathbb{R}$ and $h_{T}(\xi)$ is monotonely converging to zero for $\tilde{\mu}$-almost every $\xi$ and $h_{T}(\xi)$ is integrable for every value of $T$. Therefore by Lebesgue’s dominated convergence theorem we have,

[TABLE]

Finally note that the limit of the diagonal entries of the covariance, $\lim\limits_{T\rightarrow\infty}\mathbb{E}[|\hat{y}_{T}(\xi)|^{2}]=p_{\mu}(\xi)$ for every $\xi\in\mathbb{R}$ by the Wiener-Khintchine-Einstein Theorem [MC12]. ∎

It is well known [RW06] that for $y$ distributed as in Definition 8, the posterior distribution of $y$ given samples $t_{1},\ldots t_{s}\in[0,T]$ is also a Gaussian process. Its mean (the Bayes MMSE estimator) and its mode (the MAP estimator) coincide and are given by:

Theorem 57** (Gaussian Process Prior Signal Estimation – Finite Samples).**

Consider $y$ distributed as in Definition 8 and noise $n$ distributed as a Gaussian process covariance $\epsilon\cdot\mathbf{I}$. Given $t_{1},\ldots,t_{s}\in[0,T]$, let $\mathbf{y},\mathbf{n}\in\mathbb{R}^{s}$ be given by $\mathbf{y}(i)=y(t_{i})$ and $\mathbf{n}(t)=n(t_{i})$. Let $\mathbf{F}:\mathbb{C}^{s}\rightarrow L_{2}(\mu)$ be the operator defined by $[\mathbf{F}g](\xi)=\sum_{j=1}^{s}g(j)e^{-2\pi i\xi t_{j}}$. Both the MAP and MMSE estimates for $y$ are given by $\tilde{y}=\mathcal{F}_{\mu}^{*}\tilde{g}$ where:

[TABLE]

Proof.

Letting $\mathbf{K}=\mathbf{F}^{*}\mathbf{F}$, so $\mathbf{K}(i,j)=k_{\mu}(t_{i},t_{j})$, it is well known [RW06] that the posterior distribution of $y$ given $t_{1},\ldots,t_{s}$ is a Gaussian process with mean $\tilde{y}(t)$ given by:

[TABLE]

where $\mathbf{k}_{t}\in\mathbb{R}^{n}$ is given by $\mathbf{k}_{t}(i)=k_{\mu}(t_{i},t)$. It can be shown, analogously to the proof of Theorem 7, that $\tilde{y}=\mathcal{F}_{\mu}\tilde{g}$ where

[TABLE]

Further, since $\tilde{y}$ is the mean of the posterior distribution, it gives the Bayes MMSE estimator, and since this posterior distribution is a Gaussian process, also gives the MAP estimator. ∎

We can see that the least squares problem (10) roughly corresponds to a limit of the finite sample optimization problem of Theorem 57 as the number of samples goes to infinity. Via Theorem 3, this optimization problem can be solved approximately with $\tilde{O}(s_{\mu,\epsilon})$ samples using Algorithm 1 and the universal sampling distribution of Theorem 17. Via Claim 4 one can see that the lower bound of Section 6 (Theorem 24) extends to solving (10), even approximately, and thus our sample complexity is nearly optimal.