10. Adjoints & Differentiable FEM

PDE-constrained least squares, the discrete adjoint method, tangent vs adjoint sensitivities, unrolled vs implicit differentiation, differentiable FEM solvers, contact and the quality of gradients

Contents

1. From One-Shot Least Squares to Gradient-Based Inversion 2. The Adjoint Method: One Extra Solve for the Whole Gradient 3. Worked Example: Recovering a Stiffness Profile in a 1D Bar 4. Five Route Families to the Same Gradient 5. Differentiable FEM: Autodiff Meets the Stiffness Matrix 6. Dynamics, Contact, and the Quality of Gradients 7. Practical Recipes for a Soft Tactile Skin Interactive: Deriving the Discrete Adjoint Interactive: Adjoint Bar Inverter Interactive: Gradient-Routes Laboratory Flashcards

1. From One-Shot Least Squares to Gradient-Based Inversion

Every inversion in this guide so far, Module 7's FEMU aside, was inverted once. Tessler–Spangler shape sensing (Modules 4–6) and Tikhonov-regularized force reconstruction (Module 8) share a structural gift: the unknowns enter the observation model linearly — small-strain kinematics hold, superposition holds — so the least-squares functional is quadratic, the normal equations are assembled and solved one time, and the matrix can even be prefactorized offline. No iterations. No gradients. That gift is exactly what this module gives up.

Core Problem

Recover parameters $\mathbf{p}$ (a modulus map, a contact traction, a geometry) from measurements $\mathbf{d}$ when the predicted response $\mathbf{u}(\mathbf{p})$ is defined only implicitly, through an equilibrium residual:

\min_{\mathbf{p}} \; J(\mathbf{p}) = \tfrac{1}{2}\,\lVert \mathbf{S}\,\mathbf{u}(\mathbf{p}) - \mathbf{d} \rVert^2 + \alpha\,\Phi(\mathbf{p}) \quad \text{s.t.} \quad \mathbf{R}\big(\mathbf{u}(\mathbf{p}),\mathbf{p}\big) = \mathbf{0}

In the guide's observation model $y = \mathcal{H}(u, p, f, X_0) + b + \eta$: the response $u$ is partially measured through $\mathbf{S}$, and the unknown has moved inside the forward operator. Iterative optimization needs $\nabla_{\mathbf{p}} J$; this module is a menu of ways to get it — cheaply, exactly, or automatically.

	In this module
Measured	Displacements or strains at $m$ points — for a vision-based tactile skin, the camera-tracked marker displacement field $\mathbf{d}$
Unknown	$n_p$ parameters $\mathbf{p}$ inside the forward model: distributed contact traction (inside the operator only once the mechanics is finite-strain or the contact patch itself is unknown — a prescribed traction on a linear model enters the right-hand side and stays the one-shot problem of Module 8), a modulus map $E(\mathbf{x})$ of the elastomer, or the skin geometry
Assumed known	Mesh, boundary conditions, constitutive form, load (when material is unknown) or material (when load is unknown) — whatever is not in $\mathbf{p}$
Unobservable	Parameters with zero sensitivity — regions the forward strain or the sensor set never illuminates (Section 3 makes this exact)

Three Ways Linearity Dies

Three inverse-problem families break the one-shot structure:

Material identification. $E(\mathbf{x})$ enters the stiffness matrix, so $\mathbf{u}(\mathbf{p}) = \mathbf{K}(\mathbf{p})^{-1}\mathbf{f}$ depends nonlinearly on $\mathbf{p}$ even for linear elasticity. Module 7's FEMU already met this; here we make its gradient affordable.
Large-deformation hyperelasticity. A silicone tactile skin routinely sees 10–50% strain. Superposition fails; equilibrium $\mathbf{R}(\mathbf{u},\mathbf{p})=\mathbf{0}$ is nonlinear in $\mathbf{u}$ itself, and every objective evaluation is a Newton solve.
Unknown geometry or contact. The shape of an inclusion, the location of a contact patch — parameters that move nodes or switch constraint sets.

