From: "Jim Love" <Jim.Love@asml.com>
To: <jt@lanl.gov>
Cc: "<Brian Gough" <bjg@network-theory.co.uk>,
"<"<gsl-discuss@sources.redhat.com>
Subject: Re: Problem with Singular Value Decomposition Algorithm
Date: Wed, 19 Dec 2001 13:20:00 -0000 [thread overview]
Message-ID: <sba0e31d.048@wiltonhub.svg.com> (raw)
You have made a mistake, the equation for a z at point x,y of the "best fit" plane is:
z = (a(x-mean(x))+b(y-mean(y))-c(mean(z)))/-c
so using your values of a,b,c i get
1,1 = -0.975
1,-1 = -1.025
no need to go and calc the others - this is clearly wrong.
James A. Love
X4477
Pager 1-800-286-1188 Pin# 400659
>>> James Theiler <jt@lanl.gov> 09/13/01 03:28PM >>>
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 -----
next reply other threads:[~2001-12-19 13:20 UTC|newest]
Thread overview: 15+ messages / expand[flat|nested] mbox.gz Atom feed top
2001-12-19 13:20 Jim Love [this message]
-- strict thread matches above, loose matches on Subject: below --
2001-12-19 13:20 Jim Love
2001-12-19 13:20 ` James Theiler
2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 ` Peter Schmitteckert
2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 Jim Love
2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 ` Brian Gough
[not found] <sba09ba1.005@wiltonhub.svg.com>
2001-12-19 13:20 ` Brian Gough
2001-12-19 13:20 Jim Love
2001-12-19 13:20 ` James Theiler
2001-12-19 13:20 Jim Love
2001-12-19 13:20 Jim Love
2001-12-19 13:20 ` Brian Gough
Reply instructions:
You may reply publicly to this message via plain-text email
using any one of the following methods:
* Save the following mbox file, import it into your mail client,
and reply-to-all from there: mbox
Avoid top-posting and favor interleaved quoting:
https://en.wikipedia.org/wiki/Posting_style#Interleaved_style
* Reply using the --to, --cc, and --in-reply-to
switches of git-send-email(1):
git send-email \
--in-reply-to=sba0e31d.048@wiltonhub.svg.com \
--to=jim.love@asml.com \
--cc=bjg@network-theory.co.uk \
--cc=gsl-discuss@sources.redhat.com \
--cc=jt@lanl.gov \
/path/to/YOUR_REPLY
https://kernel.org/pub/software/scm/git/docs/git-send-email.html
* If your mail client supports setting the In-Reply-To header
via mailto: links, try the mailto: link
Be sure your reply has a Subject: header at the top and a blank line
before the message body.
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).