◆ FORGE Suite
GitHubMechanical Neuroimaging Lab · Univ. of Delaware
Skip to content

Optimization ​

Three optimization algorithms are available for iterative reconstruction. All use Armijo backtracking line search.

See the Mathematical Reference for update rule equations.

Problem Context ​

Sentinel.ForwardProblemContext Type
julia
ForwardProblemContext

Bundles all immutable problem data needed to evaluate the objective function and gradient. Passed to line search and optimizer routines to avoid long argument lists.

Fields

  • grid: Ferrite Grid

  • dh: DofHandler

  • cv_disp: CellValues for displacement interpolation

  • cv_press: CellValues for pressure (or nothing)

  • model: AbstractMaterialModel singleton

  • gp2mtrs: vector of GP2Mtr per material mesh

  • meshes: vector of MaterialMesh

  • bcs: BoundaryConditions

  • meas: measured Displacement

  • ω: angular frequency rad/s

  • ρ: density kg/m³

  • regularizations: vector of AbstractRegularization (or nothing)

  • solver: linear solver backend

  • numdispsets: number of displacement/excitation sets

  • fiber_field: per-element fiber direction vectors (or nothing)

source
Sentinel.ConvergenceHistory Type
julia
ConvergenceHistory

Records iteration-by-iteration progress of an optimization run.

Fields

  • iterations: number of iterations completed

  • objective: objective function values per iteration

  • disp_error: displacement-only error per iteration

  • mat_epsilon: relative material change per iteration

  • converged: whether the optimizer declared convergence

source

Helper Functions ​

Sentinel.copy_material Function
julia
copy_material(mat::Material) -> Material

Deep-copy a Material struct, duplicating all SingleProperty arrays.

source
Sentinel.apply_material_update! Function
julia
apply_material_update!(material::Material, delta_theta::AbstractVector,
                        alpha::Float64, grad::MaterialGradient)

Update material property values: θ += α·Δθ/scalar, using the gradient's realparam2prop reverse mapping to write each element of delta_theta back to the correct rvalue or ivalue entry.

The division by scalar matches Fortran's incrementmaterial (FEmaterial.f90): realinc = direction%rvalue(direcind) / origmaterial%prop(ii)%scalar(iii) This converts the step from physical units (produced by CG gradient scaling) back to normalized parameter units.

source
Sentinel.evaluate_objective! Function
julia
evaluate_objective!(calc::Displacement, K::SparseMatrixCSC{ComplexF64},
                     f_work::Vector{ComplexF64},
                     ctx::ForwardProblemContext, material::Material;
                     iteration::Int=1) -> Float64

Run the forward solve for all displacement sets and compute the objective.

  1. Assemble stiffness from current material

  2. For each dispset: apply BCs, solve, store in calc

  3. Return compute_objective(calc, meas, material, regularizations)

source
Sentinel.compute_full_gradient! Function
julia
compute_full_gradient!(grad::MaterialGradient, calc::Displacement,
                        adjnt::Displacement,
                        K::SparseMatrixCSC{ComplexF64},
                        ctx::ForwardProblemContext, material::Material;
                        iteration::Int=1)

Compute the full gradient including adjoint solves and regularization.

  1. For each dispset: adjoint_solve!

  2. zero!(grad), compute_gradient! (adjoint contribution)

  3. add_gradient! for each regularization term

source
Sentinel.apply_postprocessing! Function
julia
apply_postprocessing!(material::Material, regularizations)

Apply all post-processing regularizations (Van Houten, Spatial Filter, McGarry bounds) in sequence.

source

Conjugate Gradient ​

Polak-Ribiere+ nonlinear CG:

dk=−gk+βkdk−1,βk=max(0,gkT(gk−gk−1)gk−1Tgk−1)
Sentinel.conjugate_gradient! Function
julia
conjugate_gradient!(material::Material, ctx::ForwardProblemContext;
                     max_iter::Int=50, tol::Float64=1e-6,
                     alpha_init::Float64=1.0,
                     cg_threshold::Float64=0.0,
                     cg_scalar_weights::Union{Matrix{Float64}, Nothing}=nothing,
                     iteration_offset::Int=0) -> (Material, ConvergenceHistory)

Run nonlinear conjugate gradient iterations to minimize the objective.

Arguments

  • material: initial material estimate (modified in-place on convergence)

  • ctx: forward problem context

  • max_iter: maximum number of CG iterations

  • tol: convergence tolerance on material_epsilon

  • alpha_init: initial line search step size

  • cg_threshold: if > 0, cap relative material change per step via cg_update_scale

  • cg_scalar_weights: (numprop, 2) matrix of per-property CG scalar weights for gradient direction scaling. If provided, scales the gradient so the average step for each property equals weight * avg_property_value. Ports regind(7) from Fortran.

  • iteration_offset: offset for regularization weight schedule evaluation

  • step_method: :line_search (default) or :barzilai_borwein

  • bb_window: non-monotone acceptance window size (default 5, BB only)

  • bb_max_backtrack: max backtracking halvings before line search fallback (default 3, BB only)

Returns

  • (material, history, calc): updated material, convergence history, and displacement
source

Gauss-Newton ​

Diagonal Hessian preconditioned descent: Δθk=−gk/Hkk

Sentinel.gauss_newton! Function
julia
gauss_newton!(material::Material, ctx::ForwardProblemContext;
               max_iter::Int=20, tol::Float64=1e-6,
               alpha_init::Float64=1.0,
               hessian_floor::Float64=1e-30,
               iteration_offset::Int=0) -> (Material, ConvergenceHistory)

Run Gauss-Newton iterations to minimize the objective.

Arguments

  • material: initial material estimate (modified in-place)

  • ctx: forward problem context (regularizations should include Hessian-contributing regs)

  • max_iter: maximum iterations

  • tol: convergence tolerance on material_epsilon

  • alpha_init: initial line search step size

  • hessian_floor: minimum diagonal value to prevent division by zero

  • iteration_offset: offset for regularization weight schedules

Returns

  • (material, history): updated material and convergence history
source

L-BFGS ​

Limited-memory BFGS with two-loop recursion and initial scaling γ=sTy/(yTy).

Sentinel.quasi_newton! Function
julia
quasi_newton!(material::Material, ctx::ForwardProblemContext;
               max_iter::Int=50, tol::Float64=1e-6,
               memory::Int=5, alpha_init::Float64=1.0,
               iteration_offset::Int=0) -> (Material, ConvergenceHistory)

Run L-BFGS quasi-Newton iterations to minimize the objective.

Arguments

  • material: initial material estimate (modified in-place)

  • ctx: forward problem context

  • max_iter: maximum iterations

  • tol: convergence tolerance on material_epsilon

  • memory: number of L-BFGS curvature pairs to store (default 5)

  • alpha_init: initial line search step size

  • iteration_offset: offset for regularization weight schedules

Returns

  • (material, history): updated material and convergence history
source

Armijo backtracking: reduce α by factor ρ until sufficient decrease Φ(θ+αd)≤Φ(θ)+c1αgTd.

Sentinel.line_search Function
julia
line_search(ctx::ForwardProblemContext, material::Material,
            grad::MaterialGradient, direction::AbstractVector{Float64},
            obj_current::Float64;
            alpha_init::Float64=1.0, c1::Float64=1e-4,
            rho::Float64=0.5, max_iter::Int=20,
            iteration::Int=1) -> (Float64, Float64, Material)

Backtracking line search with Armijo sufficient decrease condition.

Arguments

  • ctx: forward problem context

  • material: current material state

  • grad: current gradient (must be initialized, with rvalue populated)

  • direction: search direction vector (same length as grad.rvalue)

  • obj_current: current objective value Φ(θ_k)

  • alpha_init: initial step size (default 1.0)

  • c1: sufficient decrease parameter (default 1e-4)

  • rho: backtracking reduction factor (default 0.5)

  • max_iter: maximum number of backtracking steps (default 20)

  • iteration: current optimizer iteration (for regularization weight schedules)

Returns

  • (alpha, obj_new, mat_new): step size, new objective value, updated material If no decrease is found, returns (0.0, obj_current, copy_material(material)).
source