3. Inverse Problems & Regularization

Hadamard well-posedness, least squares, SVD, discrete Picard condition, Tikhonov and TSVD, L-curve/GCV/discrepancy, Bayesian MAP, inverse problems in solid mechanics

1. Forward vs. Inverse: Which Direction Is Hard?

In a forward problem you are given the causes — loads, material parameters, geometry, boundary conditions — and you compute the effects: displacements, strains, sensor readings. Module 2 was entirely forward. The inverse problem reverses the inference: you observe noisy effects and must estimate the causes that produced them.

The asymmetry between the two directions is physical, not computational. The forward maps that matter here — loads and boundary data to remote displacement or strain — are smoothing: the underlying elliptic boundary-value problems damp high-spatial-frequency components of the input. (This is a property of these maps, not a law of continuum mechanics: hyperbolic wave propagation transports singularities rather than erasing them, and the Dirichlet-to-Neumann map amplifies high boundary frequencies instead of damping them.) By Saint-Venant-type reasoning, two statically equivalent but locally different load distributions produce nearly identical fields a short distance away. The forward map compresses distinct causes into nearly indistinguishable effects.

Concretely, in vision-based tactile sensing: a contact pressure profile on a soft skin's surface produces marker displacements measured at depth $h$ below the surface. Pressure detail finer than roughly the embedding depth is blurred away, so two distinct contact profiles map to nearly indistinguishable camera data. Inverting means un-smoothing — dividing by numbers close to zero — which amplifies measurement noise without bound.

The same structure appears everywhere:

• Image deblurring — blur (forward) is easy and stable; deblurring is unstable.
• Hearing a drum — computing the sound of a known drum is forward; inferring the drum's shape from its sound is inverse.
• Gravimetric prospecting — computing surface gravity from a known density distribution is forward; mapping subsurface density from surface gravity is inverse.

Three framings recur throughout this guide, distinguished by the kind of unknown (not by difficulty — the third row is the best-conditioned of the three, as Section 12 explains):

2. The Observation Model: One Equation, Four Families

Before any machinery, fix the vocabulary that the whole guide runs on. Every inverse method in Modules 4–12 is an instance of a single observation model, and every method is fully characterized by how it answers four questions: what is measured, what is unknown, what is assumed known, and what is unobservable.

Naming which argument of $\mathcal{H}$ is the primary target — the rest either assumed known or carried as jointly estimated nuisance unknowns, as in the all-at-once formulations of Section 12 — generates four distinct inverse families. They share the machinery of this module, but they are separate problems with different conditioning, different priors, and different literatures; none is a variant of another:

A notation collision, flagged once: in the spine equation $\mathbf{b}$ denotes systematic bias; in this module's linear algebra, following the numerical-analysis convention, $\mathbf{b}$ is the data vector in $\mathbf{A}\mathbf{x}=\mathbf{b}$. Context disambiguates — the spine's bias term will not reappear until the model-error caveat in Section 11 and the validation discussion in Section 13, where structured residuals reveal it.

Rank, Nullspace, and What "Unobservable" Means

Linearize $\mathcal{H}$ about an operating point and collect the derivatives into a matrix $\mathbf{A}$ (the sensitivity or observation matrix). Two subspaces of the unknown space now decide everything:

• rank$(\mathbf{A})$ counts the independently observable directions — the number of distinct pieces of information your sensor set can deliver, at most $\min(m,n)$.
• The nullspace, $\{\mathbf{v} : \mathbf{A}\mathbf{v} = \mathbf{0}\}$, contains the unknown components your sensors cannot see at all. Adding any nullspace vector to a solution changes nothing in the predicted data.

Mechanics supplies canonical nullspaces: strain gauges are blind to rigid-body motion — rigid translation and rotation produce exactly zero strain. Statically equivalent interior load distributions are the neighbouring but distinct case: by the Saint-Venant argument of Section 1 they produce nearly identical boundary data, so they span a near-nullspace direction, not an exact one — formally observable, practically lost in noise. Two terms of art, used consistently from here on: observability asks whether a state or field ($\mathbf{u}$, $\mathbf{X}_0$-type unknowns) is determined by the measurements; identifiability asks the same for parameters ($\mathbf{p}$-type unknowns). For linear problems both reduce to the nullspace of $\mathbf{A}$; for nonlinear parameter dependence, the rank of the Jacobian governs identifiability only locally.

Between "observable" and "unobservable" lies the practically important middle ground: directions that $\mathbf{A}$ maps to tiny but nonzero outputs. They are formally observable and practically lost in noise. Quantifying that middle ground is exactly what the SVD does in Section 5.

3. Hadamard Well-Posedness and Why Elasticity Inversion Breaks It

A problem is well-posed in the modern Hadamard sense if it passes three tests; a problem failing any one of them is ill-posed. The attribution needs a footnote: Hadamard's own 1923 wording defined a problem as “correctly set” when the accessory conditions “determine one and only one solution” — existence and uniqueness. He treated instability separately, in the same lectures, through the Cauchy problem for Laplace's equation; continuity of the inverse became the standard third condition later. The three-condition list below is the definition the field now uses, not a quotation from 1923.

Each condition fails in a concrete, mechanical way:

Existence. Noisy data $\mathbf{b}$ generally lies outside the range of $\mathbf{A}$ in an overdetermined system — no $\mathbf{x}$ satisfies $\mathbf{A}\mathbf{x}=\mathbf{b}$ exactly. Least squares (Section 4) repairs this by asking for the closest achievable fit instead.

Uniqueness. A nullspace: strain gauges cannot see rigid-body motion. (Statically equivalent interior load distributions are the near-null case of Section 2 — formally distinguishable, practically not.)

Stability. The canonical example is differentiation of noisy data. Perturb $g(x)$ by a tiny high-frequency wiggle and watch the derivative:

As $k$ grows, the data perturbation stays at size $\varepsilon$ while the answer perturbation $\varepsilon k$ blows up without bound. The same mechanism drives the backward heat equation, where recovering the initial state amplifies mode $k$ by $e^{k^2 T}$.

The Structural Reason: Compact Operators

The instability is not bad luck; it is built into the physics. When the forward map is a compact operator between infinite-dimensional spaces — the situation for a first-kind Fredholm integral equation, which is the continuous model of deconvolution and of load-to-strain maps —

a bounded inverse cannot exist, so condition (3) fails structurally, not accidentally (Engl, Hanke & Neubauer 1996). Discretization does not remove the pathology; it converts the unbounded inverse into a finite matrix whose condition number grows without bound as the mesh refines.

