Regularization

Sentinel provides 10 regularization methods in three categories: penalty methods (added to the objective), Hessian modifiers (modify the update direction), and post-processing methods (applied after each iteration).

See the Mathematical Reference for complete functional definitions.

Abstract Interface

Sentinel.AbstractRegularization Type

AbstractRegularization

Abstract supertype for all regularization methods.

source

Sentinel.WeightSchedule Type

WeightSchedule

Linear interpolation of a regularization weight from start_weight to end_weight over iterations delay+1 through max_iter.

Before iteration delay, the weight is 0. After max_iter, it stays at end_weight.

Ports the weight ramp-up logic from the Fortran regularization modules.

source

Sentinel.evaluate_weight Function

evaluate_weight(ws::WeightSchedule, iteration::Int) -> Float64

Evaluate the scheduled weight at the given iteration. Returns 0 for iterations ≤ delay, then linearly interpolates.

source

Penalty Methods

Tikhonov

$R(\mathit{\theta })=\frac{1}{2}\alpha \parallel \mathit{\theta }-{\mathit{\theta }}_{\mathrm{ref}}{\parallel }^2$

Sentinel.TikhonovReg Type

TikhonovReg <: AbstractRegularization

Tikhonov (L2) regularization penalizing deviation from reference values.

Fields

• prop_index::Int: which material property (1-based) this applies to

• weight_real::WeightSchedule: weight schedule for real part

• weight_imag::WeightSchedule: weight schedule for imaginary part

• theta_ref_real::Vector{Float64}: (npr,) reference values for real part

• theta_ref_imag::Vector{Float64}: (npi,) reference values for imaginary part

source

Total Variation

$R(\mathit{\theta })=\alpha \sum _e\sum _{\mathrm{gp}}{w}_{\mathrm{gp}}|J|\sqrt{|\mathrm{\nabla }\theta {|}^2+{\epsilon }^2}$

Sentinel.TotalVariationReg Type

TotalVariationReg <: AbstractRegularization

Total Variation regularization with smooth (Huber-like) approximation.

Fields

• prop_index::Int: which material property (1-based)

• weight_real::WeightSchedule: weight schedule for real part

• weight_imag::WeightSchedule: weight schedule for imaginary part

• epsilon::Float64: smoothing parameter δ for √(|∇θ|² + ε²)

• active::Bool: false on first iteration, true after update_weights!

• mtrl_mesh::MaterialMesh: hex8 material mesh for integration

• node_to_elements::Vector{Vector{Int}}: precomputed inverse connectivity

source

Soft Prior

Region-based Laplacian penalty: $R(\mathit{\theta })=\alpha \sum _j(L\theta {)}_j^2$

Sentinel.SoftPriorReg Type

SoftPriorReg <: AbstractRegularization

Soft prior regularization using region-based Laplacian penalty.

Fields

• prop_index::Int: which material property (1-based)

• weight_real::WeightSchedule: weight schedule for real part

• weight_imag::WeightSchedule: weight schedule for imaginary part

• region_info::RegionInfo: region assignment for material mesh nodes

source

Sentinel.RegionInfo Type

RegionInfo

Region assignment information for nodes in a material mesh. Used by soft prior regularization to penalize deviation within regions.

Fields

• node_region::Vector{Int}: (nn,) region ID per node (0 = no region)

• nregions::Int: number of distinct regions

• region_node_count::Vector{Int}: (nregions,) node count per region

• region_node_lists::Vector{Vector{Int}}: (nregions,) node lists per region

source

Exterior Constraint

$R(\mathit{\theta })=\alpha \sum _k[max(0,{\theta }_k-{\theta }_{max}{)}^2+max(0,{\theta }_{min}-{\theta }_k{)}^2]$

Sentinel.ExteriorConstraintReg Type

ExteriorConstraintReg <: AbstractRegularization

Exterior constraint penalty that enforces parameter bounds.

Penalty is zero when parameters are within bounds and grows quadratically with violation magnitude. The scalar s = prop.scalar[1] (real) or prop.scalar[2] (imaginary) converts from stored normalized values to physical units before checking bounds.

Fields

• prop_index::Int: which material property (1-based)

• weight::Float64: penalty weight

• min_real::Float64: minimum allowed physical real value

• max_real::Float64: maximum allowed physical real value

