# Efficient Evolutionary Models with Digraphons

Abhinav Tamaskar1111Corresponding Author. Email<EMAIL_ADDRESS>Bud
Mishra1
1Courant Institute of Mathematical Sciences, NYU

###### Abstract

We present two main contributions which help us in leveraging the theory of
graphons for modeling evolutionary processes. We show a generative model for
digraphons using a finite basis of subgraphs, which is representative of
biological networks with evolution by duplication[20]. We show a simple MAP
estimate on the Bayesian non parametric model using the Dirichlet Chinese
restaurant process representation, with the help of a Gibbs sampling algorithm
to infer the prior. Next we show an efficient implementation to do simulations
on finite basis segmentations of digraphons. This implementation is used for
developing fast evolutionary simulations with the help of an efficient 2-D
representation of the digraphon using dynamic segment-trees with the square-
root decomposition representation. We further show how this representation is
flexible enough to handle changing graph nodes and can be used to also model
dynamic digraphons with the help of an amortized update representation to
achieve an efficient time complexity of the update at $O(\sqrt{n}\log{n})$,
where $n$ is the number of nodes in the digraph.

##  1 Introduction

Graphical models are one of the most important tools used in machine
learning[11] and arise in most applications which involve pairwise
interactions, such as mutations in cancer evolution[26], protein networks[8,
25], hierarchical network models[5], influence in social networks[18, 21, 7],
population dynamics[1] and many more. In machine learning, there are various
techniques which forgo the use of these graphs and instead employ more
algebraic representations to take advantage of the underlying theories, such
as latent variable models[14], network or dynamic models[13], deep neural
networks (DNN)[24, 19], clustering models[2, 3]. The key advantage of the
later techniques is the reliance of the abundance of techniques developed in
linear algebra and optimized algorithms for doing fast implementations to get
efficient real world analysis.

One of the main motivations of this study is the case of evolutionary
populations, where the evolution is modeled as interactions between
individuals of the populations, such as mutations, genotypic variations and
phenotypic selection. In such cases, evolution of the gene regulatory network
(GRN) or the protein-protein interaction (PPI) network happens by specific
events, such as insertion, deletion, duplication, point mutations,
translocation and inversion[10, 6]. Of these the insertion, deletion and
duplication events have the most noticeable effect on the networks and have an
easily observable effect on the phenotype. If we think of genes as nodes in a
graph and gene interactions as edges, these events can be thought of as an
edge or node insertion, deletion or duplication.

When our network starts evolving and growing in size, we naturally think about
what the outcome of such a process would be. As our network keeps on
increasing in size, we need to extend our definition of a graph and naturally
arrive at the intuition of a _limit graph_ , i.e. a graph on an infinite
number of nodes, which we analyze not by looking at properties on non-empty
subsets of the graph. These dynamics can be represented in terms of limit
networks with the help of _digraphons_. A _digraphon_ is measurable function
$G:[0,1]^{2}\to[0,1]$. Given a digraphon $G$, there is a corresponding
countably infinite exchangeable graph $\mathcal{G}(\mathbb{N},G)$, with the
adjacency matrix $(\mathcal{G}_{ij})_{i,j\in\mathbb{N}}$ defined by the
generative model

$\displaystyle U_{i}$
$\displaystyle\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{Uniform}[0,1]~{}~{}\forall~{}i\in\mathbb{N}$
$\displaystyle\mathbb{G}_{ij}|U_{i},U_{j}$
$\displaystyle\stackrel{{\scriptstyle\text{ind}}}{{\sim}}\text{Bernoulli}(G(U_{i},U_{j}))$

Thus the digraphons are an ideal object to use for a generative model for
evolutionary networks. Any useful tool which is to be used for real world
analysis must also support _hypothesis refutability_ , which necessitates the
notion of similarity between digraphons. The space of digraphons have many
norms defined on it, such as the $l_{p}$ norm

$\left\lVert
G\right\rVert_{p}=\left(\int_{(i,j)\in[0,1]^{2}}\\!{(G_{ij})^{p}}\,\right)^{\frac{1}{p}}$

or the more interesting, _cut norm_ , which can be thought of as the maximum
dissimilarity in a bounded region,

$\left\lVert G\right\rVert_{\square}=\sup_{S,T\subseteq[0,1]}\int_{S\times
T}\\!{G}\,$

The cut metric is quite complex to calculate and can be approximated in finite
real world scenarios with the _agony_ metric which is exponentially faster to
compute[22].

These techniques for analyzing digraphons have recently been developed and
have yet to see a wider use in conjunction with the standard Bayesian
statistical tools prevalent in machine learning. The notions of limit graphs
and asymptotic behaviours of evolutionary models are very important in using
generative models from Bayesian statistics as the above model for digraphons
seems to suggest an intuitive method for reasoning an approximation of the
parameters.

The current black box learning from DNNs falls short of generating an
explainable hypothesis which is needed for refutability. Indeed such scenarios
have previously been observed such as the universal adversarial perturbations,
which have then been leveraged to design adversarial networks. But these still
fall short of an ideal answer.

