How to calculate the barycenter of a weighted graph
Sébastien Gadat
Toulouse School of Economics, Université Toulouse 1 Capitole sebastien.gadat@math.univ-toulouse.fr,
http://perso.math.univ-toulouse.fr/gadat/

arXiv:1605.04148v1 [math.PR] 13 May 2016

Ioana Gavra
Institut de Mathématiques de Toulouse, Université Toulouse 3 Paul Sabatier Ioana.Gavra@math.univ-toulouse.fr,

Laurent Risser
Institut de Mathématiques de Toulouse, CNRS lrisser@math.univ-toulouse.fr, http://laurent.risser.free.fr/

Discrete structures like graphs make it possible to naturally and flexibly model complex phenomena. Since
graphs that represent various types of information are increasingly available today, their analysis has become
a popular subject of research. The graphs studied in the field of data science at this time generally have
a large number of nodes that are not fairly weighted and connected to each other, translating a structural
specification of the data. Yet, even an algorithm for locating the average position in graphs is lacking
although this knowledge would be of primary interest for statistical or representation problems. In this work,
we develop a stochastic algorithm for finding the Fréchet mean of weighted undirected metric graphs. This
method relies on a noisy simulated annealing algorithm dealt with using homogenization. We then illustrate
our algorithm with two examples (subgraphs of a social network and of a collaboration and citation network).
Key words : metric graphs; Markov process; simulated annealing; homogeneization
MSC2000 subject classification : 60J20,60J60,90B15
OR/MS subject classification : Primary: Probability - Markov processes; secondary: Networks/graphs Traveling salesman
History :

1.

Introduction.

Numerous open questions in a very wide variety of scientific domains involve complex discrete
structures that are easily modeled by graphs. The nature of these graphs may be weighted or not,
directed or not, observed online or by using batch processing, each time implying new problems
and sometimes leading to difficult mathematical or numerical questions. Graphs are the subject of
perhaps one of the most impressive growing bodies of literature dealing with potential applications
in statistical or quantum physics (see, e.g., [21]) economics (dynamics in economy structured as
networks), biology (regulatory networks of genes, neural networks), informatics (Web understanding and representation), social sciences (dynamics in social networks, analysis of citation graphs).
We refer to [37] and [43] for recent communications on the theoretical aspects of random graph
models, questions in the field of statistics and graphical models, and related numerical algorithms.
In [34], the authors have developed an overview of numerous possible applications using graphs
and networks in the fields of industrial organization and economics. Additional applications, details
and references in the field of machine learning may also be found in [28].
In the meantime, the nature of the mathematical questions raised by the models that involve
networks is very extensive and may concern geometry, statistics, algorithms or dynamical evolution
over the network, to name a few. For example, we can be interested in the definition of suitable
random graph models that make it possible to detect specific shape phenomena observed at different scales and frameworks (the small world networks of [47], the existence of a giant connected
component for a specific range of parameters in the Erdös-Rényi random graph model [20]). A complete survey may be found in [39]. Another important field of investigation is dedicated to graph
visualization (see some popular methods in [45] and [35], for example). In statistics, a popular topic
1

2

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

deals with the estimation of a natural clustering when the graph follows a specific random graph
model (see, among others, the recent contribution [36] that proposes an optimal estimator for the
stochastic block model and smooth graphons). Other approaches rely on efficient algorithms that
analyze the spectral properties of adjacency matrices representing the networks (see, e.g, [6]). A
final important field of interest deals with the evolution of a dynamical system defined through
a discrete graph structure: this is, for instance, the question raised by gossip models that may
describe belief evolution over a social network (see [1] and the references therein).
We address a problem here that may be considered as very simple at first glance: we aim to
define and estimate the barycenter of weighted graphs. Hence, this problem involves questions
that straddle the area of statistics and the geometry of graphs. Surprisingly, as far as we know,
this question has received very little attention, although a good assessment of the location of a
weighted graph barycenter could be used for fair representation issues or for understanding the
graph structure from a statistical point of view. In particular, it could be used and extended to
produce “second order” moment analyses of graphs. We could therefore generalize a Principal
Component Analysis by extending this framework. This intermediary step was used, for example,
by [13] to extend the definition of PCA on the space of probability measures on R, and by [14]
to develop a suitable geometric PCA of a set of images. A popular strategy to define moments
in complex metric spaces is to use the variational interpretation of means (or barycenters), which
leads to the introduction of Fréchet (or Karcher) means (see Section 2.2 for an accurate definition).
This approach has been introduced in the seminal contribution [23] that makes it possible to define
p-means over any metric probability space.
The use of Fréchet means has met with great interest, especially in the field of bio-statistics and
signal processing, although mathematical and statistical derivations around this notion constitute
a growing field of interest.
• In continuous domains, many authors have recently proposed limit theorems on the empirical
Fréchet sample mean (only a sample of size n of the probability law is observed) towards its
population counterpart. These works were mostly guided by applications to continuous manifolds
that describe shape spaces introduced in [19]. For example, [38] establishes the consistency of the
population Fréchet mean and derives applications in the Kendall space. The study of [9] establishes
the consistency of Fréchet empirical means and derives its asymptotic distribution when n −→ +∞,
whereas the observations live in more general Riemannian manifolds. Finer results can be obtained
in some non-parametric restrictive situations (see, e.g., [11, 12, 15] dedicated to the so-called shape
invariant model ). Many applications in various domains involving signal processing can also be
found: ECG curve analysis [10] and image analysis [44, 2], to name a few.
• Recent works treat Fréchet means in a discrete setting, especially when dealing with phylogenetic trees that have an important hierarchical structure property. In particular, [8] proposed
a central limit theorem in this discrete case, whereas [41] used an idea of Sturm for spaces with
non-positive curvature to define an algorithm for the computation of the population Fréchet mean.
Other works deal with the averaging of discrete structure sequences such as diagrams using the
Wasserstein metric (see, e.g., [42]) or graphs [27].
Our work deviates from the above-mentioned point since we build an algorithm that recovers the
Fréchet mean of a weighted graph while observing an infinite sequence of nodes of the graph,
instead of finding the population Fréchet mean of a set of discrete structures, as proposed in
[8, 41, 42, 27]. Our algorithm uses recent contributions on simulated annealing ([3, 4]). It relies on
a continuous-time noisy simulated annealing Markov process on graphs, as well as a second process
that accelerates and homogenizes the updates of the noisy transitions in the simulated annealing
procedure.
The paper is organized as follows: Section 2 highlights the different difficulties raised by the
computation of the Fréchet mean of weighted graphs and describes the algorithm we propose. In

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

3

Section 3.1, theoretical backgrounds to understand the behavior of our method are developed on
the basis of this algorithm, and Section 3.2 states our main results. Simulations and numerical
insights are then given in Section 4. The convergence of the algorithm is theoretically established
in Section 5, whereas Section 6 describes functional inequalities in quantum graphs that will be
introduced below.

Acknowledgments.
The authors gratefully acknowledge Laurent Miclo for his stimulating discussions and helpful comments throughout the development of this work, and Nathalie Villa-Vialaneix for her interest and
advice concerning simulations. The authors are also indebted to zbMATH for making their database
available to produce numerical simulations.

2.

Stochastic algorithm on quantum graphs

We consider G = (V, E) a finite connected and undirected graph with no loop, where V = {1, . . . , N }
refers to the N vertices (also called nodes) of G, and E the set of edges that connect some couples
of vertices in G.
2.1. Undirected weighted graphs
The structure of G may be described by the adjacency matrix W that gives a non-negative weight
to each edge E (pair of connected vertices), so that W = (wi,j )1≤i,j≤N while wi,j = +∞ if there is
no direct link between node i and node j. W indicates the length of each direct link in E: a small
positive value of wi,j represents a small length of the edge {i, j }. We assume that G is undirected
(so that the adjacency matrix W is symmetric) and connected: for any couple of nodes in V , we
can always find a path that connects these two nodes. Finally, we assume that G has no self loop.
Hence, the matrix W satisfies:
∀i 6= j

wi,j = wj,i

and

∀i ∈ V

wi,i = 0

We define d(x, y) as the geodesic distance between two points (x, y) ∈ V 2 , which is the length of
the shortest path between them. The length of this path is given by the addition of the length of
traversed edges:
k−1
X
2
∀(i, j) ∈ V
d(i, j) =
min
wi` ,i`+1 .
i=i1 →i2 →...→ik+1 =j

`=1

When the length of the edges is constant and equal to 1, it simply corresponds to the number of
traversed edges. Since the graph is connected and finite, we introduce the definition of the diameter
of G:
DG := sup d(x, y).
(x,y)∈V 2

To define a barycenter of a graph, it is necessary to introduce a discrete probability distribution
ν over the set of vertices V . This probability distribution is used to measure the influence of each
node on the graph.
Example 1. Let us consider a simple scientometric example illustrated in Figure 1 and consider
a “toy” co-authorship relation that could be obtained in a subgraph of a collaboration network like
zbMATH 1 . If two authors A and B share kA,B joint papers, it is a reasonable choice to use a weight
wA,B = φ(kA,B ), where φ is a convex function satisfying φ(0) = +∞, φ(1) = 1 and φ(+∞) = 0. This
means that no joint paper between A and B leads to the absence of a direct link between A and B
1

https://zbmath.org/

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

4

1/10
ν(2) = 0.1

ν(3) = 0.1
1/20

1/5

ν(5) = 0.1
1/30

ν(1) = 0.5

ν(4) = 0.2
1
1/10

Figure 1

Example of a weighted graph between five authors, where the probability mass is {0.5, 0.1, 0.1, 0.2, 0.1}.
In this graph, Author 1 shares five publications with Author 2 and ten publications with Author 5, and
so on.

on the graph. On the contrary, the more papers there are between A and B, the closer A and B will
be on the graph. Of course, this graph may be embedded in a probability space with the additional
definition of a probability distribution over the authors that can be naturally proportional to the
number of citations of each author. This is a generalization of the Erdös graph. Note that this type
of example can also be encountered when dealing with movies and actors, leading, for example, to
the Bacon number and graph (see the website: www.oracleofbacon.org/).
2.2. Fréchet mean of an undirected weighted graph
Following the simple remark that the p-mean of any distribution ν of Rd is the point that minimizes
x 7−→ EZ∼ν [|x − Z |p ],
it is legitimate to be interested in the Fréchet mean of a graph (G, ν) where ν refers to the
probability distribution over each node. The Fréchet mean is introduced in [23] to generalize this
variational approach to any metric space. Although we have chosen to restrict our work to the
case of p = 2, which corresponds to the Fréchet mean definition, we believe that our work could be
generalized to any value of p ≥ 1. If d denotes the geodesic distance w.r.t. G, we are interested in
solving the following minimization problem:
Z
1
d2 (x, y)ν(y).
(1)
Mν := arg min Uν (x)
where
Uν (x) =
x∈E
2 G
Hence, the Fréchet mean of (G, ν), denoted Mν , is the set of all possible minimizers of Uν . Note
that this set is not necessarily a singleton and this uniqueness property generally requires some
additional topological assumptions (see [3], for example). At this point, we want to make three
important remarks about the difficulty of this optimization problem:
• The problem of finding Mν involves the minimization of Uν , which is a non-convex function
with the possibility of numerous local traps. To our knowledge, this problem cannot be efficiently
solved using either a relaxed solution or using a greedy/dynamic programing algorithm (in the
spirit of the Dijkstra method that makes it possible to compute geodesic paths [18]).

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

5

• It is thus natural to think about the use of a global minimization procedure, and, in particular,
the simulated annealing (S.A. for short) method. S.A. is a standard strategy to minimize a function
over discrete spaces and its computational cost is generally high. It relies on an inhomogeneous
Markov process that evolves on the graph with a transition kernel depending on the energy estimates Uν . In other words, it requires the computation of Uν that depends on an integral w.r.t. ν.
Since we plan to handle large graphs, this last dependency can be a very strong limitation.
• In some cases, even the global knowledge of ν may not be realistic, and the importance of
each node can only be revealed through i.i.d. sequential arrivals of new observations in E that
are distributed according to ν. This may be the case, for example, if we consider a probability
distribution over E that is cropped while gathering interactive forms on a website.

2.3. Outline of Simulated Annealing (S.A.)
The optimization with S.A. introduces a Markov random dynamical system that evolves either in
continuous time (generally for continuous spaces) or in discrete time (for discrete spaces). When
dealing with a discrete setting, S.A. is based on a Markov kernel proposition L(., .) : V × V −→
[0, 1] related to the 1-neighborhoods of the Markov chain, as introduced in [30]. It is based on
an inhomogeneous Metropolis-Hastings scheme, which is recalled in Algorithm 1. We can derive
asymptotic guarantees of the convergence towards a minimum of U as soon as the cooling schedule
is well chosen.
Algorithm 1: M.-H. Simulated Annealing
Data: Function U . Decreasing temperature sequence (Tk )k≥0
1Initialization: X0 ∈ V ;
2for k = 0 . . . + ∞ do
o
n
−1

0

0

,Xk )
Draw x0 ∼ L(Xk , .) and compute pk = 1 ∧ eTk [U (Xk )−U (x )] L(x
L(Xk ,x0 )
(
x0 with probability pk
Update Xk+1 according to Xk+1 =
Xk with probability 1 − pk

3
4
5end

6Output:

limk−→+∞ Xk

When dealing with a continuous setting, S.A. uses a drifted diffusion with a vanishing variance
(t )t≥0 over V , or an increasing drift coefficient (βt )t≥0 . We refer to [32, 40] for details and we recall
its Langevin formulation in Algorithm 2.
Algorithm 2: Langevin Simulated Annealing
Data: Function U . Increasing inverse temperature (βt )t≥0
1 Initialization: X0 ∈ V ;
2
∀t ≥ 0
dXt = −βt ∇U (Xt )dt + dBt
3Output: limt−→+∞ Xt
In both cases, we can see that S.A. with U = Uν given by (1) involves the computation of the
value of U in Line 4 of Algorithm 1, or the computation of ∇U in Line 2 of Algorithm 2. These two
computations are problematic for our Fréchet mean problem: the integration over ν is intractable
in the situation of large graphs and we are naturally driven to consider a noisy version of S.A. A
possible alternative method for this problem is to use a homogenization technique: replacing Uν in
the definition of pk by Uy (·) = d2 (·, y), where y is a value from an i.i.d. sequence (Yn )n≥0 distributed
according to ν . Such methods have been developed in [29] as a modification of Algorithm 1 with
an additional Monte-Carlo step in Line 4, when UY (x) follows a Gaussian distribution centered
around the true value of U (x) = EY ∼ν UY (x). This approach is still problematic in our case since
the Gaussian assumption on the random variable d2 (x, Y ) is unrealistic here. Another limitation
of this MC step relies on the fact that it requires a batch average of several (U (x, Yj ))1≤j≤nk where

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

6

nk is the number of observations involved at iteration k, although we also plan to develop an
algorithm that may be adapted to on-line arrivals of the observation (Yn )n≥0 . Lastly, it is important
to observe that the non-linearity of the exponential prevents the use of only one observation Yk in
the acceptation/reject ratio involved in the S.A. since it does not lead to an unbiased evaluation
of the true transition:
h −1
i
−1
0
0
EY ∼ν eTk [UY (Xk )−UY (x )] 6= eTk EY ∼ν [UY (Xk )−UY (x )] .
This difficulty does not arise in the homogenization of the simulated annealing algorithm in the
continuous case since the exponential is replaced by a gradient, i.e., the process we use is a Markov
process of the form:
dXt = −βt ∇UYt (Xt )dt + dBt ,
where Bt is a “Brownian motion” on the graph G, βt−1 refers to the inverse of the temperature, and
Yt is a continuous time Markov process obtained from the sequence (Yn )n∈N . This point motivates
the introduction of the quantum graph induced by the initial graph. Of course, dealing with a
continuous diffusion over a quantum graph G deserves special theoretical attention, which will be
given in Section 3.1.
2.4. Homogenized S.A. algorithm on a quantum graph
We now present the proposed algorithm for estimating Fréchet means. To do so, we first introduce
the quantum graph ΓG derived from G = (V, E) that corresponds to the set of points living inside
the edges e ∈ E of the initial graph. Once an orientation is arbitrarily fixed for each edge of E, the
location of a point in ΓG depends on the choice of an edge e ∈ E and on a coordinate xe ∈ [0, Le ]
where Le is the length of edge e on the initial graph. The coordinate 0 then refers to the initial
point of e and Le refers to the other extremity.
While considering the quantum graph ΓG , it is still possible to define the geodesic distance
between any point x ∈ ΓG and any node y ∈ G. In particular, when x ∈ V , we use the initial definition of the geodesic distance over the discrete graph, although when x ∈ e ∈ E with a coordinate
xe ∈ [0, Le ], the geodesic distance between x and y is:
d(x, y) = {xe + d(e(0), y)} ∧ {Le − xe + d(e(Le ), y)} .
This definition can be naturally generalized to any two points of ΓG , enabling us to consider the
metric space (ΓG , d).
Consider a positive, continuous and increasing function t 7−→ αt such that:
lim αt = +∞

t7−→+∞

∀t ≥ 0

and

βt = o(αt ).

We introduce (Ntα )t≥0 , an inhomogeneous Poisson process over R+ with intensity α. It is standard
to represent N α through a homogeneous Poisson process H of intensity 1 using the relationship:
Z t
α
∀t ≥ 0
Nt = Hh(t) , where h(t) =
αs ds.
0

Using this accelerated process N α , our optimization algorithm over ΓG is based on the Markov
process that solves the following stochastic differential equation over ΓG :
(
X0 ∈ ΓG
(2)
dXt = −βt ∇UYN α (Xt )dt + dBt .
t

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

7

Using the definition of Nt and the basic properties of a Poisson Process, it can be observed that
for all  > 0 and all t ≥ 0:


Z
Nt+ − Nt
1 t+
E
=
αs ds

 t
Hence, α should be understood as the speed of new arrivals in the sequence (Yn )n∈N .
Our definition (2) is slightly inaccurate: for any y ∈ V (G), the function Uy is C 1 except at a finite
set of points of ΓG :
• This can be the case inside an edge x ∈ e with xe ∈ (0, Le ) when at least two different geodesic
paths from x to y start in opposite directions. Then, we define ∇Uy (x) = 0.
• This can also be the case at a node x ∈ V when several geodesic paths start from x to y. In
that case, we once again arbitrarily impose a null value for the “gradient” of Uy at node x. This
prior choice will not have any influence on the behavior of the algorithm.
We will first detail below the theoretical objects involved in (2), then we will describe an efficient
discretization of (2) that makes it possible to derive our practical optimization algorithm.
Algorithm 3: Homogenized Simulated Annealing over a quantum graph
Data: Function U . Increasing inverse temperature (βt )t≥0 . Intensity (αt )t≥0
1Initialization: Pick X0 ∈ ΓG . ;
2T0 = 0 ;
3for k = 0 . . . + ∞ do
4
while Ntα = k do
5
Xt evolves as a Brownian motion, relatively to the structure of ΓG initialized at XT−k .
6
end
7
Tk+1 := inf {t : Ntα = k + 1};
8
At time t = Tk+1 , draw Yk+1 = YNtα according to ν.
9
The process Xt jumps from Xt− towards Yk+1 :
−
−−−−−−
→

Xt = Xt− + βt αt−1 Xt YNtα ,

(3)

−−−−→
where Xt YNtα represents the shortest (geodesic) path from Xt to YNtα in ΓG .
10end
11Output: limt−→+∞ Xt

Algorithm 3 could be studied following the road map of [5]. Nevertheless, this implies serious
regularity difficulties on the densities and the Markov semi-group involved. Hence, we have chosen
to consider Algorithm 3 as a natural Euler explicit discretization of our Markov evolution (2): for
a large value of k, the average time needed to travel from Tk to Tk+1 is approximately αT−1
−→ 0 as
k
k −→ +∞. On this short time interval, the drift term in (2) is the gradient of the squared geodesic
−−−→
distance between XTk and YNTα , which is approximated by our vector Xt Yk , multiplied by βTk ,
k
leading to (3). Xt now evolves as a Brownian motion over ΓG between two jump times and this
evolution can be simulated with a Gaussian random variable using a (symmetric) random walk
when the algorithm hits a node of ΓG . Figure 2 proposes a schematic evolution of (Xt )t≥0 over a
simple graph ΓG with five nodes. We will prove the following result.
Theorem 1. A constant c? (Uν ) exists such that if αt = 1 + t and βt = b log(1 + t) with b >
{c (Uν )}−1 , then (Xt )t≥0 defined in (2) converges a.s. to Mν defined in (1), when t goes to infinity.
?

A more rigorous form of this result and more details about the process defined in (2) are presented
in Section 3.2. The proof of Theorem 3 is deferred to Section 5.

3.

Inhomogeneous Markov process over ΓG

This section presents the theoretical background needed to define the Markov evolution (2).

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

8

−
Jump of size {βT α−1 } × d(XT
, YT1 )
1
1 T
1

V2
XT1 +
Brownian displacement

XT−2 +
V1

V3
+
XT−1
−
Jump of size {βT α−1 } × d(XT
, Y T2 )
2
2 T2

V5 YT2

XT2
+
V4

YT1
Figure 2

Schematic evolution of the Homogenized S.A. described in Algorithm 3 over the quantum graph. We
first observe a jump at time T1 towards YT1 = V1 and a Brownian motion on ΓG during T2 − T1 . A
second node is then sampled according to ν: here, YT2 = V5 and a jump towards YT2 occurs at time T2 .

3.1. Diffusion processes on quantum graphs
We adopt here the convention introduced in [24] and fix for any edge e ∈ E of length Le an
orientation (and parametrization se ). This means that se (0) is one of the extremities of e and
se (Le ) is the other one. By doing so, we have determined an orientation for ΓG .
Dynamical system inside one edge Following the parametrization of each edge, we can define the
second order elliptic operator Le as ∀(x, y, t) ∈ e × V (G) × R+ :
Z
1
(4)
Le f (x, y, t) = −βt ∇x Uy (x) + ∆x f (x, y) + αt [f (x, y 0 ) − f (x, y)]dν(y 0 ),
2
G
which is associated with (2) when x ∈ e: the x-component follows a standard diffusion drifted by
∇x Uy (.) inside the edge e, although the y-component jumps over the nodes of the initial graph G
with a jump distribution ν and a rate αt .
Since the drift term ∇Uy (.) is measurable w.r.t. the Lebesgue measure over e and the secondorder part of the operator is uniformly elliptic, Le uniquely defines (in the weak sense) a diffusion
process up to the first time the process hits one of the extremities of e (see, e.g., [33]), which leads
to a Feller Markov semi-group.
Dynamical system near one node We adopt the notation of [26] and write e ∼ v when a vertex
v ∈ V is an extremity of an edge e ∈ E. For any function f on ΓG , at any point v ∈ V , we can define
the directional derivative of f with respect to an edge e ∼ v according to the parametrization of
e. If nv denotes the number of edges e such that e ∼ v, we then obtain nv directional derivatives
designated as (de f (v)):


 lim f (se (h)) − f (se (0))
if
se (0) = v

h−→0+
h
de f (v) =


f (se (Le + h)) − f (se (Le ))

 lim
if
se (Le ) = v.
−
h
h−→0
It is shown in [26] that general dynamics over quantum graphs depend on a set of positive
coefficients:
(
)
X
v
A := (av , (ae,v )e∼v ) ∈ R1+n
s.t. av +
ae,v > 0 : v ∈ V .
+
e∼v

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

9

There then exists a one-to-one correspondence between A and the set of all possible continuous
Markov Feller processes on ΓG . More precisely, if the global generator Lt is defined as
∀(x, y) ∈ ΓG × V

Lt f (x, y) = Le,t (f )(x, y)

when

while f belongs to the domain:
(
D(L) :=

∀y ∈ G

x ∈ e,

)
f (., y) ∈ Cb∞ (ΓG ) : ∀v ∈ V

X

ae de f (v, y) = 0 ,

(5)

e∼v

then the martingale problem is well-posed (see [26, 22]). In our setting, the gluing conditions are
defined through the following set of coefficients A:
∀v ∈ V

av = 0

and

∀e ∼ v

ae,v =

1
.
nv

The gluing conditions defined in (5) induce the following dynamics: when the x component of
(Xt , Yt )t≥0 hits an extremity v ∈ V of an edge e, it is instantaneously reflected in one of the nv edges
connected to v (with a uniform probability distribution over the connected edges) while spending
no time on v. Using the uniform ellipticity of Le and the measurability of the drift term, Theorem
2.1 of [26] can be adapted, providing the well-posedness of the Martingale problem associated with
(L, D(L)), and the next preliminary result can then be obtained.
Theorem 2. The operator L associated with the gluing conditions A generates a Feller Markov
process on ΓG × V × R+ , with continuous sample paths on the x component. This process is weakly
unique and follows the S.D.E. (2) on each e ∈ E.
For simplicity, ∆x will refer to the Laplacian operator with respect to the x-coordinate on the
quantum graph ΓG , using our gluing conditions (5) and the formalism introduced in [26].
3.2.

Convergence of the homogenized S.A. over ΓG

As mentioned earlier, we use a homogenization technique that involves an auxiliary sequence of
random variables (Yn )n≥1 , which are distributed according to ν. More specifically, the stochastic
process (Xt , Yt )y≥0 described above is depicted by its inhomogeneous Markov generator, which can
be split into three parts:
Lt f (x, y) = L1,t f (x, y) + L2,t f (x, y).
In the equality above, L1,t is the part of the generator that acts on Yt :
Z
L1,t f (x, y) = αt [f (x, y 0 ) − f (x, y)]ν(dy 0 ),

(6)

describing the arrival of a new observation Y ∼ ν with a rate αt at time t. Concerning the action
on the x component, the generator is:
1
L2,t f (x, y) = 4x f (x, y) − βt < Ox Uy , Ox f (x, y) > .
(7)
2
Since the couple (Xt , Yt )t≥0 is Markov with a renewal of Y with ν, it can be immediately observed
that the y component is distributed at any time according to ν. We introduce the notation mt
to refer to the distribution of the couple (Xt , Yt ) at time t, and we define nt as the marginal
distribution of Xt . In the following, we will also need to deal with the conditional distribution of

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

10

Yt given the position Xt in ΓG . We will refer to this probability distribution as mt (y |x). To sum
up, we have:
Z
L (Xt , Yt ) = mt with nt (x)dx =
mt (x, y)dy and mt (y |x) := P[Y = y |Xt = x].
(8)
V

A traditional method for establishing the convergence of S.A. towards the minimum of a function
Uν consists in studying the evolution of the law of (Xt )t≥0 and, in particular, its close relationship
with the Gibbs field µβt with energy Uν and inverse temperature βt :
1
exp(−βt Uν (x)),
Zβt
Z
is the normalization factor, i.e., Zβt = e−βt Uν (x) dx.
µβt =

where Zβt

(9)

Using a slowly decreasing temperature scheme t 7−→ βt−1 , it is expected that the process (Xt )t≥0
evolves sufficiently fast (over the state space) in the ergodic sense so that its law remains close to µβt .
In addition, the Laplace method on the sequence (µβt )t≥0 ensures that the measure is concentrated
near the global minimum of Uν (see, for example, the large deviation principle associated with
(µβ )β→+∞ in [25]).
Hence, a natural consequence of the convergence “L (Xt ) −→ µβt ” and of the weak asymptotic
concentration of (µβt )t≥0 around Mν would be the almost sure convergence of the algorithm towards
the Fréchet mean:
lim P (Xt ∈ Mν ) = 1.

t−→+∞

We refer to [32] for further details. In particular, a strong requirement for this convergence can be
considered through the relative entropy of nt (the law of Xt ) with respect to µβt :


Z
nt (x)
dnt (x)
(10)
Jt := KL(nt ||µβt ) =
log
µβt (x)
ΓG
The function βt will be chosen as a C 1 function of R+ , and since µβ is a strictly positive measure over
ΓG , it implies that t 7−→ µβt (x)−1 is C 1 (R+ × ΓG ). Moreover, (t, x) 7−→ nt (x) follows the backward
Kolmogorov equation, which induces a C 1 (R+ × ΓG ) function. Since the semi-group is uniformly
elliptic on the x-component, we have:
∀t > 0

∀x ∈ ΓG

nt (x) > 0.

On the basis of these arguments, we can deduce:
Proposition 1. Assume that t 7−→ βt is C 1 , then (t, x) 7−→ nt (x) defines a positive C 1 (R+ × ΓG )
function and t 7−→ Jt is differentiable for any t > 0.
If we define:
αt = λ(t + 1) and βt = b log(t + 1),
with b a constant smaller than c? (Uν )−1 , where c? (Uν ) is the maximal depth of a well containing a
local but not global minimum of Uν , defined in (34), then our main result can be stated as follows:
Theorem 3. For any constant λ > 0 such that αt = λ(t + 1) and βt = b log(1 + t) with b >
{c (Uν )−1 }−1 , then:
lim Jt = 0
as
t −→ +∞.
?

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

11

This ensures that the process Xt will almost surely converge towards Mν and, therefore, towards
a global minimum of Uν .
The idea of the proof is to obtain a differential inequality for Jt , which implies its convergence
towards 0. It is well known that the Gibbs measure µβt is the unique invariant distribution of the
stochastic process that evolves only on the x component, whose Markov generator is given by:
1
L̂2,t (f )(x) = ∆x f − βt h∇x f, ∇x Uν i.
2

(11)

Therefore, a natural step of the proof will be to control the difference between L2,t and L̂2,t and
to use this difference to study the evolution of mt and nt . It can be seen that L̂2,t may be written
as an average action of the operator thanks to the linearity of the gradient operator:
1
L̂2,t (f ) = ∆x f − βt Ey∼ν h∇x f, ∇x Uy i.
2
When Xt = x, we know that Yt is distributed according to the distribution mt (y |x). Consequently,
the average action of L2,t on the x component is:
Z
1
1
L̃2,t (f ) = ∆x f − βt Ey∼mt (.|x) h∇x f, ∇x Uy i = ∆x f − βt h∇x f, ∇x Uy imt (y |x)dy,
(12)
2
2
V
whose expression may be close to that of L̂2,t if mt (.|x) is close to ν.
Thus, another important step is to choose appropriate values for αt and βt , i.e., to find the
balance between the increasing intensity of the Poisson process and the decreasing temperature
schedule, in order to quantify the distance between ν and mt (.|x). The main core of the proof
brings together these two aspects and is detailed in Section 5.
Another important step will be the use of functional inequalities (Poincaré and log-Sobolev
inequalities) over ΓG for the measure µβ when β = 0 and β −→ +∞. The proof of these technical
results are given in Section 6.3.
Corollary 1. Assume that βt = b log(t + 1) with b < c? (Uν )−1 and αt = λ(t + 1)γ with γ ≥ 1,
then for any neighborhood N of Mν :
lim P[Xt ∈ N ] = 1.

t−→+∞

Proof: The argument follows from Theorem 3. Consider any neighborhood N of Mν . The continuity of Uν shows that:
∃δ > 0

N c ⊂ {x ∈ ΓG : U (x) > min U + δ }.

Hence,
P[Xt ∈ N c ] ≤ P[U
Z (Xt ) > min U + δ]
=
1U (x)>min U +δ dnt (x)
ZΓG
Z
=
1U (x)>min U +δ dµβt (x) +
ΓG

1U (x)>min U +δ [nt (x) − µβt (x)]dx

ΓG

≤ µβt {U > min U + δ } + p
2dT V (nt , µβt )
≤ µβt {U > min U + δ } + 2Jt ,

where we used the variational formulation of the total variation distance and the CsiszárKullback inequality. As soon as limt−→+∞ Jt = 0, we can conclude the proof observing that
µβt {U > min U + δ } −→ 0 as βt −→ +∞.


Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

12

It can actually be proven (not shown here) that limt−→+∞ P[Xt ∈ N ] = 1 if and only if the
constant b is chosen to be lower than c? (Uν )−1 . This means that this algorithm does not allow
faster cooling schedules than the classical S.A. algorithm. Nevertheless, this is a positive result
since this homogenized S.A. can be numerically computed quickly and easily on large graphs.
Finally, we should consider this result to be theoretical. However, in practice, the simulation of this
homogenized S.A. is performed during a finite horizon time and efficient implementations certainly
deserve a specific theoretical study following the works of [17] and [46].

4.

Numerical simulation results

We now present different practical aspects of our strategy. After giving numerical details about
how to make it usable on large graphs, plus different insights into its calibration, we will present
results obtained on social network subgraphs of Facebook 2 and a citation subgraph of zbMATH.
4.1. Algorithmic complexity
It is widely recognized that the algorithmic complexity of numerical strategies dealing with graphs
can be an issue. The number of vertices and edges of real-life graphs can indeed be quite large.
For instance, the zbMATH subgraph of Section 4.3.2 has 13000 nodes and approximately 48000
edges. In this context, it is worth justifying that the motion of (Xt )t≥0 on the graph ΓG across the
iterations of Algorithm 3 is reasonably demanding in terms of computational resources. We recall
that this motion is driven by a Brownian motion (lines 4 to 6 of Algorithm 3) and an attraction
towards the vertex YNtα (line 9 of Algorithm 3) at random times (Tk )k≥0 . The number of vertices
and undirected edges with non-null weights in ΓG is also N and |E |, respectively. Note finally that
N − 1 ≤ |E | ≤ N (N − 1)/2 if ΓG has a unique connected component.
Neighborhood structure A first computational issue that can arise when moving Xt is to find all
possible neighbors of a specific vertex v. This is indeed performed every time Xt moves from one
edge to another and directly depends on how ΓG is encoded in the memory. Encoding ΓG in a list
of edges with non-null weights is common practice. In that case, the computer checks the vertex
pairs linked by all edges to find those containing v, so the algorithmic cost is 2|E |. A second classic
strategy is to encode the graph in a connectivity matrix. In this case, the computer has to go
through all the N indexes of the columns representing v to find the non-null weights. We instead
sparsely encode the graph in a list of lists: the main list has a size N and each of its elements lists
the neighbors of a specific vertex v. If the graph only contains edges with strictly positive weights,
this strategy has a computational cost N , and in all of the other cases, the cost is < N . On average,
the computational cost is 2|E |/N so that this strategy is particularly advantageous for sparse
graphs, where |E | << N (N − 1)/2, which are common for the targeted applications. For instance,
|E | = 4 ∗ 103 and N (N − 1)/2 = 1.24 ∗ 105 on the smallest Facebook subgraph of Section 4.3.1, and
|E | = 4.8 ∗ 104 and N (N − 1)/2 = 8.44 ∗ 107 on the zbMATH subgraph of Section 4.3.2.
Geodesic paths Another potential issue with our strategy is that it seeks at least one optimal
path between the vertices Xt and YNtα at each jump time (line 9 of Algorithm 3). This may be
efficiently done using a fast marching propagation algorithm (see, e.g., [18] for details), where:
(i) the distance to Xt is iteratively propagated on the whole graph ΓG until no more optimal
distance to reach Xt is updated, and
(ii) considering the shortest path between Xt and YNtα .
In this case, step (i) is particularly time-consuming and cannot be reasonably performed at each
step of Algorithm 3. Fortunately, the algorithmic cost to compute the distance between all pairs of
vertices is equivalent to that of computing step (i). We then compute these distances once and for
2

Subgraphs obtained on the Stanford Large Network Dataset Collection: https://snap.stanford.edu/data/

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

13

all at the beginning of the computations and store the result in a N × N matrix in the RAM of
the computer 3 . We can therefore very quickly use these results at each iteration of the algorithm
and, in particular, deduce the useful part of the geodesic paths involved to travel from Xt to YNtα .
The only limitation of this strategy is that the distance matrix can be memory-consuming. It can
therefore be used on small to large graphs but not on huge graphs (typically when N > 105 ) on
current desktops. An extension to deal with this scalability issue is a current subject of research,
and we will therefore not describe applications of Algorithm 3 on huge graphs in this work.

4.2. Parameter tuning
Several parameters influence the behavior of the simulated process Xt . Some of them are directly
introduced in the theoretical construction of the algorithm, i.e., the intensity of the Poisson process
or the temperature schedule. Other ones come from the practical implementation of the algorithm,
i.e., the maximal time up to which we generate Xt . The theoretical result given in Theorem 3 gives
an upper bound for the probability of Xt to be not too distant from the set of global minima. This
bound depends on βt and αt as well as on different characteristics of the graph such as its diameter
and number of nodes. We then propose an empirical strategy for parameter tuning, which we will
use later in Section 4.3.
For a graph G with N nodes and a diameter DG , we define:
∗
Tmax
= 100 + 0.1N

and

βt∗ = 2 log(t + 1)/DG .

We then choose the intensity of the Poisson process so that it will generate a reasonably large
amount of Yn at the end of the algorithm, depending on the discrete probability distribution ν.
∗
∗
. We
− 1 and Tmax
More specifically, let S ∗ be the average number of jumping times between Tmax
∗
then set S = 1000, which can be obtained by choosing:
αt∗ = λ(2t + 2)

with

λ=

S∗
∗
2Tmax
+1

.

4.3. Results
This section presents results obtained on graphs of different sizes and using different parameters.
We used the empirical method of Section 4.2 to define default parameters and altered them to
quantify the sensitivity of our algorithm to parameter variations.
Since the process Xt lives on a quantum graph, its location at time t is between two vertices, on
an edge of ΓG . However, our primary interest is to study the properties of the initial discrete graph
and therefore, the output of the algorithm will be the vertex considered as the graph barycenter.
To achieve this, we associate a frequency to each vertex. For a vertex v, this frequency is the
portion of time during which v was the closest vertex to the simulated process Xt . In our results,
we consider frequencies computed over the last 10% of the iterations of the algorithm.
4.3.1.
3

Facebook subgraphs

In our experiments, we used the all-shortest-paths function of the Python library NetworkX.

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

14

Experimental protocol We first tested our algorithm on three subgraphs of Facebook from the
Stanford Large Network Dataset Collection: (FB500) has 500 nodes and 4337 edges and contains
two obvious clusters; (FB2000) has 2000 nodes and 37645 edges and fully contains (FB500);
and (FB4000) has 4039 nodes and 88234 edges and fully contains (FB2000). For each of these
subgraphs, we considered the probability measure ν as the uniform distribution over the graph’s
vertices and used a length of 1 for all edges. We also explicitly computed the barycenter of these
graphs using an exhaustive search procedure. We are able to do that because we considered a
simplified case where the distribution over the nodes is uniform. For example, this exhaustive search
procedure required approximately 6 hours for the (FB4000) subgraph. Nevertheless, it allowed us
to compare our strategy with ground-truth results.
We used the strategy of Section 4.2 to define default parameters adapted to each subgraph. We
also tested different values for parameters β, S and Tmax in order to quantify their influence. We
repeated our algorithm 100 times for each parameter set to evaluate the algorithm stability.
In the tables representing quantitative results, Error represents the number of times, out of
100, that the algorithm converged to a node different from the ground-truth barycenter. It is a
rough indicator of the ability of the algorithm to locate the barycenter of the graph, that could be
replaced by a measure of the average distance between the last iterations of the algorithm and the
ground-truth barycenter (not shown in this work). For each subgraph and parameter set, column
Av. time contains the average times in seconds for the barycenter estimation, keeping in mind that
the Dijkstra algorithm was performed once for all before the 100 estimations. This preliminary
computation requires approximately 1, 30 and 80 seconds on a standard laptop.
Effect of β From a theoretical point of view, (βt )t≥0 should be chosen in relation to the constant
c? (Uν ), which is unknown in practice. Therefore, the practical choice of (βt )t≥0 is a real issue to
obtain a good behavior of the algorithm. Table 1 gives the results obtained on our algorithm with
different values of βt .

β
1 ∗
β
4
1 ∗
β
2

β∗
2β ∗
4β ∗
8β ∗
Table 1

FB500
Error Med. Freq. Av. time
15 %
0.6042
0.95 s
2%
0.4184
1.38 s
0%
0.8008
1.48 s
0%
0.8321
11.31 s
0%
0.8233
13.19 s
0%
0.7717
2.36 s

FB2000
Error Med. Freq. Av. time
28 %
0.3519
10.12 s
1%
0.7268
4.25 s
0%
0.9418
12.36 s
0%
0.9892
8.52 s
0%
0.9930
15.58 s
0%
0.9750
23.25 s

FB4000
Error Med. Freq. Av. time
27 %
0.3314
12.41 s
5%
0.6534
14.11 s
1%
0.8913
13.55 s
0%
0.9647
16.79 s
0%
0.9824
31.43 s
0%
0.9445
44.46 s

Experiments on the three Facebook subgraphs, while varying the values of (βt )t≥0 .

We can observe that when (βt )t≥0 is too small, then the behavior of the algorithm is deteriorated,
∗
revealing the tendancy of the process (Xt )0≤t≤Tmax
to have an excessively slow convergence rate
towards its local attractor in the graph. Roughly speaking, in such a situation, the process does
not learn fast enough. When the value of β is chosen in the range [β ∗ ; 2β ∗ ], we can observe a really
good behavior of the algorithm: it almost always locates the good barycenter in a quite reasonable
time of computation (less then 20 seconds for the largest graph). Finally, we can observe in the
column, Med. Freq., that in most of the last iterations of the algorithm (more than 80%), the
process evolves around its estimated barycenter, so that the decision to produce an estimator is
quite easy when looking at an execution of the algorithm.

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

15

Effect of S and Tmax Table 2 gives the results obtained with our algorithm while using different
values of S and Tmax . As expected, we observe that increasing the ending time of simulation always
improves the convergence rate (column Error in Table 2) of the algorithm towards the right node.
The behavior of the algorithm is also improved by increasing the value of S, which quantifies the
number of arrivals of nodes observed along the averaging procedure. Of course, the counterpart of
increasing both S and Tmax is an increasing cost of simulation (see column Av. time).

S
1 ∗
S
2

S∗
S∗
2S ∗
2S ∗
2S ∗

Tmax
∗
Tmax
∗
2Tmax
∗
4Tmax
∗
Tmax
∗
2Tmax
∗
4Tmax

Table 2

FB500
Error Med. Freq. Av. time
2%
0.7667
0.59 s
0%
0.8049
2.79 s
0%
0.8101
5.60 s
1%
0.8345
2.26 s
0%
0.8361
4.57 s
0%
0.8423
11.16 s

FB2000
Error Med. Freq. Av. time
0%
0.9344
3.22 s
0%
0.9610
9.30 s
0%
0.9677
22.21 s
0%
0.9512
8.39 s
0%
0.9586
23.30 s
0%
0.9735
11.16 s

FB4000
Error Med. Freq. Av. time
0%
0.8614
6.70 s
0%
0.8970
24.42 s
0%
0.9222
54.73 s
0%
0.9062
20.95 s
0%
0.9121
23.30 s
0%
0.9366
96.84 s

Experiments on the three Facebook subgraphs, while varying the values of S and Tmax .

As an illustration of the (small) complexity of the Facebook sub-graphs used for benchmarking
our algorithm, we provide a representation of the FB500 graph in Figure 3. This representation
has been obtained with the help of Cytoscape software and is not a result of our own algorithm.
In Figure 3, the red node is the estimated barycenter, which is also the ground-truth barycenter
located by a direct exhaustive computation. The blue nodes are the “second rank” nodes visited
by our method.
N468

N309

N356

N36

N135

N411

N443

N139

N183

N453

N206

N499

N467

N209

N447

N228
N108

N251

N197

N454
N160
N127

N336
N380
N153

N164

N485

N284

N459

N74

N472

N281

N184

N498

N320

N76

N159

N358

N117

N300
N270

N351

N424

N435

N126

N294
N178

N283

N437

N191
N384

N406
N215
N221

N486

N73

N121

N477

N27

N357
N90
N105

N10

N401

N57
N457
N45

N345
N476

N165
N88
N364

N30

N303

N260
N389
N54
N381

N50

N77

N299
N346

N39

N385

N24
N449

N69
N179

N257
N371

N480
N272

N130

N133

N3

N350

N106

N280

N236
N292
N369
N297
N231

N1
N302
N92
N362

N66
N433
N38
N94

N399

N210
N469
N441

N386

N141
N266

N329

N448

N142
N79

N21

N331
N72

N332

N420

N63

N471

N48
N415
N80
N318

N413
N495
N425

N256

N194

N198
N145

N75

N161

N204

N473
N123

N53

N276
N113

N200
N186

N122

N464

N285
N104
N254

N426

N419
N430

N16
N439

N134

N397

N462

N451
N423

N208

N322
N330

N461
N213
N180

N26

N169

N67
N393
N249

N9
N323
N188
N342
N242

N407
N392

N233
N440

N348

N341
N118
N202

N55

N445

N15

N460

N466
N442

N187
N344

N387

N379

N479
N315

N56

N203

N421

N277
N395
N166
N25

N239

N234
N148
N390
N370
N85
N290
N101
N261
N196
N404

N173

N62

N295

N431

N463

N444

N427

N455
N465

N298
N119

N98
N271

N109
N252

N363

N37

N409
N432

N456

N484

N5
N258
N172
N334
N304

N373
N274
N372

N316

N29
N335

N428
N291

N402

N238

N422

N40

N232

N374

N224

N352
N429
N146

N34

N185
N18
N199

N493
N483
N107
N414

N154
N82

N481
N100

N492

N176

N360

N207

N497

N170
N212
N163

N314
N458

N353

N354

N376
N313
N248
N156

N412
N268

N325

N132

N398

N308

N96
N408

N12
N400

N265

N452

N368
N324

N87

N450
N436

N171

N378

N375

N64

N158

N11

N355
N339

N0

N13

N128
N59

N35

N211
N223

N222

N150

N475

N217

N470

N58

N438

N403

N417

N366

N136
N487

N247
N418

N317

N494

N103
N7
N489

N31

N388

N482

N42
N235
N65

N269

N391
N219

N361

N60

N496

N338

N416

N311
N246

N189

N6

N367

N394

N33

N359

N84

N321

N434

N396
N382

N340

N168
N488
N288

N263

N250

N327

N490

N229

N22

N296
N147
N4

N195

N474

N120
N446

N410

N152
N125

N347

N491

N319

N78

N383

N143

N275

N129

N175

N405
N83

N237

N68

N95

N306
N51

N81

N365
N190
N478

N89
N102

N255

N240

N218
N181

N99

N124

N278

N349

N301

N46

N307

N61

N225
N273
N287

N328
N227

N241

N286

N23

N49
N47

N177
N377

N43

N155

N289
N131

N157
N193

N8

N267

N19
N264

N71
N86

N245

N17

N116

N140

N174

N114

N162
N91

N253

N230

N182
N115

N262

N192

N110

N293

N144

N214

N149

N259

N138

N97

N326

N137

N2

N201
N70
N52

N243

N41

N220
N312

N343
N112

N28

N205
N310
N244

N93
N167

N333

N20

N226

N111

N14

N282

N337

N216

N151

N32

N305

N44

N279

Figure 3

Region of interest presenting the results obtained on the FB500 graph.

4.3.2. zbMATH subgraph The zbMATH subgraph built from zbMATH4 has been obtained
by an iterative exploration of the co-authorship relationship in an alphabetical order. We thus naturally obtained a connected graph. This graph has been weighted by the complete number of citations obtained by each author, leading to the distribution probability ν ∝ ](number of citations).
Finally, all the edges in the graph are fixed to have a length of 1. This exploration was initialized on
the entry of the first author’s name (i.e., S. Gadat) and we stopped the process when we obtained
13000 nodes (authors) on the graph. This stopping criterion in the exploration of the zbMATH
database corresponds to a technical limitation of 40 GB memory required by the distance matrix
4

https://zbmath.org/authors/

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

16

obtained with the Dijsktra algorithm. In particular, this limitation and the starting point of the
exploration induce an important bias in the community of authors used to build the subgraph from
the zbMATH dataset: the researchers obtained in the subgraph are generally French and applied
mathematicians. The graph is more or less focused on the following research themes: probability, statistics and P.D.E. Consequently, the results provided below should be understood as an
illustration of our algorithm and not as a bibliometric study!
∗
Our experiments rely on the same choice of parameters (βt∗ )t≥0 , S ∗ , Tmax
and αt∗ as above,
indicated at the beginning of Section 4.2. Again, our algorithm produces the same outputs but of
course, in this situation, we do not know the ground-truth barycenter of the zbMath subgraph.
In Figure 4, we present a representation of the subgraph obtained with Cytoscape software, and a
zoom on a region of interest (ROI for short) in Figure 4.
Again, the main nodes visited by our algorithm are represented with a red square and the size
of the used square is larger when the node is frequently visited. According to the results obtained
on the Facebook subgraphs, we then assume that the larger red square is the barycenter of the
zbMath subgraph.
Shiryaev
Khelemski\xc4\xad

Godunov

Z\xc3\xa4hle

Ivanoff

Hawkes

Wang

Kallsen

Sudakov
St\xc3\xb3s

Klenke

Bassalygo

Onishchik

Kozlov

Beck

Goodman

Pitman

Ghomrasni
Korolyuk

Kotani
Kwa\\xc5\\x9bnicki

Raugi

Tanemura

Marckert

Pinsker

Dobrushin
Pietruska-Pa\xc5\x82uba

Popescu

Goldschmidt

Feldman
Chavel

Menard

Aleksandrov

Drapeau
Issoglio
Mekler

Ershov

Haas

Levin

L\\xc3\\xb3pez-Mimbela

Baumgarten

Sevast\xca\xb9yanov

Rosen

Utzet

Coquet

Bonichon

Kendall

Windridge

Bakhtin

Profeta

Tang

Horowitz

\\xc3\\x9cst\\xc3\\xbcnel
Melnikov

Hubalek

S\xc5\x82omi\xc5\x84ski

Schoenmakers

Eckhoff

RuizdeCh\\xc3\\xa1vez

Iscoe

vanDijk

Barrieu

Keevash

Miermont

Szulga

Norros

Rudiger

Mallows

Simatos

Zindy

Hu
Hughston

Lavrent\xca\xb9ev

Demichel

Gao

Prelov

Smirnov

Lubachevsky

Bolthausen

Fayolle

Reshetnyak
Konno

Sztonyk

Reichmann

Spehner

Berestycki

Dobri\xc4\x87

Kogan

Tribe

Zinn

Karpelevich

Aistleitner

Hill

Chern

OConnell
Gurevich
Raja

Zalgaller

Kutateladze

Robert

M\\xc3\\xa9min

Curien
Kingman

Yang

Bucklew

Doebner

Zabrodin

Siri-Jegousse

Foondun

Pecherski\xc4\xad

Andres

Kane

Warren

Taylor

Shields

Chevallier

Fujita

Hattori

Madiman

Meckes

Kabanov

Khas\\xca\\xb9minski\\xc4\\xad

Freedman

Filipovi\xc4\x87

Ney

Assaf

Biagini

Klein

Bronstein

Rhoades

Martin

Carmona

Dani

Birkner
Schweinsberg

Maslova

Marton

Jung

Sina\xc4\xad

Sayit

Stoimenov
NavarroCruz

Walsh

Feng

Zili

Durrett

Mayer-Wolf

Winter

Pra\\xc5\\x82at

Madan

Nikeghbali

Weisz

Evans

Biane

Minlos

DAristotile

Chigansky

Friedli

Rolles

Shlosman

Kulczycki

Yor

Mytnik
Morters

Pirogov

Kifer

Ayache
Bunimovich
Jakobson

Xu

Sukhov
Tsirelson

Bol\xca\xb9shev

LeGall

Privault

Lawler

Fitzsimmons

Vershik

Varadarajan

Budhiraja

Delbaen
Hinz
Hillairet
Weil
K\\xc3\\xb6rner

Panholzer
Luschgy
Reitzner

Jonasson

Miclo

Aaronson

Akhavi

Diaconis

Luczak

Berend

Taniguchi

B\\xc3\\xa1r\\xc3\\xa1ny

Pemantle

Kuba

Limic

Arguin

Cotar

Dawson

Ortgiese

Bass

Pimentel

Song

Swanson

Burdzy

Muhle-Karbe

LeNy

Beffara

Paycha

Gun

Janson

Alonso

Bordenave

Shinkar
Jolis
Gurel-Gurevich Mart\\xc3\\xadnez

Tudor

Hu

Benjamini

Zeitouni

Guiol

Newman

Garc\\xc3\\xada
Tibi
Schneider

Marsalle

Mairesse

Verd\xc3\xba

Bollob\\xc3\\xa1s
Ruci\xc5\x84ski

Baldi
Es-Sebaiy

Valencia-Pabon

Szpirglas

Ichiba

McDiarmid

Alexander

BenArous

Mahfoudh

Soundararajan

Ueltschi

Moore

Yuan

Gombani

Yadin

Rifo

VivesSanta-Eul\\xc3\\xa0ia

Coutin
Dhersin

Kozma

Stricker
Follmer
Crisan

Quer-Sardanyons

Pontier

Berger

Hoffman

Greven

Xu

Amami

Yamazaki

Cohen

Ferrari

Forde

Frohlich

Promislow

Rosenkrantz

Kira

Freire

Lacayo

Keane
Grishin

Zimmermann
Waterman

Romik

Grimmett

Cairoli

Tudor

Koralov

Ahn

Vormoor
Ross
Ouyang

Bardina

Tran

Henard

Moyal
Decreusefond

Bjork

Kurtz

Yuan

Spruill

Zajic

Alencar

Garz\xc3\xb3n

Fradon
Ruszel

Kree

Gardy
Milevsky

Puech

Stanley

Chassaing

Spencer

Panloup

Maillard

Schonmann

Maes

Guionnet

Kusuoka

Pene

Winkler
Tsur

Redig

Fernandez

Bovier

Yan
Gjerde

Smolyanov

Ramanan

Schramm

Chronopoulou
Duminil-Copin

Krylov
Salopek

Cranston

Neveu
Marchetti

Ladyzhenskaya

Wang

Brossard

Dozzi

Wilson

Angel

Liggett

Mountford

Kuptsov

Andjel

Reingold

Raimond

Phan

Gatheral

Hug

Virag

Nualart

Sidoravicius

Sortais

Zambrini

Vondra\xc4\x8dek

Zhang
Bernard

Corteel
K\\xc3\\xb6rezlio\\xc7\\xa7lu

Kurtz

Peled

denHollander
T\xc3\xb3th

Gantert
Ravishankar

Bahadoran

Ramirez
Remenik

Ub\\xc3\\xb8e

Wang

McKean

Vall\\xc3\\xa9e
Z\\xc3\\xa9mor

Szpankowski

Gillet

Noreddine

Guillotin-Plantard

Fritz

Molchanov
Huang

Auffinger

Neate

Oudjane
Ludkovski

G\xc3\xb3ra

Brightwell
Hammond
Ma\xc3\xafda

Bourgade

Xu

Quastel

Richou

Remy

Sly

Lyons

Nourdin

G\xc3\xa4rtner

Corwin

Najnudel

Hu
Unterberger

Scheicher

Deguy

Anulova

Stauffer

Peres
March
\xc4\x8cern\xc3\xbd

Revuz

Xu

Maisonneuve

Yang

Butko

Schott

Izkovsky

LeJan

Ye

Denise

Kardaras

Rolla

Caravenna

VanMoffaert

Kendall
Suryanarayana

Bottcher

Holley

Wu

H\\xc3\\xa4ggstr\\xc3\\xb6m

Jarai

Kumagai

Jacob

Kuwae

Franchi

Rubenthaler

Paccaut

Steif

Tudor
Cabezas

Kusuoka

Flajolet

Francon
Tolley

Gouyou-Beauchamps

deLima

dosSantos

Royfman
Takeda

Derguzov

Barbaresco

Ying

Bousquet-M\\xc3\\xa9lou
Thimonier

Herbertsson

Graham

Bogachev

Krawczyk

Sidorova

Fribergh

Marstrand

Fehr

Sol\\xc3\\xa9

Poisat

Ma

Cipriani

Truman

Caruana

Veretennikov

Mania

Dalang

Nardi

Deuschel

Xiang

Zhou

Krasovski\xc4\xad

L\\xc3\\xa9vy

Banderier

vanderHofstad

He

Hryniv

Housworth

Pouyanne
Dorea

Steele

Dombry

Bayer

Ponty

Rouault
Louchard

\xc5\x81uczak
Crawford

Schinazi
Fleischmann

Chen

Sagna

Sabot

Nesterenko

Randrianarimanana

Baxendale

Kesten
vandenBerg

Diehl

Ravelomanana

Yukich
Geniet

Wakolbinger

Goodman

Chen

dePanafieu

Ravanelli

Kulske

Merzbach

Poly

Heyne

Lunt

Regnier

Shwartz

Ofek

Revelle

Cherny

Barlow

Monch

Aase

Uemura

Liu

Zhang
Weiss

Vvedenskaya

Gerin

Ouerdiane

Saussol

Bansaye

Veber

Berestycki

Prodinger

Ba\xc4\xaddak

Mendez-Hernandez

Ma

Engl\\xc3\\xa4nder

Gittenberger

Trashorras

Pal

Henkel

Young

Antunes

Chung

Enriquez

Drmota

Simonyi

Kloppel

Wang

Zhou

Montseny

Gorostiza

Yan

Hsu

Lavault

Bensusan

Th\xc3\xa4le

Shkolnikov

Sun

Song

Feinberg

Bahsoun

Konig

Shen

Billera

Aldous

Ewald

Dembo
Williams

Xiong
Brydges

Davis

Mahmoud

Lew

Huang

Kupper

Kornfeld

Croydon

Reuter

Hordijk

Sen
Doring

Hattori

Khanin

Lelarge

Estrade

Donati-Martin
Port\\xc3\\xa8s

Kolmogorov

Smillie

Pesin
Bo

Kim

Cohen

Wilbertz

Yang

Kirillov

Benaych-Georges

Temme

Salhi

Brown

Gallardo

Cox
L\\xc3\\xa9ger

Rozanov

Winter

Williams

Gong

El-Nouty

Adelman

Oseledets

Hambly

Bogdan

Durran

Valchev

Balan

Zani

Braverman

Park

Barth
Santacroce
Schweizer
Cao
Singh

Koteck\\xc3\\xbd

Durand

Garcia

Vares

Roelly
Rovira

Nualart

Jansen

Zhang

Dimitroff

Logothetis

Schmitt

Lindstr\xc3\xb8m

\xc3\x87etin

Meyer

Albeverio

Lambert

Deya

Cadel

R\\xc3\\xa9veillac

Lyons

Delattre

Bugeaud

Gayrard

P\\xc3\\xa9ch\\xc3\\xa9

Pettersson

Al\\xc3\\xb2s

Yang

Vallois

Abadi

Tindel

Viens

Merola

Martinelli

Hyland

Quas

Sznitman

Stauch

Ledoux

Litterer

Zheng
Laurence

Razborov

Cass

Liu

Falconnet

Loukianova

Darses
Grorud
Bayraktar

ODonnell

Celani
B\\xc3\\xb6r\\xc3\\xb6czky

Merhav

Zhang

\xc5\xbdigo

Kannan

Runggaldier

Chesney
Kipnis

Herrmann

Carmona
Taylor

Yao

Lamberton

Fontana

Roussignol
Ouahhabi

Muller

Fernandez

Chaumont

Murr

Imkeller

Gloter

Brassesco

Ianiro

Chalendar

Ankirchner

Riedle

Castell

Hu
Cruzeiro

Najim

Locherbach

Karatzas

Nadtochiy

Liang
Khorunzhiy

Hanlon

Gu\\xc5\\xa3\\xc4\\x83

Villemonais
Jona-Lasinio

Bouleau

Jeanblanc

Varadhan
Scheutzow

Nguy\xe1\xbb\x85nHuuDu

Kleinberg

Makarychev

Guiol

Schilling
Guasoni

Presutti

Kruse

vonWeizs\\xc3\\xa4cker

Zhang

Creutzig

Peithmann
Deng

Wetzel

Zervos

Briand
Cassandro

Durrleman

Klein
Mohammed

Hein

Lorentz

Dupuis
Pokern

Roynette
Triolo

Gassiat

Davies

Bardou

Spencer

Teplyaev

Lima

Cancrini

Landim

Choulli

Parry
Frikha

Eremenko

Luccio

Olla

Spohn
Plank
Ulsamer

Straka
Frieze

DosReis

Qian

Mazliak

Grandits

Erraoui

Hennequin
Bouthemy

Leung

Zhang

Weiss

Baudoin

Aleksandrov

Delyon

Loukianov

Bertoglio

Cesi

Rozovskii

Werner

Yang
Marcuard

Hayman

Driver

Lepingle

Yang
Zhou

Mathieu

Benois

Kohatsu-Higa

ElKaroui

SanMart\\xc3\\xadn

Davis

Lacoste

Joulin
Toninelli

Mourragui

Bouten

Faggionato

Janssens

Shamir

Bonnefont

Victoir

Pavlyukevich

Blondel

Mao

Semenov

Rudolph

M\\xc3\\xa9l\\xc3\\xa9ard
Suen
Dorobantu

Stoica

Bakry

Viswanathan

Eddahbi

Ouknine
Michel

Gabrielli

DeMasi
Bertini Sellami

Fehringer
Sab

Zambotti

Dupoiron

Petunin

Evstigneev
Bally

Freidlin

Karp

Appleby

Takahashi

Denis

Bogachev

AitOuahra

Scarlatti
Swords

Kumar

Huang

deSaporta

Lakhel

Callegaro

Gundy

Thalmaier

Benassi

Aza\xc3\xafs

vanGaans
Kohlmann

Rockner

Grigoryan

Ancona

Kravitz

Rulliere

Guerbaz

Pardoux

Dereich

Aksamit

Schottstedt
Kruk

Xu

Chen

Fischer

Protter
Conforti

Sarr\xc3\xa0
Berkaoui

Elworthy
Batchelor

Lin

Horst

Gr\\xc4\\x83dinaru
Whittle

Petrov

Qian

Alon

Torres

Gobron

Bressaud

Zhang

Brown

Fang
Ledrappier

Chen

Geman

Tarr\\xc3\\xa8s

Torrecilla-Tarantino

Brunel

Shao

Galamba

LeGland

Sazonov

Ta\\xc3\\xafbi

Anderson

Friz
Ouyang

Tudor

Xu

Delmas

Tetali
Caputo

Collet

Sanz-Sole

Schneider

Richter

Allouche

Wang

ROI

Duarte

Lubetzky

Porchet
M\\xc3\\xa1rquez-Carreras

Olivieri

Dartnell

Maass

Arnaudon
Basdevant

Coulibaly-Pasquier

Bertein

Coculescu

Caderoni

Pag\\xc3\\xa8s
Jiao

Dellacherie

Metz

Duan

Galves

Galves

Hitczenko

Hobson

Zhang

Thieullen

Ob\xc5\x82\xc3\xb3j

Guivarch

\\xc3\\x89mery
Az\\xc3\\xa9ma
Michon
Aubrun
Robinson

Xing

Teichmann

Kramkov
Saada Bezerra

Li
Aida

Leonardi

Yehudayoff

Comets
Fiel

Baik

Deng

Vergne

Mossel
Hernandez-Hernandez

Picco

Mazet
Chen

Mao

Ferraz

Egami

Petz

Song

Tsagkarogiannis

Fromm

Wang

Liptser

Le\\xc3\\xb3n
Scacciatelli

Bao

Sol\\xc3\\xa9
Randriambololona

Frangos

Schapira

Altman

Cocozza-Thivent

Bascompte

Pellegrinotti
Aizenman

Gaines

Wu

Lov\\xc3\\xa1sz

Shapira

Erd\xc5\x91s

Zessin

\\xc3\\x98ksendal

Fukushima

Ren

Lau

Fricker

Chung

Meda

Fontes

Kirkpatrick

Mueller

Nadakuditi

Gouere
Gustafson

Kershaw

Gobet

Blanchet-Scalliet

Butt\xc3\xa0

Feige
Mordecki

Stoltz
Lladser

Edgar

Mathieu

Nazarov

vandenBerg

Pandurangan

Bardet

Zargari

Arnold

Feng
Gairing

Setterqvist

Cattiaux
Matoussi

Marinelli

Lueker

Elliott

Grossi

Azar

Rost

Yao

Duplantier

Krell

Beznea

Stroock

Lepeltier

Zhang

Broder

Nzi

\xc5\x81ochowski

\\xc3\\x89tor\\xc3\\xa9

DaiPra

Szatzschneider

Pietracaprina
Peleg

Mellouk

Slivkins

Zitt

Zhang

Pisanti

DasSarma

Millet

Rodkina

Russo

Ferri\\xc3\\xa8re

Schulman

Baggett

Chafa\\xc3\\xaf

Hajji

Vincenzi
Wigderson

Calka

Molla

Mercier

Robinson

Glover

Sodin

Hogele

Mauldin

Boufoussi

Hairer

Raghavan

Schachermayer

Mityagin

Carey
Sworder

Augustine

Lejay

Lapeyre
Sturm

L\\xc3\\xa9onard

LeBorgne

Cisse

Voiculescu

Masters

Montanari

Sirbu

Tanr\\xc3\\xa9
BenAlaya

Faddeev

Upfal
Karlin
Ruh

Guillin

Orlandi

Zabczyk

Graf

Borodin

Schreiber

Gross

Dyer

Nappo

Tyukov

Guerin

Hamadene

Wu

Stockbridge

Bobkov

Rutkowski

Kwapie\\xc5\\x84

Manou-Abi

Christen

Djehiche

Dwork
Albuquerque

Song

Mitzenmacher

Kumar

Marchioro
Bielecki

Maslov

Figure 4

Huet

Miller
Chistyakov

Left: General overview of the zbMath subgraph extracted for our experiments, containing approximately
13000 nodes. Right: Region of interest presenting the main results obtained on the zbMATH subgraph.

It is also possible to assert the robustness of our method with respect to several Monte Carlo
runs of our algorithm. We have produced some boxplots for each of the main authors located in the
subgraph, according to the occupation measure of the process over the last 10% of the iterations
with 10 Monte-Carlo replications. These “violin” plots are represented in Figure 5. Each execution
of the algorithm requires approximately 3 hours of computations. The algorithm seems to produce
reliable conclusions concerning the top nodes visited all along the ending iterations. Nevertheless,
it appears to be necessary to extend our investigations in order to obtain a scalable method for
handling larger graphs.

5.

Proof of the main result (Theorem 3)

We establish a differential inequality that will imply the convergence of Jt . The computations
actually lead us to a system of two differential inequalities. We therefore introduce another quantity

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

17

8
7

Frequency (%)

6
5
4
3

Figure 5

Nualart

Lions

Zeitouni

Brezis

Jaffard

Ledoux

Stroock

Yor

1

Cohen

2

“Violin” plot of the occupation measure of the algorithm on the zbMATH subgraph for 10 MC replications. The average frequency is located in the middle of the “violin” plot, although the minimal and
maximal values are shown in the extremity of the representation.

It that measures the average closeness (w.r.t. x) of the conditional law of y given x at time t to ν,
defined as:

Z
Z Z
mt (y |x)
It := KL(mt (.|x)||ν)dnt (x) =
log
mt (y |x)dy .dnt (x)
(13)
ν(y)
The next proposition links the evolution of Jt (in terms of an upper bound of ∂t Jt ) with the
spectral gap of µβt over ΓG , the diameter of the graph DG and the divergence It .
5.1.

Study of ∂t Jt

Proposition 2. Let c? (Uν ) be the term defined in Equation (34) and CΓG be the constant given
in Proposition 5. We then have:
?

∂t Jt ≤

2 0
DG
βt

2 2
+ 8DG
βt It

e−c (U )βt
Jt .
−
CΓG (1 + βt )

Proof: We compute the derivative of Jt and separately study each of the three terms:



nt (x)
∂t Jt = ∂t {log(nt (x))}dnt (x) − ∂t {log(µβt (x))}dnt (x) + log
∂t nt (x)dx
µβt
|
{z
} |
{z
} |
{z
}
Z

Z

J1,t

Study of J1,t : This term is easy to deal with:
Z
∂t {log(nt (x))}nt (x)

J1,t =
ΓG

Z

J2,t

J3,t

(14)

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

18

Z
=
ΓG

∂t {nt (x)}
nt (x) =
nt (x)


nt (x)dx = ∂t {1} = 0.

Z

Z
∂t {nt (x)}dx = ∂t

(15)

ΓG

ΓG

Study of J2,t : Using the definition of µβt given in Equation (9), we obtain for the second term:
Z

Z

J2,t = −
= βt0

∂t {log(µβt (x))}dnt (x) = −
ZΓG
Uν (x)dnt (x) +
ΓG

∂t {−βt Uν (x) − log Zβt }dnt (x)
ΓG

∂t {Zβt }
.
Zβt

According to the definition of Zβt , we have:
Z

Z
Z
e−βt Uν (x)
∂t {Zβt }
−βt Uν (x)
0
e
dx
=
−
β
U
(x)
= Zβ−1
∂
dx
=
−
∂
{
β
}
Uν (x)µβt (x)dx.
ν
t
t
t
t
t
Zβt
Zβt
ΓG
ΓG
ΓG
Hence, we obtain
J2,t = βt0

Z
Uν (x)[nt (x) − µβt (x)]dx
ΓG

Z
The graph has a finite diameter DG . We therefore have:

2
Uν (x)dnt (x) ≤ DG
. The same inequality

ΓG

holds using the measure µβt so that:

2 0
|J2,t | ≤ DG
βt .

(16)

Study of J3,t : The last term J3,t involves
Z the backward Kolmogorov equation. First, since nt is
the marginal law of Xt , we have: nt (x) = mt (x, y)dy.
Using the backward Kolmogorov equation for the Markov process (Xt , Yt )t≥0 and the Fubini
theorem, we have, for any smooth enough function ft : x ∈ ΓG 7−→ R :
Z

Z

Z



Z

Z

ft (x)∂t {nt (x)}dx =

ft (x)∂t
mt (x, y)dy dx =
V
ZΓG Z
Z
=
Lt (ft )(x)mt (x, y)dxdy =

ΓG

ΓG

V

ft (x)∂t {mt (x, y)}dxdy
ΓZ
G

ΓG

V

[L1,t + L2,t ](ft )(x)mt (x, y)dxdy,
V

where L1,t and L2,t are defined in Equations (6) and (7). Since the function ft is independent of y,
we have L1,t (ft ) = 0. For the part corresponding to L2,t , we have:
Z

Z

L2,t (ft )(x)mt (x, y)dxdy

Z 
1
=
∆x ft (x) − βt ∇x Uy (x)∇x ft (x) mt (x, y)dxdy
ZΓG V 2
Z Z
1
=
∆x ft (x)nt (x)dx − βt
∇x Uy (x)∇x ft (x)nt (x)mt (y |x)dxdy.
ΓG 2
ΓG V
ΓGZ V

because nt is the marginal distribution of Xt and mt (y |x) × nt (x) = mt (x, y).
Thus, for any smooth enough function ft (x), using the operator introduced in (12), we have:
Z
Z
ft (x)∂t nt (x)dx =
L̃2,t (ft )(x)nt (x)dx.
ΓG

ΓG

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

19

nt (x)
, we obtain:
µβt (x)


Z
Z
nt (x)
nt (x)
∂t nt (x) =
L̃2,t log
J3,t =
log
nt (dx).
µβt
µβt (x)
ΓG
ΓG

Replacing ft (x) by log

Since µβt is the invariant distribution of L̂2,t (see Equation (11)), it is natural to insert L̂2,t :


Z
nt (x)
J3,t =
L̃2,t log
nt (dx)
µβt (x) 


Z 
ZΓG

nt (x)
nt (x)
(17)
nt (dx) −
L̂2,t − L̃2,t log
nt (dx) .
=
L̂2,t log
µβt (x)
µβt (x)
ΓG
ΓG
|
{z
}
:=κt

Since L̂2,t is a diffusion operator and µβt is its invariant measure, it is well known (see, e.g., [7])
that the action of L̂2,t on the entropy is closely linked to the Dirichlet form in the following way:
(s
)!2


Z
Z
nt (x)
nt (x)
L̂2,t log
∇x
nt (dx) = −2
µβt (dx),
(18)
µβt (x)
µβt (x)
ΓG
ΓG
and therefore translates a mean reversion towards µβt in the first term of (17).
We now study the size of the difference between L̂2,t and L̃2,t and introduce the “approximation”
term of ∇x Uν (x) at time t:
Z
∇x (Uy )(x)mt (y |x)dy.

R̃t (x) =
V
√

x f
The relationship ∇x log(f ) = 2 ∇√
, the Cauchy-Schwarz inequality and 2ab ≤ a2 + b2 yield:
f



nt (x)
βt
R̃t (x) − ∇x Uν (x) ∇x log
nt (dx)
µβt (x) s
ΓG
s
(
)
Z 

nt (x)
µβt (x)
R̃t (x) − ∇x Uν (x) ∇x
nt (dx)
2βt
µβt (x)
nt (x)
ΓG
v
s
!2
sZ
uZ

2
u
nt (x)
t
2βt
R̃t (x) − ∇x Uν (x) nt (dx) ·
∇x
µβt (x)dx
µβt (x)
ΓG
ΓG
)!2
(s
Z 
Z
2
n
(x)
t
βt2
R̃t (x) − ∇x Uν (x) nt (x)dx +
µβt (dx)
∇x
µβt (x)
ΓG
Z

|κt | =

=

≤
≤





If dT V denotes the total variation distance, the first term of the right hand side leads to:
Z
∇x d2 (x, y)[mt (y |x) − ν(y)]dy
R̃t (x) − ∇x Uν (x) =
V
Z
2
≤ k∇x d (x, y)k∞
[mt (y |x) − ν(y)]dy
V

≤ 2k∇x d2 (x, y)k∞ dTs
V (mt (.|x), ν)
Z
√
mt (y |x)
2
≤ 2k∇x d (x, y)k∞
log
mt (y |x)dy
ν(y)
V

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

20

where the last line comes from the Csiszár-Kullback inequality. Since d2 (·, y) is differentiable a.e.
and its derivative is bounded for all y ∈ V by 2DG , we can use It defined in Equation (13) to obtain:
Z 
2
2
R̃t − R̂t (x)nt (x)dx ≤ 8DG
It
ΓG

Consequently, we obtain:
|κt | ≤

2 2
8DG
βt It

(s

Z
∇x

+
ΓG

nt (x)
µβt (x)

)!2
µβt (dx)

(19)

Taking the inequalities (18) and (19), we now obtain in (17):
(s
)!2
Z
n
(x)
t
2 2
J3,t ≤ 8DG
βt It −
∇x
µβt (dx)
µβt (x)
ΓG
We denote ft =

q

nt
µβt

. Since kf k22,µβ = 1, one can easily see that Entµβt (ft2 ) = Jt . Now, the Logat

rithmic Sobolev inequality on the (quantum) graph ΓG for the measure µβt stated in Proposition
5 shows that:
(s
)!2
Z
?
nt (x)
e−c (Uν )βt
∇x
µβt (dx) ≥
Jt ,
µβt (x)
CΓG (1 + βt )
ΓG
where c? (Uν ) is defined in Equation (34) and is related to the maximal depth of a well containing
a local but not global minimum. We thus obtain:


Z
?
e−c (Uν )βt
nt (x)
2 2
nt (dx) ≤ 8DG
β t It −
Jt .
J3,t = L̃2,t log
µβt (x)
CΓG (1 + βt )
The proof is concluded by regrouping the three terms.



5.2. Study of ∂t It
Proposition 3. Assume that DG ≥ 1 and βt is an increasing inverse temperature with βt ≥ 1,
then:
2
∂t It ≤ −αt It − ∂t Jt + DG
[βt0 + 6βt2 ].
(20)
Proof: We compute the derivative of It . Observing that mt (x, y) = nt (x)mt (y |x), we have:
Z
Z
 
mt (y |x)
∂t It = ∂t
nt (x)
log
mt (y |x)dy dx
ν(y)
V
Z ZΓG

mt (y |x)
= ∂t
log
mt (x, y)dxdy
ν(y)
Z ZΓG V
Z Z
Z
mt (y |x)
=
∂t {log mt (y |x)}mt (x, y)dxdy −
∂t {log ν(y)}mt (x, y)dxdy + log
∂t mt (x, y)dxdy
ν(y)
ΓG V
ΓG V
{z
}
|
{z
} |
{z
} |
:=I1,t

:=I3,t

:=I2,t

The computation of the first two terms is straightforward. For the first one we have:
Z Z
Z Z
∂t mt (y |x)
I1,t =
∂t {log mt (y |x)}mt (x, y) =
mt (x, y)dxdy

ZΓG ZV
Z ΓG VZmt (y |x)
=
∂t mt (y |x)nt (x)dxdy =
nt (x)
∂t mt (y |x)dy dx
ΓG

V

ΓG

V

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

Z
=

nt (x)∂t
ΓG




Z


| V

21






mt (y |x)dy dx = 0


{z
}
:=1

The computation of I2,t is easy since ν does not depend on t, implying that I2,t = 0.
For the third term, we use the backward Kolmogorov equation and obtain:


Z Z
Z Z
mt (y |x)
mt (y |x)
I3,t =
log
∂t mt (x, y) =
Lt log
mt (x, y)dx dy
ν(y)
ν(y)
Γ
V
Γ
V
G
G




Z Z
Z Z
mt (y |x)
mt (y |x)
=
L1,t log
L2,t log
mt (x, y)dx dy +
mt (x, y)dx dy
ν(y)
ν(y)
ΓG V
ΓG V
{z
} |
{z
}
|
1
:=I3,t

2
:=I3,t

The jump part L1,t exhibits a mean reversion on the entropy on the conditional law: applying the
Jensen inequality for the logarithmic function and the measure ν, we obtain:
1
I3,t



mt (y 0 |x)
mt (y |x)
0
0
= αt
log
− log
ν(y )dy mt (x, y)dxdy
ν(y 0 ) 
ν(y)

ZΓG ZV ZV 
mt (y 0 |x)
log
= αt
ν(y 0 )dy 0 mt (x, y)dxdy − αt It
0)
ν(y

ZΓG ZV VZ
mt (y 0 |x)
0
0
≤ αt
log
ν(y )dy mt (x, y)dxdy − αt It
0
ZV ν(y )

ZΓG ZV
0
0
≤ αt
log
mt (y |x)dy mt (x, y)dxdy − αt It
V
ZΓG ZV
≤ αt
log 1 · mt (x, y)dxdy − αt It ≤ −αt It
Z

ΓG

Z Z 

(21)

V

If we consider the action of L2,t on the entropy of the conditional law, using mt (x, y) = mt (y |x)nt (x)
yields:


Z Z
mt (y |x)
2
mt (x, y)dx dy
I3,t =
L2,t log
ν(y)


ZΓG ZV
mt (x, y)
=
L2,t log
mt (x, y)dx dy
ν(y) · nt (x)
ZΓG ZV
Z Z
=
L2,t log (mt (x, y)) mt (x, y)dxdy −
L2,t log (nt (x, y)) mt (x, y)dxdy,
(22)
ΓG

V

ΓG

V

because L2,t (ν)(y) = 0 (L2,t only involves the x component). We now study the first term of (22):
Z

Z
ΓG

V

Z Z
1
∆x [log mt (x, y)] dmt (x, y)
L2,t (log mt (x, y)) dmt (x, y) =
2 ZΓG ZV
− βt
< ∇x Uy , ∇x log mt (x, y) > dmt (x, y)
Z ΓG Z V
1
=
∂ 2 {log mt (x, y)} dmt (x, y)
2 ZΓG V x
− βt ∂x Uy ∂x log mt (x, y)dmt (x, y)

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

22

The first term of the right-hand side deserves special attention:


Z Z
Z Z
∂x {mt (x, y)}
2
∂x {log mt (x, y)} dmt (x, y) =
∂x
dmt (x, y)
mt (x, y)
ΓG V
Z Z
ZΓG ZV 2
[∂x {mt (x, y)}]2
∂x {mt (x, y)}
dmt (x, y) −
dxdy
=
ZΓG ZV mt (x, y)
Z ZΓG  V p mt (x, y)2
=
∂x2 {mt (x, y)}dxdy − 4
∂x mt (x, y) dxdy
ΓG Z V
ΓG V
X
=
[∂x mt (e(Le ), y)] − ∂x mt (e(0), y)dy
V
e∈E
Z Z  p
2
−4
∂x mt (x, y) dxdy,
ΓG

V

where {e(s), 0 ≤ s ≤ Le } is the parametrization of the edge e ∈ E introduced in Section 3.1. The
gluing conditions (5) yield:
#
Z X "X
XZ
∂x mt (e(Le ), y) − ∂x mt (e(0), y) =
−de mt (v, y) = 0.
e∈E

V

V v∈V

e∼v

As a consequence, we obtain:
Z Z
Z Z  p
2
∂x mt (x, y) dxdy
L2,t (log mt (x, y)) dmt (x, y) = −2
ΓG V
Z ΓGZ V
− βt
∂x Uy ∂x log mt (x, y)dmt (x, y)
ΓG

(23)

V

The Cauchy-Schwarz inequality applied on the second term leads to:
sZ Z
Z Z
2
βt
∂x Uy ∂x log mt (x, y)dmt (x, y) ≤ βt k∂x d2 (., .)k2,mt
(∂x log mt (x, y)) dmt (x, y)
ΓG V
sZ ΓG Z V
2

≤ βt ||∂x d2 (x, y)||∞
(∂x log mt (x, y)) mt (x, y)dxdy
ΓG V
sZ Z
 p
2
≤ 4βt DG
∂x mt (x, y) dxdy
ΓG V
Z Z  p
2
2 2
≤ 2DG βt + 2
∂x mt (x, y) dxdy
ΓG

V

Inserting this in Equation (23) leads to:
Z Z
2 2
L2,t log mt (x, y)dmt (x, y) ≤ 2DG
βt .
ΓG

(24)

V

We study the second term of (22) and use the proof of Proposition 2: the decomposition (14) with
Equations (15) and (16) yield:


Z
nt (x)
0 2
∂t Jt ≤ βt DG + L2,t log
mt (dx, dy).
µβt (x)
This implies:
Z Z
Z
2
−
L2,t (log nt (x)) mt (dx, dy) ≤ −∂t Jt + βt0 DG
− L2,t (log µβt (x)) mt (dx, dy).
ΓG

V

(25)

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

23

Using the definition of µβt , we obtain:
Z Z
Z Z
e−βt Uν (x)
−
L2,t log
=
[L2,t (βt Uν (x)) + L2,t (log Zβt )] mt (dx, dy)
Zβt
ΓG V
ΓG
Z VZ
L2,t (Uν (x)) mt (dx, dy)

= βt
ZΓG ZV
= βt
ΓG

V

1
∆x Uν (x) − βt ∇x Uy (x)∇x Uν (x)mt (dx, dy)
2

βt
2
≤ + 4βt2 DG
2
This inequality used in (25) yields:
Z
βt
2
2
− L2,t (log nt (x)) mt (x, y) = −∂t Jt + βt0 DG
+ + 4βt2 DG
2

(26)

We now use (24) and (26) in (22) and our assumptions βt ≥ 1 and DG ≥ 1 to obtain:
2
2
I3,t
≤ −∂t Jt + DG
[βt0 + 6βt2 ]

(27)

Combining (21) and (27) leads to the desired inequality given by Equation (20).



5.3. Convergence of the entropy
The use of Propositions 2 and 3 makes it possible to obtain a system of coupled differential inequalities.
2
Proof of Theorem 3: If we denote a = DG
, we can write:
(
−c? (Uν )βt
Jt0 ≤ − Ce Γ (1+βt ) Jt + aβt0 + 8aβt2 It
G

It0 ≤ −αt It − Jt0 + a(βt0 + 6βt2 )
We introduce an auxiliary function Kt = Jt + kt It , where kt is a smooth positive decreasing function
for t large enough, so that lim kt = 0. We use the system above to deduce that:
t→∞

Kt0 = Jt0 + kt0 It + kt It0
≤ Jt0 + kt It0
≤ Jt0 + kt (−αt It − Jt0 ) + akt (βt0 + 6βt2 )
≤ (1 − kt )Jt0 − kt αt It + akt (βt0 + 6βt2 ).
In the first inequality, we use the fact that kt0 ≤ 0 and It is positive. The second inequality is given
by (20). The upper bound of Proposition 2 now leads to:
Kt0 ≤ −t (1 − kt )Jt − kt αt It + 8a(1 − kt )βt2 It + aβt0 + 6akt βt2 ,
where we denoted t =

?
e−c (Uν )βt
.
CΓ (1+βt )
G

We choose the function kt to obtain a mean reversion on Kt :
kt :=

8aβt2
.
αt + 8aβt2 − t /2

(28)

Note that this function is decreasing for sufficiently large t as soon as βt = o(αt ), which is the
case according to the choices described in Theorem 3. Moreover, a straightforward consequence is
limt−→+∞ kt = 0. This ensures that a positive T0 exists such that:
∀t ≥ T0

1
0 ≤ kt ≤ .
2

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

24

Consequently, we deduce that:
∀t ≥ T0

t
Kt0 ≤ −t (1 − kt )Jt − kt It + aβt0 + 6aβt2 kt
2
t
t
t
≤ − Jt − kt It + aβt0 + 6aβt2 kt ≤ − Kt + aβt0 + 6aβt2 kt
{z
}
|
2
2
2
:=ηt

The next bound is an easy consequence of the Gronwall Lemma:
Z t
Z t
s
u
−
ds Z T
−
du
2
∀t ≥ T0
Kt ≤ KT0 e T0
+
ηs e s 2 ds,
T0

which in turn implies that Kt −→ 0 as t −→ +∞ as soon as ηt = ot∼+∞ (t ) with

R∞
0

u du = +∞.

We are now looking for a suitable choice for αt and βt . Let us assume that αt ∼ λtγ for any γ > 0.
This choice leads to:
8ab2 log(t + 1)2
kt ∼
,
λtγ
so that:
ab 48a2 b4 log(t + 1)4
ηt ∼
+
= O(t−(1∧γ) log(t)4 ).
t
λtγ
?
t−c (Uν )b
.
bCΓ log t

At the same time, we can check that t ∼

Now, our conditions on ηt and t imply that:

G

∞

Z

u du = +∞

⇐⇒

b c? (Uν ) ≤ 1.

0

At the same time:
⇐⇒

ηt = o(t )

1 ∧ γ > b c? (Uν ).

The optimal calibration of our parameters γ and b (minimal value of γ, maximal value of b) induces
the choice γ = 1 and b < c? (Uν )−1 . Up to the choices αt ∼ λt and βt ∼ b log(t), we deduce that
Kt −→ 0 when t goes to infinity. Since It , kt and Jt are positive, this also implies that Jt −→ 0
when t goes to infinity.


6.

Functional inequalities

This section is devoted to the proof of the Log-Sobolev inequality for the measure (µβ )β>0 . We dealt
with the Dirichlet form using this inequality (see Equation (18)) in Proposition 2. In particular, we
need to obtain an accurate estimate when β −→ +∞. For this purpose, we introduce the generic
notation for Dirichlet forms (see, i.e., [7] for a more in-depth description):
Z
∗
1,2
2
∀β ∈ R+ ∀f ∈ W (µβ )
Eβ (f, f ) = kOf kµβ =
|Of (x)|2 dµβ (x).
ΓG

Z
If we denote < f >µβ = µβ (f ) =

f dµβ , we are interested in showing the Poincaré inequality:
ΓG

kf − < f >µβ k22,µβ =

Z
ΓG

(f − < f >µβ )2 dµβ ≤ λ(β)Eβ (f, f )

(29)

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

25

and the Log-Sobolev Inequality (referred to as LSI below):
!2
Z
f
f 2 log
dµβ ≤ C(β)Eβ (f, f )
kf k2,µβ
ΓG

(30)

The specific feature of this functional inequality deals with the quantum graph settings and
deserves careful adaptation of the pioneering work of [31]. This technical section is split into two
parts. The first one establishes a preliminary estimate when β = 0 (i.e., when dealing with the
uniform measure on ΓG ). The second one then uses this estimate to derive the asymptotic behavior
of the LSI when β −→ +∞.
6.1. Preliminary control for µ0 on ΓG
We consider µ0 the normalized Lebesgue measure and use the standard notation for any measure
µ on ΓG :
kf kW 1,p (µ) := kf kp,µ + k∇f kp,µ
Let us establish the next elementary result:
Lemma 1. The ordinary Sobolev inequality holds on ΓG , i.e., for all measurable functions f we
have:


2
2/e 2 DG L
kf − < f >0 kp,µ0 ≤ e L
+ 1 E0 (f, f ),
(31)
2
where DG is the diameter of the graph and L its perimeter defined by:
X
Le .
L=
e∈E

Proof: First, let us remind the reader that for a given interval I in dimension 1 equipped with the
Lebesgue measure λI , the Sobolev space W 1,p (λI ) is continuously embedded in L∞ (I) (compact
injection when I is bounded). In particular, it can be shown (see, e.g., [16]) that while integrating
w.r.t. the unnormalized Lebesgue measure:
∀p ≥ 1

kf kL∞ (I) ≤ e1/e kf kW 1,p (λI ) .

We now consider f ∈ W 1,p (µ0 ). Since G has a finite number of edges we can write: ΓG =
where each ei can be seen as an interval of length `i . We have seen that for each edge ei :

Sn

i=1 ei ,

kf kL∞ (ei ) ≤ e1/e kf kW 1,p (λei ) .

We can use a union bound since ΓG represents the union of edges, and deduce that:
kf kL∞ (ΓG ) = max kf kL∞ (ei ) ≤ e1/e max kf kW 1,p (λei ) ≤ e1/e kf kW 1,p (λΓ ) ,
1≤i≤n

1≤i≤n

G

where the inequality above holds w.r.t. the unnormalized Lebesgue measure. If L denotes the sum
of the lengths of all edges in ΓG , we obtain:
kf kL∞ (ΓG ) ≤ e1/e Lkf kW 1,p (µ0 ) .

(32)

Second, we establish a simple Poincaré inequality for µ0 on ΓG . For any function f ∈ W 1,2 (µ0 ),
we use the equality:
#2
Z Z
Z Z "Z
1
1
f 0 (s)ds dµ0 (x)dµ0 (y).
kf − < f >0 k22,µ0 =
[f (x) − f (y)]2 dµ0 (x)dµ0 (y) =
2 ΓG ΓG
2 ΓG ΓG γx,y

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

26

where γx,y is the shortest path that connects x to y, parametrized with speed 1 and f 0 (s) refers to
the derivative of f w.r.t. this parametrization at time s. It should be noted that such a path exists
because the graph ΓG is connected. The Cauchy-Schwarz inequality yields:
#
Z Z "Z
1
DG L
kf − < f >0 k22,µ0 ≤
|∇f (s)|2 ds |γx,y |dµ0 (x)dµ0 (y) ≤
k∇f k22,µ0 .
(33)
2 ΓG ΓG γx,y
2
The Sobolev inequality is now an obvious consequence of the previous inequality: consider f ∈
W 1,p (µ0 ) and note that since µ0 (ΓG ) = 1, then (32) applied with p = 2 leads to:
2/q
Z
|f − < f >0 |q dµ0
∀q ≥ 1
kf − < f >0 k2q,µ0 =
ΓG

