Sparse Recovery for Orthogonal Polynomial Transforms

TL;DR

This paper develops sublinear-time algorithms for sparse recovery in orthogonal polynomial transforms, extending efficient sparse Fourier transform techniques to a broad class of classical orthogonal polynomials like Jacobi, Chebyshev, and Legendre.

Contribution

It introduces a general reduction method for sparse recovery in orthogonal polynomial transforms and provides algorithms for Jacobi polynomial-based transforms, broadening the scope beyond Fourier transforms.

Findings

Algorithms run in polynomial time in k and log N

Effective for transforms from Jacobi polynomials including Chebyshev and Legendre

High probability success for vectors with random support

Abstract

In this paper we consider the following sparse recovery problem. We have query access to a vector $\vx \in \R^N$ such that $\vhx = \vF \vx$ is $k$-sparse (or nearly $k$-sparse) for some orthogonal transform $\vF$. The goal is to output an approximation (in an $\ell_2$ sense) to $\vhx$ in sublinear time. This problem has been well-studied in the special case that $\vF$ is the Discrete Fourier Transform (DFT), and a long line of work has resulted in sparse Fast Fourier Transforms that run in time $O(k \cdot \mathrm{polylog} N)$. However, for transforms $\vF$ other than the DFT (or closely related transforms like the Discrete Cosine Transform), the question is much less settled. In this paper we give sublinear-time algorithms---running in time $\poly(k \log(N))$---for solving the sparse recovery problem for orthogonal transforms $\vF$ that arise from orthogonal polynomials. More precisely, our algorithm works for any $\vF$ that is an orthogonal polynomial transform derived from Jacobi polynomials. The Jacobi polynomials are a large class of classical orthogonal polynomials (and include Chebyshev and Legendre polynomials as special cases), and show up extensively in applications like numerical analysis and signal processing. One caveat of our work is that we require an assumption on the sparsity structure of the sparse vector, although we note that vectors with random support have this property with high probability. Our approach is to give a very general reduction from the $k$-sparse sparse recovery problem to the $1$-sparse sparse recovery problem that holds for any flat orthogonal polynomial transform; then we solve this one-sparse recovery problem for transforms derived from Jacobi polynomials.

Tables

Table 1. Table 1 : Summary of notation used in the 1-sparse Jacobi solver.

Fixed or latent parameters
$\alpha ,\beta$	parameters of the fixed Jacobi family we are considering
$𝐅$	Fixed orthogonal polynomial evaluation matrix (defined in Theorem 5.8).
$v,\ell ,{\theta }_{\ell }$	unknowns to be determined
$𝐰,𝐲$	algorithm is given query access to $𝐲=v\cdot {𝐅}^T{𝐞}_{\ell }+𝐰$
$ϵ$	fixed noise level so that ${\Vert 𝐰\Vert }_2\le ϵ\left|v\right|$
Universal constants and parameters for the fixed Jacobi family
$C^{\alpha ,\beta }$	constant in the Jacobi polynomial root distribution (Theorem 5.2)
$c$	Provides condition for Theorem 5.5 to hold (equation (42))
$c^{\prime }$	constant in the Jacobi evaluation approximation Theorem 5.5
$n,\kappa ,\gamma$	parameters in the Jacobi evaluation approximation Theorem 5.5
Constants we will fix
$\nu$	A constant given to Theorem 5.11 so that it yields a final range for ${\theta }_{\ell }$ of size $O\left(1/\nu N\right)$. Independently, we will only query ${𝒥}_j^{(\alpha ,\beta )}$ for $j\ge \nu N$, so that Theorem 5.5 holds for any ${\theta }_{\ell }$ except for a range of size $O\left(1/\nu N\right)$. It is set in Lemma 5.21.
${\delta }_0$	A constant determining the noise ${ϵ}_0$ that Theorem 5.11 sees, via equation (48). It is set in Lemma 5.21 as a function of $ϵ$.
Indices and other variables
$w$	A “blow-up” factor indexing the values $\cos \left(w{\theta }_{\ell }\right)$ that will be queried during Theorem 5.11. The algorithm will query $w\le \nu N$ to narrow the range for ${\theta }_{\ell }$ sufficiently.
$j$	Indexes the Jacobi polynomials ${𝒥}_j^{(\alpha ,\beta )}\left(x\right)$. The algorithm will query $𝐲\left[j\right]$ for $j\in [\nu N,N-\nu N]$ so that Theorem 5.5 applies.
$\Delta$	Randomly chosen value (in $[2\nu N,N-2\nu N]$) to deduce $\cos \left(w\theta \right)$ from $J_{\Delta -w},J_{\Delta },J_{\Delta +w}$ (Corollary 5.9)
$R$	Number of rounds to perform a given subroutine that has a constant probability of failure
${\rho }_{ϵ}$	A constant depending on $ϵ$ (Definition 5.10)
$\mu$	Error probability parameter of subroutines, as well as of the entire 1-sparse recovery algorithm

Equations

$${\int }_{-1}^1{p}_i(X)p_j(X)w(X)dX={\delta }_{i,j},$$

$${\mathbf{P}}_N[i,j]=p_j({\lambda }_i).$$

$${\int }_{-1}^{1}f(X)w(X)dX=\sum _{\ell =0}^{N-1}f({\lambda }_{\ell })\cdot w_{\ell }.$$

$${\mathbf{F}}_N={\mathbf{D}}_{\mathbf{w}}^{\frac{1}{2}}{\mathbf{P}}_N.$$

$${\mathbf{F}}_N^T{\mathbf{F}}_N={\mathbf{P}}_N^T{\mathbf{D}}_{\mathbf{w}}{\mathbf{P}}_N={\mathbf{I}}_N,$$

$${\mathbf{P}}_N^T{\mathbf{D}}_{\mathbf{w}}{\mathbf{P}}_N[i,j]=\sum _{k=0}^{N-1}{p}_i({\lambda }_k)w_k{p}_j({\lambda }_k)={\int }_{-1}^1{p}_i(X)p_j(X)w(X)dX={\delta }_{i,j}.$$

$${\mathbf{F}}_N^{-1}={\mathbf{F}}_N^T.$$

$$w^{(\alpha ,\beta )}(X)=(1-X{)}^{\alpha }\cdot (1+X{)}^{\beta }$$

$$\mathbf{D}[j,j]=\frac{2^{\alpha +\beta +1}}{2j+\alpha +\beta +1}\cdot \frac{\Gamma (j+\alpha +1)\Gamma (j+\beta +1)}{\Gamma (j+1)\Gamma (j+\alpha +\beta +1)}$$

