Data Unfolding with Wiener-SVD Method

TL;DR

This paper introduces a novel data unfolding method in high-energy physics that combines SVD and Wiener filtering to improve signal extraction without needing regularization parameters.

Contribution

The Wiener-SVD unfolding technique is a new approach that maximizes signal-to-noise ratio in the frequency domain, inspired by digital signal processing methods.

Findings

Effective in maximizing signal-to-noise ratio

Free from regularization parameter

Discusses advantages and limitations

Abstract

Data unfolding is a common analysis technique used in HEP data analysis. Inspired by the deconvolution technique in the digital signal processing, a new unfolding technique based on the SVD technique and the well-known Wiener filter is introduced. The Wiener-SVD unfolding approach achieves the unfolding by maximizing the signal to noise ratios in the effective frequency domain given expectations of signal and noise and is free from regularization parameter. Through a couple examples, the pros and cons of the Wiener-SVD approach as well as the nature of the unfolded results are discussed.

Figures

Tables

Table 1. Table 1 : Explanations of main symbols used in this section.

Symbols	Meaning	Dimension and Format
$A$	Constructed smearing matrix from $W$ and $V$	$n\times n$ matrix
$A_C$	Constructed smearing matrix with $C$	$n\times n$ matrix
$C$	Assisting matrix, i.e. 1st or 2nd derivative matrix	$n\times n$ matrix
$Cov$	Covariance matrix of measurement $𝐦$	$m\times m$ matrix
$Co{v}_{\hat{s}}$	Covariance matrix of unfolded result $\hat{s}$	$n\times n$ matrix
$D$	Center diagonal matrix from SVD decomposition of $R$	$m\times n$ matrix
$D_C$	Center diagonal matrix from SVD decomposition of $R\cdot C^{-1}$	$m\times n$ matrix
$d_i$	Non-negative diagonal element of $D$, $d_i$ = $D_{ii}$
$F$	regularization filter	$n\times n$ matrix
$𝐦$	Measured spectrum	$m$ vector
$M$	Measured spectrum after pre-scaling	$m$ vector
$M_U$	$U^T\cdot M$	$m$ vector
$\overline{M_U}$	Expectation of $U^T\cdot M$ based on $\overline{s}$	$m$ vector
$N$	“Noise” of measurement after pre-scaling	$m$ vector
$N_U$	$U^T\cdot N$	$m$ vector
$Q$	Lower triangular matrix from Cholesky decomposition of $Co{v}^{-1}$	$m\times m$ matrix
$𝐫$	Smearing/response function	$m\times n$ matrix
$R$	Smearing/response matrix after pre-scaling	$m\times n$ matrix
$R_{tot}$	Total transformation matrix connecting $\hat{s}$ and $𝐦$	$n\times n$ matrix
$s$	Unknown spectrum (a variable in the ${\chi }^2$ calculation)	$n$ vector
$s_{true}$	True spectrum	$n$ vector
$\overline{s}$	Expectation of true spectrum	$n$ vector
$\hat{s}$	Unfolded spectrum	$n$ vector
$T_{bias}$	Bias of unfolded result	$n$ vector
$T_{deviation}$	Deviation (square-root of variance) of unfolded result	$n$ vector
$T_{deviationj}$	$j$th element of deviation
$U$	Left U matrix from SVD decomposition of $R$	$m\times m$ matrix
$U_C$	Left matrix from SVD decomposition of $R\cdot C^{-1}$	$m\times m$ matrix
$V^T$	Right $V^T$ matrix from SVD decomposition of $R$	$n\times n$ matrix
$V_C^T$	Right matrix from SVD decomposition of $R\cdot C^{-1}$	$n\times n$ matrix
$W$	Wiener filter	$n\times n$ matrix
$W_C$	Wiener filter with $C$	$n\times n$ matrix

Equations

$$M(t^{\prime })={\int }_{-\infty }^{\infty }R(t,t^{\prime })\cdot S(t)\cdot dt+N(t^{\prime }).$$

$$R(t,t^{\prime })\equiv R(t^{\prime }-t),$$

$$M(\omega )=R(\omega )\cdot S(\omega )+N(\omega ),$$

$$\hat{S}(\omega )=\frac{M(\omega )}{R(\omega )}=S(\omega )+\frac{N(\omega )}{R(\omega )}.$$

$$\hat{S}(\omega )=\frac{M(\omega )}{R(\omega )}\cdot F(\omega ).$$

$$F(\omega )=\frac{\overline{R^2(\omega )\cdot S^2(\omega )}}{\overline{R^2(\omega )\cdot S^2(\omega )}+\overline{N^2(\omega )}},$$

$$MSE$$

$${\chi }^2\left(s\right)={\left(\mathbf{m}-\mathbf{r}\cdot s\right)}^T\cdot Co{v}^{-1}\cdot \left(\mathbf{m}-\mathbf{r}\cdot s\right).$$

$$Co{v}^{-1}=Q^T\cdot Q,$$

$${\chi }^2={\left(M-R\cdot s\right)}^T\cdot \left(M-R\cdot s\right)=\sum _i{\left(M_i-\sum _j{R}_{ij}\cdot s_j\right)}^2,$$

$$\hat{s}={\left(R^{T}R\right)}^{-1}\cdot R^T\cdot M.$$

$$\hat{s}={\left(R^{T}R\right)}^{-1}{R}^T\cdot \left(R\cdot s_{true}+N\right).$$

$$R=U\cdot D\cdot V^T,$$

$$\hat{s}$$

$$\varphi (s)=\log L(s)+\tau \Sigma (s),$$

$${\Sigma }^k\left(s(E)\right)=-\int {\left(\frac{{\text{d}}^{k}s(E)}{{\text{d}}^{k}E}\right)}^2\text{d}E,k=0,1,2,3,\cdots ,$$

$$\hat{s}$$

$$A=V\cdot F\cdot V^T,$$

$$F_{ii}=\frac{d_i^2}{d_i^2+\tau }$$

$$\hat{s}$$

