19 June 2012 07:33:39 AM SVD_DEMO: C++ version Compiled on Jun 19 2012 at 07:32:53. 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 1 1 0 9 1: 10 3 4 9 8 2: 8 1 4 4 1 3: 6 0 8 1 0 4: 4 6 8 0 3 Col: 5 6 7 Row 0: 9 7 6 1: 1 6 2 2: 4 9 8 3: 8 5 4 4: 3 9 2 The orthogonal factor U: Col: 0 1 2 3 4 Row 0: -0.421671 0.487547 0.762151 0.0590472 -0.0114523 1: -0.487891 -0.791659 0.229657 0.0328621 -0.285335 2: -0.490789 -0.0175013 -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.7103 0 0 0 0 1: 0 11.2457 0 0 0 2: 0 0 9.73453 0 0 3: 0 0 0 7.53857 0 4: 0 0 0 0 5.17902 Col: 5 6 7 Row 0: 0 0 0 1: 0 0 0 2: 0 0 0 3: 0 0 0 4: 0 0 0 The orthogonal factor V: Col: 0 1 2 3 4 Row 0: -0.448566 -0.403244 -0.19796 -0.308931 -0.140848 1: -0.15979 -0.110439 -0.0958408 0.583422 0.208195 2: -0.357014 0.0836975 -0.570569 0.285166 -0.406579 3: -0.220147 -0.608602 0.0813995 -0.280262 -0.115517 4: -0.307861 -0.145066 0.760172 0.35141 -0.209274 5: -0.350517 0.592593 0.171177 -0.192698 -0.466163 6: -0.525094 0.111491 -0.0636183 0.207607 0.530699 7: -0.32244 0.251308 0.090783 -0.447875 0.464135 Col: 5 6 7 Row 0: -0.694059 0 0 1: -0.107164 -0.705549 -0.245101 2: 0.300425 0.0603087 0.446619 3: 0.620845 -0.127205 -0.284593 4: -0.0475125 0.101693 0.355574 5: 0.0137932 -0.176028 -0.459405 6: 0.09263 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 1 1 6.63272e-15 9 1: 10 3 4 9 8 2: 8 1 4 4 1 3: 6 1.44329e-15 8 1 -3.9968e-15 4: 4 6 8 3.10862e-15 3 Col: 5 6 7 Row 0: 9 7 6 1: 1 6 2 2: 4 9 8 3: 8 5 4 4: 3 9 2 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 0 0 0 0 1: 0 0 0 0 0 2: 0 0 0 0 0 3: 0 0 0 0 0 4: 0 0 0 0 0 Col: 5 6 7 Row 0: 0 0 0 1: 0 0 0 2: 0 0 0 3: 0 0 0 4: 0 0 0 Rank R = 1 Col: 0 1 2 3 4 Row 0: 5.80876 2.06922 4.6232 2.85082 3.98668 1: 6.72098 2.39418 5.34924 3.29852 4.61276 2: 6.7609 2.4084 5.38101 3.31811 4.64016 3: 5.60125 1.9953 4.45805 2.74898 3.84426 4: 5.81116 2.07008 4.62511 2.852 3.98833 Col: 5 6 7 Row 0: 4.53907 6.79978 4.17549 1: 5.25189 7.86762 4.83121 2: 5.28309 7.91435 4.85991 3: 4.37692 6.55686 4.02632 4: 4.54094 6.80258 4.17721 Rank R = 2 Col: 0 1 2 3 4 Row 0: 3.59786 1.46371 5.0821 -0.486015 3.19131 1: 10.311 3.37739 4.6041 8.71675 5.90425 2: 6.84026 2.43014 5.36454 3.43789 4.66871 3: 4.01048 1.55963 4.78823 0.348079 3.27198 4: 5.31009 1.93285 4.72911 2.09576 3.80807 Col: 5 6 7 Row 0: 7.78814 7.41106 5.55336 1: -0.0238205 6.87504 2.59388 2: 5.16645 7.89241 4.81045 3: 6.71467 6.99669 5.01772 4: 5.27729 6.94112 4.48948 Rank R = 3 Col: 0 1 2 3 4 Row 0: 2.12916 0.752649 0.848949 0.117902 8.83116 1: 9.86839 3.16313 3.32853 8.89873 7.60369 2: 7.24363 2.62542 6.52714 3.27203 3.11977 3: 4.85309 1.96757 7.21683 0.00160564 0.0363447 4: 6.00857 2.27101 6.7423 1.80855 1.12589 Col: 5 6 7 Row 0: 9.05813 6.93906 6.22689 1: 0.358863 6.73281 2.79683 2: 4.81766 8.02204 4.62546 3: 5.98606 7.26748 4.6313 4: 4.67331 7.16559 4.16916 Rank R = 4 Col: 0 1 2 3 4 Row 0: 1.99165 1.01235 0.975885 -0.00685149 8.98759 1: 9.79186 3.30766 3.39918 8.82929 7.69074 2: 8.48008 0.290375 5.38581 4.39374 1.7133 3: 5.51654 0.714629 6.60442 0.603489 -0.718334 4: 4.15654 5.76861 8.45187 0.128384 3.23259 Col: 5 6 7 Row 0: 8.97235 7.03148 6.02753 1: 0.311125 6.78424 2.68588 2: 5.5889 7.19113 6.41801 3: 6.39989 6.82163 5.59315 4: 3.51809 8.41019 1.48416 Rank R = 5 Col: 0 1 2 3 4 Row 0: 2 1 1 6.63272e-15 9 1: 10 3 4 9 8 2: 8 1 4 4 1 3: 6 1.44329e-15 8 1 -3.9968e-15 4: 4 6 8 3.10862e-15 3 Col: 5 6 7 Row 0: 9 7 6 1: 1 6 2 2: 4 9 8 3: 8 5 4 4: 3 9 2 Original matrix A: Col: 0 1 2 3 4 Row 0: 2 1 1 0 9 1: 10 3 4 9 8 2: 8 1 4 4 1 3: 6 0 8 1 0 4: 4 6 8 0 3 Col: 5 6 7 Row 0: 9 7 6 1: 1 6 2 2: 4 9 8 3: 8 5 4 4: 3 9 2 The pseudoinverse of A: Col: 0 1 2 3 4 Row 0: -0.0289305 0.0372563 0.0159113 0.0319513 -0.0288546 1: -0.00598832 -0.000875057 -0.00984532 -0.0457147 0.07485 2: -0.0330086 0.00996223 -0.0539056 0.0742211 0.0402066 3: -0.0189293 0.053404 0.00777347 0.00564769 -0.0403379 4: 0.0606697 0.0460988 -0.0625423 -0.0210926 0.00289757 5: 0.0434277 -0.00726661 -0.0446685 0.0823752 -0.0353797 6: 0.00751516 -0.0293409 0.0624044 -0.0624723 0.0545672 7: 0.0178959 -0.0379505 0.0933326 -0.0344409 -0.0244964 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 5.82867e-16 -5.55112e-16 2.22045e-16 5.82867e-16 1: 6.245e-16 1 5.55112e-17 -1.38778e-17 3.60822e-16 2: 0 2.22045e-16 1 5.55112e-17 5.55112e-17 3: 3.46945e-16 4.16334e-16 -2.22045e-16 1 6.10623e-16 4: 2.15106e-16 3.60822e-16 -3.33067e-16 4.30211e-16 1 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.0743781 0.208513 0.430903 -0.0329765 1: -0.0743781 0.430641 0.184212 -0.0929715 0.153809 2: 0.208513 0.184212 0.706639 -0.0517415 -0.150665 3: 0.430903 -0.0929715 -0.0517415 0.517377 0.143627 4: -0.0329765 0.153809 -0.150665 0.143627 0.860968 5: 0.00957332 -0.235319 0.211651 -0.161698 0.181909 6: 0.0642907 0.3093 0.0765277 -0.0769232 0.059015 7: 0.0983159 -0.149601 -0.232075 -0.00266541 -0.122698 Col: 5 6 7 Row 0: 0.00957332 0.0642907 0.0983159 1: -0.235319 0.3093 -0.149601 2: 0.211651 0.0765277 -0.232075 3: -0.161698 -0.0769232 -0.00266541 4: 0.181909 0.059015 -0.122698 5: 0.757771 -0.0481642 0.147426 6: -0.0481642 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 07:33:39 AM