Change broadcasting for ~ macro

[trappmartin]

Currently the ~ macro automatically uses broadcasting if the left-hand side is a vector and the right-hand side a univariate distribution. It would be more convenient for testing and for the user to change it to an explicit syntax, i.e..~, which would also better match the standard Julia syntax.

@yebai What happens in the case of multivariate distributions?

[yebai]

I think ~ will check whether the dimensionality of left-side is a matrix or a vector of vectors. If so, it will treat left-side as a vector of draws from multivariate distributions.

[willtebbutt]

@trappmartin I agree with your point regarding Julia-ness, but we're quite closely tied in to the Distributions.jl conventions, and they make very specific assumptions regarding what multiple samples from a distribution look like depending upon whether it's a UnivariateDistribution, MultivarateDistribution, or a MatrixDistribution (see the Distributions.jl documentation for more info). Consequently, I think it actually makes sense not to follows the standard Julia conventions here, but rather to vectorise in the traditional sense and follow the Distributions.jl conventions.

To provide a concrete example of where we would get in trouble using a .~ notation, the Distributions convention for multiple samples from a MultivariateDistribution is to produce a matrix whose columns are observations. What, then, does the broadcasting .~ notation mean? One the one hand our intent would probably be to broadcast observation-wise, but broadcasting is already defined for AbstractMatrix in an element-wise fashion, so it would necessarily involve introducing some non-standard behaviour, which I believe would be a mistake.

[trappmartin]

Thanks, @willtebbutt! The convention of the Distributions package is quite unfortunate. I understand that with this convention, we cannot generically broadcast.

Let's keep this issue open and see if we can find a more transparent solution to this. I feel that the implicit broadcasting solution we currently have is sub-optimal.

[willtebbutt]

Agreed regarding the convention. FWIW, I had the same problem in Stheno.jl with regards to data sets. The inputs to a GP (or any other model for that matter) might be a vector of scalars, a matrix of for which each row or column should be interpreted as an observation, etc.

My solution is to require that all data sets subtype AbstractVector and define indexing behaviour. This way the semantics are unambiguous, but you can still have the data stored in whatever manner you like under the hood, and just use multiple dispatch in the usual way to specialise behaviour for different storage types. Case in point, I have a type called ColsAreObs (read: Columns are observations):

struct ColsAreObs{T, TX<:AbstractMatrix{T}} <: AbstractVector{Vector{T}}
    X::TX
    ColsAreObs(X::TX) where {T, TX<:AbstractMatrix{T}} = new{T, TX}(X)
end

designed exactly to handle the situation described above, with the interpretation is that the nth column of X is the nth observation. I then specialise a number of implementations for this type.

With this design you lose some conventions that people are used to seeing, but you gain consistency in terms of the interface for all of the data types without sacrificing any performance (provided that you specialise correctly)

[yebai]

Related discussion from Soss.jl:

cscherrer/Soss.jl#14

[xukai92]

Another reason to stick to Distributions' convention is that the vectorised version of logpdf is faster than that of the broadcasting version for some distributions, e.g. MvNormal

julia> using Distributions, BenchmarkTools

julia> d = MvNormal(zeros(2), ones(2))
DiagNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 0.0; 0.0 1.0]
)

julia> rvs = rand(d, 1000);

julia> @benchmark sum(logpdf(d, rvs))
BenchmarkTools.Trial: 
  memory estimate:  23.70 KiB
  allocs estimate:  3
  --------------
  minimum time:     4.171 μs (0.00% GC)
  median time:      5.322 μs (0.00% GC)
  mean time:        7.590 μs (22.46% GC)
  maximum time:     5.710 ms (99.87% GC)
  --------------
  samples:          10000
  evals/sample:     7

julia> rvs_col = [rvs[:,i] for i = 1:1000];

julia> @benchmark sum(logpdf.(Ref(d), rvs_col))
BenchmarkTools.Trial: 
  memory estimate:  102.63 KiB
  allocs estimate:  1022
  --------------
  minimum time:     53.432 μs (0.00% GC)
  median time:      57.578 μs (0.00% GC)
  mean time:        70.490 μs (9.19% GC)
  maximum time:     34.866 ms (99.74% GC)
  --------------
  samples:          10000
  evals/sample:     1

Ideally we'd like to call these vectorised version of logpdf in Turing instead of the broadcast version.

Although this is not very consistent in Distributions are this vectorized version for univariate version is deprecated. If you try this for Normal you will see "Warning: logpdf(d::UnivariateDistribution, X::AbstractArray) is deprecated, use logpdf.(d, X) instead."

[mkborregaard]

Is the Distributions convention set in stone?

[willtebbutt]

Is the Distributions convention set in stone?

I'm not sure what the Distributions maintainers' position is on this. FWIW my personal feeling is that we should consider adopting the convention I discussed above in Turing for consistency's sake. If we do that we could (either implicitly or explicitly) require that

1. The type of the lhs of a ~ is distributed over by the rhs, and
2. The lhs of a .~ is compatible with broadcasting (whatever that looks like)

edit: we should find out what the Distributions position is on changing their conventions for multiple sampling.