Algorithm for Functions with Short Frequency Support II

and a sampling complexity of 49 instead of a single CDFT of length1,000. Note that for this example the second method for reconstructing a function with short frequency support needs significantly less samples than the 220 samples required by the procedure described in Section 2.3. We will compare the runtimes and sampling complexities of the two methods more thoroughly in Section 2.5.

♦

2.4.2 Algorithm for Functions with Short Frequency Support II

We summarize the procedure introduced above in Algorithm 3, which deterministically finds the energetic frequencies and the Fourier coefficients of a function f ∈ C2π with bandwidth N and short frequency support if N and an upper bound B on the support length are known a priori. We will investigate its performance with respect to runtime and noisy input data in numerical experiments in Section 2.5, where we will also compare Algorithm 3 to Algorithm 2 and other sparse FFT methods.

As in Section 2.2.1, the function extended_gcd in line 13 of Algorithm 3 finds a representation of the greatest common divisor gof two integers aand bof the form

g= gcd(a, b) =ua+vb,

whereu, v∈Z. As in Algorithm 2, we always haveg= 1in line 13 by choice ofs₁, . . . , s_K. 2.4.3 Runtime and Sampling Bounds

In order to prove runtime and sample bounds for Algorithm 3, we need some estimates involving the required primes s₁, . . . , s_K. This can be done using two equivalent formu-lations of thePrime Number Theorem.

Theorem 2.27 (Prime Number Theorem) For x ∈ R define the prime-counting function

π(x) := X

p≤x pprime

1.

Then the following estimates hold.

(i) π(x) = x logx+O

x log²x

, (ii) p_l =llogl+O(llog logl).

See [MV07], Chapter 6.2, Theorem 6.9 for a proof of (i) and [HW60], Chapter 1, Theorem 8 for a proof of (ii). Recall that in line 5 of Algorithm 3 we have to compute CDFTs of lengths_k for allk∈ {1, . . . , K}. Each of them needss_k equidistant samples of f and requiresO(s_klogs_k)arithmetical operations, as already discussed in Section 1.1.2.

Thus, we have to find estimates for

K

X

k=1

s_k and

K

X

k=1

s_klogs_k,

Algorithm 3 Algorithm for Functions with Short Frequency Support II

Input: f ∈ C2π, B, N, where f has bandwidth N and a short frequency support of length at most B < N, and noise threshold ε >0.

Output: The setRof at mostB energetic frequencies off and the vectorxof estimates for their Fourier coefficients.

1: Sets₁ as the smallest prime withs₁ >2B.

2: LetK =blog_s

1Nc+ 1ands₂<· · ·< s_K theK−1smallest primes greater thans₁.

3: for kfrom1 to K do

4: a^sk ←

f

2πj sk

sk−1 j=0 5: ac^sk ←CDFT[a^sk]

6: end for

Identification of the Smallest Energetic Frequencyω1 7: for kfrom 1 to K do

8: νsk ←first support index ofac^sk

9: end for

Reconstruction of ω₁ from its Residues

10: Setk= 1,ω1=νs1 andn=s1.

11: while k < K do

12: Setk=k+ 1

13: (g, u, v)←extended_gcd(n, sk)

14: ω1 =n(((νs_k −ω1)·u) mod s_k) +ω1 15: n=s_k·n

16: end while

17: Setω1 =ω1 mod nand shiftω1 into the range

−_N

2

+ 1, . . . ,_N

2 , sinceN ≤n.

Identification of the Remaining Frequencies and Coefficients

18: for ω fromω₁ to ω₁+B−1 do

19: if

1 K

PK

k=1ac^sk_{ω mod sk}

> εthen

20: R←R∪ {ω}

21: xω← _K¹ PK

k=1ac^sk_{ω mod sk}

22: end if

23: end for Output: R,x.

2.4 Sparse FFT for Functions with Short Frequency Support II

in order to obtain bounds for the number of required samples of the function we aim to recover and the runtime of the computation of the CDFTs.

Lemma 2.28 LetB, N ∈NwithB < N and lets₁ be the smallest prime greater than 2B. SetK =blog_s1Nc+ 1and s2, . . . , sK as theK−1 smallest primes greater thans1. Then the following estimates are satisfied.

