Aligned normal_lccdf and std_normal_lccdf with the lcdf functions - #3363
Aligned normal_lccdf and std_normal_lccdf with the lcdf functions#3363jaburgoyne wants to merge 7 commits into
Conversation
Sorry, something went wrong.
|
Fyi our backend CI is migrating to a new server cluster so the jenkins will hiccup for a minute while we get everything settled |
Sorry, something went wrong.
|
It looks like your changes exposed an issue in the opencl code. Let me look into this |
Sorry, something went wrong.
|
@jaburgoyne this is interesting, for the tests it looks like there may be an autodiff issue with using the symmetry trick here (which implies there is something wrong with our code that your code is showing us). At y=-50 we are getting an error for the hessian calculation vs finite difference. ./test/unit/math/expect_near_rel.hpp:90: Failure
Value of: stan::math::is_nan(x1) && stan::math::is_nan(x2)
Actual: false
Expected: true
expect_near_rel(0, nan)
expect_near_rel; require items x1(0, 0) = x2(0, 0): hessian() Hessian for finite_diff vs fvar<var>
./test/unit/math/expect_near_rel.hpp:90: Failure
Value of: stan::math::is_nan(x1) && stan::math::is_nan(x2)
Actual: false
Expected: true
expect_near_rel(0, nan)
expect_near_rel; require items x1(0, 0) = x2(0, 0): hessian_fvar() Hessian for finite_diff vs fvar<fvar<double>>
[ FAILED ] AgradRev.mathMixScalFun_normal_lccdf (1 ms)This passed at y=50 previously and your code changes test the negative boundary which is why the error shows up now. So this issue existed before your changes. I need to dive into this a bit more, would you mind taking a look as well? Claude says it has a solution that looks simple'ish and I'll make a little branch for it |
Sorry, something went wrong.
|
Thanks for the insight and interest, @SteveBronder ! I saw the errors and wanted to have a look but haven't yet had a chance with post-conference catch-up. I think I can make some time tomorrow to try to replicate the error on my machine and start looking at solutions. I had actually already extended the range of the test beyond what is in this pull request – the current std_normal_lcdf actually matches R all the way out past -1e8 – but it is interesting to look at where and why autodiff would break down. |
Sorry, something went wrong.
|
I'm not sure why those tests were passing before, as y=50 is well into the range where even the LCCDF underflows to 0. I've pushed some new commits that do three things:
Does anybody know the history of the Hessian test? The normal_lccdf_log_test was the only one of the four that had it. It fails on my machine for the LCDFs but works for LCCDFs. Before merging, it would be good to know whether it should be there. Removing it altogether would be the most consistent with the codebase as it stood before this pull request...but perhaps we should be testing the Hessians? |
Sorry, something went wrong.
|
On a similar note, would it also be a good idea to extend the expect_ad tests in mix to cover a larger range, as I did for the tests in prim? I'm not sure where the development community lands on speed vs. completeness for the test suite. |
Sorry, something went wrong.
PR #3363 rewrote normal_lccdf as normal_lcdf(-y, -mu, sigma), which made every lccdf error report the wrong function: normal_lccdf(1.0, 0.0, -1.0) -> "normal_lcdf: Scale parameter is -1, but must be positive!" Thread the reporting name through as a leading defaulted non-type template parameter, so the reflected distributions pass their own name and the checks still run exactly once. Same shape as the existing opencl_code<const char*> kernel-string parameter in stan/math/opencl/kernel_generator/opencl_code.hpp. Putting the parameter first keeps T_y/T_loc/T_scale deducible, at the cost of breaking explicit-type-argument calls such as normal_lcdf<double, double, double>(y, mu, sigma). That idiom should not compile anyway; the two in-repo users were the deprecated *_cdf_log wrappers, which never needed the explicit arguments. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
PR #3363 rewrote the prim lccdf functions as lcdf(-y, -mu, sigma) but left the OpenCL kernels on the old erf-based algorithm, which hard-truncates the upper tail: select(scaled_diff > 8.25 * INV_SQRT_TWO, 0.0, 1.0 - erf(scaled_diff)) log(0) is -inf, so for (y - mu) / sigma > 8.25 OpenCL returned -inf while the new CPU path returned a finite value. That is the ProbDistributionsNormalLccdf CI failure: OpenCL -inf vs CPU -136.16378053699881 at z = (2, 16, -0.01). Three cases in test/unit/math/opencl/rev/normal_lccdf_test.cpp were affected (opencl_matches_cpu_small, opencl_broadcast_y, opencl_matches_cpu_big). The std_normal_lccdf tests passed only because their inputs stay in [-1, 1]; the same 8.25 cutoff was latent there. Delegating fixes it by construction: the OpenCL normal_lcdf kernel is a line-by-line transliteration of the prim implementation, so CPU/OpenCL agreement for lccdf reduces to the already-passing normal_lcdf comparison. Verified by emulating the patched kernels on CPU: values and all three gradients match prim to machine precision on every previously-failing case (0/459 gradients outside 1e-8 on the N=153 case, worst relative diff 2.2e-16). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
PR #3363 rewrites normal_lccdf as normal_lcdf(-y, -mu, sigma), which routes lccdf into argument regions of lcdf that lccdf never touched before and that lcdf's own mix test barely probes -- it stops at y = +10. Three distinct defects live in those regions, all pre-existing on develop. These tests are RED. Each fix lands in its own commit after this one. P1 exp(x2) sits in a denominator in the scaled_diff > 2.9 derivative branch. Autodiff of 1/w needs w^(k+1) representable, so the k-th derivative dies at x2 > 709.78/(k+1). Measured, at mu=0 sigma=1: d2 fvar<fvar<double>> silent 0 for y in [26.6, 37.5], NaN y >= 37.6 d2 fvar<var> NaN y >= 37.7 d3 fvar<fvar<fvar<double>>> silent 0 at y = 20, NaN y >= 22 d3 fvar<fvar<var>> NaN y >= 37.6 The silent zeros are ~1e-158 and numerically harmless; the NaNs poison whole Hessians. Note the silent-zero band is invisible to expect_ad, so it needs a direct fwd assertion. P2 erfc(|s|)^2 underflows inside log(erfc(-s)), giving a NaN third-order mixed derivative. Reachable ONLY in fvar<fvar<var>>, i.e. expect_ad's grad_hessian; fvar<var> and fvar<fvar<double>> are both clean. Window bisects to s in (-20, -19.2103473480), i.e. y in (-28.2842712475, -27.1675337575). The existing y = -20*sqrt(2) test point sits exactly on the lower edge and takes the other branch, missing it. P3 The scaled_diff < -29 residual correction adds 0.0015065154280332*x2, which grows quadratically, but d/ds log Phi grows linearly (-> 2|s|, DLMF 7.12.1). A quadratic fit cannot track a linear asymptote. Measured relative gradient error: 1.00e-4 at y=-54.46, 1.00e-3 at -66.47, 8.2e-3 at -100, 1.00e-2 at -106.06, 9.6e-2 at -300. This corrupts plain var gradients, i.e. HMC transitions, not just Hessians. Tests are placed in the tier that exhibits the defect: fwd for the pure forward orders, rev for the first-order gradient, mix for expect_ad. Reference values come from mpmath at 60 significant digits; the generating snippet is committed alongside the tables. EXPECT_FLOAT_EQ is deliberately avoided for the tail derivatives -- those magnitudes are below the float minimum, so it would compare 0 == 0 and pass vacuously. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The scaled_diff < -29 branch took the Abramowitz & Stegun negative-tail approximation and added a cubic residual fit whose leading term is 0.0015065154280332 * x2. That term grows quadratically in scaled_diff. DLMF 7.12.1 gives, for y = -scaled_diff > 0, erfc(y) ~ (exp(-y^2) / (sqrt(pi) y)) * sum_m (-1)^m (1/2)_m / y^(2m) and since dncdf_log = (2/sqrt(pi)) exp(-y^2) / erfc(y), this is dncdf_log ~ 2y / (1 - 1/(2y^2) + 3/(4y^4) - 15/(8y^6) + ...) so dncdf_log grows LINEARLY, as -2*scaled_diff. A quadratic residual fit cannot track a linear asymptote; it must eventually dominate. Measured relative error of the old branch against a 60-digit mpmath reference: scaled_diff -29 -40 -75 -100 -1000 old 9.5e-06 1.9e-04 1.0e-02 2.3e-02 6.8e-01 new 1.3e-11 1.0e-12 6.5e-15 7.0e-16 7.8e-18 At scaled_diff = -1000 the old branch is 68% wrong. In y terms (mu=0, sigma=1) the gradient passes 1e-4 relative error at y = -54.46 and 1e-2 at y = -106.06. This is a plain first-order var gradient, so it corrupts HMC transitions, not just Hessians -- and PR #3363's reflection newly exposes it on the lccdf side at large positive y. Replacing the branch with the truncated asymptotic is also exponential-free: no erf, no exp, nothing that can overflow. DLMF states the remainder is bounded by the first neglected term for real argument, which at scaled_diff = -29 is 6.5625/29^8 = 1.3e-11, matching the measured error. The seam at scaled_diff = -29 moves from 3.3e-06 to 6.2e-06 relative, still well inside the existing inter-branch jumps (~2.7e-5 at scaled_diff = 2.9). Also corrects the sqrt(2) scaling in the new rev sweep assertions: dncdf_log is d/d(scaled_diff), while the reported partial is dncdf_log/(sigma*sqrt(2)). The discrete mpmath-referenced cases were unaffected. Verified the corrected sweep is still red without this fix and green with it. All P1, P2 and P3 tests now pass across prim, rev, fwd and mix. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
PR stan-dev#3363 rewrites normal_lccdf as normal_lcdf(-y, -mu, sigma), which routes lccdf into argument regions of lcdf that lccdf never touched before and that lcdf's own mix test barely probed -- it stopped at y = +10. Three defects live in those regions, all pre-existing on develop. Verified with git diff that normal_lcdf.hpp was untouched by stan-dev#3363: the PR reaches these, it does not introduce them. Notation: s = scaled_diff, z = (y - mu) / sigma = s * sqrt(2). P1 exp(x2) sat in a denominator in the s > 2.9 derivative branch. Autodiff of 1/w needs w^(k+1) representable, so the k-th derivative dies at x2 > 709.78/(k+1). Measured at mu=0 sigma=1: d2 silently 0 for y in [26.6, 37.5] then NaN from 37.6; d3 silently 0 at y = 20, NaN from 22. Abramowitz & Stegun (1964) 7.1.26 publishes exp(-x^2) as a NUMERATOR factor; Stan had algebraically inverted it. Restoring A&S's arrangement is an exact identity -- max 6 ulp over a 200,001-point scan of s in [2.9, 26.64] -- and lets the exponential underflow instead of overflowing. R's pnorm does the same, forming the Gaussian factor only in the numerator (do_del) and guarding underflow, never overflow. P2 erfc(|s|)^2 underflowed inside log(erfc(-s)), giving a NaN third-order mixed derivative. Reachable ONLY in fvar<fvar<var>>, i.e. expect_ad's grad_hessian; fvar<var> and fvar<fvar<double>> are both clean. Window bisects to s in (-20.0000000000, -19.2103473480). Moved the crossover to R pnorm's own cutoff: it switches to the Cody tail form at y > M_SQRT_32, i.e. |x| > sqrt(32), which in scaled_diff is exactly sqrt(32)/sqrt(2) = 4. Our Cody coefficient set differs from R's p[]/q[] (ours is in s, theirs in x), so its range was measured separately: <= 9.593e-17 relative down to s = -4, degrading past -3.52. P3 The s < -29 residual correction added 0.0015065154280332 * x2, growing quadratically, but d/ds log Phi grows linearly, to 2|s| (DLMF 7.12.1). A quadratic fit cannot track a linear asymptote. Relative gradient error was 1.0e-4 at y = -54.46, 8.2e-3 at -100, and 68% at s = -1000. This corrupted plain var gradients, i.e. HMC transitions. Replaced with the truncated asymptotic Mills ratio, which needs no exponential at all: 1.3e-11 at s = -29 improving to 7.8e-18 at s = -1000. All three land at four sites together -- prim normal_lcdf, prim std_normal_lcdf, and the two OpenCL kernels, which are line-by-line transliterations whose results are compared against CPU in opencl/rev. Also in this change: - OpenCL normal_lccdf and std_normal_lccdf now reflect through lcdf like prim does. They had been left on the old erf algorithm, whose select(scaled_diff > 8.25 * INV_SQRT_TWO, 0.0, ...) makes log(0) = -inf. That is the reported CI failure: OpenCL -inf against CPU -136.16378053699881 at z = (2, 16, -0.01). - A const char* template parameter carries the caller's name into the error checks, so normal_lccdf(1, 0, -1) reports normal_lccdf rather than normal_lcdf. This breaks normal_lcdf<double, double, double>(...), an idiom that should not compile anyway; the two in-repo users were the deprecated *_cdf_log wrappers. - normal_lcdf's doc block records what each cutoff rests on: A&S 7.1.26, Cody (1969) and DLMF 7.12.1 with links, the R pnorm and SciPy log_ndtr cross-references, and the stan-dev#1411 origin of the interior Taylor cutoffs, which are original and undocumented -- the PR describes them as "original Taylor expansions ... to bridge the gap where the numerical approximations are unstable". std_normal_lcdf points at it rather than duplicating; the duplicate had already drifted. - Doxygen completed on both: \ingroup prob_dists, @return, and removal of @tparam T_loc, @tparam T_scale, @PARAM mu and @PARAM sigma from std_normal_lcdf, which does not take them. Tests, all in prim and mix: - prim: 16 rows of log Phi against references carried to 40 significant digits in comments beside the correctly-rounded doubles, spanning all three value branches out to s = -212. Mutation-checked: a 0.15% perturbation of temp_p's leading coefficient fires all six Cody rows. One row is asserted loosely and documented -- at s = 11.31 this platform's libm erfc is itself 8.0e-06 out, which passes straight through. - mix: expect_ad at every branch cutoff, bracketed rather than landing on the seam, since finite_diff_grad_hessian_auto's stencil would straddle it. - mix: per-branch derivative accuracy at each branch's own measured error, because expect_ad's blanket 1e-4 leaves the (0.8, 1.5] branch only 1.6x. - mix: order-2 and order-3 tail derivatives against 60-digit references, plus a finiteness sweep over y in [-120, 120]. These are asserted directly rather than through expect_ad because expect_ad cannot see either failure: the silent zeros are ~1e-158, invisible to a finite-difference comparison whose relative tolerance floors its denominator at 1, and test_ad.hpp stops at fvar<fvar<var>>, so fvar<fvar<fvar<double>>> is never instantiated. OpenCL cannot be executed on the development machine -- the only local device has no fp64 -- so those two kernels were verified by transliterating them back to C++ and diffing against prim over 16,000 points in y in [-400, 400]: value bit-identical, gradient worst 3.6e-16. Jenkins remains the real check. Known and deliberately untouched: the s > 2.9 branch is ~6e-5 relatively inaccurate because it uses P(t)/t where A&S 7.1.26 requires P(t)/2. This change preserves that bit-for-bit and makes no accuracy claim. Fixing it means replacing all five piecewise Taylor branches with INV_SQRT_PI * exp(-x2) / (1 - 0.5 * erfc(scaled_diff)), reusing the erfc the value branch already computes, which would shift many gradient expectations. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The four mix/prob normal tests had grown four copies of the same thing: the same d2 and d3 extractors, the same branch_ref/tail_ref/far_ref structs, the same cutoff array, and the same reference tables differing only by a sign. Adds test/unit/math/mix/prob/normal_lcdf_tail_test_helpers.hpp holding the references once and the checks as functors: d2, d3 pure-forward derivative extractors value_and_grad reverse-mode value and gradient mills_dlogphi_dy five-term DLMF 7.12.1 asymptotic expect_ad_across_cutoffs expect_ad bracketing every branch cutoff expect_ad_at_defect_inputs expect_ad at the three defects' inputs expect_branch_accuracy per-branch accuracy vs 60-digit references expect_tail_derivatives order-2 and order-3 tail references expect_derivatives_finite NaN sweep over y in [-120, 120] expect_far_tail_gradient far-tail references plus a Mills sweep All four functions share one piecewise implementation, so one set of references covers them. They are stated for normal_lcdf(y, 0, 1); a reflected function -- lccdf(y) == lcdf(-y) -- is exercised at -y with odd-order derivatives flipped. That is the `dir` argument, spelled orientation::lcdf or orientation::lccdf at each call site, which is also what makes the reflection explicit rather than buried in four hand-negated tables. Each test file drops to roughly 50 lines: its own expect_ad points, then six one-line calls. The helper is 260 lines, against about 700 of duplication it replaces. Mutation-checked after the refactor, so the move did not weaken anything: reverting the exp(-x2) rearrangement fails defect_inputs, tail_derivatives (naming y = 27, 30, 33, 37 for order 2 and 20, 22 for order 3) and derivatives_are_finite; re-centring the 1.85 series on 1.80 fails branch_cutoffs and branch_accuracy. Same failures as before the refactor. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
runTests.py --jumbo concatenates test file bodies into one .cpp, hoisting and
de-duplicating the includes (generateJumboTests). Anonymous namespaces from
every file in a folder therefore land in the same translation unit, so the
`lcdf`, `lccdf` and `dir` objects I put in `namespace { ... }` collided.
Reproduced by concatenating the four files the way generateJumboTests does:
error: redefinition of 'dir' (x3)
error: redefinition of 'lcdf'
error: redefinition of 'lccdf'
Each file now uses a namespace named after itself -- normal_lcdf_mix_test,
normal_lccdf_mix_test, std_normal_lcdf_mix_test, std_normal_lccdf_mix_test --
holding `fn` and `dir`, matching the dirichlet_test convention already in this
folder. The same concatenation now compiles clean.
Checked the prim tests for the same hazard: their reference structs are
declared inside the TEST bodies, so they are function-local and a prim jumbo
concatenates cleanly as-is.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
So in general the hessian tests are checking a finite diff for the hessian vs. Stan's AD computation of the Hessian. Does that answer your question? I was also really curious about where the cutoff points for the Taylor expansions in the adjoint calculations in normal_lcdf came from. I had claude dig backwards and found the original PR where those were written and the answer is, "Just made sense" lol. Which I think they do! They were written so that at the cutoff points we would have at least ~1e-4 relative tolerance for gradient and hessian calculations. I have a branch below that adds a lot of doc to the normal_lcdf to explain the cutoff points for the cdf and adjoint calculations. It also aligns the cdf calculations cutoff points to be more in line with R's pnorm. The main change is stan/math/prim/prob/normal_lcdf.hpp. There I added a bunch of doc explaining the cutoffs, adjusted one of the cutoffs for the cdf to be in line with R's pnorm, and fixed an overflow issue in one of the adjoint calculations that was causing the hessian failure. The tests in that are mostly to make sure we catch the hessian and gradient calculations correctly. The opencl changes are just to make sure the opencl code stays in parity with the cpu code. You can check it out locally via the below (assuming your git remote for stan-dev is named upstream) git fetch upstream
git switch -c normal-lcdf-tails upstream/fix/normal-lcdf-tails # or: git reset --hard upstream/... if he already had it
git log --oneline fix/normal_lccdf..normal-lcdf-tails # 15 commitsIf you are cool with the changes you can merge and push them to your branch via git switch fix/normal_lccdf git merge --ff-only upstream/fix/normal-lcdf-tails fyi I had claude help a good bit with this, but the docs show provenance for the actual code changes that I think make a lot of sense / align with R's pnorm. |
Sorry, something went wrong.
|
@jaburgoyne just a sanity check, did you add my branch to this branch? If not I must have fat fingered something and I can remove / revert those changes off of this branch |
Sorry, something went wrong.
|
@SteveBronder: Yes, I just added your branch! |
Sorry, something went wrong.
|
The extra documentation is a great idea! Thank you for fixing up the OpenCL bits – I had not originally realised that I needed to look there, too. I don't completely understand all of the new tests you added, but assuming they do what they say, they are certainly more rigorous than what was there before. As for the Hessian tests – my question wasn't what they were doing, but rather why they had been removed everywhere except normal_ccdf_log_test. As written, they pass for normal_lccdf and std_normal_lccdf but fail for normal_lcdf and std_normal_lcdf. Should you, I, or Claude try to make a cleaner version? |
Sorry, something went wrong.
Those tests are purely to try to 'eek out the poor hessian calculation by making sure we have better coverage on each of the cutoff points. expect_ad_across_cutoffs goes +/- 0.01 across each cutoff point to check that each one does the adjoint calculations correctly
That is a great question... maybe the first questions is whether they existed 😅 Looking at git blame for the mix tests it looks like normal_lcdf never had these checks. So it is good we are catching it now https://github.com/stan-dev/math/blame/develop/test/unit/math/mix/prob/normal_cdf_log_test.cpp#L12 |
Sorry, something went wrong.
… for the normal distribution
Right! So this is actually a much better test than what was in there before, anyway. I apologise for being so slow to piece this together – newbie growing pains! – but I think that I finally see the big picture.
Given the other tests that you added, let's just remove these. I pushed one more commit taking those Hessian tests out, and expanding the other tests slightly to:
Claude caught that second point, and so far as I can see, it was right: the suite was never checking whether autodiff on mu and sigma were working. Given the way that this pull request has developed, perhaps one other question that I originally put into #3364 but now fits better here: is it good to keep std_normal_lcdf and normal_lcdf as independent implementations, complicated branching and all, or is it worth it in this pull request (or a new one) to try to make normal_lcdf call std_normal_lcdf, with the necessary corrections for the presence of mu and sigma? I'm happy to have a go at that and now have a better sense of how to make sure the results fits better into the whole test apparatus. |
Sorry, something went wrong.
|
This PR has gotten a little big now that we found things that needed fixed. Let's merge this and do a separate PR for the normal -> std_normal |
Sorry, something went wrong.
Jenkins Console Log Machine information No LSB modules are available. Distributor ID: Ubuntu Description: Ubuntu 20.04.3 LTS Release: 20.04 Codename: focal CPU: G++: Clang: |
Sorry, something went wrong.
SteveBronder
left a comment
There was a problem hiding this comment.
Looks good to me! You should have an email that invites you as a collaborator. Pressing the accept link in that email should let you press the merge button in this PR :)
Sorry, something went wrong.
|
Thanks – great that this is ready to merge! I accepted the invitation as a collaborator, but merging is still blocked for me. Do I perhaps need to be a contributor in order to have permission do that? |
Sorry, something went wrong.
Summary
This pull request addresses #1985 by updating normal_lccdf and std_normal_lccdf to call normal_lcdf and std_normal_lcdf with the signs of x and mu reversed.
Tests
The tests have been updated to use the latest values of pnorm in R 4.6.1. These are more direct tests than the previous versions, which worked via an exp call, somewhat defeating the purpose of having separate lcdf and lccdf functions.
Side Effects
None.
Release notes
normal_lccdf and std_normal_lccdf updated to use normal_lcdf and std_normal_lcdf implementations
Checklist
Copyright holder: John Ashley Burgoyne
The copyright holder is typically you or your assignee, such as a university or company. 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
the code is written in idiomatic C++ and changes are documented in the doxygen
the new changes are tested