[avehtari]

Fixes #3153. Fixes #1578. Addresses #2508 partially.

This PR improves beta_binomial_lpmf, beta_binomial_cdf, beta_binomial_lcdf, beta_binomial_lccdf, beta_neg_binomial_lpmf, beta_neg_binomial_cdf, beta_neg_binomial_lcdf, beta_neg_binomial_lccdf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf, dirichlet_multinomial_lpmf, lkj_corr_lpdf, lkj_corr_cholesky_lpdf, lkj_cov_lpdf, student_t_lpdf, multi_student_t_lpdf, multi_student_t_cholesky_lpdf, yule_simon_lpmf, yule_simon_cdf, yule_simon_lcdf, yule_simon_lccdf, and the derivatives of lbeta.

Claude assisted with this PR.

Summary

Several distributions form a difference f(x + k) - f(x). Here f is lgamma, lbeta or digamma, x is a shape, concentration, dispersion or degrees-of-freedom parameter, and k is a count or an offset. When x is large, the two terms are large and almost equal, and the difference keeps no correct digits. In the values the terms have the size of x, so the error is about eps * x. In the gradients the terms are about log(x), and a model with the parameter on the log scale multiplies the error by x.

In the model of #3153, beta_binomial_lpmf(57 | 117, s, s) with s = exp(lc) / 2, develop returns a log probability 81 too high from lc = 42, and its gradient in lc is 0 or noise from lc = 24. Warmup then finds a flat plateau and cannot leave it. With this PR the value is within 2.3e-14 of the exact value and the gradient within 3.7e-14, for lc from 0 to 44.

lc	develop: diff to binomial limit	this PR	exact	develop d/dlc	this PR	exact
20	-4.4e-08	-1.1130e-07	-1.1130e-07	1.7e-06	1.1130e-07	1.1130e-07
28	-4.6e-05	-3.7309e-11	-3.7338e-11	0	3.7357e-11	3.7338e-11
32	+4.5e-03	-6.5937e-13	-6.8386e-13	-0.14	6.8570e-13	6.8386e-13
42	+81.098	+2.3e-14	-3.1e-17	0	0	3.1e-17

From lc = 34, the exact values are below the rounding error of any form (about 2e-14), because 57 log p + 60 log q is about -81 there.

#2508 reports the same defect: beta_binomial_lpmf(500 | 1000, s, s) for s = 0.5e10, 0.5e13, 0.5e19. Develop is wrong by 5.3e-8, 2.8e-4 and 693 there, and returns +689.47 at the largest s. This PR fixes all three points, and they are rows of the beta_binomial_lpmf test. This PR does not close #2508, because the comment there about beta_binomial_lcdf(3 | 10, 298.3, 146.75) (NaN, exact -3.93335) has a different cause: hypergeometric_3F2 does not stop at the zero term of the terminating series 3F2({1, 302.3, -6}; {5, -151.75}; 1) (exact 17.63), and a later 0 * inf gives NaN. #3437 repairs that and closes #2508.

#1578 reports wrong derivatives of neg_binomial_lpmf for large alpha. Its example (n = 0, alpha = 1e10 + 1, beta = 0.8) is already correct on develop, because develop no longer uses the Poisson approximation that the issue describes. This PR repairs the loss that remains in d/dalpha = digamma(alpha + n) - digamma(alpha) + log(beta / (1 + beta)). The loss matters when the mean alpha / beta stays fixed: the error of alpha * d/dalpha on develop reaches 3.0 at alpha = 1e15, and with this PR it is 1e-14 (rows of the new large-shape test).

The same root cause is in more functions. This PR repairs:

function	develop, worst error at large shape	this PR
beta_binomial_lpmf	value 81, log-shape gradient 6e3	rounding level
beta_neg_binomial_lpmf	value 10, gradient 13	rounding level
neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf	gradient 3, 27, 26	1e-14, 2e-14, 2e-13
neg_binomial_2_log_glm_lpmf	value 23, gradient 7.1	4e-13, 6e-14
dirichlet_multinomial_lpmf	gradient 2.8	8e-15
lkj_corr_lpdf, lkj_corr_cholesky_lpdf	value 3.9, gradient 1.5	4e-15, 2e-16
yule_simon_lccdf (and _cdf, _lcdf in relative terms)	value 7.6, gradient 2.1	2e-13, 8e-15
beta_binomial_{cdf,lcdf,lccdf}	gradient 22, lcdf inf or NaN	limited by the 3F2 series (1e-7)
beta_neg_binomial_{cdf,lccdf}	gradient 15	limited by the 3F2 series (2e-6)
student_t_lpdf	value 0.16, nu gradient 0.14	1e-11, 5e-12
multi_student_t_lpdf, multi_student_t_cholesky_lpdf	value 2e-7, nu gradient 2e-7	8e-15, 5e-15
yule_simon_lpmf	value 2e-7, gradient 3e-7	3e-14, 2e-15
lbeta (var, fvar), derivative	log-scale derivative 0.63	1e-14

