What is the time complexity of fermionic Fourier transform?

Question

Suppose $N = 2^L$ and we are interested in performing the following transformation a $\mapsto$ a_hat on arrays of $N$ complex numbers in time $O(N \log(N))$ or similar (as opposite to $O(N^2)$).

Preliminary definitions

The following construction is motivated by physics. Yet, we write it in a self-contained form to make sure no background in physics is needed to understand the question. Consider an $N$-dimensional vector space $\mathcal{H} = \mathbb{C}^N$ with $N = 2^L$ basis elements $\left|z\right>$ ($z=0,\dots,N-1$). Array a is then interpreted as a vector in $\mathcal{H}$ equal to $\sum_z a_z \left|z\right>$. For each $j=0,\dots,L-1$ we introduce maps $c^\dagger_j \colon \mathcal{H} \to \mathcal{H}$ given by $$c^\dagger_j \left|z\right> = \begin{cases} 0&\text{ if }z_j = 1,\\ (-1)^{\text{popcnt}(\lfloor z / 2^j\rfloor)} \left|z\oplus 2^j\right>&\text{ if } z_j = 0. \end{cases} \tag{1}$$ Here $\oplus$ means bitwise xor (commonly denoted by ^ in programming languages) and $z_j$ is $j$-th bit of $z$ (i.e. (z >> j) & 1). Expression $\text{popcnt}(\lfloor z / 2^j\rfloor)$ counts the number of bits set to $1$ in $z$ to the left of position $j$. The maps $c^\dagger_j$ are called "creation operators" in physics. One can then check that the composition of these maps changes sign when they change the order: $$c_j^\dagger c_l^\dagger = -c_l^\dagger c_j^\dagger. \tag{2}$$ Also, $c_j^\dagger c_j^\dagger = 0$.

With this notation we can write $$ \left|z\right> = \left(\prod_{j=L-1,\dots,0} (c^{\dagger}_{j})^{z_j}\right)\left|0\right>, \tag{3} $$ Where the product is just the composition of the maps and $(c^{\dagger}_{j})^{z_j}$ is $c^{\dagger}_{j}$ if $z_j=1$ and identity map if $z_j=0$. The ordering is such that the maps with lowest $j$ are applied first.

Definition of fermionic Fourier transform For the purpose of this question the fermionic Fourier transform is defined as follows. Take array a as input. It represents a vector $$\vec{a} = \sum_z a_z \left|z\right> = \sum_z a_z \left(\prod_{j=L-1,\dots,0} (c^{\dagger}_{j})^{z_j}\right)\left|0\right>. \tag{4}$$ In the right hand side of this expression replace every $c^{\dagger}_j$ by $\frac{1}{\sqrt{L}} \sum_{l=0}^{L-1} e^{-2\pi i jl / L} c^{\dagger}_l$ to obtain the vector $\vec{\hat{a}}$. Output elements of $\vec{\hat{a}}$ in the basis $\{\left|z\right>\}_{z=0}^{N}$.

Example

$L=4$, $N=16$, $a = [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]$ representing the vector $$\vec{a} = \left|0011_2\right> = c^\dagger_1 c^\dagger_0 \left|0\right>$$ (we write $z$ in binary notation). Then $$ \vec{\hat{a}} = \frac{1}{4} (c_0^\dagger - ic_1^\dagger - c_2^\dagger + ic_3^\dagger)(c_0^\dagger + c_1^\dagger + c_2^\dagger + c_3^\dagger)\left|0\right> = \frac{1}{4}\Bigl(-(1+i)\left|0011\right>, -2\left|0101\right>, (-1+i)\left|0110\right>, (-1+i)\left|1001\right>, 2i \left|1010\right>, (1+i)\left|1100\right>\Bigr).$$ Here in the second equality we opened the brackets, sorted $c^\dagger_l$ in decreasing order (adding $-$ signs as required by (2)), and used (3) to write the terms as basis elements.

Thus, the correct implementation of the fermionic Fourier transform would return a_hat = [0, 0, 0, -(1+i)/4, 0, -1/2, (-1+i)/4, 0, 0, (-1+i)/4, i/2, 0, (1+i)/2, 0, 0, 0].

Question

Since this is a linear operation on $N$-dimensional vectors, there is a trivial algorithm (for fixed $L$): precompute the matrix describing this linear operation, then do matrix-vector multiplication in $O(N^2)$ time. This can be improved to $O(N^2/\sqrt{\log(N)})$ by handling portions of a with different popcnt(z) separately. However, similarly to the standard FFT, one may expect that this can be done faster (possibly in $O(N \log(N))$ or, at least, $O(N \text{Poly}(\log(N)))$ time). Does such faster algorithm (for a classical computer) exist?

Answer 1

It is at most $O(2^L L^2)$

openfermion.circuits.ffft implements the fermionic Fourier transform on a gate quantum computer. One can simulate this algorithm on a classical computer at the multiplicative cost of $2^L$, which gives the overall time complexity of $O(2^L m^*(L))$, where $m^*(L)$ is the minimal size of decomposition of the classical DFT matrix $F_{jl} = \frac{1}{\sqrt{L}} e^{-2\pi jl / L}$ into 2x2 unitaries (i.e. unitary matrices which match with the identity matrix except for the intersection of two columns with rows with the same indices).

The algorithm (for classical computer) works as follows. First, notice that Eq. (4) in the question changes all creation operators $c^\dagger_j$. What if we changed just 2 of them (by applying a unitary matrix $U$ to the vector containing these two creation operators), say $c^\dagger_{j_0}$ and $c^\dagger_{j_1}$? In this case $a_z$ would remain the same if $z_{j_0} = z_{j_1} = 0$ because these terms do not involve any of the two affected creation operators. If $z_{j_0} = z_{j_1} = 1$, then $a_z$ would be multiplied by $\det(U)$. The remaining coefficients $a_z$ can be split into pairs with exactly one of $j_0$ and $j_1$-th bits set. For such pair $(z,\tilde{z})$ the coefficients $(a_z, a_{\tilde{z}})$ are multiplied by the matrix $U$ and the sign equal to $(-1)$ to the power of bits set in $z$ and $\tilde z$ between $j_0$ and $j_1$. Thus, such $2\times 2$ unitary can be applied in $O(2^L)$ time. What remains is to decompose the matrix $F$ describing the classical Fourier transform into the product of such $2\times 2$ unitaries (and, optionally, a permutation matrix which would correspond to permuting bits of $z$ and adjusting the signs of $a_z$ accordingly).

One can follow Cooley-Tukey algorithm to show that if $L=L_1\cdots L_k$ is composite then $m^*(L) \leq \sum_{j=1}^{k} m^*(L_j) L / L_j$, which gives $m^*(L) = L\log_2(L)/2$ when $L$ itself is a power of 2 (or $m^*(L) = O(L\log_2(L))$ when $L$ is a product of small primes). For prime $L$ I only can show that $m^*(L) \leq (L-1)(L+3)/4$, which gives the overall complexity of $O(L^2 2^L)$. I've posted a question on whether a decomposition into $O(L\log(L))$ 2x2 matrices is always possible (even for prime $L$) here.