text:
\section{Introduction}

Networks arise in many areas of research and applications, which come in all shapes and sizes.  The most studied and best understood are static network models. Many other network models are also in existence, but have been less studied. One such example is the multi-layer networks, which are a powerful representation of relational data, and commonly encountered in contemporary data analysis (\cite{kivela2014multilayer}). The nodes in a multi-layer network represent the entities of interest and the edges in different layers indicate the multiple relations among those entities. Examples include brain connectivity networks, world trading networks, gene-gene interactive networks and so on.
In this paper, we focus on the multi-layer networks with the same nodes set of each layer and there are no edges between two different layers.

The study on multi-layer networks has received an increasing interest. Considering the dependency among the different layers, \cite{paul2016consistent} derives consistency results for the community assignments from the maximum likelihood estimators in two models. Consistency properties of various methods for community detection under the multi-layer stochastic block model are investigated in \cite{paul2017spectral}. Three different matrix factorization-based algorithms are employed in \cite{tang2009clustering}, \cite{nickel2011three} and \cite{dong2012clustering} separately. Common community structures for multiple networks are identified via two spectral clustering algorithms with theoretical guarantee in \cite{bhattacharyya2018spectral}.  In \cite{arroyo2019inference}, authors introduce the common subspace independent-edge multiple random graph model to describe a heterogeneous collection of networks with a shared latent structure and propose a joint spectral embedding of adjacency matrices to simultaneously and consistently estimate underlying parameters for each graph. Consistency results for a least squares estimation of memberships under the multi-layer stochastic block model framework are derived in \cite{lei2019consistent}. Several literature focus on recovering the network from a collection of networks with edge contamination. The original network is estimated from multiple noisy realizations utilizing community structure in \cite{le2018estimating} and low-rank expectation in \cite{levin2019recovering}. A weighted latent position graph model contaminated via an edge weight gross error model is proposed in \cite{tang2017robust} with an estimation methodology based on robust $L_q$ estimation followed by low-rank adjacency spectral decomposition.

In applications, a random effects stochastic block model is proposed by \cite{paul2018random} for the neuroimaging data and a statistical framework with a significance and a robustness test for detecting common modules in the Drosophila melanogaster dynamic gene regulation network is proposed in \cite{zhang2017finding}.

Most of the literature about community detection in multi-layer networks is limited to consistent membership setting, which means all the layers carry information about the same community assignment. However, in reality, different layers may have different community structures. For instance, in a social network, layers related with sports (people connected with the same sport hobbies) may have different community structure with layers about movie taste (people connected with similar movie taste). Understanding the large-scale structure of multi-layer networks is made difficult by the fact that the patterns of one type of link may be similar to, uncorrelated with, or different from the patterns of another type of link. These differences from layer to layer may exist at the level of individual links, connectivity patterns among groups of nodes, or even the hidden groups themselves to which each node belongs. In \cite{de2017community}, authors pointed out that, in order to do community detection on multi-layer networks, it is crucial to know which layers have related structure and which layer are unrelated, since redundant information across layers may provide stronger evidence for clear communities than each layer would on its own. Such situation is not clearly discussed in the works mentioned above. Although, in \cite{matias2017statistical}, authors introduced community structure variety as time varying, it is hard to be applied in general multi-layer networks without time ordering.

In this paper, we introduce a general framework, i.e., mixture multi-layer stochastic block model (MMSBM), and propose a tensor-based algorithm (TWIST) to reveal both global/local memberships of nodes, and memberships of layers.
To fix ideas, we start with a simple motivating example, illustrated in Figure \ref{fig:toyexample}.  We have $L=3$ layers of networks $\{{\cal G}_1,{\cal G}_2,{\cal G}_3\}$, each containing 3 local communities. The 3 networks are of $m=2$ types: $\{{\cal G}_2\}$ and $\{{\cal G}_1,{\cal G}_3\}$. The community structure differs between  $\{{\cal G}_2\}$ and the other two networks as some members in the third community $g_{23}$ are in the second one $g_{12}$ in $\{{\cal G}_1\}$ and $\{{\cal G}_3\}$.
Viewing the 3 layers of networks together, we notice that there are 4 global communities, in which members stay in all layers throughout. Clearly, the global communities are related to, but different from the local ones in each network.  Our interest lies in detecting both local as well as global community structures, which are of great value in theory and practice.  There is an increasing literature on the global community structure as mention earlier. However, to the best of our knowledge, there is no systematic investigation into detecting local and global community structure together.

\begin{figure}
\centering
	\includegraphics[scale=.4]{figs/toy_example}
	\caption{A toy example.}
	\label{fig:toyexample}
\end{figure}

Our line of attack can be illustrated via the following diagram in Figure \ref{fig:process}. First, we pool adjacency matrices from all layers of networks to form a tensor (multi-way array), and then apply the TWIST (to be introduced later) to obtain the global community structure as well as labels of each layer. We then group the layers of networks with the same labels, which will be used to detect local community structures. Details will be unfolded next.

\begin{figure}
\centering
	\includegraphics[scale=.4]{figs/process}
	\caption{The general procedure of TWIST. }
	\label{fig:process}
\end{figure}
The main contributions of this paper are summarized as follows.

First, we propose a very general model to handle the type of problems discussed above. To be more specific, we will introduce the so-called mixture multi-layer stochastic block model (MMSBM), which can characterize the different community structures among different layers of the multi-layer network. In some way, the MMSBM resembles the relatively well studied multi-layer stochastic block model (MLSBM) \cite{paul2017spectral, han2015consistent, paul2016consistent}. However, the MMSBM is more general in that it allows the multi-layer network to contain different block structures. Thus, the MMSBM not only allows each layer to have different community structures, but also can maintain the consistent structure in the network.

Secondly, we propose a tensor-based method to study the MMSBM. The approach is referred to as the Tucker decomposition with integrated SVD transformation (TWIST). Unlike earlier approaches for multi-layer network analysis, TWIST can uncover the clusters of layers, the local and global membership of nodes simultaneously. On the theoretical front, we prove for MMSBM that TWIST can consistently recover the layer labels and global memberships of nodes under near optimal network sparsity conditions. In addition, network labels can be exactly recovered under a slightly stronger network sparsity condition.
To the best of our knowledge, this is the first systematic study on statistical guarantees about community detection in a mixture multi-layer networks using tensor decomposition. Our primary technical tool is a sharp concentration inequality of sparse tensors which might be of independent interest.

Finally, two real-world applications of the proposed methodology demonstrated to be a powerful tool in analysing multi-layer networks. The algorithm is easy to use and can help practitioners  quickly uncover interesting findings, which would otherwise be difficult by using other tools.

The rest of the paper is organized as follows. Section \ref{sec:model} introduces the mixture multi-layer stochastic block model (MMSBM) for describing the mixture structure. A new algorithm, the TWIST, is proposed in Section \ref{sec:method}. We explore the theoretical properties of the TWIST under the MMSBM in Section \ref{sec:theory}. Moreover, we make comparisons between our main results and the cutting-edge theoretical results. The advantages of the proposed method is numerically evaluated with several simulations in Section \ref{sec:simulation} and two real data examples in Section \ref{sec:realdata}. Section \ref{sec:discussion} gives concluding remarks and discussions. All the proofs are shown in the supplement.

\section{Model framework}
\label{sec:model}
\subsection{Mixture multi-layer stochastic block model (MMSBM)}

The observed data contains $L$-layers of networks on the same set of vertices: ${\cal  V}} \def\cV{{\cal  V}=[n]:=\{1,2,\cdots,n\}$:
$$
 {\cal  G}} \def\cG{{\cal  G} = \{{\cal G}_l: \ l =1,\cdots, L\}.
 $$
Assume that these networks are generated from a mixture of $m$ latent networks with probability $\pi = (\pi_1,...,\pi_m)$.  Denoting $\ell_l \in \{1,\cdots,m\}$ as a latent label of ${\cal G}_l$ with $1\le l \le L$, then $$ \PP(\ell_l=j)=\pi_j,  \quad \mbox{with}  \ \ \sum_{j=1}^m\pi_j=1$$

Assume that each of the $m$ classes of networks  satisfies the {\it stochastic block model} (SBM). More specifically, for $j\in[m]$, the $j$-th class SBM is described by the {\it membership matrix} $Z_j\in\{0,1\}^{n\times K_j}$ and the probability matrix $B_j\in[0,1]^{K_j\times K_j}$ where $K_j$ is the number of communities.  Each row of $Z_j$ has exactly one entry which is non-zero. For simplicity, we denote

\begin{itemize}
\item SBM$(Z_j, B_j) = $ the $j$-th SBM with parameter $Z_j$ and $B_j$, $j=1,...,m$.

\item ${\cal  V}} \def\cV{{\cal  V}^j_k = $ the $k$-th community in the $j$-th SBM.  So ${\cal  V}} \def\cV{{\cal  V}^j_k\subset {\cal  V}} \def\cV{{\cal  V}$ and $\cup_{k=1}^{K_j} {\cal  V}} \def\cV{{\cal  V}_k^j = {\cal  V}} \def\cV{{\cal  V}$.

\item $L_j = \#\{l: \ell_l=j, 1\leq l\leq L\} =$ the number of layers generated by SBM$(Z_j, B_j)$.
Clearly, $L=\sum_{j=1}^m L_j$.

\item $\mathring{K}=K_1+\cdots+K_m$ and $\LL=\{\ell_l\}_{l=1}^L$ and $\VV^j:=\{{\cal  V}} \def\cV{{\cal  V}^j_k\}_{k=1}^{K_j}$.
\end{itemize}

The observed adjacency matrix $A_l \in\{0,1\}^{n\times n}$ of ${\cal  G}} \def\cG{{\cal  G}_l$ obeys Bernoulli distribution:
\begin{equation}\label{eq:Al_bern}
A_{l}(i_1,i_2)\sim_{iid} {\rm Bern}\big(Z_{\ell_l}(i_1,:)B_{\ell_l}Z_{\ell_l}(i_2,:)^{\top}\big)
 \  \Longleftrightarrow  \
A_l \sim_{iid} {\rm Bern}\big(Z_{\ell_l}B_{\ell_l}Z_{\ell_l}^{\top}\big)
\end{equation}
for all $i_1\leq  i_2 \in [n]$, where $Z(i,:)$ denotes the $i$-th row of $Z$.

The resulting model is referred to as ``mixture multi-layer stochastic block model" (MMSBM).

\subsection{Adjacency tensor and its decomposition}
Observing the $L$ layers of networks, we define the adjacency tensor $\bA\in\RR^{n\times n\times L}$ so that $\bA$'s $l$-th slice
$$
A(:,:,l)=A_l,\quad \forall 1\leq l\leq L.
$$
See \cite{kolda2009tensor} for an introduction to tensor and applications.
It follows from (\ref{eq:Al_bern}) that
$$
\EE (A_l\vert \ell_l)=Z_{\ell_l}B_{\ell_l}Z_{\ell_l}^{\top},\quad \forall 1\leq l\leq L,
$$
from which we can derive the following tensor representation, whose proof is in the Appendix.

\begin{lemma}[Tensor representation]  \label{lem:EA_decomp}
We have
\begin{equation} \label{eq:barZ}
\EE(\bA |\LL)=\bB  \times_1 \bar{Z} \times_2 \bar{Z} \times_3 W,
\end{equation}
where $\LL=\{\ell_l\}_{l=1}^L$ and
\begin{itemize}
\item $\bar{Z}=(Z_1,Z_2,\cdots,Z_m)\in\{0,1\}^{n\times \mathring{K}}$ is the {\bf global} membership matrix, whereas each $Z_j$ is the {\bf local} membership matrix,

\item $W=(e_{\ell_1},e_{\ell_2},\cdots,e_{\ell_L})^{\top} \in\{0,1\}^{L\times m}$ is the network label matrix
with each row of $W$ having exactly one non-zero entry, and $e_j\in\RR^m$ being the $j$-th canonical basis vector,

\item  $\bB \in\RR^{\mathring{K}\times \mathring{K}\times m}$ is a 3-way probability tensor whose $j$-th frontal slide is
$$
B(:,:,j)={\rm} diag(0_{K_1},\cdots,0_{K_{j-1}},B_j,0_{K_{j+1}},\cdots,0_{K_m}),\quad 1\leq j\leq m
$$
with $0_K$ being a $K\times K$ zero matrix.
\end{itemize}
\end{lemma}

 \subsection{Local versus global memberships via Tucker decomposition}

 The matrix $\bar{Z}$ defined in Lemma~\ref{lem:EA_decomp} suggests the existence of {\it global community structures}.
We say that two nodes $i_1$ and $i_2$ belong to the same global community if and only if they belong to the same local community for all the $m$ classes of SBM, i.e.,
$$ \bar{Z}(i_1,:)=\bar{Z}(i_2,:). $$
Let $\bar{K}$ denote the number of global communities. Clearly, $\max_j K_j\leq \bar{K}\leq \prod_j K_j$. Denote $\bar{\VV}=\{\bar{{\cal  V}} \def\cV{{\cal  V}}_k\}_{k=1}^{\bar{K}}$ the global community clusters such that $\cup_{k=1}^{\bar{K}}\bar{{\cal  V}} \def\cV{{\cal  V}}_k={\cal  V}} \def\cV{{\cal  V}$. Therefore, for two nodes $i_1 \ne i_2$,
 \begin{equation}\label{eq:barVk}
 \{i_1, i_2\} \in\bar{{\cal  V}} \def\cV{{\cal  V}}_k \qquad \Longleftrightarrow \qquad \{i_1, i_2\} \in {\cal  V}} \def\cV{{\cal  V}^j_{k_j},   \quad k_j\in[K_j], \ \ \forall j\in[m],
 \end{equation}

 Let $r={\rm rank}(\bar{Z})$ denote the rank of $\bar{Z}$.  We hereby write the thin SVD of $\bar{Z}$ as
 \begin{equation}\label{eq:barZ_svd}
 \bar{Z}=\bar{U}\bar{D}\bar{R}^{\top}
 \end{equation}
where $\bar{U}\in\RR^{n\times r}$, $\bar{R} \in \RR^{(\mathring{K})\times r}$ have orthonormal columns, and $\bar D$ is the singular value diagonal matrix
$$
\bar{D} = diag(\sigma_1(\bar D), \cdots, \sigma_r(\bar D))  \in \RR^{r\times r},
\qquad \sigma_1(\bar D)\geq \cdots\geq \sigma_r(\bar D)>0.
$$
The global community structure can be checked by $\bar{U}$ as in Lemma~\ref{lem:barU_sep}.
 \begin{lemma}\label{lem:barU_sep}
 For $i_1\in\bar{{\cal  V}} \def\cV{{\cal  V}}_{k_1}$ and $i_2\in\bar{{\cal  V}} \def\cV{{\cal  V}}_{k_2}$ with $k_1\neq k_2$, then,
 $$
 \|\bar{U}(i_1,:)-\bar{U}(i_2,:)\|_{\ell_2}\geq \frac{1}{\sigma_1(\bar{D})}.
 $$
 \end{lemma}

 By (\ref{eq:barZ_svd}), the population adjacency tensor $\EE(\bA|\LL)$ admits the {\it Tucker decomposition} as
 \begin{equation}\label{eq:A_tucker}
 \EE(\bA|\LL)=\bar{\bC}\times_1\bar{U}\times_2\bar{U}\times_3\bar{W}
 \end{equation}
 where the core tensor $\bar \bC\in\RR^{r\times r\times m}$ is defined by
 \begin{equation}\label{eq:barC}
 \bar{\bC}=\bB\times_1(\bar{D}\bar{R}^{\top})\times_2 (\bar{D}\bar{R}^{\top})\times_3 D_L^{1/2}
 \end{equation}
 and $\bar{W}=WD_L^{-1/2} \in \RR^{L\times m}$ so that $\bar{W}^{\top}\bar{W}=I_m$, and the diagonal matrix
 $$
 D_L={\rm diag}(L_1,L_2,\cdots,L_m).
 $$
We assume that $\bar{\bC}$ has Tucker ranks $(r,r,m)$. Further assume that $m\leq r$, which is reasonable as one new type of network will introduce at least one new global community.

The decomposition (\ref{eq:A_tucker}) shows that the singular vectors of $\EE(\bA|\LL)$ contain the latent network information.
More exactly, the singular vectors in the $1$-st dimension of $\EE(\bA|\LL)$ could identify the global community structures and singular vectors in the $3$-rd dimension could identify the latent network labels.  After identifying the latent network labels, a post-processing procedure can identify the local community structures.

\section{Methodology: TWIST} \label{sec:method}

By observing the multi-layer networks $\{{\cal  G}} \def\cG{{\cal  G}_l\}_{l=1}^L$ satisfying model (\ref{eq:Al_bern}), our goals are to:
\begin{enumerate}
\item[(1)] recover the global community structures of vertices $\{\bar{{\cal  V}} \def\cV{{\cal  V}}_k\}_{k=1}^{\bar{K}}$;
\item[(2)] identify network classes $\{\ell_l\}_{l=1}^L$, and grouping networks with the same class;
\item[(3)] recover the local community structures of the vertices $\VV^j:=\{{\cal  V}} \def\cV{{\cal  V}_{k}^j:  k\in [K_j]\}$ for all $j\in[m]$.
\end{enumerate}
Note that in order to efficiently recover the local community structures, it is necessary to first identify the network classes. As a result, task (3) usually follows from task (2).

By the decomposition of oracle tensor (\ref{eq:A_tucker}), the singular vectors $\bar{U}$ contains information of global memberships since its column space comes from $\bar{Z}$. Additionally, the singular vectors $\bar{W}$ contains information of network classes. Therefore, task (1) and task (2) are both related with the Tucker decomposition of oracle tensor $\EE(\bA|\LL)$. Since the oracle is unavailable, we seek a low-rank approximation of $\bA$.

\subsection{Tucker decomposition with integrated SVD transformation (TWIST)}

In order to utilize the low rank structure of the tensor and the non-negative property of the elements, we propose a new algorithm called Tucker decomposition with integrated SVD transformation (TWIST). The general procedure is summarized below and illustrated in Figure~\ref{fig:process}.

\begin{itemize}
\item {\bf Step 1: Decomposition of adjacency tensor} \\
Apply the regularized tensor power iterations to $\bA$ to obtain its low-rank approximation. The outputs are $\widehat U$ and $\widehat W$. Details are given in Algorithm~\ref{algo:Power Iterations}.

\smallskip

\item {\bf Step 2: Global memberships} \\
Apply the standard K-means algorithm on the rows of $\widehat U$ to identify the global community memberships and output $\widehat{\bar \VV}=\{\widehat{\bar{{\cal  V}} \def\cV{{\cal  V}}}_k\}_{k=1}^{\bar{K}}$.
\smallskip

\item {\bf Step 3: Network classes}  \\
Use the rows of $\widehat W$ to identify the network classes and output the network classes: $\widehat\LL=\{\hat\ell_l\in[m]\}_{l=1}^L$. We can use either the standard K-means or the sup-norm related algorithm (Algorithm~\ref{algo:network}).

\smallskip

\item {\bf Step 4: Local memberships}  \\
We can find the local membership $\VV^j=\{{\cal  V}} \def\cV{{\cal  V}_{k}^j\}$ by focusing on networks with the same labels (\cite{lei2015consistency, rohe2011spectral}).  More precisely, for each $j\in\{1,...,m\}$, we can apply K-means either
\begin{itemize}
\item to the sum of those networks with the same label $\sum_{l: \hat\ell_l=j} A_l$, or
\item to the sub-tensor $A(:,:,\{l: \hat\ell_l=j\})$), those slides with the same labels.
\end{itemize}
Outputs are $\widehat\VV^j=\{\widehat{{\cal  V}} \def\cV{{\cal  V}}^j_{k}\}_{k=1}^{K_j}$.
\end{itemize}

\smallskip

