
% this file is called up by thesis.tex
% content in this file will be fed into the main document

%: ----------------------- introduction file header -----------------------
\chapter{Introduction}\label{C:introduction}

\section{Background and motivation}

\subsection{Introduction to complex networks}
Complex systems are ubiquitous throughout the world, both in nature and within man-made structures.  Over the past decade, large amounts of network data have become available and, correspondingly, the analysis of complex networks has become increasingly important.  In response, researchers in diverse areas have focused their attention on such problems and the study of complex systems has begun to develop into a discipline in its own right \cite{Bocca06,Brandes,Caldarelli,Ebook11,Ebook10,New10,NewBarWat03}. Mathematicians, physicists, computer scientists, biologists, and social scientists (among others) have approached these problems from various angles\cite{stareig,crofootetal,Proteomics,Proteome,Katz,dhillon}.  This large amount of research within diverse scientific communities has produced many interesting results, but there are still many fundamental questions that are not fully answered.  Additionally, as available datasets become larger and larger, efficient data structures and computational methods become paramount. 

One of the fundamental questions in network analysis is to determine the ``most important" elements in a given network.  The interpretation of what is meant by ``important" can change from application to application.  Measures of node importance are usually referred to as {\em node centralities}, and many centrality measures have been proposed, starting with the simplest of all, the node degree. This crude metric has the drawback of being too ``local", as it does not take into effect the connectivity of the immediate neighbors of the node under consideration. A number of more sophisticated centrality measures have been introduced that take into account the global connectivity properties of the network. These include various types of eigenvector centrality for both directed and undirected networks, betweenness centrality, and many others.  Overviews of various centrality measures can be found in \cite{Bocca06,Bonacich,Brandes,Ebook11,IRSurvey,New10,NewBarWat03}.  The centrality scores can be used to provide {\em rankings} of the nodes in the network. There are many different ranking methods in use (most of which depend on centrality measures), and many algorithms have been developed to compute these rankings.  Information about the many different ranking schemes can be found in \cite{BEK,Ebook11,Katz,HITS,IRSurvey,Pagerank,Whos,SALSA}.


\subsection{Outline of thesis}

This thesis investigates the use of matrix functions in the analysis of node centrality.  In the remainder of this chapter, we will provide background information from various areas of mathematics.  Sections \ref{sec:graph_theory} and \ref{sec:linear_algebra} present basic concepts from graph theory and linear algebra that will be used in this thesis.  Section \ref{sec:complex_networks} discusses properties found in many complex networks. Section \ref{sec:centrality} provides an overview of commonly used centrality and communicability measures, including HITS and PageRank while Section \ref{sec:comparing} discusses the comparison of rankings produced by different centrality measures.  Section \ref{sec:iterative_methods} provides details on iterative methods that are used in this thesis to approximate certain types of centrality scores.

Chapter \ref{C:TotalComm} examines the use of a node centrality measure based on the notion of {\em total communicability}, defined in terms of the row sums of the exponential of the adjacency matrix of the network.  We argue that this is a natural metric for ranking nodes in a network, and we point out that it can be computed very rapidly even in the case of large networks. Furthermore, we propose the total sum of node communicabilities as a useful measure of network connectivity.  Extensive numerical studies are conducted in order to compare this centrality measure with the closely related ones of subgraph centrality and Katz centrality on both synthetic and real-world networks.

In Chapter \ref{C:Robustness}, we demonstrate an analytical relationship between rankings produced by parameterized centrality measures based on the matrix exponential and resolvent with those produced by degree and eigenvector centrality.  The centrality measures considered are exponential and resolvent subgraph centrality (defined in terms of the diagonal entries of the matrix exponential and matrix resolvent, respectively), total communicability, and Katz centrality (defined in terms of the row sums of the matrix exponential and resolvent).  This analysis helps to explain the observed robustness of these rankings on many real world networks, even though the scores produced by the centrality measures are not stable.  A number of numerical experiments on both directed and undirected networks are conducted to analyze how this relationship affects node rankings.

In Chapter \ref{C:HubsAuths}, we propose an extension of the measures of subgraph centrality and communicability to directed networks, and we apply them to the problem of ranking hubs and authorities.  The extension is achieved by {\em bipartization}, i.e., the directed network is mapped onto a bipartite undirected network with twice as many nodes in order to obtain a network with a symmetric adjacency matrix.  
We explicitly determine the exponential of this adjacency matrix in terms of the adjacency matrix of the original, directed network, and we give an interpretation of centrality and communicability in this new context, leading to a technique for ranking hubs and authorities.  The matrix exponential method for computing hubs and authorities is compared to the well known HITS algorithm, both on small artificial examples 
and on more realistic real-world networks. A few other ranking algorithms are also discussed and compared with our technique. 

Finally, Chapter \ref{C:conclusions} contains concluding remarks and guidelines for future work.

The work presented in Chapters \ref{C:TotalComm} and \ref{C:Robustness} is based on \cite{BKlym13} and \cite{BEK}, respectively.

\section{Basic concepts from graph theory}
\label{sec:graph_theory}

In this section, we provide an overview of the basic concepts, definitions, and terminology from graph theory that will be used throughout the thesis.  A more comprehensive presentation can be found in \cite{Bollobas,Die00}.

A \emph{graph} $G=(V,E)$ is formed by a set of nodes (vertices) $V$ and edges (links) $E$ formed by unordered pairs of vertices.  Every network is naturally associated with a graph $G=(V,E)$ where $|V| = n$ is the number of objects in the network and $E$ ($|E|=m$) is the collection of connections between objects,  $E=\{(i,j) \, |$ there is an edge between node $i$  and node $j \}$.  Nodes represent the elements that make up the network and edges represent various types of interactions among network elements.  The {\em degree} $d_i$ of a vertex $i$ is the number of edges incident to $i$.  A graph is {\em $k$-regular} if every node has degree $d_i = k$.

A directed graph, or \emph{digraph}, $G=(V,E)$ is formed by a set of vertices $V$ and edges $E$ formed by ordered pairs of vertices.  That is, $(i,j) \in E \not\Rightarrow (j,i) \in E$.  In the case of digraphs, which model directed networks, there are two types of degree.  The {\em in-degree} of node $i$, $d_i^{in}$, is given by the number of edges which point to $i$.  The {\em out-degree}, $d_i^{out}$, is given by the number of edges pointing out from $i$. 

The first result in graph theory, which was also the start of graph theory as a proper mathematical discipline, related the number of nodes, $n$, in an undirected graph to the number of edges, $m$.  The following result, called the {\em handshaking lemma}, is from Leonhard Euler \cite{Euler} and has since been extended to directed graphs.

\begin{lemma}
\label{lemma:handshaking}
 Let $G$ be a graph with $n$ nodes and $m$ edges.  Then, 
 \begin{enumerate}[(i)]
	\item if $G$ is undirected, then $$\sum_{i=1}^n d_i = 2m,$$
	\item if $G$ is directed, then $$\sum_{i=1}^n d_i^{in} = \sum_{i=1}^n d_i^{out} = m.$$
\end{enumerate}
\end{lemma}

