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

Mby here?

https://mc-stan.org/math/d4/dd4/autodiff_test_guide.html

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);
  };

  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.

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>
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);
  }
}

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;
  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);
1 Like