It's all about groups: From Fast Fourier Transforms to QFTs

The common discrete Fourier transform

Let’s focus on the well-known discrete Fourier transform for now. As a reminder, this mathematical operation transforms a sequence \(f_1,...,f_N\) of complex numbers into another sequence of complex numbers. Sometimes the complex values are written as a function \(f(x)\) evaluated or “sampled” at regular intervals, for example using integer values \(x \in \{0,...,N-1\}\). The Fourier coefficients are then given as

where the frequencies \(k\) are the same integers as the x values. (Note that the normalisation factor is a matter of convention.) The expressions \(e^{2 \pi i \frac{k x}{N}}\) correspond to Fourier basis functions with integer-valued frequencies \(k\), and a Fourier coefficient \(\hat{f}(x_i)\) can be seen as the projection of \(f(x)\) onto the \(i\)’th basis function. The function beyond the interval \(0,..,N-1\) is thought to be “periodically continued”, which means that \(f(x_i) = f(x_i + N)\).

Let’s look at an example of such a Fourier transform:

import matplotlib.pyplot as plt
import numpy as np

N = 16

def f(x):
    """Some function on the integers 0,...,N-1."""
    x = x % N
    return 0.5*(x-4)**3

def f_hat(k):
    """Fourier coefficients of f."""
    projection = [f(x) * np.exp(2 * np.pi * 1j * k * x / N)/np.sqrt(N) for x in range(N)]
    return  np.sum(projection)

integers = np.array(range(N))
f_vec = np.array([f(x) for x in integers])
f_hat_vec = np.array([f_hat(k) for k in integers])

# plot this
fig, (ax1, ax2) = plt.subplots(2, 1)
ax1.bar(integers, np.real(f_vec), color='dimgray')  # casting to real is needed in case we perform an inverse FT
ax1.set_title(f"function f")
ax2.bar(integers + 0.05, np.imag(f_hat_vec), color='lightpink', label="imaginary part")
ax2.bar(integers, np.real(f_hat_vec), color='dimgray', label="real part")
ax2.set_title("Fourier coefficients")
plt.legend()
plt.tight_layout()
plt.show()

But why do – at least if we don’t want to incur additional headaches – the x-values have to be equidistant? Why this notion of “periodic continuation”? Why are the basis functions of exponential form? And why are the \(k\) also integers? It turns out that all of these questions have an elegant answer if we interpret the function values as a group!

More precisely, we can consider the \(x\) as elements from the set of integers \(\{0,...,N-1\}\), together with a prescription of how to combine two integers to a third from this set. Choosing “addition modulo N” for this operation (which means that \((N-1)+1 = 0\)), we get the cyclic group \(Z_N\).

This choice explains the equidistant property: integers are by nature equally spaced in \(\mathbb{R}\). It also explains the periodic continuation, as \(x = x + N\) implies \(f(x) = f(x +N)\). The Fourier basis functions are central group-theoretic objects called the characters of the group. Furthermore, the integer-valued frequencies \(k\) turn out to be elements from the so-called “dual group” \(\hat{G}\), which in this case looks exactly like the original one. We can therefore think of the \(k\) also as integers in \(\{0,...,N-1\}\) and treat Fourier space as a “mirror image” of the original space – a convenient fact that we tend to take for granted, but which does not hold for every group.

Changing the group

What happens if we exchange the cyclic group \(Z_N\) by another one? First of all, if we change to the infinite group \(\mathbb{R}\), which are just the real numbers under addition, we get the continuous Fourier transform, whose characters or basis functions look similar to the discrete ones. We could also change to a direct product of groups, such as \(Z_N \times Z_N \times ...\) or \(\mathbb{R} \times \mathbb{R} \times ...\) and get the multi-dimensional Fourier transform whose basis functions are products of characters, one for each dimension. Choosing \(N=2\) we get the group \(Z_2^n\) mentioned above, where we consider N “copies” of the binary set \(\{0,1\}\). This captures boolean logic ubiquitous in quantum algorithms. The characters for \(Z_2^n\) are of the form

where the product simply multiplies the two bits \(x_i k_i\). Each character evaluates to either -1 or 1. A character measures the parity of the subset of bits in the binary vector \(x\) that are 1 in the binary vector \(k\). A Fourier coefficient can therefore be interpreted as an average subset parity with respect to a function!

The Fourier transform with these characters is also known as the “Walsh transform”, and reads

The Fourier transform over one component group \(Z_2\) can be written as a Hadamard matrix, and the full N-dimensional version is therefore a tensor product of N Hadamards. This is exactly the reason why a quantum circuit that contains a layer of Hadamards moves into a Fourier basis as claimed above.

Let’s code up an example of a Fourier transform of a function g(x) over the group \(Z^4_2\), which transforms a function on bitstrings.

n = 4
bitstrings = [[int(c) for c in format(i, f'0{n}b')] for i in range(2**n)]

def g(x):
    """Majority function for length-n bitstrings."""
    count_ones = sum(x[i] for i in range(n))
    return count_ones/n

def character(x, k):
    """Fourier basis function for Z_2^n"""
    return np.prod([np.exp(1j * k[i] * x[i]) for i in range(n)])

def g_hat(k):
    """Fourier coefficients of f."""
    projection = [g(x) * character(x, k) for x in bitstrings]
    return  np.sum(projection)


g_vec = np.array([g(x) for x in bitstrings])
g_hat_vec = np.array([g_hat(k) for k in bitstrings])

# plot this
fig, (ax1, ax2) = plt.subplots(2, 1)
x_labels = ["".join([str(x_) for x_ in x]) for x in bitstrings]
ax1.bar(range(len(bitstrings)), np.real(g_vec), color='dimgray')  # casting to real is needed in case we perform an inverse FT
ax1.set_xticks(range(len(bitstrings)))
ax1.set_xticklabels(x_labels, rotation=45, ha='right')
ax1.set_title(f"function g")
ax2.bar(np.array(range(len(bitstrings))) + 0.05, np.imag(g_hat_vec), color='lightpink', label="imaginary part")
ax2.bar(range(len(bitstrings)), np.real(g_hat_vec), color='dimgray', label="real part")
ax2.set_xticks(range(len(bitstrings)))
ax2.set_xticklabels(x_labels, rotation=45, ha='right')
ax2.set_title("Fourier coefficients")
plt.legend()
plt.tight_layout()
plt.show()

But we can also move to any other compact group. This even includes non-Abelian groups, which are those where the composition of group elements does not commute, or \(g_1 g_2 \neq g_2 g_1\). Here, Fourier coefficients become matrix-valued objects, as the characters are replaced by matrix-valued functions called irreducible representations. (If this interests you, Parsi Diaconis’ book 2 on data analysis with group Fourier transforms is a great resource…)

The idea of Cooley-Tukey

Consider the Fourier transform of a function on \(Z_{6}\) – which for all matters and purposes one can think of as a function on the integers from 0 to 5:

The trick of the Cooley-Tukey algorithm is a change of variables such as this:

Implicitly, we are representing the set of integers \(\{0,1,2,3,4,5\}\) first as:

and then as

While doubling the amount of variables, the new variables run over fewer integers. This will allow us to “chop” the full Fourier transform into Fourier transforms over \(x_1, k_1\) only, and is the secret of the runtime reduction.

Rewriting the Fourier transform in the new variables reads

Besides some reordering, there is one slightly non-trivial identity we used above:

Essentially, the change in variables turns the Fourier basis function over period 6 into a Fourier basis function with period 2, which makes the dependency on \(k_2\) disappear. This effectively turns the Fourier transform into a sum of “smaller” Fourier transforms over \(3x_1\), namely \(\{0, 3\}\), \(\{1, 4\}\), and \(\{2, 5\}\). These are combined with an appropriate “twiddle factor” \(e^{\frac{2 \pi i}{6} (2k_2+k_1) x_2}\) that weighs the building blocks.

We explicitly implement this with the function used at the very beginning:

N = 6

def f2(x1, x2):
    """Function f from before, but with input split
       into two variables, and running from 0,...5."""
    x = 3*x1 + x2
    x = x % N
    return 0.5*(x-4)**3

def f_subgroup(x2, k1):
    """Fourier transform over a smaller group or coset."""
    projection = [f2(x1, x2) * np.exp(2*np.pi*1j * k1*x1/2)/np.sqrt(2) for x1 in range(2)]
    return  np.sum(projection)

def f_hat_fft(k1, k2):
    """Putting the smaller transforms together."""
    res = 0
    for x2 in range(3):
        twiddle_factor = np.exp(2 * np.pi * 1j * (2*k2+k1) * x2 / 6) / np.sqrt(3)
        res += twiddle_factor * f_subgroup(x2, k1)
    return res

f_hat_vec_fft = np.array([f_hat(2*k2+k1) for k1 in range(2) for k2 in range(3)])

Let’s compare the result to the DFT from before.