• min_imag::Float64: minimum allowed physical imaginary value

• max_imag::Float64: maximum allowed physical imaginary value

source

Hessian Modifiers

Joachimowicz

Scales the diagonal Hessian: $H_{kk}\leftarrow \overline{H}\cdot ϵ\cdot w$

Sentinel.JoachimowiczReg Type

JoachimowiczReg <: AbstractRegularization

Joachimowicz Hessian regularization.

Adds avg_diag · total_error · w to each diagonal entry of the Hessian for the specified property parameters.

Fields

• prop_index::Int: which material property (1-based)

• weight::WeightSchedule: weight schedule

source

Marquardt

Normalized Levenberg-Marquardt damping: normalize $H$, add $\lambda I$, rescale solution.

Sentinel.MarquardtReg Type

MarquardtReg <: AbstractRegularization

Levenberg-Marquardt Hessian regularization with adaptive weight.

Fields

• prop_index::Int: which material property (1-based)

• current_weight::Float64: current λ value

• start_weight::Float64: initial λ value (for reset)

• min_weight::Float64: minimum λ value (floor)

• delta::Float64: reduction factor (< 1.0) applied on successful steps

• adjust_threshold::Float64: |1−α| threshold for adaptation

source

Sentinel.normalize_and_regularize! Function

normalize_and_regularize!(H, grad_vec, reg::MarquardtReg, grad::MaterialGradient,
                           material::Material) -> Vector{Float64}

Perform the Marquardt normalization and regularization in-place:

1. Compute D_k = √|H[k,k]| for relevant parameters

2. Normalize: H[i,j] /= D_i·D_j, grad[k] /= D_k

3. Add λ to diagonal: H[k,k] += reg.current_weight

Returns the scaling vector D for post-solve rescaling.

source

Sentinel.rescale_solution! Function

rescale_solution!(delta_theta, D, reg::MarquardtReg, grad::MaterialGradient,
                   material::Material)

Rescale the Marquardt solution: Δθ[k] /= D_k for relevant parameters.

source

Post-Processing Methods

Van Houten

Clips power-law exponents to $[{\alpha }_{min},{\alpha }_{max}]$.

Sentinel.VanHoutenReg Type

VanHoutenReg <: AbstractRegularization

Van Houten post-processing regularization.

Clips power-law exponent values prop.rvalue[k, 2] and prop.ivalue[k, 2] to the range [-level, level] for properties with nvpp >= 3.

Fields

• level::Float64: maximum allowed power-law exponent magnitude

source

Sentinel.apply_van_houten! Function

apply_van_houten!(material::Material, reg::VanHoutenReg)

Clip power-law exponents in all properties with nvpp >= 3 to [-level, level].

source

Spatial Filter

Gaussian smoothing with widths $({\sigma }_x,{\sigma }_y,{\sigma }_z)$.

Sentinel.SpatialFilterReg Type

SpatialFilterReg <: AbstractRegularization

Spatial filter (Gaussian smoothing) post-processing regularization.

Fields

• prop_index::Int: which material property (1-based)

• sigma_real_init::Float64: initial Gaussian width for real part

• sigma_real_final::Float64: final Gaussian width for real part

• sigma_imag_init::Float64: initial Gaussian width for imaginary part

• sigma_imag_final::Float64: final Gaussian width for imaginary part

• max_iter::Int: total global iterations (for weight interpolation)

• const_reg_iters::Int: constant regularization iterations at end

• cutoff_multiplier::Float64: neighbor search radius = max_sigma × cutoff_multiplier

• mtrl_mesh::MaterialMesh: material mesh for coordinate lookup

• neighbor_indices::Vector{Vector{Int}}: precomputed neighbor lists per node

• neighbor_distances::Vector{Vector{Float64}}: precomputed distances per node

• region_info::Union{RegionInfo, Nothing}: optional region constraints

• sensitivity_mask::Union{BitVector, Nothing}: optional mask for GP-covered nodes (Fortran nodsense)

source

Sentinel.apply_spatial_filter! Function

apply_spatial_filter!(material::Material, reg::SpatialFilterReg;
                      iteration::Int=1)

Apply Gaussian spatial smoothing to material property values: θ_filtered[i] = Σ_j w_j·θ_j / Σ_j w_j where w_j = exp(−d²/(2σ²)).

Sigma is interpolated between init and final values based on iteration.