$$\overline{M_U}=U^T\cdot \overline{M}=U^T\cdot R\cdot \overline{s}=D\cdot V^T\cdot \overline{s}.$$

$$‘‘\overline{R^2\cdot S^2}"$$

$$‘‘\overline{N^2}"$$

$$W_{ik}=\frac{d_i^2\cdot {\left({\sum }_j{V}_{ij}^T\cdot {\overline{s}}_j\right)}^2}{d_i^2\cdot {\left({\sum }_j{V}_{ij}^T\cdot {\overline{s}}_j\right)}^2+1}\cdot {\delta }_{ik},$$

$${\left\{W\cdot D^{-1}\right\}}_{ij}=\frac{d_i\cdot {\left({\sum }_k{V}_{ik}^T\cdot {\overline{s}}_k\right)}^2}{\left(d_i^2\cdot {\left({\sum }_k{V}_{ik}^T\cdot {\overline{s}}_k\right)}^2+1\right)}\cdot {\delta }_{ij}.$$

$$\overline{M}=R\cdot C^{-1}\cdot C\cdot \overline{s}$$

$$R\cdot C^{-1}=U_C\cdot D_C\cdot V_C^T.$$

$$\hat{s}=C^{-1}\cdot V_C\cdot W_C\cdot V_C^T\cdot C\cdot {\left(R^{T}R\right)}^{-1}\cdot R^T\cdot M.$$

$$\hat{s}=A_C\cdot (R^{T}R{)}^{-1}\cdot R^T\cdot M,$$

$$A_C=C^{-1}\cdot V_C\cdot W_C\cdot V_C^T\cdot C.$$

$$W_{ii}=\frac{d_{Ci}^2\cdot {\left({\sum }_j{V}_{Cij}^T\cdot \left({\sum }_l{C}_{jl}\cdot {\overline{s}}_l\right)\right)}^2}{d_{Ci}^2\cdot {\left({\sum }_j{V}_{Cij}^T\cdot \left({\sum }_l{C}_{jl}\cdot {\overline{s}}_l\right)\right)}^2+1},$$

Full text

1 Introduction

Data unfolding is a common technique used in the analysis of high energy physics (HEP) experimental data. Some of the recent examples in the field of neutrino physics can be found in Refs. [1, 2, 3] and some reviews on this topic can be found in Refs. [4, 5, 6, 7]. The motivation for data unfolding is to estimate the true signal (e.g. energy spectrum) given a measurement that is affected by the detector response as well as statistical (e.g. associated with signal and backgrounds) and systematic uncertainties (e.g. associated with backgrounds, mis-modeling of detector response due to imperfect calibration or finite statistics in simulations). In many applications, data unfolding is not necessarily required. For example, in the case of a hypothesis testing problem, it is generally more advantageous to fold the detector response with the hypothesis and compare with the measurement (See Ref. [8] for more discussions). On the other hand, the data unfolding technique is helpful in many occasions where additional actions are required on the unfolded results. For example, unfolded results are convenient to compare results from different experiments that have different detector responses. Another example would be to extract the ratio of unfolded results (such as cross sections on different nuclei with different detector responses) to be compared with theoretical calculations of ratios, which are typically more precise than the calculation of individual quantities. Finally, the usage of unfolded results, which is generally closer to the true signal than the measurement, has natural advantages for the presentation purpose.

As explained by numerous reviews [4, 5, 6, 7], the main challenge to be overcome in the unfolding of data is the presence of both detector smearing and uncertainties. The random fluctuations due to the existence of statistical and systematic uncertainties would be significantly amplified by a naive inverse of the detector response matrix, which usually leads to meaningless results. This is easy to understand, as the detector smearing represents a loss of information, which in principle cannot be recovered. In HEP, there are two main data unfolding approaches to mitigate this issue. The first method is the Tikhonov regularization (or SVD unfolding) [9, 10]. In this approach, the unfolding problem is expressed as a minimization of a chi-square function comparing the measurement with the prediction. The large fluctuations (also called variance) in the unfolded results are regularized by adding a penalty term into the chi-square function. The penalty term can be chosen to regularize the strength or the curvature (second derivative) of the unfolded results, among other possible choices. A parameter commonly known as the regularization strength can be adjusted freely to control the relative size of the penalty term. A scan of the regularization strength is typically required to obtain the optimum value according to certain pre-chosen metric. A common metric is the summation of variance and bias of the unfolded results. Other choices of the metric can be found in Ref. [4]. The second method is the expectation-maximization iteration with early stopping (or Bayesian unfolding) [11]. In this approach, one would start from an initial guess of the true signal. During each iteration, the guess would be modified according to the difference between the measurement and prediction given the previous guess. Given an initial guess which is non-negative, the solution after an infinite number of iterations approaches the result of minimizing the chi-square under positivity constraints, which require all the unfolded data to be non-negative. This would again suffer from large fluctuations. To mitigate that, the regularization is achieved by stopping the iteration early before convergence. Typically, the number of iterations needs to be scanned to achieve an optimal result. Therefore, an important issue in unfolding is to find an appropriate trade-off between bias and variance of the estimators.

In both unfolding approaches, a scan of the corresponding regularization parameter is required. Inspired by the deconvolution techniques in the digital signal processing, we propose a new unfolding method based on the Wiener filter and the SVD unfolding, which takes into account both the expectation of signal and noise through maximizing the signal to noise ratios in the effective frequency domain with an orthogonal basis and avoids the scanning of any regularization parameter. In Sec. 2, we review the Wiener filter in the digital signal processing employed in Liquid Argon Time Projection Chamber (LArTPC) detectors. We then present the actual Wiener-SVD unfolding algorithm in Sec. 3. In Sec. 4 and Sec. 5, we illustrate the performance of the Wiener-SVD unfolding and compare it with the (Tikhonov) regularization method through two physics examples. The findings are summarized in Sec. 6.

2 Wiener Filter in Digital Signal Processing

