Moving function logdensity into super-lightweight package for abstract densities?

[oschulz]

@tpapp ,I've been thinking about creating a super-lightweight package AbstractDensities or so with a minimal interface for (log-)densities. At the moment, we have so many different approaches in the ecosystem and IMHO not enough compatibility.

Would you consider "donating" the function name logdensity (so have the function defined in AbstractDensities and LogDensityProblems would import/re-export if from there?

I was thinking about a really minimal non-opinionated package with just

• abstract type AbstractDensity, for direct dispatch.
• isdensity(some_density)::Bool, for Holy-trait dispatch
• logdensity(some_density, x)::Real
• logdensity(some_density)(x) == logdensity(some_density, x) (for convenient use in AD, optimization, sampling, etc.)
• LogFuncDensity(f) <: AbstractDensity as a convenient way for users to create a density from a function (e.g. a log-likelihood) that returns the log of a density (logdensity(::LogFuncDensity) would simply return f).

I think it would be great to have a common way to pass around densities that are labeled as such instead of "naked" log-density/likelihood functions. It's kinda awkward/wrong to write likelihood = some_function when that function returns a log-likelihood, but it's also awkward/asymmetric to write log_likelihood = some_function; prior = some_distribution (having "log" in the one and not in the other).

@devmotion , what would be your take on this? Maybe the Turing ecosystem (e.g. AdvancedMH) could support this too? And maybe a lightweight package like this could find a home under JuliaStats?

I would support this in BAT.jl too, of course (I'm looking for a way for users to define log-likelihoods/densities without adding a dependency on a heavy inference package).

[devmotion]

Intuitively, I would think that it is a bit limiting and opinionated if an abstract interface would enforce a specific type structure since multiple inheritance does not exist and packages (such as ones in the Turing ecosystem) already use custom type structures. I think logdensity(custom_type_or_function, x) is sufficient, hence I would assumed that a package could really just define logdensity, similar to what CommonSolve does for solve(!): https://github.com/SciML/CommonSolve.jl/blob/master/src/CommonSolve.jl.

I also think that one has to be careful about the interpretation/meaning of logdensity. E.g., I think it is too restrictive to just use it in a Bayesian setting, i.e., I think it should be fine to have densities where there is no concept of likelihoods or priors.

logdensity(some_density)(x) == logdensity(some_density, x)

Is this actually necessary to be part of the interface? One can always just use Base.Fix1(logdensity, some_density) or define it internally. One could still define a custom type such as LogDensityFunction as

struct LogDensityFunction{F}
    f::F
end

logdensity(f) = LogDensityFunction(f)
logdensity(l::LogDensityFunction, x) = l.f(x)

without a supertype that 1) would make it easy to obtain something that supports logdensity and 2) could be specialized for the type of f.

BTW (a bit unrelated): in the Bayesian setting with logprior, loglikelihood, and logjoint (as it is called in the Turing ecosystem to be clear about what density we refer to - there are algorithms such as elliptical slice sampling that don't operate with the logjoint but only with the logprior and loglikelihood separately), it also seems to be sufficient to also only use an interface that is based on functions but not on types: TuringLang/AbstractMCMC.jl#75 (comment)

[oschulz]

Intuitively, I would think that it is a bit limiting and opinionated if an abstract interface would enforce a specific type structure since multiple inheritance does not exist and packages (such as ones in the Turing ecosystem) already use custom type structures.

Good point. So how about the same as I proposed above, but minus an AbstractDensity type, so just trait-dispatch via isdensity?

I also think that one has to be careful about the interpretation/meaning of logdensity. E.g., I think it is too restrictive to just use it in a Bayesian setting

Oh sure - that was absolutely my intention. My own use cases would be mostly probability densities, but I thought the interface should just be about positive (incl. zero) mathematical densities, not even with an implication that the density can be interpreted as a probability, and nothing at all about Bayes in there.

logdensity(some_density)(x) == logdensity(some_density, x)
Is this actually necessary to be part of the interface?

I think it should be, yes - for example so that if a density is defined via a log-density function, logdensity can just return the wrapped function. It also seems natural to me to have an easy and very readable way to convert a density to a function that returns it's log value. Lastly I've encountered the need to preserve/forward "annotations" (provided via other functions, unrelated to this package), e.g. regarding variable shapes (I would support logdensity in VarableShapes.jl), auto-diff preferences, etc. I wouldn't know how to do certain things in BAT without a logdensity(some_density) that returns a LogDensityOf object (not just an anonymous closure) at part of the API, without committing (without committing type piracy, since BAT wouldn't "own" logdensity).

And we should have LogDensityFunction in the API since the need to wrap a log-density function will be very common, and users should be able to do it without depending on heavier package or rolling their own all the time.

So I think it would be very valuable to have logdensity(some_density) and LogDensityFunction as part of the API, but I agree, an AbstractDensity supertype (or other requirements regarding a supertype) is not necessary. Maybe the package should be called DensityTraits then, instead of AbstractDensities?

[devmotion]

So how about the same as I proposed above, but minus an AbstractDensity type, so just trait-dispatch via isdensity?

This seems reasonable (even though it's not completely clear to me how relevant isdensity would be - I guess you might already have some use case in mind?).

Maybe the package should be called DensityTraits then, instead of AbstractDensities?

Or maybe DensityBase, DensityCore, or DensityInterface?

[oschulz]

even though it's not completely clear to me how relevant isdensity would be

It's necessary for Holy-trait-dispatch (like other is... functions we have in traits-based interfaces, e.g. Tables.isrowtable). And so one can provide reasonable error messages to users when they supply something that's not a density where a density would be expected.

So to enable type-stable code (assuming the result of isdensity can be inferred), like

some_function(density) = some_function_impl(Val(isdensity(density)), density)

some_function_impl(::Val{true}, density) = handle_density_trait_compatible(density)
some_function_impl(::Val{false}, density) = some_other_generic_action_or_so(density)

And to implement sanity checks:

function find_max_density(density, x_init)
    isdensity(density) || throw(ArgumentError("density not compatible with the DensityTraits interface"))
    ...
end

[devmotion]

The technical reason is clear, I just don't recall any instance in the Turing ecosystem or another package where it would be needed. Even in the find_max_density example it seems it is not mandatory and one could instead just rely on the fact that some function call such as logdensity(density, x) would error if it is not a density (and one could also provide a more user-friendly error message by implementing a generic logdensity(density, x) = error("...")).

[oschulz]

True, for find_max_density it wouldn't be necessary (though it couldn't hurt, since one can throw the error earlier on in the program flow). Regarding dispatch, I think it might be a nice design reserve to have isdensity, e.g. in case applications need to support an alternate trait-based density API as well. It'll be hard to retrofit later if needed, since everybody would have to add it to their densities.

I feel it would just be nice to test if something is a density, if we don't have an abstract type for it. Also for debugging, etc. Do you think isdensity would complicate the API too much?

[devmotion]

I don't mind if it's added, I was just curious if there exists already a specific example where it is needed.

[oschulz]

Ok, maybe I make a quick draft for the package?

[tpapp]

I think that a minimal API package for log densities will make sense when the API is mature. For LogDensityProblems, this is not the case at this point, I plan various changes. So I am not in favor of a shared symbol at this point.

[oschulz]

Ok, no worries I'll choose a different name then. Sorry for spamming your repo with this longish thread, in this case.

@devmotion, Ok, if I CC you on my draft repo? I'd really like to do something that would be acceptable to the community, so Turing, etc., and I'd value your input.

[devmotion]

Sure, I can have a look.

[devmotion]

@phipsgabler This is the discussion about DensityInterface from a year ago.

Link to your PoC: https://github.com/phipsgabler/LogDensityProblems.jl/tree/phg/densityinterface

I think a bit less intrusive and breaking would be to start with just supporting the DensityInterface API on top but keeping logdensity and logdensity_and_gradient as they are without any deprecations. I.e., just adding logdensityof and DensityKind.

(Another comment regarding the sketch: I think one should refrain from making additional changes unrelated to supporting DensityInterface such as adding StaticNumbers. I'm not convinced yet that this is a good idea since most functionality is not needed, and my experience with Static.jl leads me to believe that these kinds of packages are prone to causing a lot of invalidations.)