How can I get the axes of the polarization ellipse from the Jones vector of the light? I am analyzing the polarization state of a monochromatic, coherent light source, for which I know the Jones vector of the polarization,
$$
\mathbf E
=\begin{pmatrix}E_x\\E_y\end{pmatrix}
=\begin{pmatrix}|E_x|e^{i\varphi_x}\\|E_y|e^{i\varphi_y}\end{pmatrix},
$$
and I would like to expand it in terms of a major and a minor axis of ellipticity, i.e. in the form
$$
\mathbf E=
e^{i\varphi}\left(
A \hat{\mathbf u} + i B\hat{\mathbf v}
\right)
=
e^{i\varphi}\left(
A
\begin{pmatrix}\cos(\theta)\\ \sin(\theta)\end{pmatrix}
+ i B
\begin{pmatrix}-\sin(\theta)\\ \cos(\theta)\end{pmatrix}
\right),
$$
or as shown graphically as follows:

Image source
Wikipedia provides a multi-step procedure going through the Stokes parameters, but I'm thinking there is surely a cleaner and more direct way to get $A$, $B$, $\hat{\mathbf u}$, $\hat{\mathbf v}$, $\theta$, and the components $A \hat{\mathbf u}$ and $B\hat{\mathbf v}$, from $E_x$ and $E_y$, and it's not particularly obvious from the search results I can find. What's the cleanest way to do this?

To be clear: what I think is lacking from the existing resources, and what the question is directly asking for, is an explicit set of connections, as simple as possible, for the named parameters (all of $A$, $B$, $\hat{\mathbf u}$, $\hat{\mathbf v}$, $\theta$, and the components $A \hat{\mathbf u}$ and $B\hat{\mathbf v}$), in terms of the Cartesian components $E_x$ and $E_y$. Schemes that simply send to some other set of complex manipulations are already available from Wikipedia and are not what the question is asking for.
 A: Let me try a second time. I use https://math.stackexchange.com/questions/1204131/converting-a-rotated-ellipse-in-parametric-form-to-cartesian-form as a resource.
Prior Manipulations
The physical electric field is
$$\mathbf{E}_{phys} = Re\left[\mathbf{E} e^{i\omega t}\right] = Re\left[\begin{pmatrix}|E_x|e^{i\varphi_x}\\|E_y|e^{i\varphi_y}\end{pmatrix}e^{i\omega t}\right] = \begin{pmatrix}|E_x|\cos(\omega t + \varphi_x)\\|E_y| \cos(\omega t + \varphi_y)\end{pmatrix} = \begin{pmatrix}|E_x|\left[\cos(\omega t)\cos(\varphi_x) - \sin(\omega t) \sin(\varphi_x) \right]\\|E_y| \left[\cos(\omega t)\cos(\varphi_y) - \sin(\omega t) \sin(\varphi_y) \right]\end{pmatrix}$$
This is a parametric equation for an ellipse, which is traced by the electric field.
Major axes angle
Let
$$a=|E_x|\cos(\varphi_x)$$
$$b=|E_x|\sin(\varphi_x)$$
$$c=|E_y|\cos(\varphi_y)$$
$$d=|E_y|\sin(\varphi_y)$$
Then by comparison with the linked math.SE question the accepted answer states that the major and minor axes point on the ellipse (which is centred on the origin) fulfill
$$\omega t={1\over2}\arctan{2(ab+cd)\over(a^2+c^2)-(b^2+d^2)}+{k\pi\over2}\qquad(0\leq k\leq3)\ .$$
Expansion required by the OP
Indeed this quantity is the angle $\theta$ in the expansion appearing in the OPs question. However depending on which value for $k$ is chosen in may be $\pi/2 - \theta$, a case distinction depending on the sector of the arctan is necessary.
This of course also yields $\hat{u}$ and $\hat{v}$, so the expansion can now be easily obtained by projecting the Jones vector this basis.
The formulae may be along, but constitute a close form solution to the problem up to the case distinction of choosing $k$. I do not see how a simpler solution can exist, since the formula give for the angle given does not seem algebraically reducible.
A: The cleanest way to do this is offered by Michael Berry in the paper

Index formulae for singular lines of polarization. M V Berry, J. Opt. A: Pure Appl. Opt. 6, 675–678 (2004), author eprint.

In Berry's notation, the electric field can be written as
$$
\mathbf E=\mathbf P +i \mathbf Q = e^{i\gamma} \left(\mathbf A+i\mathbf B\right),
$$
where $\mathbf P$, $\mathbf Q$, $\mathbf A$ and $\mathbf B$ are real-valued vectors, $\mathbf A$ and $\mathbf B$ are respectively the major and minor axes of the polarization ellipse, and those two are defined up to a sign by $\mathbf A\cdot\mathbf B=0$ and $|\mathbf A|\geq |\mathbf B|$. With this notation, the polarization axes and the phase are defined as
$$
\gamma = \frac12 \arg(\mathbf E\cdot\mathbf E)
\quad\text{and}\quad
\mathbf A+i\mathbf B = \frac{\sqrt{\mathbf E^*\cdot \mathbf E^*}}{\left|\sqrt{\mathbf E^*\cdot \mathbf E^*}\right|}\mathbf E,
$$
or in other words
$$
\mathbf A = \frac{1}{\left|\sqrt{\mathbf E^*\cdot \mathbf E^*}\right|}\mathrm{Re}\mathopen{}\left[\sqrt{\mathbf E^*\cdot \mathbf E^*} \: \mathbf E\right]\mathclose{}
\quad\text{and}\quad
\mathbf B = \frac{1}{\left|\sqrt{\mathbf E^*\cdot \mathbf E^*}\right|}\mathrm{Im}\mathopen{}\left[\sqrt{\mathbf E^*\cdot \mathbf E^*} \: \mathbf E\right]\mathclose{}.
$$
There is an obvious sign ambiguity, in that flipping both $\mathbf A$ and $\mathbf B$ and adding $\pi$ to $\gamma$ will not change anything (i.e. rotating the polarization ellipse by 180° is equivalent to adding a phase), which is reflected in the branch cuts of both the argument and the square root functions. These naturally mesh together so long as both branch cuts are taken on the same cut, ideally along the negative real axis.
As another tricky point, one should note that these formulas are not defined when $\mathbf E\cdot\mathbf E=0$, which corresponds to circular polarization; in this case both the polarization axes, as well as the phase $\gamma$ at the major axis, are ill-defined, so this is not a problem.
As an added bonus, this approach also naturally gives the direction of the normal to the plane of the polarization ellipse, in the form
$$
\mathbf C 
= \frac12 \mathrm{Im}\mathopen{}\left(\mathbf E^*\times\mathbf E\right)\mathclose{}
=\mathbf P\times\mathbf Q
=\mathbf A\times\mathbf B,
$$
where the cross product $\mathbf E^*\times\mathbf E$ is naturally imaginary, as its conjugate is minus itself. Of course, this will vanish if $\mathbf E$ and $\mathbf E^*$ (or $\mathbf P$ and $\mathbf Q$) are linearly dependent, which corresponds to linear polarization; in this case, $\mathbf B$ will vanish, because $\sqrt{\mathbf E^*\cdot \mathbf E^*} \: \mathbf E$ is naturally real.
Berry credits

Polarization singularities in paraxial vector fields: morphology and statistics. M R Dennis, Opt. Commun. 213,
  201–21 (2002), eprint.