The problem of data unfolding shares many common features with the digital signal processing problem, as the goal of both is to extract an estimation of signal from the data. For example, in a LArTPC, the deconvolution technique is used to “remove” the impact of field and electronics response from the measured time-series signal to recover the true signal (the time profile of the number of ionized electrons) [12, 13]. In the following, we briefly review the deconvolution technique.

Deconvolution is a mathematical technique to extract a real signal $S(t)$ from a measured signal $M(t^{\prime})$. The measured signal is modeled as a convolution integral over the real signal $S(t)$ and a given detector response function $R(t,t^{\prime})$ which gives the instantaneous portion of the measured signal at some time $t^{\prime}$ due to an element of real signal at time $t$ in addition to noises $N(t^{\prime})$:

[TABLE]

If the detector response function only depends on the relative time difference between $t$ and $t^{\prime}$,

[TABLE]

we can solve the above equation by doing a Fourier transformation on both sides of the equation:

[TABLE]

where $\omega$ is the frequency and $\hat{S}$ is an estimation of $S$. We can derive the signal in the frequency domain by taking the ratio of the measured signal and the response function:

[TABLE]

When the noise can be ignored, the real signal in the time domain can then be obtained by applying an inverse Fourier transformation to both sides of Eq. 2.4.

When the noise term $N(\omega)$ cannot be neglected, since the response function $R(\omega)$ is typically small at high frequencies due to the shaping of electronics, the noise components in those frequencies will be significantly amplified by the deconvolution $N(\omega)/R(\omega)$ leading to large fluctuations in the deconvoluted signal. Figure 1 shows an example of detector response $R(t^{\prime}-t)$ (left panel), true signal $S(t)$ (middle panel), and simulated measured signal $M(t)$ with noise (right panel).

To address the issue of noise, a filter function $F(\omega)$ is introduced to obtain the estimator of the true signal [12, 13]:

[TABLE]

Its purpose is to attenuate the problematic noise in the deconvolution. The addition of this function can be considered as an augmentation to the response function. A common choice of the filter function is the Wiener filter [14], which is constructed using the expected measured signal $\overline{R^{2}(\omega)\cdot S^{2}(\omega)}:=E[R^{2}(\omega)\cdot S^{2}(\omega)]$ and noise $\overline{N^{2}(\omega)}:=E[N^{2}(\omega)]$ in the frequency domain:

[TABLE]

with $E\left[\cdot\right]$ denotes the expectation operator, which can be understood as the average after a large amount of measurements for the quantity of interest. Note, the expectations of signal and noise squared can be viewed as prior information (e.g. previous measurements). In practice, as we will describe in Sec. 3 (Eq. 3.34), one can also estimate the expectations based on the current measurement. The functional form of Wiener filter in Eq. 2.6 is obtained by minimizing the residual [14, 15]

[TABLE]

The last step omits the $\omega$ and is obtained with the fact $E\left[N\right]\equiv 0$. The minimization of $MSE$ (mean squared error) is achieved through $\frac{\partial MSE}{\partial F}=0$. It is easy to see Eq. 2.6 is recovered.

With the construction in Eq. 2.6, the Wiener filter is expected to achieve the best signal to noise ratio. Besides digital signal processing, the Wiener filter is also widely used in other fields. For example, Wiener filters have been used in experimental astrophysics [16, 17]. Figure 2 shows the constructed Wiener filter in both the frequency and time domains given the example shown in Fig. 1.

With a suitable noise filtering model, an improved estimator for the signal $\hat{S}(t)$ in the time domain can then be found by applying an inverse Fourier transform to $\hat{S}(\omega)$. Essentially, the deconvolution replaces the real field and electronics response function ($R$) with an effective filter response function ($F$ as in the right panel of Fig. 2). For the example shown in Fig. 1, the response function in frequency domain $R(\omega)$ and the measured data in frequency domain $M(\omega)$ are shown in top left and top right panel of Fig. 3, respectively. The deconvoluted results without (left panel with Eq. (2.4)) and with (right panel with Eq. (2.5)) the Wiener filter are shown in the bottom left and the bottom right panel of Fig. 3. Without the (Wiener) filter, the noise in the measured data at high frequency is significantly amplified by dividing the small value of response function, which leads to unacceptable fluctuations in the deconvoluted results. With the Wiener filter applied, the deconvoluted results are comparable to the simulated signal truth. Since the deconvolution problem shares many common features with the data unfolding problem, it is natural to extend the application of Wiener filter technique from deconvolution to unfolding.

3 SVD Unfolding with Wiener Filter

In this section, we describe the procedure of Wiener-SVD unfolding. For clarity, Tab. 1 summarizes the symbols used in this section.

3.1 Problem definition with SVD decomposition

The data unfolding problem generally starts with a $\chi^{2}\left(s\right)$ function defined as

[TABLE]

Here, ${\bf r}$ is an $m$ (row) $\times$ $n$ (column) smearing matrix that connects the vector of measured data ${\bf m}$ (an m-dimensional vector) with an unknown vector of signal $s$ (an n-dimensional vector). Note, we later use $s_{true}$ to represent the true signal in order to differentiate from $s$ which is a variable in calculating $\chi^{2}$. This matrix is general and not limited by functional format in Eq. 2.2. We use $\hat{s}$ to represent the estimator of the true signal $s_{true}$, which is obtained after minimizing the chi-square function. We further restrict ourselves in $m\geq n$ case. The matrix $Cov$ is an $m\times m$ covariance matrix containing all statistical and systematic uncertainties associated with ${\bf m}$ and ${\bf r}$ in calculating the differences between actual measurement ${\bf m}$ and the expectation ${\bf r}\cdot s$. For example, the covariance matrix would include i) statistical uncertainties from data, ii) statistical and systematic uncertainties for the backgrounds, and iii) statistical (with Monte Carlo simulation) and systematic uncertainties associated with the detector response ${\bf r}$. 333In order to evaluate the systematic uncertainties associated with the detector response, one typically runs many Monte Carlo simulations with different detector responses to calcualte the expectated spectra. These spectra are compared with the nominal-detector-response spectrum to construct covariance matrix.

Since the covariance matrix $Cov$ is symmetric, the inverse of it ($Cov^{-1}$) is also symmetric. Hence, $Cov^{-1}$ can be decomposed with Cholesky decomposition into

[TABLE]

where $Q$ is a lower triangular matrix and $Q^{T}$ is its transpose.

Eq. 3.1 can then be rewritten as

[TABLE]

with $M:=Q\cdot{\bf m}$ and $R:=Q\cdot{\bf r}$, and this process is commonly referred to as pre-scaling or pre-whitening. The solution after minimizing Eq. 3.3 would be $M=R\cdot\hat{s}$ or

[TABLE]

In analogy to Eq. 2.3, Eq. 3.4 can be rewritten as

[TABLE]

with $N$ representing the “noise” coming from uncertainties (statistical and systematic uncertainties associated with both ${\bf m}$ and ${\bf r}$). Since $N=M-R\cdot s_{true}=Q\cdot\left({\bf m}-{\bf r}\cdot s_{true}\right)$, each term in the noise vector after pre-scaling follows a normal distribution with $\mu=0$ and $\sigma=1$, since the denominator of the chisquare function in Eq. 3.3 (i.e. square of error) is unity. Given the fact that each term in the noise vector is independent (i.e. uncorrelated), we refer the basis in this domain to be orthogonal.

Using the singular value decomposition (SVD) approach, $R$ can be decomposed as

[TABLE]

with both $U$ ($m\times m$) and $V$ ($n\times n$) being orthogonal matrices that satisfy $U^{T}\cdot U=U\cdot U^{T}=I_{(m\times m)}$ and $V^{T}\cdot V=V\cdot V^{T}=I_{(n\times n)}$. $I$ is the identity matrix and the subscript represents the dimension. $D$ is an $m\times n$ diagonal matrix with non-negative diagonal elements (known as singular values) $D_{ii}$ = $d_{i}$ arranged in descending order as $i$ increases.

Inserting Eq. 3.6 into Eq. 3.5, we have

[TABLE]

where $R_{U}:=U^{T}\cdot R$, $N_{U}:=U^{T}\cdot N$, and $M_{U}:=U^{T}\cdot M$ are transformations of the smearing matrix, the noise $N$, and the measured signal, respectively. This equation can be compared to Eq. 2.4. Note, since $U$ is an orthogonal matrix and elements of the original noise vector $N$ are uncorrelated, elements of the new noise vector $N_{U}$ are still uncorrelated. Each element follows a normal distribution with $\mu=0$ and $\sigma=1$. Thus, the basis in this new domain is still orthogonal.

Given Eq. 3.7, we can understand the large fluctuation in the unfolded results. First, due to the existence of the smearing in matrix ${\bf r}$ or $R$, the magnitude of singular values $d_{i}$ drops significantly as $i$ increases. After a certain $i$, the value of $d_{i}$ can be extremely small which leads to a gigantic value in the corresponding element in $D^{-1}$. In the case of perfect signal without any noise $N_{U}$, these gigantic diagonal elements are effectively canceled out by the small values in the signal $R\cdot s$ leading to recover the signal without any bias after data unfolding. The situation is completely changed in the presence of noise $N_{U}$. Since these noise will be significantly amplified after multiplying with $D^{-1}$, the unfolded results suffer from large fluctuations. From the above discussions, it is easy to see the similarities between Eq. 3.7 and Eq. 2.4 for the deconvolution discussion in Sec. 1. Therefore, in analogy to the results after the Fast Fourier Transformation (FFT), we refer to $M_{U}$ after the SVD transformation as the measurement in the effective frequency domain in analogy to the frequency domain in the signal processing (e.g. Eq. 2.3 and Eq. 2.4). Both these domains have orthogonal basis that are linear transformation from the original basis. While FFT requires the functional format of response function to be symmetric (i.e. Eq. 2.2), the SVD does not require this symmetric condition and is more general.

3.2 Review of traditional regularization approach

Regularization is a commonly used technique to address the problem described in the previous section. It imposes additional constraints on the estimation of true distribution $s$ by introducing a regularization function $\Sigma(s)$. The estimator can be obtained by finding the maximum of a weighted combination of log-likelihood $\log{L}$ and $\Sigma$:

[TABLE]

where $\tau$ is called regularization strength, which determines the trade-off between bias due to imposed constraints and variance due to existence of the noise in the unfolded distribution. In general, to obtain the best estimation of the signal, the log likelihood $\log{L(s)}$ as well as the regularization function $\Sigma(s)$ are required to be sufficiently well-behaved (e.g. at the very least there should not contain multiple local maxima [18]).

The so-called Tikhonov regularization technique uses the following regularization function

[TABLE]

with the spectrum $s$ depending on the variable $E$ (e.g. energy). For example, when $k=0$, we have $\Sigma^{0}\left(s(E)\right)=-\int\left(s(E)\right)^{2}\text{d}E$, which favors small values of the signal. Another commonly used example is $k=2$, in which $\Sigma^{2}\left(s(E)\right)$ represents a measure of the average curvature of distribution $s(E)$, imposing a constraint on the smoothness of the signal. When applying the regularization technique, one generally chooses a regularization function, evaluates the bias and variance of the estimator as a function of regularization strength $\tau$. The value of $\tau$ is then optimized based on a predetermined metric.

In the case of $\Sigma^{0}\left(s(E)\right)$ regularization, the unfolded results can be expressed as:

[TABLE]

Here $A$ behaves as an additional smearing matrix added to the unfolding results of Eq.3.7. It has the form of:

[TABLE]

where $F$ is an $n\times n$ diagonal matrix with elements satisfying

[TABLE]

here $\tau$ is the regularization strength. Eq. 3.10 is then changed to

[TABLE]

