In variational inference on von Mises clusters, how to find a bound for the log-Bessel function? This paper on von Mises clustering uses an upper bound on the modified log-Bessel function that I struggle to replicate. Taking results from this paper, the authors state:
$$u\frac{I'_\nu(u)}{I_\nu(u)} \le \sqrt{u^2 + v^2} $$ $$  \frac{I'_\nu(u)}{I_\nu(u)} \le \frac{\sqrt{u^2 + v^2}}{u} = \sqrt{1 + \frac{v^2}{u^2}} \tag{1}$$
The paper has a few apparent typos, but it seems to intend this bound derived by "integrating over $u > 0$," which I think means the indefinite integral:
$$\int\frac{\sqrt{u^2 + \nu^2}}{u} du = \sqrt{u^2 + \nu^2} - \nu \log(\nu (\sqrt{u^2 + \nu^2} + \nu)) + \nu \log(u) + c \tag{2}$$
The definite integral over $u >0 $ doesn't converge, but the proposed bound seems to resemble the indefinite integral of the right-hand side, less a constant $c$. Though, when I plot the indefinite integral for $c=0$, against the modified log-Bessel, the bound only holds through a few orders of $\nu$. (Image and code below.)
How might one approach integrating $(1)$ to find a bound that holds? I thought of integrating over $(a,x)$ and setting $a$ such that $I_\nu
(a) = 1$. Then,
$$\int_a^x\frac{I'_\nu(u)}{I_\nu(u)}du = \log I_\nu(u)\rvert_a^x \le \int_a^x\frac{\sqrt{u^2 + v^2}}{u} du$$
$$\log I_\nu(x) - \log I_\nu(a) \le \sqrt{u^2 + \nu^2} - \nu \log(\nu (\sqrt{u^2 + \nu^2} + \nu)) + \nu \log(u)\rvert_a^x$$
Which given that $\log I_\nu(a) = 0$ and some right-hand finagling gives:
$$\log I_\nu(x) \le \sqrt{x^2 + \nu^2} - \sqrt{a^2 + \nu^2} - \nu \log \frac{a(\sqrt{x^2 - \nu^2} + \nu)}{x(\sqrt{a^2 - \nu^2} + \nu)}\tag{3}$$
This bound passes a quick visual check for $\nu > \frac{1}{2}$, and for $\nu = \frac{1}{2}$ seems to be equal within machine precision. Is there a way to derive the result using the indefinite integral as seemingly employed in the paper, or is my bound using the definite integral correct?


"""
Created on Tue Jan 10 09:02:19 2017

@author: seaneaster
"""

import numpy as np
from matplotlib import pylab as plt
from scipy.special import iv

def upper_bound(u,v):
    sum_squares = np.power(u,2) + np.power(v,2)
    root = np.sqrt(sum_squares)
    rv = root - v * np.log(v * root + np.power(v,2)) + v * np.log(u)
    return rv

def new_upper_bound(u,v):
    a = find_intersection(lambda x: iv(v, x), lambda x: 1, np.sqrt(v))
    u_sum_squares = np.power(u,2) + np.power(v,2)
    a_sum_squares = np.power(a,2) + np.power(v,2)
    u_term_root, a_term_root = np.sqrt(u_sum_squares), np.sqrt(a_sum_squares)
    rv = u_term_root - a_term_root - v * np.log(np.divide(a * u_term_root + v, u * a_term_root + v))
    return rv, a

def find_intersection(f1, f2, x0):
    return fsolve(lambda x: f1(x) - f2(x), x0)

def plot_all(u):

    f = plt.figure()

    for i in range(9):
        v = 0.5 * (i+1)
        plt.subplot(3, 3, i + 1)
        plt.tick_params(
            which='both',      
            bottom='off',      
            top='off',
            left='off',
            labelleft='off',
            labelbottom='off')
        log_I = np.log(iv(v, u))
        bd = upper_bound(u,v)
        plt.plot(u, log_I, 'b', u, bd, 'g')
        print( (log_I <= bd).sum() )
    f.suptitle('Bound as in (1) w/ c = 0')

def plot_new(u):

    f = plt.figure()

    for i in range(9):
        v = 0.5 * (i+1)
        plt.subplot(3, 3, i + 1)
        plt.tick_params(
            which='both',      
            bottom='off',      
            top='off',
            left='off',
            labelleft='off',
            labelbottom='off')
        log_I = np.log(iv(v, u))
        bd, a = new_upper_bound(u,v)
        plt.plot(u, log_I, 'b', u, bd, 'g')
        in_range_idx = u <= a
        in_range = log_I[in_range_idx]
        perc = (in_range <= bd[in_range_idx]).sum() / len(in_range)
        print( perc)

    f.suptitle('Bound as in (3)')

if __name__ == '__main__':

    u = np.linspace(0,5,num=1000)[1:]

    plot_all(u)
    plot_new(u)

 A: It would seem that one can apply definite integral inequalities to find piecewise bounds for $\log I_\nu(x)$. The paper from which the inequalities are taken includes both lower and upper bounds:
$$
\sqrt{\frac{\nu}{\nu+1}u^2 + \nu^2} \lt u \frac{I'_\nu(u)}{I_\nu(u)} \lt \sqrt{u^2 + \nu^2} 
$$
Divide all three terms by $u$; call the the left-most side of the resulting inequality $f_\nu(u)$ and the right-most $g_\nu(u)$, $F$ and $G$ their respective anti-derivatives w.r.t. $u$. 
For any $(a, b)$:
$$
F_\nu(b) - F_\nu(a) \lt \log I_\nu(b) - \log I_\nu(a) \lt G_\nu(b) - G_\nu(a)
$$
Setting $a$ such that $I_\nu(a) = 1$ and integrating over $(x,a)$ and $(a,x)$ gives:
$$
\begin{align}
F_\nu(a) - F_\nu(x) \lt - \log I_\nu(x) \Rightarrow \log I_\nu(x) \lt F_\nu(x) - F_\nu(a) && \text{for } x \lt a
\\
\log I_\nu(x) \lt G_\nu(x) - G_\nu(a) && \text{for }  x \gt a
\end{align}
$$
These pass a manual check up to $\nu = 50$ (below). That code simply evaluates the integral, but the bounds can be expressed like so:
$$
\begin{align}
\sqrt{\frac{\nu (x^2 + \nu^2 + \nu)}{\nu+1}} - 
\sqrt{\frac{\nu (a^2 + \nu^2 + \nu)}{\nu+1}} - 
\nu\log \frac{a\big(\sqrt{(x^2+\nu^2 + \nu)} + \sqrt{\nu(\nu+1)} \big)}{x\big(\sqrt{(a^2+\nu^2 + \nu)} + \sqrt{\nu(\nu+1)}\big)}  && \text{for } x \gt a
\\
\sqrt{x^2 + \nu^2} - \sqrt{a^2 + \nu^2} - \nu \log \frac{a(\sqrt{x^2 - \nu^2} + \nu)}{x(\sqrt{a^2 - \nu^2} + \nu)}  && \text{for } x \lt a
\end{align}
$$
Deriving these from $F_\nu$ and $G_\nu$ is left to the interested reader. (I know you're out there, somewhere.)
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 17 14:08:24 2017

@author: seaneaster
"""

import numpy as np
from matplotlib import pylab as plt
from scipy.special import iv
from scipy.optimize import fsolve

def find_intersection(f1, f2, x0):
    return fsolve(lambda x: f1(x) - f2(x), x0, maxfev = 10000)

def lta_def_integral_eval(u, v, c = 0):
    first_term = np.sqrt(v / (v + 1) * (np.power(u, 2) + np.power(v, 2) + v))
    second_term = - v * np.log(v + np.power(v, 2) + np.sqrt(v * (v + 1) * (np.power(u, 2) + np.power(v, 2) + v)))
    third_term = v * np.log(u)
    return first_term + second_term + third_term + c

def gta_def_integral_eval(u, v, c = 0):
    sum_squares = np.power(u, 2) + np.power(v, 2)
    root_sum_squares = np.sqrt(sum_squares)
    return root_sum_squares - v * np.log(v * (root_sum_squares + v)) + v * np.log(u)

def bound_lta(u, v, a):
    return lta_def_integral_eval(u, v) - lta_def_integral_eval(a, v)

def bound_gta(u, v, a):
    return gta_def_integral_eval(u, v) - gta_def_integral_eval(a, v)

def find_a(v, guess = None):
    if guess == None:
        guess = np.power(v, 1.5)
    return find_intersection(lambda x: iv(v, x), lambda x: 1, guess)

def plot_all(u):

    f = plt.figure()

    for i in range(9):
        v = 0.5 * (i+1)
        a = find_a(v)
        ax = plt.subplot(3, 3, i + 1)
        plt.tick_params(
            which='both',      
            bottom='off',      
            top='off',
            left='off',
            labelleft='off',
            labelbottom='off')
        log_I = np.log(iv(v, u))
        lta_bd = bound_lta(u, v, a)
        gta_bd = bound_gta(u, v, a)
        bd = np.where(u < a, lta_bd, gta_bd)
        plt.plot(u, log_I, 'b', u, bd, 'g')
        y1, y2 = ax.get_ylim()
        plt.vlines(a, y1, y2)
        print("v: ", v, "a value: ", a, " perc: ", (log_I <= bd).sum() )
    f.suptitle('Piecewise bound')

def test_bds():
    max_order = 50
    num = 1000
    fmt_str = "Break on v = {}, count under bound (of 1000): {}, a = {}"
    breakdown_order = None
    for v in [(_ + 0.5) for _ in range(max_order + 1)]:
        a = find_a(v)
        u = np.linspace(0, a + 5, num = num + 1)[1:]
        log_I = np.log(iv(v, u))
        lta_bd = bound_lta(u, v, a)
        gta_bd = bound_gta(u, v, a)
        bd = np.where(u < a, lta_bd, gta_bd)
        count_lower = (log_I < bd).sum()
        if count_lower < num:
            print(fmt_str.format(v, count_lower, a))
    if breakdown_order == None:
        print("No break through order {}".format(max_order))

if __name__ == '__main__':

    u = np.linspace(0, 10,num=1001)[1:]
    test_bds()

