# Improved tracking of particles with highly correlated motion

Ella M. King1, Zizhao Wang2, David A. Weitz1,2, Frans Spaepen2, Michael P.
Brenner1,2 1 Physics Department, Harvard University  2 School of Engineering
and Applied Sciences, Harvard University

###### Abstract

Despite significant advances in particle imaging technologies over the past
two decades, few advances have been made in particle tracking, i.e. linking
individual particle positions across time series data. The state-of-the-art
tracking algorithm is highly effective for systems in which the particles
behave mostly independently. However, these algorithms become inaccurate when
particle motion is highly correlated, such as in dense or strongly interacting
systems. Accurate particle tracking is essential in the study of the physics
of dense colloids, such as the study of dislocation formation, nucleation, and
shear transformations. Here, we present a new method for particle tracking
that incorporates information about the correlated motion of the particles. We
demonstrate significant improvement over the state-of-the-art tracking
algorithm in simulated data on highly correlated systems.

††preprint: APS/123-QED

##  I  Introduction

Accurately tracking particle positions over time remains a major challenge in
advancing colloidal physics research. Single-particle level tracking is
critical to studies of the mechanisms of crystal nucleation[1], melting[2],
glass formation and dislocations [3], and a variety of other phenomena in
dense colloid systems. While advanced confocal microscopy and image analysis
techniques have improved the accuracy of the extraction of particle positions
from experiments, less progress has been made in accurate tracking of
particles exhibiting highly correlated particle motion, such as dense or
strongly interacting systems.

Existing tracking algorithms are based on the fundamental work of Crocker and
Grier [4], in which positions of particles undergoing Brownian motion are
linked between time steps by finding the closest candidates in successive
frames. Importantly, this method has two conditions for accurate
reconstruction of the tracks. First, the particle displacement between two
successive time steps should be much smaller than the inter-particle
separation. Second, each particle should diffuse independently. If these two
conditions are not met, the trajectory of each particle can often not be
distinguished from that of its neighbors. These conditions are violated in
systems of strongly interacting particles, which limits the accuracy of the
Crocker-Grier algorithm.

Though enhancements to Crocker and Grier’s original algorithm have been made,
such as handling uniform flow [5, 6], the fundamental requirement that
particle displacement be far less than the inter-particle separation remains a
significant limitation. This problem is particularly acute in 3D imaging due
to limitations in the frame rate.

In this work, we introduce a new algorithm that incorporates information about
correlations in particle motion, by using simulations to improve estimates of
tracking probability distributions. The algorithm achieves up to 33% more
accurate tracking than Crocker and Grier in simulated systems that go beyond
the bounds of the algorithm introduced by Crocker and Grier. We also
demonstrate robust results on noisy experimental data.

##  II Tracking Correlated Motion

###  II.1  The Crocker-Grier Algorithm

The core assumption of the Crocker-Grier tracking algorithm is that a given
particle at position $x_{0}$ will be found at a position $x$ after a time $t$
has passed with probability

$P(x,t)=\frac{1}{(4\pi Dt)^{\frac{d}{2}}}e^{-\frac{(x-x_{0})^{2}}{4Dt}}$ (1)

where $D$ is the diffusion coefficient and $d$ is the dimension. This is the
Green’s function solution to the single-particle diffusion equation [4]. The
tracking algorithm is accurate for systems in which particles diffuse
independently. However, for systems in which the motion in a group of
particles is correlated, the assumption underlying Crocker-Grier breaks down.

In particular, once multi-particle correlations are introduced, the single-
particle diffusion equation is no longer an accurate model of the system.

###  II.2  The Correlation Tracking Algorithm

This new tracking algorithm approximates multi-particle correlations using
simulated data to inform predictions. We refer to this method as “Correlation
Tracking.” Inspired by prior work showing that particle interactions can
directly tune dynamics [7], we use prior knowledge about particle interactions
to make more accurate predictions about the dynamics. By incorporating this
physical information into the algorithm, we are able to increase tracking
accuracy.

Consider first a pair of frames, $x_{0}$ and $x_{1}$. The frame $x_{i}$
contains a list of particle positions at time $i$. The aim of particle
tracking is to match particle labels between frame $x_{0}$ and frame $x_{1}$.
Because we know the interaction potential, we can simulate the motion
beginning at position $x_{0}$ and let it evolve for the same length of time as
between the pair of frames. These simulations are stochastic, and thus we run
them multiple times to extract $M$ replicates that approximate a distribution
of potential particle positions at time $x_{1}$. Using this set of potential
particle positions, we fit multivariate Gaussian distributions to subsets of
the particles.

Consider a set of $k$ particles (the $k-1$ particles nearest to a central
particle) in frame $x_{0}$. These particles are in a bath of $N-k$ particles.
To track these $k$ particles in the subsequent frame, we perform the
simulations described above. We refer to the resulting simulations as
$x_{1,\text{sim}}$. All $N$ particles are simulated in each of the $M$
replicates. The $M$ replicates provide mean positions $\mu_{i}=\langle
x_{1,\text{sim}}^{i}\rangle$ for each particle $i$ in the set of $k$
particles. Additionally, we compute correlations between each particle in the
set of $k$:

$C_{ij}=\langle(x_{1,\text{sim}}^{i}-\mu_{i})\cdot(x_{1,\text{sim}}^{j}-\mu_{j})\rangle$

Together, the computed means and correlations populate a multivariate Gaussian
distribution that describes the probability of finding a set of $k$ particles:

$P(x|x_{0},\alpha)=\frac{1}{\sqrt{(2\pi)^{d}|C|)}}\exp{-\frac{1}{2}(x-\mu)^{T}C^{-1}(x-\mu)}$
(2)

where $x$ is the set of positions of the $k$ particles we are interested in
tracking, $\alpha$ is the set of parameters describing the interaction
potential, $d$ is the dimension, and $\mu$ and $C$ are defined above.

The probability function $P(x|x_{0},\alpha)$ differs for each set of $k$
particles. To link particles in one frame to the next, we compute the
probability function for each set of $k$ particles and use it to determine
which set of $k$ particles in the next frame is most likely to be the same
set. Because it is computationally intractable to test every possible ordering
of $k$ particles in the frame, we select a set of candidate particles to test.
These candidate particles are the $p$ closest particles to the set of $k$ in
the first frame.

##  III Results on simulated data

We demonstrate our method in a several simulated systems. We focus on systems
that exhibit long-range attraction, as these are systems that do not meet the
conditions for accurate tracking with the Crocker-Grier algorithm.

In Fig. 1, we show the results of the Correlation Tracking method when applied
to a Morse pair potential of the form:

${}U(r)=D(e^{-2\alpha(r-r_{0})}-2e^{-\alpha(r-r_{0})})$ (3)

with $D=1$, $r_{0}=1$ and $\alpha=1$. In all the data shown below, the
algorithm was tested on simulated data with 100 particles undergoing Brownian
motion with periodic boundary conditions. The volume (or area) fraction of the
system was set to 0.5, with $r_{0}$ as the particle diameter. Molecular
dynamics simulations were performed using Jax-MD. When tracking sets of $k$
particles with $p$ candidate particles, we specified $k=5$ and $p=15$ for all
data. To estimate the correlation matrix for a given pair of frames, we
averaged over 10,000 independent simulations.

Figure 1: Comparison of the Correlation Tracking results to those of the
Crocker-Grier algorithm in simulated data using a Morse pair potential (eq.
3). The x-axis shows the mean particle displacement divided by the mean inter-
particle separation in units of $r_{0}$. This metric serves as a proxy for the
time and enables ready comparison between simulation and experiment. The
y-axis shows the tracking accuracy, or the fraction of correctly identified
particle labels. The uncertainties are standard errors, based on 25
independent simulations.

As shown in Fig. 1, we observe significant improvement in the tracking
accuracy for longer gaps between frames over the state-of-the-art algorithm in
both two and three dimensions.

We hypothesize that the success of our algorithm in this system is owed to the
large range of interaction. The greater the interaction range, the more
correlations can be established. To demonstrate this, we compare in Fig 2 the
accuracy of our method to that of the Crocker-Grier algorithm as a function of
the range of interaction, which is varied by tuning the parameter $\alpha$ in
the Morse potential of eq. 3. Note that $\alpha$ varies inversely with
interaction range. The figure shows that a larger interaction range indeed
leads to improved performance by the Correlation Tracking method and worse
performance by the Crocker-Grier algorithm. When the range of interaction is
sufficiently small, the accuracies of the Correlation Tracking algorithm and
the Crocker-Grier algorithm converge.

Figure 2: Comparison of the Correlation Tracking algorithm to the Crocker-
Grier algorithm as a function of interaction range. The results are computed
on simulated data using a Morse potential (eq. 3). The x-axis shows the
parameter $\alpha$ in the Morse potential which varies inversely with
interaction range. The y-axis shows the fraction of correctly identified
particle labels. Error bars are based on 25 independent simulations.

##  IV Results on Experimental Data

To verify the behavior on experimental data, we tested our tracking algorithm
on experimental hard sphere colloids. In this system, since there are no long-
range interactions and the motion is largely uncorrelated, we expect that the
algorithm should match the performance of the Crocker-Grier algorithm.

