wiener_lpdf reverse-mode derivative wrt w disagrees with finite differences of Stan Math value

Summary

stan::math::wiener_lpdf appears to return an incorrect reverse-mode adjoint for the relative starting point parameter w.

The log-density value itself appears smooth and internally consistent, but the reverse-mode derivative w.adj() disagrees with central finite differences of Stan Math’s own scalar wiener_lpdf value function.

The mismatch appears in both:

stan::math::wiener_lpdf(y, a, t0, w, v, sv)

and

stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0)

The issue also persists when sv = 0, so it does not appear to be limited to the drift-variability correction.

Environment

Stan Math version: v5.2.0
Stan Math commit: 4c5180a0d9bf0e686fdd9e71ad609fb8ffbf3b27
Compiler: g++ (GCC) 15.2.0
Compiler flags: -std=c++17 -O3 -D_REENTRANT \
  -I ~/cmdstan/stan/lib/stan_math/lib/eigen_3.4.0 \
  -I ~/cmdstan/stan/lib/stan_math/lib/boost_1.78.0 \
  -I ~/cmdstan/stan/lib/stan_math/lib/sundials_6.1.1/include \
  -I ~/cmdstan/stan/lib/stan_math/ \
  -I ~/cmdstan/stan/lib/stan_math/lib/tbb_2020.3/include \
  wiener_w_grad_bug.cpp -o wiener_w_grad_bug \
  -L ~/cmdstan/stan/lib/stan_math/lib/tbb/ \
  -ltbb
OS: NixOS 26.05 (Yarara)
CmdStan commit: 1ddbddeab633b9d0e4620d1d425325e49e3695c2

Minimal reproducer

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stan/math.hpp>

using stan::math::var;

template <typename F>
double central_diff(F&& f, double x, double h) {
  return (f(x + h) - f(x - h)) / (2.0 * h);
}

double lp5_double(double y, double a, double t0, double w, double v,
                  double sv) {
  return stan::math::wiener_lpdf(y, a, t0, w, v, sv);
}

double lp_full_double(double y, double a, double t0, double w, double v,
                      double sv, double sw, double st0) {
  return stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0);
}

void ad_5param(double y, double a, double t0, double w, double v,
               double sv) {
  var yv = y;
  var av = a;
  var t0v = t0;
  var wv = w;
  var vv = v;
  var svv = sv;

  auto lp = stan::math::wiener_lpdf(yv, av, t0v, wv, vv, svv);
  lp.grad();

  std::cout << "AD_5PARAM"
            << " lp=" << lp.val()
            << " gy=" << yv.adj()
            << " ga=" << av.adj()
            << " gt0=" << t0v.adj()
            << " gw=" << wv.adj()
            << " gv=" << vv.adj()
            << " gsv=" << svv.adj()
            << "\n";

  stan::math::recover_memory();
}

void ad_full(double y, double a, double t0, double w, double v, double sv,
             double sw, double st0) {
  var yv = y;
  var av = a;
  var t0v = t0;
  var wv = w;
  var vv = v;
  var svv = sv;
  var swv = sw;
  var st0v = st0;

  auto lp = stan::math::wiener_lpdf(yv, av, t0v, wv, vv, svv, swv, st0v);
  lp.grad();

  std::cout << "AD_FULL"
            << " lp=" << lp.val()
            << " gy=" << yv.adj()
            << " ga=" << av.adj()
            << " gt0=" << t0v.adj()
            << " gw=" << wv.adj()
            << " gv=" << vv.adj()
            << " gsv=" << svv.adj()
            << " gsw=" << swv.adj()
            << " gst0=" << st0v.adj()
            << "\n";

  stan::math::recover_memory();
}

int main() {
  std::cout << std::setprecision(17);

  const double y = 6.0;
  const double a = 10.0;
  const double t0 = 0.01;
  const double w = 0.1;
  const double v = -3.0;
  const double sv = 0.2;
  const double sw = 0.1;
  const double st0 = 0.0;

  std::cout << "=== five parameter: sv = 0.2 ===\n";
  ad_5param(y, a, t0, w, v, sv);

  for (double h : {1e-3, 1e-4, 1e-5, 1e-6, 1e-7}) {
    const double fd = central_diff(
        [&](double ww) { return lp5_double(y, a, t0, ww, v, sv); },
        w,
        h);
    std::cout << "FD_5PARAM h=" << h << " gw_fd=" << fd << "\n";
  }

  for (double ww : {0.08, 0.09, 0.10, 0.11, 0.12}) {
    std::cout << "GRID_5PARAM w=" << ww
              << " lp=" << lp5_double(y, a, t0, ww, v, sv)
              << "\n";
  }

  std::cout << "\n=== five parameter: sv = 0 control ===\n";
  ad_5param(y, a, t0, w, v, 0.0);

  for (double h : {1e-3, 1e-4, 1e-5, 1e-6, 1e-7}) {
    const double fd = central_diff(
        [&](double ww) { return lp5_double(y, a, t0, ww, v, 0.0); },
        w,
        h);
    std::cout << "FD_5PARAM_SV0 h=" << h << " gw_fd=" << fd << "\n";
  }

  std::cout << "\n=== full parameter: sw = 0.1, st0 = 0 ===\n";
  ad_full(y, a, t0, w, v, sv, sw, st0);

  for (double h : {1e-3, 1e-4, 1e-5, 1e-6, 1e-7}) {
    const double fd = central_diff(
        [&](double ww) {
          return lp_full_double(y, a, t0, ww, v, sv, sw, st0);
        },
        w,
        h);
    std::cout << "FD_FULL h=" << h << " gw_fd=" << fd << "\n";
  }

  for (double ww : {0.08, 0.09, 0.10, 0.11, 0.12}) {
    std::cout << "GRID_FULL w=" << ww
              << " lp=" << lp_full_double(y, a, t0, ww, v, sv, sw, st0)
              << "\n";
  }

  return 0;
}

Observed output

