Reproduction of NumPy's covariance and correlation matrices from scratch in Python

covariance_and_correlation_matrix.py

import numpy as np

def covariance(X):
    """
    Computes the covariance matrix of X
    :param X: Input data of shape (N, M) where N is the
              number of samples and M is the number of
              features
    :return: Output data of shape (M, M) where Output[i][j]
             is Cov(feature_i, feature_j)
    """
    mean = np.mean(X, axis=0)
    X_zero_centered = X - mean
    return X_zero_centered.T.dot(X_zero_centered)/(len(X) - 1) # Note: there is one less degree of freedom

def correlation(X):
    """
    Computes the correlation matrix of X

    :param X: Input data of shape (N, M) where N is the
              number of samples and M is the number of
              features
    :return: Output data of shape (M, M) where Output[i][j]
             is Corr(feature_i, feature_j)
    """
    mean = np.mean(X, axis=0)
    std = np.std(X, axis=0)
    X_zero_centered = X - mean
    cov = X_zero_centered.T.dot(X_zero_centered)/len(X) # Note: here it's divided by len(X) instead of len(X) - 1
    return (cov/std)/std.reshape((-1, 1))


X =[[0.2, 0.3, 1. ],
    [0. , 0.6, 1. ],
    [0.6, 0. , 3. ],
    [0.2, 0.1, 1. ]]
X = np.array(X)

assert (np.cov(X.T) - covariance(X) < 1e-15).all()
assert (np.corrcoef(X.T) - correlation(X) < 1e-15).all()

A few things to note,

• Zero centering the data has no effect on covariance and correlation since both of them subtract the features by their means. Hence f(X) == f(X - np.mean(X, axis=0)) where f could be both, covariance and correlation.
• Ideally, covariance(X_std) == correlation(X) where X_std = (X-np.mean(X, axis=0))/X.std(axis=0). However, since covariance uses len(X) - 1 instead of len(X) the results are slightly off.