The experimental system consists of core-shell colloidal particles undergoing
Brownian motion in a density-matched and refractive-index-matched solvent. The
core and shell of the particles have the same composition: poly trifluoroethyl
methacrylate (pTFEMA) and poly tert-butyl methacrylate (pTBMA). The cores are
fluorescently dyed. The core diameter is 0.8 $\mu$m, and the particle diameter
is 1.75 $\mu$m. This core-shell design separates fluorescent spots in confocal
imaging and ensures accurate particle locating in a dense colloidal
suspension. A polyelectrolyte brush is grown on the particle surface to
prevent aggregation. The particles are suspended in a mixture of urea and
formamide, with tetrabutylammonium bromide (TBAB) added as salt to screen
charges and to ensure that the particle interaction closely approximates that
of hard spheres. This colloidal solution is sealed in a glass chamber with
dimensions 20 mm x 10 mm x 0.4 mm. The particles take up 43% of the total
volume.

Figure 3: 3D rendered image of the colloidal particles in the first frame of
the tracking experiment. The front half of the sample has been removed to
reveal the red particles at the center, which are used for the tracking
analysis. The particle diameter is 1.75 $\mu$m.

We use fluorescent confocal microscopy to image the particles. Fig. 3 shows
part of the imaged volume, which measures 209 $\mu$m x 118 $\mu$m x 10 $\mu$m
and contains about 38000 particles. The imaged volume is far from the sample
boundaries to prevent edge effects. The particle centers are located and
tracked using TrackPy[6], a python package based on the Crocker-Grier
algorithm. We chose 55 particles at the center of the image for analysis,
marked in red in Fig. 3. This ensures that the analyzed particles are present
through all 100 frames. The frame rate is 1.25 seconds per 3D stack. This high
frame rate allows accurate tracking by the Crocker-Grier algorithm. We compare
frames at different time intervals to test changes in tracking accuracy.

Figure 4: Comparison of Correlation Tracking to Crocker-Grier in a three-
dimensional experimental hard sphere system. The x-axis shows the mean
particle displacement in units of the particle diameter. The y-axis shows the
fraction of correctly labeled particles. The error bars show the standard
error, computed for 10 independent tracking events that use different central
particles and different starting frames.

We approximate the hard sphere behavior in the simulations with a WCA
potential [8], which consists of the repulsive portion of a Lennard-Jones
potential. We use this potential instead of an exact hard sphere interaction
to simplify the calculation of the dynamics. We first simulate an ensemble of
trajectories starting from the positions in the first frame of the
experimental data. This ensemble is used to parameterize a probability
distribution over particle positions at a later time. We then label the
particles at this later time using the probability distribution we have
computed from the simulations. Fig 4 shows the results.

##  V Results on colored particles

Lastly, we consider systems with multiple particle types, such as colloids of
different colors. Knowledge of these types can be used to inform the tracking
algorithm and offer more accurate predictions. First, consider the case of $N$
particles interacting by a Morse potential. Suppose each particle is dyed one
of $M$ possible colors. If $M=1$, we recover the results found above. As $M$
approaches $N$, all the particles become distinguishable, and it is always
possible to track every particle accurately. In the intermediate range, more
colors lead to an effectively reduced density without interfering with
correlated particle motion. Thus, noting that introducing $N$ independent
colors is rarely experimentally feasible, we examine the tradeoff between
limiting the number of colors and the accuracy of tracking using both the
Correlated Tracking algorithm and the Crocker-Grier algorithm in Fig. 5.

To simulate adding independent colors, we added the concept of species to our
simulation. We assume that if a particle is dyed a given color, it may only be
mapped to a particle of the same color in the next frame. To enforce this
constraint, we treat each color as a different species and track it
separately, which prevents errors in labeling across colors in both the
Correlation and Crocker-Grier tracking algorithms. In this example, the
species argument is only relevant to the tracking analysis, and not to the
simulations themselves; here we assume that differently colored particles
interact by a Morse potential in the same way as particles with the same
color.

Figure 5: Comparison of the Correlation Tracking algorithm to the Crocker-
Grier algorithm in simulated data on particles with multiple colors,
interacting via a Morse potential (left) and a soft sphere potential (right).
The particles were allowed to evolve such that their distance, measured in
units of $r_{0}$ per average inter-particle separation, was 0.11 for the Morse
potential, and 0.069 for the soft sphere potential. The y-axis shows the
fraction of correctly labeled particles in that frame. Error bars are based on
25 independent simulations.

Fig. 5 shows that increasing the number of colors increases tracking accuracy
for both algorithms. Adding too many colors leads to diminishing returns in
tracking accuracy. For the Morse potential, Crocker-Grier achieves very high
accuracies with 4 colors. Correlation Tracking improves this and achieves
nearly perfect tracking accuracy with only 3 colors. For the soft sphere
potential, Crocker-Grier achieves high accuracy with 3 colors. The Correlation
Tracking algorithm matches the same level of accuracy Crocker-Grier achieves
with only 2 colors.

##  VI Discussion

Our results on experimental data indicate that the Correlation Tracking
algorithm is robust. Even though we used an imperfect description of the hard
sphere potential in the experimental system, the algorithm recovers Crocker-
Grier’s performance in the diffusive regime. Many algorithms, particularly
machine learning algorithms, require extensive modeling of noise sources in
the experimental system to achieve accurate results on experimental data [9,
10, 11]. For instance, computer vision algorithms often require modeling noise
due to the type of camera and the method of transmitting data [9]. In
contrast, we are able to accurately track experimental data without
introducing any noise modeling. Therefore, our model is robust to noise
sources in the experimental data, such as polydispersity and positional
uncertainty. Furthermore, the Correlation Tracking method relies on an
accurate description of the potential. However, we noted that our description
of the potential is imperfect, yet we are still able to track particles
accurately. As long as the potential is accurate over the range and timescale
of the measurements, Correlation Tracking will be accurate.

The increase in tracking accuracy due to the Correlation Tracking algorithm
comes at a significant computational cost. We mitigate that cost by taking
advantage of the automatic vectorization and compilation offered by the
package JAX [12, 13], so that simulations are run simultaneously. Tracking a
set of 5 particles across one pair of frames with a batch size of 10,000
simulations used to estimate the correlation matrix takes approximately 4
minutes on a K80 GPU on Google Colab, whereas performing the same analysis
with the Crocker-Grier method takes only a few seconds. However, the proof of
principle of the Correlation Tracking algorithm leaves a lot of room for
improvements in computational efficiency. Moreover, for commonly used
implementations of the Crocker-Grier algorithm, such as TrackPy, the
computational cost grows factorially when the number of particles in a
tracking set increases. This makes tracking in a dense particle system
particularly costly in computation, because many candidate particles need to
be considered in a tracking set.

Though our method increases computation time, it may reduce the overall
experimental effort. As discussed in section III, one method of increasing
tracking accuracy is to introduce multiple colors for particles in an
experimental system. Creating each color, however, comes at a high cost in
experiment time. Since the Correlation Tracking method significantly reduces
the number of colors needed to reach the highest levels of tracking accuracy,
it may save significant time.

The Correlation Tracking algorithm uses simulation in a fundamentally new way.
To date, simulations have been used to discover features and mechanisms of
material and biological systems that cannot yet be observed in experiment[14,
15], to deepen our understanding of experimental results [16, 17], or to guide
future directions for experimental exploration [18]. This algorithm, however,
intertwines simulation and experiment in the data analysis procedure. The
simulated data is used to directly increase the accuracy of particle tracking
by rooting the data analysis method in the physics of the system. There is a
rich space to explore by bringing simulation and experiment in direct contact,
a space that we anticipate will contain new lines of research on both sides.

As the Correlation Tracking algorithm increases tracking accuracy, it allows
accurate tracking for longer image acquisition intervals in experiments. This
enables a taller imaging region in 3D confocal imaging. For example, in the
experiments of this paper, the microscope can reach 40 frames per second while
maintaining the necessary image resolution in a xy-plane for accurate particle
locating. The step size in the z-direction needs to be 0.18 $\mu$m for
accurate particle locating in the z dimension. The dynamics of the colloidal
system limit the 3D imaging interval to less than 1.5 seconds to ensure
accurate particle tracking using the Crocker-Grier algorithm. These
restrictions give a maximum of 60 xy-planes or 10.8 $\mu$m at each time step.
With this limitation, only 7 particle layers can be scanned for the
1.75-$\mu$m-sized particles used in this experiment. A significant portion of
the particles moves out of the imaging region in a longer experiment. This
makes it difficult to track each particle over a long time and to study
kinetic phenomena. However, with Correlation Tracking, more particle layers
can be imaged at each time step, and long-time particle tracking becomes more
feasible.

The Correlation Tracking method promises to be of greatest use in studies with
highly correlated motion, such as studies of crystal nucleation or shear
thickening in attractive colloidal systems. The algorithm can be applied to
any experimental system for which the particle interactions are well-
characterized. Moreover, recent work has had significant success in
approximating interaction potentials from simulated data using machine
learning algorithms [19, 20]. If these interaction potentials were learned
from experimental data, the Correlation Tracking algorithm could be used to
provide accurate tracking. These approximations need only be accurate over
short distance and time scales to be useful for the Correlation Tracking
algorithm. Further, the Correlation Tracking algorithm could be inverted to
provide estimates of interaction potentials in experimental systems: given
data with accurate particle tracks, the Correlation Tracking algorithm would
provide a probability that the tracks are well explained by a given
interaction potential. One could converge on an accurate potential by
proposing a given interaction potential, evaluating the probability that it
explains the experimental particle tracks, and accordingly updating the
proposed potentials, for example with gradient-based optimization methods.

