∎

11institutetext: C.B. Scott 22institutetext: E. Mjolsness 33institutetext:
University of California, Irvine, Dept. of Computer Science, Irvine, CA, USA
Cory Scott: scottcb @ uci.edu Eric Mjolsness: emj @ uci.edu

# Novel diffusion-derived distance measures for graphs

Properties and efficient computation ††thanks: This work was supported by the
National Science Foundation NRT Award number 1633631, and also in part by the
United States Air Force under Contract No. FA8750-14-C-0011 under the DARPA
PPAML program; also by NIH grant R01 HD073179, and by the Leverhulme Trust,
and by the hospitality of the Sainsbury Laboratory Cambridge University.

Cory B. Scott  Eric Mjolsness

(Submitted September 9, 2019)

###### Abstract

We define a new family of similarity and distance measures on graphs, and
explore their theoretical properties in comparison to conventional distance
metrics. These measures are defined by the solution(s) to an optimization
problem which attempts find a map minimizing the discrepancy between two graph
Laplacian exponential matrices, under norm-preserving and sparsity
constraints. Variants of the distance metric are introduced to consider such
optimized maps under sparsity constraints as well as fixed time-scaling
between the two Laplacians. The objective function of this optimization is
multimodal and has discontinuous slope, and is hence difficult for univariate
optimizers to solve. We demonstrate a novel procedure for efficiently
calculating these optima for two of our distance measure variants. We present
numerical experiments demonstrating that (a) upper bounds of our distance
metrics can be used to distinguish between lineages of related graphs; (b) our
procedure is faster at finding the required optima, by as much as a factor of
$10^{3}$; and (c) the upper bounds satisfy the triangle inequality exactly
under some assumptions and approximately under others. We also derive an upper
bound for the distance between two graph products, in terms of the distance
between the two pairs of factors. Additionally, we present several possible
applications, including the construction of infinite “graph limits” by means
of Cauchy sequences of graphs related to one another by our distance measure.

###### Keywords:

Graph Distance Linear Assignment Problem Spectral Graph Theory Inexact Pattern
Matching

##  1 Introduction

Structure comparison, as well as structure summarization, is a ubiquitous
problem, appearing across multiple scientific disciplines. In particular, many
scientific problems (e.g. inference of molecular properties from structure,
pattern matching in data point clouds and scientific images) may be reduced to
the problem of inexact graph matching: given two graphs, compute a measure of
similarity that gainfully captures structural correspondence between the two.
Similarly, many algorithms for addressing multiple scales of dynamical
behavior rely on methods for automatically coarsening some model architecture.

In this work we present a graph distance metric, based on the Laplacian
exponential kernel of a graph. This measure generalizes the work of Hammond et
al. hammond2013graph  on graph diffusion distance for graphs of equal size;
crucially, our distance measure allows for graphs of inequal size. We
formulate the distance measure as the solution to an optimization problem
dependent on a comparison of the two graph Laplacians. This problem is a
nested optimization problem with the innermost layer consisting of
multivariate optimization subject to matrix constraints (e.g. orthogonality).
To compute this dissimilarity score efficiently, we also develop and
demonstrate the lower computational cost of an algorithm which calculates
upper bounds on the distance. This algorithm produces a
prolongation/restriction operator, $P$, which produces an optimally coarsened
version of the Laplacian matrix of a graph. Prolongation/restriction operators
produced via the method in this paper have previously been applied to
accelerate the training of machine learning algorithms in scott2018multilevel
.

###  1.1 Prior Work

Quantitative measures of similarity or dissimilarity between graphs have been
studied for decades owing to their relevance for problems in pattern
recognition including structure-based recognition of extended and compound
objects in computer vision, prediction of chemical similarity based on shared
molecular structure, and many other domains. Related problems arise in
quantitative modeling, for example in meshed discretizations of partial
differential equations and more recently in trainable statistical models of
data that feature graph-like models of connectivity such as Bayes Networks,
Markov Random Fields, and artificial neural networks. A core problem is to
define and compute how “similar” two graphs are in a way that is invariant to
a permutation of the the vertices of either graph, so that the answer doesn’t
depend on an arbitrary numbering of the vertices. On the other hand unlike an
arbitrary numbering, problem-derived semantic labels on graph vertices may
express real aspects of a problem domain and may be fair game for detecting
graph similarity. The most difficult case occurs when such labels are absent,
for example in an unstructured mesh, as we shall assume.

We mention just a few prior works to give an overview of the development of
graph distance measures over time, paying special attention to those which
share theoretical or algorithmic characteristics with the measure we
introduce. Our mathematical distinctions concern the following properties: (a)
Does the distance measure require an inner optimization loop? If so is it
mainly a discrete or continuous optimization formulation? (b) Does the
distance measure calculation naturally yield some kind of explicit map from
real-valued functions on vertices of one graph to functions on vertices of the
other? (A map from vertices to vertices would be a special case.) If we use
the term “graph signal” to mean a function $f:V(G_{1})\rightarrow S$ which
identifies each vertex of a graph $G_{1}$ with some state $s\in S$, then a
map-explicit graph distance is one whose calculation yields a second function
$g:V(G_{2})\rightarrow V(G_{1})$, with the composite function $f\circ
g:V(G_{2})\rightarrow S$. (c) Is the distance metric definable on the spectrum
of the graph alone, without regard to other data from the same graph? The
“spectrum” of a graph is a graph invariant calculated as the eigenvalues of a
matrix related to the adjacency matrix of the graph. Depending on context, the
spectrum can refer to eigenvalues of the adjacency matrix, graph Laplacian, or
normalized graph Laplacian of a graph. We will usually take the underlying
matrix to be the graph Laplacian, defined in detail in Section 1.3.
Alternatively, does it take into account more detailed “structural” aspects of
the graph? This categorization (structural vs. spectral) is similar to that
introduced in donnat2018tracking .

####  1.1.1 Graph-Edit Distance (structural, explicit, disc-opt):

The graph edit distance measures the total cost of converting one graph into
another with a sequence of local edit moves, with each type of move (vertex
deletion/addition, edge deletion/addition, edge division/contraction)
incurring a specified cost. Costs are chosen to suit the graph analysis
problem at hand; determining a cost assignment which makes the edit distance
most instructive for a certain set of graphs is an active area of research.
The distance measure is then the sum of these costs over an optimal sequence
of edits, which must be found using some optimization algorithm i.e. a
shortest-path algorithm (the best choice of algorithm may vary, depending on
how the costs are chosen). The sequence of edits may or may not (depending on
the exact set of allowable edit moves) be adaptable into an explicit map
between vertex-sets. Classic pattern recognition literature includes:
eshera1984graph  eshera1986image  gao2010survey  sanfeliu1983distance  .

####  1.1.2 Quadratic Point- and Graph- Matching (structural, explicit, cont-
opt)

Other work focuses on the construction of a point-to-point correspondence
between the vertices of two graphs. Gold et. al. gold1995new  define the
dissimilarity between two unlabelled weighted graphs with adjacency matrices
$A^{(1)}$ and $A^{(2)}$ as the solution to the following optimization problem
(for real-valued $M=[m_{ij}]$:

$\begin{array}[]{ll@{}ll}\text{minimize}&\displaystyle\sum\limits_{j=1}^{n_{2}}\sum\limits_{k=1}^{n_{1}}{\left(\sum\limits_{l=1}^{n_{2}}A^{(1)}_{jl}m_{lk}-\sum\limits_{p=1}^{n_{1}}m_{jp}A^{(2)}_{pk}\par\right)}^{2}&={\left|\left|A^{(1)}M-MA^{(2)}\right|\right|}_{F}^{2}\\\
\text{subject to}&\displaystyle\sum\limits_{i=1}^{n_{2}}m_{ij}=1,&j=1\ldots
n_{1}\\\ &\displaystyle\sum\limits_{j=1}^{n_{1}}m_{ij}=1,&i=1\ldots n_{2}\\\
&m_{ij}\geq 0&i=1\ldots n_{2}\\\ &&j=1\ldots n_{1}\\\ \end{array}$ (1)

where ${\left|\left|\cdot\right|\right|}_{F}^{2}$ is the squared Frobenius
norm. This problem is similar in structure to the optimization considered in
Section 4: a key difference being that Gold et al. consider optimization over
real-valued matchings between graph vertices, whereas we consider 0-1 valued
matchings between the eigenvalues of the graph Laplacians. In gold1995learning
and rangarajan1996novel  the authors present computational methods for
computing this optimum, and demonstrate applications of this distance measure
to various machine learning tasks such as 2D and 3D point matching, as well as
graph clustering. Gold et al. also introduce the _softassign_ , a method for
performing combinatorial optimization with both row and column constraints,
similar to those we consider.

####  1.1.3 Cut-Distance of Graphs (structural, implicit, disc-opt):

Lovász lovasz2012large  defines the _cut-distance_ of a pair of graphs as
follows: Let the $\Box$-norm of a matrix $B$ be given by:

$\displaystyle{\left|\left|B\right|\right|}_{\Box}=\frac{1}{n^{2}}\max_{S,T\subseteq
1\ldots n}\left|\sum_{i\in S,j\in T}B_{ij}\right|$ (2)

Given two labelled graphs $G_{1}$, $G_{2}$, on the same set of vertices, and
their adjacency matrices $A_{1}$ and $A_{2}$, the cut-distance
$d_{\text{cut}}(G_{1},G_{2})$ is then given by

$\displaystyle
D_{\text{cut}}(G_{1},G_{2})={\left|\left|A_{1}-A_{2}\right|\right|}_{\Box}$
(3)

(for more details, see lovasz2012large ). Computing this distance requires
combinatorial optimization (over all vertex subsets of $G_{1},G_{2}$) but this
optimization does not result in an explicit map between $G_{1}$ and $G_{2}$.

####  1.1.4 Wasserstein Earth-Mover Distance (spectral, implicit, disc-opt):

One common metric between graph spectra is the Wasserstein Earth Mover
Distance. Most generally, this distance measures the cost of transforming one
probability density function into another by moving mass under the curve. If
we consider the eigenvalues of a (possibly weighted) graph as point masses,
then the EMD measures the distance between the two spectra as the solution to
a transport problem (transporting one set of points to the other, subject to
constraints e.g. a limit on total distance travelled or a limit on the number
of ‘agents’ moving points). The EMD has been used in the past in various graph
clustering and pattern recognition contexts; see gu2015spectral . In the above
categorization, this is an optimization-based spectral distance measure, but
is implicit, since it does not produce a map from vertices of $G_{1}$ to those
of $G_{2}$ (informally, this is because the EMD is not translating one set of
eigenvalues into the other, but instead transforming their respective
histograms). Recent work applying the EMD to graph classification includes
dodonova2016kernel  and mcgregor2013sketching .

####  1.1.5 Diffusion Distance due to Hammond et al. hammond2013graph
(spectral, implicit, non-opt):

We discuss this recent distance metric more thoroughly below. This distance
measures the difference between two graphs as the maximum discrepancy between
probability distributions which represent single-particle diffusion beginning
from each of the nodes of $G_{1}$ and $G_{2}$. This distance is computed by
comparing the eigenvalues of the heat kernels of the two graphs. The
optimization involved in calculating this distance is a simple unimodal
optimization over a single scalar, $t$, representing the passage of time for
the diffusion process on the two graphs; hence we do not count this among the
“optimization based” methods we consider.

####  1.1.6 Novel Diffusion-Derived Measures (spectral, explicit, cont-opt):

In this work, we introduce a family of related graph distance measures which
compare two graphs in terms of similarity of a set of probability
distributions describing single-particle diffusion on each graph. For two
graphs $G_{1}$ and $G_{2}$ with respective Laplacians $L(G_{1})$ and
$L(G_{2})$, the matrices $e^{tL(G_{1})}$ and $e^{tL(G_{2})}$ are called the
_Laplacian Exponential Kernels_ of $G_{1}$ and $G_{2}$ ($t$ is a scalar
representing the passage of time). The column vectors of these matrices
describe the probability distribution of a single-particle diffusion process
starting from each vertex, after $t$ time has passed. The norm of the
difference of these two kernels thus describes how different these two graphs
are, from the perspective of single-particle diffusion, at time $t$. Since
these distributions are identical at very-early and very late times $t$ (we
formalize this notion in Section 2.1), a natural way to define a graph
distance is to take the supremum over all $t$. When the two graphs are the
same size, so are the two kernels, which may therefore be directly compared.
This case is the case considered by Hammond et al. hammond2013graph . However,
to compare two graphs of different sizes, we need a mapping between the column
vectors of $e^{tL(G_{1})}$ and $e^{tL(G_{2})}$. Optimization over a suitably
constrained prolongation/restriction operator between the graph Laplacians of
the two graphs is a permutation-invariant way to compare the behavior of a
diffusion process on each. The prolongation map $P$ thus calculated may then
be used to map signals (by which we mean values associated with vertices or
edges of a graph) on $G_{1}$ to the space of signals on $G_{2}$ (and vice
versa). In scott2018multilevel  we implicitly consider the weights of an
artificial neural network model to be graph signals, and use these operators
to train a hierarchy of linked neural network models. However, in that work we
do not address efficient calculation of this distance or the associated
operators, a major focus of this paper. In this work we also consider a time
conversion factor between diffusion on graphs of unequal size, and consider
the effect of limiting this optimization to sparse maps between the two graphs
(again, our case reduces to Hammond when the graphs in question are the same
size, dense $P$ and $R$ matrices are allowed, and our time-scaling parameter
is set to 1).

In this work, we present an algorithm for computing the type of nested
optimization given in our definition of distance (Equations 8 and 9). The
innermost loop of our distance measure optimization consists of a Linear
Assignment Problem (LAP, defined below) where the entries of the cost matrix
have a nonlinear dependence on some external variable. Our algorithm greatly
reduces both the count and size of calls to the external LAP solver. We use
this algorithm to compute an upper bound on our distance measure, but it could
also be useful in other similar nested optimization contexts: specifically,
nested optimization where the inner loop consists of a linear assignment
problem whose costs depend quadratically on the parameter in the outermost
loop.

###  1.2 Background

The ideal for a quantitative measure of similarity or distance on some set $S$
is usually taken to be a distance metric $d:S\times S\mapsto{\mathbb{R}}$
satisfying for all $x,y,z\in S$:

  * •

Non-negativity: $d(x,y)\geq 0$

  * •

Identity: $d(x,y)=0\iff x=y$

  * •

Symmetry: $d(x,y)=d(y,x)$

  * •

Triangle inequality: $d(x,z)\leq d(x,y)+d(y,z)$

Then $(S,d)$ is a metric space. Euclidean distance on ${\mathbb{R}}^{d}$ and
geodesic distance on manifolds satisfy these axioms. They can be used to
define algorithms that generalize from ${\mathbb{R}}^{d}$ to other spaces. A
variety of weakenings of these axioms are required in many applications, by
dropping some axioms and/or weakening others. For example if $S$ is a set of
nonempty sets of a metric space $S_{0}$, one can define the “Hausdorff
distance” on $S$ which is an extended pseudometric that obeys the triangle
inequality but not the Identity axiom and that can take values including
$+\infty$. As another example, any measure measure of distance on graphs which
is purely spectral (in the taxonomy of Section 1.1) cannot distinguish between
graphs which have identical spectra. We discuss this in more detail in Section
2.5.

Additional properties of distance metrics that generalize Euclidean distance
may pertain to metric spaces related by Cartesian product, for example, by
summing the squares of the distance metrics on the factor spaces. We will
consider an analog of this property in Section 3.3.

###  1.3 Definitions

  * •

Graph Laplacian: For an undirected graph $G$ with adjacency matrix $A$ and
vertex degrees $d_{1},d_{2}\ldots d_{n}$, we define the Laplacian of the graph
as

$\displaystyle L(G)=A-\textit{diag}(\\{d_{1},d_{2}\ldots
d_{n}\\})=A-\textit{diag}(\mathbf{1}\cdot A)=A(G)-D(G)$ (4)

$L(G)$ is sometimes instead defined as $D(G)-A(G)$; we take this sign
convention for $L(G)$ because it agrees with the standard continuum Laplacian
operator, $\Delta$, of a multivariate function $f$: $\Delta
f=\sum_{i=1}^{n}\frac{\delta^{2}f}{\delta x_{i}^{2}}$.

  * •

The squared Frobenius norm, ${\left|\left|A\right|\right|}^{2}_{F}$ of a
matrix $A$ is given by the sum of squares of matrix entries. This can
equivalently be written as $\operatorname{Tr}[A^{T}A]$.

  * •

Linear Assignment Problem (LAP): we take the usual definition of the Linear
Assignment Problem (see burkard1999linear , burkard2009assignment ): we have
two lists of items $S$ and $R$ (sometimes referred to as “workers” and
“jobs”), and a cost function $c:S\times R\rightarrow\mathbb{R}$ which maps
pairs of elements from $S$ and $R$ to an associated cost value. This can be
written as a linear program for real-valued $x_{ij}$ as follows:

$\begin{array}[]{ll@{}ll}\text{minimize}&\displaystyle\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}c(s_{i},r_{j})x_{ij}&\\\
\text{subject to}&\displaystyle\sum\limits_{i=1}^{m}x_{ij}\leq 1,&j=1\ldots
n\\\ &\displaystyle\sum\limits_{j=1}^{n}x_{ij}\leq 1,&i=1\ldots m\\\
&x_{ij}\geq 0&i=1\ldots m,j=1\ldots n\\\ \end{array}$ (5)

Generally, “Linear Assignment Problem” refers to the square version of the
problem where $|S|=|R|=n$, and the objective is to allocate the $n$ jobs to
$n$ workers such that each worker has exactly one job and vice versa. The case
where there are more workers than jobs, or vice versa, is referred to as a
Rectangular LAP or RLAP. In practice, the conceptually simplest method for
solving an RLAP is to convert it to a LAP by _augmenting_ the cost matrix with
several columns (rows) of zeros. In this case, solving the RLAP is equivalent
to solving a LAP with size $max(n,m)$. Other computational shortcuts exist;
see bijsterbosch2010solving  for details. Since the code we use to solve RLAPs
takes the augmented cost matrix approach, we do not consider other methods in
this paper.

  * •

Matching: we refer to a 0-1 matrix $M$ which is the solution of a particular
LAP as a “matching”. We may refer to the “pairs” or “points” of a matching, by
which we mean the pairs of indices $(i,j)$ with $M_{ij}=1$. We may also say in
this case that $M$ “assigns” $i$ to $j$.

  * •

Graph Lineage: we use the same definition of graph lineage as in
scott2018multilevel , namely: A graph lineage is a sequence of graphs, indexed
by $l\in\mathbb{N}=0,1,2,3\ldots$, satisfying the following:

    * –

$G_{0}$ is the graph with one vertex and one self-loop, and;

    * –

Successive members of the lineage grow roughly exponentially - that is, there
exists some base $b$ such that the growth rate as a function of level number
$l$ is $O(b^{l^{1+\epsilon}})$, for all $\epsilon>0$.

  * •

Box product $(\Box)$ of graphs: For two graphs $G$ and $H$ with vertex sets
$V(G)=\\{g_{1},g_{2}\ldots g_{n}\\}$ and $V(H)=\\{h_{1},h_{2}\ldots h_{m}\\}$,
we say the product graph $G\Box H$ is the graph with vertex set $V(G\Box
H)=V(G)\times V(H)$ and an adjacency relationship defined by:
$(g_{1},h_{1})\sim(g_{2},h_{2})$ in $G\Box H$ if and only if $g_{1}\sim g_{2}$
in $G$ and $h_{1}=h_{2}$, or $g_{1}=g_{2}$ and $h_{1}\sim h_{2}$ in $H$. Note
that the adjacency matrix of this relationship may be represented by the
following identity:

$\begin{split}A(G\Box H)=A(G)\otimes I_{m}+I_{n}\otimes A(H)\end{split}$ (6)

where $\otimes$ is the Kronecker Product of matrices (See hogben2006handbook ,
Section 11.4).

##  2 Graph Diffusion Distance Definitions

###  2.1 Diffusion Distance Definition

We generalize the diffusion distance defined by Hammond et al.
hammond2013graph . This distortion measure between two graphs $G_{1}$ and
$G_{2}$, of the same size, was defined as:

$\displaystyle
D_{\text{Hammond}}(G_{1},G_{2})=\sup_{t}{\left|\left|e^{tL_{1}}-e^{tL_{2}}\right|\right|}_{F}^{2}$
(7)

where $L_{i}$ represents the graph Laplacian of $G_{i}$. A typical example of
this distance, as $t$ is varied, can be seen in Figure 1.

Figure 1: Plot of graph distance between two graphs of the same size (a
$4\times 4$ grid and a path graph on 16 vertices), following the definition of
Hammond et al. hammond2013graph , as the parameter $t$ is varied. Note 1)
there is single global maximum of $1.592$ at $t=0.891$ and 2) similarity to
the plot of our distance measure as a function of $t$ given in Figure 2(a).

This may be interpreted as measuring the maximum divergence, as a function of
$t$, between diffusion processes starting from each vertex of each graph, as
measured by the squared Euclidean distance between the column vectors of
$e^{tL_{i}}$. Each column $v_{j}$ of $e^{tL_{i}}$ (which is called the
Laplacian Exponential Kernel) describes a probability distribution of visits
(by a random walk of duration $t$, with node transition probabilities given by
the columns of $e^{L}$) to the vertices of $G_{i}$, starting at vertex $j$.
This distance metric is then measuring the difference between the two graphs
by comparing these probability distributions; the motivation between taking
the supremum over all $t$ is that the value of the objective function at the
maximum is the most these two distributions can diverge. For further intuition
about why the peak is the most natural place to take as the distance, rather
than some other arbitrary time, note that at very early times and very late
times, the probability distribution of vertex visits is agnostic to graph
structure: at early times no diffusion has had a chance to take place, while
at very late times the distribution of vertex-visits converges to the uniform
distribution for each connected component of the graph. Hence we are most
interested in a regime of $t$-values in between these extremes, where
differences in $G_{1}$ and $G_{2}$ are apparent in their differing probability
distributions.

Our contribution generalizes this measure to allow for graphs of differing
size. We add two variables to this optimization: a _prolongation_ operator,
$P$ (represented as a rectangular matrix), and a time-scaling factor,
$\alpha$. The dissimilarity between two graphs $G_{1}$ and $G_{2}$ (with
Laplacians $L_{i}=L(G_{i})$) is then defined as:

$\displaystyle D^{2}(G_{1},G_{2})$
$\displaystyle=\sup_{t>0}\inf_{\alpha>0}\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{\frac{t}{\alpha}L_{1}}-e^{\alpha
tL_{2}}P\right|\right|}^{2}_{F}$ (8)

where $\mathcal{C}(P)$ represents some set of constraints on the matrix $P$.
It will be convenient for later calculations to introduce and assume the
concept of _transitive constraints_ \- by which we mean that for any
constraint $\mathcal{C}$, satisfaction of $\mathcal{C}$ by $P_{1}$ and $P_{2}$
implies satisfaction of $\mathcal{C}$ by their product $P_{1}P_{2}$ (when such
a product is defined). Some (non-exclusive) examples of transitive constraints
include orthogonality, particular forms of sparsity, and their conjunctions.

The simplest transitive constraint we will consider is that $P$ should be
orthogonal. Intuitively, an orthogonal $P$ represents a norm-preserving map
between nodes of $G_{1}$ and nodes of $G_{2}$, so we are measuring how well
diffusion on $G_{1}$ approximates diffusion on $G_{2}$, as projected by $P$.
Note that since in general $P$ is a rectangular matrix it is not necessarily
true that $PP^{T}=I$. We assume that $\left|G_{1}\right|=n_{1}\leq
n_{2}=\left|G_{2}\right|$; if not, the order of the operands is switched, so
that $P$ is always at least as wide as it is tall. We also briefly consider a
sparsity constraint in section 2.3 below. Since sparsity is more difficult to
treat numerically, our default constraint will be orthogonality alone. Other
constraints could include bandedness and other structural constraints (see
Section 7). We also note that because $L$ is finite-dimensional, the
exponential map is continuous and therefore we can swap the order of
optimization over $t$ and $\alpha$. The optimization procedure outlined in
this paper optimizes these variables in the order presented above (namely: an
outermost loop of maximization over $t$, a middle loop of minimization over
$\alpha$ given $t$, and an innermost loop of minimization over $P$ given $t$
and $\alpha$).

The other additional parameter, $\alpha$, controls dilation between the
passage of time in the two graphs, to account for different scales. Again, the
intuition is that we are interested in the difference between structural
properties of the graph (from the point of view of single-particle diffusion)
independent of the absolute number of nodes in the graph. As an example,
diffusion on an $n\times n$ grid is a reasonably accurate approximation of
more rapid diffusion on a $2n\times 2n$ grid, especially when $n$ is very
large. In our discussion of variants of this dissimilarity score, we will use
the notation $D^{2}(G_{1},G_{2}|x=c)$ to mean restrictions of any of our
distortion measure equations where variable $x$ is held to a constant value
$c$; In cases where it is clear from context which variable is held to a fixed
value $c$, we will write $D^{2}(G_{1},G_{2}|c)$

At very early times the second and higher-order terms of the Taylor Series
expansion of the matrix exponential function vanish, and so $e^{tL}\approx
I+tL$. This motivates the _early-time_ or “linear” version of this distance,
$\tilde{D}$:

$\displaystyle\tilde{D}^{2}(G_{1},G_{2})$
$\displaystyle=\inf_{\alpha>0}\inf_{P|\mathcal{C}(P)}{\left|\left|\frac{1}{\alpha}PL_{1}-\alpha
L_{2}P\right|\right|}^{2}_{F}$ (9)
$\displaystyle\approx\frac{1}{t^{2}}\left(\inf_{\alpha>0}\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{\frac{t}{\alpha}L_{1}}-e^{\alpha
tL_{2}}P\right|\right|}^{2}_{F}\right)$ (10)

(Note that the identity matrices cancel). The outermost optimization
(maximization over $t$) is removed for this version of the distance, as the
factor of $t$ can be pulled out:

$\displaystyle{\left|\left|\frac{t}{\alpha}PL_{1}-\alpha
tL_{2}P\right|\right|}^{2}_{F}=t^{2}{\left|\left|\frac{1}{\alpha}PL_{1}-\alpha
L_{2}P\right|\right|}^{2}_{F}$ (11)

For the exponential version of the dissimilarity score, we note briefly that
the supremum over $t$ of our objective function must exist, since for any
$G_{1},G_{2}$:

$\displaystyle D^{2}(G_{1},G_{2})\leq
D^{2}\left(G_{1},G_{2}\left|\alpha=1,P=\left[\begin{array}[]{c}I\\\
0\end{array}\right]\right.\right)$ (14)

In other words, the infimum over all $P$ and $\alpha$ is bounded above by any
particular choice of values for these variables. Since

$\displaystyle
D^{2}\left(G_{1},G_{2}\left|t=0,\alpha=1,P=\left[\begin{array}[]{c}I\\\
0\end{array}\right]\right.\right)=0,$ (17) and
$\displaystyle\lim_{t_{c}\rightarrow\infty}D^{2}\left(G_{1},G_{2}\left|t_{c},\alpha=1,P=\left[\begin{array}[]{c}I\\\
0\end{array}\right]\right.\right)=0$ (20)

this upper bound must have a supremum (possibly 0) at some
$t^{*}\in[0,\infty)$. Then

$\displaystyle
D^{2}\left(G_{1},G_{2}\left|t^{*},\alpha=1,P=\left[\begin{array}[]{c}I\\\
0\end{array}\right]\right.\right)$ (23)

must be finite and therefore so must the objective function.

###  2.2 Time-Scaled Graph Diffusion Distance

For any graphs $G_{1}$ and $G_{2}$, we define the Time-Scaled Graph Diffusion
Distance (TSGDD) as a scaled version of the linear distance given above, with
$\alpha$ fixed. Namely, let

$\displaystyle\tilde{D}^{2}_{r}(G_{1},G_{2})$
$\displaystyle={\left(n_{1}n_{2}\right)}^{-2r}\tilde{D}^{2}\left(G_{1},G_{2}\left|\alpha={\left(\frac{n_{1}}{n_{2}}\right)}^{r}\right.\right)$
(24)
$\displaystyle=\inf_{P|\mathcal{C}(P)}{\left(n_{1}n_{2}\right)}^{-2r}{\left|\left|{\left(\frac{n_{1}}{n_{2}}\right)}^{-r}PL_{1}-{\left(\frac{n_{1}}{n_{2}}\right)}^{r}L_{2}P\right|\right|}^{2}_{F}$

The intuition for this version of the distance measure is that we are
constraining the time dilation, $\alpha$, between $G_{1}$ and $G_{2}$ to be a
power of the ratio of the two graph sizes. The factor
${\left(n_{1}n_{2}\right)}^{-2r}$ is needed to ensure this version of the
distance satisfies the triangle inequality, as seen in Theorem 3.2.

###  2.3 Sparse-Diffusion Distance

We introduce the notation $\mathcal{C}(P)$ for a constraint predicate that
must be satisfied by prolongation matrix $P$, which is transitive in the sense
that:

$\displaystyle\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})\implies\mathcal{C}(P_{32}P_{21}).$
(25)

The simplest example is
$\mathcal{C}(P)=\mathcal{C}_{\text{orthog}}(P)\equiv(P^{T}P=I)$. However,
sparsity can be introduced in transitive form by
$\mathcal{C}(P)=\mathcal{C}_{\text{orthog}}(P)\wedge C_{\text{sparsity}}(P)$
where

$\mathcal{C}_{\text{sparsity}}(P)\equiv\Big{(}\max_{i,j}\text{degree}_{i,j}(P)\leq(n_{P\text{coarse}}/n_{P\text{fine}})^{s}\Big{)}$

for some real number $s\geq 0$. This predicate is transitive since

$\max_{i,j}\text{degree}_{i,j}(P_{32}P_{21})\leq\max_{i,j}\text{degree}_{i,j}(P_{32})\max_{i,j}\text{degree}_{i,j}(P_{21}),$

and since $n_{2}$ cancels out from the numerator and denominator of the
product of the fanout bounds. Here, $\text{degree}_{i,j}(M)$ is the total
number of nonzero entries in row $i$ or column $j$ of $M$.

This transitive sparsity constraint depends on a power-law parameter $s\geq
0$. When $s=0$, there is no sparsity constraint.

Another form of sparsity constraints are those which specify a pattern on
matrix entries which are allowed to be nonzero. Two simple examples (which are
also transitive) are matrices which are constrained to be upper triangular, as
well as matrices which are constrained to be of the form $A\otimes B$ where
$A$ and $B$ are themselves both constrained to be sparse. More complicated are
$n_{1}\times n_{2}$ matrices which are constrained to be banded for some
specified pattern of bands: more specifically, that there is a reordering of
the rows and columns that the number of diagonal bands (of width 1, slope
$\frac{n_{1}}{n_{2}}$) with nonzero entries is less than
${\left(\frac{n_{1}}{n_{2}}\right)}^{q}$ for some $0\leq q<1$. For example,
linear interpolation matrices between d-dimensional grids, with non-
overlapping source regions, follow this constraint.

As a final note on sparsity, we observe that any of the optimizations detailed
in this work could also be performed including a sparsity term (for example,
the $|\cdot|_{1}$-norm of the matrix $P$, calculated as
$\sum_{i}\sum_{j}|p_{ij}|$ is one possibility, as are terms which penalize $t$
or $\alpha$ far from 1), rather than explicit sparsity constraints. A
potential method of performing this optimization would be to start by
optimizing the non-sparse version of the objective function (as detailed in
Section 4.1) and then slowly increasing the strength of the regularization
term.

###  2.4 Directedness of Distance and Constraints

We note that this distance measure, as defined so far, is _directed_ : the
operands $G_{1}$ and $G_{2}$ serve differing roles in the objective function.
This additionally makes the constraint predicate $\mathcal{C}(P)$ ambiguous:
when we state that $\mathcal{C}$ represents orthogonality, it is not clear
whether we are referring to $P^{T}P=I$ or $PP^{T}=I$ (only one of which can be
true for a non-square matrix $P$). To remove this ambiguity, we will, for the
computations in the rest of this manuscript, define the distance metric to be
symmetric: the distance between $G_{1}$ and $G_{2}$ with
$\left|G_{1}\right|\leq\left|G_{2}\right|$ is always $D(G_{1},G_{2})$. $P$ is
then always at least as tall as it is wide, so of the two choices of
orthogonality constraint we select $P^{T}P=I$.

###  2.5 Variants of Distance Measure

Thus far we have avoided referring to this graph dissimilarity function as a
“distance metric”. As we shall see later, full optimization of Equations 8 and
9 over $\alpha$ and $P$ is too loose, in the sense that the distances
$D(G_{1},G_{2})$, $D(G_{2},G_{3})$, and $D(G_{1},G_{3})$ do not necessarily
satisfy the triangle inequality. The same is true for $\tilde{D}$. See Section
5.3.1 for numerical experiments suggesting a particular parameter regime where
the triangle inequality is satisfied. We thus define several
restricted/augmented versions of both $D$ and $\tilde{D}$ which are guaranteed
to satisfy the triangle inequality. These different versions are summarized in
Table 1. These variously satisfy some of the conditions necessary for
generalized versions of distance metrics, including:

  * •

Premetric: a function $d(x,y)$ for which $d(x,y)\geq 0$ and $d(x,y)=d(y,x)$
for all $x,y$.

  * •

Pseudometric: As a premetric, but additionally $d(x,z)\leq d(x,y)+d(y,z)$ for
all $x,y,z$.

  * •

$\rho$-inframetric: As a premetric, but additionally
$d(x,z)\leq\rho\left(d(x,y)+d(y,z)\right)$ and $d(x,y)=0$ if and only if
$x=y$, for all $x,y,z$.

Additionally, we note here that a distance measure on graphs using Laplacian
spectra can at best be a pseudometric, since isospectral, non-isomorphic
graphs are well-known to exist godsil1982constructing van2003graphs .
Characterizing the conditions under which two graphs are isospectral but not
isomorphic is an open problem in spectral graph theory. However, previous
computational work has led to the conjecture that “almost all” graphs are
uniquely defined by their spectra brouwer2011spectra brouwer2009cospectral
van2009developments , in the sense that the probability of two graphs of size
$n$ being isospectral but not isomorphic goes to 0 as $n\rightarrow\infty$.
Furthermore, our numerical experiments have indicated that the violation of
the triangle inequality is bounded, in the sense that
$D(G_{1},G_{3})\leq\rho*\left(D(G_{1},G_{2})+D(G_{2},G_{3})\right)$ for
$\rho\approx 2.1$. This means that even in the least restricted case our
similarity measure may be a $2.1$-infra-pseudometric on graphs (and, since
such isospectral, non-isomorphic graphs are relatively rare, it behaves like a
$2.1-$inframetric). As we will see in Section 3, in some more restricted cases
we can prove triangle inequalities, making our measure a pseudometric. In
Section 4.1, we discuss an algorithm for computing the optima in Equations (8)
and (9). First, we discuss some theoretical properties of this dissimilarity
measure.

$t$ |  $\alpha$ |  $s$ |  Classification |  Treatment in this manuscript
---|---|---|---|---
Fixed at $t_{c}<\epsilon$ |  Fixed at $\alpha_{c}=1$ |  $s=0$ |  Pseudometric |  Defined in Equation 28. Optimized by one pass of LAP solver. Triangle inequality proven in Theorem 3.1.
Fixed at $t_{c}<\epsilon$ |  Fixed at $\alpha_{c}=(\frac{n_{1}}{n_{2}})^{r}$ |  $s=0$ |  Pseudometric |  Defined in Equation (24). Optimized by one pass of LAP solver. Triangle inequality proven in Theorem 3.2.
Fixed at $t_{c}<\epsilon$ |  Optimized |  $s=0$ |  Premetric |  Defined in Equation 9. Optimized by Algorithm 1. Triangle inequality violations examined experimentally in Section 5.3.1.
Optimized |  Fixed at $\alpha_{c}=1$ |  $s=0$ |  |  When $|G_{1}|=|G_{2}|$, this is Hammond et. al’s version of graph distance.
Optimized |  Fixed at $\alpha_{c}=(\frac{n_{1}}{n_{2}})^{r}$ |  $s=0$ |  |
Optimized |  Optimized |  $s=0$ |  Premetric |  Defined in Equation 8. Optimized by Algorithm 2. Graph Product upper bound proven in Theorem 3.3. Triangle inequality violations examined experimentally in Section 5.3.1. Used to calculate graph distances in Sections 5.3.2 and 5.3.3.
Fixed at $t_{c}<\epsilon$ |  Fixed at $\alpha_{c}=1$ |  $s>0$ |  Pseudometric |  Triangle inequality proven in Theorem 3.1.
Fixed at $t_{c}<\epsilon$ |  Fixed at $\alpha_{c}=(\frac{n_{1}}{n_{2}})^{r}$ |  $s>0$ |  Pseudometric |  Triangle inequality proven in Theorem 3.2.
Fixed at $t_{c}<\epsilon$ |  Optimized |  $s>0$ |  |
Optimized |  Fixed at $\alpha_{c}=1$ |  $s>0$ |  |
Optimized |  Fixed at $\alpha_{c}=(\frac{n_{1}}{n_{2}})^{r}$ |  $s>0$ |  |
Optimized |  Optimized |  $s>0$ |  |  Discussed in Section 2.3.
Table 1: Summary of this paper’s investigation of different forms of our graph
dissimilarity measure. In this work, we systematically explore properties of
this measure given sparsity parameter $s=0$, and various regimes of $t$ (fixed
at some early time, or maximized over all $t$) and $\alpha$ (fixed at
$\alpha=1$, fixed at a constant power $r$ of the ratio of graph sizes, or
minimized over all $\alpha$. We leave exploration of nonzero values of the
sparsity parameter to future work.

##  3 Theoretical Properties of D(G1, G2)

In this section we prove several properties of various instances of our graph
dissimilarity score, including triangle inequalities for some specific
versions and an upper bound on the distance between two graph products.

###  3.1 Properties when $t$ is constant, $\alpha=1$

####  3.1.1 Triangle Inequality

######  Lemma 1.

For any matrices $M$ and $P$, with $P$ satisfying $P^{T}P=I$,
${\left|\left|PM\right|\right|}^{2}_{F}\leq{\left|\left|M\right|\right|}^{2}_{F}$
and
${\left|\left|MP\right|\right|}^{2}_{F}\leq{\left|\left|M\right|\right|}^{2}_{F}$
.

###### Proof.

Suppose without loss of generality that $P^{T}P=I$. Then:

  1. 1.

${\left|\left|PM\right|\right|}^{2}_{F}=\operatorname{Tr}[M^{T}P^{T}PM]=\operatorname{Tr}[M^{T}M]={\left|\left|M\right|\right|}^{2}_{F}$

  2. 2.

If $P^{T}P=I$, then letting $PP^{T}=\Pi$, $\Pi$ is a projection operator
satisfying $\Pi^{T}=\Pi=\Pi^{2}$. Then,

$\displaystyle{\left|\left|MP\right|\right|}^{2}_{F}$
$\displaystyle=\operatorname{Tr}[P^{T}M^{T}MP]$
$\displaystyle=\operatorname{Tr}[M^{T}MPP^{T}]=\operatorname{Tr}[M^{T}M\Pi]$

$\begin{split}{\left|\left|M\right|\right|}^{2}_{F}=\operatorname{Tr}[M^{T}M]&=\operatorname{Tr}[M^{T}M(\Pi+(I-\Pi))]\\\
&=\operatorname{Tr}[M^{T}M\Pi]+\operatorname{Tr}[M^{T}M(I-\Pi)]\\\
&=\operatorname{Tr}[M^{T}M\Pi]+\operatorname{Tr}[M^{T}M{(I-\Pi)}^{2}]\\\
&={\left|\left|MP\right|\right|}_{F}^{2}+{\left|\left|M(I-\Pi)\right|\right|}_{F}^{2}\\\
&\geq{\left|\left|MP\right|\right|}_{F}^{2}\end{split}$ (26)

∎

######  Theorem 3.1.

$\tilde{D}^{2}$ satisfies the triangle inequality for $\alpha=1$.

###### Proof.

Let $G_{1},G_{2},G_{3}$ be simple graphs, with Laplacians $L_{1},L_{2},L_{3}$.
Let

$\displaystyle
P_{31}=\arg\inf_{P|\mathcal{C}(P)}{\left|\left|PL_{1}-L_{3}P\right|\right|}_{F}^{2}$
(27)

Then

$\begin{split}\tilde{D}^{2}(G_{1},G_{3}\left|\alpha=1\right.)&={\left|\left|P_{31}L_{1}-L_{3}P_{31}\right|\right|}_{F}^{2}=\inf_{P|\mathcal{C}(P)}{\left|\left|PL_{1}-L_{3}P\right|\right|}_{F}^{2}\\\
&\leq\inf_{P_{32},P_{21}|\mathcal{C}(P_{32}P_{21})}{\left|\left|P_{32}P_{21}L_{1}-L_{3}P_{32}P_{21}\right|\right|}_{F}^{2}\\\
\end{split}$ (28)

where we write $\mathcal{C}(P_{32}P_{21})$ to signify that the product
$P_{32}P_{21}$ satisfies the original transitive constraints on $P$, e.g.
orthogonality and/or sparsity. Since the constraint predicate $\mathcal{C}(P)$
satisfies Equation (25), then so does their product, so we may write

$\begin{split}\tilde{D}(G_{1},G_{3}\left|\alpha=1\right.)&\leq\inf_{P_{32}|\mathcal{C}(P_{32})}\inf_{P_{21}|\mathcal{C}(P_{21})}{\left|\left|P_{32}P_{21}L_{1}-L_{3}P_{32}P_{21}\right|\right|}_{F}\\\
&=\inf_{P_{32}|\mathcal{C}(P_{32})}\inf_{P_{21}|\mathcal{C}(P_{21})}{\left|\left|P_{32}P_{21}L_{1}-P_{32}L_{2}P_{21}\right.\right.}\\\
&{\left.\left.\quad\quad\qquad+\quad
P_{32}L_{2}P_{21}-L_{3}P_{32}P_{21}\right|\right|}_{F}\\\
&\leq\inf_{P_{32}|\mathcal{C}(P_{32})}\inf_{P_{21}|\mathcal{C}(P_{21})}\left({\left|\left|P_{32}P_{21}L_{1}-P_{32}L_{2}P_{21}\right|\right|}_{F}\right.\\\
&\quad\quad\qquad+\quad\left.{\left|\left|P_{32}L_{2}P_{21}-L_{3}P_{32}P_{21}\right|\right|}_{F}\right)\\\
&=\inf_{P_{32}|\mathcal{C}(P_{32})}\inf_{P_{21}|\mathcal{C}(P_{21})}\left({\left|\left|P_{32}\left(P_{21}L_{1}-L_{2}P_{21}\right)\right|\right|}_{F}\right.\\\
&\quad\quad\qquad+\quad\left.{\left|\left|\left(P_{32}L_{2}-L_{3}P_{32}\right)P_{21}\right|\right|}_{F}\right)\\\
\end{split}$ (29)

By Lemma 1,

$\begin{split}\tilde{D}(G_{1},G_{3}\left|\alpha=1\right.)&\leq\inf_{P_{32}|\mathcal{C}(P_{32})}\inf_{P_{21}|\mathcal{C}(P_{21})}\left({\left|\left|P_{21}L_{1}-L_{2}P_{21}\right|\right|}_{F}\right.\\\
&\quad\quad\qquad+\quad\left.{\left|\left|P_{32}L_{2}-L_{3}P_{32}\right|\right|}_{F}\right)\\\
&=\inf_{P_{21}|\mathcal{C}(P_{21})}{\left|\left|P_{21}L_{1}-L_{2}P_{21}\right|\right|}_{F}\\\
&\quad\quad\qquad+\quad\inf_{P_{32}|\mathcal{C}(P_{32})}{\left|\left|P_{32}L_{2}-L_{3}P_{32}\right|\right|}_{F}\\\
&=\tilde{D}(G_{1},G_{2}\left|\alpha=1\right.)+\tilde{D}(G_{2},G_{3}\left|\alpha=1\right.)\end{split}$
(30)

∎

We note that in this proof we use $L_{1},L_{2},$ and $L_{3}$ (making this the
small-$t$ or linear version of the objective function), but the same argument
holds when all three are replaced with $e^{tL_{i}}$, so we also have

######  Corollary 1.

$D$ satisfies the triangle inequality for $\alpha=1$.

###### Proof.

By the same calculation as in Theorem 3.1, with all $L_{i}$ replaced by
$e^{t_{c}L_{i}}$, we have

$\displaystyle D\left(\left.G_{1},G_{3}\right|t_{c},\alpha=1\right)$
$\displaystyle\leq
D(G_{1},G_{2}\left|t_{c},\alpha=1\right)+D(G_{2},G_{3}\left|t_{c},\alpha=1\right)$
(31)

for any constant $t_{c}$. Then, letting

$\displaystyle
t_{13}=\arg\sup_{t_{c}}D\left(\left.G_{1},G_{3}\right|t_{c},\alpha=1\right)$
(32)

we have:

$\begin{split}D\left(\left.G_{1},G_{3}\right|\alpha=1\right)&=\sup_{t_{c}}D\left(\left.G_{1},G_{3}\right|t_{c},\alpha=1\right)\\\
&=D\left(\left.G_{1},G_{3}\right|t_{13},\alpha=1\right)\\\ &\leq
D(G_{1},G_{2}\left|t_{13},\alpha=1\right)+D(G_{2},G_{3}\left|t_{13},\alpha=1\right)\\\
&\leq\sup_{t_{c}}D\left(\left.G_{1},G_{2}\right|t_{c},\alpha=1\right)\\\
&\quad\quad\qquad+\quad\sup_{t_{c}}D\left(\left.G_{2},G_{3}\right|t_{c},\alpha=1\right)\\\
&=D\left(\left.G_{1},G_{2}\right|\alpha=1\right)+D\left(\left.G_{2},G_{3}\right|\alpha=1\right)\end{split}$
(33)

∎

###  3.2 Properties when $t$ is constant,
$\alpha={\left(\frac{n_{1}}{n_{2}}\right)}^{r}$

Recall our definition of the _Time-Scaled Graph Diffusion Distance_ from
Equation (24), for $\alpha={\left(\frac{n_{1}}{n_{2}}\right)}^{r}$,
$r\in\mathbb{R}$:

$\begin{split}\tilde{D}^{2}_{r}(G_{1},G_{2})&={\left(n_{1}n_{2}\right)}^{-2r}\tilde{D}^{2}\left(G_{1},G_{2}\left|\alpha={\left(\frac{n_{1}}{n_{2}}\right)}^{r}\right.\right)\\\
&=\inf_{P|\mathcal{C}(P)}{\left(n_{1}n_{2}\right)}^{-2r}{\left|\left|{\left(\frac{n_{1}}{n_{2}}\right)}^{-r}PL_{1}-{\left(\frac{n_{1}}{n_{2}}\right)}^{r}L_{2}P\right|\right|}_{F}^{2}\end{split}$
(34)

######  Theorem 3.2.

The TSGDD, as defined above, satisfies the triangle inequality.

###### Proof.

As above, let $G_{1},G_{2},G_{3}$ be three graphs with
$n_{i}=\left|G_{i}\right|$ and $n_{1}\leq n_{2}\leq n_{3}$, and let $L_{i}$ be
the Laplacian of $G_{i}$. Let $\mathcal{C}(P)$ represent a transitive
constraint predicate, also as described previously. Then, for a constant
$r\in\mathbb{R}$, we have:

$\displaystyle\tilde{D}_{r}(G_{1},G_{3})=$
$\displaystyle\qquad\inf_{P|\mathcal{C}(P)}{\left(n_{1}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{3}}\right)}^{-r}PL_{1}-{\left(\frac{n_{1}}{n_{3}}\right)}^{r}L_{3}P\right|\right|}_{F}$
$\displaystyle\qquad\leq\inf_{P_{32},P_{21}|\mathcal{C}(P_{32}P_{21})}{\left(n_{1}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{3}}\right)}^{-r}P_{32}P_{21}L_{1}-{\left(\frac{n_{1}}{n_{3}}\right)}^{r}L_{3}P_{32}P_{21}\right|\right|}_{F}$
under the assumption, as in Equation (25), that
$\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})\implies\mathcal{C}(P_{32}P_{21})$,
$\displaystyle\tilde{D}_{r}(G_{1},G_{3})\leq$
$\displaystyle\inf_{P_{32},P_{21}|\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})}{\left(n_{1}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{3}}\right)}^{-r}P_{32}P_{21}L_{1}-{\left(\frac{n_{1}}{n_{3}}\right)}^{r}L_{3}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle=\inf_{P_{32},P_{21}|\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})}{\left(n_{1}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{3}}\right)}^{-r}P_{32}P_{21}L_{1}-{\left(\frac{n_{1}n_{3}}{n_{2}^{2}}\right)}^{r}P_{32}L_{2}P_{21}\right.\right.}$
$\displaystyle\hbox{}\qquad+\qquad{\left.\left.{\left(\frac{n_{1}n_{3}}{n_{2}^{2}}\right)}^{r}P_{32}L_{2}P_{21}-{\left(\frac{n_{1}}{n_{3}}\right)}^{r}L_{3}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle\leq\inf_{P_{32},P_{21}|\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})}{\left(n_{1}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{3}}\right)}^{-r}P_{32}P_{21}L_{1}-{\left(\frac{n_{1}n_{3}}{n_{2}^{2}}\right)}^{r}P_{32}L_{2}P_{21}\right|\right|}_{F}$
$\displaystyle\hbox{}\qquad+\qquad{\left(n_{1}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}n_{3}}{n_{2}^{2}}\right)}^{r}P_{32}L_{2}P_{21}-{\left(\frac{n_{1}}{n_{3}}\right)}^{r}L_{3}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle=\inf_{P_{32},P_{21}|\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})}{\left(n_{1}n_{3}\right)}^{-r}{\left(\frac{n_{3}}{n_{2}}\right)}^{r}{\left|\left|{\left(\frac{n_{1}}{n_{2}}\right)}^{-r}P_{32}P_{21}L_{1}-{\left(\frac{n_{1}}{n_{2}}\right)}^{r}P_{32}L_{2}P_{21}\right|\right|}_{F}$
$\displaystyle\hbox{}\qquad+\qquad{\left(n_{1}n_{3}\right)}^{-r}{\left(\frac{n_{1}}{n_{2}}\right)}^{r}{\left|\left|{\left(\frac{n_{2}}{n_{3}}\right)}^{-r}P_{32}L_{2}P_{21}-{\left(\frac{n_{2}}{n_{3}}\right)}^{r}L_{3}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle=\inf_{P_{32},P_{21}|\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})}{\left(n_{1}n_{2}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{2}}\right)}^{-r}P_{32}P_{21}L_{1}-{\left(\frac{n_{1}}{n_{2}}\right)}^{r}P_{32}L_{2}P_{21}\right|\right|}_{F}$
$\displaystyle\hbox{}\qquad+\qquad{\left(n_{2}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{2}}{n_{3}}\right)}^{-r}P_{32}L_{2}P_{21}-{\left(\frac{n_{2}}{n_{3}}\right)}^{r}L_{3}P_{32}P_{21}\right|\right|}_{F}$
By Lemma 1, $\displaystyle\tilde{D}_{r}(G_{1},G_{3})$
$\displaystyle\leq\inf_{P_{32},P_{21}|\mathcal{C}(P_{32})\wedge\mathcal{C}(P_{21})}{\left(n_{1}n_{2}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{2}}\right)}^{-r}P_{21}L_{1}-{\left(\frac{n_{1}}{n_{2}}\right)}^{r}L_{2}P_{21}\right|\right|}_{F}$
$\displaystyle\hbox{}\qquad+\qquad{\left(n_{2}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{2}}{n_{3}}\right)}^{-r}P_{32}L_{2}-{\left(\frac{n_{2}}{n_{3}}\right)}^{r}L_{3}P_{32}\right|\right|}_{F}$
$\displaystyle=\inf_{P_{21}|\mathcal{C}(P_{21})}{\left(n_{1}n_{2}\right)}^{-r}{\left|\left|{\left(\frac{n_{1}}{n_{2}}\right)}^{-r}P_{21}L_{1}-{\left(\frac{n_{1}}{n_{2}}\right)}^{r}L_{2}P_{21}\right|\right|}_{F}$
$\displaystyle\hbox{}\qquad+\qquad\inf_{P_{32}|\mathcal{C}(P_{32})}{\left(n_{2}n_{3}\right)}^{-r}{\left|\left|{\left(\frac{n_{2}}{n_{3}}\right)}^{-r}P_{32}L_{2}-{\left(\frac{n_{2}}{n_{3}}\right)}^{r}L_{3}P_{32}\right|\right|}_{F}$
$\displaystyle=\tilde{D}_{r}(G_{1},G_{2})+\tilde{D}_{r}(G_{2},G_{3})$ and so
$\displaystyle\tilde{D}_{r}(G_{1},G_{3})\leq\tilde{D}_{r}(G_{1},G_{2})+\tilde{D}_{r}(G_{2},G_{3})$

for any fixed $r\in\mathbb{R}$. ∎

Note that in the proofs of Theorem 3.1, Theorem 3.2, and Corollary 1, we
assume that the constraint predicate $\mathcal{C}(P)$ includes at least
orthogonality (so that we may apply Lemma 1). However, this constraint
predicate could be more strict, e.g. include both orthogonality and sparsity.
Hence these statements also apply to the $s>0$ cases in Table 1, which we do
not otherwise consider in this work: in our numerical experiments we (for
reasons of computational simplicity) only require our optimization over $P$ be
orthogonally constrained.

###  3.3 Upper Bounds for Graph Products

We now consider the case where we want to compute the distance of two graph
box products, i.e. $D\left(\mathbf{G}_{1},\mathbf{G}_{2}\right)$ where

$\displaystyle\mathbf{G}_{1}=G_{1}^{(1)}\Box
G_{1}^{(2)}\quad\text{and}\quad\mathbf{G}_{2}=G_{2}^{(1)}\Box G_{2}^{(2)}$
(35)

and

$\begin{split}P^{(1)}=\arg\inf_{P_{c}\left|\mathcal{C}(P_{c})\right.}D\left(G_{1}^{(1)},G_{2}^{(1)}|t_{c},\alpha_{c},P_{c}\right)\\\
P^{(2)}=\arg\inf_{P_{c}\left|\mathcal{C}(P_{c})\right.}D\left(G_{1}^{(2)},G_{2}^{(2)}|t_{c},\alpha_{c},P_{c}\right)\end{split}$
(36)

are known for some $t_{c},\alpha_{c}$. Previous work scott2018multilevel
included a proof of a similar inequality for the small-$t$ (“linear”) case of
our objective function.

######  Theorem 3.3.

Let $\mathbf{G}_{1}$ and $\mathbf{G}_{2}$ be graph box products as described
above, and for a graph $G$ let $L(G)$ be its Laplacian. For fixed $t=t_{c}$,
$\alpha=\alpha_{c}$, $P^{(i)}$ as given above, for any $\lambda\in[0,1]$, we
have

$\displaystyle\inf_{P_{c}\left|\mathcal{C}(P_{c})\right.}D\left(\mathbf{G}_{1},\mathbf{G}_{2}\right)\quad\leq\quad$
(37)
$\displaystyle\qquad\lambda\left({\left|\left|e^{\frac{t_{c}}{\alpha_{c}}L(G^{(2)}_{1})}\right|\right|}_{F}+{\left|\left|e^{t_{c}\alpha_{c}L(G^{(2)}_{2})}\right|\right|}_{F}\right)D\left(G_{1}^{(1)},G_{2}^{(1)}|P^{(1)}\right)$
$\displaystyle\qquad\qquad+\left(1-\lambda\right)\left({\left|\left|e^{\frac{t_{c}}{\alpha_{c}}L(G^{(1)}_{1})}\right|\right|}_{F}+{\left|\left|e^{t_{c}\alpha_{c}L(G^{(1)}_{2})}\right|\right|}_{F}\right)D\left(G_{1}^{(2)},G_{2}^{(2)}|P^{(2)}\right)$

where all distances are evaluated at $t=t_{c}$, $\alpha=\alpha_{c}$, but we
have omitted those terms for simplicity of notation.

###### Proof.

For graph products $\mathbf{G}_{i}$, we have

$\displaystyle L(\mathbf{G}_{i})$ $\displaystyle=L(G^{(1)}_{i})\oplus
L(G^{(2)}_{i})$ (38) $\displaystyle=\left(L(G^{(1)}_{i})\otimes
I_{\left|L(G^{(2)}_{i})\right|}\right)+\left(I_{\left|L(G^{(1)}_{i})\right|}\otimes
L(G^{(2)}_{i})\right)$

(this fact can be easily verified from the formula for the adjacency matrix of
a graph box product, given in the definition in Section 1.3), and so

$\displaystyle\exp\left[{cL(\mathbf{G}_{i})}\right]$
$\displaystyle=\exp\left[c{\left(L(G^{(1)}_{i})\otimes
I_{\left|L(G^{(2)}_{i})\right|}\right)+\left(I_{\left|L(G^{(1)}_{i})\right|}\otimes
L(G^{(2)}_{i})\right)}\right].$ (39)

Because $A\otimes I_{|B|}$ and $I_{|A|}\otimes B$ commute for any $A$ and $B$,

$\displaystyle\exp\left[{cL(\mathbf{G}_{i})}\right]$
$\displaystyle=\exp\left[c{\left(L(G^{(1)}_{i})\otimes
I_{\left|L(G^{(2)}_{i})\right|}\right)}\right]\exp\left[c{\left(I_{\left|L(G^{(1)}_{i})\right|}\otimes
L(G^{(2)}_{i})\right)}\right]$ (40)
$\displaystyle=\left(\exp\left[c{L(G^{(1)}_{i})}\right]\otimes
I_{\left|L(G^{(2)}_{i})\right|}\right)\left(I_{\left|L(G^{(1)}_{i})\right|}\otimes\exp\left[c{L(G^{(2)}_{i})}\right]\right)$
$\displaystyle=\exp\left[c{L(G^{(1)}_{i})}\right]\otimes\exp\left[c{L(G^{(2)}_{i})}\right]\
$

We will make the following abbreviations:

$\begin{array}[]{ccc}\mathbf{E}_{1}=e^{\frac{t_{c}}{\alpha_{c}}L(\mathbf{G}_{1})}&E^{(1)}_{1}=e^{\frac{t_{c}}{\alpha_{c}}L(G^{(1)}_{1})}&E^{(2)}_{1}=e^{\frac{t_{c}}{\alpha_{c}}L(G^{(2)}_{1})}\\\
\mathbf{E}_{2}=e^{t_{c}\alpha_{c}L(\mathbf{G}_{2})}&E^{(1)}_{2}=e^{t_{c}\alpha_{c}L(G^{(1)}_{2})}&E^{(2)}_{2}=e^{t_{c}\alpha_{c}L(G^{(2)}_{2})}\end{array}$

Then,

$\displaystyle\inf_{P|\mathcal{C}(P)}D\left(\mathbf{G}_{1},\mathbf{G}_{2}\right)$
$\displaystyle\leq D\left(\mathbf{G}_{1},\mathbf{G}_{2}|P^{(1)}\otimes
P^{(2)}\right)$ (41) $\displaystyle={\left|\left|\left(P^{(1)}\otimes
P^{(2)}\right)\mathbf{E}_{1}-\mathbf{E}_{2}\left(P^{(1)}\otimes
P^{(2)}\right)\right|\right|}_{F}$

$\displaystyle={\left|\left|\left(P^{(1)}\otimes
P^{(2)}\right)\left(E^{(1)}_{1}\otimes
E^{(2)}_{1}\right)-\left(E^{(1)}_{2}\otimes
E^{(2)}_{2}\right)\left(P^{(1)}\otimes P^{(2)}\right)\right|\right|}_{F}$
$\displaystyle={\left|\left|\left(P^{(1)}E^{(1)}_{1}\otimes
P^{(2)}E^{(2)}_{1}\right)-\left(E^{(1)}_{2}P^{(1)}\otimes
E^{(2)}_{2}P^{(2)}\right)\right|\right|}_{F}^{2}$
$\displaystyle={\left|\left|\left(P^{(1)}E^{(1)}_{1}\otimes
P^{(2)}E^{(2)}_{1}\right)-\left(P^{(1)}E^{(1)}_{1}\otimes
E^{(2)}_{2}P^{(2)}\right)\right.\right.}$ (42)
$\displaystyle\quad+{\left.\left.\left(P^{(1)}E^{(1)}_{1}\otimes
E^{(2)}_{2}P^{(2)}\right)-\left(E^{(1)}_{2}P^{(1)}\otimes
E^{(2)}_{2}P^{(2)}\right)\right|\right|}_{F}$
$\displaystyle\leq{\left|\left|\left(P^{(1)}E^{(1)}_{1}\otimes
P^{(2)}E^{(2)}_{1}\right)-\left(P^{(1)}E^{(1)}_{1}\otimes
E^{(2)}_{2}P^{(2)}\right)\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|\left(P^{(1)}E^{(1)}_{1}\otimes
E^{(2)}_{2}P^{(2)}\right)-\left(E^{(1)}_{2}P^{(1)}\otimes
E^{(2)}_{2}P^{(2)}\right)\right|\right|}_{F}$
$\displaystyle={\left|\left|P^{(1)}E^{(1)}_{1}\otimes\left(P^{(2)}E^{(2)}_{1}-E^{(2)}_{2}P^{(2)}\right)\right|\right|}_{F}$
(43)
$\displaystyle\quad\quad+\quad{\left|\left|\left(P^{(1)}E^{(1)}_{1}-E^{(1)}_{2}P^{(1)}\right)\otimes
E^{(2)}_{2}P^{(2)}\right|\right|}_{F}$
$\displaystyle={\left|\left|P^{(1)}E^{(1)}_{1}\right|\right|}_{F}{\left|\left|P^{(2)}E^{(2)}_{1}-E^{(2)}_{2}P^{(2)}\right|\right|}_{F}$
(44)
$\displaystyle\quad\quad+\quad{\left|\left|P^{(1)}E^{(1)}_{1}-E^{(1)}_{2}P^{(1)}\right|\right|}_{F}{\left|\left|E^{(2)}_{2}P^{(2)}\right|\right|}_{F}.$

By Lemma 1,

$\begin{split}\inf_{P|\mathcal{C}(P)}D\left(\mathbf{G}_{1},\mathbf{G}_{2}\right)&\leq{\left|\left|E^{(1)}_{1}\right|\right|}_{F}{\left|\left|P^{(2)}E^{(2)}_{1}-E^{(2)}_{2}P^{(2)}\right|\right|}_{F}\\\
&\quad\quad+\quad{\left|\left|P^{(1)}E^{(1)}_{1}-E^{(1)}_{2}P^{(1)}\right|\right|}_{F}{\left|\left|E^{(2)}_{2}\right|\right|}_{F}.\end{split}$
(45) If we instead use $\left(E^{(1)}_{2}P^{(1)}\otimes
P^{(2)}E^{(2)}_{1}\right)$ as the cross term in Equation (42), we have
$\displaystyle\inf_{P}D\left(\mathbf{G}_{1},\mathbf{G}_{2}\right)\leq{\left|\left|E^{(1)}_{2}\right|\right|}_{F}{\left|\left|P^{(2)}E^{(2)}_{1}-E^{(2)}_{2}P^{(2)}\right|\right|}_{F}$
(47)
$\displaystyle\quad\quad+\quad{\left|\left|P^{(1)}E^{(1)}_{1}-E^{(1)}_{2}P^{(1)}\right|\right|}_{F}{\left|\left|E^{(2)}_{1}\right|\right|}_{F}$

A linear combination of these two bounds gives us the desired bound. ∎

This has the additional consequence that

$\displaystyle\inf_{P_{c}\left|\mathcal{C}(P_{c})\right.}D\left(\mathbf{G}_{1},\mathbf{G}_{2}\right)\quad\leq$
(48)
$\displaystyle\min\left[\left({\left|\left|e^{\frac{t_{c}}{\alpha_{c}}L(G^{(2)}_{1})}\right|\right|}_{F}+{\left|\left|e^{t_{c}\alpha_{c}L(G^{(2)}_{2})}\right|\right|}_{F}\right)D\left(G_{1}^{(1)},G_{2}^{(1)}|P^{(1)}\right),\right.$
$\displaystyle\quad\quad\left.\left({\left|\left|e^{\frac{t_{c}}{\alpha_{c}}L(G^{(1)}_{1})}\right|\right|}_{F}+{\left|\left|e^{t_{c}\alpha_{c}L(G^{(1)}_{2})}\right|\right|}_{F}\right)D\left(G_{1}^{(2)},G_{2}^{(2)}|P^{(2)}\right)\right]$

Additionally, if

$\displaystyle E^{(1)}_{i}=E^{(2)}_{i}\text{ for }i\in 1,2\quad\text{and}\quad
P^{(1)}=P^{(2)},$ (49)

This reduces further to

$\displaystyle D\left(\mathbf{G}_{1},\mathbf{G}_{2}|P^{(1)}\otimes
P^{(1)}\right)$
$\displaystyle\leq\min\left({\left|\left|E^{(1)}_{1}\right|\right|}_{F},{\left|\left|E^{(1)}_{2}\right|\right|}_{F}\right){\left|\left|P^{(1)}E^{(1)}_{1}-E^{(1)}_{2}P^{(1)}\right|\right|}_{F}$
(50)

and so

$\displaystyle D\left(\left.G_{1}^{(1)}\Box G_{1}^{(1)},G_{2}^{(1)}\Box
G_{2}^{(1)}\right|t_{c},\alpha_{c}\right)$ (51)
$\displaystyle\hbox{}\quad\quad\quad\leq\min\left({\left|\left|e^{\frac{t_{c}}{a_{c}}L(G^{(1)}_{1})}\right|\right|}_{F},{\left|\left|e^{t_{c}a_{c}L(G^{(1)}_{2})}\right|\right|}_{F}\right)D\left(G_{1}^{(1)},G_{2}^{(1)}\left.\right|t_{c},\alpha_{c}\right)$

An example of such a graph sequence is the sequence of two-dimensional square
grids, which are each the box product of two identical one-dimensional grids
i.e. path graphs: $\text{Sq}_{n}=\text{Pa}_{n}\Box\text{Pa}_{n}$.

###  3.4 Theory Summary

Triangle inequalities are proven for some members of the proposed family of
graph distortion or “distance” measures, including infinitesimal and finite
diffusion time, a power law for sparsity, and/or a power law for time
conversion factor between coarse and fine scales. However, the case of an
optimal (not power law) time conversion factor $\alpha$ needs to be
investigated by numerical experiment, and that requires new algorithms,
introduced in Section 4. We also show that in the case of distances between
graph box products, optimization over $P$ for the product graphs is bounded
above by a monotonic function of the optimum over the component graphs.

##  4 Numerical Methods for Optimal Time Conversion, $\alpha$

Optimizing the $\alpha$ parameter for conversion between coarse and fine time
scales in the proposed family of graph distance measures, in addition to
optimizing the prolongation matrix $P$ under transitive constraints
$\mathcal{C}(P)$, is a nontrivial numerical problem that in our experience
seems to require new methods. We develop such methods here and apply them to
investigate the resulting graph “distance” measure in the next section.

###  4.1 Algorithm Development

Algorithm 1 Algorithm for computing $\inf_{P,\alpha}\tilde{D}^{2}$

1:procedure LAP-Solve-Linear($\lambda^{(1)}$, $\lambda^{(2)}$, $\alpha$)

2: Compute cost matrix
$C_{ij}={\left(\frac{1}{\alpha}\lambda^{(1)}_{j}-\alpha\lambda^{(2)}_{i}\right)}^{2}$

3: Use an external solver for the LAP to find the matching matrix $M$ which
minimizes total cost $\sum_{i,j}M_{ij}C_{ij}$.

4: Return $M$.

5:end procedure

6:procedure Merge-Soln($M^{(1)}$, $M^{(2)}$,$\lambda^{(1)}$, $\lambda^{(2)}$,
$\alpha_{1}$, $\alpha_{2}$)

7: Construct $\tilde{e}^{(1)}$ and $\tilde{e}^{(2)}$ by removing
$\lambda^{(1)}_{j}$ and $\lambda^{(2)}_{i}$ whenever
$M^{(1)}_{ij}=M^{(2)}_{ij}=1$ (keeping track of the original indices of each
of the remaining values).

8: Compute $\alpha^{*}$ such that
$c\left(M^{(1)},\lambda^{(1)},\lambda^{(2)},\alpha^{*}\right)=c\left(M^{(2)},\lambda^{(1)},\lambda^{(2)},\alpha^{*}\right)$

9: if $\alpha^{*}\notin[\alpha_{1},\alpha_{2}]$, or $\tilde{e}^{(1)}$ is empty
then:

10: Mark the interval $[\alpha_{1},\alpha_{2}]$ as searched.

11: return $(M^{(1)},\alpha^{*})$ or $(M^{(2)},\alpha^{*})$.

12: else

13: Compute $M_{\text{sub}}$ as the solution given by $\textsc{LAP-Solve-
Linear}(\tilde{e}^{(2)},\tilde{e}^{(2)},\alpha^{*})$.

14: Construct a full-size solution $M$ for the original problem by starting
with $M=M^{(1)}\circ M^{(2)}$ (the elementwise product of $M^{(1)}$ and
$M^{(2)}$, which is 1 exactly where they agree and 0 otherwise) and inserting
a 1 for any pair $i,j$ which are matched by $M_{\text{sub}}$.

15: return $(M,\alpha^{*})$.

16: end if

17:end procedure

18:procedure $\tilde{D}^{2}$($G_{1},G_{2}$)

19: Compute graph Laplacians $L_{1}=L(G_{1})$ and $L_{2}=L(G_{2})$.

20: Compute $\lambda^{(1)}=\textsc{Eigenvals}(L_{1})$ and
$\lambda^{(2)}=\textsc{Eigenvals}(L_{2})$.

21: Initialize the list of known optimal matchings with $(n_{2}\times n_{1})$
matrices $\left[\begin{array}[]{c}I_{n_{1}}\\\ 0\end{array}\right]$ and
$\left[\begin{array}[]{c}0\\\ I_{n_{1}}\end{array}\right]$ and associated
$\alpha_{\text{low}}$ and $\alpha_{\text{high}}$.

22: while there are uncovered intervals within
$[\alpha_{\text{low}},\alpha_{\text{high}}]$: do

23: Order the current known pairs $(M_{i},\alpha_{i})$ by $\alpha$.

24: for each $(M_{i},\alpha_{i}),(M_{i+1},\alpha_{i+1})$ in this list: do

25: Compute $(M,\alpha^{*})=\textsc{Merge-
Soln}(M_{i},M_{i+1},\lambda^{(1)},\lambda^{(2)},\alpha_{1},\alpha_{2})$.

26: if $M=M_{i}$ or $M=M_{i+1}$ then

27: Mark the interval $[\alpha_{i},\alpha_{i+1}]$ as searched.

28: else

29: Insert $(M,\alpha^{*})$ into the appropriate place in the list of $M_{i}$.

30: end if

31: end for

32: end while

33: Find the unique global minimum of the cost of each found matching as a
function of $\alpha$.

34: Let $M^{*}$ be the matching which is minimal at its optimal $\alpha$ (out
of all matchings found).

35: Return the minimum cost $\tilde{D}^{2}(G_{1},G_{2}|M^{*},\alpha^{*})$.

36: (alternately, return the entire frontier
$\left\\{\left(M_{1},\alpha_{1}\right)\ldots\left(M_{m},\alpha_{m}\right)\right\\}$).

37:end procedure

Algorithm 2 Algorithm for computing $\sup_{t}\inf_{P,\alpha}D^{2}$

procedure LAP-Solve-Exponential($\lambda^{(1)}$, $\lambda^{(2)}$, $\alpha$, t)

Compute cost matrix
$C_{ij}={\left(e^{\frac{t}{\alpha}\lambda^{(1)}_{j}}-e^{\alpha
t\lambda^{(2)}_{i}}\right)}^{2}$

Use an external solver for the LAP to find the matching matrix $M$ which
minimizes total cost $\sum_{i,j}M_{ij}C_{ij}$.

Return $M$.

end procedure

procedure Merge-Soln($M^{(1)}$, $M^{(2)}$,$\lambda^{(1)}$, $\lambda^{(2)}$,
$\alpha_{1}$, $\alpha_{2}$, t)

Construct $\tilde{e}^{(1)}$ and $\tilde{e}^{(2)}$ by removing
$\lambda^{(1)}_{j}$ and $\lambda^{(2)}_{i}$ whenever
$M^{(1)}_{ij}=M^{(2)}_{ij}=1$ (keeping track of the original indices of each
of the remaining values).

Compute $\alpha^{*}$ such that
$c\left(M^{(1)},\lambda^{(1)},\lambda^{(2)},\alpha^{*},t\right)=c\left(M^{(2)},\lambda^{(1)},\lambda^{(2)},\alpha^{*},t\right)$

if $\alpha^{*}\notin[\alpha_{1},\alpha_{2}]$, or $\tilde{e}^{(1)}$ is empty
then:

Mark the interval $[\alpha_{1},\alpha_{2}]$ as searched.

return $(M^{(1)},\alpha^{*})$ or $(M^{(2)},\alpha^{*})$.

else

Compute $M_{\text{sub}}$ as the solution given by $\textsc{LAP-Solve-
Exponential}(\tilde{e}^{(2)},\tilde{e}^{(2)},\alpha^{*},t)$.

Construct a full-size solution $M$ for the original problem by starting with
$M=M^{(1)}\circ M^{(2)}$ (the elementwise product of $M^{(1)}$ and $M^{(2)}$,
which is 1 exactly where they agree and 0 otherwise) and inserting a 1 for any
pair $i,j$ which are matched by $M_{\text{sub}}$.

return $(M,\alpha^{*})$.

end if

end procedure

procedure T-Step($\lambda^{(1)}$, $\lambda^{(2)}$, Frontier
$=\left\\{\left(M_{1},\alpha_{1}\right),\left(M_{2},\alpha_{2}\right)\ldots\left(M_{m},\alpha_{m}\right)\right\\}$,
t)

while Frontier is still expanding do

for each $(M_{i},\alpha_{i}),(M_{i+1},\alpha_{i+1})$ in Frontier do

if $M_{i}$ and $M_{i+1}$ have already been checked then

pass

else

Compute $(M,\alpha^{*})=\textsc{Merge-
Soln}(M_{i},M_{i+1},\lambda^{(1)},\lambda^{(2)},\alpha_{i},\alpha_{i+1},t)$.

if $M=M_{i}$ or $M=M_{i+1}$ then

pass

else

Insert $(M,\alpha^{*})$ into the appropriate place in Frontier.

end if

end if

end for

end while

Return Frontier.

end procedure

procedure $D^{2}$($G_{1},G_{2}$)

Compute graph Laplacians $L_{1}=L(G_{1})$ and $L_{2}=L(G_{2})$.

