Tighter Fourier Transform Complexity Tradeoffs

The Fourier Transform is one of the most important linear transformations used in science and engineering. Cooley and Tukey's Fast Fourier Transform (FFT) from 1964 is a method for computing this transformation in time $O(n\log n)$. Achieving a matching lower bound in a reasonable computational model is one of the most important open problems in theoretical computer science. In 2014, improving on his previous work, Ailon showed that if an algorithm speeds up the FFT by a factor of $b=b(n)\geq 1$, then it must rely on computing, as an intermediate"bottleneck"step, a linear mapping of the input with condition number $\Omega(b(n))$. Our main result shows that a factor $b$ speedup implies existence of not just one but $\Omega(n)$ $b$-ill conditioned bottlenecks occurring at $\Omega(n)$ different steps, each causing information from independent (orthogonal) components of the input to either overflow or underflow. This provides further evidence that beating FFT is hard. Our result also gives the first quantitative tradeoff between computation speed and information loss in Fourier computation on fixed word size architectures. The main technical result is an entropy analysis of the Fourier transform under transformations of low trace, which is interesting in its own right.


Introduction
The (discrete) normalized Fourier transform is a complex linear mapping sending an input x ∈ C n to y = F x ∈ C n , where F is an n × n unitary matrix defined by F (k, ℓ) = n −1/2 e −i2πkℓ/n . (1.1) The Fast Fourier Transform (FFT) of Cooley and Tukey [7] is a method for computing the Fourier transform of a vector x ∈ C n in time O(n log n) using a so called linear algorithm. A linear algorithm, as defined in [12], is a sequence F 0 , F 1 , . . . , where each F i is a set of affine functions, for each i ≥ 0 F i+1 = F i ∪ {λ i f + µ i g} for some λ i , µ i ∈ C and f, g ∈ F i , and F 0 contains (projections onto) the input variables as well as constants.
As for lower bounds, it is trivial that computing the Fourier Transform requires a linear number of steps. Papadimitriou, for example, derives in [13] an Ω(n log n) lower bound for Fourier transforms in finite fields using a notion of an information flow network. It is not clear how to extend that result to the Complex field. There have also been attempts [15] to reduce the constants hiding in the upper bound of O(n log n), while also separately counting the number of additions versus the number of multiplications (by constants). In 1973, Morgenstern proved that if the modulus of the λ i 's and µ i 's is bounded by 1 then the number of steps required for computing the unnormalized Fourier transform, defined by n 1/2 F in the linear algorithm model is at least 1 2 n log 2 n. He used a potential function related to matrix determinant, which makes the technique inapplicable for deriving lower bounds for the (normalized) F . Morgenstern's result also happens to imply that the transformation √ n Id ( √ n times the identity) has the same complexity as the Fourier transform, which is not a satisfying conclusion. We also note that the unnormalized Fourier transform, as studied by Morgenstern cannot be simply implemented on a computer of fixed word length, because that transformation increases norms by √ n, and hence Ω(log n) bit words are necessary.
Ailon [1] studied the complexity of the (normalized) Fourier transform in a computational model allowing only orthogonal transformations acting on (and replacing in memory) two intermediates at each step. He showed that at least Ω(n log n) steps were required. The proof was done by defining a potential function on the matrices M (t) defined by composing the first t gates. The potential function is simply the sum of Shannon entropy of the probability distributions defined by the squared modulus of elements in the matrix rows. (Due to orthogonality, each row, in fact, thus defines a probability distribution). That result had two shortcomings: (i) The algorithm was assumed not to be allowed to use extra memory in addition to the space used to hold the input. In other words, the computation was done in place. (ii) The result was sensitive to the normalization of F , and was not useful in deriving any lower bound for γF for γ ∈ {0, ±1}.
In [2], Ailon took another step forward by showing a lower bound for computing any scaling of the Fourier transform in a stronger model of computation which we call uniformly well conditioned. At each step, the algorithm can perform a nonsingular linear transformation on at most two intermediates, as long as the matrix M (t) defining the composition of the first t steps must have condition number at most κ, for all i. We remind the reader that condition number is defined as M (t) · (M (t) ) −1 , where · is spectral norm. Equivalently, it is defined as σ 1 (M (t) )/σ n (M (t) ), where σ 1 (A) ≥ · · · ≥ σ n (A) are singular values of a matrix A. Otherwise stated, the result implies that if an algorithm computes the Fourier transform in time (n log n)/b for some b > 1, then some M (t) must have condition number at least Ω(b). This means that the computation output relies on an ill conditioned intermediate step. It is not difficult to see from the proof of the main result of [2] that, in fact most of the algorithm steps have condition number at least Ω(b). Hence, the bottleneck is ubiquitous. (We will return to the point of bottleneck ubiquity shortly.) Ill conditionedness, as claimed in [2], leads to numerical instability and compromise in robustness against noise. In this work we make a quantification of these claims. Roughly speaking, we show that a speed up of factor b implies that we either have to represent scalars with Ω(log b) bits (implying increased computation complexity in bit operations), or we'd have to sacrifice numerical errors Ω(b) times larger than those of standard FFT, for a given fixed word size architecture.