The importance of these tools has been seen in many places, such as those used
for signaling games[12], population dynamics and biomolecular networks, mesh
network topologies[4], 3-D neural imaging reconstruction[9], etc. These
techniques work directly in complement to the notions from deep neural
networks by providing an explainable AI which can then be used as a
hyperparameter in designing the DNNs by affecting their layer hierarchies,
activation functions, dropout scenarios and memory length estimates.

##  2 Background - Graphons and Digraphons

Let $\mathcal{W}$ denote the space of all symmetric, bounded, measurable
functions $\mathcal{W}:[0,1]^{2}\to\mathbb{R}$. This is defined as the set of
all _kernels_ , reminiscent of the kernels used in Support Vector Machines. If
we restrict our attention to the set of functions $W\in\mathcal{W}_{0}$ such
that $0\leq W\leq 1$, we arrive at the space of _graphons_. If we drop the
condition of symmetry, we arrive at the space of digraphons.

For our purposes, we do not distinguish between the functions which are equal
almost everywhere, as with most analytical scenarios, it is not possible to
distinguish between such functions. We will soon see that this is actually not
the only equivalence we want to put on the functions, as we will need to also
equate a larger class of functions to each other for the sake of
exchangeability, which is important for statistical modeling.

The notion of graphons is an important one where we want to look at limits of
graph sequences. Such graph sequences arise in a variety of different natural
scenarios, such as social networks, recommendation systems for advertisements,
shopping, etc., and also in biological scenarios, such as genetic mutations,
evolutionary models, population dynamics. As a general rule of thumb, any
scenario where we have a dynamic population with potential for growth in event
space is a candidate for getting sequence of graphs.

Natural questions arise in how to analyze such a sequence. Does the growth
follow a pattern? Are there noticeable features of this graph that are
preserved across its growth? Does the graph sequence converge to any
discernible end object? The theory of graphons (and digraphons) tries to
answer many such questions in a rigorous form.

For us to have a notion of convergence and similarity, we need to start with a
notion of a distance between digraphs. There are many distances defined on the
space of digraphs, the distances introduced by the $L_{p}$ norms, nuclear
norms, etc.; we will restrict our attention to the more interesting _cut
distance_.

######  Definition. 1 (Cut Distance).—

For two directed graphs $G,G^{\prime}$, on the same set of vertices $V$, the
cut distance is defined as

$d_{\square}(G,G^{\prime})=\max_{S,T\subseteq
V}\frac{e_{G}(S,T)-e_{G^{\prime}}(S,T)}{|V|^{2}}$

where $e_{G}(S,T)$ denotes the number of edges between $S$ and $T$ in the
graph $G$.

If we let $d_{1}(G,G^{\prime})$ be the $L_{1}$ distance on the adjacency
matrices of $G$ and $G^{\prime}$ we get the inequality
$d_{\square}(G,G^{\prime})\leq d_{1}(G,G^{\prime})$. Hence we see that the two
distances give differing information. As an example, for two Erdos-Renyi
graphs with $p=1/2$, we get that $E[d_{1}(G,G^{\prime})]=1/2$ while
$E[d_{\square}(G,G^{\prime})]=\theta(1/\sqrt{n})$.

For unlabeled graphs on the same set of nodes, the intuitive extension to the
cut distance is to define it via equivalence on node relabelings, which turns
out to be the correct one. Let $G,G^{\prime}$ be two graphs on same number of
nodes, then the cut distance is overloaded as

$d_{\square}(G,G^{\prime})=\min_{\phi\in\text{Hom}(G,G^{\prime})}d_{\square}(\phi(G),G^{\prime}).$

Here we minimize over all homomorphims of $G$ into $G^{\prime}$. For
generalized measurable digraphons, we first define the _cut norm_

$\left\lVert W\right\rVert_{\square}=\sup_{S,T\subseteq[0,1]}\int_{S\times
T}\\!{W}\,$

Which can then be extended to the cut distance as
$d_{\square}{W,W^{\prime}}=\left\lVert W-W^{\prime}\right\rVert_{\square}$.
Similar to the finite case, we achieve the inequalities between norms

$\left\lVert W\right\rVert_{\square}\leq\left\lVert
W\right\rVert_{1}\leq\left\lVert W\right\rVert_{p}\leq\left\lVert
W\right\rVert_{\infty}\leq 1$

Again, similar to the finite case of the cut distance, where we have
“relabelings” via measure preserving homomorphisms of $\phi:[0,1]\to[0,1]$.

$d_{\square}(W,W^{\prime})=\inf_{\phi\in\text{Hom}([0,1])}d_{\square}(\phi(W),W^{\prime})$

It is important to note that finite graphs can be represented as a specific
case of a graphon by using step functions. For example, a digraph on $[n]$
with the adjacency matrix $A_{i,j}$ can be canonically viewed as a digraphon
$W$, with $W(\frac{i}{n},\frac{j}{n})=A[i,j]$. This allows us to treat even
finite graphs as a graphons and simplify our analysis, where we no longer have
to distinguish between sequences of graphs vs sequence of graphons.