The shape is 1e15 for the first nine rows and at most 1e15 for the last four (nu decades for student_t_lpdf). "Log-shape gradient" is x * d/dx, the gradient that a sampler sees when the model parameter is log(x). All references are mpmath at 80 to 140 digits, and two independent routes agree: loggamma for the values, closed-form digamma differences and mp.diff for the gradients, and sums of the pmf for the cdfs.

The four negative binomial functions also had two defects at other arguments: the value lost phi at small phi, and the GLM partials overflowed for a linear predictor above 709.8. This PR repairs both (items 5 and 6 of "Related defects fixed on the way").

propto. In some of the new forms, a constant that propto drops cancels inside one expression against terms that depend on a parameter, so the code cannot skip it. For example, lbeta(n, alpha + 1) contains lgamma(n). To keep the propto behavior, the code adds such constants back with the opposite sign: lgamma(n) in yule_simon_lpmf, and lgamma(y + 1), phi log(phi) (data phi) and y theta (data x, alpha and beta) in neg_binomial_2_log_glm_lpmf, also in OpenCL. For the GLM with data phi, this brings back a rounding error of about eps * phi log(phi) in the propto value. A probe of all 14 changed densities, with every combination of data and parameter arguments, matches the propto and full values of develop at all 456 points on the current head, and a new test checks full - propto against the dropped terms.

Most of the added source lines are in six files:

file	added lines
prim/fun/log_rising_factorial_ratio.hpp (new helpers)	276
prim/fun/digamma_diff.hpp (new)	199
OpenCL device function log_beta_ratio (new)	168
OpenCL device function digamma_diff (new)	115
opencl/kernel_generator/elt_function_cl.hpp (registrations)	104
student_t_lpdf (series), CPU and OpenCL	89 + 41

Each of the other 26 distribution files (CPU, OpenCL and the GLM kernel) adds 3 to 61 lines.

What this PR adds

digamma_diff(x, d) = digamma(x + d) - digamma(x), a new function in prim/fun. It is public, but not exposed in the Stan language, because 19 distribution files and the derivatives of lbeta use it. log_rising_factorial can use it in a follow-up. It has three paths:

• a non-autodiff integer d from 0 to 8 (a count): sum_{j < d} 1 / (x + j). All terms are positive, and d = 1 gives exactly 1 / x;
• x < 10 and d >= 10: the plain difference, which loses only a few ulp there (the result is at least 0.72 and digamma(x) <= 2.25);
• otherwise: the recurrence shifts both arguments to at least 10, then log1p(d / y) plus the asymptotic series, with the differences y^-2i - (y + d)^-2i formed as (u - v) h_i so that nothing cancels.

Over 7 281 points, x from 1e-3 to 1e18 and d from 1e-6 to 1e7, the error is at most 4.2 ulp; the plain difference reaches 6.5e16 ulp. The function is vectorized and works with all autodiff types through generic arithmetic, with expect_ad tests. The PR also adds an OpenCL device function with the same paths and a kernel-generator registration.

internal::log_rising_factorial_ratio, internal::log_beta_ratio and their helpers in prim/fun/log_rising_factorial_ratio.hpp. They write each lgamma through lgamma_stirling_diff, so that the terms linear in the arguments cancel exactly, and they form the rest without cancellation: the x log1p(k / x) - k parts are split off for k < x, and the integer parts are added exactly. log_beta_ratio(alpha, beta, n, m) = lbeta(alpha + n, beta + m) - lbeta(alpha, beta) selects this form or the plain lbeta form by the size of their terms. The Stirling form has rounding errors of the size of n + m, and the lbeta form of the size of the smaller shape. Without the selection, points with N = 1e6 and one small shape lost up to 4e-11 where develop was correct.

The rev and fwd derivatives of lbeta. The partial in an argument of at least 10 is now -digamma_diff(a, b). Below 10 it stays digamma(a) - digamma(a + b), and the two partials share digamma(a + b) as before. lbeta itself is accurate at large arguments, but its derivative had the same cancellation. So every model that calls lbeta with a large parameter had inaccurate gradients, and so did multi_student_t_lpdf, which gets its nu gradient by autodiff through lbeta. Below 10 the plain difference has an absolute error of a few eps, which stays at the rounding level on the log scale too.

Changes per function

• beta_binomial_lpmf: value by log_beta_ratio, partials by digamma_diff.
• beta_neg_binomial_lpmf: value as -log(n) - lbeta(n, beta) + log((r)_n / (r + alpha + beta)_n) + log((alpha)_r / (alpha + beta)_r), partials by digamma_diff.
• neg_binomial_lpmf: alpha partial by digamma_diff. The beta partial term alpha / beta - alpha / (1 + beta) is now alpha / beta / (1 + beta). The value term lchoose(n + alpha - 1, n) is now -lbeta(alpha, n + 1) - log(n + alpha).
• neg_binomial_2_lpmf, neg_binomial_2_log_lpmf: phi partial by digamma_diff. The value term lchoose(n + phi - 1, n) is now -lbeta(phi, n + 1) - log(n + phi). When phi is autodiff, neg_binomial_2_lpmf evaluates log(mu + phi) once: the value and the phi partial both read it, and the lazy expression computed the log twice per element. The results are bitwise identical.
• neg_binomial_2_log_glm_lpmf: value per instance as -lbeta(phi, y + 1) - log(y + phi) - y log1p_exp(log(phi) - theta) - phi log1p_exp(theta - log(phi)), the form of neg_binomial_2_log_lpmf, instead of phi log(phi) - lgamma(phi) + lgamma(y + phi) - (y + phi) log(mu + phi) + y theta. The partials do not form exp(theta). With u = theta - log(phi) and e = exp(-|u|), mu / (mu + phi) is 1 / (1 + e) for u > 0 and e / (1 + e) otherwise, and phi / (mu + phi) is the other of the two. The theta partial is y - (y + phi) mu / (mu + phi), and the phi partial is mu / (mu + phi) - (y / phi) phi / (mu + phi) - log1p_exp(u) + digamma_diff(phi, y) (formulas from @jachymb's review). The value and the partials share e and log1p(e): log1p_exp(±u) = max(±u, 0) + log1p(e). With a scalar phi and a vector y, the terms that depend on y only (lchoose, lgamma(y + 1) under propto, digamma_diff(phi, y)) are evaluated once per distinct count, through internal::per_count_values, if the largest count is at most 4 N + 64. This is in the CPU version only.
• dirichlet_multinomial_lpmf: partials by digamma_diff.
• lkj_corr_lpdf (and lkj_corr_cholesky_lpdf, lkj_cov_lpdf through do_lkj_constant): the constant is a sum of lgamma(x + k/2) - lgamma(x). Each term is now lgamma(k/2) - lbeta(k/2, x), and the code sets the eta partial explicitly with digamma_diff.
• yule_simon_{cdf,lcdf,lccdf}: lgamma(alpha + 1) + lgamma(n + 1) - lgamma(n + alpha + 1) is now log(n) + lbeta(n, alpha + 1); partial by digamma_diff.
• yule_simon_lpmf: value log(alpha) + lbeta(n, alpha + 1); partial 1 / alpha - digamma_diff(alpha + 1, n).
• beta_binomial_{cdf,lcdf,lccdf}: the leading lbeta(nu, mu) - lbeta(alpha, beta) by log_beta_ratio, the digamma differences by digamma_diff.
• beta_neg_binomial_{cdf,lcdf,lccdf}: the leading term is the log pmf at n + 1, now computed as in beta_neg_binomial_lpmf; partials by digamma_diff.
• student_t_lpdf: unchanged for nu < 40. For nu >= 40, the terms lgamma((nu + 1) / 2) - lgamma(nu / 2) - log(nu) / 2 and digamma((nu + 1) / 2) - digamma(nu / 2) come from the asymptotic series in h = nu / 2, -log(2) / 2 - 1/(8h) + 1/(192h^3) - 1/(640h^5) + 17/(14336h^7) - 31/(18432h^9), and its derivative. Against mpmath, from h = 20 to 1e15, the errors are at most 0.9 and 2.4 eps. The series is also cheaper than two lgamma or two digamma calls.
• multi_student_t_lpdf, multi_student_t_cholesky_lpdf: the constant lgamma((nu + p) / 2) - lgamma(nu / 2) is now lgamma(p / 2) - lbeta(p / 2, nu / 2). The cholesky version sets the nu partial with digamma_diff(nu / 2, p / 2); the other version gets it through lbeta.
• OpenCL: beta_binomial_lpmf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf (prim file and kernel) and student_t_lpdf compute the same formulas.

Related defects fixed on the way

1. lkj_corr_lpdf dropped d constant / d eta at eta == 1.0. The special case for eta == 1 returned the constant as a number without a derivative. At eta = 1 with K = 3 and the matrix in the new test, develop returns d/deta = -0.165; the exact value is 1.2214. eta = exp(0) = 1 is the value of a positive parameter at the default initialization.
2. beta_binomial_lcdf computed log(1 - C), with C the complementary series. Where the cdf is small, this keeps no digits, and with the accurate C it can become log of a value at or below 0. Where C > 1/2, the function now uses the mirror relation: X <= n if and only if N - X > N - n - 1, and N - X has the distribution beta_binomial(N, beta, alpha). That is the same series with the shapes swapped, without the complement. On develop this function returned NaN at 16 of 81 test points and was wrong by up to 61. Now it has the accuracy of the series everywhere. For alpha = 1 exactly, the mirrored series is 0 / 0 at one term, so that case keeps log(1 - C).
3. OpenCL beta_binomial_lpmf did not check 0 <= n <= N. The code built the check but never computed its result, so the function never returned LOG_ZERO. Without propto the value was -inf by accident (through binomial_coefficient_log), but with a nonzero gradient. Under propto it was finite: -8.216 for n = 6, N = 5 in the new test. The CPU version returns LOG_ZERO without a gradient, and now the OpenCL version does too.
4. OpenCL neg_binomial_2_log_glm_lpmf with a scalar alpha and a vector phi failed. The kernel decided from alpha instead of phi whether to sum the phi derivative. It returned work-group sums where per-instance derivatives were necessary, and the gradient assignment threw a size mismatch. The decision now depends on phi only, and the kernel does not compute the phi derivative when phi is data. A new test covers both mixed cases: scalar alpha with vector phi, and vector alpha with scalar phi.
5. neg_binomial_lpmf, neg_binomial_2_lpmf and neg_binomial_2_log_lpmf lost phi at small phi. They computed lchoose(n + phi - 1, n) as binomial_coefficient_log(n + phi - 1, n), which rounds phi into n + phi - 1. The value then has an error of about eps * n / phi. At n = 1000, phi = 1e-8, the error is 7.9e-7 (5.0e-9 in neg_binomial_lpmf). At n = 1, phi = 1e-17, the value is -inf; the exact value is -39.1439. These three functions and the GLM, CPU and OpenCL, now use -lbeta(phi, n + 1) - log(n + phi), and the errors at these points are at most 5.3e-13. The partials were already correct, because they use the exact phi.
6. The neg_binomial_2_log_glm_lpmf partials overflowed for theta above 709.8. On develop the theta partial y - exp(theta) (y + phi) / (exp(theta) + phi) is NaN at theta = 720 (inf / inf) and -inf at theta = 709; the exact value is -1.5 at both points (y = 5 and 3, phi = 1.5). The e form in "Changes per function" has no overflow, and its second derivative does not cancel. At y = 1, phi = 0.07, the fvar<var> Hessian entry in theta is now -4.7225e-17 at theta = 35 (develop: 0) and -1.40777e-13 at theta = 27 (develop: -1.40970e-13).

Speed

For nu < 40, student_t_lpdf executes the same operations as develop plus one comparison (+2.4 % instructions for value and gradient). From nu = 40 it executes fewer: -52 % instructions for value and gradient at nu >= 100, and the scalar value takes 0.35 times as long.

Time ratio to develop for the other changed functions, at ordinary arguments (shapes below 10 unless stated). Median of 5 interleaved rounds on one whole node (AMD EPYC 9965, g++ 11.5, -O2); the spread between rounds is at most 2.9 %, except 4.3 % for one develop row.

kernel	ratio
neg_binomial, neg_binomial_2, neg_binomial_2_log gradient, counts 0 to 8	0.86 to 0.91
neg_binomial_2 gradient, phi from 10 to 1e4, counts 0 to 8	0.82
neg_binomial_2 gradient, counts 9 to 100	0.91
neg_binomial_2, 1000 counts from 0 to 20, one phi	1.04
beta_binomial gradient, N = 20	1.08
lbeta gradient, both arguments below 10	1.01
lbeta gradient, one argument from 10 to 1000	1.18
neg_binomial_2_log_glm_lpmf value, N = 1000, K = 5, counts with mean 3	0.28
same, gradient in alpha, beta	0.30
same, gradient in alpha, beta, scalar phi	0.30
same, with propto	0.39
same, with propto, all 1000 counts distinct	1.45
same, with propto, counts from 5000 to 5999	1.33
same, gradient in alpha, beta, vector phi	0.68
same, gradient in x	0.42

With a scalar phi, the GLM evaluates lchoose, lgamma(y + 1) and the digamma difference once per distinct count. Where the counts rarely repeat, the forms that repair the GLM (items 5 and 6) make it slower than develop.

Limits that remain, by a different cause

• beta_binomial_{cdf,lcdf,lccdf} stop their 3F2 series at a term below 1e-6. This leaves errors up to about 1e-5 in the gradients, and the functions return NaN or throw at many ordinary arguments, for example at beta = 1, at alpha + beta = 1 and at the point in beta_binomial_lpmf is unstable for large shape arguments #2508. Fix and speed up the terminating 3F2 series of the beta_binomial cdfs #3437 repairs these.
• beta_neg_binomial_{cdf,lcdf,lccdf} stop their series at a term below 1e-8. They throw from shape 1e4 and can take seconds per call. beta_neg_binomial_cdf and beta_neg_binomial_lcdf still form the complement 1 - ccdf, so they are wrong where the cdf is small. With the more accurate log pmf of this PR, 1 - ccdf rounds to exactly 0 where the cdf is below about 1e-15. beta_neg_binomial_cdf then returns 0 with NaN gradients (15 of 2352 points of a wide grid), where develop returns a wrong finite value with finite gradients. Compute the beta_neg_binomial cdfs from sums of the pmf #3441 replaces the series and the complement and repairs all of these.
• The log-shape gradients have a floor that no implementation of these functions can remove. The function returns d/dalpha and d/dbeta, and autodiff forms alpha d/dalpha + beta d/dbeta outside it. For tied shapes in numerical instability in beta_binomial_lpmf at large shape parameters #3153, the two products are -1.5 and +1.5, and their sum is 27 / exp(lc). Each partial is accurate to a few eps of its own terms, so the sum is accurate to a few eps of 1.5, not of its own size.
• The neg_binomial_2_log_lpmf Hessian has the quotient-rule cancellation of item 6 at large |eta| (NaN entries at eta = 720 and -400, n = 3, phi = 0.07). Develop has the same defect. A follow-up PR can use the e form of the GLM there too.
• log_rising_factorial is still lgamma(x + n) - lgamma(x). No Stan Math function calls it. A follow-up can use the new helpers there.

Tests

• test/unit/math/prim/fun/digamma_diff_test.cpp, test/unit/math/mix/fun/digamma_diff_test.cpp: special cases, exact finite sums, 17 mpmath values, the count path and the plain path (digamma_diff(x, 1) == 1 / x exactly from x = 1e-300 to 1e300); expect_ad at 8 points.
• test/unit/math/prim/prob/beta_binomial_test.cpp: values at the points of numerical instability in beta_binomial_lpmf at large shape parameters #3153 and beta_binomial_lpmf is unstable for large shape arguments #2508 and at large unequal shapes; the difference to the binomial limit at lc = 42.
• test/unit/math/rev/prob/beta_binomial_lpmf_test.cpp (new): log-shape gradients at 8 points; the lc gradient of the model in numerical instability in beta_binomial_lpmf at large shape parameters #3153 for lc = 16 ... 32. From lc = 34, develop returns exactly 0, which is within 1e-13 of the exact gradient, so a test there cannot see the defect.
• test/unit/math/rev/prob/large_shapes_test.cpp (new): 92 rows over the other functions, including student_t_lpdf, both multi_student_t functions, yule_simon_lpmf and lbeta. For each function, the rows are points where develop was wrong by more than 100 times the tolerance, plus one point at a small shape where develop was correct. The rows include the four negative binomial functions at phi = 1e-8 and 1e-17 and the GLM at theta = 720 and 709. The file also tests the lkj_corr gradient at eta = 1, and full - propto against the dropped terms for student_t_lpdf, yule_simon_lpmf, beta_neg_binomial_lpmf and three cases of the GLM.
• test/unit/math/rev/prob/neg_binomial_2_log_glm_lpmf_test.cpp: a vector call against the sum of single-instance calls, value and partials, with and without propto, once with repeated counts (per-count path) and once with a count above 4 N + 64.
• OpenCL: one absolute-reference test per changed distribution in test/unit/math/opencl/rev/ (six files), because the existing differential tests compare OpenCL with the CPU and cannot see an error that both share. The four negative binomial reference tests include the small-phi and large-theta rows. Also: a device-function test for digamma_diff; a test that beta_binomial_lpmf gives LOG_ZERO without a gradient for n outside 0...N, also under propto; a test of neg_binomial_2_log_glm_lpmf with a scalar alpha and a vector phi, and the reverse.

The tests find the defect. With the 23 develop headers swapped in (21 prim/prob files and the rev/fwd lbeta), all 4 new beta_binomial cases and both large-shape cases fail. The 3 old cases and the propto case pass (the propto case states the develop rule). With the branch headers restored, all tests pass again. In OpenCL, with the develop OpenCL files every new absolute-reference test fails (6 files), and so do the two new tests of related defects 3 and 4. Every old differential test passes, including the propto comparisons of the GLM and student_t_lpdf, except opencl_matches_cpu_big of the GLM. That test runs after the exception of defect 4 in the same binary, and the exception leaves the autodiff stack dirty. For items 5 and 6, I swapped in the negative binomial headers of the PR version before the review fixes (468ff26dbb): exactly the 10 new rows of large_shapes_test.cpp fail (13 expectations), and each of the four OpenCL negative binomial suites fails only its reference test, at the new rows.

Full suites. At 468ff26dbb, all unit tests of the touched functions pass (51 files: prim, rev, fwd, mix, including lbeta, digamma_diff, student_t, multi_student_t and lkj_cov), and so do the generated test/prob tests of beta_binomial, beta_neg_binomial, neg_binomial, neg_binomial_2, yule_simon and student_t (271 306 cases). The OpenCL suites of the six distributions, digamma_diff, lbeta and the kernel generator pass on a V100 (126 cases). The review fixes changed only the four negative binomial functions. At the current head, after the merge of develop d12f488749, these pass: the 25 unit test files of the negative binomial family (including the Laplace tests that use neg_binomial_2_log and the new glm_forward_mode_tape_test from develop), the generated test/prob/neg_binomial and test/prob/neg_binomial_2 tests (23 490 cases), and the four OpenCL negative binomial suites (V100). Doxygen 1.9.7 reports no warnings at 01c1a90f87.

Side effects

None

Release notes

Improved stability of beta_binomial_lpmf, beta_binomial_cdf, beta_binomial_lcdf, beta_binomial_lccdf, beta_neg_binomial_lpmf, beta_neg_binomial_cdf, beta_neg_binomial_lcdf, beta_neg_binomial_lccdf, neg_binomial_lpmf, neg_binomial_2_lpmf, neg_binomial_2_log_lpmf, neg_binomial_2_log_glm_lpmf, dirichlet_multinomial_lpmf, lkj_corr_lpdf, lkj_corr_cholesky_lpdf, lkj_cov_lpdf, student_t_lpdf, multi_student_t_lpdf, multi_student_t_cholesky_lpdf, yule_simon_lpmf, yule_simon_cdf, yule_simon_lcdf, yule_simon_lccdf, and the derivatives of lbeta.

Checklist

• Copyright holder: Aalto University

  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

  • unit tests pass (to run, use: ./runTests.py test/unit)
  • header checks pass, (make test-headers)
  • dependencies checks pass, (make test-math-dependencies)
  • 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

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.23	0.22	1.01	1.35% faster
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.06	0.06	0.99	-0.85% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	6.31	6.28	1.0	0.35% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.02	0.02	1.01	0.95% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	8.41	8.41	1.0	0.04% faster
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	4.46	4.45	1.0	0.32% faster
stat_comp_benchmarks/benchmarks/sir/sir.stan	169.67	168.85	1.0	0.49% faster
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.9	0.89	1.01	0.94% faster
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.63	0.62	1.02	1.53% faster
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	43.03	42.74	1.01	0.68% faster
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.59	0.59	1.01	0.63% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	20.88	20.82	1.0	0.29% faster
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.0	0.18% faster
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.19	3.18	1.0	0.21% faster
performance.compilation	389.62	386.82	1.01	0.72% faster
Mean result: 1.0052754019026933

Jenkins Console Log
Jenkins Build Stages
Commit hash: 44a50912ff3f2e71023ea29a6f91adf7c6b1cf87

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           43 bits physical, 48 bits virtual
CPU(s):                                  256
On-line CPU(s) list:                     0-255
Thread(s) per core:                      2
Core(s) per socket:                      64
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              23
Model:                                   49
Model name:                              AMD EPYC 7742 64-Core Processor
Stepping:                                0
Frequency boost:                         enabled
CPU MHz:                                 1497.237
CPU max MHz:                             3416.0681
CPU min MHz:                             1500.0000
BogoMIPS:                                4491.85
Virtualization:                          AMD-V
L1d cache:                               4 MiB
L1i cache:                               4 MiB
L2 cache:                                64 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-63,128-191
NUMA node1 CPU(s):                       64-127,192-255
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Mitigation; untrained return thunk; SMT enabled with STIBP protection
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Retpolines; IBPB conditional; STIBP always-on; RSB filling; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Not affected
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba ibrs ibpb stibp vmmcall fsgsbase bmi1 avx2 smep bmi2 cqm rdt_a rdseed adx smap clflushopt clwb sha_ni xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif v_spec_ctrl umip rdpid overflow_recov succor smca sev sev_es

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin

[jachymb]

@avehtari Hi! I was working on something related and took a similar, but a bit different approach. Here are some notes on how this can be further improved.

I have this on a branch with Hessian tests (η from −800 to 720) and the Laplace test. I can fold it into this PR or rebase it as a follow-up once this is merged, or you can paste these to your agent. lmk what is a preferred option. :)

Comments

Checked at 468ff26 with a probe (var gradients, stan::math::hessian for fvar<var> and fvar<fvar<double>>), the same probe on develop 6110e53, and mpmath as reference. The large-φ rows reproduce. These fail:

quantity	inputs	this branch	develop	exact
GLM ∂/∂φ, (theta_exp - y) / (theta_exp + phi), CPU and kernel	y=3, x=1, α=0, β=720 (data), φ=1.5	NaN	−717.2421539395109	−717.24215393951088
GLM ∂/∂β, y - theta_exp (y+φ) / (theta_exp+φ), CPU and kernel	y=5, x=1, α=0, β=720, φ=1.5	NaN	NaN	−1.5
same	β=709	−inf	−inf	−1.5
GLM ∂²/∂θ² (fvar<var> and fvar<fvar<double>>)	y=1, θ=35, φ=0.07	0	0	−4.7225e−17
same	y=1, θ=27, φ=0.07	−1.4097e−13	same	−1.40777e−13
GLM value, binomial_coefficient_log(y + phi - 1, y)	y=1, θ=0, φ=1e−17	−inf	−39.143946580898778	−39.143946580898777
neg_binomial_2_log_lpmf Hessian	n=3, φ=0.07, η=720 / η=−400	NaN entries	same	−4.4e−314 / −8.4e−173

