Fermionic fast Fourier transform

This page explains the fermionic fast Fourier transform, an efficient quantum circuit for the discrete Fourier transform on fermionic modes.

The discrete Fourier transform on fermionic creation operators is defined as

\[a^\dagger_k \mapsto \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} e^{-i 2\pi k n / N} a^\dagger_n\]

This is an orbital rotation whose matrix is the unitary DFT matrix \(\mathbf{F}_N\) with entries \((\mathbf{F}_N)_{kn} = \omega_N^{kn} / \sqrt{N}\), where \(\omega_N = e^{-i 2\pi / N}\). Rather than decomposing \(\mathbf{F}_N\) directly into Givens rotations, the fermionic fast Fourier transform (FFFT) uses the Cooley-Tukey algorithm to break it into DFTs of prime size.

Mixed-radix Cooley-Tukey algorithm

The circuit follows the general factorization described in the Variations section of the Cooley-Tukey FFT algorithm Wikipedia page. First, \(N\) is factored into primes using trial division. Writing \(N = N_1 N_2\), where \(N_2\) is the smallest prime factor, and splitting the input index as \(n = N_1 n_2 + n_1\) and the output index as \(k = N_2 k_1 + k_2\), the DFT becomes

\[X_{N_2 k_1 + k_2} = \sum_{n_1=0}^{N_1-1} \left[ e^{-\frac{2\pi i}{N_1N_2} n_1 k_2 } \right] \left( \sum_{n_2=0}^{N_2-1} x_{N_1 n_2 + n_1} e^{-\frac{2\pi i}{N_2} n_2 k_2 } \right) e^{-\frac{2\pi i}{N_1} n_1 k_1 }\]

where each inner sum is a DFT of size \(N_2\), each outer sum is a DFT of size \(N_1\).

The circuit follows this formula:

1. Perform \(N_1\) DFTs of size \(N_2\), one for each \(n_1\) (the inner sum). Since \(N_2\) is prime, each is implemented directly as an orbital rotation.

2. Multiply by the twiddle factors \(e^{-\frac{2\pi i}{N_1N_2} n_1 k_2 }\). Each twiddle factor multiplies a single mode, so it is implemented with a single-qubit phase gate.

3. Perform \(N_2\) DFTs of size \(N_1\), one for each \(k_2\) (the outer sum). These are computed by applying the same procedure recursively.

The recursion ends when the size is prime, in which case the DFT is implemented directly as an orbital rotation.

Under the Jordan-Wigner transformation, a gate implementing an orbital rotation must act on a contiguous block of qubits, because a hopping term between two modes carries a string of \(Z\) operators on all the qubits in between. The DFTs in steps 1 and 3 act on strided subsets of modes, so the circuit inserts fermionic mode permutations (themselves orbital rotations by permutation matrices) that bring each subset onto adjacent qubits before the DFTs are applied.

Quantum circuit implementation

ffsim implements the FFFT under the Jordan-Wigner transformation with the FermionicFFTSpinlessJW gate for spinless fermions and the FermionicFFTJW gate for spinful fermions.

The following code cell constructs and displays a Qiskit circuit for implementing the fermionic fast Fourier transform on \(N = 6\) spatial orbitals, so \(N_2 = 2\) and \(N_1 = 3\). In the decomposed circuit you can see the mode permutations, the three size-2 DFTs of step 1 and the two size-3 DFTs of step 3 (as orbital rotations), and the twiddle-factor phase gates of step 2.

from qiskit.circuit import QuantumCircuit, QuantumRegister

import ffsim

norb = 6

# Spinless FFFT circuit (acts on norb qubits)
qubits = QuantumRegister(norb, name="q")
circuit = QuantumCircuit(qubits)
circuit.append(ffsim.qiskit.FermionicFFTSpinlessJW(norb), qubits)

# Decompose and draw
circuit.decompose().draw("mpl", scale=0.6)

We can check that the circuit implements the DFT orbital rotation by comparing it against apply_orbital_rotation with the DFT matrix on a random state.

import numpy as np
import scipy.linalg
from qiskit.quantum_info import Statevector

nelec = 3  # spinless, so nelec is an integer
vec = ffsim.random.random_state_vector(ffsim.dim(norb, nelec), seed=1234)

