new functions: sort, rank, Kahan summation, ???

[Bob Carpenter]

[re-posted on stan-dev]

sort() is definitely useful.

If we have sort(), is there a reason to also have rank()?
(That wasn't meant to be a rhetorical question.)

I'd like to see at least one realistic example where something
like Kahan summation makes a difference before including it.
Did you work on it because of an applied problem you were
having?

Does anyone else on the list have experience with this?

Did you test gradients for summation? Writing custom
gradients is critical for speed, but the gradients need to be
tested even if they're not custom. You should just be able
to use the gradient already there for sum.

Also, which signatures did you go with for sum()? It
should probably work over the same types as sort: real[],
int[], vector and row_vector.

Could you perhaps break these into two separate pull requests
(sort/rank vs. summation).

Thanks!

- Bob


On 5/2/13 2:09 AM, Marco wrote:
> I pushed sort and rank functions to my fork now.
>
> If I have enough time, tomorrow I'll test and add them to the reference manual (as well as the float point sum precision
> functions) and then create a pull request.
>
>> There's a section on sorting and interpolation.
>> But you're right --- sort() itself isn't on the list!
>>
>> Issue 1 is which types it'll work for. Probably
>> just
>>
>> real[]
>> int[]
>> vector
>> row_vector
>>
>> We would not want it to sort in place. So you'd have
>> to copy to a std::vector, apply std::sort, then return
>> (for real[] and int[]) or use to map an Eigen vector (for
>> vector and row_vector).
>>
>> The to-do list did have this:
>>
>> -- rank(v,s): number of components of v less than v[s]
>> -- ranked(v,s): s-th smallest component of v
>>
>> These are BUGS functions.
>>
>> http://www.mrc-bsu.cam.ac.uk/bugs/winbugs/manual14.pdf
>>
>> The first can be implemented
>> in constant memory and O(s) time. The second can be
>> implemented in O(s) memory and O(n log min(s,n-s)) time.
>>
>> With sort(), I don't think we'd need these functions.
>> See this discussion on the JAGS list (JAGS uses sort()
>> instead of rank() and ranked()):
>>
>> http://sourceforge.net/p/mcmc-jags/discussion/610037/thread/64ce3a9c
>>
>> - Bob
>>
>

[Marco]

> sort() is definitely useful.
>
> If we have sort(), is there a reason to also have rank()?
> (That wasn't meant to be a rhetorical question.)

I thought about efficiency, as using sorting could be slower for big
data, but I don't know how often the need of these functions show up in
the models, probably not often as JAGS hasn't defined them (in fact JAGS
has rank(v) similar to R's one which serves a related, but different
propose).

>
> I'd like to see at least one realistic example where something
> like Kahan summation makes a difference before including it.

I think you are right and the said truth is that I couldn't come up with
a definitive case where any of them are better than sum(), so there's
no reason to include them, at least not yet. It seems that these
algorithms only help with very big data, but very big data (as my model)
takes too long to converge.

> Did you work on it because of an applied problem you were
> having?

The motivation is that model I described a few weeks ago in stan-users.
I'm still having a lot a trouble to deal with it.

>
> Does anyone else on the list have experience with this?
>
> Did you test gradients for summation? Writing custom
> gradients is critical for speed, but the gradients need to be
> tested even if they're not custom. You should just be able
> to use the gradient already there for sum.

I've done, but I sadly lack knowledge to extract useful information from
them.

>
> Also, which signatures did you go with for sum()? It
> should probably work over the same types as sort: real[],
> int[], vector and row_vector.

I tried a lot of different ways, but I decided to use just a single
eigen vector to sum.

vector [nZeros] summads;
for (i in 1:nZeros)
summands[i] <- log1p_exp(-zeta0[i]-exp(eta0[i]));
1 ~ bernoulli_logit(zeta0);

yy1 ~ poisson_log(eta1);
0 ~ bernoulli_logit(zeta1);

>
> Could you perhaps break these into two separate pull requests
> (sort/rank vs. summation).

I started them in two different branches and it looks like it was a wise
decision after all.

[Bob Carpenter]

