Paul Avery CBX 91 − 72 October 18, 1991

Applied Fitting Theory I General Least Squares Theory

Section 1: General Statements

Section 2: Simple Least Squares Fitting

Section 3: Calculation of Errors: The Covariance Matrix

Section 4: Fitting with Correlated Errors

Section 5: Fitting with Constraints

Section 6: Solving for Unknown Parameters in the Constraint Equations Section 7: Factorization of Constraints Section 8: Simple Application of Unconstrained Track Fitting Section 9: Simple Application of Constrained Track Fitting Appendix: Tricks for Inverting Large Matrices

1. General Statements

This paper is the first of four documents I am writing on general fitting theory, particularly as applied to high energy physics experiments. I have several reasons for doing this. First, I want to provide a general introduction to the subject using matrix notation for people who have not done much fitting, or who are not familiar with the very elegant matrix formulation of least squares fitting theory. Second, I want to demonstrate how constraints, even those that introduce new variables, can be incorporated into the fit (using Lagrange multipliers) and how they affect the error calculation. Third, I want to show some new techniques I have developed that greatly extend the usefulness of least squares theory. Finally, I wanted to provide a single reference for people who use least squares techniques frequently but who are unable to find practical algorithms such as the ones covered here in standard books and articles on probability and statistics.

This document replaces an earlier note, CBX 87-38, that I released in 1987. The older note served as an introduction to least squares theory, especially in the use of constraints, and was supposed to be the first of a three part series. Unfortunately, it didn’t work out that way since the followup notes were never released (although they had been written).

Since that time I have developed some new factorization algorithms for applying constraints and inverting large matrices that I felt are important enough to warrant inclusion. These methods, and the ones discussed in the original paper, have particular relevance to kinematic fitting, a favorite topic of mine. One important consequence is that it should be possible to speed up the process of kinematic fitting so that it becomes a more practical tool for high energy physics analysis. In fact, I am writing a kinematic fitting package that uses the algorithms discussed in this article. It is discussed in a separate writeup.

My plan is to provide in in the first paper some general formulas for fitting, with and without constraints. These formulas are used extensively in fitting data to histograms, determining track parameters from hit information and kinematic fitting. In the second note, I will discuss how to exploit least squares fitting to determine systematic detector effects, which are either too small to be determined in a single event by fitting, or too intertwined with other effects to be reliably estimated. In the third paper, I will present an analysis of how track parameter errors due to multiple scattering can be treated exactly, although the best possible parameters are not obtained. The technique, which I call Nonoptimal Least Squares fitting, can be applied to any problem in which correlations among large numbers of data points prevents a direct attack by traditional methods. Finally, in the fourth article I will discuss fitting techniques specific to kinematic fitting, including ways to factorize the fits and improve execution speed.

I have no intention of making this series of notes a complete reference on fitting theory. For the most part I assume Gaussian errors and do not worry about fitting low statistics histograms, which requires special techniques. Consequently, I omit discussion of how confidence intervals, upper limits and asymmetric errors are calculated.

A word on notation. Throughout this series lowercase bold letters (x, y) always refer to vector quantities and uppercase letters (A, B) represent matrix quantities. The symbol V always refers to a covariance matrix.

2. Simple Least Squares Fitting

Let’s start off with an easy problem in least squares fitting that demonstrates the basic principle and which also illustrates the compactness and power of matrix notation. Suppose we have n independently measured data values y_{l} which we believe are functions of m unknown variables α_{i}, where m ≤ n. In other words, we are hypothesizing the relation y_{l} = f_{l}(α_{i}). Note that the functions f_{l} depend on the data point in question, as can be seen easily by thinking of the index l as tracking the “independent variable” x at each point.

Since each y_{l} is a measured quantity with a standard deviation σ_{l}, the equation y_{l} = f_{l}(α_{i}) cannot be satisfied exactly if m < n. However, we can require that the equation be satisfied “as closely as possible” by defining the χ^{2} statistic, where

χ^{2} = ∑_{l} \frac{(y_{l} − f_{l}(α))^{2}}{σ_{l}^{2}} ,

and demand that the values of α_{i} be chosen so as to minimize χ^{2}. Note that χ^{2} is merely the sum of the squares of the “misses” weighted by the uncertainty σ_{l} of each data item (each of which is, for now, independent of all others) so small values are better than large ones.

Let me expand a little on what is going on here. We assume that each measurement y_{l} has a measurment error σ_{l}. If the errors are Gaussian and the expected value for the measurement is f_{l}(α) then the probability density function (pdf) for the measurements can be written

g(y_{l}) = \frac{1}{\sqrt{2π}σ_{l}} exp [−(y_{l} − f_{l})^{2}/2σ_{l}^{2}],

where the normalization is chosen to make ∫_{−∞}^{∞} g_{l}(y_{l})dy_{l} = 1. The joint pdf for n measurements is then

\begin{aligned}g(y_{1} … y_{n}) & = ∏_{l} \frac{1}{\sqrt{2π}σ_{l}} exp [− ∑_{l}(y_{l} − f_{l})^{2}/2σ_{l}^{2}] \\ & = ∏_{l} \frac{1}{\sqrt{2π}σ_{l}} exp (−χ^{2}/2) \end{aligned}

It is clear that maximizing this pdf is equivalent to minimizing χ^{2}. The adjustment of the parameters α to achieve the maximization of the pdf is called the “maximum likelihood” method, and can be carried out even when the pdf is non-Gaussian, e.g., Poisson statistics.

The requirement that the χ^{2} function be at a minimum is usually guaranteed by the m relations

\frac{∂χ^{2}}{∂α_{i}} = 0.

Unfortunately, no general method exists for solving these m equations since f_{l} and its first derivatives can in general be highly nonlinear functions of the parameters α_{i}. Often, however, the functions vary slowly enough so that one can expand to first order about an approximate solution, viz.

f_{l}(α) = f_{l}(α_{A}) + ∑_{i}(α_{i} − α_{iA})\frac{∂f_{l}(α)}{∂α_{i}}_{α_{A}} ≡ f_{lA} + ∑_{i} A_{li}η_{i}.

Under this approximation, the χ^{2} function can be written

χ^{2} = ∑_{l} \frac{(y_{l} − f_{lA} − ∑_{i} A_{li}η_{i})^{2}}{σ_{l}^{2}} ≡ ∑_{l} \frac{(Δy_{l} − ∑_{i} A_{li}η_{i})^{2}}{σ_{l}^{2}} .

Under very general conditions it can be proved that the parameters η_{i} = α_{i} − α_{iA} obtained by minimizing this linearized version of the χ^{2} function are unbiased and have minimum variance σ_{α_{i}}^{2}, i.e., one cannot find a better set of parameters.

The minimization conditions, ∂χ^{2}/∂α_{i} = 0, can be written in this linear approximation as

∑_{l} \frac{A_{li}A_{lj}}{σ_{l}^{2}} η_{j} = ∑_{l} \frac{A_{li}}{σ_{l}^{2}} Δy_{l},

where 1 ≤ i, j ≤ m. Defining the symmetric matrix on the LHS as

V_{A ij} = (∑_{l} A_{li}A_{lj}/σ_{l}^{2})^{−1} ,

we get the solution

η_{i} = ∑_{j,l} V_{A ij} \frac{A_{lj}Δy_{l}}{σ_{l}^{2}} .

Surely, you say, there must be a better way of writing these equations without worrying about all the indices. To see how this might be done, let’s go back to the linearized χ^{2} equation and write it in matrix form:

χ^{2} = (Δy − Aη)^{t}V_{y}^{−1}(Δy − Aη)

where A is the matrix of coefficients defined above and V_{y}^{−1}, η and y are given by the expressions

V_{y}^{−1} = ( \begin{matrix}1/σ_{1}^{2} & 0 & ⋯ & 0 \\ 0 & 1/σ_{2}^{2} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & 1/σ_{n}^{2} \end{matrix} ) , η = ( η_{1} ; η_{2} ; ⋮ ; η_{m} ) , y = ( y_{1} ; y_{2} ; ⋮ ; y_{n} ) .

V_{y}^{−1}, the inverse of the measurement covariance matrix, is called the “weight matrix”.

Taking the partial derivative with respect to the parameters η we obtain the m equations

A^{t}V_{y}^{−1}(Δy − Aη) = 0,

which can be written

A^{t}V_{y}^{−1}Aη = A^{t}V_{y}^{−1}Δy.

Defining as before V_{A} = (A^{t}V_{y}^{−1}A)^{−1}, we finally obtain the m equations

η = V_{A}A^{t}V_{y}^{−1}Δy,

which is a much more understandable and meaningful expression than the previous one. The parameters α can be determined from α = α_{A} + η. Note that V_{y}^{−1} has dimensions of n×n, A is of size n×m, and V_{A} is of size m×m. The advantage of the matrix notation is that we can easily put in the indices by inspection after the equation is written in compact form.

Sometimes it is useful to calculate the “residuals” r = Δy − Aη = y − f_{A}(α_{A}) − Aη, which are merely the differences of the measurements and predicted values. The χ^{2} can then be written χ^{2} = r^{t}V_{y}^{−1}r = ∑_{l} r_{l}^{2}/σ_{l}^{2}.

3. Calculation of Errors: The Covariance Matrix

Obtaining a solution for a set of parameters is not enough: we also want to know the errors of the parameters and their correlations with one another. This information is contained in the so-called covariance matrix, which is defined as

cov(x_{i}, x_{j}) ≡ V_{x ij} = ⟨(x_{i} − x̄_{i})(x_{j} − x̄_{j})⟩ ≡ ⟨δx_{i}δx_{j}⟩,

where the symbol ⟨⟩ refers to a weighted average taken over all possible values of the enclosed expression. Note that the i^{th} diagonal entry in the covariance matrix is merely σ_{x_{i}}^{2}, the square of the standard deviation of x_{i}. The covariance matrix can be written in matrix notation as V_{x} = ⟨δxδx^{t}⟩. For example, the covariance matrix for the (independent) data values y_{l} is

V_{y} = ⟨δyδy^{t}⟩ = diag(σ_{i}^{2}).

This explains how I chose the notation for the weight matrix V_{y}^{−1} in the previous section. In general the covariance matrix V is related to another matrix called the correlation matrix ρ. The elements ρ_{ij} are computed from

ρ_{ij} = \frac{V_{ij}}{\sqrt{V_{ii}V_{jj}}} = \frac{V_{ij}}{σ_{i}σ_{j}} .

The element ρ_{ij} (i ≠ j) is called the ij correlation coefficient and satisfies the inequality −1 ≤ ρ_{ij} ≤ 1.

Finding the errors in the parameters η is now straightforward. Starting with the solution η = V_{A}A^{t}V_{y}^{−1}Δy we obtain

V_{α} = V_{η} = ⟨δηδη^{t}⟩ = V_{A}A^{t}V_{y}^{−1}⟨δyδy^{t}⟩V_{y}^{−1}AV_{A} = V_{A}A^{t}V_{y}^{−1}V_{y}V_{y}^{−1}AV_{A} = V_{A},

and the errors are σ_{η_{i}} = \sqrt{V_{η ii}} = \sqrt{V_{A ii}}. I have used δΔy = δy and the definition V_{A} = (A^{t}V_{y}^{−1}A)^{−1}. The fact that V_{η} = V_{A} might seem surprising, but it turns out that the least squares technique, with or without constraints, always yields an equation of the form V^{−1}η = Cy and gives V_{η} = V.

Once a parameter error is known, it is easy to figure out for Gaussian errors the probability that the parameter will fall within ±1σ, ±2σ etc. The following table gives some typical values.

# standard deviations Probability to fall in ±nσ
0.5  0.3830
1.0  0.6826
1.5  0.8664
2.0  0.9554
2.5  0.9876
3.0  0.9974
3.5  0.9996

Of course, in real life measurement errors are rarely distributed like Gaussians and usually have some significant non-Gaussian tail. However, the Central Limit Theorem usually guarantees that the estimated parameters will be more Gaussian than the estimated measurements, provided that there are a sufficient number of measurements,

Given a set of parameters and their estimated errors, what assurance do we have that the fit is any good? To answer this question, we note that the entire least squares procedure is concerned with minimizing the χ^{2}, so the final value of χ^{2} provides a measure of how well the fit worked. Traditionally, one defines a “confidence level” which is the probability that a fit will have a χ^{2} greater than or equal to the value we found from the fit. Mathematically, we write

C.L. = ∫_{χ^{2}}^{∞}f(x)dx,

where f(χ^{2}) is the χ^{2} probability distribution function.

For Gaussian errors we can calculate f(χ^{2}) straightforwardly. The χ^{2} can be written in the compact form

χ^{2} = Δy^{t}V_{y}^{−1}Δy − η^{t}V_{η}^{−1}η,

which is the difference of two positive definite forms having n and m terms, respectively.

It is straightforward to show that if the y are distributed as Gaussians (correlated), they can be rotated and rescaled in such a way that the χ^{2} can be written in the manifestly positive form χ^{2} = w_{1}^{2} + w_{2}^{2} + ⋯ + w_{ν}^{2}, where the w_{i} are independent Gaussians having unit error and zero mean. ν = n − m is called the number of “degrees of freedom” of the fit. In that case, applying the Jacobian transformation to the joint pdf with \sqrt{χ^{2}} as the “radius” and ν − 1 angular variables over which we integrate, we find the χ^{2} probability distribution to be

f(χ^{2}) = \frac{(χ^{2})^{\frac{ν−2}{2}} e^{− \frac{1}{2} χ^{2}}}{2^{\frac{ν}{2}} Γ(ν/2)} .

This distribution has \overline{χ^{2}} = ν and σ_{χ^{2}} = \sqrt{2ν}. For large ν, the distribution of χ^{2}/ν is itself approximately Gaussian with mean 1 and standard deviation \sqrt{2/ν}, so that one can estimate the confidence level from the number of standard deviations χ^{2}/ν exceeds 1.

4. Fitting with Correlated Errors

It frequently happens that the data measurements y_{l} are not independent of each other, so that movement of a particular measurement will cause the simultaneous adjustment of other measurements. A particularly odious example of this is multiple scattering, in which a single scattering event at a particular point affects the drift distance measurements at all subsequent points. In this case the deviation of the measurements from the original track will have an independent component due to random measurement errors and a correlated component due to the presence of the multiple scattering event.

Correlated data measurements can be accounted for in a simple way in the least squares algorithm. The χ^{2} function, obtained from section 2, can be written as

χ^{2} = (Δy − Aη)^{t}V_{y}^{−1}(Δy − Aη),

where the “weight” matrix V_{y}^{−1} should now be considered to be the inverse of the covariance matrix V_{y} = ⟨δyδy^{t}⟩, and reduces to the familiar diagonal form for independent measurements. The Gauss-Markov theorem guarantees that the parameters η obtained by minimizing this χ^{2} function are both unbiased and have minimum variance, just as in the case of independent measurements. The solution for the parameters η = α − α_{A} and their covariance matrix is then exactly the same as before, η = V_{A}A^{t}V_{y}^{−1}Δy, V_{η} = V_{A}. The only difficulty is practicality: it is frequently too time consuming to invert V_{y} to get the matrix V_{y}^{−1}. In the third paper in this series I will show a way around this difficulty.

5. Fitting with Constraints

Many times it is useful to fit a set of data points to parameters whose values must satisfy constraint equations. In general, one uses constraints to impose “external knowledge”of a physical process to govern the behavior of the fit, thereby forcing the fit to conform to physical principles which are not known by the internal variables of the fit.

For example, forcing 2 tracks to come from a common vertex or to be back- to-back are common uses of constraints. In the first case, the external knowledge is the fact that the tracks had to emerge from a single space-time point (disregarding secondary vertices). In the second, knowledge of the kinematics of two body decays constrains the behavior of the fit.

When constraints are applied, the effective number of unknowns in the fit is reduced by the number of constraints. This can be seen in a trivial way by writing the chisquare equation (to be minimized) with m parameters together with the r constraint equations. By substituting the r constraints in the χ^{2} equation one is left with an expression having m − r unknowns. In the second example above, for instance, the effective number of parameters is reduced from 10 to 10 − 5 = 5 because of the 5 back-to-back constraints.

Although the elimination of unknowns by direct substitution is popular and fast, there are certain tradeoffs. First of all, the substitution can only be performed when one has relatively simple constraints. In many cases the substitution is impossible because the equations cannot be solved analytically for the variables to be substituted. Secondly, the resulting covariance matrix is calculated only for the reduced set of variables and it is sometimes difficult to extend it to the more meaningful set of original parameters by propagation of errors. In short, the somewhat reduced execution times of substitution algorithms is offset by the special computer code that must be written for each case. Moreover, the addition of new constraints is accompanied by even more program rewrites, which inevitably result in lost time and increased errors.

The method of Lagrange multipliers solves these problems by the rather counterintuitive notion of introducing more variables into the fit. Let the r constraint equations be written H(α) = 0. Expanding about α_{A} as before, we obtain the linearized equations H(α_{A}) + (α − α_{A})∂H(α_{A})/∂α ≡ Dη + d = 0, using obvious notation, where D is a r × m matrix and d is a vector of length r. The constraints are imposed by adding a new term to the χ^{2} equation,

χ^{2} = (Δy − Aη)^{t}V_{y}^{−1}(Δy − Aη) + 2λ^{t}(Dη + d),

where λ, a vector of length r, is the set of Lagrange multipliers. The solution for η is obtained by setting the partial derivatives with respect to η and λ to 0. The ∂/∂λ = 0 equation generates the constraint conditions. The solution to the equations shows that η is equal to η_{0} plus a term proportional to λ, i.e. the constraints “pull” the parameters η away from their unconstrained values η_{0}. This result suggests that the solution can be “factored” into two pieces: (1) solving the unconstrained equations for η_{0} and (2) applying the constraints to solve for η in terms of η_{0}. Since all the uncertainties in the data measurements have been absorbed into V_{η_{0}}, we don’t expect to see any explicit reference to A or V_{A} when constraints are applied. Our intuition is correct as can be seen from the fact that the full χ^{2} can be written (after a little manipulation)

χ^{2} = (Δy − Aη_{0})^{t}V_{y}^{−1}(Δy − Aη_{0}) + (η − η_{0})^{t}V_{η_{0}}^{−1}(η − η_{0}) + λ^{t}(Dη + d).

The first term is constant and is just the χ^{2} for the initial unconstrained fit. Our preferred method then is to perform the fit initially without constraints and then add the constraints after the fact using the new χ^{2}

χ^{2} = (η − η_{0})^{t}V_{η_{0}}^{−1}(η − η_{0}) + 2λ^{t}(Dη + d).

The minimization conditions ∂χ^{2}/∂η = 0 and ∂χ^{2}/∂λ = 0 yield the equations

\begin{aligned}V_{η_{0}}^{−1}(η − η_{0}) + D^{t}λ & = 0 \\ Dη + d & = 0, \end{aligned}

which, when solved, yield

λ = V_{D}(Dη_{0} + d) ≡ V_{D}α  (α = Dη_{0} + d)
η = η_{0} − V_{η_{0}}D^{t}λ
V_{D} = (DV_{η_{0}}D^{t})^{−1}
V_{λ} = V_{D}
V_{η} = V_{η_{0}} − V_{η_{0}}D^{t}V_{λ}DV_{η_{0}}.
χ^{2} = λ^{t}V_{D}^{−1}λ = α^{t}V_{D}α = λ^{t}α = α^{t}λ

The expression for the χ^{2} shows quite clearly that there are r degrees of freedom, one per constraint.

This factorization process of applying constraints to previously fitted parameters (after all the uncertainities have been absorbed in the parameter covariance matrix) will be used frequently in this series of papers. I will show later in this document that one can even factorize the constraints in the sense that the parameters fitted using one set of constraints can be used as input to a second set to obtain the same overall solution as would have been obtained from using the constraints all at once. This has application to kinematic fitting.

6. Solving for Unknown Parameters in the Constraint Equations

1. Direct solution The constrained fitting technique derived in the previous section can be generalized to

the case in which the constraint equations involve new unknown parameters which must be calculated in addition to those that appear directly in the χ^{2} function. A common example is the fitting of several tracks under the constraint that they intersect at an unknown space point.

Instead of applying the Lagrange multiplier technique directly on the original χ^{2} function, we shall use the two step procedure developed in the last section. We first perform the fit without constraints and get the solution η_{0} = V_{A}A^{t}V_{y}^{−1}Δy with the covariance matrix V_{η_{0}} = V_{A}. The new χ^{2} equation written in terms of η and including the constraints is then

χ^{2} = (η − η_{0})^{t}V_{η_{0}}^{−1}(η − η_{0}) + 2λ^{t}(Dη + Ez + d),

where z contains the undetermined parameters (q of them) and E is a r × q matrix. Minimizing this function with respect to η, λ and z yields the equations

\begin{aligned}V_{η_{0}}^{−1}(η − η_{0}) + D^{t}λ & = 0 \\ Dη + Ez + d & = 0 \\ E^{t}λ & = 0. \end{aligned}

Solving these equations is somewhat more tedious than in the previous section since more matrix manipulation is involved. It can be written

λ_{0} = V_{D}(Dη_{0} + d) ≡ V_{D}α_{0}  (α_{0} = Dη_{0} + d)
z = −V_{E}E^{t}λ_{0}
λ = λ_{0} + V_{D}Ez = V_{D}(Dη_{0} + Ez + d) ≡ V_{D}α
η = η_{0} − V_{η_{0}}D^{t}λ
χ^{2} = λ^{t}V_{D}^{−1}λ = λ_{0}^{t}V_{D}^{−1}λ_{0} − z^{t}V_{z}^{−1}z
= α^{t}V_{λ}α = λ^{t}α = α^{t}λ = λ_{0}^{t}α = λ^{t}α_{0}

where the auxiliary matrices are defined as follows

V_{D} = (DV_{η_{0}}D^{t})^{−1}
V_{z} = V_{E} = (E^{t}V_{D}E)^{−1}
V_{λ} = V_{D} − V_{D}EV_{E}E^{t}V_{D}
V_{η} = V_{η_{0}} − V_{η_{0}}D^{t}V_{λ}DV_{η_{0}}
cov(z, η) = −V_{E}E^{t}V_{D}DV_{η_{0}}

The operations should be performed in the following order (assuming that η_{0} and V_{η_{0}} are available from a previous unconstrained fit):

1. Calculate V_{D} and V_{E} 2. Calculate λ_{0}

3. Calculate z

4. Calculate λ

5. Calculate η The number of degrees of freedom is clearly r − q, the number of constraints minus the

number of unknown parameters which have to be determined.

2. Huge error method There is a clear parallel between the above equations and those derived in the previous

section without the new variables z. The only difference is that z have no prior covariance matrix and so must be determined by a special procedure. However, it is clear that we do not change the problem by adding to the χ^{2} the term z^{t}L^{−1}z, where L corresponds to huge errors in z. Since both η and z have a prior covariance matrix, the simple algorithm derived in the last section can be used to solve for η and z. This must give the same answer as the method used in this section in the limit L → ∞.

