public inbox for gsl-discuss@sourceware.org
 help / color / mirror / Atom feed
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 -----



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