# Statistical analysis for a penalized EM algorithm in high-dimensional
mixture linear regression model

Ning Wang, Xin Zhang and Qing Mai Ning Wang<EMAIL_ADDRESS>is
Assistant Professor, Center of Statistics and Data Science, Beijing Normal
University, Zhuhai, 519807, China; Xin Zhang<EMAIL_ADDRESS>is Associate
Professor, Department of Statistics, Florida State University, Tallahassee,
32312, Florida, USA; Qing Mai<EMAIL_ADDRESS>is Associate Professor,
Department of Statistics, Florida State University, Tallahassee, 32312,
Florida, USA.

(Beijing Normal University and Florida State University)

###### Abstract

The expectation-maximization (EM) algorithm and its variants are widely used
in statistics. In high-dimensional mixture linear regression, the model is
assumed to be a finite mixture of linear regression and the number of
predictors is much larger than the sample size. The standard EM algorithm,
which attempts to find the maximum likelihood estimator, becomes infeasible
for such model. We devise a group lasso penalized EM algorithm and study its
statistical properties. Existing theoretical results of regularized EM
algorithms often rely on dividing the sample into many independent batches and
employing a fresh batch of sample in each iteration of the algorithm. Our
algorithm and theoretical analysis do not require sample-splitting, and can be
extended to multivariate response cases. The proposed methods also have
encouraging performances in numerical studies.

Key Words: EM algorithm, High-dimensional regression, Mixture model.

##  1 Introduction

Consider a univariate response $Y\in\mathbb{R}$ and a $p$-dimensional
predictor $X\in\mathbb{R}^{p}$. The mixture linear regression model assumes
that

$Y=\beta_{k}^{T}X+\epsilon,\ \mathrm{for}\ k=1,\cdots,K,\ \mathrm{with}\
\mathrm{probability}\ \omega_{k}>0,$ (1)

where $K\geq 2$ is the number of mixtures, $\sum_{k=1}^{K}\omega_{k}=1$,
$\epsilon\sim N(0,\sigma^{2})$, $\sigma^{2}>0$, is independent of $X$, and
$\beta_{k}$ is the $p$-dimensional regression coefficient vector that
characterizes the linear relationship between $Y$ and $X$ in the $k$-th
mixture. By introducing a latent variable $W\in\\{1,\cdots,K\\}$, independent
of $X$, model (1) is equivalent to

$\mathbb{P}(W=k)=\omega_{k},\quad Y\mid(X,W=k)\sim
N(\beta_{k}^{T}X,\sigma^{2}).$ (2)

We consider the high-dimensional joint estimation of all the ${\beta}_{k}$’s
and provide a general estimation procedure with strong theoretical guarantees.
The latent mixtures, indicated by the latent variable $W$ in (2), make the
estimation problem much more challenging than the linear regression,
especially in high dimensions. We assume that only a subset of predictors,
indexed by $S\subset\\{1,\cdots,p\\}$, is relevant to the regression.
Therefore, $(\beta_{k})_{\mathcal{S}^{c}}=0$ for all $k$, where $S^{c}$ is the
complement of $S$. The group lasso penalty (Yuan and Lin, 2006) is naturally
applied to those $p$-dimensional $\beta_{k}$ vectors in the maximization steps
of our regularized EM algorithm to select relevant variables across all the
mixtures.

The mixture linear model and finite mixture models in general, are widely used
to account for heterogeneity in data analysis (e.g., Turner, 2000; McLachlan
and Peel, 2004; McLachlan et al., 2019). When the number of predictors is not
large, the latent mixtures and the model parameters can be estimated using the
expectation-maximization (EM, Dempster et al., 1977) algorithm. The EM
algorithm is the dominant solution for finding the maximum likelihood
estimators of mixture regression models like (1). When $\epsilon$ is not
normally distributed, the EM algorithm has been modified and extended to
robust fitting of mixture linear models under the t- and the Laplace
distributions in Yao et al. (2014) and Song et al. (2014), respectively.
Leisch (2004) provides computational and implementation details of mixture
regression models.

Although the EM algorithm has been extensively used in mixture regression
models, it is challenging to establish a rigorous theoretical characterization
of the finite-sample estimates in the iterative algorithm. Some groundbreaking
progress has been made in recent years. Balakrishnan et al. (2017) laid
theoretical foundations for quantifying the EM updates’ convergence within
statistical precision of a global optimum. For mixture linear model, strong
theoretical guarantees of the EM algorithm are often established by focusing
on the model with two equal mixtures, $K=2$ and $\omega_{1}=\omega_{2}=1/2$,
and symmetric regression coefficient vectors, $\beta_{2}=-\beta_{1}$, (e.g.,
Kwon et al., 2021, and references therein). Then the EM algorithm is
simplified substantially because $\omega_{1}=\omega_{2}=1/2$ does not require
estimation and, more importantly, the model parameters are reduced to a single
vector $\beta\equiv\beta_{1}=-\beta_{2}\in\mathbb{R}^{p}$. Many theoretical
results are established for the “sample-splitting” EM algorithm, which divides
the full data into $T$ equal batches and uses a new batch of samples in each
iteration. Without sample splitting, theoretical analysis becomes much more
challenging. This is because the function to be maximized in each EM
iteration, namely the Q-function, involves both the random samples and the
current parameter estimates, which are made independent by sample-splitting.
See Balakrishnan et al. (2017) and Klusowski et al. (2019) for recent studies
of the sample-splitting EM algorithm in mixture linear models. While the
above-mentioned works all focus on low-dimensional models and unpenalized EM
algorithm, regularized EM algorithm for high-dimensional mixture linear models
is of growing interest in recent years. A selective review is as follows.

Khalili and Chen (2007) studied the variable selection for mixture linear
regression with penalized likelihood. When $p$ is fixed and $n$ goes to
infinity, they established variable selection consistency and root-$n$
consistency for a possible local maximum. Städler et al. (2010) proposed a
lasso-type penalty (Tibshirani, 1996) on the negative-log-likelihood function
and obtained non-asymptotic convergence results for the global optimum.
However, there is no guarantee for the EM algorithm to attain either the
particular local maximum or the global maximum. To rigorously study the entire
EM iterative solution sequence in high dimensions, existing results rely on
the sample-splitting procedure. For example, convergence results based on
sample-splitting algorithms are established for a truncated EM algorithm (Wang
et al., 2014), a penalized EM algorithm (Yi and Caramanis, 2015), and a
stochastic EM algorithm (Zhu et al., 2017), all under the assumptions of
$K=2$, $\omega_{1}=\omega_{2}=1/2$, $\beta_{2}=-\beta_{1}$, and normally
distributed predictors. More recently, with the help of sample splitting,
Zhang et al. (2020) systematically studied estimation, confidence intervals,
and large-scale hypotheses testing for the mixture linear model. Sample
splitting is undoubtedly used to facilitate theoretical analyses, but is not
desirable in practice. It remains unknown how to practically choose $T$, which
is both the number of iterations and the number of sample batches, and how it
affects the estimation. Moreover, it is believed (e.g., Zhang et al., 2020)
that data splitting is unnecessary in numerical studies.

The contributions of this article are multi-fold. The first and most
significant contribution is developing a practical penalized EM algorithm for
a general high-dimensional mixture linear model and establishing its
substantial theoretical guarantees. Our algorithm is equally applicable to
random-$X$ and fixed-$X$ and multiple mixtures $K\geq 2$. In our theoretical
analysis, we do not require sample splitting and allow a relatively general
model. Specifically, we establish a non-asymptotic convergence rate for a two-
mixture linear model with unknown proportions $\omega_{1},\omega_{2}\in(0,1)$,
unrelated two regression parameter vectors
$\beta_{1},\beta_{2}\in\mathbb{R}^{p}$, and normally distributed predictors
with unknown covariance structure. To our best knowledge, we established the
first theoretical results for the high-dimensional mixture linear regression
under such a general setting without sample splitting. Compared with the high-
dimensional linear model literature, the theoretical analysis for the high-
dimensional mixture model requires bounding the supremum of random processes
and is much more challenging without sample splitting. Our general proof
strategy is related to Cai et al. (2019), which studies the penalized EM
algorithm for the Gaussian mixture model without data splitting. The
complicated relationship between the random response $Y$ and the random
predictor $X$ makes the theoretical studies of mixture linear regression even
more challenging than the Gaussian mixture model, which only involves random
$X$. Many new concentration results are needed for random processes involving
both $Y$ and $X$. For instance, unlike $X_{i}$ in the Gaussian mixture model
that is sub-Gaussian, the product term $X_{i}Y_{i}$ appears frequently in the
EM iterates and estimates and is more difficult to bound (Adamczak, 2008).
Even for the theoretical analysis of the population EM iterates, double
expectations $\mathrm{E}_{X}\\{\mathrm{E}_{Y\mid X}(\cdot)\\}$ are needed than
single expectation $\mathrm{E}_{X}(\cdot)$ in the Gaussian mixture model. With
substantial efforts, we obtain a near optimal convergence rate of
$\mathrm{log}n\sqrt{s\mathrm{log}p/n}$, with a small price $\mathrm{log}(n)$
to pay for not sample splitting.

The second contribution is our new theoretical insights on model
misspecification. Specifically, we analyze how a fixed parameter value
$\sigma^{2}$ in the penalized EM algorithm may affect the estimation of
$\beta$. To the best of our knowledge, it has not been studied in the
literature on the mixture linear regression model. In most theoretical studies
considering the EM algorithm for the mixture linear regression, the variance
$\sigma^{2}$ is often assumed to be the true parameter value
$(\sigma^{\star})^{2}$ and is a fixed constant. Kwon et al. (2021) considered
the convergence results when $\sigma^{2}$ is updated in each M step. However,
they only considered a simplified case where the mixture proportions are known
to be $1/2$ and are not updated in each M step. Besides, their theoretical
studies are limited to low dimensions. In general, when the error variance
$\sigma^{2}$ is not correctly specified in the penalized EM algorithm brings
bias for estimation. However, for the mixture model with a relatively large
signal-to-noise ratio, our theory indicates that the choice of $\sigma^{2}$
has a minor influence on the estimation, and thus an accurate estimation for
it is usually not necessary. This conclusion is further demonstrated by a
simulation study.

The third contribution is that we extend the study for the mixture linear
regression model to multiple response cases. For the mixture linear regression
model with a multivariate response, the naive approach is fitting a mixture
linear regression model separately for each element of the response. The major
drawback of doing so is each observation may be identified into different
clusters when we model each univariate response separately. We illustrated the
advantages of considering multiple responses together than handling them
separately from both theoretical and numerical aspects. The advantage of
considering multiple responses together is also demonstrated in Hyun et al.
(2020), which developed a sparse mixture linear regression model to estimate
the time-varying data sets (i.e. at each time point, the data satisfy a
mixture linear regression model). However, they are interested in developing a
time-varying model to analyze a real-world dataset, but no theoretical study
is conducted. In contrast, we rigorously characterize the advantage of
considering multiple responses simultaneously with statistical theory.

The rest of the article is organized as follows. Section 2 contains the
implementation details and discussions about the penalized EM algorithm.
Section 3 presents the theory for the penalized EM algorithm and the influence
of the choice of $\sigma^{2}$. Simulation studies are presented in Section 4.
We then extend the mixture linear regression to multiple response cases and
consider its theoretical studies in Section 5. In Section 6, we consider a
real data example followed by a short discussion in Section 7. The
supplementary material contains proofs for all the Lemmas and Theorems and
additional implementation details.

##  2 Estimation

We assume that we collect $n$ independent data points
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ from model (1). In this section, we do not
make additional assumptions on $X$ as our estimation procedure is equally
applicable to random or fixed $X$ and allows for both continuous and discrete
predictors. Let
${\theta}=\\{\omega_{1},\cdots,\omega_{k},{\beta}_{1},\cdots,{\beta}_{k}\\}$
be the unknown parameters to be estimated. In this section, we focus on the
studies of the regression coefficients and treat $\sigma^{2}$ as known. The
estimation for $\sigma^{2}$ is briefly discussed at the end of this section.
We further show in Theorem 2 that misspecification of $\sigma^{2}$ has a
relatively small impact on the final estimation.

To motivate our proposal, we first derive the standard EM algorithm and
discuss its limitations. The EM algorithm aims to maximize the log-likelihood
of $Y\mid X$ over $\theta$, by iteratively updating the sequence of solutions
$\\{\widehat{\theta}^{(t)},t=0,1,\ldots\\}$ via the Expectation-step (E-step)
and the Maximization-step (M-step). Recall that $W$ is the latent variable
representing the mixtures. Consider the $(t+1)$-th iteration with the current
value $\widehat{\theta}^{(t)}$. In the E-step, we calculate the expectation of
the log-likelihood of $W\mid(Y,X)$ at the parameter $\widehat{\theta}^{(t)}$
and the observed $(X,Y)$. This is known as the Q-function,

$\displaystyle Q(\theta\mid\widehat{\theta}^{(t)})$
$\displaystyle=-\frac{1}{2n}\sum_{i=1}^{n}\sum_{k=1}^{K}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})(Y_{i}-X_{i}^{T}\beta_{k})^{2}+\frac{1}{n}\sum_{i=1}^{n}\sum_{k=1}^{K}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})\mathrm{log}(\omega_{k}),$
(3)

where $\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})=\mathbb{P}(W_{i}=k\mid
Y_{i},X_{i},\widehat{\theta}^{(t)})$. The estimated probability
$\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})$ is given by

$\displaystyle\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})$
$\displaystyle=\frac{\widehat{\omega}_{k}^{(t)}\phi_{\sigma^{2}}(Y_{i}-X_{i}^{T}\widehat{{\beta}}_{k}^{(t)})}{\sum_{k=1}^{K}\widehat{\omega}_{k}^{(t)}\phi_{\sigma^{2}}(Y_{i}-X_{i}^{T}\widehat{{\beta}}_{k}^{(t)})},$
(4)

where $\phi_{\sigma^{2}}(u)$ is the probability density function of
$N(0,\sigma^{2})$. Then, in the M-step, we update
$\widehat{\theta}_{k}^{(t+1)}=\operatorname*{argmax}_{\theta}Q(\theta\mid\widehat{\theta}^{(t)})$
by maximizing (3).

Note that the standard EM algorithm is infeasible to high-dimensional
problems. If $p\gg n$, even when the latent random variables $W_{i}$’s are
observed, the maximizer of the Q function is not well-defined. Moreover, as
$W_{i}$ is generally latent and unobserved, we need to calculate
$\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})$ in (4), which involves the
$p$-dimensional random vector $X$ and the $p$-dimensional parameter vectors
$\beta_{k}$’s.

To estimate the linear mixture regression model in high dimensions, we modify
the standard EM algorithm by encouraging sparsity. In high-dimensional
statistics, it is often assumed that the coefficients have many elements as
zero, i.e, most elements in ${\beta}_{k}$ are zero. But we further assume that
${\beta}_{k}$ has a joint sparsity structure, in that, for most $j$, we have
$\beta_{1j}=\ldots=\beta_{Kj}=0$. The group sparsity facilitates the
interpretation, as for each $j$, if $\beta_{1j}=\ldots=\beta_{Kj}=0$, then
$X_{j}$ is unimportant for the prediction of $Y$ regardless of which mixture
the observation comes from. Moreover, the group sparsity benefits the E-step,
because, with straightforward calculation, we can rewrite (4) as

$\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})={\widehat{\omega}_{k}^{(t)}}/\big{(}\widehat{\omega}_{k}^{(t)}+\sum_{k^{\prime}\neq
k}\widehat{\omega}_{k^{\prime}}^{(t)}\mathrm{exp}\big{\\{}(\widehat{\beta}_{k^{\prime}}^{(t)}-\widehat{\beta}_{k}^{(t)})^{T}X_{i}(Y_{i}-{(\widehat{\beta}_{k}^{(t)}+\widehat{\beta}_{k^{\prime}}^{(t)})^{T}X_{i}}/{2})/\sigma^{2}\big{\\}}\big{)}.$
(5)

Equation (5) shows an advantage of the group sparsity over the individual
sparsity. It implies $X_{j}$ is unimportant for the evaluation of
$\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})$ if
$\widehat{\beta}_{k^{\prime}}^{(t)}-\widehat{\beta}_{k}^{(t)}=0$ for all
$k,k^{\prime}$. The group sparsity guarantees that such a situation happens
for most $X_{j}$ and $\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})$ is
determined by a few elements in ${X}$.

With the sparsity assumption, we modify the EM algorithm by imposing the group
lasso penalty (Yuan and Lin, 2006) on $\beta_{k}$, for $k=1,\cdots,K$. Our the
penalized EM algorithm replaces M-step by

$\widehat{\theta}_{k}^{(t+1)}=\operatorname*{argmax}_{\theta}\left\\{2Q(\theta\mid\widehat{\theta}^{(t)})-\lambda_{n}^{(t+1)}\sum_{j=1}^{p}\sqrt{\sum_{k=1}^{K}\beta_{kj}^{2}}\right\\},$
(6)

where $\lambda_{n}^{(t+1)}>0$ is the tuning parameter at the $(t+1)$-th
iteration and $\beta_{kj}$ be the $j$-th element of $\beta_{k}$. In the
E-step, we evaluate $\widehat{\theta}_{k}^{(t+1)}$ by (5), which is the same
as in the standard EM algorithm. Clearly, our penalized EM algorithm reduces
to the standard EM algorithm when $\lambda_{n}^{(t+1)}=0$ for all
$t=0,1,\ldots$.

The optimization in (6) is separable in $\omega_{k}$ and ${\beta}_{k}$. For
$\omega$, the updating equation is the same as in the standard EM algorithm,
i.e,
$\widehat{\omega}_{k}^{(t+1)}=\sum_{i=1}^{n}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})/n$.
For $\beta_{k}$, it amounts to minimizing the following objective function,

$\displaystyle\ell(\beta_{1},\cdots,\beta_{K})$
$\displaystyle=\sum_{k=1}^{K}\beta_{k}^{T}\widehat{\Sigma}_{k}^{(t+1)}\beta_{k}-2\sum_{k=1}^{K}(\widehat{\rho}_{k}^{(t+1)})^{T}\beta_{k}+\lambda_{n}^{(t+1)}\sum_{j=1}^{p}\sqrt{\sum_{k=1}^{K}\beta_{kj}^{2}},$
(7)

where
$\widehat{\rho}_{k}^{(t+1)}=\sum_{i=1}^{n}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})X_{i}Y_{i}/n$
and
$\widehat{\Sigma}_{k}^{(t+1)}=\sum_{i=1}^{n}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})X_{i}X_{i}^{T}/n$.
The convex optimization in (7) can be done efficiently by the groupwise
majorization descent algorithm (Yang and Zou, 2015). We provide implementation
details in Supplement Section S.1.

The tuning parameter $\lambda_{n}^{(t)}$ could either be fixed or varying
across iterations. For theoretical consideration, in Algorithm 1, we set
$\lambda_{n}^{(t+1)}=\kappa\lambda_{n}^{(t)}+C_{\lambda}\sqrt{\mathrm{log}(p)\mathrm{log}(n)^{2}/n}$,
where $0<\kappa<1/2$ and $C_{\lambda}$ are generic constants, for ease of
showing the statistical convergence results. Note that $\lambda_{n}^{(t)}$ is
at the order of $\sqrt{\mathrm{log}(p)\mathrm{log}(n)^{2}/n}$ when $t$ is
large. Thus, in practice, we fix $\lambda_{n}^{(t)}=\lambda$ for all $t$ and
tune $\lambda$ by the Bayesian information criterion (see our numerical
studies). For fixed $\lambda_{n}^{(t)}=\lambda$ over all iterations, our
penalized EM algorithm is maximizing

$L(\theta)-\lambda/2\sum_{j=1}^{p}\sqrt{\sum_{k=1}^{K}\beta_{kj}^{2}},$ (8)

where $L(\theta)$ is the conditional log-likelihood of $Y\mid X$. The
following lemma shows the convergence result of Algorithm 1.

######  Lemma 1.

If we set $\lambda_{n}^{(t)}=\lambda$ for all $t$ in Algorithm 1, the
objective function from (8) evaluated at $\theta^{(t+1)}$ is guaranteed to be
no less than the objective function from (8) evaluated at $\theta^{(t)}$. That
is, the sequence of iterates $\\{\theta^{(t)}\\}_{t=1}^{\infty}$ generated by
Algorithm 1 monotonically increase the value of the objective function from
(8).

In Algorithm 1, $\sigma^{2}$ is treated as a known parameter to facilitate
theoretical studies. Treating $\sigma^{2}$ as known is also common in
theoretical studies for mixture linear regression (Yi and Caramanis, 2015;
Balakrishnan et al., 2017; Zhang et al., 2020, e.g.). We leave the detailed
discussion for it in Section 3. In practice, the estimates of $\sigma^{2}$ can
be updated straightforwardly in the EM algorithm as
$n^{-1}\sum_{i=1}^{n}\sum_{k=1}^{K}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t+1)})(Y_{i}-X_{i}^{T}\widehat{\beta}_{k}^{(t+1)})^{2}$.
Similar estimates are also adopted in numerical studies of Zhang et al.
(2020).

Regularization strategies are also used by Yi and Caramanis (2015), Cai et al.
(2019), and Zhang et al. (2020) in high-dimensional EM algorithms. However,
Cai et al. (2019) considers the classification problem instead of regression
problem. Yi and Caramanis (2015); Zhang et al. (2020) studies the regression
problem with the addition of the lasso penalty (Tibshirani, 1996) instead of
the group lasso penalty. However, a key difference between our algorithm and
theirs is that they require data to be split into $T$ batches, and the
(penalized) EM algorithm iterates $T$ times, using one independent batch at
each iteration. The sample splitting is rarely performed in standard EM
algorithms on low-dimensional data, as it may decrease the computation
efficiency. Instead, the sample splitting is an attempt to circumvent
technical difficulty in proving the convergence rate. Our proposed EM
algorithm does not split sample, but we will show that it achieves a high
level of accuracy regardless. Moreover, in addition to investigating the
property of the penalized EM algorithm in estimating model (1), we also study
the effect of misspecification of $\sigma^{2}$ and the estimation of mixture
linear regression when there are multiple responses; see Theorem 2 and Section
5, respectively.

Algorithm 1 Group lasso penalized EM algorithm for model (1) in high
dimensions

  1. Input:

Initial values $\widehat{\omega}_{k}^{(0)}$, $\widehat{\beta}_{k}^{(0)}$, for
$k=1,\cdots,K$, maximum iteration number $T$, data
$\\{X_{i},Y_{i};i=1,\ldots,n\\}$, and initial tuning parameter

$\displaystyle\lambda_{n}^{(0)}=C_{1}(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{\beta}_{1}^{(0)}-\beta_{1}^{*}\|_{2}\lor\cdots\lor\quad\|\widehat{\beta}_{K}^{(0)}-\beta_{K}^{*}\|_{2})/{\sqrt{s}}+C_{\lambda}\sqrt{\mathrm{log}(n)^{2}\mathrm{log}(p)/s}$

for some positive constants $C_{1}$ and $C_{\lambda}$.

  2. Iterate:

For $t=0,\cdots,T-1$, do the following steps until convergence.

     * •

For $i=1,\cdots,n$, let

$\displaystyle\quad\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})={\widehat{\omega}_{k}^{(t)}}/\Big{(}\widehat{\omega}_{k}^{(t)}+\sum_{k^{\prime}\neq
k}\widehat{\omega}_{k^{\prime}}^{(t)}\mathrm{exp}\big{\\{}(\widehat{\beta}_{k^{\prime}}^{(t)}-\widehat{\beta}_{k}^{(t)})^{T}\cdot
X_{i}\big{(}Y_{i}-{(\widehat{\beta}_{k}^{(t)}+\widehat{\beta}_{k^{\prime}}^{(t)})^{T}X_{i}}/{2}\big{)}/\sigma^{2}\big{\\}}\Big{)}.$

     * •

For $k=1,\cdots,K$, update

$\displaystyle\widehat{\omega}_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)}),$
$\displaystyle\widehat{\rho}_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\big{(}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})X_{i}Y_{i}\big{)},$
$\displaystyle\widehat{\Sigma}_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\big{(}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})X_{i}X_{i}^{T}\big{)},$

and update $\widehat{\beta}_{k}^{(t+1)}$ by minimizing

$\displaystyle\sum_{k=1}^{K}\beta_{k}^{T}\widehat{\Sigma}_{k}^{(t+1)}\beta_{k}-2\sum_{k=1}^{K}(\widehat{\rho}_{k}^{(t+1)})^{T}\beta_{k}+\lambda_{n}^{(t+1)}\sum_{j=1}^{p}\sqrt{\sum_{k=1}^{K}\beta_{kj}^{2}}.$

with
$\lambda_{n}^{(t+1)}=\kappa\lambda_{n}^{(t)}+C_{\lambda}\sqrt{\mathrm{log}(p)\mathrm{log}(n)^{2}/n}$,
where $\kappa\in(0,1/2)$.

  3. Output:

$\widehat{\theta}^{(t+1)}=\\{\widehat{\omega}_{1}^{(t+1)},\ldots,\widehat{\omega}_{K}^{(t+1)},\widehat{\beta}_{1}^{(t+1)},\ldots,\widehat{\beta}_{K}^{(t+1)}\\}$

##  3 Theory

###  3.1 Preliminary

We begin this section with some notations. For numbers $a$ and $b$, $a\lor b$
means $\max\\{a,b\\}$. For an integer $n$, we let $[n]$ denote the set
$\\{1,\cdots,n\\}$. For a vector $x=(x_{1},\cdots,x_{p})^{T}$, $\|x\|_{0}$ is
the number of non-zero elements in $x$, $\|x\|_{1}=\sum_{i=1}^{p}|x_{i}|$, and
$\|x\|_{2}=\sqrt{\sum_{i=1}^{p}x_{i}^{2}}$. For a symmetric matrix $A$, we
denote $\lambda_{min}(A)$ and $\lambda_{max}(A)$ as the smallest and largest
eigenvalues of $A$, respectively. The Frobenius norm of a matrix $A=(a_{ij})$
is defined as $\|A\|_{F}=\sqrt{\sum_{i,j}a_{ij}^{2}}$. The $\ell_{2}$ norm of
a matrix $A$ is $\|A\|_{2}=\sqrt{\lambda_{max}(A^{T}A)}$. For a subset
$\mathcal{A}\subseteq\\{1,\cdots,p\\}$, $\mathcal{A}^{c}$ denotes its
complement. For two sequences of positive numbers $a_{n}$ and $b_{n}$,
$a_{n}=O(b_{n})$ means $a_{n}\leq cb_{n}$ for a constant $c>0$ for all $n$,
$a_{n}=o(b_{n})$ means that $a_{n}/b_{n}\to 0$ as $n\to\infty$, and $b_{n}\gg
a_{n}$ means that $a_{n}=o(b_{n})$. Let $\mathcal{S}^{p-1}$ be the unit
sphere. For a positive integer $s\leq p/2$, let set
$\Gamma_{2p}(2s)=\\{\mu\in\mathbb{R}^{2p}:\|\mu_{S^{c}}\|_{1}\leq
5\sqrt{2s}\|\mu_{S}\|_{2}+2\sqrt{2s}\|\mu\|_{2}$ for some $S\subset[2p]$ with
$|S|=4s\\}$ and $\Gamma(s)=\Gamma_{2p}(2s)_{1:p}$, where
$\Gamma_{2p}(2s)_{1:p}=\\{\mu_{1:p}:\mu\in\Gamma_{2p}(2s)\\}$. For a vector
$x$ and a symmetric matrix $A$, we define
$\|x\|_{2,s}=\sup_{\|\mu\|_{2}=1,\mu\in\Gamma(s)}\langle x,\mu\rangle$, and
$\|A\|_{2,s}=\sup_{\|\mu\|_{2}=1,\mu\in\Gamma(s)}|\mu^{T}A\mu|$.

We assume independent random predictors $X_{i}\sim N(0,\Sigma)$ in our
theoretical analysis. This is less restrictive than the assumption of
$X_{i}\sim N(0,I_{p})$ in Yi and Caramanis (2015) and Balakrishnan et al.
(2017) in that we allow the predictor to be correlated. In this section, we
consider $K=2$, which is a common assumption in theoretical analysis for high-
dimensional EM algorithm (Yi and Caramanis, 2015; Cai et al., 2019; Zhang et
al., 2020, e.g.). We re-define $\theta=\\{\omega_{1},\beta_{1},\beta_{2}\\}$
since $\omega_{2}=1-\omega_{1}$. Let ${\theta}^{*}$ be the true value of
${\theta}$, and $\widehat{{\theta}}^{(t)}$ be the estimate of ${\theta}$ at
the $t$-th step in Algorithm 1. The true parameter space we consider is

$\displaystyle{\Theta}^{*}=\\{{\theta}^{*}:\omega_{1}^{*}\in(c_{w},1-c_{w}),\|{\beta}_{k}^{*}\|_{0}\leq
s,\|{\beta}_{k}^{*}\|_{2}\leq M_{b},\ \mathrm{for}\ k=1,2\\}.$

This is a natural parameter space to consider. The condition
$\omega_{1}^{*}\in(c_{w},1-c_{w})$ guarantees the sample size from each latent
class is large enough. Condition on $\|{\beta}_{k}^{*}\|_{2}\leq M_{b}$ is
also similarly used in Yi and Caramanis (2015) under the data splitting
framework and is milder than a similar condition $\|{\beta}\|_{1}\leq M_{b}$
used in Cai et al. (2019), where the EM algorithm for Gaussian mixture model
without data splitting is studied.

In theoretical studies, we first assume that the true value of $\sigma^{2}$
denoted as $(\sigma^{\star})^{2}$ as a known parameter, namely, the input
$\sigma^{2}$ is $(\sigma^{\star})^{2}$ in Algorithm 1. Without loss of
generality, we assume $(\sigma^{\star})^{2}=1$. Treating $\sigma^{2}$ as known
is also common in state-of-the-art theoretical studies for mixture linear
regression (Yi and Caramanis, 2015; Balakrishnan et al., 2017; Zhang et al.,
2020, e.g.). Although we do not analyze ${\beta}_{k}$’s and $\sigma^{2}$
simultaneously, we investigate how the choice of $\sigma^{2}$ influence the
estimation of Algorithm 1 in Theorem 2.

Since $\sigma^{2}=1$ and $K=2$, we simplify $\mathbb{P}(W_{i}=1\mid
Y_{i},X_{i},{\theta})$ as

$\displaystyle\quad\eta_{i,1}(\theta)=1-\eta_{i,2}(\theta)=1/\big{[}1+(\omega_{2}/\omega_{1})\mathrm{exp}\\{({\beta}_{2}-{\beta}_{1})^{T}{X}_{i}\cdot(Y_{i}-{({\beta}_{1}+{\beta}_{2})^{T}{X}_{i}}/{2})\\}\big{]}.$

Then the following quantities are used repeatedly in our theoretical analysis:

$\displaystyle\widehat{\omega}_{k}({\theta})=\frac{1}{n}\sum_{i=1}^{n}\eta_{i,k}(\theta),\
\widehat{{\rho}}_{k}({\theta})=\frac{1}{n}\sum_{i=1}^{n}\eta_{i,k}(\theta){X}_{i}Y_{i},$
$\displaystyle\widehat{{\Sigma}}_{k}({\theta})=\frac{1}{n}\sum_{i=1}^{n}\eta_{i,k}(\theta){X}_{i}{X}_{i}^{T},\
{\omega}_{k}({\theta})=\mathrm{E}\big{\\{}\frac{1}{n}\sum_{i=1}^{n}\eta_{i,k}(\theta)\big{\\}},$
$\displaystyle{{\rho}_{k}}({\theta})=\mathrm{E}\big{\\{}\frac{1}{n}\sum_{i=1}^{n}\eta_{i,k}(\theta){X}_{i}Y_{i}\big{\\}},\
{\Sigma}_{k}({\theta})=\mathrm{E}\big{\\{}\frac{1}{n}\sum_{i=1}^{n}\eta_{i,k}(\theta){X}_{i}{X}_{i}^{T}\big{\\}},$

where the expectation is with respect to ${X}_{i}$ and $Y_{i}$,
$i=1,\ldots,n$. Then we let
$M({\theta})=\\{{\omega}_{k}({\theta}),{{\rho}_{k}}({\theta}),{\Sigma}_{k}({\theta}),k=1,2\\}$,
$M_{n}({\theta})=\\{\widehat{\omega}_{k}({\theta}),\widehat{{\rho}}_{k}({\theta}),\widehat{{\Sigma}}_{k}({\theta}),k=1,2\\}$,
and define $d_{2,s}\big{(}M({\theta}_{1}),M({\theta}_{2})\big{)}$ and
$d_{2}\big{(}M({\theta}_{1}),M({\theta}_{2})\big{)}$ as

$\displaystyle\max_{k=1,2}\\{|\omega_{k}({\theta}_{1})-\omega_{k}({\theta}_{2})|\lor\|{\rho}_{k}({\theta}_{1})-{\rho}_{k}({\theta}_{2})\|_{2,s}\lor\|({\Sigma}_{k}({\theta}_{1})-{\Sigma}_{k}({\theta}_{2})){\beta}_{2}^{*}\|_{2,s}\\},$
$\displaystyle\max_{k=1,2}\\{|\omega_{k}({\theta}_{1})-\omega_{k}({\theta}_{2})|\lor\|{\rho}_{k}({\theta}_{1})-{\rho}_{k}({\theta}_{2})\|_{2}\lor\|({\Sigma}_{k}({\theta}_{1})-{\Sigma}_{k}({\theta}_{2})){\beta}_{2}^{*}\|_{2}\\},$

respectively, which are distances between $M({\theta}_{1})$ and
$M({\theta}_{2})$.

Let
$\Delta=\sqrt{({\beta}_{2}^{*}-{\beta}_{1}^{*})^{T}{\Sigma}({\beta}_{2}^{*}-{\beta}_{1}^{*})}$,
which is a measure of the signal-to-noise ratio of the mixture linear
regression model. We define the contraction basin
$\mathcal{B}_{con}({\theta})$ as follows.

$\displaystyle\mathcal{B}_{con}({\theta})=\\{{\theta}:\omega_{k}\in(c_{0},1-c_{0}),\|{\beta}_{k}-{\beta}_{k}^{*}\|_{2}\leq
C_{b}\Delta,\ {\beta}_{k}-{\beta}_{k}^{*}\in\Gamma(s),\ \mathrm{for}\
k=1,2\\}.$

Intuitively, the contraction basin requires that ${\beta}_{k}$ is not far away
from the true parameter ${\beta}_{k}^{*}$. Under the technical conditions
shown later, an initialization $\widehat{{\theta}}^{(0)}$ falls in the
contraction basin can guarantee the subsequent estimators
$\widehat{{\theta}}^{(t)}$ in Algorithm 1 are all contained in the contraction
basin.

###  3.2 Main results

We first introduce some technical conditions before stating the theoretical
results.

  1. (C1)

The eigenvalues of ${\Sigma}$ satisfy that
$M_{1}\leq\lambda_{min}({\Sigma})\leq\lambda_{max}({\Sigma})\leq M_{2}$.

  2. (C2)

$n\gg s\log(p)$.

  3. (C3)

