# Improved Reconstruction for Fourier-Sparse Signals

Yeqi Gao<EMAIL_ADDRESS>The University of Washington. Zhao Song
<EMAIL_ADDRESS>Adobe Research.  Baocheng Sun<EMAIL_ADDRESS>Weizmann Institute of Science.  Omri Weinstein<EMAIL_ADDRESS>The
Hebrew University and Columbia University. Ruizhe Zhang<EMAIL_ADDRESS>Simons Institute for the Theory of Computing.

We revisit the classical problem of Fourier-sparse signal reconstruction—a
variant of the _Set Query_ problem — which asks to efficiently reconstruct (a
subset of) a $d$-dimensional Fourier-sparse signal
($\|\widehat{x}(t)\|_{0}\leq k$), from minimum _noisy_ samples of $x(t)$ in
the time domain. We present a unified framework for this problem by developing
a theory of sparse Fourier transforms (SFT) for frequencies lying on a
_lattice_ , which can be viewed as a “semi-continuous” version of SFT in
between discrete and continuous domains. Using this framework, we obtain the
following results:

  * •

Dimension-free Fourier sparse recovery We present a sample-optimal discrete
Fourier Set-Query algorithm with $O(k^{\omega+1})$ reconstruction time in one
dimension, _independent_ of the signal’s length ($n$) and
$\ell_{\infty}$-norm. This complements the state-of-art algorithm of
[Kapralov, STOC 2017], whose reconstruction time is
$\widetilde{O}(k\log^{2}n\log R^{*})$, where
$R^{*}\approx\|\widehat{x}\|_{\infty}$ is a signal-dependent parameter, and
the algorithm is limited to low dimensions. By contrast, our algorithm works
for arbitrary $d$ dimensions, mitigating the $\exp(d)$ blowup in decoding time
to merely linear in $d$. A key component in our algorithm is fast spectral
sparsification of the Fourier basis.

  * •

High-accuracy Fourier interpolation In one dimension, we design a poly-time
$(3+\sqrt{2}+\varepsilon)$-approximation algorithm for continuous Fourier
interpolation. This bypasses a barrier of all previous algorithms [Price and
Song, FOCS 2015, Chen, Kane, Price and Song, FOCS 2016], which only achieve
$c>100$ approximation for this basic problem. Our main contribution is a new
analytic tool for hierarchical frequency decomposition based on _noise
cancellation_.

###### Contents

  1. 1 Introduction
    1. 1.1 Our Results
  2. 2 Technical Overview
    1. 2.1 A general framework for Fourier set query-type problems
    2. 2.2 Our techniques for signal estimation algorithms
      1. 2.2.1 Sample-optimal and high-accuracy algorithms
      2. 2.2.2 Speed up randomized spectral sparsification
    3. 2.3 Our techniques for discrete Fourier set query
    4. 2.4 High-accuracy Fourier interpolation
      1. 2.4.1 $(9+\varepsilon)$-approximation algorithm
      2. 2.4.2 $(3+\sqrt{2}+\varepsilon)$-approximation algorithm
  3. 3 Definitions of Semi-Continuous Fourier Set Query and Interpolation
    1. 3.1 Formal definitions of Fourier set query
    2. 3.2 Formal definitions of semi-continuous Fourier interpolation
  4. 4 Preliminaries
    1. 4.1 Tools and inequalities
    2. 4.2 Basics of Fourier transformation
    3. 4.3 Facts about Lattices
    4. 4.4 Facts about importance sampling
  5. 5 Energy Bounds for Fourier Signals
    1. 5.1 Energy bound for one-dimensional signals
    2. 5.2 Energy bound for high-dimensional signals
    3. 5.3 Energy bound for discrete Fourier signals
    4. 5.4 Energy bounds imply concentrations
      1. 5.4.1 Continuous case
      2. 5.4.2 Discrete case
  6. 6 Oblivious Sketching Fourier Sparse Signals
    1. 6.1 A unified approach via uniform sampling
    2. 6.2 Weighted oblivious sketching one-dimensional signals
  7. 7 Fast Implementation of Well-Balanced Sampling Procedure
    1. 7.1 Randomized BSS implies a WBSP
    2. 7.2 Fast implementation of WBSP
    3. 7.3 Trade-off between preprocessing and query
  8. 8 Sketch Distillation for Fourier Sparse Signals
    1. 8.1 Sketch distillation for one-dimensional signals
      1. 8.1.1 Sharper bound for the energy of orthogonal part of noise
    2. 8.2 Sketch distillation for high-dimensional signals
    3. 8.3 Sketch distillation for discrete signals
  9. 9 One-dimensional Signal Estimation
    1. 9.1 Sample-optimal reduction
    2. 9.2 High-accuracy reduction
  10. 10 High-dimensional Signal Estimation
    1. 10.1 Sample-optimal reduction
    2. 10.2 Bounding the sparsity
    3. 10.3 High-accuracy reduction
  11. 11 Discrete Fourier Set Query in One Dimension
    1. 11.1 Sample-optimal set query algorithm
    2. 11.2 Composition of two WBSPs
  12. 12 High-Accuracy Fourier Interpolation Algorithm
    1. 12.1 Technical tools i@: Fourier-polynomial equivalence
    2. 12.2 Technical tools ii@: filter functions
    3. 12.3 High sensitivity frequency estimation
    4. 12.4 $(9+\varepsilon)$-approximate Fourier interpolation algorithm
    5. 12.5 Sharper error control by signal-noise cancellation effect
    6. 12.6 Technical tools iii@: HashToBins
    7. 12.7 High signal-to-noise ratio (SNR) band approximation
    8. 12.8 Ultra-high sensitivity frequency estimation
    9. 12.9 High SNR and Recoverable signals
    10. 12.10 $(3+\sqrt{2}+\varepsilon)$-approximate algorithm
  13. A High Dimensional Reduction Under Noiseless Assumption
    1. A.1 Fourier basis is linear independent on randomly sampled points
    2. A.2 Reduction
  14. B Semi-continuous Approximation
    1. B.1 Properties related to Gaussians
    2. B.2 Continuous Fourier transform
      1. B.2.1 Bounding the tails
      2. B.2.2 Bounding the convolution
    3. B.3 Semi-continuouse approximation of Fourier-sparse signals
    4. B.4 Fast optimal-sparsity Fourier sparse recovery
    5. B.5 Semi-continuous approximation with a constant frequency gap
  15. C Improving Fourier Interpolation Precision in a Smaller Range
    1. C.1 Control noise
    2. C.2 $(\sqrt{2}+\varepsilon)$-approximation ratio

##  1 Introduction

The fast Fourier transform (FFT) [CT65] is a fundamental tool in engineering,
signal processing, mathematics, and theoretical computer science, with
profound applications in theory and practice. Over the years, many variations
of FFT have been studied and developed, depending on the underlying domain and
time-invariance properties of the signal [OWN+97, Osg02, Opp11]. In this
paper, we study the _sparse Fourier transform_ (SFT), where the signal is
either discrete or continuous in time domain, but $k$-_sparse_ in the
frequency domain, i.e., $\widehat{x}$ is a discrete set of size $k$:

$\displaystyle x(t)=$
$\displaystyle~{}\sum_{j=1}^{k}v_{j}e^{2\pi\mathbf{i}\langle f_{j},t\rangle}.$

Band-limited (i.e., Fourier-sparse) signals arise in many real-world datasets
and applications, from image compression and analysis [Wat94], to compressed
sensing [Don06] and (deep) learning with frequency-invariant kernels [MMM21];
for a broader exposition of SFT and its applications, we refer the reader to
the survey [GIIS14].

A prototypical problem in this setting is _band-limited signal interpolation_
[CKPS16] and the related111The Fourier interpolation literature typically
focuses on frequency estimation followed by magnitude estimation. The
magnitude estimation can be formulated as a Set-Query problem, where we are
given a set of locations and we only try to recover the Fourier coefficients
$\widehat{x}$ at the given locations. In our setup, frequencies are assumed to
lie on a lattice, hence the problems are essentially equivalent, see Section
3. _Fourier Set-Query_ problem [Pri11], which ask how to reconstruct (a subset
of) the signal from few (ideally $\sim k$) _noisy_ samples of $x(t)$ in a time
domain $[0,T]^{d}$. In this model, the algorithm has access to samples
$x(t)+g(t)$, where the signal-to-noise ratio is guaranteed to be above a
certain constant threshold (e.g., $\|x\|_{T}\gtrsim\|g\|_{T}$, where
$\|x\|_{T}^{2}:=T^{-d}\int_{[0,T]^{d}}|x(t)|^{2}\mathrm{d}t$ is the _energy_
of the signal). In general, sparse-recovery problems have two computational
aspects: the sample complexity, i.e., number of (noisy) samples required by
the algorithm, and the reconstruction time (decoding the signal from the
measurements). This problem has a long history in signal-processing and TCS
[CT65, Rey89, ASSN08, Voe11, HIKP12a, HIKP12b, GHI+13, IK14, IKP14, Boa15,
Kap16, Kap17, KVZ19, NSW19, JLS23, SSWZ23]. A fundamental fact, pointed out in
[Moi15], is that when the frequency gap is small ($\eta:=\min_{i\neq
j\in[k]}|f_{i}-f_{j}|<1/T$), exact recovery of the signal is informational-
theoretically impossible. Complementing this negative result, [PS15] gave a
$k\cdot\mathrm{polylog}(k,FT/\delta)$-time $\delta$-error reconstruction
algorithm for one-dimensional signals where $F$ is the band-limit, assuming
the time domain satisfies $T>\Omega(\log^{2}(k/\delta)/\eta)$, and that the
frequency gap $\eta$ is known. [CKPS16] strengthened this result by showing
that even if the frequency gap is unknown, _approximate_ reconstruction of
one-dimensional signals in $\mathrm{poly}(k,\log(FT))$-samples and time is
possible, in the sense that the output signal is close to the original signal
in the time domain albeit with worse sparsity in the frequency domain222More
precisely, the error guarantee is $\|y(t)-x^{*}(t)\|_{T}\leq
O(\|g(t)\|_{T}+\delta\|x^{*}(t)\|_{T})$, where $x^{*}(t)$ is the original
signal, $y(t)$ is the reconstructed signal, and $g(t)$ is the noise
distribution.. Subsequent works [CP19b, CP19a, LLM21] have improved this
result, both in sample-complexity and decoding time. Recently, [LLM21]
improved the sparsity of the output signal from $\mathrm{poly}(k)$ to
$k\mathrm{poly}\log(k)$, settling for a somewhat weaker notion of
approximation333$\|y(t)-x^{*}(t)\|_{(1-c)T}\leq\mathrm{poly}(\log(k/c\delta))\cdot\|g(t)\|_{T}+\delta\|x^{*}\|_{T}$.
than that of [CKPS16].

In the discrete setting, [Pri11] defined and studied the Set Query problem in
the standard compressed-sensing model [Don06], where the design of the sensing
matrix is unrestricted. [Kap17] defined and studied the Set Query problem in
Fourier domain, where the sensing matrix is applied to the _Fourier transform_
of $x$ (i.e., measurements are $S\widehat{x}=S\cdot\mathsf{FFT}\cdot x$). As
such, Fourier Set Query is a more challenging problem than the former one. The
current best discrete Fourier set query algorithm in 1D is due to Kapralov
[Kap17], who gave an algorithm with near-linear sample complexity
($O(k/\varepsilon)$) and decoding time
$\widetilde{O}(\varepsilon^{-1}k\log^{2.001}n\log R^{*})$, where $n$ is the
length of signal and $R^{*}$ is (roughly) the $\ell_{\infty}$ norm of the
signal in the time domain.

Unfortunately, much less was known in higher dimensions – The “curse of
dimensionality” of the Filter function [Kap17] drastically deteriorates the
sample complexity, which grows _exponentially_ with the dimension, hence
filter-based algorithms (a-la [Kap17]) are near-optimal only for small $d$.
Indeed, this drawback was one of the principal motivations of this work.

Given the abundance of largely incomparable results and diversity of
techniques mentioned above, one might wonder whether there is a systematic,
unified framework for analyzing band-limited signal reconstruction. We propose
such a framework, which _decouples_ the band-limited signal interpolation
problem into two sub-problems:

  1. 1.

Frequency Estimation:  Find $L$ net-frequencies
$\\{f_{i}^{\prime}\\}_{i\in[L]}$ such that each $f_{i}$ is close to some
$f^{\prime}_{j}$ for $j\in[L]$.

  2. 2.

Signal Estimation:  Based on the net-frequencies, approximate the original
signal.

We note that almost all the previous works [PS15, CKPS16, LLM21] fall into
this framework, yet differ in the techniques used to solve these two sub-
problems. Our starting point is the observation that the assumption of an
$\eta$-frequency gap of the signal, is roughly equivalent to assuming that the
frequencies lie on the grid $\eta\cdot\mathbb{Z}^{d}$. It is therefore natural
to generalize this assumption by considering the case where signal frequencies
lie on a _$d$ -dimensional lattice_ ${\cal L}=\Lambda({\cal B})$, where ${\cal
B}$ is a basis of ${\cal L}$ and is given. We interchangeably call this
problem _lattice Fourier interpolation_ , or _semi-continuous signal
reconstruction_ , as lattice frequencies can be viewed as interpolating
between discrete and continuous domains. Solving this problem is our key tool
en-route to faster sparse-recovery, but it is also interesting on its own
right.

In this framework, the Fourier interpolation problem can be reduced to
sequentially solve frequency estimation and signal estimation problems. We
focus on the second step, that is, given the result of frequency estimation
for a semi-continuous signal, how to efficiently reconstruct the coefficients.
It can be formulated as a semi-continuous Fourier set-query problem, where the
queried frequency set $L$ (a set of lattice points) is obtained from Frequency
Estimation, and our goal is to recover the coefficients of these frequencies.
The generality of our framework allows us to apply it to the _discrete_
Fourier Set-Query problem as well, where the goal is to (approximately)
compute DFT only on a subset $S\subset[n]$ of $k$ coordinates. (See Remark 1.3
for more detailed discussions of the relationship between semi-continuous
signal estimation and discrete Fourier set-query.)

###  1.1 Our Results

##### Semi-continuous signal estimation

Our first main result shows that for semi-continuous signals (whose
frequencies lie on a lattice), Signal Estimation can be efficiently reduced to
Frequency Estimation:

######  Theorem 1.1 (1D Semi-continuous Signal Estimation, informal version of
Theorems 9.1 and 9.4).

Let $\Lambda(\mathcal{B})\subset\mathbb{R}$ denote the lattice
$\Lambda(\mathcal{B})=\\{z\in\mathbb{R}|z=c\eta,c\in\mathbb{Z}\\}$, and
suppose $x^{*}(t)=\sum_{j=1}^{k}v_{j}e^{2\pi\mathbf{i}f_{j}t}$ with
$f_{j}\in\Lambda({\cal B})$. Given observations of the form
$x(t)=x^{*}(t)+g(t)$ for $t\in[0,T]$ for arbitrary noise distribution $g(t)$,
let $L$ be the result of Frequency Estimation for $x(t)$ such that for each
$f_{i}$, there is an $f^{\prime}_{i}\in L$ with $|f_{i}-f^{\prime}_{i}|\leq
D/T$. Then given the set $L$ as input, we can obtain:

  * •

A Sample-optimal algorithm: There is an algorithm that takes
$O(\widetilde{k})$ samples and outputs a $\widetilde{k}$-sparse signal $y(t)$
in $O(\widetilde{k}^{\omega+1})$-time444$\omega\approx 2.373$ is the fast
matrix-multiplication exponent [Wil12, AW21], such that
$\|y-x^{*}\|_{T}^{2}\leq O(\|g\|_{T}^{2})$ holds with high probability, where
$\widetilde{k}:=O(|L|(1+D/(T\eta)))$ is the output sparsity.

  * •

A High-accuracy algorithm:  For any $\varepsilon\in(0,1)$, there is an
algorithm that takes $\widetilde{O}(\varepsilon^{-1}\widetilde{k})$ samples
and outputs a $\widetilde{k}$-sparse signal $y(t)$ in
$O(\varepsilon^{-1}\widetilde{k}^{\omega})$-time such that
$\|y-x^{*}\|_{T}^{2}\leq(1+\varepsilon)\|g\|_{T}^{2}$ holds with high
probability.

We note that our algorithms can be easily adapted to high-dimensional signal
estimation.

##### Discrete Fourier set query

Our next result provides an efficient, high-accuracy algorithm for the
discrete Fourier Set Query problem in any dimension.

######  Theorem 1.2 (Discrete Fourier Set Query, informal version of Theorems
11.1).

For any $d\geq 1$, let $n=p^{d}$ for some $p\in\mathbb{N}_{+}$. Given a vector
$x\in(\mathbb{C}^{p})^{\otimes d}$, for $k\geq 1$, any $S\subseteq[n]$,
$|S|=k$, there exists an algorithm that takes $O(\varepsilon^{-1}k)$ samples
from $x$, runs in
$O(\varepsilon^{-1}k^{\omega+1}+\varepsilon^{-1}k^{\omega-1}d)$ time, and
outputs a vector $x^{\prime}\in\mathbb{C}^{n}$ such that

$\displaystyle\|(x^{\prime}-\widehat{x})_{S}\|_{2}^{2}\leq\varepsilon\|\widehat{x}_{[n]\backslash
S}\|_{2}^{2}$

holds with probability at least $0.9$. Note that $\widehat{x}$ is the
$d$-dimensional discrete Fourier transform of $x$,
$\widehat{x}_{f}=\sum_{t\in[p]^{d}}x_{t}e^{-2\pi\mathbf{i}\langle
f,t\rangle/p},f\in[p]^{d}$.

######  Remark 1.3.

Discrete Fourier set query can be viewed as a special case of semi-continuous
signal estimation, where the queried frequencies (supported on $S$) lie on the
integer lattice $\mathbb{Z}^{d}$. And the remaining part of the signal with
frequencies outside of $S$ (i.e., $\widehat{x}_{[n]\backslash S}$) corresponds
to the noise $g$ in the semi-continuous signal estimation problem. One
difference between these two problems is that for semi-continuous signal
estimation, we assume the signal has a continuous time domain $[0,T]$, while
for discrete Fourier set query, the signal has a discrete time domain
$[p]^{d}$.

######  Remark 1.4.

In one-dimension ($d=1$), the runtime of our set query algorithm can be
simplified to $O(\varepsilon^{-1}k^{\omega+1})$. And prior to this work, the
only known result for $\ell_{2}/\ell_{2}$ Fourier set query (due to [Kap17])
achieves the same sample complexity but runs in time
$O(\varepsilon^{-1}k\log^{2.001}n\log R^{*})$. Here, $R^{*}$ is an upper bound
on the $\|\cdot\|_{\infty}$ norm of the vector, typically assumed to be
$\mathrm{poly}(n)$. We emphasize that our approach immediately yields a
$\mathrm{poly}(k)$-time algorithm, independent of $\log(n)$ and $R^{*}$, and
thus generalizes to arbitrary dimension (see next result). This, however,
comes at a price of a slower dependence on $k$, as opposed to the linear
dependence of [Kap17].

In high dimensions ($d>1$), the decoding time of [Kap17] scales as
$\log^{d}(n)$ due to the curse of dimensionality of the Filter function
used555See [NSW19] for more detailed discussions about the time complexity in
high dimensions., while our algorithm scales linearly with $d$.

##### High-accuracy Fourier interpolation

All previous algorithms for the Fourier interpolation problem provide rather
large approximation guarantees on the reconstruction error

$\|y-x^{*}\|_{T}\leq C\|g\|_{T}+\delta\|x^{*}\|_{T},$

where $C$ is an absolute constant around 100 [CKPS16], and is never better
than $3$ due to repeated applications of the _triangle inequality_ for the
errors from different sources. We develop a _sharper noise control_ technique,
together with a high sensitivity frequency estimation method and an efficient
signal estimation algorithm, which allows us to break this barrier and obtain
a $C=(1+\sqrt{2}+\varepsilon)$-approximation in
$\mathrm{poly}(k,1/\delta,1/\varepsilon)$ time. Our high-accuracy Fourier
interpolation algorithm relies on a new error analysis.

######  Theorem 1.5 (High-accuracy Fourier interpolation, informal version of
Theorem 12.41).