## Acknowledgements

We thank Peter Lu for many valuable discussions and experimental guidance.
This material is based on work supported by the National Science Foundation
Graduate Research Fellowship under Grant No. DGE1745303, the Harvard Materials
Research Science and Engineering Center (DMR 20-11754), the Office of Naval
Research (ONR N00014-17-1-3029) and the Simons Foundation.

## References

  * Rossi  _et al._ [2011] L. Rossi, S. Sacanna, W. T. Irvine, P. M. Chaikin, D. J. Pine, and A. P. Philipse, Cubic crystals from cubic colloids, Soft Matter 7, 4139 (2011).
  * Wang  _et al._ [2012] Z. Wang, F. Wang, Y. Peng, Z. Zheng, and Y. Han, Imaging the homogeneous nucleation during the melting of superheated colloidal crystals, Science 338, 87 (2012).
  * Schall  _et al._ [2007] P. Schall, D. A. Weitz, and F. Spaepen, Structural rearrangements that govern flow in colloidal glasses, Science 318, 1895 (2007).
  * Crocker and Grier [1996] J. C. Crocker and D. G. Grier, Methods of digital video microscopy for colloidal studies, Journal of colloid and interface science 179, 298 (1996).
  * Gao and Kilfoil [2009] Y. Gao and M. L. Kilfoil, Accurate detection and complete tracking of large populations of features in three dimensions, Optics express 17, 4685 (2009).
  * Allan  _et al._ [2021] D. B. Allan, T. Caswell, N. C. Keim, C. M. van der Wel, and R. W. Verweij, soft-matter/trackpy: Trackpy v0.5.0 (2021).
  * Goodrich  _et al._ [2021] C. P. Goodrich, E. M. King, S. S. Schoenholz, E. D. Cubuk, and M. P. Brenner, Designing self-assembling kinetics with differentiable statistical physics models, Proceedings of the National Academy of Sciences 118, 10.1073/pnas.2024083118 (2021).
  * Andersen  _et al._ [1971] H. C. Andersen, J. D. Weeks, and D. Chandler, Relationship between the hard-sphere fluid and fluids with realistic repulsive forces, Physical Review A 4, 1597 (1971).
  * Yue  _et al._ [2019] Z. Yue, H. Yong, Q. Zhao, L. Zhang, and D. Meng, Variational denoising network: Toward blind noise modeling and removal, arXiv preprint arXiv:1908.11314 (2019).
  * Wu  _et al._ [2018] J. Wu, S. Ruan, C. Lian, S. Mutic, M. A. Anastasio, and H. Li, Active learning with noise modeling for medical image annotation, in  _2018 IEEE 15th International Symposium on Biomedical Imaging (ISBI 2018)_ (2018) pp. 298–301.
  * Gupta and Gupta [2019] S. Gupta and A. Gupta, Dealing with noise problem in machine learning data-sets: A systematic review, Procedia Computer Science 161, 466 (2019).
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
  * Schoenholz and Cubuk [2020] S. S. Schoenholz and E. D. Cubuk, Jax m.d. a framework for differentiable physics, in  _Advances in Neural Information Processing Systems_, Vol. 33 (Curran Associates, Inc., 2020).
  * Porter  _et al._ [2021] C. L. Porter, S. L. Diamond, T. Sinno, and J. C. Crocker, Shear-driven rolling of dna-adhesive microspheres, Biophysical Journal 120, 2102 (2021).
  * Du  _et al._ [2017] C. X. Du, G. van Anders, R. S. Newman, and S. C. Glotzer, Shape-driven solid–solid transitions in colloids, Proceedings of the National Academy of Sciences 114, E3892 (2017).
  * Li  _et al._ [2021] L. Li, C. Goodrich, H. Yang, K. R. Phillips, Z. Jia, H. Chen, L. Wang, J. Zhong, A. Liu, J. Lu, _et al._ , Microscopic origins of the crystallographically preferred growth in evaporation-induced colloidal crystals, Proceedings of the National Academy of Sciences 118 (2021).
  * Audus  _et al._ [2021] D. J. Audus, S. Ali, A. M. Rumyantsev, Y. Ma, J. J. de Pablo, and V. M. Prabhu, Molecular mass dependence of interfacial tension in complex coacervation, Physical Review Letters 126, 237801 (2021).
  * Kimchi  _et al._ [2020] O. Kimchi, C. P. Goodrich, A. Courbet, A. I. Curatolo, N. B. Woodall, D. Baker, and M. P. Brenner, Self-assembly–based posttranslational protein oscillators, Science Advances 6, eabc1939 (2020).
  * Vandermause  _et al._ [2021] J. Vandermause, Y. Xie, J. S. Lim, C. J. Owen, and B. Kozinsky, Active learning of reactive bayesian force fields: Application to heterogeneous hydrogen-platinum catalysis dynamics, arXiv preprint arXiv:2106.01949 (2021).
  * Chan  _et al._ [2019] H. Chan, B. Narayanan, M. J. Cherukara, F. G. Sen, K. Sasikumar, S. K. Gray, M. K. Chan, and S. K. Sankaranarayanan, Machine learning classical interatomic potentials for molecular dynamics from first-principles training data, The Journal of Physical Chemistry C 123, 6941 (2019).


# SynGraphy: Succinct Summarisation of Large Networks via Small Synthetic
Representative Graphs

Jérôme Kunegis1, Pawan Kumar2,
Jun Sun3, Anna Samoilenko2, and Giuseppe Pirró4
1 University of Namur
2 University of Koblenz–Landau
3 University of Stuttgart
4 Institute for High Performance Computing and Networking

###### Abstract

We describe SynGraphy, a method for visually summarising the structure of
large network datasets that works by drawing smaller graphs generated to have
similar structural properties to the input graphs. Visualising complex
networks is crucial to understand and make sense of networked data and the
relationships it represents. Due to the large size of many networks,
visualisation is extremely difficult; the simple method of _drawing_ large
networks like those of Facebook or Twitter leads to graphics that convey
little or no information. While modern graph layout algorithms can scale
computationally to large networks, their output tends to a common _hairball_
look, which makes it difficult to even distinguish different graphs from each
other. Graph sampling and graph coarsening techniques partially address these
limitations but they are only able to preserve a subset of the properties of
the original graphs. In this paper we take the problem of visualising large
graphs from a novel perspective: we leave the original graph’s nodes and edges
behind, and instead summarise its properties such as the clustering
coefficient and bipartivity by generating a completely new graph whose
structural properties match that of the original graph. To verify the utility
of this approach as compared to other graph visualisation algorithms, we
perform an experimental evaluation in which we repeatedly asked experimental
subjects (professionals in graph mining and related areas) to determine which
of two given graphs has a given structural property and then assess which
visualisation algorithm helped in identifying the correct answer. Our
summarisation approach SynGraphy compares favourably to other techniques on a
variety of networks.

##  1 Introduction

Network data is nowadays ubiquitous; social networks, protein networks and the
Web are just a few examples of its spread. If network analysis is a crucial
challenge to gain insights from such wealth of data, network visualisation
represents an equally important aspect, especially because of its immediate
usefulness to both network analysts and non-expert users; displaying (portions
of) a network or structural properties (e.g. the degree distribution) makes
more accessible the specific task at hand. The visualisation of small graphs
poses no major problem; it can be achieved by graph drawing, i.e., placing
each node of a graph at a specific position on the plane, and drawing edges as
lines. The task then becomes to find suitable coordinates of nodes such that
the drawing has judicious properties for the user (Purchase et al, 1995),
including the fact that connected nodes are near to each other, lines should
not cross, points and lines should be uniformly distributed over the drawing
area, and symmetries should be highlighted.

\begin{overpic}[width=86.72267pt]{hairball-gallery/Rice31-hairball2-1.png}
\put(5.0,10.0){{{A}}} \end{overpic}

\begin{overpic}[width=86.72267pt]{hairball-gallery/IMAGE-4.png}
\put(5.0,10.0){{{B}}} \end{overpic}

\begin{overpic}[width=86.72267pt]{syngraphy.a.enron.clean}
\put(5.0,10.0){{{C}}}
\end{overpic}\begin{overpic}[width=86.72267pt]{img/syngraphy.a.prosper-loans}
\put(5.0,10.0){{{D}}} \end{overpic}

Figure 1:  Comparing usual _hairball_ -like graph drawings (A–B) with the
SynGraphy graph visualisation method introduced in this paper (C–D). (A) a
social network (Nocaj, 2015), (B) a block chain network (Harrigan, 2014), (C)
an email network (Klimt and Yang, 2004), (D) a loan network (Redmond and
Cunningham, 2013). The goal of this paper is to avoid the “hairball” effect
and make graph-theoretical features of the network apparent.

In data mining applications such as those arising on the web however, networks
are so large that graph drawing is highly problematic, for at least two
aspects. First, the computational cost associated to a certain graph drawing
method which models properties like those listed above is usually high.
Second, even if one comes up with a suitable algorithm, the drawing of very
large graphs inherently leads to the “hairball effect”, i.e., the fact that
all graphs look the same when drawn – see the article by Schulz and Hurter
(2013) for a comprehensive overview about the topic. Examples of hairball-like
plots are shown in Figure 1: even if the networks come from different domains,
it is difficult to discern any of their structural properties.

To avoid the hairball effect, other methods are typically used, such as
drawing only a part of the network (graph sampling), or aggregating its nodes
into groups to show its main overall structural features (graph coarsening).
These methods however lead to the problem of representativity, i.e., the
sampled or coarsened graphs have structural properties differing from the
original graphs. Some methods have been developed to bias graph sampling and
coarsening algorithms to preserve certain graph properties (Cem et al, 2013).
However, these approaches are not focused on taking into account specific
user’s information needs; they rely on the layout itself as a means to satisfy
all information needs.

Contributions. In this paper we investigate the problem of visualising large
graphs from the perspective of graph summarisation. Instead of applying
sampling or coarsening, _we devise an algorithm that generates synthetic
graphs whose properties match that of the original graphs_. The crucial aspect
is that the syntactic graphs generated with our algorithm do not preserve any
nodes or group of nodes from the original graphs; instead, they preserve the
original structural and graph-theoretical properties. This allows us to reduce
the problem of visualising large networks to the problem of drawing smaller,
representative graphs instead of the original network. For a given input
network, our algorithm SynGraphy performs the following steps: (i) it first
measures several (and dynamically selectable) statistical properties of a
network; (ii) scales them down to a given, small predefined number of nodes;
(iii) generates a synthetic network with the required scaled-down statistical
properties; and (iv) produces a graph drawing of the synthetic network using a
classical graph drawing algorithm. The scaling down of statistics is trained
on a large set of network datasets, in order to faithfully distinguish
particularly low and high values of each network statistic. The summarisations
of networks generated by our algorithm are evaluated in an experiment with
human volunteers that are professionals in network analysis and related
fields, asking them to estimate which of two given graphs has a higher or
lower value of a given statistical property (e.g., clustering coefficient).
SynGraphy performed well in comparison to plain graph drawing algorithms and
graph sampling algorithms.

The remainder of the paper is structured as follows. Section 2 reviews the
areas of graph drawing, sampling and network characteristics. Section 3
describes the proposed algorithm. The evaluation of our approach is discussed
in Section 4. We conclude and sketch future work in Section 5.

##  2 Related Work

The problem of graph visualisation is of primary importance in multiple fields
of research, as many types of data can be modelled as networks. In order to
visualise a graph, the most straightforward – and most common – method is to
use graph drawing, i.e., producing a plot in which each node of the graph is
represented as a point, and each edge as a line connecting a pair of points.
Other, more advanced methods may try to show only a subset of a graph, or to
summarise its properties in a specific way. The visualisation of network
datasets can serve multiple purposes that can be clearly distinguished, and
which lead to different visualisation methods:

  * •

Exploration allows to explore nodes and edges of a graph. In the simplest
case, a graph drawing allows to follow individual nodes and edges, as long as
the graph is not too large. If the graph is large, a system of navigation may
be used, in which one may click, for instance, on nodes to explore their
neighbourhood;

  * •

Structural properties – graph visualisation should allow to recognise
structural properties of the graph drawn as a whole, e.g., whether the graph
has a high clustering coefficient or not (van Ham et al, 2009). This allows to
perform other tasks such as comparing two networks, or finding an outlier in a
given set of networks.

In this work the focus is on recognising structural properties: the goal is to
allow both users and network analysts to recognise such properties in the
graph visualisation, thus leading to the task of graph summarisation – the
exploration of nodes and edges of the graph is then an orthogonal goal to
ours.

###  2.1 Graph Drawing

Graph drawing, as opposed to graph visualisation, specifically refers to the
task of laying out nodes and edges of a graph on a planar surface in order to
produce a figure that shows the input graph as a whole. Graph drawing has also
been called graph layouting and graph embedding. In general, graph drawing
algorithms are based on optimising a certain cost or utility function, which
models the quality of the graph drawing. Examples of features that can be
incorporated are: drawing nodes near to their neighbours, avoiding crossing of
edges, and avoiding multiple nodes at the same coordinates. A subset of such
algorithms performs a force-directed layout, i.e., it simulates a physical
force between nodes based on the physical model of springs connected the
nodes. In very simple cases, these can be formulated as the eigenvectors of
certain graph-characteristic matrices, in particular of the Laplacian matrix
(Koren, 2003). A well-known graph layout algorithm that is widespread in
practice is the algorithm of Fruchterman and Reingold (1991), which bases its
physical model on placing nodes near to its neighbours, but also on not
placing nodes too near to each other – this makes the problem nonlinear and
thus it cannot be expressed as eigenvectors of a matrix, but instead can be
solved efficiently by direct optimisation algorithms, even if a global minimum
of the energy function cannot always be found.

The obvious problem with pure graph drawing methods is that they don’t scale
visually – if the number of nodes is too high, the resulting drawing will
contain so many dots and lines that it will look like a single big patch of
dots, without any recognisable network structure. This is the type of drawing
that is usually called a “hairball”, and which is illustrated in Figure 1
using examples from the literature.

Another class of graph drawing algorithm approaches the problem from the point
of view of wanting to draw the network such that certain particular properties
of the network are highlighted, and thus choose to optimise the layout such
that the given properties are made visible. An example from this class of
algorithms is the hive plot (Krzywinski et al, 2012), in which nodes of the
given graph are assigned to one of three or more axes, which may be divided
into segments. Nodes are ordered on a segment based on properties such as
connectivity, density, centrality or quantitative annotations. The user then
is free to choose whatever rules fit their data and visualisation
requirements. Edges are drawn as curves. As we will see, the SynGraphy method
presented in this paper has in common with hive plots that it makes it
possible to assess network structure because both are based on network
properties, not on aesthetic layout, making visualisations of two networks
directly comparable.

###  2.2 Graph Summarisation

The inadequacy of drawing a complete graph when its size is too large leads to
several methods in which the size of the input graph is reduced to obtain a
smaller graph, which can then be drawn using any of the methods mentioned in
the previous section. Such methods can be roughly divided into two types: (i)
graph sampling techniques which choose a subset of the given graph’s nodes and
edges, and (ii) graph coarsening techniques that aggregate the graph’s nodes
into new nodes. These methods can be interpreted as graph summarisation:
Automatic summarisation refers to the generation of shorter documents that
summarise a given document, reducing it to the most important parts. Although
the term _summarisation_ has recently been applied to graphs, for instance by
Liu and colleagues (2017), the term is, as of 2020, not common in the field.

Graph sampling. Graph sampling or graph sparsification refers to reducing a
graph to a subset of its nodes and edges, with the goal of preserving certain
graph properties. Simple graph sampling techniques perform a random sampling
of nodes and edges, where each node or each edge is kept with a given
probability (Hu and Lau, 2013). In such methods, the expected properties of
the resulting network can sometimes be derived in closed form. As an example,
an edge sampling procedure that keeps each edge with probability $p$ will lead
to a graph with clustering coefficient $pc$, where $c$ is the original graph’s
clustering coefficient. On the other hand, certain graph sampling methods do
not choose nodes and edges randomly, but use heuristics. For instance, an
approach to filter a graph by removing edges in order of increasing
betweenness centrality is described by Jia and colleagues (2008). The method
preserves certain properties such as the connectivity, but there are still
some features in the resulting graph that are hidden and cannot be identified
for the entire graph.

Graph coarsening. Graph coarsening consists in aggregating multiple nodes of
the input graph into a single node of the output graph, essentially performing
vertex identification or edge contraction (Shawn et al, 2011). As a result,
high-level structural features such as connections between large communities
are preserved, while information at the node level such as the clustering
coefficient is usually lost. More examples (and drawings) of graph coarsening
methods are given by Nocaj and colleagues (2015).

###  2.3 Graph Characteristics

In order to generate a graph with similar characteristics to an input graph,
the characteristics of the latter need to be numerically quantified. This is
achieved through the use of _graph statistics_ , i.e., numerical measures of a
graph. Graph statistics comprise the important class of subgraph counts, which
count the number of occurrences of a particular small graph as a subgraph. For
instance, the number of triangles is an often used graph statistic. The number
of nodes and edges of a graph are subgraph counts too, corresponding to the
subgraph consisting of a single node and a single edge respectively. Other
graph statistics may be derived from subgraph counts, such as the clustering
coefficient, which can be expressed as $3t/s$, where $t$ is the number of
triangles and $s$ is the number of wedges, while other graph statistics are
wholly unrelated to subgraph counts (e.g., the diameter). The graph statistics
useful for the purposes of our algorithm, and described below, are summarised
in Table 1. The input graph, as well as all other graphs used in this paper
are undirected, without parallel edges, and loopless.

Table 1:  The notation used in this paper. As a general rule, symbols without an apostrophe denote the potentially large input graph, while symbols with an apostrophe denote the small output graph to be drawn.  Notation | Meaning
---|---
$G=(V,E)$ | The input graph (potentially large)
$G^{\prime}=(V^{\prime},E^{\prime})$ | The output graph (to be drawn, small)
$n$ | Number of nodes
$m$ | Number of edges
$s$ | Number of wedges
$z$ | Number of claws
$x$ | Number of crosses
$t$ | Number of triangles
$q$ | Number of squares
$d$ | Average degree
$c$ | Clustering coefficient
$y$ | 4-Clustering coefficient
$b$ | Bipartivity
$\delta$ | Diameter
$\rho$ | Degree assortativity
$\square^{\prime}$ | Scaled-down value of any statistic $\square$

Network size and volume ($n,m$). The number of nodes and edges in a network.

Number of wedges ($s$). A wedge is a subgraph consisting of a central node
connected to two other nodes. Equivalently, wedges are 2-stars or 2-paths, and
their count can be expressed as

