Fourier Neural Operators (FNOs)

The previous post built the Fourier transform: a signal is a sum of pure oscillations, and the transform recovers the amplitude and phase of each one. This post uses that machinery to build an actual architecture. The Fourier Neural Operator (FNO) is a neural operator in the sense of the neural operators post — it approximates the solution operator of a PDE, and it attaches its weights to resolution-independent objects rather than to grid points. In FNO those objects are frequency modes.

The post starts from the problem FNO exists to solve, then walks through every layer of the architecture in the order data flows through it, tracking the tensor shape at each step.

A convolution takes a short list of numbers, called the kernel, and slides it along a longer list, computing one weighted sum per position.

Where $s$ denotes the input signal as a list of numbers, $s_{i+m}$ denotes its entry at index $i+m$, $k_m$ denotes the kernel weight at offset $m$, $r$ denotes the kernel's half-width (so a three-entry kernel has $r=1$ and offsets $m\in \{-1,0,1\}$), and $o_i$ denotes the output at index $i$.

Take the kernel $k=[1,-2,1]$ and the signal $s=[2,5,3,8,1]$. To compute the output at index 2, center the kernel on $s_2=3$ and take the weighted sum of that entry and its two neighbours.

Where the three products pair the kernel weights $k_{-1}=1$, $k_0=-2$, $k_1=1$ with the signal entries $s_1=5$, $s_2=3$, $s_3=8$, and $o_2$ denotes the resulting output value at index 2.

That kernel is not arbitrary. It is Post 2's three-point curvature stencil $(u_{i+1}-2u_i+u_{i-1})/\mathrm{\Delta }{x}^2$ without the division by $\mathrm{\Delta }{x}^2$, where $u_i$ denotes the field value at grid point $i$ and $\mathrm{\Delta }x$ denotes the grid spacing. Every finite-difference update in Posts 2 through 6 was a convolution. A convolutional layer in a neural network performs the same operation, with the kernel entries learned from data instead of derived from the second derivative stencil.

A kernel of half-width $r$ combines a point with $r$ neighbours on each side — a total of $k=2r+1$ points. Track what a single output value has "seen" as more of these layers stack, working through $k=3$ (so $r=1$) one layer at a time.

Layer 0 — the raw input. Point $i$ has seen only itself.

Layer 1 — one convolution applied. Point $i$'s new value is a weighted sum of $i-1,i,i+1$. It has now seen 3 original input points.

Layer 2 — a second convolution applied to layer 1's output. Point $i$'s new value combines three layer-1 values, not three raw inputs: its own layer-1 value (which already saw $i-1,i,i+1$), its left neighbour's layer-1 value at $i-1$ (which already saw $i-2,i-1,i$), and its right neighbour's layer-1 value at $i+1$ (which already saw $i,i+1,i+2$). Taking the union of those three ranges, point $i$ has now seen $i-2$ through $i+2$ — width 5.

Generalizing this growth pattern across $L$ layers gives the receptive field: the number of input positions that can influence a single output position after $L$ layers.

Where $R(L)$ denotes the receptive-field width in grid points, $L$ denotes the number of stacked convolutional layers, and $k$ denotes the kernel width in grid points as defined above. Checking it against the two layers just worked through: $R(1)=1+1\cdot 2=3$ matches layer 1's width of 3, and $R(2)=1+2\cdot 2=5$ matches layer 2's width of 5.

Now put a number on it. For a domain discretized into 1000 grid points, connecting one end to the other requires $R(L)\ge 1000$, so $1+2L\ge 1000$, so roughly 500 layers. That is the cost of building global reach out of local, neighbour-only operations.

This matters directly for the wave equation. Post 6 established that information propagates at speed $c$, and that a stable explicit solver must respect the CFL condition $c\mathrm{\Delta }t\le \mathrm{\Delta }x$, where $c$ denotes wave speed, $\mathrm{\Delta }t$ denotes the time step and $\mathrm{\Delta }x$ denotes grid spacing. A network that has to reproduce a large time step in one forward pass needs a wavefront to move many grid cells at once, and a stack of small kernels can only deliver that with enormous depth.