Compute $\lambda^{(1)}=\textsc{Eigenvals}(L_{1})$ and
$\lambda^{(2)}=\textsc{Eigenvals}(L_{2})$.

Initialize Frontier as the frontier of
$\left\\{\left(M_{i},\alpha_{i}\right)\right\\}$ found by Algorithm 1.

Set $t_{c}=t_{\text{init}}$ (some low value e.g. $10^{-3}$).

Compute $D^{2}_{\text{max}}=\min_{i}D^{2}(G_{1},G_{2}|t_{c},\alpha_{i},M_{i})$

while $D^{2}_{\text{max}}$ is increasing do

$t_{c}+=\Delta t$

$\textsc{Frontier}=\textsc{T-Step}(\lambda^{(1)},\lambda^{(2)},\textsc{Frontier},t_{c})$

$D^{2}_{\text{max}}=\max_{i}D^{2}(G_{1},G_{2}|t_{c},\alpha_{i},M_{i})$

end while

Return $\sup_{t}\min_{i}D^{2}(G_{1},G_{2}|t_{c},\alpha_{i},M_{i})$

end procedure

In this section, we describe the algorithm used to calculate upper bounds on
graph distances as the joint optima (over $P$, $t$, and $\alpha$) of the
distance equations Equation 8 and Equation 9, under orthogonality constraints
only, i.e. the case $\mathcal{C}(P)=\\{P|P^{T}P=I\\}$. At the core of both
algorithms is a subroutine to solve the Linear Assignment Problem (LAP - see
Equation (5)) repeatedly, in order to find the subpermutation matrix which is
optimal at a particular value of $\alpha$. Namely, we are interested in
calculating $\tilde{D}$ as

$\displaystyle\tilde{D}(G_{1},G_{2})$ $\displaystyle=\min_{\alpha}f(\alpha)$
(52) where $\displaystyle f(\alpha)$
$\displaystyle=\inf_{P|P^{T}P=I}\left|\left|\frac{1}{\alpha}PL(G_{1})-\alpha
L(G_{2})P\right|\right|$ which, for orthogonality or any other compact
constraint
$\displaystyle=\min_{P|P^{T}P=I}\left|\left|\frac{1}{\alpha}PL(G_{1})-\alpha
L(G_{2})P\right|\right|.$

However, we have found that the unique structure of this optimization problem
admits a specialized procedure which is faster and more accurate than nested
univariate optimization of $\alpha$ and $t$ (where each innermost function
evaluation consists of a full optimization over $P$ at some $t,\alpha$). We
first briefly describe the algorithm used to find the optimal $P$ and $\alpha$
for $\tilde{D}^{2}$. The formal description of the algorithm is given by
Algorithm 1. In both cases, we reduce the computational complexity of the
optimization over $P$ by imposing the additional constraint that $P$ must be a
subpermutation matrix when rotated into the spectral basis (we define
subpermutations in the proof of Theorem 4.1). This constraint is compatible
with the orthogonality constraint (all subpermutation matrices are orthogonal,
but not vice versa). The tradeoff of this reduction of computational
complexity is that we can only guarantee that our optima are upper bounds of
the optima over all orthogonal $P$. However, in practice, this bound seems to
be tight: we have yet to find an example where orthogonally-constrained
optimization was able to improve in objective function value over optimization
constrained to subpermutation matrices. Therefore, we shall for the remainder
of this paper refer to the optima calculated as distance values, when strictly
they are distance upper bounds. We also note here that a distance lower bound
is also possible to calculate by _relaxing_ the constraints in
$\mathcal{C}(P)$ (for instance, by replacing the optimization over all $P$
with a less constrained matching problem - see Appendix 10.1).

(a)
(b)
(c)

Figure 2: Three plots demonstrating characteristics of distance calculation
between a $(7\times 7)$ grid and an $(8\times 8)$ grid.
(a): a plot illustrating unimodality of
$D^{2}\left(\left.\text{Sq}_{7},\text{Sq}_{8}\right|t\right)=\inf_{\alpha>0}\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{\frac{t}{a}L(\text{Sq}_{7})}-e^{t\alpha
L(\text{Sq}_{8})}P\right|\right|}_{F}^{2}$ as a function of $t$. The peak, at
$t\approx.318$, yields the distance
$D^{2}\left(\text{Sq}_{7},\text{Sq}_{8}\right)$.
(b): Plot illustrating the discontinuity and multimodality of the linear
version of distance. Each gray curve represents a function
$f_{P_{c}}(\alpha_{c})=\tilde{D}^{2}\left(\text{Sq}_{7},\text{Sq}_{8}|\alpha_{c},P_{c}\right)$.
The thicker curve is the lower convex hull of the thinner curves as a function
of $\alpha$, that is:
$\boldsymbol{f}(\alpha_{c})=\inf_{P|\mathcal{C}(P)}\tilde{D}^{2}\left(\text{Sq}_{7},\text{Sq}_{8}|\alpha_{c}\right)$.
We see that $\boldsymbol{f}(\alpha)$ is continuous, but has discontinuous
slope, as well as several local optima (marked by arrowheads). These
properties make $\tilde{D}$ difficult to optimize, necessitating the
development of Algorithm 1.
(c): As in (b), but with
$D^{2}\left(\text{Sq}_{7},\text{Sq}_{8}|t=.318\right)$ plotted instead of
$\tilde{D}^{2}$. This $t$ value is the location of the maximum in the leftmost
plot.

####  4.1.1 Optimization of $\tilde{D}^{2}$

Joint optimization of $\tilde{D}^{2}$ over $\alpha$ and $P$ is a nested
optimization problem (see onoyama2000method  and sinha2018review  for a
description of nested optimization), with potential combinatorial optimization
over $P$ dependent on each choice of $\alpha$. Furthermore, the function
$f(\alpha)=\inf_{P|\mathcal{C}(P)}\tilde{D}^{2}(G_{1},G_{2}|\alpha)$ is both
multimodal and continuous but with in general discontinuous derivative (See
Figure 2). Univariate optimization procedures such as Golden Section Search
result in many loops of some procedure to optimize over $P$, which in our
restricted case must each time compute a full solution to a LAP with
$n_{2}\times n_{1}$ weights. In our experience, this means that these
univariate methods have a tendency to get stuck in local optima. We reduce the
total number of calls to the LAP solver, as well as the size of the LAPs
solved, by taking advantage of several unique properties of the optimization
as a function of $\alpha$. When the optimal $P^{(1)}$ and $P^{(2)}$ are known
for $\alpha_{1}$ and $\alpha_{2}$, then for any $\alpha_{c}$ such that
$\min(\alpha_{1},\alpha_{2})\leq\alpha_{c}\leq\max(\alpha_{1},\alpha_{2})$,
the optimal $P^{(c)}$ at $\alpha_{c}$ must satisfy: $P^{(1)}_{ij}=1\wedge
P^{(2)}_{ij}=1\implies P^{(c)}_{ij}=1$ (see Theorem 4.2). Thus, the
optimization over $P$ at $\alpha_{c}$ is already partially solved given the
solutions at $\alpha_{1}$ and $\alpha_{2}$, and so we need only re-compute the
remaining (smaller) subproblem on the set of assignments where $P^{(1)}$ and
$P^{(2)}$ disagree. This has two consequences for our search over $\alpha$:
First, the size of LAP problems which must be solved at each step decreases
over time (as we find $P$-optima for a denser and denser set of $\alpha$).
Secondly, these theoretical guarantees mean that we can mark intervals of
$\alpha$-values as being explored (meaning we have provably found the $P$
which are optimal over the interval) and thus do not have to perform the
relatively expensive optimization over $P$ for any $\alpha$ in that interval.

####  4.1.2 Optimization of $D^{2}$

Many of the theoretical guarantees underlying our algorithm for computing
$\tilde{D}^{2}$ no longer hold for the exponential version of the distance. We
adapt our linear-version procedure into an algorithm for computing this
version, with the caveat that the lack of these guarantees means that our
upper bound on the exponential version may be looser than that on the linear
version. It is still clearly an upper bound, since the $\alpha$ and $P$ found
by this procedure satisfy the given constraints $\alpha>0$ and $P^{T}P=I$. In
particular, we have observed cases where the exponential-distance analog of
Theorem 4.2 would not hold, meaning we cannot rule out $\alpha$-intervals as
we can in the linear version. Thus, this upper bound may be looser than the
that computed for the linear objective function.

For the exponential version of the algorithm, we first compute the list of
optimal $P$ for the linear version, assuming (since $e^{tL}\approx I+L$ for
very small $t$) that this is also the list of optimal $P$ for the exponential
version of the objective function at some low $t$. We proceed to increment $t$
with some step size $\Delta t$, in the manner of a continuation method
allgower2012numerical . At each new $t$ value, we search for new optimal $P$
along the currently known frontier of optima as a function of $\alpha$. When a
new $P$ is found as the intersection of two known $P_{i},P_{i+1}$, it is
inserted into the list, which is kept in order of increasing $\alpha$. For
each $P$ in this frontier, we find the optimal $\alpha$, keeping $P$ and $t$
constant. Assuming $\inf_{P}\inf_{\alpha}D^{2}(G_{1},G_{2}|t_{c})$ is unimodal
as a function of $t_{c}$, we increase $t_{c}$ until
$\inf_{P}\inf_{\alpha}D^{2}(G_{1},G_{2}|t_{c})\geq\inf_{P}\inf_{\alpha}D^{2}(G_{1},G_{2}|t_{c}+\Delta
t)$, storing all $P$ matrices found as optima at each $t_{c}$ value. $P$ which
were on the lower convex hull at some prior value of $t$ but not the current
value are retained, as they may regain optimality for some $\alpha$-range at a
future value of $t$ (we have observed this, in practice). For this list
$P_{1},P_{2}\ldots P_{m}$, we then compute
$\sup_{t}\inf_{\alpha}\inf_{i}D^{2}(G_{1},G_{2}|P_{i})$. Since the exponential
map is continuous, and we are incrementing $t$ by very small steps, we also
propose the further computational shortcut of storing the list of optimal
$\alpha$ at time $t$ to use as starting points for the optimization at
$t+\Delta t$. In practice, this made little difference in the runtime of our
optimization procedure.

###  4.2 Algorithm Correctness Proof

######  Theorem 4.1.

For any two graphs $G_{1}$ and $G_{2}$ with Laplacians $L(G_{1})$ and
$L(G_{2})$, for fixed $\alpha$, the optimization over $P$ given in the
innermost loop of Equation 9 is upper bounded by a Linear Assignment Problem
as defined in Equation (5). This LAP is given by taking $R$ to be the
eigenvalues $\lambda^{(1)}_{j}$ of $L(G_{1})$ and $S$ to be the eigenvalues
$\lambda^{(2)}_{i}$ of $L(G_{2})$, with the cost of a pair (equivalently, one
entry of the cost matrix $C$) given by

$\displaystyle
C_{ij}=c(s_{i},r_{j})=c\left(\lambda^{(2)}_{i},\lambda^{(1)}_{j}\right)={\left(\frac{1}{\alpha}{\lambda^{(1)}_{j}}-\alpha{\lambda^{(2)}_{i}}\right)}^{2}$
(53)

###### Proof.

$L(G_{1})$ and $L(G_{2})$ are both real symmetric matrices, so they may be
diagonalized as $L(G_{i})=U_{i}\Lambda_{i}U_{i}^{T}$, where the $U_{i}$ are
rotation matrices, and the $\Lambda_{i}$ are diagonal matrices with the
eigenvalues $\lambda^{(i)}_{1},\lambda^{(i)}_{2}\ldots\lambda^{(i)}_{n_{i}}$
along the diagonal. Because the Frobenius norm is invariant under rotation, we
have:

$\displaystyle\tilde{D}^{2}(G_{1},G_{2})$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}{\left|\left|\frac{1}{\alpha}PL(G_{1})-\alpha
L(G_{2})P\right|\right|}^{2}_{F}$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}{\left|\left|\frac{1}{\alpha}U_{2}^{T}PL(G_{1})U_{1}-\alpha
U_{2}^{T}L(G_{2})PU_{1}\right|\right|}^{2}_{F}$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}{\left|\left|\frac{1}{\alpha}U_{2}^{T}PU_{1}\Lambda_{1}U_{1}^{T}U_{1}-\alpha
U_{2}^{T}U_{2}\Lambda_{2}U_{2}^{T}PU_{1}\right|\right|}^{2}_{F}$ (54)
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}{\left|\left|\frac{1}{\alpha}U_{2}^{T}PU_{1}\Lambda_{1}-\alpha\Lambda_{2}U_{2}^{T}PU_{1}\right|\right|}^{2}_{F}.$
Because the $U_{i}$ are orthogonal, the transformation
$\tilde{P}=U_{2}^{T}PU_{1}$ preserves orthogonality, so
$\displaystyle\tilde{D}^{2}(G_{1},G_{2})$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}{\left|\left|\frac{1}{\alpha}P\Lambda_{1}-\alpha\Lambda_{2}P\right|\right|}^{2}_{F}$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}{\left|\left|\frac{1}{\alpha}\Lambda_{1}\right|\right|}^{2}_{F}+{\left|\left|\frac{}{}\alpha\Lambda_{2}P\right|\right|}^{2}_{F}-2\text{Tr}\left[P^{T}\Lambda_{2}P\Lambda_{1}\right]$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}\left(\text{Tr}\left[\frac{1}{\alpha^{2}}\Lambda_{1}^{2}\right]+\text{Tr}\left[\alpha^{2}P^{T}\Lambda_{2}^{2}P\right]-2\text{Tr}\left[P^{T}\Lambda_{2}P\Lambda_{1}\right]\right)$
writing $P=[p_{ij}]$, $\displaystyle\tilde{D}^{2}(G_{1},G_{2})$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}\left(\frac{1}{\alpha^{2}}\sum\limits_{j=1}^{n_{1}}{\lambda^{(1)}_{j}}^{2}+\alpha^{2}\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}p_{ij}^{2}{\lambda^{(2)}_{i}}^{2}\right.$
(55)
$\displaystyle\left.\qquad\qquad\qquad\qquad\qquad-2\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}p_{ij}^{2}{\lambda^{(2)}_{i}}{\lambda^{(1)}_{j}}\right)$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}\left(\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}p_{ij}^{2}\left(\frac{1}{\alpha^{2}}{\lambda^{(1)}_{j}}^{2}-2{\lambda^{(2)}_{i}}{\lambda^{(1)}_{j}}+\alpha^{2}{\lambda^{(2)}_{i}}^{2}\right)\right)$
$\displaystyle=\inf_{\alpha>0}\inf_{P^{T}P=I}\left(\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}p_{ij}^{2}{\left(\frac{1}{\alpha}{\lambda^{(1)}_{j}}-\alpha{\lambda^{(2)}_{i}}\right)}^{2}\right)$
(56)

For any given $\alpha$,

$\displaystyle\inf_{P^{T}P=I}\left(\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}p_{ij}^{2}{\left(\frac{1}{\alpha}{\lambda^{(1)}_{j}}-\alpha{\lambda^{(2)}_{i}}\right)}^{2}\right)$
(57)
$\displaystyle\qquad\leq\inf_{\tilde{P}|\text{subperm}(\tilde{P})}\left(\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}\tilde{p}_{ij}^{2}{\left(\frac{1}{\alpha}{\lambda^{(1)}_{j}}-\alpha{\lambda^{(2)}_{i}}\right)}^{2}\right)\quad,$

where $\text{subperm}(\tilde{P})$ could be any other condition more strict
than the constraint $P^{T}P=I$. Here we take this stricter constraint to be
the condition that $\tilde{P}$ is a subpermutation matrix: an orthogonal
matrix (i.e. $\tilde{P}^{T}\tilde{P}=I)$ for which
$\tilde{P}\in\\{0,1\\}^{n_{2}\times n_{1}}$. Equivalently, a subpermutation
matrix is a $\\{0,1\\}$-valued matrix $[\tilde{p}_{ij}]$ such that for each
$i\in\\{1,\ldots n_{1}\leq n_{2}\\}$, exactly one $j\in\\{1,\ldots n_{2}\geq
n_{1}\\}$ takes the value 1 rather than 0 (so
$\sum_{j=1}^{n_{2}}\tilde{P}_{ji}=1$), and for each $j\in\\{1,\ldots n_{2}\geq
n_{1}\\}$, either zero or one $i\in\\{1,\ldots n_{1}\leq n_{2}\\}$ takes the
value 1 rather than 0 (so $\sum_{i=1}^{n_{1}}\tilde{P}_{ji}\leq 1$).

Furthermore, this optimization is exactly a linear assignment problem of
eigenvalues of $L(G_{1})$ to $L(G_{2})$, with the cost of a pair
$\left(\lambda^{(1)}_{j},\lambda^{(2)}_{i}\right)$ given by

$c\left(\lambda^{(1)}_{j},\lambda^{(2)}_{i}\right)={\left(\frac{1}{\alpha}{\lambda^{(1)}_{j}}-\alpha{\lambda^{(2)}_{i}}\right)}^{2}$

Note also that the same argument applies to the innermost two optimizations of
the calculation of $D^{2}$ (the exponential version of the diffusion distance)
as well as $D^{2}_{r}$. In the $D^{2}$ case the entries of the cost matrix are
instead given by

$c\left(\lambda^{(1)}_{j},\lambda^{(2)}_{i}\right)={\left(e^{\frac{1}{\alpha}{\lambda^{(1)}_{j}}}-e^{\alpha{\lambda^{(2)}_{i}}}\right)}^{2}$

If we instead loosen the constraints on $P$, we can calculate a lower bound on
the distance. See Appendix 10.1 for lower bound details.

∎

Recall that our definition of a ‘matching’ in Section 1.3 was a $P$ matrix
representing a particular solution to the linear assignment problem with costs
given as in Equation (53). For given $G_{1},G_{2}$, and some matching $M$, let

$\displaystyle
f_{M}(\alpha)=\tilde{D}^{2}(G_{1},G_{2}|\alpha,U_{2}^{T}MU_{1})$ (58)

where $U_{1},U_{2}$ diagonalize $L_{1}$ and $L_{2}$ as in Equation (54).

######  Lemma 2.

For two unique matchings $M_{1}$ and $M_{2}$ (for the same $G_{1},G_{2}$) the
equation $f_{M_{1}}(\alpha)-f_{M_{2}}(\alpha)=0$ has at most one real positive
solution in $\alpha$. This follows from the fact that when $P$ and $t$ are
fixed, the objective function is a rational function in $\alpha$ (see Equation
(55)), with a quadratic numerator and an asymptote at $\alpha=0$.

###### Proof.

By Equation (55), we have

