Source-linked AI summary
Nearly Optimal Sparse Fourier Transform
Haitham Hassanieh, Piotr Indyk, Dina Katabi, Eric Price
TL;DR
The paper asks whether sparse Fourier transforms can be computed faster than the FFT across the full sparse regime. It develops randomized sublinear-time algorithms using binning, coefficient isolation, and blockwise index recovery, and establishes a sampling lower bound. The algorithms achieve O(k log n) time for exactly k-sparse signals and O(k log n log(n/k)) time for general signals, while general algorithms require Ω(k log(n/k)/ log log n) samples.
Problem
Existing sparse Fourier algorithms did not improve over FFT for all k = o(n), and many were complex or had large big-Oh constants.
Method
The paper bins Fourier coefficients with concentrated filters, locates and estimates isolated coefficients, iteratively updates the signal, and recovers index bits in blocks for the general case.
Results
O(k log n) time is achieved for exactly k-sparse inputs and O(k log n log(n/k)) time for general inputs, with o(n log n) time for every k = o(n).
Takeaways & Limitations
These are the first known algorithms improving over FFT throughout k = o(n); the exactly sparse algorithm is optimal under an FFT-optimality assumption for k = n^Ω(1).
Takeaways & Limitations
The lower bound applies to adaptive measurements induced by the Fourier matrix, not arbitrary adaptive linear measurements.
Abstract
from arXiv · showhide
We consider the problem of computing the k-sparse approximation to the discrete Fourier transform of an n-dimensional signal. We show: * An O(k log n)-time randomized algorithm for the case where the input signal has at most k non-zero Fourier coefficients, and * An O(k log n log(n/k))-time randomized algorithm for general input signals. Both algorithms achieve o(n log n) time, and thus improve over the Fast Fourier Transform, for any k = o(n). They are the first known algorithms that satisfy this property. Also, if one assumes that the Fast Fourier Transform is optimal, the algorithm for the exactly k-sparse case is optimal for any k = n^{Ω(1)}. We complement our algorithmic results by showing that any algorithm for computing the sparse Fourier transform of a general signal must use at least Ω(k log(n/k)/ log log n) signal samples, even if it is allowed to perform adaptive sampling.
1 Introduction
Sparse Fourier algorithms target sublinear computation when Fourier outputs are sparse, addressing prior runtime and complexity limitations. This paper introduces faster randomized algorithms, supporting techniques for binning and recovery, and a sample lower bound for general signals.
- Motivation: O(n log n) is the FFT runtime for an n-dimensional DFT, while sparse outputs motivate sublinear algorithms in applications including compression and signal processing.Sparse Fourier methods exploit signals with only a small number of significant frequency coefficients.
- Prior limitations: Prior algorithms did not improve over FFT for the whole range k = o(n) and often had high complexity or large big-Oh constants.The cited prior methods improved over FFT only in narrower sparsity regimes.
- Results: O(k log n) time is achieved for exactly k-sparse inputs, while the general case takes O(k log n log(n/k)) time.Both algorithms are randomized and require n to be a power of 2.
- Results: o(n log n) time improves over FFT for every k = o(n), making these the first known algorithms with this property.Under an FFT-optimality assumption, the exactly sparse algorithm is optimal for k = n^Ω(1).
- Empirical result: A preliminary implementation beat FFTW at n = 2^22 and k ≤ 2^17, whereas a prior algorithm did so only for k ≤ 2000.The exactly sparse algorithm is described as simple and having low big-Oh constants.
- Lower bound: General-case algorithms require at least Ω(k log(n/k)/ log log n) samples, even with adaptive sampling limited to Fourier-induced measurements.The lower bound is weaker than the Ω(k log(n/k)) bound for arbitrary linear measurements.
2 Preliminaries
The paper establishes notation and assumptions for sparse Fourier recovery, including spectrum operations, filtering windows, and power-of-two signal dimensions.
- Notation: The Fourier spectrum is denoted by bx, while supp(x) and ∥x∥0 represent coordinate support and its number of nonzero entries.Indices are interpreted modulo n, enabling the paper’s convolution and coordinate-wise product notation.
- Definitions: Spectrum permutation and flat filtering windows are the two prior tools used throughout the paper.The windows have parameters B, δ, and α, with efficient constructions and evaluation procedures described later.
- Definitions: A spectrum permutation Pσ,a,b is defined when σ−1 exists modulo n.The supplied passage introduces the permutation notation but does not include its full displayed formula.
- Preliminaries: A prior lemma bounds the probability that a randomly chosen odd multiplier maps a nonzero index into a short modular interval by 4C/n.The lemma assumes n is a power of two and σ is uniformly random among odd numbers in [n].
- Assumptions: The paper assumes that n, the dimension of all vectors, is an integer power of 2.For exactly sparse inputs, Fourier coefficients are also assumed bounded by a precision parameter L=n^O(1); the general case has an additional approximation-norm condition.
3 Algorithm for the exactly sparse case
The exactly sparse algorithm repeatedly hashes Fourier coefficients into bins, identifies isolated coefficients, and updates the residual until exact recovery is achieved.
- Hashing and filtering: HASHTOBINS permutes the residual spectrum and hashes it into B bins using a flat window.When B divides n, the binned values can be computed with a B-dimensional FFT in O(B log B) time.
- Hashing and isolation: Each iteration maps frequencies into bins using a random spectrum permutation, where collisions and large offsets are the two failure events.The bin is hσ,b(i), and the offset is oσ,b(i); large offsets satisfy |oσ,b(i)| ≥ (1−α)n/(2B).
- Failure probabilities: A coefficient collides with probability at most 4|S|/B and experiences a large offset with probability at most α.Here S is the support of the current residual, and the offset bound follows from random shifting.
- Iterative recovery: The algorithm uses O(log k) iterations while progressively reducing the residual support, and the expected total runtime sums to O(k log n).The analysis bounds the probability of excessive unsuccessful iterations and derives the per-round runtime before summing across rounds.
- Main guarantee: Theorem 3.8 gives expected O(k log n) runtime and at least 2/3 success probability for recovering a k-sparse Fourier vector with bounded coefficients.The result assumes known L=n^O(1) and coefficients in {−L, …, L}.
4 Algorithm for the general case
SPARSEFFT handles general signals by iteratively locating and estimating heavy Fourier coefficients, progressively reducing the residual sparsity. Its location procedure uses blockwise phase information and non-adaptive measurements, yielding the stated recovery guarantee and runtime.
- Overview: SPARSEFFT first locates most heavy coordinates, estimates their values, and repeats on the residual with progressively smaller sparsity.Each iteration aims to make the residual approximately k/4-sparse, so later iterations run faster.
- Location: The noisy location procedure partitions frequencies into O(k/ε) bins and uses phase differences to narrow each candidate region.A t-ary search splits regions into t = Θ(log n) subregions and checks representative candidates by phase agreement.
- Location: O(log(n/k)) measurements and O(log(n/k) log n) decoding time suffice to locate frequencies, with measurements performed non-adaptively across bins.The resulting failure probability per well-hashed coordinate is 1/t^Θ(1) = o(1).
- Estimation: Median estimation makes at least 7k/8 heavy coordinates both located and well estimated with 7/8 probability, leaving a residual close to k/4-sparse.The largest k estimated coordinates contain good estimates for 3k/4 large coordinates with probability 3/4.
- Guarantee: The full algorithm succeeds with probability 2/3 and achieves the paper’s stated general-case running-time bound.The supplied runtime expression is distributed across the cited theorem and concluding runtime passages.
5 Reducing the full k-dimensional DFT to the exact k-sparse case in n dimensions
The paper reduces a full k-dimensional DFT to an exactly k-sparse DFT in n dimensions by periodically extending the input, preserving recoverability while incurring only the sparse algorithm’s runtime.
- Reduction: Assuming k divides n, periodically extending x by yi = x_i mod k produces an n-dimensional signal whose Fourier transform is k-sparse.The nonzero Fourier coefficients occur at multiples of n/k.
- Reduction: The original transform can be recovered immediately from the transformed periodic signal, so an algorithm running in T(k) yields an O(T(k)) algorithm for the k-dimensional DFT.Each requested sample of y is computed from x in constant time.
- Lower bound consequence: If the n-dimensional DFT requires Ω(n log n) time, every arbitrary-dimension k-sparse DFT algorithm requires Ω(k log k) time.This corollary follows by applying the reduction to an exactly sparse algorithm.
6 Lower Bound
The lower-bound argument shows that general sparse Fourier recovery requires many time-domain samples even under adaptive sampling. It uses a hard distribution and mutual-information bounds to establish the sample requirement.
- Statement: Any algorithm satisfying the paper’s sparse-recovery guarantee must access Ω(k log(n/k)/log log n) samples.The lower bound applies to algorithms using adaptive sampling.
- Hard distribution: The proof constructs a family of well-separated k-sparse supports and random signals supported on those sets, then adds Gaussian noise.The resulting observations form a Markov chain from the support through measurements to the recovered support.
- Information lower bound: Successful recovery implies that the recovered support identifies the planted support with sufficiently small error, enabling a mutual-information lower bound.The argument uses Fano’s inequality after relating recovery error to support error.
- Adaptive sampling: An upper bound on information per adaptive Fourier sample is combined with the lower bound I(S;S′) = Ω(k log(n/k)).The analysis conditions on previous observations and tracks the adaptively chosen measurement vectors.
- Conclusion: The theorem and corollary conclude that adaptive Fourier-coordinate sampling cannot avoid the Ω(k log(n/k)/log log n) sample requirement.For non-adaptive samples, the cited prior bound would be Ω(k log(n/k)); the paper establishes its weaker bound under adaptivity.
7 Efficient Constructions of Window Functions
The paper constructs flat window functions by combining Gaussian and box-car filters, approximating Gaussian cumulative-distribution evaluations efficiently. The resulting windows have bounded values and efficient construction and evaluation costs.
- Approximation: Gaussian CDF values can be approximated to precision ±δ in O(log(1/δ)) time.The approximation uses constant-tail behavior and an efficient expansion near zero.
- Complexity: The window G can be computed over its support in O(B log(n/δ)/α) time, while each cG′_i evaluation costs O(log(n/δ)) time.The evaluation is O(1) outside the specified transition range.
- Construction: A Gaussian convolved with a box-car filter provides the continuous prototype for the required flat window functions.The Fourier transform of the prototype is expressed as a difference of Gaussian CDF values.
- Discretization: For suitable parameters, discretization yields window values within the stated approximation error of the continuous prototype.The discretized values satisfy cG′_i = bH_i ± 6δ.
- Properties: The constructed windows satisfy cG′_i ∈ [0,1] for every i and handle the B = 1 case with a constant function.For B ≥ 2, the parameter conditions ensure the required flatness properties.
8 Open questions
The paper identifies open directions involving runtime, sample complexity, generalization, input-size assumptions, and failure probability. These include designing faster general-signal algorithms, reducing samples, extending tasks and dimensions, and avoiding amplification slowdowns.
- An O(k log n)-time algorithm for general signals remains an open design goal, or requires proving that no such algorithm exists under reasonable assumptions.
- Reducing sample complexity remains open because current sample counts are bounded only by the algorithms’ running times.
- Extending the results to related tasks, including the sparse Walsh-Hadamard Transform, is an unresolved direction.
- Extending the algorithm beyond n that is a power of 2 remains open, although some earlier algorithms work for arbitrary n.
- The algorithms currently succeed only with constant probability, and straightforward amplification incurs a log(1/p) slowdown to reach probability 1 − p.