where the pole in the last line arises from performing the energy integral. We
are left with angular integrals only, for which we write

$\displaystyle\int\mathrm{d}\Omega_{4}^{(3-2\epsilon)}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{\hat{q}_{4}})({\hat{q}_{4}}\cdot{q_{j}})}={}$
$\displaystyle\int\mathrm{d}\Omega_{4}^{(1-2\epsilon)}\int\mathrm{d}(\cos\theta)\mathrm{d}\phi(\sin\theta\sin\phi)^{-2\epsilon}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{\hat{q}_{4}})({\hat{q}_{4}}\cdot{q_{j}})}$
$\displaystyle={}$
$\displaystyle\frac{\Gamma(1-\epsilon)}{(4\pi)^{\epsilon}\Gamma(1-2\epsilon)}(2\pi)I(q_{i},q_{j})\,.$
(D.13)

Here, we introduce the auxiliary function

$\displaystyle I(q_{i},q_{j})={}$
$\displaystyle\int\mathrm{d}(\cos\theta)\int\frac{\mathrm{d}\phi}{\pi}(\sin\theta\sin\phi)^{-2\epsilon}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{\hat{q}_{4}})({\hat{q}_{4}}\cdot{q_{j}})}\,.$
(D.14)

We write it as a Laurent series in $\epsilon$, i.e.

$\displaystyle I(q_{i},q_{j})={}$
$\displaystyle\sum_{k=-1}^{2}\epsilon^{k}I^{(k)}(q_{i},q_{j})+\mathcal{O}(\epsilon^{3})\,,$
(D.15)

where the coefficients $I^{(k)}(q_{i},q_{j})$ depend on the momenta $q_{i}$
and $q_{j}$, in particular on whether they are massless or massive. All
necessary coefficients for the single-soft integrated subtraction terms,
except for the $\operatorname{\mathcal{O}}\left(\epsilon^{2}\right)$ terms,
have been obtained in Ref. [92], see Appendix A therein. For completeness, we
collect below the formulae for those $I^{(k)}(q_{i},q_{j})$ which are relevant
to our calculation.

We write the full integrated single-soft eikonal factor as

$\displaystyle\mathcal{S}^{(0)}_{ij,\mathrm{int}}={}$
$\displaystyle-\,\frac{1}{(4\pi)^{2}}\,\frac{1}{\epsilon}\,\frac{e^{\epsilon\gamma_{\mathrm{E}}}\Gamma(1-\epsilon)}{\Gamma(1-2\epsilon)}\left(\frac{\mu_{R}^{2}}{4E_{\mathrm{max}}^{2}}\right)^{\epsilon}I(q_{i},q_{j})\,.$
(D.16)

We have to distinguish the case where the emitters are two massive particles,
which we denote by a subscript $\mathrm{MM}$, and the case where emitter $i$
is massive and emitter $j$ is massless, which we denote by $\mathrm{M0}$. For
our calculation, we need the coefficients $I_{\mathrm{MM}}^{(k)}(q_{i},q_{j})$
and $I_{\mathrm{M0}}^{(k)}(q_{i},q_{j})$ for $k\in\\{-1,0,1\\}$ when the
directions of the momenta $q_{i}$ and $q_{j}$ are arbitrary. Moreover, we need
$I_{\mathrm{MM}}^{(k)}(q_{i},q_{j})$ for $k\in\\{-1,0,1,2\\}$ in two special
cases: for the case where the two momenta are equal, $q_{i}=q_{j}$ and for the
case where $q_{i}$ and $q_{j}$ are in a back-to-back configuration.

#### Two massive emitters

The reported formulae correspond to the result outlined in Eqs. (A.41) to
(A.51) of Ref. [92]. We consider two time-like momenta
$q_{i}^{\,2}=m_{i}^{\,2}$ and $q_{j}^{\,2}=m_{j}^{\,2}$. With the definition
of $v_{ij}$ from Eq. B.3 and the notation

$\displaystyle\vec{u}={}\frac{\vec{q}_{i}}{E_{i}}\,,$ and
$\displaystyle\vec{w}={}\frac{\vec{q}_{j}}{E_{j}}\,,$ (D.17)

we introduce the shorthand notations

$\displaystyle A^{2}={}$
$\displaystyle({\vec{u}}\cdot{\vec{u}})+({\vec{w}}\cdot{\vec{w}})-2({\vec{u}}\cdot{\vec{w}})\,,$
$\displaystyle X_{1}={}$
$\displaystyle({\vec{u}}\cdot{\vec{u}})-({\vec{u}}\cdot{\vec{w}})\,,$
$\displaystyle B^{2}={}$
$\displaystyle({\vec{u}}\cdot{\vec{u}})({\vec{w}}\cdot{\vec{w}})-({\vec{u}}\cdot{\vec{w}})^{2}\,,$
$\displaystyle X_{2}={}$
$\displaystyle({\vec{w}}\cdot{\vec{w}})-({\vec{u}}\cdot{\vec{w}})\,.$ (D.18)

Furthermore, we need the following arguments

$\displaystyle z_{+}={}$ $\displaystyle A+\sqrt{A^{2}-B^{2}}\,,$
$\displaystyle z_{1}={}$ $\displaystyle\sqrt{X_{1}^{\,2}+B^{2}}-X_{1}\,,$
$\displaystyle z_{-}={}$ $\displaystyle A-\sqrt{A^{2}-B^{2}}\,,$
$\displaystyle z_{2}={}$ $\displaystyle\sqrt{X_{2}^{\,2}+B^{2}}+X_{2}\,,$
(D.19)

which will be used in the function

$\displaystyle K(z)={}$
$\displaystyle-2\mathrm{Li}_{2}\left({\frac{2z_{-}(z_{+}-z)}{(z_{+}-z_{-})(z_{-}+z)}}\right)-2\mathrm{Li}_{2}\left({\frac{-2z_{+}(z_{-}+z)}{(z_{+}-z_{-})(z_{+}-z)}}\right)$
$\displaystyle-\frac{1}{2}\log^{2}\left({\frac{(z-z_{-})(z_{+}-z)}{(z+z_{-})(z_{+}+z)}}\right)\,.$
(D.20)

With these abbreviations, the coefficients of the massive-massive angular
integral read

$\displaystyle I_{\mathrm{MM}}^{(-1)}(q_{i},q_{j})={}$ $\displaystyle 0\,,$
(D.21) $\displaystyle I_{\mathrm{MM}}^{(0)}(q_{i},q_{j})={}$
$\displaystyle\frac{1}{v_{ij}}\log\left({\frac{1+v_{ij}}{1-v_{ij}}}\right)\,,$
(D.22) $\displaystyle I_{\mathrm{MM}}^{(1)}(q_{i},q_{j})={}$
$\displaystyle\frac{1-({\vec{u}}\cdot{\vec{w}})}{\sqrt{A^{2}-B^{2}}}\left({K(z_{2})-K(z_{1})}\right)\,.$
(D.23)

#### One massive and one massless emitter

The reported formulae correspond to the result outlined in Eqs. (A.22) to
(A.24) of Ref. [92]. We consider one time-like momentum
$q_{i}^{\,2}=m_{i}{\,{}^{2}}$ and one light-like momentum $q_{j}^{\,2}=0$. We
define the symbol

$\displaystyle\kappa={}\sqrt{1-\frac{m_{i}^{\,2}}{E_{i}^{\,2}}}\,.$ (D.24)

Then the coefficients of the massive-massless angular integral are

$\displaystyle I_{\mathrm{M0}}^{(-1)}(q_{i},q_{j})={}$ $\displaystyle-1\,,$
(D.25) $\displaystyle I_{\mathrm{M0}}^{(0)}(q_{i},q_{j})={}$
$\displaystyle\log\left({\frac{({\hat{q}_{i}}\cdot{\hat{q}_{j}})^{2}}{({\hat{q}_{i}}\cdot{\hat{q}_{i}})}}\right)\,,$
(D.26) $\displaystyle I_{\mathrm{M0}}^{(1)}(q_{i},q_{j})={}$
$\displaystyle-2\left\\{\frac{1}{4}\log^{2}\left({\frac{1-\kappa}{1+\kappa}}\right)+\log\left({\frac{({\hat{q}_{i}}\cdot{\hat{q}_{j}})}{1+\kappa}}\right)\log\left({\frac{({\hat{q}_{i}}\cdot{\hat{q}_{j}})}{1-\kappa}}\right)\right.$
$\displaystyle\left.\hskip
28.45274pt+\mathrm{Li}_{2}\left({1-\frac{({\hat{q}_{i}}\cdot{\hat{q}_{j}})}{1+\kappa}}\right)+\mathrm{Li}_{2}\left({1-\frac{({\hat{q}_{i}}\cdot{\hat{q}_{j}})}{1-\kappa}}\right)\right\\}\,.$
(D.27)

#### Two massive back-to-back emitters

We consider the special case of two massive emitters with the same mass,
$q_{i}^{\,2}=m^{2}$ and $q_{j}^{\,2}=m^{2}$, arranged in a back-to-back
configuration, $\vec{q}_{i}=-\vec{q}_{j}$. We then have $E_{i}=E_{j}=M_{H}/2$
and, therefore, $\kappa=\beta$. The expansion coefficients are given by

$\displaystyle I_{\mathrm{MM},\mathrm{b2b}}^{(-1)}(q_{i},q_{j})={}$
$\displaystyle 0\,,$ (D.28) $\displaystyle
I_{\mathrm{MM},\mathrm{b2b}}^{(0)}(q_{i},q_{j})={}$
$\displaystyle-\frac{1+\beta^{2}}{\beta}\log\left({\frac{1-\beta}{1+\beta}}\right)\,,$
(D.29) $\displaystyle I_{\mathrm{MM},\mathrm{b2b}}^{(1)}(q_{i},q_{j})={}$
$\displaystyle\frac{1+\beta^{2}}{\beta}\left\\{\mathrm{Li}_{2}\left({\frac{2\beta}{1+\beta}}\right)-\mathrm{Li}_{2}\left({-\frac{2\beta}{1-\beta}}\right)\right\\}\,,$
(D.30) $\displaystyle I_{\mathrm{MM},\mathrm{b2b}}^{(2)}(q_{i},q_{j})={}$
$\displaystyle\frac{1+\beta^{2}}{\beta}\left\\{-\frac{1}{3}\log^{3}\left({\frac{1-\beta}{1+\beta}}\right)+2\mathrm{Li}_{3}\left({\frac{2\beta}{1+\beta}}\right)-2\mathrm{Li}_{3}\left({-\frac{2\beta}{1-\beta}}\right)\right.$
$\displaystyle\hskip
28.45274pt\left.-\log\left({\frac{1-\beta}{1+\beta}}\right)\left[{\mathrm{Li}_{2}\left({\frac{2\beta}{1+\beta}}\right)+\mathrm{Li}_{2}\left({-\frac{2\beta}{1-\beta}}\right)}\right]\right\\}\,,$
(D.31)

where the subscript “$\mathrm{b2b}$” indicates a back-to-back configuration of
the emitters. Formulae for $I_{\mathrm{MM},\mathrm{b2b}}^{(k)}$ for
$k\in\\{-1,0,1\\}$ can be obtained from Eqs. D.21 to D.23 by taking the
relevant limit. The formula for $I_{\mathrm{MM},\mathrm{b2b}}^{(2)}$ was
calculated independently; this $\mathcal{O}(\epsilon^{2})$ term is needed in
the soft limit of the real-virtual contribution.

#### Self-correlated massive emitter

We consider the special case of a self-correlated massive emitter, $i=j$, for
a time-like momentum $q_{i}$ with $q_{i}^{\,2}=m^{2}$. Then we have

$\displaystyle I_{\mathrm{MM}}^{(-1)}(q_{i},q_{i})={}$ $\displaystyle 0\,,$
(D.32) $\displaystyle I_{\mathrm{MM}}^{(0)}(q_{i},q_{i})={}$ $\displaystyle
2\,,$ (D.33) $\displaystyle I_{\mathrm{MM}}^{(1)}(q_{i},q_{i})={}$
$\displaystyle-\frac{2}{\kappa}\log\left({\frac{1-\kappa}{1+\kappa}}\right)\,,$
(D.34) $\displaystyle I_{\mathrm{MM}}^{(2)}(q_{i},q_{i})={}$
$\displaystyle\frac{2}{\kappa}\left\\{\mathrm{Li}_{2}\left({\frac{2\kappa}{1+\kappa}}\right)-\mathrm{Li}_{2}\left({-\frac{2\kappa}{1-\kappa}}\right)\right\\}\,.$
(D.35)

The $\operatorname{\mathcal{O}}\left(\epsilon^{2}\right)$ term was calculated
independently.

###  D.3 Single-soft subtraction terms (one-loop)

We consider the soft limit of the one-loop $H\to{b\bar{b}{g}}$ amplitude which
is given in Section C.3. This factorisation involves both tree-level and one-
loop soft functions. The integration of the tree-level eikonal factors is
discussed in Section D.2, while in this section we focus on the integration of
those terms in the one-loop soft function of Section C.3 which contain
$\mathcal{R}^{(1)}_{ij,4}$.

We integrate the one-loop eikonal factors over the soft-gluon phase space,
i.e.