All three become nonlinear least squares with a PDE constraint. The optimizer (Gauss–Newton, L-BFGS — Module 3's machinery) is standard; the entire game is the gradient. The naive route, finite differences, perturbs one parameter at a time:

\text{cost of } \nabla_{\mathbf{p}} J: \quad \underbrace{(n_p+1)\,C_{\text{solve}}}_{\text{finite differences}} \quad \text{vs.} \quad \underbrace{C_{\text{solve}} + C_{\text{lin}}}_{\text{adjoint}}

For a modulus map with $n_p = 10^4$ elements, finite differences means ten thousand FEM solves per gradient evaluation, with accuracy capped near $\sqrt{\varepsilon_{\text{mach}}}$ for one-sided differences. The adjoint method delivers the same gradient — exactly, for the discrete objective — at the price of one forward solve plus one linear solve, independent of $n_p$. Even when the forward problem is itself linear the two halves of that bill are not equal: the adjoint system carries the same matrix, so it reuses the factorization the forward solve already paid for and costs a pair of triangular substitutions, not a second $C_{\text{solve}}$ (Section 2 makes this precise). For the nonlinear problems above the gap only widens — the forward solve is a Newton loop, the adjoint solve a single linear system. The routes on the menu:

(a) Finite differences — $(n_p+1)$ full FEM solves per gradient; step-size dilemma; Section 4.
(b) Direct (tangent) sensitivities — one linearized solve per parameter; wins when $n_p$ is small and outputs many; Section 2.
(c) The adjoint method — two solves per gradient regardless of $n_p$ (plus $O(n_p)$ work to assemble the components); Sections 2–3.
(d) Algorithmic differentiation of the entire FEM code — differentiable simulation, by unrolling the solver or by implicit differentiation at the fixed point; automates the tedious $\partial \mathbf{R}/\partial \mathbf{p}$ terms; Sections 4–5.

Two further families — physics-informed neural networks, which fold the PDE into a training loss and make unknown coefficients trainable, and neural operators, which amortize inversion over an offline dataset so each new inversion costs milliseconds — are deliberately not covered here. They are learning methods with their own failure modes (optimization pathologies, out-of-distribution brittleness) and their own economics (offline training bills), and they get the full treatment in Module 11: Learning-Enhanced Inversion. Every gradient method in this module runs on the physics model you already trust, with no training data — the one exception being the DiffTactile case study of Section 6, which bolts a trained optical module onto the physics.

Terminology collision — two things called "inverse FEM"

"iFEM" in the shape-sensing literature (Modules 4–6) is the Tessler–Spangler least-squares variational method: a one-shot linear reconstruction that needs no gradients, no material model, no iterations. The gradient-based, PDE-constrained inversion of this module is a different animal that happens to share the name. Conflating them wrecks every comparison: one is a fixed linear map applied at frame rate, the other is an optimization loop whose per-iteration cost is a full FEM solve. When you read "inverse FEM" in a paper, check which one is meant before comparing numbers.

Four routes from one forward model to the same gradient. Finite differences pays per parameter; direct sensitivities pay per parameter but exactly; the adjoint pays one forward solve plus one linear adjoint solve, independent of the parameter count; algorithmic differentiation automates the adjoint — if you differentiate implicitly at the fixed point rather than through the solver iterations.

The module's honest thesis

Differentiability does not repair ill-posedness. An exact gradient descends the same ill-posed landscape a finite-difference gradient stumbles down: noise amplification, null-space non-uniqueness, and zero-sensitivity regions all persist. The regularization theory of Module 3 applies verbatim to every method on this page — $\alpha\,\Phi(\mathbf{p})$ stays in the objective, and experiment design (load cases, sensor placement) stays mandatory. Gradients change what is computable, not what is identifiable.

2. The Adjoint Method: One Extra Solve for the Whole Gradient

We derive the adjoint discrete-first: differentiate the assembled residual, not the continuous PDE. The discrete route produces gradients that are exactly consistent with the discrete objective your optimizer actually sees — a property whose absence causes real line-search failures (see the collapsible below).

Setup and the Matrix Nobody Should Form

The state $\mathbf{u} \in \mathbb{R}^{n_u}$ solves $\mathbf{R}(\mathbf{u},\mathbf{p})=\mathbf{0}$; the objective is $J(\mathbf{u},\mathbf{p})$. The total derivative along the constraint manifold is

\frac{\mathrm{d}J}{\mathrm{d}\mathbf{p}} = \frac{\partial J}{\partial \mathbf{p}} + \frac{\partial J}{\partial \mathbf{u}}\,\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}\mathbf{p}}

Differentiating the constraint $\mathbf{R}(\mathbf{u}(\mathbf{p}),\mathbf{p})=\mathbf{0}$ gives $\dfrac{\partial \mathbf{R}}{\partial \mathbf{u}}\dfrac{\mathrm{d}\mathbf{u}}{\mathrm{d}\mathbf{p}} = -\dfrac{\partial \mathbf{R}}{\partial \mathbf{p}}$, so

\frac{\mathrm{d}\mathbf{u}}{\mathrm{d}\mathbf{p}} = -\left(\frac{\partial \mathbf{R}}{\partial \mathbf{u}}\right)^{\!-1}\frac{\partial \mathbf{R}}{\partial \mathbf{p}}

This is an $n_u \times n_p$ matrix, dense in general — every state's sensitivity to every parameter. (Particular problems can carry structure that thins it: in Section 3's bar, $u_k$ depends only on the elements between the support and node $k$, so the block comes out triangular. Nothing in the general case licenses that.) For a tactile-skin mesh with $n_u \sim 10^5$ and a modulus map $n_p \sim 10^4$, that is $10^9$ entries obtained through $n_p$ linear solves. Nobody should ever form it. Yet it sits inside the gradient:

\frac{\mathrm{d}J}{\mathrm{d}\mathbf{p}} = \frac{\partial J}{\partial \mathbf{p}} \;-\; \underbrace{\frac{\partial J}{\partial \mathbf{u}}\left(\frac{\partial \mathbf{R}}{\partial \mathbf{u}}\right)^{\!-1}}_{\text{a } 1 \times n_u \text{ row!}}\;\frac{\partial \mathbf{R}}{\partial \mathbf{p}}

The Associativity Trick

The whole adjoint method is a choice of parenthesization. Group the product from the left: $\frac{\partial J}{\partial \mathbf{u}}$ is $1 \times n_u$, so $\frac{\partial J}{\partial \mathbf{u}} \big(\frac{\partial \mathbf{R}}{\partial \mathbf{u}}\big)^{-1}$ is a single row vector — one linear solve with the transposed Jacobian. Name it: define the adjoint state $\boldsymbol{\lambda}$ by the adjoint equation, and the gradient collapses:

The Discrete Adjoint Method (KEY EQUATIONS)

\frac{\mathrm{d}J}{\mathrm{d}\mathbf{p}} = \frac{\partial J}{\partial \mathbf{p}} + \boldsymbol{\lambda}^{\mathsf{T}}\,\frac{\partial \mathbf{R}}{\partial \mathbf{p}}, \qquad \left(\frac{\partial \mathbf{R}}{\partial \mathbf{u}}\right)^{\!\mathsf{T}} \boldsymbol{\lambda} = -\left(\frac{\partial J}{\partial \mathbf{u}}\right)^{\!\mathsf{T}}

Recipe: (1) solve the forward problem $\mathbf{R}(\mathbf{u},\mathbf{p})=\mathbf{0}$ — possibly nonlinear; (2) solve one linear adjoint system with the transposed Jacobian; (3) assemble the gradient by cheap products $\boldsymbol{\lambda}^{\mathsf{T}}(\partial\mathbf{R}/\partial p_i)$, one per parameter. Total: one forward solve + one linear solve, independent of $n_p$.

Specialization to Linear FEM

Take $J = \tfrac12\lVert\mathbf{S}\mathbf{u}-\mathbf{d}\rVert^2$ and $\mathbf{R} = \mathbf{K}(\mathbf{p})\,\mathbf{u} - \mathbf{f}$. Then $\partial\mathbf{R}/\partial\mathbf{u} = \mathbf{K}$, $\partial J/\partial\mathbf{u} = (\mathbf{S}\mathbf{u}-\mathbf{d})^{\mathsf{T}}\mathbf{S}$, and (with $\mathbf{f}$ independent of $\mathbf{p}$) $\partial\mathbf{R}/\partial p_i = (\partial\mathbf{K}/\partial p_i)\,\mathbf{u}$:

\mathbf{K}(\mathbf{p})\,\mathbf{u} = \mathbf{f}, \qquad \mathbf{K}^{\mathsf{T}}\boldsymbol{\lambda} = -\mathbf{S}^{\mathsf{T}}\big(\mathbf{S}\mathbf{u}-\mathbf{d}\big), \qquad \frac{\mathrm{d}J}{\mathrm{d}p_i} = \boldsymbol{\lambda}^{\mathsf{T}}\,\frac{\partial \mathbf{K}}{\partial p_i}\,\mathbf{u}

Read the adjoint load: it is the negative measurement residual, back-projected to the sensor locations by $\mathbf{S}^{\mathsf{T}}$. The adjoint field $\boldsymbol{\lambda}$ is the displacement of the same structure under that fictitious residual load — a physical picture the worked example of Section 3 exploits elementwise.

Why symmetric $\mathbf{K}$ makes the adjoint nearly free

Self-adjoint elasticity gives $\mathbf{K}^{\mathsf{T}} = \mathbf{K}$, and for a properly constrained (symmetric positive definite) system the Cholesky factorization computed for the forward solve is reused, the adjoint solve using those factors instead of refactorizing. Symmetry alone is not the licence: where the tangent is symmetric but indefinite — Lagrange-multiplier contact, near-buckling tangents — the reusable factorization is $\mathbf{L}\mathbf{D}\mathbf{L}^{\mathsf{T}}$, and the argument is otherwise unchanged. What remains is not nothing: you still assemble the adjoint right-hand side, run two triangular substitutions, and contract $\boldsymbol{\lambda}$ against $\partial\mathbf{r}/\partial\mathbf{p}$ for every parameter. So the honest claim is one extra linear solve with factor reuse, plus gradient assembly — cheap next to a refactorization, and the solve count is independent of the parameter count — though assembling and storing all $n_p$ gradient components is still $O(n_p)$ work, so 'parameter-independent' describes the linear algebra, not the whole gradient. (For a nonlinear forward problem the asymmetry is larger still: the forward solve is a Newton loop, the adjoint solve is a single linear system.) This is also a trap: because $\mathbf{K}^{\mathsf{T}}=\mathbf{K}$, sloppy code that "forgets" the transpose silently works — and the habit breaks the day you meet a non-self-adjoint operator (advection, dynamics with damping, follower loads). Always derive with the transpose explicit, even when you then optimize it away.

An adjoint is a gradient mechanism, not a regularizer

A surprisingly common misreading: "adjoint inversion" is sometimes glossed as if it were a regularized variant of least squares. It is not. The adjoint machinery computes $\nabla_{\mathbf{p}} J$ — nothing more. It adds no damping, no prior, no filtering; if $J$ lacks $\alpha\Phi(\mathbf{p})$, the adjoint faithfully hands you the gradient of an unregularized, ill-posed misfit, and the optimizer will happily march into the noise-fitting regime. Regularization is a modeling decision made in $J$; the adjoint just differentiates whatever you wrote down. Never write "adjoint $\approx$ regularized LS".

Tangent Mode: the Other Direction

The direct (tangent, forward) sensitivity method solves $\big(\partial\mathbf{R}/\partial\mathbf{u}\big)\,\mathrm{d}\mathbf{u}/\mathrm{d}p_i = -\partial\mathbf{R}/\partial p_i$ once per parameter: $n_p$ linearized solves, each reusing the factorization. It is the right choice when $n_p$ is small and the outputs are many (a handful of hyperelastic constants, sensitivities of a full field wanted). The adjoint wins for the inverse-problem shape — many parameters, one scalar misfit. The two modes are exactly the forward and reverse modes of automatic differentiation, decided at the level of linear algebra.

Two Framings You Must Connect

Framing 1: the adjoint is reverse-mode AD at the equation level. Reverse-mode automatic differentiation propagates one cotangent backward through a computation; the adjoint method does precisely this, but treats "solve $\mathbf{R}=\mathbf{0}$" as a single primitive whose backward rule is the transposed linear solve. In geophysics this is the adjoint-state method: the gradient of the misfit is obtained without ever forming the Fréchet-derivative (Jacobian) matrix, at a cost essentially independent of the number of model parameters (Plessix 2006) — the property that makes large-scale seismic inversion tractable at all. In elasticity imaging — recovering a material-property distribution from measured displacement fields — the same construction is the workhorse (Oberai, Gokhale & Feijóo 2003), and Section 3's tactile-skin bar is its smallest instance.

Framing 2: discretize-then-optimize vs optimize-then-discretize. We differentiated the discrete residual (discretize-then-optimize). The alternative derives the continuous adjoint PDE first and discretizes it afterwards. The distinction sounds academic and is not:

Discretize-then-optimize vs optimize-then-discretize — why line searches care

Discretize-then-optimize (DTO): differentiate the assembled residual $\mathbf{R}$ and the discrete $J$. The resulting gradient is exactly the gradient of the number your optimizer computes — consistent to solver tolerance, by construction. This is what the discrete derivation above produces, and what autodiff of a FEM code produces when the solver is differentiated implicitly at its fixed point — unrolled AD differentiates the iterate the solver stopped at instead (Section 4).

Optimize-then-discretize (OTD): derive the adjoint PDE of the continuous problem, then discretize it — possibly with a different mesh, quadrature, or stabilization than the forward problem. The result approximates the continuous gradient but is generally not the gradient of the discrete objective. The mismatch is $O(h^k)$ — small, but a line search does not know that: it compares actual $J$ values against a predicted descent rate from an inconsistent direction, and near convergence (where the true gradient is small) the inconsistency dominates, and the search stalls mysteriously at step lengths of $10^{-8}$.

Consistency is restored only with compatible discretizations of the forward and adjoint problems. Practical rule: if you hand-derive a continuous adjoint, verify it against finite differences of the discrete objective before trusting it inside an optimizer; if the check plateaus at $10^{-3}$ instead of $10^{-7}$, you have an OTD inconsistency (or a loose solver — Section 4 shows how to tell them apart).

Foreshadow for dynamics: for time-dependent problems the adjoint marches backward in time, consuming stored (or checkpointed) forward states in reverse order. That storage bill and its checkpointing remedy are Section 6's problem.

3. Worked Example: Recovering a Stiffness Profile in a 1D Bar

Everything above in eight unknowns. A bar of length $L=1$ and area $A$, fixed at $x=0$, carries an axial load $F$ at $x=1$, discretized into $N=8$ two-node elements with per-element moduli $E_e$ — the unknowns, $n_p = 8$. You measure displacements $d_i$ at a subset $\mathcal{M}$ of nodes and minimize $J = \tfrac12\sum_{i\in\mathcal{M}}(u_i - d_i)^2$.

The worked example: an 8-element bar with unknown per-element moduli, end load $F$, and displacement sensors on a subset of nodes. Elements to the right of the last sensor will turn out to be exactly invisible.

The element stiffness and its parameter derivative — the only ingredients the adjoint recipe needs beyond the solver:

\mathbf{k}_e = \frac{E_e A}{L_e}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \qquad \frac{\partial \mathbf{k}_e}{\partial E_e} = \frac{A}{L_e}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}

$\mathbf{K}$ is tridiagonal; the forward solve is the Thomas algorithm. The three steps:

Forward: solve $\mathbf{K}\mathbf{u} = \mathbf{f}$.
Adjoint: solve $\mathbf{K}\boldsymbol{\lambda} = -\mathbf{S}^{\mathsf{T}}(\mathbf{S}\mathbf{u}-\mathbf{d})$ — $\mathbf{K}$ symmetric, so same solver, same factorization.
Assemble: $\partial\mathbf{K}/\partial E_e$ is the $E$-independent element matrix scattered into global position, so the gradient collapses to a per-element formula:

Per-Element Gradient = Forward Strain × Adjoint Strain

\frac{\mathrm{d}J}{\mathrm{d}E_e} = \frac{A}{L_e}\,\big(u_2^e - u_1^e\big)\big(\lambda_2^e - \lambda_1^e\big) = A\,L_e\;\varepsilon_e[\mathbf{u}]\;\varepsilon_e[\boldsymbol{\lambda}]

The sensitivity of the misfit to an element's stiffness is the product of the forward strain and the adjoint strain in that element — the discrete version of the sensitivity-kernel picture standard in adjoint seismology (Plessix 2006): the forward field says where the load works, the adjoint field says where the sensors look, and only their overlap contributes to the gradient. Identifiability is the stronger, residual-independent statement — a property of the sensitivity matrix $\mathbf{H} = \partial(\mathbf{S}\mathbf{u})/\partial\mathbf{E}$, its rank and its singular values, not of the gradient at the current iterate. The two coincide in the blind tail of the next subsection, whose $\mathbf{H}$ columns vanish identically so that no residual can move it; elsewhere a gradient component that happens to be zero marks a stationary point, not an unidentifiable parameter.

Identifiability, Read Off the Sensitivity Matrix

Under an end load, $u_k$ depends only on the elements between the support and node $k$. Equivalently: the adjoint bar is loaded only at sensor nodes, so elements to the right of the last sensor carry zero adjoint strain — their gradient is exactly zero, not small. No update driven by the data misfit can move them — not a smaller learning rate, not more iterations, not a better optimizer. Only a prior, or a different experiment, puts them back in play; watch the regularizer alone drag them in the interactive below. That conclusion is specific to this load path: engage the tip stop — switched on and actually in contact — and the tail stops being invisible, because the force transmitted through the contact depends on the compliance of the whole bar, so even upstream sensors start responding to tail moduli. Observability is a property of the experiment, not of the element numbering. Likewise any element with near-zero forward strain is only weakly observable regardless of sensing — a small column of $\mathbf{H}$ rather than an identically zero one, so it is noise-dominated rather than strictly unidentifiable. Its current gradient can nonetheless be exactly zero, whenever the residual happens to sit orthogonal to that column: identifiability is read off $\mathbf{H}$, never off a single gradient evaluation. The remedies are experimental, not algorithmic: multiple load cases, sensors spanning the domain, and Tikhonov regularization on $\log E$ to fill the null space with a prior.

Numerical Hygiene, Shown with Numbers

Positivity: parameterize $E_e = \exp(\theta_e)$ and descend in $\theta$. A plain gradient step on raw $E$ can drive a modulus negative and make $\mathbf{K}$ indefinite — the solver then fails in ways that look like optimizer bugs.
Verify the adjoint against central finite differences $\big(J(E_e+h) - J(E_e-h)\big)/(2h)$, sweeping $h$ over decades. The FD error curve is V-shaped — truncation error falls as $h^2$, round-off rises as $h$ shrinks — and at the bottom the match should reach relative error $\sim 10^{-6}$–$10^{-8}$ when the solver tolerance is tight. A loose forward solve caps the achievable agreement: the check plateaus at a floor set by the solver's inexactness, not by the tolerance number itself, and no amount of $h$-tuning fixes it. How the residual maps onto that floor is problem-dependent — the state error passes through the inverse tangent and then through the objective's own derivatives, so the conditioning of the tangent bounds the amplification without predicting it, and cancellation can leave the realized error far below the bound. In the laboratory of Section 4 the floor lands above the residual for some settings and below it for others, so read the plateau as a symptom of the loose solve, never as a measurement of it.

Both phenomena — the exactly-zero gradients and the V-curve — are live in the two interactives below.

Interactive Tool — Adjoint Bar Inverter

Run the full inversion on the 8-element bar: each gradient-descent step costs one forward + one adjoint solve. Watch the per-element gradient equal the product of forward and adjoint strain, watch elements the sensors cannot see barely move (hatched) — what little motion they do show comes from the smoothness prior, not from data, and use noise + $\alpha$ to replay the bias–variance trade-off of Module 3 in miniature. Hover an element for its live numbers.

