Unsupervised Liu-type Shrinkage Estimators for Mixture of Regression Models

Elsayed Ghanem†,‡, Armin Hatefi†,111Corresponding author:
Email<EMAIL_ADDRESS>and Tel: +1 (709) 864-8416 and Hamid Usefi†

†Department of Mathematics and Statistics, Memorial University of
Newfoundland, St. John’s, NL, Canada.
‡Faculty of Science, Alexandria University, Arab Republic of Egypt.

Abstract:

In many applications (e.g., medical studies), the population of interest
(e.g., disease status) comprises heterogeneous subpopulations. The mixture of
probabilistic regression models is one of the most common techniques to
incorporate the information of covariates into learning of the population
heterogeneity. Despite its flexibility, the model may lead to unreliable
estimates in the presence of multicollinearity problem. In this paper, we
develop Liu-type shrinkage methods through an unsupervised learning approach
to estimate the model coefficients in multicollinearity. The performance of
the developed methods is evaluated via classification and stochastic versions
of EM algorithms. The numerical studies show that the proposed methods
outperform their Ridge and maximum likelihood counterparts. Finally, the
developed methods are applied to analyze the bone mineral data of women aged
50 and older.

Keywords:  Multicollinearity, Maximum likelihood, Ridge penalty, Liu-type
penalty Mixture models, EM algorithm, Bone mineral data.

##  1  Introduction

As a bone metabolic disease, osteoporosis is characterized when the mineral
density of bone tissues decreases significantly, leading to various major
health problems such as skeletal fragility and osteoporotic fractures. These
can occur in different body areas, including the hip, spine and femur
(Cummings et al., 1995, Neuburger et al., 2015). Bone mineral density (BMD) is
considered one of the most reliable predictors in determining osteoporosis
status (WHO, 1994). For example, approximately every 1 out of 3 women and 1
out of 5 men aged 50 and older have osteoporosis and its fractures (Melton III
et al., 1998). Osteoporosis typically occurs without any significant symptoms.
In South Korea, for example, almost 75% of patients are unaware of their
osteoporosis problem (Lim et al., 2016). BMD score of an individual improves
until age 30 and declines as the individual ages. As the aged population is
growing, it is essential to study osteoporosis to plan the well-being and life
quality of the aged groups of the communities.

BMD values are measured through dual-energy X-ray absorptiometry imaging and
require a costly and time-consuming procedure. Unlike BMD scores, clinicians
have access to various easily attainable explanatory variables about patients,
including BMI, weight, age and test results from previous years. Linear
regression models are well-known statistical tools to investigate the impact
of a set of covariates (e.g., patients’ characteristics) on a response
variable (e.g., BMD measurement). The least squares (LS) method is a common
method for estimating the regression model coefficients; however, the LS
method can lead to extremely unreliable and misleading results in
multicollinearity when the covariates are linearly dependent. As a shrinkage
method, Ridge regression is a common solution to the multicollinearity. A
small ridge parameter is not enough to handle the ill-conditioned design
matrix when the problem is severe. On the other hand, large values of the
ridge parameter render more biases in the estimation. Liu (2003) proposed the
Liu-type (LT) shrinkage method to deal with the challenge and address the
multicollinearity. Duran et al. (2012) extended the LT method to the semi-
parametric regressions. Arashi et al. (2014) proposed the Stein-rule LT
estimators for elliptical regression models. Pearce and Hatefi (2021) used
properties of rank-based samples to improve the LT shrinkage estimators for
linear and logistic regressions.

Finite mixture models (FMMs) are probabilistic model-based tools to analyze
heterogeneous populations. The maximum likelihood (ML) method is one of the
most popular techniques for estimating the parameters of the FMMs. A part of
the popularity comes from the Expectation-Maximization (EM) algorithm
(Dempster et al., 1977) that enables computing the ML estimates of the FMMs.
The EM algorithm decomposes the ML estimation procedure into E- and M-steps
and iteratively estimates the parameters of the FMMs. Stochastic EM (SEM)
algorithm (Celeux, 1985) was developed as a stochastic teacher of the
probabilistic EM algorithm for mixture models. Also, Celeux and Govaert (1992)
proposed a classification EM (CEM) algorithm in which a classification step is
implemented in the iteration for maximization based on classified data. The
mixture of regression models (Quandt and Ramsey, 1978) incorporates the
properties of the FMMs into linear regression models. Jones and McLachlan
(1992), Hawkins et al. (2001) studied the maximum likelihood and EM algorithm
to fit a mixture of regression models. Faria and Soromenho (2010) investigated
various EM algorithms for the mixture of linear regression models. FMMs have
found many applications in the core of statistical sciences, such as times
series data (Zhang et al., 2006), sampling methods and censored data (Hatefi
et al., 2018, Wedel et al., 1998, Hatefi et al., 2015). For more details about
the theory and applications of FMMS, see McLachlan et al. (2019) and
references therein.

Despite the flexibility of the mixture of linear regressions, the ML estimates
of the mixture coefficients become unreliable in multicollinearity. Recently,
Ghanem et al. (2022) proposed shrinkage estimation methods for the mixture
logistic regressions. In the manuscript, we develop the Liu-type shrinkage
estimators for the mixture of linear regression models where the ridge EM
algorithms may not be able to address the ill-conditioned component design
metrics. Through extensive numerical studies, we show that the classification
and stochastic EM algorithms of the LT shrinkage method outperform their
counterparts and provide more reliable coefficient estimates for component
regressions. We finally apply the developed methods to analyze the bone
disorder status of women aged 50 and older.

The outline of the paper is as follows: Section 2 develops the shrinkage
methods and their EM algorithms. Section 3 evaluates the performance of the
methods by simulation studies. The methods are applied for analysis of bone
data in Section 4. Section 5 finally presents the summary and concluding
remarks.

##  2 Statistical Methods

Let ${\bf x}_{i}^{\top}=(x_{i,1},\ldots,x_{i,p})$ be the vector of $p$
explanatory variables for the $i$-th subject in a random sample of size $n$.
Let ${\bf y}=(y_{1},\ldots,y_{n})$ and ${\bf X}=({\bf
x}_{1}^{\top},\ldots,{\bf x}_{n}^{\top})^{\top}$ denote the response vector
and ($n\times p$) design matrix with $\text{rank}({\bf X})=p<n$. The
Regression model $y_{i}={\bf x}_{i}^{\top}{{\bf\beta}}+\epsilon_{i}$ for
$i=1,\ldots,n$ is one of the most common statistical methods to study the
relationship between response variable and the set of the explanatory
variables.

The mixture of regression models is a generalization of the regression model
when the underlying population comprises of heterogeneous subpopulations. The
mixture of regression models is defined by

$\displaystyle y_{i}=\left\\{\begin{array}[]{lc}{\bf
x}_{i}^{\top}{{\bf\beta}}_{1}+\epsilon_{i1}&\text{with probability
$\pi_{1}$}\\\ \vdots\\\ {\bf
x}_{i}^{\top}{{\bf\beta}}_{J}+\epsilon_{iJ}&\text{with probability
$\pi_{J}$}\\\ \end{array}\right.$ (4)

where ${\bf\beta}_{j}=(\beta_{j,1},\ldots,\beta_{J,p})$ represents the
coefficients of $p$ predictors in the $j$-th component regression for
$j=1,\ldots,J$ and ${\pi}=(\pi_{1},\ldots,\pi_{J})$ denotes the vector of
mixing proportions with $\sum_{j=1}^{M}\pi_{j}=1$ and $0<\pi_{j}<1$. Also
$\epsilon_{ij}$ are independent normal random errors from each component of
the mixture; that is $\epsilon_{ij}\overset{iid}{\sim}N(0,\sigma_{j}^{2})$.
Although we assume that the number of components $J$ in the mixture model (4)
is known throughout the paper, the component memberships of the observations
are unknown and should be estimated in an unsupervised approach. Let
${{\bf\beta}}=({\bf\beta}_{1},\ldots,{\bf\beta}_{J})$. Let
${\bf\theta}_{j}=({\bf\beta}_{j},\sigma_{j}^{2})$ denote the parameters of the
$j$-th component. Thus, we represent the vector of all unknown parameters of
the mixture (4) with
${\bf\Psi}=({\bf\pi},{\bf\theta}_{1},\ldots,{\bf\theta}_{J})$.

From regression model (4), the log-likelihood function of ${\bf\Psi}$ can be
written as

$\displaystyle\ell({\bf\Psi})=\sum_{i=1}^{n}\log\left(\sum_{j=1}^{J}\pi_{j}\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})\right),$ (5)

where $\phi_{j}({\bf x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})$ represents
the pdf of the univariate normal distribution with mean ${\bf
x}_{i}^{\top}{{\bf\beta}}_{j}$ and variance $\sigma_{j}^{2}$. We have to
maximize (5) to obtain the ML estimate of ${\bf\Psi}$. The gradient of (5) is
not tractable with respect to component parameters
${\bf\theta}_{j},j=1,\ldots,J$. We view $({\bf X},{\bf y})$ as incomplete data
and apply the expectation-maximization (EM) algorithm of (Dempster et al.,
1977) using a complete data to find $\widehat{{\bf\Psi}}_{ML}$. For each
subject $({\bf x}_{i},y_{i})$, we introduce latent variables ${\bf
Z}_{i}=(Z_{i1},\ldots,Z_{iJ})$ for $i=1,\ldots,n$ as

$\displaystyle Z_{ij}=\left\\{\begin{array}[]{lc}1&\text{if the $i$-th subject
comes from the $j$-th component,}\\\ 0&o.w.,\end{array}\right.$

where ${\bf
Z}_{i}\mathrel{\overset{iid}{\scalebox{1.5}[1.0]{$\sim$}}}\text{Multi}(1,\pi_{1},\ldots,\pi_{M})$.
From the marginal distribution of latent variables, the conditional
distribution of ${\bf Z}_{i}|y_{i}$ is given by

$\displaystyle f({\bf
z}_{i}|y_{i})=\prod_{j=1}^{J}\left\\{\frac{\pi_{j}\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})}{\sum_{j=1}^{J}\pi_{j}\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})}\right\\}^{z_{ij}}.$ (6)

From above, it is easily seen that ${\bf
Z}_{i}|y_{i}\mathrel{\overset{iid}{\scalebox{1.5}[1.0]{$\sim$}}}\text{Multi}(1,\tau_{i1}({\bf\Psi}),\ldots,\tau_{iJ}({\bf\Psi}))$
where

$\displaystyle\tau_{ij}({\bf\Psi})=\frac{\pi_{j}\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})}{\sum_{j=1}^{J}\pi_{j}\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})}.$ (7)

Let $({\bf X},{\bf y},{\bf Z})$ denote the complete data. Thus, the complete
log-likelihood function of ${\bf\Psi}$ is given by

$\displaystyle\ell_{c}({\bf\Psi})=\sum_{i=1}^{n}\sum_{j=1}^{J}z_{ij}\log(\pi_{j})+\sum_{i=1}^{n}\sum_{j=1}^{J}z_{ij}\log\\{\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})\\}.$ (8)

###  2.1 ML Estimation Method

The EM algorithm is a standard method to find the ML estimates of the mixture
model parameters. The EM algorithm employs the latent variables on top of the
observed data and decomposes the estimation procedure into an iterative
expectation (E) and maximization (M) steps. As an iterative method, EM
algorithm begins with an initial value. Let
${\bf\Psi}^{(0)}=({\bf\pi}^{(0)},{\bf\theta}_{1}^{(0)},\ldots,{\bf\theta}_{J}^{(0)})$
and ${\bf\Psi}^{(r)}$ denote the initial vector and the estimate in the $r$-th
iteration of the EM algorithm, respectively.

On the $(r+1)$-th iteration, we require to compute the conditional expectation
of (8) in the E-step. The $Q({\bf\Psi},{\bf\Psi}^{(r)})$ replaces the latent
variables by their conditional expectations as

${\bf Q}({\bf\Psi},{\bf\Psi}^{(r)})={\bf Q}_{1}({\bf\pi},{\bf\Psi}^{(r)})+{\bf
Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)}),$

where

$\displaystyle{\bf
Q}_{1}({\bf\pi},{\bf\Psi}^{(r)})=\sum_{i=1}^{n}\sum_{j=1}^{J}\tau_{ij}({\bf\Psi}^{(r)})\log(\pi_{j}),$
(9)

and

$\displaystyle{\bf
Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})=\sum_{i=1}^{n}\sum_{j=1}^{J}\tau_{ij}({\bf\Psi}^{(r)})\log\left\\{\phi_{j}({\bf
x}_{i}^{\top}{{\bf\beta}}_{j},\sigma_{j}^{2})\right\\},$ (10)

where $\tau_{ij}({\bf\Psi}^{(r)})$ is obtained by (7). In the M-step, we
maximize ${\bf Q}({\bf\Psi},{\bf\Psi}^{(r)})$ with respect to ${\bf\pi}$ and
${{\bf\theta}}$ to update ${\bf\Psi}^{(r+1)}$. One can update
${\bf\pi}^{(r+1)}$ by maximizing ${\bf Q}_{1}({\bf\pi},{\bf\Psi}^{(r)})$
subject to $\sum_{j=1}^{J}\pi_{j}=1$ as follows

$\displaystyle{\widehat{\pi}}_{j}^{(r+1)}=\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})/n;~{}~{}~{}~{}j=1,\ldots,J-1.$
(11)

The maximization of ${\bf Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})$ can be
reformulated by the weighted least squares (WLS) method as

$\displaystyle{\widehat{\bf\beta}}_{j}^{(r+1)}=\underset{{\bf\beta}_{j}}{\arg\min}~{}({\bf
y}-{\bf X}{\bf\beta})^{\top}{\bf W}_{j}({\bf y}-{\bf X}{\bf\beta})/n,$ (12)

where ${\bf W}_{j}$ is $n\times n$ diagonal matrix with diagonal elements
$(\tau_{ij}({\bf\Psi}^{(r)}),\ldots,\tau_{nj}({\bf\Psi}^{(r)}))$ for all
$j=1,\ldots,J$. One can easily update ${\widehat{\bf\beta}}_{j}^{(r+1)}$ as
the solution to (12) by

$\displaystyle{\widehat{\bf\beta}}_{j}^{(r+1)}=\left({\bf X}^{\top}{\bf
W}_{j}{\bf X}\right)^{-1}{\bf X}^{\top}{\bf W}_{j}{\bf
y},~{}~{}~{}j=1,\ldots,J.$ (13)

From (12) and following (Faria and Soromenho, 2010), we then update
${\widehat{\sigma}}_{j}^{2(r+1)}$ as follows

$\displaystyle{\widehat{\sigma}}_{j}^{2(r+1)}=\frac{({\bf y}-{\bf
X}\widehat{{\bf\beta}}^{(r+1)})^{\top}{\bf W}^{(r)}_{j}({\bf y}-{\bf
X}\widehat{{\bf\beta}}^{(r+1)})}{\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})},~{}~{}~{}j=1,\ldots,J.$
(14)

