Purpose
To compute a basis for the left and/or right singular subspace of an M-by-N matrix A corresponding to its smallest singular values.Specification
SUBROUTINE MB04XD( JOBU, JOBV, M, N, RANK, THETA, A, LDA, U, LDU,
$ V, LDV, Q, INUL, TOL, RELTOL, DWORK, LDWORK,
$ IWARN, INFO )
C .. Scalar Arguments ..
CHARACTER JOBU, JOBV
INTEGER INFO, IWARN, LDA, LDU, LDV, LDWORK, M, N, RANK
DOUBLE PRECISION RELTOL, THETA, TOL
C .. Array Arguments ..
LOGICAL INUL(*)
DOUBLE PRECISION A(LDA,*), DWORK(*), Q(*), U(LDU,*), V(LDV,*)
Arguments
Mode Parameters
JOBU CHARACTER*1
Specifies whether to compute the left singular subspace
as follows:
= 'N': Do not compute the left singular subspace;
= 'A': Return the (M - RANK) base vectors of the desired
left singular subspace in U;
= 'S': Return the first (min(M,N) - RANK) base vectors
of the desired left singular subspace in U.
JOBV CHARACTER*1
Specifies whether to compute the right singular subspace
as follows:
= 'N': Do not compute the right singular subspace;
= 'A': Return the (N - RANK) base vectors of the desired
right singular subspace in V;
= 'S': Return the first (min(M,N) - RANK) base vectors
of the desired right singular subspace in V.
Input/Output Parameters
M (input) INTEGER
The number of rows in matrix A. M >= 0.
N (input) INTEGER
The number of columns in matrix A. N >= 0.
RANK (input/output) INTEGER
On entry, if RANK < 0, then the rank of matrix A is
computed by the routine as the number of singular values
greater than THETA.
Otherwise, RANK must specify the rank of matrix A.
RANK <= min(M,N).
On exit, if RANK < 0 on entry, then RANK contains the
computed rank of matrix A. That is, the number of singular
values of A greater than THETA.
Otherwise, the user-supplied value of RANK may be changed
by the routine on exit if the RANK-th and the (RANK+1)-th
singular values of A are considered to be equal.
See also the description of parameter TOL below.
THETA (input/output) DOUBLE PRECISION
On entry, if RANK < 0, then THETA must specify an upper
bound on the smallest singular values of A corresponding
to the singular subspace to be computed. THETA >= 0.0.
Otherwise, THETA must specify an initial estimate (t say)
for computing an upper bound on the (min(M,N) - RANK)
smallest singular values of A. If THETA < 0.0, then t is
computed by the routine.
On exit, if RANK >= 0 on entry, then THETA contains the
computed upper bound such that precisely RANK singular
values of A are greater than THETA + TOL.
Otherwise, THETA is unchanged.
A (input) DOUBLE PRECISION array, dimension (LDA,N)
The leading M-by-N part of this array must contain the
matrix A from which the basis of a desired singular
subspace is to be computed.
NOTE that this array is destroyed.
LDA INTEGER
The leading dimension of array A. LDA >= max(1,M).
U (output) DOUBLE PRECISION array, dimension (LDU,*)
If JOBU = 'A', then the leading M-by-M part of this array
contains the (M - RANK) M-dimensional base vectors of the
desired left singular subspace of A corresponding to its
singular values less than or equal to THETA. These vectors
are stored in the i-th column(s) of U for which
INUL(i) = .TRUE., where i = 1,2,...,M.
If JOBU = 'S', then the leading M-by-min(M,N) part of this
array contains the first (min(M,N) - RANK) M-dimensional
base vectors of the desired left singular subspace of A
corresponding to its singular values less than or equal to
THETA. These vectors are stored in the i-th column(s) of U
for which INUL(i) = .TRUE., where i = 1,2,..., min(M,N).
Otherwise, U is not referenced (since JOBU = 'N') and can
be supplied as a dummy array (i.e. set parameter LDU = 1
and declare this array to be U(1,1) in the calling
program).
LDU INTEGER
The leading dimension of array U.
LDU >= max(1,M) if JOBU = 'A' or JOBU = 'S',
LDU >= 1 if JOBU = 'N'.
V (output) DOUBLE PRECISION array, dimension (LDV,*)
If JOBV = 'A', then the leading N-by-N part of this array
contains the (N - RANK) N-dimensional base vectors of the
desired right singular subspace of A corresponding to its
singular values less than or equal to THETA. These vectors
are stored in the i-th column(s) of V for which
INUL(i) = .TRUE., where i = 1,2,...,N.
If JOBV = 'S', then the leading N-by-min(M,N) part of this
array contains the first (min(M,N) - RANK) N-dimensional
base vectors of the desired right singular subspace of A
corresponding to its singular values less than or equal to
THETA. These vectors are stored in the i-th column(s) of V
for which INUL(i) = .TRUE., where i = 1,2,...,MIN( M,N).
Otherwise, V is not referenced (since JOBV = 'N') and can
be supplied as a dummy array (i.e. set parameter LDV = 1
and declare this array to be V(1,1) in the calling
program).
LDV INTEGER
The leading dimension of array V.
LDV >= max(1,N) if JOBV = 'A' or JOBV = 'S',
LDV >= 1 if JOBV = 'N'.
Q (output) DOUBLE PRECISION array, dimension (2*min(M,N)-1)
This array contains the partially diagonalized bidiagonal
matrix J computed from A, at the moment that the desired
singular subspace has been found. Specifically, the
leading p = min(M,N) entries of Q contain the diagonal
elements q(1),q(2),...,q(p) and the entries Q(p+1),
Q(p+2),...,Q(2*p-1) contain the superdiagonal elements
e(1),e(2),...,e(p-1) of J.
INUL (output) LOGICAL array, dimension (max(M,N))
If JOBU <> 'N' or JOBV <> 'N', then the indices of the
elements of this array with value .TRUE. indicate the
columns in U and/or V containing the base vectors of the
desired left and/or right singular subspace of A. They
also equal the indices of the diagonal elements of the
bidiagonal submatrices in the array Q, which correspond
to the computed singular subspaces.
Tolerances
TOL DOUBLE PRECISION
This parameter defines the multiplicity of singular values
by considering all singular values within an interval of
length TOL as coinciding. TOL is used in checking how many
singular values are less than or equal to THETA. Also in
computing an appropriate upper bound THETA by a bisection
method, TOL is used as a stopping criterion defining the
minimum (absolute) subinterval width. TOL is also taken
as an absolute tolerance for negligible elements in the
QR/QL iterations. If the user sets TOL to be less than or
equal to 0, then the tolerance is taken as specified in
SLICOT Library routine MB04YD document.
RELTOL DOUBLE PRECISION
This parameter specifies the minimum relative width of an
interval. When an interval is narrower than TOL, or than
RELTOL times the larger (in magnitude) endpoint, then it
is considered to be sufficiently small and bisection has
converged. If the user sets RELTOL to be less than
BASE * EPS, where BASE is machine radix and EPS is machine
precision (see LAPACK Library routine DLAMCH), then the
tolerance is taken as BASE * EPS.
Workspace
DWORK DOUBLE PRECISION array, dimension (LDWORK)
On exit, if INFO = 0, DWORK(1) returns the optimal value
of LDWORK.
LDWORK INTEGER
The length of the array DWORK.
LDWORK = max(1, LDW + max(2*P + max(M,N), LDY)), where
P = min(M,N);
LDW = max(2*N, N*(N+1)/2), if JOBU <> 'N' and M large
enough than N;
LDW = 0, otherwise;
LDY = 8*P - 5, if JOBU <> 'N' or JOBV <> 'N';
LDY = 6*P - 3, if JOBU = 'N' and JOBV = 'N'.
For optimum performance LDWORK should be larger.
If LDWORK = -1, then a workspace query is assumed;
the routine only calculates the optimal size of the
DWORK array, returns this value as the first entry of
the DWORK array, and no error message related to LDWORK
is issued by XERBLA.
Warning Indicator
IWARN INTEGER
= 0: no warning;
= 1: if the rank of matrix A (as specified by the user)
has been lowered because a singular value of
multiplicity greater than 1 was found.
Error Indicator
INFO INTEGER
= 0: successful exit;
< 0: if INFO = -i, the i-th argument had an illegal
value;
= 1: if the maximum number of QR/QL iteration steps
(30*MIN(M,N)) has been exceeded.
Method
The method used is the Partial Singular Value Decomposition (PSVD)
approach proposed by Van Huffel, Vandewalle and Haegemans, which
is an efficient technique (see [1]) for computing the singular
subspace of a matrix corresponding to its smallest singular
values. It differs from the classical SVD algorithm [3] at three
points, which results in high efficiency. Firstly, the Householder
transformations of the bidiagonalization need only to be applied
on the base vectors of the desired singular subspaces; secondly,
the bidiagonal matrix need only be partially diagonalized; and
thirdly, the convergence rate of the iterative diagonalization can
be improved by an appropriate choice between QL and QR iterations.
(Note, however, that LAPACK Library routine DGESVD, for computing
SVD, also uses either QL and QR iterations.) Depending on the gap,
the desired numerical accuracy and the dimension of the desired
singular subspace, the PSVD can be up to three times faster than
the classical SVD algorithm.
The PSVD algorithm [1-2] for an M-by-N matrix A proceeds as
follows:
Step 1: Bidiagonalization phase
-----------------------
(a) If M is large enough than N, transform A into upper
triangular form R.
(b) Transform A (or R) into bidiagonal form:
|q(1) e(1) 0 ... 0 |
(0) | 0 q(2) e(2) . |
J = | . . |
| . e(N-1)|
| 0 ... q(N) |
if M >= N, or
|q(1) 0 0 ... 0 0 |
(0) |e(1) q(2) 0 . . |
J = | . . . |
| . q(M-1) . |
| 0 ... e(M-1) q(M)|
if M < N, using Householder transformations.
In the second case, transform the matrix to the upper bidiagonal
form by applying Givens rotations.
(c) If U is requested, initialize U with the identity matrix.
If V is requested, initialize V with the identity matrix.
Step 2: Partial diagonalization phase
-----------------------------
If the upper bound THETA is not given, then compute THETA such
that precisely (min(M,N) - RANK) singular values of the bidiagonal
matrix are less than or equal to THETA, using a bisection method
[4]. Diagonalize the given bidiagonal matrix J partially, using
either QR iterations (if the upper left diagonal element of the
considered bidiagonal submatrix is larger than the lower right
diagonal element) or QL iterations, such that J is split into
unreduced bidiagonal submatrices whose singular values are either
all larger than THETA or all less than or equal to THETA.
Accumulate the Givens rotations in U and/or V (if desired).
Step 3: Back transformation phase
-------------------------
(a) Apply the Householder transformations of Step 1(b) onto the
columns of U and/or V associated with the bidiagonal
submatrices with all singular values less than or equal to
THETA (if U and/or V is desired).
(b) If M is large enough than N, and U is desired, then apply the
Householder transformations of Step 1(a) onto each computed
column of U in Step 3(a).
References
[1] Van Huffel, S., Vandewalle, J. and Haegemans, A.
An efficient and reliable algorithm for computing the singular
subspace of a matrix associated with its smallest singular
values.
J. Comput. and Appl. Math., 19, pp. 313-330, 1987.
[2] Van Huffel, S.
Analysis of the total least squares problem and its use in
parameter estimation.
Doctoral dissertation, Dept. of Electr. Eng., Katholieke
Universiteit Leuven, Belgium, June 1987.
[3] Chan, T.F.
An improved algorithm for computing the singular value
decomposition.
ACM TOMS, 8, pp. 72-83, 1982.
[4] Van Huffel, S. and Vandewalle, J.
The partial total least squares algorithm.
J. Comput. and Appl. Math., 21, pp. 333-341, 1988.
Numerical Aspects
Using the PSVD a large reduction in computation time can be gained in total least squares applications (cf [2 - 4]), in the computation of the null space of a matrix and in solving (non)homogeneous linear equations.Further Comments
NoneExample
Program Text
* MB04XD EXAMPLE PROGRAM TEXT
* Copyright (c) 2002-2017 NICONET e.V.
*
* .. Parameters ..
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D0 )
INTEGER NIN, NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER MMAX, NMAX
PARAMETER ( MMAX = 20, NMAX = 20 )
INTEGER LDA, LDU, LDV
PARAMETER ( LDA = MMAX, LDU = MMAX, LDV = NMAX )
INTEGER MAXMN, MNMIN
PARAMETER ( MAXMN = MAX( MMAX, NMAX ),
$ MNMIN = MIN( MMAX, NMAX ) )
INTEGER LENGQ
PARAMETER ( LENGQ = 2*MNMIN-1 )
INTEGER LDWORK
PARAMETER ( LDWORK = MAX( 2*NMAX, NMAX*( NMAX+1 )/2 )
$ + MAX( 2*MNMIN + MAXMN, 8*MNMIN - 5 ) )
* .. Local Scalars ..
DOUBLE PRECISION RELTOL, THETA, THETA1, TOL
INTEGER I, INFO, IWARN, J, K, LOOP, M, MINMN, N, NCOLU,
$ NCOLV, RANK, RANK1
CHARACTER*1 JOBU, JOBV
LOGICAL LJOBUA, LJOBUS, LJOBVA, LJOBVS, WANTU, WANTV
* .. Local Arrays ..
DOUBLE PRECISION A(LDA,NMAX), DWORK(LDWORK), Q(LENGQ),
$ U(LDU,MMAX), V(LDV,NMAX)
LOGICAL INUL(MAXMN)
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL MB04XD
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* .. Executable Statements ..
*
WRITE ( NOUT, FMT = 99999 )
* Skip the heading in the data file and read the data.
READ ( NIN, FMT = '()' )
READ ( NIN, FMT = * ) M, N, RANK, THETA, TOL, RELTOL, JOBU, JOBV
IF ( M.LT.0 .OR. M.GT.MMAX ) THEN
WRITE ( NOUT, FMT = 99983 ) M
ELSE IF ( N.LT.0 .OR. N.GT.NMAX ) THEN
WRITE ( NOUT, FMT = 99982 ) N
ELSE IF ( RANK.GT.MNMIN ) THEN
WRITE ( NOUT, FMT = 99981 ) RANK
ELSE IF ( RANK.LT.0 .AND. THETA.LT.ZERO ) THEN
WRITE ( NOUT, FMT = 99980 ) THETA
ELSE
READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,M )
RANK1 = RANK
THETA1 = THETA
* Compute a basis for the left and right singular subspace of A.
CALL MB04XD( JOBU, JOBV, M, N, RANK, THETA, A, LDA, U, LDU, V,
$ LDV, Q, INUL, TOL, RELTOL, DWORK, LDWORK, IWARN,
$ INFO )
*
IF ( INFO.NE.0 ) THEN
WRITE ( NOUT, FMT = 99998 ) INFO
ELSE
IF ( IWARN.NE.0 ) THEN
WRITE ( NOUT, FMT = 99997 ) IWARN
WRITE ( NOUT, FMT = 99996 ) RANK
ELSE
IF ( RANK1.LT.0 ) WRITE ( NOUT, FMT = 99996 ) RANK
END IF
IF ( THETA1.LT.ZERO ) WRITE ( NOUT, FMT = 99995 ) THETA
LJOBUA = LSAME( JOBU, 'A' )
LJOBUS = LSAME( JOBU, 'S' )
LJOBVA = LSAME( JOBV, 'A' )
LJOBVS = LSAME( JOBV, 'S' )
WANTU = LJOBUA.OR.LJOBUS
WANTV = LJOBVA.OR.LJOBVS
WRITE ( NOUT, FMT = 99994 )
MINMN = MIN( M, N )
LOOP = MINMN - 1
DO 20 I = 1, LOOP
K = I + MINMN
WRITE ( NOUT, FMT = 99993 ) I, I, Q(I), I, I + 1, Q(K)
20 CONTINUE
WRITE ( NOUT, FMT = 99992 ) MINMN, MINMN, Q(MINMN)
IF ( WANTU ) THEN
NCOLU = M
IF ( LJOBUS ) NCOLU = MINMN
WRITE ( NOUT, FMT = 99986 )
DO 40 I = 1, M
WRITE ( NOUT, FMT = 99985 ) ( U(I,J), J = 1,NCOLU )
40 CONTINUE
WRITE ( NOUT, FMT = 99991 ) NCOLU
WRITE ( NOUT, FMT = 99990 )
DO 60 I = 1, NCOLU
WRITE ( NOUT, FMT = 99989 ) I, INUL(I)
60 CONTINUE
END IF
IF ( WANTV ) THEN
NCOLV = N
IF ( LJOBVS ) NCOLV = MINMN
WRITE ( NOUT, FMT = 99984 )
DO 80 I = 1, N
WRITE ( NOUT, FMT = 99985 ) ( V(I,J), J = 1,NCOLV )
80 CONTINUE
WRITE ( NOUT, FMT = 99988 ) NCOLV
WRITE ( NOUT, FMT = 99987 )
DO 100 J = 1, NCOLV
WRITE ( NOUT, FMT = 99989 ) J, INUL(J)
100 CONTINUE
END IF
END IF
END IF
STOP
*
99999 FORMAT (' MB04XD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB04XD = ',I2)
99997 FORMAT (' IWARN on exit from MB04XD = ',I2,/)
99996 FORMAT (' The computed rank of matrix A = ',I3,/)
99995 FORMAT (' The computed value of THETA = ',F7.4,/)
99994 FORMAT (' The elements of the partially diagonalized bidiagonal ',
$ 'matrix are',/)
99993 FORMAT (2(' (',I1,',',I1,') = ',F7.4,2X))
99992 FORMAT (' (',I1,',',I1,') = ',F7.4,/)
99991 FORMAT (/' Left singular subspace corresponds to the i-th column',
$ '(s) of U for which ',/' INUL(i) = .TRUE., i = 1,...,',I1,
$ /)
99990 FORMAT (' i INUL(i)',/)
99989 FORMAT (I3,L8)
99988 FORMAT (/' Right singular subspace corresponds to the j-th colum',
$ 'n(s) of V for which ',/' INUL(j) = .TRUE., j = 1,...,',I1,
$ /)
99987 FORMAT (' j INUL(j)',/)
99986 FORMAT (' Matrix U',/)
99985 FORMAT (20(1X,F8.4))
99984 FORMAT (/' Matrix V',/)
99983 FORMAT (/' M is out of range.',/' M = ',I5)
99982 FORMAT (/' N is out of range.',/' N = ',I5)
99981 FORMAT (/' RANK is out of range.',/' RANK = ',I5)
99980 FORMAT (/' THETA must be at least zero.',/' THETA = ',F8.4)
END
Program Data
MB04XD EXAMPLE PROGRAM DATA 6 4 -1 0.001 0.0 0.0 A A 0.80010 0.39985 0.60005 0.89999 0.29996 0.69990 0.39997 0.82997 0.49994 0.60003 0.20012 0.79011 0.90013 0.20016 0.79995 0.85002 0.39998 0.80006 0.49985 0.99016 0.20002 0.90007 0.70009 1.02994Program Results
MB04XD EXAMPLE PROGRAM RESULTS The computed rank of matrix A = 3 The elements of the partially diagonalized bidiagonal matrix are (1,1) = 3.2280 (1,2) = -0.0287 (2,2) = 0.8714 (2,3) = 0.0168 (3,3) = 0.3698 (3,4) = 0.0000 (4,4) = 0.0001 Matrix U 0.8933 0.4328 -0.1209 0.2499 -0.5812 0.4913 -0.4493 0.8555 -0.2572 0.1617 -0.4608 -0.7379 -0.0079 0.2841 0.9588 -0.5352 0.1892 0.0525 0.0000 0.0000 0.0003 -0.1741 0.3389 -0.3397 0.0000 0.0000 0.0000 0.6482 0.5428 0.1284 0.0000 0.0000 0.0000 -0.4176 -0.0674 0.2819 Left singular subspace corresponds to the i-th column(s) of U for which INUL(i) = .TRUE., i = 1,...,6 i INUL(i) 1 F 2 F 3 F 4 T 5 T 6 T Matrix V -0.3967 -0.7096 0.4612 -0.3555 0.9150 -0.2557 0.2414 -0.5687 -0.0728 0.6526 0.5215 -0.2128 0.0000 0.0720 0.6761 0.7106 Right singular subspace corresponds to the j-th column(s) of V for which INUL(j) = .TRUE., j = 1,...,4 j INUL(j) 1 F 2 F 3 F 4 T