<h1>It is at most <span class="math-container">$O(2^L L^2)$</span></h1>
<p><a href="https://quantumai.google/reference/python/openfermion/circuits/ffft" rel="nofollow noreferrer">openfermion.circuits.ffft</a> implements the fermionic Fourier transform on a gate quantum computer. One can simulate this algorithm on a classical computer at the multiplicative cost of <span class="math-container">$2^L$</span>, which gives the overall time complexity of <span class="math-container">$O(2^L m^*(L))$</span>, where <span class="math-container">$m^*(L)$</span> is the minimal size of decomposition of the classical DFT matrix <span class="math-container">$F_{jl} = \frac{1}{\sqrt{L}} e^{-2\pi jl / L}$</span> 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).</p>
<p>The algorithm (for classical computer) works as follows. First, notice that Eq. (4) in the question changes all creation operators <span class="math-container">$c^\dagger_j$</span>. What if we changed just 2 of them (by applying a unitary matrix <span class="math-container">$U$</span> to the vector containing these two creation operators), say <span class="math-container">$c^\dagger_{j_0}$</span> and <span class="math-container">$c^\dagger_{j_1}$</span>? In this case <span class="math-container">$a_z$</span> would remain the same if <span class="math-container">$z_{j_0} = z_{j_1} = 0$</span> because these terms do not involve any of the two affected creation operators. If <span class="math-container">$z_{j_0} = z_{j_1} = 1$</span>, then <span class="math-container">$a_z$</span> would be multiplied by <span class="math-container">$\det(U)$</span>. The remaining coefficients <span class="math-container">$a_z$</span> can be split into pairs with exactly one of <span class="math-container">$j_0$</span> and <span class="math-container">$j_1$</span>-th bits set. For such pair <span class="math-container">$(z,\tilde{z})$</span> the coefficients <span class="math-container">$(a_z, a_{\tilde{z}})$</span> are multiplied by the matrix <span class="math-container">$U$</span> and the sign equal to <span class="math-container">$(-1)$</span> to the power of bits set in <span class="math-container">$z$</span> and <span class="math-container">$\tilde z$</span> between <span class="math-container">$j_0$</span> and <span class="math-container">$j_1$</span>. Thus, such <span class="math-container">$2\times 2$</span> unitary can be applied in <span class="math-container">$O(2^L)$</span> time. What remains is to decompose the matrix <span class="math-container">$F$</span> describing the classical Fourier transform into the product of such <span class="math-container">$2\times 2$</span> unitaries (and, optionally, a permutation matrix which would correspond to permuting bits of <span class="math-container">$z$</span> and adjusting the signs of <span class="math-container">$a_z$</span> accordingly).</p>
<p>One can follow Cooley-Tukey algorithm to show that if <span class="math-container">$L=L_1\cdots L_k$</span> is composite then <span class="math-container">$m^*(L) \leq \sum_{j=1}^{k} m^*(L_j) L / L_j$</span>, which gives <span class="math-container">$m^*(L) = L\log_2(L)/2$</span> when <span class="math-container">$L$</span> itself is a power of 2 (or <span class="math-container">$m^*(L) = O(L\log_2(L))$</span> when <span class="math-container">$L$</span> is a product of small primes). For prime <span class="math-container">$L$</span> I only can show that <span class="math-container">$m^*(L) \leq (L-1)(L+3)/4$</span>, which gives the overall complexity of <span class="math-container">$O(L^2 2^L)$</span>. I've posted a question on whether a decomposition into <span class="math-container">$O(L\log(L))$</span> 2x2 matrices is always possible (even for prime <span class="math-container">$L$</span>) <a href="https://cs.stackexchange.com/questions/162025/can-fft-for-prime-n-be-implemented-as-a-composition-of-on-logn-2-times">here</a>.</p>