The Fourier transform post gave the discrete Fourier transform.

Where $S_k$ denotes the complex coefficient of mode $k$, $s_i$ denotes the signal sample at grid index $i$, $N$ denotes the number of samples, $k$ denotes how many full oscillations the reference wave completes across the whole sampled domain, and $j$ denotes the imaginary unit with $j^2=-1$.

Every sample $s_i$ appears in that sum. A single coefficient is already a global summary of the field. Editing it and transforming back changes the field everywhere. That is the entire motivating idea of the Fourier Neural Operator: perform the layer's mixing step on coefficients rather than on grid neighbourhoods, and one layer reaches the whole domain.

The Fourier transform post established the mechanism; this section only fixes the reading of the two numbers stored per mode, because the architecture manipulates them directly.

Expanding the reference wave with Euler's formula splits the sum into two real sums.

Where $x$ denotes the real part of $S_k$, measuring how strongly the signal resembles a cosine at mode $k$, and $y$ denotes the sine similarity, which enters the coefficient with a minus sign so that $S_k=x-jy$; the other symbols are as above.

The two numbers can be read either as a pair or in polar form. The magnitude $|S_k|=\sqrt{x^2+y^2}$ is the amplitude of that oscillation in the signal, and the angle $\arg (S_k)$ is its phase — where the oscillation's peak sits relative to the origin. Both symbols are as defined above.

Check this on a signal that is exactly one cosine. Take $s=[1,0,-1,0]$ with $N=4$, which is $\cos (2\pi i/4)$ sampled at $i=0,1,2,3$, and evaluate mode $k=1$. The reference wave takes the values $1,-j,-1,j$ at those four indices.

Where $S_1$ denotes the coefficient of mode 1, the first factor in each product is the signal sample $s_i$ and the second is the reference-wave value $e^{-j2\pi i/4}$ at that index.

The result is the pair $(2,0)$: all cosine similarity, no sine similarity, phase zero. The magnitude 2 equals the signal's amplitude of 1 times $N/2$, the standard scaling of the unnormalized transform.

The concrete problem for the rest of this post is one step of the wave equation's solution operator in two spatial dimensions. Post 4 established that the wave equation needs a two-field state: the displacement field $U$ and the velocity field $V=\partial U/\partial t$, where $t$ denotes time. A single snapshot of displacement alone cannot determine the next state.

Discretize the domain on a 100 by 100 grid. The input to the network is then a tensor of shape (2, 100, 100): two channels, one per field, each holding a 100 by 100 array of grid values. A channel is a full spatial field; a grid cell holds one number per channel. At the very start, the two numbers at a given cell are that cell's displacement and velocity.

Every operation below is stated twice: what it does to the tensor shape, and what it does to the numbers sitting at one individual grid cell.

The first layer is a lifting layer: a linear map applied independently at every grid cell, which increases the channel count.

Where $a(x)$ denotes the input channel vector at grid position $x$ (here the two numbers $(U,V)$ at that cell), $P$ denotes a learnable weight matrix of shape (4, 2), $p$ denotes a learnable bias vector of length 4, and $v(x)$ denotes the resulting lifted vector of 4 numbers at the same position.

Take a cell whose values are $(U,V)=(3,-2)$ and a weight matrix

Where each row of $P$ produces one output channel, and the two entries in a row are the weights applied to the displacement and the velocity respectively. Taking the bias to be zero:

Where the input column vector holds the cell's displacement and velocity, and the output column vector holds the four lifted values at that cell.

The same matrix $P$ is applied at all 10,000 cells, so the tensor shape goes from (2, 100, 100) to (4, 100, 100).

Two channels are too narrow to carry the intermediate representation. The layers that follow mix channels together, and with only two channels there are just four mixing coefficients available per frequency mode — far too few to represent a useful learned transformation. Widening also breaks the interpretation: after lifting, the four numbers at a cell are learned features, not displacement and velocity. Nothing forces channel 1 to remain $U$. Real implementations lift to 32 or 64 channels; 4 is used here so the arithmetic stays writable.

