Efficient Computation of the Zeros of the Bargmann Transform Under Additive White Noise

We study the computation of the zero set of the Bargmann transform of a signal contaminated with complex white noise, or, equivalently, the computation of the zeros of its short-time Fourier transform with Gaussian window. We introduce the adaptive minimal grid neighbors algorithm (AMN), a variant of a method that has recently appeared in the signal processing literature, and prove that with high probability it computes the desired zero set. More precisely, given samples of the Bargmann transform of a signal on a finite grid with spacing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document}δ, AMN is shown to compute the desired zero set up to a factor of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document}δ in the Wasserstein error metric, with failure probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(\delta ^4 \log ^2(1/\delta ))$$\end{document}O(δ4log2(1/δ)). We also provide numerical tests and comparison with other algorithms.

Originally introduced as a link between configuration and phase space in quantum mechanics [7], the Bargmann transform was later recognized as a powerful tool in signal analysis [11] because it encodes the correlations between the signal f and the time-frequency shifts of the Gaussian function g(t) = 2 In the jargon of time-frequency analysis, the right-hand side of (1.2) is known as the short-time Fourier transform of f with Gaussian window, and measures the contribution to f (t) of the frequency y near t = x.
In practice, the values of the short-time Fourier transform of a signal f are only available on (a finite subset of) a grid {(δk, δj) : k, j ∈ Z}, δ > 0, (1.3) and possibly only approximately so due to numerical errors.The goal of Gabor analysis is to extract useful information about f from such limited measurements.Equivalently, by (1.2), the task is to capture the analytic function F given a limited number of its samples on a grid.This second point of view led to the most conclusive results in Gabor theory, such as the complete description of all grids (1.3) for which the Gabor transform fully retains the original analog signal f [11,21,23,24].
While Gabor signal analysis has traditionally focused on large values of the short-time Fourier transform (1.2), recent work has brought to the foreground the rich information stored in its zeros, especially when the signal is contaminated with noise.Heuristically, the zeros of the Bargmann transform of noise exhibit a rather rigid random pattern with predictable statistics, from which the presence of a deterministic signal can be recognized as a salient local perturbation [16,13,14].Remarkably, the Bargmann transform of white noise has been identified as a certain Gaussian analytic random function [5,6], and consequently the wellresearched statistics of their zero sets [19,22] can be leveraged in practice [15,Chapters 13 and 15], [5].The particular structure observed in the zeros of the Bargmann transform under even a moderate amount of white noise has also been invoked as explanation for the sparsity resulting from certain non-linear procedures to sharpen spectrograms [16], as the zeros of the Bargmann transform are repellers of the reassignment vector field [15,Chapter 12].The practical exploitation of such insights requires an effective computational link between finitely given data on the one hand and zeros of Bargmann transforms of analog signals on the other.1.2.Computation of zero sets.Suppose that the values of the Bargmann transform F of a signal f are given on a grid Λ = {δk + iδj : k, j ∈ Z}, δ > 0, (1.4) and we wish to compute an approximation of {F = 0}, the zero set of F , within the square Ω L = {x + iy : |x|, |y| ≤ L}. (1.5)More realistically, we only have access to samples of F on those grid points near the computation domain, e.g., on Λ L = {δk + iδj : k, j ∈ Z, |δk|, |δj| ≤ L}. (1.6)The inverse of the spacing of the grid, 1/δ, will be called the resolution of the data.Thresholding.The most naive approach to compute {F = 0} is thresholding: one selects all grid points λ such that |F (λ)| is below a certain threshold ε > 0: The normalizing weight e − 1 2 |λ| 2 is motivated by (1.2), as the short-time Fourier transform of a typical signal can be expected to be bounded.One disadvantage of this approach is that it requires an educated choice for the threshold ε.Moreover, computations with various reasonable choices of thresholds, such as quantiles of e − 1 2 |λ| 2 |F (λ)| calculated over all grid points λ, either fail to compute many of the zeros or capture too many points (see Figure 1).
Extrapolation.One may consider using the samples of F on the finite grid (1.6) to reconstruct the signal f , resample F at arbitrarily high density, and thus calculate more easily the zero set {F = 0}.However, computation of zeros through extrapolation may be inaccurate: while the samples of F on the infinite grid (1.4) determine F as soon as √ π • δ < 1 [11,21,24], the truncation errors involved in the approximation of F near Ω L from finite data (1.6) can only be neglected at very high resolution 1/δ.Even if the values of F are successfully extrapolated to a higher resolution grid, the remaining computation is still not trivial, as, for example, simple thresholding may fail even at high resolution (see Figure 4 and Section 5).
Minimal Grid Neighbors.A greatly effective numerical recipe for the computation of zeros of the Bargmann transform can be found in the code accompanying [13] -although not explicitly described in the text.A grid point λ is selected as a numerical approximation of a zero if e − 1 2 |λ| 2 |F (λ)| is minimal among grid neighbors, i.e., where |z| ∞ = max{|x|, |y|}.The subset of points that pass the test furnish the computation of {F = 0}.This method, which we call minimal grid neighbors (MGN), performs impressively as long as the grid resolution is moderately high.Indeed, we understand that the method is behind the simulations in [15,Chapter 15] which quite faithfully reproduce the statistics of the zeros of the Bargmann transform of complex white noise (that are known analytically [5,19]).The MGN algorithm was also used to produce the plots in [5], as pointed out in [5, Section 5.1.1.]; see also [6,Section 5], [20, Section IV], and [1].Heuristically, the test (1.8) succeeds in identifying zeros due to the analyticity of F , which implies that |F (z)|e − The MGN algorithm performs equally well when calculating the zeros of the Bargmann transform of a signal composed of a deterministic real-variable function f 1 plus complex white noise with variance σ 2 > 0. The presence of a certain amount of randomness must be behind the success of the algorithm, as, for σ = 0, the method cannot be expected to succeed.Indeed, one can engineer a deterministic signal f where the detection of zeros fails, as the value of its Bargmann transform F can be freely prescribed on any given finite subset of the computation domain [24].We are unaware of performance guarantees for MGN.1.3.Contribution.In this article we introduce a variant of MGN, called Adaptive Minimal Grid Neighbors (AMN).The algorithm is based on a comparison similar to (1.8) but incorporates an adaptive decision margin, that depends on the particular realization of F .While AMN has the same mild computational complexity and similar practical effectiveness as MGN, we are able to estimate the accuracy and confidence of the computation with AMN in terms of the grid resolution.In this way, we show that AMN is probably approximately correct for the signal model (1.9), in the sense that it computes the zero set with high probability up to the resolution of the data.
On the one hand, we present what to the best of our knowledge are the first formal guarantees for the approximate computation of zero sets of analytic functions from grid values.In fact, besides its main purpose of computation with specific data, the AMN algorithm offers a computationally attractive and provably correct method to simulate zero sets of the Gaussian entire function (2.10), by running the procedure with simulated inputs.On the other hand, our analysis is a first step towards understanding the performance of MGN.

Main Result
2.1.The adaptive minimal neighbors algorithm.We now introduce a new algorithm to compute zero sets of Bargmann transforms.Suppose again that samples of an analytic function F : C → C are given on those points of the grid (1.4) that are near the computation domain (1.5), say, on For each grid point λ strictly inside the computation domain, λ ∈ Λ L = Λ∩Ω L , we use the neighboring sample at λ + δ to compute the following comparison margin: The margin η λ therefore depends on the particular realization of F .To motivate the definition, note that, when λ = 0, the maximum in (2.2) is taken over |F (0)| and the absolute value of an incremental approximation of ∂F (0), where ∂F = ).The comparison margin thus incorporates the size and oscillation of F at z = 0.In general, η λ has a similar interpretation with respect to the covariant derivative and, indeed, η λ is defined so that The differential operator ∂ * F plays a distinguished role in the analysis of vanishing orders of Bargmann transforms [10,12] because it commutes with the translational symmetries of the space that they generate (the Bargmann-Fock shifts defined in Section 3.2).
The first step of the algorithm selects all grid points λ ∈ Λ L that pass the following comparison test: In contrast to (1.8), the comparison in (2.5) does not involve the immediate grid neighbors of λ but rather those points lying on the square centered at λ with half-side-length 2δ; see Figure 2. (In particular, the test only involves grid points µ ∈ Λ L+2δ .)Intuitively, the larger distance between λ and µ permits neglecting the error in the differential approximation (2.4).Selection step: For each grid point λ ∈ Λ L = Λ ∩ Ω L inside the target domain Ω L , we define the following comparison margin: The grid point λ is then selected if the following test is satisfied: (Note that the test (2.7)only involves grid points µ ∈ Λ L+2δ .)Let Z 1 be the set of all selected grid points.
Sieving step: Use an off-the-shelf clustering algorithm to select a subset Z ⊂ Z 1 that is 5δ separated: and maximal with respect to this property, i.e., no proper superset satisfies (2.8).
(One concrete implementation of the sieving step is described in Section 5.2.1.) Output: The set Z.
Remark 2.1.The constants 3/4 in (2.6), 2 in (2.7), and 5 in (2.8) are to some extent arbitrary, and other choices lead to similar results.These particular values are chosen to aid the exposition rather than to optimize practical performance.In fact, these choices are suboptimal at low resolutions (see Section 5).

Performance guarantees for AMN.
To study the performance of the AMN algorithm we introduce the following input model, which, as we will explain, corresponds to the Bargmann transform of an arbitrary signal contaminated with complex white noise of an arbitrary intensity.
Input model.We consider a random entire function on the complex plane where F 1 : C → C is a deterministic entire function, F 0 is a (zero mean) Gaussian analytic function with correlation kernel: z, w ∈ C, (2.10) and σ > 0 is the noise level.We assume that the deterministic function F 1 satisfies the quadratic exponential growth estimate for some constant A ≥ 0.
As for F 0 , the assumption means that for each z 1 , . . ., z n ∈ C, (F 0 (z 1 ), . . ., F 0 (z n )) is a normally distributed (circularly symmetric) complex random vector, with mean zero and covariance matrix e z k z k, .Alternatively, F 0 can be described as where (ξ n ) n≥0 are independent standard complex random variables [19].
Discussion of the model.The random function F 0 is the Bargmann transform of standard complex white noise N .As each realization of complex white noise is a tempered distribution, the computation of its Bargmann transform (1.1) is also to be understood in the sense of distributions, as in [8].(See [6] and [17, Section 5.1] for a detailed discussion on this, and alternative approaches.) We can similarly interpret F 1 as the Bargmann transform of a distribution f 1 on the real line [8].The assumption (2.11) means precisely that f 1 belongs to the modulation space M ∞ (R) consisting of distributions with Bargmann transforms bounded with respect to the standard Gaussian weight -or, equivalently, with bounded short-time Fourier transforms [9].The modulation space M ∞ (R) includes all square-integrable functions f 1 ∈ L 2 (R) and also many of the standard distributions used in signal processing.
In summary, the input model (2.9) corresponds exactly to the Bargmann transform of a random signal where f 1 ∈ M ∞ (R) and σ • N is complex white noise with standard deviation σ.
Performance analysis.We now present the following performance guarantees, pertaining to the computation domain (1.5) and the acquisition grid (2.1).To avoid immaterial technicalities, we assume that the corners of the computation domain lie on the acquisition grid.
Theorem 2.2.Fix a domain width L ≥ 1, a noise level σ > 0, and a grid spacing δ > 0 such that L/δ ∈ N. Let a realization of a random function F as in (2.9) with (2.10) and (2.11) be observed on Λ L+2δ , and let Z be the output of the AMN algorithm.
There exists an absolute constant C such that, with probability at least there is an injective map Φ : {F = 0} ∩ Ω L → Z with the following properties: • (Each zero is mapped into a near-by numerical zero) • (Each numerical zero that is away from the boundary arises in this form) A proof of Theorem 2.2 is presented in Section 4. We remark some aspects of the result.
• The AMN algorithm does not require knowledge of the noise level σ and is homogeneous in the sense that F and cF , with c ∈ C \ {0} produce the same output.
• Within the estimated success probability, the computation is accurate up to a factor of the grid spacing.
• The analysis concerns an arbitrary deterministic signal impacted by noise and is uniform over the class (2.11).As usual in such smoothed analysis, the success probability grows as the signal to noise ratio A/σ decreases, because randomness helps preclude the very untypical features that could cause the algorithm to fail [25].In fact, in the noiseless limit σ = 0, the algorithm could completely fail, since F 1 can be freely prescribed on any finite subset of the plane [24].For example, irrespectively of its values on the acquisition grid, the deterministic function F 1 could have a cluster of zeros of small diameter that would trigger a single positive minimality test.The proof of Theorem 2.2 shows that such examples are fragile, as the addition of even a moderate amount of noise regularizes the geometry of the zero set.
• Up to a small boundary effect, the guarantees in Theorem 2.2 comprise an estimate on the Wasserstein distance between the atomic measures supported on {F = 0} ∩ Ω L and on the computed set Z.More precisely, for a tolerance θ > 0 let us define the boundary-corrected Wasserstein pseudo-distance between two sets U, V ⊆ C as where the infimum is taken over all injective maps Φ : . (The definition is not symmetric in U and V , but this is not important for our purpose.)Then Theorem 2.2 reads • The presented analysis concerns a signal contaminated with complex white noise.This is a mathematical simplification; we believe that with more technical arguments a similar result can be derived for real white noise.The case of colored noise seems more challenging and will be the object of future work.

