From 8ac653aab5f793a7f4a652c1abb81e98503de944 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 3 Sep 2026 13:10:42 -0400 Subject: [PATCH 1/2] test(jacobian): add failing tests for non-square fwd-mode jacobian Regression tests for #3385. Both fail against develop, aborting inside jacobian() before reaching any assertion: Assertion failed: (rows == this->rows() && cols == this->cols() && "DenseBase::resize() does not actually allow to resize.") The fwd jacobian sizes J as n x n and fx as n from the input before evaluating f, so J.col(0) = fx_fvar.d() assigns an m-vector into an n-row column whenever m != n. Existing coverage in test/unit/math/mix/functor/autodiff_test.cpp only uses an R^2 -> R^2 functor, where m == n masks the mis-sizing entirely. These cover both directions: R^2 -> R^3 (more outputs than inputs, the case forward mode is actually cheaper for) and R^3 -> R^2, each checked against a hand-computed Jacobian. The tests live in fwd/ and include only stan/math/fwd.hpp, which is also the realistic reproduction path: with rev in scope the rev overload wins partial ordering and handles these functors correctly. --- test/unit/math/fwd/functor/jacobian_test.cpp | 75 ++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 test/unit/math/fwd/functor/jacobian_test.cpp diff --git a/test/unit/math/fwd/functor/jacobian_test.cpp b/test/unit/math/fwd/functor/jacobian_test.cpp new file mode 100644 index 00000000000..f11292ddc1e --- /dev/null +++ b/test/unit/math/fwd/functor/jacobian_test.cpp @@ -0,0 +1,75 @@ +#include +#include +#include + +using Eigen::Dynamic; +using Eigen::Matrix; + +// fun_2_3: R^2 --> R^3 | (x, y) --> [x, (x + y), (x * y)] +struct fun_2_3 { + template + inline Matrix operator()( + const Matrix& x) const { + Matrix z(3); + z << x(0), x(0) + x(1), x(0) * x(1); + return z; + } +}; + +// fun_3_2: R^3 --> R^2 | (x, y, z) --> [(x * y), (y + 2 * z)] +struct fun_3_2 { + template + inline Matrix operator()( + const Matrix& x) const { + Matrix z(2); + z << x(0) * x(1), x(1) + 2.0 * x(2); + return z; + } +}; + +TEST(FwdFunctor, jacobianMoreOutputsThanInputs) { + fun_2_3 f; + Matrix x(2); + x << 1.5, 2.0; + + Matrix fx; + Matrix J; + stan::math::jacobian(f, x, fx, J); + + EXPECT_EQ(3, fx.size()); + EXPECT_FLOAT_EQ(x(0), fx(0)); + EXPECT_FLOAT_EQ(x(0) + x(1), fx(1)); + EXPECT_FLOAT_EQ(x(0) * x(1), fx(2)); + + EXPECT_EQ(3, J.rows()); + EXPECT_EQ(2, J.cols()); + EXPECT_FLOAT_EQ(1, J(0, 0)); + EXPECT_FLOAT_EQ(0, J(0, 1)); + EXPECT_FLOAT_EQ(1, J(1, 0)); + EXPECT_FLOAT_EQ(1, J(1, 1)); + EXPECT_FLOAT_EQ(x(1), J(2, 0)); + EXPECT_FLOAT_EQ(x(0), J(2, 1)); +} + +TEST(FwdFunctor, jacobianFewerOutputsThanInputs) { + fun_3_2 f; + Matrix x(3); + x << 1.5, 2.0, -3.0; + + Matrix fx; + Matrix J; + stan::math::jacobian(f, x, fx, J); + + EXPECT_EQ(2, fx.size()); + EXPECT_FLOAT_EQ(x(0) * x(1), fx(0)); + EXPECT_FLOAT_EQ(x(1) + 2.0 * x(2), fx(1)); + + EXPECT_EQ(2, J.rows()); + EXPECT_EQ(3, J.cols()); + EXPECT_FLOAT_EQ(x(1), J(0, 0)); + EXPECT_FLOAT_EQ(x(0), J(0, 1)); + EXPECT_FLOAT_EQ(0, J(0, 2)); + EXPECT_FLOAT_EQ(0, J(1, 0)); + EXPECT_FLOAT_EQ(1, J(1, 1)); + EXPECT_FLOAT_EQ(2, J(1, 2)); +} From 9f61913945eeff9218ac5ce7aa87e6c0999a61a3 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Thu, 3 Sep 2026 13:10:42 -0400 Subject: [PATCH 2/2] fix(jacobian): size fwd-mode outputs from the range of f The fwd jacobian sized both outputs from the domain before f was ever applied, so the range dimension m was simply assumed equal to n: J.resize(x_fvar.size(), x.size()); fx.resize(x_fvar.size()); ... Matrix, Dynamic, 1> fx_fvar = f(x_fvar); J.col(0) = fx_fvar.d(); For m != n that assigns an m-vector into an n-row column. Move the J.resize below the call to f and take its rows from fx_fvar.size(), which mirrors what the rev implementation already does. The fx.resize is dropped rather than moved: fx = fx_fvar.val() is a full-matrix assignment that resizes on its own, which is why fx.size() came out correct even in the aborting runs. Moving it would only have restated that. Failure modes before the fix, both directions, with a functor that is not square: assertions on (math's own build never defines NDEBUG) - abort. -DNDEBUG, m > n - Eigen's assignment loop bounds on the destination block, so the copy truncates: J stays n x n and holds a partial Jacobian. Silent wrong answer, no out-of-bounds access. -DNDEBUG, m < n - the same destination bound reads past the source temporary. ASan reports a heap-buffer-overflow READ at jacobian.hpp:25, 8 bytes past fx_fvar. Note that this narrows the issue's description: there is no write past J's storage in either direction, so no heap corruption. No signature change, so overload resolution is untouched and the eight in-tree callers (all in rev/functor, all resolving to the rev overload, and all square) are unaffected. Verified with ASan clean in both directions, agreement with the rev overload on the same functor, the higher-order T = fvar instantiation, the fwd/functor suite, mix/functor/autodiff_test, and rev/functor/algebra_solver_fp_test. Closes #3385 --- stan/math/fwd/functor/jacobian.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/fwd/functor/jacobian.hpp b/stan/math/fwd/functor/jacobian.hpp index 030e16c1749..2bbb0788c77 100644 --- a/stan/math/fwd/functor/jacobian.hpp +++ b/stan/math/fwd/functor/jacobian.hpp @@ -14,14 +14,14 @@ inline void jacobian(const F& f, const Eigen::Matrix& x, using Eigen::Dynamic; using Eigen::Matrix; Matrix, Dynamic, 1> x_fvar(x.size()); - J.resize(x_fvar.size(), x.size()); - fx.resize(x_fvar.size()); for (int k = 0; k < x.size(); ++k) { x_fvar(k) = fvar(x(k), 0); } x_fvar(0) = fvar(x(0), 1); Matrix, Dynamic, 1> fx_fvar = f(x_fvar); + // size the outputs from the range of f, which is only known once f is applied fx = fx_fvar.val(); + J.resize(fx_fvar.size(), x.size()); J.col(0) = fx_fvar.d(); const fvar switch_fvar(0, 1); // flips the tangents on and off for (int i = 1; i < x.size(); ++i) {