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

%: ----------------------- introduction file header -----------------------
\chapter{Ranking Hubs and Authorities using Matrix Functions}\label{C:HubsAuths}

\section{Introduction}
As we have seen previously, subgraph centrality is one of the most commonly used centrality measures.  The interpretation of centrality described in \cite{NetworkProp} applies mostly to 
undirected networks. However, many important real-world networks 
(the World Wide Web, the Internet, citation networks, food webs, 
certain social networks, etc.)~are
directed. One goal of this chapter is to extend the notions of centrality 
and communicability described in 
\cite{estradahatano2,NetworkProp} to directed networks, 
with an eye towards developing new %hubs and authorities 
ranking algorithms for, e.g., document collections, web pages, and so forth.
We further compare our approach with some standard algorithms, such as HITS
(see \cite{HITS}) and a few others. 
Methods of quickly determining hub and authority 
rankings using Gauss-type quadrature rules are also 
discussed.

%The results from this chapter were published in \cite{BEK}.

\section{HITS reformulation}
\label{sec:newHITS}
One of the classical methods for ranking nodes in a directed network is the Hypertext Induced Topics Search (HITS)
algorithm, first introduced by J.~Kleinberg in \cite{HITS}, which ranks nodes as hubs and authorities.  This algorithm provides the motivation for the extension of subgraph centrality to directed graphs given in section \ref{ExtDig}.  The standard presentation of the HITS algorithm can be found in Chapter \ref{C:introduction}, Section \ref{sec:iterative_methods}.

In a digraph the adjacency matrix $A$ is generally nonsymmetric, 
however, the two matrices used in the HITS algorithm ($A^TA$ and $AA^T$) are symmetric.  Note that, setting
$$ \mathcal{A} =  \left( \begin{array}{cc}
0 & A \\
A^T & 0  \end{array} \right)\,,$$
a symmetric matrix is obtained.  Now, 
$$\mathcal{A}^2 =   \left( \begin{array}{cc} AA^T & 0 \\ 0 & A^TA  \end{array} \right); {\rm \,\,\,\,\,\,\,} \mathcal{A}^3 =   \left( \begin{array}{cc} 0 & AA^TA \\ A^TAA^T & 0  \end{array} \right). $$
In general, 
$$\mathcal{A}^{2k} =   \left( \begin{array}{cc} (AA^T)^k & 0 \\ 0 & (A^TA)^k  \end{array} \right); {\rm \,\,\,\,\,\,\,} \mathcal{A}^{2k+1} =   \left( \begin{array}{cc} 0 & A(A^TA)^k \\ (A^TA)^kA^T & 0  \end{array} \right). $$

Applying HITS to this matrix $\mathcal{A}$, $\mathcal{A}^T = \mathcal{A}$ 
so $\mathcal{A}^T\mathcal{A} = \mathcal{A}\mathcal{A}^T = \mathcal{A}^2$ 
and introducing the vector ${\bf u}^{(k)} = \left( \begin{array}{c} {\bf y}^{(k)}\\{\bf x}^{(k)}
\end{array} \right)$
for $k=1,2,3,\ldots$, equation (3.2) becomes
\begin{equation}
{\bf u}^{(k)} = \mathcal{A}^2 {\bf u}^{(k-1)} =  \left( \begin{array}{cc} AA^T & 0 \\ 0 & A^TA  
\end{array} \right) {\bf u}^{(k-1)},
\end{equation}
followed by normalization of the two vector components of ${\bf u}^{(k)}$ so that 
each has 2-norm equal to 1.
Now, if $A$ is an $n \times n$ matrix, $\mathcal{A}$ 
is $2n \times 2n$ and vector ${\bf u}^{(k)}$ is in 
$\mathbb{R}^{2n}$.  The first $n$ entries of ${\bf u}^{(k)}$ 
correspond to the hub rankings of the nodes, while 
the last $n$ entries give the authority rankings.
Under suitable assumptions (see the discussion in \cite[Chapter 11.3]{Pagerank}),
as $k\to \infty$ the
sequence $\{{\bf u}^{(k)}\}$ converges to the dominant nonnegative
eigenvector of $\cal A$, which yields the desired hub and
authority rankings. 

Hence, in HITS only information obtained from the dominant
eigenvector of $\cal A$ is used. It is natural to expect
that taking into account spectral information corresponding
to the remaining eigenvalues and eigenvectors of $\cal A$ 
may lead to improved results.

Among the limitations of HITS, we mention the possible dependence
of the rankings on the choice of the initial vectors ${\bf x}^{(0)}$,
${\bf y}^{(0)}$, see \cite{EIHITS} for examples of this.