At small values of $d_{i}$, the corresponding noise terms in $N_{U}$ are now suppressed by $d_{i}/\left(d_{i}^{2}+\tau\right)$ at finite value of $\tau$ instead of being amplified by $1/d_{i}$. Such a change would suppress the large fluctuations in the unfolded results due to the presence of the noise $N_{U}$. It’s worth noting that the regularization method is effectively introducing an additional smearing to the unfolding results, which would lead to biases on the unfolded signal. The above derivation is similar for other $\Sigma^{k}\left(s(E)\right)$ regularization schemes.

3.3 Wiener-SVD approach

In the $\Sigma^{0}\left(s(E)\right)$ regularization, the functional form of $F_{ii}=d_{i}^{2}/\left(d_{i}^{2}+\tau\right)$ is independent of the signal shape, and the choice of $\tau$ is obtained through a scan of this parameter with respect to some metrics. With the concept of the Wiener filter, we can construct $W$ (replacing $F$ in regularization) directly to optimize the signal to noise ratio. This replacement defines the Wiener-SVD approach, in which the functional form of $W$ 444We replace $F$ by $W$ for Wiener filter. also considers the expectation value of the signal in the effective frequency domain: 555In general, the $s_{true}$ is unknown, so the expectation of signal $\overline{s}$ is used.

[TABLE]