We introduce the “huge error” algorithm by adding a z covariance term (corresponding to huge errors) to the χ^{2} and then rewrite the χ^{2} in terms of new variables w:

χ^{2} = (η − η_{0})^{t}V_{η_{0}}^{−1}(η − η_{0}) + z^{t}L^{−1}z + 2λ^{t}(Dη + Ez + d)
≡ (w − w_{0})^{t}V_{w_{0}}^{−1}(w − w_{0}) + 2λ^{t}(D′w + d),

where

w = ( \begin{matrix}η \\ z \end{matrix} ) , D′ = ( D E ) and V_{w_{0}} = ( \begin{matrix}V_{η_{0}} & 0 \\ 0 & L \end{matrix} ) ,

and w_{0} = ( \begin{matrix}η_{0} \\ 0 \end{matrix} ). Since the expanded variables w have an initial covariance matrix the solution for w and V_{w} can be taken straight from Section 5:

λ′ = V_{D′}(D′η_{0} + d) ≡ V_{D′}α_{0}  (α_{0} = Dη_{0} + d)
w = w_{0} − V_{w_{0}}D^{′t}λ′
V_{D′} = (D′V_{w_{0}}D^{′t})^{−1} = (DV_{D}D^{t} + ELE^{t})^{−1}
V_{w} = V_{w_{0}} − V_{w_{0}}D^{′t}V_{D′}D′V_{w_{0}}
χ^{2} = λ^{′t}V_{D′}^{−1}λ′ = α_{0}^{t}V_{D′}α_{0} = λ^{′t}α_{0} = α_{0}^{t}λ′.

The solutions for η, z, V_{η}, cov(z, η) and V_{z} can either be extracted from the full solution above using

w = ( \begin{matrix}η \\ z \end{matrix} ) , V_{w} = ( \begin{matrix}V_{η} & cov(z, η) \\ cov(η, z) & V_{z} \end{matrix} ) ,

or by direct calculation:

η = η_{0} − V_{η_{0}}D^{t}λ′
z = −LE^{t}λ′
V_{η} = V_{η_{0}} − V_{η_{0}}D^{t}V_{D′}DV_{η_{0}}
V_{z} = L − LE^{t}V_{D′}EL
cov(z, η) = −LE^{t}V_{D′}DV_{η_{0}}.

It can be shown by expanding V_{D′} using the Woodbury inversion algorithm (see Appendix) that the huge error method in the limit L → ∞ gives exactly the same solution as the direct technique. Moreover, it is simpler to implement on a computer and it reuses the same algorithm derived for the case where no new variables were introduced in the constraints. One must be careful not to make the errors L too large, because the parameters and their covariance matrix are obtained by subtracting very large numbers, causing a significant loss of precision. One should almost certainly use double precision when applyng this method.

7. Factorization of Constraints

Earlier in this paper I showed that when constraints are applied to a fit, it is possible to factor the problem into two parts: (1) obtain the unconstrained parameters by solving the fit without constraints and (2) apply the constraints to these unconstrained parameters to get the final parameters. In this section I will show that the idea of factorization can be generalized so that constraints themselves can be applied sequentially.

The theorem to be proved can be stated as follows. Suppose we have a set of m parameters η_{0} with covariance matrix V_{η_{0}} and a set of r constraints that we want to apply (see Section 5). We divide the r constraints into disjoint subsets containing r_{1}, r_{2},… r_{n} elements. If each of these constraint subsets are applied sequentially, using the parameters and covariance matrix of the previous step and generating a new set of parameters and covariance matrix for the next step, then the final set of parameters η and its covariance matrix V_{η} is the same as would have been obtained by applying simultaneously all r = r_{1} + r_{2} + ⋯ + r_{n} constraints to the initial parameters.

In a sense, this theorem must be true because the final answer shouldn’t depend on the fact that we applied the constraints sequentially rather than simultaneously. The proof is somewhat interesting, however, so I’ll go through it.

First, I’ll divide the constraints into two pieces and apply them sequentially. The proof for two pieces is sufficient since each piece can be subdivided as needed. Write the χ^{2} as (see Section 5)

χ^{2} = (η − η_{0})^{t}V_{η_{0}}^{−1}(η − η_{0}) + 2λ^{t}(Dη + d).

where

D = \binom{D_{1}}{D_{2}} , d = \binom{d_{1}}{d_{2}} and λ = \binom{λ_{1}}{λ_{2}} .

Let’s first find the intermediate parameters η_{1} by applying the constraints D_{1}η_{1} + d_{1} = 0 to

χ_{1}^{2} = (η − η_{0})^{t}V_{η_{0}}^{−1}(η − η_{0}) + 2λ^{′t}(D_{1}η + d_{1}).

The solution from Section 5 is

λ′_{1} = V_{D_{1}}(D_{1}η_{0} + d_{1}) ≡ V_{D_{1}}α_{1}
η_{1} = η_{0} − V_{η_{0}}D_{1}^{t}λ′_{1}
V_{D_{1}} = (D_{1}V_{η_{0}}D_{1}^{t})^{−1}
V_{η_{1}} = V_{η_{0}} − V_{η_{0}}D_{1}^{t}V_{D_{1}}D_{1}V_{η_{0}}
χ_{1}^{2} = λ_{1}^{′t}V_{D_{1}}^{−1}λ′_{1} = λ_{1}^{′t}α_{1}.

The parameters η_{1} and covariance matrix V_{η_{1}} are then used as input for the constraints D_{2}η_{2} + d_{2} = 0. The solution is similar to the one above, i.e.,

