PINN or classical solver?

Section 1.4 said that for a clean forward problem a classical solver is usually faster and more accurate than a PINN, and listed the places where a PINN earns its keep. Eleven chapters later, that claim can be measured. The fair way to do it is to put each PINN of this book beside the method a numerical analyst would reach for, on the same computer, and to look at both the error and the time. Three problems: the forward problem of Chapter 9, the inverse problem of Section 8.3, and a Poisson problem in ten dimensions.

A forward problem in one dimension

Chapter 9's problem is $-u'' = (4\pi)^2\sin(4\pi x)$ on $(0, 1)$ with $u(0) = u(1) = 0$. Two classical methods: second-order finite differences, which replace $u''$ at each grid point by $(u_{j-1} - 2u_j + u_{j+1})/h^2$ and solve one tridiagonal linear system (LeVeque, 2007); and a Chebyshev spectral method, which represents $u$ by the polynomial through its values at Chebyshev points and differentiates that polynomial exactly, using the differentiation matrix of Trefethen (2000).

import math
import time
import numpy as np
from scipy.linalg import solve_banded

k = 4 * math.pi                                    # -u'' = k^2 sin(kx), u(0) = u(1) = 0; u = sin(kx)

def finite_differences(n):
    """Second-order differences on n interior points: one tridiagonal solve."""
    x = np.linspace(0, 1, n + 2)
    h = x[1] - x[0]
    bands = np.zeros((3, n))
    bands[0, 1:], bands[1], bands[2, :-1] = -1 / h**2, 2 / h**2, -1 / h**2
    u = np.zeros(n + 2)
    u[1:-1] = solve_banded((1, 1), bands, k**2 * np.sin(k * x[1:-1]))
    return x, u

def chebyshev(n):
    """Chebyshev collocation with n interior points (the differentiation matrix of Trefethen's cheb)."""
    m = n + 1
    s = np.cos(np.pi * np.arange(m + 1) / m)       # Chebyshev points on [-1, 1]
    c = np.hstack([2.0, np.ones(m - 1), 2.0]) * (-1.0) ** np.arange(m + 1)
    D = np.outer(c, 1 / c) / (s[:, None] - s[None, :] + np.eye(m + 1))
    D -= np.diag(D.sum(axis=1))
    x = (s + 1) / 2                                # map to [0, 1], so that d/dx = 2 d/ds
    u = np.zeros(m + 1)
    u[1:-1] = np.linalg.solve(-4 * (D @ D)[1:-1, 1:-1], k**2 * np.sin(k * x[1:-1]))
    return x, u

print("method               unknowns    relative L2 error at the points    milliseconds per solve")
for method, sizes in [(finite_differences, (16, 128, 1024, 8192)), (chebyshev, (16, 24, 32))]:
    for n in sizes:
        start = time.perf_counter()
        for _ in range(100):
            x, u = method(n)
        ms = (time.perf_counter() - start) * 1000 / 100
        error = np.linalg.norm(u - np.sin(k * x)) / np.linalg.norm(np.sin(k * x))
        print(f"{method.__name__:19s}  {n:6d}       {error:.1e}                          {ms:7.3f}")

method               unknowns    relative L2 error at the points    milliseconds per solve
finite_differences       16       4.7e-02                            0.046
finite_differences      128       7.9e-04                            0.064
finite_differences     1024       1.3e-05                            0.120
finite_differences     8192       2.0e-07                            0.466
chebyshev                16       1.7e-08                            0.052
chebyshev                24       3.7e-15                            0.049
chebyshev                32       9.3e-15                            0.058

With 128 unknowns, as many as the collocation points of Chapter 9's PINNs, finite differences reach an error of $7.9\times10^{-4}$ in a fraction of a millisecond. Section 9.3's PINNs took two to twelve seconds for errors between $5\times10^{-4}$ and $2\times10^{-2}$, and Levenberg–Marquardt, the most accurate optimiser of Section 9.4, eighteen seconds for $4\times10^{-6}$ to $4\times10^{-5}$. Each eightfold refinement of the grid divides the finite-difference error by about 60, the factor of $8^2$ that a second-order method promises, and 8192 unknowns give $2\times10^{-7}$ in half a millisecond: more accurate than every PINN in Chapter 9, and more than ten thousand times faster. The spectral method is in another class again: 24 unknowns reach the limit of double precision.

The gap is not an accident of our code. A classical solver uses what is known about the problem: the equation is linear, so one linear solve gives the answer; the discrete operator is sparse, so the solve is cheap; and the method's theory says how fast the error falls as the grid is refined, which is what makes a convergence study a check on the answer (Section 4.4). A PINN turns the same linear problem into a nonconvex optimisation that comes with no such guarantee, and its error has to be checked afterwards.