for this form, and that reference contains a fuller proof of how and why the decomposition works. 
This is, in fact, rather simple, once you realize that the decomposition as $\mathbf E = e^{i\gamma} \left(\mathbf A+i\mathbf B\right)$, as above, must exist, because under those conditions the dot product reduces to
$$
\mathbf E\cdot\mathbf E
=e^{2i\gamma} \left(\mathbf A+i\mathbf B\right)\cdot \left(\mathbf A+i\mathbf B\right)
=e^{2i\gamma}(\mathbf A^2-\mathbf B^2),
$$
where $\mathbf A^2-\mathbf B^2$ is real and positive, so taking the argument of both sides naturally gives the phase as $2\gamma=\arg(\mathbf E\cdot\mathbf E).$ Similarly, taking the modulus of that equation returns $\mathbf A^2-\mathbf B^2=|\mathbf E\cdot\mathbf E|$, so we can simply get the phase factor as
$$
e^{2i\gamma} 
=\frac{\mathbf E\cdot\mathbf E}{\mathbf A^2-\mathbf B^2}
=\frac{\mathbf E\cdot\mathbf E}{|\mathbf E\cdot\mathbf E|}
,\ \text{so}\ 
e^{i\gamma} 
=\frac{\sqrt{\mathbf E\cdot\mathbf E}}{|\sqrt{\mathbf E\cdot\mathbf E}|}
,\ \text{and therefore}\ 
e^{-i\gamma} 
=\frac{\sqrt{\mathbf E^*\cdot\mathbf E^*}}{|\sqrt{\mathbf E^*\cdot\mathbf E^*}|};
$$
the characterization for $\mathbf A+i\mathbf B$ then follows from $\mathbf E = e^{i\gamma} \left(\mathbf A+i\mathbf B\right)$ by simply dividing by $e^{i\gamma}$.
A: I think you can accomplish this with a singular-value decomposition.  I'll start by writing $\mathbf{E}$ in the following form
$\mathbf{E} = E_x \mathbf{\hat{x}} + E_y \mathbf{\hat{y}} = E_u \mathbf{\hat{u}} + E_v \mathbf{\hat{v}}$
where $E_x$ and $E_y$ are known complex numbers. In general $E_x E_y^* \ne 0$, but we want $E_u E_v^*=0$. We can write this in matrix form as
$\left( \begin{array}{cc}  \mathbf{\hat{x}} & \mathbf{\hat{y}} \end{array} \right)
\left( \begin{array}{c}  E_x \\ E_y \end{array} \right) = 
\left( \begin{array}{cc}  \mathbf{\hat{u}} & \mathbf{\hat{v}} \end{array} \right)
\left( \begin{array}{c}  E_u \\ E_v \end{array} \right)$
We should be able to write
$\left( \begin{array}{cc}  \mathbf{\hat{x}} & \mathbf{\hat{y}} \end{array} \right) = 
\left( \begin{array}{cc}  \mathbf{\hat{u}} & \mathbf{\hat{v}} \end{array} \right) \Lambda^T$
where $\Lambda$ is a real, orthogonal matrix which describes a rotation of coordinate systems. It is currently unknown. Since we assume that $\mathbf{\hat{x}}$ and $\mathbf{\hat{y}}$ are orthogonal, then $\mathbf{\hat{u}}$ and $\mathbf{\hat{v}}$ will also be orthogonal. And we can also write
$
\left( \begin{array}{c}  E_x \\ E_y \end{array} \right) = 
A
\left( \begin{array}{c}  1 \\ i \end{array} \right)$
where $A = \left( \begin{array}{cc}  a & b\\c & d \end{array} \right)$ is a known real matrix. We can decompose $A$ using a singular value decomposition
$A = U S V^T$
where $U$,$S$, and $V$ are real matrices. $U$ and $V$ are orthogonal matrices
and $S=\left(\begin{array}{cc}A&0\\0&B\end{array}\right)$ where $A$ and
$B$ are positive with $A\ge B$.
This all yields
$ \Lambda^T
U S V^T \left( \begin{array}{c}  1 \\ i \end{array} \right) = 
\left( \begin{array}{c}  E_u \\ E_v \end{array} \right)$
If we choose $\Lambda = U$, which gives $\mathbf{\hat{u}}$ and 
$\mathbf{\hat{v}}$, then this reduces to 
$ S V^T \left( \begin{array}{c}  1 \\ i \end{array} \right) = 
\left( \begin{array}{c}  E_u \\ E_v \end{array} \right)$
Since $V$ is an orthogonal matrix, so is $V^T$, and $S$ is 
diagonal and real, so $E_u E_v^*=0$. 