=== five parameter: sv = 0.2 ===
AD_5PARAM lp=-50.539106208790059 gy=-1.4309088946356134 ga=-3.1387340804639763 gt0=1.4309088946356134 gw=-82.160098301170336 gv=21.757018393030009 gsv=93.707129083578479
FD_5PARAM h=0.001 gw_fd=36.633174261808676
FD_5PARAM h=0.0001 gw_fd=36.632911247345135
FD_5PARAM h=1.0000000000000001e-05 gw_fd=36.632908617661997
FD_5PARAM h=9.9999999999999995e-07 gw_fd=36.632908592793001
FD_5PARAM h=9.9999999999999995e-08 gw_fd=36.632908617661997
GRID_5PARAM w=0.080000000000000002 lp=-51.286630294961604
GRID_5PARAM w=0.089999999999999997 lp=-50.908794872521725
GRID_5PARAM w=0.10000000000000001 lp=-50.539106208790059
GRID_5PARAM w=0.11 lp=-50.175601258770968
GRID_5PARAM w=0.12 lp=-49.81692817570945

=== five parameter: sv = 0 control ===
AD_5PARAM lp=-62.167447556524714 gy=-3.6469225854648935 ga=-3.9219867426130568 gt0=3.6469225854648935 gw=-73.457290943958327 gv=26.969999999999999 gsv=0
FD_5PARAM_SV0 h=0.001 gw_fd=45.335981619022192
FD_5PARAM_SV0 h=0.0001 gw_fd=45.335718604526676
FD_5PARAM_SV0 h=1.0000000000000001e-05 gw_fd=45.335715974559314
FD_5PARAM_SV0 h=9.9999999999999995e-07 gw_fd=45.335715952887767
FD_5PARAM_SV0 h=9.9999999999999995e-08 gw_fd=45.33571594578234

=== full parameter: sw = 0.1, st0 = 0 ===
AD_FULL lp=-50.058755379483443 gy=-1.5006905642810147 ga=-3.0292867403094825 gt0=1.5006905642810147 gw=-34.586923958070955 gv=21.554056380238126 gsv=91.955143887665372 gsw=8.7239803669499629 gst0=0
FD_FULL h=0.001 gw_fd=35.7193384750083
FD_FULL h=0.0001 gw_fd=35.719188331704288
FD_FULL h=1.0000000000000001e-05 gw_fd=35.719186830363014
FD_FULL h=9.9999999999999995e-07 gw_fd=35.719186815441617
FD_FULL h=9.9999999999999995e-08 gw_fd=35.719186861626895
GRID_FULL w=0.080000000000000002 lp=-50.783094978471475
GRID_FULL w=0.089999999999999997 lp=-50.418236984588823
GRID_FULL w=0.10000000000000001 lp=-50.058755379483443
GRID_FULL w=0.11 lp=-49.703547959822139
GRID_FULL w=0.12 lp=-49.351836459209828

Expected behavior

The reverse-mode adjoint for w should agree with central finite differences of the same scalar value function.

For the 5-parameter case with sv = 0.2, the finite-difference derivative is stable around:

gw ~ +36.63290859

but reverse-mode reports:

gw = -82.160098301170336

For the sv = 0 control, the finite-difference derivative is stable around:

gw ~ +45.33571595

but reverse-mode reports:

gw = -73.457290943958327

For the full-parameter case with sw = 0.1, st0 = 0, the finite-difference derivative is stable around:

gw ~ +35.719186815

but reverse-mode reports:

gw = -34.586923958070955

Possibility of a bug

For a scalar differentiable log density

$$l(w) = \log p(y\ \bar\ a, t_0, w, v, s_v),$$

the reverse-mode adjoint should satisfy

$$w.adj()\ \sim\ \frac{l(w + h) - l(w - h)}{2h}.$$

This is not a comparison against an external implementation. (However, that's how I discovered the issue, as Stan tests for a WienR reference implementation which appears to be also affected). The finite differences are computed from Stan Math's own scalaer wiener_lpdf value.

The value grids are locally monotone increasing in w, for example:

five-param:
w=0.08 lp=-51.286630294961604
w=0.10 lp=-50.539106208790059
w=0.12 lp=-49.81692817570945

and

full-param:
w=0.08 lp=-50.783094978471475
w=0.10 lp=-50.058755379483443
w=0.12 lp=-49.351836459209828

so the local derivate with respect to w should be positive. The finite differences are stable over several step sizes, but the reverse-mode adjoint is negative and of substantially different magnitude.

The sv = 0 control suggests this is not solely cause by the drift-variability term.

Location

The likely bug path is internal::wiener5_grad_w in stan/math/prim/prob/wiener5_lpdf.hpp. The full overload delegates the no-sw, no-st0 case to the 5-parameter overload, and for the full sw > 0 case it integrates internal::wiener5_grad_w for the w partial. The source page shows the full overload delegating to wiener_lpdf(..., sv, precision_derivatives) when sw == 0 && st0 == 0, and using internal::wiener5_grad_w inside the T_w partial path for the full case.

So the fix should start at the 5-parameter w derivative, not in cubature or the full wrapper.

Related tests

The existing Wiener full test appears to compare gradients against constants generated externally, for example from WienR. this issue suggests that at least the w / beta adjoint should also be tested against finite differences of Stan Math's own value function.

A possible regression test would check w.adj() against central finite differences for:

stan::math::wiener_lpdf(y, a, t0, w, v, sv)

and

stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0)

at the parameter values above.

Additional note

For the full case with sw > 0 and st == 0, there is also a simple endpoint identity for the derivative with respect to the center w0 of the uniform starting-point variability interval:

$$F(w_0) = \frac{1}{s_w} \int_{w_0 - s_w / 2}^{w_0 + s_w / 2} f(u) \,du,$$

where $f(u)$ is the 5-parameter density evaluated at starting point $u$.
Therefore

$$l(w_0) = \log F(w_0)$$

has derivative