Because the lifting map is the same at every position and involves no neighbours, it is resolution independent in the sense of the neural operators post: the same weights apply on a grid of any size.

The lifted tensor now enters the core of the architecture. An FNO block splits into two parallel paths. The spectral branch is the global one, and it runs in four steps.

Each of the 4 channels is transformed independently with a two-dimensional discrete Fourier transform. The Fourier transform post gave the two-dimensional reference wave; the transform is the same dot-product-against-a-reference-wave construction, now with a frequency index per axis.

Where $\hat{v}(k_1,k_2)$ denotes the complex coefficient of the two-dimensional mode $(k_1,k_2)$ for one channel, $v_{i,l}$ denotes that channel's value at grid cell $(i,l)$, $N_1$ and $N_2$ denote the grid size along each axis (both 100 here), $k_1$ denotes the number of oscillations across the first axis, $k_2$ the number across the second, and $j$ the imaginary unit.

A naive implementation would produce a 100 by 100 array of complex coefficients per channel. Half of that is redundant. The input field is real-valued, and for real input the coefficients obey a conjugate symmetry.

Where the overline denotes complex conjugation, meaning the real part is unchanged and the imaginary part flips sign, and $(-k_1,-k_2)$ denotes the mode with both frequency indices negated.

The reason is that a real signal has no imaginary content to encode, so the negative-frequency half carries no new information: the sine similarity at mode $-k$ is exactly minus the sine similarity at mode $+k$, and the cosine similarities are equal. This is the same statement as the sampling limit: $N$ real samples carry exactly $N/2+1$ independent frequencies, running from the constant mode $k=0$ up to the Nyquist mode $k=N/2$, the fastest oscillation the grid can represent — one full period every two grid points.

The Real-Values FFT exploits this by returning only the non-redundant half along the last axis. For a 100 by 100 field that is 51 columns, since $100/2+1=51$. The tensor is now (4, 100, 51), holding one complex number per channel per mode. This can honestly trip you up a bit, so I encourage you to make a deep dive into this topic, even though it is not necessary. Just remember, that our input value has no imaginary part. Throwing the signal into the DFT produces complex coefficients. However, the imaginary part does not contain any new information. Therefore, for the output signal to be real again, the imaginary parts of the coefficients must be conjugate symmetric for the imaginary parts to cancel out.

The next step discards most of those modes. Fix a mode count $M$ per axis and keep only the lowest $M$ frequencies; take $M=10$. The kept block has shape (4, 20, 10).

The asymmetry between 20 and 10 comes from the folding in Step 1. The second axis was folded, so its negative frequencies are already gone and keeping the 10 lowest means 10 columns. The first axis was not folded, so it still stores both signs: its 100 rows run from mode 0 upward through the Nyquist mode and then continue as the negative modes counting back toward $-1$. Keeping the 10 lowest in each direction means 10 rows from the top of the array and 10 from the bottom, hence 20.

Low modes are kept because they carry almost all of a field's shape: the coarse structure holds most of the amplitude, and coefficient magnitudes decay as frequency rises.

The more important property is what a fixed mode count buys. Mode 7 names a shape — a wave that repeats 7 times across the whole domain — not a grid location. Sample that shape with 100 points or with 400 and it is still the same 7 repetitions, just measured more finely. So a weight attached to mode 7 stays attached to the same physical wave regardless of resolution, unlike a weight attached to "grid point 7," which points somewhere different at every resolution. That is why keeping the lowest $M$ modes always costs the same $M$ weights, with no dependence on the grid size — the discretisation invariance of the neural operators post, and the reason FNO counts as a neural operator rather than a resolution-locked network.

The cost is just as direct. Keep only 10 modes, and mode 11 upward is gone, not approximated. A sharp wavefront is made of exactly those discarded modes reinforcing each other at one spot and cancelling everywhere else; remove them and the edge slumps into a slope. Post 2's heat equation was going to smooth itself anyway, so this costs it nothing. The wave equation exists to keep its edges sharp, so the same cut fights the physics directly — FNO's known weak point on hyperbolic problems.

