8. Load & Contact Force Reconstruction

Transfer matrices, ill-conditioning & SVD, Tikhonov-regularized inversion, L-curve & GCV, time-domain deconvolution, GelSlim iFEM case study, soft robot skins as distributed force sensors

Contents

1. Every Elastic Body Is a Force Sensor 2. The Transfer Matrix: Condensing an FE Model into a Sensor Model 3. Why Naive Inversion Explodes: SVD Applied to $\mathbf{H}$ 4. Tikhonov-Regularized Inversion 5. Choosing $\lambda$: L-Curve, GCV, Discrepancy 6. Time-Domain Force Identification 7. From Structures to Skin: the Tactile Inverse Problem 8. Case Study: GelSlim 2.0 and Inverse FEM 9. Soft Robot Skins as Distributed Force Sensors: Design as Conditioning 10. Practical Recipe and Pitfalls Interactive: Contact Inversion Lab Interactive: L-Curve & GCV Explorer Interactive: GelSlim Pipeline Walkthrough Flashcards

1. Every Elastic Body Is a Force Sensor

An elastic body in equilibrium encodes its applied loads in its deformation. Push on a wing, a bridge deck, or a silicone fingertip pad, and the displacement field that results is a complete — if scrambled — record of the pushing. This module is about unscrambling it: measure the deformation, know the elasticity, recover the load. It is the third of the four inverse families in this guide's taxonomy, and it behaves very differently from the shape sensing of Modules 4–6.

Core Problem

Given $m$ noisy displacement or strain measurements on an elastic body whose FE model is known, reconstruct the applied load — ideally the full traction distribution $\mathbf{t}(\mathbf{x})$ on the candidate contact surface, not merely its resultant. In the guide's observation model $y = \mathcal{H}(u, p, f, X_0) + b + \eta$: geometry $X_0$ and material $p$ are assumed known, the response $u$ is (partially) measured, and the unknown is the load $f$.

	In this module
Measured	Displacements at $m$ points (marker tracking, DIC) or strains (gauges), noisy
Unknown	$n$ load parameters $\mathbf{g}$: nodal contact forces, pressure-patch amplitudes, or a force time history $g(t)$
Assumed known	Geometry, mesh, $E$, $\nu$, boundary conditions, (for dynamics) $\mathbf{M}, \mathbf{C}$ — the entire forward model
Unobservable	Load components in the nullspace of the observation map, and spatially oscillatory traction patterns whose displacement signature is below the noise floor (Section 3)

The forward problem — given the boundary tractions, solve the elastostatic boundary-value problem for the displacements — is well-posed: the solution exists, is unique, and depends continuously on the data. The inverse load problem fails Hadamard's stability requirement (Module 3): it is the discretization of a first-kind Fredholm integral equation, the textbook setting for ill-posedness. How severe depends on the kernel: a smooth or depth-separated kernel (a load buried below the sensing plane, a thick compliant layer) smooths strongly and gives near-exponential singular-value decay, whereas the same-surface elastic kernel used later in this module is only weakly singular — Boussinesq goes like $1/r$ — and damps high spatial frequencies algebraically rather than exponentially. Both are ill-posed; they are not equally hostile, and design claims should be made against the actual kernel. Small noise in the measured displacements maps to enormous errors in the reconstructed load.

Forward vs. inverse. Elasticity destroys high-frequency load information on the way from traction to displacement ($\mathbf{K}^{-1}$ smooths); naive inversion $\mathbf{H}^{+}$ tries to re-create that information from noise and produces oscillatory garbage; Tikhonov-regularized inversion declines to invent it and returns the recoverable part.

A Decade Older Than Tactile Robotics

Load identification from measured response is a mature discipline spanning scales. Aircraft flight loads are monitored from wing strain gauges. Maincon et al. (OMAE 2008) estimated vortex-induced-vibration forces on marine risers via inverse FEM — the term "inverse FEM" for force estimation predates its tactile-robotics use by a decade. Vehicle and bridge dynamic load identification has its own literature, systematized in the reviews of Sanchez & Benaroya (2014) and Liu et al. (2022). This module's destination sits at the opposite end of the size scale: a robot fingertip whose millimeters-thick gel pad deforms under contact while a camera watches from behind.

Why Inversion Is Hard: Elasticity Smooths

The physical root of the difficulty is that elasticity is a smoothing operator. The classical Boussinesq solution makes this concrete: a normal point load $P$ on an elastic half-space produces the surface deflection

u_z(r) = \frac{(1-\nu^2)\,P}{\pi E\, r}

Read it two ways. First, the displacement footprint of a point load is wide: $1/r$ decay, no compact support — a delta function in force becomes a long-tailed bump in displacement. Second, and consequently, two nearby point loads produce almost identical displacement fields everywhere except very close to the loads. Distinguishing them from noisy data requires amplifying tiny differences — and amplifying tiny differences amplifies noise with them.

Key insight — force is a global functional

Ma et al. (ICRA 2019, Fig. 1) observed on a GelSlim sensor that markers move substantially outside the contact patch. That motion is caused by internal elastic stresses transmitted through the gel, not by any external load at those markers. The tempting per-point assumption $f_i \propto u_i$ — a diagonal stiffness — therefore produces physically wrong force maps: it paints force wherever there is motion. The force at a node depends on a coupled neighbourhood of displacements, never on that node's displacement alone: the assembled $\mathbf{K}$ is sparse but couples each node to every node it shares an element with. And the moment you invert or regularize rather than multiply, the dependence goes global — $\mathbf{K}^{-1}$, and hence $\mathbf{H} = \mathbf{S}\mathbf{K}^{-1}\mathbf{B}$, is dense. Every correct method in this module is, at bottom, machinery for honoring that coupling.

Top view of a tactile pad, redrawing the key figure logic of Ma et al. (2019). Left: the measured marker displacement field radiates far beyond the true contact patch (dashed) — driven by internal elastic stresses. Right: the reconstructed force field is concentrated in the patch, with small noise-driven vectors surviving outside it (Ma et al.'s own Fig. 9 caption reads "naturally limited in contact patch in most area, despite globally existing noises"), and near the patch boundary the force can point opposite to the local displacement. $f_i \propto u_i$ is a fallacy.

Warning — two different "inverse FEMs"

This guide contains two unrelated methods that both go by "iFEM". (1) Shape-sensing iFEM (Tessler–Spangler, Modules 4–6): measured strains $\to$ displacement field, via a variational least-squares functional; loads never appear. (2) Load-reconstruction inverse FEM (this module; Maincon et al. 2008, Ma et al. 2019): measured displacements $\to$ applied forces, via a stiffness or transfer-matrix model. Same acronym, different unknowns, different mathematics. Li et al. (2025) are the paper most often miscounted as a chain of both, so read their pipeline carefully: shape-sensing iFEM reconstructs a displacement field from measured strains, radial basis functions fit that discrete field into a continuous one, and the distributed load then falls out of the thin-plate equation of motion. There is no second, load-reconstruction iFEM in it — no transfer matrix is ever formed and inverted, which is what the authors mean when they say the method prevents ill-posedness. Do not read that as instability removed: the plate equation reaches the load through fourth spatial derivatives of the displacement field, and it is the least-squares iFEM step and the RBF fit that smooth the field first. The regularization is implicit and upstream, not absent. Always say which iFEM you mean.

The route through the module: build the linear observation model (Section 2), see exactly why naive inversion explodes (Section 3), repair it with Tikhonov regularization (Section 4) and choose the parameter honestly (Section 5), add dynamics (Section 6), then cross from structures to skin (Sections 7–9) and close with a field-tested recipe (Section 10).

2. The Transfer Matrix: Condensing an FE Model into a Sensor Model

Everything in the static linear theory follows from three matrices of bookkeeping wrapped around one FE solve. Start from the constrained static FE model of Module 2:

\mathbf{K}\mathbf{u} = \mathbf{f},\qquad \mathbf{f} = \mathbf{B}\mathbf{g},\qquad \mathbf{y} = \mathbf{S}\mathbf{u}

$\mathbf{K} \equiv \mathbf{K}_{ff}$ — the $N \times N$ free-DOF block of the stiffness matrix after boundary conditions, symmetric positive definite ($N$ = all free DOFs). Keep the unreduced $\mathbf{K}_{\mathrm{full}}$ around too: it is the one that still has rows at the constrained DOFs, and those rows are what return support reactions (trap 1 below);
$\mathbf{B}$ ($N \times n$) — the candidate-load parametrization: it maps $n$ unknown force parameters $\mathbf{g}$ (nodal forces on the candidate contact surface, or amplitudes of pressure patches) into the global load vector;
$\mathbf{S}$ ($m \times N$) — the observation operator: it selects the measured DOFs, or extracts strains via rows of the strain-displacement matrix for gauge-based sensing.

Eliminate $\mathbf{u}$ and the sensor model appears:

\mathbf{y} = \underbrace{\mathbf{S}\,\mathbf{K}^{-1}\mathbf{B}}_{\;\mathbf{H}\,\in\,\mathbb{R}^{m\times n}}\,\mathbf{g}

$\mathbf{H}$ is the transfer matrix (influence-coefficient matrix), whose units follow the bases you chose: sensor output per unit force when the columns are nodal forces (m/N for displacement sensors, strain-per-newton for strain sensors), or sensor output per unit pressure when the columns are pressure-basis amplitudes. Entry $H_{ij}$ is the response of sensor $i$ to a unit value of the load-basis coefficient $g_j$ — a unit nodal force for a nodal-force basis, a unit pressure amplitude for a pressure basis. One column of $\mathbf{H}$ per candidate load, one row per sensor.

Assembling $\mathbf{H}$ column by column. Left: candidate force nodes (open circles) and sensors (filled squares) on the meshed contact surface; a unit load at candidate node $j$. Middle: one forward FE solve. Right: the measured response at the $m$ sensors fills column $j$ of the $m \times n$ transfer matrix. $n$ forward solves, offline — or, by Maxwell–Betti, $m$ solves with unit loads at the sensors.

Practical Assembly

You never form $\mathbf{K}^{-1}$ (Module 2's rule: solve the system). Factorize $\mathbf{K}$ once, then perform $n$ forward solves with unit loads — one per candidate DOF — sampling each solution at the sensors. All offline. When $m < n$ (few sensors, many candidate loads), Maxwell–Betti reciprocity cuts the work: the displacement at $i$ due to a unit force at $j$ equals the displacement at $j$ due to a unit force at $i$ (for collocated force/displacement pairs), so you can assemble $\mathbf{H}$ row-wise from only $m$ solves with unit loads applied at the sensor DOFs.

With $m \ge n$ and $\mathbf{H}$ of full column rank — at least as many measurement channels as unknowns is necessary, never sufficient; redundant or badly placed sensors leave the rank deficient even at $m \gg n$ — the least-squares estimate is

\hat{\mathbf{g}}_{\mathrm{LS}} = \left(\mathbf{H}^{\mathsf T}\mathbf{H}\right)^{-1}\mathbf{H}^{\mathsf T}\mathbf{y}

— or, for computation, $\hat{\mathbf{g}}_{\mathrm{LS}} = \mathbf{H}^{+}\mathbf{y}$ — algebraically fine, statistically a disaster, for reasons Section 3 makes quantitative. And never form $\mathbf{H}^{\mathsf T}\mathbf{H}$ to get there: $\kappa_2(\mathbf{H}^{\mathsf T}\mathbf{H}) = \kappa_2(\mathbf{H})^{2}$, so the normal equations square an already brutal condition number and wipe out exactly the small singular directions you need to see. Use a rectangular QR or SVD of $\mathbf{H}$ itself. For harmonic and dynamic loading the same construction produces the frequency response function (FRF) matrix used in frequency-domain identification:

\mathbf{H}(\omega) = \mathbf{S}\left(\mathbf{K} - \omega^{2}\mathbf{M} + \mathrm{i}\,\omega\,\mathbf{C}\right)^{-1}\mathbf{B}

Choosing the Load Basis $\mathbf{B}$

$\mathbf{B}$ is a modeling decision, not a formality. Nodal point loads on every candidate contact node give maximal flexibility and maximal ill-posedness. Piecewise-constant pressure patches, or a low-order pressure basis over the candidate surface, shrink $n$ and build smoothness into the parametrization itself — a structural prior that does part of the regularization work before any $\lambda$ appears. A related but distinct move appears in Section 6: Li et al. (2025) compress the reconstructed displacement field with radial basis functions and then differentiate it, rather than compressing the load basis $\mathbf{B}$.

The Fully-Observed Special Case: $\mathbf{f} = \mathbf{K}\mathbf{u}$

Tactile sensing enables a degenerate case worth singling out. If all retained DOFs are measured — Ma et al. (2019) mesh the gel as a single layer of elements with a fixed base, so every retained node receives a completed displacement estimate — the tangential components recovered from the tracked image-plane motion once the projection of the normal motion is subtracted, the normal component itself supplied by the known indenter geometry — then $\mathbf{S} = \mathbf{I}$, no elimination is needed, and force recovery is a direct product:

\mathbf{f} = \mathbf{K}\,\mathbf{u}_{\mathrm{meas}}

No inversion at all: multiplication. But the pathology does not vanish — it relocates. $\mathbf{K}$ acts like a discrete second-derivative operator, so pixel-scale tracking noise in $\mathbf{u}_{\mathrm{meas}}$ is multiplied by the largest stiffness eigenvalues. Ill-posedness enters wherever the observation is partial, noisy, or the model is wrong — never only in the matrix inverse. With perfect dense data, $\mathbf{f}=\mathbf{K}\mathbf{u}$ is exact; with real data, it is a noise amplifier that still needs regularization of the kind Section 4 develops — which is the door iFEM2.0's ridge term walked through (Section 8).

Two bookkeeping traps that corrupt force maps

1. Constrained-DOF rows return reactions. Evaluating $\mathbf{f} = \mathbf{K}_{\mathrm{full}}\mathbf{u}$ on rows belonging to DOFs with imposed displacements (the bonded base, clamped edges) — rows the reduced $\mathbf{K}_{ff}$ does not even carry — yields support reactions, not applied contact loads. If your candidate contact surface touches a constrained region, applied load and reaction are inseparable there — exclude constrained DOFs from the reported force map, or report them separately as reactions.

2. Consistent nodal forces are not pointwise tractions. A nodal force is the shape-function-weighted integral of traction over the adjacent elements, $f_j = \int_{\Gamma} N_j\, t\, \mathrm{d}\Gamma$ — its magnitude scales with element size. Refine the mesh and every nodal force shrinks while the underlying traction field is unchanged. To report pressure, divide by the tributary area (lumped) or solve the boundary mass-matrix system. Plotting raw nodal forces as if they were a pressure map is a units error that mesh refinement silently exposes.

3. Why Naive Inversion Explodes: SVD Applied to $\mathbf{H}$

Module 3 built the full diagnostic machinery — SVD anatomy, the discrete Picard condition, filter factors (see Module 3, Sections 5–8). This section does not re-derive any of it; it applies it to the transfer matrix and extracts the two numbers that matter for force sensing: how fast the singular values fall, and what that decay costs you in spatial resolution.

Write $\mathbf{H} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\mathsf T} = \sum_i \sigma_i \mathbf{u}_i \mathbf{v}_i^{\mathsf T}$ with $\sigma_1 \ge \dots \ge \sigma_r > 0$, where $r = \operatorname{rank}(\mathbf{H}) \le \min(m,n)$. Only those $r$ triplets carry information; the minimum-norm least-squares solution expands as

\hat{\mathbf{g}}_{\mathrm{LS}} = \sum_{i=1}^{r}\frac{\mathbf{u}_i^{\mathsf T}\mathbf{y}}{\sigma_i}\,\mathbf{v}_i,\qquad \kappa(\mathbf{H}) = \frac{\sigma_1}{\sigma_r}

Sum through $r$, not $n$: the page's own widgets run cases with $m < n$ and with numerical rank loss, where there simply are no $n$ nonzero singular values and the $i > r$ terms would divide by zero. Read $\sigma_1/\sigma_r$ as the conditioning of the observed subspace — it is not the conditioning of load recovery, which when $r < n$ is unbounded. In that case the load is not identifiable: any component in $\operatorname{span}\{\mathbf{v}_{r+1},\dots,\mathbf{v}_n\}$ leaves no trace in $\mathbf{y}$, so the formula above returns the minimum-norm representative of a whole affine family rather than the truth.

Because elasticity smooths, the high-index right singular vectors $\mathbf{v}_i$ are spatially oscillatory force patterns whose displacement signatures $\mathbf{u}_i$ carry tiny $\sigma_i$: the singular values of discretized first-kind problems decay fast — near-exponentially for smooth or depth-separated kernels, algebraically in the high-frequency limit for the weakly singular same-surface kernel of Sections 7–9 (Hansen 1998). Additive measurement noise $\mathbf{e}$ contributes $\mathbf{u}_i^{\mathsf T}\mathbf{e}/\sigma_i$ to mode $i$. With $\sigma_r/\sigma_1 = 10^{-6}$ and 1% relative noise, the last modes contribute errors up to $10^4$ times the signal — the reconstruction is a large-amplitude oscillation bearing no resemblance to the true load. Discretized load-identification problems routinely reach $\kappa(\mathbf{H})$ of $10^6$–$10^{12}$.

Key insight — $\kappa(\mathbf{K})$ and $\kappa(\mathbf{H})$ are different diseases

The stiffness matrix's conditioning grows polynomially with mesh refinement ($\sim h^{-2}$) — a discretization artifact, treatable by preconditioning, and it caps how many digits a forward solve loses. The transfer matrix's singular-value decay is a property of the physics, not of the discretization: the elastic kernel smooths, so oscillatory loads are faint at the sensors no matter how you mesh — near-exponentially faint for smooth or depth-separated kernels, algebraically so in the high-frequency limit for the same-surface kernel of Sections 7–9. Refining the candidate-load mesh makes $\mathbf{H}$ worse, not better — you add ever-fainter columns and ever-smaller singular values. No preconditioner fixes information that never reached the sensors.

The Discrete Picard Condition on Real Data

The regularized solution is meaningful only while the data coefficients $|\mathbf{u}_i^{\mathsf T}\mathbf{y}|$ decay faster than the $\sigma_i$ (Hansen 1998). On a Picard plot of both series against mode index, noisy data coefficients track the singular values initially, then flatten at the noise floor while $\sigma_i$ keeps falling. The index where they flatten — equivalently, where the flat floor meets the still-falling $\sigma_i$ — separates recoverable modes from noise-dominated ones. Count it before promising a spatial resolution. Read the flattening, not the literal crossing point: the two series carry different units, so their vertical offset depends on how the load is scaled, and only the shape of the coefficient curve is scale-free.

Anatomy of a Picard plot for a transfer matrix. The singular values $\sigma_i$ (physics) fall steadily; the data coefficients $|\mathbf{u}_i^{\mathsf T}\mathbf{y}|$ (measurements) flatten at the marker-noise floor. Left of $i^*$, where the floor meets the $\sigma_i$ curve: recoverable traction modes, the ones any filter must keep. Right: noise-dominated modes that any sane filter must suppress. No $\lambda$ is chosen here, so no filter factor $\varphi_i$ is defined — $i^*$ is a property of the data, and picking $\lambda$ near $\sigma_{i^*}$ is what turns it into one.

The Resolution Economics of Force Sensing

How much resolution do the surviving modes buy? A model problem gives the scaling — treat this as a design estimate from standard deconvolution analysis, not a theorem about your particular sensor. For a Gaussian smoothing kernel of width $w$, the singular values behave like $\mathrm{e}^{-k^2 w^2/2}$ in spatial frequency $k$, so modes survive above noise-to-signal ratio $\eta$ only for

k_{\max} = \frac{\sqrt{2\ln(1/\eta)}}{w}

Numbers: what an SNR upgrade actually buys

Set $\eta = 1/\mathrm{SNR}$. At SNR $= 100$: $k_{\max} = \sqrt{2\ln 100}/w \approx 3.0/w$. At SNR $= 10{,}000$ — a hundred times better measurement — $k_{\max} \approx 4.3/w$: barely 40% more resolution. The square root of a logarithm is one of the cruelest functions in engineering — for this kernel. Read it as the pessimistic bracket, and note how narrow the bracket is: the square root comes from the squared exponent of the Gaussian. In general, for $\sigma(k) \sim \mathrm{e}^{-c(wk)^q}$ the threshold gives $k_{\max} \sim [\ln(\mathrm{SNR})/c]^{1/q}/w$, so $q = 2$ (Gaussian) gives $\sqrt{\ln}$ while ordinary exponential decay $\mathrm{e}^{-cwk}$ ($q = 1$) gives $k_{\max} \sim \ln(\mathrm{SNR})/(cw)$ — better, though still logarithmic. The same-surface kernel of Sections 7–9 decays only algebraically, which turns the SNR dependence into a power law. Meanwhile the kernel width $w$ sits in the denominator either way: halving it doubles $k_{\max}$ outright. Section 9 turns that into a design rule, because for an elastic layer $w$ is set by the pad thickness.

Key insight — for a Gaussian kernel, resolution grows with the square root of the logarithm of SNR

Better cameras, better marker tracking, more averaging: against a Gaussian kernel, all of it buys spatial resolution only through $\sqrt{\ln(\mathrm{SNR})}$. The square root is specific to the Gaussian's quadratic exponent, not a property of exponential decay in general: for $\sigma(k) \sim \mathrm{e}^{-c(wk)^q}$ the survival threshold gives $k_{\max} \sim [\ln(\mathrm{SNR})/c]^{1/q}/w$, so plain exponential decay ($q = 1$) already buys the full $\ln(\mathrm{SNR})$. And the exponent is a property of the kernel, not a universal law at all: when the singular values fall only algebraically, $\sigma(k) \sim k^{-p}$ — the regime of the same-surface elastic kernel of Sections 7–9 — the survival threshold $\sigma(k) \gtrsim 1/\mathrm{SNR}$ gives $k_{\max} \sim \mathrm{SNR}^{1/p}$ instead, a power law, and measurement upgrades pay off far better than the Gaussian estimate suggests. What does not change is the geometry term: thinner gel buys resolution linearly under either kernel, which is why Section 9 treats skin geometry — not sensing hardware — as the primary design knob for the inverse problem.

4. Tikhonov-Regularized Inversion

The repair (Tikhonov & Arsenin 1977) replaces the unstable least-squares problem with a penalized one:

\hat{\mathbf{g}}_{\lambda} = \arg\min_{\mathbf{g}}\;\lVert\mathbf{H}\mathbf{g}-\mathbf{y}\rVert_2^{2} + \lambda^{2}\,\lVert\mathbf{L}\mathbf{g}\rVert_2^{2}

with closed-form normal equations

\left(\mathbf{H}^{\mathsf T}\mathbf{H} + \lambda^{2}\mathbf{L}^{\mathsf T}\mathbf{L}\right)\hat{\mathbf{g}}_{\lambda} = \mathbf{H}^{\mathsf T}\mathbf{y}

Standard form $\mathbf{L} = \mathbf{I}$ penalizes force magnitude. The general form, with $\mathbf{L}$ a discrete gradient or surface Laplacian, penalizes roughness of the traction field — usually the physically right prior for contact pressure, which is smooth over the patch — and is handled via the GSVD or a standard-form transformation. In the singular basis (for $\mathbf{L} = \mathbf{I}$), Module 3's filter-factor result applies verbatim to $\mathbf{H}$:

\hat{\mathbf{g}}_{\lambda} = \sum_{i=1}^{r}\varphi_i\,\frac{\mathbf{u}_i^{\mathsf T}\mathbf{y}}{\sigma_i}\,\mathbf{v}_i = \sum_{i=1}^{r}\frac{\sigma_i}{\sigma_i^{2}+\lambda^{2}}\,\left(\mathbf{u}_i^{\mathsf T}\mathbf{y}\right)\mathbf{v}_i,\qquad \varphi_i = \frac{\sigma_i^{2}}{\sigma_i^{2}+\lambda^{2}}\;\;(\mathbf{L}=\mathbf{I})

The second form is the one to implement: it sums only over the $r$ nonzero triplets and never writes $\varphi_i/\sigma_i$ as the $0/0$ it becomes in a null direction. Read $\sigma_i/(\sigma_i^2+\lambda^2)$ there as the whole coefficient $\varphi_i/\sigma_i$, not as the filter factor — the very confusion the next paragraph warns about. And note what regularization does not do: it never resurrects $\ker(\mathbf{H})$; with $\mathbf{L} = \mathbf{I}$ the penalty simply sets those components to zero.

— a smooth roll-off passing modes with $\sigma_i \gg \lambda$ and damping modes with $\sigma_i \ll \lambda$ by the factor $(\sigma_i/\lambda)^2$. Mind the exponents: the filter factor is $\sigma_i^2/(\sigma_i^2 + \lambda^2)$, not $\sigma_i/(\sigma_i+\lambda)$ and not $\sigma_i/(\sigma_i^2+\lambda^2)$ — that last combination is the full coefficient $\varphi_i/\sigma_i$ from the display above, and promoting it to $\varphi_i$ changes the roll-off. The quadratic roll-off is what matches Tikhonov to a TSVD truncation level. Truncated SVD is precisely the hard-cutoff cousin: $\varphi_i = 1$ for $i \le k$, else $0$; with $\lambda$ matched to $\sigma_k$ the two give similar reconstructions.

What the Bias Looks Like on a Contact Patch

Small $\lambda$ leaves noise amplification (variance); large $\lambda$ over-smooths (bias). For contact problems the bias has a recognizable signature: the reconstructed patch is spread wider than the truth and the peak pressure is underestimated, while the total force — a smooth functional that lives in the well-conditioned low modes — is far less affected — less affected, not conserved: plain Tikhonov gives no exact guarantee on the resultant unless you impose one. This asymmetry is why resultant-force checks pass long after the pressure map has degraded (Section 8), and you can watch it happen live in the lab below.

Key insight — regularization declines to invent, it does not recover

No choice of $\lambda$ "restores" the information that elasticity smoothed away; regularization only trades variance for bias among the modes the data actually constrains. The resolution limit of Section 3 stands regardless — square-root-logarithmic in SNR for the Gaussian model kernel, a power law $\mathrm{SNR}^{1/p}$ when the singular values fall only algebraically, but in both cases fixed by the data and untouched by $\lambda$. Any pipeline that claims regularization "recovered" sub-kernel-width force detail is describing its prior, not its measurement.

Priors, Statistics, and Physical Constraints

Statistically, Tikhonov is MAP estimation under Gaussian noise with the Gaussian prior $\mathbf{g} \sim \mathcal{N}\!\left(\mathbf{0},\, \sigma_e^2 \lambda^{-2} (\mathbf{L}^{\mathsf T}\mathbf{L})^{-1}\right)$ — the bridge to the Bayesian and Kalman branches of load identification surveyed by Liu et al. (2022); Module 3 covers the derivation. That covariance only exists when $\mathbf{L}$ has full column rank. The standard form $\mathbf{L} = \mathbf{I}$ does; a discrete gradient or surface Laplacian does not — those annihilate constants, so $\mathbf{L}^{\mathsf T}\mathbf{L}$ is singular and $(\mathbf{L}^{\mathsf T}\mathbf{L})^{-1}$ does not exist. Read the general form instead as an improper (intrinsic) Gaussian specified by its precision $\lambda^{2}\sigma_e^{-2}\,\mathbf{L}^{\mathsf T}\mathbf{L}$: flat along $\ker(\mathbf{L})$, so the mean-pressure component is left entirely to the data. The penalized problem then still has a unique minimizer, but only under the standard general-form condition $\ker(\mathbf{H}) \cap \ker(\mathbf{L}) = \{\mathbf{0}\}$.

Physical side constraints do regularization work too, and deserve first-class status rather than afterthought status:

Unilateral contact: normal nodal forces satisfy $g_i \ge 0$ (a surface can push, not pull, absent adhesion) — solved by NNLS / active-set methods, and often worth more than a decade of $\lambda$ tuning because it forbids the oscillating negative lobes that noise amplification produces;
Friction cone: tangential traction bounded by $\lVert \mathbf{t}_t \rVert \le \mu\, p$ — couples the tangential reconstruction to the normal one;
Equilibrium / total-force equality: when a force/torque sensor is in the loop, constrain the reconstructed resultant to match its reading.

Ridge regression — standard-form Tikhonov by another name — is what iFEM2.0 (Zhao, Liu & Ma, T-RO 2025) adds to the tactile pipeline of Section 8 for noise robustness, alongside multilayer mesh constraints. Read the paper before assuming the filter factors above transfer verbatim: with a multilayer mesh the surface displacements are measured and the interior ones are not, so which unknown the ridge term penalizes — the tractions or the unobserved interior displacements — decides whether Section 4's bias signature applies as written or only by analogy.

Interactive Tool — Contact Inversion Lab

Press a 1-D compliant surface and reconstruct the contact pressure from marker displacements. The elastic layer is modeled by a Lorentzian smoothing kernel whose width scales with pad thickness $t$ — an explicit surrogate for the thin-layer Green's function, not the exact kernel. Choose the load profile, thin out or occlude the markers, add noise, then tune $\lambda$ and the physical constraints. Watch the panel: measurements, unknowns, assumptions, conditioning, and the two validation numbers that this module keeps insisting are different — resultant error vs. pressure-map error.

Load profile (truth)

Pad thickness $t$ 0.10

Bump separation $s$ 0.25

Noise $\eta$ (% of max $|y|$) 1.0

$\log_{10}\lambda$ −3.0

Method

Physical constraints

Marker spacing

markers vanish under the contact

Try this: (1) Set $\eta = 0$, marker spacing = every node, method = naive LS: inversion is essentially exact on this finite synthetic problem. Do not over-read that: ill-posedness is failure of continuous dependence, a property of the operator that is present before any noise is added. Exact data and exact arithmetic merely let you reproduce one particular truth you generated yourself; roundoff, model error, or a finer discretization expose the instability that was there all along. (With sparse markers the problem is underdetermined and even noiseless data only pins down the minimum-norm solution — and at the thick end of the $t$ slider the kernel goes numerically rank-deficient even with every node observed, so the same caveat applies there: watch the rank and nullity readouts.) (2) Restore $\eta = 1\%$, keep naive LS, and watch the clipped blow-up; then switch to Tikhonov and slide $\lambda$. (3) Two-fingertip preset: raise $t$ and watch the two bumps blur toward one broad mound, with the $\lambda$ range that still separates them shrinking as you go — the resolution crossover of Section 3, live. It is a crossover and not a floor, and the widget will show you that too: with $\eta = 0$, every node observed and $\lambda = 10^{-6}$, even $t = 0.30$ still returns two clean peaks (about 84% of the true height). Separation never becomes impossible on this exact-arithmetic synthetic problem; it becomes unaffordable, because the $\lambda$ you would need lets the noise through. (4) Occlude the markers under the contact: the load right where it is applied is now unobserved, and the resultant error jumps by about an order of magnitude — while the displayed $\sigma_1/\sigma_r$ falls. Deleting sensors deletes modes rather than making the surviving ones harder to see, so a shrinking condition number here is not good news. Watch the rank and the nullity instead: they are what actually got worse. (5) At moderate $\lambda$ (say $10^{-3}$ to $10^{-1}$), compare the two validation numbers: the resultant drifts far less than the map, which spreads while the peak drops — Section 8's partial-validation trap in miniature. Push $\lambda$ to the top of the slider and the ordering flips (resultant error $\approx 70\%$ against a map RMSE of $\approx 30\%$): that is over-smoothing past the point where any functional survives. The two percentages carry different denominators — true resultant versus true peak — so read their trends, not their difference. The friction-cone constraint is the tangential analogue of the $p \ge 0$ option; it needs a 2-D vector traction field, so here you get its scalar cousins.

5. Choosing $\lambda$: L-Curve, GCV, Discrepancy

Module 3 (Section 9) derives the three standard parameter-choice rules; here is what each one needs, and — more usefully — how each one fails on load-identification data. None of them is ground truth. Each is an estimator with a documented failure mode, and the failure modes are different, which is exactly why you run more than one.

Criterion	Chooses $\lambda$ by	Needs	Fails when
L-curve (Hansen 1992)	Corner (max curvature) of $\log\lVert\mathbf{L}\hat{\mathbf{g}}_\lambda\rVert$ vs. $\log\lVert\mathbf{H}\hat{\mathbf{g}}_\lambda-\mathbf{y}\rVert$	Nothing beyond the sweep	Very smooth exact solutions (corner disappears); extremely fast $\sigma_i$ decay
GCV (Golub, Heath & Wahba 1979)	Minimum of $G(\lambda)$ — rotation-invariant leave-one-out	No noise-level estimate	Correlated noise: minimum flattens or drifts toward undersmoothing
Discrepancy (Morozov)	Residual $= \tau\delta$, $\tau$ slightly above 1	The noise level $\delta$	$\delta$ misjudged: overestimate $\Rightarrow$ oversmoothing

G(\lambda) = \frac{m\,\lVert(\mathbf{I}-\mathbf{A}_{\lambda})\mathbf{y}\rVert_2^{2}}{\left[\operatorname{tr}\!\left(\mathbf{I}-\mathbf{A}_{\lambda}\right)\right]^{2}},\qquad \mathbf{A}_{\lambda} = \mathbf{H}\left(\mathbf{H}^{\mathsf T}\mathbf{H}+\lambda^{2}\mathbf{L}^{\mathsf T}\mathbf{L}\right)^{-1}\mathbf{H}^{\mathsf T}

The L-curve corner is operationalized as the point of maximum curvature,

\kappa(\lambda) = \frac{\rho'\eta'' - \rho''\eta'}{\left(\rho'^{2}+\eta'^{2}\right)^{3/2}},\qquad \rho = \log\lVert\mathbf{H}\hat{\mathbf{g}}_{\lambda}-\mathbf{y}\rVert_2,\;\; \eta = \log\lVert\mathbf{L}\hat{\mathbf{g}}_{\lambda}\rVert_2

and the discrepancy rule reads $\lVert\mathbf{H}\hat{\mathbf{g}}_{\lambda}-\mathbf{y}\rVert_2 = \tau\,\delta$. For tactile pipelines the GCV failure mode is not exotic: interpolating sparse marker displacements onto FE nodes correlates the noise across nodes, which is precisely the condition under which GCV drifts. Meanwhile $\delta$ is often genuinely measurable — record the unloaded sensor for a minute and read the marker-tracking repeatability — which promotes the discrepancy principle from textbook afterthought to the most reliable rule in the room.

Key insight — disagreement is a diagnostic

For $\mathbf{L} = \mathbf{I}$ the SVD is computed once ($n$ up to a few thousand is trivial), and a sweep of 100–200 log-spaced $\lambda$ values from $\sigma_r$ (the smallest nonzero singular value) to $\sigma_1$ evaluates all three criteria through filter factors at negligible cost. For a gradient or Laplacian $\mathbf{L}$ the same trick needs the GSVD of $(\mathbf{H},\mathbf{L})$ or a standard-form transformation first: the filter factors and the meaningful range of $\lambda$ are then set by the generalized singular values, which carry the relative scaling and units of both matrices, not by $\sigma_r \ldots \sigma_1$. Run all three. If the L-curve and GCV agree within an order of magnitude, proceed. If they disagree wildly, that disagreement is information: suspect correlated noise or model error before trusting either number.

Real-time practice: $\lambda$ is not chosen per frame. Freeze it offline during calibration, precompute the regularized pseudoinverse $(\mathbf{H}^{\mathsf T}\mathbf{H} + \lambda^2\mathbf{L}^{\mathsf T}\mathbf{L})^{-1}\mathbf{H}^{\mathsf T}$, and per-frame reconstruction becomes a single $n \times m$ matrix–vector product. That one precomputation is the entire gap between "inverse problem" and "sensor". Two caveats the formula hides. First, it is written that way for the algebra only: build it from the SVD/GSVD of $\mathbf{H}$ or an augmented least-squares QR, never by forming and inverting $\mathbf{H}^{\mathsf T}\mathbf{H}$ — the $\lambda^2$ term bounds the conditioning of the penalized system (for $\mathbf{L}=\mathbf{I}$, $\kappa_2 = (\sigma_1^2+\lambda^2)/(\sigma_{\min}^2+\lambda^2) \le 1 + \sigma_1^2/\lambda^2$, with $\sigma_{\min} = 0$ in a null direction — still enormous at small $\lambda$), and it cannot restore digits already lost in squaring $\mathbf{H}$ (Section 2). Second, the matrix–vector product is the unconstrained solution; impose NNLS or a friction cone and each frame becomes a small warm-started QP or SOCP over the active set instead — still real-time at these sizes, but no longer one product.

Interactive Tool — L-Curve & GCV Explorer

The fixed two-bump problem from the lab above — same 48-node discretization, $t = 0.1$, $s = 0.25$, $\eta = 1\%$ — fully observed. For that matrix $\sigma_1 \approx 0.81$ and $\sigma_{48} \approx 6.0\times10^{-7}$, so $\kappa \approx 1.4\times10^{6}$. All curves come from one precomputed SVD swept over 150 log-spaced $\lambda$ values from $10^{-6}$ to $1$ — roughly $\sigma_{\min}$ to $\sigma_1$, the range Section 5 recommends. Drag the $\lambda$ slider or click on the L-curve; compare where the corner, the GCV minimum, and the discrepancy rule land against the oracle (the truth-knowing best $\lambda$ — a luxury of simulation). Then press the failure-mode buttons: redraw the noise (the corner and GCV minimum jitter between draws) and switch to correlated noise (the GCV minimum usually drifts toward undersmoothing while the discrepancy choice stays sane).

$\log_{10}\lambda$ −3.0

correlated noise (moving-average width 5)

overlay oracle RMSE($\lambda$)

L-curve (log–log). Click to set $\lambda$.

GCV function and oracle RMSE — each curve normalized to its own range; compare the minima locations, not heights.

Reconstruction at the current $\lambda$ vs. true traction. The vertical axis is fixed to the truth's range, so an undersmoothed curve runs off-panel and is pinned to the frame — the panel says so and prints the real $\max\lvert\hat{g}\rvert$ whenever that happens.

6. Time-Domain Force Identification

Add inertia and the static picture becomes a movie:

\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{B}\,\mathbf{g}(t)

The measured response is the Duhamel convolution of the input with the impulse response function (IRF):

y(t) = \int_{0}^{t} h(t-\tau)\,g(\tau)\,\mathrm{d}\tau

Sampling with step $\Delta t$ turns this into a linear system with a lower-triangular Toeplitz operator built from IRF samples $h_k = h(k\,\Delta t)$:

y_i = \Delta t\sum_{j=1}^{i} h_{\,i-j+1}\,g_j \;\;\Longleftrightarrow\;\; \mathbf{y} = \mathbf{A}\,\mathbf{g}

Mind the index offset: with $h_k = h(k\Delta t)$ this sum places $g_j$ at $t = (j-1)\Delta t$, so the diagonal entry of $\mathbf{A}$ is $h_1 = h(\Delta t)$. That is deliberate — a displacement IRF satisfies $h(0) = 0$, so indexing the sum with $h_{i-j}$ instead would put zeros on the diagonal and make $\mathbf{A}$ singular rather than merely ill-conditioned.

— a deconvolution problem, severely ill-posed for the same structural reason as the static case, now in time: the IRF of a damped structure is smooth, so $\mathbf{A}$'s singular values decay rapidly. Light damping makes it worse — long memory, and errors that accumulate as the record grows. The canonical analysis is Jacquelin, Bennani & Hamelin (JSV 2003): impact-force reconstruction posed as regularized deconvolution,

\hat{\mathbf{g}} = \arg\min_{\mathbf{g}}\;\lVert\mathbf{A}\mathbf{g}-\mathbf{y}\rVert_2^{2} + \lambda^{2}\lVert\mathbf{g}\rVert_2^{2}\quad\left(\text{or } +\,\lambda\lVert\mathbf{g}\rVert_1\right)

comparing Tikhonov with L-curve and GCV parameter selection. Their conclusions transfer wholesale: regularization is mandatory, and the parameter choice — not the factorization — is the delicate step.

Frequency Domain, and Where It Blows Up

The alternative is per-frequency inversion, $\mathbf{G}(\omega) = \mathbf{H}(\omega)^{+}\,\mathbf{Y}(\omega)$, using the FRF matrix of Section 2. Its Achilles heel is spectral: ill-conditioning concentrates near the antiresonances of the FRF, where $|\mathbf{H}(\omega)|$ is small and per-line inversion amplifies noise violently. Transient signals add leakage and windowing artifacts on top. Hence the standard division of labor: time domain for impacts, frequency domain for stationary loading — presenting frequency-domain inversion without its antiresonance blowup is selling half a method.

The Modern Method Map

Five families now coexist. The grouping below is this guide's own synthesis, drawing on the review of Liu, Dobriban, Hou & Qian (2022) but not reproducing its organization — the review runs along frequency- and time-domain formulations, ill-posedness treatment, uncertain-parameter and statistical methods, and AI-based methods, which cuts across these rows rather than matching them (sparse regularization, for instance, is an ill-posedness remedy there, not a parallel family):

Family	Core move	Character
(i) Frequency-domain FRF inversion	Per-line pseudoinverse of $\mathbf{H}(\omega)$	Batch; stationary loads; antiresonance-fragile
(ii) Time-domain deconvolution	Tikhonov/TSVD on the Toeplitz system	Batch; impacts; the Jacquelin et al. (2003) canon
(iii) State-space recursive estimators	Augment the state with the unknown input modeled as a random walk; run a Kalman-type joint input–state filter (Lourens et al. 2012)	Online, with uncertainty; the augmented Kalman filter is the workhorse
(iv) Sparse regularization	Replace $\ell_2$ with $\ell_1$: impacts are sparse in time (or in a wavelet dictionary)	Recovers sharp impacts that Tikhonov smears
(v) Learning-based mappings	Train response $\to$ load directly (e.g., deep CNN load identification, Yang et al., MSSP 2023)	Fast inference; inherits training-distribution limits

Key insight — the prior picks the family

Tikhonov's $\ell_2$ penalty encodes "the force history is small in energy" — magnitude only, not smoothness (that needs a derivative penalty) — — the worst possible prior for a hammer blow. The $\ell_1$ family encodes "the force is zero almost always" — exactly right for impacts, useless for slowly varying service loads. The joint input–state filters encode "the force is a random walk" — right for online tracking. Choosing the method is choosing the prior; make it consciously.

Unknown Locations and Distributed Loads

When the load location is unknown, solve on a grid of candidate positions and select by minimum residual — a discrete search wrapped around the linear solver. When the load is distributed in space and varying in time, the unknown count multiplies: with sparse measurements you usually need dimension reduction, spatial–temporal regularization, or both. Whether the problem is underdetermined is a question about the number and rank of the measurements against the number of coefficients, not a law — dense enough data can support a full spatial grid, and low-order compression is one option among several. Li et al. (2025) take the compression route on a different quantity: they reconstruct distributed dynamic loads on plates in a chain that never forms a transfer matrix at all — Tessler-style iFEM turns measured strains into a discrete displacement field, Gaussian radial basis functions fit that field into a differentiable continuous one, and the thin-plate equation of motion then returns the distributed load. The RBFs compress and differentiate the displacement, not the load. That avoids forming and inverting a transfer matrix, which is the sense in which the authors describe their method as preventing ill-posedness; it does not make the recovery insensitive to noise, since the plate equation reaches the load through fourth spatial derivatives of the displacement field. The iFEM least-squares fit and the RBF fit are what smooth the field before it is differentiated — implicit regularization, not its absence (Section 1's warning box).

7. From Structures to Skin: the Tactile Inverse Problem

Now recast the problem for a vision-based tactile sensor: a millimeters-thick elastomer pad (silicone, $E$ of order 0.1–1 MPa) bonded to a rigid window, camera behind it. The measurement channels are the tracked displacements of a printed marker array (the GelSight/GelSlim lineage) and/or photometric-stereo depth. The continuous forward model is a boundary-integral relation between surface traction and surface displacement:

\mathbf{u}(\mathbf{x}) = \int_{\Gamma_{c}} \mathbf{G}(\mathbf{x}-\mathbf{x}')\,\mathbf{t}(\mathbf{x}')\,\mathrm{d}\Gamma'

Key insight — the data economics flip

Structural SHM measures a handful of strain gauges far from unknown load locations; tactile sensing measures thousands of displacement vectors on the loaded boundary itself. The kernel relating surface traction to surface displacement is far less smoothed at zero depth than at gauge distance — the singular values decay much more slowly, the conditioning is dramatically better, and this is why dense contact force reconstruction is feasible at all. Where you measure is the single biggest conditioning lever in the whole problem (Section 9).

What Still Makes It Hard

2-D projections of 3-D motion. The camera sees image-plane displacements, which conflate tangential motion with normal indentation through perspective and refraction. Ma et al. (2019) derive the correction: a marker displaced by $d$ along the normal appears shifted in-plane by $\beta = \arcsin\!\left(\frac{\sin\alpha}{\gamma}\right),\qquad \boldsymbol{\Delta} = \frac{\mathbf{r}}{\lvert\mathbf{r}\rvert}\;d\,\tan\beta$ with $\alpha$ the viewing-ray incidence angle, $\gamma$ the silicone refractive index, and $\mathbf{r}$ the radial image position (acrylic refraction neglected by index-matching design). Feed uncorrected 2-D displacements to the inverse solver and the tangential force field is biased.
Out-of-plane displacement often is not directly measured — it must come from photometric depth or from assumed contact geometry.
Soft silicones are nonlinear (hyperelastic, some viscoelastic) at deep indentation, while the pipeline assumes linear elasticity.
Marker-interpolation noise is correlated — the quiet GCV-killer of Section 5.

The GelSlim 2.0 sensing stack, drawn from the description in Ma et al. (2019): markers printed on the outer gel surface under the paint, dual-color grazing illumination, curved total-internal-reflection light guide, angled camera with its mirrored virtual position. Inset: a marker displaced $d$ along the normal appears shifted in-plane by $\Delta = d\tan\beta$ with $\beta = \arcsin(\sin\alpha/\gamma)$ — a systematic error subtracted before inversion. The ray bends toward the normal on entering the denser medium ($\beta < \alpha$, drawn here at $\sin\alpha/\sin\beta \approx 1.4$); since gel and acrylic are index-matched by design, that single bend stands in for the refraction at the camera-side boundary rather than at the gel/window interface.

Three Solution Families

Family	Idea	Strengths	Weaknesses
Analytical (GelForce-style)	Two-layer marker sensors build a dense influence matrix from Boussinesq–Cerruti half-space Green's functions and invert it globally: Kamiyama et al. assemble a constant $3m \times 3n$ block matrix $\mathbf{H}$ with $\mathbf{H}_{i,j} = \mathbf{H}(\mathbf{x}_j - \boldsymbol{\xi}_i)$, precompute the pseudoinverse $\mathbf{Q} = (\mathbf{H}^{\mathsf T}\mathbf{H})^{-1}\mathbf{H}^{\mathsf T}$ offline, and apply it per frame as one matrix–vector product	Closed-form, fast, no mesh	Semi-infinite assumption is invalid for a thin bonded pad — Ma et al. (2019) criticize precisely this, and read the semi-infinite assumption as forcing a per-marker relation — in their words, “a force vector acting on a single marker is linearly dependent only on the displacement vector of that particular marker, which is not valid in general cases.” That reading is not what GelForce does: its inversion is global, and locality does not follow from semi-infiniteness in any case, because the half-space Green’s function is itself nonlocal — a surface displacement superposes tractions from across the whole surface. The half-space assumption is the real weakness, and it fails hardest for a thin bonded pad
Model-based FE	Geometry-true stiffness model; this module's case study	Physics-consistent; calibration-light (two material constants); interpretable	Linear-elastic assumption strains at deep indentation; needs a good mesh and BCs
Learning-based	Sferrazza et al. (2019) invert the pipeline's role of FEM: offline FE simulations of the gel generate dense ground-truth force-distribution labels (resultants validated against a commercial F/T sensor); a network trained on thousands of image–label pairs reconstructs the 3-D force distribution in real time. FeelAnyForce (2024) scales the learned route across sensor types with a multi-head transformer trained on 200k+ indentations, but regresses the resultant 3-D contact force rather than a traction field — Section 8's resultant/distribution distinction again	Handles nonlinearity; simple runtime pipeline	Needs the label generator (itself an FE model) or massive data; training-distribution limits

The tradeoff in one line: model-based buys physics consistency and calibration-lightness; learned buys nonlinearity handling and pipeline simplicity. The most interesting current work (Section 8's successor discussion) imports the regularization machinery of Sections 4–5 into the model-based route rather than abandoning it.

8. Case Study: GelSlim 2.0 and Inverse FEM (Ma et al., ICRA 2019)

One paper, walked end to end — because it is the cleanest existing demonstration that a soft skin plus a calibrated FE model is a distributed force sensor. Keep the framing honest from the start: this is one bounded, material-dependent case — a specific gel, a specific geometry, linear elasticity, simple indenters — not proof of general dense 3-D force recovery. Its own authors call the validation partial; its successor exists largely to fix that.

Hardware

Markers as FE nodes: markers are printed as a regular grid directly correlated with the FE nodes, on the outer gel surface under the paint — the mesh is drawn on the sensor.
Dual-color grazing illumination: red LED from the left, green from the right (LUXEON CZ, chosen to match the Raspberry Pi Spy Camera's spectral detection bands), so contact geometry is encoded in RGB-channel differences.
Curved total-internal-reflection light guide replaces a 45° mirror — brighter image, higher frame rate, slimmer finger.
Parametric CAD scales the gel pad down to 20 × 20 mm and the finger thickness to 18 mm.

FE Model and the Key Formulation Move

The gel pad is meshed with a single layer of Hex-8 (8-node hexahedral, 24-DOF) elements — the paper reports 1964 elements and 3698 nodes — built slightly larger than the gel bounding box and then cropped. Those two counts do not quite fit together: one layer of Hex-8 elements has exactly two node planes, so 3698 nodes means 1849 per plane, while any conforming quad mesh of $Q$ faces with $B$ boundary edges has $V = Q + B/2 + 1 > Q$ vertices. A single conforming layer over 1964 elements therefore needs more than 3928 nodes. Read the pair as approximate, not as a mesh you can rebuild. Linear elasticity throughout:

\mathbf{F} = \mathbf{K}\,\mathbf{U},\qquad \mathbf{K}\in\mathbb{R}^{3n\times 3n}

\mathbf{U} = \left(\delta_{x}^{1},\delta_{y}^{1},\delta_{z}^{1},\ldots,\delta_{z}^{n}\right)^{\mathsf T},\qquad \mathbf{F} = \left(f_{x}^{1},f_{y}^{1},f_{z}^{1},\ldots,f_{z}^{n}\right)^{\mathsf T}

The only material inputs are $E$ and $\nu$, measured by tensile tests at several speeds — the authors report that the viscosity effect can be neglected, which licenses dropping viscoelasticity, not nonlinearity; linearity is a separate assumption, justified only over the small strains the sensor operates at. The paper reports $E = 147$ MPa and $\nu = 0.3223$, both flagged by the authors as sensitive to the gel preparation process. On physical grounds that modulus reads as a typo for $147$ kPa $= 0.147$ MPa (see the aside) — a hypothesis about the source, not a correction its authors state. Nothing on this page depends on the value either way, since $E$ enters as a pure gain (Section 9).

Aside: 147 MPa or 147 kPa?

The paper prints “147MPa” — plainly, with no dropped character; I checked the PDF rather than inferring a typographic accident. But $147$ MPa cannot be right for a cast silicone gel: a cast silicone sits around $0.1$–$1$ MPa, so $147$ MPa is two to three orders of magnitude too stiff — the range of a rigid thermoplastic rather than a gel (for scale, high-density polyethylene is stiffer still, around $1$ GPa) — and irreconcilable with the newton-scale forces the sensor resolves over millimetre indentations. $0.147$ MPa $= 147$ kPa is the physically consistent reading. A useful reminder that source-reading is part of modeling: a copied constant with a dropped decimal point silently rescales every reconstructed force by 1000.

The formulation move that makes it real-time: single layer of nodes + fixed bottom boundary (gel bonded to the acrylic) means every retained node has a displacement estimate — though "observed" deserves a footnote. The markers give image-plane motion, which mixes true tangential motion with the apparent shift that normal indentation produces under an oblique, refracting view (Section 7). The known indenter geometry does double duty: it supplies the normal ($z$) component, and its projection $\Delta = d\tan\beta$ is subtracted from the image-plane motion to leave the tangential ($x$, $y$) components. Ma et al. put it plainly — the observed displacements are "caused by both tangential and normal deformation of the gel pad," so the two "need to decouple." Full observation is therefore conditional on that geometric assumption, not a general property of the sensor. Granting it, this is Section 2's fully-observed special case: force recovery is the direct product $\mathbf{F} = \mathbf{K}\mathbf{U}$ — no ill-conditioned inversion online, $\mathbf{K}$ precomputed offline, estimation at camera frame rate. The pathology relocates rather than vanishes: $\mathbf{K}$ amplifies pixel-scale tracking noise through its largest eigenvalues, which is exactly the door through which iFEM2.0's ridge regularization later entered.

Displacement Pipeline

Step	What happens	Failure it guards against
1. Detect	Markers found per frame	—
2. Match	Nearest-neighbor matching with a validity check	Markers can vanish under deformation; naive matching then corrupts the field
3. Interpolate	Marker displacements interpolated onto FE nodes	— (but this step correlates the noise: Section 5)
4. Correct projection	Subtract $\Delta = d\tan\beta$, $\beta = \arcsin(\sin\alpha/\gamma)$ (Section 7)	Tangential/normal conflation through perspective and refraction
5. Reconstruct $z$	Out-of-plane displacements from contact-patch geometry for known simple indenters (sphere, cylinder)	The camera does not measure $\delta_z$ directly
6. Multiply	$\mathbf{F} = \mathbf{K}\mathbf{U}$	—

Results and Validation

The reconstructions deliver the punchline promised in Section 1: tangential forces concentrated on the contact patch even though displacements spread far outside it — the paper's own caption says the forces are “naturally limited in contact patch in most area, despite globally existing noises” — and near the patch boundary the tangential force direction can even oppose the local displacement — the global functional at work, impossible under any $f_i \propto u_i$ scheme.

Validation against an ATI Gamma F/T sensor while pushing and sliding on a sphere: standard deviation of the resultant force $(0.244,\ 0.201,\ 0.322)$ N in $(x, y, z)$ — roughly within 15% of ground truth, with the residual error attributed to calibration, marker tracking, and material nonlinearity.

Why it matters — "partial validation," in the authors' own framing

The F/T sensor checks the resultant: a handful of global linear functionals of the distribution — three force components as reported here, six once the torques of the six-axis ATI Gamma are used — all of them living in the best-conditioned modes. Infinitely many wrong pressure maps share a correct resultant — Section 4 showed regularization bias spreads the patch and cuts the peak while leaving this functional comparatively intact — a tendency, since resultants live in the well-conditioned low modes, not an exact conservation law. Resultant agreement is necessary, not sufficient; presenting it as validation of the force distribution is the field's most common overclaim. Ma et al. label their validation partial for precisely this reason. The fix requires distribution-level ground truth: Sferrazza et al. (2019) build it from a calibrated FE forward model; iFEM2.0 builds standardized simulation-based metrics.

The Successor: iFEM2.0

Zhao, Liu & Ma (IEEE T-RO 2025) upgrade the pipeline on both fronts this module has been flagging. Mechanically: a multilayer inverse FEM with mesh constraints plus ridge (Tikhonov) regularization for noise robustness — Sections 4–5 machinery brought into the tactile stack, subject to Section 4's caveat: a multilayer mesh leaves the interior displacements unobserved, so which unknown the ridge term penalizes decides how literally Section 4's bias signature transfers. Methodologically: a standardized evaluation framework with accuracy, fidelity, and noise-robustness metrics for dense 3-D contact force fields, demonstrated in simulation and experiment — a direct answer to the resultant-only validation gap.

9. Soft Robot Skins as Distributed Force Sensors: Design as Conditioning

The bridge to your own hardware: every skin design parameter is a knob on the inverse problem. Choose them with the SVD in mind and the regularization has less to repair; choose them blind and no $\lambda$ will save you.

Design knob	Effect on the inverse problem
Pad thickness $t$	The elastic kernel spreads a surface point load over a lateral scale of order $t$: thicker pads mean wider kernels, faster singular-value decay, and a characteristic spreading scale of order $t$ that sets where resolution starts to become expensive (kernel-width scaling estimate, Section 3 — a crossover, not a hard floor: what is actually recoverable also depends on the Green tensor, SNR, sampling and prior). But thicker also means larger displacement signal and more mechanical compliance for grasping — a resolution/robustness tradeoff, not a free win.
Modulus $E$	In linear elasticity with fixed geometry and $\nu$: $\mathbf{H} \propto 1/E$ — a pure gain. It does not change conditioning at all. It sets signal amplitude against the camera noise floor and the load range before nonlinearity.
Measurement density $m$	More markers average noise, but the textbook $1/m$ shrinkage per retained mode assumes independent, equal-variance, comparably informative measurements. In general a mode's noise variance is $\mathbf{u}_i^{\mathsf T}\boldsymbol{\Sigma}_e\mathbf{u}_i/\sigma_i^2$, so with the correlated marker-interpolation noise of §5 and §7 the benefit is smaller and eventually plateaus — count effective, not nominal, measurements. Either way density cannot defeat the singular-value decay — algebraic in the high-frequency limit for this bonded-layer geometry (§3) — so resolution improves with density, just slowly.
Where you measure	Sensing on the loaded surface (tactile) versus remote (structural SHM) is the single biggest conditioning lever — Section 7's data-economics flip.
Compressibility $\nu$	$\nu \to 0.5$ causes volumetric locking in standard displacement-based Hex-8 elements — the very element of Section 8's case study — requiring mixed/hybrid formulations. And note the measured GelSlim value $\nu = 0.3223$: real engineered gels can sit far from the incompressible idealization; assuming $0.5$ biases the model.

\mathbf{H} \propto \frac{1}{E}\quad(\text{fixed geometry},\ \nu),\qquad \ell_{c} \sim t

Read $\ell_c \sim t$ as the crossover scale, not a hard floor: below it, separating two features costs rapidly more SNR, denser sampling, or a stronger prior. Nothing forbids recovering detail finer than $t$ — the elastic kernel has no true cutoff — it simply stops being affordable.

Key insight — thickness is a resolution knob, modulus is a gain knob

Saying a softer or stiffer gel changes the ill-conditioning is wrong in the linear regime: $E$ scales every singular value equally and cancels out of $\kappa(\mathbf{H})$. Conditioning is controlled by geometry (thickness) and measurement placement. The corollary is diagnostic gold: a calibration error in $E$ rescales all reconstructed forces uniformly but leaves the distribution shape intact — provided the inversion is unregularized, or you rescale $\lambda$ with $E$ (nondimensionalize, and the shape really is invariant). Hold $\lambda$ fixed in absolute terms and the corollary weakens: scaling $\mathbf{H} \propto 1/E$ shifts every $\sigma_i/\lambda$, so the filter factors change mode by mode and the shape moves too. With that caveat: when your force map looks right but the magnitudes are off by a constant factor, recheck $E$; when the shape is wrong, suspect $\nu$, geometry, or boundary conditions first.

Ground truth remains the field's open sore. Point F/T sensors validate only resultants (Ma et al. 2019). Sferrazza et al. (2019) answer with a calibrated FE forward model as a dense label generator; iFEM2.0 (2025) answers with standardized simulation-based assessment metrics. Neither is a measured distributed reference; what does not yet exist is a minimally intrusive, calibrated, matched-resolution ground truth for the full 3-D vector traction field; distributed normal-pressure mats do exist and are used.

The takeaway recipe

Elastomer + camera + calibrated FE model + regularized inverse = a complete distributed force sensor — whose resolution is set by mechanics ($t$), whose accuracy is set by calibration ($E$, $\nu$, boundary conditions), and whose robustness is set by regularization.

10. Practical Recipe and Pitfalls

The Workflow

Characterize the elastomer at operating strain rates — tensile and/or indentation tests; check rate-independence as Ma et al. did to rule out viscoelasticity, and check stress–strain linearity over the operating strain range before assuming linear elasticity — rate tests alone license neither; go hyperelastic if operating strains exceed roughly 10%.
Verify the forward model first. Simulate a known indentation and compare force–depth against experiment. An inverse pipeline built on an unvalidated forward model is unfalsifiable — every wrong reconstruction can be blamed on the inversion. Hertz spherical contact is the standard sanity target — but it is a half-space result, valid only while the contact radius $a \ll t$ and the boundaries are remote. On a thin bonded pad a contact radius comparable to $t$ makes Hertz itself the wrong reference (Section 7); use a bonded-layer benchmark or a converged FE reference there instead: $p(r) = p_0\sqrt{1-\frac{r^{2}}{a^{2}}},\qquad p_0 = \frac{3F}{2\pi a^{2}}$
Choose the candidate force set and measurement set; assemble $\mathbf{H}$ column-by-column with unit-load solves — or exploit the fully-observed shortcut $\mathbf{F} = \mathbf{K}\mathbf{U}$ when every retained DOF carries a displacement estimate (and say plainly which components are tracked and which are completed from assumed geometry — Section 8).
Compute the SVD; inspect $\kappa(\mathbf{H})$ and the Picard plot with representative measured data to count recoverable modes before promising a spatial resolution.
Impose physical constraints: nonnegative normal forces (NNLS), optional friction cone, total-force equality when an F/T sensor is in the loop.
Choose $\mathbf{L}$: identity for point-load problems, surface Laplacian for smooth pressure fields.
Select $\lambda$ by the discrepancy principle wherever $\delta$ is measurable — record the unloaded sensor and read the marker-tracking repeatability — cross-checked against the L-curve corner and GCV on calibration presses. Disagreement beyond an order of magnitude is a diagnostic, not a tie to break: GCV is the one that marker-interpolation noise quietly breaks (Section 5). Then freeze it.
Precompute the regularized pseudoinverse; per-frame cost is one matrix–vector product for the unconstrained solution — this is what makes real-time possible. Keep step 5's NNLS or friction cone and each frame is instead a small warm-started QP/SOCP over the active set: still real-time at these sizes, no longer one product (Section 5).
Validate at two levels: resultants against F/T across the workspace, distributions against known indenter geometries.
Monitor drift: gel wear and aging change $E$ and thickness; recalibrate on schedule.

Why it matters — step 2 is the one people skip

Steps 3–8 all inherit their meaning from step 2. If the forward model has the wrong boundary conditions or a mischaracterized modulus, the inverse solution converges beautifully to the answer of the wrong question, and no residual, L-curve, or cross-validation will tell you. Only a forward-model check against independent physical measurement can.

The Failure Catalogue

Linear model at deep indentation — hyperelastic reality, linear pipeline (Sections 7, 8).
Wrong boundary conditions — edges assumed fixed that are not (Section 9).
Marker-interpolation noise is correlated and quietly breaks GCV (Section 5).
Regularization spreads the contact patch and underestimates peak pressure, while the resultant typically survives much better — resultants live in the low-order, well-conditioned modes that the filter barely touches. Note this is a tendency, not a conservation law: ordinary Tikhonov does not conserve $\mathbf{c}^{\mathsf T}\mathbf{g}$ exactly unless you impose it as an equality constraint — report all three numbers, not just the flattering one (Section 4).
The per-point stiffness fallacy $f_i \propto u_i$ (Section 1).
Skipping projection/refraction compensation (Section 7).
Half-space (Boussinesq–Cerruti) formulas on thin bonded pads — invalid unless contact scale $\ll$ pad thickness (Section 7).
Conflating the two iFEMs — load reconstruction vs. Tessler-style shape sensing (Section 1).
Trusting a single $\lambda$ criterion (Section 5).
Validating only resultants and claiming distribution accuracy (Section 8).
Reporting constrained-DOF entries of $\mathbf{K}\mathbf{u}$ as contact forces — they are reactions (Section 2).
Plotting consistent nodal forces as a pressure map without the tributary-area conversion (Section 2).

Flashcards

References

Tikhonov, A. N.; Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. Book (reviewed in Mathematics of Computation 32(144), 1978, and SIAM Review 21(2), 1979).
Golub, Gene H.; Heath, Michael; Wahba, Grace (1979). Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter. Technometrics 21(2), pp. 215–223. doi:10.1080/00401706.1979.10489751
Hansen, Per Christian (1992). Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review 34(4), pp. 561–580. doi:10.1137/1034115
Hansen, Per Christian (1998). Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. SIAM Monographs on Mathematical Modeling and Computation. doi:10.1137/1.9780898719697
Kamiyama, Kazuto; Vlack, Kevin; Mizota, Terukazu; Kajimoto, Hiroyuki; Kawakami, Naoki; Tachi, Susumu (2005). Vision-Based Sensor for Real-Time Measuring of Surface Traction Fields. IEEE Computer Graphics and Applications 25(1), pp. 68–75. doi:10.1109/mcg.2005.27
Jacquelin, E.; Bennani, A.; Hamelin, P. (2003). Force reconstruction: analysis and regularization of a deconvolution problem. Journal of Sound and Vibration 265(1), pp. 81–107. doi:10.1016/s0022-460x(02)01441-4
Maincon, Philippe; Barnardo, Celeste; Larsen, Carl M. (2008). VIV Force Estimation Using Inverse FEM. ASME 27th International Conference on Offshore Mechanics and Arctic Engineering (OMAE 2008), pp. 673–681. doi:10.1115/omae2008-57325
Lourens, E.; Reynders, E.; De Roeck, G.; Degrande, G.; Lombaert, G. (2012). An augmented Kalman filter for force identification in structural dynamics. Mechanical Systems and Signal Processing 27, pp. 446–460. doi:10.1016/j.ymssp.2011.09.025
Sanchez, J.; Benaroya, H. (2014). Review of force reconstruction techniques. Journal of Sound and Vibration 333(14), pp. 2999–3018. doi:10.1016/j.jsv.2014.02.025
Ma, Daolin; Donlon, Elliott; Dong, Siyuan; Rodriguez, Alberto (2019). Dense Tactile Force Estimation using GelSlim and inverse FEM. 2019 International Conference on Robotics and Automation (ICRA), pp. 5418–5424. doi:10.1109/icra.2019.8794113
Sferrazza, Carmelo; Wahlsten, Adam; Trueeb, Camill; D'Andrea, Raffaello (2019). Ground Truth Force Distribution for Learning-Based Tactile Sensing: A Finite Element Approach. IEEE Access, vol. 7, pp. 173438–173449 (arXiv:1909.04000). doi:10.1109/access.2019.2956882
Liu, Ruixue; Dobriban, Edgar; Hou, Zhichao; Qian, Kun (2022). Dynamic Load Identification for Mechanical Systems: A Review. Archives of Computational Methods in Engineering 29(2), pp. 831–863. doi:10.1007/s11831-021-09594-7
Yang, Hongji; Jiang, Jinhui; Chen, Guoping; Zhao, Jiamin (2023). Dynamic load identification based on deep convolution neural network. Mechanical Systems and Signal Processing, vol. 185, 109757. doi:10.1016/j.ymssp.2022.109757
Shahidzadeh, Amir-Hossein; Caddeo, Gabriele; Alapati, Koushik; Natale, Lorenzo; Fermuller, Cornelia; Aloimonos, Yiannis (2024). FeelAnyForce: Estimating Contact Force Feedback from Tactile Sensation for Vision-Based Tactile Sensors. arXiv:2410.02048.
Zhao, Can; Liu, Jin; Ma, Daolin (2025). iFEM2.0: Dense 3-D Contact Force Field Reconstruction and Assessment for Vision-Based Tactile Sensors. IEEE Transactions on Robotics, vol. 41, pp. 289–305 (online November 2024). doi:10.1109/tro.2024.3502197
Li, Kelu; Xiao, Longfei; Liu, Mingyue; Kou, Yufeng (2025). A distributed dynamic load identification approach for thin plates based on inverse Finite Element Method and radial basis function fitting via strain response. Engineering Structures, vol. 322, 119072. doi:10.1016/j.engstruct.2024.119072

← 7. Material Identification 9. Reference Configuration & Prestress →