$$\frac{\partial l}{\partial w_0} = \frac{f(w_0 + s_w / 2) - f(w_0 - s_w / 2)}{s_w F(w_0)}.$$

This identity agrees with the finite-difference sign and magnitude in an independent implementation, not with the reverse-mode adjoint shown above.

@SteveBronder

ping @SteveBronder @Franzi2114

The existing prim full-Wiener test still passes because it compares against the current reference constants.

A new test I am running locally within repo:

------------------------------------------------------------
test/unit/math/rev/prob/wiener_w_gradient_fd_test --gtest_output="xml:test/unit/math/rev/prob/wiener_w_gradient_fd_test.xml"
Running main() from lib/benchmark_1.5.1/googletest/googletest/src/gtest_main.cc
[==========] Running 6 tests from 1 test suite.
[----------] Global test environment set-up.
[----------] 6 tests from MathRevProbWienerWGradientFD
[ RUN      ] MathRevProbWienerWGradientFD.FiveParamI4AdDisagreesWithOwnValueFD

CASE_5 five_param_i4: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0, st0=0
AD_5 lp=-50.539106208790059 gw=-82.160098301170336
FD_5 h=0.001 gw_fd=36.633174261808676 diff_ad_minus_fd=-118.79327256297901
FD_5 h=0.0001 gw_fd=36.632911247345135 diff_ad_minus_fd=-118.79300954851547
FD_5 h=1.0000000000000001e-05 gw_fd=36.632908617661997 diff_ad_minus_fd=-118.79300691883233
FD_5 h=9.9999999999999995e-07 gw_fd=36.632908592793001 diff_ad_minus_fd=-118.79300689396334
FD_5 h=9.9999999999999995e-08 gw_fd=36.632908617661997 diff_ad_minus_fd=-118.79300691883233
GRID_5 w=0.080000000000000002 lp=-51.286630294961604
GRID_5 w=0.090000000000000011 lp=-50.908794872521725
GRID_5 w=0.10000000000000001 lp=-50.539106208790059
GRID_5 w=0.11 lp=-50.175601258770968
GRID_5 w=0.12000000000000001 lp=-49.81692817570945
test/unit/math/rev/prob/wiener_w_gradient_fd_test.cpp:219: Failure
The difference between ad.gw and fd is 118.79300689396334, which exceeds 1e-5, where
ad.gw evaluates to -82.160098301170336,
fd evaluates to 36.632908592793001, and
1e-5 evaluates to 1.0000000000000001e-05.
Stan reverse-mode gw disagrees with central finite difference of Stan Math's own wiener_lpdf(y,a,t0,w,v,sv). five_param_i4: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0, st0=0
[  FAILED  ] MathRevProbWienerWGradientFD.FiveParamI4AdDisagreesWithOwnValueFD (0 ms)
[ RUN      ] MathRevProbWienerWGradientFD.FiveParamZeroSvControlAlsoFails

CASE_5 five_param_i4_sv_zero: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0, sw=0, st0=0
AD_5 lp=-62.167447556524714 gw=-73.457290943958327
FD_5 h=0.001 gw_fd=45.335981619022192 diff_ad_minus_fd=-118.79327256298052
FD_5 h=0.0001 gw_fd=45.335718604526676 diff_ad_minus_fd=-118.793009548485
FD_5 h=1.0000000000000001e-05 gw_fd=45.335715974559314 diff_ad_minus_fd=-118.79300691851765
FD_5 h=9.9999999999999995e-07 gw_fd=45.335715952887767 diff_ad_minus_fd=-118.79300689684609
FD_5 h=9.9999999999999995e-08 gw_fd=45.33571594578234 diff_ad_minus_fd=-118.79300688974067
GRID_5 w=0.080000000000000002 lp=-63.089673159314529
GRID_5 w=0.090000000000000011 lp=-62.624325636197014
GRID_5 w=0.10000000000000001 lp=-62.167447556524714
GRID_5 w=0.11 lp=-61.717075875302008
GRID_5 w=0.12000000000000001 lp=-61.271858745773891
test/unit/math/rev/prob/wiener_w_gradient_fd_test.cpp:237: Failure
The difference between ad.gw and fd is 118.79300689684609, which exceeds 1e-5, where
ad.gw evaluates to -73.457290943958327,
fd evaluates to 45.335715952887767, and
1e-5 evaluates to 1.0000000000000001e-05.
The w-adjoint mismatch persists even at sv = 0. five_param_i4_sv_zero: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0, sw=0, st0=0
[  FAILED  ] MathRevProbWienerWGradientFD.FiveParamZeroSvControlAlsoFails (0 ms)
[ RUN      ] MathRevProbWienerWGradientFD.FullI4AdDisagreesWithOwnValueFD

CASE_FULL full_i4_sw_positive: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0.10000000000000001, st0=0
AD_FULL lp=-50.058755379483443 gw=-34.586923958070955
FD_FULL h=0.001 gw_fd=35.7193384750083 diff_ad_minus_fd=-70.306262433079254
FD_FULL h=0.0001 gw_fd=35.719188331704288 diff_ad_minus_fd=-70.306112289775243
FD_FULL h=1.0000000000000001e-05 gw_fd=35.719186830363014 diff_ad_minus_fd=-70.306110788433969
FD_FULL h=9.9999999999999995e-07 gw_fd=35.719186815441617 diff_ad_minus_fd=-70.306110773512572
FD_FULL h=9.9999999999999995e-08 gw_fd=35.719186861626895 diff_ad_minus_fd=-70.306110819697849
GRID_FULL w=0.080000000000000002 lp=-50.783094978471475
GRID_FULL w=0.090000000000000011 lp=-50.418236984588823
GRID_FULL w=0.10000000000000001 lp=-50.058755379483443
GRID_FULL w=0.11 lp=-49.703547959822139
GRID_FULL w=0.12000000000000001 lp=-49.351836459209821
test/unit/math/rev/prob/wiener_w_gradient_fd_test.cpp:253: Failure
The difference between ad.gw and fd is 70.306110773512572, which exceeds 1e-5, where
ad.gw evaluates to -34.586923958070955,
fd evaluates to 35.719186815441617, and
1e-5 evaluates to 1.0000000000000001e-05.
Stan reverse-mode gw disagrees with central finite difference of Stan Math's own full wiener_lpdf value function. full_i4_sw_positive: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0.10000000000000001, st0=0
[  FAILED  ] MathRevProbWienerWGradientFD.FullI4AdDisagreesWithOwnValueFD (1 ms)
[ RUN      ] MathRevProbWienerWGradientFD.FullI4EndpointIdentity