≤ kf − < f >0 k2L∞ (ΓG ) µ0 (ΓG )
≤ e2/e L2 kf − < f >0 k2W 1,2 (µ0 ) .

We can use the Poincaré inequality established for µ0 on ΓG in Equation (33) and obtain:




2
2/e 2 DG L
2
2/e 2 DG L
∀q ≥ 1
kf − < f >0 kq,µ0 ≤ e L
+ 1 k∇f k2,µ0 = e L
+ 1 E0 (f, f ),
2
2
which concludes the proof.

Remark 1. The constants obtained in the proof of Lemma 1 above could certainly be improved.
Nevertheless, such an improvement would have little importance for the final estimates obtained
in Proposition 5.
6.2. Poincaré Inequality on µβ
In the following, we show a Poincaré Inequality for the measure µβ for large values of β. This
preliminary estimate will be useful for deriving LSI on µβ . This functional inequality is strongly
related to the classical minimal elevation of the energy function Uν for joining any state x to any
state y.
We first introduce some useful notations. For any couple of vertices (x, y) of ΓG , and for any
path γx,y that connects them, we define h(γx,y ) as the highest value of Uν on γx,y :
h(γx,y ) = max Uν (s).
s∈γx,y

We define H(x, y) as the smallest value of h(γx,y ) obtained for all possible paths from x to y:
H(x, y) = min h(γ)
γ:x→y