\begin{algorithm}
	\caption{Regularized power iterations for sparse tensor decomposition}
	\label{algo:Power Iterations}
	\begin{algorithmic}[20]
		\Require
		$\bA_{(n\times n \times L)}$, warm initialization $\widehat U^{(0)}$ and $\widehat W^{(0)}$\\
		\hspace{1.5cm}  maximum iterations ${\rm iter}_{\max}$ and regularization parameters $\delta_1,\delta_2>0.$
		\Ensure
		$\widehat{U}$ and $\widehat W$
		\State Set counter ${\rm iter} = 0.$
        \While{${\rm iter} < {\rm iter}_{\max}$}
         \State Regularization: $\widetilde U^{({\rm iter})}\leftarrow {\cal  P}} \def\cP{{\cal  P}_{\delta_1}(\widehat{U}^{({\rm iter})})$ and $\widetilde W^{({\rm iter})}\leftarrow {\cal  P}} \def\cP{{\cal  P}_{\delta_2}(\widehat{W}^{({\rm iter})})$ by (\ref{eq:reg}).
        \State ${\rm iter} \leftarrow {\rm iter} + 1$
          \State Set $\widehat U^{({\rm iter})}$ to be the top $r$ left singular vectors of
          ${\cal  M}} \def\cM{{\cal  M}_1\big(\bA \times_{2}\widetilde{U}^{({\rm iter}-1)\top} \times_3 \widetilde{W}^{({\rm iter}-1)\top}\big).$
          \State set $\widehat W^{({\rm iter})}$ to be the top $m$ left singular vectors of ${\cal  M}} \def\cM{{\cal  M}_3\big(\bA \times_1 \widetilde{U}^{({\rm iter}-1)\top} \times_2 \widetilde{U}^{({\rm iter}-1)\top} \big)$.
        \EndWhile
        \State Return $\widehat{U}\leftarrow \widehat{U}^{({\rm iter})}$ and $\widehat{W}\leftarrow\widehat{W}^{({\rm iter})}$.
    	\end{algorithmic}
\end{algorithm}

\begin{algorithm}
\caption{Network clustering by sup-norm K-means}
\label{algo:network}
\begin{algorithmic}[20]
\Require $\widehat W$, number of clusters $m$ and threshold $\varepsilon\in(0,1)$
\Ensure  Network labels $\widehat\LL=\{\hat\ell_l\}_{l=1}^L$
\State Initiate ${\cal  C}} \def\cC{{\cal  C}\leftarrow\{1\}$, $\hat\ell_1\leftarrow1$, $k\leftarrow 1$ and $l\leftarrow 2$.
\While{$l\leq L$}
\State Compute $j\leftarrow \argmin_{j\in{\cal  C}} \def\cC{{\cal  C}} \|\widehat W(l,:)-\widehat{W}(j,:)\|$
\If{$\|\widehat W(l,:)-\widehat W(j,:)\|>\varepsilon$}
\State $k\leftarrow k+1$
\State $\hat\ell_l\leftarrow k$
\State ${\cal  C}} \def\cC{{\cal  C}\leftarrow {\cal  C}} \def\cC{{\cal  C}\cup \{l\}$
\Else
\State $\hat\ell_l\leftarrow \hat\ell_j$
\EndIf
\State $l\leftarrow l+1$
\EndWhile
\If{$k>m$ (or $k<m$)}
\State  Set $\varepsilon\leftarrow 2\varepsilon$ (or set $\varepsilon\leftarrow \varepsilon/2$); Re-run the algorithm.
\Else
\State Output $\widehat\LL=\{\hat\ell_l\}_{l=1}^L$
\EndIf
\end{algorithmic}
\end{algorithm}

\subsection{Features about TWIST}

There are several key features concerning the TWIST.

\subsubsection*{Warm starts for $\widehat U^{(0)}$ and $\widehat W^{(0)}$ in Algorithm 1}

Computing the optimal low-rank approximation of a tensor $\bA$ is NP-hard  in general; see \cite{hillar2013most}.
Algorithms with random initializations are almost always trapped in non-informative local minimals which can be nearly orthogonal to the truth, see \cite{arous2019landscape}.  To avoid these issues, tensor decomposition algorithms usually run from a warm starting point \cite{zhang2018tensor, xia2019polynomial, xia2017statistically, richard2014statistical, jain2014provable, ke2019community, zhang2019cross, sun2017provable, cai2019nonconvex, wang2018learning}.

In Section~\ref{sec:init}, we will introduce a warm initialization algorithm, obtained by applying a spectral method for initializing $\widehat{U}^{(0)}$ by summing up all the network layers.
Initialization of $\widehat W^{(0)}$ is easy whenever $\widehat U^{(0)}$ is available. We show in Lemma \ref{lem:initialization} that these initializations can indeed improve estimation accuracy.

\subsubsection*{Regularized power iterations for sparse tensor decomposition}

Adjacency matrices from some layers are often very sparse and the individual layers are even disconnected graphs. For example, in the Malaria parasite genes networks given in Section~\ref{sec:realdata},  three out of nine networks are very sparse and and disconnected.  Under these circumstances, the popular tensor power iteration algorithm, i.e., high-order orthogonal iterations (HOOI, see \cite{sheehan2007higher}) may not work. In fact, its statistical optimality was proved by \cite{zhang2018tensor} only for dense tensors, while its properties on sparse random tensors remain much more challenging.

To handle sparse random tensors, we employ a regularized tensor power iteration algorithm in Algorithm~\ref{algo:Power Iterations}, which was used in \cite{ke2019community} to deal with sparse hypergraph networks.  Regularizations to singular vectors $\widehat U^{(t)}, \widehat W^{(t)}$ are applied before each power iteration. We take $\widehat U^{(t)}$ for example as $\widehat W^{(t)}$ can be treated similarly. The regularization ${\cal  P}} \def\cP{{\cal  P}_{\delta}(U)$ is done by
\begin{equation}\label{eq:reg}
{\cal  P}} \def\cP{{\cal  P}_{\delta}(U)={\rm SVD}_r(U_{\star}) \quad {\rm where}\quad U_{\star}(i,:):=U(i,:)\cdot\frac{\min\{\delta,\|U(i,:)\|\}}{\|U(i,:)\|}\quad i\in[n],
\end{equation}
The effect of regularization is to dampen the influence of ``large" rows of $\widehat U^{(t)}$, which is due to the communities of small sizes, besides stochastic errors.
Following Lemma~\ref{lem:incoh}, the true singular vectors $\bar{U}$ are incoherent with $\max_j\|e_j^{\top}\bar{U}\|=O(\sqrt{r/n})$ if community sizes are balanced. In practice, we suggest
\begin{equation}\label{eq:delta_algo}
\hat\delta_1 =2\sqrt{r}\cdot \frac{\max_{1\leq i\leq n}{\rm deg}_i}{\sqrt{\sum_i {\rm deg}_i^2}}
\quad {\rm and}\quad\hat \delta_2=2\sqrt{m}\cdot \frac{\max_{1\leq l\leq L} {\rm neg}_l}{\sqrt{\sum_l {\rm neg}_l^2}}
\end{equation}
where the node degree ${\rm deg}_i=\sum_{j, l}A(i,j,l)$ and layer degree ${\rm neg}_l=\sum_{i,j}A(i,j,l)$.

\subsubsection*{K-means with sup-norm distance}

Clearly, the accuracy of local membership clustering (Step 4) hinges on the reliability of layer labelling. In Algorithm~\ref{algo:network},
a sup-norm version of K-means is applied to the singular vectors $\widehat W$ obtained from Algorithm~\ref{algo:Power Iterations}, and then outputs the network labels $\widehat\LL=\{\hat\ell_l\}_{l=1}^L$. The sup-norm K-means has recently been extensively investigated and shown to perform well in network community detection. See, e.g. \cite{chen2019spectral,abbe2017entrywise, kim2018stochastic} and references therein.

The rationale of Algorithm~\ref{algo:network} is that when the rows of $\bar{W}$ are well-separated (similar to Lemma~\ref{lem:barU_sep}), a row-wise screening of $\widehat W$ can immediately recover the true network labels as long as the row-wise perturbation bound of $\widehat W-W$ is small enough. As shown in Section~\ref{sec:main_theory}, Algorithm~\ref{algo:network} guarantees exact clustering of networks under weak conditions.

In Steps 2-4, one could use alternative methods other than K-means clustering, which might improve its performances. For example, we can use DBSCAN, Gaussian mixture model, the SCORE method (\cite{jin2015fast, jin2017estimating, ke2019community}).

\section{Preliminary Results} \label{sec:theory}

\subsection{Notations and definitions}

For ease of exposition, we introduce the following notations.
\begin{itemize}
\item Denote $c, c_j, C, C_j,C_j',\ j\ge 1$ as generic constants, which may vary from line to line.
\item Denote $e_k$ as the $k$-th canonical basis vector in Euclidean space (i.e., with only the $k$-th entry equal to 1 and others 0),
whose dimension depends on each context.

\item For a matrix $M=\{m_{ij}\}$, let

$\sigma_i(M) = $ the $i$-th largest singular value of $M$,

$\|M\|_{\max} = \max_{i,j} |m_{ij}|$, the maximal absolute value of all entries of $M$,

$\|M\|=\max\{\sigma_i(M^TM)\}^{1/2}$, the spectral norm (Euclidean norm for vectors).

\item For a $d_1\times d_2\times d_3$ tensor $\bT$, let ${\cal  M}} \def\cM{{\cal  M}_j(\bT)$ be a $d_j\times (d_1d_2d_3/d_j)$ matrix by unfolding $\bT$ in the $j$-th dimension.
\end{itemize}

\subsection{Signal Strengths}

Recall the low-rank decomposition of $\EE(\bA|\LL)$ in (\ref{eq:A_tucker}) with
$$
\EE(\bA|\LL)=\bar{\bC}\times_1 \bar{U}\times_2 \bar{U}\times_3 \bar{W}
$$
where $\LL=\{\ell_l\}_{l=1}^L$ and the core tensor
\begin{eqnarray*}
\bar{\bC}
&=& \bB\times_1(\bar{D}\bar{R}^{\top})\times_2 (\bar{D}\bar{R}^{\top})\times_3 D_L^{1/2} \\
&=& \bar{\bB}\times_1 \bar{D}\times_2 \bar{D}\times_3 D_L^{1/2}
\end{eqnarray*}
where
\begin{equation}\label{eq:barB}
\bar{\bB}=\bB\times_1 \bar{R}^{\top}\times_2 \bar{R}^{\top} \in\RR^{r\times r\times m}
\end{equation}
Denote the signal strengths of $\bar{\bC}$ and $\bar{\bB}$ by $\sigma_{\min}(\bar{\bC})$ and $\sigma_{\min}(\bar{\bB})$, respectively, where
\begin{equation} \label{eq:ss}
\sigma_{\min}({\bf T})=\min\big\{\sigma_{r_j}\big({\cal  M}} \def\cM{{\cal  M}_j({\bf T})\big): j=1,2,3\big\}
\end{equation}
if $\bT$ has Tucker ranks $(r_1,r_2,r_3)$.

\smallskip

We give a lower bound for the signal strength in the population adjacency tensor. The following conditions can sometimes greatly simplify our presentation.
\begin{condition}
Assume that
\begin{itemize}
     \item $(A1)$: \ $\sigma_{\min}(\bar{\bB})\geq c_1p_{\max}$, where $p_{\text{max}}=\max_{i,j,l}[\mathbb{E}A_l]_{ij}$,
     \item $(A2)$: \ $\bar{D}$ is well-conditioned, i.e. \
                         $ \sigma_1(\bar{D})\leq \kappa_0  \sigma_r(\bar{D}).$
  \item (A3): \ Minimal network balance condition:  $L_{\min} \asymp L/m$, where $L_{\min}=\min_{1\le j \le m} L_j$.
     \item (A4): \ Maximal network balance condition:  $L_{\max} \asymp L/m$, where $L_{\max}=\max_{1\le j \le m} L_j$.
\end{itemize}
\end{condition}
\begin{lemma}[signal strength] \label{lem:ss}
If conditions (A1) and (A2) hold, then we have
$$
\sigma_{\min}(\bar{\bC})\geq r^{-1}\kappa_0^{-2}\cdot m p_{\max} n \sqrt{L_{\min}}.
$$
Further if (A3) holds and $m, r, \kappa_0$ are fixed, then
$$\sigma_{\min}(\bar{\bC}) \geq C \sqrt{L} np_{\max}. $$
\end{lemma}

Condition (A1) is mild since we assumed that $\bar{\bC}$ has Tucker ranks $(r,r,m)$ implying that $\bar{\bB}$ also has the same Tucker ranks. By Lemma~\ref{lem:ss}, the signal strength of $\bar{\bC}$ is characterized by the overall network sparsity.

\smallskip

\subsection{Incoherence Property}
Theoretically, the ideal regularization parameters in Algorithm~\ref{algo:Power Iterations} are
$$
\delta_1=\max_j\|e_j^{\top}\bar{U}\|, \qquad \mbox{and} \qquad \delta_2=\max_j\|e_j^{\top}\bar{W}\|.
$$
Incoherence property ensures that singular vectors $\bar{U}$ and $\bar{W}$ are not too correlated with or incoherent to the standard basis ${e_j}$'s, as stated in the next lemma,
which the sharp convergence rates of regularized tensor power iteration algorithm rely crucially on.
\begin{lemma}[Incoherence of $\bar{U}$ and $\bar{W}$]  \label{lem:incoh}
If conditions (A1) and (A2) hold, we have
$$
\delta_1 \leq \kappa_0\sqrt{\frac{r}{n}} \quad {\rm and} \quad
\delta_2 \leq \frac{\kappa_0 r}{m\sqrt{L_{\min}}}.
$$
\end{lemma}

\smallskip

Denote $A\asymp B$ iff  $A =O(B)$ and $B=O(A)$. \
Then under conditions (A1)-(A3), it follows from Lemma \ref{lem:incoh}:
$$
\delta_2 \leq C \kappa_0r/\sqrt{mL}.
$$

\smallskip

\subsection{Tensor incoherent norms and a concentration inequality}
Given a random tensor $\mathbf{A}$, we write
$$
\bA=\EE\bA+(\bA-\EE\bA) = \mbox{the signal +  noise part}.
$$
A tensor norm is needed to measure the size of the noise. In order to deal with extremely sparse tensors, we will adopt the following definition, first introduced in \cite{yuan2017incoherent}.
\begin{definition}[Tensor incoherent norm \cite{yuan2017incoherent}] For $\delta\in(0,1]$ and $k=1,2,3$, define
$$
\Vert\mathbf{A}-\mathbb{E}\mathbf{A}\Vert_{k,\delta}:=\sup\nolimits_{\bU\in\mathcal{U}_k(\delta)}\ \langle\mathbf{A}-\mathbb{E}\mathbf{A},\bU\rangle
$$
where $\mathcal{U}_k(\delta):=\{u_1\otimes u_2\otimes u_3:\Vert u_j\Vert_{\ell_2}\le1, \forall j; \Vert u_k \Vert_{\ell_\infty}\le\delta \}$, and $\Vert u \Vert_{\ell_p}$ is the $l_p$-norm of $u$.
\end{definition}

We now present a concentration inequality for tensor incoherent norms for sparse random tensors, which is essential in proving  Theorem~\ref{thm:power_iteration} and Corollary~\ref{cor:PI} later.

\begin{theorem}[A concentration inequality for tensor incoherent norm] \label{thm:tensor_CI}
Suppose that $L\leq n$ and $Lnp_{\text{max}}\ge\log n$. Denote $n_1=n_2=n$ and $n_3=L$. Then for $k=1,2,3$, we have
\begin{equation*}
\begin{aligned}
\mathbb{P} \left\{ \Vert \mathbf{A} - \mathbb{E} \mathbf{A} \Vert_{k,\delta} \ge 3t \right\}
\le \frac{2}{n^2}+10(\log n)^2\lceil \log_2\delta^2n_k \rceil
\left[\exp\left(-\frac{t^2}{C_3p_{\text{max}}}\right)+\exp\left(-\frac{3t}{C_4\delta}\right)\right]
\end{aligned}
\end{equation*}
provided
\begin{equation*}
t \ge \max\left\{C_1 , C_2 \delta \sqrt{n_k} \log(n) \right\} \sqrt{np_{\text{max}}} \log(\delta^2n_k) \log(n).
\end{equation*}
\end{theorem}

\smallskip

We make several remarks concerning the inequality.

\begin{enumerate}
 \item  The bound in Theorem~\ref{thm:tensor_CI} is sharper than that in \cite{yuan2017incoherent}, in order to deal with extremely sparse networks. Analogous results were previously established for sparse hypergraph networks (\cite{ke2019community}), however, the dimension sizes ($n$ and $L$) in our model can be drastically different from (\cite{ke2019community}), which needs more careful treatments.

 \item   Clearly, if $\delta=1$, $\|\bA-\EE\bA\|_{k,\delta}$ reduces to the standard tensor operator norm $\|\bA-\EE\bA\|$.

From Lemma~\ref{lem:ss}, the signal strength is of the order $\sqrt{L}np_{\max}$. By using the standard operator norm to control $\bA-\EE\bA$, the success of power iterations then requires, with high probability, that
\begin{align}\label{eq:weak_cond}
\sqrt{L}np_{\max}\geq C_1'\|\bA-\EE\bA\|.
\end{align}
Under the condition $Lnp_{\max}\gg 1$, it is easy to check (by the maximum number of non-zero entries on the fibers of $\bA-\EE\bA$) that $\|\bA-\EE\bA\|\gg 1$ with high probability (see \cite[Theorem~2]{lei2019consistent}).  Consequently, by using the  standard tensor operator norm to control the size of noise part, condition (\ref{eq:weak_cond}) requires that $\sqrt{L}np_{\max}\gg 1$ rather than the optimal condition $Lnp_{\max}\gg 1$.    Indeed, if $\delta_1=O(1/\sqrt{n})$ and $\delta_2=O(1/\sqrt{L})$, Theorem~\ref{thm:tensor_CI} shows that $\|\bA-\EE\bA\|_{1,\delta_1}, \|\bA-\EE\bA\|_{3,\delta_2}=O_p(\sqrt{np_{\max}})$ up to some logarithmic factor. It can be much smaller than $\|\bA-\EE\bA\|\gg1$ (w.h.p.) when $L$ is large.

\item To apply tensor incoherent norms to analyze the convergence property of power iterations, it is necessary to prove that $\{\widehat U^{(t)}\}$ and $\{\widehat W^{(t)}\}$ are incoherent. It is possible to generalize the methods in \cite{koltchinskii2016perturbation, xia2019sup, xia2019statistical, cai2019subspace} for this purpose whose actual proof can be very involved. For simplicity, we adopt an auxiliary regularization step (\ref{eq:reg}) to truncate those singular vectors.

\end{enumerate}

\section{Main results}\label{sec:main_theory}

\subsection{Error Bound of Regularized Power Iteration}

Theorem~\ref{thm:power_iteration} states that regularized power iteration method (Algorithm~\ref{algo:Power Iterations}) works if we have a warm initialization and a strong enough signal-to-noise ratio.
These conditions are typically required (see, e.g., \cite{zhang2018tensor, xia2017statistically, ke2019community, xia2019polynomial}) and generally unavoidable (\cite{zhang2018tensor}).

For $\widehat{V}, V\in\OO_{p,r} = \{V\in\mathbb{R}^{p\times r}:V^TV=I_r \}$, the distance between their column spaces is
$$
{\rm d}(\hat{V},V):=\inf_{O\in\mathbb{O}_{r,r}}\Vert\hat{V}-VO\Vert
$$
Define
$$
{\rm Err}(t)=\max\{{\rm d}(\widehat U^{(t)}, \bar U),\ {\rm d}(\widehat W^{(t)},\bar{W}) \}.
$$
We have the following result.
\begin{theorem}[General convergence results of regularized power iterations]\label{thm:power_iteration}
Assume that
\begin{itemize}
\item the initializations $\widehat{U}^{(0)}$ and $\widehat{W}^{(0)}$ are warm, i.e., ${\rm Err}(0) \leq 1/4$,
\item the signal strength of $\bar{\bC}$ satisfies
$$
\frac{\sigma_{\min}(\bar{\bC})}{\sqrt{r\wedge 2m}}
\geq \Big(C_1 +C_2 \big((\delta_1\sqrt{n})\vee (\delta_2\sqrt{L})\big)\log n\Big)  \log\big(\delta_1^2n\vee \delta_2^2L\big) \sqrt{np_{\max}}\log(n).
$$
\end{itemize}
Then with probability at least $1-2n^{-2}$,
\begin{enumerate}
\item for all $t \le t_{\max}$, we have
$$
{\rm Err}(t) \leq \frac{1}{2}\cdot{\rm Err}({t-1})+C_3\frac{\sqrt{np_{\max}\log n}+\delta_1\delta_2\log n}{\sigma_{\min}(\bar{\bC})}
$$

\item  for $t \ge t_{\max}= C \log\big(\sigma_{\min}(\bar{\bC})/(\sqrt{np_{\max}}+\delta_1\delta_2)\big) $ iterations, we have
$$
{\rm Err}({t_{\max}}) \leq C_3\frac{\sqrt{np_{\max}\log n}+\delta_1\delta_2\log n}{\sigma_{\min}(\bar{\bC})}.
$$
\end{enumerate}
\end{theorem}

Theorem~\ref{thm:power_iteration} holds true on general tensor structures. Specializing Theorem~\ref{thm:power_iteration} to the mixture multi-layer network model (Section~\ref{sec:model}) yields the following corollary.

\begin{corollary}\label{cor:PI}
Assume that (A1)-(A3) hold. Further assume
\begin{itemize}
\item the initializations $\widehat{U}^{(0)}$ and $\widehat{W}^{(0)}$ are warm, i.e., ${\rm Err}(0) \leq 1/4$,
\item the signal strength of $\bar{\bC}$ satisfies
\begin{equation} \label{eq:sparsity_cond}
\sqrt{Lnp_{\max}} \geq \big( C_1 +C_2\kappa_0^2 (r/\sqrt{m}) \log(n) \big) r \kappa_0^2 \log(\kappa_0r) \log(n).
\end{equation}
\end{itemize}
Then with probability at least $1-2n^{-2}$, after at most $t_{\max}=O(\log(n))$ iterations,
$$
{\rm Err}({t_{\max}}) \leq C_3 \frac{\kappa_0^2 r }{\sqrt{m}} \cdot \sqrt{ \frac{\log n}{Lnp_{\max}} }.
$$
\end{corollary}
By Corollary~\ref{cor:PI}, if $\kappa_0, r, m$ are fixed and the network sparsity satisfies $Lnp_{\max} \geq C\log^4 n$, then
\begin{equation}\label{eq:rate_1}
{\rm Err}({t_{\max}}) = O_p(\sqrt{\log(n)/(Lnp_{\max}})).
\end{equation}
The convergence rate (\ref{eq:rate_1}) is optimal up to the logarithmic factor.

\subsection{Consistency of recovering global memberships}

Recall from Section~\ref{sec:model} that the global community structure is denoted as ${\cal  V}} \def\cV{{\cal  V}=\cup_{k=1}^{\bar{K}}\bar{{\cal  V}} \def\cV{{\cal  V}}_k$ with disjoint communities $\{\bar{{\cal  V}} \def\cV{{\cal  V}}_j\}_{j=1}^{\bar{K}}$ where node $i_1$ and $i_2$ belong to the same global community if and only if $(e_{i_1}-e_{i_1})^{\top}\bar{Z}=0$.
In the TWIST algorithm, after applying the K-means to the rows of $\widehat U$, we get the vertices' global membership
$\widehat{\bar \VV} = \{ \widehat{\bar{{\cal  V}} \def\cV{{\cal  V}}}_k, k\in[\bar{K}] \}$.

We measure the performance by the Hamming error of clustering:
$$
{\cal  L}} \def\cL{{\cal  L}(\widehat{\bar \VV},\bar \VV)
=\min_{\tau: \textrm{ a permutation on }[\bar{K}]} \sum_{i=1}^n {\bf 1}\big(i\in \bar{{\cal  V}} \def\cV{{\cal  V}}_k, i\notin \widehat{\bar{{\cal  V}} \def\cV{{\cal  V}}}_{\tau(k)}\big)
$$
where $\bar\VV=\{\bar{{\cal  V}} \def\cV{{\cal  V}}_k, k\in[\bar{K}]\}$.

\begin{theorem}[{\bf Consistency of global clustering}]\label{thm:global_consistency}
Assume that (A1)-(A4) hold, and that $\min_k|\bar{{\cal  V}} \def\cV{{\cal  V}}_k|\asymp n/\bar{K}$.  Then with probability at least $1-n^{-2}$, we have
$$
n^{-1} \cdot{\cal  L}} \def\cL{{\cal  L}(\widehat{\bar \VV}, \bar{\VV})
 \le C_3\kappa_0^6\frac{r^2\log n}{Lnp_{\max}}
$$
provided that the network sparsity satisfies
\begin{equation} \label{eq:global_cluster_cond1}
\sqrt{Lnp_{\max}}
\geq \big( C_1 \sqrt{\bar{K}} + C_2\kappa_0 r \log(n) \big) (\kappa_0^3 r/\sqrt m)  \log(\kappa_0r) \log(n),
\end{equation}
\end{theorem}

From Theorem~\ref{thm:global_consistency}, it follows that the relative clustering error is $O_p\big(\log^{-3}(n)\big)$ when $\bar{K}, m, \kappa_0$ are bounded. Therefore, vertices' global memberships can be consistently recovered under near optimal network sparsity conditions.

We now compare our method with some other available ones in the literature.
\begin{itemize}
\item
A special case when $m=1$ was considered by \cite{lei2019consistent}, who showed that their algorithm is able to consistently recover the communities if
$np_{\max}\sqrt{L}\gg \log^{3/2} n$. On the other hand, our result deals with more general mixture multi-layer model, is computationally more efficient, and requires weaker network sparsity: $np_{\max}L\gg \log^4(n)$. The dependence of $L$ in (\ref{eq:global_cluster_cond1}) is optimal, if we ignore the logarithmic term. This improvement is due to a sharper concentration inequality of $\bA-\EE\bA$ in terms of tensor incoherent norm.

\item
A joint matrix factorization method (Co-reg) was proposed in \cite{paul2017spectral} for a special case with $m=1$ and different $B_j$s, in which they proved that their method can consistently recover the vertices memberships if $Lnp_{\max}\gg \log n$ and the signal strengths of $B_j$s are similar. Their network sparsity condition is similar to  (\ref{eq:global_cluster_cond1}) above up to the logarithmic factor. On the other hand, our approach differs from \cite{paul2017spectral} in several aspects. Our method can perform vertices clustering and network clustering simultaneously when $m>1$.  Computationally, \cite{paul2017spectral} employed a BFGS algorithm to solve the non-convex programming, which is more computationally intensive than TWIST.
More will be discussed later in the paper after Lemma~\ref{lem:initialization}.
\end{itemize}

\subsection{Network classification}

We now show that the standard K-means algorithm on $\widehat W$ can consistently uncover the network classes of $L$ layers under the near optimal network sparsity condition (\ref{eq:global_cluster_cond1}). Further under a slightly stronger network sparsity condition (\ref{eq:network_sparsity_cond2}), we can apply Algorithm~\ref{algo:network} to exactly recover the layer labels with high probability. This shows that more layers will provide more information about layer structure and be very helpful in exact clustering of networks. Similarly, we denote
$$
{\cal  L}} \def\cL{{\cal  L}(\hat\LL, \LL)=\min_{\tau: \textrm{ permutation of }[m]}\sum_{l=1}^L {\bf 1}(\ell_l\neq \tau(\hat\ell_l)).
$$
\begin{theorem}[{\bf Consistency and exact recovery of network classes}] \label{thm:exact_recovery}
Let $\widetilde\LL=\{\tilde \ell_l\}_{l=1}^L$ be the output of the standard K-means algorithm applied to $\widehat W$.
\begin{enumerate}
\item Under the same conditions in Theorem~\ref{thm:global_consistency}, we have, with probability at least $1-n^{-2}$,
$$
L^{-1}\cdot {\cal  L}} \def\cL{{\cal  L}\big(\widetilde\LL, \LL\big)
\leq C_3\kappa_0^4\frac{(r^2/m)\log n}{Lnp_{\max}}
$$
where $\LL=\{\ell_l\}_{l=1}^L$.

\item
We further assume
\begin{equation}\label{eq:network_sparsity_cond2}
\sqrt{L}np_{\max}\geq C_1m^{-1}\kappa_0^5r^{5/2} \log(r\kappa_0) \log^{5/2}(n).
\end{equation}
Then with probability at least $1-3n^{-2}$, we have
$$
{\cal  L}} \def\cL{{\cal  L}(\widehat\LL, \LL)=0
$$
where $\widehat\LL=\{\hat{\ell}_l\}_{l=1}^L$ is the output of Algorithm~\ref{algo:network} with parameters $m$ and $\varepsilon\in[c_1, c_2] \sqrt{m/L}$.
\end{enumerate}
\end{theorem}
By Theorem~\ref{thm:exact_recovery}, in the case $r,m,\kappa_0=O(1)$, Algorithm~\ref{algo:network} is capable to exactly recover the network classes $\LL=\{\ell_l\}_{l=1}^L$ with appropriately chosen parameter $\varepsilon$ if the network sparsity (\ref{eq:network_sparsity_cond2}) satisfies $\sqrt{L}np_{\max}\gg\log^{5/2}n$.  On the other hand, consistent network clustering requires, by (\ref{eq:global_cluster_cond1}), network sparsity $Lnp_{\max}\gg \log^4n$. Therefore, condition (\ref{eq:network_sparsity_cond2}) is stronger with respect to the number of layers $L$.
\begin{remark}
Clustering of networks is essentially a problem of clustering or classifying high-dimensional data. Without exploring network structures, a simple method for clustering high-dimensional data is by spectral clustering on the left singular vectors of ${\cal  M}} \def\cM{{\cal  M}_3(\bA)$ which is a highly rectangular matrix. A simple fact is that $\|{\cal  M}} \def\cM{{\cal  M}_3(\bA-\EE\bA)\|=O_p(n\sqrt{p_{\max}})$ implying that a naive spectral clustering on ${\cal  M}} \def\cM{{\cal  M}_3(\bA)$ requires very strong condition of network sparsity for consistent network clustering. The cause of such a sub-optimality is due to the ignorance of matrix structures of rows of ${\cal  M}} \def\cM{{\cal  M}_3(\bA)$.
\end{remark}

\subsection{Consistency of local clustering}

After obtaining the network classes, it suffices to apply spectral clustering on $\sum_{l: \hat{\ell}_l=j}A_l$ for all $j\in[m]$ to recover the local memberships $\VV^j=\{{\cal  V}} \def\cV{{\cal  V}_k^j\}_{k=1}^{K_j}$. Its consistency can be directly proved by existing results in the literature (see, e.g., \cite{lei2015consistency}). We hereby omit the proof of the following theorem.
\begin{theorem}\label{thm:local_consistency}
Suppose that the conditions of Theorem~\ref{thm:exact_recovery} and (\ref{eq:network_sparsity_cond2}) hold. For all $j\in[m]$, let $\widehat\VV^j=\{\hat{\cal  V}} \def\cV{{\cal  V}_k^j\}_{k=1}^{K_j}$ denote the output of K-means algorithm on $\sum_{l: \hat{\ell}_l=j}A_l$. If $\sigma_{\min}(B_j)\geq c_1p_{\max}$ for some absolute constant $c_1>0$ and $|{\cal  V}} \def\cV{{\cal  V}_k^j|\asymp n/K_j$ for all $k\in[K_j]$, then with probability at least $1-3n^{-2}$,
$$
n^{-1}\cdot{\cal  L}} \def\cL{{\cal  L}(\widehat\VV^j, \VV^j)\leq C_1m^{-1}\kappa_0^4\cdot\frac{K_j^2\log n}{Lnp_{\max}}.
$$
\end{theorem}

\subsection{Warm initialization for regularized power iteration}\label{sec:init}
An important condition for the success of Algorithm~\ref{algo:Power Iterations} is the existence of warm initialization
$$
{\rm Err}(0) = \max\{{\rm d}(\widehat{U}^{(0)},\bar{U}), {\rm d}(\widehat{W}^{(0)},\bar{W})\}\leq 1/4.
$$
In this section, we introduce a spectral method for initializing $\widehat{U}^{(0)}$ by summing up all the layers of networks. After that, we initialize $\widehat{W}^{(0)}$ by taking the left singular vectors of ${\cal  M}} \def\cM{{\cal  M}_3(\bA)(\widetilde{U}^{(0)}\otimes\widetilde{U}^{(0)})$ where $\widetilde{U}^{(0)}={\cal  P}} \def\cP{{\cal  P}_{\delta_1}(\widehat{U}^{(0)})$. The following lemma shows that these initializations are indeed close to the truth under reasonable conditions.
\begin{lemma}[Initialization]\label{lem:initialization}
Let $\widehat{U}^{(0)}$ denote the top-$r$ left singular vectors of $\sum_{l=1}^{L}A_l$ and let $\widehat{W}^{(0)}$ be the top-$r$ left singular vectors of ${\cal  M}} \def\cM{{\cal  M}_3(\bA)(\widetilde{U}^{(0)}\otimes\widetilde{U}^{(0)})$ where $\widetilde{U}^{(0)}={\cal  P}} \def\cP{{\cal  P}_{\delta_1}(\widehat{U}^{(0)})$ with $\delta_1=\max_j\|e_j^{\top}\bar{U}\|$. Then with probability at least $1-3n^{-2}$,
\begin{equation}\label{eq:initU}
{\rm d}(\widehat{U}^{(0)},\bar{U})\leq \min\Big\{\frac{C_3\sqrt{np_{\max}}\log^2n}{\sigma_r\big(\bar{\bC}\times_3(\bfm d}   \def\bD{\bfm D}  \def\DD{\mathbb{D}_L/L)^{1/2}\big)},\ 2\Big\}
\end{equation}
where $\bfm d}   \def\bD{\bfm D}  \def\DD{\mathbb{D}_L=(L_1,\cdots,L_m)$. If $\delta_1=O(\sqrt{r/n})$ and
\begin{equation}\label{eq:init_cond}
\sigma_r\big(\bar{\bC}\times_3(\bfm d}   \def\bD{\bfm D}  \def\DD{\mathbb{D}_L/L)^{1/2}\big)\geq 4C_3\sqrt{np_{\max}}\log^2n,
\end{equation}
then with same probability,
$$
{\rm d}(\widehat{W}^{(0)},\bar{W})\leq \min\Big\{\frac{C_4\sqrt{mrnp_{\max}}\log^2(n)\log(r)}{\sigma_{\min}(\bar{\bC})},\ 2\Big\}.
$$
\end{lemma}
It is possible to improve $\log^2n$ to $\sqrt{\log n}$ in (\ref{eq:initU}) as in \cite{paul2016consistent}.
Comparing the rate of initialization (\ref{eq:initU}) and the rate after regularized power iterations in Theorem~\ref{thm:power_iteration}, Algorithm~\ref{algo:Power Iterations} improve the estimation error by a ratio of $\sigma_{\min}(\bar{\bC})$ and $\sigma_r\big(\bar{\bC}\times_3(\bfm d}   \def\bD{\bfm D}  \def\DD{\mathbb{D}_L/L)^{1/2}\big)$. In special cases, such an improvement can be significant. For instance, consider $m=r=2$ and $L_1=L_2$ and $\bar{\bC}\in\RR^{2\times 2\times 2}$ with
$$
\bar{C}(:,:,1)=\left(\begin{array}{cc}1+\varepsilon&0\\0&1+\varepsilon\end{array}\right)\quad {\rm and}\quad \bar{C}(:,:,2)=\left(\begin{array}{cc}0&1-\varepsilon\\1-\varepsilon&0\end{array}\right)
$$
for some small number $\varepsilon\in(0,1)$. It is easy to check that $\sigma_{\min}(\bar{\bC})=\sigma_2\big({\cal  M}} \def\cM{{\cal  M}_3(\bar{\bC})\big)=\sqrt{2}(1-\varepsilon)$. On the other hand,
$$
\sigma_{2}\big(\bar{\bC}\times_3(\bfm d}   \def\bD{\bfm D}  \def\DD{\mathbb{D}_L/L)^{1/2}\big)=\sqrt{2}\varepsilon.
$$
Moreover, if $\varepsilon=0$, then $\bar{\bC}\times_3(\bfm d}   \def\bD{\bfm D}  \def\DD{\mathbb{D}_L/L)^{1/2}$ is rank deficient implying that simply projecting the multi-layer networks into a graph can potentially cause serious  information loss. See more details in \cite{ke2019community} and a similar discussion in \cite{lei2019consistent}.

It is worthwhile pointing out that one could use other methods to initialize $\widehat{U}^{(0)}$ (the initialization of $\widehat{W}^{(0)}$ is easy once it is done for $\widehat{U}^{(0)}$). Examples include the HOSVD by extracting the top-$r$ left singular vectors of ${\cal  M}} \def\cM{{\cal  M}_1(\bA)$, the joint matrix factorization method in \cite{paul2017spectral}, and random projection \cite{ke2019community}.

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

We conduct several simulations to test the performance of the TWIST on the MMSBM with different choices of network sparsity, "out-in" ratio, number of layers and the size of each layer. We use K-means as the clustering algorithm.
The evaluation criterion is the mis-clustering rate. All the experiments are replicated 100 times and the average performance across the repetitions is reported.

We generate the data according to the MMSBM in Section \ref{sec:model} in the following fashion. The underlying class $\ell_l$ for each layer $l$ is generated from the multinomial distribution with $\PP(\ell_l=j)=1/m,\ j=1,\cdots,m.$  The membership $z_i^j$ for each node $i$ in layer type $j$ is generated from the multinomial distribution with $\PP(z_i^j=s)=1/K,\ s=1,\cdots,K.$ We choose the connection matrix as
$B=pI+q(\mathbf{1}\mathbf{1}^T-I),$ where $\mathbf{1}$ is a $K$-dimensional vector with all elements being $1.$ Let $\alpha=q/p$ be the out-in ratio.

\subsection{Global memberships}
First, we consider the task to detect the global memberships defined in Section~\ref{sec:model}. We compare the performance of the TWIST with Tucker decomposition with HOSVD initialization (HOSVD-Tucker), and  we also adopt a baseline method by performing spectral clustering on the sum of the adjacency matrices from all layers (Sum-Adj). Sum-Adj has been considered in literature (\cite{paul2017spectral}, \cite{dong2012clustering}, \cite{tang2009clustering}) as a simple but effective procedure \cite{kumar2011co}. The function "\textit{tucker}" from the R package "\textit{rTensor}" \cite{li2018rtensor} is adopted to apply Tucker decomposition for HOSVD-Tucker.

In Simulation 1, the networks are generated with the number of node $n=600$, the number of layers $L=20$, number of types of networks $m=3,$ number of communities of each network $K=2$ and out-in ratio of each layer $\alpha=0.4.$ The average degree $d$ of each layer varies from $2$ to $20.$

In Simulation 2, the setting is the same as in Simulation 1 except the average degree of each layer is fixed at $d=10$ and the out-in ration $\alpha$ of each layer varies from $0.1$ to $0.8.$

In Simulation 3,  the setting is the same as in Simulation 1, except that the out-in ratio $\alpha=0.6$ and the number of layers $L$ varies from $10$ to $60.$

In Simulation 4, the setting is the same as that in in Simulation 3, except that the out-in ratio is $L=20$ and the number of nodes $n$ varies from $100$ to $1200.$

\begin{figure}
	\centering
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/global_membership_d}
		\caption{\tiny{The result of Simulation 1: $n=600$, $K=2$, $m=3$, $L=20$, $\alpha=0.4$, varying $d.$ }}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/global_membership_r}
		\caption{\tiny{The result of Simulation 2: $n=600$, $K=2$, $m=3$, $L=20$, $d=10$, varying $\alpha.$} }
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/global_membership_L}
		\caption{\tiny{The result of Simulation 3: $n=600$, $K=2$, $m=3$, $d=10$, $\alpha=0.6$, varying $L.$ }}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
	\includegraphics[width=\linewidth]{figs/global_membership_n}
	\caption{\tiny{The result of Simulation 3: $K=2$, $m=3$, $d=10$, $L=20$, $\alpha=0.6$, varying $n.$ }}
    \end{subfigure}
	\caption{Overall, the TWIST and the HOSVD-Tucker perform much better than Sum-Adj. The TWIST outperforms the HOSVD-Tucker
	 when the signal is not strong enough, for instance $d<6$ in (a), $\alpha>0.5$ in (b), $L<50$ in (c) and $n<800$ in (d).}
	\label{fig:global membership}
