Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve algorithm for sampling Beta #1000

Merged
merged 11 commits into from Sep 8, 2020
Merged

Conversation

vks
Copy link
Collaborator

@vks vks commented Jul 17, 2020

Fixes #999.

cc @rasa200

@saona-raimundo
Copy link
Collaborator

saona-raimundo commented Jul 17, 2020

Great!

Do you mind sharing the reference for the algorithm you are implementing here?

Edit:

It was part of your commits and the new documentation.

Reference:

R. C. H. Cheng (1978).
Generating beta variates with nonintegral shape parameters.
Communications of the ACM 21, 317-322.
https://doi.org/10.1145/359460.359482

Thank you!!

@dhardy
Copy link
Member

dhardy commented Jul 17, 2020

That was quick @vks! Okay if I review some time this weekend? (Or is someone else wants to, I won't hog all the work 😉)

@vks
Copy link
Collaborator Author

vks commented Jul 18, 2020

@rasa200 I see you already found the source, sorry for not posting it in the PR description. The PDF is a bit tricky to access, but you could also look at the R source code for rbeta (which I personally avoided, because it is GPL code, and Rand's license is more liberal).

@dhardy Sure!

Copy link
Member

@dhardy dhardy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just starting to review, though I probably won't finish today.

@vks could you compare benchmarks before/after with a few parameter selections? Speed is not really a deciding factor in this case but it would be nice to have some idea how this performs.

rand_distr/src/gamma.rs Outdated Show resolved Hide resolved
@vks
Copy link
Collaborator Author

vks commented Jul 19, 2020

It seems like the performance of this PR is actually worse, at least for the parameters I looked into:

On master:

test distr_beta_large_param_different      ... bench:     146,886 ns/iter (+/- 2,585) = 54 MB/s
test distr_beta_large_param_similar        ... bench:     145,855 ns/iter (+/- 980) = 54 MB/s
test distr_beta_mixed_param                ... bench:     248,922 ns/iter (+/- 3,614) = 32 MB/s
test distr_beta_small_param                ... bench:     357,560 ns/iter (+/- 4,676) = 22 MB/s

This PR:

test distr_beta_large_param_different      ... bench:     252,111 ns/iter (+/- 5,843) = 31 MB/s
test distr_beta_large_param_similar        ... bench:     268,847 ns/iter (+/- 23,299) = 29 MB/s
test distr_beta_mixed_param                ... bench:     282,934 ns/iter (+/- 46,543) = 28 MB/s
test distr_beta_small_param                ... bench:     350,247 ns/iter (+/- 40,186) = 22 MB/s

@vks
Copy link
Collaborator Author

vks commented Jul 30, 2020

I think we want this new algorithm for small parameters, but it might make sense to use the old algorithm for large parameters. I also did not really optimize the algorithm in this PR.

@dhardy
Copy link
Member

dhardy commented Jul 30, 2020

Speed isn't a critical factor (though can be nice). Sorry that I didn't do the review yet.

@vks
Copy link
Collaborator Author

vks commented Jul 30, 2020

Speed isn't a critical factor (though can be nice).

Fair enough! However, the algorithm I implemented is actually two rejection-sampling algorithms: BC is used for min(alpha, beta) <= 1, BB is used for larger parameters. It would be unfortunate to switch to a less efficient algorithm for large parameters. However, the reference claims that BB is more efficient than sampling the gamma distribution twice (our current algorithm), so maybe something is off with my implementation/comparison (or maybe their results don't apply to modern hardware).

Of course, correctness is much more important, but I did not compare the algorithms for large parameters, so I don't really have an argument for using BB over the current implementation. 😅

Sorry that I didn't do the review yet.

No worries, this is not crucial or time-critical.

@vks
Copy link
Collaborator Author

vks commented Aug 4, 2020

I rebased on master, fixed the benchmarks (that were broken on master) and added them to CI (we did not test the rand_distr benchmarks, so it was not noticed that they were broken).

The last two commits could be merged separately and are arguably more urgent than the rest of this PR.

@vks vks mentioned this pull request Aug 5, 2020
3 tasks
@dhardy
Copy link
Member

dhardy commented Aug 28, 2020

Sorry that I still didn't get to this. @rasa200 would you care to do a review? The main point is to have a second person check correctness (or at least conformance to the paper).

This will have to target rand_distr v0.4, which will hopefully come out soon after rand v0.8.

@saona-raimundo
Copy link
Collaborator

Yes! I will try my best!

Copy link
Collaborator

@saona-raimundo saona-raimundo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The algorithmic part was super clear and followed the paper line-by-line. In doing so, a typo was reproduced (there is a number that should be 72 instead of 18, see comments).

The pre-simulation constants should be corrected. By making the choice of setting the parameters a and b from the beginning, before the choice of the variant BB or BC, and the introduction of the variable switched_params, there is a disparity with the formulas presented in the paper. In the paper, they set the variables a and b as the min and max respectively for the BB algorithm, but they do the other way around for the BC algorithm!! Making confusion in the formulas when coding! (You can check this in the first line of each algorihtm in the paper).

rand_distr/src/gamma.rs Outdated Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Show resolved Hide resolved
rand_distr/src/gamma.rs Outdated Show resolved Hide resolved
rand_distr/src/gamma.rs Outdated Show resolved Hide resolved
@vks
Copy link
Collaborator Author

vks commented Sep 6, 2020

Thanks for the detailed review! I overlooked that the definition for a and b was different in the two algorithms. To fix this, I decided to switch a and b in algorithm BC, so that it corresponds more closely to the reference. I also fixed the wrong switching condition.

I don't think there is a typo in the paper, as far as I can see, the division by for was absorbed in the numerator, so the denominator should not be affected.

This should be faster than the gamma variate transformation we are currently
using, and it seems to work better for parameters smaller than one.

The algorithm is also used by the R language, however I did not consult their
implementation in order to avoid licensing problems.

Reference:

R. C. H. Cheng (1978).
Generating beta variates with nonintegral shape parameters.
Communications of the ACM 21, 317-322.
https://doi.org/10.1145/359460.359482
@vks
Copy link
Collaborator Author

vks commented Sep 6, 2020

Now that the algorithm is fixed, we are doing much better with performance:

This PR:

test distr_beta_large_param_different      ... bench:      66,837 ns/iter (+/- 3,437) = 119 MB/s
test distr_beta_large_param_similar        ... bench:      73,134 ns/iter (+/- 7,125) = 109 MB/s
test distr_beta_mixed_param                ... bench:      95,295 ns/iter (+/- 7,092) = 83 MB/s
test distr_beta_small_param                ... bench:      93,504 ns/iter (+/- 5,298) = 85 MB/s

master:

test distr_beta_large_param_different      ... bench:     146,886 ns/iter (+/- 2,585) = 54 MB/s
test distr_beta_large_param_similar        ... bench:     145,855 ns/iter (+/- 980) = 54 MB/s
test distr_beta_mixed_param                ... bench:     248,922 ns/iter (+/- 3,614) = 32 MB/s
test distr_beta_small_param                ... bench:     357,560 ns/iter (+/- 4,676) = 22 MB/s

@saona-raimundo
Copy link
Collaborator

Great!! Thank you very much!!

I agree with all the changes, and sorry about the "typo" of the constants, it was my own mistake: it is right the way it is in the code and the paper.

@@ -4,6 +4,9 @@ All notable changes to this project will be documented in this file.
The format is based on [Keep a Changelog](http://keepachangelog.com/en/1.0.0/)
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).

## [Unreleased]
- Improve algorithm for sampling `Beta` (#1000)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could this be a little more specific?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I changed it to "New Beta sampling algorithm for improved performance and accuracy".

Copy link
Member

@dhardy dhardy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you both happy with the code now?

@saona-raimundo
Copy link
Collaborator

I am happy with the code :)

@vks
Copy link
Collaborator Author

vks commented Sep 7, 2020

@dhardy I'm happy with the code, too. Could you please approve the PR?

Copy link
Member

@dhardy dhardy left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sure thing. Thanks for the review, Raimundo.

@vks vks merged commit 7a1f51c into rust-random:master Sep 8, 2020
@vks vks deleted the improve-beta branch September 8, 2020 07:53
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Unexpected sample values from beta distribution for small parameters
3 participants