19 June 2012 08:48:50 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 = 5 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 3 4 Row 0: 2.000000 1.000000 1.000000 0.000000 9.000000 1: 10.000000 3.000000 4.000000 9.000000 8.000000 2: 8.000000 1.000000 4.000000 4.000000 1.000000 3: 6.000000 0.000000 8.000000 1.000000 0.000000 4: 4.000000 6.000000 8.000000 0.000000 3.000000 The orthogonal factor U: Col: 0 1 2 3 4 Row 0: -0.256248 0.569696 -0.552868 -0.505630 0.220133 1: -0.695099 0.412897 0.343245 0.270593 -0.394100 2: -0.405979 -0.190643 0.396041 -0.095763 0.795498 3: -0.355915 -0.564787 0.003267 -0.631351 -0.394621 4: -0.399600 -0.386825 -0.647812 0.513170 0.087654 The diagonal factor S: Col: 0 1 2 3 4 Row 0: 22.522623 0.000000 0.000000 0.000000 0.000000 1: 0.000000 9.778241 0.000000 0.000000 0.000000 2: 0.000000 0.000000 7.832160 0.000000 0.000000 3: 0.000000 0.000000 0.000000 3.865631 0.000000 4: 0.000000 0.000000 0.000000 0.000000 1.353373 The orthogonal factor V: Col: 0 1 2 3 4 Row 0: -0.641364 -0.121985 0.373256 -0.208728 0.625207 1: -0.228442 -0.071915 -0.384819 0.850936 0.265452 2: -0.475285 -0.629375 -0.351382 -0.194473 -0.465513 3: -0.365665 0.244289 0.597107 0.367581 -0.561213 4: -0.420546 0.723988 -0.482274 -0.243732 -0.083604 The product U * S * V': Col: 0 1 2 3 4 Row 0: 2.000000 1.000000 1.000000 0.000000 9.000000 1: 10.000000 3.000000 4.000000 9.000000 8.000000 2: 8.000000 1.000000 4.000000 4.000000 1.000000 3: 6.000000 0.000000 8.000000 1.000000 -0.000000 4: 4.000000 6.000000 8.000000 0.000000 3.000000 Frobenius Norm of A, A_NORM = 26.096 ABSOLUTE ERROR for A = U*S*V' Frobenius norm of difference A-U*S*V' = 2.65444e-14 RELATIVE ERROR for A = U*S*V': Ratio of DIF_NORM / A_NORM = 1.01718e-15 RANK_ONE_TEST: Compare A to the sum of R rank one matrices. R Absolute Relative Error Error 0 26.096 1 1 13.1807 0.505086 2 8.83841 0.338689 3 4.0957 0.156947 4 1.35337 0.0518614 5 2.65444e-14 1.01718e-15 RANK_ONE_PRINT_TEST: Print the sums of R rank one matrices. Rank R = 0 Col: 0 1 2 3 4 Row 0: 0.000000 0.000000 0.000000 0.000000 0.000000 1: 0.000000 0.000000 0.000000 0.000000 0.000000 2: 0.000000 0.000000 0.000000 0.000000 0.000000 3: 0.000000 0.000000 0.000000 0.000000 0.000000 4: 0.000000 0.000000 0.000000 0.000000 0.000000 Rank R = 1 Col: 0 1 2 3 4 Row 0: 3.701552 1.318426 2.743049 2.110386 2.427128 1: 10.040854 3.576372 7.440814 5.724647 6.583845 2: 5.864444 2.088810 4.345869 3.343528 3.845349 3: 5.141262 1.831225 3.809952 2.931216 3.371155 4: 5.772300 2.055990 4.277586 3.290993 3.784930 Rank R = 2 Col: 0 1 2 3 4 Row 0: 3.022019 0.917813 -0.762964 3.471228 6.460196 1: 9.548351 3.286021 4.899772 6.710940 9.506877 2: 6.091843 2.222871 5.519119 2.888136 2.495726 3: 5.814939 2.228386 7.285750 1.582102 -0.627154 4: 6.233705 2.328007 6.658175 2.366979 1.046471 Rank R = 3 Col: 0 1 2 3 4 Row 0: 1.405763 2.584137 0.758573 0.885663 8.548515 1: 10.551795 2.251493 3.955134 8.316175 8.210355 2: 7.249630 1.029218 4.429182 4.740278 0.999782 3: 5.824490 2.218539 7.276758 1.597380 -0.639494 4: 4.339891 4.280486 8.441004 -0.662603 3.493415 Rank R = 4 Col: 0 1 2 3 4 Row 0: 1.813737 0.920916 1.138686 0.167197 9.024907 1: 10.333463 3.141583 3.751712 8.700669 7.955409 2: 7.326898 0.714213 4.501173 4.604205 1.090008 3: 6.333904 0.141770 7.751384 0.700273 -0.044650 4: 3.925832 5.968510 8.055223 0.066576 3.009918 Rank R = 5 Col: 0 1 2 3 4 Row 0: 2.000000 1.000000 1.000000 0.000000 9.000000 1: 10.000000 3.000000 4.000000 9.000000 8.000000 2: 8.000000 1.000000 4.000000 4.000000 1.000000 3: 6.000000 0.000000 8.000000 1.000000 -0.000000 4: 4.000000 6.000000 8.000000 0.000000 3.000000 Original matrix A: Col: 0 1 2 3 4 Row 0: 2.000000 1.000000 1.000000 0.000000 9.000000 1: 10.000000 3.000000 4.000000 9.000000 8.000000 2: 8.000000 1.000000 4.000000 4.000000 1.000000 3: 6.000000 0.000000 8.000000 1.000000 0.000000 4: 4.000000 6.000000 8.000000 0.000000 3.000000 The pseudoinverse of A: Col: 0 1 2 3 4 Row 0: 0.102837 -0.165669 0.405474 -0.130873 -0.001884 1: -0.042553 -0.030585 0.121011 -0.208777 0.168883 2: -0.056738 0.094637 -0.265736 0.211215 0.006427 3: -0.163121 0.236924 -0.306959 0.095523 -0.040115 4: 0.099291 0.029699 -0.074025 0.028812 -0.019060 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 1.45e-14 A+ * A * A+ = A+ 9.52902e-16 ( A * A+ ) is MxM symmetric; 3.43958e-15 ( A+ * A ) is NxN symmetric; 3.02329e-15 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: 1.000000 -0.000000 0.000000 -0.000000 -0.000000 1: 0.000000 1.000000 -0.000000 0.000000 -0.000000 2: 0.000000 0.000000 1.000000 -0.000000 -0.000000 3: -0.000000 0.000000 0.000000 1.000000 -0.000000 4: -0.000000 0.000000 -0.000000 0.000000 1.000000 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 3 4 Row 0: 1.000000 -0.000000 -0.000000 0.000000 0.000000 1: -0.000000 1.000000 -0.000000 -0.000000 -0.000000 2: -0.000000 0.000000 1.000000 -0.000000 -0.000000 3: -0.000000 -0.000000 0.000000 1.000000 0.000000 4: 0.000000 -0.000000 -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 = 13.0767 Norm of x2 = 13.0767 Norm of residual = 6.195e-14 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 = 16.4924 Norm of x2 = 16.4924 Norm of residual = 2.71521e-13 SVD_DEMO: Normal end of execution. 19 June 2012 08:48:50 PM