Integrate_ode_bdf vs integrate_ode_rk45 for integrand that doesn't reference state variable

Developers
1

This is for stan-math used from C++; I don’t have the entire stan.

Also, a general thank-you for this very useful AD software.

I noted a difference between the behavior of the cvode bdf and boost rk45 integrator compilation when I switched to the bdf integrator for the first time today.

If I am just integrating a function, where the integrand refers only to the independent variable and the parameters, rather than the current value of the antiderivative, then there is a split between the behavior of these two functions that I wasn’t expecting.

A minimal example follows, with the differences between the two versions called out in the comments:

#include "stan/math/mix/mat.hpp"
#include <iostream>

int main(int,char*[]){
	stan::math::start_nested();
	double t(0.);
	std::vector<double> ts{2.};
	std::vector<stan::math::var> theta{1.4};
	std::vector<double> A0{0.};
	std::vector<double> xr;
	std::vector<int> xi;
    auto a([](auto t,auto&y[[maybe_unused]],auto&th,auto&,auto&,auto*){
 	 return std::vector<decltype(th[0]*t*y[0])>{th[0]*t};});
    // A= 2*1.4 = 2.8; next line will fail on line 94 of ode_system.hpp if *y[0] is
    // commented out of decltype in the lambda function above due to conversion of
    // vector<double> to vector<var>, but won't fail if using integrate_ode_rk45
	auto A(stan::math::integrate_ode_bdf(a,A0,t,ts,theta,xr,xi));
    A[0][0].grad();
    std::cout<<A[0][0]<<", "<<theta[0].adj()<<"\n";//should print: 2.8, 2
	stan::math::recover_memory_nested();
	return 0;
}

Line 94 of ode_system.hpp is in the 2nd Jacobian routine, the one that calculates the Jacobian with respect to both theta and y:
vector<var> dy_dt_var = f_(t, y_var, theta_, x_, x_int_, msgs_);

I have no idea if it makes sense to make bdf and rk45 have the same behavior for this use case. It is strange that rk45 is able to compute the same sensitivities even with *y[0] commented out.

For those who might come after, I found that on Windows MinGW GCC, before the cvode version would link, I had to go to stan-dev/math and issue the following command (discovered by reading the makefile therein):

make CC=g++ lib/cvodes_2.9.0/lib/libsundials_cvodes.a lib/cvodes_2.9.0/lib/libsundials_nvecserial.a

(and of course add these to LDLIBS in my own project makefile)

1 Like
2

This is tricky, but I think I got what is happening. So if Iget you right, then you run into problems whenever you make this change:

works with integrate_ode_bdf:

  auto a([](auto t,auto& y[[maybe_unused]],auto&th,auto&,auto&,auto*) {
      return std::vector<decltype(th[0]*t*y[0])>{th[0]*t};}
    );

fails with integrate_ode_bdf:

  auto a([](auto t,auto& y[[maybe_unused]],auto&th,auto&,auto&,auto*) {
      return std::vector<decltype(th[0]*t)>{th[0]*t};}
    );

Correct?

First: The version which fails is simply wrong in how it determines the return type. The functor must return std::vector<var> if y0 or theta is a var. The second failing version does not do this and hence is plain wrong in this regard.

You may wonder why rk45 works whereas the bdf version does not. The thing is that the rk45 version will always calculate the derivative wrt to y0 and theta. For this case your functor correctly returns a var. Now, the bdf integrator will also need in addition the Jacobian wrt to the states only. So the bdf integrator will call the functor with y0 being var while theta is just a double. The code expects a var as return in this case which it won’t get.

That’s why the two codes are different in that respect. So the wrong return type bug of the user functor is trickled when you calculate the Jacobian wrt to states only which is what the stiff solver needs.

Makes sense?

1 Like
3

That makes perfect sense, and I appreciate you taking the time to explain it to me.

1 Like