Re: SVD routine error?
Posted by Luigi Ballabio-2 on Feb 20, 2004; 8:46am
URL: http://quantlib.414.s1.nabble.com/SVD-routine-error-tp2794p2812.html
On 2004.02.20 13:07, Neil P Firth wrote:
> The SVD class was my doing. One of my early attempts ;)
> If you have some test cases or other code it would be good to get
> them into a test class in the test-suite (on my todo list..), send
> them to me and I'll see what I can do.
Neil,
I exported the SVD class to Python for ease of test. The
interface is a bit different, but the calculations come from your
class. Here's a transcript of a Python session:
>>> from QuantLib import *
>>> A = Matrix([[1,2,3,4],
... [2,0,2,1],
... [0,1,0,0]])
>>> svd = SVD(A)
U should be orthogonal:
>>> U = svd.getU()
>>> print U*transpose(U)
1.169252,0.021675,-0.000000
0.021675,1.002776,-0.000000
-0.000000,-0.000000,1.000000
Close, but no cake. The error is probably too high for numerical
applications.
V should be orthogonal too:
>>> V = svd.getV()
>>> print V*transpose(V)
1.000000,-0.000000,0.000000,0.000000
-0.000000,0.469295,-0.228611,0.443615
0.000000,-0.228611,0.901521,0.191096
0.000000,0.443615,0.191096,1.629184
Ouch. No good.
Now for the final rush:
>>> S = svd.getS()
>>> print U*S*transpose(V)
1.000000,2.027216,3.011724,4.014609
2.000000,-0.013608,1.994138,0.992695
-0.000000,0.469295,-0.228611,-0.284880
Again, close but not enough.
When I used transpose(A) to initialize the SVD, U was bogus, V was good
(V*transpose(V) == I up to 6 decimals) and U*S*transpose(V) was just as
good.
For a square matrix (A minus the rightmost column, or A plus another
row) the results were uniformly good (up to the 6 decimals I see when
printing.)
May I suggest that the test case for the test-suite might bubble
upwards a bit in your todo-list? :)
The test should:
a) create a SVD for a matrix A;
b) retrieve U and verify that U*transpose(U) == I up to some tolerance
(1e-6? 1e-8? For each element? For the L2 norm of the matrix?)
c) retrieve V and do the same
d) retrieve S and while you're at it, verify that it is diagonal and
its elements equal the singular values
e) verify that U*S*transpose(V) == A up to some tolerance
f) repeat for a few matrices (square, rectangular either way, with
elements of the same magnitude, with some elements much larger or much
smaller.)
Thoughts?
Later,
Luigi