\section{Introduction}

In many modern problems, we are interested in studying a network of
entities with multiple attributes rather than a simple univariate
attribute. For example, when an entity represents a person in a social
network, it is widely accepted that the nodal attribute is most
naturally a vector with many personal information including
demographics, interests, and other features, rather than merely a
single attribute, such as a binary vote as assumed in the current
literature of social graph estimation based on Markov random fields
\citep{Banerjee08,kolar08estimating}. In another example, when an
entity represents a gene in a gene regulation network, modern data
acquisition technologies allow researchers to measure the activities
of a single gene in a high-dimensional space, such as an image of the
spatial distribution of the gene expression, or a multi-view snapshot
of the gene activity such as mRNA and protein abundances, rather than
merely a single attribute such as an expression level, which is
assumed in the current literature on gene graph estimation based on
Gaussian graphical models~\citep{peng09partial}. Indeed, it is
somewhat surprising that existing research on graph estimation remains
largely blinded to the analysis of multi-attribute data that are
prevalent and widely studied in the network community. Existing
algorithms and theoretical analysis relies heavily on covariance
selection using graphical lasso, or penalized pseudo-likelihood. They
can not be easily extended to graphs with multi-variate nodal
attributes.

In this paper, we present a study on graph estimation from
multi-attribute data, in an attempt to fill the gap between the
practical needs and existing methodologies from the literature. Under
a Gaussian graphical model, one assumes that a $p$-dimensional random
vector $X \in \RR^p$ follows a multivariate Gaussian distribution with
the mean $\mu$ and covariance matrix $\Sigma$, with each component of
the vector corresponding to a node in the graph.  Based on $n$
independent and identically distributed observations, one can estimate
an undirected graph $G=(V,E)$, where the node set $V$ corresponds to
the $p$ variables, and the edge set $E$ describes the conditional
independence relationships among the variables, that is, variables
$X_a$ and $X_b$ are conditionally independent given all the remaining
variables if $(a,b)\notin E$. Given multi-attribute data, this
approach is clearly invalid, because it naively translates to
estimating one graph per attribute. A subsequent integration of all
such graphs to a summary graph on the entire dataset may lead to
unclear statistical interpretation.

We consider the following new setting for estimating a multi-attribute
graph.  In this setting, we consider a ``stacked'' long random vector
$X=(X^T_1,...,X^T_p)^T$ where $X_1 \in \RR^{k_1}, \ldots, X_p
\in \RR^{k_p}$ are themselves random vectors that jointly follow a
multivariate Gaussian distribution with mean $\mu = (\mu_1, \ldots,
\mu_p)^T$ and covariance matrix $\Sigma^*$, which is partitioned as
\begin{equation}
  \label{eq:vector_model}
\Sigma^* =   \left(
    \begin{array}{cccc}
     \Sigma_{11}^*  & \cdots & \Sigma_{1p}^* \\
     \vdots &   \ddots & \vdots \\
     \Sigma_{p1}^* & \cdots   & \Sigma_{pp}^* \\
    \end{array}
  \right),
\end{equation}
with $\Sigma^*_{ij} = \Cov(X_i, X_j)$.  Without loss of generality, we
assume $\mu=0$.  Let $G = (V,E)$ be a graph with the vertex set $V =
\{1,\ldots, p\}$ and the set of edges $E \subseteq V \times V$ that
encodes the conditional independence relationships among $(X_a)_{a\in
  V}$. That is, each node $a \in V$ of the graph $G$ corresponds to
the random vector $X_a$ and there is no edge between nodes $a$ and $b$
in the graph if and only if $X_a$ is conditionally independent of
$X_b$ given all the vectors corresponding to the remaining nodes,
$X_{\neg{ab}} = \{X_c\ :\ c \in V\bks\{a,b\}\}$. Such a graph is also
known as a Markov network (of Markov graph), which we shall emphasize
in this paper to compare with an alternative graph over $V$ known as
the association network, which is based on pairwise marginal
independence. Conditional independence can be read from the inverse of
the covariance matrix, as the block corresponding to $X_a$ and $X_b$
will be equal to zero. Let $\Dcal_n = \left\{x_i\right\}_{i=1}^n$ be a
sample of $n$ independent and identically distributed vectors drawn
from $\mathrm{N}(0, \Sigma)$. For a vector $x_i$, we denote $x_{i,a}
\in \RR^{k_a}$ the component corresponding to the node $a \in V$. Our
goal is to estimate the structure of the graph $G$ from the sample
$\Dcal_n$. Note that we allow for different nodes to have different
number of attributes, which is useful in many applications, e.g., when
a node represents a gene pathway in a regulatory network.

Using the standard Gaussian graphical model for univariate nodal
observations, one can estimate the Markov graph for each attribute
individually by estimating the sparsity pattern of the precision
matrix $\Omega = \Sigma^{-1}$.  This is also known as the covariance
selection problem \citep{Dempster72}. For high dimensional problems,
\citet{Meinshausen06} propose a parallel Lasso approach for estimating
Gaussian graphical models by solving a collection of sparse regression
problems. This procedure can be viewed as a pseudo-likelihood based
method. In contrast, \citet{Banerjee08}, \citet{Yuan07}, and
\citet{Friedman08} take a penalized likelihood approach to estimate
the sparse precision matrix $\Omega$. To reduce estimation bias,
\citet{Lam09}, \citet{Jalali12}, and \citet{Shen12} developed the
non-concave penalties to penalize the likelihood function. More
recently, \citet{Yuan10} and \citet{Cai11a} proposed the graphical
Dantzig selector and CLIME, which can be solved by linear programming
and are more amenable to theoretical analysis than the penalized
likelihood approach. Under certain regularity conditions, these
methods have proven to be graph estimation consistent
\citep{Ravikumar11,Yuan10,Cai11a} and scalable software packages, such
as $\mathrm{glasso}$ and $\mathrm{huge}$, were developed to implement
these algorithms \citep{Zhao12b}. However, in the case of
multi-attribute data, it is not clear how to combine estimated graphs
to obtain a single Markov graph reflecting the structure of the
underlying complex system. This is especially the case when nodes in
the graph contain different number of attributes.

In a previous work, \citet{katenka2011multi} proposed a method for
estimating association networks from multi-attribute data using
canonical correlation as a dependence measure between two groups of
attributes. However, association networks are known to confound the
direct interactions with indirect ones as they only represent marginal
associations.  In contrast, we develop a method based on partial
canonical correlation, which give rise to a Markov graph that is
better suited for separating direct interactions from indirect
confounders.  Our work is related to the literature on simultaneous
estimation of multiple Gaussian graphical models under a multi-task
setting
\citep{Guo:09,varoquaux2010brain,honorio2010multi,Chiquet2011inferring,danaher2011joint}.
However, the model given in \eqref{eq:vector_model} is different from
models considered in various multi-task settings and the optimization
algorithms developed in the multi-task literature do not extend to
handle the optimization problem given in our setting.

Unlike the standard procedures for estimating the structure of
Gaussian graphical models (e.g., neighborhood selection
\citep{Meinshausen06} or glasso \citep{Friedman08}), which infer the
partial correlations between pairs of multi-attribute nodes, our
proposed method estimates the {partial canonical correlations} between
pairs of nodes, which leads to a graph estimator over multi-attribute
nodes that bears the same probabilistic independence interpretations
as that of the graph from Gaussian graphical model over univariate
nodes. Under this new framework, the contributions of this paper
include: (i) Computationally, an efficient algorithm is provided to
estimate the multi-attribute Markov graphs; (ii) Theoretically, we
provide sufficient conditions which guarantee consistent graph
recovery; and (iii) Empirically, we apply our procedure to uncover
gene regulatory networks from gene and protein profiles, and to
uncover brain connectivity graph from functional magnetic resonance
imaging data.

\section{Methodology}

In this section, we propose to estimate the graph by estimating
non-zero partial canonical correlation between the nodes. This leads
to a penalized maximum likelihood objective, for which we develop an
efficient optimization procedure.

\subsection{Preliminaries}
\label{sec:prelim}

Let $X_a$ and $X_b$ be two multivariate random vectors.  Canonical
correlation is defined between $X_a$ and $X_b$ as
\begin{equation*}
\rho_c(X_a, X_b) = \max_{\ub \in \RR^{k_a}, \vb\in\RR^{k_b}} \Cor(\ub^TX_a, \vb^TX_b).
\end{equation*}
That is, computing canonical correlation between $X_a$ and $X_b$ is
equivalent to maximizing the correlation between two linear
combinations $\ub^TX_a$ and $\vb^TX_b$ with respect to vectors $\ub$
and $\vb$. Canonical correlation can be used to measure association
strength between two nodes with multi-attribute observations. For
example, in \citet{katenka2011multi}, a graph is estimated from
multi-attribute nodal observations by elementwise thresholding the
canonical correlation matrix between nodes, but such a graph estimator
may confound the direct interactions with indirect ones.

In this paper, we exploit the partial canonical correlation to
estimate a graph from multi-attribute nodal observations. A graph
is going to be formed by connecting nodes with non-zero partial
canonical correlation. Let $\hat\Ab = \argmin\ \EE\rbr{\norm{X_a -
  \Ab X_{\neg{ab}}}_2^2}$ and $\hat\Bb = \argmin\ \EE\rbr{\norm{X_b -
  \Bb X_{\neg{ab}}}_2^2}$, then the partial canonical correlation
between $X_a$ and $X_b$ is defined as
\begin{equation}
\label{eq:partial_correlation}
  \begin{aligned}
\rho_c(X_a, X_b; X_{\neg{ab}})
= \max_{\ub\in\RR^{k_a}, \vb\in\RR^{k_b}} \Cor\{\ub^T(X_a
- \hat\Ab X_{\neg{ab}}), \vb^T(X_b - \hat\Bb X_{\neg{ab}})\},
  \end{aligned}
\end{equation}
 that is, the partial canonical correlation between $X_a$ and $X_b$
is equal to the canonical correlation between the residual vectors of
$X_a$ and $X_b$ after the effect of $X_{\neg{ab}}$ is
removed \citep{rao_partial_1969}.