We see that any measure preserving transformation of a digraphon has zero cut-
distance to the original. An important theorem states that the only digraphons
which have cut distance zero are the ones under measure preserving
transformations of the original or of one which is equal almost everywhere,
termed as weakly isomorphic pairs.

######  Theorem. 2.1 (Weak isomorphism theorem[16]).

Let $W,W^{\prime}$ be two digraphons, then $d_{\square}(W,W^{\prime})=0$ if
and only if there exists a digraphon $Z$, such that $W=Z$ almost everywhere
and $W^{\prime}$ is a measure preserving homomorphism of $Z$222In particular
this also covers the case where $W=Z$ everywhere..

This result is one of the most important ones for the analysis of digraphons,
as this gives confidence to our sampling algorithms. Indeed, the proof of this
theorem itself relies on the canonical digraphon and sampling state
introduced. Then we generalize the distance for two digraphons $W,W^{\prime}$
by round robin chasing across the commutative diagram,

$d_{\square}(W,\text{Sample}(W))\leftrightarrow
d_{\square}(W,W^{\prime})\leftrightarrow
d_{\square}(W^{\prime},\text{Sample}(W^{\prime}))$

######  Lemma. 2.1 (Convergence in norm).

Let $W_{n},n=1,2,\ldots,$ be a sequence of digraphons such that $\left\lVert
W_{n}\right\rVert_{\square}\to 0$. Then for all dikernels $Z$, we have that
$\left\lVert W_{n}Z\right\rVert_{\square}\to 0$.

###  2.1 Sampling

One of the important parts of digraphons are the guarantees that a finite
sampling is going to converge to the underlying digraphon. The sampling works
as follows.

Given a digraphon $W$ and an ordered set
$S=(x_{1},\ldots,x_{n}),x_{i}\in[0,1]$, we create a weighted digraph $H(S,W)$
on the node set $[n]$ with the edge weights $H(i,j)=W(x_{i},x_{j})$. Now from
such an $H$ we can create a random simple unweighted digraph by trying to
sample $G$ and adding an edge $G(i,j)$ with probability $H(i,j)$.

For example, if $W$ is the uniform function with $W(i,j)=p,0\leq p\leq 1$, we
would get the standard Erdos-Renyi graphs with probability $p$. If $W=W_{G}$,
the canonical digraphon for a digraph $G$, then if we sample $k$ points from
$W_{G}$, it is the “almost” the same as calculating a random subgraph of $G$.
It is not the same as we might have sampled $x_{i},x_{j}$ from the same step
in $W_{G}$. To have an exact subgraph sampling, we need to condition on the
fact that $x_{i},x_{j}$ need to be from different steps in $W_{G}$. In
particular, we are removing sequences$(x_{1},\ldots,x_{n})$, with repetitions,
which has $\binom{k}{2}$ such sequences, and hence a measure of
$\frac{\binom{k}{2}}{n}$. This gives us that the average distance between a
randomly chosen subgraph, $R(k,G)$ and a randomly sampled
digraphon$R(k,W_{G})$ is

$d(R(k,G),R(k,W_{G}))\leq\binom{k}{2}\frac{1}{n}$

Now that we know how to sample, we can start looking at how sampling helps in
parameter estimation. For such a scenario, we want to start with a notion of
“good” parameters, which we can hope to be estimable. As it will turn out,
most of the real world scenarios are going to be good and can be estimated
using sampling. We wish to achieve some notions similar to the limit theorems
for classical statistics which will give us confidence on doing real world
analytics using EM or MAP algorithms for estimations.

A _reasonably smooth_ graph parameter is defined as a function $f(G([S]))$,
which is to say of a sampling of a digraphon, which satisfies that
$|f(G)-f(G^{\prime})|<1$, for two graphs $G,G^{\prime}$ on the same set of
nodes, whose edges differ only for a single vertex. We then achieve a sample
concentration theorem

######  Theorem. 2.2 (Sample concentration theorem for digraphons[15]).

Let $f$ be a reasonably smooth graph parameter, and let $W$ be a digraphon,
$k>1\in\mathbb{N}$. Let $f_{0}=\operatorname{\mathbf{E}}[f(R(k,W_{G}))]$, then
for all $t>0$,

$\Pr\left[f(R(k,W_{G}))>f_{0}+\sqrt{2tk}\right]<e^{-t}$

This theorem gives credibility to the fact that our intuitive sampling
algorithms are going to be working correctly on standard simulations. In fact,
we can achieve an even better result which states.

######  Theorem. 2.3 (Cut distance confidence[15]).

Let $k>1\in\mathbb{N}$ and let $G$ be a digraph on $k$ nodes. Then with
probability at least $1-\exp{(-\frac{k}{2\log k})}$, we have that

$d_{\square}(G,R(k,W_{G}))\leq\frac{20}{\sqrt{\log k}}$

Now with such a result, we can employ our standard statistical tools to carry
out prior/posterior estimations and simulations using generative models.
Indeed, we develop a maximum a posteriori (MAP) estimation algorithm for the
Dirichlet priors on digraphon generative models and due to this result, we can
be somewhat confident in the fact that we are achieving a good result, subject
to proper maximum optimizations in the log likelihood estimation.