λ′_{2} = V_{D_{2}}(D_{2}η_{1} + d_{2}) ≡ V_{D_{2}}α′_{2}
α′_{2} = α_{2} − V_{η_{0}}D_{1}^{t}V_{D_{1}}α_{1}  (α_{2} = D_{2}η_{0} + d_{2})
η_{2} = η_{1} − V_{η_{1}}D_{2}^{t}λ′_{2},
V_{D_{2}} = (D_{2}V_{η_{1}}D_{2}^{t})^{−1}
V_{η_{2}} = V_{η_{1}} − V_{η_{1}}D_{2}^{t}V_{D_{2}}D_{2}V_{η_{1}}
χ_{2}^{2} = λ_{2}^{′t}V_{D_{2}}^{−1}λ′_{2} = λ_{2}^{′t}α′_{2}.

Using the above two sets of equations we can express the solution for η_{2} directly in terms of η_{0} and V_{η_{0}}:

η_{2} = η_{0} − V_{η_{0}}D_{1}^{t}V_{D_{1}}α_{1} − (V_{η_{0}} − V_{η_{0}}D_{1}^{t}V_{D_{1}}D_{1}V_{η_{0}})D_{2}^{t}V_{D_{2}}(α_{2} − V_{η_{0}}D_{1}^{t}V_{D_{1}}α_{1})
χ^{2} = χ_{1}^{2} + χ_{2}^{2},

or, using the vector α = \binom{α_{1}}{α_{2}}

η_{2} = η_{0} −

V_{η_{0}} ( D_{1}^{t} D_{2}^{t}) ( \begin{matrix}V_{D_{1}} + V_{D_{1}}D_{1}V_{η_{0}}D_{2}^{t}V_{D_{2}}V_{η_{0}}D_{1}^{t}V_{D_{1}} & −V_{D_{1}}D_{1}V_{η_{0}}D_{2}^{t}V_{D_{2}} \\ −V_{D_{2}}D_{2}V_{η_{0}}D_{1}^{t}V_{D_{1}} & V_{D_{2}} \end{matrix} ) ( \begin{matrix}α_{1} \\ α_{2} \end{matrix} )

Compare this equation to the one derived in Section 5 for the constraints applied simultaneously. The solution can be written as

η = η_{0} − V_{η_{0}}D^{t}V_{D}α

= η_{0} − V_{η_{0}} ( D_{1}^{t} D_{2}^{t}) ( \begin{matrix}D_{1}V_{η_{0}}D_{1}^{t} & D_{1}V_{η_{0}}D_{2}^{t} \\ D_{2}V_{η_{0}}D_{1}^{t} & D_{2}V_{η_{0}}D_{2}^{t} \end{matrix} )^{−1} ( \begin{matrix}α_{1} \\ α_{2} \end{matrix} ) .

However, the inverse of the matrix in parentheses is just the large matrix in the sequential solution, as can be verified by the factorization matrix inversion algorithm discussed in the Appendix. This proof even works when there are unknown parameters z in the constraint equations, since, as we showed in Section 6, the presence of new parameters can be cast into a form without the variables using an appropriate limit process.

The fact that constraints can be applied one after the other is very useful for kinematic fitting. Consider the decay sequence

\begin{aligned}B^{−} & → D^{∗0}π^{−}, \\ D^{∗0} & → D^{0}π^{0}, \\ D^{0} & → K^{−}π^{+}, \\ π^{0} & → γγ, \end{aligned}

and assume that we want to apply mass constraints to the B^{−}, D^{∗0}, D^{0} and π^{0} and the beam energy constraint to the B^{−}. The four decay processes can be fit simultaneously to give an overall χ^{2}. A more sensible approach, however, would be to fit, in order, the π^{0}, D^{0} and D^{∗0} and then use these particles in the fit for the B^{−} with the remaining mass and beam energy constraint. The final χ^{2} is just the sum of the χ^{2}s for the individual fits. This allows us to reuse the particles in other fits and to develop standard procedures for fitting lower level particles.

8. Simple Application of Unconstrained Track Fitting

In this section I will apply some of the ideas discussed so far to a “toy” track fitting problem. Consider a set of n planar drift chambers arrayed along the x axis in a magnetic field. If the momentum is not too small, the equation of motion is y = a + bx + cx^{2}, where a, b and c are, respectively, the offset, slope and half the curvature of the track. The equation of motion is linear in these 3 parameters.

Since the problem is linear, we can set the parameters to 0 initially (α_{A} = 0, f_{A} = 0) because convergence is guaranteed in a single iteration. The solution is given by η = V_{A}A^{t}V_{y}^{−1}y. The matrices are calculated as follows.

A = ( \begin{matrix}1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ ⋮ & ⋮ & ⋮ \\ 1 & x_{n} & x_{n}^{2} \end{matrix} )

η = ( \begin{matrix}a \\ b \\ c \end{matrix} ) , y = ( y_{1} ; y_{2} ; ⋮ ; y_{n} )

V_{y}^{−1} = ( \begin{matrix}1/σ_{1}^{2} & 0 & ⋯ & 0 \\ 0 & 1/σ_{2}^{2} & ⋯ & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & ⋯ & 1/σ_{n}^{2} \end{matrix} )

V_{A} = (A^{t}V_{y}^{−1}A)^{−1} = ( \begin{matrix}∑_{l} 1/σ_{l}^{2} & ∑_{l} x_{l}/σ_{l}^{2} & ∑_{l} x_{l}^{2}/σ_{l}^{2} \\ ∑_{l} x_{l}/σ_{l}^{2} & ∑_{l} x_{l}^{2}/σ_{l}^{2} & ∑_{l} x_{l}^{3}/σ_{l}^{2} \\ ∑_{l} x_{l}^{2}/σ_{l}^{2} & ∑_{l} x_{l}^{3}/σ_{l}^{2} & ∑_{l} x_{l}^{4}/σ_{l}^{2} \end{matrix} )^{−1}

η = V_{A} ( \begin{matrix}1 & 1 & ⋯ & 1 \\ x_{1} & x_{2} & ⋯ & x_{n} \\ x_{1}^{2} & x_{2}^{2} & ⋯ & x_{n}^{2} \end{matrix} ) ( y_{1}/σ_{1}^{2} ; y_{2}/σ_{2}^{2} ; ⋮ ; y_{n}/σ_{n}^{2} ) = V_{A} ( \begin{matrix}∑_{l} y_{1}/σ_{l}^{2} \\ ∑_{l} x_{l}y_{l}/σ_{l}^{2} \\ ∑_{l} x_{l}^{2}y_{l}/σ_{l}^{2} \end{matrix} )

To get an idea what the covariance matrix looks like, let’s assume further that the measurement errors are identical (σ_{i} = σ) and that the n points are evenly distributed on the interval 0 ≤ x ≤ L. We will also assume that n is large so that the difference between n, n + 1, n + 2, etc. can be ignored. In the approximation of large n, the sums in V_{A}^{−1} can be approximately evaluated as integrals to yield

V_{A} = \frac{σ^{2}}{n} ( \begin{matrix}1 & \frac{1}{2}L & \frac{1}{3}L^{2} \\ \frac{1}{2}L & \frac{1}{3}L^{2} & \frac{1}{4}L^{3} \\ \frac{1}{3}L^{2} & \frac{1}{4}L^{3} & \frac{1}{5}L^{4} \end{matrix} )^{−1} = \frac{σ^{2}}{n} ( \begin{matrix}9 & −36/L & 30/L^{2} \\ −36/L & 192/L^{2} & −180/L^{3} \\ 30/L^{2} & −180/L^{3} & 180/L^{4} \end{matrix} ) .

In this case the errors can be obtained from the diagonal portion of V_{A} to yield:

\begin{aligned}σ_{a} & = \sqrt{\frac{9}{n}}σ \\ σ_{b} & = \sqrt{\frac{192}{n}} \frac{σ}{L} \\ σ_{c} & = \sqrt{\frac{180}{n}} \frac{σ}{L^{2}} . \end{aligned}

These results were first obtained by Gluckstern in 1964 using a somewhat different approach.

9. Simple Application of Constrained Track Fitting

Let’s take the simple model from the last section and apply constraints to it. An interesting set of constraints is to require that two tracks emerge back to back from a common interaction point. The subscripts 1 and 2 will denote the track being considered. Then η, the vector of unknowns, is 6 dimensional and has the form

η = ( \begin{matrix}a_{1} \\ b_{1} \\ c_{1} \\ a_{2} \\ b_{2} \\ c_{2} \end{matrix} )

where we want to enforce the constraints a_{2} = a_{1}, b_{2} = b_{1} and c_{2} = c_{1}. We demand that n equally spaced measurements be taken of track 1 in the interval 0 ≤ x ≤ L and the same number for track 2 be taken for −L ≤ x ≤ 0 (that’s right, folks, we’re faking a back to back 2 track event in a solenoidal detector).

The covariance matrix for the two tracks is initially block diagonal because the tracks were fit independently. We can write it as

V_{η_{0}} = \frac{σ^{2}}{n} ( \begin{matrix}9 & −36/L & 30/L^{2} & 0 & 0 & 0 \\ −36/L & 192/L^{2} & −180/L^{3} & 0 & 0 & 0 \\ 30/L^{2} & −180/L^{3} & 180/L^{4} & 0 & 0 & 0 \\ 0 & 0 & 0 & 9 & 36/L & 30/L^{2} \\ 0 & 0 & 0 & 36/L & 192/L^{2} & 180/L^{3} \\ 0 & 0 & 0 & 30/L^{2} & 180/L^{3} & 180/L^{4} \end{matrix} ) .

Some of the off-diagonal terms for the second track switch sign because the matrix value is summed over odd powers of x for these values.

The constraint equations are written Dη + d = 0. By inspection we see that d = 0 and

D = ( \begin{matrix}1 & 0 & 0 & −1 & 0 & 0 \\ 0 & 1 & 0 & 0 & −1 & 0 \\ 0 & 0 & 1 & 0 & 0 & −1 \end{matrix} ) .

We can compute V_{D} as

\begin{aligned}V_{λ} & = V_{D} = (DV_{η_{0}}D^{t})^{−1} \\ & = \frac{n}{σ^{2}} ( \begin{matrix}18 & 0 & 60/L^{2} \\ 0 & 384/L^{2} & 0 \\ 60/L^{2} & 0 & 360/L^{4} \end{matrix} )^{−1} = \frac{n}{σ^{2}} ( \begin{matrix}1/8 & 0 & −L^{2}/48 \\ 0 & L^{2}/384 & 0 \\ −L^{2}/48 & 0 & L^{4}/160 \end{matrix} ) \end{aligned}

The new track variables are computed, using the matrices we just determined, from

\begin{aligned}λ & = V_{D}(Dη_{0} + d) = V_{D} ( \begin{matrix}a_{1_{0}} − a_{2_{0}} \\ b_{1_{0}} − b_{2_{0}} \\ c_{1_{0}} − c_{2_{0}} \end{matrix} ) = \frac{n}{σ^{2}} ( \begin{matrix}(a_{1_{0}} − a_{2_{0}})/8 − (c_{1_{0}} − c_{2_{0}})L^{2}/48 \\ (b_{1_{0}} − b_{2_{0}})L^{2}/384 \\ −(a_{1_{0}} − a_{2_{0}})L^{2}/48 + (c_{1_{0}} − c_{2_{0}})L^{4}/160 \end{matrix} ) \\ η & = η_{0} − V_{η_{0}}D^{t}λ. \end{aligned}

Finally, we can compute the full covariance matrix of the new variables as

