I rewrote this answer completely, see its history for earlier versions.
Change the coordinate system in such a way that the ellipse becomes a circle of radius $1$. The following transformation achieves this:
$$ f: \begin{pmatrix}x\\y\end{pmatrix} \mapsto
\begin{pmatrix}\frac1{r_x}&0\\0&\frac1{r_y}\end{pmatrix}\cdot
\begin{pmatrix}
\cos\alpha & \sin\alpha \\
-\sin\alpha & \cos\alpha \\
\end{pmatrix}\cdot
\begin{pmatrix}x\\y\end{pmatrix}
= \begin{pmatrix}
\frac{x\cos\alpha+y\sin\alpha}{r_x} \\
\frac{y\cos\alpha-x\sin\alpha}{r_y}
\end{pmatrix}
$$
Apply this transformation to the end points to obtain new coordinates $(x_1',y_1')$ and $(x_2',y_2')$. The center of the ellipse corresponds to a point which has distance $1$ from these two points. So now you have to intersect two circles of radius $1$ around these. To do this, compute the difference vector
$$ v = \begin{pmatrix}v_x\\v_y\end{pmatrix} =
\begin{pmatrix}x_1'\\y_1'\end{pmatrix} -
\begin{pmatrix}x_2'\\y_2'\end{pmatrix}
$$
A circle of radius $1$ at the origin will intersect another circle of radius $1$ at point $v$ in two points which form a line that will pas through the midpoint between the centers, i.e. through the point $\frac12v$. A generic point on that line can be written as
$$ q = \begin{pmatrix}q_x\\q_y\end{pmatrix} =
\frac12\begin{pmatrix}v_x\\v_y\end{pmatrix} +
\lambda\begin{pmatrix}-v_y\\v_x\end{pmatrix} $$
To find those points where that line intersects the circle, compute the squared distance from the origin, and set that equal to one.
$$ \lVert q\rVert^2 = q_x^2 + q_y^2 = 1 $$
Solving this for $\lambda$ will result in two solutions, corresponding to the two points of intersection:
$$ \lambda = \pm\sqrt{\frac{1}{v_x^{2} + v_y^{2}} - \frac{1}{4}} $$
One of these will correspond to the case where the two flags agree, the other to that where they disagree. Try them both until you know which one is which. Once you have made that decision, you now have the point $q$ which corresponds to the center of your ellipse. The ellipse itself will be represented by a circle of radius $1$ around that point. The tangent to that circle through the origin will be the line
$$ q_xx + q_yy = 0 $$
You have to shift that line by $(x_2',y_2')$ and undo $f$ to obtain the equation of the tangent in your original coordinate system. To do that, you can combine all the coordinate transformations done so far into a single affine transformaton matrix.
\begin{align*}
M &=
\begin{pmatrix}
\frac1{r_x} & 0 & 0 \\
0 & \frac1{r_y} & 0 \\
0 & 0 & 1
\end{pmatrix} \cdot
\begin{pmatrix}
\cos\alpha & \sin\alpha & 0 \\
-\sin\alpha & \cos\alpha & 0\\
0 & 0 & 1
\end{pmatrix} \cdot
\begin{pmatrix}
1 & 0 & -x_2 \\
0 & 1 & -y_2 \\
0 & 0 & 1
\end{pmatrix} \\&=
\begin{pmatrix}
\frac{\cos\alpha}{r_x} &
\frac{\sin\alpha}{r_x} &
-\frac{x_2\cos\alpha + y_2\sin\alpha}{r_x} \\
-\frac{\sin\alpha}{r_y} &
\frac{\cos\alpha}{r_y} &
\frac{x_2\sin\alpha - y_2\cos\alpha}{r_y} \\
0 & 0 & 1
\end{pmatrix}
\end{align*}
Then you can take the parameters of the line and multiply them by the transpose of that matrix, or any multiple thereof.
$$ r_xr_yM^T\cdot\begin{pmatrix}q_x\\q_y\\0\end{pmatrix} = \begin{pmatrix}
-r_xq_y\sin\alpha + r_yq_x\cos\alpha \\
r_xq_y\cos\alpha + r_yq_x\sin\alpha \\
r_xx_2q_y\sin\alpha - r_xq_yy_2\cos\alpha -
r_yq_xx_2\cos\alpha - r_yq_xy_2\sin\alpha
\end{pmatrix}
$$
So the equation of your tangent in the original coordinate system will be
\begin{multline*}
\left( -r_xq_y\sin\alpha + r_yq_x\cos\alpha \right)x +
\left( r_xq_y\cos\alpha + r_yq_x\sin\alpha \right)y \\ = -
\left( r_xx_2q_y\sin\alpha - r_xq_yy_2\cos\alpha -
r_yq_xx_2\cos\alpha - r_yq_xy_2\sin\alpha \right)
\end{multline*}