Grossmann et al. (2024) made the comparison systematically, against the finite element method, for the Poisson equation in one, two and three dimensions, the Allen–Cahn equation in one and a semilinear Schrödinger equation in one and two. In their study the PINNs did not beat the finite element method. The finite element solutions were obtained one to three orders of magnitude faster for the Poisson equation, two for the Schrödinger equation in one dimension and five to six for Allen–Cahn, and they were more accurate in most cases. Where the PINNs did better was afterwards: for the Poisson problems in two and three dimensions, evaluating a trained network at new points was two to three orders of magnitude faster than interpolating the finite element solution onto them.

An inverse problem

Section 8.3 recovered the diffusivity $a$ of a rod, and its initial temperature, from 105 readings of five thermometers. Its PINN trained for 10000 epochs, about a minute on the computer that ran this book's code, and ended with $a = 0.0998$, 0.2% from the truth, and the initial temperature within $1.2\times10^{-2}$ everywhere.

The classical route wraps a solver in an optimiser. Guess $a$ and the initial temperature, solve the heat equation forward, compare the solver's temperatures at the thermometers with the readings, and let a least-squares method improve the guesses until they match. The solver is Crank–Nicolson on 99 interior grid points (LeVeque, 2007). The initial temperature is an unknown function, so it needs a representation: six sine waves, $u(x, 0) = \sum_{m=1}^{6} c_m\sin(m\pi x)$, which vanish at the ends as the boundary condition requires. That makes seven unknowns, $\log a$ and the six amplitudes. The optimiser starts from $a = 1$, as the PINN did, and from a rod at zero temperature.

from scipy.optimize import least_squares

a_true = 0.1                                       # used only to make the readings, as in Section 8.3
def truth(x, t):
    return (np.exp(-a_true * math.pi**2 * t) * np.sin(math.pi * x)
            + 0.5 * np.exp(-4 * a_true * math.pi**2 * t) * np.sin(2 * math.pi * x))

x_d, t_d = np.array([0.1, 0.3, 0.5, 0.7, 0.9]), np.linspace(0, 1, 21)
readings = truth(x_d[:, None], t_d[None, :])       # five thermometers, 21 readings each

x = np.linspace(0, 1, 101)[1:-1]                   # the solver's grid: 99 interior points
h, dt, steps = 0.01, 0.005, 10                     # ten time steps between readings
D2 = (np.diag(-2 * np.ones(99)) + np.diag(np.ones(98), 1) + np.diag(np.ones(98), -1)) / h**2
at_sensors = [int(np.argmin(abs(x - s))) for s in x_d]
sines = np.sin(np.pi * np.outer(x, np.arange(1, 7)))   # the initial temperature: six sines

def simulate(a, u0):
    """Crank-Nicolson for u_t = a u_xx with u = 0 at the ends; the temperatures at the sensors."""
    step = np.linalg.solve(np.eye(99) - 0.5 * dt * a * D2, np.eye(99) + 0.5 * dt * a * D2)
    between = np.linalg.matrix_power(step, steps)
    u, out = u0, [u0[at_sensors]]
    for _ in range(len(t_d) - 1):
        u = between @ u
        out.append(u[at_sensors])
    return np.array(out).T

simulations = 0
def misfit(p):                                     # p = (log a, the six amplitudes)
    global simulations
    simulations += 1
    return (simulate(math.exp(p[0]), sines @ p[1:]) - readings).ravel()

start = time.perf_counter()
fit = least_squares(misfit, np.r_[math.log(1.0), np.zeros(6)], method="lm")
seconds = time.perf_counter() - start
a_hat = math.exp(fit.x[0])
xs = np.linspace(0, 1, 201)
u0_hat = np.sin(np.pi * np.outer(xs, np.arange(1, 7))) @ fit.x[1:]
print(f"estimate of a: {a_hat:.6f}, {abs(a_hat - a_true) / a_true:.1e} from the truth")
print(f"largest error in the initial temperature: {np.abs(u0_hat - truth(xs, 0)).max():.1e}")
print(f"{simulations} simulations of the rod, {seconds:.2f} seconds in all")

estimate of a: 0.100013, 1.3e-04 from the truth
largest error in the initial temperature: 4.6e-05
64 simulations of the rod, 0.05 seconds in all