Main Contribution
Our first result (Theorem 5.1) is that, a speedup by a factor of b requires not only b ill-conditionedness in the usual sense, but requires a stronger notion of ill conditionedness which we define below. The theorem generalizes the main result in [2], while also providing a simpler proof.
Our second and main result (Theorem 5.3) states that, roughly speaking, in an architecture of fixed word size, speeding up FFT by a factor of b implies the existence of Ω(n) orthogonal directions in input space that are numerically unstable in that they cause overflow or underflow at Ω(n) different time steps. (The exact definition of overflow and underflow will be defined below). Avoiding these numerical instabilities could be technically achieved by allowing word size of length Ω(log b) at the offending memory word locations at particular time steps, but this incurs a bit operation complexity of Ω(n log b) because the number of such steps is Ω(n).

Fourier Transform Types
In theory and practice of engineering and computer science, there are two main types of Fourier transforms considered. The n dimensional Discrete Fourier Transform (DFT), as defined in (1.1), is a complex unitary mapping defined by the characters of the Abelian group Z/nZ. The Walsh-Hadamard transform is a real orthogonal mapping defined by the characters of the n dimensional binary hypercube (Z/2Z) n . More precisely, for n an integer power of 2, the (i, j)'th matrix element is defined as , where ·, · is dot-product, and [p] denotes (here only) the bit representation of the integer p ∈ {0, . . . , n − 1} as a vector of log n bits. Similarly to FFT for DFT, there is an O(n log n) time algorithm for computing the Walsh-Hadamard transform of an input x.
In this work we consider computation over the reals, and must hence view the n-DFT as an orthogonal mapping in R 2n in a standard way. (Each complex input coordinate is viewed as two real coordinates, one representing the real part and the other the imaginary part. Each complex matrix element of DFT is, accordingly, viewed as a two-by-two real matrix.) Throughout, we will assume n is an integer power of 2 and will use F to denote either the n−Walsh-Hadamard transform or the real representation of the n/2-DFT transform. In order for our results to work for the DFT case, however, we will have to make a slight modification to the potential function Φ defined in [2]. (This was not an issue there). We will present the modified function Φ C below and briefly explain where it is necessary, but will omit the technical details of working with Φ C , which are not much different than for Φ. For simplicity, the reader is invited to assume that F denotes the Walsh-Hadamard transform, and FFT refers to the fast Walsh-Hadamard transform.

Paper Structure
We start by defining the computational model in Section 3. In Section 4 we introduce further notation and recall (and generalizee) the matrix quasi-entropy function which was first defined in [2]. In Section 5 we present our results. The proofs, together with useful lemmas, are deferred to Section 6. One particular Lemma 6.4, a key to proving Theorem 5.3, is interesting in its own right.

Computational Model and Notation
We remind the reader of the computational model discussed in [1,2], which is a special case of the linear computational model. The machine state represents a vector in R ℓ for some ℓ ≥ n, where it initially equals the input x ∈ R n (with possible padding by zeroes, in case ℓ > n). Each step (gate) is either a rotation or a constant. A rotation applies a 2-by-2 rotation mapping on a pair of machine state coordinates (rewriting the result of the mapping to the two coordinates). We remind the reader that a 2-by-2 rotation mapping is written in matrix form as cos θ sin θ − sin θ cos θ for some real (angle) θ. A constant gate multiplies a single machine state coordinate (rewriting the result) by a nonzero constant. In case the constant equals −1, we call it a reflection gate. In case ℓ = n we say that we are in the in-place model. Any nonsingular linear mapping over R n can be decomposed into a sequence of rotation and constant gates in the in-place model, and hence our model is, in a sense, universal. FFT works in the in-place model, using rotations (and possibly reflections) only. A restricted method for dealing with ℓ > n was developed in [2], and can be applied here too in a certain sense (see Section 7 for a discussion). We focus in this work on the in-place model only.
Since both rotations and constants apply a linear transformation on the machine state, their composition is a linear transformation. If A n is an in-place algorithm for computing a linear mapping over R n , it is convenient to write it as A n = (M (0) = Id, M (1) , . . . , M (m) ) where m is the number of steps (gates), M (t) ∈ R n×n is the mapping that satisfies that for input x ∈ R n (the initial machine state), M (t) x is the machine state after t steps. (Id is the identity matrix). The matrix M (m) is the target transformation, which will typically be F in our setting. For any t ∈ [m], if the t'th gate is a rotation, then M (t) defers from M (t−1) in at most two rows, and if the t'th gate is a constant, then M (t) defers from M (t−1) in at most one row.
Both the n-Walsh-Hadamard and the real representation of the (n/2)-DFT can be implemented in-place using O(n log n) rotations.

