19 June 2012 08:48:58 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 = 8 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 Col: 5 6 7 Row 0: 9.000000 7.000000 6.000000 1: 1.000000 6.000000 2.000000 2: 4.000000 9.000000 8.000000 3: 8.000000 5.000000 4.000000 4: 3.000000 9.000000 2.000000 The orthogonal factor U: Col: 0 1 2 3 4 Row 0: -0.421671 0.487547 0.762151 0.059047 -0.011452 1: -0.487891 -0.791659 0.229657 0.032862 -0.285335 2: -0.490789 -0.017501 -0.209319 -0.530914 0.658131 3: -0.406607 0.350797 -0.437254 -0.284878 -0.662771 4: -0.421845 0.110495 -0.362461 0.795241 0.214595 The diagonal factor S: Col: 0 1 2 3 4 Row 0: 30.710280 0.000000 0.000000 0.000000 0.000000 1: 0.000000 11.245685 0.000000 0.000000 0.000000 2: 0.000000 0.000000 9.734525 0.000000 0.000000 3: 0.000000 0.000000 0.000000 7.538570 0.000000 4: 0.000000 0.000000 0.000000 0.000000 5.179020 Col: 5 6 7 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 The orthogonal factor V: Col: 0 1 2 3 4 Row 0: -0.448566 -0.403244 -0.197960 -0.308931 -0.140848 1: -0.159790 -0.110439 -0.095841 0.583422 0.208195 2: -0.357014 0.083698 -0.570569 0.285166 -0.406579 3: -0.220147 -0.608602 0.081399 -0.280262 -0.115517 4: -0.307861 -0.145066 0.760172 0.351410 -0.209274 5: -0.350517 0.592593 0.171177 -0.192698 -0.466163 6: -0.525094 0.111491 -0.063618 0.207607 0.530699 7: -0.322440 0.251308 0.090783 -0.447875 0.464135 Col: 5 6 7 Row 0: -0.694059 0.000000 0.000000 1: -0.107164 -0.705549 -0.245101 2: 0.300425 0.060309 0.446619 3: 0.620845 -0.127205 -0.284593 4: -0.047513 0.101693 0.355574 5: 0.013793 -0.176028 -0.459405 6: 0.092630 0.530362 -0.305275 7: 0.141653 -0.399713 0.478315 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 Col: 5 6 7 Row 0: 9.000000 7.000000 6.000000 1: 1.000000 6.000000 2.000000 2: 4.000000 9.000000 8.000000 3: 8.000000 5.000000 4.000000 4: 3.000000 9.000000 2.000000 Frobenius Norm of A, A_NORM = 35.327 ABSOLUTE ERROR for A = U*S*V' Frobenius norm of difference A-U*S*V' = 2.23348e-14 RELATIVE ERROR for A = U*S*V': Ratio of DIF_NORM / A_NORM = 6.32231e-16 RANK_ONE_TEST: Compare A to the sum of R rank one matrices. R Absolute Relative Error Error 0 35.327 1 1 17.4608 0.494261 2 13.3571 0.3781 3 9.14616 0.2589 4 5.17902 0.146602 5 2.23348e-14 6.32231e-16 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 Col: 5 6 7 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 3 4 Row 0: 5.808764 2.069225 4.623204 2.850824 3.986684 1: 6.720977 2.394177 5.349236 3.298520 4.612756 2: 6.760901 2.408399 5.381011 3.318114 4.640156 3: 5.601252 1.995304 4.458045 2.748982 3.844264 4: 5.811160 2.070078 4.625111 2.852000 3.988328 Col: 5 6 7 Row 0: 4.539069 6.799775 4.175487 1: 5.251888 7.867617 4.831208 2: 5.283085 7.914352 4.859906 3: 4.376916 6.556860 4.026322 4: 4.540941 6.802580 4.177209 Rank R = 2 Col: 0 1 2 3 4 Row 0: 3.597861 1.463709 5.082101 -0.486015 3.191315 1: 10.310954 3.377390 4.604097 8.716748 5.904245 2: 6.840264 2.430135 5.364538 3.437895 4.668707 3: 4.010477 1.559627 4.788228 0.348079 3.271985 4: 5.310094 1.932848 4.729113 2.095759 3.808071 Col: 5 6 7 Row 0: 7.788137 7.411060 5.553357 1: -0.023821 6.875037 2.593877 2: 5.166455 7.892409 4.810446 3: 6.714665 6.996688 5.017718 4: 5.277291 6.941118 4.489481 Rank R = 3 Col: 0 1 2 3 4 Row 0: 2.129161 0.752649 0.848949 0.117902 8.831164 1: 9.868393 3.163128 3.328531 8.898725 7.603688 2: 7.243632 2.625423 6.527142 3.272033 3.119765 3: 4.853085 1.967570 7.216831 0.001606 0.036345 4: 6.008574 2.271012 6.742303 1.808550 1.125888 Col: 5 6 7 Row 0: 9.058127 6.939064 6.226892 1: 0.358863 6.732812 2.796832 2: 4.817661 8.022039 4.625464 3: 5.986059 7.267477 4.631304 4: 4.673313 7.165588 4.169163 Rank R = 4 Col: 0 1 2 3 4 Row 0: 1.991646 1.012348 0.975885 -0.006851 8.987588 1: 9.791861 3.307661 3.399176 8.829295 7.690744 2: 8.480077 0.290375 5.385812 4.393735 1.713305 3: 5.516538 0.714629 6.604416 0.603489 -0.718334 4: 4.156538 5.768614 8.451868 0.128384 3.232585 Col: 5 6 7 Row 0: 8.972351 7.031477 6.027529 1: 0.311125 6.784243 2.685879 2: 5.588904 7.191128 6.418007 3: 6.399892 6.821626 5.593148 4: 3.518090 8.410186 1.484164 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 Col: 5 6 7 Row 0: 9.000000 7.000000 6.000000 1: 1.000000 6.000000 2.000000 2: 4.000000 9.000000 8.000000 3: 8.000000 5.000000 4.000000 4: 3.000000 9.000000 2.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 Col: 5 6 7 Row 0: 9.000000 7.000000 6.000000 1: 1.000000 6.000000 2.000000 2: 4.000000 9.000000 8.000000 3: 8.000000 5.000000 4.000000 4: 3.000000 9.000000 2.000000 The pseudoinverse of A: Col: 0 1 2 3 4 Row 0: -0.028930 0.037256 0.015911 0.031951 -0.028855 1: -0.005988 -0.000875 -0.009845 -0.045715 0.074850 2: -0.033009 0.009962 -0.053906 0.074221 0.040207 3: -0.018929 0.053404 0.007773 0.005648 -0.040338 4: 0.060670 0.046099 -0.062542 -0.021093 0.002898 5: 0.043428 -0.007267 -0.044669 0.082375 -0.035380 6: 0.007515 -0.029341 0.062404 -0.062472 0.054567 7: 0.017896 -0.037951 0.093333 -0.034441 -0.024496 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 2.60047e-14 A+ * A * A+ = A+ 1.92288e-16 ( A * A+ ) is MxM symmetric; 1.36684e-15 ( A+ * A ) is NxN symmetric; 1.73916e-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: 0.518282 -0.074378 0.208513 0.430903 -0.032977 1: -0.074378 0.430641 0.184212 -0.092971 0.153809 2: 0.208513 0.184212 0.706639 -0.051741 -0.150665 3: 0.430903 -0.092971 -0.051741 0.517377 0.143627 4: -0.032977 0.153809 -0.150665 0.143627 0.860968 5: 0.009573 -0.235319 0.211651 -0.161698 0.181909 6: 0.064291 0.309300 0.076528 -0.076923 0.059015 7: 0.098316 -0.149601 -0.232075 -0.002665 -0.122698 Col: 5 6 7 Row 0: 0.009573 0.064291 0.098316 1: -0.235319 0.309300 -0.149601 2: 0.211651 0.076528 -0.232075 3: -0.161698 -0.076923 -0.002665 4: 0.181909 0.059015 -0.122698 5: 0.757771 -0.048164 0.147426 6: -0.048164 0.616943 0.344889 7: 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 = 2.48763e-13 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 = 2.4923e-13 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.0827285 Norm of residual = 0.857524 SVD_DEMO: Normal end of execution. 19 June 2012 08:48:58 PM