P-spline mixed-model unfolding with REML smoothing selection#

The unfold_pspline_reml method is the Python analogue of the R package LMMsolver (M. P. Boer, LMMsolver: Linear Mixed Models with Sparse Matrix Methods, CRAN, 2023; M. P. Boer, Tensor product P-splines using a sparse mixed model formulation, Stat. Model. 23(5-6), 465-481, 2023): the unfolded spectrum is a P-spline (Eilers & Marx, Flexible smoothing with B-splines and penalties, Stat. Sci. 11(2), 89-121, 1996) whose smoothness is selected automatically by restricted maximum likelihood (REML) in a linear mixed-model (LMM) formulation (Wand & Ormerod, On semiparametric regression with O’Sullivan penalized splines, Aust. N. Z. J. Stat. 50(2), 179-198, 2008).

Method outline#

The spectrum is parameterised with a B-spline basis,

\[x(E) = \sum_{s=1}^{N_s} c_s\, B_s(E) \;=\; B c,\]

so the Fredholm system \(b = A x\) becomes the linear mixed model

\[b = A B c + \varepsilon, \qquad c \sim \mathcal{N}\!\left(0,\; \sigma_e^2\,\lambda\, G^{-1}\right), \qquad G = D^{(d)\,T} D^{(d)},\]

where \(D^{(d)}\) is the \(d\)-th order difference matrix (difference_matrix(), the classic P-spline penalty) and \(\lambda\) is the smoothing parameter. The penalty eigen-decomposition \(G = U \operatorname{diag}(g) U^T\) splits the coefficient space into

  • a fixed part spanning the null space of \(G\) (dimension \(d\)) — an unpenalised polynomial trend, \(X = A B U_{\mathrm{fixed}}\);

  • a random part spanning the range space with prior precision \(L = \operatorname{diag}(g_{\mathrm{random}}) > 0\), \(Z = A B U_{\mathrm{random}}\).

For a trial smoothing parameter \(\lambda\) the coefficients follow from the Henderson mixed-model equations (the system solved by LMMsolver::LMMsolve):

\[\begin{split}\begin{bmatrix} X^T W X & X^T W Z \\ Z^T W X & Z^T W Z + \lambda L \end{bmatrix} \begin{bmatrix} \hat\beta \\ \hat b \end{bmatrix} = \begin{bmatrix} X^T W y \\ Z^T W y \end{bmatrix},\end{split}\]

with per-detector weights \(W\) (uniform by default, or \(w_i = 1/b_i\) for counting statistics — weights="poisson"). The spectrum estimate is \(\hat x = X\hat\beta + Z\hat b\).

The smoothing parameter is estimated by maximising the REML profile log-likelihood of the marginal model \(y \sim \mathcal{N}(X\beta,\; \sigma_e^2 V)\), \(V = I + \lambda Z L^{-1} Z^T\):

\[\ell_R(\lambda) = -\tfrac{1}{2}\Big[ (m - p_f) \log \hat\sigma_e^2(\lambda) + \log|V| + \log|X^T V^{-1} X| \Big],\]

optimised with Brent’s method. Because \(G\) and the data term live on very different scales, the search is performed on a scale-free relative grid, \(\lambda = \lambda_{\mathrm{ref}}\,\lambda_{\mathrm{rel}}\) with \(\lambda_{\mathrm{ref}}\) equalising the average trace of both terms and \(\lambda_{\mathrm{rel}} \in [10^{-6}, 10^{6}]\); both \(\lambda\) and \(\lambda_{\mathrm{rel}}\) are reported. A fixed relative value can be forced with lam_relative.

Usage#

from bssunfold import Detector

det = Detector()
result = det.unfold_pspline_reml(readings, weights="poisson")
print(result["lam_relative"], result["ed"], result["reml_loglik"])

Diagnostics returned with the standard result dictionary: lam (absolute smoothing), lam_relative, lam_ref, sigma2 (residual variance estimate), reml_loglik, ed and ed_norm (effective dimension of the fit and its fraction of the basis size), reml_converged. The effective dimension

\[\mathrm{ed} = p_f + \operatorname{tr}\!\left[ (Z^T W Z + \lambda L)^{-1} Z^T W Z \right]\]

quantifies the flexibility actually used by the fit.

Notes#

  • Requires at least diff_order + 2 detector readings.

  • The solution is linear in the data for a fixed \(\lambda\) (like TSVD); non-negativity is enforced by clamping, so the method is best suited for smooth spectra without sharp low-energy edges.

  • Only NumPy/SciPy are used — no additional dependencies.