In elasticity specifically, the Cauchy problem — reconstructing tractions and displacements on an inaccessible boundary from overdetermined data on an accessible one — is severely ill-posed: it is the elasticity analogue of Hadamard's own counterexample, the Cauchy problem for Laplace's equation, where a boundary perturbation of spatial frequency $k$ is amplified like $e^{k}$ (Hadamard 1923). Interior parameter-field identification from full-field data is typically milder; that comparison is this guide's reading of the two literatures. The severity spectrum matters practically and returns in Section 12.

History then reversed itself: Hadamard considered ill-posed problems physically meaningless; Tikhonov's school (Tikhonov & Arsenin 1977) showed they are meaningful once prior information restores stability. The repair kit maps one-to-one onto the failures:

4. Linear Least Squares and the Normal Equations

The workhorse setting: $\mathbf{A} \in \mathbb{R}^{m\times n}$ with $m > n$ and full column rank. Since noisy $\mathbf{b}$ generically lies outside the $n$-dimensional subspace range$(\mathbf{A}) \subset \mathbb{R}^m$, minimize the squared residual instead. Setting the gradient of the quadratic $\|\mathbf{A}\mathbf{x}-\mathbf{b}\|_2^2$ to zero gives the normal equations:

Geometric reading: $\mathbf{A}\hat{\mathbf{x}}$ is the orthogonal projection of $\mathbf{b}$ onto range$(\mathbf{A})$; the residual is orthogonal to every column of $\mathbf{A}$. The solution $\hat{\mathbf{x}} = (\mathbf{A}^\mathsf{T}\mathbf{A})^{-1}\mathbf{A}^\mathsf{T}\mathbf{b} = \mathbf{A}^\dagger\mathbf{b}$ defines the Moore–Penrose pseudoinverse for full-rank $\mathbf{A}$.

Conditioning: How Many Digits Survive?

The spectral condition number measures noise amplification through a linear solve (singular values $\sigma_i$ formally defined in Section 5):

Rule of thumb: a linear solve loses about $\log_{10}\kappa$ decimal digits; double precision holds about 16. The precise perturbation facts, needed later: for consistent systems the relative error in $\mathbf{x}$ is bounded by $\kappa(\mathbf{A})$ times the relative error in $\mathbf{b}$; for least-squares problems with a large residual, an additional $\kappa(\mathbf{A})^2$ term enters even for backward-stable algorithms. So "least squares amplifies errors by $\kappa$" is only true in the small-residual regime, and only for the right algorithm — always state which algorithm and which regime.

Ill-conditioning hides in innocuous-looking matrices. The $n\times n$ Hilbert matrix — a classic severely ill-conditioned Gram matrix with clean analytic entries — has $\kappa(\mathbf{H}_5) \approx 4.8\times10^5$ and $\kappa(\mathbf{H}_{10}) \approx 1.6\times10^{13}$, so by the rule above only about three decimal digits survive a solve at $n=10$. Note that this is the severe case in Section 5's taxonomy, not the mild one: $\kappa(\mathbf{H}_n)$ grows exponentially in $n$, and by $n=12$ ($\kappa \approx 1.6\times10^{16}$) double precision really is exhausted.

Foreshadowing the rest of the module: for well-conditioned $\mathbf{A}$, least squares is numerically stable and needs no regularization for stability — which is the only thing conditioning can settle. It does not follow that shrinkage cannot pay: at low signal-to-noise a little bias can still lower the mean-squared error even on a flat spectrum, and the lab below demonstrates exactly that at $w=0.5$ ($\kappa \approx 1.7$), where raising the noise slider to 10% produces an interior minimum in the error curve. So read a flat spectrum as “no stabilizing regularization required” and decide any optional denoising on SNR or validation grounds, not from $\kappa$. Regularization earns its place as a necessity whenever the Picard diagnosis of Section 6 shows the data coefficients sinking into the noise before the singular values do — a large $\kappa$ is the warning sign, not the criterion, because a right-hand side whose small-singular-direction coefficients stay above the noise can be inverted as it stands. And which regularization depends on the spectrum's shape: a gapped spectrum calls for truncation at the numerical rank, a gradually decaying one for a parameter you must choose. That is Section 5's diagnosis.

Weighted Least Squares, Whitening, and Correlated Noise

Ordinary least squares implicitly assumes every measurement deserves equal trust. When the noise has covariance $\mathbf{C} = \operatorname{cov}(\mathbf{e})$, the statistically correct functional weights the residual by the inverse covariance (the Gaussian likelihood of Section 11 makes this exact; Kaipio & Somersalo 2005):

Equivalently, whiten: multiply both $\mathbf{A}$ and $\mathbf{b}$ by $\mathbf{C}^{-1/2}$ and solve an ordinary least-squares problem whose noise is white with unit variance. Dividing each row by its sensor channel's noise standard deviation is the diagonal special case. After whitening — and only after — the statements of Sections 5–9 that depend on whiteness hold: the Picard noise floor is flat, GCV's assumptions are met, residual-whiteness checks mean what they say. Morozov's discrepancy principle (Section 9.1) is the exception on that list and is worth separating out: it consumes a noise norm $\delta$ measured in whatever metric you are actually fitting in, so it survives correlated noise as long as you know that norm — use $\|\mathbf{e}\|$ for an unweighted fit, $\|\mathbf{C}^{-1/2}\mathbf{e}\|$ after covariance weighting. Whiteness is how $\delta$ is usually estimated, not a precondition for the principle.

5. The SVD: Anatomy of an Ill-Conditioned Problem

The singular value decomposition is the diagnostic instrument for everything in this module. It rotates the problem into a basis where the forward map is a plain diagonal scaling:

with orthonormal left/right singular vectors $\mathbf{u}_i$, $\mathbf{v}_i$ and $\sigma_1 \ge \sigma_2 \ge \cdots \ge \sigma_n \ge 0$. Write $r = \operatorname{rank}(\mathbf{A})$ for the number of nonzero singular values; the sums below run to $r$, not $n$, because $\sigma_i = 0$ is admitted (the nullspace of Section 2) and $1/\sigma_i$ is then undefined. For the full-column-rank setting of Section 4, $r = n$. In this basis the minimum-norm least-squares solution is an explicit expansion, and splitting the data $\mathbf{b} = \mathbf{b}_{\mathrm{exact}} + \mathbf{e}$ exposes exactly where inversion goes wrong:

Two facts drive everything that follows:

Fact 1 — white noise is flat in the singular basis. The exact coefficients $\mathbf{u}_i^\mathsf{T}\mathbf{b}_{\mathrm{exact}}$ decay with $i$ (Section 6), but for white noise of standard deviation $s$, every $|\mathbf{u}_i^\mathsf{T}\mathbf{e}| \approx s$ regardless of $i$ — an orthogonal rotation preserves white noise. So component $i$ of the naive solution carries a noise term of size $s/\sigma_i$: and once $\sigma_i$ is small that term is enormous. Be careful about the exact criterion, though: what decides whether component $i$ is recoverable is not $\sigma_i$ against $s$, but the exact data coefficient $|\mathbf{u}_i^\mathsf{T}\mathbf{b}_{\mathrm{exact}}|$ against $s$. Once the signal in a direction has sunk below the noise, inverting it returns amplified noise no matter how the singular value compares.

Fact 2 — small-$\sigma$ singular vectors are usually oscillatory. For discretized smoothing kernels — the deconvolution and load-to-strain maps this guide is about — the singular vectors associated with small $\sigma_i$ have increasingly many sign changes, the matrix version of "the forward map kills high frequencies". Hence the amplified garbage is high-frequency, which is why an unregularized solution looks like noise rather than like a slightly-wrong answer. Treat this as a property of that class of operators, not as a consequence of ill-conditioning: a diagonal $\mathbf{A}$ with decaying diagonal entries has exactly the gradually decaying spectrum of a discrete ill-posed problem, yet its singular vectors are the coordinate axes — no sign changes, no spatial length scale at all. Check that your own $\mathbf{v}_i$ really do order spatial frequency before reading resolution off them (Section 8).

A Taxonomy of Singular Spectra

Plot $\sigma_i$ against $i$ on a semilog axis and the problem type reads off directly (Hansen 1998; Hansen 2010):

For calibration of expectations: deconvolution with a Gaussian kernel is severely ill-posed — its spectrum decays like a Gaussian in frequency — while many elasticity identification operators are mildly to moderately ill-posed.

6. The Discrete Picard Condition

Small singular values alone do not doom a problem — what matters is how fast the data's SVD coefficients fall relative to the singular values. That comparison is the discrete Picard condition, and its plot is the single most informative diagnostic you can produce before attempting any inversion.

Reading a Picard Plot

On a semilog axis versus index $i$, plot three things: $\sigma_i$, the data coefficients $|\mathbf{u}_i^\mathsf{T}\mathbf{b}|$, and their ratio $|\mathbf{u}_i^\mathsf{T}\mathbf{b}|/\sigma_i$. With real (noisy) data the coefficients track the decaying exact coefficients only until they hit the noise floor at the per-sample noise standard deviation $s$, where they level off into a horizontal band — because $|\mathbf{u}_i^\mathsf{T}\mathbf{e}| \approx s$ for all $i$ (Section 5):

