Principal component analysis in two dimensions During my studies, I stumbled upon the following exercise:
We have the following joint probability distribution:
$$p(x,y) = p(x) p(y|x)$$
$$p(x) = \mathcal{N}(0,1), p(y \mid x) = \frac{1}{2} \delta(y -x) + \frac{1}{2} \delta(y+x)$$
where $\delta(\cdot)$ is the Direc delta function. The exercise then asks to find the principal components of $p(x,y)$. 
It is hinted that this is equivalent to finding the parameters $\theta \in [0, 2 \pi [$ that maximize the variance of the projected data: $z(\theta) = x \cos(\theta) + y \sin(\theta)$, since in linear component analysis for two-dimensional probability distributions the set of possible directions to look for in $\mathbb{R}^2$ is given by:
$\{ \begin{pmatrix} \cos(\theta) \\ \sin(\theta) \end{pmatrix}, 0 \leq \theta \leq 2\pi \}$. 
Usually I would take the Lagrangian and the derive the maximum, but I don't know how the Langrangian would look like in this case. How would I go about solving this?
 A: Since $E[z(\theta)]=0$, we have $$\operatorname{var}(z(\theta))=E[z(\theta)^2]=E[\cos(\theta)^2x^2+\sin(\theta)^2y^2+2\cos\theta\sin\theta xy]$$ Here, by definition, $E[x^2]=E[y^2]=1$. So, the expression reduces to:
$$\operatorname{var}(z(\theta))=1+\sin(2\theta) E[xy]$$
We also have $E[xy]=E[E[xy|x]]=E[xE[y|x]]=0$, which means $\operatorname{var}(z)=1$. 
This is an interesting dataset, the samples follow $y=x$ and $y=-x$ lines but we can't find these axes using this method. More importantly, $x$ and $y$ turns out to be uncorrelated. 
A: 4-fold rotational symmetry
The distribution looks like
n = 10^3
x = rnorm(n,0,1)
y = x*(rbinom(n,1,0.5)*2-1)

plot(x,y, pch=21,
     col=rgb(0,0,0,0.1),bg=rgb(0,0,0,0.1))


It has a 4-fold rotational symmetry which means that a quarter rotation leaves the distribution unchanged. That is $p(x,y) = p(-y,x)$.
For any bivariate distribution with a 4-fold rotational symmetry, there are no unique principle components.
Motivation
For any rotational transformation with an angle by $\theta$ we have that the sum of the variance of the transformed coordinates is invariant.
$$ \text{VAR}(X_\theta) + \text{VAR}(Y_\theta) = \text{VAR}(X) + \text{VAR}(Y) \tag{1}$$
Because of the symmetry we have that $$\text{VAR}(X_{\theta + 0.5 k \pi }) = \text{VAR}(X_{\theta}) \\
\text{VAR}(Y_{\theta + 0.5 k \pi }) = \text{VAR}(Y_{\theta}) \tag{2}$$
Because the axes are 90 degrees rotated from each other we also have
$$\text{VAR}(X_{\theta+ 0.5 \pi}) = \text{VAR}(Y_{\theta}) \\
\text{VAR}(Y_{\theta-0.5 \pi}) = \text{VAR}(X) \tag{3}$$
combining equations 2 and 3 it follows that for any rotation $\theta$
$$\text{VAR}(X_{\theta}) = \text{VAR}(Y_{\theta})$$
and because of equation 1, the variance must sum up to a constant, we have that after any rotation we will have that the variance of a coordinate is equal to halve the total variance.
The consequence is that we can not choose any unique direction for the components
