Parameter-Expanded ECME Algorithms for
Logistic and Penalized Logistic Regression

In a variety of estimation problems, the convergence speed of the EM algorithm can be improved by viewing the complete-data model as embedded within a larger model that has additional parameters and then performing parameter updates with respect to the larger model. An EM algorithm constructed from a larger, parameter-expanded model is typically referred to as a parameter-expanded EM algorithm (PX-EM) algorithm (Liu et al. (1998)).

To describe the steps involved in the PX-EM algorithm, we consider a general setup where one is interested in estimating a parameter vector $\mbox{\boldmath$\beta$}\in\mathcal{B}$ and one has defined a complete-data vector $\mathbf{v}$ with density function $g_{c}(\mathbf{v}|\mbox{\boldmath$\beta$})$. The complete-data vector $\mathbf{v}\in\mathcal{V}$ and observed-data vector $\mathbf{y}\in\mathcal{Y}$ are connected through $\mathbf{y}=h(\mathbf{v})$ where $h:\mathcal{V}\longrightarrow\mathcal{Y}$ is the data-reducing mapping to be used for developing the EM algorithm. As described in Liu et al. (1998), to generate a PX-EM algorithm one must first consider an expanded probability model for the complete-data vector $\mathbf{v}$, and we let $g_{PX}(\mathbf{v}|\mbox{\boldmath$\theta$},\alpha)$ denote the probability density function for this expanded probability model where $(\mbox{\boldmath$\theta$},\alpha)$ are the parameters of the expanded model and $\Theta$ denotes the parameter space for $(\mbox{\boldmath$\theta$},\alpha)$. As outlined in Liu et al. (1998) and Lewandowski et al. (2010), the expanded probability model should satisfy the following conditions in order to generate a well-defined PX-EM algorithm

1. 1.

   The observed-data model is preserved in the sense that there is a many-to-one reduction function $R:\Theta\longrightarrow\mathcal{B}$ such that

   $$g_{O}\big{\{}\mathbf{y}|\mbox{\boldmath$\beta$}=R(\mbox{\boldmath$\theta$},\alpha)\big{\}}=\int_{\mathcal{V}(\mathbf{y})}g_{PX}(\mathbf{v}|\mbox{\boldmath$\theta$},\alpha)d\mathbf{v},\quad\textrm{ for all }(\mbox{\boldmath$\theta$},\alpha)\in\Theta,$$

   where $\mathcal{V}(\mathbf{y})=\{\mathbf{v}\in\mathcal{V}:h(\mathbf{v})=\mathbf{y}\}$.

2. 2.

   There is a “null value” $\alpha_{0}$ of $\alpha$ such that the complete-data model is preserved at the null value in the sense that

   $$g_{PX}(\mathbf{v}|\mbox{\boldmath$\beta$},\alpha_{0})=g_{c}(\mathbf{v}|\mbox{\boldmath$\beta$}),\qquad\textrm{for all }\mbox{\boldmath$\beta$}\in\mathcal{B}.$$

Similar to the EM algorithm, the PX-EM algorithm consists of two steps: a “PX–E” step followed by a “PX–M” step. The PX–E step is similar to the usual “E–step” of an EM algorithm in that a “Q-function” is formed by taking the expectation of the (parameter-expanded) complete-data log likelihood conditional on the observed data and parameter values $(\mbox{\boldmath$\theta$},\alpha)=(\mbox{\boldmath$\beta$}^{(t)},\alpha_{0})$. The “PX–M” step then proceeds by finding the expanded parameter values $(\mbox{\boldmath$\theta$}^{(t+1)},\alpha^{(t+1)})$ which maximize this Q-function and computing the parameter update $\mbox{\boldmath$\beta$}^{(t+1)}$ by applying the reduction function to $(\mbox{\boldmath$\theta$}^{(t+1)},\alpha^{(t+1)})$. The two steps of the PX-EM algorithm can be summarized as:

1. 1.

   PX–E step: Compute the parameter-expanded Q-function

   $$Q_{PX}(\mbox{\boldmath$\theta$},\alpha|\mbox{\boldmath$\beta$}^{(t)},\alpha_{0})=E\Big{\{}\log g_{PX}(\mathbf{v}|\mbox{\boldmath$\theta$},\alpha)\Big{|}\mathbf{y},\mbox{\boldmath$\theta$}=\mbox{\boldmath$\beta$}^{(t)},\alpha=\alpha_{0}\Big{\}}.$$

   (7)

