A Gaussian with support on d-dimensional space will have a mean of shape(d,) and a covariance matrix of shape (d, d), which is always positive semidefinite (symmetric and all eigenvalues non-negative).

The covariance matrix can be:

1. Scalar:
   • diagonals are equal
   • off-diagonals are zero
   • represented by ()-shaped tensor
2. Diagonal:
   • diagonals can be unequal
   • off-diagonals are zero
   • represented by (d,)-shaped tensor
3. Full:
   • diagonals can be unequal
   • off-diagonals can be non-zero, but must be symmetric
   • typically represented by (d,d)-shaped tensor, but uniquely by a (d(d+1)/2,) tensor of the lower/upper triangular elements.

We can discuss here what the best approach to handling these would be. E.g., avoiding Cholesky decomposition in MultivariateNormal in cases 1 and 2 would be good.