Now, for any pair of vertices x and y, the notation γx,y will be reserved for the path that attains
the minimum in the definition of H(x, y). Such a path exists for any x, y because ΓG is connected
and possesses a finite number of paths that connect any two given vertices.
Finally, we introduce the quantity that will mainly determine the size of the spectral gap involved
in the Poincaré inequality and the constant in the LSI (see the seminal works of [25] for a LDP
probabilistic interpretation and [32] for a functional analysis point of view):
c? (Uν ) := max [H(x, y) − Uν (x) − Uν (y)] + min Uν (x)
(x,y)∈Γ2
G

x∈ΓG

In the following, every time we write min Uν we refer to min Uν (x).
x∈ΓG

(34)

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

27

Theorem 4 (Poincaré inequality for µβ ). For all measurable functions g defined on ΓG :
V arµβ (g) ≤

|E | maxe∈E Le e2DG βc? (Uν )
e
Eµβ (g, g).
2

Proof: Since the graph ΓG is connected, for any two points x and y, we can find a minimal path
γx,y that links them and that minimizes H. We denote x0 = x, x1 , · · · xl , xl+1 = y where the sequence
x1 , · · · , xl refers to the nodes included in the path γx,y .
Z Z
1
(g(y) − g(x))2 dµβ (x)dµβ (y)
V arµβ (g) =
2 ΓG ΓG
!2
Z Z
l
X
1
=
(g(xi+1 ) − g(xi ) dµβ (x)dµβ (y)
2 ΓG ΓG i=0
!2
Z Z
l Z xi+1
X
1
0
g (s)ds dµβ (x)dµβ (y)
=
2 ΓG ΓG i=0 xi
!2
Z Z
Z
1
g 0 (s)ds dµβ (x)dµβ (y)
=
2 ΓG ΓG
γx,y
Z Z Z
1
≤
|∇g(s)|2 ds|γx,y |dµβ (x)dµβ (y),
2 ΓG ΓG γx,y
The last line is implied by the Cauchy-Schwarz inequality.
We can organize the terms involved in the above upper bound in the following way: for any edge
of the graph e ∈ E, we denote Ve the set of points (x, y) such that γx,y ∩ e 6= ∅. We therefore have:
Z
Z
1X
2
|γx,y |dµβ (x)dµβ (y)
|∇g(s)| ds
V arµβ (g) ≤
2 e∈E e
Ve
Z
Z
1X
1
≤
|∇g(s)|2 µβ (s)
|γx,y |
dµβ (x)dµβ (y)ds
2 e∈E e
µ
β (s)
Ve

Z

Z
1
1X
2
|γx,y |dµβ (x)dµβ (y)
|∇g(s)| µβ (s)ds
sup
≤
2 e∈E e
s̃∈e µβ (s̃) Ve
Z
Z
1
1
≤
|∇g(s)|2 µβ (s)ds sup
|γx,y |dµβ (x)dµβ (y)
2 ΓG
s̃∈ΓG µβ (s̃) Ve
s̃
In this case, es̃ denotes the edge of the graph that contains the point s̃ and the set Ves̃ is still
the set of couples (x, y) defined above, associated with each edge es̃ . We introduce the quantity A
defined as:
Z
1
A = sup
|γx,y |dµβ (x)dµβ (y).
s̃∈ΓG µβ (s̃) Ve
s̃
Using this notation, we have obtained that for all functions g, we have the Poincaré inequality:
V arµβ (g) ≤

A
Eµ (g, g).
2 β

All that remains to be done is to obtain an upper bound of A.
Z
1
A = sup
|γx,y |dµβ (x)dµβ (y)
s̃∈ΓG µβ (s) Ve
s̃

(35)

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

28

e−β(Uν (x)+Uν (y))
= sup Zβ e
|γx,y |dxdy
Zβ2
s̃∈ΓG
Ves̃
Z
1
eβ(Uν (s̃)−Uν (x)−Uν (y)) |γx,y |dxdy
sup
=
Zβ s̃∈ΓG Ves̃
βUν (s̃)

Z

Since γx,y is the minimal path for H(x, y), we have H(x, y) = max Uν (s). Therefore:
s∈γx,y

A≤
≤
≤
≤

Z
1
eβ((H(x,y)−Uν (x)−Uν (y)) |γx,y |dxdy
sup
Zβ s̃∈ΓG ZVes̃
?
1
eβ(c (Uν )−min(Uν )) |γx,y |dxdy
sup
Zβ s̃∈ΓG Ves̃
Z
?
eβ(c (Uν )−min(Uν ))
|γx,y |dxdy
sup
Zβ
s̃∈ΓG Ve
s̃
?
eβ(c (Uν )−min(Uν ))
|E | max Le .
e∈E
Zβ

Using the definition of Zβ , we have:
eβc
A≤ Z

?

(Uν )

e

|E | maxe∈E Le

−β(Uν (x)−min(Uν ))

.

(36)

dx

ΓG




1
If x ∈ Mν is a Fréchet mean that minimizes Uν , we designate B x ,
, as the ball of center x?
β
1
and radius , for the geodesic distance d on the graph ΓG . It is easy to check that:
β


|Uν (x) − Uν (x? )| = EY ∼ν d2 (x, Y ) − d2 (x? , Y ) ≤ EY ∼ν d2 (x, Y ) − d2 (x? , Y ) ≤ 2DG × d(x, x? ).
?

?

We can then deduce a lower bound on the denominator involved in (36):
Z
Z
−β(Uν (x)−min(Uν ))
e
dx ≥  1  e−β(Uν (x)−min(Uν )) dx
Bx? ,

ΓG

β



Z
≥

1 e
Bx? , 
βZ


≥ e−2DG



−2DG βd(x,x? )



dx

1 dx.

β


Bx? ,




1
The Lebesgue measure of B x ,
may be lower bounded by β −1 since there is, at the least, one
β
path in ΓG passing by the point x? . Inserting this inequality in (36) gives:
?

A ≤ |E | max Le e2DG βeβc
e∈E

?

(Uν )

.

Using this upper bound in (35) leads to the desired Poincaré inequality.



Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

29

6.3. Sobolev Inequalities on µβ
6.3.1. Preliminary control on Dirichlet forms The next result will be useful to derive
a LSI for µβ from a Poincaré inequality on µβ (given in Theorem 4). It generalizes the Poincaré
inequality for p norms with p > 2 while using the Sobolev inequality given in Lemma 1.
First, we introduce the maximal elevation of Uν as:
M = sup Uν (x) − inf Uν (x).
x∈ΓG

x∈ΓG

(37)

2
Note that in our case, M may be upper bounded by DG
.

Proposition 4. For any p > 2 and all measurable functions f , we have:


DG L
2
2
∀β ≥ 0
kf − < f >β kp,µβ ≤ 4L
+ 1 e2/e+M β Eβ (f, f ).
2

(38)

Proof: The Jensen inequality applied to the convex function x 7−→ xp yields:
Z
|f − < f > |pβ dµβ
kf − < f >β kpp,µβ =
ZΓG
Z
p
=
f (x) −
f (y)dµβ (y) dµβ (x)
ΓG
ZΓG Z
p
=
f (x) − f (y)dµβ (y) dµβ (x)
ZΓG Z ΓG
p
≤
|f (x) − f (y)| dµβ (y)dµβ (x)
ZΓG ΓG
p
≤
(|f (x)− < f >0 | + |f (y)− < f >0 |) dµβ (y)dµβ (x)
ΓG

Again, the Jensen inequality |a + b|p ≤ 2p−1 [|a|p + |b|p ] implies that:
Z
p
p−1
kf − < f >β kp,β ≤ 2
|f (x)− < f >0 |p + |f (y)− < f >0 |p dµβ (y)dµβ (x)
Γ
Z G
≤ 2p
|f (x)− < f >0 |p dµβ (x) = 2p kf − < f >0 kpp,µβ
ΓG

We conclude that:
kf − < f >β k2p,µβ ≤ 4kf − < f >0 k2p,µβ .

Using the fact that:
µβ (x) =

(39)

e−β min Uν
e−β min Uν
≤
× Z0 µ0 (x),
Zβ
Zβ

Le−β min Uν
kf − < f >0 kpp,µ0 , since Z0 is the perimeter L of the
Zβ
graph ΓG . Consequently, we obtain:

2/p  Le−β min Uν 2/p
2
p
kf − < f >0 kp,µβ = kf − < f >0 kp,µβ
≤
kf − < f >0 k2p,µ0
Zβ
we also have: kf − < f >0 kpp,µβ ≤

Using inequality (39), the fact that Zβ ≤ Le−β min Uν and the assumption 2/p < 1, we conclude that:
kf − < f >β k2p,µβ ≤

4Le−β min Uν
kf − < f >0 k2p,µ0 .
Zβ

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

30

The Sobolev inequality (1) implies that:
kf − < f >β k2p,µβ ≤

4L3 (DG L/2 + 1) e2/e−β min Uν
E0 (f, f ).
Zβ

(40)

We now find an upper bound for the Dirichlet form E0 (f, f ) that involves Eβ (f, f ):
Z
E0 (f, f ) =
< Of, Of > dµ0
ΓGZ
e−βUν (x)
< Of, Of > eβUν (x)
dµ0
= Zβ
Zβ
ΓG 

Zβ exp β sup Uν (x) Z
x∈ΓG
≤
< Of, Of > dµβ
 Z0
 ΓG
Zβ exp β sup Uν (x)
x∈ΓG
≤
Eβ (f, f )
L

(41)

Putting (40) and (41) together concludes the proof.



6.3.2. Log-Sobolev Inequality on µβ For all probability measures µ and all measurable
functions f , we denote:


Z
f2
2
2
Entµ (f ) =
f log
dµ
kf k22,µ
ΓG
Proposition 5. The Log-Sobolev Inequality holds on ΓG . A constant CΓG exists such that for
all β ≥ 0 and all µβ -measurable functions f , we have:
?

Entµβ (f 2 ) ≤ CΓG [1 + β]ec

(Uν )β

Eβ (f, f )

Proof: We consider β ≥ 0 and a µβ measurable function f . We apply the Jensen inequality for the
logarithmic function and the measure f 2 /kf k22,µβ dµβ to obtain:
Z
ΓG

f
kf k2,µβ

!2
log

f2
kf k22,µβ

!

!2
!
f
|f |p−2
log
dµβ
p−2
kf k2,µβ
kf k2,µ
ΓG
β!
Z
f2
|f |p−2
2
≤
log
dµβ
p−2
2
p−2
ΓG kf k2,,µβ kf k2,,µβ
!
!
kf kpp,µβ
kf k2p,µβ
2
p
≤
log
=
log
.
p−2
kf kp2,µβ
p−2
kf k22,µβ

2
dµβ =
p−2

Z

Observing that for all x, δ > 0 we have log(xδ) ≤ δx, we obtain log(x) ≤ xδ + log(δ −1 ). Therefore:
Z
∀δ > 0

2

f log
ΓG

f
kf k2,µβ

!2
dµβ =
≤

kf k22,µβ

Z

pkf k22,µβ

ΓG

!2
f
log
kf k2,µβ
!
kf k2p,µβ

log
p−2
kf k22,µβ
"
#
2
pkf k2,µβ
kf k2p,µβ
1
≤
δ
+ log
.
p−2
kf k22,µβ
δ

f2
kf k22,µβ

!
dµβ

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

31

Replacing f with f − < f >β and choosing δ = e−M β leads to:

2
Z
f − < f >β
2
(f − < f >β ) log
dµβ
kf −
< f >β k2,β
ΓG
h
i
p
≤
e−βM kf − < f >β k2p,µβ + M β kf − < f >β k22,µβ .
p−2
Proposition 4 and the Poincaré inequality established in Theorem 4 yields:



Entµβ (f − < f >β )


2

f − < f >β
kf − < f >β k2,µβ

Z

!2

dµβ
(f − < f >β )2 log
ΓG


p
e2DG c? (Uν )β
2 DG L + 2 2/e
≤
4L
e + M β |E | max Le
βe
Eβ (f, f ) (42)
e∈E
p−2
2
2

=

It remains to use Rothau’s Lemma (see Lemma 5.1.4 of [7]) that states that for any measure µ and
any constant a:
Z
Entµ (g + a)2 ≤ Entµ g 2 + 2

g 2 dµ.

ΓG

Let p = 4 in Equation (42). Putting this together with Rothau’s Lemma and Theorem 4, we obtain
the following:
!
Z
2
f
Entµβ (f 2 ) =
f 2 log
dµβ
kf k22,µβ
ΓG


2
≤ Entµβ (f
i
h − < f >β ) + 2V arµβ (f ) ?
≤ Eβ (f, f ) β |E | max Le e2DG (1 + M β)ec (Uν )β + 4L2 (2 + DG L)e2/e
?

≤ CΓG [1 + β]ec

e∈E
(Uν )β

Eβ (f, f ),

where CΓG is a large enough constant (independent of β) that could be made explicit in terms of
2
constants DG , |E | and L since we trivially have M ≤ DG
.


References
[1] Acemoğlu D., Fagnani F. Ozdaglar A., Como G. 2013. Opinion fluctuations and disagreement in social
networks. Math. Oper. Res. 38(1) 1–27.
[2] Allassonnière S., Trouvé A., Kuhn E. 2010. Construction of bayesian deformable models via stochastic
approximation algorithm: A convergence study. Bernoulli 16 641–678.
[3] Arnaudon M., Miclo L. 2014. Means in complete manifolds: uniqueness and approximation. ESAIM:
Probability and Statistics 18 185–206.
[4] Arnaudon M., Miclo L. 2014. A stochastic algorithm finding generalized means on compact manifolds.
Stochastic Processes and their Applications 124 34633479.
[5] Arnaudon M., Miclo L. 2016. A stochastic algorithm finding p-means on the circle. Bernoulli, in press
.
[6] Bach F., Jordan M. 2004. Learning spectral clustering. Advances in Neural Information Processing
Systems 305–312.
[7] Bakry D., Ledoux M., Gentil I. 2014. Analysis and Geometry of Markov Diffusion Operators,
Grundlehren der mathematischen Wissenschaften, vol. 348. Springer.
[8] Barden D., Owen M., Le H. 2013. Central limit theorems for Fréchet means in the space of phylogenetic
trees. Electronic Journal of Probability 18 1–25.

32

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

[9] Bhattacharya R., Patrangenaru V. 2003. Large sample theory of intrinsic and extrinsic sample means
on manifolds. Ann. Statist. 31 1–29.
[10] Bigot, J. 2013. Fréchet means of curves for signal averaging and application to ECG data analysis.
Annals of Applied Statistics 7 1837–2457.
[11] Bigot J., Gadat S. 2010. A deconvolution approach to estimation of a common shape in a shifted curves
model. Annals of Statistics 38 224–243.
[12] Bigot J., Gendre X. 2013. Minimax properties of Fréchet means of discretely sampled curves. Annals
of Statistics 41 923–956.
[13] Bigot J., Klein T. Lopez A., Gouet R. 2015. Geodesic PCA in the wasserstein space by convex PCA.
Annales de l’Institut Henri Poincaré B: Probability and Statistics, to appear .
[14] Bigot J., Lopez A., Gouet R. 2013. Geometric PCA of images. Siam, J. Imaging Sci. 6 1851–1879.
[15] Bontemps D., Gadat S. 2014. Bayesian methods for the shape invariant model. Electronic Journal of
Statistics 8 1522–1568.
[16] Brezis, H. 1987. Analyse fonctionelle. Masson, Paris.
[17] Catoni, O. 1992. Rough large deviation estimates for simulated annealing : application to exponential
schedules. Ann. Probab. 20 1109–1146.
[18] Dijkstra, E. W. 1959. A note on two problems in connexion with graphs. Numer. Math. 1 269–271.
[19] Dryden I. L., Mardia K. V. 1998. Statistical Shape Analysis.
[20] Erdös P., Rényi A. 1960. The evolution of random graphs. Magyar Tud. Akad. Mat. Kutato Int. Kozl.
5 17–61.
[21] Estrada, E. 2015. Introduction to Complex Networks. Structure and Dynamics, chapter of Evolutionary
Equations with Applications to Natural Sciences. Lecture Notes in Mathematics, Springer.
[22] Ethier S. N., Kurtz T. 2005. Markov processes. Characterization and convergence.. Wiley Series in
Probability and Statistics, John Wiley & Sons.
[23] Fréchet, M. 1948. Les éléments aléatoires de nature quelconque dans un espace distancié. Annales de
l’Institut Henri Poincaré (B) 10 215–310.
[24] Freidlin M., Sheu S.J. 2000. Diffusion processes on graphs: stochastic differential equations, large deviation principle. Probab. Theory Relat. Fields 116 181–220.
[25] Freidlin M. I., Wentzell A. D. 1979. Random perturbations of dynamical systems, Grundlehren der
Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences] , vol. 260. 2nd ed.
Springer-Verlag, New York. Translated from the 1979 Russian original by Joseph Szücs.
[26] Freidlin M. I., Wentzell A. D. 1995. Random perturbations of Hamiltonian systems, Memoirs of the
American Mathematical Society, vol. 523. A.M.S., New York.
[27] Ginestet, C.E. 2013. Strong consistency of set-valued fréchet sample means in metric spaces. Preprint .
[28] Goldenberg A., Fienberg S.E. Airoldi E. M., Zheng A. X. 2010. A survey of statistical network models.
Found. Trends Mach. Learn. 2(2) 129–233.
[29] Gutjahr W. J., Pflug G. 1996. Simulated annealing for noisy cost functions. J. Global Optim. 8(1) 1–13.
[30] Hajeck, B. 1988. Cooling schedules for optimal annealing. Mathematics of Operations Research 13
311–329.
[31] Holley, Richard, Daniel Stroock. 1988. Simulated annealing via Sobolev inequalities. Comm. Math.
Phys. 115(4) 553–569.
[32] Holley R., Stroock D., Kusuoka D. 1989. Asymptotics of the spectral gap with applications to the theory
of simulated annealing. J. Funct. Anal. 83(2) 333–347.
[33] Ikeda N., Watanabe S. 1981. Stochastic Differential Equations and Diffusion Processes. North-Holland.
[34] Jackson, M. O. 2008. Social and Economic Networks. Princeton Univ. Press, Princeton, NJ.
[35] Kaufmann M., Wagner D. 2001. Drawing graphs: methods and models. Lecture Notes in Computer
Science, Springer.

Author: On the calculation of a graph barycenter
00(0), pp. 000–000, c 0000 INFORMS

33

[36] Klopp O., Verzelen N., Tsybakov A.B. 2016. Oracle inequalities for network models and sparse graphon
estimation. Annals of Statistics, in press .
[37] Kolaczyk, E. D. 2009. Statistical analysis of network data: methods and models. Springer Series in
Statistics, Springer New-York.
[38] Le, H. 2001. Locating Fréchet means with application to shape spaces. Adv. Appl. Probab. 33 324–338.
[39] Lovasz, L. 2012. Large networks and graph limits, vol. 60. American Mathematical Society.
[40] Miclo, Laurent. 1992. Recuit simulé sur Rn . Étude de l’évolution de l’énergie libre. Ann. Inst. H.
Poincaré Probab. Statist. 28(2) 235–266.
[41] Miller E., Provan J., Owen M. 2015. Polyhedral computational geometry for averaging metric phylogenetic trees. Advances in Applied Mathematics 68 51–91.
[42] Munch E., Bendich P. Mukherjee S. Mattingly J. Harer J., Turner K. 2015. Probabilistic Fréchet means
for time varying persistence diagrams. Electronic Journal of Statistics 9 1173–1204.
[43] Newman, M. 2010. Networks, An Introduction. Oxford University Press.
[44] Pennec, X. 2006. Intrinsic statistics on Riemannian manifolds: basic tools for geometric measurements.
J. Math. Imaging Vision 25 127–154.
[45] Shneiderman B., Aris A. 2006. Network visualization by semantic substrates. IEEE Transactions on
Visualization and Computer Graphics 12 733–740.
[46] Trouvé, A. 1993. Parallélisation massive du recuit simulé. PhD Thesis, Université d’Orsay .
[47] Watts D. J., Strogatz S. H. 1998. Collective dynamics of ’small-world’ networks. Nature 393(6684)
409–10.

