Use Eigen's igamma_der_a for the incomplete gamma gradient roots

[avehtari]

Use Eigen's igamma_der_a for the incomplete gamma gradient roots (this PR was Claude assisted)

Summary

grad_reg_inc_gamma and grad_reg_lower_inc_gamma compute the derivative of the regularized incomplete gamma function with respect to its shape parameter. They are the gradient roots under 27 tail functions in 10 distributions. Both have defects that make the gradient wrong by tens of percent, or NaN, in regions that ordinary models reach. This PR replaces both implementations with Eigen::numext::igamma_der_a, which is already vendored in lib/eigen_5.0.1. The result is more accurate everywhere measured, about four times faster, and about 120 lines shorter.

The values of the distributions do not change. Only the gradients change.

Comparison result, one table

	develop	branch
non-finite, 49 990 points	6 379	0
median rel. error, worst region	5.2e-03	6.5e-14
worst rel. error outside a≈z≈10⁴	inf	3e-11
worst rel. error at a≈z≈10⁴	NaN	5.2e-06, documented
autodiff-order spread	inf in 7/10 regions	2.2e-14
forward mode, z ≳ 700	does not terminate	finite
speed, double roots	1x	2.2x–4.4x
speed, gamma_lcdf per obs, α=50	1x	9.8x
speed, var root (unused in-repo)	1x	0.72x, explained

The defects in the current code on develop

Everything in this section and the next describes the code on develop at 5252d51d47, before this PR. All numbers were measured against an mpmath reference at 80 digits, built through two independent routes that agree to at least 13 digits, and confirmed by compiling and running the real develop headers rather than a transcription. The section "Accuracy after the change" below gives the same measurements for this PR.

1. On develop, grad_reg_inc_gamma takes an invalid branch when z is near a. The test z >= a && z >= 8 selects a fixed 10-term asymptotic expansion from DLMF 8.11.2. That expansion needs z >> a. The whole band a <= z < 2a takes it, and there the terms do not decay. The relative error at z = a is 4.4e-02 at a = 20, 5.8e-01 at a = 100 and 7.2e-01 at a = 171. The expansion also ignores the precision argument and always takes exactly 10 terms.

2. On develop, grad_reg_inc_gamma returns NaN for a above about 171.6. Every call site computes g = tgamma(a) and passes it in. Above 171.62 that is inf, so the correction term exp(-z + (a-1) log z) * S / g is inf / inf. There is no guard at any call site. The series branch fails the same way through exp(a * log z), which overflows when a log z > 709, for example at (a, z) = (150, 150).

3. On develop, grad_reg_lower_inc_gamma loses all accuracy for z well above a. Its delegation test is (a < 0.8 && z > 15) || (a < 12 && z > 30) || a < sqrt(-756 - z*z + 60*z). The third term needs -z^2 + 60z - 756 > 0, which holds only for z in (18, 42). So for z > 42 and a >= 12 the Gautschi branch always runs, and there emz * (log_z * sum_a - sum_b) is a difference of two large nearly equal sums. It returns about 1e-14 of rounding noise, with arbitrary sign, where the true value can be 1e-148. Examples: at (50, 150) the reference is -8.29e-22 and develop returns -4.70e-15; at (50, 500) the reference is -5.34e-148 and develop returns +1.47e-14.

4. On develop, the 1e-6 default precision is real but small. It sets the error of the plain series branch, 9.2e-08 at a = 0.5 and 9.0e-05 at a = 8. It is the smallest of the four problems.

On develop, the defects reach the user

Measured at the distribution level, with var, against the same reference. The CDF value is correct in every case; only the gradient is wrong. This PR brings every row below to the reference within 1e-13; see "Accuracy after the change".

function	point	reference	develop	rel err
chi_square_lccdf	nu = 200	0.0410184911365	0.0173344546475	5.8e-01
chi_square_lccdf	nu = 342	0.0311564510080	0.0086722587902	7.2e-01
chi_square_lcdf	nu = 200	-0.0388930365624	-0.0164362354566	5.8e-01
gamma_cdf	alpha = 100	-0.0399274978578	-0.0168733998162	5.8e-01
inv_gamma_lccdf	alpha = 100	-0.0777860731247	-0.0328724709132	5.8e-01