A {\em subgraph} of $G$ is a graph $H = (V', E')$ such that $V' \subseteq V$ and $E' \subseteq \{ (i,j) \in E \,|\, i,j \in V'\}$.  The subgraph $H$ is said to be an {\em induced subgraph} if $E' = \{ (i,j) \in E \,|\, i,j \in V'\}$, that is if all edges possible are present.  A subgraph $H$ is said to be {\em maximal} in regards to a property if adding any additional nodes or edges to $H$ would result in that property no longer holding.

A {\em walk} of length $k$ in $G$ is a set of nodes $i_1, i_2, \ldots i_k, i_{k+1}$ such that for all $1 \leq l \leq k$, there is an edge between $i_l$ and $i_{l+1}$ (or a directed edge from $i_l$ to $i_{l+1}$ in the case of a digraph).  A {\em closed walk} is a walk where $i_1 = i_{k+1}$.  A {\em path} is a walk with no repeated nodes.  A {\em cycle} is a closed path.  A graph $G$ is {\em walk-regular} if the number of closed walks of length $k \geq 0$ starting at a given node is the same for every node in the graph.  A graph is {\em simple} if it has no {\em loops} (edges from a node $i$ to itself), no multiple edges, and unweighted edges.

The {\em distance}, $d(i,j)$, between nodes $i$ and $j$ in a graph is given by the length of the shortest path between them.  If there is no path between two nodes, they are considered to be at an infinite distance.  The {\em diameter} of a graph, $d_G$ is given by the maximum shortest path length in the graph: $d_G = \max_{i \neq j} \, d(i,j)$.  In an undirected graph, the distance function forms a {\em metric}.  That is, the following lemma holds.
\begin{lemma}
\label{lemma:distance_metric}
Given an undirected graph $G$, the {\em distance function} $d : V \times V \to \mathbb{R}$ given by $d(i,j) = \{$the length of the shortest path between $i$ and $j \}$ satisfies the following properties and is therefore a metric: for all nodes $i$, $j$, and $k$,
\begin{enumerate}[(i)]
	\item $d(i,j) \geq 0$  and $d(i,j) = 0 \Leftrightarrow i = j$,
	\item $d(i,j) = d(j,i)$,
	\item $d(i,j) \leq d(i,k) + d(k,j)$.
\end{enumerate}
\end{lemma}
In a directed graph, there is no guarantee that distances are symmetric and, thus, the distance function does not define a metric.

An undirected graph is {\em connected} if there exists a path between every pair of nodes.  A digraph is {\em weakly connected} if the underlying undirected graph is connected.  It is {\em strongly connected} if, for every pair of vertices $i$ and $j$, there is both a directed walk from $i$ to $j$ and one from $j$ to $i$.  If an undirected graph $G$ is not connected, then $G$ can be partitioned into its {\em connected components}.  The number of connected components is the number of maximal connected subgraphs.  If $G$ is directed, then we can distinguish between the {\em strongly connected components} and the {\em weakly connected components}.  When a graph is partitioned into its (weakly) connected components, every node and edge is in exactly one component.  When a digraph is partitioned into its strongly connected components, every node is in exactly one component and every edge is in at least one component.  Often, in the analysis of complex networks we are concerned with the largest connected component of an undirected network or the largest strongly connected component of a directed network.

Often, we consider graphs with unweighted edges, but depending on the type of interaction edges can be weighted or unweighted.  Weights on edges can be used to measure the amount of traffic between two nodes or other properties of data, such as the vulnerability of the connection to attack or disruption or a physical distance between two nodes.  In the case of unweighted networks, we consider the weight of each edge to be equal to one.

Every graph is naturally associated with an {\em adjacency matrix}.  The adjacency matrix of an unweighted graph $G$ is given by the matrix $A \in \mathbb{R}^{|V| \times |V|}$ defined in the following way:
$$A = (a_{ij}); \quad a_{ij}=\left\{\begin{array}{ll}
1,& \textnormal{ if } (i,j) \textnormal{ is an edge in } G,\\
0, & \textnormal{ else. }
\end{array}\right .
$$
If $G$ is weighted, then:
$$A = (a_{ij}); \quad a_{ij}=\left\{\begin{array}{ll}
w_{ij},& \textnormal{ if } (i,j) \textnormal{ is an edge in } G,\\
0, & \textnormal{ else }
\end{array}\right .
$$
where $w_{ij}$ is the weight of edge $(i,j)$.  If $G$ contains no loops, then $A$ will have zeros along the diagonal.  If $G$ is undirected, $A$ will be a symmetric matrix.

If $A$ is the adjacency matrix of a network with unweighted edges, then $[A^k]_{ij}$ counts the number of walks of length $k$ between nodes $i$ and $j$.  The diagonal entries, $[A^k]_{ii}$ count the number of closed walks of length $k$ centered at node $i$, although some of these walks can be described as ``illogical.'' For example, the walk $i \rightarrow j \rightarrow i \rightarrow j \rightarrow i$ is a closed walk of length 4 centered at node $i$ (see \cite{Die00,NetworkProp} for more details). 


\section{Basic concepts from linear algebra}
\label{sec:linear_algebra}
This section contains an overview of some basic concepts from linear algebra.  A more comprehensive guide can be found in many places, including \cite{Meyer00,Strang}.

Given a matrix $A \in \mathbb{R}^{n \times n}$, let the {\em eigenvalues} of $A$ be given by $\lambda_1, \lambda_2, \ldots, \lambda_n$ and the {\em eigenvectors} be given by ${\bf q}_1, {\bf q}_2, \ldots {\bf q}_n$ where ${\bf q}_i$ is the eigenvector associated with eigenvalue $\lambda_i$.  If the eigenvalues of $A$ are real, we label them in non-increasing order: $\lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_n$.  If $A$ has complex eigenvalues, they are labeled with non-increasing modulus: $|\lambda_1| \geq |\lambda_2| \geq \ldots \geq |\lambda_n|$. 

The {\em spectrum} of $A$ is $\sigma(A) = \{\lambda_1, \lambda_2, \ldots, \lambda_n\}$ and the {\em spectral radius} of $A$ is given by $\rho(A) = \max_i(|\lambda_i|)$. The {\em spectral gap} of $A$ is given by the difference between the two largest eigenvalues of $A$: $|\lambda_1 - \lambda_2|$.  In the case of complex eigenvalues, it is given by the difference between the moduli of the two largest eigenvalues: $|\lambda_1| - |\lambda_2|$.  The {\em range} of $A$ is given by the span of the column space of $A$: $\mathcal{R}(A) = \{v \,| \, Ax = v$ for some $x \in \mathbb{R}^n \}$.  The {\em null space} of $A$ is given by $\mathcal{N}(A) = \{ x \, | \, Ax = 0\}.$

The matrix $A$ is said to be {\em positive} (written $A > 0$) if $a_{ij} >0$ for all $1 \leq i,j \leq n$.  It is said to be {\em non-negative} (written $A \geq 0$) if $a_{ij} \geq0$ for all $1 \leq i,j \leq n$.  The adjacency matrix of a graph $G$ is non-negative if $G$ is unweighted or if all the edge weights of $G$ are positive.
\begin{definition}
\label{def:reducible}
A matrix $A$ is said to be {\em reducible} if there exists a permutation matrix $P$ such that
\begin{displaymath}\label{Dpath}
A = P^{T}AP=  \left(\begin{array}{cc}
X & Y\\
0 & Z \\  \end{array}\right) \,
\end{displaymath}
where $X$ and $Z$ are both square \cite[p. 671]{Meyer00}.  If $A$ is not reducible, it is {\em irreducible}.
\end{definition}
The adjacency matrix of an undirected graph $G$ is irreducible if and only if $G$ is connected.  The adjacency matrix of a directed graph is irreducible if and only if $G$ is strongly connected.   

A non-negative, irreducible matrix $A$ is said to be {\em primitive} if there is only one eigenvalue $\lambda$ with $|\lambda| = \rho(A)$.  Otherwise, it is {\em imprimitive}. 

If $A$ is a non-negative, irreducible matrix, the {\em Perron-Frobenius Theorem} gives useful information on the eigenvalues and eigenvectors of $A$:
\begin{theorem}
(Perron-Forbenius Theorem, \cite[p. 673]{Meyer00}) If $A \in \mathbb{R}^{n \times n}$ is non-negative and irreducible, then the following are true:
\begin{enumerate}[(i)]
	\item $r = \rho(A) \in \sigma(A)$ and $r > 0$.
	\item the algebraic multiplicity of $r$ is 1.
	\item There exists an eigenvector $ {\bf x} > 0 $ such that $A{\bf x} = r{\bf x}$.
	\item The unique vector defined by $A{\bf p} = r{\bf p}$, $p > 0$, and $\|p\|_1=1$ is called the {\em Perron vector}.  There are no non-negative eigenvectors of $A$ except for positive multiples of ${\bf p}$, regardless of the eigenvalue. 
\end{enumerate}
\end{theorem}

If $G$ is a connected, undirected graph with adjacency matrix $A$, then $\lambda_1 > \lambda_2$ by the Perron-Frobenius theorem.   Since $A$ is a symmetric, real-valued matrix, we can decompose $A$ into $A=Q\Lambda Q^T$ where $\Lambda = {\rm diag}(\lambda_1, \lambda_2, \ldots, \lambda_n)$ with $\lambda_1 > \lambda_2 \geq \ldots \geq \lambda_n$ and $Q=[{\bf q}_1, {\bf q}_2, \ldots, {\bf q}_n]$ is orthogonal. %and ${\bf q}_i$ is the eigenvector associated with $\lambda_i$.  


If $G$ is a strongly connected digraph with adjacency matrix $A$, then, by the Perron-Frobenius theorem, $\lambda_1 = r$ is a simple eigenvalue of $A$ and both the left and right eigenvectors of $A$ associated with $\lambda_1$ are positive (as $A^T$ is also irreducible).  If $G$ is also {\em diagonalizable}, then there exists an invertible matrix $X$ such that $A=X\Lambda X^{-1}$ where $\Lambda = {\rm diag}(\lambda_1, \lambda_2, \ldots, \lambda_n)$ with $\lambda_1\geq |\lambda_i|$ for $2 \leq i \leq n$, $X=[{\bf x}_1, {\bf x}_2, \ldots, {\bf x}_n]$, and $(X^{-1})^*=[{\bf y}_1, {\bf y}_2, \ldots, {\bf y}_n]$.  The left eigenvector associated with $\lambda_i$ is ${\bf y}_i$  and the right eigenvector associated with $\lambda_i$ is ${\bf x}_i$.  In the case where $G$ is not diagonalizable, $A$ can be decomposed using its Jordan Canonical form:

\begin{displaymath}\label{Dpath}
A = XJX^{-1}= X \left(\begin{array}{cc}
\lambda_1 & {\bf 0}\\
{\bf 0} & \hat{J} \\  \end{array}\right) X^{-1} \,
\end{displaymath}
where $J$ is the Jordan matrix of $A$, except that we place the $1 \times 1$ block corresponding to $\lambda_1$ first.  More information about the eigenvalues and eigenvectors associated with graphs can be found in \cite{CvRoSi97}.

An {\em M-matrix} is a real, non-singular matrix $A = (a_{ij})$ such that $a_{ij} \leq 0$ for all $i\neq j$ and $A^{-1} \geq 0$.  Many properties of non-singular M-matrices are known (see \cite{M-matrix} for a more comprehensive overview). 
\begin{theorem}
\label{thm:M-matrix}
Let $A \in \mathbb{R}^{n \times n}$ have positive diagonal entries and non-positive off-diagonal entries.  Then, the following are equivalent  \cite{M-matrix}:
\begin{enumerate}[(i)]
	\item $A$ is a non-singular M-matrix.
	\item there exists a matrix $B \geq 0$ and a real number $r > \rho(B)$ such that $A=rI-B$.
	\item the real part of $\lambda$ is positive for all $\lambda \in \sigma(A)$.
	\item $A^{-1}$ exists and is non-negative.
\end{enumerate}
\end{theorem}

We now turn to functions of a matrix.
\begin{definition}
\label{def:matrixfunct}
Given a scalar function $f$, $f$ is said to be defined on the spectrum of $A$ if the values $$f^{(j)}(\lambda_i),\,\,\,\,\, j = 0:n_i-1,\,\,\,\,\, i=1:s$$ exist, where $s$ is the number of distinct eigenvalues of $A$ and $n_i$ is the geometric multiplicity of $\lambda_i$.  These are called the values of the function $f$ on the spectrum of $A$ (\cite[p.~3]{highambook}).
\end{definition}

Suppose $f$ is defined on the spectrum of $A$.  If $A$ is symmetric or diagonalizable, then we define $f(A) = Qf(\Lambda)Q^T$ or $f(A) = Xf(\Lambda)X^{-1}$, where $$f(\Lambda) = {\rm diag}(f(\lambda_1), f(\lambda_2), \ldots, f(\lambda_n)).$$  If $A$ is not diagonalizable, then $$f(A) = \sum_{i=1}^s \sum_{j=0}^{n_i-1} \frac{f^{(j)}(\lambda_i)}{j!} (A-\lambda_iI)^j G_i$$ where $G_i$ is the oblique projector onto $\mathcal{N}((A-\lambda_iI)^{n_i})$ along $\mathcal{R}((A-\lambda_iI)^{n_i})$ \cite[p.~603]{Meyer00}.

In this thesis, we are primarily concerned with two matrix functions: the {\em matrix exponential} and the {\em matrix resolvent}.  The (parameterized) matrix exponential is given by $e^{\beta A}$, $\beta >0$.  In many applications, the parameter $\beta$ is set to one, corresponding to the unscaled matrix exponential.  The eigenvalues of $e^{\beta A}$ are given by $e^{\beta \lambda_1}, e^{\beta \lambda_2}, \ldots, e^{\beta \lambda_n}$.  The power series expansion of $e^{\beta A}$ is given by:
\begin{equation}
\label{eq:matrix_exp}
	e^{\beta A} = I + \beta A + \frac{\beta^2A^2}{2!}+ \cdots + \frac{\beta^kA^k}{k!} + \cdots = \sum_{k=0}^{\infty} \frac{\beta^kA^k}{k!}.
\end{equation}
As $[A^k]_{ij} = [A^k]_{ji}$ counts the number of walks of length $k$ between nodes $i$ and $j$, $[e^A]_{ij}$, counts the total number of  walks from node $i$ to node $j$, penalizing longer walks by scaling walks of length $k$ by the factor $\frac{\beta^k}{k!}$.

The matrix resolvent is given by $(I-\alpha A)^{-1}$, $0 < \alpha < \frac{1}{\lambda_1}$.  It has eigenvalues $\frac{1}{1-\alpha \lambda_i}$ where $\lambda_i \in \sigma(A)$. Similarly to the matrix exponential, the entries of the matrix resolvent count the number of walks in the network, penalizing longer walks.  This can be seen by considering the power series expansion of $(I-\alpha A)^{-1}$:
\begin{equation}
\label{eq:matrix_res}
(I-\alpha A)^{-1} = I + \alpha A + \alpha^2 A^2 +\cdots + \alpha^k A^k + \cdots = \sum_{k=0}^{\infty} \alpha^kA^k.
\end{equation}
Here, $[(I-\alpha A)^{-1}]_{ij}$, counts the total number of walks from node $i$ to node $j$, weighting walks of length $k$ by $\alpha^k$.  The bounds on $\alpha$ ($0 < \alpha < \frac{1}{\lambda_1}$) ensure that the matrix $I-\alpha A$ is invertible and that the power series in (\ref{eq:matrix_res}) converges to its inverse.  The bounds on $\alpha$ also force $(I-\alpha A)^{-1}$ to be nonnegative, as $I- \alpha A$ is a nonsingular $M$-matrix.

A matrix $A$ is {\em Hermitian} if it is equal to its own conjugate transpose, i.e $A^T = A^*$.  Every real, symmetric matrix is Hermitian.  A matrix is {\em upper Hessenberg} if all the entries below the first subdiagonal are zero.

\section{Properties of complex networks}
\label{sec:complex_networks}
Networks can be used to describe and analyze many different types of interactions, from those between people (social networks), to the flow of goods across an area (transportation networks), to links between websites (the WWW graph), to thermodynamic interaction between particles (statistical mechanics), and so forth.  Every network is associated with both a graph and an adjacency matrix, allowing the use of tools from graph theory and linear algebra to be used in their analysis.  Although there is no precise definition of a complex network, there are a number of properties that complex networks tend to have in common.  More information about complex networks can be found in many places, including \cite{Ebook11,New03,New10}

Complex networks are fundamentally different from random networks generated according to the Erd\"os-R\'enyi model \cite{ErRe59,ErRe60}.  There are two closely related variants of this model: in the $G(n,m)$ model, a graph is chosen uniformly at random from the set of all graphs with $n$ nodes and $m$ edges.  In the $G(n,p)$ model, a graph with $n$ nodes is constructed by randomly connecting each pair of nodes with probability $p$. If $p$ is set to $p = {n \choose 2}m$, the two models behave similarly as $n$ tends to infinity.  Networks produced by  Erd\"os-R\'enyi random models tend to have a Poisson degree-distribution.  That is, the majority of nodes tend to have a degree close to the average degree, with relatively few outliers.

However, this is not the case for complex networks.  It has been shown that many real-world networks follow a power law degree distribution, i.e. that the fraction of nodes of degree $k$ is approximately $P(k) \approx k^{-\gamma}$ where often $2 < \gamma < 3$ \cite{prefattach,FaFaFa99}.  Others have a mixed degree distribution which follows an exponential degree distribution for low-degree nodes and changes into a power law distribution once the node degree increases sufficiently  \cite{MiQi09,ZoPeLiHeCh07}.  Networks with these types of degree-distributions are often called {\em scale-free} networks. Power law degree-distributions have been observed in many types of networks including biological networks (such as protein-protein interaction and gene-regulation networks) \cite{BK12,Proteomics,Proteome}, the internet  and world-wide-web \cite{errorandattack,prefattach,FaFaFa99}, and social networks \cite{ClNeMo04}.  More information on scale-free networks can be found in \cite{Caldarelli,Gohetal}.

The {\em clustering coefficient} of a network measures the propensity of the nodes in the network to form triangles.  The clustering coefficient of a node $i$ is given by $$C_i = \frac{2|C_3(i)|}{d_i(d_i-1)},$$ where $|C_3(i)|$ is the number of triangles in which node $i$ is involved. This measures the percentage of the paths of length two (wedges) centered at node $i$ which have closed into triangles \cite{smallw}.  The clustering coefficient of the network as a whole is given by $$\bar{C} = \frac{1}{n} = \sum_{i=1}^n C_i,$$  which gives the average node clustering coefficient across the network.  A variation on the clustering coefficient, known as the {\em global clustering coefficient} or {\em network transitivity} was introduced in \cite{New01}:  $$C = \frac{3|C_3|}{|P_2|},$$ where $|C_3|$ is the number of triangles in the network and $|P_2|$ is the number of paths of length two in the network.  Networks formed using the Erd\"os-R\'enyi random models have very low clustering coefficients while many real world networks tend to exhibit high clustering coefficients.

One property common to both Erd\"os-R\'enyi random networks and many real-world networks is the {\em small-world property} \cite{smallw}.  Informally, the existence of the small-world property in a network means that all the nodes in the network are only a relatively short distance from each other, creating a ``small world'' for the nodes to exist inside.  More formally, it states that the {\em average path length}, $$\bar{l} = \frac{1}{n(n-1)} \sum_{i,j \in V} d(i,j),$$ also known as the {\em mean shortest distance} between pairs of nodes, is small.

A variety of alternative generative models for complex networks have been proposed to address the deficiencies of the Erd\"os-R\'enyi random models.  These include the preferential attachment model and variants \cite{prefattach,KrRoRe01}, the small-world model \cite{smallw}, and the Block Two-Level Erd\"os-Renyi (BTER) model \cite{BTER}.


\section{Measures of centrality and communicability}
\label{sec:centrality}
Many measures of node centrality have been developed and used over the years.  The higher the measure of node centrality, the more ``important'' a given node is considered to be in the network.  In this section, we provide a brief overview of some of the most commonly used centrality measures.  More information about various centrality measures can be found in \cite{BE06,Ebook11,New10,vigna}

One of the simplest centrality measures is {\em degree centrality} \cite{Fr79}.  In an undirected network, the degree centrality of node $i$ is given by $C_d(i) = d_i = [A{\bf 1}]_i$, where ${\bf 1}$ is the vector of all ones.  In a directed network, it is necessary to distinguish between the in-degree centrality, $C^{in}_d(i) = d^{in}_i = [{\bf 1}^TA]_i$, and the out-degree centrality, $C^{out}_d(i) = d^{out}_i = [A{\bf 1}]_i$.  Degree centrality is a local measure of node importance, counting the number of nodes in the immediate neighborhood of a given node.

The {\em eigenvector centrality} of node $i$ in a directed network is given by $C_{ev}(i) = {\bf q}_1(i)$ \cite{Bonacich}.   It gives the limit as $k$ goes to infinity of the percentage of walks of length $k$ which start at node $i$ \cite{CvRoSi97}.  Thus, in contrast to degree-centrality, eigenvector centrality measures the global influence of nodes in the network.  Variations on eigenvector centrality have been proposed to measure slightly different aspects of node importance \cite{BL01}.  In a directed network, the eigenvector centrality can be defined by both the dominant left eigenvector and the dominant right eigenvector.  These can be used to rank nodes as authorities and hubs, respectively (see section \ref{sec:HITS} for definitions of hubs and authorities).

The {\em closeness centrality} \cite{Fr79} of node $i$ is given by $$CC(i) = \frac{n-1}{\sum_{j \neq i} d(i,j)}.$$  It measures the overall closeness of node $i$ to the rest of the network.  

Another commonly used centrality measure is {\em betweenness centrality} \cite{Fr77}.  The betweenness centrality of node $i$ is given by $$ BC(i) = \sum_{j\neq i \neq k} \frac{\sigma_{jk}(i)}{\sigma_{jk}}$$ where $\sigma_{jk}$ is the total number of shortest paths between nodes $j$ and $k$ and $\sigma_{jk}(i)$ is the number of those shortest paths which pass through node $i$.  Thus, the betweenness centrality of node $i$ measures how important $i$ is for fast communication across the network as a whole. Variations of betweenness centrality have also been proposed \cite{Br08}.

The matrix resolvent $(I-\alpha A)^{-1}$, $0 < \alpha < \frac{1}{\lambda_1}$, was first used to rank nodes in a network in the early 1950s, when Katz  used the column sums  to calculate node importance \cite{Katz}.  The Katz centrality score of node $i$ is given by $K_i(\alpha) =[(I-\alpha A)^{-1} {\bf 1}]_{i}$.  This counts the total number of walks between node $i$ and every other node in the network, weighting walks of length $k$ by a penalization factor of $\alpha^k$.  Closely related to Katz centrality is {\em resolvent subgraph centrality} \cite{NetworkProp,estradarodriguez05}.  The resolvent subgraph centrality of node $i$ is given by $RC_i(\alpha) = [(I-\alpha A)^{-1}]_{ii}$.  This counts the number of closed walks centered at node $i$, again weighting walks of length $k$ by $\alpha^k$.

A more commonly used form of subgraph centrality is {\em exponential subgraph centrality} \cite{NetworkProp,estradarodriguez05}.  The exponential subgraph centrality of node $i$ is given by $SC_i(\beta) = [e^{\beta A}]_{ii}$.  This counts the number of closed walks centered at node $i$, weighting walks of length $k$ by $\frac{\beta^k}{k!}$.  In Chapter \ref{C:TotalComm}, we introduce a new centrality measure based on exponential subgraph communicability \cite{BKlym13,estradahatano2,EHB11}.

Many centrality measures were originally developed for use on undirected networks and later needed to be extended for use on directed networks.  In Sections \ref{sec:HITS} and \ref{sec:PageRank}, we describe HITS and PageRank, two rankings methods that were developed specifically for directed networks.  Both HITS and PageRank were originally developed as web search algorithms.  However, as they provide a ranking of the nodes in directed networks, they can also be used as centrality measures in general digraphs.

The {\em communicability} between a pair of nodes $i$ and $j$ measures how well the two nodes can exchange information with each other.  A (relatively) large communicability between a pair of nodes $i$ and $j$ indicates that information flows more easily between those two nodes than between pairs of nodes with lower communicability.  In other words, a low communicability indicates that the two nodes cannot easily exchange information.  Network communicability can also be interpreted in terms of the correlations between different components of physical systems; see, e.g., \cite{EHB11}. 

Two of the most commonly used communicability measures are {\em (exponential) subgraph centrality} and {\em resolvent subgraph centrality}  \cite{estradahatano2,NetworkProp}.  The subgraph communicability between nodes $i$ and $j$ is given by $[e^{\beta A}]_{ij}$, $\beta > 0$ (note that in the case of an undirected network, $A$ is symmetric and $[e^{\beta A}]_{ij}=[e^{\beta A}]_{ji}$).  This counts the total number of walks between nodes $i$ and $j$, weighting walks of length $k$ by $\frac{1}{k!}$.  The resolvent subgraph communicability between nodes $i$ and $j$ is given by $[(I- \alpha A)^{-1}]_{ij}$, $0 < \alpha < \frac{1}{\lambda_1}$.  As in the case of exponential subgraph centrality, this counts the total number of walks between nodes $i$ and $j$, but here walks of length $k$ are weighted by $\alpha^k$.


\subsection{HITS}
\label{sec:HITS}
The classical {\em Hypertext Induced Topics Search} (HITS)
algorithm, first introduced by J.~Kleinberg in \cite{HITS} is based on the idea that in the World Wide Web, and indeed
in all document collections which can be represented by directed networks, 
there are two types of important nodes: {\em hubs} and {\em authorities}.  
Hubs are nodes which point to many nodes of the type considered important.  Authorities are these important nodes.  From this comes a circular definition: good hubs are those which point to many good authorities and good authorities are those pointed to by many good hubs.

Thus, the HITS ranking relies on an iterative method converging to a stationary solution.  Each node $i$ in the network is assigned two non-negative weights, an {\em authority weight} $x_i$ and a {\em hub weight} $y_i$.  To begin with, each $x_i$ and $y_i$ is given
an arbitrary nonzero value.  Then, the weights are updated in the following ways:
\begin{equation}
x^{(k)}_i = \sum_{j:(j,i) \in E} y^{(k-1)}_j {\,\,\,\,\,\ \rm and \,\,\,\,\,\,} y^{(k)}_i = \sum_{j:(i,j) \in E} x^{(k)}_j {\rm \,\,\,\,\, for \,\,\,\,\,} k=1,2,3...
\end{equation}
The weights are then normalized so that $\sum_j (x_j^{(k)})^2 =1$ and 
$\sum_j (y_j^{(k)})^2 =1$.

The above iterations occur sequentially and it can be shown that, under
mild conditions, both sequences of 
vectors $\{{\bf x}^{(k)}\}$ and $\{{\bf y}^{(k)}\}$ converge as $k \to \infty$.  
In practice, the iterative process is continued 
until there is
no significant change between consecutive iterates.

This iteration sequence shows the natural dependence relationship between hubs and authorities: if a node $i$ points to many nodes with large $x$-values, it receives a large $y$-value and, if it is pointed to by many nodes with large $y$-values, it receives a large $x$-value.

In terms of matrices, the equation (3.1) becomes: 
${\bf x}^{(k)} = A^T{\bf y}^{(k-1)}$ and ${\bf y}^{(k)} = A{\bf x}^{(k)}$, followed by
normalization in the 2-norm.
This iterative process can be expressed as 
\begin{equation}
{\bf x}^{(k)} = c_kA^TA{\bf x}^{(k-1)} {\rm \,\,\,\,\, and \,\,\,\,\,} {\bf y}^{(k)} = c_k'AA^T{\bf y}^{(k-1)},
\end{equation}
where $c_k$ and $c_k'$ are normalization factors.
A typical choice for the initizialization vectors ${\bf x}^{(0)}$, ${\bf y}^{(0)}$ 
would be the constant vector
$${\bf x}^{(0)} = {\bf y}^{(0)} = [1/\sqrt{n},\, 1/\sqrt{n},\, \ldots ,1/\sqrt{n}],$$ 
see \cite{EIHITS}.
Hence, HITS is just an iterative power method to 
compute the dominant eigenvector for $A^TA$ and for $AA^T$.  
The authority scores are determined by the 
entries of the dominat eigenvector of the matrix $A^TA$, 
which is called the {\em authority matrix} and the hub 
scores are determined by the entries of the
dominant eigenvector of $AA^T$, called the {\em hub matrix}. Recall that 
the eigenvalues of  both $A^TA$ and $AA^T$ are the squares of the 
singular values of $A$. Also, the eigenvectors of $A^TA$ are the right singular
vectors of $A$, and the eigenvectors of $AA^T$ are the left singular 
vectors of $A$.  See Algorithm \ref{alg:HITS} in Appendix \ref{A:algorithms} for a pseudocode implementation of the HITS algorithm.

A number of variations on the HITS algorithm have been proposed, for example \cite{BEK,LinkStructure,Dingetal,IRSurvey}


\subsection{PageRank}
\label{sec:PageRank}
The {\em PageRank} algorithm \cite{BrPa98} is perhaps the most famous (and lucrative) algorithm in the world.  It was developed to rank authoritative webpages in search engine results.  Although PageRank, and its variant Reverse PageRank, described below, were originally developed to rank webpages, they have since been used to rank nodes in many diverse types of networks \cite{BK12,Bianchini05,franceschet,KG11,MoBrHiGi05,WuLi}.

In the PageRank algorithm, the importance of a webpage is determined by ``votes,'' i.e., the number of links to that webpage. The main idea is that votes from more important sites count more than votes from less important sites and the significance of any vote is weighed by how many votes the ``voting'' webpage makes.  At its most basic, given a directed network, the PageRank of node $i$ is given by $$PR_i = \sum_{\{j \,|\, (j,i) \in E\}} \frac{PR_j}{d^{out}_j}.$$  As the PageRank of a node depends on the PageRank of its neighbors, the rankings are calculated in an iterative manner.  To begin, each node is given an initial ranking.  Often, $PR_i^{(0)} = \frac{1}{n}$.  Then, 
\begin{equation}
\label{eq:PR1}
PR^{(k)}_i = \sum_{\{j \,|\, (j,i) \in E\}} \frac{PR^{(k-1)}_j}{d^{out}_j} {\rm \,\,\,\,\, for \,\,\,\,\,} k=1,2,3...
\end{equation}
Setting $\boldsymbol\pi^{(k)T} = (PR_1^{(k)}, PR_2^{(k)}, \ldots, PR_n^{(k)})$, (\ref{eq:PR1}) becomes
\begin{equation}
\label{eq:PR2}
\boldsymbol\pi^{(k)T} = \boldsymbol\pi^{(k-1)T} P
\end{equation}
where $P = \left\{ \begin{array}{lr}
					\frac{1}{|d^{out}_i|}, {\rm if \,} (i,j) \in E  \\
					0, else.
				\end{array}
			\right.$

Equation (\ref{eq:PR2}) almost corresponds to the power method applied to a finite state space {\em Markov chain} with {\em transition probability matrix} $P$.  A Markov chain is a stochastic process that satisfies the {\em Markov property} \cite[p.~687]{Meyer00}.  That is, given a family of random variables $\{X_t\}$ in which every $X_t$ has the same range $\{S_1, S_2, \ldots, S_n\}$ (called the {\em state space}), the following holds: $$P(X_{t+1} =S_j \,|\, X_t=S_{i_t},X_{t-1} = S_{i_{t-1}}, \ldots, X_0 = S_{i_0}) = P(X_{t+1} = S_j \,|\,X_t = S_{i_t}). $$  This implies that the probability of entering state $S_j$ at step $t+1$ only depends on the state at time $t$ and not on previous states.  If the time parameter is discrete, the Markov chain is a {\em discrete-time Markov chain}. 

The {\em transition probability} $p_{ij}(t) =  P(X_{t+1} = S_j \,|\,X_t = S_{i_t})$ measures the probability of moving from state $S_i$ to $S_j$ at time $t$.  The {\em transition probability matrix} of a discrete-time Markov chain is given by $P(t) \in \mathbb{R}^{n \times n}$ with $P(t)= (p_{ij}(t))$.  If the transition probabilities are fixed, the Markov chain is said to be {\em homogenous} and $P(t) =P$.  The transition probability matrix $P$ is a {\em (row-)stochastic} matrix.  That is, $P \geq 0$ and $\sum_{j=1}^n p_{ij} = 1$ for all $1 \leq i \leq n$.

A Markov chain is {\em irreducible} if the transition probability matrix is irreducible.  Equivalently, a Markov chain is irreducible if every state can be reached from every other state.  The following theorem presents results on the convergence of irreducible Markov chains \cite[p.~693]{Meyer00}:
\begin{theorem}
\label{thm:Markov}
Let $P$ be the transition matrix for an irreducible Markov chain on states $\{S_1, S_2, \ldots, S_n\}$, i.e. $P \in  \mathbb{R}^{n \times n}$ is an irreducible stochastic matrix.  Let $\boldsymbol\pi^T$ denote the left-hand Perron vector of $P$.  The following statements are true for every initial distribution $\boldsymbol\pi^{(0)T}$:
\begin{enumerate}[(i)]
	\item The $k$th step transition matrix is $P^k$, i.e, $[P^k]_{ij}$ is the probability of moving from $S_i$ to $S_j$ in exactly $k$ steps.
	\item The $k$th step transition vector is given by $\boldsymbol\pi^{(k)T} = \boldsymbol\pi^{(0)T}P^k.$
	\item If $P$ is primitive, then $$\lim_{k \to \infty} P^k = {\bf 1}\boldsymbol\pi^T \,\,\,\,\, {\rm and} \,\,\,\,\, \lim_{k \to \infty}\boldsymbol\pi^{(k)T} = \boldsymbol\pi^T. $$
	\item If $P$ is imprimitive, then  $$\lim_{k \to \infty} \frac{I+P+\ldots +P^{k-1}}{k} = {\bf 1}\boldsymbol\pi^T \,\,\,\,\, {\rm and} \,\,\,\,\, \lim_{k \to \infty}\frac{\boldsymbol\pi^{(0)T}+\boldsymbol\pi^{(1)T}+ \ldots + \boldsymbol\pi^{(k-1)T}}{k}   = \boldsymbol\pi^T. $$
	\item For both the primitive and imprimitive cases, $\pi_j$, the $j$th component of  $\boldsymbol\pi$, represents the long-run fraction of time that the chain is in $S_j$.
	\item $\boldsymbol\pi^T$ is often called the {\em stationary distribution vector} for the chain because it is the unique distribution vector satisfying $\boldsymbol\pi^TP = \boldsymbol\pi^T$
\end{enumerate}
\end{theorem}

Equation (\ref{eq:PR2}) is not quite the power method applied to a finite state space Markov because the matrix $P$ used is not necessarily stochastic.  If node $i$ has no out-edges, then the $i$th row of $P$ will be all zeros.  This is addressed by replacing zero rows with $\frac{e^T}{n}$, forming the matrix $\bar{P}$.  In terms of web surfing, this corresponds to jumping to a random web-page when a page with no out-links is reached.  However,  $\bar{P}$ is still not guaranteed to be an irreducible matrix and, thus, the results from Theorem \ref{thm:Markov} do not guarantee converge to a unique stationary vector.  This is addressed in the PageRank algorithm by forming a matrix $E = \frac{1}{n}{\bf 1}{\bf 1}^T$ and setting 
\begin{equation}
\label{PR3}
\bar{\bar{P}} = \alpha\bar{P} + (1-\alpha)E \,\,\,\,\,\, {\rm for \,\,\,} 0 < \alpha <1.
\end{equation}
Often, the values $\alpha =0.85$ is used.  The matrix $\bar{\bar{P}}$, often referred to as the ``Google'' matrix, is stochastic and irreducible.  Therefore, by Theorem \ref{thm:Markov}, it has a unique stationary distribution vector $\boldsymbol\pi$.  The PageRank centrality of node $i$ is given by $\pi_i$, the $i$th entry of $\boldsymbol\pi$.  See Algorithm \ref{alg:PageRank} in Appendix \ref{A:algorithms} for a pseudocode implementation of the PageRank algorithm.

It was pointed out in \cite{fogaras} that applying PageRank to the digraph obtained by reversing the direction of the edges provides a natural way to rank the hubs; this is usually referred to as {\em Reverse PageRank}. In other words, authority rankings are obtained by applying PageRank to the \lq\lq Google\rq\rq\ matrix derived from $A$, and hub rankings are obtained by the same process applied to $A^T$.

A more comprehensive overview of the PageRank algorithm can be found in \cite{Pagerank}.  Many studies have been done to analyze various aspects of the PageRank algorithm \cite{Bar-Yossef,franceschet,Gleich07,Ipsen06,MiPapa,Whos}.


\section{Comparing centrality measures}
\label{sec:comparing}
There are many way to compare two sets of ranked lists.  See \cite{Whos} for an overview.  In this thesis, we use (Pearson) correlation coefficients and the intersection distance method (see \cite{Faginetal} as well as \cite{Boldi,Randomalpha}) on both the full set of nodes of a network $G$ and on partial lists of nodes to measure similarities between the rankings obtained with two methods.  The correlation coefficients are computed using lists of nodes in rank order.  The intersection distances are computed using the lists of centrality values.  

Given two ranked lists $x$ and $y$, the intersection distance between the two lists is computed in the following way: let $x_k$ and $y_k$ be the top $k$ ranked items in $x$ and $y$ respectively.  Then the top $k$ intersection distance (or {\em intersection similarity}) is given by 
\begin{equation}
\label{eq:isim}
{\rm isim}_k(x,y) := \frac{1}{k}\sum_{i=1}^k \frac{|x_i\Delta y_i|}{2i}
\end{equation}
where $\Delta$ is the symmetric difference operator between the two sets.  If the lists are identical, then ${\rm isim}_k(x,y) = 0$ for all $k$.  If the two sequences are disjoint, then ${\rm isim}_k = 1$.  Thus, a small value of ${\rm isim}_k(x,y)$ indicates a strong similarity between the two rankings.

Throughout the thesis, we denote by $\mathpzc{cc}$ the correlation coefficient between the two vector rankings, and by $\mathpzc{cc}_p$ the correlation coefficient between the top $p$\% of nodes under  two ranking systems. We denote by ${\rm isim}$ the intersection distance between two complete sets of rankings and by ${\rm isim}_{p\%}$ the intersection distance between the top $p\%$ of nodes. 


\section{Approximations of centrality scores}
\label{sec:iterative_methods}
In this section, we will provide an overview of two iterative methods used in this thesis to approximate centrality measures based on matrix functions.  The matrix exponential and resolvent are computationally expensive and exact computations are not possible for large networks.  Several approaches are available for computing the matrix exponential \cite{highambook}. A commonly used scheme is the one based on Pad\'e approximation combined with the scaling and squaring method \cite{higham,highambook}, implemented in Matlab by the
{\tt expm} function. For an $n\times n$ matrix, this method requires $O(n^2)$ storage and $O(n^3)$ arithmetic operations; current implementations are geared toward small, dense matrices. Evaluation of the matrix exponential based on diagonalization also requires$O(n^2)$ storage and $O(n^3)$ operations. Furthermore, these methods cannot be easily adapted to the case where only selected entries (e.g., the diagonal ones) of the matrix exponential are of interest.


\subsection{Approximation of diagonal entries of matrix functions}
To compute the exponential (resolvent) subgraph centralities of the nodes in a network, only the diagonal entries of the matrix exponential (resolvent) are needed.  Efficient, accurate methods for estimating (or, in some cases, bounding) 
arbitrary entries in
a matrix function $f(A)$ have been developed by Golub, Meurant and collaborators
(see \cite{golubmeurantbook} and references therein) and first applied
to problems of network analysis by Benzi and Boito in \cite{BenziBoito};
see also \cite{bonchietal}.  They have been implemented in 
the Matlab toolbox {\tt mmq}  \cite{mmq}.
Here we provide a brief description of these methods.  Additional details can be found in
 \cite{BenziBoito} and \cite{golubmeurantbook}. 
 
 Let $A$ be a real, symmetric, $n\times n$ matrix
and let $f$ be a function defined on the spectrum of $A$.
 Consider the eigendecomposition
$A=Q \Lambda Q^T$ and $f(A)=Qf(\Lambda)Q^T$, where
$Q=[{\bf q}_1, \ldots, {\bf q}_n]$ and 
$\Lambda = {\rm diag}\,(\lambda_1,\ldots ,\lambda_n)$.  For notational simplicity, here we assume
that the eigenvalues of $A$ are ordered as
$\lambda_1 \le \ldots \le \lambda_n$. For given vectors ${\bf u}$ and ${\bf v}$ we have
\begin{equation}
%u^T f(A) v=u^T Q f(\Lambda) Q^T v=w^Tf(\Lambda)z=\sum_{i=1}^n f(\lambda_i)w(i)z(i),
{\bf u}^T f(A) {\bf v}={\bf u}^T Q f(\Lambda) Q^T {\bf v}={\bf w}^Tf(\Lambda){\bf z}=\sum_{k=1}^n f(\lambda_k)w_kz_k,
\end{equation}
where ${\bf w}=Q^T{\bf u}=(w_k)$ and ${\bf z}=Q^T {\bf v}=(z_k)$.
In particular, for $f(A)=e^A$ we obtain
\begin{equation}
{\bf u}^T e^A {\bf v}= \sum_{k=1}^n e^{\lambda_k}w_kz_k.
\label{exp3}
\end{equation}
Choosing ${\bf u}={\bf v}={\bf e}_i$ (the vector with the $i$th entry equal to 1 and all the
remaining ones equal to $0$) we obtain an expression for the subgraph
centrality of node $i$:
$$ SC(i) := \sum_{k=1}^n e^{\lambda_k} {\bf q}_{k}(i)^2\,,$$
where ${\bf q}_{k}(i)$ denotes the $i$th component of vector ${\bf q}_k$.
Likewise, choosing ${\bf u}={\bf e}_i$ and ${\bf v}={\bf e}_j$ we obtain the following expression
for the communicability between node $i$ and node $j$:
$$ C(i,j) := \sum_{k=1}^n e^{\lambda_k}{\bf q}_{k}(i){\bf q}_{k}(j).$$
Analogous expressions hold for other matrix functions, such as the resolvent.
Hence, the problem is reduced to evaluating bilinear
expressions of the form ${\bf u}^T f(A) {\bf v}$. 
Such bilinear forms can be thought of as Riemann- Stieltjes integrals
with respect to a (signed) spectral measure:
$$
{\bf u}^T f(A) {\bf v}=\int_a^bf(\lambda)d\mu(\lambda),\quad
\mu(\lambda)=\left\{\begin{array}{ll}
0,& \textnormal{ if }\lambda<a=\lambda_1,\\
\sum_{k=1}^iw_kz_k, & \textnormal{ if }\lambda_{i}\leq\lambda < \lambda_{i+1},\\
\sum_{k=1}^n w_kz_k, & \textnormal{ if }b=\lambda_n\leq\lambda.
\end{array}\right .
$$
This integral can be approximated by means of a Gauss-type quadrature rule:
\begin{equation}
\int_a^bf(\lambda)d\mu(\lambda)=\sum_{j=1}^pc_jf(t_j)+\sum_{k=1}^qv_kf(\tau_k)+R[f],
\label{quadrature}
\end{equation}
where $R[f]$ denotes the error. Here the nodes $\{t_j\}_{j=1}^p$ and the 
weights $\{c_j\}_{j=1}^p$ are unknown, whereas
the nodes $\{\tau_k\}_{k=1}^q$ are prescribed. We have

\begin{itemize}
\item $q=0$ for the Gauss rule,
\item $q=1$, $\tau_1=a$ or $\tau_1=b$ for the Gauss--Radau rule,
\item $q=2$, $\tau_1=a$ and $\tau_2=b$ for the Gauss--Lobatto rule.
\end{itemize}

For certain matrix functions, including the exponential and the resolvent,
these quadrature rules can be used to obtain lower and upper bounds on the
quantities of interest; prescribing additional 
quadrature nodes leads to tighter and tighter
bounds, which (in exact arithmetic) converge monotonically to the true values 
\cite{golubmeurantbook}.
The evaluation of these quadrature rules is mathematically equivalent to the
computation of orthogonal polynomials via a three-term recurrence, or,
equivalently, to the computation of entries and spectral information of a certain
tridiagonal matrix via the Lanczos
algorithm (see Algorithm \ref{alg:Lanczos} in Appendix \ref{A:algorithms}). More information on the Lanczos algorithm can be found in \cite{Me06}. Here we briefly recall how this can be done for the case of
the Gauss quadrature rule, when we wish to estimate
the $i$th diagonal entry of $f(A)$.
It follows from (\ref{quadrature}) that the quantity
of interest has the form $\sum_{j=1}^pc_jf(t_j)$.
This can be computed from the relation \cite[p.~28]{golubmeurantbook}:
$$
\sum_{j=1}^pc_jf(t_j)={\bf e}_1^Tf(J_p){\bf e}_1,
%\sum_{j=1}^nc_jf(t_j)=\langle e_1\,|\, f(J_n)\, |\, e_1\rangle,
$$
where
\begin{displaymath}
J_p=\left(\begin{array}{ccccc}
\omega_1 & \gamma_1 &&&\\
\gamma_1 & \omega_2 & \gamma_2 &&\\
& \ddots & \ddots & \ddots &\\
&& \gamma_{p-2} & \omega_{p-1} & \gamma_{p-1}\\
&&& \gamma_{p-1} & \omega_p \end{array}\right)
\end{displaymath}
is a tridiagonal matrix whose eigenvalues are the
Gauss nodes, whereas the Gauss weights are given
by the squares of the first
entries of the normalized eigenvectors of $J_p$.
The entries of $J_p$ are computed using the Lanczos algorithm
with starting vectors
${\bf x}_{-1}={\bf 0}$ and ${\bf x}_0={\bf e}_i$. 
Note that it is not required to compute all the
components of the eigenvectors of $J_p$ if one
uses the Golub--Welsch QR algorithm; 
see \cite{golubmeurantbook}.

For small $p$ (i.e., for a small number of Lanczos steps), computing the
$(1,1)$ entry of $f(J_p)$ is inexpensive.
The main cost in estimating one entry of $f(A)$ with this approach
is associated with the sparse matrix-vector multiplies in the Lanczos
algorithm applied to the
adjacency matrix $A$. If only a small, fixed number of iterations are performed
for each diagonal element of $f(A)$, as is usually the case,
the computational cost (per node) is
at most $O(n)$ for a sparse graph, resulting
in a total cost of $O(n^2)$ for computing the subgraph centrality
of every node in the network. If only $k< n$ subgraph centralities 
are wanted, with $k$ independent of $n$, then
the overall cost of the computation will be $O(n)$ provided   
that sparsity is carefully exploited in the Lanczos algorithm
and that only a small number $p$ of iterations (independent of $n$)
is carried out. Note, however, that depending on the connectivity
characteristics of the network under consideration, the prefactor
in the $O(n)$ estimate could be large. The algorithm can be
implemented so that the storage requirements are $O(n)$ for a
sparse network---that is, a network in which the total number of links
grows linearly in the number $n$ of nodes.


Additionally, in most applications one is not so much interested in computing an exact ranking of {\em all} the nodes in a network, but only in identifying the top $k$ ranked nodes, where the integer $k$ is small compared to $n$ (for example, $k=10$ or $k=20$). Research has been done to develop methods that are capable of quickly identifying the top $k$ nodes without having to compute accurate subgraph centrality scores for each node \cite{BEK,fmrr12}.  We return to this topic in Section \ref{sec:approx}.

\subsection{Approximation of row sums of matrix functions}
\label{sec:approx_rowsums}
To compute Katz centrality or total communicability (which is introduced in Chapter \ref{C:TotalComm}), all that is needed is the row sums of the matrix resolvent or exponential: $(I-\alpha A)^{-1}{\bf 1}$ or $e^{\beta A}{\bf 1}$.  For digraphs, $A^T$ may be used in place of $A$.  In recent years, efficient algorithms have been developed for
computing the {\em action} of a matrix function on a vector, that is,
for computing the vector $f(A){\bf b}$ for a given matrix $A$ (usually large
and sparse), vector $\bf b$, and function $f$. A particularly important
case is that of the matrix exponential, since this provides a
solution method for initial value problems for first-order systems 
of linear ordinary differential equations. These algorithms, based
on variants of the Lanczos, Arnoldi or other Krylov subspace method,
access the matrix $A$ only in the form of (sparse) matrix-vector
products and have $O(n)$ storage cost for a sparse $n\times n$ matrix $A$
\cite[Chapter 13]{highambook}. When ${\bf b} = {\bf 1}$, the
vector with all its entries equal to 1, the $i$th entry
of the resulting
vector $f(A){\bf 1}$ contains the $i$th row sum of $f(A)$:
$$\left [ f(A){\bf 1} \right ]_i = \sum_{j=1}^n \left [ f(A) \right ]_{ij},
\quad 1\le i \le n\,.$$
This quantity  can be computed much faster than the diagonal entries of a matrix function using current
computational techniques.  Indeed, computing individual entries of matrix 
functions $f(A)$ is generally costly for large $A$ even with the 
best available algorithms \cite{BenziBoito,EHB11}. Of course, the same is true if the vector
${\bf 1}$ is replaced by some other vector---typically, an \lq\lq external
importance vector\rq\rq\ which can be used to take into account intrinsic,
not network-related contributions to the centrality of each node 
\cite[pp.~174--175]{New10}.

An efficient algorithm for evaluating $f(A){\bf v}$ using a restarted Krylov method has been presented in \cite{Krylov1,Krylov2}.  In this approach, the basic operation is represented by matrix-vector products with $A$.  This method has been implemented in the Matlab toolbox {\tt  funm\_kryl} 
by Stefan G\"uttel \cite{funmkryl}.

Given $A \in \mathbb{R}^{n \times n}$, ${\bf b} \in \mathbb{R}^n$, and $f$ analytic on a neighborhood containing $\sigma(A)$, the evaluation of 
\begin{equation}
\label{eq:f(A)}
f(A){\bf b}
\end{equation}
using Krylov subspace methods is based on an Arnoldi-like decomposition of  $A$ into
\begin{equation}
\label{eq:Krylov1}
 AV_m = V_{m+1} \tilde{H}_m = V_mH_m+h_{(m+1,m)}{\bf v}_{m+1}{\bf e}_m^T
\end{equation}
where $\tilde{H} = h_{ij}$ is an $(m+1) \times m$ upper Hessenberg matrix, $H_m = [I_m \,\, {\bf 0}]\tilde{H}_m$, and ${\bf e}_m \in \mathbb{R}^m$ is the $m$th unit vector.  The columns of $V_m$ form a basis of the Krylov subspace:
\begin{equation}
\label{eq:Krylov2}
\mathcal{K}_m(A, {\bf b}) = {\rm span}\{{\bf b}, A{\bf b}, \ldots, A^{m-1}{\bf b}\}.
\end{equation}
If $A$ is Hermitian, then $H_m$ in (\ref{eq:Krylov1}) is Hermitian tridiagonal and the basis $\{{\bf v}_m\}$ of $\mathcal{K}_m(A, {\bf b})$ is orthogonal.  In this case, the Arnoldi approximation of (\ref{eq:Krylov1}) is given by:
\begin{equation}
\label{eq:Krylov3}
f(A){\bf b} \approx V_m f(H_m) V_m^*{\bf b} = V_mf(H_m){\bf e}_1.
\end{equation}

The evaluation of (\ref{eq:Krylov3}) requires that all $m$ basis vectors of $\mathcal{K}_m(A, {\bf b})$ be available, which can be computationally infeasible for large $A$ and/or $m$.  A restarted Krylov subspace algorithm repeatedly generates a basis of fixed dimension $\hat{m}$, producing a sequence of approximations to $f(A){\bf b}$.  At the end of each iteration, all basis vectors but the last, ${\bf v}_{\hat{m}}$, are discarded and  ${\bf v}_{\hat{m}}$ becomes the initial vector for the next Krylov subspace.

For many network of practical interest, the cost is typically
$O(n)$, although the prefactor can vary considerably
for different types of networks.  A more comprehensive overviews of Krylov subspace methods for approximating $f(A){\bf b}$ can be found in \cite[Ch 13]{highambook}.  Additional uses of Krylov subspace methods can be found in \cite[Ch. 7]{Meyer00} and \cite[Ch. 9 and 10]{GolVan96}.

%: ----------------------- HELP: adding figures with macros
% \figuremacro{1}{2}{3}{4} calls up a series of commands formating your image.
% 1 = name of the file without extension; PNG, JPEG is ok; GIF doesn't work
% 2 = title of the figure AND the name of the label for cross-linking
% 3 = caption text for the figure
% 4 = size relative to text width which is 1; use this to reduce figures




