ALGORITHM 581 An Improved Algorithm for Computing the Singular ...
ALGORITHM 581 An Improved Algorithm for Computing the Singular Value Decomposition [F1] T O N Y F. C H A N D e p a r t m e n t of C o m p u t e r Science, Yale U n i v e r s i t y
Categories and Subject Descriptors D.3.2 [Programming Languages] Language Classifications-FORTRAN, G 1.3 [Numerical Analysis] Numerical Linear Algebra--pseudomverses, G 1 6 [Numerical Analysis]. Optlmlzations--least squares methods General Terms" Algorithms Additional Key Words and Phrases singular value decomposition
DESCRIPTION
T h e set of F O R T R A N subroutines given h e r e is a n i m p l e m e n t a t i o n of the algorithm [1] for c o m p u t i n g the Singular Value D e c o m p o s i t i o n (SVD) of a general m b y n rectangular m a t r i x A defined as A = U W V T, where U is an m × m i n ( m , n ) matrix containing the left singular vectors, W is a diagonal m a t r i x of size min(m, n) containing the singular values, a n d V is an n x m i n ( m , n) m a t r i x containing t h e right singular vectors. N o t e t h a t m is allowed to be greater t h a n or less t h a n n. F o r ease of presentation, we a s s u m e m to be greater t h a n or equal to n in the following discussion. T h e algorithm is an i m p r o v e m e n t of the G o l u b - R e i n s c h algorithm [4], which is i m p l e m e n t e d in subroutines S V D and M I N F I T in E I S P A C K [3] and in subroutine S S V D C in L I N P A C K [2]. I t should be m o r e efficient t h a n the G o l u b - R e i n s c h algorithm w h e n m is a p p r o x i m a t e l y larger tJ lan 2n, as is the case in m a n y least squares applications. T h e algorithm has a h y b r i d nature. W h e n m is a b o a t equal to n, the G o l u b - R e i n s c h algorithm is employed. W h e n the ratio m/n is larger t h a n a threshold value, which is d e t e r m i n e d b y detailed operation counts [1], the i m p r o v e d algorithm is used. T h e i m p r o v e d algorithm first c o m p u t e s the Q R factorization of A using Householder transformations, a n d t h e n uses the G o l u b - R e i n s c h a l g o r i t h m on R. A f u r t h e r i m p r o v e m e n t over the G o l u b - R e i n s c h algorithm is w h e n the left singular Received 7 October 1976; revised 25 January 1982, accepted 29 January 1982 This work was supported by NSF Grant DCR 75-13497 and NASA Ames Contract NCA 2-OR745520. The computing tune was provided by the Stanford Linear Accelerator Center (SLAC) Author's address: Department of Computer Science, Yale University, 10 Hillhouse Avenue, P O. Box 2158 Yale Station, New Haven CT 06520. © 1982 ACM 0098-3500/82/0300-0084 $00.75 ACMTransachonson MathmatlcalSoftware,Vol 8, No 1, March 1982,Pages84-88'
Algorithms
•
85
vectors are to be accumulated and saved. Here, instead of accumulating the Givens transformations (in the second phase of the algorithm where the singular values of the bidiagonal matrix are computed) on the m × n matrix containing the left singular vectors, we accumulate t h e m on a t e m p o r a r y n × n matrix. T h i s requires a small overhead in storage of an n × n matrix (small c o m p a r e d with m × n) but offers big savings in time. An additional feature of the new algorithm is t h a t it can accumulate all the left orthogonal transformations on a n u m b e r of given vectors, which can t h e n be used in computing least squares solutions. In this fashion, it is similar to the E I S P A C K routine M I N F I T . T h e r e are three main routines in the package: HYBSVD: MGNSVD: GRSVD:
This is the main routine which implements the hybrid algorithm. This performs the same thing as H Y B S V D , except t h a t it assumes m .ge. n. This is a slightly modified version of routine SVD in E I S P A C K which implements the G o l u b - R e i n s c h algorithm.
Besides, there are two utility routines: SSWAP: SRELPR:
BLAS routine for swapping two vectors. Routine for computing the machine relative precision.
T h e s e five routines must be used together. T h e y have been tested extensively on the IBM 370/168, 360/91 at the Stanford Linear Accelerator Center, and on the D E C 2060 in the C o m p u t e r Science D e p a r t m e n t at Yale. T h e y produce results t h a t agree (up to machine precision) with those produced by SVD, M I N F I T , and SSVDC. T h e y have been verified by P F O R T verifier [5] for portability. REFERENCES
l.
CHAN,
T.F. An improved algorithm for computing the Singular Value Decomposition ACM
Trans, Math. Softw. 8, 1 (Mar. 1982), 72-83.
J J., BUNCH,J.R., MOLER, C.B., AND STEWART, G.W. LINPACK User's Guide. SIAM, Philadelphia, 1979. 3. GARBOW,B.S., ET AL. Matrix E~gensystem Routmes--EISPACK Guide Extension, Lecture Notes in Computer Science Series, No. 51. Sprmger-Verlag,New York, 1977. 4. GOLUB,G.H., AND REINSCH, C Singular Value Decomposition and least squares solutions. In Handbook for Automatw Computatmn, II, Linear Algebra, J.H. Wilkinson an¢~ C. Remsch (Eds.), Springer-Ver|ag, New York, 1971. 5. RYDER,B.G. The PFORT verifier In Softw. Pract. Exper. 4 (1974) 359-377. 2. DONGARRA,
ALGORITHM
[A part of the listing is printed here. T h e complete listing is available from the ACM Algorithms Distribution Service (see page 91 for order form).]
ACM Transactions on Mathematical Software, Vol 8, No 1, March 1982
86
•
Algorithms
SUBROUTINE HYBSVD(NA, NU, NV, NZ, NB, M, N, A, W, MATU, U, MATV, * V, Z, B, IRHS, IERR, RVI) INTEGER NA, NU, NV, NZ, M, N, IRHS, IERR, MIN@ REAL A(NA, I), W(1), U(NU,I), V(NV,I), Z(NZ,I), B(NB,IRHS), RVI(1) LOGICAL MATU, MATV C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C c C C C C C C C
HYB HYB HYB HYB HYB HYB THIS ROUTINE IS A MODIFICATION OF THE GOLUB-REINSCH PROCEDURE (i) HYB HYB T HYB FOR COMPUTING THE SINGULAR VALUE DECOMPOSITION A - UWV OF A REAL M BY N RECTANGULAR MATRIX. U IS M BY MIN(M,N) CONTAINING HYB THE LEFT SINGULAR VECTORS, W IS A MIN(M,N) BY MIN(M,N) DIAGONAL HYB MATRIX CONTAINING THE SINGULAR VALUES, AND V IS N BY MIN(M,N) HYB CONTAINING THE RIGHT SINGULAR VECTORS. HYB HYB THE ALGORITHM IMPLEMENTED IN THIS HYB ROUTINE HAS A HYBRID NATURE. WHEN M IS APPROXIMATELY EQUAL TO N, HYB THE GOLUB-REINSCH ALGORITHM IS USED, BUT WHEN EITHER OF THE RATIOSHYB HYB M/N OR N/M IS GREATER THAN ABOUT 2, A MODIFIED VERSION OF THE GOLUB-REINSCH HYB ALGORITHM IS USED. THIS MODIFIED ALGORITHM FIRST TRANSFORMS A HYB T HYB HYB INTO UPPER TRIANGULAR FORM BY HOUSEHOLDER TRANSFORMATIONS L AND THEN USES THE GOLUB-REINSCH ALGORITHM TO FIND THE SINGULAR HYB HYB VALUE DECOMPOSITION OF THE RESULTING UPPER TRIANGULAR MATRIX R. WHEN U IS NEEDED EXPLICITLY IN THE CASE M.GE.N (OR V IN THE CASE HYB HYB M.LT.N), AN EXTRA ARRAY Z (OF SIZE AT LEAST HYB MIN(M,N)**2) IS NEEDED, BUT OTHERWISE Z IS NOT REFERENCED HYB AND NO EXTRA STORAGE IS REQUIRED. THIS HYBRID METHOD SHOULD BE MORE EFFICIENT THAN THE GOLUB-REINSCH ALGORITHM WHEN HYB HYB M/N OR N/M IS LARGE. FOR DETAILS, SEE (2). HYB HYB WHEN M .GE. N, HYB HYBSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST HYB SQUARES SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X-B. HYB IF M .LT. N (I.E. FOR UNDERDETERMINED SYSTEMS), THE RHS B HYB IS NOT PROCESSED. HYB h'YB NOTICE THAT THE SINGULAR VALUE DECOMPOSITION OF A MATRIX IS UNIQUE ONLY UP TO THE SIGN OF THE CORRESPONDING COLUMNS HYB RYE OF U AND V. ~B ~B THIS ROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER (3) FOR HYB ADHERENCE TO A LARGE, CAREFULLY DEFINED, PORTABLE SUBSET OF HYB AMERICAN NATIONAL STANDARD FORTRAN CALLED PFORT. HYB REFERENCES: HYB HYB HYB (i) GOLUB,G.H. AND REINSCH,C. (197@) 'SINGULAR VALUE HYB DECOMPOSITION AND LEAST SQUARES SOLUTIONS,' HYB NUMER. MATH. 14,4@3-42~, 197@. HYB HYB (2) CHAN,T.F. (1982) 'AN IMPROVED ALGORITHM FOR COMPUTING HYB THE SINGULAR VALUE DECOMPOSITION,' ACM TOMS, VOL.8, I-~YB NO. i, MARCH, 1982. HYB HYB (3) RYDER,B.G. (1974) 'THE PFORT VERIFIER,' SOFTWARE HYB PRACTICE AND EXPERIENCE, VOL.4, 359-377, 1974. HYB HYB ON INPUT: HYB
ACM Transactions on Mathematmal Software, Vol 8, No. 1, March 1982.
i@ 2@ 3@ 40 5@ 6@ 7@ 8@ 9@ i@@
11@ 12@ 13@ 14@ 15@ 16@ 17@ 18@ 19@ 2@@ 21@ 22@ 23@ 24@ 25@ 26@ 27@ 28@
29@
3@@ 31@ 32@ 33@ 34@ 35@ 36@ 37@ 38@ 39@
4@@ 41¢ 42@ 43@ 44~ 450 46@ 47@ 48~ 49@ 5@@ 51@ 52@ 53@ 540 55@ 56@ 57@ 58~ 59@ 6@@
Algorithms
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C c C c C C c C C C C C C C c C C C c C c c C
NA MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER A AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NOTE THAT NA MUST BE AT LEAST AS LARGE AS M. NU MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY U AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NU MUST BE AT LEAST AS LARGE AS H. NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NV MUST BE AT LEAST AS LARGE AS N.
•
HYB HYB HYB HYB BYB HYB HYB HYB HYB HYB HYB HYB
87
61@ 62@ 63@ 64@ 65@
66@ 67@ 68@ 69@ 7@@ 71~ 72@
HYB
730
74@ 75@ 76@ 77@
NB MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER B AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NB MUST BE AT LEAST AS LARGE AS H.
HYB HYB HYB BYE HYB BYB HYB HYB
M IS THE NUMBER OF ROWS OF A (AND U).
BYB HYB
N IS THE NUMBER OF COLUMNS OF A (AND NUMBER OF ROWS OF V).
HYB HYB
NZ MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER Z AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NOTE THAT NZ MUST BE AT LEAST " AS LARGE AS MIN(M,N).
A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED. B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED LINEAR SYSTEM A*X=B. IF IRES .GT. @ AND M .GE. N, THEN ON OUTPUT, THE FIRST N COMPONENTS OF THESE IRHS COLUMNS T WILL CONTAIN U B. THUS, TO COMPUTE THE MINIMAL LENGTH LEAST + SQUARES SOLUTION, ONE MUST COMPUTE V*W TIMES THE COLUMNS OF + + B, WHERE W IS A DIAGONAL MATRIX, W (I)=@ IF W(1) IS NEGLIGIBLE, OTHERWISE IS I/W(1). IF IRHS=@ OR M.LT.N, B IS NOT REFERENCED.
78@
79@ 8@@ 81@ 82@ 83@ 84@
85@
IIYB
86@
HYB BYB HYB HYB HYB HYB HYB HYB HYB HYB BYB HYB HYB
87@ 88@ 89@ 9@@ 91~ 92@ 93@ 94@ 95@ 96@ 97@ 98@ 99@
BYB 1@@@ IRHS IS THE NUMBER OF RIGHT-HAND-SIDES
OF THE OVERDETERMINED
HYB 1@1@
SYSTEM A*X=B. IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULAR HYB VALUE DECOMPOSITION OF A IS DESIRED. HYB HYB MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE HYB DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE. HYB HYB MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE HYB DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE. HYB HYB WHEN HYBSVD IS USED TO COMPUTE THE MINIMAL LENGTH LEAST HYB SQUARES SOLUTION TO AN OVERDETERMINED SYSTEM, MATU SHOULD HYB BE SET TO .FALSE. , AND MATV SHOULD BE SET TO .TRUE. HYB HYB
ON OUTPUT: A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V). W CONTAINS THE (NON-NEGATIVE) SINGULAR VALUES OF A (THE DIAGONAL ELEMENTS OF W). THEY ARE SORTED IN DESCENDING
HYB HYB HYB HYB HYB HYB
1@20 i@3@ 1@40 1@50 1@6@ 1@7@ 1@8@ 1090 Ii@@ iii@ 112@ 113@ 114@ 115@ 116@ 117@ 118@ 119@ 12@¢
ACM Transacttons on Mathematical Software, Vol. 8, No. 1, March 1982
88
•
Algorithms
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
ORDER. IF AN ERROR EXIT IS MADE, THE SINGULAR VALUES HYB 121@ SHOULD BE CORRECT AND SORTED FOR INDICES IERR+I,...,MIN(M,N) .HYB 122@ HYB 123@ CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE HYB 124@ DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE. IF MATU IS HYB 125@ FALSE, THEN U IS EITHER USED AS A TEMPORARY STORAGE (IF HYB 126@ M .GE. N) OR NOT REFERENCED (IF M .LT. N). HYB 127@ U MAY COINCIDE WITH A IN THE CALLING SEQUENCE. HYB 128@ IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING HYB 129@ TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT. HYB 13~@ HYB 131@ CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF HYB 132@ MATV HAS BEEN SET TO .TRUE. IF MATV IS HYB 133@ FALSE, THEN V IS EITHER USED AS A TEMPORARY STORAGE (IF HYB 134@ M .LT. N) OR NOT REFERENCED (IF M .GE. N). HYB 135@ IF M .GE. N, V MAY ALSO COINCIDE WITH A. IF AN ERROR. HYB 136@ EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF HYB 137@ CORRECT SINGULAR VALUES SHOULD BE CORRECT. HYB 138~ BYB 139@ CONTAINS THE MATRIX X IN THE SINGULAR VALUE DECOMPOSITION HYB 14@@ T HYB 141@ OF R-XSY, IF THE MODIFIED ALGORITHM IS USED. IF THE HYB 142@ GOLUB-REINSCH PROCEDURE IS USED, THEN IT IS NOT REFERENCED. HYB 143@ IF MATU HAS BEEN SET TO .FALSE. IN THE CASE M.GE.N (OR HYB 144@ MATV SET TO .FALSE. IN THE CASE M.LT.N), THEN Z IS NOT HYB 145@ REFERENCED AND NO EXTRA STORAGE IS REQUIRED. HYB 146@ HYB 147~ IERR IS SET TO HYB 148@ ZERO FOR NORMAL RETURN, HYB 149@ K BYB 15@@ IF THE K-TH SINGULAR VALUE HAS NOT BEEN DETERMINED AFTER 3@ ITERATIONS. HYB 151@ -i IF IRMS .LT. @ . HYB 152@ -2 IF M .LT. 1 .OR. N .LT. i m~B 153@ IF NA .LT. M .OR. NU .LT. M .OR. NB .LT. M. HYB 154@ -3 -4 IF NV .LT. N . HYB 155@ -5 IF NZ .LT. MIN(M,N). HYB 156@ HYB 157@ RVI IS A TEMPORARY STORAGE.ARRAY OF LENGTH AT LEAST MIN(M,N). HYB 158@ HYB 159@ PROGRAMMED BY : TONY CHAN HYB 16@@ BOX 2158, YALE STATION, HYB 161@ HYB 162@ COMPUTER SCIENCE DEPT, YALE UNIV., NEW HAVEN, CT @652@. HYB 163@ HYB 164@ LAST MODIFIED : JANUARY, 1982. HYB 165@ HYB 166@ HYBSVD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES. HYB 167@ INTERNAL GRSVD, MGNSVD, SRELPR HYB 168@ FORTRAN MIN@,ABS,SQRT,FLOAT,SIGN,AMAXI HYB 169@ BLAS SSWAP HYB 17@@ .................................................................
ACM Transactions on Mathematical Software, Vol. 8, No 1, March 1982
HYB 171@
Categories and Subject Descriptors D.3.2 [Programming Languages] Language Classifications-FORTRAN, G 1.3 [Numerical Analysis] Numerical Linear Algebra--pseudomverses, G 1 6 [Numerical Analysis]. Optlmlzations--least squares methods General Terms" Algorithms Additional Key Words and Phrases singular value decomposition
DESCRIPTION
T h e set of F O R T R A N subroutines given h e r e is a n i m p l e m e n t a t i o n of the algorithm [1] for c o m p u t i n g the Singular Value D e c o m p o s i t i o n (SVD) of a general m b y n rectangular m a t r i x A defined as A = U W V T, where U is an m × m i n ( m , n ) matrix containing the left singular vectors, W is a diagonal m a t r i x of size min(m, n) containing the singular values, a n d V is an n x m i n ( m , n) m a t r i x containing t h e right singular vectors. N o t e t h a t m is allowed to be greater t h a n or less t h a n n. F o r ease of presentation, we a s s u m e m to be greater t h a n or equal to n in the following discussion. T h e algorithm is an i m p r o v e m e n t of the G o l u b - R e i n s c h algorithm [4], which is i m p l e m e n t e d in subroutines S V D and M I N F I T in E I S P A C K [3] and in subroutine S S V D C in L I N P A C K [2]. I t should be m o r e efficient t h a n the G o l u b - R e i n s c h algorithm w h e n m is a p p r o x i m a t e l y larger tJ lan 2n, as is the case in m a n y least squares applications. T h e algorithm has a h y b r i d nature. W h e n m is a b o a t equal to n, the G o l u b - R e i n s c h algorithm is employed. W h e n the ratio m/n is larger t h a n a threshold value, which is d e t e r m i n e d b y detailed operation counts [1], the i m p r o v e d algorithm is used. T h e i m p r o v e d algorithm first c o m p u t e s the Q R factorization of A using Householder transformations, a n d t h e n uses the G o l u b - R e i n s c h a l g o r i t h m on R. A f u r t h e r i m p r o v e m e n t over the G o l u b - R e i n s c h algorithm is w h e n the left singular Received 7 October 1976; revised 25 January 1982, accepted 29 January 1982 This work was supported by NSF Grant DCR 75-13497 and NASA Ames Contract NCA 2-OR745520. The computing tune was provided by the Stanford Linear Accelerator Center (SLAC) Author's address: Department of Computer Science, Yale University, 10 Hillhouse Avenue, P O. Box 2158 Yale Station, New Haven CT 06520. © 1982 ACM 0098-3500/82/0300-0084 $00.75 ACMTransachonson MathmatlcalSoftware,Vol 8, No 1, March 1982,Pages84-88'
Algorithms
•
85
vectors are to be accumulated and saved. Here, instead of accumulating the Givens transformations (in the second phase of the algorithm where the singular values of the bidiagonal matrix are computed) on the m × n matrix containing the left singular vectors, we accumulate t h e m on a t e m p o r a r y n × n matrix. T h i s requires a small overhead in storage of an n × n matrix (small c o m p a r e d with m × n) but offers big savings in time. An additional feature of the new algorithm is t h a t it can accumulate all the left orthogonal transformations on a n u m b e r of given vectors, which can t h e n be used in computing least squares solutions. In this fashion, it is similar to the E I S P A C K routine M I N F I T . T h e r e are three main routines in the package: HYBSVD: MGNSVD: GRSVD:
This is the main routine which implements the hybrid algorithm. This performs the same thing as H Y B S V D , except t h a t it assumes m .ge. n. This is a slightly modified version of routine SVD in E I S P A C K which implements the G o l u b - R e i n s c h algorithm.
Besides, there are two utility routines: SSWAP: SRELPR:
BLAS routine for swapping two vectors. Routine for computing the machine relative precision.
T h e s e five routines must be used together. T h e y have been tested extensively on the IBM 370/168, 360/91 at the Stanford Linear Accelerator Center, and on the D E C 2060 in the C o m p u t e r Science D e p a r t m e n t at Yale. T h e y produce results t h a t agree (up to machine precision) with those produced by SVD, M I N F I T , and SSVDC. T h e y have been verified by P F O R T verifier [5] for portability. REFERENCES
l.
CHAN,
T.F. An improved algorithm for computing the Singular Value Decomposition ACM
Trans, Math. Softw. 8, 1 (Mar. 1982), 72-83.
J J., BUNCH,J.R., MOLER, C.B., AND STEWART, G.W. LINPACK User's Guide. SIAM, Philadelphia, 1979. 3. GARBOW,B.S., ET AL. Matrix E~gensystem Routmes--EISPACK Guide Extension, Lecture Notes in Computer Science Series, No. 51. Sprmger-Verlag,New York, 1977. 4. GOLUB,G.H., AND REINSCH, C Singular Value Decomposition and least squares solutions. In Handbook for Automatw Computatmn, II, Linear Algebra, J.H. Wilkinson an¢~ C. Remsch (Eds.), Springer-Ver|ag, New York, 1971. 5. RYDER,B.G. The PFORT verifier In Softw. Pract. Exper. 4 (1974) 359-377. 2. DONGARRA,
ALGORITHM
[A part of the listing is printed here. T h e complete listing is available from the ACM Algorithms Distribution Service (see page 91 for order form).]
ACM Transactions on Mathematical Software, Vol 8, No 1, March 1982
86
•
Algorithms
SUBROUTINE HYBSVD(NA, NU, NV, NZ, NB, M, N, A, W, MATU, U, MATV, * V, Z, B, IRHS, IERR, RVI) INTEGER NA, NU, NV, NZ, M, N, IRHS, IERR, MIN@ REAL A(NA, I), W(1), U(NU,I), V(NV,I), Z(NZ,I), B(NB,IRHS), RVI(1) LOGICAL MATU, MATV C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C c C C C C C C C
HYB HYB HYB HYB HYB HYB THIS ROUTINE IS A MODIFICATION OF THE GOLUB-REINSCH PROCEDURE (i) HYB HYB T HYB FOR COMPUTING THE SINGULAR VALUE DECOMPOSITION A - UWV OF A REAL M BY N RECTANGULAR MATRIX. U IS M BY MIN(M,N) CONTAINING HYB THE LEFT SINGULAR VECTORS, W IS A MIN(M,N) BY MIN(M,N) DIAGONAL HYB MATRIX CONTAINING THE SINGULAR VALUES, AND V IS N BY MIN(M,N) HYB CONTAINING THE RIGHT SINGULAR VECTORS. HYB HYB THE ALGORITHM IMPLEMENTED IN THIS HYB ROUTINE HAS A HYBRID NATURE. WHEN M IS APPROXIMATELY EQUAL TO N, HYB THE GOLUB-REINSCH ALGORITHM IS USED, BUT WHEN EITHER OF THE RATIOSHYB HYB M/N OR N/M IS GREATER THAN ABOUT 2, A MODIFIED VERSION OF THE GOLUB-REINSCH HYB ALGORITHM IS USED. THIS MODIFIED ALGORITHM FIRST TRANSFORMS A HYB T HYB HYB INTO UPPER TRIANGULAR FORM BY HOUSEHOLDER TRANSFORMATIONS L AND THEN USES THE GOLUB-REINSCH ALGORITHM TO FIND THE SINGULAR HYB HYB VALUE DECOMPOSITION OF THE RESULTING UPPER TRIANGULAR MATRIX R. WHEN U IS NEEDED EXPLICITLY IN THE CASE M.GE.N (OR V IN THE CASE HYB HYB M.LT.N), AN EXTRA ARRAY Z (OF SIZE AT LEAST HYB MIN(M,N)**2) IS NEEDED, BUT OTHERWISE Z IS NOT REFERENCED HYB AND NO EXTRA STORAGE IS REQUIRED. THIS HYBRID METHOD SHOULD BE MORE EFFICIENT THAN THE GOLUB-REINSCH ALGORITHM WHEN HYB HYB M/N OR N/M IS LARGE. FOR DETAILS, SEE (2). HYB HYB WHEN M .GE. N, HYB HYBSVD CAN ALSO BE USED TO COMPUTE THE MINIMAL LENGTH LEAST HYB SQUARES SOLUTION TO THE OVERDETERMINED LINEAR SYSTEM A*X-B. HYB IF M .LT. N (I.E. FOR UNDERDETERMINED SYSTEMS), THE RHS B HYB IS NOT PROCESSED. HYB h'YB NOTICE THAT THE SINGULAR VALUE DECOMPOSITION OF A MATRIX IS UNIQUE ONLY UP TO THE SIGN OF THE CORRESPONDING COLUMNS HYB RYE OF U AND V. ~B ~B THIS ROUTINE HAS BEEN CHECKED BY THE PFORT VERIFIER (3) FOR HYB ADHERENCE TO A LARGE, CAREFULLY DEFINED, PORTABLE SUBSET OF HYB AMERICAN NATIONAL STANDARD FORTRAN CALLED PFORT. HYB REFERENCES: HYB HYB HYB (i) GOLUB,G.H. AND REINSCH,C. (197@) 'SINGULAR VALUE HYB DECOMPOSITION AND LEAST SQUARES SOLUTIONS,' HYB NUMER. MATH. 14,4@3-42~, 197@. HYB HYB (2) CHAN,T.F. (1982) 'AN IMPROVED ALGORITHM FOR COMPUTING HYB THE SINGULAR VALUE DECOMPOSITION,' ACM TOMS, VOL.8, I-~YB NO. i, MARCH, 1982. HYB HYB (3) RYDER,B.G. (1974) 'THE PFORT VERIFIER,' SOFTWARE HYB PRACTICE AND EXPERIENCE, VOL.4, 359-377, 1974. HYB HYB ON INPUT: HYB
ACM Transactions on Mathematmal Software, Vol 8, No. 1, March 1982.
i@ 2@ 3@ 40 5@ 6@ 7@ 8@ 9@ i@@
11@ 12@ 13@ 14@ 15@ 16@ 17@ 18@ 19@ 2@@ 21@ 22@ 23@ 24@ 25@ 26@ 27@ 28@
29@
3@@ 31@ 32@ 33@ 34@ 35@ 36@ 37@ 38@ 39@
4@@ 41¢ 42@ 43@ 44~ 450 46@ 47@ 48~ 49@ 5@@ 51@ 52@ 53@ 540 55@ 56@ 57@ 58~ 59@ 6@@
Algorithms
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C c C c C C c C C C C C C C c C C C c C c c C
NA MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER A AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NOTE THAT NA MUST BE AT LEAST AS LARGE AS M. NU MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY U AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NU MUST BE AT LEAST AS LARGE AS H. NV MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER V AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NV MUST BE AT LEAST AS LARGE AS N.
•
HYB HYB HYB HYB BYB HYB HYB HYB HYB HYB HYB HYB
87
61@ 62@ 63@ 64@ 65@
66@ 67@ 68@ 69@ 7@@ 71~ 72@
HYB
730
74@ 75@ 76@ 77@
NB MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER B AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NB MUST BE AT LEAST AS LARGE AS H.
HYB HYB HYB BYE HYB BYB HYB HYB
M IS THE NUMBER OF ROWS OF A (AND U).
BYB HYB
N IS THE NUMBER OF COLUMNS OF A (AND NUMBER OF ROWS OF V).
HYB HYB
NZ MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL ARRAY PARAMETER Z AS DECLARED IN THE CALLING PROGRAM DIMENSION STATEMENT. NOTE THAT NZ MUST BE AT LEAST " AS LARGE AS MIN(M,N).
A CONTAINS THE RECTANGULAR INPUT MATRIX TO BE DECOMPOSED. B CONTAINS THE IRHS RIGHT-HAND-SIDES OF THE OVERDETERMINED LINEAR SYSTEM A*X=B. IF IRES .GT. @ AND M .GE. N, THEN ON OUTPUT, THE FIRST N COMPONENTS OF THESE IRHS COLUMNS T WILL CONTAIN U B. THUS, TO COMPUTE THE MINIMAL LENGTH LEAST + SQUARES SOLUTION, ONE MUST COMPUTE V*W TIMES THE COLUMNS OF + + B, WHERE W IS A DIAGONAL MATRIX, W (I)=@ IF W(1) IS NEGLIGIBLE, OTHERWISE IS I/W(1). IF IRHS=@ OR M.LT.N, B IS NOT REFERENCED.
78@
79@ 8@@ 81@ 82@ 83@ 84@
85@
IIYB
86@
HYB BYB HYB HYB HYB HYB HYB HYB HYB HYB BYB HYB HYB
87@ 88@ 89@ 9@@ 91~ 92@ 93@ 94@ 95@ 96@ 97@ 98@ 99@
BYB 1@@@ IRHS IS THE NUMBER OF RIGHT-HAND-SIDES
OF THE OVERDETERMINED
HYB 1@1@
SYSTEM A*X=B. IRHS SHOULD BE SET TO ZERO IF ONLY THE SINGULAR HYB VALUE DECOMPOSITION OF A IS DESIRED. HYB HYB MATU SHOULD BE SET TO .TRUE. IF THE U MATRIX IN THE HYB DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE. HYB HYB MATV SHOULD BE SET TO .TRUE. IF THE V MATRIX IN THE HYB DECOMPOSITION IS DESIRED, AND TO .FALSE. OTHERWISE. HYB HYB WHEN HYBSVD IS USED TO COMPUTE THE MINIMAL LENGTH LEAST HYB SQUARES SOLUTION TO AN OVERDETERMINED SYSTEM, MATU SHOULD HYB BE SET TO .FALSE. , AND MATV SHOULD BE SET TO .TRUE. HYB HYB
ON OUTPUT: A IS UNALTERED (UNLESS OVERWRITTEN BY U OR V). W CONTAINS THE (NON-NEGATIVE) SINGULAR VALUES OF A (THE DIAGONAL ELEMENTS OF W). THEY ARE SORTED IN DESCENDING
HYB HYB HYB HYB HYB HYB
1@20 i@3@ 1@40 1@50 1@6@ 1@7@ 1@8@ 1090 Ii@@ iii@ 112@ 113@ 114@ 115@ 116@ 117@ 118@ 119@ 12@¢
ACM Transacttons on Mathematical Software, Vol. 8, No. 1, March 1982
88
•
Algorithms
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
ORDER. IF AN ERROR EXIT IS MADE, THE SINGULAR VALUES HYB 121@ SHOULD BE CORRECT AND SORTED FOR INDICES IERR+I,...,MIN(M,N) .HYB 122@ HYB 123@ CONTAINS THE MATRIX U (ORTHOGONAL COLUMN VECTORS) OF THE HYB 124@ DECOMPOSITION IF MATU HAS BEEN SET TO .TRUE. IF MATU IS HYB 125@ FALSE, THEN U IS EITHER USED AS A TEMPORARY STORAGE (IF HYB 126@ M .GE. N) OR NOT REFERENCED (IF M .LT. N). HYB 127@ U MAY COINCIDE WITH A IN THE CALLING SEQUENCE. HYB 128@ IF AN ERROR EXIT IS MADE, THE COLUMNS OF U CORRESPONDING HYB 129@ TO INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT. HYB 13~@ HYB 131@ CONTAINS THE MATRIX V (ORTHOGONAL) OF THE DECOMPOSITION IF HYB 132@ MATV HAS BEEN SET TO .TRUE. IF MATV IS HYB 133@ FALSE, THEN V IS EITHER USED AS A TEMPORARY STORAGE (IF HYB 134@ M .LT. N) OR NOT REFERENCED (IF M .GE. N). HYB 135@ IF M .GE. N, V MAY ALSO COINCIDE WITH A. IF AN ERROR. HYB 136@ EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO INDICES OF HYB 137@ CORRECT SINGULAR VALUES SHOULD BE CORRECT. HYB 138~ BYB 139@ CONTAINS THE MATRIX X IN THE SINGULAR VALUE DECOMPOSITION HYB 14@@ T HYB 141@ OF R-XSY, IF THE MODIFIED ALGORITHM IS USED. IF THE HYB 142@ GOLUB-REINSCH PROCEDURE IS USED, THEN IT IS NOT REFERENCED. HYB 143@ IF MATU HAS BEEN SET TO .FALSE. IN THE CASE M.GE.N (OR HYB 144@ MATV SET TO .FALSE. IN THE CASE M.LT.N), THEN Z IS NOT HYB 145@ REFERENCED AND NO EXTRA STORAGE IS REQUIRED. HYB 146@ HYB 147~ IERR IS SET TO HYB 148@ ZERO FOR NORMAL RETURN, HYB 149@ K BYB 15@@ IF THE K-TH SINGULAR VALUE HAS NOT BEEN DETERMINED AFTER 3@ ITERATIONS. HYB 151@ -i IF IRMS .LT. @ . HYB 152@ -2 IF M .LT. 1 .OR. N .LT. i m~B 153@ IF NA .LT. M .OR. NU .LT. M .OR. NB .LT. M. HYB 154@ -3 -4 IF NV .LT. N . HYB 155@ -5 IF NZ .LT. MIN(M,N). HYB 156@ HYB 157@ RVI IS A TEMPORARY STORAGE.ARRAY OF LENGTH AT LEAST MIN(M,N). HYB 158@ HYB 159@ PROGRAMMED BY : TONY CHAN HYB 16@@ BOX 2158, YALE STATION, HYB 161@ HYB 162@ COMPUTER SCIENCE DEPT, YALE UNIV., NEW HAVEN, CT @652@. HYB 163@ HYB 164@ LAST MODIFIED : JANUARY, 1982. HYB 165@ HYB 166@ HYBSVD USES THE FOLLOWING FUNCTIONS AND SUBROUTINES. HYB 167@ INTERNAL GRSVD, MGNSVD, SRELPR HYB 168@ FORTRAN MIN@,ABS,SQRT,FLOAT,SIGN,AMAXI HYB 169@ BLAS SSWAP HYB 17@@ .................................................................
ACM Transactions on Mathematical Software, Vol. 8, No 1, March 1982
HYB 171@