NaN gradients, confirmed at the distribution level: chi_square_lccdf and chi_square_lcdf at nu = 360, gamma_cdf at alpha = 180, inv_gamma_lccdf at alpha = 180, scaled_inv_chi_square_lcdf at nu = 360.

gamma_lcdf and gamma_lccdf are not affected by defects 1 and 2. Their log_gamma_q_dgamma dispatch avoids that path, and they measure 1e-14 to 5e-13 in the same region. gamma_lcdf does inherit defect 3.

What this PR does

Eigen::internal::igamma_generic_impl differentiates the Cephes power series and the Cephes continued fraction term by term, stops at machine epsilon, and branches on x > 1 && x > a, which is the well-conditioned split. It is already in the tree, it is tested upstream, and it is generic in the scalar type.

Only three of its helpers are restricted to float and double: lgamma_impl, digamma_impl and cephes_helper. This PR supplies those three for Stan's autodiff scalars in two new headers, stan/math/fwd/fun/Eigen_SpecialFunctions.hpp and stan/math/rev/fun/Eigen_SpecialFunctions.hpp. igamma_num_iterations already falls back correctly and needs nothing.

Both roots then become thin wrappers, and the hand-written series, the asymptotic branch and the Gautschi branch are deleted. One algorithm now serves every scalar type and every autodiff order, so reverse mode, forward mode, fvar<fvar<double>>, fvar<var> and hessian() all agree.

Doing only the double path was tried first and rejected. It left hessian() on the old code, so gradient() and hessian() returned first derivatives that differed by 5.8e-01 at alpha = 100 and by NaN at alpha = 180. Making every order use one algorithm is the point of the design.

The signatures do not change, so no call site changes. g, dig, precision and max_steps are accepted and ignored. Removing them touches 27 call sites and is proposed separately, to keep this PR reviewable.

Accuracy after the change, this PR

Sweep of 1800 points in nine regions chosen by the branch structure, scored against the mpmath reference. Worst relative error over every region, with the count of non-finite results out of 200 per region.

region	develop grad_reg_inc_gamma	develop grad_reg_lower_inc_gamma	this PR
series, a < 1	1.5e-03	1.5e-10	8.5e-15
series, 10 < a <= 170	2.6e+29	5.8e-01	1.9e-13
asymptotic, a <= z <= 2a	5.1e-01	1.1e+07	1.7e-13
asymptotic, z >= 3a	3.7e-04	2.1e+165	3.8e-12
integer a	1.9e+15	2.6e+01	3.3e-14
a > 171.6	176 of 200 NaN	1.1e+01	3.8e-13
z >= 3a, a >= 12	9.6e-07	1.9e+127	4.3e-12

Worst over all 1800 points after the change: 4.3e-12. No non-finite result anywhere, including the 200 points above the tgamma overflow limit.

A second, larger sweep compared develop and this branch head to head: 49 990 points with a and z both in [1e-3, 1e4], with dense sampling on z ≈ a, the z = 8 seam, a = 171.62 and a·log z = 709. Develop returns NaN or inf at 6 379 of those points for grad_reg_inc_gamma. This branch is finite at every one. Median relative error on this branch is between 1e-16 and 1e-14 in every region; develop's medians reach 5e-3.

One known limit of the new implementation. Near z = a the error grows with a: worst 3e-14 for a < 10, 3e-11 for a < 1 000, 2.9e-7 for a < 10 000, and 5.2e-6 at the a = z = 10⁴ corner, where Eigen's Cephes recurrence loses digits near the mode. At that point three independent references agree to 17 digits, so this is the algorithm, not the measurement. It is finite where develop is NaN, and below the 1e-4 that expect_ad tolerates.

A reference caveat. In the far upper tail at a ≳ 10⁴ the mpmath reference is certified to 11 digits only, so errors reported there are bounded at 1e-11 by the reference, not by the code.

Autodiff orders agree on this branch to 2.2e-14 across double, var, fvar<double> and fvar<fvar<double>>. On develop, one order returns NaN where another returns a number in 7 of 10 regions.

Speed

