Dirichlet dimensionality off-by-one

[cscherrer]

As implemented, link is not bijective for Dirichlet distributions:

julia> dist = Dirichlet([1,1])
Dirichlet{Float64}(alpha=[1.0, 1.0])

julia> link(dist, [0.5,0.5])
2-element Array{Float64,1}:
 0.0
 0.0

julia> invlink(dist,[0.0,0.0])
2-element Array{Float64,1}:
 0.5
 0.5

julia> invlink(dist,[0.0,1.0])
2-element Array{Float64,1}:
 0.5
 0.5

Because values must add to one, a Dirichlet with n values only has n-1 degrees of freedom, so the bijection needs to be to ℝⁿ⁻¹. There's an implementation in TransformVariables.jl that seems to work correctly.

It seems likely this is behind the troubles described in TuringLang/Turing.jl#935

[sethaxen]

FWIW, there's another ℝⁿ⁻¹ → ℝⁿ implementation in HMCUtilities.jl.

[torfjelde]

As you've correctly pointed out this transformation is really a bijection between b: Kⁿ⁻¹→ℝⁿ⁻¹ where Kⁿ⁻¹ is a (n-1)-dim unit-simplex but both Kⁿ⁻¹ and ℝⁿ⁻¹ are considered as subsets of ℝⁿ (since a sample from Dirichlet is represented as a n-dim vector). This also means that the image of b is really only the first (n-1) components, which is why invlink(dist, [0.0, 0.0]) == invlink(dist, [0.0, 1.0]). So when you do project, invlink is only bijective on ℝⁿ⁻¹ ⊂ ℝⁿ i.e. first (n - 1) components. You can do link(dist, x; proj = Val{false}) to not project, but without projection it might not result in a simplex.

The reasoning for the current approach is that it's compatible with the Dirichlet distribution in the sense that you can do logpdf(dirichlet, invlink(x)) for any x ∈ ℝⁿ (though it's not guaranteed to be one-to-one for any x).

@mohamed82008 can probably answer this better:)

And whether or not this is the correct approach; I'm very open to a discussion as to whether or not we should it differently! Especially for the SimplexBijector as that is meant to be "independent" of the distribution and it's currently the only Bijector breaking the assumptions of a Bijector. link and invlink takes the distribution as an argument so I think it's output ought to be compatible with the Dirichlet distribution, as it is now.

[torfjelde]

Also it seems like the difference between the TransformVariables.jl and Bijectors.jl is the following:

https://github.com/tpapp/TransformVariables.jl/blob/44d9dd7add11fa02c2e57f6fbaa24756f1c789a9/src/special_arrays.jl#L120:

    y[end] = stick

vs.

Bijectors.jl/src/Bijectors.jl

Lines 209 to 213 in 45ec713

	@inbounds if proj
	y[K] = zero(T)
	else
	y[K] = one(T) - sum_tmp - x[K]
	end

That is, in Bijectors.jl it's optional to project (proj = true) onto the image of the transformation b: Kⁿ⁻¹ → ℝⁿ⁻¹ ⊆ ℝⁿ, denoted im(b). In the case when y ∈ im(b) ⊆ ℝⁿ it doesn't matter whether or not you "project", i.e. both proj = true and proj = false is going to give to the same result from invlink. But when you're no longer inside im(b), TransformVariables.jl's transform and invlink with proj = false is going to result in the same x ∈ ℝⁿ which is not going to be in the simplex Kⁿ⁻¹ thus not in the support of Dirichlet.

Because of the above, it seems to me that with proj = false and TransformVariables.jl's impls you could potentially run into issues if the HMC trajectory takes a step outside of the im(b), which is why we default to proj = true. Have you tried this in TransformVariables.jl @cscherrer ?

FWIW, there's another ℝⁿ⁻¹ → ℝⁿ implementation in HMCUtilities.jl.

So this is basically what Bijectors.jl is doing with proj = true, but we represent it using a n dim vector with the last element set to zero.

[mohdibntarek]

That Dirichlet bijector is getting a lot of attention here and on Slack and discourse, I am kind of glad! FWIW, I don't think our implementation is wrong in any conceptual way, it's confusing but it is not wrong, and we have good tests to back this claim up. But I do think we have issues with numerical stability where some scenarios can cause log of some tiny negative number. I tried to battle this in the implementation but something seems to have slipped through. @torfjelde's answers are correct.

[penelopeysm]

This was fixed some time ago

julia> using Bijectors

julia> dist = Dirichlet([1,1])
Dirichlet{Int64, Vector{Int64}, Float64}(alpha=[1, 1])

julia> link(dist, [0.5,0.5])
1-element Vector{Float64}:
 0.0