(i)

Consequently,pq−1 is the largest prime that is smaller than 2B. Theorem 2.27 (i) yields q−1 =π(2B) = 2B

Note that sK = pq+K−1. Further, by Betrand’s postulate, see [HW60], Chapter 22.3, Theorem 417, there exists at least one prime number between 2B and 2·2B. Thus, we know that 2B < s1<4B, and s1 =O(B). With K=blog_s1Nc+ 1, we obtain

Hence, the second formulation of the Prime Number Theorem (Theorem 2.27 (ii)) implies sK =pq+K−1

Using (2.14), we find for the number of required samples that

K

(ii) In order to deal with the estimates for the runtime of the CDFT computation, we first recall a property of the logarithm. For any a >0, the logarithm satisfies that

log(a·log(a)) = loga+ log log(a) =O(loga). (2.15)

Consequently, combining (2.14) and (2.15) yields

Employing the estimates proven in Lemma 2.28, we can now show the main result about the runtime and sampling complexities of Algorithm 3.

Theorem 2.29 Let B, N ∈ N with B < N and ω1 ∈

−_N

2

+ 1, . . . ,_N

2 . Let f ∈C_2π have a short frequency support of length at mostB with bandwidth N, i.e.,

f(x) =

B−1

X

k=0

cω1+k·e^i(ω1+k)x.

Then Algorithm 3 returns the energetic frequencies of f and their corresponding Fourier coefficients in

time, and has a sampling complexity of O

Proof. By construction of Algorithm 3, we know that it returns the correct frequencies and Fourier coefficients, apart from numerical errors, for exact data. We will now examine the runtime of the different parts of the algorithm using the considerations made above.

Using precomputed lists of the first, e.g., 10,000 primes, the computational costs of finding theK =blog_s1Nc+ 1smallest primes greater than2B in lines 1 and 2 are

By Section 1.1.1 and Remark 1.9, the computation of a CDFT of length M has a runtime ofO(MlogM). Consequently, it follows from Lemma 2.28 (ii) that the CDFTs in lines 3 to 6 require

O

The runtime of line 8 depends on whether we find the first support indicesν_sk of the vector ac^sk via local energies or via block search.

2.4 Sparse FFT for Functions with Short Frequency Support II

(i) If we use the local energies method to detect the first support index, computing the first local energy

we have to execute O(B+sk) additions in order to calculate all sk local energies for a fixedk. Thus, employing this method, lines 7 to 9 have a runtime of

O

(ii) If we want to identify the first support indexνs_kofac^skby looking forB+1consecutive entries that have an absolute value less than the noise threshold ε >0, we need to check each entry at most twice to also allow for blocks that are wrapped periodically around the boundary of the vector. Hence, finding the first support index ofac^sk requiresO(s_k) arithmetical operations, and lines 7 to 9 need

O

Consequently, the theoretical runtime of lines 7 to 9 is always dominated by the compu-tational effort of the CDFT computations in lines 3 to 6.

As we discussed in Section 2.2.1, the runtime of the extended Euclidean algorithm extended_gcd(n, s_k) in line 13 is O(logs_k) if k ≥ 3, since then n = Qk−1

l=1 s_l > s_k. If k = 2, we can still estimate its runtime with O(logs₂), as s₂ > n =s₁. Consequently, the CRT reconstruction procedure in lines 10 to 17, i.e., Algorithm 1, has a runtime of

O

which is again dominated by the runtime of the CDFT computations. Finally, the cal-culation of the coefficient estimates for theB possibly energetic frequencies contained in {ω₁, . . . , ω₁+B−1}in lines 18 to 23 needsO(K·B) =O(K·s_K)arithmetical operations.

Hence, we obtain that the runtime of Algorithm 3 is dominated by the runtime of the CDFT computations in lines 3 to 6, so, by (2.16), we obtain an overall runtime of

O

Further, Lemma 2.28 (i) yields that Algorithm 3 has a sampling complexity of

Im Dokument Sparse Fast Trigonometric Transforms (Seite 67-72)