\end{figure}

The results of Simulations 1-4 are given in Figure \ref{fig:global membership}.
\begin{enumerate}
\item Clearly, the mis-clustering rate of all the methods decreases as the average degree of each layer increases, the out-in ratio of each layer decreases and the number of layers increases. This is consistent with our theoretical findings.

\item TWIST and HOSVD-Tucker, both utilizing tensor structure, perform much better than Sum-Adj, which only uses the matrix structure.
The mis-clustering rate of the TWIST and the HOSVD-Tucker decreases more rapidly.

\item TWIST outperforms HOSVD-Tucker when the signal is not strong enough, e.g., for $d<6$ in Simulation 1; for $\alpha>0.5$ in Simulation 2; for $L<50$ in Simulation 3; and for $n<800$ in Simulation 4.
\end{enumerate}

\subsection{Layers' labels}
We now explore the task of clustering different types of layers. We compare the TWIST with the HOSVD-Tucker and spectral clustering applied to the mode-3 flatting of \bA (M3-SC).

In Simulation 5, the networks are generated with the number of node $n=600$, the number of layers $L=20$, number of types of networks $m=3,$ number of communities of each network $K=3$ and out-in ratio of each layer $\alpha=0.6.$ The average degree $d$ of each layer varies from $3$ to $30.$

In Simulation 6, the networks are generated as in Simulation 5, except that the average degree of each layer $d=10,$ the number of layers $L=30$ and the out-in ration $r$ of each layer varies from $0.1$ to $0.9.$

In Simulation 7, the networks are the same as in Simulation 6, except that the out-in ratio $\alpha=0.6$ and the number of layers $L$ varies from $20$ to $80.$

In Simulation 8, the networks are the same as in Simulation 7, except that the average degree of each layer $d=0.02n$ and the the size of each layer $n$ varies from $100$ to $1200.$

\begin{figure}
	\centering
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/network_membership_d}
		\caption{\tiny{The result of Simulation 4: $n=600$, $K=3$, $m=3$, $L=20$, $\alpha=0.6$, varying $d.$ }}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/network_membership_r}
		\caption{\tiny{The result of Simulation 5: $n=600$, $K=3$, $m=3$, $d=10$, $L=30$, varying $\alpha.$ }}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/network_membership_L}
		\caption{\tiny{The result of Simulation 6: $n=600$, $K=3$, $m=3$, $d=10$, $\alpha=0.6$, varying $L.$ }}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
	\includegraphics[width=\linewidth]{figs/network_membership_n}
		\caption{\tiny{The result of Simulation 7: $K=3$, $m=3$, $d=0.02n$, $\alpha=0.6$, $L=30,$ varying $n.$ }}
    \end{subfigure}
	\caption{The TWIST is the best overall, particularly when the signal is not strong enough, for instance, $d<15$ in (a), $r>0.4$ in (b), $L<50$ in (c) and $n<800$ in (d). From Simulation 7, the naive method M3-SC hardly changes as the number of layers increases.
	}
	\label{fig:network membership}