$\displaystyle f_{M_{1}}(\alpha)-f_{M_{2}}(\alpha)=$
$\displaystyle\left(\frac{1}{\alpha^{2}}\sum\limits_{j=1}^{n_{1}}{\lambda^{(1)}_{j}}^{2}+\alpha^{2}\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{1}]}_{ij}^{2}{\lambda^{(2)}_{i}}^{2}-2\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{1}]}_{ij}^{2}{\lambda^{(2)}_{i}}{\lambda^{(1)}_{j}}\right)$
(59)
$\displaystyle\quad-\left(\frac{1}{\alpha^{2}}\sum\limits_{j=1}^{n_{1}}{\lambda^{(1)}_{j}}^{2}+\alpha^{2}\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{2}]}_{ij}^{2}{\lambda^{(2)}_{i}}^{2}-2\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{2}]}_{ij}^{2}{\lambda^{(2)}_{i}}{\lambda^{(1)}_{j}}\right)$
(60)
$\displaystyle=\alpha^{2}\left(\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{1}]}_{ij}^{2}{\lambda^{(2)}_{i}}^{2}-\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{2}]}_{ij}^{2}{\lambda^{(2)}_{i}}^{2}\right)$
(61)
$\displaystyle\quad+\left(2\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{2}]}_{ij}^{2}{\lambda^{(2)}_{i}}{\lambda^{(1)}_{j}}-2\sum\limits_{i=1}^{n_{2}}\sum\limits_{j=1}^{n_{1}}{[M_{1}]}_{ij}^{2}{\lambda^{(2)}_{i}}{\lambda^{(1)}_{j}}\right)$
(62) Abbreviating the sums, we have
$\displaystyle\alpha^{2}\left(A_{1}-A_{2}\right)+\left(C_{2}-C_{1}\right)$
$\displaystyle=0$ (64) and so
$\displaystyle\alpha=\pm\sqrt{\frac{C_{2}-C_{1}}{A_{1}-A_{2}}}$ (65)

Since $A_{1},A_{2},C_{1},C_{2}$ are all nonnegative reals, at most one of
these roots is positive. ∎

We will say that a matching $M$ “assigns” $j$ to $i$ if and only if
$M_{ij}=1$.

######  Theorem 4.2.

If two matchings $M_{1}$ and $M_{3}$ which yield optimal upper bounds for the
linear distance $\tilde{D}^{2}$ (at $\alpha_{1}\leq\alpha$ and
$\alpha_{3}\geq\alpha$ respectively) agree on a set of assignments, then the
optimal $M$ at $\alpha$ must also agree with that set of assignments.

###### Proof.

We need the following lemmas:

######  Lemma 3.

If an optimal matching assigns $i$ to $m(i)$ (so that eigenvalue
$\lambda^{(1)}_{i}$ of $G_{1}$ is paired with $\lambda^{(2)}_{f(i)}$ of
$G_{2}$ in the sum of costs Equation (53)), then the sequence
$m(1),m(2),\ldots m(n_{1})$ is monotonic increasing.

###### Proof.

This follows from the fact that the two sequences of eigenvalues are monotonic
nondecreasing, so if there’s a ‘crossing’ ($i_{1}<i_{2}$ but
$m(i_{2})<m(i_{1})$) then the new matching obtained by uncrossing those two
pairs (performing a 2-opt step as defined in croes1958method ) has strictly
lesser objective function value. Hence an optimal matching can’t contain any
such crossings. ∎

######  Lemma 4.

For all positive real $\alpha^{*}\geq\epsilon>0$, let $M_{1}$ be an optimal
matching at $\alpha^{*}-\epsilon$ and $M_{2}$ be optimal at
$\alpha^{*}+\epsilon$. For $1\leq i\leq n_{1}$, let $s_{1}(i)$ and $s_{2}(i)$
be the indices of $\lambda^{(2)}$ paired with $i$ in $M_{1}$ and $M_{2}$,
respectively. Then for all $i$, $s_{1}(i)\leq s_{2}(i)$.

###### Proof.

Define a “run” for $s_{1},s_{2}$ as a sequence of consecutive indices
$l,l+1,\ldots l+k$ in $[1,n_{1}]$ such that for any $l$, $l+1$:
$\min(s_{1}(l+1),s_{2}(l+1))<\max(s_{1}(l),s_{2}(l))$. The following must be
true about a “run”:

  1. 1.

Within a run, either $s_{1}(l)<s_{2}(l)$ or $s_{1}(l)>s_{2}(l)$ for all $l$.
Otherwise, we have one or more crossings (as in Lemma 3): for some $l$ we have
$s_{1}(l)>s_{1}(l+1)$ or $s_{2}(l)>s_{2}(l+1)$. Any crossing may be uncrossed
for a strictly lower objective function value - violating optimality of
$M_{1}$ or $M_{2}$.

  2. 2.

Any pair of matchings as defined above consists of a sequence of runs, where
we allow a run to be trivial i.e. be a single index.

Next, we show that within a run, we must have $s_{1}(i)<s_{2}(i)$ for all $i$.
Let $S=\\{l,l+1,\ldots l+k\\}$ be a run. By optimality of $M_{1}$, $M_{2}$ at
$\alpha^{*}-\epsilon$ and $\alpha^{*}+\epsilon$ respectively, we have:

$\displaystyle\sum_{i\in
S}{\left(\frac{1}{\alpha^{*}-\epsilon}\lambda^{(1)}_{i}-(\alpha^{*}-\epsilon)\lambda^{(2)}_{s_{1}(i)}\right)}^{2}<\sum_{i\in
S}{\left(\frac{1}{\alpha^{*}-\epsilon}\lambda^{(1)}_{i}-(\alpha^{*}-\epsilon)\lambda^{(2)}_{s_{2}(i)}\right)}^{2}$
and $\displaystyle\sum_{i\in
S}{\left(\frac{1}{\alpha^{*}+\epsilon}\lambda^{(1)}_{i}-(\alpha^{*}+\epsilon)\lambda^{(2)}_{s_{2}(i)}\right)}^{2}<\sum_{i\in
S}{\left(\frac{1}{\alpha+\epsilon}\lambda^{(1)}_{i}-(\alpha+\epsilon)\lambda^{(2)}_{s_{1}(i)}\right)}^{2}.$

Respectively, these simplify to

$\displaystyle-\sum_{i\in
S}\left(\lambda^{(2)}_{s_{1}(i)}-\lambda^{(2)}_{s_{2}(i)}\right)\left(-2\lambda^{(i)}_{i}+{(\alpha^{*}-\epsilon)}^{2}\left(\lambda^{(2)}_{s_{1}(i)}+\lambda^{(2)}_{s_{2}(i)}\right)\right)>0$
and $\displaystyle\sum_{i\in
S}\left(\lambda^{(2)}_{s_{1}(i)}-\lambda^{(2)}_{s_{2}(i)}\right)\left(-2\lambda^{(i)}_{i}+{(\alpha^{*}+\epsilon)}^{2}\left(\lambda^{(2)}_{s_{1}(i)}+\lambda^{(2)}_{s_{2}(i)}\right)\right)>0.$

Summing these inequalities and cancelling $-2\lambda^{(i)}_{i}$, we have:

$\displaystyle\sum_{i\in
S}\left\\{{(\alpha^{*}+\epsilon)}^{2}\left({\left(\lambda^{(2)}_{s_{1}(i)}\right)}^{2}+{\left(\lambda^{(2)}_{s_{2}(i)}\right)}^{2}\right)\right.$
$\displaystyle\quad\quad\quad-\left.{(\alpha^{*}-\epsilon)}^{2}\left({\left(\lambda^{(2)}_{s_{1}(i)}\right)}^{2}+{\left(\lambda^{(2)}_{s_{2}(i)}\right)}^{2}\right)\right\\}>0.$

Summing and reducing gives us

$4\alpha^{*}\epsilon\left(\sum_{i\in
S}{\left(\lambda^{(2)}_{s_{1}(i)}\right)}^{2}-\sum_{i\in
S}{\left(\lambda^{(2)}_{s_{2}(i)}\right)}^{2}\right)>0$

and so

$\sum_{i\in S}{\left(\lambda^{(2)}_{s_{1}(i)}\right)}^{2}>\sum_{i\in
S}{\left(\lambda^{(2)}_{s_{2}(i)}\right)}^{2}.$

However, since the $\lambda^{(2)}_{j}$ are monotonic nondecreasing, this means
we cannot also have $s_{1}(i)>s_{2}(i)$ for all $i\in S$, since that would
imply

$\sum_{i=1}^{n_{1}}{\left(\lambda^{(2)}_{s_{1}(i)}\right)}^{2}<\sum_{i=1}^{n_{1}}{\left(\lambda^{(2)}_{s_{2}(i)}\right)}^{2}.$

Therefore, in a run of arbitrary length, all indices must move ‘forward’
(meaning that $s_{1}(i)<s_{2}(i)$ for all $i$ in the run), and so (since any
pair of matchings optimal at such $\alpha$ define a set of runs) we must have
$s_{1}(i)\leq s_{2}(i)$. This completes the proof of the lemma.

∎

Thus, for three matchings $M_{1},M_{2},M_{3}$ which optimal at a sequence of
$\alpha_{1}\leq\alpha_{2}\leq\alpha_{3}$, we must have $s_{1}(i)\leq
s_{2}(i)\leq s_{3}(i)$ for all $i$. In particular, if $s_{1}(i)=s_{3}(i)$, we
must also have $s_{1}(i)=s_{2}(i)=s_{3}(i)$. ∎

######  Theorem 4.3.

If two matchings $M_{1}$ and $M_{3}$ yield optimal upper bounds for the linear
distance $\tilde{D}^{2}$ at $\alpha_{1}$ and $\alpha_{3}$ respectively, and
$f_{M_{1}}(\alpha_{2})=f_{M_{2}}(\alpha_{2})$ for some $\alpha_{2}$ s.t.
$\alpha_{1}\leq\alpha_{2}\leq\alpha_{3}$, then either (1) $M_{1}$ and $M_{3}$
are optimal over the entire interval $[\alpha_{1},\alpha_{3}]$ or (1) some
other matching $M_{2}$ improves over $M_{1}$ and $M_{3}$ at $\alpha_{2}$.

###### Proof.

This follows directly from the facts that $f_{M_{1}}(\alpha)$ and
$f_{M_{2}}(\alpha)$ (as defined in Equation (58)), can only meet at one real
positive value of $\alpha$ (Lemma 2). Say that the cost curves for $M_{1}$
(known to be optimal at $\alpha=\alpha_{1}$) and $M_{3}$ (optimal at
$\alpha=\alpha_{3}$) meet at $\alpha=\alpha_{2}$, and furthermore assume that
$\alpha_{1}\leq\alpha_{2}\leq\alpha_{3}$. If some other matching $M_{2}$
improves over (meaning, has lesser obj. function value as a function of
$\alpha$) $M_{1}$ or $M_{3}$ anywhere in the interval
$[\alpha_{1},\alpha_{3}]$, it must improve over both at $\alpha=\alpha_{2}$,
since it may intersect each of these cost curves at most once on this
interval. If $M_{1}$ and $M_{3}$ are both optimal at their intersection point
(meaning no such distinct $M_{2}$ exists) then we know that no other matching
improves on either of them over the the interval $[\alpha_{1},\alpha_{3}]$ and
may therefore mark it as explored during the outermost loop (otimization over
$\alpha$) of Algorithm 1. ∎

Together, the preceeding properties verify that our algorithm will indeed find
the joint optimum over all $\alpha$ and $P$ (for fixed $t=c$, for $\tilde{D}$,
subject to subpermutation constraints on $P$): it allows us to find the entire
set of $P$ subpermutation matrices which appear on the lower convex hull of
distance as a function of alpha.

###  4.3 Implementation Details

We implement Algorithms 1 and 2 in the programming language “Python” (version
3.6.1) rossum1995python . Numerical arrays were stored using the _numpy_
package van2011numpy . Our inner LAP solver was the package _lapsolver_
heindl2018lapsolver . Univariate optimization over $t$ and $\alpha$ was
performed with the ‘bounded’ method of the _scipy.optimize_ package scipy ,
with bounds set at $[0,10.0]$ for each variable and a input tolerance of
$10^{-12}$. Laplacians were computed with the laplacian method from the
package _networkX_ hagberg2008exploring , and their eigenvalues were computed
with scipy.linalg.eigh.

Because of numerical precision issues arising during eigenvalue computation,
it can be difficult to determine when two matchings agree, using eigenvalue
comparison. In practice we ignore this issue and assume that two matchings are
only identical if they associate the same indices of the two lists of
eigenvalues. This means we may be accumulating multiple equivalent
representations of the same matching (up to multiplicity of eigenvalues)
during our sweeps through $t$ and $\alpha$. We leave mitigating this
inefficiency for future work.

##  5 Numerical Experiments

###  5.1 Graph Lineages

In this subsection we introduce several graph lineages for which we will
compute various intra- and inter-lineage distances. Three of these are well-
known lineages of graphs, and the fourth is defined in terms of a product of
complete graphs:

  * •

Path Graphs ($\text{Pa}_{n}$): 1D grid graphs of length $n$, with aperiodic
boundary conditions.

  * •

Cycle Graphs ($\text{Cy}_{n}$): 1D grid graphs of length $n$, with periodic
boundary conditions.

  * •

Square Grid Graphs ($\text{Sq}_{n}$): 2D grid graphs of dimensions $n$, with
aperiodic boundary conditions. $\text{Sq}_{n}=\text{Pa}_{n}\Box\text{Pa}_{n}$

  * •

“Multi-Barbell” Graphs ($\text{Ba}_{n}$): Constructed as $\text{Cy}_{n}\Box
K_{n}$, where $K_{n}$ is the complete graph on $n$ vertices.

Examples of each of these graph lineages can be seen in Figure 3.

2D Grids

Paths

Cycles

k-Barbell graphs

Figure 3: Graph lineages used in multiple numerical experiments in this
section.

###  5.2 Numerical Optimization Methods

We briefly discuss here the other numerical methods we have used to calculate
$\tilde{D}^{2}$ and $D^{2}$.

  * •

Nelder-Mead in Mathematica: For very small graph pairs ($n_{1}\times n_{2}\leq
100$) we are able to find optimal $P,\alpha,t$ using constrained optimization
in Mathematica 11.3 wolfram2018mathematica  using NMinimize, which uses
Nelder-Mead as its backend by default. The size limitation made this approach
unusable for any real experiments.

  * •

We also tried a variety of codes specialized for numeric optimization subject
to orthogonality constraints. These included (1) the python package PyManopt
townsend2016pymanopt , a code designed for manifold-constrained optimization;
(2) gradient descent in Tensorflow using the penalty function
$g(P)=c{\left|\left|P^{T}P-I\right|\right|}_{F}$ (with $c\ll 1$ a small
positive constant weight) to maintain orthogonality, as well as (3) an
implementation of the Cayley reparametrization method from wen2013feasible
(written by the authors of that same paper). In our experience, these codes
were slower, with poorer scaling with problem size, than combinatorial
optimization over subpermutation matrices, and did not produce improved
results on our optimization problem.

  * •

We compare in more detail two methods of joint optimization over $\alpha$ and
$P$ when $P$ is constrained to be a subpermutation matrix in the diagonal
basis for $L(G_{1})$ and $L(G_{2})$. Specifically, we compare our approach
given in Algorithm 1 to univariate optimization over $\alpha$, where each
function evaluation consists of full optimization over $P$. Figure 4 shows the
results of this experiment. We randomly sample pairs of graphs as follows:

    1. 1.

$n_{1}$ is drawn uniformly from $[5,120]$.

    2. 2.

$n_{2}$ is drawn uniformly from $[n_{1},n_{1}+60]$.

    3. 3.

$G_{1}$ and $G_{2}$ are generated by adding edges according to a Bernoulli
distribution with probability $p$. We ran 60 trials for each $p$ in { .125,
.25, .375, .5, .625, .75, .875 }.

We compute the linear version of distance for each pair. Because our algorithm
finds all of the local minima as a function of alpha, we compute the cost of
the golden section approach as the summed cost of multiple golden section
searches in alpha: one GS search starting from the initial bracket
$[0.618\alpha^{*},1.618\alpha^{*}]$ for each local minimum $\alpha^{*}$ found
by our algorithm. We see that our algorithm is always faster by at least a
factor of 10, and occasionally faster by as much as a factor of $10^{3}$. This
can be attributed to the fact that the golden section search is unaware of the
structure of the linear assignment problem: it must solve a full $n_{2}\times
n_{2}$ linear assignment problem for each value of $\alpha$ it explores. In
contrast, our algorithm is able to use information from prior calls to the LAP
solver, and therefore solves a series of LAP problems whose sizes are
monotonically nonincreasing.

Figure 4: Comparison of runtimes for our algorithm and bounded golden section
search over the same interval $[10^{-6},10]$. Runtimes were measured by a
weighted count of evaluations of the Linear Assignment Problem solver, with an
$n\times n$ linear assignment problem counted as $n^{3}$ units of cost.
Because our algorithm recovers the entire lower convex hull of the objective
function as a function of $\alpha$, we compute the cost of the golden section
search as the summed cost of multiple searches, starting from an interval
bracketing each local optimum found by our algorithm. We see that our
algorithm is much less computationally expensive, sometimes by a factor of
$10^{3}$. The most dramatic speedup occurs in the regime where $n_{1}<<n_{2}$.
Graphs were generated by drawing $n_{1}$ uniformly from $[5,120]$, drawing
$n_{2}$ uniformly from $[n_{1},n_{1}+60]$, and then adding edges according to
a Bernoulli distribution with $p$ in { .125, .25, .375, .5, .625, .75, .875 }
(60 trials each).

###  5.3 Experiments

####  5.3.1 Triangle Inequality violation of $D$ (Exponential Distance) and
$\tilde{D}$ (Linear Distance).

As stated in Section 2.5, our full graph dissimilarity measure does not
necessarily obey the triangle inequality. In this section we systematically
explore conditions under which the triangle inequality is satisfied or not
satisfied. We generate triplets $G_{1},G_{2},G_{3}$ of random graphs of sizes
$n_{i}$ for $n_{1}\in[5,30]$, $n_{2}\in[n_{1},n_{1}+30]$, and
$n_{3}\in[n_{2},n_{2}+30]$ by drawing edges from the Bernoulli distribution
with probability $p$ (we perform 4500 trials for each $p$ value in [.125, .25,
.375, .5, .625, .75, .875]). We compute the distance $\tilde{D}(G_{i},G_{k})$
(for $(i,k)\in\\{(1,3),(1,2),(2,3)\\}$). The results may be seen in Figure 5.
In this figure we plot a histogram of the “discrepancy score”

$\displaystyle\text{Disc}(G_{1},G_{2},G_{3})=\tilde{D}(G_{1},G_{3})/(\tilde{D}(G_{1},G_{2})+\tilde{D}(G_{2},G_{3})),$
(66)

which measures the degree to which a triplet of graphs violates the triangle
inequality (i.e. falls outside of the unit interval [0,1]), for approximately
$3\times 10^{4}$ such triplets. It is clear that, especially for the linear
definition of the distance, the triangle inequality is not always satisfied.
However, we also observe that (for graphs of these sizes) the discrepancy
score is bounded: no triple violates the triangle inequality by more than a
factor of approximately 1.8. This is shown by the histogram of discrepancies
in Figure 5. Additionally, the triangle inequality is satisfied in 28184
(95.2%) of cases.

