{'fraction_non_alphanumeric': 0.026144837366759686, 'fraction_numerical': 0.008783567409305092, 'mean_word_length': 4.77911434236616, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Automatic Methods for Handling Nearly Singular Covariance Structures Using the Cholesky Decomposition of an Indefinite Matrix

Linear models have found widespread use in statistical investigations. For every linear model there exists a matrix representation for which the ReML (Restricted Maximum Likelihood) can be constructed from the elements of the corresponding matrix. This method works in the standard manner when the covariance structure is non-singular. It can also be used in the case where the covariance structure is singular, because the method identifies particular non-stochastic linear combinations of the observations which must be constrained to zero. In order to use this method, the Cholesky decomposition has to be generalized to symmetric and indefinite matrices using complex arithmetic methods. This method is applied to the problem of determining the spatial size (vertex) for the Higgs Boson decay in the Higgs ->4 lepton channel. A comparison based on the Chi-Square variable from the vertex fit for Higgs signal and t-tbar background is presented and shows that the background can be greatly suppressed using the Chi-Square variable. One of the major advantages of this novel method over the currently adopted technique of b-tagging is that it is not affected by multiple interactions (pile up).


Introduction
This paper describes an original automatic method for constructing constraints necessary to establish a consistent maximum likelihood problem in the presence of a nearly singular covariance structure based on the Cholesky decomposition of an indefinite matrix. The method is applied to the problem of constructing a vertex for extrapolated tracks in the context of a high energy physics data analysis. Track extrapolations involve a nearly singular covariance structure due to observational bias which affects high transverse momentum tracks. This observational bias is inherent in all track fitting algorithms used to reconstruct charged particle trajectories in solenoidal magnetic fields. Similar effects occur in other types of tracking problems, e.g., partially observed elliptical orbits in astrodynamics, etc., and are reminiscent of spherical aberration in ray tracking. Therefore an automatic method to handle the resulting nearly singular covariance structure warrants investigation.
We make use of the generalized Cholesky decomposition of a particular indefinite work matrix, K (k by k), in order to construct a restricted maximum likelihood function (ReML) via the method described in Smith (2001a). This method is generalized to handle complex arithmetic in the indefinite case. We selected the ReML technique because our problem is well matched to minimizing residuals with respect to the fixed effects without also estimating those fixed effects. In our particular case the fixed effects correspond to the mean position of a vertex. ReML eliminates the fixed effects from the likelihood. As will be shown the leading k − 1 diagonal elements can be matched with the expansion of the determinant which normalizes the ReML gaussian distribution and the k -th diagonal element can be matched to the χ 2 argument of the exponential term.
In addition, because the problem at hand involves a nearly singular covariance structure, the corresponding near-zero diagonal terms, as compared with a tunable cut parameter, of the same Cholesky decomposition designate those rows of K which fold into vectors which form auxiliary conditions such that the corresponding rows and columns of K drop out of the primary likelihood construction and are set aside to form a set of restrictions (constraints) in the maximum likelihood problem.
The procedure we use is based on the Newton-Raphson method and requires fast computation of the first and second derivatives of the likelihood function to construct the linear system for Newton-Raphson and to proceed through enough iterations before convergence is achieved. The calculation of these derivatives requires computation of the first and second derivatives of the Cholesky decomposition as well as the auxiliary conditions (as well as Lagrange multipliers depending on the details of the implementation). In order to speed up the numerical calculations we incorporated the method of backwards differentiation. The advantage of backward differentiation is that the gradient vector (m by 1) can be calculated in one pass through function evaluation, and also the Hessian (m by m) in m passes through function evaluation. Also, backward differentiation of the Cholesky algorithm comes with a greatly reduced computer memory requirement. Smith (1995) describes how to efficiently compute both forward and backward derivatives of the Cholesky decomposition by using methods taken from automatic differentiation (Griewank, 2008).
This technique permits variance-covariance estimation by restricted maximum likelihood (ReML).
While the Cholesky decomposition and its derivatives are finding applications with ReML and in linear models typical to studies of animal breeding (Meyer and Smith, 1995), these are coming with additional innovations (Meyer, 2001;Meyer and Kirkpatrick, 2005). The computational methods have been introduced in general statistical packages, including SAS (SAS Institute, 2009;Meyer, 2007).
Outside of ReML, the Cholesky decomposition and its derivatives are finding the following applications: in spatial modeling or kriging (Toal, et al., 2009); in lattice models with an example showing non-parametric curve fitting and cross-validation (Smith, 1997); in optimization involving the inverse diffusion problem (Christianson, 1997); to differentiate a Laplace approximation of a likelihood function, thereby permitting estimations of parameters in a population model (Frimannslund, 2006); to permit Hamiltonian Markov Chain Monte Carlo in the context of non-linear regression by function factorization (Schmidt, 2009); in calculating price sensitivities associated with exposure risk of financial portfolios (Capriotti and Giles, 2010); and for optimizing several hyper-parameters within a gradient-based machine learning algorithm (Bengio, 2000). (DeHoog, et al., 2011) have introduced a general notation for treating that task involved with calculating derivatives of functions that depend on triangular matrix factors (including Cholesky's factor), and they describe several application areas too.
The prior applications consider only the Cholesky decomposition of a positive-definite (or semidefinite) matrix. This convention is not followed in (Smith, 2001a and2001b) where consideration is given to a linear state-space model, nor is it followed in the present paper. Smith used the Cholesky decomposition (and it derivatives) to estimate second moments, but in this different application the Cholesky decomposition was applied to a symmetric and indefinite matrix. In general, the Cholesky decomposition may not be possible for an indefinite matrix, in which case an attempted decomposition may lead to an interruption. Such interruption is possible when a linear model includes effects that come with a singular variance-covariance matrix structure as is possible with state-space models. It has been the case that an interrupted Cholesky decomposition permits the correct ReML likelihood calculation following (Smith, 2001a and2001b). However, Smith's approach implies that the linear model is consistent with the singular variance structure, and in general this assumption cannot be made. The purpose of the present paper is to get beyond this limitation by introducing linear constraints that guarantee consistency, and hence this enlarges the set of cases where Cholesky-based ReML can be applied. There are also some minor typographical errors in the pseudocode listed in (Smith, 2001a), and the corrected pseudocode is presented in the Appendix (Section 12) of this paper. Furthermore, it is the additional purpose of this paper to get beyond the common ReML applications and demonstrate a new Cholesky-based discrimination method that is suitable for enriching data samples that are used in the search for the Higgs Boson ("Higgs Hunting") at the Large Hadron Collider (LHC).
A singular covariance matrix V can arise in very predictable situations in linear models. For example, in multinomial or mixture problems the dependence in the covariance matrix is well anticipated. Moreover, the coefficient matrix of the normal equation may be singular because of over-parameterization, even if V is non-singular. However, the typical solution of striking out rows and columns of the matrix V, or of a coefficient matrix to the normal equations, does not work well in the present situation where reliance on the Cholesky decomposition of an indefinite matrix is made. Removing rows and columns of the indefinite matrix runs the risk of removing too much information while ignoring computed constraints that come out naturally from the Cholesky decomposition. The fact is that a singular error matrix (that is also consistent with model specifications) implies that a linear combination of the residuals either sum to zero or to a linear combination of the fixed effects, and some model specifications (or parameters in the model) may be inconsistent with computed constraints. The constraints must be checked to avoid taking model consistency for granted, and this involves a completely new and different treatment of singularities compared to past treatments. The mixed linear model and ReML are reviewed in Section 2. In Section 3, likelihood evaluation is cast in terms of residual error, and the Cholesky decomposition of a symmetric and indefinite matrix. Section 4 provides mathematical justifications for the constraints that may be needed when the Cholesky decomposition is interrupted. To treat the possibility of constraints, Section 5 describes optimization by Lagrange multipliers and Section 6 describes optimization by penalized maximum likelihood. Both methods are readily adapted to the Cholesky decomposition and its derivatives. Section 7 treats estimation of fixed effects as an auxiliary calculation to the Cholesky decomposition by solving an indefinite system of equations. Section 8 presents an example coming from experimental physics, where the ReML method is used to devise a discrimination function to reduce backgrounds and preserve signal events in order to enrich data samples used for Higgs Hunting. Numerical stability and tuning of the algorithm is discussed in Section 9 and the conclusion follows in Section 10.

Mixed Linear Models, Log-Likelihood and ReML
Notation: In the following Sections, we denote the transpose of a matrix R by R ′ . We also denote Linear models, including mixed linear models, have found widespread use in statistical investigations. The linear model, though additive, is frequently flexible enough for real situations as an approximation around the mean. Also, linear models and the associated normality assumptions are well understood. Methods as old as the analysis of variance (ANOVA) are completely consistent with mixed model methods. Furthermore, while the theory is developing in new areas, such the Gibbs Sampler or with other Bayesian methods, mixed model methods benefit as new tools come along. The mixed linear model is represented by: where y is a vector of observations, β is a vector of fixed effects, u is vector of random effects and ǫ is the observational error. The matrices X and Z are incidence matrices that relate the various effects to observations. The first moments for the random effects (their expectations) are E[u] = 0 and E[ǫ] = 0, and the variance-covariance structure is given by var[u] = G, var[ǫ] = R and cov[u, ǫ] = 0. Additional assumptions are needed to implement maximum likelihood or computer simulation, and generally y, u, and ǫ are taken as multivariate normal. As indicated in (Goldberger, 1962), the Best Linear Unbiased Prediction (BLUP) of u is found by evaluatinĝ and whereβ is the Best Linear Unbiased Estimate (BLUE) of the fixed effects obtained by the Generalized Least Squares (GLS) problem These equations can be reformulated so that the solutions can be obtained directly from the mixed model equations (Henderson, et al., 1959) However, Section 7 (see below) provides and alternative method based on the method of (Siegel, 1965).
It is well known that if a non-informative prior is used to describe the fixed effects in a The log-likelihood for the Multivariate Normal (MN) is given by The maximum likelihood estimates of β and the dispersion parameters (R and G) are found by maximizing the log-likelihood. Estimates of the dispersion parameters can be badly biased by small-sample errors induced by the estimation ofβ. This is a serious problem when the dimension of β is large relative to the information available to estimate β.
To overcome this problem (Patterson and Thompson, 1971) introduced Restricted Maximum Likelihood (ReML), where the dispersion parameters are found by maximizing whereβ is the solution obtained by GLS, Eq.
(2). ReML has the advantage of estimating away the β parameters from the likelihood. This is especially useful in cases where one wants to concentrate on minimizing the deviations from a common mean, without explicitly finding that common mean. (Wikipedia, 2011) states, "In contrast to conventional maximum likelihood estimation, ReML can produce unbiased estimates of variance and covariance parameters." (Harville, 1974) derived the likelihood in Eq. (4) to treat the "error contrasts" which are found by taking a complete set of linear combinations of the observations which are sufficient to remove the effect of β while leaving the maximal amount of information for the purpose of ReML. An early review of ReML can be found in (Harville, 1977). More reviews can be found in (Speed, 1977 and. The relevance for the particular application described in Section 8 (see below), is that the χ 2 that is determined by ReML is independent of the central vertex coordinates. It is possible to back-out the coordinates of the fitted vertex, but it is not necessary to know the coordinates in order to determine the goodness of fit.

Applying ReML to the Z = 0 Scenario
Consider the following linear model L in the case where Z = 0 where y is an independent observation vector, β is a vector of fixed effects which are to be removed using the ReML likelihood, and ǫ is a vector of random residuals. The array X is the incidence matrix that assigns fixed effects to observations. The variance-covariance matrix of the random The linear model L is now fully specified and the error matrix satisfies V = R so that the likelihood given by Eq. (4) reduces to (Smith, 2001b) associates a symmetric matrix K with the above linear model L as follows: where the 0's denote the appropriate-sized null square matrices required to fill out the rows and columns of K and X ′ is the transpose of X. Therefore, L implies the existence of a matrix K, and K implies there exists a model L. Our main interest in Eq. (6) is that R is not required to be invertible. The method we use is able to identify those linear combination of y which are associated with the zero-eigenvalues of R. These linear combinations are non-stochastic and can be eliminated from the stochastic part of the likelihood and treated by the method of Lagrange constraints. In other words, our method is able to find the natural constraints required for the maximum likelihood problem for ReML.
It is well known that the Cholesky decomposition runs to completion with any matrix that is symmetric and non-negative definite. (Smith, 2001b) shows that the Cholesky decomposition can also be performed on the matrix K, and this is most curious because while K is symmetric, it is not non-negative definite. K is classified as indefinite. However, the rows and columns of K must be first permuted, leaving the last row and column in place as required. This is enough to permit computation of the Cholesky decomposition, or the lower triangular matrix L, where K k×k = LL ′ (for some permutation involving the first k − 1 rows and columns). (Smith, 2000 and 2001a) describes the ReML function, or the likelihood that removes the impacts of β, to be a function of this particular L. This calculation is given neatly as follows: where the absolute modulus function |L i,i | transforms a possible imaginary number into a positive number, and the summation is only over the non-zero pivots. When zero pivots are encountered (when L i,i = 0 for some i) we require that the i-th column of L vanishes and this is enough to justify likelihood evaluation by the above formula. It is easy to show that this calculation agrees with the likelihood given in Eq. (5) when R is non-singular. However, in the most general case we cannot simply skip the zero pivots, and more will be said about this in the Section 4, because particular constraints must be imposed.
The likelihood function is derived from elements of the Cholesky decomposition, and so there is nothing else that is needed to perform ReML but to find the derivatives that permit the likelihood to be optimized by the iterative Newton-Raphson technique. These derivatives come automatically with the Cholesky decomposition (Smith, 2000 and2001a), and so there is little beyond the matrix K, and its decomposition, that must be considered to describe ReML. The corrected pseudocode for the differentiation of the Cholesky algorithm and directions for how to use it are found in the Appendix (see Section 12).

Justification and Additional Constraints
The Cholesky decomposition of a matrix K k×k proceeds by identifying a non-zero diagonal element (the first pivot), and then permuting rows and columns of K to reposition that diagonal to the first diagonal position. Elementary row operations are now conducted to annihilate elements below the diagonal in the first column (this is pivoting). The first row of K is now transformed into the first column of L by replacing the first diagonal by its square root and by dividing the remaining elements in the first row by that square root. The first row and column of K are now deleted to produce a smaller sub-matrix of order k − 1. This sub-matrix is called the Schur complement. An outline computation of the first Schur complement (of When K is symmetric, the Schur complement is symmetric. Therefore, the Cholesky decomposition may proceed by half-storing K. The above steps are now repeated to generate a second Schur complement, then a third and so on. The algorithm may be organized to overwrite the half-stored K with L. The operations of the Cholesky decomposition are reversible, and therefore, the Cholesky decomposition conserves information. (Smith, 2001b) generalized the Cholesky decomposition for the case when K is symmetric and indefinite by using a complex representation for the Cholesky diagonals if needed. The sub-matrix initially containing the zero entries is the non-positive definite partition, and with fill-in (generated from pivoting from the diagonals of R) the non-positive definite partition becomes more negative.
When pivoting switches over to the diagonals of the non-positive definite partition, then fill-in results in the partition initially set to R (or the non-negative definite partition) thereby making it more positive. If a few pivots are selected from the diagonals of the non-negative definite partition first, before switching over to the non-positive partition and continuing pivoting over the negative diagonals until the non-positive definite partition returns to a matrix containing zero in all its entries (expect the last diagonal), then (Smith, 2001b) referred to the particular pivot order as a standard data reduction.
After the first pivot step within the Cholesky decomposition, the lead row and column is removed from K, and this reduces the dimension of the resulting Schur complement by 1 as noted above. After a standard data reduction, the Schur complement is again reduced in dimension from K but retains the special matrix form (ignoring the last diagonal with no loss in generality): where the subscript indicates that K 1 was generated from K following the Cholesky decomposition.
We will denote K as K 0 , to maintain consistency. In shorthand, this transformation is denoted by K 0 → K 1 . (Smith, 2001b) notes that K 1 also represents a model L 1 , given by where β 1 is sub-vector of β and the variance-covariance matrix of ǫ 1 is R 1 . The model L 1 is fully determined given the Schur complement K 1 . Moreover, those observations that have been processed and eliminated from K are uncorrelated with ǫ 1 . Also, the extracted information was accounted for in the previously constructed Cholesky diagonal elements. In other words, the standard data reduction separates the data into two statistically independent parts. The first statistically independent part is used to compute that block of the log-likelihood already known for ReML and consistent with Eq. (5). The remaining part pertains to the treatment of L 1 , but because K 1 is in the form of K, the process can be iterated and the Cholesky decomposition continued to the next pivot.
The Cholesky decomposition describes a sequence of standard data reductions: And these Schur complements correspond to a sequence of models: At each transition the dimensions of the Schur complement get smaller and smaller, and all along the way information is being processed correctly to evaluate the ReML likelihood. The only question is whether the Cholesky decomposition completes and leads to a final Schur complement that folds into the likelihood function calculation. However, the last Schur complement may be one of two special forms which cannot be reduced further. In this case, the last Schur complement corresponds to the last model L which is necessarily non-stochastic (i.e., the residual vector in L vanishes because its corresponding covariance matrix is computed to be null), but consistent with the model specifications that are implied by a singular variance structure. This last model is contained in a Schur complement of the First Special Form: In the case of an incomplete Cholesky decomposition, what is made manifest are the constraints that must also be imposed on the maximization of the likelihood function. The remarkable conclusion is that all the information is contained within the Cholesky decomposition (as we will see), even when the Cholesky decomposition is unable to finish.
One might question the appropriateness of the standard data reduction, given that the pivot order in the Cholesky decomposition may be dynamic and may not follow the standard data reduction. However, (Smith, 2001b) proved that pivot orders come in equivalence classes. Any dynamic order corresponds to a sequence of standard data reductions where the ReML likelihood is treated correctly. And moreover, Schur complements are invariant to the pivot order that generates them (including the last diagonal of the Cholesky decomposition), as well as the determinant calculated as the product of pivots. The correct likelihood is calculated even when the standard data reduction is not followed, because implicit in any pivot order is a sequence of standard data reductions.
With dynamical pivoting, the Cholesky decomposition is permitted to go as far as it can while skipping zero diagonals that would otherwise be pivots. Some zero pivots may encounter fill-in during the computation, and this permits the Cholesky decomposition to continue with additional When the Cholesky decomposition is unable to finish, then one of these systems of linear equations becomes a side condition (a constraint) for the likelihood maximization exercise.
We may be uninterested in the extraneous parameters. However, we are interested in consistent models that don't contradict themselves when linear combinations are made of their components.
This concern is quite independent of the extraneous parameters as we will see. The set of nonstochastic equations that are revealed (by the First Special Form of the Schur complement) will involve extraneous parameters, and these will be ignored in as much as ReML removes their impacts from the likelihood. The revealed set of equations can imply whatever they want about β I and ReML will ignore them. This assumes that the statistical model is already consistent.
However, the non-extraneous parameters that are identified for likelihood evaluation need not conform to a consistent linear model and that's where we depart from (Smith, 2001b

Constraints: Lagrange Multiplier Method
In the event that the Cholesky decomposition encounters a zero and is unable to complete, in order to to maintain consistency one can consider appending to the ReML likelihood a Lagrange multiplier expression according to one for the following cases: where maximization proceeds with respect to the non-extraneous parameters s 1 , s 2 , s 3 , s 4 (see Section 8), β I and λ. Therefore, the side condition enforces the restriction that v is in the column space of H independent of β I .
where maximization proceeds with respect to the non-extraneous parameters s 1 , s 2 , s 3 , s 4 (see Section 8) and λ. Therefore, the side conditions enforces the restriction that v = 0.
In both cases, we find an objective function so constructed from the known elements of the unfinished Cholesky decomposition. First and second derivatives are available (Smith, 2001a) and constrained optimization is straightforward by the technique of Newton-Raphson iteration as applied in Section 8 below.
In paradoxical cases, the data may produce an inconsistent model. In which case, we recommend constrained maximization as noted above, and in this way useful information is included that would otherwise be ignored.
In case the Cholesky decomposition is unable to finish and the Schur complement is of the Second Special Form, then the particular linear combination of elements corresponding to v can be constrained to zero by adding a Lagrange multiplier term λ ′ v as indicated above. An alternative method involves adding a quadratic penalty term to the likelihood (see Section 6) can be used to adjust the constraints to zero within estimated precision.
Case 1 and Case 2 above represent possible objective functions. However, we concentrate on Case 2 because that is the one relevant to our example. The objective function to be maximized is of the form F 2 (L, λ). The three algorithms of Section 12 are used to construct the Newton-Raphson linear system. If required, the first derivatives with respect to the Lagrange multipliers (the λ vector) come directly out of the v elements of the corresponding Schur complement of the unfinished matrix L. Table 1 sketches the Cholesky decomposition which is used to calculate the components of the ReML likelihood and the possible non-stochastic linear combinations that will be subject to the constraint conditions. We initialize the array F of Table 2 to the array ∂F 2 (L, λ)/∂L ij that represents the constrained objective function including the Lagrange multipliers. We then find the mixed second derivatives directly from Q (see Table 3, which represents the forward derivative calculations.)

Constraints: Penalized Likelihood Method
There is another method of implementing constraints besides Lagrange multipliers. Instead of introducing a set of Lagrange multipliers which are treated as independent parameters to be varied in a maximum likelihood procedure, the terms which they multiply can be squared and then multiplied by a fixed constant. This is known as the "Penalized Likelihood" method. In the Penalized Likelihood method the Lagrange equations for Case 1 and Case 2 above are modified to the form where the c i are nominated constants and maximization proceeds with with respect to the nonextraneous parameters (s 1 , s 2 , s 3 , s 4 ) and β I (see Section 8).
Case 2 : where c i are nominated constants and maximization proceeds with respect to the non-extraneous parameters s 1 , s 2 , s 3 , s 4 (see Section 8). The squares of the column vectors in the last term on the RHS of the above are to be understood as vector dot products.
The main difference between the Lagrange Multiplier method and the Penalized Likelihood method is that the auxiliary constraints are imposed to the maximum machine precision by the Lagrange Multiplier method, whereas the constraints are imposed in a less extreme manner with the Penalized Likelihood method. The point of the nominated constants is to adjust the enforcement of the constraints to within reasonable limits set by the known precision of the measurements.
This allows the possibility, for example, of not demanding that the constraints be satisfied with more rigor than the measurement uncertainties can justify.
7 Estimating the Fixed Effectsβ (Siegel, 1965) described a method to compute generalized least squares (GLS) estimates by way of a system of equations with a symmetric and indefinite coefficient matrix. In our notation, this system of equations is presented below: whereβ is the GLS estimate of β, and λ are the Lagrange multipliers. This Section will follow (Smith, 2001a), and show how to calculate bothβ and λ as adjunct operations that depict backward substitution, given that the Cholesky decomposition of K is available.
Note that the both the coefficient matrix and right-hand side of Siegel's equations are submatrices of K: the coefficient matrix is the lead sub-matrix in K, and the right-hand side is fixed to the last column (or row) and remains there for all permitted row-column permutations of K.

Now rewrite Siegel's equations in simple terms:
where C and r signify the coefficient matrix and right-hand side, respectively; and b is a column vector containing λ andβ. To solve these equations, we might permute the rows and columns of C, and compute the Cholesky decomposition: TT ′ = C. To permute the rows and columns of C, and the rows of b and r are also permuted to leave Siegel's equations intact. Now multiplying both sides of Siegel's equations by the same elementary row operations that transform C into T ′ gives: When C is non-singular, a = T −1 r. The coefficient matrix in this new system of equations is upper triangular. Therefore, vector b can be solved by backward substitution.
With L computed, where K = LL ′ , note that T is the leading sub-matrix in L and a ′ is the last row vector of L (excluding the last diagonal). Therefore, having computed L we need only enter backward substitution to evaluate b. The GLS estimatesβ will be found scattered in b, noting that b is permuted.
When the i-th pivot encountered in L is zero, and the i-th column vanishes, a singularity is present and b is not estimated uniquely. We place a restriction on b: set its i-th entry to zero.
This modification is implemented when the Cholesky decomposition is unable to finish but ends with a Schur complement of the Second Special Form. The associated column of L will also be constrained to vanish. When the Cholesky decomposition ends with the First Special Form then an auxiliary estimate of β I is available from the main optimization and this estimate is used in backward substitution given that what had been computed for the rest of L is lower triangular.
The parameters on the RHS of Eq. (7)  where the sign of the distance for D is given as positive if the track circle (projection of helix into the x-y plane) contains the z-axis and negative otherwise. A typical track fitting program will produce the 5-parameters (ρ, φ, θ, dxy, dsz) along with a 5-by-5 covariance matrix which can be used to estimate the spatial error matrix in terms of the 3-by-3 correlation matrix for x(s), y(s), z(s), for example, by the method of "propagation of errors". Once these track parameters are measured, then the track can be extrapolated to any position along its trajectory by using the parametric equation of a helix as a function of s. One of the interesting properties of these spatial extrapolations is the the 3-by-3 spatial error matrix is very close to rank-2 (almost singular). Fig. (1) illustrates the reason that the spatial error matrix has a near-zero eigenvector lying almost entirely in the x-y plane. If the radius of the track circle is much larger than the radius of the tracking detector, then only a small fraction of the circumference of the track is measured. This creates a measurement or observational bias in the sample of hits along the track. Such an observational bias creates a highly eccentric error matrix for the reconstructed x and y coordinates of the centerof-curvature. This highly eccentric error matrix produces a spatial error matrix along the track helix which has a very small eigenvalue in the direction of the projection of the track tangent vector in the x-y plane. This property is independent of the method used to determine the track parameters precisely because the tracks we are most interested in have high p T and, therefore, have very large diameter track circles compared to the diameter of the tracking chamber.
In order to determine if a Higgs candidate has been selected in the data-analysis of the experiment, one of the criteria applied to the 4 selected tracks is: Are these 4 tracks consistent with originating at a common position in space? Do these 4-tracks form a "vertex" in space? If the hypothetical Higgs Boson decayed into 4 charged leptons, then each lepton would follow its own helical trajectory, but all 4 leptons would converge on a common point in space corresponding to the point of decay of the Higgs Boson. Reducible 4-lepton background events, would be inconsistent with originating from the same common location. Therefore constructing a likelihood for the hypothesis that the 4 tracks originate at the same location is a very useful algorithm to separate signal events (Higgs candidates) from reducible backgrounds.
The tracks are assumed to come from a common point by hypothesis and the likelihood is obtained by first constructing the K matrix of Eq. (6) above with the y i ( for coordinates x(s i ), y(s i ), and z(s i ) for each track (i is numbered 1 through 4) with its respective independent variable s 1 , s 2 , s 3 or s 4 . The R matrix in Eq. (6) is obtained by the "propagation of errors" method which involves Taylor expansion of the functions x(s), y(s) and z(s) about the mean positions. The overall matrix is obtained by appending each of the contributions from each track to the form the assembly y = [y 1 , y 2 , y 3 , y 4 ] and similarly for R. The incidence matrix X is similarly constructed using X i = 3-by-3 unit matrices (for i = 1, 2, 3, 4) and stacking four such matrices one-above-the-other. Given 3-by-3 unit incidence matrices for X i = I 3 (for the i-th track), then β = r c . The extraneous parameters β are the central position coordinates of the vertex which can be estimated using the method of Section 7. This relationship completes the definition of the Linear Model.
Since the extrapolated spatial error matrices have the nearly rank-2 property, it is necessary to utilize the methods described above to construct a consistent ReML for testing the hypothesis of the common spatial origin of the 4 selected leptons as well as automatically constructing the nonstochastic linear combinations which have to be constrained to zero (a la the Second Special Form mentioned above). The constraints are determined based on comparing the size of the respective diagonal element in the Cholesky decomposition to a tunable threshold value which is optimized for signal events. When such a zero is encountered, the location on the main diagonal is noted and stored for subsequent determination of the constraint itself. The ReML likelihood and collateral constraints are determined as outlined above and the system of first and second derivatives are calculated using the Tables 1, 2, and 3 of the appendix. A linear system of equations is then formed by Taylor expansion of the objective function and a step is taken by Newton-Raphson (NR) iteration to move towards the maximum likelihood position. The coordinate functions of Eq. (7), evaluated at s = 0, produce the coordinates of the helix at closest approach to the beamaxis (z-axis). Since a real Higgs would be produced very close to the colliding beam axis, then s 1 , s 2 , s 3 , s 4 should all be initialized to 0 at the beginning of the first NR step for fast convergence.
Alternatively, one could perform a binary search about s = 0 (for each track) to determine appropriate initial values for the NR method.
After convergence, the algorithm produces a set of four parameter values s 1 , s 2 , s 3 , s 4 which give the location along each track such that the assumptions of the model are satisfied. The χ 2 value is obtained as the square of the last element on the main diagonal of the resulting decomposed K matrix (χ 2 = −L 2 k,k ). This variable is only approximately χ 2 coming with degrees of freedom that are underestimated (by the number of positive pivots minus the number of negative pivots plus 1, minus 4 from estimating the position variables s 1 , s 2 , s 3 , s 4 ) and we use the empirical distribution for χ 2 in applications. ReML and χ 2 are invariant with respect to the coordinates of the central position (the vertex), because the impact of the central position was removed from the likelihood.
The crucial assumption of the the model is that there be a common origination point for the four tracks somewhere. If the χ 2 from the fit is small, the model is satisfied and the hypothesis that the 4-tracks originated from a common point is consistent with the data. This is what we expect for Higgs → 4 leptons. Reducible background events (such as tt, Zbb → 4 leptons) where the 4-leptons do not originate from the same point will not fit the hypothesis of the model and will have large χ 2 values.
In order to implement the above algorithm a sample of events was generated using the PYTHIA (Sjöstrand, et al., 2006) Event Generation program. These events were then simulated using the LDT Monte Carlo Program (Regler, et al., 2007) and (Valentan, et al., 2011) to determine the detector response and to obtain fits to the track parameters and their covariance matrix. The analysis of the fitted tracks and their error matrices was done using the Rave/Vertigo dataharvesting and vertexing environment (Waltenberger, 2011) with the above ReML algorithm inserted internally in the package.
Applying the ReML method to the case of 4 leptons from Higgs decay and also from a known major source of background (tt → 4 leptons), where the leptons can be either electrons or muons, we arrive at a comparison of signal and background shown in Fig. (2). Only events in which all four tracks had p T > 5 GeV/c were used. As can be seen Fig. (2), even though the Higgs signal events have a smaller cross section and a smaller number of expected events than tt events, they are much more peaked near zero χ 2 for the same acceptance conditions of the detector. This results in a much higher efficiency (CDF value) for detecting Higgs as compared to retaining background events as can be seen in Fig. (3).
This means that we can preferentially select Higgs candidates and reject tt background by cutting on the χ 2 value.
In the case of high luminosity collisions, there can be multiple events detected every time the detector is triggered. These multiple events are recorded by the detector in the same time window  algorithm and the 4-track vertex method will assign a large χ 2 to that event and it will be rejected.
However, if the wrong Primary Vertex is used as a reference for the b-tagging then all 4-tracks will register a large b-tagging discriminant value even though they might be matched to a common spatial position (i.e., have small χ 2 from the vertex fit) by the 4-track vertex algorithm. These events should not be discarded just because the wrong and irrelevant Primary Vertex was used as a reference marker. If the 4-tracks have a small value of χ 2 from the 4-track vertex method, then that alone is sufficient to decide if they should be further analyzed.
The 4-track vertex method involves a tunable parameter α, which determines if a pivot (diagonal element) is classified as a near-zero (constraint) or to be included in the likelihood construction.
Ideally α corresponds to the zero-precision of the machine. During the Cholesky algorithm processing steps pivots are compared with α. If the Cholesky pivot value is less than α, then the corresponding row associated with that pivot is counted as a constraint row and marked for inclusion in the Lagrange multiplier (or penalty) terms. With α tuned to maximize the number of reconstructed 4-track vertices on a Monte Carlo signal sample, the algorithm identified 2 constraints per four-track event. Other methods of vertexing based on "rules of thumb" concepts in track extrapolation apply a constraint for every track leading to 4 constraints for a four-track combination (e.g., so-called "working in the perpendicular plane" for each track). If α is set too low, then the numerical instabilities associated with attempted inversion of a nearly singular covariance structure appear. If α is set too high, then more constraints than necessary are imposed and the constraints begin to involve the fixed effects for too large a number of identified constraints corresponding to too large a value for α. In practice the distribution of pivot values for signal Monte Carlo events can be used to tune the size of alpha for good distinction between the near-zero and the non-zero pivots. Once tuned, the algorithm reconstructs all the vertices consistent with the hypothesis of the associated tracks originating from a point.
Also convergence performance of the Newton-Raphson iteration can be fortified by using binary search method to check the initialization of the algorithm. Unless tracks are constructed pathologically, or have very low values of p T and curl up in the detector (curlers), then all the vertices converge in less than 30 iteration steps and much fewer for the case of signal events.
Failures to converge are associated with a set of tracks which are inconsistent with the hypothesis that they originate from a single point in space. In the above analysis tracks were restricted to the set with p T > 5 GeV/c to avoid curlers.

Conclusion
This paper presents an automated matrix procedure for constructing the ReML likelihood function and also to identify non-stochastic linear combinations that represent quantities that must be constrained to zero in the event that the error matrix describing the data is singular. The method is applied to the problem of determining how close four tracks from a Higgs → 4 lepton decay approach a common position in space. This method can be used to discriminate between signal and background events at experiments at the LHC. It has an advantage of being independent of pile up which is not the case for b-tagging.

Acknowledgments
We are especially thankful for the assistance of Nicola 12 Appendix: Pseudocode for Differentiation of Cholesky Decomposition This Section includes corrections (shown in red) to the pseudocode which was given in (Smith, 2001a). Construction of the Cholesky decomposition for an indefinite matrix K N ×N is summarized in the following table Initializations: L is lower triangular and L ij ← ij-th element of K N ×N (x 1 , x 2 , ..., x p ).
F N ×N is a work space with elements F ij , defined respectively for i ≥ j.
S N ×N and Q N ×N are work spaces with elements S ij , Q ij , defined respectively for i ≥ j.