$\displaystyle\mathcal{R}^{(1)}_{\mathrm{int}}={}$
$\displaystyle-\int[\mathrm{d}{}q_{4}]\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)\mathcal{R}^{(1)}_{23,4}\,,$
(D.36)

where $\mathcal{S}^{(0)}_{ij,4}$ and $\mathcal{R}^{(1)}_{23,4}$ are defined in
Eqs. C.11 and C.13, respectively. Since only $\mathcal{R}^{(1)}_{23,4}$
appears, we restrict ourselves to the case of two emitters with the same non-
vanishing mass in a back-to-back configuration.

The soft-gluon phase space is parametrised using Section 2.4. Note that even
though the soft-gluon momentum factorises from the energy-momentum conserving
$\delta$-function, we restrict ourselves to the same integration region as
stated in Section 2.4; this is in accordance with the choice made for the
integrated tree-level eikonal factor. The only non-trivial integral is the
integration over the angle $\theta$ between the momentum of the soft gluon
$q_{4}$ and the $b$-quark momentum $q_{i}$. We reexpress this angle in terms
of the variable $\lambda=\cos\theta$ and use the symmetry under
$\lambda\to-\lambda$ to restrict the domain of integration to
$\lambda\in[0,1]$. We arrive at

$\displaystyle\int[\mathrm{d}{}q_{4}]$
$\displaystyle=\frac{\mu_{R}^{2\epsilon}S_{\epsilon}E_{\mathrm{max}}^{2-2\epsilon}}{(2\pi)^{3-2\epsilon}}\int\mathrm{d}\Omega_{4}^{(2-2\epsilon)}\int_{0}^{1}\mathrm{d}\xi\,\xi^{1-2\epsilon}\int_{0}^{1}\mathrm{d}\lambda\,(1-\lambda^{2})^{-\epsilon}\,.$
(D.37)

The energy dependence factorises from the integrand, which means that the
integral over $\mathrm{d}\xi$ can be solved trivially. The angular integral
$\mathrm{d}\Omega_{4}^{(2-2\epsilon)}$ is performed using Eq. 2.22.

The expressions for the expansion coefficients $R_{k}(q_{i},q_{j};q_{4})$ in
Ref. [90] are given in terms of the variables

$\displaystyle\alpha_{i}$
$\displaystyle=\frac{m_{i}^{2}}{2}\frac{({q_{j}}\cdot{q_{4}})}{({q_{i}}\cdot{q_{4}})({q_{i}}\cdot{q_{j}})}=\frac{(1-\beta^{2})(1+\beta\lambda)}{2(1+\beta^{2})(1-\beta\lambda)}\,,$
(D.38) $\displaystyle\alpha_{j}$
$\displaystyle=\frac{m_{j}^{2}}{2}\frac{({q_{i}}\cdot{q_{4}})}{({q_{j}}\cdot{q_{4}})({q_{i}}\cdot{q_{j}})}=\frac{(1-\beta^{2})(1-\beta\lambda)}{2(1+\beta^{2})(1+\beta\lambda)}\,,$
(D.39)

and contain at most classical polylogarithms with arguments composed of these
variables. We rewrite those special functions in terms of iterated integrals
of argument $\lambda$ over the alphabet

$\displaystyle\frac{\mathrm{d}\lambda}{\lambda},$
$\displaystyle\frac{\mathrm{d}\lambda}{1+\lambda},$
$\displaystyle\frac{\mathrm{d}\lambda}{1-\lambda},$
$\displaystyle\frac{\mathrm{d}\lambda}{1+\beta\lambda},$
$\displaystyle\frac{\mathrm{d}\lambda}{1-\beta\lambda}\,.$ (D.40)

As the rational coefficients in front of the iterated integrals also only
contain these letters, the integration over $\lambda$ can again be performed
in terms of iterated integrals over the same alphabet. These iterated
integrals are then evaluated at $1$ and we rewrite them in terms of harmonic
polylogarithms [66] of argument $\beta$ with up to weight four. All
manipulations of the iterated integrals are performed with the Mathematica
package HarmonicSums [93, 94, 95, 66, 96, 97, 98, 99, 100, 101, 102].

Note that the term $\mathcal{R}^{(1)}_{\mathrm{int}}$ is again a Laurent
series in $\epsilon$. Therefore, we write

$\displaystyle\mathcal{R}^{(1)}_{\mathrm{int}}={}$
$\displaystyle\frac{4C_{A}}{(4\pi)^{2}}\left({\frac{\mu_{R}^{2}}{4E_{\mathrm{max}}^{2}}}\right)^{2\epsilon}\sum_{k=-3}^{0}\mathcal{R}^{(1,k)}_{\mathrm{int}}(\beta)\,\epsilon^{k}\,.$
(D.41)

The expansion coefficients $\mathcal{R}^{(1,k)}_{\mathrm{int}}$ for
$m_{b}=4.78~{}\mathrm{GeV}$ and $M_{H}=125.09~{}\mathrm{GeV}$ evaluate to

$\displaystyle\mathcal{R}^{(1,-3)}_{\mathrm{int}}(\beta)={}$ $\displaystyle
5.52628705137596\,,$
$\displaystyle\mathcal{R}^{(1,-2)}_{\mathrm{int}}(\beta)={}$ $\displaystyle
35.3923534452863\,,$
$\displaystyle\mathcal{R}^{(1,-1)}_{\mathrm{int}}(\beta)={}$ $\displaystyle
111.992677970445\,,$
$\displaystyle\mathcal{R}^{(1,0)}_{\mathrm{int}}(\beta)={}$ $\displaystyle
245.654621810082\,.$ (D.42)

The analytical expressions can be found in computer-readable form in an
ancillary file supplied with this article.

###  D.4 Double-soft subtraction terms

For the integrated subtraction terms of the double-real contribution we have
to integrate the double-soft function, discussed in Section C.4, over the
unresolved phase space in $d$ dimensions. We define the integrated double-soft
function as

$\displaystyle\mathrm{D}{}\mathrm{Soft}^{(0)}_{ij,\mathrm{int}}(q_{2},q_{3})={}$
$\displaystyle\int[\mathrm{d}{}q_{4}][\mathrm{d}{}q_{5}]\,\mathrm{D}{}\mathrm{Soft}^{(0)}_{ij}(q_{2},q_{3};q_{4},q_{5})\,.$
(D.43)

We use the phase-space parametrisation outlined in Section 2.4, taking into
account that the gluon momenta decouple from the overall energy-momentum
conserving $\delta$-function. Note that the decoupling of soft particles means
that, in principle, their energies are unbounded unless we introduce some
constraint by hand. For simplicity, we choose to keep the bound that
explicitly appears in the formulation of the phase-space measure, i.e. we keep
$E_{45,\mathrm{max}}$ unchanged while the
$\mathrm{d}\xi_{1}\,\mathrm{d}\xi_{2}$ integration region still covers the
unit square.

It is particularly useful to keep the Born configuration fixed so that the
$b$-quark points along the $\hat{z}$-axis, i.e.

$\displaystyle
q_{2}={}\tfrac{1}{2}M_{H}\left({\hat{t}^{\mu}+\beta\hat{z}^{\mu}}\right)$ and
$\displaystyle
q_{3}={}\tfrac{1}{2}M_{H}\left({\hat{t}^{\mu}-\beta\hat{z}^{\mu}}\right)\,,$
(D.44)

where $\hat{t}^{\mu}$ and $\hat{z}^{\mu}$ are unit vectors along time and
$\hat{z}$ axes, respectively. After expressing the integrand using the phase-
space parametrisation, the $\xi_{1}$ dependence factorises and can be
integrated analytically. The remaining integrals are performed numerically.

The integrands of Eq. D.43 are still divergent in the strongly-ordered
($S_{5}$) and collinear ($C_{45}$) limits. These divergences can be handled
using the endpoint subtraction method, as discussed in Section 2.2. We write

$\displaystyle S_{45}={}$
$\displaystyle(I{}-S_{5})(I{}-C_{45})S_{45}+(I{}-S_{5})C_{45}S_{45}+S_{5}(I{}-C_{45})S_{45}+S_{5}C_{45}S_{45}\,.$
(D.45)

Note that the regularisation in the strongly-ordered limit applies only to the
$gg$ emission case of Section C.4. The relevant subtraction terms in Eq. D.45
are constructed using the factorisation formulae from Appendix C.

| $k=-3$ | $k=-2$ | $k=-1$ | $k=0$
---|---|---|---|---
$C^{(k)}_{gg,C_{A}C_{F}}$ | $-22.105148(3)$ | $-120.8071(1)$ | $-337.441(1)$ | $-613.869(4)$
$C^{(k)}_{gg,C_{F}^{2}}$ | – | $+244.3194(4)$ | $+915.818(2)$ | $+984.741(6)$
$C^{(k)}_{q\bar{q},C_{F}T_{F}}$ | – | $+7.368385(3)$ | $+39.9006(1)$ | $+124.966(1)$
Table 6: Coefficients of the $\epsilon$ expansion of the integrated double-
soft function for $gg$ channel (first two rows), and the $q\bar{q}$ channel
(the last row). We use $m_{b}=4.78~{}\mathrm{GeV}$ and
$M_{H}=125.09~{}\mathrm{GeV}$.

The integrated double-soft function can be written as a Laurent series, i.e.

$\displaystyle\mathrm{D}{}\mathrm{Soft}^{(0)}_{ij,\mathrm{int}}(q_{2},q_{3})={}$
$\displaystyle\frac{1}{(4\pi)^{4}}\,\left({\frac{\mu_{R}^{2}}{E_{\mathrm{max}}^{2}}}\right)^{2\epsilon}\sum_{k=-3}^{0}C^{(k)}_{ij}\epsilon^{k}\,.$
(D.46)

Note that the highest pole that occurs in the $gg$ channel is $\epsilon^{-3}$.
This arises from taking all singular limits in the integrand, i.e. the
strongly-ordered, collinear and the double-soft limits. In the $q\bar{q}{}$
case we only have an $\epsilon^{-2}$ pole since the strongly-ordered soft
limit, $S_{5}$, is regular. Furthermore, we split the $C^{(k)}_{ij}$
coefficients into separate colour structures,

$\displaystyle C^{(k)}_{gg}={}$ $\displaystyle
C_{A}C_{F}\,C^{(k)}_{gg,C_{A}C_{F}}+C_{F}^{\,2}\,C^{(k)}_{gg,C_{F}^{2}}\,,$
(D.47) $\displaystyle C^{(k)}_{q\bar{q}}={}$ $\displaystyle
C_{F}T_{F}\,C^{(k)}_{q\bar{q},C_{F}T_{F}}\,.$ (D.48)

We present numerical values for the integrated double-soft function
coefficients, evaluated for $m_{b}=4.78~{}\mathrm{GeV}$ and
$M_{H}=125.09~{}\mathrm{GeV}$, in Table 6.

## References

  * [1] ATLAS, CMS Collaboration, G. Aad et al., “Combined Measurement of the Higgs Boson Mass in $pp$ Collisions at $\sqrt{s}=7$ and 8 TeV with the ATLAS and CMS Experiments,” Phys. Rev. Lett. 114 (2015) 191803, arXiv:1503.07589 [hep-ex].
  * [2] ATLAS Collaboration, M. Aaboud et al., “Observation of $H\rightarrow b\bar{b}$ decays and $VH$ production with the ATLAS detector,” Phys. Lett. B786 (2018) 59–86, arXiv:1808.08238 [hep-ex].
  * [3] CMS Collaboration, A. M. Sirunyan et al., “Observation of Higgs boson decay to bottom quarks,” Phys. Rev. Lett. 121 no. 12, (2018) 121801, arXiv:1808.08242 [hep-ex].
  * [4] E. Braaten and J. P. Leveille, “Higgs Boson Decay and the Running Mass,” Phys. Rev. D22 (1980) 715.
  * [5] N. Sakai, “Perturbative QCD Corrections to the Hadronic Decay Width of the Higgs Boson,” Phys. Rev. D22 (1980) 2220.
  * [6] P. Janot, “First Order QED and QCD Radiative Corrections to Higgs Decay Into Massive Fermions,” Phys. Lett. B223 (1989) 110–118.
  * [7] M. Drees and K.-i. Hikasa, “Note on QCD corrections to hadronic Higgs decay,” Phys. Lett. B240 (1990) 455. [Erratum: Phys. Lett. B262 (1991) 497].
  * [8] A. L. Kataev and V. T. Kim, “The Effects of the massless $O(\alpha_{s}^{2})$, $O(\alpha\alpha_{s})$, $O(\alpha^{2})$ QCD and QED corrections and of the massive contributions to $\Gamma(H^{0}\to b\bar{b})$,” arXiv:hep-ph/9304282 [hep-ph].
  * [9] P. A. Baikov, K. G. Chetyrkin, and J. H. Kühn, “Scalar correlator at $O(\alpha_{s}^{4})$, Higgs decay into $b$-quarks and bounds on the light quark masses,” Phys. Rev. Lett. 96 (2006) 012003, arXiv:hep-ph/0511063 [hep-ph].
  * [10] R. Harlander and M. Steinhauser, “Higgs decay to top quarks at $O(\alpha_{s}^{2})$,” Phys. Rev. D56 (1997) 3980–3990, arXiv:hep-ph/9704436 [hep-ph].
  * [11] S. A. Larin, T. van Ritbergen, and J. A. M. Vermaseren, “The Large top quark mass expansion for Higgs boson decays into bottom quarks and into gluons,” Phys. Lett. B362 (1995) 134–140, arXiv:hep-ph/9506465 [hep-ph].
  * [12] K. G. Chetyrkin and A. Kwiatkowski, “Second order QCD corrections to scalar and pseudoscalar Higgs decays into massive bottom quarks,” Nucl. Phys. B461 (1996) 3–18, arXiv:hep-ph/9505358 [hep-ph].
  * [13] A. Primo, G. Sasso, G. Somogyi, and F. Tramontano, “Exact Top Yukawa corrections to Higgs boson decay into bottom quarks,” Phys. Rev. D99 no. 5, (2019) 054013, arXiv:1812.07811 [hep-ph].
  * [14] E. Chaubey and S. Weinzierl, “Two-loop master integrals for the mixed QCD-electroweak corrections for $H\rightarrow b\bar{b}$ through a $Ht\bar{t}$-coupling,” JHEP 05 (2019) 185, arXiv:1904.00382 [hep-ph].
  * [15] C. Anastasiou, F. Herzog, and A. Lazopoulos, “The fully differential decay rate of a Higgs boson to bottom-quarks at NNLO in QCD,” JHEP 03 (2012) 035, arXiv:1110.2368 [hep-ph].
  * [16] V. Del Duca, C. Duhr, G. Somogyi, F. Tramontano, and Z. Trócsányi, “Higgs boson decay into b-quarks at NNLO accuracy,” JHEP 04 (2015) 036, arXiv:1501.07226 [hep-ph].
  * [17] F. Caola, G. Luisoni, K. Melnikov, and R. Röntsch, “NNLO QCD corrections to associated $WH$ production and $H\to b\bar{b}$ decay,” Phys. Rev. D97 no. 7, (2018) 074022, arXiv:1712.06954 [hep-ph].
  * [18] R. Gauld, A. Gehrmann-De Ridder, E. W. N. Glover, A. Huss, and I. Majer, “Associated production of a Higgs boson decaying into bottom quarks and a weak vector boson decaying leptonically at NNLO in QCD,” arXiv:1907.05836 [hep-ph].
  * [19] R. Mondini, M. Schiavi, and C. Williams, “N3LO predictions for the decay of the Higgs boson to bottom quarks,” JHEP 06 (2019) 079, arXiv:1904.08960 [hep-ph].
  * [20] W. Bernreuther, L. Chen, and Z.-G. Si, “Differential decay rates of CP-even and CP-odd Higgs bosons to top and bottom quarks at NNLO QCD,” JHEP 07 (2018) 159, arXiv:1805.06658 [hep-ph].
  * [21] G. Ferrera, G. Somogyi, and F. Tramontano, “Associated production of a Higgs boson decaying into bottom quarks at the LHC in full NNLO QCD,” Phys. Lett. B780 (2018) 346–351, arXiv:1705.10304 [hep-ph].
  * [22] W. Astill, W. Bizoń, E. Re, and G. Zanderighi, “NNLOPS accurate associated HZ production with $H\to b\overline{b}$ decay at NLO,” JHEP 11 (2018) 157, arXiv:1804.08141 [hep-ph].
  * [23] S. Alioli, A. Broggio, S. Kallweit, M. A. Lim, and L. Rottoli, “Higgsstrahlung at NNLL${}^{\prime}+$NNLO Matched to Parton Showers in GENEVA,” arXiv:1909.02026 [hep-ph].
  * [24] F. Granata, J. M. Lindert, C. Oleari, and S. Pozzorini, “NLO QCD+EW predictions for HV and HV +jet production including parton-shower effects,” JHEP 09 (2017) 012, arXiv:1706.03522 [hep-ph].
  * [25] F. Caola, K. Melnikov, and R. Röntsch, “Nested soft-collinear subtractions in NNLO QCD computations,” Eur. Phys. J. C77 no. 4, (2017) 248, arXiv:1702.01352 [hepph].
  * [26] F. Caola, K. Melnikov, and R. Röntsch, “Analytic results for color-singlet production at NNLO QCD with the nested soft-collinear subtraction scheme,” arXiv:1902.02081 [hep-ph].
  * [27] F. Caola, K. Melnikov, and R. Röntsch, “Analytic results for decays of color singlets to $gg$ and $q\bar{q}$ final states at NNLO QCD with the nested soft-collinear subtraction scheme,” arXiv:1907.05398 [hep-ph].
  * [28] M. Czakon, “A novel subtraction scheme for double-real radiation at NNLO,” Phys. Lett. B693 (2010) 259–268, arXiv:1005.0274 [hep-ph].
  * [29] M. Czakon, “Double-real radiation in hadronic top quark pair production as a proof of a certain concept,” Nucl. Phys. B849 (2011) 250–295, arXiv:1101.0642 [hep-ph].
  * [30] M. Czakon and D. Heymes, “Four-dimensional formulation of the sector-improved residue subtraction scheme,” Nucl. Phys. B890 (2014) 152–227, arXiv:1408.2500 [hep-ph].
  * [31] M. Czakon, A. van Hameren, A. Mitov, and R. Poncelet, “Single-jet inclusive rates with exact color at $\mathcal{O}(\alpha_{s}^{4})$,” arXiv:1907.12911 [hep-ph].
  * [32] G. ’t Hooft and M. J. G. Veltman, “Regularization and Renormalization of Gauge Fields,” Nucl. Phys. B44 (1972) 189–213.
  * [33] J. F. Ashmore, “A Method of Gauge Invariant Regularization,” Lett. Nuovo Cim. 4 (1972) 289–290.
  * [34] G. M. Cicuta and E. Montaldi, “Analytic renormalization via continuous space dimension,” Lett. Nuovo Cim. 4 (1972) 329–332.
  * [35] C. G. Bollini and J. J. Giambiagi, “Dimensional Renormalization: The Number of Dimensions as a Regularizing Parameter,” Nuovo Cim. B12 (1972) 20–26.
  * [36] W. J. Marciano and A. Sirlin, “Dimensional Regularization of Infrared Divergences,” Nucl. Phys. B88 (1975) 86–98.
  * [37] T. Kinoshita, “Mass singularities of Feynman amplitudes,” J. Math. Phys. 3 (1962) 650–677.
  * [38] T. D. Lee and M. Nauenberg, “Degenerate Systems and Mass Singularities,” Phys. Rev. 133 (1964) B1549–B1562.
  * [39] T. Binoth and G. Heinrich, “An automatized algorithm to compute infrared divergent multiloop integrals,” Nucl. Phys. B585 (2000) 741–759, arXiv:hep-ph/0004013 [hep-ph].
  * [40] C. Anastasiou, K. Melnikov, and F. Petriello, “A new method for real radiation at NNLO,” Phys. Rev. D69 (2004) 076010, arXiv:hep-ph/0311311 [hep-ph].
  * [41] T. Binoth and G. Heinrich, “Numerical evaluation of phase space integrals by sector decomposition,” Nucl. Phys. B693 (2004) 134–148, arXiv:hep-ph/0402265 [hep-ph].
  * [42] S. Catani, “The Singular behavior of QCD amplitudes at two loop order,” Phys. Lett. B427 (1998) 161–171, arXiv:hep-ph/9802439 [hep-ph].
  * [43] S. Catani, S. Dittmaier, and Z. Trocsanyi, “One loop singular behavior of QCD and SUSY QCD amplitudes with massive partons,” Phys. Lett. B500 (2001) 149–160, arXiv:hep-ph/0011222 [hep-ph].
  * [44] S. M. Aybat, L. J. Dixon, and G. F. Sterman, “The Two-loop soft anomalous dimension matrix and resummation at next-to-next-to leading pole,” Phys. Rev. D74 (2006) 074004, arXiv:hep-ph/0607309 [hep-ph].
  * [45] T. Becher and M. Neubert, “Infrared singularities of QCD amplitudes with massive partons,” Phys. Rev. D79 (2009) 125004, arXiv:0904.1021 [hep-ph]. [Erratum: Phys. Rev. D80 (2009) 109901].
  * [46] M. Czakon, A. Mitov, and G. F. Sterman, “Threshold Resummation for Top-Pair Hadroproduction to Next-to-Next-to-Leading Log,” Phys. Rev. D80 (2009) 074017, arXiv:0907.1790 [hep-ph].
  * [47] A. Mitov, G. F. Sterman, and I. Sung, “The Massive Soft Anomalous Dimension Matrix at Two Loops,” Phys. Rev. D79 (2009) 094015, arXiv:0903.3241 [hep-ph].
  * [48] A. Ferroglia, M. Neubert, B. D. Pecjak, and L. L. Yang, “Two-loop divergences of massive scattering amplitudes in non-abelian gauge theories,” JHEP 11 (2009) 062, arXiv:0908.3676 [hep-ph].
  * [49] A. Mitov, G. F. Sterman, and I. Sung, “Computation of the Soft Anomalous Dimension Matrix in Coordinate Space,” Phys. Rev. D82 (2010) 034020, arXiv:1005.4646 [hep-ph].
  * [50] T. Becher and M. Neubert, “Infrared singularities of scattering amplitudes in perturbative QCD,” Phys. Rev. Lett. 102 (2009) 162001, arXiv:0901.0722 [hep-ph]. [Erratum: Phys. Rev. Lett. 111 (2013) 199905].
  * [51] S. Catani and M. H. Seymour, “A General algorithm for calculating jet cross-sections in NLO QCD,” Nucl. Phys. B485 (1997) 291–419, arXiv:hep-ph/9605323 [hep-ph]. [Erratum: Nucl. Phys. B510 (1998) 503].
  * [52] D. J. Broadhurst, N. Gray, and K. Schilcher, “Gauge invariant on-shell Z(2) in QED, QCD and the effective field theory of a static quark,” Z. Phys. C52 (1991) 111–122.
  * [53] N. Gray, D. J. Broadhurst, W. Grafe, and K. Schilcher, “Three Loop Relation of Quark (Modified) Ms and Pole Masses,” Z. Phys. C48 (1990) 673–680.
  * [54] S. Frixione, Z. Kunszt, and A. Signer, “Three jet cross-sections to next-to-leading order,” Nucl. Phys. B467 (1996) 399–442, arXiv:hep-ph/9512328 [hep-ph].
  * [55] S. Frixione, “A General approach to jet cross-sections in QCD,” Nucl. Phys. B507 (1997) 295–314, arXiv:hep-ph/9706545 [hep-ph].
  * [56] S. Catani and M. H. Seymour, “The Dipole formalism for the calculation of QCD jet cross-sections at next-to-leading order,” Phys. Lett. B378 (1996) 287–301, arXiv:hep-ph/9602277 [hep-ph].
  * [57] S. Catani, S. Dittmaier, M. H. Seymour, and Z. Trocsanyi, “The Dipole formalism for next-to-leading order QCD calculations with massive partons,” Nucl. Phys. B627 (2002) 189–265, arXiv:hep-ph/0201036 [hep-ph].
  * [58] Particle Data Group Collaboration, M. Tanabashi et al., “Review of Particle Physics,” Phys. Rev. D98 no. 3, (2018) 030001.
  * [59] L. J. Dixon, “Calculating scattering amplitudes efficiently,” in QCD and beyond. Proceedings, Theoretical Advanced Study Institute in Elementary Particle Physics, TASI-95, Boulder, USA, June 4-30, 1995, pp. 539–584.  1996\.  arXiv:hep-ph/9601359 [hep-ph].
  * [60] M. Brucherseifer, F. Caola, and K. Melnikov, “$\mathcal{O}(\alpha_{s}^{2})$ corrections to fully-differential top quark decays,” JHEP 04 (2013) 059, arXiv:1301.7133 [hep-ph].
  * [61] G. Passarino and M. J. G. Veltman, “One Loop Corrections for $e^{+}e^{-}$ Annihilation Into $\mu^{+}\mu^{-}$ in the Weinberg Model,” Nucl. Phys. B160 (1979) 151–207.
  * [62] R. K. Ellis and G. Zanderighi, “Scalar one-loop integrals for QCD,” JHEP 02 (2008) 002, arXiv:0712.1851 [hep-ph].
  * [63] S. Carrazza, R. K. Ellis, and G. Zanderighi, “QCDLoop: a comprehensive framework for one-loop scalar integrals,” Comput. Phys. Commun. 209 (2016) 134–143, arXiv:1605.03181 [hep-ph].
  * [64] J. Ablinger, A. Behring, J. Blümlein, G. Falcioni, A. De Freitas, P. Marquard, N. Rana, and C. Schneider, “Heavy quark form factors at two loops,” Phys. Rev. D97 no. 9, (2018) 094022, arXiv:1712.09889 [hep-ph].
  * [65] W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, P. Mastrolia, and E. Remiddi, “Decays of scalar and pseudoscalar Higgs bosons into fermions: Two-loop QCD corrections to the Higgs-quark-antiquark amplitude,” Phys. Rev. D72 (2005) 096002, arXiv:hep-ph/0508254 [hep-ph].
  * [66] E. Remiddi and J. A. M. Vermaseren, “Harmonic polylogarithms,” Int. J. Mod. Phys. A15 (2000) 725–754, arXiv:hep-ph/9905237 [hep-ph].
  * [67] T. Gehrmann and E. Remiddi, “Numerical evaluation of harmonic polylogarithms,” Comput. Phys. Commun. 141 (2001) 296–312, arXiv:hep-ph/0107173 [hep-ph].
  * [68] J. Gluza, A. Mitov, S. Moch, and T. Riemann, “The QCD form factor of heavy quarks at NNLO,” JHEP 07 (2009) 001, arXiv:0905.1137 [hep-ph].
  * [69] K. G. Chetyrkin, J. H. Kühn, and M. Steinhauser, “RunDec: A Mathematica package for running and decoupling of the strong coupling and quark masses,” Comput. Phys. Commun. 133 (2000) 43–65, arXiv:hep-ph/0004189 [hep-ph].
  * [70] F. Herren and M. Steinhauser, “Version 3 of RunDec and CRunDec,” Comput. Phys. Commun. 224 (2018) 333–345, arXiv:1703.03751 [hep-ph].
  * [71] K. G. Chetyrkin, “Correlator of the quark scalar currents and $\Gamma_{\text{tot}}(H\to\text{hadrons})$ at $O(\alpha_{s}^{3})$ in pQCD,” Phys. Lett. B390 (1997) 309–317, arXiv:hep-ph/9608318 [hep-ph].
  * [72] S. Catani, Y. L. Dokshitzer, M. Olsson, G. Turnock, and B. R. Webber, “New clustering algorithm for multi-jet cross-sections in $e^{+}e^{-}$ annihilation,” Phys. Lett. B269 (1991) 432–438.
  * [73] M. Cacciari, G. P. Salam, and G. Soyez, “FastJet User Manual,” Eur. Phys. J. C72 (2012) 1896, arXiv:1111.6097 [hep-ph].
  * [74] D. J. Gross and F. Wilczek, “Ultraviolet Behavior of Nonabelian Gauge Theories,” Phys. Rev. Lett. 30 (1973) 1343–1346.
  * [75] H. D. Politzer, “Reliable Perturbative Results for Strong Interactions?,” Phys. Rev. Lett. 30 (1973) 1346–1349.
  * [76] K. Melnikov and T. van Ritbergen, “The Three loop on-shell renormalization of QCD and QED,” Nucl. Phys. B591 (2000) 515–546, arXiv:hep-ph/0005131 [hep-ph].
  * [77] W. Beenakker, S. Dittmaier, M. Kramer, B. Plumper, M. Spira, and P. M. Zerwas, “NLO QCD corrections to $t\bar{t}H$ production in hadron collisions,” Nucl. Phys. B653 (2003) 151–203, arXiv:hep-ph/0211352 [hep-ph].
  * [78] M. Czakon, A. Mitov, and S. Moch, “Heavy-quark production in gluon fusion at two loops in QCD,” Nucl. Phys. B798 (2008) 210–250, arXiv:0707.4139 [hep-ph].
  * [79] S. Weinberg, “Effective Gauge Theories,” Phys. Lett. 91B (1980) 51–55.
  * [80] B. A. Ovrut and H. J. Schnitzer, “The Decoupling Theorem and Minimal Subtraction,” Phys. Lett. 100B (1981) 403–406.
  * [81] W. Wetzel, “Minimal Subtraction and the Decoupling of Heavy Quarks for Arbitrary Values of the Gauge Parameter,” Nucl. Phys. B196 (1982) 259–272.
  * [82] W. Bernreuther and W. Wetzel, “Decoupling of Heavy Quarks in the Minimal Subtraction Scheme,” Nucl. Phys. B197 (1982) 228–236. [Erratum: Nucl. Phys. B513 (1998) 758].
  * [83] W. Bernreuther, “Decoupling of Heavy Quarks in Quantum Chromodynamics,” Annals Phys. 151 (1983) 127.
  * [84] W. Bernreuther, “Heavy Quark Effects on the Parameters of Quantum Chromodynamics Defined by Minimal Subtraction,” Z. Phys. C20 (1983) 331–333.
  * [85] K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, “Decoupling relations to $O(\alpha_{s}^{3})$ and their connection to low-energy theorems,” Nucl. Phys. B510 (1998) 61–87, arXiv:hep-ph/9708255 [hep-ph].
  * [86] M. Gerlach, G. Mishima, and M. Steinhauser, “Matching coefficients in nonrelativistic QCD to two-loop accuracy,” Phys. Rev. D100 no. 5, (2019) 054016, arXiv:1907.08227 [hep-ph].
  * [87] K. Melnikov and T. v. Ritbergen, “The Three loop relation between the $\overline{\mathrm{MS}}$ and the pole quark masses,” Phys. Lett. B482 (2000) 99–108, arXiv:hep-ph/9912391 [hep-ph].
  * [88] A. Mitov and S. Moch, “The Singular behavior of massive QCD amplitudes,” JHEP 05 (2007) 001, arXiv:hep-ph/0612149 [hep-ph].
  * [89] I. Bierenbaum, M. Czakon, and A. Mitov, “The singular behavior of one-loop massive QCD amplitudes with one external soft gluon,” Nucl. Phys. B856 (2012) 228–246, arXiv:1107.4384 [hep-ph].
  * [90] M. Czakon and A. Mitov, “A simplified expression for the one-loop soft-gluon current with massive fermions,” arXiv:1804.02069 [hep-ph].
  * [91] S. Catani and M. Grazzini, “Infrared factorization of tree level QCD amplitudes at the next-to-next-to-leading order and beyond,” Nucl. Phys. B570 (2000) 287–325, arXiv:hep-ph/9908523 [hep-ph].
  * [92] S. Alioli, P. Nason, C. Oleari, and E. Re, “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX,” JHEP 06 (2010) 043, arXiv:1002.2581 [hep-ph].
  * [93] J. Ablinger, “A Computer Algebra Toolbox for Harmonic Sums Related to Particle Physics,” 2009.
  * [94] J. Ablinger, Computer Algebra Algorithms for Special Functions in Particle Physics.  PhD thesis, J. Kepler University Linz, 2012.  arXiv:1305.0687 [math-ph].
  * [95] J. A. M. Vermaseren, “Harmonic sums, Mellin transforms and integrals,” Int. J. Mod. Phys. A14 (1999) 2037–2076, arXiv:hep-ph/9806280 [hep-ph].
  * [96] J. Blümlein, “Structural Relations of Harmonic Sums and Mellin Transforms up to Weight $w=5$,” Comput. Phys. Commun. 180 (2009) 2218–2249, arXiv:0901.3106 [hep-ph].
  * [97] J. Ablinger, J. Blümlein, and C. Schneider, “Harmonic Sums and Polylogarithms Generated by Cyclotomic Polynomials,” J. Math. Phys. 52 (2011) 102301, arXiv:1105.6063 [math-ph].
  * [98] J. Ablinger, J. Blümlein, and C. Schneider, “Analytic and Algorithmic Aspects of Generalized Harmonic Sums and Polylogarithms,” J. Math. Phys. 54 (2013) 082301, arXiv:1302.0378 [math-ph].
  * [99] J. Ablinger, “The package HarmonicSums: Computer Algebra and Analytic aspects of Nested Sums,” PoS LL2014 (2014) 019, arXiv:1407.6180 [cs.SC].
  * [100] J. Ablinger, “Inverse Mellin Transform of Holonomic Sequences,” PoS LL2016 (2016) 067, arXiv:1606.02845 [cs.SC].
  * [101] J. Ablinger, “Computing the Inverse Mellin Transform of Holonomic Sequences using Kovacic’s Algorithm,” PoS RADCOR2017 (2017) 069, arXiv:1801.01039 [cs.SC].
  * [102] J. Ablinger, “Discovering and Proving Infinite Pochhammer Sum Identities,” arXiv:1902.11001 [math.CO].


# From generalized arithmetic means to geodesics to Hamilton dynamics to
Bregman divergences

Henryk Gzyl,
Centro de Finanzas IESA, Caracas, Venezuela.
<EMAIL_ADDRESS>
###### Abstract

Here we examine some connections between the notions of generalized arithmetic
means, geodesics, Lagrange-Hamilton dynamics and Bregman divergences. In a
previous paper we developed a predictive interpretation of generalized
arithmetic means. That work was more probabilistically oriented. Here we take
a geometric turn, and see that generalized arithmetic means actually minimize
a geodesic distance on ${\mathbb{R}}^{n}.$ Such metrics might result from
pull-backs of the Euclidean metric in ${\mathbb{R}}^{n}.$ We shall furthermore
see that in some cases these pull-backs might coincide with the Hessian of a
convex function. This occurs when the Hessian of a convex function has a
square root that is the Jacobian of a diffeomorphism in ${\mathbb{R}}^{n}.$ In
this case we obtain a comparison between the Bregman divergence defined by the
convex function and the geodesic distance in the metric defined by its
Hessian.

Keywords: Generalized arithmetic means, Riemannian distances, Hamilton-
Lagrange equivalence, Bregman divergence.
MSC 2000 Subject Classification 53A3,5, 53A40, 52B55, 53C21, 70H06, 70H15,
37J05.

###### Contents

  1. 1 Introduction and Preliminaries
    1. 1.1 Organization of the paper
  2. 2 Geodesics and one-dimensional changes of scale
  3. 3 Geodesics induced by a diffeomorphism
  4. 4 A detour into classical mechanics
  5. 5 Hessians of convex functions and diffeomorphisms in ${\mathbb{R}}^{n}$
    1. 5.1 From a convex function to diffeomorphisms
    2. 5.2 From diffeomorphisms to convex functions
  6. 6 A comparison result: Divergence versus geodesic distance
  7. 7 Examples
    1. 7.1 Extended generalized arithmetic means
    2. 7.2 Means defined by a flow
    3. 7.3 Harmonic means in ${\mathbb{R}}^{n}$
    4. 7.4 Multidimensional best predictors
    5. 7.5 Generalized arithmetic means and convex functions defined by the gradient of a convex function
    6. 7.6 Geodesics and distances determined by the Fenchel-Legendre conjugate

##  1 Introduction and Preliminaries

In [5] a reinterpretation of arithmetic means as best predictors was
presented. To explain, consider a strictly increasing function
$u:{\mathcal{I}}\to{\mathbb{R}},$ where ${\mathcal{I}}\subset{\mathbb{R}}$ is
an open interval. If $\\{x_{1},...,x_{n}\\}$ is a collection of points in
${\mathcal{I}},$ their generalized arithmetic mean

$c=u^{-1}\Big{(}\frac{1}{n}\sum_{i=1}^{n}u(x_{i})\Big{)}$ (1.1)

is the point at which the distance
$\Big{(}\sum(u(x_{i})-u(c))^{2}\Big{)}^{1/2}$ is minimal. This approach leads
to a notion of best predictor for random variables with extends the notion of
expected value and conditional expected value. At this point it is appropriate
to mention that this variational argument to obtain (1.1) can be traced back
at least to [2]. From a historical point of view, the generalized arithmetic
mean abstracts the notion of certain equivalent of a random cash flow for an
investor with utility function $u.$

###  1.1 Organization of the paper

The aim of this work is to look into the geometric ideas behind that
construction. For that we shall use $u(x)$ to define an appropriate metric on
${\mathcal{M}}$ and verify that the distance

$d^{2}_{u}({\boldsymbol{x}},{\boldsymbol{y}})=\Big{(}\sum(u(x_{i})-u(y_{i}))^{2}\Big{)}^{1/2}$
(1.2)

is actually the geodesic distance between ${\boldsymbol{x}}$ and
${\boldsymbol{y}}.$ This simple case is nice because it leads to a diagonal
metric in ${\mathcal{M}}$ and the realization of the metric as the Hessian of
a convex function is trivial. This is carried out in Section 2. In Section 3
we consider a more general diffeomorphism ${\boldsymbol{U}}$ defined on
${\mathcal{M}},$ which leads to a metric such that the geodesic distance in
that metric determines a multidimensional version of the generalized
arithmetic mean.

We devote Section 4 to consider the changes of variables and the integration
of the geodesic equations within the framework of Hamilton-Lagrange formalism.
The idea is to regard the integration of the equations of the geodesics from
another point of view. It is in Section 5 where we establish conditions for
the two way connection between Hessian matrices of convex functions and
diffeomorphisms in ${\mathbb{R}}^{n},$ while in Section 6 we compare the
geodesic distance to the Bregman divergence defined by the convex function. In
the last section we collect some examples and pending issues.
A word about notation, etc.
In the remainder of this paper we use ${\mathcal{M}}=({\mathcal{I}})^{n},$
where as said, ${\mathcal{I}}$ is a an interval (bounded or unbounded). We use
this setup because in many cases we want to consider convex functions of the
type $\sum u(x_{i})$ with $u(x)$ convex. As the manifold that we are
considering is a simple subset of ${\mathbb{R}}^{n}$ naturally provided with a
global chart, the standard constructs of differential geometry in this case
are very simple. For example, the tangent bundle and the cotangent bundles are
trivially identifiable with ${\mathcal{I}}\times{\mathbb{R}}^{n}.$ We use the
standard $\langle{\boldsymbol{v}}_{1},{\boldsymbol{v}}_{2}\rangle$ to denote
the usual Euclidean scalar product between the vectors
${\boldsymbol{v}}_{1},{\boldsymbol{v}}_{2}.$ And to finish, we use the
standard $h^{\prime}(x)$ to define the derivative of the function
$h(x):{\mathcal{I}}\to{\mathbb{R}},$ and $\dot{x}$ and $\ddot{x}$ to define
the first two derivatives of $x(t):(a,b)\to{\mathbb{R}}$ where $t$ is thought
of as “time.”

##  2 Geodesics and one-dimensional changes of scale

Let $u:{\mathcal{I}}\to{\mathbb{R}}.$ be a twice continuously differentiable,
strictly increasing function and $u^{\prime}(x)$ its derivative. At each
${\boldsymbol{x}}\in{\mathcal{I}}$ define the Riemannian metric (on the
tangent space to ${\mathcal{M}}$ at ${\boldsymbol{x}}$) by

$g_{i,j}({\boldsymbol{x}})=g(x_{i})\delta_{i,j}=u^{\prime}(x_{i})^{2}\delta_{i,j}.$
(2.1)

That is, not only is the coordinate system orthogonal (that is, the metric is
diagonal), but it is separable as well, that is, $g_{i,i}$ is a function of
$x_{i}$ only. The equation of the geodesic that minimizes the distance between
two points ${\boldsymbol{x}},{\boldsymbol{y}}\in{\mathcal{M}}$ is obtained
minimizing

$\int_{0}^{1}\Big{(}\sum_{i,j}g_{i,j}({\boldsymbol{x}})\dot{x}_{i}\dot{x_{j}})\Big{)}^{1/2}dt$
(2.2)

over the class of continuous functions
${\boldsymbol{x}}:[0,1]\to{\mathcal{M}},$ twice continuously differentiable on
$(0,1).$ It is a standard result, that in this case, the equations of the
geodesics are (see [7], for example):

$\frac{d}{dt}\Big{(}\frac{\partial}{\partial\dot{x}_{k}}\frac{1}{2}\sum_{i,j}g_{i,j}({\boldsymbol{x}})\dot{x}_{i}\dot{x_{j}}\Big{)}=\frac{\partial}{\partial
x_{k}}\Big{(}\frac{1}{2}\sum_{i,j}g_{i,j}({\boldsymbol{x}})\dot{x}_{i}\dot{x_{j}}\Big{)},\;\;\;k=1,...,n.$
(2.3)

In this special diagonal metric, if we put $g^{\prime}(x)=dg(x)/dx,$ the
equations of the geodesic are

$g(x)\ddot{x}_{k}+\frac{1}{2}g^{\prime}(x_{i})(\dot{x}_{k})^{2},\;\;k=1,...,n.$

Substituting $g(x_{k})=(u^{\prime}(x_{k})$ we obtain

$u^{\prime}(x_{k})\ddot{x}_{k}+u^{\prime\prime}(x_{k})\dot{x}=0,\;\;\;k=1,...,n$

This can be simply integrated. At the first step, note that it we obtain that
$u^{\prime}(x_{k}(t))\dot{x}_{k}=C_{k},$ where $C_{k}$ is some unknown
constant. This equation can be trivially integrated to obtain
$u(x_{k}(t))=u(x_{k}(0))+tC_{k}.$ To determine $C_{k}$ make use of the fact
that $x_{k}(0)=x_{k}$ and $x_{k}(1)=y_{k}$ to obtain
$C_{k}=u(y_{k})-u(x_{k}).$ To close the circle, note that the geodesic
distance between ${\boldsymbol{x}}$ and ${\boldsymbol{y}}$ is, according to
(2.2)

$\int_{0}^{1}\Big{(}\sum_{i,j}g_{i,j}({\boldsymbol{x}})\dot{x}_{i}\dot{x_{j}})\Big{)}^{1/2}dt=\int_{0}^{1}\Big{(}\sum_{i}C_{i}^{2}\Big{)}^{1/2}dt=\Big{(}\sum_{i}(u(y_{k})-u(x_{k}))^{2}\Big{)}^{1/2}$

which equals $d_{u}({\boldsymbol{x}},{\boldsymbol{y}}).$

##  3 Geodesics induced by a diffeomorphism

Now we extend the setup of the previous section to the general case. Let
${\boldsymbol{U}}:{\mathcal{M}}\to{\mathbb{R}}^{n}$ be a twice continuously
differentiable diffeomorphism between ${\mathcal{M}}$ and
${\boldsymbol{U}}({\mathcal{M}}).$ On the tangent bundle
${\boldsymbol{U}}({\mathcal{M}})\times{\mathbb{R}}^{n}$ we consider the
Euclidean metric. Its pullback to ${\mathcal{M}}\times{\mathbb{R}}^{n}$ at
${\boldsymbol{x}}\in{\mathcal{M}}$ is given by

$g_{i,j}({\boldsymbol{x}})=\sum_{k=1}^{n}U^{k}_{i}U^{k}_{j},\;\;for\;\;1\leq
i,j\leq n.$ (3.1)

Here, $U^{k}({\boldsymbol{x}})$ denotes the $k-$th component of
${\boldsymbol{U}}$ and $U^{k}_{i}$ denotes $\partial U^{k}/\partial x_{j}.$ We
denote by $V^{k}_{j}$ the inverse of $U^{k}_{j},$ that is, the Jacobian of the
inverse of ${\boldsymbol{U}}.$ Again, the geodesics minimize the distance
given by (2.2) with $g_{i,j}$ given by (3.1). If $t\to{\boldsymbol{x}}(t)$ is
a geodesic from ${\boldsymbol{x}}(0)={\boldsymbol{x}}$ to
${\boldsymbol{x}}(1)={\boldsymbol{y}},$ then it must satisfy equation (2.3).

If in (2.3) we substitute (3.1), after some simple but boring arithmetics we
obtain

$\sum_{i=1}^{n}U^{i}_{k}\Big{(}\sum_{n=1}^{n}U^{i}_{n}\ddot{x}_{n}+\sum_{k,l}U^{i}_{n,l}\dot{x}_{n}\dot{x}_{l}\Big{)}=0,\;\;\;k=1,...,n.$

Now, multiply both sides by $V^{k}_{j}$ and sum over $k=1,...,n$ to obtain

$\sum_{n=1}^{n}U^{j}_{n}\ddot{x}_{n}+\sum_{k,l}U^{j}_{n,l}\dot{x}_{n}\dot{x}_{l}=0,\;\;\;j=1,...,n.$

This can be written as

$\frac{d}{dt}\Big{(}\sum_{n=1}^{n}U^{j}_{n}\dot{x}_{n}\Big{)}=0,\;\;\;j=1,...,n.$

Which implies that

$\frac{d}{dt}U^{j}({\boldsymbol{x}}(t)=\sum_{n=1}^{n}U^{j}({\boldsymbol{x}}(t))\dot{x}_{j}(t)=C_{j},\;\;\;j=1,...,n$

for some constants $C_{j},j=1,...,n.$ Again, This implies that
$U^{j}({\boldsymbol{x}}(t))=U^{j}({\boldsymbol{x}}(0))+tC_{j}$ and that
$C_{j}=U^{j}({\boldsymbol{y}})-U^{j}({\boldsymbol{x}}),\,j=1,...,n.$ As above,
inserting this into the definition of geodesic distance between
${\boldsymbol{x}}$ and ${\boldsymbol{y}}$ we obtain

$\int_{0}^{1}\Big{(}\sum_{i,j}g_{i,j}({\boldsymbol{x}})\dot{x}_{i}\dot{x_{j}})\Big{)}^{1/2}dt=\int_{0}^{1}\Big{(}\sum_{i}C_{i}^{2}\Big{)}^{1/2}dt=\Big{(}\sum_{i}(U^{k}({\boldsymbol{y}})-U^{k}({\boldsymbol{x}}))^{2}\Big{)}^{1/2}=d_{U}({\boldsymbol{x}},{\boldsymbol{y}}).$
(3.2)

Notice that $d_{U}({\boldsymbol{x}},{\boldsymbol{y}})$ is positive, symmetric,
satisfies the triangle inequality, and since ${\boldsymbol{U}}$ is a
diffeomorphism,
$d_{U}({\boldsymbol{x}},{\boldsymbol{y}})=0\Leftrightarrow{\boldsymbol{x}}={\boldsymbol{y}},$
thus the notation is consistent. The comparison with the first case is clear:
There $U^{k}({\boldsymbol{x}})=u(x_{k})$ for $k=1,...,n.$

##  4 A detour into classical mechanics

Here we address the issue of integrating the geodesic equations from the point
of view of classical mechanics. This is to explain the trivial integration of
the geodesic equations within the context of an elegant framework. We refer
the reader to [1] for the essential notions. Cutting some edges, the general
approach to Newton’s equations, consists of specifying some manifold
(${\mathcal{M}}$ in our case), on whose tangent bundle
${\mathcal{M}}\times{\mathbb{R}}^{n}$ a Lagrangian function
$L({\boldsymbol{x}},\dot{{\boldsymbol{x}}})$ is defined. A trajectory of the
system between points ${\boldsymbol{x}}(1)$ and ${\boldsymbol{x}}(2)$ at times
$t_{1},t_{2}$ respectively, is the curve that minimizes the (action) integral:

$\int_{t_{1}}^{t_{2}}L({\boldsymbol{x}}(t),\dot{{\boldsymbol{x}}}(t),t)dt$

over the class of all (twice continuously) differentiable curves joining the
said points at the given times. It can be shown, see [1], that such trajectory
satisfies (the Euler-Lagrange) extension of the Newton’s equation of motion:

$\frac{d}{dt}\Big{(}\frac{\partial}{\partial\dot{x}_{k}}L({\boldsymbol{x}}(t),\dot{{\boldsymbol{x}}}(t),t)\Big{)}=\frac{\partial}{\partial
x_{k}}L({\boldsymbol{x}}(t),\dot{{\boldsymbol{x}}}(t),t).$ (4.1)

If we put

$L({\boldsymbol{x}}(t),\dot{{\boldsymbol{x}}}(t))=\frac{1}{2}\sum_{i,j}g_{i,j}({\boldsymbol{x}})\dot{x}_{i}\dot{x}_{j},$

we recognize the equations of the geodesics as the Euler-Lagrange equations of
some dynamical system. Notice that if we consider the curves
${\boldsymbol{y}}(t)={\boldsymbol{U}}({\boldsymbol{x}}(t))$ in
${\boldsymbol{U}}({\mathcal{M}}),$ then the Lagrangian in these coordinates
becomes

$L({\boldsymbol{y}},\dot{{\boldsymbol{y}}},t)=\frac{1}{2}\sum_{k=1}^{n}(\dot{y}_{k})^{2}$

and the Euler-Lagrange (newton) equations become $\ddot{y}_{k}=0,$ namely the
equations of straight lines. We already saw that in Section 3, hidden in the
notation. Another way of understanding that change of variables is to begin
with the Hamiltonian associated to the Lagrangian. This is defined by
introducing the momentum variables $p_{k}$ by

$p_{k}=\frac{\partial}{\partial\dot{x}_{k}}L({\boldsymbol{x}}(t),\dot{{\boldsymbol{x}}}(t),t)=\sum_{j=1}^{n}g_{j,k}({\boldsymbol{x}})\dot{x}_{j}.$

Under the assumption of invertibility of the Jacobian $U^{i}_{j},$ solving for
$\dot{{\boldsymbol{x}}}$ and substituting in the expression for $L$ we obtain:

$H({\boldsymbol{x}},{\boldsymbol{p}})=\frac{1}{2}\sum_{k,j}g^{k,j}({\boldsymbol{x}})p_{i}p_{j}.$

where we use the conventional $g^{k,j}$ to denote the inverse of $g_{k,j}.$
The analogue of the Euler-Lagrange equations are the Hamilton equations for
$({\boldsymbol{x}},{\boldsymbol{p}})$ is:

$\dot{x}_{k}=\frac{\partial H}{\partial
p_{k}}\;\;\;\dot{p}_{k}=-\frac{\partial H}{\partial x_{k}}.$ (4.2)

To complete explaining the change of variables we need to show that there
exists a mapping
$({\boldsymbol{x}},{\boldsymbol{p}})\to({\boldsymbol{y}},{\boldsymbol{\pi}})$
from ${\mathcal{M}}\times{\mathbb{R}}^{n}$ to
${\boldsymbol{U}}({\mathcal{M}})\times{\mathbb{R}}^{n},$ under which the
Hamiltonian becomes

$\tilde{H}=\frac{1}{2}\sum_{k=1}^{n}(\pi_{k})^{2},$

and that the change of variables preserves the form of the Hamiltonian
equations of motion (4.2). The mapping is:

$y^{i}=U^{i}({\boldsymbol{x}})\;\;\;\;\mbox{and}\;\;\;\pi_{i}=\sum_{k=1}p_{j}V^{j}_{i};\;\;\;i=1,...,n$

The condition for this mapping to preserve the form of the equations is that
the following holds for the given Poisson “brackets”

$[y^{i},y^{j}]=[\pi_{i},\pi_{j}]=0,\;\;\;\mbox{and}\;\;\;[y^{i},\pi_{j}]=\delta_{i,j},\;\;i,j=1,...n$

which, for any pair of continuously differentiable functions $f,g$ on
${\mathcal{M}}\times{\mathbb{R}}^{n}$ is given by:

$[f,g]=\sum_{k=1}^{n}\frac{\partial f}{\partial x_{k}}\frac{\partial
g}{\partial p_{k}}-\frac{\partial f}{\partial p_{k}}\frac{\partial g}{\partial
x_{k}}.$

As we said, we were to cut a lot of corners. See [1] for full details. That
the system is trivially integrable, means that (4.2), in the new coordinates
becomes

$\dot{y}^{i}=\frac{\partial\tilde{H}}{\partial\pi_{i}}=\pi_{i}\;\;\mbox{and}\;\;\dot{\pi}_{i}=-\frac{\partial\tilde{H}}{\partial
y^{i}}=0,$

which imply that $\pi_{i}=\rm{constant}=C_{i}$ and $y^{i}=y^{i}(0)+tC_{i}.$
That is the constants of integration that determine the solutions to the
equations of the geodesics are the constant values of the momenta. For
trajectories go from a point ${\boldsymbol{x}}$ to a point ${\boldsymbol{y}}$
in a unit of time, the momenta are related to the distance between the two
points by $y^{i}(1)-y^{i}(0)=C^{i}.$

##  5 Hessians of convex functions and diffeomorphisms in ${\mathbb{R}}^{n}$

In the two (sub)sections that come up next, we do two things. First we extend
the results in [5] that were summed up in Section 2, to a more general convex
function, after that we examine in which case a diffeomorphism
${\boldsymbol{U}}$ determines a convex function whose Hessian coincides with
the Riemann metric given by (3.1).

###  5.1 From a convex function to diffeomorphisms

Let us begin by stating some necessary properties of the convex function.

###### Assumptions 5.1

Let us now consider an open, convex subset ${\mathcal{M}}$ of
${\mathbb{R}}^{n}$ and a strictly convex function
$\Phi:{\mathcal{M}}\to{\mathbb{R}}^{n}$ satisfying
1 $\Phi$ is at least twice continuously differentiable in all variables.
Denote its Hessian matrix by $\Phi^{\prime\prime},$ that is,
$(\Phi^{\prime\prime})_{i,j}=\partial\Phi/\partial x_{i}\partial x_{j}.$ The
strict convexity of $\Phi$ implies that the Hessian $\Phi^{\prime\prime}$ is
positive definite. Let us denote by $k({\boldsymbol{x}})$ its smallest
eigenvalue and let us suppose that:
2 There exists $a>0$ such that $a\leq k(x)$ for all
${\boldsymbol{x}}\in{\mathcal{M}},$ or equivalently

$0<a\langle{\boldsymbol{\xi}},{\boldsymbol{\xi}}\rangle\leq
k({\boldsymbol{x}})\langle{\boldsymbol{\xi}},{\boldsymbol{\xi}}\rangle\leq\langle{\boldsymbol{\xi}},\Phi^{\prime\prime}{\boldsymbol{\xi}}\rangle,\;\;\;\forall\;{\boldsymbol{\xi}}\in{\mathbb{R}}^{n}.$
(5.1)

And we now have:

###### Theorem 5.1

With the notations introduced above, let us suppose that there exists a
continuously differentiable square root $S$ of $\Phi^{\prime\prime},$ that is,
$\Phi^{\prime\prime}=S^{t}S$ such that:
a $S({\boldsymbol{x}})$ is invertible at each
${\boldsymbol{x}}\in{\mathcal{M}}.$
b For every $i=1,...,n$ we have $\partial S_{i,j}/\partial x_{k}=\partial
S_{i,k}/\partial x_{j},$ for all $1\leq j,k\leq n.$
Fix some ${\boldsymbol{x}}_{0}\in{\mathcal{M}}$ and define

$U^{i}({\boldsymbol{x}})=\int S_{i,j}(\gamma(s))\dot{\gamma}(s)ds$

where $\gamma$ denotes any continuously differentiable trajectory between
${\boldsymbol{x}}_{0}$ and ${\boldsymbol{x}}.$ Then, the mapping
${\boldsymbol{U}}:{\mathcal{M}}\to{\mathbb{R}}^{n}$ with components
$U^{i}({\boldsymbol{x}})$ is well defined (up to a constant), and it is a
global diffeomorphism satisfying

$\Phi_{i,j}^{{}^{\prime\prime}}=\sum_{m=1}^{n}\frac{\partial U^{m}}{\partial
x_{i}}\frac{\partial U^{m}}{\partial x_{j}}$

Before proving the theorem we need the following result.

###### Lemma 5.1

With the notations introduced above, let $T$ denote the inverse $S^{-1}$ of
$S.$ Then $\|T\|\leq 1/a.$

Proof In item (2) of Assumption 5.1 replace ${\boldsymbol{\xi}}$ by
${\boldsymbol{V}}{\boldsymbol{\xi}}$ to obtain

$a\langle{\boldsymbol{V}}{\boldsymbol{\xi}},{\boldsymbol{V}}{\boldsymbol{\xi}}\rangle\leq
k({\boldsymbol{x}})\langle{\boldsymbol{V}}{\boldsymbol{\xi}},{\boldsymbol{V}}{\boldsymbol{\xi}}\rangle\leq\langle{\boldsymbol{V}}{\boldsymbol{\xi}},\Phi^{\prime\prime}{\boldsymbol{V}}{\boldsymbol{\xi}}\rangle=\langle{\boldsymbol{\xi}},{\boldsymbol{V}}^{t}\Phi^{\prime\prime}{\boldsymbol{V}}{\boldsymbol{\xi}}\rangle=\langle{\boldsymbol{\xi}},{\boldsymbol{V}}^{t}U^{t}U{\boldsymbol{V}}{\boldsymbol{\xi}}\rangle=\langle{\boldsymbol{\xi}},{\boldsymbol{\xi}}\rangle.$

In other words

$\|{\boldsymbol{V}}^{t}{\boldsymbol{V}}\|\leq 1/k(x)\leq
1/a\Rightarrow\|{\boldsymbol{V}}\|<K=(\sqrt{a^{-1}})$

since $\|{\boldsymbol{V}}\|$ is given by the square root of the largest
eigenvalue of ${\boldsymbol{V}}^{t}{\boldsymbol{V}}.$ $\Box$

Let us now complete the proof of Theorem 5.1.
Proof
Since ${\mathcal{M}}$ is convex, it is simply connected and assumption (b) in
the statement implies (via Stokes theorem) that ${\boldsymbol{U}}$ is a well
defined mapping on ${\mathcal{M}},$ and that its Jacobian $S$ is non vanishing
(by assumption), and thus ${\boldsymbol{U}}$ is a local diffeomorphism.
Invoking the previous lemma and Hadamard’s theorem, see, for example Theorem 2
in [8], we conclude that ${\boldsymbol{U}}$ is a diffeomorphism satisfying

$\Phi_{i,j}^{{}^{\prime\prime}}=\sum_{m=1}^{n}\frac{\partial U^{m}}{\partial
x_{i}}\frac{\partial U^{m}}{\partial x_{j}}$

which concludes the proof. $\Box$

###  5.2 From diffeomorphisms to convex functions

A convex function has a positive definite Hessian. The next result imposes a
condition on the diffeomorphism ${\boldsymbol{U}}$ so that the metric that it
defines, as explained in Section 3, can be the Hessian of a convex function.

###### Theorem 5.2

Let ${\mathcal{M}}$ be an open convex subset of ${\mathbb{R}}^{n}$ and
${\boldsymbol{U}}:{\mathcal{M}}\to{\mathbb{R}}^{n}.$ be a diffeomorphism
satisfying

$\sum_{m=1}^{n}\frac{\partial^{2}U^{m}}{\partial x_{k}\partial
x_{i}}\frac{\partial U^{m}}{\partial
x_{j}}=\sum_{m=1}^{n}\frac{\partial^{2}U^{m}}{\partial x_{i}\partial
x_{j}}\frac{\partial U^{m}}{\partial x_{k}}$ (5.2)

Then there exists a strictly convex function
$\Phi:{\mathcal{M}}\to{\mathbb{R}}^{n}$ such that

$\frac{\partial^{2}\Phi}{\partial x_{i}\partial
x_{j}}=\sum_{m=1}^{n}\frac{\partial U^{m}}{\partial x_{i}}\frac{\partial
U^{m}}{\partial x_{j}}$

Proof As before, $g_{i,j}=\sum_{m=1}^{n}\frac{\partial U^{m}}{\partial
x_{i}}\frac{\partial U^{m}}{\partial x_{j}}$ Fix an $1\leq i\leq n$ and an
${\boldsymbol{x}}(0)\in{\mathcal{M}}$ and consider the following curve joining
${\boldsymbol{x}}(0)$ to ${\boldsymbol{x}}$ defined piecewise by increasing
one coordinate at a time: At the $k-$th step, move along the $k-$th coordinate
axis from $x_{k}(0)$ to $x_{k},$ that is, along the line:

$(x_{1},x_{2},...,x_{k-1},\xi_{k},x_{k+1}(0),...,x_{n}(0))\;\;\;\mbox{with}\;\;\;x_{k}(0)\leq\xi_{k}\leq
x_{k}.$

Call this trajectory $\gamma.$ The import of condition(5.2) is to make the
integral defined below to be independent of the trajectory.

$A^{i}({\boldsymbol{x}})=A^{i}({\boldsymbol{x}}_{0})+\sum_{k=1}\int_{x_{k}(0)}^{x_{k}}g_{i,k}(\gamma)d\xi_{k}.$

where $A^{i}({\boldsymbol{x}}_{0})$ are constants of integration. Now, since
$\partial A^{i}/\partial x_{j}=g_{i,j}=g_{j,i}=\partial A^{j}/\partial x_{i},$
then the following line integral is also independent of the trajectory. Thus
integrating along the same piecewise trajectory we put:

$\Phi({\boldsymbol{x}})=\Phi({\boldsymbol{x}}(0)+\sum_{k=1}\int_{x_{k}(0)}^{x_{k}}A^{k}(\gamma)d\xi_{k}.$

We choose the special trajectory so that the verification that the Hessian of
$\Phi$ is $g$ is trivial. Thus we prove the claim. $\Box$

##  6 A comparison result: Divergence versus geodesic distance

Bregman divergences are a common measure of discrepancy. They are used to
compare how different are two objects that can be described by points in
convex subset of some many-dimensional space. The definition goes as follows:

###### Definition 6.1

Let ${\mathcal{M}}\subset{\mathbb{R}}^{n}.$ be an open convex set
$\Phi:{\mathcal{M}}\to{\mathbb{R}}$ be a strictly convex, continuously
differentiable function and put

$\delta^{2}_{\Phi}({\boldsymbol{x}},{\boldsymbol{y}})=\Phi({\boldsymbol{x}})-\Phi({\boldsymbol{y}})-\langle({\boldsymbol{x}}-{\boldsymbol{y}}),\nabla\Phi({\boldsymbol{y}})\rangle.$
(6.1)

Comment: Since the right hand is non-negative, and vanishes if and only if
${\boldsymbol{y}}={\boldsymbol{x}},$ the notation on chosen for the left hand
side is consistent, even though in general, it is not a true distance. This is
why it is called a discrimination function between ${\boldsymbol{x}}$ and
${\boldsymbol{y}}.$

This notion was introduced by Bregman in [4], and has been used in a variety
of applications. For a short list, consult with [5]. Actually, the function
$\Phi$ considered there is rather simple:

$\Phi({\boldsymbol{x}})=\sum_{i=1}^{n}\phi(x_{i})$

with $\phi:{\mathcal{I}}\to{\mathbb{R}}$ being a convex function defined on
the interval ${\mathcal{I}}$ and ${\mathcal{M}}={\mathcal{I}}^{n}.$ This is a
typical example in many applications. The thrust in [5] was to compute the
geodesic distance defined by the Riemannian metric defined by the Hessian of
$\Phi,$ and to compare it to the pseudo distance $\delta_{\Phi}$ defined in
(6.1).

In this section we generalize a comparison result previously obtained in [5]
for the separable case. Before we state the result, we shall present introduce
some notations to unclutter the typography and carry out a few elementary
calculations that form the basis of the proof of the result.

We shall use $\partial_{i}$ $\partial^{2}_{i,j}$ to denote partial (and
repeated partial) derivatives with respect to $x_{i}$ (respectively $x_{i}$
and $x_{j}$). We shall use Einstein summation convention. That is, for example
$\partial_{i}U^{m}\partial_{j}U^{m}$ stands for
$\sum_{m=1}^{n}\partial_{i}U^{m}\partial_{j}U^{m}.$ Next we present three
instances of the same computation: Once for functions, once for vector fields
and once for matrix valued functions. Just to refer to them when the time
comes up.

Again let ${\mathcal{M}}\subset{\mathbb{R}}^{n}$ be an open, convex, connected
set and $\Phi:{\mathcal{M}}\to{\mathbb{R}}$ be a twice continuously
differentiable function. Let $\gamma:[0,1]\to{\mathcal{M}}$ be a continuously
differentiable curve. We shall denote by $\gamma_{k}$ the $k-$th component of
$\gamma.$ Let ${\boldsymbol{y}},{\boldsymbol{x}}\in{\mathbb{R}}^{n}$ and
suppose that $\gamma(0)={\boldsymbol{x}}$ and $\gamma(1)={\boldsymbol{y}}.$
Below, anytime that $\gamma(t)$ appears as the argument, say of a function
defined on ${\mathcal{M}},$ we shorten it to $t.$

Starting from

$\Phi({\boldsymbol{y}})=\Phi({\boldsymbol{x}})+\int_{0}^{1}\partial_{i}\Phi(s)\dot{\gamma}_{i}(t)dt.$

As above $\dot{\gamma}$ stands for the time derivative of $\gamma.$ Apply this
same computation to the function $\partial_{i}\Phi({\boldsymbol{x}})$ that
appears under the integral sign to obtain.

$\Phi({\boldsymbol{y}})=\Phi({\boldsymbol{x}})+\int_{0}^{1}\Big{(}\partial_{i}\Phi(0)+\int_{0}^{t}\partial^{2}_{k,i}\Phi(s)\dot{\gamma}_{k}(s)dt\Big{)}\dot{\gamma}_{i}(t)dt.$

Now integrate the first term to obtain
$({\boldsymbol{y}}-{\boldsymbol{x}})_{i}\partial_{i}\Phi({\boldsymbol{x}}).$
To complete, exchange the integration over $t$ with that over $s,$ notice that
the integral over $t$ becomes $({\boldsymbol{y}}-\gamma(s))_{i}$ and the whole
identity becomes

$\Phi({\boldsymbol{y}})=\Phi({\boldsymbol{x}})+({\boldsymbol{y}}-{\boldsymbol{x}})_{i}\partial_{i}\Phi({\boldsymbol{x}})+\int_{0}^{1}\big{(}{\boldsymbol{y}}-\gamma(s)\big{)}_{i}\partial^{2}_{k,i}\Phi(s)\dot{\gamma}_{k}(s)ds.$
(6.2)

A similar identity holds componentwise for vector valued functions
${\boldsymbol{A}}:{\mathcal{M}}\to{\mathbb{R}}^{n}.$ In this case the analogue
of (6.2) becomes

$A_{m}({\boldsymbol{y}})=A_{m}({\boldsymbol{x}})+({\boldsymbol{y}}-{\boldsymbol{x}})_{i}\partial_{i}A_{m}({\boldsymbol{x}})+\int_{0}^{1}\big{(}{\boldsymbol{y}}-\gamma(u)\big{)}_{i}\partial^{2}_{k,i}A_{m}(u)\dot{\gamma}_{k}(u)du.$
(6.3)

Let us now rewrite (6.2) as

$\delta_{\Phi}^{2}({\boldsymbol{y}},{\boldsymbol{x}})=\Phi({\boldsymbol{y}})-\Phi({\boldsymbol{x}})+({\boldsymbol{y}}-{\boldsymbol{x}})_{i}\partial_{i}\Phi({\boldsymbol{x}})=\int_{0}^{1}\big{(}{\boldsymbol{y}}-\gamma(s)\big{)}_{i}\partial^{2}_{k,i}\Phi(s)\dot{\gamma}_{k}(s)ds$
(6.4)

and proceed with the right hand side as follows. Suppose that as in Section 3
that there exists a diffeomorphism ${\boldsymbol{U}}$ such that
$\partial^{2}_{i,k}\Phi=\partial_{i}U^{m}\partial_{k}U^{m}.$ Now, for a fixed
$m$ consider only $({\boldsymbol{y}}-\gamma(s)\big{)}_{i}\partial_{i}U^{m}.$
According to (6.3) in which ${\boldsymbol{x}}$ is replaced by $\gamma(s),$
this can be rewritten as

$(U^{m}({\boldsymbol{y}})-U^{m}(s))-\int_{s}^{1}({\boldsymbol{y}}-\gamma(u)\big{)}_{i}\partial^{2}_{i,j}U^{m}(u)\dot{\gamma}_{j}(u)du.$

and reinserted back in (6.5) to obtain

$\delta_{\Phi}^{2}({\boldsymbol{y}},{\boldsymbol{x}})=\int_{0}^{1}\Big{(}(U^{m}({\boldsymbol{y}})-U^{m}(s))-\int_{s}^{1}({\boldsymbol{y}}-\gamma(u)\big{)}_{i}\partial^{2}_{i,j}U^{m}(u)\dot{\gamma}_{j}(u)du\Big{)}\partial_{k}U^{m}\dot{\gamma}_{k}(s)ds.$

Therefore, the first term under the outer integral becomes

$\int_{0}^{1}\big{(}U^{m}({\boldsymbol{y}})-\gamma(s)\big{)}\partial_{k}U^{m}\dot{\gamma}_{k}(s)ds=\frac{1}{2}\|U({\boldsymbol{y}})-U({\boldsymbol{x}})\|^{2}=\frac{1}{2}d_{U}^{2}({\boldsymbol{x}},{\boldsymbol{y}})$
(6.5)

To rewrite the second term, notice that
$\partial_{k}U^{m}\dot{\gamma}_{k}(s)ds=dU^{m}(\gamma(s)).$ Now, we shall
consider a specific trajectory:
$\gamma(t)={\boldsymbol{x}}+t({\boldsymbol{y}}-{\boldsymbol{x}}).$ With this
choice we have

$({\boldsymbol{y}}-\gamma(u)\big{)}_{i}\partial^{2}_{i,j}U^{m}(u)\dot{\gamma}_{j}(u)=(1-u)({\boldsymbol{y}}-{\boldsymbol{x}})_{i}U^{m}_{i,j}(u)({\boldsymbol{y}}-{\boldsymbol{x}})_{j}$

and therefore, the second integral becomes

$\int_{0}^{1}\Big{(}\int_{s}^{1}(1-u)\big{(}({\boldsymbol{y}}-{\boldsymbol{x}})_{i}\partial^{2}_{i,j}U^{m}(u)({\boldsymbol{y}}-{\boldsymbol{x}})_{j}\big{)}du\big{)}dU^{m}.$
(6.6)

We now gather there results in the main result of this section.

###### Theorem 6.1

With the notations introduced above, suppose that the convex function
$\Phi:{\mathcal{M}}\to{\mathbb{R}}$ is at least three times continuously
differentiable, that its Hessian can be factored as
$\partial^{2}_{i,k}\Phi=\partial_{i}U^{m}\partial_{k}U^{m},$ and that the
diffeomorphism ${\boldsymbol{U}}$ is at least twice continuously
differentiable and such that the sign of
$K({\boldsymbol{\xi}})\equiv\partial^{2}_{i,j}U^{m}({\boldsymbol{\xi}})\partial_{k}U^{m}({\boldsymbol{\xi}})$
is constant over ${\boldsymbol{\xi}}\in{\mathcal{M}},$ then

$\displaystyle\delta_{\Phi}^{2}({\boldsymbol{y}},{\boldsymbol{x}})\leq\frac{1}{2}d_{U}^{2}({\boldsymbol{x}},{\boldsymbol{y}})\;\;\;\mbox{whenever}\;\;K({\boldsymbol{\xi}})\geq
0,$ (6.7)
$\displaystyle\delta_{\Phi}^{2}({\boldsymbol{y}},{\boldsymbol{x}})\geq\frac{1}{2}d_{U}^{2}({\boldsymbol{x}},{\boldsymbol{y}})\;\;\;\mbox{whenever}\;\;K({\boldsymbol{\xi}})\leq
0.$ (6.8)

The proof is contained in the computations carried out above. The conclusion
is obtained after substituting (6.5) and (6.6) into (6.2).

In the most common case in applications, when
$\Phi({\boldsymbol{x}})=\sum_{k=1}^{n}\phi(x_{i}),$ that is in the separable
case, then $\phi^{\prime\prime}(x)=(u^{\prime}(x))^{2}>0$, then
$u^{\prime}(x)>0$ and the condition upon the sign of $u^{\prime\prime}(x)$ is
equivalent to a condition upon the sign of $\phi^{\prime\prime\prime}(x).$ In
this case, the result of Theorem 7.5 was obtained in [5].

##  7 Examples

###  7.1 Extended generalized arithmetic means

Here is a direct extension of the notion of generalized arithmetic mean.

###### Theorem 7.1

Let us suppose that ${\boldsymbol{U}}({\mathcal{M}})$ is convex, and consider
a set $\\{{\boldsymbol{x}}(1),...,{\boldsymbol{x}}(M)\\}$ of points in $cM.$
There is a unique point ${\boldsymbol{c}}$ in ${\mathcal{M}}$ which minimizes
the $d_{U}$ distance to the set
$\\{{\boldsymbol{x}}(1),...,{\boldsymbol{x}}(M)\\}.$ It is given by

${\boldsymbol{c}}={\boldsymbol{U}}^{-1}\Big{(}\frac{1}{M}\sum_{m=1}^{M}{\boldsymbol{U}}({\boldsymbol{x}}(m)))\Big{)}.$
(7.1)

The proof is easy and follows the pattern of the simple one-dimensional case:
To find ${\boldsymbol{\xi}}\in{\boldsymbol{U}}({\mathcal{M}}))$ such that
$\sum_{m=1}^{M}d_{U}({\boldsymbol{x}}(m)),{\boldsymbol{\xi}})$ is easy and it
is given by
${\boldsymbol{\xi}}=\frac{1}{M}\sum_{m=1}^{M}{\boldsymbol{U}}({\boldsymbol{x}}(m)).$
Since ${\boldsymbol{\xi}}\in{\boldsymbol{U}}({\mathcal{M}}),$ let
${\boldsymbol{c}}={\boldsymbol{U}}({\boldsymbol{\xi}})$ and we are through.

###  7.2 Means defined by a flow

With the notations of Section 3 and the previous example in mind, let
${\boldsymbol{U}}(t,{\boldsymbol{x}})$ be the flow associated to the geodesics
determined by the diffeomorphism ${\boldsymbol{U}},$ and the solution to the
geodesic equation such that ${\boldsymbol{x}}(0)={\boldsymbol{x}}$ and
$\dot{{\boldsymbol{x}}}(0)={\boldsymbol{\xi}}.$ That is

${\boldsymbol{U}}(t,{\boldsymbol{x}})={\boldsymbol{U}}^{-1}\big{(}{\boldsymbol{U}}({\boldsymbol{x}})+t{\boldsymbol{\xi}}\big{)}.$
(7.2)

To verify that for any real $t,s$
${\boldsymbol{U}}(t+s,{\boldsymbol{x}})={\boldsymbol{U}}(s,{\boldsymbol{U}}(t,{\boldsymbol{x}}))$
is routine -as long as the solution to the geodesic equations is defined for
all times. With the aid of (7.2) one can construct a family of transition
kernels $\\{P_{t}({\boldsymbol{x}},A):t\geq
0,A\in{\mathcal{B}}({\mathbb{R}}^{n})\\},$ where
${\mathcal{B}}({\mathbb{R}}^{n})$ denotes the Borel-subsets of
${\mathbb{R}}^{n},$ as follows:

$P_{t}({\boldsymbol{x}},A)=I_{A}({\boldsymbol{U}}(t,{\boldsymbol{x}}))$ (7.3)

where $I_{A}$ stands for the usual indicator function of the set $A.$ To
verify that $\\{P_{t}\\}$ is indeed a (semi)group is trivial using (7.2).
Skipping a considerable amount of detail, it is intuitive that this semi group
defines a Markov process $\\{{\boldsymbol{X}}_{t}:t\geq 0\\}$ having
${\mathbb{R}}^{n}$ as state space, and $\\{P_{t}\\}$ as transition semi group.
Furthermore, this process is such that for any Borel, measurable function
$f:{\mathbb{R}}^{n}\to{\mathbb{R}},$ we have

$E[f\big{(}{\boldsymbol{X}}_{t+s}\big{)}|{\boldsymbol{X}}_{t}]=f\big{(}{\boldsymbol{U}}(s,{\boldsymbol{X}}_{t})\big{)}.$

In particular, if $X_{i}(t)$ denotes de $i-$th coordinate of
${\boldsymbol{X}}_{t}$ for $i=1,...,n$ we have
$E[X_{i}(t+s)|{\boldsymbol{X}}_{t}]=U_{i}(s,{\boldsymbol{X}}_{t}),$ or in
vector notation:

$E[{\boldsymbol{X}}(t+s)|{\boldsymbol{X}}_{t}]={\boldsymbol{U}}(s,{\boldsymbol{X}}_{t}),$

that is, the current position is the best predictor of the future values of
the position.

###  7.3 Harmonic means in ${\mathbb{R}}^{n}$

As second example, consider the inversion with respect to the unit sphere in
${\mathcal{M}}={\mathbb{R}}^{n}\setminus\\{0\\}$ The mapping
${\boldsymbol{U}}({\boldsymbol{x}})={\boldsymbol{x}}/\|{\boldsymbol{x}}\|^{2}$
is an involution of ${\mathcal{M}}.$ If
$\\{{\boldsymbol{x}}(1),...,{\boldsymbol{x}}(M)\\}$ is a finite set of points
in $cM,$ then

${\boldsymbol{c}}={\boldsymbol{U}}\Big{(}\frac{1}{M}\sum_{m=1}^{M}\frac{{\boldsymbol{x}}(m)}{\|{\boldsymbol{x}}(m)\|^{2}}\Big{)}$