The signal-to-noise ratio $\Delta>C_{1}(c_{0})$ for a constant $C_{1}(c_{0})$
only depends on $c_{0}$, and $C_{b}<C_{2}(c_{0},M_{2})$ for a constant
$C_{2}(c_{0},M_{2})$ only depends on $c_{0},M_{2}$.

  4. (C4)

The initialization
$\widehat{{\theta}}_{0}^{(0)}=(\widehat{\omega}_{1}^{(0)},\widehat{{\beta}}_{1}^{(0)},\widehat{{\beta}}_{2}^{(0)})\in\mathcal{B}_{con}({\theta}^{*})$.

Condition (C1) is a standard assumption on the covariance ${\Sigma}$ in high-
dimensional statistics (Bickel and Levina, 2008; Cai et al., 2011). Condition
(C2) is a common assumption in high dimensions on the relationship among
$(n,p,s)$ to guarantee consistent estimation (Meinshausen and Yu, 2009, e.g).
In particular, it implies that the restrictive eigenvalue condition
$\inf_{\mu\in\Gamma(s)\cap\mathcal{S}^{p-1}}\\{\mu^{T}(\sum_{i=1}X_{i}X_{i}^{T}/n)\mu\\}>\tau_{0}$
holds for a positive generic constant $\tau_{0}$ with high probability, and is
used for proving the concentration of $\widehat{{\beta}}_{k}^{(t)}$ in the
$t$-th iteration. Condition (C3) has two requirements. The first one is that
the signal-to-noise ratio is larger than a universal constant that does not
depend on $n$ and $p$ so that the two mixtures are distinguishable. This
requirement was also previously used in mixture linear model (e.g., Yi and
Caramanis, 2015; Balakrishnan et al., 2017; Zhang et al., 2020). The second
one is that, for the parameter ${\beta}_{k}$ in the contraction basin, the
distance $\|{\beta}_{k}-{\beta}_{k}^{*}\|_{2}$ is bounded by the signal-to-
noise ratio multiplied by a generic constant independent of $n$ and $p$. This
requirement makes all the ${\beta}_{k}$ in the contraction basin not too far
away from the truth ${\beta}_{k}^{*}$. Condition (C4) ensures that the
initialization is in the contraction basin. The contraction and concentration
properties shown later guarantee that the estimates in each step of Algorithm
1 stays in the contraction basin.

Next, we present two lemmas about the linear convergence of the population EM
updates and the concentration of the sample estimation to the population one
in each EM iteration. The following two lemmas together with a key technical
Lemma S.10 in the Supplement are highly non-trivial and serve as the building
blocks of the main theory for Algorithm 1.

######  Lemma 2.

Under conditions (C1) and (C3), if $\theta\in\mathcal{B}_{con}({\theta}^{*})$,
then

$d_{2}\big{(}M({\theta}),M({\theta}^{*})\big{)}\leq\kappa_{0}\big{(}|\omega_{1}({\theta})-\omega_{1}^{*}|\lor\|{\beta}_{1}-{\beta}_{1}^{*}\|_{2}\lor\|{\beta}_{2}-{\beta}_{2}^{*}\|_{2}\big{)}.$

for some $0<\kappa_{0}<\frac{1}{2\lor(64/\tau_{0})}$.

######  Lemma 3.

Suppose that ${\theta}^{*}\in{\Theta}^{*}$. Under condition (C1), there exists
a constant $C_{con}>0$, such that with probability at least $1-4p^{-1}$,

$\displaystyle\sup_{{\theta}\in\mathcal{B}_{con}({\theta}^{*})}d_{2,s}(M({\theta}),M_{n}({\theta}))\leq
C_{con}\sqrt{\frac{s\mathrm{log}(n)^{2}\mathrm{log}(p)}{n}}$

Intuitively, Lemma 2 shows the computational contraction of Algorithm 1. It
implies that, in each EM iterations, the expectations of the updated
estimators converge to the true parameters at a linear rate. On the other
hand, Lemma 3 establishes the statistical convergence rate of estimators to
their expectations in each EM iteration. When the iteration steps are large
enough, the computational error will be dominated by the statistical error,
which means that further iterations can not improve the statistical
convergence rate of the algorithm. Cai et al. (2019) also proved similar
lemmas under the Gaussian mixture model, but our proof is more challenging, as
we are interested in the mixture linear regression model. The unboundedness in
both $X$ and $Y$ makes $M(\theta)$ more complicated and $M_{n}(\theta)$ have
heavier tails.

Thanks to Lemmas 2 and 3, we can show the following result for Algorithm 1.

######  Theorem 1.

Under conditions (C1)–(C4), there exists a constant $0<\kappa<1/2$, such that
$\widehat{{\beta}}_{k}^{(t+1)}$ obtained by Algorithm 1 satisfies, with
probability $1-4p^{-1}$,

$\displaystyle\|\widehat{{\beta}}_{k}^{(t+1)}-{\beta}_{k}^{*}\|_{2}$
$\displaystyle=O\Big{(}\kappa^{t}(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{{\beta}}_{1}^{(0)}-{\beta}_{1}^{*}\|_{2}\lor\|\widehat{{\beta}}_{2}^{(0)}-{\beta}_{2}^{*}\|_{2})+\sqrt{\frac{s\mathrm{log}(n)^{2}\mathrm{log}(p)}{n}}\Big{)}.$

Consequently, for
$t\geq\\{-\mathrm{log}(\kappa)\\}^{-1}\mathrm{log}\\{n(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{{\beta}}_{1}^{(0)}-{\beta}_{1}^{*}\|_{2}\lor\|\widehat{{\beta}}_{2}^{(0)}-{\beta}_{2}^{*}\|_{2})\\}$,

$\|\widehat{{\beta}}_{k}^{(t+1)}-{\beta}_{k}^{*}\|_{2}=O\Big{(}\sqrt{\frac{s\log(n)^{2}\mathrm{log}(p)}{n}}\Big{)}.$

We make several remarks on Theorem 1. Firstly, compared with existing results
(Yi et al., 2014; Yi and Caramanis, 2015; Balakrishnan et al., 2017) requiring
${\beta}_{2}=-{\beta}_{1}$ and ${X}\sim N(0,I_{p})$, our model settings are
more general. On one hand, note that ${\beta}_{2}=-{\beta}_{1}$ is not just a
location shift of the response variable. As an illustration, when
${\beta}_{1}=(1,1,0)^{T}$ and ${\beta}_{2}=(0.5,2,0)^{T}$, a simple location
shift cannot reduce the model to the case where ${\beta}_{2}=-{\beta}_{1}$. By
removing the assumption that ${\beta}_{2}=-{\beta}_{1}$, our theory is
applicable to a larger model space. On the other hand, the predictors are not
likely to be uncorrelated in practice. Thus, it is meaningful to extend the
condition ${X}\sim N(0,I_{p})$ to ${X}\sim N(0,{\Sigma})$. Also, although we
adopt the group lasso penalty in the algorithm, the theoretical results can be
naturally extended to the penalties that are decomposable (Negahban et al.,
2012), such as the lasso penalty popular in the literature.

Moreover, unlike the extensive literature about the sample-splitting EM
algorithms for mixture linear regression (Yi et al., 2014; Yi and Caramanis,
2015; Zhang et al., 2020), the convergence result in Theorem 1 does not
require sample splitting. Sample splitting is not desirable, especially when
we have a small sample size. Splitting a limited number of observations into
$T$ batches decreases the estimation efficiency, makes the estimation less
stable and is rarely used in practice. Hence, it is meaningful to develop
theoretical results without data splitting for mixture linear regression, as
those in Theorem 1. To our best knowledge, Theorem 1 is the first theoretical
result for the high-dimensional EM algorithm of mixture linear regression
without data splitting. Also note that the convergence rate we obtain is
nearly optimal. When the latent random variables $W_{i}$’s are unknown, the
optimal rate $\sqrt{s\mathrm{log}(p)/n}$ (Ye and Zhang, 2010, e.g.), while
when $W_{i}$’s are unknown, Zhang et al. (2020) gives an estimation rate of
$\sqrt{s\mathrm{log}(p)\mathrm{log}{n}/n}$ with sample splitting. Our result
is slightly slower than these rates by the factors of $\mathrm{log}{n}$ and
$\mathrm{log}^{1/2}{n}$, respectively. The additional $\mathrm{log}(n)$ terms
are the price of no data splitting. Technically, to prove the convergence
results without data splitting we have to bound the tails of the supremum of
unbounded random processes, which is much more challenging than bounding the
random variables.

Now we turn to the effect of the misspecification of $\sigma^{2}$. In most
existing theoretical analysis including this one, $\sigma^{2}$ is usually
treated as a known parameter in theoretical studies for the mixture linear
regression model. There are two reasons for this treatment. One one hand, the
regression coefficients are of primary interest in regression models. On the
other, statistical analysis for mixture regression model with known
$\sigma^{2}$ is already challenging. However, in practice $\sigma^{2}$ is
almost never known. Statisticians often plug in an estimated value of
$\sigma^{2}$, which differs from $\sigma^{2}$. Yet it is unclear how the
misspecification affects the estimation of the mixture linear regression. In
the following theorem, we obtain a non-asymptotic convergence result for
Algorithm 1 with misspecified $\sigma^{2}$, which indicates that although
misspecified $\sigma^{2}$ brings bias to the estimation, it usually has a
minor influence on the mixture linear regression model with a large signal-to-
noise ratio.

######  Theorem 2.

Let $\widetilde{{\beta}}_{k}^{(t+1)}$ be the estimation of Algorithm 1 in the
$(t+1)$-th EM step where $\sigma^{2}$ may not equal to $(\sigma^{\star})^{2}$.
Under Conditions (C1)-(C4) and
$(\sigma^{\star})^{2}/2<\sigma^{2}<2(\sigma^{\star})^{2}$, with probability
$1-4p^{-1}$,

$\displaystyle\|\widetilde{{\beta}}_{k}^{(t+1)}-{\beta}_{k}^{*}\|_{2}$
$\displaystyle=O\Big{(}\kappa^{t}(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{{\beta}}_{1}^{(0)}-{\beta}_{1}^{*}\|_{2}\lor\|\widehat{{\beta}}_{2}^{(0)}-{\beta}_{2}^{*}\|_{2})+\sigma^{\star}(\Delta/\sigma^{\star})^{-2}|1-\frac{(\sigma^{\star})^{2}}{\sigma^{2}}|$
$\displaystyle+\sigma^{\star}\sqrt{\frac{s\mathrm{log}(n)^{2}\mathrm{log}(p)}{n}}\Big{)}.$

Consequently, for
$t\geq\\{-\mathrm{log}(\kappa)\\}^{-1}\mathrm{log}\\{n(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{{\beta}}_{1}^{(0)}-{\beta}_{1}^{*}\|_{2}\lor\|\widehat{{\beta}}_{2}^{(0)}-{\beta}_{2}^{*}\|_{2})\\}$,

$\|\widetilde{{\beta}}_{k}^{(t+1)}-\beta_{k}^{*}\|_{2}\leq
c_{1}\sigma^{\star}(\Delta/\sigma^{\star})^{-2}|1-\frac{(\sigma^{\star})^{2}}{\sigma^{2}}|+c_{2}\sigma^{\star}\sqrt{\frac{s\mathrm{log}(n)^{2}\mathrm{log}p}{n}}.$

Theorem 2 considered the case where $(\sigma^{\star})^{2}$ is misspecified as
$\sigma^{2}$ in Algorithm 1. The only difference between the convergence rates
in Theorems 1 and 2 is the bias term
$\sigma^{\star}(\Delta/\sigma^{\star})^{-2}|1-{(\sigma^{\star})^{2}}/{\sigma^{2}}|$.
It is a linear function of $|1-(\sigma^{\star})^{2}/\sigma^{2}|$ and is
exactly 0 when $\sigma^{2}=(\sigma^{\star})^{2}$. Although not using
$(\sigma^{\star})^{2}$ in Algorithm 1 may return a biased estimation, the bias
is small when $\sigma^{2}$ is close to $(\sigma^{\star})^{2}$. Note that the
bias term is also a function of $(\Delta/\sigma^{\star})^{-2}$. When the
signal-to-noise ratio $\Delta$ is large, the bias term is small regardless of
the value of $\sigma^{2}$. Specifically, if $\Delta\gg
O(n^{1/4}/[s\mathrm{log}(n)^{2}\mathrm{log}p]^{1/4})$, then the bias term is
ignorable compared with $\sqrt{s\mathrm{log}(n)^{2}\mathrm{log}p/n}$. For
large $p$, small $n$ cases, the requirement $\Delta\gg
O(n^{1/4}/[s\mathrm{log}(n)^{2}\mathrm{log}p]^{1/4})$ is usually not strict.
Intuitively, when $\Delta\to\infty$, the mixtures can be easily identified,
which makes the mixture linear regression model reduce to linear regression
models. It is well known that the choice of $\sigma^{2}$ in the linear
regression model does not affect the maximum likelihood estimation for the
coefficients. Thus, when $\Delta\to\infty$, the bias term disappears. Our
numerical studies also demonstrate that $\sigma^{2}$ has a minor influence on
the performance of Algorithm 1 for relatively large $\Delta$. Please see
Section 4.3 for details. We remark that the condition
$(\sigma^{\star})^{2}/2<\sigma^{2}<2(\sigma^{\star})^{2}$ is not necessary for
the algorithm in practice, but only a requirement in theory due to technical
reasons.

##  4 Simulation studies

###  4.1 Simulation set-up

In this section, we investigate the empirical performance of Algorithm 1. For
practical initialization, we start with lasso regression (Tibshirani, 1996) of
$Y$ on $X$ to roughly select important predictors and then apply the tensor
power method (Anandkumar et al., 2014) on the selected variables to get
initializations for ${\beta}_{k}$ and mixing proportion $\omega_{k}$,
$k=1,\cdots,K$. We use $T=20$ as the maximum number of iterations in Algorithm
1, and stop the iterations when
$\|(\widehat{\beta}^{(t+1)}_{1},\ldots,\widehat{\beta}^{(t+1)}_{K})-(\widehat{\beta}^{(t)}_{1},\ldots,\widehat{\beta}^{(t)}_{K})\|_{F}<10^{-3}$
. We use a single tuning parameter $\lambda=\lambda_{n}^{(t)}$ for all
$t=0,1,\dots,T$ and choose $\lambda$ based on the Bayesian information
criterion.

We include the following methods: 1) _Oracle_ , we fit the standard EM
algorithm on the true subset of $s$ relevant predictors. 2) _Initial_ , after
using the tensor power method to get the initialization for $\beta_{k}$,
$k=1,\cdots,K$, we plug them into (5) to get an estimation for the weight
$\widehat{\eta}_{ik}$ to the $i$-th sample, $i=1,\cdots,n$. Then, we label the
$i$-th sample with $\operatorname*{argmax}_{k=1,\cdots,K}\widehat{\eta}_{ik}$.
After that, we fit lasso regressions separately for each estimated mixtures to
estimate $\beta_{k}$. Compared with the tensor power method, this method
returns sparse estimations. 3) _GLLiM_ , the Gaussian Local Linear Mapping EM
algorithm (Deleforge et al., 2015) that is implemented in the R package xLLiM.
4) _HDEM_ , our implementation of the high-dimensional EM algorithm proposed
by Zhang et al. (2020). 5) _PSEM_ , i.e., the post-selection EM algorithm, we
fit the standard EM algorithm on the selected variables from Algorithm 1. and
finally 6) _PEM_ , the penalized EM algorithm (Algorithm 1).

We consider the following simulation models, where we first generate
independent predictor $X_{i}\sim N(0,\Sigma)$ and error $\epsilon_{i}\sim
N(0,1)$ and then $Y_{i}$ follows from the mixture linear regression (1). For
all the four simulation examples, we fix the first $s=10$ coefficients in
$\beta_{k}$ to be nonzero and consider both $p=400$ and $p=1000$ settings. The
total sample size $n$ is set to be $400$ for models (M1)–(M3), where we have
two mixtures, and $600$ for model (M4), where we have three mixtures. The
symmetric mixture assumption does not hold in any of our models,
$\beta_{1}\neq-\beta_{2}$.

  * •

(M1) Two mixtures with $\omega_{1}=\omega_{2}=0.5$, and auto-regressive
covariance structure $[{\Sigma}]_{ij}=0.3^{|i-j|}$ for $i,j=1,\ldots,p$. The
nonzero coefficients of ${\beta}_{1}$ is generated independently from
$N(0,1)$; and the nonzero coefficients of ${\beta}_{2}$ is
${\beta}_{2j}={\beta}_{1j}+2\cdot\mathrm{sgn}({\beta}_{1j})$, $j=1,\ldots,s$.

  * •

(M2) Same as (M1) but with weaker
signals:${\beta}_{2j}={\beta}_{1j}+1\cdot\mathrm{sgn}({\beta}_{1j})$,
$j=1,\ldots,s$.

  * •

(M3) Same as (M1) but with a different covariance $\Sigma$ based on Erdós-
Rényi random graph. Let $\widetilde{{\Sigma}}=(\widetilde{\sigma}_{ij})$,
where $\widetilde{\sigma}_{ij}=u_{ij}\delta_{ij}$, $\delta_{ij}$ follows
Bernoulli$(0.1)$ distribution, and
$u_{ij}\sim\mathrm{Uniform}[0.5,1]\cup\mathrm{Uniform}[-1,-0.5]$. Then let
$\widetilde{{\Sigma}}_{1}=(\widetilde{{\Sigma}}+\widetilde{{\Sigma}}^{T})/2$
and
${\Sigma}^{*}=\widetilde{{\Sigma}}_{1}+\\{\max(-\lambda_{min}(\widetilde{{\Sigma}}_{1}),0)+0.05\\}I_{p}$.
Finally, ${\Sigma}^{*}$ is standardized to have $1$’s on the diagonal.

  * •

(M4) Three mixtures with $\omega_{1}=\omega_{2}=\omega_{3}=1/3$. The nonzero
elements of ${\beta}_{1}$ and $\beta_{3}$ are $-1$ and $5$, respectively, and
the nonzero elements of ${\beta}_{2}$ are evenly spaced between $1$ and $3$.

###  4.2 Simulation results

The performances of those methods are evaluated using the following criteria.
The parameter estimation errors are defined as
$\sqrt{\sum_{k=1}^{K}\|\widehat{{\beta}}_{k}-{\beta}_{k}\|_{2}^{2}}$ for
$\beta$ and $\sum_{k=1}^{K}|\widehat{\omega}_{k}-\omega_{k}|\times 100$ for
$\omega$. The mixture estimation error is defined as
$\sum_{i=1}^{n}I(\widehat{W}_{i}\neq W_{i})/n\times 100$, where
$\widehat{W}_{i}=\operatorname*{argmax}_{k=1,\ldots,K}\widehat{\eta}_{i,k}(\widehat{\theta})$
for $i=1,\dots,n$. For methods that performs variable selection, we also
recorded the true positive and false positive rates of variable selection.

M1 | $p=400$ | $p=1000$
---|---|---
$\beta_{k}$ | $\omega_{k}$ | $W_{i}$ | $\beta_{k}$ | $\omega_{k}$ | $W_{i}$
Oracle | 0.44 (0.004) | 6.08 (0.06) | 8.43 (0.01) | 0.39 (0.001) | 6.18 (0.05) | 8.40 (0.01)
PSEM | 0.57 (0.01) | 6.31 (0.06) | 9.17 (0.04) | 0.76 (0.01) | 6.28 (0.05) | 10.50 (0.08)
GLLiM | 9.63 (0.01) | 68.18 (0.23) | 33.23 (0.12) | 9.68 (0.01) | 74.20 (0.23) | 34.24 (0.13)
Initial | 5.14 (0.01) | 32.14 (0.20) | 41.75 (0.06) | 5.34 (0.01) | 35.84 (0.22) | 42.46 (0.05)
HDEM | 1.22 (0.01) | 8.60 (0.06) | 10.71 (0.06) | 1.93 (0.02) | 10.17 (0.12) | 15.34 (0.13)
PEM | 1.04 (0.01) | 6.67 (0.06) | 9.79 (0.04) | 1.26 (0.01) | 7.08 (0.05) | 10.91 (0.08)
M2 | $p=400$ | $p=1000$
Oracle | 0.44 (0.003) | 12.29 (0.08) | 16.21 (0.02) | 0.45 (0.003) | 11.56 (0.07) | 15.97 (0.02)
PSEM | 0.64 (0.003) | 13.67 (0.09) | 17.78 (0.04) | 0.70 (0.02) | 14.06 (0.09) | 18.15 (0.02)
GLLiM | 6.74 (0.01) | 67.01 (0.20) | 33.95 (0.11) | 6.77 (0.01) | 68.60 (0.21) | 33.04 (0.10)
Initial | 2.79 (0.04) | 46.08 (0.24) | 44.03 (0.04) | 2.99 (0.04) | 55.76 (0.24) | 44.75 (0.04)
HDEM | 1.26 (0.01) | 30.87 (0.12) | 22.58 (0.06) | 1.96 (0.01) | 41.59 (0.22) | 29.65 (0.11)
PEM | 1.03 (0.004) | 23.63 (0.13) | 19.95 (0.04) | 1.39 (0.01) | 33.36 (0.23) | 23.28 (0.09)
M3 | $p=400$ | $p=1000$
Oracle | 0.56 (0.01) | 8.16 (0.06) | 10.09 (0.01) | 0.48 (0.01) | 7.60 (0.06) | 9.70 (0.02)
PSEM | 0.67 (0.01) | 7.97 (0.06) | 11.78 (0.06) | 0.65 (0.01) | 7.81 (0.06) | 11.10 (0.05)
GLLiM | 9.65 (0.01) | 68.41 (0.20) | 33.72 (0.11) | 9.54 (0.01) | 67.52 (0.22) | 33.974 (0.10)
Initial | 5.02 (0.01) | 36.47 (0.20) | 43.41 (0.05) | 5.30 (0.01) | 43.71 (0.22) | 43.79 (0.05)
HDEM | 1.54 (0.01) | 11.82 (0.07) | 14.12 (0.09) | 2.07 (0.02) | 13.39 (0.10) | 17.68 (0.13)
PEM | 1.21 (0.01) | 10.67 (0.07) | 12.27 (0.06) | 1.39 (0.01) | 10.47 (0.07) | 12.10 (0.05)
M4 | $p=400$ | $p=1000$
Oracle | 1.03 (0.03) | 5.14 (0.03) | 7.57 (0.01) | 1.35 (0.04) | 5.28 (0.05) | 7.53 (0.01)
PSEM | 1.80 (0.05) | 5.21 (0.03) | 10.82 (0.12) | 3.56 (0.07) | 6.58 (0.06) | 16.53 (0.19)
GLLiM | 17.38 (0.01) | 52.15 (0.27) | 33.79 (0.06) | 17.40 (0.01) | 54.10 (0.23) | 34.99 (0.06)
Initial | 13.16 (0.02) | 18.56 (0.10) | 49.62 (0.06) | 13.54 (0.01) | 18.48 (0.09) | 51.21 (0.06)
HDEM | 8.87 (0.09) | 8.93 (0.07) | 30.02 (0.25) | 12.08 (0.07) | 19.28 (0.13) | 43.62 (0.24)
PEM | 2.43 (0.04) | 5.08 (0.03) | 11.49 (0.12) | 3.99 (0.06) | 6.32 (0.05) | 17.10 (0.19)

Table 1:  Average estimation errors based on 100 replicates (standard errors
in parentheses).  Figure 1: The estimation error of ${\beta}_{k}$’s based on
100 replicates.

| M1 $p=400$ | M1 $p=1000$ | M2 $p=400$ | M2 $p=1000$
---|---|---|---|---
TPR | FPR | TPR | FPR | TPR | FPR | TPR | FPR
Initial | 83.2 (0.21) | 8.3 (0.05) | 76.4 (0.23) | 3.2 (0.02) | 88.3 (0.19) | 8.9 (0.05) | 79.6 (0.21) | 4.3 (0.03)
HDEM | 92.7 (0.09) | 1.4 (0.03) | 86.6 (0.15) | 1.5 (0.03) | 89.8 (0.01) | 1.9 (0.02) | 81.4 (0.25) | 1.2 (0.01)
PEM | 100 (0) | 0.9 (0.02) | 100 (0) | 0.7 (0.02) | 100 (0) | 1.2 (0.01) | 100 (0) | 0.8 (0.01)
| M3 $p=400$ | M3 $p=1000$ | M4 $p=400$ | M4 $p=1000$
Initial | 89.1 (0.16) | 0.5 (0.05) | 73.4 (0.25) | 3.9 (0.03) | 47.5 (0.31) | 4.0 (0.05) | 38.3 (0.28) | 1.7 (0.02)
HDEM | 95.0 (0.06) | 2.6 (0.05) | 80.2 (0.12) | 1.7 (0.03) | 83.4 (0.24) | 12.6 (0.15) | 64.7 (0.33) | 6.7 (0.06)
PEM | 100 (0) | 1.3 (0.02) | 100 (0) | 0.5 (0.01) | 100 (0) | 2.7 (0.11) | 100 (0) | 3.8 (0.09)

Table 2:  Average true positive rate (TPR) and false positive rate (FPR) for
variable selection based on 100 replicates (standard errors in parentheses).

The estimation results are summarized in Table 1, particularly for $\beta$’s
estimation we further plot the results from 100 replicates in Figure 1. As
expected, the _oracle_ method, i.e., the standard EM algorithm applied to the
truly relevant $s=10$ predictors, has the best performance. On the other hand,
the method that does not perform variable selection, GLLiM, failed for all the
simulation settings. The proposed method, _PEM_ , has encouraging performances
for all the simulation models; overall, it is comparable to the _oracle_
method and has a slight edge over the very recent method _HDEM_ (Zhang et al.,
2020). The advantage of our method over _HDEM_ , which updates each
$\beta_{k}$ separately via lasso penalized estimation in the M-step, can be
explained by the group-wise penalization and joint estimation of all
$\beta_{k}$’s in the M-step. In model (M4) with three mixtures, our approach
of joint estimation has even more gains in estimation accuracy.

We also note that Algorithm 1 substantially improves over the _Initial_
values. Even when the initialization is quite far from the truth, the proposed
algorithm improves over iterations in terms of all the evaluation criteria.
This is partly explained by the guaranteed monotonicity of our iterative
algorithm in the penalized log-likelihood (see Lemma S.1.). Finally, the post-
selection EM algorithm, _PSEM_ , further removes the bias in the penalized
estimation and improves our method one-step closer to the _oracle_ method.
This is because our method has excellent variable selection accuracy, as
summarized in Table 2.

###  4.3 Simulations for misspecified $\sigma^{2}$ in Algorithm 1

In this section, we show some simulation results for using different
$\sigma^{2}$ in Algorithm 1. The data is generated from (1). We set $K=2$,
$\omega_{1}=\omega_{2}=1/2$, $p=400$, $n=400$, and the true value of
$\sigma^{2}$ to be 1. The nonzero coefficients of ${\beta}_{1}$ is generated
independently from $N(0,1)$; and the nonzero coefficients of ${\beta}_{2}$ is
${\beta}_{2j}={\beta}_{1j}+\delta\cdot\mathrm{sgn}({\beta}_{1j})$,
$j=1,\ldots,p$. The signal strength $\delta$ takes value from $\\{0.75,1,2\\}$
and the input variance $\sigma^{2}$ in Algorithm 1 takes value from
$\\{0.5,0.75,1,1.5,2\\}$.

Figure 2: Coefficient estimation errors for Algorithm 1 with different input $\sigma^{2}$ based on 100 replicates.  | $\sigma^{2}=0.5$ | $\sigma^{2}=0.75$ | $\sigma^{2}=1$ | $\sigma^{2}=1.5$ | $\sigma^{2}=2$
---|---|---|---|---|---
$\delta=2$ | 1.19 (0.08) | 1.05 (0.03) | 1.05 (0.03) | 1.06 (0.03) | 1.13 (0.07)
$\delta=1$ | 1.05 (0.05) | 1.01 (0.04) | 1.03 (0.04) | 1.06 (0.03) | 1.23 (0.05)
$\delta=0.75$ | 1.50 (0.09) | 1.51 (0.10) | 1.60 (0.09) | 2.43 (0.12) | 3.28 (0.07)
Table 3:  Coefficient estimation errors for Algorithm 1 with different input
$\sigma^{2}$ based on 100 replicates. (standard errors in parentheses).

From Table 3 and Figure 2, when $\delta$ equals 1 and 2, using all the five
$\sigma^{2}$ in Algorithm 1 returns good estimations. Although when
$\sigma^{2}$ is 0.5 or 2, the performances are slightly worse, but can still
capture the mixture information and give accurate enough estimates. The
results are consistent with Theorem 2. For mixture linear regression models
with a large signal-to-noise ratio, the choice of $\sigma^{2}$ has a minor
influence on the performance of Algorithm 1. When $\delta$ is 0.75, the
signal-to-noise ratio is not large enough, the performance of Algorithm 1
becomes worse, no matter the value of $\sigma^{2}$. This is mainly caused by
the poor quality of the initialization. It deserves to notice that when
$\sigma^{2}$ equals 0.5 or 0.75, the performance of the algorithm is slightly
better than taking $\sigma^{2}=1$. This result indicates that when the signal-
to-noise ratio is relatively small, Algorithm 1 can be less sensitive to the
initialization if $\sigma^{2}$ takes a smaller value than the true one.

Algorithm 2 Group lasso penalized EM algorithm for model (9)

  1. Input:

Initial values $\widehat{\omega}_{k}^{(0)}$, $\widehat{\beta}_{k}^{(0)}$, for
$k=1,\cdots,K$, maximum iteration number $T$, data
$\\{X_{i},Y_{i};i=1,\ldots,n\\}$, and initial tuning parameter

$\displaystyle\lambda_{n}^{(0)}=C_{1}(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{\beta}_{1}^{(0)}-\beta_{1}^{*}\|_{F}\lor\cdots\lor\quad\|\widehat{\beta}_{K}^{(0)}-\beta_{K}^{*}\|_{F})/{\sqrt{s}}+C_{\lambda}\sqrt{\mathrm{log}(n)^{2}\mathrm{log}(pq)/s}$

for some positive constants $C_{1}$ and $C_{\lambda}$.

  2. Iterate:

For $t=0,\cdots,T-1$, do the following steps until convergence.

     * •

For $i=1,\cdots,n$, let

$\displaystyle\quad\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})={\widehat{\omega}_{k}^{(t)}}/\Big{(}\widehat{\omega}_{k}^{(t)}+\sum_{k^{\prime}\neq
k}\widehat{\omega}_{k^{\prime}}^{(t)}\mathrm{exp}\big{\\{}X_{i}^{T}(\widehat{\beta}_{k^{\prime}}^{(t)}-\widehat{\beta}_{k}^{(t)})^{T}{\Sigma}_{y}^{-1}\big{(}Y_{i}-{(\widehat{\beta}_{k}^{(t)}+\widehat{\beta}_{k^{\prime}}^{(t)})^{T}X_{i}}/{2}\big{)}\big{\\}}\Big{)}.$

     * •

For $k=1,\cdots,K$, update

$\displaystyle\widehat{\omega}_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)}),$
$\displaystyle\widehat{\rho}_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\big{(}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})X_{i}Y_{i}^{T}{\Sigma}_{y}^{-1}\big{)},$
$\displaystyle\widehat{\Sigma}_{k}^{(t+1)}=\frac{1}{n}\sum_{i=1}^{n}\big{(}\widehat{\eta}_{i,k}(\widehat{\theta}^{(t)})X_{i}X_{i}^{T}\big{)},$

and update $\widehat{\beta}_{k}^{(t+1)}$ by minimizing

$\displaystyle\sum_{k=1}^{K}\mathrm{tr}\big{(}{\Sigma}_{y}^{-1}\beta_{k}^{T}\widehat{\Sigma}_{k}^{(t+1)}\beta_{k}\big{)}-2\sum_{k=1}^{K}\mathrm{tr}\big{(}(\widehat{\rho}_{k}^{(t+1)})^{T}\beta_{k}\big{)}+\lambda_{n}^{(t+1)}\sum_{j=1}^{p}\sum_{l=1}^{q}\sqrt{\sum_{k=1}^{K}(\beta_{k})_{jl}^{2}}.$

with
$\lambda_{n}^{(t+1)}=\kappa\lambda_{n}^{(t)}+C_{\lambda}\sqrt{\mathrm{log}(pq)\mathrm{log}(n)^{2}/n}$,
where $\kappa\in(0,1/2)$ and $(\beta_{k})_{jl}$ is the $(j,l)$-th element of
${\beta}_{k}$.

  3. Output:

$\widehat{\theta}^{(t+1)}=\\{\widehat{\omega}_{1}^{(t+1)},\ldots,\widehat{\omega}_{K}^{(t+1)},\widehat{\beta}_{1}^{(t+1)},\ldots,\widehat{\beta}_{K}^{(t+1)}\\}$

##  5 Extension to multivariate mixture linear regression

###  5.1 Model and algorithm

In this section, we consider a generalization of the mixture linear regression
model, where $Y$ is a $q$-dimensional response. Specifically, we consider the
model

$\mathbb{P}(W=k)=\omega_{k},\quad Y\mid(X,W=k)\sim
N(\beta_{k}^{T}X,\Sigma_{y}),$ (9)

where $\beta\in\mathbb{R}^{p\times q}$ and $\Sigma_{y}\in\mathbb{R}^{q\times
q}$ is a symmetric and positive definite matrix. We will answer two questions:
Does the penalized EM algorithm still work and have similar convergence
results for multiple response cases? If so, what is the advantage of
considering $q$ responses simultaneously than $q$ mixture linear regression
problems separately?

As the mixture linear regression model, to handle the high-dimensionality of
$X$, we assume the matrix coefficients $\beta_{k}$, $k=1,\cdots,K$, to be
sparse, and the sparsity patterns are the same for all the mixtures. We
redefine $S=\\{(i,j):(\beta_{k})_{ij}\neq 0,\ i=1\cdots,p,\ j=1,\cdots,q\\}$
and $s=|S|_{0}$ be the cardinality of $S$. We allow $p\gg n$ and $q$ to grow
linearly with $s$.

In model (9), our parameter of interest is
${\theta}=\\{\omega_{1},\cdots,\omega_{k},{\beta}_{1},\cdots,{\beta}_{K}\\}$.
The error covariance $\Sigma_{y}$ is treated as a known parameter. The
penalized EM algorithm for solving (9) is analogous to Algorithm 1 and is
summarized in Algorithm 2. Similar to Algorithm 1, the formula of
$\lambda_{n}^{(t)}$ is mainly for theoretical consideration. In practice, we
fix $\lambda_{n}^{(t)}=\lambda$ over all iterations.

###  5.2 Theoretical analysis

In this section, we consider the theoretical studies for Algorithm 2. We
assume that ${\Sigma}_{y}$ is a known parameter and without loss of
generality, we assume ${\Sigma}_{y}={I}_{q}$ (If ${\Sigma}_{y}$ is not
identity, treating ${\Sigma}_{y}^{-1/2}{\beta}_{k}$ as the target parameter
will turn the problem to the case ${\Sigma}_{y}={I}_{q}$). Recall that, even
for the univariate-response mixture linear regression, theoretical analysis
for unknown variance $\sigma^{2}$ is extremely challenging and has not been
considered in high dimensions. As a result, for multiple-response mixture
linear regression, it makes sense to assume that ${\Sigma}_{y}$ is known as a
necessary simplification. As such, we focus on the coefficients ${\beta}_{k}$
as parameters of interest. To the best of our knowledge, there is little work
for multiple-response mixture regression even in this simplified context.

