Regularization & Optimization

AbstractRegularizer

A regularization functional R(σ) with value, coefficient gradient and a Gauss–Newton Hessian model, e.g. TikhonovRegularizer or TotalVariationRegularizer. Added to a data objective by RegularizedObjective.

TikhonovRegularizer(G; reference = 0)

Quadratic regularizer R(σ) = ½ (σ - σ₀)ᵀ G (σ - σ₀) for a symmetric positive semidefinite Gram matrix G and a reference conductivity σ₀ (scalar or vector). With a discretization, TikhonovRegularizer(disc; kind, reference) assembles G for the :L2, :H1semi, :H1 or :jump (piecewise constants) norm.

Interface: objective_value, value_and_gradient!, gauss_newton_hessian.

TotalVariationRegularizer(disc; ε = 1e-3)

Smoothed total variation R(σ) = TV_ε(σ) of the conductivity (see total_variation): facet jumps Σ_F |F| √((σ_K - σ_K')² + ε²) for piecewise constants, ∫ √(|∇σ|² + ε²) for continuous σ. The Gauss–Newton Hessian model is the lagged diffusivity matrix (the jump or stiffness matrix weighted by 1/√(… + ε²) at the current σ).

ε = 0 gives the exact (non-smooth) total variation. Its proximal operator (prox!) is computed exactly by the Chambolle–Pock method, but it has no gradient where σ is constant, so use it with ProximalGradient or ADMM rather than with gradient-based methods.

gauss_newton_hessian(reg, σ)

Symmetric positive semidefinite Hessian model of the regularizer at σ (sparse), used by GaussNewton: the Gram matrix for TikhonovRegularizer, the lagged diffusivity matrix H(σ) with ∇R(σ) = H(σ) σ for TotalVariationRegularizer.

RegularizedObjective(data, α₁ => R₁, α₂ => R₂, …)

J(σ) = J_data(σ) + Σₖ αₖ Rₖ(σ) for a data objective (e.g. AdjointStateObjective, KohnVogeliusObjective) and regularizers Rₖ with weights αₖ ≥ 0. The data objective must deliver coefficient gradients (gradient = CoefficientGradient()); choose the gradient representation in the optimiser instead (e.g. LBFGS(riesz = L2Gradient(mats))). GaussNewton needs a least-squares data objective with residuals and Jacobians.

prox!(z, reg, v, ρ; weights = nothing, lower = nothing, upper = nothing, tol = nothing,
      maxiter = nothing)
prox(reg, v, ρ; kwargs...)

Proximal operator argmin R(z) + ρ/2 Σᵢ wᵢ (zᵢ - vᵢ)² subject to lower ≤ z ≤ upper, in the diagonal metric w = weights (default: Euclidean; pass lumped_mass for the discrete L² metric). Exact for TikhonovRegularizer without bounds and for the non-smooth TotalVariationRegularizer with ε = 0 (Chambolle–Pock); iterative for other smooth regularizers. A ProximalMap wraps user-defined maps such as denoisers.

For the iterative TV prox, tol is an absolute tolerance on the primal–dual gap, which bounds ρ/2 ‖z - z*‖²_w (default: round-off level), and maxiter caps the iterations. The proximal methods pass tolerances tied to their own progress; other proximal maps ignore both.

ProximalMap(f!)

A regularizer given only by its proximal map f!(z, v, ρ) (write argmin R(z) + ρ/2 ‖z - v‖² into z), e.g. a denoiser for plug-and-play priors. Bounds are applied by clamping afterwards and the metric weights are ignored. Its value is reported as 0.

lumped_mass(disc)

Row sums of the mass matrix of the σ space (cell areas for piecewise constants), the diagonal metric of proximal steps. Back end contract.

AbstractOptimizer

Minimization method for minimize: GradientDescent, LBFGS, GaussNewton.

minimize(obj, σ₀, method; lower = nothing, upper = nothing, maxiter = 100, gtol = 1e-6,
         ftol = 0, ftarget = -Inf, callback = nothing, verbose = false)

Minimize the objective obj starting from σ₀ with method (GradientDescent, LBFGS, GaussNewton) subject to lower ≤ σ ≤ upper (scalars or vectors; nothing = unbounded). Returns an OptimizationState.

Stopping criteria:

• gtol: projected gradient norm ≤ gtol × its initial value;
• ftol: relative decrease of the objective in one iteration ≤ ftol;
• ftarget: objective ≤ ftarget, e.g. the discrepancy principle J ≤ τ² δ²/2 for a (whitened) noise level δ and τ slightly above 1;
• maxiter iterations; callback(state) returning true.

The objective must deliver coefficient gradients; the Riesz map (gradient representation) is an option of the method.

OptimizationState

Result and state of minimize: σ (current iterate), g (coefficient gradient), value, iteration, nevals (objective evaluations), status (:gtol, :ftol, :ftarget, :maxiter, :callback, :linesearch, :running), converged, and history (one entry (value, gnorm, step, nevals) per iterate, starting with σ₀; gnorm is the norm of the projected gradient, step the norm of the last change of σ).

GradientDescent(; riesz = CoefficientGradient())

Gradient descent σ ← P(σ - t R g) with the Riesz map riesz (e.g. L2Gradient(mats)), Barzilai–Borwein step sizes t = sᵀy / yᵀRy and projected Armijo backtracking.

LBFGS(; memory = 10, riesz = CoefficientGradient())

Limited-memory BFGS with memory curvature pairs and initial inverse Hessian γR for the Riesz map riesz. With bounds, the quasi-Newton direction is restricted to the free variables and the step is taken along the projection arc (projected L-BFGS); the memory is reset whenever the direction fails to be a descent direction.

GaussNewton(; damping = :lm, λ = nothing, scaling = :identity, linear_solver = :auto,
            cg_rtol = 1e-2, cg_itmax = 200)

Gauss–Newton method for least-squares objectives (AdjointStateObjective, possibly wrapped in a RegularizedObjective; regularizers contribute their gauss_newton_hessian).

• damping = :lm: Levenberg–Marquardt with adaptive λ (initial λ relative to the largest diagonal entry of JᵀJ + ΣαH, default 1e-3);
• damping = :linesearch: fixed damping λ (default 1e-8, relative) and projected Armijo backtracking.
• scaling: damping matrix D: :identity, :marquardt (diag(JᵀJ + ΣαH), floored), :sensitivity (diag(‖J eⱼ‖), the column norms of the Jacobian, floored: the geometric mean of the two, which damps the over-sensitive parameters near the electrodes without over-amplifying the insensitive interior) or a symmetric positive definite matrix (e.g. the σ mass matrix FEMatrices(disc).M_σ). Updated at every iterate.
• linear_solver: :dense (form the nσ × nσ matrix), :woodbury (sparse Cholesky of ΣαH + λD and an m × m system for the m residuals), :auto, or :cg: matrix-free, the Jacobian is only applied (one linearized and one adjoint solve per pattern and CG iteration) and never stored, for problems whose Jacobian does not fit into memory. cg_rtol and cg_itmax control the inexact CG solve.

ProximalGradient(α => G; weights = nothing, accelerated = true)

Proximal gradient method (FISTA) for F(σ) + α G(σ): F is the objective passed to minimize, G a regularizer with a proximal operator (prox!), e.g. the exact total variation TotalVariationRegularizer(disc; ε = 0) or a ProximalMap. Bounds are part of the proximal step. weights: diagonal metric (default Euclidean; use lumped_mass for the L² metric). accelerated: monotone FISTA with restarts.

The history records F + αG and the norm of the gradient mapping.

ADMM(α => G; ρ = 1, inner = LBFGS(), inner_maxiter = 10, inner_gtol = 1e-8,
     weights = nothing, adaptive = true)

Alternating direction method of multipliers for F(σ) + α G(σ) (see the wiki articles ADMM and Nested ADMM Reconstruction): the data step minimizes F + ρ/2 ‖x - (z - u)‖²_w with inner_maxiter iterations of the method inner (any AbstractOptimizer, bounds applied), the regularizer step is the proximal operator of G (prox!, e.g. exact TV or a ProximalMap denoiser for plug-and-play priors). adaptive: residual balancing of ρ. Returns the z iterate (feasible for G and the bounds).

The history records F(x) + α G(z), the primal residual ‖x - z‖_w as gnorm, and the change of z as step. Convergence (gtol): both residuals below gtol relative to the iterates and the scaled dual variable.

jacobian_svd(obj, θ; weights = nothing)

Singular value decomposition J = U S Vᵀ W of the Jacobian of the least-squares objective obj at θ (thin: min(m, n) singular values, in decreasing order). The columns of V are the parameter modes, orthonormal in the metric W = Diagonal(weights) (Euclidean by default), ordered by how well the data determine them. Returns a named tuple (U, s, V, r, w) with the residual r and the metric weights w (ones by default).

For a ParametrizedObjective over pixels, V[:, 1:k] is the data-optimal subspace of dimension k (see jacobian_basis).

jacobian_basis(obj, θ, k; weights = nothing)

The k leading parameter modes V[:, 1:k] of jacobian_svd: the subspace best determined by the data at θ, e.g. for SubspaceParametrization(pixels, jacobian_basis(obj, θ, k)).

TruncatedGaussNewton(; rank = nothing, rtol = nothing, weights = nothing)

Gauss–Newton method with truncated-SVD steps: at each iterate, the step uses only the leading singular modes of the Jacobian (see jacobian_svd), those with index ≤ rank and singular value ≥ rtol · s₁ (at least one of the two must be given). The step is the minimum-norm step in the metric Diagonal(weights), followed by projected backtracking; bounds are handled on the free set. weights = :sensitivity uses the column norms ‖J eⱼ‖ of the current Jacobian (floored), which keeps the steps from concentrating at the electrodes.

Truncation is the regulariser: the objective must be a pure least-squares objective (no RegularizedObjective). Stop by the discrepancy principle, ftarget = discrepancy_target(obj, noise).

sensitivity_map(obj, θ; weights = nothing)

Sensitivity of the data to every parameter at θ: the column norms ‖J eⱼ‖ of the Jacobian of the least-squares objective obj, divided by weights (e.g. pixel areas or the lumped mass, to get a density independent of the element size).

resolution_map(obj, θ; rank = nothing, rtol = nothing, λ = nothing, weights = nothing)
resolution_map(js; rank = nothing, rtol = nothing, λ = nothing)

Diagonal of the model resolution matrix R = V F Vᵀ W of the linearized reconstruction at θ, in [0, 1]: 1 where the data determine the parameter, 0 where it is left to the initial guess or prior. The filter F is either a truncation, with the modes of index ≤ rank and singular value ≥ rtol s₁ kept (as in TruncatedGaussNewton), or Tikhonov damping sᵢ² / (sᵢ² + λ s₁²) (Levenberg–Marquardt with damping matrix W). js is a precomputed jacobian_svd.

posterior_std(obj, θ; noise, prior_std, weights = nothing)
posterior_std(js; noise, prior_std)

Pointwise standard deviation of the linearized Gaussian posterior at θ, for the prior θ ~ N(θ₀, prior_std² W⁻¹) (W = Diagonal(weights), identity by default: independent parameters with standard deviation prior_std) and residual noise N(0, η² I): sqrt.(diag((JᵀJ / η² + W / prior_std²)⁻¹)). It equals the prior standard deviation times sqrt(1 - R), R the Tikhonov resolution_map with λ s₁² = η² / prior_std².

noise is η in units of the (whitened) residual, or a noise model (RelativeGaussianNoise, GaussianNoise) for an AdjointStateObjective (possibly parametrized), converted to the mean residual variance η² = 2 discrepancy_target(obj, noise; τ = 1) / n_residual(obj).