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. Here we detail several measures of graph dissimilarity, chosen by historical significance and similarity to our measure.

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:

Does the distance measure require an inner optimization loop? If so is it mainly a discrete or continuous optimization formulation?

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₁)→S which identifies each vertex of a graph G₁ with some state s ∈ S, then a map-explicit graph distance is one whose calculation yields a second function g: V(G₂)→V(G₁), with the composite function f ∘ g: V(G₂)→S.

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 [3].

For each of the graph distance variants discussed here, we label them according to the above taxonomy. For example, the two prior works by Eschera et. al. and Hammond et al (discussed in Sections 1.1.4 and 1.1.5) would be labelled as (structural, explicit, disc-opt) and (spectral, implicit, non-opt), respectively. Our distance measure would be labelled (spectral, explicit, cont-opt).

1.1.1 Quadratic matching of points and graphs (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. [4] define the dissimilarity between two unlabelled weighted graphs (with adjacency matrices A⁽¹⁾ and A⁽²⁾ and n₁ and n₂ vertices, respectively) as the solution to the following optimization problem (for real-valued M = [m_ij]:

where $$ {\parallel ·\parallel }_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 [5, 6] 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.2 Cut-distance of graphs (structural, implicit, disc-opt)

Lovász [7] defines the cut-distance of a pair of graphs as follows: Let the □-norm of a matrix B be given by:

Given two labelled graphs G₁, G₂, on the same set of vertices, and their adjacency matrices A₁ and A₂, the cut-distance d_cut(G₁, G₂) is then given by

(for more details, see [7]). Computing this distance requires combinatorial optimization (over all vertex subsets of G₁, G₂) but this optimization does not result in an explicit map between G₁ and G₂.

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 ∈ S:

Non-negativity: d(x, y)≥0

Identity: d(x, y) = 0 ⇔ x = y

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

Triangle inequality: d(x, z)≤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₀, 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 + ∞. 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.3.

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.4.

1.3 Definitions

Graph Laplacian: For an undirected graph G with adjacency matrix A and vertex degrees d₁, d₂…d_n, we define the Laplacian of the graph as

Frobenius norm: The squared Frobenius norm, $$ {\parallel A\parallel }_F^2$$ of a matrix A is given by the sum of squares of matrix entries. This can equivalently be written as Tr[A^T A].

Linear Assignment Problem (LAP): We take the usual definition of the Linear Assignment Problem (see [17, 18]): we have two lists of items S and R (sometimes referred to as “workers” and “jobs”), and a cost function $$ c:S\times R\to \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:

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 [19] 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.

Hierarchical graph sequences: A Hierarchical Graph Sequence (HGS) is a sequence of graphs, indexed by $$ l\in \mathbb{N}=0,1,2,3\dots $$, satisfying the following:

G₀ 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\left(b^{l^{1+ϵ}}\right)$$, for all ϵ > 0.

Graded graph: A graded graph is a graph along with a vertex labelling, where vertices are labelled with non-negative integers such that Δl, the difference in label over any edge, is in {−1, 0, 1}. We will refer to the Δl = 0 edges as “within-level” and the l = ±1 edges as “between-level”.

Graph lineages: A graph lineage is a graded graph with two extra conditions:

The vertices and edges with Δl = 0 form a HGS; and

the vertices and edges with Δl = ±1 form a HGS of bipartite graphs.

More plainly, a graph lineage is an exponentially growing sequence of graphs along with ancestry relationships between nodes. We will also use the term graph lineage to refer to the HGS in the first part of the definition. Some intuitive examples of graph lineages in this sense are the following:

Path graphs or cycle graphs of size bⁿ for any integer b.

More generally, grid graphs of any dimension d, of side length b, yielding a lineage which grows with size b^dn (with periodic or nonperiodic boundary conditions).

For any probability distribution p(x, y) whose support is points in the unit square, we can construct a graph by discretizing the map of p as a function of x and y, and interpreting the resulting matrix as the adjacency matrix of a graph. For a specific probability distribution p, the graphs derived this way with discretizations of exponentially increasing bin count form a graph lineage.

The triangulated mesh is a common object in computer graphics [20–22], representing a discretization of a 2-manifold embedded in R³. Finer and finer subdivisions of such a mesh constitute a graph lineage.

Several examples of graph lineages are used in the discussion of the numerical properties of Graph Diffusion Distance in Section 5.1. Additional examples (a path graph and a triangulated mesh) can be found in Figs 1 and 2.

Fig 1

The first seven levels of the graph lineage of path graphs, with ancestry relationships.

Δl = 0 edges are colored in orange, Δl = ±1 edges are colored in blue. Self-loops are not illustrated.

Fig 2

Top: subsamples of a mesh of the Utah teapot, of increasing density (each node is connected to its 8 nearest neighbors by the Δl = ±0 edges, rendered in blue). These samples form a graph lineage (Δl = ±1 edges are not illustrated). Bottom: the same set of nodes, with only Δl = ±1 edges plotted (in orange) for one node from the coarsest level and its descendants.

Box product (□) of graphs: For two graphs G and H with vertex sets V(G) = {g₁, g₂…g_n} and V(H) = {h₁, h₂…h_m}, we say the product graph G□H is the graph with vertex set V(G□H) = V(G) × V(H) and an adjacency relationship defined by: (g₁, h₁)∼(g₂, h₂) in G□H if and only if g₁ ∼ g₂ in G and h₁ = h₂, or g₁ = g₂ and h₁ ∼ h₂ in H. Note that the adjacency matrix of this relationship may be represented by the following identity:

2.1 Diffusion distance definition

We generalize the diffusion distance defined by Hammond et al. [1]. This distortion measure between two graphs G₁ and G₂, of the same size, was defined as:

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^{t{L}_i}$$. Each column v_j of $$ e^{t{L}_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. See Fig 3 for an example of a distance calculation, with a characteristic peak. 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 stationary state 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₁ and G₂ are apparent in their differing probability distributions.

Fig 3

A plot illustrating unimodality of diffusion distance.

D² was calculated between two grid graphs Sq₇ and Sq₈ of size 7 × 7 and 8 × 8, respectively. The distance is given by the formula $$ D^2({\text{Sq}}_7,{\text{Sq}}_8|t)={inf}_{\alpha >0}{inf}_{P\left|\mathcal{C}\right(P)}{\parallel P{e}^{\frac{t}{a}L\left({\text{Sq}}_7\right)}-e^{t\alpha L\left({\text{Sq}}_8\right)}P\parallel }_F^2$$ as a function of t. The peak, at t ≈.318, yields the distance D²(Sq₇, Sq₈).

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, α. The dissimilarity between two graphs G₁ and G₂ (with Laplacians L_i = L(G_i)) is then defined as:

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₁ and nodes of G₂, so we are measuring how well diffusion on G₁ approximates diffusion on G₂, 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 |G₁| = n₁ ≤ n₂ = |G₂|; 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 3.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 6). 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 α. 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 α given t, and an innermost loop of minimization over P given t and α).

The other additional parameter, α, 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 × n grid is a reasonably accurate approximation of more rapid diffusion on a 2n × 2n grid, especially when n is very large. In our discussion of variants of this dissimilarity score, we will use the notation D²(G₁, G₂|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²(G₁, G₂|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 ≈ I + tL. This motivates the early-time or “linear” version of this distance, $$ \tilde{D}$$:

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

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₁, G₂:

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

2.2 Directedness of distance and constraints

We note that this distance measure, as defined so far, is directed: the operands G₁ and G₂ serve differing roles in the objective function. This additionally makes the constraint predicate $$ \mathcal{C}\left(P\right)$$ 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₁ and G₂ with |G₁| ≤ |G₂| is always D(G₁, G₂). 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.3 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 Eqs 8 and 9 over α and P is too loose, in the sense that the distances D(G₁, G₂), D(G₂, G₃), and D(G₁, G₃) 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)≥0 and d(x, y) = d(y, x) for all x, y.

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

ρ-inframetric: As a premetric, but additionally d(x, z)≤ρ(d(x, y) + d(y, z)) and d(x, y) = 0 if and only if x = y, for all x, y, z.

Table 1

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 α (fixed at α = 1, fixed at a constant power r of the ratio of graph sizes, or minimized over all α. We leave exploration of nonzero values of the sparsity parameter to future work. Variants not explicitly called out are not considered. In the case where α and t are both optimized and s > 0, it is unclear which of the metric conditions GDD satisfies, hence the corresponding classification is left blank.

Summary of this paper’s investigation of different forms of our graph dissimilarity measure.

t	α	s	Classification	Treatment in this manuscript
Fixed at tc < ϵ	Fixed at αc = 1	s = 0	Pseudometric	Defined in Eq 18. Optimized by one pass of LAP solver. Triangle inequality proven in Theorem 2.
Fixed at tc < ϵ	Fixed at $$	s = 0	Pseudometric	Defined in Eq (24). Optimized by one pass of LAP solver. Triangle inequality proven in Theorem 4.
Fixed at tc < ϵ	Optimized	s = 0	Premetric	Defined in Eq 9. Optimized by Algorithm 1. Triangle inequality violations examined experimentally in Section 5.3.1.
Optimized	Fixed at αc = 1	s = 0	Metric	When |G1| = |G2|, this is Hammond et. al’s version of graph distance.
Optimized	Optimized	s = 0	Premetric	Defined in Eq 8. Optimized by Algorithm 2. Graph Product upper bound proven in Theorem 5. 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 tc < ϵ	Fixed at αc = 1	s > 0	Pseudometric	Triangle inequality proven in Theorem 2.
Fixed at tc < ϵ	Fixed at $$	s > 0	Pseudometric	Triangle inequality proven in Theorem 4.
Optimized	Optimized	s > 0		Discussed in Section 3.3.

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 [24, 25]. 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 [26–28], in the sense that the probability of two graphs of size n being isospectral but not isomorphic goes to 0 as n → ∞. Furthermore, our numerical experiments have indicated that the violation of the triangle inequality is bounded, in the sense that D(G₁, G₃)≤ρ*(D(G₁, G₂) + D(G₂, G₃)) for ρ ≈ 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 Eqs (8) and (9). First, we discuss some theoretical properties of this dissimilarity measure.

3.1 Triangle inequality for α = 1

Lemma 1. For any matrices M and P, with P satisfying P^T P = I, $$ {\parallel PM\parallel }_F^2\le {\parallel M\parallel }_F^2$$ and $$ {\parallel MP\parallel }_F^2\le {\parallel M\parallel }_F^2$$.

Proof. Suppose without loss of generality that P^T P = I. Then:

$$ {\parallel PM\parallel }_F^2=Tr\left[M^T{P}^{T}PM\right]=Tr\left[M^{T}M\right]={\parallel M\parallel }_F^2$$

If P^T P = I, then letting PP^T = Π, Π is a projection operator satisfying Π^T = Π = Π². Then,

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

Proof. Let G₁, G₂, G₃ be simple graphs, with Laplacians L₁, L₂, L₃. Let

Then

By Lemma 1,

We note that in this proof we use L₁, L₂, and L₃ (making this the small-t or linear version of the objective function), but the same argument holds when all three are replaced with $$ e^{t{L}_i}$$, so we also have

Corollary 3. D satisfies the triangle inequality for α = 1.

Proof. By the same calculation as in Theorem 2, with all L_i replaced by $$ e^{t_c{L}_i}$$, we have

Note that in the proofs of Theorem 2, Theorem 4, and Corollary 3, we assume that the constraint predicate $$ \mathcal{C}\left(P\right)$$ 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.2 Time-Scaled Graph Diffusion Distance

For any graphs G₁ and G₂, we define the Time-Scaled Graph Diffusion Distance (TSGDD) as a scaled version of the linear distance, with α fixed. Namely, let

The intuition for this version of the distance measure is that we are constraining the time dilation, α, between G₁ and G₂ to be a power of the ratio of the two graph sizes. The factor (n₁ n₂)^−2r is needed to ensure this version of the distance satisfies the triangle inequality, as seen in Theorem 4.

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

Proof. As above, let G₁, G₂, G₃ be three graphs with n_i = |G_i| and n₁ ≤ n₂ ≤ n₃, and let L_i be the Laplacian of G_i. Let $$ \mathcal{C}\left(P\right)$$ represent a transitive constraint predicate, also as described previously. Then, for a constant $$ r\in \mathbb{R}$$, we have:

By Lemma 1,

3.3 Sparse-diffusion distance

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

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

This transitive sparsity constraint depends on a power-law parameter s ≥ 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⊗B where A and B are themselves both constrained to be sparse. More complicated are n₁ × n₂ 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 ≤ 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 |⋅|₁-norm of the matrix P, calculated as ∑_i∑_j|p_ij| is one possibility, as are terms which penalize t or α 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.

3.4 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(G₁, G₂) where

Theorem 5. Let G₁ and G₂ be graph box products as described above, and for a graph G let L(G) be its Laplacian. For fixed t = t_c, α = α_c, P⁽ⁱ⁾ as given above, for any λ ∈ [0, 1], we have

Proof. For graph products G_i, we have

Because A⊗I_|B| and I_|A|⊗B commute for any A and B,

We will make the following abbreviations:

Then,

By Lemma 1,

If we instead use $$ (E_2^{\left(1\right)}{P}^{\left(1\right)}\otimes P^{\left(2\right)}{E}_1^{\left(2\right)})$$ as the cross term in Eq (33), we have

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

This has the additional consequence that

Additionally, if

This reduces further to

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: Sq_n = Pa_n□Pa_n.

3.5 Spectral lower bound

In Theorem 7 we will derive and make use of an upper bound on the graph distance $$ \tilde{D}(G_1,G_2)$$. This upper bound is calculated by constraining the variable P to be not only orthogonal, but also $$ P=U_{2}M{U}_1^T$$ where M is the solution (i.e. “matching”, in the terminology of that section) to a Linear Assignment problem with costs given by a function of the eigenvalues of L(G₁) and L(G₂). In this section we derive a similar lower bound on the distance.

Let G₁ and G₂ be undirected graphs with Laplacians L₁ = L(G₁) and L₂ = L(G₂), and let α > 0 be constant. By Eq (52), we have

The following upper bound on $$ \tilde{D}$$ is achieved by constraining P to be not only orthogonal, but related to a constrained matching problem between the two lists of eigenvalues:

If we discard the constraints in Eq (45), we are left with the constraint that every column of p must have unit norm.

Construct the “spectral lower bound matching” matrix P^(SLB) as follows:

For any α, this matrix is the solution to a matching problem (less constrained than the original optimization over all P) where each $$ {\lambda }_j^{\left(1\right)}$$ is assigned to the closest $$ {\lambda }_i^{\left(2\right)}$$, allowing collisions. It clearly satisfies the constraints in Eq (44), but may violate those in Eq (45). Thus, we have

Various algorithms exist to rapidly find the member of a set of points which is closest to some reference point (for example, KD-Trees [29]). For any α, the spectral lower bound can be calculated by an outer loop over alpha and an inner loop which applies one of these methods. We do not consider joint optimization of the lower bound over P and α in this work.

3.6 Regularized distance

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

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

Theorem 6. D_reg satisfies the triangle inequality.

Proof. For graphs G₁, G₂, G₃ and Laplacians L₁, L₂, L₃, for any fixed t ≥ 0, we have:

Suppose that

Then,

Since this is true for any fixed t, let

Then

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

The term “regularized” here refers to the fact that the additional terms included in D_reg and $$ {\tilde{D}}_{\text{reg}}$$ penalize α 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.

3.7 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 the time scaling factor between coarse and fine scales. However, the case of an optimal (not power law) time conversion factor α 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.

3.8 Summary of distance metric versions

Table 1 summarizes the variants of our distance metric.

4.1 Algorithm development

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

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 α and t (where each innermost function evaluation consists of a full optimization over P at some t, α). We first briefly describe the algorithm used to find the optimal P and α 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 7). 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}\left(P\right)$$ (for instance, by replacing the optimization over all P with a less constrained matching problem—see Section 3.5).

Fig 4

Two plots demonstrating characteristics of distance calculation between a (7× 7) grid and an (8 × 8) grid.

(a): Plot illustrating the discontinuity and multimodality of the linear version of distance. Each gray curve represents a function $$ f_{P_c}\left({\alpha }_c\right)={\tilde{D}}^2({\text{Sq}}_7,{\text{Sq}}_8|{\alpha }_c,P_c)$$. The thicker curve is the lower convex hull of the thinner curves as a function of α, that is: $$ f\left({\alpha }_c\right)={inf}_{P\left|\mathcal{C}\right(P)}{\tilde{D}}^2({\text{Sq}}_7,{\text{Sq}}_8|{\alpha }_c)$$. We see that f(α) 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. (b): As in (a), but with D²(Sq₇, Sq₈|t = .318) plotted instead of $$ {\tilde{D}}^2$$. This t value is the location of the maximum in Fig 3.

4.2 Algorithm correctness proof

Theorem 7. For any two graphs G₁ and G₂ with Laplacians L(G₁) and L(G₂), for fixed α, the optimization over P given in the innermost loop of Eq 9 is upper bounded by a Linear Assignment Problem as defined in Eq (5). This LAP is given by taking R to be the eigenvalues $$ {\lambda }_j^{\left(1\right)}$$ of L(G₁) and S to be the eigenvalues $$ {\lambda }_i^{\left(2\right)}$$ of L(G₂), with the cost of a pair (equivalently, one entry of the cost matrix C) given by

Proof. L(G₁) and L(G₂) are both real symmetric matrices, so they may be diagonalized as $$ L\left(G_i\right)=U_i{\Lambda }_i{U}_i^T$$, where the U_i are rotation matrices, and the Λ_i are diagonal matrices with the eigenvalues $$ {\lambda }_1^{\left(i\right)},{\lambda }_2^{\left(i\right)}\dots {\lambda }_{n_i}^{\left(i\right)}$$ along the diagonal. Because the Frobenius norm is invariant under rotation, we have:

Because the U_i are orthogonal, the transformation $$ \tilde{P}=U_2^{T}P{U}_1$$ preserves orthogonality, so

For any given α,

Furthermore, this optimization is exactly a linear assignment problem of eigenvalues of L(G₁) to L(G₂), with the cost of a pair $$ ({\lambda }_j^{\left(1\right)},{\lambda }_i^{\left(2\right)})$$ given by

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

If we instead loosen the constraints on P, we can calculate a lower bound on the distance. See Appendix 3.5 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 Eq (49). For given G₁, G₂, and some matching M, let

Lemma 8. For two unique matchings M₁ and M₂ (for the same G₁, G₂) the equation $$ f_{M_1}\left(\alpha \right)-f_{M_2}\left(\alpha \right)=0$$ has at most one real positive solution in α. This follows from the fact that when P and t are fixed, the objective function is a rational function in α (see Eq (51)), with a quadratic numerator and an asymptote at α = 0.

Proof. By Eq (51), we have

Abbreviating the sums, we have

Since A₁, A₂, C₁, C₂ 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 9. If two matchings M₁ and M₃ which yield optimal upper bounds for the linear distance $$ {\tilde{D}}^2$$ (at α₁ ≤ α and α₃ ≥ α respectively) agree on a set of assignments, then the optimal M at α must also agree with that set of assignments.

Proof. We need the following lemmas:

Lemma 10. If an optimal matching assigns i to m(i) (so that eigenvalue $$ {\lambda }_i^{\left(1\right)}$$ of G₁ is paired with $$ {\lambda }_{f\left(i\right)}^{\left(2\right)}$$ of G₂ in the sum of costs Eq (49)), then the sequence m(1), m(2), …m(n₁) 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₁ < i₂ but m(i₂)<m(i₁)) then the new matching obtained by uncrossing those two pairs (performing a 2-opt step as defined in [33]) has strictly lesser objective function value. Hence an optimal matching can’t contain any such crossings.

