\input phyzzx
\input tables
\def\A{{\rm\bf A}}
\def\B{{\rm\bf B}}
\def\C{{\rm\bf C}}
\def\D{{\rm\bf D}}
\def\E{{\rm\bf E}}
\def\G{{\rm\bf G}}
\def\H{{\rm\bf H}}
\def\R{{\rm\bf R}}
\def\r{{\rm\bf r}}
\def\V{{\rm\bf V}}
\def\d{{\rm\bf d}}
\def\f{ { {\rm\bf f}_A}}
\def\g{{\rm\bf g_0}}
\def\etab{{\bf \eta}}
\def\x{{\rm\bf x}}
\def\y{{\rm\bf y}}
\def\dy{{\rm\bf \Delta y}}
\def\z{{\rm\bf z}}
\def\w{{\rm\bf w}}
\def\cov{{\rm cov}}
\def\zero{{\rm\bf 0}}
\def\lam{{\bf \lambda}}
\def\chisq{{\chi^2}}

\def\newline{\hfil\break}

\hsize=6.50truein
\vsize=8.75truein	%Leaves enough room for page number
\hoffset=0.00truein	%Actually gives 1 inch
\voffset=0.00truein	%Exactly right

\singlespace
\rightline{\vbox{
  \hbox{Paul Avery}
  \hbox{CBX 91 $-$ 72}
  \hbox{October 18, 1991}
}}

