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
|
Fyi our backend CI is migrating to a new server cluster so the jenkins will hiccup for a minute while we get everything settled |
|
It looks like your changes exposed an issue in the opencl code. Let me look into this |
|
@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 ./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 |
|
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 |
|
I'm not sure why those tests were passing before, as
Does anybody know the history of the Hessian test? The |
|
On a similar note, would it also be a good idea to extend the |
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 The main change is 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 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 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. |
|
@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 |
|
@SteveBronder: Yes, I just added your branch! |
|
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 |
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.
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 |
… 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 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 |
2 participants
Summary
This pull request addresses #1985 by updating
normal_lccdfandstd_normal_lccdfto callnormal_lcdfandstd_normal_lcdfwith the signs ofxandmureversed.Tests
The tests have been updated to use the latest values of
pnormin R 4.6.1. These are more direct tests than the previous versions, which worked via anexpcall, somewhat defeating the purpose of having separatelcdfandlccdffunctions.Side Effects
None.
Release notes
normal_lccdfandstd_normal_lccdfupdated to usenormal_lcdfandstd_normal_lcdfimplementationsChecklist
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
./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