Imaging (BSDMM)

reconstruct_bsdmm(x_init, data, ft; kwargs...)

Block-Structured SDMM (BSDMM) image reconstruction for optical interferometry.

Solves: minimize_x f(x) + Σ_i g_i(x)

where f(x) = χ²(H(x/sum(x))) is the non-convex data fidelity term (squared visibilities, closure phases) and each g_i is a non-smooth regularizer (total variation, L2 smoothness, centering) handled via its proximal operator.

The BSDMM splits the problem into alternating steps:

1. z-updates (proximal): z_i = prox_{g_i/ρ_i}(x + u_i) — one per regularizer
2. x-update (VMLMB): minimize χ² + Σ_i (ρ_i/2)||x - z_i + u_i||² with x ≥ 0
3. Dual updates: u_i ← u_i + x - z_i

Positivity is enforced via VMLMB lower bounds. Flux normalization uses the chain-rule gradient ∇_x χ²(H(x/sum(x))) inside the x-update, plus explicit re-normalization after each step.

Supports monochromatic (nx × nx image, (1,1) data) and polychromatic (nx × nx × nwav × 1 image, (nwav,1) data) modes.

Arguments

• x_init: starting image (nx × nx for mono, nx × nx × nwav × 1 for poly)
• data: interferometric data (Matrix{OIdata})
• ft: Fourier transform plans from setup_ft

Keyword arguments (monochromatic)

• mu_reg=0.0: regularization weight (0 = no regularization)
• mu_cen=0.0: centering weight (0 = no centering)
• reg_type=:tv: regularization type (:tv or :l2smooth)

Keyword arguments (polychromatic)

• mu_tv=0.0: per-channel spatial TV/L2smooth weight
• mu_group=0.0: transspectral coupling weight (group sparsity or group TV)
• mu_cen=0.0: per-channel centering weight
• reg_type=:tv: spatial regularizer (:tv or :l2smooth)
• group_type=:sparsity: transspectral regularizer (:sparsity or :tv)

Common keyword arguments

• tv_niter=50: Chambolle iterations for TV proximal operator
• rho_init=10.0: initial ADMM penalty parameter for all blocks
• rho_min=1e-3: minimum penalty (for adaptive update)
• rho_max=1e4: maximum penalty (for adaptive update)
• adaptive=false: adaptive penalty strategy:
  • false: fixed ρ throughout
  • true or :spectral: AADMM spectral (adaptive BB + predictive dual + cross-block balancing)
  • :balanced: residual balancing (adjusts ρ based on primal/dual residual ratio)
  • :aradmm: ARADMM (spectral penalty + over-relaxation parameter γ)
• adp_freq=2: adaptive update frequency (every N outer iterations)
• adp_start=1: first iteration for adaptive updates
• maxit=300: maximum BSDMM outer iterations
• x_maxiter=50: maximum VMLMB iterations per x-update (use small values, e.g. 5)
• verb=true: print convergence info every 10 iterations
• history=false: when true, returns (image, hist) where hist is a NamedTuple with fields chi2, obj, r_norm, s_norm, rho (Dict), gamma

Returns

• Image with best objective value (same shape as x_init)
• If history=true: (image, hist) tuple

References

• Moolekamp & Melchior (arXiv:1708.09066) — BSDMM algorithm
• Xu, Taylor & Goldstein — ARADMM / AMADMM adaptive penalty
• Chambolle (2004) — TV proximal operator

See also: reconstruct, prox_tv, prox_l2smooth, prox_group_sparsity, prox_grouptv

prox_tv(y, lambda; niter=100, tau=1/8)

Proximal operator for isotropic total variation (Chambolle 2004). Solves argmin_x (1/2)||x - y||² + λ TV(x) via dual projected gradient ascent.

prox_l2smooth(y, lambda)

Proximal operator for L2 smoothness via FFT. Solves argmin_x (1/2)||x - y||² + (λ/2)||∇x||² using the discrete Laplacian eigendecomposition with periodic boundary conditions.

prox_centering!(x, y, lambda, p, q, AtA)

Proximal operator for centering: argmin_x (1/2)||x-y||² + λ[(p'x)² + (q'x)²]. Uses Woodbury formula with p = centered column indices, q = centered row indices, AtA = pre-computed [p'p p'q; q'p q'q] (2×2 matrix). Result stored in x.

tv_norm(x_img)

Isotropic total variation of a 2D image.

l2smooth_norm(x_img)

L2 smoothness norm (half squared gradient norm) of a 2D image.

prox_group_sparsity(y, lambda)

Proximal operator for group sparsity (L2,1 mixed norm) on a 3D cube (nx, ny, nwav). Block soft-thresholding: at each spatial pixel, the across-channel vector is shrunk toward zero, promoting shared support across wavelength channels.

Solves argmin_x (1/2)||x - y||² + λ Σ_{i,j} ||x[i,j,:]||₂.

prox_grouptv(y, lambda; niter=100, tau=1/8)

Proximal operator for group (vectorial) total variation on a 3D cube (nx, ny, nwav). Encourages shared edge locations across all wavelength channels via Chambolle dual projection with coupled L2,1 norm across the spectral dimension.

Solves argmin_x (1/2)||x - y||² + λ GroupTV(x).

group_sparsity_norm(x_cube)

L2,1 mixed norm of a 3D image cube: Σ_{i,j} ||x[i,j,:]||₂.

grouptv_norm(x_cube)

Group (vectorial) total variation of a 3D image cube: Σ_{i,j} √( Σ_λ [|∇_h x_{i,j,λ}|² + |∇_v x_{i,j,λ}|²] ).