Let $\Omega^*$ denote the precision matrix under the model in
\eqref{eq:vector_model}. Using standard results for the multivariate
Gaussian distribution (see also Equation (7) in
\citet{rao_partial_1969}), a straightforward calculation shows
that\footnote{Calculation given in Appendix~\ref{sec:appendix_eq_3}}
\begin{equation}
  \label{eq:comp_cov}
  \rho_c(X_a, X_b; X_{\neg{ab}}) \neq 0
  \quad
  \text{if and only if}
  \quad
   \max_{\ub\in\RR^{k_a}, \vb\in\RR^{k_b}}\ \ub^T\Omega_{ab}^*\vb \neq 0.
\end{equation}
This implies that estimating whether the partial canonical correlation
is zero or not can be done by estimating whether a block of the
precision matrix is zero or not. Furthermore, under the model in
\eqref{eq:vector_model}, vectors $X_a$ and $X_b$ are conditionally
independent given $X_{\neg{ab}}$ if and only if partial canonical
correlation is zero. A graph built on this type of inter-nodal
relationship is known as a {Markov graph}, as it captures both
local and global Markov properties over all arbitrary subsets of nodes
in the graph, even though the graph is built based on pairwise
conditional independence properties. In
$\S$\ref{sec:estimation}, we use the above observations to design
an algorithm that estimates the non-zero partial canonical correlation
between nodes from data $\Dcal_n$ using the penalized maximum
likelihood estimation of the precision matrix.

Based on the relationship given in   \eqref{eq:comp_cov}, we can
motivate an alternative method for estimating the non-zero partial
canonical correlation. Let $\bar{a} = \{b\ :\ b\in V\bks\{a\}\}$
denote the set of all nodes minus the node $a$. Then
\begin{equation*}
  \label{eq:population_block_regression}
  \EE\rbr{X_a\mid X_{\bar a} = x_{\bar a}} =
  \Sigma_{a,\ia}^*\Sigma_{\ia,\ia}^{*,-1}
  x_{\bar a}.
\end{equation*}
Since $\Omega^{*}_{a,\ia} =
-(\Sigma_{aa}^*-\Sigma_{a,\ia}^*\Sigma_{\ia,\ia}^{*,-1}\Sigma_{\ia,a}^*)^{-1}
\Sigma_{a,\ia}^*\Sigma_{\ia,\ia}^{*,-1}$, we observe that a zero block
$\Omega_{ab}$ can be identified from the regression coefficients when
each component of $X_a$ is regressed on $X_{\bar a}$. We do not build
an estimation procedure around this observation, however, we note that
this relationship shows how one would develop a regression
based analogue of the work presented in \citet{katenka2011multi}.

\subsection{Penalized Log-Likelihood Optimization}
\label{sec:estimation}

Based on the data $\Dcal_n$, we propose to minimize the penalized
negative Gaussian log-likelihood under the model in
\eqref{eq:vector_model},
\begin{equation}
  \label{eq:max_ll_opt}
  \min_{\Omegab \succ \zero}\ \Bigl\{\tr \Sbb \Omegab - \log|\Omegab|
    + \lambda \sum_{a, b} \norm{\Omegab_{ab}}_F\Bigr\}
\end{equation}
where $\Sbb = n^{-1}\sum_{i=1}^n\xb_i\xb_i^T$ is the sample covariance
matrix, $\norm{\Omegab_{ab}}_F)$ denotes the Frobenius norm of
$\Omegab_{ab}$ and $\lambda$ is a user defined parameter that controls
the sparsity of the solution $\hat\Omega$. The Frobenius norm penalty
encourages blocks of the precision matrix to be equal to zero, similar
to the way that the $\ell_2$ penalty is used in the group Lasso
\citep{Yuan06model}.  Here we assume that the same number of samples
is available per attribute. However, the same method can be used in
cases when some samples are obtained on a subset of
attributes. Indeed, we can simply estimate each element of the matrix
$\Sbb$ from available samples, treating non-measured attributes as
missing completely at random \citep[for more details
see][]{kolar12missing}.

The dual problem to \eqref{eq:max_ll_opt} is
\begin{equation}
  \label{eq:max_ll_opt_dual}
    \max_{\Sigmab} \sum_{j \in V} k_j + \log|\Sigmab| \qquad
      \text{subject to } \qquad
      \max_{a,b}\ \norm{\Sbb_{ab} - \Sigmab_{ab}}_F \leq \lambda,
\end{equation}
where $\Sigmab$ is the dual variable to $\Omegab$ and $|\Sigma|$
denotes the determinant of $\Sigma$. Note that the primal problem
gives us an estimate of the precision matrix, while the dual problem
estimates the covariance matrix. The proposed optimization procedure,
described below, will simultaneously estimate the precision matrix and
covariance matrix, without explicitly performing an expensive matrix
inversion.

We propose to optimize the objective function in \eqref{eq:max_ll_opt}
using an inexact block coordinate descent procedure, inspired by
\citet{mazumder11flexible}. The block coordinate descent is an
iterative procedure that operates on a block of rows and columns while
keeping the other rows and columns fixed. We write
\begin{equation*}
  \Omegab = \left(
  \begin{array}{cc}
    \Omegab_{aa} & \Omegab_{a,\ia} \\
    \Omegab_{\ia,a} & \Omegab_{\ia,\ia}\\
  \end{array}
  \right),
\quad
  \Sigmab = \left(
  \begin{array}{cc}
    \Sigmab_{aa} & \Sigmab_{a,\ia} \\
    \Sigmab_{\ia,a} & \Sigmab_{\ia,\ia}\\
  \end{array}
  \right),
\quad
  \Sbb = \left(
  \begin{array}{cc}
    \Sbb_{aa} & \Sbb_{a,\ia} \\
    \Sbb_{\ia,a} & \Sbb_{\ia,\ia}\\
  \end{array}
  \right)
\end{equation*}
and suppose that $(\tilde\Omegab, \tilde\Sigmab)$ are the current
estimates of the precision matrix and covariance matrix. With the
above block partition, we have $\log|\Omegab| =
\log(\Omegab_{\ia,\ia}) + \log(\Omegab_{aa} -
\Omegab_{a,\ia}(\Omegab_{\ia,\ia})^{-1}\Omegab_{\ia,a})$.  In the next
iteration, $\hat \Omegab$ is of the form
\begin{equation*}
  \hat\Omegab = \tilde\Omegab + \left(
  \begin{array}{cc}
    \Deltab_{aa} & \Deltab_{a,\ia} \\
    \Deltab_{\ia,a} & \zero
  \end{array}
\right) =
\left(
  \begin{array}{cc}
    \hat\Omegab_{aa} & \hat\Omegab_{a,\ia} \\
    \hat\Omegab_{\ia,a} & \tilde\Omegab_{\ia,\ia} \\
  \end{array}
\right)
\end{equation*}
and is obtained by minimizing
\begin{equation}
  \label{eq:partial_min}
  \tr\Sbb_{aa}\Omegab_{aa}
\! + \! 2\tr\Sbb_{a,\ia}\Omegab_{\ia,a}
  - \log|\Omegab_{aa} - \Omegab_{a,\ia}(\tilde\Omegab_{\ia,\ia})^{-1}\Omegab_{\ia,a}|
  + \lambda\norm{\Omegab_{aa}}_F + 2\lambda\sum_{b \neq
    a}\norm{\Omegab_{ab}}_F.
\end{equation}
Exact minimization over the variables $\Omegab_{aa}$ and
$\Omegab_{a,\ia}$ at each iteration of the block coordinate descent
procedure can be computationally expensive. Therefore, we propose to
update $\Omegab_{aa}$ and $\Omegab_{a,\ia}$ using one generalized
gradient step update (see \citet{beck09fast}) in each iteration.  Note
that the objective function in \eqref{eq:partial_min} is a sum of a
smooth convex function and a non-smooth convex penalty so that the
gradient descent method cannot be directly applied. Given a step size
$t$, generalized gradient descent optimizes a quadratic approximation
of the objective at the current iterate $\tilde\Omegab$, which results
in the following two updates
\begin{align}
  \label{eq:update_aa}
  \hat\Omegab_{aa} &= \argmin_{\Omegab_{aa}} \Big\{
  \tr(\Sbb_{aa} - \tilde\Sigmab_{aa})\Omegab_{aa} +
  \frac{1}{2t}\norm{\Omegab_{aa} - \tilde\Omegab_{aa}}_F^2 +
  \lambda \norm{\Omegab_{aa}}_F\Bigr\},
  \quad\text{ and} \\
  \label{eq:update_ab}
  \hat\Omegab_{ab} &= \argmin_{\Omegab_{ab}}
  \Bigl\{ \tr(\Sbb_{ab} - \tilde\Sigmab_{ab})\Omegab_{ba} +
  \frac{1}{2t}\norm{\Omegab_{ab} - \tilde\Omegab_{ab}}_F^2 +
  \lambda \norm{\Omegab_{ab}}_F\Bigr\},
  \quad \forall b \in {\bar a}.
\end{align}
If the resulting estimator $\hat\Omegab$ is not positive definite or
the update does not decrease the objective, we halve the step size $t$
and find a new update.  Once the update of the precision matrix
$\hat\Omegab$ is obtained, we update the covariance matrix
$\hat\Sigmab$.  Updates to the precision and covariance matrices can
be found efficiently, without performing expensive matrix inversion,
as we show in Appendix~\ref{appendix:a}
(see~\eqref{eq:update_aa_final}--\eqref{eq:update_sigma}).  Combining
all three steps we get the following algorithm:
\begin{enumerate}
\item
Set the initial estimator $\tilde\Omegab = \diag(\Sbb)$ and
$\tilde\Sigmab = \tilde\Omegab^{-1}$. Set the step size $t=1$.
\item
For each $a \in V$ perform the following:
\begin{tabbing}
\qquad Update $\hat\Omegab$ using
\eqref{eq:update_aa_final} and \eqref{eq:update_ab_final}.\\
\qquad If $\hat\Omegab$ is not positive definite,
  set $t \leftarrow t/2$ and repeat the update.\\
\qquad Update $\hat\Sigmab$ using \eqref{eq:update_sigma}.
\end{tabbing}
\item
  Repeat Step 2 until the duality gap
  \[
 \Bigl| \tr (\Sbb \hat \Omegab) - \log|\hat\Omegab|
  + \lambda \sum_{a, b} \norm{\hat\Omegab_{ab}}_F -
  \sum_{j \in V} k_j - \log|\Sigmab|\Bigr| \leq \epsilon,
  \]
  where $\epsilon$ is a prefixed  precision parameter (for example, $\epsilon=10^{-3}$).
\end{enumerate}
Finally, we form a graph $\hat G = (V, \hat E)$ by connecting nodes
with $\norm{\hat\Omegab_{ab}}_F \neq 0$.

Computational complexity of the procedure is given in
Appendix~\ref{appendix:b}.  Convergence of the above described
procedure to the unique minimum of the objective function in
\eqref{eq:max_ll_opt} does not follow from the standard results on the
block coordinate descent algorithm \citep{tseng01convergence} for two
reasons. First, the minimization problem in \eqref{eq:partial_min} is
not solved exactly at each iteration, since we only update
$\Omegab_{aa}$ and $\Omegab_{a,\ia}$ using one generalized gradient step
update in each iteration. Second, the blocks of variables, over which
the optimization is done at each iteration, are not completely
separable between iterations due to the symmetry of the problem. The
proof of the following convergence result is given in
Appendix~\ref{appendix:c}.

\begin{lemma}
  \label{lem:convergence}
  For every value of $\lambda > 0$, the above described algorithm
  produces a sequence of estimates $\cbr{\tilde\Omegab^{(t)}}_{t \geq
    1}$ of the precision matrix that monotonically decrease the
  objective values given in \eqref{eq:max_ll_opt}. Every element of
  this sequence is positive definite and the sequence converges to the
  unique minimizer $\hat \Omegab$ of \eqref{eq:max_ll_opt}.
\end{lemma}

\subsection{Efficient Identification of Connected Components}

When the target graph $\hat G$ is composed of smaller, disconnected
components, the solution to the problem in \eqref{eq:max_ll_opt} is
block diagonal (possibly after permuting the node indices) and can be
obtained by solving smaller optimization problems. That is, the
minimizer $\hat\Omegab$ can be obtained by solving
\eqref{eq:max_ll_opt} for each connected component independently,
resulting in massive computational gains. We give necessary and
sufficient condition for the solution $\hat\Omegab$ of
\eqref{eq:max_ll_opt} to be block-diagonal, which can be easily
checked by inspecting the empirical covariance matrix $\Sbb$.

