# Single-cell gene regulatory network analysis for mixed cell populations

Junjie Tang
School of Mathematical Sciences, Peking University
Center for Statistial Sciences, Peking University
and
Changhu Wang
School of Mathematical Sciences, Peking University
and
Feiyi Xiao
School of Mathematical Sciences, Peking University
and
Ruibin Xi
School of Mathematical Sciences, Peking University
Center for Statistial Sciences, Peking University  The authors gratefully
acknowledge the National Key Basic Research Project of China (2020YFE0204000),
the National Natural Science Foundation of China (11971039), and Sino-Russia
Mathematics Center.

###### Abstract

Gene regulatory network (GRN) refers to the complex network formed by
regulatory interactions between genes in living cells. In this paper, we
consider inferring GRNs in single cells based on single cell RNA sequencing
(scRNA-seq) data. In scRNA-seq, single cells are often profiled from mixed
populations and their cell identities are unknown. A common practice for
single cell GRN analysis is to first cluster the cells and infer GRNs for
every cluster separately. However, this two-step procedure ignores uncertainty
in the clustering step and thus could lead to inaccurate estimation of the
networks. To address this problem, we propose to model scRNA-seq by the
mixture multivariate Poisson log-normal (MPLN) distribution. The precision
matrices of the MPLN are the GRNs of different cell types and can be jointly
estimated by maximizing MPLN’s lasso-penalized log-likelihood. We show that
the MPLN model is identifiable and the resulting penalized log-likelihood
estimator is consistent. To avoid the intractable optimization of the MPLN’s
log-likelihood, we develop an algorithm called VMPLN based on the variational
inference method. Comprehensive simulation and real scRNA-seq data analyses
reveal that VMPLN performs better than the state-of-the-art single cell GRN
methods.

Keywords: Gene regulatory network; Graphical model; Precision matrix;
Variational inference; Single cell RNA sequencing; COVID-19.

##  1 Introduction

Gene regulatory network (GRN), representing the regulatory relationships
between genes, is important for understanding the complex biological system
(Arendt et al., 2016). GRNs can be inferred based on gene expression data such
as RNA sequencing (RNA-seq) data. Bulk expression data are most commonly used
for GRN inference and numerous methods have been developed (see Marbach et al.
(2012) and references therein). However, bulk data are profiled from pooled
cell populations and thus can only provide average expressions of many cells.
The recent development of single-cell RNA-seq (scRNA-seq) technologies can
measure gene expression at the single-cell level (Gohil et al., 2021; Nam et
al., 2021), thus offering unprecedented opportunity for single-cell GRN
inference.

To account for the unique features of scRNA-seq data, a number of GRN
inference methods based on scRNA-seq data have been developed (Aibar et al.,
2017; Specht and Li, 2017; Chan et al., 2017). These methods often make the
implicit assumption that all cells share the same GRN. However, single cells
in scRNA-seq data usually belong to multiple cell types and each cell type has
its own specific GRN and expression pattern. The cell identities are unknown
and have to be determined using scRNA-seq data. To infer GRNs of different
cell types, one has to first assign single cells to different cell types (e.g.
by clustering) and then estimate the GRNs using available methods. This two-
step procedure can provide accurate GRN estimation if different cell types are
well separated. If, instead, different cell types have a higher mixing degree,
a large proportion of cells cannot be confidently assigned to a cell type and
the ambiguity of the cell type assignment could seriously influence the
performance of GRN inference.

Here, we consider developing a GRN inference method based on scRNA-seq data
with mixed cell populations using mixture models. One major advantage of
mixture models is that cell types need not be predetermined before GRN
inference. Instead, mixture models can allow joint analyses of clustering and
GRN inference, and thus could give better GRN estimation when different cell
types are poorly separated. Since scRNA-seq data are count data, often rather
small count data with many zeros (the high dropout problem), graphical models
for count data would be more suitable than the widely used Gaussian graphical
model (GGM) (Meinshausen and Bühlmann, 2006; Friedman et al., 2008) or its
mixture version. Available graphical models for count data include Poisson
graphical models (PGMs) (Yang et al., 2012; Allen and Liu, 2013) and Poisson
log-normal (PLN) models (Wu et al., 2018; Chiquet et al., 2019; Silva et al.,
2019). Compared with PGMs, the PLN models can model the over-dispersion
commonly observed in scRNA-seq data (Ziegenhain et al., 2017). Therefore, we
propose to use the mixture PLN (MPLN) model for GRN inference in single cells.

A non-negative integer random vector
$\mathbf{Y}=(Y_{1},\cdots,Y_{p})^{T}\in\mathbf{R}^{p}$ follows a PLN
distribution, if conditional on a latent random vector
$\mathbf{X}=(X_{1},\cdots,X_{p})^{T}\in\mathbf{R}^{p}$ with
$\mathbf{X}\sim\mbox{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, each element
of $\mathbf{Y}$ independently follows the univariate Poisson distribution,
i.e. $Y_{j}\sim\mbox{Poisson}(\exp(X_{j}))$ ($j=1,2,\cdots,p$). Similar to the
GGM, the network of the PLN model is the precision matrix
$\boldsymbol{\Theta}={\boldsymbol{\Sigma}}^{-1}$ of the latent variable
$\mathbf{X}$. The MPLN model is a mixture of $G$ different PLN models. The
precision matrix of each component of the MPLN model represents the network of
a cell type. Assuming that the networks are sparse, we can maximize the lasso-
penalized log-likelihood of the MPLN model to estimate the networks.

We first establish the basic properties of the MPLN model: the MPLN model is
identifiable, and its Fisher information matrix is positive definite. We
further show that the network estimator by maximizing the lasso-penalized log-
likelihood of the MPLN model is consistent. As a special case, this result
also establishes the consistency result for the PLN model, which has been
lacking in the literature. Directly maximizing the lasso-penalized log-
likelihood of the MPLN is computationally intractable. We adopt the
variational inference approach (Jordan et al., 1999; Wainwright et al., 2008)
and develop an algorithm called variational mixture Poisson log-normal (VMPLN)
for simultaneous analyses of clustering and network inference.

We compare VMPLN with popular graphical methods and state-of-the-art single
cell regulatory network inference methods. Comprehensive simulation shows that
VMPLN achieves better performance especially in the scenarios that different
cell types have a high mixing degree. Benchmarking on real scRNA-seq data also
demonstrates that VMPLN can provide more accurate network estimation in most
cases. Finally, we apply VMPLN to a large scRNA-seq dataset from patients
infected with severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and
find that VMPLN identifies critical differences of regulatory networks in
immune cells between patients with moderate and severe symptoms.

The paper is organized as follows. Section 2 presents the MPLN model and its
theoretical properties. Section 3 derives the VMPLN algorithm. Simulation and
real data analyses are in Section 4 and 5, respectively. All proofs of the
theoretical properties are in supplementary material. All data used and
corresponding source codes are available and can be accessed at
https://github.com/XiDsLab/scGeneNet.

##  2 Model and Theoretical properties

###  2.1 The Mixture Poisson log-normal model for scRNA-seq data

Suppose that a scRNA-seq dataset consists of $n$ cells and $p$ genes. Let
$\mathbf{Y}_{i}=\left(Y_{i1},\cdots,Y_{ip}\right)^{T}$ be the observed
expression vector of the $i$th cell, where $Y_{ij}$’s are all non-negative
integers. Let $\mathbf{Y}=(\mathbf{Y}_{1},\cdots,\mathbf{Y}_{n})^{T}$ be the
observed count matrix. Single cells in scRNA-seq data belong to $G$ different
cell types and each cell type has its own unique mean gene expression and
regulatory network. The cell identities are unknown and have to be determined
based on the observed data $\mathbf{Y}$. The observed expression $Y_{ij}$ is a
noisy measurement of the true expression $\exp\left(X_{ij}\right)$ of the
$i$th cell at the $j$th gene. Conditional on $X_{ij}$, we assume that $Y_{ij}$
follows a Poisson distribution with a mean
$\lambda_{ij}=l_{i}\exp\left(X_{ij}\right)$, where $l_{i}$ is the library size
of the $i$th cell and can be readily estimated using available methods
(Hafemeister and Satija, 2019; Lun et al., 2016). Denote
$\mathbf{X}_{i}=\left(X_{i1},\cdots,X_{ip}\right)^{T}$ be the logarithm of the
true expression vector of the $i$th cell. We assume that the logarithm of
expressions of single cells in the $g$th cell type are normally distributed
with a mean $\boldsymbol{\mu}_{g}$ and a covariance
$\boldsymbol{\Theta}_{g}^{-1}$. Thus, given the underlying cell type $Z_{i}=g$
($g=1,\cdots,G$) of the $i$th cell, the conditional distribution of
$\mathbf{X}_{i}$ is
$\mbox{N}\left(\boldsymbol{\mu}_{g},{\boldsymbol{\Theta}_{g}}^{-1}\right)$. We
further assume that $Z_{i}$ follows a multinomial distribution
$\mbox{Multinomial}(1,\boldsymbol{\pi})$, where
$\boldsymbol{\pi}=(\pi_{1},\cdots,\pi_{G})^{T}$ is the proportion parameter
representing the composition of cell types. In summary, we have the following
MPLN model

$\displaystyle\mathbf{Y}_{i}|\mathbf{X}_{i}$
$\displaystyle\sim\prod_{j=1}^{p}\mbox{Poisson}\left[l_{i}\
\exp\left(X_{ij}\right)\right],$ (1) $\displaystyle\mathbf{X}_{i}|Z_{i}=g$
$\displaystyle\sim\mbox{N}\left(\boldsymbol{\mu}_{g},{\boldsymbol{\Theta}_{g}}^{-1}\right),\boldsymbol{\Theta}_{g}\succ
0,$ $\displaystyle Z_{i}$
$\displaystyle\sim\mbox{Multinomial}(1,\boldsymbol{\pi}),$

where $\boldsymbol{\Theta}_{g}\succ 0$ means that $\boldsymbol{\Theta}_{g}$ is
positive definite.

Denote
$\boldsymbol{\theta}=\big{(}\boldsymbol{\pi},\boldsymbol{\mu}=\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G},\boldsymbol{\Theta}=\\{\boldsymbol{\Theta}_{g}\\}_{g=1}^{G}\big{)}$
as the set of unknown model parameters. Let
$p(\mathbf{Y}_{i}|\mathbf{X}_{i})=\prod_{j=1}^{p}\left\\{\left[l_{i}\exp\left(X_{ij}\right)\right]^{Y_{ij}}\exp\left[-l_{i}\exp\left(X_{ij}\right)\right](Y_{ij}!)^{-1}\right\\}$
be the conditional probability mass function of $\mathbf{Y}_{i}$ given
$\mathbf{X}_{i}$. Suppose that
$p(\mathbf{X};\boldsymbol{\mu}_{g},\boldsymbol{\Theta}_{g})$ is the density
function of the normal distribution with mean $\boldsymbol{\mu}_{g}$ and
covariance $\boldsymbol{\Theta}_{g}^{-1}$,

$p(\mathbf{X};\boldsymbol{\mu}_{g},\boldsymbol{\Theta}_{g})=\left(2\pi\right)^{-p/2}|\boldsymbol{\Theta}_{g}|^{1/2}\exp\left[-\frac{1}{2}\left(\mathbf{X}-\boldsymbol{\mu}_{g}\right)^{T}\boldsymbol{\Theta}_{g}\left(\mathbf{X}-\boldsymbol{\mu}_{g}\right)\right].$

The conditional density function
$p(\mathbf{X}_{i}|Z_{i};\boldsymbol{\mu},\boldsymbol{\Theta})$ of
$\mathbf{X}_{i}$ given $Z_{i}$ can be written as
$p(\mathbf{X}_{i}|Z_{i};\boldsymbol{\mu},\boldsymbol{\Theta})=\prod_{g=1}^{G}\left[p(\mathbf{X}_{i};\boldsymbol{\mu}_{g},\boldsymbol{\Theta}_{g})\right]^{I\left(Z_{i}=g\right)}$.
Denote $p(Z_{i};\boldsymbol{\pi})=\prod_{g=1}^{G}\pi_{g}^{I(Z_{i}=g)}$ as the
probability mass function of the multinomial distribution. The log-likelihood
of the MPLN model (1) is

$\ell_{n}(\boldsymbol{\theta})=\sum_{i=1}^{n}\log\left(p(\mathbf{Y}_{i};\boldsymbol{\theta})\right)=\sum_{i=1}^{n}\log\iint
p(\mathbf{Y}_{i}|\mathbf{X}_{i})p(\mathbf{X}_{i}|Z_{i};\boldsymbol{\mu},\boldsymbol{\Theta})p(Z_{i};\boldsymbol{\pi})\mathrm{d}\mathbf{\mathbf{X}_{i}}\mathrm{d}Z_{i}\
,$ (2)

where $p(\mathbf{Y}_{i};\boldsymbol{\theta})$ is the marginal probability mass
function of $\mathbf{Y}_{i}$. Similar to the GGM, the precision matrix
$\boldsymbol{\Theta}_{g}$ of the $g$th cell type represents the cell-type-
specific regulatory network. Let $\Theta_{g,lm}$ be the element of the $l$th
row and the $m$th column of $\boldsymbol{\Theta}_{g}$. The regulatory networks
are sparse and can be estimated by minimizing the lasso-penalized negative
log-likelihood

$-n^{-1}\ell_{n}(\boldsymbol{\theta})+\lambda_{n}\sum_{g=1}^{G}\|{\boldsymbol{\Theta}_{g}}\|_{1,\text{
off }},$ (3)

where $\lambda_{n}>0$ is a tuning parameter and
$\|{\boldsymbol{\Theta}_{g}}\|_{1,\text{ off }}=\sum_{l\neq
m}|{\Theta_{g,lm}}|$ is the off-diagonal $l_{1}$-norm of the matrix
$\boldsymbol{\Theta}_{g}$. In the following, we first establish the
consistency of the estimator obtained by minimizing the objective function (3)
and then derive an algorithm for estimating the precision matrices based on
the variational inference method.

###  2.2 Theoretical properties

In this section, we always assume that the true means
$\boldsymbol{\mu}_{g}^{*}$ and proportions $\pi_{g}^{*}$ ($g=1,\cdots,G$) are
known. Let $\boldsymbol{\nu}_{g}={\rm vech}(\boldsymbol{\Theta}_{g})$ be the
vectorization of the precision matrix $\boldsymbol{\Theta}_{g}$ and
${\boldsymbol{\nu}}=(\boldsymbol{\nu}_{1}^{\rm
T},\ldots,\boldsymbol{\nu}_{G}^{\rm T})^{\rm T}$ (see supplementary material
for the exact definition of $\rm vech(\cdot)$). In this case, the log-
likelihood $\ell_{n}(\boldsymbol{\theta})$ can be viewed as a function of
${\boldsymbol{\nu}}$, also denoting as $\ell_{n}(\boldsymbol{\nu})$, and we
consider the estimator $\hat{\boldsymbol{\nu}}_{n}$ that minimizes
$-n^{-1}\ell_{n}(\boldsymbol{\nu})+\lambda_{n}\sum_{g=1}^{G}\|{\boldsymbol{\Theta}_{g}}\|_{1,\text{
off }}$ subject to $\boldsymbol{\Theta}_{g}\succ 0\ (g=1,\cdots,G)$. Denote
$\boldsymbol{\nu}^{*}$ as the true value of the unknown parameter
$\boldsymbol{\nu}$, $S(\boldsymbol{\nu})=\\{i|\ \boldsymbol{\nu}_{i}\neq 0\\}$
as the support of $\boldsymbol{\nu}$, and $S^{*}=S({\boldsymbol{\nu}}^{*})$.
Suppose that $\mathbf{y}$ follows the MPLN model with its log-likelihood
function $\ell(\boldsymbol{\nu},\mathbf{y})$. Denote
$\boldsymbol{\Gamma}^{*}=\boldsymbol{\Gamma}(\boldsymbol{\nu}^{*})=-{\rm
E}\left[\frac{\partial^{2}\ell(\boldsymbol{\nu},\mathbf{y})}{\partial\boldsymbol{\nu}\partial\boldsymbol{\nu}^{T}}\right]\bigg{|}_{\boldsymbol{\nu}=\boldsymbol{\nu}^{*}}$
as the Fisher information matrix of the MPLN at $\boldsymbol{\nu}^{*}$, and
${\boldsymbol{\Gamma}^{*}}_{T_{1}T_{2}}$ as the submatrix of
${\boldsymbol{\Gamma}^{*}}$ with rows and columns index by sets $T_{1}$ and
$T_{2}$, respectively. Before presenting the theoretical properties, we give
the following conditions:

  * (C1)

The eigenvalues of the precision matrices are bounded away from zero and
infinity, i.e. there are two constants $0<m<M$ such that
$m\leq\lambda_{min}(\boldsymbol{\Theta}_{g})\leq\lambda_{max}(\boldsymbol{\Theta}_{g})\leq
M$ for $g=1,\cdots,G$, where $\lambda_{min}(\boldsymbol{\Theta}_{g})$ and
$\lambda_{max}(\boldsymbol{\Theta}_{g})$ are the minimum and maximum
eigenvalue of the precision matrix $\boldsymbol{\Theta}_{g}$.

  * (C2)

The library sizes $l_{i}>0$ ($i=1,\cdots,n$) are independent and identically
distributed random variables with a bounded support.

  * (C3)

The true mean vectors $\boldsymbol{\mu}_{g}^{*}$ ($g=1,\cdots,G$) are bounded
and different from each other.

  * (C4)

The irrepresentability condition:
$||\boldsymbol{\Gamma}^{*}_{{S^{*}}^{c}S^{*}}{(\boldsymbol{\Gamma}^{*}_{S^{*}S^{*}})}^{-1}||_{\infty}<1$.

Define $\mathcal{D}=\\{{\boldsymbol{\nu}}=(\boldsymbol{\nu}_{1}^{\rm
T},\boldsymbol{\nu}_{2}^{\rm T},\ldots,\boldsymbol{\nu}_{G}^{\rm T})^{\rm
T}|~{}\boldsymbol{\nu}_{g}={\rm
vech}(\boldsymbol{\Theta}_{g}),m\leq\lambda_{min}(\boldsymbol{\Theta}_{g})\leq\lambda_{max}(\boldsymbol{\Theta}_{g})\leq
M\\}$ and $\kappa=\lambda_{min}(\boldsymbol{\Gamma}^{*})$ be the minimum
eigenvalue of the Fisher information matrix at $\boldsymbol{\nu}^{*}$.
Condition (C1-C2) are commonly used in the literature (Cai et al., 2011; Li et
al., 2020). Condition (C3) is to ensure that different components of the MPLN
model can be distinguished from each other. Under Condition (C3), the MPLN
model is identifiable and its Fisher information $\boldsymbol{\Gamma}^{*}$ is
positive definite, and thus $\kappa>0$. The irrepresentability condition (C4)
is also commonly used (Zhao and Yu, 2006; Ravikumar et al., 2011). For a GGM
with a covariance $\boldsymbol{\Sigma}$ and a known mean, its Fisher
information matrix is $\boldsymbol{\Sigma}\otimes\boldsymbol{\Sigma}$, where
$\otimes$ represents the Kronecker product. There is an intuitive explanation
for the irrepresentability condition of the GGM (Ravikumar et al., 2011).
However, the Fisher information matrix $\boldsymbol{\Gamma}^{*}$ of the MPLN
model has no closed form and we do not have an intuitive explanation for the
irrepresentability condition (C4). Based on the above conditions, we present
the theoretical properties of the MPLN and the estimator
$\hat{\boldsymbol{\nu}}_{n}$ in the following theorems.

######  Theorem 1.

Under Condition (C1-C3), the MPLN model is identifiable, and its Fisher
information matrix $\boldsymbol{\Gamma}^{*}$ at $\boldsymbol{\nu}^{*}$ is
positive definite.

Theorem 1 establishes basic properties of the MPLN model and ensures that the
MPLN model is well-behaved under rather mild conditions. The proof of this
theorem is nontrivial because the PLN distribution has no finite moment
generating function and its density function is rather complex. However, its
moments are finite and have closed forms. We use its moments to prove Theorem
1. For the identifiability, the basic idea of the proof is that
identifiability of a mixture model is equivalent to linear independence of its
components. Using moments of the PLN distribution, we can show that only a
zero vector can make the linear combination of the components of the MPLN
model as zero. To prove the positive definiteness of the Fisher information,
we also use the moments of the PLN distribution and convert the problem to
showing that a set of equations only have zero solutions. Based on this
result, we can further prove the consistency and the sign consistency of the
estimator $\hat{\boldsymbol{\nu}}_{n}$.

######  Theorem 2.

Under Conditions (C1-C3), we have

$\mathbb{P}\left[||\hat{\boldsymbol{\nu}}_{n}-\boldsymbol{\nu}^{*}||_{2}\leq\frac{3}{\kappa}\sqrt{\frac{G{p(p+1)}}{{2}}}\big{(}n^{-1}||\nabla{\ell}_{n}(\boldsymbol{\nu}^{*})||_{\infty}+2{\lambda_{n}}\big{)}\right]\rightarrow
1,\mbox{ as }n\rightarrow\infty,$

where $\nabla{\ell}_{n}(\boldsymbol{\nu}^{*})$ is the gradient of the log-
likelihood $\ell_{n}(\boldsymbol{\theta})$ at $\boldsymbol{\nu}^{*}$ and
$\lambda_{n}>0$ is the regularization parameter.

######  Theorem 3.

Under the assumptions of Theorem 2 and Condition (C4), choosing
$\lambda_{n}>0$ such that $\lambda_{n}\rightarrow 0$ and
${\sqrt{n}\lambda_{n}}\rightarrow\infty$, we have

$\mathbb{P}\left[S(\hat{\boldsymbol{\nu}}_{n})=S^{*}\right]\rightarrow 1\mbox{
as }n\rightarrow\infty.$

Theorem 3 says that if we choose $\lambda_{n}>0$ such that it does not
converge to 0 too fast, $\hat{\boldsymbol{\nu}}_{n}$ can consistently recover
the nonzero elements of $\boldsymbol{\nu}^{*}$. For example, we can choose
$\lambda_{n}=\sqrt{\log n/n}$. The consistency theory established in this
paper is only for fixed dimension $p$. Currently, we cannot prove consistency
results for the high dimensional setting. One major difficulty is that, we can
prove that the expectation of the negative log-likelihood function
$-n^{-1}{\ell}_{n}(\boldsymbol{\nu})$ at $\boldsymbol{\nu}^{*}$ is
strong-$\kappa$ convex, but we can not characterize how $\kappa$ changes as
$p\rightarrow\infty$. The PLN graphical model corresponds to the case of
$G=1$. Theorem 2 and 3 imply that, for the PLN graphical model, the network
estimated by minimizing the lasso-penalized negative log-likelihood is a
consistent estimator.

##  3 Algorithm

###  3.1 Variational inference for the MPLN model

The log-likelihood (2) of the MPLN model involves an intractable integration
and thus directly minimizing (3) is computationally very difficult. We
therefore adopt the variational inference approach to estimate the networks
(Jordan et al., 1999; Wainwright et al., 2008). We approximate the log-
likelihood by the evidence low bound (ELBO)
$\ell_{\text{E}}\left(\boldsymbol{\eta},\boldsymbol{\theta}\right)$ and
estimate $\boldsymbol{\theta}$ by minimizing
$-\ell_{\text{E}}\left(\boldsymbol{\theta},\boldsymbol{\eta}\right)+\lambda_{n}\sum_{g=1}^{G}\left\|{\boldsymbol{\Theta}_{g}}\right\|_{{1,\text{
off }}}$, where $\boldsymbol{\eta}\in\mathscr{H}$ is the parameter of the
variational distribution family
$\mathscr{L}=\\{q(\mathbf{X},\mathbf{Z};\boldsymbol{\eta}):\boldsymbol{\eta}\in\mathscr{H}\\}$.
For $\boldsymbol{\eta}\in\mathscr{H}$, the ELBO
$\ell_{\text{E}}\left(\boldsymbol{\eta},\boldsymbol{\theta}\right)$ is defined
as
$\ell_{\text{E}}\left(\boldsymbol{\eta},\boldsymbol{\theta}\right)={\mbox{E}}_{q(\mathbf{\mathbf{X},\mathbf{Z}};\boldsymbol{\eta})}\left[\log
p(\mathbf{\mathbf{X},\mathbf{Z}},\mathbf{\mathbf{Y}};\boldsymbol{\theta})-\log
q(\mathbf{\mathbf{X},\mathbf{Z}};\boldsymbol{\eta})\right]$.

For computational considerations, we consider the following variational
distribution family. Conditional on $Z_{i}=g$, this distribution family
assumes that $X_{ij}$ ($j=1,\cdots,p$) are independent normal variables with a
mean $M_{g,ij}$ and a variance $S_{g,ij}$. The distribution of $Z_{i}$ is a
multinomial distribution with probabilities
$\mathbf{P}_{i}=(P_{i1},\cdots,P_{iG})^{T}$ as proportion parameters. Denote
$\mathbf{M}_{g}=(M_{g,ij})_{n\times p}$, $\mathbf{S}_{g}=(S_{g,ij})_{n\times
p}$ and $\mathbf{P}=(\mathbf{P}_{1},\cdots,\mathbf{P}_{n})^{T}$. The
variational parameters are
$\boldsymbol{\eta}=\left\\{\boldsymbol{\eta}_{g}\right\\}_{g=1}^{G}=\left\\{\mathbf{M}_{g},\mathbf{S}_{g},\\{\mathbf{P}_{ig}\\}_{i=1}^{n}\right\\}_{g=1}^{G}$
with $\boldsymbol{\eta}\in\mathscr{H}=\left\\{\boldsymbol{\eta}\big{|}\
S_{g,ij}>0,~{}P_{ig}\geq 0,~{}\sum_{g=1}^{G}P_{ig}=1\right\\}$. Thus, the
uncentered variational distribution family is

$\mathscr{Q}=\left\\{q(\mathbf{\mathbf{X},\mathbf{Z}};\boldsymbol{\eta})=\prod_{i=1}^{n}\left[q(Z_{i};\mathbf{P}_{i})\prod_{j=1}^{p}\prod_{g=1}^{G}q\big{(}X_{ij}|Z_{i}=g;M_{g,ij},S_{g,ij}\big{)}\right],\
\boldsymbol{\eta}\in\mathscr{H}\right\\},$

where $q\big{(}X_{ij}|Z_{i}=g;M_{g,ij},S_{g,ij}\big{)}$ is the density
function of $N(M_{g,ij},S_{g,ij})$ and $q(Z_{i};\mathbf{P}_{i})$ is the
density function of $\mbox{Multinomial}(1,\mathbf{P}_{i})$.

Given two matrices $\mathbf{A}$ and $\mathbf{B}$, let
$\mathbf{A}\odot\mathbf{B}$ be the Hadamard product of $\mathbf{A}$ and
$\mathbf{B}$, $\mathbf{A}_{i\cdot}$ be the vector of the $i$th row of
$\mathbf{A}$ and $\mathbf{A}_{\cdot j}$ be the vector of the $j$th column of
$\mathbf{A}$. Given a vector $\mathbf{c}$, define $\mathbf{D}(\mathbf{c})$ as
the diagonal matrix whose diagonal elements are $\mathbf{c}$. Denote
$\boldsymbol{l}=\left(l_{1},l_{2},\cdots,l_{n}\right)^{T}$,
$\boldsymbol{\Sigma}_{gi}=(\mathbf{M}_{g,i\cdot}-\boldsymbol{\mu}_{g}){(\mathbf{M}_{g,i\cdot}-\boldsymbol{\mu}_{g})}^{T}+\mathbf{D}\left(\mathbf{S}_{g,i\cdot}\right)$,
$F_{1}\left(l_{i},M_{g,ij},S_{g,ij}\right)=\exp\left(M_{g,ij}+2^{-1}S_{g,ij}+\log
l_{i}\right)$,
$F_{2}\left(\boldsymbol{\Theta}_{g},\mathbf{M}_{g,i\cdot},\mathbf{S}_{g,i\cdot},\boldsymbol{\mu}_{g}\right)=2^{-1}\left[\log\det\boldsymbol{\Theta}_{g}-\operatorname{tr}(\boldsymbol{\Theta}_{g}\boldsymbol{\Sigma}_{gi})\right]$
and
$\boldsymbol{\theta}_{g}=\left(\boldsymbol{\pi}_{g},\boldsymbol{\mu}_{g},\boldsymbol{\Theta}_{g}\right)$
as the set of the unknown model parameters of the $g$th cell type. With the
above variational distribution family $\mathscr{L}$, the ELBO can be written
as
$\ell_{\text{E}}\left(\boldsymbol{\eta},\boldsymbol{\theta}\right)=\sum_{g=1}^{G}\ell_{\text{E}}^{(g)}\left(\boldsymbol{\eta}_{g},\boldsymbol{\theta}_{g}\right)$
with

$\ell_{\text{E}}^{(g)}\left(\boldsymbol{\eta}_{g},\boldsymbol{\theta}_{g}\right)=\mathbf{P}_{\cdot
g}^{T}\left(\boldsymbol{\Lambda}_{g}^{(1)}-\boldsymbol{\Lambda}_{g}^{(2)}+\boldsymbol{\Lambda}_{g}^{(3)}\right)\mathbf{1}_{p}+\mathbf{P}_{\cdot
g}^{T}\left(\boldsymbol{\Lambda}_{g}^{(4)}+\boldsymbol{\Lambda}_{g}^{(5)}\right)+K_{g}\left(\mathbf{Y}\right),$

where $\boldsymbol{\Lambda}_{g}^{(1)}=\mathbf{Y}\odot\mathbf{M}_{g}$,
$\boldsymbol{\Lambda}_{g}^{(2)}=\left(\boldsymbol{\Lambda}_{g,ij}^{(2)}\right)_{n\times
p}=\big{(}F_{1}\left(l_{i},M_{g,ij},S_{g,ij}\right)\big{)}_{n\times p}$,
$\boldsymbol{\Lambda}_{g}^{(3)}=2^{-1}\log\mathbf{S}_{g}$,
$\boldsymbol{\Lambda}_{g}^{(4)}=\log(\pi_{g})\mathbf{1}_{n}-\log\left(\mathbf{P}_{\cdot
g}\right)$,
$\boldsymbol{\Lambda}_{g}^{(5)}=\left(\boldsymbol{\Lambda}_{g,i}^{(5)}\right)_{n}=\big{(}F_{2}\left(\boldsymbol{\Theta}_{g},\mathbf{M}_{g,i\cdot},\mathbf{S}_{g,i\cdot},\boldsymbol{\mu}_{g}\right)\big{)}_{n}$
and
$K_{g}(\mathbf{Y})=\sum_{i,j}P_{ig}\left[-\log\left(Y_{ij}!\right)+Y_{ij}\log
l_{i}\right]$.

In real applications, we may have prior knowledge that some gene pairs cannot
have direct interactions. In this case, we can directly set the corresponding
edges as zero. Denote $E_{p}$ as the set of the edges that are priorly known
to be zero. Generally, we consider the following optimization problem

$\displaystyle\underset{\boldsymbol{\eta},\boldsymbol{\theta}}{\min}\left\\{-\ell_{\text{E}}\left(\boldsymbol{\theta},\boldsymbol{\eta}\right)+\lambda_{n}\sum_{g=1}^{G}\left\|{\boldsymbol{\Theta}_{g}}\right\|_{{1,\text{
off }}}:~{}{\boldsymbol{\Theta}}_{g}\succ 0,\
{\Theta_{g,lm}}=0~{}\mbox{for}~{}(l,m)\in
E_{p},\boldsymbol{\eta}\in\mathscr{H}\right\\}.$ (4)

###  3.2 The optimization process

We develop a block-wise descent algorithm called VMPLN to optimize (4). Let

$\displaystyle
L_{1}\left(\mathbf{M}_{g,i\cdot},\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g}\right)$
$\displaystyle=2^{-1}{\left(\mathbf{M}_{g,i\cdot}-\boldsymbol{\mu}_{g}\right)}^{T}\boldsymbol{\Theta}_{g}\left(\mathbf{M}_{g,i\cdot}-\boldsymbol{\mu}_{g}\right),$
$\displaystyle L_{2}\left(M_{g,ij},l_{i},S_{g,ij}\right)$
$\displaystyle=-\Lambda_{g,ij}^{(1)}+F_{1}\left(l_{i},M_{g,ij},S_{g,ij}^{(k+1)}\right),$
$\displaystyle\boldsymbol{\Sigma}\left(\mathbf{P}_{\cdot
g},\boldsymbol{\mu}_{g},\mathbf{M}_{g},\mathbf{S}_{g}\right)$
$\displaystyle=\left\\{\sum_{i=1}^{n}P_{ig}\left[(\mathbf{M}_{g,i\cdot}-\boldsymbol{\mu}_{g}){(\mathbf{M}_{g,i\cdot}-\boldsymbol{\mu}_{g})}^{T}+\mathbf{D}\left(\mathbf{S}_{g,i\cdot}\right)\right]\right\\}\bigg{/}\sum_{i=1}^{n}P_{ig}.$

Let
$\mathscr{N}=\left\\{1,2,\cdots,n\right\\},\mathscr{G}=\left\\{1,2,\cdots,G\right\\},\mathscr{P}=\left\\{1,2,\cdots,p\right\\}$.
Our proposed VMPLN algorithm is summarized in Algorithm 1. Given initial
values, we iteratively update $\mathbf{P}$, $\boldsymbol{\pi}$,
$\left\\{\mathbf{M}_{g}\right\\}_{g=1}^{G}$,
$\left\\{\mathbf{S}_{g}\right\\}_{g=1}^{G}$, $\boldsymbol{\mu}$ and
$\left\\{\boldsymbol{\Theta}_{g}\right\\}_{g=1}^{G}$ and terminate the
iteration if the ELBO and the network estimations only have very small changes
between two successive update steps. The step for updating $\mathbf{M}_{g}$
($g=1,\cdots,G$) is presented in the next subsection 3.3. The other steps are
straightforward and presented below.

Input: Single cell RNA-seq count data $\mathbf{Y}$, the pre-estimated library
size $\boldsymbol{l}$, the number of cell types $G$, the tunning parameter
$\lambda_{n}$, the optional prior set of zero edges $E_{p}$, the maximum
iteration number $K>0$ and the convergence thresholds
$\epsilon_{L},\epsilon_{s}$.

Output:
$\hat{\boldsymbol{\pi}},\hat{\boldsymbol{\mu}},\\{\hat{\boldsymbol{\Theta}}_{g}\\}_{g=1}^{G},\\{\hat{\mathbf{M}}_{g}\\}_{g=1}^{G},\\{\hat{\mathbf{S}}_{g}\\}_{g=1}^{G},\hat{\mathbf{P}}$.

Initialization step. Let $\delta_{L}=10^{6},\delta_{s}=p(p+1)/2,k=0$ and
initialize $\boldsymbol{\theta},\boldsymbol{\eta}$ as
$\boldsymbol{\theta}^{(0)},\boldsymbol{\eta}^{(0)}$ by equation (7).

while  _ $\left(\delta_{L}>\epsilon_{L}\ \textbf{ or }\
\delta_{s}>\epsilon_{s}\right)\ \textbf{ and }\ k\leq K$ _ do

$\mathbf{P}$-step. For each $\left(i,g\right)\in\mathscr{N}\times\mathscr{G}$,
compute
$U_{ig}^{(k)}=F_{2}\left(\hat{\boldsymbol{\Theta}}_{g}^{(k)},\hat{\mathbf{M}}_{g,i\cdot}^{(k)},\hat{\mathbf{S}}_{g,i\cdot}^{(k)},\hat{\boldsymbol{\mu}}_{g}^{(k)}\right)$
and update
${\hat{P}_{ig}}^{(k+1)}=\hat{\pi}_{g}^{(k)}\exp\left({U}_{ig}^{(k)}\right)\bigg{/}\sum_{l=1}^{G}\left[\hat{\pi}_{l}^{(k)}\exp\left({U}_{il}^{(k)}\right)\right],$

$\boldsymbol{\pi}$-step. For each $g\in\mathscr{G}$, update
${\hat{\pi}_{g}}^{(k+1)}=n^{-1}\sum_{i=1}^{n}{\hat{P}_{ig}}^{(k+1)}$

$\mathbf{M}$-step. For each $\left(i,g\right)\in\mathscr{N}\times\mathscr{G}$,
update
${\hat{\mathbf{M}}_{g,i\cdot}}^{(k+1)}=\underset{\mathbf{M}_{g,i\cdot}}{\operatorname*{\arg\min}}\left\\{L_{1}\left(\mathbf{M}_{g,i\cdot},\hat{\boldsymbol{\Theta}}_{g}^{(k)},\hat{\boldsymbol{\mu}}_{g}^{(k)}\right)+\sum_{j=1}^{p}L_{2}\left(M_{g,ij},l_{i},\hat{S}_{g,ij}^{(k)}\right)\right\\}.$
(5)

$\mathbf{S}$-step. For each
$\left(i,g,j\right)\in\mathscr{N}\times\mathscr{G}\times\mathscr{P}$, update
${\hat{\mathbf{S}}_{g,ij}}^{(k+1)}=\underset{S_{g,ij}>0}{\operatorname*{\arg\min}}\left\\{F_{1}\left(l_{i},\hat{M}_{g,ij}^{(k+1)},S_{g,ij}\right)+2^{-1}\hat{\boldsymbol{\Theta}}_{g,jj}^{(k)}S_{g,ij}-2^{-1}\log
S_{g,ij}\right\\}.$

$\boldsymbol{\mu}$-step. For each
$\left(g,j\right)\in\mathscr{G}\times\mathscr{P}$, update
$\hat{\mu}_{gj}^{(k+1)}=\sum_{i=1}^{n}\left(\hat{P}_{ig}^{(k+1)}\hat{M}_{g,ij}^{(k+1)}\right)\bigg{/}\sum_{i=1}^{n}\hat{P}_{ig}^{(k+1)}.$

$\boldsymbol{\Theta}_{g}$-step. For each $g\in\mathscr{G}$, compute
$\hat{\boldsymbol{\Sigma}}_{g}^{(k+1)}=\boldsymbol{\Sigma}\left(\hat{\mathbf{P}}_{.g}^{(k+1)},\hat{\boldsymbol{\mu}}_{g}^{(k+1)},\hat{\mathbf{M}}_{g}^{(k+1)},\hat{\mathbf{S}}_{g}^{(k+1)}\right)$
and update
$\displaystyle\hat{\boldsymbol{\Theta}}_{g}^{(k+1)}=\underset{\boldsymbol{\Theta}_{g}}{\operatorname*{\arg\min}}\bigg{\\{}$
$\displaystyle-2^{-1}\log\det\boldsymbol{\Theta}_{g}+2^{-1}\operatorname{tr}\big{(}\boldsymbol{\Theta}_{g}\hat{\boldsymbol{\Sigma}}_{g}^{(k+1)}\big{)}+\left(\lambda_{n}/\mathbf{1}_{n}^{T}\hat{\mathbf{P}}_{\cdot
g}^{(k+1)}\right)\left\|{\boldsymbol{\Theta}_{g}}\right\|_{{1,\text{ off }}}:$
$\displaystyle{\boldsymbol{\Theta}}_{g}\succ 0,\
{\Theta_{g,lm}}=0~{}\mbox{for}~{}(l,m)\in E_{p}\bigg{\\}}.$

Evaluation step. Update  $\displaystyle\delta_{s}$
$\displaystyle=\underset{g\in\mathscr{G}}{\max}\left\\{\sum_{l<m}\left|\text{sign}\left(\hat{\boldsymbol{\Theta}}_{g,lm}^{(k+1)}\right)-\text{sign}\left(\hat{\boldsymbol{\Theta}}_{g,lm}^{(k)}\right)\right|{\tbinom{p}{2}}^{-1}\right\\},$
$\displaystyle\delta_{L}$
$\displaystyle=\delta\left(\ell_{\text{E}}\left(\hat{\boldsymbol{\eta}}^{(k+1)},\hat{\boldsymbol{\theta}}^{(k+1)}\right),\ell_{\text{E}}\left(\hat{\boldsymbol{\eta}}^{(k)},\hat{\boldsymbol{\theta}}^{(k)}\right)\right),\
\text{with}\ \delta\left(a,b\right)=\left|a-b\right|/b,$ $\displaystyle k$
$\displaystyle=k+1.$

Algorithm 1 Framework of VMPLN.

The parameters $\mathbf{P}$, $\boldsymbol{\pi}$ and $\boldsymbol{\mu}$ all
have explicit updating formulas and can be efficiently calculated. The
updating form of $\hat{\boldsymbol{\mu}}_{g}$ is different from that of the
centered variational approach used by VPLN (Chiquet et al., 2019), which
involves the logarithm of the sum of $Y_{ij}$ and may be numerically unstable
when the data is of high dropout.

For the parameter $\mathbf{S}$, given all other parameters, the ELBO loss
function (4) can be decomposed into the sum of $npG$ functions, each of which
only involves one $S_{g,ij}$. Thus, in the $\mathbf{S}$-step of Algorithm 1,
the optimization problem can be decomposed into $npG$ one-dimensional convex
optimization problems. Each of these can be efficiently solved by the Newton-
Raphson algorithm.

For the network parameters $\boldsymbol{\Theta}_{g}$ ($g=1,\cdots,G$), given
all other parameters, the corresponding sub-optimization problem is equivalent
to $G$ independent Glasso problems (Meinshausen and Bühlmann, 2006; Friedman
et al., 2008), and $\boldsymbol{\Theta}_{g}$ can be updated by solving the
corresponding Glasso problem with its covariance matrix as
$\hat{\boldsymbol{\Sigma}}_{g}^{(k+1)}$ and its penality parameter as
$\lambda_{n}\big{/}\mathbf{1}_{n}^{T}\hat{\mathbf{P}}_{.g}^{(k+1)}$.

###  3.3 The optimization of the $\mathbf{M}$-step in Algorithm 1

The sub-optimization problem corresponding to $\mathbf{M}$ is

$\underset{\mathbf{M}}{\operatorname*{\arg\min}}\left\\{\sum_{g\in\mathscr{G},i\in\mathscr{N}}\mathbf{P}_{ig}\left[L_{1}\left(\mathbf{M}_{g,i\cdot},\hat{\boldsymbol{\Theta}}_{g}^{(k)},\hat{\boldsymbol{\mu}}_{g}^{(k)}\right)+\sum_{j=1}^{p}L_{2}\left(M_{g,ij},l_{i},\hat{S}_{g,ij}^{(k)}\right)\right]\right\\},$

which is equivalent to $nG$ independent optimization problems (5). The
p-dimensional optimization problem (5) can be solved by the $p$-dimensional
Newton-Raphson algorithm or coordinate descent algorithm. However, the
$p$-dimensional Newton-Raphson algorithm involves inverting $p\times p$
Heissen matrices. In each of the $nG$ optimization problems (5) and at each
step of the $p$-dimensional Newton-Raphson algorithm, the Heissen matrices are
unique and different from each other. In total, the $p$-dimensional Newton-
Raphson algorithm would involve $O(nGt)$ $p\times p$ matrix inversions, where
$t$ is the number of iterations in the Newton-Raphson algorithm, and is
computationally very expensive. The coordinate descent algorithm is also slow
because of its inferior convergence rate and the large number of optimization
problems.

Here, we instead develop a more efficient algorithm based on the alternating
direction method of multipliers (ADMM) algorithm (Boyd et al., 2011) for the
optimization problem (5). Specifically, we introduce an auxiliary matrix
$\mathbf{N}_{g}$ for $\mathbf{M}_{g}$, and denote
$\ell_{\mathbf{M}}\left(\mathbf{M}_{g,i\cdot},\mathbf{N}_{g,i\cdot}\right)=L_{1}\left(\mathbf{N}_{g,i\cdot},\hat{\boldsymbol{\Theta}}_{g}^{(k)},\hat{\boldsymbol{\mu}}_{g}^{(k)}\right)$
$+\sum_{j=1}^{p}L_{2}\left(M_{g,ij},l_{i},\hat{S}_{g,ij}^{(k)}\right)$.
Solving (5) is equivalent to solving the following problem

$\underset{\mathbf{M}_{g,i\cdot}=\mathbf{N}_{g,i\cdot}}{\operatorname*{\arg\min}}\left\\{\ell_{\mathbf{M}}\left(\mathbf{M}_{g,i\cdot},\mathbf{N}_{g,i\cdot}\right)\right\\}.$
(6)

The augmented Lagrangian of (6) is

$\ell_{\mathbf{M}}\left(\mathbf{M}_{g,i\cdot},\mathbf{N}_{g,i\cdot}\right)+\sum_{j=1}^{p}{\alpha}_{j}\left(M_{g,ij}-N_{g,ij}\right)+\rho/2\sum_{j=1}^{p}{\left(M_{g,ij}-N_{g,ij}\right)}^{2},$

where $\boldsymbol{\alpha}=\left(\alpha_{1},\cdots,\alpha_{p}\right)$ is the
Lagrangian multiplier, and $\rho$ is the step size. The corresponding ADMM
algorithm is detailed in Algorithm 2. Given initial values, we iteratively
update $\mathbf{M}_{g,i\cdot}$, $\mathbf{N}_{g,i\cdot}$ and
$\boldsymbol{\alpha}$, and terminate the iteration if the ELBO only has a very
small change between two successive update steps.

Input: Single cell RNA-seq count data $\mathbf{Y}$, the pre-estimated library
size $\boldsymbol{l}$, the current estimation of parameters
$\hat{\boldsymbol{\theta}},\hat{\boldsymbol{\eta}}$, the maximum iteration
number $T>0$ and the convergence threshold $\epsilon_{M}$.

Output: $\hat{\mathbf{M}}_{g,i\cdot}$.

Initialization step. Let $t=0,\delta_{M}=10^{6}$, and initialize
$\hat{\mathbf{M}}_{g,i\cdot}^{\left(0\right)}=\hat{\mathbf{N}}_{g,i\cdot}^{\left(0\right)}=\hat{\mathbf{M}}_{g,i\cdot}$,
$\hat{\boldsymbol{\alpha}}^{(0)}=\mathbf{0}$.

while  _$t\leq T$ and $\delta_{M}>\epsilon_{M}$_ do

Step 1. For each $j\in\mathscr{P}$, update
$\hat{M}_{g,ij}^{(t+1)}=\underset{M_{g,ij}}{\operatorname*{\arg\min}}\bigg{\\{}L_{2}\left(M_{g,ij},l_{i},\hat{S}_{g,ij}\right)+{\hat{\alpha}}_{j}^{(t)}(M_{g,ij}-\hat{N}_{g,ij}^{(t)})+\frac{\rho}{2}{(M_{g,ij}-\hat{N}_{g,ij}^{(t)})}^{2}\bigg{\\}}.$

Step 2. Update
$\hat{\mathbf{N}}_{g,i\cdot}^{(t+1)}={\left(\rho\mathbf{I}+\hat{\boldsymbol{\Theta}}_{g}\right)}^{-1}\left(\rho\hat{\mathbf{M}}_{g,i\cdot}^{(t+1)}-\hat{\boldsymbol{\alpha}}^{(t)}+\hat{\boldsymbol{\Theta}}_{g}\hat{\boldsymbol{\mu}}_{g}\right).$

Step 3. Update
${\hat{\boldsymbol{\alpha}}}^{(t+1)}={\hat{\boldsymbol{\alpha}}}^{(t)}+\rho\left(\hat{\mathbf{M}}_{g,i\cdot}^{(t+1)}-\hat{\mathbf{N}}_{g,i\cdot}^{(t+1)}\right).$

Evaluation step. Update  $\displaystyle\delta_{M}$
$\displaystyle=\delta\left(\ell_{\mathbf{M}}\left(\hat{\mathbf{M}}_{g,i\cdot}^{(t+1)},\hat{\mathbf{N}}_{g,i\cdot}^{(t+1)}\right),\ell_{\mathbf{M}}\left(\hat{\mathbf{M}}_{g,i\cdot}^{(t)},\hat{\mathbf{N}}_{g,i\cdot}^{(t)}\right)\right),\
\text{with}\ \delta\left(a,b\right)=\left|a-b\right|/b,$ $\displaystyle t$
$\displaystyle=t+1.$

Algorithm 2 ADMM algorithm for updating $\mathbf{M}_{g,i\cdot}$.

Given all other variables, the optimization problem for
$\mathbf{M}_{g,i\cdot}$ can be decomposed into $p$ independent one-dimensional
optimization problems, which can be easily optimized by the one-dimensional
Newton-Raphson algorithm (Step 1 in Algorithm 2). The optimization problem for
$\mathbf{N}_{g,i\cdot}$ has an explicit solution and involves inverting the
matrix $\rho\mathbf{I}+\hat{\boldsymbol{\Theta}}_{g}^{(k)}$. Note that in this
ADMM algorithm, we only need to calculate the matrix inversion
$(\rho\mathbf{I}+\hat{\boldsymbol{\Theta}}_{g}^{(k)})^{-1}$ once. Thus, this
ADMM algorithm is more efficient than the $p$-dimensional Newton-Raphson
algorithm. Moreover, the $\mathbf{M}$-step of Algorithm 1 involves $nG$
optimization problems (7), but we only need to calculate $G$ matrix inversions
$(\rho\mathbf{I}+\hat{\boldsymbol{\Theta}}_{g}^{(k)})^{-1}$ in the $nG$ runs
of the ADMM Algorithm 2. In total, using this ADMM algorithm to solve (7)
makes Algorithm 1 much more efficient than using the Newton-Raphson algorithm.

###  3.4 Initialization, tunning parameters selection

To initialize the parameters in the initialization step of Algorithm 1, we
perform dimension reduction using principal component analysis (PCA) on the
normalized data
$\tilde{\mathbf{Y}}\triangleq\log(\mathbf{Y}+1)-(\log\boldsymbol{l})\mathbf{1}_{p}^{T}$.
Then, we use K-means (Hartigan and Wong, 1979) to cluster single cells in the
low dimensional space. Let $\tilde{Z}_{i}\in\\{1,2,\cdots,G\\}$ be the
clustering label of the $i$th cell. The parameters are initialized as

$\displaystyle\hat{P}_{ig}^{(0)}$
$\displaystyle=I(\tilde{Z}_{i}=g),~{}\hat{M}_{g,ij}^{(0)}=\tilde{Y}_{ij},~{}\hat{S}_{g,ij}^{(0)}=10^{-5},$
(7) $\displaystyle\hat{\pi}_{g}^{(0)}$
$\displaystyle=n^{-1}\sum_{i=1}^{n}\hat{P}_{ig}^{(0)},~{}\hat{\mu}_{gj}^{(0)}=\bigg{(}\sum_{i=1}^{n}\hat{P}_{ig}^{(0)}\bigg{)}^{-1}\sum_{i=1}^{n}\hat{P}_{ig}^{(0)}\hat{M}_{g,ij}^{(0)},$
$\displaystyle\text{with}\ g=1,\cdots,G;\ i=1,\cdots,n;\ j=1,\cdots,p.$

The precision matrix ${\hat{\boldsymbol{\Theta}}_{g}}^{(0)}$ of the $g$th cell
type is initialized by the Glasso algorithm with its covariance matrix as
${\hat{\boldsymbol{\Sigma}}}_{g}^{(0)}=\boldsymbol{\Sigma}\left(\hat{\mathbf{P}}_{.g}^{(0)},\hat{\boldsymbol{\mu}}_{g}^{(0)},\hat{\mathbf{M}}_{g}^{(0)},\hat{\mathbf{S}}_{g}^{(0)}\right)$
and its penalty parameter as $10^{-6}$.

The tunning parameter $\lambda_{n}>0$ is selected for each cell type
independently by the integrated complete likelihood (ICL) criterion (Biernacki
et al., 2000), which selects the tunning parameter $\lambda_{n}>0$ for
$\boldsymbol{\Theta}_{g}$ by minimizing

$-2\
{l^{(g)}_{\text{E}}}\left(\hat{\boldsymbol{\eta}},\hat{\boldsymbol{\theta}}\right)+\log\left(\mathbf{1}_{n}^{T}\hat{\mathbf{P}}_{\cdot
g}\right)s(\hat{\boldsymbol{\Theta}}_{g}),$ (8)

where $s(\hat{\boldsymbol{\Theta}}_{g})$ denote the number of non-zero
elements in $\hat{\boldsymbol{\Theta}}_{g}$. As another choice, we can select
the tuning parameter $\lambda_{n}$ such that the estimated network has a
desired density.

##  4 Simulation

In this section, we perform simulation to evaluate the performance of VMPLN
and compare with state-of-the-art algorithms, including VPLN (Chiquet et al.,
2019), Glasso (Meinshausen and Bühlmann, 2006; Friedman et al., 2008), LPGM
(Allen and Liu, 2013), PPCOR (Kim, 2015), GENIE3 (Huynh-Thu et al., 2010) and
PIDC (Chan et al., 2017). PPCOR infers the regulatory network based on the
partial correlation coefficients. LPGM is an extension of the PGM. GENIE3, a
random-forest-based regulatory network method, is the best performing
algorithm in the DREAM4 GRN inference challenge (Greenfield et al., 2010).
PIDC is a recently developed algorithm for single cell data based on the
mutual information.

###  4.1 Simulation setups

We first generate simulation data based on the MPLN model. The number of
components is set as $G=3$ and the proportion parameter $\boldsymbol{\pi}$ is
set as $(1/3,1/3,1/3)$. The number of observations is $n=3000$. We consider 48
different simulation scenarios, which are 3 cell-population-mixing levels
(low, middle and high) $\times$ 2 dropout levels (low and high) $\times$ 2
dimension setups ($p=100$ and $300$) $\times$ 4 graph structures. In each
scenario, we generate 50 datasets. The four graph structures are as follows.

1\. Random Graph: Pairs of nodes are connected with probability 0.1. The
nonzero edges are randomly set as 0.3 or -0.3.

2\. Hub Graph: 20 $\%$ nodes are set as hub nodes. A hub node is connected
with another node with probability 0.1. Non-hub nodes are not connected with
each other. The nonzero edges are randomly set as 0.3 or -0.3.

3\. Blocked random Graph: The nodes are divided into 5 blocks of equal sizes.
Pairs of nodes within the same block are connected with probability 0.1. Nodes
in different blocks are not connected. The nonzero edges are randomly set as
0.3 or -0.3.

4\. Scale-free Graph: The Barabasi-Albert model (Barabási and Albert, 1999) is
used to generate a scale-free graph with power 1. The nonzero edges are
randomly set as 0.3 or -0.3.

Details of the data generation process are shown in supplementary material
(Section S1.2.1). To test the performance of the proposed method under the
misspecified model setting, we also generate simulation data using a mixture
multinomial log-normal distribution. In this model, the conditional Poisson
layer of the MPLN is replaced with a conditional multinomial distribution,
$\mathbf{Y}_{i}|\mathbf{X}_{i}\sim\mbox{Multinomial}(\lfloor\sum_{j=1}^{p}\lambda_{ij}\rfloor,\boldsymbol{\lambda}_{i}\big{/}\sum_{j=1}^{p}\lambda_{ij})$,
where $\lambda_{ij}=l_{i}\ \exp(X_{ij})$,
$\boldsymbol{\lambda}_{i}=\left(\lambda_{i1},\cdots,\lambda_{ip}\right)$ and
$\lfloor\sum_{j=1}^{p}\lambda_{ij}\rfloor$ represents the maximum integer that
does not exceed $\sum_{j=1}^{p}\lambda_{ij}$. The simulation data are
similarly generated.

###  4.2 Performance Comparison

For the proposed VMPLN algorithm, we use clustering results given by the
K-means algorithm as the initial value and infer the regulatory networks
jointly for all cell types. For the other algorithms, we use K-means to assign
cell types of the single cells and infer regulatory network for each cell type
separately. We compare VMPLN with other algorithms in terms of their
accuracies of network inference.

Following Pratapa et al. (2020), we define confidence scores of the predicted
edges for each method based on its reported statistics (supplementary
material, Section S1.2.2), and evaluate the algorithms by the area under the
precision-recall curve (AUPRC) ratios and the early precision (EP) ratios
based on the confidence scores. For VMPLN, VPLN and Glasso, suppose that
$\hat{\boldsymbol{\Theta}}$ is its estimation of a network, we define a
confidence score for the edge $(i,j)$ ($i\neq j$) as its absolute partial
correlation coefficient, i.e.
$\big{|}-(\hat{\Theta}_{ii}\hat{\Theta}_{jj})^{-\frac{1}{2}}\hat{\Theta}_{ij}\big{|}$.
For LPGM, we define a confidence score for each edge as its stability score.
For PPCOR, GENIE3 and PIDC, we define a confidence score for each edge as its
estimated connected weight. By varying the confidence score, we can obtain the
AUPRC for each algorithm on each simulation data. The AUPRC ratio is defined
as the ratio between the AUPRC of an algorithm and the AUPRC of the random
predictor. The EP is defined as the precision of the top $\min\\{K,s\\}$ edges
in the inferred network, where $K$ is the number of edges in the true network,
and $s$ is the number of edges in the inferred network. The EP ratio is
defined as the ratio between the EP of the algorithm and the EP of the random
predictor.

(a) AUPRC ratios

(b) EP ratios

Figure 1: The AUPRC ratios and EP ratios for the hub graph. The parameters are
set as their default values or are tuned by their default methods.

We first compare the estimated networks given by different algorithms with
their default parameters or default ways of selecting tuning parameters
(supplementary material, Section S1.2.2). Figure 1 shows the boxplots of AUPRC
ratios and EP ratios of different algorithms for the hub graph when the
parameters are selected using their default methods. The results for the
random graph, blocked random graph and scale-free graph are shown in
supplementary material (Figure S1-S3). Overall, VMPLN is the best performing
algorithm in terms of the AUPCR ratio. The advantage of VMPLN is more
pronounced when the cell types have a higher mixing level, suggesting that
compared with the two-step procedure, the joint analysis of network inference
and clustering can help to improve the network inference. In terms of the EP
ratio, VMPLN is also among the best performing algorithms. LPGM often has
larger EP ratios than VMPLN in simulation settings with low cell-type mixing
degrees. This is because LPGM uses stability to select the tuning parameter
and thus is very conservative. LPGM often only reports a few edges, which are
mostly true discoveries, leading to its high precision (as measured by the EP
ratio) but a low sensitivity (supplementary material, Figure S4). When the
cell-type mixing level is high, VMPLN also has higher EP ratio values than
other methods including LPGM.

(a) AUPRC ratios

(b) EP ratios

Figure 2: The AUPRC ratios and EP ratios for the hub graph. The confidence
score cutoffs or the tuning parameters are selected such that the network
density is 20%.

Different ways of selecting the tuning parameters can lead to network
estimations having very different densities. Dense network predictions usually
have a high sensitivity and a low specificity, but sparse network predictions
have a low sensitivity and a high specificity. Thus, dense and sparse network
predictions may not be directly comparable. To eliminate the influence of the
tuning parameter selection methods, we further compare the algorithms at the
same network density 20% (2 $\times$ the density of the true network)
(supplementary material, Section S1.2.2). Figure 2 shows the boxplots of AUPRC
ratios and EP ratios for the hub graph at the same network density. The
results for the random graph, blocked random graph and scale-free graph are
shown in supplementary material (Figure S5-S7). Similarly, VMPLN has the best
performance in most scenarios, especially in the cases with high cell-type
mixing levels. For example, in the simulation of the hub graph with $p=100$
and low-dropout, VMPLN has mean AUPRC ratios $8.63$ and $8.72$ in the high and
low cell-type mixing scenarios, respectively, about $21\%$ and $5\%$ larger
than the AUPRC ratios ($7.12$ and $8.26$) of VPLN under the same scenarios.
All algorithms tend to have decreased performances in higher dimensions or
with high dropout rates and VMPLN consistently has better performances in
these more difficult settings. Even in the low-mixing level case, where the
two-step methods should work well, VMPLN still performs better than VPLN
(supplementary material, Figure S8).

The simulation results for the data generated from the misspecified model (the
compositional model) are shown in supplementary material (Figure S9-S16).
Similarly, VMPLN also performs better in most simulation settings.

##  5 Real data analysis

###  5.1 Benchmarking on scRNA-seq data

In this section, we evaluate VMPLN and compare with other algorithms using two
real scRNA-seq datasets. LPGM is not included in this comparison because it
was unable to finish computation in a reasonable amount of time (7 days). One
dataset is the scRNA-seq of human peripheral blood mononuclear cells (PBMC)
profiled by Kang et al. (2018) (Kang data) and another dataset is the scRNA-
seq data of human PBMC cells profile by Zheng et al. (2017) (Zheng data). The
Kang data consists of two batches, the interferon Beta 1 (IFNB1)-stimulated
and control groups. The Zheng data also consists of two batches, which are
respectively sequenced by $3^{\prime}$ and $5^{\prime}$ scRNA-seq
technologies. We consider 10 cell types (7217 cells) from the Kang data and 6
cell types (5962 cells) from the Zheng data for network analysis. The other
cell types of the Kang and Zheng datasets have less than 150 cells and are not
considered for the network analysis. The number of cells in each cell type is
listed in supplementary material (Table S1).

In both datasets, we first use data from one of the two batches (construction
data) and public GRN databases (supplementary material, Table S2) to construct
a silver standard (supplementary material, Section S1.3.1). Then, we test
different algorithms using the another batches (testing data). For the testing
data, we select the top 1000 highly variable genes (HVGs) using Seurat (Stuart
et al., 2019) and cluster the cells using these HVGs. We select the top 500
HVGs as the gene sets of interest for GRN inference. The top 500 HVGs contain
39 and 19 TFs for the Kang and Zheng data, respectively. For algorithms other
than VMPLN, the cell types are first identified using clustering results of
the top 1000 HVGs and gene regulatory relationships between genes in the gene
sets of interest are inferred for each cell type. To make the comparison fair,
for VMPLN, clustering and network inference are simultaneously performed using
top 1000 HVGs and network inference focuses on genes in the gene sets of
interest.

Figure 3: The performance of the netwok inference algorithms in the Kang data.
The gene set of interest for GRN inference consists of 39 TFs and 461 non-TFs.
The colors represent the scaled values of these metrics (scaled to between 0
and 1 within each cell type) and the actual values are marked in the boxes.
Black color in the boxes: random predictor performs better.

