Can we have unnormalised densities in user functions?

General
,
1

Hi!

With reduce_sum (and map_rect) we can now parallelise Stan programs nicely. However, this comes at the price of being forced to use the standardised densities.

So while one can write this in the model block:

counts ~ poisson_log(mu);

which avoids calculating log-gamma functions needed for normalisation. Now, when you use the new parallelsiation facilities, then the log-lik must be written in a user defined function, but then users must use

log_lik_term += poisson_log(counts | mu);

… and this is wasting a lot of computational resources since the normalisation constant is being computed which is often not needed.

So can we solve this?

How about we let users define mydensity_ lpdf where they can simply use the usual ~ notation. The function would not declare a returned value which would be auto-generated from the stan parser; like:

real mydensity_lpdf(int[] counts, vector mu) {
   counts ~ poisson_log(mu);
}

This would only be compatible with the current language when we disallow a return value whenever the ~ notation is used.

Alternatively, one could turn calls to poisson_log(counts | mu) in users functions into versions which do not estimate the normalisation constant whenever the user density is used with the ~ notation.

I don’t really mind how we turn these normalisations off, but it’s a bit sad to force everyone to pay some performance penalty when you go parallel merely for the normalisation which is often not used. The alternative for users is write their pdf’s tailored for their case if that is possible (based on using sum or whatever).

Tagging some compiler experts @Bob_Carpenter, @rok_cesnovar, @nhuurre, @rybern … apologies whomever I missed and wanted a heads-up here.

Sebastian

3 Likes
2

It’s not documented but

counts ~ poisson_log(mu);

is synonymous with

target += poisson_propto_log(counts, mu);

It also works inside user-defined _lpdf/_lpmf functions but has no effect in other user-defined functions.

The infix may be changed in the future.

3 Likes
3

Sorry, that needs some qualification. _propto infix does not mean don’t normalize, it means inherit the behavior from the caller. Non-lpdf functions always normalize. Lpdf functions normalize if they’re called without propto. Model block does not normalize. Transformed data, transformed parameters and generated quantities blocks won’t even compile if there’s a propto!

reduce_sum does not pass the propto flag so anything calculated there will always normalize.

3 Likes
4

Hi!

So this is a valid user-defined function?

real mydensity_propto(int[] counts, vector mu) {
   real log_lik = 0.0;
   log_lik += poisson_propto_log(counts, mu);
   return(log_lik);
}

And it would give me a non-normalised Poisson density in my users function?

So if the function given to reduce_sum simply always works non-normalised, then that is just fine, right?

Basically, there is no magic switching between normalised / non-normalised available for user-supplied functions which hard-code the propto thing.

Is that correct?

Thanks!

Sebastian

5

Not quite. The function definition is

real mydensity_lpmf(int[] counts, vector mu) {
   real log_lik = 0.0;
   log_lik += poisson_propto_log(counts, mu);
   return(log_lik);
}

Then, if it is called like mydensity_lpmf(counts|mu) it uses normalized poisson_log but if it is called as mydensity_propto_lpmf(counts|mu) is uses unnormalized poisson_log but only if the call is in the model block or in an unnormalized lpdf function.

At C++ level lpdf functions have template argument propto__ that is passed to the nested lpdf functions.

reduce_sum(mydensity_propto_lpdf, ...) does not work because propto version is not a real function, it’s the same function with different template argument.

2 Likes
6

So what do we need to do to allow map_rect and reduce_sum so that we can pass along the propto flag?

Or would it work out to define mydensity_lpmf as you suggest and then pass to reduce_sum the function name mydensity_propto_lpmf (which would be called in the model block). Sorry, you wrote this does not work.

Or are we in the situation right now that as it stands I won’t be able to get the un-normalized density in a function passed to reduce_sum as I get you.

Do we need a reduce_sum_density which you can only give densities or can we expand what we have so that things get passed on? Any ideas?

1 Like
7

Right now this Stan code

reduce_sum(mydensity, counts, gs, mu)

creates a struct

struct mydensity_rsfunctor__ {
template <typename T3__>
typename boost::math::tools::promote_args<T3__>::type
operator()(const int& start, const int& end, const std::vector<int>& counts,
           std::ostream* pstream__, const Eigen::Matrix<T3__, -1, 1>& mu)  const 
{
return mydensity(start + 1, end + 1, counts, mu, pstream__);
}
};

and then calls

reduce_sum<mydensity_rsfunctor__>(x, 1, pstream__, mu)

For mydensity_lpmf it doesn’t compile because

struct mydensity_lpmf_rsfunctor__ {
template <bool propto__, typename T3__>
typename boost::math::tools::promote_args<T3__>::type
operator()(const int& start, const int& end, const std::vector<int>& counts,
           std::ostream* pstream__, const Eigen::Matrix<T3__, -1, 1>& mu)  const 
{
return mydensity_lpmf(start + 1, end + 1, counts, mu, pstream__);
}
};

expects a template argument propto__ (though forgets use it…) which reduce_sum doesn’t provide.

I think we could generate functors like

struct mydensity_lpmf_rsfunctor__ {
template <typename T3__>
typename boost::math::tools::promote_args<T3__>::type
operator()(const int& start, const int& end, const std::vector<int>& counts,
           std::ostream* pstream__, const Eigen::Matrix<T3__, -1, 1>& mu)  const 
{
return mydensity_lpmf<false>(start + 1, end + 1, counts, mu, pstream__);
}
};

struct mydensity_lpmf_proptorsfunctor__ {
template <typename T3__>
typename boost::math::tools::promote_args<T3__>::type
operator()(const int& start, const int& end, const std::vector<int>& counts,
           std::ostream* pstream__, const Eigen::Matrix<T3__, -1, 1>& mu)  const 
{
return mydensity_lpmf<true>(start + 1, end + 1, counts, mu, pstream__);
}
};

and then

reduce_sum(mydensity_lpmf, counts, gs, mu)
reduce_sum_propto(mydensity_lpmf, counts, gs, mu)

translates to

reduce_sum<mydensity_lpmf_rsfunctor__>(x, 1, pstream__, mu)
reduce_sum<mydensity_lpmf_proptorsfunctor__>(x, 1, pstream__, mu)

so that the first computes all normalizing constants and the second doesn’t.

2 Likes
8

Ok, I need to fix that for the release. With propto=false.

I am fine with this. What if the users names their functions with _propto from the start? So

real mydensity_propto_lpmf(int[] counts, vector mu) {
   real log_lik = 0.0;
   log_lik += poisson_propto_log(counts, mu);
   return(log_lik);
}

would that be an option to avoid more function signatures?

P.S.: @wds15 no need to list me in the compiler experts :) I am a newbie.

9

IMO propto=true is fine. It has no effect except when using these undocumented _propto_ functions. Also, an _lpmf function in reduce_sum is obviously a weird hack because the signature doesn’t make sense. It has to be

mydensity_lpmf(start| end, slice, ...);

but this implies it’s a probability distribution for the start index. For the function to make sense the slice should be the first argument and the indices must come after it. (We could allow both orderings; there’s no ambiguity and the math lib has many overloaded functions anyway.)

The contributors page disagrees :)

1 Like
10

That is a very good point.

cc @wds15 @bbbales2 : given that we havent released this yet, we can still change this. Thoughts?

11

Oof I’m struggling to catch up.

My opinion is if this is about dropping unnecessary proportionality constants, we should write a design-doc, we should demonstrate the utility (presumably with some sort of hacked up model), and then we should make a plan to implement it in three months for the next release.

Or is this about a use case where we pass an _lpdf/_lpmf function as the first argument to reduce_sum? That seems weird to me cause none of our _lpdf/_lpmfs would like start/end as a 2nd/3rd argument either.

12

I have never thought of the partial reducer function defined by the user as a function to be used on its own. As such the first three arguments were always thought to belong together.

Now, the reasons you bring up here make sense and we should consider changing it, probably…would be a bit of work, but doable (code, tests, doc, intro).

Given we changed that, would we then be able to have non normalize densities in a straightforward way (or is there any other option to have non normalized reduce sum things with what we have).

Curious what @bbbales2 and @stevebronder, @Bob_Carpenter or @mitzimorris has to say about this.

13

My question was about this, yes. Do we disallow the use of users _lpmf/_lpdf functions for this release then? Or supply propto=true/false. It has to be one or the other. C++ compile errors is not an option :)

14

The poisson distribution is such a case where it pays off a lot to not normalize. I can post probably an example model quickly.

Yup…we pass a user defined lpdf and lpmf function to reduce sum. Start/end as 1st and 2nd is probably even more weird… and this logic would not fit into the current system well whereas making start/end as 2nd and 3rd is better, i guess (sorry, i have put this current order forward…i know).

15

Dropping unnecessary proportionality constants is definitely something for 2.24.

What @nhuurre brought up and should be discussed now for this release is what to do if _lpdf/_lpmf is supplied (propto=true/false or disallow).

And whether we want to change the order of arguments. I would want to avoid duplicate signatures for lpmf functions.

16

I agree. This is why I would be ok with changing it now…even though it would be quite a bit of work.

17

Had a think about this.

I think:

reduce_sum(mydensity_lpmf, counts, gs, mu)

Should compile with propto = false (or whichever it is that includes all the constants).

And:

reduce_sum(mydensity, counts, gs, mu)

Should be an error.

Also I’m down to change the partial sum function from:

real f(int, int, real[], ...)

To:

real f(real[], int, int, ...)

It’s just a matter of swapping doc. Good catches on these @nhuurre!

1 Like
18

Great. So let’s change this such that the sliced data is in the first argument of the user supplied reducer, the start and end indices go in 2nd and 3rd place.

As a prospect, we should allow for

reduce_sum_propto

as @nhuurre points out. The propto flavour would lead to calls of the unnormalised densities as described (propto=true as template argument given to the user-functor which must be a lpdf or lpmf). That last thing would be a stanc3 only thing to handle, right? This should be done for the 2.24 release given that it is a new feature.

We probably could have caught this earlier, but I really think here the utility of a release candidate shows as now far more eyes are on this. I like it.

3 Likes
19

The whole reduce_sum_propto will be a stanc3 thing yes. We first need to settle the issue Niko posted above (change to _lupdf) and then decide on the exact formulation.

2 Likes
20

Just edit the .hpp of a few prototypes and change propto=false to propto=true to do experiments with this.

Changing the .hpp file won’t cause a rebuild cause the .stan file will be older than the .hpp.