Sensors $M$ 5

Pattern

Noise $\sigma$ (% of max $u$) 1.0

$\log_{10}\alpha$ (Tikhonov on $\log E$) −4

Step size $\eta$ 8

4. Five Route Families to the Same Gradient

Every route below converges (in exact arithmetic, at a converged forward solve) to the same vector $\nabla_{\mathbf{p}}J$. What differs is the bill: how many solves, how much memory, and how the answer degrades when the forward solve is not actually converged. Choosing a route is an engineering decision, not a mathematical one — and the two "automatic" routes are not interchangeable, which is the single most common confusion in differentiable simulation.

Route	Cost per gradient	Exactness	Memory	Breaks down when
FD one-sided	$(n_p+1)\,C_{\text{solve}}$	$O(h)$ truncation; best $\sim\sqrt{\varepsilon_{\text{mach}}}$	$O(n_u)$	$h$ ill-chosen; $n_p$ large; $J$ noisy or loosely solved
FD central	$2 n_p\,C_{\text{solve}}$	$O(h^2)$; heuristic $h^\star \sim \sqrt[3]{\varepsilon_{\text{mach}}}\,\max(\|p_i\|,s_i)$, $s_i>0$	$O(n_u)$	same, at double the price
Direct (tangent)	1 forward + $n_p$ linearized solves (factorization reused)	exact to solver tolerance	$O(n_u)$	$n_p$ large
Adjoint	1 forward + 1 transposed linear solve	exact to solver tolerance	$O(n_u)$ static	you must code $\partial\mathbf{R}/\partial\mathbf{p}$ and the transpose by hand
AD — unrolled	1 forward, backprop through every solver iteration	gradient of the unconverged iterate, not of the solution	$O(\text{iters}\cdot n_u)$ tape	many iterations; tape memory; iterate far from converged
AD — implicit	1 forward + 1 adjoint solve (registered as the solver's VJP)	exact for the converged discrete solution	$O(n_u)$	$\mathbf{R}\neq\mathbf{0}$ at the "solution" — loose tolerance corrupts it

Key insight — unrolled and implicit differentiation are different algorithms

Both are "autodiff of the solver", and papers say "differentiable FEM" for either. Unrolling records every Newton/CG iteration on the tape and differentiates the iteration path: memory grows with iteration count, and the result is the exact gradient of whatever the iterate happens to be after $k$ steps — which is the wrong object if $k$ was too small. Implicit differentiation ignores the path entirely and differentiates the fixed-point condition $\mathbf{R}(\mathbf{u}^*,\mathbf{p})=\mathbf{0}$: memory is independent of iteration count and the gradient is exact for the converged solution — but meaningless if the solve never converged. Same marketing term, different memory bills, different failure modes. Section 5 gives the implicit route its proper form.

Tighten the solve before blaming the math

Implicit differentiation (and the adjoint, which is its hand-written form) assumes $\mathbf{R}(\mathbf{u}^*,\mathbf{p})=\mathbf{0}$ holds exactly. If the forward solve stops at a residual of $10^{-4}$, the gradient inherits an error of that origin — but not of that size: the state error carries a factor $(\partial\mathbf{R}/\partial\mathbf{u})^{-1}$ and is then contracted against the objective's own derivatives, so the finite-difference plateau is the residual filtered by problem-dependent constants that the conditioning of the tangent bounds rather than predicts. Expect a floor unequal to $10^{-4}$ — it can land above it or, through cancellation, well below — and no choice of $h$ removes it. The diagnostic is free but not conclusive: an early, flat bottom in the V-curve says only that the error is $h$-independent, and a loose solve is one of several ways to produce that — a wrong $\partial\mathbf{R}/\partial\mathbf{p}$ term, a dropped sign, a missing transpose, or an OTD inconsistency (Section 2) all plateau the same way. The cheap disambiguation is to tighten the forward solve and re-run the sweep: if the floor drops, the solver was the limit; if it refuses to move, the gradient code is. Unconverged solves give unreliable gradients on every route: FD differentiates solver noise, unrolling differentiates a wrong iterate, implicit differentiation differentiates a violated identity.

The laboratory below runs these route families on the Section 3 bar — same data, same unknowns — and puts the bills side by side. Its menu has six entries for five families, and the sixth is not a sixth implementation: the implicit-AD entry runs the identical code as the hand-written adjoint, because that is the claim — a framework's custom_vjp for the solve computes what Section 2 derives by hand, and there is nothing left for a separate implementation to do. The genuinely distinct code paths are the two FD stencils, the tangent sweep, the adjoint, and the taped unrolling. The two panels at the bottom hold the module's two cautionary plots: the FD step-size V-curve (with a loose-solver toggle), and the kink that a contact active-set change leaves in an otherwise smooth objective.

Interactive Tool — Gradient-Routes Laboratory

One inverse problem, five gradient families on a six-entry menu — implicit AD is the adjoint under its autodiff name and executes the same code, by design. Pick a route and compare its per-element gradient (orange) against the exact discrete gradient (purple), with the running cost meter. Then stress it: loosen the forward solve, shrink the unrolling budget, or push the bar tip into a contact stop and watch the objective grow a kink.

Gradient route

FD step $\log_{10} h$ −6

Unroll iterations $k$ 150

Sensors $M$ (nodes 1…M) 5

Noise $\sigma$ (%) 1.0

Contact gap $g$ 0.90

loose forward solve (all routes except unrolled) contact stop at tip

FD verification V-curve (element $e = \min(3, M)$)

Objective sweep $J(s)$, $\mathbf{E} = s\,\hat{\mathbf{E}}$ — the contact kink

5. Differentiable FEM: Autodiff Meets the Stiffness Matrix

Definition — Differentiable Simulator

A simulator that exposes gradients of its outputs with respect to any of its inputs — material parameters, loads, nodal coordinates (shape), initial conditions — so the whole simulation can sit inside a gradient-based optimization or learning loop as if it were one differentiable function.

Section 2 built the adjoint by hand. Differentiable FEM is the industrialization of the same mathematics: let automatic differentiation assemble the derivative terms, and reserve your own effort for the one decision that matters — how the solve is differentiated. The two strategies from Section 4, now in their proper form:

Strategy 1 — Unrolling

Record every operation of every Newton or CG iteration on the autodiff tape and backpropagate through the lot. Correct in the limit, but memory grows with iteration count, and what you differentiate is the $k$-th iterate — an object that changes when the iteration count changes. For stiff problems needing hundreds of iterations, the tape is the dominant memory cost of the whole pipeline.

Strategy 2 — Implicit Differentiation at the Fixed Point

Apply the implicit function theorem to the converged equilibrium:

\mathbf{R}(\mathbf{u}^*,\mathbf{p}) = \mathbf{0} \;\Longrightarrow\; \frac{\mathrm{d}\mathbf{u}^*}{\mathrm{d}\mathbf{p}} = -\left(\frac{\partial \mathbf{R}}{\partial \mathbf{u}}\right)^{\!-1}\frac{\partial \mathbf{R}}{\partial \mathbf{p}}\;\Bigg|_{\mathbf{u}^*}

In reverse mode, the vector–Jacobian product the framework actually needs is exactly Section 2's adjoint solve with the transposed tangent stiffness:

\mathbf{v}^{\mathsf{T}}\,\frac{\mathrm{d}\mathbf{u}^*}{\mathrm{d}\mathbf{p}} = \boldsymbol{\lambda}^{\mathsf{T}}\,\frac{\partial \mathbf{R}}{\partial \mathbf{p}}, \qquad \left(\frac{\partial \mathbf{R}}{\partial \mathbf{u}}\right)^{\!\mathsf{T}}\boldsymbol{\lambda} = -\mathbf{v} \quad \text{(VJP = adjoint solve)}

The sign convention is Section 2's, deliberately: set $\mathbf{v} = (\partial J/\partial\mathbf{u})^{\mathsf{T}}$ and this $\boldsymbol{\lambda}$ is Section 2's $\boldsymbol{\lambda}$, and the VJP is the $\boldsymbol{\lambda}^{\mathsf{T}}\,\partial\mathbf{R}/\partial\mathbf{p}$ term of the gradient. (Autodiff libraries often define the seed without the minus sign, which flips $\boldsymbol{\lambda}$; either is fine, but carrying one symbol under two conventions across a derivation is how sign errors are born.) Memory is independent of how many iterations the forward solve took, and the gradient is exact for the converged discrete solution. In JAX this is the custom_vjp pattern: register the adjoint solve as the backward rule of the solver call, and the tape never sees the iterations. The circle closes: differentiable FEM done right is the adjoint method, automated.

What autodiff actually buys you

Not the adjoint equation — that is three lines of linear algebra. What autodiff automates is $\partial\mathbf{R}/\partial\mathbf{p}$: the derivative of the assembled residual with respect to any parameterization — a heterogeneous material field, shape derivatives via nodal coordinates, load parameters — the part that is genuinely tedious and error-prone by hand (every constitutive branch, every quadrature loop). Change the parameterization and the gradient machinery follows for free; a hand-written adjoint would need a new $\partial\mathbf{R}/\partial\mathbf{p}$ derivation each time.

Ecosystem Landmarks (verified numbers only)

DiffTaichi (Hu et al., ICLR 2020): source-code-transformation AD inside a physics DSL. Across 10 differentiable simulators, their differentiable elastic-object simulator is 4.2× shorter than the hand-engineered CUDA reference, runs as fast, and is 188× faster than a TensorFlow implementation; controllers optimize in tens of gradient iterations.
DiffPD (Du et al., ACM TOG 2021): differentiates implicit-Euler projective dynamics and reuses the prefactorized Cholesky decomposition of the PD global matrix in the backward pass — the "symmetric $\mathbf{K}$ makes the adjoint nearly free" insight, scaled up to soft-body dynamics. Reported 4–19× speedup over Newton-based differentiable-simulator baselines; applications in system identification, trajectory optimization, and closed-loop soft-robot control.
JAX-FEM (Xue et al., Computer Physics Communications 2023): open-source differentiable GPU-accelerated 3D FEM where "inverse problems are solved in a fully automatic manner without the need to manually derive sensitivities." On a 3D tensile problem with 7.7 million DOF it reports ~10× GPU acceleration over a commercial FEM code, and demonstrates topology optimization of nonlinear materials.
NVIDIA Warp (open-source): JIT-compiles Python kernels for CPU/GPU; kernels are differentiable and interoperate with PyTorch, JAX, and Paddle — the infrastructure layer for custom differentiable physics.

The idea has since reached the two communities this guide cares about. On the engineering-structures side, Wang et al. (CMAME 2025) present a differentiable finite element method with Galerkin discretization aimed at fast and accurate inverse analysis of multidimensional heterogeneous structures — differentiable FEM as a general inverse-analysis tool rather than a graphics technique. On the soft-robotics side, Zhong et al. (IEEE RA-L 2025) build an end-to-end modeling framework for pneumatic soft robots on differentiable FEM — the same solver family, pointed at exactly the compliant-body regime a tactile skin lives in.

6. Dynamics, Contact, and the Quality of Gradients

Static inversion needs one adjoint solve. Dynamic and contact-rich problems — a finger sliding across a tactile skin, a soft gripper closing — raise two genuinely new issues: where do the backward pass's states come from, and can the gradient be trusted near non-smooth events?

Memory: the Backward March

The time-dependent adjoint marches backward from the final time, and at step $t$ it needs the forward state $\mathbf{u}_t$ — which the forward pass computed long ago. Storing all $T$ steps costs $\mathcal{O}(T\,n_u)$; checkpointing trades recomputation for memory: store every $k$-th state, recompute each segment on demand during the backward sweep.

\text{memory: } \mathcal{O}(T\,n_u) \;\xrightarrow{\text{checkpoint every } k} \; \mathcal{O}\big(\big(\tfrac{T}{k} + k\big)\,n_u\big) + \text{recompute } \mathcal{O}(k) \text{ steps per segment}

Both terms are real. The $T/k$ term is the checkpoints themselves; the $+k$ term is the segment currently being replayed, which must be held in memory while the backward sweep walks it in reverse — a recomputed segment cannot be traversed backwards for free. The sum is minimized at $k \sim \sqrt{T}$, giving $\mathcal{O}(\sqrt{T}\,n_u)$: a square-root memory bill for a factor-of-two recompute. That is the result that makes checkpointing worth the bookkeeping.

Gradient Quality: Contact Is the Hard Part

Contact and friction are non-smooth physics. Penalty-based contact replaces the complementarity condition with stiff, piecewise-differentiable forces — the usual unilateral spring $k_c\max(0,\,\text{penetration})$ is smooth away from activation but its slope jumps at make/break, unless the penalty is explicitly smoothed (a softplus, or a quadratic transition through the gap function); complementarity/LCP formulations are more faithful to the physics but also only piecewise-smooth. Either way, contact make/break events and stick–slip transitions make the state trajectory non-smooth in the parameters: gradients near events can be discontinuous, biased, or high-variance. That trade-off — between versatility, computational speed, and gradient accuracy — is the axis along which the IEEE Access review of differentiable simulators (Newbury et al. 2024) organizes the field's design decisions. Its own account of the contact axis is narrower than the folklore, and worth quoting precisely: compliant (soft) models relax the step-like non-penetration and friction laws into functions with finite gradients, so derivatives can be taken directly, at the cost that the constraint is no longer strictly enforced — objects interpenetrate. Whether such simplifications help or hurt the quality of the resulting gradients the review explicitly leaves open: "it is currently not clear how these simplifications affect the differentiable simulation performance, both in terms of the quality of generated gradients and how well the simulated behavior translates to the real world" (Newbury et al. 2024, §V-A). Treat any confident soft-vs-hard gradient-quality ranking you read elsewhere as a hypothesis, not a result.

A related structural non-smoothness comes in two severities that are worth keeping apart. Remeshing genuinely swaps the discrete unknowns mid-optimization, so a difference quotient taken across it compares objectives defined on different discrete systems. The complete adaptive algorithm does still define some scalar map $\mathbf{p}\mapsto J(\mathbf{p})$ — but one carrying remeshing discontinuities, so the quotient approximates no fixed-discretization derivative and need not converge as $h\to0$. An active-set change is milder: the unknown vector persists and only the constraint equations switch branch, so one-sided derivatives exist on each side and the adjoint gradient is exact on each branch. Only a stencil that straddles the event returns garbage.

Piecewise-smooth is not a bug in your code

The kink panel in the laboratory above shows the mechanism at toy scale: when the bar tip meets the stop, the active set changes, the slope of $\mathbf{u}(s)$ changes, and $J(s)$ has a corner — the gradient jumps. On each side of the event the adjoint gradient is exact; at the event it does not exist, and a finite-difference stencil that straddles the event returns garbage that no step size fixes. In a full simulator these events are dense in time (across a tactile skin, contact nodes enter and leave the active set continually), so dynamic gradients are best treated as noisy estimates of a piecewise object. Gradient accuracy is one of the three axes — with versatility and computational speed — along which differentiable-simulator design decisions trade off (Newbury et al. 2024); a kink at an active-set change is that trade-off showing up in your plot, not a defect of your implementation.

Case Study: DiffTactile

DiffTactile (Si et al., ICLR 2024) is a physics-based differentiable tactile simulator for contact-rich robotic manipulation — and a template for this guide's readers. It models the sensing elastomer with an FEM-based soft-body model, uses a penalty-based contact model, and adds a pixel-based neural module to infer the sensor's optical response to contact. Its gradients serve two distinct masters: optimizing tactile-assisted grasping and manipulation skills, and refining the simulator's physical properties against real-world data — gradient-based system identification that narrows the sim-to-real gap. That second use is the pattern to steal for a lab-built vision-based skin: calibrate the simulated elastomer against real indentation data by gradient descent, and only then trust the simulator's gradients for force inversion and design.

Practical Guidance for Dynamic Gradients

Keep time steps implicit — stability under stiff penalty forces, and the implicit solve is exactly what implicit differentiation wants.
Average or smooth losses over windows rather than pointwise states near contact events; a windowed loss integrates over the kink instead of sitting on it.
Validate dynamic gradients by finite differences on short horizons before trusting long rollouts — the V-curve discipline of Section 3, applied per time-window.

7. Practical Recipes for a Soft Tactile Skin

The decision procedure, for the reader whose measurements are camera-observed marker displacements on a silicone skin and whose unknowns are a contact traction, a modulus map, or the skin geometry.

The gradient-based inversion loop for a vision-based tactile skin, with the classical one-shot route shown above it for contrast. The reverse arrow is the whole module: one adjoint solve per iteration prices the full parameter gradient.

Establish the linear baseline first. Classical iFEM / Tikhonov reconstruction from the earlier modules is the sanity check every fancy method must beat. If the linear regime genuinely holds, stop there — nothing below is faster than a prefactorized one-shot solve.
In the nonlinear regime, move to differentiable FEM or a hand adjoint. Silicone at 10–50% strain with $E$ of tens to hundreds of kPa is hyperelastic territory. Before optimizing anything, rescale to units in which the numbers sit near unity — kPa and mm are usually a good start, though they are still dimensional and guarantee nothing by themselves. Where parameters and residual still span decades in those units, nondimensionalize properly: divide displacements by a characteristic length $L_0$, moduli by a reference modulus $E_0$, and the residual by its own characteristic scale, so both the parameter vector and the residual come out $O(1)$.
Always verify gradients. Central differences on a handful of parameters, $h$ swept over decades, V-shaped error curve, minimum relative error $\sim 10^{-6}$–$10^{-8}$: $\frac{\mathrm{d}J}{\mathrm{d}p_i}\bigg|_{\text{FD}} = \frac{J(p_i+h)-J(p_i-h)}{2h} + \mathcal{O}(h^2), \qquad h^{\star} \sim \sqrt[3]{\varepsilon_{\text{mach}}}\,\max\big(|p_i|,\,s_i\big), \quad s_i > 0$ The step rule is a heuristic, and the nonzero scale $s_i$ is not decoration: a centred or log-parameterized coefficient routinely sits at $p_i = 0$ — the bar inverter below resets to exactly $\theta = \mathbf{0}$ — where $\sqrt[3]{\varepsilon_{\text{mach}}}\,|p_i|$ prescribes $h = 0$ and the quotient is undefined. Set $s_i$ to a typical magnitude for that parameter. The true optimum also depends on the size of $J$ and of its third derivative, which the rule does not know. If the curve bottoms out early and flat, the error is $h$-independent and the step size is not what is capping you. Tighten the forward solve first — implicit differentiation assumes $\mathbf{R}=\mathbf{0}$ holds — and re-run the sweep: a floor that drops was the solver, a floor that does not move is a wrong adjoint term or an OTD inconsistency.
Parameterize for positivity. $E = \exp(\theta)$ or softplus; a plain gradient step on raw $E$ can produce a negative modulus and an indefinite $\mathbf{K}$.
Regularize deliberately. Tikhonov on the log-field for smooth modulus maps, total variation for inclusions and defects, and multiple load cases (different indentation sites) to shrink the null space — recall the exactly-zero gradients of Section 3: $E(\mathbf{x}) = \exp\big(\theta(\mathbf{x})\big), \qquad \Phi(\theta) = \tfrac12\lVert \nabla \theta \rVert_2^2 \;\;\text{or}\;\; \lVert \nabla \theta \rVert_1 \;\;\text{(TV)}$ Both carry the single weight $\alpha$ of the objective $J = \tfrac12\lVert\mathbf{S}\mathbf{u}-\mathbf{d}\rVert^2 + \alpha\,\Phi$ — the same $\alpha$ the bar inverter's slider sets. Putting a second weight inside $\Phi$ only gives you two knobs for one degree of freedom.
Optimizer: L-BFGS with strong-Wolfe line search for smooth deterministic losses; Adam only when the loss is stochastic (minibatched loads or frames).
System identification first, inversion second. Calibrate the elastomer parameters against controlled real indentation data by gradient descent — the DiffTactile pattern (Si et al. 2024) — before trusting simulator gradients for force inversion or for optimizing skin geometry and marker layout.
If you need real-time inference, or the physics is partly unknown, the answer is not on this page: amortized neural operators and per-instance physics-informed training are Module 11's subject, where they are weighed honestly against the adjoint loop built here — which remains the certifying fallback that checks any learned answer against the actual physics.

Tool	What it is	Differentiation strategy
JAX-FEM	Open-source differentiable 3D FEM (Xue et al. 2023); automatic sensitivities, inverse design, topology optimization of nonlinear materials	implicit / custom VJP in JAX
DiffTaichi / Taichi	High-performance physics DSL with built-in AD (Hu et al. 2020); custom kernels	source-code transformation
NVIDIA Warp	Differentiable JIT-compiled Python GPU kernels; PyTorch/JAX/Paddle interop	kernel-level AD
DiffPD	Soft-body projective dynamics with fast implicit-integration gradients (Du et al. 2021)	adjoint with prefactorized Cholesky

Closing the module

Gradients change what is computable, not what is identifiable. The adjoint made a million-parameter gradient cost two solves; it did not add one bit of information to the data. Regularization and experiment design — load cases, sensor placement, the priors of Module 3 — remain the backbone of every inversion this module's machinery accelerates.

Flashcards

References

Plessix, R.-E. (2006). A review of the adjoint-state method for computing the gradient of a functional with geophysical applications. Geophysical Journal International. doi:10.1111/j.1365-246X.2006.02978.x
Oberai, Assad A.; Gokhale, Nachiket H.; Feijóo, Gonzalo R. (2003). Solution of inverse problems in elasticity imaging using the adjoint method. Inverse Problems 19(2), pp. 297–313. doi:10.1088/0266-5611/19/2/304
Hu, Y.; Anderson, L.; Li, T.-M.; Sun, Q.; Carr, N.; Ragan-Kelley, J.; Durand, F. (2020). DiffTaichi: Differentiable Programming for Physical Simulation. ICLR 2020.
Du, T.; Wu, K.; Ma, P.; Wah, S.; Spielberg, A.; Rus, D.; Matusik, W. (2021). DiffPD: Differentiable Projective Dynamics. ACM Transactions on Graphics.
Xue, Tianju; Liao, Shuheng; Gan, Zhengtao; Park, Chanwook; Xie, Xiaoyu; Liu, Wing Kam; Cao, Jian (2023). JAX-FEM: A differentiable GPU-accelerated 3D finite element solver for automatic inverse design and mechanistic data science. Computer Physics Communications 291, 108802. doi:10.1016/j.cpc.2023.108802
Wang, Xi; Yin, Zhen-Yu; Wu, Wei; Zhu, He-Hua (2025). Differentiable finite element method with Galerkin discretization for fast and accurate inverse analysis of multidimensional heterogeneous engineering structures. Computer Methods in Applied Mechanics and Engineering 437, 117755. doi:10.1016/j.cma.2025.117755
Zhong, Shaohong; Yao, Yao; Maiolino, Perla; Posner, Ingmar (2025). An End-to-End Framework for Modeling Pneumatic Soft Robots Based on Differentiable Finite Element Methods. IEEE Robotics and Automation Letters 10(12), pp. 13003–13010. doi:10.1109/LRA.2025.3625507
Si, Z.; Zhang, G.; Ben, Q.; Romero, B.; Xian, Z.; Liu, C.; Gan, C. (2024). DIFFTACTILE: A Physics-based Differentiable Tactile Simulator for Contact-rich Robotic Manipulation. ICLR 2024.
Newbury, R.; Collins, J.; He, K.; Pan, J.; Posner, I.; Howard, D.; Cosgun, A. (2024). A Review of Differentiable Simulators. IEEE Access.

← 9. Reference Configuration & Prestress 11. Learning-Enhanced Inversion →