Finally, it can be shown that $w(\vz_{0:T})$ provides an unbiased estimator of $Z_1/Z_0$, with $\mathbb{E}[w(\vz_{0:T})] = Z_1/Z_0$ [Neal, 2001].   We can thus estimate the partition function ratio using the empirical average over $K$ annealing chains,  $Z_1/Z_0 \approx \frac{1}{K} \sum w_{0:T}^{(k)}$.  We detail the complete AIS procedure in ais.

AIS is considered among the gold standard methods for estimating normalization constants.

Closely related MCMC methods involving path sampling [Gelman and Meng, 1998] include Sequential Monte Carlo [Del Moral et al., 2006], which may involve resampling steps to prioritize higher-probability $\vz_t$, or parallel tempering [Earl and Deem, 2005], which runs $T$ parallel sampling chains in order to obtain accurate samples from each $\pi_{\beta_t}(\vz)$.

§ $\PHI$-DEFORMED LOGARITHM PATHS

For a strictly positive function $\phi(v)$, Naudts, 2004, Naudts, 2009, Naudts, 2011 define the $\phi$-deformed logarithm as

\begin{align}
\log_{\phi} \of{u} = \int \limits_{1}^{u} \frac{1}{\phi(v)} dv \label{eq:log_phi}
\end{align}

The natural logarithm $\log u$ is recovered for $\phi(v) = v$, and the $q$-deformed logarithm is recovered for $\phi(v) = v^{q}$.
Note that $\log_{\phi}\of{\tpi(\vz)}$ is monotonically increasing and concave in $\tpi(\vz)$, since $\frac{d^2}{du^2}\log_{\phi}(u) = -\frac{1}{\phi(u)^2} < 0$.

In order to define the inverse function $\exp_\phi$ such that $u=\exp_\phi \of{\log_{\phi}(u) }$, it appears we need to know the indefinite integral of $1/\phi(v)$ in log_phi.   Instead, the $\phi$-exponential can be defined using the integral of an additional function $\psi(v) = \frac{d}{dv}\exp_{\phi}(v)$, where

\begin{align}
\exp_{\phi}(u) = 1 + \int \limits_0^u \psi(v)dv \label{eq:phi_exp}
\end{align}

We then have the relationships [Naudts and Zhang, 2018]

\begin{align}
\psi(u) = \phi\of{\exp_{\phi}(u)} \qquad \phi(u) = \psi \of{\log_{\phi}(u)} .
\end{align}

The $\phi$-deformed logarithmic path can now be written as

\begin{align}
\tpi_{\beta}(\vz) = \exp_{\phi}\of{(1-\beta) \log_{\phi}\of{\tpi_0(\vz)} +\beta \log_\phi\of{\tpi_1(\vz)})}
\end{align}

Its gradient with respect to $\vz$ is often necessary for taking MCMC transition steps, for example.

\begin{equation}
\frac{d}{d\vz}  \tpi_{\beta}(\vz) = \phi \of{\tpi_{\beta}} \left( (1-\beta) \frac{\tpi_0(\vz)}{\phi\of{\tpi_0(\vz)}} \frac{d\log \tpi_0(\vz)}{d\vz} + \beta \frac{\tpi_1(\vz)}{\phi\of{\tpi_1(\vz)}} \frac{d\log \tpi_1(\vz)}{d\vz} \right)
\end{equation}

It remains an interesting direction for future work to further explore special cases of the deformed logarithm [Kaniadakis and Scarfone, 2002], or even to learn the deformation function $\phi(u)$ with respect to a quantitative measure of sample quality [Syed et al., 2021].   For the latter approach, parameterizing $\phi(u)$ using a strictly positive neural network would still require integration to obtain $\log_{\phi}$.  If we were to instead directly parameterize $\log_{\phi}(u)$ using an input-convex neural network [Amos et al., 2017], evaluating (one of) the function(s) $\exp_{\phi}$, $\log_{\phi}$, or $\psi$ to evaluate $\tpi_{\beta}(\vz)$ may still be challenging.

§ PARAMETRIC BREGMAN DIVERGENCE FROM NORMALIZATION $Z_{Q}(\BETA)$

Consider a $q$-exponential family with a natural parameter vector $\boldtheta$, sufficient statistic vector $\bt$, and base measure $\expbase$,

\begin{align}
\qexpfam(\vx) &= \frac{1}{Z_q(\boldtheta)} \expbase \exp_q \{ \langle \boldtheta , \bt \rangle \} \qquad  \\
\text{where} \quad Z_q(\boldtheta) &= \int \expbase \exp_q \{ \langle \boldtheta , \bt \rangle \} d\vx . \nonumber
\end{align}

We let $\tqexpfam(\vx) = \expbase \exp_q \{ \langle \boldtheta , \bt \rangle \}$ denote the unnormalized density.

From the convexity of $\exp_q$, it can be shown that the normalization constant $Z_q(\boldtheta)$ is a convex function of the parameters $\boldtheta$, with first derivative equal to

\begin{align}
\frac{\partial Z_q(\boldtheta)}{\partial \theta_i}
%&= \int \frac{\partial }{\partial \theta_i} \expbase [1 + (1-q) \boldtheta \cdot \bt ]^{\frac{1}{1-q}} d\vx \\
&=  \int \expbase [1 + (1-q) \boldtheta \cdot \bt ]^{\frac{q}{1-q}} \cdot T_i(\vx) d\vx = \int \expbase^{1-q} \tqexpfam(\vx)^{q} \cdot T_i(\vx) d\vx
\end{align}

We show that the Bregman divergence induced by $\frac{1}{q} Z_q(\boldtheta)$, for $q > 0$, corresponds to the Amari $\alpha$-divergence between parametric unnormalized densities,

