02-Jun-2023 10:02:50 svd_test MATLAB/Octave version 5.2.0 svd() computes the singular value decomposition. svd_demo(): Test the SVD for a 3 by 3 matrix 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 Choose a "random" matrix A, with integral values. The matrix A: -15 -4 6 7 5 6 12 9 6 1 0 10 9 16 8 Orthogonal factor U: -0.3781073 0.7970853 0.2579877 -0.3800728 -0.1033291 0.3421637 0.1399026 -0.2949883 -0.0076935 -0.8810641 0.5479269 0.0070222 -0.2131509 -0.7543961 0.2918562 0.1356157 0.5183354 -0.6449996 0.4199104 0.3472573 0.6491096 0.2763356 0.6204568 0.3317338 0.0853311 Diagonal factor S: Diagonal Matrix 29.3508 0 0 0 15.3853 0 0 0 7.1988 0 0 0 0 0 0 Orthogonal factor V: 0.70252 -0.51265 -0.49361 0.63168 0.12972 0.76430 0.32779 0.84874 -0.41496 svd_product_test(): Test the value of |A-U*S*V'| The product U * X * V': -1.5000e+01 -4.0000e+00 6.0000e+00 7.0000e+00 5.0000e+00 6.0000e+00 1.2000e+01 9.0000e+00 6.0000e+00 1.0000e+00 2.2204e-15 1.0000e+01 9.0000e+00 1.6000e+01 8.0000e+00 Frobenius Norm of A, A_NORM = 33.911650 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 0 33.911650 1.000000 1 16.986178 0.500895 2 7.198832 0.212282 3 0.000000 0.000000 rank_one_print_test(): Print the sums of R rank one matrices. Rank R = 0: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Rank R = 1: -7.7964 -7.0102 -3.6377 7.0552 6.3438 3.2919 11.2980 10.1588 5.2716 2.7963 2.5144 1.3047 13.3843 12.0347 6.2450 Rank R = 2: -14.0833 -5.4195 6.7707 5.9518 6.6230 5.1188 11.2426 10.1728 5.3633 -1.2920 3.5488 8.0732 11.2048 12.5862 9.8535 Rank R = 3: -1.5000e+01 -4.0000e+00 6.0000e+00 7.0000e+00 5.0000e+00 6.0000e+00 1.2000e+01 9.0000e+00 6.0000e+00 1.0000e+00 2.2204e-15 1.0000e+01 9.0000e+00 1.6000e+01 8.0000e+00 Original matrix A: -15 -4 6 7 5 6 12 9 6 1 0 10 9 16 8 Pseudoinverse of A: -0.053300 0.023755 0.027496 0.030201 -0.036215 0.025973 -0.022775 -0.010779 -0.061191 0.082174 0.024878 0.028543 0.018793 0.067289 -0.013272 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.) A * A+: 0.844868 -0.093964 -0.256568 0.195478 0.134900 -0.093964 0.223667 0.251340 0.309187 0.077734 -0.256568 0.251340 0.345707 0.215430 0.225354 0.195478 0.309187 0.215430 0.703088 -0.168930 0.134900 0.077734 0.225354 -0.168930 0.882671 In some cases, the matrix A+ * A may be interesting (if N <= M, then it MIGHT look like the identity.) A+ * A: 1.00000 0.00000 0.00000 0.00000 1.00000 -0.00000 0.00000 0.00000 1.00000 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.038405 Norm of x2 = 13.038405 Norm of residual = 0.000000 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 = 13.038405 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 = 220.680765 Norm of x2 = 180.814518 Norm of residual = 0.000000 svd_demo(): Test the SVD for a 5 by 5 matrix 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 Choose a "random" matrix A, with integral values. The matrix A: 2 -13 -9 9 -15 -8 -9 3 -1 -9 10 8 3 15 -1 -3 -3 -12 -10 -13 3 11 -20 -4 -15 Orthogonal factor U: 0.44946020 -0.16613917 -0.75612267 0.42392649 -0.13765226 0.15441640 -0.50694667 -0.21446488 -0.80164783 -0.17471757 -0.13981885 0.63581603 -0.51159936 -0.38168002 0.41082249 0.54749438 -0.19978479 0.22781601 0.00084848 0.78002249 0.67441952 0.52079445 0.26201009 -0.17879244 -0.41631169 Diagonal factor S: Diagonal Matrix 35.5578 0 0 0 0 0 24.2814 0 0 0 0 0 22.1044 0 0 0 0 0 7.0236 0 0 0 0 0 2.2639 Orthogonal factor V: -0.038074 0.504222 -0.217601 0.413651 0.725158 -0.072422 0.746949 0.446320 -0.472539 -0.119697 -0.676636 -0.252729 -0.151422 -0.540985 0.403359 -0.179403 0.348564 -0.795807 -0.057169 -0.457975 -0.709425 0.049589 0.311787 0.556478 -0.295600 svd_product_test(): Test the value of |A-U*S*V'| The product U * X * V': 2.00000 -13.00000 -9.00000 9.00000 -15.00000 -8.00000 -9.00000 3.00000 -1.00000 -9.00000 10.00000 8.00000 3.00000 15.00000 -1.00000 -3.00000 -3.00000 -12.00000 -10.00000 -13.00000 3.00000 11.00000 -20.00000 -4.00000 -15.00000 Frobenius Norm of A, A_NORM = 48.959167 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 0 48.959167 1.000000 1 33.654800 0.687405 2 23.303656 0.475981 3 7.379422 0.150726 4 2.263864 0.046240 5 0.000000 0.000000 rank_one_print_test(): Print the sums of R rank one matrices. Rank R = 0: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Rank R = 1: -0.60849 -1.15743 -10.81387 -2.86719 -11.33788 -0.20905 -0.39765 -3.71521 -0.98505 -3.89524 0.18929 0.36005 3.36400 0.89193 3.52701 -0.74121 -1.40988 -13.17254 -3.49256 -13.81085 -0.91305 -1.73673 -16.22631 -4.30224 -17.01260 Rank R = 2: -2.64256 -4.17069 -9.79433 -4.27332 -11.53793 -6.41570 -9.59211 -0.60427 -5.27565 -4.50564 7.97371 11.89182 -0.53776 6.27323 4.29258 -3.18721 -5.03337 -11.94653 -5.18346 -14.05140 5.46314 7.70889 -19.42223 0.10556 -16.38553 Rank R = 3: 0.99434 -11.63032 -7.26353 9.02750 -16.74902 -5.38414 -11.70795 0.11356 -1.50303 -5.98370 10.43447 6.84456 1.17461 15.27268 0.76671 -4.28300 -2.78581 -12.70905 -9.19094 -12.48133 4.20289 10.29379 -20.29919 -4.50342 -14.57979 Rank R = 4: 2.22598 -13.03730 -8.87430 8.85728 -15.09212 -7.71317 -9.04734 3.15954 -1.18115 -9.11692 9.32557 8.11132 2.62486 15.42594 -0.72508 -4.28053 -2.78863 -12.71228 -9.19128 -12.47801 3.68344 10.88719 -19.61985 -4.43163 -15.27859 Rank R = 5: 2.00000 -13.00000 -9.00000 9.00000 -15.00000 -8.00000 -9.00000 3.00000 -1.00000 -9.00000 10.00000 8.00000 3.00000 15.00000 -1.00000 -3.00000 -3.00000 -12.00000 -10.00000 -13.00000 3.00000 11.00000 -20.00000 -4.00000 -15.00000 Original matrix A: 2 -13 -9 9 -15 -8 -9 3 -1 -9 10 8 3 15 -1 -3 -3 -12 -10 -13 3 11 -20 -4 -15 Pseudoinverse of A: -0.015613 -0.111759 0.127504 0.242928 -0.136369 -0.042537 0.042932 0.013471 -0.043960 0.053978 -0.058822 0.034423 0.102143 0.129014 -0.080453 0.046966 0.041535 -0.051751 -0.171636 0.080315 0.031589 -0.047842 -0.087011 -0.109901 0.031497 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.) A * A+: 1.00000 -0.00000 -0.00000 -0.00000 0.00000 0.00000 1.00000 -0.00000 -0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 -0.00000 -0.00000 1.00000 0.00000 0.00000 -0.00000 -0.00000 0.00000 1.00000 In some cases, the matrix A+ * A may be interesting (if N <= M, then it MIGHT look like the identity.) A+ * A: 1.00000 -0.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 -0.00000 0.00000 -0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 0.00000 -0.00000 1.00000 -0.00000 0.00000 0.00000 0.00000 -0.00000 1.00000 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 = 26.343880 Norm of x2 = 26.343880 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 = 236.431808 Norm of x2 = 236.431808 Norm of residual = 0.000000 svd_demo(): Test the SVD for a 8 by 8 matrix 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 Choose a "random" matrix A, with integral values. The matrix A: -1 -8 4 16 -12 -9 6 9 -4 5 -11 4 3 11 21 -17 6 4 2 -2 6 -1 -10 5 19 5 -2 7 4 14 -13 2 -5 -5 -3 -9 -5 -0 -2 -11 Orthogonal factor U: -0.052339 -0.667110 -0.619020 0.283269 -0.297989 0.770596 0.396290 -0.407213 -0.099312 -0.271018 -0.355197 0.160441 0.130167 -0.306544 -0.858593 -0.476580 0.609251 -0.471014 0.414013 0.091784 0.223939 0.031480 0.460646 0.802820 -0.303555 Diagonal factor S: Diagonal Matrix 36.3564 0 0 0 0 0 0 0 0 31.6016 0 0 0 0 0 0 0 0 23.8346 0 0 0 0 0 0 0 0 12.5518 0 0 0 0 0 0 0 0 2.6304 0 0 0 Orthogonal factor V: Columns 1 through 6: -0.4218226 0.3627335 -0.3450291 0.1694477 -0.1930617 -0.4660801 -0.0179249 0.3433037 -0.0512507 -0.4726764 -0.1630386 -0.4541064 -0.2507116 -0.2537746 0.0765142 -0.1293896 0.3038197 0.0681619 -0.0659076 -0.1717650 -0.8070802 0.0335302 -0.2890216 0.4030432 -0.0609888 0.3935383 0.1174895 -0.6289537 -0.1915423 0.5210627 0.0723584 0.5927624 -0.2363173 0.1960584 0.7011457 0.1988027 0.6922614 -0.1667070 -0.3509762 -0.3432431 0.1978911 -0.2923944 -0.5161037 -0.3501874 -0.1681114 -0.4220913 0.4391481 -0.1137697 Columns 7 and 8: 0.5262409 -0.1062991 -0.6283251 -0.1678837 0.0627737 -0.8655941 -0.2424859 -0.0941991 0.3518388 0.0060061 -0.1339178 0.0071020 0.3518513 -0.0593979 0.0037797 0.4458412 svd_product_test(): Test the value of |A-U*S*V'| The product U * X * V': Columns 1 through 6: -1.0000e+00 -8.0000e+00 4.0000e+00 1.6000e+01 -1.2000e+01 -9.0000e+00 -4.0000e+00 5.0000e+00 -1.1000e+01 4.0000e+00 3.0000e+00 1.1000e+01 6.0000e+00 4.0000e+00 2.0000e+00 -2.0000e+00 6.0000e+00 -1.0000e+00 1.9000e+01 5.0000e+00 -2.0000e+00 7.0000e+00 4.0000e+00 1.4000e+01 -5.0000e+00 -5.0000e+00 -3.0000e+00 -9.0000e+00 -5.0000e+00 5.5511e-16 Columns 7 and 8: 6.0000e+00 9.0000e+00 2.1000e+01 -1.7000e+01 -1.0000e+01 5.0000e+00 -1.3000e+01 2.0000e+00 -2.0000e+00 -1.1000e+01 Frobenius Norm of A, A_NORM = 55.253959 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 0 55.253959 1.000000 1 41.607856 0.753029 2 27.065705 0.489842 3 12.824409 0.232099 4 2.630366 0.047605 5 0.000000 0.000000 rank_one_print_test(): Print the sums of R rank one matrices. Rank R = 0: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 Rank R = 1: Columns 1 through 6: 0.802674 0.034109 0.477072 0.125414 0.116054 -0.137689 -11.817821 -0.502186 -7.023960 -1.846473 -1.708669 2.027199 5.447277 0.231476 3.237606 0.851109 0.787590 -0.934412 7.308805 0.310580 4.344012 1.141963 1.056737 -1.253734 -3.434312 -0.145937 -2.041194 -0.536593 -0.496547 0.589113 Columns 7 and 8: -1.317284 0.982079 19.394459 -14.459209 -8.939633 6.664792 -11.994625 8.942388 5.636117 -4.201912 Rank R = 2: Columns 1 through 6: -6.844385 -7.203336 5.827084 3.746522 -8.180424 -12.634160 -7.275169 3.797139 -10.202076 -3.997553 3.219763 9.450591 7.286408 1.972094 1.950919 -0.019774 2.782907 2.071011 14.292621 6.920308 -0.541984 -2.165080 8.633646 10.158895 -3.073463 0.195583 -2.293651 -0.707466 -0.105053 1.178797 Columns 7 and 8: 2.197193 8.364642 17.306722 -18.844740 -9.784871 4.889272 -15.204284 2.200128 5.470276 -4.550280 Rank R = 3: Columns 1 through 7: -1.75380 -6.44718 4.69819 15.65424 -9.91387 -9.14752 7.37552 -3.92641 4.29456 -10.94470 3.83576 2.07944 11.74422 20.71321 6.21597 1.81309 2.18830 -2.52372 3.14741 1.33784 -10.87376 18.16606 7.49567 -1.40096 6.89553 7.31466 12.81189 -11.26408 -6.86164 -0.36711 -1.45358 -9.56865 1.18490 -1.41580 1.61680 Column 8: 10.84497 -17.21309 4.36771 4.08742 -6.39603 Rank R = 4: Columns 1 through 7: -1.15133 -8.12779 4.23814 15.77346 -12.15013 -8.45043 6.15511 -4.13763 4.88377 -10.78341 3.79396 2.86345 11.49983 21.14107 5.56399 3.63179 2.68615 -2.65273 5.56742 0.58348 -9.55308 19.04661 5.03936 -2.07335 7.06978 4.04624 13.83073 -13.04778 -5.15415 -5.13018 -2.75741 -9.23077 -5.15294 0.55984 -1.84199 Column 8: 9.34421 -16.68694 5.99178 1.89398 -10.64936 Rank R = 5: Columns 1 through 6: -1.0000e+00 -8.0000e+00 4.0000e+00 1.6000e+01 -1.2000e+01 -9.0000e+00 -4.0000e+00 5.0000e+00 -1.1000e+01 4.0000e+00 3.0000e+00 1.1000e+01 6.0000e+00 4.0000e+00 2.0000e+00 -2.0000e+00 6.0000e+00 -1.0000e+00 1.9000e+01 5.0000e+00 -2.0000e+00 7.0000e+00 4.0000e+00 1.4000e+01 -5.0000e+00 -5.0000e+00 -3.0000e+00 -9.0000e+00 -5.0000e+00 5.5511e-16 Columns 7 and 8: 6.0000e+00 9.0000e+00 2.1000e+01 -1.7000e+01 -1.0000e+01 5.0000e+00 -1.3000e+01 2.0000e+00 -2.0000e+00 -1.1000e+01 Original matrix A: -1 -8 4 16 -12 -9 6 9 -4 5 -11 4 3 11 21 -17 6 4 2 -2 6 -1 -10 5 19 5 -2 7 4 14 -13 2 -5 -5 -3 -9 -5 -0 -2 -11 Pseudoinverse of A: 0.02760660 0.02005406 0.06295851 0.01819351 0.02421293 0.00191266 0.02533927 0.06640044 -0.01341369 -0.01217633 -0.03360833 -0.04008375 -0.09443263 0.00321549 -0.04365622 0.05818138 0.03975191 0.08888649 0.00452269 0.01932358 -0.00376598 0.02634689 0.08111845 -0.02136464 -0.01583650 -0.08148673 -0.06078892 -0.23264123 0.04608210 -0.07190622 -0.01852709 0.00090499 -0.06573838 -0.00976909 -0.04747685 -0.04677452 -0.05436609 -0.13069007 0.00473742 -0.08445367 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.) A * A+: 1.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 1.00000 0.00000 -0.00000 -0.00000 0.00000 0.00000 1.00000 0.00000 -0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 -0.00000 -0.00000 -0.00000 0.00000 1.00000 In some cases, the matrix A+ * A may be interesting (if N <= M, then it MIGHT look like the identity.) A+ * A: Columns 1 through 6: 0.4945403 0.1011545 -0.0932771 0.3054431 0.0583434 0.1638860 0.1011545 0.3708100 -0.0749241 0.0148500 0.4586954 0.0073260 -0.0932771 -0.0749241 0.2421603 -0.0937886 -0.0524040 0.0010032 0.3054431 0.0148500 -0.0937886 0.7698833 -0.1241290 -0.1119303 0.0583434 0.4586954 -0.0524040 -0.1241290 0.6046671 -0.0565138 0.1638860 0.0073260 0.0010032 -0.1119303 -0.0565138 0.9424931 -0.3277517 0.0783269 -0.0535713 0.1975713 0.0289176 0.1056698 -0.0076223 0.0255608 0.3934350 0.0887684 0.0552735 0.0199575 Columns 7 and 8: -0.3277517 -0.0076223 0.0783269 0.0255608 -0.0535713 0.3934350 0.1975713 0.0887684 0.0289176 0.0552735 0.1056698 0.0199575 0.7871781 -0.0081135 -0.0081135 0.7882678 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 = 17.691806 Norm of x2 = 14.985480 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 = 131.148770 Norm of x2 = 131.148770 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.161820 Norm of residual = 1.194559 svd_test(): Normal end of execution. 02-Jun-2023 10:02:50