Multiplying vectors by the covariance matrix? I thought I knew covariance but I'm starting to think that there's more to it. For example, what happens when you multiply observations by their corresponding covariance matrix? [x1,y1] * cov(x,y). I did a little experiment and am interested by the outcome.
import matplotlib.pyplot as plt
import numpy as np

mean = [0, 0]
cov = [[1, 0.9], [0.9, 1]] 

x, y = np.random.multivariate_normal(mean, cov, 1000).transpose()
plt.scatter(x, y)
plt.show()


new_x = []
new_y = []
for i in range(500):
  vec = np.array([x[i],y[i]])
  trsf = np.matmul(vec.transpose(),cov)
  new_x.append(trsf[0])
  new_y.append(trsf[1])

plt.scatter(new_x,new_y)   


As you can see, every (x,y) pair is projected onto a line w/o any deviation. I'm curious what's actually happening? My guesses are (A) these data points are projected onto the OLS line. Or (B) these data point are projected onto the leading eigenvector of the covariance matrix.
Any thoughts/ideas?
 A: They're not exactly on a line, but yes, they generally follow the main eigenvector because the projected point has the following coordinates (assuming $x$ with dimensions $2\times 1$):$$\Sigma x=\sigma_1u_1u_1^Tx+\sigma_2u_2u_2^Tx=(\sigma_1<u_1,x>)u_1+(\sigma_2<u_2,x>)u_2$$
Here, $u_i$ are eigenvectors, and $\sigma_i$ are eigenvalues. Since $\sigma_1>\sigma_2$ and for the most of the data component in the direction of $x_1$ is larger than the one in $x_2$, the coefficient before $u_1^T$ is typically much larger than the one before $u_2^T$. This causes the points align with the first eigenvector mostly.
If there were some outliers, especially nearly perpendicular to the first eigenvector, they would not be aligned so much with $u_1$'s direction.
import matplotlib.pyplot as plt
import numpy as np

mean = [0, 0]
cov = [[1, 0.9], [0.9, 1]] 

x, y = np.random.multivariate_normal(mean, cov, 1000).transpose()
x = np.hstack([x, -10])
y = np.hstack([y, 10])

plt.scatter(x, y)
plt.show()


new_x = []
new_y = []
for i in range(len(x)):
  vec = np.array([x[i],y[i]])
  trsf = np.matmul(vec.transpose(),cov)
  new_x.append(trsf[0])
  new_y.append(trsf[1])

plt.scatter(new_x,new_y)   


