I'm reading over some C++ code which does a weighted maximum likelihood fit of a 2-D normal distribution. Details aside, the code estimates the positive definite variance matrix:
$$ C=\begin{bmatrix} C_{xx} & C_{xy} \\ C_{xy} & C_{yy} \end{bmatrix},\,\, C>0 $$
Now, on line 81 the code calculates the "dominant direction" of $C$ as:
$$ \varphi = \frac{1}{2}\text{atan2}\left(-2C_{xy} , C_{yy}-C_{xx}\right) $$
Where does this formula come from (the code comments this formula as "find dominant direction via SVD")? Please derive it or point to a resource which derives it.
NB: I understand $\varphi$ to be the angle from the $x$ axis to the dominant output vector, like so:
I found the answer, buried inside the original C-code. Here's the mathematical derivation.
Suppose we want to find the dominant direction of the general 2D matrix
$$ A=\begin{bmatrix} A_{00} & A_{01} \\ A_{10} & A_{11} \end{bmatrix}\in\mathbb R^2. $$
The SVD of a 2D matrix can always be expressed as:
$$ A=USV^T=\begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix} \begin{bmatrix} \overline\sigma & 0 \\ 0 & \underline\sigma \end{bmatrix} \begin{bmatrix} \cos(\varphi) & \sin(\varphi) \\ -\sin(\varphi) & \cos(\varphi) \end{bmatrix}. $$
Since $\overline\sigma\ge\underline\sigma$, our task is basically to find $\varphi$ (then, $\begin{bmatrix} \cos(\varphi) & \sin(\varphi) \end{bmatrix}^T$ will give the dominant direction). To do so, we expand the SVD (i.e. multiply out the 3 matrices $U$, $S$ and $V^T$) and recognize that
$$ \begin{array}{rclcl} B_0 & := & A_{00}+A_{11} & = & (\overline\sigma+\underline\sigma)\cos(\varphi-\theta), \\ B_1 & := & A_{00}-A_{11} & = & (\overline\sigma-\underline\sigma)\cos(\varphi+\theta), \\ B_2 & := & A_{01}+A_{10} & = & (\overline\sigma-\underline\sigma)\sin(\varphi+\theta), \\ B_3 & := & A_{01}-A_{10} & = & (\overline\sigma+\underline\sigma)\sin(\varphi-\theta). \end{array} $$
Hence we can solve for $\varphi$,
$$ \varphi = \frac{1}{2}\left(\text{atan}\left(\frac{B_3}{B_0}\right)+\text{atan}\left(\frac{B_2}{B_1}\right)\right). $$
For the particular case of positive definite matrix $C$ in my original post, substituting in the values of $A_{ij}$, $i,j\in\{0,1\}$, gives precisely
$$ \varphi =\frac{1}{2}\text{atan}\left(\frac{2C_{xy}}{C_{xx}-C_{yy}}\right). $$
Now we can use the $\text{atan2}$ function as we desire to have $\varphi\in[-\pi,\pi]$. This is exactly the expression found in C++!