*********************************************************************** Fortran 77 Program of the EXTENDED CLASSICAL TOTAL LEAST SQUARES ALGORITHM. ----------------------------------------------------- Sabine VAN HUFFEL ESAT Laboratory, KU Leuven. Kardinaal Mercierlaan 94, 3030 Heverlee, Belgium *********************************************************************** SUBROUTINE : DTLS 1 PURPOSE: The subroutine DTLS solves the Total Least Squares (TLS) problem. The TLS problem assumes an overdetermined set of linear equations AX = B, where both the data matrix A as well as the observation matrix B are inaccurate. If the perturbations D on the data [A;B] have zero mean and their covariance matrix E(D'D), with E the expected value operator, equals the identity matrix up to an unknown scaling factor (e.g. when all errors are independent and equally sized), then the routine computes a strongly consistent estimate of the true solution of the corresponding unperturbed set. The routine also solves determined and underdetermined sets of equations by computing the minimum norm solution. It is assumed that all preprocessing measures (scaling, coordinate transformations, whitening, ... ) of the data have been performed in advance. 2 SPECIFICATION: SUBROUTINE DTLS(C, LDC, M, N, L, S, X, LDX, WRK, RANK, TOL1, TOL2, COMPRT, IERR, IWARN) INTEGER LDC, M, N, L, RANK, LDX, IERR, IWARN DOUBLE PRECISION TOL1, TOL2 DOUBLE PRECISION C(LDC,N+L), S(N+L), WRK(N+L+M), X(LDX,L) CHARACTER COMPRT*1 3 ARGUMENT LIST: 3.1 ARGUMENTS IN C - DOUBLE PRECISION array of DIMENSION (LDC,N+L). The leading M x (N + L) part of this array contains A and B as follows: The first N columns contain the data matrix A, and the last L columns the observation matrix B (right-hand sides). NOTE that this array is overwritten. LDC - INTEGER. The declared leading dimension of the array C. LDC >= max(M,(N+L)). M - INTEGER. The number of rows in the data matrix A and the observation matrix B. N - INTEGER. The number of columns in the data matrix A. L - INTEGER. The number of columns in the observation matrix B. RANK - INTEGER. The rank of the TLS approximation [A+DA;B+DB]. On entry there are two possibilities (depending on COMPRT): i) RANK is specified by the user. RANK <= min(M,N). NOTE that RANK may be overwritten, if C has multiple singular values. ii) RANK is computed by the routine. 3.2 ARGUMENTS OUT C - DOUBLE PRECISION array of DIMENSION (LDC,N+L). The first NL = N + L rows contain the NL right singular vectors of C = [A;B] S - DOUBLE PRECISION array of DIMENSION (N+L). The first min(M,N+L) entries of S contain the singular values of matrix C, arranged in descending order of magnitude. X - DOUBLE PRECISION array of DIMENSION (LDX,L). The leading N x L part of this array contains the solutions to the TLS problems specified by A and B. LDX - INTEGER. The declared leading dimension of the array X . LDX >= N. RANK - INTEGER. If not specified by the user (see COMPRT) then RANK is computed by the routine. If specified by the user (see COMPRT) then the specified RANK is changed by the routine if the RANK-th and the (RANK+1)-th singular value of C are considered to be equal or if the upper triangular matrix F (see 6 METHOD DESCRIP- TION) is singular. 3.3 WORK SPACE WRK - DOUBLE PRECISION array of DIMENSION (N+L+M). 3.4 TOLERANCES TOL1 - DOUBLE PRECISION. In computing the rank of the approximation [A+DA;B+DB] or in checking the multiplicity of singular values, TOL1 spe- cifies that S(i) and S(j) are considered to be equal if sqrt(S(i)**2 - S(j)**2) is less than or equal to TOL1 (with S(j) zero in case of rank computation). On entry there are two possibilities (depending on COMPRT): i) TOL1 is specified by the user. ii) TOL1 is computed by the routine. In this case the user must enter the value SDEV, i.e. the estimated standard deviation of the error on each element of the matrix C, as input value of TOL1. TOL1 >= machine precision. If not, DTLS corrects its value. TOL2 - DOUBLE PRECISION. In checking the singularity of the upper triangular matrix F (see description), TOL2 specifies that F is considered to be singular if the absolute value of one of its diagonal elements is smaller than or equal to TOL2. TOL2 >= machine precision. If not, DTLS corrects its value. 3.5 MODE PARAMETER COMPRT - CHARACTER*1. Specifies whether RANK and/or TOL1 is to be computed. COMPRT = 'R' or 'r': RANK ) 'T' or 't': TOL1 ) is to be computed. 'N' or 'n': neither ) any other character: both RANK and TOL1 will be computed. 3.6 ERROR INDICATORS IERR - INTEGER. On return, IERR contains 0 unless the subroutine has failed. IWARN - INTEGER. On return, IWARN contains 0 unless RANK has been lowered by the routine. 4 ERROR INDICATORS and WARNINGS: Errors detected by the routine: IERR = 0: successful completion. 1: number M of rows of array C = [A;B] is smaller than 1. 2: number N of columns of matrix A is smaller than 1. 3: number L of columns of matrix B is smaller than 1. 4: leading dimension LDC of array C is smaller than max(M,N+L). 5: leading dimension LDX of array X is smaller than N. 6: rank of the TLS approximation [A+DA;B+DB] is larger than min(M,N). 1xxx: LINPACK routine DSVDC was unable to compute all of the singular values of C. xxx is the error code from LINPACK routine DSVDC. Warnings given by the routine: IWARN = 0: no warnings. 1: the rank of matrix C has been lowered because a singular value of multiplicity > 1 has been found. 2: the rank of matrix C has been lowered because a singular upper triangular matrix F has been found. 5 EXTERNAL SUBROUTINES and FUNCTIONS: DAXPY, DDOT, DCOPY from BLAS [4]; DSVDC, DQRDC from LINPACK [5]; HOUSH, TR2 . 6 METHOD DESCRIPTION: Let [A;B] denote the matrix formed by adjoining the columns of B to the columns of A on the right. Total Least Squares (TLS) definition: Given matrices A and B, find a matrix X satisfying (A + DA) X = B + DB, where A and DA are M x N matrices, B and DB are M x L matrices, X is an N x L matrix. The solution X must be such that the Frobenius norm of [DA;DB] is a minimum and each column of B + DB is in the range of A + DA. Whenever the solution is not unique, DTLS singles out the mini- mum norm solution X. Define matrix C = [A;B] and S(i) as its i-th singular value, i = 1,...,(N+L). Let NL = N + L. The Classical TLS algorithm proceeds as follows (see [3]): Step 1: 1.a): If M >= 5*NL/3, transform C into upper triangular form R by Householder transformations. 1.b): Compute the singular value decomposition (SVD) USV' of C (or R) Step 2: If not fixed by the user, compute the rank R0 of the data [A;B] as follows: S(1) >= ... >= S(R0) > TOL1 >= ... >= S(NL). Using [2] TOL1 can be computed from the standard deviation SDEV of the errors on [A;B]: TOL1 = sqrt(2 * max(M,NL)) * SDEV. The rank R of the approximation [A+DA;B+DB] is then equal to the minimum of N and R0. Step 3: Let V2 be the matrix of the columns of V corresponding to the (NL - R) smallest singular values of C, i.e. the last (NL - R) columns of V. Compute with Householder transformations the orthogonal matrix Q such that: !VH Y! V2 x Q = ! ! !0 F! with VH being N x (N - R), Y is N x L and F is L x L upper triangular. If F is singular, then lower the rank R with the multiplicity of S(R) and repeat this step. Step 4: If F is not singular then the solutions X are obtained by solving the following equations by forward elimination: X F = -Y End Notes : - The TLS solution is unique if R = N and F is not singular and S(N) > S(N+1). - If F is singular, Golub and Van Loan (see [1]) claim that there is no TLS solution. It is however proven in [3] that for those cases, the TLS solution still makes sense and can be determined by further lowering the rank (as explained above). The obtained TLS solution remains optimal with respect to the TLS criteria in case of one ob- servation vector (L=1), as well as in case of more than one obser- vation vector (L>1) (see [3] for proofs). - DTLS is an extension [3] of the classical TLS method described by Golub and Van Loan [1]. 7 REFERENCES : [1] G.H. Golub and C.F. Van Loan, An Analysis of the Total Least Squares Problem. SIAM J. Numer. Anal., 17 (1980), 883-893. [2] J. Staar, J. Vandewalle and M. Wemans, Realization of truncated impulse response sequences with prescribed uncertainty. Proc. IFAC World Congress Kyoto, Vol.I (1981), 7-12. [3] S. Van Huffel, Analysis of the total least squares problem and its use in parameter estimation. Doctoral dissertation, Dept. of Electr. Eng., K.U.Leuven, June 1987. [4] C.L. Lawson, R.J. Hanson, F.T. Krogh and O.R. Kincaid, Basic Lin- ear Algebra Subprograms for FORTRAN Usage. ACM Trans. Math. Soft- ware, 5 (1979), 308-323. [5] J.J. Dongarra, J.R. Bunch, C.B. Moler and G.W. Stewart, LINPACK User's Guide. SIAM, Philadelphia (1979). 9 EXAMPLE: 9.1 PROGRAM TEXT PARAMETER (IN=5, IOUT=6) PARAMETER (LC=14, LM=10, LN=10, LL=4) INTEGER LDC, M, N, L, LDX, IERR, NL, RANK, IWARN DOUBLE PRECISION C(LC,LN+LL), X(LN,LL), WRK(LM+LN+LL), S(LN+LL), * SDEV, TOL2 INTEGER I, J, K CHARACTER MODE*1 LDC = LC LDX = LN READ(IN, 11) M, N, L WRITE(IOUT, 12) M, N, L NL = N + L WRITE(IOUT, 8) DO 1 I = 1, M READ(IN, 4) (C(I,J), J=1,NL) 1 CONTINUE WRITE(IOUT, 5) DO 2 I = 1, M WRITE (IOUT, 6) (C(I,J), J=1,NL) 2 CONTINUE MODE = 'X' SDEV = 1.0D-04 TOL2 = 1.0D-04 CALL DTLS(C, LDC, M, N, L, S, X, LDX, WRK, RANK, SDEV, TOL2, * MODE, IERR, IWARN) WRITE(IOUT, 7) IERR, IWARN WRITE(IOUT,14) RANK WRITE(IOUT, 9) (S(I), I=1,NL) WRITE(IOUT,13) DO 3 I = 1, L WRITE (IOUT, 10) I, (X(J,I), J=1,N) 3 CONTINUE STOP 4 FORMAT(4D15.0) 5 FORMAT(/' ',5X,'C(*,I)',8X,'C(*,I+1)',7X,'C(*,I+2)',7X, * 'C(*,I+3)') 6 FORMAT(' ',4D15.6) 7 FORMAT(/' IERR = ',I3,' IWARN = ',I3) 8 FORMAT(/' TOTAL LEAST SQUARES TEST PROGRAM'/1X,32('-')/) 9 FORMAT(/' SINGULAR VALUES S(I) OF C = '/(' ',4D15.6)) 10 FORMAT(/' X(*,'I2,') = '/(' ',4D15.6)) 11 FORMAT(3I3) 12 FORMAT(/' M = ',I3,3X,'N = ',I3,3X,'L = ',I3) 13 FORMAT(/' TLS SOLUTION :'/1X,12('*')) 14 FORMAT(/' COMPUTED RANK =',I3) END 9.2 PROGRAM DATA 6 3 1 0.80010002D+00 0.39985167D+00 0.60005390D+00 0.89999446D+00 0.29996484D+00 0.69990689D+00 0.39997269D+00 0.82997570D+00 0.49994235D+00 0.60003167D+00 0.20012361D+00 0.79011189D+00 0.90013643D+00 0.20016919D+00 0.79995025D+00 0.85002662D+00 0.39998539D+00 0.80006338D+00 0.49985474D+00 0.99016399D+00 0.20002274D+00 0.90007114D+00 0.70009777D+00 0.10299439D+01 9.3 PROGRAM RESULTS M = 6 N = 3 L = 1 TOTAL LEAST SQUARES TEST PROGRAM -------------------------------- C(*,I) C(*,I+1) C(*,I+2) C(*,I+3) 0.800100D+00 0.399852D+00 0.600054D+00 0.899994D+00 0.299965D+00 0.699907D+00 0.399973D+00 0.829976D+00 0.499942D+00 0.600032D+00 0.200124D+00 0.790112D+00 0.900136D+00 0.200169D+00 0.799950D+00 0.850027D+00 0.399985D+00 0.800063D+00 0.499855D+00 0.990164D+00 0.200023D+00 0.900071D+00 0.700098D+00 0.102994D+01 IERR = 0 IWARN = 0 COMPUTED RANK = 3 SINGULAR VALUES S(I) OF C = 0.322815D+01 0.871560D+00 0.369726D+00 0.128626D-03 TLS SOLUTION : ************ X(*, 1) = 0.500254D+00 0.800251D+00 0.299492D+00 *********************************************************************** 1988, February 15.