Lemma 11. For all positive real α* ≥ ϵ > 0, let M₁ be an optimal matching at α* − ϵ and M₂ be optimal at α* + ϵ. For 1 ≤ i ≤ n₁, let s₁(i) and s₂(i) be the indices of λ⁽²⁾ paired with i in M₁ and M₂, respectively. Then for all i, s₁(i)≤s₂(i).

Proof. Define a “run” for s₁, s₂ as a sequence of consecutive indices l, l + 1, …l + k in [1, n₁] such that for any l, l + 1: min(s₁(l + 1), s₂(l + 1)) < max(s₁(l), s₂(l)). The following must be true about a “run”:

Within a run, either s₁(l)<s₂(l) or s₁(l)>s₂(l) for all l. Otherwise, we have one or more crossings (as in Lemma 10): for some l we have s₁(l)>s₁(l + 1) or s₂(l)>s₂(l + 1). Any crossing may be uncrossed for a strictly lower objective function value—violating optimality of M₁ or M₂.

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₁(i)<s₂(i) for all i. Let S = {l, l + 1, …l + k} be a run. By optimality of M₁, M₂ at α* − ϵ and α* + ϵ respectively, we have:

Respectively, these simplify to

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

Summing and reducing gives us

However, since the $$ {\lambda }_j^{\left(2\right)}$$ are monotonic nondecreasing, this means we cannot also have s₁(i)>s₂(i) for all i ∈ S, since that would imply

Therefore, in a run of arbitrary length, all indices must move ‘forward’ (meaning that s₁(i)<s₂(i) for all i in the run), and so (since any pair of matchings optimal at such α define a set of runs) we must have s₁(i)≤s₂(i). This completes the proof of the lemma.

Thus, for three matchings M₁, M₂, M₃ which optimal at a sequence of α₁ ≤ α₂ ≤ α₃, we must have s₁(i)≤s₂(i)≤s₃(i) for all i. In particular, if s₁(i) = s₃(i), we must also have s₁(i) = s₂(i) = s₃(i).

Theorem 12. If two matchings M₁ and M₃ yield optimal upper bounds for the linear distance $$ {\tilde{D}}^2$$ at α₁ and α₃ respectively, and $$ f_{M_1}\left({\alpha }_2\right)=f_{M_2}\left({\alpha }_2\right)$$ for some α₂ s.t. α₁ ≤ α₂ ≤ α₃, then either (1) M₁ and M₃ are optimal over the entire interval [α₁, α₃] or (1) some other matching M₂ improves over M₁ and M₃ at α₂.

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

Together, the preceeding properties verify that our algorithm will indeed find the joint optimum over all α 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) [34]. Numerical arrays were stored using the numpy package [35]. Our inner LAP solver was the package lapsolver [36]. Univariate optimization over t and α was performed with the ‘bounded’ method of the scipy.optimize package [37], with bounds set at [0, 10.0] for each variable and a input tolerance of 10⁻¹². Laplacians were computed with the laplacian method from the package networkX [38], 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 α. We leave mitigating this inefficiency for future work.