Let $x^{*}(t)$ be a $k$-Fourier sparse signal with frequencies in $[-F,F]$.
Given observations $x(t)=x^{*}(t)+g(t)$ in time duration $[0,T]$, where $g$ is
arbitrary noise. For $\varepsilon,\delta\in(0,1)$, there exists an algorithm
that uses $\mathrm{poly}(k,\varepsilon^{-1},\log(1/\delta))\log(FT)$ samples
and runtime, and outputs a
$\mathrm{poly}(k,\varepsilon^{-1},\log(1/\delta))$-Fourier-sparse signal
$y(t)$ such that

$\displaystyle\|y-x^{*}\|_{T}\leq(3+\sqrt{2}+\varepsilon)\|g\|_{T}+\delta\|x^{*}\|_{T}.$

##  2 Technical Overview

The following section contains a streamlined technical overview of the main
ideas and techniques required to prove our results in Section 1.1. Section 2.1
develops a unified framework for a wide range of Fourier set query-type
problems. Section 2.2 focuses on implementing the framework for semi-
continuous signals estimation and proving Theorem 1.1. Section 2.3 focuses on
implementing the framework for discrete Fourier set query and proving Theorem
1.2. Section 2.4 shows how to apply our signal estimation algorithm and obtain
a high-accuracy Fourier interpolation algorithm as in Theorem 1.5.

###  2.1 A general framework for Fourier set query-type problems

Many Fourier-related problems can be expressed as the following set-query
problem: Given a signal $x(t)$ (either continuous or discrete), we wish to
recover the Fourier coefficients $\widehat{x}$ only at the coordinates of a
predetermined set $S$. The observed signal $x(t)$ can be accessed only through
samples, either noiseless or noisy. In the latter case, one is only allowed to
access $x(t)+g(t)$ for an arbitrary function $g(t)$.

A natural approach to such problem is to use _linear regression_ – Notice that
the observed signal $x(t)$ can be decomposed into
$x_{S}(t)+(x_{\overline{S}}(t)+g(t))$, where $x_{S}(t)$ is a portion of the
noiseless signal $x(t)$ with frequencies in $S$, and $x_{\overline{S}}(t)$ is
the remaining part such that $(x_{\overline{S}}+g)(t)$ together can be treated
as noise. In this terminology, recovering $\widehat{x}_{S}$ is equivalent to
solving a linear regression problem in the subspace spanned by the _Fourier
basis_ functions with frequencies in $S$. Implementing this approach, however,
has two substantial challenges:

  1. 1.

Sample complexity: How should one select the sample points
$\\{t_{1},\dots,t_{s}\\}$ in the time domain to solve the linear regression?
Recall we wish to use as few samples as possible.

  2. 2.

Estimation accuracy:  How can one guarantee that the solution of the linear
regression will not be corrupted by the noise? In other words, how can one
ensure the recovered signal is close to the relevant projection $x_{S}(t)$ of
the true signal?

To resolve these two issues, our framework consists of the following three
steps, and we note that this framework can be applied to both continuous and
discrete signals.

##### Step 1: Oblivious sketching

A straightforward approach for choosing the sample points for the linear
regression is to uniformly666In the one-dimensional case, it is possible to
design a more clever _nonuniform_ oblivious distribution that achieves near-
linear sample complexity [CP19a]. However, in $d\geq 2$ dimensions, the best-
known construction is a uniform oblivious sketch with polynomially many
samples; more on this in Section 2.2. sample points in the time domain. Let
$S_{0}$ denote the set of i.i.d uniform samples. Then, we need to know how
many samples are sufficient to guarantee that $S_{0}$ is a good sketch for the
signal, i.e., $\|x_{S}\|_{S_{0}}\approx\|x_{S}\|_{T}$. The size of $S_{0}$ can
be bounded by an important quantity— _the Fourier energy bound_ , which shows
how far the maximum magnitude of a signal in a family ${\cal F}$ can deviate
from its energy in a given time duration. More specifically, the energy bound
$R$ of the family ${\cal F}$ is defined as:

$\displaystyle R:=\sup_{f\in{\cal
F}}\sup_{t\in[0,T]}\frac{|f(t)|^{2}}{\|f\|_{T}^{2}},$

In [BE06, Kós08, CKPS16, CP19b], energy bounds for one-dimensional continuous
Fourier-sparse signals are proved. We further show (nearly) tight energy
bounds for high-dimensional discrete and continuous signals. Although the size
of $S_{0}$ obtained in this way may be quite large, an appealing property is
that it only depends on the Fourier sparsity instead of the actual signal
$x(t)$. In this sense, the oblivious sketching step acts as a “preconditioner"
for our regression problem.

##### Step 2: Sketch distillation

We will reduce the size of the oblivious sketching by sub-sampling a _linear-
sized_ subset $S_{1}$ from $S_{0}$, in a data-dependent fashion, such that
$S_{1}$ is still a good sketch of $x_{S}$. More specifically, suppose the
frequencies of $x_{S}$ are given. We can use a _well-balanced sampling
procedure_ defined by [CP19a] to sample a set of $s=O(k)$ points such that the
$i$-th point $t_{i}$ is sampled from a distribution $D_{i}$ supported by
$S_{0}$, and a coefficient vector $\alpha\in\mathbb{R}^{s}$. Then, this sub-
sampler satisfies the following properties:

  * •

For the weights $w_{i}:=\alpha_{i}\frac{1/|T|}{D_{i}(t_{i})}$ for $i\in[s]$,
we have $\|x\|_{S_{1},w}\approx\|x\|_{T}$ for all signals $x\in{\cal F}$,
where $\|x\|_{S_{1},w}:=\sum_{t\in S_{1}}w_{t}\cdot|x(t)|^{2}$ and ${\cal F}$
is the signal family spanned by the frequencies in $S$.

  * •

The sum of coefficients $\alpha_{i}$ is small and each distribution $D_{i}$ is
not “ill-conditioned”777Formal definition is in Definition 7.1..

Figure 1: An illustration of sketching the set query signals. The green curve
is the queried signal $x_{S}$, the red curve is the signal with remaining
frequencies $x_{\overline{S}}$, and the blue curve is the noise $g$. Suppose
$S_{1}=\\{t_{1},\dots,t_{5}\\}$ are the sampling points. Then, we need to
guarantee that $\|x_{S}\|_{S_{1},w}\approx\|x_{S}\|_{T}$ as well as
$\|x_{\overline{S}}+g\|_{S_{1},w}\lesssim\|x_{\overline{S}}+g\|_{T}$.

On the one hand, the first property guarantees that the sketch distillation
outputs a sample set $S_{1}$ of linear size and
$\|x_{S}\|_{S_{1},w}\approx\|x_{S}\|_{T}$. On the other hand, the sketch
distillation process is _robust to noise_. An easy information-theoretic
argument shows that the estimation error must be proportional to the energy of
the noise, i.e., $\|x_{\overline{S}}+g\|_{T}$. However, after the sketch
distillation the weighted energy of the noise
$\|x_{\overline{S}}+g\|_{S_{1},w}$ could be amplified (i.e.,
$\|x_{\overline{S}}+g\|_{S_{1},w}\gg\|x_{\overline{S}}+g\|_{T}$), resulting in
a large estimation error. Fortunately, the second property ensures that it
will not happen and in expectation
$\|x_{\overline{S}}+g\|_{S_{1},w}\lesssim\|x_{\overline{S}}+g\|_{T}$.

##### Step 3: Weighted linear regression

The last step is to solve a (weighted) linear regression with the noisy
samples $\\{\widetilde{x}(t_{i})\\}_{i\in[s]}$ and weight
$w\in\mathbb{R}^{s}_{>0}$. For simplicity, take the one-dimensional continuous
$k$-Fourier-sparse signal as an example. Suppose the set-query frequencies of
$x_{S}$ are $S=\\{f_{1},\dots,f_{k}\\}$. Then, we consider the following
weighted linear regression:

$\displaystyle\underset{v^{\prime}\in\mathbb{C}^{k}}{\min}~{}\left\|\sqrt{w}\circ(Av^{\prime}-b)\right\|_{2},$
(1)

where $\sqrt{w}:=(\sqrt{w_{1}},\dots,\sqrt{w_{s}})$, the coefficients matrix
$A\in\mathbb{C}^{s\times k}$ and the target vector $b\in\mathbb{C}^{s}$ are
defined as follows:

$\displaystyle A_{i,j}:=$
$\displaystyle~{}\exp(2\pi\mathbf{i}f_{j}t_{i})~{}~{}~{}\forall(i,j)\in[s]\times[k],~{}~{}\text{and}$
$\displaystyle b_{i}:=$ $\displaystyle~{}\widetilde{x}(t_{i})~{}~{}~{}\forall
i\in[s].$

Let $v^{\prime}$ denote an optimal solution of Eq. (1). Then, we output a
signal $y(t):=\sum_{i=1}^{k}v^{\prime}_{i}\exp(2\pi\mathbf{i}f_{i}t)$.

Using the error analysis in sketch distillation (Step 2), we can upper-bound
the estimation error by $\|y-x_{S}\|_{T}\lesssim\|x_{\overline{S}}+g\|_{T}$.

###  2.2 Our techniques for signal estimation algorithms

In this section, we show how to instantiate the framework and get Signal
Estimation algorithms in the semi-continuous setting (Theorem 1.1).

We first note that given the output of Frequency Estimation algorithm, a set
$L$ such that all the true frequencies are close to $L$, the Signal Estimation
problem can be formulated as a Fourier set-query problem. The idea is that by
the guarantee of $L$ and the semi-continuous assumption, there will not be too
many lattice points that are close to $L$, and we can efficiently find all of
them. Let’s denote this set of candidate frequencies by $\widetilde{S}$ of
size $\widetilde{k}$. Then, Signal Estimation problem is reduced to a set-
query for $\widehat{x^{*}}_{\widetilde{S}}$. In Section 2.2.1, we discuss how
to apply our set-query framework to obtain sample-optimal and high-accuracy
algorithms for 1-D signals and how to generalize to higher dimensions. In
Section 2.2.2, we show how to implement Sketch Distillation very efficiently
by speeding up the randomized spectral sparsification, which may be of
independent interest.

####  2.2.1 Sample-optimal and high-accuracy algorithms

Based on the aforementioned reduction, we briefly show the instantiation of
our framework for the Signal Estimation problem. For convenience, suppose
$\widetilde{k}=O(k)$ throughout this part.

##### Sample-optimal signal estimation

For one-dimensional $k$-Fourier-sparse signals, [Kós08] proved that the energy
bound is $O(k^{2})$, which implies that the uniform sketching for the signal
needs at least $\Omega(k^{3}\log k)$ samples. Although the sketch size can be
reduced to $O(k)$ via Sketch Distillation, the algorithm has to first sample
$\widetilde{O}(k^{3})$ points in $[0,T]$, which already takes $O(k^{3})$-time.
It is possible to improve this straightforward approach using _weighted
oblivious sketching_ [CP19a], namely, we sample a set $S_{0}$ of points in
$[0,T]$ from a carefully-chosen non-uniform distribution $D$, and assign each
point a weight (this works only in 1D, see last section). Using the following
distribution constructed by [CP19a]:

$\displaystyle D(t):=\begin{cases}{c}/(1-|t/T|),&\text{ for }|t|\leq
T(1-{1}/k)\\\ c\cdot k,&\text{ for }|t|\in[T(1-{1}/k),T]\end{cases}$

we only need to take $|S_{0}|=O(k\log k)$ samples to guarantee that
$\|x^{*}\|_{S_{0},w_{0}}\approx\|x^{*}\|_{T}$ holds with high probability. In
this way, oblivious sketching in Step 1 will not be a time-consuming step. In
Step 2, the sketch $S_{0}$ will be distilled to a subset $S_{1}$ of size
$O(k)$. Finally, in Step 3, we sample the signal at the points in $S_{1}$ and
solve the weighted linear regression to recover a $k$-sparse signal $y(t)$.
This gives us a _linear-sample_ reduction from Frequency Estimation to Signal
Estimation with $O(1)$-estimation error with high probability.

##### High-accuracy signal estimation

For one-dimensional semi-continuous signals, we also discover a “sample-
accuracy trade-off”. If we can use nearly-linear (i.e., $\widetilde{O}(k)$)
samples, then we can skip Step 2 and directly use $(S_{0},w_{0})$ to solve the
linear regression in Step 3.

The advantage of this approach is that the sampling procedure for
$(S_{0},w_{0})$ is well-balanced, and a sharper error analysis shows that we
can achieve much smaller errors. The main observation is that the noise can be
decomposed into $g^{\parallel}\in{\cal F}$ and $g^{\bot}$ orthogonal to ${\cal
F}$. Since we will solve a linear regression (in Step 3) in the space ${\cal
F}$, the contribution of $g^{\parallel}$ to the final estimation error will
not blow up. For the orthogonal part, we find an orthonormal basis
$\\{u_{1},\dots,u_{k}\\}$ for ${\cal F}$ and look at the weighted projection
$\langle g^{\bot},u_{i}\rangle_{S_{1},w}:=\sum_{t\in
S_{1}}w_{t}\overline{u_{i}(t)}g^{\bot}(t)$, whose magnitude indicates how much
the noise is amplified due to sketching in each direction. The well-balanced
sampling procedure gives that the total magnitudes $\sum_{i=1}^{k}|\langle
g^{\bot},u_{i}\rangle_{S_{1},w}|^{2}\leq\varepsilon\cdot\|g^{\bot}\|_{T}^{2}$
holds with high probability, which will imply an
$\varepsilon\cdot\|x_{\overline{S}}+g\|_{T}^{2}$-estimation error for any
small $\varepsilon$.

(a) The function space with norm $\|\cdot\|_{T}$ (b) The function space with
norm $\|\cdot\|_{S_{0},w_{0}}$

(c) The function space with norm $\|\cdot\|_{S_{1},w}$

Figure 2:  The composition of two WBSPs may not be a WBSP. In (a), the central
point is the ground-truth signal $x^{*}(t)$, and the yellow vector $y(t)$ is
the final output of the composited WBSPs. Our goal is to bound the estimation
error $\|y(t)-x^{*}(t)\|_{T}$. Let $y_{S_{0},w_{0}}(t)\in{\cal F}$ be the
optimal approximation of the observing signal in the basis $(S_{0},w_{0})$,
which is obtained by Step 1 of the framework. In (b), we show that by the
property of the first WBSP, the estimation error
$\|y_{S_{0},w_{0}}(t)-x^{*}(t)\|_{T}$ can be bounded by the projection of
$g^{\bot}$ in the space $({\cal F},\|\cdot\|_{S_{0},w_{0}})$, which is about
$\varepsilon\|g^{\bot}(t)\|_{T}$. In (c), the space is reduced to $({\cal
F},\|\cdot\|_{S_{1},w_{1}})$ by Step 2 of the framework, and the estimation
error $\|y_{S_{0},w_{0}}-y\|_{S_{1},w_{1}}$ is bounded by the projection of
$(x^{*}(t)+g^{\bot}-y_{S_{0},w_{0}})$, which is about
$\varepsilon\|x^{*}(t)+g^{\bot}-y_{S_{0},w_{0}}(t)\|_{S_{0},w_{0}}$. However,
these two error bounds $\varepsilon\|g^{\bot}(t)\|_{T}$ and
$\varepsilon\|x^{*}(t)+g^{\bot}-y_{S_{0},w_{0}}(t)\|_{S_{0},w_{0}}$ could not
imply any bound on $\|y(t)-x^{*}(t)\|_{T}$.

Then, in step 3, we can also decompose the final estimation error
$\|y-x^{*}\|_{T}^{2}$ into two parts: one contributed by $g^{\parallel}$ with
energy $\|g^{\parallel}\|_{T}^{2}$, and another contributed by $g^{\bot}$ with
energy at most $(1+\varepsilon)\|g^{\bot}\|_{T}^{2}$. Combining the two parts
together, we achieve a a high-accuracy guarantee:

$\displaystyle\|y-x^{*}\|_{T}^{2}\leq(1+\varepsilon)\|g\|_{T}^{2}.$

We remark that this error analysis does not work for the sample-optimal
algorithm. Roughly speaking, the weighted sketch $(S_{1},w)$ in the sample-
optimal algorithm is generated by a two-step sampling procedure: one in the
oblivious sketching (Step 1) and another in the sketch distillation (Step 2).
Even if both of them are well-balanced, when they are composited together and
considered as a single sampling procedure, they may not be well-balanced.

##### High-dimensional signal estimation

For high dimensional signals, we prove a $k^{O(d)}$-energy bound for
$k$-Fourier-sparse signals in $d$-dimensions in Section 5.2, where we also
provide a nearly-matching lower bound. We view it as a new generic tool in
Fourier analysis and are quite a tour-de-force. Thus, we basically follow the
three-step framework. One tricky thing in high-dimension is how to bound the
output signal’s sparsity, which is equivalent to the number of $d$-dimensional
lattice points close to a fixed set $L$. We further reduce it to a clean math
problem about lattices: let $\Lambda({\cal B})$ be a lattice in
$\mathbb{R}^{d}$. For $r>0$, how many lattice points can be within an
$r$-radius ball, i.e., $\sup_{x\in\mathbb{R}^{d}}|\Lambda({\cal B})\cap
B_{d}(x,r)|$. We upper-bound this quantity via different approaches, which
might be of independent interest. More details are deferred to the appendix.

####  2.2.2 Speed up randomized spectral sparsification

In Section 2.2.1, we do not discuss how to implement a well-balanced sampling
procedure. [CP19a] proved that _Randomized BSS_ algorithm in [BSS12, LS15]
yields a well-balanced sampling procedure. However, this algorithm is slow
when the sampling domain $S_{0}$ is very large. One contribution of this work
is to improve the time and space costs of the randomized BSS algorithm.

We observe that the bottleneck of each iteration in the original Randomized
BSS algorithm is to sample a point $t\in S_{0}$ from the distribution $D_{j}$
defined by $D_{j}(t)=v(t)^{\top}E_{j}v(t)$, where $v(t)$ is a $k$-dimensional
vector and $E_{j}$ is a $k$-by-$k$ positive semi-definite matrix determined by
the potential function value at the $j$-th iteration. Suppose $E_{j}$ and
$\\{v(t)\\}_{t\in S_{0}}$ have already been computed. A naive approach to
sampling from $D_{j}$ is to compute the probability $D_{j}(t)$ for each $t\in
S_{0}$, which takes $O(nk^{2})$-time per iteration, where $n=|S_{0}|$ is the
size of the sampling domain. To improve the algorithm, we consider a more
general data structure problem— _Online Quadratic-Form Sampling_. In this
problem, we are given $n$ vectors $v_{1},\dots,v_{n}\in\mathbb{R}^{k}$. In
each query, the input is a positive semi-definite matrix
$A\in\mathbb{R}^{k\times k}$ and we need to sample an $i\in[n]$ with
probability proportional to $v_{i}^{\top}Av_{i}$. (See Problem 7.4 for formal
definition). The naïve space and query time for generating a sample is
$O(nk^{2})$. We design two data structures with substantially faster time-
space tradeoff:

######  Theorem 2.1 (Online Quadratic-Form Sampling, informal version of
Theorems 7.5 and 7.8).

The Online Quadratic-Form Sampling problem admits the following two data
structures:

  * •

Data Structure 1: $O(nk^{2})$ preprocessing time of the vectors
$\\{v_{i}\\}_{i\in[n]}$ , $O(k^{2}\log n)$ query-time for generating a sample,
and $O(nk^{2})$-space.

  * •

Data Structure 2: $O(nk^{\omega-1})$ preprocessing time,
$O(k^{2}\log(n/k)+k^{\omega})$ query-time for generating a sample, and
$O(nk)$-space.

The main idea is to construct a range search tree for $\\{1,2,\dots,n\\}$, and
for a node corresponding to range $[l,r]$, stores the matrix
$\sum_{i=l}^{r}v_{i}v_{i}^{\top}$. Subsequently, for each query matrix $A$, we
traverse the root-to-leaf path in the tree (which corresponds to an element in
$[n]$), and output this element as a sample. The rule for descendin the tree
resembles a form of rejection sampling: At a node with range $[l,r]$, we know
that its left child has range $[l,m]$ and right child has right $[m+1,r]$,
where $m=\lfloor(l+r)/2\rfloor$. Then, we decide whether move to the left or
the right subtree by tossing a coin with probability:

$\displaystyle
p_{\mathsf{left}}:=\frac{\langle\sum_{i=l}^{m}v_{i}v_{i}^{\top},A\rangle}{\langle\sum_{i=l}^{r}v_{i}v_{i}^{\top},A\rangle},~{}~{}\text{and}~{}~{}p_{\mathsf{right}}:=1-p_{\mathsf{left}},$

