SUBROUTINE : ESTIM 1 PURPOSE: The subroutine ESTIM estimates an upper bound THETA using a bisec- tion method such that the bidiagonal J defined by Q and E: !Q(1) E(2) 0 ... 0 ! ! 0 Q(2) E(3) . ! J = ! . . ! ! . E(N)! ! 0 ... Q(N)! has precisely L singular values <= THETA + a given tolerance TOL1. 2 SPECIFICATION: SUBROUTINE ESTIM(Q, E, N, L, THETA, TOL1, TOL2, IWARN) INTEGER N, L, IWARN DOUBLE PRECISION Q(N), E(N), THETA, TOL1, TOL2 3 ARGUMENT LIST: 3.1 ARGUMENTS IN Q - DOUBLE PRECISION array of DIMENSION (N). Contains the diagonal elements of the bidiagonal J. E - DOUBLE PRECISION array of DIMENSION (N). E(i), i=2,...,N, contain the superdiagonal elements of the bidiagonal J, E(1) = 0.0D0. N - INTEGER. N is the dimension of the bidiagonal J defined by Q and E. L - INTEGER. L is the number of singular values of J which must be smaller than or equal to the upper bound computed by ESTIM. NOTE that L may be overwritten. THETA - DOUBLE PRECISION THETA is an initial estimate of the upper bound which is to be computed by ESTIM. If THETA < 0, then the initialization is computed by ESTIM. NOTE that THETA is overwritten. 3.2 ARGUMENTS OUT L - INTEGER. On return, L can be increased if the L-th smallest singular value of J has multiplicity > 1. In this case, L is increased by the number of singular values of J which are larger than its L-th smallest one and approach the L-th smallest singular value of J within a distance < TOL1. If L has been increased, a warning is given. THETA - DOUBLE PRECISION On return, THETA contains the upper bound computed by ESTIM such that the bidiagonal J has precisely L singular values <= THETA + TOL1. 3.4 TOLERANCES TOL1 - DOUBLE PRECISION. This parameter defines the multiplicity of singular values by considering all singular values within an interval of length TOL1 as coinciding. TOL1 is used in checking how many singu- lar values are <= THETA. Also in computing an appropriate upper bound THETA by a bisection method, TOL1 is used as stop criterion defining the minimal subinterval length. TOL2 - DOUBLE PRECISION. This parameter specifies that matrix elements Q(i) and/or E(i), which are <= TOL2 in absolute value, are considered to be zero. 3.6 ERROR INDICATORS IWARN - INTEGER. On return, IWARN = 0 unless the value of L has been increased. 4 ERROR INDICATORS and WARNINGS: Warning given by the routine: IWARN = 0: no warnings. IWARN = 1: the L-th smallest singular value of J coincides with the (L+1)-th smallest one. Therefore, the number L of singular values which must be <= the upper bound computed by ESTIM, has been increased. 5 EXTERNAL SUBROUTINES and FUNCTIONS: DAMIN, NSINGV. 6 METHOD DESCRIPTION: Let S(i), i=1,...,N, be the N singular values of the bidiagonal J in decreasing order of magnitude. ESTIM then computes an upper bound T such that S(N-L) > T >= S(N-L+1). This is done as follows (see [2]): First, if the initial estimate of THETA is not specified by the user then the subroutine initializes THETA with an estimate which is close to the requested value of THETA if S(N-L) >> S(N-L+1). Second, a bisection method (see [1, Sec.8.5]) is used which generates a sequence of shrinking intervals [Y,Z] such that (number of S(i) <= Y) < L < (number of S(i) <= Z). This bisection method is applied to an associated 2N by 2N symmetric tridiagonal matrix T" whose eigenvalues are the singular values of J and their negatives (see METHOD DESCRIPTION of the routine NSINGV). One of the starting values for the bisection method is the initial value of THETA. If this value is an upper bound, then the initial lower bound is set to zero, else the initial upper bound is computed from the Gershgorin Circle Theorem [1, Theorem 7.2-1], applied to T". The computation of the "number of S(i) <= Y (or Z)" is accomplished by the routine NSINGV and is based on applying Sylvester's Law of Inertia or equivalently Sturm sequences [1, Sec.8.5] to the associ- ated matrix T" (see METHOD DESCRIPTION of NSINGV). If at a certain moment Z-Y < TOL1, then at least 2 singular values of J lie in the interval [Y,Z] within a distance < TOL1 from each other. In this case, S(N-L) and S(N-L+1) are assumed to coincide. Then, the upper bound T is set to the value of Z , the value of L is increased and a warning is given to the user. 7 REFERENCES: [1] G.H. Golub and C.F. Van Loan, Matrix Computations. The Johns Hopkins University Press, Baltimore, Maryland (1983). [2] S. Van Huffel and J. Vandewalle, The Partial Total Least Squares Algorithm. J. Comput. and Applied Math., 21 (1988), to appear. *********************************************************************** 1988, February 15. 