ENDPOINT_IDENTITY full_i4_endpoint_identity: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0.10000000000000001, st0=0
AD_FULL gw=-34.586923958070955
FD_FULL gw=35.719186815441617
ENDPOINT gw=35.719186815203571
test/unit/math/rev/prob/wiener_w_gradient_fd_test.cpp:282: Failure
The difference between ad.gw and endpoint is 70.306110773274526, which exceeds 1e-5, where
ad.gw evaluates to -34.586923958070955,
endpoint evaluates to 35.719186815203571, and
1e-5 evaluates to 1.0000000000000001e-05.
For sw > 0 and st0 == 0, d/dw has a closed endpoint identity. Stan AD disagrees with both that identity and finite differences. full_i4_endpoint_identity: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0.10000000000000001, st0=0
[  FAILED  ] MathRevProbWienerWGradientFD.FullI4EndpointIdentity (0 ms)
[ RUN      ] MathRevProbWienerWGradientFD.ExistingFullRowsCompareWAdjointToFD
ROW_CHECK row_0: y=2, a=2, t0=1.0000000000000001e-09, w=0.10000000000000001, v=2, sv=0, sw=0, st0=0 lp=-4.2856474786661476 ad_gw=5.6712018451731829 fd_gw=5.6712018450255641 diff=1.4761880606783961e-10
ROW_CHECK row_1: y=3, a=2, t0=0.01, w=0.5, v=2, sv=0.20000000000000001, sw=0, st0=0 lp=-7.5237923514690932 ad_gw=-3.6439622109007557 fd_gw=-3.6439622110862047 diff=1.8544898949812705e-10
ROW_CHECK row_2: y=4, a=10, t0=0.01, w=0.80000000000000004, v=4, sv=0, sw=0.10000000000000001, st0=0 lp=-26.155105620924822 ad_gw=-38.777505714679201 fd_gw=-38.777505716680594 diff=2.0013928292428318e-09
ROW_CHECK row_3: y=5, a=4, t0=0.01, w=0.69999999999999996, v=3, sv=0, sw=0, st0=0.0070000000000000001 lp=-22.193913489208871 ad_gw=-14.183853955927921 fd_gw=-14.183853954463643 diff=-1.4642775880702175e-09
ROW_CHECK row_4: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0.10000000000000001, st0=0 lp=-50.058755379483443 ad_gw=-34.586923958070955 fd_gw=35.719186815441617 diff=-70.306110773512572
test/unit/math/rev/prob/wiener_w_gradient_fd_test.cpp:318: Failure
The difference between ad.gw and fd is 70.306110773512572, which exceeds 1e-4, where
ad.gw evaluates to -34.586923958070955,
fd evaluates to 35.719186815441617, and
1e-4 evaluates to 0.0001.
w-adjoint does not match central finite difference for existing Wiener test row.
Google Test trace:
test/unit/math/rev/prob/wiener_w_gradient_fd_test.cpp:292: row_4: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0.10000000000000001, st0=0
ROW_CHECK row_5: y=7, a=1, t0=0.01, w=0.90000000000000002, v=1, sv=0.20000000000000001, sw=0, st0=0.0070000000000000001 lp=-37.281726358631786 ad_gw=-10.453534568194648 fd_gw=-10.453534571297496 diff=3.1028477565087087e-09
ROW_CHECK row_6: y=8, a=3, t0=0.01, w=0.69999999999999996, v=-1, sv=0, sw=0.10000000000000001, st0=0.0070000000000000001 lp=-10.542866207943767 ad_gw=0.6795979835829058 fd_gw=0.67959798322192455 diff=3.6098124489569727e-10
ROW_CHECK row_7: y=8.8499999999999996, a=1.7, t0=0.01, w=0.92000000000000004, v=-7.2999999999999998, sv=0.69999999999999996, sw=0.01, st0=0.0089999999999999993 lp=-61.591590567424603 ad_gw=-9.9314454083420163 fd_gw=-9.9314454082843895 diff=-5.7626792226983525e-11
ROW_CHECK row_8: y=8.9000000000000004, a=2.3999999999999999, t0=0.01, w=0.90000000000000002, v=-4.9000000000000004, sv=0, sw=0, st0=0.0089999999999999993 lp=-117.23896795933561 ad_gw=2.0911720085754313 fd_gw=2.0911720071126183 diff=1.462812981856132e-09
ROW_CHECK row_9: y=9, a=11, t0=0.01, w=0.12, v=4.5, sv=0.69999999999999996, sw=0.10000000000000001, st0=0.0089999999999999993 lp=-12.578859507920518 ad_gw=-6.0858545466501583 fd_gw=-6.085851879333859 diff=-2.6673162993873234e-06
ROW_CHECK row_10: y=1, a=1.5, t0=0.10000000000000001, w=0.5, v=3, sv=0.5, sw=0.20000000000000001, st0=0 lp=-3.1448097740734999 ad_gw=-3.7487031097808319 fd_gw=-3.7487031097072077 diff=-7.3624217833412331e-11
[  FAILED  ] MathRevProbWienerWGradientFD.ExistingFullRowsCompareWAdjointToFD (15 ms)
[ RUN      ] MathRevProbWienerWGradientFD.FiveAndFullZeroValueAndWAdjointSamePath

FIVE_VS_FULL_ZERO five_vs_full_zero: y=6, a=10, t0=0.01, w=0.10000000000000001, v=-3, sv=0.20000000000000001, sw=0, st0=0
AD_5 lp=-50.539106208790059 gw=-82.160098301170336
AD_FULL_ZERO lp=-50.539106208790059 gw=-82.160098301170336
[       OK ] MathRevProbWienerWGradientFD.FiveAndFullZeroValueAndWAdjointSamePath (0 ms)
[----------] 6 tests from MathRevProbWienerWGradientFD (16 ms total)

[----------] Global test environment tear-down
[==========] 6 tests from 1 test suite ran. (16 ms total)
[  PASSED  ] 1 test.
[  FAILED  ] 5 tests, listed below:
[  FAILED  ] MathRevProbWienerWGradientFD.FiveParamI4AdDisagreesWithOwnValueFD
[  FAILED  ] MathRevProbWienerWGradientFD.FiveParamZeroSvControlAlsoFails
[  FAILED  ] MathRevProbWienerWGradientFD.FullI4AdDisagreesWithOwnValueFD
[  FAILED  ] MathRevProbWienerWGradientFD.FullI4EndpointIdentity
[  FAILED  ] MathRevProbWienerWGradientFD.ExistingFullRowsCompareWAdjointToFD

 5 FAILED TESTS
test/unit/math/rev/prob/wiener_w_gradient_fd_test --gtest_output="xml:test/unit/math/rev/prob/wiener_w_gradient_fd_test.xml" failed
exit now (05/12/26 17:58:04 CEST)
1

Can you share the test code you are using? Also if you are in the math library it would be nice to use the stan::test::expect_ad framework to show this is happening.

Of course. Sorry, I didn't know about expect_ad. I have confirmed the bug with that as well. I have written a fix and now passing tests, making a PR soon. For the time being, the tests I ran without stan::test::expect_ad:

#include <stan/math/rev.hpp>
#include <gtest/gtest.h>

#include <cmath>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

namespace {

struct Case {
  const char* name;
  double y;
  double a;
  double t0;
  double w;
  double v;
  double sv;
  double sw;
  double st0;
};

struct Ad5Result {
  double lp;
  double gy;
  double ga;
  double gt0;
  double gw;
  double gv;
  double gsv;
};

struct AdFullResult {
  double lp;
  double gy;
  double ga;
  double gt0;
  double gw;
  double gv;
  double gsv;
  double gsw;
  double gst0;
};

template <typename F>
double central_diff(F&& f, double x, double h) {
  return (f(x + h) - f(x - h)) / (2.0 * h);
}

double lp5_double(const Case& c, double w) {
  return stan::math::wiener_lpdf(c.y, c.a, c.t0, w, c.v, c.sv);
}

double lp_full_double(const Case& c, double w) {
  return stan::math::wiener_lpdf(c.y, c.a, c.t0, w, c.v, c.sv, c.sw, c.st0);
}

Ad5Result ad5(const Case& c) {
  using stan::math::var;

  var y = c.y;
  var a = c.a;
  var t0 = c.t0;
  var w = c.w;
  var v = c.v;
  var sv = c.sv;

  var lp = stan::math::wiener_lpdf(y, a, t0, w, v, sv);
  lp.grad();

  Ad5Result out{lp.val(), y.adj(), a.adj(), t0.adj(), w.adj(), v.adj(),
                sv.adj()};

  stan::math::recover_memory();
  return out;
}

AdFullResult ad_full(const Case& c) {
  using stan::math::var;

  var y = c.y;
  var a = c.a;
  var t0 = c.t0;
  var w = c.w;
  var v = c.v;
  var sv = c.sv;
  var sw = c.sw;
  var st0 = c.st0;

  var lp = stan::math::wiener_lpdf(y, a, t0, w, v, sv, sw, st0);
  lp.grad();

  AdFullResult out{lp.val(), y.adj(),  a.adj(),  t0.adj(), w.adj(),
                   v.adj(),  sv.adj(), sw.adj(), st0.adj()};

  stan::math::recover_memory();
  return out;
}

std::string case_string(const Case& c) {
  std::ostringstream ss;
  ss << std::setprecision(17)
     << c.name << ": "
     << "y=" << c.y << ", "
     << "a=" << c.a << ", "
     << "t0=" << c.t0 << ", "
     << "w=" << c.w << ", "
     << "v=" << c.v << ", "
     << "sv=" << c.sv << ", "
     << "sw=" << c.sw << ", "
     << "st0=" << c.st0;
  return ss.str();
}

void print_w_fd_trace_5(const Case& c) {
  const Ad5Result ad = ad5(c);

  std::cout << std::setprecision(17);
  std::cout << "\nCASE_5 " << case_string(c) << "\n";
  std::cout << "AD_5 lp=" << ad.lp << " gw=" << ad.gw << "\n";

  for (double h : {1e-3, 1e-4, 1e-5, 1e-6, 1e-7}) {
    const double fd = central_diff([&](double ww) { return lp5_double(c, ww); },
                                   c.w, h);
    std::cout << "FD_5 h=" << h << " gw_fd=" << fd
              << " diff_ad_minus_fd=" << (ad.gw - fd) << "\n";
  }

  for (double ww : {c.w - 0.02, c.w - 0.01, c.w, c.w + 0.01, c.w + 0.02}) {
    if (ww > 0.0 && ww < 1.0) {
      std::cout << "GRID_5 w=" << ww << " lp=" << lp5_double(c, ww) << "\n";
    }
  }
}

void print_w_fd_trace_full(const Case& c) {
  const AdFullResult ad = ad_full(c);

  std::cout << std::setprecision(17);
  std::cout << "\nCASE_FULL " << case_string(c) << "\n";
  std::cout << "AD_FULL lp=" << ad.lp << " gw=" << ad.gw << "\n";

  for (double h : {1e-3, 1e-4, 1e-5, 1e-6, 1e-7}) {
    const double fd = central_diff(
        [&](double ww) { return lp_full_double(c, ww); }, c.w, h);
    std::cout << "FD_FULL h=" << h << " gw_fd=" << fd
              << " diff_ad_minus_fd=" << (ad.gw - fd) << "\n";
  }

  for (double ww : {c.w - 0.02, c.w - 0.01, c.w, c.w + 0.01, c.w + 0.02}) {
    if (ww > 0.0 && ww < 1.0
        && (c.sw == 0.0 || (ww - c.sw / 2.0 > 0.0 && ww + c.sw / 2.0 < 1.0))) {
      std::cout << "GRID_FULL w=" << ww << " lp=" << lp_full_double(c, ww)
                << "\n";
    }
  }
}

double endpoint_identity_full_sw_only(const Case& c) {
  // Valid when sw > 0 and st0 == 0.
  //
  // L(w0) = log((1 / sw) int_{w0-sw/2}^{w0+sw/2} f(u) du)
  //
  // dL/dw0 = [f(w0 + sw/2) - f(w0 - sw/2)] / [sw * F],
  // where F is the marginal density exp(L).
  const double low = c.w - c.sw / 2.0;
  const double high = c.w + c.sw / 2.0;

  Case low_case = c;
  Case high_case = c;
  Case mid_case = c;

  low_case.sw = 0.0;
  low_case.st0 = 0.0;
  low_case.w = low;

  high_case.sw = 0.0;
  high_case.st0 = 0.0;
  high_case.w = high;

  const double f_low = std::exp(lp5_double(low_case, low));
  const double f_high = std::exp(lp5_double(high_case, high));
  const double marginal_density = std::exp(lp_full_double(mid_case, c.w));

  return (f_high - f_low) / (c.sw * marginal_density);
}

std::vector<Case> stan_existing_rows() {
  return {
      {"row_0", 2.0, 2.0, 1e-9, 0.10, 2.0, 0.0, 0.00, 0.000},
      {"row_1", 3.0, 2.0, 0.01, 0.50, 2.0, 0.2, 0.00, 0.000},
      {"row_2", 4.0, 10.0, 0.01, 0.80, 4.0, 0.0, 0.10, 0.000},
      {"row_3", 5.0, 4.0, 0.01, 0.70, 3.0, 0.0, 0.00, 0.007},
      {"row_4", 6.0, 10.0, 0.01, 0.10, -3.0, 0.2, 0.10, 0.000},
      {"row_5", 7.0, 1.0, 0.01, 0.90, 1.0, 0.2, 0.00, 0.007},
      {"row_6", 8.0, 3.0, 0.01, 0.70, -1.0, 0.0, 0.10, 0.007},
      {"row_7", 8.85, 1.7, 0.01, 0.92, -7.3, 0.7, 0.01, 0.009},
      {"row_8", 8.9, 2.4, 0.01, 0.90, -4.9, 0.0, 0.00, 0.009},
      {"row_9", 9.0, 11.0, 0.01, 0.12, 4.5, 0.7, 0.10, 0.009},
      {"row_10", 1.0, 1.5, 0.10, 0.50, 3.0, 0.5, 0.20, 0.000},
  };
}

}  // namespace