2. 2.

   PX–M step: Find

   $$(\mbox{\boldmath$\theta$}^{(t+1)},\alpha^{(t+1)})=\operatorname*{arg\,max}_{\boldsymbol{\theta,\alpha}}Q_{PX}(\mbox{\boldmath$\theta$},\alpha|\mbox{\boldmath$\beta$}^{(t)},\alpha_{0}),$$

   and update $\mbox{\boldmath$\beta$}^{(t)}$ using

   $$\mbox{\boldmath$\beta$}^{(t+1)}=R(\mbox{\boldmath$\theta$}^{(t+1)},\alpha^{(t+1)}).$$

A chief advantage of the PX-EM algorithm is that it is guaranteed to converge at least as fast as the corresponding EM algorithm. When the parameter-expanded model is constructed so that all parameter updates have a direct, closed-form, the PX-EM algorithm can have substantially faster convergence speed than the EM algorithm while maintaining the stability and computational simplicity of the original EM algorithm.

The Expectation-Conditional Maximization (ECM) algorithm (Meng & Rubin (1993)) and the Expectation-Conditional Maximization Either (ECME) algorithms (Liu & Rubin (1994)) are variations of the EM algorithm where the M-step is replaced by a sequence of conditional maximization steps. In the ECM algorithm, one breaks up the maximization of the Q-function into several conditional maximization steps where, in each conditional maximization step, one only maximizes the Q-function with respect to a subset of the parameters while keeping the remaining parameter values fixed. In the ECME algorithm, one also performs a series of conditional maximization steps. However, in the ECME algorithm, some of the conditional maximization are performed with respect to the observed-data log-likelihood function while the remaining conditional maximization steps are performed with respect to the usual Q-function.

As in the original ECME algorithm, in a parameter-expanded version of ECME (PX-ECME), one performs a sequence of conditional maximization steps for the expanded parameters $(\mbox{\boldmath$\theta$},\alpha)$ with some of the conditional maximization steps being performed with respect to the parameter-expanded Q-function while the other conditional maximization steps are performed with respect to the observed-data log-likelihood function. Specifically, when using the observed-data log-likelihood function for conditional maximization one will maximize $\ell_{o,\mathbf{s}}\{R(\mbox{\boldmath$\theta$},\alpha)|\mathbf{y}\}$ with respect to the subset of parameters being maximized in that conditional maximization step. After all of the expanded parameters have been updated through the sequence of conditional maximization steps, the update for the parameter of interest $\beta$ is found by applying the reduction function $R(\cdot)$ to the updated expanded parameters.

One example of a PX-ECME algorithm is the procedure described in Liu (1997) for estimating the parameters of a multivariate-t distribution. Here, the EM algorithm exploits the representation of the multivariate-t distribution as a scale mixture of multivariate normal distributions, and the parameter-expanded model utilized by Liu (1997) includes an expanded scale parameter that plays a role in both the covariance matrix of the normal distribution and the scale of the latent Gamma distribution. This expanded parameter is updated in each iteration through a conditional maximization step maximizing the observed-data likelihood.

One point worth noting is that a PX-ECME algorithm can often be constructed even if the expanded parameter vector $(\mbox{\boldmath$\theta$},\alpha)$ is not identifiable from the complete data. A PX-ECME algorithm can work if the subsets of parameters used in each conditional maximization step are conditionally identifiable from the observed data. Specifically, the parameters in each of the subsets should be identifiable from the observed data assuming that all the other parameters are fixed. For example, if we construct our PX-ECME algorithm so that $\theta$ is updated first followed by an update of $\alpha$, then $\theta$ should be identifiable for a fixed value of $\alpha$, and $\alpha$ should be identifiable from the observed data for a fixed value of $\theta$.

We consider the following parameter-expanded version of model (3) for the complete-data vector $\mathbf{v}=\{(y_{1},W_{1}),\ldots,(y_{n},W_{n})\}$

where the parameters in this model are $\mbox{\boldmath$\theta$}\in\mathbb{R}^{p}$ and $\alpha\in\mathbb{R}$. Under this expanded model, the observed-data model is preserved when using the reduction function $\mbox{\boldmath$\beta$}=R(\mbox{\boldmath$\theta$},\alpha)=\alpha\mbox{\boldmath$\theta$}$, and the original complete-data model (3) is preserved whenever $\alpha=\alpha_{0}=1$.

