{'fraction_non_alphanumeric': 0.030303030303030304, 'fraction_numerical': 0.007767691821343088, 'mean_word_length': 4.522758448684375, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
SVD Approach to Data Unfolding

Distributions measured in high energy physics experiments are usually distorted and/or transformed by various detector effects. A regularization method for unfolding these distributions is re-formulated in terms of the Singular Value Decomposition (SVD) of the response matrix. A relatively simple, yet quite efficient unfolding procedure is explained in detail. The concise linear algorithm results in a straightforward implementation with full error propagation, including the complete covariance matrix and its inverse. Several improvements upon widely used procedures are proposed, and recommendations are given how to simplify the task by the proper choice of the matrix. Ways of determining the optimal value of the regularization parameter are suggested and discussed, and several examples illustrating the use of the method are presented.


Introduction
In high energy physics, measurements of physical observables -spectra of invariant masses, angular distributions of particles, etc. -are usually distorted and transformed by various effects such as finite resolution and limited acceptance of the detector. For this reason, it is often impossible, or at least very difficult, to make direct comparisons of the data obtained using different detectors with each other and with various theoretical predictions.
In order to overcome this problem, some sophisticated studies of the measurement process are usually carried out. Ideally, these studies should result in a two-variable function describing the response of the detector, so that the actual measured distribution can be considered as a convolution of this function with the true one. This in general leads to an integral equation for the true distribution. Solving this equation (i.e. unfolding the true distribution) usually requires some kind of discretization, leading to a system of linear equations. The problem, however, belongs to a class of ill-posed problems, which are unstable against small variations in the initial system. Because of inevitable statistical errors in the measured distribution, the exact solution (if it exists) is usually wildly oscillating and useless (see ref. [1] for a good introduction to the subject 1 ).
For a number of reasons, in present high energy physics applications, the above approach is usually replaced by a discrete Monte Carlo simulation of the measurement process, resulting directly in a system of linear equations for the underlying true discrete distribution. In this case all the above difficulties are aggravated by statistical and possibly systematic errors in the response matrix itself.
In order to avoid these difficulties, it is sometimes advisable to fold the theoretically predicted true distribution with the estimated response matrix, and compare the thus folded theoretical spectrum with the measured one. This method is stable and may be useful in certain cases, but is clearly useless if a comparison between various experiments is required, or if the functional form of the distribution is unknown.
The problem of unfolding has been studied in various forms, giving rise to a number of independent methods described in the literature. For instance, a method based on Bayes' theorem was proposed in [2,3]. The authors manage to avoid partly the inversion difficulties by using a non-linear iterative procedure, leading asymptotically to the unfolded distribution.
Another way of overcoming the instability of unfolding is to use some kind of regularization condition, based on some a priori information about the solution. One can demand, for example, that the true solution has minimum curvature (i.e. is quite smooth [1]), or that it is strictly positive [4]. These methods usually allow the suppression of spurious oscillating components of the unfolded solution and often lead to satisfactory results, though concrete implementations happen to be quite lengthy and complicated.
In this paper, we propose a new way of analyzing unfolding problems, based on the Singular Value Decomposition (SVD) of the response matrix. A simple and transparent 2-dimensional example is used to explain and illustrate the reasons for the apparent instability of the problem, and a discrete analog of the minimum curvature condition is used to stabilize the unfolded solution. An unfolding procedure based on considerations which are similar to those presented here has been described more than a decade ago in [1], and is still widely used in experimental analyses. While there are no significant differences between the foundations of the two procedures, we believe that our formulation allows one to obtain more reliable and precise results, being at the same time much simpler and easier to implement.
Our approach is based on the extensive use of SVD, and results in a linear unfolding algorithm which is applicable to a wide range of problems. We derive also a number of important recommendations about the proper normalization of the response matrix and choice of the variables. Complete error propagation is implemented, and a vivid and reliable procedure for determining the optimal value of the regularization parameter is suggested.
The paper is organized as follows. Sect. 2 contains our conventions for the notation used throughout the paper. The problem is formulated in Sect. 3. Sect. 4 is devoted to the Singular Value Decomposition, illustrated by a simple example of a 2 × 2 matrix. The benefits and pitfalls of rescaling equations and unknowns are analysed in Sect. 5. Regularization of the system and actual unfolding is performed in Sect. 6, and the choice of the regularization parameter is discussed in Sect. 7. The step-by-step unfolding algorithm is presented in Sect. 8 and is illustrated by two distinctly different examples in Sect. 9. Conclusions are drawn out in Sect. 10.

Notation
We will adopt the following notation conventions, which will be used throughout the paper: • All one-dimensional histograms/vectors are denoted by small letters (e.g. b, z etc.).
• All two-dimensional histograms/matrices are denoted by capital letters (e.g. A, X etc.).
• The covariance matrix associated with a one-dimensional variable is denoted by the same letter in capital (i.e. the matrix W ij denotes the covariance matrix of the vector w j ).
• No implicit summation is assumed over repeated indices; i.e. no repeated index is summed unless the summation is explicitly shown.
• Upper index T stands for the transposed matrix, A T ik = A ki , so that the euclidean norm of a vector z equals √ z T z.
• Upper index −1 denotes the inverse matrix, A −1 A = AA −1 = I, where I stands for the unit matrix, I ik = δ ik .

The Problem
Let the distribution of a measured observable be stored in a vector b of dimension n b , where the ith coordinate of the vector contains the number of entries in the corresponding bin of the histogram. The measurement is affected by the finite experimental resolution and/or the limited acceptance of the detector, so that each event from the true distribution may find itself in a range of (not necessarily) adjacent bins, or nowhere at all. Let us assume, that we are able to simulate the measurement procedure of this observable (e.g. with Monte Carlo techniques). We generate the distribution x ini of dimension n x , according to some idea of the underlying physical process, and perform our detector simulation. At this stage, every entry in a measured bin (i.e. every event) can be directly traced to its origin. This gives us a well defined system of linear relations between the simulated true and measured distributions: The n b × n x matrixÂ is a probability matrix, which actually performs the simulated folding procedure. Now, withÂ and x ini given, for any vector b obtained by a real measurement using the detector described by its response matrixÂ, one can attempt to find a corresponding unfolded true distribution x. It is well known that trying to solve the linear system of equationŝ against x directly, using the exact inversion of the matrix, usually leads to completely unacceptable rapidly oscillating solutions. In the following, based on the Singular Value Decomposition (SVD) of the response matrixÂ, we analyse reasons for this behavior, locate the difficulty and propose a relatively simple and straightforward regularization method, which allows one to suppress spurious, quickly oscillating components of the solution, leaving only statistically significant terms. The above formulation of the problem may not seem general enough, but it clearly incorporates most cases of practical interest, is easily understood and interpreted, and is well-suited for the unfolding method described in this paper. Let us, however, express the discrete distributions x, b and the response matrixÂ in terms of the underlying continuous probability density functions.
Let y true be the continuous true variable under consideration, whose variation range {y true 0 ÷ y true nx } is divided into n x bins with boundaries y true j , j = 1, . . . , n x − 1. Each component of the vector x is then calculated as an integral over the true distribution function X (y true ) in the appropriate range: Analogously, let {y 0 ÷ y n b } be the variation range of the measured variable y, with bin boundaries y i , i = 1, . . . , n b − 1. Then the components of the vector b are appropriate integrals over the continuous distribution function B(y) : LetÂ(y, y true ) be the detector response function, which maps the true distribution to the observed one, according to the convolution integral: After this, the response matrixÂ can be defined as the ratio of two integrals: Each elementÂ ij equals to the probability for an event generated in the true bin j to be found in measured bin i. In high energy physics applications the response functionÂ(y, y true ) is usually not known analytically. Instead, some sophisticated detector simulation techniques are used to determine the matrix directly as explained above.

Definitions
A singular value decomposition (SVD) of a real m × n matrix A is its factorization of the form where U is an m × m orthogonal matrix, V is an n × n orthogonal matrix, while S is an m × n diagonal matrix with non-negative diagonal elements: The quantities s i are called singular values of the matrix A, and columns of U and V are called the left and right singular vectors.
The singular values contain very valuable information about the properties of the matrix. If, for example, A is itself orthogonal, all its singular values are equal to 1. On the contrary, a degenerate matrix will have at least one zero among its singular values. In fact, the rank of a matrix is the number of its non-zero singular values. If the matrix and/or the r.h.s. of a linear system is known with some level of uncertainty, and some singular values of the matrix are significantly smaller than others, the system may be difficult to solve even if formally the matrix has full rank. In many aspects such matrices behave like degenerate ones, and SVD suggests a method of treating such problems, which is common for small and exactly zero singular values.
We will assume that the singular values s i form a non-increasing sequence. This is easily achieved by swapping pairs of singular values, swapping simultaneously corresponding columns of U and V . We will assume also that m ≥ n, which means that the number of bins in the measured histogram b should not be smaller than the number of bins in the unfolded histogram x. If necessary, one can just add rows of zeroes to the initial matrix.
Comprehensive descriptions of SVD with many technical details and examples can be found in the literature (see, e.g., the excellent books [5,6]). One of the most attractive features of this procedure is that one does not really have to perform SVD by hand. A very efficient and transparent FORTRAN subroutine (called, not surprisingly, SVD) is present in the CERN program library. Some earlier implementations can be found in refs. [5,6] as well.
Once the matrix is decomposed into the form (7), its properties can be readily analyzed and it becomes very easy to manipulate, as illustrated in following subsections. This kind of analysis is extremely useful for ill-defined linear systems with almost (or even exactly) degenerate matrices, as it not only locates the difficulty, but can also suggest ways of overcoming it.

A simple example
Consider a very simple 2 × 2 example, which nevertheless incorporates most of the interesting aspects of the problem.
Let the response matrixÂ have the form with 0 ≤ ǫ ≤ 1 determining the "quality" of the detector: ǫ = 1 means an ideal detector with the response matrix equal to unity, while small ǫ ≪ 1 corresponds to a poor detector, almost unable to distinguish the two bins. Note however, that the overall efficiency is 100%, so that no event escapes detection (sum of elements in each column equals 1). The measurement process now is simulated by multiplying the matrixÂ over the true distribution x, resulting in the measured histogram b: With vector b measured and the response matrix (10) given, one can try to unfold the true distribution. Singular value decomposition of 2 × 2 matrices is very simple, involving just a single rotation from left and another from right. As the matrix (10) is symmetric, the orthogonal matrices U and V should coincide. SVD can be performed explicitly by hand, and one easily obtains: with So, the two singular values of the matrix (10) are s 1 = 1, s 2 = ǫ.

Solving a linear system using SVD
Suppose now that the apparatus described by the matrix (10) has been used to measure numbers of events in a two-bin histogram Let B be the corresponding covariance matrix, which is especially simple for purely statistical errors in independent entries b 1 and b 2 : In order to solve the system, let us use U, S and V from (13) to rotate both the unknown vector x and the r.h.s. of the system b, in order to form a diagonal system of equations where The unknown vector x can now be easily obtained by rotating z back: Expression (19) gives the exact solution of the system (11) for whatever small but finite ǫ.
Formally, SVD of the matrixÂ means a decomposition of the r.h.s. b into a series of orthogonal and normalized functions of the discrete variable i = 1, . . . , n b . The basis is given by the columns of the matrix U, and the components of the vector d form the coefficients of this decomposition. Similarly, the vector of unknowns x is also decomposed into a series of ortho-normalized functions of the discrete variable j = 1, . . . , n x , given by the columns of the matrix V , while the coefficients stored in the vector z are new unknowns. After performing these transformations, the initial system of equations (11) is reduced to the diagonal system (17) which can be easily solved: the matrix S in (13) is diagonal and can be inverted by just inverting the singular values.
The inverse matrixÂ −1 exists for any ǫ = 0: Note that the expressions (19) and (20) are exact, so that SVD and the subsequent analysis can be considered as just another method of solving well determined full rank linear systems, maybe a bit too complicated but quite capable. If all components of the rotated r.h.s. d are statistically significant and if neither of the singular values s i of the matrixÂ is too small, the system (11) can be solved without any problem using any other method like, say, Gaussian elimination. But if ǫ is small the problem becomes ill-determined, and when in addition the r.h.s. is affected by measurement errors, the exact solution usually does not make any sense. In this case clearly it is not the exact solution we are looking for, and conventional methods of solving linear systems do not work. Usually they cannot even detect the problem.
To illustrate this, let us assume that the measured event numbers b 1 and b 2 satisfy the following relation: This means that the difference b 1 − b 2 is not statistically significant, so that the first term in the exact solution (19) is in fact a random number. But if ǫ is small enough (in this case -smaller than 1/ √ b 1 + b 2 ), this first term in both x 1 and x 2 dominates over the well-behaved and statistically significant second term, leading to almost arbitrary and senseless results. This phenomenon can be easily understood: for very small ǫ the apparatus is almost "blind", and one can hardly expect to determine x 1 and x 2 separately, unless the errors in b are sufficiently small.

Locating the difficulty
In order to trace the problem back to its origin let us reconsider the rotated system (17). Under the assumption (21) about the statistical accuracy of the data, the covariance matrix D of the rotated r.h.s. d is approximately diagonal, so one has the following (almost) independent set of equations: The first equation is good and gives a sensible result for z 1 , but the second one is not too useful even if ǫ is large, as the r.h.s. is in fact a random number. However, for ǫ close to unity it is at least harmless; it becomes dangerous for small ǫ, when the random number in the r.h.s. gets strongly amplified after being divided by ǫ, and after the rotation back to x according to (19), gives a huge and senseless contribution to both components of x in the exact solution (19). This means that z 2 cannot and should not be determined from equation (24), because effectively the matrix has insufficient rank and the system is over-determined. However sensible equation (24) may look, it should be replaced by another equation which contains as much information as (24), but is much less harmful. The value of z 2 is completely arbitrary and can be determined only from some external condition. Note that each value of z 2 will lead to a new vector x. The easiest thing to do is to put z 2 = 0, which would lead to the "shortest" x: x 2 = x 2 1 + x 2 2 ≡ z 2 1 + z 2 2 , as orthogonal transformations do not change the euclidean norm of a vector. Or, alternatively, one could choose the solution minimizing the variation (x 1 − x 2 ) 2 . In our small example both these alternatives lead to the same "regularized" vector Corresponding regularized covariance matrix has the form These two equations define the best reasonable estimate for an unfolded vector x, given the condition on the statistical accuracy of the data (21) and a poor detector with ǫ ≪

Rescaling equations and normalizing unknowns
Let us now get back to the full-scale problem defined in Sect. 3 and look at the initial linear system (2) from another viewpoint. It represents the solution of the following least square problem: and is adequate if the equations are exact, or if the errors in the r.h.s. are identical. This is not generally the case, as measurement errors in the vector b vary from bin to bin, and hence, different equations have different significance. In fact, one should consider a weighted least squares problem, where the following expression is being minimized: where ∆b i is the error in b i . The general case of (29) reads where B is the covariance matrix of the measured vector b.

Normalization of the unknowns
It is well known, that the exact solution of a well-determined linear system remains unchanged, if either the equations, or the unknowns, or both are rescaled. However, in the cases under consideration (where n x ≤ n b and some singular values are small) the minimization of (30) leads to an overdetermined system which should be solved in the least-squares sense. In this case any rescaling of equations and/or unknowns changes the singular values of the system and hence the solution as well. One can suggest various ways of rescaling, and some of them may lead to a serious improvement in the system behavior. One of our most important tasks is to optimize the system by rescaling it so that significant information is not suppressed while non-significant is not enhanced. One can, for instance, divide an unknown x k everywhere in the system by a number λ, multiplying simultaneously all corresponding coefficients A ik , i = 1, . . . , n b by the same number λ. Choosing various λ's for different k's, one can obtain substantially different matrices. We argue that one particular choice of rescaling coefficients is the most suitable for our purposes, provided the probability matrixÂ is obtained using a Monte Carlo simulation procedure (see Sect. 3).
Consider a new unknown vector w j = x j /x ini j , which measures the deviation of x from the initial Monte Carlo input vector x ini . If one now multiplies each column of the probability matrixÂ ij by the corresponding number of events generated in this bin x ini j , the system becomes where A ij is no longer the probability, but rather the actual number of events, which were generated in bin j and ended up in bin i 2 . Obviously, x = x ini corresponds to all components of the vector w being equal to 1, so that b ini i = nx j=1 A ij . At the end of the unfolding procedure, in order to obtain the correctly normalized unfolded solution x j , one has to multiply the unfolded vector w by x ini : Of course, if the number of generated events is the same for each bin, x ini = const, then the probability matrixÂ and the number-of-events matrix A coincide up to an overall constant factor which is completely irrelevant for our analysis. The systems (2) and (31) are completely equivalent for any shape of x ini , if the exact solution is required, but there are two serious reasons why, for the class of problems considered here, (31) is much better suited.
The first reason is fairly simple: if the initial Monte Carlo distribution x ini is physically motivated and is reasonably close to the one being unfolded, the unknown vector w should be smooth and should have small bin-to-bin variation, thus requiring less terms in the decomposition into orthogonal functions. This in turn means that more accurate unfolding should be possible, as fewer unknowns are required in order to obtain the unfolded solution.
The second reason is more technical and is connected to the singular value analysis. Whatever high statistics is generated in order to obtain the matrix A, some of its columns and/or rows may contain very few events, and some elements A ij may have just a single entry. In the probability matrix, these elements will contain the largest possible value of 1, unjustifiably giving a high weight to that particular equation and unknown, and the fact that this element has a 100 % error is completely ignored. At the same time, highly populated columns with statistically well-determined elements usually contain values significantly smaller than 1, due to finite resolution and limited acceptance. This clearly makes the probability matrixÂ a bad choice. On the contrary, the elements of the number-of-event matrix A are large if the generated statistics is large, and vice versa, thus giving a larger weight to better determined equations and unknowns.
The latter argument is in fact based on a quite formal consideration using perturbation theorems for the singular values of a matrix [6], and is an attempt to account for the intrinsic errors in the response matrix. It can be shown that for the cases of interest, initial Monte Carlo statistics should be at least one or two orders of magnitude higher than the statistics of the measured data. If so, the error in the unfolded solution based on (31) is dominated by the measurement errors in the r.h.s. b. Note that this is not true for the system based on the probability matrix (2), where huge errors may be introduced because of the fact that scarcely populated areas of the response matrix have far larger weight than they deserve.

Rescaling equations
The very form of (29) clearly suggests the way of rescaling the equations: after dividing each equation by the corresponding error ∆b i one obtains a balanced system, where all the equations have equal weights.
If B is not diagonal, equation rescaling becomes slightly more complicated but still straightforward. Being a covariance matrix, B should be symmetric and positive-definite, so its SVD yields: Substituting B −1 into (30) one sees that after the rotation and rescaling of both the r.h.s. b and the matrix A, the expression being minimized looks very simple again, and the minimization leads to the following system: The covariance matrix of the rescaled r.h.s.,B, is now explicitly made equal to the unit matrix I, and all the equations have equal importance.

Regularization and unfolding
The transition from (28) to (35) changes the appearance of the system from (2) to (36).
The singular values of the matrix are also changed, but the main problem with small singular values still remains. The exact solution of (36) will again most certainly lead to a rapidly oscillating distribution, which may have a smaller amplitude but is still useless. This spurious oscillatory component should be suppressed, using some a priori knowledge about the solution. Technically this can be achieved by adding the regularization or stabilization term to the expression to be minimized (see [6,1,4] and references therein): Here C is a matrix which defines the a priori condition on the solution, while the value of the regularization parameter τ determines the relative weight of this condition. For example, the choice C ik = δ ik would minimize the euclidean norm of the vector w, and if τ is set to be infinitely large, this would result in a vector w j = 0 for anyÃ andb. While the optimal value of τ is very much problem-dependent and its determination is an important part of our procedure, the explicit form of the matrix C should be chosen from general considerations. The common belief is that the solution histogram w should be smooth, with small bin-to-bin variation. Let us define the "curvature" of the discrete distribution w j as the sum of the squares of its second derivatives: Then the choice will suppress solutions w having large curvatures. Minimization of (37) leads to a new linear system, which has n x additional equations: This system is clearly over-determined, and one can apply SVD to the (n b + n x ) × n x matrix in the l.h.s. in order to solve it. This is possible, but would require calling SVD for each value of τ . Fortunately, a more efficient method (called sometimes damped least squares [6] ) can be suggested, which allows to express the solution of (40) for any τ through the solution of the initial non-regularized problem corresponding to τ = 0. The first step is to make the regularization term proportional to the unit matrix I: For τ = 0 the system (41) is equivalent to (36), if the inverse C −1 exists and can be safely calculated. The "second derivative" matrix (39), however, is apparently degenerate (every column and every row sums up to zero) so some measures should be taken to make the exact inversion possible. The easiest thing to do is to add a small diagonal component, C ik ⇒ C ik + ξ · δ ik , with ξ large enough to make the inversion possible, but small enough not to change significantly the condition of minimum curvature. In most cases, ξ = 10 −3 or 10 −4 is a good choice. C now is a symmetric non-singular matrix, which can be inverted using standard techniques 3 .
Let us now solve the system (41) with τ = 0. First, one needs SVD to decompose the product of matricesÃ C −1 : Here, once again, U and V are orthogonal and S is diagonal, with non-increasing positive diagonal elements s i . The solution now proceeds as in the 2 × 2 example of section 4.2. Rotate bothb and Cw to obtain a diagonal system: The system now looks (and actually is) very simple: Note that because the covariance matrix of the r.h.s.b was made equal to the unit matrix, the orthogonality of U guarantees that the new rotated r.h.s. d also has a unit covariance matrix, i.e. the equations in (45) are completely independent and have identical unit errors in their r.h.s. Obviously, solving (45) one obtains the exact solution of the non-regularized system: and the true distribution x can be obtained by multiplying each w i by the corresponding x ini i . With τ = 0 there is no regularization, so this solution is as useless as it used to be. But the solution of the system (40) with nonzero τ can now be found very easily, 3 One has, however, to be careful, as the matrix is too close to a singular one, and some of the standard routines may not work for small ξ. E.g., RFIN, RSINV and RSFINV have failed for ξ = 10 −4 , while RINV was successful. So, it may be convenient to use SVD once again for this purpose: decompose C = U C S C V T C and then calculate using the procedure explained in detail in [6] (Chapter 25, Sect.4). In short, introducing non-zero τ is effectively equivalent to changing d i by a regularized distribution: so that the solution of the rotated system becomes One can now see how nonzero τ regularizes the singularities due to small s i 's, effectively working as a cutoff for a low-pass filter, if Fourier-transform terminology is used. Indeed, s i is small when the index i is large, which in general corresponds to quickly oscillating singular vectors (i.e. columns of U and V ) defining the new basis in the rotated space.
Continuing the analogy with Fourier analysis, one can mention that the cutoff provided by the above regularization procedure happens to be quite smooth, thus avoiding specific quasi-periodic fluctuations of the solution known as the Gibbs phenomenon.
The covariance matrices Z and W of the solutions (48) can now be easily calculated: Now in order to obtain the true unfolded distribution x and its covariance matrix X one has to multiply w and W by the initial Monte Carlo distribution x ini : It is important to note that while (51) and (52) are regularized and as such depend on the value of τ , the inverse of the covariance matrix X −1 (which should be used for any χ 2 calculation involving the unfolded distribution (51)), is regular and readily calculable: In fact, X (τ ) defined by (52) is the effective pseudoinverse of the matrix (53). This means that while the equation is valid as if X (τ ) were the true inverse of X −1 , for a different ordering one has (see [5,6]): It may be interesting to write out the exact inverse covariance matrix for the 2 × 2 example of Sect. 4: The regularized covariance matrix (26) is actually the pseudoinverse of (56). As expected, the latter is perfectly regular in the limit of zero ǫ. The standard statistical tests can be used to check whether the last n x −k components of d i are compatible with the expected normal distribution with zero mean and unit variance. If this is not the case, then the errors in the measured data (or maybe the response matrix itself) are not estimated correctly, and we recommend finding out the reason of the discrepancy before proceeding with the unfolding. If, for example,the actual measurement errors in b are under(over)estimated, then the variance of d i for i > k will be smaller (larger) than 1. Moreover, if some additional correlations exist in the measured data which are not accounted for in the covariance matrix B, then log|d i | may steadily decrease for all i, though probably for i > k the slope will be different.
All this shows that the analysis of the plot of log|d i | vs i is of great interest by itself, being able to reveal the actual level of understanding the measurement errors in the experiment described by the simulated matrix A.
Anyway, if the number of statistically significant equations is determined to be equal to k, the regularization parameter τ should be put equal to the square of the kth singular value s k ≡ S kk of the matrixÃC −1 , determined in (43): With τ given by (57), the unfolded vector x, its covariance matrix X and the inverse of the latter X −1 are completely defined by corresponding equations (51),(52) and (53) and can be easily calculated, forming the solution of the unfolding problem. Yet another (and maybe more convincing) way of determining τ is to generate a test distribution which is close to the expected true one, but still significantly different from the initial Monte Carlo distribution x ini . Then one should simulate the measurement process by applying the response matrix to it, and add corresponding random statistical errors to the thus obtained "measured" distribution. The described unfolding procedure should be applied to the latter, and the best choice for τ is the one giving the smallest χ 2 between the test and the unfolded distributions (see our second example in Sect. 9).

The algorithm
In this section we present the concise description of the complete unfolding algorithm. The algorithm is linear (i.e. contains no loops) and can be divided into three distinct parts: initialization, rescaling/rotation and actual unfolding. Each step includes references to relevant subsections and equations.

Examples
The use of the unfolding procedure described above is now illustrated with two examples.
The first one is rather academic, and we have included it only because the same example was used in [1] and [4]. The response functionÂ(y, y true ) is given bŷ Then the probability response matrix was built according to eq. (6) for 40 equidistant bins in the interval (0,2) for both y true and y. The matrix is presented in Fig. 1a. The true continuous distribution is taken to be After the convolution (5), the distribution B(y) was discretized according to (4). Simulating a counting experiment, a random normally distributed error was then added to each entry, assuming the overall initial statistics of 5000 events. The resulting "measured" distribution b is plotted in Fig. 1b by a dotted line. The distortions caused by the measurement process can be seen by comparing the latter to the true distribution (59), shown by the solid curve in Fig. 1b. The unfolding algorithm described above was then applied to the distribution b. Fig. 1c shows the plot of the rescaled and rotated r.h.s. vector d. The solid line corresponds to the actual measured histogram, and the horizontal dashed line shows the one standard deviation statistical error in d i , which is equal to 1 for each i. One can see that after i = 10 the components d i are clearly non-significant. The flatness of this distribution for i > 10 and its apparent compatibility with the expected gaussian distribution with zero mean and unit variance, is in fact a test of the gaussian random number generator used to generate the errors in the measured histogram b. If we limit ourselves to less than 10 equations, we lose some significant information. Namely, the choice k = 1 leaves effectively only one equation, and the obtained "unfolded" distribution x (1) will be nothing else but a constant. On the contrary, by taking more than 10 equations one includes rapidly oscillating components with non-significant (and large) coefficients determined by the ratio d i /s i . In this particular example taking k = 40 would result in a distribution x wildly oscillating with the amplitude of about 5000.
The shape of the distribution d suggests that the effective rank k should be put equal to 10. The dashed line on Fig. 1c shows the regularized distribution d (τ ) calculated using (47), with τ = s 2 10 . It is interesting to compare this distribution with the similar one calculated for the exact true distribution (59) by the same procedure of rescaling and rotation, but without adding the random error (the dotted histogram in Fig. 1c). One can see that the regularized distribution is quite close to the true exact one.
The obtained distribution d (τ ) is then used to calculate the unfolded histogram x (τ ) , plotted in Fig. 1b (data points). It should be compared to the true distribution (59), shown by the smooth solid curve. Note that the error bars in x (τ ) account only for the diagonal elements of the covariance matrix X, and thus underestimate the actual errors. The correlations between adjacent bins x (τ ) j are quite significant, so one should use the exact inverse of the covariance matrix X −1 for any kind of χ 2 calculation involving the unfolded vector, and the regularized covariance matrix X (τ ) for the further error propagation.   1d presents the difference between the exact distribution and the unfolded one, together with the bands showing one and two standard deviation statistical fluctuations in the true exact distribution. In this scale, one can still see some oscillations of the unfolded solution, but they are well balanced, distributed almost uniformly, and are confined inside the two standard deviation band of the true solution, thus indicating that the genuine error is about twice as large as the true statistical one would be, if the measurements were exact. The average χ 2 over the 40 bins is equal to 0.9, meaning that the unfolded distribution is quite satisfactory.
In the second example we unfold a simulated spectrum of the invariant mass of two pions, corresponding to the ρ(770)-meson mass region. An artificial two-dimensional histogram reflecting a possible detector behavior was generated as the detector response matrix A, and is shown in Fig. 2a. This time it is a number-of-event matrix resulting from some Monte Carlo simulation, as opposed to the probability matrix obtained by the integration of some analytical response function used in the previous example. The matrix is far from being diagonal, and the initial simulated distribution x ini , shown by a dotted line in Fig. 2b, is hardly a constant.
The "measured" distribution b was obtained in a way similar to the first example: a distribution x test was generated (solid line in Fig. 2b), which has a behavior distinctly different from x ini . The measurement process was then simulated by the matrix multiplication: and finally a random gaussian error was added to each entry b i , simulating statistical fluctuations.
Rescaling and rotation results in a distribution d i plotted in Fig. 2c. One sees that the effective rank of the system k is close to 9, so the parameter τ should be set to the square of the 9th singular value of the matrix AC −1 . The components d i with i > 9 are clearly compatible with zero and have variances close to 1, thus confirming that the errors in the measured data are estimated correctly.
As in the first example, the choice k = 1 would leave us effectively with only one equation, and the obtained "unfolded" distribution x (1) will be nothing else but the initial Monte Carlo distribution x ini , shown by the dotted line in Fig. 2b. As for the solution of the non-regularized system with τ = 0, it would include all non-significant components and would oscillate rapidly within the range ±(2 ÷ 3) · 10 4 . This solution depends on the machine accuracy and obviously does not make any sense.
The regularized distribution d (τ ) i is shown by a dashed line in Fig. 2c. It is to be compared with the exact distribution d test corresponding to the vector (60) after the same procedure of rescaling and rotation, but without the random error added. One can see that the regularized vector is much closer to the true exact one for large i > 9.
The resulting unfolded histogram x (τ ) i is shown by the data points in Fig. 2b. The difference of the unfolded and the exact test distributions is presented in Fig. 2d, together with one and two standard deviation bands describing the statistical errors in the test vector. Here, too, the error bars show just the diagonal elements of the error matrix, which in fact contains quite strong bin-to-bin correlations. The agreement is very good indeed, especially if one considers the four orders of magnitude variation range of the test distribution.

Conclusion
The data unfolding method developed in this paper can be used in a wide range of applications, but is especially well-suited for high energy physics, where the response matrix is usually estimated by a Monte Carlo simulation of the measurement process, using some physically motivated initial distribution of the quantity under consideration.
The extensive use of a very versatile and flexible tool -the Singular Value Decomposition of a matrix -allowed us to derive a concise loop-free algorithm for data unfolding. Our choice of the regularization term results in smooth unfolded distributions which have the smallest possible curvature among the solutions satisfying the initial linear system in the least squares sense. This is achieved by suppressing the coefficients of high-order rapidly oscillating components of the solution. The suppression factor depends on both the statistical significance of the equation and the magnitude of the corresponding singular value. The regularized solution contains as much statistically significant information from the measured data as possible, simultaneously suppressing spurious, wildly oscillating components. This suppression is a natural result of the solution process, when the initial linear system is expanded to incorporate the minimum curvature condition.
The method used to solve the regularized system is extremely simple and reliable. An easy and straightforward way to determine the optimal value of the regularization parameter is suggested, which allows at the same time to test whether the quoted measurement errors are adequate. Note that the presented solution method is quite flexible and can be used with other choices of the regularization term as well.
Obviously, as the number of statistically independent data points is usually smaller (and sometimes much smaller) than the number of bins in the unfolded histogram, the latter will probably have significant bin-to-bin correlations. In our approach, full propagation of errors from the measured distribution to the unfolded one is implemented, and both the covariance matrix of the unfolded solution and its inverse are easily calculated. This allows one to perform further error propagation and parameter fitting without any problem, so, contrary to the viewpoint expressed in [1], we do not think that one should use fewer bins and custom bin boundaries for the unfolded histogram, in order to make the covariance matrix diagonal.
Obviously, curvature minimization introduces some systematic bias into the unfolded distribution, so the method will lead to acceptable results only if the true solution is indeed smooth, if the probability response matrix is used. However, when one uses the numberof-events response matrix, the condition of minimum curvature means that the deviation of the expected distribution from the initial Monte Carlo one should be smooth enough. This clearly allows one to use our procedure in cases when the measured distribution has some structure and/or a wide variation range, provided the initial Monte-Carlo has a similar behavior. If this is the case, then even for the small effective rank of the system, when the unfolded distribution happens to be quite close to the initial Monte Carlo, the former (in conjunction with the error matrix and its inverse) is still expected to give a usable solution of the problem.
Though we have tried to present comprehensive explanations of each step of the method, one should be able to use the algorithm formulated in Sect. 8 without having to understand in full all the details of the underlying mathematics. Note however, that the mathematics involved in this paper is still much simpler than that required by other regularization methods (e.g. [1,4]) because the most difficult tasks are successfully dealt with by the SVD procedure.
In short, the method presented in this paper allows one to unfold data obtained by any measurement process, if the response matrix of the detector is known. The use of the number-of-events matrix (and not the probability matrix) allows one to minimize additional uncertainties due to statistical fluctuations in the matrix itself, when the latter is obtained using the Monte Carlo simulation process. The concise linear algorithm results in a straightforward and simple implementation with full error propagation, including the complete covariance matrix and its inverse. The method is suitable for use in a wide range of problems, and can be generalized to incorporate multi-dimensional distributions.