TEST(MathRevProbWienerWGradientFD, FiveParamI4AdDisagreesWithOwnValueFD) {
  const Case c{"five_param_i4", 6.0, 10.0, 0.01, 0.10, -3.0, 0.2,
               0.0, 0.0};

  print_w_fd_trace_5(c);

  const Ad5Result ad = ad5(c);
  const double fd = central_diff([&](double ww) { return lp5_double(c, ww); },
                                 c.w, 1e-6);

  // This is the assertion that should currently fail if the bug is present.
  EXPECT_NEAR(ad.gw, fd, 1e-5)
      << "Stan reverse-mode gw disagrees with central finite difference of "
         "Stan Math's own wiener_lpdf(y,a,t0,w,v,sv). "
      << case_string(c);
}

TEST(MathRevProbWienerWGradientFD, FiveParamZeroSvControlAlsoFails) {
  const Case c{"five_param_i4_sv_zero", 6.0, 10.0, 0.01, 0.10, -3.0,
               0.0, 0.0, 0.0};

  print_w_fd_trace_5(c);

  const Ad5Result ad = ad5(c);
  const double fd = central_diff([&](double ww) { return lp5_double(c, ww); },
                                 c.w, 1e-6);

  // This control distinguishes "sv derivative/prefactor bug" from a more
  // basic w-adjoint bug.
  EXPECT_NEAR(ad.gw, fd, 1e-5)
      << "The w-adjoint mismatch persists even at sv = 0. "
      << case_string(c);
}

TEST(MathRevProbWienerWGradientFD, FullI4AdDisagreesWithOwnValueFD) {
  const Case c{"full_i4_sw_positive", 6.0, 10.0, 0.01, 0.10, -3.0, 0.2,
               0.1, 0.0};

  print_w_fd_trace_full(c);

  const AdFullResult ad = ad_full(c);
  const double fd = central_diff(
      [&](double ww) { return lp_full_double(c, ww); }, c.w, 1e-6);

  // This is the propagated full-parameter failure.
  EXPECT_NEAR(ad.gw, fd, 1e-5)
      << "Stan reverse-mode gw disagrees with central finite difference of "
         "Stan Math's own full wiener_lpdf value function. "
      << case_string(c);
}

