19 June 2012 11:06:45.104 AM SVD_DEMO FORTRAN90 version Demonstrate the singular value decomposition (SVD) A real MxN matrix A can be factored as: A = U * S * V' where U = MxM orthogonal, S = MxN zero except for diagonal, V = NxN orthogonal. The diagonal of S contains only nonnegative numbers and these are arranged in descending order. Matrix row order M = 5 Matrix column order N = 8 Random number SEED = 123456789 (Chosen by user.) We choose a "random" matrix A, with integral values between 0 and 10. The matrix A: Col 1 2 3 4 5 Row 1 2. 1. 1. 0. 9. 2 10. 3. 4. 9. 8. 3 8. 1. 4. 4. 1. 4 6. 0. 8. 1. 0. 5 4. 6. 8. 0. 3. Col 6 7 8 Row 1 9. 7. 6. 2 1. 6. 2. 3 4. 9. 8. 4 8. 5. 4. 5 3. 9. 2. The orthogonal factor U: Col 1 2 3 4 5 Row 1 -0.421671 0.487547 0.762151 -0.590472E-01 -0.114523E-01 2 -0.487891 -0.791659 0.229657 -0.328621E-01 -0.285335 3 -0.490789 -0.175013E-01 -0.209319 0.530914 0.658131 4 -0.406607 0.350797 -0.437254 0.284878 -0.662771 5 -0.421845 0.110495 -0.362461 -0.795241 0.214595 The diagonal factor S: Col 1 2 3 4 5 Row 1 30.7103 0. 0. 0. 0. 2 0. 11.2457 0. 0. 0. 3 0. 0. 9.73453 0. 0. 4 0. 0. 0. 7.53857 0. 5 0. 0. 0. 0. 5.17902 Col 6 7 8 Row 1 0. 0. 0. 2 0. 0. 0. 3 0. 0. 0. 4 0. 0. 0. 5 0. 0. 0. The orthogonal factor V: Col 1 2 3 4 5 Row 1 -0.448566 -0.403244 -0.197960 0.308931 -0.140848 2 -0.159790 -0.110439 -0.958408E-01 -0.583422 0.208195 3 -0.357014 0.836975E-01 -0.570569 -0.285166 -0.406579 4 -0.220147 -0.608602 0.813995E-01 0.280262 -0.115517 5 -0.307861 -0.145066 0.760172 -0.351410 -0.209274 6 -0.350517 0.592593 0.171177 0.192698 -0.466163 7 -0.525094 0.111491 -0.636183E-01 -0.207607 0.530699 8 -0.322440 0.251308 0.907830E-01 0.447875 0.464135 Col 6 7 8 Row 1 -0.508569 0.374980 0.287168 2 0.343320 0.165197 0.651307 3 -0.474527E-01 -0.524943 -0.124675 4 0.661737 -0.116457 -0.176534 5 -0.269805 -0.255718 -0.290740E-01 6 0.331972 0.350015 0.975347E-01 7 0.170335E-01 0.275777 -0.553818 8 0.708900E-02 -0.528221 0.359935 The product U * S * V': Col 1 2 3 4 5 Row 1 2.00000 1.00000 1.00000 0.412257E-14 9.00000 2 10.0000 3.00000 4.00000 9.00000 8.00000 3 8.00000 1.00000 4.00000 4.00000 1.00000 4 6.00000 0.388578E-14 8.00000 1.00000 0.621725E-14 5 4.00000 6.00000 8.00000 0.249800E-14 3. Col 6 7 8 Row 1 9.00000 7.00000 6.00000 2 1.00000 6.00000 2.00000 3 4.00000 9. 8.00000 4 8.00000 5.00000 4.00000 5 3.00000 9.00000 2.00000 Frobenius Norm of A, A_NORM = 35.3270 ABSOLUTE ERROR for A = U*S*V': Frobenius norm of difference A-U*S*V' = 0.272896E-13 RELATIVE ERROR for A = U*S*V': Ratio of DIF_NORM / A_NORM = 0.772485E-15 RANK_ONE_TEST: Compare A to the sum of R rank one matrices. R Absolute Relative Error Error 0 35.3270 1.00000 1 17.4608 0.494261 2 13.3571 0.378100 3 9.14616 0.258900 4 5.17902 0.146602 5 0.272896E-13 0.772485E-15 RANK_ONE_PRINT_TEST: Print the sums of R rank one matrices. Rank R = 0 Col 1 2 3 4 5 Row 1 0. 0. 0. 0. 0. 2 0. 0. 0. 0. 0. 3 0. 0. 0. 0. 0. 4 0. 0. 0. 0. 0. 5 0. 0. 0. 0. 0. Col 6 7 8 Row 1 0. 0. 0. 2 0. 0. 0. 3 0. 0. 0. 4 0. 0. 0. 5 0. 0. 0. Rank R = 1 Col 1 2 3 4 5 Row 1 5.80876 2.06922 4.62320 2.85082 3.98668 2 6.72098 2.39418 5.34924 3.29852 4.61276 3 6.76090 2.40840 5.38101 3.31811 4.64016 4 5.60125 1.99530 4.45805 2.74898 3.84426 5 5.81116 2.07008 4.62511 2.85200 3.98833 Col 6 7 8 Row 1 4.53907 6.79978 4.17549 2 5.25189 7.86762 4.83121 3 5.28309 7.91435 4.85991 4 4.37692 6.55686 4.02632 5 4.54094 6.80258 4.17721 Rank R = 2 Col 1 2 3 4 5 Row 1 3.59786 1.46371 5.08210 -0.486015 3.19131 2 10.3110 3.37739 4.60410 8.71675 5.90425 3 6.84026 2.43014 5.36454 3.43789 4.66871 4 4.01048 1.55963 4.78823 0.348079 3.27198 5 5.31009 1.93285 4.72911 2.09576 3.80807 Col 6 7 8 Row 1 7.78814 7.41106 5.55336 2 -0.238205E-01 6.87504 2.59388 3 5.16645 7.89241 4.81045 4 6.71467 6.99669 5.01772 5 5.27729 6.94112 4.48948 Rank R = 3 Col 1 2 3 4 5 Row 1 2.12916 0.752649 0.848949 0.117902 8.83116 2 9.86839 3.16313 3.32853 8.89873 7.60369 3 7.24363 2.62542 6.52714 3.27203 3.11977 4 4.85309 1.96757 7.21683 0.160564E-02 0.363447E-01 5 6.00857 2.27101 6.74230 1.80855 1.12589 Col 6 7 8 Row 1 9.05813 6.93906 6.22689 2 0.358863 6.73281 2.79683 3 4.81766 8.02204 4.62546 4 5.98606 7.26748 4.63130 5 4.67331 7.16559 4.16916 Rank R = 4 Col 1 2 3 4 5 Row 1 1.99165 1.01235 0.975885 -0.685149E-02 8.98759 2 9.79186 3.30766 3.39918 8.82929 7.69074 3 8.48008 0.290375 5.38581 4.39374 1.71330 4 5.51654 0.714629 6.60442 0.603489 -0.718334 5 4.15654 5.76861 8.45187 0.128384 3.23259 Col 6 7 8 Row 1 8.97235 7.03148 6.02753 2 0.311125 6.78424 2.68588 3 5.58890 7.19113 6.41801 4 6.39989 6.82163 5.59315 5 3.51809 8.41019 1.48416 Rank R = 5 Col 1 2 3 4 5 Row 1 2.00000 1.00000 1.00000 0.412257E-14 9.00000 2 10.0000 3.00000 4.00000 9.00000 8.00000 3 8.00000 1.00000 4.00000 4.00000 1.00000 4 6.00000 0.388578E-14 8.00000 1.00000 0.621725E-14 5 4.00000 6.00000 8.00000 0.249800E-14 3. Col 6 7 8 Row 1 9.00000 7.00000 6.00000 2 1.00000 6.00000 2.00000 3 4.00000 9. 8.00000 4 8.00000 5.00000 4.00000 5 3.00000 9.00000 2.00000 Original matrix A: Col 1 2 3 4 5 Row 1 2. 1. 1. 0. 9. 2 10. 3. 4. 9. 8. 3 8. 1. 4. 4. 1. 4 6. 0. 8. 1. 0. 5 4. 6. 8. 0. 3. Col 6 7 8 Row 1 9. 7. 6. 2 1. 6. 2. 3 4. 9. 8. 4 8. 5. 4. 5 3. 9. 2. The pseudoinverse of A: Col 1 2 3 4 5 Row 1 -0.289305E-01 0.372563E-01 0.159113E-01 0.319513E-01 -0.288546E-01 2 -0.598832E-02 -0.875057E-03 -0.984532E-02 -0.457147E-01 0.748500E-01 3 -0.330086E-01 0.996223E-02 -0.539056E-01 0.742211E-01 0.402066E-01 4 -0.189293E-01 0.534040E-01 0.777347E-02 0.564769E-02 -0.403379E-01 5 0.606697E-01 0.460988E-01 -0.625423E-01 -0.210926E-01 0.289757E-02 6 0.434277E-01 -0.726661E-02 -0.446685E-01 0.823752E-01 -0.353797E-01 7 0.751516E-02 -0.293409E-01 0.624044E-01 -0.624723E-01 0.545672E-01 8 0.178959E-01 -0.379505E-01 0.933326E-01 -0.344409E-01 -0.244964E-01 PSEUDO_PRODUCT_TEST The following relations MUST hold: A * A+ * A = A A+ * A * A+ = A+ ( A * A+ ) is MxM symmetric; ( A+ * A ) is NxN symmetric Here are the Frobenius norms of the errors in these relationships: A * A+ * A = A 0.203432E-13 A+ * A * A+ = A+ 0.164835E-15 ( A * A+ ) is MxM symmetric; 0.192931E-14 ( A+ * A ) is NxN symmetric; 0.210915E-14 In some cases, the matrix A * A+ may be interesting (if M <= N, then it MIGHT look like the identity.) A * A+: Col 1 2 3 4 5 Row 1 1.00000 -0.277556E-15 -0.333067E-15 -0.166533E-15 -0.166533E-15 2 -0.513478E-15 1.00000 0.138778E-15 -0.444089E-15 0.631439E-15 3 0.138778E-15 0. 1.00000 -0.222045E-15 0.832667E-16 4 0.416334E-16 -0.222045E-15 0. 1.00000 0.319189E-15 5 -0.763278E-16 -0.124900E-15 -0.194289E-15 -0.555112E-15 1.00000 In some cases, the matrix A+ * A may be interesting (if N <= M, then it MIGHT look like the identity.) A+ * A: Col 1 2 3 4 5 Row 1 0.518282 -0.743781E-01 0.208513 0.430903 -0.329765E-01 2 -0.743781E-01 0.430641 0.184212 -0.929715E-01 0.153809 3 0.208513 0.184212 0.706639 -0.517415E-01 -0.150665 4 0.430903 -0.929715E-01 -0.517415E-01 0.517377 0.143627 5 -0.329765E-01 0.153809 -0.150665 0.143627 0.860968 6 0.957332E-02 -0.235319 0.211651 -0.161698 0.181909 7 0.642907E-01 0.309300 0.765277E-01 -0.769232E-01 0.590150E-01 8 0.983159E-01 -0.149601 -0.232075 -0.266541E-02 -0.122698 Col 6 7 8 Row 1 0.957332E-02 0.642907E-01 0.983159E-01 2 -0.235319 0.309300 -0.149601 3 0.211651 0.765277E-01 -0.232075 4 -0.161698 -0.769232E-01 -0.266541E-02 5 0.181909 0.590150E-01 -0.122698 6 0.757771 -0.481642E-01 0.147426 7 -0.481642E-01 0.616943 0.344889 8 0.147426 0.344889 0.591378 PSEUDO_LINEAR_SOLVE_TEST Given: b = A * x1 so that b is in the range of A, solve A * x = b using the pseudoinverse: x2 = A+ * b. Norm of x1 = 12.9228 Norm of x2 = 12.0214 Norm of residual = 0.182542E-12 Given: b = A' * x1 so that b is in the range of A', solve A' * x = b using the pseudoinverse: x2 = A+' * b. Norm of x1 = 13.3791 Norm of x2 = 13.3791 Norm of residual = 0.127106E-12 For M < N, most systems A'*x=b will not be exactly and uniquely solvable, except in the least squares sense. Here is an example: Given b is NOT in the range of A', solve A' * x = b using the pseudoinverse: x2 = A+' * b. Norm of x2 = 0.827285E-01 Norm of residual = 0.857524 COMPARE_LINPACK_LAPACK: While the singular values should be identical, the orthogonal factors may have some differences. Maximum differences: U(LAPACK) - U(LINPACK): 2.82843 S(LAPACK) - S(LINPACK): 0.141273E-13 V(LAPACK) - V(LINPACK): 3.46695 SVD_DEMO: Normal end of execution. 19 June 2012 11:06:45.108 AM