Simple and Practical Algorithm for the Sparse
Fourier Transform
Haitham Hassanieh

Piotr Indyk

Dina Katabi

Eric Price

MIT

2012-01-19

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

1 / 19

Outline

1

Introduction

2

Algorithm

3

Experiments

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

2 / 19

The Dicrete Fourier Transform

Discrete Fourier transform: given x ∈ Cn , find
X
xbi =
xj ω ij
Fundamental tool
I
I
I
I

Compression (audio, image, video)
Signal processing
Data analysis
...

FFT: O(n log n) time.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

4 / 19

Sparse Fourier Transform

Often the Fourier transform is dominated by a small number of
“peaks”
I

Precisely the reason to use for compression.

If most of mass in k locations, can we compute FFT faster?

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

5 / 19

Previous work

Boolean cube: [KM92], [GL89]. What about C?
[Mansour-92]: k c logc n.
Long list of other work [GGIMS02, AGS03, Iwen10, Aka10]
Fastest is [Gilbert-Muthukrishnan-Strauss-05]: k log4 n.
I
I
I
I

All have poor constants, many logs.
Need n/k > 40,000 or ω(log3 n) to beat FFTW.
Our goal: beat FFTW for smaller n/k in theory and practice.
Result: n/k > 2, 000 or ω(log n) to beat FFTW.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

6 / 19

Our result

Simple, practical algorithm with good constants.
√
Compute the k-sparse Fourier transform in O( kn log3/2 n) time.
Get xb0 with approximation error
kxb0 − xbk2∞ ≤

1
kxb − xbk k22
k

If xb is sparse, recover it exactly.
Caveats:
I
I

Additional kxk2 /nΘ(1) error.
n must be a power of 2.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

7 / 19

Structure of this section