To find ${\widehat{\bf\Psi}}_{ML}$, we iteratively alternate the E- and M-
steps of the EM algorithm until the stopping criterion
$|\ell({\bf\Psi}^{(r+1)})-\ell({\bf\Psi}^{(r)})|$ becomes negligible.

Classification EM Algorithm: In the above EM algorithm, we use information
from all observations (as membership probabilities) in each iteration to
estimate the parameters of the mixture model. Following Celeux and Govaert
(1992), we shall employ the classification version of EM algorithm (CEM) to
estimate ${\bf\Psi}$. The CEM algorithm incorporates a classification (C) step
between E and M steps so that the component parameters of the mixture are
updated in M-step using the classified complete data log-likelihood function.

The E-step here is identical to the E-step of the EM algorithm. In C-step, the
observations are then assigned to $J$ mutually exclusive partitions
corresponding to the $J$ components of mixture model (4). Let ${\bf
P}^{(r+1)}=(P_{1}^{(r+1)},\ldots,P_{J}^{(r+1)})$ denote the partition in the
$(r+1)$-th iteration. Each subject $({\bf x}_{i},y_{i})$ is assigned to
partition $P_{h}^{(r+1)}$ when

$\tau_{ih}({\bf\Psi}^{(r)})=\underset{j}{\arg\min}~{}\tau_{ij}({\bf\Psi}^{(r)}).$

Note when the maximum weight is not unique, the tie is broken at random. Also
the CEM algorithm is stopped and ${\bf\Psi}^{(r)}$ is returned when a
partition becomes empty or has only one observation.

In the M-step, we maximize conditional expectation function using ${\bf
P}^{(r+1)}$. From (9), the mixing proportion is updated by

$\displaystyle{\widehat{\pi}}_{j}^{(r+1)}=n_{j}/n;~{}~{}~{}~{}j=1,\ldots,J,$
(15)

where $n_{j}$ is the number of observations allocated to partition
$P_{j}^{(r+1)}$. Applying the WLS to each partition $P_{j}^{(r+1)}$, we can
update the parameters of the $j$-th component ($j=1,\ldots,J$) by

$\displaystyle{\widehat{\bf\beta}}_{j}^{(r+1)}=\left({\bf X}_{j}^{\top}{\bf
W}_{j}{\bf X}_{j}\right)^{-1}{\bf X}_{j}^{\top}{\bf W}_{j}{\bf y}_{j},$ (16)
$\displaystyle{\widehat{\sigma}}_{j}^{2(r+1)}=\frac{({\bf y}_{j}-{\bf
X}_{j}\widehat{{\bf\beta}}^{(r+1)})^{\top}{\bf W}^{(r)}_{j}({\bf y}_{j}-{\bf
X}_{j}\widehat{{\bf\beta}}^{(r+1)})}{\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})},$
(17)

where ${\bf X}_{j}$ and ${\bf y}_{j}$ represent, respectively, $(n_{j}\times
p)$ design matrix and vector of responses corresponding to $P_{j}^{(r+1)}$.
Also ${\bf W}^{(r)}_{j}$ is the diagonal weight matrix of size $n_{j}$ with
entries $(\tau_{ij}({\bf\Psi}^{(r)}),\ldots,\tau_{n_{j},j}({\bf\Psi}^{(r)}))$.
Finally, we alternate repeatedly the E-, C- and M- steps until
$|\ell({\bf\Psi}^{(r+1)})-\ell({\bf\Psi}^{(r)})|<{\epsilon}$.

Stochastic EM Algorithm: One can also apply the stochastic version of the EM
algorithm (Celeux, 1985) to fit the mixture of regression models. The
stochastic EM (SEM) algorithm implements a stochastic version of the C-step
between E- and M-steps. Although the E- and M-steps are identical to the CEM
algorithm, the SEM simulates the component membership of each observation
using conditional distribution of the latent variable given incomplete data.
On the $(r+1)$-th iteration, the S-step simulates a random allocation for each
observation using one draw out of $J$ components as

${\bf
Z}_{i}^{*}=(Z_{i1}^{*},\ldots,Z_{iJ}^{*})\mathrel{\overset{iid}{\scalebox{1.5}[1.0]{$\sim$}}}\text{Multi}(1,\tau_{i1}({\bf\Psi}^{(r)}),\ldots,\tau_{iJ}({\bf\Psi}^{(r)}))~{}~{}(i=1,\ldots,n).$

Then the $({\bf x}_{i},y_{i})$ is classified to partition $P_{j}^{(r+1)}$ if
$Z_{ij}^{*}=1;i=1,\ldots,j=1,\ldots,J$. Using the stochastic partitions, we
update the mixture parameters from (15), (16) and (17) in M-step. From
(Celeux, 1985, Faria and Soromenho, 2010), point-wise convergence in the SEM
algorithm is not guaranteed. The algorithm resembles a Markov chain where at
stationary state fluctuates around the ML estimate. Hence, we alternate the
E-, S- and M-steps until either the criterion
$|\ell({\bf\Psi}^{(r+1)})-\ell({\bf\Psi}^{(r)})|<\epsilon$ satisfies or the
chain reaches a pre-specified maximum number of iterations that is fixed for
all the algorithms for fair comparison.

###  2.2 Ridge Estimation Method

The ML method is a standard tool to estimate ${{\bf\Psi}}$; however, the ML
estimates are dramatically influenced by multicollinearity where the
covariates are linearly dependent. Ridge method (Hoerl and Kennard, 1970) is
one of the most common methods to encounter with the challenges of the LS
method. The ridge estimate for (4) can be obtained as a solution to the
penalized log-likelihood function given by

$\displaystyle\ell^{R}({\bf\Psi})=\ell({\bf\Psi})-k{\bf\beta}^{\top}{\bf\beta}/2$
(18)

where $\ell({\bf\Psi})$ comes from (5) and $k>0$ is the ridge parameter. In a
similar vein to Subsection 2.1, for each observation $({\bf
x}_{i},y_{i}),i=1,\ldots,n$, we first introduce $J$ dimensional latent vector
${\bf Z}_{i}=(Z_{i1},\ldots,Z_{iJ})$. We then develop an EM algorithm to
maximize the complete ridge log-likelihood function and obtain
${\widehat{\bf\Psi}}_{R}$.

The E-step of the ridge EM algorithm is identical to the E-step of Subsection
2.1. In the M-step, the mixing proportion are updated from (11). To update the
coefficients of the component regressions, we require to maximize ${\bf
Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})$ subject to the ridge penalty within
each component as

${\bf Q}^{R}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})={\bf
Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})-k_{j}{\bf\beta}_{j}^{\top}{\bf\beta}_{j}/2,$

where ${\bf Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})$ is from (10) and $k_{j}$ is
the ridge parameter in the $j$-th component. Like ML method, one can re-write
the maximization of ${\bf Q}^{R}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})$ as a WLS
subject to the ridge penalty by

$\displaystyle{\widehat{\bf\beta}}_{R,j}^{(r+1)}=\underset{{\bf\beta}_{j}}{\arg\min}~{}({\bf
y}-{\bf X}{\bf\beta})^{\top}{\bf W}_{j}({\bf y}-{\bf
X}{\bf\beta})+k_{j}{\bf\beta}_{j}^{\top}{\bf\beta}_{j}/2,$ (19)

where ${\bf W}_{j}$ is $n\times n$ diagonal matrix with elements
$(\tau_{ij}({\bf\Psi}^{(r)}),\ldots,\tau_{nj}({\bf\Psi}^{(r)}))$ obtained from
(7). Applying (19), ${\widehat{\bf\beta}}_{j}^{(r+1)};j=1,\ldots,J$ is updated
by

$\displaystyle{\widehat{\bf\beta}}_{R,j}^{(r+1)}=\left({\bf X}^{\top}{\bf
W}_{j}{\bf X}+k_{j}\mathbb{I}\right)^{-1}{\bf X}^{\top}{\bf W}_{j}{\bf y}.$
(20)

######  Lemma 1.

Under the assumptions of mixture of regression models (4), suppose
$\lambda_{1j},\ldots,\lambda_{pj}$ and $u_{1j},\ldots,u_{pj}$ be eigenvalues
and orthonormal eigenvectors of ${\bf X}^{\top}{\bf W}_{j}{\bf X}$ where ${\bf
W}_{j}$ is $n\times n$ diagonal matrix with entries
$(\tau_{ij}({\bf\Psi}^{(r)}),\ldots,\tau_{nj}({\bf\Psi}^{(r)}))$ under ridge
EM algorithm. Let
${{\bf\Lambda}}_{j}=\text{diag}(\lambda_{1j},\ldots,\lambda_{pj})$ and ${\bf
U}_{j}=[u_{1j},\ldots,u_{pj}]$. Then the canonical weighted ridge estimator in
each component regression is given by

${\widehat{\bf\alpha}}_{R,j}=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}{{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y},$

and

${\widehat{\bf\beta}}_{R,j}={\bf U}_{j}{\widehat{\bf\alpha}}_{R,j},$

with ${\bf V}_{1}=[v_{1j},\ldots,v_{pj}]$ where $v_{1j},\ldots,v_{pj}$ are the
orthonormal eigenvectors of ${\bf W}_{j}^{1/2}{\bf X}{\bf X}^{\top}{\bf
W}_{j}^{1/2}$.

From (19), the variance term can be updated by

$\displaystyle{\widehat{\sigma}}_{R,j}^{2(r+1)}=\frac{({\bf y}-{\bf
X}\widehat{{\bf\beta}}_{R}^{(r+1)})^{\top}{\bf W}^{(r)}_{j}({\bf y}-{\bf
X}\widehat{{\bf\beta}}_{R}^{(r+1)})}{\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})},$
(21)

where
$\widehat{{\bf\beta}}_{R}^{(r+1)}=(\widehat{{\bf\beta}}_{R,1}^{(r+1)},\ldots,\widehat{{\bf\beta}}_{R,J}^{(r+1)})$.
There are various methods available in the literature for the estimation of
$k_{j}$. Following Hoerl et al. (1975), Liu (2003), we estimate the parameter
by
$\widehat{k}_{j}=p{\widehat{\sigma}}_{ML,j}^{2}/{\widehat{\bf\beta}}_{ML,j}^{\top}{\widehat{\bf\beta}}_{ML,j}$
where ${\widehat{\sigma}}_{ML,j}^{2}$ and ${\widehat{\bf\beta}}_{ML,j}$ are
calculated from (14) and (13), respectively. The E- and M-steps are repeatedly
computed until
$|\ell^{R}({\bf\Psi}^{(r+1)})-\ell^{R}({\bf\Psi}^{(r)})|<\epsilon$.

Ridge CEM Algorithm: One can also apply the CEM algorithm to find
${\widehat{\bf\Psi}}_{R}$. Similar to the CEM algorithm in Subsection 2.1, we
require to accommodate a C-step between E- and M-steps in the ridge EM
algorithm. Here the E-step remains the same as before. Similar to the C-step
of ML method, we classify the observations to partitions ${\bf
P}^{(r+1)}=(P_{1}^{(r+1)},\ldots,P_{J}^{(r+1)})$ based on the maximum
probability of memberships; that is

$P_{j}^{(r+1)}=\\{({\bf
x}_{i},y_{i});\tau_{ij}({\bf\Psi}^{(r)})=\underset{h}{\arg\max}~{}\tau_{ih}({\bf\Psi}^{(r)})\\},~{}~{}\forall
j=1,\ldots,J.$

Based on ${\bf P}^{(r+1)}$, we use (15) to update the mixing proportions of
the mixture. The ridge parameters are estimated similar to the ridge EM
algorithm. We apply (19) to each partition $P_{j}^{(r+1)}$ and update the
coefficients and variance term of each component regression by

$\displaystyle{\widehat{\bf\beta}}_{R,j}^{(r+1)}=\left({\bf X}_{j}^{\top}{\bf
W}_{j}{\bf X}_{j}+k_{j}\mathbb{I}\right)^{-1}{\bf X}_{j}^{\top}{\bf W}_{j}{\bf
y}_{j},$ (22) $\displaystyle{\widehat{\sigma}}_{R,j}^{2(r+1)}=\frac{({\bf
y}_{j}-{\bf X}_{j}\widehat{{\bf\beta}}_{R}^{(r+1)})^{\top}{\bf
W}^{(r)}_{j}({\bf y}_{j}-{\bf
X}_{j}\widehat{{\bf\beta}}_{R}^{(r+1)})}{\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})},$
(23)

where ${\bf X}_{j}$ is $(n_{j}\times p)$ design matrix and ${\bf y}_{j}$ is
vector of responses from observations classified to $P_{j}^{(r+1)}$. ${\bf
W}^{(r)}_{j}$ is the diagonal weight matrix with entries
$\left(\tau_{ij}({\bf\Psi}^{(r)}),\ldots,\tau_{n_{j},j}({\bf\Psi}^{(r)})\right)$
from (7). Finally, the E-, C- and M- steps under ridge estimation procedure
are alternated until convergence criterion is statified.

Ridge SEM Algorithm: Ridge estimation method can be implemented via the
stochastic EM algorithm. Like SEM of ML method, the S-step determines
stochastically the component membership of observations under ridge method by
${\bf
Z}_{i}^{*}=(Z_{i1}^{*},\ldots,Z_{iJ}^{*})\mathrel{\overset{iid}{\scalebox{1.5}[1.0]{$\sim$}}}\text{Multi}(1,\tau_{i1}({\bf\Psi}^{(r)}),\ldots,\tau_{iJ}({\bf\Psi}^{(r)}));i=1,\ldots,n$
and updates ${\bf P}^{(r+1)}=(P_{1}^{(r+1)},\ldots,P_{J}^{(r+1)})$ such that
$P_{j}^{(r+1)}=\\{({\bf x}_{i},y_{i});Z_{ij}^{*}=1\\};\forall j=1,\ldots,J$.
Based on this stochastic partition of S-step, we update the mixture parameters
by (15), (22) and (23).

######  Lemma 2.

Under the assumptions of mixture of regression models (4), with component
regression models ${\bf y}_{j}={\bf X}_{j}{\bf\beta}_{j}+{\epsilon}$ based on
$n_{j}$ observations with $\text{rank}({\bf X}_{j})=p$. Suppose
$\lambda_{1j},\ldots,\lambda_{pj}$ and $u_{1j},\ldots,u_{pj}$ be eigenvalues
and orthonormal eigenvectors of ${\bf X}_{j}^{\top}{\bf W}_{j}{\bf X}_{j}$
where ${\bf W}_{j}$ is $n_{j}\times n_{j}$ diagonal matrix with entries
$(\tau_{ij}({\bf\Psi}^{(r)}),\ldots,\tau_{nj}({\bf\Psi}^{(r)}))$ under the
ridge CEM or ridge SEM algorithm. Let
${{\bf\Lambda}}_{j}=\text{diag}(\lambda_{1j},\ldots,\lambda_{pj})$ and ${\bf
U}_{j}=[u_{1j},\ldots,u_{pj}]$. Then The canonical weighted ridge estimator in
each component regression is given by