where
$\langle\sum_{i=l}^{r}v_{i}v_{i}^{\top},A\rangle:=\sum_{i=l}^{r}v_{i}^{\top}Av_{i}$
is the trace-product of matrices. Therefore, the probability
$p_{\mathsf{left}}$ equals to the _conditional_ probability $\Pr_{{\bf
i}\sim{\cal D}_{A}}[{\bf i}\in[l,m]~{}|~{}{\bf i}\in[l,r]]$. By the chain rule
of conditional probability, we get that the output distribution of this
procedure is exactly equal to ${\cal D}_{A}$. This data structure can be built
in $O(nk^{2})$ and each query only takes $O(k^{2}\log n)$-time. To store the
matrices in each node, this data structure uses $O(nk^{2})$-space. We can
further improve the space complexity to $O(nk)$ by trading-off the
preprocessing time and the query time.

By plugging-in this data structure to the Randomized BSS algorithm, we can
improve the time and space complexity by a factor of $k$ when $n$ is large.

Figure 3: An example of the outer-product range tree with $n=4$. For a query
matrix $A$, the sampling probability of the red path is
$\frac{\langle\sum_{i=1}^{2}v_{i}v_{i}^{\top},A\rangle}{\langle\sum_{i=1}^{4}v_{i}v_{i}^{\top}\rangle}\cdot\frac{\langle
v_{2}v_{2}^{\top},A\rangle}{\langle\sum_{i=1}^{2}v_{i}v_{i}^{\top},A\rangle}=\frac{v_{2}^{\top}Av_{2}}{\sum_{i=1}^{4}v_{i}^{\top}Av_{i}}$.

###  2.3 Our techniques for discrete Fourier set query

For simplicity, we consider one-dimensional discrete Fourier signal $x(t)$,
which is a length-$n$ vector such that
$x(t)=\sum_{j=1}^{n}\widehat{x}_{j}e^{2\pi\mathbf{i}jt/n}$ for any $t\in[n]$.
And the set query problem asks to recover $\widehat{x}_{S}$, for a given
$k$-subset $S\subset[n]$. Let $x_{S}$ be the part of signal with frequencies
in $S$, i.e., $x_{S}(t)=\sum_{f\in S}\widehat{x}_{f}e^{2\pi\mathbf{i}ft/n}$,
which is a $k$-Fourier-sparse signal.

The high-level idea of obtaining Theorem 1.2 is as follows. We prove that the
energy bound for discrete $k$-sparse signal is $k$, which implies that
uniformly sample a set $S_{0}\subset[n]$ of $O(k\log k)$ points can form a
good sketch of the signal $x_{S}$. Then, by Sketch Distillation, we can find a
subset $S_{1}\subset S_{0}$ of linear size together with a weight vector $w$
such that $\|x_{S}\|_{S_{1},w}\approx\|x_{S}\|_{2}/n$. Finally, we can recover
$\widehat{x}_{S}$ by solving a weighted linear regression on the samples
$\\{x(t)\\}_{t\in S_{1}}$. By a direct error analysis, it is easy to see that
this algorithm can achieve $O(1)$-estimation error. In the followings, we will
show that it can perform much better.

##### Composition of well-balanced samplers

The key step to proving a $(1+\varepsilon)$-error guarantee is to show that
the final weighted sketch $(S_{1},w)$ can be generated by a well-balanced
sampling procedure. In general, compositing two well-balanced sampling
procedures may not be well-balanced. More specifically, we can show that the
composition sampler satisfies the first property of well-balanced; that is,
the output set and weight can well-approximate the energy of every function in
the family ${\cal F}$. However, the second property about the sum of the
composition sampler’s coefficients and the condition number of the composition
sampling distribution may not hold.

The above discrete set query algorithm is a very special case of composition
in the sense that the first sampler samples each element from the same,
simplest distribution—uniform distribution over $[n]$. Then, we can prove that
in the composition sampler, each sample is equivalent to directly sampled from
the uniform distribution. Furthermore, for discrete Fourier-sparse signal, we
also have a tight energy bound of $R=k$. Using these results, we can show that
the composition sampler in our algorithm is well-balanced! Then, we are able
to apply our sharper error analysis and get that

$\displaystyle\|\widehat{y}_{S}-\widehat{x}_{S}\|_{2}\leq\varepsilon\cdot\|\widehat{x}_{\overline{S}}\|_{2}^{2}$

holds with high probability. Therefore, we obtain a linear-sample and high-
accuracy algorithm for discrete Fourier set query.

######  Remark 2.2.

This result is a fundamental departure from [CP19a] since they assumed that
the noise has a zero-mean (i.e., $\mathbb{E}[g(t)]=0$), which makes the WBSPs
composition much easier. However, we do not have such an assumption. A natural
fix is to ensure the orthogonality of the noise $g(t)$ to our basis functions
$\\{\exp(-2\pi\mathbf{i}f_{i}t/n)\\}$. Unfortunately, this only holds for the
first WBSP, but not for the second one (since sampling from the weighted
output of the first WBSP can break the orthogonality).

Our main novelty is to directly find an equivalent sampling procedure to the
combination of two WBSPs (just for the analysis). Then, using the special
properties of DFT, we prove that the composition also works in our setting
(see Section 11.2 Lemma 11.4 for more details). We believe this technique will
be useful for other set-query or active learning problems.

Algorithm 1 Discrete One-Dimensional Signal Set-Query Algorithm (Informal)

1:procedure SetQuery($x$, $n$, $k$, $S$, $\varepsilon$) $\triangleright$
Theorem 1.2

2: $\\{f_{1},f_{2},\cdots,f_{k}\\}\leftarrow S$

3:/*Step 2: Oblivious Sketching*/

4: $S_{0}\leftarrow$ $O(\varepsilon^{-2}k\log(k))$ i.i.d. samples from
$\text{Uniform}([n])$

5:/*Step 3: Sketch Distillation*/

6: Let
${\mathcal{F}}=\\{\sum_{j\in[k]}v_{j}\exp(2\pi\mathbf{i}f_{j}t/n)~{}|~{}v_{j}\in\mathbb{C}\\}$.

7:
$\\{t_{1},t_{2},\cdots,t_{s}\\},w\leftarrow\textsc{RandBSS+}(k,\mathcal{F},\mathrm{Uniform}(S_{0}),(\varepsilon/4)^{2})$
$\triangleright$ Algorithm 3

8:/*Step 4: Weighted linear regression*/

9: $A_{i,j}\leftarrow\exp(2\pi\mathbf{i}f_{j}t_{i}/n)$ for each
$(i,j)\in[s]\times[k]$

10: $b_{i}\leftarrow x(t_{i})$ for each $i\in[s]$$\triangleright$ Observe the
signal at time $t_{i}$

11:
$v^{\prime}\leftarrow\underset{v^{\prime}\in\mathbb{C}^{k}}{\arg\min}\|\sqrt{w}\circ(Av^{\prime}-b)\|_{2}$

12: return $v^{\prime}$

13:end procedure

###  2.4 High-accuracy Fourier interpolation

In this section, we introduce how to apply our four-step Fourier set-query
framework together with some new techniques to obtain a high-accuracy one-
dimensional Fourier interpolation algorithm (Theorem 1.5), which improves the
constant-accuracy algorithm by [CKPS16].

Let us briefly summarize the previous algorithm in [CKPS16]. The high-level
idea is to first find some small intervals in the frequency domain such that
each contains some significant frequencies of the signal $x^{*}$. (These
intervals are called “heavy-clusters” in their paper.) Then, they use some
filter techniques (also used in [PS15]) to reduce the problem of
reconstructing $x^{*}$, a signal with multiple heavy-clusters to several
single heavy-cluster signals. Then, for each single heavy-cluster signal,
since the band-limit is small, they can efficiently estimate its frequencies.
Finally, they reconstruct a $\mathrm{poly}(k)$-sparse signal that is close to
$x^{*}$ via a robust polynomial learning algorithm. More specifically, their
algorithm consists of the following steps:

  1. 1.

They show that the ground-truth signal
$x^{*}(t)=\sum_{j=1}^{k}v_{j}e^{2\pi\mathbf{i}f_{j}t}$ can be approximated by
$x_{S}(t)=\sum_{j\in S}v_{j}e^{2\pi\mathbf{i}f_{j}t}$, where
$S:=\\{j\in[k]:f_{j}\in\text{some heavy-cluster}~{}C_{i}\\}$ is the set of
frequencies in the heavy-clusters. This step will cause an approximation error
$E_{1}:=\|x^{*}-x_{S}\|_{T}\leq 1.2{\cal N}$888Due to the noisy observations,
not every frequency in the heavy-clusters is recoverable. This gap causes an
extra implicit error term in [CKPS16], which is about $12{\cal N}$. Section
2.4.1 has a more detailed discussion., where ${\cal
N}^{2}:=\|g\|_{T}^{2}+\delta\|x^{*}\|_{T}^{2}$ appears in the approximation
error of the Fourier interpolation problem. However, the frequency in the
heavy-cluster can’t be recovered correctly based on the method they proposed.
Due to the $S$ not having an exact property of one-cluster (See Lemma 12.34),
we will only focus on the frequency with a High signal-noise ratio, which can
correctly satisfy the one-cluster property. Furthermore, we define the
frequency we focus as $x_{S_{f}}(t)=\sum_{j\in
S_{f}}v_{j}e^{2\pi\mathbf{i}f_{j}t}$, where
$S_{f}:=\\{j\in[k]:f_{j}\in\text{some heavy-
cluster}~{}C_{i}~{}\text{and}~{}\text{ bucket}~{}j~{}\text{has high signal-
noise ratio}\\}$

  2. 2.

They solve a Frequency Estimation problem for $x_{S}$ using the filter
techniques and multiple-to-one heavy-cluster reduction, and get a list $L$ of
$O(k)$ candidate frequencies so that for each $j\in S$, $f_{j}$ is close to
some $\widetilde{f}_{p_{j}}\in L$.

  3. 3.

The signal $x_{S}(t)$ can be decomposed into
$\sum_{i=1}^{|L|}e^{2\pi\mathbf{i}\widetilde{f}_{i}t}\cdot x^{*}_{i}(t)$,
where
$x^{*}_{i}(t):=\sum_{j:p_{j}=i}e^{2\pi\mathbf{i}(f_{j}-\widetilde{f}_{i})}$ is
a one-cluster signal with small band-limit. For each $x^{*}_{i}(t)$, they
prove that there exists a low-degree polynomial $P_{i}(t)$ that can
approximate it. Let’s denote the polynomial-approximated signal
$\sum_{i=1}^{|L|}e^{2\pi\widetilde{f}_{i}t}\cdot P_{i}(t)$ by
$x_{S,\mathsf{poly}}(t)$, which has an approximation error
$E_{2}:=\|x_{S}-x_{S,\mathsf{poly}}\|_{T}\leq\delta\|x_{S}\|_{T}$.999We remark
that even if the ground-truth signal $x^{*}(t)$ can be well-approximated by a
mixed Fourier-polynomial signal $\widetilde{x}(t)$, we are unable to recover
_every_ basis of $\widetilde{x}(t)$ due to the limitation of the frequency
estimation procedure. Thus, directly applying linear regression to _partially_
reconstruct $\widetilde{x}(t)$ will not guarantee to be a
$(1+\varepsilon)$-approximation of $x^{*}(t)$.

  4. 4.

It remains to reconstruct $x_{S,\mathsf{poly}}(t)$, which is a variant of
Signal Estimation problem. They use a sampling-and-regression approach to
obtain a $\mathrm{poly}(k)$-Fourier-sparse signal $y(t)$ with an approximation
error $E_{3}:=\|y-x_{S,\mathsf{poly}}\|_{T}\leq 2200{\cal N}$.

By triangle inequality, the total approximation error is $\|y-x^{*}\|_{T}\leq
E_{1}+E_{2}+E_{3}\leq C\cdot{\cal N}$, where $C>1000$ is an absolute constant.

However, the frequency in the heavy-cluster can’t be recovered correctly based
on the method they proposed. Due to the $S$ not having an exact property of a
one-cluster (See Lemma 12.34), we will only focus on the frequency with a high
signal-to-noise ratio that can correctly satisfy the one-cluster property.
Furthermore, we will focus on $x_{S_{f}}(t)=\sum_{j\in
S_{f}}v_{j}e^{2\pi\mathbf{i}f_{j}t},$ where
$S_{f}:=\\{j\in[k]:f_{j}\in\text{some heavy-cluster}~{}C_{i}~{}\text{and the
bucket}~{}j~{}\text{has a high signal-to-noise ratio}\\}$. In the remainder of
this section, we first introduce our techniques to achieve a
$(9+\varepsilon)$-approximation error. Then, we show how to overcome the
barrier and achieve a $(3+\sqrt{2}+\varepsilon)$-approximation error.

####  2.4.1 $(9+\varepsilon)$-approximation algorithm

##### High sensitivity frequency estimation

We first improve $E_{1}$ from $1.2{\cal N}$ to $(1+\varepsilon){\cal N}$ by
proposing a high sensitivity frequency estimation method. More specifically,
to identify the heavy-clusters of the signal $x^{*}$, [CKPS16] designed a
filter function $H$ such that $H\cdot x^{*}$ has high energy in each heavy
cluster $C_{i}$; that is,

$\displaystyle\int_{C_{i}}|\widehat{H\cdot
x^{*}}(f)|^{2}\mathrm{d}f\geq\frac{T}{k}{\cal N}^{2}.$ (2)

Moreover, $H$’s frequencies are contained in a small interval of length
$\Delta$. These two properties imply that for any true frequency $f_{i}\in
C_{i}$, the signal $H\cdot x^{*}$ with frequency domain restricted to
$[f_{i}-\Delta,f_{i}+\Delta]$ is a one heavy-cluster signal with small band-
limit, which allows us to use [PS15]’s approach to estimate $f_{i}$. The
filter function $H$ in [CKPS16] is only $O(1)$-sensitive, which means it can
concentrate a constant fraction of the signal’s energy. And for those less
important frequencies, they cannot be clustered by $H$ and will be lost in the
frequency estimation procedure.

We manage to modify their filter construction and obtain a
$(1-\varepsilon)$-sensitive filter function $H$ such that the signal
$x_{S^{*}}$ consisted of the frequencies in the new heavy-clusters has about
$(1-\varepsilon)$-fraction of energy of $x^{*}$. More specifically, we have
$E^{\mathsf{new}}_{1}=\|x^{*}-x_{S^{*}}\|_{T}\leq(1+\varepsilon){\cal N}$.

To prove that we can actually estimate the frequencies in the new heavy-
clusters, we observe a subtle point: the energy condition of heavy-cluster and
the energy condition of frequency estimation are inconsistent due to the
_noise_ in observations. To be able to estimate the one heavy-cluster signal’s
frequency, it is required that

$\displaystyle\int_{C_{i}}|\widehat{H\cdot
x}(f)|^{2}\mathrm{d}f\geq\frac{T}{k}{\cal N}^{2},$ (3)

which is different from Eq. (2). In other words, not all frequencies in
$S^{*}$ are recoverable, but only most of them. Since [CKPS16] only wants a
constant approximation error, they may simply ignore this difference by losing
a constant factor in accuracy. For us, however, we need to make it precise. We
define $S$ to be a subset of $S^{*}$ containing the frequencies in the heavy-
clusters satisfying Eq. (3). We analyze the effect of $H\cdot g$ and show that
by strengthening the RHS of heavy-cluster’s energy condition (Eq. (2)) to
$\frac{4T}{k}{\cal N}^{2}$, we can bound the unrecoverable part’s energy by

$E_{1.5}:=\|x_{S^{*}}-x_{S}\|_{T}\leq(2+\varepsilon){\cal N}.$

Given the definition of $S_{f}$ as provided earlier, we also need to establish
the upper bound of

$\displaystyle
E_{S_{f}}:=\|x_{S_{f}}-x_{S}\|_{T}\leq(1+\varepsilon)\mathcal{N}.$

For the recoverable part $x_{S_{f}}$, we can just follow [CKPS16]’s approach
to estimating the frequencies in each heavy-cluster.

##### Generalized high-accuracy signal estimation

We apply our four-step Fourier set-query framework to solve the Signal
Estimation problem in the three-step of [CKPS16]’s algorithm and improve
$E_{3}$ from $2200{\cal N}$ to $(4+\varepsilon){\cal N}$. We first define the
problem more formally. By frequency estimation, we obtain a list of candidate
frequencies of $x_{S}$ and in the third step, we know that it can be
approximated by
$x_{S_{f},\mathsf{poly}}(t):=\sum_{i=1}^{|L|}e^{2\pi\mathbf{i}\widetilde{f}_{i}t}\cdot
P_{i}(t)$ with very tiny error $E_{2}$, where $P_{i}(t)$ are some degree-$d$
polynomials. We can rewrite $x_{S,\mathsf{poly}}$ in the _Fourier-monomial
mixed basis_ :

$\displaystyle
x_{S,\mathsf{poly}}(t)=\sum_{i=1}^{|L|}\sum_{j=0}^{d}v_{i,j}\cdot
e^{2\pi\mathbf{i}\widetilde{f}_{i}t}t^{j},$

where $v_{i,j}\in\mathbb{C}$ and $\widetilde{f}_{i}\in L$ are known. That is,
we need to learn $\\{v_{i,j}\\}$ given noisy observations
$x_{S_{f},\mathsf{poly}}(t)+g^{\prime}(t)$, which is a Signal Estimation
problem for the following family of signals:

$\displaystyle{\cal
F}_{\mathsf{mix}}:=\mathrm{span}\left\\{e^{2\pi\mathbf{i}\widetilde{f}_{i}t}\cdot
t^{j}~{}\bigg{|}~{}i\in[|L|],j\in\\{0,\cdots,d\\}\right\\}.$

We apply our three-step framework as follows. In Step 1, we need an energy
bound for ${\cal F}_{\mathsf{mix}}$. [CKPS16] showed that
$R_{\mathsf{mix}}:=\sup_{u(t)\in{\cal
F}_{\mathsf{mix}}}\frac{\sup_{t}|u(t)|^{2}}{\|u\|_{T}^{2}}\leq\widetilde{O}(|L|^{4}d^{4})$.
Then, we get that uniformly sample
$\widetilde{O}(|L|^{4}d^{4}\varepsilon^{-1})$ points in $[0,T]$ gives an
oblivious sketching for $x_{S,\mathsf{poly}}$. Furthermore, we can show that
this sampler is $\varepsilon$-well-balanced. In Step 2, since we aim at
achieving high accuracy, we do not distill the sketch but directly apply the
sharper error analysis to control the energy of the orthogonal part of the
noise, as we did for our high-accuracy signal estimation algorithm (Theorem
1.1). Finally, in Step 3, we solve a weighted linear regression to estimate
the coefficients and obtain a signal $y^{\prime}(t)\in{\cal F}_{\mathsf{mix}}$
such that

$\displaystyle
E^{\mathsf{new}}_{3}:=\|y^{\prime}-x_{S_{f},\mathsf{poly}}\|_{T}\leq(1+\varepsilon)\|x-x_{S_{f},\mathsf{poly}}\|_{T}\leq(4+\varepsilon){\cal
N}.$

Then, we can transform $y$ back to a $\mathrm{poly}(k)$-Fourier-sparse signal
with error $\|y-y^{\prime}\|_{T}\leq E_{2}$.

Combining them together and re-scaling $\varepsilon$ and $\delta$, we get that

$\displaystyle\|y-x^{*}\|_{T}\leq$
$\displaystyle~{}\|y-y^{\prime}\|_{T}+\|y-x_{S_{f},\mathsf{poly}}\|_{T}+\|x_{S_{f},\mathsf{poly}}-x_{S_{f}}\|_{T}+\|x_{S_{f}}-x_{S}\|_{T}+\|x_{S}-x_{S^{*}}\|_{T}+\|x_{S^{*}}-x^{*}\|_{T}$
$\displaystyle\leq$
$\displaystyle~{}E_{2}+E_{3}^{\mathsf{new}}+E_{2}+E_{S_{f}}+E_{1.5}+E_{1}^{\mathsf{new}}$
$\displaystyle\leq$
$\displaystyle~{}(9+\varepsilon)\|g\|_{T}+\delta\|x^{*}\|_{T}.$

Therefore, we obtain a Fourier interpolation algorithm with
$(9+\varepsilon)$-approximation error.

####  2.4.2 $(3+\sqrt{2}+\varepsilon)$-approximation algorithm

How can we further improve this algorithm? We observe that
$E^{\mathsf{new}}_{3}$ can be written more precisely as
$(1+\varepsilon)\|g\|_{T}+(3+\varepsilon){\cal N}$, where the first term
$(1+\varepsilon)\|g\|_{T}$ is identical to the Signal Estimation problem’s
approximation error (Theorem 1.1) and may not be further improved. On the
other hand, $E_{1}^{\mathsf{new}}$ and $E_{1.5}$ only depend on ${\cal N}$. If
we can take a smaller value for ${\cal N}^{2}$, i.e.,
$\varepsilon(\|g\|_{T}^{2}+\delta\|x\|_{T}^{2})$, then we will improve
approximation error. We show that it is possible via an _ultra-high
sensitivity frequency estimation method_.

##### Ultra-high sensitivity frequency estimation

To improve the sensitivity of the frequency estimation method, let ${\cal
N}_{1}^{2}:=\varepsilon{\cal N}^{2}$ and consider the heavy-clusters with
parameter ${\cal N}_{1}$. Let $S^{*}_{1}$ denote the set of frequencies of
$x^{*}(t)$ in the ${\cal N}_{1}$-heavy-clusters. By the same analysis as in
our previous frequency estimation approach, we have
$E_{1}^{\mathsf{new}+}:=\|x^{*}-x_{S_{1}^{*}}\|_{T}\leq(1+\varepsilon){\cal
N}_{1}$.

However, due to the inconsistent energy conditions, only those frequencies in
the heavy-clusters satisfying Eq. (3) are recoverable. Let $S_{1}$ denote the
set of such frequencies, and we need to upper bound
$\|x_{S_{1}^{*}}-x_{S_{1}}\|_{T}$. Previously, we strengthen the heavy-
cluster’s condition (Eq. (2)) and get a $E_{1.5}\leq(2+\varepsilon)\|g\|_{T}$
bound. Here, instead, we relax the RHS of Eq. (3) to
$\varepsilon\cdot\frac{T}{k}{\cal N}_{1}^{2}$. Intuitively, more frequencies
will satisfy the new frequency estimation condition; and if there is a
unrecoverable frequency $f^{*}\in S^{*}_{1}\backslash S_{1}$, it indicates
that its contribution in filtered signal $H\cdot x^{*}$ is cancelled out by
the filtered noise $H\cdot g$. Using this _signal-noise cancellation effect_ ,
we prove that:

$\displaystyle\|H(x_{S_{1}}-x_{S^{*}_{1}})\|_{T}^{2}+\|H(x-x_{S_{1}})\|_{T}^{2}\leq(1+\sqrt{\varepsilon})\|x-x_{S^{*}_{1}}\|_{T}^{2},$
(4)

which saves a factor of 2 from $E_{1.5}$ by introducing an extra term
$\|H(x-x_{S_{1}})\|_{T}$. Recall $\|x-x_{S_{1}}\|_{T}$ is related to
$E_{3}^{\mathsf{new}}$, the error of the signal estimation procedure. We can
decompose it into the “passing energy” $\|H(x-x_{S_{1}})\|_{T}$ and “filtered
energy” $\|(\mathrm{Id}-H)(x-x_{S_{1}})\|_{T}$ and bound them by:

$\displaystyle\|x-x_{S_{1}}\|_{T}\leq\|H(x-x_{S_{1}})\|_{T}+\|g\|_{T}+O(\varepsilon)\|x^{*}-x_{S_{1}}\|_{T}.$

Thus, Eq. (4) can be considered as upper-bounding $E_{1.5}$ and
$E_{3}^{\mathsf{new}}$ simultaneously. Combining them together, we get the
following error guarantee for the frequency recoverable signal $x_{S_{1}}$:

$\displaystyle\|x-x_{S_{1}}\|_{T}+\|x^{*}-x_{S_{1}}\|_{T}\leq(1+\sqrt{2}+O(\sqrt{\varepsilon}))\|g\|_{T}+O(\sqrt{\delta})\|x^{*}\|_{T}.$
(5)

(a) $\|H(t)\cdot(x_{S^{*}_{1}}-x_{S_{1}})(t)\|_{T}$ (b)
$\|H(t)\cdot(x-x_{S_{1}})(t)\|_{T}$

Figure 4:  An illustration of the signal-noise cancellation effect (Eq. (4)).
In (a), the blue region corresponds to the first term of Eq. (4), which
roughly equals the energy of $x-x_{S^{*}_{1}}$. In (b), the red and blue
regions correspond to the second term, where most of the energy is canceled.
Thus, their total energy is very close to $\|x-x_{S^{*}_{1}}\|_{T}$.

Then, by a more careful analysis of the HashToBins approach used by [CKPS16]
for Frequency Estimation, we show that $x_{S_{1}}$’s frequencies can be
efficiently approximated, which gives an ultra-high sensitivity frequency
estimation method.

The remaining part of the algorithm is almost identical to the previous one.
We run the high-accuracy signal estimation algorithm to reconstruct
$x_{S_{1}}$. Let $y(t)$ denote the output Fourier-sparse signal. By Eq. (5)
and re-scaling $\varepsilon$ and $\delta$, we have

$\displaystyle\|y-x^{*}\|_{T}\leq(3+\sqrt{2}+\varepsilon)\|g\|_{T}+\delta\|x^{*}\|_{T}.$

Therefore, we achieve a $(3+\sqrt{2}+\varepsilon)$-approximation error for the
Fourier interpolation.

## Paper Organization

The remainder of the paper is organized as follows. In Sections 3, we formally
define and study the “semi-continuous” Fourier interpolation and set query
problems over lattices. In Section 4, we provide some preliminaries on Fourier
transformation, lattices, etc.

Then, we focus on developing Fourier set query algorithms based on our three-
step framework (Section 2.1). We first build some technical components in
Sections 5 \- 8. Section 5 is for Step 1, where we prove energy bounds and
concentration properties for high-dimensional Fourier-sparse signals and
discrete Fourier-sparse signals. Section 6 is for Step 2, where we describe
fast _oblivious_ sketching methods and we use them as “preconditioners” for
continuous and discrete Fourier-sparse signals in one and higher dimensions.
Sections 7 and 8 are for Step 3. More specifically, in Sections 7, we design a
data structure for improving the time/space complexity of the Randomized BSS
algorithm (Theorem 2.1). Using this data structure, in Section 8 we describe
the _sketch distillation_ technique for different kinds of Fourier-sparse
signals and analyze its robustness to noise. Finally, these components are
wrapped up in Sections 9 \- 11. Section 9 shows sample-optimal and high-
accuracy signal estimation algorithms for one-dimensional semi-continuous
signals (Theorem 1.1). Section 10 generalizes to high-dimensional signals
using lattice theory. Section 11 gives a discrete Fourier set query algorithm
(Theorem 1.2).

In the last part of paper, in Section 12, we give a high-accuracy Fourier
interpolation algorithm for one-dimensional continuous signals, based on our
sharper error control and efficient Fourier set query algorithm.

In the appendix, Section A discusses a special case of the Signal Estimation
problem where the observation has no noise. We present a straightforward
algorithm for one-dimensional Fourier sparse signals. Section B proves that
any signal can be approximated by a semi-continuous signal with the same
sparsity and polynomially-small frequency gap101010We also show an
approximation with $O(1)$-frequency gap but slightly worse sparsity., which
implies a Fourier interpolation algorithm with optimal output-sparsity with a
different error guarantee. Finally, Section C shows that the accuracy of our
Fourier interpolation algorithm can be further improved, if we only care about
the signal in a subset of time duration.

##### Notations.

For any positive integer $n$, we use $[n]$ to denote $\\{1,2,\cdots,n\\}$. We
use $\mathbf{i}$ to denote $\sqrt{-1}$. For a complex number $z\in\mathbb{C}$
where $z=a+\mathbf{i}b$ and $a,b\in\mathbb{R}$. We use $\overline{z}$ to
denote the complex conjugate of $z$, i.e., $\overline{z}=a-\mathbf{i}b$. Then
it is obvious that $|z|^{2}=z\cdot\overline{z}=a^{2}+b^{2}$. We use $f\lesssim
g$ to denote that there exists a constant $C$ such that $f\leq Cg$, and
$f\eqsim g$ to denote $f\lesssim g\lesssim f$. We use $\widetilde{O}(f)$ to
denote $f\log^{O(1)}(f)$. We say $x(t)$ is a $k$-Fourier-sparse when
$x(t)=\sum_{j=1}^{k}v_{j}\exp(2\pi\mathbf{i}f_{j}t)$. We use $\widehat{x}(f)$
to denote the Fourier transform of $x(t)$. More specifically,
$\widehat{x}(f)=\int_{-\infty}^{\infty}x(t)\exp(-2\pi\mathbf{i}ft)\mathrm{d}t$.
We define our discrete norm as $\|g(t)\|^{2}_{W}=\frac{1}{|W|}\sum_{t\in
W}|g(t)|^{2}$ for function $g$. We define our weighted discrete norm as
$\|g(t)\|^{2}_{S,w}=\sum_{t\in S}w_{t}|g(t)|^{2}$ for function $g$. We define
the continuous $T$-norm as
$\|g(t)\|^{2}_{T}=\frac{1}{T}\int_{0}^{T}|g(t)|^{2}\mathrm{d}t$ for function
$g$.

In general, we assume $x^{*}(t)$ is our ground truth and is a $k$-Fourier-
sparse signal. We can observe function $x(t)=x^{*}(t)+g(t)$ for $g(t)$ being a
noise function. We can observe $x(t)$ in duration $[0,T]$. The ground truth
$x^{*}(t)$ has frequencies in $[-F,F]$.

##  3 Definitions of Semi-Continuous Fourier Set Query and Interpolation

In this section, we give the formal definitions of the problems studied in
this paper. In Section 3.1, we define the Fourier set query for discrete and
continuous signals. In Section 3.2, we define the Fourier interpolation
problem and its two sub-problems: frequency estimation and signal estimation.

###  3.1 Formal definitions of Fourier set query

The discrete Fourier set query problem is defined as follows:

######  Definition 3.1 (Discrete Fourier set query problem).

Let $x\in\mathbb{C}^{n}$ and $\widehat{x}$ be its discrete Fourier
transformation. Let $\varepsilon>0$. Given a set $S\subseteq[n]$ and query
access to $x$, the goal is to use a few queries to compute a vector
$x^{\prime}$ with support $\mathrm{supp}(x^{\prime})\subseteq S$ such that

$\displaystyle\|(x^{\prime}-\widehat{x})_{S}\|_{2}^{2}\leq\varepsilon\cdot\|\widehat{x}_{[n]\backslash
S}\|_{2}^{2}.$

We also define the continuous Fourier set query problem as follows:

######  Definition 3.2 (Continuous Fourier set query problem).

For $d\geq 1$, let $x^{*}(t)$ be a signal in time duration $[0,T]^{d}$. Let
$\widehat{x^{*}}(f)$ denote the continuous Fourier transformation of
$x^{*}(t)$. Let $\varepsilon>0$. Given a set $S\subseteq\mathbb{R}^{d}$ of
frequencies such that $\mathrm{supp}(\widehat{x^{*}})\subseteq S$, and
observations of the form $x(t)=x^{*}(t)+g(t)$, where $g(t)$ denotes the noise.
The goal is to output a Fourier-sparse signal $x^{\prime}(t)$ with support
$\mathrm{supp}(\widehat{x^{\prime}})\subseteq S$ such that

$\displaystyle\|x^{\prime}-x^{*}\|_{T}^{2}\leq(1+\varepsilon)\cdot\|g\|_{T}^{2}.$

###  3.2 Formal definitions of semi-continuous Fourier interpolation

In this section, we provide the following formal definition of the semi-
continuous Fourier interpolation problem, where we assume that the frequencies
of the signal are contained in a lattice.

######  Problem 3.3 (Semi-continuous Fourier interpolation problem).

Given a basis $\mathcal{B}$ of $m$ known vectors $b_{1},b_{2},\cdots
b_{m}\in\mathbb{R}^{d}$, let $\Lambda(\mathcal{B})\subset\mathbb{R}^{d}$
denote the lattice

$\displaystyle\Lambda(\mathcal{B})=\Big{\\{}z\in\mathbb{R}^{d}:z=\sum_{i=1}^{m}c_{i}b_{i},c_{i}\in\mathbb{Z},\forall
i\in[m]\Big{\\}}$

Suppose that $f_{1},f_{2},\cdots,f_{k}\in\Lambda(\mathcal{B})$, $\forall
i\in[k],|f_{i}|\leq F$. Let
$x^{*}(t)=\sum_{j=1}^{k}v_{j}e^{2\pi\mathbf{i}\langle f_{j},t\rangle}$, and
let $g(t)$ denote the noise. Given observations of the form
$x(t)=x^{*}(t)+g(t)$, $t\in[0,T]^{d}$. Let $\eta=\min_{i\neq
j}\|f_{j}-f_{i}\|_{\infty}$. There are three goals:

  1. 1.

The first goal is to design an algorithm that output
$f_{1},f_{2},\cdots,f_{k}$ _exactly_ given query access to the signal $x(t)$
for $t\in[0,T]^{d}$.

  2. 2.

The second goal is to design an algorithm that output a set $L$ of frequencies
such that, for each $f_{i}$, there is $f^{\prime}_{i}\in L$,
$\|f_{i}-f^{\prime}_{i}\|_{2}\leq D/T$.

  3. 3.

The third goal is to design an algorithm that output
$y(t)=\sum_{j=1}^{\widetilde{k}}v_{j}^{\prime}\cdot
e^{2\pi\mathbf{i}f_{j}^{\prime}t}$ such that
$\int_{[0,T]^{d}}|y(t)-x(t)|^{2}\mathrm{d}t\lesssim\int_{[0,T]^{d}}|g(t)|^{2}\mathrm{d}t$.

Then, we extract two sub-problems from Problem 3.3: Frequency Estimation and
Signal Estimation. We give their definitions below.

We first define the $d$-dimensional frequency estimation under the semi-
continuous as follows. In this problem, we want to recover each frequencies in
a small range.

######  Problem 3.4 (Frequency estimation).

Given a basis $\mathcal{B}$ of $m$ known vectors $b_{1},b_{2},\cdots
b_{m}\in\mathbb{R}^{d}$, let $\Lambda(\mathcal{B})\subset\mathbb{R}^{d}$
denote the lattice

$\displaystyle\Lambda(\mathcal{B})=\Big{\\{}z\in\mathbb{R}^{d}:z=\sum_{i=1}^{m}c_{i}b_{i},c_{i}\in\mathbb{Z},\forall
i\in[m]\Big{\\}}$

Suppose that $f_{1},f_{2},\cdots,f_{k}\in\Lambda(\mathcal{B})$. Let
$x^{*}(t)=\sum_{j=1}^{k}v_{j}e^{2\pi\mathbf{i}\langle f_{j},t\rangle}$, and
let $g(t)$ denote the noise. Given observations of the form
$x(t)=x^{*}(t)+g(t)$, $t\in[0,T]^{d}$. Let $\eta=\min_{i\neq
j}\|f_{j}-f_{i}\|_{\infty}$.

The goal is to design an algorithm that output a set $L$ of frequencies such
that, for each $f_{i}$, there is $f^{\prime}_{i}\in L$,
$\|f_{i}-f^{\prime}_{i}\|_{2}\leq D/T$.

We remark that the recovered frequencies in $L$ are not necessary to be in
$\Lambda(\mathcal{B})$, and $D$ is a parameter that can depend on $k$.

Next, we define the $d$-dimensional Signal Estimation under the semi-
continuous setting as follows. In this problem, we want to recover a signal
that can approximate the ground-truth signal in the time domain.

######  Problem 3.5 (Signal Estimation problem).

Given a basis $\mathcal{B}$ of $m$ known vectors $b_{1},b_{2},\cdots
b_{m}\in\mathbb{R}^{d}$, let $\Lambda(\mathcal{B})\subset\mathbb{R}^{d}$
denote the lattice

$\displaystyle\Lambda(\mathcal{B})=\Big{\\{}z\in\mathbb{R}^{d}:z=\sum_{i=1}^{m}c_{i}b_{i},c_{i}\in\mathbb{Z},\forall
i\in[m]\Big{\\}}$

Suppose that $f_{1},f_{2},\cdots,f_{k}\in\Lambda(\mathcal{B})$. Let
$x^{*}(t):=\sum_{j=1}^{k}v_{j}e^{2\pi\mathbf{i}\langle f_{j},t\rangle}$, and
let $g(t)$ denote the noise. Given observations of the form
$x(t)=x^{*}(t)+g(t)$, $t\in[0,T]^{d}$. Let $\eta=\min_{i\neq
j}\|f_{j}-f_{i}\|_{\infty}$.

The goal is to design an algorithm that outputs
$y(t)=\sum_{j=1}^{\widetilde{k}}v_{j}^{\prime}\cdot
e^{2\pi\mathbf{i}f_{j}^{\prime}t}$ such that

$\int_{[0,T]^{d}}|y(t)-x(t)|^{2}\mathrm{d}t\lesssim\int_{[0,T]^{d}}|g(t)|^{2}\mathrm{d}t.$

Note that outputting $y(t)=\sum_{j=1}^{\widetilde{k}}v_{j}^{\prime}\cdot
e^{2\pi\mathbf{i}f_{j}^{\prime}t}$ means outputting
$\\{v_{j}^{\prime},f_{j}^{\prime}\\}_{j\in[\widetilde{k}]}$.

######  Remark 3.6.

We note that given the solution of Frequency Estimation (Problem 3.4), Signal
Estimation (Problem 3.5) can be formulated as a Fourier set query problem
(Problem 3.2). More specifically, by Frequency Estimation, we will find a set
that contains all frequencies of the ground truth signal $x^{*}(t)$. Then, we
only need to recover the coefficients with frequencies in this set, which is
equivalent to a set query problem.

##  4 Preliminaries

This section is organized as follows. In Section 4.1, we provide some
technical tools in probability theory and linear algebra. In Section 4.2, we
review the Fourier transformation for different types of signals. In Section
4.3, we show some facts about Lattices. And in Section 4.4, we discuss the
importance sampling method.

###  4.1 Tools and inequalities

######  Lemma 4.1 (Chernoff Bound [Che52]).

Let $X_{1},X_{2},\cdots,X_{n}$ be independent random variables. Assume that
$0\leq X_{i}\leq 1$ always, for each $i\in[n]$. Let
$X=X_{1}+X_{2}+\cdots+X_{n}$ and
$\mu=\mathbb{E}[X]=\overset{n}{\underset{i=1}{\sum}}\mathbb{E}[X_{i}]$. Then
for any $\varepsilon>0$,

$\mathsf{Pr}[X\geq(1+\varepsilon)\mu]\leq\exp(-\frac{\varepsilon^{2}}{2+\varepsilon}\mu)\textit{
and
}\mathsf{Pr}[X\leq(1-\varepsilon)\mu]\leq\exp(-\frac{\varepsilon^{2}}{2}\mu).$

######  Definition 4.2 ($\varepsilon$-net).

Let $T$ be a metric space with distance measure $d$. Consider a subset
$K\subset T$ and let $\varepsilon>0$. A subset ${\cal N}\subseteq K$ is called
an $\varepsilon$-net of $K$ if every point in $K$ is within distance
$\varepsilon$ of some point of ${\cal N}$, i.e.

$\displaystyle\forall x\in K,\exists y\in{\cal
N}~{}\text{s.t.}~{}d(x,y)\leq\varepsilon.$

######  Fact 4.3 (Fast matrix multiplication).

We use ${\cal T}_{\mathrm{mat}}(a,b,c)$ to denote the time of multiplying an
$a\times b$ matrix with another $b\times c$ matrix.

We use $\omega$ to denote the exponent of matrix multiplication, i.e., ${\cal
T}_{\mathrm{mat}}(n,n,n)=n^{\omega}$. Currently $\omega\approx 2.373$ [Wil12,
LG14, AW21].

######  Fact 4.4 (Weighted linear regression).

Given a matrix $A\in\mathbb{C}^{n\times d}$, a vector $b\in\mathbb{C}^{n}$ and
a weight vector $w\in\mathbb{R}_{>0}^{n}$, it takes $O(nd^{\omega-1})$ time to
output an $x^{\prime}$ such that

$x^{\prime}=\underset{x}{\arg\min}\|\sqrt{W}(Ax-b)\|_{2}=(A^{*}WA)^{-1}A^{*}Wb.$