\end{figure}

The results are presented in Figure \ref{fig:network membership}. We make the following observations.
\begin{enumerate}
\item The mis-clustering rates of all three methods decrease as the average degree of each layer increases, the out-in ratio of each layer decreases, the number of layers increases and the size of each layer increases. This agrees with our theoretical results.

\item From Simulation 7, the naive method M3-SC shows no response to the increase of the number of layers, as might be expected.

\item Overall, the TWIST performs the best among the three methods. This can be clearly seen when the signal is not strong enough, for instance, $d<15$ in Simulation 5, $\alpha>0.4$ in Simulation 6, $L<50$ in Simulation 7 and $n<800$ in Simulation 8.
\end{enumerate}

\section{Real data analysis}\label{sec:realdata}

In this section, we apply the TWIST to two real data sets: worldwide food trading networks and Malaria parasite genes networks. The two datasets have been studied in the literature before. However, with the TWIST, we are able to make some new, interesting, and sometimes surprising findings, which the earlier methods have failed to do so.

\subsection{Malaria parasite genes networks}

The var genes of the human malaria parasite Plasmodium falciparum present a challenge to population geneticists due to their extreme diversity, which is generated by high rates of recombination. Var gene sequences are characterized by pronounced mosaicism, precluding the use of traditional phylogenetic tools. \cite{larremore2013network} identifies 9 highly variable regions (HVRs), and then maps each HVR to a complex network, see Figure \ref{fig:HVR}. They showed that the recombinational constraints of some HVRs are correlated, while others are independent, suggesting that this micromodular structuring facilitates independent evolutionary trajectories of neighboring mosaic regions, allowing the parasite to retain protein function while generating enormous sequence diversity.