${\widehat{\bf\alpha}}_{R,j}=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}{{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y},$

and

${\widehat{\bf\beta}}_{R,j}={\bf U}_{j}{\widehat{\bf\alpha}}_{R,j},$

with ${\bf V}_{1}=[v_{1j},\ldots,v_{pj}]$ where $v_{1j},\ldots,v_{pj}$ are the
orthonormal eigenvectors of ${\bf W}_{j}^{1/2}{\bf X}_{j}{\bf
X}_{j}^{\top}{\bf W}_{j}^{1/2}$.

Finally, the E-, S- and M-steps are iterated until the stopping rule is
satisfied or the algorithm reaches a pre-specified maximum number of
iterations.

###  2.3 Liu-type Estimation Method

When the design matrix is severely ill-conditioned, adding small values to the
diagonal elements by the ridge estimator may be unable to cope with the
problem. On the other side, increasing the ridge parameter may result in a
more considerable bias in the ridge estimation method. Liu (2003) proposed
Liu-type (LT) shrinkage method for the multicollinearity in estimating the
regression parameters. Like the ridge method, the LT method optimizes the
estimating equation subject to the LT penalty to control the
multicollinearity. The LT penalty is given by

$\displaystyle(-\frac{d}{k^{1/2}}){\widehat{\bf\beta}}=k^{1/2}{\bf\beta}+\epsilon^{\prime},$
(24)

where ${\widehat{\bf\beta}}$ can be any estimator of coefficients and
$d\in\mathbb{R}$ and $\lambda>0$ are two parameters of the LT method. We
develop the LT shrinkage method in estimating the unknown parameters of the
mixture model (4). We shall find the LT estimate of ${\bf\Psi}$ by maximizing
the log-likelihood function (5) subject to the LT penalty.

In a similar vein to Subsection 2.1, the penalized log-likelihood function
based on observed data is not tractable with respect to the component
parameters. We develop the LT estimation procedure through an unsupervised
approach and use the EM algorithm to estimate ${{\bf\Psi}}$. We again require
latent vectors ${\bf Z}_{i}=(Z_{i1},\ldots,Z_{iJ})$ to represent the component
membership of the $i$-th observation $({\bf x}_{i},y_{i});i=1,\ldots,n$. Let
$({\bf X},{\bf y},{\bf Z})$ denote the complete data. Then EM algorithm under
the LT method proceeds as follows.

In $(r+1)$-th iteration, the E-step remains identical to the E-step of the ML
method. The mixing proportions under the LT method are updated by (11). We
require to maximize ${\bf Q}_{2}({{\bf\theta}},{\bf\Psi}^{(r)})$ from (10)
under the LT penalty within each component to estimate the component
parameters. The LT penalized log-likelihood function can be written as a WLS
constrained on the LT penalty as

$\displaystyle{\widehat{\bf\beta}}_{LT,j}^{(r+1)}=\underset{{\bf\beta}_{j}}{\arg\min}~{}({\bf
y}-{\bf X}{\bf\beta})^{\top}{\bf W}_{j}({\bf y}-{\bf
X}{\bf\beta})+\left[(-\frac{d_{j}}{k_{j}^{1/2}}){\widehat{\bf\beta}}_{j}-k_{j}^{1/2}{\bf\beta}_{j}\right]^{\top}\left[(-\frac{d_{j}}{k_{j}^{1/2}}){\widehat{\bf\beta}}_{j}-k_{j}^{1/2}{\bf\beta}_{j}\right],$
(25)

where ${\widehat{\bf\beta}}_{j}$ can be any estimate and ${\bf W}_{j}$ is a
weight diagonal matrix with elements
$(\tau_{1j}({\bf\Psi}^{(r)}),\ldots,\tau_{nj}({\bf\Psi}^{(r)});j=1,\ldots,J$.
From (25), the coefficients and variance terms in each component regression
are updated by

$\displaystyle{\widehat{\bf\beta}}_{LT,j}^{(r+1)}=\left({\bf X}^{\top}{\bf
W}_{j}{\bf X}+k_{j}\mathbb{I}\right)^{-1}({\bf X}^{\top}{\bf W}_{j}{\bf
y}-d_{j}{\widehat{\bf\beta}}_{j}),$ (26)
$\displaystyle{\widehat{\sigma}}_{LT,j}^{2(r+1)}=\frac{({\bf y}-{\bf
X}\widehat{{\bf\beta}}_{LT}^{(r+1)})^{\top}{\bf W}^{(r)}_{j}({\bf y}-{\bf
X}\widehat{{\bf\beta}}_{LT}^{(r+1)})}{\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})},$
(27)

where
$\widehat{{\bf\beta}}_{LT}^{(r+1)}=(\widehat{{\bf\beta}}_{LT,1}^{(r+1)},\ldots,\widehat{{\bf\beta}}_{LT,J}^{(r+1)})$.
The (26) and (27) require the estimates of the LT parameters $(k_{j},d_{j})$
for each component regression. From Liu (2003), we can estimate $k_{j}$ in the
$j$-th component by ${\widehat{k}}_{LT,j}=\lambda_{1,j}-100\lambda_{p,j}/99$
where $\lambda_{1,j}$ and $\lambda_{p,j}$ are the maximum and the minimum
eigenvalues of ${\bf X}^{\top}{\bf W}_{j}{\bf X}$ on the $(r+1)$-the iteration
of the EM algorithm.

######  Lemma 3.

Under the assumptions of Lemma (1), the canonical LT estimator in the $j$-th
component regression $j=1,\ldots,J$ under EM algorithm is given by

${\widehat{\bf\alpha}}_{LT,j}=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}({{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}-d_{j}{\widehat{\bf\alpha}}_{j}),$

and

${\widehat{\bf\beta}}_{LT,j}={\bf U}_{j}{\widehat{\bf\alpha}}_{LT,j},$

where ${\widehat{\bf\alpha}}_{j}$ is the canonical estimate of
${\bf\beta}_{j}$ and ${\bf V}_{1}=[v_{1j},\ldots,v_{pj}]$ with
$v_{1j},\ldots,v_{pj}$ are the orthonormal eigenvectors of ${\bf
W}_{j}^{1/2}{\bf X}{\bf X}^{\top}{\bf W}_{j}^{1/2}$.

Following Liu (2003) and Lemma 3, the optimal $d_{j}$ can be obtained by the
next lemma within each component of the mixture of regression models.

######  Lemma 4.

Under the assumptions of Lemma 1 for $k_{j}>0$

  * i)

when ${\widehat{\bf\alpha}}_{j}={\widehat{\bf\alpha}}_{ML,j}$, then

$d_{j}=\sum_{m=1}^{p}\left((\sigma_{j}^{2}-k_{j}\alpha_{mj}^{2})/(\lambda_{mj}+k_{j})^{2}\right)/\sum_{m=1}^{p}\left((\lambda_{mj}\alpha_{mj}^{2}+\sigma_{j}^{2})/\lambda_{mj}(\lambda_{mj}+k_{j})^{2}\right),$

  * ii)

when ${\widehat{\bf\alpha}}_{j}={\widehat{\bf\alpha}}_{R,j}$, then

$d_{j}=\sum_{m=1}^{p}\left(\lambda_{mj}(\sigma_{j}^{2}-k_{j}\alpha_{mj}^{2})/(\lambda_{mj}+k_{j})^{3}\right)/\sum_{m=1}^{p}\left(\lambda_{mj}(\lambda_{mj}\alpha_{mj}^{2}+\sigma_{j}^{2})/(\lambda_{mj}+k_{j})^{4}\right),$

minimizes the $\text{MSE}({\widehat{\bf\alpha}}_{LT,j})$ within each component
(4) in the EM algorithm of the LT method.

Although Lemma 4 paves the path in estimating the optimal LT parameter $d_{j}$
within each component regression, the optimal value still depends on the
unknown quantities including $\sigma_{j}$, $k_{j}$, ${\bf\alpha}_{j}$ and
$\lambda_{m,j}$ for $m=1,\ldots,p$ and $j=1,\ldots,J$. From Lemma 4, we
propose a practical approach where $d_{j},j=1,\ldots,J$ can be updated in the
$(r+1)$-the iteration of the EM algorithm by

$\displaystyle{\widehat{d}}_{j}=\sum_{m=1}^{p}\left(\lambda_{mj}({\widehat{\sigma}}_{R,j}^{2}-{\widehat{k}}_{j}{\widehat{\alpha}}_{R,mj}^{2})/(\lambda_{mj}+{\widehat{k}}_{j})^{3}\right)/\sum_{m=1}^{p}\left(\lambda_{mj}(\lambda_{mj}{\widehat{\alpha}}_{R,mj}^{2}+{\widehat{\sigma}}_{R,j}^{4})/\lambda_{mj}(\lambda_{mj}+{\widehat{k}}_{j})^{2}\right),$
(28)

where ${\widehat{k}}_{j}={\widehat{k}}_{LT,j}$,
${\widehat{\bf\alpha}}_{R,j}=({\widehat{\alpha}}_{R,1j},\ldots,{\widehat{\alpha}}_{R,pj})$
is given by Lemma 1 and $(\lambda_{1j},\ldots,\lambda_{pj})$ are eigenvalues
of ${\bf X}_{j}^{\top}{\bf W}_{j}{\bf X}_{j}$ with
${\widehat{\sigma}}_{R,j}^{2}$ from (21). The E- and M-steps are alternated
until the stopping criterion is satisfied. As the parameters $k_{j}$ and
$d_{j}$ are updated in each iteration of the EM algorithm, the proposed LT
method is henceforth called the iterative Liu-type i.e., LT(ITR).

Unlike the iterative LT method, one can follow Hoerl et al. (1975) to estimate
the LT parameters based on ridge estimates ${\widehat{\sigma}}_{R,j}$ and
${\widehat{\bf\beta}}_{R,j}$. In other words, mixture parameter ${\bf\Psi}$
and $d_{j}$ are still iteratively updated in the EM algorithm; however,
$k_{j},j=1\ldots,J$ parameters are estimated only once throughout the EM
algorithm using the ridge estimates. Here, we estimate the parameters by
${\widehat{k}}_{LT,j}=p{\widehat{\sigma}}_{R,j}/{\widehat{\bf\beta}}_{R,j}^{\top}{\widehat{\bf\beta}}_{R,j}$.
This LT estimation method is henceforth is called HKP Liu-type i.e., LT(HKP).