$s=\sum_{u\in V}\frac{d(u)(d(u)-1)}{2},$ (1)

where $d(u)$ is the degree of node $u$.

Number of claws ($z$). A claw is a subgraph equivalent to a 3-star.

Number of triangles ($t$). A triangle is a subgraph consisting of three
mutually connected nodes.

Number of squares ($q$). A square is a group of four nodes connected to each
other in the shape of a square, i.e., a cycle of length four.

Clustering coefficient ($c$). The clustering coefficient represents the
presence of a heightened number of triangles in the network (Newman, 2003b),
and can be expressed as

$c=\frac{3t}{s}.$ (2)

The clustering coefficient can take a value between zero and one, where the
value one denotes that all possible triangles are present, and zero denotes a
triangle-free network.

4-Clustering coefficient ($y$). The 4-clustering coefficient $y$ is analogous
to the clustering coefficient, but measures the probability that a randomly
chosen path of length three is closed. It equals

$y=\frac{4q}{P_{3}},$ (3)

where $P_{3}$ is the number of paths of length three. The 4-clustering
coefficient complements the clustering coefficient, and may be used in
bipartite or near-bipartite networks (Latapy et al, 2008).

Bipartivity ($b$). We use as a measure of bipartivity the ratio of the
smallest and largest eigenvalue of the graph’s adjacency matrix $\mathbf{A}$
(Kunegis, 2015):

$b=\left|\frac{\lambda_{\min}[\mathbf{A}]}{\lambda_{\max}[\mathbf{A}]}\right|$
(4)

Average degree ($d$). The average degree measures how many edges are connected
to a node on average. The average degree is sometimes called the density. It
equals $d=2m/n$.

Diameter ($\delta$). The diameter of a network is the longest shortest
distance between any two nodes. A distance is measured by the number of edges
needed to reach one node from another.

Degree assortativity ($\rho$). In some networks, the nodes with high degree
are connected to other nodes with high degree, while nodes with low degree are
connected to nodes with low degree. This is measured by the Pearson
correlation coefficient, taken over all pairs of connected nodes (Newman,
2003a).

##  3 The SynGraphy Approach

We are now ready to introduce our algorithm. Let $G=(V,E)$ be the large input
graph that is to be summarised. Our algorithm proceeds in four steps:

  1. (i)

Compute the subgraph counts of $G$

  2. (ii)

Scale down the subgraph counts

  3. (iii)

Generate $G^{\prime}$ with the scaled down subgraph counts

  4. (iv)

Generate a drawing of $G^{\prime}$

Steps (i–iii) are described in the following subsections, while Step (iv) uses
the algorithm by Fruchterman and Reingold (1991). For Step (ii), we describe
two variants.

###  3.1 Scaling Down of Graph Statistics

The core idea of our algorithm to summarise the input graph $G$ is to generate
a new graph $G^{\prime}$ with analogous structural properties to those of $G$.
Hence, the focus is on providing a meaningful rendering of $G$’s structural
properties instead of $G$’s aesthetic. The size of $G^{\prime}$, i.e., the
number of nodes $n^{\prime}$, can be chosen arbitrarily. For the sake of
presentation, in what follows we set it to the value $n^{\prime}=80$. As for
the other statistics, their new values in the output graph $G^{\prime}$ need
to be computed before $G^{\prime}$ can be generated. Certain graph statistics,
such as the average degree $d$ and clustering coefficient $c$ are in principle
scale-invariant, and could be simply carried over from $G$ without
modification, i.e., it is possible to use $d^{\prime}=d$ and $c^{\prime}=c$.
Other statistics such as the number of triangles clearly scale with $n$ and
their values have to be adjusted. In fact, the values of $c$ and $d$ may have
to be adjusted too, as illustrated in Figure 2. In this light, we investigate
two methods for determining the new values of these graph statistics:

  * •

(NO): normalising all statistics empirically;

  * •

(SI): use of statistics that by construction are independent of size.

Empirical Normal Distribution (NO). One possible assumption is that when
normalised in certain ways, certain network statistics are size-independent.
An example for this is the clustering coefficient. However, as we can derive
empirically from actual network datasets, this is not entirely correct: the
clustering coefficient does correlate with the size of the network, when
measured in real-world networks, as shown in Figure 2. This observation leads
to the following method for scaling down graph statistics: fit a distribution
to the statistics of a set of known networks, and then read out the most
likely values of the requested graph statistics, given the target size
$n^{\prime}$.

Figure 2:  The size $n$ versus the clustering coefficient $c$ in real-world
networks. This plot shows that the clustering coefficient is not independent
of size, and gives an approximate relationship of $c\sim n^{-0.5}$. The
networks are used as provided by the KONECT project (Kunegis, 2013).

Given a set of real-world networks and a set of statistics, we compute the
values of the statistics for all networks, and then fit a multivariate normal
distribution to the data. The multivariate normal density is a generalisation
of the univariate normal density to $k\geq 2$ dimensions (Johnson and Wichern,
2007). The $k$-dimensional normal density of a $k$-dimensional variable
$\mathbf{x}$ takes the form:

$f(\mathbf{x})=\frac{1}{\sqrt{(2\pi)^{k}|\mathbf{\Sigma}|}}\exp\left\\{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu})^{\mathrm{T}}\mathbf{\Sigma}^{-1}(\mathbf{x}-\mathbf{\mu})\right\\},$
(5)

where the $k$-vector $\mathbf{\mu}$ is the mean of the distribution, the
$k\times k$ matrix $\mathbf{\Sigma}$ the covariance matrix, and
$|\mathbf{\Sigma}|$ its determinant. The conditional probability distribution
under a fixed $n=n^{\prime}$ is itself a multivariate normal distribution
whose mean is given by

$\mathbf{\mu}^{\prime}=\mathbf{\mu}_{*}-\frac{\mathbf{\mu}_{n}-n^{\prime}}{\mathbf{\Sigma}_{n,n}}\mathbf{\Sigma}_{*,n}^{\phantom{-1}},$
(6)

where the index $\square_{n}$ denotes a restriction of the given quantity to
the value corresponding to the size statistic $n$, and $\square_{*}$ denotes a
restriction to all other columns and rows. Thus, for given statistic $\alpha$,
the new value for the smaller graph is given by

$\alpha^{\prime}=\alpha-(\mu_{n}-n^{\prime})\frac{\mathbf{\Sigma}_{\alpha,n}}{\mathbf{\Sigma}_{n,n}}$
(7)

Size-Independent Statistics (SI). Under the assumption that most subgraph
counts can be expected to have values growing with the number of nodes in the
network, the number of edges $m$ for instance can be written as

$m=\frac{d}{2}n,$ (8)

where $d$ is the average degree in the network. In a social network for
instance, $d$ represents the average number of friends, and we can make the
reasonable assumption that $d$ is independent from the network size, i.e.,
people don’t have more friends because the world is larger. Note that this
assumption is not necessarily true, and the opposite case is taken into
account in the previous paragraph. Assuming that $d$ is constant, it follows
that we can use $d$ as a size-independent statistic, i.e., we can compute $d$
for the input graph, and then use $d^{\prime}=d$.

Similarly, other subgraph counts can be derived to grow as powers of the
network’s size, given an underlying graph model. In fact, this is true for
whole classes of network models; in the Erdős–Rényi model for instance, the
expected value of any subgraph count is a polynomial of the graph size
(Ruciński, 1988). For instance, the number of triangles $t$ in Erdős–Rényi
graphs grows as the cube of the number of vertices.

For each subgraph count $\alpha$, we perform the following procedure: We
derive the expected value of $\alpha$ under suitable conditions, calling it
$\hat{\alpha}$; and then derive from it the normalised statistic
$\alpha_{\mathrm{n}}^{\phantom{x}}=\alpha/\hat{\alpha}$. Finally, we use
$\alpha_{\mathrm{n}}^{\prime}=\alpha_{\mathrm{n}}^{\phantom{x}}$, giving the
new value of the plain subgraph count:

$\alpha^{\prime}=\frac{\hat{\alpha}^{\prime}}{\hat{\alpha}\phantom{{}^{\prime}}}\alpha$
(9)

We now apply this method to individual subgraph count statistics. For each
subgraph count statistic $\alpha$, we derive the expected value $\hat{\alpha}$
and the scaled down subgraph count statistic $\alpha^{\prime}$.

Number of edges ($m$): Given a fixed average degree $d^{\prime}=d=2m/n$, we
get

$\displaystyle\hat{m}=\frac{d}{2}n,\qquad
m^{\prime}=\frac{n^{\prime}}{n\phantom{{}^{\prime}}}m.$ (10)

Number of wedges ($s$): Given an average degree $d$ and an equal distribution
of degrees, the expected number of wedges can be computed by noting that the
number of wedges $s(u)$ centred at each node $u$ is given by ${d(u)\choose
2}$.

$\displaystyle\hat{s}=n{d\choose 2},\qquad
s^{\prime}=\frac{n^{\prime}}{n\phantom{{}^{\prime}}}s$ (11)

Number of claws and crosses ($z,x$): Claws and crosses are 3-stars and 4-stars
and as such, their expected values are analogous to those for 2-stars:

$\displaystyle\hat{z}=n{d\choose 3},\qquad\hat{x}=n{d\choose 4},\qquad
z^{\prime}=\frac{n^{\prime}}{n\phantom{{}^{\prime}}}z,\qquad
x^{\prime}=\frac{n^{\prime}}{n\phantom{{}^{\prime}}}x$ (12)

