Interpreting Unfolding Results#

Neutron spectrum unfolding is an ill-posed inverse problem: the number of energy groups far exceeds the number of detector readings, so many spectra fit the measurements equally well. A single spectrum vector is therefore not enough to trust a result — you need to know how the solution was obtained, which constraints and measurements drive it, and how stable it is.

bssunfold provides exactly this through the interpretation API, built on the optional pyoptexplain package (post-optimality analysis of quadratic programs). Install it with:

pip install bssunfold[interpret]

The interpretation solves the same unfolding QP used by unfold_qpsolvers / unfold_cvxpy (min 1/2 xᵀQx + cᵀx with x ≥ 0) but keeps the solver internals, dual variables and perturbation results, and turns them into a human-readable report plus quantitative metrics.

Terms. A dual variable (multiplier) \(\lambda_i\) measures how much the optimal objective would improve per unit relaxation of constraint \(i\); for a minimisation problem with constraint \(x_i \ge 0\) a positive dual means the data would benefit from \(x_i < 0\). Shadow price is the practical name of such a dual in economics and optimisation (Boyd and Vandenberghe, 2004; doi:10.1017/CBO9780511804441). The condition number of the QP matrix, the active set (constraints holding with equality at the optimum) and the regularisation sweep are defined in Parameter Selection and Stopping Criteria.

For the Bayesian counterpart of these diagnostics — global uncertainty per bin instead of local perturbations — see the CUQIpy interpretability section in CUQIpy Bayesian unfolding: uncertainty-quantified MCMC.

What the interpretation answers#

        graph TD
    A["Interpretation (pyoptexplain)"] --> B["Solve report"]
    A --> C["Shadow prices (duals)"]
    A --> D["Robustness"]
    A --> E["Detector sensitivity"]
    A --> F["Regularization sweep"]
    A --> G["Non-negativity trust"]
    A --> H["Scenarios"]

    B --> B1["solver status, objective"]
    B --> B2["active / zeroed energy groups"]
    C --> C1["price of each x ≥ 0 bound"]
    C --> C2["price of sum(x) == norm"]
    D --> D1["spectrum change for ±1…5% readings"]
    E --> E1["importance per detector sphere"]
    F --> F1["solution across a grid of α"]
    G --> G1["effect of allowing small negatives"]
    H --> H1["structured what-if cases"]

    style A fill:#4a90d9,color:#fff

Solve report — solver status, objective value, condition number and which energy groups are active (zeroed at the bound x = 0).
Shadow prices (duals) — how much the objective would change if a zeroed group were allowed to become positive. Large prices flag energy ranges the data “wants” but the non-negativity constraint keeps at zero.
Robustness — how the unfolded spectrum moves when the readings are perturbed by ±1…5%. A stable solution barely changes.
Detector sensitivity / importance — how the spectrum changes when one detector reading at a time is perturbed. This ranks the spheres by how much information they carry about the solution.
Regularization sweep — how the solution and residual change across a grid of the regularization parameter α, revealing whether the result is α-sensitive.
Non-negativity trust — what happens if small negative values are allowed, i.e. whether x ≥ 0 is forcing the solution into an artificial corner.
Scenarios — structured what-if cases over a custom scenario space.

API overview#

Level	Callable	Purpose
High-level	`Detector.unfold_interpret`	Unfold the readings and append `report` and `interpretation_metrics` to the standard result dict.
Analysis-only	`Detector.interpret_result`	Run the interpretation directly on readings without the unfolding bookkeeping; returns `report`, `metrics` and `tables`.
Low-level	`bssunfold.core.unfold_interpret.interpret_qp`	Full QP interpretation from a response matrix `A` and readings `b`.
Low-level	`bssunfold.core.unfold_interpret.build_interpretation_qp`	Build the QP (matrix `Q`, vector `c`, bounds) without solving.
Low-level	`bssunfold.core.unfold_interpret.solve_interpret`	Solve a QP and collect solver/duals/perturbation diagnostics.

Basic usage#

from bssunfold import Detector, RF_LANL

detector = Detector(RF_LANL)
readings = {"3in": 0.5, "5in": 1.2, "8in": 2.1, "12in": 3.4}

# 1) Unfold and interpret in one call
result = detector.unfold_interpret(readings, tolerance=1e-5)

spectrum  = result["spectrum"]              # unfolded flux per energy bin
report    = result["report"]                # full Markdown report
metrics   = result["interpretation_metrics"]  # JSON-friendly diagnostics

# 2) Analysis-only: same interpretation without the unfolding wrapper
ir = detector.interpret_result(readings, tolerance=1e-5)
tables = ir["tables"]  # pandas DataFrames: summary, duals, detectors, ...

# 3) L1 norm with automatic ridge for numerical stability
result = detector.unfold_interpret(
    readings, norm=1, ridge_coeff="auto",  # default; or 0.0 to disable
)

Note

About ``tolerance``: pyoptexplain’s backend can report iteration_limit on the full 11-sphere RF_LANL problem at the strictest default tolerance (1e-8). Relaxing the feasibility/optimality tolerance to 1e-5 is enough for an optimal status with a residual below 0.3 % — a good compromise for analysis-grade work. See the worked example in examples/24-interpret.ipynb.

Note

About ``norm=1`` and the ``ridge_coeff`` parameter: when the L1 penalty norm is used (norm=1), the QP matrix P = A'A can be rank-deficient if the number of detectors is much smaller than the number of energy bins. This causes OSQP to stall at the iteration limit. To cure this, build_interpretation_qp adds a small diagonal ridge to P by default (ridge_coeff="auto", which computes 1e-8 * trace(P) / n). The ridge makes P positive definite without noticeably affecting the L1 solution. You can override it with an explicit float (e.g. ridge_coeff=1e-6) or disable it entirely with ridge_coeff=0.0.

How to use it to interpret a spectrum#

Once you have result["report"], read it in this order:

Check the solve status and residual. If the solver is not optimal, or the residual is much larger than the measurement noise, the response matrix and the readings are inconsistent — check the detector set and the readings before trusting any spectrum shape.
Look at the active (zeroed) groups. Unfolding is sparse by nature: many energy groups end up pinned at x = 0. This tells you where the solution is determined by the data and where it is only constrained by the non-negativity bound. A spectrum that is concentrated in a few groups is a sparse, low-information solution; a “washed-out” spectrum with everything active usually means the data cannot resolve the energy structure.
Read the shadow prices. A large dual on a zeroed group means the data would prefer a positive flux there but the constraint forbids it. If the largest duals cluster in one energy range (e.g. a fast peak or a thermal component), that range is the least constrained part of the spectrum — the measurement set has little resolving power there, so the value is an artifact of the regularization, not of the data.
Check robustness. Compare the relative spectrum change under ±1…5 % reading perturbations. If the solution moves by more than the perturbation level, the unfolding is unstable for this detector set and you should add spheres (or a prior) before drawing quantitative conclusions.
Check non-negativity trust. If allowing 5 % negative values moves the spectrum only slightly, the x ≥ 0 bound is not distorting the solution. A large change means the physical spectrum actually needs negative components — typically a sign that the response functions or readings are miscalibrated.

How to use it to choose the detector sphere set#

The interpretation is also a sphere-selection tool. The detector_importance table ranks the spheres by how much the unfolded spectrum changes when each reading is perturbed on its own:

High-importance spheres dominate the solution — the spectrum is tightly pinned by their readings. In the RF_LANL Cf-252 example the lead-shielded and mid-size spheres (9inPb, 12inPb, 8in) carry most of the information.
Low-importance spheres (e.g. the small 3in sphere on a hard fast spectrum) barely influence the solution. Their readings are consistent with the rest of the set but add almost no information — they can be removed without changing the unfolded spectrum, saving measurement time.
High-residual spheres (from the per-detector residual table) disagree with the reconstructed spectrum. Re-examine their calibration or response function before trusting the result.

The regularization sweep complements this: if the residual stays flat over a wide range of α, the sphere set already provides enough information and the choice of regularization is not critical. If the residual changes sharply with α, the set is under-determined and more (or better chosen) spheres are needed.

Worked example#

Run the notebook examples/24-interpret.ipynb for a complete, reproducible walkthrough on the built-in RF_LANL response functions with the ISO Cf-252 reference spectrum. It produces the full Markdown report, the quantitative metrics tables and the key diagnostic plots (unfolded spectrum vs. reference, per-detector fit, detector importance, sensitivity and non-negativity).