Trouble in calculating the covariance matrix

I'm trying to calculate the covariance matrix for a dummy dataset using the following formula, but it's not matching with the actual result.

Let's say the dummy dataset contains three features, #rooms, sqft and #crimes. Each column is a feature vector, and we have 5 data points. I'm creating this dataset using the following code:

matrix = np.array([[2., 3., 5., 1., 4.],
[500., 700., 1800., 300., 1200.],
[2., 1., 2., 3., 2.]])


Let's normalize the data, so the mean becomes zero.

D = matrix.shape[0]
for row in range(D):
mean, stddev = np.mean(matrix[row,:]), np.std(matrix[row,:])
matrix[row,:] = (matrix[row,:] - mean)/stddev


Now, I can write a naive covariance calculator that looks at all possible pairs of features, and that works perfectly.

def cov_naive(X):
"""Compute the covariance for a dataset of size (D,N)
where D is the dimension and N is the number of data points"""
D, N = X.shape
covariance = np.zeros((D, D))

for i in range(D):
for j in range(i, D):
x = X[i, :]
y = X[j, :]
sum_xy = np.dot(x, y) / N
if i == j:
covariance[i, j] = sum_xy
else:
covariance[i, j] = covariance[j, i] = sum_xy
return covariance


But, if I try to implement the formula mentioned in the beginning, the result is incorrect. The method I am trying out is as follows:

def cov_naive_2(X):
"""Compute the covariance for a dataset of size (D,N)
where D is the dimension and N is the number of data points"""
D, N = X.shape
covariance = np.zeros((D, D))

for i in range(N):
x = X[:, i]
covariance += x @ x.T
return covariance / N


What am I doing wrong here?

Expected output:

array([[ 1.        ,  0.96833426, -0.4472136 ],
[ 0.96833426,  1.        , -0.23408229],
[-0.4472136 , -0.23408229,  1.        ]])


Actual output from cov_naive_2

array([[3., 3., 3.],
[3., 3., 3.],
[3., 3., 3.]])

• Can you add your expected and yielded outputs? – Romain Reboulleau Dec 14 '19 at 9:33
• @RomainReboulleau added. Please check. – Bitswazsky Dec 14 '19 at 9:41

Your approach is mathematically right, your problem comes from the fact that numpy matrix multiplication defaults to inner product when providing vectors, independently from them being transposed to row or column.

Try modifying the penultimate line like this to force outer product.

def cov_naive_2(X):
"""Compute the covariance for a dataset of size (D,N)
where D is the dimension and N is the number of data points"""
D, N = X.shape
covariance = np.zeros((D, D))

for i in range(N):
x = X[:, i]
covariance += np.outer(x, x)  # <---- here
return covariance / N
$$$$
`