\documentclass[11pt, letterpaper]{article}
\usepackage{epsfig,mst-stylefile,amssymb,amsmath, multirow, url, tcolorbox,booktabs, enumitem}

% For bibtex
%\renewcommand{\cite}{\citet*}
%\bibliographystyle{agsm}
%\citestyle{dcu}
\usepackage[ruled,vlined]{algorithm2e}

\usepackage{natbib}
\usepackage{bibentry}
\usepackage[utf8]{inputenc}
\usepackage{bbm}
\bibliographystyle{agsm}
\citestyle{dcu}
\SetKwInput{KwInput}{Input}
\usepackage{xr}
\makeatletter
\newcommand*{\addFileDependency}[1]{
  \typeout{(#1)}
  \@addtofilelist{#1}
  \IfFileExists{#1}{}{\typeout{No file #1.}}
}
\makeatother

\newcommand*{\myexternaldocument}[1]{
    \externaldocument{#1}
    \addFileDependency{#1.tex}
    \addFileDependency{#1.aux}
}
\myexternaldocument{BDWNWADC3}


%Numbered environment
%\newcounter{example}[section]
%\newenvironment{example}[1][]{\refstepcounter{example}\par\medskip
%\noindent \textbf{Example~\theexample. #1} \rmfamily}{\medskip}


%Numbered environment defined with Newtheorem
\usepackage{amsmath}
\newtheorem{SampleEnv}{Sample Environment}[section]


%adjusting margins
%\usepackage{geometry}
%\geometry{paper=letterpaper, top=1.25 in, inner=1.25 in, outer=1.25 in, bottom=1 in, footskip=.5 in}

\textwidth 6.5 in \hoffset -.8in \textheight 9 in \voffset -.8in

\newtheorem{@assumption}{\sc Assumption}[section]
\newenvironment{assumption}{\begin{@assumption}\rm}{\end{@assumption}}

\newcommand{\calW}{{\mathcal W}}
\newcommand{\bbB}{{\Bbb B}}
\newcommand{\bbX}{{\Bbb X}}
\newcommand{\bbV}{{\Bbb V}}
\newcommand{\bbW}{{\Bbb W}}
\newcommand{\bbE}{{\Bbb E}}
\newcommand{\bbY}{{\Bbb Y}}
\newcommand{\bbZ}{{\Bbb Z}}
\newcommand{\bbR}{{\Bbb R}}
\newcommand{\bbN}{{\Bbb N}}
\newcommand{\bbC}{{\Bbb C}}
\newcommand{\bbL}{{\Bbb L}}
\newcommand{\bbP}{{\Bbb P}}
\newcommand{\bbQ}{{\Bbb Q}}
\newcommand{\bbF}{{\Bbb F}}
\newcommand{\bfd}{{\bf d}}
\newcommand{\bfT}{{\bf T}}
\newcommand{\bfA}{{\bf A}}
\newcommand{\bfB}{{\bf B}}
\newcommand{\bfX}{{\bf X}}
\newcommand{\bfU}{{\bf U}}
\newcommand{\bfY}{{\bf Y}}
\newcommand{\bfy}{{\bf y}}
\newcommand{\bfV}{{\bf V}}
\newcommand{\bfW}{{\bf W}}
\newcommand{\bfZ}{{\bf Z}}
\newcommand{\bfR}{{\bf R}}
\newcommand{\bfN}{{\bf N}}
\newcommand{\bfC}{{\bf C}}
\newcommand{\bfgamma}{{\boldsymbol \gamma}}
\newcommand{\rE}{\textrm{E}}
\newcommand{\E}{\mathbb{E}}
\newcommand {\cov}{\qopname\relax n{\textrm{Cov}}}
\newcommand {\var}{\qopname\relax n{\textrm{Var}}}
\newcommand {\corr}{\qopname\relax n{\textrm{Corr}}}
\newcommand {\Le}{\mathbb{L}}
\newcommand {\B}{{\mathcal{B}}}
\newcommand {\ind} {\overset{d}{=}}
\newcommand {\inD}{\overset{\mathcal{D}}{\rightarrow}}
\newcommand {\tod} {\overset{d}{\rightarrow}}
\newcommand {\topp} {\overset{p}{\rightarrow}}
\newcommand {\xn}{\{ X_n \}}
\newcommand {\fdd}{\overset{fdd}{\longrightarrow}}
\newcommand {\indf}[1]{{1}_{\{#1 \}}}
\newcommand {\argmin}[1]{\underset{#1}{\rm{argmin}}}
\newcommand {\argmax}[1]{\underset{#1}{\rm{argmax}}}
\newcommand {\argsup}[1]{\underset{#1}{\rm{argsup}}}
\newcommand {\wpsi}{{\psi}}
\newcommand {\D}[1]{\dot{#1}}
\newcommand {\revd}{\widetilde{d}}
\newcommand {\revb}{\widetilde{b}}
\newcommand {\revx}{\widetilde{X}}
\newcommand {\fbar}{\overline{F}}
\newcommand {\xbar}{\overline{X}}
\newcommand {\ybar}{\overline{Y}}
\newcommand {\normal}{{\mathcal N}}
\newcommand {\khat}{\what{k}}
\newcommand {\that}{\what{\tau}}
\newcommand {\xhat}{\what{X}}
\newcommand {\yhat}{\what{Y}}
\newcommand {\Rhat}{\what{R}}
\newcommand {\oml}{\omega_l}
\newcommand {\lwhat}{\what{d}_{lw}}
\newcommand {\gphhat}{\what{d}_{gph}}
\newcommand {\cusum}{\mbox{CUSUM}}
\newcommand {\adjcusum}{\mbox{adjCUSUM}}
\newcommand {\hnull}{\mbox{H}_0}
\newcommand {\halt}{\mbox{H}_1}
\newcommand {\rhohat}{ \widehat{\rho}}
\newcommand {\mhat}[1]{{#1}_t^{(m)}}
\newcommand {\pc}[1]{ #1^{\scriptsize pc}  }
\newcommand{\abs}[1]{\left\vert#1\right\vert}
\setcounter{algocf}{1}       

\newcommand{\pkg}[1]{{\normalfont\fontseries{b}\selectfont #1}}
\let\proglang=\textsf
\let\code=\texttt


\title
{Supplemental Technical Appendix to ``Block Dense Weighted Networks with Augmented Degree Correction"} 
\author{Benjamin Leinwand  and Vladas Pipiras\\
University of North Carolina at Chapel Hill
}


\begin{document}
\maketitle
\appendix

This supplement discusses several topics related to the main paper, and can be treated as an extended Appendix. Appendix \ref{s:nmicim1} proposes an alternative estimation approach for Normal LSM and LSM models using Nonnegative Matrix Factorization. Appendix \ref{sa:LSM-proof} contains two concentration inequalities for LSM networks using this alternative estimation approach, followed by auxiliary lemmas. Appendix \ref{sa:repeats} discusses estimation while accounting for certain diagonal 0's for within community subnetworks, along with how to adapt the estimation procedures when there are repeated values. Appendix \ref{mod-problems} explains how modularity may fail to characterize the kinds of communities discussed in the main paper. Appendix \ref{s:greedy-alg} introduces a greedy, agglomerative approach to community detection by trying to maximize the measure $L$.  Appendix \ref{s:spatial} considers community detection as a spatial clustering problem, and Appendix \ref{s:spectral-alg} includes an algorithm built upon spectral clustering to try to maximize $L$. Appendix \ref{s:rotated} shows that the real parts of the first several eigenvectors of a normalized version of the network also appear to capture community information. Appendix \ref{more-sims} displays additional simulated networks and results. Finally, Appendix \ref{s:Hext} extends both $H$-functions and the kinds of ``errors" which may be observed in LSMs or NSMs. 

\section{Estimation for Normal LSM and LSM}\label{s:nmicim1}
\subsection{Estimation for Normal LSM}\label{s:nmicim}
Recall from Section \ref{s:cosmo} that a Normal LSM has the form 
\begin{align}
f(W_{uv}) & = \gamma + \alpha \Phi_1^{-1}(\Psi_u) + \beta \Phi_1^{-1}(\Psi_v) + \sigma \epsilon_{uv} \nonumber\\
& =  \gamma + \alpha Z_u + \beta Z_v + \sigma \epsilon_{uv}.\nonumber
\end{align}

It is for now assumed that no transformation is needed, that is, $f(W_{uv}) = W_{uv}$. We also assume for simplicity that $\gamma =0$. Then, the model (\ref{e:NLSM2}) can be expressed after exponentiation as
\begin{equation}\label{e:exp_edge}
e^{W_{uv}} = e^{\alpha Z_u}e^{\beta Z_v}e^{\sigma \epsilon_{uv}} = e^{\alpha Z_u}e^{\beta Z_v}e^{\frac{\sigma^2} {2}}\left(\frac{e^{\sigma \epsilon_{uv}}}{\E e^{\sigma \epsilon_{uv}}}\right),
\end{equation}
where the last term within the parentheses has an expected value of 1. The structure of (\ref{e:exp_edge}) enables the use of rank-one Nonnegative Matrix Factorization (NMF) to estimate the parameters of interest. NMF approximates the matrix represented by $e^{W_{uv}}$ as the decomposition $ab^\prime$ where $a = (a_u)$ and $b = (b_v)$ are column vectors with positive entries. Taking the log of this approximation yields the approximate identity
\begin{equation}\label{e:NMF}
\log(a_u) + \log(b_v) \approx \alpha Z_u +\beta Z_v+ \frac{\sigma^2}{2},
\end{equation}
which suggests the following estimators for the model parameters of interest:
\begin{itemize}
  \item[] $\widehat{\alpha}$ = SD$(\log(a))$,  \, $\widehat{\beta}$ = SD$(\log(b))$,

    \item[] $\widehat{Z}_u =  \frac{\log(a_u) - \overline{\log(a)}}{\widehat{\alpha}}$, \,  $\widehat{\Psi}_u = \Phi_1(\widehat{Z}_u)$, \,   $\widehat{Z}_v =  \frac{\log(b_v) - \overline{\log(b)}}{\widehat{\beta}}$, \, $\widehat{\Psi}_v = \Phi_1(\widehat{Z}_v)$,
\end{itemize} 
where SD stands for the standard deviation, and $\overline{\log(x)}$ indicates the sample mean of $\log(x)$. In light of (\ref{e:NLSM2}), we also set $\widehat{\sigma}$ = SD$(W_{uv} - \widehat{\alpha}  \widehat{Z}_u- \widehat{\beta} \widehat{Z}_v)$.  An adjustment for subnetworks where all nodes are in the same community is given in Appendix \ref{sa:repeats} below.
%A concentration inequality for a Normal LSM bounding the difference between the best possible estimates of the network's SC to the true generating SC process without any ``error" included is given in Appendix \ref{sa:LSM-proof} of \cite{LeinwandPipiras}.

\subsection{Extension to LSM}\label{s:nmicim2}

The difference between (\ref{e:NLSM2}) and (\ref{e:LSM}), ignoring extra subscripts, is that rather than having standard normal random variables ascribed to each node, (\ref{e:LSM}) includes random variables $h_1(\Psi_u)$ and $h_2(\Psi_v)$ with possibly different forms, albeit with identical first 2 moments. The procedure described in Appendix \ref{s:nmicim} will still apply with one exception. The relation (\ref{e:NMF}) cannot directly estimate $\widehat{Z}_u$ or $\widehat{Z}_v$, but rather $\widehat{h_1(\Psi_u)}$ and $\widehat{h_2(\Psi_v)}$. After getting these estimates, a distribution can be fit to the data points while assuming that $\Psi_u$ and $\Psi_v$ are truly distributed uniformly over the unit interval. The best fitting distribution can then be inverted to estimate $\widehat{\Psi}_u$ and $\widehat{\Psi}_v$. Finally, the parameter $\sigma$ is estimated to be the standard deviation of 
$W_{uv} - \widehat{\alpha}  \widehat{h_1(\Psi_u)}- \widehat{\beta} \widehat{h_2(\Psi_v)}$.

\section{Concentration inequality for Normal LSM}\label{sa:LSM-proof}
In Section \ref{s:intronmici}, we introduced methods for estimating sociability parameters for certain kinds of network generating models. It would be helpful to ensure these estimates are actually capturing the underlying ``error-free" network generating mechanisms, what we refer to as the network's SC. As there are several introduced models, the accuracy of the estimation results may depend on the specific functional forms of a given network. We present below a result on estimation accuracy of the SC for a Normal LSM network (\ref{e:NLSM2}). The estimation procedure is formulated through theoretical means; it remains to be seen how the procedure compares to the practical estimation approach taken in Sections \ref{s:intronmici} and \ref{s:comm-det}. The proof of the result adapts ideas from \cite{noroozi2019estimation}. 

In the theorems that follow, we assume there are true communities (with number of communities $K_*$) but the estimates $\hat{P}$ do not necessarily use these true communities, but rather chooses a partition that minimizes the resulting estimate’s Frobenius distance to the observed matrix with an added penalty term based on the number of nodes and estimated number of communities, $\hat{K}$. The penalty term is necessary in these proofs, because if no penalty were included, we could simply let each node be its own community (setting $\hat{K} =n$), then perfectly find a formula for each single edge weight that connects each pair of communities/nodes.  Having each node serve as a different community defeats the purpose of clustering nodes into communities. If we did so, we would have 2 free parameters to fit each edge. Though this procedure would make the Frobenius norm 0, this would surely be overfitting. 

 Including the penalty leads to a different clustering more in the spirit of the models introduced in the paper, which still yields a concentration inequality. We do not use the measure of clustering quality discussed in the paper as a means of estimating $\hat{P},$ instead comparing the minimizing estimate to the observed network. Establishing the theoretical consistency of the community detection approach, as well as the non-asymptotic upper bound on the proportion of the misclassified nodes are outside the scope of this paper, but are topics which merit further investigation.   

%Fortunately, the structure of a Normal LSM (\ref{e:NLSM}) network allows for a relatively direct estimation process, which in turn admits a high probability bound on the distance between the estimated network and the true edge weight generating process.  
 \subsection{Concentration for Normal LSM with known $\sigma_{max}^2$}
\begin{theorem}\label{T:conc}
Let  $A$ be a Normal LSM network such that every edge weight is generated as in (\ref{e:NLSM2}), and $P_*$ be the network such that each edge weight  has the same generating process as the corresponding edge weight in $A$, but with every $\sigma$ value set to 0. In particular, all $\{\epsilon_{uv}\}$ are i.i.d. standard normal random variables that are independent of all other model parameters, and all $\{\Psi\} \sim $ U$(0,1)$.  Also let $\hat{P}$ be the estimated network (of the form (\ref{e:NLSM2}) with $\sigma =0$)  induced by the clustering of nodes that minimizes 
\begin{equation} \label{e:frob-and-pen}
||A- \hat{P}||_F^2 + \tex{Pen}(n , \hat{K}),
\end{equation}
where $\hat{K}$ is the number of different communities in this ``best" clustering, and 
\begin{equation} \label{e:pen1}
Pen(n, K) = 6 \sigma_{max}^2 \left(C_1 nK + C_2 K^2 \log(n)\right)  + (6C_3 + \frac{2}{c}) \sigma_{max}^2(\log(n)+n\log(K)),
\end{equation}
where $n$ is the number of nodes in the network, $K$ is the number of communities in the clustering, $c$ is some constant such that $0<c<1$, $\sigma_{max}^2 < \infty$ is the largest variance parameter of any generating function for edges in network $A$, and $C_1$ and $C_2$ are as given in Lemma \ref{L:NLSM1}.
Then, for some constant $C_3$ and any $t >0$,
\begin{equation}\label{e:ineq1}
\mathbb{P}\left(||\hat{P} -P_*||_F^2 \le (1-c)^{-1}  Pen(n, K_*)+  \frac{C_3}{c}\sigma_{max}^2t\right) \ge 1-3e^{-t}.
\end{equation}
\end{theorem}

Dividing by the number of edges $n^2$, the distance from the estimate to the truth continues to shrink on a per edge basis. As an immediate consequence of the theorem, if $K\sqrt{log(n)}/n \rightarrow 0$, the estimates are consistent.   

\begin{proof}

It is first useful to define the $\hat{P}$ induced by a particular clustering. Assuming we have clustered all nodes of $A$ into $\hat{K}$ estimated communities, we denote the subnetwork including only edges connecting nodes in estimated community $\hat{k}$ to nodes in estimated community $\hat{l}$ as $A^{(\hat{k}, \hat{l})}$. % For the within community submatrix $A^{(\hat{k}, \hat{k})}$, estimate the parameters in (\ref{e:NLSM2}) as follows:
%$$\hat{\gamma}^{(\hat{k}, \hat{k})} = \frac{1}{n_{\hat{k}}(n_{\hat{k}}-1)}\sum\limits_{u,v \in \hat{k}, \, v \ne u}A_{uv}, \qquad \qquad \hat{\alpha}^{(\hat{k}, \hat{k})}= \frac{SD(A^{(\hat{k}, \hat{k})})}{\sqrt{2}}, \qquad \qquad \hat{Z}_u^{(\hat{k})} = \frac{\mathcal{A}^{(\hat{k})}_u}{2\hat{\alpha}^{(\hat{k}, \hat{k})}},$$
%where $$SD(A^{(\hat{k}, \hat{k})}) = \sqrt{\frac{\sum\limits_{u,v \in \hat{k}, \, v \ne u} \left(A_{uv} - \hat{\gamma}^{(\hat{k}, \hat{k})}\right)^2}{n_{\hat{k}}(n_{\hat{k}}-1)-1}}$$ and     $$\mathcal{A}^{(\hat{k})}_u = \frac{\Big(2\sum\limits_{u \ne v, \, v \in \hat{k}} (A_{uv} - \hat{\gamma}^{(\hat{k}, \hat{k})})\Big) - \Big(\frac{1}{n_{\hat{k}}-2} \sum\limits\limits_{u \ne v, \, v \in \hat{k}} \sum\limits_{q\in \hat{k}, \, q \ne u, v} (A_{vq} - \hat{\gamma}^{(\hat{k}, \hat{k})})\Big)}{n_{\hat{k}}-1}\,. $$
%None of these parameter estimates incorporate the certain 0's on the diagonal. $\hat{\gamma}^{(\hat{k}, \hat{k})}$ is the the average edge weight in the subnetwork, similar to $\gamma$ in (\ref{e:NLSM2}), and $\hat{\alpha}^{(\hat{k}, \hat{k})}$ represents the estimated scaling parameter of within community edge weights with respect to the estimated $\hat{Z}$ values (which are assumed to be distributed according to a standard normal). As this is a within community subnetwork, there is only a single scaling parameter, that is, $\hat{\alpha}^{(\hat{k}, \hat{k})} =\hat{\beta}^{(\hat{k}, \hat{k})}$, which also explains the \sqrt{2}$ divisor. Finally, the $\hat{Z}_u^{(\hat{k})}$ estimates can be calculated using (\ref{e:self-loop}), as if self loops were allowed, $\mathbb{E}(A_{uu}|Z) = 2\alpha Z_u$.  As the estimated $\hat{\alpha}$ parameter within communities is always non-negative, this serves as a constraint to ensure the $\hat{Z}_u$ values are identifiable. Were we to allow the within community $\hat{\alpha}$ to be possibly negative, $\hat{\alpha}$ and all $\hat{Z}_u$ values could be multiplied by $-1$ to achieve the same result. 
Letting $\mathbbm{1}_n$ be a length $n$ column vector of 1's, define $\Pi(A^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ as the projection of $A^{(\hat{k}, \hat{l})}$ of the form (\ref{e:NLSM2}) with $\sigma =0$ which minimizes the Frobenius norm to the observed subnetwork, which can be written as 
\begin{equation}\label{e:submat-projection}
\Pi(A^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})= \hat{\gamma}^{(\hat{k}, \hat{l})}\mathbbm{1}_{\hat{n}_{\hat{k}}}\mathbbm{1}_{\hat{n}_{\hat{l}}}' + \hat{\alpha}^{(\hat{k}, \hat{l})} \hat{Z}^{(\hat{k}, \hat{l})} \mathbbm{1}_{\hat{n}_{\hat{l}}}' + \hat{\beta}^{(\hat{k}, \hat{l})}\mathbbm{1}_{\hat{n}_{\hat{k}}}\hat{Z}^{(\hat{l}, \hat{k})}', 
\end{equation}
such that
\begin{equation}\label{e:submat-argmin}
\hat{\gamma}^{(\hat{k}, \hat{l})}, \hat{\alpha}^{(\hat{k}, \hat{l})}, \hat{\beta}^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{l}, \hat{k})}  = \argmin {\gamma, \alpha, \beta, Z_1, Z_2} ||A^{(\hat{k}, \hat{l})} - (\gamma\mathbbm{1}_{\hat{n}_{\hat{k}}}\mathbbm{1}_{\hat{n}_{\hat{l}}}' + \alpha Z_1 \mathbbm{1}_{\hat{n}_{\hat{l}}}' + \beta \mathbbm{1}_{\hat{n}_{\hat{k}}}Z_2')||_F. 
\end{equation}
For identifiability purposes, we also mandate the following constraints:
\begin{equation}
\sum \hat{Z}^{(\hat{k}, \hat{l})} = 0, \qquad \sum \hat{Z}^{(\hat{l}, \hat{k})}= 0,  \qquad SD(\hat{Z}_1) = 1,  \qquad SD(\hat{Z}_2) = 1, \qquad \hat{\alpha}^{(\hat{k}, \hat{l})} \ge 0,  \qquad \hat{\beta}^{(\hat{l}, \hat{k})} \ge 0.
\end{equation}

For this proof, we need not explicitly express every component quantity in (\ref{e:submat-argmin}), but note that $\Pi(A^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ has rank $\le 3$ by construction. In the between estimated community subnetworks where $\hat{k} \ne \hat{l}$, the estimates $\hat{P}^{(\hat{k}, \hat{l})} = \Pi(A^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$.

However, for within estimated community subnetworks where $\hat{k} = \hat{l}$, we must define $\Pi(A^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ slightly differently, so as to ignore any influence from certain 0's on the diagonal of $A^{(\hat{k}, \hat{k})}$. In that case, we keep the pertinent constraints while modifying (\ref{e:submat-projection}) and (\ref{e:submat-argmin}) as follows:
\begin{equation}
\Pi(A^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})= \hat{\gamma}^{(\hat{k}, \hat{k})}\mathbbm{1}_{\hat{n}_{\hat{k}}}\mathbbm{1}_{\hat{n}_{\hat{k}}}' + \hat{\alpha}^{(\hat{k}, \hat{k})} \hat{Z}^{(\hat{k}, \hat{k})} \mathbbm{1}_{\hat{n}_{\hat{k}}}' + \hat{\alpha}^{(\hat{k}, \hat{k})}\mathbbm{1}_{\hat{n}_{\hat{k}}}\hat{Z}^{(\hat{k}, \hat{k})}', 
\end{equation}
such that
\begin{equation}
\hat{\gamma}^{(\hat{k}, \hat{k})}, \hat{\alpha}^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k}, \hat{k})}  = \argmin {\gamma, \alpha, \beta, Z} ||\big(A^{(\hat{k}, \hat{k})} - (\gamma\mathbbm{1}_{\hat{n}_{\hat{k}}}\mathbbm{1}_{\hat{n}_{\hat{k}}}' + \alpha Z \mathbbm{1}_{\hat{n}_{\hat{k}}}' + \alpha \mathbbm{1}_{\hat{n}_{\hat{k}}}Z')\big)_{u>v}||_F. 
\end{equation}
This is the projection which minimizes the Frobenius norm to the original subnetwork only with respect to off-diagonal entries, i.e. where $u \ne v$. $\Pi(A^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ also has rank $\le 3$, but this projection does not necessarily have 0's on the diagonal. This gives rise to within community edge estimates, $\hat{P}^{(\hat{k}, \hat{k})}$, which is equal to $\Pi(A^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ in the off diagonal entries, but replaces the diagonal entries with 0's (which are correct by construction). We also represent this estimate as $\Pi_0(A^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$, where the 0 subscript indicates the diagonals are forced to 0. %Note that the estimates of $\hat{Z}$ do not incorporate any information from the 0's on the diagonal of $A$, so 
$\hat{P}^{(\hat{k}, \hat{k})}$ is not a projection of $A^{(\hat{k}, \hat{k})}$, but is closer to $A^{(\hat{k}, \hat{k})}$ in Frobenius norm than $\Pi(A^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$, which is itself the best projection of $A^{(\hat{k}, \hat{k})}$ with respect to the Frobenius norm for off diagonal entries.  

The projections defined on the observed network and the estimated communities differ from the projections we define on $P_*$ and the estimated communities. In defining $\Pi(P_*^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$, or $\Pi_0(P_*^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$, we treat the $\hat{Z}$ values as fixed, using the estimated values from $\Pi(A^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$, and only allowing $\hat{\gamma}^{(\hat{k}, \hat{l})}, \hat{\alpha}^{(\hat{k}, \hat{l})},$ and $\hat{\beta}^{(\hat{k}, \hat{l})}$ values to vary. 
%In general, let $\Pi(M^{(i, j)},  Z^{(i)}, Z^{(j)})$ represent the projection of subnetwork $M^{(i, j)}$ onto the (true or estimated) $Z$ values of the nodes in (true or estimated) communities $i$ and $j$. When $\hat{k} \ne \hat{l}$, we can still use estimates of $Z_u$'s for all nodes from the within community estimation step. This is then a more straightforward least squares problem, estimating $\hat{\gamma}^{(\hat{k}, \hat{l})}, \hat{\alpha}^{(\hat{k}, \hat{l})},$ and $\hat{\beta}^{(\hat{k}, \hat{l})}$ by regressing the observed submatrix values on  $\hat{Z}^{(\hat{k})}$ and $\hat{Z}^{(\hat{l})}$ values, resulting in $\Pi(M^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$. For intercommunity edges, we therefore have estimates
%$$\hat{P}^{(\hat{k}, \hat{l})}= \Pi(A^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}) = \hat{\gamma}^{(\hat{k}, \hat{l})}\mathbbm{1}_{\hat{n}_{\hat{k}}}\mathbbm{1}_{\hat{n}_{\hat{l}}}' + \hat{\alpha}^{(\hat{k},\hat{l})}\hat{Z}^{(\hat{k})}\mathbbm{1}_{\hat{n}_{\hat{l}}}' + \hat{\beta}^{(\hat{k}, \hat{l})}\mathbbm{1}_{\hat{n}_{\hat{k}}}\hat{Z}^{(\hat{l})}' $$
%where each $\hat{P}^{(\hat{k}, \hat{l})}$ also has rank $\le 3$. 

With these preliminaries in place, we follow the proof of Theorem 1 in \cite{noroozi2019estimation}. By assumption, 
$$||A- \hat{P}||_F^2 + \tex{Pen}(n , \hat{K}) \le ||A- P_*||_F^2 + \tex{Pen}(n , K_*).$$
%If we replace the diagonal entries of both $\hat{P}$ and $P_*$ by the estimate of each self loop value in $P_*$ using the estimated communities similar to in (\ref{e:self-loop}), $$\hat{P}_{uu}= P_{*uu} = \frac{\left(2\sum\limits_{\substack{u \ne v \\ v \in \hat{k}}} (P_{*uv} - \hat{\gamma}^{(\hat{k}, \hat{k})})\right) - \left(\frac{1}{n_{\hat{k}}-2} \sum\limits_{\substack{u \ne v \\ v \in \hat{k}}} \sum\limits_{\substack{q \ne u, q \ne v \\ q\in \hat{k}}} (P_{*vq} - \hat{\gamma}^{(\hat{k}, \hat{k})})\right)}{n_{\hat{k}}-1}\,, $$ WE DO THIS SO THAT I CAN JUST USE THE PROJECTION BETWEEN (A.10) and (A.11)then, the assumed inequality above still holds since we replace the diagonal 0's in each matrix with the same values.  Due to this equivalence, we move forward with these altered matrices, but do not change the notation. 
Letting 
\begin{equation}\label{Xi-def}
\Xi = A - P_*, 
\end{equation}
writing $Tr(M)$ for the trace of matrix $M$, and assuming the network has been rearranged into blocks by estimated communities, adding and subtracting $P_*$ within $||A-\hat{P}||_F^2$ on the left-hand side gives
\begin{equation}\label{e:inequality-penalty}
||\hat{P} -P_*||_F^2 \le 2Tr(\Xi'(\hat{P} -P_*))  + \tex{Pen}(n , K_*) - \tex{Pen}(n , \hat{K}).
\end{equation}
Noting $2 Tr(\Xi'(\hat{P} -P_*)) = 2 \sum\limits_{\hat{k}, \hat{l}=1}^{\hat{K}}Tr\left(\Xi^{(\hat{k}, \hat{l})}'(\hat{P}^{(\hat{k}, \hat{l})} -P_*^{(\hat{k}, \hat{l})})\right)$, adding and subtracting $\Pi_0(P_*^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ to within estimated community subnetworks and $\Pi(P_*^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ to between estimated community subnetworks yields for the trace term in (\ref{e:inequality-penalty}):
%$$2 \sum\limits_{\hat{k}, \hat{l}=1}^{\hat{K}}Tr\left(\Xi^{(\hat{k}, \hat{l})}'(\hat{P}^{(\hat{k}, \hat{l})} - \Pi(P_*^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}) + \Pi(P_*^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}) - P_*^{(\hat{k}, \hat{l})})\right) $$
\begin{align}\label{broken-trace}
&2 \sum\limits_{\hat{k}\ne \hat{l}}Tr\left(\Xi^{(\hat{k}, \hat{l})}' \Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})\right) + 2 \sum\limits_{\hat{k}=1}^{\hat{K}}Tr\left(\Xi^{(\hat{k}, \hat{k})}' \Pi_0(\Xi^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})\right) & \nonumber \\
&+2 \sum\limits_{\hat{k} \ne \hat{l}}Tr\left(\Xi^{(\hat{k}, \hat{l})}' (\Pi(P_*^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}) - P_*^{(\hat{k}, \hat{l})})\right) +2 \sum\limits_{\hat{k}=1}^{\hat{K}}Tr\left(\Xi^{(\hat{k}, \hat{k})}' (\Pi_0(P_*^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})}) - P_*^{(\hat{k}, \hat{k})})\right). &
\end{align}
We handle the four sums in (\ref{broken-trace}) separately. Note that in the between community subnetworks, $Tr\left(\Xi^{(\hat{k}, \hat{l})}'\Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})\right) = ||\Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}))||_F^2$ since $\Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ is a projection. $\Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ has rank $\le 3$ as noted above, so  
\begin{equation}\label{Theorem-6Lemma1}
2||\Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})||_F^2 \le 6||\Pi(\Xi^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})||_{op}^2 \le 6||\Xi^{(\hat{k}, \hat{l})}||_{op}^2,
\end{equation}
where $||M||_{op}$ denotes the usual spectral norm of matrix $M$. The second inequality holds because the spectral norm is sub-multiplicative, and for a projection matrix $P$, $||P||_{op} \le 1$. 

When $\hat{k}=  \hat{l}$, $$Tr\left(\Xi^{(\hat{k}, \hat{k})}' \Pi_0(\Xi^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})\right) = Tr\left(\Xi^{(\hat{k}, \hat{k})}' \Pi(\Xi^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})\right),$$
as $Tr(M'N) = \sum\limits_{u,v} M_{uv}N_{uv}$ and the diagonal of $\Xi^{(\hat{k}, \hat{k})}$ is all 0's because there are no self loops in either $A$ or $P_*$.  With this equivalence, the same argument culminating in (\ref{Theorem-6Lemma1}) also holds for within community subnetworks.

To derive a bound on the first two sums of (\ref{broken-trace}), from Lemma \ref{L:NLSM2}, we have 
\begin{equation}\label{e:Lemma1-result}
\mathbb{P}\left(\sum\limits_{\hat{k},\hat{l}=1}^{\hat{K}} ||\Xi^{(\hat{k},\hat{l})}||_{op}^2 \le \sigma_{max}^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) + C_3 (t+ \log(n) + n \log(\hat{K}))) \right) \ge 1 - e^{-t}.
\end{equation}

Moving to the last two sums in (\ref{broken-trace}), these can be rewritten in terms of the whole network instead of estimated subnetworks. Slightly abusing notation, let $\Pi(P_*, \{\hat{k}\})$ represent the full network matrix where all entries are given by the value dictated by $\Pi(P_*^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ for between estimated community edges, and by $\Pi_0(P_*^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ for within estimated community edges. The last two sums in (\ref{broken-trace}) can be represented as  
\begin{equation}\label{intro-Xi-H}
2Tr\left(\Xi'(\Pi(P_*, \{\hat{k}\}) - P_*)\right) = 2||\Pi(P_*,  \{\hat{k}\}) - P_*||_F \, |\langle\Xi, H(\{\hat{k}\})\rangle|,
\end{equation}
where%, dropping superscripts and abbreviating each $\Pi(P_*^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ as $\Pi(P_*)$, 
$$H(\{\hat{k}\}) = \frac{\Pi(P_*, \{\hat{k}\}) - P_*}{||\Pi(P_*, \{\hat{k}\}) - P_*||_F}.$$
Since for any $a, b$ and for $c >0$, $2ab \le c a^2 +\frac{b^2}{c}$,
\begin{equation}\label{second-trace-breakdown}
2Tr\left(\Xi'(\Pi(P_*, \{\hat{k}\}) - P_*)\right) \le c  ||\Pi(P_*, \{\hat{k}\}) - P_*||_F^2 + \frac{ |\langle\Xi, H(\{\hat{k}\})\rangle|^2}{c}.
\end{equation}
Denoting the set of partitions of the nodes into exactly $K$ communities as $\mathcal{G}_K$, for any fixed partition $G \in \mathcal{G}_K$, $\sum_{u, v} (H(G)_{uv})^2 = 1$, and the matrix $\Xi$ consists of independent normally distributed errors with finite variances. Representing the true communities of nodes $u$ and $v$ as $i$ and $j$, respectively, since $|\langle \Xi, H(G)\rangle| = vec(\Xi)'vec(H(G))$, observe that, 
$$ \mathbb{P}\left(|\langle\Xi, H(G)\rangle|^2 >t\right) = \mathbb{P}\left((\sum\limits_{u, v=1}^{n} \sigma_{i j} \epsilon_{uv} H(G)_{uv})^2 >t\right)$$ 
\begin{equation}
= 2\mathbb{P}\left(\sum\limits_{u, v=1}^n \sigma_{i j} \epsilon_{uv} H(G)_{uv} >\sqrt{t}\right) \le 2\mathbb{P}\left(\mathcal{N}(0, \sigma_{max}^2)>\sqrt{t}\right) \le 2e^{\left-t/2\sigma_{max}^2\right}.
\end{equation}
%The arguments from (\ref{intro-Xi-H}) to this point also hold in those subnetworks where $k=l$ when replacing $\Pi$ with $\Pi_0$. This gives 
%\begin{equation}\label{Xi-H-inequality}
 %\mathbb{P}\left(\sum\limits_{k,l=1}^{K}|\langle \Xi^{(k, l)}, H^{(k, l)}\rangle |^2 >2\sigma_{max}^2(t+ 2\log(K))\right) \le 2K^2 e^{\left(-2\log( K) - t\right)} \le 2 e^{-t}.
%\end{equation}
Applying the union bound, 
$$\mathbb{P}\left(|\langle \Xi, H(\{\hat{k}\})\rangle |^2 - 2\sigma_{max}^2(\log(n) + n\log(\hat{K})) >2\sigma_{max}^2t\right) $$
$$\le \mathbb{P}\left(\max\limits_{1 \le K \le N} \max\limits_{G \in \mathcal{G}_K} \left(|\langle \Xi, H(G)\rangle |^2 - 2\sigma_{max}^2(\log(n) + n\log(K))\right) >2\sigma_{max}^2t \right) $$
\begin{equation}\label{Xi-H-union}
}\le 2nK^n e^{-2\sigma_{max}^2(\log(n) + n\log(K)+t)/ 2\sigma_{max}^2} = 2e^{-t}.
\end{equation}
The squared Frobenius norm matrix $||\Pi(P_*, \{\hat{k}\}) - P_*||_F^2$ in (\ref{second-trace-breakdown}) can be written as 
\begin{equation}
\sum \limits_{\hat{k}\ne \hat{l}} ||\Pi(P_*^{(\hat{k}, \hat{l})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}) - P_*^{(\hat{k}, \hat{l})}||_F^2 + \sum \limits_{\hat{k}=1}^{\hat{K}} ||\Pi_0(P_*^{(\hat{k}, \hat{k})},  \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})}) - P_*^{(\hat{k}, \hat{k})}||_F^2.
\end{equation}
When $\hat{k} \ne \hat{l}$, $||\Pi(P_*^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})}) - P_*^{(\hat{k}, \hat{l})}||_F \le ||\hat{P}^{(\hat{k}, \hat{l})} - P_*^{(\hat{k}, \hat{l})}||_F$ because $\Pi(P_*^{(\hat{k}, \hat{l})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{l})})$ is the \textit{best} projection of $P_*^{(\hat{k}, \hat{l})}$ onto the estimated $\hat{Z}^{(\hat{k}, \hat{l})}$ and $\hat{Z}^{(\hat{l}, \hat{k})}$ values. When $\hat{k} = \hat{l}$, \Pi(P_*^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ is also the best projection of $P_*^{(\hat{k}, \hat{k})}$ onto the $\hat{Z}^{(\hat{k}, \hat{k})}$ values with respect to only the off diagonal entries, and the diagonal entries of both $\Pi_0(P_*^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})})$ and $P_*^{(\hat{k}, \hat{k})}$ are all 0. Therefore,  $||\Pi_0(P_*^{(\hat{k}, \hat{k})}, \hat{Z}^{(\hat{k})}, \hat{Z}^{(\hat{k})}) - P_*^{(\hat{k}, \hat{k})}||_F \le ||\hat{P}^{(\hat{k}, \hat{k})} - P_*^{(\hat{k}, \hat{k})}||_F$. Combining this point with those given in (\ref{intro-Xi-H})--(\ref{Xi-H-union}), 
\begin{equation}\label{second-trace-ineq}
\mathbb{P}\left(2Tr\left(\Xi'(\Pi(P_*, \{\hat{k}\}) - P_*)\right) \le c ||\hat{P} - P_*||_F^2  + \frac{2\sigma_{max}^2(\log(n) +n\log(\hat{K})+ t)}{c}\right)  \ge 1-2e^{-t}.
\end{equation}

Let $\Omega$ denote the set on which, for a particular $t$, both events in (\ref{e:Lemma1-result}) and (\ref{second-trace-ineq}) occur. Then $P(\Omega) \ge 1-3e^{-t}$. Therefore on $\Omega$, by using (\ref{e:inequality-penalty}),
\begin{align}\label{on-omega-result}
||\hat{P} -P_*||_F^2 \le & 6 \sigma_{max}^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) + C_3(\log(n) + n\log(\hat{K}) + t)) \nonumber \\
&+ c ||\hat{P} - P_*||_F^2  + \frac{2\sigma_{max}^2(\log(n) + n\log(\hat{K}) +t)}{c} + \tex{Pen}(n , K_*) - \tex{Pen}(n , \hat{K}) .
\end{align}
Letting $0<c < 1$ and $Pen(n, K)$ as in (\ref{e:pen1}), we derive (\ref{e:ineq1}), that is $$\mathbb{P}\big(||\hat{P} -P_*||_F^2 \le (1-c)^{-1}  Pen(n, K_*)+  \frac{C_3}{c}\sigma_{max}^2t\big) \ge 1-3e^{-t}. \qquad \qquad \qquad \QED$$ 
\end{proof}

This theorem indicates that in the Normal LSM setting, by choosing an appropriate penalty, one can choose a clustering and the number of clusters such  that, with high probability, $||\hat{P}-P||_F^2$, the squared Frobenius distance between the true generating process and the estimate thereof induced by this ``best" clustering, is bounded above by a multiple of the penalty using the true number of communities, accounting for the scale of the $\sigma$ values in the network. As might be expected, in a network with many nodes, many communities, and large $\sigma$ values, the total distance between the estimate and generating process for the network enforced by this bound can become large. One drawback is this result relies on employing a penalty which includes $\sigma_{max}^2$, a value which cannot in all cases be assumed to be known. 
    
\subsection{Case with unknown $\sigma^2$}
In practice, we do not know the value of $\sigma_{max}^2$. If instead we assume that throughout the whole network, there is a single value of $\sigma^2$, we can also estimate that value using the data.  Furthermore, letting $\hat{K} =n$ and putting each node into its own community, the estimate $\hat{P}$ will appear to be error free. For any network of interest, this is overfitting, so it is reasonable to set some restriction on the number of communities relative to the number of nodes. The following theorem extends Theorem \ref{T:conc} under these new assumptions, where we do not require prior knowledge of $\sigma_{max}^2$.
\begin{theorem} In the same setting as Theorem \ref{T:conc}, assume the choice of $\hat{K}$ is restricted such that
\begin{equation}\label{community-limit}
C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) + C_3(\log(n)+n\log(\hat{K})) \le \frac{n(n-1)}{4},
%\frac{\hat{K}^2\log n}{n^2} \rightarrow 0.
\end{equation}
where the constants $C_1, C_2, C_3$ are derived in the same manner as in Theorem \ref{T:conc}, but take on slightly different values here. Additionally, assume network $A$ has a constant variance parameter $\sigma^2$. Let $\hat{P}$ be the estimated network induced by the clustering of nodes that minimizes
\begin{equation}\label{e:frob-and-pen2}
||A- \hat{P}||_F^2 + \hat{\sigma}^2 \tex{Pen}(n , \hat{K}),
\end{equation}
where 
\begin{equation}\label{proof2-sigma-hat}
\hat{\sigma}^2 = \frac{||A-\hat{P}||_F^2}{n(n-1)},
\end{equation}
and the penalty is given by 
\begin{equation} \label{e:pen2}
Pen(n, K) = 4(C_1 nK + C_2 K^2 \log(n) + C_3(\log(n)+n\log(K))).
\end{equation}
Then, for $t >0, \epsilon \in [0, 1/2)$, 
\begin{equation}\label{proof2-result}
\mathbb{P}\left(||P_*- \hat{P}||_F^2 \le \frac{\sigma^2}{1-c}\left(C t + (1+\epsilon) \tex{Pen}(n , K_*)\right) \right) \ge 1 - 3e^{-t} -  e^{-\frac{3}{32}\epsilon^2 n(n-1)}. % e^{-\tilde{c}n^2}.$$
\end{equation}
\end{theorem}
\begin{proof}    
By assumption,  
$$||A- \hat{P}||_F^2 + \hat{\sigma}^2 \tex{Pen}(n , \hat{K}) \le ||A- P_*||_F^2 + \sigma_*^2 \tex{Pen}(n , K_*),$$
where $\hat{\sigma}^2 = \frac{||A-\hat{P}||_F^2}{n(n-1)}$ as in (\ref{proof2-sigma-hat}), and $\sigma_*^2 = \frac{||A-P_*||_F^2}{n(n-1)}$. The above is equivalent to
$$||A- \hat{P}||_F^2 \left(1+\frac{\tex{Pen}(n , \hat{K})}{n(n-1)}\right) \le ||A- P_*||_F^2 \left( 1+ \frac{\tex{Pen}(n , K_*)}{n(n-1)}\right).$ Then, following the proof of Theorem \ref{T:conc} leading to (\ref{e:inequality-penalty}), we have\\
\begin{align}
&||P_*- \hat{P}||_F^2 \left(1+\frac{\tex{Pen}(n , \hat{K})}{n(n-1)}\right)& \nonumber \\
&\le  2Tr(\Xi'(\hat{P} -P_*)) \left(1+\frac{\tex{Pen}(n , \hat{K})}{n(n-1)}\right) + \frac{||A- P_*||_F^2}{n(n-1)}\left(\tex{Pen}(n , K_*)- \tex{Pen}(n , \hat{K})\right). \nonumber
\end{align}
Dividing both sides by $1+\frac{\tex{Pen}(n , \hat{K})}{n(n-1)}$ and following the arguments culminating in (\ref{on-omega-result}) but adjusting constants as necessary, with probability $\ge 1-3e^{-t}$:
\begin{align}
&(1-c)||P_*- \hat{P}||_F^2 & \nonumber \\
&\le \sigma^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) + C_3 (\log(n)+n\log(\hat{K}))+ Ct) + \frac{\sigma_*^2 \tex{Pen}(n , K_*)}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}} - \frac{\sigma_*^2 \tex{Pen}(n , \hat{K})}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}}& \nonumber
\end{align}
\begin{equation}\label{drop-denom}
\le \sigma^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) +  C_3 (\log(n)+n\log(\hat{K})) + Ct) + \sigma_*^2 \tex{Pen}(n , K_*)- \frac{\sigma_*^2 \tex{Pen}(n , \hat{K})}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}},
\end{equation}
where the denominator from the second to last term was dropped in the last inequality. 

Note that $\sigma_*^2$ is the average of $\frac{n(n-1)}{2}$ squared $\mathcal{N}(0, \sigma^2)$ random variables, so its distribution is  $\frac{2\sigma^2}{n(n-1)} \chi_{\frac{n(n-1)}{2}}^2$. Using a result from \cite{johnstone2001chi}, for $\epsilon \in [0, 1/2)$,  
\begin{equation}\label{sigma-vs-sigma}
\mathbb{P}\big((1-\epsilon)\sigma^2 \le \sigma_*^2 \le (1+\epsilon)\sigma^2\big)  \ge 1- e^{-\frac{3}{32}n(n-1)\epsilon^2}.
\end{equation}
%$$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!######!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!$$
%Using another large deviation result (WE NEED TO JUSTIFY THIS RESULT), for $\epsilon>0$ ,with probability $\ge 1 - e^{-\tilde{c}n^2}$,
%$$(1-\epsilon)\sigma^2 \le \sigma_*^2 \le (1+\epsilon)\sigma^2.$$ 
Putting together (\ref{drop-denom}) and (\ref{sigma-vs-sigma}), for $\epsilon \in [0, .5)$,  with probability $\ge 1- 3e^{-t} - e^{-\frac{3}{32}\epsilon^2 n(n-1)}$, 
\begin{align}\label{pre-improper}
(1-c)||P_*- \hat{P}||_F^2 &\le C\sigma^2 t + (1+\epsilon)\sigma^2 \tex{Pen}(n , K_*) \nonumber \\ 
&+\sigma_*^2\left(\frac{C_1}{1-\epsilon} n\hat{K} + \frac{C_2}{1-\epsilon} \hat{K}^2 \log(n) +\frac{C_3}{1-\epsilon}(\log(n)+n\log(\hat{K})) - \frac{\tex{Pen}(n , \hat{K})}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}}\right).&
\end{align}
Using (\ref{community-limit}), (\ref{e:pen2}), and $\epsilon < 1/2$, the term in parentheses in (\ref{pre-improper}) is bounded by
$$2(C_1n\hat{K} + C_2\hat{K}^2 \log(n) +C_3(\log(n)+n\log(\hat{K}))) - \frac{\tex{Pen}(n , \hat{K})}{2}\le 0.$$
%\begin{align}\label{improper-frac}
%&\frac{\sigma_*^2}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}}\left( \frac{C_1}{1-\epsilon} n\hat{K} + \frac{C_2}{1-\epsilon} \hat{K}^2 \log(n) +\frac{C_3}{1-\epsilon}(\log(n)+n\log(\hat{K}))\right)& \nonumber \\
%&-\frac{\sigma_*^2}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}}\left(\tex{Pen}(n , \hat{K}) \left(1 - \left(\frac{C_1\hat{K}}{(1-\epsilon)(n-1)}  + \frac{C_2 \hat{K}^2 \log(n) + C_3 (\log(n)+ n\log(\hat{K}))}{(1-\epsilon)n(n-1)} \right)\right)\right).&
%\end{align}
%By the assumption (\ref{community-limit}) on the choice of $\hat{K}$, the second line of (\ref{improper-frac}) behaves asymptotically as $-\frac{\sigma_*^2 \tex{Pen}(n , \hat{K})}{1 + \frac{\tex{Pen}(n , \hat{K})}{n(n-1)}}$. 
%Using the penalty given by (\ref{e:pen2}), the second line of   (\ref{improper-frac}) becomes larger than the first line, so the whole expression becomes negative. 
Hence, with probability $\ge 1- 3e^{-t} - e^{-\frac{3}{32}n(n-1)\epsilon^2}$, $$||P_*- \hat{P}||_F^2 \le \frac{C\sigma^2 t + (1+\epsilon)\sigma^2 \tex{Pen}(n , K_*)}{(1-c)}.\qquad \qquad \QED$$
\end{proof}

\subsection{Auxiliary results}%Lemma to bound  $\sum\limits_{k,l}^K ||\Xi^{(k,l)}||_{op}^2$
The following results were used in the proof of Theorem \ref{T:conc}. 

\begin{lemma} \label{L:NLSM1} 
Let $\Xi$ be a symmetric $n\times n$ matrix with 0's on the diagonal and independent $\mathcal{N}(0, \sigma_{uv}^2)$ entries above the diagonal, where all $\sigma_{uv}^2 \le \sigma_{max}^2$. Let $\Xi$ be partitioned into $K^2$ submatrices $\Xi^{(k, l)}, k, l = 1, \ldots, K$ 
%Let $\Xi = A - P_*$, and denote the subnetwork including only edges connecting nodes in estimated community $\hat{k}$ to nodes in estimated community $\hat{l}$ as $A^{(\hat{k}, \hat{l})}$. Define $n$, $\hat{K}$ and $\sigma_{max}^2$  as in Theorem \ref{T:conc}. 
For some constants $C_1, C_2$, and $t>0$, 
\begin{equation}\label{Lemma1}
\mathbb{P}\left(\sum\limits_{k, l =1}^{K} ||\Xi^{(k, l)}||_{op}^2 \le \sigma_{max}^2 (C_1 nK + C_2 K^2 \log(n) + C_3 t) \right) \ge 1 - e^{-t}.
\end{equation}
\end{lemma}
\begin{proof}For a fixed partition, let $\xi$ and $\mu$ be vectors with entries $\xi_{k, l} = ||\Xi^{(k,l)}||_{op} $ and   \mu_{k,l} = \mathbb{E}||\Xi^{k,l}||_{op}$ and $\eta = \xi - \mu$. Then,
$$ \sum\limits_{k,l=1}^{K} ||\Xi^{(k, l)}||_{op}^2 = ||\xi||^2 \le 2||\eta||^2 + 2||\mu||^2.$$
We start by bounding $||\mu||^2$. Using Theorems 1.1 and 3.1 from \cite{Bandeira_2016}, letting $n_{k}$ and $n_{l}$ denote the number of rows and columns respectively in $\Xi^{(k, l)}$,
 $$\mu_{k, l} = \mathbb{E}||\Xi^{(k, l)}||_{op} \le (1+\epsilon)\left(\sqrt{n_{k}}\sigma_{max} + \sqrt{n_{l}}\sigma_{max} + \frac{6}{\sqrt{\log(1+\epsilon)}}\sigma_{max}\sqrt{\log(\min(n_{k}, n_{l}))}\right)$$ for any $0< \epsilon \le  1/2$. In other words, 
$$\mu_{k,l} \le C_0 \sigma_{max}(\sqrt{n_{k}} + \sqrt{n_{l}} +\sqrt{\log(\min(n_{k}, n_{l}))}\, ),$$
$$\mu_{k,l}^2 \le 3C_0^2 \sigma^2_{max}(n_{k} + n_{l} +\log(\min(n_{k}, n_{l}))),$$
\begin{equation}\label{L:mu}
||\mu||^2 \le 3 C_0^2 \sigma^2_{max}\sum\limits_{k,l=1}^{K} (n_{k} + n_{l} +\log(\min(n_k}, n_{l}))) \le 6C_0^2 \sigma^2_{max}nK + 3C_0^2 \sigma^2_{max}K^2 \log(n).
\end{equation}
Next, for $1\le k \le l \le K$, $\eta_{k,l} = \xi_{k,l} - \mu_{k, l}$ are all independent random variables since all errors are assumed to be independent. By Theorem 5.8 of \cite{boucheron2013concentration}, $$\mathbb{P}(|\eta_{k, l}| \ge t ) = \mathbb{P}(|\xi_{k,l} - \mu_{k,l}| \ge t ) \le  2e^{{\frac{-t^2}{4\sigma_{max}^2}}},$$ so $\eta_{k,l}$ is sub-gaussian. Since $\mathbb{E}(\eta_{k,l})=0$, using sub-gaussianity, from Proposition 2.5.2 of \cite{vershynin2018high}, there exists a constant $C \le 288e$ such that $$\mathbb{E}(e^{t\eta_{k,l}}) \le e^{\frac{C\sigma_{max}^2 t^2}{2}}.$$
Let $\tilde{\eta}$ be the sub-vector of $\eta$ which includes the $\eta_{k,l}$ values for $1\le k \le l \le K$. Then, Theorem 2.1 of \cite{Hsu_Kakade_Zhang_2012} ensures that, for any square matrix $M$, using the same constant $C$,

$$\mathbb{P}\left(||M\tilde{\eta}||^2 > C\sigma_{max}^2(Tr(M'M) + 2\sqrt{Tr((M'M)^2)t} + 2||M'M||_{op} t )\right) \le e^{-t}.$$
Letting $M=I_{K(K+1)/2}$, this becomes
$$\mathbb{P}\left(||\tilde{\eta}||^2 \ge C\sigma_{max}^2(K(K+1)/2 + \sqrt{2K(K+1)t} +2t\right) \le e^{-t},$$
and since $||\eta||^2 \le 2||\tilde{\eta}||^2$, 
\begin{equation}\label{L:eta}
\mathbb{P}\left(||\eta||^2 \le 2C\sigma_{max}^2K(K+1) + 6 C\sigma_{max}^2 t \right) \ge 1 - e^{-t}.
\end{equation}
Combining (\ref{L:mu}) and (\ref{L:eta}), 
$$\mathbb{P}\left( |\xi||^2 \le 12C_0^2 \sigma^2_{max}nK + 6C_0^2 \sigma^2_{max}K^2 \log(n) + 4C \sigma_{max}^2 K(K+1) + 12 C \sigma_{max}^2  t\right)\ge 1 - e^{-t}.$$
Collecting terms yields (\ref{Lemma1}). \hfill \QED
%\mathbb{P}\left(\sum\limits_{\hat{k},\hat{l}}^{\hat{K}} ||\Xi^{(\hat{k},\hat{l})}||_{op}^2 \le \sigma_{max}^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n)) + C_3\hat{K}(\hat{K}+1) + C_4 t \right) \ge 1 - e^{-t}.$$
\end{proof}\\


The following lemma is nearly identical to Lemma 6 in \cite{noroozi2019estimation}, and serves the same purpose, to translate the bound in Lemma \ref{L:NLSM1}, which is conditional on a particular partition, into an unconditional bound. 
\begin{lemma}\label{L:NLSM2}
For any t $> 0$, 
\begin{equation}
\mathbb{P}\left(\sum\limits_{\hat{k},\hat{l}=1}^{\hat{K}} ||\Xi^{(\hat{k},\hat{l})}||_{op}^2 \le \sigma_{max}^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) + C_3 (t+ \log(n) + n \log(\hat{K}))) \right) \ge 1 - e^{-t}.
\end{equation}
\end{lemma}
\begin{proof}
Denoting the set of partitions of the nodes into $K$ communities as $\mathcal{G}_K$, for any fixed partition $G \in \mathcal{G}_K$, from Lemma \ref{L:NLSM1},
$$\mathbb{P}\left(\sum\limits_{k, l =1}^{K} ||\Xi^{(k, l)}||_{op}^2 \ge \sigma_{max}^2 (C_1 nK + C_2 K^2 \log(n) + C_3 x) \right) \le  e^{-x}.$$
Taking a union bound over all possible partitions and setting $x = t+ \log(n) + n \log(K)$ 
$$\mathbb{P}\left(\sum\limits_{\hat{k},\hat{l}=1}^{\hat{K}} ||\Xi^{(\hat{k},\hat{l})}||_{op}^2 -  \sigma_{max}^2 (C_1 n\hat{K} + C_2 \hat{K}^2 \log(n) + C_3 (t + \log(n) + n \log(\hat{K}))) \ge 0 \right)$$
$$ \le \mathbb{P} \left( \max \limits_{1\le K \le n}\max \limits_{G \in\mathcal{G}_K } \sum\limits_{k, l =1}^{K} ||\Xi^{(k, l)}||_{op}^2 -  \sigma_{max}^2 (C_1 nK + C_2 K^2 \log(n) + C_3 (\log(n) + n \log(K))) \ge \sigma_{max}^2 C_3 t  \right)$$
$$\le \sum \limits_{K=1}^n \sum \limits_{G \in\mathcal{G}_K } \mathbb{P} \left(\sum\limits_{k, l =1}^{K} ||\Xi^{(k, l)}||_{op}^2 - \sigma_{max}^2 (C_1 nK + C_2 K^2 \log(n) + C_3 (\log n + n \log(K)))\ge\sigma_{max}^2 C_3 t \right) $$
$$\le nK^n e^{-t -\log(n) -n\log(K)} = e^{-t} \qquad \qquad \QED$$  





\section{Additional estimation details}\label{sa:repeats}
\subsection{Pre-estimating the diagonal for within community NMF}
When applying the methodology described in Appendix \ref{s:nmicim} to a set of edges within the same community, symmetry ensures that $a= b$, so there are no conflicts among estimators. There is, however, a need to resolve the fact that no self loops ($W_{uu} =0$) may impact the estimates of $a$ and $\alpha$ for the within community setting $i=j$. To avoid this impact, the values of the diagonal entries $W_{uu}$ should be imputed before computing the NMF. Were self loops allowed, the expected value a self loop would be $\mathbb{E}(W_{uu}|Z) = 2 \alpha Z_u$, so we can estimate this value while accounting for all diagonal zeros and estimates of other $Z$'s by setting
\begin{equation}\label{e:self-loop}
W_{uu} = \frac{\big(2\sum\limits_{v \in i , u \ne v } W_{uv}\big) - \big(\frac{1}{n_i-2} \sum\limits_{v \in i, u \ne v } \sum\limits_{q\in i, q \ne u, v} W_{vq}\big)}{n_i-1}\,,
\end{equation}
where $n_i$ is the number of nodes in community $i$. Essentially, this takes twice the average within community edge weight for a particular node $u$ and subtracts out twice the average non-diagonal edge weight connecting the other nodes in community $i$. This approximates twice the impact node $u$ while negating the collective impact of the other nodes. 

\subsection{Handling repeated edge weight values}
In Section \ref{s:nmcomp}, when estimating the uniformly distributed values of $\widehat{G}(w)$ in (\ref{e:edge-ECDF}) and $\widehat{\Psi}_u^{(j)}$ in (\ref{e:psi-est}), we make an adjustment if we observe the same value multiple times. Let $W_1 < W_2$ be two consecutive sorted values of $W_{uv}$'s in the same pair of communities, with $\widehat{G}(W_1) = \frac{k}{n+1}$ and $\widehat{G}(W_2) = \frac{k+m}{n+1}$, where $m>1$. That is, there are $m$ different edges in the set $\{W_{uv}: \, u \in i, \, v \in j, \, W_{uv} =W_2\}$. Then, we set  
\begin{equation}\label{e:reposition}
\widehat{G}(W_2) = \frac{k + \frac{m}{2} + \frac{1}{2m}}{n+1}.
\end{equation}
The purpose of (\ref{e:reposition}) is to make all the new $\{\widehat{G}(W_{uv}): \, u \in i, \, v \in j, \, W_{uv} =W_2\}$ values to be at least halfway between $\widehat{G}(W_{1})$ and the original $\widehat{G}(W_2)$ value, but where more duplicate values will bring this value down further. 
When there are no duplicate values (i.e. $W_1 < W_2 < W_3<\ldots<W_{n}$), this formula leaves $\widehat{G}(W_2)$ intact. When duplicate values are present, there will still be duplicate values after applying (\ref{e:reposition}), but they are moved to a different location between $\frac{k}{n+1}$ and $\frac{k+m}{n+1}$.  

 If the $D_{j}(u)$ values repeat in (\ref{e:local-degree}), we make an analogous adjustment to (\ref{e:psi-est}).  







\section{Additional community detection details}\label{a:comm-det2}
In this appendix, more information is provided about community detection for the kinds of networks described in the main paper, including motivating ideas, algorithms, and some discussion of why existing community detection methods may be inappropriate in this setting.  


\subsection{Limitations of modularity}\label{mod-problems}
Many clustering algorithms on networks attempt to maximize modularity, but this may not be the appropriate measure to use for dense networks. To see this, consider the network shown in Figure \ref{f:reconstruction}. That figure begins on the left with a network where $\sigma =0$ with 74 nodes in 2 communities exhibiting within community positive association and between community negative association. Within community edges are drawn from a uniform distribution with a maximum of 150, while between community edges are drawn from a uniform distribution with a maximum of 100. Node sociability parameters for both communities range from .05 to .95 in increments of .025. The $H$-function in the original network is the same as that used for the rightmost plot in Figure \ref{f:different-h-funcs}, though since the between community edges have negative association, the inputs to that $H$-function are $1 - \Psi_u$ and $1 - \Psi_v$. While the network is not strictly assortative, after accounting for node sociabilities and using the appropriate inputs, edges between nodes in the same communities have 50\% greater weights than edges between nodes in different communities. In this sense there is some notion of homophily that is absent in other regimes where modularity fails. Ordered as in the figure, one can visually identify 2 distinct communities. The second plot from the left in the figure is the estimate of the left plot using the community assignments when clustering nodes using the walktrap algorithm of \cite{pons2006computing} with four steps, which returns three communities, not two. The third plot is the estimate of the first plot using the community assignments by calculating the leading non-negative eigenvector of the modularity matrix of the graph. The fourth plot is the estimate of the first plot using the correct assignments. The modularity of the true communities on this network is in fact slightly negative. This failure of community detection algorithms will feed incorrect community labels to the estimation procedure, leading to estimates that don't preserve the structure of the original network, as seen in the figure. 
\setcounter{figure}{9} 
\begin{figure}
\begin{center}
 \includegraphics[width = \linewidth]{positive_in_group_and_negative_outside_with_estimates.png}
\end{center}
\caption{Reconstructing a network (left) based on different community assignments.}
\label{f:reconstruction}
\end{figure}

The issue is modularity tries to identify highly interconnected nodes, where edge weights within communities are expected to be higher than edge weights between communities. The methodology described above needs clusters to have a different property to work appropriately, namely that the nodes in each community should obey a kind of monotonicity. Nodes in the same community should have similar patterns of connecting to nodes in other communities, and their own community. In the network in Figure \ref{f:reconstruction}, all nodes in community 1 ``prefer" other nodes in community 1 with high sociability ($\Psi$) values, but ``prefer" nodes in community 2 with low sociability values. This shared ordering of preferences over nodes in each community is crucial for ensuring that the estimation procedure will get appropriate orderings of local sociability statistics. 

It is not difficult to construct networks where it would be useful to combine the community detection using $L$ with modularity. For example, assume communities $i$ and $j$ have positive association both within community and between the communities, but the within community edges are generally much larger than the between community edges. Also assume both $i$ and $j$ have negative association with community $k$. Clustering nodes to simply share ordering of preferences would separate nodes in community $k$ but would not distinguish nodes in community $i$ from nodes in community $j$. Subsequently employing a modularity maximization algorithm on the estimated community consisting of nodes in $i$ and $j$ would recover the true communities and lead to better estimates for all of the edges within and between $i$ and $j$. 


\subsection{Greedy algorithm for community detection} \label{s:greedy-alg}
A direct algorithm for maximizing $L$ in (\ref{e:measure}) is to try to iteratively combine nodes into communities which will greedily make $L$ larger. It is computationally impractical to test every possible partition of nodes to maximize $L$. However, based on the structure of $L$, Algorithm \ref{alg:community-greedy} tries to combine communities which are most correlated with one another. 

At the start, each node is placed into its own estimated community. The ``aggregate degree" of each estimated community is defined as the total weight of the edges emanating from any node in the estimated community, where edges are double counted if they connect two nodes in the same estimated community. In Algorithm \ref{alg:community-greedy}, the aggregate degree of an estimated community is equal to the column sum of the estimated community's corresponding column in the ``communityAggregate" matrix. When each node is in its own estimated community, the aggregate degree of the community is the same as the degree of the node. In each ``round," the algorithm orders the estimated communities by their aggregate degrees at the beginning of the round from largest to smallest. Then the algorithm visits these estimated communities in order, and for each estimated community $\hat{i}$, selects a candidate estimated community $\hat{j} \ne \hat{i}$ which maximizes
$$\mathcal{C}_{\hat{i}\hat{j}} = corr\left\{\big(\sum\limits_{u:u \in \hat{i}}W_{uv}, \sum\limits_{q:q \in \hat{j}}W_{qv}\big): v \in V\right\},$$
where $\sum_{u:u \in \hat{i}}W_{uv},  v \in V$, is a length $n$ vector where each entry has the total edge weight connecting nodes in estimated community $\hat{i}$ with each node in the network. Estimated communities $\hat{i}$ and $\hat{j}$ are merged if doing so does not decrease the measure $L$. If $\hat{i}$ and $\hat{j}$ are merged, and if $\hat{j}$ has not already been visited by the algorithm this round, the algorithm will visit the combined $\hat{i}$ and $\hat{j}$ when it would have visited $\hat{j}$. Each round completes after the algorithm has completed all scheduled visits. 

If at least two estimated communities have been merged, the algorithm proceeds to the next round. If no communities have been merged, the algorithm does a sweep, calculating $\mathcal{C}_{\hat{i}\hat{j}}$ for all estimated communities, and attempts to merge pairs of estimated communities in decreasing order of $\mathcal{C}$ values. As soon as any pair of estimated communities are merged, the algorithm stops the sweep and proceeds to the next round.  Algorithm \ref{alg:community-greedy} terminates when it goes through a full round and a sweep without merging any estimated communities, or when the number of estimated communities reaches 1.        

Two details of the algorithm should be explained further. First, the algorithm does not require merging communities to increase $L$ because at the outset, when all nodes are in their own estimated community, combining two communities \textit{cannot} increase $L$. Therefore, requiring a merger to increase $L$ would prevent the algorithm from gaining any traction. Second, why go through each round instead of constantly sweeping, or better yet just combining the estimated communities which would most increase $L$? In addition to this proposal being computationally costly, it may also lead to the initial formation of a single large but overly heterogeneous community since communities can only contribute to $L$ once they contain three or more nodes. As this is already a greedy algorithm, we view our design as a means of not overlooking any estimated community, and as a conservative precaution against premature optimization. 

\begin{algorithm}
\SetAlgoLined
\DontPrintSemicolon 
\KwResult{$\{\, \widehat{i} \,\}$}
\KwInput{$W$}
\nl labels$^{(0)} = \vec{0}_n$\;
\nl labels$^{(1)} = \{1, ..., n\}$\; 
\nl $q$=1\;
\nl communityAggregate = $W$\;
\nl \While {$labels^{(q)} \ne labels^{(q-1)}$} {
	\nl $q$++\;
	\nl labels$^{(q)}$ = labels$^{(q-1)}$\;
	\nl labelOrder = sort labels by decreasing column sum of communityAggregate\; 
	\nl \While {$length(labelOrder) > 0$}{
		\nl labels$^{(q)}$, communityAggregate, labelOrder = \proc{attemptMerge(communityAggregate, labels$^{(q)}$, labelOrder)}\;
		\nl dequeue(labelOrder(1))\;
	}
	\nl \If{$labels^{(q)} = labels^{(q-1)}$}{
		\nl labels$^{(q)}$, communityAggregate = \proc{sweep(communityAggregate, labels$^{(q)}$)}
	}
}
\nl finalClustering = labels$^{(q)}$\;
\;
  \SetKwProg{myproc}{Procedure}{}{}
  \myproc{\proc{attemptMerge(\text{communityAggregate}, \text{labels}$^{\text{(q)}}$, \text{labelOrder})}}{
		\nl newLabels = labels$^{(q)}$\;
		\nl newOrder = labelOrder\;
		\nl A = newOrder(1)\;  %#reference community
		\nl B = $\argmax{\text{B}' \ne \text{A}}$ $corr([\text{communityAggregate(A)}], [\text{communityAggregate(B')}])$\;
 \tcp{[\text{communityAggregate(A)}] represents the Ath column in the communityAggregate matrix}
		\nl mergedLabels = newLabels\;
		\nl mergedLabels(mergedLabels = B) = A\;
		\nl mergedOrder = newOrder\; 
 		\nl mergedOrder(mergedOrder = B) = A\;
		\nl \If{calculateMeasureL(\text{W}, \text{mergedLabels}) $\ge$ calculateMeasureL(\text{W}, \text{newLabels})}{
			\nl newLabels = mergedLabels\;
			\nl newOrder = mergedOrder\;
			\nl [communityAggregate(A)] += [communityAggregate(B)]\;
			\nl Remove [communityAggregate(B)] column from communityAggregate\;
			}
		  \nl \KwRet newLabels, communityAggregate, 

}
\caption{Greedy algorithm using $L$}
\label{alg:community-greedy}
\end{algorithm}
 

