Here's the note:

Fourier Transform and Fast Multiplication

When to Reach for Fourier

Here we revise when Fourier is the right tool, before how it works. The thing you have to remember: the test is about the system, not the signal.

The LTI Trigger

Fourier diagonalises every operator that is linear and shift-invariant (LTI). Complex exponentials are its eigenfunctions.

• Linear: $T(\alpha x+\beta y)=\alpha T(x)+\beta T(y)$.
• Shift-invariant: $T$ commutes with the shift $S_{\tau }$, i.e. $T(S_{\tau }x)=S_{\tau }T(x)$.
• Consequence: $T(e^{i\omega t})=H(\omega )e^{i\omega t}$. A sine goes in and the same sine comes out, only scaled and phase-shifted.

Why: the eigenvectors of the shift are exponentials, and an operator that commutes with the shift shares its eigenvectors (see Simultaneous Diagonalisation). Any LTI operator is a Convolution, and in the Fourier basis convolution becomes pointwise multiplication.

Calibration Drill

Case	Fourier?	Why
Un-blurring a photo	Yes	Blur = convolution → divide spectra; regularise where the kernel's spectrum ≈ 0 (Wiener Filter)
Singing in the shower	Yes	LTI room; hard tiles boost resonant modes and add reverb
Tides	Yes	Periodic astronomical forcing + approximately linear ocean → sum of known-period constituents
Stock prices	No	Non-stationary, no stable periods
Big-integer multiplication	Yes	Long multiplication is a convolution of digits
Heat on a ring	Yes	Mode $k$ decays as $e^{-\kappa k^{2}t}$; sharp bumps die first
A single clap	Bad lens for the clap	…but an excellent probe: impulse ≈ all frequencies → gives the room's response
JPEG	Yes	DCT; high-frequency coefficients stored crudely
Guitar vs drum	Yes, both	String modes are integer multiples (pitch); drumhead modes are not (≈1, 1.59, 2.14…)
Counting a crowd	No	Localised objects are smeared over all frequencies (except regular grids)

Anti-Triggers

• Localised in time or space (transients, edges) → Short-Time Fourier Transform, Wavelets. Uncertainty principle: narrow in time means broad in frequency.
• Nonlinear systems create new frequencies (harmonics), so no clean per-frequency description exists.
• Non-stationary rules (the system changes over time) break shift-invariance.

Convolution in Disguise

Multiplication Is Convolution

Write digits units-first, so position = power of 10: $a={\sum }_i{a}_{i}10^i$. The core fact is that exponents add when you multiply.

$$\begin{array}{rl}a\cdot b& =\sum _i\sum _j{a}_i{b}_{j}10^{i+j}=\sum _k{c}_{k}10^k,\\ c_k& =\sum _{i+j=k}{a}_i{b}_j=\sum _i{a}_i{b}_{k-i}\text{(convolution)}\end{array}$$

Example: $123\times 45$, so $a=[3,2,1]$ and $b=[5,4]$:

	$j=0$ (5)	$j=1$ (4)
$i=0$ (3)	15 → $k=0$	12 → $k=1$
$i=1$ (2)	10 → $k=1$	8 → $k=2$
$i=2$ (1)	5 → $k=2$	4 → $k=3$

Summing the anti-diagonals gives $c=[15,22,13,4]$. Carrying then gives $5535$.

Schoolbook multiplication without carrying shows exactly the convolution:

          1    2    3
      ×        4    5
          5   10   15
     4    8   12
     4   13   22   15   ← c3 c2 c1 c0

• Carrying is separate: it uses floor and mod, so it is nonlinear. Do it once, after the convolution.
• Shifted copies: $a\cdot b={\sum }_j{b}_j{S}^{j}a$. This is a weighted sum of shifted copies of $a$, structurally identical to a blur kernel.

The Flip and the ML Convolution

Sliding a kernel over "adjacent cells" (a CNN layer) is this operation. Adjacent cells = adjacent digit positions.