$${\mathbf{F}}_N[i,j]=\{\begin{array}{ll}\sqrt{\frac{1}{N}}& j=0\\ \sqrt{\frac{2}{N}}\cdot \cos \left[\frac{\pi }{N}j\left(i+\frac{1}{2}\right)\right]& j=1,\dots ,N-1.\end{array}$$

$$C_0\gamma d\le \left|\left\{{\theta }_0,\dots ,{\theta }_{d-1}\right\}\cap \left[{\theta }_i-\frac{\gamma }{2},{\theta }_i+\frac{\gamma }{2}\right]\right|\le C_1\gamma d.$$

$$\hat{\mathbf{y}}=\tilde{\mathbf{y}}+\mathbf{w},$$

$${\Vert \mathbf{w}\Vert }_2\le ϵ\left|v\right|,$$

$$\Vert \hat{\mathbf{x}}-\hat{\mathbf{z}}{\Vert }_2\lesssim \delta \Vert \hat{\mathbf{x}}{\Vert }_2,$$

$$U=\underset{i,j\in [N{]}^2}{\max }|{\mathbf{F}}_N^{-1}[i,j]|.$$

$${\mu }_0=\frac{{\mu }^2{C}_0^2{\gamma }^2}{k^2}.$$

$$\Vert \hat{\mathbf{x}}-\hat{\mathbf{z}}{\Vert }_2\le 3\delta \Vert \mathbf{x}{\Vert }_2$$

$$\Vert \hat{\mathbf{v}}-\hat{\mathbf{z}}{\Vert }_2\le 3\delta \Vert \mathbf{x}{\Vert }_2+\Vert \mathbf{w}{\Vert }_2.$$

$$\mathrm{poly}\left(\frac{k\log (1/\mu )}{\gamma \delta C_0}\right)\cdot \left(N{U}^2+Q\left(N,\frac{6\delta }{C},\delta ,\frac{{\mu }_0}{2}\right)\right)$$

$$\mathrm{poly}\left(\frac{k\log (1/\mu )}{\gamma \delta C_0}\right)\cdot \left(N{U}^2+T\left(N,\frac{6\delta }{C},\delta ,\frac{{\mu }_0}{2}\right)\right)$$

$$\left|\hat{\mathbf{z}}[i]-\hat{\mathbf{v}}[i]\right|\le Cϵ\left|\hat{\mathbf{v}}[i]\right|\text{ and }\left|\hat{\mathbf{v}}[i]\right|\ge (1-2Cϵ)\cdot \Vert \hat{\mathbf{v}}-\hat{\mathbf{z}}{\Vert }_{\infty }.$$

$$\left|\tilde{v}-\hat{\mathbf{v}}[h]\right|\le Cϵ\left|\hat{\mathbf{v}}[h]\right|\text{ and }\left|\hat{\mathbf{v}}[h]\right|\ge (1-2Cϵ)\Vert \hat{\mathbf{v}}-\tilde{\mathbf{z}}{\Vert }_{\infty }.$$

$$O\left(T_0\cdot T_2\cdot \frac{\sqrt{k}}{ϵ\gamma }\left(T_1+Q(N,6ϵ,Cϵ,{\mu }_0/2)\right)\right)$$

$$O\left(T_0\cdot T_2\cdot \frac{k^{3/2}}{ϵ\gamma }\left(T_1+T(N,6ϵ,Cϵ,{\mu }_0/2)\right)\right)$$

$$(1-ϵ)\left|\hat{\mathbf{x}}[i]\right|\le \left|\hat{\mathbf{v}}[i]\right|\le (1+ϵ)\left|\hat{\mathbf{x}}[i]\right|$$

$$(1-ϵ)\left|\hat{\mathbf{v}}[i]\right|\le \left|\hat{\mathbf{x}}[i]\right|\le (1+ϵ)\left|\hat{\mathbf{v}}[i]\right|.$$

$$|v-\hat{\mathbf{x}}[i]|\le 2ϵ(C+1)|\left|\hat{\mathbf{x}}[i]\right|.$$

$$\left|\hat{\mathbf{x}}[i]-\hat{\mathbf{v}}[i]\right|$$

$$\le \Vert \hat{\mathbf{w}}{\Vert }_2$$

$$\le (ϵ/2)\text{textsc}large(k,\hat{\mathbf{x}})$$

Full text

In this paper we give sublinear-time algorithms—running in time $\mathrm{poly}(k\log(N))$—for solving the sparse recovery problem for orthogonal transforms $\mathbf{F}$ that arise from *orthogonal polynomials. *More precisely, our algorithm works for any $\mathbf{F}$ that is an orthogonal polynomial transform derived from Jacobi polynomials. The Jacobi polynomials are a large class of classical orthogonal polynomials (and include Chebyshev and Legendre polynomials as special cases), and show up extensively in applications like numerical analysis and signal processing. One caveat of our work is that we require an assumption on the sparsity structure of the sparse vector, although we note that vectors with random support have this property with high probability.

Our approach is to give a very general reduction from the $k$-sparse sparse recovery problem to the $1$-sparse sparse recovery problem that holds for any flat orthogonal polynomial transform; then we solve this one-sparse recovery problem for transforms derived from Jacobi polynomials. Frequently, sparse FFT algorithms are described as implementing such a reduction; however, the technical details of such works are quite specific to the Fourier transform and moreover the actual implementations of these algorithms do not use the $1$-sparse algorithm as a black box. In this work we give a reduction that works for a broad class of orthogonal polynomial families, and which uses any $1$-sparse recovery algorithm as a black box.

1 Introduction

In this paper, we consider the following sparse recovery problem. Suppose that we have query access to a vector $\mathbf{x}\in\mathbb{R}^{N}$, which has the property that for a fixed orthogonal transform matrix $\mathbf{F}$, $\hat{\mathbf{x}}=\mathbf{F}\mathbf{x}$ is $k$-sparse (or approximately $k$-sparse, in the sense that $\hat{\mathbf{x}}$ is close in $\ell_{2}$ distance to a $k$-sparse vector). The goal is to recover an approximation $\hat{\mathbf{z}}$ to $\hat{\mathbf{x}}$, so that $\|\hat{\mathbf{x}}-\hat{\mathbf{z}}\|_{2}$ is small with high probability, as quickly as possible.

Variants of this problem have been studied extensively over several decades—we refer the reader to the book [Foucart2013] for many examples and references. One particularly well-studied example is the *sparse Fast Fourier Transform *(sFFT)—see the survey [FourierSurvey2014] and the references therein. In this case, the matrix $\mathbf{F}$ is taken to be the Discrete Fourier Transform (DFT) and a long line of work has produced near-optimal results: algorithms with running time $O(k\operatorname{polylog}(N))$ and sample complexity $O(k\log N)$ [HIKP12, IKP14, IK14].

We study the sparse recovery problem for a more general class of transforms $\mathbf{F}$ called *orthogonal polynomial transforms, *and in particular those that arise from *Jacobi polynomials, *a broad class of orthogonal polynomials (OPs). Jacobi polynomials include as special cases many familiar families of OPs, including Gegenbauer and in particular Chebyshev, Legendre, and Zernike111To be more precise the radial component of a Zernike polynomial is a Gegenbauer and hence, a Jacobi polynomial. polynomials, and the corresponding OP transforms appear throughout numerical analysis and signal processing.

Despite the progress on the sFFT described above, much remains unknown for general orthogonal polynomial transforms. As discussed more in Section 1.2 below, the *sample complexity *of the sparse recovery problem is well understood, and the ‘correct’ answer is known to be $\Theta(k\operatorname{polylog}(N))$ queries to $\mathbf{x}$. However, the algorithmic results that go along with these sample complexity bounds result in $\mathrm{poly}(N)$ time algorithms. Our goal in this work will be *sublinear *time algorithms as well as sublinear sample complexity. There are sublinear-time algorithms available for the special cases of Chebyshev and Legendre polynomials that work by essentially reducing to the Fourier case. For general Jacobi polynomials, such reductions are not available. We elaborate in Appendix B why reducing general Jacobi polynomials to the Fourier case does not seem easy. There are also algorithms based on Prony’s method, some of which work for quite general families of OPs [prony-general]. However these general results require exact sparsity; to the best of our knowledge versions of Prony’s method that are provably robust to noise are restricted to classes of OPs similar to the Fourier transform.

Results.

In this work, we give the first (to the best of our knowledge) sublinear-time algorithms with provable guarantees for the (approximately-)sparse recovery problem for general orthogonal transforms derived from Jacobi polynomials. We discuss our results in more detail in Section 3 and briefly summarize them here. Our algorithms run in time $\mathrm{poly}(k\log(N))$ and given query access to $\mathbf{v}=\mathbf{F}^{-1}\hat{\mathbf{v}}$, can find approximations to $\hat{\mathbf{v}}$ when $\hat{\mathbf{v}}$ is approximately $k$-sparse of an appropriate form. More precisely, we can handle vectors $\hat{\mathbf{v}}=\hat{\mathbf{x}}+\hat{\mathbf{w}}$ where $\hat{\mathbf{x}}$ is $k$-sparse with a ‘spread-out’ support (made precise in Definition 2.3), and $\hat{\mathbf{w}}$ is an adversarial noise vector with sufficiently small $\ell_{2}$ norm. We obtain guarantees of the following flavor: for any such vector $\mathbf{v}$, we can find $\hat{\mathbf{z}}$ such that $\|\hat{\mathbf{z}}-\hat{\mathbf{x}}\|_{2}\leq 0.01\|\hat{\mathbf{x}}\|_{2}$ with high probability.

We note that these results are weaker than the results for the sFFT: our sample complexity and running time are polynomially larger, and we need stronger assumptions on the sparse signals. However, we also note that the decade or so of work on the sFFT culminating in the results above began with similar results (see [GGIMS], for example, in which the dependence on $k$ is an unspecified polynomial) and we hope that this work will similarly be a first step towards near-optimal algorithms for general orthogonal polynomial transforms.

Techniques.

Our techniques follow the outline of existing algorithms for the sFFT, although as we elaborate on in Section 1.3, the situation for general Jacobi polynomials is substantially more complicated. More precisely, we first give a very general reduction, which reduces the $k$-sparse case to the $1$-sparse case. The idea of such a reduction was implicit in the sFFT literature, but previous work has relied heavily on the structure of the DFT. Our reduction applies to a broad class of OPs including Jacobi polynomials. Next, we show how to solve the $1$-sparse recovery problem for general Jacobi polynomials. The basic idea is to use known approximations of Jacobi polynomial evaluations by certain cosine evaluations [szego] in order to iteratively narrow down the support of the unknown $1$-sparse vector. We give a more detailed overview of our techniques in Section 1.3.

Organization.

For the rest of the introduction, we briefly introduce orthogonal polynomial transforms, discuss previous work, and give a high-level overview of our approach. After that we introduce the formal notation and definitions we need in Section 2, after which we state our results more formally in Section 3. Then we prove our main results: the reduction from $k$ to $1$-sparse recovery is proved in Section 4, the $1$-sparse recovery algorithm for Jacobi polynomials is presented in Section 5, and the resulting $k$-sparse recovery algorithm for Jacobi polynomials is presented in Section 6.

1.1 Orthogonal Polynomial Transforms

Orthogonal polynomials (OPs) play an important role in classical applied mathematics, mathematical physics, and the numerical analysis necessary to simulate solutions to such problems. We give more precise definitions in Section 2 but briefly a family of orthogonal polynomials $p_{0}(X),p_{1}(X),\ldots$ is a collection of polynomials defined on an interval ${\cal D}$ of $\mathbb{R}$, that are pairwise orthogonal with respect to a (non-negative) weight function $w$.

In this work we study *Jacobi polynomials *(defined formally in Section 2), which are a very general class of orthogonal polynomials. These include Chebyshev polynomials, Legendre polynomials, Zernike polynomials and more generally Gegenbauer polynomials. These OP families show up in many places. For example, Zernike polynomials are a family of orthogonal polynomials on the unit disk that permit an analytic expression of the 2D Fourier transform on the disk. They are used in optics and interferometry [Tango1977]. They can be utilized to extract features from images that describe the shape characteristics of an object and were recently used for improved cancer imaging [yu2016]. Different families of orthogonal polynomials give rise to different quadrature rules for numerical integration [R81, quadrature]. Specifically, Chebyshev polynomials are used for numerical stability (see e.g. the ChebFun package [chebfun]) as well as approximation theory (see e.g. Chebyshev approximation [cheb-approx]). Chebyshev polynomials also have certain optimal extremal properties, which has resulted in many uses in theoretical computer science, including in learning theory, quantum complexity theory, linear systems solvers, eigenvector computation, and more [musco]. Further, Jacobi polynomials form solutions of certain differential equations [jacobi]. More recent applications include Dao et al.’s [DaoDeSaRe2017] use of orthogonal polynomials to derive quadrature rules for computing kernel features in machine learning.

Orthogonal polynomials naturally give rise to (discrete) orthogonal polynomial transforms. Suppose that $\mathbf{F}$ is an $N\times N$ matrix, with each column corresponding to an orthogonal polynomial $p_{0},\ldots,p_{N-1}$ and each row an evaluation point $\lambda_{0},\ldots,\lambda_{N-1}$ in a suitable domain and suitably normalized so that it is an orthogonal matrix (Definition 2.1). A familiar example might be the DFT: in this language, the DFT matrix is defined by the polynomials $1,X,X^{2},\ldots,X^{N-1}$, evaluated at points $\lambda_{j}=\omega^{j}$ where $\omega$ is the $N$th root of unity.222We note that in this work we consider a setting slightly different than this example, where $\mathcal{D}=[-1,1]$ rather than $\mathbf{S}^{1}$. Like the Fourier Transform, it is known that all OP transforms admit ‘fast’ versions, allowing matrix-vector multiplication in time $O(N\log^{2}(N))$ [driscoll].333We note that even though the work of [driscoll] has in some sense solved the problem of computing any OP transform in near-linear time, many practical issues still remain to be resolved and the problem of computing OP transforms in near-linear time has seen a lot of research activity recently. We just mention two recent works [jacobi-1, jacobi-2] that present near-linear time algorithms for the Jacobi polynomial transforms (and indeed their notion of uniform Jacobi transform corresponds exactly to the Jacobi polynomial transform that we study in this paper). However, these algorithms inherently seem to require at least linear-time and it is not clear how to convert them into sub-linear algorithms, which is the focus of our work. Thus, our problem of sparse recovery for OP transforms is a natural extension of the sFFT problem, with applications to the areas mentioned above.

1.2 Related Work

As previously described, there has been a great deal of work on the sFFT; we refer the reader to the survey [FourierSurvey2014] for an overview. There has also been work on non-Fourier OP transforms. We break up our discussion below into discussion on the *sample complexity *(which as mentioned above is largely settled) and the *algorithmic complexity *(which remains largely open).

Sample complexity.

The sample complexity of OP transforms $\mathbf{F}$ has been largely pinned down by the *compressed sensing *literature. For example, suppose that $\mathbf{F}\in\mathbb{R}^{N\times N}$ is any orthogonal and sufficiently flat matrix, in the sense that none of the entries of $\mathbf{F}$ are too large. Then a result of Rudelson and Vershynin (and a sharpening of their result by Bourgain) shows that $m=O(k\log{k}\log^{2}{N})$ samples suffice to establish that the matrix $\Phi\in\mathbb{R}^{m\times N}$ (which is made up of $m$ sampled rows from $\mathbf{F}^{T}$) has the *Restricted Isometry Property *(RIP) [rip-B14, RV]. Finding $\hat{\mathbf{x}}=\mathbf{F}\mathbf{x}$ from samples of $\mathbf{F}$ of corresponds to the problem of finding an (approximately) $k$-sparse vector $\hat{\mathbf{x}}$ from the linear measurements $\Phi\hat{\mathbf{x}}$, which is precisely the compressed sensing problem. It is known that if $\Phi$ satisfies the RIP, then this can be solved (for example with $\ell_{1}$ minimization) in time $N^{O(1)}$. We note that very recently a result due to Błasiok et al. show that this is essentially tight, in that $O(k\log^{2}{N})$ queries (for a certain range of $k$) to $\mathbf{x}$ are not enough to compute a $k$-sparse approximation of $\mathbf{F}\mathbf{x}$ [Blasiok2019]. Bounds specific to the DFT over finite fields can be found in [Rao2019].

Foucart and Rauhut [Foucart2013] show that if the orthogonal polynomials satisfy a Bounded Orthogonal System (BOS) which are suitably flat, then if the $m$ evaluation points $\lambda_{j}$ are chosen uniformly at random proportional to the weight function $w$, then the $m\times N$ matrix $\Phi$ defined by normalizing $\mathbf{P}_{N}[i,j]=p_{j}(\lambda_{i})$ appropriately satisfies the RIP with high probability provided that $m$ has an appropriate dependence on $N,k,\epsilon,$ and the flatness of the matrix, and this again gives an $N^{O(1)}$-time algorithm to solve the sparse recovery problem.

Rauhut and Ward [RW12] show that for Jacobi polynomial transforms if the evaluation points were picked according to the Chebyshev measure, then with $O(k\,\mathrm{polylog}{N})$ random measurements, the corresponding matrix has the RIP (note that the Foucart and Rauhut sample the evaluation points according to the measure of orthogonality for the Jacobi polynomials, which in general is not the Chebyshev measure). This result again does not give a sub-linear time algorithm but was used in the result of [legendre-sparse] which we describe below.

While these approaches can give near-optimal sample complexity, they do not give sublinear-time algorithms. In fact, it is faster to compute $\hat{\mathbf{x}}$ exactly by computing $\mathbf{F}\mathbf{x}$, if we care only about the running time and not about sample complexity [driscoll]. Thus, we turn our attention to sublinear-time algorithms.

Sublinear-time algorithms for OP transforms.

There have been several works generalizing and building on the sFFT results mentioned above. One direction is to the multi-dimensional DFT (for example in [IK14, KVZ19]). Another direction is to apply the sFFT framework to orthogonal polynomials with similar structure. One example is Chebyshev polynomials and the Discrete Cosine Transform (DCT). It was observed in [legendre-sparse] (also see Appendix A) that this can be reduced to sFFT in a black box manner, solving the sparse recovery problem for Chebyshev polynomials and the DCT. A second example of OP transforms which can essentially be reduced to the sFFT is *Legendre polynomials. *Hu et al. [legendre-sparse] seek to recover an unknown $k$-term Legendre polynomial (with highest magnitude degree limited to be $N/2$), defined on $[-,1,1]$, from samples. They give a sublinear two-phase algorithm: in the first phase, they reduce $k$-sparse-Legendre to sFFT to identify a set of candidate Legendre polynomials. The second phase uses the RIP result for BOS to produce a matrix that is used to estimate the coefficients of the candidate Legendre polynomials. We note that in this work the setting is naturally continuous, while ours is discrete.

Choi et al. [Choi2018] study higher dimensions and obtain sublinear-time algorithms for more general harmonic expansions in multiple dimensions. The results of [Choi2018] complement our work. More precisely, that work shows how to use any algorithm for a univariate polynomial transform to design an algorithm for a multi-variate polynomial transform where the multi-variate polynomials are products of univariate polynomials in the individual variables. Thus our improvements for univariate polynomial transforms can be used with [Choi2018].

Finally, there are sparse OP transforms based on Prony’s method. The work [prony-general] extends Prony’s method to a very general setting, including Jacobi polynomials, and gives an algorithm that requires only $O(k)$ queries to recover exactly $k$-sparse polynomials. However, these general results work only for exact sparsity and are in general not robust to noise. There has been work extending and modifying these techniques to settings with noise (for example, [HuaSarkar1990, Potts2010]), but to the best of our knowledge the only provable results for noise are for either the sFFT or closely related polynomial families. We note that [Potts2016] presents a Prony-like algorithm for Legendre and Gegenbauer polynomials and demonstrates empirically that this algorithm is robust to noise, although they do not address the question theoretically.

1.3 Technical overview

Our technical results have two main parts. First, inspired by existing approaches to the sFFT, we provide a general reduction from the $k$-sparse recovery problem to the $1$-sparse recovery algorithm, which works for any family of OPs that is sufficiently ‘flat:’ that is, no entry of the matrix $\mathbf{F}$ is too large. Second, we provide a $1$-sparse recovery algorithm for Jacobi polynomials. We give an overview of both parts below.

For what follows, let $\mathbf{F}$ be an orthogonal matrix. For simplicity in this overview we will assume that there is no noise. That is, we want to compute the exactly $k$-sparse $\hat{\mathbf{x}}=\mathbf{F}\mathbf{x}$ given query access to $\mathbf{x}$. However, we note that our final results do work for approximately $k$-sparse vectors $\hat{\mathbf{v}}=\hat{\mathbf{x}}+\hat{\mathbf{w}}$ provided that $\|\hat{\mathbf{w}}\|_{2}$ is sufficiently small.

1.3.1 Reduction to one-sparse recovery

We give a general reduction from the $k$-sparse recovery problem to the one-sparse recovery problem, which works for a broad class of OP families defined on a finite interval.444We note that our results do not (yet) work for the case when the orthogonality is defined over an infinite interval. In particular, our reduction does not work for the Hermite and Laguerre polynomials. At a high level, the idea is as follows. Suppose that $\hat{\mathbf{x}}=\mathbf{F}\mathbf{x}$ is $k$-sparse and $\mathbf{b}\in\mathbb{R}^{N}$ is a ‘filter’: at this stage it is helpful to think of it like a boxcar filter, so $\mathbf{b}$ is $1$ on some interval $I$ and zero outside of that interval. If we choose this interval randomly, we might hope to isolate a single ‘spike’ of $\hat{\mathbf{x}}$ with $\mathbf{b}$: that is, we might hope that $\mathbf{D}_{\mathbf{b}}\hat{\mathbf{x}}$ is one-sparse, where $\mathbf{D}_{\mathbf{b}}$ is the diagonal matrix with $\mathbf{b}$ on the diagonal. Suppose that this occurs. In order to take advantage of this one-sparse vector with a black-box solution to the one-sparse recovery problem, we would need query access to the vector $\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\hat{\mathbf{x}}=\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\mathbf{F}\mathbf{x}$, while what we have is query access to $\mathbf{x}$. Thus, we would like to design $\mathbf{b}$ so that $\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\mathbf{F}$ is *row-sparse. *This would allow us to query a position of $\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\hat{\mathbf{x}}$ using only a few queries from $\mathbf{x}$.

One of our main technical contributions is showing how to design such a vector $\mathbf{b}$, so that $\mathbf{b}$ approximates a boxcar filter and so that $\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\mathbf{F}$ is row-sparse for *any *OP transform $\mathbf{F}$.

Then, given this filter, we can iteratively identify and subtract off ‘spikes’ in $\hat{\mathbf{x}}$ until we have recovered the whole thing. Of course, the actual details are much more complicated than the sketch above. First, the one-sparse solver might have a bit of error, which will get propagated through the algorithm. Second, in our analysis the vector $\hat{\mathbf{x}}$ need not be exactly $k$-sparse. Third, $\mathbf{b}$ will only approximate a boxcar filter, and this is an additional source of error that needs to be dealt with. We will see how to overcome these challenges in Section 4.

For the reader familiar with the sFFT, this approach might look familiar: most sFFT algorithms work by using some sort of filter to isolate single spikes in an approximately sparse signal. Below, we highlight some of the challenges in extending this idea beyond the Fourier transform. Some of these challenges we have overcome, and one we have not (yet) overcome. We mention this last open challenge both because it explains the assumption we have to make on the sparsity structure of $\hat{\mathbf{x}}$, and also because we hope it will inspire future work.

Challenge 1: Choice of filter.

One key difficulty in extending sFFT algorithms to general orthogonal polynomials is that the filters used in the sFFT approach are very specific to the Fourier transform. Indeed, much of the progress that has been made on that problem has been due to identifying better and better choices of filter specialized to the Fourier transform. In order to find filters that work for *any *OP family, we take a different approach and construct a filter out of low-degree Chebyshev polynomials. Then we use the orthogonality properties of the OP family to guarantee that $\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\mathbf{F}$ has the desired sparsity properties.

Challenge 2: Explicit black-box reduction.

Because our goal is generality (to as broad a class of OPs as possible), we give an explicit reduction that uses a $1$-sparse solution as a black box. To the best of our knowledge, existing work on the sFFT does not explicitly do this: a reduction of this flavor is certainly implicit in many of these works, and even explicitly given as intuition, but we are not aware of an sFFT algorithm which actually uses a $1$-sparse recovery algorithm as a black box.

Challenge 3: Equi-spaced evaluation points.

The evaluations points in DFT and the DCT are equispaced (in the angular space). This fact is crucially exploited in sFFT algorithms (as well as the reduction of DCT to DFT—see Appendix A for more details on the reduction). Unfortunately, the roots of Jacobi polynomials are no longer equally spaced—see Appendix B for why this is a barrier to reducing $k$-sparse recovery of Jacobi polynomials directly to $k$-sparse recovery for DCT or DFT. However, it is known that the roots of Jacobi polynomials are ‘spread out,’ (in a sense made below precise in Definition 2.2), and we show that this property is enough for our reduction. In fact, our reduction from $k$-sparse recovery to $1$-sparse recovery works generally for any ‘flat’ OP family with ‘spread out’ roots.

(Open) Challenge 4: Permuting the coordinates of $\hat{\mathbf{x}}$.

In the approach described above, we hoped that an interval $I$ would ‘isolate’ a single spike. In the sFFT setting, this can be achieved through a permutation of the coordinates of $\hat{\mathbf{x}}$. In our language, in the sFFT setting it is possible to define a random (enough) permutation matrix $\mathbf{P}$ so that $\mathbf{P}\hat{\mathbf{x}}$ has permuted coordinates, and so that $\mathbf{F}^{-1}\mathbf{D}_{\mathbf{b}}\mathbf{P}\mathbf{F}$ is row-sparse—this argument crucially exploits the fact that the roots of unity are equispaced in the angle space. This means that not only can we sample from the one-sparse vector $\mathbf{D}_{\mathbf{b}}\hat{\mathbf{x}}$, but also we can sample from $\mathbf{D}_{\mathbf{b}}\mathbf{P}\hat{\mathbf{x}}$, and then there is some decent probability that any given spike in $\hat{\mathbf{x}}$ is isolated by $\mathbf{b}$. However, we have not been able to come up with (an approximation to) such a $\mathbf{P}$ that works in the general OP setting. This explains why we require the assumption that the support of $\hat{\mathbf{x}}$ be reasonably ‘spread out,’ so that we can hope to isolate the spikes by $\mathbf{b}$. This assumption is made precise in Definition 2.3. We note that if such a $\mathbf{P}$ were found in future work, this would immediately lead to an improved $k$-sparse recovery result for Jacobi polynomials, which would work for arbitrary sparse signals $\hat{\mathbf{x}}$.

1.3.2 A one-sparse recovery algorithm for Jacobi polynomials

With the reduction complete, to obtain a $k$-sparse recovery algorithm for general Jacobi polynomials we need to solve the one-sparse case. We give an overview of the basic idea here. First, we note that via well-known approximations of Jacobi polynomials [szego], one can approximate the evaluation of any Jacobi polynomial at a point in $(-1,1)$ by evaluating the cosine function at an appropriate angle. Using some standard local error-correcting techniques (for example, computing $\cos(A)$ via $\cos{A}=\frac{\cos(A+B)+\cos(A-B)}{2\cos{B}}$ for a random $B$), we reduce the $1$-sparse recovery problem to computing $\theta$ from noisy values of $\cos(w\theta)$ for some integers $w\geq 1$. Since the reduction is approximate, some care has to be taken to handle some corner cases where the approximation does not hold. In particular, we have to figure out for which real numbers $y\in[0,N)$, does its orbit $\langle xy\rangle$ for $x\in\mathbb{Z}_{N}$ have small order. We give a result to handle this, which to the best of our knowledge (and somewhat surprisingly) seems to be new.555We thank Stefan Steinerberger for showing us a much simpler proof than our original more complicated proof, which also gave worse parameters. With this out of the way, our algorithm to compute the value of $\theta$ from the evaluations $\cos(w\theta)$ is based on the following idea. Assuming we already know $\cos(\theta)$ up to $\pm\epsilon$, we get a noise estimate of $\theta$ (which lives in the range $\arccos(\cos(\theta)\pm\epsilon)$) and then use the evaluations at $w>1$ to ‘dilate’ the region of $[0,\pi]$ where we know $\theta$ lies, reducing $\epsilon$. We proceed iteratively until the region of uncertainty is small enough that there are only $O(1)$ possibilities remaining, which we then prune out using the fact that $\mathbf{F}$ is orthogonal and flat, in the sense that none of its entries are too large. (We note that proving $\mathbf{F}$ is flat needs a bit of care. In particular, we need a sharper bound on Jacobi polynomials (than the cosine approximation mentioned above) in terms of Bessel functions to prove that all entries of $\mathbf{F}$ are small.) Similar ideas have been used for $1$-sparse recovery for the DFT (for example, in [HIKP12]), although our situation is more complicated than the DFT because working with cosines instead of complex exponentials means that we lose sign information about $\theta$ along the way.

2 Background and Preliminaries

2.1 Notation

We use bold lower-case letters ($\mathbf{x},\mathbf{y}$) for vectors and bold upper-case letters ($\mathbf{P},\mathbf{F}$) for matrices. Non-bold notation $x,y,U$ is used for scalars in $\mathbb{R}$. In general, if there is a given transform $\mathbf{F}$ we are considering, then the notation $\hat{\mathbf{x}}\in\mathbb{R}^{N}$ indicates $\mathbf{F}\cdot\mathbf{x}$. We use the notation $\mathbf{x}[i]$ or $\mathbf{X}[i,j]$ to index into a vector or matrix, respectively. All of our vectors and matrices are [math]-indexed, i.e. the entries of a vector $\mathbf{x}\in\mathbb{R}^{N}$ are $\mathbf{x}[0],\dots,\mathbf{x}[N-1]$. We use $[N]$ to denote the set $\{0,\ldots,N-1\}$. Given a subset $S\subset[N]$, we will denote the complement set (i.e. $[N]\setminus S$) by $S^{c}$.

Given a vector $\mathbf{x}\in\mathbb{R}^{N}$ and an integer $1\leq s\leq N$, we define $\textsc{large}(s,\mathbf{x})$ to be the magnitude of the $s$th largest value in $\mathbf{x}$ (by absolute value).

For any vector $\mathbf{u}\in\mathbb{R}^{N}$, we define $\mathbf{D}_{\mathbf{u}}\in\mathbb{R}^{N\times N}$ as the diagonal matrix with $\mathbf{u}$ on its diagonal. For a diagonal matrix $\mathbf{D}$, and any real $\alpha$ we denote $\mathbf{D}^{\alpha}$ to denote the diagonal matrix with the $(i,i)$ entry being $\left({\mathbf{D}[i,i]}\right)^{\alpha}$. Given a vector $\mathbf{x}\in\mathbb{R}^{N}$ and set $S\subseteq[N]$, $\mathbf{x}_{S}$ denotes the vector $\mathbf{x}$ where all entries out of $S$ are masked to [math]. For $\mathbf{x}\in\mathbb{R}^{N}$, $\mathrm{supp}(\mathbf{x})\subseteq[N]$ denotes the support (i.e. the set of non-zero positions) of $\mathbf{x}$.

We use $x\pm h$ to refer to either the interval $[x-h,x+h]$ or a point in this interval, whichever is clear from context. Similarly, if $S$ is an interval $[a,b]$ then $S\pm h$ is the interval $[a-h,b+h]$.

When stating algorithms, we use superscript notation to denote query access. That is $\mathcal{A}^{(\mathbf{x})}(\mathbf{z})$ takes input $\mathbf{z}$ and has query access to $\mathbf{x}$.

We use the notation $f(n)\lesssim g(n)$ to mean that there is some constant $C$ so that, for sufficiently large $n\geq n_{0}$, $f(n)\leq Cg(n)$.

2.2 Orthogonal Polynomials

For the remainder of this paper, we consider polynomials $p_{0}(X),p_{1}(X),\dots$ that form a normalized orthogonal polynomial family with respect to any compactly supported measure $w(X)$. By suitably scaling and translating $X$, we can ensure that the orthogonality is on $[-1,1]$.666See footnote 4. In particular $\deg(p_{i})=i$ and for any $i,j\geq 0$,

[TABLE]

where $\delta_{i,j}=1$ if $i=j$ and [math] otherwise.

Then for given $N$ evaluation points $\lambda_{0},\dots\lambda_{N-1}$, define the orthogonal polynomial transform $\mathbf{P}_{N}$ as follows. For any $0\leq i,j<N$, we have

[TABLE]

In other words, the rows of $\mathbf{P}_{N}$ are indexed by the evaluation points and the columns are indexed by the polynomials.

For the rest of the paper, assume $\lambda_{0}\leq\lambda_{1}\leq\cdots\leq\lambda_{N-1}$ are the roots of $p_{N}(X)$. Then it is well-known (see e.g. [szego]) that

• •

The roots lie in the support of the measure (i.e. $\lambda_{i}\in[-1,1]$) and are distinct (i.e. $\lambda_{0}<\lambda_{1}<\cdots<\lambda_{N-1}$).

• •

There exists weights Gaussian quadrature weights $w_{\ell}=\frac{1}{\sum_{j=0}^{N-1}p_{j}\left({\lambda_{\ell}}\right)^{2}},i=0,\dots,N-1$ such that for any polynomial $f(X)$ of degree at most $2N-1$,

[TABLE]

We are now ready to define the orthogonal matrix corresponding to $\mathbf{P}_{N}$ that we deal with in this paper:

Definition 2.1**.**

Let $p_{0}(X),\dots,p_{N-1}(X),\dots$ be an orthogonal polynomial family, $\lambda_{0},\dots,\lambda_{N-1}$ be the roots of $p_{N}(X)$, and $w_{0},\dots,w_{N-1}$ be the Gaussian quadrature weights. Define $\mathbf{D}_{\mathbf{w}}$ to be the diagonal matrix with $w_{0},\dots,w_{N-1}$ on its diagonal, and

[TABLE]

Note that by (1) and (2),

[TABLE]

so $\mathbf{F}_{N}$ is an orthogonal matrix. In particular,

[TABLE]

Note that since $\mathbf{F}_{N}$ is orthogonal, by definition we have

[TABLE]

2.2.1 Jacobi Polynomials and Special Cases

In this section we define Jacobi Polynomials, our main object of interest, and point out a few special cases. We note that families of named orthogonal polynomials $\{p_{i}(X)\}$ are sometimes defined through different means, hence are normalized differently up to constants. The corresponding discrete orthogonal polynomial transform (e.g. Discrete Legendre Transform) frequently refers to multiplication by $\mathbf{P}$ instead of $\mathbf{F}$. In these cases, the transform satisfies $\mathbf{P}_{N}^{T}\mathbf{D}_{\mathbf{w}}\mathbf{P}_{N}=\mathbf{D}$ for a diagonal matrix $\mathbf{D}$ corresponding to the normalization. The transform $\mathbf{F}_{N}=\mathbf{D}_{\mathbf{w}}^{\frac{1}{2}}\mathbf{P}\mathbf{D}^{-\frac{1}{2}}$ we consider (note that this matrix is indeed orthogonal) is thus equivalent up to diagonal multiplication.

Jacobi polynomials.

Jacobi polynomials are indexed by two parameters $\alpha,\beta>-1$ and these are polynomials $\left\{P_{j}^{(\alpha,\beta)}\right\}_{j\geq 0}$ that are orthogonal with respect to the measure

[TABLE]

in the range $[-1,1]$. This definition is not normalized, in the sense that we have $\mathbf{P}_{N}^{T}\mathbf{D}_{\mathbf{w}}\mathbf{P}_{N}=\mathbf{D}$, where

[TABLE]

(see [szego, Pg. 68, (4.3.3)]). We will come back to the normalization in Section 5 (cf. Corollary 5.8).

We record three well-known special cases: Chebyshev polynomials (of the first kind) are special case of $\alpha=\beta=-\frac{1}{2}$ and Legendre polynomials are the special case of $\alpha=\beta=0$ (up to potentially a multiplicative factor that could depend on the degree $j$). Another notable special case of Jacobi polynomials are the Gegenbauer or ultraspherical polynomials ($\alpha=\beta$).

Chebyshev polynomials of the 1st kind.

The Chebyshev polynomials of the 1st kind are orthogonal with respect to the weight measure $w(X)=(1-X^{2})^{-\frac{1}{2}}$.

The normalized transform $\mathbf{F}_{N}$ has the closed form

[TABLE]

This is a variant of the Discrete Cosine Transform (DCT-III, or the inverse DCT). It is well-known that the DCT-III can be ‘embedded’ into a DFT of twice the dimension, and we work out some of the details of how to use the sparse FFT to compute a sparse DCT in Appendix A.

Legendre polynomials.

Legendre polynomials are orthogonal with respect to the uniform measure i.e. $w(X)=1$ and play a critical role in multipole expansions of potential functions (whether electrical or gravitational) in spherical coordinates. They are also important for solving Laplace’s equation in spherical coordinates.

2.2.2 Roots of Orthogonal Polynomials

Since $\lambda_{i}\in[-1,1]$ for all $i$, there is a unique $\theta_{i}\in[0,\pi]$ such that $\lambda_{i}=\cos\theta_{i}$. Our reduction holds for orthogonal polynomials that have roots that are ‘well-separated’ in this angle space:

Definition 2.2**.**

Let $0<C_{0}<C_{1}$. A family of orthogonal polynomials $p_{0}(X),p_{1}(X),\ldots$ is $(C_{0},C_{1},\gamma_{0})$-dense if for all large enough $d$, the following holds.

Let $\lambda_{0},\ldots,\lambda_{d-1}$ be the roots of $p_{d}$, and $\theta_{i}=\arccos\lambda_{i}$. Then for any $i\in[d]$, for any $\gamma\geq\gamma_{0}/d$:

[TABLE]

It turns out that any family of Jacobi polynomials has the required property: their roots are spaced out such that $\theta_{\ell}$ is close to $\ell\pi/N$ (Theorem 5.2).

2.3 Sparse Recovery Problem

We will consider approximately $k$-sparse vectors $\hat{\mathbf{v}}=\hat{\mathbf{x}}+\hat{\mathbf{w}}$, where $\hat{\mathbf{x}}$ is $k$-sparse and $\|\hat{\mathbf{w}}\|_{2}$ is sufficiently small. We will require that $\hat{\mathbf{x}}$ has a ‘spread out’ support, defined as follows.

Definition 2.3**.**

Let $k\in[N]$ and $0\leq\sigma<1$. We say that a vector $\mathbf{x}\in\mathbb{R}^{N}$ is $(k,\sigma)$-sparsely separated if there are $k$ non-zero locations in $\mathbf{x}$ and any two non-zero locations are more than $\sigma N$ indices apart.

It is not hard to see that a vector $\mathbf{x}$ with random support of size $k$ is, with constant probability, $\left({k,\Omega\left({\frac{1}{k^{2}}}\right)}\right)$-sparsely separated.

In our reduction, we will reduce the $k$-sparse recovery problem to the special case of $k=1$. Next, we define some notation for the $1$-sparse case.

Definition 2.4**.**

We say that the matrix $\mathbf{F}_{N}$ has an $(N,\epsilon,\delta,\mu)$ one-sparse recovery algorithm with query complexity $Q(N,\epsilon,\delta,\mu)$ and time complexity $T(N,\epsilon,\delta,\mu)$ if there exists an algorithm $\mathcal{A}$ with the properties below:

For all $\mathbf{y}$ so that $\hat{\mathbf{y}}=\mathbf{F}_{N}\mathbf{y}$ can be decomposed as

[TABLE]

where ${\tilde{\mathbf{y}}}=v\cdot\mathbf{e}_{h}$ is $1$-sparse and

[TABLE]

we have:

$\mathcal{A}$ makes at most $Q(N,\epsilon,\delta,\mu)$ queries into $\mathbf{y}=\mathbf{F}^{-1}\left({v\cdot\mathbf{e}_{h}+\mathbf{w}}\right)$. 2. 2.

With probability at least $1-\mu$, $\mathcal{A}$ outputs $\tilde{v}\cdot\mathbf{e}_{h}$ with $|v-\tilde{v}|\leq\delta\left|v\right|$ in time $T(N,\epsilon,\delta,\mu)$.

Pre-processing time.

Our algorithm requires some pre-processing of $\mathbf{F}_{N}$. Our pre-processing step is given in Algorithm 3 and involves computing the roots $\lambda_{1},\ldots,\lambda_{N}$ of $p_{N}$ and storing them in an appropriate data structure, and additionally forming and storing some matrices that we will use in our algorithm. Finding the roots and creating the data structure can be done in time $\mathrm{poly}(N)$, and the rest of the pre-processing step also takes time $\mathrm{poly}(N)$. We note that this is an up-front cost that needs to be only paid once.

Precision.

We note that we need to make certain assumptions on size of the entries in $\hat{\mathbf{v}}$ since otherwise we would not even be able to read coefficients that are either too large or too small and need $\omega(\log{N})$ bits to represent. Towards this end we will make the standard assumption that $\left\|{\hat{\mathbf{v}}}\right\|_{{2}}=1$. In particular, this allows us to ignore any coefficients that are smaller than say $\frac{1}{N}$ since their contribution to $\left\|{\hat{\mathbf{v}}}\right\|_{{2}}$ is at most $\frac{1}{\sqrt{N}}$, which will be too small for our purposes.777More generally, we can ignore smaller coefficients as long as they are polynomially large. In particular, this implies that we only have to deal with numbers that need $O(\log{N})$ bits and as is standard in the RAM model, basic arithmetic operations on such numbers can be done in $O(1)$ time. We will implicitly assume this for the rest of the paper (except in the proof of Lemma 5.21, where we will explicitly make use of this assumption).

3 Results

In this section we state our main results. These results follow from more detailed versions which are stated with the proofs of these results.

We start off with our main result for Jacobi polynomials. We state an informal version here, and refer the reader to Corollary 6.2 for the formal result.

Theorem 3.1** (General Sparse Recovery for Jacobi Polynomial Transform, Informal).**

Fix arbitrary parameters $\alpha,\beta>-1$ for Jacobi polynomials and let $\mathbf{J}^{(\alpha,\beta)}_{N}$ be the $N\times N$ orthogonal matrix that arises from it as in Definition 2.1. Then there is an algorithm Recover that does the following. Let $\mathbf{v}=\mathbf{x}+\mathbf{w}$ where $\hat{\mathbf{x}}=\mathbf{J}^{(\alpha,\beta)}_{N}\mathbf{x}$ is $(k,C_{1}/k^{2})$-sparsely separated, and suppose that $\|\hat{\mathbf{w}}\|_{2}\lesssim\delta\min_{h\in\mathrm{supp}(\hat{\mathbf{x}})}|\hat{\mathbf{x}}[h]|.$ Then with probability at least $0.99$, Recover outputs $\hat{\mathbf{z}}$ such that

[TABLE]

with $\mathrm{poly}\left({\frac{k\log{N}}{\delta}}\right)$ queries and running time $\mathrm{poly}\left({\frac{k\log{N}}{\delta}}\right)$.

Remark 3.2**.**

The requirement on the noise term might be bad if one entry of $\hat{\mathbf{x}}$ is extremely small compared to the rest. However in this case we can decrease $k$ and add the very small entries of $\hat{\mathbf{x}}$ to the noise term $\hat{\mathbf{w}}$ resulting in a potentially better guarantee. We note that our algorithm iteratively finds the large components of $\hat{\mathbf{x}}$ and in fact has a mechanism for stopping early when all of the ‘large-enough’ entries have been found.

To prove the above result, we first reduce the $k$-sparse recovery problem to $1$-sparse recovery problem, in the presence of a small amount of noise. Next, we present an informal statement of our reduction. See Theorem 4.1 for the formal result.

Theorem 3.3** (Main Reduction, Informal).**

Let $p_{1},\ldots,p_{N}$ be a $(C_{0},C_{1},\gamma_{0})$-dense orthogonal polynomial family, and let $\mathbf{F}_{N}$ be the $N\times N$ orthogonal matrix that arises from it as in Definition 2.1. Suppose that $|\mathbf{F}^{-1}_{N}[i,j]|\lesssim 1/\sqrt{N}$ for all $i,j\in[N]$. Suppose that for some sufficiently small $\delta>0$, $\mathbf{F}_{N}$ has a $\left({N,O(\delta),\delta,O(C_{0}/k^{2})}\right)$ one-sparse recovery algorithm with query complexity $Q$ and running time $T$.

Then there is an algorithm Recover that does the following. Let $\mathbf{v}=\mathbf{x}+\mathbf{w}$ where $\hat{\mathbf{x}}=\mathbf{F}_{N}\mathbf{x}$ is $(k,C_{1}/k^{2})$-sparsely separated, and suppose that $\|\hat{\mathbf{w}}\|_{2}\lesssim\delta\min_{h\in\mathrm{supp}(\hat{\mathbf{x}})}|\hat{\mathbf{x}}[h]|.$ Then with probability at least $0.99$, Recover outputs $\hat{\mathbf{z}}$ so that

[TABLE]

with $\mathrm{poly}(k/\delta C_{0})Q$ queries and running time $\mathrm{poly}(k/\delta C_{0})T$.

The final algorithmic piece missing from the result above is the algorithm for $1$-sparse recovery. We provide this missing piece for Jacobi polynomials (see Theorem 5.1 for the formal statement):

Theorem 3.4** ($1$-Sparse Recovery for Jacobi Transform, Informal).**

There exists a universal constant $C$ such that the following holds. Consider the Jacobi transform for any fixed parameters $\alpha,\beta>-1$. There exists an $(N,\epsilon,C\cdot\epsilon,\gamma)$ $1$-sparse recovery algorithm for the Jacobi transform that makes $\mathrm{poly}\left({\log\left({\frac{N}{\gamma}}\right)\cdot\frac{1}{\epsilon}}\right)$ queries and takes time $\mathrm{poly}\left({\log\left({\frac{N}{\gamma}}\right)\cdot\frac{1}{\epsilon}}\right)$.

3.1 Open Questions

Before we dive into the proofs of the results above, we list a few questions left open by our work.

First, it is a natural to try and improve our $k$-sparse recovery algorithm to work for arbitrary $k$-sparse support, rather than ‘well-separated’ supports. One natural way to do this is to address the fourth (open) challenge in Section 1.3 for a general class of OPs. 2. 2.

Second, we could hope to handle a more general class of noise $\hat{\mathbf{w}}$ than we currently do. One could hope to handle *any *vector $\mathbf{v}$, with an error guarantee that degrades smoothly with the $\ell_{2}$ norm of the ‘tail’ of $\mathbf{v}$. 3. 3.

Third, we would like to extend our results to hold for OPs defined over infinite intervals (e.g. Hermite and Laguerre polynomials). 4. 4.

Fourth, we would like to solve the sparse recovery for $\mathbf{F}^{T}$ (where $\mathbf{F}$ is as in Definition 2.1): i.e. given query access to $\mathbf{x}$ figure out a good $k$-sparse approximation to $\mathbf{F}^{T}\mathbf{x}$ (recall that $\mathbf{F}^{-1}=\mathbf{F}^{T}$). (Note that this problem can be equivalently stated as follows: given query access to $\mathbf{F}\mathbf{y}$, compute a good $k$-sparse approximation to $\mathbf{y}$.) Currently our results do not solve this problem since we cannot show that the existence of a filter $\mathbf{b}$ such that $\mathbf{F}\mathbf{D}_{\mathbf{b}}\mathbf{F}^{T}$ is row-sparse. Note that this is not an issue for DFT since it is symmetric. 5. 5.

Finally, we would like to reduce the exponent on $k$ in our final runtime. In particular, for the case of random $k$-sparse support, the dependence on $k$ in the runtime for Jacobi transform is $k^{8}$. We note that we have not tried too hard to optimize the constants though we believe even getting a quadratic dependence on $k$ with our framework would be challenging. We would like to stress that the majority of the work in the sFFT literature has been to make the dependence on $k$ be linear and for such results, it seems very unlikely that a generic reduction from $k$-sparse recovery to $1$-sparse recovery would work. In other words, using the knowledge about the $1$-sparse recovery algorithm for DFT seems necessary to get a overall $k$-sparse FFT with running time $k\mathrm{poly}(\log{n})$.

4 Reduction to $1$-sparse case

In this section, we will prove Theorem 3.3, which shows how we can reduce the $k$-sparse recovery problem (when $\mathbf{x}$ is also ‘well-separated’) to the $1$-sparse case. Theorem 3.3 follows from the following theorem, which is the main theorem in this section. (Recall that $\textsc{large}(s,\mathbf{x})$ is the $s$’th largest value in $\mathbf{x}$ by absolute value.)

Theorem 4.1**.**

There are some constants $\delta_{0},C^{\prime}>0$ so that the following holds. Let $p_{1},\ldots,p_{N}$ be a $(C_{0},C_{1},\gamma_{0})$-dense orthogonal polynomial family, and let $\mathbf{F}_{N}$ be the $N\times N$ orthogonal matrix that arises from it as in Definition 2.1. Let

[TABLE]

Then there is an algorithm Recover that does the following.

Consider any $k\in[N]$, $\delta<\delta_{0}$, $0<\mu<1$, and $\gamma\geq\gamma_{0}/N$. Let

[TABLE]

Suppose that $\mathbf{F}_{N}$ has a $\left({N,\frac{6\delta}{C},\delta,\frac{\mu_{0}}{2}}\right)$ one-sparse recovery algorithm ${\cal A}$ for some $C\geq C^{\prime}$ with query complexity $Q\left({N,\frac{6\delta}{C},\delta,\frac{\mu_{0}}{2}}\right)$ and time complexity $T\left({N,\frac{6\delta}{C},\delta,\frac{\mu_{0}}{2}}\right)$.

Let $\hat{\mathbf{v}}=\hat{\mathbf{x}}+\hat{\mathbf{w}}$ so that $\hat{\mathbf{x}}=\mathbf{F}_{N}\mathbf{x}$ is $(k,C_{1}\gamma)$-sparsely separated, and so that $\left\|{\hat{\mathbf{w}}}\right\|_{{2}}\leq\frac{\delta}{2C}\textsc{large}(k,\hat{\mathbf{x}})$.

Then with probability at least $1-\mu$, $\textsc{Recover}^{(\mathbf{v})}(k,\delta,\mu,\gamma)$ outputs $\hat{\mathbf{z}}$ so that

[TABLE]

and so that

[TABLE]

Further, Recover makes at most

[TABLE]

queries to $\mathbf{v}$, and has running time

[TABLE]

The basic idea of Recover (Algorithm 1) is as follows. We will define an algorithm Peeler (Algorithm 2) which will iteratively ‘peel’ off the heavy hitters and store them in the approximation $\hat{\mathbf{z}}$ until the variable stop is set to true. We will show that if we stop, then every value in the residual $\hat{\mathbf{v}}-\hat{\mathbf{z}}$ will be very small, at which point we will be done. The idea for the Peeler algorithm is illustrated in Figure 2: we first use a filter to hopefully isolate a single spike of $\hat{\mathbf{x}}$, and then we use the one-sparse recovery algorithm to estimate this spike and subtract it off.

In Section 4.1 below, we prove Theorem 4.1 assuming that a suitable Peeler algorithm exists. More precisely, we will formalize what we need from Peeler in Lemma 4.2, and we will prove Theorem 4.1 assuming Lemma 4.2. Next, in Sections 4.2 and 4.3, we will state the Peeler algorithm (Algorithm 2) and prove that it works. Section 4.2 contains some useful preliminaries, and Section 4.3 contains the statement of Peeler and the proof of Lemma 4.2.

4.1 Proof of Theorem 4.1

In this section we prove Theorem 4.1 which implies Theorem 3.3. For the rest of this section, given a vector $\mathbf{x}\in\mathbb{R}^{N}$, we will denote its transform $\mathbf{F}_{N}\cdot\mathbf{x}$ by $\hat{\mathbf{x}}$.

The algorithm Recover (Algorithm 1) uses an algorithm called Peeler (Algorithm 2) multiple times. This algorithm will satisfy the following guarantee.

Lemma 4.2**.**

Consider an orthogonal polynomial family and its corresponding matrix $\mathbf{F}_{N}$ that satisfies the properties in Theorem 4.1.

There is an algorithm Peeler with the following guarantee.

Let $k\in[N]$, $\epsilon$ be sufficiently small, $0<\mu<1$, $\gamma\geq\gamma_{0}/N$, and $\mu_{0}$ as in (4). Suppose that $\mathbf{F}_{N}$ has a $(N,6\epsilon,C\epsilon,\mu_{0}/2)$ one-sparse recovery algorithm ${\cal A}$ for some sufficiently large $C$.

Suppose Peeler has query access to $\mathbf{v}=\mathbf{x}+\mathbf{w}$ such that $\hat{\mathbf{x}}$ is $\left({k,C_{1}\gamma}\right)$-sparsely separated. Consider an input $\hat{\mathbf{z}}$ with sparsity $\|\hat{\mathbf{z}}\|_{0}<k$ such that $\|\hat{\mathbf{w}}\|_{2}\leq\epsilon\|\hat{\mathbf{v}}-\hat{\mathbf{z}}\|_{\infty}$ and the following holds for every $i\in\mathrm{supp}(\hat{\mathbf{z}})$:

[TABLE]

Then, with probability at least $1-\mu$, Peeler${}^{(\mathbf{v})}(\hat{\mathbf{z}},k,\epsilon,\mu,\gamma)$ returns $(\tilde{\mathbf{z}},\text{{stop}})$, where $\tilde{\mathbf{z}}\in\mathbb{R}^{N}$ and stop is a Boolean variable, so that

• •

If stop is false, then $\tilde{\mathbf{z}}=\hat{\mathbf{z}}+\tilde{v}\cdot\mathbf{e}_{h}$ for some $h\in\mathrm{supp}(\hat{\mathbf{x}})\setminus\mathrm{supp}(\hat{\mathbf{z}})$, and (7) is satisfied for all $i\in\mathrm{supp}(\tilde{\mathbf{z}})$, i.e.

[TABLE]

• •

If stop is true, then $\mathrm{supp}(\tilde{\mathbf{z}})=\mathrm{supp}(\hat{\mathbf{z}})$, (7) is satisfied for all $i\in\mathrm{supp}(\tilde{\mathbf{z}})$, and $\|\hat{\mathbf{v}}-\tilde{\mathbf{z}}\|_{\infty}\leq\frac{\epsilon}{\sqrt{k}}\|\hat{\mathbf{v}}\|_{\infty}$.

Further, if ${\cal A}$ has query complexity $Q(N,6\epsilon,C\epsilon,\mu_{0}/2)$ and running time $T(N,6\epsilon,C\epsilon,\mu_{0}/2)$, Peeler makes at most

[TABLE]

queries to $\mathbf{x}+\mathbf{w}$, and runs in time at most

[TABLE]

where $\mu_{0},T_{0},T_{1},T_{2}$ are as defined in Algorithm 2.

The basic idea of Lemma 4.2 is that Peeler will iteratively ‘peel’ off the heavy hitters and store them in the approximation $\hat{\mathbf{z}}$ until the variable stop is set to true. As long as we are doing well enough so far, the lemma says that we will find a new heavy hitter. If we stop, then the lemma implies that every value in the residual $\hat{\mathbf{v}}-\hat{\mathbf{z}}$ is very small.

We will present Peeler and prove Lemma 4.2 in Sections 4.2 and 4.3. Before that, we prove Theorem 4.1—which says that Recover (Algorithm 1) works—assuming Lemma 4.2.

Proof of Theorem 4.1.

We will repeatedly apply Lemma 4.2 with $\epsilon=\frac{\delta}{C}$ (and preserving the remaining parameters); for the rest of this section, we use this choice for $\epsilon$. As a preliminary, we note the following useful bounds between $\hat{\mathbf{x}}$ and $\hat{\mathbf{v}}$ at the locations of the heavy hitters (or spikes), due to the condition on the noise $\hat{\mathbf{w}}$ in the statement of the theorem. More precisely, we have the following claim.

Claim 4.3**.**

For any $i\in\mathrm{supp}(\hat{\mathbf{x}})$,

[TABLE]

[TABLE]

Furthermore if $v$ satisfies $\left|v-\hat{\mathbf{v}}[i]\right|\leq C\epsilon\left|\hat{\mathbf{v}}[i]\right|$, then

[TABLE]

Proof.

For any $i\in\mathrm{supp}(\hat{\mathbf{x}})$, the noise condition says

[TABLE]

which implies (8). This also implies that

[TABLE]

which in turn implies (9).

Consequently, if $v$ satisfies $\left|v-\hat{\mathbf{v}}[i]\right|\leq C\epsilon\left|\hat{\mathbf{v}}[i]\right|$, then

[TABLE]

∎

Conditions to Lemma 4.2.

We first check that at every call to Peeler in Recover, the conditions to Lemma 4.2 are always met with this choice of $\epsilon$. First we note that the invariant (7) is vacuously satisfied at the beginning of Recover, and the guarantees of Peeler imply that it is therefore always satisfied. The only condition left to check is that the requirement on the noise $\hat{\mathbf{w}}$ in Lemma 4.2 is met: that is, that $\|\hat{\mathbf{w}}\|_{2}\leq\epsilon\|\hat{\mathbf{v}}-\hat{\mathbf{z}}\|_{\infty}$. Let $h\in\mathrm{supp}(\hat{\mathbf{x}})\setminus\mathrm{supp}(\hat{\mathbf{z}})$, which exists since Peeler is called at most $k$ iterations. By our assumption on the noise in Theorem 4.1 and (9),

[TABLE]

Establishing the error guarantee.

Next, we claim that the estimates $\hat{\mathbf{z}}$ returned by Peeler are always good in an $\ell_{2}$ sense.

Claim 4.4**.**

Let $\hat{\mathbf{z}}$ be the output of Peeler at any iteration. Let $Z$ be the support of $\hat{\mathbf{z}}$ and $X$ be the support of the largest $|Z|$ coordinates of $\hat{\mathbf{x}}$. Then

[TABLE]

Proof.

We have

[TABLE]

where we have used the guarantee from Peeler (7) as well as (10). Now suppose that $i\in Z\setminus X$ and $j\in X\setminus Z$. By definition of $X$, $\left|\hat{\mathbf{x}}[i]\right|\leq\left|\hat{\mathbf{x}}[j]\right|$. On the other hand, we have

[TABLE]

where the first and the last inequality follows from Claim 4.3 and the second inequality follows from (7) and since $\left\|{\hat{\mathbf{v}}-\hat{\mathbf{z}}}\right\|_{{\infty}}\geq\left|\hat{\mathbf{v}}[j]-\hat{\mathbf{z}}[j]\right|=\left|\hat{\mathbf{v}}[j]\right|$. Since $|Z\setminus X|=|X\setminus Z|$, the first term in (11) can be bounded by

[TABLE]

using the fact that $C$ is a sufficiently large absolute constant, and that $\epsilon$ is sufficiently small compared to $C$. Adding in the second term in (11), we have

[TABLE]

as desired.

∎

Suppose that Peeler never returns $\texttt{stop}=\textnormal{\sc true}$ in Recover. In this case, the final output $\hat{\mathbf{z}}$ has support of size $k$ so $X=\mathrm{supp}(\hat{\mathbf{x}})$ in Claim 4.4 and we get

[TABLE]

which immediately establishes the first error bound in the theorem.

On the other hand, suppose some call to Peeler returns $\texttt{stop}=\textnormal{\sc true}$. Let $\hat{\mathbf{z}}$ denote the output of Peeler on this final step, $s=\|\hat{\mathbf{z}}\|_{0}$, and let $X$ denote the support of the $s$ largest entries of $\hat{\mathbf{x}}$. Let $h$ be the index of the $(s+1)$’th largest entry in $\hat{\mathbf{x}}$, and let $h^{\prime}$ be the largest entry in $\hat{\mathbf{x}}_{\mathrm{supp}(\tilde{\mathbf{z}})^{c}}$. Note that $h^{\prime}\in\mathrm{supp}(\hat{\mathbf{x}})$ since the algorithm stopped before hitting all $k$ heavy hitters. Now

[TABLE]

In the above, the second line we have used Claim 4.3. In the second to last line we have used the guarantee from Lemma 4.2, and in the last line we use the fact that $\arg\max_{i}|\hat{\mathbf{v}}[i]|\in\mathrm{supp}(\hat{\mathbf{x}})$ which follows from the assumption on $\hat{\mathbf{w}}$ and Claim 4.3.

Thus by Claim 4.4 we have

[TABLE]

provided that $C$ is sufficiently large. This establishes the first error bound (5) in the theorem statement. Now the second error bound (6) follows from the triangle inequality.

Query and Time complexity of the Recover algorithm.

Finally, we compute the query complexity and running time. With our choices of $T_{0},T_{1}$ and $T_{2}$, it can be checked that Peeler uses query complexity

[TABLE]

and has running time

[TABLE]

Since the work of Recover is dominated by running Peeler at most $k$ times, the same query complexity and running time bounds hold for Recover (since the extra factor of $k$ can be absorbed into the $\mathrm{poly}(k)$ factor). ∎

4.2 Preliminaries to the proof of Lemma 4.2

In this section we present a few results which will be useful in the proof of Lemma 4.2. We begin by showing that $(C_{0},C_{1},\gamma_{0})$-dense orthogonal polynomial families (see Definition 2.2) and $(k,C_{1}\gamma)$-sparsely separated vectors (see Definition 2.3) work well together in the following sense.

Lemma 4.5**.**

Let $H$ be the support of any $(k,C_{1}\gamma)$-sparsely separated vector $\hat{\mathbf{x}}$ where $\gamma\geq 2\gamma_{0}/N$. Let $h\in H$. Pick $0\leq\ell<N$ uniformly at random. Then

• (i)

With probability at least $C_{0}\gamma/2$, we have

[TABLE]

and

• (ii)

Conditioned on (12), the following holds with probability $1$:

[TABLE]

Proof.

We start with part (i). Indeed note that for (12) to occur, we must have

[TABLE]

Since the orthogonal polynomial family is $(C_{0},C_{1},\gamma_{0})$-dense, the above happens for at least $\frac{C_{0}\gamma N}{2}$ values of $\ell$. Since $\ell$ was chosen randomly, the above holds with probability $\frac{C_{0}\gamma}{2}$. Once (12) holds, we show that (13) follows from $\hat{\mathbf{x}}$ being $(k,C_{1}\gamma)$-separated. Since the orthogonal polynomial family is $(C_{0},C_{1},\gamma_{0})$-dense, the number of $\ell$ such that

[TABLE]

is at most $C_{1}\gamma N$. In particular, any such $\ell$ satisfies $|\ell-h|\leq C_{1}\gamma N$. Since $\hat{\mathbf{x}}$ is $(k,C_{1}\gamma)$-separated none of these can be contained in $H\setminus\{h\}$, as desired. ∎

The Peeler algorithm relies on *boxcar polynomials, *which we define as follows (see Figure 1).

Definition 4.6**.**

Let $d\geq 1$ be an integer and $0<\gamma\leq\frac{\pi}{2},\theta\in[0,\pi],\epsilon\geq 0$ be reals. We say that $p(X)$ is a $(d,\epsilon,\theta,\gamma)$-box car polynomial if the following are true:

1. (a)

$p(X)$ has degree $d$, 2. (b)

$|p(\cos\phi)|\leq\epsilon$ for any $\phi\in[0,\pi]\setminus[\theta-2\gamma,\theta+2\gamma]$, 3. (c)

$|p(\cos\phi)-1|\leq\epsilon$ for any $\phi\in\left[{\theta-\gamma,\theta+\gamma}\right]$, and 4. (d)

$|p(\cos\phi)|\leq 1+\epsilon$ for the remaining values of $\phi$.

Fortunately for us, such polynomials exist with reasonably low degree:

Lemma 4.7**.**

For any $\theta\in[0,\pi]$, there exists an $\left({O\left({\frac{1}{\gamma\epsilon}}\right),\epsilon,\theta,\gamma}\right)$-boxcar polynomial.

To prove Lemma 4.7, we will need the following result from approximation theory [R81, jackson]:

Theorem 4.8** (Jackson’s Theorem).**

There is a constant $C$ such that: if $f:[0,2\pi]\to\mathbb{C}$ is a periodic function with Lipschitz constant $L$, then for every $d$ there is a trigonometric polynomial $T$ of degree at most $d$ such that

[TABLE]

Furthermore if $f$ is even then $T$ is even.

Recall that a trigonometric polynomial of degree $d$ has the form $a_{0}+\sum_{m=1}^{d}a_{m}\cos(mx)+\sum_{m=1}^{d}b_{m}\sin(mx)$, and it is even if all $b_{m}=0$. In particular, since $\cos(mx)=T_{m}(\cos x)$ where $T_{m}$ is the $m$th Chebyshev polynomial, an even trigonometric polynomial $T$ of degree $d$ in $x$ is a polynomial of degree $d$ in $\cos x$, in other words $T\circ\arccos$ is a polynomial of degree $d$. Now we can prove Lemma 4.7.

Proof of Lemma 4.7.

Consider the following piecewise linear function $f$ such that

[TABLE]

and such that $f$ on $[\pi,2\pi]$ is reflected from $[0,\pi]$ so that $f$ is even. Clearly $f$ has Lipschitz constant $\frac{1}{\gamma}$. Let $d$ satisfy $\frac{C}{d}\frac{1}{\gamma}\leq\epsilon$ and $T$ be a trigonometric polynomial of degree $d$ satisfying Theorem 4.8. Then (14) directly implies conditions (b), (c) and (d) for the function $p=T\circ\arccos$, which is a polynomial of degree $d$ as noted before. ∎

4.3 Proof of Lemma 4.2

In this section we prove Lemma 4.2. Fix a large enough $N$. For the rest of this section, for notational convenience we will use $\mathbf{F}=\mathbf{F}_{N}$.

We begin by describing the algorithm Peeler, which requires setting up some notation.

For every $0\leq\ell<N$, let $b_{\ell}(X)$ be an $\left({d,\frac{\epsilon}{\sqrt{k}},\theta_{\ell},\frac{\gamma}{4}}\right)$-boxcar polynomial, for some $d$ which we will choose to be

[TABLE]

Define $\mathbf{b}_{\ell}=\left({b_{\ell}(\lambda_{1}),\dots,b_{\ell}(\lambda_{N-1})}\right)$.

Given any $\tau\geq 0$, we define the following ‘truncation’ function that is defined for any $x\in\mathbb{R}$:

[TABLE]

For notational convenience, we will let $\mathbf{r}=\mathbf{v}-\mathbf{z}$ and $\hat{\mathbf{r}}=\hat{\mathbf{v}}-\hat{\mathbf{z}}$ be the residual vector.

Given this notation, the algorithm Peeler to compute $\tilde{\mathbf{z}}$ is presented in Algorithm 2, and is illustrated in Figure 2. Next, we go through the various components of Algorithm 2 and prove that they work as intended.

4.3.1 Simulating Query Access

We begin with the analysis of SimulateQueryAccess (Alg. 4), which allows Peeler (Alg. 2) to simulate query access to $\mathbf{F}^{-1}\hat{\mathbf{y}}_{\ell}$. In Peeler, this is needed both to run the one-sparse recovery algorithm, and to run the Verify algorithm.

We first prove a general property about low-degree polynomials:

Lemma 4.9**.**

Let $b(X)$ be any polynomial of degree $d$. Then for all $0\leq j<N$, the only values $0\leq i<N$ so that

[TABLE]

are those so that $|i-j|\leq d$. In particular, if we define $\mathbf{b}=\left({b(\lambda_{0}),\dots,b(\lambda_{N-1})}\right)$, then each row $j$ of the matrix $\mathbf{F}^{T}\mathbf{D}_{\mathbf{b}}\mathbf{F}$ has at most $O(d)$ non-zero values in it, at positions $i$ so that $|i-j|\leq d$.

Proof.

Since the orthogonal polynomials $\{p_{0},p_{1},p_{2},\ldots,p_{r}\}$ form a basis for polynomials of degree at most $r$ (this follows since $\deg(p_{i})=i$), it follows from orthogonality conditions that for any $r$ and for any polynomial $f(X)$ of degree strictly less than $r$, we have

[TABLE]

Then the above implies that

[TABLE]

is zero whenever $|i-j|>d$. Further, (2) implies that we have

[TABLE]

whenever $i+j+d\leq 2N-1$. Thus, (15) only holds if either $|i-j|\leq d$ or if $i+j+d\geq 2N$.

Because $i,j\in\{0,\ldots,N-1\}$, the only way that $i+j+d\geq 2N$ is if $|i-j|\leq d-2$, which is already covered by the above. ∎

With this fact out of the way, we observe that SimulateQueryAccess works as intended.

Proposition 4.10**.**

Each call to SimulateQueryAccess (Alg. 4) uses $O(d)$ queries to $\mathbf{x}+\mathbf{w}=\mathbf{F}^{-1}(\hat{\mathbf{x}}+\hat{\mathbf{w}})$, runs in time $O(kd)$, and returns $\hat{\mathbf{y}}_{\ell}[j]$.

Proof.

As in Algorithm 2, define

[TABLE]

Thus,

[TABLE]

where $\mathbf{z}=\mathbf{F}^{-1}\hat{\mathbf{z}}$ and in above we have used the definition of $\hat{\mathbf{y}}_{\ell}$ and the fact that $\mathbf{F}^{-1}=\mathbf{F}^{T}$.

Now we have

[TABLE]

which by Lemma 4.9 is only nonzero if $|i-j|\leq d$, since $d$ is the degree of $b_{\ell}$. Further, we can expand $b_{\ell}(\lambda_{t})=\sum_{r=0}^{d}b_{\ell,r}\lambda_{t}^{r}$ and observe that when (17) is nonzero, then it is equal to $\sum_{r=0}^{d}b_{r,j}\mathbf{M}_{r}[j,i]$. Thus, we have

[TABLE]

since this is exactly how we have set $\nu_{i}$ in Algorithm 4. Therefore by (16) we have

[TABLE]

which is what is returned in SimulateQueryAccess (Alg. 4).

In order to compute $\sum_{i}\nu_{i}(\mathbf{x}[i]+\mathbf{w}[i]-\mathbf{z}[i])$ we need $O(d)$ queries to $\mathbf{x}$. We also need $O(d)$ queries to $\mathbf{z}$ which we can compute directly from $\hat{\mathbf{z}}$ in time $O(k)$ per query (using the fact that $\hat{\mathbf{z}}$ is $k$-sparse, and that $\mathbf{z}=\mathbf{F}^{-1}\hat{\mathbf{z}}$). This proves the proposition. ∎

4.3.2 Correctness of the Peeler algorithm

In the rest of the proof of Lemma 4.2, we argue that $\tilde{\mathbf{z}}=\textsc{Peeler}^{(\mathbf{v})}(\hat{\mathbf{z}},k,\epsilon,\mu,\gamma)$ satisfies the required properties.

Recalling that $\hat{\mathbf{r}}=\hat{\mathbf{v}}-\hat{\mathbf{z}}$, let $h^{*}$ be the location of the largest magnitude entry in $\hat{\mathbf{r}}$. To prove the correctness of the Peeler algorithm, we will need the following two lemmas (which we will prove later in Sections 4.3.5 and 4.3.6) that reason about any iteration of the inner loop on line 11 in Algorithm 2. We first argue that a ‘good spike’ (that is, a value of $h$ so that $|\hat{\mathbf{r}}[h]|$ is large) will pass the Verify check and hence the largest heavy hitter in $\hat{\mathbf{r}}=\hat{\mathbf{v}}-\hat{\mathbf{z}}$ will be included in the set $\mathcal{L}$:

Lemma 4.11** (Good spikes get noticed).**

Let $\epsilon$ be sufficiently small compared to $C$. Let $h\in[N]$. Suppose that $\hat{\mathbf{r}}[h]$ satisfies $\left|\hat{\mathbf{r}}[h]\right|\geq\frac{1}{2}\cdot\|\hat{\mathbf{r}}\|_{\infty}$ and that $\hat{\mathbf{v}}=\hat{\mathbf{x}}+\hat{\mathbf{w}}$ so that $\left\|{\hat{\mathbf{w}}}\right\|_{{2}}\leq\epsilon\left\|{\hat{\mathbf{r}}}\right\|_{{\infty}}$. Also suppose that $\theta_{h}\in\left[{\theta_{\ell}-\frac{\gamma}{4},\theta_{\ell}+\frac{\gamma}{4}}\right]$ (line 16). Then in Peeler, with probability $1-\mu_{0}$ (over $\mathcal{A}$ (line 15) and Verify (line 17)), $(\tilde{v},h)$ gets added to $\mathcal{L}$ such that (where $v=\hat{\mathbf{r}}[h]$):

[TABLE]

Next, we will argue that ‘bad spikes’ get caught by Verify and hence they will not prevent the actual heavy hitter from being chosen:

Lemma 4.12** (Bad spikes get pruned out).**

Let $\epsilon$ be sufficiently small compared to $C$. Let $h\in[N]$. Suppose that $\hat{\mathbf{r}}[h]$ satisfies $\left|\hat{\mathbf{r}}[h]\right|<\frac{1}{2}\left\|{\hat{\mathbf{r}}}\right\|_{{\infty}}$ and that $\hat{\mathbf{v}}=\hat{\mathbf{x}}+\hat{\mathbf{w}}$ so that $\left\|{\hat{\mathbf{w}}}\right\|_{{2}}\leq\epsilon\left\|{\hat{\mathbf{r}}}\right\|_{{\infty}}$. Then with probability at least $1-\mu_{0}$, if $\left|y_{\ell}\right|\geq(1-C\epsilon)\left\|{\hat{\mathbf{v}}-\hat{\mathbf{z}}}\right\|_{{\infty}}$ then $(y_{\ell},h)$ is not added to $\mathcal{L}$.

Before we prove either of the two lemmas above, we first use them to argue the correctness of Peeler algorithm (as claimed in Lemma 4.2). Recall that by assumption,

[TABLE]

for all $i\in\mathrm{supp}(\hat{\mathbf{z}})$.

Claim 4.13**.**

In every iteration of the outer **For ** loop (line 10 in Algorithm 2), with probability at least $1-\mu/T_{2}$,

[TABLE]

and

[TABLE]

where $(\bar{v},\bar{h})\in\mathcal{L}$ is the pair picked by the Peeler algorithm after the inner For loop has completed (line 19 in Algorithm 2). (Above, $\hat{\mathbf{r}}$ is the value of $\hat{\mathbf{v}}-\hat{\mathbf{z}}$ at the beginning of the loop.)

Proof.

Let us consider a single iteration of the outer For loop.

We first note that the conditions to Lemma 4.11 and  4.12 cannot both be satisfied, and so exactly one of them is invoked in the analysis each time Verify is called. By a union bound over all $T_{0}$ iterations of the inner loop, we can say that with probability $1-T_{0}\mu_{0}$ (over all randomness in lines 15 and 17), then Lemmas 4.11 and 4.12 have the favorable outcome every time they are invoked in one iteration of the outer loop. By our choices of $T_{0},T_{2},\mu_{0}$, we have $\mu_{0}\leq\frac{\mu}{2T_{0}T_{2}}$. Indeed, plugging in our choices, we would like to show that

[TABLE]

for which it suffices to show

[TABLE]

Since $\mu_{0}=\left({\frac{\mu C_{0}\gamma}{k}}\right)^{2}$, this is true (assuming $k$ is large enough). Thus, the probability that Lemma 4.11 and Lemma 4.12 have the favorable outcome every time they are invoked in one iteration of the outer loop is at least $1-\frac{\mu}{2T_{2}}$.

Suppose that this occurs, and now consider the ‘spike’ at $h^{*}=\mathrm{argmax}_{i}|\hat{\mathbf{r}}[i]|$. By Lemma 4.5 and our choice of $T_{0}$, with probability $1-\frac{\mu}{2T_{2}}$ (over all randomness in line 12) there is some iteration $t$ in the inner loop (line 11) where $\theta_{h^{*}}\in\left[{\theta_{\ell}-\frac{\gamma}{4},\theta_{\ell}+\frac{\gamma}{4}}\right]$. More precisely, we have

[TABLE]

setting the constants in the definitions of $T_{0}$ and $\mu_{0}$ appropriately. In this iteration where $\theta_{h^{*}}\in\left[{\theta_{\ell}-\frac{\gamma}{4},\theta_{\ell}+\frac{\gamma}{4}}\right]$, $h^{*}$ satisfies the conditions of Lemma 4.11, so that $(\tilde{v},h^{*})$ gets added to $\mathcal{L}$ such that $\left|\tilde{v}\right|\geq(1-C\epsilon)\left\|{\hat{\mathbf{r}}}\right\|_{{\infty}}$. Thus with probability at least $1-\frac{\mu}{T_{2}}$, we have that $(\widetilde{v},h^{*})\in\mathcal{L}$ and Lemmas 4.11 and 4.12 have the favorable outcome if they are invoked, in every inner iteration during iteration of the outer loop.

Now by the definition of $\bar{v}$, we must have $\bar{v}\geq\widetilde{v}\geq(1-C\epsilon)\left\|{\hat{\mathbf{r}}}\right\|_{{\infty}}$. Then $|\hat{\mathbf{v}}[\bar{h}]|\geq\frac{1}{2}\left\|{\hat{\mathbf{r}}}\right\|_{{\infty}}$ or else Lemma 4.12 is contradicted. Then Lemma 4.11 implies (19), and combining with the previous equation establishes (20). ∎

By Claim 4.13 and a union bound, we have that with probability $1-\mu$ over the whole algorithm, (19) and (20) hold for every iteration of the outer loop (line 10). It remains to establish that assuming this is true, the conclusion of Lemma 4.2 holds.

We first note that if $\hat{\mathbf{z}}_{1},\hat{\mathbf{z}}_{2}$ are the values of $\hat{\mathbf{z}}$ before and after line 23 respectively (with corresponding $\hat{\mathbf{r}}_{i}=\hat{\mathbf{v}}-\hat{\mathbf{z}}_{i}$), then since $\mathbf{r}_{2}=\mathbf{r}_{1}-v\mathbf{e}_{\bar{h}}$ (19) says that

[TABLE]

In other words entries of $\hat{\mathbf{r}}$ only decrease in magnitude throughout the algorithm. Thus for all $i\in\mathrm{supp}(\hat{\mathbf{z}})$ at the beginning of the algorithm, (18) continue to hold throughout the algorithm.

Now we consider the two cases: either a $\bar{h}\not\in\mathrm{supp}(\hat{\mathbf{z}})$ is found in some iteration and Peeler returns $\texttt{stop}=\textnormal{\sc false}$, or one is never found and it returns $\texttt{stop}=\textnormal{\sc true}$. In the first case, we only need to establish that (18) holds for $\bar{h}$ as well at the time Peeler stops. It is easy to see that (19) and (20) imply (18) using the fact that $\tilde{\mathbf{z}}[\bar{h}]=\bar{v}$ and $\hat{\mathbf{r}}[\bar{h}]=\hat{\mathbf{v}}[\bar{h}]$. Further, by the assumption on $\left\|{\hat{\mathbf{w}}}\right\|_{{2}}$ implies that for small enough $\epsilon$, (18) implies that $\bar{h}\in\mathrm{supp}(\hat{\mathbf{x}})$.

Now consider the case that Peeler returns $\texttt{stop}=\textnormal{\sc true}$. That is, it completes the outer For loop and never chooses a pair $(\bar{v},\bar{h})$ so that $\bar{h}\not\in\mathrm{supp}(\hat{\mathbf{z}})$. In this case we will claim that $\hat{\mathbf{z}}$ was already very close to $\hat{\mathbf{v}}$ to begin with, in which case we will be done.

More precisely, given the choice of $T_{2}=\Theta(k\log_{1/\epsilon}(k/\epsilon))$ and the fact that $|\mathrm{supp}(\hat{\mathbf{z}})|\leq k$, if Peeler returns $\texttt{stop}=\textnormal{\sc true}$, then there is some $\bar{h}\in\mathrm{supp}(\hat{\mathbf{z}})$ that was chosen at least $\Theta(\log_{1/\epsilon}(k/\epsilon))$ times in the outer For loop. Let $\hat{\mathbf{z}}^{(0)},\hat{\mathbf{z}}^{(1)},\ldots,\hat{\mathbf{z}}^{(\log(k/\epsilon))}$ be the iterates of $\hat{\mathbf{z}}$ during all of the times that $\bar{h}$ was chosen (at the start of the outer loop); let $\hat{\mathbf{r}}^{(i)}$ be the corresponding residual vector $\hat{\mathbf{v}}-\hat{\mathbf{z}}^{(i)}$ and let $\bar{v}^{(i)}$ denote the corresponding value that was selected by Peeler. For each $i$, we have

[TABLE]

where in the second line we used (21) and in the third line we used (20). Thus, we have

[TABLE]

Iterating this $\Theta(\log_{1/\epsilon}(k/\epsilon))$ times, we see that the final $\tilde{\mathbf{z}}$ that is returned by Peeler along with $\texttt{stop}=\textnormal{\sc true}$ satisfies (assuming $\epsilon$ is sufficiently small)

[TABLE]

which is what we required from Peeler in the case that $\texttt{stop}=\textnormal{\sc true}$. Above, we are implicitly setting the constant inside the $\Theta$ in the definition of $T_{2}$ in (22). In the last inequality we used the fact that for $i\in\mathrm{supp}(\hat{\mathbf{z}})$, we have $\left|\hat{\mathbf{v}}[i]-\hat{\mathbf{z}}[i]\right|\leq C\epsilon\left|\hat{\mathbf{v}}[i]\right|\leq\left|\hat{\mathbf{v}}[i]\right|$ for small enough $\epsilon$.

This completes the proof of the error guarantee of Peeler.

4.3.3 Query and Time complexity of the Peeler algorithm

The Peeler algorithm calls the 1-sparse recovery algorithm $T_{0}\cdot T_{2}$ times, and for each of the $Q(N,6\epsilon,C\epsilon,\mu_{0}/2)$ queries that the 1-sparse recovery algorithm uses, SimulateQueryAccess uses $O(d)=O(\sqrt{k}/\epsilon\gamma)$ queries. Additionally, Verify is called $T_{0}\cdot T_{2}$ times, and uses SimulateQueryAccess $T_{1}$ times each time it is called. Thus, the total number of queries is

[TABLE]

The derivation of the running time is similar; the only overhead is that simulateQueryAccess has an additional factor of $k$ in its running time over its query complexity, leading to a running time of

[TABLE]

4.3.4 Correctness of the Verify algorithm

We begin by analyzing the estimates in Verify (which we will later use to prove Lemmas 4.11 and 4.12):

Lemma 4.14** (Estimator Lemma).**

Let $\zeta\leq 1$ be small enough. Let $\Omega$ be as chosen in Verify. Suppose that $\hat{\mathbf{u}}$ is of the form $\hat{\mathbf{u}}=\hat{\mathbf{a}}+\hat{\mathbf{q}}$, where $\left\|{\hat{\mathbf{q}}}\right\|_{{2}}\leq\zeta\left\|{\hat{\mathbf{a}}}\right\|_{{2}}$. Consider the estimate

[TABLE]

where the threshold $\tau$ satisfies ($U$ is as defined in (3))

[TABLE]

Then for any $\mu>0$, we have with probability at least $1-\mu$,

[TABLE]

for some choice of

[TABLE]

On the other hand for any $\hat{\mathbf{u}}\in\mathbb{R}^{N}$, $\tau,\zeta>0$, the estimate $\Theta$ above satisfies $\Theta\leq(1+\zeta)X^{2}$ with probability at least $1-\mu$ for any $X^{2}\geq\left\|{\hat{\mathbf{u}}}\right\|_{{2}}^{2}$ for some choice of

[TABLE]

We will use Bernstein’s inequality in the proof of above lemma, which we recall next:

Theorem 4.15** (Bernstein’s Inequality).**

Let $X_{1},\dots,X_{n}$ be independent random variables with $|X_{i}|\leq M$, then we have

[TABLE]

Proof of Lemma 4.14.

We start with the case when $\hat{\mathbf{u}}=\hat{\mathbf{a}}+\hat{\mathbf{q}}$, where $\left\|{\hat{\mathbf{q}}}\right\|_{{2}}\leq\zeta\left\|{\hat{\mathbf{a}}}\right\|_{{2}}$. Then note that

[TABLE]

where $\mathbf{q}=\mathbf{F}^{-1}\hat{\mathbf{q}}$. By definition of $U$, we have that

[TABLE]

For notational convenience, define $\mathbf{t}=\textsc{trunc}_{{\tau}}\left({{\mathbf{u}}}\right)$ and $\hat{\mathbf{t}}$ accordingly. Next, we claim that by (23), we have

[TABLE]

Indeed, we will argue that component-wise the elements of $\mathbf{z}$ are smaller (in absolute value) than $\mathbf{q}$. If $\left|\mathbf{u}[i]\right|\leq\tau$, then the claim is trivially true. If not, then $|\mathbf{u}[i]|>\tau$, $|\mathbf{t}[i]|=\tau$, and by (25) $|\mathbf{F}^{-1}\hat{\mathbf{a}}[i]|\leq\tau$. Since $\mathbf{u}[i]$ and $\mathbf{t}[i]$ have the same sign, $|\mathbf{z}[i]|=|\mathbf{t}[i]-\mathbf{F}^{-1}\hat{\mathbf{a}}[i]|\leq|\mathbf{u}[i]-\mathbf{F}^{-1}\hat{\mathbf{a}}[i]|=|\mathbf{q}[i]|$, as required.

Thus, we have $\mathbf{F}\cdot\textsc{trunc}_{{\tau}}\left({{\mathbf{u}}}\right)=\hat{\mathbf{a}}+\hat{\mathbf{z}}$ where $\left\|{\hat{\mathbf{z}}}\right\|_{{2}}\leq\left\|{\hat{\mathbf{q}}}\right\|_{{2}}$. Now by definition, we have

[TABLE]

To show that (24) holds, we first show that a similar bound holds for $\mathbf{u}$ replaced by $\mathbf{t}$ (via Bernstein’s inequality) and then show that $\left\|{\mathbf{u}}\right\|_{{2}}\approx\left\|{\mathbf{t}}\right\|_{{2}}$. Recall that the set $\Omega$ in Verify is chosen by including $T_{1}$ indices of $[N]$ uniformly at random with replacement. Suppose these indices are $j_{1},\ldots,j_{T_{1}}$, and let $X_{i}=\mathbf{t}[j_{i}]^{2}$. Thus,

[TABLE]

Further, for all $i\in[T_{1}]$, we have $|X_{i}|\leq\|\mathbf{t}\|_{\infty}^{2}$ and $\mathbb{E}X_{i}^{2}=\frac{1}{N}\sum_{i=1}^{N}\mathbf{t}[i]^{4}$. Thus by Bernstein’s inequality with $t=\frac{T_{1}}{N}\cdot\zeta\left\|{\mathbf{t}}\right\|_{{2}}^{2}$, we have

[TABLE]

Multiplying both sides of the inequality inside the probability expression by $\frac{N}{T_{1}}$ and replacing each sum by the appropriate norm expressions, we have

[TABLE]

In the above (27) follows from the fact that $\left\|{\mathbf{t}}\right\|_{{4}}\leq\sqrt{\left\|{\mathbf{t}}\right\|_{{\infty}}\cdot\left\|{\mathbf{t}}\right\|_{{2}}}$, and (28) follows from (26). Finally, (29) follows by choosing

[TABLE]

Thus, we have argued that with probability at least $1-\mu$, we have

[TABLE]

However, we wanted to prove a similar result with $\mathbf{t}$ replaced by $\mathbf{u}$. We do so next by showing that $\left\|{\mathbf{u}}\right\|_{{2}}\approx\left\|{\mathbf{t}}\right\|_{{2}}$. Indeed consider the following sequence of inequalities:

[TABLE]

where the first inequality follows from the fact that $\left\|{\mathbf{z}}\right\|_{{2}}\leq\left\|{\mathbf{q}}\right\|_{{2}}$, the second and the last inequality follow from the triangle inequality and the third inequality follows from the definition of the truncation function. Applying the bound $\left\|{\mathbf{q}}\right\|_{{2}}\leq\zeta\left\|{\hat{\mathbf{a}}}\right\|_{{2}}$, we get that

[TABLE]

In other words for $\zeta\leq 1$,

[TABLE]

Applying the above in (31), we get that with probability at least $1-\mu$

[TABLE]

Further, using the fact that $\zeta\leq 1$, implies that the RHS is upper bounded by $9\zeta\left\|{\mathbf{u}}\right\|_{{2}}^{2}$, as desired. Finally, the bound $\left\|{\mathbf{t}}\right\|_{{2}}\leq(1+\zeta)\left\|{\hat{\mathbf{a}}}\right\|_{{2}}$ along with (30) implies the required bound on $T_{1}$ in the lemma statement.

For the second part, we consider the same random variables $X_{i}$ as defined above but we use $t=\frac{T_{1}}{N}\cdot\zeta X^{2}$. Then by applying Bernstein’s inequality we get that

[TABLE]

Noting that $\left\|{\mathbf{t}}\right\|_{{2}}^{2}\leq\left\|{\mathbf{u}}\right\|_{{2}}^{2}\leq X^{2}$, the above (along with (26)) implies that

[TABLE]

where the last inequality follows from the choice of $T_{1}$ in the second part of the lemma. The proof is complete. ∎

4.3.5 Proof of Lemma 4.11

In this section, we use Lemma 4.14 to prove 4.11.

First note that $\left\|{\mathbf{w}}\right\|_{{2}}\leq\epsilon\left\|{\hat{\mathbf{v}}-\hat{\mathbf{z}}}\right\|_{{\infty}}<\left|v\right|$ for $\epsilon<\frac{1}{2}$. Hence we can assume that $h\in\mathrm{supp}(\hat{\mathbf{x}})$. Now suppose $\theta_{h}\in\left[{\theta_{\ell}-\frac{\gamma}{4},\theta_{\ell}+\frac{\gamma}{4}}\right]$. Then by part (2) of Lemma 4.5, we have that no other $\theta_{h^{\prime}}$ for $h^{\prime}\in\mathrm{supp}(\hat{\mathbf{x}})$ falls in the range

[TABLE]

This implies that we have $\hat{\mathbf{y}}_{\ell}$ (as defined in Peeler) can be expressed as

[TABLE]

where $\hat{\mathbf{q}}$ in the range $R$ only has contribution from $\hat{\mathbf{w}}$ (multiplied by a factor of at most $1+\epsilon/\sqrt{k}$) and outside the range $R$ it has contribution from $(\hat{\mathbf{v}}-\hat{\mathbf{z}})$ (multiplied by a factor of $\frac{\epsilon}{\sqrt{k}}$). (See Figure 3.)

Recalling that $\hat{\mathbf{r}}=\hat{\mathbf{v}}-\hat{\mathbf{z}}$, we can bound the noise $\hat{\mathbf{q}}$ by

[TABLE]

In the above the second line uses the properties of the boxcar polynomial $\mathbf{b}_{\ell}$ and the fact that

[TABLE]

by the fact that $\hat{\mathbf{x}}_{R\setminus\{h\}}=\mathbf{0}$ and the fact that $\mathrm{supp}(\hat{\mathbf{z}})\subseteq\mathrm{supp}(\hat{\mathbf{x}})$ (which in turn is always maintained in each application of Lemma 4.2). The third line follows from the fact that we can write

[TABLE]

The fourth line follows from the fact that $|\mathrm{supp}(\hat{\mathbf{x}})|\leq k$, and the last line uses the assumptions on $\left\|{\mathbf{w}}\right\|_{{2}}$ in the lemma statement. Finally the assumption on $\left|v\right|$ in the lemma statement yields

[TABLE]

The above implies that the ‘error’ is small enough to run our $\left({N,6\epsilon,C\epsilon,\frac{\mu_{0}}{2}}\right)$ $1$-sparse recovery algorithm. Hence, with probability at least $1-\mu_{0}/2$, we get an estimate $(y_{\ell},h)$ from the $1$-sparse recovery algorithm such that

[TABLE]

Now we apply the second part888Note that $\tau=100\left|v\right|U\geq\left|v\right|U$ and hence it satisfies (23). of Lemma 4.14 with $\hat{\mathbf{u}}=\hat{\mathbf{y}}-y_{\ell}\mathbf{e}_{h}$, $X=\frac{\left|y_{\ell}\right|}{40}$ and999We need $\left\|{\hat{\mathbf{y}}-y_{\ell}\mathbf{e}_{h}}\right\|_{{2}}\leq(6+C)\epsilon\left|v\right|\leq X$. Note that $(6+C)\epsilon\left|v\right|\leq\frac{(6+C)\epsilon\left|y_{\ell}\right|}{(1-C\epsilon)}$ and then the required bound follows for a sufficiently small choice of $\epsilon$. $\zeta=\frac{3}{5}$ to note that the estimate $g$ in Verify is, with probability at least $1-\frac{\mu_{0}}{2}$:

[TABLE]

as long as

[TABLE]

which along with the fact that we picked $\tau=100\left|y_{\ell}\right|U$, (32) implies that we would be fine if we picked:

[TABLE]

which we did.

Finally, note that by the union bound with probability at least $1-\mu_{0}$, we have that both (32) and (33) are satisfied. In other words, the 1-sparse recovery algorithm succeeds and Verify returns true, which implies that $(y_{\ell},h)$ gets added to $\mathcal{L}$. This along with (32) completes the proof.

4.3.6 Proof of Lemma 4.12

Suppose that in a given iteration, the $1$-sparse recovery solver returns $(y_{\ell},h)$. Recall that we have $\theta_{h}\in\left[{\theta_{\ell}-\frac{\gamma}{4},\theta_{\ell}+\frac{\gamma}{4}}\right]$ as otherwise $h$ will not get added to $\mathcal{L}$. If $\left|y_{\ell}\right|<\left({1-C\epsilon}\right)\left\|{\hat{\mathbf{v}}-\hat{\mathbf{z}}}\right\|_{{\infty}}$, then we have nothing to prove. So for the rest of the proof, let us assume

[TABLE]

Let us look at the potential spikes captured besides $h$,

[TABLE]

Note that by part (2) of Lemma 4.5 there can only be one such $h^{\prime}$, so this set has the form either $\{h^{\prime}\}$ for some $h^{\prime}\neq h$, or it is empty in which case for shorthand we set $h^{\prime}=\bot$.

Let $\hat{\mathbf{y}}_{\ell}$ be as defined in Peeler algorithm. Now the vector $\hat{\mathbf{y}}_{\ell}$ has two potential locations of interest at $h$ and $h^{\prime}$ and the rest can be expressed as a vector $\hat{\mathbf{q}}$

[TABLE]

where $v=\hat{\mathbf{r}}[h]$ and $v^{\prime}=\hat{\mathbf{r}}[h^{\prime}]$ Note in the case $h^{\prime}=\bot$ the second term disappears and this says $\hat{\mathbf{y}}_{\ell}=\hat{\mathbf{r}}[h]\cdot\mathbf{e}_{h}+\hat{\mathbf{q}}$, the same as in Lemma 4.11. (See Figure 4.)

First, we note that the same derivation as in proof of Lemma 4.11 works now to get an upper bound on $\|\hat{\mathbf{q}}\|_{2}$:

[TABLE]

Next we prove a lower bound on the magnitude of the vector that Verify estimates:

[TABLE]

The first line follows from the triangle inequality. The second line follows from our assumption $|\hat{\mathbf{r}}[h]|\leq\frac{1}{2}\|\hat{\mathbf{r}}\|_{\infty}$ and from (35).

Now we apply the first part of Lemma 4.14 with $\hat{\mathbf{u}}=\hat{\mathbf{y}}_{\ell}-y_{\ell}\mathbf{e}_{h}=\left[\left({\hat{\mathbf{r}}[h]-y_{\ell}}\right)\mathbf{e}_{h}+\hat{\mathbf{r}}[h^{\prime}]\mathbf{e}_{h^{\prime}}\right]+\hat{\mathbf{q}}$ and $\zeta=24\epsilon$. First, we verify that for small enough $\epsilon$

[TABLE]

and

[TABLE]

Thus, Lemma 4.14 along with (36) implies that with probability at least $1-\mu_{0}$,

[TABLE]

as long as

[TABLE]

which is how we picked $T_{1}$.

This implies that the check in Verify returns false and hence $h$ will not get added to $\mathcal{L}$.

5 One-Sparse Recovery Algorithm for Jacobi Polynomials

In this section we will present a one-sparse recovery algorithm for Jacobi polynomials [szego]. In particular, we consider a fixed family of Jacobi polynomials with parameters $\alpha,\beta>-1$ and the corresponding Jacobi transform $\mathbf{F}$ through Definition 2.1. We will prove:

Theorem 5.1**.**

There exists a universal constant $C$ such that the following holds. Fix any $\alpha,\beta>-1$. There exists an $(N,\epsilon,C\cdot\epsilon,\mu)$ $1$-sparse recovery algorithm for the Jacobi transform (with parameters $\alpha,\beta$) that makes

[TABLE]

queries and takes time

[TABLE]

In the above the Big-Oh notation hides constants that depend on $\alpha$ and $\beta$.

Recalling Definition 2.4, in other words we seek a solution to the following specific problem:

• •

Unknown: Index $\ell\in[N]$ and magnitude $v$.

• •

Query access: To a vector

[TABLE]

for some $\|\mathbf{w}\|_{2}\leq\epsilon|v|$.

• •

Output: With probability $1-\mu$, return $\ell$ and some $\widetilde{v}$ such that $|\tilde{v}-v|\leq C\epsilon|v|$.

To do so, we first provide the required background for Jacobi polynomials, including important bounds from [szego] (some of which are deferred to Appendix D). Then, we give an overview of the entire algorithm with a comprehensive summary of the notation used. There are two basic building blocks we need: (i) $\cos(\cdot)$ evaluation from Jacobi polynomial evaluations and (ii) an approximate $\arccos$ algorithm from the noisy $\cos(\cdot)$ evaluation (whose proof we cover in Section 7). We present the basic building blocks first and then delve into the details of the entire algorithm.

5.1 Background: Jacobi polynomials

Jacobi polynomials are indexed by two parameters $\alpha,\beta>-1$ and are orthogonal with respect to the measure

[TABLE]

in the range $[-1,1]$.

Traditionally, Jacobi polynomials $P_{j}^{(\alpha,\beta)}$ are defined by the following recurrence relation [szego, Sec. 4.5, (4.5.1)] for $j\geq 2$:101010$P_{0}^{(\alpha,\beta)}(X)=1$ and $P_{1}^{(\alpha,\beta)}(X)=\frac{\alpha+\beta+2}{2}\cdot X+\frac{\alpha-\beta}{2}$.

[TABLE]

However, the polynomials $P_{j}^{(\alpha,\beta)}(X)$ while being orthogonal w.r.t. the weight $w^{(\alpha,\beta)}(X)$, are not orthonormal w.r.t. it since for any $j$, we have ( [szego, Pg. 68, (4.3.3)]):

[TABLE]

where

[TABLE]

with $\Gamma(z)$ being the Gamma function.

To make the rest of the exposition simpler, we will define the Jacobi polynomial of degree $j$ with parameters $\alpha,\beta$ as

[TABLE]

Note that these polynomials are indeed orthonormal with respect to the measure $w^{(\alpha,\beta)}(X)$.

Chebyshev polynomials are special case of $\alpha=\beta=-\frac{1}{2}$ and Legendre polynomials are the special case of $\alpha=\beta=0$ (up to potentially a multiplicative factor that could depend on $j$).

We start with the following property of roots of Jacobi polynomials:

Theorem 5.2** ([szego], Theorem 8.9.1).**

Let

[TABLE]

Then for any integer $N$, let the roots of $\mathcal{J}^{({\alpha},{\beta})}_{{N}}(X)$ be given by $\cos{\theta_{\ell}}$ for $\ell\in[N]$. Then, there is a universal constant $C^{\alpha,\beta}$ such that

[TABLE]

Theorem 5.2 establishes the denseness (Definition 2.2) of Jacobi families, made explicit in Corollary 6.1. Additionally, we will make use of the following immediate corollary:

Corollary 5.3**.**

Let $[a,b]\subseteq[0,\pi]$ be an interval. Then the number of roots of $\mathcal{J}^{({\alpha},{\beta})}_{{N}}(X)$ for the form $\cos\left({\theta_{\ell}}\right)$ such that $\theta_{\ell}\in[a,b]$ is upper bounded by $O\left({\frac{(b-a)N}{\pi}}\right)$. Further, the list of such roots can be computed in time $O\left({(b-a)N}\right)$ assuming that there is a data structure that for any $i\in[N]$, stores all the $O(C^{\alpha,\beta})$ roots in the interval $\left[\frac{i\pi}{N},\frac{(i+1)\pi}{N}\right)$.

Further, we will also need the fact that all the $\theta_{\ell}$ values are bounded away from [math] and $\pi$:

Lemma 5.4**.**

Let $\alpha,\beta>-1$. Then there exists a constant $C^{\prime}>0$ such that the following holds. Let $N\geq 1$ be large enough. Let $\lambda_{0}\leq\cdots\leq\lambda_{N-1}$ be the roots of the $N$th Jacobi polynomial. Then for every $0\leq\ell<N$,

[TABLE]

Proof.

Theorem 8.1.2 in [szego] states that for large enough $N$, we have

[TABLE]

where $j^{\alpha}_{0}$ is the first positive root of $J_{{\alpha}}\left({{X}}\right)$. It is known that (see e.g. equation (5.3) in [bessel-root-lb]):

[TABLE]

Further, by Lemma D.5, we have that

[TABLE]

Thus, the lemma follows if we pick

[TABLE]

which is strictly positive since $\alpha,\beta>-1$. The proof is complete. ∎

Next, we will need an asymptotic bound on values of Jacobi polynomials111111The bounds are generally stated for $P_{j}^{(\alpha,\beta)}(X)$– we have updated it for $\mathcal{J}^{({\alpha},{\beta})}_{{j}}(X)$ by adding the required factor of $\frac{1}{\sqrt{h_{j}^{\alpha,\beta}}}$ to the approximation.:

Theorem 5.5** ([szego], Theorem 8.21.13).**

Let $\alpha,\beta>-1$. There is a constant $c>0$ such that the following holds. For any integer $j\geq 1$, let

[TABLE]

Then, there exists a constant $c^{\prime}$ such that

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Remark 5.6**.**

When $\alpha$ and $\beta$ are clear from context (which will be the case for the rest of this section), we will just use $\kappa(\theta)$ instead of $\kappa_{\alpha,\beta}\left({\theta}\right)$.

We will need the following bound on the normalization factor $h_{j}^{\alpha,\beta}$:

Lemma 5.7**.**

For large enough $j$, we have

[TABLE]

where $\Theta_{x,y}(\cdot)$ hides factors that depends on $x$ and $y$.

Proof.

We will use Stirling’s approximation for the Gamma function: for $z\geq 1$ we have

[TABLE]

By definition of $h_{j}^{\alpha,\beta}$, we have

[TABLE]

Note that for large enough $j$, the first fraction above is $\Theta_{\alpha,\beta}(1)$ so to complete the proof we will argue that the second fraction is also $\Theta_{\left|\alpha\right|,\left|\beta\right|}(1)$. Indeed, we will apply Stirling’s approximation to the Gamma function (where for notational convenience we will ignore the constant factors in the Stirling’s approximation) and get

[TABLE]

The proof is complete by noting that each of the remaining terms is $e^{O(\left|\alpha\right|+\left|\beta\right|)}$ for large enough $j$, which is $\Theta_{\left|\alpha\right|,\left|\beta\right|}(1)$, as desired, ∎

Finally, we will need the following result in our analysis.

Theorem 5.8**.**

Let $\alpha,\beta>-1$. Let $N\geq 1$ be large enough. Let $\lambda_{0},\dots,\lambda_{N-1}$ be the roots of the $N$th Jacobi polynomial. Then define $\mathbf{F}$ such that for every $0<\ell,j<N$, we have

[TABLE]

where

[TABLE]

Then for every $0\leq\ell<N$, we have

[TABLE]

and

[TABLE]

It turns out that the above result does not quite follow from Theorem 5.5 and it needs some extra work to prove (45) for the entries where equation (42) is not satisfied. In particular, we need a more general version of Theorem 5.5 (which involves the Bessel function). The proof of Theorem 5.8 is deferred to Appendix D.3 (with the more general version of Theorem 5.5 proven in Appendix D.1).

5.2 Overview of the algorithm

The main idea will be to use Theorem 5.11 to identify $\theta_{\ell}$ accurately enough (Section 5.5.3), which identifies $\ell$. Given $\ell$, $v$ can be estimated with a simple sampling argument (Algorithm 6).

In order to estimate $\theta$ accurately enough, Theorem 5.11 will require querying $\cos(w\theta_{\ell})$ for $w$ up to $\Omega(N)$, say $w\leq\nu N$ for some (universal, to be chosen later) constant $\nu$.

These queries in turn will be provided by Corollary 5.9, which computes $\cos(w\theta)$ from queries $\mathcal{J}^{({\alpha},{\beta})}_{{j}}(\cos\theta)$ for $j=\Delta-w,\Delta,\Delta+w$ for certain $\Delta$ that satisfy a technical condition (equation (46)).

• •

In order for $\mathcal{J}^{({\alpha},{\beta})}_{{j}}(\cdot)$ to be meaningful (related to $\cos$ values) via Theorem 5.5, either $j$ must be sufficiently large enough or alternatively $\theta_{\ell}$ cannot be too close to the boundaries of $[0,\pi]$ (equation (42)). By restricting the queries to $j\geq\nu N$ (enforced by only querying $\Delta\in[2\nu N,N-2\nu N]$), there will be a constant number of possible bad values of $\ell$.

• •

Because we can’t determine if (46) is true a priori, we will instead query for many random values of $\Delta$ in the allowed range. For most values of $\ell$, there are enough good values of $\Delta$ (Definition 5.15). There will be a small number of bad values of $\ell$, located in Lemma 5.16.

Finally, Algorithm 6 (with bounds given in Corollary 5.14) allows us to check if $\ell$ is any specific value. Thus the final algorithm first checks if $\ell$ is any of the possible bad values; if not, it narrows down $\ell$ using Theorem 5.11 and Corollary 5.9. This gives a final window for $\theta_{\ell}$ of size $O(1/\nu N)$, with $O(1/\nu)$ possible values of $\ell$ that can be checked with Algorithm 6.

Summary of notation

In Table 1 we summarize notation for frequently used constants and parameters in this section, including a brief description and where they are defined and used.

5.3 $\cos(\cdot)$ values from Jacobi polynomial evaluation

We use Theorem 5.5 to get the following lemma, which provides us with noisy estimates of $\cos(w\theta)$ from querying values of $\mathcal{J}^{({\alpha},{\beta})}_{{\square}}(\cos\theta)$:

Corollary 5.9**.**

Consider a fixed $w\in[N]$ and $\theta\in[0,\pi]$. Then for any integer $\Delta>w$ such that $j=\Delta-w$ satisfies (42),

[TABLE]

and we have access to values $\widetilde{{\mathcal{J}^{({\alpha},{\beta})}_{{j}}}}(\cos{\theta})$ for $j\in\left\{\Delta-w,\Delta,\Delta+w\right\}$ such that

[TABLE]

we have

[TABLE]

Proof.

This follows by applying Theorem 5.5 and the following elementary identity in a straightforward way:

[TABLE]

For notational convenience define

[TABLE]

and

[TABLE]

Now consider the following relations:

[TABLE]

as desired. In the above the first equation follows from Lemma 5.7, (5.3) follows from Theorem 5.5, (51) follows from (46), and (52) follows from (49). (Note that $\cos{A}\neq 0$ by (46).)

∎

5.4 Computing $\theta$ from noisy $\cos(\cdot)$ evaluation

We being with a definition:

Definition 5.10**.**

Let $\rho_{\epsilon_{0}}$ be such that for any $\theta\in[0,\pi]$,

[TABLE]

In Section 7, we will prove the following result:

Theorem 5.11**.**

There is an algorithm ApproxArcCos such that the following holds. Let $\theta\in(0,\pi)$. Then for any integer $\tau\geq 1$ and small enough $\epsilon_{0}$ so that $0<\rho_{\epsilon_{0}}<\pi/22$, given access to evaluations

[TABLE]

where $w\in\left\{x_{t}\cdot 2^{t-1}|\text{ for every }1\leq t\leq\tau\text{ with }x_{t}\in[1,3/2]\right\}$, then $\textsc{ApproxArcCos}(\tau,\epsilon_{0})$ (for $\tau\leq\left\lfloor\log_{2}(2N/3)\right\rfloor$) returns a range $S_{\tau}\subseteq[0,\pi]$ such that

• •

$\theta\in S_{\tau}$; and

• •

$\left|S_{\tau}\right|\leq\frac{2\rho_{\epsilon_{0}}}{2^{\tau}}$.

Finally, $\textsc{ApproxArcCos}(\tau,\epsilon_{0})$ runs in $O(\tau)$ time and makes $O(\tau)$ queries.

We would like to point out that we do not handle $\theta\in\left\{0,\pi\right\}$ in the result above. Our proof of the result above is a bit cleaner with this assumption and since we will only use the above result for $\theta=\theta_{\ell}$ where $\cos\theta_{\ell}$ is a root of the $N$th degree Jacobi polynomial this is not an issue (due to Lemma 5.4).

5.5 Details of the one-sparse recovery algorithm for Jacobi polynomials

5.5.1 A pruning step

Before we present the final $1$-sparse recovery algorithm, we need one final result that allows us to check if $\ell$ lies in a given set of candidate values:

Lemma 5.12**.**

Let $0<\epsilon<\frac{1}{200}$ be small enough. There exists an algorithm that given a subset $S\subseteq[N]$ and query access to $\mathbf{y}=v\cdot\mathbf{F}^{-1}\mathbf{e}_{\ell}+\mathbf{w}$ (such that $\left\|{\mathbf{w}}\right\|_{{2}}\leq\epsilon\left|v\right|$) does the following with probability at least $1-\mu$: if $\ell\in S$, it outputs $\ell$ and an estimate $\tilde{v}$ such that $\left|\tilde{v}-v\right|\leq 13\epsilon\left|v\right|$; otherwise it outputs fail. Further it makes $q_{\textsc{Prune}}\left({{|S|,N,\mu,\epsilon}}\right)$ queries and takes $T_{\textsc{Prune}}\left({{|S|,N,\mu,\epsilon}}\right)$ time, where

[TABLE]

and

[TABLE]

Proof.

We will solve the problem for $|S|=1$ and then just repeat the algorithm for each $\ell^{\prime}\in S$. So as long as we can solve the problem for $|S|=1$ with probability at least $1-\frac{\mu}{\left|S\right|}$, we would be fine.

Thus, the problem we want to solve is that given a target $\ell^{\prime}$, we have to decide if $\ell^{\prime}=\ell$ or not (with query access to $\mathbf{y}$). Algorithm 6 has the details on how to solve this.

Next we argue that Algorithm 6 works as it is supposed to:

Lemma 5.13**.**

Let $(b,\tilde{v})$ be the output of $\textsc{Check}^{(\mathbf{y})}\left({\ell^{\prime},\mu,\epsilon}\right)$. Then with probability at least $1-\mu$ the following is true for small enough $\epsilon$. If $\ell=\ell^{\prime}$ then $b=\textnormal{\sc true}$ and $\left|\tilde{v}-v\right|\leq 13\epsilon\cdot\left|v\right|$. Otherwise, $b=\textnormal{\sc false}$. Further, it makes $O\left({U^{2}N\cdot\frac{\log{N}}{\epsilon^{2}}\cdot\log\left(\frac{\log{N}}{\epsilon^{2}}\right)\log\left({\frac{1}{\mu}}\right)}\right)$ queries and runs in time linear in the number of queries.

Proof.

Fix any $i\in[R]$ and let $\Gamma$ be the subset of indices chosen in the $i$th iteration. We will argue that each of the following two events hold with probability at least $\frac{4}{5}$:

Let $\mathbf{w}_{\Gamma}$ be the vector $\mathbf{w}$ projected to coordinates in $\Gamma$. Then

[TABLE] 2. 2.

The following is true for any $h,h^{\prime}\in[N]$ and $v\in\mathbb{R}$:

[TABLE]

Assuming the above are true, we first complete the proof. First let us consider $u_{i}$:

[TABLE]

In the above the first containment follows from the definition of $\mathbf{y}$ and the triangle inequality while the second containment follows from (55) with $h=h^{\prime}=\ell$ and (54). Next we consider $\tilde{v}_{i}$:

[TABLE]

In the above the equality follows from the definition of $\tilde{v}_{i}$ and $\mathbf{y}$. The first containment follows from (55) (with $h=\ell$ and $h^{\prime}=\ell^{\prime}$) and Cauchy-Schwartz while the second containment follows from (54) and (55) (with $v=1$ in the latter).

Thus, by the union bound we have that with probability at least $\frac{3}{5}$ both (56) and (57) hold. Then an application of Chernoff bound implies that with probability at least $1-\mu$, we have

[TABLE]

Note that if $\ell=\ell^{\prime}$ then $\tilde{v}$ does have the required estimate. So we just need to argue that Check outputs true if $\ell=\ell^{\prime}$ and false otherwise. To see this note that the above bounds on $u$ and $\tilde{v}$ imply that

[TABLE]

Now when $\ell\neq\ell^{\prime}$, then the ratio above is at most $\frac{13\epsilon}{1-13\epsilon}\leq\frac{1}{10}$ by our choice of $\epsilon\leq\frac{1}{200}$. On the other hand, when $\ell=\ell^{\prime}$, the above ratio is at least $\frac{1-13\epsilon}{1+13\epsilon}>\frac{1}{10}$ by our choice of $\epsilon$. Thus, Check returns true or false based on if $\ell=\ell^{\prime}$ or not, as desired.

To complete the proof we argue (54) and (55). To begin with (54). Note that by definition we have

[TABLE]

Then Markov’s inequality implies that with probability at least $\frac{4}{5}$,

[TABLE]

where the second inequality follows from the assumed upper bound on $\left\|{\mathbf{w}}\right\|_{{2}}$. This proves (54).

Finally, we tackle (55): this basically follows by noting that $\mathbf{B}$ (as defined in Corollary C.2) is exactly $\sqrt{\frac{N}{s}}\mathbf{F}^{T}[\Gamma,:]$ (with $\mathbf{A}=\mathbf{F}^{T}$). Finally, note that (55) is (68) (scaled by a factor of $v$). This completes the proof. ∎

Finally, we use the Check algorithm repeatedly to get our final claimed result. In particular, for every $\ell^{\prime}\in S$, we run $\textsc{Check}^{(\mathbf{y})}\left({\ell^{\prime},\frac{\mu}{\left|S\right|},\epsilon}\right)$. If the call for $\ell^{\prime}$ returns $\left({\textnormal{\sc true},\tilde{v}}\right)$, then we return $\ell^{\prime}$ and $\tilde{v}$. Otherwise, we return fail. The correctness of this algorithm and the time complexity follows from Lemma 5.13 (and the union bound). Note that they query complexity does not have the $\left|S\right|$ factor it– this is because Corollary C.2 holds for all pairs of indices $(\ell^{\prime},\ell)$ simultaneously. Hence we can re-use the same $\Gamma$ across all $\ell^{\prime}\in S$ and this proves the query complexity, as desired. ∎

In fact, we will use Corollary 5.3 along with the above to obtain the following corollary, which will be useful in our final algorithm:

Corollary 5.14**.**

There is an algorithm Prune with the following property. If Prune has query access to $\mathbf{y}=v\cdot\mathbf{F}^{-1}\mathbf{e}_{\ell}+\mathbf{w}$ (such that $\left\|{\mathbf{w}}\right\|_{{2}}\leq\epsilon\left|v\right|$), and $\theta_{\ell}\in[a,b]$ where $[a,b]\subseteq[0,\pi]$, then $\textsc{Prune}^{(\mathbf{v})}(a,b,\mu,\epsilon)$ outputs $\ell$ with probability at least $1-\mu$ with $q_{\textsc{Prune}}\left({{\frac{(b-a)N}{\pi},N,\mu,\epsilon}}\right)$ queries and $T_{\textsc{Prune}}\left({{\frac{(b-a)N}{\pi},N,\mu,\epsilon}}\right)$. If $\theta_{\ell}\not\in[a,b]$, then with probability at least $1-\mu$, $\textsc{Prune}^{(\mathbf{v})}(a,b,\mu,\epsilon)$ outputs fail.

5.5.2 Bad potential values of $\ell$

We will also use Prune to ‘prune’ out some other ‘bad’ potential value of $\ell$, namely those for which Corollary 5.9 is difficult to apply because not enough values of $\Delta$ satisfy (46). For this we define the notion of ‘bad’ value:

Definition 5.15**.**

For any $0\leq p\leq 1$, $0\leq\nu\leq\frac{1}{2}$ and $0\leq\delta\leq 1$, we say an index $\ell\in[N]$ is $(p,\nu,\delta)$-spread if for a random $\Delta\in[\nu N,N-\nu N]$ it is the case that with probability at least $1-p$, we have:

[TABLE]

We will show that there are not that many indices that are not spread:

Lemma 5.16**.**

The number of indices that are not $\left({\frac{4\rho_{\delta}}{\pi(1-2\nu)},\nu,\delta}\right)$-spread (for any $\nu<1/2$) is $O\left({\left({\frac{1}{\rho_{\delta}^{2}}}\right)}\right)$.

Proof.

In this proof we will let $x\mod{y}$ for reals $x,y$ denote the unique value $x^{\prime}$ in $[0,y)$ such that $(x-x^{\prime})/y\in\mathbb{Z}$. For example, $x\mod{1}$ denotes its fractional part $x-\left\lfloor x\right\rfloor$.

Let $0\leq\tilde{\ell}<N$ be a real number such that

[TABLE]

Recall we want to avoid the case for $\Delta\in[\nu N,N-\nu N]$ that we have

[TABLE]

By Lemma 7.1 it is enough to avoid the case that

[TABLE]

for some $z\in\mathbb{Z}$, which is the same as

[TABLE]

Note that $\ell$ is fixed while we vary $\Delta$, which means for a fixed $\ell$, we want to figure out the probability (for a random $\Delta\in[\nu N,N-\nu N]$) that $\Delta\frac{\tilde{\ell}\pi}{N}\mod{\pi}$ is in a range $[L,R]$ of size $2\rho_{\delta}$. Note here we have $L=-\left({\frac{\alpha+\beta+1}{2}}\right)\theta_{\ell}+\left({\frac{\alpha}{2}+\frac{1}{4}}\right)\pi+\frac{\pi}{2}-\rho_{\delta}$ and $R=-\left({\frac{\alpha+\beta+1}{2}}\right)\theta_{\ell}+\left({\frac{\alpha}{2}+\frac{1}{4}}\right)\pi+\frac{\pi}{2}+\rho_{\delta}$. However, for the rest of the argument the only thing we will use about $L$ and $R$ is that $R-L=2\rho_{\delta}$.

We will bound the largest probability that for a random $\Delta\in[\nu N,N-\nu N]$, the angle $\Delta\cdot\frac{\tilde{\ell}\pi}{N}\mod{\pi}$ falls into any range inside $[0,\pi]$ of size $2\rho_{\delta}$. Further since the $\cos$ value only changes sign when then angle is increased by an integer multiple of $\pi$, we need to upper bound the size of the set:

[TABLE]

We bound the size of the set above by bounding:

[TABLE]

Or more precisely we want to figure out for how many roots $\cos\left({\theta_{\ell}}\right)$, the size of the set above is at most $2\cdot\frac{R-L}{\pi}\cdot N$. We solve this problem in Appendix E. As per the terminology in there, we consider the sequence $\frac{\tilde{\ell}}{N}$ (where $\ell\in\mathbb{Z}_{N}$ indexes the roots) and it suffices to bound the number of $\left({\frac{2\rho_{\delta}}{\pi}}\right)$-bad elements (as defined in Definition E.2) in this sequence. By Theorem 5.2, this sequence is $O(1)$-scattered (as per Definition E.1). Then Lemma E.3 implies that the number of $\tilde{\ell}$ such that the size in (58) is at least $\left({\frac{4\rho_{\delta}}{\pi}}\right)N$ is at most $O\left({\frac{1}{\rho_{\delta}^{2}}}\right)$ (call such an element bad). Note that if $\tilde{\ell}$ is not bad, then for a random $\Delta\in[\nu N,N-\nu N]$, the probability that $\Delta\frac{\tilde{\ell}\pi}{N}\in[L,R]$ is at most $\frac{1}{1-2\nu}\cdot\frac{4\rho_{\delta}}{\pi}$. This completes the proof. ∎

The above along with Lemma 5.12 (and Corollary E.8) implies the following result:

Corollary 5.17**.**

There exists an algorithm $\textsc{PruneNonSpread}(\delta,\nu,\mu,\epsilon)$ that given query access to $\mathbf{y}=v\cdot\mathbf{F}^{-1}\mathbf{e}_{\ell}+\mathbf{w}$ (such that $\left\|{\mathbf{w}}\right\|_{{2}}\leq\epsilon\left|v\right|$), with probability $1-O\left({\frac{\mu}{\rho_{\delta}^{2}}}\right)$ outputs $\ell$ if $\ell$ is not $\left({\frac{4\rho_{\delta}}{\pi(1-2\nu)},\nu,\delta}\right)$-spread. Otherwise it outputs fail. In either case, it uses $O\left({\frac{1}{\rho_{\delta}^{2}}\cdot q_{\textsc{Prune}}\left({{O(1),N,\mu,\epsilon}}\right)}\right)$ queries and $O\left({\frac{1}{\rho_{\delta}^{2}}\cdot T_{\textsc{Prune}}\left({{O(1),N,\mu,\epsilon}}\right)}\right)$ time.

5.5.3 Computing $\ell$ and $v$

The overall algorithm to compute $\ell$ and $v$ is presented in Algorithm 7.

Algorithm 7 needs access to Algorithm 8, which is based on Corollary 5.9.

Proof of correctness of Algorithm 7.

We first argue that Assertion 1-6 in Algorithm 7 and 8 hold with high (enough) probability.

Lemma 5.18**.**

With probability at least $1-\frac{\mu}{6}$, Assertion 1 holds.

Proof.

We need to have the following

[TABLE]

In other words we need (for large enough $N$),

[TABLE]

where in the above we have used that fact that $\sin(x)\approx x$ for small enough $x$. But the above is handled by the first call to Prune by making sure $C^{\prime}$ is large enough. Assertion 1 fails only if the call to Prune fails, which happens with probability at most $\frac{\mu}{6}$ (as per Corollary 5.14).

In order to further ensure that $\theta_{\ell}\geq\frac{c}{\nu N}$, it suffices to also require that $C^{\prime}\geq c$ (and hence $C^{\prime}\geq c\sqrt{\epsilon}\delta_{0}$). ∎

Lemma 5.19**.**

Conditioned on Assertion 1 being true, with probability at least $1-\frac{\mu}{6}$, Assertion 2 holds.

Proof.

By the second call to Prune we have that $\theta_{\ell}\leq\pi-\frac{c}{\nu N}$ with probability $1-\frac{\mu}{6}$ (as per Corollary 5.14). ∎

The following follows directly from Corollary 5.17:

Lemma 5.20**.**

With probability at least $1-\frac{\mu}{6}$, Assertion 3 holds.

Before arguing about Assertion 5, we first argue about correctness of Algorithm 8.

Lemma 5.21**.**

Assume Assertion 1-3 are true and let $w\leq\nu N$. Finally let $\delta_{0}=\frac{\sqrt{\epsilon}}{5000}$ and $\nu=\frac{1}{8}$. Then a call to $\textsc{QueryJacobiForCos}(w,2\nu N,R)$ returns $\frac{\cos(w\theta_{\ell})\pm O(\sqrt{\epsilon})}{1\pm O(\sqrt{\epsilon})}$ with probability $\exp\left({-\Omega(R)}\right)$. Further this call makes $3R$ queries and has time complexity $O(R)$.

Proof.

The claim on the query complexity follows from the statement of Algorithm 8. The rest of argument follows from showing that one can apply Corollary 5.9. Next, we show that for any $i\in[R]$, $v_{i}$ has the correct estimate with probability at least $\frac{2}{3}$. Then an application of the Chernoff bound would prove the claimed result.

Fix any round $i\in[R]$. We would like to show that (48) holds with high probability. Note that the value of $v_{i}$ in Algorithm 8 does not change if multiply the numerator and denominator by $\frac{1}{v\sqrt{w_{\ell}}}$. I.e. WLOG we will assume that we have query access to $\frac{\mathbf{y}[w]}{v\sqrt{w_{\ell}}}$ and the $\frac{\epsilon}{N^{3/2}}$ term in definition of $D$ would be replaced by $\frac{\epsilon}{v\sqrt{w_{\ell}}N^{3/2}}$. It suffices to show that all the conditions are met, where the queries $\widetilde{{\mathcal{J}^{({\alpha},{\beta})}_{{w}}}}(\cos{\theta_{\ell}})$ in equation (47) will be provided by $\mathbf{y}[w]/(v\sqrt{w_{\ell}})$ (due to the normalization (43).) In particular, we will show that the following conditions are satisfied:

• (42):

This follows from Assertion 2 and the fact that all of $\Delta-w,\Delta,\Delta+w$ are at least $\nu N$. This is achieved by making sure $w\leq\nu N$ (which is guaranteed by the call to ApproxArcCos having $\log(\nu N)-1$ as its first parameter) and making sure $\Delta\in[2\nu M,N-2\nu N]$ (guaranteed by the call $\textsc{QueryJacobiForCos}(w,2\nu N,R)$).

• (46):

Assertion 3 implies that this condition is met with probability at least

[TABLE]

for random $\Delta$. (In the above, the first inequality follows from Lemma 7.1 and our choice of $\nu$ while the last inequality follows from our choice of $\delta_{0}$ and the fact that $\epsilon\leq 1$.)

• (47):

We analyze the error $\widetilde{{\mathcal{J}^{({\alpha},{\beta})}_{{w}}}}(\cos{\theta_{\ell}})-\mathcal{J}^{({\alpha},{\beta})}_{{w}}(\cos{\theta_{\ell}})=\mathbf{y}[w]/(v\sqrt{w_{\ell}})-\mathcal{J}^{({\alpha},{\beta})}_{{w}}(\cos{\theta_{\ell}})=\mathbf{w}[w]/(v\sqrt{w_{\ell}})$ (by equations (39) and (43)). The expected squared error $\mathbf{w}[\Delta]^{2}$ for a random position $\Delta$ is at most $\frac{\epsilon^{2}|v|^{2}}{N(1-4\nu)}$. Thus, by Markov with probability at least $\frac{17}{18}$, the error $\mathbf{w}[\Delta]\leq\frac{\sqrt{18}\epsilon\left|v\right|}{\sqrt{N(1-4\nu)}}\leq\frac{6\epsilon\left|v\right|}{\sqrt{N}}$. The error $\mathbf{w}[\Delta]/(v\sqrt{w_{\ell}})$ is therefore $\frac{6\epsilon}{\sqrt{Nw_{\ell}}}$, which by Theorem 5.8 (more specifically (44)) is $O\left({\kappa(\theta_{\ell})\cdot\epsilon}\right)$. Similarly, for a randomly chosen $\Delta$ (and fixed $w$), the errors $\mathbf{w}[\Delta-w]$ and $\mathbf{w}[\Delta+w]$ satisfy the same bounds. By the union bound, (47) is satisfied for all three with probability at least $\frac{5}{6}$. Finally, we note that modifying the denominator changes $\mathbf{y}[\Delta]/(v\sqrt{w_{\ell}})$ by by at most another $\frac{\epsilon}{|v|\cdot N\cdot\sqrt{Nw_{\ell}}}\leq\frac{\epsilon}{\sqrt{Nw_{\ell}}}=O(\kappa(\theta_{\ell})\cdot\epsilon)$ (where in the inequality we used our assumption that $\left|v\right|\geq\frac{1}{N}$– recall the discussion on precision in Section 2.3).

Thus, we can apply (48) to conclude that with probability at least $\frac{2}{3}$:

[TABLE]

This proof is complete by noting that all the error terms are $O(\sqrt{\epsilon})$ due to Assertion 1 (and the fact that $\Delta-w\geq\nu N$ by construction) and our choice of $\delta_{0}$. ∎

We now argue that

Lemma 5.22**.**

Conditioned on Assertion 1-3 being true, with probability at least $1-\frac{\mu}{6}$, Assertion 5 holds.

Proof.

There are $O(\log{N})$ calls to $\textsc{QueryJacobiForCos}(w,2\nu N,R)$ for $R=O(\log\log{N}+\log(1/\mu))$. Thus, by Lemma 5.21 we know that each call gives back the desired output with probability at least $1-O(\mu/\log{N})$. Thus, by the union bound, all the calls work as intended with probability at least $1-\frac{\mu}{6}$. Then Assertion 5 follows from Theorem 5.11 (and noting that by Lemma 5.4 we have $\theta_{\ell}\in(0,\pi)$ as needed by Theorem 5.11). ∎

Lemma 5.23**.**

Conditioned on Assertion 5 being true, with probability at least $1-\frac{\mu}{6}$, Assertion 6 holds.

Proof.

This follows from Lemma 5.21 and Corollary 5.14. ∎

The following is a direct consequence of Lemma 5.13:

Lemma 5.24**.**

With probability at least $1-\frac{\mu}{6}$, the value $\tilde{v}$ returned by Algorithm 7 satisfies $\left|\tilde{v}-v\right|\leq 13\epsilon|v|$.

Lemma 5.18, 5.19, 5.20, 5.22, 5.23, and 5.24 implies the following result:

Lemma 5.25**.**

With probability at least $1-\mu$, Algorithm 7 correctly outputs $\ell$ and a $\tilde{v}$ such that $\left|\tilde{v}-v\right|\leq 13\epsilon|v|$.

Note that the above establishes the correctness of Algorithm 7.

Time and Query Complexity.

Lemma 5.26**.**

Algorithm 7 makes

[TABLE]

queries and takes time

[TABLE]

Proof.

We analyze the various function calls in Algorithm 7, which dominate the runtime. There are three calls to Prune each with query complexity $q_{\textsc{Prune}}\left({{\left({s,N,\frac{\mu}{6},\epsilon}\right)}}\right)$ and time complexity $T_{\textsc{Prune}}\left({{\left({s,N,\frac{\mu}{6},\epsilon}\right)}}\right)$. The first call has $s=O\left({\frac{1}{\epsilon}}\right)$ (using the fact that $c,C^{\prime},\nu$ are constants and our setting of $\delta_{0}=\Theta\left({\sqrt{\epsilon}}\right)$), the second has $s=O(1)$, and the third has $s=O(\epsilon^{1/4})$. These are dominated by the first call. The single call to PruneNonSpread, via Corollary 5.17 has query and time complexities of $O\left({\frac{1}{\sqrt{\epsilon}}\cdot q_{\textsc{Prune}}\left({{O(1),N,O(\sqrt{\epsilon}\mu),\epsilon}}\right)}\right)$ and $O\left({\frac{1}{\sqrt{\epsilon}}\cdot T_{\textsc{Prune}}\left({{O(1),N,O(\sqrt{\epsilon}\mu),\epsilon}}\right)}\right)$ respectively. Finally ApproxArcCos makes $O(\log{N})$ calls to QueryJacobiForCos with $R=O\left({\log\log{N}+\log\left({\frac{1}{\mu}}\right)}\right)$. The rest of ApproxArcCos takes $O(\log{N})$ time.121212The number of calls and time complexity follow from Theorem 5.11. Finally, by Lemma 5.21, we have that each call to QueryJacobiForCos takes time and has query complexity $O(R)$. Summing everything up, we get query complexity of

[TABLE]

(and time complexity with equal parameters but to $T_{\textsc{Prune}}\left({{}}\right)$ instead of $q_{\textsc{Prune}}\left({{}}\right)$).

Plugging the value of $q_{\textsc{Prune}}\left({{}}\right)$ from Lemma 5.12 yields a query complexity of

[TABLE]

and time complexity

[TABLE]

Theorem 5.8 (in particular, its implication that $U^{2}N=O(1)$) yields the statement. ∎

Lemma 5.25 and 5.26 implies the final result, Theorem 5.1.

6 A $k$-sparse Recovery Algorithm for Jacobi Polynomials

Finally, we use our $1$-sparse recovery algorithm (i.e. Theorem 5.1) along with Theorem 4.1 in order to obtain a $k$-sparse recovery algorithm for Jacobi polynomials. In order to apply Theorem 4.1, we first need to derive some properties of Jacobi polynomials.

Theorem 5.2 implies that that any family of Jacobi polynomials is dense as defined in Definition 2.2.

Corollary 6.1**.**

A family of Jacobi polynomials, with an associated constant $C^{\alpha,\beta}$ defined as in Theorem 5.2, is $(\frac{1}{2\pi},\frac{3}{2\pi},8C^{\alpha,\beta}\pi)$-dense.

Proof.

Consider the roots $\cos\left({\theta_{0}}\right),\dots,\cos\left({\theta_{d-1}}\right)$ of degree $d$ Jacobi polynomial. Theorem 5.2 says that for all $i\in[d]$,

[TABLE]

First, we can bound the number of roots falling in any interval $I=[a,b]$. Equation (60) implies that any root $\theta_{i}\in[a,b]$ satisfies $\frac{i}{d}\pi\in\left[a-\frac{\pi}{d}C^{\alpha,\beta},b+\frac{\pi}{d}C^{\alpha,\beta}\right]$, and therefore there are at most

[TABLE]

roots in $I$. Similarly, any $i$ such that $\frac{i}{d}\pi\in\left[a+\frac{\pi}{d}C^{\alpha,\beta},b-\frac{\pi}{d}C^{\alpha,\beta}\right]$ is a root, so there are at least

[TABLE]

roots in $I$.

Now choose $\gamma_{0}=8C^{\alpha,\beta}\pi$ in the statement of Definition 2.2. Fix any $\ell\in[d]$.

For any131313For clarity, note that this is not the constant $\gamma$ defined in Theorem 5.5, but a fresh variable in the statement of Definition 2.2. $\gamma\geq\gamma_{0}/d$, we are interested in bounding the number of roots falling in $\theta_{\ell}\pm\gamma/2$. Note that $\frac{\gamma}{2}\geq C^{\alpha,\beta}\frac{\pi}{d}$ by our choice of $\gamma_{0}$. Therefore the interval in question satisfies the containments

[TABLE]

The number of roots $r$ lying in this interval therefore satisfies

[TABLE]

Now, note that since $d\gamma\geq d\gamma_{0}\geq\frac{8C^{\alpha,\beta}}{\pi}$, we have

[TABLE]

and

[TABLE]

as desired.

∎

Finally, our following main result (which is the formal version of Theorem 3.1) follows by applying Theorem 4.1 along with Theorem 5.8 and Corollary 6.1 (with $C_{0}=\frac{1}{2\pi}$, $C_{1}=\frac{3}{2\pi}$ and $\gamma_{0}=8C^{\alpha,\beta}\pi$).

Corollary 6.2** ($k$-sparse recovery for Jacobi Transform).**

Fix parameters $\alpha,\beta>-1$. The following is true for small enough $\epsilon$ and some constant $C$. Let $\mathbf{J}^{(\alpha,\beta)}_{N}$ denote the $N\times N$ Jacobi transform (as defined in (43)) with associated constant $C^{\alpha,\beta}$ in Theorem 5.2. Then there exists an algorithm with the following property. Let $\hat{\mathbf{x}}=\mathbf{J}^{(\alpha,\beta)}_{N}\mathbf{x}$ be $(k,\frac{3}{2\pi}\gamma)$-sparsely separated for $\gamma\geq 8C^{\alpha,\beta}\pi/N$, and let $\mathbf{w}\in\mathbb{R}^{N}$ with $\hat{\mathbf{w}}=\mathbf{J}^{(\alpha,\beta)}_{N}\mathbf{w}$ be a “noise” vector so that $\left\|{\hat{\mathbf{w}}}\right\|_{{2}}\leq\frac{\delta}{2C}\textsc{large}(k,\hat{\mathbf{x}})$.

Then there exists an algorithm with the following property. for any $0<\mu<1$, the algorithm with probability at least $1-\mu$, outputs $\hat{\mathbf{z}}$ so that $\|\hat{\mathbf{x}}-\hat{\mathbf{z}}\|_{2}\leq O(\delta)\|\hat{\mathbf{x}}\|_{2}$. Further, the algorithm makes at most

[TABLE]

queries to $\mathbf{x}+\mathbf{w}$, runs in time at most

[TABLE]

7 Proof of Theorem 5.11

In this section, we complete the final missing piece by proving Theorem 5.11. We are looking for an unknown $\theta\in(0,\pi)$. The queries have the form

[TABLE]

Querying at $w=1$ gives a range (depending on the size of the noise) of possible $\theta$.

Now at beginning of every stage $t$ of the algorithm, we have a candidate interval $S_{t}$ which has width of order $2^{-t}$, which we know that $\theta$ lies in. By querying at $w=\Theta(2^{t})$, this interval is “dilated” so that the query $T_{w}(\theta)$ gives more information about $\theta$, and $S_{t}$ can be narrowed to $S_{t+1}$.

There is one obstacle in that knowing $\cos(w\theta)$ only reveals $w\theta$ up to sign. When the two possibilities are not well separated (that is, when $\cos(w\theta)$ is close to $\pm 1$), we do not receive enough information to narrow $S_{t}$. These cases can be overcome by re-querying at a different suitably chosen $w$.

ApproxArcCos

We begin with a bound on how much the noise in the query $\cos(\theta)$ affects the argument $\theta$.

Lemma 7.1**.**

If $\epsilon$ is sufficiently small, then $\rho_{\epsilon}\leq 2\sqrt{5\epsilon}$ satisfies Definition 5.10.

Proof.

Note that $\frac{1}{1\pm\epsilon}\in 1\pm 2\epsilon$ for $\epsilon<1$, so

[TABLE]

By Taylor expansion,

[TABLE]

for any $\theta\in[0,\pi]$. This completes the claim. ∎

We are now ready to prove Theorem 5.11.

Proof of Theorem 5.11.

We first note that the claim on the query locations follows from statement of Algorithm 9.

We will proceed by induction on $t$, with the inductive hypothesis that:

1. (a)

$\theta\in S_{t}$. 2. (b)

The width of $S_{t}$ is at most $|S_{t}|\leq\frac{2\rho_{\epsilon}}{2^{t}}$.

Note that the above with $t=\tau$ proves the claim on correctness of Algorithm 9.

For notational convenience, for the rest of the proof

Definition 7.2**.**

We will use $\rho$ to denote $\rho_{\epsilon}$. We will also extend the interval notation to let $\pm x\pm h$ denote $-x\pm h\cup x\pm h=[-x-h,-x+h]\cup[x-h,x+h]$.

Let $\phi(x)$ be the value $x+z\cdot 2\pi$ for some integer $z$ so that $\phi(x)\in[-\pi,\pi]$. We begin with an observation:

Observation 7.3**.**

In $\textsc{Refine}(w,S,\epsilon)$,

[TABLE]

Proof.

Since $w\theta$ can be outside of $[-\pi,\pi]$, we use the fact that $\cos(w\theta)=\cos\left({\phi(w\theta)}\right)$. If $\phi(w\cdot\theta)\in[0,\pi]$, then Lemma 7.1 gives

[TABLE]

Otherwise $\phi(w\cdot\theta)\in[-\pi,0]$ and

[TABLE]

∎

Base case.

Now we establish the base case, for $t=0$. Item (b) holds by construction of $S_{0}$, so the only things to check are that (a) holds, and also that the Assertion 1 in the pseudocode is correct.

By Observation 7.3, in Refine$(1,[0,\pi],\epsilon)$, we set

[TABLE]

Since $\theta\in[0,\pi)$ and $r\in[0,\pi]$, the sign must be positive, and we have

[TABLE]

We claim that we will set $a_{0}$ to be $a_{0}=r$. Certainly this value lives in $\bar{S}=[0,\pi]$, which ensures that $a_{0}$ will not be set to None (and hence Assertion 1 holds). Moreover, no other value $r+z\cdot 2\pi$ will be used for $a_{0}$, since for any $z\neq 0$, $r+z\cdot 2\pi$ does not live in $[0,\pi]$. Thus $a_{0}$ is as claimed.

Then $|a_{0}-\theta|\leq\rho$, which means that the choice of $S_{0}=[a_{0}-\rho,a_{0}+\rho]$ indeed satisfies (a).

Establishing the inductive hypothesis.

With the base case out of the way we proceed by induction.

Claim 7.4**.**

Suppose that the inductive hypothesis holds for $t-1$, and suppose that $w\leq\frac{3}{2}\cdot 2^{t}$, and that $\rho$ is sufficiently small (say, $\rho<\pi/4$). Then consider

[TABLE]

Then at least one of the following holds:

[TABLE]

Proof.

By Observation 7.3,

[TABLE]

which means that for some $z\in\mathbb{Z}$, and some sign $\zeta\in\{-1,1\}$,

[TABLE]

Suppose that $\zeta=+1$. Then

[TABLE]

which implies that (using the inductive hypothesis that $\theta\in S_{t-1}$) we have

[TABLE]

Moreover, this value of $z$ is the unique one with this property, since

[TABLE]

using the fact that $w\leq\frac{3}{2}\cdot 2^{t}$ and that $\rho<\pi/4$ is sufficiently small. Thus, for any other $z^{\prime}$,

[TABLE]

This implies that, in the case where $\zeta=+1$, we will choose this $z$ in our definition of $a_{t}$, and hence

[TABLE]

Similarly, if $\zeta=-1$, we will choose

[TABLE]

Since $\zeta$ is either $-1$ or $+1$, one of the two cases will hold and this proves the claim. ∎

With this claim out of the way, we will establish the inductive hypotheses for the next round. Again (b) holds by construction of $S_{t}$, so we focus on (a). Consider $a_{t},b_{t}=\textsc{Refine}(w=2^{t},S=S_{t-1},\epsilon)$, the output of the first call to Refine.

Case 0: at least one of $a_{t}$, $b_{t}$ is None. If at least one of $a_{t},b_{t}$ is None, the Claim 7.4 implies that one of them (say, $c_{t}$) is not None and that

[TABLE]

and so in particular $\theta\in S_{t}$, establishing the inductive hypothesis (a).

Case 1: both $a_{t}$ and $b_{t}$ are not None. Since they are both in $\overline{S_{t-1}}$, by the inductive hypothesis $|a_{t}-b_{t}|\leq\left|\overline{S_{t-1}}\right|\leq|S_{t-1}|+2\rho 2^{-t}\leq 2\rho{2^{-(t-1)}}+2\rho 2^{-t}=6\rho{2^{-t}}$. This means that, if $z$ is chosen in the definition of $a$ and $z^{\prime}$ for $b$, then

[TABLE]

Thus,

[TABLE]

or

[TABLE]

which implies that (since $r\in(0,\pi)$ and $(z-z^{\prime})$ is an integer),

[TABLE]

On the other hand, by Observation 7.3,

[TABLE]

so, it follows that for some integer $h$,

[TABLE]

Thus $\frac{\pi h}{2^{t}}=\theta\pm 4\frac{\rho}{2^{t}}\in S_{t-1}\pm 4\frac{\rho}{2^{t}}$. This interval has width at most $12\rho 2^{-t}$, so this integer $h$ can be found from $S_{t-1}$ assuming that $\rho<\pi/12$ is sufficiently small (and this $h$ is unique).

Assertion 2 holds.

Since $\theta\in(0,\pi)$, this implies that $h\in\{0,1,\ldots,2^{t}\}$. In fact, to prove Assertion 2, we will show that $h$ cannot be [math] or $2^{t}$.

Suppose for the sake of contradiction that $h$ is zero. Then $\theta\in[-4\rho 2^{-t},4\rho 2^{-t}]$, which means that $\theta\in(0,4\rho 2^{-t}]$, since $\theta\in(0,\pi)$. This in turn implies that

[TABLE]

where the first inclusion follows from the fact that $|S_{t-1}|\leq 4\rho 2^{-t}$.

By Observation 7.3, $r=\phi(2^{t}\cdot\theta)\pm\rho\in(0,5\rho]$ (provided $4\rho<\pi$). We note that we also used the fact that $r\in(0,\pi)$ as well as $\theta\geq 0$ to conclude that the sign in front of $\phi$ should be positive.

Note that if $14\rho<2\pi$, then $\frac{2\pi-r}{2^{t}}>\frac{9\rho}{2^{t}}$. In the call to Refine, since $\overline{S_{t-1}}\subseteq(0,9\rho 2^{-t}]$ (which in turn follows from (62)), this means $b_{t}$ would have been set to None.141414Note that $\frac{-r}{2^{t}}+z\cdot\frac{2\pi}{2^{t}}<0$ for any $z\leq 0$ (since $r>0$) and is $>\frac{9\rho}{2^{t}}$ for any $z\geq 1$ (since we showed it for $z=1$ and $z>1$ can only give a larger value) so this exhausts all possibilities. Thus the case $h=0$ could not have happened. In the case $h=2^{t}$, we similarly bound

[TABLE]

The claim on $r$ is due to the following argument. Note that $2^{t}\theta\in[2^{t}\pi-4\rho,2^{t}\pi]$, which implies $\phi\left({2^{t}\theta}\right)\in[-4\rho,0]$, which in turn implies that $-\phi\left({2^{t}\theta}\right)\in(0,4\rho]$ and this is the interval that gets ‘picked up’ when intersecting with $(0,\pi)$.

We now argue that both $\frac{r}{2^{t}}+2^{t-1}\frac{2\pi}{2^{t}}$ and $\frac{r}{2^{t}}+(2^{t-1}-1)\frac{2\pi}{2^{t}}$ are not in $\overline{S_{t-1}}$. Indeed $\frac{r}{2^{t}}+2^{t-1}\frac{2\pi}{2^{t}}>\pi$ and hence $\frac{r}{2^{t}}+2^{t-1}\frac{2\pi}{2^{t}}\not\in\overline{S_{t-1}}$. Next we argue that

[TABLE]

which implies that the LHS is not in $\overline{S_{t-1}}$ as needed. Now note that the above inequality is true if $r-2\pi<-9\rho$, which is true if $\rho<\pi/7$ (where we also used the fact that $r\leq 5\rho$). The above also implies that for $z\in\mathbb{Z}$, $\frac{r}{2^{t}}+z\cdot\frac{2\pi}{2^{t}}$ is also not in $\overline{S_{t-1}}$ since the value is $>\pi$ for $z\geq 2^{t-1}$ (since it is true for $z=2^{t-1}$ as $r>0$) and the value is $<\pi-9\rho 2^{-t}$ for $z\leq 2^{t-1}-1$ (since it is true for $z=2^{t-1}-1$). This implies that $a_{t}$ would have been None. So $h=2^{t}$ is also impossible.

Therefore $h\in\{1,\dots,2^{t}-1\}$, and write $h=2^{j}\cdot x$ where $x$ is odd and $j<t$. This establishes the Assertion 2 in the pseudocode.

Assertion 3 holds.

Now, we re-choose

[TABLE]

Since $t>j$, $w$ is an integer. Further, $2^{t}\leq w\leq\frac{3}{2}\cdot 2^{t}$, and so (because of our bound on $\tau$ in Theorem 5.11), $w$ is an integer between [math] and $N-1$ inclusive and thus is a valid query.

By (61),

[TABLE]

so we have

[TABLE]

for some odd integer $x^{\prime}$, using the fact that $w\leq\frac{3}{2}\cdot 2^{t}$.

Now consider the $r$ that is defined in the second call to Refine, with $w$ as in (63). By Observation 7.3 again,

[TABLE]

Now for any choice of $z,z^{\prime}$,

[TABLE]

where $z^{\prime\prime}\in\mathbb{N}$. This implies that there cannot both be an $a_{t},b_{t}$ in $\overline{S_{t-1}}$ as long as

[TABLE]

which holds provided $\rho<\pi/22$. Thus at least one of $a_{t},b_{t}$ in the second call to Refine is None, and by Claim 7.4 the other (call it $c_{t}$) exists and satisfies

[TABLE]

This establishes Assertion 3. Thus, when we define $S_{t}=[c_{t}-\rho/2^{t},c_{t}+\rho/2^{t}]$, we guarantee that $\theta\in S_{t}$, establishing the inductive hypothesis (a) for the next round.

Query and Runtime analysis.

It is easy to see that Algorithm 9 makes at most $2\tau$ calls to Algorithm 10, which implies that it makes at most $2\tau$ queries. It is also easy to check that other than calls to Refine, the rest of ApproxArcCos runs in time $O(\tau)$. To finish the proof, we claim that each call to Refine takes $O(1)$ time.151515This assumes computing the $\arccos$ values takes $O(1)$ time. The main operation is computing the values $a$ and $b$. To compute $a$, it suffices to find any $z\in\mathbb{Z}$ satisfying $L\leq\frac{r}{w}+z\cdot\frac{2\pi}{w}\leq R$, or any integer in $[wL/2\pi-r/2\pi,wR/2\pi-r/2\pi]$. The case for $b$ is similar. ∎

Acknowledgments

We would like to thank Mark Iwen for useful conversations. Thanks also to Stefan Steinerberger for showing us the proof of Lemma E.3 and graciously allowing us to use it, and to Clément Canonne for helpful comments on our manuscript.

ACG is partially funded by a Simons Foundation Fellowship. AR is partially funded by NSF grant CCF-1763481. MW is partially funded by NSF CAREER award CCF-1844628. AG and CR gratefully acknowledge the support of DARPA under Nos. FA87501720095 (D3M) and FA86501827865 (SDH), NIH under No. U54EB020405 (Mobilize), NSF under Nos. CCF1763315 (Beyond Sparsity) and CCF1563078 (Volume to Velocity), ONR under No. N000141712266 (Unifying Weak Supervision), the Moore Foundation, NXP, Xilinx, LETI-CEA, Intel, Google, NEC, Toshiba, TSMC, ARM, Hitachi, BASF, Accenture, Ericsson, Qualcomm, Analog Devices, the Okawa Foundation, and American Family Insurance, Google Cloud, Swiss Re, Stanford Bio-X SIG Fellowship, and members of the Stanford DAWN project: Intel, Microsoft, Teradata, Facebook, Google, Ant Financial, NEC, SAP, VMWare, and Infosys. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views, policies, or endorsements, either expressed or implied, of DARPA, NIH, ONR, or the U.S. Government.

Appendix A Chebyshev Reduction to sparse FFT

For completeness, in this section we show how to solve the sparse approximation problem for Chebyshev polynomials via a reduction to the sFFT. We note that a similar observation has been made in [legendre-sparse].

In this section, everything is [math]-indexed, and negative indices count from the end. For example, if $\mathbf{x}$ has length $N$ then $\mathbf{x}[-i]=\mathbf{x}[N-i]$.

A.1 Facts about DFT

We begin with some useful facts about the Discrete Fourier Transform (DFT). Consider the usual DFT $\hat{\mathbf{x}}=\mathcal{F}\mathbf{x}$ defined by

[TABLE]

where $\mathbf{x}\in\mathbb{C}^{N}$ and $\omega_{N}:=\exp(-2\pi i/N)$.

Shifting $\mathbf{x}$ scales $\hat{\mathbf{x}}$ elementwise, and vice versa:

Consider $\mathbf{x}\in\mathbb{C}^{N}$ and $\mathbf{y}$ such that $\mathbf{y}[\ell]=\mathbf{x}[\ell-s]$ for all $\ell\in[N]$. Then

[TABLE]

If $\mathbf{y}[\ell]=\mathbf{x}[\ell]\omega^{\ell s}$, then

[TABLE] 2. 2.

If $\mathbf{x}$ is supported only on the first half, then the even indices of $\hat{\mathbf{x}}$ can be computed by a DFT of half the size, and vice versa:

[TABLE]

[TABLE]

for $j=0,\dots,N/2-1$. 3. 3.

$\mathbf{x}$ is symmetric ($\mathbf{x}[i]=\mathbf{x}[-i]$) if and only if $\hat{\mathbf{x}}$ is symmetric.

[TABLE]

A.2 Sparse DCT Setup

Consider the Chebyshev transform $\hat{\mathbf{c}}=\mathcal{C}\mathbf{c}$

[TABLE]

where $\lambda_{j}=\cos\left(\pi\frac{2j+1}{2n}\right)$ are the Chebyshev nodes. Note that $T_{\ell}(\lambda_{j})=\cos\left(\frac{\pi}{2N}\ell(2j+1)\right)$ by properties of Chebyshev polynomials.

A.3 Reduction

We will prove the following equivalence:

Lemma A.1**.**

For any $\mathbf{c}\in\mathbb{C}^{N}$ and its transform $\hat{\mathbf{c}}=\mathcal{C}\mathbf{c}$ there exists an $\mathbf{f}\in\mathbb{C}^{2N}$ and $\hat{\mathbf{f}}=\mathcal{F}\mathbf{f}$ so that one can express each entry in $\mathbf{f}$ (and $\hat{\mathbf{f}}$) as a scalar multiple of an entry in $\mathbf{c}$ (and $\hat{\mathbf{c}}$ resp.). Further, $\mathbf{c}$ is k-sparse (or $\hat{\mathbf{c}}$ is k-sparse) if and only if $\mathbf{f}$ is 2k-sparse (or $\hat{\mathbf{f}}$ is 2k-sparse resp.).

This implies that ${\cal C}$ (and its inverse) have sparse recovery algorithms by reducing to a sparse FFT.

Proof.

We will define $\mathbf{x},\mathbf{y},\mathbf{z}\in\mathbb{C}^{4N}$ and $\mathbf{f}\in\mathbb{C}^{2N}$ as follows. Let $\omega:=\omega_{4N}$.

Define $\mathbf{x}\in\mathbb{C}^{4N}$ by

[TABLE]

Define $\mathbf{y}[j]=\mathbf{x}[j-N]$ for all $j\in[4N]$, More specifically,

[TABLE]

Note that $\mathbf{y}$ is only supported on indices $0,\dots,2N-1$. Finally, define $\mathbf{z}[j]=\omega^{j}\mathbf{y}[j]$ for all $j\in[4N]$ and $\mathbf{f}[j]=\mathbf{z}[j]$ for all $j\in[2N]$. By construction, every entry of $\mathbf{f}$ is a scalar multiple of an entry of $\mathbf{c}$. Furthermore for every $j\in[N]$ exactly two entries of $\mathbf{f}$ depend on $\mathbf{c}[j]$ (unless $j=0$, in which case one entry of $\mathbf{f}$ depends on $\mathbf{c}[0]$), so $\mathbf{c}$ is $k$-sparse if and only if $\mathbf{f}$ is $2k$-sparse.

Now we analyze $\hat{\mathbf{f}}$ and show that it has the same properties.

Claim A.2**.**

[TABLE]

for $j=0,\dots,N-1$.

Proof.

Writing $2\cos(z\pi/2N)=\omega_{4N}^{z}+\omega_{4N}^{-z}$, we have

[TABLE]

for $j=0,\dots,N-1$. In the above, the first equality follows from the definition of DCT while the last equality follows from definition of $\mathbf{x}$.

Finally, $\mathbf{x}$ is symmetric, so by fact 3 $\hat{\mathbf{x}}$ is as well. ∎

Next, by fact 1, $\hat{\mathbf{y}}[j]=\hat{\mathbf{x}}[j]\omega^{jN}$ for all $j\in[4N]$. Also by fact 1, $\hat{\mathbf{z}}[j]=\hat{\mathbf{y}}[j+1]$ for all $j\in[4N]$. Finally, by fact 2,

[TABLE]

for every $j\in[2N]$ (note here that $\hat{\mathbf{f}}$ is computed from a Fourier transform of size $2N$ instead of $4N$ since their dimensions differ). Here the last equality follows from (65).

Therefore for every $j\in[N]$, two entries of $\hat{\mathbf{f}}$ are a scalar multiple of $\hat{\mathbf{c}}[j]$. This concludes the proof of Lemma A.1.

∎

Appendix B Why not apply Jacobi to Chebyshev reduction directly to $k$-sparse recovery?

In this appendix we outline an approach to directly reduce $k$-sparse recovery of Jacobi to $k$-sparse recovery for DCT (and hence sFFT) and outline the main challenges in carrying out the reduction. The hope here is to give the reader a sense for why some obvious generalization of the ideas in Appendix A do not work for Jacobi polynomials.

The main idea is to use Theorem 5.5 to reduce evaluations of Jacobi polynomials to (approximately) evaluations of Chebyshev polynomials. One might hope then that one could directly reduce the $k$-sparse recovery problem for Jacobi to an $O(k)$-sparse recovery for Chebyshev transform. One obvious issue is that this reduction (via Theorem 5.5) is approximate and we need to handle the approximation error in our final analysis. However, for this discussion, let us assume this can be handled.

The major issue is that we now have to deal with Chebyshev polynomials evaluated at roots of Jacobi polynomials. There are two natural options here:

First observe that most roots of Jacobi polynomials are close enough to the roots of the Chebyshev polynomials. In particular, by a tighter version of Theorem 5.2 for the special case of $-1/2\leq\alpha,\beta\leq 1/2$, we know that the $\ell$th root is given by $\cos{\theta_{\ell}}$, where $\theta_{\ell}=\frac{(\ell+\tilde{\ell})\pi}{N+\frac{\alpha+\beta+1}{2}}$, for some $0<\tilde{\ell}<1$ (Theorem 6.3.2 in [szego]). So at least for the more restricted case of $-1/2\leq\alpha,\beta\leq 1/2$ (note our results hold for any fixed $\alpha,\beta>-1$) the $\ell$’th Chebyshev and Jacobi roots are fairly close. However, the crucial difference is that all the Chebyshev roots are equispaced, which allows us to permute the $\theta_{\ell}$ by just applying a permutation on $[N]$. More precisely for any permutation $\sigma$ on $[N]$, note that $\cos\left({\left({\frac{\sigma(\ell)+\frac{1}{2}}{N}}\right)\pi}\right)$ is the $\sigma(\ell)$’th root. However, since we do not have a closed form expression for Jacobi polynomial roots, it is not clear how to apply similar tricks. Further, the argument in Appendix A very crucially uses that the Chebyshev roots are equally spaced and trying to generalize those arguments to Jacobi polynomial roots breaks down immediately. In summary, we do not see how to directly reduce the Jacobi transforms to DCT or DFT. 2. 2.

Like in Section 5, we just try and do some “error-correction” to obtain cosine evaluations at regular points. Indeed, this is what the trick of $\cos{A}=\frac{\cos(A+B)+\cos(A-B)}{2\cos{B}}$ for random $B$ helped us do for the special case of $k=1$. However, this approach fails even for $k=2$ since we no longer have ‘common terms’ that cancel like they did in proof of Corollary 5.9.

Appendix C Known RIP results

We begin with some notation. Given a matrix $\mathbf{M}$, we denote $\max_{i,j}|\mathbf{M}[i,j]|$ by $\|\mathbf{M}\|_{\infty}$. Given an $N\times N$ diagonal matrix $\mathbf{D}$ and a subset $S\subseteq[N]$, $\mathbf{D}_{S}$ denotes the diagonal matrix given by

[TABLE]

Finally for any matrix $\mathbf{M}$, we define $\sqrt{\mathbf{M}}$ to be the matrix where we apply the $\sqrt{\cdot}$ operator to each entry of $\mathbf{M}$.

We will need the following well-known result:

Theorem C.1** ([RV]).**

Let $\mathbf{A}$ be an $N\times N$ orthogonal matrix (i.e. $\mathbf{A}^{T}\mathbf{A}=\mathbf{I}$) and let $\|\mathbf{A}\|_{\infty}=U$. Let $r\geq 1$ be an integer and $0<\delta<1$ be a real. Define an integer $1\leq s\leq N$ such that

[TABLE]

Pick a random subset $S\subseteq[N]$ where each element of $[N]$ is picked iid with probability $\frac{s}{N}$. Define the matrix

[TABLE]

Then the following holds:

[TABLE]

We will use the above result for $r=2$ to argue the following well-known corollary that the matrix $\mathbf{B}$ essentially preserves the inner products of columns of $\mathbf{A}$:

Corollary C.2**.**

Let $s$ be an integer such that

[TABLE]

The consider $\mathbf{A}$, $S$ and $\mathbf{B}$ be as defined in Theorem C.1 with $r=2$. Then we have for any $0<\alpha<1$,

[TABLE]

Proof.

We apply Theorem C.1 with $r=2$. Then by Markov’s inequality, we have that with probability at least $1-\alpha$, we have for any $i,j\in[N]$:

[TABLE]

and

[TABLE]

Note that

[TABLE]

Using the above we have

[TABLE]

The above along with (69) and (70) implies that for every $i,j\in[N]$,

[TABLE]

which completes the proof (since $\left(\mathbf{B}^{T}\mathbf{B}\right)[i,j]=\mathbf{e}_{i}^{T}\mathbf{B}^{T}\mathbf{B}\mathbf{e}_{j}$). ∎

Appendix D More on Jacobi polynomials

We begin by collecting more results on Jacobi polynomials and more generally orthogonal polynomials in Appendices D.1 and D.2 respectively. Finally, we prove Theorem 5.8 in Appendix D.3.

D.1 More results on Jacobi polynomials

We start off with the definition of the Bessel function:

Definition D.1**.**

The Bessel function of the first kind for parameter $\alpha\in\mathbb{R}$ is defined as follows. For every161616The Bessel function is also defined for $z\in\mathbb{C}$ but since we only need it over reals, we only state it for reals. $z\in\mathbb{R}$:

[TABLE]

We will need the following approximation and bounds for the Bessel function:

Lemma D.2**.**

Let $\alpha>-1$. For large enough $z$ (compared to $\alpha$), we have

[TABLE]

For all fixed $0<\epsilon_{0}<\sqrt{3(\alpha+1)}$ and every $z\in\left[{\epsilon_{0},\min\left\{2,\sqrt{3(\alpha+1)}\right\}}\right]$, we have

[TABLE]

For all fixed $C_{0}>0$ and every $0\leq z\leq C_{0}$, we have

[TABLE]

Finally, for all $z\geq 0$, we have

[TABLE]

Proof.

Equation (71) is equation (1.71.7) in [szego]. Further, (74) follows by combining (73) and (71). In the rest of the proof, we argue (72) and (73).

We begin with (73). We will argue that there is a large enough $\nu_{0}$ (that depends on $C_{0}$) for which there is a constant $\gamma_{0}$ (such that $|\gamma_{0}|<1$) with the following property. For every $\nu\geq\nu_{0}$, we have that

[TABLE]

where

[TABLE]

Notice (75) implies that

[TABLE]

Since $\alpha>-1$, every term in the sum is also $O(1)$ and hence (73) holds.

Next, we argue (75). Indeed,

[TABLE]

By Stirling’s approximations of the Gamma function, we have $\nu!\geq\sqrt{2\pi}\nu^{\nu+\frac{1}{2}}\cdot e^{-\nu}$ and $\Gamma(\nu+\alpha+1)\geq\Gamma(\nu)\geq\sqrt{2\pi}\nu^{\nu-\frac{1}{2}}\cdot e^{-\nu}$ (the first inequality follows since $\alpha>-1$ and the fact that $\Gamma(y)$ is an increasing function for large enough $y$). This along with the above equation implies that

[TABLE]

as desired. In the above, the last inequality follows for $\nu$ being large enough compared to $C_{0}$ and $\alpha$.

Finally, we argue (72). Note that (73) implies we only need to prove a lower bound. Consider for $\nu\geq 2$:

[TABLE]

where the inequality follows by our choice of $z$ and the fact that $\Gamma(y)$ is increasing for $y\geq 2$ (note that by our choice of $\alpha$ and $\nu$, $\nu+\alpha+1\geq 2$). Thus, for every even $\nu\geq 2$, we have:

[TABLE]

This in turn implies that

[TABLE]

as desired. In the above, the equality follows from the fact that for any $y>0$, $\Gamma(y+1)=y\Gamma(y)$; the second inequality follows from our assumption that $z\leq\sqrt{3(\alpha+1)}$ while the final inequality follows since we assumed $z\geq\epsilon_{0}$. ∎

We will use the Bessel function to present a tighter approximation of the Jacobi polynomials than Theorem 5.5. Towards this end, we begin with the following result:

Theorem D.3** ([szego], Theorem 8.21.12).**

Let $\alpha>-1$ and $\beta$ be an arbitrary real. We have for every $n\geq 1$:

[TABLE]

where

[TABLE]

Further, we define

[TABLE]

We will use the following slight re-statement of the above result:

Corollary D.4**.**

Let $\alpha>-1$ and $\beta$ be an arbitrary real. For every $n\geq 1$:

[TABLE]

for any

[TABLE]

Further, in the above $g_{{\alpha}}\left({n}\right)$ is as in (76) and

[TABLE]

Proof.

This follows from Theorem D.3 by using the definition of $\kappa_{\alpha,\beta}\left({\theta}\right)$, the fact that $\sin\theta=2\sin\frac{\theta}{2}\cos\frac{\theta}{2}$ and the observation that $\theta\sin\theta$ is $O(1)$. ∎

We also need an approximation for the range $\theta\in\left({\frac{\pi}{2},\pi}\right)$, which needs the following result:

Lemma D.5** (Equation (4.1.3) in [szego]).**

[TABLE]

We are now ready to state the equivalent of Corollary D.4 for obtuse $\theta$:

Corollary D.6**.**

Let $\beta>-1$ and $\alpha$ be an arbitrary real. For every $n\geq 1$, $P_{n}^{(\alpha,\beta)}(\cos\theta)$ is equal to

[TABLE]

for any

[TABLE]

Further, in the above $g_{{\beta}}\left({n}\right)$ is as in (76) (with $\alpha=\beta$) and $\kappa_{\alpha,\beta}\left({\theta}\right)$ is as in (77).

Proof.

Let

[TABLE]

Note that $\cos\theta=-\cos\theta^{\prime}$ and thus, by Lemma D.5, we have

[TABLE]

The above with Corollary D.4 applied with $\theta=\theta^{\prime}$ along with the observation that $\kappa_{\beta,\alpha}\left({\theta^{\prime}}\right)=\kappa_{\alpha,\beta}\left({\theta}\right)$ completes the proof. ∎

Finally, we want to argue that $g_{{\alpha}}\left({n}\right)$ is no more than a constant:

Lemma D.7**.**

For large enough $n$,

[TABLE]

Proof.

This follows from the Stirling’s approximation of the $\Gamma$ function. Indeed we have

[TABLE]

It is easy to check that each of the terms above are $\Theta_{\alpha,\beta}\left({1}\right)$, which completes the proof. ∎

We now have the final approximation for the Jacobi polynomial that we use in our proofs:

Corollary D.8**.**

Let $\alpha,\beta>-1$. For every $n\geq 1$:

[TABLE]

for any

[TABLE]

where

[TABLE]

Proof.

The proof follows from Corollaries D.4, D.6 and Lemma D.7. ∎

We will also need the following result from [szego]:

Lemma D.9** (Equation (4.21.7) in [szego]).**

For any integer $j\geq 1$:

[TABLE]

D.2 Christoffel-Darboux formula

In this section, we will only consider OPs that are orthonormal. In such a case the following result is well-known:

Theorem D.10** ([szego], Thm 3.2.1).**

Let $p_{0}(X),p_{1}(X),\dots$ be an OP family defined as in (1). Then for any $j\geq 2$, we have

[TABLE]

where

[TABLE]

We will need the following bound on $a_{j}$ for (orthonormal version of) Jacobi polynomials:

Remark D.11**.**

We note that $a_{j}=2+o(1)$ for $\mathcal{J}^{({\alpha},{\beta})}_{{j}}(X)$, where this claim follows from (40) and (41).

Next, we will need the following definition:

Definition D.12**.**

Given a family of OP $p_{0}(X),p_{1}(X),\dots$ as in (1) and an integer $n\geq 0$, define the kernel polynomial

[TABLE]

The following result is well-known (we provide its proof for the sake of completeness):

Theorem D.13** (Christoffel-Darboux formula, Theorem 3.2.2 [szego]).**

For any $\alpha\neq\beta$ and any $n\geq 1$, we have

[TABLE]

Proof.

We will first prove that for any $j\geq 0$:

[TABLE]

where for notational convenience we define $p_{i}(X)=0$ for any $i<0$. Note that summing the above for $0\leq j<n$ proves the claimed result (where we use the fact that $p_{-1}(X)=0$).

To prove (80), we first use (78) with $j+1$ instead of $j$ and evaluate all the polynomials at $\alpha$ and $\beta$ to get

[TABLE]

and

[TABLE]

Multiplying the above equations by $p_{j}(\beta)$ and $p_{j}(\alpha)$ respectively and then subtracting them one gets

[TABLE]

Re-arranging the above, one gets

[TABLE]

The above along with (79) implies (80). ∎

We will use the following corollary for the above result:

Corollary D.14** ((3.2.4) in [szego]).**

For any $\alpha$ and any $n\geq 1$, we have

[TABLE]

D.3 Proof of Theorem 5.8

We are now ready to argue that the OP transform corresponding to Jacobi polynomials for $\alpha,\beta>-1$ is flat, i.e. we prove Theorem 5.8 (which we re-state below for convenience):

Theorem D.15** (Theorem 5.8, restated).**

Let $\alpha,\beta>-1$. Let $N\geq 1$ be large enough. Let $\lambda_{0},\dots,\lambda_{N-1}$ be the roots of the $N$th Jacobi polynomial. Then define $\mathbf{F}$ such that for every $0<\ell,j<N$, we have

[TABLE]

where

[TABLE]

Then for every $0\leq\ell<N$, we have

[TABLE]

and

[TABLE]

Proof.

We will argue that for any $0\leq\ell,j<N$, we have

[TABLE]

The definition of the normalized Jacobi polynomials and Lemma 5.7 implies that

[TABLE]

The above proves both of the claimed results. Indeed the upper bound on $1/w_{\ell}$ follows by summing up the square of the above bound for all $0\leq j<N$. The above along with the bound of $O_{\left|\alpha\right|,\left|\beta\right|}\left({\frac{1}{N\kappa_{\alpha,\beta}\left({\theta_{\ell}}\right)^{2}}}\right)$ on $w_{\ell}$ from Lemma D.16 would prove our desired bound on entries of $\mathbf{F}$.

Finally, we note that (81) follows from applying (74) to the bound in Corollary D.8. ∎

Lemma D.16**.**

Let $\alpha,\beta>-1$. Let $N\geq 1$ be large enough. Let $\lambda_{0},\dots,\lambda_{N-1}$ be the roots of the $N$th Jacobi polynomial. Then

[TABLE]

The above result follows immediately from the following lemma and Lemma 5.4:

Lemma D.17**.**

Let $\alpha,\beta>-1$. Let $N\geq 1$ be large enough. Let $\lambda_{0},\dots,\lambda_{N-1}$ be the roots of the $N$th Jacobi polynomial. Then for every constant $c>0$, the following holds. Let $\theta_{\ell}$ be such that

[TABLE]

Then

[TABLE]

Proof.

We will first prove the result for every fixed $c\geq C_{0}$, where $C_{0}$ is a constant that we will fix later. Then we will prove the result for any constant $0<c<C_{0}$.

We start with the case of large $c$. We first observe that by Definition D.12, we need to lower bound $K_{N-1}(\lambda_{\ell},\lambda_{\ell})$. Then we use Corollary D.14 and the fact that $\lambda_{\ell}$ is a root of the $N$th Jacobi polynomial, to note that

[TABLE]

Remark D.11 and Lemma 5.7 imply that for large enough $N$, $\frac{1}{a_{N}\sqrt{h_{N}^{\alpha,\beta}\cdot h_{N-1}^{\alpha,\beta}}}=\Theta_{\left|\alpha\right|,\left|\beta\right|}\left({N}\right)$. Thus, to complete the proof, we need to argue that

[TABLE]

Towards this end, we recall the following identity (which appears as equation (4.5.7) in [szego]) that holds for any $n\geq 1$:

[TABLE]

Using the above for $n=N$, substituting $X=\lambda_{\ell}$ and noting that $1-\lambda_{\ell}^{2}=\sin^{2}\theta_{\ell}$, we get:

[TABLE]

This implies that

[TABLE]

where the second inequality follows from Lemma D.9 (with $j=N$).

Thus to show $E_{\ell,\alpha,\beta,N}\geq\Omega_{\left|\alpha\right|,\left|\beta\right|}\left({\kappa_{\alpha,\beta}\left({\theta_{\ell}}\right)^{2}}\right)$, it is enough to argue

[TABLE]

where the equality follows from the fact that $\kappa_{\alpha+1,\beta+1}\left({\theta}\right)=\frac{\kappa_{\alpha,\beta}\left({\theta}\right)}{\sin\frac{\theta}{2}\cos\frac{\theta}{2}}=\frac{2\kappa_{\alpha,\beta}\left({\theta}\right)}{\sin\theta}$. Corollary D.8 (with $n=N-1$ and where $\alpha\leftarrow\alpha+1$ and $\beta\leftarrow\beta+1$) implies that the above is true if

[TABLE]

where in the above we have used the fact that $N$ is large enough and where $\zeta$ is as defined in Corollary D.8. For notational simplicity, define

[TABLE]

Let $z$ be large enough so that (71) holds. I.e. there exist constant $C_{0}$ such that

[TABLE]

and

[TABLE]

We will argue that

[TABLE]

The above implies that it is enough to show that $1-O\left({\frac{1}{z}}\right)\geq\Omega_{\left|\alpha\right|,\left|\beta\right|}\left({1}\right)$, which is true if we pick $C_{0}$ to be large enough.

Now, we argue (82). Noting that $P_{N}^{(\alpha,\beta)}(\lambda_{\ell})=0$, Corollary D.8 implies that

[TABLE]

Since we have assumed $z$ is large enough and (71) holds, we have that

[TABLE]

In other words, we have

[TABLE]

where in the above we used the fact that $N\geq\Omega(z)$. The above, along with the fact that $\left|\sin A\right|\geq 1-\left|\cos A\right|$, implies (82).

If the constant $c$ in the lemma statement satisfies $c\geq C_{0}$, then we are done. So let us assume that $c<C_{0}$ and WLOG assume $\theta_{\ell}=\frac{c}{N}$. We will use the other bounds from Lemma D.2 to prove the lemma for this case. Towards that end, define

[TABLE]

and for any $0\leq j<N$,

[TABLE]

Then by our choice of $\epsilon_{0}$ and $c$, for some small enough constant $\mu>0$, there exists a subset $S\subseteq[0,N-1]$ with

[TABLE]

such that for every $j\in S$, we have $j\geq\mu N$ and:

[TABLE]

(Indeed, the above range is of size at least $\epsilon_{0}$ and every increment in $j$ increases $z_{j}$ by at most $\frac{c}{N}$, which means that are at least $\frac{\epsilon_{0}N}{c}-2$ $z_{j}$’s in the above range. We lose at most a further factor of $\mu N$ to ensure that every $j\in S$ satisfies $j\geq\mu N$.) The above along with Lemma D.2, implies that for every $j\in S$:

[TABLE]

The above with Corollary D.8 implies that for any $j\in S$ and large enough $N$:

[TABLE]

Thus, we have

[TABLE]

as desired. In the above, the second inequality follows from Lemma 5.7 and the fact that $j\geq\Omega(N)$ and the final inequality follows from (84) and the fact that $\left|S\right|\geq\Omega(N)$ (which in turn follows from (83) and our choice of parameters). The proof is complete. ∎

Appendix E A number theory problem

Let $N$ be some large enough integer. In this section we will consider sequences of reals: $0\leq y_{1}\leq y_{2}\leq\dots\leq`y_{N}<1$ that have certain nice proprieties when multiplied by integers.

We begin with some notation: given $z\in\mathbb{R}$, let $\left\{z\right\}$ denote the fractional part of $z$, i.e.

[TABLE]

We will also denote $\left\{z\right\}$ as $z\mod{1}$. We will also need the notation $\langle z\rangle$, which is its distance from the closest integer to $z$. Equivalently,

[TABLE]

We are now ready to define the kind of sequences we will encounter in our work:

Definition E.1**.**

A sequences of reals: $0\leq y_{1}\leq y_{2},\dots,y_{N}<1$ is called $s$-scattered (for any integer $s\geq 1$) if for any of the intervals $\left[\frac{i}{N},\frac{i+1}{N}\right)$ (with $i\in\mathbb{Z}_{N}$) has at most $s$ elements from the sequence in it.

We next define our notion of when a real number of good:

Definition E.2**.**

Let $0\leq\epsilon\leq 1$. We call a real $y\in[0,1)$ to be $\epsilon$-good if for any reals $0\leq\ell<r<1$, we have that

[TABLE]

If the above is not satisfied then we call $y$ to be $\epsilon$-bad.

We are interested in bounding how many bad $y_{i}$’s can be there in an $O(1)$-scattered sequence:

Lemma E.3**.**

Let $0\leq y_{1}\leq y_{2},\dots,y_{N}<1$ be an $O(1)$-scattered sequence. Define

[TABLE]

Then, we have

[TABLE]

E.1 Proof of Lemma E.3

We thank Stefan Steinerberger for showing us the following proof and kindly allowing us to use it.

We begin with Dirichlet’s approximation theorem, which implies the following result.

Lemma E.4**.**

For every $i\in\mathbb{Z}_{N}$, we have that there exists integers $0\leq p_{i}\leq q_{i}\leq N$ with $q_{i}\geq 1$ such that $\gcd(p_{i},q_{i})=1$ and

[TABLE]

The above immediately implies the following:

Corollary E.5**.**

Let $i\in\mathbb{Z}_{N}$ be such that $q_{i}\geq\left\lceil{\frac{4}{\epsilon}}\right\rceil$. Then we have for every $x\in\mathbb{Z}_{N}$

[TABLE]

The above basically says that we can essentially look at the goodness of rationals. In particular,

Lemma E.6**.**

Consider a rational $\frac{a}{b}$ with $1\leq a\leq b\leq N$ with $\gcd(a,b)=1$. Further let $b\geq\left\lceil{\frac{4}{\epsilon}}\right\rceil$. Then $\frac{a}{b}$ is $\frac{\epsilon}{2}$-good.

Proof.

Since $\gcd(a,b)=1$ (and $a\neq 0$), we have that the values $x\cdot\frac{a}{b}$ over all $x\in\mathbb{Z}_{N}$, take the values in $\frac{1}{b},\frac{2}{b},\dots\frac{b-1}{b},1$ between $\left\lfloor\frac{N}{b}\right\rfloor$ and $\left\lceil{\frac{N}{b}}\right\rceil$ times. It is not too hard to see that out of these $b$ values $(r-\ell)\cdot b\pm 2$ values can fall in the range $[\ell,r]$.171717Basically, in the worst case one can get two “extra” elements at the end points of $[\ell,r]$ for the upper bound. For the lower bound the closest two points to $\ell$ and $r$ might be just outside of $[\ell,r]$. This implies that

[TABLE]

as desired. ∎

Now Corollary E.5, Lemma E.6 and the triangle inequality implies that

Corollary E.7**.**

Let $i\in\mathbb{Z}_{N}$ be such that $q_{i}\geq\left\lceil{\frac{4}{\epsilon}}\right\rceil$. Then $y_{i}$ is $\frac{3\epsilon}{4}\leq\epsilon$-good.

We are now pretty much done. By the above result, we have that all $\epsilon$-bad $y_{i}$ have $q_{i}\leq\left\lceil{\frac{4}{\epsilon}}\right\rceil$. Further, since the sequence $\left\{y_{i}\right\}_{i\in\mathbb{Z}_{N}}$ is $O(1)$-scattered, Lemma E.4 also implies that any rational $\frac{a}{b}$ with $0\leq a\leq b\leq N$ and $\gcd(a,b)=1$, is “assigned” at most $O(1)$ many $y_{i}$.181818Indeed Lemma E.4 implies that for any $i\in\mathbb{Z}_{N}$, we have $\left|y_{i}-\frac{p_{i}}{q_{i}}\right|\leq\frac{1}{N}$ (since we must have $q_{i}\geq 1$). The number of rationals $\frac{p_{i}}{q_{i}}$ with $q_{i}\leq\left\lceil{\frac{4}{\epsilon}}\right\rceil$ is trivially at most $\left({\left\lceil{\frac{4}{\epsilon}}\right\rceil}\right)^{2}$. Thus, the overall number of $\epsilon$-bad $y_{i}$’s is at most $O\left({\frac{1}{\epsilon^{2}}}\right)$, as desired. Note that this argument implies the following result, which will be needed in our algorithms:

Corollary E.8**.**

Given a sequence $\left\{y_{i}\right\}_{i\in\mathbb{Z}_{N}}$ that is $O(1)$-scattered, all the $\epsilon$-bad $y_{i}$’s are contained in $O\left({\frac{1}{\epsilon^{2}}}\right)$ intervals, each of size at most $\frac{2}{N}$.