We see similar but even stronger results when we run the same experiment with
$D^{2}$ instead of $\tilde{D}^{2}$; these may also be seen in Figure 5. We
calculated the discrepancy score analogously, but with $D$ substituted for
$\tilde{D}$. We see similarly that the degree of violation is bounded. In this
case, no triple violated the triangle inequality by a factor of more than 5,
and in this case the triangle inequality was satisfied in 99.8% of the
triples. In both of these cases, the triangle inequality violations may be a
result of our optimization procedure finding local minima/maxima for one or
more of the three distances computed. We also repeat the above procedure for
the same triplets of graphs, but with distances computed not in order of
increasing vertex size: calculating $\text{Disc}(G_{2},G_{1},G_{3})$ and
$\text{Disc}(G_{3},G_{2},G_{1})$. All of these results are plotted in Figure
5.

Figure 5: Histograms of triangle inequality violation. These plots show the
distribution of $\text{Disc}(G_{1},G_{2},G_{3})$, as defined in the text, for
the cases (a) top: the linear or small-time version of distance and (b)
bottom: the exponential or arbitrary-time version of distance. We see that for
the sizes of graph we consider, the largest violation of the triangle
inequality is bounded, suggesting that our distance measure may be an
infra-$\rho$-pseudometric for some value of $\rho\approx 1.8$ (linear version)
or $\rho\approx 5.0$ (exponential version). See Table 1 for a summary of the
distance metric variants introduced in this paper. We also plot the same
histogram for out-of-order (by vertex size) graph sequences:
$\text{Disc}(G_{2},G_{1},G_{3})$ and $\text{Disc}(G_{3},G_{2},G_{1})$. Each
plot has a line at $x=1$, the maximum discrepancy score for which the
underlying distances satisfy the triangle inequality.

####  5.3.2 Intra- and Inter-Lineage Distances

We compute pairwise distances for sequences of graphs in the graph lineages
displayed in Figure 3. For each pair of graph families (Square Grids, Paths,
Cycles, and Multi-Barbells), we compute the distance from the $i$th member of
one lineage to the $(i+1)$-st member of each other lineage, and take the
average of the resulting distances from $i=1$ to $i=12$. These distances are
listed in Table 2. As expected, average distances within a lineage are smaller
than the distances from one lineage to another.

| Square Grids | Paths | Cycles | Multi-Barbells
---|---|---|---|---
Square Grids | 0.0096700 | 0.048162 | 0.046841 | 0.63429
Paths | 0.30256 | 0.0018735 | 0.010300 | 2.1483
Cycles | 0.27150 | 0.0083606 | 0.0060738 | 2.0357
Multi-Barbells | 0.21666 | 0.75212 | 0.72697 | 0.029317
Table 2: Mean distances between graphs in several lineages. For two lineages
$G_{1},G_{2}\ldots$ (listed at left) and $H_{!},H_{2},\ldots$ (listed at the
top), each entry shows the mean distance $D(G_{i},H_{i+1})$ (where the average
is taken over $i=1$ to $12$). As expected, we see that the distance from
elements of a graph lineage to other members of the same lineage (the diagonal
entries of the table) is smaller than distances taken between lineages.
Furthermore as expected, 1D paths are more similar (but not equal) to 1D
cycles than to other graph lineages.

Figure 6: Distances $D^{2}(G,H)$ calculated for several pairs of graphs. The
top plot shows distances where $G$ and $H$ are both chosen from
$\\{\text{Grid}_{13\times 13},P_{169},C_{169},\text{Ba}_{13}\\}$. At bottom,
distances are calculated from $G$ chosen in $\\{\text{Grid}_{12\times
12},P_{144},C_{144},\text{Ba}_{12}\\}$ to H chosen in
$\\{\text{Grid}_{13\times 13},P_{169},C_{169},\text{Ba}_{13}\\}$. As expected,
diagonal entries are smallest.

####  5.3.3 Graph Limits

Here, we provide preliminary evidence that graph distance measures of this
type may be used in the definition of a _graph limit_ \- a graphlike object
which is the limit of an infinite sequence of graphs. This idea has been
previously explored, most famously by Lovász lovasz2012large , whose
definition of a graph limit (called a _graphon_) is as follows: Recall the
definition of graph cut-distance $D_{\text{cut}}(G,H)$ from Equation 3,
namely: the cut distance is the maximum discrepancy in sizes of edge-cuts,
taken over all possible subsets of vertices, between two graphs on the same
vertex-set. A graphon is then an equivalence class of Cauchy sequences of
graphs111Here we are calling a sequence of graphs “Cauchy” if for any
$\epsilon>0$ there is some $N$ such that for all $n,m\geq N$,
$D_{\text{cut}}(G_{n},G_{m})<\epsilon$. , under the equivalence relation that
two sequences $G_{1},G_{2},\ldots$ and $H_{1},H_{2},\ldots$ are equivalent if
$D_{\text{cut}}(G_{i},H_{i})$ approaches 0 as $n\rightarrow\infty$.

We propose a similar definition of graph limits, but with our diffusion
distance substituted as the distance measure used in the definition of a
Cauchy sequence of graphs. Hammond et. al. argue in hammond2013graph  why
their variant of diffusion distance may be a more descriptive distance measure
than cut-distance. More specifically, they show that on some classes of
graphs, some edge deletions ‘matter’ much more than others: removal of a
single edge changes the diffusive properties of the graph significantly.
However, the graph-cut distance between the new and old graphs is the same,
regardless of which edge has been removed, while the diffusion distance
captures this nuance. For graph limits, however, our generalization to
_unequal-sized graphs_ via $P$ is of course essential.

We examine several sequences of graphs of increasing size for the required
Cauchy behavior (in terms of our distance measure) to justify defining “graph
limits” in this way. For each of the graph sequences defined in Section 5.1,
we examine the distance between successive members of the sequence, plotting
$D^{2}(G_{n},H_{n+1})$ for each choice of $G$ and $H$. These sequences of
distances are plotted in Figure 9.

In this figure, we see that generally distance diverges between different
graph lineages, and converges for successive members of the same lineage, as
$n\rightarrow\infty$. We note the exceptions to this trend:

  1. 1.

The distances between $n$-paths and $n+1$-cycles appear to be converging; this
is intuitive, as we would expect that difference between the two spectra due
to distortion from the ends of the path graph would decrease in effect as
$n\rightarrow\infty$.

  2. 2.

We also show analytically, under similar assumtions, that the distance between
successive path graphs also shrinks to 0 (Theorem 5.1).

We do not show that all similarly-constructed graph sequences display this
Cauchy-like behavior. We hope to address this deeper question with one or more
modified versions of the objective function (see Section 6).

In the following theorem (Theorem 5.1), we assume that the optimal value of
$t$ approaches some value $\tilde{t}$ as $n\rightarrow\infty$. We have not
proven this to be the case, but have observed this behavior for both square
grids and path graphs. Lemmas 5 and 5.1 show a related result for path graphs;
we note that the spectrum of the Laplacian (as we define it in this paper) of
a path graph of size $n$ is given by

$\displaystyle\lambda_{k}=-2+2\cos{\frac{k\pi}{n-1}}\quad\quad
k\in\\{0...n-1\\}.$ (67)

.

######  Lemma 5.

For any finite $k,t$, we have

