Calculating covariance matrices with Python

I wanted to use lower level linear algebra constructs in Python to calculate a covariance matrix.

pythonmathematics

Covariance matrices quantify the relationship between two or more variables. This post solves the problem using fundamental numpy linear algebra functions rather than relying on np.cov. The notation below follows Stat Trek’s guide to covariance matrices.

To create a variance-covariance matrix from , which is an matrix holding ordered sets of raw data, transform the raw data into deviation scores in matrix .

Where:

  • is a vector of ones.
  • is the transpose of .
  • is an matrix of deviation scores: .
  • is an matrix of raw data.

We then need to calculate which is the deviation sums of squares and cross products matrix for . Then we divide each term in the deviation sums of squares and cross products matrix by to create the variance-covariance matrix .

Where:

  • is a variance-covariance matrix.
  • is the transpose of .
  • is the deviation sums of squares and cross-products matrix.
  • n is the number of measurements/rows/scores in each column of the data matrix .

Now that we understand the mathematics, let’s look at how to implement it in Python. NumPy has the function np.cov that calculates the covariance matrix. Note that the bias parameter in np.cov normalises by if true rather than the default setting of .

X = np.array([[1, 2], [3, 6]])
np.cov(X.T, bias=True)
array([[1., 2.],
       [2., 4.]])

To verify the implementation against first principles, I’ll implement the calculation directly.

x = X - np.dot(np.ones(X.T.shape), X) * 1/X.shape[0]
var_co_var_x = np.dot(x.T, x) * 1/X.shape[0]
var_co_var_x
array([[1., 2.],
       [2., 4.]])

Both implementations produce identical results.