We compare the algorithms at the same network density ($5\%$) in terms of the
AUPRC ratio, the EP ratio and the early stability. The AUPRC ratios and EP
ratios are defined in the simulation section and are computed by comparing the
estimated network with the silver standard. Since the regulatory relationships
in the silver standard sets all involve TFs, we only consider edges involving
TFs for performance evaluation. The early stability is to measure the
stability of each algorithm by perturbing the input data. To calculate the
stability, we down-sample $90\%$ of the original data, estimate the network
using the same estimating procedure based on the down-sampled data, and
calculate the pairwise Jaccard index between the networks estimated based on
the down-sampled data. This process is repeated 100 times and the early
stability is defined as the median of the Jaccard indexes.

Figure 3 and supplementary material (Figure S17) show the heatmap of the AUPRC
ratio, the EP ratio and the early stability of the 6 methods for these two
benchmarking datasets. The algorithms are ordered by the overall AUPRC ratio
across all cell types. VMPLN has the highest AUPRC ratios and EP ratios in
most cases and is the only method that consistently performs better than the
random predictor. The stability of VMPLN is also reasonable and roughly
similar to GENIE3. Note that GENIE3 is a tree-ensemble based method that uses
data perturbation for network estimation, and thus should have good stability.

###  5.2 Application to scRNA-seq data from COVID-19 patients

In this section, we consider the scRNA-seq data of bronchoalveolar lavage
fluid macrophage cells from coronavirus disease 2019 (COVID-19) patients (Liao
et al., 2020). The data consists of 29,980 single cells from 8 patients,
including 2 patients with moderate COVID-19 infection and 6 patients with
severe infection (supplementary material, Table S3). Liao et al. (2020)
clustered the macrophages to four clusters including two classic M1-like
macrophage groups (Group1 and Group2), the alternative M2-like macrophages
(Group3) and the alveolar macrophages (Group4).

Figure 4: The GRN analysis of the COVID-19 data. (a) The inferred GRNs of
Group4 macrophages in severe and moderate patients. The size of the node
represents the weighted node degree. (b) GO enrichment analysis of TFs’ target
genes. TFs are selected as those with a large weighted degree differences
($>0.2$) between severe and moderate patients. Activation and Inhibition mean
that the regulatory relationships are positive and negative, respectively.
Left panel: p-values of GO terms. Right panel: differences of number of genes
in the GO terms between severe and moderate patients. (c) The genes regulated
by NR4A1. (d) Differential gene expression analysis between cells with and
without SARS-CoV-2 infection. MoAM2: monocyte derived aleovar macrophage 2.
TRAM: tissue resident aleovar macrophage 2. One-sided Wilcoxon rank-sum test
is used.

We perform the GRN analysis of the top 1000 HVG genes for the four macrophage
groups using VMPLN (Figure 4 a; supplementary material, Section S1.3.2 and
Figure S18-S20) and compare networks of the patients with severe and moderate
infection. We focus on the alveolar macrophages, since unlike other macrophage
groups, the proportion of the alveolar macrophages tends to be smaller in
patients with severe infection (Liao et al., 2020). A number of TFs exhibit a
large weighted degree difference between the GRNs in the moderate and severe
patients (weighted degree difference greater than 0.2; supplementary material,
Figure S21). The weighted degree of a node is defined as the summation of the
absolute partial correlation between the node and all other nodes that are
connected with it. The node degrees of TCF4, BHLHE40, JUNB, NR1H3, and JUN in
severe patients are substantially larger than in moderate patients, while the
node degrees of ELF3, NR4A1, EGR2 and ATF3 are substantially larger in
moderate patients. Gene oncology (GO) enrichment analysis of their target
genes (Figure 4 b) shows that, as expected, many of target genes are involved
in immune responses such as leukocyte migration and regulation of T cell
activation. Interestingly, we observe that a number of GO terms are only
enriched in severe patients such as response to unfolded protein (UPR) and
response to endoplasmic reticulum (ER) stress. Similar analyses of the other
macrophage groups also show that UPR is more enriched in severe patients than
in moderate patients (supplementary material, Figure S22-S24). The enrichment
of the UPR process is mainly due to the activation of target genes of NR4A1
(Figure 4 b; supplementary material, Figure S22-S24). The UPR and ER stress
processes are frequently activated in cells infected by viruses (Janssens et
al., 2014) including coronavirus (Chan et al., 2006), indicating that alveolar
macrophages might be infected by SARS-CoV-2. In fact, a recent study showed
that macrophages can be infected by SARS-CoV-2 and the infected macrophages
activate T-cells to promote alveolitis in patients with severe COVID-19 (Grant
et al., 2021). These data imply that the UPR and ER stress related genes might
be activated in SARS-CoV-2-infected macrophages through modulation of TFs such
as NR4A1. Macrophage single cell data in Grant et al. (2021) has SARS-CoV-2
infection information for each single cell and thus allows expression
comparison between cells with or without SARS-CoV-2 infection. We take the UPR
and ER stress related genes that are regulated by NR4A1 only in severe
patients (Figure 4 c) and compare their expressions in cells with and without
SARS-CoV-2 infection. We find that most of these genes are indeed
significantly differentially expressed (Figure 4 d). In light of these
findings, we reason that NR4A1 might play an important role in regulating
cellular responses to SARS-CoV-2 infection. A number of NR4A1’s target genes,
including IRF1 and HSP90, have recently been discovered to be potential
therapeutic targets for COVID-19 (Echavarría-Consuegra et al., 2021; Shaban et
al., 2021). NR4A1 and its target genes that we identified here may also serve
as potential therapeutic targets.

##  6 Discussion

In this paper, we develop a regulatory network inference method called VMPLN
for scRNA-seq data. Instead of using the two-step procedure for network
inference, VMPLN performs clustering and network inference simultaneously, and
thus are especially suitable for scRNA-seq with mixed cell types. Most of the
scRNA-seq data contain multiple cell types. We expect that VMPLN will have
many applications in single cell studies.

A potential limitation of VMPLN is that it assumes that the regulatory
relationships are linear. If the regulatory relationships are far away from
being linear, VMPLN will not perform well. Methods like GENIE3 can allow the
regulatory relationships to be nonlinear, but they require pre-clustering
before network inference. One important research direction is to develop
network methods that can account for nonlinear regulatory relationship as well
as mixed cell populations. In addition, VMPLN is developed for scRNA-seq data.
Other data and information can only be incorporated by setting the prior
edges. Currently, single cell multi-omics technologies have been developed
(Chappell et al., 2018). Developing network inference methods that can
integrate multi-omics data can help to improve sensitivity and reduce false
discoveries.

##  7 Acknowledgments

This work was supported by the National Key Basic Research Project of China
(2020YFE0204000), the National Natural Science Foundation of China (11971039),
and Sino-Russian Mathematics Center.

@columngridonene

SUPPLEMENTARY MATERIAL

##  S1 Supplementary Text

###  S1.1 Technical proofs

####  S1.1.1 Notation

For notational simplicity, we do not distinguish different lower (upper)
bounds in Condition (C1-C3), and always use $m$ to denote the lower bound and
$M$ to denote the upper bound. We always assume $M>0$. When $m$ is the lower
bound for eigenvalues of precision matrices, we always assume $m>0$. We define
two vectorization operators, ${\rm vech}$ and ${\rm vech_{2}}$. For a
symmetric matrix $\mathbf{A}\in\mathbb{R}^{p\times p}=[a_{ij}]$, ${\rm
vech}(\mathbf{A})$ is defined as

${\rm
vech}(\mathbf{A})=(a_{11},a_{12},a_{13},\ldots,a_{1p},a_{22},a_{23},\ldots,a_{2p},\ldots,a_{(p-1)(p-1)},a_{(p-1)p},a_{pp})^{\rm
T},$

and ${\rm vech_{2}}(\mathbf{A})$ is

${\rm
vech_{2}}(\mathbf{A})=(a_{11},2a_{12},2a_{13},\ldots,2a_{1p},a_{22},2a_{23},\ldots,2a_{2p},\ldots,a_{(p-1)(p-1)},2a_{(p-1)p},a_{pp})^{\rm
T}.$

Note that ${\rm vech}$ and ${\rm vech_{2}}$ only differ at off-diagonal
elements. Define ${\boldsymbol{\nu}}=(\boldsymbol{\nu}_{1}^{\rm
T},\ldots,\boldsymbol{\nu}_{g}^{\rm T},\ldots,\boldsymbol{\nu}_{G}^{\rm
T})^{\rm T},$ where $\boldsymbol{\nu}_{g}={\rm
vech}(\boldsymbol{\Theta}_{g})$. We assume the true parameter
$\boldsymbol{\nu}^{*}$ is an interior point of $\mathcal{D}$. For the PLN, we
can write the log-likelihood as follows,

$\log\,p(\mathbf{y},\boldsymbol{l};\boldsymbol{\Theta},\boldsymbol{\mu})=\sum_{i=1}^{n}\log\int
f(\mathbf{x};\boldsymbol{\Theta},\boldsymbol{\mu})p(\mathbf{y}_{i}|\mathbf{x},l_{i})d\mathbf{x}+\sum_{i=1}^{n}\log\left(p(l_{i})\right).$

where $f(\mathbf{x};\boldsymbol{\Theta},\boldsymbol{\mu})$ is the probability
density function of ${\rm N}(\boldsymbol{\Theta},\boldsymbol{\mu})$. Note that
$\sum_{i=1}^{n}\log\left(p(l_{i})\right)$ does not depend on
$\boldsymbol{\Theta}$. We only need to consider the conditional distribution
$p(\mathbf{y}_{i}|l_{i};\boldsymbol{\Theta},\boldsymbol{\mu}_{g})$. Since
$p(l)$ is independent of the unknown parameters and has a bounded support, it
can be seen from the following section that the library size $l$ does not have
essential influence on the proof. For brevity, we assume that the library size
is the constant $1$ in the proof.

In the following sections, we always use
$p(\mathbf{y};\boldsymbol{\Theta},\boldsymbol{\mu})$ and
$p\left(\mathbf{y};\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)$
to represent the density of the PLN distribution ${\rm
PLN}\left(\boldsymbol{\Theta},\boldsymbol{\mu}\right)$ and the density of the
MPLN distribution MPLN
$\left(\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)$,
respectively. Given a single sample $i$, we write the log-likelihood function
of the PLN at $\mathbf{y}_{i}$ as

$\displaystyle\ell(\boldsymbol{\Theta},\mathbf{y}_{i})$
$\displaystyle=\log\left(p(\mathbf{y}_{i};\boldsymbol{\Theta},\boldsymbol{\mu})\right)=\frac{1}{2}\log\det(\boldsymbol{\Theta})$
$\displaystyle+\log\int\det(\boldsymbol{\Theta})^{\frac{1}{2}}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right){h(\mathbf{y}_{i},\mathbf{x})}d\mathbf{x}+C(\mathbf{y}_{i}),$

where

$h(\mathbf{y}_{i},\mathbf{x})=\prod_{j=1}^{p}\exp(x_{j}y_{ij})\exp\left(-\exp(x_{j})\right)$
(S1)

and $C(\mathbf{y}_{i})=\sum_{j=1}^{p}\log y_{ij}!-2^{-1}p\log(2\pi)$. In the
following sections, we always write
$h(\mathbf{y},\mathbf{x})=\prod_{j=1}^{p}\exp(x_{j}y_{j})\exp(-\exp(x_{j}))$.
Also, we define
$\ell_{n}(\boldsymbol{\Theta})=\sum_{i=1}^{n}\ell(\boldsymbol{\Theta},\mathbf{y}_{i})$
as the log-likelihood in the PLN. For the MPLN, its log-likelihood function at
$\mathbf{y}_{i}$ is

$\displaystyle\ell(\boldsymbol{\nu},\mathbf{y}_{i})$
$\displaystyle=\log\left(\sum_{g=1}^{G}\pi_{g}p(\mathbf{y}_{i};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})\right)$
$\displaystyle=\log\left(\sum_{g=1}^{G}\pi_{g}\int\det\left(\boldsymbol{\Theta}_{g}\right)^{\frac{1}{2}}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_{g})^{\rm
T}\boldsymbol{\Theta}_{g}(\mathbf{x}-\boldsymbol{\mu}_{g})\right){h(\mathbf{y}_{i},\mathbf{x})}d\mathbf{x}\right)+C(\mathbf{y}_{i}),$

where $C(\mathbf{y}_{i})=\sum_{j=1}^{p}\log y_{j}!-2^{-1}p\log(2\pi)$. The
log-likelihood of the MPLN model is
${\ell}_{n}(\boldsymbol{\nu})=\sum_{i=1}^{n}\ell(\boldsymbol{\nu},\mathbf{y}_{i}).$
If we define

$L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})=\int\det(\boldsymbol{\Theta}_{g})^{\frac{1}{2}}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_{g})^{\rm
T}\boldsymbol{\Theta}_{g}(\mathbf{x}-\boldsymbol{\mu}_{g})\right){h(\mathbf{y},\mathbf{x})}d\mathbf{x},$

and
$L_{M}(\boldsymbol{\nu},\mathbf{y})=\sum_{g=1}^{G}\pi_{g}L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})$,
then
${\ell}(\boldsymbol{\nu},\mathbf{y}_{i})=\log\left(L_{M}(\boldsymbol{\nu},\mathbf{y}_{i})\right)+C(\mathbf{y}_{i})$.
Note that the function $L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})$ is
proportional to the density
$p(\mathbf{y};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})$. Let
$\mathcal{L}_{n}(\boldsymbol{\nu})=-n^{-1}\ell_{n}(\boldsymbol{\nu})$. The
optimization problem (3) in the main manuscript can be written as

$\hat{\boldsymbol{\nu}}_{n}\in\operatorname*{\arg\min}_{\boldsymbol{\nu}\in{\mathcal{D}}}\left\\{\mathcal{L}_{n}(\boldsymbol{\nu})+{\lambda_{n}}\mathcal{R}(\boldsymbol{\nu})\right\\},$
(S2)

where
$\mathcal{R}(\boldsymbol{\nu})=\sum_{g=1}^{G}||\boldsymbol{\Theta}_{g}||_{1,\rm
off}.$

For the PLN model, denote the derivative (the score function) and the Hessian
matrix of its log-likelihood as

$\mathcal{S}(\boldsymbol{\Theta},\mathbf{y})=\frac{\partial\ell(\boldsymbol{\Theta},\mathbf{y})}{\partial{\rm
vech}(\boldsymbol{\Theta})},{\bf
H}(\boldsymbol{\Theta},\mathbf{y})=\frac{\partial^{2}\ell(\boldsymbol{\Theta},\mathbf{y})}{\partial{\rm
vech}(\boldsymbol{\Theta})\partial{\rm vech}(\boldsymbol{\Theta})^{\rm T}}.$

For the MPLN model, we can similarly define its score function
$\mathcal{S}^{M}(\boldsymbol{\nu},\mathbf{y})$, its Hessian matrix ${\bf
F}(\boldsymbol{\nu},\mathbf{y})$, and its Fisher information matrix
$\boldsymbol{\Gamma}^{*}=\boldsymbol{\Gamma}(\boldsymbol{\nu}^{*})$.

${\bf D}(\boldsymbol{\nu})={\rm E}_{\boldsymbol{\nu}^{*}}({\bf
F}(\boldsymbol{\nu},\mathbf{y})).$ (S3)

Note that ${\bf
D}(\boldsymbol{\nu}^{*})=-\boldsymbol{\Gamma}(\boldsymbol{\nu}^{*})$. These
notations are summarized in Table S0.

Table S0: Summary of notations. Name | Definition
---|---
The density of the PLN distribution ${\rm PLN}\left(\boldsymbol{\Theta},\boldsymbol{\mu}\right)$ | $p(\mathbf{y};\boldsymbol{\Theta},\boldsymbol{\mu})$
The log-likelihood function of the PLN at $\mathbf{y}$ | $\ell(\boldsymbol{\Theta},\mathbf{y})=\log\left(p(\mathbf{y};\boldsymbol{\Theta},\boldsymbol{\mu})\right)$
The log-likelihood function of the PLN | $\ell_{n}(\boldsymbol{\Theta})=\sum_{i=1}^{n}\ell(\boldsymbol{\Theta},\mathbf{y}_{i})$
The derivative of $\ell(\boldsymbol{\Theta},\mathbf{y})$ | $\mathcal{S}(\boldsymbol{\Theta},\mathbf{y})=\frac{\partial\ell(\boldsymbol{\Theta},\mathbf{y})}{\partial{\rm vech}(\boldsymbol{\Theta})}$
The Hessian matrix of $\ell(\boldsymbol{\Theta},\mathbf{y})$ | ${\bf H}(\boldsymbol{\Theta},\mathbf{y})=\frac{\partial^{2}\ell(\boldsymbol{\Theta},\mathbf{y})}{\partial{\rm vech}(\boldsymbol{\Theta})\partial{\rm vech}(\boldsymbol{\Theta})^{\rm T}}$
The density of the MPLN distribution MPLN $\left(\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)$ | $p(\mathbf{y};\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G})$
The log-likelihood function of the MPLN at $\mathbf{y}$ | $\ell(\boldsymbol{\nu},\mathbf{y})=\log\left(p\left(\mathbf{y};\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)\right)$
The log-likelihood function of the MPLN | $\ell_{n}(\boldsymbol{\nu})=\sum_{i=1}^{n}\ell(\boldsymbol{\Theta},\mathbf{y}_{i})$
The first part of the objective function | $\mathcal{L}_{n}(\boldsymbol{\nu})=-n^{-1}\ell_{n}(\boldsymbol{\nu})$
The second part of the objective function | $\mathcal{R}(\boldsymbol{\nu})=\sum_{g=1}^{G}||\boldsymbol{\Theta}_{g}||_{1,\rm off}$
The derivative of $\ell(\boldsymbol{\nu},\mathbf{y})$ | $\mathcal{S}^{M}(\boldsymbol{\nu},\mathbf{y})=\frac{\partial\ell(\boldsymbol{\nu},\mathbf{y})}{\partial\boldsymbol{\nu}}$
The Hessian matrix of $\ell(\boldsymbol{\nu},\mathbf{y})$ | ${\bf F}(\boldsymbol{\nu},\mathbf{y})=\frac{\partial^{2}\ell(\boldsymbol{\nu},\mathbf{y})}{\partial\boldsymbol{\nu}\partial\boldsymbol{\nu}^{\rm T}}$
The expectation of ${\bf F}(\boldsymbol{\nu},\mathbf{y})$ | ${\bf D}(\boldsymbol{\nu})={\rm E}_{\boldsymbol{\nu}^{*}}({\bf F}(\boldsymbol{\nu},\mathbf{y}))$
The Fisher information matrix of MPLN | $\boldsymbol{\Gamma}^{*}=\boldsymbol{\Gamma}(\boldsymbol{\nu}^{*})=-{\bf D}(\boldsymbol{\nu}^{*})$

Finally, we denote $\mathbb{N}^{p}$ as the set of all $p$-dimensional non-
negative integer vectors. For a vector $\mathbf{a}=(a_{1},\ldots,a_{p})$, we
denote $||\mathbf{a}||_{2}=\sqrt{\sum_{j=1}^{p}a_{j}^{2}}$ as its $L_{2}$-norm
and $||\mathbf{a}||_{\infty}=\max_{j}|a_{j}|$ as its $L_{\infty}$-norm. For a
matrix $\mathbf{A}$, we denote $||\mathbf{A}||_{2}$ as its largest singular
value of $\mathbf{A}$ and $||\mathbf{A}||_{1,\infty}$ as its maximum absolute
row sum of $\mathbf{A}$. Given $\boldsymbol{\Theta}$ and $\boldsymbol{\mu}$,
we define an operator $\mathcal{T}$ that maps functions in $x$ to functions in
$y$,

$\mathcal{T}(f)=\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)f(\mathbf{x})h(\mathbf{y},\mathbf{x})d\mathbf{x}.$

We use $\mathbb{I}(\mathbf{x})\equiv 1$ as the constant function taking value
1.

####  S1.1.2 Some Lemmas

######  Lemma S1.

Let $\mathbf{y}\sim{\rm PLN}(\boldsymbol{\Theta},\boldsymbol{\mu})$. For any
$n,y\in\mathbb{N}$, we define

$\phi(y,n)=\left\\{\begin{aligned} &1&&n=0,\\\
&y(y-1)\cdots(y-n+1)&&n>0.\end{aligned}\right.$

Then, for ${\mathbf{n}}=(n_{1},\cdots,n_{p})^{\rm T}$, we have

${\rm
E}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})\right)=\exp\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}+{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}/2\right).$

######  Lemma S2 (1-dimensional dominating function).

Suppose $\theta\in[m,M]$ ($m,M>0$), and $|\mu|\leq M$. Let

$f^{1}(y,\theta,\mu)={\int\exp\left(-\frac{1}{2}\theta(x-\mu)^{2}\right)\exp\left(-\exp(x)\right)\exp(xy)dx}.$

Then, for any large enough positive integer $y$, we have

$f^{1}(y,\theta,\mu)\geq C\exp(y\log(y+1)/2),$

where $C>0$ is constant depending on $m,M$.

######  Lemma S3 ($p$-dimensional dominating function).

