A QUBO Formulation for Eigencentrality

Prosper D. Akrobotu1,2, Tamsin E. James2, Christian F. A. Negre3, Susan M.
Mniszewski2*

1 Department of Mathematical Sciences, The University of Texas at Dallas,
Richardson, TX, USA

2 Computer, Computational, and Statistical Sciences Division, Los Alamos
National Laboratory, Los Alamos, NM, USA

3 Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM, USA

These authors contributed equally to this work.

*<EMAIL_ADDRESS>
## Abstract

The efficient calculation of the centrality or “hierarchy” of nodes in a
network has gained great relevance in recent years due to the generation of
large amounts of data. The eigenvector centrality is quickly becoming a good
metric for centrality due to both its simplicity and fidelity. In this work we
lay the foundations for the calculation of eigenvector centrality using
quantum computational paradigms such as quantum annealing and gate-based
quantum computing. The problem is reformulated as a quadratic unconstrained
binary optimization (QUBO) that can be solved on both quantum architectures.
The results focus on correctly identifying a given number of the most
important nodes in numerous networks given by our QUBO formulation of
eigenvector centrality on both the D-Wave and IBM quantum computers.

## Introduction

There are several centrality measures used to identify the most influential
node or nodes within a network, each having their own benefits dependent on
the data at hand or the results desired.

For example, degree centrality [1], which is based purely on the number of
connections a node has, could be used for identifying the most popular person
within a group of people on a social media platform (number of followers).
Closeness centrality [2] is dependent on the length of the paths from one node
to all other nodes in a network, prioritizing nodes that are “closer” to all
other nodes as more central. This has been used for predicting enzyme
catalytic residues from topological descriptions of protein structures [3].
Betweenness centrality [4], is based on the number of times a node appears
when two other nodes are connected by their shortest path. This measure is
often used in biological networks, for example identifying a specific protein
that is important for information flow within a network, which could be used
in drug discovery [5]. Katz centrality [6], measures the importance of a node
through its immediate connections, and also the connections of other nodes
through the immediate neighbors. Katz centrality has also been used within a
biological setting, such as identifying disease genes [7]. PageRank centrality
[8] is a variant of eigenvector centrality designed to rank web pages by
importance based on links between pages or articles. This differs from
eigenvector centrality as it takes into account directions between nodes
(clicking from one web page to another).

It is worth noting at this point that the rankings generated by these
different centrality measures are correlated, especially for the most highly
ranked nodes [9, 10, 11, 12, 13]. Our research is concentrated on a study of
the eigenvector centrality (EC) measure [14], and using this to determine the
most important nodes in a given network. EC has been applied in many different
fields of science. For example, identifying the most important amino acid
residues in proteins undergoing an allosteric mechanism [15] and predicting
flow-paths in porous materials [16]. It is also relevant to current world
issues related to the COVID-19 pandemic, such as identifying people deemed as
“super-spreaders” and areas that are hot-spots in a pandemic (network of
people) [17], and also in the analysis of production chains in the financial
market, where micro-sectors are identified as important nodes within a chain
[18].

EC is a centrality measure which assigns to each node a value that is
proportional to the sum of the values for the node’s neighbors [19, 15]. With
this measure, the most influential node is a node that is connected to a
majority of the other important nodes in the network. The scheme designed to
implement the EC measure starts by obtaining the network’s adjacency matrix, a
matrix that describes the network’s connectivity and is defined as follows:
Let $A=[a_{ij}]_{n\times n}$ be the adjacency matrix of a network or graph
$G=(V,E)$ of $n$ nodes with entries defined by

$a_{ij}=\begin{cases}1&\\{u_{i},u_{j}\\}\in E\text{ for }u_{i},u_{j}\in V\\\
0&\text{otherwise}\end{cases}$

where $E$ is the set of edges and $V$ is the set of nodes or vertices. The
Perron-Frobenious theorem[20] states that there is a unique largest real
eigenvalue for a non-negative square matrix (here given by the adjacency
matrix), with an eigenvector solution consisting of positive elements. So let
the largest eigenvalue of $A$ be $\lambda_{1}$. Then the EC measure assigns to
each node the value

$\displaystyle\vspace{-15pt}x_{i}$
$\displaystyle=\frac{1}{\lambda_{1}}\sum_{j=1}^{n}a_{ij}x_{j}\hskip
5.0pt\,\forall\,i=1,2,\ldots,n$ (1)

where $x_{i}$ is the given centrality value for the ith node. It is
represented in matrix form by

$\displaystyle A\boldsymbol{x}=\lambda_{1}\boldsymbol{x}$ (2)

(a) Degree centrality

(b) Eigenvector centrality

Fig 1: Degree and Eigenvector Centrality of a Graph. Identifying the
importance of a node in a network based on (a) degree centrality and (b)
eigenvector centrality. The color and radius of each disk around a node is
dependent upon the centrality values. The least central nodes are colored in
purple, the mid-central nodes are colored in green with the most central nodes
colored yellow. Nodes $1,2$ and $3$ have the highest values when using the
degree centrality measure while node $0$ has the highest value when using the
EC measure.

The degree centrality can also be defined as the count of the number of walks
of length one that reach the node for which centrality is being computed. EC
instead, is a count of the number of walks of infinite length [14, 9]. A brief
demonstration of these concepts are shown in the Appendix with details also
found in [10].

Due to the number of uses of centrality measures, it appears to be a good step
forward to reformulate these problems for use on quantum computers. As
networks become larger and more complicated, such as in drug discovery, the
use of quantum computers could prove to be of great help in the future. Our
research seeks to reformulate the iterative scheme for the EC problem as a
classical optimization problem and encode it as a quadratic unconstrained
binary optimization (QUBO) problem for quantum computers.

The steady progress in the field of quantum computing since its proposal in
the early 1980s by Richard Feynman has seen researchers trying different
directions to circumvent the complexity of constructing portable physical
quantum computers. Currently, there are two major approaches for building
quantum computers: gate-based and quantum annealing. The gate-based quantum
computers are designed using quantum circuits with control and manipulative
power over the evolution of quantum states to tackle general problems arising
in nature [21, 22]. The quantum annealing approach, however, uses the natural
evolution of quantum states to tackle specific problems such as probabilistic
sampling and combinatorial optimization problems [22].

The D-Wave quantum computer is a quantum computing platform that uses quantum
annealing (QA), a heuristic search method that makes use of quantum tunneling
and quantum entanglement to solve the ground state of the Ising model
equivalence of combinatorial optimization problems. That is, any problem to be
solved by the D-Wave quantum annealer, is to be modeled as a search for the
minimum energy of the Ising Hamiltonian energy function as follows.

$\displaystyle
E(\boldsymbol{s})=\sum_{i}h_{i}s_{i}+\sum_{i<j}J_{i,j}s_{i}s_{j}$ (3)

where $s_{i}\in\\{-1,1\\}$ are magnetic spin variables subject to local fields
$h_{i}$ and nearest neighbor interactions with coupling strength $J_{ij}$ or
its Boolean equivalence obtained from the transformation
$\boldsymbol{s}=2\boldsymbol{x}-\boldsymbol{1}$ where the entries of the
vector $\boldsymbol{x}$ represent the binary variables $x_{i}\in\\{0,1\\}$ and
$\boldsymbol{1}$ is a vector of ones. The Boolean equivalence of the Ising
problem is referred to as a QUBO problem, with the following equation

