% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Ordinary least squares}

\section{Introduction}
A system of linear equations is considered {\em
overdetermined}\index{overdetermined system} if there are more equations
than unknown variables. If all equations of an overdetermined system
are linearly independent, the system has no exact solution.

A {\em linear least-squares problem} is the problem of finding an
approximate solution to an overdetermined linear system. It often arises
in applications where a theoretical model is fitted to experimental data.

\section{Linear least-squares problem}
Consider a linear system
	\begin{equation}
\ma{A}\ve{c}=\ve{b}\;,
	\end{equation}
where $\ma{A}$ is an $n\times m$ matrix, $\ve{c}$ is an $m$-component
vector of unknown variables and $\ve{b}$ is an $n$-component vector of
the right-hand side terms. If the number of equations $n$ is larger than
the number of unknowns $m$, the system is overdetermined and generally
has no solution.

However, it is still possible to find an approximate solution---the one
where $\ma{A}\ve{c}$ is only approximately equal $\ve{b}$---in the
sence that the Euclidean norm of the difference between $\ma{A}\ve{c}$
and $\ve{b}$ is minimized,
   \begin{equation}\label{eq:ls:min}
\ve{c}:\:\min_\ve{c}\|\ma{A}\ve{c}-\ve{b}\|^2 \;.
   \end{equation}
The problem (\ref{eq:ls:min}) is called the ordinary least-squares problem
and the vector $\ve{c}$ that minimizes $\|\ma{A}\ve{c}-\ve{b}\|^2$
is called the {\em least-squares solution}.

\section{Solution via QR-decomposition}
The linear least-squares problem can be solved by QR-decomposition. The
matrix $\ma{A}$ is factorized as $\ma{A}=\ma{Q}\ma{R}$, where $\ma{Q}$
is $n\times m$ matrix with orthogonal columns, $\ma{Q}^T\ma{Q}=1$, and
$\ma{R}$ is an $m\times m$ upper triangular matrix. The Euclidean norm
$\|\ma{A} \ve{c} -\ve{b}\|^{2}$ can then be rewritten as
	\begin{eqnarray}
&\|\ma{A} \ve{c} -\ve{b}\|^{2}
= \|\ma{Q}\ma{R} \ve{c} -\ve{b}\|^{2}& \nonumber\\
&= \|\ma{R} \ve{c}
- \ma{Q}^{T}\ve{b}\|^{2} + \|(1 - \ma{Q}\ma{Q}^{T})\ve{b}\|^{2}
\ge  \|(1 - \ma{Q}\ma{Q}^{T})\ve{b}\|^{2} \:.&
	\end{eqnarray}
The term $\|(1-\ma{Q}\ma{Q}^{T})\ve{b}\|^{2}$ is independent of the
variables $\ve{c}$ and can not be reduced by their variations.	However,
the term $\|\ma{R}\ve{c}-\ma{Q}^{T}\ve{b}\|^{2}$ can be reduced down
to zero by solving the $m\times m$ system of linear equations
	\begin{equation}\label{eq:ls:rc}
\ma{R}\ve{c}=\ma{Q}^{T}\ve{b} \;.
	\end{equation}
The system is right-triangular and can be readily solved by
back-substitution.  Thus the solution to the ordinary least-squares
problem (\ref{eq:ls:min}) is given by the solution of the triangular
system (\ref{eq:ls:rc}).

\section{Ordinary least-squares curve fitting}
Ordinary least-squares curve fitting is a problem of fitting $n$
(experimental) data points $\{x_i,y_{i}\pm\Delta y_{i}\}$, where $\Delta
y_{i}$ are experimental errors, by a linear combination, $F$, of $m$
functions $\{f_k(x)\:|\:k=1,\dots,m\}$,
	\begin{equation}\label{ls:eq:F}
F_\ve{c}(x)=\sum_{k=1}^{m}c_{k}f_{k}(x) \;.
	\end{equation}
The objective of the least-squares fit is to minimize the square
deviation, called $\chi^2$, between the fitting function $F(x)$ and the
experimental data,
	\begin{equation}\label{ls:eq:chi}
\chi^{2} = \sum_{i=1}^{n}\left( \frac{F(x_{i})-y_{i}}{\Delta y_{i}}
\right)^{2} \;.
	\end{equation}
Individual deviations from experimental points are weighted with their
inverse errors in order to promote contributions from the more precise
measurements.

Minimization of $\chi^2$ with respect to the coefficiendt $c_k$ in
(\ref{ls:eq:F}) is apparently equivalent to the least-squares problem
(\ref{eq:ls:min}) where
	\begin{equation}