Let

$f^{p}(\mathbf{y},\boldsymbol{\Theta},\boldsymbol{\mu})=\frac{\mathcal{T}\left((x_{1}-\mu_{1})^{4}\right)}{\mathcal{T}(\mathbb{I})},$
(S4)

where $\mathbf{y}=(y_{1},\dots,y_{p})\in\mathbb{N}^{p}$. Assuming
$\mathbf{y}\sim{\rm PLN}(\boldsymbol{\Theta}^{*},\boldsymbol{\mu})$, under
Condition (C1-C3), there exists a polynomial function $g(\mathbf{y})$ with a
constant $C_{p}$ only depending on $p,m$ and $M$,

$g(\mathbf{y})=\left(4\frac{||\mathbf{y}||_{2}}{m}\right)^{4}+C_{p}$

such that ${\rm E}_{\boldsymbol{\Theta}^{*}}(g(\mathbf{y}))<\infty$ and
$|f^{p}(\mathbf{y},\boldsymbol{\Theta},\boldsymbol{\mu})|\leq g(\mathbf{y})$
for any $\boldsymbol{\Theta},\boldsymbol{\mu}$ satisfying Condition (C1-C3).

######  Remark 1.

Applying the same proof as in Lemma S3, the polynomial $(x_{1}-\mu_{1})^{4}$
can be replaced by any polynomial with respect to $\mathbf{x}$, e.g.
$(x_{i}-\mu_{i})^{2}(x_{j}-\mu_{j})^{2}$ and
$(x_{i}-\mu_{i})(x_{j}-\mu_{j})(x_{k}-\mu_{k})^{2}$. Further, for any two
polynomial functions ${\psi}_{1}(\mathbf{x}),{\psi}_{2}(\mathbf{x})$ and

$f^{p}(\mathbf{y},\boldsymbol{\Theta},\boldsymbol{\mu})=\frac{\mathcal{T}\left({\psi}_{1}(\mathbf{x})\right)\mathcal{T}\left({\psi}_{2}(\mathbf{x})\right)}{\mathcal{T}^{2}(\mathbb{I})},$

there exists a polynomial function $g(\mathbf{y})$ with ${\rm
E}_{\Theta^{*}}(g(\mathbf{y}))<\infty$ such that
$|f^{p}(\mathbf{y},\boldsymbol{\Theta},\boldsymbol{\mu})|\leq g(\mathbf{y})$
for any $\boldsymbol{\Theta},\boldsymbol{\mu}$ satisfying Condition (C1-C3).

Proof of Lemma S1. By the property of conditional expectation, we have

${\rm E}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})\right)={\rm
E}_{\mathbf{x}}{\rm
E}_{\mathbf{y}}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})|\mathbf{x}\right).$

From the moments of the Poisson distribution, we have

${\rm
E}_{\mathbf{y}}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})|\mathbf{x}\right)=\prod_{j=1}^{p}\exp(n_{j}x_{j}).$

Further, since $\mathbf{x}\sim{\rm
N}(\boldsymbol{\mu},\boldsymbol{\Theta}^{-1})$, we have

${\rm E}_{\mathbf{x}}\left(\prod_{j=1}^{p}\exp(n_{j}x_{j})\right)={\rm
E}_{\mathbf{x}}(\exp({{\mathbf{n}}}^{\rm
T}\mathbf{x}))=\exp\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}+{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}/2\right),$

and the conclusion follows.

Proof of Lemma S2. Since for any $y\in\mathbb{N}$, $f^{1}(y,\theta,\mu)>0$, we
only need to consider $y$ large enough. Let

$g(x,y)=\exp\left(-\frac{1}{2}\theta(x-\mu)^{2}\right)\exp(-\exp(x))\exp(xy).$

Clearly, we have, for $y$ large enough,

$\int g(x,y)dx\geq\int_{\log(y)-1}^{\log(y)}g(x,y)dx.$

Let $t=x-\log(y)$, then

$\displaystyle\quad\int_{\log(y)-1}^{\log(y)}g(x,y)dx$
$\displaystyle=\int_{-1}^{0}\exp\left(-\frac{1}{2}\theta(\log(y)+t-\mu)^{2}\right)\exp(-\exp(\log(y)+t))\exp((\log(y)+t)y)dt$
$\displaystyle=\exp(y\log(y))\int_{-1}^{0}\exp\left(-\frac{1}{2}\theta(\log(y)+t-\mu)^{2}\right)\exp(-\exp(\log(y)+t))\exp(yt)dt$
$\displaystyle\geq\exp(y\log(y))\min_{t\in[-1,0]}\left\\{\exp\left(-\frac{1}{2}\theta(\log(y)-t-\mu)^{2}\right)\right\\}\int_{-1}^{0}\exp(y(t-\exp(t)))dt.$
$\displaystyle\geq\exp(y\log(y))\min_{t\in[-1,0]}\left\\{\exp\left(-\frac{1}{2}\theta(\log(y)-t-\mu)^{2}\right)\right\\}\int_{-1}^{0}\exp(y(-1-t^{2}/2))dt.$
$\displaystyle\geq\exp(y\log(y)-y)\min_{t\in[-1,0]}\left\\{\exp\left(-\frac{1}{2}\theta(\log(y)-t-\mu)^{2}\right)\right\\}\int_{-1}^{0}\exp(-t^{2}/2)dt/\sqrt{y}.$
$\displaystyle\geq\exp(y\log(y)-y)C_{2}\exp(-C_{3}\log^{2}(y)+C_{4}\log(y))/\sqrt{y},$
(S5)

where $C_{2}$, $C_{3}$ and $C_{4}$ are three constants only depending on
$m,M$. Since the leading order of (S1.1.2) is $\exp(y\log(y)-y)$ , when $y\geq
1$, there exists a constant $C$ depending on $m,M$ such that for any large
enough positive integer $y$, we have

$f^{1}(y,\theta,\mu)\geq C\exp(y\log(y+1)/2).$

Proof of Lemma S3. Similarly, we only need to consider $||\mathbf{y}||_{2}$
large enough. The proof consists of the following three steps.

Step 1. We first give a lower bound for the denominator of (S4). Note that

$\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)\geq\exp\left(-\frac{1}{2}M(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}(\mathbf{x}-\boldsymbol{\mu})\right).$

We have

$\displaystyle\quad\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}$
$\displaystyle\geq\int\exp\left(-\frac{1}{2}M(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}(\mathbf{x}-\boldsymbol{\mu})\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}$
$\displaystyle=\prod_{j=1}^{p}\int\exp\left(-\frac{1}{2}M(x_{j}-\mu_{j})^{2}\right)\exp(-\exp(x_{j}))\exp(x_{j}y_{j})dx_{j}$

By Lemma S2, we have, for any fixed $j$,
$\int\exp\left(-\frac{1}{2}M(x_{j}-\mu_{j})^{2}\right)\exp(-\exp(x_{j}))\exp(x_{j}y_{j})dx_{j}$
is greater than $C\exp(y_{j}\log(y_{j}+1)/2)$. Hence, we have

$\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}\geq
C^{p}\exp\left(\sum_{j}y_{j}\log(y_{j}+1)/2\right).$

Step 2. When $||\mathbf{x}-\boldsymbol{\mu}||_{2}\leq
A(\mathbf{y})=4\frac{||\mathbf{y}||_{2}}{m}$, we have
$(x_{1}-\mu_{1})^{4}\leq||\mathbf{x}-\boldsymbol{\mu}||_{2}^{4}\leq
A^{4}(\mathbf{y})$, where $x_{1}$ is the first element of $\mathbf{x}$. Then,
we have

$\frac{\int_{||\mathbf{x}-\boldsymbol{\mu}||_{2}\leq
A}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)(x_{1}-\mu_{1})^{4}h(\mathbf{y},\mathbf{x})d\mathbf{x}}{\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}}\leq
A^{4}(\mathbf{y}).$

Since $A^{4}(\mathbf{y})$ is a polynomial function of $\mathbf{y}$, we have
${\rm E}_{\boldsymbol{\Theta}^{*}}(A^{4}(\mathbf{y}))<\infty.$

Step 3. When $||\mathbf{x}-\boldsymbol{\mu}||_{2}>A(\mathbf{y})$, we have

$\displaystyle\quad\int_{||\mathbf{x}-\boldsymbol{\mu}||_{2}>A(\mathbf{y})}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)(x_{1}-\mu_{1})^{4}h(\mathbf{y},\mathbf{x})d\mathbf{x}$
$\displaystyle\leq\int_{||\mathbf{x}-\boldsymbol{\mu}||_{2}>A(\mathbf{y})}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)(x_{1}-\mu_{1})^{4}\exp\left(\sum_{j=1}^{p}x_{j}{y_{j}}\right)d\mathbf{x}$
$\displaystyle=\prod_{j=1}^{p}\exp(\mu_{j}{y_{j}})\int_{||{\mathbf{u}}||_{2}>A(\mathbf{y})}\exp\left(-\frac{1}{2}{\mathbf{u}}^{\rm
T}\boldsymbol{\Theta}{\mathbf{u}}\right)u_{1}^{4}\prod_{j=1}^{p}\exp(u_{j}{y_{j}})d{\mathbf{u}}\quad({\mathbf{u}}=\mathbf{x}-\boldsymbol{\mu})$
$\displaystyle\leq\prod_{j=1}^{p}\exp(\mu_{j}{y_{j}})\int_{||{\mathbf{u}}||_{2}>A(\mathbf{y})}\exp\left(-\frac{1}{2}m||{\mathbf{u}}||_{2}^{2}\right)||{\mathbf{u}}||_{2}^{4}\exp(||{\mathbf{u}}||_{2}||\mathbf{y}||_{2})d{\mathbf{u}},$

where the last inequality is by Cauchy’s inequality. Note that the area of the
$p$-dimensional sphere of radius $r$ is $M_{p}r^{p-1}$ with $M_{p}$ being a
constant only depending on $p$. Let $r=||{\mathbf{u}}||_{2}$. By the polar
decomposition, the $p$-dimensional integral can be rewritten as

$\displaystyle\quad\int_{||{\mathbf{u}}||_{2}>A(\mathbf{y})}\exp(-\frac{1}{2}m||{\mathbf{u}}||_{2}^{2})||{\mathbf{u}}||_{2}^{4}\exp(||{\mathbf{u}}||_{2}||\mathbf{y}||_{2})d{\mathbf{u}}$
$\displaystyle=\int_{r>A(\mathbf{y})}\exp(-\frac{1}{2}mr^{2})r^{4}\exp(r||\mathbf{y}||_{2})M_{p}r^{p-1}dr$
$\displaystyle=\int_{r>A(\mathbf{y})}\exp(-\frac{1}{2}r(mr-2||\mathbf{y}||_{2}))M_{p}r^{p+3}dr$
$\displaystyle\leq\int_{r>A(\mathbf{y})}\exp(-\frac{1}{2}r(2||\mathbf{y}||_{2}))M_{p}r^{p+3}dr\quad(r>4\frac{||\mathbf{y}||_{2}}{m})$
$\displaystyle\leq\int_{r>A(\mathbf{y})}\exp(-r)M_{p}r^{p+3}dr$
$\displaystyle\leq\int_{r>0}\exp(-r)M_{p}r^{p+3}dr$
$\displaystyle=C^{\prime}_{p},$

where $C^{\prime}_{p}$ is a constant only depending on $p$. Then, under
Condition (C1-C3), we have

$\displaystyle\quad\frac{\int_{||\mathbf{x}-\boldsymbol{\mu}||_{2}\leq
A}\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)(x_{1}-\mu_{1})^{4}h(\mathbf{y},\mathbf{x})d\mathbf{x}}{\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}}$
$\displaystyle\leq\frac{\prod_{j=1}^{p}\exp(\mu_{j}{y_{j}})C^{\prime}_{p}}{\exp\left(\sum_{j}y_{j}\log(y_{j}+1)/2\right)}\leq
C_{p},$

where $C_{p}$ is a constant depending on $p,m$ and $M$.

Finally, combining the results in Step 2 and 3, we get the dominating function

$g(\mathbf{y})=\left(4\frac{||\mathbf{y}||_{2}}{m}\right)^{4}+C_{p}$

which is a polynomial. By Lemma S1, we have ${\rm
E}_{\boldsymbol{\Theta}^{*}}(g(\mathbf{y}))<\infty$.

####  S1.1.3 Proof of Theorem 1, part I

To prove the the first conclusion of Theorem 1, we introduce the following
definition and give two lemmas.

######  Definition S1 (Good vector).

We call a vector $\boldsymbol{\xi}=(\xi_{1},\cdots,\xi_{G})^{\rm
T}\in\mathbb{R}^{G}$ as a good vector if one $\xi_{g}$ only appears once in
$\boldsymbol{\xi}$, i.e. $\xi_{g^{\prime}}\neq\xi_{g}$ for all $g^{\prime}\neq
g$. We call the index $g$ as a good index with respect to $\boldsymbol{\xi}$.

######  Lemma S4.

Let $\boldsymbol{\xi}=(\xi_{1},\cdots,\xi_{G})^{\rm T}$ be a good vector with
a good index $s$, $\boldsymbol{\sigma}=(\sigma_{1},\cdots,\sigma_{G})^{\rm T}$
satisfy $\sigma_{g}>0$ for $g=1,\cdots,G$ and
$\boldsymbol{\alpha}=(\alpha_{1},\cdots,\alpha_{G})^{\rm T}$. If for any
$z\in\mathbb{N}$

$\sum_{g=1}^{G}\alpha_{g}\exp(\xi_{g}z+\sigma_{g}z^{2}/2)=0,$

then $\alpha_{s}=0$.

######  Lemma S5.

For any $n>0$, let $\mathcal{M}_{i}\subset\mathbb{R}^{p},i=1,\cdots,n$ be $n$
linear proper subspaces. Then, there exists a non-negative integer vector
$\boldsymbol{\gamma}$ such that
$\boldsymbol{\gamma}\not\in\bigcup_{i=1}^{n}\mathcal{M}_{i}$.

######  Proposition S1.

$p\left(\mathbf{y};\boldsymbol{\Theta}_{1},\boldsymbol{\mu}_{1}\right),\cdots,p\left(\mathbf{y};\boldsymbol{\Theta}_{G},\boldsymbol{\mu}_{G}\right)$
are linearly independent for $\boldsymbol{\mu}_{g}$ ($g=1,\cdots,G$) that are
bounded and different from each other.

Proof of Theorem 1 part I. By Yakowitz and Spragins (1968), under Condition
(C1-C3), the identifiability of the MPLN model is equivalent to the linear
independence of the PLN components. Thus, we aim to prove Proposition S1. We
prove this by mathematical induction.

The independence for $G=1$ is trivial. Now we assume that Proposition S1 holds
for $G-1$. For any $\boldsymbol{\mu}_{1},\cdots,\boldsymbol{\mu}_{G}$ that are
bounded and different from each other, if we can prove that there exists
$\boldsymbol{\alpha}=(\alpha_{1},\cdots,\alpha_{G})^{\rm T}$ and an index
$s\in\\{1,\cdots,G\\}$ such that

$\sum_{g=1}^{G}\alpha_{g}p\left(\mathbf{y};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g}\right)=0\mbox{
and }\alpha_{s}=0,$

then by induction, we have $\boldsymbol{\alpha}=0$ and hence
$p(\mathbf{y};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})$ ($g=1,\cdots,G$)
are linearly independent. So our goal is to prove that if
$\sum_{g=1}^{G}\alpha_{g}p(\mathbf{y};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})=0$,
then we can always find an index $s$ such that $\alpha_{s}=0$.

Let ${{\mathbf{n}}}=(n_{1},\cdots,n_{p})^{\rm T}$ be any non-negative integer
vector. Then, for any positive integer $z$, by Lemma S1, there exists a
polynomial function
$p_{z}(\mathbf{y})=\prod_{j=1}^{p}\phi(zn_{j},y_{j}),z\in\mathbb{N}$ such that

$\sum_{g=1}^{G}\alpha_{g}{\rm E}_{g}\left(p_{z}(\mathbf{y})\right)=0\mbox{ and
}{\rm E}_{g}\left(p_{z}(\mathbf{y})\right)=\exp\left(z{\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{g}+z^{2}{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{g}{{\mathbf{n}}}/2\right),$

where ${\rm E}_{g}$ represents taking expectation with respect to ${\rm
PLN}\left(\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g}\right)$. Let

$\boldsymbol{\xi}=\left({\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{1},\cdots,{\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{G}\right)\mbox{ and
}\boldsymbol{\sigma}=\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{1}{{\mathbf{n}}}/2,\cdots,{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{G}{{\mathbf{n}}}/2\right).$

By Lemma S4, if there exists an ${\mathbf{n}}$ such that $\boldsymbol{\xi}$ is
a good vector with good index $s$, then $\alpha_{s}=0$ and we complete the
proof. If, on the other hand, $\boldsymbol{\xi}=\left({\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{1},\cdots,{\mathbf{n}}^{\rm T}\boldsymbol{\mu}_{G}\right)$
is not a good vector for any non-negative integer vector ${\mathbf{n}}$.
Therefore, for any ${\mathbf{n}}$, there exists $s\neq 1$ such that
${\mathbf{n}}^{\rm T}\boldsymbol{\mu}_{1}={\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{s}$. Thus, ${\mathbf{n}}$ is the solution to the linear
equation ${\mathbf{x}}^{\rm T}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{s})=0$.
We define $\mathcal{M}_{g}$ as the linear space consisting of solutions to the
linear equation ${\mathbf{x}}^{\rm
T}\left(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{g}\right)=0$ ($g\neq 1$) and
$\mathcal{M}=\cup_{g=2}^{G}\mathcal{M}_{g}$. Thus, for any non-negative
integer vector ${\mathbf{n}}$, we have ${\mathbf{n}}\in\mathcal{M}$. Since
$\boldsymbol{\mu}_{g}$ are different form each other, then
$\dim{\mathcal{M}_{g}}=p-1$ and $\mathcal{M}_{g}$ is a proper subspace of
$\mathbb{R}^{p}$. This is contradictory to Lemma S5. So there exists an
${\mathbf{n}}$ such that $\boldsymbol{\xi}=\left({\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{1},\cdots,{\mathbf{n}}^{\rm T}\boldsymbol{\mu}_{G}\right)$
is a good vector and hence for any $\boldsymbol{\mu}_{g}$ ($g=1,\cdots,G$)
that are different from each other,
$p(\mathbf{y};\boldsymbol{\Theta}_{1},\boldsymbol{\mu}_{1}),\cdots,p(\mathbf{y};\boldsymbol{\Theta}_{G},\boldsymbol{\mu}_{G})$
are linearly independent.

Proof of Lemma S4. Without loss of generality, we assume that
$(\xi_{g},\sigma_{g})$ ($g=1,\cdots,G$) are increasingly ordered (first by
$\xi$ then by $\sigma$). We say that $(\xi_{g},\sigma_{g})$ and
$(\xi_{s},\sigma_{s})$ are equivalent if
$(\xi_{g},\sigma_{g})=(\xi_{s},\sigma_{s})$. By this equivalence relationship,
$\\{(\xi_{g},\sigma_{g})\\}_{g=1}^{G}$ can be partitioned into $Q$ groups
($Q\geq 1$). Let $S_{q}$ be the index set of the $q$-th group. We have

$\sum_{q=1}^{Q}\sum_{j\in S_{q}}\alpha_{j}\exp(\xi_{j}z+\sigma_{j}z^{2}/2)=0$

for all $z\in\mathbb{N}$. Dividing $\exp(\xi_{G}z+\sigma_{G}z^{2}/2)$ on both
sides of the above equation, we get

$\sum_{q=1}^{Q-1}\sum_{j\in
S_{q}}\alpha_{j}\exp(\xi_{j}z+\sigma_{j}z^{2}/2-\xi_{G}z-\sigma_{G}z^{2}/2)+\sum_{j\in
S_{Q}}\alpha_{j}=0$ (S6)

for all $z\in\mathbb{N}$. By the choice of $\sigma_{G},\xi_{G}$, the first
summation of (S6) converges to zero when $z$ goes to infinity. So, we have
$\sum_{j\in S_{Q}}\alpha_{j}=0$. By mathematical induction, we have
$\sum_{j\in S_{q}}\alpha_{j}=0$ for $q=1,\cdots,Q$. Since $\boldsymbol{\xi}$
is a good vector with a good index $s$, ($\xi_{s},\sigma_{s}$) itself forms a
group, and hence $\alpha_{s}=0$.

Proof of Lemma S5. We prove by mathematical induction. The conclusion clearly
holds for $n=1$. Now we assume that Lemma S5 holds for $n$ and we aim to prove
that it also holds for $n+1$.

By induction hypothesis, we can take
$\boldsymbol{\alpha}\in\mathbb{N}^{p}\setminus\bigcup_{i=1}^{n}\mathcal{M}_{i}$.
If $\boldsymbol{\alpha}\not\in\mathcal{M}_{n+1}$, we have
$\boldsymbol{\alpha}\not\in\bigcup_{i=1}^{n+1}\mathcal{M}_{i}$, and the proof
is finished. Thus, we only need to consider
$\boldsymbol{\alpha}\in\mathcal{M}_{n+1}$. Similarly, we can take
${\boldsymbol{\beta}}\in\mathbb{N}^{p}\setminus\bigcup_{i=2}^{n+1}\mathcal{M}_{i}$
and ${\boldsymbol{\beta}}\in\mathcal{M}_{1}$. For any $i\neq 1$, we can prove
that there is at most one $k_{1}$ such that
$\boldsymbol{\alpha}+k_{1}\boldsymbol{\beta}\in\mathcal{M}_{i}$. In fact, if
there are $k_{1},k_{2}$ such that $k_{1}\neq k_{2}$ and
$\boldsymbol{\alpha}+k_{1}{\boldsymbol{\beta}}\in\mathcal{M}_{i},\boldsymbol{\alpha}+k_{2}{\boldsymbol{\beta}}\in\mathcal{M}_{i}$,
then ${\boldsymbol{\beta}}\in\mathcal{M}_{i}$, which is contradictory to the
fact that
${\boldsymbol{\beta}}\in\mathbb{N}^{p}\setminus\bigcup_{i=2}^{n+1}\mathcal{M}_{i}$.
Furthermore, there is no $k\in\mathbb{N}$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\in\mathcal{M}_{1}$. If otherwise,
there exists a $k\in\mathbb{N}$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\in\mathcal{M}_{1}$. Then, we have
$\boldsymbol{\alpha}\in\mathcal{M}_{1}$, which is also a contradiction. So we
could find at most $n$ positive integers $k$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\in\bigcup_{i=1}^{n+1}\mathcal{M}_{i}$.
Since there are infinitely many non-negative numbers, we prove that there
exists $k\in\mathbb{N}$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\not\in\bigcup_{i=1}^{n+1}\mathcal{M}_{i}$,
and Lemma S5 is proved.

####  S1.1.4 Proof of Theorem 1, part II

To prove this conclusion, we need to give the explicit formula for the score
functions and the Fisher information matrices of the PLN and MPLN. The Hessian
matrix ${\bf H}\left(\boldsymbol{\Theta},\mathbf{y}\right)$ of the PLN is a
$p(p+1)/2\times p(p+1)/2$ matrix. For notational convenience, we let ${\bf
H}_{(i,j)(i^{\prime},j^{\prime})}=\frac{\partial^{2}\ell(\boldsymbol{\Theta},\mathbf{y})}{\partial\Theta_{i^{\prime}j^{\prime}}\partial\Theta_{ij}}$,
$(i\leq j,i^{\prime}\leq j^{\prime})$ as the element at the $(2p-i+1)i/2-p+j$
row and $(2p-i^{\prime}+1)i^{\prime}/2-p+j^{\prime}$ column of the Hessian
matrix. It is clear that the densities of the PLN and MPLN satisfy the
regularity conditions in Shao (2003). Then, we can calculate the score
function and the Fisher information as follows. The score function of the PLN
can be written as

$\mathcal{S}(\boldsymbol{\Theta},\mathbf{y})=\frac{1}{2}{\rm
vech_{2}}\left(\boldsymbol{\Theta}^{-1}\right)-\frac{1}{2}\frac{\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right){\rm
vech_{2}}\left((\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}}{\int\exp\left(-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\boldsymbol{\Theta}(\mathbf{x}-\boldsymbol{\mu})\right)h(\mathbf{y},\mathbf{x})d\mathbf{x}}.$

Especially, at the true parameter $\boldsymbol{\Theta}^{*}$ , we have
$\mathcal{S}(\boldsymbol{\Theta}^{*},\mathbf{y})=\frac{1}{2}{\rm
E}_{\mathbf{x}}\left({\rm
vech_{2}}\left(\boldsymbol{\Theta}^{*-1}-(\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\right)|\mathbf{y}\right)$. We use $(i,j),i\leq j$ to index $\mathcal{S}$.
Using the operator $\mathcal{T}$, the element of the score function at
$(2p-i+1)i/2-p+j$ can be rewritten as

$\mathcal{S}_{(i,j)}\left(\boldsymbol{\Theta},\mathbf{y}\right)=\frac{1}{2}{\rm
vech_{2}}\left(\boldsymbol{\Theta}^{-1}\right)_{(i,j)}-\frac{1}{2}\frac{\mathcal{T}\left({\rm
vech_{2}}\left((\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\right)_{(i,j)}\right)}{\mathcal{T}\left(\mathbb{I}\right)}.$

Using the operator $\mathcal{T}$, the Fisher information matrix can be written
as follows. Let $\Sigma=\boldsymbol{\Theta}^{-1}.$

When $i=j,i^{\prime}=j^{\prime}$,

$\displaystyle{\bf
H}_{(i,i)(i^{\prime},i^{\prime})}\left(\boldsymbol{\Theta},\mathbf{y}\right)$
$\displaystyle=-\frac{1}{2}\Sigma_{ii^{\prime}}\Sigma_{i^{\prime}i}+\frac{1}{4}\frac{\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i}^{2}({\mathbf{x}-\boldsymbol{\mu}})_{i^{\prime}}^{2}\right)}{\mathcal{T}\left(\mathbb{I}\right)}$
$\displaystyle-\frac{1}{4}\frac{\mathcal{T}\left(({\mathbf{x}-\boldsymbol{\mu}})_{i^{\prime}}^{2}\right)\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i}^{2}\right)}{\mathcal{T}^{2}\left(\mathbb{I}\right)}.$

When $i\neq j,i^{\prime}=j^{\prime}$,

$\displaystyle{\bf
H}_{(i,j)(i^{\prime},i^{\prime})}\left(\boldsymbol{\Theta},\mathbf{y}\right)$
$\displaystyle=-\Sigma_{ii^{\prime}}\Sigma_{i^{\prime}j}+\frac{1}{2}\frac{\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i}{({\mathbf{x}-\boldsymbol{\mu}})}_{j}({\mathbf{x}-\boldsymbol{\mu}})_{i^{\prime}}^{2}\right)}{\mathcal{T}\left(\mathbb{I}\right)}$
$\displaystyle-\frac{1}{2}\frac{\mathcal{T}\left(({\mathbf{x}-\boldsymbol{\mu}})_{i^{\prime}}^{2}\right)\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i}{({\mathbf{x}-\boldsymbol{\mu}})}_{j}\right)}{\mathcal{T}^{2}\left(\mathbb{I}\right)}.$

When $i=j,i^{\prime}\neq j^{\prime}$,

$\displaystyle{\bf
H}_{(i,i)(i^{\prime},j^{\prime})}\left(\boldsymbol{\Theta},\mathbf{y}\right)$
$\displaystyle=-\Sigma_{ii^{\prime}}\Sigma_{j^{\prime}i}+\frac{1}{2}\frac{\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i^{\prime}}{({\mathbf{x}-\boldsymbol{\mu}})}_{j^{\prime}}({\mathbf{x}-\boldsymbol{\mu}})_{i}^{2}\right)}{\mathcal{T}\left(\mathbb{I}\right))}$
$\displaystyle-\frac{1}{2}\frac{\mathcal{T}\left(({\mathbf{x}-\boldsymbol{\mu}})_{i}^{2}\right)\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i^{\prime}}{({\mathbf{x}-\boldsymbol{\mu}})}_{j^{\prime}}\right)}{\mathcal{T}^{2}(\mathbb{I})}.$

When $i\neq j,i^{\prime}\neq j^{\prime}$,

$\displaystyle{\bf
H}_{(i,j)(i^{\prime},j^{\prime})}\left(\boldsymbol{\Theta},\mathbf{y}\right)$
$\displaystyle=-(\Sigma_{ii^{\prime}}\Sigma_{j^{\prime}j}+\Sigma_{ij^{\prime}}\Sigma_{i^{\prime}j})+\frac{\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i^{\prime}}{({\mathbf{x}-\boldsymbol{\mu}})}_{j^{\prime}}({\mathbf{x}-\boldsymbol{\mu}})_{i}({\mathbf{x}-\boldsymbol{\mu}})_{j}\right)}{\mathcal{T}\left(\mathbb{I}\right)}$
$\displaystyle-\frac{\mathcal{T}\left(({\mathbf{x}-\boldsymbol{\mu}})_{i}({\mathbf{x}-\boldsymbol{\mu}})_{j}\right)\mathcal{T}\left({({\mathbf{x}-\boldsymbol{\mu}})}_{i^{\prime}}{({\mathbf{x}-\boldsymbol{\mu}})}_{j^{\prime}}\right)}{\mathcal{T}^{2}\left(\mathbb{I}\right)}.$

######  Lemma S6.

Assume $\mathbf{y}\sim{\rm PLN}(\boldsymbol{\Theta},\boldsymbol{\mu})$. Under
Condition (C1-C3), there exist two polynomial functions
$K_{1}(\mathbf{y}),K_{2}(\mathbf{y})$ with ${\rm E}[K_{1}(y)]<\infty$ and
${\rm E}[K_{2}(y)]<\infty$ such that for any $i,j,i^{\prime},j^{\prime}$,
$|\mathcal{S}_{(i,j)}(\boldsymbol{\Theta},\mathbf{y})|\leq K_{1}(\mathbf{y})$,
$|{\bf H}_{(i,i)(i^{\prime},i^{\prime})}(\boldsymbol{\Theta},\mathbf{y})|\leq
K_{2}(\mathbf{y})$.

Now we consider the score function and the Fisher information matrix of the
MPLN. The score function of the MPLN can be written as

$\mathcal{S}^{M}(\boldsymbol{\nu},\mathbf{y})=L_{M}(\boldsymbol{\nu},\mathbf{y})^{-1}\left({\pi_{1}\frac{\partial
L_{1}(\boldsymbol{\nu}_{1},\mathbf{y})}{\partial\boldsymbol{\nu}_{1}}},\ldots,{\pi_{G}\frac{\partial
L_{G}(\boldsymbol{\nu}_{G},\mathbf{y})}{\partial\boldsymbol{\nu}_{G}}}\right):=\left(\mathcal{S}^{M}_{1}(\boldsymbol{\nu},\mathbf{y}),\cdots,\mathcal{S}^{M}_{G}(\boldsymbol{\nu},\mathbf{y})\right).$

Similarly, we use $(g,i,j),i\leq j$ to index $\mathcal{S}^{M}$. Recall that
$\boldsymbol{\nu}_{g}={\rm vech}(\boldsymbol{\Theta}_{g})$. We use
$(g,i,j),(g^{\prime},i^{\prime},j^{\prime}),i\leq j,i^{\prime}\leq j^{\prime}$
to index the element at the $(g-1)p(p+1)/2+(2p-i+1)i/2-p+j$ row and
$(g^{\prime}-1)p(p+1)/2+(2p-i^{\prime}+1)i^{\prime}/2-p+j^{\prime}$ column of
$F(\boldsymbol{\nu},\mathbf{y})$, respectively.

When $g=g^{\prime}$

$\displaystyle{\bf
F}_{(g,i,j,g,i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y})$
$\displaystyle=\frac{\partial}{{\partial\boldsymbol{\Theta}_{g,i^{\prime}j^{\prime}}}}\mathcal{S}^{M}_{(g,i,j)}(\boldsymbol{\nu},\mathbf{y})$
$\displaystyle=\frac{\pi_{g}L_{g}(\boldsymbol{\Theta}_{g},\mathbf{y})}{L_{M}(\boldsymbol{\nu},\mathbf{y})}\left(H_{(i,j)(i^{\prime},j^{\prime})}(\boldsymbol{\Theta}_{g},\mathbf{y})+\mathcal{S}_{(i^{\prime},j^{\prime})}(\boldsymbol{\Theta}_{g},\mathbf{y})\mathcal{S}_{(i,j)}(\boldsymbol{\Theta}_{g},\mathbf{y})\right)$
$\displaystyle-\frac{(\pi_{g}L_{g}(\boldsymbol{\Theta}_{g},\mathbf{y}))^{2}}{L_{M}(\boldsymbol{\nu},\mathbf{y})^{2}}\mathcal{S}_{(i^{\prime},j^{\prime})}(\boldsymbol{\Theta}_{g},\mathbf{y})\mathcal{S}_{(i,j)}(\boldsymbol{\Theta}_{g},\mathbf{y}).$
(S7)

When $g\neq g^{\prime}$, we have

$\displaystyle{\bf
F}_{(g,i,j,g^{\prime},i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y})$
$\displaystyle=\frac{\partial}{{\partial\boldsymbol{\Theta}_{g^{\prime},i^{\prime}j^{\prime}}}}\mathcal{S}^{M}_{(g,i,j)}(\boldsymbol{\nu},\mathbf{y})$
$\displaystyle=-\frac{\pi_{g}\pi_{g}^{\prime}L_{g}(\boldsymbol{\Theta}_{g},\mathbf{y})L_{g^{\prime}}(\boldsymbol{\Theta}_{g^{\prime}},\mathbf{y})}{L_{M}(\boldsymbol{\nu},\mathbf{y})^{2}}\mathcal{S}_{(i^{\prime},j^{\prime})}(\boldsymbol{\Theta}_{g^{\prime}},\mathbf{y})\mathcal{S}_{(i,j)}(\boldsymbol{\Theta}_{g},\mathbf{y}).$
(S8)

Recall the definition (S3) of ${\bf D}(\boldsymbol{\nu})$. Then,

$\boldsymbol{\Gamma}(\boldsymbol{\nu}^{*})=-{\bf D}(\boldsymbol{\nu}^{*})={\rm
E}\left(\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y})\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y})^{\rm
T}\right)$

is the Fisher information matrix of the MPLN at $\boldsymbol{\nu}^{*}$.

######  Lemma S7.

Assume $\mathbf{y}\sim{\rm MPLN}(\boldsymbol{\nu},\boldsymbol{\mu})$. Under
Condition (C1-C3), there exists a polynomial function $K(\mathbf{y})$ with
${\rm E}[K(y)]<\infty$ such that for any
$g,i,j,g^{\prime},i^{\prime},j^{\prime}$, $|{\bf
F}_{(g,i,j,g^{\prime},i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y})|\leq
K(\mathbf{y})$.

In addition, we require the following two lemmas.

######  Lemma S8.

Let $\mathbf{y}$ and $\mathbf{x}$ be random variables as in the PLN model (1)
in the main manuscript of the paper and $\phi(y,n)$ is the same as Lemma S1.
Let ${{\mathbf{n}}}=(n_{1},\cdots,n_{p})^{\rm T}$ and $\mathbf{T}$ be a
$p\times p$ matrix. $\mathbf{T}\in\mathbb{R}^{p\times p}$. We have

$\displaystyle\quad{\rm
E}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})\operatorname{tr}\left(\mathbf{T}\left(\boldsymbol{\Theta}^{-1}-({\mathbf{x}}-\boldsymbol{\mu})({\mathbf{x}}-\boldsymbol{\mu})^{\rm
T}\right)\right)\right)$ $\displaystyle=\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}\mathbf{T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}\right)\exp\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}+{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}/2\right).$

######  Lemma S9.

For any $n>0$, let $\mathcal{M}_{i}\subset\mathbb{R}^{p},i=1,\cdots,n$ be $n$
linear proper subspaces. Let $A\neq{\bf 0}$ be a $p\times p$ symmetric matrix
and $\mathcal{V}=\\{\mathbf{y}:\mathbf{y}^{\rm T}A\mathbf{y}=0\\}$. Then,
there exists a non-negative integer vector $\boldsymbol{\gamma}$ such that
$\boldsymbol{\gamma}\not\in\bigcup_{i=1}^{n}\mathcal{M}_{i}$ and
$\boldsymbol{\gamma}\not\in\mathcal{V}$.

Proof of Theorem 1, part II. Note that

$\boldsymbol{\Gamma}(\boldsymbol{\nu}^{*})={\rm
E}_{\mathbf{y}}(\mathcal{S}_{M}(\boldsymbol{\nu}^{*},\mathbf{y})\mathcal{S}_{M}(\boldsymbol{\nu}^{*},\mathbf{y})^{\rm
T}).$

If there exists a non-zero vector
$\mathbf{t}=(\mathbf{t}_{1},\cdots,\mathbf{t}_{G})^{\rm T}$ such that ${\rm
E}(\mathbf{t}^{\rm
T}\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y})\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y})^{\rm
T}\mathbf{t})=0$, then we aim to prove that $\mathbf{t}={\bf 0}$. Since y is a
discrete random variable, it follows that for any $\mathbf{y}$,
$\mathbf{t}^{\rm T}\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y})=0$. Then,
we have

${L_{M}(\boldsymbol{\nu},\mathbf{y})}^{-1}\sum_{g=1}^{G}\mathbf{t}_{g}{\pi_{g}\frac{\partial
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})}{\partial\boldsymbol{\nu}_{g}}}=0.$

Since $L_{M}(\boldsymbol{\nu},\mathbf{y})\neq 0$, we have

$\sum_{g=1}^{G}\mathbf{t}_{g}^{\rm T}{\pi_{g}\frac{\partial
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})}{\partial\boldsymbol{\nu}_{g}}}=0.$

Let $\psi(\mathbf{y})=\prod_{j=1}^{p}\phi(zn_{j},y_{j}),z\in\mathbb{N}^{p}$.
Then

$\sum_{g=1}^{G}\mathbf{t}_{g}^{\rm T}{\pi_{g}\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})}{\partial\boldsymbol{\nu}_{g}}}L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})\psi(\mathbf{y})=0.$

Since $L_{r}(\boldsymbol{\nu}_{r},\mathbf{y})$ is proportional to the density
of the PLN with parameters $\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g}$, the
above equation can be rewritten as

$\sum_{g=1}^{G}\mathbf{t}_{g}^{\rm T}{\pi_{g}\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})}{\partial\boldsymbol{\nu}_{g}}}p({\mathbf{y}};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})\psi(\mathbf{y})=0.$

Summing over $\mathbf{y}$, we get

$\sum_{\mathbf{y}}\sum_{g=1}^{G}\mathbf{t}_{g}^{\rm
T}{\pi_{g}\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})}{\partial\boldsymbol{\nu}_{g}}}p({\mathbf{y}};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})\psi(\mathbf{y})=0.$

By Fubini ’s Theorem, we get

$\sum_{g=1}^{G}\sum_{\mathbf{y}}\mathbf{t}_{g}^{\rm
T}{\pi_{g}\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y})}{\partial\boldsymbol{\nu}_{g}}}p({\mathbf{y}};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})\psi(\mathbf{y})=0.$
(S9)

Then, let ${\mathbf{n}}=(n_{1},\cdots,n_{p})^{\rm T}$ and $\mathbf{T}_{g}$ be
the symmetric matrix such that ${\rm vech}(\mathbf{T}_{g})=\mathbf{t}_{g}$.
For a fixed $g$, we have

$\displaystyle\quad\sum_{{\mathbf{y}}}\mathbf{t}_{g}^{\rm
T}{\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},{\mathbf{y}})}{\partial\boldsymbol{\nu}_{g}}}p({\mathbf{y}};\boldsymbol{\Theta}_{g},\boldsymbol{\mu}_{g})\psi({\mathbf{y}})={\rm
E}_{{\mathbf{y}}_{g}}\left(\psi(\mathbf{y}_{g})\mathbf{t}_{g}^{\rm
T}{\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y}_{g})}{\partial\boldsymbol{\nu}_{g}}}\right)$
$\displaystyle=2^{-1}{\rm
E}_{\mathbf{x}_{g},\mathbf{y}_{g}}\left(\psi(\mathbf{y}_{g})\operatorname{tr}\left\\{\mathbf{T}_{g}\boldsymbol{\Theta}^{*-1}_{g}-\mathbf{T}_{g}(\mathbf{x}_{g}-\boldsymbol{\mu})(\mathbf{x}_{g}-\boldsymbol{\mu})^{\rm
T}\right\\}\right),$

where $\mathbf{y}_{g}$ follows the PLN distribution with parameters
$\boldsymbol{\Theta}_{g}$ and $\boldsymbol{\mu}_{g}$, $\mathbf{x}_{g}\sim{\rm
N}\left(\boldsymbol{\mu}_{g},\boldsymbol{\Theta}_{g}^{-1}\right)$ is the
corresponding latent variable. By Lemma S8, we get