Each kept coefficient is multiplied by a learnable complex weight.

Where $x+yj$ denotes the coefficient entering the layer, $a+bj$ denotes the learnable weight at that mode, and $j$ denotes the imaginary unit with $j^2=-1$. Both sides are pairs of real numbers, so in practice this is stored and computed as two real tensors — no complex data type is needed.

Multiplying by a complex number rescales the coefficient's magnitude and rotates its phase, the two quantities defined in Section 2: the oscillation's amplitude changes and its peak shifts, while its frequency stays fixed.

So essentially such complex weight lets the network turn one frequency up or down, and slide that whole oscillation sideways — all of its repeating peaks move together, since a phase shift is just the same wave pattern translated along the grid — independently of every other frequency. Give it one such weight per mode, and it can learn to keep the frequencies that matter for reproducing the dynamics and suppress the ones that don't and by that, essentially move the fields forward in time.

Channels mix here the same way as in a normal convolutional layer: each output channel is a weighted combination of all input channels. The only difference is that every kept mode gets its own separate weights.

where $d$ denotes the output channel, $c$ denotes an input channel, and $(k_x,k_y)$ denotes one kept mode.

This makes the weight tensor's shape (4, 4, 20, 10). The coefficient tensor stays (4, 20, 10).

Post 7's widget runs this weight step in miniature: one channel instead of four, sine shapes on a string with fixed ends instead of complex waves, and every weight set to the exact solution operator for one time step $T$, where $T$ denotes the length of that step. Under the heat equation each weight only turns its mode down. Under the wave equation it acts on two numbers per mode, the displacement coefficient and the velocity coefficient, and turns that pair by a fixed angle: the same rescale and rotate as the complex weight above, applied to a different pair of numbers. In its rollout experiment, opened below, every output becomes the next input, so the per-mode weights alone carry the field forward in time.

The multiplied block has shape (4, 20, 10), but the real-valued inverse Fast Fourier Transform (RFFT) expects a full coefficient array of shape (4, 100, 51). As we know, we intentionally discarded the high modes, so the missing entries are simply filled back up with zeros as coefficients.

This padding happens in frequency space, before the inverse transform. It does not restore the discarded content and therefore states that the model's output has no content at those frequencies. The inverse transform then returns to the grid, producing shape (4, 100, 100). The output of the spectral branch is a bandlimited field: it has no content above the cutoff frequency, and therefore cannot reproduce sharp edges. Note how this might become a problem for hyperbolic PDEs, which are known to maintain sharp wavefronts rather than smoothing them out. The local branch is introduced to address this problem.

Running alongside the spectral branch is a second path applied to the same lifted input. It is a linear map applied at each grid cell independently, with no neighbours involved — a 1 by 1 convolution.

Illustrate with two channels for readable arithmetic. Take the weight matrix and a cell holding $(3,5)$:

Where $W$ denotes the learnable local weight matrix, the input column vector holds the channel values at one grid cell, and the output column vector holds that cell's transformed channel values. In the running example $W$ is 4 by 4, holding 16 weights.

The same $W$ is used at every one of the 10,000 positions, and the shape stays (4, 100, 100).

This branch reads the untruncated input. The spectral branch saw only 200 modes, whereas the local branch sees the field as it actually is at each cell, including whatever high-frequency content the truncation in the spectral branch removed. It is therefore the one route by which fine detail can reach the block's output. There is a major limitation to this though: being pointwise, it can only pass a cell's own information forward. Independently of it's neighbours. Because of this, the local branch is literally blind to shapes, slopes and curvatures.

The obvious question is why this branch is not a 3 by 3 convolution, which would at least see neighbours. The answer is resolution invariance. A 3 by 3 kernel spans a fixed number of grid cells, not a fixed physical distance. On a unit-length domain with 100 points, $\mathrm{\Delta }x\approx 0.01$ and the kernel reaches about $0.02$ in physical units; on the same domain with 400 points it reaches about $0.005$, a quarter of the footprint, where $\mathrm{\Delta }x$ denotes grid spacing. The same learned weights would then implement a different physical operation at each resolution, which is exactly the property FNO's mode-indexed weights were built to avoid.