On 5/2/13 4:00 PM, Marco wrote:
>
>> sort() is definitely useful.
>>
>> If we have sort(), is there a reason to also have rank()?
>> (That wasn't meant to be a rhetorical question.)
> I thought about efficiency, as using sorting could be slower for big data, but I don't know how often the need of these
> functions show up in the models, probably not often as JAGS hasn't defined them (in fact JAGS has rank(v) similar to R's
> one which serves a related, but different propose).
>
>
>>
>> I'd like to see at least one realistic example where something
>> like Kahan summation makes a difference before including it.
> I think you are right and the said truth is that I couldn't come up with a definitive case where any of them are better
> than sum(), so there's no reason to include them, at least not yet. It seems that these algorithms only help with very
> big data, but very big data (as my model) takes too long to converge.
>
>
>> Did you work on it because of an applied problem you were
>> having?
> The motivation is that model I described a few weeks ago in stan-users. I'm still having a lot a trouble to deal with it.

Are you sure it's a problem with summation precision per se?

I want to understand where you need to add numbers of
vastly different scales AND the lower-order numbers
matter. That is, when do you do (1e-30 + 1e10) and
the 1e-30 term matters in the sum? We try to handle
the gradients separately, and in the usual ways to do
sums (+ and sum()), the gradients are handled independently,
so you get the right answer there even if the 1e-30 term
drops out of the sum due to underflow.

So what I'm wondering is where higher-precision summation
will make a difference.

The integrator we use to simulate the Hamiltonian dynamics
has a higher error than floating point precision.

>>
>> Does anyone else on the list have experience with this?
>>
>> Did you test gradients for summation? Writing custom
>> gradients is critical for speed, but the gradients need to be
>> tested even if they're not custom. You should just be able
>> to use the gradient already there for sum.
> I've done, but I sadly lack knowledge to extract useful information from them.

You basically want to implement a complex summation
so that the gradients are done in the same way as for
sum(), but the value is calculated to higher precision.
In a sum(), the partial is just 1 for all of the components,
so there's really no more precision to be gained there.

>> Also, which signatures did you go with for sum()? It
>> should probably work over the same types as sort: real[],
>> int[], vector and row_vector.
> I tried a lot of different ways, but I decided to use just a single eigen vector to sum.
>
> vector [nZeros] summads;
> for (i in 1:nZeros)
> summands[i] <- log1p_exp(-zeta0[i]-exp(eta0[i]));
> 1 ~ bernoulli_logit(zeta0);
>
> yy1 ~ poisson_log(eta1);
> 0 ~ bernoulli_logit(zeta1);

What I meant is that if you implement kahan_sum(), you probably
want to implement it for real[], vector, and row-vector. Forget
what I said about int[], because you won't get any extra precision.

>> Could you perhaps break these into two separate pull requests
>> (sort/rank vs. summation).
> I started them in two different branches and it looks like it was a wise decision after all.

Great. Send us the pull request for sort() and rank() when
you have them. Just please make sure to run fairly extensive
tests. Testing these functions is usually harder than writing
them, so it's not doing us much of a favor to write the function without
extensive testing. For sort, you'd want to test zero size arrays/vectors
as input, and cases where all the numbers are the same, and some random
ones. There's no issue with gradients here because neither rank() nor
sort() have gradients to calculate.

- Bob

[Marco]

> Are you sure it's a problem with summation precision per se?
>
> I want to understand where you need to add numbers of
> vastly different scales AND the lower-order numbers
> matter. That is, when do you do (1e-30 + 1e10) and
> the 1e-30 term matters in the sum? We try to handle
> the gradients separately, and in the usual ways to do
> sums (+ and sum()), the gradients are handled independently,
> so you get the right answer there even if the 1e-30 term
> drops out of the sum due to underflow.
>
> So what I'm wondering is where higher-precision summation
> will make a difference.
>
> The integrator we use to simulate the Hamiltonian dynamics
> has a higher error than floating point precision.

Thanks very much for helping me understand the point. I was ignoring the
importance of gradients when I started creating those functions. When I
started reading the files in src/prob/distribution yesterday, however, I
began to get a better idea of how Stan works to calculate log prob and
gradients. And as my model derivatives are all in closed form, I'll try
to write a distribution function for it directly and check if I can get
interesting results :)

> Great. Send us the pull request for sort() and rank() when
> you have them. Just please make sure to run fairly extensive
> tests. Testing these functions is usually harder than writing
> them, so it's not doing us much of a favor to write the function without
> extensive testing. For sort, you'd want to test zero size arrays/vectors
> as input, and cases where all the numbers are the same, and some random
> ones. There's no issue with gradients here because neither rank() nor
> sort() have gradients to calculate.

Sure, I tested them today and they are working fine. I just sent the
pull request.