#!/usr/bin/env python
from __future__ import division
from numpy import *
from numpy.random import randn
import numpy.linalg as la

def C(x):
	return sum(x*x.reshape(-1,1,N),axis=-1)/N

N = 1000000
x = randn(2,N)
x = vstack((x,x[0:1]+x[1:]))
Cx = C(x)
w,v = la.eig(Cx)
print 'Cx\n',Cx
print 'Cx evals\n',w
print 'Cx evecs\n',v
M = (v.T)[:2] #drop third column because zero eigenvalue
y = dot(M,x)
Cy = C(y)
print 'Cy\n',Cy
w,v = la.eig(Cy)
print 'Cy evals\n',w
print 'Cy evecs\n',v