Numerical experiments.
In Section 5, we report on numerical experiments that compare the AMN and MGN algorithms.We also include a modified version of thresholding (ST), that uses a thereshold proportional to the grid spacing and incorporates a sieving step as in AMN (while standard thresholding without sieving performs extremely poorly, as seen in Figure 1).
The performance of AMN, MGN, and ST is first tested indirectly, by using these algorithms to simulate the zero sets of the random functions in the input model (2.9).We then compare theoretically derived statistics of the zeros of (2.9) to empirical statistics obtained from the output of AMN and MGN under various simulated realizations of (2.9).
Second, we perform a consistency experiment that aims at estimating the probability of computing a low-distortion parametrization of the zero set of F , as in Theorem 2.2.Specifically, we simulate a realization of the random input F sampled at high-resolution and use the output of AMN or MGN as a proxy for the ground truth {F = 0}.We then test the extent to which this set is captured by the output of AMN, MGN, or ST from lower resolution subsets of the same simulated data.
The performance of AMN and MGN is almost identical, although the minimal resolution at which MGN starts to perform well is slightly lower than that for AMN.(This is to be expected, as the constants 2 and 5 used in (2.7) and (2.8) are not adequate for low resolutions, cf.Remark 2.1.)Both AMN and MGN significantly outperform ST.See also Figure 4 for an illustration.
The favorable performance of AMN is interesting also when the input is just noise, as it gives a fast and provably accurate method to simulate the zeros of the Gaussian entire function (2.10).(The simulations that we present in Section 5 use certain heuristic shortcuts to accelerate the simulation of the input (2.10) -see Section 5.1; although we do not formally analyze these, they are implicitly validated, as the simulated point process reproduces the expected theoretical statistics.)All numerical experiments can be reproduced with openly-accessible software and our code is available at https://github.com/laescudero/discretezeros.Our implementation of the Bargmann transform uses [4].

AMN and MGN ST
2.4.Organization.Section 3 introduces the notation and basic technical tools about analytic functions, Bargmann-Fock shifts, and their applications to random functions and their zeros.Theorem 2.2 is proved in Section 4, while numerical experiments are presented in detail in Section 5. Conclusions and outlook on future directions are discussed in Section 6.

Preliminaries
3.1.Notation.For a complex number z = x + iy, we use the notation |z| ∞ = max{|x|, |y|}, while |z| denotes the usual absolute value.The zero set of F is denoted by {F = 0}.The differential of the (Lebesgue) area measure on the plane will be denoted for short dm, while the measure of a set E is |E|.With a slight abuse of notation, we also denote the cardinality of a finite set Z by |Z|.Squares on the complex plane are denoted by Q r (z) = {w ∈ C : |z − w| ∞ ≤ r}.For two non-negative functions f, g, we write f g if there exists an absolute constant C such that f (x) ≤ Cg(x), for all x.We write f g if f g and g f .

