Total least squares

The total least squares approximation of the data is generically equivalent to the best, in the Frobenius norm, low-rank approximation of the data matrix.

Linear model

Background

In the least squares method of data modeling, the objective function to be minimized, S, is a quadratic form:

${\displaystyle S={\mathbf{r}}^{\mathbf{T}}\mathbf{W}\mathbf{r},}$

where r is the vector of residuals and W is a weighting matrix. In linear least squares the model contains equations which are linear in the parameters appearing in the parameter vector ${\displaystyle \mathit{\beta }}$, so the residuals are given by

${\displaystyle \mathbf{r}\mathbf{=}\mathbf{y}\mathbf{-}\mathbf{X}\mathit{\beta }.}$

There are m observations in y and n parameters in β with m>n. X is a m×n matrix whose elements are either constants or functions of the independent variables, x. The weight matrix W is, ideally, the inverse of the variance-covariance matrix ${\displaystyle {\mathbf{M}}_y}$ of the observations y. The independent variables are assumed to be error-free. The parameter estimates are found by setting the gradient equations to zero, which results in the normal equations

${\displaystyle {\mathbf{X}}^{\mathbf{T}}\mathbf{W}\mathbf{X}\mathit{\beta }\mathbf{=}{\mathbf{X}}^{\mathbf{T}}\mathbf{W}\mathbf{y}.}$

Allowing observation errors in all variables

Now, suppose that both x and y are observed subject to error, with variance-covariance matrices ${\displaystyle {\mathbf{M}}_x}$ and ${\displaystyle {\mathbf{M}}_y}$ respectively. In this case the objective function can be written as

${\displaystyle S={\mathbf{r}}_{\mathbf{x}}^{\mathbf{T}}{\mathbf{M}}_{\mathbf{x}}^{\mathbf{-}\mathbf{1}}{\mathbf{r}}_{\mathbf{x}}\mathbf{+}{\mathbf{r}}_{\mathbf{y}}^{\mathbf{T}}{\mathbf{M}}_{\mathbf{y}}^{\mathbf{-}\mathbf{1}}{\mathbf{r}}_{\mathbf{y}},}$

where ${\displaystyle {\mathbf{r}}_x}$ and ${\displaystyle {\mathbf{r}}_y}$ are the residuals in x and y respectively. Clearly these residuals cannot be independent of each other, but they must be constrained by some kind of relationship. Writing the model function as ${\displaystyle \mathbf{f}\mathbf{(}{\mathbf{r}}_{\mathbf{x}}\mathbf{,}{\mathbf{r}}_{\mathbf{y}}\mathbf{,}\mathit{\beta }\mathbf{)}}$, the constraints are expressed by m condition equations.

${\displaystyle \mathbf{F}\mathbf{=}\mathbf{\Delta }\mathbf{y}\mathbf{-}\frac{\mathrm{\partial }\mathbf{f}}{\mathrm{\partial }{\mathbf{r}}_{\mathbf{x}}}{\mathbf{r}}_{\mathbf{x}}\mathbf{-}\frac{\mathrm{\partial }\mathbf{f}}{\mathrm{\partial }{\mathbf{r}}_{\mathbf{y}}}{\mathbf{r}}_{\mathbf{y}}\mathbf{-}\mathbf{X}\mathbf{\Delta }\mathit{\beta }\mathbf{=}\mathbf{0}.}$

Thus, the problem is to minimize the objective function subject to the m constraints. It is solved by the use of Lagrange multipliers. After some algebraic manipulations, the result is obtained.

${\displaystyle {\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-}\mathbf{1}}\mathbf{X}\mathbf{\Delta }\mathit{\beta }\mathbf{=}{\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-}\mathbf{1}}\mathbf{\Delta }\mathbf{y},}$

or alternatively ${\displaystyle {\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-}\mathbf{1}}\mathbf{X}\mathit{\beta }\mathbf{=}{\mathbf{X}}^{\mathbf{T}}{\mathbf{M}}^{\mathbf{-}\mathbf{1}}\mathbf{y},}$ where M is the variance-covariance matrix relative to both independent and dependent variables.

${\displaystyle \mathbf{M}\mathbf{=}{\mathbf{K}}_{\mathbf{x}}{\mathbf{M}}_{\mathbf{x}}{\mathbf{K}}_{\mathbf{x}}^{\mathbf{T}}\mathbf{+}{\mathbf{K}}_{\mathbf{y}}{\mathbf{M}}_{\mathbf{y}}{\mathbf{K}}_{\mathbf{y}}^{\mathbf{T}}\mathbf{;}{\mathbf{K}}_{\mathbf{x}}\mathbf{=}\mathbf{-}\frac{\mathrm{\partial }\mathbf{f}}{\mathrm{\partial }{\mathbf{r}}_{\mathbf{x}}}\mathbf{,}{\mathbf{K}}_{\mathbf{y}}\mathbf{=}\mathbf{-}\frac{\mathrm{\partial }\mathbf{f}}{\mathrm{\partial }{\mathbf{r}}_{\mathbf{y}}}.}$

Example

When the data errors are uncorrelated, all matrices M and W are diagonal. Then, take the example of straight line fitting.

${\displaystyle f(x_i,\beta )=\alpha +\beta x_i}$

in this case

${\displaystyle M_{ii}={\sigma }_{y,i}^2+{\beta }^2{\sigma }_{x,i}^2}$

showing how the variance at the ith point is determined by the variances of both independent and dependent variables and by the model being used to fit the data. The expression may be generalized by noting that the parameter ${\displaystyle \beta }$ is the slope of the line.

${\displaystyle M_{ii}={\sigma }_{y,i}^2+{\left(\frac{dy}{dx}\right)}_i^2{\sigma }_{x,i}^2}$

An expression of this type is used in fitting pH titration data where a small error on x translates to a large error on y when the slope is large.

Algebraic point of view

As was shown in 1980 by Golub and Van Loan, the TLS problem does not have a solution in general. The following considers the simple case where a unique solution exists without making any particular assumptions.

The computation of the TLS using singular value decomposition (SVD) is described in standard texts. We can solve the equation

${\displaystyle XB\approx Y}$

for B where X is m-by-n and Y is m-by-k.

That is, we seek to find B that minimizes error matrices E and F for X and Y respectively. That is,

${\displaystyle {\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{B,E,F}\Vert [EF]{\Vert }_F,(X+E)B=Y+F}$

where ${\displaystyle [EF]}$ is the augmented matrix with E and F side by side and ${\displaystyle \Vert \cdot {\Vert }_F}$ is the Frobenius norm, the square root of the sum of the squares of all entries in a matrix and so equivalently the square root of the sum of squares of the lengths of the rows or columns of the matrix.

This can be rewritten as

${\displaystyle [(X+E)(Y+F)]\left[\begin{array}{c}B\\ -I_k\end{array}\right]=0.}$

where ${\displaystyle I_k}$ is the ${\displaystyle k\times k}$ identity matrix. The goal is then to find ${\displaystyle [EF]}$ that reduces the rank of ${\displaystyle [XY]}$ by k. Define ${\displaystyle [U][\mathrm{\Sigma }][V{]}^{\ast }}$ to be the singular value decomposition of the augmented matrix ${\displaystyle [XY]}$.

${\displaystyle [XY]=[U_X{U}_Y]\left[\begin{array}{cc}{\mathrm{\Sigma }}_X& 0\\ 0& {\mathrm{\Sigma }}_Y\end{array}\right]{\left[\begin{array}{cc}{V}_{XX}& V_{XY}\\ V_{YX}& V_{YY}\end{array}\right]}^{\ast }=[U_X{U}_Y]\left[\begin{array}{cc}{\mathrm{\Sigma }}_X& 0\\ 0& {\mathrm{\Sigma }}_Y\end{array}\right]\left[\begin{array}{cc}{V}_{XX}^{\ast }& V_{YX}^{\ast }\\ V_{XY}^{\ast }& V_{YY}^{\ast }\end{array}\right]}$

where V is partitioned into blocks corresponding to the shape of X and Y.

Using the Eckart–Young theorem, the approximation minimising the norm of the error is such that matrices ${\displaystyle U}$ and ${\displaystyle V}$ are unchanged, while the smallest ${\displaystyle k}$ singular values are replaced with zeroes. That is, we want

${\displaystyle [(X+E)(Y+F)]=[U_X{U}_Y]\left[\begin{array}{cc}{\mathrm{\Sigma }}_X& 0\\ 0& 0_{k\times k}\end{array}\right]{\left[\begin{array}{cc}{V}_{XX}& V_{XY}\\ V_{YX}& V_{YY}\end{array}\right]}^{\ast }}$

so by linearity,

${\displaystyle [EF]=-[U_X{U}_Y]\left[\begin{array}{cc}{0}_{n\times n}& 0\\ 0& {\mathrm{\Sigma }}_Y\end{array}\right]{\left[\begin{array}{cc}{V}_{XX}& V_{XY}\\ V_{YX}& V_{YY}\end{array}\right]}^{\ast }.}$

We can then remove blocks from the U and Σ matrices, simplifying to

${\displaystyle [EF]=-U_Y{\mathrm{\Sigma }}_Y{\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]}^{\ast }=-[XY]\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]{\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]}^{\ast }.}$