# Apply the FFFT circuit to the state
qiskit_vec = ffsim.qiskit.ffsim_vec_to_qiskit_vec(vec, norb=norb, nelec=nelec)
result = Statevector(qiskit_vec).evolve(circuit)
result = ffsim.qiskit.qiskit_vec_to_ffsim_vec(np.array(result), norb=norb, nelec=nelec)

# Apply the DFT matrix as an orbital rotation
dft_mat = scipy.linalg.dft(norb, scale="sqrtn")
expected = ffsim.apply_orbital_rotation(vec, dft_mat, norb=norb, nelec=nelec)

np.testing.assert_allclose(result, expected)

For spinful fermions, FermionicFFTJW acts on \(2N\) qubits and applies the FFFT independently to the spin \(\alpha\) orbitals (the first \(N\) qubits) and the spin \(\beta\) orbitals (the last \(N\) qubits). The following code cell applies it to the Hartree-Fock state.

norb = 4
nelec = (2, 2)  # 2 alpha electrons, 2 beta electrons

qubits = QuantumRegister(2 * norb, name="q")
circuit = QuantumCircuit(qubits)
circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)
circuit.append(ffsim.qiskit.FermionicFFTJW(norb), qubits)

circuit.decompose().draw("mpl", scale=0.6)

Cost compared with a direct orbital rotation

The DFT could also be implemented as a single orbital rotation by passing the DFT matrix to OrbitalRotationSpinlessJW, which decomposes it directly into Givens rotations. The following code cell compiles both circuits to a basis of CZ, \(R_Z\), \(\sqrt{X}\), and \(X\) gates. For each one it reports the number of CZ gates, the number of \(R_Z\) rotations whose angle is not a multiple of \(\pi/2\) (that is, the non-Clifford rotations), and how many distinct angles those rotations use.

import math

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

basis_gates = ["cz", "id", "rz", "sx", "x"]
pass_manager = generate_preset_pass_manager(
    optimization_level=1, basis_gates=basis_gates
)


def cost(gate) -> tuple[int, int, int]:
    """Return CZ count, non-Clifford Rz count, and number of distinct angles."""
    circuit = QuantumCircuit(gate.num_qubits)
    circuit.append(gate, range(gate.num_qubits))
    transpiled = pass_manager.run(circuit)
    angles = []
    for instruction in transpiled.data:
        if instruction.operation.name == "rz":
            angle = float(instruction.operation.params[0]) % (2 * math.pi)
            remainder = angle % (math.pi / 2)
            if min(remainder, math.pi / 2 - remainder) > 1e-8:
                angles.append(round(angle, 8))
    return transpiled.count_ops().get("cz", 0), len(angles), len(set(angles))


print(f"{'':>3}{'CZ gates':^18}{'non-Clifford Rz':^18}{'distinct angles':^18}")
print(f"{'N':>3}" + f"{'FFFT':>9}{'direct':>9}" * 3)
# List powers of 2 first, then sizes with an odd prime factor
for n in [4, 8, 16, 6, 12]:
    ffft = cost(ffsim.qiskit.FermionicFFTSpinlessJW(n))
    direct = cost(
        ffsim.qiskit.OrbitalRotationSpinlessJW(n, scipy.linalg.dft(n, scale="sqrtn"))
    )
    print(f"{n:>3}" + "".join(f"{f:>9}{d:>9}" for f, d in zip(ffft, direct)))

        CZ gates      non-Clifford Rz   distinct angles
  N     FFFT   direct     FFFT   direct     FFFT   direct
  4       14       12        8       18        1        6
  8       72       56       26      105        3       49
 16      328      240       74      465        7      175
  6       36       30       32       56        6       28
 12      174      132       80      253        7      104

The FFFT needs far fewer non-Clifford rotations than the direct orbital rotation, and those rotations use only a handful of distinct angles. This is most pronounced when \(N\) is a power of 2. Then every Givens rotation is either a mode permutation (angle \(\pi/2\)) or a DFT of size 2 (angle \(\pi/4\)), and the only other angles come from twiddle factors, which are multiples of \(2\pi/N\). The direct orbital rotation, by contrast, needs a generic angle for almost every rotation. On fault-tolerant hardware, where each arbitrary-angle rotation must be synthesized from many expensive non-Clifford gates, this makes the FFFT much cheaper.

The price is more CZ gates, since the mode permutations add two-qubit gates that the direct decomposition doesn’t need. When CZ count is the main cost, as on current noisy hardware, the direct orbital rotation is the better choice.