Incorporate relative + absolute tolerance into grad_2F1 #2175
Conversation
…to grad_2F1 and increasing the max number of default steps. Overall it should terminate sooner for very large gradients . New tests are referenced vs. a quick discrete diff in Mathematica.
…erence interval in the test of the neg_binomial_2_ccdf was too narrow so bumping it higher.
…4.1 (tags/RELEASE_600/final)
|
This works locally for the failing test so if there's a standard way to tickle the test system I'd appreciate any help. |
|
I am able to replicate locally with clang 10.0 on Ubuntu. Not sure on the solution though. |
|
@rok-cesnovar thanks for checking the bug happens for you, I'm a little burned out on reading C++ errors at the moment but I'm sure I can figure it out. I'm on arch linux so it's usually one of the newest compiler versions and going to Ubuntu is like getting in a time machine. |
There was a problem hiding this comment.
I had a closer look. I'm not sure relative error makes sense here. As far as tests, the thing that seemed to matter is the increased number of steps of the algorithm, not the error test change.
If I print out the error on convergence, I'm always seeing really small numbers which makes me think log_t_new is an increment (and so just making sure it is below some value is already a relative error check).
It looks like from the line log_t_new += log(fabs(p)) + log_z; that log_t_new is being incremented to some value, but I think this is actually the log of the terms in the power series of the function and so the terms are getting smaller in the sequence and the log is getting more negative (and so the log_t_new accumulation is only negative, and isn't really converging itself to anything, cuz ideally it goes to -inf which is zero error).
The variable p at the top of the loop looks like incrementally computing new factors to be included in each term in the series (look at the definition at the bottom of the Notation section here -- I wonder if p was for Pochhammer: https://en.wikipedia.org/wiki/Generalized_hypergeometric_function#Notation)
|
@bbbales2 just for context: I re-wrote and re-organized the hypergeometric functions and their gradients last time I worked significantly on Stan, so I'm familiar with the calculations and their challenges. It would take another significant chunk of work to improve them over the current status and that is not what I'm trying to address in the PR. |
|
Yeah that's fine. I don't think we need any major work here. I'm happy just multiplying max_steps by 10 and calling it a day. Since the relative tolerance was new, I tried to verify that tit was behaving well at the switch and I wasn't able to do that. The tests I ran ( |
|
oh this is definitely not worth your time to try and figure out the hypergeometric function notation! It is different everywhere and I standardized our notation across the different functions last time I touched them but the reference are still weird! Your comments are the nudge I needed to look at the parts that were bothering me yesterday so thanks! You got it right that log_t_new is based on the increment (it's a lagged part of each new term in the sum) so I need to actually include |
|
Thanks again for the review, it really helps especially with these obscure calculations... even if I'm the one who re-wrote them... :/ |
|
Looks like there were some failures from the test system that aren't code-related so I upped precision in one test and kicked it off again... |
|
BTW: upshot is that bumping the max_steps to 1e6 is still required to get the boundary value down right so the PR is adding relative tolerance to avoid excessive computation and also bumping max_steps to avoid a failure in Stan models. Downstream tests failed (out of memory) so I re-started those. |
|
Upstream tests fails again, anybody know what to do with this? It's been too long for me. |
|
That is an issue with upstream tests that we are looking into. Will make sure to restart these tests once the problem is fixed. Sorr for the inconvenience. |
Jenkins Console Log Machine informationProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010CPU: G++: Clang: |
|
Wow those ARMA/GARCH performance results |
|
Got it on my list for tomorrow! |
|
I think these convergence checks should be on I tried out: T g_a1_prev = g_a1;
T g_b1_prev = g_b1;
g_a1 += log_g_old_sign[0] > 0 ? exp(log_g_old[0]) : -exp(log_g_old[0]);
g_b1 += log_g_old_sign[1] > 0 ? exp(log_g_old[1]) : -exp(log_g_old[1]);
if (std::abs(value_of_rec(g_a1) - value_of_rec(g_a1_prev)) <=
std::max(std::abs(value_of_rec(g_a1)) * precision, precision) &&
std::abs(value_of_rec(g_b1) - value_of_rec(g_b1_prev)) <=
std::max(std::abs(value_of_rec(g_b1)) * precision, precision)) {
std::cout << "k: " << k << std::endl;
std::cout << "log(fabs(term)): " << log(fabs(term)) << std::endl;
std::cout << "log_t_new: " << value_of(log_t_new) << std::endl;
std::cout << "log_g_old[0]: " << value_of(log_g_old[0]) << std::endl;
std::cout << "log_g_old[1]: " << value_of(log_g_old[1]) << std::endl;
std::cout << "g_old[0]: " << exp(value_of(log_g_old[0])) << std::endl;
std::cout << "g_old[1]: " << exp(value_of(log_g_old[1])) << std::endl;
std::cout << "abs(g_a1) * precision: " << std::setprecision(17) << std::abs(value_of(g_a1)) * precision << std::endl;
std::cout << "abs(g_b1) * precision: " << std::setprecision(17) << std::abs(value_of(g_b1)) * precision << std::endl;
std::cout << "precision: " << exp(log_precision) << std::endl;
std::cout << "g_a1: " << std::setprecision(17) << value_of(g_a1) << std::endl;
std::cout << "g_b1: " << std::setprecision(17) << value_of(g_b1) << std::endl;
return;
}Running The problem is that the thing converging here is If I change to the checks above I get: In this case the gradients are off in the 11th digit, but we're stopping the calculation way sooner (in this case when we're incrementing numbers smaller than |
|
I think I agree, I doubt I can get to this today but I'll do it tomorrow (if the world isn't on fire) |
|
Just to bump this along I went ahead and made the changes. Have a look and revert them or change them if they're not good. |
…4.1 (tags/RELEASE_600/final)
Jenkins Console Log Machine informationProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010CPU: G++: Clang: |
stan/math/prim/fun/grad_2F1.hpp
Outdated
| @@ -86,13 +86,16 @@ void grad_2F1(T& g_a1, T& g_b1, const T& a1, const T& a2, const T& b1, | |||
| log_g_old[1] = log_t_new + log(fabs(term)); | |||
| log_g_old_sign[1] = term >= 0.0 ? log_t_new_sign : -log_t_new_sign; | |||
|
|
|||
| T g_a1_prev = g_a1; | |||
There was a problem hiding this comment.
@bbbales2 thanks for prodding me along, shouldn't we be storing these as just doubles? Is there any performance penalty involved in storing the type?
There was a problem hiding this comment.
Yeah I guess since we're just looking at the values we could do value_of_rec. Will push that up in a sec.
|
So... that increment you calculate as the convergence condition, that is |
|
@sakrejda that makes sense. I switched it to logs. The convergence numbers to compare to above are: |
I think this is a third thing, I should have time later today to update this to what I think it should be and I'll ask you to check it then, ok? I think we're getting there, it's just a weird calculation... |
Sure! Btw I don't think the red X is anything to worry about. Sounds like Jenkins got sad over the weekend. It'll probably fix itself with the next push you make. |
|
I would also merge recent develop in. There were a couple of fixes in the last few days for the random false-positive failures, like #2187. |
…-grad-2F1-precision
…b.com/stan-dev/math into feature/issue-2158-grad-2F1-precision
Jenkins Console Log Machine informationProductName: Mac OS X ProductVersion: 10.11.6 BuildVersion: 15G22010CPU: G++: Clang: |
|
@bbbales2 you have a "change request" pending on this but I just read the current code again and I think it does look right now. The tests are passing so I think this is good to merge. Can you take a look again and I'd love to get it merged. |
|
@sakrejda in! |
Summary
Fixes issue #2158 by introducing combined absolute/relative tolerance to grad_2F1 and increasing the max number of default steps. Overall it should terminate sooner for very large gradients, the increased max number of steps is required to get the stated precision in the tests. New tests are referenced vs. a quick discrete diff in Mathematica.
The absolute/relative tolerance calculation is borrowed from the unit test formulation of the same problem:
Tests
Test case uses values from a user's failing model and reference values from mathematica. The gradients are in the tens of millions and the original absolute tolerance calculation was excessive there.
Side Effects
Nothing noticeable. Most likely improved performance for models that hit the boundaries in neg_binomial_cdf and anything else that relies on grad_2F1. Fixes a user-reported model bug as a side-effect (linked above).
Release notes
Make gradients for negative binomial and 2F1 function more robust for boundary values.
Checklist
Math issue User-reported that grad_2F1 bumps into internal counter limit #2158
Copyright holder: Krzysztof Sakrejda
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
./runTests.py test/unit)make test-headers)make test-math-dependencies)make doxygen)make cpplint)the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested