The multiplication between a matrix and its transposed is not symmetric and psd

I want to compute the covariance matrix of a batch of feature vectors (replaced by torch.randn here) using the code below:

import torch

a = torch.randn(64, 512)
cov_mat = (a.T @ a) / (a.size(0) - 1)

The produced cov_mat is symmetric and positive semidefinite. So I changed the clause above by adding cov_mat = (cov_mat + cov_mat.T) / 2

Now the covariance matrix is symmetric, but still not positive semidefinite. I've tried torch.mm but still didn't work. Does it have something to do with the precision?

Here is the specific version:
Python 3.7
PyTorch 1.8.1

I got similar issues, and the conclusion was indeed that it was a precision issue. When I executed your code (without your modification), I observed that cov_mat is almost symmetric, with an error around 1e-6. This is actually expected, as explained here:

> Floating point precision limits the number of digits you can accurately rely on.
torch.float32 has precision between 6 and 7 digits regardless of exponent, i.e. precision issues are expected for values < 1e-6. Similarly torch.float64 will have precision issues for values < 1e-15.

If you really need more precision, you can use torch.float64 instead:

# with the default float32:
a = torch.randn(64, 512)
c = (a.T @ a) / (a.size(0) - 1)
((c - c.T) < 1e-7).all()    >>> tensor(False)
((c - c.T) < 1e-6).all()    >>> tensor(True)

# with float64:
a = torch.randn(64, 512, dtype=torch.float64)
c = (a.T @ a) / (a.size(0) - 1)
((c - c.T) == 0).all()      >>> tensor(True)

While your workaround using cov_mat = (cov_mat + cov_mat.T) / 2 also achieves similar precision using float32 (in the sense that the matrix is exactly symmetric), it will make the matrix "less" psd.

Unfortunately, whatever you do, you will always suffer from precision errors (1e-6 for float32, 1e-15 for float64). As such, we can see that the eigenvalues of the covariance matrix may be very slightly below zero. I tried manually setting the eigenvalues to 0, but there still are errors.

TL;DR My final advice would be:

• if you can afford smaller errors, then just go for float64
• if you need the matrix to be psd to use specific properties (like the fact that x.T @ C @ x > 0) then clamp the resulting value to be above 0
• if you absolutely need the matrix to be psd, you can do: c += 1e-12 * torch.eye(a.size(1)) (I don't recommend it if you don't need it because it modifies the matrix)

Hope that helps you!