Develop and this branch compiled from one source and run back to back on the same exclusive node: Xeon E5-2680 v3 at 2.50 GHz, performance governor, 1024-point input block, with an accumulator-dependent term added to every input so the inner loop is not loop-invariant, a baseline subtracted, and every result checked against a plausible cycle count. Two rounds each; run-to-run spread under 3 %. The input grid stays where develop terminates and is finite, a in [0.5, 60] and z in [0.1, 40], which favours develop.

Function level, net ns per call:

kernel	develop	this PR	ratio
grad_reg_inc_gamma, double, incl. tgamma+digamma	805	370	2.2x
grad_reg_inc_gamma, fvar	1 749	592	3.0x
grad_reg_inc_gamma, fvar<fvar>	3 230	1 320	2.4x
grad_reg_lower_inc_gamma, double	1 135	256	4.4x
grad_reg_lower_inc_gamma, fvar	1 766	468	3.8x
grad_reg_lower_inc_gamma, fvar<fvar>	2 536	1 133	2.2x
grad_reg_inc_gamma, var, value + grad	2 690	3 727	0.72x
grad_reg_lower_inc_gamma, var, value + grad	3 921	3 618	1.1x

The one slower row is the var instantiation of the upper root. Eigen's series and continued fraction take more operations than develop's 10-term asymptotic branch on the part of this grid where that branch runs, and with var every operation is a tape node. No function in the library calls the roots with var; every distribution passes double partials. The row is included for completeness.

Distribution level, which is what a model pays: N = 1000 observations, one shared var shape, value plus reverse sweep timed together, ns per observation:

function	alpha	develop	this PR	ratio
gamma_cdf	5	561	353	1.6x
gamma_lcdf	5	1 331	366	3.6x
gamma_lccdf	5	862	393	2.2x
chi_square_lccdf	5	1 072	845	1.3x
inv_gamma_lcdf	5	1 087	905	1.2x
gamma_cdf	50	1 546	497	3.1x
gamma_lcdf	50	4 974	509	9.8x
gamma_lccdf	50	1 449	544	2.7x
chi_square_lccdf	50	1 934	941	2.1x
inv_gamma_lcdf	50	840	755	1.1x

Forward mode did not terminate on develop

stan/math/fwd/fun/gamma_q.hpp on develop held its own inlined copy of the series with while (fabs(delta) > 1e-6) and no iteration cap. For z above about 700 the alternating terms (-z)^k / k! overflow to inf before they decay, delta stays inf, and the loop never exits. The root's own asymptotic branch calls gamma_q(a, z) with fvar arguments, so grad_reg_inc_gamma reached it too. In the comparison runs a single call at a = 0.168, z = 4 748 did not return in 15 minutes. The existing mix test for gamma_q carries a commented-out expect_value(f, 8.01006, 2.47579e+215) with the note "this still fails forward mode", which is the same defect.

This PR removes that code, so the loop is gone. Because a test that hangs is not usable in CI, the regression test mathMixScalFun.gammaQ_fwd_shape_derivative_large_z uses points where the old series terminates but is wrong: at (a, z) = (1.5, 300) develop returns 1.2e+112 for a derivative whose value is 5.7e-129. On develop the test fails; on this PR it passes.

Testing

New file test/unit/math/prim/fun/grad_reg_inc_gamma_accuracy_test.cpp, with fixed references from mpmath at 80 digits. expect_ad compares against finite differences and cannot see errors of this size, so these are absolute references rather than autodiff comparisons.

New case ProbInternalMath.gradRegIncGamma_gradient_matches_hessian asserts that the first derivative of gamma_cdf is the same whether it comes from reverse mode or from hessian(). It needs no external reference.

Two new cases in test/unit/math/mix/fun/gamma_q_test.cpp: gammaQ_fwd_shape_derivative_large_z (the indirect non-termination test described above) and gammaQ_fwd_z_derivative_large_a, which checks the second-argument derivative at (a, z) = (200, 150) against the closed form -z^(a-1) e^(-z) / Γ(a); on develop that returned NaN because both pow and tgamma overflow. Installing the develop versions of the three changed files makes both cases fail and the rest of the file pass; restoring this PR makes all four cases pass, with the file checksums back to their starting values.

Evidence that the tests find the defects, by installing the develop headers, keeping the new tests, and restoring afterwards:

stage	suites failed	cases failed
this PR	0 of 4	0
develop headers	2 of 4	10
this PR again	0 of 4	0

The header checksums in the third stage equal the first, so the swap is clean.

Three constants in test/unit/math/mix/fun/grad_reg_inc_gamma_test.cpp had to be corrected. They expected 0.38984156 for d/da Q(0.5, 1.0), whose true value is 0.38983726432851057. The old constant encoded two errors together: the rounded support values that the test passes in, since g = 1.77245 is 2.2e-06 away from tgamma(0.5), and the 1e-6 series truncation.

Two cases in that file were not touched and still pass: gradRegIncGamma_fv_1stderiv and gradRegIncGamma_fv_2ndderiv assert the second and third derivatives at (0.5, 1.0) to 1e-6 against constants that predate this work. That is independent confirmation of the higher-order derivatives.

The existing grad_reg_lower_inc_gamma_test checks 9600 fixed values from Mathematica with a 1e-10 absolute tolerance. Eigen satisfies all 9600. That grid covers a in [0.0001, 19.75] and z in [0.0001, 29.75], and contains no point with a >= 12 and z >= 3a, which is exactly the region where defect 3 lives. The grid stopped one step short of it.

On develop, a duplicated copy of the same series in fwd/fun/gamma_q.hpp

On develop, stan/math/fwd/fun/gamma_q.hpp carried its own inlined copy of the series in grad_reg_inc_gamma: the same while (fabs(delta) > 1e-6), the same (1.0 - u) * (dig - l) + exp(x1 * l) * S / g, the same hard-coded tolerance, and without the second branch. This PR replaces both copies with a call to grad_reg_inc_gamma, which is the function they were copied from. stan/math/fwd/fun/gamma_p.hpp was already written this way, and calls grad_reg_lower_inc_gamma.

The duplicate had been hiding the defect. In test/prob/chi_square, the generated third-order case compares the distribution's analytic partials against autodiff through the reference formula log(gamma_q(nu * 0.5, y * 0.5)). Both routes used the same inaccurate series, so their errors cancelled and the comparison passed. With an accurate root and a stale copy, the same comparison fails by 2.4e-03. After the copy is removed, it passes with the new root and with the old one, because the two routes are then consistent either way.

One item in that file is left for a follow-up: der2 is -exp(-x2) * pow(x2, x1 - 1) / tgamma(x1), which overflows for x1 above 171.6 and silently returns -0, in the same way as defect 2 above. gamma_p already uses the log-space form.

Dependency on Eigen

The two new headers specialize lgamma_impl, digamma_impl and cephes_helper inside Eigen::internal. Stan Math already specializes inside that namespace, in stan/math/rev/core/Eigen_NumTraits.hpp. The cost is that an Eigen upgrade which renames those three helpers breaks the build. That failure is loud and immediate, not silent, and the alternative was to copy about 300 lines of the Cephes algorithm into Stan Math and maintain it here.

The headers belong to the fun layer, not core, because they need digamma and lgamma. Putting them beside Eigen_NumTraits.hpp in core breaks the include order of every mix target.

Not included

OpenCL is unaffected: no OpenCL code calls these roots, and the five affected families have OpenCL implementations only for their _lpdf, which does not use the incomplete gamma.

The incomplete beta gradient roots, inc_beta_dda, inc_beta_ddb and grad_reg_inc_beta, have their own defects: they stop on an absolute threshold of 1e-10 and reach 12 % relative error at (a, b, z) = (2500, 2, 0.999), a case their own docstring names as tested. Boost supplies only ibeta_derivative, which is the derivative with respect to x, so there is no vendored replacement. The published algorithm is Boik and Robinson-Cox (1998), Journal of Statistical Software 3(1). That work is separate from this PR.

Removing the now-unused g, dig, precision and max_steps arguments is also a separate PR.

Release notes

Use Eigen's igamma_der_a for the incomplete gamma gradient roots

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

  • new 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

[avehtari]

I'll be away for two weeks, but did create this PR to avoid overlap in work. Similar PR should be made for incomplete beta gradient roots (mentioned in PR text)