We now turn to describing how to use the parameter-expanded model (8) to construct a PX-ECME algorithm for maximizing a penalized weighted log-likelihood function with penalty function $P_{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$})$. A PX-ECME algorithm for maximizing a penalized log-likelihood is very similar to that of maximizing a log-likelihood function the only difference being that a penalty function will be subtracted from either the parameter-expanded Q-function or the observed-data log-likelihood function. Throughout the remainder of this paper, we will assume the penalty function has the form $P_{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$})=\sum_{j=1}^{p}P_{\boldsymbol{\eta},j}(\beta_{j})$, where $\mbox{\boldmath$\eta$}\in\mathbb{R}^{q}$ is a vector of tuning parameters. The observed-data weighted penalized log-likelihood function $p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}$ that we seek to maximize is then

and the complete-data penalized weighted log-likelihood corresponding to model (8) is

where $C$ is a constant that does not depend on $(\mbox{\boldmath$\theta$},\alpha)$. The penalized, parameter-expanded Q-function $Q_{PX}^{\boldsymbol{\eta}}(\mbox{\boldmath$\theta$},\alpha|\mbox{\boldmath$\beta$}^{(t)},\alpha_{0})$ is obtained by just subtracting $P_{\boldsymbol{\eta}}(\alpha\mbox{\boldmath$\theta$})$ from the usual parameter-expanded Q-function (7). For model (8), this is given by

where $\mathbf{W}(\mbox{\boldmath$\beta$}^{(t)})$ is the $n\times n$ weight matrix defined in (5).

One feature of the parameter-expanded model (8) is that the parameter vector $(\mbox{\boldmath$\theta$},\alpha)$ is not identifiable from the complete-data vector $\mathbf{v}$. This may be seen by noting that, for a positive scalar c, the parameters $(\alpha,\mbox{\boldmath$\theta$})$ and $(\alpha/c,c\mbox{\boldmath$\theta$})$ both lead to the same complete-data penalized log-likelihood. Nevertheless, $(\mbox{\boldmath$\theta$},\alpha)$ are conditionally identifiable meaning that, for a fixed $\theta$, $\alpha$ is identifiable, and for a fixed $\alpha$, $\theta$ is identifiable. Hence, we can first update $\mbox{\boldmath$\theta$}^{(t)}$ by fixing $\alpha=\alpha^{(t)}$ and maximizing $Q_{PX}^{\boldsymbol{\eta}}(\mbox{\boldmath$\theta$},\alpha|\mbox{\boldmath$\beta$}^{(t)},\alpha_{0})$ with respect to $\theta$

We then update $\alpha^{(t)}$ by maximizing the observed-data penalized log-likelihood

The PX-ECME parameter update of $\mbox{\boldmath$\beta$}^{(t)}$ is given by $\mbox{\boldmath$\beta$}^{(t+1)}=R(\mbox{\boldmath$\theta$}^{(t+1)},\alpha^{(t+1)})=\alpha^{(t+1)}\mbox{\boldmath$\theta$}^{(t+1)}$.

It is worth noting that one can express $\mbox{\boldmath$\theta$}^{(t+1)}$ in terms of the EM parameter update as $\mbox{\boldmath$\theta$}^{(t+1)}=\mbox{\boldmath$\beta$}^{(t+1),EM}/\alpha^{(t)}$ where the EM update $\mbox{\boldmath$\beta$}^{(t+1),EM}$ is obtained from $\mbox{\boldmath$\beta$}^{(t)}$ by maximizing $Q_{PX}^{\boldsymbol{\eta}}(\mbox{\boldmath$\theta$},1|\mbox{\boldmath$\beta$}^{(t)},1)$ with respect to $\theta$. Hence, we can express the PX-ECME update as

where $\rho^{(t+1)}=\operatorname*{arg\,max}_{\rho\in\mathbb{R}}p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\rho\mbox{\boldmath$\beta$}^{(t+1),EM}|\mathbf{y})$. In other words, the PX-ECME update is found by simply first computing the EM update $\mbox{\boldmath$\beta$}^{(t+1),EM}$ of the regression coefficients and then multiplying $\mbox{\boldmath$\beta$}^{(t+1),EM}$ by the scaling factor $\rho^{(t+1)}$, where $\rho^{(t+1)}$ is the scalar maximizing the weighted observed-data penalized log-likelihood function among regression coefficients of the form $\rho\mbox{\boldmath$\beta$}^{(t+1),EM}$. The complete PX-ECME procedure is summarized in Algorithm 1.