This provides E and F so that

${\displaystyle [(X+E)(Y+F)]\left[\begin{array}{c}{V}_{XY}\\ V_{YY}\end{array}\right]=0.}$

Now if ${\displaystyle V_{YY}}$ is nonsingular, which is not always the case (note that the behavior of TLS when ${\displaystyle V_{YY}}$ is singular is not well understood yet), we can then right multiply both sides by ${\displaystyle -V_{YY}^{-1}}$ to bring the bottom block of the right matrix to the negative identity, giving

${\displaystyle [(X+E)(Y+F)]\left[\begin{array}{c}-V_{XY}{V}_{YY}^{-1}\\ -V_{YY}{V}_{YY}^{-1}\end{array}\right]=[(X+E)(Y+F)]\left[\begin{array}{c}B\\ -I_k\end{array}\right]=0,}$

and so

${\displaystyle B=-V_{XY}{V}_{YY}^{-1}.}$

A naive GNU Octave/MATLAB implementation of this is:

<span id="mw7A"></span><span class="k" id="mw7Q">function</span><span class="w" id="mw7g"> </span>B<span class="w" id="mw7w"> </span><span class="p" id="mw8A">=</span><span class="w" id="mw8Q"> </span><span class="nf" id="mw8g">tls</span><span class="p" id="mw8w">(</span>X, Y<span class="p" id="mw9A">)</span> <span class="w" id="mw9Q"> </span><span class="p" id="mw9g">[</span><span class="n" id="mw9w">m</span><span class="w" id="mw-A"> </span><span class="n" id="mw-Q">n</span><span class="p" id="mw-g">]</span><span class="w" id="mw-w"> </span><span class="p" id="mw_A">=</span><span class="w" id="mw_Q"> </span><span class="nb" id="mw_g">size</span><span class="p" id="mw_w">(</span><span class="n" id="mwAQA">X</span><span class="p" id="mwAQE">);</span><span class="w" id="mwAQI"> </span><span class="c" id="mwAQM">% n is the width of X (X is m by n)</span> <span class="w" id="mwAQQ"> </span><span class="n" id="mwAQU">Z</span><span class="w" id="mwAQY"> </span><span class="p" id="mwAQc">=</span><span class="w" id="mwAQg"> </span><span class="p" id="mwAQk">[</span><span class="n" id="mwAQo">X</span><span class="w" id="mwAQs"> </span><span class="n" id="mwAQw">Y</span><span class="p" id="mwAQ0">];</span><span class="w" id="mwAQ4"> </span><span class="c" id="mwAQ8">% Z is X augmented with Y.</span> <span class="w" id="mwARA"> </span><span class="p" id="mwARE">[</span><span class="n" id="mwARI">U</span><span class="w" id="mwARM"> </span><span class="n" id="mwARQ">S</span><span class="w" id="mwARU"> </span><span class="n" id="mwARY">V</span><span class="p" id="mwARc">]</span><span class="w" id="mwARg"> </span><span class="p" id="mwARk">=</span><span class="w" id="mwARo"> </span><span class="nb" id="mwARs">svd</span><span class="p" id="mwARw">(</span><span class="n" id="mwAR0">Z</span><span class="p" id="mwAR4">,</span><span class="w" id="mwAR8"> </span><span class="mi" id="mwASA">0</span><span class="p" id="mwASE">);</span><span class="w" id="mwASI"> </span><span class="c" id="mwASM">% find the SVD of Z.</span> <span class="w" id="mwASQ"> </span><span class="n" id="mwASU">VXY</span><span class="w" id="mwASY"> </span><span class="p" id="mwASc">=</span><span class="w" id="mwASg"> </span><span class="n" id="mwASk">V</span><span class="p" id="mwASo">(</span><span class="mi" id="mwASs">1</span><span class="p" id="mwASw">:</span><span class="n" id="mwAS0">n</span><span class="p" id="mwAS4">,</span><span class="w" id="mwAS8"> </span><span class="mi" id="mwATA">1</span><span class="o" id="mwATE">+</span><span class="n" id="mwATI">n</span><span class="p" id="mwATM">:</span><span class="k" id="mwATQ">end</span><span class="p" id="mwATU">);</span><span class="w" id="mwATY"> </span><span class="c" id="mwATc">% Take the block of V consisting of the first n rows and the n+1 to last column</span> <span class="w" id="mwATg"> </span><span class="n" id="mwATk">VYY</span><span class="w" id="mwATo"> </span><span class="p" id="mwATs">=</span><span class="w" id="mwATw"> </span><span class="n" id="mwAT0">V</span><span class="p" id="mwAT4">(</span><span class="mi" id="mwAT8">1</span><span class="o" id="mwAUA">+</span><span class="n" id="mwAUE">n</span><span class="p" id="mwAUI">:</span><span class="k" id="mwAUM">end</span><span class="p" id="mwAUQ">,</span><span class="w" id="mwAUU"> </span><span class="mi" id="mwAUY">1</span><span class="o" id="mwAUc">+</span><span class="n" id="mwAUg">n</span><span class="p" id="mwAUk">:</span><span class="k" id="mwAUo">end</span><span class="p" id="mwAUs">);</span><span class="w" id="mwAUw"> </span><span class="c" id="mwAU0">% Take the bottom-right block of V.</span> <span class="w" id="mwAU4"> </span><span class="n" id="mwAU8">B</span><span class="w" id="mwAVA"> </span><span class="p" id="mwAVE">=</span><span class="w" id="mwAVI"> </span><span class="o" id="mwAVM">-</span><span class="n" id="mwAVQ">VXY</span><span class="w" id="mwAVU"> </span><span class="o" id="mwAVY">/</span><span class="w" id="mwAVc"> </span><span class="n" id="mwAVg">VYY</span><span class="p" id="mwAVk">;</span> <span class="k" id="mwAVo">end</span>

The way described above of solving the problem, which requires that the matrix ${\displaystyle V_{YY}}$ is nonsingular, can be slightly extended by the so-called classical TLS algorithm.

Computation

The standard implementation of classical TLS algorithm written by Van Huffel is available through Netlib, see also. All modern implementations based, for example, on solving a sequence of ordinary least squares problems, approximate the matrix ${\displaystyle B}$ (denoted ${\displaystyle X}$ in the literature), as introduced by Van Huffel and Vandewalle. It is worth noting, that this solution is, however, only the TLS solution when ${\displaystyle \mathrm{rank}(W^{(q,e)})=e}$ and ${\displaystyle \mathrm{rank}(V_{12}^{(-e)})=d-e}$. In other cases it represents a ${\displaystyle B_{\text{T-TLS}}}$, a solution of the truncated total least squares problem (T-TLS). The T-TLS has a regularization property that real TLS does not have.