\begin{figure}
	\centering
	\begin{subfigure}[b]{.6\linewidth}
		\includegraphics[width=\linewidth]{figs/HVR}
		\caption{Extracted highly variable regions.}
	\end{subfigure}
\\
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_1}
		\caption{HVR1.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_2}
		\caption{HVR2.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_3}
		\caption{HVR3.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_4}
		\caption{HVR4.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_5}
		\caption{HVR5.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_6}
		\caption{HVR6.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_7}
		\caption{HVR7.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_8}
		\caption{HVR8.}
	\end{subfigure}
	\begin{subfigure}[b]{.18\linewidth}
		\includegraphics[width=\linewidth]{figs/Net_HVR_9}
		\caption{HVR9.}
	\end{subfigure}
\\
	\begin{subfigure}[b]{.6\linewidth}
		\includegraphics[width=\linewidth]{figs/Genes_net_9}
		\caption{The $212\times 212\times 9$ mixture multi-layer tensor.}
	\end{subfigure}
	\caption{9 extracted highly variable regions and corresponding networks. HVR2, HVR3 and HVR4 are quite sparse and not mainly connected with only some small components.}
	\label{fig:HVR}
\end{figure}

Despite the innovative network approach, there are still some drawbacks in \cite{larremore2013network}.
\begin{enumerate}
	\item Even though 9 HVRs have been identified, only 6 HVRs have been used in the analysis, while the other three HRVs are discarded due to their sparse structures, as seen in Figure \ref{fig:HVR} (c,d,e). However, these sparse networks still contain valuable information, which would be of great interest to researchers and practitioners.
	\item Community structures are identified individually for each network, and then compared with each other to identify similar structures. This is not only very demanding and tedious computationally, but also involves much human intervention. This becomes increasing undesirable as the number of networks grows bigger.
\end{enumerate}

Here we propose to employ the TWIST to the problem, in order to overcome the above difficulties.
The data under investigation are the 9 highly variable regions (HVRs) used in \cite{larremore2013network}. Each network is  derived from the same set of 307 genetic sequences from var genes of malaria parasites. A node represents a specific gene and an edge is generated by comparing sequences pair-wisely within each HVR. More information about the data and data pre-processing could be found in \cite{larremore2013network}.
In our study, we consider 212 nodes which appear on all 9 layers. This results in a $212\times 212\times 9$ mixture multi-layer tensor, as shown in Figure \ref{fig:HVR} (k).

