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 think there may be some confusion about transposes; see below.

jt

On Thu, 13 Sep 2001, Jim Love wrote:

] 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

isn't it 0.698103, -0.017900, -0.715774
ie, the third column and not the third row?

note also, if you use the columns of V
and not the rows, then the sign ambiguity
doesn't matter.

]
] 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 mean -1,1 ??

let's see z = zo - [a(x-xo)+b(y-yo)]/c
if i use

a = 0.698103
b = -0.017900
c = -0.715774

then i get

 1  1    0.97530415
 1 -1    1.0253199
-1 -1   -0.92530415
-1  1   -0.97531993

which also seems reasonable to me.



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

---------------------------------------------
James Theiler                     jt@lanl.gov
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----




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