If xb is k-sparse with known support S, find xbS exactly in
O(k log2 n) time.
√
e nk ) time.
In general, estimate xb approximately in O(

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform

2012-01-19

9 / 19

Intuition
Original signal (time)

Original signal (freq)

n-dimensional DFT

Cutoff signal (time)

Cutoff signal (freq)

n-dimensional DFT
of first B terms.

Cutoff signal (time)

Cutoff signal, subsampled (freq)

B-dimensional DFT
of first B terms.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

10 / 19

Framework
Cutoff signal (time)

Cutoff subsampled signals (freq)

“Hashes” into B buckets in B log B time.
Issues:
I

“Hashing” needs a random hash function
F

I

Collisions
F

I
I

Access xt0 = ω −at xσt , so xb0 t = xbσ−1 t+a [GMS-05]
Have B > 4k , repeat O(log n) times and take median. [Count-Sketch,
CCF02]

Leakage
Finding the support. [Porat-Strauss-12], talk at 9:45am.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

11 / 19

Leakage
Cutoff signal (time)

Cutoff, subsampled signals (freq)



1 i <B
be the “boxcar” filter. (Used in
0 otherwise
[GGIMS02,GMS05])
Observe
Let Fi =

b ∗xb, B).
DFT(F ·x, B) = subsample(DFT(F ·x, n), B) = subsample(F
b of boxcar filter is sinc, decays as 1/i.
DFT F
Need a better filter F !
Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

12 / 19

Filters
Filter (time)

2.0

25

Filter (freq)

20

1.5

15
1.0
10
0.5

0.00

5

20

40

60

80

100

0

Bin

b ∗ xb, B) in O(B log B) time.
Observe subsample(F
Needs for F :
I
I
I

supp(F ) ∈ [0, B]
b | < δ = 1/nΘ(1) except “near” 0.
|F
b ≈ 1 over [−n/2B, n/2B].
F

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

13 / 19

Filters
Filter (time)

1.0

18

Filter (freq)

16
0.8

14
12

0.6

10
8

0.4

6
4

0.2

2
0.00

20

40

60

80

100

0

Bin

b ∗ xb, B) in O(B log B) time.
Observe subsample(F
Needs for F :
I
I
I

supp(F ) ∈ [0, B]
b | < δ = 1/nΘ(1) except “near” 0.
|F
b ≈ 1 over [−n/2B, n/2B].
F

Gaussians:
I
I
I
I

p
Standard deviation σ
p= B/ log n
DFT has σ
b = (n/B) log n
Nontrivial leakage into O(log n) buckets.
But likely trivial contribution to correct bucket.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

13 / 19

Filters
Filter (time)

1.0

80

Filter (freq)

70
0.8

60
50

0.6

40
0.4

30
20

0.2

10
0.00

50

100

150

200

0

Bin

b ∗ xb, B) in O(B log B) time.
Observe subsample(F
Needs for F :
I
I
I

supp(F ) ∈ [0, B log n]
b | < δ = 1/nΘ(1) except “near” 0.
|F
b
F ≈ 1 over [−n/2B, n/2B].

Gaussians:
I
I
I
I

p
p 

B/
Standard deviation σ = 
log n B · log n
p
p 

DFT has σ
b =
(n/B)
 log n (n/B)/ log n
Nontrivial leakage into 0 buckets.
But likely trivial contribution to correct bucket.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

13 / 19

Filters
Filter (time)

1.0

80

Filter (freq)

70
0.8

60
50

0.6

40
0.4

30
20

0.2

10
0.00

50

100

150

200

0

Bin

p
Let G be Gaussian with σ = B log n
H be box-car filter of length n/B.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

14 / 19

Filters
Filter (time)

0.10

1.2

0.08

1.0

0.06

0.8

0.04

0.6

0.02

0.4

0.00

0.2

0.020

0.0

50

100

150

200

Filter (freq)

Bin

p
Let G be Gaussian with σ = B log n
H be box-car filter of length n/B.
b ∗ H.
b =G
Use F
I
I
I

b so supp(F ) ⊂ [0, B log n].
F = G · H,
b
|F | < 1/nΘ(1) outside −n/B, n/B.
b | = 1 ± 1/nΘ(1) within n/2B, n/B.
|F

Hashes correctly to one bucket, leaks to at most 1 bucket.
Replace Gaussians with “Dolph-Chebyshev window functions”:
factor 2 improvement.
Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

14 / 19

Algorithm to estimate xbS

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Algorithm to estimate xbS

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Algorithm to estimate xbS

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Algorithm to estimate xbS

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Algorithm in general

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Algorithm in general

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.
To find S: choose all that map to the top 2k values.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Algorithm in general

For O(log n) different permutations of xb, compute
b ∗ xb, B).
subsample(F
Estimate each xi as median of values it maps to.
To find S: choose all that map to the top 2k values.
nk /B candidates to update at each iteration: total
(

√
nk
+ B log n) log n = nk log3/2 n
B

time.
Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

15 / 19

Empirical Performance: runtime
Run Time vs Signal Size (k=50)
10

22

Run Time vs Sparsity (n=2 )

sFFT 1.0
sFFT 2.0

10

FFTW
FFTW OPT

1

Run Time (sec)

Run Time (sec)

AAFFT 0.9

0.1

0.01

1

sFFT 1.0

0.1

sFFT 2.0
FFTW

0.001

FFTW OPT
AAFFT 0.9

0.01

2

14

2

15

2

16

2

17

2

18

2

19

2

20

2

21

2

22

2

23

2

24

2

25

2

26

100

Signal Size (n)

1000

Number of Non-Zero Frequencies (k)

Compare to FFTW, previous best sublinear algorithm (AAFFT).
e 1/3 k 2/3 ).
Offer a heuristic that improves time to O(n
I
I

Filter from [Mansour ’92].
Can’t rerandomize, might miss elements.

Faster than FFTW for n/k > 2,000.
Faster than AAFFT for n/k < 1,000,000.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

17 / 19

Empirical Performance: noise

22

Robustness vs SNR (n=2 , k=50)
sFFT 1.0

Average L1 Error per Enrty

1

sFFT 2.0
AAFFT 0.9

0.1
0.01
0.001
0.0001
1e-05
1e-06
1e-07
-20

0

20

40

60

80

100

SNR (dB)

Just like in Count-Sketch, algorithm is noise tolerant.

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

18 / 19

Conclusions

Roughly: fastest algorithm for n/k ∈ [2 × 103 , 106 ].
Recent improvements [HIKP12b?]
I
I
I

O(k log n) for exactly sparse xb
O(k log kn log n) for approximation.
Beats FFTW for n/k > 400 (in the exact case).

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

19 / 19

Hassanieh, Indyk, Katabi, and Price (MIT) Simple and Practical Algorithm for the Sparse Fourier Transform 2012-01-19

20 / 19