E230: Constrained and Unconstrained Linear Least Squares Fitting

Author(s): W. Hart, W. Matt Library: KERNLIB
Submitter: Submitted: 01.01.1975
Language: Fortran Revised: 04.02.1986

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 tex2html_wrap_inline281 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

displaymath285

to be solved such that the Euclidian norm tex2html_wrap_inline287 is minimised. Instead of determining x from the Normal Equation tex2html_wrap_inline289 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

displaymath291

we have that tex2html_wrap_inline293 and the vector x is obtained by backward substitution in tex2html_wrap_inline295 . As a byproduct, the sum of squares of residuals is directly calculated as tex2html_wrap_inline297 .

Now consider A and b to be composed of tex2html_wrap_inline299 constraints to be satisfied exactly, followed by tex2html_wrap_inline301 equations to be minimised. Writing

displaymath303

then tex2html_wrap_inline305 has to be minimized subject to tex2html_wrap_inline307 .

This problem is solved by eliminating tex2html_wrap_inline309 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 tex2html_wrap_inline311 . Thus it is possible to terminate the calculation whenever tex2html_wrap_inline313 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 tex2html_wrap_inline315 and tex2html_wrap_inline317 matrices with the pivoting strategy applied to the first column of b.

The triangular form of tex2html_wrap_inline319 allows the error matrix, E, of the fitted parameters to be derived directly from tex2html_wrap_inline321 itself without inverting. The equation is

displaymath323

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

displaymath327

Note that these residuals do not have to be calculated to find the fitted tex2html_wrap_inline329 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
The parameter IER returns the number of parameters solved for, or else -1001 if either tex2html_wrap_inline331 , tex2html_wrap_inline333 or A has rank less than N.
Constrained Least Squares Fitting 
    CALL TLSC(A,B,AUX,IPIV,EPS,X)
A
(REAL) The combined constraint / derivative matrix of dimension tex2html_wrap_inline335 , the upper M1 rows being the constraints.
B
(REAL) The combined constraint / measurement matrix of dimension tex2html_wrap_inline337 , the upper M1 rows being the constraints.
X
(REAL) The matrix of dimension tex2html_wrap_inline339 returning the L least squares solutions.
AUX
(REAL) Working array of length tex2html_wrap_inline341 . 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 tex2html_wrap_inline347 . Typically tex2html_wrap_inline349 .
Subroutines called: TLSMSQ, TLSWOP, TLUK, TLSTEP.
When constraint equations are present, the full pivoting strategy cannot be adopted and so all parameters are solved for, i.e., IER returns the value N or -1001. Under these circumstances EPS is used to reduce the amount of pivoting to those cases where it is felt to be absolutely necessary.

Unconstrained Least Squares Fitting 

    CALL TLS(A,B,AUX,IPIV,EPS,X)
A
(REAL) tex2html_wrap_inline351 derivative matrix.
B
(REAL) tex2html_wrap_inline353 matrix of measurements.
X
(REAL) tex2html_wrap_inline355 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:
tex2html_wrap_inline357 Termination when r.m.s. residual tex2html_wrap_inline359 ,
tex2html_wrap_inline361 Termination when the reduction in the residual tex2html_wrap_inline363 ,
tex2html_wrap_inline365 Unconditionally solve for all tex2html_wrap_inline367
Subroutines called: TLSMSQ, TLSWOP, TLUK, TLSTEP, TLPIV.
As previously indicated, full pivoting is possible without constraints, hence the allowance for premature exit.
Fitted Error Matrix 
    CALL TLERR(A,E,AUX,IPIV)
The parameter and subroutine arguments defined previously in COMMON /TLSDIM/ require the output values from a call to TLS or TLSC. E is an tex2html_wrap_inline369 matrix which, upon return, will contain the unnormalised covariance matrix of the fitted parameters, tex2html_wrap_inline371 . A may be overwritten by E and the routine may be called independently from TLS/TLSC by setting IER to zero.
Subroutines called: TLUK, TLSTEP.
Fitted Residuals 
    CALL TLRES(A,B,AUX)
All the arguments and common variables require the output values from a call to TLS or TLSC. Upon return, B will give the matrix of residuals, i.e., for each set of least squares equations the column vector tex2html_wrap_inline373 .
Subroutine called: TLSTEP.

Notes:

  1. The pivoting and exit criteria of TLS are calculated using the first vector of measurements; therefore it is wise to have tex2html_wrap_inline375 if tex2html_wrap_inline377 .
  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 tex2html_wrap_inline379 or tex2html_wrap_inline381 .
  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 tex2html_wrap_inline383 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 tex2html_wrap_inline385 . The first step is to perform the Choleski decomposition of G, which is positive semi-definite (see TR (F112)), such that tex2html_wrap_inline387 , U being upper triangular. The problem is then reduced to minimising tex2html_wrap_inline389 , where tex2html_wrap_inline391 and tex2html_wrap_inline393 , 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.
tex2html_wrap_inline395

Michel Goossens Wed Jun 5 03:21:20 METDST 1996