\begin{algorithm}
\DontPrintSemicolon 
\SetAlgoLined
\setcounter{AlgoLine}{28}
	\SetKwProg{myproc}{Procedure}{}{}
	\myproc{\proc{sweep(communityAggregate, labels$^{(q)}$)}}{
		\nl correlationOrder = sort non-matching label pairs by decreasing correlations of columns in communityAggregate\;
		\nl \While{$length(correlationOrder) > 0$}{
			\nl newLabels = labels$^{(q)}$\;
			\nl A = correlationOrder(1, 1)\;  %#reference community
			\nl B = correlationOrder(1, 2)\;  %#candidate community for merger\;
			\nl mergedLabels = newLabels\; 
			\nl mergedLabels(mergedLabels = B) = A\;
			\nl \eIf{calculateMeasureL(W, mergedLabels) $\ge$ calculateMeasureL(W, newLabels)}{
				\nl newLabels = mergedLabels\;
				\nl [communityAggregate(A)] += [communityAggregate(B)]\;
				\nl Remove [communityAggregate(B)] column from communityAggregate\;
				\nl dequeueAll(correlationOrder)
			}{
			\nl dequeue(correlationOrder(1))
			}
		}
\nl \KwRet newLabels, communityAggregate
}
 \end{algorithm}


\subsection{Spatial clustering}\label{s:spatial}
For another perspective on community detection, we represent the original network as a matrix, and think of each row vector as a point in $n$-dimensional space corresponding to a particular node, where $n$ is the number of nodes in the network. Spatial clustering of the points then gives us communities, which is motivated as follows. Consider two nodes $u$ and $v$ in the same community with similar $\Psi$ values. For any third node $x \in j$, as $u \in i$ and $v \in i$, $H_{ij}(\Psi_u, \Psi_x)$ should be similar to $H_{ij}(\Psi_v, \Psi_x)$, thus $W_{ux}$ should be close to $W_{vx}$. Therefore, nodes in the same community with similar $\Psi$ values are expected to have similar edge weights, so entries in their corresponding rows should be similar in $n-2$ dimensions (the exceptions being due to no self loops). Therefore, in this $n$-dimensional space, two nodes with the same community and similar $\Psi$ values should be neighbors. 

Spatial clustering can be achieved using many existing clustering algorithms. For certain methods, such as those based on distances, one may need to account for the zeros on the diagonal by only calculating the distance between two nodes on the remaining $n-2$ dimensions. However, as with other spatial clustering problems, there is no algorithm which correctly clusters the nodes for every possible network.

\subsection{Spectral clustering algorithm incorporating $L$}\label{s:spectral-alg}
Appendix \ref{s:spatial} justifies the application of spatial clustering techniques for our networks of interest. In this section, we focus on one spatial clustering algorithm, spectral clustering, which we employ in tandem with the measure $L$ in (\ref{e:measure}), used to choose both a particular number of estimated communities as well as the best clustering for that number of communities. If the number of communities $K$ is known, we can follow the methodology of \cite{ng_jordan_weiss}. Measure the distance between every column in our network using a radial basis function kernel. Then construct a neighbors graph based on these distances, and a Laplacian based on this neighbors graph.  Finally, run K-means on the Laplacian to get clusterings. This can all be done using the existing \code{specc} function from the \proglang{R} package  \pkg{kernlab} just by specifying the number of centers, as \code{specc} will automatically select a scale parameter for the kernel. However, as this is not a deterministic algorithm, it can be helpful to run several replicates and take clustering that maximizes $L$.    