source

Bounds (McGarry)

Enforces Poisson ratio bounds via clipping.

Sentinel.BoundsReg Type

BoundsReg <: AbstractRegularization

McGarry bounds post-processing regularization.

Fields

• max_poisson::Float64: maximum allowed Poisson's ratio (e.g., 0.485) for Model 3

• property_bounds::Vector{PropertyBounds}: per-property min/max constraints

source

Sentinel.apply_mcgarry_bounds! Function

apply_mcgarry_bounds!(material::Material, reg::BoundsReg)

Enforce material property bounds (post-processing, no gradient contribution).

1. General min/max clamping: for each constrained property, clamp scalar * rvalue to [min, max] (matching Fortran mcgarrybounds.f90).

2. Model 3 Poisson ratio: clamp λ so that ν ≤ ν_max.

source

Force Balance

Physics-consistency penalty: the equilibrium residual of the measured displacement field under the current property estimate, $R(\mathit{\theta })=\frac{1}{2}w\parallel MK(\mathit{\theta }){\hat{u}}_{\mathrm{meas}}{\parallel }^2$ over non-Dirichlet displacement DOFs. Model 1 (isotropic incompressible) only.

Weights are configured as fractions of the initial data misfit (relative mode, the default for config-driven runs): the fraction is resolved to an absolute weight once at solve start via $w=\mathrm{frac}\cdot D_0/S_0$, which makes the same number portable across meshes and datasets. MGH calibration (benchmark/results/mgh_fb_calibration/): fractions ≤ 0.2 are gentle (data misfit at or slightly better than baseline, per-voxel correlation ≥ 0.99 against the unregularized reconstruction); ≥ 5 over-regularizes. Enabling force balance roughly doubles per-iteration cost (~2.05× on the MGH batched-GPU benchmark, down from 2.26× via exact affine line-search caches for the penalty and gradient — see fb_begin_line_search!); the dominant remaining cost is the per-iteration gradient evaluations, with bounds-clamped line-search trials falling back to full evaluation.

macOS: relative-weight resolution is a one-time iterative solve

Resolving the relative weight (w = frac · D0/S0) needs one global forward solve for D0, once at solve start. On macOS this uses a condensed, factorization-free iterative solve (seconds to a few minutes) instead of a direct factorization, which avoids an upstream Apple-libmalloc crash; the resolved value is accurate and Linux is unaffected (it keeps the faster direct solve). A slow ForceBalanceReg: relative weights resolved, D0/S0 = … step on macOS is the expected one-time cost — see macOS FB / xzone Crash.

Choosing the fraction — noise matters (ground-truth phantom arbitration, benchmark/results/fb_phantom_arbitration/): on equilibrium-consistent (clean / high-SNR / denoised) data, strong force balance genuinely improves accuracy (frac = 1.0: −21 % RMSE vs truth, better contrast recovery). On realistically noisy data it becomes a bias — the equilibrium pull drives the properties to "explain" measurement noise (frac = 1.0 at 2 % noise: +25 % RMSE, flattened contrast). Keep the default frac = 0.1 (neutral) for noisy in-vivo data; reserve larger fractions for high-SNR or denoised acquisitions.

Sentinel.ForceBalanceReg Type

ForceBalanceReg <: AbstractRegularization

Force-balance (equilibrium residual on measured displacements) regularization.

Fields

• weight::WeightSchedule: penalty weight schedule

• relative::Bool: when true, weight holds fractions of the natural scale D0/S0 (initial data misfit over raw equilibrium penalty at the initial material) rather than absolute weights. run_inverse_solve! resolves relative weights to absolute ones once at solve start (one extra global forward solve); a relative instance cannot be bound directly.

• bound context fields (set by bind_context!; untyped because ForwardProblemContext is defined after the regularization includes)

Construct with ForceBalanceReg(weight::WeightSchedule; relative=false) or ForceBalanceReg(w::Float64; relative=false); the context is bound later per zone.

The raw penalty ½‖M·K(θ)·û_meas‖² has force² units and its magnitude is mesh- and dataset-dependent (MGH calibration 2026-07-12: raw penalty 1.3e5 vs data misfit 5.2e-2 — a natural scale of ~4e-7). Relative mode makes the configured number portable: w = 0.2, relative = true means "an equilibrium penalty worth 20% of the initial data misfit", on any mesh and any dataset.

