# Unit test for differentiation

**URL:** <https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812>\
**Category:** Developers\
**Created:** [January 26, 2020, 7:54pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812 "2020-01-26T19:54:45Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [January 26, 2020, 7:54pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/1 "2020-01-26T19:54:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevebronder](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stevebronder/32/17242_2.png) [@stevebronder](https://discourse.mc-stan.org/u/stevebronder)\
**Post date:** [January 26, 2020, 8:03pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/2 "2020-01-26T20:03:24Z")

</div>

Mby here?

[https://mc-stan.org/math/d4/dd4/autodiff\_test\_guide.html](https://mc-stan.org/math/d4/dd4/autodiff_test_guide.html)

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [January 26, 2020, 8:10pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/3 "2020-01-26T20:10:53Z")

</div>

That looks spot on. Trying it out.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [January 26, 2020, 10:27pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/4 "2020-01-26T22:27:05Z")

</div>

Sooo… here’s my [prototype function](https://github.com/stan-dev/math/blob/try-efficient_hmm_gradient/stan/math/prim/prob/hmm_marginal_lpdf.hpp) and my [prototype unit test](https://github.com/stan-dev/math/blob/try-efficient_hmm_gradient/test/unit/math/prim/prob/hmm_marginal_test.cpp), 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);
  };

  stan::test::expect_ad(hmm_functor, 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
      ./test/unit/math/test_ad.hpp:433:35)>' requested here
  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
      'stan::test::internal::expect_all_throw<(lambda at ./test/unit/math/test_ad.hpp:433:35)>' requested here
    internal::expect_all_throw(h(0), x);
              ^
./test/unit/math/test_ad.hpp:642:13: note: in instantiation of function template specialization
      'stan::test::internal::expect_ad_helper<(lambda at test/unit/math/prim/prob/hmm_marginal_test.cpp:202:22), (lambda at
      ./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
      'stan::test::internal::expect_ad_vvv<(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
  internal::expect_ad_vvv(tols, f, x1, x2, x3);
            ^
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
  stan::test::expect_ad(tols, hmm_functor, log_omegas, Gamma, rho);

```

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.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [January 27, 2020, 7:31pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/5 "2020-01-27T19:31:14Z")

</div>

The compilation error above is fixed by replacing `value_of` with `value_of_rec`, see [this post](https://discourse.mc-stan.org/t/value-of-vs-value-of-rec-in-stan-dev-math/12807).

(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>
void expect_ad_derivatives(const ad_tolerances& tols, const G& g,
                           const Eigen::VectorXd& x, bool test_ad_hessian = 1) {
  double gx = g(x);
  test_gradient(tols, g, x, gx);
  test_gradient_fvar(tols, g, x, gx);
  if (test_ad_hessian) {
    test_hessian(tols, g, x, gx);
    test_hessian_fvar(tols, g, x, gx);
    test_grad_hessian(tols, g, x, gx);
  }
}

```

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [January 28, 2020, 1:21am UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/6 "2020-01-28T01:21:15Z")

</div>

> [@charlesm93](#):
>
> We could add a knob

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.

---

<div class="post-metadata">

**Author:** ![charlesm93](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charlesm93/32/4634_2.png) [@charlesm93](https://discourse.mc-stan.org/u/charlesm93)\
**Post date:** [January 28, 2020, 3:12pm UTC](https://discourse.mc-stan.org/t/unit-test-for-differentiation/12812/7 "2020-01-28T15:12:48Z")

</div>

Ok. The following works:

```
  stan::test::ad_tolerances tols;
  double infinity = std::numeric_limits<double>::infinity();
  tols.hessian_val_ = infinity;
  tols.hessian_grad_ = infinity;
  tols.hessian_hessian_ = infinity;
  tols.hessian_fvar_val_ = infinity;
  tols.hessian_fvar_grad_ = infinity;
  tols.hessian_fvar_hessian_ = infinity;
  tols.grad_hessian_val_ = infinity;
  tols.grad_hessian_hessian_ = infinity;
  tols.grad_hessian_grad_hessian_ = infinity;

  stan::test::expect_ad(tols, hmm_functor, log_omegas, Gamma, rho);

```