$\displaystyle
E(\boldsymbol{x})=\boldsymbol{x}^{T}Q\boldsymbol{x}=\sum_{i=1}^{N}x_{i}Q_{ii}+\sum_{i\neq
j}x_{i}x_{j}Q_{ij}.$ (4)

D-Wave quantum annealers include the current 2000Q and the new Advantage [23].
The quantum processing unit (QPU) of the D-Wave 2000Q has up to 2048 qubits
and 6061 couplers sparsely connected as a Chimera graph $C_{16}$ [23]. While
the Advantage consists of more than 5000 superconducting qubits connected with
35,000 couplers on a Pegasus graph [23]. The sparse connectivity of the
chimera graph of the D-Wave 2000Q requires a “minor embedding” of the Ising
model connectivity of the problem onto the hardware. This results in chains of
physical qubits representing logical qubits leading to a maximum capacity of
64 fully connected logical qubits/variables [24]. The D-Wave quantum computers
possess the ability to sample degenerate ground state solutions and have been
utilized in solving several problems such as quantum isomer search [25], graph
partitioning [26], community detection [27], binary clustering [28], graph
isomorphism [29] and machine learning [30, 31]. It has also been used in
solving physical problems related to atomistic configuration stability [32],
job-shop scheduling [33], airport and air traffic management [34, 35].

The gate-based quantum computers use unitary operations defined on a quantum
circuit to transform input data into a desired output data [21, 36]. This
computational mechanism is employed in the design of IBM quantum computers
which are made available through a cloud-based platform called IBM Quantum
(IBM-Q) Experience [37]. The unitary operations are designed to process data
with high fidelity and to tackle both combinatorial optimization problems and
non-combinatorial problems like prime factorization [21]. IBM-Q consists of
both quantum hardware and simulators. The basic steps to follow in carrying
out any experiment on the IBM-Q Experience requires first to specify a quantum
circuit via a graphic interface called composer or a text-based editor (the
cloud version is called quantum lab), then run the circuit on a simulator to
verify specifications, and finally execute the circuit on the quantum
processor for a number $N$ of shots with $N=8192$ being the maximum allowed on
current devices [36].

In this paper, our main goal is to reformulate the EC problem as a QUBO
problem that correctly identifies a desired number of the most important nodes
of an undirected graph when solved on quantum annealing and gate-based quantum
computers. As a secondary objective, we probe the possibility of using our
formulation to generate a ranking of the nodes that correspond to the rankings
obtained via the usual iterative method of solving the EC problem.

The paper is organized as follows. Section 2 elaborates on our extraction of
an optimization problem out of the EC problem classically. We then construct a
QUBO formulation for identifying a number of the most important nodes of the
graph by the EC measure to be implemented on quantum computing devices.
Section 3 presents the software tools, implementations and results obtained
from the classical and quantum computers such as D-Wave 2000Q and IBM-Q. It
also stipulates ways of defining a rank from the results obtained from solving
the QUBO on D-Wave 2000Q and provides discussions on the conclusions derived
from the results. In section 4, we provide a summary of the results and
conclude with suggestions of future directions to be considered.

## Methods

In this section we present the mathematical formulation of the EC problem as
an optimization problem. A careful examination of the scheme in Eq. (1) and
Eq. (2) shows that the EC is a problem of determining the eigenvector
corresponding to the leading eigenvalue of the adjacency matrix of the
network. To this end, we recall some useful properties for the leading
eigenvalue of a symmetric matrix.

  * •

$\frac{x^{T}Ax}{x^{T}x}\leq\lambda_{1}\leq d_{max}$, where $d_{max}$ is the
maximum degree of the graph.

Hence it is understood that the maximum of the set of numbers
$\\{\frac{x^{T}Ax}{x^{T}x}\\}$ coincides with the leading eigenvalue of the
adjacency matrix and thus one can construct a maximization problem from this
property as follows:

$\displaystyle\max_{\boldsymbol{x}\,\in\,\mathbb{R}^{n}}\boldsymbol{x}^{T}A\boldsymbol{x}\hskip
10.0pt\mathrm{s.t.}\,\hskip 5.0pt\parallel\boldsymbol{x}\parallel=1\hskip
5.0pt$ (5)

This maximization problem is a constrained optimization problem which is
equivalent to the following unconstrained minimization problem:

$\displaystyle\min_{\boldsymbol{x}\,\in\,\mathbb{R}^{n}}\left[-\boldsymbol{x}^{T}A\boldsymbol{x}+P\left(\sum_{i=1}^{n}x_{i}^{2}-1\right)^{2}\right]$
(6)

where $P$ is the Lagrange multiplier or simply a penalty constant and we have
used the fact that
$1=\parallel\boldsymbol{x}\parallel^{2}=\sum_{i=1}^{n}x_{i}^{2}$. The goal
here is that the argument of a typical solution to the unconstrained
minimization problem should preserve the ranking on the nodes in the network
when the usual iterative scheme Eq. (2) is used in determining the importance
of a node. That is, we are more interested in the rank assigned than the
centrality values assigned to each node in the network.

### Eigenvector centrality as a QUBO problem

Recall that EC measures the influence of a node in a network. We therefore
seek a QUBO problem whose solution classifies the most influential nodes in a
network when the EC measure is used. That is to say, the proposed QUBO should
directly determine the most influential nodes that would have been identified
from the eigvenvector in Eq. (1).

The word “binary” suggests that our search domain must be a binary field,
however the original problem requires a search space of real vector spaces.
Thus to efficiently transition into binary fields, we will consider simply
splitting the set of nodes into two categories: most central and least
central, where a value of $1$ denotes a node is most central and a value of
$0$ denotes a node is least central. To this end, we define $\tau$ as the
number of most central nodes we wish to be identified from the set of nodes
(i.e. how many nodes are to be assigned the value of 1). The problem of
splitting into two categories is binary as we only have two categories: 0 and
1, or high and low. The value $\tau$ therefore must be chosen with
consideration of factors such as the size of the network and how many
important nodes you wish to identify. This definition will require a slight
modification to our unconstrained minimization model, Eq. (6) and a need for
more constraints. Now, consider the second term in Eq. (6), since the search
field is now binary, we have for each $i$, $x_{i}^{2}=x_{i}$ and based on the
definition of $\tau$, $\sum_{i=1}^{n}x_{i}=\tau$, we adapt the following
modification to the penalty term.

$\displaystyle P\left(\sum_{i=1}^{n}x_{i}^{2}-1\right)^{2}\mapsto
P\left(\sum_{i=1}^{n}x_{i}-\tau\right)^{2}$ (7)

We now write the modified penalty term in matrix notation.

$\displaystyle\left(\sum_{i=1}^{n}x_{i}-\tau\right)^{2}$
$\displaystyle=\left(\sum_{i=1}^{n}x_{i}\right)^{2}-2\tau\sum_{i=1}^{n}x_{i}+\tau^{2}$
$\displaystyle=\sum_{i=1}x_{i}^{2}+2\sum_{i<j}x_{i}x_{j}-2\tau\sum_{i=1}^{n}x_{i}+\tau^{2}$
$\displaystyle=(1-2\tau)\sum_{i=1}^{n}x_{i}^{2}+2\sum_{i<j}x_{i}x_{j}+\tau^{2}$
$\displaystyle=\boldsymbol{x}^{T}\left[(1-2\tau)I+U\right]\boldsymbol{x}+\tau^{2}$
$\displaystyle=\boldsymbol{x}^{T}C\boldsymbol{x}+\tau^{2}$ (8)