$\displaystyle\lim_{n\rightarrow\infty}n{\left(e^{t(-2+2\cos(\frac{\pi
k}{n}))}-e^{t(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}=0$ (68)

###### Proof.

Clearly for finite $k,t$

$\displaystyle\lim_{n\rightarrow\infty}{\left(e^{t(-2+2\cos(\frac{\pi
k}{n}))}-e^{t(-2+2\cos(\frac{\pi k}{n+1}))}\right)}=0$ (69)

Then,

$\displaystyle\lim_{n\rightarrow\infty}n{\left(e^{-2+2\cos(\frac{\pi
k}{n})}-e^{-2+2\cos(\frac{\pi k}{n+1})}\right)}$ (70)
$\displaystyle\quad\quad\quad=\lim_{n\rightarrow\infty}\frac{{\left(e^{-2+2\cos(\frac{\pi
k}{n})}-e^{-2+2\cos(\frac{\pi k}{n+1})}\right)}}{\frac{1}{n}}$ (71)

Evaluating this expression requires applying L’Hôpital’s rule. Hence, we have:

$\displaystyle\lim_{n\rightarrow\infty}\frac{{\left(e^{-2+2\cos(\frac{\pi
k}{n})}-e^{-2+2\cos(\frac{\pi k}{n+1})}\right)}}{\frac{1}{n}}$ (72)
$\displaystyle\quad\quad=\lim_{n\rightarrow\infty}\frac{2\pi
kt\left(\frac{\sin\left(\frac{\pi k}{n}\right)e^{2t\left(\cos\left(\frac{\pi
k}{n}\right)-1\right)}}{n^{2}}-\frac{\sin\left(\frac{\pi
k}{n+1}\right)e^{2t\left(\cos\left(\frac{\pi
k}{n+1}\right)-1\right)}}{(n+1)^{2}}\right)}{\frac{-1}{n^{2}}}$ (73)
$\displaystyle\quad\quad=2\pi
kt\lim_{n\rightarrow\infty}\left(\frac{n^{2}\sin\left(\frac{\pi
k}{n+1}\right)e^{2t\left(\cos\left(\frac{\pi
k}{n+1}\right)-1\right)}}{(n+1)^{2}}-\sin\left(\frac{\pi
k}{n}\right)e^{2t\left(\cos\left(\frac{\pi k}{n}\right)-1\right)}\right).$
(74)

Since both of the limits

$\displaystyle\lim_{n\rightarrow\infty}\left(\frac{n^{2}\sin\left(\frac{\pi
k}{n+1}\right)e^{2t\left(\cos\left(\frac{\pi
k}{n+1}\right)-1\right)}}{(n+1)^{2}}\right)$ (75) and
$\displaystyle\lim_{n\rightarrow\infty}\left(-\sin\left(\frac{\pi
k}{n}\right)e^{2t\left(\cos\left(\frac{\pi k}{n}\right)-1\right)}\right)$ (77)

exist (and are 0),

$\displaystyle 2\pi
kt\lim_{n\rightarrow\infty}\left(\frac{n^{2}\sin\left(\frac{\pi
k}{n+1}\right)e^{2t\left(\cos\left(\frac{\pi
k}{n+1}\right)-1\right)}}{(n+1)^{2}}-\sin\left(\frac{\pi
k}{n}\right)e^{2t\left(\cos\left(\frac{\pi k}{n}\right)-1\right)}\right)=0$
(78) and therefore
$\displaystyle\lim_{n\rightarrow\infty}n{\left(e^{t(-2+2\cos(\frac{\pi
k}{n}))}-e^{t(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}=0$ (79)

∎

######  Theorem 5.1.

If
$\lim_{n\rightarrow\infty}\arg\sup_{t}D^{2}\left(\left.\text{\emph{Pa}}_{n},\text{\emph{Pa}}_{n+1}\right|t\right)$
exists, then:

$\displaystyle\lim_{n\rightarrow\infty}D^{2}\left(\text{Pa}_{n},\text{Pa}_{n+1}\right)=0.$
(80)

###### Proof.

Assume that
$\lim_{n\rightarrow\infty}\arg\sup_{t}D^{2}\left(\left.\text{Pa}_{n},\text{Pa}_{n+1}\right|t\right)=\tilde{t}$.
Then, we must have

$\displaystyle\lim_{n\rightarrow\infty}D^{2}\left(\text{Pa}_{n},\text{Pa}_{n+1}\right)$
$\displaystyle\leq\lim_{n\rightarrow\infty}D^{2}\left(\left.\text{Pa}_{n},\text{Pa}_{n+1}\right|\tilde{t}\right)$
(81)

Hence, it remains only to prove that

$\displaystyle\lim_{n\rightarrow\infty}D^{2}\left(\left.\text{Pa}_{n},\text{Pa}_{n+1}\right|t\right)=0$
(82)

for any finite $t$ (which will then include $\tilde{t}$). First, for any
particular $(n+1)\times n$ subpermutation matrix $S$, note that

$\displaystyle D^{2}\left(\left.\text{Pa}_{n},\text{Pa}_{n+1}\right|t\right)$
$\displaystyle=\inf_{\alpha>0}\inf_{P|\mathcal{C}(P)}D^{2}\left(\left.\text{Pa}_{n},\text{Pa}_{n+1}\right|t,P,\alpha\right)$
(83) $\displaystyle\leq
D^{2}\left(\left.\text{Pa}_{n},\text{Pa}_{n+1}\right|t,\alpha=1,S\right)$ (84)

If at each $n$ we select $S$ to be the subpermutation
$S=\left[\begin{array}[]{c}I\\\ 0\end{array}\right]$, then the objective
function simplifies to:

$\displaystyle D^{2}\left(\text{Pa}_{n},\text{Pa}_{n+1}\right|t,P=S,\alpha=1)$
(85)
$\displaystyle\quad={\left|\left|Se^{c\Lambda_{\text{Pa}_{n}}}-e^{c\Lambda_{\text{Pa}_{n+1}}}S\right|\right|}_{F}^{2}$
(86) $\displaystyle\quad=\sum_{k=0}^{n-1}{\left(e^{c(-2+2\cos(\frac{\pi
k}{n}))}-e^{c(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}$ (87)
$\displaystyle\quad\leq\max_{0\leq k\leq n-1}n{\left(e^{c(-2+2\cos(\frac{\pi
k}{n}))}-e^{c(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}$ (88)

By Lemma 5, for any finite $k,t$, we have

$\displaystyle\lim_{n\rightarrow\infty}n{\left(e^{t(-2+2\cos(\frac{\pi
k}{n}))}-e^{t(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}=0$ (89)

So for any $\epsilon>0$, $\exists N$ such that when $n\geq N$, for any $c,k$,

$\displaystyle n{\left(e^{c(-2+2\cos(\frac{\pi
k}{n}))}-e^{c(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}<\epsilon$ (90)

But then

$\displaystyle\sum_{k=0}^{n-1}{\left(e^{c(-2+2\cos(\frac{\pi
k}{n}))}-e^{c(-2+2\cos(\frac{\pi k}{n+1}))}\right)}^{2}<\epsilon$ (91)

as required. Thus, the Cauchy condition is satisfied for the lineage of path
graphs $\text{Pa}_{n}$ ∎

Figure 7: Limiting behavior of $D$ and two parameters as path graph size
approaches infinity. All distances were calculated between $\text{Path}_{n}$
and $\text{Path}_{n+1}$. We plot the value of the objective function, as well
as the optimal values of $\alpha$ and $t$, as $n\rightarrow\infty$. Optimal
$\alpha$ rapidly approach 1 and the optimal distance tends to 0. Additionally,
the optimal $t$ value approaches a constant $(t\approx.316345)$, providing
experimental validation of the assumption we make in proving Theorem 5.1.

Given a graph lineage which consists of levelwise box products between two
lineages, it seems natural to use our upper bound on successive distances
between graph box products to prove convergence of the sequence of products.
As an example, the lineage consisting of square grids is the levelwise box
product of the lineage of path graphs with itself. However, in this we see
that this bound may not be very tight. Applying Equation (51) from Theorem
3.3, we have (for any $t_{c},\alpha_{c}$):

$\displaystyle D\left(\left.\text{Sq}_{n},\text{Sq}_{n+1}\right.\right)$
$\displaystyle\leq
D\left(\left.\text{Sq}_{n},\text{Sq}_{n+1}\right|t_{c},\alpha_{c}\right)$ (92)
$\displaystyle\leq
D\left(\text{Pa}_{n+1},\text{Pa}_{n+1}\left.\right|t_{c},\alpha_{c}\right)\left({\left|\left|e^{\frac{t_{c}}{a_{c}}L(\text{Pa}_{n})}\right|\right|}_{F}\right.$
$\displaystyle\quad\quad\quad\quad+\left.{\left|\left|e^{t_{c}a_{c}L(\text{Pa}_{n+1})}\right|\right|}_{F}\right)$

As we can see in Figure 8, the right side of this inequality seems to be
tending to a nonzero value as $n\rightarrow\infty$, whereas the actual
distance (calculated by our optimization procedure) appears to be tending to
0.

Figure 8: Comparison of the distance $D(\text{Sq}_{n},\text{Sq}_{n+1})$ as a
function of $n$, to the upper bound calculated as the optimum of distance
between $\text{Pa}_{n}$ and $\text{Pa}_{n+1}$. We see that the upper found
converges to some constant $D\approx 0.01782$, whereas the actual distance
appears to be converging to 0 as $n\rightarrow\infty$. Figure 9: Cauchy-like
behavior of graph distance as a function of sequence index, $n$. The distance
between successive square grids and all other graph sequences appears to
diverge (the same behavior is seen for k-barbells). Notably, the distance
between $\text{Grid}_{n\times n}$ and $\text{Grid}_{(n+1)\times(n+1)}$ does
not appear to converge, until much higher values of $n$ $(n>100)$ than the
other convergent series. This may be because the distances calculated are an
upper bound, and may be converging more slowly than the ‘true’ optima.

##  6 Regularized Distance

We can add a regularization term to the graph diffusion distance, as follows:
define

$\displaystyle D_{\text{reg}}(G_{1},G_{2})$
$\displaystyle=\sup_{t}\inf_{P|\mathcal{C}(P)}\inf_{\alpha>0}\left\\{{\left|\left|Pe^{\frac{t}{\alpha}L_{1}}-e^{t\alpha
L_{2}}P\right|\right|}_{F}+{\left|\left|e^{\frac{t}{\alpha}L_{1}}-e^{tL_{1}}\right|\right|}_{F}\right.$
$\displaystyle\quad\quad\quad+\left.{\left|\left|e^{tL_{2}}P-e^{t\alpha
L_{2}}P\right|\right|}_{F}\right\\}$ (93)

We can show analytically that this distance satisfies the triangle inequality:

######  Theorem 6.1.

$D_{\text{reg}}$ satisfies the triangle inequality.

###### Proof.

For graphs $G_{1},G_{2},G_{3}$ and Laplacians $L_{1},L_{2},L_{3}$, for any
fixed $t\geq 0$, we have:

$\displaystyle D_{\text{reg}}(G_{1},G_{3}|t)$
$\displaystyle=\inf_{P|\mathcal{C}(P)}\inf_{\alpha>0}\left\\{{\left|\left|Pe^{\frac{t}{\alpha}L_{1}}-e^{t\alpha
L_{3}}P\right|\right|}_{F}+{\left|\left|e^{\frac{t}{\alpha}L_{1}}-e^{tL_{1}}\right|\right|}_{F}\right.$
$\displaystyle\quad\quad\quad+\left.{\left|\left|e^{tL_{3}}P-e^{t\alpha
L_{3}}P\right|\right|}_{F}\right\\}$ (94) $\displaystyle\leq
D_{\text{reg}}(G_{1},G_{3}|t,\alpha=1)$ (95)
$\displaystyle=\inf_{P|\mathcal{C}(P)}\left\\{{\left|\left|Pe^{tL_{1}}-e^{tL_{3}}P\right|\right|}_{F}+{\left|\left|e^{tL_{1}}-e^{tL_{1}}\right|\right|}_{F}\right.$
$\displaystyle\quad\quad\quad+\left.{\left|\left|e^{tL_{3}}P-e^{tL_{3}}P\right|\right|}_{F}\right\\}$
(96)
$\displaystyle=\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{tL_{1}}-e^{tL_{3}}P\right|\right|}_{F}$
(97)

Suppose that

$\displaystyle\alpha_{32},P_{32}$
$\displaystyle=\arg\inf_{a>0}\inf_{P|\mathcal{C}(P)}\left\\{{\left|\left|Pe^{\frac{t}{\alpha}L_{2}}-e^{t\alpha
L_{3}}P\right|\right|}_{F}+{\left|\left|e^{\frac{t}{\alpha}L_{2}}-e^{tL_{2}}\right|\right|}_{F}\right.$
$\displaystyle\quad\quad\quad+\left.{\left|\left|e^{tL_{3}}P-e^{t\alpha
L_{3}}P\right|\right|}_{F}\right\\}$ (98) $\displaystyle\alpha_{21},P_{21}$
$\displaystyle=\arg\inf_{a>0}\inf_{P|\mathcal{C}(P)}\left\\{{\left|\left|Pe^{\frac{t}{\alpha}L_{1}}-e^{t\alpha
L_{2}}P\right|\right|}_{F}+{\left|\left|e^{\frac{t}{\alpha}L_{1}}-e^{tL_{1}}\right|\right|}_{F}\right.$
$\displaystyle\quad\quad\quad+\left.{\left|\left|e^{tL_{2}}P-e^{t\alpha
L_{2}}P\right|\right|}_{F}\right\\}$ (99)

Then,

$\displaystyle\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{tL_{1}}-e^{tL_{3}}P\right|\right|}_{F}$
$\displaystyle\leq{\left|\left|P_{32}P_{21}e^{tL_{1}}-e^{tL_{3}}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{tL_{1}}-e^{tL_{3}}P\right|\right|}_{F}$
$\displaystyle\leq{\left|\left|P_{32}P_{21}e^{tL_{1}}-P_{32}P_{21}e^{\frac{t}{\alpha_{21}}L_{1}}\right.\right.}+P_{32}P_{21}e^{\frac{t}{\alpha_{21}}L_{1}}$
$\displaystyle\quad-
P_{32}e^{t\alpha_{21}L_{2}}P_{21}+P_{32}e^{t\alpha_{21}L_{2}}P_{21}-P_{32}e^{tL_{2}}P_{21}$
$\displaystyle\quad+P_{32}e^{tL_{2}}P_{21}-P_{32}e^{\frac{t}{\alpha_{32}}L_{2}}P_{21}+P_{32}e^{\frac{t}{\alpha_{32}}L_{2}}P_{21}$
$\displaystyle\quad-e^{t\alpha_{32}L_{3}}P_{32}P_{21}{\left.\left.+e^{t\alpha_{32}L_{3}}P_{32}P_{21}-e^{tL_{3}}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle\leq{\left|\left|P_{32}P_{21}e^{tL_{1}}-P_{32}P_{21}e^{\frac{t}{\alpha_{21}}L_{1}}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|P_{32}P_{21}e^{\frac{t}{\alpha_{21}}L_{1}}-P_{32}e^{t\alpha_{21}L_{2}}P_{21}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|P_{32}e^{t\alpha_{21}L_{2}}P_{21}-P_{32}e^{tL_{2}}P_{21}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|P_{32}e^{tL_{2}}P_{21}-P_{32}e^{\frac{t}{\alpha_{32}}L_{2}}P_{21}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|P_{32}e^{\frac{t}{\alpha_{32}}L_{2}}P_{21}-e^{t\alpha_{32}L_{3}}P_{32}P_{21}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|e^{t\alpha_{32}L_{3}}P_{32}P_{21}-e^{tL_{3}}P_{32}P_{21}\right|\right|}_{F}$
by Lemma 1,
$\displaystyle\inf_{P|\mathcal{C}(P)}{\left|\left|Pe^{tL_{1}}-e^{tL_{3}}P\right|\right|}_{F}$
$\displaystyle\leq{\left|\left|e^{tL_{1}}-e^{\frac{t}{\alpha_{21}}L_{1}}\right|\right|}_{F}+{\left|\left|P_{21}e^{\frac{t}{\alpha_{21}}L_{1}}-e^{t\alpha_{21}L_{2}}P_{21}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|e^{t\alpha_{21}L_{2}}P_{21}-e^{tL_{2}}P_{21}\right|\right|}_{F}$
$\displaystyle+{\left|\left|e^{tL_{2}}-e^{\frac{t}{\alpha_{32}}L_{2}}\right|\right|}_{F}+{\left|\left|P_{32}e^{\frac{t}{\alpha_{32}}L_{2}}-e^{t\alpha_{32}L_{3}}P_{32}\right|\right|}_{F}$
$\displaystyle\quad+{\left|\left|e^{t\alpha_{32}L_{3}}P_{32}-e^{tL_{3}}P_{32}\right|\right|}_{F}$
$\displaystyle=D_{\text{reg}}(G_{1},G_{2}|t=c)+D_{\text{reg}}(G_{2},G_{3}|t=c)$

Since this is true for any fixed $t$, let

$\displaystyle t^{*}=\arg\sup_{t}D_{\text{reg}}(G_{1},G_{3}|t)$

. Then

$\displaystyle D_{\text{reg}}(G_{1},G_{3})$
$\displaystyle=\sup_{c}D_{\text{reg}}(G_{1},G_{3}|t)$
$\displaystyle=D_{\text{reg}}(G_{1},G_{3}|t^{*})$ $\displaystyle\leq
D_{\text{reg}}(G_{1},G_{2}|t^{*})+D_{\text{reg}}(G_{2},G_{3}|t=t^{*})$
$\displaystyle\leq\sup_{t_{21}}D_{\text{reg}}(G_{1},G_{2}|t_{21})+\sup_{t_{32}}D_{\text{reg}}(G_{2},G_{3}|t_{32})$
$\displaystyle=D_{\text{reg}}(G_{1},G_{2})+D_{\text{reg}}(G_{2},G_{3})$

∎

We can construct a similar regularized version of the linear objective
function:

$\displaystyle\tilde{D}_{\text{reg}}(G_{1},G_{2})$
$\displaystyle={\left|\left|\frac{1}{\alpha}PL_{1}-\alpha
L_{2}P\right|\right|}+{\left|\left|\frac{1}{\alpha}L_{1}-L_{1}\right|\right|}+{\left|\left|\frac{\hbox{}}{\hbox{}}PL_{2}-\alpha
L_{2}P\right|\right|}$

The term “regularized” here refers to the fact that the additional terms
included in $D_{\text{reg}}$ and $\tilde{D}_{\text{reg}}$ penalize $\alpha$
distorting the respective Laplacians far from their original values. In
practice, many of the theoretical guarantees provided earlier in this
manuscript may not apply to optimization of the augmented objective function.
Hence, a major area of future work will be modification of our optimization
procedure to compute this form of distance.

##  7 Applications and Future Work

We briefly discuss possible applications of both our distance metric and our
procedure for calculating the relevant minima.

###  7.1 Algebraic Multigrid

The need for prolongation/restriction operators arises naturally in the
Algebraic MultiGrid (AMG) context, where a hierarchy of progressively coarser
meshes are constructed, with the goal of speeding convergence of a model with
local update (“smoothing”) rules. A model with modes of behavior at
wavelengths which are much larger than the neighborhood of one update will
take many update steps to converge. Thus, the goal in AMG is to iteratively
construct a series of coarsened meshes, so that update steps at the coarser
scales can address coarser modes of behavior. A fine-scale model state is
translated into a coarse-scale state via a “restriction” operator. After a
coarse-scale smoothing step, the new coarse state is translated back to the
fine-scale by “prolonging” it. Our procedure for calculating $P$ could be
incorporated as a preprocessing step, in the case where the series of meshes
are known in advance; otherwise, the $P$ from the previous round of coarsening
could be used as the initial conditions to a modified version of our solver.
In either case, the matrix $P$ is a natrual choice of prolongation/restriction
operator for this type of coarsening scheme, since it optimally transforms the
Laplacian of one graph into another.

###  7.2 Graph Limits

In this work we briefly introduce a new definition of graph limits based on
the diffusion distance, which raises several natural questions: What does the
“limit” of a sequence of graphs under diffusion distance look like? Are there
pairs of sequences that converge to the same such object, as in the example of
path graphs and cycle graphs? Can we separate graph sequences into equivalence
classes based on which of these they converge to? We hope to address these
questions in future work.

###  7.3 Graph Convolutional Networks

Graph convolutional networks (GCNs) are a variant of the _convolutional neural
networks_ (CNNs) widely used in machine vision. In the same way that CNNs
learn a set of trained image filters and apply them across multiple spatial
locations in an image, GCNs learn a set of filters which are applied to local
neighborhoods of a graph. One implementation of GCNs due to Kipf and Welling
kipf2016semi  uses a Chebyshev polynomial of the Laplacian matrix as an
approximation of the graph fourier transform, demonstrating comparable results
to the full transform but far fewer multiplication operation needed. However,
construction of pooling operators for GCNs is still an area of open research.
Since our $P$ is a restriction operator that preserves information about the
Laplacian, it is natural to use it as a pooling operator in this type of
model.

###  7.4 Graph Clustering

We can also use the diffusion distance and its variants to compare graphs or
neighborhoods of graphs for structural similarity, independent of graph size.
This is similar to the approach of gold1995new  for comparing point clouds in
2D and 3D, in the sense that both approaches optimize an objective function
based to a matching between elements of the two graphs. This type of
similarity measure may then be used to convert a dataset of graphs to a
distance-to-cluster-centers representation, or for any other of the typical
methods used in machine learning for converting sets of pairwise distances
into fixed-length feature vectors ($k$-medoids, kernel methods,
multidimensional scaling, etc.). In this setting, our distance measure has an
additional benefit: since computing it yields an explicit projection operator
between the nodes of the graphs, we may use the set of $P$ we compute to
project signals (e.g. labels on the vertices of each graph in the dataset) to
a common space.

##  8 Conclusion

In this work, we present a novel generalization of graph diffusion distance
which allows for comparison of graphs of inequal size. We consider several
variants of this distance measure to account for sparse maps between the two
graphs, and for maps between the two graphs which are optimal given a fixed
time-dilation factor $\alpha$. We prove several important theory properties of
distances in this family of measures, including triangle inequalities in some
cases and Cauchy-like behavior of some graph sequences. We present a new
procedure for optimizing the objective function defined by our distance
measure, prove the correctness of this procedure, and demonstrate its
efficiency in comparison to univariate search over the dilation parameter,
$\alpha$. Numerical experiments suggest that this dissimilarity score
satisfies the triangle inequality up to some constant $\rho\approx 2.1$. We
demonstrate that this measure of graph distance may be used to compare graph
lineages (families of exponentially-growing graphs with shared structure), and
additionally that certain lineages display Cauchy-sequence like behavior as
the graph size approaches infinity. We suggest several possible applications
of our distance measure to scientific problems in the contexts of pattern
matching and machine learning.

##  9 Acknowledgements

This work was supported by U.S. National Science Foundation NRT Award number
1633631, U.S. National Institute of Aging grant AG059602, Human Frontiers
Science Program grant HFSP - RGP0023/ 2018, U.S. National Institutes for
Health grant R01HD073179, USAF/DARPA FA8750-14-C-0011, and by the Leverhulme
Trust and and the hospitality of the Sainsbury Laboratory at Cambridge
University, as well as the hospitality of the Center for Non-Linear Studies
(CNLS) at Los Alamos National Laboratory.

## References

  * (1) Allgower, E.L., Georg, K.: Numerical continuation methods: an introduction, vol. 13.  Springer Science & Business Media (2012)
  * (2) Bentley, J.L.: Multidimensional binary search trees used for associative searching.  Communications of the ACM 18(9), 509–517 (1975)
  * (3) Bijsterbosch, J., Volgenant, A.: Solving the Rectangular Assignment Problem and Applications.  Annals of Operations Research 181(1), 443–462 (2010)
  * (4) Brouwer, A.E., Haemers, W.H.: Spectra of Graphs.  Springer Science & Business Media (2011)
  * (5) Brouwer, A.E., Spence, E.: Cospectral Graphs on 12 Vertices.  The Electronic Journal of Combinatorics 16(1), N20 (2009)
  * (6) Burkard, R.E., Cela, E.: Linear Assignment Problems and Extensions.  In: Handbook of Combinatorial Optimization, pp. 75–149. Springer (1999)
  * (7) Burkard, R.E., Dell’Amico, M., Martello, S.: Assignment Problems.  Springer (2009)
  * (8) Croes, G.A.: A method for solving traveling-salesman problems.  Operations research 6(6), 791–812 (1958)
  * (9) Dodonova, Y., Belyaev, M., Tkachev, A., Petrov, D., Zhukov, L.: Kernel Classification of Connectomes Based on Earth Mover’s Distance between Graph Spectra.  arXiv preprint arXiv:1611.08812 (2016)
  * (10) Donnat, C., Holmes, S.: Tracking Network Dynamics: A Survey of Distances and Similarity Metrics.  arXiv preprint arXiv:1801.07351 (2018)
  * (11) Eshera, M., Fu, K.S.: A Graph Distance Measure for Image Analysis.  IEEE Transactions on Systems, Man, and Cybernetics SMC-14(3), 398–408 (1984)
  * (12) Eshera, M., Fu, K.S.: An Image Understanding System using Attributed Symbolic Representation and Inexact Graph-Matching.  IEEE Transactions on Pattern Analysis and Machine Intelligence PAMI-8(5), 604–618 (1986)
  * (13) Gao, X., Xiao, B., Tao, D., Li, X.: A Survey of Graph Edit Distance.  Pattern Analysis and Applications 13(1), 113–129 (2010)
  * (14) Godsil, C.D., McKay, B.: Constructing Cospectral Graphs.  Aequationes Mathematicae 25(1), 257–268 (1982)
  * (15) Gold, S., Lu, C.P., Rangarajan, A., Pappu, S., Mjolsness, E.: New Algorithms for 2d and 3d Point Matching: Pose Estimation and Correspondence.  In: Advances in neural information processing systems, pp. 957–964 (1995)
  * (16) Gold, S., Rangarajan, A., Mjolsness, E.: Learning with preknowledge: Clustering with Point and Graph Matching Distance Measures.  In: Advances in Neural Information Processing Systems, pp. 713–720 (1995)
  * (17) Gu, J., Hua, B., Liu, S.: Spectral Distances on Graphs.  Discrete Applied Mathematics 190, 56–74 (2015)
  * (18) Hagberg, A., Swart, P., S Chult, D.: Exploring Network Structure, Dynamics, and Function using NetworkX.  Tech. rep., Los Alamos National Lab.(LANL), Los Alamos, NM (United States) (2008)
  * (19) Hammond, D.K., Gur, Y., Johnson, C.R.: Graph Diffusion Distance: A Difference Measure for Weighted Graphs based on the Graph Laplacian Exponential Kernel.  In: 2013 IEEE Global Conference on Signal and Information Processing, pp. 419–422. IEEE (2013)
  * (20) Heindl, C.: lapsolver: Fast Linear Assignment Problem (LAP) Solvers for Python Based on c-extensions.  https://github.com/cheind/py-lapsolver (2018)
  * (21) Hogben, L.: Handbook of Linear Algebra.  CRC Press (2006)
  * (22) Inc., W.R.: Mathematica, Version 11.3.  Champaign, IL, 2018
  * (23) Eric Jones, Travis Oliphant, Pearu Peterson, et al.: SciPy: Open Source Scientific Tools for Python (2001–).  URL http://www.scipy.org/.  [Online; accessed June 20 2019]
  * (24) Kipf, T.N., Welling, M.: Semi-supervised classification with graph convolutional networks.  arXiv preprint arXiv:1609.02907 (2016)
  * (25) Lovász, L.: Large Networks and Graph Limits, vol. 60.  American Mathematical Soc. (2012)
  * (26) McGregor, A., Stubbs, D.: Sketching Earth-Mover Distance on Graph Metrics.  In: Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pp. 274–286. Springer (2013)
  * (27) Onoyama, T., Kubota, S., Oyanagi, K., Tsuruta, S.: A Method for Solving Nested Combinatorial Optimization Problems-a Case of Optimizing a Large-Scale Distribution Network.  In: SMC 2000 Conference Proceedings. 2000 IEEE International Conference on Systems, Man and Cybernetics. ’Cybernetics Evolving to Systems, Humans, Organizations, and their Complex Interactions’ (cat. no. 0, vol. 1, pp. 340–345. IEEE (2000)
  * (28) Rangarajan, A., Gold, S., Mjolsness, E.: A Novel Optimizing Network Architecture with Applications.  Neural Computation 8(5), 1041–1060 (1996)
  * (29) Rossum, G.: Python Reference Manual.  Technical Report (1995)
  * (30) Sanfeliu, A., Fu, K.S.: A Distance Measure between Attributed Relational Graphs for Pattern Recognition.  IEEE Transactions on Systems, Man, and Cybernetics SMC-13(3), 353–362 (1983)
  * (31) Scott, C., Mjolsness, E.: Multilevel Artificial Neural Network Training for Spatially Correlated Learning (in press).  SIAM SISC Copper Mountain Special Issue (2019)
  * (32) Sinha, A., Malo, P., Deb, K.: A Review on Bilevel Optimization: from Classical to Evolutionary Approaches and Applications.  IEEE Transactions on Evolutionary Computation 22(2), 276–295 (2018)
  * (33) Townsend, J., Koep, N., Weichwald, S.: Pymanopt: A Python Toolbox for Optimization on Manifolds using Automatic Differentiation.  The Journal of Machine Learning Research 17(1), 4755–4759 (2016)
  * (34) Van Dam, E.R., Haemers, W.H.: Which Graphs are Determined by their Spectrum?  Linear Algebra and its applications 373, 241–272 (2003)
  * (35) Van Dam, E.R., Haemers, W.H.: Developments on Spectral Characterizations of Graphs.  Discrete Mathematics 309(3), 576–586 (2009)
  * (36) Van Der Walt, S., Colbert, S.C., Varoquaux, G.: The NumPy Array: a Structure for Efficient Numerical Computation.  Computing in Science & Engineering 13(2), 22 (2011)
  * (37) Wen, Z., Yin, W.: A Feasible Method for Optimization with Orthogonality Constraints.  Mathematical Programming 142(1-2), 397–434 (2013)

##  10 Appendices

###  10.1 Spectral Lower Bound

In Theorem 4.1 we derived an upper bound on the graph distance