Description
If
Sigma is the diagonal
m by
n matrix whose
(i,i) entry is the
i-th element of
S, then
M = U Sigma Vt. This is the singular value decomposition of
M. The entries of
S are (up to roundoff error) the eigenvalues of the Hermitian matrix
M * (conjugate transpose M)
M may also be a
MutableMatrix in which case the returned values
U and
Vt are also
mutable matrices.
If
M is over
CC, then
U and
Vt are unitary matrices over
CC. If
M is over
RR,
U and
Vt are orthogonal over
RR.
i1 : printingPrecision=2;
|
i2 : M = map(RR^3, RR^5, (i,j) -> (i+1)^j * 1.0)
o2 = | 1 1 1 1 1 |
| 1 2 4 8 16 |
| 1 3 9 27 81 |
3 5
o2 : Matrix RR <-- RR
53 53
|
i3 : (S,U,V) = SVD(M)
o3 = ({88 }, | -.016 -.4 -.92 |, | -.014 -.038 -.11 -.32 -.94 |)
{3.9} | -.21 -.89 .4 | | -.28 -.41 -.57 -.59 .29 |
{.8 } | -.98 .2 -.068 | | -.74 -.41 .066 .51 -.15 |
| .084 .33 -.81 .47 -.081 |
| -.61 .74 .077 -.27 .062 |
o3 : Sequence
|
i4 : M' = (transpose U) * M * (transpose V)
o4 = | 88 -8.6e-15 1.1e-14 -1.8e-15 -1.7e-15 |
| 7e-15 3.9 -5.6e-15 -5.8e-16 -2.3e-15 |
| -4.3e-16 1.2e-14 .8 2.6e-16 7.2e-16 |
3 5
o4 : Matrix RR <-- RR
53 53
|
We can clean the small entries from the result above with
clean.
i5 : e = 1e-10;
|
i6 : clean_e M'
o6 = | 88 0 0 0 0 |
| 0 3.9 0 0 0 |
| 0 0 .8 0 0 |
3 5
o6 : Matrix RR <-- RR
53 53
|
i7 : clean_e norm (1 - U * transpose U)
o7 = 0
o7 : RR (of precision 53)
|
Alternatively, if the only issue is display of the matrix, we may set the printing accuracy.
i8 : printingAccuracy = 2
o8 = 2
|
i9 : M'
o9 = | 88 -0 0 -0 -0 |
| 0 3.9 -0 -0 -0 |
| -0 0 .8 0 0 |
3 5
o9 : Matrix RR <-- RR
53 53
|
Now let's try the divide and conquer algorithm and compare answers.
i10 : (S', U', V') = SVD(M, DivideConquer => true)
o10 = ({88 }, | -.02 -.4 -.92 |, | -.01 -.04 -.11 -.32 -.94 |)
{3.9} | -.21 -.89 .4 | | -.28 -.41 -.57 -.59 .29 |
{.8 } | -.98 .2 -.07 | | -.74 -.41 .07 .51 -.15 |
| .08 .33 -.81 .47 -.08 |
| -.61 .74 .08 -.27 .06 |
o10 : Sequence
|
i11 : norm \ ({S', U', V'}-{S, U, V})
o11 = {0, 0, 0}
o11 : List
|
The SVD routine calls on the SVD algorithms in the lapack and eigen libraries.