Number of triangles ($t$): For the number of triangles, we assume a constant
clustering coefficient $c$. The clustering coefficient can be expressed by
definition as $c=3t/s$. This gives

$\displaystyle\hat{t}=\frac{c}{3}\hat{s}=\frac{c}{3}{d\choose 2}n,\qquad
t^{\prime}=\frac{n^{\prime}}{n\phantom{{}^{\prime}}}t.$ (13)

Number of squares ($q$): This case is defined analogously to the number of
triangles.

$\displaystyle\hat{q}=\frac{y}{4}m(d-1)^{2}=\frac{y}{4}\frac{d}{2}(d-1)^{2}n,\qquad
q^{\prime}=\frac{n^{\prime}}{n\phantom{{}^{\prime}}}q$ (14)

###  3.2 Graph Generation

As previously mentioned, graphs are generated by taking the fixed size of the
network ($n^{\prime}=80$), and corresponding scaled down values of six
statistics including the count of edges, wedges, claws, crosses, triangles,
and squares as obtained by applying the two methods (i.e., NO and SI)
discussed above. The layout is drawn using the force-directed algorithm of
Fruchterman and Reingold (1991).

The general form of our graph generation algorithm applying to any set of
network statistics is given in Algorithm 1. The algorithm works by taking as a
starting point an Erdős–Rényi graph with the correct number of nodes and
edges, and then modifying the graph iteratively until the resulting graph is
near enough to the target. At each step of the iteration, we need to consider
a certain number of possible changes in the graph, and choose the change which
leads to the lowest error measure. In order to compute the changes in the
statistics efficiently, individual changes that are considered should be
small, such that the changes in the statistic values can be easily computed.
The smallest change that we can make in a graph (without changing the node
set) is to add or remove an edge. In fact, the change in subgraph count
statistics for the addition and removal of an edge can be expressed in terms
of the immediate properties of the two involved nodes. Furthermore, in order
to optimise the algorithm, the computation of the changes in statistics over
all changes considered in one step makes use of common subexpression. Thus, we
consider at each step the additions and removals of edges connected to one
single node. As we show in the next section, this leads to efficient
expressions for the change in the subgraph count statistics. We also note that
the described algorithm, while performing the changes with best reduction in
error at each step, is not a pure greedy algorithm; indeed, the algorithm
continues to iterate in the case that the error $e$ cannot be further reduced
by any single edge switch.

Algorithm 1 The small graph generation algorithm in non-vectorised form.

0: a size $n^{\prime}$

0: a set of statistic labels $\mathcal{S}=\\{m,s,z,x,t,q\\}$ and their values
$\mathbf{x}$ for the input graph

0: a convergence parameter $\epsilon>0$

0: a graph $G^{\prime}=(V^{\prime},E^{\prime})$

$G^{\prime}=\textrm{ErdősRényi}(n^{\prime},m^{\prime}/{n^{\prime}\choose 2})$

for all $S\in\mathcal{S}$ do

$\mathbf{y}^{S}=S(G^{\prime})$

end for

repeat

Choose a node $u\in V^{\prime}$ at random

for all $S\in\mathcal{S}$ do

for all $w\in V^{\prime}\setminus\\{u\\}$ do

$\Delta^{S}_{w}=S(G^{\prime}\pm\\{u,w\\})-S(G^{\prime})$

end for

end for

$v=\mathrm{argmin}_{w\in
V^{\prime}\setminus\\{u\\}}\sum_{S\in\mathcal{S}}((\mathbf{y}^{S}+\Delta^{S}_{w}-\mathbf{x}^{S})/\mathbf{x}^{S})^{2}$

$G^{\prime}=G^{\prime}\pm\\{u,v\\}$

for all $S\in\mathcal{S}$ do

$\mathbf{y}^{S}=\mathbf{y}^{S}+\Delta^{S}_{v}$

end for

$e=\sum_{S\in\mathcal{S}}((\mathbf{y}^{S}-\mathbf{x}^{S})/\mathbf{x}^{S})^{2}$

until $e$ has not achieved a new minimum value in the last
$(-n^{\prime}\log\epsilon)$ iterations

The function $\textrm{ErdősRényi}(n^{\prime},p)$ generates an Erdős–Rényi
graph with $n^{\prime}$ vertices and individual edge probability $p$.
$G^{\prime}\pm\\{u,v\\}$ denotes the graph $G^{\prime}$ in which the state of
the edge $\\{u,v\\}$ has been switched, i.e. removed or added depending on
whether $\\{u,v\\}$ is present or not. The convergence parameter $\epsilon$
ensures that the expected ratio of nodes $u$ that were not visited since the
last new minimum value $e$ was found equals $\epsilon$. In all experiments, we
use a value of $\epsilon=0.01$.

Algorithm 1 requires to compute the difference in statistics between the
current graph $G^{\prime}$ and the graph $G^{\prime}$ in which one edge is
added or removed:

$\Delta^{S}_{w}=S(G^{\prime}\pm\\{u,w\\})-S(G^{\prime})$ (15)

In order to speed up the algorithm, this calculation is performed in a
vectorised way. The existence of closed-form expressions for $\Delta^{S}_{w}$
decides whether a particular statistic can be used in our algorithm. Since
$\Delta^{S}_{w}$ must be computed at each step for all nodes $w\in V$ (except
for $w=v$), we give vectorised expressions that give a vector $\Delta^{S}$
containing the value $\Delta^{S}_{w}$ for all $w\in V$. The individual value
computed for $w=v$ is then simply ignored. During the run of the algorithm,
the graph $G^{\prime}$ is always represented by its symmetric adjacency matrix
$\mathbf{A}\in\\{0,1\\}^{n\times n}$.

In what follows we provide expressions for the vectors $\Delta^{S}$ measuring
the change in statistic $S$ expressed as functions of the graph’s adjacency
matrix $\mathbf{A}$, for each statistic $S\in\\{m,s,z,x,t,q\\}$. These
expressions make use of the degree vector $\mathbf{d}$, which is updated along
with the matrix $\mathbf{A}$. $\mathbf{u}\circ\mathbf{v}$ denotes the entry-
wise product between two vectors $\mathbf{u}$ and $\mathbf{v}$, and
$\mathbf{A}_{:u}$, the $u$th column of $\mathbf{A}$.

Number of edges ($m$): The number of edges will always increase or decrease by
one, depending on the previous state of the edge $\\{u,w\\}$.

$\Delta^{m}=-2\mathbf{A}_{:u}+1$ (16)

Number of wedges ($s$): When adding an edge, the number of wedges increases by
the sum of the degrees of the two connected nodes. When removing an edge, the
number of wedges decreases by the sum of the degrees of the two nodes, minus
two.

$\Delta^{s}=\Delta^{m}\circ(\mathbf{d}+\mathbf{d}_{u})+2\mathbf{A}_{:u}$ (17)

Number of claws ($z$): When the two nodes $u$ and $w$ are not connected, the
number of added claws equals
$\mathbf{d}_{u}(\mathbf{d}_{u}-1)+\mathbf{d}_{w}(\mathbf{d}_{w}-1)$. When the
two nodes are connected, the number of removed wedges can be computed in the
same way, but based on the degrees after the removal.

$\displaystyle\Delta^{z}=\frac{1}{2}\Delta^{m}\circ[(\mathbf{d}-\mathbf{A}_{:u})\circ(\mathbf{d}-1-\mathbf{A}_{:u})$
(18)
$\displaystyle\qquad\qquad\qquad\qquad\qquad+(-\mathbf{A}_{:u}+\mathbf{d}_{u})\circ(-\mathbf{A}_{:u}-1+\mathbf{d}_{u})]$

Number of crosses ($x$): The expression for the number of $k$-stars for higher
$k$ follows the same pattern as for $k=2$ and $k=3$. Due to the asymmetry
between the addition and the removal of edges, the resulting expressions get
increasingly complex.

$\displaystyle\Delta^{x}=\frac{1}{6}[(\mathbf{d}-1)\circ(\mathbf{d}-2)\circ(\mathbf{A}_{:u}\circ(-2\mathbf{d}+3)+\mathbf{d})$
(19)
$\displaystyle\qquad\qquad\;+(\mathbf{d}_{u}-1)(\mathbf{d}_{u}-2)((3-2\mathbf{d}_{u})\mathbf{A}_{:u}+\mathbf{d}_{u})]$

Number of triangles ($t$): When adding an edge between two nodes, the number
of added triangles equals the number of common neighbours between the two
nodes. Likewise when removing an edge, the number of removed triangles equals
the number of common neighbours of the two nodes. We thus get the following
expression for the change in the number of triangles.

$\Delta^{t}=(\mathbf{A}\mathbf{A}_{:u})\circ\Delta^{m}$ (20)

We thus need to perform one sparse matrix-vector multiplication for each
iteration step.

Number of squares ($q$): To compute the number of squares added or removed, we
count the number of paths of length three added or removed between $u$ and
$w$. In principle, this can be achieved by using the expression for
$\Delta^{t}$ and multiplying $\mathbf{A}_{:u}$ once more by $\mathbf{A}$.
However, this will also include the number of paths of length three that
include an edge $\\{u,x\\}$ or $\\{w,x\\}$ multiple times, or that include the
edge $\\{u,w\\}$ if it is present. Thus, these cases must be subtracted to get
the correct number of squares added or removed.