We first define the following additional notations. For a matrix
${A}\in\mathbb{R}^{p\times q}$, $\|{A}\|_{F}$ is its Frobenius norm and
$\|{A}\|_{F,s}=\sup_{{\mu}\in\mathbb{R}^{p\times
q},\|{\mu}\|_{F}=1,\mathrm{vec}({\mu})\in{\Gamma}(s)}\langle{A},{\mu}\rangle_{F}$.
The true parameter space we consider now is

${\Theta}^{*}=\\{{\theta}^{*}:\omega_{1}^{*}\in(c_{w},1-c_{w}),\|\mathrm{vec}({\beta}_{k}^{*})\|_{0}\leq
s,\|{\beta}_{k}^{*}\|_{F}\leq M_{b},\mathrm{for}\ k=1,2\\},$

and the signal-to-noise ratio is re-defined as
${\Delta}=\sqrt{\mathrm{tr}(({\beta}_{2}^{*}-{\beta}_{1}^{*})^{T}{\Sigma}({\beta}_{2}^{*}-{\beta}_{1}^{*}))}/q$.
Then we define $d_{F,s}\big{(}M({\theta}_{1}),M({\theta}_{2})\big{)}$ and
$d_{F}\big{(}M({\theta}_{1}),M({\theta}_{2})\big{)}$ as

$\displaystyle\max_{k=1,2}\\{|\omega_{k}({\theta}_{1})-\omega_{k}({\theta}_{2})|\lor\|{\rho}_{k}({\theta}_{1})-{\rho}_{k}({\theta}_{2})\|_{F,s}\lor\|({\Sigma}_{k}({\theta}_{1})-{\Sigma}_{k}({\theta}_{2})){\beta}_{2}^{*}\|_{F,s}\\},$
$\displaystyle\max_{k=1,2}\\{|\omega_{k}({\theta}_{1})-\omega_{k}({\theta}_{2})|\lor\|{\rho}_{k}({\theta}_{1})-{\rho}_{k}({\theta}_{2})\|_{F}\lor\|({\Sigma}_{k}({\theta}_{1})-{\Sigma}_{k}({\theta}_{2})){\beta}_{2}^{*}\|_{F}\\},$

respectively, and the new constriction basin $\mathcal{B}_{con}({\theta})$ as

$\mathcal{B}_{con}({\theta})=\\{{\theta}:\omega_{k}\in(c_{0},1-c_{0}),\|{\beta}_{k}-{\beta}_{k}^{*}\|_{F}\leq
C_{b}\Delta,\ \mathrm{vec}({\beta}_{k}-{\beta}_{k}^{*})\in\Gamma(s),\
\mathrm{for}\ k=1,2\\}.$

We require the same technical conditions as the mixture linear regression with
a modification to Condition (C3):

  1. (C3’)

The new signal-to-noise ratio
$\Delta=\sqrt{\mathrm{tr}(({\beta}_{2}^{*}-{\beta}_{1}^{*})^{T}{\Sigma}({\beta}_{2}^{*}-{\beta}_{1}^{*}))}/q>C_{1}(c_{0})$
for a constant $C_{1}(c_{0})$ only depends on $c_{0}$, and
$C_{b}<C_{2}(c_{0},M_{2})$ for a constant $C_{2}(c_{0},M_{2})$ only depends on
$c_{0},M_{2}$.

The new signal-to-noise ratio reveals a fundamental difference of considering
$q$ responses together than $q$ mixture linear regressions separately. To see
this, define $\beta_{k,j}$ as the $j$-th column of $\beta_{k}$, i.e, the
regression coefficient of $(Y_{j},{X})$. Further define the individual signal-
to-noise-ratio as
$\Delta_{j}=\sqrt{({\beta}_{1,j}-{\beta}_{2,j})^{T}{\Sigma}({\beta}_{1,j}-{\beta}_{2,j})}$.
If we apply Algorithm 1 to $(Y_{j},X)$ for $j=1,\ldots,q$, then we need all of
$\Delta_{j}$ to be bounded below, as discussed in Section 3. In other words,
when some $\Delta_{j}$ are small and the two mixtures are not well separated
on this coordinate, Algorithm 1 may not be able to estimate the corresponding
coefficients. However, if we consider all the responses simultaneously in
multivariate mixture linear regression, we only require
$\Delta=\sqrt{\sum_{j=1}^{q}\Delta_{j}^{2}}/q$, the average of $q$ signal-to-
noise ratios, to be greater than a constant. We can allow some responses to
have small signal-to-noise ratios without rendering Algorithm 2 inapplicable.
This is because all the response shares the common latent structure. The
responses with large signal-to-noise ratios can help enhance the
identification of the mixtures for the responses with small signal-to-noise
ratios.

As a brief numerical illustration, we generated simulated datasets from (9)
with $K=2$, $n=400$, $p=100$, $s=10$, $q=2$, ${\Sigma}_{y}={I}_{q}$ and
$[{\Sigma}]_{ij}=0.3^{|i-j|}$. The first 5 rows of ${\beta}_{k}$, $k=1,2$, are
set to be nonzero. More specifically, the nonzero coefficients of the first
and second columns of ${\beta}_{1}$ are generated from $N(0,1)$; and the
nonzero coefficients of ${\beta}_{2}$ are
$({\beta}_{2})_{j1}=({\beta}_{1})_{j1}+2\cdot\mathrm{sgn}(({\beta}_{1})_{j1})$
and
$({\beta}_{2})_{j2}=({\beta}_{1})_{j2}+\delta\cdot\mathrm{sgn}(({\beta}_{1})_{j2})$.
The signal strength is 2 for the first response and $\delta$ for the second
one. We vary $\delta$ from 0.5 to 2 with increments of 0.5.

Figure 3: The reported results for each signal strength are medians of 100
replicates. In the left panel, ${\beta}_{:j}(M)$ represents the estimation
error for the $j$-th column of $\beta$ using Algorithm 2; ${\beta}_{:j}$
represents the estimation error for the $j$-th column of $\beta$ using
Algorithm 1 separately for each response. In the right panel, MMLR represents
the estimation error for the mixtures using Algorithm 2; MLR-Response1 and
MLR-Response2 represent the estimation error for the mixtures using Algorithm
1 separately for the first and second element of the response, respectively.

From Figure 3, we can see that the coefficient estimation errors for Algorithm
2 remain at a low level for all signal strengths and the mixture estimation
error gains slightly improvement when the signal strength of the second
response increases. When the signal strength of the second response is weak,
the strong signal of the first response assists the recovery of the mixtures
and helps to increase the coefficient estimation accuracy of the second
coefficient substantially. As a comparison, when using Algorithm 1 separately
for each response, the estimation errors for the second coefficient and
mixture increase quickly when the signal strength decreases. In terms of the
mixture estimation error, it even performs like random guessing when the
signal strength is 0.5.

The following lemmas and theorem rigorously show the statistical convergence
for Algorithm 2.

######  Lemma 4.

Under conditions (C1),(C3’), if $q=O(s)$ and
$\theta\in\mathcal{B}_{con}({\theta}^{*})$, then

$d_{F}\big{(}M({\theta}),M({\theta}^{*})\big{)}\leq\kappa_{0}\big{(}|\omega_{1}({\theta})-\omega_{1}^{*}|\lor\|{\beta}_{1}-{\beta}_{1}^{*}\|_{F}\lor\|{\beta}_{2}-{\beta}_{2}^{*}\|_{F}\big{)}.$

for some $0<\kappa_{0}<\frac{1}{2\lor(64/\tau_{0})}$.

######  Lemma 5.

Suppose that ${\theta}^{*}\in{\Theta}^{*}$. Under condition (C1), there exists
a constant $C_{con}>0$, such that with probability at least $1-4(pq)^{-1}$,

$\displaystyle\sup_{{\theta}\in\mathcal{B}_{con}({\theta}^{*})}d_{F,s}(M({\theta}),M_{n}({\theta}))\leq
C_{con}\sqrt{\frac{s\mathrm{log}(n)^{2}\mathrm{log}(pq)}{n}}$

Lemmas 4 & 5 can be interpreted similarly to Lemmas 1 and 2, respectively.
However, since we are now interested in multivariate response, we further
assume condition $q=O(s)$ in addition to the technical conditions (C1) to
(C4). The assumption on $q$ results from the signal-to-noise ratio. Note that
$\Delta=O(s/q)$. In order for the condition $\Delta>C_{1}$ to hold, we need
$q=O(s)$. However, when this assumption holds, $q$, i.e, the dimensionality of
${Y}$, has a relatively small effect on the estimation error, as we only have
an additional factor of $\mathrm{log}{q}$ in Lemma 5 in comparison to Lemma 3.

Combine Lemmas 4 & 5 and we have the following theorem.

######  Theorem 3.

Under conditions (C1), (C2), (C3’), and(C4) and $q=O(s)$, there exists a
constant $0<\kappa<1/2$, such that $\widehat{{\beta}}_{k}^{(t+1)}$ obtained by
Algorithm 2 satisfies, with probability greater than $1-4(pq)^{-1}$,

$\displaystyle\|\widehat{{\beta}}_{k}^{(t+1)}-{\beta}_{k}^{*}\|_{F}$
$\displaystyle=O\Big{(}\kappa^{t}(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{{\beta}}_{1}^{(0)}-{\beta}_{1}^{*}\|_{F}\lor\|\widehat{{\beta}}_{2}^{(0)}-{\beta}_{2}^{*}\|_{F})+\sqrt{\frac{s\mathrm{log}(n)^{2}\mathrm{log}(pq)}{n}}\Big{)}.$

Consequently, for
$t\geq\\{-\mathrm{log}(\kappa)\\}^{-1}\mathrm{log}\\{n(|\widehat{\omega}_{1}^{(0)}-\omega_{1}^{*}|\lor\|\widehat{{\beta}}_{1}^{(0)}-{\beta}_{1}^{*}\|_{F}\lor\|\widehat{{\beta}}_{2}^{(0)}-{\beta}_{2}^{*}\|_{F})\\}$,

$\|\widehat{{\beta}}_{k}^{(t+1)}-{\beta}_{k}^{*}\|_{F}=O\Big{(}\sqrt{\frac{s\log(n)^{2}\mathrm{log}(pq)}{n}}\Big{)}.$

Again, the convergence rate for multivariate mixture linear regression model
is similar to that for the mixture linear regression model with an additional
$\mathrm{log}(q)$ term. The convergence rate grows slowly as a function of
$q$. However, theoretical requirement $\Delta>C_{1}$ and practical
consideration about the estimation for $\Sigma_{y}$ restrict $q$ to grow
linearly with $s$ and slower than $n$. It would be interesting to develop a
method that allows $q$ to grow faster, which we leave a challenging future
topic.

##  6 Real data analysis

The cancer cell line encyclopedia (CCLE) dataset contains 8-point dose-
response curves for 24 chemical compounds across over 400 cell lines, which is
a publicly available dataset at www.broadinstitute.org/ccle and is also
considered in Li et al. (2019). Because the cell lines are not consistent for
different chemical compounds, we consider chemical compounds: Lapatinib,
AZD6244, and PD-0325901, three popular chemical compounds for cancer
treatment. Analogous to Li et al. (2019), we use the area under the dose-
response curve, also known as the activity area, to measure the sensitivity of
a drug for each cell line. Besides the drug information, the data also
contains the expression data of 18,926 genes for each cell line. Aiming to
identify the genes sensitive to the chemical compounds, we treat the active
area as the response and the gene expressions as the predictor. After
identifying cell lines treated using all three chemical compounds, we get the
sample size $n=490$. We keep $p=500$ genes that are highly correlated with
three responses (for each gene, we calculated the sum of its absolute
correlations with the three responses). Due to the complexity of cancer, we
expect the data to be heterogeneous.