position:   0  3  2  1  0       (123 units-first, zero-padded)
k=0:       [4  5]           → 15
k=1:          [4  5]        → 22
k=2:             [4  5]     → 13
k=3:                [4  5]  → 4

	Math convolution	ML "convolution"
Formula	${\sum }_i{a}_i{b}_{k-i}$	${\sum }_i{a}_{k+i}{w}_i$
Kernel	Flipped	Not flipped (cross-correlation)
Why it doesn't matter	—	Kernel is learned, so the network learns the flip

Circular Convolution and Padding

The DFT indexes modulo $N$, so it computes circular convolution: high positions wrap onto low ones. With $N=3$, $c_3$ wraps into $c_0$, giving $15+4=19$ ✗.

Pad to $N\ge \text{len}(a)+\text{len}(b)-1$.

Same phenomenon as zero-padding vs circular padding in CNNs.

Other Disguised Convolutions

• Polynomial multiplication (coefficients).
• PMF of a sum of independent random variables → Characteristic Function.
• Moving averages, discretised diffusion.
• Wildcard string matching via FFT.

From Convolution to Pointwise Product

Evaluate → Multiply → Interpolate

Replace 10 by $x$: $A(x)=3+2x+x^2$, $B(x)=5+4x$, and $C(x)=A(x)B(x)$ for every $x$.

• $x=2$: $A=11$, $B=13$, $C=143=11\cdot 13$.
• $x=-1$: $A=2$, $B=1$, $C=2$.

At a single point, the whole convolution collapses to one scalar multiplication.

How many points: $\mathrm{deg}C+1=\text{len}(a)+\text{len}(b)-1$. Two distinct polynomials of degree $\le d$ agreeing on $d+1$ points would have a nonzero difference of degree $\le d$ with $d+1$ roots, which is impossible.

Worked interpolation with points $\{0,1,-1,2\}$, $C$-values $\{15,54,2,143\}$:

$$\begin{array}{rl}{c}_0& =C(0)=15\\ C(1)+C(-1)& =2c_0+2c_2=56\Rightarrow c_2=13\\ C(1)-C(-1)& =2c_1+2c_3=52\\ C(2)& =15+2c_1+52+8c_3=143\Rightarrow c_3=4,c_1=22\end{array}$$

	Logarithm tables	Fourier
Transform	$\log$	Evaluate at points
Easy operation	Add	Multiply pointwise, $O(n)$
Inverse	$\exp$	Interpolate

Link to Karatsuba and Toom–Cook

Evaluating at small points is itself a real algorithm family. Karatsuba Algorithm: split into halves (degree-1 polynomials), evaluate at $\{0,1,\infty \}$, so 3 products instead of 4, giving $O(n^{1.585})$. Toom–3: points $\{0,1,-1,-2,\infty \}$. The FFT is the limit of this idea, with $N$ points chosen so the transform itself is fast.

The FFT

Even/Odd Split

$$A(x)=E(x^2)+xO(x^2)$$

For $A=[3,2,1,0]$: $E(y)=3+y$ and $O(y)=2+0y=2$. The $x$ is factored out of the odd terms, so $O$ is constant here.

Butterfly (from one pair of child values, get two parent values):

$$\begin{array}{rl}A(x)& =E(x^2)+xO(x^2)\\ A(-x)& =E(x^2)-xO(x^2)\end{array}$$

Going up = evaluating the nested expression from the inside out: $A(x)=[3+x^2\cdot 1]+x[2+x^2\cdot 0]$. The brackets depend only on $x^2$, so they are shared by $\pm x$.

Why Roots of Unity

The regrouping holds for any $x$. Roots are needed so that the set of points halves at every level.

Point set ($N=8$)	Distinct points per level	Saving
Random	8 → 8 → 8	None
$\pm 1,\pm 2,\pm 3,\pm 4$	8 → 4 → 4 → 4	One level only (squares are positive)
8th roots of unity	8 → 4 → 2 → 1	Every level

