Inefficient modified bessel function and von_mises_lpdf derivative calculation

[venpopov]

Description

The current implementation of the von_mises_lpdf can be sped up ~4 times with a minor change to the derivative calculation that doesn't change the numerical output in any way.

The relevant code is:

edge<2>(ops_partials).partials_
        = cos_mu_minus_y
          - modified_bessel_first_kind(-1, kappa_val)
                / modified_bessel_first_kind(0, kappa_val);

For integer order, $I_n(x) = I_{-n}(x)$ (see 10.27.1: https://dlmf.nist.gov/10.27). Thus modified_bessel_first_kind(-1, kappa_val) is the same as modified_bessel_first_kind(1, kappa_val).

modified_bessel_first_kind() calls boost::math::cyl_bessel_i(), which has a special optimizations for order 0 and 1 (see. When cyl_bessel_i() is called with order -1, it uses a much more inefficient calculation.

Calculating the modified_bessel_function is the most expensive part of calculating the von_mises_lpdf, taking ~80% of the time.

Example

The following speed test shows that both functions are equal, but modified_bessel_first_kind(1,x) is ~10 times faster:

library(rstan)
library(microbenchmark)

x <- rgamma(1000,10,0.5)
stanFunction('modified_bessel_first_kind',nu=0, x=x)
microbenchmark(modified_bessel_first_kind(1,x),
               modified_bessel_first_kind(-1,x), check = 'equal')

Unit: microseconds
                              expr      min        lq       mean   median       uq      max neval
  modified_bessel_first_kind(1, x)   95.118   96.0355   96.87414   96.406   97.112  112.300   100
 modified_bessel_first_kind(-1, x) 1240.216 1246.2575 1251.92232 1251.828 1253.942 1318.924   100

Expected Output

Change modified_bessel_first_kind(-1, kappa_val)tomodified_bessel_first_kind(1, kappa_val) here

Current Version:

v4.8.0

[andrjohns]

Fantastic catch, thanks! I'll put together a PR for this now so it will be in the 2.35 release

[venpopov]

Awesome, thanks! I was just about to try to figure out how to do that myself, but never having contributed before, it was gonna take some time :)

[andrjohns]

No worries! Did you want me to leave it for you still? Or just tag you once the PR is up so you can see what would be needed for next time?

[andrjohns]

And to clarify, what would be needed for a PR would be to:

• Change the modified_bessel_first_kind call in von_mises_lpdf that you mentioned
• Update the von_mises mix test to use our expect_ad gradient-testing framework
  • For an example of this, see the gamma_lcdf test

[venpopov]

Yes, please do it and tag me! I opened a second issue for a more general application for other models that use the modified bessel function of order 0 (which is what actually led me to discover this speed up), so I can perhaps tackle that after seeing how you do it here and reading the docs. Although in the new issue I guess first the devs have to agree on the corrects implementation.