# Mysterious failure to generate independent set of random variables

After getting my answer in this:

Making a well conditioned orthonormal basis

I am running into a problem which I do not understand. I have n dependent gaussian random variables that are related by a known n x n covariance matrix C. The solution I was given was to take the positive definite square root of the inverse of C - which means I put the n independent random variables (call them the vector x) through a linear transformation (D^-.5*Q)*x, where D is a square diagonal matrix of the eigenvalues (of C) and Q's rows are the corresponding eigenvectors.

The trouble is, the covariance matrix isn't necessarily positive definite under certain values of the parameters of the problem. I don't want a complex matrix, taking the negative eigenvalues to the -.5 power, so my idea was, I keep just the positive eigenvalues. So I look at D, and if it's a 100 x 100 matrix but only 67 of the eigenvalues are positive, I throw out the negative ones, and use only a 67 x 67 matrix for D, which I can take to the -.5 power, and then Q is no longer a square matrix but a 67 x 100 matrix, and the 100-length vectors are the normalized eigenvectors of C corresponding to the eigenvalues that I kept. And lo and behold, when I find this matrix A=(D^-.5*Q), and I put x under the linear transformation Ax, so that the covariance matrix of the result SHOULD be a 67 x 67 matrix, for 67 random variables that are actually variance-1 and independent and constructed out of the original 100 random variables that were too highly correlated for their own good, and I calculate what SHOULD be the new covariance matrix ACA', I get that it is very nearly the 67 x 67 entry identity matrix. GREAT.

Unfortunately....

it does not hold up to actual observation. When I actually try this transformed set of random variables, they are not even CLOSE to being independent. But apparently, only when I throw out some negative eigenvalues, or when it's really close to not being positive definite. In other words, if the parameters of the problem are such that all 100 eigenvalues are quite positive, it works - Ax produces 100 random variables that appear, by monte carlo simulations, to truly be independent. But as soon as I throw away any eigenvalues - and I don't even have to do that, as soon as the parameters of the problem put it CLOSE to being at the point where I would be getting negative eigenvalues, suddenly different entries in Ax start exhibiting absurdly high correlations. 0.6, 0.7, as high as 0.9. Even though ACA' is very nearly an identity matrix, albeit of fewer dimensions than the original C. I thought at first that maybe I wasn't doing enough tests, but that thought is gone - I try it with tens of thousands of trials and the correlation of something that was 0.709 before maybe changes to 0.706 with the completely new data. They ARE highly correlated. I thought for a while that it was a result of small numerical errors in dealing with large matrices. But that doesn't explain why it suddenly happens when the parameters of the problems contrive to make the lowest eigenvalue close to crossing below 0.

Am I missing something? Isn't ACA' ALWAYS supposed to be the covariance matrix of Ax if C is the covariance of x? Is it ever not? Because ACA' here is very nearly an identity matrix and Ax's elements are flagrantly dependent. And if not, is there another method of producing as many independent equal variance random variables as possible besides throwing out the negative eigenvalues as I was doing?

Firstly, let $n$ be the number of rows and $p$ be the number of columns (variables) of your data matrix $\mathbf{X}$. You never said how many records $n$ you were using, and only stated, for example, the size of your $100 \times 100$ covariance matrix $\mathbf{C}$.
You need to look at the Marcenko-Pastur (MP) law on the density of eigenvalues of $\mathbf{C}$ for a random White Wishart matrix, which is essentially an $n \times p$ data matrix, $\mathbf{X}$, with independent columns. In other words, a White Wishart matrix is filled with i.i.d. standard normal variates.
By definition, if $p>n$, then there will be $p-n$ zero eigenvalues of $\mathbf{C}$. In addition, as $n\rightarrow p$, there will be near-zero eigenvalues of $\mathbf{C}$.
Now with regard to the MP Law, the value $\gamma=p/n$ determines the ratio of the dimensions to the size of the $\mathbf{X}$ matrix. The MP law provides the lower $(\lambda^-)$ and upper bound $(\lambda^+)$ of the eigenvalues for $\mathbf{C}$ derived from an $n \times p$ White Wishart $\mathbf{X}$ matrix. Therefore, knowing $n$ and $p$ of any data matrix (i.e., your data matrix, $\mathbf{X}$), you can first use the MP law to calculate $\lambda^-$ and $\lambda^+$, and any of your eigenvalues from $\mathbf{C}$ which exceed $\lambda^+$ are considered to be the signal eigenvalues and those below $\lambda^+$ are believed to be in the noise region. (eigenvalues in the interval $(\lambda^-, \lambda^+)$ correspond to eigenvalues of a White Wishart matrix having the same dimensions as your data $\mathbf{X}$ matrix.