where we have used the fact that $x_{i}=x_{i}^{2}$ and $C=(1-2\tau)I+U$, I is
the $n\times n$ identity matrix and $U=[u_{ij}]$ is an $n\times n$ matrix with
entries $u_{ij}$ defined by

$\displaystyle u_{ij}$ $\displaystyle=\begin{cases}1&i\neq j\\\
0&\text{otherwise}\end{cases}$ (9)

We now attempt to build a problem Hamiltonian from the eigenvector equation,
Eq. (2). Note that

$\displaystyle\left(A\boldsymbol{x}-\lambda_{1}\boldsymbol{x}\right)^{2}$
$\displaystyle=[(A-\lambda_{1}I)\boldsymbol{x}]^{T}[(A-\lambda_{1}I)\boldsymbol{x}]$
$\displaystyle=\boldsymbol{x}^{T}(A-\lambda_{1}I)^{T}(A-\lambda_{1}I)\boldsymbol{x}$
$\displaystyle=\boldsymbol{x}^{T}A^{T}A\boldsymbol{x}-2\lambda_{1}\boldsymbol{x}^{T}A\boldsymbol{x}+\lambda_{1}^{2}\boldsymbol{x}^{T}\boldsymbol{x}$
$\displaystyle=\boldsymbol{x}^{T}A^{T}A\boldsymbol{x}-2d_{max}\boldsymbol{x}^{T}A\boldsymbol{x}+d_{max}^{2}\boldsymbol{x}^{T}\boldsymbol{x}+\mathrm{error}$
$\displaystyle=\boldsymbol{x}^{T}A^{2}\boldsymbol{x}-2d_{max}\boldsymbol{x}^{T}A\boldsymbol{x}+d_{max}^{2}\boldsymbol{x}^{T}\boldsymbol{x}+\mathrm{error}$
(10)

where we have used $A^{T}=A$ for undirected graphs. We wish to obtain the
ground state of Eq. (Eigenvector centrality as a QUBO problem) and determine
whether there is any meaningful information leading to the identification of
the most influential node of the graph. Our investigations showed no
conclusive information embedded in the ground state for identifying the most
central node. However, conclusive information could be draw from the first
excited state. This then suggested a need to modify the above objective
function Eq. (Eigenvector centrality as a QUBO problem). Therefore motivated
by the fact that EC is also a measure of walks of infinite length, we proposed
a modified objective function whose symmetric matrix is defined by Eq. (11).
This form was inspired by the form of the quantum search Hamiltonian of Childs
and Goldstone’s continuous time quantum walk algorithm described in [38].

$\displaystyle
Q=-P_{0}A^{2}\boldsymbol{\hat{d}}\boldsymbol{\hat{d}}^{T}A-P_{0}A\boldsymbol{\hat{d}}\boldsymbol{\hat{d}}^{T}A^{2}+P_{1}C$
(11)

where $P_{0},P_{1}$ are penalty constants such that $P_{1}>P_{0}$, and

$\displaystyle\boldsymbol{d}$
$\displaystyle=\sum_{i}^{n}d_{i}\boldsymbol{e_{i}}\hskip
5.0pt\mathrm{where}\hskip
5.0ptd_{i}=\sum_{j=1}^{n}\boldsymbol{e_{i}}^{T}A\boldsymbol{e_{j}},\hskip
10.0pt\mathrm{and}\hskip
5.0pt\boldsymbol{\hat{d}}=\frac{\boldsymbol{d}}{\parallel\boldsymbol{d}\parallel}$
(12)

Here the vectors $\boldsymbol{e_{i}}\in\mathbb{R}^{n}$ are the canonical basis
vectors of $\mathbb{R}^{n}$. We minimize over all the state vectors of the
nodes, $\boldsymbol{x}\in\\{0,1\\}^{n}$, of the graph $G$ and verify using
experimental results on several graphs the proposition:

###### Proposition 1

Let $G=(V,E)$ be an undirected graph with $n=|V|$ nodes, adjacency matrix $A$
and degree sequence
$\boldsymbol{d}=(d_{1},d_{2},\ldots,d_{n})\in\mathbb{R}^{n}$. Using the EC
measure, the most central node of the graph $G$ is the ground state of the
QUBO problem.

$\displaystyle\min_{\boldsymbol{x}\in\\{0,1\\}^{n}}\boldsymbol{x}^{T}Q\boldsymbol{x}$
$\displaystyle Q$
$\displaystyle=-P_{0}A^{2}\boldsymbol{\hat{d}}\boldsymbol{\hat{d}}^{T}A-P_{0}A\boldsymbol{\hat{d}}\boldsymbol{\hat{d}}^{T}A^{2}+P_{1}C$
$\displaystyle C$ $\displaystyle=(1-2\tau)I+U$

with $\tau=1$, where the matrix $U$ is defined in Eq. (9).

## Results

### Tools and Implementation

Using Python packages NumPy [39], SciPy [40], D-Wave Ocean [41], IBM Qiskit
[42], NetworkX [43], and Matplotlib [44] and the D-Wave 2000Q and IBM Q
hardware, we were able to verify proposition 1 (with $\tau\geq 1$) for
different graphs (see S1 Fig. and S2 Fig.) by comparing our D-Wave 2000Q/IBM-Q
output with centrality rankings obtained from the NetworkX EC tool. The QUBO
constructed for the EC problem was implemented on the D-Wave 2000Q_LANL
machine at Los Alamos National Laboratory [45] and also on IBM Q using IBM
Qiskit QASM simulator or real quantum devices available on the IBM Q, in
particular ibmq_manhattan [42]. At the front-end of the D-Wave platform, we
use D-Wave Ocean tools to submit instructions for the optimal ground state
solution for the problem Hamiltonian/QUBO with specified parameters such as
the anneal time, chain strength, post-processing method and the number of
samples to be collected. The front-end then sends the instructions to the
2000Q_LANL solver chip for processing. Once the problem Hamiltonian is
successfully embedded unto the chip, the annealer solves the QUBO for the
minimum energy solution which is a bit string that minimally violates the
constraint. Note that the bit string returned has $\tau$ number of $1$’s whose
corresponding index denote the the top $\tau$ influential nodes of the graph.
To solve the QUBO problem on the IBM-Q Experience platform, we employed
Qiskit’s CPLEX tools [46] in generating a quadratic program that is converted
to a QUBO/Ising operator for building quantum instances on available QASM
simulators or real quantum devices available on the IBM Q. The ground state of
the QUBO Hamiltonian is then solved using a Minimum Eigen Solver [47] such as
quantum approximate optimization algorithm (QAOA) [48]. Due to qubit
limitation, the QUBO can be implemented for graphs with at most $65$ nodes on
IBM-Q’s Manhattan which can encode at most $65$ qubits.

### Results