We apply the TWIST to this $212\times 212\times 9$ tensor with the core tensor $15\times 15\times 3.$ The embedding of each layer is plotted in Figure \ref{fig:genes_layer_embedding}. We make the following comments.
\begin{enumerate}
	\item The 9 HVRs fall into 4 groups (Figure \ref{fig:genes_layer_embedding} (a)):
	$
	\{1, 2, 3, 4, 5, 6\}, \{7\}, \{8\}, \{9\}.
	$

	\smallskip

	By comparison, \cite{larremore2013network} found that the 6 HVRs fall into 4 groups (without layers 2-4):
	$\{1, 5, 6\}, \{7\}, \{8\}, \{9\}.$  The two findings are consistent.

	\smallskip

	\item  TWIST places sparse networks of layers 2-4 to the same group as layers 1 and 5.

	\smallskip

	By comparison, the sparse layers 2-4 had to be discarded in \cite{larremore2013network}. The new result implies that the sequences remains mostly unchanged in the beginning (HVRs 1-6), and start to diversity from HVR 7 onward.

	\smallskip

	\item Hierarchical structure of the 9 HVRs.

	\smallskip

	If we zoom in the mini group $\{1, 2, 3, 4, 5, 6\}$ (Figure \ref{fig:genes_layer_embedding}(b)), we notice that the first 5 layers are more tied together, so we have a finer partition: $\{1, 2, 3, 4, 5\}, \{6\}$. This operation can be repeated. Therefore, TWIST can be easily used to form a hierarchical structure of the 9 HVRs (Figure \ref{fig:genes_layer_embedding} (c)).

	\smallskip

	\item Computational ease of the TWIST

	\smallskip

	TWIST can easily cluster layers and nodes using $K$-means. This is much easier than the procedure in \cite{larremore2013network}, which first finds the community structure for each layer, and then computes their similarities.

	\smallskip

	\item Better community structure is obtained by combining information from similar layers.

	\smallskip

	TWIST is applied to the first 6 similar layers $\{1, 2, 3, 4, 5, 6\}$ to identify their common local structure, while spectral clustering is applied to HVR 6 to find its community structure, as was done in \cite{larremore2013network}, see Figure \ref{fig:genes_embedding} (a)-(b).
	Clearly, the 4 local communities are much more separated in Figure \ref{fig:genes_embedding}(a) than in (b).

\end{enumerate}

\begin{figure}
	\centering
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/network_embedding_9}
		\caption{Embedding of each network.}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
		\includegraphics[width=\linewidth]{figs/network_embedding_6}
		\caption{Zoom center area in (a).}
	\end{subfigure}
	\begin{subfigure}[b]{.4\linewidth}
	\includegraphics[width=\linewidth]{figs/hierarchical}
	\caption{Hierarchical structure of network classes.}
\end{subfigure}
	\caption{Embedding of each layer in malaria parasite genes networks and hierarchical structure of network classes.} \label{fig:genes_layer_embedding}
\end{figure}

\begin{figure}
	\centering
	\begin{subfigure}[b]{.50\linewidth}
		\includegraphics[width=\linewidth]{figs/genes_tucker_community_123456}
		\caption{TWIST applied to HVRs 1-6}
	\end{subfigure}
	\begin{subfigure}[b]{.44\linewidth}
		\includegraphics[width=\linewidth]{figs/genes_spectral_community_6}
		\caption{Spectral clustering applied to HVR 6}
	\end{subfigure}
	\caption{Nodes embedding using TWIST and spectral decompositions with $K=4$.}
	\label{fig:genes_embedding}
\end{figure}

\subsection{Worldwide food trading networks}