Hybrids that attach a U-Net to the FNO block do recover local detail and measurably improve sharp-feature accuracy, but they lose resolution invariance at their downsampling and upsampling steps. The direction a fix would take is to define the local operation as a continuous kernel function evaluated at the actual grid spacings $\mathrm{\Delta }x$ and $\mathrm{\Delta }y$. However, this is still an open problem, and not part of this blog post yet. In here, we discuss the original FNO architecture, which is resolution invariant through the mode truncation and 1x1 kernel convolutions.

The two branches produce tensors of identical shape (4, 100, 100), and are combined by elementwise addition followed by a pointwise nonlinearity. That is the complete FNO block.

Where $v_{\ell }(x)$ denotes the channel vector at position $x$ entering block $\ell$, $v_{\ell +1}(x)$ the vector leaving it, $W$ the local branch's weight matrix, $b$ its bias, $\mathcal{F}$ and ${\mathcal{F}}^{-1}$ the forward and inverse Fourier transforms, $R$ the mode-indexed complex weights including the truncation and zero-padding, and $\sigma$ a pointwise nonlinearity such as GELU.

Numerically, at one grid cell: suppose the spectral branch outputs $(1.4,-0.3,2.0,0.5)$ and the local branch outputs $(-0.4,-0.7,-1.0,0.5)$. Adding channel by channel gives $(1.0,-1.0,1.0,1.0)$, and applying a ReLU nonlinearity gives $(1.0,0.0,1.0,1.0)$. The same computation runs at all 10,000 cells, and the shape stays (4, 100, 100) throughout.

A single block applies one global mixing step and one nonlinearity. A full FNO stacks $L$ of them, typically 4, where $L$ denotes the number of blocks. Each block holds its own independent parameters in both branches: its own mode-indexed weights $R$ and its own local matrix $W$, both as defined above. Nothing is shared between blocks. The channel count and grid size stay at (4, 100, 100) from the first block's output to the last.

Depth here buys something different from depth in a CNN. A CNN needs depth to grow its receptive field; an FNO already reaches globally in one block. Depth instead buys composition — repeated alternation between global linear mixing and pointwise nonlinearity — which is what allows the model to represent nonlinear solution operators rather than a single linear reweighting of modes.

https://ar5iv.labs.arxiv.org/html/2010.08895

The last block's output has 4 channels of learned features. The projection layer maps those back to the 2 physical fields.

Where $v_L(x)$ denotes the 4-channel vector at grid position $x$ leaving the final block, $Q$ denotes a learnable weight matrix of shape (2, 4), $q$ denotes a learnable bias vector of length 2, and $u(x)$ denotes the predicted displacement and velocity at that position. In practice the projection is often two linear maps with a nonlinearity between them, widening to an intermediate size such as 128 before narrowing to 2, which gives the readout enough capacity to be more than a fixed linear combination.

The shape returns to (2, 100, 100): displacement and velocity on the original grid, one time step ahead.

The full path through the network, in shapes:

(2, 100, 100) → lifting → (4, 100, 100) → block × $L$ → (4, 100, 100) → projection → (2, 100, 100)

Each block splits into a spectral branch — transform to (4, 100, 51), truncate to (4, 20, 10), multiply by mode-indexed complex weights, zero-pad, inverse transform to (4, 100, 100) — and a pointwise local branch at the same shape, with the two summed and passed through a nonlinearity.

Real configurations are of the same order. Grids run from 64 to a few hundred points per side. The mode count does not scale with the grid: a model trained at 64 points per side with 12 modes is evaluated at 256 points per side with the same 12 modes and the same weights. That is what discretisation invariance means in practice.

The weakness is however, that truncation to $M$ modes obviously deletes high frequency details, resulting in a smoothed output. For hyperbolic PDEs such as the wave equation, this is particularly problematic. The natural next question is whether an architecture can keep FNO's resolution independence without paying for it by losing high-frequency content.