Our investigation considered fabricated and famous graphs such as shown in S1
Fig. and S2 Fig.. These graphs were created using the NetworkX graph generator
algorithms [49]. The graphs in S2 Fig. include scale-free networks, that is
networks with power law or scale-free degree distribution. The Barabasi-Albert
graph, $G_{4}$ (see S2 Fig. d), is an example of a scale-free network which
integrates two essential concepts in real networks: growth (increasing number
of nodes in the network) and preferential attachment (highly connected nodes
have a maximum likelihood of obtaining new connections) [50]. Other well known
graphs considered are social networks such as the Davis Southern women social
network (see S2 Fig. b) – the network of a Southern women social club made up
of $18$ women who attended $14$ different events, and the Karate Club graph
(see S2 Fig. f) – network of a university karate club, Bull graph, $G_{10}$
(see S1 Fig. d), complete graphs $K_{n}$ (e.g., see S1 Fig. f), complete
bipartite graphs $K_{m,n}$ (e.g., see S1 Fig. e), Sedgewick-Maze graph (see S1
Fig. a), Barbell graph – two complete graphs joined together by a path graph,
(e.g., the Lollipop graph, $G_{7}$ (see S2 Fig. f)), Tutte graph $G_{3}$ (see
S2 Fig. c) – a cubic polyhedral graph with $46$ nodes and 69 edges. The
fabricated graphs were mostly tree graphs– connected acyclic undirected graphs
such as graphs $G_{1}$ (see S2 Fig. a), and $G_{8}$ (see S2 Fig. b). The tree
graphs $G_{1}$ and $G_{8}$ were fabricated mainly to test and show that the
QUBO formulation is correctly executing EC rather than a degree centrality
measure.

Table 1: Table of Networkx and QUBO results with $\tau=1$ and $\tau=5$.

Graph $G=(V,E)$ Features | Most central node by EC | Top 5 most central nodes
---|---|---
Name |  $|V|$ |  $|E|$ |  Density $\rho(G)$ |  Nx |  QA |  Nx |  QA
Bull |  5 |  5 |  0.5 |  1 |  1 |  0,1,2,3,4 |  0,1,2,3,4
$K_{3,5}$ |  8 |  15 |  0.54 |  0,1,2 |  2 |  0,1,2,3,4 |  0,1,2,3,4
$K_{5}$ |  5 |  10 |  1 |  - |  0 |  - |  -
$K_{6}+S_{6}$ |  12 |  18 |  0.32 |  5 |  5 |  0,1,2,3,4,5 |  0,1,2,3,5
Florentine Family $F$ (S1 Fig. h) |  15 |  20 |  0.19 |  Medici |  Medici |  Guadagni, Medici, Ridolfi, Strozzi, Tomabuoni |  Guadagni, Medici, Ridolfi, Strozzi, Tomabuoni
Tutte $G_{3}$ |  46 |  69 |  0.07 |  - |  7 |  - |  12, 30, 36, 44, 45
Lollipop $L_{10,20}=G_{7}$ |  30 |  65 |  0.15 |  9 |  9 |  0,1,2,4,9 |  0,1,2,4,9
Barabasi-Albert $BA(50,5)=G_{4}$ |  50 |  225 |  0.18 |  7 |  7 |  5,6,7,8,11 |
Davis southern women $G_{2}$ |  32 |  89 |  0.18 |  E8 |  E8 |  E7,E8,E9, Evelyn Jefferson, Theresa Anderson |  E7,E8,E9, Evelyn Jefferson, Theresa Anderson
Karate club |  34 |  78 |  0.14 |  33 |  33 |  0,1,2,32,33 |  0,1,2,32,33
Sedgewick-Maze |  8 |  10 |  0.36 |  4 |  4 |  0,3,4,5,7 |  0,3,4,5,7
$G_{8}$ |  16 |  15 |  0.13 |  0 |  0 |  0,1,2,3,4 |  0,1,2,3,8
$G_{1}$ |  20 |  19 |  0.1 |  1 |  1 |  0,1,2,3,5 |  0,1,2,5,6
$M$ (S1 Fig. i) |  13 |  15 |  0.19 |  2 |  2 |  0,1,2,6,9 |  0,1,2,6,9

The first 4 columns describe some basic graph features such as the name of the
graph, number of nodes $|V|$, number of edges $|E|$ and the density of the
graph $\rho=\frac{2|E|}{|V|(|V|-1)}$. Columns 5 and 6 describe the most
central node corresponding to the EC NetworkX result and D-Wave QA result for
the QUBO. The last two columns describe the top 5 most central nodes using the
EC measure in NetworkX and by implementing the QUBO, $Q$ on D-Wave (QA). A“-”
signifies no node was identified as most central or least central, all nodes
have the same centrality.

NetworkX was used to obtain initial EC measures and node rankings for
comparison with results from the quantum computations on small graphs.

### Discussion

We first analyzed Eq. (6), using classical optimization tools from SciPy [40].
The analysis is aimed at verifying the possibility of constructing an
unconstrained optimization problem whose solution determines the EC values of
a graph exactly or within a small margin of error and concurrently preserves
the node rankings of the graphs before we went on to formulate the quantum
version.

Classical minimization using this equation preserved the ranking of the nodes
within a small error (see Table 2). This is of small concern however as these
nodes have extremely similar values. Therefore, the results obtained using
classical minimization in Eq. (6) were in support of our hypothesis with the
worst case showing the invariance of the rankings in the top 25% of the most
central nodes.

Table 2: NetworkX EC values and node rankings for graph $G_{8}$.

Rank | Node | NetworkX | Node | Scipy NM
---|---|---|---|---
1 | 0 | 0.4629 | 0 | 0.4625
2 | 1 | 0.4082 | 3 | 0.4116
3 | 2 | 0.4082 | 2 | 0.4075
4 | 3 | 0.4082 | 1 | 0.4068
5 | 4 | 0.1543 | 12 | 0.1580
6 | 5 | 0.1543 | 15 | 0.1563
7 | 6 | 0.1543 | 14 | 0.1556
8 | 7 | 0.1543 | 4 | 0.1552
9 | 8 | 0.1543 | 8 | 0.1546
10 | 9 | 0.1543 | 10 | 0.1545
11 | 10 | 0.1543 | 9 | 0.1539
12 | 11 | 0.1543 | 13 | 0.1538
13 | 12 | 0.1543 | 5 | 0.1536
14 | 13 | 0.1543 | 6 | 0.1522
15 | 14 | 0.1543 | 7 | 0.1513
16 | 15 | 0.1543 | 11 | 0.1504

The result was obtained using the SciPy minimizer with the Nelder-Mead method
in solving Eq.(6). The column titled “Rank”, defines the rank on each node in
the columns titled “Nodes” using the centrality values obtained from both
methods.

(a) NetworkX ranking

(b) Ranking obtained from IBM-Q and D-Wave for $\tau=1$

(c) Ranking obtained on both IBM-Q and D-Wave using $\tau=5$

Fig 2: Most Central Nodes Using QUBO and NetworkX’s EC. The graph $G_{8}$
showing the most central nodes using the (a) NetworkX EC algorithm, and QUBO
with (b) $\tau=1$ and (c) $\tau=5$ on IBM-Q and D-Wave quantum computers for
$P_{0}=\frac{1}{\sqrt{n}}$ and $P_{1}=5n$ where $n$ is the number of nodes of
the graphs. For (a), the most central node(s) is(are) of brighter colors and
are encircled by larger circles. In (b) and (c), the bright colored (yellow)
nodes are the most central nodes whiles the dark colored (purple) nodes are
the least central nodes relative to $\tau$.

With these results, we moved forward to the construction of a quantum
minimization problem:

$\min_{\boldsymbol{x}\in\\{0,1\\}^{n}}\boldsymbol{x}^{T}Q\boldsymbol{x}$

with the objective of identifying the top $\tau$ most important nodes within
our graphs, as confirmed by the results from NetworkX and the classical
minimization mentioned previously. Note that here, the goal is to capture the
top $\tau$ most important nodes of the graph irrespective of their order.

Experimenting with the QUBO Eq. (11) on both the D-Wave 2000Q and IBM-Q
devices (QASM simulator and ibmq_manhattan for the Karate club graph), we
obtained results showing that the formulation indeed captures EC and not
degree centrality ranking when compared with the NetworkX output. Fig. 2b
shows the output for a search for the most central node (colored yellow) which
in this case requires $\tau=1$ of the graph $G_{8}$ in S1 Fig.. Compared to
the NetworkX output in Fig. 2a, the quantum computing scheme supports
proposition 1, the claim that the ground state of the QUBO encapsulates the
information about the nodes with high ranked EC values and not those of high
ranked degree centrality values. Fig. 2c, shows the result obtained for a
search for the top $\tau=5$ most important nodes (in yellow) of the graph
$G_{8}$. For the NetworkX results, the most central nodes are identified by
the size and brightness of the color of the disk around the nodes; the larger
and brighter the disk the more central the node. For the graph $G_{8}$, the
yellow node ($0$) is the most central, followed by nodes $1,2,3$ and the nodes
$4,5,6,\ldots,14$ are all of the same centrality value and are the least
central nodes. It was observed that it is sufficient to choose the penalty
constants $P_{0}=\frac{1}{\sqrt{n}}$, and $P_{1}>P_{0}$ (in our case
$P_{1}=5n$ worked well).

To obtain optimal results using the D-Wave 2000Q, it is best to set the chain
strength to the maximum possible, $1000$ in this case. It was observed that
very low chain strength values resulted in broken chains affecting the
probability of the QA settling into a global minimum solution. Setting the
post-processing method to “optimization” and the number of samples to the
maximum, $10,000$ boosted the chances of obtaining the optimal solution.
However as the graph gets larger more samples are required to increase the
probability of observing the ground state solution. Here is where some
inconsistencies in obtaining the lowest energy output for larger graphs was
highlighted: the first (or second, etc.) run may not result in the expected
output. Running multiple times would eventually result in the correct output
occurring once, but due to the noisy and quantum nature of the quantum
machine, we do not know the exact number of runs required to give the global
minimum solution. This was mostly observed in the Karate Club graph in Fig. 3e
and the Davis Southern women network. With a little bit of luck the result can
be obtained in the first run or in a few runs. This behavior is not surprising
since an increase in the problem size decreases the probability of finding an
optimal solution due to annealing error and imperfect hardware [51].

The quadratic solver QBSolv in Ocean served as a benchmark for determining the
correct minimum of the D-Wave annealing output. QBSolv provided expected
results of the minimization on most occasions with degenerate ground state
solutions in some instances. Degeneracy here refers to the outputs with the
same minimum energy value of the QUBO due to multiple nodes having the same
centrality value. Whenever there is degeneracy in the QBSolv output, one of
the solutions correctly identifies the top $\tau$ nodes. For example, in Fig.
2c, we have 12 degenerate solutions corresponding to the 12 leaves of the tree
graph. The solution graphed included node $14$ in the top $5$ important nodes,
however node $4,5,6,7,8,9,10$, $11,12,13$ or $15$ are valid replacements for
node $14$ in this output since column 3 of Table 2 shows that all these nodes
have the same centrality values.

Solving the problem for the QUBO in Eq. (11) for the graph in Fig. 1 and graph
$G_{8}$ in S2 Fig. correctly identify the most important nodes for any given
$\tau$, including selecting node $0$ as the node with the highest EC (see Fig.
2). The penalty weights used here were $P_{0}=\frac{1}{\sqrt{n}}$ and
$P_{1}=5n$ for $n$ number of nodes of the graphs. The same results were
obtained from experiments on IBM-Q’s QASM simulator using QAOA. However the
ground state solution for the graphs with smaller numbers of nodes ($n\leq
16$) were obtained on a single run with the penalty constants
$P_{0}=\frac{1}{\sqrt{n}}$ and $P_{1}=5n$. When considering graphs with larger
numbers of nodes, the program had to be run multiple times using the same
penalty constants above to be able to capture the global minimum solution.

Solving for the ground state solution of the QUBO in Eq. (11) for the graph
$G_{1}$ and Karate Club graph in S2 Fig. was quite challenging. On most
occasions, the solution for the Karate Club graph with $\tau\geq 3$ required
multiple runs before settling on the global minimum solution. In other words,
this graph required more samples to be able to output the global minimum. For
$3\leq\tau\leq 5$, the QUBO couldn’t capture the ordering that matched that of
NetworkX rankings when using the penalty constants $P_{0}=\frac{1}{\sqrt{n}}$
and $P_{1}=5n$ for the graph $G_{1}$. The nodes $0,3,4$ were always skipped in
the search for the top $\tau=3,4,5$. From Fig. 3d, we see that the output for
the top $6$ most important nodes excludes node $4$ and includes node $2$ which
shouldn’t be the case since from Fig. 3a, node $4$ is more central than node
$2$. The exact reason for this occurrence for this particular graph is
unknown, however it seems the program picks only one of the degenerate nodes
$3$, or $4$ and moves on to select the next central node $2$.

(a) NetworkX result for the graph $G_{1}$

(b) NetworkX result for the Karate Club graph

(c) NetworkX result for the Barabasi-Albert, $BA(50,5)=G_{4}$

(d) D-Wave_2000Q result for the graph $G_{1}$, $\tau=6$

(e) D-Wave_2000Q result for the Karate Club graph, $\tau=5$

(f) D-Wave_2000Q result for the Barabasi-Albert graph $BA(50,5)=G_{4}$ with
$\tau=3$

(g) D-Wave_2000Q result for the graph $G_{1}$, $\tau=1$

(h) D-Wave_2000Q result for the Karate Club graph, $\tau=1$

(i) D-Wave_2000Q result for the Barabasi-Albert graph $BA(50,5)=G_{4}$ with
$\tau=1$

Fig 3: Results for challenging graphs. NetworkX results for the challenging
graphs encountered; (a) graph $G_{1}$, (b) Karate Club graph, and (c)
Barabasi-Albert graph, $BA(50,5)=G_{4}$. D-Wave results obtained for the QUBO
in Eq. 11 using $P_{0}=\frac{1}{\sqrt{n}}$ and $P_{1}=5n$ for (d) graph
$G_{1}$ with $\tau=6$, (e) Karate Club graph with $\tau=5$, (f)
$BA(50,5)=G_{4}$ with $\tau=3$, (g) $G_{1}$ with $\tau=1$, (h) Karate Club
graph with $\tau=1$ and (i) $BA(50,5)=G_{4}$ with $\tau=1$. The D-Wave results
colors the top $\tau$ most important nodes yellow and the least central nodes
purple. The NetworkX result identifies most central nodes with larger and
brighter circles (yellow being the top) and least central nodes with smaller
and darker circles (least being purple).