where
$\sqrt{W}:=\mathrm{diag}(\sqrt{w_{1}},\dots,\sqrt{w_{n}})\in\mathbb{R}^{n\times
n}$, and $\omega\approx 2.373$ is the exponent of matrix multiplication
[Wil12, LG14, AW21].

######  Fact 4.5.

For any $x\in(0,1)$, we have $\cos(x)\leq\exp(-x^{2}/2)$.

###  4.2 Basics of Fourier transformation

The definition of high dimensional Fourier transform is as follows:

$\displaystyle\widehat{x}(f)=\int_{(-\infty,\infty)^{d}}{x}(t)\exp(-2\pi\mathbf{i}\langle
f,t\rangle)\mathrm{d}t,\text{~{}where~{}}f\in\mathbb{R}^{d},$

and the definition of high dimensional inverse Fourier transform is as
follows:

$\displaystyle
x(t)=\int_{(-\infty,\infty)^{d}}\widehat{x}(f)\exp(2\pi\mathbf{i}\langle
f,t\rangle)\mathrm{d}f,\text{~{}where~{}}t\in\mathbb{R}^{d}.$

Note that when we replace $d=1$ in the definition of high dimensional Fourier
transform and inverse Fourier transform above, we get the definition of one-
dimensional Fourier transform and inverse Fourier transform.

The definition of discrete Fourier transform is as follows:

$\displaystyle\widehat{x}_{f}=\sum_{t=1}^{n}{x}_{t}\exp(-2\pi\mathbf{i}ft/n),\text{~{}where~{}}f\in[n],$

and the definition of discrete inverse Fourier transform is as follows:

$\displaystyle
x_{t}=\frac{1}{n}\sum_{f=1}^{n}\widehat{x}_{f}\exp(2\pi\mathbf{i}ft/n),\text{~{}where~{}}t\in[n].$

A continuous $k$-Fourier sparse signal
$x(t):\mathbb{R}^{d}\rightarrow\mathbb{C}$ can be represented as follows:

$\displaystyle x(t)=\sum_{j=1}^{k}v_{j}\exp(2\pi\mathbf{i}\langle
f_{j},t\rangle),~{}v_{j}\in\mathbb{C},f_{j}\in\mathbb{R}^{d},~{}\forall
j\in[k].$

Thus, $\widehat{x}(f)$ is:

$\displaystyle\widehat{x}(f)=\sum_{j=1}^{k}v_{j}\delta(t-f_{j}).$

A discrete $k$-Fourier sparse signal $x\in\mathbb{C}^{n}$ can be represented
as follows:

$\displaystyle x_{t}=\sum_{j\in
S}v_{j}\exp(2\pi\mathbf{i}jt/n),~{}S\subseteq[n],|S|=k,~{}v_{j}\in\mathbb{C},\forall
j\in S.$

So, $\widehat{x}_{f}$ is:

$\displaystyle\widehat{x}_{f}=\left\\{\begin{array}[]{ll}v_{j}&,j\in S\\\
0&,\text{o.w.}\end{array}\right.$

###  4.3 Facts about Lattices

######  Definition 4.6 (Lattice).

A lattice ${\cal L}$ in $\mathbb{R}^{d}$ is defined as follows:

$\displaystyle{\cal
L}:=\Big{\\{}\sum_{i=1}^{k}\lambda_{i}b_{i}:\lambda_{1},\dots,\lambda_{k}\in\mathbb{Z}\Big{\\}},$

where $b_{1},\dots,b_{k}\in\mathbb{R}^{d}$ are linearly independent vectors.
And we denote the matrix
$B:=\begin{bmatrix}b_{1}&\cdots&b_{k}\end{bmatrix}\in\mathbb{R}^{n\times k}$
as the basis of the lattice ${\cal L}$.

######  Definition 4.7 (Fundamental parallelepiped).

For a lattice ${\cal L}$ with basis $B$, its fundamental parallelepiped is
defined to be:

$\displaystyle{\cal P}(B):=\\{Bx~{}|~{}x\in[0,1)^{d}\\}.$

######  Fact 4.8.

For any lattice with basis $B$, we have

$\displaystyle\mathrm{vol}({\cal P}(B))=\sqrt{\det(B^{\top}B)}.$

In particular, if $B$ is full-rank, $\mathrm{vol}({\cal P}(B))=|\det(B)|$.

######  Lemma 4.9 (The number of lattice points within a ball 111111We thank
Thomas Rothvoss for providing the proof of this bound. ).

Let ${\cal L}$ be any lattice with basis $B$ such that the spectral norm
$\|B\|\leq\ell$. Then, the number of lattice points inside a ball centered at
0 with radius $R$ is upper bounded by:

$\displaystyle|{\cal L}\cap
B_{d}(0,R)|\leq(1+\frac{\sqrt{k}\ell}{R})^{k}\cdot\frac{\mathrm{vol}(B_{k}(0,R))}{\mathrm{vol}({\cal
P}(B))}.$

###### Proof.

We first show that for two different lattice points $x,y\in{\cal L}\cap
B_{d}(0,R)$, the translations of ${\cal P}$ at $x$ and at $y$ are disjoint,
i.e., $(x+{\cal P})\cap(y+{\cal P})=\emptyset$.

Suppose $(x+{\cal P})\cap(y+{\cal P})\neq\emptyset$ for some $x,y\in{\cal
L}\cap B_{d}(0,R)$. Then, we have

$\displaystyle
x+\sum_{i=1}^{k}\lambda_{i}b_{i}=y+\sum_{i=1}^{k}\tau_{i}b_{i},$

where $\lambda_{i},\tau_{i}\in[0,1)$. It gives

$\displaystyle x-y=\sum_{i=1}^{k}(\tau_{i}-\lambda_{i})b_{i}.$

Note that $x-y\in{\cal L}$, which means

$\displaystyle\sum_{i=1}^{k}(\tau_{i}-\lambda_{i})b_{i}=\sum_{i=1}^{k}c_{i}b_{i},$

where $c_{1},\dots,c_{k}\in\mathbb{Z}$. Since $\tau_{i}-\lambda_{i}\in(-1,1)$,
we get that $c_{i}=0$ for all $i\in[k]$. Thus, $x=y$.

Then, for any point $y\in x+{\cal P}$, where $x\in{\cal L}\cap B_{d}(0,R)$, we
have

$\displaystyle\|y\|_{2}=\|x+z\|_{2}\leq\|x\|_{2}+\|z\|_{2}\leq
R+\sqrt{k}\ell,$ (6)

where $z\in{\cal P}$ and the last step follows from $x\in B_{d}(0,R)$ and
$\|z\|_{2}=\|B\lambda\|_{2}\leq\|B\|\|\lambda\|_{2}\leq\sqrt{k}\ell$, for some
$\lambda\in[0,1)^{k}$.

Then, we have

$\displaystyle|{\cal L}\cap B_{d}(0,R)|\leq$ $\displaystyle~{}\sum_{x\in{\cal
L}\cap B(0,R)}\frac{\mathrm{vol}(x+{\cal P})}{\mathrm{vol}({\cal P}(B))}$
$\displaystyle\leq$
$\displaystyle~{}\frac{\mathrm{vol}(B_{k}(0,R+\sqrt{k}\ell))}{\mathrm{vol}({\cal
P}(B))}$ $\displaystyle\leq$
$\displaystyle~{}(1+\frac{\sqrt{k}\ell}{R})^{k}\cdot\frac{\mathrm{vol}(B_{k}(0,R))}{\mathrm{vol}({\cal
P}(B))},$

where the second step follows from the disjointness of translations and the
bound on the total width (Eq. (6)).

The lemma is then proved. ∎

We define the shortest vector problem (SVP) as follows:

######  Definition 4.10 (SVP).

Let ${\cal L}$ denote a Lattice. We define $\mathrm{SVP}({\cal L})$,

$\displaystyle\mathrm{SVP}({\cal L}):=\min\\{\|x\|_{2}~{}|~{}x\in{\cal
L}\backslash\\{{\bf 0}\\}\\}.$

Given a basis of ${\cal L}$, the goal is to compute $\mathrm{SVP}({\cal L})$.

In fact, compute SVP (or even approximations of SVP) is an NP-hard problem.
The following theorem shows a well-known lower bound for the shortest vector
length.

######  Theorem 4.11 (Theorem 1.10 in [Rot16], a lower bound on shortest
vector).

Let ${\cal L}$ denote a lattice with basis ${\cal B}$. Let
$(b_{1}^{*},\cdots,b_{n}^{*})$ be its Gram-Schmit orthogonalization. Then

$\displaystyle\mathrm{SVP}({\cal L})\geq\min_{i\in[n]}\|b_{i}^{*}\|_{2}$

######  Fact 4.12.

The Gram–Schmidt process takes a finite, linearly independent set of vectors
$S=\\{v_{1},\cdots,v_{k}\\}$ for $k\leq n$, runs $O(nk^{2})$ time, and
generates an orthogonal set $S^{\prime}=\\{u_{1},\cdots,u_{k}\\}$ that spans
the same $k$-dimensional subspace of $\mathbb{R}^{n}$ as $S$.

###  4.4 Facts about importance sampling

Important sampling try to estimate a statistic value in one distribution by
taking samples in another distribution. In particular, [CP19a] considered the
importance sampling for estimating the norm of functions in a linear family
${\cal F}$.

In this followings, we first provide some basic definitions about linear
function family.

######  Definition 4.13 (Condition number of sampling distribution).

Let $G$ be any domain and ${\cal F}$ is a linear function family from $G$ to
$\mathbb{C}$. Let $D$ be an arbitrary distribution over $G$. Then the
condition number of $D$ with respect to ${\cal F}$ is defined as follows:

$K_{D}:=\sup_{t\in
G}\sup_{f\in{\mathcal{F}}}\frac{|f(t)|^{2}}{\|f\|_{D}^{2}},$

where

$\displaystyle\|f\|_{D}^{2}:=\int_{G}D(t)\cdot|f(t)|^{2}\mathrm{d}t.$

######  Definition 4.14 (Orthonormal basis for linear function family).

Let $G$ be any domain. Given a linear function family $\mathcal{F}$ from $G$
to $C$, and a probability distribution $D$ over $G$. We say
$\\{v_{1},\ldots,v_{d}\\}$ form an orthonormal basis of $\mathcal{F}$ with
respect to $D$, if they satisfy the following properties:

  * •

for any $i,j\in[d]$, $\int_{G}D(t)v_{i}(t)\overline{v_{j}(t)}\mathrm{d}t={\bf
1}_{i=j}$, and

  * •

for any $f\in{\cal F}$, $f\in\mathrm{span}\\{v_{1},\dots,v_{d}\\}$.

######  Fact 4.15.

Let $\\{v_{1},\dots,v_{k}\\}$ be an orthonormal basis of ${\cal F}$ with
respect to $D$. For any function $f\in\mathcal{F}$, let $\alpha(f)$ denote the
coefficients under the basis $\\{v_{1},\ldots,v_{d}\\}$, i.e.,
$h=\sum_{i=1}^{d}\alpha(h)_{i}\cdot v_{i}$. Then,

$\|\alpha(h)\|_{2}=\|h\|_{D}.$

For an unknown function $f\in{\cal F}$, the goal of importance sampling is to
estimate $\|f\|_{D}$, given samples from another distribution $D^{\prime}$.
The following definition introduces the importance sampling procedure and
condition number of the importance sampling distribution.

######  Definition 4.16 (Definition 3.1 of [CP19a]).

For any unknown distribution $D^{\prime}$ over the domain $G$ and any function
$f\in{\cal F}$, let
$f^{(D^{\prime})}(t):=\sqrt{\frac{D(t)}{D^{\prime}(t)}}\cdot f(t)$ be the
importance sampling function for some known distribution $D$ such that

$\underset{t\sim
D^{\prime}}{\operatorname*{{\mathbb{E}}}}\left[|f^{(D^{\prime})}(t)|^{2}\right]=\underset{t\sim
D^{\prime}}{\operatorname*{{\mathbb{E}}}}\left[\frac{D(t)}{D^{\prime}(t)}|f(t)|^{2}\right]=\underset{t\sim
D}{\operatorname*{{\mathbb{E}}}}\left[|f(t)|^{2}\right].$

Then, we can use samples from $D^{\prime}$ to estimate
$\|f^{(D^{\prime})}\|_{D^{\prime}}$, which gives an estimate of $\|f\|_{D}$.

When the family $\mathcal{F}$ and $D$ is clear, we use
$K_{\mathsf{IS},D^{\prime}}$ to denote the condition number of importance
sampling from $D^{\prime}$:

$\displaystyle
K_{\mathsf{IS},D^{\prime}}=\underset{t}{\sup}\left\\{\underset{f\in\mathcal{F}}{\sup}\left\\{\frac{|f^{(D^{\prime})}(t)|^{2}}{\|f^{D^{\prime}}\|_{D^{\prime}}^{2}}\right\\}\right\\}=\underset{t}{\sup}\bigg{\\{}\frac{D(t)}{D^{\prime}(t)}\cdot\underset{f\in\mathcal{F}}{\sup}\big{\\{}\frac{|f(t)|^{2}}{\|f\|_{D}^{2}}\big{\\}}\bigg{\\}}.$
(7)

From Definition 4.16, we know that the efficiency of importance sampling
depends on how many samples we need to estimate
$\|f^{D^{\prime}}\|_{D^{\prime}}$. The following lemma provide a criteria for
judging whether a set of samples gives a good estimation for the norm of
function.

######  Lemma 4.17 (Lemma 4.2 in [CP19a]).

For any $\varepsilon\in(0,1)$, let $S=\\{t_{1},\dotsc,t_{s}\\}$ and the weight
vector $w\in\mathbb{R}_{>0}^{s}$. Define a matrix $A\in\mathbb{R}^{s\times d}$
be the $s\times d$ matrix defined as $A_{i,j}=\sqrt{w_{i}}\cdot v_{j}(t_{i})$,
where $\\{v_{1},\dots,v_{d}\\}$ is an orthonormal basis for ${\cal F}$. Then

$\|h\|^{2}_{S,w}:=\sum_{j=1}^{s}w_{j}\cdot|h(x_{j})|^{2}\in[1\pm\varepsilon]\cdot\|h\|_{D}^{2}\text{
\quad for every }h\in{\mathcal{F}}$

if and only if the eigenvalues of $A^{*}A$ are in
$[1-\varepsilon,1+\varepsilon]$.

The following lemma shows that the sample complexity depends on the condition
number $K_{\mathsf{IS},D^{\prime}}$:

######  Lemma 4.18 (Lemma 6.6 in [CP19a]).

Let $D^{\prime}$ be an arbitrary distribution over $G$ and let
$K_{\mathsf{IS},D^{\prime}}$ be the condition number of importance sampling
from $D^{\prime}$ (defined by Eq. (7)). There exists an absolute constant $C$
such that for any $\varepsilon\in(0,1)$ and $\delta\in(0,1)$, let
$S=\\{t_{1},\dots,t_{s}\\}$ be a set of i.i.d. samples from the distribution
$D^{\prime}$ and let $w$ be the weight vector defined by
$w_{j}=\frac{D(t_{j})}{s\cdot D^{\prime}(t_{j})}$ for each $j\in[s]$. Then, as
long as

$\displaystyle s\geq\frac{C}{\varepsilon^{2}}\cdot
K_{\mathsf{IS},D^{\prime}}\log\frac{d}{\delta},$

the $s\times d$ matrix $A_{i,j}=\sqrt{w_{i}}\cdot v_{j}(t_{i})$ satisfies

$\|A^{*}A-I\|_{2}\leq\varepsilon\text{ with probability at least }1-\delta.$

##  5 Energy Bounds for Fourier Signals

The energy bound shows that the maximum value of a Fourier sparse signal in a
certain interval can be bounded by its energy on the interval. One interesting
fact is that the approximation ratio in the energy bound is only relate to the
sparsity $k$, and have no relationship with time duration $T$ and band-limit
$F$. An application of energy bound is preserving the norm, that is what is
the least size of set $S$, such that $\|f\|_{S}=\|f\|_{T}$, for any function
$f$ in a certain function family. The relationship between energy bound and
norm preserving can be build by Chernoff bound.

[BE06, Kós08, CKPS16, CP19b] proved energy bounds for sparse Fourier signal
under one-dimensional continuous Fourier transform. We further generalize
these results to discrete Fourier sparse signal under discrete Fourier
transform and high-dimensional Fourier sparse signal under continuous Fourier
transform.

This section is organized as follows:

  * •

Section 5.1 reviews previous results for one-dimensional continuous Fourier-
sparse signals.

  * •

Section 5.2 proves a new energy bound for high dimensional Fourier-sparse
signals, and also gives a nearly matching lower bound.

  * •

Section 5.3 proves energy bound for discrete Fourier-sparse signals.

  * •

Section 5.4 builds the connection between energy bound and the concentration
property.

###  5.1 Energy bound for one-dimensional signals

In this section, we review the energy bound proved in prior work [BE06, Kós08,
CKPS16, CP19b].

[Kós08] proved the following energy bound:

######  Theorem 5.1 ([Kós08, CKPS16]).

Define a family of $F$-band-limit, $k$-sparse Fourier signals:

$\displaystyle{\cal F}:=\Big{\\{}x(t)=\sum_{j=1}^{k}v_{j}\cdot
e^{2\pi\mathbf{i}f_{j}t}~{}\Big{|}~{}f_{j}\in\mathbb{R}\cap[-F,F]\Big{\\}}$

Then, for any $t\in(-1,1)$,

$\displaystyle\underset{x\in{\cal
F}}{\sup}\frac{|x(t)|^{2}}{\|x\|_{D}^{2}}\lesssim k^{2}.$

[BE06] also proved a time-dependent energy bound for one-dimensional signal:

######  Theorem 5.2 ([BE06, CP19a]).

Define a family of $F$-band-limit, $k$-sparse Fourier signals:

$\displaystyle{\cal F}:=\Big{\\{}x(t)=\sum_{j=1}^{k}v_{j}\cdot
e^{2\pi\mathbf{i}f_{j}t}~{}\Big{|}~{}f_{j}\in\mathbb{R}\cap[-F,F]\Big{\\}}$

Then, for any $t\in(-1,1)$,

$\displaystyle\underset{x\in{\cal
F}}{\sup}\frac{|x(t)|^{2}}{\|x\|_{D}^{2}}\lesssim\frac{k}{1-|t|}.$

###  5.2 Energy bound for high-dimensional signals

The goal of this section is to prove Theorem 5.3, which gives an energy bound
for $d$-dimensional Fourier signal. It can be viewed as a $d$-dimensional
version of [CKPS16, Lemma 5.1]. We also prove a lower bound in Lemma 5.5.

######  Theorem 5.3 (Energy bound in $d$-dimensional).

For any $d$-dimensional $k$-Fourier-sparse signal
$x(t):\mathbb{R}^{d}\rightarrow\mathbb{C}$ and any duration $T$, we have

$\displaystyle\max_{t\in[0,T]^{d}}|x(t)|^{2}\leq k^{O(d)}\|x\|_{T}^{2},$

where $\|x\|_{T}^{2}=\frac{1}{T^{d}}\int_{[0,T]^{d}}|x(t)|^{2}\mathrm{d}t$.

###### Proof.

Without loss of generality, we fix $T=1$. Then
$\|x\|_{T}^{2}=\int_{[0,1]^{d}}|x(t)|\mathrm{d}t$. Because $\|x\|_{T}^{2}$ is
the average over the interval $[0,T]^{d}$, if the maximizer
$t^{*}=\arg\max_{t\in[0,T]^{d}}|x(t)|^{2}$ is not $0^{d}$ or $T=1$, we can
scale the two intervals $[0^{d},t^{*}]$ and $[t^{*},T^{d}]$ to $[0,1]$ and
prove the desired property separately. Hence we assume that
$|x(0)|^{2}=\max_{t\in[0,T]}|x(t)|^{2}$ in the proof.

In the next a few paragraphs, we show how to use Lemma 5.4 to prove Theorem
5.3.

We use $0^{d}$ to denote a length-$d$ vector with $0$ everywhere. Due to Lemma
5.4, we can choose $t_{0}=0^{d}$ such that $\forall\tau\in\mathbb{R}_{>0}^{d}$
there exist $C_{1},\cdots,C_{m}\in\mathbb{C}$, and

$\displaystyle x(0^{d})=\sum_{j\in[m]}C_{j}\cdot x(j\cdot\tau).$

By the Cauchy-Schwarz inequality, it implies that for any $\tau$,

$\displaystyle|x(0^{d})|^{2}\leq
m\sum_{j\in[m]}|C_{j}|^{2}|x(j\cdot\tau)|^{2}$

