The Complexity of Computing a Fourier Perturbation
Nir Ailon\footnote{nailon@cs.technion.ac.il
The complexity of computing the Fourier transform is a longstanding open problem. Very recently, Ailon (2013, 2014, 2015) showed in a collection of papers that, roughly speaking, a speedup of the Fourier transform computation implies numerical ill-condition. The papers also quantify this tradeoff. The main method for proving these results is via a potential function called quasi-entropy, reminiscent of Shannon entropy. The quasi-entropy method opens new doors to understanding the computational complexity of the important Fourier transformation. However, it suffers from various drawbacks. This paper, motivated by one such drawback, eliminates it by extending the method. The argument goes as follows: If instead of computing the Fourier transform $Fx$ of input $x\in \R^n$ we were to compute a Fourier $\eps$-perturbation, defined as $(\Id+\eps F)x$, then the quasi-entropy method in the well-conditioned regime would, without any adjustments, lead to a linear algebraic lower bound of $\Omega(\eps^2 n \log n)$ many operations (counting additions and multiplications on scalar variables). Had this bound been matched by an algorithm, then we would have been able to extract $Fx$ in time $O(\eps^2 n\log n) $ by first computing $(\Id+\eps F)x$, then subtracting $x$ from the output and dividing the result by $\eps$. By taking $\eps=\Theta(1/\sqrt{ \log n)$ we could artificially drive the computation time toward the trivial linear time lower bound. Such a scheme would suffer on a real computer from numerical errors, but this could be avoided by extending the computer word size by only $\Theta(\log \eps^{-1) = \Theta(\log\log n)$ bits. The end result is a Fourier algorithm in running time $\tilde O(n\log\log n)$ (counting logical bit operations, and using fast integer multiplication). We generalize the quasi-entropy method so as to show that driving $\eps$ down does not allow such a free ride in case of the Walsh-Hadamard Fourier transform, and that the linear algebraic lower bound is, in fact, $\Omega((n\log n)/\log \eps^{-1)$. This exactly `cancels out' the numerical accuracy overhead. It also strengthens our belief that, roughly speaking, Fourier computation requires $\Omega(n\log n)$ time in a computational model that takes into account numerical accuracy and logical bit operations.
low_bound_fourier_perturbation.tex
% ----------------------------s------------------------------------
% AMS-LaTeX Paper ************************************************
% **** -----------------------------------------------------------

\documentclass[11pt]{article}
\usepackage{geometry}
\geometry{verbose,tmargin=3cm,bmargin=3cm,lmargin=1in,rmargin=1in}
\usepackage{color}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{graphicx}
% ----------------------------------------------------------------
\vfuzz2pt % Don't report over-full v-boxes if over-edge is small
\hfuzz2pt % Don't report over-full h-boxes if over-edge is small
% THEOREMS -------------------------------------------------------
\newtheorem{thm}{Theorem}[section]
\newtheorem{cor}[thm]{Corollary}
\newtheorem{conj}[thm]{Conjecture}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\theoremstyle{definition}
\newtheorem{defn}[thm]{Definition}
\newtheorem{fact}[thm]{Fact}
\theoremstyle{remark}
\newtheorem{rem}[thm]{Remark}
\numberwithin{equation}{section}
% MATH -----------------------------------------------------------
\newcommand{\norm}[1]{\left\Vert#1\right\Vert}
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\Real}{\mathbb R}
\newcommand{\eps}{\varepsilon}
\newcommand{\To}{\longrightarrow}
\newcommand{\BX}{\mathbf{B}(X)}
\newcommand{\A}{\mathcal{A}}

\newcommand{\yl}{{\check{y}}}
\newcommand{\yh}{{\hat{y}}}
\newcommand{\bl}{{\check{b}}}
\newcommand{\bh}{{\hat{b}}}

\newcommand{\NA}[1]{{\color{blue} #1}}
\newcommand{\GY}[1]{{\color{red} #1}}
\def\id{\operatorname{Id}}
\def\supp{\operatorname{supp}}
%\def\span{\operatorname{span}}
\def\sign{\operatorname{sign}}
\def\CT{{\cal C}}
% ----------------------------------------------------------------
\def\ent12{\hat \Phi_{|1,2}}
\def\R{\mathbb R}
\def\eps{\varepsilon}
\def\polylog{\operatorname{polylog}}
\def\errmat{{\mathcal E}}
\def\Ibig{{I_{\operatorname{big}}}}
\def\Ismall{{I_{\operatorname{small}}}}
\renewcommand{\cos}{\operatorname{cos}}
\renewcommand{\sin}{\operatorname{sin}}
\title{The Complexity of Computing a Fourier Perturbation}
\author{Nir Ailon\footnote{\texttt{nailon@cs.technion.ac.il}} \hspace{1cm}  Gal Yehuda \footnote{\texttt{gal2016@gmail.com}} \\ Department of Computer Science \\ Technion Israel Institute of Technology \\ Haifa, Israel }

%-----------------------------------------------------------------
\begin{document}

\maketitle

% ----------------------------------------------------------------
%\begin{abstract}
%\end{abstract}
\def\sign{\operatorname{sgn}}
\def\dim{n}
\def\C{\mathbb C}
\def\R{\mathbb R}
\def\P{\mathcal P}
\def\trace{\operatorname{tr}}
\def\diag{\operatorname{diag}}
\def\rank{\operatorname{rank}}
\def\F{{\mathcal F}}
\def\Id{\operatorname{Id}}
\def\f{{\hat f}}
\def\A{{\mathcal A}}
\def\Z{{\mathbb Z}}
% ----------------------------------------------------------------

\begin{abstract}

The complexity of computing the Fourier transform is a longstanding open problem.  Very recently, Ailon (2013, 2014, 2015) showed in a collection of papers that, roughly speaking, a speedup of the Fourier transform computation  implies numerical ill-condition.  The papers also quantify this tradeoff. 
The main method for proving these results is via a potential function called quasi-entropy, reminiscent of Shannon entropy.    The quasi-entropy method opens new doors to understanding the computational complexity of the important  Fourier transformation.  However, it suffers from various drawbacks.   
This paper, motivated by one such drawback, eliminates it by extending the method.   

 The argument goes as follows: If instead of computing the Fourier transform $Fx$ of input $x\in \R^n$ we were to compute a Fourier $\eps$-perturbation, defined as $(\Id+\eps F)x$, then the quasi-entropy method in the well-conditioned regime would, without any adjustments, lead to a linear algebraic lower bound of $\Omega(\eps^2 n \log n)$ many operations (counting additions and multiplications on scalar variables).  Had this bound been matched by an algorithm, then we would have been able to extract $Fx$ in time $O(\eps^2 n\log n) $ by first  computing $(\Id+\eps F)x$, then  subtracting $x$ from the output and dividing
the result by $\eps$.  By taking  $\eps=\Theta(1/\sqrt{ \log n})$  we could artificially  drive the  computation time toward the trivial linear time lower bound.  Such a scheme would suffer on a real computer from numerical errors, but this could be avoided by extending the computer word size by   only $\Theta(\log \eps^{-1}) = \Theta(\log\log n)$ bits.  The end result is a Fourier algorithm in running time $\tilde O(n\log\log n)$ (counting logical bit operations, and using fast integer multiplication).

We generalize the quasi-entropy method  so as to show that driving $\eps$ down does not allow such a free ride in case of the  Walsh-Hadamard Fourier transform, and that the linear algebraic lower bound is, in fact, $\Omega((n\log n)/\log \eps^{-1})$.  This exactly `cancels out' the numerical accuracy  overhead.  It also  strengthens our belief that, roughly speaking,  Fourier computation requires $\Omega(n\log n)$ time in a computational model that takes into account numerical accuracy and logical bit operations.
\end{abstract}

\pagebreak

\section{Introduction}\label{sec:intro}
The Fourier transform is one of the most important  linear transformations in science and engineering.
The (normalized) Discrete Fourier Transform (DFT) $\hat x\in\C^n$ for input signal $x\in \C^n$ is
defined by
$$ \hat x_i = n^{-1/2}\sum_{j=1}^n e^{-2\pi \iota (i-1)(j-1)/n} x_j\ ,$$
where $\iota = \sqrt{-1}$.   
 DFT has applications in many fields, including fast polynomial multiplication \cite[chapter~30]{Cormen:2009:IAT:1614191}, fast
integer multiplication \cite{Furer:2007:FIM:1250790.1250800}, fast large scale linear algebra and matrix sketching \cite{Avron:2010:BSL:1958627.1958633,sarlos2006improved},
signal processing \cite[chapters~ 6-9]{elliott2013handbook}, \cite{mallat1999wavelet,vetterli1984simple} and more.
From a theoretical perspective, the DFT is a special case  of the more general  Fourier transform
on  abelian groups,  with respect to the group  $\Z/n\Z$.   Another special case, known as the Walsh-Hadamard
transform, is defined over the group $(\Z/2\Z)^{\log_2 n}$ (for integer $\log_2 n$).
The Walsh-Hadamard transform $\hat x^{WH}$is given by 
$$ \hat x^{WH}_i = n^{-1/2} \sum_{j=1}^n (-1)^{\langle i-1, j-1\rangle} x_j\ ,$$
where $\langle a, b \rangle$ here is dot product of two bit vectors corresponding to the base-$2$ representation of
integers $a,b$.
The WH transform has applications in coding theory and digital image processing \cite{bayer1985image,berrouche2014improved,decker1970hadamard, lum2015fast,vetterli1984multi} as well as  in fast large scale
linear algebra and matrix sketching \cite{boutsidis2013improved,martinsson2010randomized}.  It is also instrumental in analysis of boolean functions
\cite{blais2009testing,blais2012testing,fischer2002testing} and more generally in theoretical computer science and learning theory \cite{de2008brief,mansour1994learning}.
%The main result in this work applies only to the WH transform.

From a computational point of view, an $O(n\log n)$ algorithm is known for both the DFT and the WH
transform.  For the DFT case, this was discovered by Cooley and Tukey in 1965 \cite{1965-cooley} in a seminal
paper.   For the WH case, the corresponding WH transform has been discovered in 1976 \cite{10.1109/TC.1976.1674569}.
Both algorithms run in a linear algebraic computational model (more on that in Section~\ref{sec:notation_model}).

The complexity of computing the Fourier transform is, on the other hand a longstanding open problem.  
An $\Omega(n)$ bound is trivial due to the necessity to consider all input coordinates.
The gap between $\Omega(n)$ and $O(n\log n)$ may seem small.  Nevertheless, owing to the
importance of the Fourier transform (both DFT and WH), it is crucial to close it in a reasonable model of computation.  Some early
results \cite{Morgenstern:1973:NLB:321752.321761} proved a lower bound of the \emph{unnormalized} Fourier transform,
defined as a scaling up of the Fourier transform by $n^{1/2}$,
using a potential function that is related to matrix determinant.  The result used the fact
that scaling up an $n$-by-$n$ matrix by a factor of $\alpha$ increases its potential by a factor
of $\alpha^n$, while any linear algebraic operation involving numbers of at most constant modulus (as used in FFT) can
increase the potential by a constant factor.  
 This result, though shedding light on an important problem, unfortunately does not explain why the normalized (unitary)  Fourier transform has complexity $\Omega(n\log n)$ and, conversely, does not explain
why the Fourier transform is computationally more complex than a scaling of the identity matrix by a
factor $n^{1/2}$.
A result by Papadimitriou \cite{Papadimitriou:1979:OFF:322108.322118} provides a lower bound
of computing the Fourier transform in finite fields, in a computational model that is not comparable
with ours. 

 Recently
Ailon \cite{DBLP:journals/cjtcs/Ailon13,DBLP:conf/icalp/Ailon15,DBLP:journals/toct/Ailon16} showed in a collection of papers that speedup of Fourier computation (both DFT and WH) implies
ill-conditioned computation (see also Section~\ref{sec:notation_model} below for a precise definition).   The result uses a potential 
reminiscent of Shannon entropy function on probability vectors, except that it is applied to 
any real vector (including e.g. negative numbers).   
%The method can also be used  over the complex
%field, but we prefer to avoid this complication by working over the reals only.  


%Henceforth,
%we will only consider the WH transform, which will be denoted by $F$.
%By FFT we mean, the fast WH transform.
%By DFT we shall refer to its embedding as a real $2n$-by-$2n$ matrix, acting on input vectors $x$
%in $\R^{2n}$ in  a standard way.

\subsection{The Quasi-Entropy Method}
In this work we work over the reals $\R$, and leave the DFT (which is complex) to future work.
We briefly remind the reader of matrix quasi-entropy:  For a nonsingular real $n$-by-$n$ matrix $M$, the matrix quasi-entropy $\Phi(M)$ is given as
\begin{equation}\label{eq:defphi}
 \Phi(M) = -\sum_{i=1}^n\sum_{j=1}^n M(i,j)M^{-T}(i,j)\log | M(i,j)M^{-T}(i,j)|\ .
 \end{equation}
where $M^{-T}$ is shorthand inverse-transpose.
The main result in \cite{DBLP:journals/toct/Ailon16} is, that as long as the condition number of $M$ is at most $\kappa$, a \emph{planar rotation} operation applied to a pair of rows of $M$ can change entropy by at most $O(\kappa)$. A  planar rotation matrix, or simply a rotation matrix for brevity, is obtained by multiplying
$M$ on  the left with the $(i,i')$-plane $\Theta$-rotation matrix $R_{i,i',\Theta}\in \R^{n\times n}$ defined by
the identity $\left ( \begin{matrix} R_{i,i',\Theta}(i,i)& R_{i,i',\Theta}(i,i') \\R_{i,i',\Theta}(i',i) & R_{i,i',\Theta}(i',i') \end{matrix}\right)=\left ( \begin{matrix} \cos\Theta & \sin \Theta \\ -\sin\Theta & \cos \Theta \end{matrix}\right)$, the remaining  diagonal elements  $1$, and then the remaining elements  $0$.

The quasi-entropy of the $n$-by-$n$ identity matrix  $\Id_n$ is $0$, that of the Fourier matrix
(both DFT and WH) is $\Omega(n \log n)$.  Therefore, in the uniformly $\kappa$-well conditioned
model (see Section~\ref{sec:notation_model} for an exact definition) a lower bound of
the number of steps is $\Omega(\kappa^{-1}n\log n)$. Equivalently, a speedup of FFT by factor of
$\kappa>1$ implies $\kappa$-ill conditioned computation.  In fact, as argued in followup work \cite{DBLP:conf/icalp/Ailon15},
the $\kappa$-ill condition implication is very strong in the sense that at $\Omega(n)$ pairwise orthogonal directions in input space, at some point along the computation,  $\Omega(\log \kappa)$ bits of information about each direction are lost due to overflow or underflow.   We refer the reader to  \cite{DBLP:conf/icalp/Ailon15} for the precise statement, however, we need to use here an extension of the quasi-entropy function that was defined there.  For two matrices $A,B$ with $n$ rows and $n'$ columns each, the $(A,B)$-preconditioned quasi entropy of $M$, denoted $\Phi_{A,B}(M)$, is defined as
\begin{equation}\label{eq:defprecond}
\Phi_{A,B}(M) = -\sum_{i=1}^n \sum_{j=1}^n (MA)(i,j) (M^{-T} B)(i,j) \log |(MA)(i,j) (M^{-T} B)(i,j)|\ .
\end{equation}
The reason we work with the preconditioned quasi-entropy function is that we will often track the change
of the entropy as $M$ evolves in the computation, while $A,B$ remain fixed.

%In followup work \cite{} it was shown that
%$R_{i,j,\Theta}(i,i)=R_{i,j,\Theta}(j,j)=\cos\Theta$, $
%$$ R_{i,j,\Theta} = \left ( 
%\begin{matrix}
%\Id_{i-1}     &                &                           &                      &  \\
 %                         &   \cos\Theta    &         &                    \sin\Theta &   \\
  %                        &                  &      \Id_{j-i-1}                   &                          &   &                      &   \\
 %                         &   -\sin\Theta    &         &                    \cos\Theta &   \\             
 %                         &                            &         &                                          & \Id_{n-j}       \\
%\end{matrix}\ ,
%\right )\ ,
%$$

\subsection{Problems with the Quasi-Entropy Method}

\begin{itemize}
\item[(P1)] {\bf  Suboptimality of speedup vs condition number lower bound}.
The main criticism of the method is the apparent sub-optimality of the derived lower bounds.
We provide an intuitive explanation to this point here, refer to \cite{DBLP:conf/icalp/Ailon15} for further details.
Roughly speaking,  the difficulty with a $\kappa$-conditioned algorithm for large $\kappa$ is, that
a fixed computer word size cannot accommodate the numbers needed for accurately carrying on the
computation.  In that sense, the results in \cite{DBLP:conf/icalp/Ailon15} tell us that speedup of FFT results in accuracy issues.
  But these issues  could be dealt with by computing with words of size $\Theta(\log \kappa)$, 
incurring a \emph{bitwise} time complexity cost of $ \Omega(\log \kappa)$ per operation.  By bitwise time
complexity we mean that we count logical binary operations, as measured in fast integer
multiplication algorithm design.   (See \cite{harvey2014even} for a groundbreaking recent result in that field, almost matching
the trivial $\Omega(\log \kappa)$ bound.)  If the lower bound of \cite{DBLP:journals/toct/Ailon16} were tight for $\kappa = \Theta(\log n)$, we could
speed up FFT by a factor of $\kappa$ to obtain a linear time algorithm (in the sense of counting linear algebraic steps), compensating for the loss of accuracy by emulating words
of size $\Theta(\log \log n)$ bits, resulting in bitwise running time of $\tilde O(n \log \log n)$.\footnote{Here, $\tilde O$ hides $o(\log \log n)$ factors, arising from the state-of-the-art integer-integer multiplicaiton algorithm \cite{harvey2014even}.}\footnote{We say \emph{emulating}
because, in a standard computer architecture the word size is constant, but we can emulate ``big numbers'' in software.}  Consequently, if  one believes (as we do) that the true bit complexity lower bound  of FFT computation should be $\Omega(n \log n)$, then the lower bound in \cite{DBLP:journals/toct/Ailon16} is probably not tight.
A major open conjecture in this line of work, which we will come back to in Section~\ref{sec:future} is the following tighter speedup vs condition number bound:
\begin{conj}[Informal]\label{conj:main} Speedup of FFT by factor of $\kappa>1$ (counting linear algebraic computation) can only be obtained using  $e^{\Omega(\kappa)}$-conditioned computation.
\end{conj}
\item[(P2)] {\bf Fourier perturbation complexity.}
%In this work we address a different shortcoming that might seem of lesser importance at first.
%Addressing this shortcoming, nevertheless, reveals that there is probably more `juice' in the quasi-entropy
%method than we currently understand.
   %Consider the following issue.
Instead of computing the Fourier
transform $Fx$ of input $x\in \R^n$,  assume we compute a Fourier $\eps$-perturbation thereof,
which we define as 
\begin{equation}\label{eq:defperturbation}
(\Id+\eps F) x\ .
\end{equation}
The corresponding matrix $(\Id+\eps F)$ is not orthogonal, but if $\eps$ is small then it is almost orthogonal,
in the sense that it has condition number $1+O(\eps)$.  In fact, it can be shown (see Appendix~\ref{sec:can_compute}) that a uniformly $(1+O(\eps))$-well conditioned 
algorithm (as defined in detail in Section~\ref{sec:notation_model})  \emph{can} be used to compute $(\Id+\eps F)$
and, for that matter, any $(1+O(\eps))$-conditioned matrix.
Now, it is easy to check that the quantity $\Phi(\Id+\eps F)$ is $-\Omega(\eps^2 n \log n)$.  
Hence, using the arguments in \cite{DBLP:journals/toct/Ailon16} without any change,
we get (in the WH case) a lower bound of 
\begin{equation}\label{eq:weakbound}
\Omega(\eps^2 n \log n)
\end{equation}
linear algebraic steps for computing a Fourier $\eps$-perturbation
in an algorithm of condition number $(1+O(\eps))$.
The problem is, that it is unlikely that (\ref{eq:weakbound}) is tight.   Indeed, assume that a matching algorithm existed.
Then we could compute $Fx$ by first computing $(\Id+\eps F)x$ in time $O(\eps^2n\log n)$
linear algebraic operations, subtracting $x$ from the output (after having stored a copy somewhere) and dividing the result by $\eps$.
Now we can drive $\eps$ down to $\Theta(1/\sqrt{\log n})$ so that the running time would match the trivial
lower bound of linear time.  To avoid numerical accuracy problems, we could enlarge the computer word
by $\Theta(\log 1/\eps) = \Theta(\log\log n)$ bits, incurring a $\tilde \Theta(\log \log n)$ complexity
for each addition or multiplication, taking into account logical bit operations.\footnote{Notation $\tilde \Theta$ here hides $\log\log \log n$ factors, coming from the current state-of-the-art machinery for fast integer multiplication \cite{furer2009faster}}.  The bottom line is an algorithm of running time $\tilde O(n\log\log n)$
(counting logical bit operations) for computing the Fourier transform.  We don't believe this is possible.
\end{itemize}

Notice that both problem (P1) and problem (P2) have a similar flavor. They both illustrate
that the quasi-entropy method without any adjustment is not strong enough to explain our
 belief that speeding up Fourier computation by some factor (counting linear
algebraic operations) must be offset by the overhead necessasry to maintain numerical accuracy.
Hence, we hope that solving one problem would shed light on the other.
\subsection{Our results}
We address problem (P2) in the WH case.
The main result in this work (Theorem~\ref{thm:main}) is an exponentially (in $\eps^{-1}$) improvement (over (\ref{eq:weakbound})) of the form
\begin{equation}\label{eq:mainineq}\Omega\left (\frac{n\log n}{\log \eps^{-1}}\right)\ .
\end{equation}
Note that  the term $\log \eps^{-1}$ in (\ref{eq:mainineq}) exactly offsets  the numerical overhead arising
from working with large numbers, as explained in (P2) above.
This leads us to believe that the improvement (\ref{eq:mainineq}) is optimal.
 The main technique used for deriving the improved bound
is a new extension of the quasi-entropy method.  Our  hope is that these techniques (in conjunction with others) could
be used to make progress with Conjecture~\ref{conj:main} in (P1).



% The main idea for `future work' item in this paper is to
%try to address problem (P1) by further understanding the quasi-entropy method.



\section{Notation and Computational Model}\label{sec:notation_model}
The results in this work apply to the WH transform only, which will henceforth be denoted by $F$.
By FFT we mean, the fast WH transform.
For the DFT case we believe that a similar method (or extension) should work -  see Section~\ref{sec:future}.
%The main result in this work applies both to WH and to the (real embedding of the) DFT, but for simplicity
%we will only assume the WH case.
%$F\in \R^{n\times n}$ will either denote the (real embedding of the)
 %$(n/2)$-DFT, or the WH transform.
%Both transformations have a uniform bound of $O(n^{-1/2})$ on the absolute
%value of their matrix elements.
For the WH case, all matrix elements are \emph{exactly} in $\{\pm n^{-1/2}\}$.
 Additionally, the WH transform is symmetric.   Hence $F=F^T, F^TF = F^2=\Id$. These algebraic facts make it easier to work with.
%Accordingly, FFT refers to the fast WH transform.
%
%Accordingly, FFT will either refer to the fast DFT transform (its usual meaning), or the fast WH transform.
%It is also an orthogonal matrix, and symmetric,  Hence, $F^2=\Id$. 

We recall the computational model in \cite{DBLP:journals/toct/Ailon16,DBLP:conf/icalp/Ailon15}.  We assume a linear algebraic computational model,
which means that the state of the machine at each step is a linear transformation of the initial state,
which is the input.
Therefore,  an algorithm running in $m$ steps is expressed as $\A=\{M^{(0)}=\Id, M^{(1)}\dots M^{(m)}\}$, 
where $M^{(t)}$ is the transformation taking the input $x$ to the $t$'th state of the machine.  
Two consecutive matrices
$M^{(t)}$ and $M^{(t-1)}$ differ in either one row indexed $i_t$ (a constant gate) or two rows indexed $i_t, i'_t$ (a rotation).  In case of a constant gate, 
$M^{(t)}$ is obtained by multiplying the  $i_t$'th row of $M^{(t)}$ by a constant $c_t\neq 0$.
In case of a rotation, $M^{(t+1)}$ is given as $R_{i_t, i'_t, \Theta_t} M^{(t)}$ for some angle constant
$\Theta_t$, where $R_{\cdot, \cdot, \cdot}$
was defined in Section~\ref{sec:intro}.
%rows $i_t, j_t$ of $ M^{(t)}$ are obtained by the corresponding rows in $M^{(t-1)}$ by a rotation of angle
%$\phi_t$.  [Exact definition of  a rotation]
We say  ``$\A$ computes $X$'' if $M^{(m)}=X$.

As argued in \cite{DBLP:conf/icalp/Ailon15,DBLP:journals/toct/Ailon16}, the computational
model is equivalent via simple reductions, to the so-called linear straight line computation (see
\cite{Ben-Or:1983:LBA:800061.808735} and references therein for a deifinition), except that
we  assume a \emph{no extra memory} scenario here.  This means that the machine
state at each moment is a vector in $\R^n$, represented using exactly $n$ computer words.  In particular,
 we are not allowed to pad the input with  zeros, somehow taking advantage of  the extra space throughout the computation.  Note that FFT  work in this model.
We refer the reader to \cite{DBLP:journals/toct/Ailon16} for an initial treatment of the \emph{additional memory} regime.  This paper shows that there is still much to uncover in the no extra memory regime, and we leave the extra memory
regime to a separate research branch.

Throughout $M^T$ is transpose of matrix $M$, and for invertible $M$,  $M^{-T}$ abbreviates $(M^{(-1)})^T$.
We also borrow from Matlab notation to create submatrices, e.g. $M([1\ 3\ 10], [5\ 6])$ is the submatrix obtained
from rows $1,3,10$ and columns $5,6$. Colon is used to specifiy the full range of indices, so e.g $M([1\ 3], :)$ is the
submatrix obtained by stacking rows $1$ and $3$ on top of each other.
For matrices $A,B$ with same number of rows, $[A,\,B]$ is the matrix obtained by stacking $B$ to the
right of $A$.
All logarithms are assume in base $2$.

The condition number $\kappa(M)$ of a nonsingular matrix $M$ is the ratio between its top and bottom singular values, equivalently, the ratio between its sepctral norm and that of its inverse.
Condition number is a standard measure used for quantifying numerical robustness of linear
algebraic computations, see e.g. chapter 11 of \cite{GolubV2013}.

\begin{defn} An algorithm $\A$ is (uniformly) $\kappa$-well conditioned (for $\kappa\geq 1$)  if for all $t\in [m]$ the condition number $\kappa(M^{(t)})$ of
the $t$'th matrix is bounded by $\kappa$.
\end{defn}

%Throughout we will assume that $\A$ is $(1+\eps)$-well conditioned.

%For any two matrices of same size $M$, $L$ (not necessarily square) we define the matrix bi-entropy function as follows:
%$$ \Phi(M, L) = -\sum M(i,j)  L (i,j) \log |M(i,j)L(i,j)|\ ,$$

%For  an invertible matrix $M$, the matrix pseudo-entropy of $M$  is defined \cite{} as $\Phi(M, M^{-T}) $, which is equivalently $$ -\sum M(i,j) (M^{-T})(i,j) \log |M(i,j) (M^{-T})(i,j)|\ .$$
%Abusing notation, and for the sake of compatibility with prior work,  we let $\Phi(M)$ be shorthand for the matrix pseudo-entropy of $M$.

%We also need the following useful definition \cite{}: For any two matrices of same size $M,L$ and another  pair $P,Q$ such
%that $MP$ and $LQ$ are defined (the number of columns of $M,L$ equals the number of rows of $P,Q$),
%$$
% \Phi_{P,Q}(M, L)  := \Phi(MP, LQ)\ .
%$$
%The definition is used for ease of notation, because often $P,Q$ will remain fixed while we focus our attention on the
%evolution of a sequence of pairs $(M,L)$.  The pair $P,Q$ will be often referred to as \emph{preconditioners}.
%If $M$ is invertible, then similarly to above:
%$$  \Phi_{P,Q}(M) := \Phi_{P,Q}(M, M^{-T})\ .$$

Recall the definitions of the quasi-entropy function $\Phi(M)$ and the preconditioned version $\Phi_{A,B}(M)$ from Section~\ref{sec:intro}.
\begin{thm}\label{thm:mainfromprev}\cite{DBLP:conf/icalp/Ailon15,DBLP:journals/toct/Ailon16}
If $M^{(t)}$ is obtained from $M^{(t-1)}$ by a constant gate, $\Phi(M^{(t)}) = \Phi(M^{(t-1)})$.
If $M^{(t)}$ is obtained from $M^{(t-1)}$ by a rotation acting on rows $i_t, i'_t$, then
\begin{equation}\label{eq:ineq}
 |\Phi_{P,Q}(M^{(t)}) - \Phi_{P,Q}(M^{(t-1)})| \leq \|(M^{(t-1)}P)([i_t\  i'_t], :)\|_F \cdot \|((M^{(t-1)})^{-T}Q)([i_t\ i'_t], :)\|_F\ ,
\end{equation}
 where $\|\cdot\|_F$ is Frobenius norm.
\end{thm}

Note that in particular, if $M^{(t-1)}$ (and hence, also $M^{(t)}$) is $\kappa$-well conditioned and $P,Q$ have both spectral norm $O(1)$, then
$$ |\Phi(M^{(t)}) - \Phi(M^{(t-1)})| =  O(\kappa)\ .$$

For $\eps<1/2$, a  Fourier $\eps$-perturbation is defined as the matrix $(\Id+\eps F)$.  
The inverse of the $\eps$-Fourier perturbation  is given by $\Id - \eps F + Z$, where $Z$ has spectral norm at most $O(\eps^2)$.
Using a simple entropy approximation lemma (deferred to Appendix~\ref{sec:applemma}),
it is not difficult to see that 
   the quasi entropy of a Fourier $\eps$-perturbation is
   \begin{equation}\label{eq:EpsPertPotential}\Phi(\Id+\eps F) = - \Omega(\eps^2 n\log n)\ .
   \end{equation}
%   For the DFT case, the potential is, in fact smaller in absolute value, but we shall not show this computation because the bound in the next Section is stronger, and applies to both WH and DFT.
   
Note that the potential is negative, but this will not matter for the purpose of establishing computational
lower bounds.
We will assume throughout that $1/\eps = n^{o(1)}$.  (Otherwise, computing with numerical accuracy $O(\eps)$ probably requires $\Omega(\log n)$ bits per word, which is not very interesting because then 
any multiplication or addition requires $\Omega(\log n)$ logical bit operations).


The condition number of a Fourier $\eps$-perturbation is clearly $1+O(\eps)$.  This means that
a necessary condition for an algorithm computing a Fourier $\eps$-perturbation is that it is not
$\kappa$ well conditioned for $\kappa = 1+o(\eps)$.  It turns out that
there exists an algorithm that is  $(1+O(\eps))$-well conditioned that computes $(\Id +\eps F)x$ given
input $x$. Refer to Appendix~\ref{sec:can_compute} for details of this simple fact.
By Theorem~\ref{thm:mainfromprev} we conclude:
\begin{prop}\label{prop:epssquare}
%In case $F$ is the WH transform,
 If algorithm $\A$ computes $(\Id+\eps F)$
 and $\A$ is  $(1+O(\eps))$-well conditioned, then $\A$ makes $\Omega(\eps^2 n \log n)$ steps.
\end{prop}

%\emph{Remark:}  The last warmup analysis applies to WH only; for DFT some adjustments are necessary.  In any case, the improved bounds below apply to both WH and DFT, as promised.
%We will need a slight generalization of the last theorem.  The generalization is obtained by noticing that the proof in \cite{}
%used Cauchy-Schwarz to derive (\ref{eq:ineq}) from the following stronger inequality:
%\begin{eqnarray}
 %|\Phi_{P,Q}(M^{(t)}) - \Phi_{P,Q}(M^{(t-1)})| &\leq& \sum_{k=1}^n \|(M^{(t)}P) ([i_t\ j_t], k)\|_2\cdot \|((M^{(t)})^{-T}Q) ([i_t\ j_t], k)\|_2 \ . \nonumber  \\
%\end{eqnarray}
%Taking a slightly different path from \cite{}, we get
%\begin{eqnarray}
 %|\Phi_{P,Q}(M^{(t)}) - \Phi_{P,Q}(M^{(t-1)})|&\leq& \sum_{k=1}^n \|(M^{(t)}P) ([i_t\ j_t], k)\|_1 \cdot \|((M^{(t)})^{-T}Q) ([i_t\ j_t], k)\|_1 \nonumber \\
%&\leq & 2 \|(M^{(t-1)}P)([i_t\  j_t], :)\|_1 \cdot \|((M^{(t-1)})^{-T}Q)t([i_t\ j_t], :)\|_\infty \label{eq:ineq1}
%\end{eqnarray}

%where the first inequality is obtained from the trivial  $\ell_1$ vs $\ell_2$ comparison, and the second  is the $(1,\infty)$ case of H\"older.  (Note that  we could use any $(p,q)$-H\"older inequality for dual $p,q$ norms, but
%we just need the trivial case here.  The more general inequality could be used in future work.)


%We will need a 
\section{From $\Omega(\eps^2 n \log n)$ to $\Omega(\eps n \log n)$} \label{bound1}

Assume  an algorithm $\A=(\Id=M^{(0)}, M^{(1)}, \dots, M^{(m)}=(\Id+\eps F))$ computes the Fourier $\eps$-perturbation, and that $\A$ is $(1+O(\eps))$-well conditioned.
The lower bound of $\Omega(\eps^2 n \log n)$  from Proposition~\ref{prop:epssquare} was obtained by tracking the following sequence of matrix quasi-entropies:

$$0 = \Phi(M^{(0)}), \Phi(M^{(1)}), \dots, \Phi(M^{(m)}) = \Phi(\Id + \eps F) = -\Omega(\eps^2 n \log n)\ .$$

A first step toward improving this bound is by identifying a pair of preconditioning matrices $A,B$
of at most constant spectral norm, for which
$$ \Phi_{A,B}(M^{(0)}) = \Phi_{A,B}(\Id) = \pm o(\eps  n\log n)\ \ \ \Phi_{A,B}(M^{(m)}) = \Phi_{A,B}(\Id+\eps F) = \pm\Omega(\eps n \log n)\ .$$

Using Theorem~\ref{thm:mainfromprev} (together with the bound on the spectral norm of $A,B$) implies that for all $t\in [m]$:
$$ \left | \Phi_{A,B}(M^{(t)}) - \Phi_{A,B}(M^{(t+1)}) \right | = O(1)\ .$$
Combining, the implication is that $m = \Omega(\eps n \log n)$.
%In order to get $\Omega(\eps n \log n)$, we track the following sequence:
It turns out that this can be achieved by taking $A=\Id$ and $B=F$. It is easy to check that
$$ \Phi_{\Id, F}(M^{(0)}) = \Phi_{\Id, F}(\Id) = O(\sqrt n \log n)\ .$$
%\begin{equation}\label{eq:chain}
%0=\Phi_{\Id, F}(M^{(0)}), \Phi_{\Id, F}(M^{(1)}), \dots, \Phi_{\Id, F}(M^{(m)}) = \Phi_{\Id, F}(\Id+\eps F)\ . %\Omega(\eps n\log n) .
%\end{equation}

We will now prove that $\Phi_{\Id, F} (\Id + \eps F) = \Omega(\eps n \log n) $ using the definition of $\Phi_{\Id, F} ()$ and a straightforward calculation. As before, recall  that $(\Id+\eps F)^{-1} = \Id-\eps F + Z$ where $Z$ has spectral norm $O(\eps^2)$.
%First, note that up to factor of $\eps^2$, it holds that $(\Id + \eps F)^{-T} = (\Id - \eps F)$. \NA{GAL: You need more details or notastion.  you can't just say ``up to factor $\eps^2$'' in a paper.  in a talk yes, in a paper no.}Hence: 
\begin{eqnarray}
\Phi_{\Id,F}(\Id + \eps F) &=& - \sum_{i = 1}^{n} \sum_{j=1}^n (\Id + \eps F)(i,j)\cdot (F - \eps F^2 + ZF)(j,i) \label{eq:doublesum} \\
& & \hspace{3cm} \log | (\Id + \eps F)(i,j)\cdot (F - \eps F^2 + ZF)(j,i)|  \ .\nonumber
\end{eqnarray}
%\Phi(\Id + \eps F, (\Id - \eps F)F) = \Phi(\Id + \eps F, F - \eps \Id)\ .$$

The diagonal summands in the last RHS are numbers in the set $\{x\log|x|: |x|=O(\eps + 1/\sqrt n)\}$.
Hence their total contribution to the sum is (in absolute value) 

\begin{equation}\label{eq:diagonal}O(n(\eps+1/\sqrt n)\log n)\ .\end{equation}

As for the non-diagonal elements, an intuitive analysis uses the fact that the dominant term
in $(\Id + \eps F)(i,j)\cdot (F + \eps F^2 + ZF)(j,i)$ is $\eps F(i,j)F(i,j)$.  All other terms have order $O(\eps^2)$.
For an exact analysis, refer to Appendix~\ref{sec:applemma}  
with $\ell=n-1$ to  bound  the inner sum of the RHS of (\ref{eq:doublesum}) (excluding the diagonal  $j=i$)
by $\Omega(\eps \log n)$.  Combining,
\begin{eqnarray}
\Phi_{\Id,F}(\Id + \eps F) &=& \Omega(\eps n \log n)\ ,
\end{eqnarray}
as required.
%
%We will denote $A = \Id + \eps F$, and $ B = F - \eps \Id$. By the definition of $\Phi$ we get :
%\begin{align}
%\Phi_{\Id, F}(\Id + \eps F, F - \eps \Id) &= - \sum_{i = 1}^{n} \sum _{j \neq i}A_{i,j}B_{i,j} \log |A_{i,j}B_{i,j}|   \nonumber \\
%&= - \sum_{i = 1}^{n} \sum _{j \neq i}\eps F_{i,j}^2 \log |\eps F_{i,j}^2| = - \sum_{i = 1}^{n}\eps \log (\frac{\eps}{n}) \nonumber \\
%&= \eps n \log (\frac{n}{\eps}) \nonumber \\ 
%& = \Omega(\eps n \log n)\ . \nonumber
%\end{align}
The method described in this section can probably not be used to improve on the lower bound
of $\Omega(\eps n \log n)$.  We need new ideas.


\section{From $\Omega(\eps n \log n)$ to $\Omega((n\log n)/\log \eps^{-1})$}

We state and prove the main result in this paper.
\begin{thm}\label{thm:main}
 If algorithm $\A$ computes $(\Id+\eps F)$
 and $\A$ is  $(1+O(\eps))$-well conditioned, then $\A$ makes $\Omega( (n \log n)/\log \eps^{-1})$ steps.
\end{thm}
\begin{proof}
In order to achieve the $\Omega((n \log n)/\log \eps^{-1})$ bound, we will use a different version of the preconditioned potential function. 
For a nonsingular matrix $M$ and two $n$-by-$2n$ matrices $P,Q$, we define
\begin{eqnarray}
\hat \Phi_{P,Q}(M) &=& -\sum_{i=1}^n \sum_{j=1}^n ((MP)(i,j)(M^{-T}Q)(i,j) +(MP)(i,j+n)(M^{-T}Q)(i,j+n) ) \nonumber \\
& & \hspace{3cm}+ \log |(MP)(i,j)(M^{-T}Q)(i,j) +(MP)(i,j+n)(M^{-T}Q)(i,j+n)| \nonumber 
%\hat \Phi_{P,Q}(X,Y) &=& \hat \Phi(XP, YQ)\ .
\end{eqnarray}  

\noindent
We will use $\hat \Phi$ in conjunction with the preconditioners:
\begin{align*}
P &= [\Id,\,-F] \\
Q &= [F,\,\Id]\ . \\
\end{align*}



%We will need the 
%For two $n \times 2n$ matrices $X, Y$ we define a variant of the quasi-entropy function, $\hat \Phi$, as:
%Given the above definition, the potential function is defined as follows:
%\begin{eqnarray}
%\Phi_{P^{(1)},P^{(2)}}(M^{(t)}) = \Phi(M^{(t)}[Id; -F], M^{(t)^{-T}}[F; Id]).
%\end{eqnarray}

\noindent
We aim to prove that a rotation can change the potential by at most $O(\eps \log \eps^{-1})$.
Let $t$ be such that $M^{(t+1)}$ is obtained from $M^{(t)}$ by a rotation. (It is clear to see that
a constant gate does not change the potential.)
Therefore we can assume that  $M^{(t+1)} = R_{i, i', \Theta}M^{(t)}$ for some row indices $i,i'$ and
angle $\Theta$.  Without loss of generality we can assume that $i=1, i'=2$, in other words that 
the rotation affects rows $i=1,2$ only.
First, one can notice that since $\kappa(M^{(t)}) = 1 + O(\eps)$, we can write 
\begin{align}
M^{(t)} &= U+\Delta \label{eq:MUDelta} \\
{(M^{(t)})}^{-T} &= U + \Gamma  \label{eq:MUGamma}
\end{align}
where $U$ is an orthogonal matrix, and $\Delta,\Gamma$ have spectral norm $O(\eps)$.

%$M^{(t)}$ as $U+\Delta$ where $U$ is a unitary matrix, and $\Delta$ is a matrix with spectral norm $O(\eps)$. With this notation,
%$$ M^{(t)}^{-1} = (U+\Delta)^{-1} $$
% $M^{(t)^{-T}} = U - \Delta$. Hence we get:
%\begin{eqnarray}
%\hat \Phi_{P,Q}(M^{(t)}, {(M^{(t)}}^{-T}) &=&\hat  \Phi((U+\Delta)[Id\ \  -F], (U+\Gamma)[F\ \  Id]) \nonumber \\
%&=& \hat \Phi\left([(U+\Delta)\ \ -(U+\Delta)F], [(U+\Gamma)F\ \ (U+\Gamma)]\right ) \label{eq:entprime}\ .
%\end{eqnearray}
%For simplicity of notation, let $A = U+\Delta$, $B = -(U+ \Delta)F$, $C = (U+\Delta')F$, $D = U + \Delta'$. 

%If $M^{(t+1)}$ is obtained from $M^{(t)}$ by a constant gate, then it is east to check that the potential doesn't change, namely that
%$\hat \Phi_{P,Q}(M^{(t)}) = \hat \Phi_{P,Q}(M^{(t+1)})$.

For the purpose of tracking the change in potential $\hat \Phi_{P,Q}$ we can ignore the contribution  to the potential coming from rows
$i>2$.
Let $V = UF, \Xi = \Delta F,  \Lambda = \Gamma F$, and $L(x) = x\log|x|$ for any real $x$.
If $\ent12$ denotes the contribution to $\hat \Phi_{P,Q}(M^{(t)})$ coming from rows $i=1,2$, then we can now express it as:
%Additionally for $i=1,2$ and $k=1..n$ let
%\begin{align}
%a(i,k) &= U(i,k) \alpha(i,k) &= \Delta(i,k) \\
%b(i,k) = (UF)(i,k) \beta(i,k) &= \Delta'(i,j) \\
%c(i,k)
%can write the contribution of these rows to the entropy $\hat \Phi_{P,Q}(M^{(t)})$ as
\begin{eqnarray}
&-&\sum_{j=1}^n L((U(1,j)+\Delta(1,j))(V(1,j) + \Lambda(1,j)) - (V(1,j)+\Xi(1,j))(U(1,j)+\Gamma(1,j))) \nonumber\\
&-& 
\sum_{j=1}^n L((U(2,j)+\Delta(2,j))(V(2,j) + \Lambda(2,j)) - (V(2,j)+\Xi(2,j))(U(2,j)+\Gamma(2,j)))\ . \nonumber 
\end{eqnarray}
Notice now that  the term $U(1,k)V(1,k)$ disappears from the first row, and $U(2,k)V(2,k)$ from the second.
Hence, %the last expression, which we now denote by $\ent12$, equals:

\begin{eqnarray}
\ent12=  &-&\sum_{j=1}^n L(  U(1,j)\Lambda(1,j) + V(1,j)\Delta(1,j) + \Delta(1,j)\Lambda(1,j) \label{eq:phi12}\\
& & \hspace{24ex}-    V(1,j)\Gamma(1,j) - U(1,j)\Xi(1,j)  - \Gamma(1,j)\Xi(1,j)) \nonumber\\
&-& 
\sum_{j=1}^n L(  U(2,j)\Lambda(2,j) + V(2,j)\Delta(2,j) + \Delta(2,j)\Lambda(2,j)  \nonumber \\
& & \hspace{24ex} -    V(2,j)\Gamma(2,j) - U(2,j)\Xi(2,j)  - \Gamma(2,j)\Xi(2,j))\ .\nonumber
\end{eqnarray}

Now let 
\begin{align}
r(j) &= \sqrt{U(1,j)^2 + V(1,j)^2 + U(2,j)^2 + V(2,j)^2} \label{eq:r}\\
\rho(j) &= \sqrt{\Delta(1,j)^2 + \Xi(1,j)^2 + \Delta(2,j)^2 + \Xi(2,j)^2 +\Gamma(1,j)^2 + \Lambda(1,j)^2 + \Gamma(2,j)^2 + \Lambda(2,j)^2 }\ . \label{eq:rho}
\end{align}

\noindent
Note that by our construction we have:
\begin{eqnarray}
\sum_{j=1}^n r(j)^2 = 4  & &
\sum_{j=1}^n \rho(j)^2 = O(\eps) \label{eq:boundrrho} 
\end{eqnarray}
where the left  identity is by orthogonality of $U,V$ and the right bound is by the spectral bound of $O(\eps)$
on $\Delta, \Gamma, \Delta, \Xi$ and $\Lambda$.
Dividing and multiplying by $r^2(j)$,  (\ref{eq:phi12}) now  becomes
\begin{eqnarray}
\ent12&=&\sum_{j=1}^n L\left (r^2(j)\left ( \frac{ U(1,j)\Lambda(1,j)}{r^2(j)} + \frac{V(1,j)\Delta(1,j)}{r^2(j)} + \frac{\Delta(1,j)\Lambda(1,j)}{r^2(j)} \right .  \right . \label{eq:dividemult}\\
& & \hspace{24ex} \left . \left .-    \frac{V(1,j)\Gamma(1,j)}{r^2(j)} - \frac {U(1,j)\Xi(1,j)}{r^2(j)}  - \frac{\Gamma(1,j)\Xi(1,j)}{r^2(j)} \right )\right ) \nonumber\\
&+& 
\sum_{j=1}^n L\left ( r^2(j) \left (\frac{U(2,j)\Lambda(2,j)}{r^2(j)} + \frac{V(2,j)\Delta(2,j)}{r^2(j)} + \frac{\Delta(2,j)\Lambda(2,j)}{r^2(j)}\right .  \right .  \nonumber \\
& & \hspace{24ex} \left . \left . -   \frac{ V(2,j)\Gamma(2,j)}{r^2(j)} - \frac{U(2,j)\Xi(2,j)}{r^2(j)}  - \frac{\Gamma(2,j)\Xi(2,j)}{r^2(j)}\right )\right ) \nonumber
\end{eqnarray}

\noindent
For simplicity of notation, for rows $i=1,2$  and any column $j$ we define:
\begin{eqnarray}
X(i,j) &=&  \frac{ U(i,j)\Lambda(i,j)}{r^2(j)} + \frac{V(i,j)\Delta(i,j)}{r^2(j)} + \frac{\Delta(i,j)\Lambda(i,j)}{r^2(j)} \label{eq:Xij}\\
 &-&    \frac{V(i,j)\Gamma(i,j)}{r^2(j)} - \frac {U(i,j)\Xi(i,j)}{r^2(j)}  - \frac{\Gamma(i,j)\Xi(i,j)}{r^2(j)}. \nonumber
\end{eqnarray}

Now recall that $M^{(t+1)}$ is obtained from $M^{(t)}$ by a  rotation matrix $R_{1,2,\Theta}$, affecting the first two rows only.  Therefore, from (\ref{eq:MUDelta}) and (\ref{eq:MUGamma}) we have
\begin{align}
M^{(t+1)} &= U'+\Delta' \nonumber \\ % R_{1,2,\Theta}U+R_{1,2,\Theta}\Delta \label{eq:MUDelta'} \\
{(M^{(t+1)})}^{-T} &= U' + \Gamma' \nonumber %R_{1,2,\Theta}U + R_{1,2,\Theta}\Gamma  \label{eq:MUGamma'}
\end{align}
where $U'=R_{1,2,\Theta}U$, $\Delta' = R_{1,2,\Theta}\Delta$, $\Gamma' = R_{1,2,\Theta}\Gamma$.
Similarly, we define $V' = R_{1,2,\Theta}V$, $\Xi' = R_{1,2,\Theta}\Xi$, $\Lambda' = R_{1,2,\Theta}\Lambda$
as the `post-rotation' version of the corresponding variables.  We can  also define $r'(j), \rho'(j)$ 
similarly to (\ref{eq:r}) and (\ref{eq:rho}), but with the  post-rotation variables.  However  clearly
$r(j) = r'(j)$ and $\rho(j)=\rho'(j)$ because rotation is an isometry, so $U(1,j)^2+U(2,j)^2 = U'(1,j)^2+U'(2,j)^2$ and similarly
for the other components in the RHS's of (\ref{eq:r}) and (\ref{eq:rho}).
%we will consider the case where the rotation is performed on each matrix separately, and denote by $U'(1,k)$ the matrix after the rotation (same for all matrices). 
%We will now analyze how rotation influences on the potential. 
%First, it can be seen that $r(k)$ remains the same: we can think of 
The ultimate goal is to compare the corresponding potentials $\ent12$ and $\ent12'$, where $\ent12'$
is obtained as in (\ref{eq:dividemult}), but using the post-rotation variables.

%$(U(1,k),U(2,k)), (V(1,k),V(2,k))$ as two vectors in $\R^2$, hence their norm will not change after rotation, which means $U(1,k)^2 + U(2,k)^2 = U'(1,k)^2 + U'(2,k)^2$ and $V(1,k)^2 + V(2,k)^2 = V'(1,k)^2 + V'(2,k)^2$, so $r'(k)=r(k)$. 
Now consider the expression $X(1,j) + X(2,j)$ for fixed $j$.  This expression can be viewed as sum of (scaled) inner products.  For example, the first inner product is 
$$\frac{\left \langle (U(1,j),U(2,j)),\ (\Lambda(1,j),\Lambda(2,j))\right \rangle}{r^2(j)}\ .$$ 
Hence, 
\begin{equation}\label{eq:equivXs} X(1,j)+X(2,j) = X'(1,j)+X'(2,k)\ ,\end{equation} where $X'(i,j)$ is obtained as in (\ref{eq:Xij}), but using the post-rotation variables.  Indeed planar inner products are not affected by a planar rotation.  
We are now finally in a position to compare $\ent12$ with $\ent12'$.
%is not affected by rotation: $X(1,k) + X(2,k) = X'(1,k) + X'(2,k)$. Therefore :
\begin{eqnarray}
|\ent12 - \ent12'| &=& \left |\sum_{j=1}^n (r(j)^2X(1,j)\log(|r(j)^2X(1,j)|)+ r(j)^2X(2,j)\log(|r(j)^2X(2,j))|) \right . \\
& &- \left .  \sum_{j=1}^n (r(j)^2X'(1,j)\log(r(j)^2X'(1,j))+ r(j)^2X'(2,j)\log(r(j)^2X'(2,j)) )\right |\nonumber \\
&=&\left  |\sum_{j=1}^n (r(j)^2\log(|r(j)^2|) (X(1,j)+X(2,j)-X'(1,j)-X'(2,j)) \right .
 \\   
& &+ \left .\sum_{j=1}^n (r(j)^2 (L(X(1,j)) + L(X(2,j)) - L(X'(1,j)) - L(X'(2,j)))\right | \nonumber \\
&=& \left |\sum_{j=1}^n (r(j)^2 (L(X(1,j)) + L(X(2,j)) - L(X'(1,j)) - L(X'(2,j)))\right | \nonumber \\
&\leq& \sum_{j=1}^n r(j)^2\left |(L(X(1,j)) + L(X(2,j)) - L(X'(1,j)) - L(X'(2,j))\right | \label{eq:last} \ ,
\end{eqnarray}
where the first equality is application of the rule $\log(xy)=\log x + \log y$, the second is from (\ref{eq:equivXs}) and the inequality is application of the triangle inequality.
%equality 
Clearly, $|U(i,j)|, |V(i,j)| \leq r(j)$ and $|\Lambda(i,j)|, |\Delta(i,j)|, |\Gamma(i,j)|, |\Xi(i,j)| \leq \rho(j)$ for $i=1,2$
and $j\in[n]$.  Therefore by definition of $X(i,j)$,  $|X(i,j)| \leq 6 \rho(j)/r(j)$.  A similar bound
holds for $X'(i,j)$.  For fixed $j$, the inner sum of  (\ref{eq:last}) is hence bounded above by
$$ r(j)^2 \left (2\max_{|x| \leq 6\frac{\rho(j)}{r(j)} } L(x) - 2\min_{|y| \leq 6\frac{\rho(j)}{r(j)} } L(y) \right ) \ .$$
%% $L(X(1,k)) \leq L(C\frac{\rho(k)}{r(k)})$ for some constant $C$. 
%Hence, the above expression is bounded by order of :
%\begin{eqnarray} \label{eq:sum1}
%\sum_{k=1}^n r(k)^2 (\max (L(\frac{\rho(k)}{r(k)})) - \min (L(\frac{\rho(k)}{r(k)}))) 
%\end{eqnarray}\label{eq:sum1}

Let $e$ be the natural logarithm basis.
 For nonnegative $a$ such that $a \leq e^{-1}$, 
$\max_{x:|x|\leq a} L(x) = -|a|\log|a|$ and $\min_{x:|x|\leq a} L(x) =|a|\log|a|$. For $a > e^{-1}$,   $\max_{|x|\leq a} L(x) \leq a(3+\log a)$ and $\min_{|x|\leq a} L(x) \geq -a(3+\log a)$.
We will now define the following subsets of column indices, indexed by an integer $h$:
\begin{align*}
J_h = \begin{cases} 
\left  \{j: 6\frac{\rho(j)}{r(j)} \leq \eps\right \}   & h=0 \\
\left  \{j : 2^{h-1}\eps \leq 6\frac{\rho(j)}{r(j)} \leq \min\{2^{h}\eps, e^{-1}\} \right  \} & 
1\leq h \leq \lceil -\log (\eps e)\rceil  \\
\left  \{j : \max\{e^{-1}, 2^{h-1}\eps\} \leq 6\frac{\rho(j)}{r(j)} \leq 2^{h}\eps \right  \} & 
 h >  \lceil -\log (\eps e)\rceil \\
\end{cases}\ .
\end{align*}
Splitting  the sum (\ref{eq:last}) and applying our analysis of the function $L(x)$, we get
\begin{eqnarray} 
|\ent12 - \ent12'| &\leq&  \sum_{h\geq 0}\sum_{j\in J_h} r(j)^2 \left (2\max_{|x| \leq 6\frac{\rho(j)}{r(j)} } L(x) - 2\min_{|y| \leq 6\frac{\rho(j)}{r(j)} } L(y) \right ) \nonumber \\ 
&\leq & -4\sum_{j\in J_0}r(j)^2\eps\log \eps - 4\sum_{h=1}^{\lceil -\log (\eps e)\rceil}\sum_{j\in J_h} r(j)^2  2^h \eps\log (2^h\eps) \nonumber \\
& & +4\sum_{h > \lceil -\log (\eps e)\rceil}\sum_{j\in J_h} r(j)^2\, 2^h \eps (3+\log (2^h \eps))\ .\label{eq:last2}
\end{eqnarray}

Let us now bound $\sum_{j\in J_h} r(j)^2$ for all $h$.  For $h=0$ we can simply recall  (\ref{eq:boundrrho}) to get the trivial 
\begin{equation}\label{eq:boundsumrj1}
\sum_{j\in J_0} r(j)^2\leq 4\ .
\end{equation}
%This is obvious from the definition (\ref{eq:r}), together with the fact
%that $U,V$ are orthogonal matrices. 
 As for $h\geq 1$,
from the definition of $J_h$, we get that for all $j\in J_h$: 
$$r^2(j) \leq \frac{36\rho^2(j)}{2^{2h-2}\eps^2}\ .$$ 
Therefore,
\begin{equation}\label{eq:boundsumrj2}\sum_{j\in J_h} r^2(j) \leq \sum_{j=1}^n \frac{36\rho^2(j)}{2^{2h-2}\eps^2} \leq C2^{-2h}\ ,\end{equation}
where the rightmost bound is by (\ref{eq:boundrrho}) and $C$ is proportional to the constant hiding in the $O()$-notation
there.  Plugging the bounds (\ref{eq:boundsumrj1}) and
(\ref{eq:boundsumrj2}) in (\ref{eq:last2}) gives:
\begin{eqnarray} 
|\ent12 - \ent12'| &\leq& -4\eps\log \eps -4 \sum_{h=1}^{\lceil -\log (\eps e)\rceil} C\,2^{-2h}  2^h \eps\log  (2^h \eps)  \nonumber\\
& & + 4\sum_{h > \lceil -\log (\eps e)\rceil}C\,2^{-2h} 2^h \eps (3+\log (2^h \eps)) \nonumber \\
&=& O(-\eps \log \eps)\ .
\end{eqnarray}
In order to complete the proof, we will show that $\hat{\Phi}_{P,Q}(Id)=0$ and $\hat{\Phi_{P,Q}(Id+\eps F)} = O(\eps n \log n)$.
\begin{eqnarray}
\hat \Phi_{P,Q}(Id) &=&-\sum_{i=1}^n \sum_{j=1}^n ((P)(i,j)(Q)(i,j) +(P)(i,j+n)(Q)(i,j+n) ) \nonumber \\
& & \hspace{3cm}+ \log |(P)(i,j)(Q)(i,j) +(P)(i,j+n)(Q)(i,j+n)| \nonumber 
\end{eqnarray}
By the definition of $P$ and $Q$, it holds that $P_{i,j}Q_{i,j} = -P_{i,j+n}Q_{i,j+n}$ and hence $\hat{\Phi}_{P,Q}(Id) = 0$. 
We now prove that  ${\hat \Phi_{P,Q}(Id +\eps F)} = \Omega(\eps n \log (n))$ . We remind the reader that $(Id + \eps F)^{-1} = \Id - \eps F + Z$ where $Z$ has spectral norm $O(\eps^2)$.
\begin{eqnarray}
& &\hat \Phi_{P,Q}(\Id + \eps F) = \label{eq:lastpotential}\\
& & \hspace{1cm} -\sum_{i=1}^n \sum_{j=1}^n ((\Id+\eps F)(i,j)(F-\eps F^2 +Z^TF)(i,j) +(-(F+\eps F^2))(i,j)(\Id - \eps F +Z^T)(i,j) ) \nonumber \\
& & \hspace{1cm}+ \log |((\Id+\eps F)(i,j)(F-\eps F^2+Z^TF)(i,j) +(-(F+\eps F^2))(i,j)(\Id - \eps F +Z^T)(i,j) )| \nonumber 
\end{eqnarray}
The diagonal elements in the last sum are of the form $\{x\log |x|: |x| = O(\eps + \frac{1}{\sqrt{n}})  \}$. Therefore their contribution to the sum is $O(n(\eps + \sqrt{n}^{-1})\log n)$. 

The main thing to note about the off-diagonal elements in the sum (\ref{eq:lastpotential}), as complicated as it may seem, is that products of the form $\Id(i,j)F(i,j)$ always appear  together with both a $+$ and a $-$ sign, hence cancel each other out.  The bottleneck is elements multiplied by $\eps$ exactly once. (Note that our bound on the spectral norm of $Z$ involves $\eps^2$,  hence elements such as $Z(i,j)$ or $(ZF)(i,j)$ are not part of the bottleneck.)  
These bottleneck elements give the theorem's
required bound.  The remaining elements are `noise', which can be accounted for using 
Lemma~\ref{lem:applemma} in Appendix~\ref{sec:applemma} for each row $i$.  To apply the lemma, take $\ell = n-1$ (accounting for all elements in  the row except diagonal), $x$'s coordinates are the elements
$2\eps F(i,j)^2$ (for $j\neq i$) and the vector $y$ is given by collecting the remaining elements, namely:
$$ -\eps^2 F(i,j) (F^2)(i,j) + \eps F(i,j)(Z^TF)(i,j) -F(i,j)Z(i,j) - \eps^2 (F^2)(i,j)F(i,j) - \eps(F^2)(i,j)Z^T(i,j)\ .$$
%To use the lemma at each row $i$, take $x$ (from the lemma's notation) to be the vector $(F(i,j))_{j\neq i} \in \R^{\ell }$,  $\ell$ to be $n-1$ and $y,z$ 
\noindent
The resulting estimation is:
\begin{eqnarray}
\hat \Phi_{P,Q}(\Id + \eps F) = \Omega(\eps n \log (n/\eps))\ . 
\end{eqnarray}
\noindent
(Recall that we assumed $1/\eps$ is $n^{o(1)}$, hence $\log(n/\eps) = \Theta(\log n)$.)
Combining this with the bound  $O(-\eps\log \eps)$ on the change of $\hat \Phi$ at each step
concludes the proof of the theorem.
This concludes the proof of the theorem.
\end{proof}

\section{Future Work}\label{sec:future}

The main technique used here was extending the quasi-entropy function.  The extension was `taylored'
for the Fourier (WH) $\eps$-perturbation, but we believe that a similar analysis is doable for DFT as well.
In fact, one can significantly extend the method giving rise to much
 freedom in proving bounds for other interesting problems.
A reasonable way to further extend the quasi-entropy function is to define

\begin{eqnarray*}
\Phi^{(k)}_{[A_1,\dots,A_k],\ [B_1,\dots, B_k]}(M)&=& \sum_{i=1}^n\sum_{j=1}^n \left(\sum_{p=1}^k (MA_p)(i,j)(M^{-T}B_p)(i,j)\right)  \\
& & \hspace{2cm} \log \left | \sum_{p=1}^k (MA_p)(i,j)(M^{-T}B_p)(i,j)\right | \ ,
\end{eqnarray*}
where $M$ is nonsingular and the $A_p$'s, $B_p$'s are square $n\times n$. The function $\hat \Phi$ we used corresopnds exactly to to $k=2$.

 The main hope is that this
technique, or some further extension, could be used to make progress on Conjecture~\ref{conj:main}.
Another open problem is to match the lower bound of Theorem~\ref{thm:main}, using a $(1+O(\eps))$-well conditioned algorithm (without extra memory).  We show in the appendix that a quadratic time algorithm exists, leaving quite a large gap.  




%In addition, from the definition of $\rho(k)$ we know that: $$\sum_{k=0}^n \rho^2(k) = 8\eps^2$$ Notice now that for each $k$, $\rho^2(k) \geq 0$, therefore each $\rho^2(k) \leq 8\eps^2$, which means that each $r^2(k) \leq \frac{8}{2^{2h}}$. Hence \ref{sum3} equals to: 
%\begin{eqnarray}
%\sum_{k\in I_{h}^1} \frac{8}{2^{2h}} (2^h\eps \log(2^{-h}) + 2^h\eps \log(\eps^{-1})) = 8\eps \sum_{h=0}^{\infty}\frac{2^h\log(2^{-h})}{2^{2h}} + 8\eps \log (\eps^{-1}) \sum_{h=0}^{\infty}\frac{2^h}{2^{2h}} \nonumber
%\end{eqnarray}
%\sum_{h\geq 0} 
%r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)})) \label{sums:1} \\
%+ \sum_{k\in I_{h}^2} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)})) \label{sums:2} \\
%+\sum_{k\in I_{h}^3} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)}))
%\end{eqnarray}
%

%J^0 &= \left  \{j: 6\frac{\rho(j)}{r(j)} \leq \eps\right \}   &
%J_{h}^1 &=\left  \{j : 2^{h}\eps \leq 6\frac{\rho(k)}{r(k)} \leq \min\{2^{h+1}\eps, e^{-1}\} \right  \}, h \leq \lceil -\log \eps e\rceil-1  \nonumber \\
%I_{h}^2 &= \{k : \frac{1}{e} \leq 2^{h}\eps \leq \frac{\rho(k)}{r(k)} \leq 2^{h+1}\eps  \} &
%I_{h}^3 &= \{k : 2^{h}\eps \leq \frac{\rho(k)}{r(k)} \leq 2^{h+1}\eps, 2^{h}\eps < \frac{1}{e}, 2^{h+1}\eps > %\frac{1}{e}  \} \nonumber 
%\end{align}
%% Therefore, it makes sense to define $$\Ibig = \left \{j\in[n]: 6\frac{\rho(j)}{r(j)} > e^{-1}\right \},\ \  \Ismall = [n] \setminus \Ibig\ .$$

%Let us consider the function $L(x)= x \log (|x|)$.
%, where we formally define $L(0) = 0$. 
%On a closed interval $[a,2a]$ such that $0 \leq a \leq \frac{1}{e}$, the function is monotonically decreasing, thus the maximum is obtained at $x = a$ and the minimum is obtained at $x = 2a$. It follows that on such interval $[a,2a]$: 
%\begin{eqnarray} \label{b1}
%\max (L(x)) - \min(L(x)) = O(a\log(a^{-1})). \label{b1}
%\end{eqnarray} \label{b1}
%On a interval $[a,2a]$ such that $a > \frac{1}{e}$, the function is monotonically increasing, hence the maximum is obtained at $x = 2a$ and the minimum is obtained at $x = a$, which means:
%\begin{eqnarray} 
%\max (L(x)) - \min(L(x)) = O(a\log(a)). \label{eq:bound2}
%\end{eqnarray}
%The last case is where $a < \frac{1}{e}$ but $2a > \frac{1}{e}$. In that case, the minimum is $\frac{1}{e}$, and the maximum is $\max \{a\log(a), 2a\log(2a)\}$, which is $O(a\log(a))$. \\


%Each $h \in \mathbb{N}$ can be thought as a \textit{layer}, such that for each layer $h$, $I_{h}^i$ (where $i = 1,2,3$) represents the $k's$ in that layer. Instead of summing over $k's$, we can sum over \textit{layers}, and \ref{eq:sum1} becomes:
%\begin{eqnarray} \label{sums}
%\sum_{k\in I_{h}^1} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)})) \label{sums:1} \\
%+ \sum_{k\in I_{h}^2} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)})) \label{sums:2} \\
%+\sum_{k\in I_{h}^3} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)}))
%\end{eqnarray}

%Let us bound \ref{1}. From  \ref{b1} we know that: 
%\begin{eqnarray}\label{sum3}
%\sum_{k\in I_{h}^1} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)})) &\leq & \sum_{k\in I_{h}^1} r^2(k) (2^h\eps \log(\frac{1}{2^h\eps})) \nonumber \\
%&\leq & \sum_{k\in I_{h}^1} r^2(k) (2^h\eps \log(2^{-h}) + 2^h\eps \log(\eps^{-1})). \label{sum3}
%\end{eqnarray}\label{sum3}
%From the definition of $I_{h}^{1}$, we get that $r^2(k) \leq \frac{\rho^2(k)}{2^{2h}\eps^2}$. In addition, from the definition of $\rho(k)$ we know that: $$\sum_{k=0}^n \rho^2(k) = 8\eps^2$$ Notice now that for each $k$, $\rho^2(k) \geq 0$, therefore each $\rho^2(k) \leq 8\eps^2$, which means that each $r^2(k) \leq \frac{8}{2^{2h}}$. Hence \ref{sum3} equals to: 
%\begin{eqnarray}
%\sum_{k\in I_{h}^1} \frac{8}{2^{2h}} (2^h\eps \log(2^{-h}) + 2^h\eps \log(\eps^{-1})) = 8\eps \sum_{h=0}^{\infty}\frac{2^h\log(2^{-h})}{2^{2h}} + 8\eps \log (\eps^{-1}) \sum_{h=0}^{\infty}\frac{2^h}{2^{2h}} \nonumber
%\end{eqnarray}
%The above sums convergence, and we get that (4.10) = $O(\eps \log (\eps^{-1}))$. \\
%Using the same technique, we will bound (4.11): 
%\begin{eqnarray}
%\sum_{k\in I_{h}^2} r^2(k) (\max(\frac{\rho(k)}{r(k)})-\min(\frac{\rho(k)}{r(k)})) &\leq & \sum_{k\in I_{h}^2} r^2(k) (2^h\eps \log(2^h\eps)) \nonumber \\
%&\leq & \sum_{k\in I_{h}^1} r^2(k) (2^h\eps \log(2^{h}) + 2^h\eps \log(\eps)). \nonumber
%\end{eqnarray}
%From the definition of $I_{h}^2$, we know that $\eps \geq \frac{1}{2^he}$, which means that (4.11) is bounded by: 
%\begin{eqnarray}
%\sum_{h=0}^{\infty} (\frac{\log(2^h)}{2^{2h}} + \frac{\log(2^{-h})}{2^{2h}}) = C.
%\end{eqnarray}
%For some constant $C$. \\
%$I_{h}^3$ can be bounded by exactly the same methods. \\
%The above leads to the conclusion that the entropy can change by at most $O(\eps \log \eps^{-1})$ after rotation. 




%Delta(1,k)^2 + V(1,k)^2 + U(2,k)^2 + V(2,k)^2}
%Let $A_{i,k} = a_{i,k} + \alpha_{i,k}, B_{i,k} = b_{i,k} + \beta_{i,k}, C_{i,k} = c_{i,k} + \gamma_{i,k}, D_{i,k} = d_{i,k} + \delta_{i,k}$, where $a_{i,k}, b_{i,k},c_{i,k}, d_{i,k} $ comes from the unitary part, and $\alpha_{i,k}, \beta_{i,k},\gamma_{i,k}, \delta_{i,k} $ comes from the part with the norm of $\eps$. By the construction, $a_{i,k}c_{i,k}+b_{i,k}d_{i,k} = 0$ for $i = 1,2$ and $k=0...n$.   



\bibliographystyle{plain}

\bibliography{low_bound_fft}


\appendix
\section{A Useful Lemma on Entropy  with Noise}\label{sec:applemma}
%\begin{lem}\label{lem:applemma}
%For any $C\geq 1$ there exists $\ell_0$ such that for all $\ell>\ell_0$ the following holds:
%For any $x\in \R^{\ell}$  such that $\|x\|_2=1$ and $|x(i)| \leq C/\sqrt \ell$ for
%all $i\in[\ell]$ and any $y\in \R^{\ell}$  such that $\|y\|_2^2\leq 1/(30C)^2$\ ,
%\begin{equation}\label{eq:applemma} -\sum_{i=1}^\ell  x_i(x_i+y_i) \log|x_i (x_i+y_i)| = \Omega( \log\ell)\ .\end{equation}
%\end{lem}
%\begin{proof}
%Let $\Ibig$ denote $\{i\in [\ell]: |yi| \geq |x_i|/2\}$, and let $\Ismall$ denote $[\ell]\setminus \Ibig$.
%By the monotonicity of the function $-x\log|x|$ in the range $x\in[0,e^{-1}]$, and by choosing $\ell$ large enough so that $|x_i(x_i+y_i)| \leq e^{-1}$ for $i\in \Ismall$,
%\begin{eqnarray*}
%-\sum_{i\in \Ismall} x_i(x_i+y_i)\log|x_i(x_i+y_i)|  &\geq& -\sum_{i\in \Ismall} \frac{x_i^2}2\log \frac{x_i^2}{2}\\
%&\geq & -\sum_{i\in \Ismall} \frac{x_i^2}2\log (C/\sqrt \ell) \\
%& = & \sum_{i\in \Ismall} \frac{x_i^2}4 \log \ell - \sum_{i\in \Ismall} \frac{x_i^2}2\log C\ .
%\end{eqnarray*}
%We now  estimate $\sum_{i\in \Ismall} x_i^2 = 1 -\sum_{i\in \Ibig} x_i^2 $. By definition of $\Ibig$,  $\sum_{i\in \Ibig} x_i^2 \leq  \sum_{i\in\Ibig} 4y_i^2 \leq 4\|y\|^2 \leq 1/2$.  Therefore, $\sum_{i\in \Ismall} x_i^2 \geq 1/2$ and consequently 
%\begin{equation}\label{eq:Ismall} -\sum_{i\in \Ismall} x_i(x_i+y_i)\log|x_i(x_i+y_i)|  \geq \frac 1 8 \log \ell - \frac 1 4 \log C\ .\end{equation}
%
%We now bound the contribution from $\Ibig$.  For each $i\in \Ibig$, by definition $|x_i(x_i+y_i)|\leq (C/\sqrt \ell)3|y_i|$.  Therefore, $$Z:=\sum_{i\in\Ibig} |x_i(x_y+y_i)| \leq 3(C/\sqrt \ell)\sum |y_i| \leq 3(C/\sqrt \ell)(\sqrt \ell / (30C)) = 1/10\ ,$$
%where the rightmost inequality is from the well known bound $\|y\|_1\leq \sqrt \ell \|y\|_2$.
%It is folklore (from information theory that for a vector $p_1\dots p_\ell$ of nonnegative numbers summing up to one, one has
%$-\sum p_i\log p_i \leq \log \ell$.  Therefore,
%\begin{eqnarray}
%\left | \sum_{i\in \Ibig} -x_i(x_i+y_i)\log|x_i(x_i+y_i)|\right | &\leq &
%\sum_{i\in \Ibig} -|x_i(x_i+y_i)|\log|x_i(x_i+y_i)| \nonumber \\
% &=& -\sum_{i\in\Ibig}Z\left(\frac{|x_i(x_i+y_i)|}{Z}\right)\log Z \left(\frac{|x_i(x_i+y_i)|}{Z}\right ) \nonumber \\
%&\leq &   Z\log \ell - Z\log Z\ . \nonumber \\
%&\leq & \frac 1 {10}\log \ell - \frac 1{10}\log {\frac 1{10}}\ . \label{eq:Ibig}
%\end{eqnarray}
%
%\noindent
%Combining (\ref{eq:Ismall}) with (\ref{eq:Ibig}) gives the required result.
%\end{proof}
%
%


\begin{lem}\label{lem:applemma}
%There exists a global constant $\gamma>0$ such that for all  integer $\ell$ the following holds.
%For any $C\geq 1$ there exists $\ell_0=\ell_0(C)$ such that for all $\ell>\ell_0$ the following holds:
For any large enough integer $\ell$, for any $x\in (\R^+)^{\ell}$  such that $\|x\|_1\leq 1$,  $\|x\|_\infty \leq 4\|x\|_1/\ell$
and for any $y\in \R^{\ell}$ such that $\|y\|_1 \leq C\|x\|_1$  (for some global $C$):
\begin{equation}\label{eq:applemma} -\sum_{i=1}^\ell  (x_i+y_i) \log|x_i+y_i| \geq \|x\|_1 \log \frac {\ell} {\|x\|_1} - 10\ .\end{equation}
\end{lem}
\begin{proof}
Let $\Ibig$ denote $\{i\in [\ell]: |yi| \geq x_i/2\}$, and let $\Ismall$ denote $[\ell]\setminus \Ibig$.
By the monotonicity of the function $-x\log|x|$ in the range $x\in[0,e^{-1}]$, and by choosing $\ell$ large enough so that $|x_i + y_i| \leq e^{-1}$ for $i\in \Ismall$,\footnote{By the construction, such $\ell$ can 
be chosen independently of $x,y$}
\begin{eqnarray*}
-\sum_{i\in \Ismall} (x_i+y_i)\log|x_i+y_i|  &\geq& -\sum_{i\in \Ismall} \frac{x_i}2\log \frac{x_i}{2}\\
&\geq & -\sum_{i\in \Ismall} \frac{x_i}2\log (4\|x\|_1/\ell)\\
%& = & -\|x\|_1\log(4\|x\|_1/\ell)\ .
\end{eqnarray*}
We now  estimate $\sum_{i\in \Ismall} x_i = 1 -\sum_{i\in \Ibig} x_i $. By definition of $\Ibig$,  $\sum_{i\in \Ibig} x_i \leq  \sum_{i\in\Ibig} 2y_i \leq 2C\|x\|_1 $ which is at most $\|x\|_1/2$ (assuming $C$ small enough).  Therefore, $\sum_{i\in \Ismall} x_i \geq \|x\|_1/2$.
 and consequently 
\begin{equation}\label{eq:Ismall} -\sum_{i\in \Ismall} (x_i+y_i)\log|x_i+y_i|  \geq \frac 1 4\|x\|_1 \log (\ell/(4\|x\|_1))\ .\end{equation}

We now bound the contribution from $\Ibig$.  For each $i\in \Ibig$, by our construction, $|x_i+y_i| \leq 3|y_i|$.  Therefore,
\begin{eqnarray} 
\left  |-\sum_{i\in \Ibig} (x_i+y_i)\log|x_i+y_i| \right |&\leq & -\sum_{i\in \Ibig} |x_i+y_i|\log|x_i+y_i| \nonumber\\
&\leq & -\sum_{i\in \Ibig} 3|y_i|\log|3 y_i| \nonumber\\
&\leq& 3\|y\|_1\log \ell - 3\|y\|_1\log(3\|y\|_1) \nonumber\\
&\leq& 3C\|x\|_1\log \ell - 3C\|x\|_1\log(3C\|x\|_1) \nonumber \\
&=&   3C\|x\|_1\log (\ell / (3C\|x\|_1))\ ,   \label{eq:Ibig}
 \end{eqnarray}
 where the third inequality is by well known properties of the Shannon entropy function for nonnegative vectors, and the fourth inequality assumes that $C$ is  small enough so that $3C\|x\|_1 \leq e^{-1}$ (recall
 that $-x\log x$ is monotone  increasing in $[0,e^{-1}]$).
 Choosing $C$ small enough, and combining (\ref{eq:Ismall}) with (\ref{eq:Ibig}) gives the desired result.
 \end{proof}










%
%
%\begin{lem}\label{lem:applemma}
%%There exists a global constant $\gamma>0$ such that for all  integer $\ell$ the following holds.
%%For any $C\geq 1$ there exists $\ell_0=\ell_0(C)$ such that for all $\ell>\ell_0$ the following holds:
%For any large enough integer $\ell$, for any $x\in \R^{\ell}$  such that $\|x\|_2=1$ and $\|x(i)\|_\infty \leq4/\sqrt \ell$,  and  for any $y\in \R^{\ell}$  such that $\|y\|_2 \leq C$ (for some global $C$)\ ,
%\begin{equation}\label{eq:applemma} -\sum_{i=1}^\ell  x_i(x_i+y_i) \log|x_i (x_i+y_i)| \geq \frac 1 8 \log \ell -  3\ .\end{equation}
%%where $\gamma$ is a universal constant.
%\end{lem}
%\begin{proof}
%Let $\Ibig$ denote $\{i\in [\ell]: |yi| \geq |x_i|/2\}$, and let $\Ismall$ denote $[\ell]\setminus \Ibig$.
%By the monotonicity of the function $-x\log|x|$ in the range $x\in[0,e^{-1}]$, and by choosing $\ell$ large enough so that $|x_i(x_i+y_i)| \leq e^{-1}$ for $i\in \Ismall$,
%\begin{eqnarray*}
%-\sum_{i\in \Ismall} x_i(x_i+y_i)\log|x_i(x_i+y_i)|  &\geq& -\sum_{i\in \Ismall} \frac{x_i^2}2\log \frac{x_i^2}{2}\\
%&\geq & -\sum_{i\in \Ismall} \frac{x_i^2}2\log (16/\ell) \\
%& = & \sum_{i\in \Ismall} \frac{x_i^2}4 \log\ell - \sum_{i\in \Ismall}2 x_i^2 \ .
%\end{eqnarray*}
%We now  estimate $\sum_{i\in \Ismall} x_i^2 = 1 -\sum_{i\in \Ibig} x_i^2 $. By definition of $\Ibig$,  $\sum_{i\in \Ibig} x_i^2 \leq  \sum_{i\in\Ibig} 4y_i^2 \leq 4C^2 \leq 1/2$ (assuming $C$ large enough).  Therefore, $\sum_{i\in \Ismall} x_i^2 \geq 1/2$.
% and consequently 
%\begin{equation}\label{eq:Ismall} -\sum_{i\in \Ismall} x_i(x_i+y_i)\log|x_i(x_i+y_i)|  \geq \frac 1 8 \log \ell - 2\ .\end{equation}
%
%We now bound the contribution from $\Ibig$.  For each $i\in \Ibig$, by our construction, $|x_i(x_i+y_i)|\leq 12|y_i|/\sqrt \ell$.  Therefore,
% $$Z:=\sum_{i\in\Ibig} |x_i(x_y+y_i)| \leq (12/\sqrt \ell)\sum |y_i| \leq (12/\sqrt \ell)(\sqrt \ell C)) = 12C \ ,$$
%where the rightmost inequality is from the well known bound $\|y\|_1\leq \sqrt \ell \|y\|_2$.
%It is folklore (from information theory that for a vector $p_1\dots p_\ell$ of nonnegative numbers summing up to one, one has
%$-\sum p_i\log p_i \leq \log \ell$.  Therefore, also using the fact that $|x_i(x_i+y_i)| \leq 1$ for all $i$ (assuming $\ell$ sufficiently large),
%\begin{eqnarray}
%\left | \sum_{i\in \Ibig} -x_i(x_i+y_i)\log|x_i(x_i+y_i)|\right | &\leq &
%\sum_{i\in \Ibig} -|x_i(x_i+y_i)|\log|x_i(x_i+y_i)| \nonumber \\
% &=& -\sum_{i\in\Ibig}Z\left(\frac{|x_i(x_i+y_i)|}{Z}\right)\log Z \left(\frac{|x_i(x_i+y_i)|}{Z}\right ) \nonumber \\
%&\leq &   12C\log \ell  - (12 C)\log (12 C)\ . \label{eq:Ibig} 
%\end{eqnarray}
%Taking $C= \frac 1 {192}$ gives an upper bound of $\frac 1 {16}\log \ell + 1$ on the last RHS.
%Combining (\ref{eq:Ismall}) with (\ref{eq:Ibig}) now gives the required result.
%\end{proof}
%
%
%\noindent
%Using similar techniques, the following can also be quite easily shown:
%
%\begin{lem}\label{lem:applemma2}
%For any large enough integer $\ell$, for any $x\in \R^{\ell}$  such that $\|x\|_2=1$ and $\|x(i)\|_\infty \leq4/\sqrt \ell$,  and  for any $y,z\in \R^{\ell}$  such that $\max\{\|z\|_2,\|y\|_2\} \leq C$ (for some global $C$)\ ,
%\begin{equation}\label{eq:applemma} -\sum_{i=1}^\ell  (x_i+z_i)(x_i+y_i) \log|(x_i+z_i) (x_i+y_i)| \geq \frac 1 {16} \log \ell -  6\ .\end{equation}
%%where $\gamma$ is a universal constant.
%\end{lem}



\section{Fourier $\eps$-Perturbation Can be Computed by a $(1+O(\eps))$-Well Conditioned Algorithm} \label{sec:can_compute}

For completeness, we prove the simple fact stated in the section title.  By the SVD theorem, the matrix $(\Id+\eps F)$ can be written as a product of three matrices $U\Sigma V^T$, where $U$ and $V$ are real orthogonal and $\Sigma$ is diagonal nonnegative, with the elements on the diagonal in the range $[1-\eps, 1+\eps]$.
Therefore, to compute $(\Id+\eps F)x$ we can first compute $V^T x$, using the well known fact that any real orthogonal matrix is a composition of $O(n^2)$  rotations (see Chapter 5 on \emph{Givens}  rotations in \cite{GolubV2013}).  
Continuing from there, we can compute $\Sigma V^T x$ using constant gates, one per coordinate. 
Finally, we get $U\Sigma V^T x$ be decomposing $U$ as $O(n^2)$ rotations.
Clearly this computation is $(1+O(\eps))$-well conditioned.



\end{document}

low_bound_fourier_perturbation.bbl
\begin{thebibliography}{10}

\bibitem{DBLP:journals/cjtcs/Ailon13}
Nir Ailon.
\newblock A lower bound for fourier transform computation in a linear model
  over 2x2 unitary gates using matrix entropy.
\newblock {\em Chicago J. Theor. Comput. Sci.}, 2013, 2013.

\bibitem{DBLP:conf/icalp/Ailon15}
Nir Ailon.
\newblock Tighter {F}ourier transform lower bounds.
\newblock In {\em Automata, Languages, and Programming - 42nd International
  Colloquium, {ICALP} 2015, Kyoto, Japan, July 6-10, 2015, Proceedings, Part
  {I}}, pages 14--25, 2015.

\bibitem{DBLP:journals/toct/Ailon16}
Nir Ailon.
\newblock An omega((n log n)/r) lower bound for {F}ourier transform computation
  in the r-well conditioned model.
\newblock {\em {TOCT}}, 8(1):4, 2016.

\bibitem{Avron:2010:BSL:1958627.1958633}
Haim Avron, Petar Maymounkov, and Sivan Toledo.
\newblock Blendenpik: Supercharging {L}{A}{P}{A}{C}{K}'s least-squares solver.
\newblock {\em SIAM J. Sci. Comput.}, 32(3):1217--1236, April 2010.

\bibitem{bayer1985image}
Bryce~E Bayer.
\newblock Image processing method using a collapsed {W}alsh-{H}adamard
  transform, October~22 1985.
\newblock US Patent 4,549,212.

\bibitem{Ben-Or:1983:LBA:800061.808735}
Michael Ben-Or.
\newblock Lower bounds for algebraic computation trees.
\newblock In {\em Proceedings of the Fifteenth Annual ACM Symposium on Theory
  of Computing}, STOC '83, pages 80--86, 1983.

\bibitem{berrouche2014improved}
Yaakoub Berrouche and Ra{\"\i}s~El’hadi Bekka.
\newblock Improved multiple description wavelet based image coding using
  {H}adamard {T}ransform.
\newblock {\em AEU-International Journal of Electronics and Communications},
  68(10):976--982, 2014.

\bibitem{blais2009testing}
Eric Blais.
\newblock Testing juntas nearly optimally.
\newblock In {\em Proceedings of the forty-first annual ACM symposium on Theory
  of computing}, pages 151--158. ACM, 2009.

\bibitem{blais2012testing}
Eric Blais.
\newblock {\em Testing properties of Boolean functions}.
\newblock PhD thesis, US Army, 2012.

\bibitem{boutsidis2013improved}
Christos Boutsidis and Alex Gittens.
\newblock Improved matrix algorithms via the subsampled randomized {H}adamard
  transform.
\newblock {\em SIAM Journal on Matrix Analysis and Applications},
  34(3):1301--1340, 2013.

\bibitem{1965-cooley}
James~W. Cooley and John~W. Tukey.
\newblock {An algorithm for the machine calculation of complex {F}ourier
  series}.
\newblock {\em Mathematics of Computation}, 19:297--301, 1965.

\bibitem{Cormen:2009:IAT:1614191}
Thomas~H. Cormen, Charles~E. Leiserson, Ronald~L. Rivest, and Clifford Stein.
\newblock {\em Introduction to Algorithms, Third Edition}.
\newblock The MIT Press, 3rd edition, 2009.

\bibitem{de2008brief}
Ronald De~Wolf.
\newblock A brief introduction to {F}ourier analysis on the boolean cube.
\newblock {\em Theory of Computing, Graduate Surveys}, 1:1--20, 2008.

\bibitem{decker1970hadamard}
JA~Decker.
\newblock {H}adamard--transform image scanning.
\newblock {\em Applied optics}, 9(6):1392--1395, 1970.

\bibitem{elliott2013handbook}
Douglas~F Elliott.
\newblock {\em Handbook of digital signal processing: engineering
  applications}.
\newblock Academic press, 2013.

\bibitem{10.1109/TC.1976.1674569}
B.J. Fino and V.R. Algazi.
\newblock Unified matrix treatment of the fast {W}alsh-{H}adamard transform.
\newblock {\em IEEE Transactions on Computers}, 25(11):1142--1146, 1976.

\bibitem{fischer2002testing}
Eldar Fischer, Guy Kindler, Dana Ron, Shmuel Safra, and Alex Samorodnitsky.
\newblock Testing juntas [combinatorial property testing].
\newblock In {\em Foundations of Computer Science, 2002. Proceedings. The 43rd
  Annual IEEE Symposium on}, pages 103--112. IEEE, 2002.

\bibitem{Furer:2007:FIM:1250790.1250800}
Martin F\"{u}rer.
\newblock Faster integer multiplication.
\newblock In {\em Proceedings of the Thirty-ninth Annual ACM Symposium on
  Theory of Computing}, STOC '07, pages 57--66, New York, NY, USA, 2007. ACM.

\bibitem{GolubV2013}
G.~H. Golub and C.~F. van Loan.
\newblock {\em Matrix Computations}.
\newblock Johns Hopkins University Press, Baltimore, 4th edition, 2013.

\bibitem{harvey2014even}
David Harvey, Joris Van Der~Hoeven, and Gr{\'e}goire Lecerf.
\newblock Even faster integer multiplication.
\newblock {\em arXiv preprint arXiv:1407.3360}, 2014.

\bibitem{lum2015fast}
Daniel~J Lum, Samuel~H Knarr, and John~C Howell.
\newblock Fast {H}adamard transforms for compressive sensing of joint systems:
  measurement of a 3.2 million-dimensional bi-photon probability distribution.
\newblock {\em Optics express}, 23(21):27636--27649, 2015.

\bibitem{mallat1999wavelet}
St{\'e}phane Mallat.
\newblock {\em A wavelet tour of signal processing}.
\newblock Academic press, 1999.

\bibitem{mansour1994learning}
Yishay Mansour.
\newblock Learning boolean functions via the {F}ourier transform.
\newblock In {\em Theoretical advances in neural computation and learning},
  pages 391--424. Springer, 1994.

\bibitem{martinsson2010randomized}
Gunnar Martinsson, Adrianna Gillman, Edo Liberty, Nathan Halko, Vladimir
  Rokhlin, Sijia Hao, Yoel Shkolnisky, Patrick Young, Joel Tropp, Mark Tygert,
  et~al.
\newblock Randomized methods for computing the singular value decomposition
  ({S}{V}{D}) of very large matrices.
\newblock In {\em Workshop on Algorithms for Modern Massive Data Sets, Palo
  Alto}, 2010.

\bibitem{Morgenstern:1973:NLB:321752.321761}
Jacques Morgenstern.
\newblock Note on a lower bound on the linear complexity of the fast {F}ourier
  transform.
\newblock {\em J. ACM}, 20(2):305--306, April 1973.

\bibitem{Papadimitriou:1979:OFF:322108.322118}
Christos~H. Papadimitriou.
\newblock Optimality of the fast {F}ourier transform.
\newblock {\em J. ACM}, 26(1):95--102, January 1979.

\bibitem{sarlos2006improved}
Tamas Sarlos.
\newblock Improved approximation algorithms for large matrices via random
  projections.
\newblock In {\em Foundations of Computer Science, 2006. FOCS'06. 47th Annual
  IEEE Symposium on}, pages 143--152. IEEE, 2006.

\bibitem{vetterli1984multi}
Martin Vetterli.
\newblock Multi-dimensional sub-band coding: some theory and algorithms.
\newblock {\em Signal processing}, 6(2):97--112, 1984.

\bibitem{vetterli1984simple}
Martin Vetterli, Henri~J Nussbaumer, et~al.
\newblock Simple {F}{F}{T} and {D}{C}{T} algorithms with reduced number of
  operations.
\newblock {\em Signal processing}, 6(4):267--278, 1984.

\end{thebibliography}
