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.]])
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)  # &lt;---- here
    return covariance / N
```