##  3 Evolution by Duplication

In many scenarios, preferential attachment happens at a larger scale than a
single vertex, where multiple parts, aka clusters, of the network get
duplicated. This process is hard to model in a one shot setting where we have
to duplicate the parts one at a time. For added robustness, we expand the
definition of _preferential attachment_ to allow duplicating subgraphs.

We specifically want to model graphs for evolutionary events, such as
graphical models and causal networks. Each node represents an event that can
take place, while each edge (potentially weighted) represents the _influence_
from an event to another. We call this the _event graph_ and represent this as
a weighted adjacency matrix. An example of such a network that is the Suppes
Bayes causal network (SBCN).

We define evolution by duplication as an extension of the preferential
attachment model, where we allow extending the graph by duplicating a larger
subgraph. This subsumes the original case where preferential attachment of a
single node and allows for a larger, more robust evolutionary model.

######  Definition. 2 (Evolution by Duplication).—

Let $G$ be a directed graph, and let
$X=[a_{i},\ldots,a_{j}]\times[b_{k},\ldots,b_{l}]\subset[1,\ldots,n]^{2}$. A
preferential attachment of $X$ on $G$, with the weight function $\theta$, is
the new digraphon $G^{\prime}$, defined by

$G^{\prime}=(1-\theta(X))\cdot G+\theta(X)\cdot G|_{X}$

Typically, we want to randomly select segments which are going to be attached.
This is carried out with the help of a weight function $\theta$, which is
inversely proportional to the measure of the attached segment,
$\theta(X)\propto\frac{1}{\int_{X}^{~{}}\\!{G}\,}=\frac{1}{\sum_{\begin{subarray}{c}x\in[a_{i},\ldots,a_{j}]\\\
y\in[b_{k},\ldots,b_{l}]\end{subarray}}w(a_{x},b_{y})}$. This weighting
implies that smaller segments are easier to attach than larger segments, which
mimics the biological characteristics, wherein it is easier to duplicate
smaller, simpler networks over complex multi-path networks.

Figure 1: Segment basis for a digraph.
The colored rectangles denote _segments_ for a digraph, which can be
duplicated by _preferential attachment_. It is not necessary for the segments
to be non-intersecting, the overlapped regions will be attached multiple times
through different segments. Each segment also has an associated weight, which
influences the evolutionary priority of duplicating that segment.

Due to the exponentially large size of the segment space, it is preferred to
have a smaller finite collection of segments of interest. In biological
networks, such as those for cancer somatic mutations, we typically want to
restrict our attention to either driver genes or those known to interact with
cancer mutations up to some extent. For example, The Cancer Genome Atlas
(TCGA) program or the Catalogue of Somatic Mutations in Cancer (COSMIC)
databases are the main places which give information about important genes.

######  Definition. 3 (Segmented Digraph).—

A _segmented_ digraph, fig. 1, is a digraph with a finite collection of
weighted segments $(X_{i},w_{i})$, $X_{i}=[a_{i},b_{i}]\times[c_{i},d_{i}]$,
with weights $w_{i}\in\mathbb{R}_{+},\sum w_{i}=1$.

We can also allow dynamic sized vertices by allowing an attachment to
introduce a node, thereby adding a row or a column to our adjacency matrix.
For our current analysis, we restrict ourselves to static sized attachments,
which only affect weights of the graph.

It is important to note that this is an exchangeable process; relabeling the
segments, $X_{i}$, does not lead to a different simulation, we only depend
upon the weights of the segments. Thus we can leverage all the theoretical
properties of digraphons, such as asserting that any such evolutionary
process, wherein we do simulations either by attachments to make larger
graphs, or by performing weighted boosting, is going to converge to a final
digraphon.

abcdefghabcdefgh Figure 2: Evolution by subgraph attachment via duplication.
Preferential attachment of the segment $[d,f,g,h]^{2}$. The same effect can be
seen for more than one segment attachment, as evidenced by preferentially
attaching $[d,f,g]\times[f,g,h]$.

##  4 Finite Modeling and Implementations

To model such an evolutionary process is non trivially complex, as the graphs
can have large number of nodes and edges. Naive algorithms for doing
_preferential attachment_ of a segment $[i,j]\times[k,l]$, have a time
complexity of $(j-i)\cdot(k-l)=O(n^{2})$. In addition to that, to do a
weighted sample, we need to get the current weight of a segment, which would
also take $O(n^{2})$. This strategy is feasible in cases of small segments,
but has an undesirable asymptotic behaviour, which we improve upon. In this
section, we detail an efficient data structure for simulating the preferential
attachment framework on a digraph.

There are two key operations that we want our model to access:

  1. 1.

Get the current weight of a segment $X$.

  2. 2.

_Preferentially Attach_ a segment, $X$, with a particular weight, $\theta(X)$.

We note that we can perform (2) in two steps - (a) Multiply the whole graph by
$(1-\theta(X))$, (b) Multiply the segment $X$ by $\frac{1}{1-\theta(x)}$.

Hence, our data structure needs to support the following operations:

  1. 1.

Multiply a segment $[a,b]\times[c,d]$ by some value $c$.

  2. 2.

Return the sum of all elements in a segment $[a,b]\times[c,d]$.
For a finite model implementation, we assume that the size of the digraph is
fixed.

$T_{0}$$T_{1}$$T_{2}$ Figure 3: Graph data structure for preferential
attachment.
We break the graph into $\sqrt{n}$ groups along the rows and create a segment
tree, $T_{i}$, for each of these groups. We see that any preferential
attachment that is done can be broken down into $\sqrt{n}$ parts, where we do
updates for each tree, with a total of $\sqrt{n}$ updates. The overlaps are
handled efficiently using lazy propagation in the segment tree, which allow us
to do the weight updates in $O(\log n)$ for each individual tree.

We represent the digraph $G$ using the square-root decomposition along the
rows and a segment tree along the columns for each group of $\sqrt{n}$ rows,
fig. 3.
Each $T_{i}$ represents a collection of $\sqrt{n}$ rows with a segment tree.
An update of a segment $X=[a,b]\times[c,d]$ can span across segment trees. The
square root decomposition facilitates fast updates in $O(\sqrt{n}\log{n})$
time.

For each segment tree, it is possible to do _range sum_ together with _range
multiply_ using _lazy propagation_ in $O(\log{n})$, algorithm 1 and algorithm
2.

To prove that the whole data structure works in the time complexity described
above, we break the analysis into two parts. First, we show that in an
individual segment tree, we can solve range sum and range multiplication in
$O(\log n)$. Then we show how to extend this across rows by grouping multiple
trees together.

###  4.1 Segment Tree - Lazy Multiply and Sum

The per row simplification of this problem boils down to

######  Problem. 4.1 (Lazy Multiply and Sum).

Given an array of numbers, $[{a}_{1},\ldots,{a}_{n}]$, perform the following
operations in $O(\log n)$

  1. 1.

Find the sum of all numbers in a contiguous range $[i,j]$.

  2. 2.

Multiply all numbers in a contiguous range $[i,j]$ by some value $k$.

$[{a}_{1},\ldots,{a}_{n}]-\\{sum=s_{0},mult=m_{0}\\}$$[{a}_{1},\ldots,{a}_{n/2}]-\\{s_{1},m_{1}\\}$$[{a}_{n/2+1},\ldots,{a}_{n}]-\\{s_{2},m_{2}\\}$$\vdots$[]
Figure 4: Segment tree for lazy propagation.
For each node in the segment tree, we store the subtree range along with two
additional parameters, current sum and propagated multiplicand. The
multiplicand is not propagated fully until a subrange is updated in a
subsequent query or until a subrange weight is queried. And even in those
cases, we only propagate it optimally with lazy propagation.

We solve this by creating a segment tree of nodes, fig. 4 ,with each node
containing the following metadata:

  * •

_low, high_ \- the lower and upper end points of the range of the node

  * •

_sum_ \- current sum of all values in the range of the node

  * •

_mult_ \- current multiplicand not yet propagated to lower nodes

  * •

_left, right_ \- left and right children nodes of current node

Algorithm 1 Segment Tree - Lazy Multiply

1:function LazyMultiply(Segment Tree Node $node$, int $l$, int $r$, int $c$)

2: # Multiplies the part of the range $[l,r]$ contained inside $node$ by $c$

3: # Time complexity = $O(\log{n})$

4: if $node.low>r$ or $node.high<l$ then

5: # No intersection, nothing to do

6: return

7: elseif $node.low\geq l$ and $node.high\leq r$ then

8: # Fully contained, multiply current multiplicand by $c$

9: $node.mult\mathrel{*}=c$

10: $node.sum\mathrel{*}=c$

11: else

12: LazyMultiply($node.left$, $node.low$, $node.high$, $node.mult$)

13: LazyMultiply($node.right$, $node.low$, $node.high$, $node.mult$)

14: $node.mult=1$

15: LazyMultiply($node.left$, $l$, $r$, $c$)

16: LazyMultiply($node.right$, $l$, $r$, $c$)

17: $node.sum=node.left.sum+node.right.sum$

18: end if

19:end function

The lazy multiplication algorithm, algorithm 1, works by manipulating the
$mult$ parameter to keep track of the accumulated multiplication and only
propagating it when intersecting ranges are updated.

######  Theorem. 4.1 (Lazy Multiply).

The running time of algorithm 1 is $O(\log{n})$ for all ranges $[l,r]$.

###### Proof.

Notice that the recursion stops at any node at which the range of the node is
full contained inside the range to be updated. For example, if our tree is on
$[1,\ldots,8]$, we will have $15$ total nodes,

$[1,8],[1,4],[5,8],[1,2],[3,4],[5,6],[7,8],[1,1],\ldots,[8,8]$