Our first result follows immediately from the Karush-Kuhn-Tucker
conditions for the optimization problem \eqref{eq:max_ll_opt} and
states that if $\hat\Omegab$ is block-diagonal, then it can be
obtained by solving a sequence of smaller optimization problems.
\begin{lemma}
  If the solution to \eqref{eq:max_ll_opt} takes the form
  $\hat\Omega={\rm diag}(\hat\Omega_1,
  \hat\Omega_2,\ldots,\hat\Omega_l)$, that is, $\hat\Omega$ is a block
  diagonal matrix with the diagonal blocks $\hat\Omega_1, \ldots,
  \hat\Omega_l$, then it can be obtained by solving
  \begin{equation*}
  \min_{\Omegab_{l'} \succ \zero}\ \Bigl\{ \tr \Sbb_{l'} \Omegab_{l'} - \log|\Omegab_{l'}|
    + \lambda \sum_{a, b} \norm{\Omegab_{ab}}_F\Bigr\}
  \end{equation*}
  separately for each $l'=1,\ldots,l$, where $\Sbb_{l'}$ are
  submatrices of $\Sbb$ corresponding to $\Omegab_{l'}$.
\end{lemma}
Next, we describe how to identify diagonal blocks of
$\hat\Omegab$. Let $\Pcal = \{P_1, P_2, \ldots, P_l\}$ be a partition
of the set $V$ and assume that the nodes of the graph are ordered in
a way that if $a \in P_j$, $b \in P_{j'}$, $j < j'$, then $a < b$.
The following lemma states that the blocks of $\hat\Omega$ can be
obtained from the blocks of a thresholded sample covariance matrix.
\begin{lemma}
  \label{lem:nec_suff_block_diag}
  A necessary and sufficient conditions for $\hat\Omegab$ to be block
  diagonal with blocks $P_1, P_2, \ldots, P_l$ is that
  $\norm{\Sbb_{ab}}_F \leq \lambda$ for all $a \in P_j$, $b \in
  P_{j'}$, $j \neq j'$.
\end{lemma}

Blocks $P_1, P_2, \ldots, P_l$ can be identified by forming a $p\times
p$ matrix $\Qb$ with elements $q_{ab} = \ind\{\norm{\Sbb_{ab}}_F >
\lambda\}$ and computing connected components of the graph with
adjacency matrix $\Qb$. The lemma states also that given two penalty
parameters $\lambda_1 < \lambda_2$, the set of unconnected nodes with
penalty parameter $\lambda_1$ is a subset of unconnected nodes with
penalty parameter $\lambda_2$. The simple check above allows us to
estimate graphs on datasets with large number of nodes, if we are
interested in graphs with small number of edges. However, this is
often the case when the graphs are used for exploration and
interpretation of complex systems.
Lemma~\ref{lem:nec_suff_block_diag} is related to existing results
established for speeding-up computation when learning single and
multiple Gaussian graphical models \citep{witten2011new,
  mazumder2011exact, danaher2011joint}. Each condition is different,
since the methods optimize different objective functions.

\section{Consistent Graph Identification}
\label{sec:theory}

In this section, we provide theoretical analysis of the estimator
described in $\S$\ref{sec:estimation}. In particular, we provide
sufficient conditions for consistent graph recovery.  For simplicity
of presentation, we assume that $k_a = k$, for all $a \in V$, that is,
we assume that the same number of attributes is observed for each
node. For each $a=1,\ldots,kp$, we assume that
$(\sigma_{aa}^*)^{-1/2}X_a$ is sub-Gaussian with parameter $\gamma$,
where $\sigma_{aa}^*$ is the $a$th diagonal element of
$\Sigmab^*$. Recall that $Z$ is a sub-Gaussian random variable if
there exists a constant $\sigma \in (0, \infty)$ such that
\begin{equation*}
  \EE\rbr{\exp(tZ)} \leq \exp(\sigma^2t^2),\ \text{for all
  } t \in \RR.
\end{equation*}

Our assumptions involve the Hessian of the function $f(\Ab) =
\tr\Sbb\Ab - \log|\Ab|$ evaluated at the true $\Omegab^*$, $\Hcal =
\Hcal(\Omega^{*}) = (\Omega^{*})^{-1}\otimes(\Omega^{*})^{-1} \in
\RR^{(pk)^2\times(pk)^2}$, with $\otimes$ denoting the Kronecker
product, and the true covariance matrix $\Sigmab^*$. The Hessian and
the covariance matrix can be thought of as block matrices with blocks
of size $k^2 \times k^2$ and $k \times k$, respectively. We will make
use of the operator $\Ccal(\cdot)$ that operates on these block
matrices and outputs a smaller matrix with elements that equal to the
Frobenius norm of the original blocks. For example, $\Ccal(\Sigmab^*)
\in \RR^{p\times p}$ with elements $\Ccal(\Sigmab^*)_{ab} =
\norm{\Sigmab^*_{ab}}_F$. Let $\Tcal = \{(a,b): \norm{\Omegab_{ab}}_F
\neq 0\}$ and $\Ncal = \{(a,b): \norm{\Omegab_{ab}}_F = 0\}$. With
this notation introduced, we assume that the following
{irrepresentable} condition holds. There exists a constant $\alpha \in
[0, 1)$ such that
\begin{equation}
\label{eq:assum:irrepresentable}
\opnorm{
\Ccal\left(
  \Hcal_{\Ncal\Tcal}(\Hcal_{\Tcal\Tcal})^{-1}
\right)}{\infty} \leq 1 - \alpha,
\end{equation}
where $\opnorm{A}{\infty} = \max_i \sum_j|A_{ij}|$.
We will also need the following quantities to specify the results $
\kappa_{\Sigmab^*} = \opnorm{\Ccal(\Sigmab^*)}{\infty}$ and
$\kappa_{\Hcal} =
\opnorm{\Ccal(\Hcal_{\Tcal\Tcal}^{-1})}{\infty}$. These conditions
extend the conditions specified in \citet{Ravikumar11} needed for
estimating graphs from single attribute observations.

We have the following result that provides sufficient conditions for
the exact recovery of the graph.
\begin{proposition}
\label{prop:sparsistency}
Let $\tau > 2$. We set
the penalty parameter $\lambda$ in \eqref{eq:max_ll_opt} as
\[
\lambda = 8k\alpha^{-1}
\rbr{128(1+4\gamma^2)^2(\max_a(\sigma_{aa}^*)^2)n^{-1}(2\log(2k) +
  \tau\log(p))}^{1/2}.
\]
If $
  n > C_1s^2k^2(1+8\alpha^{-1})^2(\tau\log p + \log 4 + 2\log k),$
where $s$ is the maximal degree of nodes in $G$,
$
C_1 =
(48\sqrt{2}(1+4\gamma^2)(\max_a \sigma_{aa}^*)
\max(\kappa_{\Sigmab^*}\kappa_{\Hcal},
\kappa_{\Sigmab^*}^3\kappa_{\Hcal}^2))^2
$
 and
\[
\min_{(a,b)\in\Tcal, a\neq b} \norm{\Omegab_{ab}}_F >
16\sqrt{2}(1+4\gamma^2)(\max_a \sigma_{aa}^*)
          (1+8\alpha^{-1})\kappa_\Hcal k
          \rbr{\frac{\tau\log p + \log 4 + 2\log k}{n}}^{1/2},
\]
then $ \PP\rbr{\hat G = G} \geq 1 - p^{2-\tau}$.
\end{proposition}

The proof of Proposition~\ref{prop:sparsistency} is given in
Appendix~\ref{appendix:c}. We extend the proof of \citet{Ravikumar11}
to accommodate the Frobenius norm penalty on blocks of the precision
matrix. This proposition specifies the sufficient sample size and a
lower bound on the Frobenius norm of the off-diagonal blocks needed
for recovery of the unknown graph.  Under these conditions and
correctly specified tuning parameter $\lambda$, the solution to the
optimization problem in \eqref{eq:max_ll_opt} correctly recovers the
graph with high probability. In practice, one needs to choose the
tuning parameter in a data dependent way. For example, using the
Bayesian information criterion.  Even though our theoretical analysis
obtains the same rate of convergence as that of \citet{Ravikumar11},
our method has a significantly improved finite-sample performance
(More details will be provided in $\S$\ref{sec:simulation}.). It
remains an open question whether the sample size requirement can be
improved as in the case of group Lasso \citep[see, for
example,][]{lounici10oracle}.  The analysis of \cite{lounici10oracle}
relies heavily on the special structure of the least squares
regression. Hence, their method does not carry over to the more
complicated objective function as in \eqref{eq:max_ll_opt}.

\section{Interpreting Edges}
\label{sec:interpretation}

We propose a post-processing step that will allow us to quantify the
strength of links identified by the method proposed in
$\S$\ref{sec:estimation}, as well as identify important attributes that
contribute to the existence of links.

For any two nodes $a$ and $b$ for which $\Omegab_{ab}\neq0$, we define
$\Ncal(a,b) = \{ c \in V \bks \{a,b\}\ :\ \Omegab_{ac} \neq 0 \text{
  or } \Omegab_{bc} \neq 0 \}$, which is the Markov blanket for the set
of nodes $\{ \Xb_a, \Xb_b \}$. Note that the conditional distribution
of $(\Xb_a^T, \Xb_b^T)^T$ given $\Xb_{\neg{ab}}$ is equal to the
conditional distribution of $(\Xb_a^T, \Xb_b^T)^T$ given
$\Xb_{\Ncal(a,b)}$.  Now,
\begin{equation*}
  \begin{aligned}
\rho_c(\Xb_a, \Xb_b; \Xb_{\neg{ab}})
&= \rho_c(\Xb_a, \Xb_b; \Xb_{\Ncal(a,b)})\\
&= \max_{\wb_a\in\RR^{k_a}, \wb_b\in\RR^{k_b}} \Cor(\ub^T(\Xb_a
- \tilde\Ab\Xb_{\Ncal(a,b)}), \vb^T(\Xb_b - \tilde\Bb\Xb_{\Ncal(a,b)})),
  \end{aligned}
\end{equation*}
where $\tilde\Ab = \argmin\ \EE\rbr{\norm{\Xb_a -
  \Ab\Xb_{\Ncal(a,b)}}_2^2}$ and $\tilde\Bb = \argmin\ \EE\rbr{\norm{\Xb_b
  - \Bb\Xb_{\Ncal(a,b)}}_2^2}$. Let $\bar\Sigmab(a,b) = \Var(\Xb_a,
\Xb_b\mid\Xb_{\Ncal(a,b)})$. Now we can express the partial canonical
correlation as
\[
\rho_c(\Xb_a, \Xb_b; \Xb_{\Ncal(a,b)}) =
\max_{\wb_a\in\RR^{k_a}, \wb_b\in\RR^{k_a} }\
\frac{\wb_a^T\bar\Sigma_{ab}\wb_b}
{\rbr{\wb_a^T\bar\Sigma_{aa}\wb_a}^{1/2}\rbr{\wb_b^T\bar\Sigma_{bb}\wb_b}^{1/2}}
\]
where
\[
\bar\Sigmab(a,b) = \left(
\begin{array}{cc}
\bar\Sigmab_{aa} & \bar\Sigmab_{ab} \\
\bar\Sigmab_{ba} & \bar\Sigmab_{bb}
\end{array}
\right).
\]
The weight vectors $\wb_a$ and $\wb_b$ can be easily found by solving
the system of eigenvalue equations
\begin{equation}
  \label{eq:eigenvalue_system}
\left\{
\begin{array}{c}
  \bar\Sigmab_{aa}^{-1}\bar\Sigmab_{ab}\bar\Sigmab_{bb}^{-1}\bar\Sigmab_{ba}\wb_a
  = \phi^2\wb_a \\
  \bar\Sigmab_{bb}^{-1}\bar\Sigmab_{ba}\bar\Sigmab_{aa}^{-1}\bar\Sigmab_{ab}\wb_b
  = \phi^2\wb_b
\end{array}
\right.
\end{equation}
with $\wb_a$ and $\wb_b$ being the vectors that correspond to the
maximum eigenvalue $\phi^2$. Furthermore, we have $\rho_c(\Xb_a,
\Xb_b; \Xb_{\Ncal(a,b)}) = \phi$. Following \citet{katenka2011multi},
the weights $\wb_a$, $\wb_b$ can be used to access the relative
contribution of each attribute to the edge between the nodes $a$ and
$b$. In particular, the weight $(w_{a,i})^2$ characterizes the
relative contribution of the $i$th attribute of node $a$ to
$\rho_c(\Xb_a, \Xb_b; \Xb_{\Ncal(a,b)})$.

Given an estimate $\hat\Ncal(a, b) = \{ c \in V \bks \{a,b\}\ :\
\hat\Omegab_{ac} \neq 0 \text{ or } \hat\Omegab_{bc} \neq 0 \}$ of the
Markov blanket $\Ncal(a, b)$, we form the residual vectors
\[
\rb_{i, a} = \xb_{i,a} - \check{\Ab}\xb_{i, \hat \Ncal(a,b)},
\qquad
\rb_{i, b} = \xb_{i,b} - \check{\Bb}\xb_{i, \hat \Ncal(a,b)},
\]
where $\check{\Ab}$ and $\check{\Bb}$ are the least square estimators
of ${\tilde{\Ab}}$ and ${\tilde{\Bb}}$. Given the residuals, we form
$\check\Sigmab(a,b)$, the empirical version of the matrix $\bar
\Sigmab(a,b)$, by setting
\[
{\check \Sigma}_{aa} = {\Cor}\rbr{\{\rb_{i,a}\}_{i \in [n]}},\quad
{\check \Sigma}_{bb} = {\Cor}\rbr{\{\rb_{i,b}\}_{i \in [n]}},\quad
{\check \Sigma}_{ab} = {\Cor}\rbr{\{\rb_{i,a} \}_{i \in [n]}, \{\rb_{i,a} \}_{i \in [n]}}.
\]
Now, solving the eigenvalue system in \eqref{eq:eigenvalue_system}
will give us estimates of the vectors $\wb_a$, $\wb_b$ and the partial
canonical correlation.

Note that we have described a way to interpret the elements of the
off-diagonal blocks in the estimated precision matrix. The elements of
the diagonal blocks, which correspond to coefficients between
attributes of the same node, can still be interpreted by their
relationship to the partial correlation coefficients.

\section{Simulation Studies}
\label{sec:simulation}

In this section, we perform a set of simulation studies to illustrate
finite sample performance of our method. We demonstrate that the
scalings of $(n,p,s)$ predicted by the theory are sharp. Furthermore,
we compare against three other methods: 1) a method that uses
the ${\rm glasso}$ first to estimate one graph over each of the $k$
individual attributes and then creates an edge in the resulting graph
if an edge appears in at least one of the single attribute graphs, 2)
the method of \cite{Guo:09} and 3) the method of
\cite{danaher2011joint}.  We have also tried applying the ${\rm
  glasso}$ to estimate the precision matrix for the model in
\eqref{eq:vector_model} and then post-processing it, so that an edge
appears in the resulting graph if the corresponding block of the
estimated precision matrix is non-zero. The result of this method is
worse compared to the first baseline, so we do not report it here.

All the methods above require setting one or two tuning parameters
that control the sparsity of the estimated graph. We select these
tuning parameters by minimizing the Bayesian information criterion,
which balances the goodness of fit of the model and its complexity,
over a grid of parameter values. For our multi-attribute method, the
Bayesian information criterion takes the following form
\[
{\rm BIC}(\lambda) = \tr(\Sbb\hat\Omegab) - \log|\hat\Omegab|
+ \sum_{a < b} \ind\{\hat\Omegab_{ab} \neq \zero\}k_ak_b\log(n).
\]
Other methods for selecting tuning parameters are possible, like
minimization of cross-validation or Akaike information criterion.
However, these methods tend to select models that are too dense.

Theoretical results given in $\S$\ref{sec:theory} characterize the
sample size needed for consistent recovery of the underlying graph.
In particular, Proposition~\ref{prop:sparsistency} suggests that we
need $n = \theta s^2k^2\log(pk)$ samples to estimate the graph
structure consistently, for some $\theta > 0$. Therefore, if we plot
the hamming distance between the true and recovered graph
against $\theta$, we expect the curves to reach zero distance for
different problem sizes at a same point. We verify this on randomly
generated chain and nearest-neighbors graphs.

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/glassoNormal}
    \caption{glasso procedure}
    \label{fig:normal:glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/multiTaskNormal}
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:normal:mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/guoNormal}
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:normal:guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.95\textwidth}
    \centering
    \includegraphics[width=\textwidth]{plots/multiAttNormal}
    \caption{Multi-attribute procedure}
    \label{fig:normal:ma}
  \end{subfigure}%
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent
    runs. Off-diagonal blocks are full matrices.}
  \label{fig:simulation:normal}
\end{figure}

We generate data as follows. A random graph with $p$ nodes is created
by first partitioning nodes into $p/20$ connected components, each
with $20$ nodes, and then forming a random graph over these $20$
nodes. A chain graph is formed by permuting the nodes and connecting
them in succession, while a nearest-neighbor graph is constructed
following the procedure outlined in \citet{li06gradient}. That is, for
each node, we draw a point uniformly at random on a unit square and
compute the pairwise distances between nodes. Each node is then
connected to $s=4$ closest neighbors.  Since some of nodes will have
more than $4$ adjacent edges, we  randomly remove edges from nodes that
have degree larger than $4$ until the maximum degree of a node in a
network is $4$. Once the graph is created, we construct a
precision matrix, with non-zero blocks corresponding to edges in the
graph. Elements of diagonal blocks are set as $0.5^{|a-b|}$, $0\leq
a,b\leq k$, while off-diagonal blocks have elements with the same
value, $0.2$ for chain graphs and $0.3/k$ for nearest-neighbor
networks. Finally, we add $\rho I$ to the precision matrix, so that
its minimum eigenvalue is equal to $0.5$.  Note that $s=2$ for the
chain graph and $s=4$ for the nearest-neighbor graph. Simulation
results are averaged over 100 replicates.

Figure~\ref{fig:simulation:normal} shows simulation results.  Each row
in the figure reports results for one method, while each column in the
figure represents a different simulation setting. For the first two
columns, we set $k=3$ and vary the total number of nodes in the
graph. The third simulation setting sets the total number of nodes
$p=20$ and changes the number of attributes $k$. In the case of the
chain graph, we observe that for small sample sizes the method of
\cite{danaher2011joint} outperforms all the other methods. We note
that the multi-attribute method is estimating many more parameters,
which require large sample size in order to achieve high
accuracy. However, as the sample size increases, we observe that
multi-attribute method starts to outperform the other methods. In
particular, for the sample size indexed by $\theta = 13$ all the graph
are correctly recovered, while other methods fail to recover the
graph consistently at the same sample size. In the case of
nearest-neighbor graph, none of the methods recover the graph well for
small sample sizes. However, for moderate sample sizes,
multi-attribute method outperforms the other methods.  Furthermore,
as the sample size increases none of the other methods recover the
graph exactly. This suggests that the conditions for consistent
graph recovery may be weaker in the multi-attribute setting.

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoEasy2}
    \end{subfigure}%
    \caption{glasso procedure}
    \label{fig:easy:glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskEasy2}
    \end{subfigure}%
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:easy:mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoEasy2}
    \end{subfigure}%
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:easy:guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttEasy1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttEasy2}
    \end{subfigure}%
    \caption{Multi-attribute procedure}
    \label{fig:easy:ma}
  \end{subfigure
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent
    runs. Blocks $\Omegab_{ab}$ of the precision matrix $\Omegab$ are
    diagonal matrices. }
  \label{fig:simulation:easy}
\end{figure}

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoHard2}
    \end{subfigure}%
    \caption{glasso procedure}
    \label{fig:hard:glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskHard2}
    \end{subfigure}%
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:hard:mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoHard2}
    \end{subfigure}%
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:hard:guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttHard1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttHard2}
    \end{subfigure}%
    \caption{Multi-attribute procedure}
    \label{fig:hard:ma}
  \end{subfigure}%
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent runs.
    Off-diagonal blocks $\Omegab_{ab}$ of the precision matrix
    $\Omegab$ have zeros as diagonal elements.}
  \label{fig:simulation:hard}
\end{figure}

\begin{figure}[p]
  \centering

  \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/glassoRS2}
    \end{subfigure}%
    \caption{glasso procedure}
    \label{fig:rs::glasso}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiTaskRS2}
    \end{subfigure}%
    \caption{Procedure of \cite{danaher2011joint}}
    \label{fig:rs::mt}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/guoRS2}
    \end{subfigure}%
    \caption{Procedure of \cite{Guo:09}}
    \label{fig:rs::guo}
  \end{subfigure}%

 \begin{subfigure}[b]{0.9\textwidth}
    \centering
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttRS1}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.4\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/multiAttRS2}
    \end{subfigure}%
    \caption{Multi-attribute procedure}
    \label{fig:rs::ma}
  \end{subfigure}%
  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent runs.
    Off-diagonal blocks $\Omegab_{ab}$ of the precision matrix
    $\Omegab$ have elements uniformly sampled from $[-0.3, -0.1]
    \bigcup [0.1,0.3]$.  }
  \label{fig:simulation:RS}
\end{figure}

\subsection{Alternative Structure of Off-diagonal Blocks}

In this section, we investigate performance of different estimation
procedures under different assumptions on the elements of the
off-diagonal blocks of the precision matrix.

First, we investigate a situation where the multi-attribute method does
not perform as well as the methods that estimate multiple graphical
models. One such situation arises when different attributes are
conditionally independent. To simulate this situation, we use the data
generating approach as before, however, we make each block
$\Omegab_{ab}$ of the precision matrix $\Omega$ a diagonal
matrix. Figure \ref{fig:simulation:easy} summarizes results of the
simulation. We see that the methods of \cite{danaher2011joint} and
\cite{Guo:09} perform better, since they are estimating much fewer
parameters than the multi-attribute method. ${\rm glasso}$ does not
utilize any structural information underlying the estimation problem
and requires larger sample size to correctly estimate the graph than
other methods.

A completely different situation arises when the edges between nodes
can be inferred only based on inter-attribute data, that is, when a
graph based on any individual attribute is empty. To generate data
under this situation, we follow the procedure as before, but with the
diagonal elements of the off-diagonal blocks $\Omegab_{ab}$ set to
zero. Figure~\ref{fig:simulation:hard} summarizes results of the
simulation. In this setting, we clearly see the advantage of the
multi-attribute method, compared to other three methods. Furthermore,
we can see that ${\rm glasso}$ does better than multi-graph methods of
\cite{danaher2011joint} and \cite{Guo:09}. The reason is that ${\rm
  glasso}$ can identify edges based on inter-attribute relationships
among nodes, while multi-graph methods rely only on intra-attribute
relationships. This simulation illustrates an extreme scenario where
inter-attribute relationships are important for identifying edges.

So far, off-diagonal blocks of the precision matrix were constructed
to have constant values. Now, we use the same data generating
procedure, but generate off-diagonal blocks of a precision matrix in a
different way. Each element of the off-diagonal block $\Omega_{ab}$ is
generated independently and uniformly from the set $[-0.3, -0.1]
\bigcup [0.1,0.3]$. The results of the simulation are given in
Figure~\ref{fig:simulation:RS}. Again, qualitatively, the results are
similar to those given in Figure~\ref{fig:simulation:normal}, except
that in this setting more samples are needed to recover the graph
correctly.

\subsection{Different Number of Samples per Attribute}

In this section, we show how to deal with a case when different number
of samples is available per attribute. As noted in
$\S$\ref{sec:estimation}, we can treat non-measured attributes as
missing completely at random \citep[see][for more
details]{kolar12missing}.

Let $R = (r_{il})_{i\in\{1,\ldots,n\},l\in\{1,\ldots,pk\}}\in \RR^{n \times pk}$ be an
indicator matrix, which denotes for each sample point $x_i$ the components
that are observed. Then the sample covariance matrix $S =
(\sigma_{lk}) \in \RR^{pk \times pk}$ is estimated as
$\sigma_{lk} = \rbr{\sum_{i=1}^n r_{i,l}r_{i,k}}^{-1}
{\sum_{i=1}^n r_{i,l}r_{i,k}x_{i,l}x_{i,k}}.$
This estimate is plugged into the objective in \eqref{eq:max_ll_opt}.

\begin{figure}[b]
  \centering

  \includegraphics[width=\textwidth]{plots/biometrikaDiffSample}

  \caption{Average hamming distance plotted against the rescaled
    sample size. Results are averaged over 100 independent
    runs. Additional samples available for the first attribute. }
  \label{fig:simulation:diff_sample}
\end{figure}

We generate a chain graph with $p = 60$ nodes, construct a precision
matrix associated with the graph and $k=3$ attributes, and generate $n
= \theta s^2k^2\log(pk)$ samples, $\theta > 0$. Next, we generate
additional $10\%$, $30\%$ and $50\%$ samples from the same model, but
record only the values for the first attribute. Results of the
simulation are given in Figure~\ref{fig:simulation:diff_sample}.
Qualitatively, the results are similar to those presented in
Figure~\ref{fig:simulation:normal}.

\section{Illustrative Applications to Real Data}

In this section, we illustrate how to apply our method to data arising
in studies of biological regulatory networks and Alzheimer's disease.

\subsection{Analysis of a Gene/Protein Regulatory Network}
\label{sec:real_data:sup}

We provide illustrative, exploratory analysis of data from the
well-known NCI-60 database, which contains different molecular
profiles on a panel of 60 diverse human cancer cell
lines.
Data set consists of protein profiles (normalized reverse-phase lysate
arrays for 92 antibodies) and gene profiles (normalized RNA microarray
intensities from Human Genome U95 Affymetrix chip-set for $> 9000$
genes). We focus our analysis on a subset of 91 genes/proteins for
which both types of profiles are available. These profiles are
available across the same set of $60$ cancer cells. More detailed
description of the data set can be found in \citet{katenka2011multi}.

We inferred three types of networks: a network based on protein
measurements alone, a network based on gene expression profiles and a
single gene/protein network. For protein and gene networks we use the
\texttt{glasso}, while for the gene/protein network, we use our
procedure outlined in $\S$\ref{sec:estimation}. We use the stability
selection~\citep{stability:10} procedure to estimate stable
networks. In particular, we first select the penalty parameter
$\lambda$ using cross-validation, which over-selects the number of
edges in a network. Next, we use the selected $\lambda$ to estimate
100 networks based on random subsamples containing 80\% of the
data-points. Final network is composed of stable edges that appear in
at least 95 of the estimated networks. Table~\ref{tab:net_stats}
provides a few summary statistics for the estimated
networks. Furthermore, protein and gene/protein networks share $96$
edges, while gene and gene/protein networks share $104$ edges. Gene
and protein network share only $17$ edges. Finally, $66$ edges are
unique to gene/protein network. Figure~\ref{fig:node_degree} shows
node degree distributions for the three networks.  We observe that the
estimated networks are much sparser than the association networks in
\citet{katenka2011multi}, as expected due to marginal correlations
between a number of nodes.  The differences in networks require a
closer biological inspection by a domain scientist.

\begin{table}[b]

 \caption{Summary statistics for protein, gene, and gene/protein
   networks ($p = 91$).}

 \begin{tabular}{lrrr}
 ~ & { protein network} & {gene network} & {gene/protein network} \\

Number of edges & 122&   214 &   249 \\
Density & 0.03 & 0.05 & 0.06 \\
Largest connected component & 62 & 89 & 82 \\
Avg Node Degree  & 2.68 &  4.70 & 5.47\\
Avg Clustering Coefficient  & 0.0008 &    0.001 &   0.003\\
 \end{tabular}

\label{tab:net_stats}
\end{table}

\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{nodeDegree}
  \caption{Node degree distributions for protein, gene and
    gene/protein networks.}
  \label{fig:node_degree}
\end{figure}

\begin{figure}[t]
  \centering
  \includegraphics[width=\columnwidth]{nodeClass}

\caption{Edge and node classification based on $w_p^2$.}
  \label{fig:node_class}
\end{figure}

We proceed with a further exploratory analysis of the gene/protein
network. We investigate the contribution of two nodal attributes to the
existence of an edges between the nodes. Following
\citet{katenka2011multi}, we use a simple heuristic based on the weight
vectors to classify the nodes and edges into three classes. For an
edge between the nodes $a$ and $b$, we take one weight vector, say
$\wb_a$, and normalize it to have unit norm. Denote $w_{p}$ the
component corresponding to the protein attribute. Left plot in
Figure~\ref{fig:node_class} shows the values of $w_p^2$ over all
edges. The edges can be classified into three classes based on the
value of $w_p^2$. Given a threshold $T$, the edges for which $w_p^2
\in (0, T)$ are classified as gene-influenced, the edges for which
$w_p^2 \in (1-T, 1)$ are classified as protein influenced, while the
remainder of the edges are classified as mixed type. In the left plot
of Figure~\ref{fig:node_class}, the threshold is set as
$T=0.25$. Similar classification can be performed for nodes after
computing the proportion of incident edges. Let $p_1$, $p_2$ and $p_3$
denote proportions of gene, protein and mixed edges, respectively,
incident with a node. These proportions are represented in a simplex
in the right subplot of Figure~\ref{fig:node_class}. Nodes with mostly
gene edges are located in the lower left corner, while the nodes with
mostly protein edges are located in the lower right corner. Mixed
nodes are located in the center and towards the top corner of the
simplex. Further biological enrichment analysis is possible (see
\citet{katenka2011multi}), however, we do not pursue this here.

\subsection{Uncovering Functional Brain Network}
\label{sec:real_data}

We apply our method to the Positron Emission Tomography dataset, which
contains 259 subjects, of whom 72 are healthy, 132 have mild cognitive
Impairment and 55 are diagnosed as Alzheimer's \& Dementia.  Note that
mild cognitive impairment is a transition stage from normal aging to
Alzheimer's \& Dementia. The data can be obtained from {
 \it http://adni.loni.ucla.edu/}.  The preprocessing is done in the same
way as in \cite{huang2009learning}.