Finding $\rho^{(t+1)}$ only involves solving a one-dimensional root-finding problem, namely, $\partial p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\rho\mbox{\boldmath$\beta$}^{(t+1),EM}|\mathbf{y})/\partial\rho=0$, of which there are a wide range of available numerical methods, and in such cases, robust root-finding methods such as Brent’s method (Brent (2013)) may be directly used to find $\rho^{(t+1)}$.

An attractive feature of the PX-ECME algorithm (Algorithm 1) is that it is, like the EM algorithm, a monotone algorithm. The monotonicity property helps to ensure that the procedure stable and robust in practice. The monotonicity of PX-ECME directly follows from the fact that the EM algorithm is monotone with respect to the penalized log-likelihood (i.e., $p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$}^{(t+1),EM}|\mathbf{y})\geq p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$}^{(t)}|\mathbf{y})$) and that $p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\rho^{(t+1)}\mbox{\boldmath$\beta$}^{(t+1),EM}|\mathbf{y})\geq p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$}^{(t+1),EM}|\mathbf{y})$.

The PX-ECME algorithm accelerates convergence by multiplying the EM parameter update by a single scalar, and there are likely other straightforward modifications of the EM algorithm that can produce substantial speed-ups in convergence without sacrificing the monotonicity property.

An attractive class of methods with straightforward parameter updates are those that simply take a linear combination of the previous two EM parameter updates with the weights in the linear combination determined adaptively. One can ensure that this scheme is monotone by comparing the value of the objective function for a parameter update versus that of the current parameter value. One method for finding the weights in this linear combination is the order-1 Anderson acceleration (Walker & Ni (2011)). Applying the order-1 Anderson acceleration to the EM algorithm for penalized logistic regression leads to the following parameter updating scheme

where $\gamma^{(t)}$ is the scalar $\gamma^{(t)}=\mathbf{v}_{t}^{T}\mathbf{r}_{t}/\mathbf{v}_{t}^{T}\mathbf{v}_{t}$ and where $\mathbf{r}_{t}=\mbox{\boldmath$\beta$}^{(t+1),EM}-\mbox{\boldmath$\beta$}^{(t)}$ and $\mathbf{v}_{t}=\mbox{\boldmath$\beta$}^{(t+1),EM}-\mbox{\boldmath$\beta$}^{(t)}+\mbox{\boldmath$\beta$}^{(t-1)}-\mbox{\boldmath$\beta$}^{(t),EM}$. Because the underlying EM algorithm is monotone, only accepting the linear combination $\mbox{\boldmath$\beta$}^{(t+1)}$ if $p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$}^{(t+1)}|\mathbf{y})\geq p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$}^{(t+1),EM}|\mathbf{y})$ ensures that the above iterative scheme is monotone.

The PX-ECME updating scheme outlined in Algorithm 1 implicitly defines a mapping $G:\mathbb{R}^{p}\longrightarrow\mathbb{R}^{p}$ which performs the PX-ECME update of $\mbox{\boldmath$\beta$}^{(t)}$, i.e., $\mbox{\boldmath$\beta$}^{(t+1)}=G_{PX}(\mbox{\boldmath$\beta$}^{(t)})$. Likewise, the EM algorithm also defines a mapping $\mbox{\boldmath$\beta$}^{(t+1),EM}=G_{EM}(\mbox{\boldmath$\beta$}^{(t),EM})$ which, from (6), has the form $G_{EM}(\mbox{\boldmath$\beta$})=(\mathbf{X}^{T}\mathbf{S}\mathbf{W}(\mbox{\boldmath$\beta$})\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{S}\mathbf{u}$. Because the PX-ECME parameter update is a scalar multiple of the EM update, we can express the mapping $G_{PX}$ in terms of $G_{EM}$ as

where $h:\mathbb{R}^{p}\longrightarrow\mathbb{R}$ is the function defined as $h(\mbox{\boldmath$\beta$})=\operatorname*{arg\,max}_{\rho}[p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\rho\mbox{\boldmath$\beta$}|\mathbf{y})]$.

