# Unit test for differentiation

I recalled @Bob_Carpenter worked on ways to automate the testing of our autodiff against finite diff. After a quick search, I couldn’t find any doc on on discourse, the Stan wiki, or the test functions I looked at (for normal and multinormal). So I’m posting this here.

I’m familiar with the old fashion way, which means coding the call to the `grad()` function and manually coding finite diff. Nothing too bad, but if there’s a slicker solution, I’d prefer using it.

1 Like
3 Likes

That looks spot on. Trying it out.

Sooo… here’s my prototype function and my prototype unit test, with the relevant piece of code:

``````  auto hmm_functor = [](const auto& log_omegas,
const auto& Gamma,
const auto& rho) {
return hmm_marginal_lpdf(log_omegas, Gamma, rho);
};

``````

But this doesn’t compile with a long error message that repeats itself three times (I’m guessing once for each argument):

``````lib/eigen_3.3.3/Eigen/src/Core/functors/AssignmentFunctors.h:24:104: error: assigning to 'double' from incompatible type
'const stan::math::fvar<stan::math::var>'
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(DstScalar& a, const SrcScalar& b) const { a = b; }
^
lib/eigen_3.3.3/Eigen/src/Core/AssignEvaluator.h:637:15: note: in instantiation of member function 'Eigen::internal::assign_op<double,
stan::math::fvar<stan::math::var> >::assignCoeff' requested here
m_functor.assignCoeff(m_dst.coeffRef(index), m_src.coeff(index));
^
lib/eigen_3.3.3/Eigen/src/Core/AssignEvaluator.h:497:14: note: in instantiation of member function
'Eigen::internal::generic_dense_assignment_kernel<Eigen::internal::evaluator<Eigen::Matrix<double, -1, -1, 0, -1, -1> >,
Eigen::internal::evaluator<Eigen::Matrix<stan::math::fvar<stan::math::var>, -1, -1, 0, -1, -1> >,
Eigen::internal::assign_op<double, stan::math::fvar<stan::math::var> >, 0>::assignCoeff' requested here
kernel.assignCoeff(i);
^
lib/eigen_3.3.3/Eigen/src/Core/AssignEvaluator.h:741:34: note: in instantiation of member function
'Eigen::internal::dense_assignment_loop<Eigen::internal::generic_dense_assignment_kernel<Eigen::internal::evaluator<Eigen::Matrix<double,
-1, -1, 0, -1, -1> >, Eigen::internal::evaluator<Eigen::Matrix<stan::math::fvar<stan::math::var>, -1, -1, 0, -1, -1> >,
Eigen::internal::assign_op<double, stan::math::fvar<stan::math::var> >, 0>, 1, 0>::run' requested here
dense_assignment_loop<Kernel>::run(kernel);
^
lib/eigen_3.3.3/Eigen/src/Core/AssignEvaluator.h:879:5: note: in instantiation of function template specialization
'Eigen::internal::call_dense_assignment_loop<Eigen::Matrix<double, -1, -1, 0, -1, -1>,
Eigen::Matrix<stan::math::fvar<stan::math::var>, -1, -1, 0, -1, -1>, Eigen::internal::assign_op<double,
stan::math::fvar<stan::math::var> > >' requested here
call_dense_assignment_loop(dst, src, func);
^
lib/eigen_3.3.3/Eigen/src/Core/AssignEvaluator.h:836:46: note: in instantiation of member function
'Eigen::internal::Assignment<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<stan::math::fvar<stan::math::var>, -1, -1,
0, -1, -1>, Eigen::internal::assign_op<double, stan::math::fvar<stan::math::var> >, Eigen::internal::Dense2Dense, void>::run'
requested here
Assignment<ActualDstTypeCleaned,Src,Func>::run(actualDst, src, func);
^
lib/eigen_3.3.3/Eigen/src/Core/PlainObjectBase.h:728:17: note: (skipping 5 contexts in backtrace; use -ftemplate-backtrace-limit=0 to
see all)
internal::call_assignment_no_alias(this->derived(), other.derived(), internal::assign_op<Scalar,typename OtherDerived::S...
^
./test/unit/math/test_ad.hpp:354:3: note: in instantiation of function template specialization
'stan::test::internal::expect_throw<stan::math::fvar<stan::math::fvar<stan::math::var> >, (lambda at
expect_throw<fvar<fvar<var>>>(f, x, "fvar<fvar<var>>");
^
./test/unit/math/test_ad.hpp:442:15: note: in instantiation of function template specialization
internal::expect_all_throw(h(0), x);
^
./test/unit/math/test_ad.hpp:642:13: note: in instantiation of function template specialization
./test/unit/math/test_ad.hpp:637:13), Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>,
Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
internal::expect_ad_helper(tols, f, g1, serialize_args(x1), x1, x2, x3);
^
./test/unit/math/test_ad.hpp:1173:13: note: in instantiation of function template specialization
-1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
^
test/unit/math/prim/prob/hmm_marginal_test.cpp:210:15: note: in instantiation of function template specialization
'stan::test::expect_ad<(lambda at test/unit/math/prim/prob/hmm_marginal_test.cpp:202:22), Eigen::Matrix<double, -1, -1, 0, -1,
-1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >' requested here
``````

The function I wrote uses `operands_and_partials`, which if my understanding is correct works for rev mode, but `expect_ad` also checks the forward mode. Questions: (i) should I expect the function to be forward differentiable and (ii) if not, is there a routine I can call to only test reverse mode (but the looks of `test_ad.hpp`, yes.

The compilation error above is fixed by replacing `value_of` with `value_of_rec`, see this post.

(Code compiles and runs, though the unit tests seem to fail for the hessians.)

EDIT: The code doesn’t compute Hessians, because all the adjoint calculations are done using Eigen. Stan functions usually incur an overhead that slows down computation, and since the primary use of the function is gradient-based algorithms, I’m not keen to get a slow down.

We could add a knob to only test first-order derivatives. Something like that (in `test_ad.hpp`, specifically `expect_ad_derivatives`):

``````template <typename G>
const Eigen::VectorXd& x, bool test_ad_hessian = 1) {
double gx = g(x);
test_hessian(tols, g, x, gx);
test_hessian_fvar(tols, g, x, gx);
}
}
``````

There are lots of tests and lots of entry points. The easiest thing to do would be to generalize `ad_tolerances` to allow individual tests to be turned off the same way their tolerances are set now. In fact, setting a tolernace of +infinity could signal to not do a test.

1 Like

Ok. The following works:

``````  stan::test::ad_tolerances tols;
double infinity = std::numeric_limits<double>::infinity();
tols.hessian_val_ = infinity;