Levels and Twiddle Factors

Points are fixed by the level, not chosen by hand ($N=8$, $\omega =e^{2\pi i/8}$):

Level	Sub-polynomial size	Points
0	1	$\{1\}$
1	2	$\{1,-1\}$
2	4	$\{1,i,-1,-i\}$
3	8	every $45°$

• At each level, the $x$ values are the points in the upper half of the circle, $[0°,180°)$, and $-x$ is the antipode.
• The butterfly at $x={\omega }^k$ needs child values at ${\omega }^{2k}$, which were computed one level below.
• The multiplier $x$ is the twiddle factor.

Worked Example, $N=4$

$A=[3,2,1,0]$ at $\{1,-1,i,-i\}$:

$$\begin{array}{rl}\text{L1:}& E(\pm 1)=3\pm 1=4,2O(\pm 1)=2\pm 0=2,2\\ \text{L2:}& A(\pm 1)=4\pm 2=6,2A(\pm i)=2\pm 2i\end{array}$$

Complexity

$$T(N)=2T(N/2)+O(N)=O(N\log N),\#\text{butterflies}=\frac{N}{2}{\log }_{2}N$$

For $N=8$: 3 levels × 4 butterflies = 12.

Inverse FFT

Evaluation at the roots is the Vandermonde Matrix $V_{kj}={\omega }^{jk}$. Orthogonality of sampled waves:

$$\sum _{k=0}^{N-1}{\omega }^{k(j-l)}=N{\delta }_{jl}\Rightarrow V^{-1}=\frac{1}{N}\overline{V}$$

So interpolation is the same FFT with ${\omega }^{-1}$, divided by $N$. With arbitrary points, inverting $V$ is expensive and ill-conditioned.

Where the Fourier Part Is

Evaluation at Roots = DFT

$$A({\omega }^k)=\sum _j{a}_j{e}^{2\pi ijk/N},e^{i\theta }=\cos \theta +i\sin \theta$$

The powers of ${\omega }^k$ are a sampled wave with $k$ cycles. So $A({\omega }^k)$ = (sequence · probe wave), i.e. how much of frequency $k$ is in the sequence.

• $A(1)={\sum }_j{a}_j$: frequency 0 (DC, the mean).
• $A(-1)={\sum }_j(-1{)}^j{a}_j$: the highest frequency (Nyquist).

Standard DFT convention uses $e^{-2\pi i}$: same thing, traversing the circle the other way.

Algebra view	Signal view
$C(x)=A(x)B(x)$ at every point	Convolution theorem: $\mathcal{F}(a\ast b)=\mathcal{F}a\cdot \mathcal{F}b$
Evaluate at the roots of unity	Measure each frequency
Even/odd split + halving lemma	Cooley and Tukey's FFT (1965), anticipated by Gauss

Eigenvector View

Convolution with $b$ is a Circulant Matrix $C_b$, and the DFT diagonalises it:

$$C_b=F^{-1}\mathrm{diag}(Fb)F$$

The algebra explains why it's fast; the eigenvector view explains why it works for any LTI system.

Big-Integer Multiplication: Full Algorithm

Pipeline

1. Digits units-first; pad to a power of 2 with $N\ge \text{len}(a)+\text{len}(b)-1$.
2. FFT both, multiply pointwise, inverse FFT.
3. Round to nearest integer (removes floating-point noise).
4. Carry, $O(n)$.

import cmath

