text:
\section{Introduction}

A fundamental problem in graph mining is to
identify the most central nodes in a graph.
Numerous centrality measures have been proposed,
including
degree centrality,
closeness centrality~\cite{closenesscentrality},
betweenness centrality~\cite{betweennesscentrality},
random-walk centrality~\cite{randomwalkcentrality},
Katz centrality~\cite{katzcentrality}, and
PageRank~\cite{PageRank}.

In the interest of robustness
many centrality measures use random walks:
while the shortest-path distance between two nodes can
change dramatically by inserting or deleting a single edge,
distances based on random walks account for multiple paths and offer a
more global view of the connectivity between two nodes.
In this spirit, the random-walk centrality of {\it one} node with
respect to {\it all nodes} of the graph is defined as the
expected time needed to come across this node in a random walk that
starts in any other node of the graph~\cite{randomwalkcentrality}.

In this paper, we consider a measure that generalizes random-walk
centrality for a {\it set} of nodes $C$
with respect to a set of \emph{query nodes} $Q$.
Our centrality measure is defined as the expected length
of a random walk that starts from any node in $Q$ until
it reaches any node in $C$ --- at which point the random walk
is {\it ``absorbed''} by $C$. Moreover, to allow for adjustable
importance of query nodes in the centrality measure,
we consider random walks {\it with restarts}, that occur with a fixed
probability $\alpha$ at each step of the random walk.
The resulting computational problem
is to find a set of $k$ nodes $C$ that optimizes this measure with respect to nodes $Q$,
which are provided as input.
We call this measure  $k$ {\em absorbing random-walk centrality}
and the corresponding optimization problem {{{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}}.

To motivate the {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem, let us consider the scenario of
searching the Web graph and summarizing the search results.
In this scenario, nodes of the graph correspond to webpages,
edges between nodes correspond to links between pages, and the set of
query nodes $Q$ consists of all nodes that match a user query,
i.e.,  all webpages that satisfy a keyword search.
Assuming that the size of $Q$ is large, the goal is to find the $k$
most central nodes with respect to $Q$,
and present those to the user.

It is clear that ordering the nodes of the graph by their individual
random-walk centrality scores and taking the top-$k$ set does not solve the
{{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem,
as these nodes may all be located in the
same ``neighborhood'' of the graph, and thus,
may not provide a good absorbing set for the query.
On the other hand,
as the goal is to minimize the expected absorption time for walks starting at $Q$,
the optimal solution to the {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem will be a set of $k$, both
{\em centrally-placed} and {\em diverse}, nodes.

This observation has motivated researchers in the informa\-tion-retrieval
field to consider random walks with absorbing states in order
to diversify web-search results~\cite{Zhu:2007tj}.
However, despite the fact that similar problem definitions and algorithms have
been considered earlier, the {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem has not been formally
studied and there has not been a theoretical analysis of its properties.

Our key results in this paper are the following:
we show that the {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem is {\ensuremath{\mathbf{NP}}}-hard, and
we show that the $k$ absorbing random-walk centrality measure is
{\em monotone} and {\em supermodular}.
The latter property allows us to quantify the approximation guarantee
obtained by a natural greedy algorithm,
which has also been considered by previous work~\cite{Zhu:2007tj}.
Furthermore, a na\"{i}ve implementation of the greedy algorithm requires
many expensive matrix inversions, which make the algorithm particularly slow.
Part of our contribution is to show how to make use of the
Sherman-Morrison inversion formula to implement the greedy algorithm
with only one matrix inversion and more efficient
matrix$\,\times\,$vector multiplications.

Moreover, we explore the performance of faster, heuristic algorithms,
aiming to identify methods that are faster than the greedy approach
without significant loss in the quality of results. The heuristic
algorithms we consider include the personalized PageRank
algorithm~\cite{PageRank,langville2005survey}
as well as algorithms based on spectral clustering~\cite{von2007tutorial}. We find that, in
practice, the personalized PageRank algorithm offers a very good trade-off
between speed and quality.

The rest of the paper is organized as follows.
In Section~\ref{sec:related},
we overview previous work and discuss how it compares to this paper.
We define our problem in Section~\ref{section:problem-definition}
and provide basic background results on absorbing random walks in
Section~\ref{section:absorbing}.
Our main technical contributions are given in Sections~\ref{section:absorbing} and~\ref{sec:properties},
where we characterize the complexity of the problem, and provide the details of the
greedy algorithm and the heuristics we explore.
We evaluate the performance of algorithms in Section~\ref{sec:experiments},
over a range of real-world graphs,
and Section~\ref{sec:conclusions} is a short conclusion.
Proofs for some of the theorems shown in the paper are
provided in the Appendix.

\section{Related work}
\label{sec:related}

Many works in the literature explore ways to quantify
the notion of node centrality on graphs~\cite{boldi2014axioms}.
Some of the most commonly-used measures include the following:
($i$) {\em degree centrality}, where the centrality of a node is
simply quantified by its degree;
($ii$) {\em closeness centrality} \cite{leavitt1951some,closenesscentrality}, defined
as the average distance of a node from all other nodes on the
graph;
($iii$) {\em betweenness centrality}~\cite{betweennesscentrality},
defined as the number of shortest paths between pairs of
nodes in the graph that pass through a given node;
($iv$) {\em eigenvector centrality}, defined as the stationary probability
that a Markov chain on the graph visits a given node, with
Katz centrality~\cite{katzcentrality} and PageRank~\cite{PageRank}
being two well-studied variants; and
($v$) {\em random-walk centrality}~\cite{randomwalkcentrality}, defined as
the expected first passage time of a random walk from a given node,
when it starts from a random node of the graph.
The measure we study in this paper
generalizes the notion of {\em random-walk centrality} to a set
of absorbing nodes.

Absorbing random walks have been used in previous work
to select a {\it diverse} set of nodes from a graph.
For example,
an algorithm proposed by Zhu et~al.~\cite{Zhu:2007tj} selects nodes in the following manner:
($i$) the first node is selected based on its PageRank value and is set as absorbing;
($ii$) the next node to be selected is the node that maximizes the expected first-passage time from the already selected absorbing nodes.
Our problem definition differs considerably from the one considered in that work,
as in our work the expected first-passage times are always computed from the set of
query nodes that are provided in the input, and not from the nodes that participate in the solution so far.
In this respect, the greedy method proposed by Zhu et~al.\ is not associated with a crisp problem definition.

Another conceptually related line of work aims to
select a diverse subset of query results, mainly
within the context of document retrieval~\cite{Agrawal:2009bc,Angel:2011be,Vieira:2011bx}.
The goal, there, is to select $k$ query results to
optimize a function that quantifies the trade-off between relevance
and diversity.

Our work is also remotely related to the problem studied by Leskovec et al.\ on
{\it cost-effective outbreak detection} \cite{leskovec2007cost}.
One of the problems discussed there is to select
nodes in the network so that the detection time for a set of cascades is minimized.
However, their work differs from ours on the fact that they consider as input a set of {\it cascades},
each one of finite size,
while in our case the input consists of a set of query {\it nodes}
and we consider a probabilistic model that generates random walk paths,
of possibly infinite size.

\section{Problem definition}
\label{section:problem-definition}

We are given a graph $G=(V,E)$ over a set of nodes $V$ and set of
undirected edges $E$. The number of nodes $|V|$ is denoted by $n$ and the
number of edges $|E|$ by $m$. The input also includes a subset of
nodes $Q\subseteq V$, to which we refer as the {\em query nodes}.
As a special case, the set of query nodes $Q$ may be equal to the
whole set of nodes, i.e., $Q=V$.

Our goal is to find a set $C$ of $k$
nodes that are {\em central} with respect to the query nodes $Q$.
For some applications it makes sense to restrict the central nodes to
be only among the query nodes, while in other cases,
the central nodes may include any node in $V$.
To model those different scenarios, we consider a set of candidate nodes $D$,
and require that the $k$ central nodes should belong in this candidate set,
i.e., $C\subseteq D$.
Some of the cases include $D=Q$, $D=V$, or $D=V\setminus Q$, but
it could also be that $D$ is defined in some other way
that does not involve $Q$.
In general, we assume that $D$ is given as input.

The {\em centrality} of a set of nodes $C$ with respect to query nodes $Q$
is based on the notion of absorbing random-walks and their expected
length.
More specifically, let us consider a random walk on the nodes $V$ of
the graph, that proceeds at discrete steps: the walk starts from a
node $q\in Q$ and, at each step
moves to a different node, following edges in $G$,
until it arrives at some node in $C$.
The {\em starting} node $q$ of the walk is chosen according to a probability distribution ${\ensuremath{\mathbf{s}}}$.
When the walk arrives at a node $c\in C$ for the first time,
it terminates, and we say that the random walk is {\em absorbed} by that node~$c$.
In the interest of generality, and to allow for adjustable
importance of query nodes in the centrality measure,
we also allow the random walk to restart. Restarts occur with a
probability $\alpha$ at each step of the random walk, where $\alpha$ is
a parameter that is specified as input to the problem.
When restarting, the walk proceeds to a query node selected randomly
according to ${\ensuremath{\mathbf{s}}}$.
Intuitively, larger values of $\alpha$ favor nodes that are closer to nodes $Q$.

We are interested in the expected length (i.e., number of steps)
of the walk that starts from a query node $q\in Q$ until it gets
absorbed by some node in~$C$, and we denote this expected length
by ${\ensuremath{\mathrm{ac}}}^q_{_Q}(C)$.
We then define the {\em absorbing random-walk centrality} of a set of
nodes $C$ with respect to query nodes~$Q$,~by
\[
{\ensuremath{\mathrm{ac}}}_Q(C) = \sum_{q\in Q} {\ensuremath{\mathbf{s}}}(q) \, {\ensuremath{\mathrm{ac}}}^q_{_Q}(C).
\]

The problem we consider in this paper is the following.
\begin{problem}
{\em (}{{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}{\em )}
\label{problem:k-ac}
We are given a graph $G=(V,E)$,
a set of query nodes $Q \subseteq V$,
a set of candidate nodes $D\subseteq V$,
a starting probability distribution ${\ensuremath{\mathbf{s}}}$ over $V$ such that ${\ensuremath{\mathbf{s}}}(v)=0$ if $v \in V \setminus Q$,
a restart probability~$\alpha$,
and an integer $k$.
We ask to
find a set of $k$ nodes $C\subseteq D$
that {\em minimizes} ${\ensuremath{\mathrm{ac}}}_Q(C)$, i.e.,
the expected length of a random walk that starts from $Q$ and proceeds
until it gets absorbed in some node in $C$.
\end{problem}
In cases where we have no reason to distinguish among the query nodes,
we consider the uniform starting probability distribution
${\ensuremath{\mathbf{s}}}(q) = 1/|Q|$.
In fact, for simplicity of exposition, hereinafter we focus on the
case of uniform distribution. However, we note that all our definitions
and techniques generalize naturally, not only to general starting
probability distributions ${\ensuremath{\mathbf{s}}}(q)$, but also to
{\em directed} and {\em weighted} graphs.

\section{Absorbing random walks}
\label{section:absorbing}

In this section we review some relevant background on absorbing random
walks. Specifically, we discuss how to calculate the objective function
${\ensuremath{\mathrm{ac}}}_Q(C)$ for Problem~\ref{problem:k-ac}.

Let  $\mathbf{P}$ be the transition matrix for a random walk,
with $\mathbf{P}(i,j)$ expressing the probability that the random
walk will move to node~$j$ given that it is currently at node $i$.
Since random walks can only move to absorbing nodes $C$, but not away
from them, we set
$\mathbf{P}(c,c)=1$ and
$\mathbf{P}(c,j)=0$, if $j\neq c$,
for all absorbing nodes $c\in C$.
The set $T = V \setminus C$ of non-absorbing nodes is called {\em transient}.
If $N(i)$ are the neighbors
of a node $i \in T$ and $d_i = |N(i)|$ its degree, the transition
probabilities from node $i$ to other nodes are
\begin{equation}
\mathbf{P}(i,j) = \left\{
\begin{array}{ll}
\alpha \, {\ensuremath{\mathbf{s}}}(j) & \text{ if } j \in Q \setminus N(i), \\
(1-\alpha) / d_i+ \alpha \, {\ensuremath{\mathbf{s}}}(j)  & \text{ if } j \in N(i).
\end{array}
\right.
\end{equation}
Here, ${\ensuremath{\mathbf{s}}}$ represents the starting probability vector.
For example, for the uniform distribution over query nodes we have
${\ensuremath{\mathbf{s}}}(i) = 1/|Q|$ if $i\in Q$ and $0$ otherwise.
The transition matrix of
the random walk can be
written as follows
\begin{equation}
\mathbf{P} = \left(
\begin{array}{cc}
\mathbf{P}_{TT} & \mathbf{P}_{TC} \\
\mathbf{0} & \mathbf{I}
\end{array}
\right).
\end{equation}
In the equation above, $\mathbf{I}$ is an $(n - |T|)\times(n - |T|)$
identity matrix and $\mathbf{0}$ a matrix with all its entries equal
to $0$;
$\mathbf{P}_{TT}$ is the $|T| \times |T|$ sub-matrix of $\mathbf{P}$
that contains the transition probabilities between transient nodes; and
$\mathbf{P}_{TC}$ is the $|T|\times|C|$ sub-matrix of $\mathbf{P}$
that contains the transition probabilities from transient to absorbing nodes.

The probability of the walk being on node $j$ at exactly
$\ell$ steps having started at node~$i$,
is given by the $(i,j)$-entry of the matrix $\mathbf{P}_{TT}^\ell$.
Therefore, the expected total number of times that the random walk visits
node $j$ having started from node~$i$ is given by the $(i,j)$-entry of
the $|T|\times |T|$ matrix
\begin{equation}
\mathbf{F} = \sum_{\ell=0}^{\infty} \mathbf{P}_{TT}^\ell
= \left( \mathbf{I} - \mathbf{P}_{TT}\right)^{-1},
\end{equation}
which is known as the {\em fundamental matrix} of the
absorbing random walk.
Allowing the possibility to start the random walk at an absorbing node
(and being absorbed immediately), we see that the expected length of a
random walk that starts from node $i$ and gets absorbed by the set $C$
is given by the $i$-th element of the following $n\times 1$ vector
\begin{equation}
\mathbf{L} = \mathbf{L}_C = \left(
\begin{array}{c}
\mathbf{F} \\
\mathbf{0}
\end{array}
\right) \mathbf{1},
\end{equation}
where $\mathbf{1}$ is an $T\times 1$ vector of all 1s.
We write $\mathbf{L}=\mathbf{L}_C$ to emphasize the dependence on the
set of absorbing nodes $C$.

The expected number of steps when starting from a node in $Q$ and
until being absorbed by some node in $C$ is then obtained by summing
over all query nodes, i.e.,
\begin{equation}
\label{eq:inversion}
{\ensuremath{\mathrm{ac}}}_Q(C) =  {\ensuremath{\mathbf{s}}}^T\, \mathbf{L}_C.
\end{equation}

\subsection{Efficient computation of absorbing centrality}
\label{section:fasteval}

Equation~(\ref{eq:inversion}) pinpoints the difficulty of the problem we consider:
even computing the objective function ${\ensuremath{\mathrm{ac}}}_Q(C)$ for a
candidate solution $C$ requires an expensive matrix inversion;
$\mathbf{F} = \left( \mathbf{I} - \mathbf{P}_{TT}\right)^{-1}$.
Furthermore,
searching for the optimal set $C$
involves an exponential number of candidate sets,
while evaluating each one of them requires a matrix inversion.

In practice, we find that we can compute ${\ensuremath{\mathrm{ac}}}_Q(C)$ much faster approximately, as shown in Algorithm~\ref{algo:appox_ac}.
The algorithm follows from the infinite-sum expansion of Equation~(\ref{eq:inversion}).
\begin{eqnarray*}
{\ensuremath{\mathrm{ac}}}_Q(C) =  {\ensuremath{\mathbf{s}}}^T\, \mathbf{L}_C = {\ensuremath{\mathbf{s}}}^T\,
\left(
\begin{array}{c}
\mathbf{F} \\
\mathbf{0}
\end{array}
\right) \mathbf{1}
= {\ensuremath{\mathbf{s}}}^T\,
\left(
\begin{array}{c}
\sum_{\ell=0}^{\infty}\mathbf{P}_{TT}^\ell \\
\mathbf{0}
\end{array}
\right) \mathbf{1} \\
= {\ensuremath{\mathbf{s}}}^T\,
\sum_{\ell=0}^{\infty}
\left(
\begin{array}{c}
\mathbf{P}_{TT}^\ell \\
\mathbf{0}
\end{array}
\right) \mathbf{1}
= \left(\sum_{\ell=0}^{\infty}  {\ensuremath{\mathbf{s}}}^T\,
\left(
\begin{array}{c}
\mathbf{P}_{TT}^\ell \\
\mathbf{0}
\end{array}
\right)\right) \mathbf{1} \\
= \left(\sum_{\ell=0}^{\infty}  \mathbf{x}_\ell \right) \mathbf{1}
= \sum_{\ell=0}^{\infty}  \mathbf{x}_\ell \mathbf{1},
\end{eqnarray*}
with
\begin{equation}
\mathbf{x}_0 = {\ensuremath{\mathbf{s}}}^{^T} \;\text{ and }\; \mathbf{x}_{\ell+1} = \mathbf{x}_\ell \left(
\begin{array}{c}
\mathbf{P}_{TT} \\
\mathbf{0}
\end{array}
\right).
\end{equation}
Note that computing each vector $\mathbf{x}_\ell$ requires time~$\mathcal{O}(n^2)$.
Algorithm~\ref{algo:appox_ac} terminates
when the increase of the sum
due to the latest term
falls below a pre-defined threshold $\epsilon$.

\begin{algorithm}[t]
\caption{{{\sf\small  Approximate\-AC}}}
\label{algo:appox_ac}
\begin{algorithmic}
\STATE {\bf Input}: Transition matrix $\mathbf{P}_{TT}$, threshold $\epsilon$,
\\ starting probabilities ${\ensuremath{\mathbf{s}}}$
\STATE {\bf Output}: Absorbing centrality ${\ensuremath{\mathrm{ac}}}_Q$
\STATE $\mathbf{x_0} \leftarrow {\ensuremath{\mathbf{s}}}^{^T}$
\STATE $\delta \leftarrow \mathbf{x_0}\cdot\mathbf{1}$
\STATE ${\ensuremath{\mathrm{ac}}} \leftarrow  \delta$
\STATE $\ell \leftarrow 0$
\WHILE {$\delta < \epsilon$}
\STATE $\mathbf{x_{\ell + 1}} \leftarrow \mathbf{x_\ell} \left(
\begin{array}{c}
\mathbf{P}_{TT} \\
\mathbf{0}
\end{array}
\right)$
\STATE $\delta \leftarrow \mathbf{x_{\ell + 1}}\cdot \mathbf{1}$
\STATE ${\ensuremath{\mathrm{ac}}} \leftarrow {\ensuremath{\mathrm{ac}}} + \delta$
\STATE $\ell \leftarrow \ell + 1$
\ENDWHILE
\STATE \textbf{return} {\ensuremath{\mathrm{ac}}}
\end{algorithmic}
\end{algorithm}

\section{Problem characterization}
\label{sec:properties}

We now study the {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem in more detail.
In particular, we show that the function ${\ensuremath{\mathrm{ac}}}_Q$
is monotone and supermodular, a property that is used later to provide
an approximation guarantee for the greedy algorithm.
We also show that {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ is  {\ensuremath{\mathbf{NP}}}-hard.

Recall that a function $f:2^V\rightarrow\mathbb{R}$ over subsets of
a ground set $V$ is {\em submodular} if it has the
{\em diminishing returns} property
\begin{equation}
f(Y \cup \{u\}) - f(Y) \le f(X \cup \{u\}) - f(X),
\end{equation}
for all $X\subseteq Y\subseteq V$ and $u\not\in Y$.
The function $f$ is {\em supermodular} if $-f$ is submodular.
Submodularity (and supermodularity) is a very useful property for
designing algorithms.
For instance, minimizing a submodular function is a polynomial-time solvable problem,
while the maximization problem is typically amenable to approximation algorithms,
the exact guarantee of which depends on other properties of the function and requirements of the problem,
e.g., monotonicity, matroid constraints, etc.

Even though the objective function ${\ensuremath{\mathrm{ac}}}_Q(C)$ is
given in closed-form by Equation~(\ref{eq:inversion}),
to prove its properties we find it more convenient to work with its
descriptive definition, namely,
${\ensuremath{\mathrm{ac}}}_Q(C)$ being the expected length
for a random walk starting at nodes of $Q$
before being absorbed at nodes of~$C$.

For the rest of this section we consider that the set of query nodes
$Q$ is fixed, and for simplicity we write
${\ensuremath{\mathrm{ac}}}={\ensuremath{\mathrm{ac}}}_Q$.

\begin{proposition}[Monotonicity]
\label{proposition:monotonicity}
For all $X\subseteq Y\subseteq V$ it is
${\ensuremath{\mathrm{ac}}}(Y)\le{\ensuremath{\mathrm{ac}}}(X)$.
\end{proposition}

The proposition states that absorption time decreases with more absorbing
nodes. The proof is given in the Appendix.

Next we show that the absorbing random-walk centrality measure
${\ensuremath{\mathrm{ac}}}(\cdot)$ is supermodular.

\begin{proposition}[Supermodularity]
\label{proposition:supermodularity}
For all sets $X\subseteq Y\subseteq V$ and $u \not\in Y$ it is
\begin{equation}
\label{eq:supermodularity}
{\ensuremath{\mathrm{ac}}}(X) - {\ensuremath{\mathrm{ac}}}(X \cup \{u\}) \ge
{\ensuremath{\mathrm{ac}}}(Y) - {\ensuremath{\mathrm{ac}}}(Y \cup \{u\}) .
\end{equation}
\end{proposition}

\begin{IEEEproof}
Given an instantiation of a random walk, we define the following propositions
for any pair of nodes $i, j\in V$, non-negative integer $\ell$,
and set of nodes $Z$:
\begin{description}
\item[$A_{i,j}^\ell(Z)$:] The random walk started at node $i$ and
  visited node~$j$ after exactly $\ell$ steps, without visiting  any
  node in set~$Z$.
\item[$B_{i,j}^\ell(Z, u)$:] The random walk started at node $i$ and
  visited node~$j$ after exactly $\ell$ steps, having previously visited
  node~$u$ but without visiting any node in the set~$Z$.
\end{description}
It is easy to see that the set of random walks for which
$A_{i,j}^\ell(Z)$ is {{\tt true}}\ can be partitioned into
those that visited $u$ within the first $\ell$ steps and those
that did not.
Therefore, the probability that proposition $A_{i,j}^\ell(Z)$
is {{\tt true}}\ for any instantiation of a random walk generated by our model is equal to
\begin{equation}
\label{eq:split_paths_probability}
\Pr \left[ A_{i,j}^\ell(Z) \right] =
\Pr \left[A_{i,j}^\ell(Z \cup \{u\})\right] + \Pr \left[B_{i,j}^\ell(Z,u)\right].
\end{equation}
Now, let ${\ensuremath{\mathbf{\Lambda}}}(Z)$ be the number of steps for a random walk to reach the nodes in $Z$.
${\ensuremath{\mathbf{\Lambda}}}(Z)$ is a random variable and its expected value
over all random walks generated by our model is equal
to ${\ensuremath{\mathrm{ac}}}(Z)$.
Note that the proposition ${\ensuremath{\mathbf{\Lambda}}}(Z) \geq \ell + 1$
is {{\tt true}}\ for a given instantiation of a random walk only
if there is a pair of nodes $q\in Q$ and $j\in V\setminus Z$,
for which the proposition $A^\ell_{q,j}(Z)$ is {{\tt true}}.
Therefore,
\begin{equation}
\Pr \left[ {\ensuremath{\mathbf{\Lambda}}}(Z) \geq \ell + 1 \right]
= \sum_{q \in Q} \sum_{j \in V \setminus Z} \Pr \left[ A_{q,j}^\ell(Z)\right].
\end{equation}
From the above, it is easy to calculate ${\ensuremath{\mathrm{ac}}}(Z)$ as
\begin{eqnarray}
{\ensuremath{\mathrm{ac}}}(Z) & = & E[{\ensuremath{\mathbf{\Lambda}}}(Z)] \nonumber  \nonumber \\
& = & \sum_{\ell = 0}^{\infty} \ell \, \Pr \left[ {\ensuremath{\mathbf{\Lambda}}}(Z) = \ell \right]  \nonumber \\
& = & \sum_{\ell = 1}^{\infty} \Pr \left[ {\ensuremath{\mathbf{\Lambda}}}(Z) \geq \ell \right]  \nonumber \\
& = & \sum_{\ell = 0}^{\infty} \Pr \left[ {\ensuremath{\mathbf{\Lambda}}}(Z) \geq \ell + 1 \right]  \nonumber \\
& = & \sum_{\ell = 0}^{\infty} \sum_{q \in Q} \sum_{j \in V \setminus Z} \Pr \left[ A_{q,j}^\ell(Z)\right] .
\label{equation:expectation}
\end{eqnarray}
The final property we will need is the observation that,
for $X \subseteq Y$,
$B_{i,j}^\ell (Y,u)$ implies $B_{i,j}^\ell (X,u)$  and thus
\begin{equation}
\label{eq:Beta_event_probability}
\Pr \left[ B_{i,j}^\ell (X,u) \right] \ge \Pr \left[ B_{i,j}^\ell (Y,u) \right].
\end{equation}
By using Equation~(\ref{equation:expectation}),
the Inequality~(\ref{eq:supermodularity}) can be rewritten as
\begin{align}
\nonumber
\sum_{\ell=0}^{\infty} \sum_{q \in Q} & \sum_{j \in V \setminus X} \Pr
\left[ A_{q,j}^\ell(X)\right] - \\
\nonumber
& \sum_{\ell=0}^{\infty} \sum_{q \in Q} \sum_{j \in V \setminus \{X
  \cup \{u\}\}} \Pr \left[ A_{q,j}^\ell(X \cup \{u\})\right] \\
\nonumber
\ge \sum_{\ell=0}^{\infty} \sum_{q \in Q} & \sum_{j \in V \setminus Y}
\Pr \left[ A_{q,j}^\ell(Y)\right] - \\
& \sum_{\ell=0}^{\infty} \sum_{q \in Q} \sum_{j \in V \setminus \{Y \cup \{u\}\}} \Pr \left[ A_{q,j}^\ell(Y \cup \{u\})\right].
\end{align}
We only need to show that the inequality holds for an arbitrary value of $\ell$ and $q \in Q$,
that is
\begin{align}
\nonumber
&\sum_{j \in V \setminus X} \Pr \left[ A_{q,j}^\ell(X)\right]
- \sum_{j \in V \setminus \{X \cup \{u\}\}} \Pr \left[ A_{q,j}^\ell(X \cup \{u\})\right]  \ge \\
& \sum_{j \in V \setminus Y} \Pr \left[ A_{q,j}^\ell(Y)\right]
- \sum_{j \in V \setminus \{Y \cup \{u\}\}} \Pr \left[ A_{q,j}^\ell(Y \cup \{u\})\right].
\end{align}
Notice that $\Pr \left[ A_{i,u}^\ell(Y \cup \{u\}) \right] = 0$,
so we can rewrite the above inequality as
\begin{align}
\nonumber
&\sum_{j \in V \setminus X} \Pr \left[ A_{q,j}^\ell(X)\right]
- \sum_{j \in V \setminus X } \Pr \left[ A_{q,j}^\ell(X \cup \{u\})\right] \ge  \\
&\sum_{j \in V \setminus Y} \Pr \left[ A_{q,j}^\ell(Y)\right]
- \sum_{j \in V \setminus Y } \Pr \left[ A_{q,j}^\ell(Y \cup \{u\})\right].
\end{align}
To show the latter inequality we start from the left hand side and use
Inequality (\ref{eq:Beta_event_probability}).
We have
\begin{align*}
\nonumber
\sum_{j \in V \setminus X} \Pr & \left[ A_{i,j}^\ell(X)\right]
- \sum_{j \in V \setminus X } \Pr \left[ A_{i,j}^\ell(X \cup \{u\})\right] \\
= & \sum_{j \in V \setminus X} \Pr \left[ B_{i,j}^\ell(X, u)\right]  \\
\ge & \sum_{j \in V \setminus Y} \Pr \left[ B_{i,j}^\ell(Y, u)\right] \\
= & \sum_{j \in V \setminus Y} \Pr \left[ A_{i,j}^\ell(Y)\right]
- \sum_{j \in V \setminus Y } \Pr \left[ A_{i,j}^\ell(Y \cup \{u\})\right],
\end{align*}
which completes the proof.
\end{IEEEproof}

\smallskip
Finally, we establish the hardness of $k$ absorbing centrality,
defined in Problem~\ref{problem:k-ac}.

\begin{theorem}
The {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem is {\ensuremath{\mathbf{NP}}}-hard.
\end{theorem}
\begin{IEEEproof}
We obtain a reduction from the {{\sc VertexCover}} problem~\cite{GJ}.
An instance of the {\sc VertexCover}\ problem is specified by a graph
$G = (V, E)$ and an integer $k$, and asks whether there exists a
set of nodes $C\subseteq V$ such that $|C|\le k$ and $C$ is a
vertex cover,
(i.e., for every $(i,j) \in E$ it is $\{i,j\} \cap C\ne\emptyset$).
Let $|V|=n$.

Given an instance of the {{\sc VertexCover}} problem, we construct an instance of
the decision version of {{{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}} by taking the same graph $G=(V,E)$ with
query nodes $Q=V$ and asking whether there is a set of absorbing nodes
$C$ such that $|C| \le k$ and
${\ensuremath{\mathrm{ac}}}_{_Q}(C) \le 1 - \frac{k}{n}$.

We will show that $C$ is a solution for {{\sc VertexCover}}
if and only if ${\ensuremath{\mathrm{ac}}}_{_Q}(C) \le 1 - \frac{k}{n}$.

Assuming first that $C$ is a vertex cover.
Consider a random walk starting uniformly at random from a node $v\in Q=V$.
If $v\in C$ then the length of the walk will be 0,
as the walk will be absorbed immediately.
This happens with probability $|C|/|V|=k/n$.
Otherwise, if  $v\not\in C$ the length of the walk will be 1,
as the walk will be absorbed in the next step
(since $C$ is a vertex cover all the neighbors of $v$ need to
belong in $C$).
This happens with the rest of the probability $1-k/n$.
Thus, the expected length of the random walk is
\begin{equation}
{\ensuremath{\mathrm{ac}}}_{_Q}(C)
=  0\cdot \frac{k}{n} + 1\cdot\left(1 - \frac{k}{n}\right)
= 1 - \frac{k}{n} \\
\end{equation}
Conversely, assume that $C$ is not a vertex cover for $G$.
Then, there should be an uncovered edge $(u,v)$.
A random walk that starts in $u$ and then goes to $v$
(or starts in $v$ and then goes to $u$)
will have length at least 2,
and this happens with probability at least
$\frac{2}{n}\frac{1}{d_{\max}} \ge \frac{2}{n^2}$.
Then, following a similar reasoning as in the previous case, we have
\begin{eqnarray}
\nonumber
{\ensuremath{\mathrm{ac}}}_{_Q}(C) &= &
\sum_{k=0}^\infty k\, \Pr  \left( \text{absorbed in exactly } k \text{ steps}\right) \\
\nonumber
& =  & \sum_{k=1}^\infty \Pr \left( \text{absorbed after at least } k \text{ steps}\right) \\
& \ge & \left(1 - \frac{k}{n}\right) +  \frac{2}{n^2} > 1 - \frac{k}{n}.
\end{eqnarray}
\end{IEEEproof}

\section{Algorithms}
\label{section:algorithms}

This section presents algorithms to solve the {{$k$-{\sc arw}-{\sc Central}\-{\sc ity}}}\ problem.
In all cases, the set of query nodes $Q\subseteq V$ is given as input,
along with a set of candidate nodes $D\subseteq V$ and the restart probability
$\alpha$.

\subsection{Greedy approach}
\label{section:k-ac-algorithms}

The first algorithm is a standard greedy algorithm, denoted {{\sf\small Greedy}},
which exploits the
supermodularity of the absorbing random-walk centrality
measure. It starts with the result set $C$ equal to the empty set,
and iteratively adds a node from the set of candidate nodes $D$,
until $k$ nodes are added.
In each iteration the node added in the set $C$ is the one that brings
the largest improvement to ${\ensuremath{\mathrm{ac}}}_Q$.

As shown before, the objective function to be
minimized, i.e., ${\ensuremath{\mathrm{ac}}}_Q$, is supermodular and monotonically decreasing.
The {{\sf\small Greedy}}\ algorithm is not an approximation
algorithm for this minimization problem. However, it can be shown to
provide an approximation guarantee for maximizing the
{\em absorbing centrality gain} measure, defined below.
\begin{definition}[Absorbing centrality gain]
Given a graph $G$,
a set of query nodes $Q$,
and a set of candidate nodes $D$,
the absorbing centrality gain of a set of nodes
$C\subseteq D$ is defined as
\[
{\ensuremath{\mathrm{acg}}}_Q(C) = {\ensuremath{m}}_Q - {\ensuremath{\mathrm{ac}}}_Q(C),
\]
where
${\ensuremath{m}}_Q = {\min}_{v\in D}\{{\ensuremath{\mathrm{ac}}}_Q(\{v\})\}$.
\end{definition}

\para{Justification of the gain function.}
The reason to define the absorbing centrality gain is
to turn our problem into a submodular-maxi\-mization problem
so that we can apply standard approximation-theory results
and show that the greedy algorithm provides a constant-factor approximation guarantee.
The \emph{shift} ${\ensuremath{m}}_Q$ quantifies the absorbing centrality of the best
single node in the candidate set.
Thus, the value of ${\ensuremath{\mathrm{acg}}}_Q(C)$
expresses how much we gain in expected random-walk length
when we use the set $C$ as absorbing nodes
compared to when we use the best single node.
Our goal is to maximize this gain.

Observe that the gain function ${\ensuremath{\mathrm{acg}}}_Q$ is not non-negative everywhere.
Take for example any node $u$
such that ${\ensuremath{\mathrm{ac}}}_Q(\{u\}) > {\ensuremath{m}}_Q$. Then,
${\ensuremath{\mathrm{acg}}}_Q(\{u\}) < 0$. Note also that we could have obtained a
non-negative gain function
by defining gain with respect to the \emph{worst}
single node, instead of the best.
In other words, the gain function
${\ensuremath{\mathrm{acg}}}'_Q(C) = {\ensuremath{M}}_Q - {\ensuremath{\mathrm{ac}}}_Q(C)$,
with
${\ensuremath{M}}_Q = {\max}_{v\in D}\{{\ensuremath{\mathrm{ac}}}_Q(\{v\})\}$,
is non-negative everywhere.

Nevertheless, the reason we use the gain function ${\ensuremath{\mathrm{acg}}}_Q$
instead of ${\ensuremath{\mathrm{acg}}}'_Q$ is that ${\ensuremath{\mathrm{acg}}}'_Q$ takes much larger
values than ${\ensuremath{\mathrm{acg}}}_Q$, and thus, a multiplicative approximation
guarantee on ${\ensuremath{\mathrm{acg}}}'_Q$ is a
weaker result than a multiplicative approximation guarantee on ${\ensuremath{\mathrm{acg}}}_Q$.
On the other hand,
our definition of ${\ensuremath{\mathrm{acg}}}_Q$ creates a technical difficulty with
the approximation guarantee,
that is defined for non-negative functions.
Luckily, this difficulty can be overcome easily
by noting that, due to the monotonicity of ${\ensuremath{\mathrm{acg}}}_Q$,
for any $k>1$,
the optimal solution of the function ${\ensuremath{\mathrm{acg}}}_Q$,
as well as the solution returned by {{\sf\small Greedy}},
are both non-negative.

\para{Approximation guarantee.}
The fact that the {{\sf\small Greedy}}\ algorithm gives an approximation guarantee
to the problem of maximizing absorbing centrality gain is a standard
result from the theory of submodular functions.

\begin{proposition}
The function ${\ensuremath{\mathrm{acg}}}_Q$ is monotonically increasing, and
submodular.
\end{proposition}

\begin{proposition}
Let $k>1$.
For the problem of finding a set $C\subseteq D$ with $|C|\le k$,
such that ${\ensuremath{\mathrm{acg}}}_Q(C)$ is maximized, the {{\sf\small Greedy}}\ algorithm gives a
$\left(1-\frac{1}{e}\right)$-approximation guarantee.
\end{proposition}

We now discuss the complexity of the {{\sf\small Greedy}}\ algorithm.
A na\"ive implementation
requires computing the absorbing
centrality ${\ensuremath{\mathrm{ac}}}_Q(C)$ using Equation~(\ref{eq:inversion})
for each set $C$ that needs to be evaluated during the execution of
the algorithm.
However, applying Equation~(\ref{eq:inversion}) involves a matrix
inversion, which is a very expensive operation.
Furthermore, the number of times that we need to evaluate
${\ensuremath{\mathrm{ac}}}_Q(C)$ is ${\cal O}(k|D|)$, as for each iteration of
the greedy we need to evaluate the improvement over the current set of
each of the ${\cal O}(|D|)$ candidates.
The number of candidates can be very large, e.g., $|D|=n$, yielding an
${\cal O}(kn^4)$ algorithm, which is prohibitively expensive.

We can show, however, that we can execute {{\sf\small Greedy}}\ significantly more
efficiently. Specifically, we can prove the following two propositions.
\begin{proposition}
Let $C_{i-1}$ be a set of $i-1$ absorbing nodes,
$\mathbf{P}_{i-1}$ the corresponding
transition matrix, and let $\mathbf{F}_{i-1} = (\mathbf{I}-\mathbf{P}_{i-1})^{-1}$.
Let $C_i = C_{i-1} \cup \{u\}$.
Given  $\mathbf{F}_{i-1}$ the value ${\ensuremath{\mathrm{ac}}}_Q(C_i)$ can be computed in $\mathcal{O}(n^2)$.
\label{prop:fast-updates-next-step}
\end{proposition}
\begin{proposition}
Let $C$ be a set of absorbing nodes, $\mathbf{P}$ the corresponding
transition matrix, and $\mathbf{F} = (\mathbf{I}-\mathbf{P})^{-1}$.
Let $C' = C - \{v\}  \cup \{u\}$, $u, v \in C$.
Given  $\mathbf{F}$ the value ${\ensuremath{\mathrm{ac}}}_Q(C')$ can be computed in time $\mathcal{O}(n^2)$.
\label{prop:fast-updates-same-step}
\end{proposition}
The proofs of these two propositions can be found in the Appendix.
Proposition~\ref{prop:fast-updates-next-step} implies that in order to compute~${\ensuremath{\mathrm{ac}}}_Q(C_i)$
for absorbing nodes $C_i$ in $\mathcal{O}(n^2)$, it is enough to maintain the matrix $\mathbf{F}_{i-1}$,
computed in the previous step of the greedy algorithm for absorbing nodes $C_{i-1}$.
Proposition~\ref{prop:fast-updates-same-step}, on the other hand, implies that we can compute
the {{absorbing centrality}}\ of each set of absorbing nodes of a fixed size~$i$ in $\mathcal{O}(n^2)$,
given the matrix $\mathbf{F}$, which is computed for one arbitrary set of absorbing nodes $C$
of size $i$. Combined, the two propositions above yield a greedy algorithm that runs in
${\cal O}(kn^3)$ and offers the approximation guarantee discussed above. We outline it as
Algorithm~\ref{algo:greedy}.

\begin{algorithm}[t]
\caption{{{\sf\small Greedy}}}
\label{algo:greedy}
\begin{algorithmic}
\STATE {\bf Input}: graph $G$, query nodes $Q$, candidates $D$, $k\geq 1$
\STATE {\bf Output}: a set of $k$ nodes $C$
\STATE Compute ${\ensuremath{\mathrm{ac}}}_Q(\{v\})$ for arbitrary $v\in D$
\STATE For each $u\in (D - \{v\})$, use Prop.\ref{prop:fast-updates-same-step} to compute ${\ensuremath{\mathrm{ac}}}_Q({u})$
\STATE Select $u_1\in D$ s.t.\ $u_1 \leftarrow \arg\max_{u\in D} {\ensuremath{\mathrm{ac}}}_Q({u})$
\STATE Initialize solution $C \leftarrow \{u_1\}$
\FOR{$i = 2 .. k$}
\STATE For each $u\in D$, use Prop.\ref{prop:fast-updates-next-step} to compute ${\ensuremath{\mathrm{ac}}}_Q(C\cup\{u\})$
\STATE Select $u_i \in D$ s.t.\ $u_i \leftarrow \arg\max_{u_i\in (D-C)} {\ensuremath{\mathrm{ac}}}_Q(C\cup\{u\})$
\STATE Update solution $C \leftarrow C \cup \{u_i\}$
\ENDFOR
\STATE \textbf{return} $C$
\end{algorithmic}
\end{algorithm}

\para{Practical speed-up.}
We  found that the following heuristic lets us speed-up {{\sf\small Greedy}}\ even further,
with no significant loss in the quality of results.
To select the first node for the solution set $C$ (see Algorithm~\ref{algo:greedy}),
we calculate the {\it PageRank} values of all nodes in $D$ and evaluate ${\ensuremath{\mathrm{ac}}}_Q$
only for the $t << k$ nodes with highest PageRank score,
where $t$ is a fixed parameter.
In what follows, we will be using this heuristic version of {{\sf\small Greedy}}, unless explicitly stated otherwise.

\subsection{Efficient heuristics}

Even though {{\sf\small Greedy}}\ runs in polynomial time, it can
be quite inefficient when employed on moderately sized
datasets (more than some tens of thousands of nodes).
We thus describe algorithms that we study as
efficient heuristics for the problem.
These algorithms do not offer guarantee for their performance.

\spara{Spectral}
methods have been used extensively for the problem of
graph partitioning.
Motivated by the wide applicability of this family of algorithms,
here we explore three spectral algorithms:
{{\sf\small  SpectralQ}}, {{\sf\small  SpectralC}}, and
{{\sf\small  SpectralD}}.
We start by a brief overview of the spectral method;
a comprehensive presentation can be found in the
tutorial by von Luxburg~\cite{von2007tutorial}.

The main idea of spectral approaches is to project the original graph
into a low-dimensional Euclidean space so that distances between
nodes in the graph correspond to Euclidean distances between the
corresponding projected points.
A standard spectral embedding method, proposed by Shi and
Malik~\cite{Shi:2000gf},
uses the ``random-walk'' {\em Laplacian} matrix
$\mathbf{L}_G=\mathbf{I}-\mathbf{D}^{-1}\mathbf{A}$
of a graph $G$, where $\mathbf{A}$ is the adjacency matrix of the graph,
and forms the matrix
$\mathbf{U} =[u_2,\ldots,u_{d+1}]$
whose columns are the
eigenvectors of  $\mathbf{L}_G$
that correspond to the smallest eigenvalues $\lambda_2\le\ldots\le\lambda_{d+1}$,
with $d$ being the target dimension of the projection.
The spectral embedding is then defined by mapping
the $i$-th node of the graph to a point in~$\mathbb{R}^d$,
which is the $i$-row of the matrix~$\mathbf{U}$.

\iffalse
\begin{itemize}[itemsep=-1pt]
\item
Given a graph $G=(V,E)$ with adjacency matrix $\mathbf{A}$,
we compute the ``random-walk'' {\em Laplacian} matrix
$\mathbf{L}_G=\mathbf{I}-\mathbf{D}^{-1}\mathbf{A}$ of $G$, where
$\mathbf{D}$ is the diagonal  matrix whose $d_{ii}$ entry is equal to
the degree of the $i$-th node of~$G$.
\item
Given a number $\mu$, we compute the $n\times \mu$
matrix $\mathbf{U} =[u_2,\ldots,u_{\mu+1}]$ whose columns are the
eigenvectors of  $\mathbf{L}_G$
that correspond to the smallest eigenvalues
$\lambda_2\le\ldots\le\lambda_{\mu+1}$.
The eigenvector that corresponds to the smallest eigenvalue
$\lambda_1=0$ is the vector of all $1$'s and is not used.
\item
The matrix $\mathbf{U}$ is then used to define the embedding
$\phi:V\rightarrow\mathbb{R}^\mu$:
the $i$-th node of the graph is mapped to a vector in the Euclidean
space $\mathbb{R}^\mu$, which is the $i$-row of the
matrix~$\mathbf{U}$.
\end{itemize}
\fi

The algorithms we explore are adaptations of the spectral method.
They all start by computing the spectral embedding $\phi:V\rightarrow\mathbb{R}^d$,
as described above, and then, proceed as follows:

\sparanbf{{{\sf\small  SpectralQ}}}
performs $k$-means clustering on the embeddings
of the {\it query nodes}, where $k$ is the desired size of the result set.
Subsequently, it selects {\it candidate nodes} that are close to the
computed centroids.
Specifically, if $s_i$ is the size of the $i$-th cluster,
then $k_i$ candidate nodes are selected whose embedding is the nearest
to the $i$-th centroid. The number $k_i$ is selected so that
$k_i \propto s_i$ and $\sum k_i = k$.

\sparanbf{{{\sf\small  SpectralC}}}
is similar to {{\sf\small  SpectralQ}},
but it performs the  $k$-means clustering on the embeddings
of the {\it candidate nodes},
instead of the query nodes.

\sparanbf{{{\sf\small  SpectralD}}}
performs $k$-means clustering on the embeddings
of the {\it query nodes}, where $k$ is the desired result-set size.
Then, it selects the $k$ candidate nodes whose embeddings minimize
the sum of squared $\ell_2$-distances from the centroids,
with no consideration of the relative sizes of the clusters.

\spara{Personalized Pagerank} ({{\sf\small  PPR}}).
This is the standard Pagerank~\cite{PageRank} algorithm
with a damping factor equal to the restart probability $\alpha$ of the random walk
and personalization probabilities equal to the start probabilities  ${\ensuremath{\mathbf{s}}}(q)$.
Algorithm {{\sf\small  PPR}}\ returns the $k$ nodes with highest PageRank values.

\spara{Degree and distance centrality.}
Finally, we consider the standard degree and distance centrality measures.

\sparanbf{{{\sf\small  Degree}}} returns the $k$ highest-degree nodes.
Note that this baseline is oblivious to the query nodes.

\sparanbf{{{\sf\small  Distance}}} returns the $k$ nodes with highest distance
centrality with respect to $Q$.
The distance centrality
of a node $u$
is defined as
${\ensuremath{\mathrm{dc}}}(u) = \left( \sum_{v\in Q} d(u,v)\right)^{-1}$.

\section{Experimental evaluation}
\label{sec:experiments}

\subsection{Datasets}
\label{section:datasets}

We evaluate the algorithms described in Section~\ref{section:algorithms}
on two sets of real graphs:
one set of small graphs that allows us to compare the performance of
the fast heuristics against the greedy approach;
and one set of larger graphs,
to compare the performance of the heuristics against
each other on datasets of larger scale.
Note that the bottleneck of the computation lies in the evaluation
of centrality.
Even though the technique we describe in
Section~\ref{section:fasteval} allows it to scale
to datasets of tens of thousands of nodes on a single processor,
it is still prohibitively expensive for massive graphs.
Still, our experimentation allows us to discover the traits of the different
algorithms and understand what performance to anticipate when they are
employed on graphs of massive size.

The datasets are listed in Table~\ref{table:datasets}. Small graphs are obtained from
Mark Newman's repository\footnote{{http://www-personal.umich.edu/\%7Emejn/netdata/}},
larger graphs from SNAP.\footnote{\url{http://snap.stanford.edu/data/index.html}}
For {{\tt kddCoauthors}}, {{\tt livejournal}}, and {{\tt roadnet}}\
we use samples of the original datasets.
In the interest of repeatability, our code and datasets are
made publicly available.\footnote{{https://github.com/harrymvr/absorbing-centrality}}

\begin{table}[t]
\centering
\caption{Dataset statistics}
\begin{tabular}{lrr}
\toprule
Dataset   & $|V|$ & $|E|$ \\
\midrule
{{\tt karate}} & $34$ & $78$ \\
{{\tt dolphins}} & $62$ & $159$ \\
{{\tt lesmis}} & $77$ & $254$ \\
{{\tt adjnoun}} & $112$ & $425$ \\
{{\tt football}} & $115$ & $613$ \\ \hline
{{\tt kddCoauthors}} & $2\,891$   & $2\,891$\\
{{\tt livejournal}} & $3\,645$ & $4\,141$ \\
{{\tt ca-GrQc}} & $5\,242$ & $14\,496$\\
{{\tt ca-HepTh}} & $9\,877$ & $25\,998$\\
{{\tt roadnet}} & $10\,199$ & $13\,932$\\
{{\tt oregon-1}} & $11\,174$ & $23\,409$\\
\bottomrule
\end{tabular}
\label{table:datasets}
\end{table}

\subsection{Evaluation Methodology}

Each experiment in our
evaluation framework is defined by
a graph~$G$,
a set of query nodes~$Q$,
a set of candidate nodes~$D$, and
an algorithm to solve the problem.
We evaluate all algorithms presented in Section~\ref{section:algorithms}.
For the set of candidate nodes~$D$, we consider two cases:
it is equal to either the set of query nodes, i.e., $D = Q$,
or the set of all nodes, i.e., $D = V$.

Query nodes $Q$ are selected randomly,
using the following process:
First, we select a set $S$ of $\ensuremath{s}$ seed nodes, uniformly at random
among all nodes.
Then, we select a ball $B(v,r)$
of predetermined radius $r = 2$,
around each seed $v\in S$.\footnote{For the planar
{{\tt roadnet}}\ dataset we use $r = 3$.}
Finally, from all balls,
we select a set of query nodes $Q$ of predetermined size $q$,
with $q = 10$ and $q = 20$,
respectively, for the small and larger datasets.
Selection is done uniformly at random.

Finally,
the restart probability $\alpha$ is set to $\alpha=0.15$
and the starting probabilities ${\ensuremath{\mathbf{s}}}$ are uniform over $Q$.

\subsection{Implementation}
All algorithms are implemented in Python using the
NetworkX package~\cite{hagberg-2008-exploring},
and were run on an Intel Xeon 2.83GHz with 32GB RAM.

\subsection{Results}

\begin{figure*}[t]
        \centering
        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{karate_quality.png}
                \caption{{{\tt karate}}}

        \end{subfigure}%
        ~

        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{dolphins_quality.png}
                \caption{{{\tt dolphins}}}

        \end{subfigure}

        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{lesmis_quality.png}
                \caption{{{\tt lesmis}}}

        \end{subfigure}
        ~
        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{adjnoun_quality.png}
                \caption{{{\tt adjnoun}}}

        \end{subfigure}

        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{football_quality.png}
                \caption{{{\tt football}}}

        \end{subfigure}
        \caption{Results on small datasets for varying $k$ and $\ensuremath{s} = 2$.}
        \label{fig:small}
\end{figure*}

\begin{figure*}[h]
        \centering
        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{ca-GrQc_quality.png}
                \caption{{{\tt ca-GrQc}}}

        \end{subfigure}%

        ~
        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{ca-HepTh_quality.png}
                \caption{{{\tt ca-HepTh}}}

        \end{subfigure}

        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{livejournal_1000_5_quality.png}
                \caption{{{\tt livejournal}}}

        \end{subfigure}
        ~
        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{oregon1_010526_quality.png}
                \caption{{{\tt oregon-1}}}

        \end{subfigure}

        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{roadnet_ca_10000_0_quality.png}
                \caption{{{\tt roadnet}}}

        \end{subfigure}
        ~
        \begin{subfigure}[b]{0.49\textwidth}
                \includegraphics[width=\textwidth]{author-author-wpapers_KDD_all_quality.png}
                \caption{{{\tt kddCoauthors}}}

        \end{subfigure}
        \caption{Results on large datasets for varying $k$ and $\ensuremath{s} = 5$.}
        \label{fig:large_quality}
\end{figure*}

Figure~\ref{fig:small} shows the centrality scores achieved by different algorithms
on the small graphs for varying $k$ (note: lower is better).
We present two settings: on the left, the candidates are all nodes ($D = V$),
and on the right, the candidates are only the query nodes ($D = Q$).
We observe that {{\sf\small  PPR}}\ tracks well
the quality of solutions returned by {{\sf\small Greedy}},
while {{\sf\small  Degree}}\ and {{\sf\small  Distance}}\ often come close to that.
Spectral algorithms do not perform that well.

Figure~\ref{fig:large_quality} is similar to Figure~\ref{fig:small}, but
results on the larger datasets are shown, not including {{\sf\small Greedy}}.
When all nodes are candidates,
{{\sf\small  PPR}}\ typically has the best performance,
followed by {{\sf\small  Distance}},
while {{\sf\small  Degree}}\ is unreliable.
The spectral algorithms typically perform worse than {{\sf\small  PPR}}.

When only query nodes are candidates,
all algorithms demonstrate similar performance,
which is most typically worse than the performance of
{{\sf\small  PPR}}\ (the best performing algorithm) in the previous setting.
Both observations can be explained by the fact that
the selection is very restricted by the requirement $D = Q$,
and there is not much flexibility for the best
performing algorithms
to produce a better solution.

In terms of running time on the larger graphs, {{\sf\small  Distance}}\ returns
within a few minutes
(with observed times between 15 seconds to 5 minutes)
while {{\sf\small  Degree}}\ returns within seconds
(all observed times were less than 1 minute).
Finally, even though {{\sf\small Greedy}}\ returns within 1-2 seconds for the small datasets,
it does not scale well for the larger datasets
(running time is orders of magnitude worse than the heuristics and not included in the experiments).

Based on the above, we conclude that {{\sf\small  PPR}}\ offers the best trade-off of quality
versus running time for datasets of at least moderate size (more than $10\,\text{k}$ nodes).

\section{Conclusions}
\label{sec:conclusions}

In this paper, we have addressed the problem of
finding central nodes in a graph with respect to a set of query nodes~$Q$.
Our measure is based on absorbing random walks:
we seek to compute $k$ nodes that minimize
the expected number of steps that a random walk will need to
reach at (and be ``absorbed'' by) when it starts from the query nodes.
We have shown that the problem is {\ensuremath{\mathbf{NP}}}-hard and
described an $\mathcal{O}(kn^3)$ greedy algorithm to solve it approximately.
Moreover, we experimented with heuristic algorithms to solve the problem on
large graphs.
Our results show that, in practice, personalized PageRank
offers a good combination of quality and speed.

\bibliographystyle{abbrv}

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Given an undirected hypergraph $H=(V,E)$, the \emph{$k$-way hypergraph partitioning problem}
is to find disjoint subsets of its vertex set, $V_1, \dots, V_k$, called \emph{blocks}, such that the
blocks have roughly equal size and an objective function involving the cut hyperedges is minimized, e.g.\ , the sum of the weights of those hyperedges that connect multiple blocks.
The hypergraph partitioning problem has many important applications in practice such as
scientific computing or VLSI design~\cite{Papa2007}.
In particular VLSI design is a field where small improvements can lead to significant savings~\cite{Wichlund98}.
Hence, our focus in this work is on \emph{solution quality}.

It is well known that the hypergraph partitioning problem (HGP) is NP-hard~\cite{Lengauer:1990} so that mostly heuristic algorithms are used in practice.
A successful heuristic to partition large hypergraphs is the \emph{multilevel} approach~\cite{SPPGPOverviewPaper}.
Here, the hypergraph is recursively \emph{contracted} to obtain smaller hypergraphs which should reflect the same basic structure as the input. After applying an \emph{initial partitioning} algorithm to the smallest hypergraph, contraction is undone and, at each level, a
\emph{local search} method is used to improve the partitioning induced by the coarser level.
The intuition behind this approach is that a good partition at one level of the hierarchy will also be a good partition on the next finer level. Hence,
depending on the definition of the neighborhood, local search algorithms are able explore local solution spaces very effectively in this setting.
However, they are also prone to get stuck in local optima~\cite{hMetisKway}.
The multilevel paradigm helps to some extent, since
local search has a more global view on the problem on the coarse levels and a very fine-grained view on the fine levels of the multilevel hierarchy.
In addition, as with many other metaheuristics, multilevel HGP gives better
results if several repeated runs are made with some measures taken to diversify
the search.

Still even a large number of repeated executions can only scratch the surface of the huge space of possible partitionings.
In order to explore the global solution space extensively we need more sophisticated metaheuristics.
This is where memetic algorithms (MAs), i.e.~, genetic algorithms combined with local search~\cite{KimHKM11},
come into play. Memetic algorithms allow for effective exploration (global search) and exploitation (local search) of the solution space.
The general idea behind genetic algorithms is to use mechanisms inspired by biological evolution such as selection, mutation, recombination and survival of the fittest.
A genetic algorithm (GA) starts with a population of individuals (in our case partitions of the hypergraph) and evolves the population over several generational cycles~(rounds).
In each round, the GA uses a selection rule based on the fitness of the individuals of the population to select good individuals and combines them to obtain improved offspring~\cite{goldbergGA89}.
When an offspring is generated an eviction rule is used to select a member of the population to be replaced by the new offspring.
For an evolutionary algorithm it is of major importance to preserve diversity in the population~\cite{baeckEvoAlgPHD96}, i.e., the individuals should not become too similar in order to avoid a premature convergence of the~algorithm.
This is usually achieved by using mutation operations and by using eviction rules that take similarity of individuals into account.

Several genetic and memetic hypergraph partitioning algorithms have already been proposed in the literature~\cite{Areibi00anintegrated,AreibiY04,ArmstrongGAD10,BuiMoon94,Kim2004}. However
\emph{none of them} is considered to be truly competitive with state-of-the-art tools~\cite{Cohoon2003}. We believe that this is due to the fact that all of them employ
 \emph{flat} (i.e., non-multilevel) partitioning algorithms to drive the exploitation of the local solution space.

 Our \emph{main contribution} in this paper therefore is a technique that integrates a memetic algorithm with a \emph{multilevel} hypergraph partitioner.
To this end, we present sophisticated recombination and mutation operators as well as a replacement rule that uses a problem specific similarity measure.
In contrast to previous work~\cite{Areibi00anintegrated,AreibiY04,ArmstrongGAD10,BuiMoon94,Kim2004}, which only considered small and outdated~\cite{MCNCoutdated,ISPD98} ACM/SIGDA benchmark instances~\cite{MCNC} (dating back to the late 1980s), we perform extensive experiments on a large benchmark set containing hypergraphs from several application areas.
Our experiments indicate that our algorithm is able to compute partitions of very high quality and scales well to large networks.
It performs better than KaHyPar, which seems to be the current method of choice among the available hypergraph partitioning tools unless speed is more important than quality~\cite{hs2017sea},
and the state-of-the-art HGP tools hMetis~\cite{hMetisRB,hMetisKway} and PaToH~\cite{PaToH}.
In a setting where competing algorithms get the same fairly large amount of time to compute a solution,
   our new algorithm computes the best result on $648$ out of the $700$  benchmark instances.
This is in contrast to previous \emph{non-multilevel} evolutionary algorithms for the problem, which are not considered to be competitive with state-of-the-art tools~\cite{Cohoon2003}.

\section{Preliminaries}
\subsubsection*{Notation and Definitions.}
An \textit{undirected hypergraph} $H=(V,E,c,$ $\omega)$ is defined as a set of $n$ vertices $V$ and a
set of $m$ hyperedges/nets $E$ with vertex weights $c:V \rightarrow \mathbb{R}_{>0}$ and net
weights $\omega:E \rightarrow \mathbb{R}_{>0}$, where each net is a subset of the vertex set $V$ (i.e., $e \subseteq V$). The vertices of a net are called \emph{pins}.
We extend $c$ and $\omega$ to sets, i.e., $c(U) :=\sum_{v\in U} c(v)$ and $\omega(F) :=\sum_{e \in F} \omega(e)$.
A vertex $v$ is \textit{incident} to a net $e$ if $v \in e$. $\mathrm{I}(v)$ denotes the set of all incident nets of $v$.
The set $\Gamma(v) := \{ u~|~\exists~e \in E : \{v,u\} \subseteq e\}$ denotes the neighbors of $v$.
The \textit{size} $|e|$ of a net $e$ is the number of its pins.
A \emph{$k$-way partition} of a hypergraph $H$ is a partition of its vertex set into $k$ \emph{blocks} $\mathrm{\Pi} = \{V_1, \dots, V_k\}$
such that $\bigcup_{i=1}^k V_i = V$, $V_i \neq \emptyset $ for $1 \leq i \leq k$ and $V_i \cap V_j = \emptyset$ for $i \neq j$.
We use $b[v]$ to refer to the block of vertex $v$.
We call a $k$-way partition $\mathrm{\Pi}$ \emph{$\mathrm{\varepsilon}$-balanced} if each block $V_i \in \mathrm{\Pi}$ satisfies the \emph{balance constraint}:
$c(V_i) \leq L_{\max} := (1+\varepsilon)\lceil \frac{c(V)}{k} \rceil$ for some parameter $\mathrm{\varepsilon}$.
Given a $k$-way partition $\mathrm{\Pi}$, the number of pins of a net $e$ in block $V_i$ is defined as
$\mathrm{\Phi}(e,V_i) := |\{v \in V_i~|~v \in e \}|$.
For each net $e$, $\mathrm{\Lambda}(e) := \{V_i~|~ \mathrm{\Phi}(e, V_i) > 0\}$ denotes the \emph{connectivity set} of $e$.
The \emph{connectivity} of a net $e$ is the cardinality of its connectivity set: $\mathrm{\lambda}(e) := |\mathrm{\Lambda}(e)|$.
A net is called \emph{cut net} if $\mathrm{\lambda}(e) > 1$.
The \emph{$k$-way hypergraph partitioning problem} is to find an $\varepsilon$-balanced $k$-way partition $\mathrm{\Pi}$ of a hypergraph $H$ that
minimizes an objective function over the cut nets for some $\varepsilon$.
Several objective functions exist in the literature~\cite{Alpert19951,Lengauer:1990}.
The most commonly used cost functions are the \emph{cut-net} metric $\text{cut}(\mathrm{\Pi}) := \sum_{e \in E'} \omega(e)$ and the
\emph{connectivity} metric $(\mathrm{\lambda} - 1)(\mathrm{\Pi}) := \sum_{e\in E'} (\mathrm{\lambda}(e) -1)~\omega(e)$, where $E'$ is the set of all cut nets~\cite{UMPa,donath1988logic}.
Optimizing both objective functions is known to be NP-hard \cite{Lengauer:1990}.
In this paper, we use the connectivity-metric $(\mathrm{\lambda} - 1)(\mathrm{\Pi})$.
\emph{Contracting} a pair of vertices $(u, v)$ means merging $v$ into $u$.
The weight of $u$ becomes $c(u) := c(u) + c(v)$.  We connect $u$ to the former neighbors $\Gamma(v)$ of $v$ by replacing
$v$ with $u$ in all nets $e \in \mathrm{I}(v) \setminus \mathrm{I}(u)$ and remove $v$ from all nets $e \in \mathrm{I}(u) \cap \mathrm{I}(v)$.
\emph{Uncontracting} a vertex $u$ reverses the contraction.

\vspace*{-.25cm}
\subsection{Related Work}
\label{s:related}
\paragraph{Overview.}Driven by applications in VLSI design and scientific computing, HGP has evolved into a broad research area since the 1990s.
We refer to \cite{Alpert19951,DBLP:conf/dimacs/2012,Papa2007,trifunovic2006parallel} for an extensive overview.
In the following, we focus on issues closely related to the contributions of our paper.
Memetic algorithms~(MAs) were introduced in~\cite{MAs} and formalized in~\cite{RadcliffeS94} as an extension to the concept of genetic algorithms~(GAs)~\cite{Holland:1975}.
While GAs effectively explore the \emph{global} solution space, MAs additionally allow for exploitation of the \emph{local} solution space by incorporating local search
methods into the genetic framework. We refer to \cite{Moscato2010} for an introduction to memetic algorithms.
While several genetic and memetic flat (i.e., non-multilevel) hypergraph partitioning algorithms have been proposed in the literature, \emph{none of them} is considered to be truly competitive with state-of-the-art tools~\cite{Cohoon2003}.
Well-known multilevel HGP software packages with certain distinguishing characteristics include PaToH~\cite{PaToH} (originating from scientific computing),
hMetis~\cite{hMetisRB,hMetisKway} (originating from VLSI design), KaHyPar~\cite{ahss2017alenex,hs2017sea,KaHyPar-R} (general purpose, $n$-level), Mondriaan~\cite{Mondriaan} (sparse matrix
partitioning), MLPart~\cite{MLPart} (circuit partitioning), Zoltan~\cite{Zoltan}, Parkway~\cite{Parkway2.0}, and SHP~\cite{SHP}  (distributed),
UMPa~\cite{DBLP:conf/dimacs/CatalyurekDKU12} (directed hypergraph model, multi-objective), and kPaToH (multiple constraints, fixed vertices)~\cite{Aykanat:2008}.
\vfill\pagebreak
\paragraph{Evolutionary Hypergraph Partitioning.}
Saab and Rao~\cite{SaabR89} present an evolution-based approach for solving a $k$-way multi-objective, multi-constraint hypergraph partitioning problem.
Since the algorithm only works with \emph{one} individual, it does not use any recombination operators. Instead, the solution initially
generated via  bin packing is evolved using a randomized algorithm that moves vertices to different blocks if their
gain is greater than some random value. Hulin~\cite{Hulin1991} provides a GA that uses a coding scheme specifically tailored to circuit bipartitioning along with
crossover and mutation operations that respect the coding.
Bui and Moon~\cite{BuiMoon94} present a steady-state MA for ratio cut bipartitioning of hypergraphs, which uses a weak variation of the FM algorithm~\cite{FMAlgorithm}
as local search engine. To improve the performance of the crossover operation, a preprocessing step re-indexes the vertices by the visiting
order of a weighted depth first search on the clique-representation~\cite{HuMoerder85} of the hypergraph.
Areibi~\cite{Areibi00anintegrated} present a memetic algorithm that combines a GA with a modified version of Sanchis' $k$-way FM algorithm~\cite{HypergraphKFM}.
Areibi and Yang~\cite{AreibiY04} enhance
the MA presented in~\cite{Areibi00anintegrated} with a preprocessing step that clusters and contracts vertices to reduce the complexity of the hypergraphs.
Furthermore, the initial population contains both random as well as good solutions generated using the GRASP heuristic~\cite{FeoRS94}.
Armstrong et al.~\cite{ArmstrongGAD10} propose a $k$-way MA that performs crossover, mutation and local search on multiple individuals in parallel.
The traditional FM algorithm~\cite{FMAlgorithm} and Sanchis' $k$-way FM version~\cite{HypergraphKFM} are used for local search.
Kim et al.~\cite{Kim2004} present a steady-state MA for hypergraph bipartitioning, which uses a modified FM algorithm that works with lock-gains~\cite{KimM04}.
Note that \emph{none} of these algorithms makes use of the multilevel~paradigm.

\paragraph{Evolutionary Graph Partitioning.}
We refer to the survey of Kim et~al.\  \cite{KimHKM11} for a general overview and more material on genetic approaches for graph partitioning.
Soper et al. \cite{soper2004combined} provide the first algorithm that combined an evolutionary search algorithm within a multilevel graph partitioner. Here, crossover and mutation operators compute edge biases based on the input individuals.
A similar approach based on perturbations of edge weights has been used by Delling et~al.\ \cite{delling2010graph}.
Benlic et al. \cite{benlichao2010} provide a multilevel memetic algorithm for balanced graph partitioning.
PROBE \cite{ChardaireBM07} is a metaheuristic which can be viewed as a genetic algorithm without selection. It outperforms other metaheuristics, but it is restricted to the case $k=2$ and $\varepsilon=0$.
KaHIP~\cite{sandersschulz2013} contains  KaFFPaE~\cite{kaffpaE}, which has a general recombine operator framework based on a multilevel algorithm.

\subsection{$k$-way Hypergraph Partitioning using KaHyPar}
Since our memetic algorithm builds on top of the KaHyPar framework, we briefly review its core components.
While traditional multilevel HGP algorithms contract matchings or clusterings and therefore
work with a coarsening hierarchy of $\Oh{\log n}$ levels, KaHyPar instantiates the multilevel paradigm in the
extreme $n$-level version, removing only a \emph{single} vertex between two levels.
Vertex pairs $(u,v)$ to be contracted are determined using the heavy-edge rating function $r(u,v) := \sum_{e \in E'}  \omega(e)/(|e| - 1)$, where $E' := \{\mathrm{I}(u) \cap \mathrm{I}(v)\}$.
The coarsening process stops
as soon as the number of vertices drops below a certain threshold or no more contractions are possible. The framework currently contains
two coarsening algorithms.
The first algorithm~\cite{KaHyPar-R} contracts vertices in decreasing rating score order using a priority queue to store and update the ratings.
 The second algorithm~\cite{ahss2017alenex} immediately contracts each vertex with its highest-rated neighbor in \emph{random} order.
After coarsening, a portfolio of simple algorithms is used to create an initial partition of the coarsest hypergraph. During uncoarsening,
strong localized local search heuristics based on the FM algorithm~\cite{FMAlgorithm,HypergraphKFM} are used to refine
the solution by moving vertices to other blocks in the order of improvements in the optimization objective.
The framework provides a recursive bisection
algorithm to optimize the cut-net metric (KaHyPar-R~\cite{KaHyPar-R}) as well as a direct $k$-way algorithm to optimize the $(\lambda-1)$ metric (KaHyPar-K~\cite{ahss2017alenex}).
Recently,
Heuer and Schlag~\cite{hs2017sea} integrated an improved coarsening scheme into KaHyPar-K that incorporates global information
about the structure of the hypergraph into the coarsening process. It uses community detection in a preprocessing step and prevents
inter-community contractions during coarsening. This version is referred
to as KaHyPar-CA. Unless mentioned otherwise, we use the default
configurations provided by~the~authors\footnote{\url{https://github.com/SebastianSchlag/kahypar/tree/master/config}}.

\vspace*{-.35cm}
\section{Memetic Multilevel Hypergraph Partitioning}
\vspace*{-.25cm}
We now explain the components of our memetic multilevel hypergraph partitioning algorithm.
Given a hypergraph $H$ and a time limit $t$, the algorithm starts by creating an initial
population of $\mathcal{P}$ \emph{individuals}, which in our case correspond to $\varepsilon$-balanced $k$-way partitions of $H$.
The population size $|\mathcal{P}|$ is determined dynamically by first measuring the time $t_\text{I}$ spend to create one individual.
Then $\mathcal{P}$ is chosen such that the time to create $|\mathcal{P}|$ individuals is a certain percentage $\delta$ of the total running time $t$:
$|\mathcal{P}| := \max(3,\min(50,\delta\cdot(t/t_I)))$, where $\delta$ is a tuning parameter.
The lower bound on the population size is chosen to ensure a certain minimum of diversity, while
the upper bound is used to ensure convergence.
In contrast to previous approaches~\cite{Areibi00anintegrated,ArmstrongGAD10,BuiMoon94,Hulin1991,Kim2004} the population is not
filled with randomly generated individuals, but \emph{high-quality} solutions computed by KaHyPar-CA.

To judge the \emph{fitness} of an individual
we use the connectivity $(\lambda -1)(\mathrm{\Pi})$ of its partition $\mathrm{\Pi}$.
The initial population is evolved over several generational cycles using the \emph{steady-state} paradigm~\cite{EvoComp}, i.e., we generate
only \emph{one} offspring per generation.
Our two-point and multi-point recombination operators described in Section~\ref{sec:combine} improve the average quality of the population
by effectively combining different solutions to the HGP problem.
In order to sufficiently explore the global search space and to prevent premature convergence, it is important to keep
the population diverse~\cite{baeckEvoAlgPHD96}. This becomes even more relevant in our case, since with KaHyPar-CA we use powerful
heuristics to exploit the local solution space.
Previous work on evolutionary algorithms for HGP~\cite{AreibiY04,ArmstrongGAD10,BuiMoon94,Hulin1991,Kim2004} used simple mutations that
change the block of each vertex uniformly at random with a small probability. In contrast to these simple, problem agnostic operators,
we propose  mutation operators based on V-cycles~\cite{WalshawVcycle} that
exploit knowledge of the problem domain and create offspring solutions in the \emph{vicinity} of
the current population.
Furthermore in Section~\ref{sec:replace} we propose a replacement strategy which considers
fitness \emph{and} similarity to determine the individual to be evicted from the population.

\subsection{Recombination Operators}\label{sec:combine}
\vspace*{-.125cm}
The evolutionary algorithms for HGP presented in Section~\ref{s:related}
use simple multi-point crossover operators which split the parent partitions into several parts and then combine these parts to form new offspring (see Figure~\ref{fig:combine} (a)).
Since these operators do not take the structure of the hypergraph into account, offspring solutions may have
considerably worse fitness than their parents. By generalizing the recombine operator framework presented in \cite{kaffpaE} from
graphs to hypergraphs, our two-point recombine operators described in this section assure that the fitness of the offspring is
\emph{at least as good as the best of both parents}. The edge frequency based multi-point recombination operator described afterwards
gives up this property, but still generates good offspring.

\paragraph{Two-Point Recombine.}
The operator starts with selecting parents for recombination using binary tournament selection (without replacement)~\cite{BlickleT96}.
Two individuals $I_1$ and $I_2$ are chosen uniformly at random from $\mathcal{P}$ and
the individual with better fitness (i.e., lower $(\lambda-1)$ objective) becomes the first parent~$P_1$. This process is then repeated
to determine the second parent $P_2$. A tournament size of two is chosen to keep the selection
pressure low and to avoid premature convergence, since all our individuals already constitute high-quality solutions.
Both individuals/partitions are then used as input of a modified multilevel partitioning scheme as follows:
During coarsening, two vertices $u$ and $v$ are only allowed to be contracted if \emph{both parents agree
on the block assignment of both vertices}, i.e., if $b_1[u] = b_1[v] \wedge b_2[u] = b_2[v]$. This
is a generalization from multilevel evolutionary GP, i.e.~\cite{kaffpaE}, where \emph{edges running between two blocks are not eligible} for contraction
and therefore remain in the graph. In other words, our generalization allows two vertices of the same cut net to be contracted as long as the input individuals agree that they belong to the same block.
For HGP, this restriction ensures that cut nets $e$ remain in the coarsened hypergraph
\emph{and} maintain their connectivity $\lambda(e)$ regarding \emph{both} partitions. This modification is important for our optimization objective,
because it allows us to use the partition of the \emph{better} parent as initial partition of the offspring.
Since we can skip the initial partitioning phase and therefore do not need a sufficiently large number of vertices in the coarsest hypergraph to compute a good initial partition~\cite{hMetisKway}, we alter the stopping criterion of the coarsening phase such that it stops when no more contractions are possible.
The high quality solution of the coarsest hypergraph contains two different classes of vertices: Those for which both parent partitions agree on a block assignment
and those for which they don't (see Figure~\ref{fig:combine} (b) for an example).
During the uncoarsening phase, local search algorithms can then use this initial partitioning to (i) exchange good parts of the solution on the coarse
levels by moving few vertices and (ii) to find the best block assignment for those vertices, for which the parent partitions disagreed.
Since KaHyPar's refinement algorithms guarantee nondecreasing solution quality, the fitness of offspring solutions generated using
this kind of recombination is always \emph{at least as good as the better of both parents}.

\begin{figure}[t!]
    \vspace*{-.25cm}
  \centering
  \includegraphics[width=.8\textwidth]{img/combine}
  \caption{(a) Traditional, problem agnostic crossover operation to combine parent partitions $P_1$ and $P_2$ to offspring $o$.
    (b) Recombination using modified multilevel coarsening to combine two partitions (dashed red line and solid green line).
    Each cut net $e$ remains in the coarse hypergraph and maintains its connectivity $\lambda(e)$ regarding both partitions.}\label{fig:combine}
    \vspace*{-.5cm}
\end{figure}

\paragraph{Edge-Frequency Multi-Recombine.}
The operator described previously is restricted to recombine $p=2$ partitions to improved offspring of nondecreasing quality.
Sanders and Schulz~\cite{kaffpaE} specifically restrict their operators to this case, arguing that in the course of the algorithm
a series of two-point recombine operations to some extend emulates a multi-point recombine.
We here present a reasonable multi-point recombine operation to partially evaluate this hypothesis in our experimental evaluation.
Our recombine operator uses the concept of (hyper)edge frequency~\cite{Wichlund98} to pass information about the cut nets of the $t$ \emph{best individuals} in the population on to new offspring.
The frequency $f(e)$ of a net $e$ hereby refers to the number of times it appears in the cut in the $t$ best solutions: $f(e) := |\{I \in t~|~\lambda(e) > 1\}|$. We use $t = \lceil \sqrt{|\mathcal{P}|} \rceil$,
which is a common value in evolutionary algorithms~\cite{delling2010graph}.
Our multi-recombine operator then uses this information to create a \emph{new} individual in the following way. The
coarsening algorithm is modified to prefer to contract vertex pairs $(u,v)$ which share a large number of small, low-frequency nets. This is achieved
by replacing the standard heavy-edge rating function of KaHyPar with the rating function~\cite{Wichlund98} shown in Eq.~\ref{ef}:

\begin{equation}\label{ef}
r(u,v) := \frac{1}{c(v) \cdot c(u)}~\sum \limits_{e \in \{\mathrm{I}(v) \cap \mathrm{I}(u)\}}  \frac{\exp(-\gamma f(e))}{|e|}.
\end{equation}

This rating function disfavors the contraction of vertex pairs incident to cut nets with high frequency, because these nets
are likely to appear in the cut of high quality solutions. The tuning parameter $\gamma$ is used as a damping factor.
After coarsening stops, we compute an initial partition of the coarsest hypergraph using KaHyPar's initial partitioning algorithms and refine
it during the uncoarsening and local search phase.

\subsection{Mutation Operations and Diversification}\label{sec:mutate}
We define two mutation operators based on V-cycles.
All operators are applied to a random individual $I$ of the current~population.
The main idea of V-cycle based mutation operators is to reuse an already computed partition as input for the multilevel approach and to iterate
coarsening and local search phases several times using different seeds for randomization. This approach has been applied successfully
in evolutionary GP~\cite{kaffpaE}, therefore we also adopt it for HGP.
Similar to the recombine operator described in Section~\ref{sec:combine}, the quality of the solution is maintained by only contracting vertex pairs $(u,v)$
belonging to the same block ($b[u] = b[v]$). By distinguishing two possibilities for initial partitioning, we define two different
mutation operators: The first one uses the current partition of the individual as initial partition of the coarsest hypergraph and guarantees nondecreasing
solution quality. The second one employs KaHyPar's portfolio of initial partitioning algorithms to compute a new solution for the coarsest hypergraph.
During uncoarsening, local search algorithms improve the solution quality and thereby further mutate the individual.
Since the second operator computes a new initial partition which might be different from the original partition of $I$, the fitness of offspring generated by this operator can be worse than the fitness of $I$.

\subsection{Replacement Strategy}\label{sec:replace}
All recombination and mutation operators create one new offspring $o$.
In order to keep the population diverse, we evict the individual \emph{most similar} to the offspring among all individuals whose fitness is equal to or worse than $o$.
Previous work on bipartitioning~\cite{BuiMoon94,Kim2004} used the Hamming distance as a metric to measure the similarity between partitions.
We propose a more sophisticated similarity measure that takes into account the connectivity $\lambda(e)$ of each cut net $e$.
For each individual, we compute the multi-set $D :=  \{(e, m(e)) : e \in E \} $, where $m(e) := \lambda(e) - 1$ is the
multiplicity (i.e.~ number of occurrences) of $e$. Thus each cut net $e$ is represented $\lambda(e) - 1$ times in $D$.
The difference of two individuals $I_1$ and $I_2$ is the computed as $d(I_1, I_2) := | D_1 \ominus D_2 |$, where  $\ominus$ is the symmetric difference.

\section{Experimental Evaluation}\label{sec:eval}
\subsubsection*{System and Methodology.}\label{Methodology}
We implemented the memetic algorithm described in the previous section using the latest version of the KaHyPar framework.
The code is written in C++ and compiled using g++-5.2 with flags \texttt{-O3} \texttt{-mtune=native} \texttt{-march=native}.
We refer to the algorithm presented in this paper as EvoHGP.
All experiments are performed on a cluster with $512$ nodes, where each node has two Intel Xeon E5-2670 Octa-Core (Sandy Bridge) processors
clocked at $2.6$ GHz, $64$~GB main memory, $20$ MB L3- and 8x256 KB L2-Cache and runs RHEL 7.4.

We compare EvoHGP with two different configurations of KaHy-Par-CA~\cite{hs2017sea},
as well as to the $k$-way (hMetis-K) and the recursive bisection variant (hMetis-R) of hMetis 2.0 (p1)~\cite{hMetisRB,hMetisKway}, and to PaToH 3.2~\cite{PaToH}.
These HGP libraries were chosen because they provide the best solution quality~\cite{hs2017sea,ahss2017alenex}.
The first configuration of KaHyPar-CA  corresponds to the default configuration as described in \cite{hs2017sea}.
Since it is known that global search strategies are more effective than plain restarts~\cite{KahipWFCycles}, we
augment KaHyPar-CA with V-cycles (in a similar fashion as the first mutation operator) using
a maximum number of $100$ V-cycle iterations per partitioner call. This \emph{new} enhanced version of KaHyPar-CA
constitutes the second configuration and is referred to as KaHyPar-CA-V.
hMetis and PaToH are configured as described in \cite{hs2017sea}.
Since PaToH ignores the random seed if configured to use the quality preset, we only report result of the default configuration (PaToH-D).
For all five algorithms we perform repeated runs using different random seeds for each run.

To evaluate the impact of different algorithmic components of our algorithm in Section~\ref{subs:config} each EvoHGP configuration gets \emph{two} hours time
\emph{per} test instance to compute a solution. For the final evaluation in Section~\ref{subsec:eval} all algorithms get \emph{eight} hours time per test instance.
In both cases, we perform five repetitions with different seeds for each test instance
and algorithm. Due to the large amount of computing time necessary to perform these experiments, we always partition $16$ instances in
parallel on a single node.
We use the \emph{arithmetic mean} when averaging over solutions of the same instance and the \emph{geometric mean} when averaging over different instances
in order to give every instance a comparable influence on the final result.
In order to compare EvoHGP with the different algorithms, we present two kinds of plots:
\emph{Convergence plots}~\cite{kaffpaE} show the evolution of solution quality over time normalized by instance size, while \emph{performance plots}~\cite{KaHyPar-R} are used to compare
the best solutions of all algorithms on a per-instance basis.

\paragraph{Convergence Plots.}
We start by explaining how to compute the data for a single instance $I$, i.e., a $k$-way partition of a hypergraph $H$.
Whenever an algorithm computes a partition that improves the solution quality, it
reports a pair ($t$, $(\lambda - 1)$), where the timestamp $t$ is the currently elapsed time.
For $r$ repetitions with different seeds $s$, these $r$ sequences $T_s^I$ of
pairs are merged into one sequence  $T^I$ of triples $(t, s, (\lambda - 1))$, which is sorted by the timestamp $t$.
Since we are interested in the \emph{evolution} of the solution quality, we compute the sequence $T_{\text{min}}^I$ representing \emph{event-based} average values.
We start by computing the average connectivity $\overline{c}$ and the average time $\overline{t}$ using the first pair ($t$, $(\lambda - 1)$) of all $r$ sequences $T_s^I$ and insert
$(\overline{t}, \overline{c})$ into $T_{\text{min}}^I$. We then sweep through the remaining entries $(t, s, (\lambda - 1))$ of $T^I$. Each entry corresponds
to a partition computed at timestamp $t$ using seed $s$ that improved the solution quality to $(\lambda - 1)$. For each entry we therefore replace the old connectivity value of seed $s$
that took part in the computation of $\overline{c}$ with the new value $(\lambda - 1)$, recompute $\overline{c}$ and insert a new pair $(t, \overline{c})$ into $T_{\text{min}}^I$.
$T_{\text{min}}^I$ therefore represents the evolution of the average solution quality $\overline{c}$ for instance $I$ over time.
In a final step, we create the \emph{normalized} sequence $N_{\text{min}}^I$, where each entry $(t, \overline{c})$ in $T_{\text{min}}^I$ is replaced by
$(t_n, \overline{c})$ where $t_n := t/t_I$ and $t_I$ is the average time that KaHyPar-CA needs to compute a $k$-way partition of $H$.
Average values over \textit{multiple instances} are then obtained as follows:
All sequences $N^I_{\text{min}}$ of pairs $(t_n,\overline{c})$ are merged into a sequence $N_{\text{min}}$ of triples  $(t_n, \overline{c}, I)$,
which is then sorted by $t_n$. The final sequence $S_{\mathcal{G}}$ presenting event-based \emph{geometric} averages values is then computed as follows:
We start by computing the average normalized time $\overline{t_n}$ and the geometric mean connectivity $\mathcal{G}$ over all instances $I$ using the first value of all $N^I_\text{min}$ and
insert $(\overline{t_n}, \mathcal{G})$ into $S_{\mathcal{G}}$. We then sweep through the remaining entries of $N_{\text{min}}$.
For each entry $(t_n, \overline{c}, I)$, we replace the old connectivity value of $I$ that took part in the computation of $\mathcal{G}$ with the new value $\overline{c}$,
recompute $\mathcal{G}$ and insert $(t_n, \mathcal{G})$ into $S_{\mathcal{G}}$. The sequence $S_{\mathcal{G}}$ therefore represents the evolution of the solution quality
averaged over all instances~and~repetitions.

\paragraph{Performance Plots.}
These plots relate the smallest minimum connectivity of all algorithms to the corresponding connectivity produced by each algorithm on a per-instance basis.
For each algorithm, these ratios are sorted in increasing order. The plots use a cube root scale for both axes to reduce right skewness~\cite{st0223}
and show $1-(\text{best}/\text{algorithm})$ on the y-axis to highlight the instances were each partitioner performs badly.
A point close to one indicates that the partition produced by the corresponding algorithm was considerably worse than the
partition produced by the best algorithm. A value of zero therefore indicates that the corresponding algorithm produced the best solution.
Points above one correspond to infeasible solutions that violated the balance constraint.
Thus an algorithm is considered to outperform another algorithm if its corresponding ratio values are below those of the other algorithm.
In order to include instances with a cut of zero into the results, we set the corresponding cut values to \emph{one} for ratio computations.

\vspace*{-.5cm}
\subsubsection*{Benchmark Instances.} \label{Instances}
We evaluate our algorithm on a representative subset of $100$ hypergraphs from the benchmark set of Heuer and Schlag~\cite{hs2017sea}\footnote{The benchmark set was downloaded from \url{http://algo2.iti.kit.edu/schlag/sea2017/}.}, which
contains instances from four benchmark sets: the ISPD98 VLSI Circuit Benchmark Suite~\cite{ISPD98},
the DAC 2012 Routability-Driven Placement Contest~\cite{DAC2012},
the University of Florida Sparse Matrix Collection~\cite{FloridaSPM}, and the international SAT Competition 2014~\cite{SAT14Competition}.
Sparse matrices are translated into hypergraphs using the row-net model~\cite{PaToH}, i.e., each row is treated as a net and each column as a vertex.
SAT instances are converted to three different representations: For literal hypergraphs, each boolean \emph{literal} is mapped to one vertex and each clause constitutes a net~\cite{Papa2007}, while in the \emph{primal} model each variable is represented by a vertex and each clause is represented by a net.
In the \emph{dual} model the opposite is the case~\cite{FormulaPartitioning14}.
The latter two models are more common in the SAT solving community than the literal model proposed in~\cite{Papa2007}.
All hypergraphs have unit vertex and net weights.
An overview of our benchmark sets is given in Tables~\ref{tbl:instancessmall} and \ref{tbl:instanceslarge} in Appendix~\ref{app:hypergraphs}.
To compare EvoHGP with the best competing algorithms, all $100$ hypergraphs are partitioned into $k \in \{2,4,8,16,32,64,128\}$ blocks with $\varepsilon = 0.03$.
For each hypergraph $H$ and each value of $k$, a $k$-way partition $H$ is considered to be \emph{one} test instance, resulting in a total of $700$ instances.
In order to save running time, we choose a subset of 25 hypergraphs shown in Table~\ref{tbl:instancessmall}, $k=32$, and $\varepsilon = 0.03$ to evaluate the impact of different algorithmic components of our algorithm
 (recombine and mutation operations) before we run the algorithms on the large benchmark set.

\subsection{Influence of Algorithmic Components}\label{subs:config}
All configurations determine their population size $\mathcal{P}$ dynamically such that $\delta=15\%$ of the total time is spent to create the initial population.
According to the results of Wichlund and Aas~\cite{Wichlund98},
\begin{wrapfigure}{r}{8.5cm}
\centering
\vspace*{-.5cm}
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}

{\centering \includegraphics[width=.5\textwidth]{./plots/algo_conf-1}

}
\end{knitrout}
\caption{Influence of EvoHGPs algorithmic components.}\label{fig:config}
\vspace*{-.75cm}
\end{wrapfigure}
the damping factor $\gamma$ used for edge frequency calculations is set to $\gamma=0.5$.
We use a naming scheme to refer to different configurations of our algorithm. All configuration names start with EvoHGP followed by abbreviations
for the added recombine and mutation operations (multiple abbreviations are used to add multiple operations).
Abbreviation +C refers to EvoHGP using two-point recombine operations, +ER refers to EvoHGP using multi-recombine operations,
and finally +M$_{x}$ adds mutation operations with a mutation chance of $x$.
Whenever a mutation operation is performed, both operators have a 50 percent change of being chosen.
Figure~\ref{fig:config} compares different configurations of EvoHGP.
Of all configurations, EvoHGP+ER, which relies only on multi-point recombine operations,
performs worst, being only slighthly better than KaHyPar-CA-V.
Comparing its performance with EvoHGP+C (which uses only two-point recombine operations), we can see that it is indeed beneficial to guarantee nondecreasing solution quality for
combine operations. However combining both recombination operators results in a performance similar to EvoHGP+C+M$_{0.1}$. This can be explained by the fact that
multi-recombines also act as mutation operator in that they don't guarantee nondecreasing performance.
Due to the fact that the strong multilevel local search engine KaHyPar-CA computes high quality solutions, we see that a significant amount of mutations
is necessary to ensure diversity in the population. While EvoHGP+C+M$_{0.1}$ ($10\%$ mutation chance performed best for evolutionary graph partitioning in~\cite{kaffpaE}) performs
equally well as EvoHGP+C+ER, increasing the mutation rate to $50\%$ (EvoHGP+C+M$_{0.5}$) improves the overall performance of the algorithm. Moreover, we see that using both recombination operators
and mutations (EvoHGP+C+ER+M$_{0.5}$) also performs well.
Since EvoHGP+C+M$_{0.5}$ and EvoHGP+C+ER+M$_{0.5}$ show the best convergence behavior, we restrict ourselves to these configurations for the remaining experiments performed in the paper.

\vspace*{-.25cm}
\subsection{Evaluation}\label{subsec:eval}
\vspace*{-.25cm}
\begin{figure}[t!]
\centering
\begin{knitrout}
\definecolor{shadecolor}{rgb}{0.969, 0.969, 0.969}\color{fgcolor}

{\centering \includegraphics[width=\textwidth]{./plots/convergence-1}

}
\caption{Convergence plots for all instances and for different values of $k$.}\label{fig:convergence}
\end{knitrout}
\vspace*{-.5cm}
\end{figure}

\begin{table}
  \centering
\caption{Average improvement in connectivity metric over both KaHyPar-CA and the new version KaHyPar-CA-V.}~\label{tbl:imp}
\begin{tabular}{@{\extracolsep{2pt}}rcccc@{}}
  \toprule
  & \multicolumn{2}{c}{KaHyPar-CA vs. EvoHGP} & \multicolumn{2}{c}{KaHyPar-CA-V vs. EvoHGP}                 \\
  \cline{2-3}\cline{4-5}
  $k$    & +C+M$_{0.5}$                    & +C+ER+M$_{0.5}$ & +C+M$_{0.5}$ & +C+ER+M$_{0.5}$\\
  \midrule
  all & 3.3\%                           & 3.4\%           & 2.3\%        & 2.4\%           \\
2     & 0.9\%                           & 0.9\%           & 0.3\%        & 0.4\%           \\
4     & 1.3\%                           & 1.4\%           & 0.8\%        & 1.0\%           \\
8     & 2.7\%                           & 2.9\%           & 1.9\%        & 2.0\%           \\
16    & 3.5\%                           & 3.6\%           & 2.5\%        & 2.6\%           \\
32    & 4.3\%                           & 4.6\%           & 3.2\%        & 3.5\%           \\
64    & 4.9\%                           & 5.0\%           & 3.5\%        & 3.6\%           \\
128   & 5.4\%                           & 5.4\%           & 3.7\%        & 3.7\%           \\
\bottomrule
\end{tabular}
\end{table}

\begin{figure}[t]
  \centering
   \includegraphics[width=.45\textwidth]{./plots/evo_comparison_performance-1}
\includegraphics[width=.45\textwidth]{./plots/performance_all_instances-1}
\caption{Performance plots comparing the two final configurations of our evolutionary algorithm (left) and the best configuration of EvoHGP to the strongest non-evolutionary algorithms (right).
EvoHGP+C+ER+M$_{0.5}$ outperforms the non-evolutionary algorithms on almost all instances.}~\label{fig:performance}
\end{figure}

We now switch to our large benchmark set to evaluate the performance of the different algorithms under consideration.
Table~\ref{tbl:imp} and Figure~\ref{fig:convergence} compare the performance of our memetic algorithms with repeated executions of
KaHyPar-CA and KaHyPar-CA-V. When looking at convergence plots, note that KaHyPar-CA-V starts later than all other algorithms and has an initially better solution quality.
This is due to the fact it uses up to $100$ V-cycles before reporting the first solution.
The improvements of our memetic algorithms increase with increasing $k$. This is expected as the search space of possible partitionings
increases with the number of blocks.
Looking at Table~\ref{tbl:imp}, we see that both memetic algorithms on average outperform the best partitioner currently available (KaHyPar-CA), culminating in an improvement of $5.4\%$
for $k=128$. Furthermore both EvoHGP+C+M$_{0.5}$ and EvoHGP+C+ER+M$_{0.5}$ are able to improve upon the \emph{new} V-cycling version KaHyPar-CA-V for all values of $k$ and performs
$3\%$ better on average than KaHyPar-CA-V for $k \geq 32$. While the difference in solution quality between both memetic algorithms is small on average,
a Wilcoxon matched pairs signed rank test~\cite{wilcoxon} (using a 1\% significance level) reveals that the improved solution quality of EvoHGP+C+ER+M$_{0.5}$
is statistically significant ($Z=-2.992857$, $p=0.002763795$). This is also confirmed by the performance plot shown in Figure~\ref{fig:performance} (left).

Looking at the performance plot that compares the strongest non-evolutionary algorithms with the strongest memetic configuration in Figure~\ref{fig:performance} (right),
we see that EvoHGP+C+ER+M$_{0.5}$ performs significantly better than \emph{all} other algorithms. It produces the best partitions for $648$ of the $700$ instances.
It is followed by KaHyPar-CA-V (107), KaHyPar-CA (79), hMetis-R (55), hMetis-K (46), and PaToH-D (35).
Note that for some instances, multiple partitioners computed the same solution.
Comparing the best solutions of EvoHGP+C+ER+M$_{0.5}$ to each partitioner individually, it produced better partitions than PaToH-D, hMetis-K, hMetis-R, KaHyPar-CA,
KaHyPar-CA-V in $661$, $644$, $636$, $609$, $585$ cases, respectively.

This shows that even a large number of repeated executions helps only partially to explore the huge space of possible partitionings.
By combining effective exploration (global search) with exploitation (in our case using powerful $n$-level HGP algorithms)
our memetic algorithm can effectively help to break out of local minima and hence explore the global solution space more extensively.

\vspace*{-.25cm}
\section{Conclusion and Future Work}
\vspace*{-.25cm}
EvoHGP is the first \emph{multilevel} memetic algorithm to tackle the balanced hypergraph partitioning problem.
Key components of our contribution are new effective multilevel recombine and mutation operations that incorporate information about the best solutions in the coarsening process
and provide a large amount of diversity.
Experiments comparing EvoHGP with a V-cycling version of KaHyPar-CA and the well known HGP tools hMetis and PaToH indicate that our evolutionary algorithm computes by far the best solutions on almost all instances.
This confirms our conjecture that previous attempts to solve the HGP problem using memetic algorithms failed to be competitive with
state-of-the-art tools because (i) only flat partitioning algorithms were used to drive the exploitation phase and (ii) evolutionary
operators were problem agnostic and thus did not incorporate enough structural information into the algorithm.
We therefore believe that EvoHGP is helpful in a wide area of application areas in which solution quality is of major importance.
In the future, it would be interesting to apply EvoHGP in such application areas and to try other domain specific recombine operators that offer more specific knowledge of the application domain.
In addition, it may be worth to investigate shared-memory parallelization as in~\cite{ArmstrongGAD10} or a distributed memory parallelization based on islands as in~\cite{kaffpaE}.
Lastly, we plan to integrate our algorithm in the KaHyPar framework.

\vfill\pagebreak