Algorithmic Foundations for the Diffraction Limit

Introduction

where $J_{1}$ is a Bessel function of the first kind. Under Feynman’s path integral formalism, $I(x)$ is precisely the pdf of the distribution over where the photon is detected (see Appendix A). The physical properties of the optical system, namely its numerical aperture and the wavelength of light being observed, determine $\sigma$ which governs the amount by which each point source gets blurred.

Intuitively, when point sources are closer together it seems harder to resolve them. However, despite considerable interest over the years [Abb73, Ray79, Sch04, Spa16, Hou27, Bux37], our understanding of what exactly can and cannot be resolved has never risen above heuristic arguments. In 1879 Lord Rayleigh [Ray79] proposed a criterion for assessing the resolving power of an optical system, which is still widely-used today, of which he wrote:

“This rule is convenient on account of its simplicity and it is sufficiently accurate in view of the necessary uncertainty as to what exactly is meant by resolution.”

Over the years, many researchers have proposed alternative criteria and offered arguments about why some are more appropriate than others. For example, in 1916 Carroll Sparrow proposed a new criterion [Spa16] that bears his name, which he justified as follows:

“It is obvious that the undulation condition should set an upper limit to the resolving power $\dots$ The effect is observable both in positives and in negatives, as well as by direct vision $\dots$ My own observations on this point have been checked by a number of my friends and colleagues.”

Even more resolution criteria were proposed, both before and after, by Ernst Abbe [Abb73], Sir Arthur Schuster [Sch04], William Houston [Hou27], etc. Their popularity varies depending on the application area and research community. Many researchers have also pushed back on the idea that there is a fundamental diffraction limit at all. In his 1964 Lectures on Physics [FLS11, Section 30-4], Richard Feynman writes:

“$\dots$ it seems a little pedantic to put such precision into the resolving power formula. This is because Rayleigh’s criterion is a rough idea in the first place. It tells you where it begins to get very hard to tell whether the image was made by one or by two stars. Actually, if sufficiently careful measurements of the exact intensity distribution over the diffracted image spot can be made, the fact that two sources make the spot can be proved even if $\theta$ is less than $\lambda/L$.”

“Mathematics cannot set any lower limit for the distance of two resolvable points.”

Our goal in this work is to remedy this gap in the literature and place the notion of the diffraction limit on rigorous statistical foundations by drawing new connections to recent work in theoretical computer science on provably learning mixture models, as we will describe next. First we remark that the way the diffraction limit is traditionally studied is in fact a mixture model. In particular we assume that, experimentally, we can measure photons that are sampled from the true diffracted image. However we only observe a finite number of them because our experiment has finite exposure time, and indeed as we will see in some settings the number of samples needed to resolve closely-spaced objects can explode and be essentially impossible just from statistical considerations. Moreover we may only be able to record the location of observed photons up to some finite accuracy, which can also be thought of as being related to sampling error. The main question we will be interested in is:

How many samples (i.e. photons) are needed to accurately estimate the centers and relative intensities of a mixture (i.e. superposition) of two or more Airy disks, as a function of their separation and the parameters of the optical system?

This is a central question in optics. Fortunately, there are many parallels between this question and the problem of provably learning mixture models that surprisingly seem to have gone undiscovered. In particular, let us revisit Sparrow’s argument that resolution is impossible when the density function becomes unimodal. In fact there are already counter-examples to this claim, albeit not for mixtures of Airy disks. It is known that there are algorithms for learning mixtures of two Gaussians that take a polynomial number of samples and run in polynomial time. These algorithms work even when the density function is unimodal, and require just that the overlap between the components can be bounded away from one. Moreover when there are $k$ components it is known that there is a critical separation above which it is possible to learn the parameters accurately with a polynomial number of samples, and below which accurate learning requires a superpolynomial number of samples information-theoretically [RV17]. Thus a natural way to formulate what the diffraction limit is, so that it can be studied rigorously, is to ask:

At what critical separation does the sample complexity of learning mixtures of $k$ Airy disks go from polynomial to exponential?

In this work we will give algorithms whose running time and sample complexity are polynomial in $k$ above some critical separation, and prove that below some other critical separation the sample complexity is necessarily exponential in $k$. These bounds will be within a universal constant, and thus we approximately locate the true diffraction limit. There will also be some surprises along the way, such as the fact that the Abbe limit, which has long been postulated to be the true diffraction limit, is not actually the correct answer!

Before we proceed, we also want to emphasize that there is an important conceptual message in our work. First, for mixtures of Gaussians the model was only ever supposed to be an approximation to the true data generating process. For example, Karl Pearson introduced mixtures of Gaussians in order to model various physical measurements of the Naples crabs. However mixtures of Gaussians always have some chance of producing samples with negative values, but Naples crabs certainly do not have negative forehead lengths! In contrast, for mixtures of Airy disks the model is an extremely accurate approximation to the observations in many experimental setups because it comes from first principles. It is particularly accurate in astronomy where for all intents and purposes the lens is spherical and the star is so far away that it is a point source, and the question itself is highly relevant because it arises when we want to locate double-stars [Fal67] .

Furthermore we believe that there ought to be many more examples of inverse problems in science and engineering where tools and ideas from the literature on provably learning mixture models ought to be useful. Indeed both mixtures of Gaussians and mixtures of Airy disks can be thought of as inverse problems with respect to simple differential equations, for the heat equation and a modified Bessel equation respectively. While this is a well-studied topic in applied mathematics, usually one makes some sort of smoothness assumption on the initial data. What is crucial to both the literature on learning mixtures of Gaussians and our work is that we have a parametric assumption that there are few components. Thus we ask: Are there provable algorithms for other inverse problems, coming from differential equations, under parametric assumptions? Even better: Could techniques inspired by the method of moments play a key role in such a development?

It is often the case that heuristic arguments, despite being quite far from a rigorous proof, predict the correct thresholds for a wide range of statistical problems. However here there will be a surprise. In a seminal work in 1873, Ernest Abbe formulated what is now called the Abble limit. Since then it has been widely accepted in the optics literature as the critical distance below which diffraction makes resolution impossible for classical optical systems. In the mixture model formalism outlined above, it corresponds to a separation of $\pi\sigma$ between any pair of Airy disk centers $\bm{\mu}_{i},\bm{\mu}_{j}$. This distance arises naturally because it corresponds to the radius of the support of the Fourier transform of the Airy disk kernel $A_{\sigma}:\bm{x}\mapsto\frac{1}{\pi\sigma^{2}}\left(\frac{J_{1}(\bm{x}/\sigma)}{\bm{x}/\sigma}\right)^{2}$ (see Appendix B.4 for further discussion).

One of the main results of this work is to show that resolution is statistically hard even above the Abbe limit! Specifically, we show that even for mixtures of Airy disks whose centers have a pairwise separation that is a constant factor larger than the Abbe limit, the problem of recovering their locations can require $\exp(\Omega(\sqrt{k}))$ samples. The main challenge is that no configuration where the Airy disk centers are all on the same line can beat the Abbe limit. Instead we construct a new, natural lower bound instance.

Let $\underline{\gamma}\triangleq\sqrt{4/3}\approx 1.155$. For any $0<\epsilon<1$, there exist two superpositions of $k$ Airy disks $\rho,\rho^{\prime}$ which are both $\underline{\gamma}\cdot(1-\epsilon)\cdot\pi\sigma$-separated and such that 1) $\rho$ and $\rho^{\prime}$ have noticably different sets of centers, and yet 2) it would take at least $\exp(\Omega(\epsilon\sqrt{k}))$ samples to distinguish whether the samples came from $\rho$ or from $\rho^{\prime}$.

On the other hand, we also show that when the Airy disks have separation that is a small constant factor larger than this critical distance, there is an algorithm for recovering the centers that takes a polynomial number of samples and runs in polynomial time.

Define the absolute constant $\overline{\gamma}=\frac{2j_{0,1}}{\pi}=1.530\ldots$, where $j_{0,1}$ is the first positive zero of the Bessel function $J_{0}$. Let $\rho$ be a $\overline{\gamma}\cdot\pi\sigma$-separated superposition of $k$ Airy disks where every disk has relative intensity at least $\lambda$. Then for any target error $\epsilon>0$, there is an algorithm with time and sample complexity $N={\mathsf{poly}}\left(k,1/\Delta,1/\lambda,1/\epsilon\right)$ which outputs an estimate for the centers and relative intensities of $\rho$ which incurs error $\epsilon$ with probability at least $9/10$. Furthermore, this holds even when there is granularity in the photon detector, as long as it is at most some inverse polynomial in all parameters.

The main open question of our work is to prove matching upper and lower bounds that pin down the true diffraction limit. However, as we will discuss, this is a challenging problem in harmonic analysis, despite being connected to areas where there has been considerable recent progress. Moreover this phase transition for resolution is actually more dramatic than what happens for mixtures of Gaussians [RV17]. Even ignoring the issue of computational complexity, for spherical Gaussian mixtures it is known that at separation $o(\sqrt{\log k})$, super-polynomially many samples are needed, while at separation $\Omega(\sqrt{\log k})$, polynomially many suffice.

We now say a word about the techniques that go into proving Theorem 1.1 and Theorem 1.2. It turns out that both are closely related to proving a modified version of an Ingham-type estimate [KL05]:

What is the smallest $\Delta$ for which the quantity

In particular, the main technical step for showing Theorem 1.2 is to show that the critical $\Delta$ in Question 1.1 is at most $2j_{0,1}/\pi$. This can be obtained via the following extremal function. A ball minorant, is a function $F$ satisfying the properties that

$F(x)\leq\mathop{\mathds{1}}[x\in B]$ and

$\widehat{F}$ is supported on the ball of radius $\Delta$

In [HV+96, CCLM17, Gon18] it was shown that such a ball minorant exists for $\Delta=2j_{0,1}/\pi$ (interestingly, this paved the way to some recent progress on Montgomery’s famous pair correlation conjecture for the Riemann zeta function [CCLM17]). One can use property $(1)$ to pass from integrating against the function $\mathop{\mathds{1}}[x\in B]$ to integrating against $F$. And because by property $(2)$ $F$ is localized in the frequency domain, the latter integral is large. In fact the one-dimensional analogue of Question 1.1 was resolved in [Moi15] using the univariate analogue of $F$, namely the Beurling-Selberg minorant. However the algorithmic approach only made sense in one-dimension. In our case, we employ the tensor generalization of the matrix pencil method, originally introduced in [HK15]. We defer the details of this to Section 4.3.

For the lower bound in Theorem 1.1, we need to answer a variant of Question 1.1.

The connection to Theorem 1.1 is straightforward: By Plancherel’s and smoothness properties of $A_{\sigma}$, one can upper bound the $L_{1}$ distance between the mixture of Airy disks given by parameters $\{\lambda_{j}\},\{\bm{\mu}_{j}\}$ and the mixture given by $\{\lambda^{\prime}_{j}\},\{\bm{\mu}^{\prime}_{j}\}$ in terms of the left-hand side of (3). So if one can construct a set of $\Delta$-separated centers $\{\bm{\mu}_{j}\},\{\bm{\mu}^{\prime}_{j}\}$ for which (3) fails to hold but for which the collection $\{\bm{\mu}_{j}\}$ is separated from $\{\bm{\mu}^{\prime}_{j}\}$ but the resulting mixtures of Airy disks are $o(1/{\mathsf{poly}}(k))$-close in total variation distance it implies that resolution is statistically impossible with a polynomial number of samples. This is the recipe used in known lower bounds [MV10, HP15, RV17] for learning mixtures of Gaussians.

For our purposes, it turns out that “tensoring” one-dimensional lower bounds does not work because it would not beat the Abbe limit [Moi15]. Morally, this is because tensoring the unit interval with itself would give us the unit square, which corresponds to separation in the $L_{\infty}$ distance rather than the $L_{2}$ distance, and the $L_{2}$ distance is the right distance in optics because it is rotationally invariant. The main technical contribution in our lower bound is to give a more sophisticated construction given by interleaving two triangular lattices and placing the centers at points on these lattices (see Figure 5). The analysis is rather delicate, and we defer the details to Section 2 and Section 5.

To complete the picture, we show that there is no diffraction limit when the number of Airy disks is a constant. In particular we show that for any constant number of Airy disks there is an algorithm that takes a polynomial number of samples and runs in polynomial time that learns the parameters to any desired accuracy regardless of the separation.

Let $\rho$ be a $\Delta$-separated superposition of $k$ Airy disks where every disk has relative intensity at least $\lambda$. Then for any target error $\epsilon>0$ and failure probability $\delta>0$, there is an algorithm which draws $N={\mathsf{poly}}\left((k\sigma/\Delta)^{k^{2}},1/\lambda,1/\epsilon,\log(1/\delta)\right)$ samples from $\rho$, runs in time $O(N)$, and outputs an estimate for the centers and relative intensities of $\rho$ which incurs error $\epsilon$ with probability at least $1-\delta$. Furthermore, this holds even when there is granularity in the photon detector, as long as it is at most some inverse polynomial in $N$.

This result turns out to be simple in retrospect, and comes from assembling a few standard tools from the literature on provably learning mixture models. Nevertheless it underscores an important point that existing tools can already have important implications for inverse problems the sciences. Our approach is to first estimate the Fourier transform $\widehat{\rho}$ from samples and then pointwise divide by $\widehat{A}_{\sigma}$. In this way we can simulate noisy access to the Fourier transform of the mixture of delta functions at $\bm{\mu}_{1},\ldots,\bm{\mu}_{k}$. However $\widehat{A}_{\sigma}$ has compact support, so we can only access frequencies with bounded $L_{2}$ norm. Now we can reduce to the one-dimensional case [Moi15] by projecting $\rho$ along two nearby directions, solving each resulting univariate problem, and then solving an appropriate linear system to recover the centers and relative intensities. This method is reminiscent of [KMV10, MV10], which gives algorithms for learning high-dimensional mixtures of Gaussians based on reducing to a series of one-dimensional problems and stitching together these estimates carefully.

2 Related Work

We have already mentioned that our work is closely related to the vast literature on learning mixture models and, in particular, on learning mixtures of Gaussians [Das99, DS00, AK01, VW02, AM05, BV08, KMV10, MV10, BS15, HP15, HK13, GHK15, RV17, HL18, KS17, DKS18]. Here we mention some other connections to work on recovering spike trains from noisy, band-limited Fourier measurements.

