Whitening transformation does NOT return a unit covariance matrix

Question

For this question, I am using the following Wiki definition of Matrix whitening:

Suppose $X$ is a random (column) vector with non-singular covariance matrix $\Sigma$ and mean 0. Then the transformation $Y=WX$ with a whitening matrix $W$ satisfying the condition $W^TW=\Sigma^{-1}$ yields the whitened random vector $Y$ with unit diagonal covariance.

From the definition, I expect the covariance matrix of $Y$ to be the identity matrix. However, this is far from the truth!

Here is the reproduction:

import numpy as np
# random matrix
dim1 = 512 # dimentionality_of_features
dim2 = 100 # no_of_samples

X = np.random.rand(dim1, dim2)
# centering to have mean 0
X = X - np.mean(X, axis=1, keepdims=True)

# covariance of X
Xcov = np.dot(X, X.T) / (X.shape[1] - 1)

# SVD decomposition
# Eigenvecors and eigenvalues
Ec, wc, _ = np.linalg.svd(Xcov)
# get only the first positive ones (for numerical stability)
k_c = (wc > 1e-5).sum()
# Diagonal Matrix of eigenvalues
Dc = np.diag((wc[:k_c]+1e-6)**-0.5)
# E D ET should be the whitening matrix
W = Ec[:,:k_c].dot(Dc).dot(Ec[:,:k_c].T)

# SVD decomposition End

Y = W.dot(X)
# Now apply the same to the whitened X
Ycov = np.dot(Y, Y.T) / (Y.shape[1] - 1)
print(Ycov)

>> [[ 0.19935189 -0.00740203 -0.00152036 ...  0.00133161 -0.03035149
      0.02638468]  ...

It seems that it won't give me a unit diagonal matrix, unless, dim2 >> dim1.

If I take dim2=1 then I get a vector (although in the example I get an error due to division by 0), and by the Wikis definition, it is incorrect?

Answer 1

As the comments notice, you may have some confusion in your head between covariance and sample covariance. However, that's not what causes your error.

First of all, forget about getting the covariance of $Y$ to be the identity matrix ("unit diagonal"). For rank reasons, it's possible only if $XX^T$ is full rank, and that is not the case because you are in an unusual setup in which the number of samples is smaller than the dimensionality: the rank of $XX^T$ is at most dim2, and so is the rank of $YY^T=(WX)(WX)^T$.

I don't know where you found that algorithm for the whitening matrix, but I find it very confusing. There are some steps in it to take care of rank deficiency (that whole k_c business), but I'm not sure if they are actually helpful; it seems to me that whenever k_c is not equal to dim1 you will get a wrong result: if I am not mistaken, the formulas simplify to $YY^T/(n-1) = EDE^T$, where $E$ is orthogonal, and $D$ is a diagonal that has k_c ones and the rest zeros; when $D$ contains zeros, this matrix is not diagonal (unless $E$ is very special).

However, even if $XX^T$ does not have full rank, you can still get $Y$ to have diagonal covariance equal to $D$, with ones and zeros on the diagonal. To do this, you just use a similar formula but without the first factor $E$:

Dc = np.zeros_like(Ec)
Dc[:k_c,:k_c] = np.diag((wc[:k_c])**-0.5)
W = Dc.dot(Ec.T)

I have also removed the +1e-6, which seems superfluous (and even harmful) because you are already truncating small eigenvalues anyway.

As a final note, as someone who works in numerical linear algebra, I can't avoid pointing out that svd($XX^T$) is not the most numerically stable way to compute the singular values and vectors of $X$, but again that's not the main issue here. Use Ec, wc, _ = np.linalg.svd(X); wc = wc / math.sqrt(X.shape[1]-1) instead.