Scientific Machine Learning:
Embedding Physical Knowledge into Neural Networks
In many scientific and engineering problems, purely data-driven learning is insufficient. Physical systems are governed by underlying structure: conservation laws, symmetries, boundary conditions, causality. A network that does not know about this structure must rediscover it from data, and it will only ever rediscover it approximately. Incorporating the structure directly improves data efficiency, generalization and physical consistency. This notebook is a systematic catalogue of how that is done — first the general choice between penalizing a violation and making it impossible, then the standard constraints one at a time, ending with the two symmetry requirements that are easiest to confuse with each other.
Everything below runs in plain PyTorch on ordinary supervised data. Nothing in this notebook solves a differential equation or differentiates a network output — each experiment fits a network to samples and then measures whether the constraint holds. The point is always the same: compare a network with the constraint built in against the same network without it, and read off what the constraint bought.
Learning Objectives¶
After studying this notebook, you should be able to
- explain why purely data-driven learning is often insufficient for physical systems, and what is gained by building known structure into the network;
- distinguish a soft constraint, imposed as a penalty in the loss, from a hard constraint, imposed by the architecture, and state precisely what each does and does not guarantee;
- decide between the two from the certainty of the prior knowledge, the cost of a violation, and the tractability of an architectural encoding;
- enforce positivity, negativity and monotonicity exactly by transforming the raw network output, and write the corresponding penalty when no exact encoding is available;
- write boundary and initial conditions as residual penalties, and say why the resulting weights are hard to balance and why a boundary error is not a local error;
- impose the same conditions exactly instead, by the masking parameterization $\hat{v} = g + s\,v$, and construct a masking function for a given domain;
- enforce a global integral constraint by an additive shift, and positivity together with normalization by the continuous softmax;
- define equivariance $f(Tq) = T f(q)$ with a transformation operator acting on both input and output, and verify that the symmetrized network $\tilde{f}(q) = \tfrac12\bigl[f(q) + R\,f(Rq)\bigr]$ satisfies it exactly;
- define invariance $f(Tq) = f(q)$, and say when each of the two is the right requirement;
- compare the three standard routes to a symmetry — structured layers, averaging over the group, and data augmentation — in terms of guarantee and cost;
- implement each of these constructions in a few lines of PyTorch, and measure what the constraint bought: how far an unconstrained network drifts outside the admissible set, and how much accuracy the constrained one gave up in exchange.
Constraints on the output of a network arise from many sources: physical laws, mathematical requirements, conservation principles, symmetry properties. Whatever the origin, one question comes first — how strongly must the constraint hold? The answer divides every technique in this notebook into two families.
It helps to think of the question in two parts: what the prior knowledge is, and at what level it is enforced. When the knowledge is exact and certain — a boundary value that is prescribed, a total mass that is fixed — hard enforcement is appropriate. When it is approximate, uncertain, or too intricate to encode in the architecture, soft enforcement is the flexible alternative. The choice is therefore not a matter of implementation convenience. It records how much we trust the constraint and how much a violation would cost.
Setup
NumPy, Matplotlib and PyTorch.
%matplotlib inline
import numpy as np
import torch
import torch.nn as nn
from matplotlib import pyplot as plt
1.1.1. The Penalty Term¶
A soft constraint enters through the loss. Instead of forcing the network to satisfy a requirement, we add a term that grows whenever the requirement is violated, and let training push it down along with everything else.
If the output is to satisfy some requirement $\mathcal{C}$, the loss becomes
$$ \mathcal{L}(\theta) = \mathcal{L}_{\mathrm{data}}(\theta) + \lambda \, \mathcal{R}(\theta) $$
where $\mathcal{R}(\theta) \geq 0$ measures the degree of violation and $\lambda > 0$ sets how strongly it is enforced. The penalty is almost always written so that it is zero on the feasible set and smooth elsewhere, because a term that is merely an indicator of infeasibility carries no gradient and teaches the optimizer nothing.
1.1.2. What a Penalty Does and Does Not Guarantee¶
The appeal of this route is its generality: any differentiable measure of violation can be used, and nothing about the architecture has to change.
The limitation is equally plain. Minimizing $\mathcal{L}_{\mathrm{data}} + \lambda\mathcal{R}$ trades one term against the other, so at the minimum $\mathcal{R}$ is small but not zero. The constraint is satisfied approximately, and the residual violation depends on $\lambda$, on the optimizer, and on how much the data term resists. Raising $\lambda$ shrinks the violation but stiffens the optimization problem, and in the limit $\lambda \to \infty$ the loss surface becomes ill-conditioned and training slows down or stalls. There is no value of $\lambda$ at which the constraint becomes exact.
A hard constraint is enforced through the architecture, so that the output satisfies it exactly for every value of the parameters $\theta$ — at initialization, during training, and at inference. Rather than penalizing a violating output, the network is redesigned so that producing one is structurally impossible. There are two standard mechanisms.
1.2.1. Output Parameterization¶
The raw output $v(x)$ is passed through a fixed transformation that lands in the constraint-satisfying set. If the output must be non-negative,
$$ \hat{v}(x) = \exp\!\bigl(v(x)\bigr) $$
is non-negative whatever $v$ does. The transformation is part of the model, so gradients flow through it and the network learns in the transformed coordinates.
1.2.2. Projection¶
A correction is applied after the raw output is computed, mapping it onto the constraint-satisfying set. If the output must have zero mean, $\int_\Omega \hat{v}(x)\,dx = 0$, subtract the mean of the raw output:
$$ \hat{v}(x) = v(x) - \frac{1}{|\Omega|}\int_\Omega v(y)\,dy $$
This guarantees $\int_\Omega \hat{v}(x)\,dx = 0$ exactly, whatever $v$ produces.
The advantage of both mechanisms is exact satisfaction, everywhere and always. Two costs come with it. Not every constraint has a simple architectural encoding — that is the main limitation. And the transformation changes the shape of the loss surface, so a parameterization that is mathematically correct can still be numerically awkward; the exponential of Section 1.2.1 is an example, since it can produce very large values before training has settled.
1.3.1. Three Deciding Factors¶
Certainty of the prior knowledge. If the constraint is exact and known with certainty — a Dirichlet boundary condition, a conservation law — hard enforcement is appropriate. If it is approximate, domain-dependent, or derived from uncertain physical assumptions, soft enforcement is more robust, because a hard encoding of a slightly wrong constraint is a wrong model that cannot be corrected by data.
Cost of violation. In safety-critical or physically interpretable settings, even a small violation may be unacceptable: a negative predicted density is not a small error but a meaningless one. Hard enforcement removes the risk entirely. In exploratory or data-rich settings, approximate satisfaction is often sufficient.
Architectural tractability. Not every constraint admits a clean hard encoding. Nonlinear invariants, inequality constraints coupling several outputs, or conditions that depend on the input in a complicated way may be difficult or impossible to impose exactly. Then soft enforcement is the only practical option.
1.3.2. Combining the Two¶
The two are not exclusive, and in practice they are usually combined: hard constraints wherever an exact encoding is available, soft constraints for the remainder. This keeps the flexibility of penalty-based training while guaranteeing the structural requirements that matter most.
There is a second benefit worth noting. Every constraint moved from the loss into the architecture is one penalty weight $\lambda$ that no longer has to be tuned, and one competing term removed from the objective. Hard constraints therefore tend to make the remaining optimization problem easier, not merely more correct.
| soft constraint | hard constraint | |
|---|---|---|
| imposed by | a penalty term in the loss | the architecture |
| satisfaction | approximate | exact, for every $\theta$ |
| tuning | requires a weight $\lambda$ | none |
| applicability | any differentiable measure of violation | only constraints with a known encoding |
| failure mode | residual violation; ill-conditioning as $\lambda$ grows | no encoding exists, or the encoding distorts training |
A positivity constraint requires the output to be non-negative everywhere,
$$ f(x) \geq 0 \qquad \forall x . $$
This arises naturally throughout engineering. Stress intensity factors, energy densities, material stiffnesses, absolute temperatures and probability densities are all non-negative by definition, and a prediction that violates the bound is not merely inaccurate — it is outside the set of physically meaningful values.
2.1.1. Output Transformation¶
Pass the raw network output $v(x)$ through a strictly positive function:
$$ \hat{f}(x) = \exp\bigl(v(x)\bigr) $$
or the softplus variant
$$ \hat{f}(x) = \log\bigl(1 + \exp(v(x))\bigr) . $$
Both are positive for every real $v$, so the constraint holds by construction.
The two differ in how they behave away from the origin. The exponential grows without bound and can produce enormous values early in training, when $v$ is unconstrained; softplus is linear for large positive arguments and decays gently for large negative ones, so it is the more stable choice in practice. Where the output may legitimately be exactly zero, softplus is also the better model, since it approaches zero smoothly rather than requiring $v \to -\infty$.
2.1.2. Soft Enforcement¶
If a transformation is inconvenient — for instance when the same output must also satisfy other conditions that the transformation would break — add a penalty that acts only on the negative part:
$$ \mathcal{L}_{\mathrm{pos}}(\theta) = \frac{1}{N} \sum_{j=1}^{N} \bigl[ \max\bigl(0,\, -f(x_j)\bigr) \bigr]^2 $$
The term is zero when the output is non-negative at every sampled point and grows quadratically with the size of a violation. Two limitations follow from the formula itself. It is evaluated at a finite sample $\{x_j\}$, so it says nothing about points that were never sampled; and, being a penalty, it is only approximately satisfied even there.
2.2.1. Reversing the Sign¶
A negativity constraint requires the output to be non-positive,
$$ f(x) \leq 0 \qquad \forall x , $$
which arises when the quantity is bounded above by zero: a compressive stress component, a dissipation rate written as a negative quantity, a signed distance function restricted to the interior of a domain.
Everything above carries over by negating the output. The transformations become
$$ \hat{f}(x) = -\exp\bigl(v(x)\bigr) \qquad\text{or}\qquad \hat{f}(x) = -\log\bigl(1 + \exp(v(x))\bigr) , $$
and the penalty acts on the positive part instead:
$$ \mathcal{L}_{\mathrm{neg}}(\theta) = \frac{1}{N} \sum_{j=1}^{N} \bigl[ \max\bigl(0,\, f(x_j)\bigr) \bigr]^2 . $$
2.2.2. Which Strategy to Use¶
For a sign constraint the output transformation is almost always the right choice. It is exact, it costs one elementwise function evaluation, and it needs no tuning. The penalty is worth using only when the transformation cannot be applied — typically because the output already carries another hard constraint that the transformation would destroy.
Remark (a bound that is not at zero). A constraint $f(x) \geq a$ is the same construction shifted, $\hat{f} = a + \exp(v)$, and a two-sided bound $a \leq f(x) \leq b$ is obtained from the sigmoid, $\hat{f} = a + (b - a)\,\sigma(v)$. The pattern is general: to impose a set constraint exactly, find a smooth surjection from $\mathbb{R}$ onto that set and compose it with the network.
2.3.1. An Experiment¶
The target is $e^{-3x}$, which is positive everywhere. The network is trained on $[0, 1]$ and then evaluated on $[0, 2]$, so that half of the test range lies outside the training data. This is where an unconstrained network is free to do whatever its extrapolation happens to do — including going negative.
x_train = np.linspace(0, 1.0, 26) # training range
y_train = np.exp(-3*x_train) # the target, positive everywhere
x_test = np.linspace(0, 2.0, 201) # test range, twice as wide
xt = torch.tensor(x_train, dtype = torch.float32).reshape(-1, 1)
yt = torch.tensor(y_train, dtype = torch.float32).reshape(-1, 1)
xa = torch.tensor(x_test, dtype = torch.float32).reshape(-1, 1)
The two models differ in one line: the transformation applied to the raw output. Passing it through the identity is the unconstrained network; passing it through the softplus is the hard constraint of Section 2.1.1.
# the unconstrained network: the raw output is the prediction
torch.manual_seed(0)
net_plain = nn.Sequential(nn.Linear(1, 20), nn.Tanh(),
nn.Linear(20, 20), nn.Tanh(),
nn.Linear(20, 1))
opt = torch.optim.Adam(net_plain.parameters(), lr = 0.01)
for k in range(3000):
loss = torch.mean((net_plain(xt) - yt)**2)
opt.zero_grad()
loss.backward()
opt.step()
loss_plain = loss.item()
y_plain = net_plain(xa).detach().numpy().ravel()
$$ \hat{f}(x) = -\log\bigl(1 + \exp(v(x))\bigr) $$
# the constrained network: the raw output passes through a softplus, which is always positive
torch.manual_seed(0)
net_soft = nn.Sequential(nn.Linear(1, 20), nn.Tanh(),
nn.Linear(20, 20), nn.Tanh(),
nn.Linear(20, 1))
opt = torch.optim.Adam(net_soft.parameters(), lr = 0.01)
for k in range(3000):
z = net_soft(xt)
pred = torch.log(1 + torch.exp(z)) # softplus
loss = torch.mean((pred - yt)**2)
opt.zero_grad()
loss.backward()
opt.step()
loss_soft = loss.item()
y_soft = torch.log(1 + torch.exp(net_soft(xa))).detach().numpy().ravel()
print('%-18s %10s %14s %s' % ('', 'train loss', 'min prediction', 'negative points'))
print('%-18s %10.2e %+14.4f %3d / %d'
% ('plain network', loss_plain, y_plain.min(), np.sum(y_plain < 0), len(x_test)))
print('%-18s %10.2e %+14.4f %3d / %d'
% ('softplus network', loss_soft, y_soft.min(), np.sum(y_soft < 0), len(x_test)))
plt.figure(figsize = (12, 4))
plt.subplot(1, 2, 1)
plt.plot(x_test, np.exp(-3*x_test), 'k', linewidth = 2, alpha = 0.3)
plt.plot(x_test, y_plain, 'r--', linewidth = 1.5)
plt.plot(x_test, y_soft, 'b--', linewidth = 1.5)
plt.plot(x_train, y_train, 'ko', markersize = 3)
plt.plot([0, 2], [0, 0], 'k', linewidth = 0.8)
plt.legend(['exact', 'plain network', 'softplus network', 'training data'])
plt.grid(True)
plt.xlim(0.0, 2.0)
plt.xlabel('x', fontsize = 12)
plt.ylabel('f(x)', fontsize = 12)
plt.title('the whole range', fontsize = 12)
plt.subplot(1, 2, 2)
plt.plot(x_test, np.exp(-3*x_test), 'k', linewidth = 2, alpha = 0.3)
plt.plot(x_test, y_plain, 'r--', linewidth = 1.5)
plt.plot(x_test, y_soft, 'b--', linewidth = 1.5)
plt.plot([1.0, 2], [0, 0], 'k', linewidth = 0.8)
plt.grid(True)
plt.xlim(1.0, 2.0)
plt.ylim(-0.015, 0.06)
plt.xlabel('x', fontsize = 12)
plt.title('the extrapolation range, magnified')
plt.show()
Both networks fit the training range. Beyond it the unconstrained one crosses zero and stays below it, while the constrained one cannot: the softplus has no negative values in its range, so no choice of parameters and no behaviour of the extrapolation can produce one. Note also that the training loss of the constrained model is not worse — the constraint has cost nothing here, because the target satisfied it in the first place.
A monotonicity constraint requires the output to change consistently in one direction as the input increases. Monotone increasing means
$$ x_1 \leq x_2 \implies f(x_1) \leq f(x_2) , $$
as for the deflection of a beam under a growing load; monotone decreasing means
$$ x_1 \leq x_2 \implies f(x_1) \geq f(x_2) , $$
as for the fatigue life of a material against the amplitude of the applied stress.
The two cases are not separate problems. A monotone decreasing function is $f(x) = -h(x)$ with $h$ monotone increasing, so every strategy below is described for the increasing case and the decreasing case follows by negation.
3.1.1. Non-negative Derivative¶
Monotonicity is a statement about the derivative,
$$ \frac{\partial f}{\partial x} \geq 0 , $$
so one exact route is to model the derivative rather than the function, force it to be non-negative by the transformation of Chapter 2, and integrate:
$$ f(x) = f(x_0) + \int_{x_0}^{x} \exp\bigl(g(t)\bigr)\,dt $$
where $g$ is an unconstrained network. Since $\exp(g(t)) > 0$ everywhere, the integral is strictly increasing in $x$, for any $g$ and any parameters.
This gives the strongest guarantee and leaves $g$ completely free, so no expressive power is lost. Its cost is the integral, which must be evaluated numerically at every forward pass, and the constant $f(x_0)$, which becomes an additional learnable parameter.
The two cases are drawn below. What the constraint asks for is visible in the tangent at every point: it may never tip the other way.
Monotone increasing (left) and monotone decreasing (right). The dashed line is the tangent at a point: for an increasing function its slope is never negative, for a decreasing one never positive. A beam deflects more as the load grows; the fatigue life of a material falls as the stress amplitude grows.
3.1.2. Monotone Architecture¶
Constrain every weight in the network to be non-negative — in practice by storing $w$ and using $\exp(w)$ or $\mathrm{softplus}(w)$ — and use only non-decreasing activations such as ReLU or the sigmoid. A composition of non-decreasing functions is non-decreasing, so the network is monotone in its input by construction, with no integral to evaluate.
The price is expressive power. Non-negative weights mean the network can never subtract one feature from another, and such a network cannot represent functions that are monotone overall but built from competing internal terms. It is the simplest construction and the most restrictive.
3.1.3. Soft Enforcement¶
Sample pairs $x_1^i \leq x_2^i$ and penalize the violations:
$$ \mathcal{L}_{\mathrm{mono}}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \bigl[ \max\bigl(0,\, f(x_1^i) - f(x_2^i)\bigr) \bigr]^2 $$
Zero when monotonicity holds at the sampled pairs, quadratic in the size of a violation elsewhere. As always with a penalty, the guarantee extends only to what was sampled, and only approximately.
3.2.1. Reversing the Sign¶
The same three strategies apply with the sign reversed. The integral construction becomes
$$ f(x) = f(x_0) - \int_{x_0}^{x} \exp\bigl(g(t)\bigr)\,dt , $$
the architectural route constrains all weights to be non-positive — or, equivalently and more simply, applies the monotone increasing architecture to $-f(x)$ — and the penalty reverses its argument:
$$ \mathcal{L}_{\mathrm{mono}}(\theta) = \frac{1}{N} \sum_{i=1}^{N} \bigl[ \max\bigl(0,\, f(x_2^i) - f(x_1^i)\bigr) \bigr]^2 . $$
3.2.2. Which Strategy to Use¶
The non-negative derivative construction gives the strongest guarantee and keeps the network expressive; it is the default when the integral is affordable. The monotone architecture is simpler and cheaper but limits what can be represented. Soft enforcement is the easiest to add to an existing model and the only one of the three with no exact guarantee.
Remark (monotone in one input only). In several dimensions the constraint is usually required in one variable at a time — the output must increase with the load, but may do anything with the temperature. The architectural route then restricts the weights only on the path from the monotone input, and the integral route integrates only in that variable, treating the others as parameters of $g$.
3.3.1. An Experiment¶
The data are a monotone increasing curve with noise added. The noise destroys monotonicity in the data itself — it decreases at $46$ of the $100$ intervals of the grid — so a network that fits the data closely will inherit some of that. The question is whether it has to.
n = 101
x = np.linspace(0, 1, n)
dx = x[1] - x[0]
np.random.seed(0)
y_true = 1/(1 + np.exp(-10*(x - 0.5))) # a monotone increasing curve
y_data = y_true + 0.05*np.random.randn(n) # the noise destroys monotonicity
xt = torch.tensor(x, dtype = torch.float32).reshape(-1, 1)
yt = torch.tensor(y_data, dtype = torch.float32).reshape(-1, 1)
plt.figure(figsize = (6, 4))
plt.plot(x, y_data, 'ko', markersize = 3, alpha = 0.3)
plt.xlim(0, 1)
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('f(x)', fontsize = 12)
plt.show()
The unconstrained network is fitted first.
# the unconstrained network
torch.manual_seed(0)
net_plain = nn.Sequential(nn.Linear(1, 30), nn.Tanh(),
nn.Linear(30, 30), nn.Tanh(),
nn.Linear(30, 1))
opt = torch.optim.Adam(net_plain.parameters(), lr = 0.01)
for k in range(4000):
loss = torch.mean((net_plain(xt) - yt)**2)
opt.zero_grad()
loss.backward()
opt.step()
loss_plain = loss.item()
y_plain = net_plain(xt).detach().numpy().ravel()
The monotone model is the construction of Section 3.1.1, with the integral evaluated as a cumulative sum on the same grid. The network produces the rate $g$; the exponential makes it positive; the cumulative sum integrates it. The starting value $f(x_0)$ is one extra learnable parameter.
$$ f(x) = f(x_0) - \int_{x_0}^{x} \exp\bigl(g(t)\bigr)\,dt , $$
# the monotone network: the network gives the rate, the curve is its running integral
torch.manual_seed(0)
g = nn.Sequential(nn.Linear(1, 30), nn.Tanh(),
nn.Linear(30, 30), nn.Tanh(),
nn.Linear(30, 1))
f0 = torch.zeros(1, requires_grad = True) # the starting value, also learned
opt = torch.optim.Adam(list(g.parameters()) + [f0], lr = 0.01)
for k in range(4000):
rate = torch.exp(g(xt)) # exp makes the rate positive
f = f0 + dx*torch.cumsum(rate, 0) - dx*rate[0]
loss = torch.mean((f - yt)**2)
opt.zero_grad()
loss.backward()
opt.step()
loss_mono = loss.item()
rate = torch.exp(g(xt))
y_mono = (f0 + dx*torch.cumsum(rate, 0) - dx*rate[0]).detach().numpy().ravel()
d_data = np.diff(y_data) # the step between neighbouring points
d_plain = np.diff(y_plain)
d_mono = np.diff(y_mono)
print('%-18s %10s %16s %s' % ('', 'train loss', 'decreasing steps', 'most negative step'))
print('%-18s %10s %6d / %-7d %+.2e'
% ('the noisy data', '-', np.sum(d_data < 0), n - 1, d_data.min()))
print('%-18s %10.2e %6d / %-7d %+.2e'
% ('plain network', loss_plain, np.sum(d_plain < 0), n - 1, d_plain.min()))
print('%-18s %10.2e %6d / %-7d %+.2e'
% ('monotone network', loss_mono, np.sum(d_mono < 0), n - 1, d_mono.min()))
print('')
print('error against the true monotone curve')
print(' plain %.2e' % np.mean((y_plain - y_true)**2))
print(' monotone %.2e' % np.mean((y_mono - y_true)**2))
plt.figure(figsize = (6, 4))
plt.plot(x, y_data, 'ko', markersize = 3, alpha = 0.3)
plt.plot(x, y_plain, 'r--', linewidth = 1.5)
plt.plot(x, y_mono, 'b--', linewidth = 1.5)
plt.legend(['noisy data', 'plain network', 'monotone network'])
plt.xlim([0, 1])
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('f(x)', fontsize = 12)
plt.show()
Read the two error figures together. The unconstrained network decreases at twelve of the hundred intervals, following the noise; the constrained one decreases nowhere, because it cannot. And the error against the true curve is the same to three digits for both. The constraint removed every violation and cost no accuracy at all — which is what one should expect whenever the prior knowledge is correct.
For a differential equation, the solution must satisfy prescribed conditions on the boundary of the domain and, for time-dependent problems, at the initial instant. Without them the equation alone does not determine a solution, so these conditions are not decoration: they are what makes the problem well posed.
Both routes of Chapter 1 are available. The conditions can be added to the loss as residuals, which is what a physics-informed network does by default, or they can be built into the output so that they hold identically. This chapter presents the penalty first, because it is the standard arrangement and because its difficulties are what motivate the construction that follows.
4.1.1. The Boundary and Initial Residuals¶
In a physics-informed setting the loss is already a sum of residuals, and the conditions simply join it as two more terms:
$$ \mathcal{L}(\theta) = \mathcal{L}_{\mathrm{PDE}}(\theta) + \lambda_{b}\,\mathcal{L}_{b}(\theta) + \lambda_{0}\,\mathcal{L}_{0}(\theta) , $$
where the boundary term is the mean squared violation of $v = g$ at points sampled on $\partial\Omega$,
$$ \mathcal{L}_{b}(\theta) = \frac{1}{N_{b}} \sum_{j=1}^{N_{b}} \bigl( \hat{v}(x_j) - g(x_j) \bigr)^{2} , \qquad x_j \in \partial\Omega , $$
and the initial term is the same thing at $t = 0$,
$$ \mathcal{L}_{0}(\theta) = \frac{1}{N_{0}} \sum_{j=1}^{N_{0}} \bigl( \hat{v}(x_j, 0) - v_0(x_j) \bigr)^{2} . $$
Nothing about the network changes. The conditions are sampled like any other data, and the optimizer is asked to reconcile them with the equation residual.
4.1.2. What the Penalty Costs¶
This arrangement works, and for a domain whose geometry admits no convenient parameterization it is the only practical option. But it has four costs that are worth naming, because each of them disappears under the construction of Section 4.2.
The conditions compete with the equation. At the minimum of a weighted sum, no term is zero. The network settles at a compromise in which the boundary is slightly wrong and the equation is slightly wrong, and the balance point is set by $\lambda_b$ and $\lambda_0$ rather than by the physics.
The weights are hard to choose. The three terms have different units and different typical magnitudes, and the right ratio changes during training. Too small a $\lambda_b$ and the network ignores the boundary; too large and the boundary term dominates the early gradients, the network fits the edges of the domain first, and the interior is learned slowly against a stiff objective.
A boundary error is not a local error. For an elliptic problem the solution at every interior point depends on the whole boundary, so a small violation at the edge is not confined to the edge — it propagates inward and shifts the entire field. The same is true in time: an error in the initial condition is carried along the whole trajectory.
Only sampled points are constrained. The penalty says nothing about the parts of the boundary that were not sampled, which matters most exactly where the geometry is complicated and the sampling is thin.
The conditions of a well-posed problem are, however, known exactly — they are prescribed, not inferred. By the criteria of Section 1.3.1 they are the clearest candidates in this notebook for hard enforcement, and the rest of the chapter shows how.
4.2.1. The Masking Construction¶
Suppose the output must satisfy $v(x) = g(x)$ on the boundary $\partial\Omega$. Construct the network output as
$$ \hat{v}(x) = g(x) + s(x)\,v(x) $$
where $v$ is the raw network output and $s$ is a masking function with
$$ s(x) = 0 \ \text{ on } \partial\Omega , \qquad s(x) > 0 \ \text{ inside } \Omega . $$
On the boundary the second term vanishes and $\hat{v} = g$ exactly, whatever the network produces. Inside the domain $s > 0$, so the network retains full freedom, merely rescaled.
Compare this with Section 4.1. The boundary condition is now satisfied at initialization, so training never spends capacity learning it; it remains satisfied at every later iteration, so it cannot be traded against the equation residual; it holds at every boundary point, not only the sampled ones; and the weight $\lambda_b$ has disappeared from the objective altogether. The loss is left with the equation residual alone.
4.2.2. Choosing a Masking Function¶
The masking function depends on the geometry of the domain. For the unit interval $\Omega = (0, 1)$ the simplest choice is
$$ s(x) = x(1 - x) , $$
which vanishes at both endpoints and is positive in between. With homogeneous conditions $g = 0$ the parameterization reduces to
$$ \hat{v}(x) = x(1 - x)\,v(x) . $$
For a product domain the masks multiply. On the unit square $\Omega = (0,1)^2$,
$$ s(x_1, x_2) = x_1(1 - x_1)\,x_2(1 - x_2) , $$
which vanishes on all four edges by construction.
Remark (what the mask costs). A mask that decays quickly near the boundary forces $\hat{v}$ to be small there, so the network must produce large values of $v$ to represent a solution with a steep boundary layer. The construction is exact for any admissible $s$, but the conditioning of training is not indifferent to the choice; masks that vanish linearly, as above, are generally better behaved than higher-order ones. For a domain with no simple algebraic description, an approximate distance function to $\partial\Omega$ serves the same purpose.
4.3.1. Temporal Masking¶
For a time-dependent problem the output must satisfy $v(x, 0) = v_0(x)$. The same idea applies with a mask in time:
$$ \hat{v}(x, t) = v_0(x) + \tau(t)\,v(x, t) $$
where $\tau(0) = 0$ and $\tau(t) > 0$ for $t > 0$. At $t = 0$ the second term vanishes and $\hat{v}(x, 0) = v_0(x)$ exactly.
The simplest choice is $\tau(t) = t$. Another common one is
$$ \tau(t) = 1 - e^{-\alpha t} , \qquad \alpha > 0 , $$
which also vanishes at $t = 0$ but saturates instead of growing, so the network contribution is not scaled up at late times. The parameter $\alpha$ controls how quickly the network is switched on relative to the initial condition, and is worth choosing on the time scale of the problem.
4.3.2. What the Parameterization Achieves¶
In both cases the output has been split into a constraint-carrying part — the known boundary or initial value — and a free part, the network. The network is unconstrained on the interior and for $t > 0$; the condition holds by construction where it is prescribed.
The same pattern will appear again in Chapter 5 and Chapter 6: identify what must hold exactly, build a form in which it holds identically, and let the network occupy whatever freedom is left. That is the general recipe behind every hard constraint in this notebook.
In practice the two routes of this chapter are often mixed, as Section 1.3.2 described. A Dirichlet condition on a simple part of the boundary is absorbed into the mask, while a condition on an irregular part, or a Neumann condition for which no convenient mask is at hand, is left in the loss as a residual. Each condition moved into the parameterization is one weight fewer to tune and one competing term fewer in the objective.
4.4.1. An Experiment¶
The target vanishes at both ends of $(0, 1)$, and the training data are drawn from the interior only: there is no sample at $x = 0$ or $x = 1$. This is the situation the masking parameterization is built for — the condition is known exactly even though nothing in the data says so.
np.random.seed(0)
x_train = np.linspace(0.08, 0.92, 22) # no sample at x = 0 or x = 1
y_train = np.sin(np.pi*x_train)*(1 + 0.3*x_train) + 0.02*np.random.randn(22)
x_test = np.linspace(0, 1, 201)
u_test = np.sin(np.pi*x_test)*(1 + 0.3*x_test) # the exact solution, zero at both ends
xt = torch.tensor(x_train, dtype = torch.float32).reshape(-1, 1)
yt = torch.tensor(y_train, dtype = torch.float32).reshape(-1, 1)
xa = torch.tensor(x_test, dtype = torch.float32).reshape(-1, 1)
plt.figure(figsize = (6, 4))
plt.plot(x_train, y_train, 'ko', markersize = 4, alpha = 0.3)
plt.plot([0, 1], [0, 0], 'k', linewidth = 0.8)
plt.xlim([0, 1])
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('u(x)', fontsize = 12)
plt.show()
The mask is the simplest admissible one for this domain, $s(x) = x(1-x)$, and the boundary values are homogeneous, so the parameterization is $\hat{v} = s\,v$.
# the unconstrained network
torch.manual_seed(0)
net_plain = nn.Sequential(nn.Linear(1, 20), nn.Tanh(),
nn.Linear(20, 20), nn.Tanh(),
nn.Linear(20, 1))
opt = torch.optim.Adam(net_plain.parameters(), lr = 0.01)
for k in range(3000):
loss = torch.mean((net_plain(xt) - yt)**2)
opt.zero_grad()
loss.backward()
opt.step()
loss_plain = loss.item()
y_plain = net_plain(xa).detach().numpy().ravel()
$$ \hat{v}(x) = x(1 - x)\,v(x) . $$
# the masked network: the output is multiplied by s(x) = x(1 - x), which vanishes at both ends
torch.manual_seed(0)
net_mask = nn.Sequential(nn.Linear(1, 20), nn.Tanh(),
nn.Linear(20, 20), nn.Tanh(),
nn.Linear(20, 1))
opt = torch.optim.Adam(net_mask.parameters(), lr = 0.01)
for k in range(3000):
pred = xt*(1 - xt)*net_mask(xt)
loss = torch.mean((pred - yt)**2)
opt.zero_grad()
loss.backward()
opt.step()
loss_mask = loss.item()
y_mask = (xa*(1 - xa)*net_mask(xa)).detach().numpy().ravel()
print('%-16s %10s %10s %10s %s' % ('', 'train loss', 'u(0)', 'u(1)', 'max error on [0,1]'))
print('%-16s %10.2e %+.3e %+.3e %.2e'
% ('plain network', loss_plain, y_plain[0], y_plain[-1], np.max(np.abs(y_plain - u_test))))
print('%-16s %10.2e %+.3e %+.3e %.2e'
% ('masked network', loss_mask, y_mask[0], y_mask[-1], np.max(np.abs(y_mask - u_test))))
plt.figure(figsize = (6, 4))
plt.plot(x_test, y_plain, 'r--', linewidth = 1.5)
plt.plot(x_test, y_mask, 'b--', linewidth = 1.5)
plt.plot(x_train, y_train, 'ko', markersize = 4, alpha = 0.3)
plt.legend(['plain network', 'masked network', 'training data'])
plt.xlim([0, 1])
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('u(x)', fontsize = 12)
plt.show()
The masked network returns exactly zero at both ends — not approximately, and not as a result of training, but because $s$ vanishes there and the product with anything finite is zero. The unconstrained network misses by $3 \times 10^{-2}$ at one end and $7 \times 10^{-2}$ at the other, and those two numbers are the largest errors it makes anywhere: the worst part of the solution is precisely the part the prior knowledge already determined. Removing it also halves the maximum error over the whole interval.
A conservation law requires that some quantity computed from the output remains fixed. Unlike the conditions of Chapter 4, which act pointwise, a conservation law is a single scalar condition on the output as a whole — which is what makes it awkward for a network that produces its output one point at a time, and what makes the corrections below worth knowing.
5.1.1. Soft Enforcement¶
The most common form is a global integral constraint,
$$ \int_\Omega v(x)\,dx = C , $$
where $C$ is a prescribed constant such as a total mass or a total energy. The penalty is the squared violation,
$$ \mathcal{L}(\theta) = \mathcal{L}_{\mathrm{data}}(\theta) + \lambda\left(\int_\Omega \mathcal{G}_\theta(u)(x)\,dx - C\right)^{2} , $$
with $\mathcal{G}_\theta(u)$ the predicted field. The integral itself is evaluated by quadrature, so the accuracy of the constraint is limited by the accuracy of that quadrature as well as by $\lambda$.
A control volume is the physical picture behind such a constraint. Whatever crosses the boundary inwards must be balanced by what crosses it outwards, so the quantity held inside is fixed no matter how the field inside rearranges itself.
A conserved quantity as a control volume. At steady state the net flow in equals the net flow out, so the contents are fixed. The constraint is one number for the whole region, not a condition at each point — which is why it is imposed on an integral of the output rather than on the output itself.
5.1.2. Hard Enforcement by an Additive Shift¶
Because the constraint is a single scalar condition and the correction has a single degree of freedom, it can be imposed exactly by shifting the raw output by a constant:
$$ \hat{v}(x) = v(x) + \frac{C - \int_\Omega v(y)\,dy}{|\Omega|} $$
Integrating both sides over $\Omega$ gives $\int_\Omega \hat{v} = \int_\Omega v + \bigl(C - \int_\Omega v\bigr) = C$ exactly, for any $v$. This is the projection of Section 1.2.2: the raw output is mapped onto the affine set of functions with the prescribed integral, along the direction of the constants.
The shift changes the mean of the output and nothing else, so all the structure the network has learned is preserved. Its one requirement is that the integral of $v$ be available at every forward pass, which again means quadrature — here it enters the model rather than the loss.
5.2.1. The Continuous Softmax¶
When the conserved quantity must also be non-negative — a probability density, a concentration, a mass distribution — the additive shift of Section 5.1.2 will not do, because a shift can push values below zero. Both conditions are imposed at once by the continuous softmax:
$$ \hat{v}(x) = \frac{\exp\bigl(v(x)\bigr)}{\int_\Omega \exp\bigl(v(y)\bigr)\,dy} $$
The numerator is positive by the argument of Section 2.1.1, and dividing by the integral of the numerator makes the result integrate to one. Hence $\hat{v}(x) > 0$ everywhere and $\int_\Omega \hat{v}(x)\,dx = 1$, exactly, for any $v$.
This is the discrete softmax of a classification network with the sum replaced by an integral, and it inherits the same practical caution: the exponential should be evaluated after subtracting $\max_x v(x)$, so that the largest exponent is zero and the computation does not overflow.
Remark (normalizing to a constant other than one). For $\int_\Omega \hat{v} = C$ with $C \neq 1$, multiply the softmax by $C$. The construction does not extend as cleanly to several simultaneous integral constraints, which in general require a projection onto the intersection of the corresponding affine sets — a small linear solve rather than a closed-form correction.
5.3.1. An Experiment¶
A density is learned from noisy measurements. The true field is non-negative and integrates to one, but the measurements are not: the noise carries several of them below zero in the tail, where the density is nearly exhausted. The same network is trained four times — with no constraint, with the penalty of Section 5.1.1, with the additive shift of Section 5.1.2 applied to the output at every step, and with the continuous softmax of Section 5.2.1. The last two make the integral exact by construction, so the question is what each of them does to the sign.
n = 101
x = np.linspace(0, 1, n)
dx = x[1] - x[0]
np.random.seed(0)
p = np.exp(-60*(x - 0.35)**2) + 0.7*np.exp(-120*(x - 0.72)**2)
def integral(f):
return (0.5*(f[1:] + f[:-1])*dx).sum() # trapezoid rule, for arrays and for tensors
p = p/integral(p) # the true density, with integral 1
idx = np.arange(0, n, 4) # 26 noisy measurements
y_obs = p[idx] + 0.1*np.random.randn(len(idx))
xt = torch.tensor(x, dtype = torch.float32).reshape(-1, 1)
xo = xt[idx]
yo = torch.tensor(y_obs, dtype = torch.float32).reshape(-1, 1)
plt.figure(figsize = (6, 4))
plt.plot(x, p, 'k', linewidth = 2, alpha = 0.2)
plt.plot(x[idx], y_obs, 'ko', markersize = 4, alpha = 0.5)
plt.plot([0, 1], [0, 0], 'k', linewidth = 0.8)
plt.xlim([0, 1])
plt.legend(['true density', 'measurements'])
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('v(x)', fontsize = 12)
plt.show()
# no constraint: the raw output is the density
torch.manual_seed(0)
net = nn.Sequential(nn.Linear(1, 30), nn.Tanh(),
nn.Linear(30, 30), nn.Tanh(),
nn.Linear(30, 1))
opt = torch.optim.Adam(net.parameters(), lr = 0.01)
for k in range(4000):
loss = torch.mean((net(xo) - yo)**2)
opt.zero_grad()
loss.backward()
opt.step()
v_plain = net(xt).detach().numpy().ravel()
$$ \mathcal{L}(\theta) = \mathcal{L}_{\mathrm{data}}(\theta) + \lambda\left(\int_\Omega \mathcal{G}_\theta(u)(x)\,dx - C\right)^{2} , $$
# soft: the squared violation of the integral is added to the loss
torch.manual_seed(0)
net = nn.Sequential(nn.Linear(1, 30), nn.Tanh(),
nn.Linear(30, 30), nn.Tanh(),
nn.Linear(30, 1))
opt = torch.optim.Adam(net.parameters(), lr = 0.01)
for k in range(4000):
v = net(xt)
loss = torch.mean((net(xo) - yo)**2) + 10.0*(integral(v) - 1)**2
opt.zero_grad()
loss.backward()
opt.step()
v_soft = net(xt).detach().numpy().ravel()
$$ \hat{v}(x) = v(x) + \frac{C - \int_\Omega v(y)\,dy}{|\Omega|} $$
# hard: the output is shifted at every step, so the integral is exact throughout training
torch.manual_seed(0)
net = nn.Sequential(nn.Linear(1, 30), nn.Tanh(),
nn.Linear(30, 30), nn.Tanh(),
nn.Linear(30, 1))
opt = torch.optim.Adam(net.parameters(), lr = 0.01)
for k in range(4000):
v = net(xt)
v = v + (1 - integral(v))
loss = torch.mean((v[idx] - yo)**2)
opt.zero_grad()
loss.backward()
opt.step()
v = net(xt)
v_shift = (v + (1 - integral(v))).detach().numpy().ravel()
$$ \hat{v}(x) = \frac{\exp\bigl(v(x)\bigr)}{\int_\Omega \exp\bigl(v(y)\bigr)\,dy} $$
# hard: the continuous softmax gives the integral and the sign at once
torch.manual_seed(0)
net = nn.Sequential(nn.Linear(1, 30), nn.Tanh(),
nn.Linear(30, 30), nn.Tanh(),
nn.Linear(30, 1))
opt = torch.optim.Adam(net.parameters(), lr = 0.01)
for k in range(4000):
e = torch.exp(net(xt))
v = e/integral(e)
loss = torch.mean((v[idx] - yo)**2)
opt.zero_grad()
loss.backward()
opt.step()
e = torch.exp(net(xt))
v_softmax = (e/integral(e)).detach().numpy().ravel()
print('%-22s %9s %10s %9s %s' % ('', 'integral', 'min value', 'negative', 'error vs true'))
for name, v in [('plain network', v_plain), ('soft penalty', v_soft),
('additive shift', v_shift), ('continuous softmax', v_softmax)]:
print('%-22s %9.5f %+10.5f %9d %.4f'
% (name, integral(v), v.min(), np.sum(v < 0), np.sqrt(np.mean((v - p)**2))))
# the two corrections applied to the same unconstrained output, for comparison
v_a = v_plain + (1 - integral(v_plain)) # additive shift
v_b = v_plain/integral(v_plain) # rescaling
print('')
print('distance from the unconstrained output, among the corrections with the right integral')
print(' additive shift %.5f' % np.sqrt(integral((v_a - v_plain)**2)))
print(' rescaling %.5f' % np.sqrt(integral((v_b - v_plain)**2)))
plt.figure(figsize = (11, 4))
plt.subplot(1, 2, 1)
plt.plot(x, p, 'k', linewidth = 2, alpha = 0.2)
plt.plot(x[idx], y_obs, 'ko', markersize = 4, alpha = 0.5)
plt.plot(x, v_plain, 'r--', linewidth = 1.5)
plt.plot(x, v_softmax, 'b--', linewidth = 1.5)
plt.plot([0, 1], [0, 0], 'k', linewidth = 0.8)
plt.legend(['true density', 'measurements', 'plain network', 'continuous softmax'], fontsize = 10)
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('v(x)', fontsize = 12)
plt.subplot(1, 2, 2)
plt.plot(x, p, 'k', linewidth = 2, alpha = 0.2)
plt.plot(x, v_plain, 'r--', linewidth = 1.5)
plt.plot(x, v_shift, 'g--', linewidth = 2)
plt.plot(x, v_softmax, 'b--', linewidth = 1.5)
plt.plot([0, 1], [0, 0], 'k', linewidth = 0.8)
plt.xlim(0.82, 1.0)
plt.ylim(-0.12, 0.25)
plt.legend(['true density', 'plain network', 'additive shift', 'continuous softmax'], fontsize = 10)
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.title('the tail, magnified', fontsize = 12)
plt.show()
The unconstrained network misses the prescribed integral by three percent, and the penalty brings it to three parts in a thousand — closer, but never exact, and what remains is the price of a finite $\lambda$. The shift and the softmax are exact to every digit printed, because in their case the integral is not something the optimizer is asked to achieve.
What separates the two is the sign. The shift subtracts one constant from the whole field, and the constant it must subtract here is positive, so the correction pushes the tail further below zero: the plain network is negative at two grid points, the shifted one at five, and its lowest value is twice as far below zero. The softmax has no negative values to lose, and it satisfies both requirements at once, no further from the true density than the shift is. The last two lines confirm the claim of Section 5.1.2 along the way — among the corrections that give the right integral, the shift is the smaller one — which is exactly the point: a correction can be optimal in one norm and still be inadmissible, and when the quantity is a density the constraint to be satisfied is the pair.
Many networks are asked to learn mappings that respect a symmetry of the underlying physics. If the input is rotated, translated or reflected, the output should transform accordingly. Building that structure into the network improves data efficiency and guarantees consistent predictions on transformed inputs, which no amount of training on untransformed data can.
6.1.1. Definition¶
Let $f$ be a network and $T$ a transformation — a rotation, a translation, a reflection. The network is equivariant with respect to $T$ when it commutes with it:
$$ f(Tq) = T\,f(q) $$
Transforming the input and then applying the network gives the same result as applying the network and then transforming the output.
Two points about the notation are worth making explicit, because they are the source of most confusion later. First, the same symbol $T$ appears on both sides but acts on different spaces: on the left it transforms the input of the network, on the right it transforms the output. Second, $T$ acts on whole functions, not on numbers. If the input is a field $q$ defined on a domain, then $Tq$ is another field, and the action of $T$ on the output is whatever the corresponding transformation of that output is — a rotation of a velocity field rotates both the positions and the vectors.
6.1.2. Examples¶
A network predicting the velocity field of a fluid should be equivariant to rotation: rotating the whole flow configuration by ninety degrees should rotate the predicted velocity field by ninety degrees.
A network predicting a temperature distribution from a heat source should be equivariant to translation: shifting the source to the right should shift the entire predicted field by the same amount.
In both cases the network that is equivariant responds correctly to a transformed input without ever having been trained on one. A network that is not equivariant may respond to the rotated configuration with a field that is not the rotation of anything it predicted before — and no inspection of its loss on the training set would reveal that.
The figure below collects the transformations that occur most often. In each case the same $T$ acts on the input and on the output, which is the content of $f(Tx) = T f(x)$.
Transformations that commute with an equivariant network. Translation, rotation, reflection and permutation of labelled objects are the four that appear most often. Equivariance is the requirement that moving the input by $T$ moves the output by the same $T$.
6.2.1. Group-Equivariant Layers¶
Build the symmetry into the layers, by constructing weight-sharing patterns that commute with the transformation. The canonical example is the convolution layer, which is equivariant to translation: the same filter is applied at every position, so shifting the input shifts the output by the same amount. Group convolutional networks extend the idea to rotations by sharing weights across rotated copies of each filter.
This gives exact equivariance at one forward pass, which is the strongest of the three routes. It requires designing the layer for the group in question, and not every group has a convenient one.
6.2.2. Output Symmetrization¶
Leave the architecture alone and average the network over the group. For a finite set of transformations $\{T_1, T_2, \ldots, T_n\}$,
$$ \tilde{f}(q) = \frac{1}{n} \sum_{i=1}^{n} T_i^{-1} f(T_i q) $$
Each term transforms the input by $T_i$, applies the network, and then undoes the transformation on the output with $T_i^{-1}$. Every term is therefore a prediction expressed in the original frame, and averaging them gives a result that no longer depends on which frame was used.
The construction is exact for any $f$, needs no change to the architecture, and costs $n$ forward passes per input.
6.2.3. Data Augmentation¶
Show the network transformed versions of each training example. From a training pair $(q_i, y_i)$, generate
$$ (T q_i,\; T y_i) $$
for various $T$ — note that the label is transformed as well, since that is what equivariance asserts.
This requires no change to the architecture or to inference, and is by far the easiest route. It is also the weakest: the network is encouraged to be equivariant on the transformations and inputs it has seen, and may still be measurably inequivariant on unseen ones.
Of the three, group-equivariant layers give the strongest guarantee, data augmentation the weakest. Output symmetrization sits in between: exact and architecture-free, at the price of $n$ forward passes.
6.3.1. The Reflection Operator¶
Consider a clamped-clamped beam on $x \in [-L/2,\, L/2]$, with the origin at the centre. Given a distributed load $q(x)$, the network predicts the deflection, $v = f(q)$.
The relevant transformation is the reflection about the centre of the beam. Written as an operator acting on a function,
$$ (Rq)(x) = q(-x) , $$
so $Rq$ is the load obtained by flipping the beam left to right. The same operator acts on deflections, $(Rv)(x) = v(-x)$. Reflecting twice returns the original function, so
$$ R^2 = I \qquad\Longrightarrow\qquad R^{-1} = R , $$
and $\{I, R\}$ is a group with two elements.
The physics is equivariant under $R$. Flipping the load flips the deflection and does nothing else, because the beam and its supports are themselves symmetric about $x = 0$. Writing $f^{\star}$ for the exact solution operator,
$$ f^{\star}(Rq) = R\,f^{\star}(q) . $$
A network $f$ trained without any symmetry constraint does not satisfy this. Nothing in a data-fitting loss forces it to.
The figure below is that statement drawn once. The top row sends a load through the network and obtains a deflection; the bottom row sends the mirrored load through the same network, and the deflection that comes out is the mirror of the first.
Equivariance of the clamped-clamped beam. Top: a load $q$ and the deflection $v = f(q)$ predicted from it. Bottom: the mirrored load $Rq$, whose deflection is the mirrored deflection $R\,f(q)$. The beam and its supports are symmetric about $x = 0$, so nothing distinguishes the two situations except the reflection itself.
6.3.2. The Symmetrized Network¶
Apply the construction of Section 6.2.2 to the group $\{I, R\}$. With $n = 2$, $T_1 = I$ and $T_2 = R$, and using $R^{-1} = R$,
$$ \tilde{f}(q) = \frac{1}{2}\Bigl[\, f(q) + R\,f(Rq) \,\Bigr] $$
Read it as three operations. Pass the load through the network. Pass the flipped load through the same network and flip the resulting deflection back. Average the two.
Written pointwise, with $v = f(q)$ and $\bar{v} = f(Rq)$,
$$ \tilde{v}(x) = \frac{1}{2}\bigl[\, v(x) + \bar{v}(-x) \,\bigr] . $$
The construction is exactly equivariant, and the verification is two lines. Using $R^2 = I$,
$$ \tilde{f}(Rq) = \frac{1}{2}\Bigl[ f(Rq) + R\,f(R^{2}q) \Bigr] = \frac{1}{2}\Bigl[ f(Rq) + R\,f(q) \Bigr] , $$
while applying $R$ to the output gives
$$ R\,\tilde{f}(q) = \frac{1}{2}\Bigl[ R\,f(q) + R^{2}f(Rq) \Bigr] = \frac{1}{2}\Bigl[ R\,f(q) + f(Rq) \Bigr] . $$
The two right-hand sides are the same, so
$$ \tilde{f}(Rq) = R\,\tilde{f}(q) $$
holds identically — for every load, and for every value of the parameters, even if the network $f$ has learned nothing about the symmetry at all. The cost is two forward passes instead of one.
The construction is easiest to read as a picture.
The two terms of the average. The upper beam carries the original load and gives $f(q)$; the lower one carries the flipped load, and its deflection is flipped back to give $R\,f(Rq)$. Both are now expressed in the same frame, so averaging them is meaningful — and the average is equivariant whatever the network has learned.
6.3.3. A Symmetric Load as a Special Case¶
Equivariance is a statement about all loads. The familiar special case follows from it in one line.
A symmetric load is one that the reflection leaves unchanged — a fixed point of $R$:
$$ Rq = q , \qquad\text{that is}\qquad q(-x) = q(x) . $$
Apply the equivariance established in Section 6.3.2 to such a load. Since $\tilde{f}(Rq) = R\,\tilde{f}(q)$ holds for every load, and $Rq = q$ for this one,
$$ R\,\tilde{f}(q) = \tilde{f}(Rq) = \tilde{f}(q) , $$
which written pointwise is
$$ \tilde{v}(-x) = \tilde{v}(x) . $$
The predicted deflection is even. This is what one expects physically — a beam that is symmetric about its centre, loaded symmetrically, deflects symmetrically — but it is worth seeing that nothing extra had to be imposed to obtain it. Symmetry of the response is a consequence of equivariance, applied to a load that happens to be a fixed point of $R$.
What the construction reduces to. For a symmetric load the two forward passes of Section 6.3.2 are calls to the network on the same input, since $Rq = q$. The symmetrized prediction
$$ \tilde{f}(q) = \frac{1}{2}\Bigl[\, f(q) + R\,f(Rq) \,\Bigr] = \frac{1}{2}\Bigl[\, f(q) + R\,f(q) \,\Bigr] $$
is therefore obtained in three steps and one forward pass:
- Pass the load $q$ through the network to obtain the raw deflection $v = f(q)$.
- Flip the output: $Rv$, that is $v(-x)$.
- Average the two.
$$ \tilde{v}(x) = \frac{v(x) + v(-x)}{2} $$
In words, the symmetrized prediction is the even part of the raw prediction. Whatever antisymmetric component the network produced — and an unconstrained network will produce one, since nothing in a data-fitting loss forbids it — is discarded, and $\tilde{v}(-x) = \tilde{v}(x)$ holds by construction.
6.4.1. An Experiment¶
The beam of Section 6.3 is the test case. A unit point load acts at position $c$ on the beam, and the network reads $c$ and returns the deflection at the $m = 21$ interior grid points. For a point load the exact deflection is available in closed form: with $a = c + L/2$ the distance of the load from the left support and $b = L - a$,
$$ v(s) = \frac{b^{2} s^{2}\,\bigl[\,3aL - (3a + b)\,s\,\bigr]}{6\,EI\,L^{3}} \quad (s \le a) , $$
and the same expression with $a$ and $b$ exchanged and $s$ replaced by $L - s$ for $s \ge a$. It serves as the truth. The reflection $R$ sends a load at $c$ to a load at $-c$ and reverses the order of the grid values of the deflection, so the network is asked to satisfy $f(-c) = R\,f(c)$.
The symmetry is withheld from the data on purpose: every training load lies on the left half of the beam, $c \in [-0.35, -0.05]$, while the test loads lie anywhere, $c \in [-0.40, 0.40]$. Three of the routes of Section 6.2 are compared against the unconstrained network — data augmentation, the symmetrization of Section 6.2.2 applied to a network that has already been trained, and the same symmetrization built in from the start, so that what the optimizer sees is $\tilde{f}$ rather than $f$.
L = 1.0
m = 21 # interior grid points
x = np.linspace(-L/2, L/2, m + 2)[1:-1] # the beam, origin at the centre
def deflection(c): # exact deflection for a unit point load at x = c
a = c + L/2 # distance of the load from the left support
b = L - a
s = x + L/2 # distance of each grid point from the left support
v = np.zeros(m)
left = s <= a
v[left] = b**2*s[left]**2*(3*a*L - (3*a + b)*s[left])/(6*L**3)
r = L - s[~left]
v[~left] = a**2*r**2*(3*b*L - (3*b + a)*r)/(6*L**3)
return v
R = np.eye(m)[::-1].copy() # reflection: reverses the order of the entries
print('R is an involution || R R - I || = %.2e' % np.linalg.norm(R.dot(R) - np.eye(m)))
print('the exact operator is equivariant || f(Rq) - R f(q) || = %.2e'
% np.linalg.norm(deflection(-0.2) - R.dot(deflection(0.2))))
print('centre deflection %.6e P L^3 / 192 EI = %.6e' % (deflection(0.0).max(), 1/192))
rng = np.random.default_rng(0)
c_train = rng.uniform(-0.35, -0.05, 300) # every training load sits on the left half
rng = np.random.default_rng(1)
c_test = rng.uniform(-0.40, 0.40, 300) # the test loads sit anywhere
V_train = np.zeros((300, m))
V_test = np.zeros((300, m))
for i in range(300):
V_train[i] = deflection(c_train[i])
V_test[i] = deflection(c_test[i])
scale = 1/np.abs(V_train).max() # the deflections are of order 1e-3
ct = torch.tensor(c_train.reshape(-1, 1), dtype = torch.float32)
vt = torch.tensor(V_train*scale, dtype = torch.float32)
# the unconstrained network
torch.manual_seed(0)
net = nn.Sequential(nn.Linear(1, 64), nn.Tanh(),
nn.Linear(64, 64), nn.Tanh(),
nn.Linear(64, m))
opt = torch.optim.Adam(net.parameters(), lr = 0.002)
for k in range(3000):
loss = torch.mean((net(ct) - vt)**2)
opt.zero_grad()
loss.backward()
opt.step()
print('plain network train loss = %.2e' % loss.item())
# data augmentation: each load and its mirror, with the label mirrored as well
c_aug = np.concatenate([c_train, -c_train])
V_aug = np.vstack([V_train, V_train[:, ::-1]])
ca = torch.tensor(c_aug.reshape(-1, 1), dtype = torch.float32)
va = torch.tensor(V_aug*scale, dtype = torch.float32)
torch.manual_seed(0)
net_aug = nn.Sequential(nn.Linear(1, 64), nn.Tanh(),
nn.Linear(64, 64), nn.Tanh(),
nn.Linear(64, m))
opt = torch.optim.Adam(net_aug.parameters(), lr = 0.002)
for k in range(3000):
loss = torch.mean((net_aug(ca) - va)**2)
opt.zero_grad()
loss.backward()
opt.step()
print('data augmentation train loss = %.2e' % loss.item())
# the symmetrized model, trained as such: the optimizer never sees the raw network alone
torch.manual_seed(0)
net_sym = nn.Sequential(nn.Linear(1, 64), nn.Tanh(),
nn.Linear(64, 64), nn.Tanh(),
nn.Linear(64, m))
opt = torch.optim.Adam(net_sym.parameters(), lr = 0.002)
for k in range(3000):
pred = 0.5*(net_sym(ct) + torch.flip(net_sym(-ct), [1]))
loss = torch.mean((pred - vt)**2)
opt.zero_grad()
loss.backward()
opt.step()
print('symmetrized network train loss = %.2e' % loss.item())
Two numbers are reported for each model. The equivariance error is the relative discrepancy $\lVert f(-c) - R\,f(c)\rVert / \lVert f(c)\rVert$ averaged over the test loads; it is a property of the map alone and never looks at the exact solution. The test RMSE compares the prediction with the exact deflection, and is the accuracy.
def predict(network, c):
with torch.no_grad():
return network(torch.tensor(c.reshape(-1, 1), dtype = torch.float32)).numpy()/scale
def symmetrize(network, c): # 0.5*[ f(q) + R f(Rq) ]
return 0.5*(predict(network, c) + predict(network, -c)[:, ::-1])
print('%-26s %18s %14s' % ('', 'equivariance error', 'test RMSE'))
for name, P, PR in [('plain network', predict(net, c_test), predict(net, -c_test)),
('data augmentation', predict(net_aug, c_test), predict(net_aug, -c_test)),
('symmetrized, trained', symmetrize(net_sym, c_test), symmetrize(net_sym, -c_test))]:
error = np.mean(np.linalg.norm(PR - P[:, ::-1], axis = 1)/np.linalg.norm(P, axis = 1))
print('%-26s %18.2e %14.3e' % (name, error, np.sqrt(np.mean((P - V_test)**2))))
c_one = np.array([0.25]) # a load on the right half, never seen in training
v_exact = deflection(c_one[0])
v_raw = predict(net, c_one)[0]
v_mirror = predict(net, -c_one)[0][::-1]
v_sym = symmetrize(net_sym, c_one)[0]
plt.figure(figsize = (12, 4))
plt.subplot(1, 2, 1)
plt.plot([c_one[0], c_one[0]], [0, 1], 'k', linewidth = 2)
plt.plot([-c_one[0], -c_one[0]], [0, 1], 'k--', linewidth = 1.5)
plt.xlim(-0.5, 0.5)
plt.legend(['load at c', 'mirrored load'], fontsize = 11)
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('q(x)', fontsize = 12)
plt.subplot(1, 2, 2)
plt.plot(x, v_exact, 'k', linewidth = 3, alpha = 0.3)
plt.plot(x, v_sym, 'b--', linewidth = 1.5)
plt.legend(['exact', 'symmetrized network'], fontsize = 11)
plt.grid(True)
plt.xlabel('x', fontsize = 12)
plt.ylabel('v(x)', fontsize = 12)
plt.show()
The unconstrained network is not equivariant by any margin — the discrepancy between the two sides of the identity is as large as the prediction itself — and on loads placed where it was never trained its error is two thirds of the deflection it is meant to produce. It learned the left half of the problem and nothing about the right.
Augmentation repairs most of that. The equivariance error falls to about one percent and the accuracy becomes good, because the mirrored copies are in effect the missing half of the training set. What it does not deliver is the identity itself: one percent is small, but it is not zero, and nothing keeps it small on a load unlike those it was shown.
Both symmetrized rows are exactly equivariant, to every digit printed, and the difference between them is the point of the experiment. Symmetrizing a network that has already been trained buys the identity immediately and little else — the error falls only from $1.5\times10^{-3}$ to $1.1\times10^{-3}$, because the average of one good prediction and one useless one is a half-useless prediction, and no algebraic identity can supply information the network never saw. Training the symmetrized model instead gives both at once: the identity for free, and an error of $1.4\times10^{-5}$ from data covering half the beam, because the construction supplies the other half before the optimizer ever sees it.
The figure is the same statement for a single load. The two curves $f(q)$ and $R\,f(Rq)$ are what equivariance requires to be equal; for the unconstrained network the first is wrong everywhere while the second, which reaches the network through the half of the beam it does know, lies on the exact solution. Averaging them is what the symmetrized model does, and training with that average in place is what makes the average right.
Equivariance is the right requirement when the output carries direction or position. When it does not — when the output is a single number that should not change at all under the transformation — the requirement is invariance, and the constructions change in one small but decisive way.
7.1.1. Definition¶
A network is invariant with respect to a transformation $T$ when
$$ f(Tq) = f(q) . $$
This is the special case of equivariance in which $T$ acts trivially on the output.
A network predicting the total energy of a molecule should return the same value however the molecule is rotated or translated, because total energy is a scalar. A network classifying an image as "cat" or "not cat" should return the same label when the image is shifted a few pixels.
7.1.2. Equivariance or Invariance?¶
The distinction is decided by the output, not by the transformation.
If the output is a field or a vector — something with a direction or a position — equivariance is the requirement, and the transformation must be undone on the output. If the output is a single number with no directional content, invariance is the requirement, and there is nothing to undo.
The two symmetrization formulas differ by exactly that factor:
$$ \tilde{f}(q) = \frac{1}{2}\bigl[\, f(q) + R\,f(Rq) \,\bigr] \qquad\text{(equivariant)} , $$
$$ \tilde{f}(q) = \frac{1}{2}\bigl[\, f(q) + f(Rq) \,\bigr] \qquad\text{(invariant)} . $$
7.2.1. Invariant Input Features¶
Transform the raw input into features that are unchanged by the transformation, and pass only those to the network. For rotational invariance, pairwise distances between points serve, since distances do not change under rotation. The network then operates entirely on invariant quantities and its output is invariant automatically.
This is the strongest and most efficient of the three, because the invariance is built into the representation and costs one forward pass. It requires knowing a complete set of invariant features; an incomplete set discards information the network may need.
The idea is shown below for two interacting particles.
A complete invariant for two particles. The positions $(x_1, y_1)$ and $(x_2, y_2)$ change under rotation and translation, but the distance $r_{12} = \lVert r_1 - r_2\rVert$ between them does not. A network fed $(q_1, q_2, r_{12})$ rather than the raw coordinates returns an energy $U$ that is invariant by construction.
7.2.2. Pooling over Transformations¶
Evaluate the network on all transformed versions of the input and aggregate:
$$ \tilde{f}(q) = \frac{1}{n} \sum_{i=1}^{n} f(T_i q) $$
Because the sum runs over the whole group, transforming the input merely permutes the terms and the result is unchanged. A maximum can replace the average,
$$ \tilde{f}(q) = \max_{i}\, f(T_i q) , $$
with the same argument. This is exact and architecture-free, at the price of $n$ forward passes.
7.2.3. Data Augmentation¶
Train on transformed versions of each example. From a pair $(q_i, y_i)$ with a scalar label, generate
$$ (T q_i,\; y_i) . $$
Unlike the equivariant case of Section 6.2.3, the label is left unchanged, which is precisely the difference between the two requirements. The route is easy and approximate: the network is encouraged towards invariance but is not forced into it.
7.3.1. Why the Energy Is Invariant¶
Take the clamped-clamped beam of Section 6.3 again, on $x \in [-L/2,\, L/2]$ with load $q(x)$, but let the network predict the total strain energy
$$ U = \int_{-L/2}^{L/2} \frac{1}{2}\, EI \left(\frac{d^{2}v}{dx^{2}}\right)^{2} dx , $$
where $EI$ is the bending stiffness.
Flipping the load flips the deflection, by the equivariance established in Section 6.3.1 — but the energy is an integral of the squared curvature over the whole beam, and reflecting the integrand merely reverses the order in which it is integrated. The energy is therefore unchanged:
$$ U[Rq] = U[q] . $$
The deflection field is equivariant; the scalar computed from it is invariant. Both statements describe the same physics, and which one applies to a given network is settled by what that network outputs.
Accordingly, a network predicting the energy should satisfy
$$ f(Rq) = f(q) , $$
and if this is not enforced it will generally predict two different energies for $q$ and $Rq$ — two different answers to what is physically one question.
7.3.2. The Three Strategies Applied¶
Invariant input features. Symmetrize the load before it reaches the network:
$$ \tilde{q}(x) = q(x) + q(-x) . $$
Since $R\tilde{q} = \tilde{q}$, the network only ever sees symmetric inputs and its output is invariant. Note what this costs: the antisymmetric part of the load has been discarded, and it does contribute to the energy. Here the feature is invariant but not complete, which is the general hazard of this route.
Pooling over transformations. Evaluate the network on both the original and the flipped load and average:
$$ \tilde{f}(q) = \frac{1}{2}\bigl[\, f(q) + f(Rq) \,\bigr] . $$
No inverse transformation appears, because the output is a scalar. The result satisfies $\tilde{f}(Rq) = \tilde{f}(q)$ identically, whatever $f$ has learned, at the cost of two forward passes. This is the invariant counterpart of the construction of Section 6.3.2, and comparing the two formulas side by side is the clearest way to hold the distinction.
Data augmentation. Train on both $(q_i, U_i)$ and $(Rq_i, U_i)$, with the label unchanged. The network is encouraged to return the same energy for both, but may still differ slightly on unseen loads. It is the easiest to implement, requiring no change to the architecture or to inference, and the only one of the three with no exact guarantee.
7.4.1. An Experiment¶
Two point charges in the plane interact with the Coulomb energy
$$ U = \frac{k\,q_1 q_2}{r_{12}} , \qquad r_{12} = \lVert r_1 - r_2 \rVert , $$
which depends on the separation alone and not on where the pair sits or how it is turned. Two networks are trained on the same energies: one reads the four raw coordinates, the other reads only $r_{12}$. To make the difference visible, the training configurations are restricted to a quarter turn of orientation, while the test configurations use every orientation.
k = 1.0
q1 = 1.0
q2 = -1.0
m = 2000
# training pairs: the orientation covers a quarter turn only
rng = np.random.default_rng(0)
r_train = rng.uniform(0.5, 2.0, m) # separation
th = rng.uniform(0, np.pi/4, m) # orientation
c = rng.uniform(-1, 1, (m, 2)) # centre of the pair
d = 0.5*np.stack([r_train*np.cos(th), r_train*np.sin(th)], axis = 1)
X_train = np.hstack([c + d, c - d]) # the four raw coordinates
U_train = k*q1*q2/r_train
# test pairs: every orientation
rng = np.random.default_rng(1)
r_test = rng.uniform(0.5, 2.0, m)
th = rng.uniform(0, 2*np.pi, m)
c = rng.uniform(-1, 1, (m, 2))
d = 0.5*np.stack([r_test*np.cos(th), r_test*np.sin(th)], axis = 1)
X_test = np.hstack([c + d, c - d])
U_test = k*q1*q2/r_test
# the network that reads the four raw coordinates
torch.manual_seed(0)
net_raw = nn.Sequential(nn.Linear(4, 40), nn.Tanh(),
nn.Linear(40, 40), nn.Tanh(),
nn.Linear(40, 1))
a = torch.tensor(X_train, dtype = torch.float32)
b = torch.tensor(U_train, dtype = torch.float32).reshape(-1, 1)
opt = torch.optim.Adam(net_raw.parameters(), lr = 0.005)
for i in range(4000):
loss = torch.mean((net_raw(a) - b)**2)
opt.zero_grad()
loss.backward()
opt.step()
p_train_raw = net_raw(a).detach().numpy().ravel()
p_test_raw = net_raw(torch.tensor(X_test, dtype = torch.float32)).detach().numpy().ravel()
# the network that reads the separation only
torch.manual_seed(0)
net_inv = nn.Sequential(nn.Linear(1, 40), nn.Tanh(),
nn.Linear(40, 40), nn.Tanh(),
nn.Linear(40, 1))
a = torch.tensor(r_train.reshape(-1, 1), dtype = torch.float32)
b = torch.tensor(U_train, dtype = torch.float32).reshape(-1, 1)
opt = torch.optim.Adam(net_inv.parameters(), lr = 0.005)
for i in range(4000):
loss = torch.mean((net_inv(a) - b)**2)
opt.zero_grad()
loss.backward()
opt.step()
p_train_inv = net_inv(a).detach().numpy().ravel()
p_test_inv = net_inv(torch.tensor(r_test.reshape(-1, 1), dtype = torch.float32)).detach().numpy().ravel()
print('%-20s %10s %s' % ('', 'train RMSE', 'test RMSE (all orientations)'))
print('%-20s %10.4f %10.4f'
% ('raw coordinates', np.sqrt(np.mean((p_train_raw - U_train)**2)),
np.sqrt(np.mean((p_test_raw - U_test)**2))))
print('%-20s %10.4f %10.4f'
% ('distance r12 only', np.sqrt(np.mean((p_train_inv - U_train)**2)),
np.sqrt(np.mean((p_test_inv - U_test)**2))))
The aggregate numbers already tell the story, but the direct test is better: take one configuration, rotate it through a full turn without changing the separation, and watch the prediction.
r0 = 1.2
angles = np.linspace(0, 2*np.pi, 73)
pred_raw = np.zeros(73)
for i in range(73): # the same pair, turned through a full circle
dx1 = 0.5*r0*np.cos(angles[i])
dy1 = 0.5*r0*np.sin(angles[i])
pair = np.array([[dx1, dy1, -dx1, -dy1]])
pred_raw[i] = net_raw(torch.tensor(pair, dtype = torch.float32)).detach().item()
pred_inv = net_inv(torch.tensor([[r0]], dtype = torch.float32)).detach().item()
print('one configuration at r = %.1f, rotated through a full turn' % r0)
print(' exact energy %+.4f' % (k*q1*q2/r0))
print(' raw-coordinate network : min %+.4f max %+.4f spread %.4f'
% (pred_raw.min(), pred_raw.max(), pred_raw.max() - pred_raw.min()))
print(' distance-only network : %+.4f (constant, the input never changes)' % pred_inv)
plt.figure(figsize = (6, 4))
plt.plot(angles, pred_raw, 'r--', linewidth = 1.5)
plt.plot(angles, pred_inv*np.ones(73), 'b--', linewidth = 1.5)
plt.plot(angles, (k*q1*q2/r0)*np.ones(73), 'k', linewidth = 4, alpha = 0.3)
plt.legend(['raw coordinates', 'distance only', 'exact'], fontsize = 12)
plt.grid(True)
plt.xlabel('orientation of the pair', fontsize = 12)
plt.ylabel('predicted energy', fontsize = 12)
plt.show()
The network given the raw coordinates fits its training set to $0.005$ and then misses the test set by $1.37$ — a factor of about $260$ — and its prediction for one physical configuration varies by $2.3$ as that configuration is merely turned around. The network given $r_{12}$ has the same accuracy on both sets, and its prediction is exactly constant under rotation — its flat line lies on top of the exact one in the figure — because the number it receives does not change. The invariance was not learned; it was made unnecessary to learn.
Key takeaways
| Constraint | Hard enforcement | Soft enforcement | | :--- | :--- | :--- | | $f \geq 0$ | $\hat{f} = \exp(v)$ or $\log(1 + e^{v})$ | $\bigl[\max(0, -f)\bigr]^2$ | | $f \leq 0$ | $\hat{f} = -\exp(v)$ | $\bigl[\max(0, f)\bigr]^2$ | | monotone increasing | $f(x_0) + \int_{x_0}^{x}\exp(g)\,dt$; non-negative weights | $\bigl[\max(0, f(x_1) - f(x_2))\bigr]^2$ | | $v = g$ on $\partial\Omega$ | $\hat{v} = g + s\,v$, $\;s|_{\partial\Omega} = 0$ | $\bigl(\hat{v} - g\bigr)^2$ on sampled boundary points | | $v(x,0) = v_0$ | $\hat{v} = v_0 + \tau(t)\,v$, $\;\tau(0) = 0$ | $\bigl(\hat{v}(\cdot,0) - v_0\bigr)^2$ | | $\int_\Omega v = C$ | $\hat{v} = v + \bigl(C - \int_\Omega v\bigr)/|\Omega|$ | $\bigl(\int_\Omega v - C\bigr)^2$ | | $v > 0$ and $\int_\Omega v = 1$ | continuous softmax | both penalties | | equivariance | group layers; $\tfrac12[f(q) + R f(Rq)]$ | augmentation with transformed labels | | invariance | invariant features; $\tfrac12[f(q) + f(Rq)]$ | augmentation with unchanged labels |
- A soft constraint is a term in the loss; a hard constraint is a property of the architecture. The first applies to anything differentiable and is satisfied approximately; the second is exact for every parameter value but exists only where an encoding is known. Raising the penalty weight does not turn one into the other — it only stiffens the optimization.
- The recipe behind every hard constraint is the same. Identify what must hold exactly, write a form in which it holds identically, and let the network occupy the freedom that remains: $\exp(v)$ for a sign, $g + s\,v$ for a boundary value, an additive shift for an integral, a group average for a symmetry.
- Sign and bound constraints are essentially free. A smooth surjection from $\mathbb{R}$ onto the admissible set, composed with the network, costs one elementwise evaluation and removes a penalty weight from the objective.
- Monotonicity is a statement about the derivative, so the exact constructions either model the derivative and integrate it, or restrict the weights so that the composition cannot decrease. The second is simpler and strictly less expressive.
- Boundary and initial conditions are the clearest case for hard enforcement: they are known exactly, they hold from initialization onwards, and they can never be traded against the equation residual during training. As penalties they compete with the equation, need two weights that are hard to balance, and leave an error at the edge that an elliptic problem carries into the whole interior.
- Equivariance means $f(Tq) = T f(q)$, with the transformation acting on the input on one side and on the output on the other. Averaging over a group with the transformation undone on the output, $\tilde{f}(q) = \tfrac12\bigl[f(q) + R f(Rq)\bigr]$ for a reflection, is exactly equivariant for any network, at the price of one extra forward pass.
- Invariance means $f(Tq) = f(q)$, and its symmetrization is the same formula with the output transformation removed. Which of the two is wanted is decided by the output: a field or a vector calls for equivariance, a scalar for invariance.
- The symmetry of a particular solution is a consequence, not the requirement. A symmetric load is a fixed point of $R$, so $R\tilde{f}(q) = \tilde{f}(Rq) = \tilde{f}(q)$ and the predicted deflection is even. Building only for that case would leave every asymmetric load unconstrained.
- The constructions cost little and the measurements say so. In the experiments of this notebook the constrained model matched the unconstrained one on training loss and on error against the truth, while removing every violation: the softplus network never went negative where the plain one did over a fifth of the test range; the monotone network decreased nowhere where the plain one decreased at twelve intervals; the masked network was exactly zero at both ends, where the plain network's largest error was at those ends. A correct prior is not a restriction that has to be paid for.
- An exact construction is exact for a reason that can be checked. Commuting with the reflection to machine precision, returning zero at a masked boundary, integrating to the prescribed constant — each is one line of verification, and each fails loudly when the construction is subtly wrong.
Looking ahead. Every construction in this notebook fixes what the network cannot get wrong; none of them says anything about what it should get right in between. That is the work of the training objective, and the two meet in physics-informed learning, where the governing equation itself becomes the loss and the constructions above decide which of its conditions are paid for by a penalty weight and which are imposed for free. The same distinction reappears at the level of the architecture in operator learning, where a translation-equivariant spectral layer is chosen for exactly the reason given in Section 6.2.1 — a symmetry that is built in is a symmetry that never has to be learned from data.