\begin{align}
D_{\frac{1}{q} Z_q}[\bthetaprime: \boldtheta] &= \frac{1}{q} Z_q(\bthetaprime) - \frac{1}{q} Z_q(\boldtheta) - \langle \nabla \frac{1}{q} Z_q(\boldtheta), \bthetaprime - \boldtheta \rangle \\
&=\frac{1}{q} Z_q(\bthetaprime) - \frac{1}{q} Z_q(\boldtheta) - \int \expbase^{1-q} \qexpfam(\vx)^{q} \Big( \langle  \bthetaprime, \bt \rangle  - \langle  \btheta, \bt \rangle \Big)  d\vx \nonumber \\
&\overset{(1)}{=} \frac{1}{q} Z_q(\bthetaprime) - \frac{1}{q} Z_q(\boldtheta) - \frac{1}{q} \int \expbase^{1-q} \qexpfam(\vx)^{q} \Big( \log_q \frac{\tilde{\pi}^{(q)}_{\bthetaprime}(\vx)}{\expbase} -  \log_q \frac{\tqexpfam(\vx)}{\expbase} \Big)  d\vx \nonumber \\
&= \frac{1}{q} Z_q(\bthetaprime) - \frac{1}{q} Z_q(\boldtheta) - \frac{1}{q (1-q)} \int \expbase^{1-q} \qexpfam(\vx)^{q} \Big( \frac{\tilde{\pi}^{(q)}_{\bthetaprime}(\vx)}{\expbase}^{1-q} -   \frac{\tqexpfam(\vx)}{\expbase}^{1-q} \Big)  d\vx \nonumber \\
&= \frac{1}{q} Z_q(\bthetaprime) - \frac{1}{q} Z_q(\boldtheta) - \frac{1}{q (1-q)} \int \Big(  \tilde{\pi}^{(q)}_{\bthetaprime}(\vx)^{1-q} \qexpfam(\vx)^{q} d\vx + \frac{1}{q(1-q)} \int \tqexpfam(\vx)   d\vx \nonumber \\
&= \frac{1}{q}\int \tilde{\pi}^{(q)}_{\bthetaprime}(\vx) d\vx + \frac{1}{1-q} \int \tqexpfam(\vx) d\vx - \frac{1}{q(1-q)} \int   \tilde{\pi}^{(q)}_{\bthetaprime}(\vx)^{1-q} \qexpfam(\vx)^{q} d\vx \nonumber \\
&= D_{A}^{(q)}[\tilde{\pi}^{(q)}_{\bthetaprime} : \tqexpfam ] \nonumber
\end{align}

where in $(1)$ we use the fact that $\log_q \frac{\tqexpfam(\vx)}{\expbase} = \langle \btheta, \bt \rangle$.

For the special case of our parametric interpretation of the $q$-annealing path, we have the (normalized) deformed exponential family

\begin{align}
\pi_{\beta}^{(q)}(\vx) &= \frac{1}{Z_q(\beta)} \tpi_0(\vz) \exp_q \{ \beta \cdot \log_q \frac{\tpi_1(\vx) }{\tpi_0(\vx)} \} .
\end{align}

In particular, between $\beta = 0$ and $\beta = 1$ we have

\begin{align}
D_{\frac{1}{q} Z_q}[0: 1] = D_{A}^{(q)}[\tpi_0 : \tpi_1 ].
\end{align}

§.§ Exponential Family Case

For $q=1$ and the exponential family, we have

\begin{align}
Z(\boldtheta) &= \int \tpi_{\boldtheta}(\vx) d\vx = \int \expbase \exp \{ \langle \boldtheta, \bt \rangle \} d\vx \\
\nabla Z(\boldtheta) &= \int \tpi_{\boldtheta}(\vx) \bt d\vx = \int \expbase \exp \{ \langle \boldtheta, \bt \rangle \} \bt d\vx  \nonumber ,
\end{align}

which leads to the Bregman divergence

\begin{align}
D_{Z}[\bthetaprime : \boldtheta]&= Z(\bthetaprime) - Z(\boldtheta) - \langle \nabla Z(\boldtheta), \bthetaprime - \boldtheta \rangle  \\
&= Z(\bthetaprime) - Z(\boldtheta) - \int \tpi_{\boldtheta}(\vx) \langle \bthetaprime, \bt \rangle d\vx - \int \tpi_{\boldtheta}(\vx) \langle \boldtheta, \bt \rangle d\vx \nonumber \\
&=  \int \tpi_{\bthetaprime}(\vx)d\vx - \int \tpi_{\boldtheta}(\vx)d\vx - \int \tpi_{\boldtheta}(\vx) \Big( \log \tpi_{\bthetaprime}(\vx) - \log \tpi_{\boldtheta}(\vx) \Big) d\vx \nonumber \\
&= \DKL[\tpi_{\boldtheta}(\vx) : \tpi_{\bthetaprime}(\vx)]. \nonumber
\end{align}

By contrast, the divergence generated by the log partition function $\log Z(\boldtheta)$ yields the kl divergence between normalized distributions.  Since

\begin{align}
\nabla \log Z(\boldtheta) = \frac{1}{Z(\boldtheta)} \int \tpi_{\boldtheta}(\vx)\bt d\vx = \int \expfam(\vx) \bt d\vx,  \nonumber
\end{align}

we have

\begin{align}
D_{\log Z}[\bthetaprime : \boldtheta]&= \log Z(\bthetaprime) - \log Z(\boldtheta) - \langle \nabla \log Z(\boldtheta), \bthetaprime - \boldtheta \rangle  \\
&= \log Z(\bthetaprime) - \log Z(\boldtheta) -
\int \expfam(\vx) \langle  \bthetaprime - \boldtheta, \bt \rangle d\vx \nonumber \\
&= \medmath{ \log \int \tpi_{\bthetaprime}(\vx)d\vx - \log \int \tpi_{\boldtheta}(\vx)d\vx - \int \pi_{\boldtheta}(\vx) \Big( \log \tpi_{\bthetaprime}(\vx) - \log \tpi_{\boldtheta}(\vx) \Big) d\vx} \nonumber \\
&= \DKL[\pi_{\boldtheta}(\vx) : \pi_{\bthetaprime}(\vx)] \nonumber
\end{align}

§ LIMITING BEHAVIOR OF RHO-TAU DIVERGENCES

In this section, we show that the limiting behavior of the functions generating the $\alpha$-divergence as a rho-tau Bregman divergence recovers the familiar generators for the divergence.
Similar techniques can be used to derive the limiting behavior for the $\beta$-divergence choice of $ \rho, \tau, f,$ and $f^*$ in limiting_kl.

$f, \rho$ Limiting Behavior as $\mathbf{q\rightarrow 0}$:
Recall the choices of $\rho$ and $f$ which generate the $\alpha$-divergence,

\begin{align}
\rho_{\tpi} &=  \log_q \tpi(\vz) =  \frac{1}{1-q} \tpi(\vz)^{1-q} - \frac{1}{1-q}   \qquad \qquad  \\ % \log_q \tpi(\vz) = \frac{1}{1-q} \big( \tpi(z)^{1-q} - 1 \big)  \\%\qquad
f(\rho) &=\frac{1}{q} \exp_q \bracket{ \rho} - \frac{1}{q} \rho - \frac{1}{q} = \frac{1}{q} [1+(1-q)\rho]_+^{\frac{1}{1-q}} - \frac{1}{q} \rho - \frac{1}{q} .  \nonumber %\label{eq:alpha_rho_app}
\end{align}

Considering the limiting behavior of $f(\rho)$ as $q \rightarrow 0$, note that both the denominator $d(q) = q$ and numerator $n(q) = [1+(1-q)\cdot \rho]^{\frac{1}{1-q}} - \rho - 1$ have $\lim_{q\rightarrow0} d(q) = 0$, $\lim_{q\rightarrow0} n(q) = 0$ and are well defined away from $q=0$.   We can thus use L'Hôpital's rule to find

\begin{align}
\lim \limits_{q\rightarrow 0} f(\rho) = \lim \limits_{q\rightarrow 0} \frac{n^\prime(q)}{d^\prime(q)} =  \frac{1}{1} \cdot n^{\prime}(q) \Big|_{q=0}
\end{align}

Differentiating $n(q) = [1+(1-q)\cdot \rho]^{\frac{1}{1-q}} - \rho - 1$ with respect to $q$ using $\frac{d}{dq}( f(q)^{g(q)} ) = \frac{d}{dq} e^{g(q) \log f(q)}  = f(q)^{g(q)} \left( f^{\prime}(q) \frac{g(q)}{f(q)} + g^{\prime}(q) \log f(q) \right)$, we have

\begin{align}
n^{\prime}(q)\Big|_{q=0} &= \medmath{[1+(1-q) \rho]^{\frac{1}{1-q}} \left(-\rho \frac{\frac{1}{1-q}}{1+(1-q)\rho} + \frac{1}{(1-q)^2}\log \big(1+(1-q) \rho\big)\right) \Big|_{q=0}} \nonumber \\
& = ( 1+ \rho)  \left( - \frac{\rho}{1+\rho} + \log (1 + \rho) \right) \nonumber \\[1.25ex]
% & = (1+\rho) \log (1+\rho) -  \rho \nonumber
\implies \lim \limits_{q\rightarrow 0} f(\rho) &= (1+\rho) \log (1+\rho) - \rho
\end{align}

where we have $\rho_{\tpi}(\vz) = \log_{0}\tpi(\vz) = \tpi - 1$ and $1 + \rho = \tpi(\vz)$.
The function $f(\rho(\tpi)) = \tpi(\vz) \log \tpi(\vz) - \tpi(\vz) + 1$ thus matches the Bregman generator for the negative Shannon entropy and divergence $\DKL[\firstarg:\secondarg]$.
Indeed, the rho-tau Bregman divergence simplifies as

\begin{align}
\rhobreg[&\rho_{\firstarg} : \rho_{\secondarg} ] = \int f\of{\rho_{\firstarg}(\vz)}- f\of{\rho_{\secondarg}(\vz)} - \left( \rho_{\firstarg}(\vz) - \rho_{\secondarg}(\vz) \right) \tau_{\secondarg}(\vz) d\vz \nonumber \\
&=\int \medmath{ \firstarg(\vz) \log \firstarg(\vz) -\firstarg(\vz) - \secondarg(\vz)\log \secondarg(\vz) + \secondarg(\vz) - \Big( \firstarg(\vz) - \secondarg(\vz) \Big) \log {\secondarg}(\vz) d\vz}  \nonumber \\
&= \int \firstarg(\vz) \log \frac{\firstarg(\vz)}{\secondarg(\vz)} d\vz - \int \firstarg(\vz) d\vz + \int \secondarg(\vz) d\vz = \DKL[\firstarg:\secondarg] \label{eq:rho_kl}.
\end{align}

$f, \rho$ Limiting Behavior as $\mathbf{q\rightarrow 1}$:  For $q \rightarrow 1$, we directly reason that $\rho_{\tpi}(\vz) =  \log \tpi(\vz)$ and $f(\rho) = \exp\{\rho\} - \rho - 1$.   This leads to the generator $f(\rho_\tpi(\vz)) = -\log \tpi(\vz) + \tpi(\vz) - 1$, which matches the $f$-divergence generator of the divergence.

$f^*, \tau$ Limiting Behavior as $\mathbf{q\rightarrow 1}$:
Recall the choices of $\tau$ and $f^*$ which generate the $\alpha$-divergence,

\begin{align}
\tau_{\tpi}(\vz) &=  \log_{1-q} \tpi(\vz)  = \frac{1}{q} \, \tpi(\vz)^q - \frac{1}{q} \quad \,\,\\
f^{*}(\tau) &= \frac{1}{1-q} \exp_{1-q} \bracket{ \tau} - \frac{1}{1-q} \tau - \frac{1}{1-q} = \frac{1}{1-q} [1+q \, \tau]_+^{\frac{1}{q}} - \frac{1}{1-q} \tau - \frac{1}{1-q} \nonumber
\end{align}

For $q\rightarrow 1$, we have $\tau_{\tpi}(\vz) = \tpi(\vz) - 1$.

To calculate $\lim_{q\rightarrow 1} f^*(\tau)$, we apply L'Hôpital's rule using similar reasoning as above for $d(q) = \frac{1}{1-q}$ and $n(q) = [1+q\tau]^{\frac{1}{q}} - \tau - 1$.  We differentiate using $\frac{d}{dq}( f(q)^{g(q)} ) = \frac{d}{dq} e^{g(q) \log f(q)}  = f(q)^{g(q)} ( f^{\prime}(q) \frac{g(q)}{f(q)} + g^{\prime}(q) \log f(q) )$ with $g(q) = \frac{1}{q}$, $f(q) = 1 + q\tau$ to obtain

\begin{align}
\lim \limits_{q\rightarrow 1} f^*(\tau) = \lim \limits_{q\rightarrow 1} \frac{n^\prime(q)}{d^\prime(q)} =  \frac{1}{-1} \cdot n^{\prime}(q) \Big|_{q=1} &=  {-
%  \left(
[1+q \tau]^{\frac{1}{q}} \left( \tau \frac{\frac{1}{q}}{1+q\tau} - \frac{1}{q^2}\log \big(1+q \tau \big)
% \Big) - 1
\right) \Big|_{q=1}} \nonumber \\
& = - ( 1+ \tau)  \left( \frac{\tau}{1+\tau} - \log (1 + \tau) \right)\nonumber \\[1.25ex]
% & = (1 + \tau) \log (1+\tau) - \tau  \nonumber
\implies \lim \limits_{q\rightarrow 1} f^*(\tau) &= (1+\tau) \log (1+\tau) - \tau
\end{align}

Using $\tau_{\tpi}(\vz) = \tpi(\vz) - 1$ for $q=1$, we have $f^{*}(\tau_{\tpi}(\vz)) = \tpi(\vz) \log \tpi(\vz) - \tpi(\vz) + 1$, which matches the pointwise negative Shannon entropy or the $f$-divergence generator for $D_f[\firstarg:\secondarg] = \DKL[\firstarg:\secondarg]$.
Identical derivations as in rho_kl confirm that the rho-tau divergence recovers $\taubreg[\tau_{\firstarg}:\tau_{\secondarg}] =\DKL[\firstarg:\secondarg]$.

This observation is indicative of the representational duality with respect to the functions $\rho(u) = \log_q u$ and $\tau(t) = \log_{1-q}t$.   Clearly, switching $q \rightarrow 1-q$ and $1-q \rightarrow q$ switches the role of $f, \rho$ and $f^*, \tau$, so that we obtain the same divergence using $D_{f,\rho}^{(q)}[\firstarg:\secondarg]$ and $ D_{f^*,\tau}^{(1-q)}[\firstarg:\secondarg]$.   For example, we have seen that using $q = 0$ for $\rho(u)$ and $q=1$ for $\tau(t)$ recover the same divergence.

$f^*, \tau$ Limiting Behavior as $\mathbf{q\rightarrow 0}$:  For $q \rightarrow 0$, we similarly find that  $\tau_{\tpi}(\vz) =  \log \tpi(\vz)$ and $f^*(\tau) = \exp\{\tau\} -\tau - 1$.   This leads to the generator $f^*(\tau_\tpi(\vz)) = -\log \tpi(\vz) + \tpi(\vz) - 1$, which matches the $f$-divergence generator of the divergence, with $\DKL[\secondarg:\firstarg] = D_f[\firstarg:\secondarg] =  D_{f,\rho}^{(q=1)}[\firstarg:\secondarg] = D_{f^*,\tau}^{(q=0)}[\firstarg:\secondarg]$.

§ INFORMATION GEOMETRY FROM DIVERGENCE FUNCTIONS

The `Eguchi relations', from the seminal work of Eguchi, 1983, Eguchi, 1985,

show that a statistical divergence $D[\firstargn:\secondargn]$ naturally

induces a Riemannian metric and pair of conjugate affine connections on the tangent space of a manifold $\mathcal{M}$ of probability distributions or densities.
We first review the Eguchi relations for the parametric case, where the parameters
$\boldtheta(\pi): \mathcal{M} \mapsto U \subset \mathbb{R}^N$

provide a coordinate system for points $\pi_{\btheta} \in \mathcal{M}_{\btheta}$ on the manifold, before reviewing the nonparametric information geometry of Zhang, 2013, where each point $\tpi \in \mathcal{M}$ on the manifold is represented by its (infinite-dimensional) density function.

§.§ Parametric Information Geometry

In coordinate notation with

basis vectors of the tangent space $\partial_i = \frac{\partial}{\partial \theta^i}$, the metric

= \langle \partial_i, \partial_j \rangle$ specifies a bilinear form at a point indexed by $\boldtheta$.

An affine connection $\nabla$ defines notions of curvature (of the manifold),
covariant differentiation (of vector fields), and parallel transport (of tangent vectors along curves).

In particular, the covariant derivative can be expressed using the (scalar) Christoffel symbols $\Gamma_{ij,k}(\btheta) = \langle \nabla_{\partial_i} \partial_j, \partial_k \rangle$.

We refer

to Amari and Nagaoka, 2000 or Nielsen, 2020 for detailed background.

For a given divergence function, taking the second and third order differentials yield the following metric and conjugate pair of affine connections

\begin{align}
g_{ij}(\btheta) = - ( {\partial_j})_{\firstargparam} ( {\partial_k} )_\secondargparam  &D[\firstargparam:\secondargparam] \Big|_{\firstargparam=\secondargparam} \label{eq:metric_param}\\
\Gamma_{ij,k}(\btheta) = - ( {\partial_i})_\firstargparam ( {\partial_j})_\firstargparam ( {\partial_k})_\secondargparam &D[\firstargparam:\secondargparam] \Big|_{\firstargparam=\secondargparam} \label{eq:gamma_param} \\
\Gamma^*_{ij,k}(\btheta) = - ( {\partial_i})_\secondargparam ( {\partial_j})_\secondargparam ( {\partial_k})_\firstargparam &D[\firstargparam:\secondargparam] \Big|_{\firstargparam=\secondargparam} \label{eq:gammastar_param}
\end{align}

where $({\partial_j})_{\firstargparam}$ indicates partial differentiation with respect to the parameter $\theta_{a}^j$ of the first argument.  Note that conjugacy of connections amounts to $\Gamma_{ij,k} + \Gamma^*_{ik,j}=\partial_i g_{jk}$ and ensures the inner product of tangent vectors
is preserved by parallel translation of each vector according to the respective connection (Amari and Nagaoka, 2000 Sec. 3.1).

It can be shown that any $f$-divergence [Eguchi, 1983], including the kl divergence in either direction, Amari's $\alpha$-divergence, and the Jensen-Shannon divergence[The JSD, rescaled by $\frac{1}{\beta(1-\beta)}$ to ensure the standard form that $f^{\prime\prime}(1)=1$, is also an $f$-divergence for $   f_{\beta}(u) = \frac{1}{\beta(1-\beta)} \left((1-\beta) u \log u + \big((1-\beta) u + \beta\big)\log \big((1-\beta) u + \beta \big)\right)$.], yields the Fisher Information metric, which is the unique metric that is invariant to reparameterization [Chentsov, 1982],

\begin{align}
g_{ij}(\btheta) =  \int \pi_{\btheta}(\vx) \frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^i} \frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^j} &d\vx
= \int \frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^i} \frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^j} d\vx
\label{eq:metric_log}\\
&= \int  \frac{\partial \log_q \pi_{\btheta}(\vx)}{\partial \theta^i} \frac{\partial \log_{1-q} \pi_{\btheta}(\vx)}{\partial \theta^j} d\vx , \nonumber
\end{align}

where the expression in terms of $q$-representations reflects a duality of representation functions [Amari and Nagaoka, 2000, Zhang, 2004, Nielsen, 2020] which we will study in alpha_div.
However, the divergences we consider will differ in their induced dual affine connections (Eguchi, 1983, Zhang, 2013;connections).

Bregman Divergence and KL Divergence
It is well-known that any Bregman divergence induces a dually flat space,

in which there exists a coordinate system such that either $\Gamma_{ij,k}$ or $\Gamma^{*}_{ij,k}$ vanishes everywhere (Amari and Nagaoka, 2000 Ch. 3,
Nielsen, 2020 Sec. 3.7,
Amari, 2016 Ch. 6).
It is easy to confirm that $\Gamma_{ij,k}(\btheta)=0$ for the parametric Bregman divergence $D_{\psi}[\btheta_a : \btheta_b] = \DKL[\secondargparam : \firstargparam]$ (exp_fam_div associated with the exponential family.

For arbitrary input distributions, differentiating the reverse kl divergence ($\alpha \rightarrow 1$) using $-(\partial_i)_{a} (\partial_j)_a (\partial_k)_b \, \DKL[\secondargn:\firstargn]$, we obtain

\begin{align}
\Gamma_{ij,k}(\btheta) = \Gamma^{(1)}_{ij,k}(\btheta) &\coloneqq \int
\frac{\partial^2 \log \pi_{\btheta}(\vx)}{\partial \theta^i \partial \theta^j}\frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^k} d\vx\label{eq:e}\\
\Gamma^*_{ij,k}(\btheta) = \Gamma^{(0)}_{ij,k}(\btheta) &\coloneqq \int \frac{\partial^2  \pi_{\btheta}(\vx)}{\partial \theta^i \partial \theta^j}\frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^k} d\vx. \label{eq:m}
%   \\[1.25] &=
%   \int
%   \medmath{\pi_{\btheta}(\vx)
%   \left( \frac{\partial^2 \log \pi_{\btheta}(\vx)}{\partial \theta^i \partial \theta^j} + \frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^i}\frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^j}\right)\frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^k}  d\vx} \nonumber
\end{align}

These dual connections are well-known as the $e$- and $m$-connections, respectively [Amari and Nagaoka, 2000].

Similarly, the forward kl divergence ($\alpha \rightarrow 0$) induces $\Gamma_{ij,k} = \Gamma^{(0)}_{ij,k}$ and $\Gamma^*_{ij,k} = \Gamma^{(1)}_{ij,k}$.

$\alpha$-Connection and $f$-Divergence
More generally, Amari, 1982 introduced the family of $\alpha$-connections which, among other interpretations, may be viewed as an arithmetic mixture

of $\Gamma^{(1)}_{ij,k}$ and $\Gamma^{(0)}_{ij,k}$,

\begin{align}
\Gamma^{(\alpha)}_{ij,k}(\btheta) & =
(1-\alpha) \, \Gamma^{(0)}_{ij,k} + \alpha \Gamma^{(1)}_{ij,k} \eqqcolon  \langle \nabla^{(\alpha)}_{\partial_i} \partial_j, \partial_k \rangle  \label{eq:alpha_connection}
%   \qquad \,
%   \Gamma^{*(\alpha)}_{ij,k}(\btheta)=
%   \alpha \Gamma^{(0)}_{ij,k} + (1-\alpha)  \Gamma^{(1)}_{ij,k}
\\[1.25ex]
&= \int \frac{1}{\pi_{\btheta}(\vx)} \left( \frac{\partial^2  \pi_{\btheta}(\vx)}{\partial \theta^i \partial \theta^j} \frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^k} - \alpha \frac{1}{\pi_{\btheta}(\vx)} \frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^i} \frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^j}\frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^k} \right) d\vx \nonumber
%\\[1.25ex]  &= \int (1-\alpha) \frac{\partial^2  \pi_{\btheta}(\vx)}{\partial \theta^i \partial \theta^j}\frac{\partial \log \pi_{\btheta}(\vx)}{\partial \theta^k} + \alpha  \frac{\partial^2 \log \pi_{\btheta}(\vx)}{\partial \theta^i \partial \theta^j}\frac{\partial \pi_{\btheta}(\vx)}{\partial \theta^k} d\vx
\end{align}

where the second line can be derived by simplifying from e-(<ref>).
Adapted to our notation, it can be shown that any $f$-divergence induces a pair of conjugate connections
$\Gamma_{ij,k} =\Gamma^{(\alpha)}_{ij,k}$ and $\Gamma^{*}_{ij,k}= \Gamma^{(1-\alpha)}_{ij,k}$

where $\alpha = -f^{\prime\prime\prime}(1)-1$ (see connections).

In particular, we see that both the $\alpha$-divergence and the (weighted, scaled) Jensen-Shannon divergence induce the $\alpha$-connection.  Although the limiting behavior as $\alpha \rightarrow \{0,1\}$ recovers the kl divergence for both divergences, note that the $f$-divergence generator for the $\alpha$-divergence in amari_alpha is a function of the representation $\rho_q(u) = \log_q(u)$ with $\alpha = q$.   By contrast, the JSD can be derived as either a Jensen diversity or $f$-divergence as a function of the the mixture parameter $\beta$.

§.§ Information Geometry of Rho-Tau Divergence

Parametric Case
Using the Eguchi relations in metric_param-(<ref>), Zhang, 2004, Zhang, 2013 derive the metric and affine connections induced by the representational $\alpha$-divergence or scaled Jensen diversity in scaled_breg.    In particular, the metric becomes

\begin{equation}
% \resizebox{\textwidth}{!}{
% $
\hspace*{-.2cm}  g_{ij}(\boldtheta) = \int \frac{\partial \rho\of{\expfam(\vx)}}{\partial \theta^i} \frac{\partial \tau\of{\expfam(\vx)}}{\partial \theta^j} d\vx = \int
{\rho^{\prime}\of{\expfam(\vx)} \tau^{\prime}\of{\expfam(\vx)}}\frac{\partial \expfam(\vx)}{\partial \theta^i}\frac{\partial \expfam(\vx)}{\partial \theta^j} d\vz, \label{eq:rhotau_param_metric}
% $}
\end{equation}

which generalizes metric_log beyond $\rho(\pi) = \log_q(\pi)$ and $\tau(\pi) = \log_{1-q}(\pi)$.

Zhang, 2004 show that the divergence $\frac{1}{\beta(1-\beta)} \, \BregInfoarg{\potential_{f}}{\rho} \of{\bracket{\tpi_0, \tpi_1},  \bracket{1-\beta, \beta}}$ in zhang_div induces the following analogue of the $\alpha$-connection as its primal connection (gamma_param)

\begin{equation}
\begin{aligned}
\GammaRT{\beta}_{ij,k}(\btheta)&= \int (1-\beta) \frac{\partial^{2} \tau\of{\expfam(\vz)}}{\partial \theta^i \partial \theta^j}\frac{\partial \rho\of{\expfam(\vz)}}{\partial \theta^k}d\vz + \beta \frac{\partial^{2} \rho\of{\expfam(\vz)}}{\partial \theta^i \partial \theta^j}\frac{\partial \tau\of{\expfam(\vz)}}{\partial \theta^k}d\vz . \nonumber
%_{\expfam}(\vx) \medmath
\\ &= {\int  {\rho^{\prime}_{\expfam} \tau^{\prime}_{\expfam}} \left( \frac{\partial^2 \expfam}{\partial \theta^i \partial \theta^j} \frac{\partial \expfam}{\partial \theta^k} + \big( (1-\beta) \frac{\tau^{\prime\prime}_{\expfam}}{\tau^{\prime}_{\expfam}} + \beta \frac{\rho^{\prime\prime}_{\expfam}}{\rho^{\prime}_{\expfam}}  \big) \frac{\partial \expfam}{\partial \theta^i }\frac{\partial \expfam}{\partial \theta^j } \frac{\partial \expfam}{\partial \theta^k}  \right)d\vx} \nonumber
\end{aligned}
\label{eq:rhotau_param_connection}
\end{equation}

which is a mixture $\GammaRT{\beta}_{ij,k}(\btheta) = (1-\beta) \GammaRT{0}_{ij,k}(\btheta) + \beta \GammaRT{1}_{ij,k}(\btheta)$ of the conjugate connections induced by the rho-tau Bregman divergences at $\beta =0$ and $\beta = 1$.
This generalizes the expressions for the $\alpha$-connection $\Gamma^{(\alpha)}_{ij,k}$ in e-(<ref>), which are recovered for $\tau(\pi) = \pi$ and $\rho(\pi) = \log \pi$ (geo_kl).

Nonparametric Statistical Manifolds
Since the target density in MCMC applications is usually  unnormalized and can not be represented using a parametric model, we would like to consider a nonparametric information geometry over the space of
unnormalized measures.
In this case, the approach is to construct a Banach manifold $\mathcal{M}_{\nu_0}$ of all densities which are absolutely continuous with respect to a base measure $\nu_0$ [Pistone and Sempi, 1995, Loaiza and Quiceno, 2013].   Similarly to our likelihood ratio interpretation in q_lref,  the coordinate mappings correspond to $q$-logarithmic likelihood ratios, while the

tangent space can be identified with all one-dimensional
$q$-exponential families with base measure $\tpi_0$ [Loaiza and Quiceno, 2013].   For detailed constructions over normalized nonparametric distributions, see Loaiza and Quiceno, 2013, Loaiza and Quiceno, 2013 or Pistone and Sempi, 1995, Gibilisco and Pistone, 1998, Grasselli, 2010.  Extension to a nonparametric manifold of unnormalized densities should be straightforward, but remains for future work.

Assuming the manifold and tangent spaces are well-defined, Zhang, 2004, Zhang, 2013 obtain the following analogues of the Fisher metric and $\alpha$-connection.   For general tangent vectors $u(\vx), v(\vx), w(\vx)$ at a point $\tpi(\vx)$,

\begin{align}
g_{u,v}({\tpi})&= \langle u, v \rangle
%= \int \frac{1}{\tpi(\vz)} u(\vz) v(\vz) d\vz
%{\rho^{\prime}\of{\tpi(\vx)} \tau^{\prime}\of{\tpi(\vx)}}
\rho^{\prime}_{\tpi}{(\vx)} \tau^{\prime}_{\tpi}{(\vx)} \,
u(\vz) v(\vz) d\vz \label{eq:rhotau_metric} \\
% \end{align}
% \begin{align}
\GammaRT{\beta}_{wu,v}(\tpi) &= \langle \nabla^{(\beta)}_w u, v \rangle    \label{eq:rhotau_connection}  \\
%  =  \int \rho^{\prime} \tau^{\prime}
%  % \of{\tpi(\vx)}
%   (d_w u) v + \big( (1-\alpha) \tau^{\prime\prime} \rho^{\prime}
%   +\alpha \rho^{\prime\prime} \tau^{\prime} \big)  u w v d\vz
%  &=  \int \rho^{\prime}_{\tpi}{(\vx)} \tau^{\prime}_{\tpi}{(\vx)}
%   %\of{\tpi(\vx)}
%   \, d_w u(\vx) v(\vx) d\vx \label{eq:rhotau_connection} \\
%   &\phantom{==} + \int \Big( (1-\alpha) \tau^{\prime\prime}_{\tpi}{(\vx)} \rho^{\prime}_{\tpi}{(\vx)}
%   +\alpha \rho^{\prime\prime}_{\tpi}{(\vx)} \tau^{\prime}_{\tpi}{(\vx)} \Big)  u(\vx) v(\vx) d\vx . \nonumber \\
&=  \medmath{ \int \rho^{\prime}_{\tpi}{(\vx)} \tau^{\prime}_{\tpi}{(\vx)} \left(
\Big( d_w u(\vx) \Big) v(\vx)
+  \Big( (1-\beta)
\frac{\tau^{\prime\prime}_{\tpi}(\vx)}{\tau^{\prime}_{\tpi}(\vx)}
+\beta \frac{\rho^{\prime\prime}_{\tpi}{(\vx)}}{\rho^{\prime}_{\tpi}{(\vx)}}\Big)
%\frac{\tau^{\prime\prime}\of{\tpi{(\vx)}}}{\tau^{\prime}\of{\tpi{(\vx)}} } +\alpha \frac{\rho^{\prime\prime}\of{\tpi{(\vx)}}}{\rho^{\prime}\of{\tpi{(\vx)}}} \big)
u(\vx) w(\vx) v(\vx) \right) d\vx . }
\nonumber
% &\phantom{==}
%   &=  \int \rho^{\prime}\of{\tpi(\vx)} \tau^{\prime}\of{\tpi(\vx)}
%   %\of{\tpi(\vx)}
%   \big(d_w u\big)(\vx) v(\vx) \\
%   &\phantom{==} + \big( (1-\beta) \tau^{\prime\prime}\of{\tpi(\vx)} \rho^{\prime}\of{\tpi(\vx)}
%   +\beta \rho^{\prime\prime}\of{\tpi(\vx)} \tau^{\prime}\of{\tpi(\vx)} \big)  u(\vx) w(\vx) v(\vx) d\vx
\end{align}

where $d_w u$ is the directional derivative of $u$ in the direction of $w$ and $\rho^{\prime}_{\tpi}{(\vx)} = \rho^{\prime}\of{\tpi(\vx)}$.
The parametric expression above can be recovered using, for example, $u(\vx) = \frac{\partial}{\partial \theta^i}\expfam(\vx)$, and we also note the similar form of rhotau_connection to the expression in alpha_connection.
Again, the connections $\GammaRT{\beta}_{wu,v} = (1-\beta) \GammaRT{0}_{wu,v} + \beta \GammaRT{1}_{wu,v}$ are obtained as mixtures of the conjugate connections induced by the rho-tau Bregman divergence.

§ GEODESICS FOR THE RHO-TAU BREGMAN DIVERGENCE

In this section, we show that the quasi-arithmetic mean in the $\rho(\tpi)$ representation of densities is a geodesic with respect to the primal connection induced by the rho-tau divergence $\rhobreg[\rho(\tpi_1):\rho(\tpi_0)]$.

Recall from Zhang, 2013 Sec. 2.3 (Eq. 88) that the $\alpha$-connection associated with the rho-tau Bregman divergence has the form

\begin{align}
\nabla^{(\alpha)}_{\dot{\gamma}} \dot{\gamma} &= (d_{\dot{\gamma}} \dot{\gamma})_{\gamma_t} + \frac{d}{d\gamma} \big( \alpha \log \rho^{\prime}(\gamma_t) + (1-\alpha) \log \tau^{\prime}(\gamma_t) \big) \cdot \dot{\gamma}^2 \label{eq:affine_alpha} \\
&= (d_{\dot{\gamma}} \dot{\gamma})_{\gamma_t} +  \left( \alpha \frac{ \rho^{\prime\prime}(\gamma_t)}{\rho^{\prime}(\gamma_t)} + (1-\alpha) \frac{ \tau^{\prime\prime}(\gamma_t)}{\tau^{\prime}(\gamma_t)} \right) \cdot \dot{\gamma}^2 \nonumber
\end{align}

where $\dot{\gamma} = \frac{d\gamma}{dt}$ and the parameter $\alpha$ plays the role of the convex combination or mixture parameter $\beta$.    In this section, we use the notation $\nabla^{(\alpha)}$ to represent the affine connection, instead of the Christoffel symbol notation from e.g. rhotau_connection, with $\Gamma^{(\alpha)}_{wu,v}(\tpi) = \langle \nabla^{(\alpha)}_w u, v \rangle$.

We are interested in the Bregman divergence and $\nabla^{(1)}$ connection for $\alpha =1$.   Using abstract_affine, we need to show that the geodesic equation $ \nabla^{(1)}_{\dot{\gamma}} \dot{\gamma} = 0$ holds (Nielsen, 2020 Sec. 3.12) for curves which are linear in the $\rho$-representation.

We simplify each of the terms in the geodesic equation, where we rewrite the desired geodesic equation in geodesic_eq to match Eq. 88 of Zhang, 2013,

\begin{align}
\nabla^{(1)}_{\dot{\gamma}} \dot{\gamma}  =
d_{\dot{\gamma}} \dot{\gamma}
%\frac{d \dot{\gamma}}{d\gamma} \cdot \dot{\gamma}
+ \big(\dot{\gamma}\big)^{2} \cdot  \left(\frac{d}{d\gamma} \log {\rho^{\prime}(\gamma)}
%\log \rho^{\prime}(\gamma)
\right) = 0
\end{align}

First, note the particularly simple expression for $\frac{d\rho(\gamma_t)}{dt} = \rho(\tpi_1) - \rho(\tpi_0)$ given the definition $\rho(\gamma(t)) = (1-t) \rho(\tpi_0) + t \, \rho(\tpi_1)$.  Noting the chain rule $\frac{d\rho_t}{dt} = \frac{d\rho(\gamma_t)}{d\gamma} \frac{d\gamma_t}{dt}$, we can rearrange to obtain an expression for $\dot{\gamma}(t)$

\begin{align}
\dot{\gamma}(t) = \frac{d\gamma(t)}{dt} = \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \frac{d\rho(\gamma_t)}{dt} =  \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \left( \rho(\tpi_1) - \rho(\tpi_0) \right)  \label{eq:dot_gamma}
%\frac{1}{\rho(\tpi_1) - \rho(\tpi_0)} \frac{d\rho(\gamma(t))}{dt}
\end{align}

Taking the directional derivative $d_{\dot{\gamma}}\dot{\gamma} = \frac{d\dot{\gamma}}{d\gamma}\cdot \dot{\gamma}$,

\begin{align}
d_{\dot{\gamma_t}}\dot{\gamma_t}  &= \dot{\gamma_t} \cdot \frac{d}{d\gamma} \left[ \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \left( \rho(\tpi_1) - \rho(\tpi_0) \right) \right] \\
&=  - \dot{\gamma_t}\left( \rho(\tpi_1) - \rho(\tpi_0) \right) \frac{d\rho(\gamma_t)}{d\gamma}^{-2} \frac{d^2\rho(\gamma_t)}{d\gamma^2} \nonumber
\end{align}

Rewriting the final term in geodesic_eq, we have

\begin{align}
\frac{d}{d\gamma} \big(  \log \rho^{\prime}(\gamma_t) \big) &= \frac{d}{d\gamma} \big(  \log \frac{d \rho(\gamma_t)}{d\gamma}  \big) =  \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \frac{d^2  \rho(\gamma_t) }{d\gamma^2}.
\end{align}

Putting it all together, we have

\begin{align}
\nabla^{(1)}_{\dot{\gamma_t}} \dot{\gamma_t}  &= -\frac{d \dot{\gamma_t}}{d\gamma} \cdot \dot{\gamma_t} + \big(\dot{\gamma_t} \big)^{2} \cdot \frac{d}{d\gamma} \big(  \log \rho^{\prime}(\gamma_t) \big) \\
&= -\dot{\gamma_t} \left( \rho(\tpi_1) - \rho(\tpi_0) \right) \frac{d\rho(\gamma_t)}{d\gamma}^{-2} \frac{d^2\rho(\gamma_t)}{d\gamma^2} + (\dot{\gamma_t})^2 \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \frac{d^2  \rho(\gamma_t) }{d\gamma^2} \nonumber \\
\intertext{Noting that $\frac{d\rho(\gamma_t)}{d\gamma}^{-1} = \dot{\gamma}(\rho(\tpi_1) - \rho(\tpi_0))$ from \myeq{dot_gamma}, we have }
\nabla^{(1)}_{\dot{\gamma_t}} \dot{\gamma_t} &= -\dot{\gamma_t} \cancel{\frac{\rho(\tpi_1) - \rho(\tpi_0)}{\rho(\tpi_1) - \rho(\tpi_0)}} \cdot \dot{\gamma_t}  \frac{d\rho(\gamma_t)}{d\gamma}^{-1} \frac{d^2\rho(\gamma_t)}{d\gamma^2} + (\dot{\gamma_t})^2 \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \frac{d^2  \rho(\gamma_t) }{d\gamma^2} \nonumber \\
&= -(\dot{\gamma_t})^2  \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \frac{d^2\rho(\gamma_t)}{d\gamma^2} + (\dot{\gamma_t})^2 \left(\frac{d\rho(\gamma_t)}{d\gamma}\right)^{-1} \frac{d^2  \rho(\gamma_t) }{d\gamma^2} \nonumber \\
&=0 \nonumber
\end{align}

which proves the proposition.

We can define an information manifold $(\mathcal{M}_+, G_{u,v}^{(D^{(\alpha)}_{f,\rho})}, \nabla^{(\alpha)}, \nabla^{(-\alpha)})$ from any divergence $D^{(\alpha)}_{f,\rho}$ or $D^{(\beta, q)}$ [Eguchi, 1985].   As shown in Proposition 2 and Corollary 3 of Zhang, 2013, the metric for this manifold has a common form for any rho-tau divergence, which depends on the embedding functions $\rho, \tau$.   Consider two elements of the tangent space to $\mathcal{M}_{+}$ at a point $\tpi$,  $u(\vz), v(\vz) \in \mathcal{T}_{\tpi}\mathcal{M}_+$, where $u$ and $v$ are themselves are functions over $\vz$ (Zhang, 2013 Sec. 2.1).  The metric is given as

\begin{align}
g_{u,v}(\tpi) = \int \rho^{\prime}\big(\tpi(\vz)\big) \tau^{\prime}\big(\tpi(\vz)\big) u(\vz) v(\vz) d\vz \, .
\end{align}

See Zhang, 2013 Sec. 3 (e.g. Proposition 9) for analysis of the parametric case, where the metric above recovers to the Fisher Information.

Recall that, as $\beta \rightarrow 1$,  $D^{(\beta, q)}$ approaches Amari's $\alpha$-divergence, which has the form of a rho-tau Bregman divergence.  A family of affine $\alpha$-connections (here, denoted as $\beta$) can be associated with a rho-tau Bregman divergences.
For representations $\rho, \tau$,

Zhang, 2013 Sec. 2.3 show that the $\alpha$-connection has the form

\begin{align}
\nabla^{(\beta)}_{\dot{\gamma}} \dot{\gamma} = (d_{\dot{\gamma}} \dot{\gamma})_{\gamma_t} + \frac{d}{d\gamma} \big( \beta \log \rho^{\prime}(\gamma(t)) + (1-\beta) \log \tau^{\prime}(\gamma(t)) \big) \cdot \dot{\gamma}^2 \label{eq:affine_alpha}
\end{align}

where $ (d_{\dot{\gamma}} \dot{\gamma})_{\gamma_t} = \frac{d \dot{\gamma(t)}}{d\gamma_t} \cdot \dot{\gamma}$.

We proceed to show that our $q$-path forms a geodesic with respect to the affine connection $\nabla^{(1)}$ in the information manifold induced by $D^{(\beta=1, q)} = D^{(\alpha=q)}_{A}$.

Let $\gamma_t = \gamma(t) = \tpi_{t,q}(z)$ indicate the $q$-path distribution with mixing parameter $t$,

with $\dot{\gamma} = \frac{d \gamma(t)}{dt}$.   We need to show that (Nielsen, 2020 Sec. 3.12)

\begin{align}
\nabla_{\dot{\gamma}} \dot{\gamma} = 0.
\end{align}

Plugging $\beta= 1$ into (<ref>) to express the desired covariant derivative,
we need to show that $\nabla_{\dot{\gamma}} \dot{\gamma} = 0$, or

\begin{align*}
\nabla^{(1)}_{\dot{\gamma}} \dot{\gamma} = \frac{d \dot{\gamma}}{d\gamma} \cdot \dot{\gamma} + \big(\dot{\gamma}\big)^{2} \cdot \frac{d}{d\gamma} \big(  \log \rho^{\prime}(\gamma(t)) \big)
\end{align*}

We will first derive several useful quantities.
First, we obtain the desired $\frac{d\gamma(t)}{dt}$ using the chain rule $\frac{d}{dt} \rho\big( \gamma(t) \big) = \frac{d \rho(\gamma) }{d\gamma(t)} \frac{d\gamma(t)}{dt}$, where it is easy to differentiate $\rho(\gamma(t)) = (1- t) \rho(\tpi_0) + t \, \rho(\tpi_1)$ with respect to $\gamma$,

\begin{align}
\dot{\gamma} = \frac{d\gamma(t)}{dt} &= \bigg( \frac{d \rho(\gamma_t) }{d\gamma(t)} \bigg)^{-1} \, \frac{d \rho(\gamma_t) }{dt} \label{eq:dgamma_dt_precursor}
% \\[1.4ex]
%  &= \gamma^{\alphaaneg} \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) \label{eq:dgamma_dt} %\\[1.4ex]
% \implies  \dot{\gamma} &=\gamma_t^}\alphaaneg \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) \quad (\text{for $\alpha = 1$}) \label{eq:dgamma_dt}
\end{align}

We can now explicitly derive the above quantities for the case of $q$-paths and $\rho(u) = \log_q u$, with

\begin{align}
\frac{d \rho(\gamma_t) }{d\gamma_t}  &= \frac{d}{d \gamma_t} \alphaainv \big( \gamma_t^\alphaa - 1 \big) = \gamma_t^{\alphaa - 1} =  \gamma_t^{-\alphaaneg} \\
\text{and } \quad \frac{d \rho(\gamma_t) }{dt}  &= \frac{d}{dt} \big( (1-t) \rho(\tpi_0) + t \, \rho(\tpi_1) \big) = \rho(\tpi_1) - \rho(\tpi_0).
\end{align}

Plugging these back into dgamma_dt_precursor

\begin{align}
\dot{\gamma}= \frac{d\gamma(t)}{dt}  &= \gamma^{\alphaaneg} \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) \label{eq:dgamma_dt}
\end{align}

We now differentiate the $\frac{d}{d\gamma} \big(  \log \rho^{\prime}(\gamma(t)) \big) $ term in the expression for the affine connection in affine_alpha .   From above, $d \rho(\gamma_t) / d\gamma(t) = \gamma_t^{-q}$, so that

\begin{align}
\frac{d}{d\gamma} \big(  \log \frac{d \rho(\gamma(t))}{d\gamma(t)}  \big) &=   \frac{d}{d\gamma}   \log \big( \gamma_t^{-q} \big) \\
&= -q \cdot \frac{d}{d\gamma}   \log \gamma_t \\
%&\phantom{\frac{d}{d\gamma}   \log \big( \gamma_t^{-q} \big)} = -q \cdot \gamma_t^{-1} \\
% &= \gamma_t^{\frac{(1+\alpha)}{2}} \cdot \frac{d}{d\gamma} \gamma_t^{\frac{(1+\alpha)}{2}} \\
% % &=
% % \bigg( \frac{d}{d\gamma} \rho(\gamma_t) \bigg)^{-1} \cdot \frac{d}{d\gamma}  \rho(\gamma_t) \\
% &= \gamma^{\frac{(1+\alpha)}{2}} \cdot  - \frac{(1+\alpha)}{2} \gamma^{-\frac{(1+\alpha)}{2}-1} \\
% &= -\gamma_t \cdot \gamma_t^{-2} \\
\implies \frac{d}{d\gamma} \big(  \log \rho^{\prime}(\gamma(t)) \big) &=  -q \cdot \gamma_t^{-1} \label{eq:db_dgamma}
\end{align}

Putting it all together, and with dot products turning into scalar products for our one-dimensional curve $t$, we have

\begin{align}
\nabla^{(1)}_{\dot{\gamma}} \dot{\gamma} &=  \frac{d \dot{\gamma}_t}{d\gamma_t} \cdot \dot{\gamma}_t  + \big(\dot{\gamma}\big)^{2} \frac{d}{d\gamma} \big(  \log \rho^{\prime}(\gamma(t)) \big) \\
&=  \frac{d }{d\gamma_t} \bigg( \gamma_t^{q} \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) \bigg) \cdot \dot{\gamma}_t  -  q \cdot (\dot{\gamma}_t)^{2} \cdot \gamma_t^{-1}
\qquad \quad \hfill \text{(using \eqref{eq:dgamma_dt} and \eqref{eq:db_dgamma})} \\
&= q \cdot \gamma_t^{q-1}  \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) \cdot \dot{\gamma}_t -  q \cdot (\dot{\gamma})^{2} \cdot \gamma_t^{-1} \\ %\qquad \text{(rearranging \eqref{eq:dgamma_dt})} \\
&= q \cdot \gamma_t^{q-1}  \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) -  q \cdot \dot{\gamma}_t \cdot \gamma_t^{-1}  \qquad \qquad \hfill \text{(Factor out  } \dot{\gamma}_t) \\
&= q \cdot \gamma_t^{q-1}  \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) -  q \cdot \gamma_t^{q}  \cdot \big( \rho(\tpi_1) - \rho(\tpi_0) \big) \cdot \gamma_t^{-1}  \quad \hfill \text{(Substituting \eqref{eq:dgamma_dt})} \nonumber \\
\end{align}


#  Efficient and accurate computation to the $\varphi$-function and its action
on a vector

Siyu Yang  Dongping Li<EMAIL_ADDRESS>Department of Mathematics, Changchun
Normal University, Changchun 130032, PR China  Department of Mathematics,
Jilin University, Changchun 130012, PR China

###### Abstract

In this paper, we develop efficient and accurate algorithms for evaluating
$\varphi(A)$ and $\varphi(A)b$, where $A$ is an $N\times N$ matrix, $b$ is an
$N$ dimensional vector and $\varphi$ is the function defined by
$\varphi(z)\equiv\sum\limits^{\infty}_{k=0}\frac{z^{k}}{(1+k)!}$. Such matrix
function (the so-called $\varphi$-function) plays a key role in a class of
numerical methods well-known as exponential integrators. The algorithms use
the scaling and modified squaring procedure combined with truncated Taylor
series. The backward error analysis is presented to find the optimal value of
the scaling and the degree of the Taylor approximation. Some useful techniques
are employed for reducing the computational cost. Numerical comparisons with
state-of-the-art algorithms show that the algorithms perform well in both
accuracy and efficiency.

###### keywords:

$\varphi$-function , Truncated Taylor series , Scaling and modified squaring
method , Backward error, Paterson-Stockmeyer method

###### MSC:

[2010] 65L05 , 65F10, 65F30

left

##  1 Introduction

In this work, we consider numerical methods for approximating the first matrix
exponential related function and its action on a vector, that is,

$\varphi(A)~{}~{}\text{and}~{}~{}\varphi(A)b,$ (1)

where

$\varphi(z)=\sum\limits^{\infty}_{k=0}\frac{z^{k}}{(1+k)!},~{}A\in\mathbb{C}^{N\times
N},~{}b\in\mathbb{C}^{N}.$ (2)

The $\varphi$-function satisfies the recursive relation

$\varphi(z)=\frac{e^{z}-1}{z}.$ (3)

The problem of numerically approximating such matrix function is of great
importance and is commonly encountered in the solution of constant
inhomogeneous linear system of ordinary differential equations and in the
exponential integrators for solving semi-linear problems. For example, the
well-known exponential Euler method for solving the autonomous semi-linear
problems of the form

$y^{\prime}(t)=Ay(t)+N(y(t)),~{}~{}y(t_{n})=y_{n}$ (4)

yields

$y_{n+1}=e^{hA}y_{n-1}+h\varphi(hA)N(y_{n}).$ (5)

If Eq. (4) has a constant inhomogeneous term, i.e., $N(y(t))\equiv b$, then
the scheme (3) is the exact solution of (4). Utilizing the relationship (3),
it can be shown that (5) is equivalent to

$y_{n+1}=h\varphi(hA)(AN(y_{n})+y_{n-1}).$ (6)

The main cost in the scheme (6) originates from the need to accurately solve
the $\varphi$-function at each time step. For a detailed overview on
exponential integrators, see [13, 17].

Over the past few years, there has been a tremendous effort to develop
efficient approaches to deal with such matrix functions, see, e.g., [2, 3, 4,
12, 25, 16, 14, 20, 26]. These methods are generally divided into two classes.
The first class of methods compute $\varphi(A)$ explicitly. Among them, the
scaling and modified squaring method combined with Padé approximation [12, 26]
is perhaps the most popular choice for small and medium sized $A$. The method
is a variant of the well-known scaling and squaring approach for computing the
matrix exponential [1, 9]. An alternative computation is based on the formula
[22]:

$\begin{aligned}
e^{\mathbb{A}}=\left(\begin{tabular}[]{cccccc}$e^{A}$&$\varphi(A)$\\\
$0$&$I$\end{tabular}\right),\end{aligned}~{}\text{where}~{}\begin{aligned}
\mathbb{A}=\left(\begin{tabular}[]{cccccc}$A$&$I$\\\
$0$&$0$\end{tabular}\right)\in\mathbb{C}^{2N\times 2N}.\end{aligned}$ (7)

Thus the computation of $\varphi(A)$ can be reduced to that of the matrix
exponential. The effective evaluation of matrix exponential, which arise in
many areas of science and engineering, have been extensively investigated in
the literature; see, e.g., [1, 6, 7, 9, 10, 18, 22, 24, 27] and the references
given therein.

In some applications, it requires the computation of matrix-function vector
product $\varphi(A)b$ rather than the single $\varphi(A)$. When $A$ is very
large, it is prohibitive to explicitly compute $\varphi(A)$ and then form the
the product with vector $b$. The second class of methods enable evaluation of
$\varphi(A)b$ using matrix-vector products and avoids the explicit computation
of the generally dense matrix $\varphi(A)$. This type of methods is especially
well-suited to large and sparse $A$. We mention two typical strategies in such
an approach: Krylov subspace methods [25, 20] and the scaling-and-squaring
method [2]. The former are iterative and difficult to determine a reasonable
convergence criterion to guarantee a sufficiently accurate approximation. The
latter evaluate $\varphi(A)b$ by computing the action of a matrix exponential
$e^{\mathbb{A}}$ of dimension $N+1$ on a vector. The method is numerical
stable and can achieve a machine accuracy in exact arithmetic.

In the present paper we focus on the direct approach and develop the scaling
and modified squaring method in combination with Taylor series to efficiently
and accurately evaluate $\varphi(A)$ and $\varphi(A)b,$ respectively. The
backward error are used to determine the scaling value $s$ and the Taylor
degree $m$. Numerical experiments with other state-of-the-art MATLAB routines
illustrate that a straight implementation of the scaling and modified squaring
algorithm may be the most efficient.

This paper is organized as follows. Section 2 presents two algorithms for
computing $\varphi(A)$. Section 3 deals with algorithm for evaluating
$\varphi(A)b.$ Numerical experiments are given to illustrate the benefits of
the algorithms in Section 4. Finally, conclusions are given in Section 5.

Throughout the paper, we use $\|\cdot\|$ to denote an induced matrix norm, and
in particular $\|\cdot\|_{1}$, the 1-norm. Let $I$ be the identity and 0 be
the zero matrix or vector whose dimension are clear from the context. $e_{i}$
denotes the $i$-th coordinate vector with appropriate size. $\lfloor x\rfloor$
denotes the largest integer not exceeding $x$ and $\lceil x\rceil$ denotes the
smallest integer not less than $x$. Standard MATLAB notations are used
whenever necessary.

##  2 Computing $\varphi(A)$

For a given matrix $A\in\mathbb{C}^{N\times N},$ the scaling and modified
squaring method exploits the identity [12, 26]

$\varphi(A)=\frac{1}{2}\varphi(\frac{1}{2}A)(e^{\frac{1}{2}A}+I).$ (8)

Applying recursively (8) $s$ times yields

$\displaystyle\varphi(A)=(\frac{1}{2})^{s}\varphi(X)(e^{X}+I)(e^{2X}+I)\ldots(e^{2^{s-1}X}+I),~{}~{}s\in\mathbb{N},$
(9)

where $X=2^{-s}A$. The $\varphi(A)$ then can be evaluated using rational
polynomial to approximate $\varphi(X)$ and $e^{X}$ and employing the following
coupled recurrences:

$\left\\{\begin{array}[]{l}\varphi(2X)=\frac{1}{2}\varphi(X)(e^{X}+I),\vspace{1ex}\\\
e^{2X}=e^{X}\cdot e^{X}.\end{array}\right.$ (10)

The scaling parameter $s$ is chosen such that $\|X\|$ is sufficiently small
and the method can achieve a prescribed accuracy.

In our algorithm, we use the truncated Taylor series $T_{m}(X)$ to approximate
$\varphi(X)$, i.e.,

$T_{m}(X):=\sum\limits^{m}_{k=0}\frac{X^{k}}{(1+k)!}.$ (11)

Then, the approximation to $e^{X}$ is naturally chosen as

$\tilde{T}_{m}:=XT_{m}(X)+I.$ (12)

Here $e^{X}=\tilde{T}_{m}(X)+\mathcal{O}(X^{m+2})$. The computation of
$\tilde{T}_{m}$ only requires one matrix multiplication and one matrix
summation.

In practical the truncated Taylor series $T_{m}(X)$ in (11) can be computed
efficiently by using the Paterson-Stockmeyer method [21]. The expression is

$\tilde{T}_{m}=\sum\limits^{r}_{0}B_{k}\cdot(X^{q})^{k},~{}~{}~{}r=\lfloor
m/q\rfloor,$ (13)

where $q$ is a positive integer and

$B_{k}=\left\\{\begin{array}[]{l}\sum\limits^{q-1}_{i=0}\frac{1}{(1+qk+i)!}X^{i},~{}~{}~{}k=0,1,\ldots,r-1,\vspace{1ex}\\\
\sum\limits^{m-qr}_{i=0}\frac{1}{(1+qr+i)!}X^{i},~{}~{}~{}k=r.\\\
\end{array}\right.$ (14)

Applying Horner’s method to (13), then the number of matrix multiplications
for computing $T_{m}(X)$ is minimized by $q$ either $\lfloor\sqrt{m}\rfloor$
or $\lceil\sqrt{m}\rceil$, and both choices yield the same computational cost.

To obtain a more accurate approximation to $\varphi(X)$, we compute $T_{m}$
only when $m$ belongs to the index sequence
$\mathbb{M}=\\{2,4,6,9,12,16,20,25,30,36,\ldots\\}$. Assume that $m_{i}$ is
the $i$-th element of the set $\mathbb{M}$, it is shown in [10, Table 4.1],
[23] that the number of matrix multiplications for computing $T_{m_{i}}(X)$ is
the same amount as $T_{k}(X)$ of $m_{i-1}<k<m_{i}$. Then the number of matrix
multiplications for computing $T_{m}(X)$ is

$\pi_{m}=\lceil\sqrt{m}\rceil+\lfloor m/\lceil\sqrt{m}\rceil\rfloor-2.$ (15)

Table 1 lists the corresponding number of matrix multiplications $\pi_{m}$ to
evaluate $T_{m}$ for the first 12 values of $m$ belonging to $\mathbb{M}$. A
brief sketch of the algorithm for solving $\varphi(A)$ is given in Algorithm
1.

Table 1: Number of matrix multiplications $\pi_{m}$ required to evaluate $T_{m}$ for the first 12 optimal values of $m$.  $m$ | $2$ | $4$ | $6$ | $9$ | $12$ | $16$ | $20$ | $25$ | $30$ | $36$ | $42$ | $49$
---|---|---|---|---|---|---|---|---|---|---|---|---
$\pi_{m}$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$ | $9$ | $10$ | $11$ | $12$
Algorithm 1 Given $A\in\mathbb{C}^{N\times N},$ this algorithm computes
$\varphi(A)$ by the scaling and modified squaring based on Taylor series.

1:Select optimal values of $m$ and $s.$

2:Compute $X=2^{-s}A.$

3:Compute $T=\sum\limits^{m}_{k=0}\frac{X^{k}}{(1+k)!}$ by PS method.

4:Compute $\tilde{T}:=XT+I.$

5:Compute $Y=\tilde{T}+I.$

6:for $i=1:s-1$ do

7: Compute $\tilde{T}=\tilde{T}^{2}.$

8: Compute $Y=\frac{1}{2}Y(\tilde{T}+I).$

9:end for

10:Compute $Y=\frac{1}{2}TY.$

11: $Y$

Now we consider the concrete choice of $m$ and $s$. We formulate two
approaches to choose the scaling value $s$ and the Taylor degree $m$, which
were similarly introduced in [1, 2]. Define the function
$h_{m+2}(X)=\log(e^{-X}\tilde{T}_{m}(X))$, then

$\tilde{T}_{m}(X)=e^{X+h_{m+2}(X)}$ (16)

and

$\displaystyle\begin{aligned}
\varphi(A)&\approx(\frac{1}{2})^{s}T_{m}(\tilde{T}_{m}+I)(\tilde{T}_{m}^{2}+I)\cdots(\tilde{T}_{m}^{2^{s-1}}+I)\\\
&=(\frac{1}{2})^{s}\frac{\tilde{T}_{m}-I}{X}(\tilde{T}_{m}+I)(\tilde{T}_{m}^{2}+I)\cdots(\tilde{T}_{m}^{2^{s-1}}+I)\\\
&=\frac{e^{2^{s}X+2^{s}h_{m+2}(X)}-I}{2^{s}X}\\\ &=\frac{e^{A+\Delta
A}-I}{A},\end{aligned}$ (17)

where $\Delta A=2^{s}h_{m+2}(X)$ is the backward error resulting from the
approximation of $\varphi(A).$

Let $X\in\Omega_{m}:=\\{X\in C^{N\times
N}:~{}~{}\rho(e^{-X}\tilde{T}_{m}-I)<1\\},$ then the function $h_{m+2}(X)$ has
a power series expansion

$h_{m+2}(X)=\sum\limits^{\infty}_{k=m+2}c_{k}X^{k}.$ (18)

By Theorem 4.2 of [1] we have

$\frac{\|\Delta
A\|}{\|A\|}=\frac{\|h_{m+2}(X)\|}{\|X\|}\leq\tilde{h}_{m+2}(2^{-s}\alpha_{p}(A)),~{}~{}p(p-1)\leq
m+2,$ (19)

where $\tilde{h}_{m+2}(x)=\sum\limits^{\infty}_{k=m+2}|c_{k}|x^{k-1}$ and
$\alpha_{p}(A)=\max(\|A^{p}\|^{1/p},\|A^{p+1}\|^{1/(p+1)}).$ Given a tolerance
Tol, one can computes

$\theta_{m}=max\\{\theta:{\tilde{h}_{m+2}(\theta)}\leq\text{Tol}\\}.$ (20)

Table 2 presents the maximal values $\theta_{m}$ satisfying the backward error
bound (20) of $\text{Tol}=2^{-53}$ for the first 12 values of $m$ in
$\mathbb{M}$. Thus, once the scaling $s$ is chosen such that

$2^{-s}\alpha_{p}(A)\leq\theta_{m},~{}~{}p(p-1)\leq m+2,$ (21)

it follows that

$\|\Delta A\|\leq\|A\|\cdot\text{Tol}.$ (22)

A straightforward computation of inequality (21) yields

$s\geq\lceil\log_{2}(\alpha_{p}(A)/\theta_{m})\rceil.$ (23)

We naturally choose the smallest $s$ so that the inequality (21) holds. The
total number of matrix multiplications $C_{m}$ to evaluate $\varphi(A)$ then
is

$C_{m}=\pi_{m}+2s=\lceil\sqrt{m}\rceil+\lfloor\sqrt{m}\rfloor-2+2\max(\lceil\log_{2}(\alpha_{p}(A)/\theta_{m})\rceil,~{}0).$
(24)

Table 2: Maximal values $\theta_{m}$ such that the backward error bound (20) does not exceed $\text{tol}=2^{-53}$ for the first 12 optimal values of $m$. $m$ | $2$ | $4$ | $6$ | $9$ | $12$ | $16$ | $20$ | $25$ | $30$ | $36$ | $42$ | $49$
---|---|---|---|---|---|---|---|---|---|---|---|---
$\theta_{m}$ | $1.39\text{e-}5$ | $2.40\text{e-}3$ | $2.38\text{e-}2$ | $1.44\text{e-}1$ | $4.00\text{e-}1$ | $9.31\text{e-}1$ | $1.62$ | $2.64$ | $3.77$ | $5.22$ | $6.73$ | $8.55$

In Figure 1 we have plotted $C_{m}$ as a function of $m$ for ten different
values of $\alpha_{p}(A)$. We see the location of the first optimal value of
$m$, that is, the first value that minimizes $C_{m}$ is no more than 25. Thus
we consider $m$ with $m\in\\{2,4,6,9,12,16,20,25\\}$ in the remainder of the
section.

Figure 1: $m$ versus cost $C_{m}$ with different $\alpha_{p}.$

In order to get the optimal value of $m$, we consider the following two
strategies:

$\bullet$ Choose the first $m\in\\{2,4,6,9,12,16,20,25\\}$ such that
$\eta_{m}\leq\theta_{m}$, where $\eta_{m}=\min\\{\alpha_{p}(A),p(p-1)\leq
m+2\\}$, and set $s=0$. When $\eta_{25}>\theta_{25}$, set $m=25$ and
$s=\lceil\log_{2}(\eta_{25}/\theta_{25})\rceil$. To reduce the computational
cost, in practical implementation the bound $\|A^{p}\|^{1/p}$ are estimated
using the products of bounds or norms of matrices that have been computed. The
details of the process are summarized in Algorithm 2.

$\bullet$ Select the parameters $m$ and $s$ such that the total computational
cost (23) is the lowest. This requires pre-evaluating the first six 1-norm of
matrix power, i.e., $\|A^{k}\|$, $k=1,2,\cdots,6$. The full procedure is given
in Algorithm 3.

Algorithm 2 Given $A\in\mathbb{C}^{N\times N},$ this algorithm computes the
parameters $m,~{}s$ and $A_{i}=A^{i}.$

1:$s=0.$

2:$A_{1}=A,$ $A_{2}=A^{2},$ $d_{1}=\|A_{1}\|_{1},$ $d_{2}=\|A_{2}\|_{1},$
$d_{3}=d_{1}d_{2}.$

3:$\alpha_{1}=d_{1},$ $\alpha_{2}=\max(d_{2}^{1/2},d_{3}^{1/3}),$
$\eta_{1}=\alpha_{2}.$

4:if $\eta_{1}<=\theta_{2},$ then $m=2,$ return .

5:end if

6:if $\eta_{1}<=\theta_{4},$ then $m=4,$ return .

7:end if

8:$A_{3}=A_{1}A_{2},$ $d_{3}=\|A_{3}\|_{1},$ $d_{4}=\min(d_{1}d_{3},d_{4}),$
$\alpha_{2}=\max(d_{2}^{1/2},d_{3}^{1/3}),$
$\alpha_{3}=\max(d_{3}^{1/3},d_{4}^{1/4}).$

9:$\eta_{2}=\min(\alpha_{2},\alpha_{3}).$

10:if $\eta_{2}<=\theta_{6},$ then $m=6,$ return .

11:end if

12:if $\eta_{2}<=\theta_{9},$ then $m=9,$ return .

13:end if

14:$A_{4}=A_{2}^{2},$ $d_{4}=\|A_{4}\|_{1},$
$d_{5}=\min(d_{1}d_{4},d_{2}d_{3}),$
$\alpha_{3}=\max(d_{3}^{1/3},d_{4}^{1/4}),$
$\alpha_{4}=\max(d_{4}^{1/4},d_{5}^{1/5}).$

15:$\eta_{3}=\min(\alpha_{2},\alpha_{3},\alpha_{4}).$

16:if $\eta_{3}<=\theta_{12},$ then $m=12,$ return .

17:end if

18:if $\eta_{3}<=\theta_{16}$ then $m=16,$ return .

19:end if

20:$A_{5}=A_{1}A_{4},$ $d_{5}=\|A_{5}\|_{1},$
$d_{6}=\min(d_{1}d_{5},d_{2}d_{4},d_{3}^{2}),$
$\alpha_{4}=\max(d_{4}^{1/4},d_{5}^{1/5}),$
$\alpha_{5}=\max(d_{5}^{1/5},d_{6}^{1/6}).$

21:$\eta_{4}:=\min(\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5}).$

22:if $\eta_{4}<=\theta_{20}$ then $m=20,$ return .

23:end if

24:if $\eta_{4}<=\theta_{25}$ then $m=25,$ return .

25:end if

26:$m=25,$ $s=\lceil\log_{2}(\eta_{4}/\theta_{25})\rceil,$

27:$A_{i}=2^{-is},$ $i=1,\cdots,5.$

Algorithm 3 Given $A\in\mathbb{C}^{N\times N},$ this algorithm computes the
parameters $m$ and $s$ based on the number of matrix-matrix products.

1:$M=[2,~{}4,~{}6,~{}9,~{}12,~{}16,~{}20,~{}25].$

2:$p_{max}=5,$ $m_{max}=8.$

3:$A_{1}=A,$ $d_{1}=\|A\|_{1}.$

4:for $p=2:p_{max}+1$ do

5: $c=\text{normest}(A,p).$

6: $d_{p}=c^{1/p}.$

7:end for

8:$\alpha_{1}=d_{1}.$

9:for $p=2:p_{max}$ do

10: $\alpha_{p}=\max(d_{p},d_{p+1}).$

11:end for

12:$\eta_{1}=\alpha_{2}.$

13:for $p=2:p_{\max}$ do

14: $\eta_{p}=\min(\eta_{p-1},\alpha_{p}).$

15:end for

16:for $m=[2,~{}4,~{}6,~{}9,~{}12,~{}16,~{}20,~{}25]$ do

17: if $m=2$ then

18: $s_{m}=\max(\lceil\log_{2}(\eta_{2}/\theta_{m})\rceil,0).$

19: else if $m=[4,6,9]$ then

20: $s_{m}=\max(\lceil\log_{2}(\eta_{3}/\theta_{m})\rceil,0).$

21: else if $m=[12,16]$ then

22: $s_{m}=\max(\lceil\log_{2}(\eta_{4}/\theta_{m})\rceil,0).$

23: else

24: $s_{m}=\max(\lceil\log_{2}(\eta_{5}/\theta_{m})\rceil,0).$

25: end if

26: $q_{m}=\sqrt{m},$ $C_{m}=\lceil q_{m}\rceil+\lfloor
q_{m}\rfloor-2+2*s_{m}.$

27:end for

28:$m=\text{argmin}_{m\in M}C_{m},$ $s=s_{m},$ $q=\lceil q_{m}\rceil.$

29:$A_{i}=A_{i-1}A,$ $i=2:q.$

30:$A_{i}=2^{-is}A_{i},$ $i=2:q.$

##  3 Computing $\varphi(A)b$

We now focus our attention on accurately and efficiently evaluating
$\varphi(A)b$ for sparse and large matrix $A$. Following an idea of Al-Mohy
and Higham [2], we will use the scaling part of the scaling and modified
squaring method in combination with truncated Taylor series to approximate the
function. The computational cost of the method is dominated by matrix-vector
products.

We start by recalling the following general recurrence [26]:

$\varphi((\alpha+\beta)z)=\frac{1}{(\alpha+\beta)}[\beta e^{\alpha
z}\varphi(\beta z)+\alpha\varphi(\alpha
z)],~{}~{}~{}\alpha,\beta\in\mathbb{R},~{}z\in\mathbb{C}.$ (25)

As a special case of (25), we have

$\varphi(sz)=\frac{1}{s}[e^{(s-1)z}\varphi(z)+(s-1)\varphi((s-1)z)],~{}~{}s\in\mathbb{N}.$
(26)

Taking $Y=\frac{1}{s}A,$ and using (26) it follows that

$\displaystyle\begin{aligned}
\varphi(A)=&\frac{1}{s}[e^{(s-1)Y}\varphi(Y)+(s-1)\varphi((s-1)Y)]\\\
=&\frac{1}{s}[e^{(s-1)Y}\varphi(Y)+e^{(s-2)Y}\varphi(Y)+(s-2)\varphi((s-2)Y)]\\\
&\vdots\\\
=&\frac{1}{s}\varphi(Y)[e^{(s-1)Y}+e^{(s-2)Y}+\cdots+e^{2Y}+e^{Y}+I].\end{aligned}$
(27)

Choose the integers $m$ and $s$ such that $\varphi(Y)$ and $e^{Y}$ can be
well-approximated by the truncated Taylor series $T_{m}(Y)$ and
$\tilde{T}_{m}(Y):=YT_{m}(Y)+I$ defined by (11) and (12). Then $\varphi(A)b$
can be approximated by firstly evaluating the recurrence

$\displaystyle
b_{1}=T_{m}(Y)b~{}~{}~{}\text{and}~{}~{}~{}b_{i+1}=\tilde{T}_{m}(Y)b_{i},~{}~{}i=1,2,\cdots,s-1,$
(28)

and then computing $\frac{1}{s}$ times the sum of $b_{i},i=1,2,\cdots,s-1.$
This process requires $s$ multiplications of matrix polynomial with a vector,
$s$ vector additions, and 1 scalar multiplication. The number of matrix-vector
products for evaluating $\varphi(A)b$ by recurrence (28) is $C_{m}=s(m+1)-1.$

Algorithm 4 Given $A\in\mathbb{C}^{N\times N},~{}b\in\mathbb{C}^{N\times
n_{0}},$ this algorithm computes $\varphi(A)b$ by the scaling and modified
squaring based on Taylor series.

1:Select optimal values of $m$ and $s.$

2:Compute $Y=A/s.$

3:Compute $b_{1}=\sum\limits^{m}_{k=0}\frac{Y^{k}}{(1+k)!}b$ based on matrix-
vector products.

4:Compute $f=b_{1}.$

5:for $i=1:s-1$ do

6: Compute $b_{i+1}=\sum\limits^{m+1}_{k=0}\frac{Y^{k}}{k!}b_{i}$ based on
matrix-vector products.

7: Compute $f=f+b_{i}.$

8:end for

9:Compute $f=\frac{1}{s}f.$

10: $f$

The procedure described above mentions two key parameters: the degree $m$ of
the matrix polynomial $T_{m}(Y)$ and the scaling parameter $s$. We use the
backward error analysis combined with the computational cost to choose an
optimal parameters $m$ and $s$. The backward error analysis of the method is
exactly the same as the above section. The only difference is the form of
their scaling coefficients. The former is $2^{-s}$ and the latter is
$\frac{1}{s}.$ The relative backward error of the method satisfies

$\frac{\|\Delta A\|}{\|A\|}\leq\tilde{h}_{m+2}(\frac{1}{s}\alpha_{p}(A)),$
(29)

where $\tilde{h}_{m+2}(x)$ and $\alpha_{p}(A)$ are defined exactly as in (15).
Given a tolerance Tol and integer $m,$ the parameter $s$ is chosen so that
$s^{-1}\alpha_{p}(A)\leq\theta_{m},$ i.e.,

$s\geq\lceil\alpha_{p}(A)/\theta_{m}\rceil.$ (30)

And the cost of the algorithm in matrix-vector products is

$C_{m}=(m+1)\lceil\alpha_{p}(A)/\theta_{m}\rceil-1.$ (31)

Let $p_{max}$ denote the largest positive integer $p$ such that $p(p-1)\leq
m_{max}+2.$ Then the optimal cost is

$C_{m^{*}}=\min\\{(m+1)\lceil\alpha_{p}(A)/\theta_{m}\rceil-1:~{}~{}2\leq
p\leq p_{max},p(p-1)-2\leq m\leq m_{max}\\},$ (32)

where $m^{*}$ denotes the smallest value of $m$ at which the minimum is
attained. The optimal scaling parameter is $s=C_{m^{*}}/m^{*}.$

The cost of computing $\alpha_{p}(A)$ for $p=2:p_{max}$ is approximately
$2lp_{max}(p_{max}+3),~{}~{}l=1~{}\text{or}~{}2$. If the cost $C_{m_{max}}$
matrix-vector products of evaluating $\varphi(A)b$ with $m$ determined by
using $\|A\|_{1}$ in place of $\alpha_{p}(A)$ in (31) is no larger than the
cost of computing the $\alpha_{p}(A)$, i.e.

$\|A\|_{1}<=[4*p_{max}*(p_{max}+3)+1]*[\theta_{m_{max}}/(m_{max}+1)].$ (33)

Then we should certainly use $\|A\|_{1}$ in place of the $\alpha_{p}(A)$. The
details of the method is summarized in Algorithms 4 and 5.

Algorithm 5 Given $A\in\mathbb{C}^{N\times N},~{}b\in\mathbb{C}^{N\times
n_{0}},$ $p_{max}$ and $m_{max}.$ this algorithm computes the parameters $m$
and $n$ based on the number of matrix-vector products.

1:$M=1:m_{max}.$

2:$d_{1}=\|A\|_{1}.$

3:if $d_{1}<=\theta_{m_{max}}*(4*p_{max}*(p_{max}+3)+1)/m_{max}$ then

4: $m=\arg\min\\{(m+1)\lceil d_{1}/\theta_{m}\rceil-1,~{}1\leq m\leq
m_{max}\\}.$

5: $n=\lceil d_{1}/\theta_{m}\rceil,$ return .

6:end if

7:for $p=2:p_{max}+1$ do

8: $c=\text{normest}(A,p).$

9: $d_{p}=c^{1/p}.$

10:end for

11:$\alpha_{1}=d_{1}.$

12:for $p=2:p_{max}$ do

13: $\alpha_{p}=\max(d_{p},d_{p+1}).$

14:end for

15:$[m,p]=\arg\min\\{(m+1)\lceil\alpha_{p}/\theta_{m}\rceil-1:2\leq p\leq
p_{max},p(p-1)-2\leq m\leq m_{max}\\}.$

16:$s=\max(\lceil\alpha_{p}/\theta_{m}\rceil,1).$

17: $m,s.$

##  4 Numerical experiments

In this section we perform two numerical experiments to test the performance
of the approach that has been presented in the previous sections. All tests
are performed under Windows 10 and MATLAB R2018b running on a laptop with an
Intel Core i7 processor with 1.8 GHz and RAM 8 GB.

We use Algorithm 1 in combination with Algorithm 2 and Algorithm 3 to evaluate
$\varphi(A),$ and Algorithm 4 combined with Algorithm 5 to compute
$\varphi(A)b.$ The three combined algorithms are denoted as phitay1, phitay2
and phimv, respectively.

######  Experiment 1.

In this experiment we compare algorithms phitay1 and phitay2 with existing
MATLAB routine phipade13 from [26]. The function phipade13 employs scaling and
modified squaring method based on [13/13] padé approximation to evaluate
$\varphi(A).$ We use a total of 201 matrices divided into two sets to test
these algorithms. The two test sets are described as follows:

$\bullet$ The first test set contains 62 $8\times 8$ test matrices as in [9]
and [24, sec. 4.1]. The first 48 matrices are obtained from the subroutine
matrix111The subroutine matrix can generate fifty-two matrices. Matrices 17,
42, 44, 43 are excluded the scope of the test as the first three overflow in
double precision and the last is repeated as matrix 49. in the matrix
computation toolbox [11]. The other fourteen test matrices of dimension $2-20$
come from [5, Ex. 2], [7, Ex. 3.10], [15, p. 655], [19, p. 370], [27, Test
Cases 1–4].

$\bullet$ The second test set is essentially the same tests as in [6], which
consists of 139 matrices of dimension $n=128.$ The first 39 matrices are
obtained from MATLAB routine matrix of the Matrix Computation Toolbox [11].
The remaining 100 matrices are generated randomly, half of which are
diagonalizable and half non diagonalizable matrices.

In this implementation, we evaluated the relative errors in the 1-norm of the
computed solutions $Y$, i.e.,

$Error=\frac{\|Y-\varphi(A)\|_{1}}{\|\varphi(A)\|_{1}}.$ (34)

The ”exact” $\varphi(A)$ is computed using MATLAB build-in function expm of
[9, 1] by evaluating the augmented matrix exponential (7) at 100-digit
precision using MATLAB’s Symbolic Math Toolbox.

a. Normwise relative errors

b. Performance of errors
c. Ratio of execution times d. Performance of execution times

Figure 2: Results for test matrix set 1 for Experiment 1.

In Figs. 2 and 3, we present the relative errors, the performances of relative
errors, the ratio of execution times and the performances of execution times
for each test. Figs. 2(a) and 3(a) display the relative error of the
algorithms in our test sets, sorted by decreasing relative condition number of
$\varphi(A)$ at $A$. The solid black line represents the unit roundoff
multiplied by the relative condition number, which is estimated by the MATLAB
routine funm condest1 in the Matrix Function Toolbox [11]. Figs. 2(b) and 3(b)
show the performance profiles of the three solvers on the same error data. For
a given $\alpha,$ the corresponding value of $p$ on each performance curve is
the probability that the algorithm has a relative error lower than or equal to
$\alpha$ times the smallest error over all the methods involved [8]. The
results show that the methods based on Taylor series are more accurate than
the implementation based on Padé series, and phitay2 is slightly more accurate
than phitay1. Figs. 2(c) and 3(c) show the the ratio of execution times of the
three solvers with respect to phipade13. The performance on the execution
times of three methods is compared in Figs. 2(d) and 3(d). We notice that
phitay1 and phitay2 have lower execution times than phipade13, and the
execution time of phitay1 is slightly lower than the execution time of phitay2
on test set 1 but the opposite is true on test set 2. This can be attributed
to the choice of the key parameters $m$ and $s$ in the methods. The phitay2
uses a minimum amount of computational costs to determine the optimal
parameters by evaluating exactly the 1-norm of $\|A^{k}\|_{1}^{1/k}$ for a few
values of $k$ using a matrix norm estimator. Although this requires some extra
calculations in estimating the norm of matrix power, the computational
advantages will gain as the dimension of the matrix increases.

a. Normwise relative errors

b. Performrnce of errors
c. Ratio of execution times d. Performance of execution times

Figure 3: Experimental results for test matrix set 2 for Experiment 1.

######  Experiment 2.

This experiment uses the same tests as [20]. There are four different sparse
matrices test matrices. The matrices details are

$\bullet$ The first matrix orani678 is an unsymmetric sparse matrix of order
$N=2,529$ with $nnz=90,158$ nonzero elements and its 1-norm is 1.04e+03.

$\bullet$ The second matrix bcspwr10 is a symmetric Hermitian sparse matrix of
order $N=5,300$ with $nnz=21,842$ nonzero elements and its 1-norm is 14.

$\bullet$ The third matrix gr 30 30 is an symmetric sparse matrix of order
$N=900$ with $nnz=7,744$ nonzero elements and its 1-norm is 16.

$\bullet$ The fourth matrix helm2d03 is a sparse matrix of order $N=392,257$
has $nnz=2,741,935$ nonzero elements and its 1-norm is 10.72.

We use our algorithm phimv with other two popular MATLAB routines phiv of [25]
and phipm of [20] to evaluate $\varphi(tA)b$ and
$\varphi_{0}(tA)b_{0}+t\varphi_{1}(tA)b_{1}$ with $t=10$ for the first test
matrix and $t=2$ for the other three, respectively. As in [20], we choose the
vectors $b=b_{0}=b_{1}=[1,1,...,1,1]^{T}$ except for the second test matrix
with $b=[1,0,...,0,1]^{T}.$ The MATLAB routines phiv and phipm are run with
their default parameters and the uniform convergence tolerance Tol=eps (eps is
the unit roundoff) in our experiments.

In this tests, we assess the accuracy of the computed solution $y$ by the
relative errors

$Error=\frac{\|y-y_{exa}\|_{2}}{\|y_{exa}\|_{2}}.$ (35)

where $y_{exa}$ is a reference solution obtained by computing the action of
the augmented matrix exponential using MATLAB routine expmv of [2]. We measure
the average ratio of execution times ($t_{ratio}$) of each of the codes
relative to phikmv by running the comparisons 100 times. Tables 3 and 4 show
the numerical results. All three methods deliver almost the same accuracy, but
phikmv performs to be the fastest.

Table 3: Comparisons of the average speedup of phiv, phip and phipm with respect to phimv and the relative errors for solving $\varphi(A)b$. method | orani678 | bcspwr10 | gr 30 30 | helm2d03
---|---|---|---|---
| error | $t_{ratio}$ | error | $t_{ratio}$ | error | $t_{ratio}$ | error | $t_{ratio}$
phiv | 1.1966e-15 | 1.38 | 4.8688e-16 | 45.90 | 1.8016e-14 | 9.92 | 1.0842e-13 | 20.98
phipm | 1.7612e-15 | 0.79 | 6.6655e-16 | 7.4 | 3.7391e-15 | 5.54 | 1.3550e-13 | 1.59
phimv | 1.1682e-15 | 1 | 3.6051e-16 | 1 | 1.2622e-15 | 1 | 6.2692e-14 | 1
Table 4: Comparisons of the average speedup of phiv, phip and phipm with respect to phimv and the relative errors for solving $\varphi_{0}(tA)b_{0}+t\varphi_{1}(tA)b_{1}$. method | orani678 | bcspwr10 | gr 30 30 | helm2d03
---|---|---|---|---
| error | $t_{ratio}$ | error | $t_{ratio}$ | error | $t_{ratio}$ | error | $t_{ratio}$
phiv | 4.7206e-16 | 1.06 | 7.4407e-16 | 43.63 | 2.1790e-15 | 15.58 | 1.5239e-14 | 20.73
phipm | 1.0408e-15 | 0.63 | 1.2216e-15 | 8.73 | 7.4383e-15 | 5.58 | 8.7267e-15 | 1.77
phimv | 1.8024e-15 | 1 | 7.6561e-16 | 1 | 8.7257e-16 | 1 | 8.7682e-15 | 1

##  5 Conclusion

The computation of $\varphi$-functions can lead to a large computational
burden of exponential integrators. In this work three accurate algorithms
phitay1, phitay2 and phimv have been developed to compute the first
$\varphi$-function and its action on a vector. The first two are used for
solving $\varphi(A)$ and the last one is used for $\varphi(A)b$. These
algorithms employ the scaling and modified squaring procedure based on the
truncated Taylor series of the $\varphi$-function and are backward stable in
exact arithmetic. For phitay1 and phitay2, the optimal Horner and Paterson-
Stockmeyer’s technique has been applied to reduce the computational cost. The
main difference of both is the estimation of matrix powers
$\|A^{k}\|_{1}^{1/k}$ for a few values of $k$ and phitay2 allows to determine
the optimal values of scaling and the degree of the Taylor approximation by
the minimum amount of computational costs. The phimv takes the similar
approach as phitay2 to determine the key parameters. The computational costs
mostly focused on computing matrix-vector products which is especially well-
suited large sparse matrix. Numerical comparisons with other state-of-the-art
MATLAB routines illustrate that the methods proposed are efficient and
reliable. In the future we hope to further generalize these methods to the
general exponential related functions and their linear combination.

## Acknowledgements

This work was supported in part by the Jilin Scientific and Technological
Development Program (Grant Nos. 20200201276JC and 20180101224JC) and the
Natural Science Foundation of Jilin Province (Grant No. 20200822KJ), and the
Scientific Startup Foundation for Doctors of Changchun Normal University
(Grant No. 002006059).

## References

  * [1] A.H. Al-Mohy, N.J. Higham, _A new scaling and squaring algorithm for the matrix exponential,_ SIAM J. Matrix Anal. Appl., 31 (3) (2009), 970-989.
  * [2] A. Al-Mohy and N. Higham, _Computing the action of the matrix exponential, with an application to exponential integrators,_ SIAM J. Sci. Comput., 33 (2011), pp. 488-511.
  * [3] G. Beylkin, J.M. Keiser, L. Vozovoi, _A new class of time discretization schemes for the solution of nonlinear PDEs,_ J. Comput. Phys., 147 (1998), 362–387.
  * [4] M. Caliari, M. Vianello, L. Bergamaschi, _Interpolating discrete advection-diffusion propagators at Leja sequences_ , J. Comput. Appl. Math. 172 (1)(2004), pp.79-99.
  * [5] I. Davies and N. J. Higham, _A Schur-Parlett algorithm for computing matrix functions,_ SIAM J. Matrix Anal. Appl., 25 (2003), pp. 464-485.
  * [6] E. Defez, J. Ibáñez, J. Sastre, J. Peinado and P. Alonso, _A new efficient and accurate spline algorithm for the matrix exponential computation,_ J. Comput. Appl. Math., 337 (2018), pp. 354-365.
  * [7] L. Dieci and A. Papini, _Padé approximation for the exponential of a block triangular matrix,_ Linear Algebra Appl., 308 (2000), pp. 183-202.
  * [8] E. D. Dolan and J. J. Moré _Benchmarking optimization software with performance profiles,_ Math. Program., 91 (2002), pp. 201-213.
  * [9] N. J. Higham, _The scaling and squaring method for the matrix exponential revisited,_ SIAM J. Matrix Anal. Appl., 26 (2005), pp. 1179-1193.
  * [10] N. J. Higham, _Functions of matrices: theory and computation,_ SIAM, Philadelphia, 2008.
  * [11] N. J. Higham, _The Matrix Computation Toolbox,_ http://www.ma.man.ac.uk/ higham/mctoolbox.
  * [12] M. Hochbruck, C. Lubich, H. Selhofer, _Exponential integrators for large systems of differential equations._ SIAM J. Sci. Comput. 19 (1998), pp. 1552-1574.
  * [13] M. Hochbruck,and A. Ostermann, _Exponential Integrators,_ Acta Numer., 19 (2010), pp. 209-286.
  * [14] A.K. Kassam and L.N. Trefethen, _Fourth-order time stepping for stiff PDEs,_ SIAM J. Sci. Comput., 26 (2005), pp. 1214-1233.
  * [15] C. S. Kenney and A. J. Laub, _A Schur-Fréchet algorithm for computing the logarithm and exponential of a matrix,_ SIAM J. Matrix Anal. Appl., 19 (1998), pp. 640-663.
  * [16] Y.Y. Lu, _Computing a matrix function for exponential integrators,_ J. Comput. Appl. Math. 161 (1) (2003), pp. 203–216.
  * [17] B.V. Minchev and W.M. Wright, _A review of exponential integrators for first order semi-linear problems,_ Tech. report 2/05, Department of Mathematics, NTNU, 2005.
  * [18] C. Moler, C.V. Loan, _Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,_ SIAM Review, 45 (2003), pp. 3-49.
  * [19] I. Najfeld and T. F. Havel, _Derivatives of the matrix exponential and their computation,_ Adv. in Appl. Math., 16 (1995), pp. 321-375.
  * [20] J. Niesen, W. Wright, _Algorithm 919: A Krylov subspace algorithm for evaluating the phi- functions appearing in exponential integrators_ , ACM Trans. Math. Software, 38 (3) (2012), Article 22.
  * [21] M.S. Paterson, L.J. Stockmeyer, _On the number of nonscalar multiplications necessary to evaluate polynomials,_ SIAM J. Comput. 2 (1) (1973), pp.60-66.
  * [22] Y. Saad, _Analysis of some Krylov subspace approximations to the matrix exponential operator,_ SIAM J. Numer. Anal., 29 (1992), pp. 209-228.
  * [23] J. Sastre, J. Ibáñez, E. Defez, P. Ruiz, _Efficient orthogonal matrix polynomial based method for computing matrix exponential,_ Appl. Math. Comput, 217 (14) (2011), pp. 6451-6463.
  * [24] J. Sastre, J. Ibáñez, E. Defez, and P. Ruiz, _New Scaling-Squaring Taylor Algorithms for Computing the Matrix Exponential,_ SIAM J. Sci. Comput., 37 (1) (2015), pp. 439-455.
  * [25] R.B. Sidje, _Expokit: A software package for computing matrix exponentials,_ ACM Trans. Math. Softw., 24 (1998), pp. 130-156.
  * [26] B. Skaflestad and W.M. Wright, _The scaling and modified squaring method for matrix functions related to the exponential,_ Applied Numerical Mathematics, 59 (2009), pp. 783-799.
  * [27] R. C. Ward, _Numerical computation of the matrix exponential with accuracy estimate,_ SIAM J. Numer. Anal., 14 (1977), pp. 600-610.