Numerical Recipes Forum  

Go Back   Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 2, 11, and 18

Thread Tools Display Modes
Old 09-22-2008, 08:11 PM
Traveler Traveler is offline
Registered User
Join Date: Sep 2008
Posts: 1
Misapplication of SVDCMP?

I am attempting to do least-squares data reduction on data which describes a circular orbit in 3 dimensions. I'm trying to determine the center of orbit and orbit radius. Gander et al, in their paper "Least-Squares Fitting of Circles and Ellipses" (available on CiteSeer) in section 3, describe a particularly simple residual equation for this problem in 2 dimensions. It is easily extended to 3 dimensions.

Because the Jacobian is non-linear in the parameters, Gander et al suggest the use of Gauss-Newton to solve this problem. The highly recommended NR routine "mrqmin" is not readily applicable since there is no easily identifiable independent data vector, "x". With only the 3D samples to work with, it isn't possible to provide a predictor function "funcs" for the "mrqmin" routine.

Gauss-Newton is a particularly simple non-linear least squares solver. The Jacobian of the residual equation(s) is used to calculate both the gradient and an approximation of the Hessian for the problem. In turn, those are used to write "normal equations". Instead of the normal equations generating a solution to the problem, they generate an update vector, which is used to iteratively improve the parameter vector.

Typically an algorithm like Gauss-Jordan (gaussj) or LU decomposition with forward/backward substitution (ludcmp and lubksb) is/are used to solve the normal equations. When I tried Gauss-Jordan on a simple problem with synthetic data, it complained that the pseudo-Hessian matrix was singular.

When I reverted from 3D to 2D, the singularity problem went away. The issue is that a circle is a planar object in an essential way. The singularity appears to be an essential part of the problem, not an attribute of poorly selected basis functions as suggested in NR.

When I substituted "svdfit" for "gaussj", I was able to eliminate the singularity. As you might guess, the null value on the diagonal of the "w" matrix is in the same column as the column vector in the "v" matrix, which is normal to the plane of the circle.

When I tried the modified algorithm on synthetic data from a randomly oriented circle, the solution vector headed towards the correct center and radius of the circle, but stopped far short.

Next, I oriented the normal vector along the Z axis, and I presented the same data to both the 2D version of the algorithm (using gaussj) and the 3D version of the algorithm (using svdfit) (though in the 3D case, all Z coordinates were set to 0). Both algorithms converged to the correct answer in the same number of iterations, with all of the same intermediate results (except for the extra zero Z coordinates and null Z vectors).

I apologize for being so long winded... Am I applying svdcmp incorrectly? What would cause it to fail to produce useable results when the null space isn't conveniently aligned with one of the axes?
Reply With Quote

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is Off
HTML code is Off

Forum Jump

All times are GMT -5. The time now is 10:22 AM.

Powered by vBulletin® Version 3.8.6
Copyright ©2000 - 2014, Jelsoft Enterprises Ltd.