#1




Misapplication of SVDCMP?
I am attempting to do leastsquares 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 "LeastSquares 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 nonlinear in the parameters, Gander et al suggest the use of GaussNewton 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. GaussNewton is a particularly simple nonlinear 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 GaussJordan (gaussj) or LU decomposition with forward/backward substitution (ludcmp and lubksb) is/are used to solve the normal equations. When I tried GaussJordan on a simple problem with synthetic data, it complained that the pseudoHessian 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? 
Thread Tools  
Display Modes  