TEST(MathRevProbWienerWGradientFD, FullI4EndpointIdentity) {
  const Case c{"full_i4_endpoint_identity", 6.0, 10.0, 0.01, 0.10, -3.0,
               0.2, 0.1, 0.0};

  ASSERT_GT(c.sw, 0.0);
  ASSERT_EQ(c.st0, 0.0);

  const AdFullResult ad = ad_full(c);
  const double fd = central_diff(
      [&](double ww) { return lp_full_double(c, ww); }, c.w, 1e-6);
  const double endpoint = endpoint_identity_full_sw_only(c);

  std::cout << std::setprecision(17)
            << "\nENDPOINT_IDENTITY " << case_string(c) << "\n"
            << "AD_FULL gw=" << ad.gw << "\n"
            << "FD_FULL gw=" << fd << "\n"
            << "ENDPOINT gw=" << endpoint << "\n";

  // This checks that the independent endpoint identity agrees with finite
  // differences of the value function.
  EXPECT_NEAR(endpoint, fd, 1e-5);

  // This is expected to fail if the current Stan adjoint is wrong.
  EXPECT_NEAR(ad.gw, endpoint, 1e-5)
      << "For sw > 0 and st0 == 0, d/dw has a closed endpoint identity. "
         "Stan AD disagrees with both that identity and finite differences. "
      << case_string(c);
}

TEST(MathRevProbWienerWGradientFD, ExistingFullRowsCompareWAdjointToFD) {
  const auto rows = stan_existing_rows();

  for (const auto& c : rows) {
    SCOPED_TRACE(case_string(c));

    // Avoid cases where finite differencing w would violate Stan's support
    // constraint w - sw/2 > 0 and w + sw/2 < 1.
    const double h = 1e-6;
    if (!(c.w - h > 0.0 && c.w + h < 1.0)) {
      continue;
    }
    if (c.sw > 0.0
        && !(c.w - h - c.sw / 2.0 > 0.0
             && c.w + h + c.sw / 2.0 < 1.0)) {
      continue;
    }

    const AdFullResult ad = ad_full(c);
    const double fd = central_diff(
        [&](double ww) { return lp_full_double(c, ww); }, c.w, h);

    std::cout << std::setprecision(17)
              << "ROW_CHECK " << case_string(c)
              << " lp=" << ad.lp
              << " ad_gw=" << ad.gw
              << " fd_gw=" << fd
              << " diff=" << (ad.gw - fd)
              << "\n";

    EXPECT_NEAR(ad.gw, fd, 1e-4)
        << "w-adjoint does not match central finite difference for existing "
           "Wiener test row.";
  }
}

TEST(MathRevProbWienerWGradientFD, FiveAndFullZeroValueAndWAdjointSamePath) {
  const Case c{"five_vs_full_zero", 6.0, 10.0, 0.01, 0.10, -3.0, 0.2,
               0.0, 0.0};

  const Ad5Result ad5_result = ad5(c);
  const AdFullResult ad_full_result = ad_full(c);

  std::cout << std::setprecision(17)
            << "\nFIVE_VS_FULL_ZERO " << case_string(c) << "\n"
            << "AD_5 lp=" << ad5_result.lp
            << " gw=" << ad5_result.gw << "\n"
            << "AD_FULL_ZERO lp=" << ad_full_result.lp
            << " gw=" << ad_full_result.gw << "\n";

  EXPECT_NEAR(ad5_result.lp, ad_full_result.lp, 1e-12);
  EXPECT_NEAR(ad5_result.gw, ad_full_result.gw, 1e-12);
}

The tests I ran using expect_ad:

#include <stan/math/mix.hpp>
#include <test/unit/math/test_ad.hpp>
#include <gtest/gtest.h>

TEST(MathMixProbWienerLpdf, fiveParamWGradientExpectAd) {
  auto f = [](const auto& w) {
    return stan::math::wiener_lpdf(6.0, 10.0, 0.01, w, -3.0, 0.2);
  };

  stan::test::ad_tolerances tols;
  tols.gradient_grad_ = 1e-5;
  stan::test::expect_ad(tols, f, 0.1);
}

TEST(MathMixProbWienerLpdf, fiveParamZeroSvWGradientExpectAd) {
  auto f = [](const auto& w) {
    return stan::math::wiener_lpdf(6.0, 10.0, 0.01, w, -3.0, 0.0);
  };

  stan::test::ad_tolerances tols;
  tols.gradient_grad_ = 1e-5;
  stan::test::expect_ad(tols, f, 0.1);
}

TEST(MathMixProbWienerLpdf, fullParamWGradientExpectAd) {
  auto f = [](const auto& w) {
    return stan::math::wiener_lpdf(6.0, 10.0, 0.01, w, -3.0, 0.2, 0.1,
                                   0.0);
  };

  stan::test::ad_tolerances tols;
  tols.gradient_grad_ = 1e-5;
  stan::test::expect_ad(tols, f, 0.1);
}

TEST(MathMixProbWienerLpdf, existingFullRowsWGradientExpectAd) {
  struct Case {
    const char* name;
    double y;
    double a;
    double t0;
    double w;
    double v;
    double sv;
    double sw;
    double st0;
  };

  const std::vector<Case> cases = {
      {"row_0", 2.0, 2.0, 1e-9, 0.10, 2.0, 0.0, 0.00, 0.000},
      {"row_1", 3.0, 2.0, 0.01, 0.50, 2.0, 0.2, 0.00, 0.000},
      {"row_2", 4.0, 10.0, 0.01, 0.80, 4.0, 0.0, 0.10, 0.000},
      {"row_3", 5.0, 4.0, 0.01, 0.70, 3.0, 0.0, 0.00, 0.007},
      {"row_4", 6.0, 10.0, 0.01, 0.10, -3.0, 0.2, 0.10, 0.000},
      {"row_5", 7.0, 1.0, 0.01, 0.90, 1.0, 0.2, 0.00, 0.007},
      {"row_6", 8.0, 3.0, 0.01, 0.70, -1.0, 0.0, 0.10, 0.007},
      {"row_7", 8.85, 1.7, 0.01, 0.92, -7.3, 0.7, 0.01, 0.009},
      {"row_8", 8.9, 2.4, 0.01, 0.90, -4.9, 0.0, 0.00, 0.009},
      {"row_9", 9.0, 11.0, 0.01, 0.12, 4.5, 0.7, 0.10, 0.009},
      {"row_10", 1.0, 1.5, 0.10, 0.50, 3.0, 0.5, 0.20, 0.000},
  };

  stan::test::ad_tolerances tols;
  tols.gradient_grad_ = 1e-4;

  for (const auto& c : cases) {
    SCOPED_TRACE(c.name);

    auto f = [c](const auto& w) {
      return stan::math::wiener_lpdf(c.y, c.a, c.t0, w, c.v, c.sv, c.sw,
                                     c.st0);
    };

    // The row sweep is intended to check the reverse-mode w adjoint against
    // finite differences. Some full-Wiener rows are not stable enough for
    // higher-order mixed-mode finite-difference checks
    stan::test::expect_ad<true>(tols, f, c.w);
  }
}