\normalspace
\fourteenpoint
\centerline{\bf Applied Fitting Theory I}
\centerline{\bf General Least Squares Theory}
\twelvepoint
\vskip 0.1 truein
%
\singlespace
{\obeylines
{\bf Section 1:} General Statements
{\bf Section 2:} Simple Least Squares Fitting
{\bf Section 3:} Calculation of Errors: The Covariance Matrix
{\bf Section 4:} Fitting with Correlated Errors
{\bf Section 5:} Fitting with Constraints
{\bf Section 6:} Solving for Unknown Parameters in the Constraint Equations
{\bf Section 7:} Factorization of Constraints
{\bf Section 8:} Simple Application of Unconstrained Track Fitting
{\bf Section 9:} Simple Application of Constrained Track Fitting

{\bf Appendix:}  Tricks for Inverting Large Matrices
}
\normalspace
\parindent 20pt
%
\chapter {\bf 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 Non-optimal 
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.

\chapter {\bf 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 $\alpha_i$, where $m\leq n$. In
other words, we are hypothesizing the relation $y_l = f_l(\alpha_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 $\sigma_l$,
the equation $y_l = f_l(\alpha_i)$ cannot be satisfied exactly if 
$m < n$. However, we can require that the equation be satisfied ``as 
closely as possible'' by defining the $\chisq$ statistic, where
$$\chisq = \sum_l { {(y_l - f_l(\alpha))}^2 \over {\sigma_l^2} }, $$
and demand that the values of $\alpha_i$ be chosen so as to minimize
$\chisq$. Note that $\chisq$ is merely the sum of the squares
of the ``misses'' weighted by the uncertainty $\sigma_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 $\sigma_l$. If the errors are
Gaussian and the expected value for the measurement is $f_l(\alpha)$
then the probability density function (pdf) for the measurements
can be written
$$
   g(y_l) = {1 \over \sqrt{2\pi}\sigma_l}\exp{[-(y_l - f_l)^2/2\sigma_l^2]},
$$
where the normalization is chosen to make
$\int_{-\infty}^{\infty}g_l(y_l)dy_l = 1$. The joint pdf for $n$ measurements
is then
$$\eqalign{
   g(y_1\dots y_n) &= \prod_l{1 \over \sqrt{2\pi}\sigma_l}
            \exp{[-\sum_l(y_l - f_l)^2/2\sigma_l^2]}  \cr
                   &= \prod_l{1 \over \sqrt{2\pi}\sigma_l}
            \exp{(-\chisq/2)}  \cr
}$$
It is clear that maximizing this pdf is equivalent to minimizing $\chisq$.
The adjustment of the parameters $\alpha$ 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 $\chisq$ function be at a minimum is usually 
guaranteed by the $m$ relations
$$ { {\partial \chisq} \over \partial \alpha_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 $\alpha_i$. Often, however, the 
functions vary slowly enough so that one can expand to first order
about an approximate solution, viz.
$$ f_l(\alpha) = f_l(\alpha_A) + \sum_i (\alpha_i - \alpha_{iA})
{ {\partial f_l(\alpha)} \over   \partial\alpha_i}\vrule_{\alpha_A}  
\equiv f_{lA} + \sum_i A_{li}\eta_i. $$
Under this approximation, the $\chisq$ function can be written
$$\eqalign{
\chisq = \sum_l { {(y_l - f_{lA} - \sum_i A_{li}\eta_i)}^2 \over 
        {\sigma_l^2} }
 \equiv \sum_l { {(\Delta y_l - \sum_i A_{li}\eta_i)}^2 \over 
        {\sigma_l^2} } \cr
}.$$
Under very general conditions it can be proved that the parameters 
$\eta_i = \alpha_i - \alpha_{iA}$ obtained by minimizing this 
linearized version of the $\chisq$ function are {\it unbiased} and have
{\it minimum} variance $\sigma_{\alpha_i}^2$, i.e., one cannot find a 
better set of parameters.

The minimization conditions, $\partial \chisq / \partial \alpha_i = 0$,
can be written in this linear approximation as
$$ { \sum_l { {A_{li}A_{lj}} \over \sigma_l^2} }\eta_j = 
 \sum_l { {A_{li}} \over \sigma_l^2} \Delta y_l,$$
where $ 1 \leq i,j \leq m$. Defining the symmetric matrix on the LHS as
$$V_{A\;ij} = \left(\sum_l A_{li}A_{lj}/\sigma_l^2 \right)^{-1},$$
we get the solution
$$ \eta_i = \sum_{j,l} V_{A\;ij} { {A_{lj} \Delta y_l} \over\sigma_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 $\chisq$ equation and write it in
matrix form:
$$ \chisq = (\dy - \A\etab)^t \V_y^{-1} (\dy - \A\etab)$$
where $\A$ is the matrix of coefficients defined above and
$\V_y^{-1}$, $\etab$ and $\y$ are given by the expressions
$$ \V_y^{-1} =
\left( \matrix{
1/\sigma_1^2  &      0       &      \cdots  &      0      \cr
     0        & 1/\sigma_2^2 &      \cdots  &      0      \cr
   \vdots     &    \vdots    &      \ddots  &    \vdots   \cr
     0        &      0       &      \cdots  &    1/\sigma_n^2     \cr
}\right),
\quad \etab =
\left( \matrix{
 \eta_1 \cr
 \eta_2 \cr
 \vdots \cr
 \eta_m \cr
}\right),
\quad \y =
\left( \matrix{
 y_1 \cr
 y_2 \cr
 \vdots \cr
 y_n \cr
}\right).
$$
$\V_y^{-1}$, the inverse of the measurement covariance matrix,
is called the ``weight matrix''.

Taking the partial derivative
with respect to the parameters $\etab$ we obtain the $m$ equations
$$ \A^t \V_y^{-1} (\dy - \A\etab) = \zero,$$
which can be written
$$ \A^t \V_y^{-1} \A \etab = \A^t \V_y^{-1} \dy.$$
Defining as before $\V_A = (\A^t \V_y^{-1} \A)^{-1}$, we finally obtain 
the $m$ equations
$$ \etab = \V_A \A^t \V_y^{-1} \dy,$$
which is a much more understandable and meaningful expression than the 
previous one. The parameters $\alpha$ can be
determined from $\alpha = \alpha_A + \etab$. Note that 
$\V_y^{-1}$ has dimensions of $n\times n$, \A\ is of size $n\times m$, 
and $\V_A$ is of size $m\times 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 = \dy - \A \etab = \y - \f(\alpha_A) - \A \etab$, which are merely
the differences of the measurements and predicted values. The $\chisq$
can then be written
$\chisq = \r^t \V_y^{-1} \r = \sum_l r_l^2/\sigma_l^2$.

\chapter {\bf 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 {\it covariance} matrix,
which is defined as
$$ \cov(x_i,x_j) \equiv V_{x\;ij} = \langle(x_i - \bar x_i) (x_j - \bar x_j)
     \rangle
\equiv \langle\delta x_i \delta x_j\rangle,$$
where the symbol $\langle\rangle$ 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 $\sigma_{x_i}^2$, the square of the standard deviation of
$x_i$. The covariance matrix can be written in matrix notation as
$ \V_{\x} = \langle\delta\x \delta\x^t\rangle$.
For example, the covariance matrix for the (independent) data values $y_l$ is
$$ \V_y = \langle\delta \y \delta \y^t\rangle = {\rm diag}(\sigma_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
{\it correlation} matrix $\rho$. The elements $\rho_{ij}$ are computed from
$$ \rho_{ij} = { {V_{ij}} \over {\sqrt{V_{ii}V_{jj}} } } = 
   { {V_{ij}} \over {\sigma_i\sigma_j} }.$$
The element $\rho_{ij}$ ($i\neq j$) is called the $ij$ correlation coefficient
and satisfies the inequality $-1 \leq \rho_{ij} \leq 1$.

Finding the errors in the parameters $\etab$ is now straightforward. 
Starting with the
solution $ \etab = \V_A \A^t \V_y^{-1} \dy$ we obtain
$$ \V_{\alpha} = \V_{\eta} = 
\langle\delta\etab\delta\etab^t\rangle = 
\V_A \A^t \V_y^{-1} \langle\delta\y\delta\y^t\rangle \V_y^{-1} \A \V_A = 
\V_A \A^t \V_y^{-1} \V_y \V_y^{-1} \A \V_A = \V_A,$$
and the errors are $\sigma_{\eta_i} = \sqrt{V_{\eta\;ii}} = \sqrt{V_{A\;ii}}$.
I have used $\delta \dy = \delta \y$ and the definition
$\V_A = (\A^t \V_y^{-1} \A)^{-1}$. 
The fact that $\V_{\etab} = \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}\etab = \C\y$ and 
gives $\V_{\etab} = \V$.

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

\begintable
\# standard deviations  | Probability to fall in $\pm n \sigma$ \nr
    0.5  |  0.3830 \nr
    1.0  |  0.6826 \nr
    1.5  |  0.8664 \nr
    2.0  |  0.9554 \nr
    2.5  |  0.9876 \nr
    3.0  |  0.9974 \nr
    3.5  |  0.9996
\endtable
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 $\chisq$,
so the final value of $\chisq$ 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 $\chisq$ greater than or equal to the
value we found from the fit. Mathematically, we write
$$
    {\rm C.L.} = \int_{\chisq}^{\infty} f(x) dx,
$$
where $f(\chisq)$ is the $\chisq$ probability distribution function.

For Gaussian errors we can calculate $f(\chisq)$ straightforwardly.
The $\chisq$ can be written in the compact form
$$
\chisq = \dy^t \V_y^{-1} \dy - \etab^t \V_{\etab}^{-1} \etab,
$$
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 $\chisq$ can be written in the manifestly positive
form $\chisq = w_1^2 + w_2^2 + \cdots + w_{\nu}^2$, where the $w_i$
are independent Gaussians having unit error and zero mean.
$\nu = 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{\chisq}$
as the ``radius'' and $\nu-1$ angular variables over which we integrate,
we find the $\chisq$ probability distribution to be
$$
  f(\chisq) = { (\chisq)^{\nu-2\over2} e^{-\half\chisq} \over
     2^{\nu\over 2} \Gamma({\nu/2})}.
$$
This distribution has
$\overline{\chisq} = \nu$ and $\sigma_{\chisq} = \sqrt{2\nu}$. For large
$\nu$, the distribution of $\chisq/\nu$ is itself approximately Gaussian
with mean $1$ and standard deviation $\sqrt{2/\nu}$, so that one can estimate
the confidence level from the number of standard deviations
$\chisq/\nu$ exceeds 1.

\chapter {\bf 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 $\chisq$ function, obtained from section 2,
can be written as
$$ \chisq = (\dy - \A\etab)^t \V_y^{-1} (\dy - \A\etab),$$
where the ``weight'' matrix $\V_y^{-1}$ should now be considered to be 
the inverse of the covariance matrix
$ \V_y = \langle\delta\y \delta\y^t\rangle$, and reduces to the familiar
diagonal form for independent measurements. The Gauss-Markov theorem
guarantees that the parameters $\etab$ obtained by minimizing this 
$\chisq$ function are both unbiased and have minimum variance, just as in
the case of independent measurements. The solution
for the parameters $\etab = \alpha - \alpha_A$ and their covariance
matrix is then exactly the same as before,
$ \etab = \V_A \A^t \V_y^{-1} \dy$,
$\V_{\etab} = \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.

\chapter {\bf 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 $\chisq$ 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(\alpha) = \zero$. Expanding about $\alpha_A$ as before, we
obtain the linearized equations                                                 
$\H(\alpha_A) + (\alpha - \alpha_A) \partial \H(\alpha_A) / \partial\alpha
 \equiv \D\etab + \d = \zero$, using obvious notation, where
$\D$ is a $r \times m$ matrix and $\d$ is a vector of length $r$. 
The constraints are imposed by adding a new term to the $\chisq$ equation,
$$
\chisq = (\dy - \A\etab)^t \V_y^{-1} (\dy - \A\etab) + 2\lam^t(\D\etab + \d),
$$
where $\lam$, a vector of length $r$, is the set of Lagrange multipliers.
The solution for $\etab$ is obtained by setting the partial derivatives with
respect to $\etab$ and $\lam$ to 0. The $\partial/\partial\lam = \zero$
equation generates the constraint conditions.
The solution to the equations shows that $\etab$ is equal to $\etab_0$
plus a term proportional to $\lam$, i.e. the
constraints ``pull'' the parameters $\etab$ away from
their unconstrained values $\etab_0$. This result suggests that
the solution can be ``factored'' into two pieces: (1) solving the
unconstrained equations for $\etab_0$ and (2) applying the constraints to
solve for $\etab$ in terms of $\etab_0$.
Since all the uncertainties in the data measurements have been
absorbed into $\V_{\etab_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 $\chisq$
can be written (after a little manipulation)
$$
 \chisq = (\dy - \A\etab_0)^t \V_y^{-1} (\dy - \A\etab_0) + 
       (\etab - \etab_0)^t \V_{\etab_0}^{-1} (\etab - \etab_0) + 
       \lam^t(\D\etab + \d).
$$
The first term is constant and is just the $\chisq$ for the initial
unconstrained fit. Our preferred method then is to perform the fit initially
{\it without} constraints and then add the constraints after the fact using
the new $\chisq$
$$ \chisq = (\etab - \etab_0)^t \V_{\etab_0}^{-1}
   (\etab - \etab_0) + 2\lam^t (\D \etab + \d).
$$
The minimization conditions $\partial\chisq / \partial\etab = \zero$ and
$\partial\chisq / \partial\lam = \zero$ yield the equations
$$\eqalign{
 \V_{\etab_0}^{-1} (\etab - \etab_0) + \D^t \lam &= \zero    \cr
 \D \etab + \d &= \zero,     \cr
}$$
which, when solved, yield
$$\eqalign{
 \lam &= \V_D (\D\etab_0 + \d) \equiv \V_D\alpha 
                           {\quad\quad (\alpha = \D\etab_0 + \d)}  \cr
 \etab &= \etab_0 - \V_{\etab_0} \D^t \lam    \cr
 \V_D &= (\D \V_{\etab_0} \D^t)^{-1}    \cr
 \V_{\lam} &= \V_D     \cr
 \V_{\etab} &= \V_{\etab_0} - \V_{\etab_0} 
              \D^t \V_{\lam} \D \V_{\etab_0}. \cr
 \chisq &= \lam^t \V_D^{-1} \lam = \alpha^t \V_D \alpha = \lam^t\alpha = 
              \alpha^t\lam \cr
}$$
The expression for the $\chisq$ 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.

\chapter {\bf Solving for Unknown Parameters in the Constraint Equations}
\noindent{\bf 1. Direct solution}\newline
\indent 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 $\chisq$ 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 $\chisq$ function, we shall use the two step procedure developed
in the last section. We first perform the fit without constraints and
get the solution $\etab_0 = \V_A \A^t \V_y^{-1} \dy $ with the
covariance matrix $\V_{\etab_0} = \V_A$. 
The new $\chisq$ equation written in terms of $\etab$ and including the 
constraints is then
$$ \chisq = (\etab - \etab_0)^t \V_{\etab_0}^{-1}
   (\etab - \etab_0) + 2\lam^t (\D \etab + \E\z + \d), $$
where $\z$ contains the undetermined parameters ($q$ of them) 
and $\E$ is a $r\times q$ matrix.
Minimizing this function with respect to $\etab$, $\lam$ and $\z$
yields the equations
$$\eqalign{
 \V_{\etab_0}^{-1} (\etab - \etab_0) + \D^t \lam &= \zero    \cr
 \D \etab + \E\z + \d &= \zero     \cr
 \E^t \lam &= \zero.    \cr
}$$

Solving these equations is somewhat more tedious than in the previous 
section since more matrix manipulation is involved. It can be written
$$\eqalign{
\lam_0 &= \V_D (\D\etab_0 + \d) \equiv \V_D\alpha_0
                           {\quad\quad (\alpha_0 = \D\etab_0 + d)}   \cr
\z &= -\V_E \E^t \lam_0     \cr
\lam &= \lam_0 + \V_D \E \z = \V_D(\D\etab_0 + \E\z + \d) \equiv \V_D\alpha \cr
\etab &= \etab_0 - \V_{\etab_0} \D^t \lam   \cr
\chisq &= \lam^t \V_D^{-1} \lam = \lam_0^t \V_D^{-1} \lam_0 -
                 \z^t \V_{\z}^{-1} \z \cr
       &= \alpha^t \V_{\lam} \alpha  = \lam^t\alpha = \alpha^t\lam
          = \lam_0^t\alpha = \lam^t\alpha_0   \cr
}$$
where the auxiliary matrices are defined as follows
$$\eqalign{
\V_D &= (\D \V_{\etab_0} \D^t)^{-1}      \cr
\V_{\z} &= \V_E = (\E^t \V_D \E)^{-1}      \cr
\V_{\lam} &= \V_D - \V_D \E \V_E \E^t \V_D   \cr
\V_{\etab} &= \V_{\etab_0} - \V_{\etab_0} \D^t 
      \V_{\lam} \D \V_{\etab_0}    \cr
\cov(\z,\etab) &= -\V_E\E^t\V_D\D\V_{\etab_0}  \cr
}$$

The operations should be performed in the following order (assuming that
$\etab_0$ and $\V_{\etab_0}$ are available from a previous
unconstrained fit):
\item{1.\ }Calculate $\V_D$ and $\V_E$
\item{2.\ }Calculate $\lam_0$
\item{3.\ }Calculate $\z$
\item{4.\ }Calculate $\lam$
\item{5.\ }Calculate $\etab$

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.

\noindent{\bf 2. Huge error method}\newline
\indent 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 $\chisq$ the term
$\z^t L^{-1} \z$, where $L$ corresponds
to huge errors in $\z$. Since both $\etab$ and $\z$ have a prior covariance
matrix, the simple algorithm derived in the last section can be used to solve
for $\etab$ and $\z$. This must give the same answer as the method used in
this section in the limit $L \to \infty$.

We introduce the ``huge error'' algorithm by adding a $\z$ covariance term
(corresponding to huge errors) to the $\chisq$ and then rewrite the $\chisq$
in terms of new variables $\w$:
$$\eqalign{
 \chisq &= (\etab - \etab_0)^t \V_{\etab_0}^{-1} (\etab - \etab_0)
     + \z^t L^{-1} \z + 2\lam^t (\D \etab + \E\z + \d)  \cr
    &\equiv (\w - \w_0)^t \V_{\w_0}^{-1} (\w - \w_0)
         + 2\lam^t (\D' \w + \d),  \cr
}$$
where $$
\w = \left( \matrix{\etab \cr \z \cr} \right), \quad
\D' = \left( \matrix{\D & \E} \right) \quad {\rm and} \quad
\V_{\w_0} = \left( \matrix{\V_{\etab_0} & 0  \cr  0 & L \cr} \right),
$$
and $\w_0 = \left(\matrix{\etab_0  \cr  \zero}\right)$. Since the expanded
variables $\w$ have an initial covariance
matrix the solution for $\w$ and $\V_{\w}$ can be taken
straight from Section 5:
$$\eqalign{
 \lam' &= \V_{D'} (\D'\etab_0 + \d) \equiv \V_{D'}\alpha_0
                           {\quad\quad (\alpha_0 = \D\etab_0 + \d)}  \cr
 \w &= \w_0 - \V_{\w_0} \D'^t \lam'   \cr
 \V_{D'} &= (\D' \V_{\w_0} \D'^t)^{-1} = (\D \V_D \D^t + \E L \E^t)^{-1}  \cr
 \V_{\w} &= \V_{\w_0} - \V_{\w_0} \D'^t \V_{D'} \D' \V_{\w_0} \cr
 \chisq &= \lam'^t \V_{D'}^{-1} \lam' = \alpha_0^t \V_{D'} \alpha_0 = 
\lam'^t\alpha_0 = \alpha_0^t\lam'. \cr
}$$

The solutions for $\etab$, $\z$, $\V_{\etab}$, $\cov(\z,\etab)$ and
$\V_{\z}$ can either be extracted from the full solution above using
$$
  \w = \left( \matrix{\etab \cr \z \cr} \right), \quad
  \V_{\w} = \left( \matrix{\V_{\etab}        & \cov(\z,\etab) \cr
                            \cov(\etab,\z)   & \V_{\z}        \cr} \right),
$$
or by direct calculation:
$$\eqalign{
        \etab &= \etab_0 - \V_{\etab_0} \D^t \lam'   \cr
         \z   &= -L \E^t \lam'    \cr
   \V_{\etab} &= \V_{\etab_0} - \V_{\etab_0} D^t \V_{D'} D \V_{\etab_0} \cr
       V_{\z} &= L - L \E^t \V_{D'} \E L  \cr
\cov(\z,\etab) &=  -L \E^t \V_{D'} \D \V_{\etab_0}.  \cr
}$$

It can be shown by expanding $\V_{D'}$ using the Woodbury inversion
algorithm (see Appendix) that the huge error method in the limit
$L \to \infty$ 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.

\chapter {\bf 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 $\etab_0$ with covariance matrix $\V_{\etab_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$,\dots $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 $\etab$ and its covariance matrix $\V_{\etab}$
is the same as would have been obtained by applying simultaneously all
$r=r_1 + r_2 + \cdots + r_n$ constraints to the initial parameters.

In a sense, this theorem {\sl 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 $\chisq$ as (see Section 5)
$$ \chisq = (\etab - \etab_0)^t \V_{\etab_0}^{-1}
   (\etab - \etab_0) + 2\lam^t (\D \etab + \d).
$$
where 
$$\D = \left({\D_1\atop \D_2}\right) ,\quad \d = \left({\d_1 \atop \d_2}\right)
             {\quad\rm and\quad} \lam = \left({\lam_1\atop \lam_2}\right).
$$
Let's first find the intermediate parameters $\etab_1$ by applying
the constraints $\D_1\etab_1 + \d_1 = \zero$ to
$$ \chi^2_1 = (\etab - \etab_0)^t \V_{\etab_0}^{-1}
   (\etab - \etab_0) + 2\lam'^t (\D_1 \etab + \d_1).
$$

The solution from Section 5 is
$$\eqalign{
 \lam'_1 &= \V_{D_1} (\D_1\etab_0 + \d_1) \equiv \V_{D_1}\alpha_1  \cr
 \etab_1 &= \etab_0 - \V_{\etab_0} \D_1^t \lam'_1   \cr
 \V_{D_1} &= (\D_1 \V_{\etab_0} \D_1^t)^{-1}    \cr
  \V_{\etab_1} &= \V_{\etab_0} - \V_{\etab_0} 
              \D_1^t \V_{D_1} \D_1 \V_{\etab_0} \cr
    \chi^2_1 &= \lam'^t_1 \V_{D_1}^{-1} \lam'_1 = \lam'^t_1 \alpha_1. \cr
}$$
The parameters $\etab_1$ and covariance matrix $\V_{\etab_1}$ are
then used as input for the constraints $\D_2\etab_2 + \d_2 = \zero$.
The solution is similar to the one above, i.e.,
$$\eqalign{
 \lam'_2 &= \V_{D_2} (\D_2\etab_1 + \d_2)  \equiv \V_{D_2}\alpha'_2  \cr
    \alpha'_2 &= \alpha_2 - \V_{\etab_0}\D_1^t\V_{D_1}\alpha_1
            {\qquad (\alpha_2 = \D_2\etab_0 + \d_2)}  \cr
 \etab_2 &= \etab_1 - \V_{\etab_1} \D_2^t \lam'_2,  \cr
 \V_{D_2} &= (\D_2 \V_{\etab_1} \D_2^t)^{-1}    \cr
  \V_{\etab_2} &= \V_{\etab_1} - \V_{\etab_1} 
              \D_2^t \V_{D_2} \D_2 \V_{\etab_1} \cr
  \chi^2_2 &= \lam'^t_2 \V_{D_2}^{-1} \lam'_2 = \lam'^t_2 \alpha'_2. \cr
}$$

Using the above two sets of equations we can express the solution
for $\etab_2$ directly in terms of $\etab_0$ and $\V_{\etab_0}$:
$$\eqalign{
\etab_2 &= \etab_0
   - \V_{\etab_0} \D_1^t \V_{D_1}\alpha_1 -
 (\V_{\etab_0} - \V_{\etab_0}\D_1^t \V_{D_1} \D_1 \V_{\etab_0})
  \D_2^t \V_{D_2} (\alpha_2 - \V_{\etab_0}\D_1^t\V_{D_1}\alpha_1)  \cr
   \chisq &= \chi^2_1 + \chi^2_2,  \cr
}$$
or, using the vector
$\alpha = \left({\alpha_1 \atop \alpha_2}\right)$
$$\displaylines{
\etab_2 = \etab_0 \ - \hfill \cr
\V_{\etab_0}
\left(
\matrix{\D_1^t & \D_2^t \cr}
\right)
\left(
\matrix{\V_{D_1} + \V_{D_1}\D_1\V_{\eta_0}\D_2^t
           \V_{D_2}\V_{\etab_0}\D_1^t\V_{D_1}   &
       -\V_{D_1} \D_1 \V_{\etab_0} \D_2^t \V_{D_2}  \cr
       -\V_{D_2} \D_2 \V_{\etab_0} \D_1^t \V_{D_1}  & 
       \V_{D_2}    \cr}
\right)
\left(
\matrix{\alpha_1  \cr  \alpha_2  \cr}
\right)  \cr
}$$

Compare this equation to the one derived
in Section 5 for the constraints applied simultaneously. The
solution can be written as
$$\eqalign{
  \etab &= \etab_0 - \V_{\etab_0} \D^t \V_D \alpha \cr
        &= \etab_0 - \V_{\etab_0}
\left(
\matrix{\D_1^t & \D_2^t  \cr}
\right)
\left(
\matrix{
\D_1\V_{\etab_0}\D_1^t  &  \D_1\V_{\etab_0}\D_2^t  \cr
\D_2\V_{\etab_0}\D_1^t  &  \D_2\V_{\etab_0}\D_2^t  \cr }
\right)^{-1}
\left(
\matrix{\alpha_1  \cr  \alpha_2  \cr}
\right).
}$$
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
$$\eqalign{
    B^-   &\to  D^{*0} \pi^-,   \cr
  D^{*0}  &\to  D^0 \pi^0,      \cr
  D^0     &\to  K^-\pi^+,       \cr
  \pi^0   &\to  \gamma\gamma,   \cr
}$$
and assume that we want to apply mass constraints to the
$B^-$, $D^{*0}$, $D^0$ and
$\pi^0$ and the beam energy constraint to the $B^-$.
The four decay processes can be fit simultaneously to give an overall
$\chisq$. A more sensible approach, however, would be to fit, in order, the
$\pi^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
$\chisq$ is just the sum of the $\chisq$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.

\chapter {\bf 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
($\alpha_A = \zero,\;\f = \zero$) because
convergence is guaranteed in a single iteration. The solution is given
by $\etab = \V_A \A^t \V_y^{-1} \y$. The matrices are calculated
as follows.
$$ \A = \left( \matrix{
  1     &      x_1       &      x_1^2      \cr
  1     &      x_2       &      x_2^2      \cr
\vdots  &    \vdots      &    \vdots       \cr
  1     &      x_n       &      x_n^2      \cr
}\right)$$
$$ \etab = \left( \matrix{
  a \cr
  b \cr
  c \cr
}\right), \qquad  \y = \left( \matrix{
  y_1  \cr
  y_2  \cr
\vdots \cr
  y_n  \cr
}\right)$$
$$ \V_y^{-1} = \left( \matrix{
1/\sigma_1^2  &      0       &      \cdots  &      0      \cr
     0        & 1/\sigma_2^2 &      \cdots  &      0      \cr
   \vdots     &    \vdots    &      \ddots  &    \vdots   \cr
     0        &      0       &      \cdots  &    1/\sigma_n^2     \cr
}\right)$$
$$ \V_A = (\A^t \V_y^{-1} \A)^{-1} = \left( \matrix{
\sum_l 1/\sigma_l^2 & \sum_l x_l/\sigma_l^2 & \sum_l x_l^2/\sigma_l^2  \cr
\sum_l x_l/\sigma_l^2 & \sum_l x_l^2/\sigma_l^2 & \sum_l x_l^3/\sigma_l^2  \cr
\sum_l x_l^2/\sigma_l^2 & \sum_l x_l^3/\sigma_l^2 & \sum_l x_l^4/\sigma_l^2 \cr
}\right)^{-1}$$
$$ \etab = \V_A 
\left( \matrix{
   1    &    1    & \cdots  &   1      \cr
  x_1   &   x_2   & \cdots  &  x_n     \cr
  x_1^2 &   x_2^2 & \cdots  &  x_n^2   \cr
}\right)
\left( \matrix{
   y_1/\sigma_1^2     \cr
   y_2/\sigma_2^2     \cr
      \vdots          \cr
   y_n/\sigma_n^2     \cr
}\right) = \V_A
\left( \matrix{
   \sum_l y_1/\sigma_l^2           \cr
   \sum_l x_l y_l/\sigma_l^2       \cr
   \sum_l x_l^2 y_l/\sigma_l^2     \cr
}\right)$$

To get an idea what the covariance matrix looks like, let's assume
further that the measurement errors are identical ($\sigma_i = \sigma$)
and that the $n$ points are evenly distributed on the interval
$0 \leq x \leq 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 =  {\sigma^2\over n}
\left( \matrix{
    1         &    {1\over 2}L     & {1\over 3}L^2  \cr
{1\over 2}L   &    {1\over 3}L^2   & {1\over 4}L^3  \cr
{1\over 3}L^2 &    {1\over 4}L^3   & {1\over 5}L^4  \cr
}\right)^{-1} = {\sigma^2\over n}
\left( \matrix{
    9   &    -36/L     & 30/L^2   \cr
-36/L   &    192/L^2   &-180/L^3  \cr
 30/L^2 &   -180/L^3   & 180/L^4  \cr
}\right).
$$
In this case the errors can be obtained from the diagonal portion of
$\V_A$ to yield:
$$\eqalign{
            \sigma_a &= \sqrt { {9\over n} } \sigma             \cr
            \sigma_b &= \sqrt { {192\over n} } {\sigma\over L}   \cr
            \sigma_c &= \sqrt { {180\over n} } {\sigma\over L^2}. \cr
} $$
These results were first obtained by Gluckstern in 1964 using a somewhat
different approach.

\chapter {\bf 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 $\etab$, the vector
of unknowns, is 6 dimensional and has the form
$$ \etab = \left( \matrix{
  a_1 \cr
  b_1 \cr
  c_1 \cr
  a_2 \cr
  b_2 \cr
  c_2 \cr
}\right)$$
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 \leq x \leq L$ and the same number for track 2 be taken for
$-L \leq x \leq 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_{\etab_0} = {\sigma^2\over n}
\left( \matrix{
    9   &    -36/L     &  30/L^2  &    0     &     0        &    0      \cr
-36/L   &    192/L^2   &-180/L^3  &    0     &     0        &    0      \cr
 30/L^2 &   -180/L^3   & 180/L^4  &    0     &     0        &    0      \cr
  0     &     0        &    0     &    9     &     36/L     &  30/L^2  \cr
  0     &     0        &    0     &   36/L   &    192/L^2   & 180/L^3  \cr
  0     &     0        &    0     &   30/L^2 &    180/L^3   & 180/L^4  \cr
}\right).      
$$
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\etab + \d = \zero$. By inspection
we see that $\d = \zero$ and
$$ \D = \left( \matrix{
   1  &  0  &  0  & -1  &  0  &  0  \cr
   0  &  1  &  0  &  0  & -1  &  0  \cr
   0  &  0  &  1  &  0  &  0  & -1  \cr
}\right).
$$

We can compute $\V_D$ as 
$$\eqalign{
 \V_{\lam} &= \V_D = (D \V_{\etab_0} \D^t)^{-1}   \cr
           &= {n\over \sigma^2}
\left( \matrix{
    18       &   0       &   60/L^2    \cr
     0       &  384/L^2  &    0        \cr
    60/L^2   &   0       &  360/L^4    \cr
}\right)^{-1} = {n \over \sigma^2}
\left( \matrix{
      1/8    &   0     & -L^2/48   \cr
       0     & L^2/384 &    0      \cr
     -L^2/48 &   0     &  L^4/160  \cr
}\right)
}$$

The new track variables are computed, using the matrices we just determined,
from
$$\eqalign{
 \lam &= \V_D (\D\etab_0 + \d) = \V_D
\left( \matrix{
    a_{1_0} - a_{2_0}    \cr
    b_{1_0} - b_{2_0}    \cr
    c_{1_0} - c_{2_0}    \cr
} \right) =   {n\over \sigma^2}
\left( \matrix{
    (a_{1_0}-a_{2_0})/8 - (c_{1_0}-c_{2_0})L^2/48        \cr
            (b_{1_0}-b_{2_0})L^2/384                     \cr
   - (a_{1_0}-a_{2_0})L^2/48 + (c_{1_0}-c_{2_0})L^4/160  \cr
} \right)                    \cr
 \etab &= \etab_0 - \V_{\etab_0} \D^t \lam.    \cr 
} $$

Finally, we can compute the full covariance matrix of the new variables
as
$$\eqalign{
 \V_{\etab} &= \V_{\etab_0} - \V_{\etab_0} 
              \D^t \V_{\lam} \D \V_{\etab_0}      \cr
   &=  {\sigma^2\over n}
\left( \matrix{
  9/8     &   0    &  -15/8L^2  &    9/8    &     0     &  -15/8L^2 \cr
   0      &  3/2L^2&    0       &     0     &    3/2L^2 &     0     \cr
 -15/8L^2 &   0    &   45/8L^4  &  -15/8L^2 &     0     &   45/8L^4 \cr
  9/8     &   0    &  -15/8L^2  &    9/8    &     0     &  -15/8L^2 \cr
   0      &  3/2L^2&    0       &     0     &    3/2L^2 &     0     \cr
 -15/8L^2 &   0    &   45/8L^4  &  -15/8L^2 &     0     &   45/8L^4 \cr
}\right)      \cr
} $$
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.

\vskip 0.4in
\centerline{\bf\fourteenpoint Appendix: Tricks for Inverting Large Matrices}
\vskip 0.4in
\noindent{\bf 1. Factorization method}\hfil\newline
\indent 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
$$ \left( \matrix{
xxx    &0      &0      &\cdots &yyy    \cr
0      &xxx    &0      &\cdots &yyy    \cr
0      &0      &xxx    &\cdots &yyy    \cr
\vdots &\vdots &\vdots &\ddots &\vdots \cr
yyy    &yyy    &yyy    &\cdots &zzz    \cr
}\right)$$
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} =
\left( \matrix{
V_{11} & V_{12}  \cr
V_{21} & V_{22}  \cr
}\right)^{-1} = 
\left( \matrix{
A & C  \cr
B & D  \cr
}\right)
$$
where $V_{11}$ and $V_{22}$ are square matrices and assume that $V_{11}$
is invertible. Then the equation
$$
\left( \matrix{
A & C  \cr
B & D  \cr
}\right)
\left( \matrix{
V_{11} & V_{12}  \cr
V_{21} & V_{22}  \cr
}\right) =
\left( \matrix{
1 & 0  \cr
0 & 1  \cr
}\right)
$$
must be satisfied. Solving the implied four equations for $A$, $B$, $C$
and $D$ yields
$$
V^{-1} =
\left( \matrix{
V_{11} & V_{12}  \cr
V_{21} & V_{22}  \cr
}\right)^{-1} = 
\left( \matrix{
V_{11}^{-1} +
V_{11}^{-1}V_{12}S^{-1}V_{21}V_{11}^{-1} & -V_{11}^{-1}V_{12}S^{-1}  \cr
        -S^{-1}V_{21}V_{11}^{-1}         & S^{-1}  \cr
}\right)
$$
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} =
\left( \matrix{
V_{11} & V_{12}  \cr
V_{21} & V_{22}  \cr
}\right)^{-1} = 
\left( \matrix{
    S'^{-1}                 &  -S'^{-1}V_{12}V_{22}^{-1}  \cr
-V_{22}^{-1}V_{21}S'^{-1}  &  V_{22}^{-1} + V_{22}^{-1}V_{21}S'^{-1}V_{12}
                                 V_{22}^{-1}  \cr
}\right)
$$
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.

\noindent{\bf 2. Woodbury formula}\hfil\newline
\indent 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\times n$ matrix and $U$, $V$ are $n\times 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\times p$.

The general expansion of the inverse of a sum of matrices $A + \delta A$ is
$$
    (A+\delta A)^{-1} = A^{-1} - A^{-1}\delta A A^{-1} +
           A^{-1}\delta A A^{-1} \delta A A^{-1} + \dots
$$
In our case $\delta 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} - \dots
$$
Looking at the expansion we can pick out a repeating unit
$X = V^tA^{-1}U$ which is of size $p\times p$. The equation now becomes
$$\eqalign{
    (A+UV^t)^{-1} &= A^{-1} - A^{-1}UV^t A^{-1} +
           A^{-1}UXV^t A^{-1} -
           A^{-1}UX^2V^t A^{-1} + \dots  \cr
                  &= A^{-1} - A^{-1}U(1-X+X^2-X^3+\dots)V^t A^{-1}  \cr
                  &= A^{-1} - A^{-1}U(1+X)^{-1}V^t A^{-1}  \cr
                  &= A^{-1} - A^{-1}U(1+V^tA^{-1}U)^{-1}V^t A^{-1}  \cr
}$$
where $(1 + V^t A^{-1} U)$ is a $p\times p$ matrix. This expression is called
the Woodbury formula. It is stated without proof in Numerical Recipes, p.68.
\end