The seminal work of [DS89, Don92] was one of the first to put this question on rigorous footing. Donoho studied the modulus of continuity for this problem on a grid as the grid width goes to zero. Later Candes and Fernandez-Granda [CFG14] gave a practical algorithm based on $L_{1}$ minimization over a continuous domain. There has been a long line of work on this problem which it would also be impossible to survey fully, so we refer the reader to [CFG13, TBSR13, FG13, Lia15, Moi15, FG16, KPRvdO16, MC16] and references therein.

We remark that essentially all works on super-resolution in high dimensions focus on the case where measurements are $L_{\infty}$-band-limited rather than $L_{2}$-band-limited. Given the prevalence of Airy disks and circular apertures in statistical optics, one upshot of our work is that, technical issues related to the so-called box (aka $L_{\infty}$ ball) minorant problem notwithstanding, the $L_{2}$ setting may be the more practically relevant one to consider anyways.

There are also connections to the extensive literature on the sparse Fourier transform, which can be interpreted in some sense as the “agnostic” version of the super-resolution problem where the goal is to compete with the error of the best $k$-sparse approximation to the discrete Fourier transform, even in the presence of noise, using few measurements [GGI+02, GMS05, HIKP12, GIIS14, IKP14, Kap16]. When the $k$ spikes need not be at discrete locations and the low-frequency measurements are randomly chosen, this is the problem of compressed sensing off the grid introduced by [TBSR13], for which recovery is possible with far fewer measurements. This can be thought of as the one-dimensional case of the setting of [HK15]. To our knowledge, the only work that addresses the continuous, high-dimensional version of the sparse Fourier transform is the very recent work of [JLS20]. The emphasis in this literature is primarily on obtaining sample complexity near-linear in $k$, whereas our guarantees are only polynomial in $k$. Consequently the results in the sparse Fourier transform literature lose log factors in the level of separation they require, whereas in our setting the emphasis is primarily on the level of separation needed to get polynomial-time and -sample algorithms.

3 Visualizing the Diffraction Limit

In this short section we provide some figures to help conceptualize our results. Figure 2 illustrates the basic notion that separation is information-theoretically unnecessary for parameter learning of superpositions of Airy disks. We compare the discretized empirical distribution of samples from two diffraction patterns whose components have separation well below the diffraction limit and thus well below what conventional wisdom in optics suggests is resolvable. While the differences in the diffraction patterns are minute, they do indeed become statistically significant with enough samples. Eventually it becomes possible to conclude that the gray diffraction pattern is generated by one point source and the red diffraction pattern is generated by two.

Next, we present a striking visual representation of the statistical barrier imposed by the diffraction limit when the number of components is large. Recall that the upshot of Theorems 1.1 and 1.2 is that $k$ plays a leading role in determining when resolution is and is not feasible: slightly above the Abbe limit, the sample (and computational) complexity is polynomial in $k$, and anywhere beneath the Abbe limit, the sample complexity becomes exponential in $k$. This helps clarify why in some domains like astronomy, where there are only ever a few tightly spaced point sources, there is evidently no diffraction limit. Yet in other domains like microscopy where there are a large number of tightly spaced objects, the diffraction limit is indeed a fundamental barrier, at least in the classical physical setup. This helps explain why different communities have settled on different beliefs about whether there is or is not a diffraction limit.

In Figure 3 we experimentally investigate this phenomenon and illustrate how the total variation distance scales as we vary the number of disks and the separation in our earlier constructions. It is evident from these plots that for any superposition of a few Airy disks, there is no sharp dividing line between what is and is not possible to resolve. But when the number of Airy disks becomes large, with any reasonable number of samples, it is feasible to resolve the superposition if and only if their separation is at least as large as the Abbe limit.

4 Roadmap

In Section 2 we give a preview of our lower bound proof by providing a self-contained answer to Question 1.2. In Section 3, we give an overview of our probabilistic model, some notation, and other mathematical preliminaries. In Section 4, we prove the algorithmic results in Theorems 1.3 and 1.2. In Section 5 we complete the proof of our lower bound from Theorem 1.1. In Section 6 we conclude with some directions for future work. In Appendix A, we overview previous attempts in the optics literature to put the diffraction limit on rigorous footing. In Appendix B, we describe and motivate our model and also define the various resolution criteria which have appeared in the literature. In Appendix C, we catalogue quotations from the literature that are representative of the points of view addressed in the introduction. In Appendix D, we complete some deferred proofs. Lastly, in Appendix E, we give details on how Figure 3 was generated.

Lower Bound Preview

In this section we give a self-contained proof of one of the main technical ingredients in the proof of our main result, Theorem 1.1. Before proceeding, it will be convenient to introduce a bit of notation; any outstanding notation we will present Section 3, e.g. our convention for the Fourier transform. Recalling that $\underline{\gamma}\triangleq\sqrt{4/3}$, define

for any small constant $\epsilon>0$ so that the critical level of separation for which Theorem 1.1 applies is $\Delta\triangleq 2/m=\underline{\gamma}\cdot(1-\epsilon)\cdot\pi\sigma$.This $\pi\sigma$ scaling is not important to the result in this section but is the natural choice of scaling for Airy disks, so it will be convenient to work with this when we apply the results of this section to prove Theorem 1.1. Additionally, let $k$ be an odd square and define

This construction is illustrated in Figure 5: there, similarly colored points correspond to centers in the same mixture, and our choice of $\{\nu_{j_{1},j_{2}}\}$ ensures that the level of separation between any two points in a particular mixture is $\Delta$, which is slightly less than $\sqrt{4/3}$ times the Abbe limit of $\pi\sigma$. As such, the following tells us that the answer to Question 1.2 is surprisingly at least $\sqrt{4/3}$, rather than 1 as the Abbe limit would suggest:

for all $\lVert x\rVert\leq 1/\pi\sigma$. Furthermore, $\mathop{\text{sgn}}(u_{j_{1},j_{2}})=(-1)^{j_{1}+j_{2}}$, and

We need the following ingredient from the proof of the one-dimensional lower bound in [Moi15].

where $h_{j_{1},j_{2}}=\alpha^{\prime}_{j_{1}}\alpha^{\prime}_{j_{2}}\geq 0$, where $\alpha^{\prime}_{j}\triangleq\alpha_{j}\cdot\mathop{\mathds{1}}[j=0]+m\alpha_{j}\cdot\mathop{\mathds{1}}[j\neq 0]$. We will take

Observe that $\mathop{\text{sgn}}(u_{j_{1},j_{2}})=(-1)^{j_{1}+j_{2}}$ as desired.

By taking the inverse Fourier transform of $\hat{H}$, we get that

To complete our proof, it therefore suffices to show that $H(\bm{x})\leq\exp(-\Omega(\epsilon\sqrt{k}))$ for all $\lVert x\rVert\leq 1/\pi\sigma$.

We claim that $S$ contains the origin-centered $L_{\infty}$ ball of radius $\epsilon/2\sqrt{2}$. Note that $S$ is given by translating the four connected components of $(\times)\backslash B_{1}$, which is nonempty because $B_{1}$ consists of points $(x_{1},x_{2})$ satisfying

In particular, for $x_{1},x_{2}\in$ satisfying $|x_{1}-1/2|,|x_{2}-1/2|>(1-\epsilon)/2$, observe that the left-hand quantity in (14) satisfies

where the last step follows by our choice of $\underline{\gamma}=\sqrt{4/3}$. We conclude that $S$ contains the origin-centered $L_{\infty}$ ball of radius $\epsilon/2$ as claimed.

The last step is just to scale $u$ so that (7) holds. First note that by substituting $x=0$ into (13), we have that

Preliminaries

In this section we explain the terminology and notation that we will adopt in this work and also provide some technical preliminaries that will be useful later.

We first formally define the family of distributions we study in this work.

Note that the factor of $\frac{1}{\pi\sigma^{2}}$ in the definition of $A_{\sigma}$ is to ensure that $A_{\sigma}(\cdot)$ is a probability density.

It will be straightforward to extend the above model to take into account error stemming from the fact that the photon detector itself only has finite precision.

Given discretization parameter $\varsigma>0$, we say $\bm{x}$ is a $\varsigma$-granular sample from $\rho$ if it is produced via the following generative process: 1) a point $\bm{x}^{\prime}$ is sampled from $\rho$, 2) $\bm{x}$ is obtained by moving $\bm{x}^{\prime}$ an arbitrary distance of at most $\varsigma$.

$\hat{A}_{\sigma}[\omega]=\frac{2}{\pi}(\arccos(\pi\sigma\lVert\omega\rVert)-\pi\sigma\lVert\omega\rVert\sqrt{1-\pi^{2}\sigma^{2}\lVert\omega\rVert^{2}}$.

It is enough to show this for $\sigma=1$. Let $G(\bm{x})\triangleq J_{1}(\lVert\bm{x}\rVert)/\lVert\bm{x}\rVert$. It is a standard fact that the zeroeth-order Hankel transform of the function $r\mapsto J_{1}(r)/r$ is the indicator function of the interval $$. Using our convention for the Fourier transform (see (18)), this implies that$\hat{G}[\omega]=2\pi\cdot\mathop{\mathds{1}}[\lVert\omega\rVert\in[0,1/2\pi]]$. Because$A_{1}=G^{2}/\pi$, by the convolution theorem we conclude that$\hat{A}_{1}$is$\frac{1}{\pi}$times the convolution of$\hat{G}$with itself, which is just$4\pi^{2}$times the convolution of the indicator function of the unit disk of radius$1/2\pi$with itself. By elementary Euclidean geometry one can compute this latter function to be$\omega\mapsto\frac{1}{2\pi^{2}}\cdot\left(\arccos(\pi\lVert\omega\rVert)-\pi\lVert\omega\rVert\sqrt{1-\pi^{2}\lVert\omega\rVert^{2}}\right)$, from which the claim follows. ∎

In optics, the two-dimensional Fourier transform of the point-spread function is called the optical transfer function, a term we will occasionally use in the sequel.

Now note that by Fact 2, $\hat{A}_{\sigma}$ is supported only over the disk of radius $\frac{1}{\pi\sigma}$ centered at the origin in the frequency domain. In the spatial domain, this corresponds to a separation of $\pi\sigma$; this is the definition of the Abbe limit. We will need the following elementary estimate for $\widehat{A}[\omega]$:

For all $\lVert\omega\rVert_{2}\leq 1$, $\widehat{A}[\omega]\geq(1-\lVert\omega\rVert_{2})^{2}$.

As the algorithms we give will be scale-invariant, we will assume that $\sigma=1/\pi$ in the rest of this work and refer to $A_{1/\pi}$ as $A$.

The following terminology formalizes what it means for an algorithm to return an accurate estimate for the parameters of a superposition of Airy disks.

$\left(\{\lambda^{*}_{i}\}_{i\in[k]},\{\bm{\mu}^{*}_{i}\}_{i\in[k]}\right)$is an $(\epsilon_{1},\epsilon_{2})$-accurate estimate for the parameters of a superposition of $k$ Airy disks $\rho$ with centers $\{\bm{\mu}_{i}\}$ and relative intensities $\{\lambda_{i}\}$ if there exists a permutation $\tau$ for which

Given matrices $M,N$, we will denote by $(M,N)$ the generalized eigenvalue problem $Mx=\lambda Nx$. In any solution $(\lambda,x)$ to this, $\lambda$ is called a generalized eigenvalue and $x$ is called a generalized eigenvector.

In Section 5, we need the following estimate for $J_{\nu}(z)$:

for all $(i_{1},i_{2},i_{3})\in[m^{\prime}_{1}]\times[m^{\prime}_{2}]\times[m^{\prime}_{3}]$.

Learning Superpositions of Airy Disks

In this section we present the technical details of our algorithmic results. In Sections 4.2 and 4.4, we prove the following formal version of Theorem 1.3.

Let $\rho$ be a $\Delta$-separated superposition of $k$ Airy disks with minimum mixing weight $\lambda_{\min}$ and such that $\lVert\bm{\mu}_{i}\rVert\leq\mathcal{R}$ for all $i\in[k]$.

For any $\epsilon_{1},\epsilon_{2}>0$, there is some $\alpha={\mathsf{poly}}\left(\log 1/\delta,1/\lambda_{\min},1/\epsilon_{1},1/\epsilon_{2},\mathcal{R},(k\sigma/\Delta)^{k^{2}}\right)^{-1}$ for which there exists an algorithm with time and sample complexity ${\mathsf{poly}}(1/\alpha)$ which, given $\varsigma={\mathsf{poly}}(\alpha)$-granular sample access to $\rho$, outputs an $(\epsilon_{1},\epsilon_{2})$-accurate estimate for the parameters of $\rho$ with probability at least $1-\delta$.

Specifically, in Section 4.2, we show how one can use the matrix pencil method to recover the parameters for $\rho$ given oracle access to the optical transfer function, i.e. the two-dimensional Fourier transform of $\rho$, up to some small additive error. In Section 4.4, we show how to implement this approximate oracle.

In Section 4.3, we also use the oracle of Section 4.4 to prove the following formal version of Theorem 1.2.

Let $\rho$ be a $\Delta$-separated superposition of $k$ Airy disks with minimum mixing weight $\lambda_{\min}$ and such that $\lVert\bm{\mu}_{i}\rVert\leq\mathcal{R}$ for all $i\in[k]$. Let

where $j_{0,1}$ is the first positive zero of the Bessel function of the first kind $J_{0}$. For any $\Delta>\overline{\gamma}\cdot\pi\cdot\sigma$, the following holds:

For any $\epsilon_{1},\epsilon_{2}>0$, there is some $\alpha=1/{\mathsf{poly}}\left(k,\mathcal{R},\sigma/\Delta,1/\lambda_{\min},1/\epsilon_{1},1/\epsilon_{2},1/(\Delta-\overline{\gamma})\right)$ for which there exists an algorithm with time and sample complexity ${\mathsf{poly}}(1/\alpha)$ which, given ${\mathsf{poly}}(\alpha)$-granular sample access to $\rho$, outputs an $(\epsilon_{1},\epsilon_{2})$-accurate estimate for the parameters of $\rho$ with probability at least $4/5$.

In this section we reduce the problem of learning superpositions of Airy disks to the problem of learning a convex combination of Dirac deltas given the ability to make noisy, band-limited Fourier measurements.

Formally, suppose we had access to the following oracle:

As we will see in Section 4.4, $\mathcal{O}$ will be constructed by sampling some number of points from $\rho$ and computing empirical averages. The number $m$ and accuracy $\eta$ of queries that $\mathcal{O}$ can answer dictates the sample complexity of this procedure. As we will see in the proofs of Lemma 4.6 and Lemma 4.12 below, the $m$ that we need to take will be small, so the reader can ignore $m$ and pretend it is unbounded for most of this section.

where for $\omega=(r\cos\theta,r\sin\theta)$, we have by Fact 2 that

In particular, $\widehat{A}[\omega]$ only depends on $r=\lVert\omega\rVert$ (because $A(\cdot)$ is radially symmetric), so henceforth regard $\widehat{A}$ as a function merely of $r$.

This is a trigonometric polynomial to which we have noisy pointwise access using $\mathcal{O}$:

Let $0<r<1$. With an $\eta$-approximate OTF oracle $\mathcal{O}$, on input $\omega\in B^{2}(r)$ we can produce an estimate of $F(\omega)$ to within $\eta/\widehat{A}[r]$ additive error.

By dividing by $\widehat{A}[\omega]$ on both sides of (22), we get that

where the last step uses the fact that $\widehat{A}[\cdot]$ is decreasing on the interval $$. ∎

So concretely, given an $\eta$-approximate OTF oracle, we have reduced the problem of learning superpositions of Airy disks to that of recovering the locations of $\{\bm{\mu}_{j}\}$ given the ability to query $F(\omega)$ at arbitrary frequencies $\omega$ for which $\lVert\omega\rVert_{2}<1$ to witin additive accuracy $\eta/\widehat{A}[\lVert\omega\rVert_{2}]$.

Lastly, for reasons that will become clear in subsequent sections (see e.g. (64)), it will be convenient to assume that $\mathcal{R}\leq 1/3$. This is without loss of generality, as otherwise, we can scale the data down by a factor of $3\mathcal{R}$ so that they are now i.i.d. samples from the superposition of Airy disks with density $\rho^{\prime}(\bm{x})\triangleq\sum^{k}_{j=1}\lambda_{j}\cdot A_{1/\mathcal{R}}\left(\bm{x}-\bm{\mu}_{j}/\mathcal{R}\right)$. Define the rescaled centers $\bm{\mu}^{\prime}_{j}\triangleq\bm{\mu}_{j}/\mathcal{R}$ and note that by assumption, $\|\bm{\mu}^{\prime}_{j}\|_{2}\leq 1/3$ for all $j\in[k]$.

The Fourier transform of $\rho^{\prime}$ is then given by $\hat{\rho}^{\prime}(\omega)=\hat{A}_{1/\mathcal{R}}[\omega]\sum^{k}_{j=1}\lambda_{j}e^{-2\pi i\langle\bm{\mu}^{\prime}_{j},\omega\rangle}$, so by the proof of Lemma 4.1 we conclude that with an $\eta$-approximate OTF oracle for $\rho$, for any $0<r<1$ on input $\omega\in B^{2}(r\cdot\mathcal{R})$ we can produce an estimate of $\sum^{k}_{j=1}\lambda_{j}e^{-2\pi i\langle\bm{\mu}^{\prime}_{j},\omega\rangle}$ to within $\eta/\hat{A}[r]$ additive error. Recovering the centers $\{\bm{\mu}^{\prime}_{j}\}$ to within additive error $\epsilon$ then translates to recovering the centers $\{\bm{\mu}_{j}\}$ to within additive error $3\mathcal{R}\epsilon$. For this reason, we will henceforth assume that $\mathcal{R}\leq 1/3$.

2 Learning via the Optical Transfer Function

By the discussion at the end of Section 4.1, we may assume $\|\bm{\mu}_{i}\|_{2}\leq 1/2$ for all $i\in[k]$, so $\lVert\bm{\mu}_{i}-\bm{\mu}_{j}\rVert_{2}\leq 1$ for all $i\neq j$. For $j\in[k]$, let $m_{j}=\langle\bm{\mu}_{j},v\rangle$ and $\alpha_{j}=e^{2\pi i\cdot(m_{j}/4k)}$. In this section we will assume that $m_{j}\neq 0$ for all $j\in[k]$

Consider the generalized eigenvalue problem $(VD_{\lambda}V^{\top},VD_{\lambda}D_{\alpha}V^{\top})$ where

The following standard facts are key to the matrix pencil method:

The generalized eigenvalues of $(VD_{\lambda}V^{\top},VD_{\lambda}D_{\alpha}V^{\top})$ are exactly $\alpha_{1},...,\alpha_{k}$.

so we must instead work with the generalized eigenvalue problem $(VD_{\lambda}V^{\top}+E,VD_{\lambda}D_{\alpha}V^{\top}+F)$, where the $(i,j)$-th entry of $E$ (resp. $F$) is the noise $\eta^{\prime}_{i+j-2}$ (resp. $\eta^{\prime}_{i+j-1}$) in the observation of $v_{i+j-2}$ (resp. $v_{i+j-1}$).

If $V$ is well-conditioned, one can apply standard perturbation bounds to argue that the solutions to this generalized eigenvalue problem are close to those of the original $(VD_{\lambda}V^{\top},VD_{\lambda}D_{\alpha}V^{\top})$. Moreover, given approximations $\widehat{\alpha}_{1},...,\widehat{\alpha}_{k}$ to these generalized eigenvalues, we can find approximations $\widehat{\lambda}_{1},...,\widehat{\lambda}_{k}$ to $\lambda_{1},...,\lambda_{k}$ by solving the system of equations $\bm{v}=\widehat{V}\bm{\lambda}$, where $\bm{v}=(v_{0},...,v_{k-1})$, $\bm{\lambda}=(\widehat{\lambda}_{1},...,\widehat{\lambda}_{k})$, and

The formal specification of the matrix pencil method algorithm ModifiedMPM that we use is given in Algorithm 1.

The following theorem, implicit in the proof of Theorem 2.8 in [Moi15], makes the above reasoning precise. Henceforth, let $\kappa(\Delta^{\prime})$ and $\sigma_{\min}(\Delta^{\prime})$ respectively denote the condition number and minimum singular value of $V$ when $\frac{m_{i}}{4k},\frac{m_{j}}{4k}$ have minimum separation $\Delta^{\prime}$ for all $i\neq j$, and define $\lambda_{\min}=\min_{i}\lambda_{i}$, $\lambda_{\max}=\max_{i}\lambda_{i}$.

Then if $\lVert E\rVert+\lVert F\rVert<\sigma_{\min}(\Delta^{\prime})^{2}\lambda_{\min}$ and $\gamma<\Delta^{\prime}/4$, ModifiedMPM produces estimates $\{\widehat{\lambda}_{i}\}$ for the mixing weights and estimates $\{\widehat{m}_{i}\}$ for the projected centers such that for some permutation $\tau$:

Note that the guarantees of Theorem 33 are stated in [Moi15] in terms of wraparound distance on the interval $[-1/2,1/2]$, but because $\frac{m_{i}}{4k}\in[-1/4,1/4]$ for all $j\in[k]$, $\frac{m_{1}}{4k},....,\frac{m_{k}}{4k}$ have pairwise separation $\Delta^{\prime}$ both in absolute and wraparound distance.

In other words, the output of ModifiedMPM converges to the true values for $\{\langle\bm{\mu}_{1},v\rangle\}_{j\in[k]}$ and $\{\lambda_{j}\}_{j\in[k]}$ at a rate polynomial in the noise rate, condition number of $V$, and relative intensity of the Airy disks, provided $\sigma_{\min}(\Delta^{\prime})$ is inverse polynomially large and $\kappa(\Delta^{\prime})$ is polynomially small in those parameters.

To complete the argument, we must establish these bounds on $\sigma_{\min}$ and $\kappa$. Henceforth, let

where in the first step we used the standard fact that the absolute value of the determinant of a square matrix is equal to the product of its singular values, in the second step we used the standard identity for the determinant of a Vandermonde matrix, in the third step we used the angular separation of the $\alpha_{i}$’s, and in the final step we used the elementary inequality $\cos(\Delta^{\prime})\leq 1-\Delta^{\prime 2}/2$. We may thus naively lower bound $\sigma_{\min}(V)$ by $\frac{\Delta^{\prime k(k-1)/2}}{k^{k-1}}$, from which the lemma follows. ∎

This yields the following consequence for ModifiedMPM.

then ModifiedMPM produces estimates $\{\widehat{\lambda}_{i}\}$ for the mixing weights and estimates $\{\widehat{m}_{i}\}$ for the projected centers such that for some permutation $\tau$:

From Lemma 35, we have that $\sigma_{\min}(\Delta^{\prime})^{2}\geq(\Delta^{\prime k}/k^{2})^{k-1}$. Then because $\kappa(\Delta^{\prime})^{2}\leq k^{2}/\sigma_{\min}(\Delta^{\prime})^{2}$, we would like to conclude by Theorem 33 that $|m_{i}-\widehat{m}_{\tau(i)}|\leq 8\gamma$, where

where in the last step we use that the vector $\eta^{\prime}$ has length $O(k)$ and satisfies $\lVert\eta^{\prime}\rVert_{\infty}\leq O(\eta)$ by (28). To do so, we just need to verify that $\lVert E\rVert+\lVert F\rVert<\sigma_{\min}(\Delta^{\prime})^{2}\lambda_{\min}$ and $\gamma<\Delta^{\prime}/4$. The latter clearly follows from the bound (37) for sufficiently small $\Cr{etagap}$. For the former, note that

Finally, Theorem 33 also implies that $|\lambda_{i}-\widehat{\lambda}_{i}|\leq\zeta$, where

2.2 Combining Directions

We first show a “pairing lemma” stating that if $\omega_{1}$ is chosen randomly and $\omega_{2}$ is chosen to be close to $\omega_{1}$, then if one sorts the centers $\bm{\mu}_{1},...,\bm{\mu}_{k}$, first in terms of their projections in the $\omega_{1}$ direction, and then in terms of their projections in the $\omega_{2}$ direction, the corresponding elements in these two sorted sequences will correspond to the same centers.

We require the following elementary fact.

Then with probability at least $1-\frac{k(k-1)\theta}{\pi}$, for every $i\neq j$ the following are equivalent: I) $m_{i}>m_{j}$, II) $m^{\prime}_{i}>m^{\prime}_{j}$, III) $\hat{m}_{i}>\hat{m}_{j}$, and IV) $\hat{m}^{\prime}_{i}>\hat{m}^{\prime}_{j}$.