${\rm
E}_{\mathbf{x}_{g},\mathbf{y}_{g}}\left(\psi(\mathbf{y}_{g})\mathbf{t}_{g}^{\rm
T}{\frac{\partial\log
L_{g}(\boldsymbol{\nu}_{g},\mathbf{y}_{g})}{\partial\boldsymbol{\nu}_{g}}}\right)=z^{2}\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{g}\mathbf{T}_{g}\boldsymbol{\Theta}^{-1}_{g}{{\mathbf{n}}}\right)\exp\left(z{{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}_{g}+z^{2}{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}_{g}^{-1}{{\mathbf{n}}}/2\right).$

Then, (S9) can be rewritten as, for all $z\in\mathbb{N}$,

$\sum_{g=1}^{G}\pi_{g}z^{2}\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{g}\mathbf{T}_{g}\boldsymbol{\Theta}^{-1}_{g}{{\mathbf{n}}}\right)\exp\left(z{{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}_{g}+z^{2}{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}_{g}^{-1}{{\mathbf{n}}}/2\right)=0.$

In order to show $\mathbf{T}_{1}=0$, similar to the proof of the first
conclusion, we define $\mathcal{M}_{g}$ as linear space consisting of
solutions to the linear equation ${\mathbf{x}}^{\rm
T}(\boldsymbol{\mu}_{1}-\boldsymbol{\mu}_{g})=0$ ($g=2,\cdots,G$) and
$\mathcal{M}=\cup_{g=2}^{G}\mathcal{M}_{g}$. For any
${\mathbf{n}}\not\in\mathcal{M}$, then $({\mathbf{n}}^{\rm
T}\boldsymbol{\mu}_{1},\ldots,{\mathbf{n}}^{\rm T}\boldsymbol{\mu}_{G})$ is a
good vector with a good index $1$. Since $\pi_{1}>0$, we must have
${{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{1}\mathbf{T}_{1}\boldsymbol{\Theta}^{-1}_{1}{{\mathbf{n}}}=0$.
By Lemma S9, if
$\boldsymbol{\Theta}^{-1}_{1}\mathbf{T}_{1}\boldsymbol{\Theta}^{-1}_{1}$ is
not a zero matrix, then there exists an ${{\mathbf{n}}}$ such that
${\mathbf{n}}\not\in\mathcal{M}$ and ${{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{1}\mathbf{T}_{1}\boldsymbol{\Theta}^{-1}_{1}{{\mathbf{n}}}\neq
0$, which is contradictory to the fact that ${\mathbf{n}}\not\in\mathcal{M}$
implies ${{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}_{1}\mathbf{T}_{1}\boldsymbol{\Theta}^{-1}_{1}{{\mathbf{n}}}=0$.
Hence, we must have
$\boldsymbol{\Theta}^{-1}_{1}\mathbf{T}_{1}\boldsymbol{\Theta}^{-1}_{1}={\bf{0}}$
and thus $\mathbf{T}_{1}={\bf{0}}$. Similarly, we get
$\mathbf{T}_{g}={\bf{0}}$ for all $g=1,\ldots,G$. It follows that
$\mathbf{t}={\bf 0}$, and we compete the proof.

Proof of Lemma S6.  We only prove that there is a dominating function for
$\mathcal{S}_{(i,j)}(\boldsymbol{\Theta},\mathbf{y})$. Others can be proved
similarly. Since
$m\leq\lambda_{min}(\boldsymbol{\Theta})\leq\lambda_{max}(\boldsymbol{\Theta})\leq
M$, we have $\lambda_{max}(\boldsymbol{\Sigma})\leq 1/m$ and thus ${\rm
vech_{2}}(\boldsymbol{\Theta}^{-1})_{(i,j)}$ is bounded by a constant.
Further, since ${\rm
vech_{2}}((\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{\rm
T})_{(i,j)}$ is a polynomial function of $\mathbf{x}$, by Remark 1, we have

$\frac{1}{2}\frac{\mathcal{T}\left({\rm
vech_{2}}\left((\mathbf{x}-\boldsymbol{\mu})(\mathbf{x}-\boldsymbol{\mu})^{\rm
T}\right)_{(i,j)}\right)}{\mathcal{T}\left(\mathbb{I}\right)}$

can be bounded by an integrable polynomial function and we prove the existence
of $K_{1}(\mathbf{y})$.

Proof of Lemma S7. If $g=g^{\prime}$, since
$\pi_{g}L_{g}(\boldsymbol{\nu},\mathbf{y})/L_{M}(\boldsymbol{\nu},\mathbf{y})\leq
1$, then we have,

$|{\bf
F}_{(g,i,j,g,i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y})|\leq|{\bf
H}_{(i,j,i^{\prime}j^{\prime})}(\boldsymbol{\Theta}_{g},\mathbf{y})|+2|\mathcal{S}_{(i,j)}(\boldsymbol{\Theta}_{g},\mathbf{y})\mathcal{S}_{(i^{\prime},j^{\prime})}(\boldsymbol{\Theta}_{g},\mathbf{y})|.$

The, by Lemma S6, we know that, there exists a function $K(\mathbf{y})$ such
that $|{\bf
F}_{(g,i,j,g,i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y})|\leq
K(\mathbf{y})$ and ${\rm E}(K(\mathbf{y}))<\infty$. The same proof can be
applied to the $g\neq g^{\prime}$ case.

Proof of Lemma S8. By Lemma S1

${\rm
E}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})\operatorname{tr}\left(\mathbf{T}\boldsymbol{\Theta}^{-1}\right)\right)=\operatorname{tr}\left(\mathbf{T}\boldsymbol{\Theta}^{-1}\right)\exp\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}+{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}/2\right).$

Similar to the proof of Lemma S1, by the moment generating function of the
normal distribution, we have

$\displaystyle\quad{\rm
E}\left(\prod_{j=1}^{p}\phi({n_{j}},y_{j})\operatorname{tr}\left(\mathbf{T}({\mathbf{x}}-\boldsymbol{\mu})({\mathbf{x}}-\boldsymbol{\mu})^{\rm
T}\right)\right)$ $\displaystyle={\rm
E}_{\mathbf{x}}\left(\exp({\mathbf{n}}^{\rm
T}\mathbf{x})\operatorname{tr}\left(\mathbf{T}({\mathbf{x}}-\boldsymbol{\mu})({\mathbf{x}}-\boldsymbol{\mu})^{\rm
T}\right)\right)$
$\displaystyle=\left\\{\operatorname{tr}\left(\mathbf{T}\boldsymbol{\Theta}^{-1}\right)+{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}\mathbf{T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}\right\\}\exp\left({{\mathbf{n}}}^{\rm
T}\boldsymbol{\mu}+{{\mathbf{n}}}^{\rm
T}\boldsymbol{\Theta}^{-1}{{\mathbf{n}}}/2\right).$

Lemma S8 follows from the above two equations.

Proof of Lemma S9. By the proof of Lemma S5, there exists
$\boldsymbol{\alpha}\not\in\bigcup_{i=1}^{n}\mathcal{M}_{i}$. We can assume
$\boldsymbol{\alpha}\in\mathcal{V}$. Otherwise, we complete the proof. Since
$A$ is not a zero matrix, we can take
${\boldsymbol{\beta}}\not\in\mathcal{V}$. On the one hand, since
${\boldsymbol{\beta}}\not\in\mathcal{V}$, there are at most two integers
$k\in\mathbb{N}$ satisfying the quadratic equation
$(\boldsymbol{\alpha}+k{\boldsymbol{\beta}})^{\rm
T}A(\boldsymbol{\alpha}+k{\boldsymbol{\beta}})=0$. On the other hand, for any
$\mathcal{M}_{i},i=1,\ldots,n$, there is at most one integer $k$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\in\mathcal{M}_{i}$. Otherwise, if
there exist $k_{1}\neq k_{2}$ and $i$ satisfying

$\boldsymbol{\alpha}+k_{1}{\boldsymbol{\beta}}\in\mathcal{M}_{i},\boldsymbol{\alpha}+k_{2}{\boldsymbol{\beta}}\in\mathcal{M}_{i},$

then we have $(k_{2}-k_{1}){\boldsymbol{\beta}}\in\mathcal{M}_{i}$. It follows
that $\boldsymbol{\alpha}\in\mathcal{M}_{i}$, which is contradictory to the
fact that $\boldsymbol{\alpha}\not\in\bigcup_{i=1}^{n}\mathcal{M}_{i}$. Hence,
there are at most $(n+2)$ $k$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\in\bigcup_{i=1}^{n}\mathcal{M}_{i}\bigcup\mathcal{V}$.
Since there are infinitely many non-negative integers, there exists an integer
$k$ such that
$\boldsymbol{\alpha}+k{\boldsymbol{\beta}}\not\in\bigcup_{i=1}^{n}\mathcal{M}_{i}\bigcup\mathcal{V}$.

####  S1.1.5 Proof of Theorem 2

Since we only prove the positive definiteness of the Fisher information, we
can only get the local strong convexity. In order to prove the convergence
rate, we first give the consistency of MLE. The proof is based on the
M-estimator theory. For convenience, we introduce some notations in the
M-estimator theory. We let
$m_{\boldsymbol{\nu}}(\mathbf{y}_{i})=\ell(\boldsymbol{\nu},\mathbf{y}_{i})$,
$M_{n}(\boldsymbol{\nu})=n^{-1}\sum_{i=1}^{n}m_{\boldsymbol{\nu}}(\mathbf{y}_{i})=-\mathcal{L}_{n}(\boldsymbol{\nu})$
and $M(\boldsymbol{\nu})={\rm
E}_{\boldsymbol{\nu}^{*}}[m_{\boldsymbol{\nu}}(\mathbf{y}_{i})]$. Let
$\mathcal{D}_{0}=\\{\boldsymbol{\nu}_{0}\in{\mathcal{D}}\ |\ {\rm
E}[m_{\boldsymbol{\nu}_{0}}]=\sup_{\boldsymbol{\nu}}{\rm
E}[m_{\boldsymbol{\nu}}]\\}$. By Jessen’s inequality and the identifiablity of
the MPLN, we get that $\mathcal{D}_{0}$ only contains one element
$\boldsymbol{\nu}^{*}$. Then we need the following two conditions and Lemma
S10. The proof of Lemma S10 can be found in Van der Vaart (2000).

  * (SC1)

$\limsup_{\boldsymbol{\nu}_{n}\rightarrow\boldsymbol{\nu}}m_{\boldsymbol{\nu}_{n}}(\mathbf{y})\leq
m_{\boldsymbol{\nu}}(\mathbf{y})\ \mbox{ for all }\boldsymbol{\nu}\mbox{ and
a.s. }\mathbf{y}$.

  * (SC2)

For all sufficiently small ball $U\subset{\mathcal{D}}$,
$\mathbf{y}\mapsto\sup_{\boldsymbol{\nu}\in
U}m_{\boldsymbol{\nu}}(\mathbf{y})$ is measurable and satisfies ${\rm
E}(\sup_{\boldsymbol{\nu}\in U}m_{\boldsymbol{\nu}}(\mathbf{y}))<\infty.$

######  Lemma S10 (Wald’s consistency).

Assume that Condition (SC1-SC2) hold for $m_{\boldsymbol{\nu}}(\mathbf{y})$.
Suppose that $\hat{\boldsymbol{\nu}}_{n}$ is any sequence of random vectors
such that $M_{n}(\hat{\boldsymbol{\nu}}_{n})\geq
M_{n}(\boldsymbol{\nu}_{0})-o_{p}(1)$ for some
$\boldsymbol{\nu}_{0}\in\mathcal{D}_{0}$. Then for any $\epsilon>0$, and every
compact set $K\subset{\mathcal{D}}$, as $n\rightarrow\infty$, we have

$\mathbb{P}(d(\hat{\boldsymbol{\nu}}_{n},\mathcal{D}_{0})\geq\epsilon\wedge\hat{\boldsymbol{\nu}}_{n}\in
K)\rightarrow 0,$

where
$d(\hat{\boldsymbol{\nu}}_{n},\mathcal{D}_{0})=\inf_{\boldsymbol{\nu}_{0}\in\mathcal{D}_{0}}||\hat{\boldsymbol{\nu}}_{n}-\boldsymbol{\nu}_{0}||_{2}$.

######  Lemma S11.

Assume that the estimator $\hat{\boldsymbol{\nu}}_{n}$ minimizes (S2) in
parameter space ${\mathcal{D}}$ and $\lambda_{n}$ goes to zero. Then, for any
$\epsilon>0$, as $n\rightarrow\infty$, we have:

$\mathbb{P}(||\hat{\boldsymbol{\nu}}_{n}-\boldsymbol{\nu}^{*}||_{2}\geq\epsilon)\rightarrow
0,\mbox{ as }n\rightarrow\infty.$

######  Lemma S12 (Uniform law of large numbers).

For all $g,i,j,g^{\prime},i^{\prime},j^{\prime}$ we have

$\mathbb{P}\left\\{\lim_{n\rightarrow\infty}\sup_{\boldsymbol{\nu}\in{\mathcal{D}}}\left|\frac{1}{n}\sum_{k=1}^{n}{\bf
F}_{(g,i,j)(g^{\prime},i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y}_{k})-{\bf
D}_{(g,i,j)(g^{\prime},i^{\prime},j^{\prime})}(\boldsymbol{\nu})\right|=0\right\\}=1,$

where ${\bf D}(\boldsymbol{\nu})={\rm E}_{\boldsymbol{\nu}^{*}}({\bf
F}(\boldsymbol{\nu},\mathbf{y})).$ Furthermore, we have

$\mathbb{P}\left\\{\lim_{n\rightarrow\infty}\sup_{\boldsymbol{\nu}\in{\mathcal{D}}}\left\|\frac{1}{n}\sum_{k=1}^{n}{\bf
F}(\boldsymbol{\nu},\mathbf{y}_{k})-{\bf
D}(\boldsymbol{\nu})\right\|_{2}=0\right\\}=1.$

######  Remark 2.

By the proof of Lemma S12, we know that there exists a function
$F_{0}(\mathbf{y})$ such that for any $\boldsymbol{\nu}$, $|{\bf
F}_{(g,i,j,g^{\prime},i^{\prime},j^{\prime})}(\nu,\mathbf{y})|\leq
F_{0}(\mathbf{y})$ and ${\rm E}(F_{0}(\mathbf{y}))<\infty$ . Then, by the
dominated convergence theorem, ${\bf D}(\boldsymbol{\nu})={\rm
E}_{\boldsymbol{\nu}^{*}}({\bf F}(\boldsymbol{\nu},\mathbf{y}))$ is
continuous.

Let $\Delta_{n}=\hat{\boldsymbol{\nu}}_{n}-\boldsymbol{\nu}^{*}$ . We define

$\delta\mathcal{L}_{n}=\mathcal{L}_{n}(\boldsymbol{\nu}^{*}+\Delta_{n})-\mathcal{L}_{n}(\boldsymbol{\nu}^{*})-<\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*}),{\Delta_{n}}>,$

where
$<\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*}),{\Delta_{n}}>=\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*})^{\rm
T}{\Delta_{n}}$. According to Lemma S12, we can prove the following lemma.

######  Lemma S13.

Under Condition (C1-C3), with high probability, we have

$\delta\mathcal{L}_{n}\geq\frac{\kappa}{3}||{\Delta_{n}}||_{2}^{2},$

where $\kappa=\lambda_{min}(\boldsymbol{\Gamma}^{*})$.

Proof of Theorem 2. Define

$\mathcal{F}(\Delta)=\mathcal{L}_{n}(\boldsymbol{\nu}^{*}+\Delta)-\mathcal{L}_{n}(\boldsymbol{\nu}^{*})+\lambda_{n}(\mathcal{R}(\boldsymbol{\nu}^{*}+\Delta)-\mathcal{R}(\boldsymbol{\nu}^{*})).$

Since $\mathcal{F}(0)=0$ and $\hat{\boldsymbol{\nu}}_{n}$ minimizes (S2), we
must have $\mathcal{F}(\Delta_{n})\leq 0.$ Further, by Lemma S13, with high
probability,
$\delta\mathcal{L}_{n}\geq\frac{\kappa}{3}||\Delta_{n}||_{2}^{2}$. Then,
combining with Cauchy’s inequality and the convexity of the lasso penalty,
with high probability, we get

$\displaystyle\mathcal{F}(\Delta_{n})$
$\displaystyle=\delta\mathcal{L}_{n}+<\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*}),\Delta_{n}>+{\lambda_{n}}(\mathcal{R}(\boldsymbol{\nu}^{*}+\Delta_{n})-\mathcal{R}(\boldsymbol{\nu}^{*}))$
$\displaystyle\geq\frac{\kappa}{3}||\Delta_{n}||_{2}^{2}-(||\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*})||_{\infty}+{2\lambda_{n}})||\Delta_{n}||_{1}$
$\displaystyle\geq\frac{\kappa}{3}||\Delta_{n}||_{2}^{2}-\sqrt{Gp(p+1)/2}(||\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*})||_{\infty}+{2\lambda_{n}})||\Delta_{n}||_{2}.$

When

$||\Delta_{n}||_{2}>\frac{3}{\kappa}\sqrt{Gp(p+1)/2}(||\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*})||_{\infty}+{2\lambda_{n}}),$

we have $\mathcal{F}(\Delta_{n})>0$. By the fact that
$\mathcal{F}(\Delta_{n})\leq 0$. we obtain

$||\Delta_{n}||_{2}\leq\frac{3}{\kappa}\sqrt{Gp(p+1)/2}(||\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*})||_{\infty}+{2\lambda_{n}}).$

Proof of Lemma S11. We use Lemma S10 to prove this lemma. We first check
Condition (SC1-SC2). Note that
$m_{\boldsymbol{\nu}}(\mathbf{y})=\ell(\boldsymbol{\nu},\mathbf{y})=\log(L_{M}(\boldsymbol{\nu},\mathbf{y}))+C(\mathbf{y})$
where $C(\mathbf{y})=\sum_{j=1}^{p}\log y_{j}!-2^{-1}p\log(2\pi)$ and
$L_{M}(\boldsymbol{\nu},\mathbf{y})$ is continuous at all
$\boldsymbol{\nu}\in{\mathcal{D}}$ for any fixed $\mathbf{y}$. Condition (SC1)
thus follows.

For Condition (SC2), we first check the measurability of
$\sup_{\boldsymbol{\nu}\in U}m_{\boldsymbol{\nu}}(\mathbf{y})$ for any small
ball $U$. Let $Q_{U}=\\{\boldsymbol{\nu}|\boldsymbol{\nu}\ \mbox{ is a
rational point },\boldsymbol{\nu}\in U\\}$. Then $Q_{U}$ has a countable
number of elements. From the measurability of
$m_{\boldsymbol{\nu}}(\mathbf{y})$, we get that $\sup_{\boldsymbol{\nu}\in
Q_{U}}m_{\boldsymbol{\nu}}(\mathbf{y})$ is measurable. On the other hand, from
the continuity of $m_{\boldsymbol{\nu}}(\mathbf{y})$ in $\boldsymbol{\nu}$, we
get

$\sup_{\boldsymbol{\nu}\in
Q_{U}}m_{\boldsymbol{\nu}}(\mathbf{y})=\sup_{\boldsymbol{\nu}\in
U}m_{\boldsymbol{\nu}}(\mathbf{y}),$

and hence $\sup_{\boldsymbol{\nu}\in U}m_{\boldsymbol{\nu}}(\mathbf{y})$ is
measurable. Finally, we prove ${\rm E}[\sup_{\boldsymbol{\nu}\in
U}m_{\boldsymbol{\nu}}(\mathbf{y})]<\infty$. We have

$\begin{array}[]{ll}&\quad\log(L_{M}(\boldsymbol{\nu},\mathbf{y}))\\\
&=\log\left(\sum_{g=1}^{G}\pi_{g}\int{\exp\left(\sum_{j=1}^{p}\left(x_{j}y_{j}-\exp\left(x_{j}\right)\right)-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Theta}_{g}(\mathbf{x}-\boldsymbol{\mu}_{g})\right)}\det(\boldsymbol{\Theta}_{g})^{1/2}d\mathbf{x}\right)\\\
&\leq\log\left(\sum_{g=1}^{G}\pi_{g}\int{\exp\left(\mathbf{x}^{T}\mathbf{y}-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu}_{g})^{T}\boldsymbol{\Theta}_{g}(\mathbf{x}-\boldsymbol{\mu}_{g})\right)}\det(\boldsymbol{\Theta}_{g})^{1/2}d\mathbf{x}\right)\\\
&=\log\left(\sum_{g=1}^{G}\pi_{g}\exp\left(\frac{1}{2}\mathbf{y}^{T}\boldsymbol{\Theta}_{g}^{-1}\mathbf{y}+\mathbf{y}^{T}\mu_{g}\right)\right)+2^{-1}p\log(2\pi).\end{array}$

and
$\sup_{\boldsymbol{\nu}\in{\mathcal{D}}}\frac{1}{2}\mathbf{y}^{T}\boldsymbol{\Theta}_{g}^{-1}\mathbf{y}+\mathbf{y}^{T}\boldsymbol{\mu}_{g}\leq\|\mathbf{y}\|_{2}^{2}/(2M)+M\|\mathbf{y}\|_{1}$.
Then, we have

$\sup_{\boldsymbol{\nu}\in{\mathcal{D}}}\log(L_{M}(\boldsymbol{\nu},\mathbf{y}))\leq\|\mathbf{y}\|^{2}_{2}/(2M)+M\|\mathbf{y}\|_{1}+2^{-1}p\log(2\pi).$

Also, we have

$C(\mathbf{y})=\sum_{j=1}^{p}\log
y_{j}!-2^{-1}p\log(2\pi)\leq\sum_{j=1}^{p}y_{j}\log(y_{j}+1)-2^{-1}p\log(2\pi).$

