This is a partial answer. I'm going to derive the equations satisfied by the center of desired ellipses. As of this moment,, solving the equations by anything non-numerical is out of reach. The main result is
The center lies on the intersection of two algebraic curves, one cubic and another quartic.
To avoid conflict with other uses of the variable $a,b$. I will use $\alpha,\beta$ to denote the semi-major/semi-minor axis for the desired ellipse.
main result - center lies on intersection of a cubic and a quartic curve.
First, let us determine the condition for the origin to be such a center. If $\theta$ is the angle between the semi-major axis with the $x$-axis.The equation of the ellipse will be
$$\begin{align}
& {\small \frac{(x\cos\theta + y\sin\theta)^2}{\alpha^2} + \frac{(-x\sin\theta+y\cos\theta)^2}{\beta^2} = 1}\\
\iff & {\small \left(\frac{\cos^2\theta}{\alpha^2} + \frac{\sin^2\theta}{\beta^2}\right)x^2
+ \left(\frac{\sin^2\theta}{\alpha^2} + \frac{\cos^2\theta}{\beta^2}\right)y^2 +
\left(\frac{1}{\alpha^2}-\frac{1}{\beta^2}\right)\sin(2\theta)xy = 1}\\
\iff & {\small \frac12\left(\frac{1}{\alpha^2} + \frac{1}{\beta^2}\right)(x^2 + y^2)
+ \frac12\left(\frac{1}{\alpha^2}-\frac{1}{\beta^2}\right)\left((x^2-y^2)\cos(2\theta) + 2xy\sin(2\theta)\right) = 1 }
\end{align}
$$
Let $\epsilon = \frac{\sqrt{\alpha^2-\beta^2}}{\alpha}$ be the eccentricity of the ellipse and define $\sigma$ and $\lambda$ by
$$\frac{1}{\sigma^2} = \frac12\left(\frac1{\alpha^2} + \frac{1}{\beta^2}\right)
\quad\text{ and }\quad
\lambda = \frac{\alpha^2-\beta^2}{\alpha^2+\beta^2} = \frac{\epsilon^2}{2-\epsilon^2}$$
$\sigma$ and $\lambda$ will be an alternative measure of the size and eccentricity of the ellipse. In particular, for the conic above to be an ellipse, we need $0 \le \epsilon < 1 \iff 0 \le \lambda < 1$. In terms of them, above equation becomes
$$\frac{1}{\sigma^2}(x^2+y^2) - \frac{\lambda}{\sigma^2}\left[(x^2-y^2)\cos(2\theta)+2xy\sin(2\theta)\right] = 1$$
For any $p = (x,y), q = (x', y') \in \mathbb{R}^2$, define
- $\tilde{p} = (-y,x)$, the image of $p$ after a $90^\circ$ counterclockwise rotation.
- $p \cdot q = xx' + y y'\;$ be ordinary 2-d dot product.
- $p \times q \stackrel{def}{=} \tilde{p}\cdot q = x y' - y x'\;$ be the "cross product" for 2-d numbers.
- $U(p) = \begin{bmatrix} x^2 + y^2\\ x^2 - y^2\\ 2xy\end{bmatrix}$
and $V(p) = U(\tilde{p}) = \begin{bmatrix} x^2 + y^2\\ -x^2 + y^2\\ -2xy\end{bmatrix}$
In terms of them, the equation for 3 points $A,B,C$ to lie on the ellipse centered at origin is
$$\Lambda \cdot U(A) = \Lambda \cdot U(B) = \Lambda U(C) = 1\quad\text{ where }\quad
\Lambda = \frac{1}{\sigma^2}\begin{bmatrix}
1\\
-\lambda\cos(2\theta)\\
-\lambda\sin(2\theta)
\end{bmatrix}$$
With a little bit of vector algebra, we can solve above three equations to get
$$\Lambda = \frac{U(A)\times U(B) + U(B)\times U(C) + U(C)\times U(A)}{U(A)\cdot (U(B) \times U(C))}$$
It is not hard to verify the $U(\cdot), V(\cdot)$ satisfy following identities
- $U(P)\cdot U(Q) = V(P)\cdot V(Q) = 2(P\cdot Q)^2$.
- $U(P)\cdot V(Q) = V(P)\cdot U(Q) = 2(P\times Q)^2$.
- $U(P)\times U(Q) = 2(P\times Q)\left[V\left(\frac{P+Q}{2}\right)-V\left(\frac{P-Q}{2}\right)\right]$
Using these identities, the denominator simplifies to
$$\small \begin{align} & U(A)\cdot(U(B)\times U(C))\\
= & 2(B\times C)\,U(A)\times\left[V\left(\frac{B+C}{2}\right) - V\left(\frac{B-C}{2}\right)\right]\\
= & 2(B\times C)\left[ 2\left(A \times \frac{B+C}{2}\right)^2 - 2\left(A \times \frac{B-C}{2}\right)^2\right]\\
= & -4(B \times C)(C\times A)(A\times B)
\end{align}$$
Let $\Delta$ be the area of $\triangle ABC$ and $u,v,w$ be the barycentric coordinates of the origin with respect to $\triangle ABC$. We have
$$2\Delta u = B \times C,\quad 2\Delta v = C \times A\quad\text{ and }\quad 2\Delta w = A \times B$$
We find the denominator equals to $-32\Delta^3 uvw$.
Doing similar thing to the numerator and let $D, E, F$ be the midpoint of $BC, CA$ and $AB$, the numerator becomes
$$\small 4\Delta \left[ u (V(D) - V(E-F)) + v(V(E) - V(F-D)) + w(V(F)-V(D-E)) \right]$$
This is a little bit clumsy to write down. Let us use the notation $\sum_{cyc}$ to indicate a cyclic sum of over parameter set
$$u \to v \to w,\quad A \to B \to C\quad\text{ and }\quad D \to E \to F$$
In this new notation, the equation for $\Lambda$ becomes
$$\Lambda = - \frac{1}{8\Delta^2 uvw} \sum_{cyc} u(V(D) - V(E-F))$$
Switch to other coordinate system where the center, let's call it $Z$, is no longer the origin, the equation for the center becomes
$$\Lambda = -\frac{1}{8\Delta^2 uvw}\sum_{cyc}u(V(Z-D) - V(E-F))\tag{*1}$$
To proceed further, we will do two things.
We will switch to a coordinate system where $N$, the nine-point center of $\triangle ABC$ is the origin. In this coordinate system, $D,E,F$ will be lying on a circle centered at origin $N$ with radius $\frac{R}{2}$ ($R$ is the circumradius of $\triangle ABC$).
We will identify the Euclidean plane with the complex plane, we will use the lower case letter to denote the complex number corresponds to a point. e.g. $Z$ becomes $z$ and $D, E, F$ become $d, e, f$.
After this, we can repress $(*1)$ as two equations
$$\begin{align}
\Lambda_1 = \frac{1}{\sigma^2}
&= -\frac{1}{8\Delta^2 uvw} \sum_{cyc} u(|z-d|^2 - |e - f|^2)\tag{*2a}\\
\Lambda_2 + i\Lambda_3 = -\frac{\lambda}{\sigma^2} e^{2i\theta}
&= +\frac{1}{8\Delta^2 uvw} \sum_{cyc} u((z-d)^2 - (e - f)^2)\tag{*2b}
\end{align}
$$
Let $o = d + e + f$ (the point corresponding to it is the circumcenter $O$ of $\triangle ABC$) and $u',v',w'$ the barycentric coordinate of $Z$ with respect to $\triangle DEF$. We have
$$u = \frac{1-u'}{2},\cdots
\implies \sum_{cyc} ud = \frac12\sum_{cyc}(1-u')d = \frac12(o-z)$$
Notice
$$\begin{align}
|z-d|^2 - |e-f|^2 &= |z-d|^2 + |e+f|^2 - 2|e|^2 - 2|f|^2\\
&= |z-d|^2+|o-d|^2 - 4|d|^2\\
&= |z|^2 -d(\overline{z+o}) - \bar{d}(z+o) + |o|^2 - 2|d|^2
\end{align}$$
Multiply by $u$ and take cyclic sum, we get
$$\begin{align} \sum_{cyc}u(|z-d|^2 - |e-f|^2)
= & |z|^2 - \frac{(o - z)(\overline{z+o})}{2} - \frac{(\overline{o-z})(z+o)}{2} + |o|^2 - 2|d|^2\\
= & 2|z|^2 - \frac{R^2}{2}\end{align}$$
Equation $(*2a)$ becomes
$$\bbox[border:1px solid blue;padding: 1em;]{z\bar{z} - \frac{R^2}{4} + \frac{4\Delta^2 uvw}{\sigma^2} = 0} \tag{*3a}$$
Since $u,v,w$ are linear in $(x,y) = (\Re z,\Im z)$. This describe a cubic curve
in the Euclidean plane and it is the locus of the center when the size parameter $\sigma$ is held fixed.
Let $\Omega = \frac{def}{|d|^2} = \frac{4def}{R^2}$. By a similar procedure, we have
$$\begin{align}(z-d)^2 - (e-f)^2
&= (z-d+e+f)(z-d-e-f) = (z-o+2e)(z-o+2f)\\
&= (z-o)^2+2(e+f)(z-o) + 4ef\\
&= z^2 - o^2 - 2d(z-o) + 4\Omega \bar{d}
\end{align}
$$
Multiply by $u$ and take cyclic sum, we get
$$\begin{align}\sum_{cyc} u((z-d)^2 - (e-f)^2)
&= z^2 - o^2 - (o-z)(z-o) + 2\Omega(\bar{o} - \bar{z})\\
&= 2\left(z(z-o) - \Omega(\bar{z} - \bar{o})\right)\end{align}$$
Equations $(*2b)$ becomes
$$z(z-o) - \Omega(\bar{z} - \bar{o})
= -\lambda e^{2i\theta}\frac{4\Delta^2 uvw}{\sigma^2}$$
Compare this with equation $(*3a)$, we obtain
$$z(z-o) - \Omega(\bar{z} - \bar{o}) = \lambda e^{2i\theta} \left(z\bar{z} - \frac{R^2}{4}\right)$$
Taking absolute value and square, we obtain a quartic curve
$$\bbox[border:1px solid blue;padding: 1em]{|z(z-o) - \Omega(\bar{z}-\bar{o})|^2 = \lambda^2 \left|z\bar{z} - \frac{R^2}{4}\right|^2}\tag{*3b}
$$
This is the locus of the center when the eccentricity $\epsilon$ is help fixed.
observation I - locus for $\sigma$.
All locus of $\sigma$ passes through $D,E,F$ and the three foots of $\triangle DEF$.
When orthocenter $H$ of $\triangle ABC$ belongs to interior of $\triangle DEF$. $\sigma$ takes a local minimum at $H$. This should corresponds to a special circumellipse associated with $\triangle ABC$. I'm unable to figure out what that is. Anyone has any idea?
observation II - locus for $\lambda$.
When $\epsilon = 0 \implies \lambda = 0$, the locus reduces to a single point $O$. The circumcenter of $\triangle ABC$.
When $\epsilon = 1 \implies \lambda = 1$, the quartic terms in Eq. ($*3b$) cancel out. The corresponding curve reduces a union of $3$ lines through the midpoints. i.e. the lines $DE$, $EF$ and $FD$.
These $3$ lines split the plane into $7$ regions. When one varies $\epsilon$ (and hence $\lambda$) from $0$ to $1$. The locus sweep across the interior of
$4$ out of the $7$ regions ( $\triangle DEF$ or the $3$ cone with apex at $D, E, F$). This means in order for a point to be center for an ellipse, it need to be
either $D, E, F$ or belongs to the interior of above $4$ regions.
Infinitely many loci of $\lambda$ passing through $D, E, F$. Combine with the observation $I_1$, $D, E, F$ are centers for infinitely many circumellipses of $A,B,C$.
To be continued???
Best Answer
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
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) $
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} $
Update the estimate of the solution by the Newton-Raphson recurrence
$$ r_{0}^{(k+1)} = r_{0}^{(k)} - J^{-1} y^{(k)} $$
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) $