How to Sort Covariance Matrix with Python

Question

I have a covariance matrix $\Sigma$ that corresponds to measurements taken at times $t$. So $\Sigma_{ij}$ is the covariance between the measurement taken at time $t_i$ and $t_j$. The trouble is that $t$ is that is unsorted. I need to use $\Sigma$ but where it's entries are sorted with time. How would I sort the entries of $\Sigma$ with increasing time using NumPy?

import numpy as np


# Make random values in the shape of a covariance matrix (the values aren't important for this problem)
covmat = np.random.randn(10, 10)
times = np.random.randn(10,1).squeeze()
indices = np.argsort(times)
# Now what do I do with these indices?

Here's a working example done with loops. I think there is probably a more elegent solution:

import numpy as np


A = np.array([[1,  2,  4,  6,  9],
              [2,  3,  5,  8, 11],
              [4,  5,  7, 10, 13],
              [6,  8, 10, 12, 14],
              [9, 11, 13, 14, 15]])

sort_indices = [2, 3, 4, 1, 0]
sorted_A = np.zeros_like(A)
for i in range(A.shape[0]):
    for j in range(A.shape[1]):
        sorted_A[i,j] = A[sort_indices[i], sort_indices[j]]
print(sorted_A)

>[[ 7 10 13  5  4]
 [10 12 14  8  6]
 [13 14 15 11  9]
 [ 5  8 11  3  2]
 [ 4  6  9  2  1]]

Answer 1

The covariance matrix has by definition the entries $\Sigma_{ij} = \operatorname{cov}(x(t_i), x(t_j))$ now, you are interested in the covariance matrix of the reordered time series $\tilde t_k = t_{\pi(k)}$, where $\pi$ is a given permutation (e.g. the permutation that *un*sorts the time series). Then correspondingly, $$ \tilde \Sigma_{ij} = \operatorname{cov}(x(t_{\pi(i)}), x(t_{\pi(j)})) = \Sigma_{\pi(i), \pi(j)}$$

You have 2 options: to permute the rows/columns of a matrix, you can multiply by permutation matrices from the left/right: $$\tilde \Sigma = P^T \Sigma P \iff \Sigma = P\tilde \Sigma P^T$$

where $P$ is the permutation matrix induced by $\pi$. Alternatively, you can use the indices obtained by the p = np.argsort(p) directly to avoid having to do the actual matrix multiplication via indexing S[p][:,p]

short demo:

import numpy as np
N = 5
t_org = np.arange(N)
pi = np.random.permutation(N)
x_org = np.array([np.random.randn(100)*k for k in range(N)])
S_org = np.cov(x_org)
print("Covariance of sorted time steps", S_org, sep="\n")

t_obs = t_org[pi]
x_obs = x_org[pi]
S_obs = np.cov(x_obs)
print("Covariance of unsorted time steps", S_obs, sep="\n")

#%% Using indexing S[p][:,p]
p = np.argsort(t_obs)
print("Reconstruction equals original:", S_obs[p][:,p] == S_org, sep="\n")

#%% Alternative using Permutation matrix
P = np.eye(N)[p]
print("Permutation matrix", P, sep="\n")
print("Reconstruction equals original:", P@[email protected] == S_org, sep="\n")