This formula rescales a standard ellipse $(\cos(\theta + d/2), \cos(\theta - d/2))$ (inscribed within the unit square) by the diagonal matrix $(x_0, y_0)$. By symmetry, the values $\theta = 0$ and $\theta = \pi/2$ correspond to vertices of this standard ellipse, allowing us to find their coordinates (and thus the lengths of the semi-axes), whence we easily deduce its equation is $x^2 + y^2 - 2 \rho x y = 1 - \rho^2$. Applying the diagonal matrix gives the conventional implicit form
$$\left(\frac{x}{x_0}\right)^2 + \left(\frac{y}{y_0}\right)^2 - 2 \rho \frac{x}{x_0} \frac{y}{y_0} = 1 - \rho^2 \text{.}$$
From here you can look up anything you want.
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
The minor axis will be perpendicular to the line of the major axis.
To find the length of the semi-major axis , find the distance between a point on intersection of the line above with the ellipse and the origin . Likewise the length of the semi-minor axis will be the distance between the point of intersection of the second line with the ellipse and the origin. I obtained this value by finding the maximum of x^2 + y^2 on the ellipse, ( this can be done in many ways, one is the method of lagrange multipliers.)
I'll include my sketch. You can take it down off the screen for a closer look. In regard to rotation of the axis, find the angle the red line makes with the x-axis , this will give you a good idea of what to try in terms of rotating your coordinate system.