Module 4 of 14 · one week
Two other ways to put a network in front of a PDE, and what each one runs into. DeepONet factorizes the operator into a learned basis and learned coefficients, and pays for it with a hard information limit that no amount of training removes. A PINN fits one solution by penalising the equation, and fails for a reason about optimization, not capacity. The hybrids are what the field built after each pure form met its limit.
Module 0 separated three objectives. This module supplies the architectures behind the two that Modules 2 and 3 did not cover. DeepONet is the older answer to the surrogate objective and still the default in the DeepXDE ecosystem, so any literature search will return it. The PINN is the constraint objective in its purest form. And the hybrids matter because they are where the two lines cross, which is where most current work sits.
There is also a specific claim to settle. The earlier episode says at 00:11:58 that "PINNs don't work everywhere", which is true and, in the episode, unexplained. Krishnapriyan et al. explain it with numbers and a diagnosis, and that diagnosis is the load-bearing fact for everything in Module 11. Getting it right here saves a mistaken design later.
Chen and Chen proved in 1995 that a network with one hidden layer can approximate any continuous operator between function spaces on a compact set. That result sat unused for two decades because it says nothing about how to build such a network. DeepONet (Lu, Jin, Karniadakis, 2019) is the first architecture designed around it, and the design reads the theorem's proof almost literally.
The input function \(u\) is read at \(m\) fixed sensor locations \(x_1,\dots,x_m\). A branch network turns the vector of readings into \(p\) coefficients. A trunk network turns a query point \(y\) into \(p\) values. The prediction is their dot product (Equation 2 of the paper):
\[G_\theta(u)(y) = \sum_{k=1}^p b_k\big(u(x_1),\dots,u(x_m)\big)\,t_k(y) + b_0 .\]
That is a basis expansion whose basis is learned. The trunk supplies functions \(t_k\) of the output variable; the branch supplies the coefficients of the answer in that basis. The dependence on the input is nonlinear, because the branch is a network, but the reconstruction from coefficients to output function is linear, and that asymmetry is the architecture's defining property. Two variants exist: a "stacked" version with one branch per coefficient, and an "unstacked" version with a single branch producing all \(p\) at once, which has fewer parameters and generalizes better (Section 4.1, research_notes/web-papers.md:50-58).
The paper's numbers are the ones to keep, because they set expectations for the whole family (:60-80). On the antiderivative operator, 100 sensors and 10,000 training functions give a test error 10 to 20 times lower than a fully connected network's \(10^{-4}\) mean squared error. On the forced pendulum, error falls like \(4.6^{-m}\) in the sensor count \(m\) until it saturates near 10 sensors, and a longer horizon needs more: about 25 sensors at \(T=5\). Against dataset size the error falls exponentially below roughly \(10^4\) samples and algebraically above it, with generalization error scaling like \(N^{-1}\), which the paper flags as surprising against the classical \(N^{-1/2}\).
The two families are not rivals with separate theories; one contains the other. JMLR Proposition 5, on page 25, states that "a neural operator with a point-wise parameterized first kernel and discretized integral operators yields a DeepONet", and the proof constructs it: choose \(W_0=0\), make the first kernel separable as \(\kappa^{(0)}_j(y,z) = \mathbf 1(y)w_j(z)\), replace the inner integral by its Monte Carlo estimate on the sensor points, and the layer collapses to \(\sum_k G_k(\tilde v)\varphi_k(x)\), which is Equation 2 above.
That construction also shows the defect. The branch weights \(\tilde w_j\) are tied to the fixed sensor points, so refining the mesh changes what the parameters mean and the parameter count has to grow. This is why DeepONet scores "no" on discretization invariance in Module 1's Table 1: not from carelessness, but because the Monte Carlo step froze the discretization into the weights. The paper's own repair, which it calls DeepONet-Operator (Equation 33), replaces the finite inner product \(\langle\tilde w_j,\tilde a\rangle\) with the function-space inner product \(\langle w_j,a\rangle\), and the result is discretization invariant (research_notes/pdf-fno-jmlr.md:430-436). Module 3's LNO is that same structure arrived at from the cost side.
The deeper limit, which the repair does not fix. JMLR page 26 points out that any architecture with a fixed trunk basis is a linear approximation method: whatever the branch computes, the output is confined to the span of \(t_1,\dots,t_p\). The best possible error of a linear method over a set of output functions is the Kolmogorov \(n\)-width of that set, and for advection-dominated maps the \(n\)-width decays slowly, so no training improves it (:437-445). Module 3 gave the mirror image of this from the spectrum side: low-rank works when the operator smooths. Here it is the output set rather than the operator, and the same conclusion follows.
Before any of the approximation theory, there is a limit that comes from the sampling alone and applies to every model that reads its input at fixed points, DeepONet included. The proof is two lines and it puts a floor under a whole design family.
Claim. Let \(G_\theta\) be any model whose only access to the input is the values at a sensor set \(S=\{x_1,\dots,x_m\}\). Suppose \(u_1\) and \(u_2\) agree at every point of \(S\), while \(\|G^\dagger(u_1)-G^\dagger(u_2)\|_Y = \delta\). Then at least one of the two prediction errors is at least \(\delta/2\).
Proof. The two inputs produce identical sensor vectors, so \(G_\theta(u_1)=G_\theta(u_2)\); call that common output \(w\). By the triangle inequality, \[\delta = \|G^\dagger(u_1)-G^\dagger(u_2)\| \le \|G^\dagger(u_1)-w\| + \|w-G^\dagger(u_2)\|,\] so the two errors sum to at least \(\delta\) and the larger is at least \(\delta/2\). \(\square\)
The claim is only interesting if such a pair \(u_1, u_2\) exists, so construct one with numbers. Put \(m=16\) sensors at \(x_j = j/16\) on \((0,1)\), and take the antiderivative operator \(G^\dagger u(x) = \int_0^x u\). Define a bump living entirely in the gap between the sensors at 0.5 and 0.5625:
\[\varphi(x) = \begin{cases}\sin^2\big(16\pi(x-0.5)\big), & 0.5 \le x \le 0.5625,\\ 0, & \text{otherwise.}\end{cases}\]
It is continuous with a continuous derivative, and it vanishes at both ends of the gap because \(\sin(16\pi \cdot 0.0625) = \sin\pi = 0\), hence at every sensor. So \(u\) and \(u+\varphi\) are the same input as far as the model can tell, for any \(u\) at all. Their images are not the same. Since \(\int \sin^2 = \tfrac12\) times the interval length, \[\int_0^1 \varphi = \tfrac12 \cdot 0.0625 = 0.03125,\] so \(G^\dagger(u+\varphi) - G^\dagger(u)\) is zero up to \(x=0.5\), rises across the gap, and equals 0.03125 for every \(x \ge 0.5625\). Its \(L^2(0,1)\) norm is at least \(0.03125\sqrt{1-0.5625} = 0.0207\). By the claim, one of the two predictions is wrong by at least 0.0103 in \(L^2\), for every set of weights, after any amount of training.
Two features of that number decide how seriously to take it. It scales like \(1/m\) for a bump of fixed height, since the gap width does, so doubling the sensors halves the floor. And it assumed the bump has height one, which is a statement about the input distribution, not about the sensors: if the inputs are drawn from a measure whose high-frequency content is small, a bump that fits between sensors is small too, and the floor drops with it.
Lanthaler, Mishra and Karniadakis turn that observation into the exact bound. The encoding error of \(m\) sensors is bounded below by the tail of the input measure's covariance spectrum, \(\sqrt{\sum_{k>m}\lambda_k} \le \hat{\mathcal E}_{\mathcal E}\) (Theorem 3.8, research_notes/web-papers.md:337-343). Compare that with Module 0's Fact 2 calculation on Lab 01's inputs, where \(\lambda_k = 625/(k^2+25)^2\): the tail beyond \(m\) behaves like \(m^{-3}\), so the sensor floor falls fast for that measure and would not fall at all for white noise.
One more prerequisite is easy to skip. Point evaluation \(u \mapsto u(x_j)\) is not defined on \(L^2\) at all, by Module 0's Fact 1, so the sensor reading is meaningless unless the input space is one where evaluation is bounded: \(H^s(D)\) with \(s > d/2\), or \(C(\bar D)\). Every DeepONet result silently assumes this, and it is the same assumption as route 1 of Module 0's three routes.
With the floor established, the rest of the theory splits the remaining error into three named pieces (Lanthaler, Mishra, Karniadakis, Theorem 3.3, :331-336). The encoding error asks how well \(m\) sensor values capture the input; the bound above is its lower limit. The approximation error asks how well the finite branch and trunk networks approximate the induced map between finite-dimensional spaces. The reconstruction error asks how well \(p\) trunk functions span the outputs actually needed, and it too is bounded below by a covariance tail, this time of the output measure, \(\sqrt{\sum_{k>p}\lambda_k}\) (Theorem 3.6). For Lipschitz operators into \(H^s(\mathbb T^n)\), reconstruction costs at most \(Cp^{-s/n}\) with a trunk of size at most \(Cp(1+\log^2p)\) (Theorem 3.7).
The middle term is where the curse of dimensionality would normally live. A generic Lipschitz map \(\mathbb R^m \to \mathbb R\) needs a ReLU network of size about \(\varepsilon^{-m/2}\) (Remark 3.4, after Yarotsky), which at \(m = 100\) sensors is hopeless. The paper's four PDE examples nonetheless achieve algebraic rates, and the reason is that their solution operators are smooth in a specific sense and not merely Lipschitz (:348-363). That escape clause is the same shape as every positive result in this course, and Module 5 puts it next to the FNO version.
Raissi, Perdikaris and Karniadakis take a different route entirely (research_notes/web-papers.md:105-175). A network \(u_\theta(t,x)\) is trained so that its own automatic derivatives satisfy \(u_t + \mathcal N[u] = 0\). The loss has two parts: squared error against \(N_u\) known initial and boundary values, plus the squared PDE residual evaluated at \(N_f\) collocation points where nothing is known. L-BFGS optimizes both. Part II makes the equation's unknown coefficients trainable parameters of the same loss, so the method recovers them from scattered noisy observations.
| Example (Part I) | Setup | Relative L2 |
|---|---|---|
| Burgers, \(u_t + uu_x - (0.01/\pi)u_{xx}=0\) | continuous time, \(N_u=100\), \(N_f=10{,}000\), 9 layers × 20 | 6.7e-4 |
| Burgers, discrete time | \(q=500\) RK stages, one step \(\Delta t=0.8\), \(N_n=250\), 4 × 50 | 8.2e-4 |
| Nonlinear Schrödinger | \(N_0=N_b=50\), \(N_f=20{,}000\), 5 × 100 tanh | 1.97e-3 |
| Allen-Cahn, discrete time | \(q=100\), \(\Delta t=0.8\), \(N_n=200\), 4 × 200 | 6.99e-3 |
One number in that table deserves attention because it isolates the error source. The discrete-time model embeds a \(q\)-stage implicit Runge-Kutta scheme, whose truncation error is \(O(\Delta t^{2q})\). At \(q=500\) and \(\Delta t=0.8\) that is \(0.8^{1000} \approx 10^{-97}\), which is zero for every practical purpose. So the \(8.2\times10^{-4}\) that the method actually achieves is entirely network approximation and optimizer error, with no contribution from the time integrator (:149-152). The authors are candid about the standing of the whole approach, writing that PINNs are "not... replacements of classical numerical methods" (:159-162).
Part II's inverse results carry a lesson of their own. On Burgers with no noise, the recovered coefficients are off by 0.096 percent for \(\lambda_1\) and 0.469 percent for \(\lambda_2\); at 10 percent noise, \(\lambda_2\)'s error rises to 6.39 percent while \(\lambda_1\) stays much better (:200-228). The asymmetry is not an accident. \(\lambda_2\) multiplies the diffusion term, and additive noise at every point looks locally like diffusion, so the two are hard to tell apart. An inverse problem's conditioning depends on which coefficient you are after.
The structural fact, which no tuning changes. A PINN produces one function for one problem instance. A new initial condition means a new optimization from scratch. That is why PDEBench has to train ten separate PINNs per dataset only to report error bars, cannot run a PINN on Darcy flow at all because a coefficient-to-solution map is not a single solution to fit, and skips PINNs on 2D and 3D compressible flow for memory reasons (research_notes/pdf-mgn-pdebench.md:227, :243). Comparing a PINN with an operator on equal terms is therefore hard by construction, and any benchmark table that puts them in the same column is comparing two different objectives.
Krishnapriyan, Gholami, Zhe, Kirby and Mahoney ran the experiment that settles the mechanism (NeurIPS 2021, research_notes/web-papers.md:585-653). Their test problem is convection, \(u_t + \beta u_x = 0\) on \([0,2\pi]\) with \(u(x,0)=\sin x\) and periodic boundaries, chosen because the exact solution \(\sin(x-\beta t)\) is known for every \(\beta\) and because \(\beta\) turns difficulty into a dial.
The result is a wall rather than a slope. A standard PINN is accurate at small \(\beta\) and reaches roughly 100 percent relative error once \(\beta\) passes about 10. A reaction problem, \(u_t - \rho u(1-u)=0\), shows about 98 percent error at \(\rho=5\). Reaction-diffusion with \(\rho=5\) and \(\nu\) from 2 to 6 lands between 50 and 96 percent. Crucially, the same network handed the correct solution as an initialisation holds it, so the capacity was always there. What changes with \(\beta\) is the loss landscape: viewed along the top two Hessian eigenvectors it is smooth at \(\beta=1\) and "complex and non-symmetric" at high \(\beta\), and the condition number grows like \((\beta N)^2\) for convection and \((\nu N^2)^2\) for diffusion. The authors' diagnosis, in their words: "The failure is due to optimization difficulties associated with the PINN's soft PDE constraint."
Both of their repairs follow from that diagnosis rather than from architecture, and both work.
| Fix | Mechanism | Case | Before | After | Factor |
|---|---|---|---|---|---|
| Curriculum regularization | train at \(\beta=1\), then warm-start progressively harder \(\beta\) | convection, \(\beta=30\) | 8.97e-1 | 2.02e-2 | 44× |
| Sequence to sequence | march in time windows \(\Delta t\), same collocation budget | reaction-diffusion, \(\nu=\rho=5\), \(\Delta t=0.1\) | 9.35e-1 | 2.39e-2 | 39× |
Those two factors are the "1 to 2 orders of magnitude" of the abstract, and neither one adds a parameter. Both change only the sequence of optimization problems the network is asked to solve. Module 11 takes the curriculum idea as the template for using physics as a data engine instead of a loss term.
Each pure form has a complementary weakness. A PINN gets a dense training signal from the equation but solves one instance at a time. An operator amortizes across instances but needs paired data, which means running the expensive solver you were trying to avoid. The obvious move is to put the residual into an operator's loss, and two papers do it from the two directions.
PINO (Li et al. 2021, arXiv 2111.03794) trains an FNO with \(\mathcal J_{\text{data}} + \lambda\mathcal J_{\text{pde}}\), and the interesting engineering is that the two terms live on different grids. Data supervision goes on a coarse grid, where labelled solutions are affordable; the residual goes on a fine grid, where it costs nothing but a forward pass. That split is possible only because the model is an operator, so one parameter set evaluates on both grids, and the derivatives in the residual are computed exactly in Fourier space, not by autodiff.
The numbers say the combination pays (research_notes/web-papers.md:466-501). On Kolmogorov flow, a PINN reaches 18.7 percent error in 4577 seconds while PINO reaches 0.9 percent in 536 seconds, which is 20 times better in a fifth of the time. On Burgers, PINO with data reaches 1.22 percent, PINO with the PDE loss and no paired data at all reaches 1.50 percent, and a DeepONet reaches 6.97 percent. On Darcy, 11×11 supervision plus a 61×61 residual gives 1.56 percent at four times the training resolution, against 9.46 percent for the same data with no residual. And zero-shot across resolutions PINO holds 6.04, 6.02 and 6.01 percent at 1×, 2× and 4× where a plain FNO gives 8.28 to 8.30.
Wang, Wang and Perdikaris (2021, arXiv 2103.10974) add a residual term computed by autodiff at collocation points to a DeepONet's loss, and their headline is stronger in kind, not in degree: the model trains with no paired input-output data at all, using only initial and boundary values. Reported errors are \(3.25\times10^{-3}\) on the antiderivative operator, 0.45 percent on diffusion-reaction, 1.38 percent on Burgers and \(4.22\times10^{-3}\) on the Eikonal equation (:546-573).
The Burgers number comes with a footnote. The 1.38 percent uses a hand-tuned loss weight \(\lambda = 20\), selected from the set \(\{1,5,10,20,50,100\}\). The paper has no bound tying generalization error to the residual, and it poses "what are appropriate weights?" as an open question (:559-573). PINO's own Section 5 lists comparable limits: the FFT backbone does not extend to very high dimensions, and gradient-based fine-tuning converges more slowly than refining a numerical grid (:503-507). Both papers are the Krishnapriyan diagnosis in a different form. The residual term reshapes the loss landscape, sometimes helpfully, and its weight is a hyperparameter with no theory behind it. Module 11's lab sweeps that weight on the thermal problem for exactly this reason.
LAB Not yet run. The course uses PyTorch and neuraloperator only, and DeepXDE would be a whole dependency for one experiment. Three substitutes keep the content, and all three are small.
neuraloperator ships BurgersEqnLoss and ICLoss in neuralop/losses/equation_losses.py (research_notes/repos.md:57). Add a WeightedSumLoss of data plus residual to Lab 01 run A and sweep the weight over \(\{0, 0.1, 1, 10\}\). This is the PINO recipe on the bundled data, and it is the warm-up for Module 11's lab, where the same sweep runs on the thermal problem.research_notes/web-papers.md:10-101. The five-author journal listing there is inferred, not fetched (:17-28).:105-228.:311-381. The exact algebraic exponents for the four PDE examples were not extracted into the notes; pull them from the PDF if a page needs them (:711-715).:452-519. Physics-informed DeepONet arXiv 2103.10974: :523-581. Krishnapriyan et al. arXiv 2109.01050 v2: :585-653.21-1524.pdf p.25; notes at research_notes/pdf-fno-jmlr.md:430-445.research_notes/pdf-mgn-pdebench.md:227, :243.You can now state a lower bound on the error of any fixed-sensor model and build the input pair that realises it, explain why a PINN's difficulty is in the optimizer and not the network, and say what a residual term does and does not buy when it is added to an operator's loss.