Positive definiteness error due to import of pyro

Hello,

I am trying to create a Multivariate normal distribution from a mean and covariance matrix obtained with SWAG.

When using pyro.distributions to create the distribution I got a positive definiteness error such as the following:

ValueError: Expected parameter covariance_matrix (Tensor of shape (2570, 2570)) of distribution MultivariateNormal(loc: torch.Size([2570]), covariance_matrix: torch.Size([2570, 2570])) to satisfy the constraint PositiveDefinite(), but found invalid values:
tensor([[ 4.8379e-01, -1.0609e-02,  2.0726e-01,  ..., -3.8059e+00,
          4.2559e+00, -2.2484e+00],
        [-1.0609e-02,  1.2820e-03, -1.3205e-02,  ...,  1.9880e-01,
         -1.9384e-01,  1.0893e-01],
        [ 2.0726e-01, -1.3205e-02,  6.3967e-01,  ..., -3.7263e+00,
          3.7177e+00, -2.1286e+00],
        ...,
        [-3.8059e+00,  1.9880e-01, -3.7263e+00,  ...,  1.4046e+02,
         -7.0172e+01,  3.8460e+01],
        [ 4.2559e+00, -1.9384e-01,  3.7177e+00,  ..., -7.0172e+01,
          1.5973e+02, -4.2080e+01],
        [-2.2484e+00,  1.0893e-01, -2.1286e+00,  ...,  3.8460e+01,
         -4.2080e+01,  4.6706e+01]], device='cuda:0')

I have checked the constraints for positive definitness and all of the them are fulfilled by my covariance matrix.
Here is a sample of the checks I performed:

import numpy as np
cov = torch.load('cov_mat')
cov = cov.cpu().numpy()
np.allclose(cov,cov.T)   # Returned true
np.linalg.eig(cov)         #Returned only positive eigenvalues
np.linalg.cholesky(cov)       #Was feasible

I thus tried to use torch.distributions.MultivariateNormal instead of the corresponding pyro implementation, which from what I have read is simply a wrapper class of the corresponding pytorch class.

Doing so, I had absolutely no issue creating the distribution and could sample from it without any issues.

I then tried to import pyro and then redo the same with pyro.distributions.MultivariateNormal, but it failed again.

In short after a bit of testing, I noticed that simply importing pyro without using it would make torch.distributions.MultivariateNormal throw the same error as pyro.distributions.MultivariateNormal.

So the following code only works with a fresh start of python if “import pyro” is commented out:

import torch
import pyro

cov = torch.load('/cov_mat')
dis = torch.distributions.MultivariateNormal(torch.zeros(len(cov)).to('cuda'), cov)
print(dis.sample())

I’ve tried running the following tests for positive definiteness which are the ones run in pytorch:

print(torch.distributions.constraints._Symmetric().check(cov))
print(torch.linalg.cholesky_ex(cov).info.eq(0))

both return True, meaning that my covariance matrix is indeed positive definite.

Does anyone know why simply importing pyro without using it makes the torch code fail? And why the pyro implementation does not work, while the pytorch one does (without the import of pyro)?

I am thinking that precision might be the issue. So could it be that a simple import of pyro somehow changes the precision?

This is strange. Could you please create a github issue with small reproducible code?

Sure, but how should I go about transferring the matrix I use as a covariance matrix? (It has dimensions 2570 x 2570). I tried taking very small subsections of it to easily be able to write them in some code, but then it does not throw the error anymore.

Maybe you can try to enable a debugger and check what happens when around the error message?

I have just done so and have found the reason for the difference/error.

I have been running the following with the debugger:

import torch
import pyro

cov = torch.load('cov_mat')
dis = torch.distributions.MultivariateNormal(torch.zeros(len(cov)), cov)
print(dis.sample())

with a break point on the line with “dis = torch …”.

When commenting out the import of pyro and running the script, the PositiveDefinite() function from torch.distributions.constraints is being called to check for positive definiteness.
The function there checks if the matrix is symmetric and if a cholesky factorization is feasible.

Here is the function:

class _PositiveDefinite(_Symmetric):
    """
    Constrain to positive-definite matrices.
    """
    def check(self, value):
        sym_check = super().check(value)
        if not sym_check.all():
            return sym_check
        return torch.linalg.cholesky_ex(value).info.eq(0)

When importing pyro on the other hand, the function PositiveDefinite() from pyro.distributions.torch_patch gets called instead. It checks if all eigenvalues are positive and is defined as follows:

@patch_dependency("torch.distributions.constraints._PositiveDefinite.check")
def _PositiveDefinite_check(self, value):
    matrix_shape = value.shape[-2:]
    batch_shape = value.shape[:-2]
    flattened_value = value.reshape((-1,) + matrix_shape)
    return torch.stack(
        [torch.linalg.eigvalsh(v)[:1] > 0.0 for v in flattened_value]
    ).view(batch_shape)

In the first scenario, the matrix passes both the symmetry test and the cholesky factorization test, while in the second it does not pass the test because some of the eigenvalues are negative.

I am using version 1.11.0 of pytorch but have checked and the function is still unchanged in the newest version of pytorch on github.

Since my matrix does contain some elements with small positive values in its diagonal, I am guessing that a higher computational error with one method than the other could maybe explain why one positive definiteness check fails and not the other. Or maybe a higher tolerance for one than the other could also be the reason?

From the algorithm I know for eigenvalue computations, they should be computed recursively. Since my matrix is rather big (2570 x 2570) I am assuming that the propagation of computational error is accumulating to a significant degree and since some positive values in the diagonal of my covariance matrix are rather small, that probably leads to a big error on the last eigenvalues computed. Problem that maybe doesn’t arise when doing a cholesky factorization? (I haven’t looked into the algorithms used for it). That would maybe also explain why I did not get the error when only taking subsections of the covariance matrix.

a possible work around might be to initialize the MultivariateNormal with
scale_tril = torch.linalg.cholesky(...) instead of covariance_matrix=cov

Thanks for debugging, really helpful for me to understand the issue. Do you want to submit a PR to remove the patch? It was used to address some upstream pytorch issue back then, but unnecessary now.

Thanks for the tip that does indeed solve it :slight_smile:

You are welcome. Thanks for helping on the forum and helping develop Pyro :slight_smile: .
I created a pull request, it should be ready for review.