From Worden

Jump to: navigation, search

Manuscript from the Laplacian Clustering project.

The material is here:

Laplacianoid - definition of the Laplacianoid and analysis of the Laplacianoid ODE
Laplacian Experiments/Laplacianoid ODEs - numerical exploration of the Laplacianoid ODE
Laplacian Experiments/Stochastic - the fully stochastic Laplacianoid process
Laplacian Paper/Figures - specific code for the figures for this paper

Following some discussion on Google+ about this manuscript, I am considering alternative names.

Mega-Laplacian? Meta-Laplacian? Double Laplacian? Diplacian? Dilaplacian?
- (I have looked into defining an operation on the graph Laplacian matrix rather than on the adjacency matrix, which would make it more of a meta-Laplacian, but it didn't seem to work so well.)
- Laplacienne?
- Maybe I'll have a naming contest after the draft is more written...

Some things that probably won't make it into the paper:

Analysis for directed graphs (I've done a fair bit of this, but there are fewer results)
Laplacian deformation of continuous metric spaces (I haven't done this but would like to)
Is there a Laplacian-derived measure of clustering

... after it's written, consider making a presentation of it as sequential art?

[log, pdf]Reflexive graph operations for clustering

Lee Worden, Junling Ma, Jonathan Dushoff

Abstract.

We use a graph’s Laplacian operator in a symmetrical way to transform the graph’s own structure, analogous to the use of the graph Laplacian operator in diffusion of a quantity along graph edges, and the use of the 2-D discrete Laplacian in image blurring. We look at the ability of this “matrix Laplacian” operation to add local clustering to a graph in three ways: by inducing continuous changes in the adjacency matrix, and interpreting it as a weighted graph; by making continuous changes in the adjacency matrix and interpreting it as a matrix of edge-presence probabilities; and by inducing discrete changes in the adjacency matrix using a stochastic process.

1 Introduction

In many disciplines there is a need to generate random graphs, in order to study how outcomes depend on graph structure, from dynamics of ecological food webs to spread of behaviors on social networks, the impact of contact network structure on disease transmission, behavior of networks of metabolic pathways, or search algorithms on the worldwide web. Exactly what aspects of a graph’s structure are salient can be a vexing question, and while degree distribution is the most often considered measure, clustering coefficient, degree assortativity, mean path lengths, and other graph characteristics can have a significant impact on outcomes as well.

High clustering coefficients are common to many measured real-world networks from various disciplines [], and highly clustered networks can have significantly different dynamical characteristics from otherwise similar graphs []. There are two popular methods for generating graphs with desired degree distributions and clustering coefficients, edge swapping techniques [] and exponential random graph models [], both of which are effective but computationally demanding and mathematically complex. In this paper we introduce a novel method of generating graphs with a desired clustering coefficient, which has the potential to be simpler though it may turn out to be less flexible. The process we introduce may turn out to have theoretical interest as well, because of its close relationship to the graph Laplacian operator, and the way in which it opens up the general subject of applying linear operations directly to graph structures.

The Laplacian operator originates in the study of the heat equation, in which it describes the temperature difference between any point and its surroundings, which drives the flow of heat from warmer to cooler locations. It is useful in describing diffusion of any number of quantities, from liquids and gases to crowds of random walkers, across any of a variety of spaces, Euclidean and nonlinear, continuous and discrete. The discrete Laplacian, on a lattice, approximates the continuous Laplacian on $\mathbb{R}^{2}$ , and is used in numerical approximations to PDEs and image processing. In image processing, Laplacian operations are used to add blurring to pixelated images, and in feature detection and image sharpening. The more general graph Laplacian, on arbitrary graph structures rather than lattices, can be used to model heat, random walks, and other diffusion processes on graphs.

In this paper we investigate the possibility of using a graph’s Laplacian reflexively on the structure of the graph itself, as a potential way of adding clustering to an existing graph. By analogy to blurring of an image, or diffusion of heat or random walkers across a graph, we use the Laplacian to make the graph’s edges themselves diffuse to neighboring locations, adding triangles and higher-order polygons where they were previously not. Once the properties of this operation are well understood, it may be useful for generating graphs with a desired combination of degrees and clustering. It may also have other unforeseen applications.

1.1 Laplacian diffusion of graph edges

Consider a graph $G$ with adjacency matrix $\mathbf{e}$ , so that $e_{{ij}}$ is 1 if there is an arc from vertex $i$ to vertex $j$ , and 0 otherwise. Let $\phi$ be a vector recording the temperature, or other diffusing quantity, at each vertex of $G$ . Then the values of $\phi$ flow across arcs of the graph, from higher to lower concentrations:

\frac{d}{dt}\phi _{i}=\sum _{j}e_{{ij}}(\phi _{j}-\phi _{i}).

This is equivalent to writing

\frac{d}{dt}\phi=\mathbf{L}(\mathbf{e})\phi

where $\mathbf{L}(\mathbf{e})$ is the Laplacian operator associated with $G$ . The Laplacian is a linear operator, so its action can be written as multiplication by a matrix¹¹Note that this matrix is distinct from the graph Laplacian matrix associated with $G$ , which for traditional reasons is generally taken to be the negative of this matrix. Here we are concerned with the Laplacian operator and its matrix representation.:

\left[\mathbf{L}(\mathbf{e})\right]_{{ij}}=\begin{cases}e_{{ij}}&\text{if }i\not=j\\ -\sum _{k}e_{{ik}}&\text{if }i=j.\end{cases}

This operator implements diffusion of a quantity across the graph $G$ . For instance, if initially $\phi$ is concentrated at one vertex, with $\phi _{i}=1$ and $\phi _{j}=0$ at all other vertices, under the influence of this diffusion dynamics the concentration at vertices neighboring $i$ will become positive, and the concentration at $i$ will decrease, until ultimately it has spread to all vertices reachable from $i$ .

To apply a similar dynamic to the structure of $G$ itself, where initially there is an edge connecting vertex $i$ to vertex $j$ , we would like that connectedness to “spill over” to neighboring locations, so that neighbors of $i$ become connected to $j$ , and $i$ becomes connected to neighbors of $j$ . The $j$ th column of $\mathbf{e}$ records which edges are connected to vertex $j$ , so naively one might like to make this happen by applying the Laplacian directly to the matrix $\mathbf{e}$ , which is equivalent to operating on each of its columns in just the same way as to $\phi$ :

\frac{d}{dt}\mathbf{e}=\mathbf{L}(\mathbf{e})\mathbf{e},

or

\frac{d}{dt}e_{{ij}}=\sum _{k}e_{{ik}}(e_{{kj}}-e_{{ij}}).

In this system, if initially there is a connection from vertex $i$ to $j$ , with $e_{{ij}}=1$ , and $e_{{kj}}=0$ for other $k$ , this “concentration” of connectedness will diffuse to neighbors of $i$ , making $e_{{kj}}>0$ for those $k$ and $e_{{ij}}<1$ .

However, this operation has the drawback that it makes no connections to neighbors of $j$ : it only operates on arrow tails, leaving arrow heads alone. To make it operate symmetrically, it’s necessary to operate on both the rows and columns of $\mathbf{e}$ , that is, to operate on the matrix from both the left and the right:

\frac{d}{dt}e_{{ij}}=\sum _{k}\left[e_{{ik}}(e_{{kj}}-e_{{ij}})+(e_{{ik}}-e_{{ij}})e_{{kj}}\right].

This double application of the graph Laplacian is a linear operator that acts on square matrices, rather than on the vectors on which the Laplacian acts. We may call this the graph’s matrix Laplacian operator, $\mathbf{L}_{m}$ : using this notation we can rewrite the above equation as

\frac{d}{dt}\mathbf{e}=\mathbf{L}_{m}(\mathbf{e})\mathbf{e}.

For any square matrix $\mathbf{x}$ the matrix Laplacian of $\mathbf{x}$ is a linear operator $\mathbf{L}_{m}(\mathbf{x})$ on the space of matrices of the same size as $\mathbf{x}$ :

\left[\mathbf{L}_{m}(\mathbf{x})\mathbf{y}\right]_{{ij}}=\sum _{k}\left[x_{{ik}}(y_{{kj}}-y_{{ij}})+(y_{{ik}}-y_{{ij}})x_{{kj}}\right],

or

\mathbf{L}_{m}(\mathbf{x})\mathbf{y}=\mathbf{L_{o}}(\mathbf{x})\mathbf{y}+\mathbf{y}\mathbf{L_{i}}(\mathbf{x}),

where $\mathbf{L_{o}}(\mathbf{x})$ is the “out-degree Laplacian,” with vertices’ out degrees on the diagonal (equal to the $\mathbf{L}(\mathbf{x})$ defined above), and $\mathbf{L_{i}}(\mathbf{x})$ is the corresponding in-degree Laplacian:

\left[\mathbf{L_{i}}(\mathbf{e})\right]_{{ij}}=\begin{cases}e_{{ij}}&\text{if }i\not=j\\ -\sum _{k}e_{{kj}}&\text{if }i=j.\end{cases}

This operation “diffuses” any nonzero edge $e_{{ij}}$ to neighbors of both $i$ and $j$ . Note that the notation $\mathbf{L}_{m}(\mathbf{x})\mathbf{y}$ refers only to application of a linear operator $\mathbf{L}_{m}(\mathbf{x})$ , not to matrix multiplication.

The matrix Laplacian is a linear operator that preserves nonnegativity and symmetry of its matrix arguments, that is, when both $\mathbf{x}$ and $\mathbf{y}$ represent undirected graphs, the result $\mathbf{L}_{m}(\mathbf{x})\mathbf{y}$ does as well. The above matrix Laplacian differential equation transforms unweighted graphs — those with adjacency matrices containing only 0 and 1 — into weighted graphs whose adjacency matrix entries are real numbers in $[0,1]$ .

We study what the reflexive matrix Laplacian operation $\mathbf{L}_{m}(\mathbf{e})\mathbf{e}$ does to a graph’s structure over time, in three ways. First we use the above ordinary differential equation as is, which produces a continuously changing weighted graph. We consider the properties of this weighted graph itself, and then look at two ways of producing unweighted graphs from it. Then we turn to a stochastic process on unweighted graphs that is constructed from the same reflexive Laplacian operation on the graph’s adjacency matrix, and give the same consideration to the effects of this process on the graph’s structure.

In what follows, all the graph measures we use are defined in terms of the graph’s adjacency matrix, for mathematical clarity. Given adjacency matrix $\mathbf{e}$ for graph $G$ , we define the density of $G$ as $\sum _{{i,j}}e_{{ij}}/{N^{2}}$ , where $N$ is the number of vertices in the graph and the indices of summation range from 1 to $N$ . This is equal to the standard measure of density, with loops — edges from a vertex to itself — included. Transitivity is defined as

\frac{\sum _{{i,j,k}}e_{{ij}}e_{{jk}}e_{{ik}}}{\sum _{{i,j,k}}e_{{ij}}e_{{jk}}},

that is, the proportion of all the paths of length two whose transitive edge is present in the graph. This is equivalent to the conventional clustering coefficient for undirected graphs,

\frac{3\cdot\text{\# of triangles}}{\text{\# of connected triples of vertices}},

except that it explicitly counts triangles and 2-paths that include loops, and can be used on directed graphs. The matrix Laplacian operations tend to create loops, so this formulation seems appropriate, in addition to being simple to compute in terms of an adjacency matrix.

2 Methods

2.1 Deterministic diffusion of weighted graph edges

To begin with, we will describe the behavior of the matrix ordinary differential equation

\frac{d}{dt}\mathbf{e}(t)=\mathbf{L}_{m}(\mathbf{e}(t))\mathbf{e}(t).

(1)

As we will see, this ODE is well-behaved in that when edge weights are contained in the interval $[0,1]$ , they remain in that interval, and in that undirected graphs remain undirected, that is, the space of symmetric matrices is an invariant region of the dynamics. The behavior of this system is straightforward and instructive.

2.2 Generation of unweighted graphs by Bernoulli sampling

The above differential equation defines a continuous-time process in which the weights of the graph’s edges, which are initially all either 1 or 0, become smoothly changing real values between 0 and 1. It is often desirable to generate graphs with unweighted edges, and one way to do this is to construe the weighted edges generated by the ODE as probabilities, and create unweighted random graphs by defining each edge to be a Bernoulli random variable (independent of the other edges) with the probability given by its counterpart in the weighted graph. As we will see, this process can be used to generate graphs with nontrivial clustering, but it does not seem to be as powerful as the fully stochastic process we describe next.

2.3 Stochastic diffusion of unweighted graph edges

After examining the above differential equation, which produces smoothly varying graph edge weights, we turn to a closely related model, in which edge weights are only 0 and 1, and flip on and off stochastically, at rates given by the graph’s matrix Laplacian applied reflexively.

We define this model by converting the above ODE’s rate of change

\frac{d}{dt}e_{{ij}}=\left[\mathbf{L}_{m}(\mathbf{e})\mathbf{e}\right]_{{ij}}=\sum _{k}\left[e_{{ik}}(e_{{kj}}-e_{{ij}})+(e_{{ik}}-e_{{ij}})e_{{kj}}\right]

to an average rate of change. When all matrix entries are 0 or 1, this quantity can be seen to be nonnegative when $e_{{ij}}$ is 0, and nonpositive when $e_{{ij}}$ is 1. Therefore we can simply consider a stochastic model in which each adjacency matrix entry $e_{{ij}}$ flips its state from 0 to 1 or vice versa at rate $|\left[\mathbf{L}_{m}(\mathbf{e})\mathbf{e}\right]_{{ij}}|$ . We restrict this operation to undirected graphs by changing symmetric pairs of matrix entries as a unit.

We will show that the absorbing states of this model satisfy the same conditions as the equilibria of the ODE: all the edges among any given triangle of vertices must be equal — either all zero or all one — implying a union of disjoint cliques. Unlike the ODE, this system can and often does come to rest at a disconnected graph. Whether connected or not, of course, this clique structure maximizes transitivity.

3 Results

3.1 Deterministic diffusion of weighted graph edges

In the matrix differential equation (1), given an initial graph that connects all the graph’s vertices, the initial graph gradually transforms into one in which all possible edges are present with equal weight, and the total weight is equal to the total weight of the initial graph’s edges. First the transitive edges implied by the initial graph’s edges begin to appear, then the transitive edges generated by those, and so on, until the entire transitive closure of the graph is visible, and finally the edge weights all equalize. See below for a mathematical analysis of this process.

Figure 1: Snapshots of a continuously changing weighted graph following equation 1.

Figure 2: Edge weights vs. time in the same graph shown in figure 1

3.1.1 Linear algebra of $\mathbf{L}_{m}$

$\mathbf{L}_{m}(\mathbf{e})$ is a linear operator on the $n^{2}$ -dimensional space of $n\times n$ matrices, $\mbox{Mat}_{{n{\times}n}}$ , and has $n^{2}$ eigenvalues and the same number of linearly independent eigenmatrices (rather than eigenvectors). If $\mathbf{v}_{i}$ , $i=1,\ldots,n$ are the $n$ eigenvectors of $\mathbf{L}(\mathbf{e})$ with eigenvalues $\lambda _{i}$ , the $n^{2}$ eigenmatrices of $\mathbf{L}_{m}(\mathbf{e})$ are the matrices $\mathbf{v}_{i}\mathbf{v}^{{\mathrm{T}}}_{j}$ , with eigenvalues $\lambda _{{ij}}=\lambda _{i}+\lambda _{j}$ :

	$\displaystyle\mathbf{L}_{m}(\mathbf{e})\mathbf{v}_{i}\mathbf{v}_{j}^{\mathrm{T}}$	$\displaystyle=\mathbf{L}(\mathbf{e})\mathbf{v}_{i}\mathbf{v}_{j}^{\mathrm{T}}+\mathbf{v}_{i}\mathbf{v}_{j}^{\mathrm{T}}\mathbf{L}(\mathbf{e})$
		$\displaystyle=(\lambda _{i}+\lambda _{j})\mathbf{v}_{i}\mathbf{v}_{j}^{\mathrm{T}}.$

When $\mathbf{e}$ is symmetric, its matrix Laplacian operator $\mathbf{L}_{m}(\mathbf{e})$ maps symmetric matrices to symmetric matrices, so that it can also be considered as a linear operator on $\mbox{Sym}_{n}$ , the $n(n+1)/2$ -dimensional space of $n\times n$ symmetric matrices. The eigenmatrices of this restricted Laplacian operator are the $n(n+1)/2$ symmetrized combinations of the above matrices, $(\mathbf{v}_{i}\mathbf{v}_{j}^{\mathrm{T}}+\mathbf{v}_{j}\mathbf{v}_{i}^{\mathrm{T}})$ .

In this section we consider only symmetric $\mathbf{e}$ (equivalently, assume that the graph is undirected), so that $\mathbf{L}(\mathbf{e})=\mathbf{L_{o}}(\mathbf{e})=\mathbf{L_{i}}(\mathbf{e})$ is also symmetric.

When $\mathbf{e}$ is symmetric, $\mathbf{L}(\mathbf{e})$ has nonnegative eigenvalues, so the only zero eigenvalues of $\mathbf{L}_{m}(\mathbf{e})$ come from adding zero eigenvalues of $\mathbf{L}(\mathbf{e})$ . If the graph $G$ is connected, the only null vector of the Laplacian is the vector of all ones, $\mathbf{1}$ , so the only null matrix of $\mathbf{L}_{m}(\mathbf{e})$ is $\mathbf{1}\text{ }\mathbf{1}^{{\mathrm{T}}}$ .

If the graph is not connected, then it is a union of disjoint subgraphs $G_{1},\ldots,G_{k}$ with adjacency matrices $\mathbf{e}_{1},\ldots,\mathbf{e}_{k}$ , with $\mathbf{e}=\mathbf{e}_{1}+\cdots+\mathbf{e}_{k}$ . The nullspace of $\mathbf{L}(\mathbf{e})$ contains $k$ orthogonal eigenvectors $\mathbf{w}_{i}$ with $(\mathbf{w}_{i})_{j}=1$ if vertex $j$ is part of subgraph $G_{i}$ , $0$ otherwise. It follows that the nullspace of $\mathbf{L}_{m}(\mathbf{e})$ in $\mbox{Mat}_{{n\times n}}$ is $k^{2}$ -dimensional and spanned by the eigenmatrices $\mathbf{w}_{i}\mathbf{w}_{j}^{\mathrm{T}}$ constructed by these nullvectors. The nullspace of $\mathbf{L}_{m}(\mathbf{e})$ in $\mbox{Sym}_{n}$ is the restriction of the above nullspace to $\mbox{Sym}_{n}$ , the span of the $k(k+1)/2$ symmetric combinations $(\mathbf{w}_{i}\mathbf{w}_{j}^{\mathrm{T}}+\mathbf{w}_{j}\mathbf{w}_{i}^{\mathrm{T}})$ of nullvectors.

It is sometimes helpful to write the Laplacian operator of a graph $\mathbf{e}$ in matrix form,

\mathbf{L}(\mathbf{e})_{{ij}}=e_{{ij}}-\sum _{k}e_{{ik}}\delta _{{ij}}

where $\delta _{{ij}}$ is $1$ when $i=j$ and $0$ otherwise. The corresponding expression for the matrix Laplacian is

\mathbf{L}_{m}(\mathbf{e})_{{ij,k\ell}}=e_{{ik}}\delta _{{\ell j}}+e_{{\ell j}}\delta _{{ik}}-\sum _{p}(e_{{ip}}+e_{{pj}})\delta _{{ik}}\delta _{{\ell j}}.

This can be used directly to express the matrix Laplacian operation:

	$\displaystyle(\mathbf{L}_{m}(\mathbf{e})\mathbf{x})_{{ij}}$	$\displaystyle=\sum _{{k\ell}}\mathbf{L}_{m}(\mathbf{e})_{{ij,k\ell}}x_{{k\ell}}$
		$\displaystyle=\sum _{k}e_{{ik}}x_{{kj}}+\sum _{{\ell}}e_{{\ell j}}x_{{i\ell}}-\sum _{p}(e_{{ip}}+e_{{pj}})x_{{ij}}$

In this way we can think of the matrix Laplacian as a sort of hypermatrix — an object with four indices that operates on matrices the way a matrix operates on vectors. This notation will be useful in working with the Jacobian of the matrix Laplacian dynamics, below.

3.1.2 Equilibria of the matrix Laplacian differential equation

The differential equation $\dot{\mathbf{e}}=\mathbf{L}_{m}(\mathbf{e})\mathbf{e}$ , i.e.,

	$\displaystyle\dot{e}_{{ij}}$	$\displaystyle=\sum _{k}\left[e_{{ik}}(e_{{kj}}-e_{{ij}})+(e_{{ik}}-e_{{ij}})e_{{kj}}\right]$
		$\displaystyle=\sum _{k}\left[2e_{{ik}}e_{{kj}}-(e_{{ik}}+e_{{kj}})e_{{ij}}\right],$

describes a square matrix whose entries change continuously, or equivalently, a weighted graph changing continuously in time on a fixed set of vertices. Equilibria of this system are matrices corresponding to particular graph structures. Simulation shows that given a starting graph that is connected, this system tends to relax to a matrix in which $e_{{ij}}=k$ for all $i,j$ , that is, to a graph whose edge weights are all equal. Analysis reveals that the full set of equilibria is a bit larger.

An equilibrium solution of the ODE (1) is a matrix that is in the nullspace of its own matrix Laplacian. This provides the general symmetric equilibrium solution of (1): it is a matrix $\mathbf{e}^{*}$ which on the one hand corresponds to an undirected graph $G$ partitioned into one or more disjoint subgraphs $G_{i}$ , and on the other hand can be expressed as a linear combination of the basis matrices of the nullspace of its own matrix Laplacian,

\mathbf{e}^{*}=\sum _{{s\leq t}}\alpha _{{st}}(\mathbf{w}_{s}\mathbf{w}_{t}^{\mathrm{T}}+\mathbf{w}_{t}\mathbf{w}_{s}^{\mathrm{T}})

where $\mathbf{w}_{s}$ and $\mathbf{w}_{t}$ range over the nullvectors of the Laplacian of $\mathbf{e}^{*}$ , as above. For $G$ to be partitioned, we must have $\alpha _{{st}}=0$ when $s\not=t$ . Thus an equilibrium graph is a linear combination of $m$ disjoint complete graphs,

\mathbf{e}^{*}=\sum _{{s=1}}^{m}k_{{s}}\mathbf{w}_{s}\mathbf{w}_{s}^{\mathrm{T}}.

The $i,j$ entry of this matrix is

	$\displaystyle\mathbf{e}^{*}_{{ij}}$	$\displaystyle=\sum _{{s}}k_{s}(\mathbf{w}_{s})_{i}(\mathbf{w}_{s})_{j}$
		$\displaystyle=\sum _{s}k_{s}I_{{G_{s}}}(i)I_{{G_{s}}}(j),$

where $I_{{G_{s}}}(i)$ is the indicator function for membership of vertex $i$ in subgraph $G_{s}$ . This gives us an $m$ -dimensional linear space of equilibrium solutions.

However, we can show that (1) conserves the sum of all edge weights, given that $\mathbf{e}$ is symmetric:

	$\displaystyle\frac{d}{dt}\sum _{{i,j}}e_{{ij}}$	$\displaystyle=\sum _{{i,j}}\left[\mathbf{L}_{m}(\mathbf{e}(t))\mathbf{e}(t)\right]_{{ij}}$
		$\displaystyle=\sum _{{i,j,k}}\left[2e_{{ik}}e_{{kj}}-e_{{ij}}(e_{{ik}}+e_{{kj}})\right]$
		$\displaystyle=\sum _{{i,j,k}}\left[2e_{{ik}}e_{{kj}}-e_{{ji}}e_{{ik}}-e_{{kj}}e_{{ji}}\right]$
		$\displaystyle=0.$

Therefore the total edge weight at equilibrium must be the same as in the initial graph. This reduces the space of attainable equilibria to $m-1$ dimensions.

Also, as mentioned above, when all the matrix entries (edge weights) $e_{{ij}}$ are in $[0,1]$ , it follows directly that $de_{{ij}}/dt\geq 0$ when $e_{{ij}}=0$ and $de_{{ij}}/dt\leq 0$ when $e_{{ij}}=1$ , which is sufficient to guarantee that entries remain within $[0,1]$ for all future times, including at equilibrium.

3.1.3 Stability of equilibria

To assess the stability of the ODE’s equilibria, we need to look at its Jacobian. Since

\dot{e}_{{ij}}=\sum _{k}\left[2e_{{ik}}e_{{kj}}-(e_{{ik}}+e_{{kj}})e_{{ij}}\right],

we have

\frac{\partial\dot{e}_{{ij}}}{\partial e_{{ij}}}=2e_{{ii}}+2e_{{jj}}-2e_{{ij}}-\sum _{k}(e_{{ik}}+e_{{kj}}),

and for $k\not=i$ and $\ell\not=j$ ,

	$\displaystyle\frac{\partial\dot{e}_{{ij}}}{\partial e_{{i\ell}}}$	$\displaystyle=2e_{{\ell j}}-e_{{ij}},$
	$\displaystyle\frac{\partial\dot{e}_{{ij}}}{\partial e_{{kj}}}$	$\displaystyle=2e_{{ik}}-e_{{ij}},$
	$\displaystyle\frac{\partial\dot{e}_{{ij}}}{\partial e_{{k\ell}}}$	$\displaystyle=0.$

We can encapsulate this by writing (for all $k,\ell$ )

	$\displaystyle\frac{\partial\dot{e}_{{ij}}}{\partial e_{{k\ell}}}$	$\displaystyle=(2e_{{ik}}-e_{{ij}})\delta _{{\ell j}}+(2e_{{\ell j}}-e_{{ij}})\delta _{{ik}}-\sum _{p}(e_{{ip}}+e_{{pj}})\delta _{{ik}}\delta _{{\ell j}}$
		$\displaystyle=\mathbf{L}_{m}(\mathbf{e})_{{ij,k\ell}}+(e_{{ik}}-e_{{ij}})\delta _{{\ell j}}+(e_{{\ell j}}-e_{{ij}})\delta _{{ik}}.$

As with a regular vector dynamical system, if $\mathbf{e}^{*}$ is an equilibrium, at a point near the equilibrium $\mathbf{e}=\mathbf{e}^{*}+\mathbf{x}$ the derivative is approximated by the linearization, $\dot{e}_{{ij}}=\dot{x}_{{ij}}\approx\sum _{{k,\ell}}\frac{\partial\dot{e}_{{ij}}}{\partial e_{{k\ell}}}x_{{k\ell}}$ , and we want to know if there is any such $\mathbf{x}$ that grows rather than shrinking. Of course, this really is a regular vector dynamical system, notwithstanding the fact that the state variables are written in a square rather than a column, so all the linear theory applies, including the fact that $\mathbf{e}^{*}$ is hyperbolically stable if and only if all the eigenvalues of the Jacobian are negative.

We analyze the eigenvalues by examining the quadratic form associated with the Jacobian, because if we can show that the quadratic form is never positive (that is, if it is a negative definite or negative semidefinite form) we can conclude that the eigenvalues of the Jacobian are all negative or non-positive. Because the ODE is restricted to symmetric matrices as long as the initial value is symmetric, we can restrict ourselves to the action of the Jacobian on symmetric perturbations $\mathbf{x}$ . The quadratic form is

	$\displaystyle Q(\mathbf{x})$	$\displaystyle=\sum _{{i,j,k,\ell}}x_{{ij}}\frac{\partial\dot{e}_{{ij}}}{\partial e_{{k\ell}}}x_{{k\ell}}$
		$\displaystyle=\sum _{{i,j,k,\ell}}x_{{ij}}[(2e_{{\ell j}}-e_{{ij}})\delta _{{ik}}+(2e_{{ik}}-e_{{ij}})\delta _{{\ell j}}-\sum _{p}(e_{{ip}}+e_{{pj}})\delta _{{ik}}\delta _{{\ell j}}]x_{{k\ell}}$
		$\displaystyle=\sum _{{i,j,\ell}}(2e_{{\ell j}}-e_{{ij}})x_{{ij}}x_{{i\ell}}+\sum _{{i,j,k}}(2e_{{ik}}-e_{{ij}})x_{{ij}}x_{{kj}}-\sum _{{i,j,p}}(e_{{ip}}+e_{{pj}})x_{{ij}}^{2}$
		$\displaystyle=\sum _{{c,b,a}}2e_{{ab}}x_{{cb}}x_{{ca}}-\sum _{{a,b,c}}e_{{ab}}x_{{ab}}x_{{ac}}+\sum _{{a,c,b}}2e_{{ab}}x_{{ac}}x_{{bc}}$
		$\displaystyle\quad\quad\mbox{}-\sum _{{a,b,c}}e_{{ab}}x_{{ab}}x_{{cb}}-\sum _{{a,c,b}}e_{{ab}}x_{{ac}}^{2}-\sum _{{c,b,a}}e_{{ab}}x_{{cb}}^{2}$
		$\displaystyle=\sum _{{a,b,c}}e_{{ab}}\left(2x_{{cb}}x_{{ca}}-x_{{ab}}x_{{ac}}+2x_{{ac}}x_{{bc}}-x_{{ab}}x_{{cb}}-x_{{ac}}^{2}-x_{{cb}}^{2}\right),$

and for symmetric $\mathbf{x}$ ,

	$\displaystyle Q(\mathbf{x})$	$\displaystyle=\sum _{{a,b,c}}e_{{ab}}\left(4x_{{ac}}x_{{bc}}-x_{{ab}}x_{{ac}}-x_{{ab}}x_{{bc}}-x_{{ac}}^{2}-x_{{bc}}^{2}\right)$
		$\displaystyle=\sum _{{a,b,c}}e_{{ab}}\left(-(x_{{ac}}-x_{{bc}})^{2}-x_{{ac}}(x_{{ab}}-x_{{bc}})-x_{{bc}}(x_{{ab}}-x_{{ac}})\right).$

Near the connected equilibrium,

\displaystyle Q(\mathbf{x})

\displaystyle=\sum _{{a,b,c}}k\left(-(x_{{ac}}-x_{{bc}})^{2}-x_{{ac}}(x_{{ab}}-x_{{bc}})-x_{{bc}}(x_{{ab}}-x_{{ac}})\right)

but since we are considering only symmetric $\mathbf{x}$ , we can write

	$\displaystyle\sum _{{a,b,c}}x_{{ac}}(x_{{ab}}-x_{{bc}})$	$\displaystyle=\sum _{{a,b,c}}x_{{ac}}x_{{ab}}-\sum _{{a,b,c}}x_{{ac}}x_{{bc}}$
		$\displaystyle=\sum _{{a,b,c}}x_{{ac}}x_{{ab}}-\sum _{{a,b,c}}x_{{ca}}x_{{cb}}$
		$\displaystyle=0;$

and the third term of the above expression cancels in the same way as the second, leaving us with

\displaystyle Q(\mathbf{x})

\displaystyle=-k\sum _{{a,b,c}}(x_{{ac}}-x_{{bc}})^{2},

which shows that $Q$ is negative semidefinite. $Q(\mathbf{x})$ is zero when $x_{{ac}}=x_{{bc}}$ for all $a,b,c$ ; for symmetric $\mathbf{x}$ this condition implies that all entries of $\mathbf{x}$ are equal: $x_{{ij}}=x_{{\ell j}}=x_{{j\ell}}=x_{{k\ell}}$ . These are the eigenmatrices of the Jacobian corresponding to a zero eigenvalue, and all its other eigenvalues are negative.

Consequently this $\mathbf{e}^{*}$ is a hyperbolically stable fixed point, except that it is neutral to the addition of a constant to all entries of the matrix. The ODE conserves the sum of all matrix entries, so in fact this fixed point is hyperbolically stable.

When $\mathbf{e}^{*}$ is made of disjoint cliques, the result is different:

\displaystyle Q(\mathbf{x})

\displaystyle=\sum _{s}\sum _{{a,b\in G_{s},c\in G}}k_{s}\left(-(x_{{ac}}-x_{{bc}})^{2}-x_{{ac}}(x_{{ab}}-x_{{bc}})-x_{{bc}}(x_{{ab}}-x_{{ac}})\right)

If $\mathbf{x}$ has an edge $x_{{ij}}>0$ from one clique to another, we can have a positive term at $(a,b,c)=(i,i,j)$ .

The linearized system is

	$\displaystyle\dot{x}_{{ij}}$	$\displaystyle=\sum _{{k,\ell}}\frac{\partial\dot{e}_{{ij}}}{\partial e_{{k\ell}}}x_{{k\ell}}$
		$\displaystyle=\sum _{{k,\ell}}[(2e_{{\ell j}}-e_{{ij}})\delta _{{ik}}+(2e_{{ik}}-e_{{ij}})\delta _{{\ell j}}-\sum _{p}(e_{{ip}}+e_{{pj}})\delta _{{ik}}\delta _{{\ell j}}]x_{{k\ell}}$
		$\displaystyle=\sum _{\ell}(2e_{{\ell j}}-e_{{ij}})x_{{i\ell}}+\sum _{k}(2e_{{ik}}-e_{{ij}})x_{{kj}}-\sum _{p}(e_{{ip}}+e_{{pj}})x_{{ij}},$

and for this $\mathbf{e}^{*}$ this has two cases: if $i$ and $j$ are in the same $G_{s}$ , i.e. if $e_{{ij}}>0$ ,

	$\displaystyle\dot{x}_{{ij}}$	$\displaystyle=\sum _{{\ell\in G_{s}}}k_{s}x_{{i\ell}}-\sum _{{\ell\not\in G_{s}}}k_{s}x_{{i\ell}}+\sum _{{k\in G_{s}}}k_{s}x_{{kj}}-\sum _{{k\not\in G_{s}}}k_{s}x_{{kj}}-\sum _{{p\in G_{s}}}2k_{s}x_{{ij}}$
		$\displaystyle=\sum _{{k\in G_{s}}}k_{s}(x_{{ik}}+x_{{kj}}-2x_{{ij}})-\sum _{{k\not\in G_{s}}}k_{s}(x_{{ik}}+x_{{kj}}),$

and if $i\in G_{s}$ and $j\in G_{t}$ ( $s\not=t$ ),

	$\displaystyle\dot{x}_{{ij}}$	$\displaystyle=\sum _{{\ell\in G_{t}}}2k_{t}x_{{i\ell}}+\sum _{{k\in G_{s}}}2k_{s}x_{{kj}}-\sum _{{p\in G_{s}}}k_{s}x_{{ij}}-\sum _{{p\in G_{t}}}k_{t}x_{{ij}}$
		$\displaystyle=\sum _{{k\in G_{s}}}k_{s}(2x_{{kj}}-x_{{ij}})+\sum _{{k\in G_{t}}}k_{t}(2x_{{ik}}-x_{{ij}})$
		$\displaystyle=2\sum _{{k\in G_{s}}}k_{s}x_{{kj}}+2\sum _{{k\in G_{t}}}k_{t}x_{{ik}}-(k_{s}\|G_{s}\|+k_{t}\|G_{t}\|)x_{{ij}}.$

We wish to construct an eigenmatrix with a positive eigenvalue, that is, a symmetric matrix $\mathbf{x}$ such that $\dot{\mathbf{x}}=\lambda\mathbf{x}$ with $\lambda>0$ . This will demonstrate that a union of multiple disjoint cliques is not a stable equilibrium.

3.1.4 A positive eigenmatrix

Given an equilibrium matrix $\mathbf{e}^{*}$ made of multiple disjoint cliques, let us choose two distinct cliques $G_{s}$ and $G_{t}$ , and let $x_{{ij}}=0$ except that $x_{{ij}}$ is $z$ for every edge connecting $G_{s}$ and $G_{t}$ (and vice versa), $u$ for every edge within $G_{s}$ and $v$ for every edge within $G_{t}$ .

Then for an edge from $G_{s}$ to $G_{t}$

\displaystyle\dot{x}_{{ij}}

\displaystyle=(|G_{s}|k_{s}+|G_{t}|k_{t})z,

and the value is the same for an edge from $G_{t}$ to $G_{s}$ , while for an edge within $G_{s}$

\displaystyle\dot{x}_{{ij}}

\displaystyle=-2|G_{t}|k_{s}z

and within $G_{t}$

\displaystyle\dot{x}_{{ij}}

\displaystyle=-2|G_{s}|k_{t}z.

For all other edges $\dot{x}_{{ij}}=0$ .

So for this to be an eigenmatrix we must have

	$\displaystyle(\|G_{s}\|k_{s}+\|G_{t}\|k_{t})z$	$\displaystyle=\lambda z$
	$\displaystyle-2\|G_{t}\|k_{s}z$	$\displaystyle=\lambda u$
	$\displaystyle-2\|G_{s}\|k_{t}z$	$\displaystyle=\lambda v.$

We can satisfy these equalities by setting

\lambda=|G_{s}|k_{s}+|G_{t}|k_{t},

and given this we set $u=-2|G_{t}|k_{s}z/\lambda$ and $v=-2|G_{s}|k_{t}z/\lambda$ , and the matrix is fully determined for any choice of $z$ .

Notice that the sum of all entries of $\mathbf{x}$ is

	$\displaystyle\sum _{{ij}}x_{{ij}}$	$\displaystyle=\|G_{s}\|^{2}u+2\|G_{s}\|\|G_{t}\|z+\|G_{t}\|^{2}v$
		$\displaystyle=(-2\|G_{s}\|^{2}\|G_{t}\|k_{s}/\lambda+2\|G_{s}\|\|G_{t}\|-2\|G_{s}\|\|G_{t}\|^{2}k_{t}/\lambda)z$
		$\displaystyle=-2\|G_{s}\|\|G_{t}\|z\frac{\|G_{s}\|k_{s}-(\|G_{s}\|k_{s}+\|G_{t}\|k_{t})+\|G_{t}\|k_{t}}{\|G_{s}\|k_{s}+\|G_{t}\|k_{t}}$
		$\displaystyle=0,$

confirming that the growth of $\mathbf{x}$ takes place within the invariant manifold of matrices with sum of entries conserved.

(Note that $\mathbf{x}$ is not itself a graph’s adjacency matrix, it is a perturbation to $\mathbf{e}^{*}$ , so it can have negative entries. It adds weight to edges that are zero, and subtracts from edges that are positive.)

3.2 Generation of unweighted graphs by Bernoulli sampling

Figures 3–7 present the results of taking the real-valued edge weights given by the matrix Laplacian differential equation as probabilities for randomly sampled edges. In the early stages of the ODE system’s evolution, this sampling indeed adds structure to the graph, producing configurations that are similar to the initial graph but have nonzero transitivity. Later, however, the transitivity declines. This is not mysterious when one considers that in the final state of the ODE, all edge weights are equal. This implies that the final state of the random sampling process is an Erdős-Rényi random graph, which is known to have vanishing transitivity when the graph’s size is large relative to its average degree.

Figure 3: Graphs generated by sampling edges using the probabilities given by the edge weights shown in figure 1.

Figure 4: Graph transitivity vs. time in the above randomly sampled graphs

Figure 5: Graph transitivity vs. time in an ensemble of

n

randomly sampled graphs made by applying equation 1 to random power-law graphs of 100 vertices.

Figure 6: Graph density vs. time in the same ensemble as figure 5.

Figure 7: Degree distribution vs. time in the ensemble of figure 5 (upper left:

t=0

; upper right:

t=0.01

; lower left:

t=0.02

; lower right:

t=0.09

).

3.3 Stochastic diffusion of unweighted graph edges

Figure 8: Snapshots of a stochastically changing unweighted graph following equation 2 (left:

t=0

; center:

t=0.1

; right:

t\to\infty

).

Figure 9: Graph transitivity vs. time in the system of figure 8.

Figure 10: Graph transitivity vs. time in an ensemble of

n

randomly sampled graphs made by applying equation 2 to random power-law graphs of 100 vertices.

Figure 11: Graph density vs. time in the same ensemble as figure 10.

Figure 12: Degree distribution vs. time in the ensemble of figure 10 (upper left:

t=0

; upper right:

t=0.01

; lower left:

t=0.02

; lower right:

t=0.09

).

The stochastic matrix Laplacian system is defined by a Poisson rate of change for each possible graph edge $(i,j)$ :

	$\displaystyle\mbox{rate}(e_{{ij}}\to 1-e_{{ij}})$	$\displaystyle=\|\left[\mathbf{L}_{m}(\mathbf{e})\mathbf{e}\right]_{{ij}}\|$		(2)
		$\displaystyle=\|\sum _{k}(2e_{{ik}}e_{{kj}}-(e_{{ik}}+e_{{kj}})e_{{ij}})\|.$

At that rate, each edge $(i,j)$ is added if it is not present and removed if it is present. Given an undirected graph, corresponding to a symmetric matrix $\mathbf{e}$ , this rate is a symmetric function of $i$ and $j$ , so it does not matter in which order they are considered.

If edge $(i,j)$ is present, the signed rate of change $[\mathbf{L}_{m}(\mathbf{e})\mathbf{e}]_{{ij}}$ is negative or zero, and if the edge is not present this rate is positive or zero. As a consequence, the instantaneous expected change in any given matrix entry $e_{{ij}}$ corresponds to the matrix Laplacian differential equation:

\displaystyle\lim _{{\Delta t\to 0}}\mathbb{E}\left[\Delta e_{{ij}}\right]

\displaystyle=\left[\mathbf{L}_{m}(\mathbf{e})\mathbf{e}\right]_{{ij}}\Delta t.

It follows that the number of edges $\sum _{{i,j}}e_{{ij}}$ is a martingale:

	$\displaystyle\lim _{{\Delta t\to 0}}\mathbb{E}\left[\sum _{{i,j}}e_{{ij}}(t+\Delta t)\right]$	$\displaystyle=\sum _{{i,j}}\left[e_{{ij}}(t)+\left(\mathbf{L}_{m}(\mathbf{e}(t))\mathbf{e}(t)\right)_{{ij}}\Delta t\right]$
		$\displaystyle=\sum _{{i,j}}e_{{ij}},$

and consequently the graph density is also a martingale.

This process’s absorbing states are the equilibria of the ODE, complete graphs and unions of disjoint complete graphs. These are graphs with transitivity of 1, and in simulations (e.g., figures 8, 9, and 10) the transitivity tends to increase quickly from 0 and approach 1 asymptotically.

It follows from the martingale that the graph’s initial density is an upper bound on the probability of ending at a single complete graph:

		$\displaystyle=\mathbb{E}[\text{final density}]$
		$\displaystyle=\mathbb{P}[\text{fixation at complete graph}]$
		$\displaystyle\quad\quad\mbox{}+\sum _{{\text{other final states }e^{}}}\text{(density of }e^{}\text{)}\mathbb{P}[\text{fixation at }e^{*}]$
		$\displaystyle\geq\mathbb{P}[\text{fixation at complete graph}].$

4 Discussion

The stochastic operation can probably be used for a short time to produce a graph that has a moderate amount of clustering and looks something like the original graph. The laplacian operation distorts the degree distribution (you can’t really see in the figures, but it adds a whole lot of zero-degree vertices), but maybe it’s possible to “tune” the transitivity and degree distribution at the same time?

These procedures may be useful for someone. Or they may inspire other useful ideas. Maybe the general idea of linear operations on entire graphs is a productive one. Maybe it’s worth generalizing from graphs to other metric spaces…

Some readers might wonder what these changes do to the spectrum of the graph’s Laplacian as the graph changes?

My main question right now is what does this manuscript need before I upload it to arXiv. I also want to send it to a journal, but that’ll be a month or more from now. I think i need to improve the intro, and add more citations, and obviously I need to say something reasonable here in the conclusions. I would also like to write something more about the stochastic process, and the figures should probably be improved and discussed more. What do you think?

5 Acknowledgements

[To do: Mention funding here.]

Supplemental material, including the source code for the above experiments and further unpublished results, is available at the first author’s online open lab notebook [1].

References

L. Worden (2011)
Reflexive Laplacian operations on graphs,

Lee Worden Research Wiki.

Note: \hrefhttp://lalashan.mcmaster.ca/theobio/worden/index.php/Laplacian_Paperhttp://lalashan.mcmaster.ca/theobio/worden/index.php/Laplacian_Paper

Cited by: 5.

ms.bib

\usepackage{hyperref}
@Article{Worden-wiki,
  author =    {Worden, Lee},
  title =     {Reflexive {Laplacian} operations on graphs},
  journal = {{Lee} {Worden} {Research} {Wiki}},
  year =      {2011--2012},
  note =      {\href{http://lalashan.mcmaster.ca/theobio/worden/index.php/Laplacian\_Paper}{http://lalashan.mcmaster.ca/theobio/worden/index.php/Laplacian\_Paper}}}

Makefile

ifndef RESOURCES
include ../.workingwiki/makefile
else
 
default: ms.pdflatex.pdf
 
paper-only :
	$(MAKE) -f $(RESOURCES)/makefile PAPER_ONLY=1 ms.pdflatex.pdf
 
.PHONY: default paper-only
 
ifndef PAPER_ONLY
figures/% : $(FIGURES)/% figures
	cp -f $(FIGURES)/$* $@
 
$(FIGURES)/% : /proc/uptime
	@$(MAKE) -C $(FIGURES) -f $(RESOURCES)/makefile $*
 
figures:
	mkdir $@
endif
 
.SECONDARY: $(patsubst %,$(FIGURES)/%,deterministic.snapshot.0.png deterministic.snapshot.0.2.png deterministic.snapshot.10.png deterministic.weights-vs-time.png sampling.snapshot.0.png sampling.snapshot.0.2.png sampling.snapshot.10.png sampling.transitivity.png sampling.ensemble.density.png sampling.ensemble.degree.distribution.0.png sampling.ensemble.degree.distribution.0.01.png sampling.ensemble.degree.distribution.0.02.png sampling.ensemble.degree.distribution.0.09.png stochastic.snapshot.0.png stochastic.snapshot.0.1.png stochastic.snapshot.0.99.png stochastic.transitivity.png stochastic.ensemble.density.png stochastic.ensemble.transitivity.png stochastic.ensemble.degree.distribution.0.png stochastic.ensemble.degree.distribution.0.01.png stochastic.ensemble.degree.distribution.0.02.png stochastic.ensemble.degree.distribution.0.09.png)
 
endif

	$\displaystyle\sum _{{ij}}x_{{ij}}$	$\displaystyle=\|G_{s}\|^{2}u+2\|G_{s}\|\|G_{t}\|z+\|G_{t}\|^{2}v$
		$\displaystyle=(-2\|G_{s}\|^{2}\|G_{t}\|k_{s}/\lambda+2\|G_{s}\|\|G_{t}\|-2\|G_{s}\|\|G_{t}\|^{2}k_{t}/\lambda)z$
		$\displaystyle=-2\|G_{s}\|\|G_{t}\|z\frac{\|G_{s}\|k_{s}-(\|G_{s}\|k_{s}+\|G_{t}\|k_{t})+\|G_{t}\|k_{t}}{\|G_{s}\|k_{s}+\|G_{t}\|k_{t}}$
		$\displaystyle=0,$

Laplacian Paper

From Worden

[log, pdf]Reflexive graph operations for clustering

Abstract.

1 Introduction

1.1 Laplacian diffusion of graph edges

2 Methods

2.1 Deterministic diffusion of weighted graph edges

2.2 Generation of unweighted graphs by Bernoulli sampling

2.3 Stochastic diffusion of unweighted graph edges

3 Results

3.1 Deterministic diffusion of weighted graph edges

3.1.1 Linear algebra of $\mathbf{L}_{m}$

3.1.2 Equilibria of the matrix Laplacian differential equation

3.1.3 Stability of equilibria

3.1.4 A positive eigenmatrix

3.2 Generation of unweighted graphs by Bernoulli sampling

3.3 Stochastic diffusion of unweighted graph edges

4 Discussion

5 Acknowledgements

References

Views

Personal tools

Navigation

Search

Projects

Toolbox

help links

Laplacian Paper

From Worden

[log, pdf]Reflexive graph operations for clustering

Abstract.

1 Introduction

1.1 Laplacian diffusion of graph edges

2 Methods

2.1 Deterministic diffusion of weighted graph edges

2.2 Generation of unweighted graphs by Bernoulli sampling

2.3 Stochastic diffusion of unweighted graph edges

3 Results

3.1 Deterministic diffusion of weighted graph edges

3.1.1 Linear algebra of Lm\mathbf{L}_{m}

3.1.2 Equilibria of the matrix Laplacian differential equation

3.1.3 Stability of equilibria

3.1.4 A positive eigenmatrix

3.2 Generation of unweighted graphs by Bernoulli sampling

3.3 Stochastic diffusion of unweighted graph edges

4 Discussion

5 Acknowledgements

References

Views

Personal tools

Navigation

Search

Projects

Toolbox

help links

3.1.1 Linear algebra of $\mathbf{L}_{m}$