\section{Subgraph centralities and communicabilities}\label{SC}
In \cite{NetworkProp}, the authors review the use of (exponential) subgraph centrality and communicability for ranking nodes in undirected networks.  However, these measures are less meaningful when applied to directed networks.  Although the matrix exponential is certainly well-defined for any
matrix, whether symmetric or not, the {\em interpretation} of the   
notion of subgraph centrality 
for directed networks can be problematic. 
To see this, consider the directed path graph consisting of $n$ nodes,
with edge set $E = \{(1,2)\,, (2,3)\,, \ldots , (n-1,n)\}$ and
adjacency matrix
%\begin{displaymath}\label{AP}
\begin{equation}\label{AP}
A=\left(\begin{array}{ccccc}
0 & 1 & 0 &\cdots & 0\\
0 & 0 & 1 & \cdots & 0\\
\vdots & \ddots & \ddots & \ddots &\vdots \\
0 & 0 & 0 & \cdots & 1\\
0 & 0 & 0 & \cdots & 0 \end{array}\right) \,.
%\end{displaymath}
\end{equation}
The entries of $e^A$ are given by
$$[e^A]_{ij}=\left\{\begin{array}{ll}
1/(j-i)!,& \textnormal{ if } j\ge i,\\
0, & \textnormal{ else. }
\end{array}\right .
$$
In particular, the diagonal entries of $e^A$ are all equal to 1. 
Therefore, it is impossible to distinguish any of the nodes from
the others on the basis of this centrality measure; yet, it
is clear that the first and last node are rather special, and
certainly more \lq\lq peripheral\rq\rq\ (less \lq\lq central\rq\rq)   
than the other nodes.
Part of the
problem, of course, is that the path digraph contains no closed
walks. In the next section we show one way to extend the notion
of subgraph centrality to digraphs that is immune from such 
shortcomings, and correctly differentiates between nodes in the example
above.
(On the other hand, it is interesting to note that the interpretation
of the off-diagonal
entries of $e^A$ in terms of communicabilities is straightforward
for the directed path. All entries of $e^A$ below the
main diagonal are zero, reflecting the fact that information can only flow
from a node to higher-numbered nodes. Also, the entries of $e^A$ decay
rapidly away from the main diagonal, reflecting 
the fact that the ``ease" of communication 
between a node and a higher numbered one decreases rapidly with the distance.)  Although the diagonal entires of $e^A$ do not provide a meaningful ranking of the nodes in the network, here the row and column sums are able to identify nodes as hubs and authorities.

Another issue when extending the notions of subgraph centrality and
communicability to directed graphs is that
computational difficulties may arise. While the
computations involved do not pose a problem for small 
networks, many real-world networks are large enough that 
directly computing the exponential of the adjacency matrix 
is prohibitive.  In \cite{BenziBoito}, 
techniques for bounding and estimating individual entries of the 
matrix exponential using Gaussian quadrature rules are discussed; see also
\cite{bonchietal} and section \ref{sec:approx} below.  
The ability to find upper and lower bounds for the 
entries requires that the matrix be symmetric, 
thus these bounds cannot be directly computed using the 
adjacency matrix of a directed network.
Again, these difficulties can be
circumvented using the approach discussed in the 
next section.


\section{An extension to digraphs}\label{ExtDig}
Although the techniques described in 
\cite{BenziBoito} cannot be 
directly applied to non-symmetric matrices, setting
\begin{equation}\label{calA}
\mathcal{A} =  \left( \begin{array}{cc}
0 & A \\
A^T & 0  \end{array} \right)
\end{equation}
produces a symmetric matrix $\mathcal{A}$ and, thus, upper 
and lower bounds of individual entries of $e^{\mathcal{A}}$ 
can be computed.  In Proposition \ref{prop1} 
below we relate $e^{\mathcal{A}}$ to the underlying hub 
and authority structure of 
the original digraph. By $B^{\dagger}$ we denote the Moore--Penrose 
generalized inverse of 
matrix $B$.

\begin{proposition}\label{prop1}
Let $\mathcal{A}$ be as described in equation (\ref{calA}).  Then,
$$e^{\mathcal{A}} =  \left( \begin{array}{cc}
\cosh{\left(\sqrt{AA^T}\right)} & A\left(\sqrt{A^TA }\right)^{\dagger}\sinh{\left(\sqrt{A^TA }\right)} \\
\sinh{\left(\sqrt{A^TA}\right)}\left(\sqrt{A^TA }\right)^{\dagger}A^T & \cosh{\left(\sqrt{A^TA}\right)}  \end{array} \right).$$
\end{proposition}

\begin{proof}
Let $A = U\Sigma V^T$ be the SVD of the original (non-symmetric) adjacency matrix $A$.  Then, $\mathcal{A}$ can be decomposed as
 $\mathcal{A} = \left( \begin{array}{cc} U & 0 \\ 0 & V  \end{array} \right) \left( \begin{array}{cc} 0 & \Sigma \\ \Sigma & 0  \end{array} \right) \left( \begin{array}{cc} U^T & 0 \\ 0 & V^T  \end{array} \right)$.  Hence, 
\begin{equation} \label{expA}
e^{\mathcal{A}} = \left( \begin{array}{cc} U & 0 \\ 0 & V  \end{array}\right) \exp{\left( \begin{array}{cc} 0 & \Sigma \\ \Sigma & 0  \end{array} \right)} \left( \begin{array}{cc} U^T & 0 \\ 0 & V^T  \end{array} \right).
\end{equation}
Now,
$$\exp{\left( \begin{array}{cc} 0 & \Sigma \\ \Sigma & 0  \end{array} \right)} = \cosh{\left( \begin{array}{cc} 0 & \Sigma \\ \Sigma & 0  \end{array} \right)} + \sinh{\left( \begin{array}{cc} 0 & \Sigma \\ \Sigma & 0  \end{array} \right)}$$ $$= {\left( \begin{array}{cc} \cosh(\Sigma) & 0 \\ 0 & \cosh(\Sigma)  \end{array} \right)} + {\left( \begin{array}{cc} 0 & \sinh(\Sigma) \\ \sinh(\Sigma) & 0  \end{array} \right)}. $$
Thus,
\begin{equation} \label{expS}
\exp{\left( \begin{array}{cc} 0 & \Sigma \\ \Sigma & 0  \end{array} \right)} = {\left( \begin{array}{cc} \cosh(\Sigma) & \sinh(\Sigma) \\ \sinh(\Sigma) & \cosh(\Sigma) \end{array} \right)}.
 \end{equation}
Putting together equations (\ref{expA}) and (\ref{expS}),
$$e^{\mathcal{A}} =   \left( \begin{array}{cc} U & 0 \\ 0 & V  \end{array}\right)  {\left( \begin{array}{cc} \cosh(\Sigma) & \sinh(\Sigma) \\ \sinh(\Sigma) & \cosh(\Sigma) \end{array} \right)} \left( \begin{array}{cc} U^T & 0 \\ 0 & V^T  \end{array} \right)$$ $$=  \left( \begin{array}{cc}
\cosh{\left(\sqrt{AA^T}\right)} & A\left(\sqrt{A^TA }\right)^{\dagger}\sinh{\left(\sqrt{A^TA }\right)} \\
\sinh{\left(\sqrt{A^TA}\right)}\left(\sqrt{A^TA}\right)^{\dagger}A^T & \cosh{\left(\sqrt{A^TA}\right)}  \end{array} \right).$$
The identities involving the off-diagonal blocks can be easily checked using
the SVD of $A$.
\end{proof}


\subsection{Interpretation of diagonal entries}
\label{sec:interp}
In the context of undirected networks, 
the interpretation of the entries of the matrix exponential 
in terms of subgraph centralities and communicabilities is well-established,
see e.g.~\cite{NetworkProp}.  
In the case of directed networks and $e^{\mathcal{A}}$, things are not as clear.  The network behind $\mathcal{A}$ can be thought of as follows: take the vertices from the original network $A$ and make two copies of them, $V$ and $V'$.  Then, undirected edges exist between the two sets based on the following rule: $E' = \{(i,j') |$there is a directed edge from $i$ to $j$ in the original network$\}$.  This creates a bipartite graph with $2n$ nodes: $1, 2, \ldots ,n,n+1,n+2, \ldots ,2n$. We denote by $V(\mathcal{A})$ this set of
nodes. The use of bipartization to treat rectangular and structurally
unsymmetric matrices is of course standard in numerical linear algebra.

In the undirected case, each node had only one role to play in the network: any information that came into the node could leave by any edge.  In the directed case, there are two roles for each node: that of a hub and that of an authority.  It is unlikely that a high ranking hub will also be a high ranking authority, but each node can still be seen as acting in both of these roles.  In the network $\mathcal{A}$, the two aspects of each node are separated.  Nodes 
$1, 2,\ldots ,n$ in $V(\mathcal{A})$ represent the original nodes 
in their role as hubs and nodes $n+1, n+2,\ldots, 2n$ in $V(\mathcal{A})$ represent the original nodes in their role as authorities.

Given a directed network, an {\em alternating walk} of length $k$, starting with an out-edge, from node $v_1$ to node $v_{k+1}$  is a list of nodes $v_1, v_2, ..., v_{k+1}$ such that there exists edge $(v_i, v_{i+1})$ if $i$ is odd and edge $(v_{i+1}, v_i)$ if $i$ is even:
$$v_1 \rightarrow v_2 \leftarrow v_3 \rightarrow \cdots $$
An {\em alternating walk} of length $k$, starting with an in-edge, from node $v_1$ to node $v_{k+1}$ in a directed network is a list of nodes $v_1, v_2, ..., v_{k+1}$ such that there exists edge $(v_{i+1}, v_i)$ if $i$ is odd and edge $(v_i, v_{i+1})$ if $i$ is even:
$$v_1 \leftarrow v_2 \rightarrow v_3 \leftarrow \cdots $$
From graph theory (see also \cite{croftetal}), it is known that $[AA^TA \ldots]_{ij}$ (where there are $k$ matrices being multiplied) counts the number of alternating walks of length $k$, starting with an out-edge, from node $i$ to node $j$, 
whereas  $[A^TAA^T\ldots]_{ij}$ (where there are $k$ matrices being multiplied) counts the number of alternating walks of length $k$, starting with an in-edge, from node $i$ to node $j$.  That is, $[(AA^T)^k]_{ij}$ and $[(A^TA)^k]_{ij}$ count the number of alternating walks of length $2k$.

In the original network $A$, if node $i$ is a good hub, it will point to many good authorities, which will in turn be pointed at by many hubs.  These hubs will also point to many authorities, which will again be pointed at by many other hubs.  Thus, if $i$ is a good hub, it will show up many times in the sets of hubs described above.  That is, there should be many even length alternating walks, starting with an out-edge, from node $i$ to itself.  Giving a walk of length $2k$ a weight of 
$\frac{1}{(2k)!}$, these walks can be counted using the $(i,i)$ entry of the matrix
$$I + \frac{AA^T}{2!} + \frac{AA^TAA^T}{4!} + \cdots + \frac{(AA^T)^{k}}{(2k)!}+ \cdots$$
Letting $A = U\Sigma V^T$ be the SVD of $A$, this becomes:
%$$U\left(I + \frac{AA^T}{2!} + \frac{AA^TAA^T}{4!} + \cdots + 
%\frac{(AA^T)^{k}}{(2k)!}+ \cdots\right)U^T$$ $$= U\cosh(\Sigma)U^T 
%= \cosh(\sqrt{AA^T})$$
$$U\left(I + \frac{\Sigma ^2}{2!} + \frac{\Sigma ^4}{4!} + \cdots +
\frac{\Sigma^{2k}}{(2k)!}+ \cdots\right)U^T$$ $$= U\cosh(\Sigma)U^T
= \cosh(\sqrt{AA^T})\,.$$

The {\em hub centrality} of node $i$ (in the original network) 
is thus given by $$[e^{\mathcal{A}}]_{ii} = [\cosh(\sqrt{AA^T})]_{ii}.$$  
This measures how well node $i$ transmits information to the 
authoritative nodes in the network.

Similarly, if node $i$ is a good authority, there will be many even length alternating walks, starting with an in-edge, from node $i$ to itself.  Giving a walk of length $2k$ a weight of $\frac{1}{(2k)!}$, these walks can be counted using the $(i,i)$ entry of $\cosh(\sqrt{A^TA})$.  

Hence, the {\em authority centrality} of node $i$ is given by 
$$[e^{\mathcal{A}}]_{n+i,n+i}= [\cosh(\sqrt{A^TA})]_{ii}.$$  
It measures how well node $i$ receives information from the hubs in the network.

Note that the traces of the two diagonal blocks in $e^{\cal A}$ are identical,
so each accounts for half of the Estrada index of the bipartite graph.
Also, recalling the well-known fact that the eigenvalues of $\cal A$
are $\pm \sigma_i$ where $\sigma_i$ denotes the singular values of $A$,
we have 
$${\rm Tr}\,(e^{\cal A}) = 
\sum_{i=1}^n e^{\sigma_i} + \sum_{i=1}^n e^{-\sigma_i} 
 = 2\sum_{i=1}^n \cosh\,(\sigma_i),$$
an identity that can also be obtained directly from the
expression for $e^{\cal A}$ given in Proposition \ref{prop1}.

Returning to the example of the directed path graph with adjacency
matrix $A$ given by (\ref{AP}), one finds that using the diagonal entries of
$e^{\cal A}$ to rank the nodes gives node 1 as the least authoritative
node, and node $n$ as the one with the lowest hub ranking, with all
the other nodes being tied. Thus we see that, while $e^A$ fails to 
differentiate between the nodes of this graph, using $e^{\cal A}$
yields a very reasonable hub/authority ranking of the nodes.

\subsection{Interpretation of off-diagonal entries}
Although not used in the remainder of this paper, for the sake of
completeness we give here
an interpretation of the off-diagonal entries
of $e^{\cal A}$. As we will see, this interpretation is
rather different from the one usually given for the off-diagonal
entries of $e^A$, and provides information of a different nature
on the structure of the underlying graph. 

In discussing the off-diagonal entries of $\mathcal{A}$, there are three 
blocks to consider.  First, there are the off-diagonal entries of the upper-left block, $\cosh(\sqrt{AA^T})$, then there are the off-diagonal entries of the lower-right block, $\cosh(\sqrt{A^TA})$.  Finally, there is the off-diagonal block, 
$ A\left(\sqrt{A^TA }\right)^{\dagger}\sinh \left(\sqrt{A^TA }\right)$
(the fourth block in $e^{\cal A}$ being its transpose).

From section \ref{sec:interp}, 
$[e^{\mathcal{A}}]_{ij} = [\cosh(\sqrt{AA^T})]_{ij}$, $1 \leq i,j \leq n$, counts the number of even length alternating walks, starting with an out-edge, from node $i$ to node $j$, weighting walks of length $2k$ by a factor of $\frac{1}{(2k)!}$.  When $i \neq j$, these entries measure how similar nodes $i$ and $j$ are as hubs.  That is, if nodes $i$ and $j$ point to many of the same nodes, there will be many short even length alternating walks between them.

The {\em hub communicability} between nodes $i$ and $j$, $1 \leq i,j \leq n$, is given by 
$$[e^{\mathcal{A}}]_{ij} = [\cosh(\sqrt{AA^T})]_{ij}$$
This measures how similar nodes $i$ and $j$ are in their roles as hubs.  That is, a large value of hub communicability between nodes $i$ and $j$ indicates that they point to many of the same authorities.  In other words, they point to nodes which are authorities on the same subjects.

Similarly, $[e^{\mathcal{A}}]_{n+i,n+j} = [\cosh(\sqrt{A^TA})]_{ij}$, $1 \leq i,j \leq n$, counts the number of even length alternating walks, starting with an in-edge, from node $i$ to node $j$, also weighing walks of length $2k$ by a factor of $\frac{1}{(2k)!}$.  When $i \neq j$, these entries measure how similar the two nodes are as authorities.  If $i$ and $j$ are pointed at by many of the same hubs, there will be many short even length alternating walks between them.

The {\em authority communicability} between nodes $i$ and $j$, $1 \leq i,j, \leq n$, is given by 
$$[e^{\mathcal{A}}]_{i+n,j+n} = [\cosh(\sqrt{A^TA})]_{ij}$$
This measures how similar nodes $i$ and $j$ are in their roles as authorities.  That is, a large value of authority communicability between nodes $i$ and $j$ means that they are pointed to by many of the same hubs and, as such, are likely to contain information on the same subjects.

Let us now consider the off-diagonal blocks of $\mathcal{A}$. 
Here, $[\sinh(\sqrt{A^TA})]_{ij}$ counts the number of odd length alternating walks, starting with an out-edge, from node $i$ to node $j$, weighing walks of length 
$2k+1$ by $\frac{1}{(2k+1)!}$.  This measures the communicability between node $i$ as a hub and node $j$ as an authority.

The {\em hub-authority communicability} between nodes $i$ and $j$ (that is, the communicability between node $i$ as a hub and node $j$ as an authority) is given by:
$$[e^{\mathcal{A}}]_{i,n+j} = [A\left(\sqrt{A^TA }\right)^{\dagger}\sinh \left(\sqrt{A^TA }\right)]_{ij} $$ $$= [\sinh{\left(\sqrt{A^TA}\right)}\left(\sqrt{A^TA}\right)^{\dagger}A^T]_{ji} = [e^{\mathcal{A}}]_{n+j, i}.$$
A large hub-authority communicability between nodes $i$ and $j$ means that they are likely in the same \lq\lq part\rq\rq\ of the directed network: node $i$ tends to point to nodes that contain information similar to that on which node $j$ is an authority.

\subsection{Relationship with HITS}\label{hits_comp}
As described in \ref{sec:newHITS}, the HITS ranking of nodes 
as hubs and authorities uses only information from the dominant 
eigenvector of $\mathcal{A}$. Here we show that when using the
diagonal of $e^{\mathcal A}$, we exploit information contained
in all the eigenvectors of $\mathcal{A}$; moreover, the HITS rankings
can be regarded as an approximation of those given by the
diagonal entries of $e^{\mathcal A}$.

Assume the eigenvalues of $\mathcal{A}$ can be ordered as
$\lambda_{1} > \lambda_{2} \geq \lambda_{3} \geq \cdots \geq 
\lambda_{2n}$.
Then, $\mathcal{A}$ can be written as $\mathcal{A} = \sum_{i=1}^{2n} 
\lambda_i{\bf u}_i{\bf u}_i^T$ where ${\bf u}_1, {\bf u}_2,\ldots, {\bf u}_{2n}$ are the 
normalized eigenvectors of $\mathcal{A}$.  Taking the exponential of 
$\mathcal{A}$, we get: 
$$e^{\mathcal{A}} = \sum_{i=1}^{2n} 
e^{\lambda_i}{\bf u}_i{\bf u}_i^T = e^{\lambda_{1}}{\bf u}_1{\bf u}_1^T + \sum_{i=2}^{2n} 
e^{\lambda_i}{\bf u}_i{\bf u}_i^T.$$

Now, the hub and authority rankings come from the diagonal entries 
of $e^{\mathcal{A}}$: 
$${\rm diag}\,(e^{\mathcal{A}}) = e^{\lambda_{1}}
{\rm diag}\,({\bf u}_1{\bf u}_1^T) + \sum_{i=2}^{2n}e^{\lambda_i}{\rm diag}\,({\bf u}_i{\bf u}_i^T).$$ 
Rescaling the hub and authority scores by $e^{\lambda_{1}}$ does not alter
the rankings; hence, we can instead consider
$${\rm diag}\,(e^{-\lambda_{1}}e^{\mathcal{A}}) = 
{\rm diag}\,(e^{\mathcal{A}-\lambda_{1}I}) =
{\rm diag}\,({\bf u}_1{\bf u}_1^T) + \sum_{i=2}^{2n}e^{\lambda_i-\lambda_{1}}{\rm diag}\,({\bf u}_i{\bf u}_i^T).$$
Now, the diagonal entries of the rank-one matrix ${\bf u}_1{\bf u}_1^T$ are 
just the squares of the (nonnegative) entries of the dominant eigenvector of
$\mathcal{A}$; hence,
the rankings provided by the first
term in the expansion of $e^{\mathcal{A}}$ in the eigenbasis of $\cal A$ are 
precisely those given by HITS.

It is also clear that if $\lambda_{1} \gg \lambda_{2}$, 
then the rankings provided
by the diagonal entries of $e^{\mathcal{A}}$ are unlikely to differ much from
those of HITS, since the weights $e^{\lambda_i-\lambda_{1}}$ will be tiny, for
all $i=2,\ldots ,2n$. 
 Conversely, if the gap between $\lambda_{1}$ and the rest of the
spectrum is small ($\lambda_{1}\approx \lambda_{2}$), then the contribution 
from the remaining eigenvectors,
$\sum_{i=2}^{2n}e^{\lambda_i-\lambda_{1}}{\rm diag}\,({\bf u}_i{\bf u}_i^T)$, may be 
non-negligible relative to the first term and therefore the resulting rankings
could differ significantly from those obtained using HITS. %This is true, 
%{\em a fortiori}, if the dominant eigenvalue of $\cal A$ is not simple. 
In section \ref{sec:applications} we will see examples of real networks illustrating
both scenarios.

Summarizing, use of the matrix exponential for ranking hubs and authorities
amounts to using the (squared) entries of {\em all} the eigenvectors of $\cal A$,
weighted by the exponential of the corresponding eigenvalues. 
Of course, in place of the exponential, a number of other functions
could be used; see the discussion in the next section. 
As shown above, the HITS ranking scheme 
uses the leading term only, corresponding to the approximation 
$e^{\mathcal A}\approx e^{\lambda_1}{\bf u}_1{\bf u}_1^T$. Between these two extremes
one could also use approximations of the form
\begin{equation}\label{k_terms}
e^{\mathcal A} \approx \sum_{i=1}^k e^{\lambda_i}{\bf u}_i{\bf u}_i^T\,,
\end{equation}
where $1< k < n$; indeed, in most cases of practical interest a modest
value of $k$ ($\ll n$)
 will be sufficient for a very good approximation, since the
eigenvalues of $\cal A$ are often observed to decay rapidly from a certain
index $k$ onward. We return on this topic in section \ref{sec:approx}.


\section{Other ranking schemes}
In this section we discuss a few other schemes that have been proposed
in the literature, and compare them with the hub and authority centrality
measures based on the exponential of $\cal A$.

\subsection{Resolvent-based measures}
As already noted, another function that has been successfully
used to define centrality and communicability measures for an undirected
network is the matrix
resolvent, $(I-\alpha A)^{-1}$. This approach was pioneered early
on by Katz \cite{Katz}, and variants thereof have since been used 
by numerous authors; see, e.g., 
\cite{bonchietal,Brandes,EHB11,NetworkProp,franceschet,vigna}. It is
hardly necessary to mention the close relationship existing between the
resolvent and the exponential function, which can be expressed via the
Laplace transform. For the adjacency matrix $\cal A$ of a bipartite
graph given by (\ref{calA}), the resolvent
is easily determined to be
\begin{equation}\label{resolv}
(I - \alpha {\cal A})^{-1} = 
\left( \begin{array}{cc} (I - \alpha^2AA^T)^{-1} & \alpha A(I - \alpha^2A^TA)^{-1}\\
                         \alpha(I - \alpha^2A^TA)^{-1}A^T     & (I - \alpha^2A^TA)^{-1}
                         \end{array} \right).
\end{equation}
The condition on $\alpha$ can be expressed as $0 < \alpha < 1/\sigma_1$, 
where $\sigma_1=\|A\|_2$
denotes the largest singular value of $A$, the adjacency matrix of the undirected 
network. This ensures that the matrix in (\ref{resolv}) is well-defined and
 nonnegative, with
positive diagonal entries. The diagonal entries of $(I - \alpha^2AA^T)^{-1}$
provide the hub scores, those of $(I - \alpha^2A^TA)^{-1}$ the authority scores. 
A drawback of this approach is the need to select
the parameter $\alpha$, and the fact that different values of $\alpha$ may lead
to different rankings. We have performed numerical experiments with this
approach and we found that for certain values of $\alpha$, particularly 
those close to the upper limit $1/\sigma_1$, the hub and authority
rankings obtained with the resolvent function are not too different
from those obtained with the matrix exponential, which is not surprising in light of the results from Chapter \ref{C:Robustness}. Moreover,
not surprisingly, as the value of $\alpha$ is reduced, one obtains hub and
authority rankings that
are strongly correlated with the out- and in-degree of the nodes,
respectively.
 Overall, because the resolvent tends to weigh short walks more
heavily than the exponential, and since longer walks contribute relatively
little to the centrality scores, it is fair to say that the exponential
is less \lq\lq degree biased\rq\rq\ than the resolvent function. Also,
since (for $\beta = 1$) the exponential rankings do not depend on a tunable parameter,
they provide unambiguous rankings.


We note that
``Katz'' authority and hub scores may also be obtained by considering 
the column and row sums of the (nonsymmetric) matrix resolvent
$(I - \alpha A)^{-1}$, where $A$ is the adjacency
matrix of the original digraph and $\alpha>0$ is again assumed to
be small enough for the corresponding Neumann series to converge. 
Indeed, the row sums of $(I - \alpha A)^{-1}$ count the number of
(weighted) walks out of each node, while the column sums count the
number of (weighted) walks into each node. Denoting by $\bf 1$ the
vector of all ones, hub and authority rankings can be obtained by solving the
two linear systems
\begin{equation}\label{katz}
(I - \alpha A){\bf y} = {\bf 1} \quad {\rm and} \quad (I - \alpha A^T){\bf x} = {\bf 1}\,,
\end{equation}
respectively. Here the parameter $\alpha$ must satisfy $0 < \alpha < 1/\rho(A)$,
where $\rho(A)$ denotes the spectral radius of $A$.
The results of numerical experiments comparing the Katz scores with those
based on the exponential of $\cal A$ are given in section \ref{sec:applications}.
Here we observe that these Katz scores are also dependent on the choice
of the parameter $\alpha$, and similar considerations to those made for
$(I - \alpha {\cal A})^{-1}$ apply.

A natural analogue to this approach is the use 
of row and column sums of the
exponential $e^A$ to rank hubs and authorities. Some results obtained with 
this approach are discussed in section \ref{sec:applications}. We note that this method 
is different from the {\em Exponentiated Inputs HITS Method} of \cite{EIHITS}.  This method is a modification to HITS, which was developed in order to correct the issue of non unique results in certain networks.  If the dominant eigenvalue of $A^TA$ (and, consequently of $AA^T$ ) is not simple, then the corresponding eigenspace is multidimensional. This means that the choice of the initial vector affects the convergence of the HITS algorithm and different hub and authority vectors can be produced using different initial vectors. This can occur when $A^{T} A$ is reducible, that is when the original network is not strongly connected.    In \cite{EIHITS}, Farahat {\em et. al.}  propose a modification to the HITS algorithm which amounts to replacing $A$ and $A^T$ with $e^A-I$ and $(e^A-I)^T$ in the HITS iteration.  They prove that, as long as the original network is weakly connected, the dominant eigenvalue of $(e^A-I)^T (e^A - I)$ is simple. Thus, HITS with this exponentiated input produces unique hub and authority rankings.  However, a result of this replacement is that nodes with 0 in-degree (or a low in-degree) are less important in the calculation of authority scores than nodes with a high in-degree.   When there are many nodes with 0 in-degree or whose edges point to only a few other nodes, dropping these edges can greatly affect the HITS rankings.

An obvious disadvantage of this algorithm is its cost, since it requires iterating with a matrix exponential and its transpose. It can be implemented using only matrix- vector products involving $\tilde{A}$ and $\tilde{A}^T$ by means of techniques, like Krylov subspace methods, for evaluating the action of a matrix function on a given vector; see, e.g., \cite[Chapter 13]{highambook}. This approach leads to a nested iteration scheme, with HITS as the outer iteration and the Krylov method as the inner one.

\subsection{PageRank and Reverse PageRank}
As mentioned in Section \ref{sec:PageRank}, the (now) classical PageRank algorithm 
provides a means of finding
the authoritative nodes in a
digraph.  Details about the PageRank algorithm can be found in Chapter \ref{C:introduction}, Section \ref{sec:iterative_methods}. 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$. 
 Like HITS, PageRank and Reverse PageRank are
eigenvector-based ranking algorithms that do not take into account
information about the network contained in the non-dominant eigenvectors.
As already mentioned,
it has been argued \cite{MiPapa} that eigenvector-based algorithms
tend to be degree-biased. 
Furthermore, like the Katz-type algorithms, the rankings obtained
depend on the choice of a tuneable damping parameter.
While the success of PageRank in finding authoritative nodes is
well known and very well documented, much less is known about the
effectiveness of Reverse PageRank in identifying hubs; some references
are \cite{Bar-Yossef,crofootetal,WuLi,WuLietal}. We present
the results of a few numerical experiments with PageRank and Reverse PageRank
in section \ref{sec:applications}.


\section{Examples on small digraphs}\label{expes}
In this section and the next we illustrate the proposed method on some simple networks
of small size, as well as on some larger data sets corresponding
to real networks. We also compare our approach with HITS and other
rankings schemes, including Katz, PageRank and Reverse PageRank. Here, we compare out and in-degree counts, HITS, and
our proposed method to obtain hub and authority rankings in a 
few small digraphs. The purpose of this section is mostly
pedagogical.

\subsection{Example 1}
Consider the small directed network in Fig.~\ref{graphs} (left panel).
The adjacency matrix is given by
$$A =  \left( \begin{array}{cccc} 0 & 1 & 1 & 0 \\ 1 & 0 & 1 &0 \\ 0 & 1 & 0 & 1 \\ 0 & 1 & 0 & 0 \end{array} \right).$$

\begin{figure}[t!]
\centering
\includegraphics[width=0.30\textwidth]{digraph.eps}
\hspace{0.3in}
\includegraphics[width=0.22\textwidth]{newbipartite1.eps}
\caption{The original directed network from Example 1, with adjacency matrix $A$ (left) 
and the bipartite network with adjacency matrix $\mathcal{A}$ (right).}
\label{graphs}
\end{figure}

%\begin{figure}[ht]
%\centering
%\subfigure[The original directed network, with adjacency matrix $A$.]{
%%\includegraphics[width=0.25\textwidth]{digraph.eps}
%\includegraphics[width=0.25\textwidth]{digraph}
%\label{fig:subfig1}
%}
%%\hspace{.3in}
%%\subfigure[The bipartite network with adjacency matrix $\mathcal{A}$]{
%\includegraphics[width=0.18\textwidth]{bipartite.eps}
%\includegraphics[width=0.18\textwidth]{bipartite}
%\label{fig:subfig2}
%}
%\label{fig:subfigureExample}
%\caption[Optional caption for list of figures]{Networks behind $A$ and $\mathcal{A}$.}
%\end{figure}

The corresponding bipartite graph is shown in  Fig.~\ref{graphs} (right panel).
If hubs and authorities are determined simply using in-degree and out-degree 
counts, the result is as follows: 

\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & out-degree & in-degree \\ \hline
1 & 2 & 1\\
2 & 2 & 3\\
3 & 2 & 2\\
4 & 1 & 1\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

Under this ranking, the hub ranking of the nodes is: 
$\{1, 2, 3 {\rm \,(tie)}; 4\}$.  The authority ranking 
of the nodes is: $\{2;3; 1,4 {\rm \, (tie)}\}$. %These rankings
%are not satisfactory as they result in many nodes being tied.
We obtain somewhat different results
using the HITS algorithm. The 
eigenvectors of $AA^T$ and $A^TA$ corresponding to the 
largest eigenvalue $\lambda_{\max} \approx 3.9563$, which is
simple, yield the
following rankings for hubs and authorities:

\vspace{0.12in}


\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & hub rank & authority rank \\ \hline
1 & .3383 & .0965 \\
2 & .1729 & .4618\\
3 & .2798 & .2854\\
4 & .2091 & .1562\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

Here, the ranking for hubs is: $\{1;3;4;2\}$. 
The ranking for authorities is: $\{2;3;4;1\}$.
Note that node 2, which was given a top
hub score by looking just at the out-degrees, is judged
by HITS as the node with the lowest hub score. 

Using $e^{\mathcal{A}}$ as described above, the 
rankings for hub centralities and authority centralities are:

\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & hub centrality = $[e^{\mathcal{A}}]_{ii}$ & authority centrality 
= $[e^{\mathcal{A}}]_{4+i,4+i}$  \\ \hline
1 & 2.3319 & 1.5906\\
2 & 2.2289 & 3.0209\\
3 & 2.2812 & 2.2796\\
4 & 1.6414 & 1.5922\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}


With this method, the hub ranking of the nodes is: $\{1;3;2;4\}$.  
The authority ranking is: $\{2;3;4;1\}$. On this example, our method produces  
the same authority ranking as HITS. The hub ranking, however, is slightly
different: both methods identify node 1 as the one with the highest hub
score, followed by node 3; however, our method assigns the lowest hub score
to node 4 rather than node 2. This is arguably a more meaningful
ranking. 

\subsection{Example 2}
Consider the small directed network in Fig.~\ref{graphs2} (left panel). The adjacency matrix is given by 
$$A =  \left( \begin{array}{cccc} 0 & 0 & 1 & 0 \\ 1 & 0 & 0 &1 \\ 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 \end{array} \right).$$

\begin{figure}[t!]
\centering
\includegraphics[width=0.30\textwidth]{digraph2.eps}
\hspace{0.3in}
\includegraphics[width=0.21\textwidth]{newbipartite2.eps}
\caption{The original directed network from Example 2, with adjacency matrix $A$ (left) 
and the bipartite network with adjacency matrix $\mathcal{A}$ (right).}
\label{graphs2}
\end{figure}

The corresponding bipartite graph is shown in Fig.~\ref{graphs2} (right panel).  If hubs and authorities are determined only using in-degrees and out-degrees, the result is as follows:

\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & out-degree & in-degree \\ \hline
1 & 1 & 1\\
2 & 2 & 2\\
3 & 1 & 1\\
4 & 1 & 1\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

Under this criterion, the hub and authority rankings are
both %given by
$\{2; 1, 3, 4\, {\rm (tie)}\}$.
While it is intuitive that 
node 2 should be given a high score (both as an authority
and as a hub), just looking at the degrees does not allow
one to distinguish the remaining nodes.

Consider now the use of HITS.  The largest eigenvalue 
of $AA^T$ (and $A^TA$) is $\lambda_{\max} =2$ and it has multiplicity two.  
Thus, different starting vectors for the HITS algorithm may produce different rankings,
as discussed in \cite{EIHITS}.  
Starting from a constant authority vector $x^{(0)}$,
as suggested in \cite{HITS},
produces the following scores:

\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & hub rank & authority rank \\ \hline
1 & .0000 & .3333\\
2 & .5000 & .3333\\
3 & .2500 & .0000\\
4 & .2500 & .3333\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

\noindent The ranking for hubs is: $\{2; 3, 4\, {\rm (tie)}; 1\}$.  The ranking for 
authorities is the following: 
$\{1, 2, 4\, {\rm (tie)}; 3\}$. 
%The fact that node 2 has a higher in-degree than any 
%of the other nodes does not 
%affect its authority ranking under the HITS algorithm, which shows a clear
%drawback of HITS.

If the ranking is determined using $e^{\mathcal{A}}$ as described above, 
the resulting scores are:

\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & hub centrality = $[e^{\mathcal{A}}]_{ii}$ & authority 
centrality = $[e^{\mathcal{A}}]_{4+i,4+i}$  \\ \hline
1 & 1.5431 & 1.5891\\
2 & 2.1782 & 2.1782\\
3 & 1.5891 & 1.5431\\
4 & 1.5891 & 1.5891\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

\noindent With this method, the hub ranking of the nodes is the same as in HITS: 
$\{2; 3,4\, {\rm (tie)}; 1\}$.  However, in the authority ranking, 
node 2 is the clear winner rather than being part of a three-way 
tie for first place: $\{2; 1, 4\, {\rm (tie)}; 3\}$. %Identical rankings
%are obtained by the exponentiated input HITS algorithm. 
In this example, the method based on the matrix exponential is
able to identify a top authority node by making
use of additional spectral information. 

\subsection{Example 3}
\label{example3}
Let $G$ be the small directed network in Fig.~\ref{graphs3}.  The adjacency matrix is given by
$$A =  \left( \begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 1 & 0 \end{array} \right).$$

\begin{figure}[t!]
\centering
\includegraphics[width=0.35\textwidth]{example3.eps}
\hspace{0.3in}
\includegraphics[width=0.2\textwidth]{example3b.eps}
\caption{The original directed network from Example 3, with adjacency matrix $A$ (left)
and the bipartite network with adjacency matrix $\mathcal{A}$ (right).}
\label{graphs3}
\end{figure}

If hubs and authorities are determined using only in-degrees and out-degrees, the result is:

\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & out-degree & in-degree \\ \hline
1 & 0 & 4\\
2 & 1 & 1\\
3 & 1 & 1\\
4 & 1 & 1\\
5 & 1 & 1\\
6 & 4 & 0\\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

The hub ranking of the nodes using degrees is: \{6; 2,3,4,5 (tie); 1\}.  The authority ranking is \{1; 2,3,4,5 (tie); 6\}.

If the HITS algorithm is used, the resulting rankings are similar, 
but not exactly the same.  Starting with a constant authority vector 
%$x^{(0)} = [\frac{1}{6}, \frac{1}{6}, \frac{1}{6}, \frac{1}{6}, 
%\frac{1}{6}, \frac{1}{6}]$, the rankings produced are:
$x^{(0)}$, the results are:


\vspace{0.12in}

\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & hub rank & authority rank \\ \hline
1 & .000 & .200 \\
2 & .125 & .200 \\
3 & .125 & .200 \\
4 & .125 & .200 \\
5 & .125 & .200 \\
6 & .500 & .000 \\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

The hub ranking of the nodes is: $\{6; 2,3,4,5\, {\rm (tie)}; 1\}$.  
The authority ranking is: \{1,2,3,4,5 (tie); 6\}.  Here, HITS does
not differentiate
between node 1 and nodes 2, 3, 4, and 5 in 
terms of the authority score, even though node 1 has by far the highest 
in-degree.
This appears as a failure of HITS, since it is intuitive that node 1
should be regarded as very authoritative.
%If HITS with exponentiated input is used, node 1 does get a 
%higher authority score than all the other nodes:

%\vspace{0.12in}

%\begin{center}
%\begin{tabular}{|c|c|c|}
%\hline
%node & Exp.~input HITS hub rank & Exp.~input HITS authority rank \\ \hline
%1 & .0000 & .4472 \\
%2 & .1382 & .1382 \\
%3 & .1382 & .1382 \\
%4 & .1382 & .1382 \\
%5 & .1382 & .1382 \\
%6 & .4472 & .0000 \\
%\hline
%\end{tabular}
%\end{center}

%\vspace{0.12in}

%Under HITS with exponentiated input, the hub ranking is the
%same one obtained by HITS, namely, \{6; 2,3,4,5 (tie); 1\},  
%while the authority ranking is: \{1; 2,3,4,5 (tie); 6\}.

When $e^{\mathcal{A}}$ is used to calculate the hub and authority scores, node 1 does get a higher authority ranking than all the other nodes:

\vspace{0.12in}


\begin{center}
\begin{tabular}{|c|c|c|}
\hline
node & hub centrality = $[e^{\mathcal{A}}]_{ii}$ & authority 
centrality = $[e^{\mathcal{A}}]_{6+i,6+i}$  \\ \hline
1 & 1.0000 & 3.7622 \\
2 & 1.6905 & 1.6905 \\
3 & 1.6905 & 1.6905 \\
4 & 1.6905 & 1.6905 \\
5 & 1.6905 & 1.6905 \\
6 & 3.7622 & 1.0000 \\
\hline
\end{tabular}
\end{center}

\vspace{0.12in}

Note that, if desired, the value 1 can be subtracted from these 
scores since it does not affect the relative ranking of the nodes.
The hub ranking is \{6; 2,3,4,5 (tie); 1\}, and 
the authority ranking is: \{1; 2,3,4,5 (tie); 6\}.  

%One may ask under which conditions, if any, the exponentiated
%input HITS method and the one based on $e^{\cal A}$ result in
%different rankings.
%As we will see in the next example, when there 
%are many nodes with 0 in-degree, the rankings using  
%$e^{\mathcal{A}}$ start to differ from those using HITS 
%with exponentiated input.

%\subsection{Example 4}
%\label{example4}
%Let $G$ be the 10-node directed network in Fig.~\ref{graphs4}.  The adjacency matrix of $G$ is given by:
%$$A =  \left( \begin{array}{cccccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0
%& 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\  0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \end{array} \right).$$


%\begin{figure}[t!]
%\centering
%\includegraphics[width=0.6\textwidth]{newexample4.eps}
%\hspace{0.3in}
%\includegraphics[width=0.14\textwidth]{example4b.eps}
%\caption{The original directed network from Example 4, with adjacency matrix $A$ (left)
%and the bipartite network with adjacency matrix $\mathcal{A}$ (right).}
%\label{graphs4}
%\end{figure}

%Using just in-degree and out-degree, the hub and authority 
%rankings of $G$ are as follows:

%\vspace{0.12in}

%\begin{center}
%\begin{tabular}{|c|c|c|}
%\hline
%node & out-degree & in-degree \\ \hline
%1 & 0 & 4\\
%2 & 1 & 1\\
%3 & 1 & 1\\
%4 & 1 & 1\\
%5 & 1 & 5\\
%6 & 4 & 0\\
%7 & 1 & 0\\
%8 & 1 & 0\\
%9 & 1 & 0\\
%10 & 1 & 0\\
%\hline
%\end{tabular}
%\end{center}

%\vspace{0.12in}


%The hub ranking is: \{6; 2,3,4,5,7,8,9,10 (tie); 1\}.  
%The authority ranking is \{5; 1; 2,3,4 (tie); 6,7,8,9,10 (tie)\}.
%Using HITS, different rankings are produced.  
%Starting with authority vector $x_0 = (\frac{1}{10},
%\frac{1}{10}, \frac{1}{10}, \frac{1}{10}, \frac{1}{10}, 
%\frac{1}{10}, \frac{1}{10}, \frac{1}{10}, \frac{1}{10}, \frac{1}{10})$, 
%the rankings produced are:
%Starting with a constant authority vector $x^{(0)}$,
%the rankings produced are:

%\vspace{0.12in}

%\begin{center}
%\begin{tabular}{|c|c|c|}
%\hline
%node & hub rank & authority rank \\ \hline
%1 & .0000 & .0008\\
%2 & .0003 & .1665\\
%3 & .0003 & .1665\\
%4 & .0003 & .1665\\
%5 & .0003 & .4996\\
%6 & .3330 & .0000\\
%7 & .1665 & .0000\\
%8 & .1665 & .0000\\
%9 & .1665 & .0000\\
%10 & .1665 & .0000\\
%\hline
%\end{tabular}
%\end{center}

%\vspace{0.12in}

%The hub ranking is: \{6; 7,8,9,10 (tie); 2,3,4,5 (tie); 1\}.  The authority ranking is: \{5; 2,3,4
%(tie); 1; 6,7,8,9,10 (tie)\}.  Since nodes 2, 3, 4, and 5 are given such low hub scores, node 1 is
%deemed an extremely low authority, even though it has an in-degree of 4.

%By comparison, when using HITS with exponentiated input, node 1 is considered very authoritative:

%\vspace{0.12in}
%\begin{center}
%\begin{tabular}{|c|c|c|}
%\hline
%node & Exp.~input HITS hub rank & Exp.~input HITS authority rank \\ \hline
%1 & .0000 & .4196\\
%2 & .0849 & .1095\\
%3 & .0849 & .1095\\
%4 & .0849 & .1095\\
%5 & .0849 & .2519\\
%6 & .2871 & .0000\\
%7 & .0934 & .0000\\
%8 & .0934 & .0000\\
%9 & .0934 & .0000\\
%10 & .0934 & .0000\\
%\hline
%\end{tabular}
%\end{center}

%\vspace{0.12in}


%The hub ranking is: \{6; 7,8,9,10 (tie); 2,3,4,5 (tie); 1\}, identical to
%HITS.  
%The authority ranking is: \{1; 5; 2,3,4 (tie); 6,7,8,9,10 (tie)\}.  
%Here, even though node 5 has a higher in-degree than node 1, node 1
%is given a higher authority rank.

%On the other hand,
%using the rankings from $e^{\mathcal{A}}$, node 5 is considered a 
%higher authority than node 1:

%\vspace{0.12in}

%\begin{center}
%\begin{tabular}{|c|c|c|}
%\hline
%node & hub centrality = $[e^{\mathcal{A}}]_{ii}$ & 
%authority centrality = $[e^{\mathcal{A}}]_{n+i,n+i}$  \\ \hline
%1 & 1.0000 & 3.7622\\
%2 & 1.6905 & 1.6974\\
%3 & 1.6905 & 1.6974\\
%4 & 1.6905 & 1.6974\\
%5 & 1.6905 & 4.9203\\
%6 & 4.0063 & 1.0000\\
%7 & 1.7516 & 1.0000\\
%8 & 1.7516 & 1.0000\\
%9 & 1.7516 & 1.0000\\
%10 & 1.7516 & 1.0000\\
%\hline
%\end{tabular}
%\end{center}

%\vspace{0.12in}

%The hub ranking is the same as with HITS and its
%exponentiated input variant. % \{6; 7,8,9,10 (tie); 2,3,4,5 (tie); 1\}.  
%The authority ranking is \{5; 1; 2,3,4 (tie); 6,7,8,9,10 (tie)\}.
%Nodes 1 and 5 swap place in the authority ranking, with everything
%else staying the same. Whether this ranking is \lq\lq better\rq\rq\
%than the one provided by HITS with exponentiated input is open to debate. %Node 1
%is pointed to by less nodes than 5, but one of the nodes pointing
%toit, node 5, is quite authoritative. Node 5 has more nodes pointing
%to it, but none of them is authoritative. 
%What is clear is that both methods 
%identify nodes 1 and 5 as the most
%authoritative ones by a considerable margin, whereas HITS with a uniform
%starting vector completely
%fails in identifying node 1 as authoritative.



\section{Application to web graphs}\label{sec:applications}
Similarly to HITS, and in analogy to subgraph centrality for
undirected networks,  the rankings produced by the values 
on the diagonal of $e^{\mathcal{A}}$ can be used to 
rank websites as hubs and authorities in web searches 
(many other applications are of course also possible).
 Three of the data sets considered here are small web graphs consisting of
web sites on various topics and can be found at \cite{datasets} 
along with the website associated with each node;
see also \cite{LinkStructure}.  The experiments for this paper were run on the ``Expanded" version of the data sets.
Each data set is named after the corresponding
topic.\footnote{It should be noted, 
however, that in the node list for the adjacency matrix, 
the node labeling begins with 1 and in the list of websites 
associated with the nodes found at \cite{datasets}, node 
labeling begins at 0.  Thus, node $i$ in the adjacency 
matrix is associated with website $i-1$.} 
In addition, we include results for the {\em wb-cs-stanford}
data set from the University of Florida Sparse Matrix Collection
\cite{UFL}. This digraph represents a subset of the
 Stanford University web. In this section, 
the hub and authority rankings obtained from 
$e^{\mathcal{A}}$ are compared with those from HITS, 
%implemented using the Matlab code described in \cite{Pagerank}.  
Katz (using (\ref{katz}) with $\alpha=1/(\rho(A) + 0.1)$),
the row and column sums of the exponential $e^A$ of
the nonsymmetric matrix $A$, and PageRank/Reverse PageRank. 
For the latter we use the standard value $\alpha = 0.85$ for the damping
parameter.
All experiments are performed using Matlab Version 7.9.0 
(R2009b) on a MacBook Pro running 
OS X Version 10.6.8, a 2.4 GHZ Intel Core i5 processor 
and 4 GB of RAM. For the purpose of these tests we use the
built-in Matlab function {\tt expm} to compute the matrix
exponentials, and backslash to compute the
Katz scores. Other approximations of $e^{\mathcal{A}}$ are 
discussed in section \ref{sec:approx}.

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\textwidth]{evaluesofnewabortion.eps}
\vspace{-.25in}
\caption{Plot of the eigenvalues of the expanded abortion matrix $\mathcal{A}$. }
\label{abortionevalues}
\end{figure}

\subsection{Abortion data set}

The {\em abortion} data set contains $n=2293$ nodes and $m=9644$ 
directed edges.  The expanded matrix $\mathcal{A} =  
\left( \begin{array}{cc} 0 & A \\ A^T & 0  \end{array} \right)$ 
has order $N=2n=4586$ and contains $2m=19288$ nonzeros.  
The maximum eigenvalue of $\mathcal{A}$ is $\lambda_{N} \approx 
31.91$ and the second largest eigenvalue is 
$\lambda_{N-1} \approx 26.04$.  
In this matrix, the largest eigenvalue is fairly 
well-separated from the second largest so that one would 
expect the HITS rankings (which only use information 
from the dominant eigenpair of $\cal A$) to be reasonably close to the 
rankings from $e^{\mathcal{A}}$ (which use information from 
all of the eigenvalues and corresponding eigenvectors). 
A plot of the eigenvalues of the expanded abortion data set 
matrix can be found in Fig.~\ref{abortionevalues}. 
Note the high multiplicity of the zero eigenvalue in this matrix, as well as in the adjacency matrices of the computational complexity and death penalty data sets considered below. Due to this, the numerical rank of the matrix is very low and, as such, the diagonal values of $e^A$ can be estimated using only a few eigenvalues (see section \ref{sec:approx} for further discussion on this).

%\begin{table}[t!]
%\centering
%\caption{Top 10 hubs of the abortion
%web graph, ranked using $[e^{\mathcal{A}}]_{ii}$ and HITS.}
%\label{abortionhubs}
%\begin{tabular}{|c|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & $[e^{\mathcal{A}}]_{ii}$ score & 
%HITS nodes & HITS score  \\
% \hline
%48  & 3.3049e11 & 48  & .0.008968 \\
% \hline
% 1021 & 3.2080e11 & 1006  & .0.008836 \\
%\hline
% 1007 & 3.2080e11 & 1007  & .0.008836 \\
%\hline
%1006 & 3.2080e11 & 1021 & .0.008836 \\
%\hline
%1053 & 3.1958e11 &  1053 &0.008819 \\
%\hline
%1020 & 3.1956e11 & 1020 &0.008819 \\
%\hline
%987 & 3.1891e11 & 960 & 0.008810 \\
%\hline
%990 & 3.1891e11 & 968 &  0.008810  \\
%\hline
%985 & 3.1891e11 & 969 & 0.008810 \\
%\hline
%989 & 3.1891e11 & 970 & 0.008810 \\
%\hline
%\end{tabular}
%\end{table}

%\begin{table}[t!]
%\centering
%\caption{Top 10 authorities of the abortion
%web graph, ranked using $[e^{\mathcal{A}}]_{ii}$ and HITS.}
%\label{abortionauthorities}
%\begin{tabular}{|c|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & $[e^{\mathcal{A}}]_{ii}$ score & 
%HITS nodes & HITS score  \\
% \hline
% 967 & 4.0249e12 & 939  & 0.108667\\
% \hline
% 958 & 4.0249e12 & 958 & 0.108667 \\
%\hline
%939 & 4.02492e12 & 967 & 0.108667 \\
%\hline
%962 & 4.0012e12 & 961& 0.108347 \\
%\hline
%963 & 4.0012e12 & 962 & 0.108347 \\
%\hline
%964 &4.0012e12  & 963 & 0.108347\\
%\hline
%962 & 4.0012e12 & 964 & 0.108347 \\
%\hline
%965 & 4.0012e12 & 965 & 0.108347 \\
%\hline
%966 & 4.0012e12 & 966 & 0.108347 \\
%\hline
%587 & 1.3828e10 & 1582 & 0.002888 \\
%\hline
%\end{tabular}
%\end{table}

\begin{table}[t!]
\centering
\caption{Top 10 hubs of the abortion
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz, 
$e^A$ row sums and Reverse PageRank with $\alpha = 0.85$.}
\label{abortionhubs}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ &
HITS & Katz & $e^A$ rs & RPR \\
 \hline
48  &  48  & 80  & 80 & 125  \\
 \hline
 1021 & 1006  & 1431  & 1431  & 2184  \\
\hline
 1007 & 1007  & 1432 & 1432  & 79  \\
\hline
1006 & 1021 & 1387 & 1426 & 81  \\
\hline
1053 &  1053 & 1388  & 1425 & 48  \\
\hline
1020 & 1020 & 1389  & 1415  & 1424  \\
\hline
987 & 960 & 1397  & 1388  & 1447  \\
\hline
990 & 968 & 1425  & 1389  &  78  \\
\hline
985 & 969 & 1426  & 1397  & 134  \\
\hline
989 & 970 & 1415  & 1387  & 1445  \\
\hline
\end{tabular}
\end{table}

\begin{table}[t!]
\centering
\caption{Top 10 authorities of the abortion
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz,
$e^A$ column sums and PageRank with $\alpha=0.85$.}
\label{abortionauthorities}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ & 
HITS & Katz & $e^A$ cs  & PR  \\
 \hline
 967 & 939 & 1430  & 1430  & 1609 \\
 \hline
 958 & 958 & 1387  & 1387  & 1941 \\
\hline
 939 & 967 & 1425  & 1425  & 1948  \\
\hline
 962 & 961 & 1426  & 1426  & 1608 \\
\hline
 963 & 962 & 1429  & 1417  & 587 \\
\hline
 964 & 963 & 1396  & 1409  & 1610 \\
\hline
 961 & 964 & 1405 & 1429  & 2045  \\
\hline
 965 & 965 & 1406  & 1406 &  317  \\
\hline
966 & 966 & 1409  & 1396 & 2191  \\
\hline
587 & 1582 & 1417  & 1405  & 753  \\
\hline
\end{tabular}
\end{table}



%\begin{table}[t!]
%\centering
%\caption{Top 10 hubs of the abortion web graph, 
%ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, and HITS with exponentiated input.}
%\label{abortionhubs}
%\begin{tabular}{|c|c|c|c|c|}
%\begin{tabular}{|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & HITS nodes & Exp.~input HITS nodes  \\
% \hline
% 48 & 48 &  80 \\
% \hline
%1021 & 1006 & 1431 \\
%\hline
%1007 &  1007 & 1432 \\
%\hline
%1006 & 1021 & 1426 \\
%\hline
%1053 & 1053 & 1425 \\
%\hline
%1020 & 1020 &  1415 \\
%\hline
%987 & 960 & 1388 \\
%\hline
%990 & 968 & 1389 \\
%\hline
%985 & 969 & 1397 \\
%\hline
%989 & 970 & 1387 \\
%\hline
%\end{tabular}
%\end{table}


%\begin{table}[h!]
%\centering
%\caption{Top 10 authorities of the abortion web graph, 
%ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, and HITS with exponentiated input.}
%\label{abortionauthorities}
%\begin{tabular}{|c|c|c|c|c|}
%\begin{tabular}{|c|c|c|}
%\hline
%  $[e^{\mathcal{A}}]_{ii}$ nodes & HITS nodes & Exp.~input HITS nodes \\
% \hline
% 967 &  939 & 1430 \\
% \hline
% 958 & 958 & 1387 \\
%\hline
%939 & 967 & 1425 \\
%\hline
%962 & 961 & 1426 \\
%\hline
%963 & 962 & 1409 \\
%\hline
%964 & 963 & 1417 \\
%\hline
%961 & 964 & 1429 \\
%\hline
%965 & 965 & 1406 \\
%\hline
%966 & 966 & 1396 \\
%\hline
%587 & 1582  & 1405 \\
%\hline
%\end{tabular}
%\end{table}




%\begin{figure}[ht]
%\centering
%\subfigure[Hub rankings]{
%\includegraphics[width=0.8\textwidth]{abortionhubsemilogyrank.eps}
%}
%\subfigure[Authority rankings]{
%\includegraphics[width=0.8\textwidth]{abortionauthsemilogyrank.eps}
%}
%\caption[Optional caption for list of figures]{A semilogy plot of 
%the two normalized hub and authority rankings for each node 
%in the abortion dataset. The rankings using $e^{\mathcal{A}}$ 
%are in blue and the rankings using HITS are in red.}
%\label{abortionplot}
%\end{figure}

The top 10 hubs and authorities of the abortion data set, as determined 
using the diagonal entries of
$e^{\mathcal{A}}$, HITS with constant initial vector, 
the row/column sums of $(I - cA)^{-1}$ (\lq\lq Katz\rq\rq), the 
row/column sums of $e^A$ and Reverse PageRank/PageRank are shown in
Tables \ref{abortionhubs} and \ref{abortionauthorities}. 
We observe that there is a good deal of agreement between the $e^{\cal A}$
rankings
and the HITS ones: indeed, both methods identify the websites 
labeled 48, 1021, 1007, 1006, 1053, 1020 as the top 6 hubs, and
both pick web site $48$ as the top one. Also, there are 7 
web sites identified by both methods as being among the top 10 
authorities. The top authority identified by HITS is ranked third
by $e^{\mathcal{A}}$, and conversely the top authority identified by 
$e^{\mathcal{A}}$ is third in the HITS ranking. 
%In contrast, the exponentiated input HITS algorithm returns 
%completely different rankings. 
%Node 80, which is deemed the best
%hub by the exponentiated input HITS, is ranked 1236 by HITS and
%851 using $e^{\mathcal{A}}$. 
%Exponentiated input HITS 
%ranks node 1430 as the top authority. This node is ranked 731
%by HITS and 429 by $e^{\mathcal{A}}$.
%The odd behavior of
%exponentiated input HITS is due to the fact that the 
%web graph abortion data set contains many nodes with no
%in-links (at least, none that are included in the data set).
%As seen in section \ref{example4}, nodes with no in-links are
%given less value by this algorithm and thus will not be
%reported as top hubs (nor will nodes pointed to by many
%nodes with no in-links be reported as top authorities).
%Since exponentiated input HITS behaved unreliably on the remaining 
%two data sets as well, we do not show
%the corresponding results.
The Katz rankings and those based on $e^A$ show considerable
agreement with one another, but are very different 
from the HITS ones and from those based on $e^{\cal A}$. 
Node 48, which is the top-ranked hub according to HITS
and $e^{\cal A}$, is now not even among the top 100. Conversely,
node 80, which is ranked the top hub by Katz and $e^A$, is not
in the top 100 nodes according to HITS or to $e^{\cal A}$.
This is not too surprising, since the metrics based on $A$ and those
based on $\cal A$ are obtained by counting 
rather different types of graph walks. Finally, for this
network Reverse PageRank
and PageRank return rankings with almost no overlap with any of
the other methods. 

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\textwidth]{evaluesofnewCompComplex.eps}
\vspace{-0.25in}
\caption{Plot of the eigenvalues of the expanded computational 
complexity matrix $\mathcal{A}$. }
\label{compcomplexevalues}
\end{figure}

\subsection{Computational complexity data set}

The {\em computational complexity\/} data set contains 
$n=884$ nodes and $m= 1616$ directed edges.  
The expanded matrix $\mathcal{A}$ has order $N = 2n =1768$ and 
contains $2m=2232$ nonzeros.  The maximum eigenvalue of 
$\mathcal{A}$ is $\lambda_N \approx 10.93$ and the second 
largest eigenvalue is $\lambda_{N-1} \approx 9.86$.  Here, 
the (relative) spectral gap between the first and the second
eigenvalue is smaller than in the previous example; 
consequently, we expect the rankings produced 
using $e^{\mathcal{A}}$ and HITS to be less similar than for 
the abortion data set.  A plot of the eigenvalues of the 
expanded computational complexity data set matrix can be found 
in Fig.~\ref{compcomplexevalues}.

%\begin{table}[t!]
%\centering
%\caption{Top 10 hubs of the computational complexity 
%web graph, ranked using $[e^{\mathcal{A}}]_{ii}$ and HITS.}
%\label{comphubs}
%\begin{tabular}{|c|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & $[e^{\mathcal{A}}]_{ii}$ score & 
%HITS nodes & HITS score  \\
% \hline
% 57 & 2.6518e04 &  57 & 0.276589\\
% \hline
% 17 & 7.2059e02 &  634 & 0.035592 \\
%\hline
%644 & 6.6561e02 &  644 & 0.020557 \\
%\hline
%643 & 6.1256e02 &  721 & 0.018340\\
%\hline
%634 & 5.5558e02 &  643 & 0.017880\\
%\hline
%106 & 4.7486e02 &  554 & 0.014191\\
%\hline
%119 & 4.2791e02 &  632 & 0.013106 \\
%\hline
%529 & 3.8451e02 &  801 & 0.012383 \\
%\hline
%86 & 3.6528e02 & 640 & 0.011566\\
%\hline
%162 & 3.5502e02 & 639 & 0.010893\\
%\hline
%\end{tabular}
%\end{table}

%\begin{table}[h!]
%\centering
%\caption{Top 10 authorities of the computational complexity 
%web graph, ranked using $[e^{\mathcal{A}}]_{ii}$ and HITS.}
%\label{compauthorities}
%\begin{tabular}{|c|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & $[e^{\mathcal{A}}]_{ii}$  score &  
%HITS nodes & HITS score  \\
% \hline
% 1 & 4.8958e03 &  719 & 0.012155\\
% \hline
% 315 & 1.4747e03 &  717 & 0.011501 \\
%\hline
%673 & 8.0015e02 & 727 & 0.009972 \\
%\hline
%148 & 7.3093e02 & 723 & 0.009131\\
%\hline
%719 & 6.6746e02 & 808 & 0.008828\\
%\hline
%717 & 5.8437e02 & 735 & 0.008785\\
%\hline
%2 & 5.5637e02 &  737 & 0.008721 \\
%\hline
%45 & 4.0969e02 & 1 & 0.008550 \\
%\hline
%727 & 4.0315e02 & 722 & 0.008491\\
%\hline
%534 & 3.4473e02 & 770 & 0.008491\\
%\hline
%\end{tabular}
%\end{table}

The top 10 hubs and authorities of the computational complexity 
data set, determined by the various ranking methods, can be found in Tables 
\ref{comphubs} and \ref{compauthorities}.
%We also report the actual scores obtained for these nodes.
 As expected, we see
less agreement between HITS and the diagonals of $e^{\cal A}$. 
Concerning the
hubs, both methods agree that the web site labelled 57 is by far the most
important hub on the topic of computational complexity. However,
the method based on $e^{\cal A}$ identifies as the
second most important hub the web site corresponding to node 17, 
which is ranked only 39th by HITS. The two methods
agree on the next three hubs, but after that they return completely
different results. The difference is even more pronounced for the
authority rankings. The method based on $e^{\cal A}$ clearly
identifies web site 1 as the most authoritative one, whereas
HITS relegates this node to 8th place. The top authority
acording to HITS, web site 719, places 5th in the ranking 
obtained by $e^{\cal A}$. The two methods agree on only two
other web sites as being in the top 10 authorities (717 and 727). 
The Katz rankings and those based on $e^A$ show little overlap
for this data set, although node 57 is clearly considered an
important hub by all measures. A natural question is how much 
these results are affected by the choice of the parameter $c$
used to compute the Katz scores. We found experimentally that,
in contrast to the situation for the other data sets, small
changes in the value of $\alpha$ can significantly affect the Katz
ranking for this particular data set. Changing the value of $c$
to $\alpha = 1/(\rho(A) + 0.3)$ results in hub and authority rankings
that are much closer to those given by the column/row sums of
$e^A$. The potential sensitivity to $\alpha$ is a clear drawback of the Katz-based
approach compared to the methods based on the matrix exponential.
Coming to (Reverse) PageRank, it is interesting to note that for
this data set it provides rankings that are at least in partial
agreement with some of the other measures, especially those
based on $e^A$. Looking at the authority scores, we also notice
a good degree of overlap among all methods, except HITS. Due to 
the small spectral gap, HITS is probably the least reliable of all
ranking methods on this particular data set.

\begin{table}[t!]
\centering
\caption{Top 10 hubs of the computational complexity
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz, 
$e^A$ row sums and Reverse PageRank with $\alpha=0.85$.}
\label{comphubs}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ &
HITS & Katz & $e^A$ rs & RPR \\
 \hline
57  &  57  & 56  & 57  & 57  \\
 \hline
17 & 634  & 709  & 56 & 56  \\
\hline
644 & 644  & 57  & 17  & 17  \\
\hline
643 & 721 & 697  & 51  & 51  \\
\hline
634 & 643 & 705  & 634  & 21  \\
\hline
106 & 544 & 690  & 21  & 11  \\
\hline
119 & 632 & 714  & 255  & 255  \\
\hline
529 & 801 & 708   & 173  &  12  \\
\hline
86 & 640 & 712  & 709 & 13  \\
\hline
162 & 639 & 715  & 45 & 45  \\
\hline
\end{tabular}
\end{table}

\begin{table}%[t!]
\centering
\caption{Top 10 authorities of the computational complexity
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz,
$e^A$ column sums and PageRank with $\alpha=0.85$.}
\label{compauthorities}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ & 
HITS & Katz & $e^A$ cs  & PR  \\
 \hline
 1 & 719 & 688 & 673  & 673 \\
 \hline
 315 & 717 & 685  & 1  & 664 \\
\hline
 673 & 727 & 673  & 664  & 534  \\
\hline
 148 & 723 & 690 & 534  & 45 \\
\hline
 719 & 808 & 56  & 45  & 2 \\
\hline
 717 & 735 & 686  & 473  & 1 \\
\hline
 2 & 737 & 664 & 315  & 376  \\
\hline
 45 & 1 & 1 & 376 &  341  \\
\hline
727 & 722 & 45  & 688  & 50  \\
\hline
534 & 770 & 534  & 599  & 51  \\
\hline
\end{tabular}
\end{table}

%\begin{figure}[ht]
%\centering
%\subfigure[Hub rankings]{
%\includegraphics[width=0.8\textwidth]{compcomplexityhubsemilogyrank.eps}
%}
%\subfigure[Authority rankings]{
%\includegraphics[width=0.8\textwidth]{compcomplexityauthsemilogyrank.eps}
%}
%\caption[Optional caption for list of figures]{A semilogy plot of the two normalized hub and authority rankings for each node in the computational complexity 
%data set. The rankings using $e^{\mathcal{A}}$ are in blue and the rankings using HITS are in red.}
%\label{compcomplexityplot}
%\end{figure}


\begin{figure}[t!]
\centering
\includegraphics[width=0.8\textwidth]{evaluesofnewdeathpenalty.eps}
\vspace{-0.25in}
\caption{Plot of the eigenvalues of the expanded death penalty matrix $\mathcal{A}$. }
\label{deathpenaltyevalues}
\end{figure}

\subsection{Death penalty data set}

The {\em death penalty} data set contains $n=1850$ and $m=7363$ 
directed edges.  The expanded matrix $\mathcal{A}$ has order $N=2n=3700$ 
and contains $m=14726$ nonzeros.  The maximum eigenvalue of 
$\mathcal{A}$ is $\lambda_N \approx 28.02$ and the second largest 
eigenvalue $\lambda_{N-1} \approx 17.68$.  In this case, the largest 
and second largest eigenvalues are quite far apart, and the relative 
gap is larger than in the previous examples.  
A plot of the eigenvalues of the
expanded death penalty matrix can be found in Fig.~\ref{deathpenaltyevalues}.

Due to the presence of a large spectral gap, 
much of the information used in forming the 
rankings of $e^{\mathcal{A}}$ is also used in the HITS ranking,
and we expect the two methods to produce similar results; see
section \ref{hits_comp}.
Indeed, as shown in Table \ref{deathhubs} (hubs) and
Table \ref{deathauthorities} (authorities),    
in this case the top 10 rankings produced by the
two methods are actually identical.  

%\begin{table}%[t!]
%\centering
%\caption{Top 10 hubs of the death penalty web graph, 
%ranked using $[e^{\mathcal{A}}]_{ii}$ and HITS.}
%\label{deathhubs}
%\begin{tabular}{|c|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & $[e^{\mathcal{A}}]_{ii}$ score &  
%HITS nodes & HITS score  \\
% \hline
% 210 & 6.7731e10 &  210 & 0.017562\\
% \hline
% 637 & 3.5805e10 &  637 & 0.012769 \\
%\hline
%413 & 3.2347e10 & 413 & 0.012137 \\
%\hline
%1586 & 1.9676e10 & 1586 & 0.009466\\
%\hline
%552 & 1.7611e10 & 552 & 0.008955\\
%\hline
%462 & 1.1998e10 & 462 & 0.007392\\
%\hline
%930 & 1.1772e10 & 930 & 0.007322 \\
%\hline
%542 & 1.1713e10 & 542 & 0.007303 \\
%\hline
%618 & 1.1163e10 & 618 & 0.007130\\
%\hline
%1275 & 1.0556e10 & 1275 & 0.006933\\
%\hline
%\end{tabular}
%\end{table}

%\begin{table}%[h!]
%\centering
%\caption{Top 10 authorities of the death penalty web graph, 
%ranked using $[e^{\mathcal{A}}]_{ii}$ and HITS.}
%\label{deathautorities}
%\begin{tabular}{|c|c|c|c|}
%\hline
% $[e^{\mathcal{A}}]_{ii}$ nodes & $[e^{\mathcal{A}}]_{ii}$ score &  
%HITS nodes & HITS  score  \\
% \hline
% 4 & 6.4460e10 &  4 & 0.023556\\
% \hline
% 1 & 5.4816e10 &  1 & 0.021723 \\
%\hline
%6 & 3.8451e10 &   6 & 0.018193 \\
%\hline
%7 & 3.0091e10 &  7 & 0.016094\\
%\hline
%10 & 2.8139e10 & 10 & 0.015564\\
%\hline
%16 & 2.6059e10 & 16 & 0.014977\\
%\hline
%2 & 2.5742e10 &  2 & 0.014886 \\
%\hline
%53 & 2.4149e10 &  3 & 0.014418 \\
%\hline
%4 & 1.9763e10 & 44 & 0.013043\\
%\hline
%27 & 1.8591e10 & 27 & 0.012651\\
%\hline
%\end{tabular}
%\end{table}

%\begin{figure}[ht]
%\centering
%\subfigure[Hub rankings]{
%\includegraphics[width=0.8\textwidth]{deathpenaltyhubsemilogyrank.eps}
%}
%\subfigure[Authority rankings]{
%\includegraphics[width=0.8\textwidth]{deathpenaltyauthsemilogyrank.eps}
%}
%\caption[Optional caption for list of figures]{A semilogy plot of the two normalized hub and authority rankings for each node in the death penalty data set. The rankings using $e^{\mathcal{A}}$ are in blue and the rankings using HITS are in red.}
%\label{deathpenaltyplot}
%\end{figure}

Looking at the Katz scores and those based on $e^A$, we see in
this case a great deal of overlap between these two, but almost
completely different rankings compared to HITS and $e^{\cal A}$ 
(although node 210 is clearly an important hub by any
standard). Note that node 1632 is both the top hub and the
top authority according to Katz and to $e^A$. PageRank and
Reverse PageRank show a limited amount of overlap with the other 
measures; nevertheless, nodes 210 and 1632 are also found 
to be important hubs and nodes 1632, 1 and 4 are found to be
authoritative, in agreement with some of the other measures.


\subsection{Stanford web graph}
The {\em wb-cs-stanford} data set from the University of Florida sparse 
matrix collection contains $n=9914$ nodes and $m=36854$ directed edges.  
The expanded matrix $\mathcal{A}$ has order $N=2n=19828$ and contains $m=73708$ nonzeros.  The maximum
eigenvalue of $\mathcal{A}$ is $\lambda_N \approx 38.38$ and the 
second largest is $\lambda_{N-1} \approx 32.12$, hence there
is a sizeable gap.
Tables \ref{stanfordhubs}-\ref{stanfordauthorities} 
report the results obtained with the various ranking schemes.


\begin{table}[t!]
\centering
\caption{Top 10 hubs of the death penalty
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz, 
$e^A$ row sums and Reverse PageRank with $\alpha=0.85$.}
\label{deathhubs}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ &
HITS & Katz & $e^A$ rs & RPR \\
 \hline
210  & 210  & 1632  & 1632  & 210  \\
 \hline
637 & 637  & 133  & 133  & 1632  \\
\hline
413 & 413  & 1671  & 1671 & 70  \\
\hline
1586 & 1586 & 552  & 552 & 95  \\
\hline
552 & 552 & 1651  & 1651  & 135  \\
\hline
462 & 462 & 1673  & 210 & 133  \\
\hline
930 & 930 & 1328  & 1673  & 55  \\
\hline
542 & 542 & 1653  & 1653  & 958  \\
\hline
618 & 618 & 210  & 1328  & 1077  \\
\hline
1275 & 1275  & 1709  & 1709 & 315  \\
\hline
\end{tabular}
\end{table}

\begin{table}%[t!]
\centering
\caption{Top 10 authorities of the death penalty
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz,
$e^A$ column sums and PageRank with $\alpha=0.85$.}
\label{deathauthorities}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ & 
HITS & Katz & $e^A$ cs  & PR  \\
 \hline
 4 & 4 & 1632   & 1632  & 993 \\
 \hline
 1 & 1 & 1662 & 1662  & 667 \\
\hline
 6 & 6 & 1697  & 1697  & 3  \\
\hline
 7 & 7 & 1689  & 1689  & 736 \\
\hline
 10 & 10 & 1653  & 1653  & 735 \\
\hline
 16 & 16 & 1671   & 1671  & 1632 \\
\hline
 2 & 2 & 1675  & 1675 & 42  \\
\hline
 3 & 3 & 1684  & 1684 &  1  \\
\hline
44 & 44 & 798 & 789  & 4  \\
\hline
27 & 27 & 1652  & 1654  & 1212  \\
\hline
\end{tabular}
\end{table}


The first thing to observe is the remarkable agreement
between the HITS, $e^{\cal A}$, Katz, and $e^A$ rankings
of both hubs and authorities. This in stark contrast with
the results for the previous three data sets. Moreover,
many of the nodes that are ranked highly as hubs are also ranked
highly as authorities.  A plausible
explanation of these observations is that the adjacency matrix $A$
for this digraph is much closer to being symmetric than
in the other cases. Indeed, the percentage of \lq\lq bidirectional\rq\rq\
edges in the wb-cs-stanford graph is $47.63\%$; the
corresponding percentages for the abortion, computational
complexity and death penalty graphs are just $2.72\%$, $2.97\%$
and $4.02\%$, respectively. 

Interestingly, the (Reverse) PageRank results are now drastically
different fron the ones provides by all the other measures in nearly
all cases. 
The only (partial) exception is that PageRank finds nodes 6837, 6839 and 6840 
to be among the top 10 authorities; these three
nodes are identified as the three most authoritative ones by the
remaining methods.

\begin{table}[t!]
\centering
\caption{Top 10 hubs of the wb-cs-stanford
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz, 
$e^A$ row sums and Reverse PageRank with $\alpha=0.85$.}
\label{stanfordhubs}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ &
HITS & Katz & $e^A$ rs & RPR \\
 \hline
6562  & 6562 & 6562 & 6562 & 251  \\
 \hline
6838 & 6838  & 6837 &  6837 & 252  \\
\hline
6840 & 6837  & 6838  & 6838 & 253  \\
\hline
6837 & 6839 & 6839  & 6839 & 254  \\
\hline
6839 & 6840 & 6840  & 6840 & 271  \\
\hline
6616 & 6616 & 6669  & 6669 & 2240  \\
\hline
6765 & 6615 & 6668  & 6668  & 2241  \\
\hline
6615 & 6765 & 6670  & 6670  & 2242  \\
\hline
6669 & 6669 & 6616  & 6616  & 2243   \\
\hline
6731 & 6731  & 6615  & 6615 & 348  \\
\hline
\end{tabular}
\end{table}

\begin{table}%[t!]
\centering
\caption{Top 10 authorities of the wb-cs-stanford
web graph, ranked using $[e^{\mathcal{A}}]_{ii}$, HITS, Katz,
$e^A$ column sums and PageRank with $\alpha=0.85$.}
\label{stanfordauthorities}
\begin{tabular}{|c|c|c|c|c|}
\hline
 $[e^{\mathcal{A}}]_{ii}$ & 
HITS & Katz & $e^A$ cs  & PR  \\
 \hline
6837 & 6837 & 6837   & 6837  & 2264 \\
 \hline
 6840 & 6839 & 6839  & 6839  & 8226 \\
\hline
 6839  & 6840 & 6840  & 6840  & 8059  \\
\hline
 6838 & 6838 & 6838  & 6838   & 8057 \\
\hline
 6617 & 6617 & 6573  & 6573  & 4485 \\
\hline
 6615 & 6615 & 6574  & 6575  & 5707 \\
\hline
 6766 & 6614  & 6575  & 6576 & 8225  \\
\hline
 6764 & 6616 & 6576  & 6577 & 6837  \\
\hline
 6616 & 6764 & 6577  & 6578  & 6839  \\
\hline
 6614 & 6766 & 6578  & 6579  & 6840  \\
\hline
\end{tabular}
\end{table}

\section{Approximating the matrix exponential}
\label{sec:approx}

For the purpose of ranking hubs and authorities in a directed network,
only the main diagonal
of $e^{\cal A}$ is required.  As was discussed in Chapter \ref{C:introduction}, Section \ref{sec:iterative_methods}, this can be done without having to
compute {\em all} the entries in $e^{\cal A}$.
When applying the approach based on Gauss quadrature rules
to the $2n\times 2n$ matrix $\cal A$, only matrix-vector
products with $A$ and its transpose are required, 
just like in the HITS algorithm. 
If only the hub scores are wanted, it is also possible to apply the
described techniques to the symmetric matrix $AA^T$ using the
function $f(\lambda) = \cosh(\sqrt{\lambda})$; the same applies
if only the authority scores are wanted, working this time with the matrix
$A^TA$. The problem with this approach is that only estimates (rather
than increasingly accurate lower and upper bound) can be obtained,
due to the fact that the function $f(\lambda) = \cosh(\sqrt{\lambda})$
is not strictly completely monotonic on the positive real axis. We refer
to \cite{BenziGolub} for details. In our experiments we always work
with the matrix $\cal A$, since we are interested in computing both
hub and authority scores. 


\subsection{Test results}
\label{sec:tests}
Accurate evaluation of {\em all} 
the diagonal entries of $e^{\mathcal{A}}$ using quadrature
rules may be too expensive for
truly large-scale graphs. In most applications,
fortunately, it is not necessary to rank all the nodes in
the network: only the top few hubs and authorities are likely
to be of interest. 
When using quadrature rules, the number of quadrature nodes 
(Lanczos iterations) required to correctly rank 
the nodes as hubs or authorities varies and depends on 
both the eigenvalues of $e^{\mathcal{A}}$ and how close 
the diagonal entries are in value.  If the rankings of 
the nodes are
very close, it can take many iterations for the ordering 
to be exactly determined.  However, since estimates for diagonal 
entries are calculated individually, once the top 10 (say) 
nodes have been identified, additional iterations can be 
performed only on these nodes in order to determine their exact ranking.

\begin{table}[t!]
\centering
\caption{The number of iterations necessary for the top 10 
hubs or authorities to be determined (not necessarily in 
the correct order).\\}
%along with the time it took to calculate the bounds.}
\label{intop10}
\vspace{-0.2in}
\subfloat{
\begin{tabular}{|c|c|c|}
\hline
 Dataset & hub (lower bound) & hub (upper bound)  \\
 \hline
Abortion & $>40$ & $>40$ \\
%Comp. Complex. & 3 (1.09s) & 3 (1.09s) \\
Comp. Complex. & 3  & 3  \\
%Death Penalty & 5 (6.10s) & 3 (3.28s) \\
Death Penalty & 5 & 3  \\
Stanford      & 8 & 8 \\
\hline
\end{tabular}
} \\
\subfloat{
\vspace{1in}
\begin{tabular}{|c|c|c|}
\hline
 Dataset & authority (lower bound) & authority (upper bound)  \\
 \hline
%Abortion & 2 (1.57s) & 2 (1.57s)  \\
Abortion & 2  & 2  \\
%Comp. Complex. & 4 (1.36s) & 5 (2.23s)  \\
Comp. Complex. & 4 & 5   \\
%Death Penalty & 4 (3.71 s) & 2 (1.41s) \\
Death Penalty & 4 & 2  \\
Stanford      & 7 & 8 \\
\hline
\end{tabular}
}
\end{table}

Our approach exploits the monotonicity of the Gauss-Radau bounds:
as soon as the lower bound for node $i$ is above the upper bounds for
other nodes, we know that node $i$ will be ranked higher than those
othe nodes. This observation leads to a simple algorithm for
identifying the top-$k$ nodes.
The number of Lanczos iterations per node necessary %(and the time required) 
to identify the top $k=10$ hubs and authorities, using Gauss-Radau 
lower and upper bounds, for the four data sets from section 
\ref{sec:applications} is given in Table \ref{intop10}.  
Our implementation is based on Meurant's Matlab code \cite{mmq},
From the table it can be seen that, in most cases, only 
2-5 iterations per node are needed. An exception is the determination
of the top 10 
hubs of the abortion data set, for which the number 
of iterations is large ($>40$).  
This is due to a cluster of nodes (nodes 960 and 968-990) 
that have nearly identical hub rankings. These nodes' 
scores agree to 15 significant digits.
%the rankings range from 3.189134930290351e11 to 3.189134930290353e11.  
However, for  most applications, if
a subset of nodes are so closely ranked, their 
exact ordering may not be so important.  Table 
\ref{intop30} reports the number of Lanczos iterations neeeded for the 
top $k=10$ hubs and authorities to be ranked at least in the top 30.  
Here, the number of iterations per node needed is never more than 7.
The total cost is thus $O(n)$ Lanczos iterations, again leading to an
$O(n^2)$ overall complexity.
%Using a simple implementation based on Meurant's Matlab code \cite{mmq},
%which has not been optimized,
%the running times 
%to estimate the hub or authority score for all the nodes in the 
%data set using the number of Lanczos iterations listed range from 
%0.58s for the computational
%complexity data set
%to 6.41s for the abortion data set. 
Various enhancements can be used
to reduce costs,  
including the use of sparse-sparse mat-vecs
in the Lanczos iteration, and the
exclusion of nodes with zero out-degree (for hub computations)
and zero in-degree (for authority computations) from the 
top-$k$ calculations. It is also safe to assume that in most 
cases of interest, one can also
exclude nodes with in- and out-degree 1 from the computations,
leading to further savings.

\begin{table}[t!]
\centering
\caption{The number of iterations necessary for the top 10 hubs or 
authorities to be ranked in the top 30.\\}
%, along with the time it took 
%to calculate the bounds.}
\label{intop30}
\vspace{-0.2in}
\subfloat{
\begin{tabular}{|c|c|c|}
\hline
 Dataset & hub (lower bound) & hub (upper bound)  \\
 \hline
%Abortion & 5 (8.54s) & 4 (6.41s) \\
Abortion & 5 & 4 \\
%Comp. Complex. & 2 (.58s) & 2 (.58s) \\
Comp. Complex. & 2 & 2  \\
%Death Penalty & 2 (1.75s) & 2 (1.75s) \\
Death Penalty & 2  & 2  \\
Stanford      & 7 & 4 \\
\hline
\end{tabular}
}\\
\subfloat{
\begin{tabular}{|c|c|c|}
\hline
 Dataset & authority (lower bound) & authority (upper bound)  \\
 \hline
%Abortion & 2 (1.57s) & 2 (1.57s)  \\
Abortion & 2 & 2   \\
%Comp. Complex. & 4 (1.63s) & 2 (.59s)  \\
Comp. Complex. & 4 & 2   \\
%Death Penalty & 2 (1.41s) & 2(1.41s) \\
Death Penalty & 2  & 2 \\
Stanford      & 2 & 4 \\
\hline
\end{tabular}
}
\end{table}

\section{Conclusions and outlook}
In this chapter we have presented a new approach to ranking hubs and
authorities in directed networks using functions of matrices. 
Bipartization is used to transform the original directed network
into an undirected one with twice the number of nodes. The adjacency
matrix of the bipartite graph is symmetric, and this allows the
use of subgraph centrality (and communicability) measures for undirected
networks. We showed that the diagonal entries of the matrix exponential  
provide hub and authority rankings, and we gave an interpretation
for the off-diagonal entries (communicabilities). Unlike HITS, the results are
independent of any starting vectors; and unlike the Katz-based ranking
schemes, there is no dependency on an arbitrary parameter.
% and the proposed method appears
%to be superior to the
%variant of HITS known as \lq\lq exponentiated 
%input HITS\rq\rq, which is unreliable in the presence of a large 
%number of nodes with 0 in-degree (a common occurrence in many real-life
%networks).

Several examples, both synthetic and corresponding to real data sets,
have been used to demonstrate the effectiveness of the proposed
ranking algorithms relative to HITS and to other ranking schemes
based on the matrix resolvent and on the exponential of the
adjacency matrix  of the original digraph. 
Our experiments indicate that our method
results in rankings that are frequently different from
those computed by HITS, at least in the absence of large gaps
between the dominant singular value of the adjacency matrix
$A$ and the remaining ones. This is to be expected,
since our method uses information from all the singular spectrum of the
network, not just the dominant singular triplets.

As usual in this field, there is no simple way to compare different
ranking schemes, and therefore it is impossible to state with certainty
that a ranking scheme will give \lq\lq better\rq\rq\ results than a 
different scheme in practice.
 It is, however, certainly the case that the method
based on the exponential of $\cal A$ takes into account more spectral
information than HITS does; moreover, the rankings so obtained are
unambiguous, in that they do not depend on the choice of an initial
guess or on a tuneable parameter. As we saw, the latter is a weak spot
of Katz-like methods, and a similar case can be made for 
PageRank and Reverse PageRank. 

Compared to HITS, the new technique has
a higher computational cost.
We showed how Gaussian quadrature rules
can be used to quickly identify the top ranked
hubs and authorities for networks involving
thousands of nodes. We note that such schemes
require a symmetric input matrix and are
not readily applicable to nonsymmetric matrices,
since in this case one can only hope for estimates 
instead of lower and upper bounds.  

Future work should include a more efficient implementation
and tests on larger networks. It is likely
that the proposed approach based on Gaussian quadrature will prove to be
too expensive for truly large-scale networks with millions of nodes. We
hope to explore techniques similar to those presented in \cite{bonchietal,fmrr12}
and \cite{dhillon} in order to extend our methodology to truly large-scale
networks. Another relevant question is the study of
the rate of convergence of the
Lanczos algorithm for estimating bilinear forms associated with adjacency
matrices of graphs of different types. %This work will be presented
%elsewhere. 


% the code below specifies where the figures are stored

%: ----------------------- HELP: references
% References can be links to figures, tables, sections, or references.
% For figures, tables, and text you define the target of the link with \label{XYZ}. Then you call cross-link with the command \ref{XYZ}, as above
% Citations are bound in a very similar way with \cite{XYZ}. You store your references in a BibTex file with a programme like BibDesk.

%\section{Motivation}
%%%%%%%%%%%%%%%%%%%%%%%%
%In this thesis we are interested in efficient methods for the numerical solution of saddle point systems arising from finite element discretization of two important classes of fluid dynamic problems, namely the \emph{incompressible Navier-Stokes equations} and the \emph{Brinkman equations}. 

%: ----------------------- 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