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}$.