>> svd_demo ( 5, 8 ) 14-Sep-2006 11:55:58 SVD_DEMO (MATLAB 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 = 8 We choose a "random" matrix A, with integral values between 0 and 10. The random number seed used to generate the matrix A is 1067729193 The matrix A: 4.000000 3.000000 7.000000 10.000000 6.000000 5.000000 7.000000 4.000000 4.000000 1.000000 8.000000 2.000000 9.000000 6.000000 8.000000 6.000000 9.000000 3.000000 4.000000 2.000000 6.000000 9.000000 1.000000 7.000000 2.000000 2.000000 7.000000 4.000000 7.000000 6.000000 9.000000 9.000000 5.000000 5.000000 4.000000 3.000000 1.000000 5.000000 6.000000 1.000000 The orthogonal factor U: -0.4767 0.3835 -0.5181 -0.5840 -0.1271 -0.4945 0.1123 0.4514 0.2288 -0.6977 -0.4267 -0.8670 0.0034 -0.2426 0.0855 -0.5078 0.2922 0.3992 0.1003 0.6981 -0.2971 -0.0565 -0.6070 0.7332 0.0492 The diagonal factor S: 33.9639 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 9.3302 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 7.8150 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.9855 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 3.2359 0.0000 0.0000 0.0000 The orthogonal factor V: -0.3011 -0.5914 -0.3165 0.0526 -0.2743 -0.6038 -0.0229 -0.1276 -0.1680 -0.1111 -0.4261 0.3241 0.2533 0.1684 0.2752 0.7104 -0.4046 0.2073 0.0465 0.0816 -0.3232 0.0906 -0.7569 0.3144 -0.2807 0.3563 -0.5753 -0.6553 0.1373 -0.0034 0.0335 -0.1179 -0.4040 0.0106 0.4045 -0.2939 -0.4924 0.1091 0.5298 0.2272 -0.4041 -0.4010 -0.0629 0.1077 0.1181 0.6761 -0.0358 -0.4337 -0.4143 0.5366 -0.0079 0.5620 0.0594 -0.2314 0.2370 -0.3334 -0.3748 -0.1381 0.4664 -0.2057 0.6909 -0.2764 -0.1075 0.1242 The product U * S * V' 4.000000 3.000000 7.000000 10.000000 6.000000 5.000000 7.000000 4.000000 4.000000 1.000000 8.000000 2.000000 9.000000 6.000000 8.000000 6.000000 9.000000 3.000000 4.000000 2.000000 6.000000 9.000000 1.000000 7.000000 2.000000 2.000000 7.000000 4.000000 7.000000 6.000000 9.000000 9.000000 5.000000 5.000000 4.000000 3.000000 1.000000 5.000000 6.000000 1.000000 Frobenius Norm of A, A_NORM = 36.565011 ABSOLUTE ERROR for A = U*S*V': Frobenius norm of difference A-U*S*V' = 0.000000 RELATIVE ERROR for A = U*S*V': Ratio of DIF_NORM / A_NORM = 0.000000 RANK_ONE_TEST: Compare A to the sum of R rank one matrices. R Absolute Relative Error Error 1 13.544504 0.370423 2 9.818386 0.268519 3 5.943618 0.162549 4 3.235944 0.088498 5 0.000000 0.000000 RANK_ONE_PRINT_TEST: Print the sums of R rank one matrices. Rank R = 1 4.875035 2.720119 6.551500 4.544171 6.541873 6.542116 6.708725 6.067748 5.056714 2.821490 6.795656 4.713520 6.785671 6.785922 6.958741 6.293876 4.363400 2.434642 5.863920 4.067260 5.855304 5.855521 6.004644 5.430937 5.193117 2.897599 6.978966 4.840665 6.968712 6.968970 7.146450 6.463651 3.038323 1.695291 4.083166 2.832115 4.077166 4.077317 4.181155 3.781671 Rank R = 2 2.758917 2.322671 7.293214 5.819110 6.579622 5.107418 8.628688 5.573805 4.436953 2.705087 7.012887 5.086919 6.796727 6.365733 7.521053 6.149212 9.147907 3.333264 4.186914 1.184644 5.769955 9.099350 1.663640 6.547735 3.580622 2.594741 7.544157 5.812176 6.997477 5.875721 8.609474 6.087263 3.350092 1.753847 3.973888 2.644278 4.071605 4.288692 3.898286 3.854444 Rank R = 3 4.040333 4.047785 7.104733 8.148596 4.941788 5.362298 8.660874 3.685281 3.320616 1.202213 7.177087 3.057531 8.223566 6.143688 7.493013 7.794444 9.139511 3.321960 4.188149 1.169380 5.780687 9.097680 1.663429 6.560109 2.593405 1.265695 7.689365 4.017514 8.259283 5.679359 8.584678 7.542203 4.851386 3.774973 3.753065 5.373482 2.152733 4.587307 3.935995 1.641868 Rank R = 4 3.887144 3.104209 6.867024 10.056489 5.797437 5.048579 7.024447 4.284213 3.380636 1.571907 7.270221 2.310017 7.888322 6.266604 8.134167 7.559782 9.075867 2.929946 4.089392 1.962025 6.136172 8.967343 0.983566 6.808939 2.619717 1.427767 7.730195 3.689807 8.112313 5.733244 8.865757 7.439328 5.043716 4.959634 4.051509 2.978119 1.078464 4.981183 5.990530 0.889908 Rank R = 5 4.000000 3.000000 7.000000 10.000000 6.000000 5.000000 7.000000 4.000000 4.000000 1.000000 8.000000 2.000000 9.000000 6.000000 8.000000 6.000000 9.000000 3.000000 4.000000 2.000000 6.000000 9.000000 1.000000 7.000000 2.000000 2.000000 7.000000 4.000000 7.000000 6.000000 9.000000 9.000000 5.000000 5.000000 4.000000 3.000000 1.000000 5.000000 6.000000 1.000000 Original matrix A: 4.000000 3.000000 7.000000 10.000000 6.000000 5.000000 7.000000 4.000000 4.000000 1.000000 8.000000 2.000000 9.000000 6.000000 8.000000 6.000000 9.000000 3.000000 4.000000 2.000000 6.000000 9.000000 1.000000 7.000000 2.000000 2.000000 7.000000 4.000000 7.000000 6.000000 9.000000 9.000000 5.000000 5.000000 4.000000 3.000000 1.000000 5.000000 6.000000 1.000000 The pseudoinverse of A: 0.0055 0.0405 0.0488 -0.0883 0.0344 -0.0219 -0.0632 0.0032 0.0384 0.0868 0.0142 0.0845 -0.0267 -0.0532 0.0058 0.1281 -0.0845 0.0057 0.0024 -0.0493 0.0331 0.1220 0.0056 -0.0851 -0.0787 -0.0239 -0.0231 0.0402 0.0179 0.0285 -0.0398 0.0250 -0.0704 0.0467 0.0846 -0.0344 -0.1277 0.0460 0.1700 -0.0519 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.000000 A+ * A * A+ = A+ 0.000000 ( A * A+ ) is MxM symmetric; 0.000000 ( A+ * A ) is NxN symmetric; 0.000000 In some cases, the matrix A * A+ may be interesting (if M <= N, then it MIGHT look like the identity.) 1.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 1.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 1.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 1.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 1.0000 In some cases, the matrix A+ * A+ may be interesting (if N <= M, then it MIGHT look like the identity.) 0.6186 0.1987 0.0775 -0.0163 0.1070 0.3520 -0.1768 -0.1535 0.1987 0.3913 -0.0303 0.0751 -0.3256 0.2041 0.2106 -0.0121 0.0775 -0.0303 0.3200 0.0628 0.3197 0.0481 0.3052 -0.0954 -0.0163 0.0751 0.0628 0.9850 0.0094 -0.0477 -0.0480 0.0173 0.1070 -0.3256 0.3197 0.0094 0.6558 0.0438 -0.0246 0.0589 0.3520 0.2041 0.0481 -0.0477 0.0438 0.3536 0.0203 0.2368 -0.1768 0.2106 0.3052 -0.0480 -0.0246 0.0203 0.7791 0.0029 -0.1535 -0.0121 -0.0954 0.0173 0.0589 0.2368 0.0029 0.8967 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 = 16.093477 Norm of x2 = 13.992648 Norm of residual = 0.000000 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.431677 Norm of x2 = 16.431677 Norm of residual = 0.000000 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' * x1 = b using the pseudoinverse: x2 = A+ * b. Norm of x2 = 0.103199 Norm of residual = 0.111918 SVD_DEMO: Normal end of execution. 14-Sep-2006 11:56:02 >>