Then, we obtain

$\displaystyle|x(0^{d})|^{2}=$
$\displaystyle~{}m^{d}\int_{[0,1/m]^{d}}|x(0^{d})|^{2}\mathrm{d}\tau$
$\displaystyle\lesssim$
$\displaystyle~{}m^{d}\cdot\int_{[0,1/m]^{d}}\left(m\sum_{j=1}^{m}|x(j\cdot\tau)|^{2}\right)\mathrm{d}\tau$
$\displaystyle=$
$\displaystyle~{}m^{d+1}\cdot\sum_{j=1}^{m}\int_{[0,1/m]^{d}}|x(j\cdot\tau)|^{2}\mathrm{d}\tau$
$\displaystyle=$
$\displaystyle~{}m^{d+1}\cdot\sum_{j=1}^{m}\frac{1}{j^{d}}\int_{[0,j/m]^{d}}|x(\tau)|^{2}\mathrm{d}\tau$
$\displaystyle\leq$
$\displaystyle~{}m^{d+1}\cdot\sum_{j=1}^{m}\frac{1}{j^{d}}\int_{[0,1]^{d}}|x(\tau)|^{2}\mathrm{d}\tau$
$\displaystyle\lesssim$ $\displaystyle~{}m^{d+1}\log m\cdot\|x\|_{T}^{2}$
$\displaystyle\leq$ $\displaystyle~{}k^{O(d)}\|x\|_{T}^{2},$ (8)

where the third step follows by moving $m$ outside of the integral and
swapping the integration and the summation, the fourth step follows by
replacing $j\tau$ by $\tau$, the fifth step follows by $j/m\leq 1$, the sixth
step follows by $\sum_{j=1}^{m}1/j^{d}\leq\sum_{j=1}^{m}1/j=O(\log m)$ and the
definition of $\|x\|_{T}^{2}$, and the last step follows from Lemma 5.4 that
$m=\mathrm{poly}(k)$.

Thus, we have the desired bound. ∎

The following lemma shows that each point of the signal can be expressed as a
linear combination of about $k^{2}$ equally spaced signal points.

######  Lemma 5.4 ($d$-dimensional signal interpolation).

For any $k$ and $d$, there exists $m=O(k^{2}\log k)$ such that for any
$d$-dimensional $k$-Fourier-sparse signal $x(t)$, any
$t_{0}\in\mathbb{R}^{d}_{\geq 0}$ and $\tau\in\mathbb{R}_{>0}^{d}$, there
always exist $C_{1},C_{2},\cdots,C_{m}\in\mathbb{C}$ such that the following
properties hold,

$\displaystyle\mathrm{Property~{}\@slowromancap i@}$
$\displaystyle\quad~{}|C_{j}|\leq 11~{}\text{~{}for~{}all~{}}j\in[m],$
$\displaystyle\mathrm{Property~{}\@slowromancap ii@}$
$\displaystyle\quad~{}x(t_{0})=\sum_{j\in[m]}C_{j}\cdot x(t_{0}+j\cdot\tau).$

###### Proof.

Consider a specific signal
$x(t):=\sum_{i=1}^{k}v_{i}e^{2\pi\mathbf{i}f_{i}^{\top}t}$ for
$t\in\mathbb{R}^{d}$, where $f_{i}\in\mathbb{R}^{d}$ are given. We fix
$t_{0}\in\mathbb{R}^{d}$ and $\tau\in\mathbb{R}^{d}$, and then rewrite
$x(t_{0}+j\cdot\tau)$ as a polynomial of $b_{i}:=v_{i}\cdot
e^{2\pi\mathbf{i}f_{i}^{\top}t_{0}}$ and
$z_{i}:=e^{2\pi\mathbf{i}f_{i}^{\top}\tau}$ for each $i\in[k]$.

$\displaystyle x(t_{0}+j\cdot\tau)=$
$\displaystyle~{}\sum_{i=1}^{k}v_{i}e^{2\pi\mathbf{i}f_{i}^{\top}(t_{0}+j\tau)}$
$\displaystyle=$
$\displaystyle~{}\sum_{i=1}^{k}v_{i}e^{2\pi\mathbf{i}f_{i}^{\top}t_{0}}e^{2\pi\mathbf{i}f_{i}^{\top}j\tau}$
$\displaystyle=$ $\displaystyle~{}\sum_{i=1}^{k}b_{i}\cdot z_{i}^{j}.$

where the last step follows from the definition of $b_{i}$ and $z_{i}$.

Given $k$ and $z_{1},\cdots,z_{k}$, let $P(z)=\sum_{j=0}^{m}c_{j}z^{j}$ be the
degree $m$-polynomial in [CKPS16, Lemma 5.4].

$\displaystyle\sum_{j=0}^{m}c_{j}x(t_{0}+j\tau)=$
$\displaystyle~{}\sum_{j=0}^{m}c_{j}\sum_{i=1}^{k}b_{i}\cdot z_{i}^{j}$
$\displaystyle=$ $\displaystyle~{}\sum_{i=1}^{k}b_{i}\sum_{j=0}^{m}c_{j}\cdot
z_{i}^{j}$ $\displaystyle=$ $\displaystyle~{}\sum_{i=1}^{k}b_{i}P(z_{i})$
$\displaystyle=$ $\displaystyle~{}0,$

where the last step follows by Property i@ of $P(z)$ in [CKPS16, Lemma 5.4].

By Property ii@ and iii@ in [CKPS16, Lemma 5.4], we have
$x(t_{0})=-\sum_{j=1}^{m}c_{j}x(t_{0}+j\tau)$. ∎

The energy bound in Theorem 5.3 for $d$-dimensional signals is nearly optimal
due to a $k^{\Omega(d)}$ lower bound as follows.121212The proof is due to Yang
P. Liu.

######  Lemma 5.5.

Given $d\geq 1,\delta\in(0,0.1),k\in\mathbb{Z}$ such that $k\geq
O(d^{1+\delta})$. Then, there is a $d$-dimensional $k$-Fourier-sparse signal
$x(t):\mathbb{R}^{d}\rightarrow\mathbb{C}$ and a duration $T$ such that,

$\displaystyle\max_{t\in[0,T]^{d}}|x(t)|^{2}\geq k^{\Omega(\delta
d)}\|x(t)\|_{T}^{2}.$

###### Proof.

We consider the following construction of $x(t)$:

$\displaystyle x(t):=2^{-k}(1+e^{2\pi\mathbf{i}\langle f_{0},t\rangle})^{k},$

where $f_{0}={\bf 1}/(100dT)\in\mathbb{R}^{d}$.

It is easy to see that $x$ is a $k$-Fourier sparse signal, and

$\displaystyle|x(t)|^{2}=$ $\displaystyle~{}2^{-2k}|1+e^{2\pi\mathbf{i}\langle
f_{0},t\rangle}|^{2k}$ $\displaystyle=$
$\displaystyle~{}2^{-2k}\big{(}(1+e^{2\pi\mathbf{i}\langle
f_{0},t\rangle})(1+e^{-2\pi\mathbf{i}\langle f_{0},t\rangle})\big{)}^{k}$
$\displaystyle=$ $\displaystyle~{}2^{-2k}\big{(}2+e^{2\pi\mathbf{i}\langle
f_{0},t\rangle}+e^{-2\pi\mathbf{i}\langle f_{0},t\rangle}\big{)}^{k}$
$\displaystyle=$ $\displaystyle~{}2^{-2k}\cdot 2^{k}(1+\cos(2\pi\langle
f_{0},t\rangle))^{k}$ $\displaystyle=$
$\displaystyle~{}\Big{(}\frac{1+\cos(2\pi\langle
f_{0},t\rangle)}{2}\Big{)}^{k}$ $\displaystyle=$
$\displaystyle~{}\cos(\pi\langle f_{0},t\rangle)^{2k},$

where the first step follows from the definition, the second step follows from
$|z|^{2}=z\overline{z}$, the third step is straightforward, the fourth step
follows from $e^{ia}+e^{-ia}=2\cos(a)$, the fifth step is straightforward, the
last step follows from $(\cos(a)+1)/2=\cos(a/2)^{2}$.

Then, we know that

$\displaystyle\max_{t\in[0,T]^{d}}|x(t)|^{2}=|x({\bf 0})|^{2}=1.$ (9)

It remains to upper bound

$\displaystyle\|x(t)\|_{T}^{2}=T^{-d}\int_{[0,T]^{d}}|x(t)|^{2}\mathrm{d}t=T^{-d}\int_{[0,T]^{d}}\cos(\pi\langle
f_{0},t\rangle)^{2k}\mathrm{d}t.$

Let $r$ be a parameter. We have

$\displaystyle\int_{[0,T]^{d}}\cos(\pi\langle
f_{0},t\rangle)^{2k}\mathrm{d}t=$
$\displaystyle~{}\int_{W_{d}(r)}\cos(\pi\langle
f_{0},t\rangle)^{2k}\mathrm{d}t+\int_{[0,T]\backslash W_{d}(r)}\cos(\pi\langle
f_{0},t\rangle)^{2k}\mathrm{d}t$ $\displaystyle\leq$
$\displaystyle~{}\int_{W_{d}(r)}1\cdot\mathrm{d}t+T^{d}\cdot\max_{t\in[0,T]\backslash
W_{d}(r)}~{}\cos(\pi\langle f_{0},t\rangle)^{2k},$

where $W_{d}(r):=\\{t\in[0,T]^{d}~{}|~{}t_{1}+\cdots+t_{d}\leq r\\}$.

We first bound the second term:

$\displaystyle\cos(\pi\langle
f_{0},t\rangle)^{2k}\leq\exp\big{(}-\pi^{2}k\langle
f_{0},t\rangle^{2}\big{)}\leq\exp(-\Omega(kr^{2}/(dT)^{2})),$

where the first step follows from Fact 4.5, and the second step follows from
the definition of $f_{0}$ and $t\notin W_{d}(r)$. Hence,

$\displaystyle T^{d}\cdot\max_{t\in[0,T]\backslash W_{d}(r)}~{}\cos(\pi\langle
f_{0},t\rangle)^{2k}\leq T^{d}\cdot\exp(-\Omega(kr^{2}/(dT)^{2})).$ (10)

Next, we bound the first term, which is equal to the volume of $W_{d}(r)$.
Note that $W_{d}(r)$ is contained in the following simplex:

$\displaystyle
P_{d}(r):=\\{t\in\mathbb{R}_{+}^{d}~{}|~{}t_{1}+\cdots+t_{d}\leq r\\}.$

Thus, we have

$\displaystyle\mathrm{Vol}(W_{d}(r))\leq$
$\displaystyle~{}\mathrm{Vol}(P_{d}(r))\leq r^{d}\cdot\mathrm{Vol}(P_{d}(1))$
$\displaystyle\leq$ $\displaystyle~{}r^{d}/d!$ $\displaystyle=$
$\displaystyle~{}O(r/d)^{d},$ (11)

where the first step is straightforward, the second step follows from the
scaling of the volume, the third step follows from a well-known fact on the
volume of a $d$-dimensional simplex $P_{d}(1)=1/d!$ (see e.g. [Ste66]), and
the last step follows from Stirling’s approximation.

Combining Eqs. (10) and (5.2) together, we get that

$\displaystyle T^{-d}\|x(t)\|_{T}^{2}\leq
O(\frac{r}{dT})^{d}+\exp(-\Omega(kr^{2}/(dT)^{2})).$

By taking $r=Tdk^{-\delta/6}$, we have

$\displaystyle O(\frac{r}{dT})^{d}=O(k^{-\delta/6})^{d}=k^{-\Omega(\delta
d)},$

and

$\displaystyle\exp(-\Omega(kr^{2}/(dT)^{2}))=$
$\displaystyle~{}\exp(-\Omega(k^{1-\delta/3}))$ $\displaystyle=$
$\displaystyle~{}\exp(-\Omega(d^{(1+\delta)(1-\delta/3)}))$
$\displaystyle\leq$ $\displaystyle~{}\exp(-\Omega(\delta d\log k))$
$\displaystyle=$ $\displaystyle~{}k^{-\Omega(d)},$

where the first step is straightforward, the second step follows from
$k=O(d^{1+\delta})$, and the last step follows from $\delta\in(0,1)$.

Therefore, we have

$\displaystyle T^{-d}\|x(t)\|_{T}^{2}\leq
k^{-\Omega(d)}=k^{-\Omega(d)}\cdot\max_{t\in[0,T]^{d}}~{}|x(t)|^{2},$

where the last step follows from Eq. (9).

The lemma is then proved.

∎

###  5.3 Energy bound for discrete Fourier signals

Recall the definition of one-dimensional discrete sparse Fourier signal: for
$t\in\\{0,1,\dots,n-1\\}$,

$\displaystyle
x_{t}=\frac{1}{n}\sum_{i=1}^{k}\widehat{x}_{f_{i}}\exp\left(\frac{2\pi\mathbf{i}}{n}f_{i}t\right),~{}f_{i}\in[n]$
(12)

More generally, for $d\geq 1$, the $d$-dimensional discrete sparse Fourier
signal can be defined as follows. Let $n=p^{d}$ where both $p$ and $d$ are
positive integers. Recall the definition of high-dimensional discrete sparse
Fourier signal (see e.g. [NSW19]):

$\displaystyle
x_{t}=\frac{1}{n}\sum_{i=1}^{k}\widehat{x}_{f_{i}}\exp\left(\frac{2\pi\mathbf{i}}{p}\langle
f_{i},t\rangle\right)~{}~{}~{}\forall t\in[p]^{d},$ (13)

where each $f_{i}\in[p]^{d}$.

In this section, we prove the following discrete Fourier signals energy bound
that works for any dimension:

######  Theorem 5.6 (Discrete $d$-dimensional Fourier energy bound).

For $d\geq 1$ and any discrete $d$-dimensional $k$-sparse Fourier signal
$\\{x_{i}\\}_{i\in[n]}$, we have

$\displaystyle\|x\|_{\infty}^{2}\leq k\cdot\frac{\|x\|_{2}^{2}}{n}$

###### Proof.

$\displaystyle\|x\|_{\infty}^{2}\leq\frac{1}{n^{2}}\|\widehat{x}\|_{1}^{2}\leq\frac{k}{n^{2}}\|\widehat{x}\|_{2}^{2}=\frac{k}{n}\|x\|_{2}^{2},$

where the first step follows from Claim 5.7, the second step follows from
Cauchy-Schwarz inequality, and the last step follows from Theorem 5.8. ∎

######  Claim 5.7.

For any $d\geq 1$ and any discrete $d$-dimensional Fourier signal $x$,

$\displaystyle\|x\|_{\infty}\leq\frac{1}{n}\|\widehat{x}\|_{1}.$

###### Proof.

By triangle inequality,

$\displaystyle|x(t)|\leq\frac{1}{n}\left|\sum_{i=1}^{k}\widehat{x}_{i}\exp\left(\frac{2\pi\mathbf{i}}{p}\langle
f_{i},t\rangle\right)\right|\leq\frac{1}{n}\sum_{i=1}^{k}|\widehat{x}_{i}|=\frac{1}{n}\|\widehat{x}\|_{1}.$

∎

######  Theorem 5.8 (Parseval’s theorem).

For any $d\geq 1$ and any discrete $d$-dimensional Fourier signal (Eq. (13)),

$\displaystyle\|x\|_{2}^{2}=\frac{1}{n}\|\widehat{x}\|_{2}^{2}.$

###  5.4 Energy bounds imply concentrations

By using Chernoff bound, we prove the following lemma to show the performance
of uniformly sampling.

####  5.4.1 Continuous case

######  Lemma 5.9.

Let $d\in\mathbb{Z}_{+}$. Let $R$ be a parameter. Given any function
$x(t):\mathbb{R}^{d}\rightarrow\mathbb{C}$ with
$\underset{t\in[0,T]^{d}}{\max}|x(t)|^{2}\leq R\|x(t)\|_{T}^{2}$. Let $S$
denote a set of points chosen uniformly at random from $[0,T]^{d}$. We have
that

$\displaystyle\mathsf{Pr}\left[\left|\frac{1}{|S|}\sum_{i\in
S}|x(t_{i})|^{2}-\|x(t)\|_{T}^{2}]\right|\geq\varepsilon\|x(t)\|_{T}^{2}\right]$
$\displaystyle\leq$ $\displaystyle\exp(-\Omega(\varepsilon^{2}|S|/R)),$

where $\|x(t)\|_{T}^{2}=\frac{1}{T^{d}}\int_{[0,T]^{d}}|x(t)|^{2}\mathrm{d}t$.

###### Proof.

Let $M$ denote $\underset{t\in[0,T]^{d}}{\max}|x(t)|^{2}$. Replacing $X_{i}$
by $\frac{|x(t_{i})|^{2}}{M}$ and $n$ by $|S|$ in Lemma 4.1, we obtain that

$\displaystyle\mathsf{Pr}[|X-\mu|>\varepsilon\mu]\leq
2\exp(-\frac{\varepsilon^{2}}{3}\mu)$

The above equation implies

$\displaystyle\mathsf{Pr}\left[\left|\sum_{i\in
S}\frac{|x(t_{i})|^{2}}{M}-|S|\frac{\|x(t)\|_{T}^{2}}{M}\right|>\varepsilon|S|\frac{\|x(t)\|_{T}^{2}}{M}\right]\leq
2\exp(-\frac{\varepsilon^{2}}{3}\mu)$

Multiplying $M$ on the both sides

$\displaystyle\mathsf{Pr}\left[\left|\frac{1}{|S|}\sum_{i\in
S}|x(t_{i})|^{2}-\|x(t)\|_{T}^{2}\right|\geq\varepsilon\|x(t)\|_{T}^{2}\right]\leq
2\exp(-\frac{\varepsilon^{2}}{3}\mu)$

Applying bound on $\mu$

$\displaystyle\mathsf{Pr}\left[\left|\frac{1}{|S|}\sum_{i\in
S}|x(t_{i})|^{2}-\|x(t)\|_{T}^{2}\right|\geq\varepsilon\|x(t)\|_{T}^{2}\right]\leq
2\exp(-\frac{\varepsilon^{2}}{3}|S|\frac{\|x(t)\|_{T}^{2}}{M})$

which is less than $2\exp(-\frac{\varepsilon^{2}}{3}|S|/R)$, thus completes
the proof. ∎

####  5.4.2 Discrete case

######  Lemma 5.10.

Let $R$ be a parameter. Given any function $x\in\mathbb{C}^{n}$ with
$\|x\|_{\infty}^{2}\leq R\|x\|_{2}^{2}/n$. Let $S$ denote a set of points
chosen uniformly at random from $[n]$. We have that

$\displaystyle\mathsf{Pr}\left[\left|\frac{1}{|S|}\sum_{t\in
S}|x_{t}|^{2}-n^{-1}\|x\|_{2}^{2}]\right|\geq\varepsilon
n^{-1}\|x\|_{2}^{2}\right]$ $\displaystyle\leq$
$\displaystyle\exp(-\Omega(\varepsilon^{2}|S|/R)).$

###### Proof.

Let $M$ denote $\underset{t\in[n]}{\max}|x_{t}|^{2}$. Replacing $X_{i}$ by
$\frac{|x_{t}|^{2}}{M}$ and $n$ by $|S|$ in Lemma 4.1, we obtain that

$\displaystyle\mathsf{Pr}[|X-\mu|>\varepsilon\mu]\leq
2\exp(-\frac{\varepsilon^{2}}{3}\mu)$

The above equation implies that

$\displaystyle\mathsf{Pr}\left[\left|\sum_{t\in
S}\frac{|x_{t}|^{2}}{M}-|S|\frac{\|x\|_{2}^{2}}{nM}\right|>\varepsilon|S|\frac{\|x\|_{2}^{2}}{nM}\right]\leq
2\exp(-\frac{\varepsilon^{2}}{3}\mu)$

Multiplying the normalization factor on both sides,

$\displaystyle\mathsf{Pr}\left[\left|\frac{1}{|S|}\sum_{t\in
S}|x_{t}|^{2}-n^{-1}\|x\|_{2}^{2}\right|\geq\varepsilon
n^{-1}\|x\|_{2}^{2}\right]\leq 2\exp(-\frac{\varepsilon^{2}}{3}\mu)$

Applying bound on $\mu$

$\displaystyle\mathsf{Pr}\left[\left|\frac{1}{|S|}\sum_{t\in
S}|x_{t}|^{2}-n^{-1}\|x\|_{2}^{2}\right|\geq\varepsilon
n^{-1}\|x\|_{2}^{2}\right]\leq
2\exp(-\frac{\varepsilon^{2}}{3}|S|\frac{\|x\|_{2}^{2}}{nM})$

which is less than $2\exp(-\frac{\varepsilon^{2}}{3}|S|/R)$, thus completes
the proof. ∎

##  6 Oblivious Sketching Fourier Sparse Signals

In this section, we show an intermediate step in the reduction from Frequency
estimation to Signal estimation: constructing a small sketching subset $S$ of
the time domain _obliviously_ (without making any query to the signal), so
that any signal in the family ${\cal F}$ discretized by $S$ has norm close to
the original continuous signal. More formally, we define the _oblivious
sketching Fourier signal problem_ as follows:

######  Problem 6.1 (Oblivious sketching Fourier sparse signal problem).

Let ${\cal F}$ be a family of Fourier-sparse signals. Let
$\varepsilon\in(0,0.1)$ denote the accuracy parameter. The goal is to find a
set $S=\\{t_{1},\dots,t_{s}\\}$ in the time domain of size $s$ such that, with
high probability, for any signal $x\in{\cal F}$, it holds that

$\displaystyle(1-\varepsilon)\|x\|_{T}\leq\|x\|_{S}\leq(1+\varepsilon)\|x\|_{T},$

where

$\displaystyle\|x\|_{T}^{2}:=\frac{1}{T^{d}}\int_{[0,T]^{d}}|x(t)|^{2}\mathrm{d}t,~{}~{}\mathrm{and}~{}~{}\|x\|^{2}_{S}:=\frac{1}{|S|}\sum_{i\in[s]}|x(t_{i})|^{2}.$

We remark that the concentration inequalities in Section 5.4 do not give a
small sketching directly. To apply concentration inequalities for all signals
in ${\cal F}$, we need an $\varepsilon$-net, which will increase the sketching
size by a factor of $\mathrm{poly}(k)$.

In Section 6.1, we show a unified approach to obtain oblivious sketching using
uniform samples for general signals. Then, in Section 6.2, we show how to
sketch one-dimensional signals with nearly-optimal size using importance
sampling.

###  6.1 A unified approach via uniform sampling

In this sub-section, we prove a general result of oblivious sketching for
discrete and continuous Fourier-sparse signals. The technical tools we employ
are Fourier energy bounds and importance sampling in Section 4.4.

######  Lemma 6.2 (Oblivious sketching via uniform sampling).

Let ${\cal F}\subset\\{x(t):{\cal G}\rightarrow\mathbb{C}\\}$ be a family of
signals of dimension $k$. Suppose the energy bound $R\geq\sup_{t\in{\cal
G}}\sup_{x\in{\cal F}}\frac{|x(t)|}{\|x(t)\|_{T}}$ holds. For any
$\varepsilon,\rho\in(0,1)$, let $S=\\{t_{1},t_{2},\dots,t_{s}\\}$ be a set of
i.i.d. samples chosen uniformly at random over the time domain ${\cal G}$ of
size

$\displaystyle s\geq\Omega\left(\varepsilon^{-2}R\log(k/\rho)\right).$

Then, with probability at least $1-\rho$, it holds that for all $x\in{\cal
F}$,

$\displaystyle(1-\varepsilon)\|x\|_{T}\leq\|x\|_{S}\leq(1+\varepsilon)\|x\|_{T}.$

###### Proof.

To show that sampling from the distribution $D=\mathrm{Uniform}({\cal G})$
gives a good sketch, we apply Lemma 4.18 with $D^{\prime}=D$, $d=k$,
$\delta=\rho$, $w_{i}=1/s$, which implies that the matrix
$A\in\mathbb{C}^{s\times k}$ defined by $A_{i,j}:=\sqrt{w_{i}}\cdot
v_{j}(t_{i})$ satisfies:

$\displaystyle\|A^{*}A-I\|_{2}\leq\varepsilon$ (14)

with probability at least $1-\rho$, as long as
$s=\Omega(\varepsilon^{-2}K_{\mathsf{IS},D}\log(k/{\rho}))$.

Under the condition that Eq. (14) holds, Lemma 4.17 shows that

$\displaystyle(1-\varepsilon)\|x\|_{T}^{2}\leq\|{x}\|^{2}_{S}\leq(1+\varepsilon)\|x\|_{T}^{2}~{}~{}~{}\forall
x\in{\cal F},$

which proves the oblivious sketching guarantee.

It remains to upper-bound the condition number $K_{\mathsf{IS},D}$. By
Definition 4.16, we have

$\displaystyle K_{\mathsf{IS},D}=$ $\displaystyle~{}\underset{t\in{\cal
G}}{\sup}~{}\underset{x\in\mathcal{F}}{\sup}~{}\frac{|x(t)|^{2}}{\|x\|_{T}^{2}}\leq
R.$

Thus, we get that

$\displaystyle s=\Omega\left(\varepsilon^{-2}R\log(k/\rho)\right).$

The lemma is then proved. ∎

Using the $d$-dimensional Fourier-sparse signals’ energy bound (Theorem 5.3)
and the discrete Fourier-sparse signals’ energy bound (Theorem 5.6), we
immediately obtain the following corollaries:

######  Corollary 6.3 (Oblivious sketching high-dimensional continuous
signal).

Let $d>1$ be the dimension of the signal. Let
$f_{1},\dots,f_{k}\in\mathbb{R}^{d}$. For any $\varepsilon\in(0,1)$, let
$S_{d}$ be a set of i.i.d. samples chosen uniformly at random over $[0,T]^{d}$
of size $|S_{d}|\geq\varepsilon^{-2}k^{O(d)}\log(k/\rho)$. Let
$V:=\\{\exp(2\pi\mathbf{i}\langle f_{i},t\rangle)~{}|~{}i\in[k]\\}$. Then,
with probability at least $1-\rho$, for all $x\in\mathrm{span}\\{V\\}$, we
have

$(1-\varepsilon)\|x\|_{T}\leq\|{x}\|_{S_{d}}\leq(1+\varepsilon)\|{x}\|_{T}.$

######  Corollary 6.4 (Oblivious sketching discrete signal).

For $d\geq 1$, let $n=p^{d}$ for some positive integer $p$. Let
$k\in\mathbb{N}_{+}$ and $f_{1},\dots,f_{k}\in[p]^{d}$. Define
$V:=\left\\{(e^{2\pi\mathbf{i}\langle
f_{i},t\rangle/p})_{t\in[p]^{d}}~{}{|}~{}\forall
i\in[k]\right\\}\subseteq\mathbb{C}^{[p]^{d}}$. For any
$\varepsilon,\rho\in(0,1)$, let $S$ be a set of i.i.d. samples chosen
uniformly at random over $[n]$ of size $|S|\geq
O(\varepsilon^{-2}k\log(k/\rho))$. Then, with probability at least $1-\rho$,
for all $u\in\mathrm{span}\\{V\\}$, we have

$(1-\varepsilon)\|u\|^{2}_{2}\leq
n\|{u}\|^{2}_{S}\leq(1+\varepsilon)\|{u}\|^{2}_{2},$

where $\|u\|_{S}^{2}=\sum_{i\in S}|u_{i}|^{2}/|S|$.

###  6.2 Weighted oblivious sketching one-dimensional signals

For one-dimensional signals, using uniform sampling (Lemma 6.2) and the energy
bound (Theorem 5.1) gives a sketching set of size at least
$\varepsilon^{-2}k^{2}$. Motivated by [CP19a], we show a more efficient
sketching method using only $\varepsilon^{-2}k$ samples by assigning different
weights to each sample point. In the following lemma, we give a sketch for any
one-dimensional Fourier sparse signal with nearly-optimal size:

######  Lemma 6.5 (Nearly-optimal weighted sketch for one-dimensional
signals).

For $k\in\mathbb{N}_{+}$, define a probability distribution ${\cal D}(t)$ as
follows:

$\displaystyle{\cal D}(t):=\begin{cases}{c}/(1-|t/T|),&\text{ for }|t|\leq
T(1-{1}/k)\\\ c\cdot k,&\text{ for }|t|\in[T(1-{1}/k),T]\end{cases}$ (15)

where $c=\Theta(T^{-1}\log^{-1}(k))$ is a normalization factor such that
$\int_{-T}^{T}{\cal D}(t)\mathrm{d}t=1$.

For any $f_{1},\cdots,f_{k}\in[-F,F]$, let ${\cal
F}:=\Big{\\{}x(t)=\sum_{j=1}^{k}v_{j}\cdot
e^{2\pi\mathbf{i}f_{j}t}~{}\Big{|}~{}v_{j}\in\mathbb{C}\Big{\\}}$. For any
$\varepsilon,\rho\in(0,1)$, let $S_{{\cal D}}=\\{t_{1},\cdots,t_{s}\\}$ be a
set of i.i.d. samples from ${\cal D}(t)$ of size $s\geq
O(\varepsilon^{-2}k\log(k)\log(k/\rho))$. Let the weight vector
$w\in\mathbb{R}^{s}$ be defined by $w_{i}:=2/(Ts{\cal D}(t_{i}))$ for
$i\in[s]$. Then with probability at least $1-\rho$, we have

$(1-\varepsilon)\|x\|_{T}\leq\|x\|_{S_{\cal
D},w}\leq(1+\varepsilon)\|x\|_{T},$

where $\|x\|_{T}^{2}:=\frac{1}{2T}\int_{-T}^{T}|x(t)|^{2}\mathrm{d}t$.131313We
use time duration $[-T,T]$ for convenience. It is easy to transform to $[0,T]$
by shifting and re-scaling.

###### Proof.

Let $\\{v_{1}(t),v_{2}(t),\cdots,v_{k}(t)\\}$ be an orthonormal basis for
${\cal F}$ with respect to the distribution ${\cal D}$, i.e.,

$\displaystyle\int_{0}^{T}{\cal D}(t)\cdot
v_{i}(t)\overline{v_{j}(t)}\mathrm{d}t=~{}{\bf 1}_{i=j},\quad\forall
i,j\in[k].$

We first prove that the distribution ${\cal D}$ is well-defined. By the
condition that $\int_{-T}^{T}{\cal D}(t)\mathrm{d}t=1$, we have

$\displaystyle
2\int_{0}^{T(1-{1}/({k}))}\frac{c}{(1-|t/T|)}\mathrm{d}t+2\int_{T(1-{1}/({k}))}^{T}c\cdot
k^{2}k\mathrm{d}t=1,$

which implies that

$\displaystyle c^{-1}=$
$\displaystyle~{}2\int_{0}^{T(1-{1}/k)}\frac{1}{(1-|t/T|)}\mathrm{d}t+2\int_{T(1-{1}/k)}^{T}k^{2}\mathrm{d}t$
$\displaystyle=$ $\displaystyle~{}2T{\log k}+2T$ $\displaystyle=$
$\displaystyle~{}\Theta(T\log(k)).$

Thus, we get that $c=\Theta(T^{-1}\log^{-1}(k))$.

To show that sampling from distribution ${\cal D}$ gives a good weighted
sketch, we will use some technical tools in Section 4.4. Applying Lemma 4.18
with $D^{\prime}={\cal D}$, $D=\mathrm{Uniform}([-T,T])$, $d=k$,
$\delta=\rho$, we have that, with probability at least $1-\rho$, the matrix
$A\in\mathbb{C}^{s\times k}$ defined by $A_{i,j}:=\sqrt{w_{i}}\cdot
v_{j}(t_{i})$ satisfying

$\displaystyle\|A^{*}A-I\|_{2}\leq\varepsilon,$

as long as
$s\geq\Omega(\varepsilon^{-2}K_{\mathsf{IS},D^{\prime}}\log(k/{\rho}))$, where
$w_{i}=\frac{D(t_{i})}{s\cdot D^{\prime}(t_{i})}=\frac{2}{Ts{\cal D}(t_{i})}$.
Then, by Lemma 4.17, it implies that for every $x\in{\mathcal{F}}$,

$\displaystyle(1-\varepsilon)\|x\|_{T}^{2}\leq\|{x}\|^{2}_{S_{\cal
D},w}\leq(1+\varepsilon)\|x\|_{T}^{2}.$

It remains to upper-bound the condition number $K_{\mathsf{IS},D^{\prime}}$
(see Definition 4.16):

$\displaystyle K_{\mathsf{IS},D^{\prime}}:=$
$\displaystyle~{}\underset{t}{\sup}\\{\frac{D(t)}{D^{\prime}(t)}\cdot\underset{f\in\mathcal{F}}{\sup}\\{\frac{|f(t)|^{2}}{\|f\|_{D}^{2}}\\}\\}$
$\displaystyle=$ $\displaystyle~{}\underset{t}{\sup}\\{\frac{1}{2T{\cal
D}(t)}\cdot\underset{f\in\mathcal{F}}{\sup}\\{\frac{|f(t)|^{2}}{\|f\|_{T}^{2}}\\}\\}$
$\displaystyle\leq$ $\displaystyle~{}\underset{t}{\sup}\\{\frac{1}{2T{\cal
D}(t)}\cdot{\min}\\{\frac{k}{1-|t/T|},k^{2}\\}\\}$ $\displaystyle\leq$
$\displaystyle~{}{\max}\\{\frac{(1-|t/T|)}{2cT}\frac{k}{1-|t/T|},\frac{1}{2cTk}k^{2}\\}$
$\displaystyle=$ $\displaystyle~{}\frac{k}{2cT}$ $\displaystyle=$
$\displaystyle~{}O(k\log k),$

where the first step follows from the definition, the second step follows from
$D(t)=\mathrm{Uniform}([-T,T])(t)=\frac{1}{2T}$, the third step follows from
Theorem 5.1 and Theorem 5.2, and the remaining steps follow from direct
calculations. Thus, we get that

$\displaystyle s\geq\Omega\left(\varepsilon^{-2}k\log(k)\log(k/\rho)\right).$

The lemma is then proved. ∎

##  7 Fast Implementation of Well-Balanced Sampling Procedure

Well-balanced sampling procedure was first defined in [CP19a] to study the
active linear regression problem. Our signal estimation algorithm will call it
as a sub-procedure. In this section, we give a fast implementation of well-
balanced sampling procedure based on the Randomized BSS algorithm [BSS12,
LS15].

First, we restate the definition of well-balanced sampling procedure in
[CP19a].

######  Definition 7.1 (Well-balanced sampling procedure (WBSP), [CP19a]).

Given a linear family $\mathcal{F}$ and underlying distribution $D$, let $P$
be a random sampling procedure that terminates in $m$ iterations ($m$ is not
necessarily fixed) and provides a coefficient $\alpha_{i}$ and a distribution
$D_{i}$ to sample $x_{i}\sim D_{i}$ in every iteration $i\in[m]$.

We say $P$ is an $\varepsilon$-WBSP if it satisfies the following two
properties:

  1. 1.

With probability $0.9$, for weight
$w_{i}=\alpha_{i}\cdot\frac{D(x_{i})}{D_{i}(x_{i})}$ of each $i\in[m]$,

$\sum_{i=1}^{m}w_{i}\cdot|h(x_{i})|^{2}\in\left[1-10\sqrt{\varepsilon},1+10\sqrt{\varepsilon}\right]\cdot\|h\|_{D}^{2}\quad\forall
h\in\mathcal{F}.$

  2. 2.

The coefficients always have $\sum_{i=1}^{m}\alpha_{i}\leq\frac{5}{4}$ and
$\alpha_{i}\cdot K_{\mathsf{IS},D_{i}}\leq\frac{\varepsilon}{2}$ for all
$i\in[m]$.

This definition describes a general sampling procedure that uses a few samples
to represent the whole continuous signal, and the sampling procedure should
satisfy two properties: one guarantees that the norm of any function in a
function family is preserved, and another guarantees that the norm of noise is
also preserved.

In Section 7.1, we review some results in [CP19a] and show that WBSP can be
implemented via randomized spectral sparsification. In Section 7.2, we design
a data structure and improve the time efficiency of the WBSP. In Section 7.3,
we discover a tradeoff between the preprocessing cost and the query cost,
which can improve the space complexity.

###  7.1 Randomized BSS implies a WBSP

In this section, we review the result of [CP19a], which shows that the
Randomized BSS algorithm [BSS12, LS15] implies a well-balanced sampling
procedure.

######  Lemma 7.2 (Lemma 5.1 in [CP19a]).

Let $G$ be any domain. Given any dimension $d$ linear function family
$\mathcal{F}$ of function $f:G\rightarrow\mathbb{C}$,

$\displaystyle{\mathcal{F}}=\\{f(t)=\sum_{j=1}^{d}v_{j}u_{j}(t)|v_{j}\in\mathbb{C}\\},$

where $u_{j}:G\rightarrow\mathbb{C}$. Given any distribution $D$ over $G$, and
any $\varepsilon>0$, there exists an efficient procedure (Algorithm 2) that
runs in $O(\varepsilon^{-1}d^{3}|G|+\varepsilon^{-1}d^{\omega+1})$ time and
outputs a set $S\subseteq G$ and weight $w$ such that

  * •

$|S|=O(d/\varepsilon),~{}w\in\mathbb{R}^{|S|}$,

  * •

the procedure is an $\varepsilon$-WBSP,

holds with probability $1-\frac{1}{200}$.

Algorithm 2 A well-balanced sampling procedure based on Randomized BSS (see
[CP19a])

1:procedure RandBSS($d,\mathcal{F},D,\varepsilon$)

2: Find an orthonormal basis $v_{1},\ldots,v_{d}$ of $\mathcal{F}$ under $D$

3: Set $\gamma\leftarrow\sqrt{\varepsilon}/3$ and
$\textsf{mid}\leftarrow\frac{4d/\gamma}{1/(1-\gamma)-1/(1+\gamma)}$

4: $j\leftarrow 0,B_{0}\leftarrow 0$

5: $l_{0}\leftarrow-2d/\gamma,u_{0}\leftarrow 2d/\gamma$

6: while $u_{j+1}-l_{j+1}<8d/\gamma$ do

7:
$\Phi_{j}\leftarrow\mathrm{tr}[(u_{j}I-B_{j})^{-1}]+\mathrm{tr}[(B_{j}-l_{j}I)^{-1}]$
$\triangleright$ The potential function at iteration $j$.

8: Set the coefficient
$\alpha_{j}\leftarrow\frac{\gamma}{\Phi_{j}}\cdot\frac{1}{\textsf{mid}}$

9: Set $v(x)\leftarrow\big{(}v_{1}(x),\ldots,v_{d}(x)\big{)}$

10: for $x\in\mathrm{supp}(D)$ do

11: Set the distribution  $D_{j}(x)\leftarrow
D(x)\cdot\bigg{(}v(x)^{\top}(u_{j}I-B_{j})^{-1}v(x)+v(x)^{\top}(B_{j}-l_{j}I)^{-1}v(x)\bigg{)}/\Phi_{j}$

12: end for

13: Sample $x_{j}\sim D_{j}$ and set a scale
$s_{j}\leftarrow\frac{\gamma}{\Phi_{j}}\cdot\frac{D(x_{j})}{D_{j}(x_{j})}$

14: $B_{j+1}\leftarrow B_{j}+s_{j}\cdot v(x_{j})v(x_{j})^{\top}$

15: $u_{j+1}\leftarrow u_{j}+\frac{\gamma}{\Phi_{j}(1-\gamma)},\quad
l_{j+1}\leftarrow l_{j}+\frac{\gamma}{\Phi_{j}(1+\gamma)}$

16: $j\leftarrow j+1$

17: end while

18: $m\leftarrow j$

19: Assign the weight $w_{j}\leftarrow s_{j}/\textsf{mid}$ for each $x_{j}$

20: return $\\{x_{1},x_{2},\cdots,x_{m}\\},w$

21:end procedure

###  7.2 Fast implementation of WBSP

In this section, we give a fast implementation of Algorithm 2:

######  Theorem 7.3 (Fast implementation of WBSP).

Let $G$ be any domain. Given any dimension $d$ linear function family
$\mathcal{F}$ of function $f:G\rightarrow\mathbb{C}$,

$\displaystyle{\mathcal{F}}=\\{f(t)=\sum_{j=1}^{d}v_{j}u_{j}(t)|v_{j}\in\mathbb{C}\\},$

where $u_{j}:G\rightarrow\mathbb{C}$. Given any distribution $D$ over $G$, and
any $\varepsilon>0$, there exists an efficient procedure (Algorithm 3) that
runs in
$O(d^{2}|G|+\varepsilon^{-1}d^{3}\log|G|+\varepsilon^{-1}d^{\omega+1})$ time
and outputs a set $S\subseteq G$ and weight $w\in\mathbb{R}^{|S|}$ such that
the following properties hold with probability at least 0.995:

  * •

$|S|=O(d/\varepsilon)$,

  * •

the procedure is an $\varepsilon$-WBSP.

Algorithm 3 Our fast implementation of well-balanced sampling procedure

1:procedure RandBSS+($d,\mathcal{F},D,\varepsilon$)$\triangleright$ Theorem
7.3

2: /*Preprocessing*/

3: Find an orthonormal basis $v_{1},\ldots,v_{d}$ of $\mathcal{F}$ under $D$

4: $\gamma\leftarrow\sqrt{\varepsilon}/3$ and
$\textsf{mid}\leftarrow\frac{4d/\gamma}{1/(1-\gamma)-1/(1+\gamma)}$

5: $j\leftarrow 0,B_{0}\leftarrow 0$

6: $l_{0}\leftarrow-2d/\gamma,u_{0}\leftarrow 2d/\gamma$

7: $\delta\leftarrow 1/\mathrm{poly}(d)$

8: $\triangleright$ Let
$v(x)=\big{(}v_{1}(x),\ldots,v_{d}(x)\big{)}\in\mathbb{R}^{d}$

9:
$\mathrm{DS}.\textsc{Init}(|D|,d,\\{v(x_{1}),\cdots,v(x_{|D|})\\}\subset\mathbb{R}^{d},\\{D(x_{1}),\dots,D(x_{|D|})\\}\subset\mathbb{R})$
$\triangleright$ Algorithm 4

10: /*Iterative step*/

11: while $u_{j+1}-l_{j+1}<8d/\gamma$ do

12:
$\Phi_{j}\leftarrow\mathrm{tr}[(u_{j}I-B_{j})^{-1}]+\mathrm{tr}[(B_{j}-l_{j}I)^{-1}]$
$\triangleright$ The potential function at iteration $j$.

13: $\alpha_{j}\leftarrow\frac{\gamma}{\Phi_{j}}\cdot\frac{1}{\textsf{mid}}$

14: $E_{j}\leftarrow(u_{j}I-B_{j})^{-1}+(B_{j}-l_{j}I)^{-1}$

15: $q\leftarrow\mathrm{DS}.\textsc{Query}(E_{j}/\Phi_{j})$ $\triangleright$
$q\in[|D|]$, Algorithm 4

16: ${\sf x}_{j}\leftarrow x_{q}$ and set a scale
$s_{j}\leftarrow\frac{\gamma}{v({\sf x}_{j})^{\top}E_{j}v({\sf x}_{j})}$

17: $B_{j+1}\leftarrow B_{j}+s_{j}\cdot v({\sf x}_{j})v({\sf x}_{j})^{\top}$

18: $u_{j+1}\leftarrow u_{j}+\frac{\gamma}{\Phi_{j}(1-\gamma)},\quad
l_{j+1}\leftarrow l_{j}+\frac{\gamma}{\Phi_{j}(1+\gamma)}$

19: $j\leftarrow j+1$

20: end while

21: $m\leftarrow j$

22: Assign the weight $w_{j}\leftarrow s_{j}/\textsf{mid}$ for each $x_{j}$

23: return $\\{{\sf x}_{1},{\sf x}_{2},\cdots,{\sf x}_{m}\\},w$

24:end procedure

Our algorithm is based on a data structure for solving the _online quadratic-
form sampling problem_ defined as follows:

######  Problem 7.4 (Online Quadratic-Form Sampling Problem).

Given $n$ vectors $v_{1},\dots,v_{n}\in\mathbb{R}^{d}$ and $n$ coefficients
$\alpha_{1},\dots,\alpha_{n}$, for any PSD matrix $A\in\mathbb{R}^{d\times
d}$, output a sample $i\in[n]$ from the following distribution ${\cal D}_{A}$:

$\displaystyle\Pr_{{\cal D}_{A}}[i]:=\frac{\alpha_{i}\cdot
v_{i}^{\top}Av_{i}}{\sum_{j=1}^{n}\alpha_{j}\cdot
v_{j}^{\top}Av_{j}}~{}~{}~{}\forall i\in[n].$ (16)

######  Theorem 7.5.

There is a data structure (Algorithm 4) that uses $O(nd^{2})$ spaces for the
Online Quadratic-Form Sampling Problem with the following procedures:

  * •

Init$(n,d,\\{v_{1},\dots,v_{n}\\}\subset\mathbb{R}^{d},\\{\alpha_{1},\dots,\alpha_{n}\\}\subset\mathbb{R})$:
the data structure preprocesses in time $O(nd^{2})$.

  * •

Query$(A\in\mathbb{R}^{d\times d})$: Given a PSD matrix $A$, the Query
operation samples $i\in[n]$ exactly from the probability distribution ${\cal
D}_{A}$ defined in Problem 7.4 in $O(d^{2}\log n)$-time.

###### Proof.

The pseudo-code of the algorithm is given as Algorithm 4. The idea is to build
a binary tree such that each node has an interval in
$[l,\dots,r]\subset[1,\dots,n]$ and stores a matrix
$\sum_{i=l}^{r}\alpha_{i}\cdot v_{i}v_{i}^{\top}$. For each internal node with
interval $[l,\dots,r]$, its left child node has interval
$[l,\dots,\lfloor(l+r)/2\rfloor]$, and its right child node has interval
$[\lfloor(l+r)/2\rfloor+1,\dots,r]$.

We first prove the correctness. Suppose the output of Query is $i\in[n]$. We
compute its probability. Let $u_{0}=\mathsf{root},u_{1},\dots,u_{t}$ be the
path from the root of the tree to the leaf with $\mathrm{id}=i$. Then, we have

$\displaystyle\Pr[u_{t}]=\prod_{j=1}^{t}\Pr[u_{j}|u_{j-1}]=\prod_{j=1}^{t}\frac{\sum_{k=l_{j}}^{r_{j}}\alpha_{k}\cdot
v_{k}^{\top}Av_{k}}{\sum_{k=l_{j-1}}^{r_{j-1}}\alpha_{k}\cdot
v_{k}^{\top}Av_{k}}=\frac{\alpha_{i}\cdot
v_{i}^{\top}Av_{i}}{\sum_{k=1}^{n}\alpha_{k}\cdot v_{k}^{\top}Av_{k}},$

where $[l_{j},\dots,r_{j}]$ is the range of the node $u_{j}$, the first step
follows from the conditional probability, the second step follows from Line 34
in Algorithm 4, and the last step follows from the telescoping products.
Hence, we get that

$\displaystyle\Pr[\textsc{Query}(A)=i]=\Pr_{{\cal D}_{A}}[i]~{}~{}~{}\forall
i\in[n].$

Hence, the sampling distribution is the same as the Online Quadratic-Form
Sampling Problem’s distribution.

For the running time, in the preprocessing stage, we build the binary tree
recursively. It is easy to see that the number of nodes in the tree is $O(n)$
and the depth is $O(\log n)$. For a leaf node, we take $O(d^{2})$-time to
compute the matrix $\alpha_{i}\cdot v_{i}v_{i}^{\top}\in\mathbb{R}^{d\times
d}$. For an internal node, we take $O(d^{2})$-time to add up the matrices of
its left and right children. Thus, the total preprocessing time is
$O(nd^{2})$.

In the query stage, we walk along a path from the root to a leaf, which has
$O(\log n)$ steps. In each step, we compute the inner product between $A$ and
the current node’s matrix, which takes $O(d^{2})$-time. And we compute the
inner product between $A$ and its left child node’s matrix, which also takes
$O(d^{2})$-time. Then, we toss a coin and decide which subtree to move. Hence,
each query takes $O(d^{2}\log n)$-time.

The theorem is then proved. ∎

Algorithm 4 Quadratic-form sampling data structure

1:structure Node

2: $V\in\mathbb{R}^{d\times d}$

3: $\mathrm{left},\mathrm{right}$ $\triangleright$ Point to the left/right
child in the tree

4:end structure

5:data structure DS

6:members

7: $n\in\mathbb{N}$$\triangleright$ The number of vectors

8: $v_{1},\dots,v_{n}\in\mathbb{R}^{d}$ $\triangleright$ $d$-dimensional
vectors

9: $\alpha_{1},\dots,\alpha_{n}\in\mathbb{R}$ $\triangleright$ Coefficients

10: $\mathsf{root}$: Node $\triangleright$ The root of the tree

11:end members

12:procedure BuildTree($l,r$)$\triangleright$ $[l,\dots,r]$ is the range of
the current node

13: ${\sf p}\leftarrow\mathbf{new}$ Node

14: if $l=r$ then$\triangleright$ Leaf node

15: ${\sf p}.V\leftarrow\alpha_{l}\cdot v_{l}v_{l}^{\top}$ $\triangleright$ It
takes $O(d^{2})$-time

16: else$\triangleright$ Internal node

17: $mid\leftarrow\lfloor(l+r)/2\rfloor$

18: ${\sf p}.\mathrm{left}\leftarrow\textsc{BuildTree}(l,mid)$

19: ${\sf p}.\mathrm{right}\leftarrow\textsc{BuildTree}(mid+1,r)$

20: ${\sf p}.V\leftarrow({\sf p}.\mathrm{left}).V+({\sf p}.\mathrm{right}).V$
$\triangleright$ It takes $O(d^{2})$-time

21: end if

22: return ${\sf p}$

23:end procedure

24:procedure
Init($n,d,\\{v_{i}\\}_{i\in[n]}\subseteq\mathbb{R}^{d},\\{\alpha_{i}\\}_{i\in[n]}\subseteq\mathbb{R}$)

25: $v_{i}\leftarrow v_{i}$, $\alpha_{i}\leftarrow\alpha_{i}$ for $i\in[n]$

26: $\mathsf{root}\leftarrow\textsc{BuildTree}(1,n)$

27:end procedure

28:procedure Query($A\in\mathbb{R}^{d\times d}$)

29: ${\sf p}\leftarrow{\sf root},~{}l\leftarrow 1,~{}r\leftarrow n$

30: $s\leftarrow 0$

31: while $l\neq r$ do $\triangleright$ There are $O(\log n)$ iterations

32: $w\leftarrow\langle{\sf p}.V,A\rangle$$\triangleright$ It takes
$O(d^{2})$-time

33: $w_{\ell}\leftarrow\langle({\sf p}.\mathrm{left}).V,A\rangle$

34: Sample $c$ from $\mathrm{Bernoulli}(w_{\ell}/w)$

35: if $c=0$ then

36: ${\sf p}\leftarrow{\sf p}.\mathrm{left}$,
$r\leftarrow\lfloor(l+r)/2\rfloor$

37: else

38: ${\sf p}\leftarrow{\sf p}.\mathrm{right}$,
$l\leftarrow\lfloor(l+r)/2\rfloor+1$

39: end if

40: end while

41: return $l$

42:end procedure

43:end data structure

######  Lemma 7.6 (Running time of Procedure RandBSS+ in Algorithm 3).

Algorithm 3 runs in

  * •

$O(|D|d^{2})$-time for preprocessing,

  * •

$O(d^{2}\log(|D|)+d^{\omega})$-time per iteration, and

  * •

$O(\varepsilon^{-1}d)$ iterations.

Thus, the total running time is,

$\displaystyle O(|D|d^{2}+\varepsilon^{-1}d\cdot(d^{2}\log|D|+d^{\omega})).$

###### Proof.

In each call of the Procedure RandBSS+ in Algorithm 3,

  * •

Finding orthonormal basis takes $O(|D|d^{2})$.

  * •

In the line 9, it runs $O(|D|d^{2})$ times.

  * •

The while loop repeat $O(\varepsilon^{-1}d)$ times.

    * –

Line 14 is computing $(u_{j}I-B_{j})\in\mathbb{C}^{d\times d}$,
$(u_{j}I-B_{j})^{-1}$. This part takes $O(d^{\omega})$ time141414Note that
this step seems to be very difficult to speed up via the Sherman-Morrison
formula since $u_{j}$ changes in each iteration and the update is of high
rank..

    * –

Note that line 15 of Procedure RandBSS+ in Algorithm 3 runs $O(d^{2}\log|D|)$
times.

So, the time complexity of Procedure RandBSS+ in Algorithm 3 is

$\displaystyle O(|D|d^{2}+\varepsilon^{-1}d\cdot(d^{2}\log|D|+d^{\omega})).$

∎

######  Lemma 7.7 (Correctness of Procedure RandBSS+ in Algorithm 3).

Given any dimension $d$ linear space ${\cal F}$, any distribution $D$ over the
domain of ${\cal F}$, and any $\varepsilon>0$,
$\textsc{RandBSS+}(d,\mathcal{F},D,\varepsilon)$ is an $\varepsilon$-WBSP that
terminates in $O(d/\varepsilon)$ rounds with probability $1-1/200$.

###### Proof.

We first claim that, for each $j\in[m]$, ${\sf x}_{j}$ has the same
distribution as $D_{j}$, where

$\displaystyle
D_{j}(x)=D(x)\cdot(v(x)^{\top}E_{j}v(x))/\Phi_{j}~{}~{}~{}\forall x\in D$

Notice that sampling from distribution $D_{j}$ can be reformulated as an
Online Quadratic-Form Sampling Problem: the vectors are $\\{v(x)\\}_{x\in D}$
, the coefficients are $\\{D(x)\\}_{x\in D}$, and the query matrix is
$E^{\prime}_{j}:=E_{j}/\Phi_{j}$. Then, we have $D_{j}={\cal
D}_{E^{\prime}_{j}}$ defined in Problem 7.4. Hence, by Theorem 7.5, we can use
the data structure (Algorithm 4) to efficiently sample from $D_{j}$.

Therefore, the sample ${\sf x}_{j}$ in each iteration is generated from the
same distribution as the original randomized BSS algorithm (Algorithm 2).
Then, the WBSP guarantee and the number of iterations immediately follow from
the proof of [CP19a, Lemma 5.1].

The proof of the lemma is then completed. ∎

###### Proof of Theorem 7.3.

The running time of the algorithm follows from Lemma 7.6, and the correctness
follows from Lemma 7.7. ∎

###  7.3 Trade-off between preprocessing and query

In this section, we consider the preprocessing and query trade-off in the data
structure for quadratic form sampling problem. In the following theorem, we
give a new data structure that takes less time in preprocessing and more time
for each query than Theorem 7.5, and the space complexity is also reduced from
$O(nd^{2})$ to $O(nd)$.

######  Theorem 7.8.

There is a data structure (Algorithms 5 and 6) that uses $O(nd)$ spaces for
the Online Quadratic-Form Sampling Problem with the following procedures:

  * •

Init$(n,d,\\{v_{1},\dots,v_{n}\\}\subset\mathbb{R}^{d},\\{\alpha_{1},\dots,\alpha_{n}\\}\subset\mathbb{R})$:
the data structure preprocesses in time $O(nd^{\omega-1})$.

  * •

Query$(A\in\mathbb{R}^{d\times d})$: Given a PSD matrix $A$, the Query
operation samples $i\in[n]$ exactly from the probability distribution ${\cal
D}_{A}$ defined in Problem 7.4 in $O(d^{2}\log(n/d)+d^{\omega})$-time.

###### Proof.

The time and space complexities follow from Lemma 7.9. And the correctness
follows from Lemma 7.10. ∎

Algorithm 5 Quadratic-form sampling with preprocessing-query trade-off:
Preprocessing

1:structure Node

2: $V_{1},V_{2}\in\mathbb{R}^{d\times d}$

3: $\mathrm{left},\mathrm{right}$ $\triangleright$ Point to the left/right
child in the tree

4:end structure

5:data structure DS+$\triangleright$ Theorem 7.8

6:members

7: $n\in\mathbb{N}$$\triangleright$ The number of vectors

8: $m\in\mathbb{N}$$\triangleright$ The number of blocks

9: $v_{1},\dots,v_{n}\in\mathbb{R}^{d}$ $\triangleright$ $d$-dimensional
vectors

10: $\mathsf{root}$: Node $\triangleright$ The root of the tree

11:end members

12:procedure BuildTree($l,r$)$\triangleright$ $[l,\dots,r]$ is the range of
the current node

13: ${\sf p}\leftarrow\mathbf{new}$ Node

14: if $l=r$ then$\triangleright$ Leaf node

15: ${\sf
p}.V_{2}\leftarrow\begin{bmatrix}v_{(l-1)d+1}&\cdots&v_{ld}\end{bmatrix}$

16: ${\sf p}.V_{1}\leftarrow({\sf p}.V_{2})\cdot({\sf p}.V_{2})^{\top}$
$\triangleright$ It takes $O(d^{\omega})$-time

17: $\triangleright$ ${\sf
p}.\mathrm{mat1}=\sum_{i=(l-1)d+1}^{ld}v_{i}v_{i}^{\top}$

18: else$\triangleright$ Internal node

19: $mid\leftarrow\lfloor(l+r)/2\rfloor$

20: ${\sf p}.\mathrm{left}\leftarrow\textsc{BuildTree}(l,mid)$

21: ${\sf p}.\mathrm{right}\leftarrow\textsc{BuildTree}(mid+1,r)$

22: ${\sf p}.V_{1}\leftarrow({\sf p}.\mathrm{left}).V_{1}+({\sf
p}.\mathrm{right}).V_{1}$ $\triangleright$ It takes $O(d^{2})$-time

23: end if

24: return ${\sf p}$

25:end procedure

26:procedure
Init($n,d,\\{v_{i}\\}_{i\in[n]}\subseteq\mathbb{R}^{d},\\{\alpha_{i}\\}_{i\in[n]}\subseteq\mathbb{R}$)

27: $v_{i}\leftarrow v_{i}\cdot\sqrt{\alpha_{i}}$ for $i\in[n]$

28: $m\leftarrow n/d$$\triangleright$ We assume that $n$ is divisible by $d$

29: Group $\\{v_{i}\\}_{i\in[n]}$ into $m$ blocks
$B_{1},\dots,B_{m}$$\triangleright$ $B_{i}=\\{v_{(i-1)d+1},\dots,v_{id}\\}$
for $i\in[m]$

30: $\mathsf{root}\leftarrow\textsc{BuildTree}(1,m)$

31:end procedure

32:end data structure

Algorithm 6 Quadratic-form sampling with preprocessing-query trade-off: Query

1:data structure DS+$\triangleright$ Theorem 7.8

2:members

3: $n\in\mathbb{N}$$\triangleright$ The number of vectors

4: $m\in\mathbb{N}$$\triangleright$ The number of blocks

5: $v_{1},\dots,v_{n}\in\mathbb{R}^{d}$ $\triangleright$ $d$-dimensional
vectors

6: $\mathsf{root}$: Node $\triangleright$ The root of the tree

7:end members

8:procedure BlockSampling(${\sf p}$, $l\in\mathbb{N}$,
$A\in\mathbb{R}^{d\times d}$)$\triangleright$ ${\sf p}$ is a leaf node with
index $l$

9: $U\leftarrow({\sf p}.V_{2})^{\top}\cdot A\cdot({\sf
p}.V_{2})$$\triangleright$ It takes $O(d^{\omega})$-time

10: Define a distribution ${\cal D}_{l}$ over $[d]$ such that ${\cal
D}_{l}(i)\propto U_{i,i}$

11: Sample $i\in[d]$ from ${\cal D}_{l}$$\triangleright$ It takes $O(d)$-time

12: return $(l-1)d+i$

13:end procedure

14:procedure Query($A\in\mathbb{R}^{d\times d}$)

15: ${\sf p}\leftarrow{\sf root},~{}l\leftarrow 1,~{}r\leftarrow m$

16: $s\leftarrow 0$

17: while $l\neq r$ do $\triangleright$ There are $O(\log m)$ iterations

18: $w\leftarrow\langle{\sf p}.V_{1},A\rangle$$\triangleright$ It takes
$O(d^{2})$-time

19: $w_{\ell}\leftarrow\langle({\sf p}.\mathrm{left}).V_{1},A\rangle$

20: Sample $c$ from $\mathrm{Bernoulli}(w_{\ell}/w)$

21: if $c=0$ then

22: ${\sf p}\leftarrow{\sf p}.\mathrm{left}$,
$r\leftarrow\lfloor(l+r)/2\rfloor$

23: else

24: ${\sf p}\leftarrow{\sf p}.\mathrm{right}$,
$l\leftarrow\lfloor(l+r)/2\rfloor+1$

25: end if

26: end while

27: return BlockSampling(${\sf p}$, $l$, $A$)

28:end procedure