Least-squares fit of a nonlinear (polar) system I want to determine the six unknown coefficients (uppercase letters) of the model 
$$x=X_c+(Au+B)\cos(Cv+D),\\y=Y_c+(Au+B)\sin(Cv+D)$$
given a set of data $(x_k,y_k,u_k,v_k)$, by least-squares. As you see, this model is a (modified) polar-to-Cartesian transform. The number of points is moderate, say between $10$ and $100$.
I know that such a fitting can be done by Levenberg–Marquardt, but I wonder if a simpler ad-hoc approach is possible. For LM, I also need a good starting approximation anyway.
Do you see a way to make the problem more tractable, possibly more "linear", or suggest a way to get an approximate solution ? For instance, first finding the polar origin without looking for the other coefficients.

Update:
I am now considering the following approach: one can choose points vertically aligned in the $(u,v)$ plane and relate them to nearby known data points by some form of interpolation (f.i. linear between three points). Then the corresponding points are estimated in the $(x,y)$ plane by applying the same interpolation.
In this plane, the reconstructed points will be approximately aligned on radial lines, from which the polar origin can be estimated. (Similarly, points reconstructed from horizontals will be aligned on polar circles.)
Knowing $(X_c,Y_c)$, by applying the inverse polar transform you recompute quadruples $(x_k,y_k,\rho_k,\theta_k)$ which allow you to estimate $A,B,C,D$ by two linear fits.
 A: $$\begin{cases} x-X_c=(Au+B)\cos(Cv+D)\\y-Y_c=(Au+B)\sin(Cv+D)\end{cases} \quad\to\quad (x-X_c)^2+(y-Y_c)^2=(Au+B)^2$$
FIRST PART :
In 3D. Cartesian system $(x,y,u)$, this equation is the equation of a cone with vertical axe and center $(X_c\:,Y_c\:,-\frac{B}{A})$.
So, the problem can be interpreted as the fitting of a cone to a given set of points.
A cone is a particular quadratic surface. The problem is a particular case of the general problem of regression of quadratic surface. This problem is treated page 17 in the paper : https://fr.scribd.com/doc/14819165/Regressions-coniques-quadriques-circulaire-spherique (in French). A rough translation (adapted to the present case of cone) is given below.
The equation can be written as :
$$x^2+y^2=A^2u^2+2X_cx+2Y_cy+2ABu-X_c^2-Y_c^2+B^2$$
Let $\quad \begin{cases}R=x^2+y^2 \\ S=A^2 \\ T=2X_c \\ U=2Y_c \\ V=2AB \\ W=-X_c^2-Y_c^2+B^2 \end{cases} \to\quad R(x,y,u)=Su^2+Tx+Uy+Vu+W$
Given the data $(x_k,y_k,u_k)$ one computes $R_k=X_k^2+y_k^2$ from $k=1$ to $n$. 
Then, the linear regression for $S,T,U,V,W$ $$R_k \simeq Su_k^2+Tx_k+Uy_k+Vu_k+W\qquad (n\text{ points})$$
leads to approximates of $X_c\:,Y_c\:,A,B$.
These approximates are good only if the scatter of the data is small, due to the not-matching number of parameters.
SECOND PART :
Using the approximates $X_c\:,Y_c\:,A,B$ computed above, the system of equation can be rewritten as :
$$\begin{cases} \cos^{-1}\left(\frac{x-X_c}{Au+B}\right)=Cv+D\\ \sin^{-1}\left(\frac{y-Y_c}{Au+B}\right)=Cv+D\end{cases}$$
Let $Y_k= \cos^{-1}\left(\frac{x-X_c}{Au+B}\right)$  and $Y_{n+k}=\sin^{-1}\left(\frac{y-Y_c}{Au+B}\right)$ from $k=1$ to $n$ which gives $2n$ values of $Y_k$.
Then, the linear regression for $C,D$ 
$$Y_k \simeq Cv_k+D\qquad (2n\text{ points})$$
leads to approximates of $C$ and $D$.
COMMENTS :
The above approximates of $X_c\:,Y_c\:,A,B,C,D$ are obtained according to a criteria of fitting which is not exactly the least mean squares of distances between the cone and the points. If you definitively want the mean least squares, a non linear regression has to be used, for example with the Levenberg–Marquardt method. It would requires an iterative process starting from initial guess of the six parameters. This is a difficult guess in general. The above direct method provides some values of the parameters which can be used as initial guess.
