SSR unfolding: Sign-Simplicity-Regression solver (sisireg port)#

The unfold_ssr method is the Python port of the R package sisireg 1.2.1 (L. Metzner, sisireg: Sign-Simplicity-Regression- Solver, CRAN, GPL>=2). The Sign-Simplicity-Regression (SSR) model of Metzner (2020, Trendbasierte Prognostik; 2021, Adaequates Maschinelles Lernen) is a nonparametric regression built on the signs of the residuals: it seeks the most parsimonious (fewest extrema) regression function that is statistically adequate with respect to two sign criteria.

Sign criteria#

For residuals r_i = y_i - mu_i and their signs s_i = sign(r_i) the model must satisfy

  • the partial sum criterion: for every interval length k the maximum absolute sum of consecutive residual signs stays below the 95% quantile of the partial sums,

    \[\max_t \Big| \sum_{j=t}^{t+k-1} s_j \Big| \;\le\; F(n, k) = \min\!\big(\sqrt{1 + 2.33 \ln(n)\, k},\; k\big);\]

  • the maximum run criterion: with the 95% quantile of the maximum run length of equal signs,

    \[k_{\mathrm{run}}(n) = \lfloor 3.3 + 1.44 \ln(n) \rfloor,\]

    no window of \(k_{\mathrm{run}} + 1\) consecutive residuals may sum to more than \(k_{\mathrm{run}}\) in absolute value.

The regression function itself is computed with the quantised Gauss-Seidel (QSOR) iteration of the package: every interior point is replaced by the simplicitic linear interpolation of its neighbours and the update is reverted whenever it would violate the partial sum criterion with threshold \(h\) (fn). Iterating this projection yields the most parsimonious statistically adequate function. The port reproduces the R/C semantics exactly, including the rolling-median start values, the truncating integer conversions of the .C interface and both the L1 (ssr/ssr_ne) and the L2 standardised-stencil (ssr(..., funk=2)) solvers.

Unfolding scheme#

The SSR model regresses directly observed data, while unfolding only observes the folded readings \(b = A x\). solve_ssr therefore alternates

  1. a data-fidelity step — one multiplicative MLEM update of the spectrum (positivity preserving), and

  2. an SSR parsimony step — a non-equidistant SSR QSOR sweep of the current spectrum over the energy grid, which removes sign-inadequate wiggles while bounding the deviation from the data-fit in every window.

The strength of the parsimony step is the partial sum threshold \(h\). Following Metzner’s minimum statistic, fn="auto" starts from

\[h_{\mathrm{start}} = \lfloor 0.66\, F(n, k_{\mathrm{run}}) \rfloor,\]

decreases the threshold by one and stops at the first candidate whose folded residuals violate the data-space sign adequacy (partial sum / maximum run test applied to \(b - A x\)) or whose parsimony (number of extrema) exceeds that of the start model; the last adequate candidate is returned — the analogue of the R ssr_ne_minR loop transplanted to the unfolding system. A fixed threshold can be forced with fn=<int>.

Usage#

from bssunfold import Detector

det = Detector()
result = det.unfold_ssr(readings)             # fn="auto"
print(result["fn"], result["n_extrema"],
      result["ps_valid_data"], result["run_valid_data"])

# fixed parsimony threshold
result = det.unfold_ssr(readings, fn=3)

The pure regression building blocks are exported as well and can be used standalone:

from bssunfold.core import ssr, ssr_ne, ssr_min_statistic_ne
from bssunfold.core import max_run_quantile, partial_sum_valid
from bssunfold.core import ssr_predict

Diagnostics returned with the standard result dictionary: fn (threshold used), fn_start (ladder start), k_run (run length quantile), n_extrema (parsimony of the unfolded spectrum), ps_valid_data / run_valid_data (data-space adequacy of the folded residuals), ssr_converged and, from solve_ssr_full, the full fn_ladder trail with per-candidate adequacy results.

Notes#

  • Requires at least 8 energy bins (rolling-median windows) and 3 detector readings; no upper limit — the QSOR sweeps are O(n) sequential passes.

  • Non-negativity is preserved by construction (interpolation of a non-negative spectrum), so the method suits spectra with strong low-energy components.

  • For fewer than 25 detector readings the partial sum test is trivially satisfied (its longest tested interval is m // 5); the selection is then driven by the maximum run test and the parsimony criterion.

  • Pure NumPy — no additional dependencies. The R package sisireg is GPL (>= 2), compatible with the GPL-3 license of bssunfold.

Spatial regression: ssr3d#

The 1-D solvers generalise to scattered data \(z(x, y)\) on the plane (R ssr3d, R/ssr3d.R + src/ssr3d.c). The regression surface is the most parsimonious minimal surface that is statistically adequate with respect to the residual signs in the k-quadrant-neighbourhood of every observation: the k nearest observations in each of the four quadrants north-east, north-west, south-west and south-east of a reference point. Adequacy is tested over the 4k + 1 nearest-neighbourhoods (the partial sum criterion with threshold fn); the Gauss-Seidel iteration replaces every point by the exponential-distance-weighted mean of its quadrant neighbours and reverts the update to the observation whenever the new value violates the criterion.

from bssunfold.core import ssr3d, ssr3d_predict, ps_statistic_3d

model = ssr3d(coords, dat)          # coords (n, 2), dat (n,)
z = ssr3d_predict(model, queries)   # exponential weights
z_ms = ssr3d_predict(model, queries, ms=True)  # 4-point minimal surface

# partial sums vs. 95% quantiles (computation behind R psplot3d)
ps, fn = ps_statistic_3d(coords, dat, model.mu)

Neighbourhood construction (near_neighbors_quadrant, near_neighbors), the partial sum statistic (ps_max_3d, ps_statistic_3d), the weighted means (wmean, wmean_exp, wmean_ms) and the model dataclass SSR3DModel are exported as well. The default k = max_run_quantile(n) / 2 reproduces the R true division including its truncating conversions. The port was verified to reproduce the original R output to machine precision (same values of mu for scattered and jittered-grid data); on exactly regular grids, distance ties make the neighbourhood composition depend on the last bits of the floating-point distances, so tiny deviations of the tie handling are possible there — as between any two R builds.

Neural regression: ssrMLP#

ssrMLP (R ssrMLP.R) is a two-hidden-layer perceptron with sigmoid activations and linear output whose SGD training criterion is not the least squares residual but Metzner’s partial sum criterion: the sum of the positive parts of the absolute input-neighbourhood residual sign sums above the threshold, optionally combined with squared residuals (opt='ps_lse'), with a curvature penalty of the learned function over the input neighbourhoods (opt='ps_l1'), with plain least squares (opt='lse') or with user-supplied error and factor functions (opt='ext'). The port reproduces the R semantics exactly, including two quirks of the original: the hidden layer receives the transposed weight update (t(t(temp) %*% O1)), which confines the configuration to square hidden layers hl[0] == hl[1], and the fn / alpha arguments of ssrmlp_train are accepted but never forwarded, so the criterion defaults fn = 4 / alpha = 1e-4 apply (custom values are possible through opt='ext').

from bssunfold.core import ssrmlp_train, ssrmlp_predict, fii_model

model = ssrmlp_train(X, Y, opt="ps", max_iter=1000, rng=0)
yp = ssrmlp_predict(X, model)
importance = fii_model(model)   # per-input (bias included), sums to 1

The criterion building blocks (check_ps, err_ps/fac_ps, err_lse/fac_lse, err_ps_lse/fac_ps_lse, err_ps_l1/fac_ps_l1), the forward pass calc_out and the per-sample factor importance fii_prediction are exported as well. Two training epochs reproduce the original R result bit-exactly (fixture-tested); longer runs can diverge in the last bits because the sign criterion is discontinuous — a residual within BLAS noise of zero can flip its sign — which is inherent to any cross-platform reimplementation of the R algorithm.