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?