CC BY-SA 4.0
2023-09-14T00:48:01.787Z
2023-09-14T00:48:01.787Z
fiktor
<p>Suppose <span class="math-container">$N = 2^L$</span> and we are interested in performing the following transformation <code>a</code> <span class="math-container">$\mapsto$</span> <code>a_hat</code> on arrays of <span class="math-container">$N$</span> complex numbers in time <span class="math-container">$O(N \log(N))$</span> or similar (as opposite to <span class="math-container">$O(N^2)$</span>).</p>
<p><strong>Preliminary definitions</strong></p>
<p>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 <span class="math-container">$N$</span>-dimensional vector space <span class="math-container">$\mathcal{H} = \mathbb{C}^N$</span> with <span class="math-container">$N = 2^L$</span> basis elements <span class="math-container">$\left|z\right&gt;$</span> (<span class="math-container">$z=0,\dots,N-1$</span>). Array <code>a</code> is then interpreted as a vector in <span class="math-container">$\mathcal{H}$</span> equal to <span class="math-container">$\sum_z a_z \left|z\right&gt;$</span>. For each <span class="math-container">$j=0,\dots,L-1$</span> we introduce maps <span class="math-container">$c^\dagger_j \colon \mathcal{H} \to \mathcal{H}$</span> given by
<span class="math-container">$$c^\dagger_j \left|z\right&gt; = \begin{cases}
0&amp;\text{ if }z_j = 1,\\
(-1)^{\text{popcnt}(\lfloor z / 2^j\rfloor)} \left|z\oplus 2^j\right&gt;&amp;\text{ if } z_j = 0.
\end{cases} \tag{1}$$</span>
Here <span class="math-container">$\oplus$</span> means bitwise xor (commonly denoted by <code>^</code> in programming languages) and <span class="math-container">$z_j$</span> is <span class="math-container">$j$</span>-th bit of <span class="math-container">$z$</span> (i.e. <code>(z &gt;&gt; j) &amp; 1</code>).
Expression <span class="math-container">$\text{popcnt}(\lfloor z / 2^j\rfloor)$</span> counts the number of bits set to <span class="math-container">$1$</span> in <span class="math-container">$z$</span> to the left of position <span class="math-container">$j$</span>.
The maps <span class="math-container">$c^\dagger_j$</span> are called &quot;creation operators&quot; in physics. One can then check that the composition of these maps changes sign when they change the order:
<span class="math-container">$$c_j^\dagger c_l^\dagger = -c_l^\dagger c_j^\dagger. \tag{2}$$</span>
Also, <span class="math-container">$c_j^\dagger c_j^\dagger = 0$</span>.</p>
<p>With this notation we can write
<span class="math-container">$$
\left|z\right&gt; = \left(\prod_{j=L-1,\dots,0} (c^{\dagger}_{j})^{z_j}\right)\left|0\right&gt;,
\tag{3}
$$</span>
Where the product is just the composition of the maps and <span class="math-container">$(c^{\dagger}_{j})^{z_j}$</span> is <span class="math-container">$c^{\dagger}_{j}$</span> if <span class="math-container">$z_j=1$</span> and identity map if <span class="math-container">$z_j=0$</span>. The ordering is such that the maps with lowest <span class="math-container">$j$</span> are applied first.</p>
<p><strong>Definition of fermionic Fourier transform</strong>
For the purpose of this question the fermionic Fourier transform is defined as follows. Take array <code>a</code> as input. It represents a vector
<span class="math-container">$$\vec{a} = \sum_z a_z \left|z\right&gt; = \sum_z a_z \left(\prod_{j=L-1,\dots,0} (c^{\dagger}_{j})^{z_j}\right)\left|0\right&gt;. \tag{4}$$</span>
In the right hand side of this expression replace every <span class="math-container">$c^{\dagger}_j$</span> by <span class="math-container">$\frac{1}{\sqrt{L}} \sum_{l=0}^{L-1} e^{-2\pi i jl / L} c^{\dagger}_l$</span> to obtain the vector <span class="math-container">$\vec{\hat{a}}$</span>. Output elements of <span class="math-container">$\vec{\hat{a}}$</span> in the basis <span class="math-container">$\{\left|z\right&gt;\}_{z=0}^{N}$</span>.</p>
<p><strong>Example</strong></p>
<p><span class="math-container">$L=4$</span>, <span class="math-container">$N=16$</span>, <span class="math-container">$a = [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]$</span> representing the vector
<span class="math-container">$$\vec{a} = \left|0011_2\right&gt; = c^\dagger_1 c^\dagger_0 \left|0\right&gt;$$</span>
(we write <span class="math-container">$z$</span> in binary notation). Then
<span class="math-container">$$ \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&gt; = \frac{1}{4}\Bigl(-(1+i)\left|0011\right&gt;, -2\left|0101\right&gt;, (-1+i)\left|0110\right&gt;, (-1+i)\left|1001\right&gt;, 2i \left|1010\right&gt;, (1+i)\left|1100\right&gt;\Bigr).$$</span>
Here in the second equality we opened the brackets, sorted <span class="math-container">$c^\dagger_l$</span> in decreasing order (adding <span class="math-container">$-$</span> signs as required by (2)), and used (3) to write the terms as basis elements.</p>
<p>Thus, the correct implementation of the fermionic Fourier transform would return
<code>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]</code>.</p>
<p><strong>Question</strong></p>
<p>Since this is a linear operation on <span class="math-container">$N$</span>-dimensional vectors, there is a trivial algorithm (for fixed <span class="math-container">$L$</span>): precompute the matrix describing this linear operation, then do matrix-vector multiplication in <span class="math-container">$O(N^2)$</span> time. This can be improved to <span class="math-container">$O(N^2/\sqrt{\log(N)})$</span> by handling portions of <code>a</code> with different <code>popcnt(z)</code> separately. However, similarly to the standard FFT, one may expect that this can be done faster (possibly in <span class="math-container">$O(N \log(N))$</span> or, at least, <span class="math-container">$O(N \text{Poly}(\log(N)))$</span> time). Does such faster algorithm (for a classical computer) exist?</p>

<time-complexity><fourier-analysis>
What is the time complexity of fermionic Fourier transform?
fiktor
# What is the time complexity of fermionic Fourier transform?

(Asked by: username_0 on 2023-04-10)

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?

**Tags:** time-complexity, fourier-analysis

---

## 1 answers

### Answer by username_0 on 2023-09-14. Score: 0

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

[openfermion.circuits.ffft](https://quantumai.google/reference/python/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](https://cs.stackexchange.com/questions/162025/can-fft-for-prime-n-be-implemented-as-a-composition-of-on-logn-2-times).


---