The Wirtinger derivative of a function
).When we need to stress on which variable the derivative is taken we write subindices, e.g., ∂ w F (z, w).
A Gaussian entire function (see [19,Ch. 2] and [22]) is a random function F : C → C that is almost surely entire, and such that for every z 1 , . . ., z n ∈ C, F (z 1 ), . . ., F (z n ) is a circularly symmetric complex normal vector.We will be only concerned with the random function F given in (2.9).We also use the notation (1.4), (1.5), (1.6), possibly for distinct values of L.
3.2.Bargmann-Fock shifts and stationarity of amplitudes.The analysis of the AMN algorithm is more transparent when formulated in terms of the Bargmann-Fock shifts.For a function F : C → C we let The amplitude of an entire function F is defined as the weighted magnitude and satisfies The comparison margin of the AMN algorithm (2.6) can be expressed in terms of Bargmann-Fock shifts as and leads to the approximation (2.4) because (Here and throughout we write F λ (0) for ∂[F λ ](0).)Similarly, in terms of amplitudes, the test (2.7) reads With respect to the input model (2.9) we note that, if F 0 is the (zero mean) Gaussian entire function with correlation kernel (2.10), then the Bargmann-Fock shifts F 0 → F 0 w preserve the stochastics of F 0 , as they leave its covariance kernel invariant.As a consequence, for any w ∈ C, F 0 w (0), F 0 w (0) are independent standard complex normal random variables (with zero mean and variance 1).Indeed, by the mentioned invariance it suffices to consider w = 0, and, in this case, F 0 w (0), F 0 w (0) are the coefficients ξ 0 and ξ 1 in (2.12).
Hence, the maximum principle for subharmonic functions together with (3.7) implies that log H(z) and therefore |F (z As F is non-vanishing on D, it follows that ∂ z F − zF = 0 on D, and therefore on D. This contradiction shows that F must vanish on D.

Linearization.
In what follows, we derive basic facts about the input model (2.9), and always assume that (2.10) and (2.11) hold.
The following is a strengthened version of [19,Lemma 2.4.4].Lemma 3.2.Let F be as in (2.9).Then there exists an absolute constant C > 0 such that for all L ≥ 1 and t ≥ A, Proof.We consider the Taylor expansion of F : where E 2 can be bounded in terms of the amplitude (3.2) as We also note that for |z| ≤ 10, Hence, for |z| ≤ 10, We apply the previous bounds to F w , note that, by (3.3), G w (z) = G(z + w), and obtain that for |z| ≤ 10, Let G 0 and G 1 be the amplitudes corresponding to F 0 and F 1 , respectively.Then by (2.11), Hence, To conclude, we claim that the following excursion bound holds: where C > 0 is an absolute constant.For L ≤ 1/4 this follows for example from [19,Lemma 2.4.4].In general, we cover the domain with L 2 squares of the form w + [−1/4, 1/4] 2 , apply the previously mentioned bound to G 0 w (z) = G 0 (z + w), and use a union bound.This completes the proof.

Almost multiple zeros.
It is easy to see that, almost surely, the random function (2.9) has no multiple zeros.In the analysis of the AMN algorithm, we will also need to control the occurrence of zeros that are multiple up to a certain numerical precision, in the sense that F and its derivative are simultaneously small.The following lemma is a first step in that direction, as it controls the probability of finding a grid point that is an almost multiple zero.Lemma 3.3.Let F be as in (2.9) and α, β > 0. Then the probability that for some grid point λ ∈ Λ L the following occurs: , where C is an absolute constant.
3.6.First intensity of zeros.The following proposition is not used in the proof of Theorem 2.2, but rather as a benchmark in the numerical experiments (Section 5).Proposition 3.4.Let F be as in (2.9).Then for every Borel set where Proof.The set of zeros and thus ρ 1 (ζ) does not change if we scale F by a fixed constant.Hence, by considering the function 1 σ F in place of F , we can assume that σ = 1.The expected number of points {z ∈ B : F (z) = 0} of a Gaussian random field F is given by Kac-Rice's formula: where p F (ζ) (0) is the probability density of F (ζ) at 0; see, e.g., [3,Th. 6.2].

Proof of Theorem 2.2
We present the proof of Theorem 2.2 in several steps.The strategy is two-fold: (i) to show that computed zeros are close to true ones, we relate the comparison test (2.7) to a similar property involving non-grid points and apply the minimum principle from Lemma 3.1; (ii) to show that true zeros do trigger a detection, we show that the test (2.7) is satisfied by linearly approximating the input function.The two objectives are in tension: while a large comparison margin η λ would facilitate (i) by absorbing possible oscillations between a grid and a close-by nongrid point, a small margin makes the comparison test easier to satisfy and thus facilitates (ii).The core of the proof consists in showing that the adaptive margin (2.6) strikes the desired balance with high probability.
Initially we bound the Hausdorff distance between the exact and computed zero sets (showing that each of the sets lies in a small neighborhood of the other).We then refine this conclusion to a bound on the Wasserstein distance by analyzing the sieving step.4.1.Preparations.Let L, σ, δ, and F satisfy the assumptions of the theorem, and denote by Z 1 the set produced by the AMN algorithm after the selection step.Recall that L ≥ 1.By choosing a sufficiently large constant in (2.14), we can assume that δ ≤ 1/5; otherwise, the success probability would be trivial.For the same reason, we can assume that 4.2.Excluding bad events.We let γ = 8σ and wish to apply Lemma 3.2 with t = γ log(1/δ).By (4.1), Hence t ≥ A and we can apply Lemma 3.2 to conclude that for all w ∈ Ω L+1 and |z| ≤ 10, except for an event of probability at most 2 − A 2 , we further have Second, we select a large absolute constant κ > 1 to be specified later, and use Lemma 3.3 with both), (4.4) except for an event with probability at most L 2 log 2 (1/δ)δ 4 .
Overall we have excluded events with total probability In what follows, we show that under the complementary events the conclusions of Theorem 2.2 hold.

Numerical Experiments
In this section we perform a series of tests of the AMN algorithm and compare its performance with MGN and thresholding supplemented with a sieving step (ST).5.1.Simulation.We first discuss how to simulate samples from the input model (2.9).To make simulations tractable, we introduce a fast method to draw samples of the Gaussian entire function F 0 given by (2.10) on the finite grid (1.6).The method is based on the relation between the Bargmann transform and the shorttime Fourier transform (1.2) and amounts to discretizing the underlying signal f .We fix L > 0, T > 0, and δ > 0. For convenience, we further let σ = 1 and assume that T δ −1 is an integer.Recall that we also assumed that Lδ −1 is an integer.
To model a discretization of N , we take i.i.d.noise samples in the interval [−T − L, T + L] ⊆ R spaced by a distance δ.More specifically, we consider a random vector w = (w −(T +L)δ −1 , . . ., w (T +L)δ −1 ), where the elements w s ∼ N C (0, δ) are independent, i.e., E[w s w s ] = δ, and E[w s w s ] = 0 for s = s .Here, w s can be interpreted as an integration of N over the interval [δs, δ(s + 1)].
Let f 1 : R → C, ϕ = g| [−T,T ] the restriction of g(t) = ( 2π ) 1 4 e −t 2 to the compact support [−T, T ] and define for k, j ∈ {−Lδ −1 , . . ., Lδ −1 }.The mean of H is given by and approximates the integral x)e −2iyt dt = e −ixy e − 1 2 (x 2 +y 2 ) F 1 (z), with x = δk and y = δj.Furthermore, the covariance of H is For small δ and sufficiently large T , this is an approximation of the integral ∞ −∞ g t − δk g t − δk e −2it(j−j )δ dt = e − u 2 +v 2 +x 2 +y 2 2 e i(uv−xy) e (x−iy)(u+iv) , (5.2) with x = δk, y = δj, u = δk , and v = δj .Therefore, if we take T large enough so that we can ignore the numerical error introduced by the truncation of the normalized Gaussian window g, we obtain in (5.1) a random Gaussian vector whose covariance structure approximates the right-hand side of (5.2) on the grid Λ L , provided that δ is small.
To obtain a vector whose covariance structure approximates (2.10) we proceed as follows.By conjugating z in (5.1) and multiplying by the deterministic factor e −ixy , we obtain an approximate sampling of (2.9) with weight e − 1 2 |z| 2 : (5. 3) e − 1 2 |z| 2 F (z) ≈ e −ixy H(z) for z = δk + iδj.We carry out all computations with the weighted function (5.3), as the unweighted version can lead to floating point arithmetic problems.Note that, for a grid point λ, the comparison margin (2.6) can be expressed in terms of 5.2.Specifications for the experiments.