A_{ik} = \frac{f_{k}(x_{i})}{\Delta y_{i}}\;,\; b_{i} =
\frac{y_{i}}{\Delta y_{i}}\;.
	\end{equation}
If $\ma{Q}\ma{R}=\ma{A}$ is the QR-decomposition of the matrix $\ma{A}$,
the formal least-squares solution to the fitting problem is
   \begin{equation}\label{ls:eq:c}
\ve{c} = R^{-1}Q^{T}\ve{b}\;.
   \end{equation}
However, in practice one has to back-substitute the system
$\ma{R}\ve{c}=\ma{Q}^{T}\ve{b}$.

\subsection{Variances and correlations of fitting parameters} Suppose
$\delta y_{i}$ is a small deviation of the measured value of the physical
observable at hand from its exact value. The corresponding deviation
$\delta c_{k}$ of the fitting coefficient is then given as
	\begin{equation}\label{ls:eq:o}
\delta c_{k}=\sum_i\frac{\partial c_{k}}{\partial y_{i}}\delta y_{i}\;.
	\end{equation}
In a good experiment the deviations $\delta y_{i}$ are statistically
independent and distributed normally with the standard deviations $\Delta
y_{i}$. The deviations (\ref{ls:eq:o}) are then also distributed normally
with {\em variances}
	\begin{equation}
\langle\delta c_{k}\delta c_{k}\rangle = \sum_{i}\left(\frac{\partial
c_{k}}{\partial y_{i}}\Delta y_{i}\right)^2 = \sum_{i}\left(\frac{\partial
c_{k}}{\partial b_{i}}\right)^2\;.
	\end{equation}
The standard errors in the fitting coefficients are then given as the
square roots of variances,
	\begin{equation}
\Delta c_{k} = \sqrt{\langle\delta c_{k}\delta c_{k}\rangle}
=\sqrt{\sum_{i}\left(\frac{\partial c_{k}}{\partial b_{i}}\right)^2}\;.
	\end{equation}

The variances are diagonal elements of the {\em covariance matrix},
$\Sigma$, made of {\em covariances},
	\begin{equation}\label{ls:eq:S}
\Sigma_{kq} \equiv \langle\delta c_{k}\delta c_{q}\rangle
=
\sum_{i}
\frac{\partial c_{k}}{\partial b_{i}}
\frac{\partial c_{q}}{\partial b_{i}}\;.
	\end{equation}
Covariances $\langle\delta c_{k}\delta c_{q}\rangle$ are measures of
to what extent the coefficients $c_k$ and $c_q$ change together if the
measured values $y_i$ are varied. The normalized covariances,
	\begin{equation}
\frac{\langle\delta c_{k}\delta c_{q}\rangle} {\sqrt{\langle\delta
c_{k}\delta c_{k}\rangle \langle\delta c_{q}\delta c_{q}\rangle}}
	\end{equation}
are called {\em correlations}.

Using (\ref{ls:eq:S}) and (\ref{ls:eq:c}) the covariance matrix can be
calculated as
	\begin{equation}
\Sigma = \left(\frac{\partial\ve{c}}{\partial\ve{b}}\right)
\left(\frac{\partial\ve{c}}{\partial\ve{b}}\right)^T
= R^{-1} (R^{-1})^T = (R^{T}R)^{-1} = (A^{T}A)^{-1}\;.
	\end{equation}
The square roots of the diagonal elements of this matrix provide the
estimates of the errors $\Delta \ve{c}$ of the fitting coefficients,
	\begin{equation}
\Delta c_k = \sqrt{\Sigma_{kk}} \;\Big|_{k=1\dots m}\;,
	\end{equation}
and the (normalized) off-diagonal elements provide the estimates of their
correlations.

A C implementation of the ordinary least squares fit with QR-decomposition
is shown in Table~(\ref{ls:t:f}).

	\begin{table}[htb]
\caption{C implementation of ordinary least squares fit with
QR-decomposition.}
	\label{ls:t:f}
\lstinputlisting{c/least.squares/lsfit.c}
   \end{table}

An illustration of a fit is shown on Figure~(\ref{ls:f:p}) where a
polynomial is fitted to a set of data.

	\begin{figure}[htb]
\centerline{\includegraphics{c/least.squares/fit.eps}}
\caption{Ordinary least squares fit of $F_\ve{c}(x)=c_1+c_2x+c_3x^2$
to a set of data. Shown are fits with optimal coefficiens $\ve{c}$
as well as with $\ve{c}+\Delta\ve{c}$ and $\ve{c}-\Delta\ve{c}$.}
  \label{ls:f:p}
	\end{figure}

