19 June 2012 08:48:42 PM SVD_DEMO: C++ version Compiled on Jun 19 2012 at 20:48:18. 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 the user.) We choose a "random" matrix A, with integral values between 0 and 10. The matrix A: Col: 0 1 2 Row 0: 2.000000 1.000000 1.000000 1: 10.000000 3.000000 4.000000 2: 8.000000 1.000000 4.000000 3: 6.000000 0.000000 8.000000 4: 4.000000 6.000000 8.000000 The orthogonal factor U: Col: 0 1 2 3 4 Row 0: -0.122437 -0.045288 0.140845 -0.229960 -0.954064 1: -0.552266 -0.468282 0.415306 0.550237 0.021788 2: -0.447998 -0.400116 -0.056977 -0.757087 0.250556 3: -0.485998 0.124504 -0.821077 0.242781 -0.123271 4: -0.493068 0.776574 0.360930 -0.110611 0.106358 The diagonal factor S: Col: 0 1 2 Row 0: 19.303064 0.000000 0.000000 1: 0.000000 6.203906 0.000000 2: 0.000000 0.000000 4.111361 3: 0.000000 0.000000 0.000000 4: 0.000000 0.000000 0.000000 The orthogonal factor V: Col: 0 1 2 Row 0: -0.737695 -0.664259 0.120690 1: -0.268643 0.452811 0.850172 2: -0.619384 0.594746 -0.512485 The product U * S * V': Col: 0 1 2 Row 0: 2.000000 1.000000 1.000000 1: 10.000000 3.000000 4.000000 2: 8.000000 1.000000 4.000000 3: 6.000000 0.000000 8.000000 4: 4.000000 6.000000 8.000000 Frobenius Norm of A, A_NORM = 20.6882 ABSOLUTE ERROR for A = U*S*V' Frobenius norm of difference A-U*S*V' = 2.03228e-14 RELATIVE ERROR for A = U*S*V': Ratio of DIF_NORM / A_NORM = 9.82342e-16 RANK_ONE_TEST: Compare A to the sum of R rank one matrices. R Absolute Relative Error Error 0 20.6882 1 1 7.44256 0.35975 2 4.11136 0.19873 3 2.03228e-14 9.82342e-16 RANK_ONE_PRINT_TEST: Print the sums of R rank one matrices. Rank R = 0 Col: 0 1 2 Row 0: 0.000000 0.000000 0.000000 1: 0.000000 0.000000 0.000000 2: 0.000000 0.000000 0.000000 3: 0.000000 0.000000 0.000000 4: 0.000000 0.000000 0.000000 Rank R = 1 Col: 0 1 2 Row 0: 1.743481 0.634916 1.463864 1: 7.864136 2.863849 6.602895 2: 6.379395 2.323157 5.356275 3: 6.920497 2.520208 5.810596 4: 7.021170 2.556869 5.895123 Rank R = 2 Col: 0 1 2 Row 0: 1.930113 0.507694 1.296763 1: 9.793926 1.548354 4.875054 2: 8.028272 1.199154 3.879950 3: 6.407417 2.869963 6.269983 4: 3.820907 4.738418 8.760483 Rank R = 3 Col: 0 1 2 Row 0: 2.000000 1.000000 1.000000 1: 10.000000 3.000000 4.000000 2: 8.000000 1.000000 4.000000 3: 6.000000 0.000000 8.000000 4: 4.000000 6.000000 8.000000 Original matrix A: Col: 0 1 2 Row 0: 2.000000 1.000000 1.000000 1: 10.000000 3.000000 4.000000 2: 8.000000 1.000000 4.000000 3: 6.000000 0.000000 8.000000 4: 4.000000 6.000000 8.000000 The pseudoinverse of A: Col: 0 1 2 3 4 Row 0: 0.013663 0.083436 0.058289 -0.018860 -0.053710 1: 0.027523 0.059386 -0.034751 -0.153936 0.138178 2: -0.017969 -0.078940 -0.016880 0.129878 0.045278 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 9.74976e-15 A+ * A * A+ = A+ 7.33729e-17 ( A * A+ ) is MxM symmetric; 1.23645e-15 ( A+ * A ) is NxN symmetric; 3.51083e-16 In some cases, the matrix A * A+ may be interesting (if M <= N, then it MIGHT look like the identity.) A * A+: Col: 0 1 2 3 4 Row 0: 0.036879 0.147319 0.064947 -0.061779 0.076036 1: 0.147319 0.696764 0.411118 -0.130901 0.058545 2: 0.064947 0.411118 0.364041 0.214692 -0.110391 3: -0.061779 -0.130901 0.214692 0.925862 0.039965 4: 0.076036 0.058545 -0.110391 0.039965 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: 0 1 2 Row 0: 1.000000 0.000000 -0.000000 1: -0.000000 1.000000 0.000000 2: -0.000000 0.000000 1.000000 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 = 1.94793e-14 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.17232 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 = 8.59287e-14 SVD_DEMO: Normal end of execution. 19 June 2012 08:48:42 PM