08-Jan-2022 10:10:04 svd_test MATLAB/Octave version 9.8.0.1380330 (R2020a) Update 2 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: 5 -13 -13 18 -4 30 -23 3 7 9 36 -1 3 28 7 Orthogonal factor U: -0.2973 -0.0928 -0.4423 -0.6968 -0.4710 0.1368 0.9509 0.1114 -0.0967 -0.2353 0.0019 -0.2284 0.8448 -0.2130 -0.4344 0.7353 -0.1834 -0.2675 0.2466 -0.5416 0.5935 -0.0377 0.0815 -0.6316 0.4907 Diagonal factor S: 48.6940 0 0 0 36.2264 0 0 0 26.9544 0 0 0 0 0 0 Orthogonal factor V: 0.1916 0.5560 -0.8088 0.9531 -0.3020 0.0182 0.2341 0.7744 0.5878 svd_product_test(): Test the value of ||A-U*S*V'|| The product U * X * V': 5.0000 -13.0000 -13.0000 18.0000 -4.0000 30.0000 -23.0000 3.0000 7.0000 9.0000 36.0000 -1.0000 3.0000 28.0000 7.0000 Frobenius Norm of A, A_NORM = 66.407831 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 66.407831 1.000000 1 45.154085 0.679951 2 26.954384 0.405892 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: -2.7739 -13.7979 -3.3894 1.2763 6.3484 1.5595 0.0175 0.0869 0.0213 6.8605 34.1253 8.3827 5.5380 27.5470 6.7668 Rank R = 2: -4.6427 -12.7828 -5.9924 20.4280 -4.0547 28.2355 -4.5818 2.5852 -6.3849 3.1674 36.1314 3.2387 4.7779 27.9600 5.7080 Rank R = 3: 5.0000 -13.0000 -13.0000 18.0000 -4.0000 30.0000 -23.0000 3.0000 7.0000 9.0000 36.0000 -1.0000 3.0000 28.0000 7.0000 Original matrix A: 5 -13 -13 18 -4 30 -23 3 7 9 36 -1 3 28 7 Pseudoinverse of A: 0.0107 0.0118 -0.0288 0.0081 -0.0007 -0.0053 -0.0052 0.0025 0.0157 0.0120 -0.0131 0.0234 0.0136 -0.0062 0.0038 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.2926 -0.1782 -0.3530 -0.0832 -0.2090 -0.1782 0.9353 -0.1228 -0.1036 0.0544 -0.3530 -0.1228 0.7659 -0.1828 0.0786 -0.0832 -0.1036 -0.1828 0.6458 0.4215 -0.2090 0.0544 0.0786 0.4215 0.3604 In some cases, the matrix A+ * A may be interesting (if N <= M, then it MIGHT look like the identity.) A+ * A: 1.0000 0.0000 -0.0000 0.0000 1.0000 0.0000 0.0000 -0.0000 1.0000 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 = 15.165751 Norm of x2 = 15.165751 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 = 15.165751 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 = 209.523268 Norm of x2 = 7.366483 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: 3 -8 14 3 11 -8 -29 -17 -9 11 9 14 -1 0 -9 -11 3 -2 -2 1 -11 -8 3 6 -12 Orthogonal factor U: 0.1029 -0.8516 0.2824 0.0131 0.4292 0.9002 0.1514 -0.2658 0.1774 0.2541 -0.4099 0.2326 -0.3036 0.3539 0.7487 0.0533 0.2539 0.1519 -0.8576 0.4173 0.0906 0.3650 0.8569 0.3281 0.1286 Diagonal factor S: 40.8709 0 0 0 0 0 21.9804 0 0 0 0 0 19.6003 0 0 0 0 0 11.3032 0 0 0 0 0 1.2251 Orthogonal factor V: -0.2977 -0.3859 -0.5539 0.6750 -0.0097 -0.8131 0.1601 -0.2654 -0.4859 -0.0816 -0.3251 -0.6434 0.5634 -0.0431 0.4013 -0.1800 -0.1017 0.4121 0.1881 -0.8672 0.3349 -0.6334 -0.3682 -0.5205 -0.2831 svd_product_test(): Test the value of ||A-U*S*V'|| The product U * X * V': 3.0000 -8.0000 14.0000 3.0000 11.0000 -8.0000 -29.0000 -17.0000 -9.0000 11.0000 9.0000 14.0000 -1.0000 0.0000 -9.0000 -11.0000 3.0000 -2.0000 -2.0000 1.0000 -11.0000 -8.0000 3.0000 6.0000 -12.0000 Frobenius Norm of A, A_NORM = 51.643005 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 51.643005 1.000000 1 31.568501 0.611283 2 22.659082 0.438764 3 11.369382 0.220153 4 1.225085 0.023722 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: -1.2514 -3.4183 -1.3669 -0.7567 1.4080 -10.9516 -29.9157 -11.9622 -6.6220 12.3221 4.9863 13.6206 5.4464 3.0150 -5.6102 -0.6489 -1.7726 -0.7088 -0.3924 0.7301 -1.1028 -3.0124 -1.2046 -0.6668 1.2408 Rank R = 2: 5.9712 -6.4160 10.6764 1.1469 13.2641 -12.2358 -29.3827 -14.1036 -6.9604 10.2140 3.0133 14.4394 2.1566 2.4950 -8.8489 -2.8026 -0.8788 -4.2999 -0.9600 -2.8052 -4.1987 -1.7275 -6.3668 -1.4827 -3.8413 Rank R = 3: 2.9051 -7.8851 13.7953 3.4282 11.2260 -9.3505 -28.0002 -17.0386 -9.1072 12.1319 6.3089 16.0185 -1.1958 0.0429 -6.6582 -4.4515 -1.6689 -2.6226 0.2669 -3.9012 -13.5016 -6.1851 3.0965 5.4390 -10.0251 Rank R = 4: 3.0051 -7.9571 13.7890 3.4560 11.1489 -7.9970 -28.9746 -17.1249 -8.7300 11.0881 9.0089 14.0748 -1.3681 0.7954 -8.7403 -10.9950 3.0417 -2.2051 -1.5567 1.1447 -10.9985 -7.9872 2.9368 6.1366 -11.9554 Rank R = 5: 3.0000 -8.0000 14.0000 3.0000 11.0000 -8.0000 -29.0000 -17.0000 -9.0000 11.0000 9.0000 14.0000 -1.0000 0.0000 -9.0000 -11.0000 3.0000 -2.0000 -2.0000 1.0000 -11.0000 -8.0000 3.0000 6.0000 -12.0000 Original matrix A: 3 -8 14 3 11 -8 -29 -17 -9 11 9 14 -1 0 -9 -11 3 -2 -2 1 -11 -8 3 6 -12 Pseudoinverse of A: 0.0036 0.0069 0.0227 -0.0637 -0.0127 -0.0412 -0.0378 -0.0511 0.0078 -0.0334 0.1728 0.0633 0.2316 0.1365 0.0541 -0.2942 -0.1872 -0.5297 -0.3079 -0.0697 -0.0797 -0.0589 -0.1937 -0.0667 -0.0707 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.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 -0.0000 -0.0000 0.0000 1.0000 -0.0000 0.0000 0 0 0.0000 1.0000 In some cases, the matrix A+ * A may be interesting (if N <= M, then it MIGHT look like the identity.) A+ * A: 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 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 = 20.639767 Norm of x2 = 20.639767 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 = 167.928556 Norm of x2 = 167.928556 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: 11 -11 9 -2 16 2 26 -4 15 24 -8 -2 -8 -12 -7 -18 1 -6 -14 14 7 -11 2 8 -15 7 -14 3 8 1 -1 -9 -7 -2 5 2 -2 7 -19 1 Orthogonal factor U: -0.6921 0.5041 -0.0565 -0.2354 -0.4564 0.6595 0.6995 -0.0614 0.0737 -0.2580 -0.0878 -0.0014 0.7832 0.5605 -0.2544 0.1927 -0.0538 0.5807 -0.7890 -0.0167 0.2029 -0.5037 -0.2059 -0.0505 -0.8125 Diagonal factor S: Columns 1 through 7 43.7319 0 0 0 0 0 0 0 35.1881 0 0 0 0 0 0 0 28.7893 0 0 0 0 0 0 0 20.8980 0 0 0 0 0 0 0 10.1372 0 0 Column 8 0 0 0 0 0 Orthogonal factor V: Columns 1 through 7 -0.0485 0.5789 -0.2789 0.5390 -0.3163 0.3399 -0.2628 0.5696 0.3377 -0.0373 -0.2118 0.1837 0.1508 0.4169 -0.2735 -0.0797 -0.6996 0.0114 -0.2279 -0.3310 0.5127 -0.0041 -0.1022 0.4353 0.2729 -0.3756 0.3207 0.6447 -0.3619 0.0863 0.3518 -0.3179 -0.5453 -0.1581 -0.1115 -0.1537 -0.3112 -0.3075 -0.4145 -0.0713 0.7748 -0.0778 -0.6136 0.5067 0.1341 -0.1802 0.4819 0.1132 0.2412 -0.2592 -0.4160 0.0752 0.5335 0.3721 0.1055 0.0449 Column 8 0.0981 0.5338 0.0695 -0.2570 0.5496 0.0412 -0.1141 0.5638 svd_product_test(): Test the value of ||A-U*S*V'|| The product U * X * V': Columns 1 through 7 11.0000 -11.0000 9.0000 -2.0000 16.0000 2.0000 26.0000 15.0000 24.0000 -8.0000 -2.0000 -8.0000 -12.0000 -7.0000 1.0000 -6.0000 -14.0000 14.0000 7.0000 -11.0000 2.0000 -15.0000 7.0000 -14.0000 3.0000 8.0000 1.0000 -1.0000 -7.0000 -2.0000 5.0000 2.0000 -2.0000 7.0000 -19.0000 Column 8 -4.0000 -18.0000 8.0000 -9.0000 1.0000 Frobenius Norm of A, A_NORM = 67.223508 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 67.223508 1.000000 1 51.054074 0.759468 2 36.990711 0.550265 3 23.226923 0.345518 4 10.137201 0.150798 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 7 1.4671 -17.2420 8.2777 0.1246 10.9554 4.6508 18.5731 -1.3979 16.4294 -7.8875 -0.1188 -10.4391 -4.4316 -17.6977 0.1861 -2.1873 1.0501 0.0158 1.3898 0.5900 2.3562 -0.4085 4.8012 -2.3050 -0.0347 -3.0506 -1.2951 -5.1718 -0.4300 5.0540 -2.4264 -0.0365 -3.2113 -1.3633 -5.4442 Column 8 7.8464 -7.4766 0.9954 -2.1849 -2.3000 Rank R = 2: Columns 1 through 7 11.7342 -11.2527 6.8636 -1.6876 12.4860 -0.8689 27.5610 12.8498 24.7409 -9.8499 -2.6336 -8.3151 -12.0913 -5.2251 0.1579 -2.2038 1.0540 0.0208 1.3856 0.6052 2.3315 -1.5044 4.1619 -2.1541 0.1587 -3.2140 -0.7059 -6.1312 -10.6903 -0.9313 -1.0132 1.7745 -4.7408 4.1527 -14.4261 Column 8 0.4681 -17.7155 1.0157 -1.3974 5.0734 Rank R = 3: Columns 1 through 7 12.1878 -11.1920 8.0015 -2.3956 11.9138 -0.3687 27.3430 13.3425 24.8068 -8.6136 -3.4027 -8.9367 -11.5480 -5.4620 -6.1297 -3.0447 -14.7211 9.8348 9.3171 -6.3280 5.3540 -6.1665 3.5384 -13.8510 7.4356 2.6670 -5.8467 -3.8900 -9.0370 -0.7102 3.1348 -0.8060 -6.8264 5.9758 -15.2209 Column 8 0.3459 -17.8483 2.7104 -0.1408 4.6278 Rank R = 4: Columns 1 through 7 9.5364 -10.1500 7.9455 -3.7379 13.4773 1.6703 28.2296 14.1726 24.4805 -8.5961 -2.9825 -9.4262 -12.1864 -5.7396 0.1842 -5.5262 -14.5878 13.0313 5.5938 -11.1838 3.2428 -15.0536 7.0311 -14.0386 2.9364 7.9077 0.9879 -0.9184 -9.6055 -0.4867 3.1228 -1.0939 -6.4911 6.4131 -15.0308 Column 8 -2.2783 -17.0267 8.9597 -8.9370 4.0650 Rank R = 5: Columns 1 through 7 11.0000 -11.0000 9.0000 -2.0000 16.0000 2.0000 26.0000 15.0000 24.0000 -8.0000 -2.0000 -8.0000 -12.0000 -7.0000 1.0000 -6.0000 -14.0000 14.0000 7.0000 -11.0000 2.0000 -15.0000 7.0000 -14.0000 3.0000 8.0000 1.0000 -1.0000 -7.0000 -2.0000 5.0000 2.0000 -2.0000 7.0000 -19.0000 Column 8 -4.0000 -18.0000 8.0000 -9.0000 1.0000 Original matrix A: 11 -11 9 -2 16 2 26 -4 15 24 -8 -2 -8 -12 -7 -18 1 -6 -14 14 7 -11 2 8 -15 7 -14 3 8 1 -1 -9 -7 -2 5 2 -2 7 -19 1 Pseudoinverse of A: 0.0178 0.0213 0.0149 -0.0266 0.0175 -0.0100 0.0100 -0.0125 0.0089 -0.0161 0.0147 0.0016 -0.0125 -0.0153 0.0231 0.0116 0.0075 0.0286 -0.0008 0.0278 0.0344 0.0083 0.0155 0.0183 0.0390 0.0065 -0.0075 -0.0174 0.0094 0.0127 -0.0030 -0.0124 -0.0121 0.0052 -0.0492 -0.0248 -0.0199 0.0076 -0.0197 -0.0269 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.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.) A+ * A: Columns 1 through 7 0.8058 0.0059 0.2404 0.0856 -0.0294 -0.2878 0.0361 0.0059 0.5185 -0.2009 -0.1799 -0.2230 -0.1064 -0.0567 0.2404 -0.2009 0.6227 -0.2065 -0.0333 0.2935 -0.0783 0.0856 -0.1799 -0.2065 0.4155 0.2639 -0.1878 -0.2211 -0.0294 -0.2230 -0.0333 0.2639 0.6605 0.0912 0.1075 -0.2878 -0.1064 0.2935 -0.1878 0.0912 0.3919 -0.0643 0.0361 -0.0567 -0.0783 -0.2211 0.1075 -0.0643 0.9160 -0.0793 -0.3356 -0.0273 0.0821 -0.2881 -0.1015 0.0415 Column 8 -0.0793 -0.3356 -0.0273 0.0821 -0.2881 -0.1015 0.0415 0.6690 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 = 29.631065 Norm of x2 = 25.992079 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 = 219.544984 Norm of x2 = 219.544984 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.192984 Norm of residual = 2.197066 svd_test(): Normal end of execution. 08-Jan-2022 10:10:05