Improved numerical stability of binomial_coefficient_log #1614

martinmodrak · 2020-01-14T12:27:40Z

Summary

Reimplemented binomial_coefficient_log. Previously, this was done via calls to lgamma or normal approximation to the binomial coefficient. The current version delegates to lbeta for most cases, using the fact that binomial_coefficient_log(N,n) == -lbeta(N - n + 1, n + 1) - log(N + 1) . This in turn relies on the stability improvements made to lbeta in #1612 (this PR depends on #1612 - the branch has been merged here as otherwise the tests don't past).

Additionally, better handling of edge cases has been improved. The function is now valid to call whenever calling lgamma(N + 1) - lgamma(n + 1) - lgamma(N + 1 - n) would be, i.e. for n >= -1, N >= -1, N + 1 >= n. With the cases when equality holds return 0 or infinity with correct sign.

Tests

Test against values and derivatives precomputed in Mathematica.
- For small arguments, I needed to increase the test tolerance for derivatives to 1e-7, but I think it is not a big deal.
Test against the successor identity ( binomial_coefficient_log(n, k) == binomial_coefficient_log(n - 1, k - 1) + log(n) - log(k)). The code could be improved once [WIP] Tools for testing identities #1677 is done.
Test for lbeta against a few identities. It is a bit weaker than it should be, because I triggerred a bug Derivatives of log_sum_exp are wrong for large arguments #1679 which causes false positives in the test. The code can be simplified once [WIP] Tools for testing identities #1677 is done.
Test for edge and error cases

Side Effects

I don't think so, although the code might slightly change performance in either direction.

Checklist

Math issue binomial_coefficient_log produces wrong results when N >> n #1592
Copyright holder: Institute of Microbiology of the Czech Academy of Sciences

By submitting this pull request, the copyright holder is agreeing to the license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)
- Documentation: CC-BY 4.0 (https://creativecommons.org/licenses/by/4.0/)
the basic tests are passing
- unit tests pass (to run, use: ./runTests.py test/unit)
- header checks pass, (make test-headers)
- docs build, (make doxygen)
- code passes the built in C++ standards checks (make cpplint)
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested

…omial_coefficient_log

…gs/RELEASE_500/final)

…omial_coefficient_log

…s' into bugfix/1592-binomial_coefficient_log

…gs/RELEASE_500/final)

…s' into bugfix/1592-binomial_coefficient_log

…gs/RELEASE_500/final)

…mial_coefficient_log

…omial_coefficient_log

…mial_coefficient_log

…gs/RELEASE_500/final)

…stable/2017-11-14)

…_log' into bugfix/1592-binomial_coefficient_log

…stable/2017-11-14)

martinmodrak · 2020-03-09T20:53:41Z

