From a PDE solve to operator learning¶
A familiar problem is to solve a differential equation for one forcing function:
$$ \boxed{\mathcal L u=f},\qquad \mathcal L=-\frac{d^2}{dx^2},\qquad u(0)=u(1)=0. $$
Given $f(x)$, we find $u(x)$. For example, $f(x)=1$ gives $u(x)=x(1-x)/2$. A different forcing requires a different solve.
Operator learning asks for the rule that performs all these solves. With a fixed differential operator $\mathcal L$ and boundary conditions, write
$$ \boxed{\mathcal G:\ f\longmapsto u,\qquad u=\mathcal G(f)=\mathcal L^{-1}f.} $$
The distinction matters: $\mathcal L$ takes a solution to a forcing; its inverse $\mathcal G$ takes a forcing to a solution. Given many observed pairs $(f^{(i)},u^{(i)})$, a neural operator ordinarily learns $\mathcal G$. Learning the unknown differential law $\mathcal L$ itself is a different problem: operator identification.
What makes the inverse a genuine map? For every admissible forcing, the PDE must have one well-defined solution. Here the two boundary conditions guarantee uniqueness. Without suitable boundary conditions, $-u''=f$ has infinitely many solutions differing by an affine function.
A neural network only processes finitely many numbers. So the first obstacle is not the neural network: how can finite numbers represent a whole input function?
1. Can finite sensors describe a function?¶
Place $m$ sensors at positions $x_1,\ldots,x_m$. The whole input curve is replaced by a vector:
$$ \mathbf s(f)=\big[f(x_1),\ldots,f(x_m)\big]\in\mathbb R^m. $$
Connect the sampled values to make a reconstructed curve $f_m(x)$. We can measure what is lost using the supremum norm:
$$ \boxed{\|f-f_m\|_\infty=\max_{x\in[0,1]}|f(x)-f_m(x)|.} $$
This is the greatest vertical separation between the two curves, anywhere on the interval. Saying $\|f-f_m\|_\infty<\varepsilon$ means every point is reconstructed to within our tolerance $\varepsilon$.
Can a finite number of sensors guarantee that accuracy for every continuous function? Try the spike below. It is placed exactly between two sensors. Narrow it until every measurement reads zero, even though the true curve still reaches one.
interactive(
controls=[
dict(key='m', label='Sensors m', min=4, max=22, step=2, value=10),
dict(key='w', label='Spike width', min=0, max=26, step=0.1, value=5, play='bounce', duration=7,
hint='Play: widen and narrow the spike'),
dict(key='eps', label='Tolerance ε', min=0.05, max=0.85, step=0.05, value=0.20),
],
panels=[
dict(title='A curve and its sensor reconstruction', xlabel='position x', ylabel='function value'),
dict(title='The error is the largest vertical gap', xlabel='position x', ylabel='absolute error'),
],
draw='''
var m = P.m, xs = linspace(0, 1, m), h = 1 / (m - 1);
var j = (m - 2) / 2, center = (xs[j] + xs[j + 1]) / 2; // m is even, so the spike sits at x = 1/2
var halfwidth = (0.04 + 0.045 * P.w) * h;
function spike(x) { return max(1 - abs(x - center) / halfwidth, 0); }
var X = grid(1201, xs.concat([center - halfwidth, center, center + halfwidth]));
var f = X.map(spike), samples = xs.map(spike), recon = interp(X, xs, samples);
var difference = X.map(function (x, i) { return abs(f[i] - recon[i]); }), k = argmax(difference), err = difference[k];
panel[0].draw({ylim: [-0.07, 1.13], series: [
{type: 'line', x: X, y: f, color: 's1', width: 2.5, name: 'f(x)'},
{type: 'line', x: X, y: recon, color: 's2', dash: '6 4', name: 'fₘ(x)'},
{type: 'dots', x: xs, y: samples, color: 's2', name: 'sensors'}]});
panel[1].draw({ylim: [-0.07, 1.15], series: [
{type: 'line', x: X, y: difference, color: 's3', width: 2.5, name: '|f − fₘ|'},
{type: 'hline', y: P.eps, color: 's2', width: 1.6, dash: '6 4', name: 'ε'},
{type: 'dots', x: [X[k]], y: [err], color: 's3'}]});
readout('‖f − fₘ‖∞ ≈ ' + fmt(err, 3) + ' ' + (err >= P.eps ? '(above tolerance)' : '(below tolerance)'));
''')
Straight lines, staircases, and a genuine jump¶
The sensor vector does not dictate how we reconstruct the curve. Two simple choices are to connect samples with straight lines or assign each point its nearest sensor value, making a staircase. Both use exactly the same measurements.
Change the input below. For a smooth curve, more sensors make both approximations better. A fixed narrow spike may disappear on a coarse grid and reappear as sensors are added. In the earlier figure, the spike could move to avoid any fixed grid. A true step is different: it jumps abruptly, so it does not belong to our compact family of continuous functions. Unless the jump's location is known or captured, fixed-grid samples cannot guarantee a uniformly small error near it.
interactive(
controls=[
dict(key='family', label='Input curve', options=['Smooth wave', 'Fixed spike', 'Step (jump)'], value='Smooth wave'),
dict(key='m', label='Sensors m', min=5, max=61, step=2, value=11, play='loop', duration=8),
],
panels=[
dict(title='One sensor vector, two reconstructions', xlabel='position x', ylabel='function value'),
dict(title='Largest error anywhere', xlabel='position x', ylabel='absolute error'),
],
draw='''
function example_input(x) {
if (P.family === 'Smooth wave') return 0.65 * sin(PI * x) + 0.17 * sin(7 * PI * x);
if (P.family === 'Fixed spike') return max(1 - abs(x - 0.463) / 0.026, 0); // fixed location and width
return x >= 0.467 ? 1 : 0; // a genuine jump
}
var m = P.m, sx = linspace(0, 1, m), ss = sx.map(example_input), cuts = [];
for (var i = 0; i + 1 < m; i++) { var b = (sx[i] + sx[i + 1]) / 2; cuts.push(b - 1e-9, b); }
var X = grid(1601, sx.concat(cuts, [0.437, 0.463, 0.489, 0.467 - 1e-9, 0.467]));
function nearest_value(x) { return ss[Math.min(m - 1, Math.floor(x * (m - 1) + 0.5))]; } // nearest sensor
var yy = X.map(example_input), linear = interp(X, sx, ss), steps = X.map(nearest_value);
var eLinear = X.map(function (x, i) { return abs(yy[i] - linear[i]); });
var eStep = X.map(function (x, i) { return abs(yy[i] - steps[i]); });
panel[0].draw({ylim: [-0.15, 1.25], series: [
{type: 'line', x: X, y: yy, color: 's1', width: 2.8, name: 'true f(x)'},
{type: 'line', x: X, y: linear, color: 's2', dash: '6 4', name: 'straight-line fₘ'},
{type: 'line', x: X, y: steps, color: 's3', width: 1.7, name: 'staircase fₘ'},
{type: 'dots', x: sx, y: ss, color: 'ink', r: 3.5, name: 'sensor values'}]});
panel[1].draw({ylim: [-0.04, 1.1], series: [
{type: 'line', x: X, y: eLinear, color: 's2', name: 'linear: ' + fmt(eLinear[argmax(eLinear)], 3)},
{type: 'line', x: X, y: eStep, color: 's3', name: 'stairs: ' + fmt(eStep[argmax(eStep)], 3)}]});
readout({'Fixed spike': 'A fixed spike can hide on a coarse grid but becomes visible as sensors get closer.',
'Step (jump)': 'The jump is discontinuous. Without its location, grid sampling cannot guarantee small uniform error.',
'Smooth wave': 'A smooth input can be reconstructed more accurately as sensor spacing shrinks.'}[P.family]);
''')
At its narrowest, the two inputs $f(x)=0$ and $f(x)=\text{spike}(x)$ can produce exactly the same sensor vector, yet $\|f-0\|_\infty=1$. No predictor receiving only the sensor vector can distinguish those functions. Adding any fixed number of sensors does not help if I am allowed to place an arbitrarily narrow spike between them.
This tells us what we need: a family of allowable inputs whose behavior between measurements is controlled. That leads to compactness.
2. Compactness: what do we restrict?¶
Start with an ordinary interval. Square brackets include endpoints; parentheses exclude them:
$$ [0,1]:\ 0\le x\le1,\qquad (0,1):\ 0<x<1. $$
Both are bounded, but only $[0,1]$ is compact. Follow the sequence $1/2,1/4,1/8,\ldots$. Its destination is zero. That destination belongs to the first interval but has been left out of the second. Compactness prevents a sequence from converging out of the set. In ordinary Euclidean space, closed and bounded is equivalent to compact.
interactive(
controls=[dict(key='n', label='Halvings n', min=1, max=18, step=1, value=1, play='loop', duration=9)],
panels=[dict(axes=False, height=150)],
draw='''
var pos = Math.pow(2, -P.n), trail = [], series = [];
for (var k = 1; k < P.n; k++) trail.push(Math.pow(2, -k));
[{y: 1, closed: true, text: '[0,1] includes zero'}, {y: 0, closed: false, text: '(0,1) excludes zero'}].forEach(function (r) {
series.push({type: 'line', x: [0, 1], y: [r.y, r.y], color: 's1', width: 4});
series.push({type: 'dots', x: [0, 1], y: [r.y, r.y], color: 's1', r: 6, hollow: !r.closed});
series.push({type: 'dots', x: trail, y: trail.map(function () { return r.y; }), color: 's2', r: 3.5, opacity: 0.35});
series.push({type: 'dots', x: [pos], y: [r.y], color: 's2', r: 7});
series.push({type: 'text', x: 1.055, y: r.y, text: r.text, size: 13});
});
panel[0].draw({xlim: [-0.08, 1.65], ylim: [-0.5, 1.5], series: series});
readout('2<sup>−' + P.n + '</sup> = ' + pos.toFixed(5) + ' approaches 0 as n grows');
''')
In a function space, the points are entire curves, and distance means the maximum separation between curves. A set of curves is compact if every sequence of curves contains a subsequence that converges uniformly to a curve in the set.
This is stricter than simply bounding their height. Consider the functions
$$ f_n(x)=\sin(2\pi n x). $$
Each stays between $-1$ and $1$, but higher $n$ creates finer and finer oscillations. The curves cannot settle toward one continuous curve uniformly. Try increasing $n$: the height is bounded, but the horizontal detail is not.
interactive(
controls=[dict(key='n', label='Frequency n', min=1, max=25, step=1, value=2, play='loop', duration=8)],
panels=[dict(title='fₙ(x) = sin(2πnx): bounded height, unbounded frequency', xlabel='x', ylabel='f_n(x)', height=250)],
draw='''
var xx = linspace(0, 1, 2400), y = xx.map(function (x) { return sin(2 * PI * P.n * x); });
panel[0].draw({ylim: [-1.2, 1.2], series: [
{type: 'hline', y: 1, color: 's2', width: 1.2, dash: '6 4'},
{type: 'hline', y: -1, color: 's2', width: 1.2, dash: '6 4'},
{type: 'line', x: xx, y: y, color: 's1', width: 2.5, name: 'fₙ(x)', legend: false}]});
''')
The missing ingredient is that one horizontal scale must work for every curve. Given a small height tolerance $\varepsilon$, there must be a distance $\delta$ so that any two points closer than $\delta$ have function values within $\varepsilon$, regardless of which curve we picked:
$$ |x-z|<\delta\quad\Longrightarrow\quad |f(x)-f(z)|<\varepsilon\qquad\text{for all }f\in K. $$
This is equicontinuity. Together with a uniform height bound, it characterizes relative compactness in $C([0,1])$ (Arzelà–Ascoli); including the uniform limits makes the family compact. The oscillating sines fail equicontinuity: every horizontal scale eventually misses a rapid oscillation.
A convenient concrete family is
$$ \boxed{K=\{f:\ |f(x)|\le M,\quad |f(x)-f(z)|\le L|x-z|\}}. $$
The shared slope bound $L$ rules out arbitrarily narrow spikes. The word shared is the point: a different bound for each curve is not enough.
From compactness to a sensor error bound¶
Put evenly spaced sensors a distance $h=1/(m-1)$ apart, and connect measurements with straight lines to get $f_m$. Between adjacent sensors, write $x=x_i+th$, $0\le t\le1$. Then
$$ \begin{aligned} f_m(x)&=(1-t)f(x_i)+tf(x_{i+1}),\\ |f(x)-f_m(x)|&\le 2Lht(1-t)\le \frac{Lh}{2}. \end{aligned} $$
We have bounded the error for every location and every input function:
$$ \boxed{\sup_{f\in K}\|f-f_m\|_\infty\le\frac{Lh}{2}.} $$
The outer $\sup_{f\in K}$ means take the worst function in the whole family. To guarantee error below $\varepsilon$, choose enough sensors that $Lh/2<\varepsilon$. Vary the number of sensors and the input phase below. The measured error changes, but the bound always holds.
interactive(
controls=[
dict(key='m', label='Sensors m', min=5, max=81, step=2, value=13),
dict(key='phase', label='Curve (phase)', min=0, max=6.28, step=0.01, value=0.8, play='loop', duration=10),
dict(key='eps', label='Target ε', min=0.04, max=0.30, step=0.02, value=0.16),
],
panels=[
dict(title='Input function and linear interpolation', xlabel='position x'),
dict(title='Actual error and guaranteed bound', xlabel='position x'),
],
draw='''
var L = 0.30 * 2 * PI + 0.12 * 6 * PI; // one bound for every phase in this example
function smooth_family(x) { return 0.30 * sin(2 * PI * x) + 0.12 * sin(6 * PI * x + P.phase); }
var xs = linspace(0, 1, P.m), sample = xs.map(smooth_family), X = grid(1601, xs);
var truth = X.map(smooth_family), recon = interp(X, xs, sample);
var e = X.map(function (x, i) { return abs(truth[i] - recon[i]); }), err = e[argmax(e)], bound = L / (2 * (P.m - 1));
panel[0].draw({ylim: [-0.55, 0.55], series: [
{type: 'line', x: X, y: truth, color: 's1', width: 2.7, name: 'f(x)'},
{type: 'line', x: X, y: recon, color: 's2', dash: '6 4', name: 'fₘ(x)'},
{type: 'dots', x: xs, y: sample, color: 's2', r: 3.5}]});
panel[1].draw({ylim: [0, max(0.12, bound, P.eps) * 1.25], series: [
{type: 'line', x: X, y: e, color: 's3', width: 2.4, name: 'actual |f − fₘ|'},
{type: 'hline', y: bound, color: 'ink', width: 1.7, dash: '0.5 4', name: 'Lh/2 (guarantee)'},
{type: 'hline', y: P.eps, color: 's2', width: 1.5, dash: '6 4', name: 'ε (target)'}]});
readout('Observed max error: ' + fmt(err, 4) + ' Universal bound Lh/2: ' + fmt(bound, 4) + ' ' +
(bound < P.eps ? 'guaranteed below target' : 'not yet guaranteed by this bound'));
''')
Compactness does not mean that a function has finitely many degrees of freedom. It means that, at any prescribed accuracy, finitely many measurements can approximate every function in the family. As we demand smaller error, we may need more sensors.
Now we can return to the mathematical problem. We know how to make an input function accessible to a network. But does a small error in the input cause a small error in the solution?
3. What does the operator actually do?¶
Consider the same boundary-value problem throughout:
$$ \boxed{ \begin{aligned} \mathcal L u(x)&=f(x),\qquad\mathcal L=-\frac{d^2}{dx^2},\quad 0<x<1,\\ u(0)&=u(1)=0, \end{aligned}} \qquad \boxed{u=\mathcal G(f)=\mathcal L^{-1}f.} $$
Each input forcing $f$ produces an output solution $u$. The operator $\mathcal G$ is the entire rule connecting the two families, not one particular solution.
To see it working, take two forcing modes:
$$ f(x)=a\sin(\pi x)+b\sin(3\pi x). $$
Change $a$ and $b$. The two plots use different vertical scales so the solution remains visible.
interactive(
controls=[
dict(key='a', label='a · slow mode', min=-2, max=2, step=0.05, value=1.2),
dict(key='b', label='b · fast mode', min=-2, max=2, step=0.05, value=1.0),
],
panels=[
dict(title='Input f(x): forcing', xlabel='position x', ylabel='forcing amplitude'),
dict(title='Output u(x) = 𝒢(f)(x): solution', xlabel='position x', ylabel='solution amplitude'),
],
draw='''
function forcing(x) { return P.a * sin(PI * x) + P.b * sin(3 * PI * x); }
function solution(x) { return P.a * sin(PI * x) / (PI * PI) + P.b * sin(3 * PI * x) / (9 * PI * PI); }
var X = linspace(0, 1, 801);
panel[0].draw({ylim: [-3.3, 3.3], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: X.map(forcing), color: 's2', width: 2.8, name: 'f(x)', legend: false}]});
panel[1].draw({ylim: [-0.23, 0.23], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: X.map(solution), color: 's1', width: 2.8, name: 'u(x)', legend: false}]});
''')
A tiny ripple in the solution can become a large ripple in the forcing¶
The previous figure showed the inverse map, forcing $f$ $\longmapsto$ solution $u$. Look at the same equation in the other direction: $\mathcal L u=f$. Differentiation acts on every sine mode according to
$$ -\frac{d^2}{dx^2}\sin(n\pi x)=(n\pi)^2\sin(n\pi x). $$
Take a smooth-looking solution with a tiny high-frequency ripple:
$$ \boxed{u(x)=\sin(\pi x)+0.01\sin(10\pi x).} $$
Differentiating twice gives
$$ \boxed{f(x)=-u''(x)=\pi^2\sin(\pi x)+\pi^2\sin(10\pi x).} $$
The high-frequency term is 1% of the main term in $u$, but 100% of its amplitude in $f$. Change the frequency and ripple amplitude below. The lower panels isolate the ripple, so you can see it even when it is almost invisible in the full solution.
interactive(
controls=[
dict(key='n', label='Frequency n', min=2, max=20, step=1, value=10, play='bounce', duration=8),
dict(key='ripple', label='Ripple α', min=0, max=0.04, step=0.0005, value=0.01),
],
panels=[
dict(title='Solution u: the ripple is almost invisible', xlabel='x', ylabel='solution', height=250, basis='calc(50% - 11px)', left=64),
dict(title='Forcing f = −u″: ripple amplified', xlabel='x', ylabel='forcing', height=250, basis='calc(50% - 11px)', left=64),
dict(title='Isolated ripple in u', xlabel='x', ylabel='u_high', height=170, basis='calc(50% - 11px)', left=64),
dict(title='Isolated ripple in f', xlabel='x', ylabel='f_high', height=170, basis='calc(50% - 11px)', left=64),
],
draw='''
var n = P.n, ripple = P.ripple, k2 = (n * PI) * (n * PI), X = linspace(0, 1, 2001);
var low = X.map(function (x) { return sin(PI * x); });
var high = X.map(function (x) { return ripple * sin(n * PI * x); });
var uu = low.map(function (v, i) { return v + high[i]; });
var flow = low.map(function (v) { return PI * PI * v; });
var fhigh = high.map(function (v) { return k2 * v; });
var ff = flow.map(function (v, i) { return v + fhigh[i]; });
var F = 1.15 * (PI * PI + ripple * k2), cu = max(0.004, 1.17 * ripple), cf = max(0.25, 1.17 * ripple * k2);
panel[0].draw({ylim: [-0.12, 1.15], series: [
{type: 'line', x: X, y: uu, color: 's1', width: 2.8, name: 'u(x)'},
{type: 'line', x: X, y: low, color: 'ink2', width: 1.5, dash: '5 4', name: 'slow mode only'}]});
panel[1].draw({ylim: [-F, F], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: ff, color: 's2', width: 2.7, name: 'f(x)'},
{type: 'line', x: X, y: flow, color: 'ink2', width: 1.5, dash: '5 4', name: 'slow mode only'}]});
panel[2].draw({ylim: [-cu, cu], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: high, color: 's1', width: 2.3, name: 'u_high', legend: false, digits: 5}]});
panel[3].draw({ylim: [-cf, cf], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: fhigh, color: 's2', width: 2.3, name: 'f_high', legend: false}]});
readout('High-frequency amplitude: in u: ' + fmt(ripple, 3) + '; in f: ' + fmt(ripple * k2, 3) +
' (relative amplification: n² = ' + n * n + ')<br>' +
'Slow-mode amplitude: u = 1, f = π² = ' + fmt(PI * PI, 3) + '. High/low amplitude: u = ' + fmt(100 * ripple, 1) +
'%, f = ' + fmt(100 * ripple * n * n, 1) + '%.');
''')
Both curves here are infinitely differentiable: a rapidly oscillating sine is not less differentiable than a slowly oscillating one. What changes is the size of the derivatives and the weight of high-frequency modes. In Sobolev regularity, solving this Poisson equation can gain two derivatives under the usual conditions. That is why the solution operator $\mathcal G=\mathcal L^{-1}$ smooths while the differential operator $\mathcal L$ amplifies high frequencies.
4. The matrix hidden inside the operator¶
Why does changing the fast sine mode barely affect the solution? Differentiate any sine mode twice:
$$ -\frac{d^2}{dx^2}\sin(k\pi x)=(k\pi)^2\sin(k\pi x). $$
The equation $\mathcal Lu=f$ therefore divides the amplitude of mode $k$ by $(k\pi)^2$:
$$ \boxed{u(x)=\frac{a}{\pi^2}\sin(\pi x)+\frac{b}{9\pi^2}\sin(3\pi x).} $$
Now forget the curves and look at their coefficients. The entire solution operator, on this two-mode subspace, is a matrix:
$$ \boxed{ \underbrace{\begin{bmatrix}a/\pi^2\\b/(9\pi^2)\end{bmatrix}}_{\text{output coefficients}} =\underbrace{\begin{bmatrix}1/\pi^2&0\\0&1/(9\pi^2)\end{bmatrix}}_{\mathcal G\text{ in the sine basis}} \underbrace{\begin{bmatrix}a\\b\end{bmatrix}}_{\text{input coefficients}}. } $$
The third sine mode is attenuated nine times more than the first. Here the map between function spaces became a simple diagonal scaling of two numbers. For a nonlinear PDE, the coefficient map will generally be nonlinear; nevertheless, its input and output can still be finite coefficient vectors once we choose finite approximations.
This is the geometric insight behind operator learning: choose coordinates for functions, then learn how the coordinates transform.
5. Continuity: will input errors ruin the solution?¶
Imagine perturbing the forcing by a sine wave:
$$ \Delta f(x)=\delta\sin(n\pi x). $$
Our solution operator transforms that perturbation into
$$ \boxed{\Delta u(x)=\frac{\delta}{(n\pi)^2}\sin(n\pi x).} $$
Change $\delta$ or the frequency. Watch the solution respond smoothly. Higher frequency perturbations have smaller effect on the output, even for the same input amplitude.
interactive(
controls=[
dict(key='delta', label='Change δ', min=0, max=1, step=0.01, value=0.8),
dict(key='n', label='Frequency n', min=1, max=11, step=1, value=3, play='bounce', duration=7),
],
panels=[
dict(title='Change in forcing Δf', xlabel='position x', ylabel='input change'),
dict(title='Change in solution Δu', xlabel='position x', ylabel='output change'),
],
draw='''
var X = linspace(0, 1, 1201), c = P.delta / ((P.n * PI) * (P.n * PI));
var du = X.map(function (x) { return P.delta * sin(P.n * PI * x); }); // change in the forcing
var dv = X.map(function (x) { return c * sin(P.n * PI * x); }); // change in the solution
var U = 1.15 * max(c, 1e-6); // the axis follows the amplitude of the response
panel[0].draw({ylim: [-1.15, 1.15], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: du, color: 's2', width: 2.6, name: 'Δf(x)', legend: false}]});
panel[1].draw({ylim: [-U, U], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'line', x: X, y: dv, color: 's1', width: 2.6, name: 'Δu(x)', legend: false, digits: 6}]});
readout('‖Δf‖∞ = ' + fmt(P.delta, 4) + ' ‖Δu‖∞ = ' + fmt(c, 5) + ' attenuation = 1/(' + P.n + 'π)²');
''')
For the operator to be continuous, nearby input functions must produce nearby output functions. We measure nearness with the same maximum-gap distance:
$$ \|f-g\|_\infty\text{ small}\quad\Longrightarrow\quad \|\mathcal G(f)-\mathcal G(g)\|_\infty\text{ small}. $$
For this Poisson equation we can do better than saying "small": the operator obeys the quantitative bound
$$ \boxed{\|\mathcal G(f)-\mathcal G(g)\|_\infty\le \frac18\|f-g\|_\infty.} $$
Where does $1/8$ come from? If a forcing is at most one in magnitude everywhere, the largest possible solution amplitude is bounded by the solution for the constant forcing $f=1$:
$$ -u''=1\quad\Longrightarrow\quad u(x)=\frac{x(1-x)}2, \qquad \max_x u(x)=\frac18. $$
More generally, the positive Green's function for these Dirichlet boundary conditions gives the same bound for positive or negative forcing. Thus an input error no greater than $\eta$ everywhere produces an output error no greater than $\eta/8$ everywhere. The operator cannot amplify a small input error into a large output error.
Now combine continuity with our finite sensors. If $f_m$ is the piecewise-linear sensor reconstruction and $K$ has slope bound $L$, then
$$ \boxed{ \|\mathcal G(f)-\mathcal G(f_m)\|_\infty \le \frac18\|f-f_m\|_\infty \le \frac{Lh}{16}. } $$
We have connected a finite list of input measurements to an error bound on the entire output function. That is the mathematical reason finite measurements can work for a continuous solution operator.
Continuity alone does not promise this particular constant; for a general continuous operator on a compact input family, the guarantee is instead: a sufficiently small input error gives a uniformly small output error. If the PDE lacks a unique or stable solution, there may be no such well-behaved solution map to learn.
6. One last reduction: from an output function to coefficients¶
The network still needs to produce a whole output curve. We can represent that curve using finitely many basis functions:
$$ \boxed{\mathcal G(f)(y)\approx\sum_{j=1}^{p} c_j(f)\,\phi_j(y).} $$
Here $\phi_j(y)$ describes an output shape and $c_j(f)$ determines how much of it the current forcing produces. We have met this idea already. For our two-mode example,
$$ \underbrace{\mathcal G(f)(y)}_{\text{solution value at }y} =\underbrace{\begin{bmatrix}a/\pi^2&b/(9\pi^2)\end{bmatrix}}_{\text{output coefficients}} \underbrace{\begin{bmatrix}\sin(\pi y)\\\sin(3\pi y)\end{bmatrix}}_{\text{basis values at }y}. $$
Slide the query point $y$. The coefficients stay fixed for that input function; the basis values change with the query. The sum gives the solution at that point. The separation between coefficients and where to evaluate them is the branch–trunk idea.
interactive(
controls=[
dict(key='y', label='Query y', min=0, max=1, step=0.005, value=0.35, play='bounce', duration=8),
dict(key='a', label='a · slow mode', min=-2, max=2, step=0.05, value=1.2),
dict(key='b', label='b · fast mode', min=-2, max=2, step=0.05, value=1.0),
],
panels=[
dict(title='The output function: u(y)', xlabel='position y', ylabel='solution value', xname='y'),
dict(title='Dot product at the chosen location', xlabel='contribution to u(y)', crosshair=False, left=72),
],
draw='''
var coeff = [P.a / (PI * PI), P.b / (9 * PI * PI)], y = P.y;
var basis = [sin(PI * y), sin(3 * PI * y)];
var terms = [coeff[0] * basis[0], coeff[1] * basis[1]], val = terms[0] + terms[1];
var X = linspace(0, 1, 801), u = X.map(function (x) { return coeff[0] * sin(PI * x) + coeff[1] * sin(3 * PI * x); });
panel[0].draw({ylim: [-0.23, 0.23], series: [
{type: 'hline', y: 0, color: 'muted', width: 1},
{type: 'vline', x: y, color: 's3', width: 1.2, dash: '5 4'},
{type: 'line', x: X, y: u, color: 's1', width: 2.6, name: 'u(y)', legend: false},
{type: 'dots', x: [y], y: [val], color: 's3', r: 6}]});
panel[1].draw({xlim: [-0.23, 0.23], categories: ['c_1φ_1(y)', 'c_3φ_3(y)'], series: [
{type: 'hbars', values: terms, colors: ['s1', 's3'], labels: false},
{type: 'text', x: -0.22, y: 0.96, text: 'coefficient vector = (' + fmt(coeff[0], 3) + ', ' + fmt(coeff[1], 3) + ')'}]});
readout('u(' + fmt(y, 2) + ') = (' + fmt(coeff[0], 4) + ')(' + fmt(basis[0], 4) + ') + (' + fmt(coeff[1], 4) + ')(' +
fmt(basis[1], 4) + ') = ' + fmt(val, 5));
''')
7. DeepONet's universal approximation theorem, without a leap¶
The branch network receives a finite sensor vector $\mathbf s(f)=[f(x_1),\ldots,f(x_m)]$ and predicts output coefficients. The trunk network takes the query location $y$ and returns basis values. Their dot product gives the predicted solution:
$$ \boxed{ \mathcal G_\theta(f)(y) =\sum_{j=1}^{p} \underbrace{b_{\theta,j}\big(f(x_1),\ldots,f(x_m)\big)}_{\text{branch: coefficients}} \underbrace{t_{\theta,j}(y)}_{\text{trunk: basis values}}. } $$
Our figure used coefficients calculated from the known equation, not coefficients learned from training. DeepONet learns that coefficient rule and can learn the basis shapes too.
Why does the theorem work? Suppose $K$ is a compact family of continuous input functions, and $\mathcal G:K\to C(Y)$ is a continuous operator on a compact output domain $Y$.
- Compact inputs: a sufficiently fine finite sensor grid captures all functions in $K$ to any prescribed accuracy. The spike example showed why the restriction matters.
- Continuity of the map: input approximation can be made accurate enough that output errors are uniformly small. The Poisson example gave a concrete bound.
- Compact outputs: because $\mathcal G$ is continuous and $K$ is compact, the output family $\mathcal G(K)$ is compact too. Its curves can therefore be approximated by finitely many output basis functions—for instance, piecewise-linear interpolation functions.
- Ordinary neural-network approximation: what remains is a continuous finite-dimensional map from sampled inputs to output coefficients. Ordinary UAT provides neural networks that approximate it. For a general operator, a rigorous construction uses a continuous finite-dimensional surrogate (rather than assuming reconstruction stays inside $K$).
Taken together, we can make the total error smaller than any $\varepsilon>0$:
$$ \boxed{ \sup_{f\in K}\sup_{y\in Y} \left|\mathcal G(f)(y)-\mathcal G_\theta(f)(y)\right|<\varepsilon. } $$
What the theorem does not say: how many sensors $m$, output modes $p$, or training examples are needed. Compactness guarantees approximation exists, not that it is cheap. A family of moving sharp fronts can be compact yet require many fixed output basis functions.
The central question in practical operator learning is therefore not merely can we approximate the operator? It is can we find coordinates in which the map from one function to another is simple? Our Poisson operator happened to be diagonal in the sine basis. For more complicated physics, discovering good coordinates is much of the research problem.