Analysis

Conditioning and identifiability diagnostics for identified parameters (pinnrl.analysis). Two questions are kept apart: how well the data determine the parameters (sensitivities, Fisher information, the Cramér–Rao bound; exact from a closed form or from the trained surrogate by the implicit-function theorem) and how well-conditioned training is (loss-Hessian and tangent-kernel spectra).

exact_sensitivities

exact_sensitivities(
    pde, parameters, x, t=None, values=None
)

du/dtheta of a closed-form solution at (x, t), shape (N, k, p).

The compiled expression is evaluated with a per-point copy of each parameter that requires grad, so the derivative comes from autograd and works for every function the expression language accepts (abs, heaviside, sigmoid, ...). It is computed in double precision.

Parameters:

Name	Type	Description	Default
pde	Any	A CustomPDE with an expression-valued exact_solution.	required
parameters	Sequence[str]	Names of the parameters to differentiate with respect to.	required
x	Tensor	Spatial coordinates, (N, d).	required
t	Tensor | None	Time, (N, 1); None for a steady problem.	None
values	dict[str, float] | None	Parameter values to evaluate at; defaults to the PDE's true values (pde.parameters), not the current trainable guesses.	None

Raises:

Type	Description
ValueError	If the PDE has no expression-valued exact solution, or a name is not one of its parameters.

Source code in pinnrl/analysis/sensitivity.py

def exact_sensitivities(
    pde: Any,
    parameters: Sequence[str],
    x: torch.Tensor,
    t: torch.Tensor | None = None,
    values: dict[str, float] | None = None,
) -> torch.Tensor:
    """``du/dtheta`` of a closed-form solution at ``(x, t)``, shape ``(N, k, p)``.

    The compiled expression is evaluated with a per-point copy of each
    parameter that requires grad, so the derivative comes from autograd and
    works for every function the expression language accepts (``abs``,
    ``heaviside``, ``sigmoid``, ...). It is computed in double precision.

    Args:
        pde: A ``CustomPDE`` with an expression-valued ``exact_solution``.
        parameters: Names of the parameters to differentiate with respect to.
        x: Spatial coordinates, ``(N, d)``.
        t: Time, ``(N, 1)``; ``None`` for a steady problem.
        values: Parameter values to evaluate at; defaults to the PDE's true
            values (``pde.parameters``), not the current trainable guesses.

    Raises:
        ValueError: If the PDE has no expression-valued exact solution, or a
            name is not one of its parameters.
    """
    compiled = getattr(pde, "_exact", None)
    if not compiled:
        raise ValueError(
            f"{pde.config.name}: exact_sensitivities needs an expression-valued exact_solution"
        )
    unknown = [p for p in parameters if p not in pde.parameters]
    if unknown:
        raise ValueError(
            f"{unknown} are not parameters of {pde.config.name}: {list(pde.parameters)}"
        )
    params = dict(pde.parameters)
    params.update(values or {})
    n = x.shape[0]
    like = torch.ones(n, 1, dtype=torch.float64)
    inputs: dict[str, Any] = {}
    for i, c in enumerate(pde.spatial_coordinates):
        inputs[c] = x[:, i : i + 1].detach().double()
    if t is not None:
        inputs["t"] = t.detach().double()
    # Per-point copies of the parameters to differentiate; the others are constants.
    thetas = {
        name: torch.full((n, 1), float(params[name]), dtype=torch.float64, requires_grad=True)
        for name in parameters
    }
    constants = {name: float(v) for name, v in params.items() if name not in thetas}
    columns = []
    for field in pde.fields:
        expr = compiled.get(field)
        if expr is None:
            columns.append(torch.zeros(n, len(parameters), dtype=torch.float64))
            continue
        u = expr({**constants, **thetas, **inputs}, like)
        per_param = []
        for name in parameters:
            if name not in expr.parameters:
                per_param.append(torch.zeros(n, dtype=torch.float64))
                continue
            (g,) = torch.autograd.grad(u.sum(), thetas[name], retain_graph=True, allow_unused=True)
            per_param.append(
                g.reshape(-1) if g is not None else torch.zeros(n, dtype=torch.float64)
            )
        columns.append(torch.stack(per_param, dim=-1))  # (N, p)
    return torch.stack(columns, dim=1).to(x.dtype)  # (N, k, p)

fisher_information

fisher_information(
    jacobian, sigma, parameters, rank_tol=None
)

Fisher information J^T Sigma^{-1} J for Gaussian noise and its Cramér–Rao bound.

Parameters:

Name	Type	Description	Default
jacobian	Tensor	Sensitivities, (N, p) or (N, k, p) (flattened over k).	required
sigma	float | Tensor	Noise standard deviation: a scalar, one value per observation point (N,), one per field (k,) or one per entry (N, k).	required
parameters	Sequence[str]	Names for the p axis.	required
rank_tol	float | None	Eigenvalues of F at or below rank_tol * lambda_max count as zero (default: p times the double-precision epsilon).	None

Raises:

Type	Description
ValueError	If sigma has none of the accepted shapes.

Source code in pinnrl/analysis/sensitivity.py

def fisher_information(
    jacobian: torch.Tensor,
    sigma: float | torch.Tensor,
    parameters: Sequence[str],
    rank_tol: float | None = None,
) -> FisherSummary:
    """Fisher information ``J^T Sigma^{-1} J`` for Gaussian noise and its Cramér–Rao bound.

    Args:
        jacobian: Sensitivities, ``(N, p)`` or ``(N, k, p)`` (flattened over ``k``).
        sigma: Noise standard deviation: a scalar, one value per observation
            point ``(N,)``, one per field ``(k,)`` or one per entry ``(N, k)``.
        parameters: Names for the ``p`` axis.
        rank_tol: Eigenvalues of ``F`` at or below ``rank_tol * lambda_max``
            count as zero (default: ``p`` times the double-precision epsilon).

    Raises:
        ValueError: If ``sigma`` has none of the accepted shapes.
    """
    if jacobian.dim() == 2:
        jacobian = jacobian[:, None, :]
    n, k, p = jacobian.shape
    J = jacobian.reshape(n * k, p).double()
    if isinstance(sigma, torch.Tensor):
        s = sigma.double()
        if s.numel() == 1:
            s = s.reshape(1).expand(n * k)
        elif tuple(s.shape) == (n, k):
            s = s.reshape(-1)
        elif tuple(s.shape) == (n,):
            s = s.repeat_interleave(k)
        elif tuple(s.shape) == (k,):
            s = s.repeat(n)
        else:
            raise ValueError(
                f"sigma must be a scalar or have shape ({n},), ({k},) or ({n}, {k}); "
                f"got {tuple(s.shape)}"
            )
        Jw = J / s[:, None]
    else:
        Jw = J / float(sigma)
    F = Jw.T @ Jw
    eigenvalues, eigenvectors = torch.linalg.eigh(F)
    order = torch.argsort(eigenvalues, descending=True)
    eigenvalues, eigenvectors = eigenvalues[order], eigenvectors[:, order]
    lam_max = float(eigenvalues[0]) if p else 0.0
    rel_tol = rank_tol if rank_tol is not None else p * torch.finfo(torch.float64).eps
    tol = rel_tol * lam_max
    keep = eigenvalues > tol
    inv = torch.where(keep, 1.0 / eigenvalues.clamp_min(1e-300), torch.zeros_like(eigenvalues))
    cov = eigenvectors @ torch.diag(inv) @ eigenvectors.T
    null = eigenvectors[:, ~keep]
    if null.shape[1]:
        identifiable = (null**2).sum(dim=1) <= 1e-12
    else:
        identifiable = torch.ones(p, dtype=torch.bool)
    std = cov.diagonal().clamp_min(0).sqrt()
    std = torch.where(identifiable, std, torch.full_like(std, math.inf))
    cov = cov.clone()
    cov.diagonal()[~identifiable] = math.inf
    denom = std[:, None] * std[None, :]
    corr = cov / denom
    corr[~torch.isfinite(denom)] = math.nan
    corr.diagonal()[identifiable] = 1.0
    smallest = float(eigenvalues[-1]) if p else 0.0
    cond = lam_max / smallest if smallest > tol and smallest > 0 else math.inf
    return FisherSummary(
        parameters=tuple(parameters),
        information=F,
        covariance=cov,
        std=std,
        correlation=corr,
        identifiable=identifiable,
        eigenvalues=eigenvalues,
        eigenvectors=eigenvectors,
        condition_number=cond,
    )

closed_form_bound

closed_form_bound(pde, parameters, x, t, sigma)

The Cramér–Rao bound for parameters from observations at (x, t) with noise sigma.

A convenience around :func:exact_sensitivities and :func:fisher_information.

Source code in pinnrl/analysis/sensitivity.py

def closed_form_bound(
    pde: Any,
    parameters: Sequence[str],
    x: torch.Tensor,
    t: torch.Tensor | None,
    sigma: float | torch.Tensor,
) -> FisherSummary:
    """The Cramér–Rao bound for ``parameters`` from observations at ``(x, t)`` with noise ``sigma``.

    A convenience around :func:`exact_sensitivities` and :func:`fisher_information`.
    """
    J = exact_sensitivities(pde, parameters, x, t)
    return fisher_information(J, sigma, parameters)

FisherSummary dataclass

FisherSummary(
    parameters,
    information,
    covariance,
    std,
    correlation,
    identifiable,
    eigenvalues,
    eigenvectors,
    condition_number,
)

Fisher information for a parameter set and what it implies.

Attributes:

Name	Type	Description
parameters	tuple[str, ...]	Parameter names, in the order of every matrix axis.
information	Tensor	The Fisher matrix F, shape (p, p).
covariance	Tensor	The Cramér–Rao bound F^{-1} on the identifiable subspace; inf on the diagonal of a parameter with a component in the null space of F.
std	Tensor	Lower bound on the standard deviation of each parameter (inf when the parameter is not identifiable from these observations).
correlation	Tensor	Correlation matrix of the bound (nan where undefined).
identifiable	Tensor	Per parameter, whether it lies outside the null space of F.
eigenvalues	Tensor	Eigenvalues of F, descending.
eigenvectors	Tensor	Columns are the corresponding parameter directions.
condition_number	float	Ratio of the largest to the smallest eigenvalue (inf when F is singular).

relative_std

relative_std(true_values)

std / |true value| per parameter (inf when not identifiable).

Source code in pinnrl/analysis/sensitivity.py

def relative_std(self, true_values: dict[str, float]) -> dict[str, float]:
    """``std / |true value|`` per parameter (``inf`` when not identifiable)."""
    return {
        name: float(self.std[i]) / abs(float(true_values[name]))
        for i, name in enumerate(self.parameters)
    }

surrogate_sensitivities

surrogate_sensitivities(
    pde,
    model,
    parameters,
    x,
    t,
    *,
    n_collocation=2048,
    damping=0.001,
    cg_iterations=300,
    cg_tol=1e-05,
    seed=0,
)

du/dtheta at (x, t) from the trained surrogate, by the implicit-function theorem.

Pass the same pde the model was trained with: the Hessian is taken of that object's physics loss (its loss weights, boundary and initial point counts, loss function and device), at the current values of its parameters. Parameters that are not registered as trainable are made trainable for the duration of the call, so a forward-trained problem works as it is.

Parameters:

Name	Type	Description	Default
pde	Any	The problem the model was trained on.	required
model	Module	The trained network, a minimum of the physics loss for the parameters' current values.	required
parameters	Sequence[str]	Names of the PDE parameters to differentiate with respect to.	required
x	Tensor	Spatial coordinates of the observation points, (N, d).	required
t	Tensor | None	Time, (N, 1); None for a steady problem (the network then receives t = 0, as the trainer does).	required
n_collocation	int	Collocation points on which the Hessian is evaluated.	2048
damping	float	Tikhonov damping lambda added to H_ww (the Hessian of a trained network is singular along flat directions). When the damped Hessian still shows negative curvature, which happens when the network is not at an exact minimum, the solve restarts with ten times the damping, up to six times.	0.001
cg_iterations	int	Maximum conjugate-gradient iterations per parameter.	300
cg_tol	float	Relative residual at which the solve stops.	1e-05
seed	int	Seed of the collocation sample (drawn in a forked RNG, so the caller's random state is untouched).	0

Raises:

Type	Description
ValueError	If a name is not a parameter of pde, or a solve diverges.

Source code in pinnrl/analysis/surrogate.py

def surrogate_sensitivities(
    pde: Any,
    model: torch.nn.Module,
    parameters: Sequence[str],
    x: torch.Tensor,
    t: torch.Tensor | None,
    *,
    n_collocation: int = 2048,
    damping: float = 1e-3,
    cg_iterations: int = 300,
    cg_tol: float = 1e-5,
    seed: int = 0,
) -> SurrogateSensitivity:
    """``du/dtheta`` at ``(x, t)`` from the trained surrogate, by the implicit-function theorem.

    Pass the same ``pde`` the model was trained with: the Hessian is taken of
    that object's physics loss (its loss weights, boundary and initial point
    counts, loss function and device), at the current values of its
    parameters. Parameters that are not registered as trainable are made
    trainable for the duration of the call, so a forward-trained problem works
    as it is.

    Args:
        pde: The problem the model was trained on.
        model: The trained network, a minimum of the physics loss for the
            parameters' current values.
        parameters: Names of the PDE parameters to differentiate with respect to.
        x: Spatial coordinates of the observation points, ``(N, d)``.
        t: Time, ``(N, 1)``; ``None`` for a steady problem (the network then
            receives ``t = 0``, as the trainer does).
        n_collocation: Collocation points on which the Hessian is evaluated.
        damping: Tikhonov damping ``lambda`` added to ``H_ww`` (the Hessian of
            a trained network is singular along flat directions). When the
            damped Hessian still shows negative curvature, which happens when
            the network is not at an exact minimum, the solve restarts with
            ten times the damping, up to six times.
        cg_iterations: Maximum conjugate-gradient iterations per parameter.
        cg_tol: Relative residual at which the solve stops.
        seed: Seed of the collocation sample (drawn in a forked RNG, so the
            caller's random state is untouched).

    Raises:
        ValueError: If a name is not a parameter of ``pde``, or a solve diverges.
    """
    unknown = [p for p in parameters if p not in getattr(pde, "parameters", {})]
    if unknown:
        raise ValueError(f"{unknown} are not parameters of {pde.config.name}")
    device = next(model.parameters()).device
    was_training = model.training
    model.eval()
    try:
        with _temporarily_trainable(pde, parameters) as thetas:
            result = _sensitivities(
                pde,
                model,
                parameters,
                thetas,
                x,
                t,
                device,
                n_collocation,
                damping,
                cg_iterations,
                cg_tol,
                seed,
            )
    finally:
        if was_training:
            model.train()
    return result

SurrogateSensitivity dataclass

SurrogateSensitivity(
    jacobian,
    parameters,
    cg_iterations,
    cg_residuals,
    damping,
    physics_loss,
)

The Jacobian du/dtheta of a trained surrogate and how the solve went.

Attributes:

Name	Type	Description
jacobian	Tensor	(N, k, p) sensitivities at the observation points.
parameters	tuple[str, ...]	The parameter names along the last axis.
cg_iterations	list[int]	Conjugate-gradient iterations used per parameter.
cg_residuals	list[float]	Relative residual of each solve.
damping	list[float]	Damping each solve ended with (raised from the requested value when the damped Hessian showed negative curvature).
physics_loss	float	The physics loss at which the Hessian was taken.

loss_hessian_spectrum

loss_hessian_spectrum(
    loss_fn, parameters, iterations=50, seed=0
)

Extreme eigenvalues of d^2 loss / d parameters^2 by Lanczos with full reorthogonalization.

Parameters:

Name	Type	Description	Default
loss_fn	Callable[[], Tensor]	Returns the scalar loss as a function of parameters (it is called once; the graph is kept for Hessian-vector products).	required
parameters	list[Parameter]	The tensors to differentiate with respect to.	required
iterations	int	Lanczos steps, hence the number of Ritz values. With 50 steps the few largest and smallest eigenvalues are converged to several digits for typical PINN Hessians.	50
seed	int	Seed of the starting vector (a forked RNG, so the caller's state is untouched).	0

Source code in pinnrl/analysis/spectrum.py

def loss_hessian_spectrum(
    loss_fn: Callable[[], torch.Tensor],
    parameters: list[torch.nn.Parameter],
    iterations: int = 50,
    seed: int = 0,
) -> HessianSpectrum:
    """Extreme eigenvalues of ``d^2 loss / d parameters^2`` by Lanczos with full reorthogonalization.

    Args:
        loss_fn: Returns the scalar loss as a function of ``parameters`` (it is
            called once; the graph is kept for Hessian-vector products).
        parameters: The tensors to differentiate with respect to.
        iterations: Lanczos steps, hence the number of Ritz values. With 50
            steps the few largest and smallest eigenvalues are converged to
            several digits for typical PINN Hessians.
        seed: Seed of the starting vector (a forked RNG, so the caller's state
            is untouched).
    """
    loss = loss_fn()
    grads = torch.autograd.grad(loss, parameters, create_graph=True, allow_unused=True)
    flat_g = torch.cat(
        [
            (g if g is not None else torch.zeros_like(p)).reshape(-1)
            for g, p in zip(grads, parameters, strict=True)
        ]
    )

    def hvp(v: torch.Tensor) -> torch.Tensor:
        hv = torch.autograd.grad(flat_g @ v, parameters, retain_graph=True, allow_unused=True)
        return torch.cat(
            [
                (h if h is not None else torch.zeros_like(p)).reshape(-1)
                for h, p in zip(hv, parameters, strict=True)
            ]
        )

    n = flat_g.numel()
    m = min(iterations, n)
    with torch.random.fork_rng(devices=[]):
        torch.manual_seed(seed)
        q = torch.randn(n, device=flat_g.device, dtype=flat_g.dtype)
    q = q / q.norm()
    Q = [q]
    alphas, betas = [], []
    for k in range(m):
        z = hvp(Q[k]).detach()
        alpha = float(Q[k] @ z)
        z = z - alpha * Q[k] - (betas[-1] * Q[k - 1] if k else 0.0)
        # Full reorthogonalization keeps the Ritz values honest at this size.
        for qq in Q:
            z = z - (qq @ z) * qq
        beta = float(z.norm())
        alphas.append(alpha)
        if beta < 1e-10 or k == m - 1:
            break
        betas.append(beta)
        Q.append(z / beta)
    T = torch.diag(torch.tensor(alphas, dtype=torch.float64))
    if betas:
        off = torch.tensor(betas[: len(alphas) - 1], dtype=torch.float64)
        T += torch.diag(off, 1) + torch.diag(off, -1)
    ritz = torch.linalg.eigvalsh(T).flip(0)
    return HessianSpectrum(eigenvalues=ritz, iterations=len(alphas), loss=float(loss.detach()))

HessianSpectrum dataclass

HessianSpectrum(eigenvalues, iterations, loss)

Ritz values of the loss Hessian from Lanczos iteration.

Attributes:

Name	Type	Description
eigenvalues	Tensor	Ritz values, descending; the extremes converge first.
iterations	int	Lanczos steps taken.
loss	float	The loss the Hessian belongs to.

largest property

largest

The largest eigenvalue (the curvature that limits the learning rate).

smallest property

smallest

The smallest Ritz value (negative when the point is not a minimum).

condition_number

condition_number(floor=1e-12)

Largest over smallest positive Ritz value.

Source code in pinnrl/analysis/spectrum.py

def condition_number(self, floor: float = 1e-12) -> float:
    """Largest over smallest positive Ritz value."""
    positive = self.eigenvalues[self.eigenvalues > floor]
    return float(positive[0] / positive[-1]) if positive.numel() else float("inf")

tangent_kernel_spectrum

tangent_kernel_spectrum(model, values_fn)

Spectrum of the tangent kernel of a per-point quantity of the network.

Parameters:

Name	Type	Description	Default
model	Module	The network.	required
values_fn	Callable[[Module], Tensor]	Maps the network to the per-point values of the term, shape (N,) or (N, k): the residual on collocation points, the boundary mismatch, the network output itself, ... Derivatives inside it must be taken with create_graph=True so that the values stay differentiable with respect to the weights.	required

The kernel is K = J J^T with J = d values / d weights of shape (N k, P), assembled row by row with one backward pass per value, so keep N k in the hundreds.

Source code in pinnrl/analysis/spectrum.py

def tangent_kernel_spectrum(
    model: torch.nn.Module,
    values_fn: Callable[[torch.nn.Module], torch.Tensor],
) -> KernelSpectrum:
    """Spectrum of the tangent kernel of a per-point quantity of the network.

    Args:
        model: The network.
        values_fn: Maps the network to the per-point values of the term, shape
            ``(N,)`` or ``(N, k)``: the residual on collocation points, the
            boundary mismatch, the network output itself, ... Derivatives
            inside it must be taken with ``create_graph=True`` so that the
            values stay differentiable with respect to the weights.

    The kernel is ``K = J J^T`` with ``J = d values / d weights`` of shape
    ``(N k, P)``, assembled row by row with one backward pass per value, so
    keep ``N k`` in the hundreds.
    """
    params = [p for p in model.parameters() if p.requires_grad]
    values = values_fn(model).reshape(-1)
    rows = []
    for i in range(values.numel()):
        grads = torch.autograd.grad(values[i], params, retain_graph=True, allow_unused=True)
        rows.append(
            torch.cat(
                [
                    (g if g is not None else torch.zeros_like(p)).reshape(-1)
                    for g, p in zip(grads, params, strict=True)
                ]
            )
        )
    J = torch.stack(rows).double()
    K = J @ J.T
    eig = torch.linalg.eigvalsh(K).flip(0).clamp_min(0)
    return KernelSpectrum(eigenvalues=eig, trace=float(eig.sum()))

KernelSpectrum dataclass

KernelSpectrum(eigenvalues, trace)

Eigenvalues of the neural tangent kernel of a loss term on a batch.

Attributes:

Name	Type	Description
eigenvalues	Tensor	Descending eigenvalues of J_w J_w^T (one per point).
trace	float	Their sum, the total rate at which the term is fitted.

condition_number

condition_number(floor=1e-12)

Largest over smallest eigenvalue above floor.

Source code in pinnrl/analysis/spectrum.py

def condition_number(self, floor: float = 1e-12) -> float:
    """Largest over smallest eigenvalue above ``floor``."""
    positive = self.eigenvalues[self.eigenvalues > floor]
    return float(positive[0] / positive[-1]) if positive.numel() else float("inf")

perturb_observations

perturb_observations(
    obs, kind, magnitude, *, seed=0, sensors=None
)

A disturbed copy of {"x", "t", "u"} observations.

Parameters:

Name	Type	Description	Default
obs	Observations	Observation tensors of shape (N, d), (N, 1), (N, k).	required
kind	str	One of "noise" (add Gaussian noise of standard deviation magnitude), "bias" (add the constant magnitude to the values), "dropout" (remove a fraction magnitude of the points at random), "delay" (the values are attributed to times t + magnitude: a sensor whose timestamps lag the physics), "outliers" (replace a fraction magnitude of the values by draws of ten times the data's standard deviation).	required
magnitude	float	Size of the disturbance, in the units the kind names.	required
seed	int	Seed of the random draws.	0
sensors	Tensor | None	Optional boolean mask (N,) selecting the points a "bias" or "delay" applies to (a single faulty probe); all points by default.	None

Raises:

Type	Description
ValueError	For an unknown kind.

Source code in pinnrl/analysis/disturbance.py

def perturb_observations(
    obs: Observations,
    kind: str,
    magnitude: float,
    *,
    seed: int = 0,
    sensors: torch.Tensor | None = None,
) -> Observations:
    """A disturbed copy of ``{"x", "t", "u"}`` observations.

    Args:
        obs: Observation tensors of shape ``(N, d)``, ``(N, 1)``, ``(N, k)``.
        kind: One of
            ``"noise"`` (add Gaussian noise of standard deviation ``magnitude``),
            ``"bias"`` (add the constant ``magnitude`` to the values),
            ``"dropout"`` (remove a fraction ``magnitude`` of the points at random),
            ``"delay"`` (the values are attributed to times ``t + magnitude``:
            a sensor whose timestamps lag the physics),
            ``"outliers"`` (replace a fraction ``magnitude`` of the values by
            draws of ten times the data's standard deviation).
        magnitude: Size of the disturbance, in the units the kind names.
        seed: Seed of the random draws.
        sensors: Optional boolean mask ``(N,)`` selecting the points a
            ``"bias"`` or ``"delay"`` applies to (a single faulty probe);
            all points by default.

    Raises:
        ValueError: For an unknown kind.
    """
    gen = torch.Generator(device="cpu").manual_seed(seed)
    x, t, u = obs["x"].clone(), obs["t"].clone(), obs["u"].clone()
    n = u.shape[0]
    mask = torch.ones(n, dtype=torch.bool) if sensors is None else sensors.bool().cpu()
    if kind == "noise":
        u = u + magnitude * torch.randn(u.shape, generator=gen).to(u)
    elif kind == "bias":
        u[mask] = u[mask] + magnitude
    elif kind == "dropout":
        keep = torch.rand(n, generator=gen) >= magnitude
        x, t, u = x[keep], t[keep], u[keep]
    elif kind == "delay":
        t[mask] = t[mask] + magnitude
    elif kind == "outliers":
        hit = torch.rand(n, generator=gen) < magnitude
        scale = 10.0 * float(u.std())
        u[hit] = scale * torch.randn((int(hit.sum()), u.shape[1]), generator=gen).to(u)
    else:
        raise ValueError(
            f"unknown disturbance {kind!r}; use noise, bias, dropout, delay or outliers"
        )
    return {"x": x, "t": t, "u": u}

linear_response

linear_response(jacobian, sigma, delta_y)

First-order parameter shift for a data disturbance delta_y.

delta_theta = (J^T Sigma^-1 J)^-1 J^T Sigma^-1 delta_y with Sigma = sigma^2 I, the response of a weighted least-squares estimate. A (N, k, p) Jacobian is flattened over k; delta_y is (N,) or (N, k). Returns (p,).

Source code in pinnrl/analysis/disturbance.py

def linear_response(
    jacobian: torch.Tensor, sigma: float | torch.Tensor, delta_y: torch.Tensor
) -> torch.Tensor:
    """First-order parameter shift for a data disturbance ``delta_y``.

    ``delta_theta = (J^T Sigma^-1 J)^-1 J^T Sigma^-1 delta_y`` with
    ``Sigma = sigma^2 I``, the response of a weighted least-squares estimate.
    A ``(N, k, p)`` Jacobian is flattened over ``k``; ``delta_y`` is ``(N,)``
    or ``(N, k)``. Returns ``(p,)``.
    """
    J = jacobian.reshape(-1, jacobian.shape[-1]).double()
    d = delta_y.reshape(-1).double()
    if isinstance(sigma, torch.Tensor):
        s = sigma.reshape(-1).double()
        if s.numel() == 1:
            s = s.expand(J.shape[0])
        Jw, dw = J / s[:, None], d / s
    else:
        Jw, dw = J / float(sigma), d / float(sigma)
    F = Jw.T @ Jw
    return torch.linalg.solve(F, Jw.T @ dw)

stability_ratio

stability_ratio(pde, model, n_points=8192, seed=0)

error / residual of a trained surrogate on a dense independent sample of the domain.

The residual is the equation applied to the network at fresh uniform points (not the training points), so it includes the quadrature gap; the error is against the closed form. Requires pde.has_exact_solution.

Source code in pinnrl/analysis/disturbance.py

def stability_ratio(
    pde: Any, model: torch.nn.Module, n_points: int = 8192, seed: int = 0
) -> StabilityRatio:
    """``error / residual`` of a trained surrogate on a dense independent sample of the domain.

    The residual is the equation applied to the network at fresh uniform
    points (not the training points), so it includes the quadrature gap;
    the error is against the closed form. Requires ``pde.has_exact_solution``.
    """
    if not pde.has_exact_solution:
        raise ValueError(f"{pde.config.name}: stability_ratio needs an exact solution")
    was_training = model.training
    model.eval()
    gen = torch.Generator(device="cpu").manual_seed(seed)
    d = max(int(pde.dimension), 1)
    lo = torch.tensor([a for a, _ in pde.domain[:d]], dtype=torch.float32)
    hi = torch.tensor([b for _, b in pde.domain[:d]], dtype=torch.float32)
    t0, t1 = pde.time_domain
    x = (lo + (hi - lo) * torch.rand(n_points, d, generator=gen)).to(pde.device)
    t = (t0 + (t1 - t0) * torch.rand(n_points, 1, generator=gen)).to(pde.device)
    residual = pde.compute_residual(model, x, t).detach()
    with torch.no_grad():
        u = model(torch.cat([x, t], dim=1))
        exact = pde.exact_solution(x, t)
        u = u[:, : exact.shape[1]]
        err = float((u - exact).norm() / exact.norm().clamp_min(1e-12))
        # Boundary: both faces of the first spatial axis at the sampled times.
        xb = torch.cat([x.clone(), x.clone()])
        xb[:n_points, 0], xb[n_points:, 0] = lo[0], hi[0]
        tb = torch.cat([t, t])
        ub = model(torch.cat([xb, tb], dim=1))[:, : exact.shape[1]]
        b_err = float((ub - pde.exact_solution(xb, tb)).pow(2).mean().sqrt())
        i_err = 0.0
        if pde.time_dependent and pde.initial_condition is not None:
            ti = torch.full((n_points, 1), float(t0), device=pde.device)
            ui = model(torch.cat([x, ti], dim=1))[:, : exact.shape[1]]
            i_err = float((ui - pde.initial_condition(x, ti)).pow(2).mean().sqrt())
    res = float(residual.pow(2).mean().sqrt())
    if was_training:
        model.train()
    return StabilityRatio(
        error=err,
        residual=res,
        boundary_error=b_err,
        initial_error=i_err,
        ratio=err / res if res > 0 else float("inf"),
        points=n_points,
    )

StabilityRatio dataclass

StabilityRatio(
    error,
    residual,
    boundary_error,
    initial_error,
    ratio,
    points,
)

Error, residual and their ratio on a dense independent grid.

Attributes:

Name	Type	Description
error	float	Relative L2 error of the field against the closed form.
residual	float	Root-mean-square residual on the grid (all equations).
boundary_error	float	RMS mismatch of the field against the closed form on the spatial boundary (zero when the exact solution is unavailable there).
initial_error	float	RMS mismatch against the initial condition at t = 0.
ratio	float	error / residual, the empirical stability constant.
points	int	Grid points used.

field_response

field_response(
    model_base, model_perturbed, x, t, reference_delta=None
)

Change of the trained field between a base and a perturbed problem, on (x, t).

Returns the relative L2 norm of u_perturbed - u_base and, when the reference solution's change reference_delta is given, the amplification ratio ||delta u_PINN|| / ||delta u_ref|| and the relative error of the PINN's response.

Source code in pinnrl/analysis/disturbance.py

def field_response(
    model_base: torch.nn.Module,
    model_perturbed: torch.nn.Module,
    x: torch.Tensor,
    t: torch.Tensor,
    reference_delta: torch.Tensor | None = None,
) -> dict[str, float]:
    """Change of the trained field between a base and a perturbed problem, on ``(x, t)``.

    Returns the relative L2 norm of ``u_perturbed - u_base`` and, when the
    reference solution's change ``reference_delta`` is given, the
    amplification ratio ``||delta u_PINN|| / ||delta u_ref||`` and the
    relative error of the PINN's response.
    """
    inputs = torch.cat([x, t], dim=1)
    with torch.no_grad():
        ub = model_base(inputs)[:, :1]
        up = model_perturbed(inputs)[:, :1]
    delta = up - ub
    out = {"response": float(delta.norm() / ub.norm().clamp_min(1e-12))}
    if reference_delta is not None:
        ref = reference_delta.reshape(delta.shape)
        out["amplification"] = float(delta.norm() / ref.norm().clamp_min(1e-12))
        out["response_error"] = float((delta - ref).norm() / ref.norm().clamp_min(1e-12))
    return out