A Linear Iterative Unfolding Method

A frequently faced task in experimental physics is to measure the probability distribution of some quantity. Often this quantity to be measured is smeared by a non-ideal detector response or by some physical process. The procedure of removing this smearing effect from the measured distribution is called unfolding, and is a delicate problem in signal processing, due to the well-known numerical ill behavior of this task. Various methods were invented which, given some assumptions on the initial probability distribution, try to regularize the unfolding problem. Most of these methods definitely introduce bias into the estimate of the initial probability distribution. We propose a linear iterative method, which has the advantage that no assumptions on the initial probability distribution is needed, and the only regularization parameter is the stopping order of the iteration, which can be used to choose the best compromise between the introduced bias and the propagated statistical and systematic errors. The method is consistent: “binwise” convergence to the initial probability distribution is proved in absence of measurement errors under a quite general condition on the response function. This condition holds for practical applications such as convolutions, calorimeter response functions, momentum reconstruction response functions based on tracking in magnetic field etc. In presence of measurement errors, explicit formulae for the propagation of the three important error terms is provided: bias error, statistical error, and systematic error. A trade-off between these three error terms can be used to define an optimal iteration stopping criterion, and the errors can be estimated there. We provide a numerical C library for the implementation of the method, which incorporates automatic statistical error propagation as well.

💡 Research Summary

The paper addresses the classic unfolding problem in experimental physics, where a measured probability density function (pdf) g(y) is a smeared version of the true underlying pdf f(x) due to detector response or physical processes. Mathematically this is expressed as a linear integral equation g = A_ρ f, where A_ρ is the folding operator defined by the response function ρ(y|x). Because A_ρ is continuous while its inverse is generally discontinuous, direct inversion is ill‑posed and leads to wildly oscillating, non‑physical solutions when applied to noisy data.

The authors propose a linear iterative unfolding algorithm that requires no prior assumptions about the shape of f. The method is motivated by the Neumann series and Landweber iteration from functional analysis, but is specially adapted to probability‑theoretic settings where the usual convergence conditions do not hold. The key steps are:

1. Normalization constant
   K_ρ = ∬ ρ(y|z) ρ(y|x) dy dz, assumed finite for the response under consideration (holds for convolutions, calorimeter responses, tracking‑based momentum reconstructions, etc.).

2. Initial estimate
   f₀ = K_ρ⁻¹ A_ρᵀ g, where A_ρᵀ is the transpose (or adjoint) of the folding operator.

3. Iterative update
   f_{n+1} = f_n + ( f₀ – K_ρ⁻¹ A_ρᵀ A_ρ f_n ).

The iteration proceeds until a stopping criterion is met. The authors prove a convergence theorem: for any compact region S, the average deviation of f_n from the orthogonal projection of the true pdf onto the complement of the kernel of A_ρ tends to zero as n → ∞, provided K_ρ < ∞. If A_ρ is invertible, the full pdf is recovered; otherwise the algorithm reconstructs the maximal information possible, i.e., f – P_{Ker(A_ρ)} f.

Three distinct error contributions are treated analytically:

• Bias error – the deterministic distance between the current iterate and the best achievable reconstruction. Theorem 2 gives an upper bound that decays roughly as 1/√(bin volume) with increasing iteration order, showing that bias diminishes with more iterations.

• Statistical error – arising from the finite‑sample fluctuations in the measured histogram g. By factorizing the covariance matrix C = EEᵀ, the same iterative scheme is applied to the matrix E, yielding a propagated covariance C_n = E_n E_nᵀ at each step (Theorem 3). This exact propagation is a direct consequence of the linearity of the method.

• Systematic error – due to uncertainties in the response function ρ. Because the algorithm is linear in ρ, a perturbation δρ can be propagated through the same iteration, allowing a quantitative estimate of the systematic contribution.

These three components combine to a total error estimate TotalError(N) = Bias(N)² + Stat(N)² + Syst(N)², which can be evaluated as a function of the iteration number N. The optimal stopping criterion is defined as the N that minimizes this total error, providing a principled balance between decreasing bias and increasing propagated statistical/systematic uncertainties.

The paper also discusses practical aspects:

• The method works with histograms, but does not require a specific binning; any density estimator (e.g., kernel density) can be used, and the binning effect itself can be folded into ρ and unfolded simultaneously.
• For discrete implementations, A_ρ and A_ρᵀ become matrices; the algorithm reduces to simple matrix‑vector operations.
• A C library implementing the algorithm is supplied, featuring automatic statistical error propagation and example programs. Users can supply custom implementations of A_ρ and A_ρᵀ to handle complex detector responses.

In summary, the authors deliver a transparent, mathematically rigorous unfolding technique that:

1. Guarantees convergence under a broad class of realistic response functions.
2. Provides explicit formulas for bias, statistical, and systematic error propagation.
3. Requires only a single regularization parameter—the iteration count—facilitating an objective stopping rule.
4. Is implemented in a lightweight, portable C library suitable for integration into existing analysis frameworks.

Compared with traditional non‑linear Bayesian iterative unfolding, the proposed linear method retains positivity and integral conservation only implicitly (through the stopping rule) but gains the crucial advantage of exact error propagation and a clear convergence proof. This makes it a compelling alternative for high‑energy physics, nuclear physics, and any domain where deconvolution of smeared probability distributions is required.