Best fit of ellipse with given size and tilt to a set of points

The best fit ellipse has the equation

$ (r - r_0)^T Q (r - r0) = 1 $

where

$ r = [x,y]^T $ is the position vector of each of the data points, and $r_0=[h, k]^T$ is the center of the best fit ellipse, and $Q$ is given and is equal to,

$ Q = \begin{bmatrix} \dfrac{1}{a^2} \cos^2 \theta + \dfrac{1}{b^2} \sin^2 \theta && \cos \theta \sin \theta \bigg( \dfrac{1}{b^2} - \dfrac{1}{a^2} \bigg ) \\ \cos \theta \sin \theta \bigg( \dfrac{1}{b^2} - \dfrac{1}{a^2} \bigg) && \dfrac{1}{a^2} \sin^2 \theta + \dfrac{1}{b^2} \cos^2 \theta \end{bmatrix} $

The error function to be minimized is the sum of the squared errors from the above ellipse model and is given by

$ E = \displaystyle \sum_i \bigg( (r_i - r_0)^T Q (r - r_0 ) - 1 \bigg)^2 = \sum_i \bigg( r_i^T Q r_i - 2 r_i^T Q r_0 + r_0^T Q r_0 - 1 \bigg)^2 $

The gradient of $E$ with respect to $r_0$ is given by

$ \dfrac{\partial E}{\partial r_0} = \displaystyle 2 \sum_i \bigg( r_i^T Q r_i - 2 r_i^T Q r_0 + r_0^T Q r_0 - 1 \bigg) \bigg( - 2 Q (r_i - r_0) \bigg) = 0 $

The above equation is a set of $2$ equations in the two unknowns $r_{0x} = h $ and $r_{0y}= k $

Since these equations are nonlinear, we can use the well-known multivariate Newton-Raphson method to find the solution.

Simply select an initial guess of $r_0$, let's call it $r_0^{(0)} $, for $k = 1, 2, 3, \dots $ you can generate the next guess as follows

1. Set $ y^{(k)} = \displaystyle \sum_i \bigg( r_i^T Q r_i - 2 r_i^T Q r_0^{(k)} + {r_0^{(k)}}^T Q {r_0^{(k)}} - 1 \bigg) \bigg( Q (r_i - r_0^{(k)}) \bigg) $

2. Calculate the Jacobian matrix $J = [J_{ij} ]= \bigg[ \dfrac{\partial y_i}{\partial x_j} \bigg]$, where $x_1 = r_{0x}$, and $x_2 = r_{0y} $

3. Update the estimate of the solution by the Newton-Raphson recurrence

$$ r_{0}^{(k+1)} = r_{0}^{(k)} - J^{-1} y^{(k)} $$

4. Stop when $\| y^{(k)} \| \le \epsilon $

To calculate the Jacobian, define

$ v_i = r_i - {r_0}^{(k)}$

$ f_i = {v_i}^T Q {v_i} - 1 $

Then

$ y^{(k)} = \displaystyle \sum_i f_i Q v_i $

Then by direct differentiation, we have

$ J =\displaystyle \sum_i \bigg( -2 Q v_i v_i^T Q - f_i Q \bigg) $