$\Delta^{q}=(\mathbf{A}^{2}\mathbf{A}_{:u})\circ\Delta^{m}+\mathbf{A}_{:u}\circ(\mathbf{d}+\mathbf{d}_{u}-1)$
(21)

We thus need to perform two matrix-vector multiplications in this step.

##  4 Experimental Evaluation

In this section, we describe the user experiments conducted in order to
evaluate the two proposed graph visualisation methods against baseline
methods, as well as evaluate the runtime of SynGraphy empirically, and perform
an experiment to determine the influence of the size $n^{\prime}$ of the
generated graphs on the drawings.

FR LA SU SN NO SI

Figure 3:  The same network dataset, the Reactome network (Joshi-Topé et al,
2005), visualised using the six algorithms evaluated in this paper’s
experiment. The images are those used on the experimental website. The
abbreviations of the graph visualisation methods are those given in Table 2.
(FR–SN) Baseline methods. (NO–SI) SynGraphy methods.

###  4.1 Experimental Setup

The human experiment was conducted via a website111The site was available at
http://nohairball.west.uni-koblenz.de/, and is offline as of this writing. on
which practitioners in the fields of graph mining, network analysis, and
related fields were asked to participate. After showing the participants the
definitions of the graph statistics

$\displaystyle d,c,\delta,q,z,s,b$

as defined in Section 2.3, each participant is shown, in turn, twelve pages
organised as follows. After choosing one graph visualisation method at random
(see list in Table 2), two networks are chosen at random from a set of 18
predefined networks, and visualised using the chosen method. Then, a graph
statistic is chosen at random from the set of the seven mentioned graph
statistics; the participant is asked to determine which of the two shown
graphs has a higher or lower value for the given graph statistic. The layout
of the page (which graph is shown left and right), and whether the question is
formulated as _higher_ or _lower_ is chosen randomly. The user may also
determine that she cannot determine the difference by clicking on a button for
that purpose. The users’ answers are logged, and evaluated in the following.
The participants in the experiments were mostly from a computer science and
network science background, partly associated with the authors’ institution,
and partly external.

The experiment’s web page was available in mid 2016, and was disseminated
using various channels. An anonymised screenshot of the experiment running in
a browser is shown in Figure 4. For reference, Figure 3 shows a single network
dataset visualised with all six algorithms used in the experiments.

Table 2:  The six graph visualisation methods evaluated in the experiments. The last two are the variants of the SynGraphy algorithm described in this paper.  | Method
---|---
FR | Fruchterman–Reingold (1991)
LA | Laplacian embedding (see e.g. Koren, 2003)
SU | Uniform vertex sampling1
SN | Node sampling1
NO | SynGraphy with empirical statistic distribution1
SI | SynGraphy with size-independent statistics1
1These methods are followed by the use of the Fruchterman–Reingold method for
drawing the resulting small graph.
Figure 4:  An anonymised screenshot of the human experiment. The screenshot
shows a page of the experiment in which the human is tasked to determine
visually which of the two shown networks (left or right) has a higher degree
assortativity. The human subject must click on the network to answer the
question.  Figure 5:  The experimental results showing the success rate and
confusion rate of the six evaluated graph visualisation methods. The rightmost
two, NO and SI, are the variants of the SynGraphy algorithm described in this
paper.  Figure 6:  The experimental results showing the success rate and
confusion rate for each of the network statistics employed in the experiments.
In particular, the plots shows that as measured over all graph visualisation
methods, certain network statistics are very difficult to discern.  Figure 7:
The experimental results broken down both by graph visualisation method and
network statistic. The grey bars represent the confusion rate.

The experiment was divided into two parts. In the first part, we show to the
experimental subject the definitions of the various network properties and
show two example graphs with high and low values of each properties. In the
second part, some random graphs (generated using the methods discussed above)
were shown to the users. The user had to identify the graph with a low or high
value of a certain property. The six graph visualisation methods used in the
experiment are given in Table 2. For the uniform vertex sampling method, we
chose a random subset of the nodes, and keep all edges connected to at least
one of them (Clauset, 2014; Hu and Lau, 2013). For node sampling, we chose a
random subset of all nodes, and keep all edges connecting two of them.

###  4.2 Results

From the logs of the experiments, we determine for each page view whether the
person was correct or incorrect in their assessment of the given network
statistic for the two given graphs visualised with the given method. In total,
3,049 questions where answered by 79 unique identified participants, plus an
unknown number of anonymous participants.222Participants were given the
possibility to optionally enter their name into the system. Each question is
thus answered in exactly one of three ways: correct, incorrect, or skipped.
Thus, we collect three counts:

  * •

$N^{+}$, the number of correct answers,

  * •

$N^{-}$, the number of incorrect answers,

  * •

$N^{0}$, the number of skipped answers.

The total number of questions shown is then $N^{+}+N^{-}+N^{0}$. From these
numbers, we compute the success rate $R^{s}$ and the confusion rate $R^{c}$ as
follows:

$\displaystyle R^{s}=N^{+}/(N^{+}+N^{-}+N^{0}),\qquad
R^{c}=N^{0}/(N^{+}+N^{-}+N^{0})$ (22)

The experimental results are shown in figures 5, 6 and 7.

Success Rate by Visualisation Method. Figure 5 shows the overall success rate
for each graph visualisation method used in the experiment, including the two
SynGraphy variants. The result shows that the overall success rate of the
SynGraphy variant based on normalisation (NO) is significantly higher than all
the other methods used in the experiment ($\chi^{2}$ test, $p=0.000287$). The
SynGraphy variant based on statistics that are by construction independent of
size (SI) does not show statistically better results than the node sampling
(SN) method ($\chi^{2}$ test, $p=0.447$) and the uniform vertex sampling (SU)
method ($\chi^{2}$ test, $p=0.984$). The fact that the empirical normalisation
method (NO) performs much better than the one using predetermined normalised
statistics (SI) is also enlightening, as it suggests that statistics such as
the clustering coefficient, which are size-independent in theory by
construction, are not necessarily so in practice, and that in practice it is
more effective to normalise subgraph count statistics using available
empirical data from actual network datasets.

Success Rate by Network Statistic. Comparing the individual network statistics
with each other, it appears that the task of identifying structural network
properties is, as a general rule, a difficult one. As seen in Figure 6, only a
single statistic reaches an overall success rate (over all visualisation
methods combined) of more than 50%, this being the clustering coefficient.
This may be explained by the better knowledge of the clustering coefficient of
many participants, who, despite the explanations given at the start of the
experiments, may not be familiar with less well-known network statistics such
as the degree assortativity, and may thus have a harder time recognising it in
graph drawings. It can also be observed that non-normalised statistics such as
the raw number of wedges ($s$) have low success rates – justifying the need
for normalised statistics such as the clustering coefficient.

Confusion Rate. In Figure 5, we may also compare the confusion rate of the
different graph visualisation methods. The interpretation of the confusion
rate is ambivalent: A high confusion rate may indicate that a particular graph
visualisation method is confusing for the user, but at the same time a graph
visualisation method should avoid misleading users – it is better if the user
is undecided, than for her to be wrong. Thus, the confusion rate should not be
taken as an indicator of the quality of a graph visualisation method, but
rather as a measure of how convincing the resulting plots look. Thus, a graph
visualisation method with a low confusion rate may lead to the user often
thinking that they are able to recognise the properties of the graph,
regardless of whether they are right or wrong. That being said, the smallest
confusion rates are achieved by the two sampling methods (SU, SN), and the
highest by the direct drawing methods (FR, LA). The two SynGraphy variants
then achieve an intermediate confusion rate. To be precise, the size-
independent SynGraphy (SI) method has a significantly lower confusion rate
than the Fruchterman–Reingold (FR) method ($\chi^{2}$ test,
$p=2.78{\times}10^{-7}$), and also than the four baseline methods combined
($\chi^{2}$ test, $p=0.00356$).

Individual Results. The results shown in Figure 7 allow us to make more
detailed remarks about the performance of the individual graph visualisation
methods with respect to individual network statistics. We observe that the
SynGraphy variant based on empirically derived normalisation (NO) gives the
highest success rate for the number of claws and density, while for other
network statistics, other visualisation methods than SynGraphy perform better,
although no single one is better for many network statistics. In particular,
the plain Fruchterman–Reingold method has the highest success rate for the
number of wedges and squares, the Laplacian embedding has the highest success
rate for the diameter, uniform vertex sampling has the highest success rate
for the clustering coefficient and assortativity, and finally node sampling
has the highest success rate for the bipartivity coefficient. Thus, it seems
that if a graph visualisation is needed in which a particular network property
should be emphasized, then we cannot claim that the SynGraphy method is the
best choice in all cases. It should also be noted that if a particular
property such as bipartivity is to be emphasized, then other, more specific
methods may be chosen.

###  4.3 Runtime Experiments

In order to assess the runtime of the method empirically, we show the runtime
of the SynGraphy algorithm in Figure 8. The two variants of SynGraphy have,
for all practical purposes, the same runtime, as they differ only in the step
of scaling down the network statistics, whose runtime is constant and very
small, i.e., it does not depend on the network size. The amount of memory
taken up by a sparse graph is proportional to the number of edges rather than
the number of nodes, and thus we measure the runtime in function of the number
of edges, i.e., the graph’s volume $m$.

