This is the mail archive of the
gsl-discuss@sources.redhat.com
mailing list for the GSL project.
Re: Problem with Singular Value Decomposition Algorithm
- To: <bjg at network-theory dot co dot uk>
- Subject: Re: Problem with Singular Value Decomposition Algorithm
- From: "Jim Love" <Jim dot Love at asml dot com>
- Date: Thu, 13 Sep 2001 14:52:30 -0400
- Cc: "<"<gsl-discuss at sources dot redhat dot com>
I guess the best way to explain this is via a physical example. If I have a plane that is tilted some angle about the Y-axis and I make 4 measurements of the distance from the xy plane to the surface at 4 points.
1,1 I get 1
1,-1 I get 1
-1, -1 I get -0.9
-1, 1 I get -1
You can then process the data to set up the least square fit to find the "best fit" plane by making a matrix:
x1-mean(x) y1-mean(y) z1-mean(z)
x2-mean(x) y2-mean(y) z2-mean(z)
x3-mean(x) y3-mean(y) z3-mean(z)
x4-mean(x) y4-mean(y) z4-mean(z)
so in this example you would have
1.000000 1.000000 0.975000
1.000000 -1.000000 0.975000
-1.000000 -1.000000 -0.925000
-1.000000 1.000000 -1.025000
find the svd of this matrix.
The solution of the problem is the vector from v the corresponds to the smallest eigenvalue
in our case:
min eigenvalue = 0.035791
the vector from v is -0.698332 -0.000000 0.715774
the equation of a plane is
a(x-xo) + b(y-yo) + c(z-zo) = 0
in this case:
a = -0.698332
b = -0.000000
c = 0.715774
therefore
-0.698332*x + -0.000000*y + 0.715774*z = 0
substituting back our initial points we get:
1,1 = 1.000633
1,-1 = 1.000633
-1, -1 = -0.950633
1,-1 = -0.950633
You can see by inspection that this is correct.
Now the solution that you API provides is:
a = -0.698332
b = -0.000000
c = -0.715774
-0.698332*x + -0.000000*y + -0.715774*z = 0
1,1 = -0.950633
1,-1 = -0.950633
-1,-1 = 1.000633
1,-1 = 1.000633
This is clearly not the "best fit" plane to the data and is clearly wrong.
The sign is not arbitrary for this case or any use of SVD.
James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659
>>> Brian Gough <bjg@network-theory.co.uk> 09/13/01 12:21PM >>>
Columns 2 and 3 have the opposite sign, but this is the arbitrary
minus-sign factor referred to earlier. The results satisfy m =
u*diag(s)*v' and u'*u = I, v'*v = I numerically --- so the
decomposition looks ok. Let me know if there's something I missing
here.
regards
Brian Gough
Jim Love writes:
> This code fixes the order problem of the S vector and the other matrix, but their is still a sign problem. Using this matrix for A:
>
> 1.000000 1.000000 0.975000
> 1.000000 -1.000000 0.975000
> -1.000000 -1.000000 -0.925000
> -1.000000 1.000000 -1.025000
>
> The modified code produces:
>
> s:
> 2.793961 2.000000 0.035791
>
> This is correct!
>
> For V:
>
> -0.715538 -0.025633 0.698103
> 0.018347 -0.999671 -0.017900
> -0.698332 -0.000000 -0.715774
>
> This is NOT correct!
>
> The correct answer for V is:
>
> -0.7155 0.0256 -0.6981
> 0.0183 0.9997 0.0179
> -0.6983 -0.0000 0.7158
>
> U is also wrong: the program outputs:
>
> -0.493230 -0.512652 -0.493875
> -0.506363 0.487019 0.506368
> 0.480733 0.512652 -0.506047
> 0.518861 -0.487019 0.493554
>
> The correct U is:
>
> -0.4932 0.5127 0.4939
> -0.5064 -0.4870 -0.5064
> 0.4807 -0.5127 0.5060
> 0.5189 0.4870 -0.4936
>
> Note last column missing for both solutions for U.
>