Assigning eigenvectors of a covariance matrix to the variables it was generated from

Question

Preface: This is the follow-up question according to the insights I got from MvG on my earlier post.

Think of an ellipsoid in the n-dimensional space defined by $$ell: (x-\mu)'A(x-\mu)=1.$$ Then one could calculate $A$ from the covariance matrix of a data set, such that $$A = cov(dataset)^{-1}.$$ Say, each of the $N$ rows of $dataset$ holds an observation of 3 variables, namely $x_{1}, x_{2},$ and $x_{3}$. Then $dataset$ is a matrix with dimension $N\times3$, and $A$ is then a matrix of $3\times3$.

Refering to the eigenvectors of $A$ by $eVecs$ and the corresponding eigenvalues by $eVals$, the eigenvector with the largest eigenvalue will correspond to the direction along which $dataset$ has the maximum variance see here for details.

What I would like to know: I would be happy to get advise how to sort $eVecs$ and $eVals$ to correspond to the variables $x_{1}, x_{2},$ and $x_{3}$. See the example below for better understanding, there, the variables $X$ are drawn from uniform distributions where the upper limit of the distribution is increases from 1 over 10 to 100. The eigenvalues are then sorted by software R with respect to their magnitude (no option to avoid sorting). Nevertheless, the first eigenvalue (12.7) would correspond to the last variable (limit = 100). Is there an algorithm to resort $eVals$ in less clear situations? Is sorting them according to the variances of the $X$ enough ?(I guess not since "direction" does not mean "variable" when the variables are correlated - right?)

  > dataset = cbind(runif(1000,0,1),runif(1000,0,10),runif(1000,0,100))
  > A = solve(cov(dataset))
  > eigen(A)$values
  > 12.685907528  0.129176140  0.001206275

What I need to do: I want to draw uniformly from $ell$, for which $eVecs$ and $eVals$ are needed. Here is a matlab code for that. Nevertheless (and having no idea of matlab), this algorithm seems to produce unreasonable results in R, and I read (can't verify) that matlab does not sort $eVecs$ in the same way R does. Tests look to me - no definite conclusion - that sorting of the $eVecs$ produces the problems in R.

Thanks a lot for any help, Johannes

Answer 1

I'll concentrate on this part of your question:

What I need to do: I want to draw uniformly from ell, for which eVecs and eVals are needed.

The function eigen basically provides you with a diagonalization of the matrix. So you get

$$A=Q\cdot \Lambda\cdot Q^T$$

where $\Lambda$ is a diagonal matrix of eigenvalues, and the columns of $Q$ are the eigenvectors. Since the eigenvectors are of unit length, these columns will form an orthonormal system. Hence $Q^{-1}=Q^T$ which simplifies things a lot.

Assume that the ellipse is centered around the origin, i.e. your $\mu$ is zero. Plugging a vector $x$ into your equation becomes

$$x^T\cdot A\cdot x = x^T\cdot Q\cdot\Lambda\cdot Q^T\cdot x = (Q^T\cdot x)^T\cdot\Lambda\cdot(Q^T\cdot x) = y^T\cdot\Lambda\cdot y = 1$$

So instead of plugging your original vector $x$, from your original coordinate system, into the original matrix $A$, you can as well plug the transformed vecotr $y=Q^T\cdot x$ into the diagonal matrix $\Lambda$, which is a lot simpler. Since $\Lambda$ is diagonal, it describes an ellipsoid which is aligned with the coordinate axes. The eigenvalues represent reciprocal squared radii:

$$y^T\cdot\Lambda\cdot y = \lambda_1y_1^2 + \lambda_2y_2^2 + \dots + \lambda_ny_n^2 = \left(\frac{y_1}{r_1}\right)^2 + \left(\frac{y_2}{r_2}\right)^2 + \dots + \left(\frac{y_n}{r_n}\right)^2 = 1 \\ r_i = \frac1{\sqrt{\lambda_i}} $$

So here is what I'd do: Start by sampling uniformly from the shere. By the way, do you mean surface or interior of the sphere? For the Surface of $S^2$, the normal sphere as embedded in 3d space, I once heard that combining three normally distributed coordinates and normalizing the result to unit length would yield uniform distribution on the surface. I don't have a reference for this just now.

Next, scale each coordinate separately by the corresponding radius $r_i$ as obtained above. The result will be sampled from the ellipsoid in the new transformed coordinate system. If you were sampling from the interior of ell, then I'd say the sampling would still be unique. If you are sampling from the surface, then probably not. Then transform that result back into the old coordinate system via the transormation $x=Q\cdot y$. Since this is essentially a rotation, the sampling stays uniform if it was before, i.e. for the interior. If your $\mu$ is non-zero, then you'd add that as a last step.

how to sort eVecs and eVals to correspond to the variables $x_1,x_2$ and $x_3$

There is no one-to-one correspondence between eigenvectors and data variables. Instead, each eigenvector represents a linear combination of variables, some with positive and some with negative coefficient. This is a bit akin to how the main axes of an ellipsoid don't exactly map to coordinate axes if that ellipsoid is rotates somehow. I elaborated this point in mor detail in this post, which I also mentioned in a comment.