Code for computing diffusion distance, both with our algorithm and with naive univariate optimiztion, may be found in the S1 Data associated with this paper, as well as a maintained GitHub repository [39].

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 (Pa_n): 1D grid graphs of length n, with aperiodic boundary conditions.

Cycle Graphs (Cy_n): 1D grid graphs of length n, with periodic boundary conditions.

Square Grid Graphs (Sq_n): 2D grid graphs of dimensions n, with aperiodic boundary conditions. Sq_n = Pa_n□Pa_n

“Multi-Barbell” Graphs (Ba_n): Constructed as Cy_n□K_n, where K_n is the complete graph on n vertices.

These familes are all illustrated in Fig 5.

Fig 5

Graph lineages used in multiple numerical experiments in the main text.

Additionally, some examples distances between elements of these graph lineages are illustrated in Fig 6. In these tables we see that in general intra-lineage distances are small, and inter-lineage distances are large.

Fig 6

Distances D²(G, H) calculated for several pairs of graphs.

The top plot shows distances where G and H are both chosen from {Grid_13×13, P₁₆₉, C₁₆₉, Ba₁₃}. At bottom, distances are calculated from G chosen in {Grid_12×12, P₁₄₄, C₁₄₄, Ba₁₂} to H chosen in {Grid_13×13, P₁₆₉, C₁₆₉, Ba₁₃}. As expected, diagonal entries are smallest.

5.2 Numerical optimization methods

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

Nelder-Mead in Mathematica For very small graph pairs (n₁ × n₂ ≤ 100) we are able to find optimal P, α, t using constrained optimization in Mathematica 11.3 [40] using NMinimize, which uses Nelder-Mead as its backend by default. The size limitation made this approach unusable for any real experiments.

Orthogonally Constrained Opt. We also tried a variety of codes specialized for numeric optimization subject to orthogonality constraints. These included (1) the python package PyManopt [41], a code designed for manifold-constrained optimization; (2) gradient descent in Tensorflow using the penalty function g(P) = c||P^T P − I||_F (with c ≪ 1 a small positive constant weight) to maintain orthogonality, as well as (3) an implementation of the Cayley reparametrization method from [42] (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.

Black-Box Optimization Over α. We compare in more detail two methods of joint optimization over α and P when P is constrained to be a subpermutation matrix in the diagonal basis for L(G₁) and L(G₂). Specifically, we compare our approach given in Algorithm 1 to univariate optimization over α, where each function evaluation consists of full optimization over P. Fig 7 shows the results of this experiment. We randomly sample pairs of graphs as follows:

n₁ is drawn uniformly from [5, 120].

n₂ is drawn uniformly from [n₁, n₁ + 60].

G₁ and G₂ 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 }.

Fig 7

Comparison of runtimes for our algorithm and bounded golden section search over the same interval [10⁻⁶, 10].

Runtimes were measured by a weighted count of evaluations of the Linear Assignment Problem solver, with an n × n linear assignment problem counted as n³ units of cost. Because our algorithm recovers the entire lower convex hull of the objective function as a function of α, 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³. The most dramatic speedup occurs in the regime where n₁ ≪ n₂. Graphs were generated by drawing n₁ uniformly from [5, 120], drawing n₂ uniformly from [n₁, n₁ + 60], and then adding edges according to a Bernoulli distribution with p in {.125, .25, .375, .5, .625, .75, .875 } (60 trials each).

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α*, 1.618α*] for each local minimum α* 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³. 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₂ × n₂ linear assignment problem for each value of α 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.

5.3 Experiments

5.3.1 Triangle inequality violation of D (exponential distance) and $$ \tilde{D}$$ (linear distance)

As stated in Section 2.3, 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₁, G₂, G₃ of random graphs of sizes n_i for n₁ ∈ [5, 30], n₂ ∈ [n₁, n₁ + 30], and n₃ ∈ [n₂, n₂ + 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)∈{(1, 3), (1, 2), (2, 3)}). The results may be seen in Fig 8. In this figure we plot a histogram of the “discrepancy score”

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 × 10⁴ 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 Fig 8. Additionally, the triangle inequality is satisfied in 28184 (95.2%) of cases.

