This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Re: Problem with Singular Value Decomposition Algorithm


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.
 > 


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]