@bbbales2 Can I ask for a review? This is a prerequisite for the neg. binomial PR (#1497).

bbbales2 · 2020-03-09T21:06:44Z

Sure yeah. Did we do a pull related to this for the last release? I went to the closed pull reqs. and didn't find it.

bbbales2 · 2020-03-09T21:07:48Z

Found it: #1612

stan-buildbot · 2020-03-09T21:12:54Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.85	4.85	1.0	0.06% faster
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.02	2.03% faster
eight_schools/eight_schools.stan	0.09	0.09	1.03	2.45% faster
gp_regr/gp_regr.stan	0.22	0.22	1.01	0.52% faster
irt_2pl/irt_2pl.stan	6.47	6.49	1.0	-0.39% slower
performance.compilation	88.71	86.23	1.03	2.8% faster
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	7.53	7.52	1.0	0.07% faster
pkpd/one_comp_mm_elim_abs.stan	20.6	21.24	0.97	-3.12% slower
sir/sir.stan	95.43	95.89	1.0	-0.48% slower
gp_regr/gen_gp_data.stan	0.05	0.05	0.99	-0.61% slower
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.96	2.96	1.0	0.04% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.32	0.31	1.03	2.61% faster
arK/arK.stan	1.75	1.74	1.0	0.22% faster
arma/arma.stan	0.66	0.66	1.0	-0.1% slower
garch/garch.stan	0.51	0.52	1.0	-0.26% slower
Mean result: 1.00413980921

Jenkins Console Log
Blue Ocean
Commit hash: 47c7962

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

bbbales2

Looks good, just a couple explanations so I can double check yah on the last couple bits.

test/unit/math/rev/fun/binomial_coefficient_log_test.cpp

test/unit/math/rev/fun/lbeta_test.cpp

bbbales2 · 2020-03-10T19:47:37Z

stan/math/prim/fun/binomial_coefficient_log.hpp

+      }
+    }
+    if (!is_constant_all<T_k>::value) {
+      if (k_ == 0 && n_ == -1.0) {


Could you explain the gradient logic for k here?

I think I kinda followed the ones for n but these ones confused me.

Not sure what to add. The only fact I use is that lim x-> 0 digamma(x) from above is negative infinity. Mentioned that in the comments, if that is still confusing, let me know (it is also possible I am missing some of the edge cases).

Ah yeah, that makes sense.

bbbales2 · 2020-03-10T19:49:29Z

Just for clarity just explain here on Github then we can throw them in the code if that seems worth the 8 hours of testing or just staple it on the next pull.

mcol

I left some style comments. As @bbbales2 is actually reviewing this, consider my comments optional unless Ben says otherwise.

stan/math/prim/fun/binomial_coefficient_log.hpp

mcol · 2020-03-10T20:17:05Z

stan/math/prim/fun/binomial_coefficient_log.hpp

+        if (k_ == 0) {
+          ops_partials.edge1_.partials_[0] = 0;
+        } else {
+          ops_partials.edge1_.partials_[0] = stan::math::NEGATIVE_INFTY;


This could be ops_partials.edge1_.partials_[0] = k_dbl == 0 ? 0 : NEGATIVE_INFTY.

Yes, I just thought the ternary operator is discouraged. But if that is the Stan style, I can happily use it.

I don't think it's discouraged, I've been suggested to use it in similar cases. It's a matter of preference, if you find your version clearer, don't change it. :) For me, given that there are a series of nested ifs, it helped reducing the number of things I had to keep in mind while reading.

…stable/2017-11-14)

martinmodrak · 2020-03-12T10:00:24Z

Thanks @mcol and @bbbales2 for the suggestions I implemented most of them (the others I commented and left as unresolved).

stan-buildbot · 2020-03-12T17:26:32Z

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
gp_pois_regr/gp_pois_regr.stan	4.84	4.93	0.98	-1.95% slower
low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	0.98	-1.96% slower
eight_schools/eight_schools.stan	0.09	0.09	0.94	-5.94% slower
gp_regr/gp_regr.stan	0.22	0.22	0.98	-1.66% slower
irt_2pl/irt_2pl.stan	6.44	6.45	1.0	-0.11% slower
performance.compilation	87.51	88.01	0.99	-0.57% slower
low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	7.53	7.53	1.0	0.07% faster
pkpd/one_comp_mm_elim_abs.stan	21.55	20.83	1.03	3.35% faster
sir/sir.stan	95.48	94.15	1.01	1.39% faster
gp_regr/gen_gp_data.stan	0.05	0.05	1.01	0.59% faster
low_dim_gauss_mix/low_dim_gauss_mix.stan	2.96	2.95	1.0	0.11% faster
pkpd/sim_one_comp_mm_elim_abs.stan	0.31	0.31	1.0	-0.19% slower
arK/arK.stan	1.75	1.74	1.0	0.28% faster
arma/arma.stan	0.66	0.66	1.0	0.32% faster
garch/garch.stan	0.51	0.52	0.99	-0.92% slower
Mean result: 0.995619320051

Jenkins Console Log
Blue Ocean
Commit hash: 6aed316

Machine information

ProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010

CPU:
Intel(R) Xeon(R) CPU E5-1680 v2 @ 3.00GHz

G++:
Configured with: --prefix=/Applications/Xcode.app/Contents/Developer/usr --with-gxx-include-dir=/usr/include/c++/4.2.1
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

Clang:
Apple LLVM version 7.0.2 (clang-700.1.81)
Target: x86_64-apple-darwin15.6.0
Thread model: posix

bbbales2

Awesome!

bbbales2 · 2020-03-13T00:08:46Z

stan/math/prim/fun/binomial_coefficient_log.hpp

+      }
+    }
+    if (!is_constant_all<T_k>::value) {
+      if (k_ == 0 && n_ == -1.0) {


Ah yeah, that makes sense.

martinmodrak and others added 5 commits January 14, 2020 11:36

Failing test

8a8ea67

Fixes stan-dev#1592

7991f40

Test for derivatives, docs

77d1e85

Merge branch 'bugfix/1611-lbeta-large-arguments' into bugfix/1592-bin…

8bef133

…omial_coefficient_log

[Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (ta…

118adb2

…gs/RELEASE_500/final)

martinmodrak mentioned this pull request Jan 14, 2020

Fixing negative binomial phi cutoff #1497

Merged

5 tasks

martinmodrak and others added 19 commits January 14, 2020 13:37

Lint, test sources

3578ff6

Mismatched define

e03205d

Merge branch 'bugfix/1611-lbeta-large-arguments' into bugfix/1592-bin…

e44bd90

…omial_coefficient_log

Improved code style

035212d

Merge remote-tracking branch 'origin/bugfix/1611-lbeta-large-argument…

e7194bf

…s' into bugfix/1592-binomial_coefficient_log

[Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (ta…

7922dd5

…gs/RELEASE_500/final)

Merge remote-tracking branch 'origin/bugfix/1611-lbeta-large-argument…

3fd0f3f

…s' into bugfix/1592-binomial_coefficient_log

Fixed lint error

6d6ac35

[Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (ta…

3a320df

…gs/RELEASE_500/final)

Merge remote-tracking branch 'stan-dev/develop' into bugfix/1592-bino…

37760cd

…mial_coefficient_log

Removed problematic constexpr

6ebef3e

Merge branch 'bugfix/1611-lbeta-large-arguments' into bugfix/1592-bin…

3df81a2

…omial_coefficient_log

Merge branch 'bugfix/1611-lbeta-large-arguments' into bugfix/1592-bin…

d0f7f45

…omial_coefficient_log

Merge remote-tracking branch 'stan-dev/develop' into bugfix/1592-bino…

c01d2be

…mial_coefficient_log

Updated and cleaned up test generating code

74d2824

[Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (ta…

dd59b06

…gs/RELEASE_500/final)

Improved precomputed tests, first pass at formula test

5bda4f5

Work towards formula tests (failing)

cd3381a

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

d014464

…stable/2017-11-14)

martinmodrak mentioned this pull request Feb 4, 2020

neg_binomial_2_log_lpmf is not stable for large phi #1495

Closed

martinmodrak and others added 4 commits March 7, 2020 18:10

Merge branch 'develop' into bugfix/1592-binomial_coefficient_log

9a5ee11

Passing identity tests, ad test

51ef491

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

215cad0

…stable/2017-11-14)

Resolved int and fwd issues

71cbd40

martinmodrak and others added 5 commits March 7, 2020 23:07

Merge remote-tracking branch 'origin/bugfix/1592-binomial_coefficient…

cfa5b81

…_log' into bugfix/1592-binomial_coefficient_log

Attempt to handle edge cases for derivatives

4fa9c52

Moved to operands_and_partials, fixed edge cases

eab19c3

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

75a6f73

…stable/2017-11-14)

Fixed lint

47c7962

bbbales2 self-assigned this Mar 9, 2020

bbbales2 requested changes Mar 10, 2020

View reviewed changes

mcol reviewed Mar 10, 2020

View reviewed changes

martinmodrak and others added 3 commits March 12, 2020 10:56

Handling review comments

cf05d40

Merge commit 'f2a3c1a6f8de1d5bd380c265f0bc92472dce2d42' into HEAD

5b93278

[Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/…

c84a259

…stable/2017-11-14)

One more edge case-test+fix

6aed316

bbbales2 approved these changes Mar 13, 2020

View reviewed changes

bbbales2 merged commit fe3a41c into stan-dev:develop Mar 13, 2020

mcol linked an issue Mar 13, 2020 that may be closed by this pull request

binomial_coefficient_log produces wrong results when N >> n #1592

Closed

SteveBronder mentioned this pull request Apr 16, 2020

Stan Math 3.2 release #1826

Closed

bbbales2 mentioned this pull request Apr 20, 2020

Stanc3 release for Cmdstan 2.23 stan-dev/stanc3#498

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Improved numerical stability of binomial_coefficient_log #1614

Improved numerical stability of binomial_coefficient_log #1614

martinmodrak commented Jan 14, 2020 •

edited

Loading

martinmodrak commented Mar 9, 2020

bbbales2 commented Mar 9, 2020

bbbales2 commented Mar 9, 2020

stan-buildbot commented Mar 9, 2020

bbbales2 left a comment

bbbales2 Mar 10, 2020

martinmodrak Mar 12, 2020

bbbales2 Mar 13, 2020

bbbales2 commented Mar 10, 2020

mcol left a comment

mcol Mar 10, 2020

martinmodrak Mar 12, 2020

mcol Mar 12, 2020

martinmodrak commented Mar 12, 2020

stan-buildbot commented Mar 12, 2020

bbbales2 left a comment

bbbales2 Mar 13, 2020

Improved numerical stability of binomial_coefficient_log #1614

Improved numerical stability of binomial_coefficient_log #1614

Conversation

martinmodrak commented Jan 14, 2020 • edited Loading

Summary

Tests

Side Effects

Checklist

martinmodrak commented Mar 9, 2020

bbbales2 commented Mar 9, 2020

bbbales2 commented Mar 9, 2020

stan-buildbot commented Mar 9, 2020

bbbales2 left a comment

Choose a reason for hiding this comment

bbbales2 Mar 10, 2020

Choose a reason for hiding this comment

martinmodrak Mar 12, 2020

Choose a reason for hiding this comment

bbbales2 Mar 13, 2020

Choose a reason for hiding this comment

bbbales2 commented Mar 10, 2020

mcol left a comment

Choose a reason for hiding this comment

mcol Mar 10, 2020

Choose a reason for hiding this comment

martinmodrak Mar 12, 2020

Choose a reason for hiding this comment

mcol Mar 12, 2020

Choose a reason for hiding this comment

martinmodrak commented Mar 12, 2020

stan-buildbot commented Mar 12, 2020

bbbales2 left a comment

Choose a reason for hiding this comment

bbbales2 Mar 13, 2020

Choose a reason for hiding this comment

martinmodrak commented Jan 14, 2020 •

edited

Loading