Numerical Architecture
The in-place model implicitly assumes representation of a vector in R n in memory using n words. A typical computer word represents a coordinate (with respect to some fixed orthogonal basis) in the range [−1, 1] to within some accuracy ε. The accuracy level ε is Θ(1). For sake of simplicity, it should be thought of as 2 −31 or 2 −63 in modern computers of 32 or 64 bit words, respectively.
To explain the difficulties in speeding up FFT on computers of fixed precision words in the in-place model, we need to understand whether (and in what sense) standard FFT is at all suitable on such machines. First, we must restrict the domain of inputs. Clearly this domain cannot be R n , because computer words can only represent coordinates in the range [−1, 1], by our convention. We could restrict the input to come from the cube [−1, 1] n . Initializing the machine with such input is easy, but the output (as well as intermediate steps) of the standard FFT map the cube to orthogonal transformations thereof, which might exceed the cube. The extreme case is an orthogonal transformation that translates the unit vector (1/ √ n, . . . , 1/ √ n) to (0, . . . , 1). Under such transformation, an input of (1, . . . , 1) from the cube is mapped to (0, . . . , 0, √ n), overflowing in the last coordinate. It makes more sense to consider input from an n-ball of radius Θ( √ n), which we denote ). An n-ball is invariant to orthogonal transformations. Additionally, a coordinate of a random input from B(Θ( √ n)) is of typical magnitude Θ(1), and hence typically does not overflow. The problem is that, using well known principles from high dimensional geometry, it is very likely for some coordinate to be of magnitude Θ( √ log n). Hence, paradoxically, it is likely that some coordinate would overflow. Perhaps the reason that this issue is hardly every considered in literature on Fourier transform complexity (including the author's prior work) is because a number of magnitude Θ( √ log n) can be accommodated using Θ(log log n) bit words. It is unlikely to encounter an application with n = Ω 2 2 32 (let alone Ω 2 2 64 ), and hence this is really a non-issue for applications. Nevertheless, it's an issue for complexity theoretical purposes. So is the true running time of standard FFT O(n log n) × (time to multiply two Θ(log log n) bit words) , owing to a tiny minority of coordinates that overflow? The answer is no, because on expectation the number of bits required to encode a coordinate x(i) of a random vector x from B(Θ( √ n)) to within accuracy ε = Θ(1), is Θ(1). 1 In fact, using state-of-the-art integer multiplication technology [8,9], the expected time required to multiply such a random coordinate with a number of Θ(1) bits is Θ(1) when taking logical bit operations into account. However, we prefer to avoid dealing with bit operations in this work, and leave it to future work discussion (see Section 7). Instead, we shall take an information theoretical approach, and think of a word in memory, at a particular time step of the algorithm, as a channel that is used for encoding a random coordinate. If the coordinate is drawn from a distribution that requires O(1) bits on expectation to encode to within the fixed accuracy of ε, then no overflow occurs. If ω(1) bits are necessary, then we define this as overflow. 2 By this definition, standard FFT for input drawn uniformly from B(Θ( √ n)) does not overflow at all, because any coordinate of the machine state at any step is tightly concentrated (in absolute values) around Θ(1).
It will be easier to replace the uniform distribution from the ball with the multivariate Gaussian N (0, Id(Θ(n))), which is a good approximation of the former for large n. With this assumption, any coordinate of the standard FFT machine state at any step follows the law N (0, Θ(1)). By simple integration against the Gaussian measure, one can verify that the expected number of bits required to encode such a random variable (to within fixed accuracy ε) is Θ(1), hence no overflow occurs. This input assumption together with the no-overflow guarantee will serve as our benchmark.
Two notes regarding the numerical architecture: 1. Our definition of overflow is counterintuitive, because we are used to thinking about overflow as an offending machine state at a particular step of the algorithm execution for a particular input, while our definition is stochastic. The reason we use this definition is from the combination of (a) our desire to work with a spherically symmetric input and (b) avoiding measuring complexity in the granularity of logical bit operations, stemming from the varying word length typically arising even in the standard FFT benchmark. 3 2. We are completely ignoring the fact that the machine state after t steps on input x in any computer is not only a result of quantizing the vector M (t) x. Rather, errors are accumulated from the effects of quantizing at earlier steps. Taking accumulated errors into account should affect the numerical accuracy of both the benchmark and of any speed-up, but quantifying this effect seems extremely difficult. Again, we take an information theoretical approach by saying that the machine state after t steps contains at most the information in the quantization of 1 By "encoding" here we simply mean the base-2 representation of the integer ⌊x(i)/ε⌋. 2 By encoding of a random number X here we simply mean the integer representation of its quantization ⌊X/ε⌋. 3 It is possible to somewhat practically justify the stochastic definition of overflow by thinking of running FFT on a large number L of iid inputs, which can be thought of as being stacked as columns of an input matrix X. The larger L is, the more concentrated the total number of bits required to encode each row of the matrix around Θ(L) will be. For L polynomial in n, the probability of requiring more than Θ(L) bits per row becomes exponentially small at any step of standard FFT. We could have defined this as our benchmark, but we do not think that anyone would want to compute over a tightly compressed representation of matrix rows into arrays of words. Hence we've decided to settle for the stochastic definition of overflow. the coordinates of M (t) x, and whatever information lost (due to this quantization) cannot be later recovered because the machine state encodes everything that is known about the input (and output) at any given step.
Our chosen approach (and compromises) enable us to obtain what we believe are important insights into inherent tradeoffs of information version computation speed.

The Matrix Quasi-Entropy Function
The set {1, . . . , q} is denoted by [q]. By R a×b we formally denote matrices of a rows and b columns. Matrix transpose is denoted by (·) T . We use (·) −T as shorthand for ( We slightly abuse notation and extend the definition of the quasi-entropy function Φ(M ) defined on nonsingular matrices M from [2], as follows. Given two matrix arguments A, B ∈ R a×b for some This extends naturally to vectors, namely for u, v ∈ R a , Φ(u, v) is as above by viewing R a as R a×1 . If A, B ∈ R a×b and a, b are even, then we define the complex quasi-entropy function Φ C (A, B) to be: The function Φ C can be used for proving our results for the real representation of the complex DFT, which we omit from this manuscript for simplicity. The reason we need this modification to Φ for DFT is explained in the proof of Lemma 6.4, needed by Theorem 5.3 below. Elsewhere, we will work (for convenience and brevity) only with Φ. Abusing notation, and following [2], we define for any nonsingular matrix M : It is easy to see that Φ(F ) = n log n for the Walsh-Hadamard transform, because all matrix elements are ±1/ √ n. If F is a real representation of the (n/2)-DFT, then it is easy to see that Φ C (F ) = n log(n/2), because all matrix elements of the (complex representation of the) (n/2)-DFT are complex unit roots divided by n/2. It will be also useful to consider a generalization of the potential of a nonsingular matrix M , by allowing linear operators acting on the rows of M and M −T , respectively. More precisely, we will let Φ P,Q (M ) be shorthand for where P, Q ∈ R n×a are some mappings. (We will only be working with projection matrices P, Q here). Similarly, Φ C P,Q (M, M −T ) := Φ C (M P, M −T Q). Finally, for any matrix A ∈ R n×n , let σ 1 (A), . . . , σ n (A) denote its singular values, where we use the convention σ 1 (A) ≥ · · · ≥ σ n (A). If A is nonsingular, then the condition number κ(A) is defined by σ 1 (A)/σ n (A). For any matrix A, we let A denote its spectral norm and A F its Frobenius norm. If x is a vector, hence, x = x 2 = x F . Let B denote the Euclidean unit ball in R n 5 Results

FFT Speedup ⇒ Bottleneck
The result for this section should interpreted as follows: If an in-place algorithm A n = (M (0) = Id, . . . , M (m) = F ) speeds up FFT by a factor of b ≥ 1, then there exists an intermediate step M (t) that is ill conditioned (in a generalized sense, to be explained).
Theorem 5.1. Fix n, and let A n = {Id = M (0 , . . . , M (m) } be an in-place algorithm computing some linear function in R n and let P, Q ∈ R n×n be two matrices. For any t ∈ [m], let {i t , j t } denote the set of at most two indices that are affected by the t'th gate (if the t'th gate is a constant gate, then i t = j t , otherwise it's a rotation acting on indices i t , j t ). Then for Additionally, if R = 1 then the t'th gate can be assumed to be a rotation.
In particular, if M (m) = F and m = (n log n)/b for some b ≥ 1 ("A n speeds up FFT by a factor of b") and P = Q = Id, then For the main result in this paper in the next section, we will only need the case R = 1 of the theorem. It is worthwhile, however, to state the following corollary for general R, which is a simple consequence of the fact that for any matrix X and integer q, q i=1 σ 2 q (X) ≥ X(I, :) 2 F whenever |I| ≤ q.
This repeats the main result in [2], for the case R = 1 (our proof is also simper). The bound (5.3) is stronger in the sense that it imposes a stronger restriction on the spectrum of M (t) for larger values of R. The interesting case is, say, R = 2 b/2 . With this choice, (5.3) becomes: One extreme case of (5.4) is when σ 1 (M (t) ) = Ω(2 b/2 ) and σ 2 (M (t) ) = · · · = σ n (M (t) ) = 1. In this case, the condition number is exponentially higher than the bound achieved in [2]. The other extreme is the case σ 1 (M (t) ) = · · · = σ 2 b/2+1 (M (t) ) = α and σ n−2 b/2+1 +1 (M (t) ) = · · · = σ n (M (t) ) = 1, for some global α > 1. In other words, α-ill-conditionedness is exhibited by 2 b/2 disjoint pairs of eigenvalues of M (t) . We do not continue the discussion of this extended notion of ill condition, because we only need the R = 1 case for our main result in what follows.

Implications of Theorem 5.1 on Numerical Accuracy
We discuss the implication of Theorem 5.1, in case R = 1, P = Q = Id. The theorem implies that an algorithm with m = (n log n)/b must exhibit an intermediate matrix M (t) and a pair of indices i t , j t such that the t'th gate is a rotation acting on i t , j t and additionally: Hence, either Case (i). We can assume wlog that Then, letting x T over = M (t) (i t , :)/ M (t) (i t , :) ∈ R n (x over is the normalized i t 'th row of M (t) , transposed), we get that the i t 'th coordinate of M (t) x over is at least √ b. But this implies that, for input x = x over , the t'th step acts on a pair of numbers (M (t) x over )(i t ), (M (t) x over )(j t ), at least one of which requires representation using O(log b) in the numerical architecutre described in Section 3. Now, for input x ∼ N (0, Θ(n)), the i t 'th coordinate just before the application of the t'th gate equals M (t) (i t , :) x T x over , and is hence distributed Using (5.5), this distribution is N (0, Ω(b)). If b = ω(1), then by our definition we reach overflow.
Note also that it is possible as a preprocessing step to replace x with x − (x T x over )x over (eliminating the overflow component), and then to reintroduce the offending component by adding (x T x over )F x over as a postprocessing step. In the next section, however, we shall show that, in fact, there must be Ω(n) pairwise orthonormal directions (in input space) that cause such an overflow at Ω(n) different time steps, so such a simple "hack" cannot not work.
Case (ii). This scenario, as the reader guesses, should be called underflow. In case (ii), wlog Now define x T under = (M (t) ) −T (i t , :)/ (M (t) ) −T (i t , :) ∈ R n , and consider the orthonormal basis u 1 , . . . u n ∈ R n so that u 1 = x under . For any t ′ ∈ [m] (and in particular for t ′ = t): Now notice that the i t 'th coordinate of (x T over (M (t) ) −1 ) has magnitude at least b/2 by (5.6) and the construction of x over . Also notice that for all i = i t , the row M (t) (i, :) is orthogonal to x over , by matrix inverse definition. This means that coordinate i = i t of M (t) x contains no information about g 1 . All the information in g 1 is hence contained in (M (t) x)(i t ). More precisely, g 1 is given by where e is a random variable independent of g 1 . But |(( is known only up to an additive error of ε, due to our assumptions on quantization in the numerical architecture. This means that g 1 can only be known up to an additive error of at least ε b/2, for any value of e. It is important to note that this uncertainty cannot be "recovered" later by the algorithm, because at any step the machine state contains all the information about the input (aside from the input distribution prior). In other words, any information forgotten at any step cannot be later recalled. The uncertainty interval of size ≥ ε b/2 is not very severe, in the sense that it affects just one out of n input components. Indeed, at step 0 the input vector coordinates x(1), . . . , x(n) are represented in individual words, each of which gives rise to an uncertainty interval of width ε. This means that merely storing the input in memory in the standard coordinate system implies knowing its location up to an uncertainty n-cube with side ε, and of diameter ε √ n. 4 An uncertainty interval of size ε b/2 = O(ε √ log n) in a single direction is therefore hardly a deterrent for someone who wants to asymptotically beat FFT. The next section tells us, however, that the situation is much more severe, in the sense that underflow occurs in Ω(n) pairwise orthogonal input directions. (i) (Severe Overflow) There exists an orthonormal system v 1 , . . . , v n ′ ∈ R n , integers t 1 , . . . , t n ′ ∈ [m] and i 1 , . . . , i n ′ ∈ [n] with n ′ = Ω(n) such that for all j ∈ [n ′ ],

7)
where P j is projection onto the space orthogonal to v 1 , . . . , v j−1 .
In both cases of severe overflow and severe underflow, the gates at time t 1 , . . . t n ′ are all rotations, for all j ∈ [n ′ ] the index i j is one of the two indices affected by the corresponding rotation. Additionally, the set {t 1 , . . . , t n ′ } is of size at least n ′ /2.
The proof is deferred to Section 6.2. We discuss its implications, continuing the discussion following Theorem 5.1 in Section 5.2.
In the severe overflow case, Theorem 5.3 tells us that there exists an orthonormal collection v 1 , . . . , v n ′ (with n ′ = Ω(n)) in the input space, such that each v i behaves like x over (from Section 5.2). This means that, if the speedup factor b is ω(1), we have overflow caused by a linear number of independent in put components, occurring at Ω(n) different time steps (by the last sentence in the theorem). In the extreme case of speedup b = Θ(log n) (linear number of gates), this means that in a constant fraction of time steps overflow occurs.
The severe underflow is more interesting to analyze. An intuitive analysis is as follows: For each orthogonal direction v j , for j ∈ [n ′ ], we lose Θ(log b) bits of information about x T v j just before applying the t j 'th gate because, roughly speaking, we are trying to squeeze the value x T v j /γ j ∼ N (0, Θ(1/γ 2 j )) in the i j 'th word while not storing any information about x T v j in any other word, and γ j ≥ b/2. This lost information cannot be recovered later.
Although this explanation can be made rigorous using information theoretical tools, we prefer to offer a geometric interpretation. The theorem tells us that there exists an orthonormal collection u 1 , . . . , u n ′ in the input space that is bad in the following sense. For each j ∈ [n ′ ], redefine g j = u T j x to be the input component in direction u j . Again, the variables g 1 , . . . , g n ′ are iid N (0, Θ(1)). The first element in the series, u 1 , can be analyzed as x under (Section 5.2), whereby it was argued that before the t 1 'th step, the component g 1 = u T 1 x can only be known to within an interval of width Ω(γ 1 ε), independently of information from components orthogonal to u 1 . We remind the reader that by this we mean that the width of the interval is independent, but the location of the interval depends smoothly (in fact, linearly) on information from orthogonal components of x (see Figure 1).
As for u 2 , . . . , u n ′ : For each j ∈ [n ′ ], let z j := (M (t j ) ) −T (i j , :). Therefore u 1 = z 1 / z 1 and by (5.8), for j > 1 we can write z j = γ j u j + h j where h j ∈ span{u 1 , . . . , u j−1 }. Treating z j / z j as x under in Section 5.2, we conclude that the component (z j / z j ) T x of x can only be known to within an interval of size Ω(ε z j ), given any value of the projection of input x onto the space orthogonal to z. We extend the list of vectors z 1 , . . . , z n ′ , orthonormal vectors u 1 , . . . , u n ′ , numbers γ 1 , . . . , γ n ′ and projections Q 1 , . . . , Q n ′ to size n as follows. Having defined z j , u j , Q j , γ j for some j ≥ n ′ , we inductively define Q j+1 as projection onto the space orthogonal to span{z 1 , . . . , z j } = span{u 1 , . . . , u j } and z j+1 to be a standard basis vector such that Q j+1 z j+1 2 ≥ 1 − j/n. (Such a vector exists because there must exist an index i 0 ∈ [n] such that j j ′ =1 u j ′ (i 0 ) 2 ≤ j/n, by orthonormality of the collection u 1 , . . . , u j ; Now set z j+1 to have a unique 1 at coordinate i 0 and 0 at all other coordinates.) We let u j+1 be Q j+1 z j+1 / Q j+1 z j+1 , that is, a normalized vector pointing to the component of z j+1 that is orthogonal to span{z 1 , . . . , z j } = span{u 1 , . . . , u j }. The number γ j+1 is defined as Q j+1 z j+1 . By construction, γ j+1 ≥ 1 − j/n.
The above extends the partial construction arising from the severe underflow to a full basis, with the following property: Proposition 5.4. For any j ∈ n, even given exact knowledge of the exact projectionx of x onto the space orthogonal to z j , the quantity x T (z j / z j ) upon termination of the algorithm can only be known to within an interval of the form [s, s + ε z j ] where s depends smoothly (in fact, linearly) onx.
The proposition is simply a repetition of the analysis done for x under in Section 5.2, for j ≤ n ′ . For j > n ′ it is a simple consequence of the fact that upon initialization of the algorithm with input x, each coordinate of x (and in particular x T z j ) is stored in a single machine word, while all other machine words store information independent of x T z j . Hence the uncertainty of width ε z j = ε.
What do we know about x upon termination of the algorithm? As stated earlier, any information that was lost during execution, cannot be later recovered. Let I denote the set of possible inputs, given the information the we are left with upon termination. Consider the projection Q 2 onto the space orthogonal to u 1 = z 1 / z 1 , as a function defined over I. Let I 2 = Q 2 I denote its image. The preimage of any point w ∈ I 2 must contain a line segment of length at least εγ 1 parallel to u 1 , due to the uncertainty in x T u 1 . Hence the volume of I is at least εγ 1 times the (n − 1)-volume of I 2 . 5 Continuing inductively, we lower bound the (n − j + 1)-volume of I j := Q j I = Q j I j−1 for j > 2. Consider the projection Q j as a function operating on I j−1 , and any point w in the image I j . By definition of Q j , there existsŵ ∈ I such that Q jŵ = w. By proposition 5.4, the intersection of the line L = {ŵ + ηz j : η ∈ R} with I must contain a segment ∆ of size ε z j . The projection Q j ∆ of this segment is contained in the line Q j L = {w + ηu j : η ∈ R}. The size of the segment is ε Q j z j = εγ j . This means that the (n − j + 1)-volume of I j+1 is at least εγ j times the (n − j)-volume of I j+1 = Q j+1 I j .
Concluding, we get that the volume of I is at least n j=1 γ j . From the construction immediately preceding Proposition 5.4, we get (using the fact that n ′ = Ω(n)): The bound (5.9) tells us that the volume of uncertainty in the input (and hence, the output) of a b-speedup of FFT in the in-place model is at least b Ω(n) times the volume of uncertainty incurred simply by storing the input in memory.

Proofs
We start with some useful lemmas. Lemma 6.1. Let x, y ∈ R a for some integer a, with x 2 = y 2 = 1. Then − log a ≤ Φ(x, y) ≤ log a . (6.1) The proof is a simple done by a simple analysis of the function Φ(x, y) under the stated constraints using, say, Lagrange multipliers. Lemma 6.2. Let A, B ∈ R a×n . Let U ∈ R a×a be orthogonal. Then rs · x(j)y(j) log |rs · x(j)y By orthogonality of U , we have that a j=1 x(j)y(j) = a j=1 (U x)(j)(U y)(j). Also for the same reason we have U x 2 = U y 2 = 1. Using Lemma 6.1, we conclude Proof. Let U, V ∈ R a be orthogonal and Σ ∈ R a diagonal (and nonsingular) so that D = U ΣV .
(Such a composition exists by standard SVD theory.) as required. (We used Lemma 6.2 twice in the second inequality, and the orthogonality of U for the last derivation. The reason the middle term in the RHS of the first inequality is null is by properties of Φ that are trivial to check.) The following is the most important technical lemma in this work. Roughly speaking, it tells us that application of operators that are close to Id to the rows of F and F −T does not reduce the corresponding potential by much. Similarly, assuming that P, Q are PSD with spectral norm at most 1, applying these transformations to the rows of Id does not increase the corresponding potential by much. Lemma 6.4. Let P, Q ∈ R n×n be two matrices. LetP = Id −P,Q = Id −Q.Then If, additionally, P and Q are positive semi-definite contractions (the spectral norm of P, Q is at most 1), then The proof is deferred to Section 6.3. The bounds (6.8)-(6.9) use the smoothness of the matrices F and Id (that is, almost all matrix elements have exactly the same magnitude). This is the reason we needed to modify Φ and work with Φ C for the complex case: If F were the real representation of the n/2-DFT matrix, then it is not smooth in this sense. It does hold though that for any i ∈ [n] and j ∈ [n/2]: F (i, 2j − 1) 2 + F (i, 2j) 2 = 2/n, in other words, the matrix is smooth only in the sense that all pairs of adjacent elements have the same norm (viewed as vectors in R 2 ).

Proof of Theorem 5.1
We directly prove the less general (5.3). The more general bound (5.1) is shown similarly, but with more notation. Fix R ∈ [⌊n/2⌋]. Let m ′ be the smallest integer divisible by R satisfying m ′ ≥ m. If m ′ > m, then "pad" the algorithm A n by defining M (m+1) . . . M (m ′ ) = M (m) = F . By the triangle inequality, Combining (6.11) with (6.10), we get For any matrix A ∈ R n×n and any subset I ⊆ [n], we have A(I, : (this can be seen e.g. using the SVD theorem). Therefore, But Φ(M (m ′ ) ) = Φ(F ) = n log n and Φ(M (0) ) = 0, hence, there must exists 2Rn log n 2((n log n)/b) log 2R = Rb log 2R .

Proof of Theorem 5.3
Let V = {v 1 , . . . , v k } ⊆ R n be an orthonormal set satisfying the properties described in case (i) of the theorem. Similarly, let U = {u 1 , . . . , u ℓ } ⊆ R n be an orthonormal set satisfying the properties described in case (ii) of the theorem. We show that as long as k + ℓ is at most O(n), then we can extend one of the two sets by one element. Let P denote the projection onto (span V ) ⊥ and Q the projection onto (span U ) ⊥ . Using Lemma 6.4, we have that for some global C > 0, whereP = Id −P,Q = Id −Q (the orthogonal projections). By known properties of projection matrices, trP = P 2 F = k, trQ = Q 2 F = ℓ. Therefore, for some global C ′ > 0. This implies, using Theorem 5.1 (with R = 1) that for some t ∈ {0, . . . , m} and i ∈ [n]: Hence, as long as k + ℓ ≤ n/(2C ′ ): This implies that either M (t) P ⊥ (i, :) ≥ b/2 or that (M (t) ) −T Q ⊥ (i, :) ≥ b/2. In the former case we can extend the set V by adding v k+1 = M (t) P (i, :)/ M (t) P (i, :) , which is orthogonal to v 1 , . . . , v k by construction. In the latter case we can extend the set U by adding u k+1 = (M (t) ) −T Q(i, :)/ (M (t) ) −T Q(i, :) which is again, orthogonal to u 1 , . . . , u ℓ by construction.
This process of augmenting V and U can continue until k + ℓ ≥ n/2C ′ , which implies that either k ≥ n/4C ′ (establishing extreme overflow of the theorem) or ℓ ≥ n/4C ′ (establishing extreme underflow).
Assume k ≥ n/4C ′ and let n ′ = k. For j ∈ [n ′ ] let t j denote the time step of the overflow corresponding to direction v j ∈ R n and let i j denote the coordinate at which the overflow occurs. It is clear from the construction that i j is one of the at most two coordinates affected by the t j 'th step. We show that there exist no 1 ≤ j < j ′ ≤ n ′ such that (t j , i j ) = (t j ′ , i j ′ ). Indeed, note that M (t j ) P j ′ (i j , :) must be null by construction, contradicting the fact that . The conclusion is that for any t ∈ [m] there can be at most two indices j, j ′ ∈ [n ′ ] such that t j = t j ′ = t, and therefore the cardinality of the set {t 1 , . . . , t n ′ } is at least n/2. A similar argument is done for the extreme underflow case, concluding the proof.
For A, we just need to notice that: We now bound B.
This concludes the proof of (6.8).

Future Work
Taking into account bit operation complexity, and using state-of-the-art integer multiplication algorithms [8,9] it can be quite easily shown that both severe overflow and severe underflow could be resolved by allowing flexible word size, accommodating either large numbers (in the overflow case) or increased accuracy (in the underflow case). In fact, allowing O(log b)-bit words at the time steps at which overflow (or underflow) occur, of which there are Ω(n) many by Theorem 5.3, suffice. Hence, this work does not rule out the possibility of (in the extreme case of b = Θ(log n)) a Fourier transform algorithm in the in-place model using a linear number of gates, in bit operation complexity ofΩ(n log log n), whereÕ() here hides log log log n factors arising from fast integer multiplication algorithms . We conjecture that such an algorithm does not actually exist, and leave this as the main open problem.
Another problem that was left out in this work is going beyond the in-place model. In the more general model, the algorithm works in space R ℓ for ℓ > n, where the (ℓ − n) extra coordinates can be assumed to be initialized with 0, and the first n are initialized with the input x ∈ R n . The final matrix M (m) of Fourier transform algorithm A n = {Id = M (0) , . . . , M (m) } contains F as a sub matrix, so that the output F x can simply be extracted from a subset of n coordinates of M (m) x, which can be assumed to be the first. The matrix M (m) (and its inverse-traspose) therefore contains (ℓ − n) extra rows. The submatrix defined by the extra rows (namely, the last ℓ − n) and the first n columns were referred to in [2] as the "garbage" part of the computation. To obtain an Ω(n log n) computational lower bound in the model assumed there, 6 it was necessary to show that Φ P,P (M (m) )P = Ω(n log n), where P ∈ R ℓ×ℓ is projection onto the space spanned by the first n standard basis vectors. 7 To that end, it was shown that such a potential lower bound held as long as spectral norm of the "garbage" submatrices was properly upper bounded. That result, in fact, can be deduced as a simple outcome of Lemma 6.4 that was developed here. What's more interesting is how to generalize Theorem 5.3 to the non in-place model, and more importantly how to analyze the numerical accuracy implications of overflow and underflow to the non in-place model. Such a generalization is not trivial and is an immediate open problem following this work.
The closing question of Section 5 is another avenue of possible research stemming from this work, and seems (to the author) as an achievable goal in the quest to solve the half decade old problem of understanding the complexity of Fourier transform. In view of this work and of the preceding [1,2], it seems that the correct way to understand the hardness of this important transformation is by way of trading off time versus accuracy, for a given machine architecture.
Another interesting possible avenue is to study the complexity of Fourier transform on input x for which some prior knowledge is known. The best example is when F x is assumed sparse, for which much interesting work on the upper bound side has been recently done by Indyk et al. (see [10] and references therein).
Many algorithms use the Fourier transform as a subroutine. In certain cases (fast polynomial multiplication, fast integer multiplication [8,9], fast Johnson-Lindenstrauss transform for dimensionality reduction [3,4,5,11] and the related restricted isometry property (RIP) matrix construction [14,6,11]) the Fourier transform subroutine is the algorithm's bottleneck . Can we use the techniques developed here to derive lower bounds (or rather, time-accuracy tradeoffs) for those algorithms as well? Moreover, we can ask how the implications of speeding up the Fourier transform subroutine (as derived in this work) affect the numerical outcome of these algorithms, assuming they insist on using Fourier transform as a black box.