integrate_1d_gauss_kronrod performance fix

.gitattributes

@@ -5,3 +5,7 @@
 *.tex text diff=tex
 lib/** binary
 lib/eigen_5.0.1/STAN_CHANGES.md -binary text diff=markdown
+lib/boost_1.87.0/STAN_CHANGES -binary text
+# Stan-patched vendored source: keep it diffable so the patch stays reviewable
+# and so a future Boost upgrade that drops it is visible rather than silent.
+lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp -binary text diff=cpp

lib/boost_1.87.0/STAN_CHANGES

@@ -0,0 +1,9 @@
+This file documents changes done for the stan-math project
+
+- Added an `abs_tol` parameter to
+  `boost/math/quadrature/gauss_kronrod.hpp`'s public
+  `gauss_kronrod<Real, N>::integrate()`, and forwarded it to the five
+  `recursive_adaptive_integrate()` entry points in place of the hard-coded
+  `Real(0)`.
+
+  Used by `stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp`.

lib/boost_1.87.0/boost/math/quadrature/gauss_kronrod.hpp

@@ -1227,7 +1227,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
 
 public:
    template <class F>
-   static auto integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon<Real>(), Real* error = nullptr, Real* pL1 = nullptr)->decltype(std::declval<F>()(std::declval<Real>()))
+   static auto integrate(F f, Real a, Real b, unsigned max_depth = 15, Real tol = tools::root_epsilon<Real>(), Real* error = nullptr, Real* pL1 = nullptr, Real abs_tol = Real(0))->decltype(std::declval<F>()(std::declval<Real>()))
    {
       typedef decltype(f(a)) K;
       static_assert(!std::is_integral<K>::value,
@@ -1248,7 +1248,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
                return res;
             };
             recursive_info<decltype(u)> info = { u, tol };
-            K res = recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+            K res = recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, abs_tol, error, pL1);
             return res;
          }
 
@@ -1263,7 +1263,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
                return res;
             };
             recursive_info<decltype(u)> info = { u, tol };
-            K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+            K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, abs_tol, error, pL1);
             if (pL1)
             {
                *pL1 *= 2;
@@ -1280,7 +1280,7 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
                return f(b - arg) * z * z;
             };
             recursive_info<decltype(v)> info = { v, tol };
-            K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, Real(0), error, pL1);
+            K Q = Real(2) * recursive_adaptive_integrate(&info, Real(-1), Real(1), max_depth, abs_tol, error, pL1);
             if (pL1)
             {
                *pL1 *= 2;
@@ -1297,9 +1297,9 @@ class gauss_kronrod : public detail::gauss_kronrod_detail<Real, N, detail::gauss
             recursive_info<F> info = { f, tol };
             if (b < a)
             {
-               return -recursive_adaptive_integrate(&info, b, a, max_depth, Real(0), error, pL1);
+               return -recursive_adaptive_integrate(&info, b, a, max_depth, abs_tol, error, pL1);
             }
-            return recursive_adaptive_integrate(&info, a, b, max_depth, Real(0), error, pL1);
+            return recursive_adaptive_integrate(&info, a, b, max_depth, abs_tol, error, pL1);
          }
       }
       return static_cast<K>(policies::raise_domain_error(function, "The domain of integration is not sensible; please check the bounds.", a, Policy()));

stan/math/fwd/functor/integrate_1d_gauss_kronrod.hpp

@@ -12,7 +12,7 @@ namespace stan {
 namespace math {
 
 /**
- * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21)
+ * Return the integral of f from a to b using adaptive Gauss-Kronrod (G10,K21)
  * quadrature, with tangents computed via finite differences over the
  * integrand parameters.
  *
@@ -64,7 +64,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_gauss_kronrod_tol(
 }
 
 /**
- * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21)
+ * Return the integral of f from a to b using adaptive Gauss-Kronrod (G10,K21)
  * quadrature, with tangents computed via finite differences over the
  * integrand parameters.
  *

stan/math/prim/functor/integrate_1d_gauss_kronrod.hpp

@@ -29,7 +29,7 @@ constexpr int INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH = 15;
 
 /**
  * Integrate a single variable function f from a to b using Boost's adaptive
- * Gauss-Kronrod (G21,K21) quadrature, with QUADPACK-style mixed convergence
+ * Gauss-Kronrod (G10,K21) quadrature, with QUADPACK-style mixed convergence
  * criterion. The integration succeeds (returns the Boost estimate Q)
  * whenever
  *   error <= max(relative_tolerance * L1, absolute_tolerance)
@@ -41,9 +41,18 @@ constexpr int INTEGRATE_1D_GAUSS_KRONROD_MAX_DEPTH = 15;
  * checking accumulated floating-point round-off against itself (this
  * happens routinely in nested integrate_1d_gauss_kronrod calls when the
  * outer integration probes the deep tail of the integrand and every
- * inner evaluation sees an essentially-zero integrand). Setting it to
- * zero (the default) reproduces the strict pure-relative-tolerance
- * behaviour of integrate_1d.
+ * inner evaluation sees an essentially-zero integrand).
+ *
+ * absolute_tolerance is applied twice, in the same units: as a floor on
+ * refinement inside Boost's adaptive recursion (a panel whose error already
+ * sits below the floor is not bisected, which is what bounds the work in the
+ * round-off regime above) and as the floor on the convergence test below.
+ *
+ * Setting it to zero (the default) reproduces the strict
+ * pure-relative-tolerance behaviour of integrate_1d. Note that zero is
+ * Boost's sentinel for "derive the refinement budget from the root panel's
+ * own relative target", so a positive but negligible absolute_tolerance
+ * removes that derived budget and can refine slightly MORE than zero does.
  *
  * The signature for f should be:
  *   double f(double x, double xc)
@@ -86,7 +95,7 @@ inline double integrate_gk(const F& f, double a, double b,
   const unsigned int depth
       = max_depth < 0 ? 0u : static_cast<unsigned int>(max_depth);
   double Q = gauss_kronrod<double, INTEGRATE_1D_GAUSS_KRONROD_ORDER>::integrate(
-      f_wrap, a, b, depth, relative_tolerance, &error, &L1);
+      f_wrap, a, b, depth, relative_tolerance, &error, &L1, absolute_tolerance);
 
   // QUADPACK-style mixed convergence: throw only if the Boost error
   // exceeds both the relative-tolerance target (rel_tol * L1) and the
@@ -106,7 +115,7 @@ inline double integrate_gk(const F& f, double a, double b,
 
 /**
  * Compute the integral of the single variable function f from a to b to within
- * a specified relative tolerance using adaptive Gauss-Kronrod (G21,K21)
+ * a specified relative tolerance using adaptive Gauss-Kronrod (G10,K21)
  * quadrature. a and b can be finite or infinite.
  *
  * @tparam F type of function to integrate
@@ -147,7 +156,7 @@ inline double integrate_1d_gauss_kronrod_tol(const F& f, double a, double b,
 
 /**
  * Compute the integral of the single variable function f from a to b using
- * adaptive Gauss-Kronrod (G21,K21) quadrature. a and b can be finite or
+ * adaptive Gauss-Kronrod (G10,K21) quadrature. a and b can be finite or
  * infinite.
  *
  * The signature for f should be:

stan/math/rev/functor/integrate_1d_gauss_kronrod.hpp

@@ -16,7 +16,7 @@ namespace stan {
 namespace math {
 
 /**
- * Return the integral of f from a to b using adaptive Gauss-Kronrod (G21,K21)
+ * Return the integral of f from a to b using adaptive Gauss-Kronrod (G10,K21)
  * quadrature.
  *
  * @tparam F Type of f
@@ -57,7 +57,7 @@ inline return_type_t<T_a, T_b, Args...> integrate_1d_gauss_kronrod_tol(
 
 /**
  * Compute the integral of the single variable function f from a to b using
- * adaptive Gauss-Kronrod (G21,K21) quadrature. a and b can be finite or
+ * adaptive Gauss-Kronrod (G10,K21) quadrature. a and b can be finite or
  * infinite.
  *
  * f should be compatible with reverse mode autodiff and have the signature:

test/unit/math/mix/functor/integrate_1d_gauss_kronrod_test.cpp

@@ -25,3 +25,20 @@ TEST(mixFunctor, integrate1DGaussKronrod) {
   // a correctness issue in the integrate_1d_gauss_kronrod wrapper, and so
   // the NaN-input case is intentionally omitted here.
 }
+
+// A positive absolute_tolerance must not disturb higher-order autodiff. The
+// integrand is scaled to 1e-6 so the refinement floor is actually reached,
+// which is the regime where abs_tol changes how many panels are evaluated.
+TEST(mixFunctor, integrate1DGaussKronrodPositiveAbsoluteTolerance) {
+  auto f = [](const auto& theta, const auto& lb, const auto& ub) {
+    auto func = [](const auto& x, const auto& xc, std::ostream* msgs,
+                   const auto& theta) {
+      return 1e-6 * theta * stan::math::sin(7.0 * x);
+    };
+    std::ostringstream* msgs = nullptr;
+    return stan::math::integrate_1d_gauss_kronrod_tol(func, lb, ub, 1e-8, 1e-10,
+                                                      5, msgs, theta);
+  };
+
+  stan::test::expect_ad(f, 0.75, 0.0, 1.0);
+}

test/unit/math/prim/functor/integrate_1d_gauss_kronrod_test.cpp

@@ -444,3 +444,161 @@ TEST(StanMath_integrate_1d_gk_prim, abs_tol_argument_smoke) {
   EXPECT_NEAR(Q0, 1.0423499493102901, 1e-8);
   EXPECT_NEAR(Q1, Q0, 1e-12);
 }
+
+// ---------------------------------------------------------------------------
+// absolute_tolerance during refinement
+//
+// absolute_tolerance is applied in two places, in the same units: as a floor
+// on refinement inside Boost's adaptive recursion, and as the floor on the
+// convergence test. The first of those depends on the abs_tol parameter added
+// to Boost's gauss_kronrod::integrate by a Stan-local patch
+// (lib/boost_1.87.0/STAN_CHANGES).
+//
+// ---------------------------------------------------------------------------
+
+// The motivating case: an integrand so small that the relative-tolerance test
+// degenerates into comparing accumulated round-off against itself. A positive
+// abs_tol stops the pointless refinement.
+//
+// ALSO the regression guard on the vendored-Boost patch -- see above.
+TEST(StanMath_integrate_1d_gk_prim,
+     positive_abs_tol_reduces_work_on_negligible_integrand) {
+  constexpr double scale = 1e-12;
+  constexpr double frequency = 127.0;
+  constexpr double absolute_tolerance = 1e-14;
+  const double expected = scale * (1.0 - std::cos(frequency)) / frequency;
+
+  auto run = [](double abs_tol, int &evaluations) {
+    auto integrand = [&evaluations](double x, double xc, std::ostream *msgs) {
+      ++evaluations;
+      return scale * std::sin(frequency * x);
+    };
+    return stan::math::integrate_1d_gauss_kronrod_tol(
+        integrand, 0.0, 1.0, 1e-12, abs_tol, 5, integrate_1d_gk_test::msgs);
+  };
+
+  int relative_evaluations = 0, absolute_evaluations = 0;
+  const double relative_result = run(0.0, relative_evaluations);
+  const double absolute_result = run(absolute_tolerance, absolute_evaluations);
+
+  // The contract: within the absolute tolerance requested. Asserting anything
+  // tighter would assert an accident of how much more accurate K21 happens to
+  // be than the caller asked for.
+  EXPECT_NEAR(absolute_result, expected, absolute_tolerance);
+  EXPECT_NEAR(relative_result, expected, absolute_tolerance);
+  EXPECT_NEAR(absolute_result, relative_result, absolute_tolerance);
+  EXPECT_LT(std::abs(absolute_result - expected) / std::abs(expected), 1e-10);
+
+  // Materially less work. Stated as a ratio so it does not pin the exact
+  // recursion counts of the quadrature.
+  EXPECT_LT(absolute_evaluations, relative_evaluations);
+  EXPECT_LT(2 * absolute_evaluations, relative_evaluations);
+}
+
+// Regression guard for unbounded refinement. x^{-0.9} has an endpoint
+// singularity Gauss-Kronrod cannot resolve, so no attainable tolerance is met
+// and every panel looks "not yet good enough". Refinement stays panel-local,
+// so the cost stays at the level of the abs_tol == 0 call.
+TEST(StanMath_integrate_1d_gk_prim,
+     positive_abs_tol_bounds_work_on_unresolvable_integrand) {
+  auto run = [](double absolute_tolerance, int &evaluations) {
+    auto integrand = [&evaluations](double x, double xc, std::ostream *msgs) {
+      ++evaluations;
+      return std::pow(x, -0.9);
+    };
+    EXPECT_THROW(stan::math::integrate_1d_gauss_kronrod_tol(
+                     integrand, 1e-300, 1.0, 1e-12, absolute_tolerance, 15,
+                     integrate_1d_gk_test::msgs),
+                 std::domain_error);
+  };
+
+  int relative_evaluations = 0, absolute_evaluations = 0;
+  run(0.0, relative_evaluations);
+  run(1e-14, absolute_evaluations);
+
+  EXPECT_LE(absolute_evaluations, relative_evaluations);
+  EXPECT_LT(absolute_evaluations, 50000);
+}
+
+// Among POSITIVE absolute tolerances, raising the floor never increases work.
+// Zero is excluded because it is the derived-budget sentinel above.
+TEST(StanMath_integrate_1d_gk_prim, work_is_monotone_in_positive_abs_tol) {
+  auto evaluations_for = [](double absolute_tolerance) {
+    int evaluations = 0;
+    auto integrand = [&evaluations](double x, double xc, std::ostream *msgs) {
+      ++evaluations;
+      return std::exp(-x * x) * std::cos(30 * x);
+    };
+    try {
+      stan::math::integrate_1d_gauss_kronrod_tol(integrand, 0.0, 3.0, 1e-12,
+                                                 absolute_tolerance, 12,
+                                                 integrate_1d_gk_test::msgs);
+    } catch (const std::domain_error &) {
+      // Convergence failure is irrelevant here; only the cost is.
+    }
+    return evaluations;
+  };
+
+  int previous = evaluations_for(1e-300);
+  for (double absolute_tolerance : {1e-16, 1e-12, 1e-8, 1e-4, 1e-1}) {
+    const int current = evaluations_for(absolute_tolerance);
+    EXPECT_LE(current, previous) << "abs_tol = " << absolute_tolerance;
+    previous = current;
+  }
+}
+
+// A positive abs_tol must not paper over a genuinely unconverged result.
+TEST(StanMath_integrate_1d_gk_prim, positive_abs_tol_still_throws_when_needed) {
+  auto integrand = [](double x, double xc, std::ostream *msgs) {
+    return std::pow(x, -0.9);
+  };
+  EXPECT_THROW(
+      stan::math::integrate_1d_gauss_kronrod_tol(
+          integrand, 1e-300, 1.0, 1e-12, 1e-14, 15, integrate_1d_gk_test::msgs),
+      std::domain_error);
+}
+
+// max_depth = 0 means "one panel, no bisection"; a positive abs_tol must not
+// disturb that (the refinement floor is never consulted).
+TEST(StanMath_integrate_1d_gk_prim, positive_abs_tol_with_max_depth_zero) {
+  auto integrand
+      = [](double x, double xc, std::ostream *msgs) { return std::exp(x); };
+  const double Q = stan::math::integrate_1d_gauss_kronrod_tol(
+      integrand, 0.0, 1.0, 1e-10, 1e-12, 0, integrate_1d_gk_test::msgs);
+  EXPECT_NEAR(Q, std::exp(1.0) - 1.0, 1e-12);
+}
+
+// A positive abs_tol must agree with the abs_tol == 0 answer under every
+// change of variable, not just the finite one.
+TEST(StanMath_integrate_1d_gk_prim, positive_abs_tol_domain_transformations) {
+  constexpr double relative_tolerance = 1e-12;
+  constexpr double absolute_tolerance = 1e-10;
+  constexpr int max_depth = 8;
+  const double infinity = std::numeric_limits<double>::infinity();
+
+  auto check_integral = [&](const auto &integrand, double lower, double upper,
+                            double expected) {
+    const double legacy_result = stan::math::integrate_1d_gauss_kronrod_tol(
+        integrand, lower, upper, relative_tolerance, 0.0, max_depth,
+        integrate_1d_gk_test::msgs);
+    const double absolute_result = stan::math::integrate_1d_gauss_kronrod_tol(
+        integrand, lower, upper, relative_tolerance, absolute_tolerance,
+        max_depth, integrate_1d_gk_test::msgs);
+    EXPECT_NEAR(absolute_result, expected, absolute_tolerance);
+    EXPECT_NEAR(absolute_result, legacy_result, absolute_tolerance);
+  };
+
+  auto increasing_exponential
+      = [](double x, double xc, std::ostream *msgs) { return std::exp(x); };
+  auto decreasing_exponential
+      = [](double x, double xc, std::ostream *msgs) { return std::exp(-x); };
+  auto gaussian_kernel = [](double x, double xc, std::ostream *msgs) {
+    return std::exp(-x * x);
+  };
+
+  check_integral(increasing_exponential, 0.0, 1.0, std::exp(1.0) - 1.0);
+  check_integral(decreasing_exponential, 0.0, infinity, 1.0);
+  check_integral(increasing_exponential, -infinity, 0.0, 1.0);
+  check_integral(gaussian_kernel, -infinity, infinity,
+                 std::sqrt(stan::math::pi()));
+}

test/unit/math/rev/functor/integrate_1d_gauss_kronrod_test.cpp

@@ -481,4 +481,52 @@ TEST_F(AgradRev, StanMath_integrate_1d_gk_rev_TestUniform) {
   EXPECT_FLOAT_EQ(1, 1 + g[1]);
 }
 
+// Reverse-mode counterpart of
+// StanMath_integrate_1d_gk_prim.positive_abs_tol_reduces_work_on_negligible_
+// integrand. The gradient is computed by integrating the adjoint, so the
+// refinement floor has to hold for that pass too -- and in practice that is
+// where it matters most, because a parameter derivative is routinely orders
+// of magnitude smaller than the value it differentiates.
+TEST_F(AgradRev, StanMath_integrate_1d_gk_rev_positive_abs_tol_reduces_work) {
+  using stan::math::var;
+
+  constexpr double scale = 1e-12;
+  constexpr double frequency = 127.0;
+  constexpr double absolute_tolerance = 1e-14;
+  const double expected_adjoint
+      = scale * (1.0 - std::cos(frequency)) / frequency;
+  const double theta_value = 1.7;
+
+  auto run = [&](double abs_tol, int &evaluations, double &value,
+                 double &adjoint) {
+    stan::math::nested_rev_autodiff nested;
+    auto integrand = [&evaluations](double x, double xc, std::ostream *msgs,
+                                    const auto &theta) {
+      ++evaluations;
+      return theta[0] * scale * stan::math::sin(frequency * x);
+    };
+    std::vector<var> theta{theta_value};
+    var integral = stan::math::integrate_1d_gauss_kronrod_tol(
+        integrand, 0.0, 1.0, 1e-12, abs_tol, 5, msgs, theta);
+    integral.grad();
+    value = integral.val();
+    adjoint = theta[0].adj();
+  };
+
+  int relative_evaluations = 0, absolute_evaluations = 0;
+  double relative_value = 0, relative_adjoint = 0;
+  double absolute_value = 0, absolute_adjoint = 0;
+  run(0.0, relative_evaluations, relative_value, relative_adjoint);
+  run(absolute_tolerance, absolute_evaluations, absolute_value,
+      absolute_adjoint);
+
+  EXPECT_NEAR(absolute_value, theta_value * expected_adjoint,
+              absolute_tolerance);
+  EXPECT_NEAR(absolute_adjoint, expected_adjoint, absolute_tolerance);
+  EXPECT_NEAR(absolute_value, relative_value, absolute_tolerance);
+  EXPECT_NEAR(absolute_adjoint, relative_adjoint, absolute_tolerance);
+
+  EXPECT_LT(absolute_evaluations, relative_evaluations);
+}
+
 }  // namespace integrate_1d_gk_test