By Lemma 4.3 and a union bound we have that with probability $1-\frac{k(k-1)\theta}{\pi}$, $|m_{i}-m_{j}|>\Delta\sin\theta$ for all $i\neq j$. Fix any $i\neq j$ and suppose that $m_{i}>m_{j}$. Then by triangle inequality and Cauchy-Schwarz, we have that

where the final inequality follows by the definition of $\upsilon$. So I) implies II) and by symmetry we can show II) implies I). We also have that

so I) implies III) and by symmetry we can show II) implies IV).

It is enough to show that III) implies I). Suppose $\hat{m}_{i}>\hat{m}_{j}$. Then

so it must be the case that $m_{i}-m_{j}>0$ given that $|m_{i}-m_{j}|>\Delta\sin\theta$. ∎

We now show that we can combine these projected center estimates to approximately recover the two-dimensional centers by solving a linear system. The specification of this algorithm, which we call PreConsolidate, is given in Algorithm 2.

so it remains to bound $\sigma_{\min}(A)$. Without loss of generality we may assume $\omega_{1}=(1,0)$ and $\omega_{2}=(x,\sqrt{1-x^{2}})$ for $x\triangleq 1-\upsilon^{2}/2$, in which case $\sigma_{\min}(A)=\upsilon\sqrt{1-\upsilon^{2}/4}$, and the claim follows. ∎

Finally, we show how to boost the success probability via the following naive clustering-based algorithm Select (Algorithm 3), whose guarantees we establish below.

We can now give the full specification of our algorithm LearnAiryDisks (see Algorithm 4).

Let $\rho$ be a $\Delta$-separated superposition of $k$ Airy disks. For any $\epsilon_{1},\epsilon_{2},\delta>0$, let

Without loss of generality suppose $\epsilon_{1}<3\Delta/8$. Then the output $(\lambda^{*}_{1},\lambda^{*}_{2},\bm{\mu}^{*}_{1},\bm{\mu}^{*}_{2})$ of LearnAiryDisks, given $\epsilon_{1},\epsilon_{2},\delta$ and access to an $\eta$-approximate, $O(\log(1/\delta))$-query OTF oracle $\mathcal{O}$ for $\rho$, satisfies

for some permutation $\tau$ with probability at least $1-\delta$. Furthermore, the runtime of LearnAiryDisks is dominated by the time it takes to invoke the OTF oracle $O(\log(1/\delta))$ times.

Suppose we are given a valid $\eta$-approximate OTF oracle $\mathcal{O}$. By taking $\theta=\frac{\pi}{3k^{2}(k-1)}$ and invoking Lemmas 4.1 and 4.5, we ensure that a single run of PreConsolidate in an iteration of the loop in Step 6 of LearnAiryDisks will yield, with probability at least $1-\frac{1}{3k}$, an $(\epsilon^{\prime}_{1},\epsilon^{\prime}_{2})$-accurate estimate, where

In this case we say that such an iteration of the loop in LearnAiryDisks “succeeds.” Note that if we take

then we can ensure that $\epsilon^{\prime}_{1}=\epsilon_{1}/3$ and $\epsilon^{\prime}_{2}=\epsilon_{2}$. The bound in (46) then follows from the elementary inequality $\sin\theta\geq\theta/2$ for $0\leq\theta\leq 1$, together with our choice of $\theta=\frac{\pi}{3k^{2}(k-1)}$.

Each iteration of the loop in Step 6 of LearnAiryDisks individually succeeds with probability at least $1-\frac{1}{3k}$. So by a Chernoff bound, by taking $T=\Omega(\log(1/\delta))$, we conclude that with probability at least $1-\delta$, at least $1-\frac{1}{2k}$ fraction of these iterations will succeed. So of the $k\cdot T$ elements in $\mathcal{S}$, at most $T/2$ correspond to failed iterations.

Observe that $G$ is $k$-partite because every vertex in $V$ is $3\epsilon^{\prime}_{1}$-close to some center of $\rho$, but two vertices which are $3\epsilon^{\prime}_{1}$-close to $\bm{\mu}_{i}$ and $\bm{\mu}_{j}$ respectively for $i\neq j$ must be distance at least $\Delta-6\epsilon^{\prime}_{1}>2\epsilon^{\prime}_{1}$ apart. We conclude that with high probability, Select will output $3\epsilon^{\prime}_{1}=\epsilon_{1}$-accurate estimates for the centers of $\rho$.

Finally, note that in each iteration of the main loop of LearnAiryDisks, $\mathcal{O}$ is invoked exactly six times. Furthermore, other than these invocations of $\mathcal{O}$, the remaining steps of LearnAiryDisks all require constant time. So the runtime of LearnAiryDisks is indeed dominated by the $O(\log(1/\delta))$ calls to $\mathcal{O}$. ∎