Rows 1 and 6 are regressions of this branch (inf/inf; y + phi is rounded to ulp(y) before the call). The rest is exp(theta) on the tape, as on develop. The second-derivative error (quotient rule, relative error ≈ eps·e^θ/φ) makes block_matrix_sqrt reject a rank-deficient W in the embedded Laplace approximation: a Laplace test of mine (user likelihood calling the GLM, e^θ/φ ≈ 6e9) throws on this branch and on develop.

All of it goes away by never forming exp(theta): with $u = \theta - \log\phi$ and $e = e^{-|u|}$,

$$s=\frac{e^{\theta}}{e^{\theta}+\phi}=\begin{cases}1/(1+e) & u>0\ e/(1+e) & u\le 0\end{cases},\qquad 1-s=\begin{cases}e/(1+e) & u>0\ 1/(1+e) & u\le 0\end{cases}$$ $$\frac{\partial\log p}{\partial\theta}=y-(y+\phi)\,s,\qquad \frac{\partial\log p}{\partial\phi}=s-\frac{y}{\phi}\,(1-s)-\log\left(1+e^{u}\right)+\psi(y+\phi)-\psi(\phi)$$

where log(1+e^u) is your log1p_exp_theta_m_log_phi and ψ(y+φ)−ψ(φ) your digamma_diff(phi, y), both from the same e (log1p_exp(±u) = max(±u, 0) + log1p(e)); same in the kernel and in neg_binomial_2_log_lpmf. For the value, -lbeta(phi, y + 1) - log(y + phi) is what binomial_coefficient_log evaluates for y + φ ≥ 10, with φ unrounded.