The remaining issue is how to choose the number of clusters, $K$, which is a priori unknown. However, we can use the introduced measure $L$ as a measure of clustering success, so we can impose a simple stopping rule, which is shown in Algorithm \ref{alg:community}. 

%We start at $K$*=1 and increment $K$* by 1. When $K >1$, K-means is not deterministic, so we run $specc$ several times for $K$* clusters, and keep only the clustering that performs best according to our measure for that value of $K$*. As soon as we see a drop in the measure from the best clustering with $K$* communities to the best clustering with $K$*+1 communities, we stop and take our final clustering to be the best clustering with $K$=$K$* communities.     

\begin{algorithm}
\SetAlgoLined
\KwResult{$\{\, \widehat{i} \,\}$}
\KwInput{$W$, replicates} \tcp{replicates is the number of times to run \code{specc} for each value of K}
\nl $\Delta$ = 1; $K$ =0\; 
	\nl \While {$\Delta > 0$}{
	\nl $K = K++$\;
	\nl clusterGoodness($K$) = $-\infty$\;
	\nl index $= 0$\;
    \nl \While{$\text{index}$ $<$ $\text{replicates}$}{
		\nl index$++$\;
		\nl candidateClustering = \code{specc}($W$, centers = $K$)\;
		\nl candidateGoodness =calculateMeasureL($W$, candidateClustering)\; 
			\nl \If{$candidateGoodness > clusterGoodness(K)$}{
    			\nl clusterGoodness($K$) = candidateGoodness\;
				\nl clustering($K$) = candidateClustering\;
			}
	}
\nl $\Delta$ = clusterGoodness($K$) $-$ clusterGoodness($K-1$) 
}
\nl finalClustering = clustering($K-1$)
 \caption{Community detection using spectral clustering and stopping criterion}
\label{alg:community}
 \end{algorithm}




\subsection{Clustering using normalized network}\label{s:rotated}

Another approach to community detection which seems to return interesting results is as follows: take the network as a whole and normalize each row in the original matrix by taking $N_{uv} = \frac{A_{uv}- \bar{A}_{u\bullet}}{SD(A_{u\bullet})}$, where $A_{u\bullet}$ represents the row in $A$ corresponding to node $u$. Next, calculate the eigenvalues $\lambda_1, \ldots, \lambda_n$ of the normalized network $N$, and note the index  corresponding to the largest difference in the absolute values of the real parts of these successive eigenvalues, $\argmax {1\le i \le n-1} |Re(\lambda_{i})| -  |Re(\lambda_{i+1})|$. Keep only the real parts of the first several eigenvectors corresponding to eigenvalues $\lambda_1, \ldots, \lambda_i$, and look for clusters in this lower dimensional space. 

A normalized version of the network presented in Figure \ref{f:reconstruction} is shown in the left plot of Figure \ref{f:normalized}. In the right plot of Figure \ref{f:normalized}, each point plots a row of the real parts of the first 2 eigenvectors of this row-normalized network, where the color of the point represents the true community of the corresponding node.

Plotting each row of these eigenvalues appears to give clearly separate clusters. Normalizing the network so all nodes have degree 1 doesn't give the same results. The first eigenvector of the original matrix often represents degree information, so the normalization to calculate $N$ should discard that degree information, but in the process, it also seems have some kind of downstream effect on other eigenvectors which helps them to capture the underlying communities. 
\begin{figure}
\begin{center}
    \includegraphics[width = 5in]{row_normalized_positive_in_group_and_negative_outside_combined2.png}
\end{center}
\caption{Normalized network with implied clustering.}
\label{f:normalized}
\end{figure}
This phenomenon is not  specific to the constructed network in Figure \ref{f:reconstruction}. Figure \ref{f:rot-eig} shows a network generated via an NSM with an $H$-function combining two exponential random variables, in the style of the third plot in Figure \ref{f:different-h-funcs}, where within community connections show positive association and between community connections show Simpson association, along with the eigenvectors from the normalized version of that network. Again, the real parts of the eigenvectors of the normalized network show a clear separating plane between the two communities, but it is not along one of the 2 dimensions, but rather along a combination of them.
\begin{figure}
\begin{center}
    \includegraphics[width = 4in]{exponential_NSM_with_neg_off_diag_combined3.png}
\end{center}
\caption{Network generated via NSM, with the clustering from its normalized version.}
\label{f:rot-eig}
\end{figure}

There is some intuition for why clustering based on normalizing the network might work. Ignoring the diagonal, normalizing the matrix makes each column sum to 0, and have variance 1. Since this is a dense weighted network, there is no true concept of a “hub,” but instead just a node that has greater edge weights. For this reason, we have no reason to value one node/dimension over another, and we can avoid having to find complicated kernels with different variances across different dimensions. If we don't normalize, the distance in only one dimension can totally dominate. In the case where all associations are positive, when doing community detection, we want to ignore degree in favor of preference, so it makes sense to divide by a standard deviation. In the case where some associations are negative (or Simpson), extreme $\Psi$ values will indicate greater variance, but again, we still want to ignore this for the purposes of community detection, so it still makes sense to normalize. The reason to normalize is since each node is connecting to all (or most) others, we want to ensure that we are accounting for a given node's average edge weight and the variance of the edge weights. In this way, we are still looking for patterns of preferences across other nodes. By normalizing the columns, we can compare them to one another on an apples-to-apples basis. A negative weight in the normalized matrix means that the edge between the reference node and another node is less than the average weight emanating from the reference node. For the purposes of community detection in this model, this is really what we care about, that is, patterns across nodes in a given community that show preferences across nodes throughout the whole network.

