Source: Exercise 2.3.18 (p.54) from Convex functions: constructions, characterizations and counterexamples, J.M. Borwein & J.D. Vanderwerff (2010).
Consider $E:=\{(x,y):x^2/a^2+y^2/b^2=1\}$ in standard form. Show that the best approximation is:
$$P_E\,(u,v)=\left(\frac{a^2u}{a^2-t},\frac{b^2v}{b^2-t}\right)$$
where $t$ solves $\frac{a^2u^2}{(a^2-t)^2}+\frac{b^2v^2}{(b^2-t)^2}=1$.
More explicitly, we have the decomposition
$$\begin{pmatrix}\cos\,t\\\cos\,t+\sin\,t\end{pmatrix}=\begin{pmatrix}\cos\,\lambda&-\sin\,\lambda\\\sin\,\lambda&\cos\,\lambda\end{pmatrix}\cdot\begin{pmatrix}\sqrt{1+\phi}\cos(t+\eta)\\\sqrt{2-\phi}\sin(t+\eta)\end{pmatrix}$$
where $\tan\,\lambda=\phi$, $\tan\,\eta=1-\phi$, and $\phi=\dfrac{1+\sqrt{5}}{2}$ is the golden ratio. You can check that both your original parametric equations and the new decomposition both satisfy the Cartesian equation $2x^2-2xy+y^2=1$. What the decomposition says is that your curve is an ellipse with axes $\sqrt{1+\phi}$ and $\sqrt{2-\phi}$, with the major axis inclined at an angle $\lambda$.
If we take the linear algebraic viewpoint, as suggested by Robert in the comments, what the decomposition given above amounts to is the singular value decomposition (SVD) of the shearing matrix; i.e.,
$$\begin{pmatrix}1&0\\1&1\end{pmatrix}=\begin{pmatrix}\cos\,\lambda&-\sin\,\lambda\\\sin\,\lambda&\cos\,\lambda\end{pmatrix}\cdot\begin{pmatrix}\sqrt{1+\phi}&\\&\sqrt{2-\phi}\end{pmatrix}\cdot\begin{pmatrix}\cos\,\eta&\sin\,\eta\\-\sin\,\eta&\cos\,\eta\end{pmatrix}^\top$$
The SVD is in fact an excellent way to look at how a matrix transformation geometrically affects points: the two orthogonal matrices on the left and right can be thought of as rotation matrices, reflection matrices, or products thereof, and the diagonal matrix containing the singular values amounts to nothing more than a scaling about the axes of your coordinate system.
Best Answer
There doesn't seem to be a nice closed-form solution to this, but the following might be helpful..
For an ellipse $\frac {x^2}{a^2}+\frac {y^2}{a^2}=1$, the equation of the normal at point $P(a \cos\theta, b\sin\theta)$ on the ellipse is
$$\frac {a\sin\theta}{b\cos\theta}x-y=a\left(\frac{a^2-b^2}{ab}\right)\sin\theta\tag{1}$$
A hyperbola which just touches or does not intersect at all with the ellipse above has the equation $xy=\frac {mab}2$ where $m\ge 1$. The equation of the normal at point $Q (v,\frac {mab}{2v})$ on the hyperbola is $$\frac {2v^2}{mab}x-y=v\left(\frac {2v^2}{mab}-\frac {mab}{2v^2}\right)\tag{2}$$
The minimum distance between the ellipse and hyperbola is the distance $PQ$ when $(1)=(2)$, i.e. both normals are the same line.
As the coefficients of $y$ are the same in both $(1),(2)$, equating coefficients of $x$ in $(1),(2)$ gives $$v=a\sqrt \frac{m\tan\theta}{2}\tag{3}$$ This relationship ensures that the tangents and normals at $P,Q$ are parallel to each other respectively (but the normals are not necessarily the same line).
Putting $(3)$ in $(2)$ gives
$$\left(\frac ab \tan\theta\right)x-y=a\sqrt{\frac{m\tan\theta} 2}\left(\frac ab\tan\theta-\frac ba\cot\theta\right)\tag{4}$$
To ensure that both normals are the same line, we need to equate RHS of $(1),(4)$. This gives $$\left(\frac{a^2-b^2}{ab}\right)\sin\theta=\sqrt{\frac{m\tan\theta}2}\left(\frac ab \tan\theta-\frac ba\cot\theta\right)$$ which is equivalent to $$(a^2-b^2)\sin\theta=\sqrt{\frac{m\tan\theta}2}\left(\frac{a^2\sin^2\theta-b^2\cos^2\theta}{\sin\theta\cos\theta}\right)\tag{4}$$
Solve numerically $\theta$ in $(4)$, find corresponding value $v$ in $(3)$, then calculate $PQ$. This should give the minimum distance between the ellipse and hyperbola.
See desmos implementation here.
In the trivial case where $a=b$ (i.e. ellipse is a circle), then $\theta=\frac \pi 4$ and $v=a\sqrt{\frac m2}$ . This gives $P=\left(\frac a{\sqrt{2}}, \frac a{\sqrt{2}}\right)$ and $Q=\left(a\sqrt{\frac m2}, a\sqrt{\frac m2}\right)$ and the distance $PQ=a(\sqrt m-1)$.