Mathematical Formulation#
This section collects the common variational, iterative, entropy-based,
Bayesian, parametric, and first-order-optimisation formulations used
throughout bssunfold.
The package currently ships 93 public unfolding methods in
bssunfold.core, which this page groups into eleven method families.
For each family we state the optimisation objective, the constraint
structure, the convergence facts that are actually guaranteed by the
literature, and the stopping criteria used in practice. Terms are
defined where they first appear; every quantitative claim is backed by a
reference with a DOI in References.
General Inverse Problem#
The neutron spectrum unfolding problem is modelled as a discretised Fredholm integral equation of the first kind:
Here \(\mathbf{N}\in\mathbb{R}^M\) denotes the detector readings,
\(\mathbf{R}\in\mathbb{R}^{M\times K}\) the response matrix, and
\(\boldsymbol{\Phi}\in\mathbb{R}^K_{\ge 0}\) the non-negative
spectrum to be reconstructed. The reference RF_GSF detector of this
package uses \(K = 60\) logarithmically spaced energy bins spanning
1e-9 to 631 MeV and \(M = 10\) Bonner spheres; multisphere systems in
general use between roughly \(M = 8\) and \(M = 18\) spheres.
Because \(M \ll K\), the discretised problem is severely
underdetermined (there exist infinitely many non-negative spectra that
reproduce the readings exactly), and because the singular values of
\(\mathbf{R}\) decay exponentially fast (Hansen, 1994), it is
ill-posed: small perturbations of \(\mathbf{N}\) can produce
arbitrarily large perturbations of \(\boldsymbol{\Phi}\). Every
practical solver therefore combines one of the formulations below with
regularisation, positivity constraints, or early stopping.
Terms used throughout this page.
Response matrix \(\mathbf{R}\) — the \(M\times K\) table of detector sensitivities: entry \(\mathbf{R}_{ij}\) is the count rate that sphere \(i\) registers per unit fluence in energy bin \(j\) (IAEA, 2001).
Regularisation parameter \(\alpha\) — the weight that trades data fit against smoothness (or another prior belief). Too small: the solver fits noise (over-fitting, oscillating spectrum). Too large: the solver ignores the data (over-smoothing, biased fluence and dose). Selection rules are discussed in Parameter Selection and Stopping Criteria.
Condition number \(\kappa(\mathbf{R}) = \sigma_{\max}/\sigma_{\min}\) — the ratio of extreme singular values; it lower-bounds the relative error amplification \(\|\delta\boldsymbol{\Phi}\|/\|\boldsymbol{\Phi}\| \lesssim \kappa(\mathbf{R})\,\|\delta\mathbf{N}\|/\|\mathbf{N}\|\) (Tikhonov et al., 1995). Bonner-sphere matrices are ill-conditioned precisely because many columns (energy responses of adjacent sphere sizes) are nearly collinear.
Prior spectrum \(\boldsymbol{\Phi}_0\) — a physically plausible initial guess (flat, Watt fission, or a catalogue shape) used by ratio, entropy, and Bayesian methods. A good prior speeds up convergence and stabilises ill-conditioned directions, but a wrong prior biases the result towards itself — the classic bias-variance trade-off of regularisation (Reginatto, 2010; Engl, Hanke and Neubauer, 1996).
Noise model. Counting statistics make the readings Poisson distributed,
Most solvers in bssunfold work with the standard Gaussian
approximation of this model, \(\boldsymbol{\varepsilon}\sim
\mathcal{N}(\mathbf{0}, \boldsymbol{\Sigma})\) with diagonal
\(\boldsymbol{\Sigma}\), where \(\sigma_i^2 \approx N_i\) for
pure counting statistics; an additional relative-uncertainty floor (a
systematic component, e.g. sigma_factor in unfold_maxed) is
added in practice to account for calibration and response-matrix
uncertainties. The Gaussian likelihood underlies all
\(\chi^2\)-based objectives in this page, while the Poisson
likelihood is used explicitly by the EM family (unfold_mlem and
relatives).
Semi-convergence. Iterative methods applied to this problem share one important property: after an initial phase in which the smooth, data-dominated components of the solution are recovered, further iterations start amplifying noise associated with the smallest singular values of \(\mathbf{R}\). This semi-convergence behaviour is not restricted to Krylov methods — it is characteristic of Landweber and CGLS iterations, ART/SART, MLEM and its Bayesian (D’Agostini) relative alike (Engl, Hanke and Neubauer, 1996; Hansen, 1994). Consequently the iteration number itself acts as a regularisation parameter, and every iterative family in this package offers or recommends an early-stopping rule, most commonly the discrepancy principle or GCV (see Parameter selection).
Parameter Selection and Stopping Criteria#
Several families share the same parameter-selection machinery, so the definitions are collected once here.
Discrepancy principle (Morozov). Let \(\delta\) be a bound on the data noise, \(\|\boldsymbol{\varepsilon}\|_2 \le \delta\) (e.g. \(\delta^2 = \sum_i \sigma_i^2\)). The regularisation parameter \(\alpha\) (or the stopping index \(k\)) is chosen as the smallest value for which the residual satisfies
with a safety factor \(\tau > 1\), typically \(\tau \in
[1.01,\,1.1]\) (Morozov, 1984; Engl, Hanke and Neubauer, 1996). When the
\(M\) residuals are scaled by their known standard deviations
\(\sigma_i\), the reduced criterion is
\(\chi^2 \approx M\): for \(M\) independent bins with correctly
known variances, \(\chi^2\) is asymptotically
\(\chi^2_M\)-distributed, with expectation \(M\) and variance
\(2M\), so values outside roughly \(M \pm \sqrt{2M}\) indicate
over-fitting (chi-square too small, noise fitted) or over-regularisation
(chi-square too large). This \(M \pm \sqrt{2M}\) band is used by
several solvers (e.g. the automatic smoothing-weight bisection in
unfold_ferdor and the entropy-weight calibration of unfold_maxed
derivatives) as a rule of thumb, not as a rigorous test.
Generalised cross-validation (GCV). For a filtering or iterative method whose effective influence matrix is \(\mathbf{A}_k\), choose \(k\) minimising
which requires no noise-level estimate (Golub, Heath and Wahba, 1979).
L-curve. For a family of regularised solutions \(\boldsymbol{\Phi}_\alpha\), plot \(\log\|\mathbf{L}\boldsymbol{\Phi}_\alpha\|_2\) against \(\log\|\mathbf{R}\boldsymbol{\Phi}_\alpha - \mathbf{N}\|_2\) and select the point of maximum curvature (the corner) as the best compromise between residual size and solution norm (Hansen, 1992).
KKT residuals. For the QP-based solvers (unfold_cvxpy,
unfold_qpsolvers, unfold_interpret), the Karush-Kuhn-Tucker
residuals reported at termination measure first-order optimality of
the numerical solver — they certify that the convex program has been
solved to the requested tolerance. They are solver diagnostics, not a
statement about statistical convergence of the unfolding itself; the
statistical quality of the unfolded spectrum is still governed by the
regularisation choice and the noise model.
Variational (Tikhonov-Type) Methods#
Methods: unfold_cvxpy, unfold_qpsolvers, unfold_tikhonov_tv,
unfold_tikhonov_legendre, unfold_statreg, unfold_reconst,
unfold_epic, unfold_scipy_direct_method, unfold_cs,
unfold_nnksvd, unfold_docplex, unfold_scip, unfold_smt,
unfold_qubo, unfold_interpret.
The generic penalised objective is
where \(\mathbf{W}\) encodes the measurement covariances, \(\mathbf{L}_0 = \mathbf{I}\), \(\mathbf{L}_1\) (first difference) or \(\mathbf{L}_2\) (second difference / Legendre-basis derivatives) selects the roughness penalty, and \(p \in \{1, 2\}\). For \(p = 2\) the unconstrained objective is quadratic; if \(\ker(\mathbf{R}) \cap \ker(\mathbf{L}) = \{\mathbf{0}\}\), the Hessian \(\mathbf{R}^T\mathbf{W}^T\mathbf{W}\mathbf{R} + \alpha\,\mathbf{L}^T\mathbf{L}\) is positive definite and the minimiser is unique (Tikhonov et al., 1995; Engl, Hanke and Neubauer, 1996). This is the typical situation for \(\mathbf{L} = \mathbf{L}_1\) or \(\mathbf{L}_2\) on the log-spaced grid of a Bonner sphere system.
For total-variation regularisation the penalty becomes
TV regularisation is convex but not strictly convex: for
\(p = 1\) uniqueness of the minimiser is not guaranteed, and
solutions may be non-unique along flat segments. The statement
“unique for \(\alpha > 0\)” therefore applies only to the
\(p = 2\) (quadratic) case above. unfold_tikhonov_tv uses the
standard smoothed-TV approximation to retain differentiability.
Closed-form statistical regularisation. unfold_statreg
(Turchin’s method, cf. the RECONST algorithm ported by
unfold_reconst) solves the linear system directly:
implemented as solve(RᵀΣ⁻¹R + αΩ, RᵀΣ⁻¹N) with a fallback to
lstsq. This closed form is exact only for the unconstrained
problem: the moment a non-negativity constraint
\(\boldsymbol{\Phi}\ge 0\) is enforced, the problem becomes a
quadratic program and the closed-form solution no longer applies (it
is then merely the solution of the equality-constrained relaxation).
Sparse formulations. unfold_cs (compressed sensing with OMP /
K-SVD style dictionaries), unfold_nnksvd (non-negative K-SVD
dictionary learning) and the QUBO/Ising formulation unfold_qubo
express the solution as sparse over-complete expansions; they are
solved with the corresponding combinatorial or convex optimisation
backends (unfold_docplex, unfold_scip, unfold_smt).
Convergence guarantees in this family are those of the underlying
numerical backend, plus the general regularisation theory above.
Spectral Filtering and Krylov Methods#
Methods: unfold_tsvd, unfold_lanczos, unfold_cgls,
unfold_gks, unfold_hybrid_gmres.
The TSVD solution is defined by the truncated singular value expansion
where \(\mathbf{u}_i\), \(\mathbf{v}_i\) and \(\sigma_i\)
are the singular vectors and values of \(\mathbf{R}\). The
truncation index \(k\) is chosen by GCV, the discrepancy principle
or the L-curve (see Parameter Selection and Stopping Criteria). In the implementation of
unfold_tsvd the small singular-value contributions are discarded
and any remaining negative bins are clipped,
\(\boldsymbol{\Phi}\leftarrow\max(0, \boldsymbol{\Phi})\). This
post-hoc clipping is a heuristic projection applied after the
filtering step; it is not equivalent to solving the non-negativity-
constrained filtering problem exactly.
CGLS (unfold_cgls, Hestenes and Stiefel, 1952) and LSQR-type
recurrences (Paige and Saunders, 1982) applied to the regularised
least-squares problem converge towards the (Tikhonov-regularised)
solution while implicitly filtering the singular spectrum; the
Krylov-subspace methods unfold_gks and unfold_hybrid_gmres
build the Golub-Kahan bidiagonalisation
(Golub and Kahan, 1965) and solve the projected problem with an explicit Tikhonov regularisation term whose parameter is re-selected at each iteration — the hybrid strategy (Gazzola, Hansen and Nagy, 2018). Hybrid methods delay but do not eliminate semi-convergence; the inner regularisation parameter is chosen by GCV or discrepancy rules. All methods in this section exhibit the semi-convergence behaviour described in General Inverse Problem and must be stopped early.
Algebraic Iterative Methods#
Methods: unfold_landweber, unfold_kaczmarz,
unfold_randomized_kaczmarz, unfold_sart, unfold_fista,
unfold_doroshenko.
The Landweber iteration with positivity projection is
which converges for \(0 < \omega < 2/\sigma_{\max}^2\) where
\(\sigma_{\max}\) is the largest singular value of
\(\mathbf{R}\) (unfold_landweber uses the conservative step
\(\omega = 1/\sigma_{\max}^2\)). Kaczmarz/ART-type updates
(unfold_kaczmarz) project onto the hyperplane defined by one
detector row at a time; the randomized Kaczmarz variant
(unfold_randomized_kaczmarz) selects rows with probability
proportional to \(\|\mathbf{R}_i\|_2^2\) and converges linearly in
expectation at the rate
\(\mathbb{E}\,\|\boldsymbol{\Phi}^{(k)} -
\boldsymbol{\Phi}_\star\|_2^2 \le (1 - \sigma_{\min}^2 /
\|\mathbf{R}\|_F^2)^k\,\|\boldsymbol{\Phi}^{(0)} -
\boldsymbol{\Phi}_\star\|_2^2\) (Strohmer and Vershynin, 2009).
For SART the implementation in unfold_sart follows the classical
row-and-column-normalised scheme (Andersen and Kak, 1984):
i.e. the residual is normalised by the forward-projected unit spectrum and the update by the back-projected unit sensitivity, with relaxation \(\lambda(n)\) (constant or sequence; convergence of the exact SART recursion is guaranteed for \(0 < \lambda \le 1\), and more generally the algebraic family requires \(0 < \lambda < 2\); Jiang and Wang, 2003). Note that this is the simultaneous scheme — all rows contribute in each sweep — which is distinct from the purely row-by-row ART/Kaczmarz recursion.
For FISTA the accelerated proximal step reads
with \(f\) the smooth least-squares term, \(h\) the
non-smooth penalty (non-negativity or \(\ell_1\)) and the step
\(\gamma \le 1/\|\mathbf{R}^T\mathbf{R}\|_2 = 1/\sigma_{\max}^2\).
FISTA improves the worst-case objective-value rate from
\(\mathcal{O}(1/k)\) to \(\mathcal{O}(1/k^2)\) (Beck and
Teboulle, 2009); the iterate sequence itself is non-monotone, so the
objective may temporarily increase — monotone behaviour requires
modified variants (e.g. MFISTA). unfold_fista follows the IRtools
implementation (Gazzola, Hansen and Nagy, 2018).
unfold_doroshenko applies a related coordinate-update relaxation of
the same least-squares objective.
All methods in this section are stopped by semi-convergence-aware rules: a relative-change threshold, a discrepancy-principle test, or a maximum iteration count.
First-Order Optimisation Methods (MIPT Course Port)#
Methods: unfold_pgd, unfold_frank_wolfe, unfold_mirror_descent,
unfold_admm, unfold_lbfgsb, unfold_coordinate_descent,
unfold_subgradient, unfold_extragradient,
select_regularization_1d.
This family ports the algorithmic core of the MIPT course “Optimization Methods in Machine Learning” (lectures 7-15; the course repository contains lecture notes and homework notebooks). Each method is a classical first-order scheme applied to the unfolding objective; the pedagogical value is that every design choice has a provable effect, stated below.
Projected gradient descent (PGD) — unfold_pgd (lecture 9 /
homework 14). The plain gradient step
\(\boldsymbol{\Phi}^{(k+1)} = \boldsymbol{\Phi}^{(k)} -
\gamma\nabla f(\boldsymbol{\Phi}^{(k)})\) is followed by the Euclidean
projection \(\mathcal{P}_{\mathcal{C}}\) onto the constraint set
\(\mathcal{C}\) — the non-negative orthant, a box, or the
fluence simplex \(\{\boldsymbol{\Phi}\ge 0,\ \sum_j \Phi_j = F\}\).
For a fixed total fluence the simplex projection preserves the physical
normalisation exactly at every iterate — unlike iterative solvers,
whose fluence drifts during semi-convergence. Convergence for smooth
convex \(f\) is
\(f(\boldsymbol{\Phi}^{(k)}) - f^\star = \mathcal{O}(1/k)\) for
\(\gamma \le 1/L\) with \(L = \|\mathbf{R}^T\mathbf{R}\|_2\)
(Beck and Teboulle, 2009); Armijo backtracking (halving
\(\gamma\) until the Armijo descent condition holds) makes the step
adaptive to the local curvature and removes the need to estimate
\(L\) manually (Nocedal and Wright, 2006). The reported
NNLS duality gap — the difference between the primal objective and a
dual-feasible lower bound reconstructed from multipliers — is an
optimality certificate: when it reaches zero the returned spectrum is
provably the constrained optimum, not merely a stationary point
(Nocedal and Wright, 2006; Boyd and Vandenberghe, 2004).
Frank-Wolfe (conditional gradient) — unfold_frank_wolfe
(lecture 9). Instead of a gradient step plus projection, the method
minimises the linearisation of \(f\) over the fluence simplex —
a trivial linear minimisation oracle that picks the single most
descent-deserving energy vertex — and moves a fraction
\(\gamma \in (0,1]\) towards it. The Frank-Wolfe gap
\(g(\boldsymbol{\Phi}^{(k)}) =
\nabla f(\boldsymbol{\Phi}^{(k)})^T(\boldsymbol{\Phi}^{(k)} -
\mathbf{s}^{(k)})\) with \(\mathbf{s}^{(k)}\) the oracle vertex is
always an upper bound on the sub-optimality, so it doubles as a
natural stopping certificate (Frank and Wolfe, 1956;
Jaggi, 2013). Wolfe’s away steps let iterates move away from
previously chosen vertices, curing the classical zig-zagging of
boundary solutions (Wolfe, 1970). Because every iterate is a convex
combination of simplex vertices, the total fluence is preserved to
machine precision — the method of choice when the fluence
normalisation is itself a measured quantity.
Mirror descent — unfold_mirror_descent (lecture 10 / homework
16). Generalises projected gradient descent by replacing the
Euclidean geometry with a Bregman divergence generated by a
mirror map \(\psi\) (Nemirovski and Yudin, 1983;
Beck and Teboulle, 2003). The entropy map
\(\psi(\boldsymbol{\Phi}) = \sum_j \Phi_j \ln \Phi_j\) yields the
multiplicative update
\(\Phi_j^{(k+1)} \propto \Phi_j^{(k)}
\exp(-\gamma \nabla_j f)\) — the same recursion family as MLEM, GRAVEL
and SAND-II, derived here from first principles. Multiplicative
updates preserve positivity and act relatively (a bin with large
fluence is corrected proportionally more), which is why
entropy-geometry solvers dominate Poisson-count unfolding. The
\(\ell_2\) and p-norm maps give additive updates with different
bias profiles: the \(\ell_2\) geometry distributes corrections
uniformly and is more robust to near-zero bins that dominate the
log-space relative error. A per-iteration golden-section line search
along the mirror trajectory keeps the step inside the orthant interior
(Nemirovski and Yudin, 1983).
Consensus ADMM — unfold_admm (lecture 11 / homework 18). The
splitting
\(\min f(\boldsymbol{\Phi}) + g(\mathbf{z})\) s.t.
\(\boldsymbol{\Phi} = \mathbf{z}\) turns the L1/TV-regularised
problem into alternating (i) an exact non-negative least squares
x-update on the augmented system, (ii) a soft-thresholding z-update
that realises the proximal operator of the \(\ell_1\) penalty, and
(iii) a scaled dual ascent step (Boyd et al., 2011; Gabay and Mercier,
1976). The primal/dual residual stopping test
\(\|\mathbf{R}\boldsymbol{\Phi}^{(k)} - \mathbf{z}^{(k)}\| \to 0\),
\(\rho(\mathbf{z}^{(k+1)}-\mathbf{z}^{(k)}) \to 0\) certifies
consensus between the data-fidelity and penalty variables. The
adaptive penalty \(\rho\) (Boyd et al., 2011, sec. 3.4.1)
balances the two residual norms each outer iteration; it makes the
solver insensitive to the absolute count scale of the problem — the
single most common failure mode of fixed-\(\rho\) ADMM on
Poisson-counted data.
L-BFGS-B quasi-Newton — unfold_lbfgsb (lecture 7 / homework
10). A limited-memory quasi-Newton method with box bounds: curvature
pairs \((\mathbf{s}_i, \mathbf{y}_i)\) from the last
\(m\) iterates approximate the Hessian inverse, giving
superlinear-type local convergence on the smooth Tikhonov objective at
\(\mathcal{O}(n\,m)\) memory (Byrd et al., 1995; Nocedal and
Wright, 2006). The optional second-difference (curvature) penalty
\(\|\mathbf{D}_2\boldsymbol{\Phi}\|^2\) penalises oscillations;
increasing it visibly suppresses high-frequency noise but flattens
sharp spectral peaks — the parameter should be raised only until the
residual chi-square leaves the \(M \pm \sqrt{2M}\) band of
Parameter Selection and Stopping Criteria.
Coordinate descent — unfold_coordinate_descent (lecture 15).
Minimises the NNLS objective one coordinate at a time in closed form,
\(\Phi_j \leftarrow \max\!\bigl(0,\,
(\mathbf{a}_j^T\mathbf{r} + \|\mathbf{a}_j\|^2\Phi_j -
\lambda_1)/(\|\mathbf{a}_j\|^2 + \lambda_2)\bigr)\), with the residual
\(\mathbf{r}\) updated in \(\mathcal{O}(M)\) per coordinate
(Wright, 2015; Luo and Tseng, 1992). Each coordinate update
monotonically decreases the objective, which makes the method
extremely robust on ill-conditioned systems — at the cost of slower
progress on strongly coupled bins. Cyclic order is deterministic;
seeded random order escapes the worst-case cycling of badly scaled
columns (Wright, 2015).
Projected subgradient — unfold_subgradient (lecture 8 /
homework 12). For the nonsmooth L1/TV objective the gradient is
replaced by a subgradient — any vector in the subdifferential
(Nemirovski and Yudin, 1983; Shor, 1985). The step-size policies have
sharply different guarantees: the Polyak step
\(\gamma_k = (f(\boldsymbol{\Phi}^{(k)}) -
f^\star)/\|\mathbf{g}^{(k)}\|^2\) converges when the optimal value
\(f^\star\) is estimated well (Polyak, 1967), square-summable
diminishing steps (\(\sum_k \gamma_k = \infty\),
\(\sum_k \gamma_k^2 < \infty\)) converge for any fixed
\(f^\star\) but slowly, and fixed steps converge only to a
neighbourhood of the optimum whose radius is proportional to
\(\gamma\). The best-iterate-by-objective return is the standard
remedy for subgradient non-monotonicity.
Extragradient (Korpelevich) — unfold_extragradient
(lecture 13 / homework 20). The robust saddle formulation
\(\min_{\boldsymbol{\Phi}\ge 0}\ \max_{\|\mathbf{y}\|\le 1}
\frac{1}{2}\|\mathbf{R}\boldsymbol{\Phi}-\mathbf{N}\|^2 +
\delta\,\mathbf{y}^T(\mathbf{R}\boldsymbol{\Phi}-\mathbf{N})\) is
equivalent to least squares against an adversary that may add any
bounded perturbation of L2 norm \(\delta\) to the readings — the
uncertainty-aware reformulation of noise-robust unfolding. The plain
projected gradient method on saddle problems oscillates; Korpelevich’s
extra prediction-correction gradient step restores convergence at the
rate \(\mathcal{O}(1/k)\) for monotone operators
(Korpelevich, 1976; Facchinei and Pang, 2003). Larger
noise_level widens the guaranteed robustness radius at the price
of a more conservative (smoother) spectrum.
1D parameter search — select_regularization_1d (lecture 1 /
homework 1) sweeps \(\log_{10}\alpha\) with golden-section /
dichotomy / Brent minimisation of a GCV, Morozov or predictive-risk
score, replacing hand-tuned regularisation grids (Brent, 1973;
Golub, Heath and Wahba, 1979).
When to choose which. Use PGD/FISTA-type methods for the fastest smooth-objective descent; Frank-Wolfe when the fluence must be conserved exactly; mirror descent (entropy map) for Poisson-counted data; ADMM when an exact L1/TV penalty is required with certified consensus; L-BFGS-B for the largest smooth problems; coordinate descent for maximal per-iteration robustness; subgradient methods only when a nonsmooth objective cannot be split; extragradient when robustness to bounded reading errors must be explicit.
Iterative Ratio Methods (SAND-II Family)#
Methods: unfold_sandii, unfold_gravel, unfold_bunki,
unfold_bunkiut, unfold_rebunki, unfold_nsduaz.
The SAND-II family corrects a trial spectrum multiplicatively from the
measured-to-calculated count-rate ratios. In unfold_sandii
(McElroy et al., 1967; Griffin, Kelly and VanDenburg, 1994) the update
is the weighted geometric mean
and the iteration stops when the chi-square of the fit drops to the
number of detectors or the maximum relative change of the spectrum
falls below the tolerance. GRAVEL (unfold_gravel; Matzke, 2003)
uses the same ratio logic with count-based weights
\(W_{ij} \propto N_i\,\mathbf{R}_{ij}\,\Phi_j/\hat N_i\), an
optional regularisation term on the log-spectrum, and a stop test on
the relative change of the spectrum; its objective is the weighted sum
of squared logarithmic ratio deviations, which dampens the influence
of high-count channels. The BUNKI/BUNKI-UT variants
(unfold_bunki, unfold_bunkiut) implement the SPUNIT and BON31G
ratio recursions, unfold_rebunki adds iteration-dependent
relaxation, and unfold_nsduaz replaces the initial guess by a
catalogue of analytic shapes (fission Watt, evaporation, thermal — see
Parametric Models) before applying ratio corrections.
These methods require strictly positive measurements (zero or negative channels are dropped), preserve positivity exactly, and — as all iterative schemes here — display semi-convergence; they are typically run for a modest number of iterations with chi-square monitoring (Machado, García-Baonza and Vega-Carrillo, 2024).
Maximum Entropy and Information Divergence#
Methods: unfold_maxed, unfold_amaxed, unfold_amaxed_regularization,
unfold_imaxed, unfold_directed_divergence.
The entropy-regularised reconstruction minimises the relative entropy (Kullback and Leibler, 1951) with respect to a prior spectrum \(\boldsymbol{\Phi}_0 > 0\):
In the canonical MAXED formulation (Reginatto and Goldhagen, 1999; Reginatto, Goldhagen and Neumann, 2002) one solves the constrained problem
whose Lagrangian dual yields the solution \(\boldsymbol{\Phi}(\theta)\) parametrised by the multiplier \(\theta > 0\) attached to the chi-square constraint; \(\theta\) is then adjusted (bisection on \(\Omega\)) until \(\chi^2 \approx \Omega\), with the usual choice \(\Omega = M\) (or the \(M \pm \sqrt{2M}\) band) for known counting variances. The dual functional \(\Psi(\boldsymbol{\lambda})\) is strictly concave in the dual variables, and its Hessian, \(\mathbf{R}^T\operatorname{diag}(\boldsymbol{\Phi})\mathbf{R} + \operatorname{diag}(\sigma_i^2/2\theta)\), is positive definite for every \(\theta > 0\) thanks to its diagonal term — no additional rank condition on \(\mathbf{R}\) is required.
The implementation of unfold_maxed in this package minimises the
equivalent primal objective directly,
in log-space (\(y_j = \ln\Phi_j\), which enforces positivity
exactly) with L-BFGS-B; the relative weight of the chi-square term is
controlled by sigma_factor rather than by an explicit bisection on
\(\Omega\). unfold_amaxed, unfold_amaxed_regularization
and unfold_imaxed are adaptive/iterative refinements of the same
primal scheme (Wong, 2024). unfold_directed_divergence minimises
the directed (Kullback) divergence
\(\sum_j \Phi_j \ln(\Phi_j/\Phi_{0,j})\) with multiplicative
updates — the same functional used inside the N-spline solver
unfold_nspline (see Parametric Models).
Because the entropy functional is strictly convex on
\(\mathbb{R}^K_{>0}\) (for the Kullback-Leibler form above, the
Hessian is \(\operatorname{diag}(1/\Phi_j) \succ 0\)), the primal
problem is well posed once the feasible set is non-empty; positivity of
the solution is automatic since \(S\) is defined only for
\(\boldsymbol{\Phi} > 0\). Semi-convergence does not arise for the
strictly convex entropy-plus-quadratic objective solved to optimality,
but the noise-fitting problem reappears when \(\Omega\) (or
sigma_factor) is chosen too small.
Poisson-Likelihood (EM) Family#
Methods: unfold_mlem, unfold_mlem_stop, unfold_mlem_odl,
unfold_odl_advanced, unfold_osem, unfold_osem_anlm,
unfold_mapem, unfold_bsrem.
The Poisson log-likelihood maximisation problem is
where \(V\) is an optional roughness penalty (MAP estimation). The classical MLEM update is multiplicative,
preserves positivity exactly, and increases the likelihood
monotonically (Shepp and Vardi, 1982; Richardson, 1972; Lucy, 1974).
Because the likelihood is ill-posed, unconstrained MLEM converges to
the maximum-likelihood fit of the noise — the semi-convergence problem
again — so regularised or early-stopped variants are preferred:
ordered-subset acceleration (unfold_osem), MAP penalties
(unfold_mapem), block-sequential regularised EM (unfold_bsrem;
De Pierro, 1995) whose step sizes satisfy the stochastic-approximation
conditions \(\sum_k \alpha_k = \infty\) and
\(\sum_k \alpha_k^2 < \infty\) (Robbins and Monro, 1951), and the
proximal ODL variants (unfold_mlem_odl, unfold_odl_advanced).
OSEM with asymptotic non-local means. unfold_osem_anlm
(Jamaati et al. 2026) interleaves the ordered-subset update with the
two-stage asymptotic non-local means (ANLM) filter applied to the
intermediate spectrum after every subset update (optionally only once
after the last update, anlm_mode='post'). Stage 1 applies NLM with
the uniform parameter \(h_1 = 0.5\,\sigma\); stage 2 applies the
point-wise parameter of the article’s eq. 6,
i.e. the noise standard deviation smoothed by the initial NLM weights \(w(i,j)\). For one-dimensional spectra the 2D windows become index windows on the energy grid (search window \(N\), Gaussian-weighted similarity window \(\nu\)), and the filter by default operates in log space, where the relative EM noise is approximately additive and a single noise-level estimate matches every bin; the automatic estimate uses the MAD of the second differences.
MLEM with J-factor stopping. unfold_mlem_stop implements the
early-stopping criterion of Montgomery et al. (2020): at iteration
\(k\) the indicator
is computed and the iteration is stopped once \(J(k)\) falls below
a threshold. In the reference implementation the threshold is selected
automatically from the count level (the cps_crossover parameter,
default 30000 counts, switches between two threshold regimes); the
statistically motivated check \(J(k) \le M + \sqrt{2M}\) corresponds
to the discrepancy band discussed in Parameter Selection and Stopping Criteria. The point
estimate returned is the spectrum at the stopping index, which acts as
the regularisation parameter.
Quadratic Methods with Full Covariance Treatment#
Methods: unfold_staysl, unfold_ferdor.
This family solves weighted least-squares problems with full propagation of measurement covariances — it is not entropy-based, although historically these codes are discussed alongside the MaxEnt codes in the Bonner sphere literature (Reginatto, 2010).
STAY’SL (unfold_staysl; Perey, 1977) performs least-squares
dosimetry unfolding with input covariance matrices for the measured
responses, the response functions and the prior spectrum, and returns
the updated spectrum together with its full covariance matrix. The
solution minimises the quadratic form
subject to linear constraints, solved through the corresponding normal equations with the appropriate covariance weighting.
FERDOR (unfold_ferdor; Burrus, ORNL-4154, 1965) seeks the
spectrum that reproduces the measurements within their uncertainties
while being as smooth as possible, realised here as
where \(\boldsymbol{\Sigma}\) is the diagonal measurement-covariance matrix and \(\mathbf{D}_2\) the second difference operator. The smoothing weight \(\alpha\) is not fixed a priori: it is tuned iteratively (bisection) so that the reduced chi-square of the fit meets the discrepancy band of Parameter Selection and Stopping Criteria.
Bayesian and Stochastic Methods#
Methods: unfold_bayes, unfold_bayes_spline_regularization,
unfold_bayesian_parametric, unfold_mcmc, unfold_eki.
The posterior density combines the likelihood and prior terms,
with a Poisson or Gaussian likelihood as described in General Inverse Problem.
Bayesian iterative unfolding (D’Agostini). unfold_bayes
implements the Bayesian iterative scheme of D’Agostini (1995). With
the response normalised column-wise,
\(P(E_j \mid C_i) = \mathbf{R}_{ij}\Phi_j^{(k)} /
\sum_m \mathbf{R}_{im}\Phi_m^{(k)}\), the update is
or, equivalently, in multiplicative form
i.e. exactly the MLEM recursion without the sensitivity normalisation \(1/\sum_i \mathbf{R}_{ij}\) (the missing factor is absorbed by the column normalisation of \(\mathbf{R}\)). The two expressions above are algebraically identical; in particular the update contains the same factor \(\Phi_j^{(k)}\) only once, and no additional power of the current spectrum appears. The implementation works in effective-count space and rescales the result to physical units.
MCMC. unfold_mcmc samples the posterior with the NUTS
Hamiltonian sampler (PyMC), placing a smoothness (Ornstein-Uhlenbeck)
prior on the log-spectrum anchored at a data-driven centre. Chain
convergence is assessed with the Gelman-Rubin statistic and effective
sample sizes as computed by ArviZ. For \(m\) chains of length
\(n\), with within-chain variance \(W\) and between-chain
variance \(B\),
(Gelman and Rubin, 1992; Brooks and Gelman, 1998); the threshold \(\hat R < 1.05\) used by the package is a common practical choice, while the rank-normalised split-\(\hat R\) recommended today (Vehtari et al., 2021) — and implemented by ArviZ — tightens this to \(\hat R < 1.01\).
Ensemble Kalman inversion. unfold_eki (Iglesias, Law and
Stuart, 2013) propagates an ensemble of spectra through the Kalman
update of the forward model. In the limit \(N_{\mathrm{ens}} \to
\infty\) and \(t \to \infty\) the ensemble collapses onto the
least-squares solution of the data-misfit problem; useful regularised
solutions are obtained either by early stopping of the ensemble
evolution (the iteration index playing the role of
\(\alpha^{-1}\), with the stopping level selected by a discrepancy
or GCV rule) or by the Tikhonov-type regularisation term and prior
covariance supported by the implementation (regularization,
inflation, ensemble size n_ensemble). Convergence is monitored
through the ensemble spread and the data misfit; the spread stabilises
when the ensemble has collapsed onto the attractor of the update.
Parametric Models#
Methods: unfold_parametric, unfold_parametric2,
unfold_hybrid_parametric, unfold_crystal_ball, unfold_nspline,
unfold_zfit, unfold_lmfit, unfold_express, unfold_fruit_like,
unfold_rfsp_jul, unfold_mystic, unfold_genetic, unfold_maeo.
The spectrum is represented by a low-dimensional parameter vector \(\boldsymbol{\theta}\) and the corresponding data misfit is minimised:
Analytic building blocks used by the shape catalogues (e.g.
unfold_nsduaz, unfold_express) must be written for the fluence
density. The fission (Watt) spectrum is
with the common \({}^{252}\mathrm{Cf}\) values \(a = 1.025\) MeV, \(b = 2.926\ \mathrm{MeV}^{-1}\) (Watt, 1952); the evaporation and thermal Maxwellian fluence components enter with the same functional form \(\Phi(E) \propto E\,e^{-E/kT}\) at different temperatures, with a \(1/v\)-controlled tail in the thermal limit. Note that the form \(\sqrt{E}\,e^{-E/T}\) sometimes seen in the literature is the Maxwellian flux density convention, not the fluence density; the shapes above are the fluence forms used by the code.
The N-spline model of unfold_nspline represents
\(\ln\Phi(E)\) piecewise by splines of the form
\(\exp(a_k + q_k\ln E + r_k E)\) on segments between free knots
(Islamgulov and Lartsev, 2008); the node-continuity constraints are
enforced through KKT elimination and the directed divergence is
minimised at each iteration — see N-spline unfolding (Islamgulov & Lartsev, 2008) for the full
formulation and preset knot layouts.
These models are solved with nonlinear least squares (L-BFGS-B, least-squares
and probabilistic backends) or hybrid local/global optimisers:
unfold_lmfit and unfold_zfit wrap the corresponding fitting
libraries, while unfold_mystic, unfold_genetic and
unfold_maeo provide global and multiobjective search. For
bound-constrained smooth problems the package relies on L-BFGS-B (Byrd
et al., 1995), which converges globally to first-order stationary
points under standard line-search conditions; its observed rate is
typically linear, and no superlinear convergence is claimed for the
limited-memory variant.
Meta-Strategies and Ensembles#
Methods: unfold_ensemble, unfold_combined, unfold_composite,
unfold_cascade, unfold_iterative_refinement, unfold_binned.
Meta-methods do not introduce new variational principles; they combine
base solvers. unfold_ensemble averages (weighted mean, median,
trimmed mean, or best-residual selection) several base solutions;
unfold_combined chains methods in a pipeline; unfold_composite
and unfold_cascade run coarse-to-fine sequences;
unfold_iterative_refinement performs a two-pass unfold with an
auto-selected blending factor; and unfold_binned selects the
per-bin best method from a pre-computed benchmark lookup table. The
statistical behaviour of a meta-method is inherited from its
constituents: it can only average out solver-specific artefacts, not
replace an informed regularisation choice.
References#
Andersen, A. C., Kak, A. C. (1984). Simultaneous algebraic reconstruction technique (SART): a superior implementation of the ART algorithm. Ultrasonic Imaging 6, 81-94. doi:10.1177/016173468400600107 |
Beck, A., Teboulle, M. (2009). A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM J. Imaging Sci. 2, 183-202. doi:10.1137/080716542 |
Beck, A., Teboulle, M. (2003). Mirror descent and nonlinear projected subgradient methods for convex optimization. Operations Research Letters 31, 167-175. doi:10.1016/S0167-6377(02)00231-6 |
Boyd, S., Parikh, N., Chu, E., Peleato, B., Eckstein, J. (2011). Distributed optimization and statistical learning via the alternating direction method of multipliers. Foundations and Trends in Machine Learning 3, 1-122. doi:10.1561/2200000016 |
Boyd, S., Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press. doi:10.1017/CBO9780511804441 |
Brent, R. P. (1973). Algorithms for Minimization Without Derivatives. Prentice-Hall, Englewood Cliffs. |
Brooks, S. P., Gelman, A. (1998). General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434-455. doi:10.1080/10618600.1998.10474787 |
Burrus, W. R. (1965). Utilization of a priori information in the statistical interpretation of measured distributions. Report ORNL-4154, Oak Ridge National Laboratory. |
Byrd, R. H., Lu, P., Nocedal, J., Zhu, C. (1995). A limited memory algorithm for bound constrained optimization. SIAM J. Sci. Comput. 16, 1190-1208. doi:10.1137/0916069 |
D’Agostini, G. (1995). A multidimensional unfolding method based on Bayes’ theorem. Nucl. Instrum. Meth. A 362, 487-498. doi:10.1016/0168-9002(95)00274-X |
De Pierro, A. R. (1995). A modified expectation maximization algorithm for penalized likelihood estimation. IEEE Trans. Med. Imaging 14, 132-137. doi:10.1109/42.370409 |
Engl, H. W., Hanke, M., Neubauer, A. (1996). Regularization of Inverse Problems. Kluwer, Dordrecht. doi:10.1007/978-94-009-1740-8 |
Facchinei, F., Pang, J.-S. (2003). Finite-Dimensional Variational Inequalities and Complementarity Problems, vol. II. Springer, New York. doi:10.1007/b97411 |
Frank, M., Wolfe, P. (1956). An algorithm for quadratic programming. Naval Research Logistics Quarterly 3, 95-110. doi:10.1002/nav.3800030109 |
Gabay, D., Mercier, B. (1976). A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Computers & Mathematics with Applications 2, 17-40. doi:10.1016/0898-1221(76)90003-1 |
Gazzola, S., Hansen, P. C., Nagy, J. G. (2018). IR Tools: a MATLAB package of iterative regularization methods and large-scale test problems. Numer. Algorithms 81, 773-811. doi:10.1007/s11075-018-0570-7 |
Gelman, A., Rubin, D. B. (1992). Inference from iterative simulation using multiple sequences. Statist. Sci. 7, 457-472. doi:10.1214/ss/1177011136 |
Golub, G. H., Heath, M., Wahba, G. (1979). Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics 21, 215-223. doi:10.1080/00401706.1979.10489751 |
Golub, G. H., Kahan, W. (1965). Calculating the singular values and pseudo-inverse of a matrix. J. SIAM Numer. Anal. B 2, 205-224. doi:10.1137/0702016 |
Griffin, P. J., Kelly, D. G., VanDenburg, S. W. (1994). User’s manual for SNL-SAND-II code. Report SAND92-2357 / OSTI 10149711, Sandia National Laboratories. doi:10.2172/10149711 |
Hansen, P. C. (1992). Analysis of discrete ill-posed problems by means of the L-curve. SIAM Review 34, 561-580. doi:10.1137/1034115 |
Hansen, P. C. (1994). Regularization tools: a Matlab package for analysis and solution of discrete ill-posed problems. Numer. Algorithms 6, 1-35. doi:10.1007/BF02149761 |
Hestenes, M. R., Stiefel, E. (1952). Methods of conjugate gradients for solving linear systems. J. Res. NBS 49, 409-436. doi:10.6028/jres.049.044 |
Iglesias, M. A., Law, K. J. H., Stuart, A. M. (2013). Ensemble Kalman methods for inverse problems. Inverse Problems 29, 045001. doi:10.1088/0266-5611/29/4/045001 |
IAEA (2001). Compendium of Neutron Spectra and Detector Responses for Radiation Protection Purposes. Technical Reports Series No. 403, Vienna. |
Islamgulov, D. G., Lartsev, A. V. (2008). N-spline reconstruction of neutron spectra on the basis of activation measurements. Atomic Energy 104, 387-397. doi:10.1007/s10512-008-9045-6 |
Jaggi, M. (2013). Revisiting Frank-Wolfe: projection-free sparse convex optimization. Proc. ICML 2013, 427-435. |
Jiang, M., Wang, G. (2003). Convergence of the simultaneous algebraic reconstruction technique (SART). IEEE Trans. Image Process. 12, 957-961. doi:10.1109/TIP.2003.815295 |
Korpelevich, G. M. (1976). The extragradient method for finding saddle points and other problems. Ekonomika i Matematicheskie Metody 12, 747-756. |
Kullback, S., Leibler, R. A. (1951). On information and sufficiency. Ann. Math. Statist. 22, 79-86. doi:10.1214/aoms/1177729694 |
Lucy, L. B. (1974). An iterative technique for the rectification of observed distributions. Astron. J. 79, 745-754. doi:10.1086/111605 |
Luo, Z.-Q., Tseng, P. (1992). On the linear convergence of the coordinate descent method for convex differentiable minimization. J. Optim. Theory Appl. 72, 7-35. doi:10.1007/BF00939952 |
Machado, G. V., García-Baonza, A. M., Vega-Carrillo, H. R. (2024). Stopping criteria for neutron spectrum unfolding algorithms. Appl. Radiat. Isot. 212, 111456. doi:10.1016/j.apradiso.2024.111456 |
Matzke, M. (2003). Unfolding procedures. Radiat. Prot. Dosimetry 107, 155-174. doi:10.1093/oxfordjournals.rpd.a006384 |
McElroy, W. N., Berg, S., Crockett, T., Hawkins, R. G. (1967). A computer-automated iterative method for neutron flux spectra determination by foil activation. Report AFWL-TR-67-41 (vols. I-IV), Air Force Weapons Laboratory. |
Montgomery, L. et al. (2020). A novel MLEM stopping criterion for unfolding neutron fluence spectra in radiation therapy. Nucl. Instrum. Meth. A 957, 163400. doi:10.1016/j.nima.2020.163400 |
Morozov, V. A. (1984). Methods for Solving Incorrectly Posed Problems. Springer, New York. doi:10.1007/978-1-4612-5280-1 |
Nemirovski, A., Yudin, D. (1983). Problem Complexity and Method Efficiency in Optimization. Wiley, New York. |
Nocedal, J., Wright, S. J. (2006). Numerical Optimization, 2nd ed. Springer, New York. doi:10.1007/978-0-387-40065-5 |
Paige, C. C., Saunders, M. A. (1982). LSQR: an algorithm for sparse linear equations and sparse least squares. ACM Trans. Math. Software 8, 43-71. doi:10.1145/355984.355989 |
Perey, F. G. (1977). Least squares dosimetry unfolding: the program STAY’SL. Report ORNL/TM-6062, Oak Ridge National Laboratory. |
Polyak, B. T. (1967). A general method for solving extremal problems. Doklady Akademii Nauk SSSR 174, 33-36. |
Reginatto, M. (2010). Overview of spectral unfolding techniques and uncertainty estimation in neutron spectrometry. Radiat. Meas. 45, 1323-1329. doi:10.1016/j.radmeas.2010.06.016 |
Reginatto, M., Goldhagen, P. (1999). MAXED, a computer code for maximum entropy deconvolution of multisphere neutron spectrometry measurements. Health Phys. 77, 579-583. doi:10.1097/00004032-199911000-00012 |
Reginatto, M., Goldhagen, P., Neumann, S. (2002). Spectrum unfolding, sensitivity analysis and propagation of uncertainties with the MAXED unfolding algorithm. Nucl. Instrum. Meth. A 476, 242-246. doi:10.1016/S0168-9002(01)01439-5 |
Richardson, W. H. (1972). Bayesian-based iterative method of image restoration. J. Opt. Soc. Am. 62, 55-59. doi:10.1364/JOSA.62.000055 |
Robbins, H., Monro, S. (1951). A stochastic approximation method. Ann. Math. Statist. 22, 400-407. doi:10.1214/aoms/1177729586 |
Shepp, L. A., Vardi, Y. (1982). Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imaging 1, 113-122. doi:10.1109/TMI.1982.4307558 |
Shor, N. Z. (1985). Minimization Methods for Non-Differentiable Functions. Springer, Berlin. doi:10.1007/978-3-642-82118-9 |
Strohmer, T., Vershynin, R. (2009). A randomized Kaczmarz algorithm with exponential convergence. J. Fourier Anal. Appl. 15, 262-278. doi:10.1007/s00041-008-9030-4 |
Thomas, D. J., Alevra, A. V. (2002). Bonner sphere spectrometers — a critical review. Nucl. Instrum. Meth. A 476, 12-20. doi:10.1016/S0168-9002(01)01379-1 |
Tikhonov, A. N., Goncharsky, A. V., Stepanov, V. V., Yagola, A. G. (1995). Numerical Methods for the Solution of Ill-Posed Problems. Kluwer, Dordrecht. doi:10.1007/978-94-015-8480-7 |
Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., Bürkner, P.-C. (2021). Rank-normalization, folding, and localization: an improved \(\hat R\) for assessing convergence of MCMC. Bayesian Analysis 16, 667-718. doi:10.1214/20-BA1221 |
Watt, B. E. (1952). Energy spectrum of neutrons from thermal fission of U-235. Phys. Rev. 87, 1037-1041. doi:10.1103/PhysRev.87.1037 |
Wolfe, P. (1970). Convergence theory in nonlinear programming. In: Integer and Nonlinear Programming (J. Abadie, ed.), North-Holland, Amsterdam. |
Wong, O. (2024). Modernising neutron spectrum unfolding for fusion applications. PhD Thesis, Sheffield Hallam University. shura.shu.ac.uk/36014 |
Wright, S. J. (2015). Coordinate descent algorithms. Math. Program. 151, 3-34. doi:10.1007/s10107-015-0892-3 |
Integration Notes#
The mathematical section is referenced from the overview page through
labelled cross-references (:ref: targets math-tikhonov,
math-krylov, math-iterative, math-ratio, math-maxent,
math-em, math-quadratic, math-bayesian, math-parametric,
math-meta) and is rendered with MathJax via sphinx.ext.mathjax
using the mathjax3_config entry in conf.py (Sphinx >= 4
semantics). Inline mathematics in reStructuredText sources must use
the :math: role; the MathJax inlineMath configuration only
affects the browser-side rendering of pre-existing dollar-sign input,
not the Sphinx parsing. The optional diagram extension
sphinxcontrib-mermaid is imported defensively in conf.py and
requires the package (and, for image output, the mermaid CLI) to be
present in the documentation build environment.
Integration Notes#
The mathematical section is referenced from the overview page through
labelled cross-references (:ref: targets math-tikhonov,
math-krylov, math-iterative, math-optcourse, math-ratio,
math-maxent, math-em, math-quadratic, math-bayesian,
math-parametric, math-meta) and is rendered with MathJax via
sphinx.ext.mathjax
using the mathjax3_config entry in conf.py (Sphinx >= 4
semantics). Inline mathematics in reStructuredText sources must use
the :math: role; the MathJax inlineMath configuration only
affects the browser-side rendering of pre-existing dollar-sign input,
not the Sphinx parsing. The optional diagram extension
sphinxcontrib-mermaid is imported defensively in conf.py and
requires the package (and, for image output, the mermaid CLI) to be
present in the documentation build environment.