Liu-type CEM Algorithm: Like previous subsections, the CEM algorithm
partitions the observations in the C-step and then update the parameters with
in each partition. In the $(r+1)$-th iteration of the CEM algorithm, the
E-step remains the same as before. The C-step classifies the observations into
partition ${\bf P}^{(r+1)}=(P_{1}^{(r+1)},\ldots,P_{J}^{(r+1)})$ where
$P_{j}^{(r+1)}=\\{({\bf
x}_{i},y_{i});\tau_{ij}({\bf\Psi}^{(r)})=\underset{h}{\arg\max}~{}\tau_{ih}({\bf\Psi}^{(r)})\\}$
with $(\tau_{i1}({\bf\Psi}^{(r)},\ldots,\tau_{iJ}({\bf\Psi}^{(r)})$ are
obtained from (7). Using ${\bf P}^{(r+1)}$, we update the mixing proportions
from (15). We then require to estimate the LT parameters $(k_{,}d_{j})$ in
each iteration of the CEM algorithm, Like Liu-type EM algorithm, we propose
${\widehat{k}}_{LT,j}=\lambda_{1,j}-100\lambda_{p,j}/99$ where $\lambda_{1,j}$
and $\lambda_{p,j}$ are the maximum and the minimum eigenvalues of ${\bf
X}_{j}^{\top}{\bf W}_{j}{\bf X}_{j}$.

Figure 1: The median (M), lower (L) and upper (U) bounds of 95% CIs for
$\text{SSE}(\widehat{\bf\beta})$ of the estimators when the population is a
mixture of two regression models with $n=60$.

######  Lemma 5.

Under the assumptions of Lemma (1), the canonical LT estimator in the $j$-th
component regression $j=1,\ldots,J$ under the CEM algorithm is given by

${\widehat{\bf\alpha}}_{LT,j}=\left({{\bf\Lambda}}_{j}+k_{j}\right)^{-1}({{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}-d_{j}{\widehat{\bf\alpha}}_{j}),$

and

${\widehat{\bf\beta}}_{LT,j}={\bf U}_{j}{\widehat{\bf\alpha}}_{LT,j},$

where ${\widehat{\bf\alpha}}_{j}$ is the canonical estimate of
${\bf\beta}_{j}$ and ${\bf V}_{1}=[v_{1j},\ldots,v_{pj}]$ with
$v_{1j},\ldots,v_{pj}$ are the orthonormal eigenvectors of ${\bf
W}_{j}^{1/2}{\bf X}_{j}{\bf X}_{j}^{\top}{\bf W}_{j}^{1/2}$.

From Lemma 5 and Lemma 4, one can estimate parameter $d_{j}$ based on
partition $P_{j}^{(r+1)}$ from (28) where $(\lambda_{1j},\ldots,\lambda_{pj})$
are eigenvalues of ${\bf X}_{j}^{\top}{\bf W}_{j}{\bf X}_{j}$ and
${\widehat{\sigma}}_{R,j}^{2}$ from (23). To estimate the regression
parameters, we implement a WLS based on the LT penalty as

$\displaystyle{\widehat{\bf\beta}}_{LT,j}^{(r+1)}=\underset{{\bf\beta}_{j}}{\arg\min}~{}({\bf
y}_{j}-{\bf X}_{j}{\bf\beta})^{\top}{\bf W}_{j}({\bf y}_{j}-{\bf
X}_{j}{\bf\beta})+\left[(-\frac{\widehat{d}_{j}}{{\widehat{k}}_{LT,j}^{1/2}}){\widehat{\bf\beta}}_{j}-{\widehat{k}}_{LT,j}^{1/2}{\bf\beta}_{j}\right]^{\top}\left[(-\frac{\widehat{d}_{j}}{{\widehat{k}}_{LT,j}^{1/2}}){\widehat{\bf\beta}}_{j}-{\widehat{k}}_{LT,j}^{1/2}{\bf\beta}_{j}\right],$
(29)

where ${\bf y}_{j}$ and ${\bf X}_{j}$ are response vector and design matrix
under $P_{j}^{(r+1)}$ and ${\widehat{\bf\beta}}_{j}$ can be any estimate for
${{\bf\beta}}_{j}$. Also, ${\bf W}_{j}$ is a weight diagonal matrix with
entries
$(\tau_{1j}({\bf\Psi}^{(r)}),\ldots,\tau_{n_{j},j}({\bf\Psi}^{(r)});j=1,\ldots,J$.
One can easily find the solution to (25) and update the regression parameters
by

$\displaystyle{\widehat{\bf\beta}}_{LT,j}^{(r+1)}=\left({\bf X}_{j}^{\top}{\bf
W}_{j}{\bf X}_{j}+{\widehat{k}}_{LT,j}\mathbb{I}\right)^{-1}({\bf
X}_{j}^{\top}{\bf W}_{j}{\bf y}_{j}-\widehat{d}_{j}{\widehat{\bf\beta}}_{j}),$
(30) $\displaystyle{\widehat{\sigma}}_{LT,j}^{2(r+1)}=\frac{({\bf y}_{j}-{\bf
X}_{j}\widehat{{\bf\beta}}_{LT}^{(r+1)})^{\top}{\bf W}^{(r)}_{j}({\bf
y}_{j}-{\bf
X}_{j}\widehat{{\bf\beta}}_{LT}^{(r+1)})}{\sum_{i=1}^{n}\tau_{ij}({\bf\Psi}^{(r)})},$
(31)

with
$\widehat{{\bf\beta}}_{LT}^{(r+1)}=(\widehat{{\bf\beta}}_{LT,1}^{(r+1)},\ldots,\widehat{{\bf\beta}}_{LT,J}^{(r+1)})$.
The E-, C- and M-steps are repeatedly computed until the convergence criterion
is satisfied.

Figure 2: The median (M), lower (L) and upper (U) bounds of 95% CIs for
$\text{SSE}(\widehat{\pi})$ of the estimators when the population is a mixture
of two regression models with $n=60$.

Unlike the iterative Liu-type CEM algorithm, one may estimate the parameters
of the mixture model via the HKP Liu-type CEM algorithm where the LT
parameters $(k_{j},d_{j}),j=1,\ldots,J$ are updated once throughout the
algorithm. From Hoerl et al. (1975), Liu (2003), we propose to estimate
${\widehat{k}}_{LT,j}=p{\widehat{\sigma}}_{R,j}/{\widehat{\bf\beta}}_{R,j}^{\top}{\widehat{\bf\beta}}_{R,j}$
and ${\widehat{d}}_{j}$ from (28) where ${\widehat{\bf\beta}}_{R,j}$ and
${\widehat{\sigma}}_{R,j}$ come from (22) and (23), respectively.

Liu-type SEM Algorithm: Similar to the SEM algorithms, the S-step partition
the observations stochastically from
$\text{Multi}(1,\tau_{i1}({\bf\Psi}^{(r)}),\ldots,\tau_{iJ}({\bf\Psi}^{(r)}))$
for $i=1,\ldots,n$. Once the partition established, the rest of the Liu-type
SEM algorithms are implemented similar to the Liu-type CEM algorithms.

##  3 Simulation Studies

In this section, we examine the performance of the ML, Ridge and Liu-type (LT)
methods when the underlying population is a mixture of linear regression
models with multicollinearity problem. We present two simulations enabling us
to study the effect of the sample size, multicollinearity levels and the
number of components of mixture models on the estimation and prediction of the
proposed methods.

Figure 3: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\text{SSE}(\widehat{\sigma}^{2})$ of the estimators when the population is a mixture of two regression models with $n=60$.  Table 1: The median (M) and the length (L) of 95% CIs for the RMSEP of the ML, ridge and LT methods in predicting the mixture of two regression models when $n=60$. | $\rho$ |  | 0.88 |  | 0.90 |  | 0.95 |  | 0.97 |  | 0.99
---|---|---|---|---|---|---|---|---|---|---|---
Method | Algorithm |  | M | L |  | M | L |  | M | L |  | M | L |  | M | L
ML | EM |  | 16.3 | 10.8 |  | 16.7 | 11.1 |  | 17.7 | 12.5 |  | 18.1 | 12.1 |  | 18.3 | 12.3
| CEM |  | 16.5 | 10.7 |  | 16.7 | 11.1 |  | 17.7 | 12.0 |  | 18.0 | 12.0 |  | 18.2 | 12.3
| SEM |  | 16.3 | 11.3 |  | 16.7 | 11.4 |  | 17.6 | 11.9 |  | 17.9 | 11.9 |  | 18.2 | 12.4
Ridge | EM |  | 16.5 | 11.0 |  | 16.7 | 11.2 |  | 17.9 | 12.1 |  | 18.1 | 12.0 |  | 18.3 | 12.0
| CEM |  | 16.4 | 10.8 |  | 16.7 | 11.1 |  | 17.8 | 11.6 |  | 17.9 | 11.9 |  | 18.3 | 12.7
| SEM |  | 16.4 | 11.0 |  | 16.7 | 11.5 |  | 17.7 | 11.6 |  | 18.0 | 11.8 |  | 18.3 | 12.3
LT(HKP) | EM |  | 16.4 | 10.6 |  | 16.6 | 11.1 |  | 17.7 | 11.7 |  | 18.1 | 12.5 |  | 18.3 | 12.0
| CEM |  | 16.4 | 11.3 |  | 16.6 | 11.2 |  | 17.6 | 11.9 |  | 18.1 | 11.6 |  | 18.3 | 12.2
| SEM |  | 16.4 | 10.9 |  | 16.6 | 11.2 |  | 17.8 | 12.0 |  | 18.0 | 12.1 |  | 18.3 | 12.1
LT(ITE) | EM |  | 16.4 | 10.8 |  | 16.8 | 10.9 |  | 17.3 | 11.4 |  | 17.7 | 11.8 |  | 18.1 | 12.0
| CEM |  | 16.5 | 10.9 |  | 16.7 | 10.9 |  | 17.2 | 11.7 |  | 17.7 | 11.8 |  | 18.1 | 11.9
| ESM |  | 16.4 | 11.1 |  | 16.7 | 11.0 |  | 17.2 | 11.7 |  | 17.7 | 12.3 |  | 18.1 | 11.9

In the first study, we simulate data from a population corresponding to a
mixture of two regression models with four covariates $(x_{1},\ldots,x_{4})$.
Following Inan and Erdogan (2013), we use $\rho$ denoting the correlation
between covariates to simulate the multicollinearity in the mixture model. To
do so, we first generate random numbers
$\\{w_{ij};i=1,\ldots,n;j=1,\ldots,5\\}$ from the standard normal
distribution. The covariates are then generated by

$x_{ij}=(1-\rho^{2})^{1/2}w_{ij}+\rho w_{i,5},~{}~{}~{}j=1,\ldots,4,$

where we set $\rho=\\{0.88,0.9,0.95,0.97,0.99\\}$ to simulate the
multicollinearity levels in the mixture of regression models. The responses
are then generated from mixture model (4) whose true parameters are given by
${\bf\Psi}_{0}=(\pi_{0},{\bf\beta}_{01},{\bf\beta}_{02},\sigma^{2}_{01},\sigma^{2}_{02})$
with $\pi_{0}=0.7$, ${\bf\beta}_{01}=(1,3,4,5,6)$ ,
${\bf\beta}_{02}=(-1,-1,-2,-3,-5)$ and $\sigma^{2}_{01}=\sigma^{2}_{02}=1$. We
measured the estimation performance of the ML, Ridge and LT methods by sum of
squared errors (SSE) in estimating $(\pi,{\bf\beta},\sigma^{2})$ and computed
$\text{SSE}(\widehat{\bf\beta})=\left[(\widehat{\bf\beta}-{\bf\beta}_{0})^{\top}(\widehat{\bf\beta}-{\bf\beta}_{0})\right]$,
$\text{SSE}({\widehat{\pi}})=(\widehat{\pi}-\pi_{0})^{2}$ and
$\text{SSE}(\widehat{\sigma}^{2})=\left[(\widehat{\sigma}^{2}-{\sigma_{0}^{2}})^{\top}(\widehat{\sigma}^{2}-{\sigma_{0}^{2}})\right]$,
where $\widehat{\bf\beta}=(\widehat{\bf\beta}_{1},\widehat{\bf\beta}_{2})$,
$\widehat{\sigma}^{2}=({\widehat{\sigma}}_{1}^{2},{\widehat{\sigma}}_{2}^{2})$,
$\widehat{\bf\beta}_{0}=(\widehat{\bf\beta}_{01},\widehat{\bf\beta}_{02})$ and
${\sigma_{0}^{2}}=({\sigma_{01}^{2}},{\sigma_{02}^{2}})$. We also used the
root mean squared errors of prediction (RMSEP) to evaluate the prediction
performance of the methods. To do so, we first compute the RMSEP of the $m$-th
replicate through $K$-fold cross-validation by

$\text{RMSEP}^{(m)}=\left(1/n\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i}^{(m)})^{2}\right)^{1/2},$

where $\widehat{y}_{i}^{(m)}$ is the predicted response of the $i$-th
observation in the $m$-th replicate.

We computed the estimation and prediction measures for the ML, Ridge and LT
methods as follows. We first generated a sample of size $n=\\{60,100\\}$ from
the underlying mixture of regression models as described above. We then used
the EM, CEM and SEM algorithms to estimate the parameters of the mixture
population via ML, ridge and LT methods. We applied the idea of $K=5$ cross-
validation to assess the prediction performance of the methods. To this end,
we divided the sample into $K$ folds of equal sizes. We used $K-1$ folds for
training and the remaining fold for prediction. We repeated the procedure for
all $k=1,\ldots,K$ to compute the $\text{RMSEP}^{(m)}$. Eventually, we
replicated the entire procedures $m=2000$ times and computed the median and
95% confidence interval (CI) for the SSE and RMSEP measures. The lower and
upper bounds of the CI correspond to 2.5 and 97.5 percentiles of 2000
replications, respectively.

Figure 4: The median (M), lower (L) and upper (U) bounds of 95% CIs for
$\text{SSE}(\widehat{\pi})$ of the estimators when the population is a mixture
of three regression models with $n=60$.

We show the results of the simulation study in estimating
$(\pi,{\bf\beta},{\bf\sigma}^{2})$ in Figures 1-7. It is observed that the ML
methods estimate slightly better the mixing proportion than the ridge and LT
methods. This happens because the ridge and LT estimators are biased shrinkage
methods where a slight bias is incorporated into the estimation to encounter
the multicollinearity problem. We observe that the multicollinearity
significantly impacts the ML estimates of the coefficients and results in
extremely unreliable estimates for all EM, CEM and SEM algorithms. Unlike ML
estimates, a significant improvement is seen in the performance of the
shrinkage methods in estimating the coefficients of the component regressions.
The LT methods appear more reliable than their ridge counterparts in the
multicollinearity. From a comparison between LT(ITR) and LT(HKP), we see that
LT(HKP) provides more reliable estimates for ${\bf\sigma}^{2}$. Among the
LT(HKP) estimators, the CEM algorithm almost always outperforms its EM ad SEM
counterparts. Tables 1-3 show the median and the length of 95% CI of the RMSEP
for all the developed methods. The tables clearly show that the prediction
performances of all the methods and EM algorithms are almost identical. This
finding is consistent with Inan and Erdogan (2013), Ghanem et al. (2022) that
multicollinearity seriously affects the estimation of the methods while
prediction levels stay almost the same.

Table 2: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\sqrt{\text{SSE}}$ of the methods in the analysis of bone mineral data with sample size $n=60$. |  |  | CEM |  | SEM |  | EM
---|---|---|---|---|---|---|---
Methods | $\boldsymbol{\Psi}$ |  | M | L | U |  | M | L | U |  | M | L | U
ML | $\beta$ |  | .010 | .002 | .165 |  | .019 | .003 | .213 |  | .018 | .003 | .134
| $\pi$ |  | .333 | .100 | .366 |  | .333 | .183 | .366 |  | .218 | .015 | .365
| $\sigma^{2}$ |  | .003 | .000 | .014 |  | .006 | ..000 | .014 |  | .004 | .000 | .014
Ridge | $\beta$ |  | .009 | .002 | .165 |  | .013 | .002 | .166 |  | .012 | .002 | .118
| $\pi$ |  | .333 | .100 | .366 |  | .333 | .166 | .366 |  | .214 | .019 | .366
| $\sigma^{2}$ |  | .003 | .000 | .014 |  | .006 | .000 | .014 |  | .004 | .000 | .014
LT(HKP) | $\beta$ |  | .009 | .002 | .165 |  | .010 | .003 | .183 |  | .010 | .003 | .067
| $\pi$ |  | .333 | .100 | .366 |  | .333 | .150 | .366 |  | .205 | .013 | .372
| $\sigma^{2}$ |  | .003 | .000 | .014 |  | .006 | .000 | .014 |  | .004 | .000 | .016
LT(ITE) | $\beta$ |  | .009 | .002 | .010 |  | .009 | .006 | .011 |  | .009 | .007 | .010
| $\pi$ |  | .300 | .100 | .366 |  | .350 | .116 | .566 |  | .575 | .032 | .599
| $\sigma^{2}$ |  | .002 | .000 | .014 |  | .005 | .000 | .014 |  | .003 | .000 | .009

In the second simulation study, we evaluate the performance of the estimators
when the underlying mixture model consists of three component regressions with
two covariates. We set $\rho=\\{0.9,0.92,0.95,0.97,0.99\\}$ to simulate the
multicollinearity in the mixture. As described earlier, we generated the
covariates and responses from the mixture model with parameters
$\pi_{0}=(\pi_{01},\pi_{02},\pi_{03})=(0.3,0.4,0.3)$,
${\bf\beta}_{01}=(1,3,4)$, ${\bf\beta}_{02}=(-1,-1,-2)$,
${\bf\beta}_{03}=(-3,1,-4)$ and $\sigma^{2}_{0}=(0.25,1,0.09)$. Similar to the
settings of the first simulation study, we replicated 2000 times all the
estimation and prediction procedures under the EM, CEM and SEM algorithms and
computed the median and 95% CI for the SSEs and RMSEP for size sizes
$n=\\{60,100\\}$. The estimation and prediction results are demonstrated in
Figures 4-12 and Tables 4-5, respectively. Here, we also observe that the ML
method slightly better estimates the mixing proportions; however, the ML
methods result in extremely unreliable estimates for the coefficients of
component regressions. It is easy to see that shrinkage estimators do better
in estimating the component regression parameters. In addition, the LT(HKP)
almost always outperforms other methods and provides a more reliable estimate
of the mixture of regression models. Therefore, the LT(HKP) method based on
the CEM algorithm is recommended to fit the mixture of linear regression
models in multicollinearity.

##  4 Bone Data Analysis

As a bone metabolic disease, osteoporosis occurs when the bone mineral
architecture of the body deteriorates. This deterioration results in skeletal
fragility and a high risk of osteoporotic fractures in different body areas
such as the hip and femur. Osteoporosis and its related diseases significantly
impact the patient’s health and survival. For example, one out of every two
patients with osteoporotic hip fractures can no longer live independently, and
one out of three may die within one year after the medical complication of the
broken bone (Bliuc et al., 2009, Neuburger et al., 2015). The bone mineral
density (BMD) of individuals improves until age 30 and then decreases as
individuals age. The BMD score of an individual is compared with a BMD norm to
determine the bone disorder status. The BMD norm is computed by the mean BMD
scores of healthy adults aged 20-30. The bone status of an individual is
diagnosed as osteoporosis when the BMD score is less than -2.5 SD from the BMD
norm of the population.

The BMD scores are obtained vis an expensive and time-consuming procedure.
Despite this, the researchers have access to various easy-to-measure patients’
characteristics, such as age, weight, BMI, and test results from earlier
surveys (Kim et al., 2012, Felson et al., 1993). Regression models are among
the most common methods to investigate the impact of a set of patients’
characteristics on the BMD responses. The impact of the characteristics may
differ at different BMD levels. Hence, the inference on BMD measurements can
be handled as a problem of the mixture of linear regression models.

The bone mineral data in this section were obtained from the National Health
and Nutritional Examination Survey (NHANES III) administrated by the CDC on
more than 33999 American adults between years 1988 to 1994. One hundred
eighty-two white women aged 50 and older participated in all two bone
examinations during the survey. Owing to the significant impact of
osteoporosis on older women, we treated these 182 women as our underlying
population. We considered the femur BMD from the second bone examination as
the response variable. We also used two easily attainable physical
characteristics as the explanatory variables of the regression. These physical
characteristics include arm and bottom circumferences. The association between
the two explanatory variables is $\rho=0.81$, indicating the multicollinearity
in the regression. Using the ML method based on the information of all
individuals in the population, the BIC criterion suggests that a mixture of
two regression models is the best fit with parameters
$\pi_{0}=(\pi_{01},\pi_{02})=(0.64,0.36)$, ${\bf\beta}_{01}=(0.004,0.006)$,
${\bf\beta}_{02}=(0.011,0.005)$ and $\sigma^{2}_{0}=(0.0081,0.0144)$. These
parameters were then considered as the true parameters of the mixture
population. We replicated 2000 times the ML, ridge and LT methods under EM,
CEM and SEM algorithms in estimating the parameters of the bone mineral
population based on sample sizes $n=\\{60,100\\}$. We then computed the
estimation and prediction measures
$\sqrt{\text{SSE}(\widehat{\bf\beta})},\sqrt{\text{SSE}(\widehat{\pi})},\sqrt{\text{SSE}(\widehat{\sigma^{2}})}$
and MRSEP using 5-fold cross-validation as described in Section 3.

Tables 2-7 report the median and 95% CI for the above measures in estimating
and predicting the bone mineral population. The CI’s lower (L) and upper (U)
bounds correspond to 2.5 and 97.5 percentiles, respectively. Although all ML,
Ridge and LT methods almost perform identically in estimating the mixing
proportion and component variances, $\widehat{\bf\beta}_{ML}$ become
considerably unreliable. Unlike ML methods, the LT and ridge shrinkage methods
could appropriately handle the multicollinearity in estimating the
coefficients of component regressions. Comparing the shrinkage methods, we
observe that the LT estimators appear more reliable than their ridge
counterparts in estimating the parameters of the bone mineral population.

##  5 Summary and Concluding Remarks

In medical and environmental research (e.g., osteoporosis research), linear
regression models are used as standard statistical methods to investigate the
relationship between a set of covariates with the response variable. When the
underlying population is heterogeneous, the impact of the covariates on the
response may change in different subpopulations. A mixture of linear
regression models can be considered a solution to the problem. The maximum
likelihood (ML) method is a common technique to fit a mixture of regression
models; however, the ML estimates become unreliable in multicollinearity. We
investigated shrinkage methods, including ridge and Liu-type (LT) estimators
based on the EM, CEM and SEM algorithms in estimating the parameters of the
mixture of regression models. We showed that LT estimators outperformed their
ridge and ML counterparts through extensive numerical studies. Finally, we
applied the developed methods to analyze the bone mineral data of women aged
50 and older. While the methods were only applied to an osteoprosis research,
the developed methods are generic and can be applied to other medical and
environmental studies.

## Acknowledgment

Armin Hatefi and Hamid Usefi acknowledge the research support of the Natural
Sciences and Engineering Research Council of Canada (NSERC).

## References

  * Arashi et al. (2014) Mohammad Arashi, BM Golam Kibria, Mina Norouzirad, and Saralees Nadarajah.  Improved preliminary test and stein-rule liu estimators for the ill-conditioned elliptical linear regression model.  _Journal of Multivariate Analysis_ , 126:53–74, 2014.
  * Bliuc et al. (2009) Dana Bliuc, Nguyen D Nguyen, Vivienne E Milch, Tuan V Nguyen, and John A Eisman.  Mortality risk associated with low-trauma osteoporotic fracture and subsequent fracture in men and women.  _Jama_ , 301(5):513–521, 2009.
  * Celeux (1985) Gilles Celeux.  The sem algorithm: a probabilistic teacher algorithm derived from the em algorithm for the mixture problem.  _Computational statistics quarterly_ , 2:73–82, 1985.
  * Celeux and Govaert (1992) Gilles Celeux and Gérard Govaert.  A classification em algorithm for clustering and two stochastic versions.  _Computational statistics & Data analysis_, 14(3):315–332, 1992.
  * Cummings et al. (1995) Steven R Cummings, Michael C Nevitt, Warren S Browner, Katie Stone, Kathleen M Fox, Kristine E Ensrud, Jane Cauley, Dennis Black, and Thomas M Vogt.  Risk factors for hip fracture in white women.  _New England journal of medicine_ , 332(12):767–774, 1995.
  * Dempster et al. (1977) Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Duran et al. (2012) Esra Akdeniz Duran, Wolfgang Karl Härdle, and Maria Osipenko.  Difference based ridge and liu type estimators in semiparametric regression models.  _Journal of Multivariate Analysis_ , 105(1):164–175, 2012.
  * Faria and Soromenho (2010) Susana Faria and Gilda Soromenho.  Fitting mixtures of linear regressions.  _Journal of Statistical Computation and Simulation_ , 80(2):201–225, 2010.
  * Felson et al. (1993) David T Felson, Yuqing Zhang, Marian T Hannan, and Jennifer J Anderson.  Effects of weight and body mass index on bone mineral density in men and women: the framingham study.  _Journal of Bone and Mineral Research_ , 8(5):567–573, 1993.
  * Ghanem et al. (2022) Elsayed Ghanem, Armin Hatefi, and Hamid Usefi.  Liu-type shrinkage estimators for mixture of logistic regressions: An osteoporosis study.  _Liu-type Shrinkage Estimators for Mixture of Logistic Regressions: An Osteoporosis Study_ , Submitted:1–21, 2022.
  * Hatefi et al. (2015) Armin Hatefi, Mohammad Jafari Jozani, and Omer Ozturk.  Mixture model analysis of partially rank-ordered set samples: Age groups of fish from length-frequency data.  _Scandinavian Journal of Statistics_ , 42(3):848–871, 2015.
  * Hatefi et al. (2018) Armin Hatefi, Nancy Reid, Mohammad Jafari Jozani, and Omer Ozturk.  Finite mixture modeling, classification and statistical learning with order statistics.  _Statistica Sinica_ , pages 1–50, 2018.
  * Hawkins et al. (2001) Dollena S Hawkins, David M Allen, and Arnold J Stromberg.  Determining the number of components in mixtures of linear models.  _Computational Statistics & Data Analysis_, 38(1):15–48, 2001.
  * Hoerl and Kennard (1970) Arthur E Hoerl and Robert W Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  _Technometrics_ , 12(1):55–67, 1970.
  * Hoerl et al. (1975) Arthur E Hoerl, Robert W Kannard, and Kent F Baldwin.  Ridge regression: some simulations.  _Communications in Statistics-Theory and Methods_ , 4(2):105–123, 1975.
  * Inan and Erdogan (2013) Deniz Inan and Birsen E Erdogan.  Liu-type logistic estimator.  _Communications in Statistics-Simulation and Computation_ , 42(7):1578–1586, 2013.
  * Jones and McLachlan (1992) PN Jones and Geoffrey J McLachlan.  Fitting finite mixture models in a regression context.  _Australian Journal of Statistics_ , 34(2):233–240, 1992.
  * Kim et al. (2012) Sang Jun Kim, Won-Gyu Yang, Eun Cho, and Eun-Cheol Park.  Relationship between weight, body mass index and bone mineral density of lumbar spine in women.  _Journal of bone metabolism_ , 19(2):95–102, 2012\.
  * Lim et al. (2016) Hee-Sook Lim, Soon-Kyung Kim, Hae-Hyeog Lee, Dong Won Byun, Yoon-Hyung Park, and Tae-Hee Kim.  Comparison in adherence to osteoporosis guidelines according to bone health status in korean adult.  _Journal of bone metabolism_ , 23(3):143–148, 2016.
  * Liu (2003) Kejian Liu.  Using liu-type estimator to combat collinearity.  _Communications in Statistics-Theory and Methods_ , 32(5):1009–1020, 2003.
  * McLachlan et al. (2019) Geoffrey J McLachlan, Sharon X Lee, and Suren I Rathnayake.  Finite mixture models.  _Annual review of statistics and its application_ , 6:355–378, 2019.
  * Melton III et al. (1998) L Joseph Melton III, Elizabeth J Atkinson, Michael K O’connor, W Michael O’fallon, and B Lawrence Riggs.  Bone density and fracture risk in men.  _Journal of Bone and Mineral Research_ , 13(12):1915–1923, 1998.
  * Neuburger et al. (2015) Jenny Neuburger, Colin Currie, Robert Wakeman, Carmen Tsang, Fay Plant, Bianca De Stavola, David A Cromwell, and Jan van der Meulen.  The impact of a national clinician-led audit initiative on care and mortality after hip fracture in england: an external evaluation using time trends in non-audit data.  _Medical care_ , 53(8):686, 2015.
  * Pearce and Hatefi (2021) Andrew David Pearce and Armin Hatefi.  Multiple observers ranked set samples for shrinkage estimators.  _arXiv preprint arXiv:2110.07851_ , 2021.
  * Peel and MacLahlan (2000) DAVID Peel and G MacLahlan.  Finite mixture models.  _John & Sons_, 2000.
  * Quandt and Ramsey (1978) Richard E Quandt and James B Ramsey.  Estimating mixtures of normal distributions and switching regressions.  _Journal of the American statistical Association_ , 73(364):730–738, 1978.
  * Wedel et al. (1998) Michel Wedel, Frenkel Ter Hofstede, and Jan-Benedict EM Steenkamp.  Mixture model analysis of complex samples.  _Journal of Classification_ , 15(2):225–244, 1998\.
  * WHO (1994) WHO.  Assessment of fracture risk and its application to screening for postmenopausal osteoporosis: report of a who study group [meeting held in rome from 22 to 25 june 1992].  1994\.
  * Zhang et al. (2006) Zhiqiang Zhang, Wai Keung Li, and Kam Chuen Yuen.  On a mixture garch time-series model.  _Journal of Time Series Analysis_ , 27(4):577–597, 2006.

##  6 Appendix

###  6.1  Proof of Lemma 1

The positive eigenvalues of ${\bf X}^{\top}{\bf W}_{j}{\bf X}$ and ${\bf
W}_{j}^{1/2}{\bf X}{\bf X}^{\top}{\bf W}_{j}^{1/2}$ must be the same. Hence,
the eigenvalue of ${\bf W}_{j}^{1/2}{\bf X}{\bf X}^{\top}{\bf W}_{j}^{1/2}$
are given by $\lambda_{1j},\ldots,\lambda_{pj}$ and the other $(n-p)$ values
must be zero. From singular value decomposition, it s easy to see ${\bf
V}_{1}={\bf W}_{j}^{1/2}{\bf X}{\bf U}_{j}{{\bf\Lambda}}_{j}^{-1/2}$ and
${{\bf\Lambda}}_{j}^{1/2}={\bf V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf X}{\bf
U}_{j}$. From the definition of ${\bf V}_{1}$ and ${{\bf\Lambda}}_{j}^{1/2}$,
we can show

$\displaystyle{\bf W}_{j}^{1/2}{\bf X}={\bf V}_{1}{\bf V}_{1}^{\top}{\bf
W}_{j}^{1/2}{\bf X}{\bf U}_{j}{\bf U}_{j}^{\top}={\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top}.$ (32)

From (32), we can write the canonical form of the regression by

$\displaystyle{\bf W}_{j}^{1/2}{\bf y}={\bf W}_{j}^{1/2}{\bf
X}{{\bf\beta}}_{j}+{\bf W}_{j}^{1/2}\epsilon={\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top}{{\bf\beta}}_{j}+{\bf
W}_{j}^{1/2}\epsilon={\bf V}_{1}{{\bf\Lambda}}_{j}^{1/2}{{\bf\alpha}}_{j}+{\bf
W}_{j}^{1/2}\epsilon.$ (33)

From (33), we can derive the canonical form of the weighted ridge estimator in
each component by

$\displaystyle{\widehat{\bf\alpha}}_{R,j}$ $\displaystyle=\left(({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2})^{\top}({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2})+k_{j}\mathbb{I}\right)^{-1}({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2})^{\top}{\bf W}_{j}^{1/2}{\bf y}$
$\displaystyle=\left({{\bf\Lambda}}_{j}^{1/2}{\bf V}_{1}^{\top}{\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}+k_{j}\mathbb{I}\right)^{-1}{{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}$
$\displaystyle=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}{{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}.$ $\displaystyle{\bf
U}_{j}{\widehat{\bf\alpha}}_{R,j}=$ $\displaystyle{\bf
U}_{j}\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}{\bf
U}_{j}^{\top}{\bf U}_{j}{{\bf\Lambda}}_{j}^{1/2}{\bf V}_{1}^{\top}{\bf
W}_{j}^{1/2}{\bf y}$ $\displaystyle=$ $\displaystyle\left({\bf
U}_{j}{{\bf\Lambda}}_{j}^{1/2}{{\bf\Lambda}}_{j}^{1/2}{\bf
U}_{j}^{\top}+k_{j}\mathbb{I}\right)^{-1}({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top})^{\top}{\bf W}_{j}^{1/2}{\bf
y}$ $\displaystyle=$ $\displaystyle\left(({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top})^{\top}({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf
U}_{j}^{\top})+k_{j}\mathbb{I}\right)^{-1}({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top})^{\top}{\bf W}_{j}^{1/2}{\bf
y}$ $\displaystyle=$ $\displaystyle\left(({\bf W}_{j}^{1/2}{\bf
X})^{\top}({\bf W}_{j}^{1/2}{\bf X})+k_{j}\mathbb{I}\right)^{-1}({\bf
W}_{j}^{1/2}{\bf X})^{\top}{\bf W}_{j}^{1/2}{\bf y}$ $\displaystyle=$
$\displaystyle\left({\bf X}^{\top}{\bf W}_{j}{\bf
X}+k_{j}\mathbb{I}\right)^{-1}{\bf X}^{\top}{\bf W}_{j}{\bf y}$

$\square$

###  6.2  Proof of Lemma 2

The lemma can be proved in a similar vein to Lemma 1. $\square$

###  6.3  Proof of Lemma 3

From Lemma 1, (32) and (33), we can write the canonical form of the weighted
LT estimator in each component regression as

$\displaystyle{\widehat{\bf\alpha}}_{LT,j}$ $\displaystyle=\left(({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2})^{\top}({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2})+k_{j}\mathbb{I}\right)^{-1}\left(({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2})^{\top}{\bf W}_{j}^{1/2}{\bf
y}-d_{j}{\widehat{\bf\alpha}}_{j}\right)$
$\displaystyle=\left({{\bf\Lambda}}_{j}^{1/2}{\bf V}_{1}^{\top}{\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}+k_{j}\mathbb{I}\right)^{-1}\left({{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}-d_{j}{\widehat{\bf\alpha}}_{j}\right)$
$\displaystyle=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}\left({{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}-d_{j}{\widehat{\bf\alpha}}_{j}\right).$
$\displaystyle{\bf U}_{j}{\widehat{\bf\alpha}}_{R,j}=$ $\displaystyle{\bf
U}_{j}\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}\left({{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}-d_{j}{\widehat{\bf\alpha}}_{j}\right)$
$\displaystyle=$ $\displaystyle\left({\bf
U}_{j}{{\bf\Lambda}}_{j}^{1/2}{{\bf\Lambda}}_{j}^{1/2}{\bf
U}_{j}^{\top}+k_{j}\mathbb{I}\right)^{-1}\left(({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top})^{\top}{\bf W}_{j}^{1/2}{\bf
y}-d_{j}{\bf U}_{j}{\widehat{\bf\alpha}}_{j}\right)$ $\displaystyle=$
$\displaystyle\left(({\bf V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf
U}_{j}^{\top})^{\top}({\bf V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf
U}_{j}^{\top})+k_{j}\mathbb{I}\right)^{-1}\left(({\bf
V}_{1}{{\bf\Lambda}}_{j}^{1/2}{\bf U}_{j}^{\top})^{\top}{\bf W}_{j}^{1/2}{\bf
y}-d_{j}{\widehat{\bf\beta}}_{j}\right)$ $\displaystyle=$
$\displaystyle\left(({\bf W}_{j}^{1/2}{\bf X})^{\top}({\bf W}_{j}^{1/2}{\bf
X})+k_{j}\mathbb{I}\right)^{-1}\left(({\bf W}_{j}^{1/2}{\bf X})^{\top}{\bf
W}_{j}^{1/2}{\bf y}-d_{j}{\widehat{\bf\beta}}_{j}\right)$ $\displaystyle=$
$\displaystyle\left({\bf X}^{\top}{\bf W}_{j}{\bf
X}+k_{j}\mathbb{I}\right)^{-1}\left({\bf X}^{\top}{\bf W}_{j}{\bf
y}-d_{j}{\widehat{\bf\beta}}_{j}\right).$

$\square$

###  6.4  Proof of Lemma 4

i) Since ${\widehat{\bf\alpha}}_{ML,j}={{\bf\Lambda}}_{j}^{-1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}$, it is easy to show that

$\displaystyle{\widehat{\bf\alpha}}_{LT,j}=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}({{\bf\Lambda}}_{j}-d_{j}\mathbb{I}){\widehat{\bf\alpha}}_{ML,j}.$
(34)

From (34), the bias and covariance of ${\widehat{\bf\alpha}}_{LT,j}$ are
computed by

$\displaystyle\text{Bias}({\widehat{\bf\alpha}}_{LT,j})$
$\displaystyle=\mathbb{E}({\widehat{\bf\alpha}}_{LT,j})-{{\bf\alpha}}_{j}$
$\displaystyle=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}({{\bf\Lambda}}_{j}-d_{j}\mathbb{I}){{\bf\alpha}}_{j}-{{\bf\alpha}}_{j}$
$\displaystyle=-\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}(k_{j}+d_{j}){{\bf\alpha}}_{j}.$
$\displaystyle\text{cov}({\widehat{\bf\alpha}}_{LT,j})$
$\displaystyle=\sigma_{j}^{2}({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}({{\bf\Lambda}}_{j}-d_{j}\mathbb{I}){{\bf\Lambda}}_{j}^{-1}({{\bf\Lambda}}_{j}-d_{j}\mathbb{I})({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}.$

Following Liu (2003), we then can find the
$\text{MSE}({\widehat{\bf\alpha}}_{LT,j})$ using the bias and covariance as
follows

$\displaystyle\text{MSE}({\widehat{\bf\alpha}}_{LT,j})$
$\displaystyle=||\text{Bias}({\widehat{\bf\alpha}}_{LT,j})||^{2}+\text{tr}(\text{cov}({\widehat{\bf\alpha}}_{LT,j}))$
$\displaystyle=\sum_{m=1}^{p}(d_{j}+k_{j})^{2}\alpha_{m}^{2}/(\lambda_{m}+k_{j})^{2}+\sigma_{j}^{2}\sum_{m=1}^{p}(d_{j}-\lambda_{m})^{2}/\lambda_{m}(\lambda_{m}+k_{j})^{2}.$

Differentiating $\text{MSE}({\widehat{\bf\alpha}}_{LT,j})$ with resect to
$d_{j}$, it is easy to obtain

$d_{opt,j}\sum_{m=1}^{p}\left((\lambda_{m}\alpha_{m}^{2}+\sigma_{j}^{2})/\lambda_{m}(\lambda_{m}+k_{j})^{2}\right)=\sum_{m=1}^{p}\left((\sigma_{j}^{2}-k_{j}\alpha_{m}^{2})/(\lambda_{m}+k_{j})^{2}\right).$

ii) From Ridge method, one can easily show
${\widehat{\bf\alpha}}_{R,j}=({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}{{\bf\Lambda}}_{j}^{1/2}{\bf
V}_{1}^{\top}{\bf W}_{j}^{1/2}{\bf y}$. Thus,

$\displaystyle{\widehat{\bf\alpha}}_{LT,j}=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}\left({{\bf\Lambda}}_{j}+(k_{j}-d_{j})\mathbb{I}\right){\widehat{\bf\alpha}}_{R,j}.$
(35)

From (35), the bias and covariance of ${\widehat{\bf\alpha}}_{LT,j}$ are
computed by

$\displaystyle\text{Bias}({\widehat{\bf\alpha}}_{LT,j})$
$\displaystyle=\mathbb{E}({\widehat{\bf\alpha}}_{LT,j})-{{\bf\alpha}}_{j}$
$\displaystyle=\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}({{\bf\Lambda}}_{j}+(k_{j}-d_{j})\mathbb{I})\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}{{\bf\alpha}}_{j}-{{\bf\alpha}}_{j}$
$\displaystyle=\left[\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}({{\bf\Lambda}}_{j}+(k_{j}-d_{j})\mathbb{I})\left({{\bf\Lambda}}_{j}+k_{j}\mathbb{I}\right)^{-1}-\mathbb{I}\right]{{\bf\alpha}}_{j}.$
$\displaystyle\text{cov}({\widehat{\bf\alpha}}_{LT,j})$
$\displaystyle=\sigma_{j}^{2}({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}\left({{\bf\Lambda}}_{j}+(k_{j}-d_{j})\mathbb{I}\right)({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}{{\bf\Lambda}}_{j}({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}\left({{\bf\Lambda}}_{j}+(k_{j}-d_{j})\mathbb{I}\right)({{\bf\Lambda}}_{j}+k_{j}\mathbb{I})^{-1}.$

Following Liu (2003), we then can find the
$\text{MSE}({\widehat{\bf\alpha}}_{LT,j})$ using the bias and covariance as
above.

Similar to part (i), differentiating
$\text{MSE}({\widehat{\bf\alpha}}_{LT,j})$ with resect to $d_{j}$, it is easy
to obtain

$d_{opt,j}\sum_{m=1}^{p}\left(\lambda_{m}(\sigma_{j}^{2}+\lambda_{m}\alpha_{m}^{2})/(\lambda_{m}+k_{j})^{4}\right)=\sum_{m=1}^{p}\left(\lambda_{m}(\sigma_{j}^{2}-k_{j}\alpha_{m}^{2})/(\lambda_{m}+k_{j})^{3}\right).$

$\square$

###  6.5  Proof of Lemma 5

The lemma can be proved in a similar vein to Lemma 3. $\square$

Figure 5: The median (M), lower (L) and upper (U) bounds of $95\%$ CIs for $\text{SSE}(\widehat{\bf\beta})$ of the estimators when the population is a mixture of two regression models with $n=100$.  Figure 6: The median (M), lower (L) and upper (U) bounds of $95\%$ CIs for $\text{SSE}(\widehat{\pi})$ of the estimators when the population is a mixture of two regression models with $n=100$.  Figure 7: The median (M), lower (L) and upper (U) bounds of $95\%$ CIs for $\text{SSE}({\widehat{\sigma}}^{2})$ of the estimators when the population is a mixture of two regression models with $n=100$.  Figure 8: The median (M), lower (L) and upper (U) bounds of $95\%$ CIs for $\text{SSE}(\widehat{\pi})$ of the estimators when the population is a mixture of three regression models with $n=100$.  Figure 9: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\text{SSE}({\widehat{\sigma}}^{2})$ of the estimators when the population is a mixture of three regression models with $n=100$.  Figure 10: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\text{SSE}(\widehat{\bf\beta})$ of the estimators when the population is a mixture of three regression models with $n=100$.  Figure 11: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\text{SSE}(\widehat{\pi})$ of the estimators when the population is a mixture of three regression models with $n=100$.  Figure 12: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\text{SSE}({\widehat{\sigma}}^{2})$ of the estimators when the population is a mixture of three regression models with $n=100$.  Table 3: The median (M) and the length (L) of 95% CIs for the RMSEP of the ML, ridge and LT methods in predicting the mixture of two regression models when $n=100$. | $\rho$ |  | 0.88 |  | 0.90 |  | 0.95 |  | 0.97 |  | 0.99
---|---|---|---|---|---|---|---|---|---|---|---
Method | Algorithm |  | M | L |  | M | L |  | M | L |  | M | L |  | M | L
ML | EM |  | 16.6 | 8.2 |  | 16.9 | 8.8 |  | 17.7 | 9.0 |  | 18.1 | 8.9 |  | 18.5 | 9.3
| CEM |  | 16.5 | 8.6 |  | 16.8 | 8.7 |  | 17.6 | 9.4 |  | 18.1 | 9.0 |  | 18.5 | 9.9
| SEM |  | 16.6 | 8.5 |  | 16.9 | 8.9 |  | 17.7 | 9.0 |  | 18.0 | 9.4 |  | 18.5 | 9.8
Ridge | EM |  | 16.6 | 8.8 |  | 16.8 | 8.6 |  | 17.7 | 9.1 |  | 18.1 | 9.2 |  | 18.6 | 9.4
| CEM |  | 16.5 | 8.6 |  | 16.8 | 8.6 |  | 17.6 | 9.4 |  | 18.0 | 9.4 |  | 18.4 | 9.3
| SEM |  | 16.5 | 8.7 |  | 16.7 | 8.5 |  | 17.8 | 9.0 |  | 18.0 | 9.3 |  | 18.5 | 9.4
LT(HKP) | EM |  | 16.5 | 8.6 |  | 16.8 | 8.7 |  | 17.6 | 9.0 |  | 18.1 | 8.9 |  | 18.5 | 9.9
| CEM |  | 16.4 | 8.8 |  | 16.8 | 8.5 |  | 17.6 | 8.9 |  | 18.1 | 9.4 |  | 18.4 | 9.6
| SEM |  | 16.5 | 8.6 |  | 16.9 | 8.8 |  | 17.7 | 8.9 |  | 18.0 | 9.2 |  | 18.4 | 9.6
LT(ITE) | EM |  | 16.6 | 8.5 |  | 16.9 | 8.4 |  | 17.5 | 9.3 |  | 17.9 | 9.3 |  | 18.2 | 9.0
| CEM |  | 16.5 | 8.6 |  | 16.8 | 8.8 |  | 17.3 | 9.0 |  | 17.8 | 9.2 |  | 18.2 | 9.3
| ESM |  | 16.6 | 8.4 |  | 16.8 | 8.5 |  | 17.4 | 8.8 |  | 17.8 | 9.7 |  | 18.2 | 9.5
Table 4: The median (M) and the length (L) of 95% CIs for the RMSEP of the ML, ridge and LT methods in predicting the mixture of three regression models when $n=60$. | $\rho$ |  | 0.90 |  | 0.92 |  | 0.95 |  | 0.97 |  | 0.99
---|---|---|---|---|---|---|---|---|---|---|---
Method | Algorithm |  | M | L |  | M | L |  | M | L |  | M | L |  | M | L
ML | EM |  | 6.0 | 4.0 |  | 6.1 | 3.9 |  | 6.2 | 4.1 |  | 6.4 | 4.1 |  | 6.4 | 4.3
| CEM |  | 6.1 | 3.9 |  | 6.2 | 4.0 |  | 6.4 | 4.3 |  | 6.5 | 4.3 |  | 6.6 | 4.4
| SEM |  | 6.1 | 3.9 |  | 6.2 | 4.1 |  | 6.3 | 4.2 |  | 6.4 | 4.2 |  | 6.5 | 4.3
Ridge | EM |  | 6.0 | 4.0 |  | 6.1 | 3.9 |  | 6.2 | 4.1 |  | 6.3 | 4.2 |  | 6.4 | 4.5
| CEM |  | 6.1 | 4.0 |  | 6.2 | 4.1 |  | 6.3 | 4.2 |  | 6.5 | 4.2 |  | 6.6 | 4.3
| SEM |  | 6.0 | 4.0 |  | 6.2 | 4.2 |  | 6.4 | 4.1 |  | 6.5 | 4.2 |  | 6.5 | 4.3
LT(HKP) | EM |  | 5.9 | 3.9 |  | 6.0 | 3.9 |  | 6.2 | 4.2 |  | 6.3 | 4.2 |  | 6.4 | 4.1
| CEM |  | 6.0 | 4.1 |  | 6.1 | 4.0 |  | 6.3 | 4.3 |  | 6.5 | 4.3 |  | 6.5 | 4.4
| SEM |  | 6.0 | 4.0 |  | 6.1 | 4.0 |  | 6.3 | 4.3 |  | 6.5 | 4.2 |  | 6.6 | 4.4
LT(ITE) | EM |  | 6.5 | 6.9 |  | 6.1 | 3.9 |  | 6.1 | 4.1 |  | 6.3 | 4.0 |  | 6.3 | 4.1
| CEM |  | 6.3 | 4.4 |  | 6.3 | 4.0 |  | 6.4 | 4.2 |  | 6.5 | 4.2 |  | 6.6 | 4.3
| SEM |  | 6.2 | 4.5 |  | 6.2 | 4.0 |  | 6.3 | 4.1 |  | 6.4 | 4.2 |  | 6.5 | 4.3
Table 5: The median (M) and the length (L) of 95% CIs for the RMSEP of the ML, ridge and LT methods in predicting the mixture of three regression models when $n=100$. | $\rho$ |  | 0.90 |  | 0.92 |  | 0.95 |  | 0.97 |  | 0.99
---|---|---|---|---|---|---|---|---|---|---|---
Method | Algorithm |  | M | L |  | M | L |  | M | L |  | M | L |  | M | L
ML | EM |  | 5.9 | 3.0 |  | 6.1 | 3.2 |  | 6.3 | 3.1 |  | 6.4 | 3.4 |  | 6.5 | 3.4
| CEM |  | 6.1 | 3.0 |  | 6.3 | 3.1 |  | 6.4 | 3.3 |  | 6.5 | 3.3 |  | 6.6 | 3.4
| SEM |  | 6.1 | 3.0 |  | 6.1 | 3.2 |  | 6.4 | 3.3 |  | 6.5 | 3.1 |  | 6.5 | 3.4
Ridge | EM |  | 6.0 | 3.1 |  | 6.1 | 3.1 |  | 6.3 | 3.3 |  | 6.4 | 3.2 |  | 6.5 | 3.3
| CEM |  | 6.1 | 3.0 |  | 6.2 | 3.0 |  | 6.4 | 3.2 |  | 6.5 | 3.4 |  | 6.6 | 3.4
| SEM |  | 6.0 | 3.0 |  | 6.2 | 3.1 |  | 6.4 | 3.2 |  | 6.5 | 3.3 |  | 6.6 | 3.5
LT(HKP) | EM |  | 5.9 | 3.1 |  | 6.0 | 3.0 |  | 6.2 | 3.2 |  | 6.4 | 3.2 |  | 6.5 | 3.3
| CEM |  | 6.0 | 3.2 |  | 6.2 | 3.2 |  | 6.3 | 3.2 |  | 6.5 | 3.4 |  | 6.6 | 3.4
| SEM |  | 6.0 | 3.1 |  | 6.1 | 3.1 |  | 6.9 | 3.1 |  | 6.5 | 3.2 |  | 6.5 | 3.4
LT(ITE) | EM |  | 6.4 | 6.2 |  | 6.1 | 3.0 |  | 6.2 | 3.1 |  | 6.3 | 3.0 |  | 6.3 | 3.5
| CEM |  | 6.3 | 3.3 |  | 6.3 | 3.2 |  | 6.4 | 3.3 |  | 6.5 | 3.4 |  | 6.6 | 3.6
| ESM |  | 6.2 | 4.2 |  | 6.2 | 3.0 |  | 6.3 | 3.1 |  | 6.3 | 3.2 |  | 6.4 | 3.4
Table 6: The median (M), lower (L) and upper (U) bounds of 95% CIs for $\sqrt{\text{SSE}}$ of the methods in the analysis of bone mineral data with sample size $n=100$. |  |  | CEM |  | SEM |  | EM
---|---|---|---|---|---|---|---
Methods | $\boldsymbol{\Psi}$ |  | M | L | U |  | M | L | U |  | M | L | U
ML | $\beta$ |  | .010 | .002 | .126 |  | .014 | .003 | .202 |  | .014 | .002 | .112
| $\pi$ |  | .350 | .100 | .380 |  | .360 | .210 | .380 |  | .222 | .016 | .370
| $\sigma^{2}$ |  | .005 | .000 | .014 |  | .006 | ..000 | .014 |  | .003 | .000 | .014
Ridge | $\beta$ |  | .009 | .002 | .123 |  | .011 | .003 | .165 |  | .009 | .002 | .086
| $\pi$ |  | .350 | .100 | .380 |  | .360 | .210 | .380 |  | .220 | .019 | .370
| $\sigma^{2}$ |  | .005 | .000 | .014 |  | .006 | .000 | .014 |  | .003 | .000 | .014
LT(HKP) | $\beta$ |  | .009 | .002 | .123 |  | .010 | .003 | .133 |  | .010 | .002 | .047
| $\pi$ |  | .350 | .100 | .380 |  | .360 | .190 | .380 |  | .207 | .019 | .370
| $\sigma^{2}$ |  | .004 | .000 | .014 |  | .005 | .000 | .014 |  | .003 | .000 | .014
LT(ITE) | $\beta$ |  | .009 | .002 | .010 |  | .009 | .007 | .010 |  | .009 | .007 | .009
| $\pi$ |  | .310 | .100 | .380 |  | .360 | .150 | .580 |  | .584 | .040 | .600
| $\sigma^{2}$ |  | .004 | .000 | .014 |  | .005 | .000 | .014 |  | .002 | .000 | .007
Table 7: The median (M), the lower (L) and the upper (U) bounds of 95% CIs for the RMSEPof the ML, ridge and LT methods in predicting the Bone mineral population when $n=\\{60,100\\}$. | $n$ |  | 60 |  | 100
---|---|---|---|---|---
Method | Algorithm |  | M | L | U |  | M | L | U
ML | EM |  | .139 | .105 | .195 |  | .135 | .110 | .173
| CEM |  | .149 | .114 | .220 |  | .141 | .113 | .190
| SEM |  | .140 | .104 | .234 |  | .137 | .109 | .205
Ridge | EM |  | .137 | .104 | .186 |  | .133 | .109 | .172
| CEM |  | .148 | .115 | .207 |  | .139 | .113 | .189
| SEM |  | .140 | .103 | .223 |  | .136 | .109 | .194
LT(HKP) | EM |  | .135 | .104 | .182 |  | .132 | .109 | .168
| CEM |  | .148 | .115 | .208 |  | .139 | .113 | .188
| SEM |  | .139 | .104 | .221 |  | .137 | .109 | .194
LT(ITE) | EM |  | .125 | .101 | .155 |  | .124 | .104 | .145
| CEM |  | .153 | .119 | .193 |  | .146 | .117 | .181
| SEM |  | .143 | .110 | .186 |  | .140 | .113 | .171


First we investigate the properties of the GFs near $R^{-1}\\{0\\}$,
especially when it is near $Q^{*}$. Recall that
$Q^{*}=\bigcup_{r=m_{0}}^{m}\left(\bigcup_{P,\pi}Q_{P,\pi}^{r}\right)$, whence
by an appropriate coordinate transformation the GF near some $Q_{P,\pi}^{r}$
becomes system (20). Since $Q_{P,\pi}^{r}$ has the same geometry as
$Q_{P}^{r}$, we will only prove the cases for $Q_{P}^{r}$’s; the results then
apply to $Q_{P,\pi}^{r}$ for all permutation $\pi$’s on $\\{1,...,m\\}$.

######  Proposition 6.2 (Local properties of GF).

For $\theta^{*}\in R^{-1}\\{0\\}$, let $\gamma$ denote any GF converging to
$\theta^{*}$. Let $\pi_{s},\pi_{c},\pi_{p}$ be defined as above. Given
$m_{0}\leq r\leq m$ and $n$ Type-II separating samples of $R$. The following
results hold.

  * (a)

Suppose that $n\geq m+(m+m_{0}-r)d$ and $\frac{m_{0}+m}{2}\leq r\leq m$. If
$P$ has deficient number $l=2r-m-m_{0}$, then for almost all $\theta^{*}\in
Q_{P}^{r}$, $\gamma$ converges to $\theta^{*}$ at linear rate and satisfies
$|\pi_{p}(\gamma(t)-\theta^{*})|=O(|\pi_{s}(\gamma(t)-\theta^{*})|)$ as
$t\to+\infty$. When $n\leq(d+1)m_{0}$, the same result holds for any
$\theta^{*}\in R^{-1}\\{0\\}$ sufficiently close to $Q^{*}$.

  * (b)

Suppose that $n\geq m+(m+m_{0}-r)d$. If $m_{0}\leq r<\frac{m+m_{0}}{2}$ or the
deficient number $l$ of $P$ satisfies $l<2r-m-m_{0}$, then for almost all
$\theta^{*}\in Q_{P}^{r}$, if

$l(\gamma)(t):=\int_{t}^{\infty}\sqrt{|\pi_{c}(\gamma(\zeta)-\theta^{*})|^{2}+|\pi_{s}(\gamma(\zeta)-\theta^{*})|^{2}}d\zeta=\Omega(e^{-\mu
t})$

for some $\mu$ greater than the largest negative eigenvalue of
$\operatorname{Hess}R(\theta^{*})$, then
$|\pi_{s}(\gamma(t)-\theta^{*})|=O(|\pi_{c}(\gamma(t)-\theta^{*})|^{2})$ as
$t\to+\infty$.

  * (c)

In particular, (a) holds for all $\theta^{*}\in Q_{P}^{m}$ and (b) holds for
all $\theta^{*}\in Q_{P}^{m_{0}}$.

###### Proof.

  * (a)

First assume that $n\geq m+(m+m_{0}-r)d$ and $\frac{m_{0}+m}{2}\leq r\leq m$.
For almost all $\theta^{*}\in Q_{P}^{r}$ we have

$s=\operatorname{rank}\left(\operatorname{Hess}R(\theta^{*})\right)=\operatorname{codim}\,Q_{P}^{r}.$

Fix any such $\theta^{*}$. Find some $U\ni\theta^{*}$ open, such that for any
$\tilde{\theta}^{*}\in U\cap Q_{P}^{r}$ we have
$\operatorname{rank}\operatorname{Hess}R(\tilde{\theta}^{*})=\operatorname{rank}\left(\operatorname{Hess}R(\theta^{*})\right)$.
Let $\tau:U\to\mathbb{R}^{s}\times\mathbb{R}^{(d+1)m-s}$ be a (smooth)
coordinate transformation such that $\tau(U\cap
Q_{P}^{r})\subseteq\\{0\\}\times\mathbb{R}^{p}$, where $p=(d+1)m-s$, and let
$(0,z^{*}):=\tau(\theta^{*})$, where $0\in\mathbb{R}^{s}$ and
$z^{*}\in\mathbb{R}^{p}$. In this way we obtain the system (20) with $c=0$ and
$U_{p}=\tau(U)$. By Remark 6.1, $\tau(\gamma(t))=(y(t),z(t))$ converges to
$(0,z^{*})$ at linear rate and by Lemma 6.1 (a), $|z(t)-z^{*}|=O(|y(t)|^{2})$.
By transforming back to the original coordinate system, we see that
$\gamma(t)$ converges at linear rate and
$|\pi_{p}(\gamma(t)-\theta^{*})|=O(|\pi_{s}(\gamma(t)-\theta^{*})|)$.

Now assume that $n\leq(d+1)m_{0}$. By Lemma 5.4 (a), there is some open
$U\subseteq\mathbb{R}^{(d+1)m}$ such that $Q^{*}\subseteq U\cap R^{-1}\\{0\\}$
and $R$ is Morse–Bott at $U\cap R^{-1}\\{0\\}$. In particular this means by
applying a coordinate transformation we obtain the system (20) with $c=0$.
Arguing in the same way as above, we can see that $\gamma(t)\to\theta^{*}$ at
linear rate and
$|\pi_{p}(\gamma(t)-\theta^{*})|=O(|\pi_{s}(\gamma(t)-\theta^{*})|)$, as
$t\to+\infty$.

  * (b)

Let $s$ be the maximum of
$\operatorname{rank}\left(\operatorname{Hess}R(\theta^{*})\right)$ for
$\theta^{*}\in Q_{P}^{r}$. By Lemma 5.4, for almost all $\theta^{*}\in
Q_{P}^{r}$ we have
$\operatorname{rank}\left(\operatorname{Hess}R(\theta^{*})\right)=s$. Fix any
such $\theta^{*}$. Find some $U\ni\theta^{*}$ open, such that for any
$\tilde{\theta}^{*}\in U\cap Q_{P}^{r}$,
$\operatorname{rank}\left(\operatorname{Hess}R(\tilde{\theta}^{*})\right)=s$.
Since $r<\frac{m_{0}+m}{2}$, by Lemma 5.4 $s<\operatorname{codim}\,Q_{P}^{r}$,
so $c:=\operatorname{codim}\,Q_{P}^{r}-s>0$. Also define $p:=\dim Q_{P}^{r}$.

At each $\tilde{\theta}^{*}\in Q_{P}^{r}$ there is a coordinate transformation
$\tau_{\tilde{\theta}^{*}}:U\cap(\tilde{\theta}^{*}+(Q_{P}^{r})^{\bot})\to\mathbb{R}^{c}\times\mathbb{R}^{s}$,
parametrized smoothly by $\tilde{\theta}^{*}$, such that

    * i)

$\tau(U\cap(\tilde{\theta}^{*}+V_{s}))=\tau(\tilde{\theta}^{*})+\\{0\\}\times\mathbb{R}^{s}$,
where $V_{s}\subseteq\mathbb{R}^{(d+1)m}$ is the largest subspace such that
for any $v\in V_{s}\setminus\\{0\\}$, $\langle
v,\operatorname{Hess}R(\tilde{\theta}^{*})v\rangle\neq 0$.

    * ii)

$\tau(U\cap(\tilde{\theta}^{*}+V_{c}))=\tau(\tilde{\theta}^{*})+\mathbb{R}^{c}\times\mathbb{R}^{s}$,
where $V_{c}=(V_{s}+Q_{P}^{r})^{\bot}$.

Since the eigenvectors of $\operatorname{Hess}R(\theta)$ depends smoothly on
$\theta$, these $\tau_{\tilde{\theta}^{*}}$’s can be combined into an
embedding $\tau:U\to\mathbb{R}^{c}\times\mathbb{R}^{s}\times\mathbb{R}^{p}$.
In this way we obtain system (20) with $U_{P}=\tau(U)$. By Lemma 6.1 (b) and
Remark 6.2, we thus have $|y(t)|=O(|x(t)|^{2})$. Now apply
$\tau^{-1}:\operatorname{ran}\tau\to U$ to see that
$|\pi_{s}(\gamma(t)-\theta^{*})|=O(|\pi_{c}(\gamma(t)-\theta^{*})|^{2})$.

  * (c)

When $r\in\\{m_{0},m\\}$ and $n\geq r+(m+m_{0}-r)d$, the rank of
$\operatorname{Hess}R$ is constant on $Q_{P}^{r}$.

∎

Since $R$ is continuous, $\lim_{t\to+\infty}R(\gamma(t))$ exists for any GF
$\gamma$ sufficiently close to $R^{-1}\\{0\\}$. We then investigate the
convergence rates of $R(\gamma)$. In particular, we show that even if $\gamma$
does not converge at linear rate, $R(\gamma)$ decreases to $0$ much more
quickly compared to the convergence of $\gamma$.

######  Corollary 6.1 (Local properties of loss).

With the notations in Proposition 6.2, given $m_{0}\leq r\leq m$ and $n$ Type-
II separating samples of $R$. The following results hold.

  * (a)

Suppose that $n\geq m+(m+m_{0}-r)d$ and $\frac{m_{0}+m}{2}\leq r\leq m$ and
$P$ has deficient number $l=2r-m-m_{0}$. For almost all $\theta^{*}\in
Q_{P}^{r}$, we have $R(\gamma(t))\to R(\theta^{*})=0$ at linear rate as
$t\to+\infty$. When $n\leq(d+1)m_{0}$, the same result holds for any
$\theta^{*}\in R^{-1}\\{0\\}$ sufficiently close to $Q^{*}$.

  * (b)

Suppose $m_{0}\leq r<\frac{m+m_{0}}{2}$ and the deficient number $l$ of $P$
satisfies $l<2r-m-m_{0}$. For almost all $\theta^{*}\in Q_{P}^{r}$, there is
some $\beta>0$ such that $R(\gamma(t))=O(|\gamma(t)-\theta^{*}|^{4}+e^{-\beta
t})$ for all sufficiently large $t$.

###### Proof.

  * (a)

Since $R$ is analytic, it is in particular locally Lipschitz, so for any
bounded open $U\subseteq\mathbb{R}^{(d+1)m}$ containing $\theta^{*}$, there is
a constant $c>0$ with

$R(\theta)=|R(\theta)-R(\theta^{*})|\leq c|\theta-\theta^{*}|.$

By Proposition 6.2, for almost all $\theta^{*}\in Q_{P}^{r}$ we have
$|\gamma(t)-\theta^{*}|\leq c^{\prime}e^{-\beta t}$ for some
$c^{\prime},\beta>0$. Therefore,

$R(\gamma(t))\leq c|\gamma(t)-\theta^{*}|\leq cc^{\prime}e^{-\beta t}.$

In other words, $R(\gamma(t))\to 0$ at linear rate as $t\to+\infty$.

Now assume that $n\leq(d+1)m_{0}$. Then there is some open
$U\subseteq\mathbb{R}^{(d+1)m}$ such that $Q^{*}\subseteq U\cap R^{-1}\\{0\\}$
and $R$ is Morse–Bott at $U\cap R^{-1}\\{0\\}$. Thus, for any $\theta^{*}\in
U\cap R^{-1}\\{0\\}$ we can apply Proposition 6.2 and argue in the same way as
above to deduce that $R(\gamma(t))\to 0$ at linear rate as $t\to+\infty$.

  * (b)

By Proposition 6.2 (b), we actually have for almost all $\theta^{*}\in
Q_{P}^{r}$,
$|\pi_{s}(\gamma(t)-\theta^{*})|=O(|\pi_{c}(\gamma(t)-\theta^{*})|^{2}+e^{-\mu
t})$ as $t\to+\infty$. for some $\mu>0$. Using the Talyor expansion of $R$
near $\theta^{*}$, we can write

$\displaystyle R(\theta)=R(\theta^{*})$ $\displaystyle+\langle\nabla
R(\theta^{*}),\theta-\theta^{*}\rangle$
$\displaystyle+\frac{1}{2}\langle\operatorname{Hess}R(\theta^{*})(\theta-\theta^{*}),\theta-\theta^{*}\rangle+O(|\theta-\theta^{*}|^{3}).$

Since $R$ is analytic and $R\geq 0$, we can further write and simplify it as

$\displaystyle R(\theta)$
$\displaystyle=\frac{1}{2}\langle\operatorname{Hess}R(\theta^{*})(\theta-\theta^{*}),\theta-\theta^{*}\rangle+O(|\theta-\theta^{*}|^{4})$
$\displaystyle=\frac{1}{2}\langle\operatorname{Hess}R(\theta^{*})\pi_{s}(\theta-\theta^{*}),\pi_{s}(\theta-\theta^{*})\rangle+O(|\theta-\theta^{*}|^{4})$
$\displaystyle\leq\frac{1}{2}\left\|\operatorname{Hess}R(\theta^{*})\right\||\pi_{s}(\theta-\theta^{*})|^{2}+O(|\theta-\theta^{*}|^{4}).$

Since $\lim_{t\to+\infty}\gamma(t)=\theta^{*}$, it follows that for
sufficiently large $t$ we have

$\displaystyle R(\gamma(t))$
$\displaystyle\leq\frac{1}{2}\left\|\operatorname{Hess}R(\theta^{*})\right\||\pi_{s}(\gamma(t)-\theta^{*})|^{2}+O(|\gamma(t)-\theta^{*}|^{4})$
$\displaystyle=\frac{1}{2}\left\|\operatorname{Hess}R(\theta^{*})\right\|O((|\pi_{c}(\gamma(t)-\theta^{*})|^{2}+e^{-\mu
t})^{2})+O(|\gamma(t)-\theta^{*}|^{4})$
$\displaystyle=\frac{1}{2}\left\|\operatorname{Hess}R(\theta^{*})\right\|O((|\gamma(t)-\theta^{*}|^{2}+e^{-\mu
t})^{2})+O(|\gamma(t)-\theta^{*}|^{4})$
$\displaystyle=O(|\gamma(t)-\theta^{*}|^{4}+e^{-\beta t}),$

where $\beta>0$ is some constant depending on $\mu$. This proves the desired
result.

∎

##  7 Conclusion and Discussion

###  7.1 A Detailed Summary of Results

This part serves as a more detailed version of Section 3.2 (Main results). The
geometry of the loss landscape near $Q^{*}$ is summarized in the following
table.

Sample size | Separation of $Q^{*}$ | Geometry of $Q^{*}$
---|---|---
$\leq(d+1)m_{0}$ | Nowhere separated |  A disjoint union of $Q_{P,\pi}^{r}$’s and each $\overline{Q_{P,\pi}^{r}}$ is an affine subspace in $\mathbb{R}^{(d+1)m}$3
$\geq m+m_{0}d$ |  $Q_{P,\pi}^{m}$ separated1
$\vdots$ | $\vdots$
$\geq r+(m+m_{0}-r)d$ |  $Q_{P,\pi}^{r}$ separated1
$\vdots$ | $\vdots$
$\geq m_{0}+md$ |  $Q_{P,\pi}^{m_{0}}$ separated1
$\geq(d+1)m$ |  $Q^{*}$ is separated2
| 1\. Lemma 5.4
---
2\. Lemma 5.5
3\. Proposition 5.1
Table 1: How sample size determines $Q^{*}$ as a subset of $R^{-1}\\{0\\}$,
under Assumptions 3.1 and 3.2

In the table, the first column lists the important turning points of sample
size. The second and third columns list the separation of the
$Q_{P,\pi}^{r}$’s, and geometry of $Q^{*}$ as a whole, respectively. For
example, when $n\leq(d+1)m_{0}$, no $Q_{P,\pi}^{r}$ is separated. The smallest
$n$ that makes some $Q_{P,\pi}^{r}$ become separated is $m+m_{0}d$; at this
point $Q_{P,\pi}^{m}$ is separated. In fact, the larger the $r$ is, the fewer
samples are needed to make it separated. Finally, when $n\geq(d+1)m$,
$Q^{*}=R^{-1}\\{0\\}$ by Lemma 3.6, so every $Q_{P,\pi}^{r}$ is separated.

Then we summarize the functional properties of $R$ and gradient flows (GF)
near $Q^{*}$ in the following table.

Sample size | $R$ (near $Q_{P,\pi}^{r}$) | Behavior of $\gamma$
---|---|---
Strong Conv. | Weak Conv, | Gen. Stab.
$\leq(d+1)m_{0}$ |  Morse–Bott1 | Yes, linear conv.2 | Yes, linear conv.2 |  No1
|  $\geq r+(m+m_{0}-r)d$,
---
$r\geq\frac{m+m_{0}}{2}$ and $l$ largest
Morse–Bott a.e.1 | Yes4
|  $\geq r+(m+m_{0}-r)d$,
---
$r\leq\frac{m+m_{0}}{2}$ or $l$ not largest
Nowhere Morse–Bott1 |  | Yes, but may not
---
be linear conv.2
| Yes,
---
$R(\gamma)=O(|\gamma-\theta^{*}|^{4}+e^{-\beta t})$3
| 1\. Lemma 5.4
---
2\. Theorem 6.2 and Proposition 6.2
3\. Corollary 6.1
4\. Remark after Definition 6.1
Table 2: How sample size determines $R$ and the GF near $Q^{*}$, under
Assumptions 3.1 and 3.2

In the table, the first column lists several cases of sample sizes. The second
column describes the degeneracy of $\operatorname{Hess}R$. The rest describe
local behaviors of the GF, $\gamma$, including the convergence and convergence
rate of $\gamma$ and $R(\gamma)$ (“Strong Conv.” and “Weak Conv.” in table,
and generalization stability (“Gen. Stab.” in table). When $n\leq(d+1)m_{0}$,
$R$ is a Morse–Bott function near every $Q_{P,\pi}^{r}$, namely,
$\text{rank}\,\operatorname{Hess}R=\text{codim}R^{-1}\\{0\\}$ near every
$Q_{P,\pi}^{r}$. At this sample size, $\gamma(t)\to R^{-1}\\{0\\}$,
$R(\gamma(t))\to 0$ (as $t\to+\infty$) both at exponential rate, but no
$Q_{P,\pi}^{r}$ has generalization stability. When $r\geq(m+m_{0})/2$ and
$n\geq r+(m+m_{0}-r)d$, $R$ is Morse–Bott near $Q_{P,\pi}^{r}$. Since
$Q_{P,\pi}^{r}$ is separated, $\gamma(t)\to Q_{P,\pi}^{r}$ (almost surely) and
$R(\gamma(t))\to 0$ are both linear convergences, and $Q_{P,\pi}^{r}$ has
generalization stability (at which it is separated). Intuitively, it means
near $Q_{P,\pi}^{r}$, the training of the student network is quick, finds the
optimal solution, and stable under perturbation. When $r<(m+m_{0})/2$ and
$n\geq r+(m+m_{0}-r)d$, $R$ is not Morse–Bott near $Q_{P,\pi}^{r}$, and
$\gamma(t)\to Q_{P,\pi}^{r}$ (almost surely) but (by Lemma 5.2 and Remark 5.3)
the convergence may not be linear. Since $\gamma(t)\to Q_{P,\pi}^{r}\subseteq
R^{-1}\\{0\\}$, $R(\gamma(t))\to 0$; in particular $R(\gamma(t))\to
O(|\gamma(t)-\theta^{*}|^{4})$ at exponential rate, where $\theta^{*}$ is the
limit of $\gamma$. Intuitively, it means even if $\operatorname{Hess}R$ is
very degenerate, the training loss still decreases quickly in the training
dynamics, thus, a short training time would often give a significantly small
loss. Note that in this case $Q_{P,\pi}^{r}$ also has generalization
stability.

###  7.2 Conclusion and Future Work

In this paper, we investigate the structure of loss landscape of two-layer NNs
near global minima, determine the set of parameters which give perfect
generalization, and fully characterize the gradient flows around it. Our work
uncovers some simple aspects of the complicated loss landscape and reveals how
model, target function, samples and initialization affect the training
dynamics differently. We also show that when and when not the model is
generalization stable. Surprisingly, these all admit a hierarchical structure
behind them. The results are the fruit of combination of “static” function
theory, the optimization theory, and generalization theory. Meanwhile, we deal
with the problems in novel ways, by developing and introducing tools such as
the approach to prove linear independence of (certain type of) functions, the
Morse—Bott function theory, and a parametrized version of central manifold
theory. The generality of such methods means they could produce more fruitful
results in future theoretical study of machine learning models, including
neural networks.

While we have detailly characterized the loss landscape and training dynamics
near the global minima of two-layer NNs, a “truly” global characterization of
either object is still challenging for us. Indeed, we need to consider not
just the structure of global minima, but also that of local minima, saddles,
plateaus, etc. Besides, unlike the situation in this paper, even the
convergence of gradient flow may not be guaranteed. To solve such problems, we
suggest future works start by characterizing the loss function and gradient
flows near local minima, as we believe they have they are similar to the
global minima in certain aspects (e.g., their geometry).

Another interesting problem is to characterize the loss landscape and training
dynamics near the global minima of arbitrarily deep NNs. The challenge for
this problem mainly lies in the difficulty in determining and proving linear
independence of “NN-type functions” – or more generally, the linear
independence results for each layer of an NN. Indeed, the interplay of non-
linearity from each layer in the composition by activation function is
extremely complicated. However, this result is important because, similar as
our remark for good activation, the feature-maps of the model do not eliminate
information while propagating it from one layer to the next. Inspired by the
paper, we expect these claims (see also [15]) to be true for arbitrarily deep
NNs with most analytic activations. Future work is needed to mathematically
characterize it.

## References

  * [1] P. Absil  “Convergence of Iterates of Descent Methods for Analytic Cost Functions”  In _Society for Industrial and Applied Mathematics_ 16.2, 2005, pp. 531–547
  * [2] J. Carr  “Applications of Center Manifold Theory”, Applied Mathematical Science  Springer-Verlag New York, 1981
  * [3] P. Cheridito, A. Jentzen and F. Rossmannek  “Landscape Analysis for Shallow Neural Networks: Complete Classification of Critical Points for Affine Target Functions”  In _Journal of Nonlinear Science_ Springer, 2022  DOI: https://doi.org/10.1007/s00332-022-09823-8
  * [4] R. Cohen, K. Iga and P. Norbury  “Topics in Morse Theory: Lecture Notes”  Standford, 2006
  * [5] K. Fukumizu and S. Amari  “Local Minima and Plateaus in Hierarchical Structures of Multilayer Perceptrons”  In _Neural Networks_ 13, 2000, pp. 317–327
  * [6] A. Jentzen and A. Riekert  “A Proof of Convergence for the Gradient Descent Optimization Method with Random Initializations in the Training of Neural Networks with ReLU Activation for Piecewise Linear Target Functions”  In _Journal of Machine Learning Research_ 23, 2022, pp. 1–50
  * [7] S.G. Krantz and H.R. Parks  “A Primer of Real Analytic Functions”, Birkhauser Advanced Texts  Springer Science+Business Media New York, 2002
  * [8] J.M. Lee  “Introduction to Smooth Manifolds”, Graduate Texts in Mathematics  Springer Science+Business Media New York, 2013
  * [9] Bo Liu  “Understanding the loss landscape of one-hidden-layer ReLU networks”  In _Knowledge-Based Systems_ 220, 2021, pp. 106923
  * [10] P.Feehan  “Resolution of Singularities and Geometric Proofs of the Lojasiewicz Inequality”  In _Geometry and Topology_ 23.7, 2019, pp. 3273–3313
  * [11] B. Simsek et al.  “Geometry of the Loss Landscape in Overparametrized Neural Networks: Symmetry and Invariances”  In _PMLR_ 139, 2021
  * [12] L. Venturi and J. A.S.  “Spurious Valleys in One-hidden-layer Neural Network Optimization Landscapes”  In _Journal of Machine Learning Research_ 20, 2019
  * [13] Chiyuan Zhang et al.  “Understanding Deep Learning Requires Rethinking Generalization”  In _International Conference on Learning Representations_ , 2017
  * [14] Y. Zhang, Z. Zhang, T. Luo and Z-Q. Xu  “Emebdding Principle of Loss Landscape of Deep Neural Networks”  In _Conference on Neural Information Processing Systems_ , 2021
  * [15] Y. Zhang et al.  “Emebdding Principle: a Hierarchical Structure of Loss Landscape of Deep Neural Networks”  In _Journal of Machine Learning_ , 2022, pp. 60–113