3 Learning Airy Disks Above the Diffraction Limit

In this subsection we present the proof of Theorem 4.2. Recall that we are assuming that $\sigma=1/\pi$ and $\Delta>\overline{\gamma}$, where $\overline{\gamma}$ is defined in (21). Let $c\triangleq\frac{1}{2}(\Delta+\overline{\gamma})$ and define $R\triangleq\frac{\overline{\gamma}}{2c}$ and $r\triangleq 1/2-R$.

We will use the following Algorithm 5 that we call TensorResolve. While this is only a slight modification of the tensor decomposition algorithm of [HK15] for high-dimensional superresolution, our analysis is novel and obtains sharper results in low dimensions by using certain extremal functions [Gon18, HV+96, CCLM17] arising in the study of de Branges spaces (see Theorem 4.4.

and note that it admits a low-rank decomposition as

The following is a consequence of the stability of Jennrich’s algorithm.

The setting of parameters in [HK15] is slightly different from ours, so we provide a self-contained proof of Lemma 4.7 in Appendix D.

We will also need the following basic lemma about the stability of solving for $\hat{\lambda}$ in Step 10 in TensorResolve.

By triangle inequality and definition of $\hat{\lambda}$, $\lVert\hat{V}(\hat{\lambda}-\lambda)\rVert_{2}\leq 2\epsilon\lVert\lambda\rVert_{2}+2\epsilon^{\prime}$, so $\lVert\hat{\lambda}-\lambda\rVert_{2}\leq\frac{2\epsilon\lVert\lambda\rVert_{2}+2\epsilon^{\prime}}{\sigma_{\min}(\hat{V})}$. The lemma follows because $\sigma_{\min}(V^{\prime})\geq\sigma_{\min}(V)-\epsilon$. ∎

It remains to show the following condition number bound.

For any $\delta>0$, if $m=\Theta\left(\frac{k^{2}\log(/\delta)}{(\Delta-\overline{\gamma})\wedge 1}\right)$, then $\kappa(V)\leq O\left(k\vee\frac{k}{\sqrt{\Delta-\overline{\gamma}}}\right)$ and $\sigma_{\min}(V)\geq\Omega\left(k^{2}\log(/\delta)\right)$ with probability at least $1-\delta$.

$\kappa(V)\leq\sqrt{2k}\cdot\kappa(V^{*})$.

So by matrix Hoeffding applied to the random variables $m\cdot{V^{*}}^{\dagger}_{1}V^{*}_{1},.\ldots,m\cdot{V^{*}}^{\dagger}_{m}V^{*}_{m}$, each of which is upper bounded in spectral norm by $m\cdot k$ based on (55), we conclude that

Lemma 2 below allows us to bound the quadratic form given by the expectation term evaluated at any $\lambda$. Taking $t=O(\sqrt{\log k/\delta})$ and $m=\Theta\left(\frac{k^{2}\log(k/\delta)}{(\Delta-\overline{\gamma})\wedge 1}\right)$ in (56) and applying Lemma 2, we conclude that with probability at least $1-\delta$,

from which it follows that with this probability, $\kappa(V^{*})\leq O\left(\frac{k}{(\Delta-\overline{\gamma})\wedge 1}\right)^{1/2}$, from which the lemma follows by Lemma 4.10. ∎

It remains to show Lemma 2 below, the key technical ingredient of this section. We will require the following special case of a result of [Gon18], which essentially follows from results of [CCLM17, HV+96]. This can be thought of as the high-dimensional generalization of the well-known Beurling-Selberg minorant (see, e.g., [Vaa85] for a discussion of the one-dimensional case).

The upper bound follows by (55). We now show the lower bound. By Theorem 4.4 applied to $d=2$, for any $\overline{\gamma}/2<r<\frac{j_{1,1}}{\pi}$ there is a function $M$ which minorizes the indicator function of $B^{2}(1)$ and has Fourier transform supported in $B^{2}(r)$. Take $r=\{\Delta R\wedge\frac{\overline{\gamma}/2+j_{1,1}/\pi}{2}\}$ which satisfies $\overline{\gamma}/2<r<\frac{j_{1,1}}{\pi}$. This implies that the function $M^{\prime}(\omega)\triangleq\frac{1}{\pi R^{2}}\cdot\frac{1}{R}\cdot M(\omega/R)$ minorizes the density $\psi(\omega)$, has Fourier transform supported in $B^{2}(r)\subseteq B^{2}(\Delta)$, and satisfies

where in the last step we used that $R<1/2$. We can lower bound (58) by

where the last step follows by (59) and the fact that $\widehat{M^{\prime}}[\mu_{j}-\mu_{j^{\prime}}]=0$ for all $j\neq j^{\prime}$. The lemma follows from noting that $4r-2\overline{\gamma}>\{2\overline{\gamma}(\Delta/c-1)\}\wedge\left\{\frac{2j_{1,1}}{\pi}-\overline{\gamma}\right\}\geq O(\Delta-\overline{\gamma}\wedge 1)$. ∎

Putting everything together, we have the following guarantee:

Let $\rho$ be a $\Delta$-separated superposition of $k$ Airy disks. For any $\epsilon_{1},\epsilon_{2},\delta>0$, let

Without loss of generality suppose $\epsilon_{1}<1/6$. Then the output $(\lambda^{*}_{1},\lambda^{*}_{2},\bm{\mu}^{*}_{1},\bm{\mu}^{*}_{2})$ of TensorResolve, given $\epsilon_{1},\epsilon_{2},\delta$ and access to an $\eta$-approximate, $m$-query OTF oracle $\mathcal{O}$ for $\rho$, satisfies

for some permutation $\tau$ with probability at least $1-\delta$. Furthermore, the runtime of LearnAiryDisks is polynomial in $k$, the number of OTF oracle queries, and the time it takes to make those queries.

and because of the elementary inequality $|e^{-2\pi ix}-1|\geq 2|x|$ for any $|x|\leq 2/3$ and the fact that

To show that the mixing weights are $\epsilon_{2}$-close to the true mixing weights, we can apply Lemma 4.8 to conclude that

so, possibly by modifying $\epsilon_{1}$ to be $\frac{\epsilon_{2}}{k^{2}\log(1/\delta)}$, we get that the estimates $\{\hat{\lambda}_{j}\}$ for the mixing weights are $\epsilon_{2}$-close to the true mixing weights. ∎

Note that we can also amplify the success probability of TensorResolve by running Select from Section 4.2, but we do not belabor this point here.

4 Approximating the Optical Transfer Function

In this section, we show that the following algorithm DFT is a valid implementation of an approximate OTF oracle. We begin by showing that when the samples have granularity $\varsigma=0$, DFT can achieve arbitrarily small error with polynomially many samples.

By a union bound, it suffices to show that for any single $j\in[m]$, $|u_{j}-\widehat{\rho}[\omega_{j}]|\leq\eta$ with probability at least $1-\beta/m$. Note that

where the last step follows by the fact that $\widehat{\rho}$ is real-valued (by circular symmetry of $A$). Furthermore, the summands in $\sum^{N}_{i=1}\cos(2\pi\cdot\langle\omega_{j},\bm{x}_{i}\rangle)$ are $$-valued, so by Chernoff,

from which the lemma follows by our choice of $N$. ∎

We now show that for general granularity $\varsigma>0$, the output of DFT still achieves error $\eta+O(\varsigma)$.

Take any collection of 0-granular samples $\bm{x}^{\prime}_{1},...,\bm{x}^{\prime}_{N}$ for which the averages $u^{\prime}_{1},...,u^{\prime}_{m}$ computed by DFT would be $\eta$-accurate. If DFT were instead passed $\varsigma$-granular samples $\bm{x}_{1},...,\bm{x}_{N}$ for which $\lVert\bm{x}^{\prime}_{i}-\bm{x}_{i}\rVert_{2}\leq\varsigma$ for each $i\in[N]$, then by triangle inequality, the averages $u_{1},...,u_{m}$ computed by DFT with these samples would satisfy $|u_{j}-u^{\prime}_{j}|\leq\eta+O(\varsigma\cdot\lVert\omega_{j}\rVert_{2})$ for each $j\in[m]$, as claimed. ∎

Finally, with Lemma 4.6 and Lemma 4.13, we can complete the proof of Theorem 4.1.

By Lemma 4.6, it suffices to produce an $\eta$-approximate, $m$-query OTF oracle for $\eta$ defined in (46) and $m=O(\log 1/\delta)$. By Corollary 4.2, this can be done using

samples of granularity $\eta/2$ with probability at least $1-\delta$. Theorem 4.1 then follows by a union bound over the failure probabilities of LearnAiryDisks and DFT, and replacing $2\delta$ with $\delta$ and absorbing constant factors. Finally, note that the dependence on $\mathcal{R}$ follows by the discussion at the end of Section 4.1. ∎

By Lemma 4.12, it suffices to produce an $\eta$-approximate, $m$-query OTF oracle for $\eta$ defined in (62) and $m=\Theta\left(\frac{k^{2}}{(\Delta-\overline{\gamma})\wedge 1}\right)$. By Corollary 4.2, this can be done with probability $9/10$ using

samples of granularity $\eta/2$ with probability at least $1-\delta$. Theorem 4.1 then follows by a union bound over the failure probabilities of TensorResolveCorrect and DFT. As in the proof of Theorem 4.1, the dependence on $\mathcal{R}$ follows by the discussion at the end of Section 4.1. ∎

Information Theoretic Lower Bound

for some $0<\sigma<1$. Note that under this setting of parameters, the Abbe limit corresponds to separation $\pi\sigma$.

Let $\underline{\gamma}\triangleq\sqrt{4/3}$. There exists a choice of $\{\bm{\mu}_{i}\}$, $\{\bm{\mu}^{\prime}_{i}\}$, $\{\lambda_{i}\}$, $\{\lambda^{\prime}_{i}\}$ such that the minimum separation among centers of $\rho$ and among centers of $\rho^{\prime}$ is $\Delta=(1-\epsilon)\underline{\gamma}\pi\sigma$, and $d_{\text{TV}}(\rho,\rho^{\prime})\leq\exp(-\Omega(\epsilon\sqrt{k}))$.

We first bound $\lVert\rho-\rho^{\prime}\rVert_{L^{2}}$. By Plancherel’s,

where the inequality follows by the elementary bound

Recall now the construction in Lemma 7 (see the beginning of Section 2). As the entries of the vector $u$ constructed in Lemma 7 satisfy $\mathop{\text{sgn}}(u_{j_{1},j_{2}})=(-1)^{j_{1},j_{2}}$, let $I_{0}$ (resp. $I_{1}$) denote the elements $\bm{j}=(j_{1},j_{2})\in\mathcal{J}\times\mathcal{J}$ for which $j_{1}+j_{2}$ is even (resp. odd), and for every $\bm{j}\in I_{0}$ (resp. $\bm{j}\in I_{1}$), define $\lambda_{\bm{j}}$ and $\bm{\mu}_{\bm{j}}$ (resp. $\lambda^{\prime}_{\bm{j}}$ and $\bm{\mu}^{\prime}_{\bm{j}}$) by $u_{\bm{j}}$ and $\left(\frac{j_{1}}{m},\frac{\sqrt{3}j_{2}}{m}\right)$. This construction is illustrated for $k=25$ in Figure 5. By design, $\{\lambda_{\bm{j}}\}_{\bm{j}\in I_{0}}$ and $\{\lambda^{\prime}_{\bm{j}}\}_{\bm{j}\in I_{1}}$ consist solely of nonnegative scalars and respectively sum to 1, so $\rho,\rho^{\prime}$ are valid superpositions of Airy disks. Furthermore, by design,

so we may now bound (71) to get the desired $L_{2}$ bound

We are now ready to show that $d_{\text{TV}}(\rho,\rho^{\prime})$ is small. The following is a generic $L^{1}$ bound for functions whose univariate restrictions have bounded $L^{2}$ mass, whose derivatives inside some region $\Omega$ are bounded, and which decay sufficiently quickly outside of $\Omega$.

For all $y\in[-T,T]$, $\max_{x\in[-T,T]}|f^{\prime}(x,y)|\leq C$,

$f(-T,y)\leq\delta$ for all $y\in[-T,T]$.

By the triangle inequality and condition 3, it is enough to verify that

Note that for a fixed $y\in[-T,T]$, we have by the fundamental theorem of calculus and conditions 2 and 3 that for any $x\in[-T,T]$,

Define $g(y)\triangleq\int^{T}_{-T}f(t,y)^{2}dt$ and note that $\int^{T}_{-T}g(y)dy\leq\lVert f\rVert^{2}_{L^{2}}$. Then

where the penultimate inequality follows from the measure-theoretic generalization of Jensen’s inequality. ∎

As we will see below, by triangle inequality it will suffice to verify certain properties of $A_{\bm{a},y}$.

For $T>\Delta(\sqrt{k}-1)/4$, we have that $f(-T,y)\leq\Omega\left((T/\sigma)^{-8/3}\right)$ for all $y\in[-T,T]$.

By linearity, it suffices to show that for any $y\in[-T,T]$ and any $j_{1},j_{2}\in\mathcal{J}$, the claimed bound holds for $A_{\nu_{j_{1},j_{2}},y}(-T)$. By Theorem 3.1, we know that

where in the last step we used the fact that $j_{1}/m\leq\Delta(\sqrt{k}-1)/4<T$. ∎

For $T>\Delta(\sqrt{k}-1)/2$, we have that $\int_{\Omega^{c}}|f|\leq O(T^{-2/3}\sigma^{8/3})$, where $\Omega=[-T,T]^{2}$.

By linearity and the fact that $\lVert\bm{x}-(j_{1}/m,\sqrt{3}j_{2}/m)\rVert_{2}\geq T-j_{1}/m\geq T/2$ for every $\bm{x}\not\in\Omega$, it suffices to show that for any $j_{1},j_{2}\in\mathcal{J}$, the claimed bound holds for $\int_{B_{0}(T/2)^{c}}|A(x,y)|dx\,dy$. Expressing this as a polar integral, we have

Take $T=\Theta(\lVert f\rVert^{-1/5}_{L^{2}})\geq\exp(-\Omega(\epsilon\sqrt{k}))$ which we may assume without loss of generality, by scaling $\sigma,\Delta$ appropriately, is greater than $\Delta\sqrt{k}$. By Lemma 5.1 and Lemmas 5.2, 5.3, and 5.4, we have that for $f$ defined by (74),