Here, the expectation of signal is assumed to be known. We will come back to this point later in this section on how to obtain the expectation of signal. The construction of $W$ is based on the Wiener filter ${\rm``}\overline{R^{2}\cdot S^{2}}{\rm"}/\left({\rm``}\overline{R^{2}\cdot S^{2}}{\rm"}+{\rm``}\overline{N^{2}}{\rm"}\right)$ as in Eq. 2.6. Taking Eq. 3.14, at bin $i$, we have

[TABLE]

resulting in a Wiener filter of

[TABLE]

replacing $F$ in Eq. 3.11. Here, Eq. 3.16 is obtained, since each element of noise $N_{U}$ follows a normal distribution with $\mu=0$ and $\sigma=1$. We have

[TABLE]

The small value of $d_{i}$ is balanced by the finite value of the expectation value of $\overline{N^{2}}\equiv 1$. From Eq. 3.17, the construction of the Wiener filter takes into account both the strengths of signal and noise expectations and is free from regularization strength $\tau$.

A few comments should be made regarding the Wiener-SVD approach:

• •

Generalized Wiener-SVD approach:

As shown in Ref. [9], the regularization can be applied on the curvature of the spectrum instead of the strength of the spectrum. This involves an additional matrix $C_{2}$. This can also be achieved in the Wiener-SVD approach:

[TABLE]

by including an additional matrix $C$ that has the commonly used regularization forms, such as the first and second derivatives. Since the effective frequency domain is determined by the smearing matrix $R$, the inclusion of $C$ would alter the basis of the effective frequency domain. In this case, the SVD decomposition becomes

[TABLE]

The final solution of the regularization would become

[TABLE]

or

[TABLE]

where

[TABLE]

The corresponding Wiener filter would be

[TABLE]

where $C_{jl}$, $V^{T}_{Cij}$, and $d_{Ci}$ are matrix elements of matrices $C$, $V_{C}$, and $D_{C}$, respectively.

• •

Covariance matrix of unfolded results:

Since the unfolded results are a linear transformation of the measurement, we can easily evaluate the uncertainties associated with them. Eq. 3.22 can be rewritten into

[TABLE]

with

[TABLE]

Then, the covariance matrix of $\hat{s}$ can be deduced from the $Cov$ (the covariance matrix of $M$) as

[TABLE]

• •

Variance:

The variances of the unfolded data can also be easily calculated given that their origin $N$ in Eq. 3.21 is well understood. Defining $N\left(i\right)$ as a vector with the $i$th element being 1 and the rest of elements being 0, we can calculate the variance in $s$ due to $i$th element in $N$ as:

[TABLE]

with $T_{deviation}\left(i\right)$ being a vector. The variance of the $j$th element of $T_{deviation~{}j}$ can thus be written as:

[TABLE]

after summing the contribution from each independent noise source. The square of $T_{deviation~{}j}$ corresponds to the $j$th diagonal element of covariance matrix in Eq. 3.27.

• •

Bias:

Given Eq. 3.22, we can understand the entire process of unfolding as to "remove" the effect of $R$ through multiplying $\left(R^{T}R\right)^{-1}\cdot R^{T}$ and then replace it with a new smearing matrix $A_{C}$. Therefore, it is straightforward to estimate the bias on the unfolded results:

[TABLE]

with $I$ being identity matrix and $\overline{s}$ being the expectation of the signal.

Given only the measurement m, an alternative approach given in Ref. [4] can be used to estimate the bias. The bias for bin $i$ now is defined as:

[TABLE]

where $\hat{M}=R\cdot\hat{s}$. Using Eq. 3.22 and  3.31, we have

[TABLE]

As can be seen, in Eq. • ‣ 3.3 the bias is directly calculated if $\overline{s}$ is known, while for Eq. 3.32 the bias is estimated by replacing the $\overline{s}$ in Eq. • ‣ 3.3 with $\hat{s}$. For Tikhonov regularization, the bias estimation of Eq. 3.32 deviates from that of Eq. • ‣ 3.3 at large $\tau$ values when the unfolded spectrum significantly deviates from the true spectrum. At very small values of $\tau$, the unfolded spectrum suffers from large fluctuations leading to a significant overestimation of bias with Eq. 3.32. Moreover, since $M$ has fluctuations, the variance of $T_{bias}^{\prime}$ can also be calculated. The variance in $T_{bias}^{\prime}$ due to $i$th element in $N_{U}$ has the form of:

[TABLE]

• •

Expectation of signal:

We should note that Eq. 3.15 only considers one particular model of signal $s$. In reality, the expectation of signal should cover a range of possible signals $\overline{s(k)}$ (e.g. prior models) that are compatible with existing observations

[TABLE]

Here, $\chi_{k}^{2}$ is the chi-square representing the compatibility between the prediction $s(k)$ and the measurement. When there is no prior models available, one can construct these models using general functions (e.g. Legendre polynomials or spline functions).

• •

Regularization interpretation of Wiener-SVD approach:

We show that the Wiener-SVD unfolding method is equivalent to a regularization which attempts to maximize the signal to noise ratio in the effective frequency domain $M_{U~{}i}=D_{i}\cdot\sum_{j}V^{\intercal}_{ij}s_{j}$. Recall Eq. 3.8, one now has:

[TABLE]

with the expectation of noise square being $1$ in the effective frequency domain. Using the procedure detailed by [18], by maximizing $\phi(s)$ one obtains the following estimator

[TABLE]

where

[TABLE]

With $X$ and $Y$ evaluated at the expectation of $\overline{s}$ and $\overline{M}$, Eq.3.36 can be rewritten as

[TABLE]

where

[TABLE]

Therefore, we have

[TABLE]

with

[TABLE]

One recovers Eq.3.17.

4 Data Unfolding Example: Cross Section Extraction

In this example, we apply the Wiener-SVD unfolding on a neutrino cross section extraction problem. As introduced in Sec. 1, the data unfolding technique can be useful for this problem when the ratio of cross sections is desired or when the comparison of cross section measurements from different experiments is needed.

Experiments that engage in neutrino cross-section measurements generally consist of two parts: a neutrino beam produced by bombarding a target with a proton beam, and a detector (or a series of detectors) located a few hundred meters away from the target to detect neutrino interactions. The neutrino beam composition and energy distribution are generally well-understood. Depending on the detector technology, the neutrino energy can be reconstructed via calorimetry or from the kinematics of final state particles. Due to the smallness of neutrino cross-sections, the signal statistics are typically low in such measurements. Depending on the beam configuration (on-axis or off-axis) one can have a broad or narrow neutrino energy spectrum. Some neutral hadrons produced by neutrino interactions, neutrons in particular, could leave undetectable for calorimetric or tracking detectors, therefore the reconstructed visible energy tends to be smaller than the true neutrino energy. For simplicity, we neglect the neutrino flux uncertainties and only consider the reconstructed neutrino energy spectrum $M(E_{r})$. Figure 4 shows the true energy spectrum $s_{true}$ , detector smearing matrix ${\bf r}$, Asimov spectrum [19] $\overline{\textbf{m}}={\textbf{r}}\cdot s_{true}$, measured spectrum m which is the Asimov spectrum with random Poisson statistical fluctuation, and covariance matrix with statistical uncertainty only. A Gaussian true spectrum is assumed; detector smearing matrix is mocked up such that reconstructed energy is skewed towards energy lower than the true neutrino energy. No systematic uncertainty or background is considered in this toy experiment.

In order to illustrate the performance of the Wiener-SVD method, we compare the unfolded results with those from the Tikhonov regularization described in Sec. 3.2. Three choices of $C$ matrices are used for comparison. They are:

[TABLE]

which correspond to the k=0, k=1 (first-order derivative), and k=2 (second-order derivative, or curvature) cases in Eq. 3.9, respectively. Since $C^{-1}$ is needed to construct the Wiener filter $W$, to make $C$ invertible, a very small value of $\epsilon=10^{-8}$ is added to the diagonal elements of $C_{2}$ matrix.

In addition, it should be noted that one has the freedom to normalize the unfolded distribution to that of the measured distribution. This is equivalent to imposing a constraint on the total number of events in Eq. 3.8:

[TABLE]

where $\lambda$ is a Lagrange multiplier, and $\partial\phi/\partial\lambda=0$. This normalization can be of particular importance for the $C_{0}$ case with low statistics or large systematic uncertainties.

Given a choice of the matrix $C$, the unfolded results with the Wiener-SVD method can be obtained through Eq. 3.21. For Tikhonov regularization, the unfolded results also depend on the regularization strength $\tau$. For the results shown in this section, the optimal regularization strength is determined by minimizing the following Mean Square Error (MSE) [4]:

[TABLE]

where $\sigma^{2}=\sum T_{deviation~{}i}^{2}$ is the total variance and $b^{2}=\sum T_{bias~{}i}^{2}$ is the total bias square. Here $i$ represents the $i$th bin and the definition of $T_{deviation}$ and $T_{bias}$ can be found in Eq. 3.28 and Eq. • ‣ 3.3, respectively.

Figure 5 shows the Wiener filter and regularization filter in the effective frequency domain for the $C_{0}$ and $C_{2}$ cases. To construct the Wiener filter in this toy example, the expectation of signal $\overline{s}$ is taken to be the true signal $s_{true}$. We use Eq. 3.34 to construct the Wiener filter in the example described in next section. The Wiener and regularization filters assign different weights to each effective frequency bin: regularization is scaled by the eigenvalue corresponding to each bin, whereas the Wiener filter is scaled by the “signal/noise-weighted”eigenvalue. Both the Wiener and regularization filters suppress high frequency bins, and therefore reduce the impact from random fluctuation at small $d_{i}$ values. In order to emphasize the importance of normalization, the left panel of Fig. 5 assumes that the total systematic uncertainty is $\sqrt{10}$ times as the statistical uncertainty for the $C_{0}$ case. When the statistics are low, or the systematic uncertainty is high, the Wiener and regularization filters both yield greater suppression. Since the filter is multiplied on the measurement, the normalization is often needed to further reduce the bias, especially when the $C_{0}$ matrix is used.

Figure 6 shows the unfolded results based on the $C_{2}$ case for the regularization method (left panels) and the Wiener-SVD method (right panels). The unfolded spectra (top panel) and residual (the third panel to top) are similar between the two methods. The additional smearing matrix (defined by Eq. 3.11 and shown in the second panel to top) from regularization is more local than that of Wiener-SVD. This is straightforward to understand, as the $C_{2}$ regularization constrains on the smoothness and Wiener-SVD constrains on the signal to noise ratio in the effective frequency domain. The bottom panels of Fig. 6 show the covariance matrices of the unfolded results. In comparison to the diagonal covariance matrix of the measurement, the unfolded covariance matrix is no longer diagonal due to the application of an additional smearing matrix.

Figure 7 shows the quantitative comparisons of the results from Tikhonov regularization and Wiener-SVD unfolding in this example. In the left panel, the bias squared v.s. variance is plotted. For the Tikhonov regularization method, the regularization strength $\tau$ is scanned from 0 to 1. As shown, at fixed variance (bias), the bias (variance) of the Wiener-SVD result is smaller than those of the Tikhonov regularization method. For both methods, the variance and bias of unfolded results with $C_{2}$ applied are better than those of $C_{1}$ and $C_{0}$. The right panel of Fig. 7 shows the MSE as a function of the regularization strength. We see that the MSEs of Wiener-SVD are smaller than the corresponding ones of Tikhonov regularization.

5 Data Unfolding Example: Reactor Neutrino Flux

In the previous section, we illustrated the Wiener-SVD method with a low-statistics neutrino cross section extraction example, in which the problem is simplified by using a diagonal covariance matrix with only statistical uncertainties. In this section, we show a high-statistics example by constructing a toy reactor neutrino experiment to extract the reactor antineutrino energy spectrum from the measured visible energy spectrum using the Wiener-SVD approach. In particular, we will implement a more realistic covariance matrix including both statistical and systematic uncertainties and illustrate how to construct the Wiener filter with a group of theoretically well-motivated models.

In a typical reactor antineutrino experiment such as the Daya Bay experiment [20], the antineutrinos are detected through the inverse beta decay (IBD) process $\bar{\nu}_{e}+p\to e^{+}+n$. The positron gives a prompt signal including its kinetic energy and the two 511 keV annihilation gamma-rays, whereas the neutron after thermalization gets captured in the detector and yields a delayed signal. Since the energy carried away by the recoil neutron is small, the neutrino energy $E_{\nu}$ can be approximately calculated from the prompt energy $E_{p}$ by $E_{p}\simeq E_{\nu}$ - 0.8 MeV. The measurement of $E_{p}$ can be affected by a variety of systematic effects. For instance, in the Daya Bay experiment where the liquid scintillator is used as a calorimeter to determine the particle energy, the response from a particle’s true energy to its visible energy is nonlinear. The nonlinearity is caused by both the quenching effect of the scintillator and the additional photons produced by the Cerenkov radiation. In addition, particles could lose energy in the non-scintillating materials, which further alters the visible energy. Various electronics nonlinear response can also occur and impact the total visible energy. The resolution of $E_{p}$, typically $\sim$8%, is mainly determined by the fluctuation of photoelectrons that follows the Poisson distribution. The gain variation, dark noise, and detector non-uniformity further add to the energy resolution. In order to construct the detector energy response matrix and the associated uncertainties, typically a comprehensive detector calibration campaign and data-Monte-Carlo comparison is necessary to fully understand these detector effects.

We generated a $\sim$50k events toy reactor neutrino experiment using the Huber and Mueller reactor models[21, 22] with a typical commercial reactor fission fractions for the four main isotopes: 235U, 238U, 239Pu, and 241Pu. We used the detector energy response reported by the Daya Bay experiment [23] to resemble a realistic situation. The covariance matrix used for the “measured” prompt spectrum in this toy study is generated based on the covariance matrix given in Ref. [23] with the corresponding statistics. Figure 8 shows the inputs used for this toy study. The left panel shows the “measured” prompt spectrum (blue) that includes the detector smearing effect and fluctuation due to uncertainties, the true neutrino spectrum (black), and the detector smearing matrix (in the inset). The smearing matrix includes all detector response effects as mentioned above. The right plot of Fig. 8 shows the covariance matrix for the prompt spectrum, which includes both the statistical and systematic uncertainties.

As discussed in Sec. 3, the first step of unfolding using SVD method is to do a pre-scaling to normalize and remove correlations of uncertainties among bins. Figure 9 shows the prompt spectrum (left plot) and smearing matrix (right plot) after the pre-scaling. As can be seen, both the pre-scaled prompt spectrum and the smearing matrix are quite different from the original ones (i.e. Fig. 8).

In practice, the “true” model is always unknown, so it is not directly available for the purpose of constructing the Wiener filter $W$. Instead, the Wiener filter can be constructed through a group of theoretical models using Eq. 3.34. In this example, we consider a variety of reactor flux models generated from the linear combinations of the calculations in Ref. [24, 21, 22]. The $\chi^{2}_{k}$ for model $s_{k}$ is then constructed by comparing the spectra of the prediction ${\bf m_{k}}$ $:=$ ${\bf r}\cdot s_{k}$ and the measurement ${\bf m}$.

[TABLE]

where $M_{k}$ is ${\bf m_{k}}$ after pre-scaling: $M_{k}$ $:=$ $Q\cdot{\bf m_{k}}$.

Figure 10 compares the Wiener filters $W$ and Tikhonov regularization filters $F$ in the effective frequency domain. They are constructed with the $C_{0}$ and $C_{2}$ matrices. For the regularization filters, the value of the regularization strength $\tau$ is chosen by minimizing the MSE defined in Eq. 4.3. As discussed previously, since the regularization filters only consider the $R$ in the effective frequency domain, the larger the “frequency” the more the suppression. As can be seen in both panels of Fig. 10, all of them behave as monotonically decreasing functions. On the other hand, the Wiener-SVD method considers not only $R$ but also the signal to noise ratio of each bin in the effective frequency domain. Therefore, the shapes of the Wiener filters do not necessarily behave monotonically decreasing. For instance, The Wiener filter constructed with the $C_{0}$ (left panel of Fig. 10) has very large suppressions for the two medium frequency bins (bin 6 and 7). Compared with the Wiener filter, the monotonicity of the regularization filter dictates that it inevitably will keep more noise at these medium frequency bins and remove more signal at some of higher frequency bins even when they are not very noisy. For filters constructed with $C_{2}$ (right panel of Fig. 10), the Wiener filter has a similar shape as the regularization filter. Nevertheless, the details of suppression at high frequency bins are still different.

The unfolded results of the Wiener-SVD and regularization methods can be seen in Fig. 11. The regularization unfolding uses the regularization strength $\tau$ = 2.4$\times$10*-5* through minimizing the MSE in Eq. 4.3. For simplicity, only results with the using of a $C_{2}$ matrix are shown. Both methods produce reasonable unfolded results. To compare the Wiener-SVD and regularization unfolded results with different values of $\tau$, we plot the total variance versus total bias square in the left panel of Fig. 12. The bias is calculated using Eq. • ‣ 3.3 with $\overline{s}$ set to the model that has the smallest $\chi^{2}$ value (see Eq. 5.1). The black curve is from the regularization method with a wide range of $\tau$. The red square is from the Wiener-SVD method described previously. Similar to the cross section example in the previous section, at the same variance (bias), the Wiener-SVD method has a smaller bias (variance). The right panel of Fig. 12 shows the corresponding MSE values. Again, similar to the previous section, the Wiener-SVD unfolded result has a smaller MSE than any unfolded result from the regularization method. To illustrate the necessity of using $\bar{s}$ to represent the unknown $s_{true}$, we also show an unfolded result (blue triangle) with a Wiener filter constructed using an improper expectation (quite different from the $\bar{s}$). In this case, the MSE of the unfolded results with the Wiener filter is no longer the smallest. In practice, it is crucial to use Eq. 3.34 to evaluate the signal expectation in order to achieve the optimal bias and variance in the Wiener-SVD approach. This is the common issue for all kinds of unfolding approaches, and in fact the best result of regularization method (minimum point in Fig. 12 right panel) would also be altered by comparing the unfolded result with the improper signal expectation.

6 Discussions and Recommendations

Based on the examples in previous two sections, we make the following recommendations regarding data unfolding:

• •

The SVD-based unfolding methods (traditional regularization filter or Wiener filter) are equivalent to replacing the detector smearing matrix with a new smearing matrix $A_{C}$. The application of this new smearing matrix is crucial to suppress the large oscillation (high variance) of the direct matrix inversion unfolded results. In evaluating bias, the $A_{C}$ is applied to the true spectrum, which leads to bias.

• •

We recommend to report this new smearing matrix $A_{C}$ in the publication together with the unfolded results to enable a more direct comparison of expectations (e.g. from new theoretical calculations) with the unfolded results. In practice, the new smearing matrix should be applied to the theoretical calculation before comparing to the unfolded results. Through reporting the new smearing matrix $A_{C}$, one can avoid including the bias due to unfolding in the final uncertainties.

• •

For SVD-based approach, the $C_{2}$ typically yields a better result than those of $C_{1}$ and $C_{0}$.

• •

The covariance matrix, variance, and the bias of the unfolded results can be calculated with Eq. 3.27, Eq. 3.28, and Eq. • ‣ 3.3, respectively.

• •

The Wiener filter should be constructed by Eq. 3.34, which takes into account the measurement given a range of prior expectations.

From these two toy examples shown in Sec. 4 and Sec. 5, we can conclude the following pros and potential cons for the Wiener-SVD approach in comparison to the Tikhonov regularization approach:

• •

The Wiener-SVD approach is free from a regularization parameter, which is required to be optimized for the regularization approach. This is achieved by utilizing the expectations of signal and noise, which can be viewed as a direct determination of an effective regularization parameter. At a fixed variance (bias), unfolded results from the Wiener-SVD method gives a better bias (variance) than the regularization method. This presumably is due to the optimized signal to noise ratio in the effective frequency domain for the Wiener filter. When evaluated with the MSE (metric defined in Eq. 4.3), the unfolded results based on the Wiener filter is comparable and sometimes better than the best result from the regularization approach.

• •

For $C_{1}$ and $C_{2}$, the traditional regularization method pulls the estimator towards smoothness, even though the true distribution is not necessarily so. The Wiener filter considers the signal to noise ratio in the effective frequency domain, which does not require smoothness and is more general.

• •

To construct the Wiener filter, an estimation of the true model is required. The unfolded results in the Wiener-SVD method is strictly model dependent. Such dependence is reduced when the estimation of the true model takes into account the actual measurement as illustrated in Eq. 3.34. In addition, the different choices of the true model only affects the construction of the new smearing matrix. By reporting the new smearing matrix, the model dependence of the results can be avoided, as the new smearing matrix can be applied to the other expectations to be tested.

• •

As shown in Fig. 6, the new smearing matrix of the Wiener-SVD is less localized than those of the regularization method. This could be a potential disadvantage of the Wiener-SVD approach, but can again be mitigated by reporting the new smearing matrix.

The derivation of the Wiener filter construction as shown in Eq. 3.24 is based on the covariance matrix and SVD decomposition of the smearing matrix $R\cdot C^{-1}$. In the case when the uncertainties cannot be simply expressed via a covariance matrix (i.e. Gaussian approximation), the Wiener filter can still be constructed if the smearing matrix can be constructed. In this case, the expectation of signal square can still be constructed using Eq. 3.34. The expectation of noise square would be obtained through a Monte Carlo approach.

7 Summary

Inspired by the deconvolution technique employed in the digital signal processing, we introduce a new unfolding technique based on the Wiener filter and SVD technique for HEP data analysis. Through maximizing the signal to noise ratios in the effective frequency domain, the Wiener-SVD unfolding avoids the scanning of any regularization parameter in the traditional approaches. Through a couple examples, we show that the unfolded results from the Wiener-SVD method generally have a smaller bias (variance) at fixed variance (bias) than the unfolded results from the Tikhonov regularization method. The overall MSE averaging the total bias and variance is also generally smaller for the Wiener-SVD method. These features support the Wiener-SVD method as an attractive option for the data unfolding problem. An implementation of the Wiener-SVD method can be found in Ref. [25].

Acknowledgments

We thank Tom Junk, Clark McGrew for fruitful discussions and David Caratelli for carefully reading of the manuscript. This work is supported by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, and Early Career Research Program under contract number DE-SC0012704.