B-spline MLEM unfolding (MLEM-BS)#

The unfold_mlem_bs method implements the MLEM-BS algorithm of

V. Mazankova, L. Torokova, D. Trunec, Z. Kopecky, Z. Matej, Experimental Measurement of Neutron Flux and Its Mathematical Data Processing, Proceedings of the 5th International Conference CNDGS’2026, Brno, Czech Republic, 7–10 September 2026, ISSN 2538-8959, https://doi.org/10.47459/cndcgs.2026.61.

Method outline#

The measurement is described by the discrete Fredholm integral equation of the first kind with Poisson-distributed counts,

\[n_i = \sum_j R_{ij}\, x_j, \qquad i = 1,\dots,I. \qquad (2)\]

Instead of solving for the \(n\) bin values \(x_j\) directly, the paper represents the sought spectrum in a B-spline basis,

\[x(E) = \sum_{s=1}^{N_s} b_s\, B_s(E),\]

so that the effective system matrix becomes \(\widetilde{R} = R\,B\) (build_bspline_basis() assembles the design matrix \(B\); clamped uniform or logarithmic knot grids are supported, selected automatically from the span of the energy grid). The B-spline coefficients \(b_s\) are found with the regularized MLEM iteration (Eq. 4 of the paper):

\[b_s^{(k+1)} = \frac{b_s^{(k)}}{\sum_i (RB)_{is} + \beta\, \partial P/\partial b_s} \sum_i (RB)_{is}\, \frac{n_i}{\sum_{s'} (RB)_{is'}\, b_{s'}^{(k)}},\]

where the penalty (Eq. 5) suppresses the noise influence,

\[P(b) = \|D^{(2)} b\|_2^2, \qquad \frac{\partial P}{\partial b} = 2\, (D^{(2)})^T D^{(2)} b,\]

with \(D^{(2)}\) the second-derivative (second finite difference) matrix (second_difference_matrix()).

The solution is restricted to the B-spline sieve — the cone of non-negative coefficients — following Z. Szkutnik, B-splines and discretization in an inverse problem for Poisson processes, Journal of Multivariate Analysis 93, 198–221 (2005) (reference [7] of the paper). Positivity of the spectrum is preserved automatically by the multiplicative MLEM update started from a non-negative sieve projection of the initial spectrum. In the overdetermined regime of the paper (thousands of measured bins vs \(N_s\) splines) the sieve restriction also provides a strong noise-suppressing regularisation by itself, even with \(\beta = 0\).

Penalty strength#

The absolute penalty parameter \(\beta\) of Eq. 4 is problem-scale dependent (the paper uses \(\beta = 1.0 \times 10^{-17}\) for its counting setup). For convenience the method additionally accepts beta_relative, which defines the effective penalty as \(\beta = \beta_{rel}\, \bar{s}\) with \(\bar{s}\) the mean column sum of the effective matrix \(RB\), i.e. the penalty gradient is measured relative to the MLEM sensitivity. Exactly one of beta / beta_relative may be provided; with neither, the iteration reduces to a pure sieve MLEM on the B-spline parameterisation.

Automatic parameter selection#

Following the paper, the number of iterations and the parameters \(N_s\) (B-spline space dimension) and \(\beta\) (penalty strength) are selected by minimization of the goodness-of-fit statistic (Eq. 6):

\[K_S = \left| \frac{\sum_i \left(n_i - \sum_s (RB)_{is}\, b_s^{(k)}\right)^2} {\sum_i \sum_s (RB)_{is}\, b_s^{(k)}} - 1 \right|,\]

which approaches zero when the residuals of the fit are consistent with Poisson noise. With auto_params=True the solver scans a grid of candidate \(N_s\) and relative penalty strengths (AUTO_BETA_RELATIVE_GRID), runs the iteration for each candidate, and keeps the combination with the smallest \(K_S\); the selected values are reported in the auto_selection key of the result. Even in manual mode the iteration tracks the \(K_S\) history and keeps the best iterate (early stop with ks_patience iterations without improvement).

Confidence intervals#

The paper estimates confidence intervals of the unfolded spectrum with a Poisson bootstrap (Eqs. 7–9): the backward-reconstructed counts \(n^{(0)} = R\, x^{(0)}\) are resampled as \(n^{*(b)}_i \sim \mathrm{Poisson}(n^{(0)}_i)\), each replicate is unfolded with the same settings, and the \(100(1-\alpha)\%\) interval is formed from the \(\alpha/2\) and \(1-\alpha/2\) quantiles of the bootstrap distribution (a 95% CI corresponds to \(\alpha = 0.05\)). Enable it with bootstrap_ci=True; the result then contains ci_low, ci_high, bootstrap_mean and bootstrap_std arrays.

Note

The count-resampling bootstrap captures the statistical (Poisson) uncertainty. In strongly underdetermined BSS systems (few spheres, many bins) the systematic model bias dominates and the intervals become optimistic; smaller \(N_s\), auto_params=True or an informative initial_spectrum make the intervals more conservative (see the module docstring of bssunfold.core.unfold_mlem_bs).

Usage#

from bssunfold import Detector

detector = Detector()
readings = {"3in": 0.053, "5in": 0.184, "10in": 0.172, "18in": 0.034}

result = detector.unfold_mlem_bs(
    readings,
    n_basis=None,          # None -> automatic (log knots for wide grids)
    beta_relative=None,    # None -> pure sieve MLEM (beta = 0)
    max_iterations=500,
    bootstrap_ci=True,     # Poisson bootstrap CI (Eqs. 7-9)
    n_bootstrap=100,
    random_state=42,
)

spectrum = result["spectrum"]      # unfolded fluence spectrum
ks = result["ks_final"]            # K_S statistic of Eq. 6
ci = (result["ci_low"], result["ci_high"])

# fully automatic selection of N_s, beta and the iteration count
result = detector.unfold_mlem_bs(readings, auto_params=True)
print(result["auto_selection"]["chosen"])

Diagnostics returned by the method: ks_history (per-iteration \(K_S\)), ks_final, chi2_pearson, n_basis, interior_knots, beta_effective, coefficients (non-negative B-spline coefficients), and, with auto_params=True, auto_selection with the full candidate table.