Define the crossing index $i^*$: the index at which the measured coefficients stop tracking the exact ones and flatten onto the noise plateau — equivalently, where $|\mathbf{u}_i^\mathsf{T}\mathbf{b}_{\mathrm{exact}}|$ sinks to $s$. Read $i^*$ as a budget, not a partition. Below it most components carry recoverable signal; above it most are noise, the ratio grows like $s/\sigma_i$, and inverting those directions amplifies noise — but $i^*$ is the centre of a transition band, not a hard boundary, and individual exact coefficients do resurface above the floor past it (Section 8's selective-filtering caveat). In the lab configuration quoted at the end of this paragraph the exact coefficients cross the noise floor five times before settling below it for good. One caveat on any such count: the lab's singular values come in degenerate cosine/sine pairs, and within a repeated pair the split of the energy between the two coefficients depends on which basis of that 2-D subspace you happen to hold. The crossing count and the count above the floor are therefore basis-dependent at the margin; what is invariant is the energy of the whole pair, $\|\mathbf{U}_{\text{pair}}^\mathsf{T}\mathbf{b}\|_2$, against a two-dimensional noise threshold. Quote the counts as a reading of a particular plot, not as an operator invariant. Note also what $i^*$ is not: the index where $\sigma_i$ itself drops below $s$. That crossing is a different (and later) event — in the figure above the coefficients flatten at $i^* \approx 16$ while $\sigma_i$ only reaches the noise level near $i = 25$ — and using it would overstate the information you actually hold. Recoverability is decided by the data coefficients, not by the singular values alone. $i^*$ is the honest measure of the information content of the experiment — often brutally small, and, once the unknown grid is fine enough to resolve what the sensors can see, stable under further refinement of the unknowns. It is not a discretization invariant in the strict sense: changing $n$ changes $\mathbf{A}$, its left singular vectors and the coefficients $\mathbf{u}_i^\mathsf{T}\mathbf{b}$, so while the grid is still too coarse the count can move by an index or two. What refinement never does is create new experimental information. The count depends on the kernel, the true signal and the noise level, so quote it with its configuration: in the 64-point periodic Gaussian deconvolution of the lab below (PSF width $w = 1.5$ grid units, 1% relative noise), only about two dozen of the 64 exact-data coefficients sit above the noise floor.

The experiment-design reading: adding sensors, moving them closer to the source, or reducing noise raises $i^*$. Adding unknowns does not — past the point where the grid resolves the measurement-supported subspace, refinement leaves $i^*$ where it was.

This diagnosis motivates both regularization families: TSVD keeps components $i \le k \approx i^*$ (Section 8); Tikhonov rolls components off smoothly around $\sigma_i \approx \lambda$ (Section 7).

7. Tikhonov Regularization

Tikhonov's repair (Tikhonov & Arsenin 1977): stop asking for the best data fit, and minimize a penalized functional that trades fidelity against a prior notion of "reasonable":

The seminorm matrix $\mathbf{L}$ and reference $\mathbf{x}_0$ encode the prior. Standard form: $\mathbf{L}=\mathbf{I}$, $\mathbf{x}_0=\mathbf{0}$ — penalize magnitude. Common general forms: $\mathbf{L}$ = discrete first-derivative (bidiagonal) or Laplacian (tridiagonal) operator — penalize roughness, appropriate when the unknown field (e.g. a contact pressure distribution under a smooth indenter) is expected smooth. Setting the gradient to zero gives the regularized normal equations, uniquely solvable provided $\mathbf{L}$ actually covers what $\mathbf{A}$ cannot see, i.e. $\ker\mathbf{A}\cap\ker\mathbf{L}=\{\mathbf{0}\}$. That kernel condition buys uniqueness, and by itself nothing more. Conditioning is a separate question with a separate answer: for $\mathbf{L}=\mathbf{I}$ the eigenvalues of the left-hand side are $\sigma_i^2+\lambda^2$, so for $\lambda$ not too small the condition number drops to $(\sigma_1^2+\lambda^2)/(\sigma_n^2+\lambda^2)$ — a genuine and large improvement. For a general $\mathbf{L}$ it depends on the scaling of $\mathbf{L}$ and on the generalized spectrum of the pair $(\mathbf{A},\mathbf{L})$, and it can come out far worse than $\mathbf{A}^\mathsf{T}\mathbf{A}$ alone: a badly scaled $\mathbf{L}$ with a trivial kernel is enough to do it, since the penalty then adds its own dynamic range to the matrix. Compute $\kappa$ of the regularized system; do not assume it improved. The kernel condition itself is automatic for $\mathbf{L}=\mathbf{I}$ but not for derivative penalties: a first-difference $\mathbf{L}$ has the constants in its nullspace, and what a discrete Laplacian has depends entirely on its boundary rows — a boundary-free 1-D second-difference stencil (the $(n-2)\times n$ operator) has constants and linears, a Neumann or graph Laplacian has constants only, and a Dirichlet Laplacian is nonsingular. If $\mathbf{A}$ is also blind to whatever that nullspace turns out to contain, no $\lambda$ will pin those modes down.

Filter Factors: Tikhonov in the Singular Basis

For $\mathbf{L}=\mathbf{I}$ and $\mathbf{x}_0=\mathbf{0}$ — standard form, assumed for the rest of this section and for the bias–noise split below — the SVD diagonalizes everything, and the Tikhonov solution is a filtered version of the least-squares expansion:

The filter factors interpolate smoothly between keep and kill: $f_i \approx 1$ for $\sigma_i \gg \lambda$ (component kept), $f_i \approx \sigma_i^2/\lambda^2 \to 0$ for $\sigma_i \ll \lambda$ (component suppressed). Tikhonov is a low-pass filter in the singular basis, with $\lambda$ the cutoff singular value. The two forms above are algebraically identical wherever $\sigma_i > 0$; the right-hand one is what you implement, since it never divides by a small (or zero) $\sigma_i$.

With a nonzero reference the same algebra applies to the shifted problem — filter the residual data $\mathbf{b}-\mathbf{A}\mathbf{x}_0$ and add the reference back:

and every statement below then reads with $\mathbf{x}_{\mathrm{exact}}$ replaced by $\mathbf{x}_{\mathrm{exact}}-\mathbf{x}_0$: the bias is bias towards the reference, not towards zero.

The Bias–Noise Trade-Off Is a U-Shaped Curve

Substituting $\mathbf{b} = \mathbf{A}\mathbf{x}_{\mathrm{exact}} + \mathbf{e}$ into the filtered expansion splits the error into two opposing terms:

Numerical Practice

Never form $\mathbf{A}^\mathsf{T}\mathbf{A}$ (Section 4). Solve the augmented least-squares system by QR instead — it is an ordinary least-squares problem with the penalty rows stacked beneath the data rows:

Caveat for general $\mathbf{L}$: the SVD filter-factor formulas above hold for $\mathbf{L}=\mathbf{I}$ only. For derivative-type $\mathbf{L}$ the correct analysis uses the generalized SVD (GSVD), and components in null$(\mathbf{L})$ — the constant mode for a first-derivative $\mathbf{L}$ — pass through completely unregularized. Verify that $\mathbf{A}$ constrains them, or the "regularized" solution can still contain an arbitrarily wild unpenalized component. The same caveat decides where $\lambda$ lives: the $[\sigma_n,\sigma_1]$ sweep range of Section 9 is calibrated to $\mathbf{L}=\mathbf{I}$, because replacing $\mathbf{L}$ by $c\mathbf{L}$ rescales the useful $\lambda$ by $1/c$ while leaving every $\sigma_i(\mathbf{A})$ untouched — $\mathbf{x}_\lambda$ for the pair $(\lambda,\mathbf{L})$ is exactly $\mathbf{x}_{\lambda/c}$ for $(\lambda/c,c\mathbf{L})$. For a general $\mathbf{L}$, bracket the sweep with the generalized singular value ratios $\alpha_i/\beta_i$ of the pair $(\mathbf{A},\mathbf{L})$, or widen it adaptively until the residual and the penalty have each gone flat at their respective ends.

8. Truncated SVD and Friends

Truncated SVD (TSVD) is the blunt-scissors alternative: keep the first $k$ terms of the SVD expansion and discard the rest.

The truncation index $k$ and the parameter $\lambda$ are analogous controls — each sets where the expansion stops being trusted — and $\lambda \approx \sigma_k$ is the usual rule of thumb for translating one into the other (Hansen 2010). Read it as a hint, not an identity: at $\lambda = \sigma_k$ Tikhonov keeps that component at $f_k = 1/2$ while TSVD keeps it whole, so wherever component $k$ carries real weight the two solutions part company. The pairing is dependable only across a genuine spectral gap, $\sigma_{k+1} \ll \lambda \ll \sigma_k$, where the ambiguous component is negligible either way; under gradual decay there is no one-to-one mapping, and how close the two land also depends on the data coefficients, not on spectral decay alone. In the lab below at $w=3$, TSVD at $k=24$ and Tikhonov at $\lambda=\sigma_{24}$ differ by 64% of $\|\mathbf{x}_\lambda\|$. Read that as a warning about degeneracy rather than as a test of the safe regime, because the safe regime is exactly what this configuration does not satisfy. The lab's periodic kernel has doubly degenerate singular values — one cosine and one sine per frequency — so $\sigma_{24} = \sigma_{25}$ exactly ($1.9389\times10^{-3}$), $k=24$ splits that degenerate pair, and $\sigma_{k+1} \ll \lambda$ fails outright. Truncating the whole block instead ($k=25$) barely helps: still 64%, because consecutive blocks here sit only a factor $\approx 3$ apart, which is nowhere near the $\ll$ the rule asks for. Fast decay is not the same thing as a gap.

Cost. A full SVD is $O(mn^2)$ — trivial for the $n \lesssim 10^3$ unknowns typical of per-frame reconstruction problems. For large 3D field identification ($n \sim 10^{5+}$) one uses iterative methods (CGLS/LSQR), which exploit semi-convergence: early iterations pick up the large-$\sigma$ components first, so the iteration count itself acts as a regularization parameter. Stop early, using the same criteria as Section 9 — running the iteration to convergence faithfully reproduces the noise-dominated naive solution.

Selective-filtering caveat. TSVD assumes the coefficient decay is roughly monotone. If the Picard plot shows an exact-coefficient spike beyond $k$, the hard cutoff deletes real signal that Tikhonov's tail would have partially kept.

Resolution: What the Filtered Solution Can Show

Apply any filtered method to noise-free data and the algebra of Section 7 gives its reconstruction of the truth as

— the true solution seen through a low-pass filter. For TSVD this is exactly the orthogonal projection onto $\mathrm{span}\{\mathbf{v}_1,\dots,\mathbf{v}_k\}$: features of $\mathbf{x}_{\mathrm{exact}}$ outside the first $k$ singular vectors are gone from this reconstruction. That is a property of the truncation you chose, not of the physics — as long as $\sigma_i > 0$ the data do contain $\mathbf{u}_i^\mathsf{T}\mathbf{b}_{\mathrm{exact}} = \sigma_i\,\mathbf{v}_i^\mathsf{T}\mathbf{x}_{\mathrm{exact}}$, and raising $k$ recovers the component whenever the data are clean enough to support it. Only exact-nullspace directions, $\sigma_i = 0$, are invisible no matter how good the data (Section 2). What the noise level fixes is where raising $k$ stops paying: past $i^*$ the recovered component is mostly amplified noise, so at a given noise level the spatial resolution of an honest reconstruction is capped near the length scale of $\mathbf{v}_{i^*}$ — the Picard index again, now read as a resolution limit. That last reading assumes the singular vectors of your $\mathbf{A}$ genuinely order spatial frequency (Fact 2, Section 5); on an operator where they do not, $i^*$ still counts recoverable components but no longer converts into a length.

9. Choosing the Regularization Parameter

The bias–noise curve of Section 7 must be located blind — $\mathbf{x}_{\mathrm{exact}}$ is not available. Three standard criteria exist, each with an honest failure profile. None is ground truth, and — worth being exact about — none of them actually estimates the parameter-space error you would like to minimize: Morozov matches a residual to the noise level, the L-curve locates a corner in a trade-off plot, and GCV estimates predictive risk in data space. They are computable proxies for the bottom of the U, and they land on it only to the extent that their proxy tracks it.

9.1 Morozov's Discrepancy Principle

If the noise level $\delta = \|\mathbf{e}\|_2$ is known, choose $\lambda$ so the residual matches it — fitting the data below the noise level means fitting noise:

The residual norm is monotone increasing in $\lambda$, so this is a 1-D root find. Read $\tau\delta$ as a selection target — the level you deliberately push the residual up to — and not as a prediction of what a well-fitted residual should equal: a correct model leaves an expected residual below $\delta$, because the fit absorbs part of the noise (Section 13 gives the expression). Pushing the residual up to $\tau\delta$ is precisely why the principle over-smooths. Behavior: theoretically convergent as $\delta \to 0$ (Engl, Hanke & Neubauer 1996); tends to slightly over-smooth; its quality stands or falls with the $\delta$ estimate. It is the best default in engineering settings where sensor noise is characterizable.

9.2 The L-Curve

Plot the trade-off itself (Hansen 1992): a parametric log-log plot of $(\|\mathbf{A}\mathbf{x}_\lambda-\mathbf{b}\|_2,\; \|\mathbf{L}(\mathbf{x}_\lambda-\mathbf{x}_0)\|_2)$ as $\lambda$ sweeps — the ordinate must be the penalty term actually being traded against the residual, which reduces to $\|\mathbf{L}\mathbf{x}_\lambda\|_2$ only in the common case $\mathbf{x}_0=\mathbf{0}$ assumed for the rest of this section. Under-regularized solutions populate a steep vertical branch (solution norm explodes while the residual barely improves); over-regularized ones populate a flat horizontal branch (norm barely drops while the residual grows). Pick the corner — operationalized by Hansen & O'Leary (1993) as the point of maximum curvature of the log-log curve. It needs no noise estimate and is visually informative even when not used for selection. (The figure below marks its corner with a cheaper nearest-corner heuristic; curvature-based selection is the standard implementation, and the two can disagree on a rounded L.)

A documented failure, not anecdote — but read its scope carefully. Hanke (1996) constructs an infinite-dimensional example in which the L-curve-selected parameter vanishes too rapidly as the noise-to-signal ratio in the data goes to zero, so the computed reconstructions do not converge to the true solution at all; he then gives numerical examples showing that similar behaviour appears in discrete ill-posed problems when the exact solution is smooth. That is one proved counterexample plus numerical evidence, not a theorem covering every smooth solution. Note also that this concerns the decay of the solution coefficients, a different quantity from the singular values. Separately, the corner degenerates when the singular values decay slowly: the corner is the transition between the noise-dominated vertical branch and the bias-dominated horizontal branch, and slow (algebraic) decay smears that transition across many indices into a smooth arc. Fast decay does the opposite — it makes the corner sharp, the limiting case being a genuine spectral gap, where every $\lambda$ between $\sigma_{r+1}$ and $\sigma_r$ leaves both norms essentially unchanged and the curve turns a literal right angle. And a computational trap: the corner is defined on the log-log curve — drawn or differentiated in linear axes, the corner moves and often disappears entirely.

9.3 Generalized Cross-Validation

GCV (Golub, Heath & Wahba 1979) chooses $\lambda$ to minimize a rotation-invariant leave-one-out criterion — "how well does the model predict measurements it did not see":

where $\mathbf{A}_\lambda^{\#}$ is the regularized inverse mapping $\mathbf{b} \mapsto \mathbf{x}_\lambda$. The trace costs nothing once the filter factors are known: $\operatorname{trace} = m - \sum_i f_i$. GCV estimates predictive risk without needing a noise estimate. Its optimality result is narrower than the usual one-line quotation of it: Golub, Heath & Wahba (1979) bound the expectation inefficiency — and the parameter it is evaluated at is a minimizer $\tilde\lambda$ of the expected GCV functional $\mathbb{E}V(\lambda)$, not the random minimizer of the one realized GCV curve you compute from your data. Their bound is on $\mathbb{E}T(\tilde\lambda)/\mathbb{E}T(\lambda^0)$, the ratio of the risk at that $\tilde\lambda$ to the smallest attainable risk, and it tends to 1 only when the normalized traces of the influence matrix $\mathbf{H}_\lambda = \mathbf{A}\mathbf{A}_\lambda^{\#}$, namely $\tfrac{1}{m}\operatorname{trace}\mathbf{H}_\lambda$ and $(\operatorname{trace}\mathbf{H}_\lambda)^2/(m\operatorname{trace}\mathbf{H}_\lambda^2)$, both vanish as the problem grows. White Gaussian noise alone does not deliver that: it is a limit statement about a sequence of problems, carrying its own conditions on the spectrum, and never a finite-sample guarantee for the matrix on your desk. Its failure profile: positively correlated noise typically drives GCV to under-smooth badly, because the smooth part of the noise reads as signal — GCV is not assumption-free; it assumes white (uncorrelated, equal-variance) noise, exactly what camera-based measurements with spatial post-processing violate (Section 4). And its minimum is often very flat, demanding careful 1-D minimization rather than eyeballing.

Practice

Sweep 50–100 log-spaced $\lambda$ values across $[\sigma_n, \sigma_1]$ — padded a decade or two beyond each end, so that all three criteria are actually bracketed rather than pinned to an endpoint — then evaluate them on the sweep and inspect the selected solutions side by side. That range is calibrated to $\mathbf{L}=\mathbf{I}$; for a general $\mathbf{L}$ the useful $\lambda$ moves with the scaling of $\mathbf{L}$ (Section 7's caveat), so bracket it from the generalized singular values of $(\mathbf{A},\mathbf{L})$ or by widening until both ends go flat. If a criterion's condition is met nowhere in the range, say so; do not report the boundary value as a pick. If the lower end of the sweep had to be clamped, disclose the clamp and say why: $\sigma_n$ having genuinely reached roundoff is a different situation from $\sigma_n$ merely dropping under some floor your code imposes, and only the first is a statement about the problem. Agreement is evidence; disagreement is diagnosis — a GCV pick far below the discrepancy pick suggests correlated noise; a wandering L-curve corner suggests slowly decaying singular values — a poorly separated signal-to-noise transition. The criteria estimate the optimum, they do not certify it.

10. Nonlinear Problems: Gauss–Newton and Levenberg–Marquardt

The machinery of Sections 4–9 assumed a linear forward map (Sections 1–3 did not: they are written for a general $G$ and $\mathcal{H}$). Material identification breaks that assumption: the response is generally nonlinear in the parameters even when the physics is linear (Section 12). The standard treatment is iterative — linearize, solve a linear least-squares problem for the update, repeat.

Gauss–Newton. At iterate $\mathbf{m}_k$, with Jacobian $\mathbf{J}_k = \partial G/\partial \mathbf{m}$ evaluated at $\mathbf{m}_k$, the step solves the linearized problem:

Levenberg–Marquardt damps the step with a parameter $\mu$, which in normal-equations form reads

— formally a Tikhonov term with $\mathbf{L}=\mathbf{I}$ applied to the update. Large $\mu$ shortens the step toward gradient descent; small $\mu$ recovers Gauss–Newton. In practice $\mu$ is shrunk as iterations succeed and grown when they fail.

The connection to the rest of the module is direct: the Jacobian plays the role of $\mathbf{A}$, so each iteration is a (regularized) linear least-squares problem, and the SVD of $\mathbf{J}$ at the solution governs local identifiability — small singular values of $\mathbf{J}$ mark parameter combinations the data barely constrains. The Section 4 warning applies verbatim: forming $\mathbf{J}^\mathsf{T}\mathbf{J}$ squares the condition number; solve the augmented system $[\mathbf{J};\,\mu\mathbf{I}]$ by QR instead.

11. The Bayesian View: Regularization as a Prior

Reframe the whole construction statistically (Kaipio & Somersalo 2005; Calvetti & Somersalo 2018). Model the noise as Gaussian, $\mathbf{e}\sim\mathcal{N}(\mathbf{0},\sigma^2\mathbf{I})$, and encode prior belief as a Gaussian prior on the unknown. Bayes' theorem multiplies likelihood and prior; taking $-\log$ of the posterior and dropping constants yields exactly the Tikhonov functional.

Three consequences worth internalizing:

• $\lambda$ stops being a fudge factor. $\lambda^2$ is the noise-to-prior variance ratio, so $\lambda$ itself is the ratio of standard deviations — "how much do I trust the sensor relative to my prior about the field". A strong prior (small $\gamma$) means large $\lambda$.
• $\mathbf{L}$ is a statement about spatial correlation. A first-derivative $\mathbf{L}$ says increments are independent Gaussians — smooth fields are a-priori probable. That is precisely right for a pressure field under a smooth indenter, and wrong (over-smoothing) for sharp contact edges — which motivates non-Gaussian alternatives: a Laplace prior gives an $\ell_1$ penalty, a sparsity-promoting MAP that is no longer linear in $\mathbf{b}$.
• The posterior comes with error bars. In the linear-Gaussian setting the posterior is fully Gaussian, so regularized inversion silently carries an uncertainty estimate (details in the collapsible below).

Caveats, stated bluntly: MAP = Tikhonov requires Gaussian noise and a Gaussian prior; a Laplace/sparsity prior breaks it outright, replacing the quadratic penalty with an $\ell_1$ one. A nonlinear forward model and a positivity constraint do not break it — they give nonlinear Tikhonov (Section 10) and constrained Tikhonov respectively, with the same objective. Keep two things apart for the constrained case: with a linear $\mathbf{A}$, positivity-constrained Tikhonov still minimizes a quadratic (it is a bound-constrained quadratic program, solved by projected gradient / active-set methods rather than one factorization), but its solution is no longer a linear function of $\mathbf{b}$, so closed-form filter factors and the Gaussian posterior go, while the objective stays exactly the one above. What linearity is needed for is the Gaussian posterior: only then is the posterior mean the MAP point and the closed-form covariance below valid. And the posterior quantifies uncertainty only within the assumed model: model error — wrong $E$, wrong geometry, the bias term $\mathbf{b}$ of Section 2's spine — is not included in those error bars.

12. A Taxonomy of Inverse Problems in Solid Mechanics

Organize the field by which ingredient of the boundary-value problem is unknown, everything else given. The four-class cut below is this guide's own synthesis, not a scheme lifted from a source: Bonnet & Constantinescu's (2005) review, the standard entry point for classes 2 and 3, states its own scope as the identification of distributions of elastic moduli, of model parameters, and of buried objects such as cracks. Read the four classes alongside Section 2's spine, not as a refinement of it: the two partitions cut differently. Classes 1 and 2 map cleanly onto the spine's $\mathbf{f}$ and $\mathbf{p}$; class 3 lumps together $\mathbf{u}$-type shape sensing, the reference-configuration problem $\mathbf{X}_0$ of Module 9, and defect geometry (which is not an argument of $\mathcal{H}$ at all); and class 4's unknown boundary tractions are $\mathbf{f}$ again, but restricted to an inaccessible surface, which is what makes it the severe case.

Class 1 — Load Identification

Geometry, material, and boundary conditions known; applied forces or pressures unknown. For linear elasticity the response is linear in the loads, so discretization gives

with $\mathbf{A}$ an influence/transfer matrix assembled from the FE model (columns = responses to unit loads) — directly a linear discrete ill-posed problem to which Sections 4–9 apply verbatim. Ill-posedness worsens with sensor-to-load distance (Saint-Venant smoothing, Section 1). Dynamic variants add time deconvolution on top; Liu et al. (2022) review the frequency- and time-domain machinery, the ill-posedness handling, and the ML variants. Contact-force-map reconstruction in tactile skin is squarely this class.

Class 2 — Material / Parameter Identification

Constitutive parameters or spatial modulus fields unknown. The response is generally nonlinear in the parameters even when the physics is linear, so the problem is solved iteratively — each Gauss–Newton step (Section 10) is a regularized linear least-squares problem in the parameter update. For full-field data the method families are FEMU (finite element model updating), the constitutive equation gap, the virtual fields method, the equilibrium gap, and the reciprocity gap (Avril et al. 2008). Römer et al. (2025) organize much of this work with a reduced (a forward solver in the loop) versus all-at-once (field and parameters jointly unknown) distinction. Take that as an organizing lens, not an exhaustive partition: they claim it structures “a large portion” of the literature, and they say explicitly that some methods do not fit exactly into either class — the VFM, which optimizes over the parameter domain alone while substituting measured data for the unknown state, is their named example, with the equilibrium gap method treated as a particular case of it. PINNs land on both sides of that cut in their treatment: the conventional "inverse PINN", which identifies material and network parameters simultaneously, belongs to the all-at-once family, while a parametric PINN — one taking the material parameters as extra network inputs — is a surrogate forward model and therefore supports a reduced formulation. Elastography / elasticity imaging is the field-identification version of this class.

Class 3 — Shape, Configuration, and Defect Identification

Unknown geometry. The hard direction: cracks, voids, inclusions — nonlinear, handled with adjoint fields, level sets, and topological derivatives. The benign direction: shape sensing — reconstructing the deformed configuration from discrete strain measurements. This is where iFEM lives (Module 4): a variational least-squares state-reconstruction problem, much better conditioned than the other classes — but still an inverse problem, with observability that depends on sensor placement.

Class 4 — Boundary-Condition Identification

Unknown tractions or displacements on an inaccessible boundary, from overdetermined data on the accessible part — the Cauchy problem of elasticity. Severely ill-posed: it is the elasticity analogue of Hadamard's Cauchy problem for Laplace's equation, whose boundary perturbations of spatial frequency $k$ are amplified like $e^{k}$ (Hadamard 1923; Section 3). Regularization is indispensable, not optional.

13. Practical Workflow: Diagnosing and Solving Your Inverse Problem

A field recipe consolidating the module. Follow it in order; each step is cheap compared to debugging its omission.

1. Scale first. Whiten the rows: multiply $\mathbf{A}$ and $\mathbf{b}$ by $\mathbf{C}^{-1/2}$ (Section 4). Dividing each row by its sensor channel's noise standard deviation is the diagonal special case — it equalizes variances and leaves every correlation in place, so for camera-derived data with spatially correlated noise it is not whitening and does not license the clean statements that follow. Non-dimensionalize the columns by the change of variables $\mathbf{x} = \mathbf{D}\mathbf{z}$, $\mathbf{D} = \operatorname{diag}(c_1,\dots,c_n)$ with $c_j$ the characteristic magnitude of unknown $j$: solve for $\mathbf{z}$ in $\mathbf{A}\mathbf{D}$ — i.e. multiply column $j$ by $c_j$, not divide — then recover $\mathbf{x}=\mathbf{D}\mathbf{z}$, and transform the penalty operator to $\mathbf{L}\mathbf{D}$ to match. Conditioning and the SVD are not scale-invariant: a $\kappa = 10^9$ that drops to $10^3$ after sensible scaling tells you most of that number was units, not physics. It does not by itself certify the problem well-posed — a discretized compact operator still has a spectrum that decays, and a $\kappa$ that keeps climbing as you refine the mesh is ill-posedness however tidy the units are. Inspect the scaled spectrum and its refinement behaviour before deciding.
2. Compute the SVD of the scaled $\mathbf{A}$ and classify the spectrum (Section 5): flat → ordinary least squares, no stabilizing regularization needed (at low SNR still ask whether shrinkage lowers your validation error); gap → rank-deficient, TSVD at the gap; gradual decay → discrete ill-posed, continue.
3. Picard plot with real data (Section 6); read off the crossing index $i^*$ = number of recoverable components. If $i^*$ is far below the number of unknowns, consider re-parameterizing coarser — fewer basis functions for the pressure field — instead of fighting with heavy regularization.
4. Choose $\mathbf{L}$ from prior physics (Section 7): identity for magnitude, gradient/Laplacian for smoothness; check that null$(\mathbf{L})$ is observable through $\mathbf{A}$.
5. Sweep $\lambda$ over 50–100 log-spaced values spanning $[\sigma_n, \sigma_1]$ padded a decade or two beyond each end (so a criterion cannot silently return an endpoint), and remember that this range assumes $\mathbf{L}=\mathbf{I}$ — a general $\mathbf{L}$ moves the useful $\lambda$ with its own scaling (Section 7). Evaluate discrepancy (if $\delta$ is known), L-curve corner, and GCV (Section 9); compare the selected $\lambda$'s and inspect the corresponding solutions side by side. Report any pick that landed on an endpoint as unbracketed, and disclose any clamp you put on the low end.
6. Validate honestly. (a) Residual check: after whitening, the residual should look like sensor noise — but compare it against what the fit itself does to the noise, not against raw white noise. With hat matrix $\mathbf{H}_\lambda$ the residual is $\mathbf{r} = (\mathbf{I}-\mathbf{H}_\lambda)\mathbf{b}$, so $\operatorname{cov}(\mathbf{r}) = \sigma^2(\mathbf{I}-\mathbf{H}_\lambda)(\mathbf{I}-\mathbf{H}_\lambda)^\mathsf{T} \neq \sigma^2\mathbf{I}$ and $$\mathbb{E}\|\mathbf{r}\|_2^2 = \underbrace{\left\|(\mathbf{I}-\mathbf{H}_\lambda)\mathbf{b}_{\mathrm{exact}}\right\|_2^2}_{\text{regularization bias}} + \sigma^2\operatorname{trace}\!\left[(\mathbf{I}-\mathbf{H}_\lambda)(\mathbf{I}-\mathbf{H}_\lambda)^\mathsf{T}\right], \qquad \operatorname{trace}[\cdot] = m - 2\operatorname{trace}\mathbf{H}_\lambda + \operatorname{trace}\mathbf{H}_\lambda^2$$ for symmetric $\mathbf{H}_\lambda$ — equivalently $\sigma^2\big[(m-n) + \sum_i (1-f_i)^2\big]$ in filter-factor terms. Use that expression, not $\sigma^2(m-\text{effective DOF})$: the "effective degrees of freedom" shortcut, $\sigma^2(m-\operatorname{trace}\mathbf{H}_\lambda)$, is exact only when $\mathbf{H}_\lambda$ is an orthogonal projection (idempotent). That class is wider than full-rank least squares: it covers ordinary least squares of any rank, where $\mathbf{H} = \mathbf{A}\mathbf{A}^\dagger = \mathbf{U}_r\mathbf{U}_r^\mathsf{T}$ and the shortcut gives $\sigma^2(m-r)$ — the familiar $\sigma^2(m-n)$ when $r=n$ — and it covers TSVD, a regularizer on this page, whose $\mathbf{H}_k = \mathbf{U}_k\mathbf{U}_k^\mathsf{T}$ is a projection too. What breaks idempotency is a fractional filter: ordinary Tikhonov's $f_i = \sigma_i^2/(\sigma_i^2+\lambda^2)$ is strictly between 0 and 1, so $\mathbf{H}_\lambda^2 \neq \mathbf{H}_\lambda$, and at usable $\lambda$ the shortcut can be wrong by tens of percent — enough to turn a healthy residual into a false model-error alarm. Note also what this is not: Morozov's $\tau\delta$ (Section 9.1) is a rule for choosing $\lambda$, not a prediction of the post-fit residual. With the bias term small the fitted residual sits below $\delta$, by the amount the fit absorbed; it climbs back above $\delta$ only when the bias term grows, which is what over-smoothing and model error both look like. Structure beyond what the expression above predicts means model error (the bias term $\mathbf{b}$ of Section 2's spine). (b) Withheld data: hold out sensors or whole loading cases from the fit and score the reconstruction on predicting them — a solution that only explains the data it was fitted to has explained nothing. (c) No inverse crimes (box below). (d) Report the bias: at the $\lambda$ you actually selected, Tikhonov attenuates narrow peaks and sharpens nothing — state this wherever peak values (e.g. peak contact pressure) are the quantity of interest. (Under-regularized solutions err the other way, ringing and overshooting; check which regime you are in before reporting a direction.)
7. For repeated real-time inversion (fixed $\mathbf{A}$, streaming $\mathbf{b}$): precompute the SVD, or the QR/Cholesky factorization of the augmented system at the chosen $\lambda$; the per-frame cost drops to two matrix–vector products, $O(mn)$.

Tools note: Hansen's Regularization Tools — the companion package to Hansen (2010) — implement Picard plots, filter factors, the L-curve, GCV, and the standard test problems of this literature (including the deconvolution family used in this page's interactives).

Flashcards

References

• Hansen, P. C.; O’Leary, D. P. (1993). The Use of the L-Curve in the Regularization of Discrete Ill-Posed Problems. SIAM Journal on Scientific Computing 14(6), 1487–1503. doi:10.1137/0914086
• Hadamard, Jacques (1923). Lectures on Cauchy's Problem in Linear Partial Differential Equations. Yale University Press (Silliman Lectures), New Haven.
• Tikhonov, A. N.; Arsenin, V. Y. (1977). Solutions of Ill-Posed Problems. V. H. Winston & Sons / Wiley, Washington DC (Scripta Series in Mathematics).
• Golub, Gene H.; Heath, Michael; Wahba, Grace (1979). Generalized Cross-Validation as a Method for Choosing a Good Ridge Parameter. Technometrics 21(2), 215–223. doi:10.1080/00401706.1979.10489751
• Hansen, Per Christian (1990). The discrete Picard condition for discrete ill-posed problems. BIT Numerical Mathematics 30, 658–672. doi:10.1007/BF01933214
• Hansen, Per Christian (1992). Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review 34(4), 561–580. doi:10.1137/1034115
• Engl, Heinz W.; Hanke, Martin; Neubauer, Andreas (1996). Regularization of Inverse Problems. Kluwer Academic Publishers, Dordrecht (Mathematics and Its Applications, vol. 375). doi:10.1007/978-94-009-1740-8
• Hanke, Martin (1996). Limitations of the L-curve method in ill-posed problems. BIT Numerical Mathematics 36, 287–301. doi:10.1007/BF01731984
• Hansen, Per Christian (1998). Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion. SIAM (Society for Industrial and Applied Mathematics). doi:10.1137/1.9780898719697
• Bonnet, Marc; Constantinescu, Andrei (2005). Inverse problems in elasticity. Inverse Problems 21(2), R1–R50 (topical review). doi:10.1088/0266-5611/21/2/R01
• Kaipio, Jari P.; Somersalo, Erkki (2005). Statistical and Computational Inverse Problems. Springer, Applied Mathematical Sciences (vol. 160). doi:10.1007/b138659
• Avril, Stéphane; Bonnet, Marc; Bretelle, Anne-Sophie; Grédiac, Michel; Hild, François; Ienny, Patrick; Latourte, Félix; Lemosse, Didier; Pagano, Stéphane; Pagnacco, Emmanuel; Pierron, Fabrice (2008). Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements. Experimental Mechanics 48(4), 381–402. doi:10.1007/s11340-008-9148-y
• Hansen, Per Christian (2010). Discrete Inverse Problems: Insight and Algorithms. SIAM, Philadelphia (Fundamentals of Algorithms series). doi:10.1137/1.9780898718836
• Calvetti, Daniela; Somersalo, Erkki (2018). Inverse problems: From regularization to Bayesian inference. WIREs Computational Statistics 10(3), e1427. doi:10.1002/wics.1427
• Liu, Ruixue; Dobriban, Edgar; Hou, Zhichao; Qian, Kun (2022). Dynamic Load Identification for Mechanical Systems: A Review. Archives of Computational Methods in Engineering 29, 831–863. doi:10.1007/s11831-021-09594-7
• Römer, Ulrich; Hartmann, Stefan; Tröger, Jendrik-Alexander; Anton, David; Wessels, Henning; Flaschel, Moritz; De Lorenzis, Laura (2025). Reduced and All-at-Once Approaches for Model Calibration and Discovery in Computational Solid Mechanics. Applied Mechanics Reviews 77(4), 040801. doi:10.1115/1.4066118
• Dasgupta, Agnimitra; Ramaswamy, Harisankar; Murgoitio-Esandi, Javier; Foo, Ken Y.; Li, Runze; Zhou, Qifa; Kennedy, Brendan F.; Oberai, Assad A. (2025). Conditional score-based diffusion models for solving inverse elasticity problems. Computer Methods in Applied Mechanics and Engineering 433, 117425. doi:10.1016/j.cma.2024.117425