3 ms·

Turing.jl is in an interesting spot because it is essentially a DSL-free probabilistic programming language. While it technically has a DSL of sorts given by th

Turing.jl is in an interesting spot because it is essentially a DSL-free probabilistic programming language. While it technically has a DSL of sorts given by the `@model` macro, anything that is AD-compatible can be used in this macro and since Julia's AD tools work on things written in the Julia language, this means that you can just throw code from other Julia packages into Turing and just expect AD-compatible things to work with Hamiltonian Monte Carlo and all of that. So things like DifferentialEquations.jl ODEs/SDEs/DAEs/DDEs/etc. work quite well with this, along with other "weird things for a probabilistic programming language to support" like nonlinear solving (via NLsolve.jl) or optimization (via Optim.jl, and I mean doing Bayesian inference where a value is defined as the result of an optimization). If you are using derivative-free inference methods, like particle sampling methods or variants of Metropolis-Hastings, then you can throw pretty much any existing Julia you had as a nonlinear function and do inference around it. So while it's in some sense similar to PyMC3 or Stan, there's a huge difference in the effective functionality that you get by supporting a language-wide infrastructure vs the more traditional method of one-by-one adding features and documenting them. So while PyMC3 ran a Google Summer of Code to get some ODE support (https://docs.pymc.io/notebooks/ODE_API_introduction.html https://docs.pymc.io/notebooks/ODE_API_introduction.html) and Stan has 2 built-in methods you're allowed to use (https://mc-stan.org/docs/2_19/stan-users-guide/ode-solver-chapter.html https://mc-stan.org/docs/2_19/stan-users-guide/ode-solver-ch...), with Julia you get all of DifferentialEquations.jl just because it exists (https://docs.sciml.ai/latest/ https://docs.sciml.ai/latest/). This means that Turing.jl doesn't document or doesn't have to document most of its features, but they exist due to composibility. That's quite different from a "top down" approach to library support. This explains why Turing has been able to develop so fast as well, since it's developer community isn't just "the people who work on Turing", but it's pretty much the whole ecosystem of Julia. Its distributions are defined by Distributions.jl (https://github.com/JuliaStats/Distributions.jl https://github.com/JuliaStats/Distributions.jl), its parallelism is given by Julia's base parallelism work + everything around it like CuArrays.jl and KernelAbstractions.jl (https://github.com/JuliaGPU/KernelAbstractions.jl https://github.com/JuliaGPU/KernelAbstractions.jl), derivatives come from 4 libraries, ODEs from etc. the list keeps going. So bringing it back to deep learning, Turing currently has 4 modes for automatic differentiation (https://turing.ml/dev/docs/using-turing/autodiff https://turing.ml/dev/docs/using-turing/autodiff), and thus supports any library that's compatible with those. It turns out that Flux.jl is compatible with them, so therefore Turing.jl can do Bayesian deep learning. In that sense it's like Edward or Pyro, but supporting "anything that AD's with Julia AD packages" (which soon will allow multi-AD overloads via ChainRules.jl) instead of "anything on TensorFlow graphs" or "anything compatible with PyTorch". As for performance and robustness, I mentioned in a SciML ecosystem release today that our benchmarks pretty clearly show Turing.jl as being more robust than Stan while achieving about a 3x-5x speedup in ODE parameter estimation (https://sciml.ai/2020/05/09/ModelDiscovery.html https://sciml.ai/2020/05/09/ModelDiscovery.html). However, that's utilizing the fact that Turing.jl's composibility with packages gives it top notch support (I want to work with Stan developers so we can use our differential equation library with their samplers to better isolate differences and hopefully improve both PPLs, but for now we have what we have). If you isolate it down to just "Turing.jl itself", it has wins and losses against Stan (https://github.com/TuringLang/Turing.jl/wiki https://github.com/TuringLang/Turing.jl/wiki). That said, there's some benchmarks which indicate using the ReverseDiff AD backend will give about 2 orders of magnitude performance increases in many situations (https://github.com/TuringLang/Turing.jl/issues/1140 https://github.com/TuringLang/Turing.jl/issues/1140, note that ThArrays is benchmarking PyTorch AD here) which would then probably tip the scales in Turing's favor. As for benchmarking against Pyro or Edward, it would probably just come down to benchmarking the AD implementations.

sgt101 6y ago
Thank you for this very helpful comment.
cuban_MMA 6y ago
Hi Chris, as a heads up, Stan actually has had three built-in methods for a while now. There is a non-stiff Adams-Moulton solver introduced in 2018. It unfortunately was only just exposed in the Stan 2.23 documentation: https://mc-stan.org/docs/2_23/stan-users-guide/ode-solver-chapter.html https://mc-stan.org/docs/2_23/stan-users-guide/ode-solver-ch.... Certainly, the devs have been talking about adding more solvers for a while, including SDE and DDE solvers, and your DifferentialEquations.jl ecosystem is an excellent model; it is an area that we know Stan has been lacking in. I think Steve Gronder will be trying to work with you regarding benchmarking.
- ChrisRackauckas 6y ago
  Awesome. If we can get FFI with Stan I'd like to connect DifferentialEquations.jl to it and poke at it with various problems to see how well it does on a few things. We can provide custom gradients if there's an interface for it, but I couldn't figure out how to do it without modifying the Stan source itself.
  - ivanyashchuk 6y ago
    In order to connect Stan with DifferentialEquations.jl the steps would be: 1. Create "diffeqcpp", a C++ interface to DifferentialEquations.jl (that would be similar to diffeqpy, diffeqr) possibly using CxxWrap.jl 2. Make it possible to evaluate vector-Jacobian products (VJP) with "diffeqcpp". Probably that would require ODE RHS to be coded as a string of Julia code, to make Julia AD libraries compatible with it. At this point, it should be possible to call Julia solvers from C++ and evaluate the derivatives. In Stan, there is stan::math::adj_jac_apply that makes it possible to define custom functions with custom VJP without having to deal with Stan autodiff types, it works for example with Eigen::Matrix<double>. https://discourse.mc-stan.org/t/adj-jac-apply/5163 https://discourse.mc-stan.org/t/adj-jac-apply/5163 3. Make a class (let's call it JuliaODESolver) that implements two methods: operator() // calls Julia solver for the given input multiply_adjoint_jacobian() // evaluates VJP for the given vector 4. In .stan file add a custom function in "functions {}" block, and write a header file that implements that custom function. That would probably be one line return stan::math::adj_jac_apply<JuliaODESolver>(ode_solver_inputs); More info on using external C++ code is in Section 4.5 CmdStan Manual. 5. Modify cmdstan/main.cpp to initialize and finalize Julia context to be able to call Julia functions. This is probably the only place where Stan source itself needs to be modified. I don't know what would be needed to make forward mode, and higher-order derivatives to work. I think it would be much better for a fair benchmarking if there was a convenient and documented interface to Stan algorithms to use with user-provided log-density function, similar to DynamicHMC.jl and AdvancedHMC.jl libraries. It would be then easy to call it from Julia/Python/R/C++ or anything else.

sgt101 6y ago

Thank you for this very helpful comment.

cuban_MMA 6y ago

Hi Chris, as a heads up, Stan actually has had three built-in methods for a while now. There is a non-stiff Adams-Moulton solver introduced in 2018. It unfortunately was only just exposed in the Stan 2.23 documentation: https://mc-stan.org/docs/2_23/stan-users-guide/ode-solver-chapter.html https://mc-stan.org/docs/2_23/stan-users-guide/ode-solver-ch.... Certainly, the devs have been talking about adding more solvers for a while, including SDE and DDE solvers, and your DifferentialEquations.jl ecosystem is an excellent model; it is an area that we know Stan has been lacking in. I think Steve Gronder will be trying to work with you regarding benchmarking.

ChrisRackauckas 6y ago
Awesome. If we can get FFI with Stan I'd like to connect DifferentialEquations.jl to it and poke at it with various problems to see how well it does on a few things. We can provide custom gradients if there's an interface for it, but I couldn't figure out how to do it without modifying the Stan source itself.
- ivanyashchuk 6y ago
  In order to connect Stan with DifferentialEquations.jl the steps would be: 1. Create "diffeqcpp", a C++ interface to DifferentialEquations.jl (that would be similar to diffeqpy, diffeqr) possibly using CxxWrap.jl 2. Make it possible to evaluate vector-Jacobian products (VJP) with "diffeqcpp". Probably that would require ODE RHS to be coded as a string of Julia code, to make Julia AD libraries compatible with it. At this point, it should be possible to call Julia solvers from C++ and evaluate the derivatives. In Stan, there is stan::math::adj_jac_apply that makes it possible to define custom functions with custom VJP without having to deal with Stan autodiff types, it works for example with Eigen::Matrix<double>. https://discourse.mc-stan.org/t/adj-jac-apply/5163 https://discourse.mc-stan.org/t/adj-jac-apply/5163 3. Make a class (let's call it JuliaODESolver) that implements two methods: operator() // calls Julia solver for the given input multiply_adjoint_jacobian() // evaluates VJP for the given vector 4. In .stan file add a custom function in "functions {}" block, and write a header file that implements that custom function. That would probably be one line return stan::math::adj_jac_apply<JuliaODESolver>(ode_solver_inputs); More info on using external C++ code is in Section 4.5 CmdStan Manual. 5. Modify cmdstan/main.cpp to initialize and finalize Julia context to be able to call Julia functions. This is probably the only place where Stan source itself needs to be modified. I don't know what would be needed to make forward mode, and higher-order derivatives to work. I think it would be much better for a fair benchmarking if there was a convenient and documented interface to Stan algorithms to use with user-provided log-density function, similar to DynamicHMC.jl and AdvancedHMC.jl libraries. It would be then easy to call it from Julia/Python/R/C++ or anything else.