The rate of convergence of a fixed-point iteration of the form $\mbox{\boldmath$\beta$}^{(t+1)}=G(\mbox{\boldmath$\beta$}^{(t)})$ is typically defined as the maximum eigenvalue of the Jacobian matrix of $G$ when the Jacobian is evaluated at the convergence point $\mbox{\boldmath$\beta$}^{*}$ of the algorithm (see, e.g., McLachlan & Krishnan (2007)). As shown in Durante & Rigon (2018), the Jacobian matrix $\mathbf{J}_{EM}(\mbox{\boldmath$\beta$})$ associated with the mapping $G_{EM}(\mbox{\boldmath$\beta$})=(\mathbf{X}^{T}\mathbf{S}\mathbf{W}(\mbox{\boldmath$\beta$})\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{S}\mathbf{u}$ is given by

where $\mathbf{R}(\mbox{\boldmath$\beta$})$ is the $n\times n$ diagonal matrix whose $i^{th}$ diagonal element is $m_{i}\pi(\mathbf{x}_{i}^{T}\mbox{\boldmath$\beta$})[1-\pi(\mathbf{x}_{i}^{T}\mbox{\boldmath$\beta$})]$ and $\pi(u)=1/\{1+\exp(-u)\}$. As we show in Theorem 1 below, the maximum eigenvalue of the Jacobian matrix of $G_{PX}$ at $\mbox{\boldmath$\beta$}^{*}$ can be no larger than the maximum eigenvalue of $\mathbf{J}_{EM}(\mbox{\boldmath$\beta$}^{*})$, and hence the rate of convergence of PX-ECME is no slower than the EM algorithm.

For many penalty functions, the maximization in (10) does not yield an easily-computable solution and may require the use of a time-consuming iterative procedure to compute $\mbox{\boldmath$\theta$}^{(t+1)}$. To better handle such cases, we consider a generalized PX-ECME algorithm where $\mbox{\boldmath$\theta$}^{(t+1)}$ maximizes a more manageable, surrogate parameter-expanded Q-function. The more manageable surrogate Q-functions $\tilde{Q}_{PX}^{\boldsymbol{\eta}}(\mbox{\boldmath$\theta$},\alpha|\mbox{\boldmath$\beta$}^{(t)},\alpha_{0})$ are assumed to have the form

where $\mathbf{H}^{(t)}$ is a positive semi-definite $p\times p$ matrix. As in the PX-ECME algorithm for penalized logistic regression, to update $\mbox{\boldmath$\beta$}^{(t)}$ one first finds $\mbox{\boldmath$\theta$}^{(t+1)}$ and $\alpha^{(t+1)}$ by performing two conditional maximization steps with respect to $\tilde{Q}_{PX}^{\boldsymbol{\eta}}$ and then sets $\mbox{\boldmath$\beta$}^{(t+1)}=\alpha^{(t+1)}\mbox{\boldmath$\theta$}^{(t+1)}$. We will refer to an algorithm which updates $\mbox{\boldmath$\beta$}^{(t)}$ in this manner as a generalized parameter-expanded ECME (GPX-ECME) algorithm for penalized weighted logistic regression.

The steps involved in a GPX-ECME are outlined in Algorithm 2. Notice in this Algorithm that, as in the case of a PX-ECME algorithm, the process of first computing $\mbox{\boldmath$\theta$}^{(t+1)}$ and $\alpha^{(t+1)}$ and then setting $\mbox{\boldmath$\beta$}^{(t+1)}=\alpha^{(t+1)}\mbox{\boldmath$\theta$}^{(t+1)}$ may be equivalently expressed as $\mbox{\boldmath$\beta$}^{(t+1)}=\rho^{(t+1)}\mbox{\boldmath$\beta$}^{(t+1),GEM}$, where $\mbox{\boldmath$\beta$}^{(t+1),GEM}$ maximizes $\tilde{Q}_{PX}^{\boldsymbol{\eta}}(\mbox{\boldmath$\beta$},1|\mbox{\boldmath$\beta$}^{(t)},\alpha_{0})$ and $\rho^{(t+1)}$ maximizes $p\ell_{o,\mathbf{s}}^{\boldsymbol{\eta}}(\rho\mbox{\boldmath$\beta$}^{(t+1),GEM}|\mathbf{y})$ with respect to $\rho\in\mathbb{R}$. As stated in the following proposition, as long as $\mathbf{H}^{(t)}$ is a positive semi-definite matrix for each $t$, each iterate generated from a GPX-ECME algorithm will increase the value of the penalized log-likelihood.