[avehtari]

@jachymb thanks for the review and suggestions. I'll check these

[avehtari]

@jachymb Rows 1 and 6 are regressions of this PR, and I could reproduce them. Your formulas are correct, and I used them. In the other PR checking your comments lead to additional findings, so I did ask Claude to check your comments, too, but I guess here the fix is the same as you had. Below there is one request for you

• neg_binomial_2_log_glm_lpmf (CPU and OpenCL kernel): the partials use 1 / (1 + e) and e / (1 + e) with e = exp(-|theta - log(phi)|) and do not form exp(theta). This fixes rows 1 to 3. It also fixes the GLM second derivatives of rows 4 and 5, because the derivative of e / (1 + e) does not cancel: the fvar<var> Hessian entry is now -4.7225e-17 at theta = 35 (before: 0) and -1.40777e-13 at theta = 27 (before: -1.40970e-13).
• The value term lchoose(y + phi - 1, y) is now -lbeta(phi, y + 1) - log(y + phi). The defect is wider than row 6: binomial_coefficient_log(y + phi - 1, y) rounds phi into y + phi - 1, so the value has an error of about eps * y / phi. Develop has this form in neg_binomial_2_lpmf, neg_binomial_2_log_lpmf and neg_binomial_lpmf (with alpha), and this PR had brought it into the GLM. At y = 1000, phi = 1e-8 the error was 7.9e-7 (5.0e-9 in neg_binomial_lpmf), and at phi = 1e-17 the value was -inf. The four functions, CPU and OpenCL, now use the lbeta form; at these points the errors are at most 5.3e-13.
• New test rows with mpmath references in large_shapes_test.cpp and in the four OpenCL reference tests: the GLM at theta = 720 and 709, and the four functions at phi = 1e-8 and 1e-17. With the previous head of this PR exactly these rows fail. With the new commit all unit tests of the negative binomial family, the Laplace tests that use neg_binomial_2_log, the generated test/prob tests (23 490 cases) and the OpenCL suites (V100) pass.

Row 7, the neg_binomial_2_log_lpmf Hessian, has the same quotient-rule cancellation, and develop has it too. To keep this PR reviewable, could you (@jachymb) make it a follow-up PR after this one is merged, together with your Hessian tests and the Laplace test? Your e form applies there in the same way.

I also saw #3446 for the inv_logit derivatives. The GLM uses the explicit e / (1 + e) form, so it does not depend on #3446.

[WardBrian]

@avehtari could you resolve the conflicts in this PR to let jenkins run?

[avehtari]

conflict resolved (it was only in the include block)

[stan-buildbot]

Name	Old Result	New Result	Ratio	Performance change( 1 - new / old )
stat_comp_benchmarks/benchmarks/gp_regr/gen_gp_data.stan	0.04	0.04	0.99	-1.21% slower
stat_comp_benchmarks/benchmarks/gp_regr/gp_regr.stan	0.16	0.17	0.96	-4.08% slower
stat_comp_benchmarks/benchmarks/garch/garch.stan	0.87	0.86	1.01	1.24% faster
stat_comp_benchmarks/benchmarks/low_dim_corr_gauss/low_dim_corr_gauss.stan	0.01	0.01	0.94	-6.19% slower
stat_comp_benchmarks/benchmarks/eight_schools/eight_schools.stan	0.11	0.11	1.02	2.42% faster
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix/low_dim_gauss_mix.stan	7.97	7.98	1.0	-0.11% slower
stat_comp_benchmarks/benchmarks/gp_pois_regr/gp_pois_regr.stan	3.69	3.75	0.98	-1.71% slower
stat_comp_benchmarks/benchmarks/pkpd/one_comp_mm_elim_abs.stan	32.92	32.96	1.0	-0.1% slower
stat_comp_benchmarks/benchmarks/pkpd/sim_one_comp_mm_elim_abs.stan	0.38	0.38	1.01	0.59% faster
stat_comp_benchmarks/benchmarks/irt_2pl/irt_2pl.stan	6.79	6.83	0.99	-0.61% slower
stat_comp_benchmarks/benchmarks/sir/sir.stan	136.74	136.78	1.0	-0.03% slower
stat_comp_benchmarks/benchmarks/arK/arK.stan	3.21	3.23	0.99	-0.77% slower
stat_comp_benchmarks/benchmarks/arma/arma.stan	0.6	0.61	0.99	-0.78% slower
stat_comp_benchmarks/benchmarks/low_dim_gauss_mix_collapse/low_dim_gauss_mix_collapse.stan	25.57	25.26	1.01	1.23% faster
performance.compilation	387.23	379.43	1.02	2.01% faster
Mean result: 0.9950799155442522

Jenkins Console Log
Jenkins Build Stages
Commit hash: 29e32c4215464fd5550540ce9a448016d22be56f

Machine information

Distributor ID:	Ubuntu
Description:	Ubuntu 20.04.3 LTS
Release:	20.04
Codename:	focal

CPU:

Architecture:                            x86_64
CPU op-mode(s):                          32-bit, 64-bit
Byte Order:                              Little Endian
Address sizes:                           52 bits physical, 57 bits virtual
CPU(s):                                  192
On-line CPU(s) list:                     0-191
Thread(s) per core:                      2
Core(s) per socket:                      48
Socket(s):                               2
NUMA node(s):                            2
Vendor ID:                               AuthenticAMD
CPU family:                              25
Model:                                   17
Model name:                              AMD EPYC 9474F 48-Core Processor
Stepping:                                1
Frequency boost:                         enabled
CPU MHz:                                 1497.526
CPU max MHz:                             4114.4229
CPU min MHz:                             1500.0000
BogoMIPS:                                7189.39
Virtualization:                          AMD-V
L1d cache:                               3 MiB
L1i cache:                               3 MiB
L2 cache:                                96 MiB
L3 cache:                                512 MiB
NUMA node0 CPU(s):                       0-47,96-143
NUMA node1 CPU(s):                       48-95,144-191
Vulnerability Gather data sampling:      Not affected
Vulnerability Indirect target selection: Not affected
Vulnerability Itlb multihit:             Not affected
Vulnerability L1tf:                      Not affected
Vulnerability Mds:                       Not affected
Vulnerability Meltdown:                  Not affected
Vulnerability Mmio stale data:           Not affected
Vulnerability Old microcode:             Not affected
Vulnerability Reg file data sampling:    Not affected
Vulnerability Retbleed:                  Not affected
Vulnerability Spec rstack overflow:      Mitigation; Safe RET
Vulnerability Spec store bypass:         Mitigation; Speculative Store Bypass disabled via prctl
Vulnerability Spectre v1:                Mitigation; usercopy/swapgs barriers and __user pointer sanitization
Vulnerability Spectre v2:                Mitigation; Enhanced / Automatic IBRS; IBPB conditional; STIBP always-on; PBRSB-eIBRS Not affected; BHI Not affected
Vulnerability Srbds:                     Not affected
Vulnerability Tsa:                       Mitigation; Clear CPU buffers
Vulnerability Tsx async abort:           Not affected
Vulnerability Vmscape:                   Mitigation; IBPB before exit to userspace
Flags:                                   fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush mmx fxsr sse sse2 ht syscall nx mmxext fxsr_opt pdpe1gb rdtscp lm constant_tsc rep_good amd_lbr_v2 nopl xtopology nonstop_tsc cpuid extd_apicid aperfmperf rapl pni pclmulqdq monitor ssse3 fma cx16 pcid sse4_1 sse4_2 x2apic movbe popcnt aes xsave avx f16c rdrand lahf_lm cmp_legacy svm extapic cr8_legacy abm sse4a misalignsse 3dnowprefetch osvw ibs skinit wdt tce topoext perfctr_core perfctr_nb bpext perfctr_llc mwaitx cpb cat_l3 cdp_l3 hw_pstate ssbd mba perfmon_v2 ibrs ibpb stibp ibrs_enhanced vmmcall fsgsbase bmi1 avx2 smep bmi2 erms invpcid cqm rdt_a avx512f avx512dq rdseed adx smap avx512ifma clflushopt clwb avx512cd sha_ni avx512bw avx512vl xsaveopt xsavec xgetbv1 xsaves cqm_llc cqm_occup_llc cqm_mbm_total cqm_mbm_local user_shstk avx512_bf16 clzero irperf xsaveerptr rdpru wbnoinvd amd_ppin cppc arat npt lbrv svm_lock nrip_save tsc_scale vmcb_clean flushbyasid decodeassists pausefilter pfthreshold avic v_vmsave_vmload vgif x2avic v_spec_ctrl vnmi avx512vbmi umip pku ospke avx512_vbmi2 gfni vaes vpclmulqdq avx512_vnni avx512_bitalg avx512_vpopcntdq la57 rdpid overflow_recov succor smca fsrm flush_l1d debug_swap

G++:

g++ (Ubuntu 9.4.0-1ubuntu1~20.04) 9.4.0
Copyright (C) 2019 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

Clang:

clang version 10.0.0-4ubuntu1 
Target: x86_64-pc-linux-gnu
Thread model: posix
InstalledDir: /usr/bin