Source-linked AI summary
A regularizing iterative ensemble Kalman method for PDE-constrained inverse problems
Marco A. Iglesias
TL;DR
PDE-constrained inverse problems can be unstable, while many iterative regularization methods require derivatives that black-box forward models may not expose. This paper develops a derivative-free ensemble Kalman regularization method and finds numerically that sufficiently large ensembles can inherit regularizing behavior, including in a Darcy example.
Problem
PDE-constrained inverse problems are ill-posed, and derivative requirements can limit the use of iterative regularization when forward-model derivatives are unavailable.
Source: §1 · p. 2, §2.1 · p. 5Method
The method approximates the regularizing LM scheme with ensemble covariances, an LM-like regularization-parameter choice, and discrepancy-based early termination.
Source: §1 · p. 1, (3) Analysis step. Update each… · p. 15Results
With sufficiently large ensembles, the method numerically inherits LM regularizing properties; for Darcy flow, Ne = 150 yielded stable, reasonably accurate estimates in about 12 iterations with ρ = 0.7 and M = 100 measurements.
Source: §6 · pp. 40, 42Takeaways & Limitations
The derivative-free method is demonstrated on Darcy and EIT identification and shape identification using a level-set parameterization.
Source: §6 · pp. 41–42Convergence theory for the proposed method remains open, and stability and accuracy depend on selecting an appropriate initial ensemble.
Source: §1.1 · p. 4, §6 · p. 42
Abstract
from arXiv · showhide
We introduce a derivative-free computational framework for approximating solutions to nonlinear PDE-constrained inverse problems. The aim is to merge ideas from iterative regularization with ensemble Kalman methods from Bayesian inference to develop a derivative-free stable method easy to implement in applications where the PDE (forward) model is only accessible as a black box. The method can be derived as an approximation of the regularizing Levenberg-Marquardt (LM) scheme [14] in which the derivative of the forward operator and its adjoint are replaced with empirical covariances from an ensemble of elements from the admissible space of solutions. The resulting ensemble method consists of an update formula that is applied to each ensemble member and that has a regularization parameter selected in a similar fashion to the one in the LM scheme. Moreover, an early termination of the scheme is proposed according to a discrepancy principle-type of criterion. The proposed method can be also viewed as a regularizing version of standard Kalman approaches which are often unstable unless ad-hoc fixes, such as covariance localization, are implemented. We provide a numerical investigation of the conditions under which the proposed method inherits the regularizing properties of the LM scheme of [14]. More concretely, we study the effect of ensemble size, number of measurements, selection of initial ensemble and tunable parameters on the performance of the method. The numerical investigation is carried out with synthetic experiments on two model inverse problems: (i) identification of conductivity on a Darcy flow model and (ii) electrical impedance tomography with the complete electrode model. We further demonstrate the potential application of the method in solving shape identification problems by means of a level-set approach for the parameterization of unknown geometries.
1. Introduction
The paper proposes a derivative-free regularizing ensemble Kalman method for PDE-constrained identification, combining ensemble approximations with iterative regularization. It investigates the method numerically while noting that its convergence theory remains open.
1. Introduction
- PDE-constrained inverse problems are ill-posed: small data perturbations can have uncontrolled effects on parameter estimates, motivating stable approximations controlled by the noise level.Iterative regularization aims for estimates that converge to a solution as the noise tends to zero.Source: §1 · p. 2
1.1. Contribution of this work
- Many iterative regularization methods require forward-map derivatives and adjoints, which may be unavailable or inaccessible in modular black-box simulation software.The proposed framework uses the forward map in black-box fashion.Source: §1.1 · p. 2
- The method replaces forward-map derivatives and adjoints with ensemble covariances, selects a regularization parameter in an LM-like way, and stops using a discrepancy-principle criterion.It is presented as both an approximation of the regularizing LM scheme and a regularizing version of standard iterative Kalman methods.Source: §1.1 · pp. 2–3
- The work tests stability and accuracy on Darcy-flow conductivity identification and electrical impedance tomography, and explores shape identification with a level-set parameterization.For the shape application, the paper combines the ensemble method with the level-set approach of.Source: §1.1 · p. 3
- The paper reports that convergence theory for the proposed method remains an open problem, and presents numerical investigation of its convergence and regularization properties.The investigation varies ensemble size, number of measurements, initial ensemble, and tunable parameters.Source: §1.1 · p. 4
2. PDE-constrained inverse problems
The section formulates PDE-constrained identification as fitting observations with a forward operator, then introduces an ensemble Kalman scheme with regularization and discrepancy-based stopping. It describes Darcy-flow and EIT test models and gives preliminary evidence that regularization and early stopping improve estimates.
2.1. General framework
- The identification framework maps PDE parameters to observable quantities in a finite-dimensional observation space, while allowing an infinite-dimensional parameter space.The stated application scales are about 10^3–10^4 observations and parameter dimensions greater than 10^6 after discretization.Source: §2.1 · p. 4
- PDE parameter identification seeks u such that G(u) matches noisy observations, but compact forward operators can make the inverse problem unstable under unregularized minimization.Many iterative regularization methods require derivatives that may not be available in commercial forward-simulation software.Source: §2.1 · p. 5
2.2. The regularizing ensemble Kalman method
- The proposed algorithm updates each member of an initial ensemble and terminates using a discrepancy principle applied to the average of the ensemble's model predictions.That averaged prediction approximates G of the ensemble mean to first order, but is not the estimate's data misfit.Source: (1) Prediction step. Evaluate · p. 6
2.3. Test Models
- The Darcy test identifies log hydraulic conductivity from noisy pointwise head measurements, while the EIT test identifies conductivity from electrode voltage measurements with known contact impedances.The EIT forward model is the Complete Electrode Model.Source: §2.3.1 · p. 7, §2.3.2 · pp. 7–8
2.4. Regularizing properties of the ensemble Kalman method
- For nonlinear forward models, the standard Kalman choice α_n = 1 can lead to instabilities, whereas the proposed method selects α_n and stops according to regularization principles.The preliminary Darcy experiment uses the same initial ensemble for the regularized and unregularized methods.Source: §2.4 · pp. 8–9
- For the Darcy example, the unregularized estimate becomes uncontrolled, while the regularized estimate is more accurate when stopped near the noise level.The regularized method decreases data misfit and relative error slowly; continuing beyond the stopping point can let error increase.Source: §2.4 · p. 9, Figure 1 · p. 10
3. Properties and computational aspects of the regularizing ensemble Kalman method
The proposed method's ensemble mean remains in the span of the initial ensemble, making initialization a design choice that shapes the estimates the method can produce.
3.1. The initial ensemble
- At every iteration, the ensemble mean lies in the subspace generated by the initial ensemble members.This invariance makes the starting ensemble a constraint on the estimate's possible values.Source: §3.1 · p. 10
- The initial ensemble is therefore crucial to performance; prior knowledge of admissible solutions can guide its construction.The paper suggests truncated bases or samples with desired regularity, including Gaussian or Besov priors.Source: §3.1 · p. 11
3.2. The Regularizing LM scheme
This section presents the regularizing Levenberg–Marquardt scheme and its parameter and stopping conditions, then gives a finite-observation-space reformulation useful for deriving the ensemble approximation.
3.2. The Regularizing LM scheme
- The LM update uses the Fréchet derivative and its adjoint, with regularization weighted by C^-1, an a priori operator enforcing regularity on the parameter space.In the Bayesian framework, C^-1 is the inverse covariance operator of a prior distribution.Source: §3.2 · p. 12
- The LM regularization parameter must satisfy a residual inequality, and the scheme terminates when the discrepancy condition is met with τ > 1/ρ.The parameter ρ is selected a priori in (0, 1).Source: §3.2 · p. 12
- When the observation space is finite-dimensional and derivative functionals are linearly independent, the LM update can be rewritten with an inversion in observation space.After discretization, this inversion can have negligible cost when the number of observations is small.Source: §3.2 · pp. 12–13
3.3. The proposed method as an approximation of the regularizing LM scheme
This section derives the ensemble update by replacing derivative-based terms with ensemble-based approximations. It culminates in a formula for the ensemble mean and the proposed ensemble Kalman scheme.
3.3. The proposed method as an approximation of the regularizing LM scheme
- Terms in the update are replaced by approximations derived from the ensemble covariance construction.The supplied passages introduce a first-order approximation around the ensemble mean before describing the replacements.Source: §3.3 · p. 13
- The proposed scheme combines covariance-based approximations with an ensemble-mean update and corresponding choices of regularization and stopping criteria.The section concludes by presenting the resulting method as Algorithm 1.Source: §3.3 · pp. 13–14
3.4. The proposed method as a sequence of linear inverse problems
The method is recast using augmented variables and projection operators. The resulting scheme is organized into prediction, discrepancy-principle, and covariance steps.
3.4. The proposed method as a sequence of linear inverse problems
- In augmented variables, the algorithm uses projections onto the data and parameter spaces and proceeds through prediction, discrepancy checking, and covariance construction.The passages identify H and H⊥ as the projections and list the prediction step, discrepancy principle, and sample covariance.Source: §3.4 · pp. 14–15
- The augmented formulation uses the product space Z = X × Y and represents the algorithm in terms of its projections onto the two component spaces.The initial ensemble is specified in X, with ρ ∈ (0, 1) and τ > 1/ρ.Source: §3.4 · p. 14
(3) Analysis step. Update each ensemble member as follows
The analysis step has a Tikhonov-regularization interpretation in an ensemble-defined subspace. Its regularization parameter is selected by a discrepancy principle, while ρ and τ shape stabilization and stopping.
(3) Analysis step. Update each ensemble member as follows
- For a finite linearly independent ensemble, Cn is positive definite on its range, allowing the analysis update to be expressed through the regularized inverse problem.The proposition states the linear-independence assumption for the ensemble at each iteration.Source: (3) Analysis step. Update each… · pp. 15–16
- The ensemble mean is the Tikhonov-regularized solution of an artificial linear inverse problem in the ensemble subspace Zn.The regularization uses the norm induced by the sample covariance Cn.Source: (3) Analysis step. Update each… · pp. 15–16
- The artificial inverse problem's noise level matches that of the original nonlinear problem, and the discrepancy rule applies to that artificial problem at every iteration.This interpretation is stated for iterations before convergence.Source: (3) Analysis step. Update each… · p. 16
- At each iteration, α is selected so the ensemble mean's artificial inverse problem honors the discrepancy principle.The analysis identifies this α-selection rule as central to the regularizing method.Source: (3) Analysis step. Update each… · pp. 15–16
- Choosing ρ near 1 yields a large α and can stabilize the estimate, whereas smaller ρ calls for earlier stopping to prevent data overfitting.The passage cautions that ρ near 1 may also be detrimental to performance; τ is chosen just above 1/ρ.Source: (3) Analysis step. Update each… · p. 17
3.5. The connection with the Bayesian framework
The method is connected to Bayesian ensemble smoothing, which gives a posterior approximation in the linear-Gaussian case. For nonlinear forward maps, the smoother is only an approximation, motivating the paper's regularizing deterministic approach.
3.5. The connection with the Bayesian framework
- With a linear forward operator, the updated ensemble approximates the Bayesian posterior as ensemble size tends to infinity, and its mean tends to the maximum a posteriori estimate.The stated result assumes Gaussian noise and a Gaussian prior.Source: §3.5 · p. 17
- For nonlinear forward maps, the posterior is generally non-Gaussian, and the ensemble smoother is only an approximation with no convergence theory known to the authors.The passage also notes that straightforward application can yield inaccurate posterior approximations.Source: §3.5 · p. 17
- The paper adapts an intrinsically Bayesian method with regularizing properties for deterministic inverse problems using a derivative-free computational framework.This frames the method's role relative to the Bayesian ensemble smoother.Source: §3.5 · p. 18
3.6. Small ensemble size ill-conditioning
When the ensemble is smaller than the number of measurements, the covariance matrix is singular and the method can be ill-conditioned. The authors motivate iterative regularization as a way to stabilize ensemble Kalman methods beyond standard ad-hoc fixes.
3.6. Small ensemble size ill-conditioning
- When the ensemble size is smaller than the number of measurements, the covariance matrix is singular.The passage identifies this setting as Ne < M.Source: §3.6 · p. 18
- Even with a strictly positive definite regularized inverse, some choices of Γ can produce a large conditioning number.The passage compares this ill-conditioning to the ill-posedness of the Fréchet derivative in the regularizing LM scheme.Source: §3.6 · p. 18
- The authors argue that iterative regularization can address this ill-conditioning, while standard fixes lack a mathematical framework for general implementation guidelines.Covariance localization and inflation are cited as standard ad-hoc fixes.Source: §3.6 · p. 18
3.7. Computational aspects of the scheme
Under the assumption that the number of measurements is small relative to the discretized parameter-space dimension, forward-model evaluations dominate the scheme's computational cost. Feasibility depends on obtaining reasonable accuracy with a few ensemble members and iterations, though some cases require larger ensembles.
3.7. Computational aspects of the scheme
- The main computational cost is n⋆Ne forward-model evaluations, with n⋆ determined by the stopping criterion.Construction of Cuw_n and inversion of the measurement-space matrix are negligible for the measurement counts considered.Source: §3.7 · p. 18
- The scheme is computationally feasible if a few ensemble members achieve reasonable accuracy in only a few iterations.The authors note that some cases require relatively large ensembles for stable computations with reasonable accuracy.Source: §3.7 · p. 18
4. Numerical Experiments
Numerical experiments on Darcy flow and EIT examine how ensemble size, measurement count, stopping and parameter choices, and the initial ensemble affect the method’s regularizing behavior and stability. The EIT experiments further assess how initial-ensemble correlation length relates to accuracy and representation of the truth’s spatial structure.
4.2.2. The selection of the initial ensemble
- Because estimates lie in the initial ensemble’s span, different ensembles can produce different estimates; the EIT study examined five estimates for a prescribed truth unrelated to the ensemble distribution.The EIT examples used five initial ensembles with Ne = 100 and ρ = 0.5.Source: §4.1 · p. 20, §4.2 · p. 29, §4.2.1 · p. 31
4.1.1. Effect of the ensemble size Ne.
- For Darcy flow, the critical ensemble size was Ne = 125; at Ne ≥ 125, stopping near η/ρ generally preceded error growth, while larger ensembles improved accuracy.For Ne ≥ 200, relative error did not increase even when the misfit fell below η/ρ, though a slow increase appeared below the noise level.Source: §4.1.1 · pp. 21–22, Figure 5 · p. 22
4.1.2. Effect of the number of measurements.
- For M ≤ 100, the experiments indicated Ne > M was needed for stabilization, whereas for M > 100 some Ne < M choices, including Ne = 169 for M = 225, were regularized.For M = 400, Ne = 196 and Ne = 225 also avoided error growth before the misfit reached η/ρ.Source: §4.1.2 · pp. 23–24, Figure 7 · p. 25, Figure 8 · p. 26
4.1.3. The tunable parameters ρ.
- For Darcy flow, ρ = 0.7 provided a reasonable accuracy-cost balance; lower values led to rapid error growth below η/ρ, while increasing ρ slightly improved accuracy at greater cost.These experiments used Ne = 150 and 100 measurements, with stopping set near τ = 1/ρ.Source: §4.1.3 · pp. 24, 27, Figure 9 · p. 27
4.1.4. Regularizing properties of the proposed method in the small-noise limit
- At the smallest tested Darcy noise level, 0.5%, Ne = 150 was needed to stabilize computations, supporting inherited regularization for sufficiently large ensembles.The algorithm used the discrepancy stopping criterion with τ = 1/ρ.Source: §4.1.4 · pp. 27–28, Figure 10 · p. 28
4.1.5. Comparison with the regularizing LM scheme.
- In Darcy-flow experiments, the ensemble Kalman method approximated regularizing LM constrained to the initial-ensemble subspace and minimized data misfit stably.The comparison used 40 initial ensembles of size Ne = 150.Source: §4.1.5 · p. 28, Figure 11 · p. 29
4.2. Test Model II: EIT
- For EIT with M = 240, a critical ensemble size of Ne = 75 was sufficient for stabilization, below the measurement count and consistent with addressing the small-ensemble effect.The experiments used ρ = 0.5 and stopping near τ = 1/ρ.Source: §4.2 · pp. 29–30, Figure 14 · p. 31
4.2.1. Ensemble size Ne, tunable parameter ρ and small noise limit.
- For EIT, ρ = 0.5 balanced accuracy and computational cost, while the small-noise experiments showed stable convergence as η → 0.With Ne = 100, stopping near τ = 1/ρ was essential because relative error rose abruptly after misfit fell below η/ρ.Source: §4.2.1 · p. 31, Figure 16 · p. 32
4.2.2. The selection of the initial ensemble
- The study varies L, which controls sample spatial correlation and regularity, while holding θ fixed; smaller L produces less spatially correlated initial-ensemble samples.Figure 18 illustrates Gaussian-prior draws for several L values.Source: §4.2.2 · p. 32, Figure 18 · p. 33
- Across 40 experiments, smaller correlation lengths generally stabilized relative error and reduced data misfit but degraded accuracy; at L = 1 and L = 0.05, error rose before misfit crossed η/ρ.The initial-ensemble choice therefore affects the scheme’s regularizing properties.Source: §4.2.2 · pp. 33–34, Figure 19 · p. 34
- Because updates remain in the initial ensemble’s invariance subspace, low-correlation members cannot represent the truth’s ~0.3-diameter high-conductivity inclusions, which smoother fields with similar correlation identify better.The authors suggest further research to identify initial-ensemble regularity parameters that better characterize the truth.Source: §4.2.2 · p. 34
5. Applications to geometric inverse problems
The paper applies a regularizing ensemble Kalman method to shape identification, parameterizing unknown interfaces with level-set functions. In synthetic geologic-facies and EIT examples, results depend on ensemble size and initialization choices, including the level-set correlation length.
5. Applications to geometric inverse problems
- The target is a sharp interface separating regions with different parameter values, represented as an unknown region within the PDE domain.The motivating subsurface example has conductivity values that may differ by several orders of magnitude across geologic facies.Source: §5 · p. 35
- The level-set ensemble method estimates shapes without shape derivatives or level-set evolution, and its estimates do not become flat.The update with controlled regularization moves the shape; the authors connect avoidance of flattening to the initial ensemble and the method's invariant-subspace property.Source: §5 · p. 36
5.1. Estimation of geologic facies
- The geologic-facies initialization uses a covariance with an intrinsic, nonadjustable correlation length.The example assumes the two conductivity values are known and treats the interface between them as unknown.Source: §5.1 · p. 36
- For N_e > 150, average data misfit and conductivity error decrease; with N_e = 150, five conductivity estimates visually agree with the truth.The experiments use 40 initial ensembles for the ensemble-size comparison; the five estimates use initial ensembles generated from the same Gaussian measure.Source: §5.1 · pp. 37–38, Figure 24 · p. 39
5.2. EIT. Sharp interfaces.
- For EIT, 0.06 ≤ L ≤ 0.1 gives the lowest average truth-relative error; L ≥ 0.2 produces larger data-misfit fluctuations.The authors associate the larger fluctuations with high conductivity appearing close to the electrodes.Source: §5.2 · pp. 39–40
- When the optimal EIT correlation length is unknown, estimating L within the method is beyond this manuscript's scope.The authors report that varying L changes the spatial correlation of the interface and the resulting estimates.Source: §5.2 · pp. 39–40
6. Conclusions
The conclusions report stable estimates for sufficiently large ensembles across Darcy flow and EIT, while emphasizing the method’s derivative-free operation and its application to level-set shape identification. Stability and accuracy depend on ensemble size and, for EIT, the initial ensemble’s spatial correlation length.
6. Conclusions
- With Ne=150 and M=100, Darcy estimates were stable and reasonably accurate on average within 12 iterations at ρ=0.7, using 1,800 forward simulations.The authors describe this computational cost as comparable with gradient-based techniques.Source: §6 · p. 40
- The method requires no derivatives of the forward map, and its computational cost and accuracy with large enough ensembles are described as competitive with variational iterative regularization methods.The Darcy example used 1,800 forward simulations.Source: §6 · pp. 40–41
- For the 16-electrode EIT setting, Ne=100 was sufficient for stability across a wide class of tunable parameters, but accuracy depended on the initial ensemble’s spatial correlation length.The authors report better approximations when the initial samples’ correlation length resembles that of the true conductivity.Source: §6 · pp. 41–42
- The regularized ensemble method was successfully used for level-set shape identification, producing stable computations of unknown high-conductivity regions with a reasonable ensemble size.The approach induces controlled shape motion without solving the level-set equation for shape evolution.Source: §6 · p. 42
- The numerical investigation reports that sufficiently large ensembles let the method inherit the regularizing properties of the LM scheme and address instability associated with small ensembles.The authors present this as numerical evidence that iterative-regularization ideas can stabilize ensemble methods for inverse estimates.Source: §6 · pp. 40, 42
7. Appendix
The appendix proves Lemma 3.1 for an infinite-dimensional parameter space by expressing the finite set of observation derivatives through functionals and representers, then deriving the update expression.
7.1. Proof of Lemma 3.1
- Each observation derivative D_mG(u_n) acts as a scalar functional on X and is represented through the inner product with DG*_m(u_n).The proof uses these representers to handle the infinite-dimensional case where the finite-dimensional matrix argument does not apply.Source: §7.1 · p. 43
- In the representer decomposition, the component b is orthogonal to every r_j, so D_mG(u_n)b vanishes for each observation.This follows from the stated orthogonality and the functional representation.Source: §7.1 · p. 43
- The proof uses linear independence of the observation derivatives to obtain an invertible matrix and then derives expression (24) from the minimizer of (69).The final step writes u_{n+1}−u_n=v and uses expressions (64)–(65) and (70).Source: §7.1 · p. 44