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