source

Sentinel.bind_context! Function

bind_context!(reg::AbstractRegularization, ctx; kwargs...)

Bind a ForwardProblemContext to a regularization that needs FE-level data (mesh, DOFs, boundary conditions, measured displacements) beyond the Material passed to compute_penalty/add_gradient!. Called from solve_zone! after the context is built. No-op for most regularizations.

source

bind_context!(regs, ctx; kwargs...)

Bind a context to every regularization in a collection (nothing allowed).

source

bind_context!(reg::ForceBalanceReg, ctx;
              gradient_backend=nothing, ka_edata=nothing,
              ka_workspace=nothing)

Bind a zone's ForwardProblemContext to the regularizer so that compute_penalty/add_gradient! can assemble the equilibrium residual. Called from solve_zone!. Throws for unsupported material models.

source

Sentinel.resolve_fb_weights! Function

resolve_fb_weights!(regularizations, ctx, material; K=nothing) -> Union{Float64, Nothing}

Resolve every relative-weighted ForceBalanceReg in regularizations to an absolute weight, anchored at the current (initial) material:

w_abs = w_frac · D0 / S0

where S0 = ½‖M·K(θ₀)·û_meas‖² is the raw equilibrium penalty (computed with a temporary weight-1 instance bound to the global ctx) and D0 is the data misfit of the forward solution at θ₀ (one global forward solve — the only nontrivial cost, comparable to a single global iteration's solve; pass a preallocated global stiffness via K to avoid an extra allocation).

On macOS (default), D0 is computed via a condensed, factorization-free iterative solve (_fb_d0_condensed_iterative) to avoid the xzone-malloc crash; the K preallocation is used only by the direct path (Linux, or macOS with SENTINEL_ALLOW_MACOS_FB_RELATIVE=1).

Returns the scale D0/S0, or nothing if no unresolved relative instance is present (in which case nothing is computed). Called by run_inverse_solve! before the first zone solve; safe to call directly for custom drivers.

MGH calibration (2026-07-12, 12-iter batched-GPU sweep): fractions ≤ 0.2 are gentle (data misfit slightly better than baseline, per-voxel corr ≥ 0.99 vs unregularized), 1.0 gave the best misfit with substantive reconstruction changes, and 5.0 over-regularized. See benchmark/results/mgh_fb_calibration/.

source

Sentinel.fb_begin_line_search! Function

fb_begin_line_search!(regs, material, delta_theta, grad) -> nothing

Arm the line-search caches on every bound ForceBalanceReg in regs, for a search from material (θ0) along delta_theta (the direction apply_material_update! applies). Both FB terms are exact in the step α because the Model-1 element stiffness is affine in μ (K_uu = μ·G − ω²ρ·M; K_up, K_pp are μ-independent):

• penalty: ½w(c0 + 2α·c1 + α²·c2) from residuals at θ0 and θ0 + 1·d;

• gradient: g(α) = g0 + α·g_d, where g0 reuses the CG-point gradient (add_gradient! remembers its last full evaluation keyed by the μ values, and the CG driver always evaluates at the search base), so arming costs the two residual loops plus a single extra gradient evaluation.

Trial-point compute_penalty/add_gradient! then cost O(nμ)/O(nparams). Only activates when the direction touches property 1 (μ) exclusively — a direction with κ content breaks affinity (K_pp ∝ 1/κ) and leaves the cache off (full evaluation per trial, exact as before). Trials knocked off the cached line (e.g. by bounds clamping) are detected by an α-projection consistency guard and evaluated fully. Call fb_end_line_search! after the search commits.

source

Sentinel.fb_end_line_search! Function

fb_end_line_search!(regs) -> nothing

Disarm the line-search quadratic cache on every ForceBalanceReg in regs (the buffers are kept for reuse by the next search).

source

CG Update Scaling

Sentinel.cg_update_scale Function

cg_update_scale(alpha, material_orig, material_new, threshold) -> Float64

Scale the CG step size alpha so that the relative material change does not exceed threshold.

Computes eps = material_epsilon(orig, new), then:

• If eps > threshold: return alpha * threshold / eps

• Otherwise: return alpha unchanged

source

Specification API

The high-level regularization spec types passed to InversionProblem (regularization = [...]). These are declarative descriptions that the solver lowers to the runtime regularizers documented above; they are distinct from the legacy *Reg classes.

Weights on TotalVariation, SpatialFilter, Tikhonov, Joachimowicz, SoftPrior, Marquardt, ExteriorConstraint, and CGResidualScaling may be per property: pass a vector aligned with props (e.g. weight = [1e-16, 2e-16]). Legacy .dat runfiles that set distinct weights per property now survive load_inversion and the TOML round trip without collapsing to the first property: per-property weights are preserved for every regularizer (no weight collapses). The one remaining scalar collapse is TotalVariation's delta (the Charbonnier smoothing epsilon), which the spec carries as a single value; a legacy runfile with per-property/complex delta collapses to the first value with a warning.

Sentinel.Tikhonov Type

Tikhonov(; props=:all, weight=1e-18)

Tikhonov (L2) penalty toward the initial property values (indicator 3).

source

Sentinel.TotalVariation Type

TotalVariation(; props=:all, weight=5e-16, delta=1e-19)

Total-variation penalty (runfile indicator 1).

source

Sentinel.SpatialFilter Type

SpatialFilter(; props=:all, sigma=0.003 => 0.0015)

Gaussian spatial filtering of the property update (post-processing smoothing, runfile indicator 2). sigma is the Gaussian width in meters, not a penalty weight.

source

Sentinel.Marquardt Type

Marquardt(; props=:all, start_weight=100.0, min_weight=1e-11,
            delta=0.5, adjust_threshold=0.25)

Marquardt adaptive damping (indicator 4).

source

Sentinel.Joachimowicz Type

Joachimowicz(; props=:all, weight=5.0)

Joachimowicz parameter-magnitude penalty (indicator 5).

source

Sentinel.ExteriorConstraint Type

ExteriorConstraint(; props=:all, weight=1e-14, bounds=nothing)

Soft penalty for property values outside bounds (indicator 6). bounds is one (min_re, max_re, min_im, max_im) tuple per property; nothing uses the standard EMPIRE bounds (3-property layout only).

source

Sentinel.CGResidualScaling Type

CGResidualScaling(; weight=0.04 => 0.01)

Conjugate-gradient residual scaling (indicator 7). Applies to all properties.

source

Sentinel.VanHouten Type

VanHouten(level=1.2)

Van Houten update-limiting post-processing (indicator 8).

source

Sentinel.McGarryBounds Type

McGarryBounds(; bounds=nothing)

Hard clamping of property values to bounds after each iteration (indicator 9). Bound semantics as in ExteriorConstraint; when both are present they must agree (the runfile format shares one bounds table).

source

Sentinel.SoftPrior Type

SoftPrior(; props=:all, weight=1e-11, delay=0)

Region-based soft-prior penalty (indicator 10).

Warning

Not yet wired into the solver: nothing reads the region stack into a RegionInfo, so InversionProblem warns and the term is ignored at setup. The type exists so the TOML schema is stable for when the loader lands.

source

Sentinel.ForceBalance Type

ForceBalance(; weight=0.1, relative=true)

Force-balance (equilibrium residual on the measured displacement field) penalty — a physics-consistency term for Model-1 (isotropic incompressible) reconstructions. It has no Fortran indicator slot; it is a Julia-side extension carried alongside the runfile regularizers.

In relative mode (the default) weight is a fraction of the initial data misfit D0/S0, resolved to an absolute weight once at solve start — the same number is portable across meshes and datasets. With relative=false it is an absolute weight (mesh- and dataset-dependent). weight may be a scalar or an init => final ramp.

Calibrated guidance (benchmark/results/mgh_fb_calibration/, benchmark/results/fb_phantom_arbitration/): weight ≤ 0.2 is gentle and safe on noisy in-vivo data; larger fractions improve accuracy only on clean/high-SNR/denoised data.

source

Sentinel.default_regularization Function

default_regularization() -> Vector{RegularizationSpec}

The regularization stack InversionProblem applies when the regularization keyword is not given:

RegularizationSpec[SpatialFilter(), CGResidualScaling(), VanHouten(), McGarryBounds()]

Gaussian spatial filtering of the property update, conjugate-gradient residual scaling, Van Houten update limiting, and McGarry hard bounds — each with its constructor defaults. Returns a fresh vector on every call, so the result can be modified and passed back via the regularization keyword.

source