Noisy data and identifiability

Section 8.3 recovered a diffusivity from exact readings. Real thermometers add noise, and a single training run then gives one estimate among the many that equally good data could have produced. This section asks how much the estimate moves with the noise, whether a PINN does as well as the data allow, and what happens when the data cannot determine the unknown at all.

The functions are Section 8.3's, repeated so that this section runs on its own, with the solution written for any diffusivity:

import math
import numpy as np
import torch
import torch.nn as nn
from scipy.optimize import curve_fit

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

def solution(x, t, a):
    return (torch.exp(-a * math.pi**2 * t) * torch.sin(math.pi * x)
            + 0.5 * torch.exp(-4 * a * math.pi**2 * t) * torch.sin(2 * math.pi * x))

X, T = torch.meshgrid(torch.tensor([0.1, 0.3, 0.5, 0.7, 0.9]), torch.linspace(0, 1, 21), indexing="ij")
x_d, t_d = X.reshape(-1, 1), T.reshape(-1, 1)
u_d = solution(x_d, t_d, 0.1)                       # exact readings; the true diffusivity is 0.1

def train_inverse(u_d, x_d=x_d, t_d=t_d, a_start=1.0, seed=0, epochs=5000):
    g = torch.Generator().manual_seed(seed)
    torch.manual_seed(seed)
    net = nn.Sequential(nn.Linear(2, 32), nn.Tanh(), nn.Linear(32, 32), nn.Tanh(),
                        nn.Linear(32, 32), nn.Tanh(), nn.Linear(32, 1))
    u = lambda x, t: net(torch.cat([2 * x - 1, 2 * t - 1], dim=1))
    log_a = torch.tensor([math.log(a_start)], requires_grad=True)
    x = torch.rand(1000, 1, generator=g).requires_grad_(True)
    t = torch.rand(1000, 1, generator=g).requires_grad_(True)
    t_b = torch.rand(200, 1, generator=g)
    x_b = (torch.rand(200, 1, generator=g) < 0.5).float()
    opt = torch.optim.Adam([{"params": net.parameters()}, {"params": [log_a], "lr": 1e-2}], lr=1e-3)
    for epoch in range(epochs):
        opt.zero_grad()
        u_c = u(x, t)
        r = d(u_c, t) - log_a.exp() * d(d(u_c, x), x)
        loss = (r ** 2).mean() + ((u(x_d, t_d) - u_d) ** 2).mean() + (u(x_b, t_b) ** 2).mean()
        loss.backward()
        opt.step()
    return log_a.exp().item()

Noise

Each reading gets independent noise with standard deviation $\sigma$, a fraction of the initial peak temperature of about 1. For a yardstick, the same noisy readings are also fitted by a method that knows far more than the PINN: the exact form of the solution, two decaying sine waves, with only their two amplitudes and $a$ unknown. That fit, an ordinary nonlinear least-squares problem, shows how well these readings can determine $a$ at best.

def noisy(sigma, draw):
    g = torch.Generator().manual_seed(1000 + draw)
    return u_d + sigma * torch.randn(u_d.shape, generator=g)

def fit_known_form(u_obs):
    model = lambda xt, A, B, a: (A * np.exp(-a * np.pi**2 * xt[1]) * np.sin(np.pi * xt[0])
                                 + B * np.exp(-4 * a * np.pi**2 * xt[1]) * np.sin(2 * np.pi * xt[0]))
    xt = (x_d.double().numpy().ravel(), t_d.double().numpy().ravel())
    params, _ = curve_fit(model, xt, u_obs.double().numpy().ravel(), p0=[1.0, 1.0, 1.0])
    return params[2]

print("noise     PINN estimates (three noise draws)     fit of the known form (same draws)")
for sigma in (0.01, 0.05, 0.1):
    pinn = [train_inverse(noisy(sigma, k)) for k in range(3)]
    fits = [fit_known_form(noisy(sigma, k)) for k in range(3)]
    print(f"{sigma:4.2f}      " + "  ".join(f"{a:.4f}" for a in pinn) + "           "
          + "  ".join(f"{a:.4f}" for a in fits))

noise     PINN estimates (three noise draws)     fit of the known form (same draws)
0.01      0.0987  0.1000  0.0997           0.0994  0.1006  0.1003
0.05      0.0963  0.1026  0.1007           0.0972  0.1031  0.1013
0.10      0.0935  0.1059  0.1021           0.0942  0.1063  0.1027

Two things stand out. The spread grows in proportion to the noise: about $\pm0.5\%$ of $a$ for noise of 0.01, $\pm3\%$ for 0.05 and $\pm6\%$ for 0.1. And the PINN's estimates follow the fit's, draw by draw, a little below them. The fit knows the exact form of the solution and the PINN does not, yet the PINN moves with each noise draw just as the fit does: the spread comes from the data, not from the method, and no method could do much better with these readings. The small, steady offset of about 0.6% is the PINN's own error after 5000 epochs, about what Section 8.3's run had reached at the same point.

How much a reading says about $a$

Whether data can determine an unknown comes down to sensitivity: how much the readings change when the unknown changes, compared with the noise. Here is the largest change in any of the five readings at a given time, when $a$ changes by 1%:

x_s = torch.tensor([[0.1], [0.3], [0.5], [0.7], [0.9]])
print("time t:                     " + "  ".join(f"{t:6.2f}" for t in (0.0, 0.1, 0.25, 0.5, 1.0)))
changes = [(solution(x_s, torch.full_like(x_s, t), 0.101) - solution(x_s, torch.full_like(x_s, t), 0.1)).abs().max().item()
           for t in (0.0, 0.1, 0.25, 0.5, 1.0)]
print("change for a 1% change in a: " + "  ".join(f"{c:6.4f}" for c in changes))

time t:                       0.00    0.10    0.25    0.50    1.00
change for a 1% change in a: 0.0000  0.0020  0.0033  0.0037  0.0037

At $t = 0$ the readings do not depend on $a$ at all. Later, a 1% change in $a$ moves them by at most 0.004, and by less early on: information about $a$ builds up as the temperature has time to diffuse. A change of 0.004 is smaller than the noise in a single reading even at $\sigma = 0.01$. The estimate is still good to better than 1% because a hundred readings carry the change and their noise averages out (Exercise 2). At $\sigma = 0.1$ the noise in each reading is about 25 times the effect of a 1% change, and the spread of several per cent above is what remains after the averaging.

When the data cannot decide

Now give the PINN only the five readings at $t = 0$, and train twice, from two different starting values of $a$:

first = t_d.flatten() == 0
for a_start in (1.0, 0.01):
    a_hat = train_inverse(u_d[first], x_d[first], t_d[first], a_start=a_start)
    print(f"readings at t = 0 only, starting from a = {a_start}:  estimate {a_hat:.4g}")

readings at t = 0 only, starting from a = 1.0:  estimate 0.008841
readings at t = 0 only, starting from a = 0.01:  estimate 0.00364

Both runs end far from the true 0.1, and far from each other: 0.0088 and 0.0036. Readings at a single time say nothing about how fast the temperature changes (Exercise 1), so the data term is equally content with any $a$. The residual is not: a temperature that stays as it is, with $u_t = 0$, satisfies the equation exactly when $a = 0$. With nothing in the data to stop it, the training drifts towards that trivial answer at a speed set by the optimiser, and the estimate is wherever it happens to be when training stops. A single run would have reported one of these numbers as the answer.

Beyond a single number

Repeating the training with fresh noise, as above, gives a spread but no proper statement of uncertainty, and it costs one training run per draw. Bayesian PINNs (Yang et al., 2021) treat the weights and the unknowns as random and return a whole distribution for $a$, from a single dataset, at a higher computational price. The question they answer, how sure can we be, is the right one to ask of every inverse problem; this section's checks are the cheapest way to start answering it.

Exercises

1. The sensitivity table shows the readings at $t = 0$ do not depend on $a$ at all. Explain why from the equation.
2. With noise of 0.1 in every reading, the PINN's estimates still came out within about 6% of the truth, although a 1% change in $a$ moves no reading by as much as 0.004. How can that be?
3. Section 4.2 showed that running the heat equation backwards amplifies fine detail without limit. Suppose that instead of $a$ you wanted the initial temperature, from readings at $t = 1$ only. What would you expect a PINN to recover, and what not?

Recap

• With noisy readings, an estimate is one draw from a spread. Measure the spread by repeating with fresh noise, and compare it with the best a simpler method can do with the same data.
• Identifiability is sensitivity against noise: readings that do not change when the unknown changes cannot determine it, and a PINN trained on them returns an answer that depends on where it started.
• Bayesian PINNs turn the spread into a distribution, at a price.