With these interesting results, we further examined the possibility of
defining a hierarchy of nodes in the network from the QUBO results obtained
from the D-Wave 2000Q and IBM Q. By hierarchy we mean a ranking that orders
the nodes based on importance or influence using EC. The hierarchy can then be
used to identify for example super-spreaders of disease. We compare the
hierarchy of nodes obtained using our QUBO formulation with that obtained from
NetworkX when using the EC algorithm. The result obtained for graph $M$ is
shown in Table 3. To determine the node rank, we consider the set of $\tau$
nodes obtained from the QUBO results for each value of $0<\tau\leq n$ and
compute the symmetric difference. The $i$th rank is determined by taking the
difference of the QUBO result for $\tau=i$ and $\tau=i-1$. For example, to
rank the nodes for graph $M$, the 1st most important node is obtained by
running the QUBO for $\tau=1$. The 2nd most important node is determined by
taking the difference between the QUBO results for $\tau=2$, $\\{0,2\\}$ and
the result $\\{2\\}$ of $\tau=1$. Thus $\\{0,2\\}\setminus\\{2\\}=\\{0\\}$
implies that node $0$ is the second most important node. Continuing in this
manner the importance of the nodes can be ranked. Computing the difference
$\\{0,1,2,3,4,5,6,7,9,11\\}\setminus\\{0,1,2,3,4,5,6,9,11\\}=\\{7\\}$ for
$\tau=10$ and $\tau=9$ we see that the $10$th important node is node $7$.

Table 3: Node rank of graph $M$ using QUBO results.

$\tau$ |  (QA/QAOA) QUBO results for top $\tau$ most central nodes of graph $M$ |  Node Rank of (QA/QAOA) QUBO results for $M$ |  Node Rank of NetworkX EC result for $M$
---|---|---|---
1 |  2 |  2 |  2
2 |  0, 2 |  0 |  0
3 |  0, 2, 9 |  9 |  9
4 |  0, 1, 2, 9 |  1 |  1
5 |  0, 1, 2, 6, 9 |  6 |  6
6 |  0, 1, 2, 3, 6, 9 |  3 |  3
7 |  0, 1, 2, 3, 5, 6, 9 |  5 |  5
8 |  0, 1, 2, 3, 5, 6, 9, 11 |  11 |  11
9 |  0, 1, 2, 3, 4, 5, 6, 9, 11 |  4 |  7
10 |  0, 1, 2, 3, 4, 5, 6, 7, 9, 11 |  7 |  8
11 |  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11 |  8 |  4
12 |  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12 |  12 |  10
13 |  0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 |  10 |  12

Determining node rank from the QUBO result obtained from D-Wave (QA) and IBM-Q
(QAOA) using symmetric difference of results for $\tau=1$ to $\tau=n$. In
column 4, the ith in rank is determined by taking the difference of the result
for $\tau=i$ and $\tau=i-1$ e.g. the 1st (most central) node $2$ is determined
by implementing the QUBO for $\tau=1$, the 2nd is determined by taking the
difference between the result of $\tau=2$, $\\{0,2\\}$ and the result of
$\tau=1$,$\\{2\\}$, i.e. $\\{0,2\\}\setminus\\{2\\}=\\{0\\}$, implying node
$0$ is ranked second in importance.

(a) NetworkX ranking

(b) Ranking obtained from IBM-Q and D-Wave for $\tau=1$

(c) Ranking obtained from IBM-Q and D-Wave using $\tau=10$

Fig 4: Results for graph $M$. The graph $M$ showing the most central nodes
using (a) the NetworkX EC algorithm, and ((b) and (C)) QUBO on the IBM-Q and
D-Wave quantum computers for (b) $\tau=1$ and (c) $\tau=10$ with
$P_{0}=\frac{1}{\sqrt{n}}$ and $P_{1}=5n$ where $n$ is the number of nodes of
the graphs. For (a), the most central node(s) is(are) of brighter colors and
are encircled by larger circles. In (b) and (c), the bright colored (yellow)
nodes are the most central nodes whiles the dark colored (purple) nodes are
the least central nodes relative to $\tau$.

## Conclusion

We have formulated and shown a QUBO problem for the identification of most
important nodes or hot-spots in a graph based on the EC measure is possible.
Using quantum computing algorithms such as quantum annealing on the D-Wave
2000Q and QAOA on IBM-Q, our formulation was able to correctly identify all
$\tau<n$ most important nodes for graphs with less than $17$ nodes ($n<17$).
For graphs with more than 16 nodes, the quantum computing algorithm always
identified the top $\tau\leq 6$ most important nodes for all the graphs
considered correctly except for the Davis Southern women network and the tree
graph $G_{1}$ whose outputs for some values of $1<\tau\leq 6$ showed some
marginal inconsistencies. Marginal because the differences are negligible. For
example, for the graph $G_{1}$, the node $3$ was not selected since it’s
centrality value is the same as that of node $4$. Despite this challenge, the
results obtained from all graphs considered for $\tau=1$ verified Proposition
1. We have also demonstrated the feasibility of defining a hierarchy of nodes
in a graph from our formulation using the QUBO results from the D-Wave 2000Q
and IBM-Q’s QASM simulator. Given that the current quantum resources (D-Wave
2000Q and IBM-Q’s Manhattan) at our disposal limit us on the size of graphs to
explore and in the presence of uncontrollable noise which affects the
probabilities of obtaining quality results, we were unable to experiment with
real life data. Therefore in the future when powerful and less noisy quantum
computers are made available for our perusal, we wish to test our hypothesis
further to establish a more generalized formulation that works for all $\tau$
on all graphs. That is, we want to verify

###### Claim 1

For any graph $G=(V,E)$, with adjacency matrix $A$ and degree sequence $D$.
The indices of the nonzero elements of the ground state solution to the QUBO
problem

$\displaystyle\min_{\boldsymbol{x}\in\\{0,1\\}^{n}}\boldsymbol{x}^{T}Q\boldsymbol{x}$
$\displaystyle Q$
$\displaystyle=-P_{0}A^{2}\boldsymbol{\hat{d}}\boldsymbol{\hat{d}}^{T}A-P_{0}A\boldsymbol{\hat{d}}\boldsymbol{\hat{d}}^{T}A^{2}+P_{1}C$
$\displaystyle C$ $\displaystyle=(1-2\tau)I+U$

where the matrix $U$ is defined in Eq. (9), corresponds to the $\tau$ most
central nodes via EC measure of the graph $G$ for any $\tau\leq n=|V|$.

For this claim we wish to investigate both directed and undirected graphs,
“will replacing $A^{2}$ by $AA^{T}$ or $A^{T}A$ in $Q$ still work for directed
graphs or will it require a modified QUBO?”

## Supporting information

#### S1 Appendix.

Computing Degree Centrality and Eigencentrality from Exponential Function.

Consider the exponential function defined by

$\displaystyle f(x)=e^{x}=\sum_{k=0}^{\infty}\frac{x^{k}}{k!}$ (13)

Since $f$ is continuous and analytic with radius of convergence infinity, we
observe that for a primitive matrix $A$ which is mostly the case for the
adjacency matrix of simply connected undirected graphs

$\displaystyle f(\gamma A)=e^{\gamma A}=\sum_{k=0}^{\infty}\frac{(\gamma
A)^{k}}{k!}$ (14)

Then $f(\gamma A)$ is a matrix and the components of the vector $f(\gamma
A)\boldsymbol{1}$ where $\boldsymbol{1}$ is a vector of ones, counts the walks
of infinite length centered at each node. Let
$\\{\boldsymbol{e_{i}}\,\in\,\mathbb{R}^{n}\\}$ be the set of canonical basis
vectors for $\mathbb{R}^{n}$ whose only nonzero element is the ith component.

$\displaystyle v_{i}(\gamma)$
$\displaystyle=\sum_{j=0}^{n-1}\boldsymbol{e_{i}}^{T}f(\gamma
A)\boldsymbol{e_{j}}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{j=0}^{n-1}\boldsymbol{e_{i}}^{T}\frac{(\gamma
A)^{k}}{k!}\boldsymbol{e_{j}}$
$\displaystyle=\sum_{k=0}^{\infty}\sum_{j=0}^{n-1}a_{k}\gamma^{k}\boldsymbol{e_{i}}^{T}A^{k}\boldsymbol{e_{j}}$
$\displaystyle=a_{0}\sum_{j=0}^{n-1}\boldsymbol{e_{i}}^{T}\boldsymbol{e_{i}}+a_{1}\gamma\sum_{j=0}^{n-1}\boldsymbol{e_{i}}^{T}A\boldsymbol{e_{j}}+\ldots$
$\displaystyle=a_{0}+a_{1}\gamma
d_{i}+\sum_{k=0}^{\infty}\sum_{j=0}^{n-1}a_{k}\gamma^{k}\boldsymbol{e_{i}}^{T}A^{k}\boldsymbol{e_{j}}$
$\displaystyle c_{i}(\gamma)=\frac{v_{i}(\gamma)-a_{0}}{a_{1}\gamma}$
$\displaystyle=d_{i}+\sum_{k=2}^{\infty}\sum_{j=0}^{n-1}\frac{a_{k}}{a_{1}}\gamma^{k-1}\boldsymbol{e_{i}}^{T}A^{k}\boldsymbol{e_{j}}$
(15)

where $d_{i}$ is the degree of the $i$th node. In the limit
$\gamma\longrightarrow 0^{+}$, Eq. (15) converges to the degree centrality and
by expanding in the eigenbasis, it converges to EC for
$\gamma\longrightarrow\infty$ [10].

#### S1 Fig.

Small graphs considered. (a) Sedgewick-Maze graph, (b) Graph $G_{8}$, (c)
Graph $G_{9}$, (d) Bull graph $G_{10}$, (e) $K_{3,5}$, (f) $K_{5}$, (g)
$K_{6}+S_{6}$, (h) Florentine families graph $F$, (i) Graph $M$.

(a) Sedgewick-Maze graph

(b) Graph $G_{8}$

(c) Graph $G_{9}$

(d) Bull graph $G_{10}$

(e) $K_{3,5}$

(f) $K_{5}$

(g) Graph $K_{6}+S_{6}$

(h) Florentine families graph $F$

(i)  Graph $M$

#### S2 Fig.

Large graphs considered. (a) Graph $G_{1}$, (b) Davis southern women network,
$G_{2}$ (c) Tutte graph $G_{3}$, (d) Barabasi-Albert graph $G_{4}$, (e)
Lollipop $L_{10,20}$ graph $G_{7}$, (f) Karate Club graph.

(j) Graph $G_{1}$

(k) Davis southern women network, $G_{2}$

(l) Tutte graph $G_{3}$

(m) Barabasi-Albert graph $G_{4}$

(n) Lollipop $L_{10,20}$ graph $G_{7}$

(o) Karate Club graph

## Acknowledgments

Research presented in this article was supported by the Laboratory Directed
Research and Development (LDRD) program of Los Alamos National Laboratory
(LANL) under project number 20200056DR. LANL is operated by Triad National
Security, LLC, for the National Nuclear Security Administration (NNSA) of the
U. S. Department of Energy (contract no. 89233218CNA000001). TJ acknowledges
support from the U.S. Department of Energy (DOE) through a quantum computing
program sponsored by the Los Alamos National Laboratory (LANL) Information
Science & Technology Institute. This research was also supported by the U. S.
Department of Energy (DOE) National Nuclear Security Administration (NNSA)
Advanced Simulation and Computing (ASC) program at LANL. We acknowledge the
ASC program at LANL for use of their Ising D-Wave 2000Q quantum computing
resource. Quantum resources from the IBM-Q Hub are also acknowledged.
Assigned: Los Alamos Unclassified Report LA-UR-21-24030.

## References

  * 1. Freeman LC.  Centrality in social networks conceptual clarification.  Social Networks. 1978-1979;1(3):215–239.
  * 2. Sabidussi G.  The centrality index of a graph.  Psychometrika. 1966;31:581–603.
  * 3. Chea E, Livesay DR.  How accurate and statistically robust are catalytic site predictions based on closeness centrality?  Bmc Bioinformatics. 2007;8(1):1–14.
  * 4. Freeman LC.  A Set of Measures of Centrality Based on Betweenness.  Sociometry. 1977;40(1):35–41.
  * 5. Yu H, Kim PM, Sprecher E, Trifonov V, Gerstein M.  The importance of bottlenecks in protein networks: correlation with gene essentiality and expression dynamics.  PLoS Comput Biol. 2007;3(4):e59.
  * 6. Katz L.  A New Status Index Derived from Sociometric Index.  Psychometrika. 1953; p. 39–43.
  * 7. Zhao J, Yang TH, Huang Y, Holme P.  Ranking candidate disease genes from gene expression and protein interaction: a Katz-centrality based approach.  PloS one. 2011;6(9):e24306.
  * 8. Brin S, Page L.  The anatomy of a large-scale hypertextual Web search engine.  Computer Networks and ISDN Systems. 1998;30(1-7):107–117.
  * 9. Benzi M, Klymko C.  Total communicability as a centrality measure.  Journal of Complex Networks. 2013;1(2):124–149.
  * 10. Benzi M, Klymko C.  On the Limiting Behavior of Parameter-Dependent Network Centrality Measures.  SIAM Journal on Matrix Analysis and Applications. 2015;36(2):686–706.  doi:10.1137/130950550.
  * 11. Schoch D, Valente TW, Brandes U.  Correlations among centrality indices and a class of uniquely ranked graphs.  Social Networks. 2017;50:46–54.  doi:https://doi.org/10.1016/j.socnet.2017.03.010.
  * 12. Ronqui JRF, Travieso G.  Analyzing complex networks through correlations in centrality measurements.  Journal of Statistical Mechanics: Theory and Experiment. 2015;2015(5):P05030.  doi:10.1088/1742-5468/2015/05/p05030.
  * 13. Mihail M, Papadimitriou C.  On the Eigenvalue Power Law.  In: Randomization and Approximation Techniques in Computer Science. RANDOM 2002. Springer,Berlin, Heidelberg; 2002. p. 254–262.  Available from: https://doi.org/10.1007/3-540-45726-7_20.
  * 14. Bonacich P.  Power and Centrality: A Family of Measures.  American Journal of Sociology. 1987;92(5):1170–1182.
  * 15. Negre C, Morzan U, Hendrickson H, Pal R, Lisi G, Loria J, et al.  Eigenvector centrality for characterization of protein allosteric pathways.  Proc Natl Acad Sci USA. 2018;115(52):E12201–E12208.
  * 16. Jimenez-Martinez J, Negre CFA.  Eigenvector centrality for geometric and topological characterization of porous media.  Phys Rev E. 2017;96:013310.  doi:10.1103/PhysRevE.96.013310.
  * 17. Azad S, Devi S.  Tracking the spread of COVID-19 in India via social networks in the early phase of the pandemic.  Journal of Travel Medicine. 2020;27(8).
  * 18. Navaretti GB, Calzolari G, Dossena A, Lanza A, Pozzolo AF.  In and out lockdowns: Identifying the centrality of economic activities.  Covid Economics. 2020;17:189–204.
  * 19. Newman M.  Networks.  OUP Oxford; 2018.  Available from: https://books.google.com/books?id=YdZjDwAAQBAJ.
  * 20. Perron O.  Zur Theorie der Matrices.  Math Ann. 1907;64:248–263.
  * 21. Quantum computing in a nutshell; 2021.  Available from: https://qiskit.org/documentation/qc_intro.html#quantum-computing-in-a-nutshell.
  * 22. Welcome to D-Wave; 2020.  https://docs.dwavesys.com/docs/latest/c_gs_1.html.
  * 23. D-Wave QPU Architecture: Chimera; 2020.  Available from: https://docs.dwavesys.com/docs/latest/c_gs_4.html [cited 2020 Dec 25].
  * 24. Minor-Embedding a Problem onto the QPU; 2020.  Available from: https://docs.dwavesys.com/docs/latest/c_gs_7.html#em [cited 2020 Dec 21].
  * 25. Terry JP, Akrobotu PD, Negre CFA, Mniszewski SM.  Quantum isomer search.  PLoS ONE. 2020;15(1).  doi:https://doi.org/10.1371/journal.pone.0226787.
  * 26. Ushijima-Mwesigwa H, Negre CFA, Mniszewski SM.  Graph Partitioning Using Quantum Annealing on the D-Wave System.  In: Proceedings of the Second International Workshop on Post Moores Era Supercomputing. PMES’17. New York, NY, USA: ACM; 2017. p. 22–29.  Available from: http://doi.acm.org/10.1145/3149526.3149531.
  * 27. Negre CFA, Ushijima-Mwesigwa H, Mniszewski SM.  Detecting multiple communities using quantum annealing on the D-Wave system.  PLoS ONE. 2020;15(2).  doi:https://doi.org/10.1371/journal.pone.0227538.
  * 28. Bauckhage C, Brito E, Cvejoski K, Ojeda C, Sifa R, Wrobel S.  Ising Models for Binary Clustering via Adiabatic Quantum Computing.  In: Pelillo M, Hancock E, editors. Energy Minimization Methods in Computer Vision and Pattern Recognition. Cham: Springer International Publishing; 2018. p. 3–17.
  * 29. Calude CS, Dinneen MJ, Hua R.  QUBO formulations for the graph isomorphism problem and related problems.  Theoretical Computer Science. 2017;701:54–69.  doi:https://doi.org/10.1016/j.tcs.2017.04.016.
  * 30. Pudenz KL, Lidar DA.  Quantum adiabatic machine learning.  Quantum Inf Process. 2013;12:2027–2070.  doi:https://doi.org/10.1007/s11128-012-0506-4.
  * 31. Adachi SH, Henderson MP.  Application of Quantum Annealing to Training of Deep Neural Networks.  arXiv e-prints. 2015; p. arXiv:1510.06356.
  * 32. Carnevali V, Siloi I, Di Felice R, Fornari M.  Vacancies in graphene: an application of adiabatic quantum optimization.  Phys Chem Chem Phys. 2020;22:27332–27337.  doi:10.1039/D0CP04037A.
  * 33. Venturelli D, Marchand DJJ, Rojo G. Quantum Annealing Implementation of Job-Shop Scheduling; 2016.
  * 34. Stollenwerk T, Lobe E, Jung M.  Flight Gate Assignment with a Quantum Annealer.  In: Feld S, Linnhoff-Popien C, editors. Quantum Technology and Optimization Problems. Cham: Springer International Publishing; 2019. p. 99–110.
  * 35. Stollenwerk T, O’Gorman B, Venturelli D, Mandrà S, Rodionova O, Ng H, et al.  Quantum Annealing Applied to De-Conflicting Optimal Trajectories for Air Traffic Management.  IEEE Transactions on Intelligent Transportation Systems. 2020;21(1):285–297.  doi:10.1109/TITS.2019.2891235.
  * 36. Michielsen K, Nocon M, Willsch D, Jin F, Lippert T, De Raedt H.  Benchmarking gate-based quantum computers.  Computer Physics Communications. 2017;220:44–55.  doi:https://doi.org/10.1016/j.cpc.2017.06.011.
  * 37. IBM. IBM Quantum Experience is quantum on the cloud; 2021.  Available from: https://www.ibm.com/quantum-computing/experience/.
  * 38. Childs AM, Goldstone J.  Spatial search by quantum walk.  Phys Rev A. 2004;70:022314.  doi:10.1103/PhysRevA.70.022314.
  * 39. Numpy; 2020.  Available from: https://www.numpy.org/doc/stable [cited 2020 Dec 25].
  * 40. SciPy; 2020.  Available from: https://docs.scipy.org/doc/scipy/reference/ [cited 2020 Dec 25].
  * 41. D-Wave’s Ocean Software; 2020.  Available from: https://ocean.dwavesys.com/ [cited 2020 Dec 25].
  * 42. Qiskit; 2020.  Available from: https://qiskit.org/ [cited 2020 Dec 25].
  * 43. Hagberg AA, Schult DA, Swart PJ.  Exploring Network Structure, Dynamics, and Function using NetworkX.  In: Proceedings of the 7th Python in Science Conference (SciPy 2008). SciPy 2008. ACM; 2008. p. 11–16.
  * 44. Hunter JD.  Matplotlib: A 2D graphics environment.  Computing in Science & Engineering. 2007;9(3):90–95.  doi:10.1109/MCSE.2007.55.
  * 45. D-Wave. Los Alamos National Laboratory Upgrades to D-Wave 2000Q™Quantum Computer; 2019.  https://www.dwavesys.com/press-releases/los-alamos-national-laboratory-upgrades-d-wave-2000q%E2%84%A2-quantum-computer.
  * 46. Converters for Quadratic Programs; 2020.  Available from: https://qiskit.org/documentation/tutorials/optimization/2_converters_for_quadratic_programs.html.
  * 47. Minimum Eigen Optimizer; 2020.  Available from: https://qiskit.org/documentation/tutorials/optimization/3_minimum_eigen_optimizer.html.
  * 48. Farhi E, Goldstone J.  A Quantum Approximate Optimization Algorithm.  arXiv preprint arXiv:14114028. 2014;.
  * 49. Developers N. Graph generators; 2020.  Available from: https://networkx.org/documentation/stable/reference/generators.html [cited 2021 March 9].
  * 50. Wikipedia contributors. Barabási–Albert model — Wikipedia, The Free Encyclopedia; 2021.  https://en.wikipedia.org/w/index.php?title=Barab%C3%A1si%E2%80%93Albert_model&oldid=998323343.
  * 51. Pudenz KL, Albash T, Lidar DA.  Quantum annealing correction for random Ising problems.  Physical Review A. 2015;91(4):042302.