Result of in-tree AD test:

test/unit/math/mix/prob/wiener_lpdf_ad_test --gtest_output="xml:test/unit/math/mix/prob/wiener_lpdf_ad_test.xml"
Running main() from lib/benchmark_1.5.1/googletest/googletest/src/gtest_main.cc
[==========] Running 4 tests from 1 test suite.
[----------] Global test environment set-up.
[----------] 4 tests from MathMixProbWienerLpdf
[ RUN      ] MathMixProbWienerLpdf.fiveParamWGradientExpectAd
./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 118.79300689265223, which exceeds tol_val, where
x1 evaluates to 36.632908591481893,
x2 evaluates to -82.160098301170336, and
tol_val evaluates to 0.00059396503446326125.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient() grad for finite diff vs auto diff for finite_diff vs var

./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 118.79300689265223, which exceeds tol_val, where
x1 evaluates to 36.632908591481893,
x2 evaluates to -82.160098301170336, and
tol_val evaluates to 0.0059396503446326118.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient_fvar() grad for finite_diff vs fvar<double>

[  FAILED  ] MathMixProbWienerLpdf.fiveParamWGradientExpectAd (0 ms)
[ RUN      ] MathMixProbWienerLpdf.fiveParamZeroSvWGradientExpectAd
./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 118.7930068939788, which exceeds tol_val, where
x1 evaluates to 45.335715950020479,
x2 evaluates to -73.457290943958327, and
tol_val evaluates to 0.00059396503446989407.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient() grad for finite diff vs auto diff for finite_diff vs var

./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 118.7930068939788, which exceeds tol_val, where
x1 evaluates to 45.335715950020479,
x2 evaluates to -73.457290943958327, and
tol_val evaluates to 0.0059396503446989398.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient_fvar() grad for finite_diff vs fvar<double>

[  FAILED  ] MathMixProbWienerLpdf.fiveParamZeroSvWGradientExpectAd (1 ms)
[ RUN      ] MathMixProbWienerLpdf.fullParamWGradientExpectAd
./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 70.306110772440491, which exceeds tol_val, where
x1 evaluates to 35.719186814369543,
x2 evaluates to -34.586923958070955, and
tol_val evaluates to 0.00035153055386220248.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient() grad for finite diff vs auto diff for finite_diff vs var

./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 70.306110772440491, which exceeds tol_val, where
x1 evaluates to 35.719186814369543,
x2 evaluates to -34.586923958070955, and
tol_val evaluates to 0.0035153055386220248.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient_fvar() grad for finite_diff vs fvar<double>

[  FAILED  ] MathMixProbWienerLpdf.fullParamWGradientExpectAd (5 ms)
[ RUN      ] MathMixProbWienerLpdf.existingFullRowsWGradientExpectAd
./test/unit/math/expect_near_rel.hpp:36: Failure
The difference between x1 and x2 is 70.306110772440491, which exceeds tol_val, where
x1 evaluates to 35.719186814369543,
x2 evaluates to -34.586923958070955, and
tol_val evaluates to 0.0035153055386220248.
expect_near_rel_finite in: expect_near_rel; require items x1(0) = x2(0): gradient() grad for finite diff vs auto diff for finite_diff vs var

Google Test trace:
test/unit/math/mix/prob/wiener_lpdf_ad_test.cpp:67: row_4
[  FAILED  ] MathMixProbWienerLpdf.existingFullRowsWGradientExpectAd (44 ms)
[----------] 4 tests from MathMixProbWienerLpdf (50 ms total)

[----------] Global test environment tear-down
[==========] 4 tests from 1 test suite ran. (50 ms total)
[  PASSED  ] 0 tests.
[  FAILED  ] 4 tests, listed below:
[  FAILED  ] MathMixProbWienerLpdf.fiveParamWGradientExpectAd
[  FAILED  ] MathMixProbWienerLpdf.fiveParamZeroSvWGradientExpectAd
[  FAILED  ] MathMixProbWienerLpdf.fullParamWGradientExpectAd
[  FAILED  ] MathMixProbWienerLpdf.existingFullRowsWGradientExpectAd

 4 FAILED TESTS
test/unit/math/mix/prob/wiener_lpdf_ad_test --gtest_output="xml:test/unit/math/mix/prob/wiener_lpdf_ad_test.xml" failed
exit now (05/12/26 19:51:50 CEST)
1

Which is the same I have experienced using an external implementation.

Submitted a PR.

Navigation Menu

Repository navigation

`wiener_lpdf` reverse-mode derivative wrt `w` disagrees with finite differences of Stan Math value #3322

Summary

Environment

Minimal reproducer

Observed output

Expected behavior

Possibility of a bug

Location

Related tests

Additional note

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

wiener_lpdf reverse-mode derivative wrt w disagrees with finite differences of Stan Math value #3322

Description

Summary

Environment

Minimal reproducer

Observed output

Expected behavior

Possibility of a bug

Location

Related tests

Additional note

Activity

WardBrian commented on May 12, 2026

martonaronvarga commented on May 12, 2026

SteveBronder commented on May 12, 2026

martonaronvarga commented on May 12, 2026

martonaronvarga commented on May 12, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

`wiener_lpdf` reverse-mode derivative wrt `w` disagrees with finite differences of Stan Math value #3322