Could you please also add a Cholesky example?

Two comments on the interface:

1. It would be clearer to provide two alternative arguments sigma=None and chol_sigma=None and derive one from the other, rather than provide a single sigma and a flag is_cholesky=False. This is also more consistent with other distributions classes like Bernoulli and Categorical which can take either probs=None or logits=None and derive one from the other internally. I suppose the more common should be first, but I don't know which is more common.
2. Going forward we're trying to maximize compatibility with Tensorflow Distributions, so it would be preferable to rename mu -> loc, sigma -> covariance_matrix, and sigma_cholesky -> scale_tril. For reference, see https://www.tensorflow.org/versions/r1.0/api_docs/python/tf/contrib/distributions/MultivariateNormalCholesky and https://www.tensorflow.org/versions/r1.0/api_docs/python/tf/contrib/distributions/MultivariateNormalFull

This normalized flag is nice to have. Right now we don't have the plumbing to support it as an argument to batch_log_pdf(), but we can support it as an argument to the __init__ method. Could you move it there and make an instance attribute .normalized?

I guess the gradients will be incorrect if sigma_cholesky.requires_grad=True and normalized=False. Do you think it would be reasonable to either (1) raise a NotImplementedError in this case (to ensure users set sigma_cholesky.requires_grad=False) or (2) compute a numerically-stable quantity whose value is zero but whose gradient agrees with the gradient of the normalization factor (this is a research question 😉 )?

I see the problem there, of course the gradient of the normalization coefficient is not taken into account at the moment. I am however not completely sure how much of a problem this is, as I would presume that the normalization might not be of real importance for many applications. I must admit that I have not thought too much about this problem until now, but I will look into it. Approach (2) seems better to me though, as (1) would disable differentiation for those who don't care about normalization too.

If (2) is difficult, I think it's fine to merge (1) for now and implement (2) later once we know how. I'd like to merge your PR soon since @fehiepsi is already building a Gaussian Process tutorial #650 and it would be nice if he could base that on your work.

Regarding (1) I have added a warning about the incorrect gradients to batch_log_pdf. This seems to be the best option to me, as not being able to differentiate if normalized is set to false would seriously limit the usefulness of this feature in my opinion.

Approach (2) might also not be as difficult as it seems. The gradient of the determinant of X apparently is 2*X.inverse() - torch.diag(torch.diag(X.inverse())), a rather nice expression (see (141) in http://www.math.uwaterloo.ca/~hwolkowi//matrixcookbook.pdf) especially since the inverse is already being computed. I am however not sure how to inject that gradient into autograd, do you have any suggestions as to how that might be achieved efficiently?

I think the easiest way is to create an autograd function:

from torch.autograd import Function

class _NonnormalizedDeterminate(Function):
    @staticmethod
    def forward(ctx, matrix):
        ctx.save_for_backward(matrix)
        return matrix.new([1.0])  # A bogus value.

    @staticmethod
    def backward(ctx, grad_output):
        matrix, = ctx.saved_variables
        inv = matrix.inverse()
        grad = 2 * inv - torch.diag(torch.diag(inv))
        return grad_output * grad

You can then use this as

normalization_factor = torch.log(_NonnormalizedDeterminant.apply(sef.sigma)) + ...

Also feel free to open an issue on this so we don't lose this extended discussion.

nit: You can decorate with @copy_docs_from(Distribution) and then omit trivial method docstrings. See other distributions for usage.

I think the easiest interface to use would be to require one of the args to be specified (as with Bernoulli):

def __init__(self, loc, covariance_matrix=None, scale_tril=None, ...):
    if covariance_matrix is None and scale_tril is None:
        raise ValueError('At least one of covariance_matrix or scale_tril must be specified')
    if scale_tril is None:
        covariance_matrix = ...
    elif covariance_matrix is None:
        scale_tril = ...

Okay, I am going to update that. That will make using the Nones in the test obsolete too.

Could you omit the None arguments below:

{'loc': [2.0, 1.0], 'covariance_matrix': [[1.0, 0.5], [0.5, 1.0]],
 'test_data': [[2.0, 1.0], [9.0, 3.4]]},
{'loc': [2.0, 1.0], 'scale_tril': [[1.0, 0.5], [0, 3900231685776981/4503599627370496]],
 'test_data': [[2.0, 1.0], [9.0, 3.4]]},

I like your potri_compat() function above. Could you factor this line out as a similar potrf_compat()?

No problem, I am going to include that in the next commit.

As I understand, it's generally a good idea to avoid computing an inverse when solving an equation: if you have a cholesky you should use a Cholesky solver, or if you have a matrix you should use a general matrix solver. From what you know of PyTorch, is this not yet possible?

It looks like you could do this with torch.gesv(): http://pytorch.org/docs/master/torch.html#torch.gesv

I was thinking that it is easier to use torch.inverse() because torch.gesv() is not batch solver.

Yeah, that is really a problem. It looks like using torch.gesv() only brings performance gains if the batch size is considerably smaller than that of the covariance matrix.

nit: Could you ensure that docstrings are at most 80 characters wide (here and in .sample() below)? This is important for users of Jupyter Notebooks so that help(MultivariateNormal) prints correctly.

Your review comment was helpful in understanding how we'll update this code in the future. Could you insert that comment here as a # TODO It may be useful to switch between matrix_inverse_compat() and linear_solve_compat() ...?

I'm inclined to raise an error rather than warn because I do not understand what is being computed in this case. Do you have a clear understanding of what is being computed and whether there are any statistical applications for computing this new quantity?

Of course the gradients would only be a kind of rather bad estimate of the true gradient, but it looked to me like it was good enough to do some things like fitting the parameters of a kernel for a gaussian process. But I did investigate that very thoroughly so I will change the code to raise an error for the moment. It looks like (2) can be implemented very soon anyway.

Thanks, raising an error seems like the right thing to do, and you've convinced me that a correct solution should be fairly cheap.