19 June 2012 11:03:13.856 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 = 3 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 Row 1: 2.00000 1.00000 1.00000 2: 10.0000 3.00000 4.00000 3: 8.00000 1.00000 4.00000 4: 6.00000 0.00000 8.00000 5: 4.00000 6.00000 8.00000 The orthogonal factor U: Col 1 2 3 4 5 Row 1: -0.122437 0.452879E-01 0.140845 0.268239 -0.944017 2: -0.552266 0.468282 0.415306 0.468609 0.289209 3: -0.447998 0.400116 -0.569765E-01 -0.782546 -0.153558 4: -0.485998 -0.124504 -0.821077 0.272025 0.118521E-01 5: -0.493068 -0.776574 0.360930 -0.148586 0.383248E-01 The diagonal factor S: Col 1 2 3 Row 1: 19.3031 0.00000 0.00000 2: 0.00000 6.20391 0.00000 3: 0.00000 0.00000 4.11136 4: 0.00000 0.00000 0.00000 5: 0.00000 0.00000 0.00000 The orthogonal factor V: Col 1 2 3 Row 1: -0.737695 0.664259 0.120690 2: -0.268643 -0.452811 0.850172 3: -0.619384 -0.594746 -0.512485 The product U * S * V': Col 1 2 3 Row 1: 2.00000 1.00000 1.00000 2: 10.0000 3.00000 4.00000 3: 8.00000 1.00000 4.00000 4: 6.00000 0.177636E-14 8.00000 5: 4.00000 6.00000 8.00000 Frobenius Norm of A, A_NORM = 20.6882 ABSOLUTE ERROR for A = U*S*V': Frobenius norm of difference A-U*S*V' = 0.110783E-13 RELATIVE ERROR for A = U*S*V': Ratio of DIF_NORM / A_NORM = 0.535492E-15 RANK_ONE_TEST: Compare A to the sum of R rank one matrices. R Absolute Relative Error Error 0 20.6882 1.00000 1 7.44256 0.359750 2 4.11136 0.198730 3 0.110783E-13 0.535492E-15 RANK_ONE_PRINT_TEST: Print the sums of R rank one matrices. Rank R = 0 Col 1 2 3 Row 1: 0.00000 0.00000 0.00000 2: 0.00000 0.00000 0.00000 3: 0.00000 0.00000 0.00000 4: 0.00000 0.00000 0.00000 5: 0.00000 0.00000 0.00000 Rank R = 1 Col 1 2 3 Row 1: 1.74348 0.634916 1.46386 2: 7.86414 2.86385 6.60289 3: 6.37939 2.32316 5.35627 4: 6.92050 2.52021 5.81060 5: 7.02117 2.55687 5.89512 Rank R = 2 Col 1 2 3 Row 1: 1.93011 0.507694 1.29676 2: 9.79393 1.54835 4.87505 3: 8.02827 1.19915 3.87995 4: 6.40742 2.86996 6.26998 5: 3.82091 4.73842 8.76048 Rank R = 3 Col 1 2 3 Row 1: 2.00000 1.00000 1.00000 2: 10.0000 3.00000 4.00000 3: 8.00000 1.00000 4.00000 4: 6.00000 0.177636E-14 8.00000 5: 4.00000 6.00000 8.00000 Original matrix A: Col 1 2 3 Row 1: 2.00000 1.00000 1.00000 2: 10.0000 3.00000 4.00000 3: 8.00000 1.00000 4.00000 4: 6.00000 0.00000 8.00000 5: 4.00000 6.00000 8.00000 The pseudoinverse of A: Col 1 2 3 4 5 Row 1: 0.136627E-01 0.834365E-01 0.582892E-01 -0.188605E-01 -0.537102E-01 2: 0.275234E-01 0.593865E-01 -0.347508E-01 -0.153936 0.138178 3: -0.179694E-01 -0.789400E-01 -0.168804E-01 0.129878 0.452783E-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.857176E-14 A+ * A * A+ = A+ 0.122872E-15 ( A * A+ ) is MxM symmetric; 0.111472E-14 ( A+ * A ) is NxN symmetric; 0.105909E-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: 0.368794E-01 0.147319 0.649473E-01 -0.617791E-01 0.760358E-01 2: 0.147319 0.696764 0.411118 -0.130901 0.585450E-01 3: 0.649473E-01 0.411118 0.364041 0.214692 -0.110391 4: -0.617791E-01 -0.130901 0.214692 0.925862 0.399650E-01 5: 0.760358E-01 0.585450E-01 -0.110391 0.399650E-01 0.976453 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 Row 1: 1.00000 0.277556E-15 0.222045E-15 2: 0.333067E-15 1.00000 0.222045E-15 3: -0.277556E-15 -0.333067E-15 1.00000 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 = 9.84886 Norm of x2 = 9.84886 Norm of residual = 0.367925E-13 For N < M, 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.172320 Norm of residual = 0.475409 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 = 10.3441 Norm of x2 = 8.44769 Norm of residual = 0.512380E-13 COMPARE_LINPACK_LAPACK: While the singular values should be identical, the orthogonal factors may have some differences. Frobenius differences: U(LAPACK) - U(LINPACK): 2.12500 S(LAPACK) - S(LINPACK): 0.758860E-14 V(LAPACK) - V(LINPACK): 2.00000 SVD_DEMO: Normal end of execution. 19 June 2012 11:03:13.858 AM