\begin{aligned}V_{η} & = V_{η_{0}} − V_{η_{0}}D^{t}V_{λ}DV_{η_{0}} \\ & = \frac{σ^{2}}{n} ( \begin{matrix}9/8 & 0 & −15/8L^{2} & 9/8 & 0 & −15/8L^{2} \\ 0 & 3/2L^{2} & 0 & 0 & 3/2L^{2} & 0 \\ −15/8L^{2} & 0 & 45/8L^{4} & −15/8L^{2} & 0 & 45/8L^{4} \\ 9/8 & 0 & −15/8L^{2} & 9/8 & 0 & −15/8L^{2} \\ 0 & 3/2L^{2} & 0 & 0 & 3/2L^{2} & 0 \\ −15/8L^{2} & 0 & 45/8L^{4} & −15/8L^{2} & 0 & 45/8L^{4} \end{matrix} ) \end{aligned}

This covariance matrix (which took me 2 hours to compute, so I hope someone is reading this) exhibits all the properties we would expect from a constraint in which track parameters are forced to be equal: (1) The correlation between any two parameters in the same track is the same as that between the same two parameters in the other track or between both tracks; (2) The curvature variance (parameter 3) is 32 times smaller than it was for the unconstrained fit. This makes sense because we can think of the 2 tracks forming a single long track with twice as many measurements and twice the length of a single track. The factor of 1/32 comes from the fact that the variance for the curvature scales like 1/nL^{4}. This same argument cannot be applied to the other 2 parameters because their accuracy depends on where the measurements take place, i.e., in the unconstrained case they are determined externally to the measurement region whereas when constraints are applied they are determined within the measurement region of the “single” track.

Appendix: Tricks for Inverting Large Matrices

1. Factorization method Many square matrices are “almost” invertible in the sense that they can be written as

the sum of a large invertible piece plus a piece that spoils the inversion. An example which comes up frequently is the band matrix

( \begin{matrix}xxx & 0 & 0 & ⋯ & yyy \\ 0 & xxx & 0 & ⋯ & yyy \\ 0 & 0 & xxx & ⋯ & yyy \\ ⋮ & ⋮ & ⋮ & ... & ⋮ \\ yyy & yyy & yyy & ⋯ & zzz \end{matrix} )

where “xxx”, “yyy” and “zzz” are matrices (xxx and zzz are square). Although there are no doubt published ways of inverting such matrices, I came up with an exact algorithm which solves the problem neatly.

Let’s consider the general case first, i.e. not assume any symmetry. We write the matrix we want to invert and its inverse as

V^{−1} = ( \begin{matrix}V_{11} & V_{12} \\ V_{21} & V_{22} \end{matrix} )^{−1} = ( \begin{matrix}A & C \\ B & D \end{matrix} )

where V_{11} and V_{22} are square matrices and assume that V_{11} is invertible. Then the equation

( \begin{matrix}A & C \\ B & D \end{matrix} ) ( \begin{matrix}V_{11} & V_{12} \\ V_{21} & V_{22} \end{matrix} ) = ( \begin{matrix}1 & 0 \\ 0 & 1 \end{matrix} )

must be satisfied. Solving the implied four equations for A, B, C and D yields

V^{−1} = ( \begin{matrix}V_{11} & V_{12} \\ V_{21} & V_{22} \end{matrix} )^{−1} = ( \begin{matrix}V_{11}^{−1} + V_{11}^{−1}V_{12}S^{−1}V_{21}V_{11}^{−1} & −V_{11}^{−1}V_{12}S^{−1} \\ −S^{−1}V_{21}V_{11}^{−1} & S^{−1} \end{matrix} )

where S = V_{22} − V_{21}V_{11}^{−1}V_{12}. Since the only matrix that must be inverted is S (V_{11}^{−1} is already known), the dimensions of V_{22}, and hence S, should be chosen to be as small possible for maximum effectiveness. If V is symmetric then V_{21} = V_{12}^{t}.

If the matrix is such that only the element V_{22} can be inverted, the sequence of steps that led top the previous solution can be modified to give

V^{−1} = ( \begin{matrix}V_{11} & V_{12} \\ V_{21} & V_{22} \end{matrix} )^{−1} = ( \begin{matrix}S^{′−1} & −S^{′−1}V_{12}V_{22}^{−1} \\ −V_{22}^{−1}V_{21}S^{′−1} & V_{22}^{−1} + V_{22}^{−1}V_{21}S^{′−1}V_{12}V_{22}^{−1} \end{matrix} )

where S′ = V_{11} − V_{12}V_{22}^{−1}V_{21}. The same statements about maximum effectiveness apply here: the dimensions of V_{11} should be chosen to be as small as possible.

2. Woodbury formula Instead of factoring the matrix into 4 pieces, it is sometimes possible to write it as

A + UV^{t}, where A is invertible. A is a n × n matrix and U, V are n × p matrices, where p < n. Our strategy, as in the previous section, is to find a way to invert the total matrix by inverting another smaller matrix of size p × p.

The general expansion of the inverse of a sum of matrices A + δA is

(A + δA)^{−1} = A^{−1} − A^{−1}δAA^{−1} + A^{−1}δAA^{−1}δAA^{−1} + …

In our case δA = UV^{t} so the expansion becomes

(A + UV^{t})^{−1} = A^{−1} − A^{−1}UV^{t}A^{−1} + A^{−1}UV^{t}A^{−1}UV^{t}A^{−1} − …

Looking at the expansion we can pick out a repeating unit X = V^{t}A^{−1}U which is of size p × p. The equation now becomes

(A + UV^{t})^{−1} = A^{−1} − A^{−1}UV^{t}A^{−1} + A^{−1}UXV^{t}A^{−1} − A^{−1}UX^{2}V^{t}A^{−1} + …

= A^{−1} − A^{−1}U(1 − X + X^{2} − X^{3} + …)V^{t}A^{−1} = A^{−1} − A^{−1}U(1 + X)^{−1}V^{t}A^{−1} = A^{−1} − A^{−1}U(1 + V^{t}A^{−1}U)^{−1}V^{t}A^{−1}

where (1 + V^{t}A^{−1}U) is a p × p matrix. This expression is called the Woodbury formula. It is stated without proof in Numerical Recipes, p.68.