Since any polynomial of a MPLN random variable $\mathbf{y}$ is integrable, we
prove ${\rm E}(\sup_{\boldsymbol{\nu}\in
U}m_{\boldsymbol{\nu}}(\mathbf{y}))<\infty.$ In addition, we have

$\begin{array}[]{ll}M_{n}(\hat{\boldsymbol{\nu}}_{n})&\geq\lambda_{n}\sum_{g=1}^{G}\|\hat{\boldsymbol{\Theta}}_{g}\|_{1,\rm
off}+M_{n}(\boldsymbol{\nu}^{*})-\lambda_{n}\sum_{g=1}^{G}\|\boldsymbol{\Theta}_{g}^{*}\|_{1,\rm
off}\\\ &\geq
M_{n}(\boldsymbol{\nu}^{*})-\lambda_{n}\sum_{g=1}^{G}\|\boldsymbol{\Theta}_{g}^{*}\|_{1,\rm
off}\\\ &\geq M_{n}(\boldsymbol{\nu}^{*})-C\lambda_{n}\\\
&=M_{n}(\boldsymbol{\nu}^{*})-o(1).\end{array}$

Thus, all conditions in Lemma S10 are satisfied. Finally, note that
$\mathcal{D}_{0}$ only contains one element $\boldsymbol{\nu}^{*}$. Taking
$K={\mathcal{D}}$, we get, for any $\epsilon>0$,

$\mathbb{P}(||\hat{\boldsymbol{\nu}}_{n}-\boldsymbol{\nu}^{*}||_{2}\geq\epsilon)\rightarrow
0,\mbox{ as }n\rightarrow\infty.$

Proof of Lemma S12. By Theorem 16(a) in Ferguson (2017), for any
$g,i,j,g^{\prime},i^{\prime},j^{\prime}$ we only need to verify there exists a
function $F_{0}(\mathbf{y})$ such that $|{\bf
F}_{(g,i,j)(g^{\prime},i^{\prime},j^{\prime})}(\boldsymbol{\nu},\mathbf{y})|\leq
F_{0}(\mathbf{y})$ and ${\rm
E}_{\boldsymbol{\nu}^{*}}(F_{0}(\mathbf{y}))<\infty$. The existence of such
function is guaranteed by Lemma S7 .

Proof of Lemma S13. By the definition of $\mathcal{L}_{n}$ and Taylor
expansion, we have

$\displaystyle\delta\mathcal{L}_{n}$
$\displaystyle=\mathcal{L}_{n}(\boldsymbol{\nu}^{*}+\Delta_{n})-\mathcal{L}_{n}(\boldsymbol{\nu}^{*})-<\nabla\mathcal{L}_{n}(\boldsymbol{\nu}^{*}),\Delta_{n}>$
$\displaystyle=-\frac{1}{n}\sum_{i=1}^{n}\Delta_{n}^{\rm T}{\bf
F}(\check{\boldsymbol{\nu}},\mathbf{y}_{i})\Delta_{n}$
$\displaystyle=\Delta_{n}^{\rm T}({\bf
D}(\check{\boldsymbol{\nu}})-\frac{1}{n}\sum_{i=1}^{n}F(\check{\boldsymbol{\nu}},\mathbf{y}_{i})){\Delta_{n}}+{\Delta_{n}}^{\rm
T}({\bf D}(\boldsymbol{\nu}^{*})-{\bf
D}(\check{\boldsymbol{\nu}})){\Delta_{n}}+\Delta_{n}^{\rm T}(-{\bf
D}(\boldsymbol{\nu}^{*})){\Delta_{n}}$

where
$\check{\boldsymbol{\nu}}=\boldsymbol{\nu}^{*}+\theta\Delta_{n},0\leq\theta\leq
1.$ By Lemma S12, we get, for any $\epsilon$, there exists $N$ such that when
$n>N$, with probability $1-\epsilon$,

$\Delta_{n}^{\rm T}({\bf
D}(\check{\boldsymbol{\nu}})-\frac{1}{n}\sum_{i=1}^{n}{\bf
F}(\check{\boldsymbol{\nu}},\mathbf{y}_{i}))\Delta_{n}\geq-\frac{\kappa}{3}||\Delta_{n}||_{2}^{2}.$
(S10)

Also, by the continuity of ${\bf D}(\boldsymbol{\nu})$ (Remark 2) and
Wielandt-Hoffman Theorem (Bhatia, 2013), there exists a constant $\tau$ such
that when $||\Delta_{n}||_{2}<\tau$, we have

$\lambda_{min}({\bf D}(\boldsymbol{\nu}^{*})-{\bf
D}(\check{\boldsymbol{\nu}}))\geq-\frac{\kappa}{3}.$

It follows that

$\Delta_{n}^{\rm T}({\bf D}(\boldsymbol{\nu}^{*})-{\bf
D}(\check{\boldsymbol{\nu}}))\Delta_{n}\geq-\frac{\kappa}{3}||\Delta_{n}||_{2}^{2}.$
(S11)

Finally, by $\kappa=\lambda_{min}(\boldsymbol{\Gamma}^{*})$, we have

$\Delta_{n}^{\rm T}(-{\bf
D}(\boldsymbol{\nu}^{*}))\Delta_{n}\geq\kappa||\Delta_{n}||_{2}^{2}.$ (S12)

Combining the above inequalities (S10-S12), we have when
$||\Delta_{n}||_{2}<\tau$, with high probability,

$\delta\mathcal{L}_{n}\geq\frac{\kappa}{3}||\Delta_{n}||_{2}^{2}.$

Finally, since $\tau$ is a constant, by the consistency of the MPLN, we have
$\mathbb{P}(||{\Delta}_{n}|||_{2}<\tau)\rightarrow 1\mbox{ as
}n\rightarrow\infty$. Then, we get
$\delta\mathcal{L}_{n}\geq\frac{\kappa}{3}||{\Delta}_{n}|||_{2}^{2}.$ Thus we
complete the proof.

####  S1.1.6 Proof of Theorem 3

In this subsection, we simplify the notation and use $S$ and $S^{c}$ to denote
$S(\boldsymbol{\nu}^{*})$ and $S^{c}(\boldsymbol{\nu}^{*})$, respectively.
Define $\hat{\boldsymbol{\Gamma}}^{*}_{i}=-{\bf
F}(\boldsymbol{\nu}^{*},\mathbf{y}_{i})$ and
$\hat{\boldsymbol{\Gamma}}^{*}=\frac{1}{n}\sum_{i=1}^{n}\hat{\boldsymbol{\Gamma}}^{*}_{i}$,
which is an estimator of $\boldsymbol{\Gamma}^{*}$. We write
$\alpha=1-||\boldsymbol{\Gamma}^{*}_{S^{c}S}{(\boldsymbol{\Gamma}^{*}_{SS})}^{-1}||_{1,\infty}$.
By Condition (C4), $\alpha>0$.

######  Lemma S14.

For any $\epsilon>0$, with high probability,
$\hat{\boldsymbol{\Gamma}}_{SS}^{*}$ is invertible and
$||\hat{\boldsymbol{\Gamma}}^{*}_{S^{c}S}{(\hat{\boldsymbol{\Gamma}}^{*}_{SS})}^{-1}||_{1,\infty}\leq
1-\alpha/2.$

######  Lemma S15.

For any $\lambda_{n}>0$, with high probability, the solution $\hat{\nu}_{n}$
to the optimization problem (S2) is characterized by

$-n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\hat{\boldsymbol{\nu}}_{n},\mathbf{y}_{i})+\lambda_{n}{\hat{\mathbf{K}}}=0,$

where ${\hat{\mathbf{K}}}$ is the subdifferential of
$\mathcal{R}(\boldsymbol{\nu})$ at $\hat{\boldsymbol{\nu}}_{n}$.

Next, we construct the primal-dual witness solution
$(\widetilde{\boldsymbol{\nu}}_{n},{\widetilde{\mathbf{K}}})$. Let
$\widetilde{\boldsymbol{\nu}}_{n}$ be the solution to the restricted
optimization problem

$\widetilde{\boldsymbol{\nu}}_{n}\in\arg\min_{\boldsymbol{\nu}_{S^{c}}=0,\boldsymbol{\nu}\in{\mathcal{D}}}\mathcal{L}_{n}(\boldsymbol{\nu})+{\lambda_{n}}\mathcal{R}(\boldsymbol{\nu}).$

Denote

$\widetilde{\mathbf{K}}=\lambda_{n}^{-1}\left\\{n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\widetilde{\boldsymbol{\nu}}_{n},\mathbf{y}_{i})\right\\}.$
(S13)

Then, we have

$-n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\widetilde{\boldsymbol{\nu}}_{n},\mathbf{y}_{i})+\lambda_{n}{\widetilde{\mathbf{K}}}=0.$
(S14)

Define
$\widetilde{\Delta}_{n}=\boldsymbol{\nu}^{*}-\widetilde{\boldsymbol{\nu}}_{n}$.
We rewrite (S14) as

$-n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\widetilde{\boldsymbol{\nu}}_{n},\mathbf{y}_{i})+n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y}_{i})-n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y}_{i})+\hat{\boldsymbol{\Gamma}}^{*}\widetilde{\Delta}_{n}-\hat{\boldsymbol{\Gamma}}^{*}\widetilde{\Delta}_{n}+\lambda_{n}{\widetilde{\mathbf{K}}}=0.$

Let

${{\bf
R}}=-n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\widetilde{\boldsymbol{\nu}}_{n},\mathbf{y}_{i})+n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y}_{i})+\hat{\boldsymbol{\Gamma}}^{*}\widetilde{\Delta}_{n},$

and
${\mathbf{W}}=-n^{-1}\sum_{i=1}^{n}\mathcal{S}^{M}(\boldsymbol{\nu}^{*},\mathbf{y}_{i})$.
Then, we have

${{\bf
R}}-\hat{\boldsymbol{\Gamma}}^{*}\widetilde{\Delta}_{n}+{\mathbf{W}}+\lambda_{n}{\widetilde{\mathbf{K}}}=0.$
(S15)

Since we restrict the solution to the set of true support, similarly to
Theorem 1, we can proceed analogously to the proof of

$\mathbb{P}\left[||\widetilde{\nu}_{n}-\nu^{*}||_{2}\leq\frac{3}{\kappa}\sqrt{\frac{G{p(p+1)}}{{2}}}\left(n^{-1}||\nabla{\ell}_{n}(\boldsymbol{\nu}^{*})||_{\infty}+2{\lambda_{n}}{}\right)\right]\rightarrow
1,\mbox{ as }n\rightarrow\infty.$ (S16)

Considering the local convex property of loss function at
$\boldsymbol{\nu}^{*}$ in Lemma S13, if we verify that with high probability,
the strict dual feasibility condition $||\widetilde{\mathbf{K}}||_{\infty}\leq
2$ holds, then we prove that with high probability
$\widetilde{\boldsymbol{\nu}}_{n}$ is equal to $\hat{\boldsymbol{\nu}}_{n}$.
Then, the model can recover all zeros. Similarly to Lemma S15, by the fact
that $\widetilde{\boldsymbol{\nu}}_{n}$ is the restricted construction and the
definition of $\widetilde{\mathbf{K}}$ in (S13), we have

$||\widetilde{\mathbf{K}}_{S}||_{\infty}\leq 2.$ (S17)

Therefore, we only need to show that
$||\widetilde{\mathbf{K}}_{S^{c}}||_{\infty}\leq 2$. Note that the infinity
norm of $\widetilde{\mathbf{K}}_{S^{c}}$ is less than 2 instead of 1. The
reason is
$\mathcal{R}(\boldsymbol{\nu})=\sum_{g=1}^{G}||\boldsymbol{\Theta}_{g}||_{1,\rm
off}=2\sum_{g=1}^{G}\sum_{i<j}|\boldsymbol{\Theta}_{g,ij}|$. Hence, we aim to
verify the strict dual feasibility.

######  Lemma S16 (Strict dual feasibility).

Under Condition (C1-C4), suppose that $\hat{\boldsymbol{\Gamma}}_{SS}^{*}$ is
invertible and

$||{\mathbf{W}}||_{\infty}+||{{\bf
R}}||_{\infty}<\frac{\alpha\lambda_{n}}{4},\,\left\|\hat{\boldsymbol{\Gamma}}^{*}_{S^{c}S}{\left(\hat{\boldsymbol{\Gamma}}^{*}_{SS}\right)}^{-1}\right\|_{1,\infty}\leq
1-\alpha/2.$

Then, the matrix $\widetilde{\mathbf{K}}$ satisfies

$||\widetilde{\mathbf{K}}_{S^{c}}||_{\infty}\leq 2.$

######  Lemma S17.

$\frac{||{{\bf R}}||_{\infty}}{||\widetilde{\Delta}_{n}||_{2}}\rightarrow 0$,
in probability as $n\rightarrow\infty$.

Applying Chebyshev’s inequality, the following lemma is clear.

######  Lemma S18.

Let $a_{n}$ be any sequence such that
${a_{n}}{\sqrt{n}}\rightarrow\infty,a_{n}>0$ and $a_{n}\rightarrow 0$. Then,
we have

$\mathbb{P}(||{\mathbf{W}}||_{\infty}\leq a_{n})\rightarrow 1,\mbox{ as
}n\rightarrow\infty.$

Proof of Theorem 3. Note that
${\mathbf{W}}=-n^{-1}\nabla\ell_{n}(\boldsymbol{\nu}^{*})$ and
$\widetilde{\Delta}_{n}=\boldsymbol{\nu}^{*}-\widetilde{\boldsymbol{\nu}}_{n}$.
Applying Lemma S17, we have $||{{\bf
R}}||_{\infty}=o_{p}(1)||\widetilde{\Delta}_{n}||_{2}.$ By (S16), we have

$||{{\bf
R}}||_{\infty}=o_{p}(1)\left(\frac{3}{\kappa}\sqrt{\frac{G{p(p+1)}}{{2}}}\left(||{\mathbf{W}}||_{\infty}+2{\lambda_{n}}\right)\right).$

It follows that

$||{\mathbf{W}}||_{\infty}+||{{\bf
R}}||_{\infty}=||{\mathbf{W}}||_{\infty}+o_{p}(1)\left(\frac{3}{\kappa}\sqrt{\frac{G{p(p+1)}}{{2}}}\left(||{\mathbf{W}}||_{\infty}+2{\lambda_{n}}\right)\right)=(1+o_{p}(1))||{\mathbf{W}}||_{\infty}+o_{p}(1)\lambda_{n}.$

In order to show $||{\mathbf{W}}||_{\infty}+||{{\bf
R}}||_{\infty}<\frac{\alpha\lambda_{n}}{4},$ we only need to show that with
high probability,

$||{\mathbf{W}}||_{\infty}\leq(1+o_{p}(1))^{-1}\left(\frac{\alpha}{4}-o_{p}(1)\right)\lambda_{n}.$

By the choice of $\lambda_{n}$ such that $\lambda_{n}\rightarrow 0$ and
${\sqrt{n}\lambda_{n}}\rightarrow\infty$, applying Lemma S18, we have
$\mathbb{P}(||{\mathbf{W}}||_{\infty}\leq\left(1+o_{p}(1))^{-1}\left(\frac{\alpha}{4}-o_{p}(1)\right)\lambda_{n}\right)\rightarrow
1.$ By Lemma S14, we have with high probability
$\hat{\boldsymbol{\Gamma}}_{SS}^{*}$ is invertible and
$\left\|\hat{\boldsymbol{\Gamma}}^{*}_{S^{c}S}{\left(\hat{\boldsymbol{\Gamma}}^{*}_{SS}\right)}^{-1}\right\|_{1,\infty}\leq
1-\alpha/2.$ Then, by Lemma S16, with high probability, we have the strict
dual feasibility condition holds. Hence, the witness solution
$\widetilde{\boldsymbol{\nu}}_{n}$ is equal to the original solution
$\hat{\boldsymbol{\nu}}_{n}$. Since $\widetilde{\boldsymbol{\nu}}_{n}$ is the
restricted solution, with high probability, $\widetilde{\boldsymbol{\nu}}_{n}$
can recover all zeros. It follows that $\hat{\boldsymbol{\nu}}_{n}$ can
recover all zeros. Finally, since
$||\hat{\boldsymbol{\nu}}_{n}-\boldsymbol{\nu}^{*}||_{2}\rightarrow 0$ in
probability, with high probability, $\hat{\boldsymbol{\nu}}_{n}$ can recover
all non zeros.

Proof of Lemma S15. Since $\boldsymbol{\nu}^{*}$ is an interior point of
$\mathcal{D}$ , we only need to prove with high probability, the maximum value
will not be taken at $\partial{\mathcal{D}}$. By Jessen’s inequality and the
identifiablity of the MPLN, we have
$\sup_{\boldsymbol{\nu}\in\partial{\mathcal{D}}}{\rm
E}_{\boldsymbol{\nu}^{*}}\left[\log\frac{p\left(\mathbf{y}_{i};\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)}{p\left(\mathbf{y}_{i};\boldsymbol{\nu}^{*},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)}\right]<0$.
By uniform law of large numbers, we have with high probability,

$\sup_{\boldsymbol{\nu}\in\partial{\mathcal{D}}}\left\\{n^{-1}\sum_{i=1}^{n}\log\,p\left(\mathbf{y}_{i};\boldsymbol{\nu},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)-n^{-1}\sum_{i=1}^{n}\log\,p\left(\mathbf{y}_{i};\boldsymbol{\nu}^{*},\\{\boldsymbol{\mu}_{g}\\}_{g=1}^{G}\right)\right\\}<0.$

Thus we complete the proof.

Proof of Lemma S16. We split $\widetilde{\Delta}_{n}$ into
$\widetilde{\Delta}_{nS}$ and $\widetilde{\Delta}_{nS^{c}}$. By
$\widetilde{\Delta}_{nS^{c}}=0$, (S15) can be rewritten as two blocks of
linear equations

${{\bf
R}}_{S}-\hat{\boldsymbol{\Gamma}}_{SS}^{*}\widetilde{\Delta}_{nS}+{\mathbf{W}}_{S}+\lambda_{n}{\widetilde{\mathbf{K}}}_{S}=0,$
(S18) ${{\bf
R}}_{S^{c}}-\hat{\boldsymbol{\Gamma}}_{S^{c}S}^{*}\widetilde{\Delta}_{nS}+{\mathbf{W}}_{S^{c}}+\lambda_{n}{\widetilde{\mathbf{K}}}_{S^{c}}=0.$
(S19)

From (S18), we have

$\widetilde{\Delta}_{nS}={\left(\hat{\boldsymbol{\Gamma}}^{*}_{SS}\right)}^{-1}\left({{\bf
R}}_{S}+{\mathbf{W}}_{S}+\lambda_{n}{\widetilde{\mathbf{K}}}_{S}\right).$

Substituting this expression into (S19), with high probability, we have

${\widetilde{\mathbf{K}}}_{S^{c}}=\lambda_{n}^{-1}\left\\{-{{\bf
R}}_{S^{c}}+\hat{\boldsymbol{\Gamma}}_{S^{c}S}^{*}{\left(\hat{\boldsymbol{\Gamma}}^{*}_{SS}\right)}^{-1}\left({{\bf
R}}_{S}+{\mathbf{W}}_{S}+\lambda_{n}{\widetilde{\mathbf{K}}}_{S}\right)-{\mathbf{W}}_{S^{c}}\right\\}.$

Let $\mathbf{A}$ be a matrix and $\mathbf{a}$ be a vector. By the fact that
$||\mathbf{A}\mathbf{a}||_{\infty}\leq||\mathbf{A}||_{1,\infty}||\mathbf{a}||_{\infty}$,
we have

$\displaystyle||{\widetilde{\mathbf{K}}}_{S^{c}}||_{\infty}$
$\displaystyle\leq{\lambda_{n}}^{-1}\left(1+\bigg{|}\bigg{|}\hat{\boldsymbol{\Gamma}}_{S^{c}S}^{*}{\left(\hat{\boldsymbol{\Gamma}}^{*}_{SS}\right)}^{-1}\bigg{|}\bigg{|}_{1,\infty}\right)\left(||{\mathbf{W}}||_{\infty}+||{{\bf
R}}||_{\infty}\right)+\bigg{|}\bigg{|}\hat{\boldsymbol{\Gamma}}_{S^{c}S}^{*}{\left(\hat{\boldsymbol{\Gamma}}^{*}_{SS}\right)}^{-1}\bigg{|}\bigg{|}_{1,\infty}||{\widetilde{\mathbf{K}}}_{S}||_{\infty}$
$\displaystyle\leq
2(1-\alpha/2)+\frac{(2-\alpha/2)}{\lambda_{n}}(||{\mathbf{W}}||_{\infty}+||{{\bf
R}}||_{\infty})$ $\displaystyle\leq 2.$

Here we use the inequality that $||{\widetilde{\mathbf{K}}}_{S}||_{\infty}\leq
2$ in (S17). Thus we complete the proof.