def f(x):
    """Function f from before, but restricted to the integers 0,...,5."""
    x = x % N
    return 0.5*(x-4)**3

def f_hat(k):
    """Fourier coefficients of f."""
    projection = [f(x) * np.exp(2 * np.pi * 1j * k * x / N)/np.sqrt(N) for x in range(N)]
    return  np.sum(projection)

f_hat_vec = np.array([f_hat_fft(k1, k2) for k1 in range(2) for k2 in range(3)])
print("FFT and DFT coincide:", np.allclose(f_hat_vec, f_hat_vec_fft))

FFT and DFT coincide: True

Quantifying the runtime gains

The FFT algorithm devises to compute \(\hat{f}(k) = \hat{f}(2k_2+k_1)\) in two steps:

1. First, compute the terms in the inner sum, \(\tilde{f}(x_2, k_1) = \sum_{x_1=0}^1 e^{\frac{2 \pi i}{2} x_1 k_1 } f(3 x_1 + x_2)\) as done in f_subgroup. There are 9 such terms (as we have 3 options to pick \(k_1\) and 3 options to pick \(x_2\)). Each term requires summing over 2 expressions (as \(x_1\) runs over 2 values). Overall, there are 18 terms involved.

2. Second, combine the 9 terms using the outer sum \(\hat{f}(k_2, k_1) = \sum_{x_2=0}^2 e^{\frac{2 \pi i}{6} (2k_2+k_1) x_2} \tilde{f}(k_1, x_2)\). This computes \(2 \cdot 3=6\) Fourier coefficients, and each has to sum over 3 terms (as \(x_2\) runs over 3 values). Overall, there are 12 terms involved.

Together, the FFT algorithm includes 12 + 18 = 30 terms. The naive sum, instead, would compute the 6 Fourier coefficients with a sum over 6 terms each, using \(6 \cdot 6=36\) terms altogether.

Of course, in this “mini” example the savings are not very impressive. The overall scaling however reduces from \(O(N^2)\) to \(O(N \mathrm{log}N)\), which is a quadratic speedup, and made it possible for the FFT to become one of the most important algorithms in science and engineering.

A group-theoretic interpretation

As hinted at before, the change of variables that was the core idea behind the Cooley-Tukey implementation of the Fast Fourier Transform is in fact a decomposition of the original group – here \(Z_{6}\) – into a subgroup with the elements \(\{0,3\}\), and its cosets or “shifted copies”, \(\{1,4\}\) and \(\{2,5\}\) (see 3, 4). The subgroup is “isomorphic to” \(Z_{2}\), which loosely speaking means that it “behaves like” \(Z_{2}\).

The new variable \(x_2\) selects between cosets, while \(x_1\) picks within a coset. For example, for \(5=3 \cdot 1 + 2\) we have that \(x_2=2\) selects the third coset of the subgroup, while \(x_1=1\) selects the second element within the coset. It is a fundamental result in group theory that “shifted” copies of a subgroup tile the entire group, which is why the new coordinates can always be used to express any element \(x \in G\).

The splitting of the variable \(k\) is related to the concept of “restricting a character” to the subgroup. Essentially, it allows us to turn the characters, or Fourier basis functions \(e^{\frac{2 \pi i}{6} x k }\) related to the original group, into characters of the subgroup, \(e^{\frac{2 \pi i}{2} x_1 k_1 }\) by changing the period. Surprisingly, such a trick generalises to much more complicated groups, including those where there is no “cyclic” notion of ordered integers.

Conclusion

Congratulations if you made it to the end, this was a lot of material to cover! Let’s summarise everything we learned as three important insights:

1. The Fourier transform can be defined with respect to functions over groups. The usual discrete Fourier transform, for example, uses the cyclic group of integers.

2. The Fast Fourier Transform exploits that for some groups, a Fourier transform can be decomposed into “smaller” Fourier transforms that consider subgroups (and their copies). This provides a speedup to nearly-linear runtime in the size of the group.

3. The Quantum Fourier Transform calculates these smaller blocks in “quantum parallel”, thereby realising a runtime that is logarithmic in the size of the group.

Group structure and Fourier transforms are important building blocks of quantum algorithms, and a fascinating reason for why quantum computers might offer unique ways of information processing. Of course, the caveat is that quantum implementations of Fourier transforms cannot compute Fourier coefficients directly. After a QFT, the Fourier coefficients are the amplitudes of a quantum state. All we can do is sample from the final state, or manipulate it in interesting ways – creating a wonderful playground for quantum algorithm design.