The classical fit finds $a$ with a relative error of $1.3\times10^{-4}$, set by the solver's grid and time step rather than by the fit, and the initial temperature within $5\times10^{-5}$, after 64 simulations of the rod that took a twentieth of a second in all. That is fifteen times closer than the PINN in $a$, over two hundred times closer in the initial temperature, and a thousand times faster.

What matters more is what the classical route needed that the PINN did not:

• A solver for the forward problem. Twenty lines here, for a linear equation on an interval. For a new equation on a three-dimensional geometry, writing and verifying one can be weeks of work. The PINN for a new equation is the old one with a new residual.
• A representation of the unknown function. Six sines is a choice, and it is the classical method's regularisation. With exact readings, any number from two to nine recovers $a$ to about one part in ten thousand; with noisy readings, more sines let the initial temperature follow the noise (Exercise 2). A PINN makes the equivalent choice implicitly, through its architecture and its training.
• Derivatives with respect to the unknowns. least_squares estimated them by finite differences, one extra simulation per unknown, which is why seven unknowns cost 64 simulations. For an unknown that varies across a domain and has thousands of values, that is hopeless. The classical answer is an adjoint solver, which gives the gradient of the misfit with respect to all the unknowns for the cost of about one extra solve, but has to be derived and written for each equation. A PINN gets the same gradient from automatic differentiation, whatever the unknown.

This is the case that reviews of the field make for physics-informed learning (Karniadakis et al., 2021): classical inverse problems with hidden physics can be prohibitively expensive, and need their own formulations and elaborate codes. For a small linear problem with a solver to hand, like this one, the classical code is still the better tool. The balance tips towards the PINN as the forward solver gets harder to write and the unknowns more numerous.

Many dimensions

Everything so far was in one dimension, where a grid is cheap. The cost of a grid grows exponentially with the dimension. Finite differences need about twenty points in each direction for an error under 1% on the solution below (0.75% in one dimension, Exercise 1, and a little less in two and three), and in $d$ dimensions that makes $20^d$ unknowns: 400 in two dimensions, 8000 in three, 3.2 million in five, and $10^{13}$ in ten, eighty terabytes for the solution alone, before any matrix. This is the curse of dimensionality. Equations in many dimensions are not exotic: in quantum mechanics, finance and control, the dimension counts particles, assets or agents.

A PINN builds no grid. Its points are drawn at random, and each one tests the equation in every direction at once. Here is the Poisson equation $-\nabla^2 u = f$ on the ten-dimensional cube $[0, 1]^{10}$, with $f$ and the values on the faces chosen so that

$$u = \frac{1}{\sqrt{10}}\sum_{i=1}^{10}\sin(2\pi x_i),$$

a full sine wave in every direction, with mean zero and a spread of about 0.7 over the cube. Each step draws 500 fresh points inside the cube and 500 on its faces (Section 7.1), the residual is divided by its natural scale $4\pi^2$ (Section 4.5), and the learning rate decays from $10^{-3}$ (Section 9.2):

import torch
import torch.nn as nn

def d(y, x):
    return torch.autograd.grad(y, x, torch.ones_like(y), create_graph=True)[0]

dim = 10
exact = lambda p: torch.sin(2 * math.pi * p).sum(dim=1, keepdim=True) / math.sqrt(dim)   # -lap u = 4 pi^2 u
g = torch.Generator().manual_seed(0)

def inside(n):
    return torch.rand(n, dim, generator=g).requires_grad_(True)

def on_faces(n):                                        # each point on one of the 20 faces, chosen at random
    p = torch.rand(n, dim, generator=g)
    p[torch.arange(n), torch.randint(0, dim, (n,), generator=g)] = (torch.rand(n, generator=g) < 0.5).float()
    return p

torch.manual_seed(0)
net = nn.Sequential(nn.Linear(dim, 64), nn.Tanh(), nn.Linear(64, 64), nn.Tanh(),
                    nn.Linear(64, 64), nn.Tanh(), nn.Linear(64, 1))
u = lambda p: net(2 * p - 1)

def residual(p):
    grad_u = d(u(p), p)
    laplacian = sum(d(grad_u[:, i:i + 1], p)[:, i:i + 1] for i in range(dim))
    return (-laplacian - 4 * math.pi**2 * exact(p)) / (4 * math.pi**2)