5.2.1.
Implementation of the sieving step in AMN.In order to fully specify the AMN algorithm we need to fix an implementation of the sieving step, which provides a subset Z ⊆ Z 1 satisfying (2.8), and such that no proper superset Z 1 ⊇ Z Z satisfies (2.8).We choose an implementation that uses knowledge of the input F to decide which points are to be discarded.We assume that Z 1 is non-empty, otherwise Z is trivial.
Algorithm S1: Obtain a maximal subset that is 5δ separated.
Input: Values of a function F on a grid Λ L .A discrete non-empty set Z 1 ⊆ Λ L .
Step 1: Copy the set Z 1 to Z aux 1 .
Step 6: If the set Z aux 1 is not empty, repeat Steps 3-6.If the set Z aux 1 is empty, then the algorithm ends.
Output: The set Z.
The resulting set Z ⊆ Z 1 always satisfies (2.8).Moreover, any superset Z Z included in Z 1 must contain some of the discarded points µ ∈ Z 1 , which by construction satisfy (5.5) for some λ ∈ Z, and therefore Z is not 5δ separated.Thus, Z is indeed maximal with respect to (2.8).
The choice of λ ∈ Z aux 1 in Step 3 of S1 is not essential.Our particular choice is motivated by finding the zeros of F ; however, we did not observe any significant performance difference when using other algorithms than S1 as the sieving step of AMN.5.2.2.Specification of the compared algorithms.Given the values of a function F : C → C on the grid Λ L , we consider the following three algorithms to compute an approximation of {F = 0} ∩ Λ L−1 .
• AMN: the AMN algorithm run with domain length L − 1 and with sieving step S1 implemented as described in Section 5.2.1, • MGN: outputs the set of all grid points λ ∈ Λ L−1 such that • ST: outputs the set of grid points λ ∈ Λ L−1 obtained as the result of applying the sieving algorithm S1 to Note that each of the algorithms relies only on the samples of F on Λ L+2δ−1 .The use of a common input grid Λ L simplifies the notation when considering various grid spacing parameters δ.

5.2.3.
Varying the grid resolution.In the numerical experiments, we start with a small minimal spacing value δ = δ Hi , that provides a high resolution approximation in (5.1), and simulate F as in Section 5.1.We then incrementally double δ to produce coarser grid resolutions and subsample F accordingly.More precisely, each element of the grid Λ L can be written as for adequate M , N > 0. If F is given on Λ L , we subsample it by setting S(F )(λ k,l ) := F λ 2k+i2l , (5.8) for values (k, l) such that the indices 2k + i2l are valid.5.3.Faithfulness of simulation of zero sets.As a first test, we simulate random inputs from the model (2.9), as specified in Section 5.1, apply the abovedescribed three different algorithms, and test whether this process faithfully simulates the zero sets of the random function (2.9).To this end, we estimate first or second order statistics on the computed zero sets by averaging over several realizations of (2.9), and compare them to the corresponding expected values concerning the zero sets of (2.9).

5.3.1.
No deterministic signal.We first consider the case F 1 ≡ 0 and σ = 1 in (2.9).Let F δ Hi 1 , . . ., F δ Hi R be R independent realizations of samples of (2.9) on a grid Λ L with resolution δ = δ Hi , simulated as in Section 5.1.These are then subsampled with (5.8) yielding F δ k r = S (k) (F δ Hi r ) and used as input for AMN, MGN, and ST, as specified in Section 5.2.2.The corresponding output sets are denoted Z δ r where we omit the dependence on the method to simplify the notation.These sets should approximately correspond to {F r = 0} ∩ Λ L−1 , for R independent realizations of (2.9).We now put that statement to test.
The expected number of zeros of the random function F on a Borel set Θ ⊆ C is see, e.g., [19,Section 2.4].We define the following empirical estimator for the first intensity ρ 1 = 1/π: If the computed set Z δ r were replaced by {F = 0} in (5.10), the estimator would be unbiased.The mean of the estimation error ρ(Θ, r, δ) − 1/π thus measures the quality of the algorithm used to compute Z δ r , as it should be close to zero when the algorithm is faithful.In Table 1, we present the empirical means and the empirical standard deviations of the estimation error over R = 1000 independent realizations F δ r for L = 7, Θ = Ω L−1 , T = 6, and various grid sizes δ.To derive a benchmark for the empirical standard deviation of ρ(Θ, r, δ) − 1/π, we express the variance of |{F = 0} ∩ Θ|/|Θ| in terms of the second intensity function ρ 2 (ζ, ζ ) of {F = 0} as follows: A formula for ρ 2 (ζ, ζ ) is provided in [18] and numerical integration over Θ = Ω L−1 results in Var[|{F = 0} ∩ Θ|/|Θ|] ≈ 0.01165.We see in Table 1 that the methods AMN and MGN almost perfectly match the expected mean and standard deviation while ST does not.We only test first order statistics of the computed zero sets.The benchmark is provided by Proposition 3.4: the expected number of zeros of F in Θ is (5.12) where ρ 1 is given by (3.11) (with σ = 1).For each of the tested algorithms, we define an estimator for the error resulting from replacing {F = 0} in (5.12) by the computed set Z δ r (for 1 ≤ r ≤ R): As before, we simulate R = 100 realizations of F = F 0 + F 1 on a grid with a certain spacing δ.The empirical average of β(Θ, r, δ) over all realizations is denoted β R (Θ, δ).As ρ 1 is not constant when F 1 = 0, this time we calculate β R (Θ, δ) on Θ = Ω L 1 for several values of L 1 .
The results for δ = 2 −9 are depicted in Figure 5.We see that the performance of AMN and MGN is indistinguishable, while ST may perform poorly even at such high resolution.Lower grid resolutions yield similar results.5.4.Failure probabilities and consistency as resolution decreases.Having tested the statistical properties of the computed zero sets under the input model (2.9) we now look into the accuracy of the computation for an individual realization F .We aim to test the existence of a map as in Theorem 2.2, that assigns true zeros to computed ones with small distortion and almost bijectively.As a proxy for the Table 2. Functions f 1 and their Bargmann transforms F 1 = B(f ).(unavailable) ground truth {F = 0} we will use the output of AMN from data at very high resolution (computations with MGN yield indistinguishable results).We thus conduct a consistency experiment, where the zero set of the same realization of F is computed from samples on grids of different resolution, and the existence of a map as in Theorem 2.2 between both outputs is put to test.
Suppose that samples of a function F are simulated on a high-resolution grid Λ L with spacing δ = δ Hi and restricted to the low-resolution grid Λ L with spacing δ = δ Lo by subsampling.We compute Z δ Hi ⊆ Ω L−1 from the high-resolution data using AMN, and Z δ Lo ⊆ Ω L−1 from the low-resolution data, using one of the algorithms described in Section 5.2.2.
Second we construct a set U ⊆ Z δ Hi and a map φ : U → Z δ Lo r with the following greedy procedure:

Construction of U and φ
Input: Two subsets of Ω L−1 : Z δ Hi and Z δ Lo .
Step 1: Choose a total order on Z δ Hi .Let U and U be empty sets.If Z δ Hi is empty, output U = ∅ and φ = ∅.Otherwise proceed to Step 2.

Output:
The set U and the map φ.
The resulting function φ is injective and satisfies We say that the computation of Z δ   ∩ Ω (L−1)−2δ Lo ⊆ φ(U ).(5.14)In this case, the map φ satisfies properties analogous to the ones in Theorem 2.2.Conceivably, other such maps may exist even if the one constructed in the greedy fashion fails to satisfy (5.14).We define the following computation certificate: M( Z δ Hi , Z δ Lo ) = 0 if (5.14) holds 1 otherwise.
analysis is the assumption that grid samples of the Bargmann transform are exactly given, while, more realistically, acquired data corresponds to averages of the signal values resulting from analog to digital conversion and numerical integration.Second, we considered complex-valued white noise, while in practice noise may also be colored or real-valued.We understand that the techniques used to prove Theorem 2.2 are general enough to allow for a refinement of the result in these directions.Similarly, we expect to be able to adapt our analysis of AMN to other ensembles of analytic functions, which are relevant in connection to other signal transforms.A more challenging open direction is the investigation of rigorous performance guarantees for MGN, which remains the algorithm of choice in practice.

Figure 1 .
Figure 1.Calculation of zero sets by thresholding: Values below the same threshold (marked by circles) fail to detect some zeros and at the same time cannot clearly separate other zeros.
(a) One of the immediate neighbors is used to calculate a comparison margin.(b) The weighted values of F are compared against points in the larger box.

Figure 2 .
Figure 2. The selection step of AMN.
(a) A point passes the selection test.(b) The same zero causes a second detection.

Figure 4 .
Figure 4.A realization of e − 1 2 |z| 2 |(F 0 (z) + F 1 (z))| with F 1 the Bargmann transform of f 1 .The deterministic functions are scaled to obtain the prescribed A. Zeros computed with AMN, MGN, and ST are calculated from grid samples with δ = 2 −9 .Zeros from AMN and MGN coincide (circle), while ST (cross) fails either by detecting false zeros (left) or by not capturing all of them (right).

1 , 2 − 9 ) f 1 1 , 2 − 9 ) f 1 Figure 5 .
Figure 5. Empirical mean of β(Θ, r, δ) for different choices of f 1 and A, increasing domain Θ = Ω L 1 for L 1 < L, and the three methods.Note the different scale in the bottom right plot illustrating a systematic error in the ST method.

r
was certified to be accurate ifZ δ Hi r ⊆ U and Z δ Lo r