For this question, I am using the following Wiki definition of Matrix whitening:
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?