*********************************************************************** Fortran 77 Program of the PARTIAL SINGULAR VALUE DECOMPOSITION ALGORITHM. --------------------------------------------------- Sabine VAN HUFFEL ESAT Laboratory, KU Leuven. Kardinaal Mercierlaan 94, 3030 Heverlee, Belgium. *********************************************************************** SUBROUTINE : PSVD 1 PURPOSE: The subroutine PSVD computes in an efficient and reliable way a basis for the left and/or right singular subspace of a matrix corresponding to its smallest singular values. The dimension of the desired sub- space may be given or may depend on a given upper bound for those smallest singular values. 2 SPECIFICATION: SUBROUTINE PSVD(A, LDA, M, N, RANK, THETA, U, LDU, V, LDV, Q, INUL, WRK, TOL1, TOL2, MODE, IERR, IWARN) INTEGER LDA, M, N, RANK, LDU, LDV, MODE, IERR, IWARN DOUBLE PRECISION THETA, TOL1, TOL2 DOUBLE PRECISION A(LDA,N), U(LDU,M), V(LDV,N),Q(min(M,N)+min(M+1,N)), WRK(*) LOGICAL INUL(max(M,N)) 3 ARGUMENT LIST: 3.1 ARGUMENTS IN A - DOUBLE PRECISION array of DIMENSION (LDA,N). The leading M x N part of this array contains the matrix A for which the basis of a desired singular subspace must be computed. NOTE that A is destroyed by PSVD. LDA - INTEGER. LDA is the leading dimension of the array A (LDA >= M). M - INTEGER. M is the number of rows of the matrix A. N - INTEGER. N is the number of columns of the matrix A. RANK - INTEGER. Specifies on entry whether the rank is given or is to be computed. i) RANK < 0: the rank is to be computed (see THETA). ii) RANK >= 0: specifies the given rank. NOTE that RANK may be overwritten if A has multiple singular values. THETA - DOUBLE PRECISION. On entry, there are two possibilities (depending on RANK): i) RANK < 0: THETA specifies an upper bound on the smallest singular values of A corresponding to the singular sub- space to be computed. THETA >= 0. THETA allows to compute the rank of A. ii) RANK >= 0: THETA is an initial estimate for computing an upper bound on the min(M,N) - RANK smallest singular values of A. If not available, assign a negative value (< 0) to THETA. NOTE that THETA is overwritten. 3.2 ARGUMENTS OUT RANK - INTEGER. If not specified by the user, then RANK is computed by the routine. If specified by the user, then the specified RANK is changed by the routine if the RANK-th and the (RANK+1)-th singular value of A are considered to be equal. THETA - DOUBLE PRECISION. If RANK >= 0, then THETA specifies the computed upper bound such that precisely RANK singular values of A are greater than THETA + TOL1. U - DOUBLE PRECISION array of DIMENSION (LDU,M). The leading M x S part of this array (where S = min(M,N) or S = M, depending on MODE) contains the S - RANK M-dimensional base vectors of the desired left singular subspace of A cor- responding to its singular values <= THETA. These vectors are stored in those columns of U whose index i corresponds with INUL(i) = .TRUE. If 10 <= MODE < 20, then S = M. MODE >= 20, then S = min(M,N). 10 > MODE , then U is not referenced. U may not be identified with A in the subroutine call. LDU - INTEGER. LDU is the leading dimension of the array U (LDU >= M). V - DOUBLE PRECISION array of DIMENSION (LDV,N). The leading N x S part of this array (where S = min(M,N) or S = N, depending on MODE) contains the S - RANK N-dimensional base vectors of the desired right singular subspace of A cor- responding to its singular values <= THETA. These vectors are stored in those columns of V whose index i corresponds with INUL(i) = .TRUE. If MODE mod 10 = 1, then S = N, else S = min(M,N). If MODE mod 10 = 0, then V is not referenced. V may not be identified with A in the subroutine call. LDV - INTEGER. LDV is the leading dimension of the array V (LDV >= N). Q - DOUBLE PRECISION array of DIMENSION (min(M,N)+min(M+1,N)). Returns the partially diagonalized bidiagonal computed from A, at the moment that the desired singular subspace has been found. The first p = min(M,N) entries of Q contain the diagonal elements q(1),...,q(p) and the entries Q(p+2),...,Q(p+s) (with s = min(M+1,N)) contain the superdiagonal elements e(2),...,e(s). Q(p+1) = 0. INUL - LOGICAL array of DIMENSION (max(M,N)). The indices of the elements of INUL 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, if computed. They also equal the indices of the diago- gonal entries of the subbidiagonals in Q, corresponding to the computed singular subspaces. 3.3 WORK SPACE WRK - DOUBLE PRECISION array of DIMENSION (t). If M < 5*N/3 or if no basis of a left singular subspace is requested (i.e. if MODE < 10), then t equals M + N, else t equals M + N + N*(N+1)/2 in order to provide N*(N+1)/2 extra storage locations. 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. According to the numerical properties of the SVD, TOL1 must be >= !!A!! x EPS where !!A!! denotes the L2-norm and EPS is the machine precision. TOL2 - DOUBLE PRECISION. Working precision for the computation of the desired singular subspaces of A. This parameter specifies that elements of matrices used in the computation, which are <= TOL2 in abso- lute value, are considered to be zero. 3.5 MODE PARAMETER MODE - INTEGER. MODE controls the computation of the desired singular sub- space. It has the decimal expansion AB with the following meaning: A = 0, do not compute the left singular subspace. A = 1, return the M - RANK base vectors of the desired left singular subspace in U. A >= 2, return the first min(M,N) - RANK base vectors of the desired left singular subspace in U. B = 0, do not compute the right singular subspace. B = 1, return the N - RANK base vectors of the desired right singular subspace in V. B >= 2, return the first min(M,N) - RANK base vectors of the desired right singular subspace in V. 3.6 ERROR INDICATORS IERR - INTEGER. On return, IERR contains 0 unless the routine 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 A smaller than 1. 2: number N of columns of array A smaller than 1. 3: leading dimension LDA of array A smaller than M. 4: leading dimension LDU of array U smaller than M. 5: leading dimension LDV of array V smaller than N. 6: rank of matrix A (RANK) larger than min(M,N). 7: the parameters RANK and THETA are both negative (< 0). 8: tolerance TOL1 is negative. 9: tolerance TOL2 is negative. 10: maximum number of QR/QL iteration steps (50) exceeded. 11: parameter MODE out of range. Warnings given by the routine. IWARN = 0: no warnings. 1: the rank of matrix A, specified by the user, has been lowered because a singular value of multiplicity > 1 has been found. 5 EXTERNAL SUBROUTINES and FUNCTIONS: DCOPY form BLAS [5]; DQRDC from LINPACK [6]; BIDIAG, INIT, CANCEL, QRQL, RESTOR. 6 METHOD DESCRIPTION: PSVD is an efficient method (see [1]), computing the singular sub- space 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. First, the Householder transformations of the bidiagonalization need only to be applied on the base vectors of the desired singular sub- spaces. Second, the bidiagonal needs only to be partially diagonalized. Third, the convergence rate of the iterative diagonalization can be improved by an appropriate choice between QL and QR iterations. Depending on the gap, the desired numerical accuracy and the dimen- sion of the desired singular subspace, PSVD can be 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 ----------------------- 1.a): If M >= 5*N/3, transform A into upper triangular form R. 1.b): Transform A (or R) into bidiagonal form: !q(1) e(2) 0 ... 0 ! (0) ! 0 q(2) e(3) . ! J = ! . . ! ! . e(N)! ! 0 ... q(N)! if M >= N, or !q(1) e(2) 0 ... 0 0 ! (0) ! 0 q(2) e(3) . . ! J = ! . . . ! ! . e(M) . ! ! 0 ... q(M) e(M+1)! if M < N, using Householder transformations. 1.c): If U is requested, initialize U with the identity matrix. If V is requested, initialize V with the identity matrix. 1.d): If M < N, then cancel e(M+1), and reduce the bidiagonal to M x M. Accumulate the Givens rotations in V (if V is desired). Step 2: Partial diagonalization phase ----------------------------- If the upper bound THETA is not given, then compute THETA such that precisely p - RANK singular values (p=min(M,N)) of the bidiagonal are <= THETA, using a bisection method [4]. Diagonalize the given bidiagonal J partially, using either QL itera- tions (if the upper left diagonal element of the considered subbi- diagonal > the lower right diagonal element) or QR iterations, such that J is splitted into unreduced subbidiagonals 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 ------------------------- 3.a): Apply the Householder transformations of step 1.b) onto the columns of U and/or V associated with the subbidiagonals with all singular values <= THETA, (if U and/or V is desired). 3.b): If M >= 5*N/3 and U is desired, then apply the Householder transformations of step 1.a) onto each computed column of U in step 3.a). NOTE. If M > N (resp.,M < N), then the base vectors of the orthogonal complement of the column (resp.,row) space of the M by N matrix A can also be computed if desired (see MODE) by applying step 3 onto the last M - N (resp.,N - M) columns of U (resp.,V). 7 REFERENCES: [1] S. Van Huffel, J. Vandewalle and A. Haegemans, An efficient and reliable algorithm for computing the singular subspace of a matrix associated with its smallest singular values. J. Comput. and Applied Math., 19 (1987), 313 - 330. [2] 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. [3] T.F.Chan, An improved algorithm for computing the singular value decomposition. ACM Trans. Math. Software, 8 (1982), 72-83. [4] S. Van Huffel and J. Vandewalle, The Partial Total Least Squares Algorithm. J. Comput. and Applied Math., 21 (1988), to appear. [5] 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. [6] J.J. Dongarra, J.R. Bunch, C.B. Moler and G.W. Stewart, LINPACK User's Guide. SIAM, Philadelphia (1979). 8 NUMERICAL ASPECTS: Using 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. 9 EXAMPLE: 9.1 PROGRAM TEXT PARAMETER (IN = 5, IOUT = 6) PARAMETER (LDA = 30, LDU = 30, LDV = 11, LN = 11, LM = 30, * LQ = 22, LW = 107, LINUL = 30) INTEGER M, N, RANK, MODE, IERR, IWARN DOUBLE PRECISION A(LDA,LN), U(LDU,LM), V(LDV,LN), Q(LQ), * WRK(LW), * TOL1, TOL2, THETA INTRINSIC MIN LOGICAL INUL(LINUL) INTEGER I, J, K READ(IN,12) M, N, RANK, THETA WRITE(IOUT,13) M, N, RANK, THETA WRITE(IOUT,9) DO 1 I = 1, M READ(IN,5) (A(I,J), J = 1, N) 1 CONTINUE WRITE(IOUT,6) DO 2 I = 1, M WRITE(IOUT,7) (A(I,J), J = 1, N) 2 CONTINUE MODE = 11 TOL1 = 1.0D-08 TOL2 = 1.0D-10 CALL PSVD(A, LDA, M, N, RANK, THETA, U, LDU, V, LDV, Q, INUL, * WRK, TOL1, TOL2, MODE, IERR, IWARN) WRITE(IOUT,8) IERR, IWARN WRITE(IOUT,18) RANK, THETA K = MIN(M,N) WRITE(IOUT,10) (Q(I), I = 1, K) WRITE(IOUT,11) (Q(K+I), I = 2, K) WRITE(IOUT,14) K = 0 DO 3 I = 1, M IF (INUL(I)) THEN K = K + 1 WRITE(IOUT,15) K, I WRITE(IOUT,7) (U(J,I), J = 1, M) END IF 3 CONTINUE WRITE(IOUT,16) K = 0 DO 4 I = 1, N IF (INUL(I)) THEN K = K + 1 WRITE(IOUT,17) K, I WRITE(IOUT,7) (V(J,I), J = 1, N) END IF 4 CONTINUE STOP 5 FORMAT(4D15.0) 6 FORMAT(/' ',5X,'A(*,I)',8X,'A(*,I+1)',7X,'A(*,I+2)',7X, * 'A(*,I+3)') 7 FORMAT(' ',4D15.6) 8 FORMAT(/' IERR = ',I3,' IWARN = ',I3) 9 FORMAT(/' PARTIAL SVD (PSVD) TEST PROGRAM'/' ',32('-')/) 10 FORMAT(/' DIAGONAL OF THE PARTIALLY DIAGONALIZED BIDIAGONAL=' * /(1X,4D15.6)) 11 FORMAT(/' SUPERDIAGONAL OF THE PARTIALLY DIAGONALIZED ', * 'BIDIAGONAL=',/(1X,4D15.6)) 12 FORMAT(3I3,D12.5) 13 FORMAT(/' M =',I3,' N =',I3,' RANK =',I3,' THETA =',D12.5) 14 FORMAT(/' BASIS OF THE COMPUTED LEFT SINGULAR SUBSPACE :'/1X, * 44('-')) 15 FORMAT(' THE ',I2,'-TH BASE VECTOR U(*,',I2,') =') 16 FORMAT(/' BASIS OF THE COMPUTED RIGHT SINGULAR SUBSPACE :'/1X, * 45('-')) 17 FORMAT(' THE ',I2,'-TH BASE VECTOR V(*,',I2,') =') 18 FORMAT(/' COMPUTED RANK =',I3,4X,' COMPUTED BOUND THETA = ', * D12.5) END 9.2 PROGRAM DATA 6 4 -1 1.00000D-03 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 PR0GRAM RESULTS: M = 6 N = 4 RANK = -1 THETA = 0.10000D-02 PARTIAL SVD (PSVD) TEST PROGRAM -------------------------------- A(*,I) A(*,I+1) A(*,I+2) A(*,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 COMPUTED BOUND THETA = 0.10000D-02 DIAGONAL OF THE PARTIALLY DIAGONALIZED BIDIAGONAL= 0.322802D+01 0.871401D+00 0.369809D+00 0.128626D-03 SUPERDIAGONAL OF THE PARTIALLY DIAGONALIZED BIDIAGONAL= 0.286516D-01 0.167536D-01 -0.739580D-22 BASIS OF THE COMPUTED LEFT SINGULAR SUBSPACE : -------------------------------------------- THE 1-TH BASE VECTOR U(*, 4) = 0.269797D+00 0.153118D+00 -0.536944D+00 -0.186820D+00 0.642075D+00 -0.410236D+00 THE 2-TH BASE VECTOR U(*, 5) = -0.578307D+00 -0.456351D+00 0.180389D+00 0.336878D+00 0.552879D+00 -0.748493D-01 THE 3-TH BASE VECTOR U(*, 6) = 0.484175D+00 -0.742503D+00 0.646079D-01 -0.334913D+00 0.115913D+00 0.290665D+00 BASIS OF THE COMPUTED RIGHT SINGULAR SUBSPACE : --------------------------------------------- THE 1-TH BASE VECTOR V(*, 4) = -0.355483D+00 -0.568663D+00 -0.212821D+00 0.710606D+00 *********************************************************************** 1988, February 15.