Multiscale High-Dimensional Sparse Fourier Algorithms for Noisy Data

TL;DR

This paper introduces a multiscale sparse Fourier algorithm that is both efficient and robust against noise, improving upon previous methods for high-dimensional signals that are sparse or nearly sparse in the Fourier domain.

Contribution

It presents a novel multiscale algorithm capable of handling noisy high-dimensional Fourier data efficiently, extending prior work to practical noisy scenarios.

Findings

Algorithm is robust against noise in high-dimensional Fourier data.

Achieves efficient runtime comparable to prior sparse Fourier methods.

Effectively approximates nearly sparse signals in noisy environments.

Abstract

We develop an efficient and robust high-dimensional sparse Fourier algorithm for noisy samples. Earlier in the paper ``Multi-dimensional sublinear sparse Fourier algorithm" (2016), an efficient sparse Fourier algorithm with $\Theta(ds \log s)$ average-case runtime and $\Theta(ds)$ sampling complexity under certain assumptions was developed for signals that are $s$-sparse and bandlimited in the $d$-dimensional Fourier domain, i.e. there are at most $s$ energetic frequencies and they are in $ \left[-N/2, N/2\right)^d\cap \mathbb{Z}^d$. However, in practice the measurements of signals often contain noise, and in some cases may only be nearly sparse in the sense that they are well approximated by the best $s$ Fourier modes. In this paper, we propose a multiscale sparse Fourier algorithm for noisy samples that proves to be both robust against noise and efficient.

Figures

Equations

$$f(\mathbf{x})=\sum _{j=1}^s{a}_j{e}^{2\pi i{\mathbf{w}}_j\cdot \mathbf{x}}$$

$$h(\mathbf{t}):=f(g(\mathbf{t})).$$

$${\mathbf{h}}_p^{\tilde{k}}[\ell ]:=h\left(\frac{\ell }{p}{\mathbf{e}}_{\tilde{k}}\right)\text{and}{\mathbf{h}}_{p;ϵ}^{\tilde{k};k}[\ell ]:=h\left(\frac{\ell }{p}{\mathbf{e}}_{\tilde{k}}+ϵ{\mathbf{e}}_k\right)$$

$$\mathcal{F}\left({\mathbf{h}}_p^{\tilde{k}}\right)[m]$$

$$\mathcal{F}\left({\mathbf{h}}_{p;ϵ}^{\tilde{k};k}\right)[m]$$

$$w_{\bar{j},k}^{\prime }=\frac{1}{2\pi ϵ}\mathrm{Arg}\left(\frac{\mathcal{F}\left({\mathbf{h}}_{p;ϵ}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{h}}_p^{\tilde{k}}\right)[m]}\right)\text{and}{a}_{\bar{j}}=\frac{\mathcal{F}\left({\mathbf{h}}_p^{\tilde{k}}\right)[m]}{p},$$

$$\left|\frac{\mathcal{F}\left({\mathbf{h}}_{p;ϵ}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{h}}_p^{\tilde{k}}\right)[m]}\right|=1,$$

$${\mathbf{r}}_p^{\tilde{k}}[\ell ]:={\mathbf{h}}_p^{\tilde{k}}[\ell ]+n_{\ell }=h\left(\frac{\ell }{p}{\mathbf{e}}_{\tilde{k}}\right)+n_{\ell }$$

$$\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]=\mathcal{F}\left({\mathbf{h}}_p^{\tilde{k}}\right)[m]+\sum _{\ell =0}^{p-1}{n}_{\ell }{e}^{-2\pi im\frac{\ell }{p}}$$

$$\mathbb{E}\left[\sum _{\ell =0}^{p-1}{n}_{\ell }{e}^{-2\pi im\frac{\ell }{p}}\right]=0$$

$$\text{Var}\left[\sum _{\ell =0}^{p-1}{n}_{\ell }{e}^{-2\pi im\frac{\ell }{p}}\right]=p{\sigma }^2,$$

$$\mathbb{E}\left[\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]\right]=\mathcal{F}\left({\mathbf{h}}_p^{\tilde{k}}\right)[m]$$

$$\text{Var}\left[\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]\right]=p{\sigma }^2.$$

$$\mathbb{E}\left[\mathcal{F}\left({\mathbf{r}}_{p;ϵ}^{\tilde{k};k}\right)[m]\right]=\mathcal{F}\left({\mathbf{h}}_{p;ϵ}^{\tilde{k};k}\right)[m]$$

$$\text{Var}\left[\mathcal{F}\left({\mathbf{r}}_{p;ϵ}^{\tilde{k};k}\right)[m]\right]=p{\sigma }^2.$$

$$\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]=p{a}_{\bar{j}}+\mathcal{O}(\sigma \sqrt{p})$$

$$\mathcal{F}\left({\mathbf{r}}_{p;ϵ}^{\tilde{k};k}\right)[m]=p{a}_{\bar{j}}{e}^{2\pi i{w}_{\bar{j},k}^{\prime }ϵ}+\mathcal{O}(\sigma \sqrt{p})$$

$$\frac{\mathcal{F}\left({\mathbf{r}}_{p;ϵ}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]}=e^{2\pi i{w}_{\bar{j},k}^{\prime }ϵ}+\mathcal{O}\left(\frac{\sigma }{a_{\bar{j}}\sqrt{p}}\right)$$

$${\hat{w}}_{\bar{j},k}:=\frac{1}{2\pi ϵ}\mathrm{Arg}\left(\frac{\mathcal{F}\left({\mathbf{r}}_{p;ϵ}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]}\right),$$

$$\Vert \mathrm{Arg}(\gamma +\nu )-\mathrm{Arg}(\gamma ){\Vert }_{2\pi \mathbb{Z}}={\Vert \mathrm{Arg}\left(1+\frac{\nu }{\gamma }\right)\Vert }_{2\pi \mathbb{Z}}\le \frac{\pi }{2}\left|\frac{\nu }{\gamma }\right|,$$

$${\Vert \mathrm{Arg}\left(\frac{\mathcal{F}\left({\mathbf{r}}_{p;ϵ}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]}\right)-2\pi w_{\bar{j},k}^{\prime }ϵ\Vert }_{2\pi \mathbb{Z}}\le \mathcal{O}\left(\frac{\sigma }{|a_{\min }|\sqrt{p}}\right).$$

$${\Vert {\hat{w}}_{\bar{j},k}-w_{\bar{j},k}^{\prime }\Vert }_{\mathbb{Z}}\le \mathcal{O}\left(\frac{\sigma }{2\pi ϵ|a_{\min }|\sqrt{p}}\right),$$