epochs, w_b = 6000, 10.0
opt = torch.optim.Adam(net.parameters(), lr=1e-3)
schedule = torch.optim.lr_scheduler.CosineAnnealingLR(opt, epochs, eta_min=1e-5)
x_test = torch.rand(10000, dim, generator=torch.Generator().manual_seed(1))
start = time.time()
for epoch in range(1, epochs + 1):
    p, q = inside(500), on_faces(500)                   # fresh points at every step
    opt.zero_grad()
    loss = (residual(p) ** 2).mean() + w_b * ((u(q) - exact(q)) ** 2).mean()
    loss.backward()
    opt.step()
    schedule.step()
    if epoch % 1000 == 0:
        with torch.no_grad():
            error = ((u(x_test) - exact(x_test)).norm() / exact(x_test).norm()).item()
        print(f"epoch {epoch}:  relative L2 error at 10000 random points {error:.1e}   {time.time() - start:4.0f} s")

epoch 1000:  relative L2 error at 10000 random points 8.0e-02     25 s
epoch 2000:  relative L2 error at 10000 random points 5.1e-02     50 s
epoch 3000:  relative L2 error at 10000 random points 4.2e-02     76 s
epoch 4000:  relative L2 error at 10000 random points 4.1e-02    100 s
epoch 5000:  relative L2 error at 10000 random points 3.6e-02    125 s
epoch 6000:  relative L2 error at 10000 random points 3.5e-02    149 s

After 6000 steps, two and a half minutes, the error at 10000 random points of the cube is 3.5%, and it is still falling slowly. That is far from the accuracy of the one-dimensional runs of Chapter 9, and typical of PINNs in many dimensions: the training sets the error, and more of it would help. The weight of 10 on the boundary term is there because the twenty faces of the cube are a large part of this problem. An error left on them spreads into the interior as a function whose Laplacian is zero, which the residual cannot see, so the weight tips Section 10.1's compromise towards the faces; with equal weights the same run ended at 4.9%. But the comparison that matters is with a grid, and there is none to make. A grid with as many points as one of these batches would have fewer than two points in each direction.

Three cautions keep this honest. The solution is a sum of one-dimensional functions, which suits a network, and a classical method that knew it could exploit it too: sparse grids and low-rank tensor methods are built for structure of this kind. Monte Carlo methods, which estimate the solution at one point by averaging over random walks that start there, also need no grid, and their cost grows only slowly with the dimension; they are the classical competitor here. And the error above was measured against an exact solution that a real problem would not have. Without one, the checks of Sections 12.1 to 12.3 are all there is.

Two early demonstrations that neural networks can solve equations in many dimensions were the deep BSDE method of Han et al. (2018), for parabolic equations including the nonlinear Black–Scholes, Hamilton–Jacobi–Bellman and Allen–Cahn equations in a hundred dimensions, and the Deep Galerkin Method of Sirignano & Spiliopoulos (2018), which, like a PINN, trains on randomly sampled points instead of a mesh, and solved free-boundary problems in up to 200 dimensions. Of all the reasons to prefer a mesh-free neural method, the dimension is the one for which grids have no answer at all.

Choosing

The measurements of this section and the studies above add up to a short guide:

None of these rows is a law; each is where the evidence points today. Whichever tool you choose, a PINN is only as good as its checks: the residual at fresh points, several seeds, and an independent reference wherever one can be had. The worked problems of Part 5 do exactly that, one problem at a time.

In the Studio

The Studio's domains are a line, a rectangle and a disc, with time as an option: one to three dimensions, where a classical solver is the benchmark to beat. To compare the two on a problem without an exact solution, set the Report block's Rows to Regular grid, export the field, and solve the same problem with a classical method on the same grid, as Section 13.2 does with finite elements. Then report the difference between the two, and the time each took, as the box above describes.

Exercises

1. The error of second-order finite differences for $\sin(kx)$ is close to $(kh)^2/12$, with $h$ the grid spacing. Check this against the table. Use it to estimate the error for $\sin(2\pi x)$ with 20 interior points, and the number of unknowns that would give $10^{-6}$ on Chapter 9's problem.
2. In the inverse problem, add noise with standard deviation 0.01 to the readings and fit again with 2, 6 and 9 sines. Which changes more with the number of sines, the estimate of $a$ or the initial temperature, and why?
3. The ten-dimensional PINN used 500 interior points at every step, yet the cost of a step grows with the dimension. Why? Roughly how would it change from ten dimensions to a hundred?

Recap

• For a forward problem in one to three dimensions, a classical solver is faster by orders of magnitude and more accurate: here finite differences beat every PINN of Chapter 9 in both, and a spectral method reached rounding error.
• A PINN's advantages lie elsewhere: one code for forward and inverse problems with any unknown, equations in many dimensions, and fast, differentiable evaluation once trained.
• Compare at equal accuracy or equal time, against the method a specialist would use, and count the training.