1 Introduction↩︎

Classical information geometry is agnostic to the geometry of the sample space, which is chosen to be \(\mathbb{R}^n\) in this article. This sample space is equipped with the standard inner product and the corresponding Euclidean distance function. The Kullback-Leibler divergence (KL-divergence), defined for probability measures on \(\mathbb{R}^n\), plays a central role in information geometry [1], [2]. It is, however, not coupled with the standard geometry of \(\mathbb{R}^n\). For instance, if we consider two Dirac measures concentrated in distinct points, then their KL-divergence is infinite and therefore not sensitive to the actual distance between these points. we derive an explicit closed-form expression for the WKL-divergence and study its geometric and continuity properties.

Based on previous work [3], [4], a general construction of a canonical divergence between two probability measures \(\mu\) and \(\nu\) was introduced. It is defined by the integral expression \[\begin{align} D(\mu \| \nu) = \int_0^1 t \lVert \dot{\gamma}(t) \rVert_{\gamma(t)}^2 dt, \end{align}\] where \(\gamma\) is the so-called e-geodesic (with respect to the exponential affine connection) connecting \(\mu\) and \(\nu\) and the norm is induced by a chosen Riemannian metric. The affine connection defining the e-geodesic is chosen to be the dual of the mixture connection where duality between the affine connections is again defined with respect to the chosen Riemannian metric. When the Fisher\(-\)Rao metric is used, the canonical divergence in 8 reduces to the KL-divergence with reversed order, i.e., \(D(\mu \| \nu)=D_{\rm KL}(\nu \| \mu)\). Building on this idea, a new canonical divergence has been proposed in [5], replacing the Fisher\(-\)Rao metric with the Otto metric [6], a Riemannian structure on the space of probability densities which underlies optimal transport theory. Since the Otto metric induces the Wasserstein distance as its Riemannian distance, the resulting divergence can be seen as the analogue of the KL-divergence within the Wasserstein geometric framework [7]. Since this represents the Wasserstein version of the classical KL-divergence, we therefore refer to it as the Wasserstein KL-divergence, abbreviated by \({\rm WKL}\).

2 Information Geometry of Gaussian Measures with the Otto Metric↩︎

In this section, we introduce the information geometric constructions leading to the definition of the Wasserstein KL divergence as the canonical divergence. We briefly review the preliminaries and notation used throughout the paper. We follow the development from [5] and apply these ideas to the setting of Gaussian probability measures on \(\mathbb{R}^n\). For further details, see [5]. On one hand, the results in the present paper can be seen as a special case for instantiation of the results from [5] to the case of finite dimensional Gaussian measures. On the other hand, [5] assumes a compact manifold structure which excludes the Gaussian measures defined on \(\mathbb{R}^n\).

2.1 Tangent and Cotangent Spaces of the Gaussian Manifold↩︎

Let \(\lambda\) denote the Lebesgue measure on \(\mathbb{R}^n\), and let \(C^\infty(\mathbb{R}^n)\) denote the space of smooth real-valued functions on \(\mathbb{R}^n\). Let the space of finite signed measures on \(\mathbb{R}^n\) be denoted by \(\mathcal{M}(\mathbb{R}^n)\) which is a Banach space endowed with the total variation norm. Let \[\begin{align} \mathcal{S}^\infty := \left\{ f\lambda \in \mathcal{M}(\mathbb{R}^n) \;\middle|\; f \in C^\infty(\mathbb{R}^n), \quad \int_{\mathbb{R}^n} |f(x)|\, d\lambda(x) < \infty \right\}. \end{align}\] be the space of finite signed measures on \(\mathbb{R}^n\) admitting a smooth density with respect to Lebesgue measure. We further define the linear subspace \[\begin{align} \mathcal{S}_0^\infty := \left\{ \mu \in \mathcal{S}^\infty \;\middle|\; \mu(\mathbb{R}^n)=0 \right\}. \end{align}\] Finally define \[\begin{align} \mathcal{P}_+^\infty := \left\{ \mu \in \mathcal{S}^\infty \;\middle|\; \mu(\mathbb{R}^n)=1, \quad \mu > 0 \right\}, \end{align}\] to be the space of smooth positive probability measures on \(\mathbb{R}^n\).

Let \(S_n^+\) denote the set of symmetric positive-definite matrices in \(\mathbb{R}^{n\times n}\). Define the parametrization \[\begin{align} \varphi : \mathbb{R}^n \times S_n^+ \longrightarrow \mathcal{P}_+^\infty, \quad (m,\Sigma) \longmapsto p(\;\cdot\;;m,\Sigma)\lambda, \end{align}\] where \[\begin{align} p(x;m,\Sigma) := \frac{1}{(2\pi)^{n/2}\det(\Sigma)^{1/2}} \exp\left( -\frac{1}{2}(x-m)^\top \Sigma^{-1}(x-m) \right) \end{align}\] is the multivariate Gaussian density with mean \(m \in \mathbb{R}^n\) and covariance matrix \(\Sigma \in S_n^+\). The image of \(\varphi\) is denoted by \(\mathcal{G} := \varphi(\mathbb{R}^n \times S_n^+) \subset \mathcal{P}_+^\infty,\) and is referred to as the manifold of Gaussian measures. We sometimes denote a Gaussian measure with mean \(m\) and covariance \(\Sigma\) by \(\mathcal{N}(m,\Sigma)\).

Let us now introduce the tangent space at a point \(\mu \in \mathcal{G}\). To this end, consider a smooth curve \((-\varepsilon,\varepsilon)\ni t \mapsto (m_t,\Sigma_t) \in \mathbb{R}^n \times S_n^+\) and the \(\varphi\) image of this curve in \(\mathcal{G}\) given by \[\begin{align} (-\varepsilon,\varepsilon)\ni t \mapsto \mu_t := \varphi(m_t,\Sigma_t) = p(\;\cdot\; ;m_t,\Sigma_t) \lambda \in \mathcal{G}. \end{align}\] First note that \[\begin{align} &\frac{d}{dt} p(x;m_t,\Sigma_t) \\ &= p(x;m_t,\Sigma_t) \frac{d}{dt}\left(\ln p(x;m_t,\Sigma_t) \right)\\ &=p(x;m_t,\Sigma_t)\left(\frac{1}{2}(x-m_t)^\top \Sigma_t^{-1}\dot{\Sigma}_t\Sigma_t^{-1}(x-m_t) + \dot{m}_t^\top\Sigma_t^{-1}(x-m_t)-\frac{1}{2}\textrm{Tr}\left({\Sigma_t^{-1}\dot{\Sigma}_t}\right)\right)\\ &=p(x;m_t,\Sigma_t)\cdot g(x;m_t,\Sigma_t,\dot{m}_t,\dot{\Sigma}_t), \end{align}\] where the quadratic function \(g\) is defined for \((m,\Sigma)\in \mathbb{R}^n\times S_n^{+}\) and \((m_v,\Sigma_v)\in \mathbb{R}^n\times S_n\) \[\begin{align} g(x;m,\Sigma,m_v,\Sigma_v)=\frac{1}{2}(x-m)^\top \Sigma^{-1}\Sigma_v\Sigma^{-1}(x-m) + m_v^\top\Sigma^{-1}(x-m)-\frac{1}{2}\textrm{Tr}\left({\Sigma^{-1}\Sigma_v}\right). \end{align}\] Let us compute the velocity of the curve \(\mu_t\) at \(t=0\) as \[\begin{align} \frac{d}{dt}\Big|_{t=0} \mu_t = \lim_{t\rightarrow 0}\frac{1}{t} \left(\mu_t -\mu_0\right) &= \lim_{t\rightarrow 0}\frac{1}{t}\left(\big(p(\;\cdot\; ;m_t,\Sigma_t) - p(\;\cdot\; ; m_0,\Sigma_0) \big)\lambda\right)\\ &= \left(\frac{d}{dt}\Big|_{t=0} p(\;\cdot\; ;m_t,\Sigma_t)\right)\lambda \\ &=g(\;\cdot\; ; m_0,\Sigma_0,\dot{m}_0,\dot{\Sigma}_0) \mu_0, \end{align}\] where the derivative is understood in the norm topology induced by the ambient space \(\mathcal{M}(\mathbb{R}^n)\). It can be verified that \[\begin{align} \left(\frac{d}{dt}\Big|_{t=0} \mu_t\right)(\mathbb{R}^n)=\int_{\mathbb{R}^n}g(\;\cdot\; ; m_0,\Sigma_0,\dot{m}_0,\dot{\Sigma}_0) \mu_0 dx=\mathbb{E}_{\mu_0}\left[g(\;\cdot\; ; m_0,\Sigma_0,\dot{m}_0,\dot{\Sigma}_0)\right]=0, \end{align}\] i.e., \(\left(\frac{d}{dt}\Big|_{t=0} \mu_t\right) \in \mathcal{S}_0^{\infty}\). Since a tangent vector at a point \(\mu \in \mathcal{G}\) is an equivalence classes of curves \(\gamma:(-\varepsilon,\varepsilon) \rightarrow \mathcal{G}\) with \(\gamma(0)=\mu\) with the equivalence relation being defined by equating velocities at \(t=0\), we can identify the tangent vectors with the velocities at \(t=0\). Thus, \(\frac{d}{dt}\Big|_{t=0} \mu_t \in \mathcal{S}_0^{\infty}\) is a tangent vector with the tangent space at a point \(\mu=\mathcal{N}(m,\Sigma) \in \mathcal{G}\) given by \[\begin{align} T_{\mu}\mathcal{G}=\left\{g(\;\cdot\; ; m,\Sigma,m_v,\Sigma_v)\mu \in \mathcal{S}_0^{\infty}\Big| (m_v,\Sigma_v) \in \mathbb{R}^n\times S_n\right\}. \end{align}\] Note that \(T_{\mu}\mathcal{G}\) is a finite dimensional subspace of \(\mathcal{S}_0^{\infty}\) of dimension \(d=n+\frac{n(n+1)}{2}\) which can also be identified with \(\mathbb{R}^n\times S_n\). Observe that \(g\) is a quadratic function of \(x\). With this motivation, define the set of quadratic functions \[\begin{align} \mathcal{Q}=\left\{f\in C^{\infty}(\mathbb{R}^n)\Big| f(x)=x^\top A x + b^\top x + c, A \in S_n, b\in \mathbb{R}^n, c \in \mathbb{R}\right\}. \end{align}\] Now consider an equivalence relation between two functions \(f,g\in\mathcal{Q}\), where \(f \sim g\) if \(f-g\) is the constant function. Let \(f+\mathbb{R}\) denote the equivalence class of functions differing from \(f\) by a constant function. Define \[\begin{align} \mathcal{Q}/\mathbb{R}=\left\{f+\mathbb{R}\;\big|\; f\in \mathcal{Q}\right\}= \left\{f+\mathbb{R} \;\Big|\; f(x)=x^\top A x + b^\top x, A \in S_n, b\in \mathbb{R}^n\right\}. \end{align}\] For a given measure \(\mu=\mathcal{N}(m,\Sigma) \in \mathcal{G}\), define \[\begin{align} \mathcal{Q}_{\mu}^0&=\left\{g(\;\cdot\; ;m,\Sigma,m_v,\Sigma_v)\in \mathcal{Q}\Big| (m_v,\Sigma_v)\in \mathbb{R}^n\times S_n\right\}. \end{align}\] Since for \(\mu=\mathcal{N}(m,\Sigma) \in \mathcal{G}\), \(\mathbb{E}_{\mu}\left[g(\;\cdot\; ;m,\Sigma,m_v,\Sigma_v)\right]=0\) for all \((m,\Sigma)\in \mathbb{R}^n\times S_n^{+}\) and \((m_v,\Sigma_v)\in \mathbb{R}^n\times S_n\), it can be verified that \[\begin{align} \mathcal{Q}_{\mu}^0&=\left\{f\in \mathcal{Q}\Big| \mathbb{E}_\mu[f]=0\right\}. \end{align}\] The tangent space at \(\mu \in \mathcal{G}\) thus has the compact description \[\begin{align} T_{\mu}\mathcal{G}= \left\{f \mu \in \mathcal{S}_0^{\infty}\Big| f \in \mathcal{Q}_{\mu}^0\right\}. \end{align}\]

Since tangent vectors at \(\mu \in \mathcal{G}\) are signed measures of the form \(a=f\mu\) with \(f \in \mathcal{Q}_\mu^0\), quadratic functions naturally define linear functionals on \(T_\mu\mathcal{G}\) through integration. Moreover, this pairing is well defined because for Gaussian \(\mu\), every polynomial function is \(\mu\)-integrable. Since \(T_\mu\mathcal{G}\) is finite-dimensional, every linear functional on it is continuous. Equivalently, if \(a=\sum_{j=1}^d a_j e_j\mu\) in a fixed basis, then \[\left|\int_{\mathbb{R}^n} f\,da\right| = \left|\sum_{j=1}^d a_j \int_{\mathbb{R}^n} f\,d(e_j\mu)\right| \le \sum_{j=1}^d |a_j| \left|\int_{\mathbb{R}^n} f e_j\,d\mu\right| \le C_f \|a\|_{TV},\] for some constant \(C_f>0\), where the last inequality uses equivalence of norms on the finite-dimensional space \(T_\mu\mathcal{G}\). This suggests the natural pairing \[(f,a)\mapsto \int_{\mathbb{R}^n} f\,da.\] Furthermore, since every tangent vector \(a\in T_\mu\mathcal{G}\) satisfies \(a(\mathbb{R}^n)=0\), constant functions act trivially on \(T_\mu\mathcal{G}\) under integration. Hence the pairing depends only on the class of a quadratic function modulo constants. Let \[\mathcal{Q}/\mathbb{R}:=\mathcal{Q}\big/\{c\mathbf{1}\mid c\in\mathbb{R}\}, \qquad [f]:=f+\mathbb{R},\] where \(\mathbf{1}\) denotes the constant function equal to \(1\). Then the pairing descends to \[\mathcal{Q}/\mathbb{R}\times T_\mu\mathcal{G}\ni (f+\mathbb{R},a)\longmapsto \int_{\mathbb{R}^n} f\,da,\] which is well defined because if \(f'=f+c\), then \[\int_{\mathbb{R}^n} f'\,da = \int_{\mathbb{R}^n} f\,da + c\,a(\mathbb{R}^n) = \int_{\mathbb{R}^n} f\,da.\] Moreover, every class \(f+\mathbb{R}\in \mathcal{Q}/\mathbb{R}\) has a unique representative in \(\mathcal{Q}_\mu^0\), namely \[f^0:=f-\mathbb{E}_\mu[f],\] so the canonical map \[\mathcal{Q}/\mathbb{R}\longrightarrow \mathcal{Q}_\mu^0,\qquad f+\mathbb{R}\longmapsto f-\mathbb{E}_\mu[f]\] is an isomorphism, and \(\mathcal{Q}_\mu^0\cong \mathcal{Q}/\mathbb{R}\).

Since \(\mathcal{Q}_\mu^0\cong \mathcal{Q}/\mathbb{R}\), showing that \(T_\mu^*\mathcal{G}\cong \mathcal{Q}_{\mu}^0\) completes the proof. Since \(\mathcal{Q}_\mu^0\) is finite dimensional with dimension \(d=n+\frac{n(n+1)}{2}\), there exists a basis \(\{e_1,\dots,e_d\},\) of \(\mathcal{Q}_{\mu}^0\) and therefore \(\{e_1\mu,\dots,e_d\mu\}\) is a basis of \(T_\mu\mathcal{G}\). In particular, every tangent vector \(a \in T_\mu\mathcal{G}\) can be written uniquely as \[a = \sum_{i=1}^d a_i\, e_i\mu, \qquad c_i \in \mathbb{R}.\] Let \(\ell\in (T_\mu\mathcal{G})^*\) be given. Define the coefficients \[l_j:=\ell(e_j\mu), \qquad j=1,\dots,d.\] We now seek \(f\in\mathcal{Q}_\mu^0\) of the form \[f=\sum_{i=1}^d f_i e_i\] such that \[\ell(a)=\int_{\mathbb{R}^n} f\, da \qquad \text{for all } a\in T_\mu\mathcal{G}.\] For each basis vector \(e_j\mu\), this requirement becomes \[l_j = \ell(e_j\mu) = \int_{\mathbb{R}^n} f\, d(e_j\mu) = \int_{\mathbb{R}^n} f e_j\, d\mu.\] Substituting \(f=\sum_i f_i e_i\), we obtain the linear system \[l_j=\sum_{i=1}^d f_i \int_{\mathbb{R}^n} e_i e_j\, d\mu, \qquad j=1,\dots,d.\] Equivalently, writing \[G_{ji}:=\int_{\mathbb{R}^n} e_i e_j\, d\mu,\] we have \[\begin{bmatrix} l_1 \\ \vdots \\ l_d \end{bmatrix} = \begin{bmatrix} G_{11} & \ldots & G_{1d}\\ \vdots & \ddots & \vdots \\ G_{d1} & \ldots & G_{dd} \end{bmatrix}\begin{bmatrix} f_1 \\ \vdots \\ f_d \end{bmatrix}\] where \(G\) is the Gram matrix of the basis \(\{e_i\}\). Since the basis elements are linearly independent, the Gram matrix \(G\) is invertible. Hence there exists a unique \[f=\sum_{i=1}^d f_i e_i\in \mathcal{Q}_\mu^0,\] such that \(\ell(e_j\mu) = \int_{\mathbb{R}^n} f\, d(e_j\mu)\). Now let \(a=\sum_{j=1}^d a_j e_j\mu\in T_\mu\mathcal{G}\). Then \[\begin{align} \int_{\mathbb{R}^n} f\, da &= \int_{\mathbb{R}^n} f \left(\sum_{j=1}^d a_j e_j\right) d\mu \\ &= \sum_{j=1}^d a_j \int_{\mathbb{R}^n} f e_j\, d\mu \\ &= \sum_{j=1}^d a_j \ell(e_j\mu) = \ell(a). \end{align}\] Thus \(\ell(a)=\int f\, da\) for all \(a\in T_\mu\mathcal{G}\).

We therefore identify the cotangent space as \[T_\mu^*\mathcal{G}\cong \mathcal{Q}/\mathbb{R}.\]

2.2 Affine Vector Fields and Their Induced Flows↩︎

Let \(\mathcal{T}(\mathbb{R}^n)\) denote the space of smooth vector fields on \(\mathbb{R}^n\), namely \[\begin{align} \mathcal{T}(\mathbb{R}^n) = \left\{ X : \mathbb{R}^n \to \mathbb{R}^n \;\middle|\; X \in C^\infty(\mathbb{R}^n,\mathbb{R}^n) \right\}. \end{align}\] In this work, we are particularly interested in the class of affine vector fields. A vector field \(X \in \mathcal{T}(\mathbb{R}^n)\) is called affine if it is of the form \(X(x) = Ax + b\), where \(A \in \mathbb{R}^{n \times n}\) and \(b \in \mathbb{R}^n\). The set of affine vector fields is denoted by \(\mathcal{T}_a(\mathbb{R}^n)\), which is naturally isomorphic to \(\mathbb{R}^{n\times n} \times \mathbb{R}^n\).

Given a vector field \(X \in \mathcal{T}(\mathbb{R}^n)\), let \(\varphi_t : \mathbb{R}^n \to \mathbb{R}^n\) denote its associated flow map, defined as the solution to the ordinary differential equation \[\frac{d}{dt}\varphi_t(x) = X(\varphi_t(x)), \qquad \varphi_0(x)=x. \label{eq:flow95equation}\tag{1}\] For affine vector fields, the flow \(\varphi_t\) forms a one-parameter family of affine transformations on \(\mathbb{R}^n\). The flow induces a natural action on measures through push-forward. Given a measure \(\mu\) on \(\mathbb{R}^n\), the push-forward measure \(\mu_t := (\varphi_t)_*(\mu)\) is defined by \[(\varphi_t)_*(\mu)(A) = \mu(\varphi_t^{-1}(A)) \label{eq:pushforward95measure}\tag{2}\] for every measurable set \(A \subseteq \mathbb{R}^n\). Equivalently, for every \(\mu\)-integrable test function \(f : \mathbb{R}^n \to \mathbb{R}\), \[\begin{align} \label{eq:pushforward95prop95measures} \int_{\mathbb{R}^n} f(x)\, d\mu_t(x) = \int_{\mathbb{R}^n} f(\varphi_t(x))\, d\mu(x). \end{align}\tag{3}\]

Fix \(\mu=\mathcal{N}(m,\Sigma)\in\mathcal{G}\). For an affine vector field \[X(x)=Ax+b, \qquad A\in\mathbb{R}^{n\times n},\;\;b\in\mathbb{R}^n,\] let \(\varphi_t\) denote its flow and define the evolved measure \[\mu_t := (\varphi_t)_*\mu.\] Since \(\mu\) has a smooth strictly positive density with respect to Lebesgue measure, each \(\mu_t\) is again absolutely continuous with respect to \(\mu\). This allows us to define the \(\mu\)-divergence of \(X\) by \[\begin{align} \label{eq:div95mu95defn} \operatorname{div}_\mu(X) :=- \left.\frac{d}{dt}\right|_{t=0}\frac{d\mu_t}{d\mu}. \end{align}\tag{4}\] In particular, \(\operatorname{div}_\mu(X)\) is a smooth function on \(\mathbb{R}^n\), and it satisfies \[\int_{\mathbb{R}^n} \operatorname{div}_\mu(X)\, d\mu = 0,\] so that \(\operatorname{div}_\mu(X)\in \mathcal{Q}_\mu^0\).

The flow of the affine vector field \(X(x)=Ax+b\) is given by \[\varphi_t(x)=e^{tA}x+\int_0^t e^{(t-s)A}b\,ds.\] Hence, if \(\mu=\mathcal{N}(m,\Sigma)\), then \(\mu_t=(\varphi_t)_*\mu\) is again Gaussian, with \[\mu_t=\mathcal{N}(m_t,\Sigma_t),\qquad m_t=e^{tA}m+\int_0^t e^{(t-s)A}b\,ds,\qquad \Sigma_t=e^{tA}\Sigma e^{tA^\top}.\] Differentiating at \(t=0\) gives \[\dot{m}_0=Am+b, \qquad \dot{\Sigma}_0=A\Sigma+\Sigma A^\top.\] By the formula computed in the previous subsection for the derivative of a Gaussian density with respect to its parameters, \[\left.\frac{d}{dt}\right|_{t=0}\frac{d\mu_t}{d\mu} = \left.\frac{d}{dt}\right|_{t=0}\frac{p(\cdot;m_t,\Sigma_t)}{p(\cdot;m,\Sigma)} = g(\cdot\,;m,\Sigma,\dot{m}_0,\dot{\Sigma}_0).\] Therefore \[\operatorname{div}_\mu(X) = -g(\cdot\,;m,\Sigma,Am+b,A\Sigma+\Sigma A^\top).\] Since this is of the form \(-g(\cdot\,;m,\Sigma,m_v,\Sigma_v)\), it lies in \(\mathcal{Q}_\mu^0\). Accordingly, we obtain a linear map \[\operatorname{div}_\mu : \mathcal{T}_a(\mathbb{R}^n)\longrightarrow \mathcal{Q}_\mu^0, \qquad X\longmapsto \operatorname{div}_\mu(X),\] which associates to each affine vector field its induced centered quadratic score function. Note that \[A\Sigma+\Sigma A^\top \in S_n.\] Since \(\Sigma\in S_n^+\) is symmetric positive definite, the Lyapunov operator \[\mathcal{L}_\Sigma : S_n \to S_n, \qquad B \longmapsto B\Sigma+\Sigma B,\] is an isomorphism. Hence there exists a unique matrix \(A_{\mathrm{symm}}\in S_n\) such that \[A_{\mathrm{symm}}\Sigma+\Sigma A_{\mathrm{symm}}^\top = A\Sigma+\Sigma A^\top.\] Define \[A^\perp := A-A_{\mathrm{symm}}, \qquad b^\perp := -A^\perp m, \qquad b_{\mathrm{symm}} := b+A^\perp m.\] Then \[X_{\mathrm{grad}}(x):=A_{\mathrm{symm}}x+b_{\mathrm{symm}}, \qquad X_{\mathrm{divfree}}(x):=A^\perp x+b^\perp=A^\perp(x-m),\] so that \[X=X_{\mathrm{grad}}+X_{\mathrm{divfree}}.\] Moreover, one checks that \[\operatorname{div}_\mu(X_{\mathrm{divfree}})=0 \qquad \text{and} \qquad \operatorname{div}_\mu(X_{\mathrm{grad}})=\operatorname{div}_\mu(X).\]

This decomposition is in fact an orthogonal decomposition with respect to the \(L^2(\mu)\) inner product on vector fields, as summarized in the next proposition.

The earlier discussion already shows the unique existence of the vector fields \(X_{\mathrm{divfree}}\) and \(X_{\mathrm{grad}}\). So we only need to prove that \(\langle\!\langle X_{\mathrm{divfree}},X_{\mathrm{grad}}\rangle\!\rangle_{\mu}^{O} =0\). To this end, observe that \[\begin{align} \langle\!\langle X_{\mathrm{divfree}},X_{\mathrm{grad}}\rangle\!\rangle_{\mu}^{O}&=\mathbb{E}_{\mu}\!\left[ X_{\mathrm{divfree}}(x)^\top X_{\mathrm{grad}}(x) \right] \\ &=\mathbb{E}_{\mu}\!\left[\left(A^\perp(x-m)\right)^\top \left(A_{\mathrm{symm}}x+b_{\mathrm{symm}}\right)\right]\\ &=\mathbb{E}_{\mu}\!\left[(x-m)^\top \left(A^\perp\right)^\top A_{\mathrm{symm}}x\right]\\ &=\mathbb{E}_{\mu}\!\left[(x-m)^\top \left(A^\perp\right)^\top A_{\mathrm{symm}}\left(x-m\right)\right]\\ &=\textrm{Tr}\left({\left(A^\perp\right)^\top A_{\mathrm{symm}}\Sigma}\right)=-\textrm{Tr}\left({\left(A^\perp\right)^\top A_{\mathrm{symm}}\Sigma}\right)\\ &=0, \end{align}\] where we have used the cyclic property of the trace operator and \(A^\perp\Sigma+\Sigma \left(A^\perp\right)^\top=0\).

Consequently, it is enough to consider gradient vector fields, which can be written in the form \[X_{\mathrm{grad}}(x)=\nabla f(x), \qquad f\in \mathcal{Q}.\] Since two functions differing by a constant have the same gradient, we may restrict attention to \(\mathcal{Q}/\mathbb{R}\) and arrive at the following proposition.

Let \(f(x)=\frac{1}{2} x^\top A x+b^\top x\) with \(A\in S_n\). Then \[\nabla f(x)=Ax+b,\] and Proposition 2 gives \[\Delta_\mu(f+\mathbb{R}) = \operatorname{div}_\mu(\nabla f) = -g(\cdot\,;m,\Sigma,Am+b,A\Sigma+\Sigma A).\]

For surjectivity, let \(g(\cdot\,;m,\Sigma,m_v,\Sigma_v)\in\mathcal{Q}_\mu^0\) be arbitrary. Since \(\Sigma\in S_n^+\), the Lyapunov operator \[\mathcal{L}_\Sigma:S_n\to S_n,\qquad A\mapsto A\Sigma+\Sigma A\] is an isomorphism. Let \(A=\mathcal{L}_\Sigma^{-1}(\Sigma_v)\) and define \[b:=m_v-Am.\] Then with \[f(x)=\frac{1}{2} x^\top A x+b^\top x\] we obtain \[\Delta_\mu(f+\mathbb{R}) = -g(\cdot\,;m,\Sigma,m_v,\Sigma_v).\] Injectivity follows because if \(\Delta_\mu(f+\mathbb{R})=0\), then \(A=0\) and \(b=0\), so \(f\) is constant modulo \(\mathbb{R}\).

2.3 Riemannian Metrics on the Gaussian Manifold↩︎

2.3.1 Fisher–Rao Metric↩︎

We now introduce the Fisher–Rao metric by translating the \(L_2(\mu)\) product defined for functions from the cotangent space to the tangent space. Consider the \(L^2(\mu)\) product for functions \(f,g \in \mathcal{Q}_\mu^0\) defined by \[\langle\!\langle f,g\rangle\!\rangle_{\mu}^{FR} = \int_{\mathbb{R}^n} f(x)g(x)\,\mu(dx),\] which can be translated to functions in \((f+\mathbb{R}),(g+\mathbb{R}) \in \mathcal{Q}/\mathbb{R}\) as \[\begin{align} \langle\!\langle f+\mathbb{R},g+\mathbb{R}\rangle\!\rangle_{\mu}^{FR} = \int_{\mathbb{R}^n} \left(f(x)-\mathbb{E}_{\mu}[f(x)]\right)\left(g(x)-\mathbb{E}_{\mu}[g(x)]\right)\,\mu(dx). \end{align}\] Note that since \[\begin{align} \langle\!\langle f+\mathbb{R},g+\mathbb{R}\rangle\!\rangle_{\mu}^{FR}= \int_{\mathbb{R}^n} \left(f-\mathbb{E}_{\mu}[f]\right)\left(g-\mathbb{E}_{\mu}[g]\right)d\mu=\int_{\mathbb{R}^n} f da, \end{align}\] where \(a=\left(g-\mathbb{E}_{\mu}[g]\right)\mu\), we have that \(\langle\!\langle \;\cdot\;,g+\mathbb{R}\rangle\!\rangle_{\mu}^{FR}\) acts on \(f+\mathbb{R}\) in the same way as the measure \(a\) acts on \(f+\mathbb{R}\). We thus have the isomorphism \[\phi_\mu : \mathcal{Q}/\mathbb{R} \to T_\mu \mathcal{G}, \qquad f+\mathbb{R} \mapsto \left(f-\mathbb{E}_{\mu}[f]\right)\mu,\] with the inverse \[\phi_\mu^{-1} : T_\mu \mathcal{G} \to \mathcal{Q}/\mathbb{R}, \qquad a \mapsto \frac{da}{d\mu}+\mathbb{R}.\] We can use \(\phi_\mu\) to translate the FR inner product to the tangent space from the cotangent space: \[\begin{align} \langle\!\langle a,b\rangle\!\rangle_{\mu}^{FR} :&= \left\langle\!\left\langle \phi_{\mu}^{-1}(a), \phi_{\mu}^{-1}(b) \right\rangle\!\right\rangle_{\mu}^{FR} = \left\langle\!\left\langle \frac{da}{d\mu}+\mathbb{R}, \frac{db}{d\mu} +\mathbb{R}\right\rangle\!\right\rangle_{\mu}^{FR}\\ &= \int_{\mathbb{R}^n} \left(\frac{da}{d\mu}\right) \left(\frac{db}{d\mu}\right)\, d\mu\\ &= \mathbb{E}_\mu\!\left[\left(\frac{da}{d\mu}\right)\left(\frac{db}{d\mu}\right)\right]\\ \end{align}\]

2.3.2 Otto Metric↩︎

We now follow a similar recepie for the Otto metric by translating the \(L_2(\mu)\) product defined for gradient vector fields from the cotangent space to the tangent space. For \(f,g \in \mathcal{Q}/\mathbb{R}\), let \[\langle\!\langle f+\mathbb{R},g+\mathbb{R}\rangle\!\rangle_{\mu}^{O} = \int_{\mathbb{R}^n} \nabla f(x)^\top \nabla g(x)\,\mu(dx).\] Since \[\begin{align} {2} \langle\!\langle f+\mathbb{R},g+\mathbb{R}\rangle\!\rangle_{\mu}^{O} &= \langle\!\langle \nabla f,\nabla g\rangle\!\rangle_{\mu}^{O} = \int_{\mathbb{R}^n} \nabla f^{\top}\nabla g \, d\mu &\\ &= \int_{\mathbb{R}^n} \left.\frac{d}{dt}\right|_{t=0} (f\circ \varphi_t)\, d\mu &\text{(where } \varphi_t \text{ is the flow of } \nabla g)\\ &= \left.\frac{d}{dt}\right|_{t=0}\int_{\mathbb{R}^n} (f\circ \varphi_t)\, d\mu &\text{(because f\circ \varphi_t is quadratic)} \\ &= \left.\frac{d}{dt}\right|_{t=0}\int_{\mathbb{R}^n} f\, d\mu_t &\text{(because of \eqref{eq:pushforward95prop95measures})}\\ &= \left.\frac{d}{dt}\right|_{t=0}\int_{\mathbb{R}^n} f\,\left(\frac{d\mu_t}{d\mu}\right) d\mu & \\ &= \int_{\mathbb{R}^n} f\, \left.\frac{d}{dt}\right|_{t=0} \left(\frac{d\mu_t}{d\mu}\right)\, d\mu &\text{(because \frac{d}{dt}\Big|_{t=0} \left(\frac{d\mu_t}{d\mu}\right) is quadratic)}\\ &= \int_{\mathbb{R}^n} f\bigl(-\operatorname{div}_{\mu}(\nabla g)\bigr)\, d\mu & \text{(using definition \eqref{eq:div95mu95defn})}\\ &= \int_{\mathbb{R}^n} f\bigl(-\Delta_{\mu}(g+\mathbb{R})\bigr)\, d\mu & \text{(using definition \eqref{eq:Laplace95op95defn})}\\ &= \int_{\mathbb{R}^n} fda \qquad &\text{(where a= -\Delta_{\mu}(g+\mathbb{R})\mu)},\\ &= \int_{\mathbb{R}^n} (f+\mathbb{R})da \qquad &\text{since a\in \mathcal{S}_0^{\infty}}. \end{align}\] We have used the fact that polynomials of any degree are in \(L_1(\mu)\), which make all the integrals in the above chain of equalities well defined and justify moving \(\frac{d}{dt}\) in and out of the integral.

Thus, \(\langle\!\langle \;\cdot\;,g+\mathbb{R}\rangle\!\rangle_{\mu}^{O}\) acts on \(f+\mathbb{R}\) in the same way as the measure \(a\) acts on \(f+\mathbb{R}\). With this motivation, define the isomorphism \(\phi_{\mu}\) and its inverse as \[\begin{align} \phi_{\mu} : \mathcal{Q}/\mathbb{R} &\to T_{\mu}\mathcal{G}, & g+\mathbb{R} &\mapsto -\Delta_{\mu}(g+\mathbb{R})\,\mu, \\ \phi_{\mu}^{-1} : T_{\mu}\mathcal{G} &\to \mathcal{Q}/\mathbb{R}, & a &\mapsto -\Delta_{\mu}^{-1}\!\left(\frac{da}{d\mu}\right)+\mathbb{R}, \end{align}\] respectively. We can now translate the inner product from the cotangent to the tangent space via the isomorphisms as \[\begin{align} \langle\!\langle a,b\rangle\!\rangle_{\mu}^{O} :&=\left\langle\!\left\langle \phi_{\mu}^{-1}(a), \phi_{\mu}^{-1}(b) \right\rangle\!\right\rangle_{\mu}^{O}\\ &= \left\langle\!\left\langle -\Delta_{\mu}^{-1}\!\left(\frac{da}{d\mu}\right)+\mathbb{R}, -\Delta_{\mu}^{-1}\!\left(\frac{db}{d\mu}\right)+\mathbb{R} \right\rangle\!\right\rangle_{\mu}^{O} \\ &= \int_{\mathbb{R}^n} \nabla\!\left(\Delta_{\mu}^{-1}\!\left(\frac{da}{d\mu}\right)\right)^{\top} \nabla\!\left(\Delta_{\mu}^{-1}\!\left(\frac{db}{d\mu}\right)\right)\, d\mu \\ &= \int_{\mathbb{R}^n} (A_a x+b_a)^{\top}(A_b x+b_b)\, d\mu\\ &=\operatorname{Tr}(A_aA_b\Sigma) + (A_a m+b_a)^\top(A_b m+b_b), \end{align}\] where \(A_a\), \(b_a\), \(A_b\) and \(b_b\) are solutions to equations \[\begin{align} A_a\Sigma+\Sigma A_a^{\top}=\dot{\Sigma}_a,\quad A_a m+b_a&=\dot{m}_a,\\ A_b\Sigma+\Sigma A_b^{\top},=\dot{\Sigma}_b, \quad A_b m + b_b&=\dot{m}. \end{align}\]

2.4 Exponential Affine Connections Associated with the Fisher–Rao and Otto Metrics↩︎

In this subsection, we define the \(e_0\)- and \(e_1\)-connections as in [5]. The construction is the same in both cases: identify a tangent vector \(a\) with its cotangent representative \(f_a+\mathbb{R}=\phi_{\mu}^{-1}(a)\), transport that representative trivially as \[\begin{align} T_{\mu}^*\mathcal{G}\ni (\mu,f_a+\mathbb{R}) \mapsto (\nu,f_a+\mathbb{R}) \in T_{\nu}^*\mathcal{G}, \end{align}\] and map \(f_a+\mathbb{R}\) at \(\nu\) back to the tangent vector \(\phi_{\nu}(f_a+\mathbb{R})\) at \(\nu\).

2.4.1 The \(e_0\)-Connection for the Fisher–Rao Metric↩︎

Following the above described construction for the Fisher–Rao metric gives the parallel transport map \[\begin{align} \Pi_{\mu,\nu}^{(e_0)}: T_{\mu}\mathcal{G}\ni (\mu,a) \mapsto (\nu,(\phi_{\nu}\circ \phi_{\mu}^{-1})(a)) \in T_{\nu}\mathcal{G}, \end{align}\] where \[\begin{align} (\phi_{\nu}\circ \phi_{\mu}^{-1})(a)=\left(\frac{da}{d\mu}-\mathbb{E}_{\nu}\left[\frac{da}{d\mu}\right]\right)\nu. \end{align}\] This can be explicitly written for \(\mu=\mathcal{N}(m_{\mu},\Sigma_{\mu})\) and \(\nu=\mathcal{N}(m_{\nu},\Sigma_{\nu})\) in coordinates as \[\begin{align} \Pi_{\mu,\nu}^{(e_0)}: (\mu,g(x;m_{\mu},\Sigma_{\mu},m_v,\Sigma_v)\mu) \mapsto \Big(\nu,\big(g(x;m_{\mu},\Sigma_{\mu},m_v,\Sigma_v)-\mathbb{E}_{\nu}[g(x;m_{\mu},\Sigma_{\mu},m_v,\Sigma_v)]\big)\nu\Big). \end{align}\] The corresponding \(e_0-\)geodesic \(\gamma_0\) can is described by \[\begin{align} \label{eq:e0-geodesic} \dot{\gamma}_0(t)=\Pi_{\mu,\gamma_0(t)}^{(e_0)}a=\left(\frac{da}{d\mu}-\mathbb{E}_{\gamma_0(t)}\left[\frac{da}{d\mu}\right]\right)\gamma_0(t). \end{align}\tag{5}\] In order to write the geodesic equations in coordinates, consider the curve \(t\mapsto \mu_t\) as defined in Section 2.1 and let \((m_a,\Sigma_a)\in \mathbb{R}^n\times S_n\) represent the initial velocity. The geodesic equation leads to \[\begin{align} g(x;m_t,\Sigma_t,\dot{m}_t,\dot{\Sigma}_t)\cdot \mu_t=\left(g(x;m_0,\Sigma_0,m_a,\Sigma_a)-\mathbb{E}_{\mu_t}[g(x;m_0,\Sigma_0,m_a,\Sigma_a]\right)\mu_t. \end{align}\] Equating the quadratic and linear terms in \(x\) gives \[\begin{align} \Sigma_t^{-1}\dot{\Sigma}_t\Sigma_t^{-1}&=\Sigma_0^{-1}\Sigma_a\Sigma_0^{-1},\\ -\Sigma_t^{-1}\dot{\Sigma}_t\Sigma_t^{-1}m_t+\Sigma_t^{-1}\dot{m}_t &= -\Sigma_0^{-1}\Sigma_a\Sigma_0^{-1}m_0+\Sigma_0^{-1}m_a, \end{align}\] which can be rearranged to obtain \[\begin{align} \dot{\Sigma}_t&=\Sigma_t\Sigma_0^{-1}\Sigma_a\Sigma_0^{-1}\Sigma_t,\\ \dot{m}_t &= \Sigma_t(\Sigma_0^{-1}\Sigma_a\Sigma_0^{-1}m_t-\Sigma_0^{-1}\Sigma_a\Sigma_0^{-1}m_0)+\Sigma_t \Sigma_0^{-1}m_a. \end{align}\] These equations can be solved by transforming to the "natural" exponential family coordinates \(P_t=\Sigma_t^{-1}\) and \(\eta_t=\Sigma_t^{-1}m_t\) to get \[\begin{align} {2} \dot{P}_t&=-P_0\Sigma_a P_0 \quad &\implies& \quad P_t=P_0-t\cdot P_0\Sigma_aP_0,\\ \dot{\eta}_t&=P_0m_a-P_0\Sigma_a\eta_0 \quad &\implies& \quad \eta_t=\eta_0+t\cdot (P_0m_a-P_0\Sigma_a\eta_0), \end{align}\] which are straight-line trajectories in the natural coordinates.

If one starts instead with two points \(\mu,\nu\in \mathcal{G}\), then it is possible to solve for initial velocities \((m_a,\Sigma_a)\) such that the geodesics connect \(\mu_0=\mathcal{N}(m_0,\Sigma_0)\) to \(\mu_1=\mathcal{N}(m_1,\Sigma_1)\) in unit time. The resulting geodesic equations are \[\begin{align} m_t&=\Sigma_t\left((1-t)\Sigma_0^{-1}m_0+t\Sigma_1^{-1}m_1\right),\\ \Sigma_t&=\left((1-t)\Sigma_0^{-1}+t\Sigma_1^{-1}\right)^{-1}. \end{align}\] In coordinate free notation, the \(e_0\)-geodesic can be written as \[\begin{align} \gamma_0(t)=\frac{\left(\frac{d\nu}{d\mu}\right)^t}{Z(t)} \mu=\rho_t \cdot \mu, \textrm{ where }Z(t)=\int_{\mathbb{R}^n} \left(\frac{d\nu}{d\mu}\right)^t d\mu. \end{align}\]

2.4.2 The \(e_1\)-Connection for the Otto Metric↩︎

We now investigate what happens when we use the Otto metric instead of the Fisher–Rao metric.

The parallel transport map defined on the tangent space is \[\begin{align} \Pi_{\mu,\nu}^{(e_1)}: T_{\mu}\mathcal{G}\ni (\mu,a) \mapsto (\nu,(\phi_{\nu}\circ \phi_{\mu}^{-1})(a)) \in T_{\nu}\mathcal{G}, \end{align}\] where \[\begin{align} (\phi_{\nu}\circ \phi_{\mu}^{-1})(a)=\left(\left(\Delta_{\nu} \circ \Delta_{\mu}^{-1}\right)\left(\frac{da}{d\mu}\right)\right)\nu. \end{align}\] This can be explicitly computed for \(\mu=\mathcal{N}(m_{\mu},\Sigma_{\mu})\) and \(\nu=\mathcal{N}(m_{\nu},\Sigma_{\nu})\) in coordinates. To this end, let \(g(\;\cdot\; ;m_{\mu},\Sigma_{\mu},m_v,\Sigma_v)\in \mathcal{Q}_{\mu}^0\). Recall that \[\begin{align} \Delta_{\mu}^{-1}(g(\;\cdot\; ;m_{\mu},\Sigma_{\mu},m_v,\Sigma_v))=\underbrace{\left(x\mapsto \frac{1}{2}x^\top A_{\mu} x+b_{\mu}\right)}_{f}+\mathbb{R}, \end{align}\] where \(A_{\mu}\) is the unique symmetric solution to the Lyapunov equation \(A_{\mu}\Sigma_{\mu}+\Sigma_{\mu}A_{\mu}=\Sigma_v\) and \(b_{\mu}=m_v-A_{\mu}m_{\mu}\). Finally, \[\begin{align} \Delta_{\nu}(f+\mathbb{R})=g(\;\cdot\; ;m_{\nu},\Sigma_{\nu},A_{\mu}m_{\nu}+b_{\mu},A_{\mu}\Sigma_{\nu}+\Sigma_{\nu}A_{\mu}). \end{align}\] Therefore, the transport map in coordinates can be written as \[\begin{align} \Pi_{\mu,\nu}^{(e_1)}: (\mu,g(\;\cdot\; ;m_{\mu},\Sigma_{\mu},m_v,\Sigma_v)\mu) \mapsto \Big(\nu,g(\;\cdot\; ;m_{\nu},\Sigma_{\nu},A_{\mu}m_{\nu}+b_{\mu},A_{\mu}\Sigma_{\nu}+\Sigma_{\nu}A_{\mu})\nu\Big), \end{align}\] where \(A_{\mu},b_{\mu}\) satisfy \[\begin{align} \label{eq:A95mu95b95mu} A_{\mu}\Sigma_{\mu}+\Sigma_{\mu}A_{\mu}&=\Sigma_v, \\ A_{\mu}m_{\mu}+b_{\mu}&=m_v. \end{align}\tag{6}\]

The corresponding \(e_1-\)geodesic \(\gamma_1\) can is described by \[\begin{align} \label{eq:e1-geodesic} \dot{\gamma}_1(t)=\Pi_{\mu,\gamma_1(t)}^{(e_1)}a=\left(\left(\Delta_{\gamma_1(t)} \circ \Delta_{\mu}^{-1}\right)\left(\frac{da}{d\mu}\right)\right)\gamma_1(t). \end{align}\tag{7}\] In order to write the geodesic equations in coordinates, consider the curve \(t\mapsto \mu_t\) as defined in Section 2.1 and let \((m_a,\Sigma_a)\in \mathbb{R}^n\times S_n\) represent the initial velocity. The geodesic equation leads to \[\begin{align} g(x;m_t,\Sigma_t,\dot{m}_t,\dot{\Sigma}_t) \mu_t=\left(g(\;\cdot\; ;m_t,\Sigma_t,A_{\mu}m_t+b_{\mu},A_{\mu}\Sigma_t+\Sigma_tA_{\mu})\right)\mu_t. \end{align}\] Equating the quadratic and linear terms in \(x\) gives \[\begin{align} \dot{\Sigma}_t&=A_{\mu}\Sigma_t+\Sigma_tA_{\mu},\\ \dot{m}_t &= A_{\mu}m_t+b_{\mu}, \end{align}\] where \(A_{\mu}\) and \(b_{\mu}\) are given in 6 . These linear differential equations can be solved to obtain \[\begin{align} \Sigma_t&=e^{A_{\mu}t}\Sigma_0 e^{A_{\mu}t},\\ m_t&=e^{A_\mu t}m_0+\int_0^t e^{A_\mu(t-s)}b_\mu\,ds. \end{align}\] If one starts instead with two points \(\mu,\nu\in \mathcal{G}\), then it is possible to solve for initial velocities \((m_a,\Sigma_a)\) such that the geodesics connect \(\mu_0=\mathcal{N}(m_0,\Sigma_0)\) to \(\mu_1=\mathcal{N}(m_1,\Sigma_1)\) in unit time. This is analyzed later in Lemma 2.

2.5 Canonical Divergence Function↩︎

Recall the general construction of a canonical divergence alluded to in the introduction. The canonical divergence between two probability measures \(\mu\) and \(\nu\) is given by \[\begin{align} \label{eq:KL95energy95formula} D(\mu \| \nu) = \int_0^1 t \lVert \dot{\gamma}(t) \rVert_{\gamma(t)}^2 dt, \end{align}\tag{8}\] where \(\gamma\) is an e-geodesic (with respect to one of the two exponential affine connections discussed above) connecting \(\mu\) and \(\nu\) and the norm is induced by the corresponding chosen Riemannian metric.

For the \(e_{0}\)-geodesic given by 5 , the canonical divergence 8 reduces to the KL-divergence, i.e., \[\begin{align} \label{eq:KL95energy95formula95specialize95to95KL} D^{(\rm e_0)}(\mu \| \nu) &:=\int_0^1 t \langle\dot{\gamma}_0(t),\dot{\gamma}_0(t) \rangle_{\gamma_0(t)}^{\rm FR} dt= D_{\rm KL}(\nu \| \mu). \end{align}\tag{9}\] We show this in Appendix 5. For the \(e_{1}\)-geodesic given by 7 , the canonical divergence 8 reduces to \[\begin{align} D^{(\rm e_1)}(\mu \| \nu):=\int_0^1 t \langle\dot{\gamma}_1(t),\dot{\gamma}_1(t) \rangle_{\gamma_1(t)}^{\rm O} dt=\int_{\mathbb{R}^n} \int_0^1 \left(f \circ \varphi_1 - f \circ \varphi_t\right) \, dt \, d\mu. \end{align}\] This is shown in Appendix 7. This defines the Wasserstein analogue of the KL-divergence, which we call the Wasserstein KL-divergence, abbreviated by \({\rm WKL}\). \[\label{eq:contrast95func} D_{\rm WKL}(\mu \| \nu) \; := \; \int_{\mathbb{R}^n} \int_0^1 \left(f \circ \varphi_1 - f \circ \varphi_t\right) \, dt \, d\mu.\tag{10}\]

We verify that \(D_{\rm WKL}\) satisfies the usual properties of a divergence. Since the dynamics are governed by a gradient flow, observe that \[\begin{align} \frac{d}{dt} f(\varphi_t(x))=\langle \nabla{f}(\varphi_t(x)), \frac{d}{dt}\varphi_t(x)\rangle=\langle \nabla{f}(\varphi_t(x)),\nabla{f}(\varphi(t))\rangle \geq 0, \end{align}\] for all \(t \in [0,1]\) and for all \(x\in \mathbb{R}^n\). Hence for every \(x\in \mathbb{R}^n\), the function \(t \mapsto f (\varphi_t(x))\) is non-decreasing, ensuring \(f \circ \varphi_1 - f \circ \varphi_t\geq 0\) for all \(x\in \mathbb{R}^n\). Hence it follows that \(D_{\rm WKL}(\mu \| \nu) \geq 0\). Moreover, \(D_{\rm WKL}(\mu \| \nu)=0\) implies that the statement \[\begin{align} f \circ \varphi_1 = f \circ \varphi_t\textrm{ for all } t \in [0,1] \end{align}\] holds \(\mu\) almost surely. Therefore, we have that for all \(t \in [0,1]\), \[\begin{align} \frac{d}{dt} f(\varphi_t(x))=\langle \nabla{f}(\varphi_t(x)), \frac{d}{dt}\varphi_t(x)\rangle=\langle \frac{d}{dt}\varphi_t(x),\frac{d}{dt}\varphi_t(x)\rangle = 0 \end{align}\] holds \(\mu\) almost surely. Therefore we conclude that \(\frac{d}{dt}\varphi_t(x)=0\) for all \(t\in [0,1]\) and for almost all \(x\). Therefore, \(\varphi_t\) is the identity map \(\mu\) almost everywhere which implies that \(\mu=\nu\). Therefore, \(D_{\rm WKL}\) satisfies the usual properties of a divergence.

Clearly, the outlined definition of the WKL-divergence is rather implicit and requires the knowledge about a potential function \(f\) that induces the transport of \(\mu\) to \(\nu\) in terms of the gradient flow of \(f\). This construction can be made explicit when restricting attention to \(\mathcal{G}\). As a result, we provide an explicit formula for the WKL-divergence in terms of the means \(m_0, m_1\) and covariance matrices \(\Sigma_0,\Sigma_1\) of the respective Gaussian distributions \({\mathcal{N}}(m_0, \Sigma_0)\) and \({\mathcal{N}}(m_1, \Sigma_1)\) belonging to \(\mathcal{G}\). We compare the WKL-divergence with the classical KL-divergence and show that the WKL-divergence is indeed nicely coupled with the geometry of the sample space, that is \(\mathbb{R}^n\).

2.5.0.1 Notation

We use \(I\) to denote the identity matrix of appropriate size. The transpose of a matrix \(A\) is represented by \(A^T\), and its Moore–Penrose pseudoinverse is denoted by \(A^{\dagger}\). The projection matrix onto the null space of \(A\) is denoted by \(A^{\perp} = I - A A^{\dagger}\). The gradient of a function \(f: \mathbb{R}^n\to \mathbb{R}\) evaluated at a point \(x\) is the vector \(\nabla{f}(x)\in \mathbb{R}^n\). We use \(A \succ 0\) (\(A\succeq 0\)) to indicate that \(A\) is symmetric positive definite (semi-definite). The cone of real symmetric positive definite \(n \times n\) matrices is denoted by \(S_n^{+}\), and its closure, the cone of positive semi-definite matrices, is denoted by \(\mathrm{cl}(S_n^{+})\). The matrix exponential is denoted by \(e^A\). For any \(A\succ 0\), the matrix logarithm is denoted by \(\log(A)\) and for any \(A\succeq 0\), the symmetric positive semi-definite square root of \(A\) is denoted by \(\sqrt{A}\) (or \(A^{\frac{1}{2}}\)). Finally, the Frobenius norm of \(A\) is represented by \(\lVert A \rVert_F\) and \(\textrm{Tr}\left({A}\right)\) denotes the trace of \(A\).

3 Main Results↩︎

Let \(f:\mathbb{R}^n\rightarrow \mathbb{R}\) be the quadratic potential function \[\begin{align} \label{eq:potential} f(x)=\frac{1}{2}x^TAx + b^Tx \end{align}\tag{11}\] where \(A\in\mathbb{R}^{n\times n}\) is symmetric¹ and \(b\in \mathbb{R}^{n}\). Consider the gradient flow dynamics given by \[\begin{align} \label{eq:dynamics} \dot{x}(t)=\nabla{f} (x(t)), \quad x(0)=x_0 \end{align}\tag{12}\] and let \(\varphi_t:x_0 \mapsto x(t)\) denote its flow map. Our first Lemma provides a formula for the inner integral on the right hand side of 10 in terms of the \(A\) and \(b\) that define the potential function \(f\).

The gradient flow dynamics 12 can be explicitly solved to obtain \[\begin{align} x(t)= e^{At}x_0+ \left(\int_0^t e^{A(t-\tau)}d\tau\right) b &=e^{At}\underbrace{\left(x_0+A^{\dagger}b\right)}_y+t(I-AA^{\dagger})b-A^{\dagger}b \\ &=e^{At}y+tA^{\perp}b-A^{\dagger}b, \end{align}\] where we have defined \(y=x_0+A^{\dagger}b\) for convenience and have used properties² of \(A^{\perp}\), \(A^{\dagger}\) and \(e^A\).

We now compute \(f\circ \varphi_t\) \[\begin{align} f(x(t))&=\frac{1}{2}x(t)^TAx(t) + b^Tx(t)=\frac{1}{2} \left(y^T Ae^{2At}y-b^TA^{\dagger}b\right)+tb^TA^{\perp}b+y^TA^{\perp}b \end{align}\] Therefore, \[\begin{align} (f\circ\varphi_1) (x_0) - (f \circ \varphi_t) (x_0) =\frac{1}{2} y^T A\left(e^{2A}-e^{2At}\right)y + (1-t)b^TA^{\perp}b. \end{align}\] Integrating with respect to time, we get, \[\begin{align} \int_0^1 \left(f \circ \varphi_1 - f \circ \varphi_t\right)(x_0) dt&=\int_0^1\left(\frac{1}{2} y^T A\left(e^{2A}-e^{2At}\right)y + (1-t) b^TA^{\perp}b \right) dt \\ &=\frac{1}{4} y^T \left(2Ae^{2A}-e^{2A}+I\right)y+ \frac{1}{2} b^TA^{\perp}b. \end{align}\] Plugging in the definitions of \(y\) and \(M\), we get the desired identity. The next lemma considers two Gaussian distributions \(\mu\) and \(\nu\) and provides a function \(f\) with the desired property that \(\nu\) is the image of \(\mu\) with respect to its gradient flow after a unit of time.

We have seen in the proof of Lemma 1, that dynamics 12 can be solved to obtain \[\begin{align} x(t)=e^{At}x_0+\left(e^{At}A^{\dagger}+tA^{\perp}-A^{\dagger}\right)b. \end{align}\] Using properties of Gaussian random variables, it can be shown that \[\begin{align} x_0 \sim \mathcal{N}(m_0,\Sigma_0) \implies x(t) \sim \mathcal{N}\left(\underbrace{e^{At}{m}_0+\left(e^{At}A^{\dagger}+tA^{\perp}-A^{\dagger}\right)b}_{m_t},\underbrace{e^{At}\Sigma_0 e^{At}}_{\Sigma_t}\right). \end{align}\] Observe that with the prescribed \(A\) and \(b\) given in ?? and ?? , we indeed get \[\begin{align} \label{eq:riccatti} \Sigma_1&=e^A\Sigma_0e^A,\\ m_1&=e^{A}{m}_0+\left((e^{A}-I)A^{\dagger}+A^{\perp}\right)b. \end{align}\tag{13}\] Note that with a change of variable \(X=e^A\), 13 reduces to the matrix equation \(X\Sigma_0 X = \Sigma_1\) which is a special case of the algebraic Riccati equation and has been extensively studied in control theory (see [8] for example). Since \(e^A \succ 0\) for any symmetric \(A\), we are interested in the positive definite solutions \(X\) to the matrix equation \(X\Sigma_0 X = \Sigma_1\). Under the constraints that \(\Sigma_0\succ 0\) and \(\Sigma_1\succ 0\), uniqueness of the positive definite solution follows from [8] proving the uniqueness of the solution \(A\) to 13 . Uniqueness of \(b\) is obtained immediately since \(\left((e^{A}-I)A^{\dagger}+A^{\perp}\right)\) is non-singular. This proves the final statement. We now present the main result of the paper which provides a formula for the WKL-divergence between two Gaussian distributions.

Note that the outer integral in the right-hand-side of 10 corresponds to taking an expectation with respect to \(\mu\). Using Lemma 2, we define \(A\) and \(b\) according to ?? and ?? to obtain the property that \(x(0) \sim \mu\) implies \(x(1) \sim \nu\). Using Lemma 1 along with the well-known properties of the expectation and trace operators³, we get that \[\begin{align} D_{\rm WKL}(\mu \| \nu)&=\int_{\mathbb{R}^n}\int_0^1 \left(f \circ \varphi_1 - f \circ \varphi_t\right)(x_0) dt d\mu(x_0)\nonumber \\ &=\mathbb{E}\left[\frac{1}{4}(x_0+A^{\dagger}b)^TM(x_0+A^{\dagger}b) + \frac{1}{2} b^TA^{\perp}b\right]\nonumber \\ &=\frac{1}{4}\left(\textrm{Tr}\left({(M\Sigma_0)}\right)+(m_0+A^{\dagger}b)^TM(m_0+A^{\dagger}b)\right) + \frac{1}{2} b^TA^{\perp}b \label{eq:D95of95mu95nu} \end{align}\tag{14}\] where \(A\) and \(b\) are given by ?? and ?? . We now substitute the expression ?? for \(b\) into 14 term by term. Let us first compute and simplify \(m_0 + A^{\dagger}b\) as \[\begin{align} m_0 + A^{\dagger} b & = m_0 + A^{\dagger}\left( A^{\dagger}(e^{A}-I)+A^{\perp} \right)^{-1}(m_1-e^Am_0) \nonumber \\ & = m_0 + (e^{A}-I)^{\dagger}(m_1-e^Am_0) \nonumber = (e^{A}-I)^{\perp}m_0+ (e^{A}-I)^{\dagger}(m_1-m_0). \nonumber \end{align}\] Using the fact that \(M(e^{A}-I)^{\perp}=0\), we get that \[\begin{align} (m_0+A^{\dagger}b)^TM(m_0+A^{\dagger}b) &=(m_1-m_0)^T (e^{A}-I)^{\dagger}M (e^{A}-I)^{\dagger}(m_1-m_0).\label{eq:term2} \end{align}\tag{15}\] Since \(A^{\perp}\left( A^{\dagger}(e^{A}-I)+A^{\perp} \right)^{-1}=A^{\perp}\) and \(A^{\perp}e^{A}=A^{\perp}\), the last term can be computed as \[\begin{align} b^TA^{\perp}b=(m_1-e^Am_0)^TA^{\perp}(m_1-e^Am_0) &=(m_1-m_0)^TA^{\perp}(m_1-m_0).\label{eq:term3} \end{align}\tag{16}\] Finally, the cyclic property of the trace operator gives us \[\begin{align} \textrm{Tr}\left({M\Sigma_0}\right) &=\textrm{Tr}\left({\left(\log(R^2)R^2-R^2+I\right)\Sigma_0}\right) \nonumber\\ &=\textrm{Tr}\left({\log(R^2)R^2\Sigma_0-\Sigma_1+\Sigma_0}\right) \label{eq:term1}. \end{align}\tag{17}\] Plugging in 15 , 16 and 17 in 14 and using \(M=2Ae^{2A}-e^{2A}+I=\log(R^2)R^2-R^2+I\) as well as \(A=\log(R)\) gives the desired result. We now show that \[\begin{align} \label{eq:mat95in95mean95norm} Q+2\log(R)^{\perp}=(e^{A}-I)^{\dagger}M (e^{A}-I)^{\dagger}+2A^{\perp} \succ 0 \end{align}\tag{18}\] ensuring the square root is well defined. Consider the orthogonal eigenvalue decomposition \(A=U\Lambda U^T\) where \(\Lambda\) is a diagonal matrix consisting of eigenvalues \(\lambda_i\), \(i \in \{1,2,\cdots,n\}\) of \(A\) and \(U\) is an orthogonal matrix consisting of eigenvectors as columns. Since \(A\) and \(e^{A}\) have the same eigenvectors, observe that \(U^TMU=2\Lambda e^{2\Lambda}-e^{2\Lambda}+I\), which is a diagonal matrix consisting of entries \(2\lambda_i e^{2\lambda_i}-e^{2\lambda_i}+1\), \(i \in \{1,2,\cdots,n\}\) on the diagonal. Observe that the function \(h:\lambda\mapsto (e^{-2\lambda}-1)\) is strictly convex, which implies that \(e^{-2\lambda}-1=h(\lambda)> h(0)+h'(0)\lambda=-2\lambda\) for all \(\lambda \in \mathbb{R}\setminus \{0\}\) which implies that \(2\lambda e^{2\lambda}-e^{2\lambda}+1> 0\) for all \(\lambda \in \mathbb{R}\setminus \{0\}\). Hence each diagonal entry of \(U^T M U\) is nonnegative, with equality if and only if \(\lambda_i = 0\). That directly implies \(M \succeq 0\) and \(M\) has an eigenvalue at \(0\) if and only if \(A\) has an eigenvalue at \(0\). Therefore, \(Q=(e^{A}-I)^{\dagger}M (e^{A}-I)^{\dagger}\succeq 0\). Furthermore, note that \(z^TQz=0\) only if \(z\) belongs to the null space of \(A\) which implies that \(z^TA^{\perp}z=\lVert z\rVert^2\). This shows that \(z^T\left(Q+2\log(R)^{\perp}\right)z=0\) if and only if \(z=0\) thereby showing that \(Q+2\log(R)^{\perp}\succ 0\). This ensures the square root is well defined.

Finally, when \(\Sigma_0 \Sigma_1 = \Sigma_1 \Sigma_0\), the first term can be reformulated using \(M(e^{A}-I)^{\perp}=0\) as \[\begin{align} \textrm{Tr}\left({M\Sigma_0}\right) &= \textrm{Tr}\left({(e^{A}-I)(e^{A}-I)^{\dagger}M(e^{A}-I)^{\dagger}(e^{A}-I)\Sigma_0}\right) \\ &=\left\lVert \sqrt{Q}\left(\sqrt{\Sigma_1}-\sqrt{\Sigma_0}\right) \right\rVert_F^2. \end{align}\]

Let us now show that the WKL-divergence ?? is the sum of two non-negative terms. Furthermore, the first term is zero if and only if \(\Sigma_1=\Sigma_0\). This can be seen by observing that the first term \(\frac{1}{4}\textrm{Tr}\left({\Sigma_0-\Sigma_1+\Sigma_0R^2\log(R^2)}\right)=\frac{1}{4}\textrm{Tr}\left({M\Sigma_0}\right)\) (see 17 ) where \(M=2Ae^{2A}-e^{2A}+I\succeq 0\). Therefore, \(\textrm{Tr}\left({M \Sigma_0}\right)=\textrm{Tr}\left({\sqrt{\Sigma_0}M \sqrt{\Sigma_0}}\right)=0\) implies \(\sqrt{\Sigma_0}M \sqrt{\Sigma_0}=0\) and thus, \(M=0\). Recall from the proof of Theorem 5 that \(M\) and has eigenvalues \(2\lambda_i e^{2\lambda_i}-e^{2\lambda_i}+1\) where \(\lambda_i\), \(i\in\{1,2,\cdots,n\}\) are the eigenvalues of \(A\). Thus \(M=0\) implies that \(2\lambda_i e^{2\lambda_i}-e^{2\lambda_i}+1=0\) which in turn implies that \(\lambda_i=0\) for all \(i\in\{1,2,\cdots,n\}\) (see proof of Theorem 5). Therefore, \[\begin{align} A=\log \left(\Sigma_0^{-\frac{1}{2}}\left(\Sigma_0^{\frac{1}{2}}\Sigma_1\Sigma_0^{\frac{1}{2}}\right)^{\frac{1}{2}}\Sigma_0^{-\frac{1}{2}}\right)=0 \end{align}\] which gives us \(\Sigma_1=\Sigma_0\). Similarly, owing to the fact that \(Q+2\log(R)^{\perp}\succ 0\) (see proof of Theorem 5), we see that the second term \(\frac{1}{4}\left\lVert \sqrt{Q+2\log(R)^{\perp}} (m_1-m_0) \right\rVert^2\) is zero if and only if \(m_1=m_0\). Finally, note that when \(\Sigma_1=\Sigma_0\), we get that \(R=I\) and \(Q=0\). This gives us that \(D_{\rm WKL}(\mu \| \nu) = \frac{1}{2}\left\lVert (m_1-m_0) \right\rVert^2\) showing that the WKL-divergence indeed captures the geometry of the sample space. Finally, let us compare the WKL-divergence to the KL-divergence between Gaussian distributions given by \[\begin{align} \label{eq:D95main95formula95KL} D_{\rm KL}(\nu \| \mu) &= \frac{1}{2}\left(\log \left(\frac{\det \Sigma_0}{\det \Sigma_1}\right)-n + \textrm{Tr}\left({\Sigma_0^{-1}\Sigma_1}\right)\right) +\frac{1}{2}\left\lVert \Sigma_0^{-\frac{1}{2}} (m_1-m_0) \right\rVert^2, \end{align}\tag{19}\] which is clearly distinct from \(D_{\rm WKL}(\mu \| \nu)\).

3.1 Gaussian transport examples↩︎

To illustrate the geometry underlying the Wasserstein KL-divergence, we consider three explicit Gaussian transport examples in dimension two. In each case, the transport is generated by an affine vector field, so the associated \(e_1\)-geodesic is represented by a Gaussian flow whose covariance and mean evolve in closed form. This is depicted in Figure 1.

3.2 Relation to Optimal Transport or Wasserstein Distance↩︎

The 2-Wasserstein distance \(W_2\) is a Riemannian distance on the space of probability distributions (more precisely, on the density manifold [9]), arising as the geodesic distance induced from the Otto metric [6]. In this sense, it plays a role analogous to that of the Fisher\(-\)Rao distance in information geometry. While our proposed WKL-divergence is not a distance, it is closely related to Wasserstein geometry, and therefore invites comparison with \(W_2\), especially in the Gaussian case where both quantities admit closed-form expressions that have a close resemblance.

We review the case of zero-mean Gaussian distributions \(\mu = \mathcal{N}(0, \Sigma_0)\) and \(\nu = \mathcal{N}(0, \Sigma_1)\) and point the reader to [10] for further details. The 2-Wasserstein distance between them arises as the minimal cost in an optimal transport problem. Specifically, the task is to find a linear transport map represented by a matrix \(R \in \mathbb{R}^{n \times n}\) such that if \(x \sim \mu\), then \(Rx \sim \nu\). This requirement leads to the algebraic condition \(\Sigma_1 = R \Sigma_0 R^T\), which generally admits multiple solutions. To ensure uniqueness, we look at the linear maps \(R\) that minimize the cost \(\mathbb{E}_{x \sim \mu}\left[\lVert x-Rx\rVert ^2\right]\) subject to the constraint \(\Sigma_1 = R \Sigma_0 R^T\). This leads to the map \[\begin{align} \label{eq:brenier95map} R=\Sigma_0^{-\frac{1}{2}}\left(\Sigma_0^{\frac{1}{2}}\Sigma_1\Sigma_0^{\frac{1}{2}}\right)^{\frac{1}{2}}\Sigma_0^{-\frac{1}{2}}, \end{align}\tag{20}\] which is the unique optimal transport map in this setting. The square of the associated 2-Wasserstein distance is then given by the minimal cost under this map. An alternative and equivalent strategy to ensure uniqueness of solutions to \(\Sigma_1 = R \Sigma_0 R^T\) is to restrict attention to linear maps \(R\) that can be described as the gradient of a convex (quadratic) potential which corresponds to requiring \(R\succeq 0\) yielding a unique solution given by 20 . For general Gaussian distributions (not necessarily zero mean) \(\mu = \mathcal{N}(m_0, \Sigma_0)\) and \(\nu = \mathcal{N}(m_1, \Sigma_1)\), the squared 2-Wasserstein distance is given by \[\begin{align} W_2(\mu,\nu)^2=\left \lVert m_0-m_1\right \rVert^2 + \textrm{Tr}\left({\Sigma_0 + \Sigma_1-2\left(\Sigma_0^{\frac{1}{2}}\Sigma_1\Sigma_0^{\frac{1}{2}}\right)^{\frac{1}{2}}}\right) \end{align}\] where the second term reduces to \(\left \lVert \sqrt{\Sigma_0} - \sqrt{\Sigma_1}\right \rVert_F^2\) when \(\Sigma_0\) and \(\Sigma_1\) commute.

Comparing this with the expression for \(D_{\rm WKL}\) given in Theorem 5, we observe a close structural resemblance: both quantities separate into a quadratic mean difference term and a term involving square roots of the covariance matrices. However, the WKL-divergence involves an additional factor in both terms. A more direct connection emerges when comparing the transport maps involved in both constructions. Comparing the map in 20 with the expression ?? , we see that \(e^A\) plays the role of \(R\). Furthermore, as shown in Lemma 2, the expression for \(A\) is chosen so that \(e^A\) transports \(\mu\) to \(\nu\) via a gradient flow, and symmetry of \(A\) ensures that this flow arises from a quadratic (not necessarily convex) potential. The matrix exponential \(e^A\) is positive definite irrespective of \(A\) and hence is uniquely determined without assuming convexity of the potential. In contrast, the uniqueness of the map \(R\) in optimal transport is ensured by restricting to the class of convex quadratic potentials.

3.3 Univariate Gaussian Distributions↩︎

We now take a closer look at the simple case of univariate Gaussian distributions, i.e., \(n=1\). The main formula is given in the following corollary which is a direct application of Theorem 5.

Substituting \(A\leftarrow a\), \(\Sigma_0\leftarrow \sigma_0^2\), \(\Sigma_1\leftarrow \sigma_1^2\), we get \[\begin{align} R=\frac{\sigma_1}{\sigma_0} \quad \textrm{ and } \quad Q=\left(\sigma_0^2-\sigma_1^2+\sigma_1^2\log\left(\frac{\sigma_1^2}{\sigma_0^2}\right) \right)\left(\left(\sigma_1-\sigma_0\right)^2\right)^{\dagger}. \end{align}\] This leads to \[\begin{align} &D_{\rm WKL}\left(\mu || \nu\right) = \frac{Q}{4} \left((\sigma_1-\sigma_0)^2 + (m_1-m_0)^2\right) + \frac{1}{2}\log(R)^{\perp} (m_1-m_0)^2 \end{align}\] which gives the desired result. The final limit can be shown via successive applications of L’Hopital’s rule.

In contrast to the WKL-divergence, the KL-divergence given in 19 for the univariate Gaussian distributions simplifies to \[\begin{align} D_{\rm KL}\left(\mu || \nu\right) = \log \left(\frac{\sigma_0}{\sigma_1}\right)+\frac{\sigma_1^2}{2\sigma_0^2}+\frac{(m_1-m_0)^2}{2\sigma_0^2} - \frac{1}{2}. \end{align}\] Therefore, if we consider two univariate Gaussian distributions \(\mu=\mathcal{N}(m_0,\sigma^2)\) and \(\nu=\mathcal{N}(m_1,\sigma^2)\) with equal variance \(\sigma^2\) and possibly different means \(m_0\) and \(m_1\), the WKL-divergence and the classical KL-divergence between these distributions can be computed as \[\begin{align} D_{\rm WKL}(\mu\lVert\nu)= \frac{(m_1-m_0)^2}{2} \quad\quad \textrm{and}\quad\quad D_{\rm KL}(\nu\lVert\mu)= \frac{(m_1-m_0)^2}{2\sigma^2}. \end{align}\] Note that since our WKL-divergence corresponds to the dual of the KL-divergence [5], a direct comparison requires reversing the order of the distributions. It is clear from this example that as \(\sigma \rightarrow 0\), the KL-divergence diverges to infinity whereas the WKL-divergence remains constant since it depends solely on the distance between the means. This is illustrated in Figure 2 (left) where we set \(m_0=0\), \(m_1=2\) and the standard deviations \(\sigma_0=\sigma_1=\sigma\) is varied. This illustrates the coupling of the divergence with the geometry of the sample space which in this case is \(\mathbb{R}\). This allows us to approximate the divergence between two Dirac measures such that it is proportional to the distance in the sample space.

We now look at the local curvature of the KL-divergence and the WKL-divergence around an optimum which has a strong influence on performance of optimization algorithms. The curvature corresponds to the Hessian evaluated at the optimum and it can be shown that \[\begin{align} \frac{\partial^2}{\partial \sigma^2} D_{\rm WKL}(\mu\lVert\nu)\bigg|_{\sigma=\sigma_{\textrm{opt}}}= 1 \quad\quad \textrm{and}\quad\quad \frac{\partial^2}{\partial \sigma^2} D_{\rm KL}(\nu\lVert\mu)\bigg|_{\sigma=\sigma_{\textrm{opt}}}= \frac{2}{\sigma^2_{\textrm{opt}}}. \end{align}\] Observe that local curvature of the WKL-divergence is independent of \(\sigma_{\textrm{opt}}\) whereas the local curvature of the KL-divergence blows up to infinity as \(\sigma_{\textrm{opt}}\) approaches zero. Figure 2 (right) plots the divergences for \(m=0\) and \(\sigma_{\textrm{opt}}\in\{1,3\}\) and it can be observed that the local curvature of \(D_{\rm WKL}\) is independent of \(\sigma_{\textrm{opt}}\) whereas the local curvature of \(D_{\rm KL}\) is high for \(\sigma_{\textrm{opt}}=1\) than that for \(\sigma_{\textrm{opt}}=3\).

Figure 3 displays the sub-level sets (divergence balls) for both divergences, each centered at different reference distributions. We observe that the sub-level sets of the KL-divergence are sensitive to changes in \(\sigma_0\), as evidenced by the significant variation in their size when \(\sigma_0\) varies. In contrast, the sizes of the WKL-divergence sub-level sets exhibit much less sensitivity to \(\sigma_0\). Finally, while the KL-divergence consistently produces convex sub-level sets regardless of the reference parameters, the WKL-divergence sub-level sets become non-convex when \(\sigma_0\) is small.

Figure 4 depicts gradient descent trajectories in the \((m_1,\sigma_1)\)-plane, initialized from a uniform grid around the reference parameters \((m_0,\sigma_0)\). The plot illustrates the local curvature structure of \(D_{\rm KL}\) and \(D_{\rm WKL}\) and highlights the basins of attraction induced by their gradient flows. For \(D_{\rm KL}\), all trajectories converge reliably to the global minimizer at \((m_0,\sigma_0)\). In contrast, for \(D_{\rm WKL}\) most trajectories converge to the global minimizer, but certain initializations are drawn toward the boundary \(\sigma_1=0\). This behavior indicates that, while the WKL-divergence inherits desirable geometric properties, its optimization landscape may admit boundary-attracting trajectories. Understanding this phenomenon and developing robust optimization methods in the WKL framework constitute interesting directions for future research.

3.4 Continuity Properties of \(D_{\rm WKL}\)↩︎

In this subsection we collect the main continuity and boundary behaviour of the WKL-divergence \(D_{\mathrm{WKL}}\) on the Gaussian family. In particular, we (i) derive a uniform lower bound that links \(D_{\mathrm{WKL}}\) to the Euclidean distance between the means, (ii) describe how the divergence behaves when one or both covariance matrices approach singularity, and (iii) contrast these phenomena with the classical KL-divergence.

We start by deriving a uniform lower bound on \(D_{\rm WKL}\) in the next proposition.

Recall that we have already seen in the discussion after Theorem 5 that the WKL-divergence is the sum of two non-negative terms. The second term thus already provides a lower bound on the WKL-divergence, that is, \[\begin{align} D_{\rm WKL}(\mu \| \nu) \geq \frac{1}{4}\left\lVert \sqrt{Q+2\log(R)^{\perp}} (m_1-m_0) \right\rVert^2. \end{align}\] Note that if \(Q+2\log(R)^{\perp}-I\succeq 0\), then we have that \[\begin{align} \frac{1}{4}\left\lVert \sqrt{Q+2\log(R)^{\perp}} (m_1-m_0) \right\rVert^2&= \frac{1}{4}(m_1-m_0)^T\left(Q+2\log(R)^{\perp}\right)(m_1-m_0)\\ &\geq \frac{1}{4}(m_1-m_0)^T(m_1-m_0) =\frac{1}{4}\left\lVert m_1-m_0 \right\rVert^2 \end{align}\] which gives us the first statement. In order to see that \(Q+2\log(R)^{\perp}-I\succeq 0\) is true, consider the orthogonal eigenvalue decomposition \(R=P\Lambda P^T\) and multiply the matrix \(Q+2\log(R)^{\perp}\) from the left and right by \(P^T\) and \(P\), respectively, to diagonalize it. This gives us \[\begin{align} P^T \left(Q+2\log(R)^{\perp}\right) P&=(\Lambda-I)^{\dagger}\left(\log(\Lambda^2)\Lambda^2-\Lambda^2+I\right)(\Lambda-I)^{\dagger}+2\log(\Lambda)^{\perp} \end{align}\] which is a diagonal matrix with the \(i^{\textrm{th}}\) diagonal entry given by \[\begin{align} g(\lambda_i):=\begin{cases} \frac{\log(\lambda_i^2)\lambda_i^2 - \lambda_i^2+1}{(\lambda_i-1)^2}\textrm{ if } \lambda_i \neq 1,\\2\textrm{otherwise}, \end{cases} \end{align}\] where \(\lambda_i\) is the \(i^{\textrm{th}}\) diagonal entry of \(\Lambda\). We next show that \(g(\lambda_i)\geq 1\) which implies that the eigenvalues of \(\left(Q+2\log(R)^{\perp}\right)\) are greater than or equal to \(1\) proving \(Q+2\log(R)^{\perp}-I\succeq 0\). Observe that convexity of the function \(h:t\mapsto t\log t\) implies that \[\begin{align} \lambda_i\log \lambda_i =h(\lambda_i)\geq h(1)+h'(1)(\lambda_i-1)=\lambda_i-1. \end{align}\] Multiplying both sides by \(2\lambda_i\) (note that \(\lambda_i> 0\) since \(R\succ 0\)), then adding \(1\) and subtracting \(\lambda_i^2\) from both sides gives \[\begin{align} \lambda_i^2 \log (\lambda_i^2)-\lambda_i^2+1=2\lambda_i^2\log \lambda_i+1-\lambda_i^2 \geq \lambda_i^2-2\lambda_i+1=(\lambda_i-1)^2. \end{align}\] This completes the proof of the first statement.
Let \(\Sigma_0\succ 0\) be fixed and consider a sequence of symmetric positive definite matrices \(\Sigma_1^{(1)}, \Sigma_1^{(2)}, \cdots\) that converges to the zero matrix. Hence \[\begin{align} R^{(i)}&=\Sigma_0^{-\frac{1}{2}}\left(\Sigma_0^{\frac{1}{2}}\Sigma_1^{(i)}\Sigma_0^{\frac{1}{2}}\right)^{\frac{1}{2}}\Sigma_0^{-\frac{1}{2}} \rightarrow 0,\\ Q^{(i)}&=\left(R^{(i)}-I\right)^{\dagger}\left(\log\left(\left(R^{(i)}\right)^2\right)\left(R^{(i)}\right)^2-\left(R^{(i)}\right)^2+I\right)\left(R^{(i)}-I\right)^{\dagger} \rightarrow I. \end{align}\] Furthermore, note that there exists a natural number \(N\), such that for \(i > N\), \(\lVert R^{(i)}\rVert <1\) which implies that the largest eigenvalue of \(R^{(i)}\) is less than \(1\). Furthermore, since \(R^{(i)}\succ 0\), we get that \(\log\left(R^{(i)}\right)^{\perp}=0\) for all \(i>N\). This gives us that \[\begin{align} \lim_{\Sigma_1 \rightarrow 0} D_{\rm WKL}\left(\mathcal{N}(m_0,\Sigma_0) \| \mathcal{N}(m_1,\Sigma_1)\right) &= \frac{1}{4}\textrm{Tr}\left({\Sigma_0}\right)+ \frac{1}{4}\left\lVert m_1-m_0 \right\rVert^2. \end{align}\] Taking the limit as \(\Sigma_0\) converges to the zero matrix gives us the desired result ?? . In contrast to the limit obtained in ?? , we get a different limit if we approach the origin along the diagonal, i.e., when \(\Sigma_0=\Sigma_1\). To see this, note that if \(\Sigma_0=\Sigma_1=\Sigma\succ 0\), then we have that \(R=I\), \(Q=0\) and \(\log (R)^{\perp}=I\). Plugging this in ?? gives us \[\begin{align} D_{\rm WKL}\left(\mathcal{N}(m_0,\Sigma) \| \mathcal{N}(m_1,\Sigma)\right) &= \frac{1}{2}\left\lVert m_1-m_0 \right\rVert^2\label{eq:limit2} \end{align}\tag{21}\] which is independent of \(\Sigma\). Therefore, \[\begin{align} \lim_{\Sigma \rightarrow 0} D_{\rm WKL}\left(\mathcal{N}(m_0,\Sigma) \| \mathcal{N}(m_1,\Sigma)\right) \neq \lim_{\Sigma_0 \rightarrow 0}\left(\lim_{\Sigma_1 \rightarrow 0} D_{\rm WKL}\left(\mathcal{N}(m_0,\Sigma_0) \| \mathcal{N}(m_1,\Sigma_1)\right)\right). \end{align}\] In contrast to the finite limits described above, the limit as \(\Sigma_0\) approaches singularity (when \(\Sigma_1\succ 0\) is fixed) is not finite. To see this, let \(\Sigma_0=\sigma_0^2 I\) which gives us \(R=\frac{1}{\sigma_0}\Sigma_1\) which grows unbounded as \(\sigma_0\) approaches \(0\). Hence the WKL-divergence grows unbounded as \(\sigma_0\) approaches \(0\). We summarize these observations as follows:

1. If \(\Sigma_1\to 0\) and then \(\Sigma_0\to 0\), the WKL-divergence converges to \(\frac{1}{4}\|m_1-m_0\|^2\).

2. Along \(\Sigma_0=\Sigma_1=\Sigma\), the WKL-divergence equals \(\frac{1}{2}\|m_1-m_0\|^2\) and is independent of \(\Sigma\)

3. If \(\Sigma_1\succ 0\) is fixed and \(\Sigma_0\to0\), then the WKL-divergence diverges to \(\infty\).

The first two items show that the WKL-divergences can be naturally extended to give a finite squared-distance between Dirac measures. The following theorem formally analyzes the continuity properties of the WKL-divergence.

We divide the argument into several steps.
Step 1. Continuity on the interior. Recall that the maps \[\begin{align} S_n^{+} \times S_n^{+} \ni (X,Y) &\mapsto XYX \in S_n^{+}, \\ S_n^{+} \ni X &\mapsto X^{\frac{1}{2}} \in S_n^{+},\\ S_n^{+} \ni X &\mapsto X^{-1} \in S_n^{+} \end{align}\] are continuous and compositions of continuous maps are continuous. Therefore, the map \[\begin{align} (\Sigma_0,\Sigma_1) \mapsto R=\Sigma_0^{-\frac{1}{2}}\left(\Sigma_0^{\frac{1}{2}}\Sigma_1\Sigma_0^{\frac{1}{2}}\right)^{\frac{1}{2}}\Sigma_0^{-\frac{1}{2}} \end{align}\] is continuous on \(S_n^{+} \times S_n^{+}\). Consequently, the map \[\begin{align} (\Sigma_0,\Sigma_1) \mapsto \frac{1}{4}\textrm{Tr}\left({\Sigma_0-\Sigma_1+\Sigma_0R^2\log(R^2)}\right) \end{align}\] is continuous which proves the continuity of the first term in ?? .
We now show that the map \[\begin{align} R \mapsto (R-I)^{\dagger}\left(\log(R^2)R^2-R^2+I\right)(R-I)^{\dagger} + 2\log(R)^{\perp} \label{eq:R95map} \end{align}\tag{22}\] is continuous on \(S_n^{+}\). To this end, consider a sequence of matrices \(R_i\), \(i \in \{1,2,\cdots\}\) converging to \(R\) with a spectral decomposition \(R = P \Lambda P^T\). Standard perturbation theory for symmetric matrices allows us to choose spectral decompositions \(R_i = P_i \Lambda_i P_i^T\) such that \(P_i \to P\) and \(\Lambda_i \to \Lambda\). Then \[\begin{align} &(R_i-I)^{\dagger}\left(\log(R_i^2)R_i^2-R_i^2+I\right)(R_i-I)^{\dagger}+2\log(R)^{\perp} \nonumber \\ &=P_i\left((\Lambda_i-I)^{\dagger}\left(\log(\Lambda_i^2)\Lambda_i^2-\Lambda_i^2+I\right)(\Lambda_i-I)^{\dagger}+2\log(\Lambda)^{\perp}\right)P_i^T, \label{eq:sequence95EVD} \end{align}\tag{23}\] where the central matrix is a diagonal matrix with entries of the form \[\begin{align} g(\lambda)= (\lambda-1)^{\dagger}(\lambda^2 \log (\lambda^2) - \lambda^2+1)(\lambda-1)^{\dagger} + 2\log(\lambda)^{\perp}, \label{eq:diag95scalar95term} \end{align}\tag{24}\] where \(\lambda>0\) is a place-holder for the eigenvalues of the matrix \(\Lambda_i\). Here \((\lambda-1)^{\dagger}\) denotes the scalar pseudoinverse, i.e. \[(\lambda-1)^{\dagger}=\begin{cases} \frac{1}{\lambda-1}, & \lambda\neq 1,\\ 0, & \lambda=1. \end{cases}\] Thus \(g(\lambda)\) equals \[g(\lambda)= \begin{cases} \frac{\lambda^2\log(\lambda^2)-\lambda^2+1}{(\lambda-1)^2}, & \lambda\neq 1,\\[6pt] 2, & \lambda=1, \end{cases}\] which is continuous for all \(\lambda>0\). A direct calculation shows \[\begin{align} \lim_{\lambda \to 1} \frac{\lambda^2 \log (\lambda^2) - \lambda^2+1}{(\lambda-1)^2} = 2 \end{align}\] which implies that the \(g\) extends continuously to all \(\lambda \in (0,\infty)\) and therefore, the central matrix in 23 converges to \[\begin{align} (\Lambda-I)^{\dagger}\left(\log(\Lambda^2)\Lambda^2-\Lambda^2+I\right)(\Lambda-I)^{\dagger}+2\log(\Lambda)^{\perp} \end{align}\] whenever \(\Lambda_i\) converges to \(\Lambda\). Furthermore, the left and right factors in 23 converge to \(P\) and \(P^T\), respectively. Hence, the map 22 is continuous on \(S_n^{+}\). Finally, since the map \[\begin{align} (\mathbb{R}^n\times S_n^{+}) \ni (m,X) \mapsto m^T X m \in [0,\infty) \end{align}\] is continuous, it follows that the second term of ?? is continuous. Thus \(\mathcal{L}\) is continuous on \((\mathbb{R}^n\times S_n^{+}) \times (\mathbb{R}^n\times S_n^{+})\).
Step 2. Extension when \(\Sigma_1\) approaches the boundary. Fix \(\Sigma_0 \succ 0\) and let \(\Sigma_1\) approach the boundary of \(S_n^{+}\), i.e., let one of the eigenvalues of \(\Sigma_1\) approach \(0\). It follows that since \(\Sigma_0\succ 0\), one eigenvalue of \(R\) also converges to \(0\). Nevertheless, since \[\begin{align} \lim_{\lambda \to 0} \lambda \log(\lambda) = 0, \end{align}\] the map \(R\mapsto R^2\log(R^2)\) can be continuously extended to \(S_n^+\). The same reasoning along with a diagonalization step shows that the map \[\begin{align} R \mapsto (R-I)^{\dagger}\left(\log(R^2)R^2-R^2+I\right)(R-I)^{\dagger} + 2\log(R)^{\perp} \end{align}\] can also be continuously extended to \(\mathrm{cl}(S_n^{+})\). Therefore, \(\mathcal{L}\) can be continuously extended to \(\mathbb{R}^n\times S_n^{+} \times \mathbb{R}^n\times \mathrm{cl}(S_n^{+})\).
Step 3. Extension when \(\Sigma_0\) approaches the boundary. Let us now consider the case when \(\Sigma_1\succ 0\) is fixed and \(\Sigma_0\) approaches the boundary of \(S_n^{+}\), i.e., at least one of the eigenvalues of \(\Sigma_0\) approaches \(0\). Let \(\lambda_{\textrm{min}}\left(\Sigma_0\right)\) be the smallest eigenvalue of \(\Sigma_0\) with the corresponding eigenvector \(v\), i.e., \(\Sigma_0 v= \lambda_{\textrm{min}}(\Sigma_0) \cdot v\). Then, we get that \[\begin{align} v^T R v =\frac{1}{\lambda_{\textrm{min}}(\Sigma_0)} v^T\left(\Sigma_0^{\frac{1}{2}}\Sigma_1\Sigma_0^{\frac{1}{2}}\right)^{\frac{1}{2}}v &\geq \frac{1}{\lambda_{\textrm{min}}(\Sigma_0)}\sqrt{\lambda_{\textrm{min}}\left(\Sigma_0\right)\lambda_{\textrm{min}}\left(\Sigma_1\right)}\lVert v \rVert^2\\ &=\sqrt{\frac{\lambda_{\textrm{min}}\left(\Sigma_1\right)}{\lambda_{\textrm{min}}\left(\Sigma_0\right)}}\lVert v \rVert^2. \end{align}\] Thus, as the smallest eigenvalue of \(\Sigma_0\) approaches \(0\), the largest eigenvalue of \(R\) diverges to \(\infty\). Therefore we may extend \(\mathcal{L}\) continuously by setting its value to \(\infty\) whenever \(\Sigma_0\) is singular and \(\Sigma_1 \succ 0\).
Step 4. Boundary case when both \(\Sigma_0,\Sigma_1\) are singular. Finally, on the boundary when both \(\Sigma_0\) and \(\Sigma_1\) are singular, \(\mathcal{L}\) cannot be continuously extended. For instance, if \(\Sigma_0 = \Sigma_1 = 0\), then \[\begin{align} &\lim_{s\to 0} \lim_{t\to 0} \left( (\Sigma_0+sI)^{-\frac{1}{2}} \left((\Sigma_0+sI)^{\frac{1}{2}}(\Sigma_1+tI)(\Sigma_0+sI)^{1/2}\right)^{\frac{1}{2}} (\Sigma_0+sI)^{-1/2} \right)\\ &= \lim_{s\to 0} \lim_{t\to 0}\sqrt{\frac{t}{s}}=0, \end{align}\] whereas reversing the order of limits gives \(\infty\). Thus continuity may fail at such points. Nevertheless, using the uniform lower bound ?? , we can define a lower semi-continuous extension given in ?? . By construction, this extension is lower semi-continuous. This completes the proof.

Let us now specialize Theorem 7 to the case of univariate Gaussian distributions. Let \(\mathcal{L}: \mathbb{R} \times (0,\infty) \times \mathbb{R} \times (0,\infty) \to [0,\infty]\) be defined by \[\begin{align} \mathcal{L}\big(m_0,\sigma_0^2,m_1,\sigma_1^2\big) = D_{\rm WKL}\big(\mathcal{N}(m_0,\sigma_0^2) \,\|\, \mathcal{N}(m_1,\sigma_1^2)\big), \end{align}\] where \(D_{\rm WKL}\) is given in ?? . Applying Theorem 7, we can continuously extend \(\mathcal{L}\) to \(\big(\mathbb{R} \times (0,\infty) \times \mathbb{R} \times [0,\infty)\big) \;\cup\; \big(\mathbb{R} \times [0,\infty) \times \mathbb{R} \times (0,\infty)\big)\) as \[\begin{align} \Bar{\mathcal{L}}\left(m_0,\sigma_0^2,m_1,\sigma_1^2\right)&=\begin{cases} \mathcal{L}\left(m_0,\sigma_0^2,m_1,\sigma_1^2\right)&\textrm{ if } \sigma_0> 0,\, \sigma_1> 0,\\[6pt] \lim_{t \to 0}\mathcal{L}\left(m_0,\sigma_0^2,m_1,t\right)&\textrm{ if } \sigma_0> 0,\, \sigma_1 = 0,\\[6pt] \lim_{s \to 0}\mathcal{L}\left(m_0,s,m_1,\sigma_1^2\right)&\textrm{ if } \sigma_0= 0,\,\sigma_1> 0 \end{cases}\nonumber \\ &=\begin{cases} \mathcal{L}\left(m_0,\sigma_0^2,m_1,\sigma_1^2\right)&\textrm{ if } \sigma_0> 0,\, \sigma_1> 0,\\[6pt] \frac{1}{4}\left(\sigma_0^2 + \left(m_1-m_0 \right)^2\right)&\textrm{ if } \sigma_0> 0,\, \sigma_1 = 0,\\[6pt] \infty&\textrm{ if } \sigma_0= 0,\,\sigma_1> 0. \label{eq:lsc95extension95univariate} \end{cases} \end{align}\tag{25}\] Furthermore, Theorem 7 also shows that \(\Bar{\mathcal{L}}\) admits a lower semi-continuous extension if we define the value of the function for \(\sigma_1=\sigma_0=0\) to be equal to \(\frac{1}{4}(m_1-m_0)^2\).

This is depicted in Figure 5 and is summarized in the next corollary.

We finally note that \(D_{\rm WKL}\) is continuous at \(\Sigma_0=\Sigma_1=\Sigma\succ 0\) and takes a value independent of \(\Sigma\) (see 21 ). This allows us to take limit as \(\Sigma\) tends to \(0\) to approximate the divergence between Dirac measures concentrated at \(m_0\) and \(m_1\) while providing a finite value that is proportional to the squared distance \(\lVert m_0-m_1\rVert^2\). In contrast, the KL-divergence for \(\Sigma_0=\Sigma_1=\Sigma\) gives \(D_{\rm KL}(\nu\lVert \mu)=\frac{1}{2}\left\lVert \Sigma^{-1/2}(m_1-m_0)\right\rVert^2\) which diverges to \(\infty\) as \(\Sigma\) approaches singularity. Figure 6 shows the surface plots for the KL-divergence and the WKL-divergence for univariate Gaussian distributions.

4 Outlook↩︎

In this work, we introduced the WKL-divergence for multivariate Gaussian distributions, building on the framework of [5] and analyzed its continuity properties. Promising directions for future research include exploring its role in information-theoretic tasks such as maximum likelihood estimation and statistical inference. Another important direction is a systematic comparison of its empirical performance against classical divergences on real-world machine learning tasks. We believe the WKL-divergence provides a robust alternative to KL-based methods, with potential to enrich both theoretical developments and practical applications in information geometry and beyond.

5 Reduction of the Canonical Divergence to KL-Divergence↩︎

In this appendix we prove that, when the Fisher–Rao metric and the \(e_0\)-geodesic are used, the canonical divergence reduces exactly to the KL-divergence, i.e., \[\begin{align} D^{(\rm e_0)}(\mu \| \nu) &:=\int_0^1 t \langle\dot{\gamma}_0(t),\dot{\gamma}_0(t) \rangle_{\gamma_0(t)}^{\rm FR} dt= D_{\rm KL}(\nu \| \mu). \end{align}\] Recall that \[\begin{align} \gamma_0(t)=\frac{\left(\frac{d\nu}{d\mu}\right)^t}{Z(t)} \mu=\rho_t \cdot \mu, \textrm{ where }Z(t)=\int_{\mathbb{R}^n} \left(\frac{d\nu}{d\mu}\right)^t d\mu. \end{align}\] Note that if \(\mu=f_{\mu}\lambda\) and \(\nu=f_{\nu}\lambda\), then \[\begin{align} Z(t)=\int_{\mathbb{R}^n} \left(\frac{d\nu}{d\mu}\right)^td\mu=\int_{\mathbb{R}^n} \left(\frac{f_{\nu}}{f_{\mu}}\right)^tf_{\mu}d\lambda=\int_{\mathbb{R}^n} \left(\frac{f_{\nu}}{f_{\mu}}\right)^tf_{\mu}d\lambda=\int_{\mathbb{R}^n} f_{\nu}^tf_{\mu}^{(1-t)}d\lambda. \end{align}\] Note that \(f_{\nu}^tf_{\mu}^{(1-t)}\) is differentiable with respect to \(t\) and \[\begin{align} \frac{d}{dt} \left(f_{\nu}^tf_{\mu}^{(1-t)}\right)=\left(f_{\nu}^tf_{\mu}^{(1-t)}\right)\ln \! \left(\frac{f_{\nu}}{f_{\mu}}\right). \end{align}\] It can be shown that \(f_{\nu}^tf_{\mu}^{(1-t)}\) is a Gaussian density for any \(t\in[0,1]\) and it decays exponentially as \(\lVert x \rVert\) grows whereas \(\ln \left(\frac{f_{\nu}}{f_{\mu}}\right)\) has at most polynomial growth. This justifies differentiating under the integral sign yielding \[\begin{align} \frac{d}{dt}Z(t)=\int_{\mathbb{R}^n} \left(f_{\nu}^tf_{\mu}^{(1-t)}\right)\ln \left(\frac{f_{\nu}}{f_{\mu}}\right)d\lambda=\int_{\mathbb{R}^n} \left(\frac{f_{\nu}}{f_{\mu}}\right)^t\ln \left(\frac{f_{\nu}}{f_{\mu}}\right)d\mu=\int_{\mathbb{R}^n} \left(\frac{d\nu}{d\mu}\right)^t\ln \left(\frac{d\nu}{d\mu}\right)d\mu. \end{align}\] This gives us \[\begin{align} \frac{1}{Z(t)}\frac{d}{dt}Z(t) &=\frac{1}{\int_{\mathbb{R}^n} \left(\frac{d\nu}{d\mu}\right)^td\mu} \int_{\mathbb{R}^n} \left(\frac{d\nu}{d\mu}\right)^t \ln \left(\frac{d\nu}{d\mu}\right)d\mu =\int_{\mathbb{R}^n} \ln \left(\frac{d\nu}{d\mu}\right)d\gamma_0(t). \end{align}\] Differentiating \(\gamma_0\) with respect to time and using the above expression for \(\frac{\dot{Z}(t)}{Z(t)}\), we get that \[\begin{align} \dot{\gamma}_0(t)=\left(\frac{\left(\frac{d\nu}{d\mu}\right)^t \ln \left(\frac{d\nu}{d\mu}\right)}{Z(t)}-\frac{\left(\frac{d\nu}{d\mu}\right)^t \dot{Z}(t)}{Z(t)^2} \right)\mu&=\left(\ln \left(\frac{d\nu}{d\mu}\right)-\frac{\dot{Z}(t)}{Z(t)} \right)\gamma_0(t)\\ &=\left(\ln \left(\frac{d\nu}{d\mu}\right)-\int_{\mathbb{R}^n} \ln \left(\frac{d\nu}{d\mu}\right)d\gamma_0(t) \right)\gamma_0(t)\\ &=\left(\ln \left(\frac{d\nu}{d\mu}\right) - \mathbb{E}_{\gamma_0(t)}\left[\ln \left(\frac{d\nu}{d\mu}\right)\right] \right) \gamma_0(t). \end{align}\] Plugging this into 8 and using the definition of the Fisher\(-\)Rao metric gives \[\begin{align} D^{(\rm e_0)}(\mu \| \nu) &= \int_0^1 t \langle \dot{\gamma}_0(t), \dot{\gamma}_0(t)\rangle^{\rm FR}_{\gamma_0(t)} dt \nonumber \\ &= \int_0^1 t \left(\int_{\mathbb{R}^n} \left(\ln \left(\frac{d\nu}{d\mu}\right) - \mathbb{E}_{\gamma_0(t)}\left[\ln \left(\frac{d\nu}{d\mu}\right)\right] \right)^2 d\gamma_0(t)\right) dt \nonumber \\ &= \int_0^1 t \cdot \mathbb{E}_{\gamma_0(t)}\left[\left(\ln \left(\frac{d\nu}{d\mu}\right) - \mathbb{E}_{\gamma_0(t)}\left[\ln \left(\frac{d\nu}{d\mu}\right)\right] \right)^2\right] dt. \nonumber \\ \end{align}\] Another justified application of the differentiation under the integral sign yields the standard identity \[\begin{align} \frac{d}{dt}\left(\mathbb{E}_{\gamma_0(t)}\left[\ln \frac{d\nu}{d \mu}\right]\right)=\mathbb{E}_{\gamma_0(t)}\left[\left(\ln \left(\frac{d\nu}{d\mu}\right) - \mathbb{E}_{\gamma_0(t)}\left[\ln \left(\frac{d\nu}{d\mu}\right)\right] \right)^2\right]. \end{align}\] This, together with the application of integration by parts gives us \[\begin{align} \label{eq:KL95energy95formula95specialize95to95KL95intermediate} D^{(\rm e_0)}(\mu \| \nu) &= \int_0^1 t \cdot \frac{d}{dt}\left(\mathbb{E}_{\gamma_0(t)}\left[\ln \frac{d\nu}{d \mu}\right]\right) dt \nonumber \\ &= \mathbb{E}_{\gamma_0(1)}\left[\ln \frac{d\nu}{d \mu}\right]-\int_0^1 \left(\mathbb{E}_{\gamma_0(t)}\left[\ln \frac{d\nu}{d \mu}\right]\right) dt. \end{align}\tag{26}\] Finally, using the fact that \(\gamma_0(0)=\mu\) and \(\gamma_0(1)=\nu\) and \(\frac{\dot{Z}(t)}{Z(t)}=\int_{\mathbb{R}^n} \ln \left(\frac{d\nu}{d\mu}\right)d\gamma_0(t)\), we get that \[\begin{align} \int_0^1 \left(\mathbb{E}_{\gamma_0(t)}\left[\ln \frac{d\nu}{d \mu}\right]\right) dt &= \int_0^1 \frac{\dot{Z}(t)}{Z(t)} dt= \int_0^1 \frac{d}{dt}\ln (Z(t)) dt=\ln(Z(1))-\ln (Z(0))=0 \end{align}\] which finally gives us \[\begin{align} D^{(\rm e_0)}(\mu \| \nu) &= \mathbb{E}_{\nu}\left[\ln \frac{d\nu}{d \mu}\right]=\int_{\mathbb{R}^n}\ln \frac{d\nu}{d \mu} d\nu = D_{\rm KL}(\nu \| \mu). \end{align}\]

6 Characterization of Quadratic Potentials for Tangent Vectors↩︎

In this appendix, we show that for any tangent vector \(v\) in \(T_{\mu}\mathcal{G}\), there exists a unique quadratic function \(f_v \in C^{\infty}(\mathbb{R}^n)\), uniquely determined up to an additive constant, such that \(v={\rm div_{\mu}} (\nabla{f_v})\mu\). We restrict attention to quadratic functions \(f_v:\mathbb{R}^n\rightarrow \mathbb{R}\) of the form \(f_v(x)=\frac{1}{2}x^TAx+b^Tx\) such that its gradient vector field is given by \(\nabla{f_v}=Ax+b\). Since \(\nabla{f_v}\) is unchanged under \(f_v \mapsto f_v + c\), the quadratic potential is unique up to an additive constant. Using properties of Gaussian random variables, it can be shown that \[\begin{align} x_0 \sim \mu=\mathcal{N}(m_0,\Sigma_0) \implies x(t) \sim \mu_t=\mathcal{N}(\underbrace{e^{At}{m}_0+\left(e^{At}A^{\dagger}+tA^{\perp}-A^{\dagger}\right)b}_{m_t},\underbrace{e^{At}\Sigma_0 e^{At}}_{\Sigma_t}). \end{align}\] This gives us the density \({\rm div_{\mu}} (\nabla{f_v})\) of \(\frac{d}{dt}\mu_t\big|_{t=0} \in T_{\mu} \mathcal{G}\) with respect to \(\mu\) in coordinates as \[\begin{align} \frac{d}{dt}m_t\big|_{t=0}&=Am_0+(AA^{\dagger}+A^{\perp})b=Am_0+b,\\ \frac{d}{dt}\Sigma_t|_{t=0}&=A\Sigma_0+\Sigma_0A. \end{align}\] Now consider an arbitrary tangent vector in \(T_{\mu} \mathcal{G}\) described in coordinates by \((\dot{m},\dot{\Sigma}) \in \mathbb{R}^n\times S_n\). The vector field \(f_v\) whose gradient flow realizes the desired tangent vector can be obtained by solving \[\begin{align} A\Sigma_0+\Sigma_0A&=\dot{\Sigma},\\ Am_0+b&=\dot{m} \end{align}\] for \(A\in S_n\) and \(b \in \mathbb{R}^n\). The Lyapunov equation \[\begin{align} A\Sigma_0+\Sigma_0A&=\dot{\Sigma} \end{align}\] possesses a unique symmetric solution if \(\Sigma_0\succ 0\). To see this explicitly, write \(\Sigma_0=U\Lambda U^T\) (the orthogonal eigenvalue decomposition of \(\Sigma_0\)). Multiplying the above equation from the left by \(U^T\) and from the right by \(U\), we obtain the equation \[\begin{align} \underbrace{U^TAU}_B\Lambda+\Lambda \underbrace{U^TAU}_B&=\underbrace{U^T\dot{\Sigma} U}_C. \end{align}\] This yields the explicit solution \(B_{ij}=\frac{1}{\lambda_i+\lambda_j}C_{ij}\) which is well-defined since \(\lambda_i + \lambda_j>0\) owing to the positive definiteness of \(\Sigma_0\). This uniquely determines \(A=UBU^T\). Once \(A\) is determined, \(b=\dot{m} - Am_0\) follows directly from the mean constraint. Thus, for any tangent vector in \(T_{\mu} \mathcal{G}\), there exists a unique quadratic function \(f_v\) (unique up to an additive constant) such that \(v={\rm div_{\mu}} (\nabla{f_v})\mu\).

7 Reduction of the Canonical Divergence to WKL-Divergence↩︎

In this appendix we prove that, when the Otto metric and the \(e_1\)-geodesic are used, the canonical divergence reduces exactly to the Wasserstein KL-divergence formula 10 , i.e., \[\begin{align} D^{(\rm e_1)}(\mu \| \nu):=\int_0^1 t \langle\dot{\gamma}_1(t),\dot{\gamma}_1(t) \rangle_{\gamma_1(t)}^{\rm O} dt=\int_{\mathbb{R}^n} \int_0^1 \left(f \circ \varphi_1 - f \circ \varphi_t\right) \, dt \, d\mu. \end{align}\] Differentiating \(\gamma_1(t)=\rho_t \mu\) with respect to time, we get that \[\begin{align} \dot{\gamma}_1(t)=\dot{\rho}_t \mu = -{\rm div_{\gamma_1(t)}}(\nabla{f})\rho_t \mu = -{\rm div_{\gamma_1(t)}}(\nabla{f})\gamma_1(t). \end{align}\] By the definition of the Otto metric, \[\begin{align} \langle \dot{\gamma}_1(t),\dot{\gamma}_1(t) \rangle_{\gamma_1(t)}^{\rm O} :&= \int_{\mathbb{R}^n} \langle -\nabla{f}, -\nabla{f} \rangle d\gamma_1(t)= \int_{\mathbb{R}^n} \langle \nabla{f}, \nabla{f} \rangle d\gamma_1(t). \end{align}\] Recall that the pushforward relation \(\gamma_1(t) = (\varphi_t)_* \mu\) implies that \[\begin{align} \int_{\mathbb{R}^n}h(y)d\gamma_1(t)(y)=\int_{\mathbb{R}^n}h(\varphi_t(x))\mu(dx) \end{align}\] holds for any measurable function \(h\). Hence, \[\begin{align} \langle \dot{\gamma}_1(t),\dot{\gamma}_1(t) \rangle_{\gamma_1(t)}^{\rm O} = \int_{\mathbb{R}^n} \langle \nabla{f}, \nabla{f} \rangle d\gamma_1(t)&= \int_{\mathbb{R}^n} \langle \nabla{f}(\varphi_t(x)),\nabla{f}(\varphi_t(x))\rangle \mu(dx)\\ &=\int_{\mathbb{R}^n} \langle \nabla{f}(\varphi_t(x)),\frac{d}{dt}\varphi_t(x)\rangle \mu(dx)\\ &=\int_{\mathbb{R}^n} \frac{d}{dt}f(\varphi_t(x)) \mu(dx) \end{align}\] where we have used gradient flow dynamics \(\tfrac{d}{dt}\varphi_t(x) = \nabla{f}(\varphi_t(x))\) and the chain rule of differentiation. Plugging this into 8 gives us \[\begin{align} \label{eq:KL95energy95formula95specialize95to95WKL} D^{(\rm e_1)}(\mu \| \nu) &= \int_0^1 t \langle \dot{\gamma}_1(t), \dot{\gamma}_1(t)\rangle^{\rm O}_{\gamma_1(t)} dt \nonumber \\ &= \int_{\mathbb{R}^n} \left(\int_0^1 t \frac{d}{dt}(f(\varphi_t(x))) dt\right) \mu(x) \nonumber \\ &= \int_{\mathbb{R}^n} \int_0^1 \left(f \circ \varphi_1 - f \circ \varphi_t\right) \, dt \, d\mu, \end{align}\tag{27}\] where we have exchanged the order of integration. Since \(f\) is quadratic and \(\mu\) is Gaussian, all integrals are finite, and Fubini’s theorem justifies this. This completes the reduction, showing that the canonical divergence with the Otto metric recovers the WKL-divergence formula 10 .

References↩︎

---

1. Note that the assumption on the symmetry of \(A\) is without loss of generality since replacing \(A\) by \(\frac{1}{2}(A+A^T)\) keeps the function unchanged.↩︎

2. We use the fact that \(A^{\perp}\), \(A^{\dagger}\) and \(e^A\) can be simultaneously diagonalized by orthogonal matrices owing to the symmetry of \(A\) which leads to a number of useful properties such as commutativity. These are used throughout the paper.↩︎

3. We mainly use the linearity and the cyclic property of the trace operators↩︎