def fft(a, inv=False):
    n = len(a)
    if n == 1:
        return a
    E, O = fft(a[0::2], inv), fft(a[1::2], inv)   # even / odd split
    s = -1 if inv else 1
    out = [0] * n
    for k in range(n // 2):
        w = cmath.exp(s * 2j * cmath.pi * k / n)   # twiddle factor
        out[k] = E[k] + w * O[k]                   # A(x)
        out[k + n // 2] = E[k] - w * O[k]          # A(-x)
    return out

def multiply(x, y):
    a = [int(d) for d in str(x)[::-1]]             # units first
    b = [int(d) for d in str(y)[::-1]]
    n = 1
    while n < len(a) + len(b) - 1:                 # avoid wrap-around
        n *= 2
    A = fft(a + [0] * (n - len(a)))
    B = fft(b + [0] * (n - len(b)))
    c = fft([p * q for p, q in zip(A, B)], inv=True)
    c = [round((v / n).real) for v in c]           # 1/N, kill float noise
    return sum(ck * 10**k for k, ck in enumerate(c))   # carry

assert multiply(123, 45) == 5535

State of the Art

• Number-Theoretic Transform: replace $\mathbb{C}$ with $\mathbb{Z}/p$ (roots of unity exist mod suitable primes), so the arithmetic is exact.
• Schönhage and Strassen's algorithm (1971): rings $\mathbb{Z}/(2^m+1)$, where 2 is a root of unity. $O(n\log n\log \log n)$.
• Harvey and van der Hoeven's algorithm (2019): $O(n\log n)$, conjectured optimal.
• GMP: schoolbook → Karatsuba → Toom–Cook → FFT as the size grows.

Drawbacks

• Constant factors: the FFT wins only for large $n$.
• Floating-point precision: rounding fails at large $n$ or large digit bases → NTT.
• Wrap-around if padding is forgotten.

CNN Layers Through the Fourier Lens

Conv Layer as Per-Frequency Matrices

Without the activation, a conv layer is LTI. In the Fourier domain:

$${\hat{y}}_{c_{\text{out}}}(\xi )=\sum _{c_{\text{in}}}{\hat{w}}_{c_{\text{out}},c_{\text{in}}}(\xi ){\hat{x}}_{c_{\text{in}}}(\xi )$$

So a conv layer is one $C_{\text{out}}\times C_{\text{in}}$ matrix per frequency, with no mixing across frequencies.

• Fast conv via FFT for large kernels; 3×3 kernels use Winograd instead.
• Small kernels give blunt filters: compact support forces smooth frequency responses (uncertainty principle). Depth multiplies responses, which sharpens them.
• Exact singular values = singular values of the per-frequency matrices → Lipschitz control (Sedghi et al., ICLR 2019).
• Fourier Neural Operator (Anandkumar's group, ICLR 2021): learn the per-frequency weights directly; keep low modes.

Where the Picture Breaks

Component	Breaks	Effect
ReLU	Linearity	Creates new frequencies; can't Fourier the whole net
Stride	Shift-invariance	Downsampling → Aliasing; fixed by blur-before-subsample (Richard Zhang, ICML 2019)
Zero padding	Circularity	Boundary effects; FFT assumes wrap-around

Connections

Adjacent and Distant Ideas

• Logarithms: the original transform → cheap operation → inverse trick.
• Probability: characteristic functions turn sums of independent RVs into products.
• Heat equation: Fourier's original problem; each mode decays independently at a rate $\propto k^2$.
• Graph Fourier Transform: on a graph, Laplacian eigenvectors replace sines; this is the basis of spectral GNNs.
• Consensus dynamics: averaging on a ring or lattice is LTI, so disagreement decomposes into modes; long-wavelength disagreement decays slowest. This is relevant when thinking about how norms spread over a lattice of agents (norm-emergence).

Checklist

• State the LTI test and why exponentials are eigenfunctions of shift-commuting operators.
• Classify a new case as Fourier / not Fourier, asking about the system.
• Derive $c_k={\sum }_{i+j=k}{a}_i{b}_j$ from multiplying powers of 10.
• Explain the kernel flip and why ML layers skip it.
• Padding rule $N\ge \text{len}(a)+\text{len}(b)-1$ and what goes wrong without it.
• Evaluate → multiply → interpolate; why $\mathrm{deg}+1$ points suffice.
• Even/odd split, the butterfly, and why $\pm$ real pairs fail at level 2.
• Which points each level uses, and the twiddle factor.
• Count: $\frac{N}{2}{\log }_{2}N$ butterflies.
• Inverse FFT via $V^{-1}=\frac{1}{N}\overline{V}$.
• Interpret $A({\omega }^k)$ as frequency content; $A(1)$ = DC, $A(-1)$ = Nyquist.
• Conv layer = per-frequency channel-mixing matrices; where ReLU, stride and padding break it.
• Why the FFT needs NTT at scale; Schönhage–Strassen vs Harvey–van der Hoeven.

I ran the Python and it matches Python's built-in multiplication on random 60-digit numbers. The note also includes the inverse FFT and the A(1) answer, which we hadn't reached in the walkthrough.

Fourier's heat problem — Step 1: Heat on a ring of cells

Fourier studied how temperature evolves in a solid, including a metal ring. Start with a ring cut into cells. At each time step, every cell moves a little toward the average of its two neighbours:

u_j ← u_j + α·(u_{j−1} − 2u_j + u_{j+1})

A cell hotter than its neighbours cools down, and a colder one warms up. Written as a kernel, the update is a sliding window with weights [α, 1−2α, α]: the same kernel applied at every cell, wrapping around the ring.

Example: 4 cells, α = 0.25, starting at [4, 0, 0, 0], so all the heat is in one spot.

• u₀ = 4 + 0.25·(0 − 8 + 0) = 2
• u₁ = 0 + 0.25·(4 − 0 + 0) = 1
• u₂ = 0 + 0.25·(0 − 0 + 0) = 0
• u₃ = 0 + 0.25·(0 − 0 + 4) = 1

After one step the heat is [2, 1, 0, 1]: the peak has spread to its neighbours, and the total is still 4.

Check (two parts):

1. Is this update linear and shift-invariant? What does that tell you, given everything above?
2. Compute the next step from [2, 1, 0, 1].

Old Note

The algorithm has been the same, some ideas are in Fourier Series, but architectures change, which means there are new ways to make this algorithm even faster.

Example of transforms

We have learned in Algebra lineare numerica, Cambio di Base that linear transforms are usually a change of basis. They are matrix vector multiplications (additions and multiplications by constants). The optimizations are based on what sorts of transforms we have (e.g. Sparse Matrix Vector Multiplication, or dense versions). The same idea applies also for Fourier transforms.

Discrete Fourier Transforms

Primitive roots

It uses primitive roots of $n$ as coefficients. Discrete Fourier transform computes:

$$\begin{array}{rl}{y}_k& =\sum _{l=0}^{n-1}{x}_l\cdot e^{-2\pi i\frac{kl}{n}}\\ & =\sum _{l=0}^{n-1}{x}_l\cdot {\omega }^{kl}\end{array}$$

Where $\omega =e^{-2\pi i/n}$ is the primitive root of $n$.

The matrix $DF{T}_n$ is a square matrix of size $n$ and the entries are given by:

$DF{T}_n[k][l]={\omega }^{kl}=e^{-2\pi i\frac{kl}{n}}$

Example of DFT matrices

At each index we will have the $ij$ th unit of root. For example if $n=2$ we know that the root is $-1$, then the matrix would be:

$$DF{T}_2=\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)$$

This has two adds (one add and sub), it would be quite inefficient to write the multiplication part. This is usually called butterflies in signal processing, due to the nature of the circuit.

for $n=4$, the primitive root is $i$ we have:

$$DF{T}_4=\left(\begin{array}{cccc}1& 1& 1& 1\\ 1& i& -1& -i\\ 1& -1& 1& -1\\ 1& -i& -1& i\end{array}\right)$$

If we assume complex everything, it has 12 additions and 4 multiplications. This matrix can be factorized, then we can add it as a series of matrix vector products.

Decomposing into sparse MVM

If $T=T_1{T}_2...T_k$, we can write:

$$\begin{array}{rl}{t}_1& =T_{1}x\\ t_2& =T_2{t}_1\\ ⋮& \\ t_k& =T_k{t}_{k-1}\\ y& =t_k\end{array}$$

This makes sense if $T_i$ are sparse (see Sparse Matrix Vector Multiplication, and not too many $m$ steps).

You can use Group Theory to find this kind of transformations.

The above representation can be rewritten using matrix algebra:

$$DF{T}_4=(DF{T}_2\otimes I_2)\cdot \text{diag}(1,1,1,i)\cdot (I_2\otimes DF{T}_2)L_2^4$$

You can see that it is also quite easy to implement it in hardware.

It can be similar for size $16$ and radix $4$.

Other transforms

IMDCT and DCT are useful transforms for MPEG and JPEG format, while DFT and its inverse can be considered universal transformation tools.

FFT transformations

General Radix Recursive Cooley-Tukey FFT

This is the most general version of DFT, it has a similar structure compared to the previous ones.

$$\begin{array}{rl}DF{T}_n& =(DF{T}_{n/r}\otimes I_r)\cdot \text{diag}(1,1,...,{\omega }^0,{\omega }^1,...,{\omega }^{r-1})\cdot (I_{n/r}\otimes DF{T}_r)L_{n/r}^n\end{array}$$

The diagonal matrix has the roots of unity. This is called the decimation in time frequency. It is possible also to use the decimation in frequency, which is just the transpose. We can decompose everything until we have size two butterfly operations (bunch of adds).

Cost of DFT

Independently of recursion (the radix), we have a number of

• $n{\log }_2(n)$ additions and $n{\log }_2(n)/2$ multiplications, but to be real we would need to consider real additions, so we have to multiply everything by two, to get the adds. For multiplication, we have four multiplication and two additions (image it as a matrix) (in reality you use the structure of roots of unity, so in reality you have less adds). Usually best is around radix 8
• You will have about $\le 3n{\log }_2(n)$ real additions and $\le 2n{\log }_2(n)$ complex multiplications, so upper bound of 5 operations. For a long time radix 2 iterative version was the best version, after superscalar processors like Skylake Microprocessor it has better versions. There is a practical lower bound analysis, but that is like Strassen Fast Linear Algebra, which is not much useful.

Recursive vs iterative FFT

The recursive version is good for re-use, has better locality. Triple loop is like the blocked MMM in Fast Linear Algebra. If we have big stages computed, we can continue with the next stage immediately.

Iterative similar, the op count is the same, is slower because of the (log n passes through the data). The code change is quite difficult. If we have much data, then at each iteration we need to load the whole data which could be expensive.

Decimation in Time and Frequency

TODO

Other FFTs

Pease FFT

Highly parallelized with maps and shuffles.

Stockham FFT

It is good for vector computers of the 80s, where we have long vectors

Six-Step FFT

Multi-core FFT

Good for Sharing memory.

Fastest Fourier Transform of the West

Also known as FFTW, we say.

Ideas for FFTW

Merging Steps

1. In the original version, the original permutation could be very expensive, to load everything. There is no need to keep it as a sperate step, we can fuse it in the next step.
   1. This fusion would need more data, because it needs to write to another place (no-inplace computation!).
2. The scaling can also be fused, with the previous or next step (just change the kernel).

Strides will change as we recurse for the DFTrec, and we can also interchange input and output to save space. Because of the power of two strides, often it is very difficult to reach peak performance.

Constants

We precompute the constants of the sins and cosines to save time (just need what roots of unity you need in order to compute it). They just depend on the input size, so they precompute it at the beginning of the function.

Block un-rolling

We un-roll the basic blocks at the end to get longer codes, less recursion (size 32 is good basecase for example).

Program Generator

We have DAG generator (for many types of FFT), a simplifier for the DAG (specific for the FFT computation).

• An example of a simplification is when you have $x-y$ and $y-x$, an optimization would be flipping one after you have computed one.
• Also constants come in pairs (minus and positive version).
  • Reduces register pressure.