. If we wish to do a range multiplication on $[4,7]$, the algorithm will stop
at the top most nodes possible which unify up to our desired range, which in
this case - $[4,4],[5,6],[7,7]$, which is strictly smaller than $4$, the range
to be updated.

Notice that we will only ever update a maximum of two nodes of the same
length, and they will never be neighbours, as the recursion would instead stop
at the parent node. Hence, to represent our range $[l,r]$, as a unification of
$k$ ranges - $[l,r]=\cup_{i=1}^{k}[x_{i},y_{i}]$, we can only have $k<2\log
n$, as the maximum size of a node’s subrange is $n/2$.

This argument shows that we will only ever visit $O(\log n)$ nodes, proving
our running time proposition. ∎

We can reach the desired worst case of $\sim 2\log n$ for a tree of range
$[1,2^{n}]$ and a range update for $[2,2^{n}-1]$.

Algorithm 2 Segment Tree - Lazy Sum

1:function LazySum(Segment Tree Node $node$, int $l$, int $r$)

2: # Returns the sum of the part of the range $[l,r]$ in $node$

3: # Time complexity = $O(\log{n})$

4: if $node.low>r$ or $node.high<l$ then

5: # No intersection

6: return 0

7: elseif $node.low\geq l$ and $node.high\leq r$ then

8: # Fully contained

9: return $node.sum$

10: else

11: LazyMultiply($node.left$, $node.low$, $node.high$, $node.mult$)

12: LazyMultiply($node.right$, $node.low$, $node.high$, $node.mult$)

13: $node.mult=1$

14: return LazySum($node.left$, $l$, $r$) + LazySum($node.right$, $l$, $r$)

15: end if

16:end function

######  Theorem. 4.2 (Lazy Sum).

The running time of algorithm 2 is $O(\log n)$ all ranges $[l,r]$.

The proof is identical to that of the lazy multiplication, where we only look
at the top level nodes.

###  4.2 Lazy Attach

To utilize the previous algorithms for a 2-D structure, we have to make
certain modifications.
We first create a segment tree for each row and then group them together in
groups of size $\sqrt{n}$. For each group, create a parent segment tree, with
the same metadata, except that for each node, store the parent.sum =
$\sum_{i\in\text{rows}}$ row.sum. The parent tree stores the aggregate
information across all rows, which allows us to do range queries across the
whole group in $O(\log n)$.
It is vital to note that this is only for the whole group and not a subset of
the group. To query for a subset of the group, we have to go to each
individual row and query its segment tree. Due to the fact that each group is
of size $\sqrt{n}$, we are guaranteed that each operations only touches a
maximum of $\sqrt{n}$ groups, out of which only two groups need to ever do
individual row updates, as show in algorithm 3.

Algorithm 3 Preferential Attachment on a Digraph

1:function PreferentialAttach(Digraph $G$, Segment $X$, Weight k)

2: # Modifies $G$ in place by attaching the 2-D segment $X=[a,b]\times[c,d]$

3: # Time complexity = $O(\sqrt{n}\log{n})$

4: for $T_{i}\in[1,\sqrt{n}]$

5: LazyMultiply($T_{i}.root,1,n,1-k$)

6: end for

7: for $T_{i}\in[a,b]$

8: # These trees are fully inside the range and can be updated as a group

9: LazyMultiply($T_{i}.root$, $c$, $d$, k)

10: end for

11: for boundary trees $T_{a}$ and $T_{b}$

12: # These two trees have partial intersection with the range $[a,b]$ and
must be updated manually

13: # Each individual row is also represented as a segment tree

14: for $R_{i}\in T_{a},T_{b}$ and $i\in[a,b]$

15: LazyMultiply($R_{i}.root$, $c$, $d$, 2)

16: end for

17: end for

18:end function

####  4.2.1 Real World Optimizations

There are some factors that can be considered for optimizations.

  * •

If it is known that segment sizes to be updated are within a certain bounded
width $w$, then it is possible to create segment groups of size $w$ instead of
$\sqrt{n}$, this implies, that there will only ever be $O(w\log n)$ maximum
time.

  * •

For really small width segments, it can also be possible to use data
structures such as a Quad Tree or a k-D tree.

  * •

If working with large segments, but a small number of such segments, it is
beneficial to look at binary space partition (BSP) trees and pre-partition
segments into non intersecting parts. This can quickly become complex, with a
growth in number of segments. BSPs are used in computer vision to do
segmentation, which allows us to use highly optimized implementations, if such
a scenarios is feasible.

##  5 Learning and Inference

Sampling of a digraphon is done by the Chinese Restaurant Process(CRP), which
is a staple tool to model the Dirichlet distribution[23, 17]. The CRP process
aims to model how people are assigned tables when sitting at a shared seating
restaurant. With higher probability, people wish to sit next to others for a
more pleasant experience, while with a lower probability, they wish to get a
new table. This is formalized as follows.

Let $\alpha\in(0,1)$ be a hyperparameter. At any point of time $n$, let us
have $k$ groups of size $[g_{1},\ldots,g_{k}]$, which are a partition of
$[1,\ldots,n]$. At time $n+1$, we wish to assign a group to the element $n+1$
which is done as follows:

  * •

With probability $\frac{(1-\alpha)|g_{i}|}{n+1}$, $n+1$ is assigned to group
$i$.

  * •

With probability $\frac{\alpha}{n+1}$, $n+1$ is assigned to a new group
$g_{k+1}$

As we scale $n\to\infty$, we achieve a distribution of the set $\mathbb{N}$,
which is the Dirichlet distribution with scaling parameter $\alpha$. The CRP
representation is very useful convenient when performing finite sampling and
parameter estimation. Another advantage of the CRP representation is the ease
of computation and the fast simulations, which are important for larger
models.

###  5.1 Sampling a Digraphon

Similar to Bayesian statistical models, we need a generative model for a
digraphon to be able to get insights using parameter estimation. The most
common generative models used are the Dirichlet prior, which do leverage using
the CRP model.
Let $\alpha$ be the hyper parameter affecting the CRP model to get the
clustering assignment of the vertices and let $\beta$ be another hyper
parameter for the standard Dirichlet process, such as the gamma
representation.
The generative model for the digraphon is as follows, where we generate a
digraph on $n$ vertices:

  1. 1.

Draw clustering assignments, $\zeta$, for each vertex,

$\zeta\sim CRP(\alpha)$

  2. 2.

Draw weights for the edges, for each pair of groups $r\neq s$

$\eta_{r,s}|\beta\sim\text{Dirichlet}(\beta)$

  3. 3.

Set the edge using the measure of the partition

$\mathcal{G}_{ij}=\text{Categorical}(\eta_{\zeta_{i},\zeta_{j})}$

We see that this is an exchangeable process as well, as the clustering
assignments are generated irrespective of the actual labels of the vertices.
This reasoning implies that as the number of vertices $n\to\infty$, the
generated digraph converges to the digraphon.

###  5.2 MAP inference

Let $\mathcal{G}$ be the final digraph generated from the digraphon generative
model. Our aim is to infer the weights $\eta_{r,s}$.

The likelihood that $\mathcal{G}$ is sampled is given by

$\Pr(\mathcal{G}|\zeta)=\prod_{r\neq s}(\zeta_{r,s})^{m_{r,s}}$

where $m_{r,s}$ denotes the number of edges from cluster $r$ to cluster $s$
and we assume no self loops for simplicity.

$\Pr(\zeta|\alpha)=\frac{\prod_{r\neq
s}(\zeta_{r,s})^{(\alpha-1)}}{\mathbf{B}(\alpha)}$

where
$\mathbf{B}(\omega)=\frac{\prod_{i}\Gamma(\omega_{i})}{\Gamma(\sum_{i}\omega_{i})}$
is the multivariate beta function.

$\displaystyle\Pr(\mathcal{G}|\alpha)$
$\displaystyle=\int\\!{\Pr(\mathcal{G}|\zeta)\cdot\Pr(\zeta|\alpha)}\,$
$\displaystyle=\int\\!{\frac{\prod_{r\neq
s}(\zeta_{r,s})^{m_{r,s}}\prod_{r\neq
s}(\zeta_{r,s})^{(\alpha-1)}}{\mathbf{B}(\alpha)}}\,$
$\displaystyle=\frac{1}{\mathbf{B}(\alpha)}\prod_{r\neq
s}\int\\!{\zeta_{r,s}^{m_{r,s}+\alpha-1}}\,\mathrm{d}\zeta_{r,s}$
$\displaystyle=\frac{1}{\mathbf{B}(\alpha)}\prod_{r\neq
s}\mathbf{B}(m_{r,s}+\alpha)$

We can remove the constant $\frac{1}{\mathbf{B}(\alpha)}$ and maximize the
negative log likelihood
$\sum_{r\neq s}\log\mathbf{B}(m_{r,s}+\alpha)$ at

$\zeta_{r,s}=\frac{m_{r,s}+\alpha}{\displaystyle\binom{n}{2}\cdot\alpha+\sum_{r\neq
s}m_{r,s}}$

##  6 Discussion

As final concluding remarks, we have seen how to use digraphons as generative
models for directed graphs for evolutionary models. We have further developed
a robust modeling data structure for fast simulations which is easily
extendable for other evolutionary mechanisms, as the data structure allows for
generic operations which can be used for other world models. This opens up
further testing grounds for hypothesis checking by comparing simulations to
real world dynamics.

There are many extensions that have yet to be explored using the theory of
digraphons, of which an important one is the use of the various distance
metrics on the digraphon space. There are many metrics, such as the ones
induced by the $L_{p}$ norm , nuclear norm, and the more interesting cut-
distance, which are used. Some of the more interesting questions are presented
below.

(1) Can we use the metric on the digraphon space to measure similarity of
models for two distinct populations?
An example of such a scenario would be, given data for two distinct
populations, we wish to know if they have evolved from a common ancestor
population. If so, how far back did they diverge? Can we quantify the
divergence of populations? Even more thoroughly, can we find the evolutionary
pathway used by the populations to reach the current state?

