Tue Oct 19 17:12:08 2021 svd_test() Python version: 3.6.9 svd() carries out the singular value decomposition. 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 = #d 3 Matrix column order N = #d 3 We choose a "random" matrix A, with integral values between 0 and 10. The matrix A: Col: 0 1 2 Row 0 : 4 0 1 1 : 9 1 0 2 : 1 7 8 The orthogonal factor U: Col: 0 1 2 Row 0 : -0.252401 -0.310811 -0.916346 1 : -0.495446 -0.77194 0.398298 2 : -0.83116 0.554531 0.0408489 The diagonal factor S: Col: 0 1 2 Row 0 : 11.2775 0 0 1 : 0 9.21961 0 2 : 0 0 0.904073 The orthogonal factor V: Col: 0 1 2 Row 0 : -0.558616 -0.828254 -0.0440835 1 : -0.559838 0.3373 0.756842 2 : -0.611988 0.447464 -0.65211 The product U * S * V Col: 0 1 2 Row 0 : 4 2.68805e-16 1 1 : 9 1 -9.67486e-16 2 : 1 7 8 Frobenius Norm of A, A_NORM = 14.594520 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 9.263826 0.634747 2 0.904073 0.061946 3 0.000000 0.000000 rank_one_print_test: Print the sums of R rank one matrices. Rank R = 1 Col: 0 1 2 Row 0 : 1.59007 1.59355 1.74199 1 : 3.1212 3.12803 3.41941 2 : 5.23612 5.24758 5.7364 Rank R = 2 Col: 0 1 2 Row 0 : 3.96348 0.627001 0.459764 1 : 9.01587 0.727469 0.234818 2 : 1.00163 6.97205 8.02408 Rank R = 3 Col: 0 1 2 Row 0 : 4 2.68805e-16 1 1 : 9 1 -9.67486e-16 2 : 1 7 8 Original matrix A: Col: 0 1 2 Row 0 : 4 0 1 1 : 9 1 0 2 : 1 7 8 The pseudoinverse of A: Col: 0 1 2 Row 0 : 0.0851064 0.0744681 -0.0106383 1 : -0.765957 0.329787 0.0957447 2 : 0.659574 -0.297872 0.0425532 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+: Col: 0 1 2 Row 0 : 1 0 -4.23273e-16 1 :-3.33067e-16 1 1.38778e-16 2 :-8.88178e-16 4.44089e-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 Row 0 : 1 -1.21431e-16 -1.52656e-16 1 :-2.22045e-16 1 2.10942e-15 2 :-2.08167e-17 -2.47719e-15 1 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 = 8.246211 Norm of x2 = 8.246211 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 = 9.539392 Norm of x2 = 9.539392 Norm of residual = 0.000000 svd_test: Normal end of execution. Tue Oct 19 17:12:09 2021