Let's have spatio-temporal dataset ($Y \in \mathbb{R}^{L \times T}$). Where $L$ stands for spatial grid points and $T$ for time. Now let's say that the noise of the system follows a multivatiate gaussian distribution. $y^\ell_t = \epsilon^\ell_t$, where $\epsilon^\ell_t \sim \mathcal{N}(0, \Sigma)$. Standing $\Sigma$ for a the covariance matrix $L \times L$.
Now, we have three different blocs (close spatial points) where the spatial points of each bloc covary together, but independently of the others bloc. Then, our $\Sigma$ is a diagonal block matrix. For example
import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
Sigma = np.zeros((8*30,30*8))
Sigma[:40, :40] = 0.3 # First bloc
Sigma[88:136, 88:136] = 0.3 # Second block
Sigma[200:, 200:] = 0.3
np.fill_diagonal(Sigma,1)
plt.imshow(Sigma)
Now we create a random data from that $\Sigma$ covariance matrix. And perform pca on top of it. I would expect the first three components to stand for each region where noise covary. For the shake of interpretability, we will assume that our spatio-temporal model lives in a $8 \times 30$ space, that is $L = 240$. In this spatio temporal model there are three blocs of covariant points corresponding to the contiguous grid-points 0:44, 88:136 and 200:240. According to $\Sigma$. Then,
Y = np.random.multivariate_normal(np.zeros(8*30), Sigma, size=50000)
pca = PCA(n_components=3)
pca.fit(Y)
for i in range(3):
# Reshaping the components to a spatial representation allows us to plot
# the components as the contribution of each grid point to each component.
# The spatial dimensions are 30 x 8.
plt.imshow(pca.components_[i, :].reshape((30,8)))
plt.colorbar()
plt.show()
plt.close()
The fact is, that I see that they are correlated. And I'm not really sure why, given that I would expect that the principal components of PCA to be orthogonal between them. I also plot the correlation matrix of the first three components.
components = pca.components_
corr_2 = np.corrcoef(components)
plt.imshow(corr_2)
plt.colorbar()
Shouldn't the elements outside of the diagonal of that matrix be 0?
I've also try Varimax rotation, where the three components are now independent.
I guess I am missing some fundamental understanding of what means that the components are orthogonal. Could someone bring light on why the correlation matrix of the 3 comnponents has non-zeros outside of the diagonal?