We consider the dataset on the worldwide food trading networks, which is collected by \cite{de2015structural}, and is available at \url{http://www.fao.org}. The data contains an economic network in which layers represent different products, nodes are countries and edges at each layer represent trading relationships of a specific food product among countries.

We focus on the trading data in 2010 only. We convert the original directed networks to undirected ones by ignoring the directions. We delete the links with weight less than 8 (the first quartile) and abandon the layers whose largest component consists less than 150 nodes. These are done to to filter out the less important information. Finally we extract the intersections of the largest components of the remaining layers.

After data preprocessing, we obtained a 30-layers network with 99 nodes at each layer. Each layer represents trading relationships between 99 countries/regions worldwide with respect to 30 different food products. Together they form a mixture multi-layer tensor of dimension $99\times99\times30$.

We first apply Algorithm 1 in the TWIST procedure to the mixture multi-layer tensor, which results in a tensor decomposition with a core tensor of dimension $20\times 20\times 2$. The resulting two clusters of layers are listed in Table \ref{food_cluster}. We then apply Algorithm 2 in the TWIST procedure to each cluster separately (here we have two clusters) to find the community structures for each cluster, in order to obtain the clustering result of countries. This time, we take the core tensor of dimension $4\times 4\times 1.$ The embedding of 99 countries with clustering results from K-means are shown in Fig. \ref{fig:country_embedding}.
For the two types of networks, we plot the sum of adjacency matrices with nodes arranged according the community labels in Figure \ref{fig:food_networks} to have a glance of different community structures of two network types.

\begin{table}
\begin{center}
	\begin{tabular}{ | l  p{10cm} |}
		\hline
		Food cluster 1: &
		Beverages non alcoholic, Food prep nes,       Chocolate products nes, Crude materials, Fruit prepared nes, Beverages distilled alcoholic, Coffee green, Pastry, Sugar confectionery, Wine, Tobacco unmanufactured    \\ \hline
		Food cluster 2: &
	    Cheese whole cow milk, Cigarettes, Flour wheat, Beer of barley, Cereals breakfast, Milk skimmed dried, Juice fruit nes, Maize, Macaroni, Oil palm, Milk whole dried, Oil essential nes, Rice milled, Sugar refined, Tea, Spices nes, Vegetables preserved nes, Waters ice etc, Vegetables fresh nes  \\
		\hline
	\end{tabular}
\caption{List of two clusters of food.}
\label{food_cluster}
\end{center}
\end{table}

\begin{figure}
	\centering
	\begin{subfigure}[b]{.7\linewidth}
		\includegraphics[width=\linewidth]{figs/food_clustering1_country}
		\caption{Embedding of countries for networks in cluster 1}
	\end{subfigure}
	\begin{subfigure}[b]{.7\linewidth}
		\includegraphics[width=\linewidth]{figs/food_clustering2_country}
		\caption{Embedding of countries for networks  in cluster 2}
	\end{subfigure}
	\caption{Embedding of countries on two different types of food trading networks.}
	\label{fig:country_embedding}
\end{figure}

\begin{figure}
	\centering
	\begin{subfigure}[b]{.46\linewidth}
		\includegraphics[width=\linewidth]{figs/food_clustering1_heatmap}
		\caption{Heat map of networks in cluster 1.}
	\end{subfigure}
	\begin{subfigure}[b]{.46\linewidth}
		\includegraphics[width=\linewidth]{figs/food_clustering2_heatmap}
		\caption{Heat map of networks in cluster 2.}
    \end{subfigure}
	\caption{Heat maps of two types of networks.}
	\label{fig:food_networks}
\end{figure}

We make the following remarks from Table \ref{food_cluster}, Figures \ref{fig:country_embedding} and  \ref{fig:food_networks}.
\begin{enumerate}
\item
Trading patterns of food are different for unprocessed and processed foods.

\smallskip

Specifically, cluster 1 consists mainly of raw or unprocessed food (e.g., crude materials, coffee green, unmanufactured tobacco), while cluster 2 is mainly made of processed food (e.g., such as cigarettes, flour wheat, essential oil, milled rice, refined sugar).

\smallskip

\item  For unprocessed food, global trading is the more dominant trading pattern than regional one. Some countries have closer trading ties with countries across the globe.

\smallskip

From cluster 1, a small number of countries, such as China, Canada, United Kingdom, United States, France, Germany, are very active in trading with others as well as amongst themselves. This small group of countries is called a hub community.
This implies reflects the fact that these large countries import unprocessed food from, and/or export unprocessed food to a great number of other countries worldwide.

\smallskip

\item For processed foods, regional trading is very dominant. In fact, the world trading map is striking similar to the world geography map in Figure \ref{fig:country_embedding} (b).

\smallskip

I cluster 2, countries are mainly clustered by the geographical location, i.e.,  countries in the same continent have closer trading ties. Examples of these clusters include countries in America (United State, Canada, Mexico, Brazil, Chile), in Asia and Africa (China, Japan, Singapore, Thailand, Indonesia, Philippines, India), and in Europe (Germany, Italy, Poland, Spain, Denmark, Switzerland). Regional trading of processed food can have many advantages, e.g., keeping the food cost low due to lower transportation cost, and keeping food refresh due to faster delivery.

\smallskip

There are some interesting "outliers" as well. For instance, United Kingdom has closer trading ties with African and Middle Eastern countries than its European neighbor, which might be interesting to delve into further.

\end{enumerate}

\section{Conclusion and discussion}\label{sec:discussion}

In this paper, we have proposed a novel mixture multi-layer stochastic block model (MMSBM) to capture the intrinsic local as well as global community structures. A tensor-based algorithm, TWIST, was proposed to conduct community detection on multi-layer networks and shown to have near optimal error bounds under weak conditions in the MMSBM framework. In particular, the method allows for very sparse networks in many layers. The proposed method outperforms other state of the art methods both in nodes community detection and layers clustering by extensive simulation studies. We also applied the algorithm to two real dataset and found some interesting results.

A number of future directions are worth exploring. As a natural extension, one can generalize the tensor-based representation to account for adjacency matrices capturing the degree heterogeneity of nodes.  The layers of networks could have the spacial and temporal structures of networks in many real applications, one could incorporate these into the model. On a more theoretical level, it is of interest to explore  theoretical properties in more sparse scenario. It is also important to develop scalable algorithms which can handle millions nodes with thousands layers in this big data era.

\bibliographystyle{plainnat}

redpajama_set_name: RedPajamaArXiv


text:
Q: boost hash returning same value for different inputs I have two objects, Account and Transaction where Transaction is the unique pair of Account and an incrementing id number.  I want to use boost::hash to get unique values for these and have overloaded the hash_value method per the instructions: http://www.boost.org/doc/libs/1_53_0/doc/html/hash/custom.html
class Account {
  ...
};

class Transaction
{
    Account account;
    unsigned int id;
};

Account's hash_value method works correctly, and the value returned is always unique for a given account, however to make the unique pair, Transaction's method needs to use hash
_combine (per boost's instructions ):
inline std::size_t hash_value( const Account& acct )
{
    boost::hash<int> hasher;
    size_t rval = hasher( acct.id() ); //just an int. guaranteed to be unique
    return rval;
}

inline std::size_t hash_value( const Transaction& t )
{
    std::size_t seed = 0;
    boost::hash_combine( seed, t.account );
    boost::hash_combine( seed, t.id );

    return seed;
}

This returns the same values for different inputs sometimes.  Why??  I only have a few thousand accounts, and the id number only goes up to a few hundred thousand.  this doesn't seem like an upper bound issue.
Does anyone know if this is a bug, or if I need to seed boost hash?
Thanks

A: Look up perfect hashing, and the birthday paradox, and for completeness's sake the pigeonhole principle.
What it boils down to is hash functions generally do produce collisions,unless what you're hashing has very specific properties you've taken advantage of.  Your chances of seeing a hash collision for any given set of keys is going to be counterintuitively high just because that's one of the mathematical realities we're not wired for: with a 1/365 chance of getting any particular hash, your odds of a collision are 50/50 given just 23 keys.

A: Boost provides good generic hash functions because it makes no/few assumptions about the input and tries to be fast. In most cases you can make specific assumptions about the input to create a far better hash function than what you get from boost. For example you can optimize a string hash function by assuming the string contains english text. By using assumptions you can make far better hash functions (as in: far less collisions). For example if you need to merge two hash values that are each integers between 1 and 1000 it's obvious that you will not get collisions is you multiply one of them by 1000 and then add the other.
Be very careful when writing custom hash functions because there is a clear disadvantage beyond getting it wrong: Code robustness always suffers.
Example 1: You optimize a UTF-8 string hash for english language strings. Suddenly the application gets chinese language strings.
Example 2: You assume an ID is always small because the IDs start at 1, increase by one each time one is assigned and there are never more than a few thousand assigned. Now someone changes the id to be a random GUID.

redpajama_set_name: RedPajamaStackExchange


As a Collaborative Divorce Coach I work as a neutral member of Collaborative Teams to help couples in the Collaborative divorce process to have a healthier divorce for themselves and their children. I believe that families don't belong in court. I work as a peacemaker to coach the parties to improve communication and co-parenting skills and help to manage the process of respectful resolution to conflict. If you are considering divorce or separation, consider the Collaborative process to avoid the adversarial approach. I also work as a trainer to provide basic and advanced trainings in the Collaborative Process.
Comments: I help to resolve disputes collaboratively. My passion is to work as a peacemaker for families. I help to bring respect and peace to families in crisis in order to resolve family problems and avoid court. I believe that the Collaborative Practice offers a healthy alternative to the traditional method of resolving conflict. I am always excited to help to facilitate a healthy process of resolving disputes for the benefit of all parties, including children. Healthy communication helps to move all parties to a place where interests can be met and hopes and fears can be addressed.
October 25, 2014: A Practical Approach to Handling Coercive Control or Violence Issues in CP, Lori Gephart, Lynn Hutchinson, Patricia Rogers, Lisa Alexander, Mary Pence, 1.5 hours.
November 3, 2014: Strategies for streamlining the Collaborative Process, Derek Smith, Paula Hopkins, Lori Gephart, 1 hour.
November 3, 2014: How to set up and Manage a Profitable Collaborative Practice - Panel Mark Byelich (Bucks), Peter Bort (Philadelphia), Mary Louise Parker, Diane Hitzemann (East Central), Lori Gephart (Pgh), 1 hour.
February 6, 2015: Domestic Violence and Collaborative Practice 1 hour.
April 23-25, 2015: "Three Day Interdisciplinary Basic Training in Collaborative Practice" by Lea Anderson, Marc Criss, Loretta Gephart, M.A., Patricia Henk and Paula Hopkins. 18 hours.
November 6, 2015: ACBA Collaborative Law Committee Panel Discussion mental health issues in Collaborative practice by Loretta Gephart, M.A. and Patricia Boswell 1 hour.
Barbara A. Burr, JD; Lisa Herrick, PhD; Kate Sharff, MSW, LCSW-C, LICSW; Marshall Yoder, MA, JD.
Linda Solomon, LMFT and Barbara Kelly, PhD.
Jane Iannuzzelli, M.A., Eve Orlow, M.S., Steven Cohen, Ph.D.
Jane Iannuzzelli, M.Ed., M.A., Eve Orlow, Ed.D., M.S., Steven Cohen, Ph.D.
Karen Levitt, JD; Lynda Robbins, JD and Sanford Portnoy, Ph.D.


text:
Poppy is a FABULOUS example of what Bridge can achieve with the determination and drive.
Poppy, 24 from Hexham joined Bridge Northumberland for support in Feb 2021 being referred after a visit to her GP and speaking with the Social Prescriber at the surgery.
Poppy left university in 2018 after completing a HR Business degree, and found work in the family business until it was sold in May 2019. Becoming unemployed triggered past anxiety and mental wellbeing issues for Poppy which were further impacted by the pandemic and the numerous lockdowns. Poppy became increasingly agoraphobic and even the simplest of day to day things were impossible for her to do, she never left her home. The situation continued for months until in Feb 2021 Poppy decided enough was enough and she went to her GP to ask for help.
The GP referred her to the clinics Social Prescriber Alex, who was able to give Poppy time to explain her difficulties and work with her in finding different solutions to support her. Alex suggested Poppy join Bridge Northumberland to assist her in her journey to finding work.
On joining Poppy met with Sheena her Bridge Worker. It was evident to Sheena that Poppy had skills and abilities and a degree but she really needed to rebuild her confidence before looking for work. After the financial assessment was completed with Northumberland Citizens Advice, Sheena started with arranging coaching sessions weekly with Cygnus Support.
Poppy says the coaching sessions with Cygnus were different to mental health support she had had previously, Cygnus really listened and focused it all about you. The regular sessions helped Poppy rebuild her confidence to the point she felt able to start a volunteering role when Sheena suggested it in May.
Working with Lisa at Northumberland CVA Poppy started volunteering weekly with Full Circle Food Project, in a admin role. A job and tasks she was familiar with and able to do. The weekly sessions really helped rebuild Poppy's confidence and then she progressed into volunteering at Children's cooking classes, at Full Circle.
In July Sheena suggested Poppy start the Skills U Need course run by Northumberland CVA, the course covers writing the perfect CV, carrying out effective job searches, complete application forms and learn how to be successful at interview. The training also helps participants prepare for work and manage on a budget. Poppy was happy to attend the weekly course to support herself in searching for work.
In August at one of the weekly sessions at Skills U Need, Michelle the course tutor told Poppy of a job opportunity that was being advertised at Northumberland CVA. Michelle suggested she apply for the job as she had all the skills and experience the role required. Poppy put in an application and soon heard she had an interview in early September. The interview went well and Poppy was delighted when she was offered the role of Receptionist & Admin Assistant at Northumberland CVA. Poppy accepted and started her new job in Oct.
Sarah at Full Circle Food Project says: "Poppy was quite nervous coming to us initially as an Admin. Volunteer, but she had all of the intelligence and skills that helped us. She gradually developed confidence that meant that she not only helped in the office, but also helped with children's cooking classes as well. We are so pleased that she has got the job she deserves, but we will certainly miss her! "
The decisive action taken back in February 2021 to where she is now 8 months later, Poppy has succeeded by being driven and determined to improve her life and make lasting changes. She has taken all the support and help on offer to enable her to reach her goal.
Poppy says: "You can't change until you do it yourself!"
Sheena says of Poppy's journey: "Poppy is a FABULOUS example of what Bridge can achieve with the determination and drive. Poppy was lacking confidence and felt far-removed from the job market, but with the support of Bridge and The Full Circle Food Project, together with her hard work and motivation, she has succeeded in securing a great new role and we are all so proud."

redpajama_set_name: RedPajamaCommonCrawl


text:
Among the Karo people of North Sumatra, Indonesia, the oukup is a traditional form of sauna in which the bather is wrapped in blankets with spice-scented steam coming from a boiling cauldron. The modern oukup is now done using steam piped into a small room. A similar bath, known as mar-tup was practised by the Batak Toba people. Both tribes traditionally performed the bath after childbirth, but it is now used as a general restorative or for relaxation, and oukup parlours typically also offer massage services.

Many modern oukup parlours can be found in Medan and other cities in North Sumatra, since the opening of the first modern oukup parlour in 1992 in Medan. The bath is scented using fresh ginger, turmeric, pepper, star anise, temu kunci, coriander, celery leaves, kencur, and many other spices and herbs. Oukup parlours can now also be found in other parts of Indonesia, such as Jakarta.

References

Batak Karo
Bathing
Sauna

redpajama_set_name: RedPajamaWikipedia


Deanna deBara is a freelance writer living in Portland, OR. Free from alcohol for nearly eight years, Deanna writes regularly about addiction, recovery, and mental health. You can learn more about her at deannadebara.com.


Here at CSSI, the BackOffice Team serves our clients by providing a wide range of services. Whether your firm is looking to outsource operations due to unexpected employee turnover, internal staff transitions, or just trying to reduce cost CSSI can help. Through secure encrypted connectivity our BackOffice team is able to provide temporary, specialized, or permanent services remotely.
What Do CSSI's BackOffice Services include?
Don't see something your firm needs? Just ask! CSSI prides ourselves in our ability to meet client's ever evolving needs.
to learn more about how CSSI can assist your firm with any and all of your BackOffice needs.
Advent has been providing portfolio management and accounting, straight through processing and trade order management software, services and data to the world's leading financial professionals. CSSI holds a current Advent Axys® license.
Founded in 1999, ByAllAccounts delivers reconciliation-ready data to investment managers and advisory firms. Using innovative expert system technology, ByAllAccounts provides unparalleled access to financial account data from thousands of custodians and enables efficient and timely account reconciliation, performance analysis and consolidated reporting.
GIPS® is a trademark owned by CFA Institute.
Portfolio Accounting is the backbone of an investment firms operation and reporting process. With CyGems outsourced reconciliation service you can concentrate on the big picture of servicing your customers while CyGems professional team manages the process.
Tamarac was founded in 2000 as a Software as a Service (SaaS) company with an emphasis on service. Tamarac's focus on innovation and best-in-breed technology allows us to incorporate user feedback and add significant enhancements to our platform; ongoing field testing keeps us ahead of an ever-changing industry. These enhancements help ensure that your firm is always using the most cutting edge reporting, portfolio management, and client engagement technology.
Tamarac is an open platform with a large and growing list of supported integration partners — providing your firm with the flexibility and support you need to automate more of your business.


Part 1 of 2. The efforts of rural Hartsville, S.C., to provide a quality education to children are chronicled in this documentary, which spotlights parents, students, teachers and administrators at two elementary schools. Many students come from families at or near the poverty line, with more than half qualifying for free or reduced-price lunches. Yet, the school district has a 92-percent graduation rate, which is higher than the norm for low-income students.
Conclusion. The efforts of rural Hartsville, S.C., to provide a quality education to children are chronicled in this documentary, which spotlights parents, students, teachers and administrators at two elementary schools. Many students come from families at or near the poverty line, with more than half qualifying for free or reduced-price lunches. Yet, the school district has a 92-percent graduation rate, which is higher than the norm for low-income students.
There are no featured audience reviews yet. Click the link below to see what others say about 180 Days: Hartsville: Season 1!


text:
Q: Laravel - LIKE operator to search for encrypted values I am trying to implement a search module in my Laravel app that could filter users by name. In order to protect users, the 'name' column rows are encrypted on DB. The problem that I am facing is that query below always returns 0 results.
I am encrypting the search input value before searching into DB.
$patients = DB::select( DB::raw("SELECT * FROM patient WHERE name LIKE '%".Crypt::encrypt($searchText)."%';"));

What am I doing wrong here ?

A: Crypt::encrypt("Text")

The above will rarely generate the same value each time you call it, therefore you will not have the same value stored in your DB.
Suggestion
Don't encrypt users names unless you really have to. If you do have to encrypt the names then try the following.
$patients = DB::select("SELECT * FROM patient WHERE name = '".Crypt::encrypt($searchText)."'");

Your encrypted value will be case sensitive so change the column collation to latin1_general_cs

A:
What am I doing wrong here ?

Laravel's encryption is properly implemented, which means that each message is randomized!
Randomization introduces a challenge for implementing searchable encryption: You never get the same ciphertext twice. Fortunately, there's a workaround for making encryption searchable.

*

*Encrypt each record using authenticated encryption.

*Separately, store HMAC-SHA2/BLAKE2 of the plaintext (with a separate key than the encryption key). You may truncate the hash function output. This is called a blind index.

*

*Additionally, you may store multiple hash functions of various slices or transformations of the plaintext. (e.g. last 4 digits, first initial + last name, etc.)

*Each blind index should have a distinct key.

*You can combine multiple fields into a single blind index.

*For low-entropy inputs (e.g. boolean fields), never use them in isolation.

*You may also use a slow KDF function (PBKDF2, Argon2id) for highly sensitive fields, although this is not strictly necessary.

*Recalculate the blind indexes when doing SELECT queries, rather than trying to select based on the ciphertext itself.

This may sound like a lot, but in practice, CipherSweet implements all the low-level machinery and should be easy enough for most PHP developers to use. (If this turns out to not be the case, file a bug on Github; we want it to be the case!)