is the point closest to $\\{{\boldsymbol{x}}(1),...,{\boldsymbol{x}}(M)\\}$ is
the $d_{U}$ distance. Thus we have a variational interpretation of the
$n-$dimensional extension of the notion of harmonic mean.

###  7.4 Multidimensional best predictors

Here we extend the situation considered in the second example. Let
$(\Omega,{\mathcal{F}},{\mathbb{P}})$ be a probability space, that is, a set
$\Omega,$ a $\sigma-$algebra of subsets of $\Omega,$ and a probability
${\mathbb{P}}.$ By the customary $E_{\mathbb{P}}$ we shall denote expectation
with respect to ${\mathbb{P}}.$ All the random variables ${\boldsymbol{X}}$ in
this section will take values in ${\boldsymbol{U}}({\mathcal{M}})$ and be such
that $E_{\mathbb{P}}[{\boldsymbol{U}}({\boldsymbol{X}})^{2}]<\infty.$ Let
${\mathcal{G}}\subset{\mathcal{F}}$ be a sub-$\sigma-$algebra. We have

###### Theorem 7.2

Define the ${\boldsymbol{U}}$ distance between any two random variables
${\boldsymbol{X}},{\boldsymbol{Y}}$ by

$\delta_{U}({\boldsymbol{X}},{\boldsymbol{Y}})=\Big{(}E_{\mathbb{P}}[d_{U}({\boldsymbol{X}},{\boldsymbol{Y}})^{2}]\Big{)}^{1/2}.$

Then, there is a unique, ${\mathcal{G}}-$measurable, square integrable random
variable ${\boldsymbol{X}}_{\mathcal{G}}$ which satisfies

${\boldsymbol{X}}_{\mathcal{G}}=argmin\\{\delta_{U}({\boldsymbol{X}},{\boldsymbol{Y}})|{\boldsymbol{Y}}\;\;\mbox{measurable
with respect to}\;\;{\mathcal{G}}\\}.$

It is given by

${\boldsymbol{X}}_{\mathcal{G}}={\boldsymbol{U}}^{-1}\Big{(}E_{\mathbb{P}}[{\boldsymbol{U}}(X)|{\mathcal{G}}]\Big{)}.$

This result was obtained in [6]. The aim there was to present generalized
arithmetic means as best predictors, that is, as the solutions to the
variational problem stated in Theorem 7.2. A particular case of this result is
contained in the following result.

###### Theorem 7.3

Let ${\mathcal{G}}=\\{\emptyset,\Omega\\}$ be the trivial sigma algebra, and
Let ${\boldsymbol{X}}$ be a random variable taking finitely many values
$\\{{\boldsymbol{x}}_{1},...,{\boldsymbol{x}}_{M}\\}$ with probabilities
$P({\boldsymbol{X}}={\boldsymbol{x}}_{i})=p_{i},i=1,...,M.$ The rest of the
notations are as before. The best predictor of $X$ given no information is
given by the generalized arithmetic mean. As noted at the outset, it happens
to coincide with the notion of certainty equivalent.

$\langle{\boldsymbol{X}}\rangle={\boldsymbol{U}}^{-1}\Big{(}E_{\mathbb{P}}[{\boldsymbol{U}}(X)]\Big{)}={\boldsymbol{U}}^{-1}\Big{(}\sum_{k=1}^{M}{\boldsymbol{U}}({\boldsymbol{x}}_{i})P_{i}\Big{)}.$
(7.4)

###  7.5 Generalized arithmetic means and convex functions defined by the
gradient of a convex function

A particular case of the results in sections 3 and 5 is provided by strictly
convex (or concave) functions, that is functions whose Hessian is strictly
positive (or negative). Suppose $\Psi({\boldsymbol{x}})$ is a strictly
positive function on ${\mathbb{R}}^{n},$ then $\nabla\Psi$ is locally
invertible. Let us suppose that its inverse is global. If we consider the
metric introduced in (3.1), that is,

$g_{i,j}({\boldsymbol{x}})=\Psi_{i,k}({\boldsymbol{x}})\Psi_{j,k}({\boldsymbol{x}})$

using Einstein’s summation convention, then the generalized mean of the set
$\\{{\boldsymbol{x}}_{1},...{\boldsymbol{x}}_{M}\\}$ that it defines is given
by

$\big{(}\nabla\Psi\big{)}^{-1}\Big{(}\frac{1}{M}\sum_{k=1}^{M}\nabla\Psi({\boldsymbol{x}}_{k})\Big{)}.$

If the condition mentioned in Section 5.2 hold, namely:

$\mbox{For each
fixed}\;\;i=1,...,n\;\;\Psi_{i,m,k}\Psi_{j,k}=\Psi_{i,j,k}\Psi_{m,k}$

then there exists a convex function $\Phi({\boldsymbol{x}})$ on
${\mathbb{R}}^{n}$ such that

$\Phi_{i,j}({\boldsymbol{x}})=\Psi_{i,k}({\boldsymbol{x}})\Psi_{j,k}({\boldsymbol{x}})$

to which the comparison results established in Section 6 apply, whenever the
condition mentioned in Theorem 7.5 holds, namely:

$\mbox{The sign
of}\;\;\;\Psi_{i,j,m}({\boldsymbol{x}})\Psi_{k,m}({\boldsymbol{x}})\;\;\mbox{is
constant over}\;\;{\mathbb{R}}^{n}\;\;\mbox{independently of}\;\;i,j,k\;\;.$

###  7.6 Geodesics and distances determined by the Fenchel-Legendre conjugate

Consider the simplest possible case in which
$\Phi:{\mathbb{R}}^{n}\to{\mathbb{R}}$ has a strictly positive Hessian and
that $\nabla\Phi:{\mathbb{R}}^{n}\to{\mathbb{R}}^{n}$ is a diffeomorphism
whose range is ${\mathbb{R}}^{n}.$ Recall that:

###### Definition 7.1

The Fenchel-Lagrange dual $\Phi^{*}$ is defined by

$\Phi^{*}({\boldsymbol{\xi}})=\sup\\{\langle{\boldsymbol{\xi}},{\boldsymbol{x}}\rangle-\Phi({\boldsymbol{x}})|{\boldsymbol{x}}\in{\mathbb{R}}^{n}\\}.$

The basic properties of this definition can be seen, for example, in Borwein
and Lewis [3]. To begin with, we have:

###### Lemma 7.1

With the assumptions made at the outset of the section, a simple computation
shows that:

$\Phi({\boldsymbol{\xi}})=\langle{\boldsymbol{\xi}},\big{(}\nabla\Phi\big{)}^{-1}({\boldsymbol{\xi}})\rangle-\Phi\big{(}\big{(}\nabla\Phi\big{)}^{-1}({\boldsymbol{\xi}})\big{)}.$

And more importantly

$\nabla\Phi^{*}({\boldsymbol{\xi}})=\big{(}\nabla\Phi\big{)}^{-1}({\boldsymbol{\xi}}).$
(7.5)

And also (using Einstein’s summation convention):

$\partial^{2}_{i,k}\phi^{*}({\boldsymbol{\xi}})\partial^{2}_{k,j}\Phi\big{(}\big{(}\nabla\Phi\big{)}^{-1}({\boldsymbol{\xi}})\big{)}=\delta_{i,j}.$
(7.6)

The following result is known. See [9] for example. Its proof drops out of a
computation.

###### Theorem 7.4

With the notations introduced above, let $\\{{\boldsymbol{x}}(t):0\leq t\leq
1\\}$ be a geodesic in ${\mathbb{R}}^{n}$ between the points
${\boldsymbol{x}}_{1}={\boldsymbol{x}}(0)$ and
${\boldsymbol{x}}_{2}={\boldsymbol{x}}(1)$ with respect to the metric given by
the Hessian matrix of $\Phi.$ Let us put
${\boldsymbol{\xi}}(t)=\big{(}\nabla\Phi\big{)}({\boldsymbol{x}}(t))$ with
${\boldsymbol{\xi}}_{1}=\big{(}\nabla\Phi\big{)}({\boldsymbol{x}}_{1})$ and
${\boldsymbol{\xi}}_{2}=\big{(}\nabla\Phi\big{)}({\boldsymbol{x}}_{2}).$ Then
${\boldsymbol{\xi}}(t)$ is a geodesic between ${\boldsymbol{x}}_{1}$ and
${\boldsymbol{\xi}}_{2}$ in the metric given by the Hessian of $\Phi^{*}.$

Proof To verify the assertion it suffices invoke (7.5) and (7.6) to verify
that

$\partial^{2}_{i,j}\phi^{*}({\boldsymbol{\xi}})\dot{{\boldsymbol{\xi}}}_{i}(t)\dot{{\boldsymbol{\xi}}}_{j}(t)=\partial^{2}_{i,j}\phi({\boldsymbol{x}})\dot{{\boldsymbol{x}}}_{i}(t)\dot{{\boldsymbol{x}}}_{j}(t).$

This is left for the reader to carry out. $\Box$

Let us now verify that if the Hessian of $\Phi$ can be factored as in Section
3, then the Hessian of $\Phi^{*}$ can be factored as well.

###### Theorem 7.5

Let $\Phi$ and $\Phi^{*}$ be as above. Suppose that there is a diffeomorphism
${\boldsymbol{U}}$ of ${\mathbb{R}}^{n}$ such that
$\partial^{2}_{ij}\Phi({\boldsymbol{x}})=\partial_{i}U_{k}({\boldsymbol{x}})\partial_{j}U_{k}({\boldsymbol{x}}).$
Then
${\boldsymbol{U}}^{*}({\boldsymbol{\xi}})={\boldsymbol{U}}(\big{(}\nabla\Phi\big{)}^{-1}({\boldsymbol{\xi}}))$
is a factorization of the Hessian of $\Phi^{*}.$

Proof Observe that invoking (7.6) we obtain:

$\partial_{j}U^{*}_{k}({\boldsymbol{\xi}})=\partial_{l}U_{k}\big{(}\big{(}\nabla\Phi\big{)}^{-1}({\boldsymbol{\xi}})\big{)}\partial^{2}_{n}i\Phi*({\boldsymbol{\xi}}).$

From this and (7.6) it follows that

$\partial_{i}U^{*}_{k}({\boldsymbol{\xi}})\partial_{j}U^{*}_{k}({\boldsymbol{\xi}})=\partial_{i,j}\Phi^{*}({\boldsymbol{\xi}}).$

Thus concludes the proof. $\Box$
To finish, let us verify that the distances along the geodesics in the
Hessians of $\Phi^{*}$ and $\Phi$ coincide.

###### Corollary 7.1

With the notations introduced above, let ${\boldsymbol{\xi}}(t)$ and
${\boldsymbol{x}}(t)$) be geodesics described in Theorem 7.5. Then

$\delta_{\Phi^{*}}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2})=\delta_{\Phi}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2}).$

Proof We saw in Section 3 that
$\delta^{2}_{\Phi^{*}}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2})=\|{\boldsymbol{U}}^{*}({\boldsymbol{\xi}}_{1})-{\boldsymbol{U}}^{*}({\boldsymbol{\xi}}_{2})\|^{2}.$
From the definition of ${\boldsymbol{U}}^{*}$ and since
${\boldsymbol{\xi}}(t)=\nabla\Phi({\boldsymbol{x}})$ we obtain

$\delta^{2}_{\Phi^{*}}({\boldsymbol{\xi}}_{1},{\boldsymbol{\xi}}_{2})=\|{\boldsymbol{U}}^{*}({\boldsymbol{\xi}}_{1})-{\boldsymbol{U}}^{*}({\boldsymbol{\xi}}_{2})\|^{2}=\|{\boldsymbol{U}}({\boldsymbol{x}}_{1})-{\boldsymbol{U}}({\boldsymbol{x}}_{2})\|^{2}=\delta^{2}_{\Phi}({\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2}).$

And thus the assertion is verified. $\Box$

Acknowledgments: I would like to thank Soumalya Mukhopadhyay for bringing
references [2] up to my attention.

## References

  * [1] Arnold, V.I. (1978). Mathematical Methods of Classical Mechanics, Springer, New York.
  * [2] Berger, R.L. and Casella, G. (1992). Deriving Generalized Means as Least Squares and Maximum Likelihood Estimates, The American Statistician, 46, 279-282.
  * [3] Borwein, J. M. and Lewis, A.S. (2000). Convex Analysis and Non linear Optimization, CMS Books in Mathematics, Springer, New York.
  * [4] Bregman, L. (1967). The relaxation method of finding common points of convex sets and its application to the solution of problems in convex programming, Comp. Math. Phys., USSR, 7, 2100-217.
  * [5] Gzyl, H.(2020). Prediction in Riemannian metrics derived from divergence functions, Comm. in Stats., https://doi.org/10.1080/03610926.2020.1752384.
  * [6] Gzyl, H. (2020). A predictive approach to generalized arithmetic means, http://arxiv.org/abs/2005.07770.
  * [7] Hicks, N. (1965). Notes on Differential Geometry, Van Nostrand, Princeton,.
  * [8] Miller, J.D. (1984). Some Global Inverse Function Theorems, Journal of Mathematical Analysis and its Applications, 100, 375-384.
  * [9] Nielsen, F. (2018). An elementary introduction to information geometry, https://arxiv.org/abs/1808.08271
  * [10] Sternberg, S. (2012). Curvature in Mathematics and Physics, Dover Pubs. Inc., New York.