(2) Do the metrics induce a EM-type convergence for learning algorithms?
Given multi-dimensional data for a population, we wish to stratify it into
sub-populations with individuals having a similar evolutionary patterns. This
problem is reminiscent of k-means clustering where we wish to use the
digraphon metric to perform the clustering. There has been extensive work in
the analysis of Euclidean norms for showing convergence (if only to a local
minima for the error function), which would be important to translate to the
digraphon spaces. The work by Lovász, et al.,[16] has shown many convergence
results which may be helpful for such scenarios.

## References

  * [1] Rob J de Boer.  Modeling population dynamics: A graphical approach, 2012.
  * [2] Åke Brännström and David JT Sumpter.  The role of competition and clustering in population dynamics.  Proceedings of the Royal Society B: Biological Sciences, 272(1576):2065–2072, 2005.
  * [3] Chibiao Chen, Eric Durand, Florence Forbes, and Olivier François.  Bayesian clustering algorithms ascertaining spatial population structure: a new computer program and a comparison study.  Molecular Ecology Notes, 7(5):747–756, 2007.
  * [4] Vin De Silva and Robert Ghrist.  Coverage in sensor networks via persistent homology.  Algebraic & Geometric Topology, 7(1):339–358, 2007.
  * [5] Walter Dempsey, Brandon Oselio, and Alfred Hero.  Hierarchical network models for exchangeable structured interaction processes.  Journal of the American Statistical Association, pages 1–43, 2021\.
  * [6] Douglas H Erwin and Eric H Davidson.  The evolution of hierarchical gene regulatory networks.  Nature Reviews Genetics, 10(2):141–148, 2009.
  * [7] Amit Goyal, Francesco Bonchi, and Laks VS Lakshmanan.  Learning influence probabilities in social networks.  In Proceedings of the third ACM international conference on Web search and data mining, pages 241–250, 2010.
  * [8] Marco Grzegorczyk.  Extracting protein regulatory networks with graphical models.  Proteomics, 7(S1):51–59, 2007.
  * [9] Xiaoling Hu, Li Fuxin, Dimitris Samaras, and Chao Chen.  Topology-preserving deep image segmentation.  arXiv preprint arXiv:1906.05404, 2019.
  * [10] David Juan, Florencio Pazos, and Alfonso Valencia.  Co-evolution and co-adaptation in protein networks.  FEBS letters, 582(8):1225–1230, 2008.
  * [11] Daphne Koller and Nir Friedman.  Probabilistic graphical models: principles and techniques.  MIT press, 2009.
  * [12] Travis LaCroix.  Evolutionary explanations of simple communication: Signalling games and their models.  Journal for General Philosophy of Science, 51(1):19–43, 2020.
  * [13] Simon A Levin.  Population dynamic models in heterogeneous environments.  Annual review of ecology and systematics, pages 287–310, 1976.
  * [14] John C Loehlin.  Latent variable models.  hillsdale, nj: erlbaum, 1987.
  * [15] László Lovász.  Large networks and graph limits, volume 60.  American Mathematical Soc., 2012.
  * [16] László Lovász and Katalin Vesztergombi.  Nondeterministic graph property testing.  arXiv preprint arXiv:1202.5337, 2012.
  * [17] Steven N MacEachern and Peter Müller.  Estimating mixture of dirichlet process models.  Journal of Computational and Graphical Statistics, 7(2):223–238, 1998.
  * [18] Winter A Mason, Frederica R Conrey, and Eliot R Smith.  Situating social influence processes: Dynamic, multidirectional flows of influence within social networks.  Personality and social psychology review, 11(3):279–300, 2007.
  * [19] Risto Miikkulainen, Jason Liang, Elliot Meyerson, Aditya Rawal, Daniel Fink, Olivier Francon, Bala Raju, Hormoz Shahrzad, Arshak Navruzyan, Nigel Duffy, et al.  Evolving deep neural networks.  In Artificial intelligence in the age of neural networks and brain computing, pages 293–312. Elsevier, 2019.
  * [20] Susumu Ohno.  Evolution by gene duplication.  Springer Science & Business Media, 2013.
  * [21] Garry Robins and Philippa Pattison.  Random graph models for temporal processes in social networks.  Journal of Mathematical Sociology, 25(1):5–41, 2001.
  * [22] Nikolaj Tatti.  Faster way to agony.  In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pages 163–178. Springer, 2014.
  * [23] Yee Whye Teh.  Dirichlet process., 2010.
  * [24] Suneetha Uppu and Aneesh Krishna.  Tuning hyperparameters for gene interaction models in genome-wide association studies.  In International Conference on Neural Information Processing, pages 791–801. Springer, 2017.
  * [25] Adriano V Werhli, Marco Grzegorczyk, and Dirk Husmeier.  Comparative evaluation of reverse engineering gene regulatory networks with relevance networks, graphical gaussian models and bayesian networks.  Bioinformatics, 22(20):2523–2531, 2006.
  * [26] Haitao Zhao and Zhong-Hui Duan.  Cancer genetic network inference using gaussian graphical models.  Bioinformatics and biology insights, 13:1177932219839402, 2019.