Fig 8

Histograms of triangle inequality violation.

These plots show the distribution of Disc(G₁, G₂, G₃), 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-ρ-pseudometric for some value of ρ ≈ 1.8 (linear version) or ρ ≈ 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: Disc(G₂, G₁, G₃) and Disc(G₃, G₂, G₁). Each plot has a line at x = 1, the maximum discrepancy score for which the underlying distances satisfy the triangle inequality.

We see similar but even stronger results when we run the same experiment with D² instead of $$ {\tilde{D}}^2$$; these may also be seen in Fig 8. 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 Disc(G₂, G₁, G₃) and Disc(G₃, G₂, G₁). All of these results are plotted in Fig 8.

5.3.2 Intra- and inter-lineage distances

We compute pairwise distances for sequences of graphs in the graph lineages displayed in Fig 5. For each pair of graph families (Square Grids, Paths, Cycles, and Multi-Barbells), we compute the distance from the ith 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.

Table 2

For two lineages G₁, G₂… (listed at left) and H_!, H₂, … (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.

Mean distances between graphs in several lineages.

	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

We note here that the idea of computing intra- and inter- lineage distances is similar to recent work [43] computing distances between graph ensembles: certain classes of similarly-generated random graphs. Graph diffusion distance has been previously shown (in [43]) to capture key structural information about graphs; for example, GDD is known to be sensitive to certain critical transitions in ensembles of random graphs as the random parameters are varied. This is also true for our time dilated version of GDD. More formally: let G_p and $$ G_p^{{}^{\prime }}$$ represent random graphs on n vertices, drawn from the Erdős-Renyi distribution with edge probability p. Then $$ D(G_p,G_p^{{}^{\prime }})$$ has a local maximum at $$ p=\frac{1}{n}$$, representing the transition between disconnected and connected graphs. This is true for our distance as well as the original version due to Hammond.

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 [7], whose definition of a graph limit (called a graphon) is as follows: Recall the definition of graph cut-distance D_cut(G, H) from Eq 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 graphs, under the equivalence relation that two sequences G₁, G₂, … and H₁, H₂, … are equivalent if D_cut(G_i, H_i) approaches 0 as n → ∞. Here we are calling a sequence of graphs “Cauchy” if for any ϵ > 0 there is some N such that for all n, m ≥ N, D_cut(G_n, G_m)<ϵ.

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 [1] 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. Furthermore, previous work [44] on sparse graph limits has shown that in the framework of Lovász all sequences of sparse graphs converge (in the infinite-size limit) to the zero graphon. Graph convergence results specific to sparse graphs include the Benjamini-Schramm framework [45], in which graph sequences are compared using the distributional limits of subgraph frequencies. These two graph comparison methods both have the characteristic that the “limit object” of a sequence of graphs is rigorously defined. In this section we attempt to show empirically that such a limit object of graph sequences under GDD may exist, and therefore merit further investigation.

We examine several sequences of graphs of increasing size for the required Cauchy behavior (in terms of our distance measure) to justify this variant definition of a “graph limit”. For each of the graph sequences defined in Section 5.1, we examine the distance between successive members of the sequence, plotting D²(G_n, H_n+1) for each choice of G and H. These sequences of distances are plotted in Fig 9.

Fig 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 Grid_n×n and Grid_(n+1)×(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.

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

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 → ∞.

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

We do not show that all similarly-constructed graph sequences display this Cauchy-like behavior. We hope to address this deeper question, as well as a more formal exploration of the limit object, with one or more modified versions of the objective function (see Section 3.6).

5.3.4 Limit of path graph distances

In this section, we demonstrate analytically that the sequence of path graphs of increasing size is Cauchy in the sense described by the previous section. In the following theorem (Theorem 14), we assume that the optimal value of t approaches some value $$ \tilde{t}$$ as n → ∞. We have not proven this to be the case, but have observed this behavior for both square grids and path graphs (see Fig 10 for an example of this behavior). Lemmas 13 and 14 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

Fig 10

Limiting behavior of D and two parameters as path graph size approaches infinity.

All distances were calculated between Path_n and Path_n+1. We plot the value of the objective function, as well as the optimal values of α and t, as n → ∞. Optimal α rapidly approach 1 and the optimal distance tends to 0. Additionally, the optimal t value approaches a constant (t ≈.316345), providing experimental validation of the assumption we make in proving Theorem 14.

Lemma 13. For any finite k, t, we have

Proof. Clearly for finite k, t

Then,

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

Since both of the limits

and

exist (and are 0),

and therefore

Theorem 14. If lim_n→∞ arg sup_t D² (Pa_n, Pa_n+1|t) exists, then:

Proof. Assume that $$ {lim}_{n\to \infty }arg{sup}_t{D}^2({\text{Pa}}_n,{\text{Pa}}_{n+1}|t)=\tilde{t}$$. Then, we must have

Hence, it remains only to prove that

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

Here, U_n and U_n+1 are the matrices which diagonalize L(Pa_n) and L(Pa_n+1) respectively (note also that a diagonalizer of a matrix L also diagonalizes e^L). If at each n we select S to be the subpermutation $$ S=\left[\begin{array}{c}I\hfill \\ 0\hfill \end{array}\right]$$, then (using the same argument as in Theorem 7) the objective function simplifies to:

By Lemma 13, for any finite k, t, we have

So for any ϵ > 0, ∃N such that when n ≥ N, for any c, k,

But then

as required. Thus, the Cauchy condition is satisfied for the lineage of path graphs Pa_n

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 Eq (41) from Theorem 5, we have (for any t_c, α_c):

As we can see in Fig 11, the right side of this inequality seems to be tending to a nonzero value as n → ∞, whereas the actual distance (calculated by our optimization procedure) appears to be tending to zero.

Fig 11

Comparison of the distance D(Sq_n, Sq_n+1) as a function of n, to the upper bound calculated as the optimum of distance between Pa_n and Pa_n+1.

We see that the upper found converges to some constant D ≈ 0.01782, whereas the actual distance appears to be converging to 0 as n → ∞.

5.3.5 Shape analysis for discretized meshes

In this section we demonstrate that graph diffusion distance captures structural properties of 3D point clouds. Ten 3D meshes (see Fig 12 for an illustration of the meshes used) were chosen to represent an array of objects with varying structural and topological properties. Not all of the mesh files chosen are simple manifolds: for example, the “y-tube” is an open-ended cylinder with a fin around its equator. Each mesh was used to produce multiple graphs, via the following procedure:

Subsampling the mesh to 1000 points;

Performing a clustering step on the new point cloud to identify 256 cluster centers;

Connecting each cluster center to its 16 nearest neighbors in the set of cluster centers.

Fig 12

3D meshes used in the shape analysis experiment.

Each mesh was used to produce several sampled discretizations, which were then compared using GDD.

Since each pass of this procedure (with different random seeds) varied in Step 1, each pass produced a different graph. We generated 20 graphs for each mesh, and compared the graphs using GDD.

The results of this experiment can be seen in Fig 13. This Figure shows the three first principal components of the distance matrix of GDD on the dataset of graphs produced as described above. Each point represents one graph in the dataset, and is colored according to the mesh which was used to generate it. Most notably, all the clusters are tight and do not overlap. Close clusters represent structurally similar objects: for example, the cluster of graphs from the tube mesh is very close to the cluster derived from the tube with an equatorial fin. This synthetic dataset example demonstrates that graph diffusion distance is able to compare structural information about point clouds and meshes.

Fig 13

Embedding of pairwise distances between mesh discretizations.

We see that GDD clusters each category of mesh tightly, and furthermore that clusters are nearby when they are structurally similar meshes, and distant otherwise. Axes represent the three principal components of the distance matrix and are thus unitless.

6.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.

6.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.

6.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 [46] 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. In [47] we make use of the optimization procedure described in this manuscript to find optimal (in the sense of Graph Diffusion Distance) coarsenings of a graph representing a protein nanotube. Since our procedure yields a P matrix which maps between the coarse and fine graphs, we are then able to define a machine learning model which learns to reproduce the protein’s energetic dynamics at multiple spatial scales. As a result of this multiscale construction, this multiscale model is more computationally efficient by an order of magnitude.

6.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. Section 5.3.5 showcases an example of doing this with a small synthetic database of 3D point clouds derived from mesh files of several well-known objects. This is similar to the approach of [4] 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.