\begin{figure}[!b]
  \centering

    \begin{subfigure}[b]{0.33\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/normGraph1}
      \caption{Healthy subjects}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.33\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/mciGraph1}
      \caption{Mild Cognitive Impairment}
    \end{subfigure}%
    \hfill
    \begin{subfigure}[b]{0.33\textwidth}
      \centering
      \includegraphics[width=\textwidth]{plots/adGraph1}
      \caption{Alzheimer's \& Dementia}
    \end{subfigure}%

  \caption{Brain connectivity networks}
  \label{fig:brain}
\end{figure}

Each Positron Emission Tomography image contains $91 \times 109 \times
91 = 902,629$ voxels. The effective brain region contains $180,502$
voxels, which are partitioned into $95$ regions, ignoring the regions
with fewer than $500$ voxels. The largest region contains $5,014$
voxels and the smallest region contains $665$ voxels. Our
preprocessing stage extracts $948$ representative voxels from these
regions using the $K$-median clustering algorithm. The parameter $K$ is
chosen differently for each region, proportionally to the initial
number of voxels in that region.  More specifically, for each category
of subjects we have an $n \times (d_1 + \ldots + d_{95})$ matrix, where
$n$ is the number of subjects and $d_1 + \ldots + d_{95} = 902,629$ is
the number of voxels. Next we set $K_i = \lceil d_i / \sum_j d_j
\rceil$, the number of representative voxels in region $i$,
$i=1,\ldots,95$. The representative voxels are identified by running
the $K$-median clustering algorithm on a sub-matrix of size $n \times
d_i$ with $K = K_i$.

We inferred three networks, one for each subtype of subjects using the
procedure outlined in $\S$\ref{sec:estimation}.  Note that for
different nodes we have different number of attributes, which
correspond to medians found by the clustering algorithm.  We use the
stability selection~\citep{stability:10} approach to estimate stable
networks.  The stability selection procedure is combined with our
estimation procedure as follows.  We first select the penalty
parameter $\lambda$ in \eqref{eq:max_ll_opt} using cross-validation,
which overselects the number of edges in a network.  Next, we create
$100$ subsampled data sets, each of which contain $80\%$ of the data
points, and estimate one network for each dataset using the selected
$\lambda$.  The final network is composed of stable edges that appear
in at least $95$ of the estimated networks.

We visualize the estimated networks in Figure~\ref{fig:brain}.
Table~\ref{tab:brain_net_stats} provides a few summary statistics for
the estimated networks. Appendix~\ref{appendix:d} contains names of
different regions, as well as the adjacency matrices for networks.
From the summary statistics, we can observe that in normal subjects
there are many more connections between different regions of the
brain.  Loss of connectivity in Alzheimer's \& Dementia has been
widely reported in the literature \citep{greicius2004default,
  hedden2009disruption, andrews2007disruption, wu2011altered}.

Learning functional brain connectivity is potentially valuable for
early identification of signs of Alzheimer's disease.
\citet{huang2009learning} approach this problem using exploratory data
analysis. The framework of Gaussian graphical models is
used to explore functional brain connectivity. Here we point out that
our approach can be used for the same exploratory task, without the
need to reduce the information in the whole brain to one number. For
example, from our estimates, we observe the loss of connectivity in
the cerebellum region of patients with Alzheimer's disease, which has
been reported previously in \citet{sjobeck2001alzheimer}. As another
example, we note increased connectivity between the frontal lobe and
other regions in the patients, which was linked to compensation for
the lost connections in other regions \citep{stern2006cognitive,
  gould2006brain}.

\begin{table}[t]

 \caption{Summary statistics for protein, gene, and gene/protein
   networks ($p = 91$)}

 \begin{tabular}{lrrr}
 ~ & { Healthy }
 & {Mild Cognitive}
& {Alzheimer's \& } \\

& subjects & Impairment & Dementia \\

Number of edges & 116 & 84  &   59 \\
Density & 0.030 & 0.020 & 0.014 \\
Largest connected component & 48 & 27 & 25 \\
Avg Node Degree  & 2.40 & 1.73 &  1.2 \\
Avg Clustering Coefficient  & 0.001  & 0.0023 &    0.0007 \\
 \end{tabular}

\label{tab:brain_net_stats}
\end{table}

\subsection*{Acknowledgments}

{\small
We thank Eric D. Kolaczyk and Natallia V. Katenka for sharing
preprocessed data used in their study with us. Eric P. Xing is
partially supported through the grants NIH R01GM087694 and AFOSR
FA9550010247.  The research of Han Liu is supported by NSF grant
IIS-1116730.
 }


\section{Introduction}

In the year 2012 we celebrate the $100^\mathrm{th}$ anniversary of the
discovery of cosmic rays. But even after a hundred years of research
many questions about the origin, acceleration, and composition of
ultra-high energy cosmic rays remain unanswered. The multi-messenger
approach, combining the information gained from electromagnetic
radiation from radio to TeV photons, charged cosmic rays, and
neutrinos promises to resolve these problems. Neutrinos are ideal
messengers in the sense that they are undeflected by magnetic fields
during their propagation and that they rarely interact, preserving
their initial direction and energy until detected at Earth.

Ultra-high energy (UHE; $E_\nu \gtrsim 100$~PeV) neutrinos, offer a
very rich physics program, including astrophysics, cosmology, particle
physics, and physics beyond the Standard Model:

\begin{itemize}
\item Cosmogenic neutrinos are produced in the interactions of charged
  cosmic rays at ultra-high energies with the cosmic microwave
  background \cite{Berezinsky:1969}, typically within a distance of a
  few ten Mpc of the source \cite{Greisen:1966}. Thus, for very far
  sources, they allow for a good pointing towards the source. The flux
  of UHE cosmogenic neutrinos is very sensitive to the chemical
  composition of the charged cosmic rays (e.g.~\cite{Kotera:2010,
    Ahlers:2012}).
\item Since UHE neutrinos reach us from very high redshifts, their
  flux is also sensitive to the evolution of cosmic ray sources in the
  earlier universe (e.g.~\cite{Kotera:2010}). Resonant $Z$ boson
  production could reveal the cosmic neutrino background and allow us
  to determine the neutrino masses \cite{Ringwald:2008}.
\item Measuring the neutrino flux with different mass overburden,
  e.g. at different zenith angles with an underground detector, will
  allow us to determine the neutrino absorption in the Earth and thus
  probe the neutrino nucleon cross section at high center-of-mass
  energies \cite{Connoly:2006}.
\item There are many theoretical models of physics beyond the Standard
  Model which predict large deviations of the neutrino nucleon cross
  section from the Standard Model at high energies
  (e.g.~\cite{Anchordoqui:2003, Mattingly:2010}). UHE neutrinos will
  allow us to test many of these models.
\end{itemize}

To achieve all these goals we need to measure UHE neutrinos with
reasonable statistics and good energy and direction resolution. This
requires a detector with a volume $\ge 100$~km$^3$. There are
different experimental techniques to build such large scale detectors
which are currently pursued either in running experiments or as
feasibility studies.

Radio detection experiments are looking for short radio pulses in the
hundreds of MHz to GHz frequency range emitted by the electromagnetic
cascade generated in a neutrino interaction. These experiments can be
embedded in radio-transparent, homogeneous media like ice
\cite{Besson:2012} or salt, or use balloons or satellites to observe
large natural ice volumes. Also, observations of the Moon with radio
telescopes are employed to look for neutrino interactions in the lunar
regolith \cite{Seckel:2012}.

Extensive air shower experiments can detect UHE neutrinos either as
highly inclined, ``young'' air showers, where the primary neutrino has
penetrated deep into the atmosphere before interacting, or as up-going
air showers from Earth-skimming neutrinos. The HiRes detector
\cite{Abbasi:2008} and the Pierre Auger Observatory \cite{Abreu:2012}
have used these methods to set upper limits on the flux of UHE
neutrinos.

Finally, acoustic neutrino detectors are searching for ultrasonic
pressure pulses generated in the instantaneous heating and expansion
of the medium induced by electromagnetic and hadronic cascades. Water
\cite{Graf:2012}, ice, salt, and permafrost soil \cite{Nahnhauer:2008}
have been discussed as detection media. In this work we review the
development, status and perspectives of acoustic neutrino detection in
ice.

\section{A brief history of acoustic neutrino detection in ice}

Acoustic neutrino detection in liquids is based on the thermo-acoustic
model \cite{Askaryan:1977, Bowen:1977}: When a neutrino of any flavor
interacts via a charged- or neutral current interaction, a hadronic
and/or electromagnetic cascade develops at the interaction vertex,
which carries a significant amount of the neutrino energy. In a dense
medium this energy is dissipated in a volume of typically $10$~m in
length and a few centimeters in diameter. This leads to an
instantaneous heating of the cascade volume. The corresponding rapid
expansion of the volume propagates as an ultrasonic shock wave
perpendicular to the cascade axis and can be measured as a short
(i.e.~with a broad frequency spectrum), bipolar pressure pulse with a
duration of several ten microseconds. The details of the pulse depend
on the material properties of the medium and on the modeling of the
cascade energy deposition density.

The first ideas about acoustic detection of particles in liquids date
back to the 1950s \cite{Askaryan:1957}. It was then revived in the
1970s and studied in great detail in the context of the DUMAND
project, leading to detailed calculations of the expected acoustic
signals from the thermo-acoustic model \cite{Learned:1979,
  Askaryan:1979} and first measurements with a proton beam from an
accelerator dumped in water \cite{Sulak:1979}.

With the design and construction of the AMANDA optical {Cherenkov}
neutrino telescope at South Pole interest in acoustic neutrino
detection in ice began. Ice, in contrast to water, allows for the
propagation of longitudinal (pressure, p) sound waves and transverse
(shear, s) waves. The formalism of the thermo-acoustic model can be
expanded to the case of solid media and predicts the excitation of
mainly pressure waves by neutrino interactions
(cf.~e.g.~\cite{Salomon:2007}); shear waves can be generated at
impurities in the crystal structure of the medium.

\begin{figure}
  \includegraphics[width=0.44\textwidth]{arena2012_karg_fig01}
  \caption{Laser-induced acoustic pulses in water (top) and ice
    (bottom) generated with identical laser pulses (from
    \cite{Heinen:2011}).\label{fig:laser}}
\end{figure}

Acoustic signals predicted by the thermo-acoustic model scale, for
equal energy deposition densities, with the thermo-elastic properties
of the detection medium. Due to the nearly equal matter densities of
water and ice the energy deposition density from a neutrino
interaction is very similar. From detailed calculations, taking into
account the elastic properties of ice, the amplitudes of
neutrino-induced thermo-acoustic signals in ice are expected to have
amplitudes which are larger by a factor of about four compared to
water \cite{Salomon:2007}. This is supported by measurements in the
Aachen Acoustic Laboratory \cite{Heinen:2012} where laser-induced
acoustic pulses in water and ice have been studied
\cite{Heinen:2011}. Figure~\ref{fig:laser} shows that the scaling of the
signal amplitude from water to ice is compatible with expectations. It
can also be seen that in ice, due to the larger speed of sound, higher
frequency signals are generated.

Phenomenological studies of the ice acoustic properties predicted an
acoustic attenuation length of several kilometers \cite{Price:1993,
  Price:2006} and low background noise \cite{Price:2006}, which would
allow for very large, sparsely instrumented detection
volumes. Subsequent studies favored the radio technique as being more
sensitive than acoustics \cite{Price:1996}. In the same article
(Ref.~\cite{Price:1996}) the possibility of hybrid detection, using
several complementary techniques (radio and optical) is discussed. In
the following years experimental limits on the flux of high energy
neutrinos became more stringent and theoretical flux predictions
decreased accordingly. It became clear that detector volumes $\ge
100$~km$^3$ are required which are difficult to achieve with optical
Cherenkov detectors. Since radio and acoustic signals were expected to
have similar attenuation lengths, hybrid radio-acoustic detectors were
discussed and simulation studies showed very promising results
\cite{Vandenbroucke:2006}.

However, it was clear that the predicted acoustic properties of the
ice need to be tested by in-situ measurements. Different piezoelectric
sensors for use in ice were developed and characterized
\cite{Boser:2006} and sound generation in ice by an accelerator proton
beam was studied \cite{Nahnhauer:2003}. These efforts led to the
construction and deployment of the South Pole Acoustic Test Setup
(SPATS) that will be discussed later in this work.

\section{Measuring acoustic waves in ice}

Building a large acoustic detector with reasonable energy and
direction resolution requires to fully understand the sensor response
in-situ. Sensor sensitivity is not just a single number to convert
incident pressure to output voltage measured at the sensor, but is a
function of incident wave direction, wave mode, temperature, and
possibly other environmental parameters. Since no pre-calibrated
sensors for ice are commercially available that can be used for
relative calibration, extensive studies have been performed to use the
reciprocity calibration method, that does not require a reference
receiver, in ice \cite{Semburg:2011}. The in-situ calibration of
sensors deep in natural glacial ice is even more challenging due to
the limited possibilities of access to the detectors.

In the SPATS project (cf.~next section) it has been tried to factorize
the problem in the laboratory: SPATS sensors have been absolutely
calibrated in water at $0^\circ$C before deployment \cite{Abdou:2011}
and the angular response of the sensor has been determined at
different frequencies \cite{Fischer:2006}. It is not obvious whether the
calibration results can be transferred to operation conditions, where
the sensors are frozen in the deep ice at
South Pole. There, they are subject to low temperatures of
approx.~$-50^\circ$C, increased static pressure, and a different
sensor-medium interface (ice to steel). The different effects have
been studied separately in the laboratory:

\begin{figure}
  \includegraphics[width=0.44\textwidth]{arena2012_karg_fig02}
  \caption{Ratio of the sensitivities measured in water and ice of two
    SPATS sensor channels (from
    \cite{Paul:2011}).\label{fig:sensitivity}}
\end{figure}

\begin{itemize}
\item In air at constant pressure, a sensor has been cooled down from
  $0^\circ$C to $-50^\circ$C and its response to the signal from an
  external transmitter, kept at constant temperature, has been used as
  an estimator for the receiver sensitivity. It has been found that
  the sensor sensitivity increases by a factor of $1.5 \pm 0.2$ when
  the temperature was lowered from $0^\circ$C to $-50^\circ$C
  \cite{Abdou:2011}.
\item In a pressure vessel, filled with an emulsion of water and oil,
  a sensor has been exposed to static pressure up to $100$~bar. A
  transmitter placed outside the pressure vessel and transmitting
  through the steel vessel has been used to determine possible changes
  in sensitivity. The results indicate that the variation of the
  receiver sensitivity is less than $30$\% between $1$ and $100$~bar
  static pressure \cite{Abdou:2011}.
\item The Aachen Acoustic Laboratory \cite{Heinen:2012} allows for the
  production of volumes up to $3$~m$^3$ of clear ice with temperatures
  down to $-25^\circ$C which are ideal for the study of changes in
  sensor sensitivity when going from the water to the ice
  phase. Figure~\ref{fig:sensitivity} shows the ratio of the
  sensitivities of a SPATS sensor measured in water and ice using the
  reciprocity calibration method. It can be seen that the ratio is
  compatible with unity within its errors, indicating that the
  sensitivity does not change in the frequency range relevant for
  acoustic neutrino detection (approx.~from $10$ to $50$~kHz) when the
  sensor is frozen into bulk ice.
\end{itemize}

Assuming that the influences of the environmental effects are
independent, it has been concluded that for the SPATS sensors the
sensitivity in ice will be increased by a factor of $1.5 \pm 0.4$
compared to the pre-deployment calibration in water
\cite{Abbasi:2012}. This factor takes into account the uncertainties
from the temperature and pressure measurements.

It is under investigation how naturally occurring transient noise
events and artificial calibration transmitters can be utilized for
sensor relative calibration and angular response measurements in-situ
\cite{Berdermann:2012}.

\section{Site exploration}

\begin{figure}
  \includegraphics[width=0.65\textwidth]{arena2012_karg_fig03}
  \caption{Overview of the South Pole Acoustic Test Setup (SPATS)
    frozen into the upper 500~m of holes drilled for the IceCube
    neutrino telescope (from \cite{Abdou:2011}).\label{fig:spats_overview}}
\end{figure}

Another important step towards a large scale detector embedded in a
natural detection medium is the full understanding of the signal
propagation properties and backgrounds therein. For an acoustic
experiment this means the determination of the sound speed depth
profile, the attenuation length, the noise level, and possible
transient backgrounds. The sound speed profile determines possible
refraction during the signal propagation that impedes accurate vertex
reconstruction up the existence of multiple solutions. The attenuation
length and noise level will determine the detector geometry required
to achieve a given neutrino energy threshold. Transient noise sources
need to be identified and characterized to separate them from neutrino
induced events. The measurement of the sound speed profile and
transient backgrounds are easier to accomplish in the sense that they
only rely on time information. The attenuation length and noise level
measurement depend on amplitude information and are thus subject to
the calibration challenge discussed in the previous section.

\subsection{SPATS -- Hardware}

To carry out these measurements in the Antarctic ice at the South
Pole, the site of the IceCube neutrino observatory, the South Pole
Acoustic Test Setup (SPATS) has been designed and is successfully
operated since January 2007 \cite{Abdou:2011}. SPATS consists of four
vertical strings that are deployed in the upper $500$~m part of
IceCube bore holes after the installation of the optical IceCube
string. Horizontal baselines between $125$~m and $543$~m are
covered. Each SPATS string is instrumented with seven stages, each
containing an acoustic receiver and a transmitter. The SPATS sensor is
made from a steel housing with three piezoceramic disks pressed to the
inner wall at $120^\circ$ separation for full azimuthal coverage. The
signals are amplified in the sensor module and the differential
analogue signal is transmitted via twisted pair cable to the surface
where it is digitized and time stamped in a String-PC. The data from
all four strings are collected by a Master-PC housed in the IceCube
counting house and are prepared for satellite transmission to the
IceCube central data storage. At two positions an alternative sensor
type, HADES, is installed, where the piezoceramic element is cast in
resin and mounted below the steel housing. HADES is used for
systematic studies of the sensor medium coupling. The SPATS
transmitter consists of a piezoceramic ring cast in resin and frozen
directly into the ice. It is connected to a high voltage pulser which
is protected in a steel housing and steered by the String-PC. A
schematic overview of the SPATS hardware is shown in
Fig.~\ref{fig:spats_overview}.

SPATS is complemented by a mobile acoustic transmitter, called
``pinger'', which can be lowered into freshly drilled water filled
IceCube holes while continuously emitting acoustic pulses with high
stability. The pinger is retrieved from the hole after operation.

\subsection{SPATS -- Results}

\begin{figure}
  \includegraphics[width=0.5\textwidth]{arena2012_karg_fig04a}
  \hspace{0.5cm} \includegraphics[width=0.5\textwidth]{arena2012_karg_fig04b}
  \caption{Acoustic attenuation coefficient for 30~kHz (left) and
    45~kHz (right) measured with the SPATS retrievable pinger in
    different sensor channels (from
    \cite{Abdou:2012}).\label{fig:attenuation}}
\end{figure}

SPATS has measured the speed of sound depth profile in the Antarctic
ice at the South Pole in the depth range from $80$~m to $500$~m using
the retrievable pinger over horizontal baselines of $125$~m
\cite{Abbasi:2010}. The pinger emits longitudinal waves that propagate
through the water column in the drill hole and are then transmitted
into the ice. When the incidence on the water-ice interface is
non-normal, part of the waves energy is transferred into a shear
wave. Thus, the sound speed profile for pressure and shear waves could
be determined. The speed of sound is found to be increasing in the top
$200$~m of the ice where a gradual transition from a snow/air mixture
occurs (firn layer) and is found to be constant below that depth. The best
fit values \cite{Abbasi:2010} for the sound speed $v$ at $375$~m depth
and its gradient $g$ at this depth are for pressure waves

\begin{eqnarray*}
  v_p & = & (3878 \pm 12) \, \mathrm{m} \, \mathrm{s}^{-1} \\
  g_p & = & (0.087 \pm 0.13) \, \frac{\mathrm{m} \,
    \mathrm{s}^{-1}}{\mathrm{m}}
\end{eqnarray*}

\noindent and for shear waves

\begin{eqnarray*}
  v_s & = & (1975.8 \pm 8.0) \, \mathrm{m} \, \mathrm{s}^{-1} \\
  g_s & = & (0.067 \pm 0.086) \, \frac{\mathrm{m} \,
    \mathrm{s}^{-1}}{\mathrm{m}}
\end{eqnarray*}

\noindent Since the gradient is compatible with zero only very little
refraction is expected below the firn layer.

Pinger data have also proven very valuable to determine the signal
amplitude attenuation length \cite{Abbasi:2011} in the frequency range
from $10$ to $30$~kHz. This analysis requires the comparison of the
signal observed at different distances. To reduce systematic
uncertainties from the sensor absolute calibration and angular
response, attenuation lengths are derived for each sensor channel. For
this the pinger, which produces highly reproducible pulses, is
deployed in different IceCube drill holes at increasing distances, but
aligned in direction, so that only a small range of the angular
response of the sensor is probed. Averaging over all available sensor
channels leads to an signal amplitude attenuation length
\cite{Abbasi:2011} of

\begin{equation*}
\langle \lambda \rangle = 312^{+68}_{-47} \, \mathrm{m}
\end{equation*}

\noindent This value is much smaller than the several kilometers
initially expected from theoretical calculations \cite{Price:2006}. To
study this discrepancy a modified pinger has been constructed that
emits gated sine bursts at different frequencies ($30$, $45$, and
$60$~kHz). It has been used to study the frequency dependence of the
attenuation length. Absorption, if it would be the main attenuation
mechanism, is expected to be frequency independent, whereas the
scattering coefficient would increase $\propto f^4$, where $f$ is the
frequency of the signal \cite{Price:2006}. The data from the modified
pinger are analyzed with the same methods as described above
\cite{Abdou:2012}. Figure~\ref{fig:attenuation} shows the attenuation
coefficient measured in the different sensor channels for $30$~kHz and
$45$~kHz as well as the mean value and spread of the data points. The
signal-to-noise ratio for $60$~kHz has been too poor to extract an
attenuation length since the transmit response of the pinger
piezoelectric element is small at this frequency. It can be seen that
the two values are compatible with a frequency independent attenuation
length and that an $f^4$ frequency dependence is hard to reconcile
with the data.

SPATS unbiased noise data, that are recorded every hour for $0.1$~s
have been used to study the background noise level. It has been shown
that the noise is Gaussian and that the RMS is very stable over time
\cite{Abbasi:2012}. Using the corrections on the sensor calibration
discussed in the previous section, the absolute noise level below the
firn layer has been estimated to be $14$~mPa, corresponding to the
signal expected from a $10^{11}$~GeV neutrino interacting at a
distance of $1000$~m to the sensor \cite{Abbasi:2012}. However, the
systematic uncertainty on the noise level is still large, making
additional measurements with sensors pre-calibrated in ice desirable.

\begin{figure}
  \includegraphics[width=0.475\textwidth]{arena2012_karg_fig05}
  \caption{Upper limits (90\% C.L.) on the all-flavor flux of
    ultra-high energy neutrinos set by acoustic neutrino detection
    experiments in water (SAUND II, ACoRNE) and ice (SPATS12). Several
    models on for the flux of cosmogenic neutrinos and the most
    stringent limits set by radio experiments are shown for comparison
    (adapted from \cite{Abbasi:2012}).\label{fig:limits}}
\end{figure}

An analysis of transient events triggering all four strings of the
SPATS detector revealed only man-made sources: all events were
reconstructed to either re-freezing IceCube drill holes or Rodriguez
Wells, large caverns in the ice used to circulate the water for the
IceCube hot water drill system \cite{Abbasi:2012}. The absence of
unidentified transient noise sources in the deep ice allowed setting
an upper limit on the flux of ultra-high energy
neutrinos. Figure~\ref{fig:limits} compares the limit set by SPATS to
other experiments and the expected cosmogenic neutrino flux. It has to
be kept in mind that none of the acoustic experiments shown there were
designed as neutrino detectors, but were built for site-exploration
(SPATS) or parasitically use existing military hydrophone arrays to
search for neutrino induced signals (SAUND II, ACoRNE).

\section{A path forwards}

The SPATS results show that acoustic neutrino detection in the glacial
ice at the South Pole is feasible. The acoustic attenuation length is
shorter than for radio signals
but of comparable magnitude which opens the possibility to design a
hybrid radio/acoustic detector. Due to the intrinsically higher energy
threshold of the acoustic technique a possible scenario is to use the
radio sub-array to trigger the acoustic sub-array up to energies where
the acoustic detector becomes fully efficient by itself. In this case
a single in-time hit in an acoustic sensor can already be highly
significant since all known backgrounds produce either only radio
emission or only sound emission.

The construction of such a detector at the South Pole will require the
deployment of a few hundred strings per $100$~km$^2$ instrumented area
(estimating a spacing of $500$~m between strings based on the measured
attenuation lengths for radio and acoustic) reaching below the firn
layer. A design which is scalable in size is desirable so that, once
the magnitude of the flux of ultra-high energy neutrinos is measured,
the detector can be expanded in size to accumulate sufficient event
statistics to answer the physics questions discussed in the
introductory section.

\begin{figure}
  \includegraphics[width=0.35\textwidth]{arena2012_karg_fig06}
  \caption{Schematic of an autonomous drilling and deployment probe
    based on the IceMole concept (cf.~text). The cable and sensors are
    stored in and deployed from the probe, so that the hole is allowed
    to refreeze immediately after the passage of the probe (from
    \cite{Laihem:2012}).\label{fig:drilling}}
\end{figure}

A strong R\&D program will be required to realize a new large area
hybrid detector. Apart from South Pole there are several other sites
in Antarctica where scientific infrastructure exists and which are
worth evaluating:\\
- On the Ross Ice Shelf the ARIANNA radio neutrino
detector\footnote{http://arianna.ps.uci.edu/} is currently under
construction.\\
- Concordia Station at Dome C is located on top of more than $3$~km of
very cold ice which is favorable for acoustic and radio signal
propagation

The installation of a new detector consisting of several hundred to a
few thousand strings will also require new techniques for drilling and
deployment. One possibility that is being discussed is the use of
autonomous drilling- and/or melting-probes similar to the IceMole
\cite{Dachwald:2011} prototype. To minimize human intervention in the
drilling and deployment procedure, the sensors and cable would be
stored in and deployed from the disposable icecraft which would remain
in the ice at the bottom of the string after deployment. This allows
for the hole to refreeze immediately after the passage of the probe. A
schematic of the procedure is shown in Fig.~\ref{fig:drilling}.

A large area detector at a remote site will also require new concepts
for calibration, communication and power supply. The large overall
extent and the large distances between the components prohibit a fully
cabled design. Wireless data transmission and power generation at the
site of the component by e.g.~wind or solar power are
required. Ideally one would have a combination of wind and solar power
since beyond the polar circles solar power is unavailable half of the
year. Valuable lessons can be learned from large area experiments
which are already in operation, like the Pierre Auger Observatory
\cite{Abraham:2004}.

\section{Conclusions}

Ultra-high energy neutrinos offer a vast physics program covering
astrophysics, cosmology, particle physics, and physics beyond the
Standard Model. The acoustic neutrino detection technique has made
large advances over the last few years: sensors have been designed and
their behavior in ice is largely understood. The acoustic properties
of the South Pole glacier have been measured and found to be suitable
for neutrino detection. It is expected that acoustic can play an
important part in a future hybrid neutrino telescope in ice. To
realize this a strong R\&D program has to be established and first
promising studies have already been presented.

\begin{theacknowledgments}
  T.K.~is supported by the ``Helmholtz Alliance for Astroparticle
  Physics HAP'' funded by the Initiative and Networking Fund of the
  Helmholtz Association.
\end{theacknowledgments}

\bibliographystyle{aipproc}