These last points also indicate the zeros on diagonal can bring up issues.  Imagine a situation where every other edge weight emanating from a node is extremely large but low variance. The 0 will shrink the average and increase the standard deviation by a lot, but it's purely artificial. If we use a distance between nodes/columns to cluster, it is imperative to ignore both the rows corresponding to those nodes. If there is actually a relatively large edge weight between them, that edge weight will be subtracted and squared twice. This is also purely artificial due to no self loops.
 

\section{Additional simulations}\label{more-sims}
\subsection{$H$-Normal LSM with $\sigma =0$} 
Figure \ref{f:no-noise-LSM} shows an $H$-Normal LSM with $\sigma =0$, where the value for within community edges are equal to $5+ 3\Phi_1^{-1}(\Psi_u) + 3\Phi_1^{-1}(\Psi_v)$. The between community edges, assuming $u$ is in community 1 and $v$ is in community 2, are equal to $8 - 3\Phi_1^{-1}(\Psi_u) + 1.5\Phi_1^{-1}(\Psi_v)$. The original network is shown on the left. In the middle is the reconstructed estimated network, which looks nearly identical to the original network. $\widehat{\sigma}$ is essentially 0 for this network, so the MSE (as defined in (\ref{e:MSE})) for each subnetwork is shown on the right. 
 \begin{figure}[h!]
 \includegraphics[width=\linewidth]{LSM_reconstruct2.png}
\caption{$H$-Normal LSM with $\sigma =0$, its estimate, and the subnetwork level MSE (\ref{e:MSE}).}
\label{f:no-noise-LSM}
\end{figure}

\subsection{Varied network}
Figure \ref{f:tricky} adds more complexity and departs even further from an $H$-Normal NSM, this time with 200 nodes and four communities of possibly different sizes. In this case, nodes are not assigned $\Psi$ values but i.i.d. Gamma(shape = 5, scale = 10) random values $\Gamma_u^{(1)}$ and $\Gamma_u^{(2)} = 1/\Gamma_u^{(1)}$.  However, within each community, $\Gamma_u^{(1)}$ and $\Gamma_u^{(2)}$ are normalized by dividing by $\sum_{u \in i} \Gamma_u^{(1)}$ and $\sum_{u \in i} \Gamma_u^{(2)}$, respectively. Each pair of communities is also assigned independently a random Gamma(shape =50$\sqrt{10}$, scale = 50$\sqrt{10}$) value $\Gamma_{ij}$. There is no mathematical significance to these parameters, they were chosen only to create a striking image. Within community edges are distributed according to a negative binomial distribution where the target number of successful trials is $2\times\Gamma_{ii}$ and the probability of success in each trial is $(1-\Gamma_u^{(1)})(1- \Gamma_v^{(1)})$. ``Adjacent" communities in the graph are distributed according to a negative binomial where the target number of successful trials is $1.5 \times\Gamma_{ij}$ and the probability of success in each trial is $(1-\Gamma_u^{(2)})(1- \Gamma_v^{(2)})$. Connections between communities 1 and 3 or between communities 2 and 4 are distributed according to a Poisson distribution with parameter $100\times\Gamma_{ij}\Gamma_u^{(1)}\Gamma_v^{(1)}$. Finally, connections between 1 and 4 are normally distributed with mean $10000\times(\Gamma_u^{(1)}+\Gamma_v^{(1)})$ and variance 1. Though the original network looks noisy, the estimate seems to capture a smooth approximation. In fact, in this network, $\hat{\sigma}$ is indistinguishable from 0 everywhere, so the MSE (\ref{e:MSE}) of each subnetwork, shown in the right plot of Figure \ref{f:tricky}, would be used for the bootstrap. 
\begin{figure}[h!]
 \includegraphics[width=\linewidth]{Tricky_example_reconstruction2.png}
\caption{From left to right: original, reordered, estimate, MSE (\ref{e:MSE}) of each subnetwork.}
\label{f:tricky}
\end{figure}


\subsection{Correlation matrix}
In this network, we first generate 3 ``lodestar" series $\mathcal{L}_1, \mathcal{L}_2, \mathcal{L}_3$ of length 1000 composed entirely of i.i.d. U$(0,1)$ random variables. We create a network with 200 nodes where the first 50 nodes each get a series which are (to different degrees) positively correlated with the first lodestar series, and negatively correlated with the second lodestar series. This is done by calculating the series for node $u$ at time $t$ as $$u(t) = \beta_1(u)\mathcal{L}_1(t) +  \beta_2(u)\mathcal{L}_2(t) +\beta_3(u)\mathcal{L}_3(t) + \big(1-\sqrt{\beta_1^2(u)+\beta_2^2(u) +\beta_3^2(u)}\big)\epsilon_u(t),$$ where $\epsilon_u(t) \sim$ U$(0, .35)$ is an idiosyncratic term for each node at each time step. For the first 50 nodes, the $\beta_1(u)$ values are positive and increasing with $u$, while the $\beta_2(u)$ values are negative and getting more negative with $u$, and $\beta_3(u) =0$. For the second 50 nodes, the $\beta_2(u)$ values are positive and increasing with $u$, while the $\beta_1(u)$ values are negative and getting more negative with $u$. The 101st to 150th nodes get series which are increasingly positively correlated to $\mathcal{L}_3$ series and increasingly negatively correlated to $\mathcal{L}_2$. Finally, the last 50 nodes have series which are increasingly negatively correlated to $\mathcal{L}_3$. Taking the correlations of the 200 series, we get the correlation matrix on the left of Figure \ref{f:corr-less}. Though the subnetwork between the first and fourth communities may look disordered, because nodes within each community are relatively correlated with each other, there is a discernible ordering in that between community subnetwork. Even so, from a practical perspective, edge weights only take on a narrow range of values near 0 in that subnetwork.  
\begin{figure}[h!]
 \includegraphics[width=\linewidth]{Correlation_mat_reconstruction_less_noise2.png}
\caption{Correlation network. From left to right: original, estimate, $\hat{\sigma}$ for each subnetwork, MSE (\ref{e:MSE}) of each subnetwork.}
\label{f:corr-less}
\end{figure}


\subsection{Edge weights with injected noise}\label{s:noisyedgesim}
In plot A of Figure \ref{f:noisyedgesim}, we start with the same underlying network as in Figure \ref{f:no-noise-NSM} but add Gaussian noise $\zeta_{uv}$ centered at 0 with a variance of 36 \textit{to the final edge weights}, not in ``normal" space, so the edge weights in the network can go below 0 and above 150. In plot B, the added Gaussian noise $\zeta_{uv}$ has variance 100 for within community edges and variance 225 for between community edges. In both cases, the reconstructed networks recover the underlying pattern, but the noisier network is estimated more coarsely.  
\begin{figure}[h!]
 \includegraphics[width=\linewidth]{HNSM_with_external_noise.png}
\caption{The network in plot A of Figure \ref{f:no-noise-NSM} with different values for $\sigma_{\zeta}^2$}
\label{f:noisyedgesim}
\end{figure}

Finally, in plot C, the same network is taken with $\sigma = .05$ everywhere and external noise is included by adding $\zeta_{uv}$ with variance 36 to the final edge weights. Even in this last case, the underlying signal is broadly recovered.    