Figure 8:  The runtime of the SynGraphy algorithm for each tested network
dataset, as a function of the number of edges in the network. For small
networks (number of edges $m<10^{4}$) the runtime is dominated by the constant
runtime of the scaling down, small graph generation and graph drawing steps of
the SynGraphy algorithm, while for larger networks ($m>10^{4}$), the runtime
is empirically quadratic due to the subgraph count computation step, in
accordance with the theoretical model described in the text.

We can derive the runtime of the SynGraphy algorithm in function of the input
graph’s volume $m$ as follows. The different steps of the algorithm are in
turn: (i) computation of subgraph counts, (ii) scaling down of subgraph
counts, (iii) generation of the synthetic graph, and (iv) generation of the
graph layout using the Fruchterman–Reingold algorithm. Steps (ii–iv) have
runtime independent of the graph’s number of edges $m$, as they are performed
on data of fixed size, and thus have a formal runtime of $O(1)$. Step (i), the
counting of subgraphs, is the only part of the algorithm depending on $m$, and
its actual runtime is dominated by the computation of the counting of squares
(4-cycle subgraphs), which is quadratic. Thus, the overall runtime of
SynGraphy is expected on theoretical grounds to be $O(m^{2})$. The empirical
data as shown in Figure 8 is consistent with that runtime of $O(m^{2})$.

###  4.4 Varying $n^{\prime}$

In order to assess the influence of the size $n^{\prime}$ of the generated
small synthetic graphs on the quality of the summarisation, we generate graph
summarisations for varying $n^{\prime}$ using a single network, the email
network of Enron (Klimt and Yang, 2004) already used in Figure 1, visualised
using the SynGraphy variant based on the empirically learned normal
distribution of statistics (NO). The resulting graphs are shown in Figure 9,
for varying $n^{\prime}$ between a value of $n^{\prime}=10$ up to a value of
$n^{\prime}=1000$.

$n^{\prime}=10$ $15$ $20$ $30$ $50$ $70$ $100$ $150$ $200$ $300$ $500$ $700$
$1000$

Figure 9:  The email dataset from (Klimt and Yang, 2004) summarised with a
varying number of nodes $n^{\prime}$ in the resulting graph, using SynGraphy
with the empirical statistic distribution (NO). When $n^{\prime}$ is very
small, the generated graphs are uninformative, while for too large
$n^{\prime}$, the drawings converge to a “hairball”-like drawing. For
intermediate values of $n^{\prime}$ in the range $70\leq n^{\prime}\leq 100$,
the drawings strike a balance between enough detail and enough space to make
the network structure evident.

We observe that for very small $n^{\prime}$, the small number of drawn nodes
and edges is not enough to convey any meaningful information about the
network, and thus the summarisation of the original graph is not a useful
representation of the original network. For too large values of $n^{\prime}$,
the drawn graphs approach a “hairball”-like image, and thus they too are not
useful representations of the graph to be summarised. Intermediate values in
the approximate range $70\leq n^{\prime}\leq 100$ give an acceptable trade-off
between the two conditions, leading us to recommend a value of 80 for
$n^{\prime}$. This is also the value used in the previous experiments, as
explained in the previous section.

##  5 Discussion

We have seen that in the overall evaluation the proposed SynGraphy variant
based on the empirical distribution of graph statistics (NO) has the highest
success rate of all tested visualisation methods, and we may thus recommend it
as a first step in the analysis of networked data. However, if certain
specific graph properties are to be emphasized, then traditional graph drawing
methods may be best instead – our experiments showed for instance that a
Laplacian embedding is well suited to making evident the diameter of a graph.
We must of course stress that this method is only useful if general,
structural properties are to be summarised – an exploration of the actual
nodes and edges of the network is, by construction, not possible with graph
summarisation methods. The problem of visualising the individual nodes and
edges of a given large network is a separate problem that, in the end, cannot
be solved without showing only a small subset of the given network, giving
rise to the problem of letting the user search the nodes of a graph, for
instance using a search engine or even a recommender system, if the nodes of
the graph are amenable to that – that problem is however clearly outside of
the scope of graph summarisation.

Finally, we should address the question whether graph visualisation and
summarisation are a necessary part of the graph mining process. If certain
specific properties such as the bipartivity or assortativity of a given
network are to be detected, it is usually more productive to simply compute
and show the corresponding numerical statistics. Graph visualisation and
summarisation however are usually used for a higher-level purpose: As recently
shown (Matejka and Fitzmaurice, 2017), a general-purpose visualisation of a
given dataset is often a simple and effective way of starting an in-depth
analysis. Network datasets are of course not exempt from this, and thus a
graph summarisation method such as SynGraphy can be recommended to fill that
role, and, as we showed in our experimental evaluation, give practitioners a
better idea of a given network dataset’s overall structural properties.

An evident future extension of our method is to apply it to further graph
properties, although with more elaborate graph statistics, the task of graph
generation becomes more runtime-intensive, as the vectorised expressions as
given in Section 3.2 are then not possible anymore.

## References

  * Cem et al (2013) Cem E, Tozal ME, Sarac K (2013) Impact of sampling design in estimation of graph characteristics. In: Proc. Perform. Comput. and Commun. Conf., pp 1–10
  * Clauset (2014) Clauset A (2014) Network analysis and modeling, small worlds and navigable networks. Lecture notes
  * Fruchterman and Reingold (1991) Fruchterman TMJ, Reingold EM (1991) Graph drawing by force-directed placement. Software: Pract and Exp 21(11):1129–1164
  * van Ham et al (2009) van Ham F, Schulz HJ, Dimicco JM (2009) Honeycomb: Visual analysis of large scale social networks. In: Proc. Conf. Human-Computer Interaction, pp 429–442
  * Harrigan (2014) Harrigan M (2014) A network analyst’s view of the block chain. http://www.coindesk.com/network-analysts-view-block-chain/
  * Hu and Lau (2013) Hu P, Lau WC (2013) A survey and taxonomy of graph sampling. CoRR abs/1308.5865
  * Jia et al (2008) Jia Y, Hoberock J, Garland M, Hart JC (2008) Visualization of social and other scale-free networks. Trans Visualization and Computer Graphics 14(6):1285–1292
  * Johnson and Wichern (2007) Johnson RA, Wichern DW (2007) Applied Multivariate Statistical Analysis, 6th edn. Pearson
  * Joshi-Topé et al (2005) Joshi-Topé G, Gillespie M, Vastrik I, D’Eustachio P, Schmidt E, de Bono B, Jassal B, Gopinath G, Wu G, Matthews L, Lewis S, Birney E, Stein L (2005) Reactome: A knowledgebase of biological pathways. Nucleic Acids Res 33(Database):D428–D432
  * Klimt and Yang (2004) Klimt B, Yang Y (2004) The Enron corpus: A new dataset for email classification research. In: Proc. Eur. Conf. on Mach. Learn., pp 217–226
  * Koren (2003) Koren Y (2003) On spectral graph drawing. In: Proc. Int. Comput. and Comb. Conf., pp 496–508
  * Krzywinski et al (2012) Krzywinski M, Birol I, Jones SJM, Marra MA (2012) Hive plots – rational approach to visualizing networks. Brief Bioinform 13(5):627–644
  * Kunegis (2013) Kunegis J (2013) KONECT – The Koblenz Network Collection. In: Proc. Int. Conf. on World Wide Web Companion, pp 1343–1350
  * Kunegis (2015) Kunegis J (2015) Exploiting the structure of bipartite graphs for algebraic and spectral graph theory applications. Internet Math 11(3):201–321
  * Latapy et al (2008) Latapy M, Magnien C, Vecchio ND (2008) Basic notions for the analysis of large two-mode networks. Soc Netw 30(1):31–48
  * Liu et al (2017) Liu Y, Dighe A, Safav T, Koutra D (2017) Graph summarization: A survey. ACM Comput Surv In review
  * Matejka and Fitzmaurice (2017) Matejka J, Fitzmaurice G (2017) Same stats, different graphs: Generating datasets with varied appearance and identical statistics through simulated annealing. In: Proc. Conf. on Human Factors in Comput. Syst., pp 1290–1294
  * Newman (2003a) Newman MEJ (2003a) Mixing patterns in networks. Phys Rev E 67:026,126
  * Newman (2003b) Newman MEJ (2003b) The structure and function of complex networks. SIAM Rev 45(2):167–256
  * Nocaj (2015) Nocaj A (2015) Untangling networks: Focus on less to see more. PhD thesis, Univ. of Konstanz
  * Nocaj et al (2015) Nocaj A, Ortmann M, Brandes U (2015) Untangling the hairballs of multi-centered, small-world online social media networks. J of Graph Algorithms and Appl 19(2):595–618
  * Purchase et al (1995) Purchase HC, Cohen RF, James M (1995) Validating graph drawing aesthetics. In: Proc. Graph Drawing Symp., pp 435–446
  * Redmond and Cunningham (2013) Redmond U, Cunningham PD (2013) A temporal network analysis reveals the unprofitability of arbitrage in the Prosper marketplace. Expert Syst with Appl: An Int J 40(9)
  * Ruciński (1988) Ruciński A (1988) When are small subgraphs of a random graph normally distributed? Prob Th Rel Fields 78:1–10
  * Schulz and Hurter (2013) Schulz HJ, Hurter C (2013) Workshop: Grooming the hairball – how to tidy up network visualizations? In: Proc. Inf. Vis. Conf.
  * Shawn et al (2011) Shawn M, Brown WM, Klavans R, Boyack KW (2011) OpenOrd: An open-source toolbox for large graph layout. In: Proc. Vis. and Data Anal.