\section{Introduction}

The calculation of the power spectrum of cosmological fluctuations
for stochastic initial conditions is an important and challenging task.
It provides the basis for analyzing data on the large-scale structure of
the universe and can lead to constraints on the parameters of
the cosmological model. The spectrum can be computed through
$N$-body simulations. However, these are CPU-intensive, while they
do not provide
an intuitive understanding of the underlying physical processes.
The alternative option is an analytical treatment,
which amounts to solving the collisionless Vlasov-Poisson
equation perturbatively in the density contrast
for appropriate classes of initial conditions \cite{bernardeau0}.
The use of perturbative techniques in this context is complicated by the
growth of non-linearities at short length scales.
Within the strongly non-linear regime, velocity dispersion
and virialization becomes important.

The inability of analytical methods to describe the short scales reliably suggests
to reformulate cosmological perturbation theory
based on an effective description, applicable only above
some length scale.
The description should contain effective
parameters which
absorb the effect of short-scale perturbations that are `integrated out'~\cite{baumann,effective,Pietroni:2011iz}.
In the conventional description of dark matter
as a pressureless ideal fluid, only the lowest moments
of the phase-space distribution are taken into account.
Enlarging this framework by including higher moments results
in viscous transport coefficients that parametrize a non-ideal stress tensor
for the dark matter.
In the following we present heuristic arguments
that support the suggestion that an efficient fluid description of dark matter
must allow for such non-ideal terms. The role of the effective viscosity and sound velocity is to account for
the interaction of large-scale fluctuations (with
wavenumbers $k < k_m$) with the
short-scale ones ($k > k_m$) that are not followed explicitly
in a fluid description limited to $k<k_m$.
In our analysis the UV contributions are incorporated in the large-scale theory
at the level of a one-loop approximation.
Going beyond one-loop can be achieved through
renormalization-group techniques \cite{Max1,floerchinger}.
The approximation that we review here provides an intuitive and rather simple framework, while it still
resolves the
main deficiency of standard perturbation theory (SPT), namely
the strong dependence on the short-scale dynamics that are out of the reach of
analytical techniques.
We review the calculation of the density power
spectrum presented in ref. \cite{blas} in the context of the coarse-grained
theory.
We also present results for the velocity and cross spectra within
the same framework.

\section{The relevant scales}

Our approach to the problem of cosmological matter perturbations assumes
the presence of two scales: the momentum scale $k_\Lx$ at which a fluid
description becomes feasible, and the scale $k_m < k_\Lx$ at which the
description can be based on a small number of
effective parameters. In practice we
can take $k_\Lx \sim 1-3 ~h/$Mpc, corresponding to length scales
$\sim 3-10$ Mpc, and $k_m\sim 0.5-1~h/$Mpc, corresponding to length scales
$\sim 10-20$ Mpc. At the scale $k_m$ the density contrast is of order 1, while
above the scale $k_\Lx$ it is much larger than 1 and the
dynamics involves shell crossing, making a fluid description unfeasible.
The description at $k_m$ must include viscosity terms arising either
through the coarse-graining of the scales $k>k_m$ (effective viscosity), or
fundamental dark matter interactions (fundamental viscosity).
The scales $k>k_\Lx$ correspond to virialized structures, which are
expected to have a negligible effect on long-distance dynamics \cite{Peebles:1980}. It seems
reasonable then to expect that the effective viscosity results mainly
from the integration of the modes with $k_m<k<k_\Lx$. The form of the linear
power spectrum (which falls off roughly as $k^{-3}$) is consistent
with this expectation: the
effective viscosity at $k_m$ is dominated by the modes with $k$ slightly
above this value, while the scales $k>k_\Lx$ give a negligible contribution.
This efficient decoupling of deep UV modes is also supported by numerical
evidence from $N$-body simulations as well as analytical arguments based on
non-perturbative relations of response functions \cite{nishimichi,Garny:2015oya}.

There are several reasons for introducing the scale $k_m$
and distinguishing it from $k_\Lx$:
\begin{itemize}
\item
In the absence of dark-matter interactions,
the form of the effective terms beyond the perfect-fluid
description (such as effective viscosity) at the scale $k_m$ is in a certain sense universal.
They do not depend strongly
on the initial conditions at $k_\Lx$, but only on $k_m$.
If fundamental interactions exist, the fundamental viscosity should
vary only little between $k_\Lx$ and $k_m$, so that it
can be estimated at $k_m$ without loss of accuracy.
\item
The main role of the scale $k_m$ is to act as an UV cutoff for perturbative
corrections in the large-scale theory, thus eliminating
contributions of dubious validity, arising from the modes near and above $k_\Lx$. This effect is seen clearly in the good convergence of our results,
in contrast to standard perturbation theory \cite{blas}.
\item
The scale $k_m$ can be lowered even further than the range that we are
considering, if one would like to focus on the effective theory for length scales
much larger than the range of baryon acoustic oscillations (BAO).
In such a case, the renormalization-group improvement of the
effective description is crucial \cite{floerchinger}.
\end{itemize}

\section{Estimate of the effective viscosity}

The validity of the fluid description for cold dark matter cannot be established
similarly to systems close to thermal equilibrium. The crucial elements supporting this description are the small dark-matter velocity and
the finite age of the universe, during which dark matter particles can
drift over a finite distance, much smaller than the Hubble radius.
A more quantitative discussion of this point is given in ref. \cite{baumann}.
In the absence of dark-matter interactions, one starts from the Vlasov-Poisson
equations
and writes the particle phase-space density as
\be
f({\bf x},{\bf p},\tau)=f_0(p)[1+\delta_f({\bf x},p,\hat{\bf p},\tau)],
\label{phase} \ee
with $\tau$ the conformal time.
The Fourier modes of the perturbation are then expanded in terms of
Legendre polynomials $P_n$:
\be
\delta_f({\bf k},p,\hat{\bf p},\tau)=
\sum_{n=0}^\infty (-i)^n (2n+1)\,\delta_f^{[n]}({\bf k},p,\tau)
\, P_n(\hat{\bf k}\cdot \hat{\bf p}),
\label{legendre} \ee
with unit vectors $\hat{\bf k}$, $\hat{\bf p}$.
The Vlasov-Boltzmann equation implies a hierarchy of evolution equations for the
higher moments:
\be
\frac{d\delta_f^{[n]}}{d\tau}=k v_p \left[\frac{n+1}{2n+1}\delta_f^{[n+1]}
-\frac{n}{2n+1}\delta_f^{[n-1]} \right],
~~~~~~~~~n\geq2,
\label{momev} \ee
with $v_p=p/{am}$ the particle velocity.
If the typical value of
$v_p$ is sufficiently small, the higher moments can be neglected.

A quantitative constraint can be obtained if one assumes that the
particle velocity
is of the order of the fluid velocity $v$ at small length scales.
We emphasize that this assumption is valid only for cold dark matter,
for which the particle velocity is induced by the gravitational fields
generated by the growth of large-scale structure.
The assumption is supported by the estimated typical velocities
of dark matter particles in galaxy haloes.\footnote{Due to the cosmological expansion, the massive particles
that constitute the cold dark matter have negligible velocity
originating in their thermal motion before freeze-out.
Warm dark-matter particles would have much larger velocities, which
makes the application of the fluid description to this case problematic.}
At the comoving momentum scale $k$ we expect $[\theta/\Hc]^2\sim k^3 P^L(k)$,
with $P^L(k)$ the linear power spectrum evaluated on the growing mode,
$\theta=\vec{\bf k}\cdot \vec{\bf v}$ the velocity divergence, and $\Hc=(1/a)(da/d\tau)$.
The linear spectrum scales roughly as $k^{-3}$ above a scale $\sim k_m$,
so that $k^3 P^L(k)$ is roughly constant, with a value of order 1 today.
Its time dependence is given by $D_L^2$,
 with $D_L$ the linear growth factor.
 We expect then that the maximal
particle velocity can be identified with the fluid velocity at the
scale $k_m$, and is roughly
\be
v_p\sim \frac{\Hc}{k_m}\, D_L.
\label{partvel} \ee
The time $\tau$
available for the higher moments to grow
is $\sim 1/\Hc$.
Therefore, the dimensionless factor characterizing the
growth of higher moments is $k v_p/\Hc\sim D_L\, k/k_m$.
Scales with $k\gg k_m/D_L$ require the use
of the whole Boltzmann hierarchy and cannot be treated through a fluid
description that takes into account only the lowest moments.
In practice, the validity
of the fluid description extends beyond $k_m$ at all times,
as the initial values of the
higher moments are much smaller than 1, while
$D_L\leq 1$.
Therefore, on large scales, cold dark matter can be described as a cosmological
fluid, not because microscopic interaction rates are sufficiently large to maintain local thermal equilibrium, but because the lifetime
of the universe is too short for dark matter
to deviate strongly from local equilibrium.

The above picture can lead to an estimate of the effective viscosity
$\eta$ that
must be attributed to
the dark-matter fluid because of the particle velocities
that deviate from the uniform motion.
In general, we expect $\eta/(\rho+p)\sim l_{\rm free}v_p$, with
$l_{\rm free}$ the mean free path. For the effective viscosity we
can estimate $l_{\rm free}\sim v_p /H$, with $H=\Hc/a$.
In this way we obtain for the kinematic viscosity of the
dark-matter fluid
\be
\nu_{\rm eff} \Hc =\frac{\eta_{\rm eff}}{\rho_m \, a}\Hc
\sim l_{\rm free}v_p H\sim \frac{\Hc^2}{k_m^2} D^2_L,
\label{etaeff} \ee
where we have neglected the pressure.
The numerical value is small, of the order of $10^{-6}$ today,
and even smaller at higher redshifts. However, viscosity plays a significant role in the context of the coarse-grained theory.

We can also estimate the viscosity induced by possible fundamental interactions between the dark-matter particles.
For interacting dark matter, with number density $n$, mass $m$ and
cross section $\sigma$, the mean free path is $l_{\rm free}\sim 1/(n\sigma)$.
We expect then fundamental viscosity
\be
\nu_{\rm fund} \Hc =\frac{\eta_{\rm fund}}{\rho_m \, a}\Hc
\sim l_{\rm free}v_p H \sim \frac{1}{n\sigma} \frac{H^2}{k_m}a D_L
\sim \frac{m}{\sigma} \frac{8\pi G_N}{3k_m} \frac{a D_L }{\Omega_m}.
\label{etafund} \ee
The time dependence, given by the last factor, indicates that the fundamental
viscosity drops very quickly at early times.
We point out that the above estimate is valid only for cold dark matter,
for which the typical particle velocity is induced by the late-time
growth of structure. For warm dark matter, the typical velocity
and the resulting viscosity would evolve differently with time.

We note that that the fundamental viscosity of cold dark matter does not
grow indefinitely for decreasing $\sigma$. The mean free path
of dark matter particles cannot exceed $v_p/ H\sim a D_L /k_m$.
As a result, for
values of $\sigma$ such that $1/(n\sigma)$ exceeds $a D_L /k_m$,
the particle scattering has no effect.
This implies that the fundamental viscosity
cannot take values larger than roughly $(\Hc^2/k_m^2) D_L^2$, and
is always subleading or comparable to the effective viscosity.

As a final remark, we would like to emphasize that the constraints
derived in refs. \cite{skordis1,kunz,skordis2} for the non-ideal properties
of dark matter are not applicable to the case of cold dark matter that we
are considering. In these works the viscosity and sound velocity are
assumed to stay constant or even grow with redshift, in contrast to
the redshift-dependence indicated by our eqs. (\ref{etaeff}) and (\ref{etafund}).
In our picture, the particle velocity is induced by the growth of structure
and is a late-time phenomenon. The assumptions of refs. \cite{skordis1,kunz,skordis2} can by justified only within a warm-dark-matter
scenario, in which the particles have velocities originating in their
primordial motion.

\section{One-loop determination of the effective viscosity}

For the construction of the coarse-grained theory we follow the
Wilsonian approach. A detailed description, including a discussion of the
 renormalization-group flow of the coarse-grained theory, is given
in ref. \cite{floerchinger}. The main idea is that the
effective parameters of the large-scale theory result from the
integration of the UV modes. For the problem at hand, the effective
viscosity and sound velocity of the
dark-matter fluid are given by one-loop expressions that depend on
the dimensionful scale
\be
\sigma^2_{dm}({\etat})=\frac{4\pi}{3} \int_{k_m}^\infty dq\, P^L(q,\etat)
=\frac{4\pi}{3}  D^2_L(\etat) \int_{k_m}^\infty dq\, P^L(q,0),
\label{sigmadm} \ee
with $\etat=\ln a=-\ln(1+z)$ and $D_L(\etat)$ the linear growth factor:
$\delta(\bk,\etat)=D_L(\etat)\,\delta(\bk,0)$ on the growing mode.
For a spectrum that scales as $1/k^3$, the integral is dominated by the region near $k_m$.

A precise definition of the effective couplings can be given through
 the propagator of theory.
One first defines the doublet
\begin{equation}
 \left(
\begin{array}{c}
\varphi_{1}(\textbf{k},\etat)\\ \\ \varphi_{2}(\textbf{k},\etat)
\end{array}
\right)
=\left(
\begin{array}{c}
\delta_\bk(\etat)\\ \\-\dfrac{\theta_{\bk}( \etat)}{\mathcal{H}}
\end{array}
\right),
\label{oneloopprop}
\end{equation}
and the propagator
\be
G_{ab}(\bk, \tilde\eta, \tilde\eta^\prime)\, \delta^{(3)}(\bk - \bk^\prime) = \left\langle  \frac{\delta \varphi_a(\bk,\tilde\eta)}{\delta \varphi_b(\bk^\prime, \tilde\eta^\prime)} \right\rangle\, .
\label{propagator}
\ee
The power spectrum is defined as
\be
 \langle\varphi_{a}(\bk, \etat)\varphi_{b}(\bk', \etat)\rangle
= \delta^{(3)}\,({\bk+\bk'}) \, P_{ab}({\bf{k}}, \etat).
\label{spectrum}
\ee

We concentrate on an Einstein-de Sitter (EdS) background, for which
$D_L(\tilde\eta)=\exp(\tilde\eta)$ and the linear propagator for the growing mode in
the perfect-fluid theory is
\be
g_{ab}(\tilde\eta-\tilde\eta^\prime) = \frac{e^{\tilde\eta-\tilde\eta^\prime}}{5} \begin{pmatrix} 3 && 2 \\ 3 && 2 \end{pmatrix}.
\label{linprop} \ee
The {\it one-loop} contribution to the low-energy
propagator with wavenumber $k \to 0$, coming from modes with $q\geq k_m$, is given by
\be
\delta g_{ab}(\bk, \tilde\eta,\tilde\eta^\prime)
=  - k^2 e^{\tilde\eta-\tilde\eta^\prime} \sigma_{dm}^2(\tilde\eta)
\begin{pmatrix} \frac{61}{350} && \frac{61}{525} \\ \frac{27}{50} && \frac{9}{25} \end{pmatrix}\,
\simeq - k^2 \, e^{\eta} \, \sigma_{dm}^2(\etat) \frac{3}{50}
\begin{pmatrix}3 \cdot 0.968 && 2\cdot 0.968 \\ 9 && 6 \end{pmatrix}\, ,
\label{onlooppro}
\ee
where we have taken the limit $\etat\gg \etat'$ for the second equality.
On the other hand, the {\it linear} propagator of a theory with small kinematic
viscosity and sound velocity of the form
\be
\nu\Hc=\frac{\eta}{\rhoz\, a}\Hc=\frac{3}{4}\, \beta_\nu \, e^{2 \etat}\, \frac{\Hc^2}{k^2_m}\, ,
~~~~~~~~~~
c^2_s = \frac{\dpp}{\drho}= \beta_s\,  e^{2 \etat}\, \frac{\Hc^2}{k^2_m}\, ,
\label{deff} \ee
contains a contribution
\be
\delta g_{ab}(\bk,\etat) = -\frac{k^2}{k^2_m} (\beta_\nu+\beta_s)
\, \frac{e^{3 \etat}}{45} \,
\begin{pmatrix} 3 && 2 \\ 9 && 6 \end{pmatrix}
\label{corrprop} \ee
in addition to the perfect-fluid term (\ref{linprop}),
where we have kept only the leading contribution for $\etat\gg\etat'$.

We can now derive precise expressions for the effective parameters of the
large-scale theory, by identifying the
{\it linear} contribution to the propagator (\ref{corrprop})
for the {\it viscous} theory,
with the {\it one-loop} correction (\ref{onlooppro}) of the {\it perfect-fluid} propagator. This can be achieved with 1\% accuracy, and we
find
\be
\beta_s+\beta_\nu=\frac{27}{10} \, k^2_m\sigma^2_{dm}(0).
\label{matching} \ee
It must be emphasized that there are
no free parameters in our approach, as the effective couplings are uniquely
determined. Moreover, the deep UV region, which is out of the reach of
perturbation theory, gives a negligible contribution to the integral of eq.
(\ref{sigmadm}), consistently with the expectation of decoupling of virialized
structures. The dominant contribution comes from the region slightly above
$k_m$, for which perturbation theory is expected to give reliable results.

By matching heuristically viscous fluid dynamics to
the one-loop contribution (\ref{corrprop}), effective couplings were explored in \cite{blas}. This procedure
can be obtained as a limiting case from the renormalization-group formulation of Ref.~\cite{floerchinger}.
The latter approach also allows one to break the degeneracy between $\beta_\nu$ and $\beta_s$~\cite{floerchinger}.
The numerical results that we present below depend only on the sum $\beta_\nu+\beta_s$ to a very good approximation.

\begin{figure*}
\centering
\includegraphics[width=6cm]{A-PoverNbodyVisc}\qquad
\includegraphics[width=6cm]{A-PoverNbodySPT}
\caption{Density power spectrum at redshift $z=0$,
normalized to the $N$-body power spectrum~\cite{Kim:2011ab}. The left panel
shows the results obtained in the viscous
description, for various values of the scale
$k_m$.
The coloured lines show the linear spectra,
the open and filled symbols the one- and two-loop
spectra.
For comparison, the right panel
shows the corresponding one- and two-loop results
in standard perturbation theory. Here $\Lambda$ is a sharp cutoff in $k$-space. The grey shaded region shows the estimated uncertainty of the $N$-body simulation. Figure taken from Ref.~\cite{blas}. }
\label{fig-2}
\end{figure*}

\section{Density and velocity power spectra}

The procedure that we outlined in the previous section for an EdS background
can be generalized to the $\Lambda$CDM case with an appropriate change of
variables \cite{blas}.
We present here results for the density, velocity and mixed power spectra for
a large-scale viscous theory with parameters determined through
eq. (\ref{matching}).
The non-linear spectra for the low-energy theory can be computed
through perturbation theory, with the crucial modification that (i) all
momentum integrals possess a UV cut off $k_m$, and (ii) the Euler equation
includes pressure and viscosity contributions. The presence of the cutoff
prevents the inclusion of spurious effects from the deep UV, so that
fast convergence is expected. It must be emphasized that the precise value
of the scale $k_m$ should not influence the final results for physical
quantities, such as the power spectra: the explicit dependence of the effective
pressure and viscosity on $k_m$ compensates the $k_m$-dependence introduced
by the cutoff. In some sense,
the scale $k_m$ may be viewed as
a technical device that guarantees the good convergence of our scheme.
Nevertheless, because of the perturbative truncation we expect a residual dependence
that is a measure of the remaining theoretical uncertainty.

In order to test the accuracy of our results we compare them
to the power spectrum extracted from $N$-body simulations \cite{Kim:2011ab}.
In the left panel of fig.\,\ref{fig-2}
we show the density power spectrum at redshift $z=0$,
obtained from solving the viscous fluid equations
and averaging over an initially Gaussian random distribution of the density field.
The different symbols show our results obtained for the four values
$k_m=0.4, 0.6, 0.8, 1.0\, h/$Mpc. In addition, open symbols are the one-loop results,
and the filled symbols show the power spectrum up to two-loops. The loop expansion is
perturbative in the density contrast. However, the effect of the time-dependent
pressure and viscosity on the propagation of density and velocity perturbations is
taken into account non-perturbatively.
We normalized all curves to the result of the $N$-body simulation.
The grey-shaded area corresponds to an estimate of the error of the simulation.
The agreement between $N$-body data and the two-loop results is at the percent level for $k\lesssim 0.2\, h/$Mpc.
In addition, the variation of the results when changing the scale $k_m$ is mild, at the $\pm 2 \%$ level for
$k\lesssim 0.2\, h/$Mpc. Our results can be compared to the ones obtained in
standard perturbation theory (SPT) \cite{bernardeau0}, depicted in the
right panel of \ref{fig-2}.  For SPT we show the one- and two-loop contributions obtained when imposing a
cutoff $\Lambda$ in $k$-space. The cutoff dependence in SPT is sizable at two loops, and amounts to $\pm 5\%\ (\pm 10\%)$
for $0.4\, h/$Mpc$<\Lambda<1\, h/$Mpc ($0.4\, h/$Mpc$<\Lambda\lesssim 5\, h/$Mpc).
It is apparent that, as compared to SPT, the viscous
description reduces the uncertainty generated by
the treatment of UV modes significantly (left panel in fig.\,\ref{fig-2}).
The remaining uncertainty within the viscous framework is consistent with the expected impact of effects not captured by the fluid dynamical
description employed here \cite{Pueblas:2008uv}.

We also consider
the power spectrum for the velocity divergence and the cross power spectrum.
The velocity spectrum at $z=0$ is shown in fig.\,\ref{fig-3}, where we again compare the viscous description (left panel) with the result obtained
in SPT (right panel) for various values of $k_m$ and $\Lambda$, respectively.
All curves are normalized to the density power spectra corresponding to the same approximation.
Fig.\,\ref{fig-4} shows a similar comparison for the cross power spectrum.
The extraction of power spectra for the velocity divergence from $N$-body simulations is less robust than
for the density, because of complications arising when extracting a continuous velocity field from a
discrete set of particles. We show the results obtained from two different methods \cite{Jennings:2012ej, Hahn:2014lca} by the two thick red lines
in figs.\,\ref{fig-3} and \ref{fig-4}. The results are sensitive to certain parameters of the simulation, especially the resolution,
and we show the uncertainty quoted in refs. \cite{Jennings:2012ej, Hahn:2014lca} as shaded regions. Within this relatively large uncertainty,
the results obtained in the viscous description (left panels in figs.\,\ref{fig-3} and \ref{fig-4}) at two loops are again in agreement up to about $k\lesssim 0.2\, h/$Mpc. In addition, the variation when changing
$k_m$ is much smaller than the corresponding cutoff dependence in SPT (right panels in figs.\,\ref{fig-3} and \ref{fig-4}).
The residual $k_m$-dependence yields a quantitative estimate of the theoretical uncertainty.
The difference of the perturbative prediction and the $N$-body results, as well
as the theoretical uncertainty, increase quickly with $k$.
Given that the predictions within the scheme discussed here do not require to fit any free parameters, the level of agreement between the two approaches for $k\lesssim 0.2\, h/$Mpc is remarkable.

\section{Conclusions}

We conclude that the description of dark matter as
a viscous fluid discussed here yields a robust framework for predicting power spectra for $k\lesssim 0.2\, h/$Mpc at $z=0$
without the need to adjust any free parameters. The framework can be
constructed and extended through a formal analysis
based on the renormalization group.
We refer the reader to ref. \cite{floerchinger} for details.

\begin{figure*}
\centering
\includegraphics[width=6cm]{A-P22OverP11Viscz0}\qquad
\includegraphics[width=6cm]{A-P22OverP11SPTz0}
\caption{Power spectrum of the velocity-divergence, normalized to the density power spectrum, for the viscous theory (left)
and standard perturbation theory (right). The various lines show the dependence on $k_m$ and on $\Lambda$, respectively, as
in fig.\,\ref{fig-2}. The red lines are $N$-body results ((J12) \cite{Jennings:2012ej}, (HAA14) \cite{Hahn:2014lca}), and the shaded
region is the quoted uncertainty of the velocity power spectra extracted from $N$-body data.
}
\label{fig-3}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[width=6cm]{A-P12OverP11Viscz0}\qquad
\includegraphics[width=6cm]{A-P12OverP11SPTz0}
\caption{Cross density-velocity power spectrum,
corresponding to the same
approximations as in fig.\,\ref{fig-3}.}
\label{fig-4}
\end{figure*}


\section{INTRODUCTION}\label{S:Introduction}

	A network of $n$ objects can be denoted by $N=(V,E)$, where $V=\{v_1,v_2,...,v_n \}$ is the set of nodes, and $E$ is the set of edges between nodes. In this article, we focus on symmetric weighted networks represented by an $n \times n$ adjacency matrix, $A$, where the element $A_{ij}$ measures the strength of the relationship between nodes $v_i$ and $v_j$.

	Traditionally, statistical network analysis focuses on modeling \textit{observed} network structure (e.g.,  highway systems or electrical transmission grids).  In this situation, nodes are well defined and the physical links between these nodes can be directly observed \cite{Hiller01,Newman11}. However, in some fields of research the explicit network structure may not be observable.  This is particularly true of the social sciences.  In these fields, the observable data are groups of individuals and it is presumed that the groups are the result of some latent network. The fundamental task is to construct a network from such data.

	A familiar example is provided by Wasserman and Faust \cite{Wasserman94} who introduce the problem of inferring the relationships among a collection of children based on their attendance at birthday parties. In this introductory dataset, the children represent nodes in the network and the birthday parties represent subsets of children.

	In our paper, a collection of nodes observed in the same sample is called a \textit{group} and a dataset of these observations is referred to as \textit{grouped data}.  In Wasserman and Faust's example, each party defines a group and the set of all parties is the grouped data. Two individuals are said to \textit{co-occur} if they appear in the same group.

	One common technique used to estimate $A$ from grouped data is to count the number of times that a pair of nodes appears in the same group \cite{Zachary77,Freeman89,Wasserman94,Kolaczyk09,Brent11}.  Frequently, a threshold is applied to this count to create a symmetric unweighted undirected adjacency matrix; however, it has been shown that the characteristics of networks inferred by this technique are sensitive to the choice of threshold \cite{Choudhury10}.  Therefore, we adopt a generalized version of the  inter-citation frequency \cite{Kolaczyk09} which measures the number of times a pair of nodes is observed to co-occur in the dataset.  We refer to this measure as the \textit{co-occurrence matrix}.

	An alternative technique, called the \textit{half weight index} \cite{Cairns87,Bejder98}, estimates $A$ by the frequency that two nodes co-occur given that one of them is observed.  This addresses a shortcoming of the co-occurrence matrix in which nodes that appear rarely can be estimated to have a weak relationship even though the relationship is quite strong \cite{Cairns87,Voelkl11}.

	The co-occurrence matrix and half weight index both have probabilistic interpretations. The co-occurrence matrix estimates the probability that two nodes will be observed together.  The half weight index estimates the probability that two nodes will be observed together given that one of them is observed. However, the probability or conditional probability of co-occurrences is not equivalent to the probability of connection of nodes. In fact, neither of these techniques describe the process which leads to the generation of the observed groups. It is unclear how the descriptive statistics relate to the grouped data in these methods.

	We propose a model-based approach to infer the latent network from grouped data. We refer to this model of group formation as the \textit{Hub Model} because each observed group is assumed to be brought together by a hub node.  A group generated by the Hub Model is illustrated in Figure \ref{F:Sample Hub}.

	\begin{figure}[H]
		\begin{center}
			\fbox{
				\includegraphics{Example_Star-eps-converted-to}
			}
		\end{center}
		\caption{The generating mechanism of the Hub Model is demonstrated on a group of 10 nodes. In the observed sample, nodes $v_1,\dots,v_6$ are members of the group while nodes $v_7,\dots,v_{10}$ are not members of the group.  The observed group is the result of the hub node, $v_1$, bringing together nodes $v_2,\dots,v_6$.}
		\label{F:Sample Hub}
	\end{figure}

	The Hub Model belongs to the family of finite mixture models. Such mixture models  and their variants have been applied in many different situations including text classification  \cite{Carreira00}, topic models \cite{Anandkumar15}, fingerprint identification \cite{Vretos12}, and product recommendation \cite{Colace15}.

	One of the advantages of Hub Models is that we are able to give the strength of the relationship $A_{ij}$ a definition which is both mathematically clear and at the same time practical to researchers.  In the Hub Model, $A_{ij}$ is defined as the probability that node $v_i$ will include node $v_j$ when $v_i$ is the hub node of  a group. The formal definition of the Hub Model will be given in Section \ref{S:Model}.

	As a graphical introduction to the performance of Hub Models, consider Figure \ref{F:Shining_True} where relationship strength is represented by the width of the links between nodes.  In this example there is a pair of nodes, $v_1$ and $v_2$, which never directly bond to each other.  Despite the aversion between these two nodes, they share a number of strong relationships with common nodes and therefore often co-occur.  In Figure \ref{F:Shining_O}, we can see that the co-occurrence matrix mistakenly assigns a relatively strong relationship to nodes $v_1$ and $v_2$. In Figure \ref{F:Shining_H}, the half weight index arrives at a very similar conclusion as the co-occurrence matrix.  In both Figures \ref{F:Shining_O} and \ref{F:Shining_H}, the non-existent relationship between nodes $v_1$ and $v_2$ is actually estimated to be stronger than all other relationships. By contrast, results of the Hub Model in Figure \ref{F:Shining_A} clearly capture the latent network of the population.

	\begin{figure}[H]
		\centering
		\begin{subfigure}[b]{0.2\textwidth}
			\includegraphics[width=\textwidth]{Shining_True-eps-converted-to}
			\caption{\small Truth}
			\label{F:Shining_True}
		\end{subfigure
		~
		\begin{subfigure}[b]{0.2\textwidth}
			\includegraphics[width=\textwidth]{Shining_O-eps-converted-to}
			\caption{\small Co-occurrence}
			\label{F:Shining_O}
		\end{subfigure}
		~
		\begin{subfigure}[b]{0.2\textwidth}
			\includegraphics[width=\textwidth]{Shining_H-eps-converted-to}
			\caption{\small Half Weight Index}
			\label{F:Shining_H}
		\end{subfigure}
		~
		\begin{subfigure}[b]{0.2\textwidth}
			\includegraphics[width=\textwidth]{Shining_A-eps-converted-to}
			\caption{\small Hub Model}
			\label{F:Shining_A}
		\end{subfigure}
		\caption{Comparison of Estimation Techniques}\label{F:Shining}
	\end{figure}

	To the best of our knowledge, there have been very limited attempts to infer latent networks by model-based approaches. Rabbat et al.\cite{Rabbat08} provides an application for telecommunication networks. Rabbat et al. modeled the process of group formation as a random walk from a source node to a terminal node.  The nodes along the path were subjected to an unknown permutation to account for the lack of order information.  Treating permutations as missing data, Rabbat et al. employed a \textit{Monte Carlo EM} algorithm based on importance sampling to estimate the parameters of the model. This model assumed a distinctly different process of group formation.

	In the following sections we present a formal description of the grouped data structure, review existing techniques, and define Hub Models.  Then we study Hub Model identifiability and prove that a symmetry constraint on the adjacency matrix is a sufficient condition for identifiability. We propose an EM algorithm to solve the maximum likelihood estimator of Hub Model.  We then evaluate the performance of the proposed method by simulation studies. We apply the Hub Model to infer the relationships between United States Senators in the 110$^{th}$ Congress based on their co-sponsorship of legislation, the relationships among the characters of the 18$^{th}$ century Chinese novel, \textit{Dream of the Red Chamber}, and apply the Hub Model to floral dispersion over North America.  We close with a discussion that returns to the subject of identifiability along with a phenomenon of ``self-sparsity'' of the Hub Model estimators.

	\section{GROUPED DATA}

	\subsection{Data Structure}

	For a population of $n$ individuals, $V=\{v_1,\dots, v_n\}$, we observe  $T$ subsets of the global population, $\{V^{(t)}|V^{(t)} \subseteq V, t=1,...,T \}$.
	Each observed subset $V^{(t)}$ can be coded as an $n$ length row vector $G^{(t)}$ where:

	\[ G_i^{(t)} = \left\{
	\begin{array}{l l}
	1 & \quad \textnormal{if $v_i\in V^{(t)}$}\\
	0 & \quad \textnormal{if $v_i \notin V^{(t)}$}
	\end{array} \right.\]

	The full set of observations is denoted by a $T \times n$ matrix, $G$.  The $t^{th}$ row of $G$ is $G^{(t)}$.

	\subsection{Existing Methods}

	Inferring latent networks from grouped data relies on descriptive statistics which count the number of times that two nodes are observed together.  We focus on two popular techniques which estimate probabilities of individual behavior.

	An intuitive and computationally simple measure of grouped data is the \textit{co-occurrence matrix}.  Versions of this technique appear throughout the literature under many names and notations including: \textit{capacity matrix} \cite{Zachary77}, \textit{sociomatrix} \cite{Wasserman94}, \textit{inter-citation frequency} \cite{Kolaczyk09}, \textit{cocitation matrix} \cite{Newman11}, and \textit{strength} \cite{Brent11}.

	A co-occurrence matrix, $O$, is an $n\times n $ symmetric matrix, defined by:
	\begin{equation}
	O=\frac{G' G}{T},
	\end{equation}

	which estimates the frequency that the nodes $v_i$ and $v_j$ are observed in the same group.

	One shortcoming of the co-occurrence matrix is that it estimates the probability that two nodes \textit{will be observed} to co-occur in a given observation.  That is, if two nodes have a strong relationship, but appear in the dataset infrequently, the co-occurrence matrix will estimate a low probability that the two nodes \textit{will be observed} to co-occur.

	As an example, consider four nodes $v_1,\dots, v_4$ and the grouped data represented in Table \ref{T:O_vs_H_Example}.
	{\linespread{1.0}
		\begin{table}[hbt!]
			\begin{center}
				\begin{tabular}{|c  | c c c c  |}
					\hline

					&\multicolumn{4}{c|}{Node} \\
					\hline
					Event & $v_1$ & $v_2$ & $v_3$ & $v_4$  \\
					\hline
					1 & 1 & 0 & 0 & 0 \\
					2 & 1 & 1 & 0 & 0 \\
					3 & 1 & 1 & 0 & 0 \\
					4 & 1 & 0 & 1 & 1 \\
					5 & 0 & 1 & 1 & 1 \\
					\hline

				\end{tabular}
			\end{center}
			\caption{Notional Grouped Data}
			\label{T:O_vs_H_Example}
		\end{table}
	}
	For this dataset, both $O_{1,2}=\frac{2}{5}$ and $O_{3,4}=\frac{2}{5}$.  However, notice that every time node $v_3$ is present node $v_4$ is also present.  A researcher may conclude that there is some aspect of the relationship between nodes $v_3$ and $v_4$ which has been understated.

	As an alternative, the \textit{half weight index} estimates the probability that two nodes will be observed to co-occur given that one of them is observed \cite{Cairns87}.

	The half weight index has been introduced in a number of equivalent forms \cite{Dice45}.  Computationally, the most direct form is:
	\begin{equation}\label{E:HWI1}
	H_{ij} = \frac{2\sum_t G_i^{(t)} G_j^{(t)}}{\sum_t G_i^{(t)}+\sum_t G_j^{(t)}}.
	\end{equation}

	The half weight index estimates the frequency that the nodes $v_i$ and $v_j$ are observed in the same group given that one of them is observed.

	Returning to the example in Table \ref{T:O_vs_H_Example}, we can see that $H_{1,2}=\frac{4}{7}$ while $H_{3,4}=\frac{4}{4}$.  Therefore, the half weight index provides different information about the latent network than the co-occurrence matrix.

	\section{HUB MODELS}\label{S:Model}

	\subsection{Generating Mechanism}

	Hub Models (HM) assume that each group is a star subgraph on the global population. The hub node of the star subgraph connecting $G^{(t)}$ is represented by an $n$ length row vector, $S^{(t)}$, where

	\[ S_i^{(t)} = \left\{
	\begin{array}{l l}
	1 & \quad \text{if $v_i$ the hub node of sample $t$},\\
	0 & \quad \text{otherwise}.
	\end{array} \right.\]

	There is one and only one element of $S^{(t)}$ that is equal to 1.

	Under the Hub Model, each group $G^{(t)}$ is independently generated by a two step process.

	\begin{enumerate}
		\item The hub node is drawn from a multinomial distribution with parameter $\rho=(\rho_1,...,\rho_n)$, i.e., $\rho_i=\mathbb{P}(S_i^{(t)}=1)$. The following constraint applies to $\rho_i$:
		\begin{equation}\label{E:USM Constraint}
		\sum_i \rho_i=1.
		\end{equation}

		\item The hub node, $v_i$, chooses to include $v_j$ in the group with probability $A_{ij}$, i.e., $A_{ij}=	\mathbb{P}(G_j^{(t)}=1|S_i^{(t)}=1)$.
	\end{enumerate}

	In most practical applications, the hub node of each group $S^{(t)}$ is missing data. We focus on this case in this article. When necessary,we  refer to the model with both $S$ and $G$ being observed as the Known Hub Model (KHM).

	Since many existing techniques for inferring networks produce a symmetric adjacency matrix, we assume $A_{ij}=A_{ji}$. The symmetry condition will be shown to ensure the identifiability of the Hub Model when $S$ is unobserved(Section \ref{S:Ident}).

	Further, we assume that the hub node will always include itself in the group, i.e. $A_{ii}=1$ for all $i$.

	This generating mechanism implies that each observed group is independent of every other observed group.  In particular, $G^{(t)}$ is not a transformation of $G^{(t-1)}$ and the order in which groups are observed contains no information about the relationships between group members.  Researchers often collect data in such a way to ensure this property \cite{Bejder98}.

	\subsection{Likelihood of the Hub Model}

	From now on, when we mention the likelihood of the Hub Model, we mean the marginal likelihood of the group observations $G$.

	Under the HM, the probability of $G^{(t)}$ has the form of a finite mixture model with $n$ components:
	\begin{equation}\label{E:HM ProbOfGt}
	\mathbb{P}(G^{(t)}|A,\rho)=\sum_{i=1}^n \rho_i G_i^{(t)} \prod_{j}{A_{ij}^{G_j^{(t)}} (1-A_{ij})^{1-G_j^{(t)}}}.
	\end{equation}

	By taking the log of the product of individual observed groups, the log likelihood function for the full set of observations is:
	\begin{equation}\label{E:HM LL}	\mathcal{L}(G|A,\rho)=\sum_t\textnormal{log}\Big[\sum_{i=1}^n \rho_i G_i^{(t)} \prod_{j}{A_{ij}^{G_j^{(t)}} (1-A_{ij})^{1-G_j^{(t)}}}\Big].
	\end{equation}

	Solving the MLE of HM is an optimization problem with the equality constraints $\sum_i \rho_i=1$, and $A_{ij}=A_{ji}$ for all $i$ and $j$. From \eqref{E:HM LL}, we denote the log likelihood function as $\mathcal{L}(G|A,\rho)$.  This gives the following Lagrange function:
	\begin{equation}\label{E:HM Objective Function}
	\Lambda(G|A,\rho)=\mathcal{L}(G|A,\rho)-\lambda_o  [(\sum_i \rho_i)-1  ]-\sum_{i<j}\lambda_{ij}(A_{ij}-A_{ji}).
	\end{equation}

	The log likelihood does not have a closed-form solution for the MLE.  Instead we will derive estimating equations which can be incorporated into an Expectation Maximization algorithm.  Before doing so we investigate the identifiability problem of the Hub Model.

	\subsection{Ensuring Identifiability}\label{S:Ident}

	A basic requirement for any model is \textit{identifiability}.  For Hub Models, this means for any two sets of parameters $\{ A,\rho \}$ and $\{A^*, \rho^*\}$:
	\begin{align}
	\mathbb{P}(G=g|A,\rho)=\mathbb{P}(G=g|A^*,\rho^*) \hspace{2mm} \forall g \implies A=A^*, \rho=\rho^* .
	\end{align}

	We first show that if no constraint is put on $A$, the model is unidentifiable.

	This can be demonstrated by the following simple counterexample.  Consider a network of size $n=4$ with $\{A,\rho\}$ defined in Table \ref{T:IdentEx}. Notice that nodes $v_1$ and $v_2$ both always produce the same group while nodes $v_3$ and $v_4$ produce a different group.

	{\linespread{1.0}
		\begin{table}[hbt!]
			\begin{center}
				\begin{tabular}{| l | c | l l l l|}
					\cline{2-6}
					\multicolumn{1}{c}{}  &  \multicolumn{5}{|c|}{$A_{ij}$}\\
					\cline{3-6}
					\multicolumn{1}{c}{}  & \multicolumn{1}{|c}{}    & \multicolumn{4}{|c|}{$j$}\\
					\hline
					$\rho_i$ &  $i$ & 1 & 2 & 3 & 4\\
					\hline
					0.25 & 1 & 1 & 1 & 1 & 0 \\
					0.25 & 2 & 1 & 1 & 1 & 0 \\
					0.25 & 3 & 0 & 1 & 1 & 1 \\
					0.25 & 4 & 0 & 1 & 1 & 1 \\
					\hline
				\end{tabular}
			\end{center}
			\caption{Example of a Set of Parameters Which are Not Identifiable}
			\label{T:IdentEx}
		\end{table}
	}
	The probability of $G$ has the form:

	\[ \mathbb{P}(G=g|A,\rho) = \left\{
	\begin{array}{l l}
	\frac{1}{2} & \quad  g=\{1,1,1,0\}, \\
	\frac{1}{2} & \quad  g=\{0,1,1,1\}, \\
	0 & \textnormal{otherwise} .
	\end{array} \right.\]

	There are an infinite number of parameters yielding the same distribution, but a simple alternative is as follows.  Let $\rho^*=(0.5, 0 , 0.5,0)$, leave the first and third rows of $A$ unchanged but let all other components  of $A^*$ assume arbitrary values. Obviously, we have $	\mathbb{P}(G=g|A,\rho)=\mathbb{P}(G=g|A^*,\rho^*)$ for all $g$. This counterexample demonstrates that the model requires an additional condition to be identifiable.

	The following theorem shows that symmetry of $A$ is a sufficient condition for identifiability.

	\newtheorem{HM}{Theorem}

	\begin{HM}\label{HM}
		Let $A$ and $A^*$ be symmetric adjacency matrices with $A_{ii}=A_{ii}^*=1$ for all $i$, $A_{ij}<1$ and $A_{ij}^*<1$ for all ${i\neq j}$. If  $\mathbb{P}(g|A,\rho)=\mathbb{P}(g|A^*,\rho^*)$ for all $g$, then $\{A,\rho\}=\{A^*,\rho^*\}$.
	\end{HM}

	Let $g^x$ and $g^y$ denote the singleton groups which consist only of nodes $v_x$ and $v_y$, respectively.  Further, let $g^{xy}$ denote the group representing the pair of $v_x$ and $v_y$.

	From \eqref{E:HM ProbOfGt} the probability of the singletons is:
	\begin{align}
	\mathbb{P}(g^x|A,\rho)&=\rho_x (1-A_{xy}) \prod_{j\ne\{x,y\}}(1-A_{xj})\label{E:Prob_gx}\\
	\mathbb{P}(g^y|A,\rho)&=\rho_y (1-A_{xy}) \prod_{j\ne\{x,y\}}(1-A_{yj})\label{E:Prob_gy}.
	\end{align}

	In \eqref{E:Prob_gy} we have taken advantage of the symmetry of $A$ to replace $A_{yx}$ with $A_{xy}$.

	Now, we consider the probability of $g^{xy}$.
	\begin{align}
	\mathbb{P}(g^{xy}|A,\rho)& =\rho_x A_{xy} \prod_{j\ne\{x,y\}}(1-A_{xj})
	+ \rho_y A_{xy} \prod_{j\ne\{x,y\}}(1-A_{yj}) \notag \\
	& =A_{xy} \Big[\rho_x \prod_{j\ne\{x,y\}}(1-A_{xj})
	+ \rho_y  \prod_{j\ne\{x,y\}}(1-A_{yj})\Big]  \notag \\
	& =A_{xy} \Big[ \frac{\mathbb{P}(G=g^x|A,\rho)}{(1-A_{xy})}
	+ \frac{\mathbb{P}(g^y|A,\rho)}{(1-A_{xy})}\Big]  \notag \\
	&	=\frac{A_{xy}}{(1-A_{xy})} \Big[ \mathbb{P}(g^x|A,\rho)
	+ \mathbb{P}(g^y|A,\rho)\Big], \label{E:Prob_gxy}
	\end{align}
	which implies that:
	\begin{equation}\label{E:HM_MoM}
	A_{xy}=\frac{\mathbb{P}(g^{xy}|A,\rho)}{\mathbb{P}(g^x|A,\rho)
		+ \mathbb{P}(g^y|A,\rho)+\mathbb{P}(g^{xy}|A,\rho)}.
	\end{equation}

	Therefore, $A_{xy}=A_{xy}^*$ for all $x$ and $y$.

	To complete the proof, consider an arbitrary node $v_x$ which appears as a singleton represented by $g^x$:
	\begin{equation}
	\mathbb{P}(g^x|A,\rho)=\rho_x  \prod_{j \ne x }(1-A_{xj}).
	\end{equation}
	If $A_{xy}=A_{xy}^*$ for all $x$ and $y$ and $\mathbb{P}(g|A,\rho)=\mathbb{P}(g|A^*,\rho^*)$ for all $g$, then:
	\begin{equation}
	\rho_x  \prod_{j \ne x }(1-A_{xj})=\rho_x^*  \prod_{j \ne x }(1-A_{xj})
	\end{equation}
	and it is easy to see that $\rho_x=\rho_x^*$ for all $x$.

	$\square$

	\paragraph{Remarks}

	Before proceeding we would like to make three remarks about Hub Models.

	\begin{enumerate}

		\item The generating mechanism for Hub Models is equivalent to a finite mixture model of multivariate Bernoulli random variables. In general, such a model is not identifiable \cite{Teicher61}. This shortcoming does not prevent such models from being useful in many applications. For example, when dealing with classification problems where the researcher only has to identify which component density an observation came from, this type of mixture can be effectively used \cite{Carreira00}.  In such a situation, the individual parameters of the multivariate Bernoulli random variables are not of interest.  However, the issue of identifiability presents a challenge in network inference because we are specifically interested in the individual parameters of $A$.

		\item Equation \eqref{E:HM_MoM} suggests a method of moments estimator for $A_{xy}$ based on frequencies of doubletons and singletons. However,  this estimator requires that the probability of doubletons and singletons be estimated accurately so this technique would be very inefficient for most real cases, because small groups appear infrequently in many datasets. Therefore, we will continue to consider the MLE which presumably uses all available information.

		\item It is worth noticing that even though symmetry of $A$ is a natural assumption in network analysis, it is only a sufficient condition for identifiability according to Theorem \ref{HM}. For future work, we will explore other assumptions to ensure identifiability and ultimately find a necessary and sufficient condition.

	\end{enumerate}

	\subsection{Estimating Equations}

	The maximum likelihood estimator of HM does not have a closed-form solution for the parameters. In this section, we derive estimating equations for the conditions that the MLE $\{\hat{A}, \hat{\rho}\}$ must satisfy. Then we will show that solving these equations iteratively is equivalent to an EM algorithm. The details of the EM algorithm will be given in the next section.

	We begin by taking the derivative of \eqref{E:HM LL} with respect to $A_{xy}$ and $A_{yx}$.

	\begin{align}
	\frac{\partial \Lambda(G|A,\rho)}{\partial A_{xy}}& =\frac{\partial \mathcal{L}(G|A,\rho)}{\partial A_{xy}}-\lambda_{xy}=0 \mbox{ if } x<y, \\
	\frac{\partial \Lambda(G|A,\rho)}{\partial A_{yx}}& =\frac{\partial \mathcal{L}(G|A,\rho)}{\partial A_{yx}}+\lambda_{xy}=0 \mbox{ if } x>y.
	\end{align}
	Therefore,
	\begin{equation}\label{E:Symmetric Constraint}
	\frac{\partial \mathcal{L}(G|A,\rho)}{\partial A_{xy}}=-\frac{\partial \mathcal{L}(G|A,\rho)}{\partial A_{yx}}.
	\end{equation}

	We now focus on the derivative of the log likelihood function of \eqref{E:HM LL}:
	\begin{equation}\label{E:HM Derivative of LL1}
	\sum_t \frac{\rho_x G_x^{(t)}\Big[\prod_{j \ne y}{A_{xj}^{G_j^{(t)}} (1-A_{xj})^{1-G_j^{(t)}}}\Big]\frac{\partial}{\partial A_{xy}} \Big( A_{xy}^{G_y^{(t)}}(1-A_{xy})^{1-G_y^{(t)}} \Big) }{\sum_{i=1}^n \rho_i G_i^{(t)} \prod_{j}{A_{ij}^{G_j^{(t)}} (1-A_{ij})^{1-G_j^{(t)}}}}.
	\end{equation}

	Note that the derivative in the numerator of \eqref{E:HM Derivative of LL1} is equal to 1 if node $v_y$ is in observation $G^{(t)}$, and $-1$ if $v_y$ is not in the observation.  We represent this by the function:

	\[ \gamma(G_y^{(t)}) = \left\{
	\begin{array}{l l}
	1 & \quad \textnormal{if $G_y^{(t)}=1$},\\
	-1 & \quad \textnormal{if $G_y^{(t)}=0$}.
	\end{array} \right.\]
	Therefore,
	\begin{equation}\label{E:HM Derivative of LL2}
	\frac{\partial}{\partial A_{xy}} \mathcal{L}(G|A,\rho)=\sum_t  \frac{\Big[ \rho_x G_x^{(t)}\prod_{j \ne y}{A_{xj}^{G_j^{(t)}} (1-A_{xj})^{1-G_j^{(t)}}}\Big]\gamma(G_y^{(t)})}{\sum_{i=1}^n \rho_i G_i^{(t)} \prod_{j}{A_{ij}^{G_j^{(t)}} (1-A_{ij})^{1-G_j^{(t)}}}}.
	\end{equation}

	The denominator of \eqref{E:HM Derivative of LL2} is simply the probability of $G^{(t)}$ (see \eqref{E:HM ProbOfGt}).  In addition, the term in brackets can be made equal to $\mathbb{P}(G^{(t)},S_x=1)$ by multiplying $A_{xy}^{G_y^{(t)}}(1-A_{xy})^{(1-G_y^{(t)})}$. To conserve space, we suppress $\{A,\rho\}$ going forward.  This gives:
	\begin{equation}\label{E:HM Derivative of LL3}
	\frac{\partial}{\partial A_{xy}} \mathcal{L}(G)=\sum_t \frac{\gamma(G_y^{(t)}) \mathbb{P}(G^{(t)}, S_x^{(t)}=1)}{A_{xy}^{G_y^{(t)}}(1-A_{xy})^{(1-G_y^{(t)})} \mathbb{P}(G^{(t)})}.
	\end{equation}

	This equation can be further simplified by observing that  $\frac{\mathbb{P}(G^{(t)},S_x^{(t)}=1)}{\mathbb{P}(G^{(t)})}$ is equivalent to $\mathbb{P}(S_x^{(t)}=1|G^{(t)})$:
	\begin{equation}\label{E:HM Derivative of LL Final}
	\frac{\partial}{\partial A_{xy}} \mathcal{L}(G)=\sum_t \frac{\gamma(G_y^{(t)}) \mathbb{P}(S_x^{(t)}=1|G^{(t)})}{A_{xy}^{G_y^{(t)}}(1-A_{xy})^{(1-G_y^{(t)})}}.
	\end{equation}

	Plugging \eqref{E:HM Derivative of LL Final} into \eqref{E:Symmetric Constraint}, we get:
	\begin{equation}
	\sum_t \frac{\gamma(G_y^{(t)}) \mathbb{P}(S_x^{(t)}=1|G^{(t)})}{A_{xy}^{G_y^{(t)}}(1-A_{xy})^{(1-G_y^{(t)})}}=
	- \sum_t \frac{\gamma(G_x^{(t)}) \mathbb{P}(S_y^{(t)}=1|G^{(t)})}{A_{yx}^{G_x^{(t)}}(1-A_{yx})^{(1-G_x^{(t)})}}.
	\end{equation}

	By applying symmetry and breaking the summations, this becomes:
	\begin{align}
	\sum_{t:G_y^{(t)}=1} \frac{\mathbb{P}(S_x^{(t)}=1|G^{(t)})}{A_{xy}}
	&- \sum_{t:G_y^{(t)}=0} \frac{\mathbb{P}(S_x^{(t)}=1|G^{(t)})}{1-A_{xy}}\notag\\
	= - \sum_{t:G_x^{(t)}=1} \frac{\mathbb{P}(S_y^{(t)}=1|G^{(t)})}{A_{xy}}
	&+ \sum_{t:G_x^{(t)}=0} \frac{\mathbb{P}(S_y^{(t)}=1|G^{(t)})}{1-A_{xy}}.
	\end{align}
	With some simple algebra, it is easy to see that:
	\begin{equation}\label{E:HM A_hat}
	\hat{A}_{xy}=\frac{\sum_t G_y^{(t)} \mathbb{P}(S_x=1|G^{(t)})+\sum_t G_x^{(t)}\mathbb{P}(S_y=1|G^{(t)})}{\sum_t \big[\mathbb{P}(S_x=1|G^{(t)})+\mathbb{P}(S_y=1|G^{(t)})\big]}.
	\end{equation}

	It is worth repeating that \eqref{E:HM A_hat} is not a closed form solution for $\hat{A}_{xy}$.  This is because the right hand side of the equation depends on $\hat{A}_{xy}$.

	We next derive the estimating equation for $\hat{\rho}$.
	By taking the derivative of \eqref{E:HM Objective Function} with respect to $\rho_x$, we get the following:
	\begin{align}
	\frac{\partial}{\partial \rho_x} \Lambda(G)&=\sum_t \frac{ G_x^{(t)} \prod_{j}{A_{xj}^{G_j^{(t)}} (1-A_{xj})^{1-G_j^{(t)}}}}{\mathbb{P}(G^{(t)})}-\lambda_o\notag\\
	&=\sum_t \frac{ \mathbb{P}(G^{(t)},S_x^{(t)}=1)}{\rho_x \mathbb{P}(G^{(t)})}-\lambda_o\notag\\
	&=\frac{1}{\rho_x}\sum_t \mathbb{P}(S_x^{(t)}=1|G^{(t)})-\lambda_o.
	\end{align}
	Solving this equation for zero, we obtain:
	\begin{equation}\label{E:HM rho with lambda}
	\rho_x=\frac{1}{\lambda_o}\sum_t \mathbb{P}(S_x^{(t)}=1|G^{(t)}).
	\end{equation}
	Using the constraint on $\rho$ \eqref{E:USM Constraint}, we get:
	\begin{equation}\label{E:HM rho_hat}
	\hat{\rho}_x=\frac{\sum_{t=1}^T \mathbb{P}(S_x^{(t)}=1|G^{(t)})}{T}.
	\end{equation}

	\section{EM ALGORITHM}

	The previous section derived estimating equations which depended on the probability $\mathbb{P}(S_x^{(t)}=1|G^{(t)})$. This implies a fairly intuitive algorithm updating $\{\hat{A},\hat{\rho}\}$ and $\mathbb{P}(S_x^{(t)}=1|G^{(t)}) $ iteratively, which can be fitted into the general framework of EM algorithm.

	The key technique of any EM algorithm is to formulate a complete data model then solve the model as if some data is observed and other data is missing.  In this case, the Known Hub Model serves as the complete data model, $G$ is the observed data, and $S$ is the missing data. Each iteration of the EM algorithm consists of an expectation step followed by a maximization step \cite{McLachlan08}.

	\vspace{4mm}
	\textbf{E-Step}
	\vspace{4mm}

	Since the log likelihood function of the complete data model is linear in the unobserved data, $S_i^{(t)}$, the E-Step (on the $(m+1)^{th}$ iteration) simply requires calculating the current conditional expectation of $S_i^{(t)}$ given the observed data, $G^{(t)}$ (see \cite{McLachlan08} for detailed explanation).
	\begin{align}\label{E:ProbKisCenter}
	E[S_x^{(t)}|G^{(t)}]&=\mathbb{P}(S_x^{(t)}=1|G^{(t)})\notag\\
	&=
	\frac{ \rho_x G_x^{(t)} \prod_{j}{A_{xj}^{G_j^{(t)}} (1-A_{xj})^{1-G_j^{(t)}}}}{\sum_{i=1}^n \rho_i G_i^{(t)} \prod_{j}{A_{ij}^{G_j^{(t)}} (1-A_{ij})^{1-G_j^{(t)}}}}
	\end{align}

	\vspace{4mm}
	\textbf{M-Step}
	\vspace{4mm}

	The M-Step replaces $\mathbb{P}(S_x^{(t)}=1|G^{(t)})$ on the right hand side of \eqref{E:HM A_hat} and \eqref{E:HM rho_hat} with $E[S_x^{(t)}|G^{(t)}]$ from \eqref{E:ProbKisCenter}.

	\vspace{4mm}

	\vspace{4mm}
	\textbf{Algorithm}
	\vspace{4mm}

	Algorithm \ref{A:EM Algorithm} illustrates the details of the algorithm for the Hub Model.

	Several standard techniques are used to improve the performance of the EM algorithm.  Firstly, we run the EM algorithm 10 times with different staring points and choose the solution with the highest likelihood. Secondly, we limit the number of iterations applied to a starting point. This second treatment is based in part on the observation that when this algorithm has a bad starting point, it will take a very long time to converge and the point that it converges to is not close to the maximum.  As a final step, we treat any $\hat{A}_{xy}\le 10^{-4}$ as $\hat{A}_{xy}=0$.  We apply this finishing step to remove clutter from the returned solutions.

	\begin{spacing}{0.6}
		\begin{algorithm}[H]
			\KwData{G}
			\KwResult{$\hat{A}, \hat{\rho}$}
			Initialize:\\
			$\mathcal{L}(G|\hat{A})=-\infty$\\
			\For{rep=1 to 10}{
				Initialize:\\
				$\hat{A}_{ij}^{(0)}=unif(0,1) \hspace{2 mm} \forall \{i,j\}$\\
				$X_i=unif(0,1) \hspace{2 mm} \forall i$\\
				$\hat{\rho}_i^{(0)}=\frac{X_i}{\sum_k X_k}$\\
				$\Delta \mathcal{L}(G|A^{(0)})=10^4$\\
				counter=1\\
				\While{$\vert \frac{\Delta \mathcal{L}(G|A^{(m+1)})}{ \mathcal{L}(G|A^{(m)})} \vert >10^{-4} \textnormal{and counter}<100$}{
					\textbf{E-Step}\\
					$\>$ Update $\mathbb{P}(S_k^{(t)}=1|G^{(t)})$ by Equation \ref{E:ProbKisCenter}\\
					\textbf{M-Step}\\
					$\>$ Update $A^{(m+1)}$ by Equation \ref{E:HM A_hat}\\
					$\>$ Update $\rho^{(m+1)}$ by Equation \ref{E:HM rho_hat}\\
					$\Delta \mathcal{L}(G|A^{(m+1)})=\mathcal{L}(G|A^{(m+1)})-\mathcal{L}(G|A^{(m)})$\\
					counter=counter$+1$\\
				}
				\If{$\mathcal{L}(G|A^{(m+1)})>\mathcal{L}(G|\hat{A})$}{
					\eIf{$\hat{A}_{ij} \le 10^{-4}$}{
						$\hat{A}_{ij}=0$
					}
					{$\hat{A}_{ij}=A_{ij}^{(m+1)}$}
				}
			}
			\caption{Expectation Maximization Algorithm for the Hub Model}
			\label{A:EM Algorithm}
		\end{algorithm}
	\end{spacing}

	\section{SIMULATION}\label{S:Simulation}

	In order to perform simulations, we generate $\{A,\rho\}$ using the following techniques.

	For $\rho$, we select $n$ i.i.d. random numbers uniformly, $X_i$, and divide each random number by the sum of all $X_i$'s.  That is, $\rho_i=\frac{X_i}{\sum_i X_i}$.

	We use a two step process to generate $A$.  First, we create a symmetric unweighted undirected random graph on $n$ nodes using the configuration model with a power law degree distribution.  We refer to this unweighted graph as the \textit{structure} of the network.  Each edge in the graph is then assigned a relationship strength with a beta distribution,
	\[ A_{ij} = \left\{
	\begin{array}{l l}
	Beta(\alpha, \beta) & \textnormal{if there is an edge between $v_i$ and $v_j$}\\
	0 & \textnormal{otherwise}
	\end{array} \right.\]
	We simply let $A_{ji}=A_{ij}$ since we assume $A$ is symmetric.

	In Table \ref{T:Converge}, we consider four different networks with $n=\{10,20,30,50\}$.  The latent network for each scenario was generated with power $2$, and beta distribution $\alpha=1$, and $\beta=4$. For each network, we generated 100 datasets with $T=50,000$ and estimated the latent network using subsets of the dataset ranging from 100 observations to the full dataset.

	For each combination of $n$ and $T$, we report four different measures of performance.

	The first performance measure is the average run time to estimate $\{A,\rho\}$ on a Intel Pentium CPU G2030 at 3.00 GHz with 4.00GB of RAM.

	Next, we measured the ability of $\hat{A}$ to correctly identify the structure of the latent network.  To do this we define true positives and true negatives as follows:
	\begin{align}
	TP_{ij} &=I(  A_{ij} =1 \mbox{ and } \hat{A}_{ij} =1) \\
	TN_{ij} &=I ( A_{ij} =0 \mbox{ and } \hat{A}_{ij} =0).
	\end{align}

	Recall that we apply a threshold to $\hat{A}_{ij}$ in the EM algorithm such that $\hat{A_{ij}}$ is zero if it is less than $10^{-4}$.

	Using this notation, we assess the accuracy of the estimated network structure using:
	\begin{equation}
	Accuracy = \frac{\sum_{i<j} TP_{ij}+TN_{ij}}{\binom{n}{2}}
	\end{equation}

	\cite{Han11}.

	To measure the difference between $\{A,\rho\}$ and $\{\hat{A},\hat{\rho}\}$, we calculate the mean absolute error (MAE):
	\begin{align}
	MAE(A)&=\frac{1}{\binom{n}{2}} \sum_{i<j} |\hat{A}_{ij}-A_{ij}|\\
	MAE(\rho)&=\frac{1}{n} \sum_{i} |\hat{\rho}_{i}-\rho_{i}|.
	\end{align}

	The first observation from Table \ref{T:Converge} is simply that as the number of observations increases, the average error in the estimates tends to decline and the variability of the estimate also decreases.

	A more interesting observation concerns the runtime of the estimates.  For increasing values of $n$, the average runtime for a single estimate also increases.  However, for all but the case of $n=50$, as the number of observations increases, the run time actually decreases.  This is  due to reductions in the number of iterations necessary to achieve convergence as the dataset gets larger.

	The most interesting observation is that while the number of parameters is of order $O(n^2)$ the accuracy of the estimates is high even when there is very little data available.  This is due to a property of the model which we call \textit{self-sparsity} and discuss in more detail in Section \ref{S:Discussion}.

	{\linespread{1.0}
		\begin{table}[H]
			\begin{center}
				\resizebox{5in}{!}{
					\begin{tabular}{|c | c c | c c | c c |}
						\hline
						& \multicolumn{6}{c|}{$n=10$}\\
						\hline
						& Avg Run& Avg & Avg & StDev & Avg & StDev \\
						Obs &  Time (sec) & Accuracy & MAE($\rho$) & MAE($\rho$) & MAE($A$) & MAE($A$)\\
						\hline
						100   & 0.0953 & 0.9329 & 0.0254 & 0.0073 & 0.0132 &  0.0035\\
						200   & 0.0901 & 0.9656 & 0.0181 & 0.0051 & 0.0089 &  0.0022 \\
						500   & 0.0624 & 0.9933 & 0.0115 & 0.0030 & 0.0056 &  0.0015\\
						1000  & 0.0573 & 0.9998 & 0.0081 & 0.0021 & 0.0039 &  0.0011 \\
						2000  & 0.0542 & 1.0000 & 0.0057 & 0.0014 & 0.0027 &  0.0007\\
						5000  & 0.0504 & 1.0000 & 0.0036 & 0.0009 & 0.0017 &  0.0004 \\
						10000 & 0.0494 & 1.0000 & 0.0025 & 0.0007 & 0.0012 &  0.0003\\
						20000 & 0.0500 & 1.0000 & 0.0018 & 0.0005 & 0.0008 &  0.0002\\
						50000 & 0.0504 & 1.0000 & 0.0011 & 0.0003 & 0.0005 &  0.0001\\
						\hline
					\end{tabular}
				}
				\vspace{2mm}
				\resizebox{5in}{!}{
					\begin{tabular}{|c | c c | c c | c c | }
						\hline
						& \multicolumn{6}{c|}{$n=20$} \\
						\hline
						& Avg Run& Avg & Avg & StDev & Avg & StDev \\
						Obs &  Time (sec) & Accuracy & MAE($\rho$) & MAE($\rho$) & MAE($A$) & MAE($A$) \\
						\hline
						100   & 0.2190 & 0.9468 & 0.0194 & 0.0031 & 0.0157 &  0.0054\\
						200   & 0.2121 & 0.9635 & 0.0135 & 0.0024 & 0.0103 &  0.0038\\
						500   & 0.1766 & 0.9756 & 0.0079 & 0.0016 & 0.0048 &  0.0013\\
						1000  & 0.1557 & 0.9861 & 0.0055 & 0.0010 & 0.0030 &  0.0006\\
						2000  & 0.1400 & 0.9945 & 0.0040 & 0.0007 & 0.0021 &  0.0004\\
						5000  & 0.1445 & 0.9997 & 0.0026 & 0.0004 & 0.0013 &  0.0002\\
						10000 & 0.1464 & 1.0000 & 0.0017 & 0.0003 & 0.0009 &  0.0001\\
						20000 & 0.1898 & 1.0000 & 0.0012 & 0.0002 & 0.0006 &  0.0001\\
						50000 & 0.2166 & 1.0000 & 0.0008 & 0.0001 & 0.0004 &  0.000\\
						\hline
					\end{tabular}
				}
				\vspace{2mm}

				\resizebox{5in}{!}{
					\begin{tabular}{|c | c c | c c | c c |}
						\hline
						& \multicolumn{6}{c|}{$n=30$} \\
						\hline
						& Avg Run& Avg & Avg & StDev & Avg & StDev \\
						Obs &  Time (sec) & Accuracy & MAE($\rho$) & MAE($\rho$) & MAE($A$) & MAE($A$) \\
						\hline

						100   & 0.4226 & 0.9577 & 0.0153 & 0.0021 & 0.0128 &  0.0025 \\
						200   & 0.4306 & 0.9720 & 0.0109 & 0.0017 & 0.0088 &  0.0020 \\
						500   & 0.3419 & 0.9884 & 0.0066 & 0.0010 & 0.0045 &  0.0006 \\
						1000  & 0.3003 & 0.9943 & 0.0047 & 0.0007 & 0.0030 &  0.0004 \\
						2000  & 0.2479 & 0.9960 & 0.0033 & 0.0005 & 0.0021 &  0.0003 \\
						5000  & 0.2286 & 0.9981 & 0.0021 & 0.0004 & 0.0013 &  0.0002 \\
						10000 & 0.2326 & 0.9994 & 0.0014 & 0.0002 & 0.0009 &  0.0001 \\
						20000 & 0.2440 & 0.9999 & 0.0010 & 0.0002 & 0.0006 &  0.0001 \\
						50000 & 0.2481 & 1.0000 & 0.0007 & 0.0001 & 0.0004 &  0.0001 \\

						\hline
					\end{tabular}
				}
				\vspace{2mm}
				\resizebox{5in}{!}{
					\begin{tabular}{|c | c c | c c | c c | }
						\hline
						& \multicolumn{6}{c|}{$n=50$}\\
						\hline
						& Avg Run& Avg & Avg & StDev & Avg & StDev \\
						Obs &  Time (sec) & Accuracy & MAE($\rho$) & MAE($\rho$) & MAE($A$) & MAE($A$) \\
						\hline

						100   & 0.7410 & 0.9457 & 0.0121 & 0.0012 & 0.0192 &  0.0029\\
						200   & 0.9810 & 0.9562 & 0.0089 & 0.0011 & 0.0142 &  0.0025\\
						500   & 1.2460 & 0.9760 & 0.0054 & 0.0006 & 0.0069 &  0.0014\\
						1000  & 1.2914 & 0.9828 & 0.0037 & 0.0004 & 0.0040 &  0.0004\\
						2000  & 1.3745 & 0.9838 & 0.0025 & 0.0003 & 0.0026 &  0.0002\\
						5000  & 1.8633 & 0.9850 & 0.0016 & 0.0002 & 0.0016 &  0.0001\\
						10000 & 2.5726 & 0.9874 & 0.0011 & 0.0001 & 0.0011 &  0.0001\\
						20000 & 10.9601 & 0.9941 & 0.0008 & 0.0001 & 0.0008 &  0.0001\\
						50000 & 21.0469 & 0.9976 & 0.0005 & 0.0001 & 0.0005 &  0.0000\\

						\hline

					\end{tabular}
				}

			\end{center}
			\caption{Average and Standard Deviation of Mean Absolute Error as Observations Increase}
			\label{T:Converge}
		\end{table}
	}
	As $n$ increases for a fixed number of observations, the estimates of $\rho$ appear to become more accurate. This counterintuitive result is a consequence of the way that $\rho$ is distributed in this set of examples.  As $n$ increases, the mean and variance of $\rho_i$ decreases and leads to the perception of less error.

	Mean absolute error of the matrix is only a measure of overall estimator performance and cannot give the details of each element of the adjacency matrix.  To explore the manner in which individual elements of the estimate change as the number of observations increases, we use the same latent network with 10 nodes which was used above and generate datasets of increasing size. It can be seen from Figure \ref{F:NodeConverge} that the difference between the parameter and estimates is quickly reduced to very small values and appears to converge to 0.

	\begin{figure}[H]
		\centering
		\includegraphics[width=3in]{NodeConverge-eps-converted-to}
		\caption{Convergence of Individual Node Estimates}
		\label{F:NodeConverge}
	\end{figure}

	\section{DATA ANALYSIS}\label{S:Data}

	\subsection{Introduction}

	In this section, we perform data analysis on three datasets.

	The first dataset records co-sponsorship of legislation in the Senate of the 110$^{th}$ United States Congress.  The rules of the Senate require that each piece of legislation have a unique sponsor; however, other members may co-sponsor the bill \cite{Fowler06a}.  These rules mean that the data conform to the assumption of the Hub Model.

	The second dataset is from the 18$^{th}$ century Chinese novel, \textit{Dream of the Red Chamber}.  The observed groups in this dataset do not necessarily conform to the Hub Model assumption.  However, we will show that even without this assumption being explicitly valid, important information about the relationships between individuals can be estimated.

	The final dataset has been extracted from the USDA plant database.  Unlike the first two datasets, this one does not deal with ``social'' data, but with ``spatial'' data.  Each observation represents a single species of plant along with each North American state or territory in which the plant is observed to grow.  For this analysis, states and territories represent nodes. As with the second dataset, we find that the Hub Models return meaningful information about underlying structure.

	\subsection{Senate of the 110$^{th}$ United States Congress}\label{S:Senate}

	The United States Senate is a chamber in the bicameral legislature of the United States, and together with the U.S. House of Representatives makes up the U.S. Congress.  A key function of both chambers of Congress is to originate legislation.  Each piece of legislation can have only one originating sponsor; however, since the mid-1930s, Senators have had an opportunity to express support for a piece of legislation by signing it as a co-sponsor \cite{Fowler06a}.

	The 110$^{th}$ United States Congress occurred between January 3, 2007 and January 3, 2009.  The Democratic Party controlled a majority in both chambers for the first time since 1995 with a voting share of 50.5 \% of the Senate membership.

	The United States Senate consists of 100 members with each state represented by 2 Senators at any time.  During this session of Congress, there were a total of 102 individuals who served in the Senate.  One original member died and a second resigned to become a lobbyist.  Both members were replaced by appointed state representatives.

	Data for legislative co-sponsorship are available in the Library of Congress Thomas legislative database. This database includes more than 280,000 pieces of legislation proposed in the U.S. House and Senate with over 2.1 million co-sponsorship signatures. Most bills do not pass, and cosponsors need not invest time and resources crafting legislation; so co-sponsorship is a relatively costless way to signal one's position on issues important to constituents and fellow legislators. For the purposes of this study, we include all forms of legislation including all available resolutions, public and private bills, and amendments \cite{Fowler06b}.  During the 110$^{th}$ Congress, the Senate initiated 10,327 pieces of legislation.

	It is a trivial task to apply the KHM to this dataset when we treat the sponsor as known; therefore, we focus on the case where the sponsor is unknown.  That is, we intentionally confound sponsors with co-sponsors so that the only data that we have is $G$. We would like to investigate whether the HM can provide a meaningful estimate of the latent network even when the information of hub nodes is missing. The average difference between edges estimated by KHM and HM is 0.03, which suggests that the HM estimate is very accurate even when we confound the hub nodes.

	In Figure \ref{F:Senate}, we plot the co-occurrence matrix, half-weight index, and the adjacency matrix of the Hub Model using the force directed graph drawing technique of Fruchterman-Reingold.  Each Senator is represented as a node where the color of the node represents the Senator's official political party. Red nodes represent Republicans while blue nodes represent Democrats. We omit the links between nodes in the three sub-figures, since the networks constructed by the co-occurrence matrix and half-weight index are almost complete graphs. The estimate $\hat{A}$ of HM in Figure \ref{F:Senate_A} shows a clearly clustered pattern of Republicans and Democrats. By contrast, the figures of co-occurrence matrix and half-weight index show no pattern.

	\begin{figure}[H]
		\centering
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\textwidth,trim = 20mm 20mm 20mm 20mm, clip]{plot_Senate_O-eps-converted-to}
			\caption{Co-occurrence}
			\label{F:Senate_O}
		\end{subfigure
		~
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\textwidth,trim = 20mm 20mm 20mm 20mm, clip]{plot_Senate_H-eps-converted-to}
			\caption{Half Weight Index }
			\label{F:Senate_H}
		\end{subfigure}
		~
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\textwidth,trim = 20mm 20mm 20mm 20mm, clip]{plot_Senate_A-eps-converted-to}
			\caption{Hub Model }
			\label{F:Senate_A}
		\end{subfigure}
		\caption{Comparison of Estimation Techniques for the 110$^{th}$ Senate}\label{F:Senate}
	\end{figure}

	This implies that the HM provides more meaningful information about the community structure of this dataset than classical measures.  This is further supported by the normalized cut value corresponding to the Senator's official party membership \cite{Shi00} .
	\begin{equation}\label{E:Normalized Cut}
	\frac{\sum\limits_{i\in C_1,j\in C_2} A_{ij}}{\sum\limits_{i,j\in C_1} A_{ij}}+\frac{\sum\limits_{i\in C_1,j\in C_2} A_{ij}}{\sum\limits_{i,j\in C_2} A_{ij}}.
	\end{equation}

	\eqref{E:Normalized Cut} gives the normalized cut value for two communities $C_1$ and $C_2$.  Lower normalized cut values indicate stronger community differentiation.  For the Senate data, Table \ref{T:Normalized Cut} presents the values for each inferred network. The normalized cut value for HM is lower than the co-occurrence matrix and half weight index, which strengthens the visual intuition that the estimates from HM provide better distinction between communities.
	{\linespread{1.0}
		\begin{table}[hbt!]
			\begin{center}
				\begin{tabular}{|c| c|}
					\hline
					& Normalized Cut Value\\
					\hline

					Co-Occurrence Matrix    &  0.837   \\
					\hline

					Half Weight Index &  0.823  \\
					\hline
					Hub Model & 0.757 \\

					\hline
				\end{tabular}
			\end{center}
			\caption{Normalized Cut Ratios of 110$^{th}$ Senate for Different Inference Techniques}
			\label{T:Normalized Cut}
		\end{table}
	}

	\subsection{Dream of the Red Chamber}

	As noted by Kolaczyk \cite{Kolaczyk09}, a significant challenge with parameter estimation in the context of implicit networks is that for a real world dataset there is usually no way to verify the extent to which the estimate matches reality.  That is, there is no so-called ``ground truth'' or ``golden standard'' to evaluate the performance of the estimated results against.  Therefore, to test the performance of Hub Models, it is useful to analyze data about which there is some qualitative knowledge of the relationships between nodes.

	To this end, we construct a dataset of  characters from the 18th century Chinese novel \textit{Dream of the Red Chamber}, also known as \textit{The Story of the Stone}.  Since novels contain a qualitative social structure that is familiar to readers, the results of quantitative analysis can be compared to this standard.  This novel is used here because the relationships between the characters are subtle and complex. In this way, the story presents a challenge to the task of estimating the relationships between the characters.

	Traditional approaches to building datasets from novels require carefully reading the text and identifying dyadic interactions between characters based on criteria established by the researchers, e.g., character $A$ has a conversation with character $B$ \cite{MacCarron13}. This method may construct high quality datasets; however, in order to identify the dyadic interaction, it requires readers who are familiar with the languages used in the novel and who have time to build the datasets.  Since \textit{Dream of the Red Chamber} is written in classical Chinese and the English translation runs over 2,600 pages, directly generating the dataset would be excessively time consuming.

	To overcome the size of the novel, we built the dataset by employing an automatic text mining technique. We define a group as characters who co-occur in the same paragraph. The novel can be text-mined automatically since character names and paragraph markers are easily identified.  Paragraphs with no characters named in them are ignored.

	We analyze the relationships of 29 important characters. The character names presented here are based on the original pinyin pronunciations and the David Hawkes translation \cite{Hawkes74}.  A Chinese version of the novel was used for text-mining.

	This complete novel contains 120 chapters, but we focus on the first 80 chapters because it is commonly believed that the last 40 chapters are written by a different author and may not reflect the original themes of the novel.  The resulting dataset has 1,389 observations of groups containing at least one of the 29 characters.

	Figure \ref{F:DRC} presents the results of the three techniques for estimating social structure discussed throughout this paper.  In this section, we employ an alternative visualization technique for the relationships between individual nodes.  The adjacency matrix, $A$ is represented as an $n \times n$ grid where the $i^{th} \times j^{th}$ cell represents the relationship $A_{ij}$.  The strength of a relationship is represented by the cell's color.  Nodes with weak relationships have light cells while nodes with strong relationship have dark cells.  Cells representing relationships of intermediate strength are shaded along the gray scale.

	This visualization demonstrates another difference in the performance of the techniques.  The co-occurrence matrix estimates all relationships as being very weak and it is difficult to differentiate strong relationships from the absence of a relationship.  The half-weight index presents a much stronger set of relationships but there is evidence of relationships which have been imputed transitively. In general, HM returns a much sparser network where the strength of relationships demonstrates higher contrast.  This tendency towards sparsity is discussed in more detail in Section \ref{S:Discussion}.

	\begin{figure}[H]
		\centering
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\textwidth]{DRC_O_names-eps-converted-to}
			\caption{Co-Occurrence Matrix}
			\label{F:DRC_O_names}
		\end{subfigure}
		~
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\textwidth]{DRC_H_names-eps-converted-to}
			\caption{Half Weight Index}
			\label{F:DRC_H_names}
		\end{subfigure}
		\\
		\begin{subfigure}[b]{0.3\textwidth}
			\includegraphics[width=\textwidth]{DRC_A_names-eps-converted-to}
			\caption{HM Adjacency Matrix}
			\label{F:DRC_A_names}
		\end{subfigure}
		\caption{Comparison of Estimation Techniques for \textit{Dream of the Red Chamber}}\label{F:DRC}
	\end{figure}

	The EM algorithm of HM provides very stable solutions. By selecting multiple starting points, we find that the adjacency matrix (Figure \ref{F:DRC_A_names}) is repeatedly returned as the most likely parameter of the observed data.

	The standard deviation of the parameters of the Hub Model was estimated using the bootstrap technique.  In general, the standard deviation was low. This was particularly true for $\hat{\rho}$ where the maximum standard deviation was 0.0173. Table \ref{T:DRC_Percentiles} presents the standard deviation of $\hat{A}$ at different percentiles.
	{\linespread{1.0}
		\begin{table}[hbt!]
			\begin{center}
				\resizebox{5in}{!}{
					\begin{tabular}{|c| c c c c c c c|}
						\hline
						Percentile & Max & 95 \% & 75 \% & Med & 25 \% & 5 \% & Min \\
						\hline
						StDev &         0.2696 & 0.1025 & 0.0374 & 0.0100 & 0.0000 & 0.0000 &  0.0000 \\
						\hline
					\end{tabular}
				}
			\end{center}
			\caption{Percentiles of Standard Deviation in $\hat{A}$ estimated by HM for \textit{Dream of the Red Chamber}}
			\label{T:DRC_Percentiles}
		\end{table}
	}

	One of the main themes of the \textit{Dream of the Red Chamber} is the love story surrounding the protagonist Jia Baoyu (1st character in Figure \ref{F:DRC_A_names}) and two potential fiancees.  These are the sickly Lin Daiyu (2nd character) and the ``ideal'' Xue Baochai (3rd character).  Although Jia Baoyu shares a special bond with Lin Daiyu and has no significant emotional connection to Xue Baochai, he is ultimately tricked into marrying Xue Baochai.  In Table \ref{T:DRC_Main}, we present the relationships between these two girls and the other characters as estimated by the co-occurrence matrix, half weight index, and HM.

	From the novel, Lin Daiyu is a sensitive and sickly girl who prefers to be alone. By contrast, Xue Baochai is a social and calculating girl. She is extremely good at interpersonal communication especially with the protagonist's mother (Lady Wang) and grandmother (Grandmother Jia). These significantly different personalities are clearly represented by the HM estimator while the other estimators do not identify this difference.

	Xue Baochai generally has much stronger relationships with other characters, except for three: Jia Baoyu (the protagonist), Miaoyu (a nun with a very similar personality to Lin Daiyu) and Zijuan (a maid of Lin Daiyu). The co-occurrence matrix and half weight index fail to identify such a clear pattern.
	{\linespread{1.0}
		\begin{table}[H]
			\begin{center}
				\resizebox{6in}{!}{
					\begin{tabular}{| c | c c | c c | c c |}
						\hline
						& \multicolumn{2}{c|}{Co-Occurrence Matrix ($O$)} & \multicolumn{2}{c|}{Half Weight Index ($H$)} & \multicolumn{2}{c|}{Hub ($\hat{A}$)}\\
						\hline
						& Lin   & Xue      & Lin   & Xue      & Lin & Xue  \\
						& Daiyu & Baochai  & Daiyu & Baochai  & Daiyu & Baochai\\
						\hline
						Jia Baoyu &   0.1728 &	0.1274 &	0.4563 &	0.3587 &	0.3113 &	0.2258 \\
						Lin Daiyu &   1.0000 &	0.1109 &  1.0000 &	0.4866 &	1.0000 &	0.4072 \\
						Xue Baochai & 0.1109 &	1.0000 &	0.4866 &	1.0000 &	0.4072 &	1.0000 \\
						Jia Yuanchun &0.0072 &	0.0050 &	0.0531 &	0.0449 &	0.0156 &	0.0228 \\
						Jia Tanchun  & 0.0439 &	0.0533 &	0.2490 &	0.3482 &	0.0915 &	0.4848 \\
						Shi Xiangyun & 0.0590 &	0.0490 &	0.3273 &	0.3119 &	0.2194 &	0.2365 \\
						Miaoyu       & 0.0072 &	0.0036 &	0.0552 &	0.0337 &	0.0597 &	0\\
						Jia Yingchun & 0.0252 &	0.0274 &	0.1667 &	0.2141 &	0 &	0.2846 \\
						Jia Xichun   & 0.0187 &	0.0202 &	0.1313 &	0.1692 &	0.0102 &	0.2461 \\
						Wang Xifeng  & 0.0497 &	0.0526 &	0.1840 &	0.2131 & 	0.0317 &	0.0697 \\
						Jia Qiaojie  & 0.0022 &	0.0022 &	0.0170 &	0.0208 &	0 &	0.0348 \\
						Li Wan       & 0.0367 &	0.0482 &	0.2086 &	0.3160 &	0.0580 &	0.3384 \\
						Qin Keqing   & 0.0007 &	0.0007 &	0.0052 &	0.0062 &	0 &	0 \\
						\hline
						Grandmother Jia    & 0.0655 &	0.0648 &	0.2725 &	0.2985 &	0.1925 &	0.2820 \\
						\hline
						Jia She      & 0.0065 &	0.0043 &	0.0449 &	0.0357 &	0 &	0 \\
						Jia Zheng    & 0.0122 &	0.0144 &	0.0701 &	0.0952 &	0.0143 &	0.0174 \\
						Jia Lian     & 0.0072 &	0.0036 &	0.0423 &	0.0245 &	0.0002 &	0.0073 \\
						Xiangling    & 0.0180 &	0.0252 &	0.1185 &	0.1961 &	0.0741 &	0.2344 \\
						Ping'er      & 0.0122 &	0.0209 &	0.0668 &	0.1306 &	0.0016 &	0.1643 \\
						Xue Pan      & 0.0043 &	0.0101 &	0.0292 &	0.0809 &	0 &	0  \\
						Granny Liu   & 0.0072 &	0.0050 &	0.0493 &	0.0411 &	0.0101 &	0.0113 \\
						\hline
						Lady Wang    & 0.0490 &	0.0590 &	0.2248 &	0.3037 &	0.0224 &	0.2065 \\
						\hline
						Aunt Xue     & 0.0302 &	0.0396 &	0.1806 &	0.2750 &	0.0479 &	0.1657 \\
						Hua Xiren    & 0.0403 &	0.0389 &	0.1938 &	0.2105 &	0.0283 &	0.1469 \\
						Qingwen      & 0.0166 &	0.0115 &	0.1020 &	0.0829 &	0.0155 &	0.0886 \\
						Yuanyang     & 0.0086 &	0.0101 &	0.0556 &	0.0763 &	0  &	0.0430 \\
						Mingyan      & 0.0007 &	0.0007 &	0.0053 &	0.0064 &	0  &	0  \\
						Zijuan       & 0.0317 &	0.0108 &	0.2184 &	0.0888 &	0.1775 &	0.0376 \\
						Concubine Zhao & 0.0050 &	0.0058 &	0.0361 &	0.0495 &	0 &	0.0338 \\
						\hline
					\end{tabular}
				}
			\end{center}
			\caption{Relationships of Lin Daiyu and Xue Baochai to other characters in \textit{Dream of the Red Chamber}}
			\label{T:DRC_Main}
		\end{table}
	}

	Figure \ref{F:DRC} and Table \ref{T:DRC_Main} show that even when the model assumption is not necessarily valid, important distinctions can be drawn from the Hub Model which are not possible with the other techniques.

	\subsection{North American Flora}\label{S:Flora}

	In the previous examples, we have worked with datasets which are essentially from the social sciences.  However, we believe that Hub Models are useful in other situations where observations are the result of nodes coalescing around a single node or observations are the result of some resource dispersing outward from a single node to multiple nodes.

	As a demonstration of how this kind of data can be used to estimate the relationship between different regions, we use a dataset from the University of California Irvine Machine Learning Repository which had been extracted and encoded from the USDA plants database \cite{Hamalainen08}.  34,781 plant species or geneses are included in the dataset.  For each plant, the dataset indicates which of 68 areas the plant is found in.  These areas include all United States states, Canadian provinces and territories, along with the Virgin Islands, Puerto Rico, Greenland, and St. Pierre and Miquelon (islands off the northeast coast of Canada).  For simplicity, we will refer to these areas as \textit{states}.

	We would expect that contiguous states would tend to have many flora in common while states which are far apart would be less likely to share common flora.  For example, Connecticut and Massachusetts are small states which share a common boarder; therefore, we would expect them to appear together in the regions of many plants.  Conversely, California and Greenland are very far apart and at different latitudes; therefore, we would expect a weak relationship.

	Of course, the map of North America is well known and our objective here is not to compare the Hub Model to spacial modeling.  Instead, we are using the regions of North America as a proxy for a system of distribution.

	To demonstrate the ability of the Hub Model to capture the connections between states, we split them into 13 different regions.  The United States is identified by the 9 divisions of the US Census Bureau.  The Canadian provinces are identified by three regions.  The final region includes islands in the Atlantic ocean which are not included in any other region.

	\begin{figure}[H]
		\centering
		\fbox{
			\begin{subfigure}[b]{0.2\textwidth}
				\includegraphics[width=\textwidth,trim = 85mm 40mm 73mm 40mm, clip]{plot_Plant_O-eps-converted-to}
				\caption{Co-Occurrence}
				\label{F:Plant_O}
			\end{subfigure}
		}
		~
		\fbox{
			\begin{subfigure}[b]{0.2\textwidth}
				\includegraphics[width=\textwidth,trim = 85mm 40mm 73mm 40mm, clip]{plot_Plant_H-eps-converted-to}
				\caption{Half Weight Index}
				\label{F:Plant_H}
			\end{subfigure}
		}
		~
		\fbox{
			\begin{subfigure}[b]{0.2\textwidth}
				\includegraphics[width=\textwidth,trim = 86mm 40mm 73mm 40mm, clip]{plot_Plant_A-eps-converted-to}
				\caption{Hub Model}
				\label{F:Plant_A}
			\end{subfigure}
		}
		\\
		\begin{subfigure}[b]{.7\textwidth}
			\includegraphics[width=\textwidth]{plot_Plant_legend}
			\label{F:Plant_ledgend}
		\end{subfigure}
		\caption{Adjacency Matrix Estimates for North American Flora Data}\label{F:Plant}
	\end{figure}

	The HM graph is striking in how closely states in the same region are grouped.  Southern states are generally on the right side of Figure \ref{F:Plant_A} while northern states are generally to the left.  Eastern states are generally at the bottom of the figure while western states are generally at the top.  Atlantic islands are on the outer edge of the plot.

	In Figures \ref{F:Plant_O} and \ref{F:Plant_H}, there is almost no distinction between the organization of the states.  This suggests that even in situations where the data does not clearly conform to the Hub Model assumption that valuable information about the relationship between nodes can be identified.

	\section{DISCUSSION}\label{S:Discussion}

	\subsection{Self-sparsity}

	In Section \ref{S:Simulation}, we introduced a property of the Hub Model estimators which we refer to as \textit{self-sparsity}.  When $T$ is small relative to $n$, the model tends to produce a sparse adjacency matrix. Rabbat et al.\cite{Rabbat08} observed similar behavior in their research.  This sparsity in $A$ is achieved without any penalty in the log-likelihood, hence the name.

	To begin, observe that the true probability of co-occurrence is related to $\{A,\rho\}$ by the following equation:
	\begin{equation}\label{E:Sparsity}
	\mathbb{P}(v_i \textnormal{ and } v_j \textnormal{ co-occur})=\sum_{k=1}^n \rho_k A_{ki} A_{kj}.
	\end{equation}

	Suppose that there is a pair of nodes, $\{v_i,v_j\}$, for which the probability of co-occurrence is exactly zero. Equation \eqref{E:Sparsity} implies that:
	\begin{align*}
	\rho_k A_{ki} A_{kj} &= 0 \quad \forall k.
	\end{align*}

	Hence, for every $k$, at least one of the following is true: $\rho_k=0$, $A_{ki}=0$, or $A_{kj}=0$.  At a minimum, this requires that there be $n$ elements of the parameters $\{A,\rho\}$ which are exactly equal to zero for every pair of nodes which fails to co-occur.

	Clearly, $O_{ij}=0$ implies $A_{ij}=0$.  However, self-sparsity shows that the absence of co-occurrence contains even more information than just the relationship between two nodes.  Absence of co-occurrence means that no member of the population chooses to simultaneously interact with both nodes.

	However, the question remains as to why self-sparsity occurs in the estimation of $\{A,\rho\}$ when $T$ is small relative to $n$.  We observe that sparsity in $\hat{A}$ or $\hat{\rho}$ is a consequence of the EM-algorithm of HM.

	First, it is easy to check that zero is an absorbing state for $\hat{A}$ and $\hat{\rho}$ in the EM-algorithm.  If at the $m^{th}$ iteration, $\hat{\rho}_i^{(m)}=0$, then by  \eqref{E:ProbKisCenter}  $\mathbb{P}(S_i=1|G^{(t)})=0$ for all $i$, and by \eqref{E:HM rho_hat} $\hat{\rho}_i^{(m+1)}=0$.  Therefore, $\hat{\rho}_i=0$ is an absorbing state.  For a similar reason, $\hat{A}_{ij}=0$ is also an absorbing state.

	Note that in the EM-algorithm we set $\hat{A}_{ij}^{(m)}=0$ when it is below a certain threshold.  But why the estimate approaches zero is not fully understood.  This aspect of the model will be explored in future work.

	\subsection{Identifiability}

	Recall from Section \ref{S:Ident} that when we allow $A$ to be asymmetric, the model is not identifiable.  That is, the following condition is not satisfied:
	\begin{equation}
	\mathbb{P}(G=g|A,\rho)=\mathbb{P}(G=g|A^*,\rho^*) \hspace{2mm} \forall g \implies \{A,\rho\}=\{A^*,\rho^*\}.
	\end{equation}

	In Section \ref{S:Ident}, we showed that symmetry ($A_{ij}=A_{ji}$) is an identifiability condition.  In this subsection we explore identifiability in more detail.

	We give a simple counterexample in Section \ref{S:Ident} to demonstrate that without any constraints on $A$ and $\rho$, the model is not identifiable. Here we will randomly select parameters to explore the issue in general. We randomly generate an asymmetric adjacency matrix $\{A,\rho\}$ with $n=4$ (see Table \ref{T:A12}).  Note that for the asymmetric case, $n=4$ is the smallest population size where the  number of possible groups exceeds the number of parameters to estimate.  The number of observed groups is set high ($T=100,000$) to ensure good performance of the algorithm.  We ran Algorithm 1 100 times, and obtained 100 different estimators with the same or very close likelihoods. Figure \ref{F:A12} gives a scatterplot of the one hundred pairs of $\{\hat{A}_{1,2},\hat{A}_{2,1}\}$ indicated by blue circles.  This plot clearly demonstrates a non-linear relationship between these two values where increases in one are associated with decreases in the other.
	{\linespread{1.0}
		\begin{table}[H]
			\begin{center}
				\begin{tabular}{| l | c | l l l l |}
					\cline{2-6}
					\multicolumn{1}{c}{}  &  \multicolumn{5}{|c|}{$A_{ij}$}\\
					\cline{3-6}
					\multicolumn{1}{c}{}  & \multicolumn{1}{|c}{}    & \multicolumn{4}{|c|}{j}\\
					\hline
					$\rho_i$ &  i & 1 & 2 & 3 & 4 \\
					\hline
					0.5499 & 1 & 1.0000 & 0.7854 & 0.9063 &  0.7957 \\
					0.3269 & 2 & 0.7032 & 1.0000 & 0.8324 &  0.5885 \\
					0.1016 & 3 & 0.9464 & 0.8817 & 1.0000 &  0.9334 \\
					0.0216 & 4 & 0.7452 & 0.8594 & 0.9478 &  1.0000 \\
					\hline
				\end{tabular}
			\end{center}
			\caption{True Adjacency Matrix in Identifiability Example}
			\label{T:A12}
		\end{table}
	}
	\begin{figure}[H]
		\centering
		\includegraphics[width=0.4\textwidth]{A12-eps-converted-to}
		\caption{Nonlinear Relationship of Symmetric Elements for $\hat{A}$}
		\label{F:A12}
	\end{figure}

	We now derive the theoretical relationship between symmetric elements of the adjacency matrix. As in Section \ref{S:Ident}, let $g^x$ and $g^y$ represent the observed groups which contain only node $v_x$ and $v_y$ respectively.  Further, let $g^{xy}$ represent the group that is observed to contain only the pair $v_x$ and $v_y$.

	Under the Hub Models $\mathbb{P}(G=g^x)$ and $\mathbb{P}(G=g^x)$ are given by \eqref{E:Prob_gx} and \eqref{E:Prob_gy}.  In the asymmetric case, the probability of the pair is:
	\begin{equation}\label{E:gxy}
	\mathbb{P}(G=g^{xy})=\rho_x A_{xy} \prod_{j\ne \{x,y\}}(1-A_{xj}) + \rho_y A_{yx} \prod_{j\ne \{x,y\}}(1-A_{yj}).
	\end{equation}

	By simply reordering the terms of Equations \ref{E:Prob_gx} and \ref{E:Prob_gy}, they can be plugged back into Equation \ref{E:gxy} to give:
	\begin{equation}\label{E:gxy2}
	\mathbb{P}(G=g^{xy})=\frac{A_{xy}}{(1-A_{xy})} \mathbb{P}(G=g^x)  +
	\frac{A_{yx}}{(1-A_{yx})} \mathbb{P}(G=g^y)
	\end{equation}

	By some simple algebra, we find the following relationship between the elements of the adjacency matrix:
	\begin{equation}\label{E:SymRule}	A_{xy}=\frac{\mathbb{P}(G=g^{xy})-A_{yx}\big[\mathbb{P}(G=g^{y})+\mathbb{P}(G=g^{xy})\big]}{\big[\mathbb{P}(G=g^{x})+\mathbb{P}(G=g^{xy})\big]-A_{yx}\big[\mathbb{P}(G=g^{x})+\mathbb{P}(G=g^{y})+\mathbb{P}(G=g^{xy})\big]}.
	\end{equation}

	Using \eqref{E:SymRule}, we can calculate the relationship between $A_{1,2}$ and $A_{2,1}$ from the example above.  This is represented in Figure \ref{F:A12} by the solid line.  Clearly the observed solutions are falling along this theoretical curve.

	\section{CONCLUSION}

	To the best of our knowledge, Hub Models introduce an innovative approach to the task of implicit network inference.  By defining a model-based generating mechanism to link the latent network to observed grouped data and applying an EM algorithm, we are able to estimate the network using this model.

	Not only are the estimators easy to calculate in a reasonable amount of time, but the estimators have a practical interpretation.  The parameter $\rho$ measures the probability that a member of a population will form a group.  $A_{ij}$ measures the probability that a member of the population will be included in a group formed by node $v_i$.

	The Hub Models compare favorably against existing techniques.  Since the co-occurrence matrix and half weight index lack a generating mechanism to connect them to the observed grouped data, these measures often cannot detect important features of a network.

	By applying the Hub Model to the 110$^{th}$ United States Congress dataset, the 18$^{th}$ century Chinese novel \textit{Dream of the Red Chamber}, and a dataset of flora in North America, we demonstrate that the HM is able to detect important features in the relationships between nodes in complex situations.  We conclude by providing some initial insights into self-sparsity and the identifiability of HM. To fully understand these properties is an intriguing research topic and will be explored in our future works.

	\par

	\vskip 14pt
	\noindent {\large\bf Acknowledgements}

	This work is partially supported by NSF DMS 1513004.
	\par