$$a_{\bar{j}}=\frac{1}{p}\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]+\mathcal{O}\left(\frac{\sigma }{\sqrt{p}}\right).$$

$$w_{j,k}^0:=\frac{1}{2\pi {ϵ}_0}\mathrm{Arg}\left(\frac{\mathcal{F}\left({\mathbf{r}}_{p;{ϵ}_0}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]}\right),$$

$$b_{k,\alpha }:=\frac{1}{2\pi }\mathrm{Arg}\left(\frac{\mathcal{F}\left({\mathbf{r}}_{p;{ϵ}_{\alpha }}^{\tilde{k};k}\right)[m]}{\mathcal{F}\left({\mathbf{r}}_p^{\tilde{k}}\right)[m]}\right)$$

$$w_{j,k}^{\alpha }:=w_{j,k}^{\alpha -1}+\frac{\left(b_{k,\alpha }-{ϵ}_{\alpha }{w}_{j,k}^{\alpha -1}\right)\left(\mathrm{mod}\left[-\frac{1}{2},\frac{1}{2}\right)\right)}{{ϵ}_{\alpha }}$$

$$b_{k,\alpha }\approx {ϵ}_{\alpha }{w}_{j,k}^{\prime }\left(\mathrm{mod}\left[-\frac{1}{2},\frac{1}{2}\right)\right),$$

$${ϵ}_{\alpha }(w_{j,k}^{\prime }-w_{j,k}^{\alpha -1})$$

$$\approx (b_{k,\alpha }-{ϵ}_{\alpha }{w}_{j,k}^{\alpha -1})\left(\mathrm{mod}\left[-\frac{1}{2},\frac{1}{2}\right)\right)$$

$$\mathcal{O}\left(\frac{\sigma }{{ϵ}_{\alpha }{a}_{\min }\sqrt{p}}\right)$$

$$=w_{j,k}^{\prime }-\left(w_{j,k}^{\alpha -1}+\frac{\left(b_{k,\alpha }-{ϵ}_{\alpha }{w}_{j,k}^{\alpha -1}\right)\left(\mathrm{mod}\left[-\frac{1}{2},\frac{1}{2}\right)\right)}{{ϵ}_{\alpha }}\right)$$

$$=w_{j,k}^{\prime }-w_{j,k}^{\alpha }.$$

Full text

**Keywords ** Higher dimensional sparse FFT $\cdot$ Fast Fourier algorithms $\cdot$ Fourier analysis $\cdot$ Multiscale algorithms

**Mathematics subject classification ** 65T50 $\cdot$ 68W25

1 Introduction

Sparsity or compressibility in large data sets appears in many applications. Efficient algorithms taking advantage of these properties have been developed in different areas [1, 2, 5, 13] in order to reduce sampling and/or runtime complexities. Compressive sensing [6, 7] demonstrates that an appropriate small number of samples are sufficient to solve under-determined linear systems under certain well known conditions. More precisely, assuming $\boldsymbol{x}\in\mathbb{C}^{L}$ is $s$-sparse then $A\boldsymbol{x}=\boldsymbol{y}$ where $A\in\mathbb{C}^{I\times L}$, $\boldsymbol{y}\in\mathbb{C}^{I}$ and $I<L$ can be solved under appropriate conditions on $A$ and $I$. Compressive sensing has led to an explosion in the study of sparsity and algorithms that can take advantage of sparsity to drastically reduce both runtime and sample complexities. One example of the algorithms in this type is the sparse Fourier transform, which finds the energetic Fourier modes for a signal that is $s$-sparse in the Fourier domain quickly with reduced number of samples.

Several algorithms with varying approaches have been proposed to achieve sublinear runtime and sampling complexity for sparse Fourier transform, for both one-dimensional and higher-dimensional settings. In the one-dimension, the first sparse Fourier transform was a randomized algorithm introduced in [9] having $\mathcal{O}(s^{2}\log^{c}N)$ runtime and $\mathcal{O}(s^{2}\log^{c}N)$ sampling complexity with a small positive constant $c$. The constant $c$ controls the balance between accuracy and efficiency. In the follow up work [10] both runtime and sampling complexity were improved to $\mathcal{O}(s\log^{c}N)$. The randomized algorithms in [11, 12] have $\mathcal{O}(s\log N\log N/s)$ average-case runtime complexity. The first deterministic algorithm is introduced in [14], which uses techniques from number theory and combinatorics and has $\mathcal{O}(s^{2}\log^{4}N)$ runtime and sampling complexity. In [15], an improved deterministic algorithm was introduced, and the extension to higher dimensional function was suggested but it suffered the exponential dependence of runtime complexity on the dimension. Another deterministic algorithm was introduced in [17] which uses the similar idea in the frequency recovery through phaseshift and works for noiseless samples from exactly $s$-sparse functions. It has $\Theta(s\log s)$ average-case runtime and $\Theta(s)$ sampling complexity. In [4] the method from [17] was extended to incorporate a multiscale technique that works robustly with noisy samples and has $\Theta(s\log s\log N/s)$ average-case runtime and $\Theta(s\log N/s)$ sampling complexity.

Extension of one-dimensional spare Fourier transform to the multi-dimensional problem setting is in general not straightforward. One simple way of extension is to unwrap the multi-dimensional signal to one-dimensional signal. However, this method suffers from the exponentially large runtime complexity due to the curse of dimensionality [15]. The first randomized algorithm for two-dimensional problem was introduced in [8] through the use of parallel projections of frequencies. It has $\mathcal{O}(s)$ sampling and $\mathcal{O}(s\log s)$ runtime complexity on average for the exactly sparse signals, and $\mathcal{O}(s\log N)$ sampling and $\mathcal{O}(s\log^{2}N)$ runtime complexity on average for the approximately sparse signals. In [20] a general $d$-dimensional sparse Fourier algorithm was developed to achieve $\mathcal{O}(ds^{2}N)$ samples and $\mathcal{O}(ds^{3}+ds^{2}N\log(sN))$ runtime complexity. The algorithm uses rank-1 lattices and it finds energetic frequencies in a dimension-incremental fashion. While it is a deterministic algorithm it can also be modified into a randomized algorithm with $\mathcal{O}(ds+dN)$ sample complexity and $\mathcal{O}(ds^{3})$ runtime complexity. A randomized algorithm introduced in [16] requires $2^{\mathcal{O}(d^{2})}(s\log N^{d}\mathop{\rm loglog}N^{d})$ samples and $2^{\mathcal{O}(d^{2})}s\log^{d+3}N^{d}$ runtime. In [19], two deterministic sampling sets $\mathcal{O}(s^{2}d^{2}N)$ and $\mathcal{O}((s^{2}d^{3}N\log N)$ are constructed and the corresponding algorithms have $\mathcal{O}(s^{2}d^{2}N^{d})$ and $\mathcal{O}(s^{2}d^{3}N^{2}\log N)$ runtime complexities respectively under the assumption that $N$ is a prime number. In [3] we have developed an algorithm combining phaseshift from [15] and various transformations including parallel projection was given with average-case runtime complexity of $\mathcal{O}(sd\log s)$ and sampling complexity $\mathcal{O}(sd)$ under certain assumptions.

In this paper, we introduce a multi-dimensional sparse Fourier algorithm for noisy samples. Our algorithm uses the techniques from our algorithm for noiseless samples in [3] and the multiscale technique from [4] to overcome some of the challenges. Let $f:[0,1)^{d}\rightarrow\mathbb{C}$ be defined by $f(\boldsymbol{x})=\sum_{j=1}^{s}a_{j}e^{2\pi i\boldsymbol{w}_{j}\cdot\boldsymbol{x}}+n(\boldsymbol{x})$ where $\boldsymbol{w}_{j}\in[-N/2,N/2)^{d}\cap\mathbb{Z}^{d}$, $s\ll N^{d}$ and $n(\boldsymbol{x})$ represents noise. The algorithms from [3] are not robust to noise since it needs to compute the ratio of discrete Fourier transforms of samples from $f$ at shifted and unshifted points, which is not robust when the shift is small or the noise level $\sigma$ is high. To overcome this we adopt a multiscale approach similar to the one introduced in [4]. This allows us to progressively approximate the significant modes while controlling the influence of noise. In higher dimensions there are some additional challenges we are able to overcome. Details will be discussed later in the paper.

Our algorithm assumes that we have access to an underlying continuous function $f$, i.e., we can sample at anywhere we want. However, samples are sometimes given at the beginning and getting extra samples can be very expensive. Accordingly, it is necessary to do approximation of extra samples using given samples. A fully discrete sparse Fourier transform was introduced in [18] combining periodized Gaussian filters and one-dimensional sparse Fourier transform in continuous setting such as [15] and [4]. For future work we are hopeful that similar approaches can help us to adopt the multiscale high-dimensional algorithm to the case of fully discrete samples.

The rest of this paper is organized as follows. In Section 2 we introduce our problem setting, necessary notation and our noise model, and review briefly about high-dimensional sparse Fourier algorithm from [3]. Section 3 introduces a multiscale method for the high-dimensional sparse Fourier algorithm from our previous work. This modifies the algorithm to be able to recover noisy signals. In Section 4, parameters that determine the performance of our algorithm are introduced and the pseudocode is given with the description and analysis. The results of the numerical experiments are shown in Section 5 and the conclusion is in Section 6.

2 Preliminaries

2.1 Notation and Review

In this section, we introduce the notation used throughout the rest of this paper and the brief review of high-dimensional sparse Fourier algorithm for samples without noise in [3]. Let $s$, $d$ and $N$ be natural numbers, $s\ll N^{d}$, and $D:=[0,1)^{d}$. We consider a function $f:D\rightarrow\mathbb{C}$ which is $s$-sparse in the $d$-dimensional Fourier domain as follows

[TABLE]

where each $\boldsymbol{w}_{j}\in[-N/2,N/2)^{d}\cap\mathbb{Z}^{d}$ and $a_{j}\in\mathbb{C}$. We note that $f$ can be regarded as a periodic function defined on $\mathbb{R}^{d}$. The aim of sparse Fourier algorithms is to rapidly reconstruct a function $f$ using small number of its samples. In [3], the methods were introduced using several different transformations and parallel projections along coordinate axes of frequencies in order to exploit the one-dimensional sparse Fourier algorithm from [17]. These transformations such as partial unwrapping and tilting methods are introduced in order to change the locations of energetic frequencies when the current energetic frequencies are hard or impossible to find through the parallel projections directly. Through those manipulations, each frequency vector $\boldsymbol{w}_{j}$ is recovered in an entry-wise fashion. In this way, the linear dependence of runtime and sampling complexities on the dimension $d$ could be shown empirically, which is a great improvement when compared to the $d$-dimensional FFT with the exponential dependence on $d$. The transformations and projections occur in the physical domain, which provides the separation of the frequency vectors in the Fourier domain. That is, each $\boldsymbol{w}_{j}$ is transformed to $\boldsymbol{w}_{j}^{\prime}$ and then projected onto several lines, and these can be done by manipulating the sampling points in the physical domain $D$. Let $g:D^{\prime}\rightarrow D$ represent those transformations where $D^{\prime}$ is $d$ or less dimensional space and is determined by each transformation. We assume that $D^{\prime}$ has $d^{\prime}$ dimensions. A new function $h:D^{\prime}\rightarrow\mathbb{C}$ is defined as a composition of $f$ and $g$, i.e.,

[TABLE]

We note that $h$ is still $s$-sparse in the $d^{\prime}$-dimensional Fourier domain, $[-N^{\prime}/2,N^{\prime}/2)^{d^{\prime}}\cap\mathbb{Z}^{d^{\prime}}$, whose bandwidth $N^{\prime}$ depends on each transformation. For example, consider a 4-dimensional function $f$ with the Fourier domain, $[-N/2,N/2)^{4}\cap\mathbb{Z}^{4}$. If a partial unwrapping is applied to $f$ which unwraps each $\boldsymbol{w}_{j}=(w_{j,1},w_{j,2},w_{j,3},w_{j,4})$ to $\boldsymbol{w}_{j}^{\prime}=(w_{j,1}+Nw_{j,2},w_{j,3}+Nw_{j,4})$ then it implies that $g:(t_{1},t_{2})\rightarrow(t_{1},Nt_{1},t_{2},Nt_{2})$ and $h$ has the Fourier domain, $[-N^{2}/2,N^{2}/2)^{2}\cap\mathbb{Z}^{2}$ where accordingly $N^{\prime}=N^{2}$ and $d^{\prime}=2$. This is one example and there are variations of partial unwrapping methods and tilting methods which can be found in [3]. Using samples of $h$(or $f$), the transformed frequency vectors $\boldsymbol{w}_{j}^{\prime}$ and corresponding Fourier coefficients $a_{j}$ are found through the parallel projection method and $\boldsymbol{w}_{j}^{\prime}$ are transformed back to $\boldsymbol{w}_{j}$. Now, we introduce how $\boldsymbol{w}_{j}^{\prime}=(w_{j,1},w_{j,2},\cdots,w_{j,d^{\prime}})$ can be recovered element-wisely using the parallel projection method and the ideas from one-dimensional sparse Fourier algorithm. Let $p$ be a prime number greater than a constant multiple of $s$, i.e., $p>cs$ for some constant $c$, $j$ be a fixed integer among $\{1,2,\cdots,s\}$, $k$ be fixed among the set $\{1,2,\cdots d^{\prime}\}$ and $\boldsymbol{e}_{k}$ be a vector with all zero entries but $1$ at the index $k$. Furthermore, $\epsilon$ is defined as a positive number $\leq 1/N^{\prime}$. To recover $k$-th element of each $\boldsymbol{w}_{j}^{\prime}$, we use two sets of $p$-length equispaced samples $\mathbf{h}_{p}^{\widetilde{k}}$ and $\mathbf{h}_{p;\epsilon}^{\widetilde{k};k}$ as follows,

[TABLE]

where $\ell=0,1,\cdots,p-1$ and $\widetilde{k}\in\{1,2,\cdots,d^{\prime}\}$ is the index of coordinate axis where a particular $\boldsymbol{w}^{\prime}_{\bar{j}}$ has $\widetilde{k}$-th element, $w^{\prime}_{\bar{j},\widetilde{k}}$, different from the $\widetilde{k}$-th elements of any other energetic frequency vectors. In this case, we refer to the above phenomenon as “no collision from projection”. At the same time, if there is no collision modulo $p$, i.e., $w^{\prime}_{\bar{j},\widetilde{k}}$ has the unique remainder modulo $p$ from others then the discrete Fourier transform of each sample set is

[TABLE]

respectively. The equations above in (2.1) give a unique entry for the $k^{\rm th}$ element assuming there does not exist a collision modulo $p$ of the vetor projected onto the $k^{\rm th}$ axis. Hence, in the above equations, when the second equalities hold, we can recover $a_{\bar{j}}$ and $w^{\prime}_{\bar{j},k}$ as follows,

[TABLE]

where the function ${\rm Arg}(z)$ is defined to be the argument of $z$ in the branch $[-\pi,\pi)$. The right choice of the branch and the shift size $\epsilon\leq\frac{1}{N^{\prime}}$ make it possible to find the correct $w^{\prime}_{\bar{j},k}$. Algorithmically, the two kinds of collisions are guaranteed not too happen using the following test:

[TABLE]

which is inspired by the test used in [17]. Practically, we put some threshold $\tau>0$ so that if the difference between that the left and right-hand sides is less than $\tau$, then we conclude that there are no collisions of both kinds. In this case, each $w_{\bar{j},k}$ for $k=1,2,\cdots,d^{\prime}$ can be recovered using a pair of sets $\mathbf{h}_{p}^{\widetilde{k}}$ and $\mathbf{h}_{p;\epsilon}^{\widetilde{k};k}$, respectively. Otherwise, we take another prime number for sample length, switch the index of coordinate axis for the projection, update the samples by eliminating the influence from previously found fourier modes and repeat our procedure as before. Switching the coordinate axis and updating samples reduce the occasions of collision from projection. In [17], moreover, it is proved that the probability is very low when $d^{\prime}$ is large that all remaining frequency vectors have collisions from projection onto all coordinate axes, which we call the “worst case scenario”. In the “worst case scenario”, the parallel projections we just used do not work and thus, we need a rotation mapping, $g^{\prime}$, defining another function $h^{\prime}=f\circ g^{\prime}$.

2.2 Noise Model

In this section, we introduce a model system which contains noise. We use this model system to quantify the behavior of the algorithm in the presence of noise. The algorithm we introduced in Section 2.1 works well when $h$ (or $f$) is exactly $s$-sparse and the samples from $h$ (or $f$) are not noisy. The high-dimenional sparse Fourier transform, described in the previous section, is not robust to noise since in order to find entries of energetic frequency vectors we compute the fraction {\mathcal{F}\left(\mathbf{h}_{p;\epsilon}^{\widetilde{k};k}\right)[m]}\big{/}{\mathcal{F}\left(\mathbf{h}_{p}^{\widetilde{k}}\right)[m]} which is sensitive to noise. The model we consider here is,

[TABLE]

where $h$ is the tranformed $f$ defined in the previous section, $\ell=0,1,\cdots,p-1$, and $\boldsymbol{n}:=(n_{0},n_{1},\cdots,n_{p-1})$ is a complex Gaussian random variable with mean $\boldsymbol{0}$ and variance $\sigma^{2}I$. If we apply DFT to this sample set, we get

[TABLE]

Since $n_{\ell}$ are i.i.d Gaussian variables, the expectation and variance of the second term in (2.4) are

[TABLE]

and

[TABLE]

respectively. Accordingly,

[TABLE]

and

[TABLE]

For noisy shifted sample set $\mathbf{r}_{p;\epsilon}^{\widetilde{k};k}:=\mathbf{h}_{p;\epsilon}^{\widetilde{k};k}+\widetilde{\boldsymbol{n}}$ with an i.i.d Gaussian random vector $\widetilde{\boldsymbol{n}}$, we have likewise

[TABLE]

and

[TABLE]

In the case of $w^{\prime}_{\bar{j},\widetilde{k}}$ not having collisions both from projection and modulo $p$, and $w^{\prime}_{\bar{j},\widetilde{k}}\equiv m~{}(\bmod~{}p)$,

[TABLE]

and

[TABLE]

for each $k=1,2,\cdots,d^{\prime}$. As a result, we get

[TABLE]

and note that if there were no noise in samples, we only have the first term on the right side of (2.6) which makes it possible to recover $w^{\prime}_{\bar{j},k}$ by taking its argument and dividing it by $2\pi\epsilon$ as (2.2). With noisy samples, however, it is corrupted with noise which is a multiple of $\frac{\sigma}{a_{\bar{j}}\sqrt{p}}$. Defining

[TABLE]

we want to see how far $\hat{w}_{\bar{j},k}$ is from $w^{\prime}_{\bar{j},k}$. For this purpose, we introduce the Lee norm associated with a lattice $\mathcal{L}$ in $\mathbb{R}$ as $\|z\|_{\mathcal{L}}:=\min_{y\in\mathcal{L}}|z-y|$ for $z\in\mathbb{R}$ and the related property that under the Lee norm associated with the lattice $2\pi\mathbb{Z}$,

[TABLE]

where $|\gamma|\geq|\nu|$ with $\gamma,\nu\in\mathbb{C}$. By choosing the sample length $p$ large enough depending on the least magnitude nonzero $a_{min}$ and the noise level $\sigma$, (2.8) can be applied to (2.6) as follows,

[TABLE]

Consequently,

[TABLE]

which implies that the error of our estimate $\hat{w}_{\bar{j},k}$ to $w^{\prime}_{\bar{j},k}$ is controlled by the size of $\frac{\sigma}{\epsilon~{}|a_{\min}|\sqrt{p}}$. Thus, $p$ needs to be chosen carefully depending on $\frac{\sigma}{\epsilon~{}|a_{\min}|}$. On the other hand, from (2.5), we can approximate the corresponding coefficient $a_{\bar{j}}$ with the error of size $\mathcal{O}(\sigma/\sqrt{p})$ as follows

[TABLE]

3 Multiscale Method

Here we introduce the multiscale approach for recovering frequencies in the noisy setting. The method was introduced in [4]. The method is part of the overall sublinear algorithm described in Section 4. The basic idea of the algorithm is to find the most significant bits which are the least susceptible to noise. Then the algorithm subtracts the leading bits and shifts the remaining bits to the most significant digits. As we describe here in Section 3, this will decrease the impact of the noise by the use of larger shifting size $\epsilon_{\alpha}$ recovering the next most significant bits. This is repeated to recover the entire entries of the frequency vector.

In [4], rounding and multiscale methods for the one-dimensional sparse Fourier algorithm for noisy data were introduced. Both methods use the fact that the peaks of DFT are robust to relatively high noise, i.e., $\ell$ can be correctly found such that $w\equiv\ell(\bmod~{}p)$ for an energetic frequency $w$. The rounding method is efficient when $\sigma$ is relatively small, which approximates such $w:=bp+\ell$ for some $b$ up to $p/2$ error and rounds a multiple of $p$ in order to get the correct $b$. It was shown that $p\geq\max\{c_{1}s,c_{2}(\frac{\sigma}{\epsilon a_{\min}})^{2/3}\}$ for some constants $c_{1}$ and $c_{2}$ makes it possible to correctly find $w$ through the rounding method. On the other hand, the multiscale method was introduced for relatively large $\sigma$. It prevents $p$ from becoming too large, which happens for large $\sigma$ in the rounding method. In this section, we focus on extending the multiscale method to recover high dimensional frequencies by gradually fixing each entry estimation with several shifts $\epsilon_{\alpha}$.

3.1 Description of Frequency Entry Estimation

In this section, we give an overview of how the multiscale method works. Let ${j}\in\{1,2,\cdots,s\}$ and $k\in\{1,2,\cdots,d^{\prime}\}$ be fixed. The target frequency entry $w^{\prime}_{j,k}$ is assumed not to have collisions both from projection and modulo $p$. We start with a coarse estimation $w^{0}_{j,k}$ of $w^{\prime}_{j,k}$ defined by

[TABLE]

where $\epsilon_{0}\leq 1/N^{\prime}$. Then $w^{0}_{j,k}\equiv w^{\prime}_{j,k}(\bmod~{}p)$ even though it is not guaranteed that $w^{0}_{j,k}=w^{\prime}_{j,k}$. Thus, we need to improve the approximation. With each correction, the solution is improved by $l$ digits where $l$ depends on the parameters that are chosen in the method as well as noise. Each correction term is calculated with a choice of growing $\epsilon_{\alpha}>1/N^{\prime}$, i.e., $\epsilon_{\alpha-1}<\epsilon_{\alpha}$ for all $\alpha\geq 1$. With the initialization $b_{k,0}:=\epsilon_{0}w^{0}_{j,k}$, the correction terms are calculated in the following way,

[TABLE]

and

[TABLE]

for $\alpha\geq 1$, where $x~{}(\bmod~{}[-1/2,1/2))$ is defined to be the value $a\in[-1/2,1/2)$ such that $x\equiv a~{}(\bmod~{}1)$. Using the fact,

[TABLE]

the error $w^{\prime}_{j,k}-w^{\alpha-1}_{j,k}$ can be estimated as

[TABLE]

and thus, in a similar manner to (2.9),

[TABLE]

As the correction is repeated with larger $\epsilon_{\alpha}$, the error of the estimate decreases. In other words, we approximate $w^{\prime}_{j,k}$ by its most significant bits and the next significant bits repeatedly. The performance of this multiscale method is shown in detail in the next section.

3.2 Analysis of Multiscale Method

In this section, we establish the multiscale algorithm recovering a fixed number of additional bits of the frequency with each iteration, and we further establish the rate of reconstruction. The following theorem shows how the correction term $c_{k,\alpha}/\epsilon_{\alpha}$ is constructed in each iteration and how large the error of estimate $w^{\alpha}_{j,k}$ after $M$ iterations is in the multiscale frequency entry estimation procedure.

Theorem 1**.**

Let $j\in\{1,2,\cdots,s\},k\in\{1,2,\cdots,d^{\prime}\}$ be fixed and $w^{\prime}_{j,k}\in\left[-\frac{N^{\prime}}{2},\frac{N^{\prime}}{2}\right)$. Let $0<\epsilon_{0}<\epsilon_{1}<\cdots<\epsilon_{M}$ and $b_{k,0},b_{k,1},\cdots,b_{k,M}\in\mathbb{R}$ such that

[TABLE]

where $0<\delta<\frac{1}{4}$. Assume that $\epsilon_{0}\leq\frac{1-2\delta}{N^{\prime}}$ and $\beta_{\alpha}:=\frac{\epsilon_{\alpha}}{\epsilon_{\alpha-1}}\leq\frac{1-2\delta}{2\delta}$. Then there exist $c_{k,0},c_{k,1},\cdots,c_{k,M}\in\mathbb{R}$, each computable from $\{\epsilon_{\alpha}\}$ and $\{b_{k,\alpha}\}$ such that

[TABLE]

Proof.

The proof is the same as the proof of Theorem 4.2 in [4] since each entry of frequency vectors is corrected in the same way as each one-dimensional frequency $w$ in [4] is. Each $c_{k,\alpha}$ is defined as

[TABLE]

for $\alpha\geq 1$ where $\lambda_{k,\alpha}=c_{k,\alpha}/\epsilon_{\alpha}$ for $\alpha\geq 0$. ∎

Corollary 1**.**

Assume that we let $\beta_{\alpha}=\beta$ in Theorem 1 where $\beta\leq(1-2\delta)/(2\delta)$, i.e., $\epsilon_{\alpha}=\beta^{\alpha}\epsilon_{0}$ for all $\alpha\geq 1$. Let $p>0$ and $M\geq\left\lfloor\log_{\beta}\frac{2\delta}{\epsilon_{0}}\right\rfloor+1$. Then,

[TABLE]

Proof.

This is straightforward corollary of Theorem 1. ∎

Corollary 1 looks very similar to Corollary 4.3 in [4]. However, it is different in that the iteration number is increased in order to make the error between $\widetilde{w}_{j,k}$ and $w^{\prime}_{j,k}$ is less than $1/2$ instead of $p/2$. In [4], the remainder $m$ modulo $p$ of each energetic $w$ is known by sorting out the $s$ largest DFT components of $p$-length unshifted sample vector so that $p/2$ error bound is enough to guarantee the exact recovery of $w$. On the other hand, in the high-dimensional setting of this paper, the remainders of all entries of $w^{\prime}_{j,k}$ are not known. Instead, the remainder of $w^{\prime}_{j,\widetilde{k}}$ is only known where $\widetilde{k}$ is the index of coordinate axis where the frequencies are projected. Thus, we decrease the error further by enlarging the iteration number $M$ and are able to recover the exact $w^{\prime}_{j,k}$ by rounding when the error is less $1/2$.

Remark 1**.**

Similar to Theorem 4.4 in [4], the admissible size of $\delta$ can be estimated as follows

[TABLE]

under the assumption $\epsilon_{0}\leq\frac{1-2\delta}{N^{\prime}}$ of Corollary 1 together with the assumption $\beta\leq\frac{1-2\delta}{2\delta}$ of Theorem 1.

4 Algorithm

In this section, we present the overall multiscale high-dimensional sublinear sparse FFT which combines the multiscale method from Section 3 with our previous work in [3]. In addition, we describe how key parameters of the algorithm are chosen. As will be discussed below, the choice of parameters are affected by the noise level, $\sigma$. It is important to know how the frequency entry estimation works in the algorithm and how the collision detection tests are modified from the tests in [3] in order to make them tolerant of noise. Furthermore, we present the analysis of the average-case runtime and sampling complexity under assumption that the worst case scenario does not happen.

4.1 Choice of $p$

In this section, we establish the length of subsampling vector, $p$. The sample length $p$ affects the total runtime complexity since the discrete Fourier transform is applied to all sample sets taken to recover the frequency entries. Due to this, we want to make it as small as possible. At the same time, however, we can see from (2.9) that the error between the target entry and its approximation becomes smaller if a larger $p$ is taken. Thus, as discussed in Section 3, if $p$ is large enough, then the rounding method instead of multiscale method can recover the exact frequency by rounding (2.7) to the nearest integer of the form $pv+m$ with an integer $v$. In this case, $p$ is large enough to diminish the influence of $\sigma$, and therefore if $\sigma$ is large, so is $p$. Instead, the multiscale method makes it possible to enlarge $p$ moderately. From Theorem 1, we get

[TABLE]

and (3.3) implies

[TABLE]

By putting the right side of (4.2) as $c_{\sigma}\frac{\sigma}{\epsilon_{\alpha}~{}|a_{\min}|\sqrt{p}}$ with some constant $c_{\sigma}$ and equating both right sides of (4.1) and (4.2), $\beta:=\epsilon_{\alpha+1}/\epsilon_{\alpha}$ can be estimated as

[TABLE]

$\beta$ determines the choice of $\epsilon_{\alpha}$ for each iteration, i.e., $\epsilon_{\alpha}=\beta^{\alpha}\epsilon_{0}$ where $\epsilon_{0}$ is chosen to be less than $1/N^{\prime}$. Combining (3.4) and (4.3), the sample length $p$ can be calculated as

[TABLE]

when $\epsilon_{0}=\frac{1}{2N^{\prime}}$ and $\beta>1$, which implies $\delta=\frac{1}{2\beta+2}$. Eventually, the sample length $p$ for the multiscale method needs to satisfy

[TABLE]

where $p>c_{1}s$ ensures the sample is long enough so that the $90\%$ of all energetic frequencies are not collided modulo $p$ on average which comes from the pigeonhole argument in [17].

4.2 Collision Detection Tests

As mentioned in Section 2.1, frequencies are recovered only when there are no collisions from the projection and the modulo $p$ division. These conditions are satisfied if

[TABLE]

for $k=1,\cdots,d^{\prime}$ and some small $\tau>0$ which are the practical tests of (2.3). In our noisy setting, Equation (2.6) implies that the left hand side of (4.5) is bounded above by $\mathcal{O}(\frac{\sigma}{a_{\min}\sqrt{p}})$. Thus, we set our threshold $\tau$ as a constant multiple of $\frac{\sigma}{a_{\min}\sqrt{p}}$. Moreover, since we iteratively update the estimates $w^{\alpha}_{j,k}$ for $\alpha=0,1,\cdots,M$, we reject the estimate after $M$ iterations if the tests fail for more than $\eta(M+1)$ times for each $k$-th entry with $k=1,2,\cdots,d^{\prime}$ where $\eta$ is a fraction$<1$. Numerical experiments indicate $\eta=\frac{1}{4}$ is a good number.

4.3 Number of Iterations

In this section, we give a specific choice for the number of iterations in our multiscale algorithms. From Corollary 1, $M\geq\left\lfloor\log_{\beta}\frac{2\delta}{\epsilon_{0}}\right\rfloor+1$ guarantees we get the approximation error, $\left|\widetilde{w}_{j,k}-w^{\prime}_{j,k}\right|<\frac{1}{2}$ which is required to recover the exact $w^{\prime}_{j,k}$ by rounding $\widetilde{w}_{j,k}$ to the nearest integer. With our choice of $\epsilon=\frac{1}{2N^{\prime}}$ and the fact that $\delta<1$, $M=\left\lfloor\log_{\beta}N^{\prime}\right\rfloor+1$ suffices to satisfy the $1/2$ error bound. For example, if each $\boldsymbol{w}_{j}\in\left[-\frac{N}{2},\frac{N}{2}\right)^{d}\cap\mathbb{Z}^{d}$ is partially unwrapped to some value in $\left[-\frac{N^{d_{1}}}{2},\frac{N^{d_{1}}}{2}\right)^{d_{2}}\cap\mathbb{Z}^{d_{2}}$ where $d_{1}$ and $d_{2}$ are positive integers satisfying $d=d_{1}d_{2}$, then $d_{2}$ entries of each energetic frequency are recovered element-wisely after $M=\mathcal{O}(d_{1}\log N)$ iterations.

4.4 Description of Our Pseudocode

In this section, we explain the multiscale high-dimensional sparse Fourier transform whose pseudocode is provided in Algorithm 1. The set $R$ contains the identified Fourier frequencies and their corresponding coefficients, and it is an empty set initially. Parameter $i$ is the counting number determining the index $\widetilde{k}$ of the coordinate axes where the frequencies are projected in line 6. Parameters $p,\tau$ and $M$ are determined as discussed in the previous sections. Function $q$ in line 7 is a function constructed from the previously found Fourier modes which is used in updating our samples in lines 9 and 17. In line 9, the unshifted samples $\boldsymbol{r}^{\widetilde{k}}_{p}$ corrupted by random Gaussian noise $n_{\ell}$ are taken and in line 11, DFT is applied to these samples and the transformed vector is sorted in the descending order of magnitude. Only its $s^{\ast}$ largest components are taken into account under the $s^{\ast}$-sparsity assumption. In the loop from line 12 through 39, the entries of frequencies corresponding to these $s^{\ast}$ components are estimated iteratively. The shift size $\epsilon_{\alpha}$ is updated in line 14 and we get the $d^{\prime}$ number of length $p$ samples at the points shifted by $\epsilon_{\alpha}$ along each axis. Each length-$p$ sample will be used to approximate each entry. Similar to line 11, DFT is applied to each length-$p$ sample and the transformed vector is sorted again following the index order of the sorted DFT of unshifted samples in line 19. In line 22, we check whether the $d^{\prime}$ tests are passed at the same time or not. If not, the $vote$ is increased by 1. From lines 24 through 34, the estimate $w_{j,k}^{\prime}$ for $k$th entry is updated. Except when $\alpha=0$, $w_{j,k}^{\prime}$ is improved by adding the correction term shown in line 29 in each iteration. In the last iteration when $\alpha=M$, $w_{j,k}^{\prime}$ is rounded to the nearest integer in order to recover the exact entry, as guaranteed by Corollary 1. Whether this estimate is stored in the set $R$ or not is determined by checking if the $vote$ after $M$ iterations is less than $\eta(M+1)$. If it is less, this implies that the failure rate of the collision detection tests is less than $\eta$. Accordingly, we estimate $a_{\boldsymbol{w}^{\prime}_{j}}$ from the DFT of unshifted samples and store $(\boldsymbol{w}^{\prime}_{j},a_{\boldsymbol{w}^{\prime}_{j}})$ in $R$. The entire while loop repeats until $s$ energetic Fourier modes are all found switching the projection coordinate. Once we find all $s$ frequency vectors, each $\boldsymbol{w}^{\prime}_{j}\in R$ is transformed to $d$-dimensional $\boldsymbol{w}_{j}=g^{-1}(\boldsymbol{w}^{\prime}_{j})$.

4.5 Runtime and Sampling Complexity for the Average Case Signals Under No Worst-case Scenario Assumption

In this section, we explain the performance of the multiscale high-dimensional sparse Fourier algorithms. In particular, we will restrict ourselves to the average-case analysis under the assumption that the wort-case scenario does not happen.

Theorem 2**.**

Let $f^{z}(\boldsymbol{x})=f(\boldsymbol{x})+z(\boldsymbol{x})$, where $\hat{f}(\boldsymbol{w})$ is $s$-sparse with each frequency $\boldsymbol{w}_{j}\in[-N/2,N/2)^{d}\cap\mathbb{Z}^{d}$ corresponding to the nonzero Fourier coefficient for $j\in\{1,2,\cdots,s\}$ and not forming any worst case scenario, and $z$ is complex i.i.d. Gaussian noise of variance $\sigma^{2}$. Moreover, suppose that $s>C(\beta(\beta+1)c_{\min}\sigma)^{2}$ for some constant $C$. Algorithm 1, given $N,d,s,\beta$ with $N>5s$ and access to $f^{z}(\boldsymbol{x})$ returns a list of $s$ pairs $(\hat{\boldsymbol{w}},a_{\hat{\boldsymbol{w}}})$ such that (i) each $\hat{\boldsymbol{w}}=\boldsymbol{w}_{j}$ for some $j\in\{1,2,\cdots,s\}$ and (ii) for each $\hat{\boldsymbol{w}}$, $|a_{\boldsymbol{w}_{j}}-a_{\hat{\boldsymbol{w}}}|\leq C\sigma/\sqrt{{s}}$. The average-case runtime and sampling complexity are

[TABLE]

respectively, over the class of random signals.

Proof.

The difference of Algorithm 1 in this paper from Algorithm 2 in [3] appears in lines 13 through 39. Algorithm 1 has the multiscale frequency entry estimation. Thus, the average-case runtime and sampling complexity of Algorithm 2 from [3] is increased by a factor of $M$ which is the number of the repetition in the multiscale frequency entry estimation from Section 4.3. Corollary 1 ensures that the returned frequency vectors $\hat{\boldsymbol{w}}$ are correct, and the coefficient $a_{\hat{\boldsymbol{w}}}$ has the desired error bound from (2.10). ∎

5 Empirical Evaluation

In this section, we show the empirical evaluation of the multiscale high-dimensional sparse Fourier algorithm. The empirical evaluation was done for test functions $f(\boldsymbol{x})$ which consist of $a_{j}$ randomly chosen from a unit circle in $\mathbb{C}$ and $\boldsymbol{w}_{j}$ randomly chosen from $[-N/2,N/2)^{d}\cap\mathbb{Z}^{d}$. Sparsity $s$ varied from $1$ to $2^{10}=1024$ by factor of 2. The noise term added to each sample of $f$ came from the Gaussian distribution $\mathcal{N}(0,\sigma^{2})$. Standard deviation $\sigma$ varied from $0.001$ to $0.512$ by factor of 2. The dimension $d$ was chosen to be $100$ and $1000$, and $N$ was chosen as $20$. The transformation $g$ was the one for partial unwrapping that was used in [3], i.e., every $5$-dimensional subvector of each frequency vector was unwrapped to a one-dimensional vector and therefore each $100$ and $1000$-dimensional function $f$ was unwrapped to $20$ or $200$-dimensional function $f\circ g$ whose Fourier domain is $[-20^{5}/2,20^{5}/2)^{20}\cap\mathbb{Z}^{20}$ or $[-20^{5}/2,20^{5}/2)^{200}\cap\mathbb{Z}^{200}$, respectively. The input parameters $c_{1}=2$, $c_{\sigma}=6$, $\eta=1/4$ and $\beta=2.5$ were empirically chosen to balance the runtime and accuracy as in [4]. The initial shift size $\epsilon_{0}$ was set to $\frac{1}{2\cdot 20^{5}}$. All experiments are performed in MATLAB.

The three plots in Figure 1 show the average over 10 trials of the $\ell_{1}$ error, the number of samples, and the runtime in seconds as the noise level $\sigma$ changes. These values are in logarithm in the plots. Dimension $d$ and sparsity $s$ are fixed with $100$ and $256$, respectively. On the other hand, the other three plots in Figure 2 show the average over 10 trials of the $\ell_{1}$ error, the number of samples, and the runtime in seconds as the sparsity $s$ changes when $d=100$ and $1000$. These values are in logarithm in the plots, and the noise level is fixed to $0.512$.

5.1 Accuracy

In this section, we do numerical experiments to investigate the accuracy of the algorithm. In Figure 1(a) and Figure 2(a), the $\ell_{1}$ errors of Fourier coefficient vectors are given under various parameter changes. Throughout all trials conducted in these experiments frequencies were always recovered exactly even for the noise level $\sigma=0.512$ which is relatively large compared to the true coefficients from the unit circle in $\mathbb{C}$. Thus, we can observe the errors only from coefficients whose size is a constant multiple of $\frac{\sigma}{\sqrt{p}}$ from (2.10). Due to the characteristic of the multiscale method which uses less samples compared to the rounding method, $\ell_{1}$ errors are relatively large in nature. From Figure 1(a), $\ell_{1}$ error looks increasingly linear as $\sigma$ increases, which meets our expectation. In Figure 2(a), the plot does not look exactly linear, but between $\log_{2}s=5$ and $6$, there is a transition of slope. This is because the sample length $p$ from (4.4) changes from $c_{1}s$ to $\left(\frac{\beta(\beta+1)a_{\min}c_{\sigma}\sigma}{\pi}\right)^{2}$ during this transition.

5.2 Sampling complexity

In this section, we numerically explore the sampling complexity of the algorithm. Sample numbers along $\sigma$ changes in Figure 1(b) seems irregular at first sight. Looking at the scale of vertical axis, however, we can see that the difference between maximum and minimum is less than $0.3$. Therefore, sampling complexity is not very affected by noise level. In Figure 2(b), the red graph shows the average sample numbers as the sparsity increases when $d=100$ and the blue graph shows the ones when $d=1000$. Since our multiscale algorithm recovers each frequency entry iteratively using $\log N$ sets of $\mathcal{O}(s)$-length samples, the average-case sampling complexity is indeed $\mathcal{O}(sd\log N)$ when the worst-case scenario does not happen. Two graphs in Figure 2(b) look close to be linear excluding the transition between $\log_{2}s=5$ and $6$, which again is caused by the change of $p$ from (4.4). Moreover the difference between the values of the red and blue graphs are close to $3$, which implies that the sampling number depends linearly on $d$. The $d$-dimensional FFT whose sampling complexity is $\mathcal{O}(N^{d})$ cannot deal with our high-dimensional problem computationally, whereas our algorithm uses only millions to billions of samples for reconstruction.

5.3 Runtime complexity

In this section, we consider the average-case runtime complexity of our multiscale high-dimensional sparse Fourier transform. Figures 1(c) and 2(c) demonstrate the average-case runtime complexity of the algorithm. The time for evaluating the samples from functions is excluded when measuring the runtime. For the main algorithm, we demonstrated that it is $\mathcal{O}(sd\log s\log N)$ because for each entry recovery, DFT with $\mathcal{O}({s\log s})$ runtime complexity is applied in $d\log N$ iterations. In Figures 1(c), the runtimes in seconds look irregular but the scale of vertical axis is less than $0.3$ so that we can conclude that similar to sample numbers the runtime is not affected by $\sigma$ very much. Overall, it took less than a second on average. In Figure 2(c), the red graph represents the runtimes as the sparsity changes when $d=100$, and the blue graph represents the ones when $d=1000$. Those graphs do not look linear, but considering the average slope we can see that the runtime is increased by around $2^{8}$ while $s$ is increased by $2^{10}$. On the other hand, the difference between two graphs implies that the runtime complexity is linear in $d$. Compared to the FFT with runtime complexity of $\mathcal{O}(N^{d}\log N^{d})$ which is impossible to be practical in high-dimensional problem, our algorithm is quite effective, taking only a few seconds.

6 Conclusion

In this paper, we developed a multiscale high-dimensional sparse Fourier algorithm recovering a few energetic Fourier modes using noisy samples. As the estimation error is controlled by the noise level $\sigma$ and the sample length $p$, larger $p$ reduces the error. Rather than recovering the frequencies in a single step, however, we choose multiscale approach in order to make the sample length $p$ increase moderately by improving the estimate iteratively through correction terms determined by a sequence of shifting sizes $\epsilon_{\alpha}$. We showed that a finite number of correction terms are enough to make the error smaller than $1/2$ so that we can reconstruct each integer frequency entry by rounding. As a result, the algorithm has $\mathcal{O}(sd\log N)$ sampling complexity and $\mathcal{O}(sd\log s\log N)$ runtime complexity on average under the assumption that there is no worst case scenario happening by combining the result from [3] and [4]. In the numerical experiment we ran, with a noise of $\sigma=0.512$ we were able to recover 100% of the frequencies up to dimension 1000. The methods introduced either in [3] and in this paper assume that we can get the measurement at any sample point. However, this is not always the case in practice. Our future work will be a modification of the algorithm to make it work for the given discrete signals using the idea of filtering from [18].

ACKNOWLEDGEMENTS We would like to thank Mark Iwen for his valuable advice. This research is supported in part by AFOSR grants FA9550-11-1-0281, FA9550-12-1-0343 and FA9550-12-1-0455, NSF grant DMS-1115709, and MSU Foundation grant SPG-RG100059, as well as Hong Kong Research Grant Council grants 16306415 and 16317416.