E230 Constrained and Unconstrained Linear Least Squares Fitting

The TL package finds the least squares solution to a set of unweighted linear equations, possibly subject to a set of equality constraints. The solution is found by Householder triangularisation (see Ref. 1 for details) with parameter elimination if constraints are present. This write-up ends with a few words on generalised least squares fitting (unequal weighting) which is a simple application of the TL package.

All matrices are assumed to be stored row-wise and without gaps, contrary to the Fortran convention, i.e., if the Fortran statement DIMENSION A(NJ,NI) reserves memory for the matrix A the element $A$_ij is found in word A(J,I).

Structure:

SUBROUTINE subprograms
User Entry Names: TLSC, TLS, TLERR, TLRES
Internal Entry Names: TLSMSQ, TLSWOP, TLUK, TLSTEP, TLPIV

Usage:

General Description
Consider the set of M linear equations

$\sum$_j=1^NA_ijx_j= b_i&sp;(i=1,2,...,M with N≤M)

to be solved such that the Euclidian norm $||Ax\; -\; b||$₂= S² is minimised. Instead of determining x from the Normal Equation $x\; =\; (A\text{'}A)-1A\text{'}b$

it is found by applying successive Householder transformations ( Q) which reduce A to upper triangular form without changing the norm of the columns of A or the vector b. This is beneficial from the point of view of stability and flexibility of application. Writing

$QA\; =\; R\; =\; [R$₁O]

$N\; rows$

$M-N\; rows$

₁y₂]

$N\; rows$

$M-N\; rows$

$$

we have that $||Rx\; -\; y||$₂= ||Ax - b||₂

and the vector x is obtained by backward substitution in $R$₁x = y₁ . As a byproduct, the sum of squares of residuals is directly calculated as $S2=||y$₂||₂ .

Now consider A and b to be composed of $M$₁ constraints to be satisfied exactly, followed by $M-M$₁ equations to be minimised. Writing

$A\; =\; [A$₁A₂]

$M$1 rows

$M-M$1 rows

₁b₂]

$M$1 rows

$M-M$1 rows

$$

then $||A$₂x-b₂||₂= S² has to be minimized subject to $A$₁x-b₁= 0 .

This problem is solved by eliminating $M$₁ parameters and then evaluating the reduced set of parameters (see Ref. 2 for details).

An attractive feature of the unitary Householder transformations is that when each parameter is eliminated ("solved for") column pivoting allows the selection of that parameter which gives the maximum reduction in the current value of $S2$ . Thus it is possible to terminate the calculation whenever $S2$ or its current reduction become acceptably small. This can be exploited when iterating. If there is more than one RHS vector, then x and b become $N\; xL$ and $M\; xL$ matrices with the pivoting strategy applied to the first column of b.

The triangular form of $R$₁ allows the error matrix, E, of the fitted parameters to be derived directly from $R$₁

itself without inverting. The equation is

$E=R$₁^-1(R₁^-1)'.

Moreover, the vector of fitted residuals is most easily computed by applying the inverse Householder transformation to $y$₂ , i.e.

$Ax-b=Q-1[O\; y$₂].

Note that these residuals do not have to be calculated to find the fitted $S2$ which is output from the fitting routines.

In all routines described below, the dimensionality of the problem is transmitted via the common block

    COMMON /TLSDIM/ M1,M,N,L,IER

    CALL TLSC(A,B,AUX,IPIV,EPS,X)

A

( REAL) The combined constraint / derivative matrix of dimension $M\; xN$ , the upper M1 rows being the constraints.

B

( REAL) The combined constraint / measurement matrix of dimension $M\; xL$ , the upper M1 rows being the constraints.

X

( REAL) The matrix of dimension $N\; xL$

returning the L least squares solutions.

AUX

( REAL) Working array of length $N+max(N,L)$ . On output AUX(J),(J=1,L) contain the minimised sum of squares.

IPIV

( INTEGER) Working array of length N which holds the exchange information (column pivoting is employed if necessary).

EPS

( REAL) Parameter specifying a pivoting criterium. There is no exchange of columns I and 1 unless $EPS*PIVOT(I)\; >\; PIVOT(1)$ . Typically $EPS\; 0.1$ .

Unconstrained Least Squares Fitting

    CALL TLS(A,B,AUX,IPIV,EPS,X)

A

( REAL) $M\; xN$ derivative matrix.

B

( REAL) $M\; xL$ matrix of measurements.

X

( REAL) $N\; xL$ parameter solution matrix.

AUX

( REAL) Working array as for TLSC.

IPIV

( INTEGER) Working array as for TLSC.

EPS

( REAL) Input parameter used for prematurely terminating the calculation:
$>\; 0:$ Termination when r.m.s. residual ,
Termination when the reduction in the residual ,
$=\; 0:$ Unconditionally solve for all $X$_j.

    CALL TLERR(A,E,AUX,IPIV)

    CALL TLRES(A,B,AUX)

Notes:

1. The pivoting and exit criteria of TLS are calculated using the first vector of measurements; therefore it is wise to have $EPS=0$ if $L>1$ .
2. TLERR and/or TLRES may be called in any order after TLS or TLSC.
3. TLS or TLSC may be used for solving simultaneous linear equations by setting $M=N$ or $M1=N$ .
4. Useful examples in the application of these routines can be found in the HYDRA Geometry / Kinematics processors.

Generalized Least Squares Fitting
The problem is to minimise $(Ax-b)\text{'}G(Ax-b)$ where G, the weight matrix, is the inverse of the error matrix of the measurement vector b. Once again Householder triangularisation offers an attractive alternative to the Normal Equation solution $x=(A\text{'}GA)-1A\text{'}Gb$ . The first step is to perform the Choleski decomposition of G, which is positive semi-definite (see TR (F112)), such that $G=U\text{'}U$ , U being upper triangular. The problem is then reduced to minimising $||A1x-b1||$₂ , where $A1=UA$ and $b1=Ub$ , which is just the unweighted case previously described. This has the feature that if A has already been triangularised then the product UA remains triangular and only back substitution is necessary to find the weighted least squares solution.

References:

1. G. Golub, Numerical methods for solving linear least squares problems, Numer. Math. 7 (1965) 206--216.
2. Å. Björck and G. Golub, Iterative refinement of linear least square solutions by Householder transformation, BIT 7 (1967) 322--337.

E250