We consider $K=2,3,4,5,6$, and $8$ for the proposed algorithms and denote
individual mixture linear regression using Algorithm 1 as MLR and multivariate
mixture linear regression using Algorithm 2 as MMLR. In addition, we consider
the methods used in the simulation studies including HDEM, Initial method, and
PSEM and LASSO regression as competitors. We repeatedly split the whole data
set into $1:4$ ratios, and use $1/5$ of the data as the testing samples and
the rest as the training samples. The process is repeated 100 times. We report
the mean squared prediction errors. Please see Table 4 for the overall
prediction errors and Figures 4 and 5 for the boxplots. More detailed
comparison results for each response are presented in Section S.8 of the
Supplement.

When $K=5$, the overall mean squared prediction errors for MMLR is the
smallest. However, from Figure 4, we can see that MMLR is insensitive to $K$.
When $K=4,5$ and $6$, the mean squared prediction errors are very close. As a
comparison, MLR is more sensitive to $K$ and achieves the best prediction
error when $K=4$. Two possible reasons make MMLR perform better than MLR. The
first reason is that considering the three responses together results in a
larger signal-to-noise ratio, which facilitates the mixture identification and
coefficients estimation. The second reason is that, unlike MMLR, the estimated
mixtures are different for the three responses in MLR, which may reduce the
estimation accuracy. Although MLR performs slightly worse than MMLR, it still
outperforms the other competitors, especially the LASSO regression. It implies
that the data set is quite heterogeneous, the proposed mixture linear
regression approach improves the performance of a single linear regression
significantly.

When $K=4$ and $6$, we also consider the estimated label for each observation
based on the full data set. The number of observations in each estimated
mixture is reported in Table 5. Based on this table, we recommend using $K=2$,
3 and 3 for MLR to the three responses, respectively, and using $K=5$ for
MMLR. Except from the first response, the recommended $K$’s make MLR and MMLR
achieve the lowest or among the lowest prediction error. For the first
response, MLR achieves the smallest prediction error when $K=6$. However, it
only separates the observations into 2 mixtures. Note that the estimated label
for the $i$th observation is
$\operatorname*{argmax}_{k=1,\cdots,K}\widehat{\eta}_{i,k}(\widehat{\theta})$,
it can happen that no observation is classified into some mixtures. For the
first response, the two mixtures are more balanced with the increase of $K$,
which may be one reason makes the prediction error smaller. Another reason is
that the relationship between the first response and the predictors may be
complicated and non-linear, which can not be captured by a mixture linear
regression model. Hence, MLR tends to use more mixtures to approximate the
true relationship. Another phenomena we observe is that MMLR uses more
mixtures than MLR. This happens since the true mixtures can be different for
the three responses. By separating the samples into more mixtures, MMLR
becomes more accurate and robust.

K | MLR | MMLR | HDEM | Initial | PSEM | LASSO
---|---|---|---|---|---|---
2 | 0.413 (0.005) | 0.388 (0.004) | 0.420 (0.004) | 0.612 (0.009) | 1.461 (0.020) | 0.948 (0.026)
3 | 0.275 (0.004) | 0.276 (0.004) | 0.347 (0.004) | 0.493 (0.010) | 1.376 (0.021)
4 | 0.254 (0.005) | 0.238 (0.004) | 0.324 (0.004) | 0.443 (0.008) | 1.079 (0.020)
5 | 0.276 (0.007) |  0.237 (0.004) | 0.315 (0.005) | 0.416 (0.008) | 0.569 (0.012)
6 | 0.302 (0.007) | 0.239 (0.004) | 0.310 (0.004) | 0.395 (0.007) | 0.412 (0.008)
8 | 0.347 (0.006) | 0.243 (0.004) | 0.310 (0.005) | 0.357 (0.006) | 0.375 (0.007) |

Table 4:  Mean squared prediction errors based on 100 replicates (Standard
errors are given in parenthesis).

Figure 4: CCLE data: Boxplot for the mean squared prediction error  Figure 5: Overall prediction error for MLR and MMLR. The red boxes represent results for MLR. | K=4 | K=6
---|---|---
MLR (Lapatinib) | 0 | 0 | 147 | 343 | 0 | 0 | 0 | 0 | 162 | 328
MLR (AZD6244) | 0 | 87 | 162 | 241 | 0 | 0 | 0 | 85 | 187 | 218
MLR (PD-0325901) | 0 | 105 | 170 | 215 | 0 | 0 | 0 | 101 | 180 | 209
MMLR | 76 | 128 | 138 | 148 | 0 | 48 | 73 | 104 | 114 | 151
Table 5:  Numbers of samples in each estimated mixtures.

##  7 Discussion

The paper studies a group lasso penalized EM algorithm for high-dimensional
mixture linear regression. We obtained an encouraging non-asymptotic
convergence rate without data splitting and under a general model setting.
Although our theory is for normally distributed predictors and two-mixtures
regression, the penalized EM algorithm is applicable to essentially any
predictors and to more than two mixtures. We then extended the mixture linear
regression model and the penalized EM algorithm to the multivariate response
case and established its non-asymptotic theory.

The theoretical study is currently for the model with two mixtures. It
provides insights and foundations to the theoretical studies for cases with
$K\geq 3$. We leave it as a future work. Besides, the error is assumed to
follow a normal distribution, it is interesting to considering more general
error assumptions, such as the student t distribution, to achieve robustness.
Finally, it is also a meaningful future work to extend the studies to
generalized mixture linear regression models.

## References

  * (1)
  * Adamczak (2008) Adamczak, R. (2008), ‘A tail inequality for suprema of unbounded empirical processes with applications to markov chains’, Electronic Journal of Probability 13, 1000–1034.
  * Anandkumar et al. (2014) Anandkumar, A., Ge, R., Hsu, D., Kakade, S. M. and Telgarsky, M. (2014), ‘Tensor decompositions for learning latent variable models’, Journal of machine learning research 15, 2773–2832.
  * Balakrishnan et al. (2017) Balakrishnan, S., Wainwright, M. J. and Yu, B. (2017), ‘Statistical guarantees for the em algorithm: From population to sample-based analysis’, The Annals of Statistics 45(1), 77–120.
  * Bickel and Levina (2008) Bickel, P. J. and Levina, E. (2008), ‘Covariance regularization by thresholding’, The Annals of Statistics 36(6), 2577–2604.
  * Cai et al. (2011) Cai, T., Liu, W. and Luo, X. (2011), ‘A constrained l1 minimization approach to sparse precision matrix estimation’, Journal of the American Statistical Association 106(494), 594–607.
  * Cai et al. (2019) Cai, T. T., Ma, J. and Zhang, L. (2019), ‘Chime: Clustering of high-dimensional gaussian mixtures with em algorithm and its optimality’, The Annals of Statistics 47(3), 1234–1267.
  * Deleforge et al. (2015) Deleforge, A., Forbes, F. and Horaud, R. (2015), ‘High-dimensional regression with gaussian mixtures and partially-latent response variables’, Statistics and Computing 25(5), 893–911.
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M. and Rubin, D. B. (1977), ‘Maximum likelihood from incomplete data via the em algorithm’, Journal of the Royal Statistical Society: Series B (Methodological) 39(1), 1–22.
  * Hyun et al. (2020) Hyun, S., Cape, M. R., Ribalet, F. and Bien, J. (2020), ‘Modeling cell populations measured by flow cytometry with covariates using sparse mixture of regressions’, arXiv preprint arXiv:2008.11251 .
  * Khalili and Chen (2007) Khalili, A. and Chen, J. (2007), ‘Variable selection in finite mixture of regression models’, Journal of the american Statistical association 102(479), 1025–1038.
  * Klusowski et al. (2019) Klusowski, J. M., Yang, D. and Brinda, W. (2019), ‘Estimating the coefficients of a mixture of two linear regressions by expectation maximization’, IEEE Transactions on Information Theory 65(6), 3515–3524.
  * Kwon et al. (2021) Kwon, J., Ho, N. and Caramanis, C. (2021), On the minimax optimality of the em algorithm for learning two-component mixed linear regression, in ‘International Conference on Artificial Intelligence and Statistics’, PMLR, pp. 1405–1413.
  * Leisch (2004) Leisch, F. (2004), ‘Flexmix: A general framework for finite mixture models and latent class regression in r’, Journal of Statistical Software 11(8), 1–18.
  * Li et al. (2019) Li, Q., Shi, R. and Liang, F. (2019), ‘Drug sensitivity prediction with high-dimensional mixture regression’, PloS one 14(2), e0212108.
  * McLachlan et al. (2019) McLachlan, G. J., Lee, S. X. and Rathnayake, S. I. (2019), ‘Finite mixture models’, Annual review of statistics and its application 6, 355–378.
  * McLachlan and Peel (2004) McLachlan, G. J. and Peel, D. (2004), Finite Mixture Models, John Wiley & Sons.
  * Meinshausen and Yu (2009) Meinshausen, N. and Yu, B. (2009), ‘Lasso-type recovery of sparse representations for high-dimensional data’, The annals of statistics 37(1), 246–270.
  * Negahban et al. (2012) Negahban, S. N., Ravikumar, P., Wainwright, M. J. and Yu, B. (2012), ‘A unified framework for high-dimensional analysis of $m$-estimators with decomposable regularizers’, Statistical science 27(4), 538–557.
  * Song et al. (2014) Song, W., Yao, W. and Xing, Y. (2014), ‘Robust mixture regression model fitting by laplace distribution’, Computational Statistics & Data Analysis 71, 128–137.
  * Städler et al. (2010) Städler, N., Bühlmann, P. and Van De Geer, S. (2010), ‘$\ell$1-penalization for mixture regression models’, Test 19(2), 209–256.
  * Tibshirani (1996) Tibshirani, R. (1996), ‘Regression shrinkage and selection via the lasso’, Journal of the Royal Statistical Society: Series B (Methodological) 58(1), 267–288.
  * Turner (2000) Turner, T. R. (2000), ‘Estimating the propagation rate of a viral infection of potato plants via mixtures of regressions’, Journal of the Royal Statistical Society: Series C (Applied Statistics) 49(3), 371–384.
  * Wang et al. (2014) Wang, Z., Gu, Q., Ning, Y. and Liu, H. (2014), ‘High dimensional expectation-maximization algorithm: Statistical optimization and asymptotic normality’, arXiv preprint arXiv:1412.8729 .
  * Yang and Zou (2015) Yang, Y. and Zou, H. (2015), ‘A fast unified algorithm for solving group-lasso penalize learning problems’, Statistics and Computing 25(6), 1129–1141.
  * Yao et al. (2014) Yao, W., Wei, Y. and Yu, C. (2014), ‘Robust mixture regression using the t-distribution’, Computational Statistics & Data Analysis 71, 116–127.
  * Ye and Zhang (2010) Ye, F. and Zhang, C.-H. (2010), ‘Rate minimaxity of the lasso and dantzig selector for the lq loss in lr balls’, The Journal of Machine Learning Research 11, 3519–3540.
  * Yi and Caramanis (2015) Yi, X. and Caramanis, C. (2015), ‘Regularized em algorithms: A unified framework and statistical guarantees’, Advances in Neural Information Processing Systems 28, 1567–1575.
  * Yi et al. (2014) Yi, X., Caramanis, C. and Sanghavi, S. (2014), Alternating minimization for mixed linear regression, in ‘International Conference on Machine Learning’, pp. 613–621.
  * Yuan and Lin (2006) Yuan, M. and Lin, Y. (2006), ‘Model selection and estimation in regression with grouped variables’, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 49–67.
  * Zhang et al. (2020) Zhang, L., Ma, R., Cai, T. T. and Li, H. (2020), ‘Estimation, confidence intervals, and large-scale hypotheses testing for high-dimensional mixed linear regression’, arXiv preprint arXiv:2011.03598 .
  * Zhu et al. (2017) Zhu, R., Wang, L., Zhai, C. and Gu, Q. (2017), High-dimensional variance-reduced stochastic gradient expectation-maximization algorithm, in ‘International Conference on Machine Learning’, PMLR, pp. 4180–4188.