\subsection{Multiplying sociability parameters}
It may be worth considering alternative models which can generate networks similar to LSMs and NSMs. For example, if edge weights are generated by multiplying sociability parameters, some surprising things happen. Ignoring diagonal zeros, if we take the values from .01 to 1 by .01 on each axis and let the value in the matrix equal the product of the axes, that results in the leftmost plot in Figure \ref{f:mult}.  In other words, if all values are positive, the result looks like the 4th plot in Figure \ref{f:different-h-funcs}. But if we center each axis to be mean 0 and recalculate, that produces the second plot in Figure \ref{f:mult}. This is not the intended setting for the models discussed in this paper, but if we were to treat these networks as such, in the first case, every node would be put into a single community. In the second case, those nodes with negative values would be put in a separate community from nodes with positive values. Using this split and reordering based on within community degree gives the rightmost plot in Figure \ref{f:mult}. 

In principle, were noise added to the edge weights of this network, knowing the true generating model type might improve estimation, as one may be able to smooth noise out over more observations by keeping all nodes in one large community. However, even using the ``wrong" communities, our estimation procedure appears to replicate the underlying network. Even though the generating process for this network is the same across both estimated communities, separating the second network into two communities is therefore a reasonable choice, especially since the two communities can be so easily defined. This kind of pattern only arises when multiplying nodes with positive sociabilities and others with negative sociabilities, not when all sociabilities have the same sign. While multiplying sociabilities hints at the idea of negative association, networks generated from these models are still restricted to symmetric contours of the type seen in the 4th plot of Figure \ref{f:different-h-funcs}. However, simply multiplying sociability parameters can’t generate networks that have contours of the type seen in the first 3 plots of Figure \ref{f:different-h-funcs}, nor can it give positive association patterns within a community but negative association patterns between communities when there are more than two communities.    
\begin{figure}
\begin{center}
 \includegraphics[width = \linewidth]{multiplicative_appendix.png}
\end{center}
\caption{From left to right: Network generated by multiplying sociability parameters. Network generated by demeaning sociability parameters then multiplying. Reordering of second network based on within estimated community degree.}
\label{f:mult}
\end{figure}   

\section{Higher dimensional $H$-functions and ``failure"}\label{s:Hext}
Thus far, we have defined $H$-functions as a class of functions that take in two uniform random variables and output another uniform random variable. This model can be extended to include a broader class of functions.   

\begin{definition}($d$-dimensional $H$-function.) A function $H:(0, 1)^d\rightarrow (0,1)$ is a $d$-dimensional $H$-function if the inputs are $d$ uniform random variables, and the output is a uniform random variable which is monotonic in each argument. 
\end{definition}

Note the $H$-Normal NSM (\ref{e:HNSM}) can be expressed in terms of a 3-dimensional $H$-function as:
\begin{equation}\label{e:HNSM2}  
W_{u v} = G^{-1}\left (H(\Psi_{u}, \Psi_{v}, \eta_{u v})\right)
\end{equation}
with $\eta_{uv} = \Phi_1(\epsilon_{uv}) \sim U(0,1)$, where
\begin{equation}\label{e:3DH} 
H(x, y, \eta) =\Phi_1\left(\frac{1}{\sqrt{1+\sigma^2}} \Phi_1^{-1}(h(x,y)) + \frac{\sigma}{\sqrt{1+\sigma^2}} \Phi_1^{-1}(\eta)\right),
\end{equation}
where $h(x,y)$ is a 2-dimensional $H$-function. The first two plots in Figure \ref{f:noise-failure} show the values along the $x$ and $y$ axes of functions which have the form of (\ref{e:3DH}) where $h(x,y)$ is in the schema of (\ref{e:HCDF}) such that $F_1$ and $F_2$ are Exponential distributions, and $F_{1, 2}$ is a Gamma distribution. As ``errors" are injected at each ($x$, $y$) pair, this can also be seen as a subnetwork generated via an $H$-Normal NSM where $H(x,y)$ in (\ref{e:HNSM}) is that of the third plot in Figure \ref{f:different-h-funcs}, and $\eta$ is injected at the edge level.

\begin{figure}
 \includegraphics[width=\linewidth]{additive_noise_and_failure_in_exponential_dist3.png}
\caption{Discretized 3-dimensional $H$-functions where axes show $x$ and $y$ but $\eta$ is drawn randomly at each ($x$, $y$) pair. The 2-dimensional $H$-function with arguments $x$ and $y$ is the same as in the third plot of Figure \ref{f:different-h-funcs}. From left to right: $H$-Normal NSM as in (\ref{e:HNSM}) with $\sigma = .05$. $H$-Normal NSM with $\sigma=.5$. Failure with $\alpha = .95$.  Failure with $\alpha = .5$.}
\label{f:noise-failure}
\end{figure}
The specific $H$-Normal NSM (\ref{e:H3-with-noise}) is a simple example which can be expressed in closed form as a 3-dimensional $H$-function as in (\ref{e:3DH}) where
\begin{equation}\label{e:3DH2}
H(x, y, \eta) =\Phi_1\left(\frac{1}{\sqrt{1+\sigma^2} \sqrt{1+\rho^2}} \Phi_1^{-1}(x) + \frac{\rho}{\sqrt{1+\sigma^2 }\sqrt{1+\rho^2}} \Phi_1^{-1}(y) + \frac{\sigma}{\sqrt{1+\sigma^2}} \Phi_1^{-1}(\eta)\right).
\end{equation}
In the network context, $\rho$ defines the relative influence of each of the node sociabilities, while  $\sigma$ controls the ``signal-to-noise" ratio of this 3-dimensional $H$-function. If one imagines observing several instances of the same network given by (\ref{e:HNSM2}) with idiosyncratic $\eta$ values, then increasing $\sigma$ would increase the variance of the individual edge weights from one instance to another. 
The 3-dimensional $H$-function (\ref{e:3DH2}) can be viewed as a composition of two 2-dimensional $H$-functions as follows:
$$ H(x, y, \eta) = H_{\sigma^2}(H_{\rho^2}(x,y), \eta), $$
where $H_{\rho^2}$ and $H_{\sigma^2}$ are both of the form (\ref{e:H-normal-rho}) but with different variance parameters, as indicated by their subscripts. 
All the observations about (\ref{e:3DH2}) do not depend on using 3-dimensional $H$-functions built from normal distributions or even $H$-Normal NSMs, but rather one can use any $d$-dimensional functions similar to (\ref{e:HCDF}) where $d>2$, with suitable adjustments based on the chosen distributions $F_1, F_2$, and $F_{1,2}$. 

One general method of creating higher dimensional $H$-functions is by chaining together lower dimensional $H$-functions as follows:
\begin{equation}\label{H-chain} 
H(x, y, \eta) = F_{1,2}\left(  F_1^{-1} \left( F_{3,4}(F_3^{-1}(x), F_4^{-1}(y)) \right), F_2^{-1}(\eta)   \right).
\end{equation}  
In general, the inner functions $F_3, F_4$ do not need to share the same form as the outer functions $F_1, F_2$. In the left plot of Figure \ref{f:Exp-Unif}, $F_3$, $F_4$ and $F_{3,4}$ correspond to the third plot in Figure \ref{f:different-h-funcs}, $\eta \sim$ U$(0,1)$ for each edge, and $F_1$, $F_2$ and $F_{1,2}$ correspond to the rightmost plot in Figure \ref{f:different-h-funcs}. The right plot of Figure \ref{f:Exp-Unif} swaps the roles of the inner and outer 2-dimensional $H$-functions in the left plot of Figure \ref{f:Exp-Unif}. While building higher dimensional $H$-functions in this way provides a lot of flexibility, this method may not give simple closed form expressions, and may not guarantee identifiability.  

\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth]{Exp_Unif_3d_H_functions.PNG}
\end{center}
\caption{3-dimensional $H$-functions of the form (\ref{H-chain}), where axes represent values of $x$ and $y$.}
\label{f:Exp-Unif}
\end{figure}


While the $H$-Normal NSM is a model with additive ``error," (\ref{e:HNSM2}) is more general. The following can be used to generate a different kind of ``error," one we shall call \textit{failure}. In this context, let 

\begin{equation}\label{e:failure}
H(x, y, \eta) = ( h(x, y))^\alpha \delta^{1-\alpha},
\end{equation}
where $\alpha$ is a value between 0 and 1, $h$ is a 2-dimensional $H$-function, and $\delta^{1-\alpha}$ is given by $F_\alpha^{-1}(\eta)$ where 
 $$F_\alpha(x) =
  \begin{cases}
 0, & x<0, \\
  \frac{x}{1-\alpha}, & 0 \le x \le 1-\alpha, \\
1, & x >  1-\alpha.
   \end{cases}
$$
Seen another way, 
$$\delta =
  \begin{cases}
 1 & \text{with probability } \alpha, \\
  \eta^{\frac{1}{1-\alpha}} & \text{with probability } 1-\alpha.
   \end{cases}
$$ 

One can see that (\ref{e:failure}) satisfies the definition of a 3-dimensional $H$-function by taking the Laplace-Stieltjes (LS) transform of the log of (\ref{e:failure}) with uniform inputs and recognizing that it matches the LS transform of the log of a uniform distribution. For $\eta \sim$ U$(0,1)$, the model (\ref{e:failure}) can also be written as
$$ H(x, y, \eta) \overset{d}{=}
  \begin{cases}
 h(x, y)^\alpha, & \text{with probability } \alpha, \\
   h(x, y)^\alpha \varepsilon, & \text{with probability } 1-\alpha,
   \end{cases}
$$ 
where $\varepsilon$ is also a uniform random variable. When $\alpha$ = 1, there is no ``error;" when $\alpha$ = 0, there is no degree correction, and increasing $\alpha$ increases the ``signal-to noise-ratio." As a contrast to the additive ``error" regime, in the failure case, when $\alpha$ is large, several different instances of the same network would share many of the \textit{exact} same edge weights, as in expectation, $100\times(1-\alpha)\%$ of the edge weights are given precisely by $h(x, y)^\alpha$.  The reason to call this kind of error ``failure" is that rather than defining a distribution that is concentrated near $h(x, y)$ with relatively small variation, even when $\alpha$ is large, there are infrequent occasions where the value will fall far below the modal value of $h(x, y)^\alpha$. The injected error will never raise the value greater than $h(x, y)^\alpha$, which also accounts for why the modal value lies at $h(x, y)^\alpha$ rather than $h(x, y)$. This kind of variation is reminiscent of each component in a system possessing a particular capacity, but occasional component failures cause that capacity to not be met. Figure \ref{f:noise-failure} depicts different levels of $\sigma$ and $\alpha$ being injected into the third plot of Figure \ref{f:different-h-funcs}. 

%\begin{figure}
% \includegraphics[width = \linewidth]{Collection_of_3d_H_functions.PNG}
%\caption{3-dimensional $H$-functions where axes represent values of $x$ and $y$, with the histogram of resulting values for each function below. The y-axes of the histograms are not aligned.}
%\label{f:too-many-plots}
%\end{figure}

%The $x$ and $y$ vectors used in Figure \ref{f:noise-failure} have the same length, although they can have different lengths $n_x$ and $n_y$ in general. When $\eta$ varies at the edge level, it must have length $n_x \times n_y$. This is sensible for incorporating idiosyncratic ``errors" into each edge in a network, but there are other use cases for higher dimensional $H$-functions. 
Consider the example of a road network, where the vertices are geographic locations, edges are roads, and edge weights are the number of cars that travel along the road each day. In this case, there should be degree correction, as there should be heavier traffic between certain locations than others. However, in addition to random variation for travelers along each road (which could be represented by additive ``error"), on some days, whether due to accidents, construction, or some other issue, the traffic along certain roads may fall dramatically. This latter case would be an example of failure. In this case, it may be better to model the network using a 4-dimensional $H$-function of the form
$$H(x,y, \eta_1, \eta_2) =   (F_{1,2,3}\left(F_1(x)+ F_2(y) + F_3(\eta_1)\right))^{\alpha} \delta_2^{1-\alpha},$$
which would incorporate additive ``error" through $\eta_1$, and failure through $\delta_2^{1-\alpha}$, which is a function of $\eta_2$. While this would be a valid 4-dimensional $H$-function in theory, in practice there may be issues arising from dependence between inputs, such as the probability of failure being correlated with the additive ``error."  Furthermore, one may want to include other covariates in the $H$-function which are not completely random, but rather systematic features, like time, that are different from node sociabilities, for use with tensors rather than matrices. There are many other kinds of $H$-functions that can incorporate various covariates and errors to model specific phenomena. Our goal, however, is not to catalog these possibilities, but to illustrate the richness and generality of the class of $H$-functions, particularly for generating random networks. 

\bibliography{yourbibfile}


\end{document}




