Skip to content

Re-evaluate ODE default max_num_steps values #1969

Open
@bbbales2

Description

@bbbales2

Description

Practically I think the default max_num_steps arguments to both solvers are too high.

Currently we use 1e6 for the RK45 solver and 1e8 for the BDF solver.

This is the number of steps between output points in an ODE solve, and the max_num_steps is expected to throw an error if the solver is spending too much time.

The reason I think max_num_steps is too large that if we set up a simple ODE example like this:

struct CosArg1 {
  template <typename T0, typename T1, typename T2>
  inline Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1>
  operator()(const T0& t, const Eigen::Matrix<T1, Eigen::Dynamic, 1>& y,
             std::ostream* msgs, const T2& a) const {
    Eigen::Matrix<stan::return_type_t<T1, T2>, Eigen::Dynamic, 1> out(1);
    out << stan::math::cos(a * t);
    return out;
  }
};

TEST(StanMathOde_ode_bdf_tol, too_much_work) {
  Eigen::Matrix<var, Eigen::Dynamic, 1> y0
      = Eigen::VectorXd::Zero(1).template cast<var>();
  var t0 = 0;
  std::vector<var> ts = {0.45, 1e10};

  var a = 1.0;

  EXPECT_THROW_MSG(stan::math::ode_rk45(CosArg1(), y0, t0, ts, nullptr, a),
                   std::domain_error,
                   "Failed to integrate to next output time");
}

the rk45 solver takes about a second to give up. The bdf would be slower per-step, and also has a much higher default so it will be slower.

If we assume a model takes 50000 leapfrog steps to sample, and should hopefully be done within an hour, that leaves 70ms per leapfrog step. It takes roughly a second for the current RK45 defaults to fail on a simple problem like this. If we cut the default max_num_steps to 1e5, we're looking at roughly 100ms to fail, or with 1e4 that's 10ms to fail, which is probably more in line with what we'd expect a useful default to do.

People with particularly slow models can increase the max_num_steps after they've exhausted their other options (switching solvers, adjusting priors, reparameterizing the model, and changing tolerances).

Current Version:

v3.2.0

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions