8
$\begingroup$

I'm attempting to find a precise as possible approximation to the logarithm of the Modified Bessel Function of the Second Kind:

$$\log K_{\alpha}(x) = \log\Big[\frac{1}{2}\int_0^{∞} t^{\alpha-1} \exp\{-\frac{1}{2}x(t+t^{-1})\}dt\Big].$$

The problem is that $K_{\alpha}(x)$ diverges both as $\alpha\rightarrow\infty$ and as $x\rightarrow0$, which is why I'm taking the logarithm in the first place. But of course the log trick only works if I can use it to break up the expression in some way, which I can't do in this case because I can't pass it through the integral.

There is an asymptotic approximation for large $\alpha$, given by:

$$K_{\alpha}(x)\sim\sqrt{\frac{\pi}{2\alpha}}\Big(\frac{ex}{2\alpha}\Big)^{-\alpha}.$$

Clearly this approximation can be broken up by the log, unfortunately it is not sufficiently precise (for my purposes) for values of $\alpha$ and $x$ which cause overflow.

Does anyone have another way to get a more precise estimate for $\log K_{\alpha}(x)$?

If you want you can assume that $\alpha=k-\frac{1}{2}$ for $k\in{\bf Z}_+$.

$\endgroup$
1
  • $\begingroup$ the integral will be dominated by the stationary point of the integrand at $x\approx 2 (\alpha+1)$. For large enough $\alpha$,$\exp(-1/2x)$ is almost constant in this region an might be taylor expanded. What remains are $\Gamma$ type integrals which can be tackled by stirlings approximation. Can you take it from here ? wolframalpha.com/input/…*Exp%5B-t%2F2%5D%7D,%7Bt,0,40%7D%5D $\endgroup$
    – tired
    Oct 10, 2016 at 23:57

1 Answer 1

6
$\begingroup$

I am currently struggling with the same problem, getting a bit further but not all the way. If you solved it in the meantime, I'd be glad to hear about it.

EDIT 2: Improved answer some more as my technique improved some more.

EDIT 3: The pull request (linked below) now has somewhat better implementation then discussed here that reduces the error a bit further, but still breaks for some values.

C++ code (written for Stan math library) can be found at https://github.com/stan-dev/math/pull/1121

Case 0: $x$ small relative to $\alpha$

I use this when $\alpha > 15$ AND $10x < \sqrt{\alpha + 1}$. The formula can be found on Wikipedia on Bessel K and also as formula 1.10 of Temme, Journal of Computational Physics, vol 19, 324 (1975). It has:

$$ K_\alpha(x) \simeq \frac{\Gamma(\alpha)}{2} \left( \frac{2}{x} \right)^\alpha $$

Case 1: Small to medium $\alpha$ and $x$

My main workhorse is using the logarithm of Equation 26 of Rothwell: Computation of the logarithm of Bessel functions of complex argument and fractional order.

which has (maintaining their notation)

$$ K_\nu(z) = \frac{\sqrt{\pi}}{\Gamma(\nu + \frac{1}{2})}(2z)^{-\nu}e^{-z} \int_0^1 [\beta e^{-u^\beta}(2z + u^\beta)^{\nu-1/2} u^{n-1} + e^{-1/u} u^{-2\nu-1}(2zu + 1)^{\nu-1/2}] \mathrm{d}u $$

Where $\beta = \frac{2n}{2\nu+1}$ and the authors suggest using $n = 8$.

You still can't push the log through the integral, but for large areas of the parameter space, it is OK to just compute the integral and take the log afterwards.

The integral ceases to be numerically tractable around roughly $\nu > 50$ or when $\nu > 0.5$ AND $log(z) > \frac{300}{\nu - 0.5} - log(2)$ (one of the summands turn to infinity).

Fruther, the formula ceases to be very accurate for small $z$ and small $\nu$, where below roughly $z < 10^-4$ the relative error starts climbing from around $10^{-30}$ to around $10{^-2}$ for $z < 10^{-7}$

There is a C++ implementation of the Rothwell formula at https://github.com/stan-dev/stan/wiki/Stan-Development-Meeting-Agenda/0ca4e1be9f7fc800658bfbd97331e800a4f50011

Case 2: $\alpha \gg x$

Here the log of the approximate formula you've written works OK, but I got a bit better results with the same formula as in Case 0.

Case 3: $x \gg \alpha$

Asymptotic formula such as 10.4.2 here: https://dlmf.nist.gov/10.40 can be used. Once again, you need to compute the sum on non-log scale (it has negative elements), and take the log afterward, but this has worked great for me. Summing up to 10 first elements worked OK in my case.

Remaining cases For large $x$ with comparably large $\alpha$, I still haven't found a reliable solution.

Here is a plot of error with respect to recursive formula for consecutive neighbours with the approach I've described (ratio is the relative error of the actual value, I also test the gradient of this function wrt. both parameters as computed by Stan's autodiff - the gradient is a tougher problem, but probably not relevant here):

plot of relative error of the approximation

The biggest relative error shown is 3.8e-02 for $\alpha = 148$ and $x = 105$.

$\endgroup$
3
  • 1
    $\begingroup$ When you say 'compute the integral' do you mean numerically compute it with quadrature? ..it seems that would be relatively very slow compared to just evaluating it with the built-in functions in whatever computing environment you're using? $\endgroup$
    – sonicboom
    Jun 24, 2019 at 18:56
  • 1
    $\begingroup$ I used the built-in functions for numerical quadrature. I don't think you can solve the integral analytically. My actual code (within Stan math library) can be found at github.com/stan-dev/math/pull/1121 and is yours to inspect/use. $\endgroup$ Jun 25, 2019 at 12:38
  • $\begingroup$ Ok cheers, I'm interested in integer orders for the moment but at a later stage I'll be moving to fractional orders and I think I will run into this issue so its good to see some work on it. $\endgroup$
    – sonicboom
    Jun 26, 2019 at 14:09

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .