Use boost beta, simplify overloads by andrjohns · Pull Request #3212 · stan-dev/math

andrjohns · 2025-06-25T14:23:29Z

Summary

This PR replaces our manual implementation of the beta function with the one from Boost which is more accurate over a broader domain

I've also combined the multiple overloads into a single overload per prim, rev, fwd now that we have if constexpr.

Tests

N/A

Side Effects

N/A

Release notes

Use Boost's implementation of beta for improved stability

Checklist

Copyright holder: Andrew Johnson

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

andrjohns · 2025-06-25T14:26:49Z

@SteveBronder what do you think of this structure for combining overloads with if constexpr? Just figuring out a good consistent structure to follow

SteveBronder · 2025-10-20T16:11:21Z

stan/math/fwd/fun/beta.hpp

+  const FvarInnerT beta_val = beta(a_val, b_val);
+  const FvarInnerT digamma_ab = digamma(a_val + b_val);
+  FvarInnerT beta_d(0);
+  if constexpr (!is_constant<Ta>::value) {


Suggested change

if constexpr (!is_constant<Ta>::value) {

if constexpr (is_autodiff_v<Ta>) {

SteveBronder · 2025-10-20T16:12:01Z

stan/math/fwd/fun/beta.hpp

+  if constexpr (!is_constant<Ta>::value) {
+    beta_d += (digamma(a_val) - digamma_ab) * beta_val * a.d_;
+  }
+  if constexpr (!is_constant<Tb>::value) {


Suggested change

if constexpr (!is_constant<Tb>::value) {

if constexpr (is_autodiff_v<Tb>) {

SteveBronder · 2025-10-20T16:36:50Z

stan/math/prim/fun/beta.hpp

+template <
+    typename T1, typename T2, require_any_container_t<T1, T2>* = nullptr,
+    require_t<math::disjunction<
+        is_arithmetic<return_type_t<T1, T2>>, is_fvar<return_type_t<T1, T2>>,
+        is_std_vector<T1>, is_std_vector<T2>>>* = nullptr>


Little confused. Wouldn't is_fvar<return_type> conflict with the fvar version of beta?

SteveBronder · 2025-10-20T16:38:56Z

stan/math/rev/fun/beta.hpp

+inline auto beta(const T1& a, const T2& b) {
+  arena_t<ref_type_t<T1>> arena_a = a;
+  arena_t<ref_type_t<T2>> arena_b = b;


Suggested change

inline auto beta(const T1& a, const T2& b) {

arena_t<ref_type_t<T1>> arena_a = a;

arena_t<ref_type_t<T2>> arena_b = b;

inline auto beta(T1&& a, T2&& b) {

arena_t<ref_type_t<T1>> arena_a = std::forward<T1>(a);

arena_t<ref_type_t<T2>> arena_b = std::forward<T2>(b);

SteveBronder · 2025-10-20T16:40:09Z

stan/math/rev/fun/beta.hpp

+  const auto& beta_val = beta(value_of(arena_a), value_of(arena_b));
+  using return_type_t = return_var_matrix_t<decltype(beta_val), T1, T2>;
+  arena_t<return_type_t> res(beta_val);


Suggested change

const auto& beta_val = beta(value_of(arena_a), value_of(arena_b));

using return_type_t = return_var_matrix_t<decltype(beta_val), T1, T2>;

arena_t<return_type_t> res(beta_val);

using inner_ret_t = decltype(beta(value_of(arena_a), value_of(arena_b)));

using return_type_t = return_var_matrix_t<inner_ret_t, T1, T2>;

arena_t<return_type_t> res(beta(value_of(arena_a), value_of(arena_b)));

SteveBronder · 2025-10-20T16:45:23Z

stan/math/rev/fun/beta.hpp

+    auto&& b_array = as_array_or_scalar(arena_b);
+    const auto& res_array = as_array_or_scalar(res);
+    const auto& digamma_ab = digamma(value_of(a_array) + value_of(b_array));
+    const auto& adj_val = res_array.adj() * res_array.val();


This produces a new expression objet so we don't need a const& of it

Suggested change

const auto& adj_val = res_array.adj() * res_array.val();

auto adj_val = res_array.adj() * res_array.val();

Anywhere that we will produce a new Eigen expression it should just be an object.

Same for several places here

SteveBronder · 2025-10-20T16:48:02Z

stan/math/rev/fun/beta.hpp

+  reverse_pass_callback([arena_a, arena_b, res]() mutable {
+    auto&& a_array = as_array_or_scalar(arena_a);
+    auto&& b_array = as_array_or_scalar(arena_b);
+    const auto& res_array = as_array_or_scalar(res);


You can use decltype(auto) here since as_array_or_scalar. decltype(auto) is allowed to be either a reference to an object or an object. The function as_array_or_scalar can either return a reference, for instance for a scalar or if an array is passed in as input and then given as output. It can also return a new object, like when we are return an array view of a matrix. That array view is an actual object.

SteveBronder · 2025-10-20T16:48:27Z

stan/math/rev/fun/beta.hpp

-  auto digamma_ab = (digamma(b.val()) - digamma(a + b.val())) * beta_val;
-  return make_callback_var(beta_val, [b, digamma_ab](auto& vi) mutable {
-    b.adj() += vi.adj() * digamma_ab;
+    if constexpr (!is_constant<T1>::value) {


It would be nice to use the new is_autodiff_v everywhere you use !is_constant

SteveBronder · 2025-10-20T16:49:05Z

test/unit/math/prim/fun/beta_test.cpp

+  Eigen::VectorXd in1 = Eigen::VectorXd::Random(6).cwiseAbs();
+  Eigen::VectorXd in2 = Eigen::VectorXd::Random(6).cwiseAbs();


Why did this have to change?

andrjohns and others added 2 commits June 25, 2025 22:14

Use boost beta, simplify overloads

e5e2f13

[Jenkins] auto-formatting by clang-format version 10.0.0-4ubuntu1

544a5e8

andrjohns added 4 commits June 25, 2025 22:35

Update doxygen

2c81206

Expression fixes

87f0309

arena for ret

f691622

Merge branch 'develop' into simplify-beta

0cbdf91

SteveBronder requested changes Oct 20, 2025

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Use boost beta, simplify overloads#3212

Use boost beta, simplify overloads#3212
andrjohns wants to merge 6 commits intodevelopfrom
simplify-beta

andrjohns commented Jun 25, 2025

Uh oh!

andrjohns commented Jun 25, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

SteveBronder Oct 20, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

	if constexpr (!is_constant<Ta>::value) {
	if constexpr (is_autodiff_v<Ta>) {

	if constexpr (!is_constant<Tb>::value) {
	if constexpr (is_autodiff_v<Tb>) {

	const auto& adj_val = res_array.adj() * res_array.val();
	auto adj_val = res_array.adj() * res_array.val();

		Eigen::VectorXd in1 = Eigen::VectorXd::Random(6).cwiseAbs();
		Eigen::VectorXd in2 = Eigen::VectorXd::Random(6).cwiseAbs();

Uh oh!

Conversation

andrjohns commented Jun 25, 2025

Summary

Tests

Side Effects

Release notes

Checklist

Uh oh!

andrjohns commented Jun 25, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants