Non-orthogonal technique analogous to PCA Suppose I have a 2D point dataset and I want to detect the directions of all the local maxima's of variance in the data, for example:

PCA does not help in this situation as it is an orthogonal decomposition and therefore cannot detect both the lines I indicated in blue, rather its output may look like the one shown by green lines. 
Please recommend any technique which might be suitable for this purpose. Thanks.
 A: There are PCA-like procedures for the so-called "oblique" case. In stat-software like SPSS (and possibly also in its freeware clone) PSPP one finds the equivalently called "oblique rotations", and instances of them named as "oblimin","promax" and something more. If I understand things correctly the software tries to "rectangularize" the factor-loadings by re-calculating their coordinates in an orthogonal, euclidean space (as for instance shown in your picture) into coordinates of a space whose axes are non-orthogonal perhaps with some technique known from multiple regression. Moreover I think this works only iteratively and consumes one or more degrees of freedom in the statistical testing of the model.
See wikipedia for rotation-methods in factor analysis
An article with an example of comparision PCA and oblique rotation
The reference-manual of SPSS (at the IBM-site) for oblique-rotations contains even formulae for the computation.     
[Update] (Upps, sorry, just checked that PSPP does not provide "rotations" of the oblique type)
A: I don't have much experience with it, but Vidal, Ma, and Sastry's Generalized PCA was made for a very similar problem.
A: The other answers have already given some useful hints about techniques you can consider, but nobody seems to have pointed out that your assumption is wrong: the lines shown in blue on your schematic picture are NOT local maxima of the variance.
To see it, notice that the variance in direction $\mathbf{w}$ is given by $\mathbf{w}^\top\mathbf{\Sigma}\mathbf{w}$, where $\mathbf{\Sigma}$ denotes covariance matrix of the data. To find local maxima we need to put the derivative of this expression to zero. As $\mathbf{w}$ is constrained to have unit length, we need to add a term $\lambda(\mathbf{w}^\top\mathbf{w}-1)$ where $\lambda$ is a Lagrange's multiplier. Differentiating, we obtain the following equation: $$ \mathbf{\Sigma}\mathbf{w} - \lambda \mathbf{w} = 0.$$
This means that $\mathbf{w}$ should be an eigenvector of the covariance matrix, i.e. one of the principal vectors. In other words, PCA gives you all local maxima, there are no others.
A: Independent Component Analysis should be able to provide you with s good solution. It is able to decompose non-orthogonal components (like in your case) by assuming that your measurements result from a mixture of statistically independent variables.
There are plenty of good tutorials in Internet, and quiet a few freely available implementations to try out (for example in scikit or MDP).
When does ICA not work?
As other algorithms, ICA is optimal when the assumptions for which it was derived apply. Concretely,


*

*sources are statistically independent

*the independent components are non-Gaussian

*the mixing matrix is invertible


ICA returns an estimation of the mixing matrix and the independent components.
When your sources are Gaussian then ICA cannot find the components. Imagine you have two independent components, $x_{1}$ and $x_{2}$, which are $N(0,I)$. Then,
$$
p(x_{1}, x_{2}) = p(x_{1})p(x_{2}) = \frac{1}{2\pi}\exp \left( -\frac{x_{1}^{2}+x_{2}^{2}}{2} \right) = \frac{1}{2\pi}\exp -\frac{||\mathbf{x}||^{2}}{2}
$$
where $||.||$. is the norm of the two dimensional vector. If they are mixed with an orthogonal transformation (for example a rotation $R$), we have, $||R\mathbf{x}|| = ||\mathbf{x}||$, which means that the probability distribution does not change under the rotation. Hence, ICA cannot find the mixing matrix from the data.
