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.