The equation of an ellipse in the plane is given by
$ (r - E)^T Q (r - E) = 1 \hspace{20pt} (1)$
where $E$ is the center, and $Q$ is a $2 \times 2$ symmetric matrix.
At any point on the ellipse the gradient vector which is the normal (perpendicular) vector to the ellipse curve is given by
$ g = 2 Q (r - E) \hspace{20pt} (2) $
If the tangency point on $AB$ is $r_1$ then the gradient vector at $r_1$ will be parallel to the normal to the line segment $AB$, that is,
$ Q (r_1 - E) = \alpha n_1 \hspace{20pt} (3)$
where $n_1$ is known (it is the normal vector to $AB$). From this it follows that
$ r_1 - E = \alpha Q^{-1} n_1 \hspace{20pt} (4)$
Since $r_1$ is on the ellipse we can plug this expression into the equation of the ellipse, to deduce that
$\alpha = \dfrac{1}{\sqrt{n_1^T Q^{-1} n_1 } } \hspace{20pt} (5)$
On the other hand, premultiplying $(4)$ by $n_1^T$ and using $(5)$
$n_1^T (r_1 - E) = \sqrt{n_1^T Q^{-1} n_1 } $
The equation of the line segment $AB$ is $n_1^T (r - A) = 0 $, and we have
$n_1^T (r_1 - E) = n_1^T (r_1 - A + A - E) = n_1^T (A - E )$
because $n_1^T (r_1 - A) = 0$. Thus we have
$n_1^T (A - E) = \sqrt{ n_1^T Q^{-1} n_1 } \hspace{20pt} (6.1)$
Similar equations can be written for the other three line segments $BC, CD, DA$,
whose normals are $n_2, n_3, n_4$:
$n_2^T (B - E) = \sqrt{ n_2^T Q^{-1} n_2 } \hspace{20pt}(6.2) $
$n_3^T (C - E) = \sqrt{n_3^T Q^{-1} n_3 } \hspace{20pt} (6.3) $
$n_4^T (D - E) = \sqrt{n_4 ^T Q^{-1} n_4 } \hspace{20pt} (6.4)$
Only $3$ equations out of $4$ are needed to find the matrix $Q^{-1}$. Using the theorem mentioned in the question, a matrix $Q$ satisfying any three of the equations $(6.1) - (6.4)$ will satisfy the fourth equation if the center $E$ is on the line segment between the midpoint of $AC$ and the midpoint of $BD$. ( I don't have a proof of that but I verified it numerically).
Solving any three of equations $(6.1) - (6.4)$ can be done quite easily. Let matrix
$M = Q^{-1} = \begin{bmatrix} M_{11} && M_{12} \\ M_{12} && M_{22} \end{bmatrix} $
Then from $(6.1)$,
$n_1^T Q^{-1} n_1 = \left( n_1^T (A - E) \right)^2 \hspace{20pt} (7) $
The right hand side of $(7)$ is known, while the left hand side is linear in the entries of $Q^{-1} $, namely,
$n_1^T Q^{-1} n_1 = M_{11} n_{1x}^2 + M_{22} n_{1y}^2 + 2 M_{12} n_{1x} n_{1y} \hspace{20pt} (8) $
Two more equations can be written creating a $3 \times 3$ in the unknowns $M_{11}, M_{22}, M_{12} $, and solved.
Once $Q^{-1}$ is found, the ellipse specification is complete. Finding all the parameters of the ellipse follows directly by diagonalizing matrix $Q$ as follows:
$Q = R D R^T \hspace{20pt} (9)$
It can be assumed that $D_{11} \le D_{22}$. The semi-major axis $a = \dfrac{1}{\sqrt{D_{11}}}$ and the semi-minor axis $b = \dfrac{1}{\sqrt{D_{22}}} $
The focii are along the major axis, and are given by
$F_1 = E + a e R_1 $ and $F_2 = E - a e R_1 \hspace{20pt} (10) $
where $R_1$ is the first column vector of the matrix $R$ specified in $(9)$, and $e$ is the eccentricy,
$e = \sqrt{ 1 - \left( \dfrac{b}{a} \right)^2 }$
The figure below shows an example with the above method.