Computation for Latent Variable Model Estimation: A Unified Stochastic Proximal Framework

Abstract

Latent variable models have been playing a central role in psychometrics and related fields. In many modern applications, the inference based on latent variable models involves one or several of the following features: (1) the presence of many latent variables, (2) the observed and latent variables being continuous, discrete, or a combination of both, (3) constraints on parameters, and (4) penalties on parameters to impose model parsimony. The estimation often involves maximizing an objective function based on a marginal likelihood/pseudo-likelihood, possibly with constraints and/or penalties on parameters. Solving this optimization problem is highly non-trivial, due to the complexities brought by the features mentioned above. Although several efficient algorithms have been proposed, there lacks a unified computational framework that takes all these features into account. In this paper, we fill the gap. Specifically, we provide a unified formulation for the optimization problem and then propose a quasi-Newton stochastic proximal algorithm. Theoretical properties of the proposed algorithms are established. The computational efficiency and robustness are shown by simulation studies under various settings for latent variable model estimation.

Latent variable models have been playing a central role in psychometrics and related fields. Commonly used latent variable models include item response theory models (Embretson & Reise 2000; Reckase 2009), latent class models (Clogg, 1995; Rupp et al., 2010; von Davier & Lee, 2019), structural equation models (Bollen, 1989), error-in-variable models (Carroll et al., 2006), random-effects models (Hsiao, 2014), and models for missing data (Little & Rubin, 1987), where latent variables have different interpretations, such as hypothetical constructs, ‘true’ variables measured with error, unobserved heterogeneity, and missing data. We refer the readers to Rabe-Hesketh & Skrondal (2004) and Bartholomew et al. (2011) for a comprehensive review of latent variable models.

A latent variable model contains unobserved latent variables and unknown parameters. For example, an item response theory model contains individual-specific latent traits as latent variables and item-specific parameters as model parameters. Comparing with models without latent variables, such as linear regression and generalized linear regression, the estimation of latent variable models is typically more involved. This estimation problem can be viewed from three perspectives: (1) fixed latent variables and parameters, (2) random latent variables and fixed parameters, and (3) random latent variables and parameters.

The first perspective, i.e., fixed latent variables and parameters, leads to the joint maximum likelihood (JML) estimator. This estimator can often be efficiently computed, for example, by an alternating minimization algorithm (Birnbaum, 1968; Chen et al., 2019; 2020). Unfortunately, however, the JML estimator is typically statistically inconsistent (Neyman & Scott, 1948; Andersen, 1973; Haberman, 1977; Ghosh, 1995), except under some high-dimensional asymptotic regime that is suitable for large-scale applications (Chen et al., 2019; 2020; Haberman, 1977; 2004). Treating both latent variables and parameters as random variables, the third perspective leads to a full Bayesian estimator, for which many Markov chain Monte Carlo (MCMC) algorithms have been developed (e.g., Béguin & Glas, 2001; Bolt & Lall, 2003; Dunson, 2000; 2003; Edwards, 2010).

The second perspective, i.e., random latent variables and fixed parameters, essentially follows an empirical Bayes (EB) approach (Robbins, 1956; C.-H. Zhang, 2003). This perspective is the most commonly adopted one (Rabe-Hesketh & Skrondal, 2004). Throughout the paper, we refer to estimators derived under this perspective as EB estimators. Both the full-information marginal maximum likelihood (MML) estimator (Bock & Aitkin, 1981) and the limited-information composite maximum likelihood (CML) estimator (Jöreskog & Moustaki, 2001; Vasdekis et al., 2012) can be viewed as special cases. Such estimators involve optimizing an objective function with respect to the fixed parameters, while the objective function is often intractable due to an integral with respect to the latent variables. The most commonly used algorithm for this optimization problem is the expectation-maximization (EM) algorithm (Dempster et al., 1977; Bock & Aitkin, 1981). This algorithm typically requires to iteratively evaluate numerical integrals with respective to the latent variables, which is often computationally unaffordable when the dimension of the latent space is high.

A high-dimensional latent space is not the only challenge to the computation of EB estimators. Penalties and constraints on parameters may also involve in the optimization, further complicating the computation. In fact, penalized estimators have become increasingly more popular in latent variable analysis for learning sparse structure, with applications to restricted latent class analysis, exploratory item factor analysis, variable selection in structural equation models, differential item functioning analysis, among others (Chen et al., 2015; Sun et al., 2016; Chen et al., 2018; Lindstrøm & Dahl, 2020; Tutz & Schauberger, 2015; Jacobucci et al., 2016; Magis et al., 2015). The penalty function is often non-smooth (e.g., Lasso penalty, Tibshirani, 1996), for which many standard optimization tools (e.g., gradient descent methods) are not applicable. In addition, complex inequality constraints are also commonly encountered in latent variable estimation, for example, in structural equation models (Van De Schoot et al., 2010) and restricted latent class models (e.g., de la Torre, 2011, Xu, 2017). Such complex constraints further complicate the optimization.

In this paper, we propose a quasi-Newton stochastic proximal algorithm that simultaneously tackles the computational challenges mentioned above. This algorithm can be viewed as an extension of the stochastic approximation (SA) method (Robbins & Monro, 1951). Comparing with SA, the proposed method converges faster and is more robust, thanks to the use of Polyak–Ruppert averaging (Polyak & Juditsky, 1992; Ruppert, 1988). The proposed method can also be viewed as a stochastic version of a proximal gradient descent algorithm (Chapter 4, Parikh & Boyd, 2014), in which constraints and penalties are handled by a proximal update. As will be illustrated by examples later, the proximal update is easy to evaluate for many commonly used penalties and constraints, making the proposed algorithm computationally efficient. Theoretical properties of the proposed method are established, showing that the proposed one is almost optimal in its convergence speed.

The proposed method is closely related to the stochastic-EM algorithm (Celeux, 1985; Ip, 2002; Nielsen, 2000; S. Zhang et al., 2020b) and the MCMC stochastic approximation algorithms (Cai, 2010a; b; Gu & Kong, 1998), two popular methods for latent variable model estimation. Although these methods perform well in many problems, they are not as powerful as the proposed one. Specifically, the MCMC stochastic approximation algorithms cannot handle complex inequality constraints or non-smooth penalties, because they rely on stochastic gradients which do not always exist when there are complex inequality constraints or non-smooth penalties. In addition, as will be discussed later, both the stochastic-EM algorithm and the MCMC stochastic approximation algorithms are computationally less efficient than the proposed method, even for estimation problems without complex constraints or penalties.

The proposed method is also closely related to a perturbed proximal gradient algorithm proposed in Atchadé et al. (2017). The current development improves upon that of Atchadé et al. (2017) from two aspects. First, the proposed method is a Quasi-Newton method, in which the second-order information (i.e., second derivatives) of the objective function is used in the update. Although this step may only change the asymptotic convergence speed by a constant factor (when the number of iterations grows to infinity), our simulation study suggests that the new method converges much faster than that of Atchadé et al. (2017) empirically. Second, the theoretical analysis of Atchadé et al. (2017) only considers a convex optimization setting, while we consider a non-convex setting which is typically the case for latent variable model estimation. Note that the analysis is much more involved when the objective function is non-convex. Therefore, our proof of sequence convergence is different from that of Atchadé et al. (2017). Specifically, the convergence theory is established by analyzing the convergence of a set-valued generalization of an ordinary differential equation (ODE).

The rest of the paper is organized as follows. In Sect. 1, we formulate latent variable model estimation as a general optimization problem which covers many commonly used estimators as special cases. In Sect. 2, a quasi-Newton stochastic proximal algorithm is proposed. Theoretical properties of the proposed algorithm are established in Sect. 3, suggesting that the proposed algorithm achieves the optimal convergence rate. The performance of the proposed algorithm is demonstrated and compared with other estimators by simulation studies in Sect. 4. We conclude with some discussions in Sect. 5. An R package has been developed that can be found on https://github.com/slzhang-fd/lvmcomp2.

1. Estimation of Latent Variable Models

1.1. Problem Setup

We consider the estimation of a parametric latent variable model. We adopt a general setting, followed by concrete examples in Sects. 1.2 and 1.3. Let $Y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Y}}$$\end{document} be a random object representing observable data and let $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}$$\end{document} be its realization. For example, in item factor analysis (IFA), $Y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Y}}$$\end{document} represents (categorical) responses to all the items from all the respondents. A latent variable model specifies the distribution of $Y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Y}}$$\end{document} by introducing a set of latent variables $\xi \in \Xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi } \in \Xi $$\end{document} , where $\Xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Xi $$\end{document} denotes the state space of the latent vector $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} . For example, in item factor analysis, $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} consists of the latent traits of all the respondents and $\Xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Xi $$\end{document} is a Euclidean space. Let $\beta ={({\beta }_1,\dots ,{\beta }_p)}^{\top }\in B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\beta } = (\beta _1,\ldots , \beta _p)^\top \in {\mathcal {B}}$$\end{document} be a set of parameters in the model, where $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}$$\end{document} denotes the parameter space. The goal is to estimate $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\beta }$$\end{document} given observed data $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}$$\end{document} .

We consider an EB estimator which takes the form

(1) $\begin{array}{c}l\left(\beta \right)=log\left({\int }_{\Xi },f,(y,\xi \mid \beta ),d,\xi \right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} l(\varvec{\beta }) = \log \left( \int _{\Xi } f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta })d \varvec{\xi }\right) , \end{aligned}$$\end{document}

where $f(y,\xi \mid \beta )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta })$$\end{document} is a complete-data likelihood/pseudo-likelihood function that has an analytic form. We assume that the objective function $l\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l(\varvec{\beta })$$\end{document} is finite for any $\beta \in B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}\in {\mathcal {B}}$$\end{document} and is also smooth in $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\beta }$$\end{document} .

The estimator is given by solving the following optimization problem

(2) $\begin{array}{c}\hat{\beta }=\underset{\beta \in B}{\text{arg min}}-l\left(\beta \right)+R\left(\beta \right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\hat{{\varvec{\beta }}}} = \mathop {\hbox {arg min}}\limits _{\varvec{\beta } \in {\mathcal {B}}} -l(\varvec{\beta }) + R(\varvec{\beta }), \end{aligned}$$\end{document}

where $R\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R(\varvec{\beta })$$\end{document} is a penalty function that has an analytic form, such as Lasso, ridge, or elastic net regularization functions. Note that the penalty function often depends on tuning parameters. Throughout this paper, we assume these tuning parameters are fixed and thus do not explicitly indicate them in the objective function (2). In practice, tuning parameters are often unknown and need to be chosen by cross-validation or certain information criterion. We point out that many commonly used estimators take the form of (2), including the MML estimator, the CML estimator, and regularized estimators based on the MML and CML. We also point out that despite its general applicability to latent variable estimation problems, the proposed method is more useful for complex problems that cannot be easily solved by the classical EM algorithm. For certain problems, such as the estimation of linear factor models and simple latent class models, both the E- and M-step of the EM algorithm have closed-form solutions. In that situation, the classical EM algorithm may be computationally more efficient, though the proposed method can still be used.

1.2. High-dimensional Item Factor Analysis

Item factor analysis models are commonly used in social and behavioral sciences for analyzing categorical response data. For exposition, we focus on binary response data and point out that the extension to ordinal response data is straightforward. Consider N individuals responding to J binary-scored items. Let $Y_{\mathrm{ij}}\in \{0,1\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{ij}\in \{0, 1\}$$\end{document} be a random variable denoting person i’s response to item j and let $y_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$y_{ij}$$\end{document} be its realization. Thus, we have $Y={\left(Y_{\mathrm{ij}}\right)}_{N\times J}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Y}}= (Y_{ij})_{N\times J}$$\end{document} and $y={\left(y_{\mathrm{ij}}\right)}_{N\times J}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}= (y_{ij})_{N\times J}$$\end{document} , where $Y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Y}}$$\end{document} and $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}$$\end{document} are the generic notations introduced in Sect. 1.1 for our data. A comprehensive review of IFA models and their estimation can be found in Chen & Zhang (2020a).

It is assumed that the dependence among an individual’s responses is driven by a set of latent factors, denoted by $\xi ={\left({\xi }_{\mathrm{ik}}\right)}_{N\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi } = (\xi _{ik})_{N\times K}$$\end{document} , where ${\xi }_{\mathrm{ik}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\xi _{ik}$$\end{document} represents person i’s kth factor. Recall that $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} is our generic notation for the latent variables in Sect. 1.1 and here the state space $\Xi =R^{N\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Xi = {\mathbb {R}}^{N\times K}$$\end{document} . Throughout this paper, we assume the number of factors K is known.

An IFA model makes the following assumptions:

1. 1. ${\xi }_i={({\xi }_{i1},\dots ,{\xi }_{\mathrm{iK}})}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i = (\xi _{i1},\ldots , \xi _{iK})^\top $$\end{document} , $i=1,\dots ,N$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1,\ldots , N$$\end{document} , are independent and identically distributed (i.i.d.) random vectors, following a multivariate normal distribution $N(0,\Sigma )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N({\mathbf {0}}, \varvec{\Sigma })$$\end{document} . The diagonal terms of $\Sigma ={\left({\sigma }_{k{k}^{\prime }}\right)}_{K\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma } = (\sigma _{kk'})_{K\times K}$$\end{document} are set to one for model identification. As $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} is a positive semi-definite matrix, it is common to reparametrize $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} by Cholesky decomposition,

   $\begin{array}{c}\Sigma =B{B}^{\top },\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\Sigma } = {\mathbf {B}}{\mathbf {B}}^\top , \end{aligned}$$\end{document}
   where $B={\left(b_{k{k}^{\prime }}\right)}_{K\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}} = (b_{kk'})_{K\times K}$$\end{document} is a lower triangular matrix. Let $b_k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {b}}_k$$\end{document} be the kth row of $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} . Then $\Vert b_k\Vert =1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert {\mathbf {b}}_k \Vert = 1$$\end{document} , $k=1,\dots ,K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k=1,\ldots , K$$\end{document} , since the diagonal terms of $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} are constrained to value 1.

2. 2. $Y_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{ij}$$\end{document} given ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} follows a Bernoulli distribution satisfying

   (3) $\begin{array}{c}P(Y_{\mathrm{ij}}=1\mid {\xi }_i)=\frac{exp(d_j+a_j^{\top }{\xi }_i)}{1+exp(d_j+a_j^{\top }{\xi }_i)},\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i) = \frac{\exp (d_j + {\mathbf {a}}_j^\top \varvec{\xi }_i)}{1+\exp (d_j + {\mathbf {a}}_j^\top \varvec{\xi }_i)}, \end{aligned}$$\end{document}
   where $d_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j$$\end{document} and $a_j={(a_{j1},\dots ,a_{\mathrm{jK}})}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {a}}_j = (a_{j1},\ldots , a_{jK})^\top $$\end{document} are item-specific parameters. The parameters $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} are often known as the loading parameters.

3. 3. $Y_{i1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{i1}$$\end{document} ,..., $Y_{\mathrm{iJ}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{iJ}$$\end{document} are assumed to be conditionally independent given ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} , which is known as the local independence assumption.

Note that we consider the most commonly used logistic model in (3). It is worth pointing out that the proposed algorithm also applies to the normal ogive (i.e., probit) model which assumes that $P(Y_{\mathrm{ij}}=1\mid {\xi }_i)=\Phi (d_j+a_j^{\top }{\xi }_i)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$ {\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i) =\Phi (d_j + {\mathbf {a}}_j^\top \varvec{\xi }_i)$$\end{document} . Under the current setting and using the reparametrization for $\Sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document} , our model parameters are $\beta =\{B,d_j,a_j,j=1,\dots ,J\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}= \{{\mathbf {B}}, d_j, {\mathbf {a}}_j, j =1,\ldots , J\}$$\end{document} . The marginal likelihood function takes the form

(4) $\begin{array}{c}l\left(\beta \right)=\sum _{i=1}^{N}log\left({\int }_{x\in R^K},\prod _{j=1}^J,\frac{exp\left[y_{\mathrm{ij}}(d_j+a_j^{\top }x)\right]}{1+exp(d_j+a_j^{\top }x)},\varphi ,(x\mid B),d,x\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} l({\varvec{\beta }}) = \sum _{i=1}^N \log \left( \int _{{\mathbf {x}}\in {\mathbb {R}}^{K}} \prod _{j=1}^J \frac{\exp [y_{ij}(d_j + {\mathbf {a}}_j^\top {\varvec{x}})]}{1+\exp (d_j + {\mathbf {a}}_j^\top {\varvec{x}})} \phi ({\mathbf {x}}\mid {\mathbf {B}}) d{\mathbf {x}}\right) , \end{aligned}$$\end{document}

where $\varphi (x\mid B)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\phi ({\mathbf {x}}\mid {\mathbf {B}})$$\end{document} is the density function for multivariate normal distribution $N(0,B{B}^{\top })$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N({\mathbf {0}}, {\mathbf {B}}{\mathbf {B}}^\top )$$\end{document} . The K-dimensional integrals involved in (4) cause a high computational burden for a relatively large K (e.g., $K\ge 5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K \ge 5$$\end{document} ).

IFA models are commonly used for both exploratory and confirmatory analyses. In exploratory IFA, an important problem is to learn a sparse loading matrix ${\left(a_{\mathrm{ij}}\right)}_{J\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(a_{ij})_{J\times K}$$\end{document} from data, which facilitates the interpretation of the factors. One approach is by the $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} -regularized estimator (Sun et al., 2016) which takes the form

(5) $\begin{array}{c}\begin{array}{cc}\hat{\beta }=& \underset{\beta \in B}{\text{arg min}}-l\left(\beta \right)+R\left(\beta \right),\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned} {\hat{{\varvec{\beta }}}} =&\mathop {\hbox {arg min}}\limits _{{\varvec{\beta }}\in {\mathcal {B}}} - l({\varvec{\beta }}) + R({\varvec{\beta }}), \end{aligned} \end{aligned}$$\end{document}

where the parameter space

$\begin{array}{c}B=\{\beta \in R^p:b_{k{k}^{\prime }}=0,1\le k<k^{\prime }\le K,\sum _{k^{\prime }=1}^K{b}_{k{k}^{\prime }}^2=1,k=1,\dots ,K\},\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathcal {B}} = \{{\varvec{\beta }}\in {\mathbb {R}}^{p}: b_{kk'} = 0, 1\le k < k' \le K, \sum _{k' = 1}^K b_{kk'}^2 = 1, k = 1,\ldots , K\}, \end{aligned}$$\end{document}

and the penalty term

(6) $\begin{array}{c}R\left(\beta \right)=\lambda \sum _{j=1}^J\sum _{k=1}^K\left|a_{\mathrm{jk}}\right|.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} R({\varvec{\beta }}) = \lambda \sum _{j = 1}^J\sum _{k=1}^K \vert a_{jk}\vert . \end{aligned}$$\end{document}

In $R\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R({\varvec{\beta }})$$\end{document} , $\lambda >0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda > 0$$\end{document} is a tuning parameter assumed to be fixed throughout this paper. This regularized estimator resolves the rotational indeterminacy issue in exploratory IFA, as the $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} penalty term is not rotational invariant. Consequently, under mild regularity conditions, the loading matrix can be consistently estimated only up to a column swapping. Note that only the $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} matrix has constraints, as reflected by the parameter space $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}$$\end{document} . Here $b_{k{k}^{\prime }}=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$b_{kk'} = 0$$\end{document} is due to that $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} is a lower triangle matrix and ${\sum }_{k^{\prime }=1}^K{b}_{k{k}^{\prime }}^2=1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{k' = 1}^K b_{kk'}^2 = 1$$\end{document} is due to that the diagonal terms of $\Sigma =B{B}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma } = {\mathbf {B}}{\mathbf {B}}^\top $$\end{document} are all 1. We remark that it is possible to replace the $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} penalty in $R\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R({\varvec{\beta }})$$\end{document} by other penalty functions for imposing sparsity, such as the elastic net penalty (Zou & Hastie, 2005)

(7) $\begin{array}{c}R\left(\beta \right)={\lambda }_1\sum _{j=1}^J\sum _{k=1}^K{a}_{\mathrm{jk}}^2+{\lambda }_2\sum _{j=1}^J\sum _{k=1}^K\left|a_{\mathrm{jk}}\right|,\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} R({\varvec{\beta }}) = \lambda _1 \sum _{j = 1}^J\sum _{k=1}^K a_{jk}^2 + \lambda _2 \sum _{j = 1}^J\sum _{k=1}^K \vert a_{jk}\vert , \end{aligned}$$\end{document}

where ${\lambda }_1,{\lambda }_2>0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _1, \lambda _2 > 0$$\end{document} are two tuning parameters.

In confirmatory IFA, zero constraints are imposed on loading parameters, based on prior knowledge about the measurement design. More precisely, these zero constraints can be coded by a binary matrix $Q={\left(q_{\mathrm{jk}}\right)}_{J\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}} = (q_{jk})_{J\times K}$$\end{document} . If $q_{\mathrm{jk}}=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$q_{jk} = 0$$\end{document} , then item j does not load on factor k and $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} is set to 0. Otherwise, $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} is freely estimated. These constraints lead to parameter space $B=\{\beta :b_{k{k}^{\prime }}=0,1\le k<k^{\prime }\le K;{\sum }_{k^{\prime }=1}^K{b}_{k{k}^{\prime }}^2=1,a_{\mathrm{jk}}=0\text{for}{q}_{\mathrm{jk}}=0,j=1,\dots ,J,k=1,\dots ,K\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}} = \{{\varvec{\beta }}: b_{kk'} = 0, 1\le k < k' \le K; \sum _{k' = 1}^K b_{kk'}^2 = 1, a_{jk} = 0 ~\text{ for }~ q_{jk} = 0, j = 1,\ldots ,J, k = 1,\ldots , K\}$$\end{document} . The MML estimator for confirmatory IFA is then given by

(8) $\begin{array}{c}\begin{array}{cc}\hat{\beta }=& \underset{\beta \in B}{\text{arg min}}-l\left(\beta \right).\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned} {\hat{{\varvec{\beta }}}} =&\mathop {\hbox {arg min}}\limits _{{\varvec{\beta }}\in {\mathcal {B}}} - l({\varvec{\beta }}). \end{aligned} \end{aligned}$$\end{document}

Besides parameter estimation, another problem of interest in confirmatory IFA is to make statistical inference, for which it is required to compute the asymptotic variance of $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} . The estimation of the asymptotic variance often requires to compute the Hessian matrix of $l\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l({\varvec{\beta }})$$\end{document} at $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} , which also involves intractable K-dimensional integrals. As we will see in Sect. 2.1, this Hessian matrix, as well as quantities taking a similar form, can be easily obtained as a by-product of the proposed algorithm.

1.3. Restricted Latent Class Model

Our second example is restricted latent class models which are also widely used in social and behavioral sciences. For example, they are commonly used in education for cognitive diagnosis (von Davier & Lee, 2019). These models differ from IFA models in that they assume discrete latent variables. Here, we consider a setting for cognitive diagnosis when both data and latent variables are binary. Consider data taking the same form as that for IFA, denoted by $Y={\left(Y_{\mathrm{ij}}\right)}_{N\times J}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Y}}= (Y_{ij})_{N\times J}$$\end{document} and $y={\left(y_{\mathrm{ij}}\right)}_{N\times J}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}= (y_{ij})_{N\times J}$$\end{document} . In this context, $Y_{\mathrm{ij}}=1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{ij} = 1$$\end{document} means that item j is answered correctly and $Y_{\mathrm{ij}}=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{ij} = 0$$\end{document} means an incorrect answer.

The restricted latent class model assumes that each individual is characterized by a K-dimensional latent vector ${\xi }_i={({\xi }_{i1},\dots ,{\xi }_{\mathrm{iK}})}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i = (\xi _{i1},\ldots , \xi _{iK})^\top $$\end{document} , $i=1,\dots ,N$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1,\ldots , N$$\end{document} , where ${\xi }_{\mathrm{ik}}\in \{0,1\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\xi _{ik} \in \{0, 1\}$$\end{document} . Thus, the latent variables are $\xi ={\left({\xi }_{\mathrm{ik}}\right)}_{N\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi } = (\xi _{ik})_{N\times K}$$\end{document} , whose state space $\Xi ={\{0,1\}}^{N\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Xi = \{0,1\}^{N\times K}$$\end{document} contains all $N\times K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N\times K$$\end{document} binary matrices. Each dimension of ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} represents a skill, and ${\xi }_{\mathrm{ik}}=1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\xi _{ik} = 1$$\end{document} indicates that person i has mastered the kth skill and ${\xi }_{\mathrm{ik}}=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\xi _{ik} = 0$$\end{document} otherwise.

The restricted latent class model can be parametrized as follows.

1. 1. The person-specific latent vectors ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} , $i=1,\dots ,N$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1,\ldots , N$$\end{document} , are i.i.d., following a categorical distribution satisfying

   $\begin{array}{c}P({\xi }_i=\alpha )=\frac{exp\left({\nu }_{\alpha }\right)}{{\sum }_{{\alpha }^{\prime }\in {\{0,1\}}^K}exp\left({\nu }_{{\alpha }^{\prime }}\right)},\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbb {P}}(\varvec{\xi }_i = \varvec{\alpha }) = \frac{\exp (\nu _{\varvec{\alpha }})}{\sum _{\varvec{\alpha }' \in \{0, 1\}^K}\exp (\nu _{\varvec{\alpha }'})}, \end{aligned}$$\end{document}
   where $\alpha \in {\{0,1\}}^K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha } \in \{0, 1\}^K$$\end{document} represents an attribute profile representing the mastery status on all K attributes, and we set ${\nu }_{{\alpha }^{\prime }}=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu _{\varvec{\alpha }'} = 0$$\end{document} as the baseline, for ${\alpha }^{\prime }={(0,\dots ,0)}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }' = (0,\ldots , 0)^\top $$\end{document} .

2. 2. $Y_{\mathrm{ij}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{ij}$$\end{document} given ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} follows a Bernoulli distribution, satisfying

   $\begin{array}{c}P(Y_{\mathrm{ij}}=1\mid {\xi }_i=\alpha )=\frac{exp\left({\theta }_{j,\alpha }\right)}{1+exp\left({\theta }_{j,\alpha }\right)},\alpha \in {\{0,1\}}^K.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i = \varvec{\alpha }) = \frac{\exp (\theta _{j, \varvec{\alpha }})}{1+\exp (\theta _{j, \varvec{\alpha }})}, ~\varvec{\alpha } \in \{0, 1\}^K. \end{aligned}$$\end{document}

3. 3. Local independence is still assumed. That is, $Y_{i1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{i1}$$\end{document} ,..., $Y_{\mathrm{iJ}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$Y_{iJ}$$\end{document} are conditionally independent given ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} .

The above model specification leads to a marginal likelihood function

(9) $\begin{array}{c}l\left(\beta \right)=\sum _{i=1}^{N}log\left(\sum _{\alpha \in {\{0,1\}}^K},\frac{exp\left({\nu }_{\alpha }\right)}{{\sum }_{{\alpha }^{\prime }\in {\{0,1\}}^K}exp\left({\nu }_{{\alpha }^{\prime }}\right)},\prod _{j=1}^J,\frac{exp\left(y_{\mathrm{ij}}{\theta }_{j,\alpha }\right)}{1+exp\left({\theta }_{j,\alpha }\right)}\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} l({\varvec{\beta }}) = \sum _{i=1}^N \log \left( \sum _{\varvec{\alpha }\in \{0,1\}^K} \frac{\exp (\nu _{\varvec{\alpha }})}{\sum _{\varvec{\alpha }' \in \{0, 1\}^K}\exp (\nu _{\varvec{\alpha }'})} \prod _{j=1}^J \frac{\exp (y_{ij}\theta _{j, \varvec{\alpha }})}{1+\exp (\theta _{j, \varvec{\alpha }})} \right) , \end{aligned}$$\end{document}

where $\beta =\{{\nu }_{\alpha },{\theta }_{j,\alpha },\alpha \in {\{0,1\}}^K,j=1,\dots ,J\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}= \{\nu _{\varvec{\alpha }}, \theta _{j, \varvec{\alpha }}, \varvec{\alpha }\in \{0, 1\}^K, j = 1,\ldots , J\}$$\end{document} .

We consider a confirmatory setting where there exists a design matrix, similar to the $Q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}}$$\end{document} -matrix in confirmatory IFA. With slight abuse of notation, we still denote $Q={\left(q_{\mathrm{jk}}\right)}_{J\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}} = (q_{jk})_{J\times K}$$\end{document} , where $q_{\mathrm{jk}}\in \{0,1\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$q_{jk} \in \{0, 1\}$$\end{document} . Here, $q_{\mathrm{jk}}=1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$q_{jk} = 1$$\end{document} indicates that solving item j requires the kth skill and $q_{\mathrm{jk}}=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$q_{jk} = 0$$\end{document} otherwise. As will be explained below, this design matrix leads to equality and inequality constraints in model parameters.

Denote $q_j={(q_{j1},\dots ,q_{\mathrm{jK}})}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {q}}_{j} = (q_{j1},\ldots , q_{jK})^\top $$\end{document} as the design vector for item j. For $\alpha ={({\alpha }_1,\dots ,{\alpha }_K)}^{\top }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha } = (\alpha _1,\ldots , \alpha _K)^\top $$\end{document} , we write

$\begin{array}{c}\alpha ⪰q_j,\text{if}{\alpha }_k\ge q_{\mathrm{jk}}\text{for}\text{all}k\in \{1,\dots ,K\},\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\alpha } \succeq {\mathbf {q}}_{j}, ~\text{ if }~ \alpha _k \ge q_{jk} ~\text{ for } \text{ all }~ k \in \{1,\ldots , K\}, \end{aligned}$$\end{document}

and write

$\begin{array}{c}\alpha ⪰̸q_j,\text{if}\text{there}\text{exists}k\text{such}\text{that}{\alpha }_k<q_{\mathrm{jk}}.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\alpha } \nsucceq {\mathbf {q}}_{j}, ~\text{ if } \text{ there } \text{ exists } k \text{ such } \text{ that }~ \alpha _k < q_{jk}. \end{aligned}$$\end{document}

That is, $\alpha ⪰q_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha } \succeq {\mathbf {q}}_{j}$$\end{document} if profile $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }$$\end{document} has all the skills needed for solving item j and $\alpha ⪰̸q_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha } \nsucceq {\mathbf {q}}_{j}$$\end{document} if not. The design information leads to the following constraints:

1. 1. $P(Y_{\mathrm{ij}}=1\mid {\xi }_i=\alpha )=P(Y_{\mathrm{ij}}=1\mid {\xi }_i={\alpha }^{\prime })$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i = \varvec{\alpha }) = {\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i = \varvec{\alpha }')$$\end{document} , if both $\alpha ,{\alpha }^{\prime }⪰q_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }, \varvec{\alpha }' \succeq {\mathbf {q}}_{j}$$\end{document} . That is, individuals who have mastered all the required skills have the same chance of answering the item correctly.

2. 2. $P(Y_{\mathrm{ij}}=1\mid {\xi }_i=\alpha )\ge P(Y_{\mathrm{ij}}=1\mid {\xi }_i={\alpha }^{\prime })$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i = \varvec{\alpha }) \ge {\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i = \varvec{\alpha }')$$\end{document} if $\alpha ⪰q_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }\succeq {\mathbf {q}}_{j}$$\end{document} and ${\alpha }^{\prime }⪰̸q_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }'\nsucceq {\mathbf {q}}_{j}$$\end{document} . That is, students who have mastered all the required skills have a higher chance of answering the item correctly than those who do not.

3. 3. $P(Y_{\mathrm{ij}}=1\mid {\xi }_i=\alpha )\ge P(Y_{\mathrm{ij}}=1\mid {\xi }_i=0)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {P}}(Y_{ij} = 1\mid \varvec{\xi }_i=\varvec{\alpha })\ge {\mathbb {P}}(Y_{ij}=1\mid \varvec{\xi }_i={\varvec{0}})$$\end{document} for all $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }$$\end{document} . That is, students who have not mastered any skill have the lowest chance of answering correctly.

We refer the readers to Xu (2017) for more discussions on these constraints which are key to the identification of this model. Under these constraints, the MML estimator is given by

(10) $\begin{array}{c}\begin{array}{cc}\hat{\beta }=& \underset{\beta \in B}{\text{arg min}}-l\left(\beta \right),\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned} {\hat{{\varvec{\beta }}}} =&\mathop {\hbox {arg min}}\limits _{{\varvec{\beta }}\in {\mathcal {B}}} - l({\varvec{\beta }}), \end{aligned} \end{aligned}$$\end{document}

where

$\begin{array}{c}\begin{array}{c}B=\{\beta :\underset{\alpha ⪰q_j}{max}{\theta }_{j,\alpha }=\underset{\alpha ⪰q_j}{min}{\theta }_{j,\alpha }\ge {\theta }_{j,{\alpha }^{\prime }}\ge {\theta }_{j,0},\text{if}{\alpha }^{\prime }⪰̸q_j,{\nu }_0=0\}.\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned} {\mathcal {B}} = \{{\varvec{\beta }}: \max _{\varvec{\alpha }\succeq {\mathbf {q}}_{j}} \theta _{j, \varvec{\alpha }} = \min _{\varvec{\alpha }\succeq {\mathbf {q}}_{j}} \theta _{j, \varvec{\alpha }} \ge \theta _{j, \varvec{\alpha }'} \ge \theta _{j,{\varvec{0}}}, ~\text{ if }~ \varvec{\alpha }^\prime \nsucceq \mathbf {q}_j, \,\, \nu _{{\varvec{0}}} = 0\}. \end{aligned} \end{aligned}$$\end{document}

When K is relatively large, the computation for solving (10) becomes challenging, due to both the summation over $2^K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$2^K$$\end{document} possible values of $\alpha$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha }$$\end{document} in $l\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l({\varvec{\beta }})$$\end{document} , and the large number of inequality constraints.

2. Stochastic Proximal Algorithm

In this section, we propose a quasi-Newton stochastic proximal algorithm for the computation of (2). The description in this section will focus on the computation aspect, without emphasizing the regularity conditions needed for its convergence. A rigorous theoretical treatment will be given in Sect. 3. In what follows, we describe the algorithm in its general form in Sect. 2.1, followed by details for two specific models in Sects. 2.2 and 2.3, and finally comparisons with related algorithms in Sect. 2.4.

2.1. General Algorithm

For ease of exposition, we introduce some new notation. We write the penalty function as the sum of two terms, $R\left(\beta \right)=R_1\left(\beta \right)+R_2\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R({\varvec{\beta }}) = R_1({\varvec{\beta }}) + R_2({\varvec{\beta }})$$\end{document} , where $R_1\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_1({\varvec{\beta }})$$\end{document} is a smooth function and $R_2\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_2({\varvec{\beta }})$$\end{document} is non-smooth. In the example of regularized estimation for exploratory IFA, $R_1\left(\beta \right)=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_1({\varvec{\beta }}) = 0$$\end{document} and $R_2\left(\beta \right)=\lambda {\sum }_{j=1}^J{\sum }_{k=1}^K\left|a_{\mathrm{jk}}\right|$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_2({\varvec{\beta }}) = \lambda \sum _{j = 1}^J\sum _{k=1}^K \vert a_{jk}\vert $$\end{document} , when $R\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R({\varvec{\beta }})$$\end{document} is an $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} penalty as in (6). When an elastic net penalty is used as in (7), $R_1\left(\beta \right)={\lambda }_1{\sum }_{j=1}^J{\sum }_{k=1}^K{a}_{\mathrm{jk}}^2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_1({\varvec{\beta }}) = \lambda _1 \sum _{j = 1}^J\sum _{k=1}^K a_{jk}^2 $$\end{document} and $R_2\left(\beta \right)={\lambda }_2{\sum }_{j=1}^J{\sum }_{k=1}^K\left|a_{\mathrm{jk}}\right|$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_2({\varvec{\beta }}) = \lambda _2 \sum _{j = 1}^J\sum _{k=1}^K \vert a_{jk}\vert $$\end{document} .

The optimization problem can be reexpressed as

(11) $\begin{array}{c}\underset{\beta }{min}h\left(\beta \right)+g\left(\beta \right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \min _{\varvec{\beta }}~ h({\varvec{\beta }}) + g({\varvec{\beta }}), \end{aligned}$$\end{document}

where $h\left(\beta \right)=-l\left(\beta \right)+R_1\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }}) = -l(\varvec{\beta }) + R_1(\varvec{\beta })$$\end{document} and $g:R^p\to R\cup \{+\infty \}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g: {\mathbb {R}}^{p} \rightarrow {\mathbb {R}} \cup \{+\infty \}$$\end{document} is a generalized function taking the form $g\left(\beta \right)=R_2\left(\beta \right)+I_B\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }}) = R_2({\varvec{\beta }}) + I_{{\mathcal {B}}}({\varvec{\beta }})$$\end{document} , where

(12) $\begin{array}{c}{I}_B\left(\beta \right)=\left(\begin{array}{cc}0,& \text{if}\beta \in B,\\ +\infty ,& \text{otherwise.}\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} I_{{\mathcal {B}}}({\varvec{\beta }}) = \left\{ \begin{array}{ll} 0, &{} ~\text{ if }~ {\varvec{\beta }}\in {\mathcal {B}},\\ +\infty , &{} ~\text{ otherwise. }~ \end{array}\right. \end{aligned}$$\end{document}

Note that since both $l\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l({\varvec{\beta }})$$\end{document} and $R_1\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_1({\varvec{\beta }})$$\end{document} are smooth in $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} , $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} is still smooth in $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} . The second term $g\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }})$$\end{document} is non-smooth in $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} , unless it is degenerate (i.e., $g\left(\beta \right)\equiv 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }}) \equiv 0$$\end{document} ). We further write

(13) $\begin{array}{c}H(\xi ,\beta )=-logf(y,\xi \mid \beta )+R_1\left(\beta \right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} H(\varvec{\xi }, {\varvec{\beta }}) = - \log f({\mathbf {y}}, \varvec{\xi }\mid {\varvec{\beta }}) + R_1({\varvec{\beta }}), \end{aligned}$$\end{document}

which can be viewed as a complete-data version of $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} that will be used in the algorithm.

The algorithm relies on a scaled proximal operator (Lee et al., 2014) for the g function, defined as

$\begin{array}{c}{\text{Prox}}_{\gamma ,g}^D\left(\beta \right)=\underset{x\in R^p}{\text{arg min}}\left(g,\left(x\right),+,\frac{1}{2\gamma },{\Vert x-\beta \Vert }_D^2\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \text {Prox}_{\gamma , g}^{{\mathbf {D}}}({\varvec{\beta }}) = \mathop {\hbox {arg min}}\limits _{{\mathbf {x}}\in {\mathbb {R}}^{p}} \left\{ g({\mathbf {x}}) + \frac{1}{2\gamma } \Vert {\mathbf {x}}- {\varvec{\beta }}\Vert ^2_{{\mathbf {D}}} \right\} , \end{aligned}$$\end{document}

where $\gamma >0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma >0$$\end{document} , $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}$$\end{document} is a strictly positive definite matrix, and ${\Vert ·\Vert }_D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert \cdot \Vert _{{\mathbf {D}}}$$\end{document} is a norm defined by ${\Vert x\Vert }_D^2={⟨x,x⟩}_D=x^{\top }Dx$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert {\mathbf {x}}\Vert _{{\mathbf {D}}}^2 = \langle {\varvec{x}},{\varvec{x}}\rangle _{{\mathbf {D}}}= {\mathbf {x}}^\top {\mathbf {D}} {\mathbf {x}}$$\end{document} . The choices of $\gamma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma $$\end{document} , $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}$$\end{document} , and the intuition behind the proximal operator will be explained in the sequel.

Our general algorithm is described in Algorithm 1, followed by implementation details. The proposed algorithm is an extension of a perturbed proximal gradient algorithm (Atchadé et al., 2017). The major difference is that the proposed algorithm makes use of second-order information from the smooth part of the objective function, which can substantially speed up its convergence. See Sect. 2.4 for further comparison.

Algorithm 1

(Stochastic Proximal Algorithm)

• Input Data $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}$$\end{document} , initial parameters ${\beta }^{\left(0\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(0)}$$\end{document} , a sequence of step size ${\gamma }_s,s=1,2,\dots$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _s, s= 1, 2,\ldots $$\end{document} , pre-specified tuning parameters $c_2\ge c_1>0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_{2} \ge c_{1} > 0$$\end{document} , and burn-in size $\varpi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varpi $$\end{document} .

• Update At tth iteration where $t\ge 1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t \ge 1$$\end{document} , we perform the following two steps:

  1. 1. Stochastic step Sample $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} from the conditional distribution of $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} given $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}$$\end{document} ,

     $\begin{array}{c}\psi \left(\xi \right)\propto f(y,\xi \mid {\beta }^{(t-1)}),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \psi (\varvec{\xi }) \propto f({\mathbf {y}}, \varvec{\xi }\mid {\varvec{\beta }}^{(t-1)}), \end{aligned}$$\end{document}
     and obtain ${\xi }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }^{(t)}$$\end{document} . The sampling can be either exact or approximated by MCMC.

  2. 2. Proximal step Update model parameters by

     (14) $\begin{array}{c}{\beta }^{\left(t\right)}={\text{Prox}}_{{\gamma }_t,g}^{D^{\left(t\right)}}({\beta }^{(t-1)}-{\gamma }_t{\left(D^{\left(t\right)}\right)}^{-1}{G}^{\left(t\right)}),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\varvec{\beta }}^{(t)} = \text {Prox}_{\gamma _t, g}^{{\mathbf {D}}^{(t)}}\big ({\varvec{\beta }}^{(t-1)} - \gamma _t({\mathbf {D}}^{(t)})^{-1}{\mathbf {G}}^{(t)}\big ), \end{aligned}$$\end{document}
     where
     $\begin{array}{c}{G}^{\left(t\right)}=\frac{\partial H({\xi }^{\left(t\right)},\beta )}{\partial \beta }{|}_{\beta ={\beta }^{(t-1)}}.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbf {G}}^{(t)} = \frac{\partial H(\varvec{\xi }^{(t)}, {\varvec{\beta }})}{\partial {\varvec{\beta }}} \bigg \vert _{{\varvec{\beta }}= {\varvec{\beta }}^{(t-1)}}. \end{aligned}$$\end{document}
     $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} is a diagonal matrix with diagonal entries
     $\begin{array}{c}{\delta }_i^{\left(t\right)}=\frac{t-1}{t}{\delta }_i^{(t-1)}+\frac{1}{t}T\left({\tilde{\delta }}_i^{\left(t\right)},;,c_1,,,c_2\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \delta _i^{(t)} = \frac{t-1}{t} \delta _i^{(t-1)} + \frac{1}{t}T\left( {\tilde{\delta }}_i^{(t)}; c_{1}, c_{2}\right) , \end{aligned}$$\end{document}
     where $T(x;c_1,c_2)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$T(x;c_1,c_2)$$\end{document} is a truncation function defined as
     (15) $\begin{array}{c}T(x;c_1,c_2)=\left(\begin{array}{cc}{c}_1,& \text{if}x<c_1,\\ x,& \text{if}x\in [c_1,c_2],\\ c_2,& \text{if}x>c_2.\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} T(x;c_1,c_2) = \left\{ \begin{array}{ll} c_1, &{} ~\text{ if }~ x < c_1, \\ x, &{} ~\text{ if }~ x \in [c_1,c_2],\\ c_2, &{} ~\text{ if }~ x > c_2. \end{array}\right. \end{aligned}$$\end{document}
     Here ${\tilde{\delta }}_i^{\left(t\right)}={\tilde{\delta }}_{i,1}^{\left(t\right)}+{\left({\tilde{\delta }}_{i,2}^{\left(t\right)}\right)}^2,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{\delta }}_i^{(t)} = {\tilde{\delta }}_{i,1}^{(t)} + ({\tilde{\delta }}_{i,2}^{(t)})^2,$$\end{document} where
     $\begin{array}{cc}{\tilde{\delta }}_{i,1}^{\left(t\right)}=& (1-{\gamma }_t){\tilde{\delta }}_{i,1}^{(t-1)}+{\gamma }_t\left(\frac{{\partial }^{2}H({\xi }^{\left(t\right)},\beta )}{\partial {\beta }_i^2},{|}_{\beta ={\beta }^{(t-1)}},,-,,{\left(\frac{\partial H({\xi }^{\left(t\right)},\beta )}{\partial {\beta }_i}\right)}^2,{|}_{\beta ={\beta }^{(t-1)}}\right),\\ {\tilde{\delta }}_{i,2}^{\left(t\right)}=& (1-{\gamma }_t){\tilde{\delta }}_{i,2}^{(t-1)}+{\gamma }_t{\left(\frac{\partial H({\xi }^{\left(t\right)},\beta )}{\partial {\beta }_i}\right)}^2{|}_{\beta ={\beta }^{(t-1)}}.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\tilde{\delta }}_{i,1}^{(t)}= & {} (1\!-\!\gamma _t){\tilde{\delta }}_{i,1}^{(t-1)} \!+\! \gamma _t \left( \frac{\partial ^2 H(\varvec{\xi }^{(t)}, {\varvec{\beta }})}{\partial \beta _i^2}\bigg \vert _{{\varvec{\beta }}= {\varvec{\beta }}^{(t-1)}} \!-\! \left( \frac{\partial H(\varvec{\xi }^{(t)},{\varvec{\beta }})}{\partial \beta _i}\right) ^2\bigg \vert _{{\varvec{\beta }}={\varvec{\beta }}^{(t-1)}}\right) ,\\ {\tilde{\delta }}_{i,2}^{(t)}= & {} (1-\gamma _t){\tilde{\delta }}_{i,2}^{(t-1)} + \gamma _t \left( \frac{\partial H(\varvec{\xi }^{(t)},{\varvec{\beta }})}{\partial \beta _i}\right) ^2 \bigg \vert _{{\varvec{\beta }}={\varvec{\beta }}^{(t-1)}}. \end{aligned}$$\end{document}
  Iteratively perform these two steps until a stopping criterion is satisfied and let n be the last iteration number.

• Output ${\overline{\beta }}_n={\sum }_{t=\varpi +1}^n{\beta }^{\left(t\right)}/(n-\varpi )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n = {\sum _{t=\varpi +1}^n {\varvec{\beta }}^{(t)}}/(n-\varpi )$$\end{document} .

In what follows, we make a few remarks to provide some intuitions about the algorithm.

Remark 1

(Connection with stochastic gradient descent) To provide some intuitions about the proposed method, we first make a connection between the proposed method and the stochastic gradient descent (SGD) algorithm. In fact, when the sampling of $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} is exact in the stochastic step, then $G^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}^{(t)}$$\end{document} is a stochastic gradient of the smooth part of our objective function, in the sense that $E(G^{\left(t\right)}\mid y,{\beta }^{(t-1)})=\nabla h\left(\beta \right){|}_{\beta ={\beta }^{(t-1)}}.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {E}}({\mathbf {G}}^{(t)}\mid {\mathbf {y}}, {\varvec{\beta }}^{(t-1)}) = \nabla h({\varvec{\beta }})\vert _{{\varvec{\beta }}= {\varvec{\beta }}^{(t-1)}}.$$\end{document} If, in addition, there is no constraint or non-smooth penalty, i.e., $g\left(\beta \right)\equiv 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }}) \equiv 0$$\end{document} , then the proximal step degenerates to an SGD update ${\beta }^{\left(t\right)}={\beta }^{(t-1)}-{\gamma }_t{\left(D^{\left(t\right)}\right)}^{-1}{G}^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)} = {\varvec{\beta }}^{(t-1)} - \gamma _t({\mathbf {D}}^{(t)})^{-1}{\mathbf {G}}^{(t)}$$\end{document} . In that case, the proposed method becomes a version of SGD.

Remark 2

(Proximal step) We provide some intuitions about the proximal step. We start with two special cases. First, as mentioned in Remark 1, if there is no constraint or non-smooth penalty, then the proximal step is nothing but a stochastic gradient descent step. This is because, the scaled proximal operator degenerates to an identity map, i.e., ${\text{Prox}}_{\gamma ,g}^D\left(\beta \right)=\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\text {Prox}_{\gamma , g}^{{\mathbf {D}}}({\varvec{\beta }}) = {\varvec{\beta }}$$\end{document} . Second, when the g function involves constraints but does not contain a non-smooth penalty, then the proximal step is a projected stochastic gradient descent step. That is, one first performs a stochastic gradient descent update ${\tilde{\beta }}^{\left(t\right)}={\beta }^{(t-1)}-{\gamma }_t{\left(D^{\left(t\right)}\right)}^{-1}{G}^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{{\varvec{\beta }}}}^{(t)} = {\varvec{\beta }}^{(t-1)} - \gamma _t({\mathbf {D}}^{(t)})^{-1}{\mathbf {G}}^{(t)}$$\end{document} . Then ${\tilde{\beta }}^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{{\varvec{\beta }}}}^{(t)}$$\end{document} is projected back to the feasible region $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}$$\end{document} by the scaled proximal operator:

$\begin{array}{c}\hat{\beta }=\underset{\beta \in B}{\text{arg min}}{\Vert \beta -{\tilde{\beta }}^{\left(t\right)}\Vert }_D,\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\hat{{\varvec{\beta }}}} = \mathop {\hbox {arg min}}\limits _{{\varvec{\beta }}\in {\mathcal {B}}} \Vert {\varvec{\beta }}- {\tilde{{\varvec{\beta }}}}^{(t)}\Vert _{{\mathbf {D}}}, \end{aligned}$$\end{document}

which is a projection under the norm ${\Vert ·\Vert }_D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert \cdot \Vert _{{\mathbf {D}}}$$\end{document} . When $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}$$\end{document} is an identity matrix as in the vanilla (i.e., non-scaled) proximal operator, then the projection is based on the Euclidean distance.

More generally, when the g function involves non-smooth penalties, then the proximal step can be viewed as minimizing the sum of $g\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }})$$\end{document} and a quadratic approximation of $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} at ${\beta }^{(t-1)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t-1)}$$\end{document} ; see Lee et al. (2014) for more explanations. We provide an example to facilitate the understanding. Suppose that

$\begin{array}{c}g\left(\beta \right)=\lambda \sum _{i=1}^p\left|{\beta }_i\right|\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} g({\varvec{\beta }}) = \lambda \sum _{i=1}^p \vert \beta _i\vert \end{aligned}$$\end{document}

is the Lasso penalty, and $D=diag({\delta }_1,\dots ,{\delta }_p)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}} = diag(\delta _1,\ldots , \delta _p)$$\end{document} is a diagonal matrix, where $\lambda ,{\delta }_i>0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda , \delta _i > 0$$\end{document} , $i=1,\dots ,p$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1,\ldots , p$$\end{document} . Then ${\text{Prox}}_{\gamma ,g}^D\left({\tilde{\beta }}^{\left(t\right)}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\text {Prox}_{\gamma , g}^{{\mathbf {D}}}({\tilde{{\varvec{\beta }}}}^{(t)})$$\end{document} involves solving p optimization problems separately, each of which takes the form

(16) $\begin{array}{c}{\hat{\beta }}_i=\underset{x}{\text{arg min}}\frac{1}{2}{(x-{\tilde{\beta }}_i^{\left(t\right)})}^2+\frac{\lambda \gamma }{{\delta }_i}\left|x\right|.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\hat{\beta }}_i = \mathop {\hbox {arg min}}\limits _{x} \frac{1}{2} (x - {\tilde{\beta }}_i^{(t)})^2 + \frac{\lambda \gamma }{\delta _i} |x|. \end{aligned}$$\end{document}

It is well known that (16) has a closed-form solution given by soft-thresholding (see Chapter 3, Friedman et al., 2001):

$\begin{array}{c}{\hat{\beta }}_i=\left(\begin{array}{cc}{\tilde{\beta }}_i^{\left(t\right)}-\frac{\lambda \gamma }{{\delta }_i},& \text{if}{\tilde{\beta }}_i^{\left(t\right)}>\frac{\lambda \gamma }{{\delta }_i},\\ {\tilde{\beta }}_i^{\left(t\right)}+\frac{\lambda \gamma }{{\delta }_i},& \text{if}{\tilde{\beta }}_i^{\left(t\right)}<-\frac{\lambda \gamma }{{\delta }_i},\\ 0,& \text{otherwise}.\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\hat{\beta }}_i = \left\{ \begin{array}{ll} {\tilde{\beta }}_i^{(t)} - \frac{\lambda \gamma }{\delta _i}, &{} ~\text{ if }~ {\tilde{\beta }}_i^{(t)} > \frac{\lambda \gamma }{\delta _i}, \\ {\tilde{\beta }}_i^{(t)} + \frac{\lambda \gamma }{\delta _i}, &{} ~\text{ if }~ {\tilde{\beta }}_i^{(t)} < - \frac{\lambda \gamma }{\delta _i}, \\ 0, &{} ~\text{ otherwise }. \end{array}\right. \end{aligned}$$\end{document}

Remark 3

(Role of $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} ) Our proximal step is a quasi-Newton proximal update proposed in Lee et al. (2014) under a non-stochastic optimization setting. As shown in Lee et al. (2014), quasi-Newton proximal methods converge faster than first-order proximal methods under the non-stochastic setting. Here, the diagonal matrix $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} is used to approximate the Hessian matrix of $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} at ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} . When ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} converges to $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} , then ${\delta }_i^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta _i^{(t)}$$\end{document} , the ith diagonal term of $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} converges to $T\left(\frac{{\partial }^{2}h}{\partial {\beta }_i^2},{|}_{\beta =\hat{\beta }},;,c_1,,,c_2\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$T\left( \frac{\partial ^2 h}{\partial \beta _i^2}\vert _{{\varvec{\beta }}= {\hat{{\varvec{\beta }}}}}; c_{1}, c_{2}\right) $$\end{document} where T is the truncation function defined in (15); see Remark 8 for more explanations.

In the proposed update, we choose $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} to be a diagonal matrix for computational convenience. Specifically, as discussed in Remark 2, the proximal step is in a closed form when $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} is a diagonal matrix. In addition, the proximal step requires to calculate the inverse of $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} , whose complexity is much lower when $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} is diagonal.

We point out that using a diagonal matrix to approximate the Hessian matrix is a popular and effective trick in numerical optimization (e.g., Chapter 5, Bertsekas et al., 1992; Becker & Le Cun, 1988), especially for large-scale optimization problems. In principle, it is possible to allow $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} to be non-diagonal. In fact, it is not difficult to generalize the BFGS updating formula for $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} given in Lee et al. (2014) to a stochastic version.

Our choice of $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} guarantees its eigenvalues to be constrained in the interval $[c_1,c_2]$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$[c_{1}, c_{2}]$$\end{document} . It rules out the singular situation when $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} is not strictly positive definite. In the implementation, we set $c_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_{1}$$\end{document} ’s to be a sufficiently small constant and set $c_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_{2}$$\end{document} ’s to be a sufficiently large constant. According to simulation, the algorithm tends to be insensitive to these choices.

We further provide some remarks regarding the implementation details.

Remark 4

(Choices of step size) As will be shown in Sect. 3, the convergence of the proposed method requires the step size to satisfy ${\sum }_{t=1}^{\infty }{\gamma }_t=\infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{t=1}^{\infty } \gamma _t = \infty $$\end{document} and ${\sum }_{t=1}^{\infty }{\gamma }_t^2<\infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{t=1}^{\infty } \gamma _t^2 < \infty $$\end{document} . This requirement is also needed in the Robbins–Monro algorithm. Here, we choose the step size ${\gamma }_t=\mu t^{-\frac{1}{2}-\epsilon }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t = \mu t^{-\frac{1}{2} - \varepsilon }$$\end{document} so that the above requirement is satisfied, where $\mu$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mu $$\end{document} is a positive constant and $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} is a small positive constant. As will be shown in Sect. 3, with sufficiently small $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} , ${\overline{\beta }}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n$$\end{document} is almost optimal in terms of its convergence speed. We point out that $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} is needed to prove the convergence of ${\overline{\beta }}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n$$\end{document} , under our non-convex setting. It is not needed, if the objective function (2) is convex; see Atchadé et al. (2017). The requirement of $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} may be an artifact due to our proof strategy. Simulation results show that the algorithm converges well even if we set $\epsilon =0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon = 0$$\end{document} . For the numerical analysis in this paper, we set $\epsilon ={10}^{-2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon = 10^{-2}$$\end{document} .

We point out that our choice of step size is very different from the step size in the Robbins–Monro algorithm, for which asymptotic results (Fabian, 1968) suggest that the optimal choice of step size satisfies ${\gamma }_t=O(1/t)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t = O(1/t)$$\end{document} .

Remark 5

(Starting point) As the objective function (2) is typically non-convex for most latent variable models, the choice of the starting point ${\beta }^{\left(0\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(0)}$$\end{document} matters. The algorithm is more likely to converge to the global optimum given a good starting point. One strategy is to run the proposed algorithm with multiple random starting points and then choose the best-fitting solution. Alternatively, one may find a good starting point using less accurate but computationally faster estimators, such as the constrained joint maximum likelihood estimator (Chen et al., 2019; 2020) or spectral methods (H. Zhang et al., 2020a). Moreover, to further avoid convergence to local optima, one may also use multiple random starting points and choose the one with the smallest objective function value.

Remark 6

(Sampling in stochastic step) As mentioned in Remark 1, when the latent variables $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} can be sampled exactly in the stochastic step, then $G^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}^{(t)}$$\end{document} is a stochastic gradient of $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} . Unfortunately, exact sampling is only possible under some situations such as restricted latent class analysis. In most cases, we only have approximate samples from an MCMC algorithm. For example, as discussed below, the latent variables in IFA can be sampled by a block-wise Gibbs sampler. With approximate samples, $G^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}^{(t)}$$\end{document} is only approximately unbiased. As we show in Sect. 3, such $G^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}^{(t)}$$\end{document} may still yield convergence of ${\overline{\beta }}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n$$\end{document} .

Remark 7

(Stopping criterion) In the implementation of Algorithm 1, we stop the iterative update by monitoring a window of successive differences in ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} . More precisely, we stop the iteration if all differences in the window are less than a given threshold. Unless otherwise stated, the numerical analysis in this paper uses a window size 3. The same stopping criterion is also adopted by the Metropolis–Hasting Robins–Monro algorithm proposed by Cai (2010a).

Finally, as we explain in Remark 8, certain quantities, including the Hessian matrix of $l\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l({\varvec{\beta }})$$\end{document} , can be obtained as a by-product of the proposed algorithm.

Remark 8

(By-product) It is often of interest to compute quantities of the form

(17) $\begin{array}{c}\hat{M}=E\left(m,(,y,,,\xi ,\mid ,\beta ,),\mid ,y,,,\beta \right){|}_{\beta =\hat{\beta }},\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\hat{M}} = {\mathbb {E}}\left[ m({\mathbf {y}}, \varvec{\xi }\mid {\varvec{\beta }})\mid {\mathbf {y}}, {\varvec{\beta }}\right] \big \vert _{{\varvec{\beta }}= {\hat{{\varvec{\beta }}}}}, \end{aligned}$$\end{document}

where $m(y,\xi \mid \beta )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m({\mathbf {y}}, \varvec{\xi }\mid {\varvec{\beta }})$$\end{document} is a given function with an analytic form and the conditional expectation $E\left(·,\mid ,y,,,\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {E}}\left[ \cdot \mid {\mathbf {y}}, {\varvec{\beta }}\right] $$\end{document} is with respect to the conditional distribution of $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} given $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {y}}$$\end{document} . The quantity (17) is intractable due to the high-dimensional integral with respect to $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} . One such example is the Hessian matrix of $l\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l(\varvec{\beta })$$\end{document} at $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} as discussed in Sect. 1.2 that is a key quantity for the statistical inference of $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} . In fact, by Louis’ formula (Louis, 1982),

$\begin{array}{c}\begin{array}{cc}\frac{{\partial }^{2}l\left(\beta \right)}{\partial \beta \partial {\beta }^{\top }}=& E\left(\left(\frac{{\partial }^2(logf(y,\xi \mid \beta ))}{\partial \beta \partial {\beta }^{\top }},+,\frac{\partial (logf(y,\xi \mid \beta \left)\right)}{\partial \beta },{\left(\frac{\partial (logf(y,\xi \mid \beta \left)\right)}{\partial \beta }\right)}^{\top }\right),y,,,\beta \right)\\ & -E\left(\left(\frac{\partial (logf(y,\xi \mid \beta \left)\right)}{\partial \beta }\right),y,,,\beta \right){\left(E,\left(\left(\frac{\partial (logf(y,\xi \mid \beta \left)\right))}{\partial \beta }\right),y,,,\beta \right)\right)}^{\top }.\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned} \frac{\partial ^{2} l(\varvec{\beta })}{\partial \varvec{\beta } \partial \varvec{\beta }^{\top }}=&{\mathbb {E}}\left[ \left. \frac{\partial ^{2} (\log f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta }))}{\partial {\varvec{\beta }}\partial {\varvec{\beta }}^{\top }} + \frac{\partial (\log f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta }))}{\partial {\varvec{\beta }}}\left[ \frac{\partial (\log f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta }))}{\partial {\varvec{\beta }}}\right] ^{\top } \right| {\mathbf {y}}, {\varvec{\beta }}\right] \\&- {\mathbb {E}}\left[ \left. \frac{\partial (\log f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta }))}{\partial {\varvec{\beta }}} \right| {\mathbf {y}}, {\varvec{\beta }}\right] \left( {\mathbb {E}}\left[ \left. \frac{\partial (\log f({\mathbf {y}}, \varvec{\xi }\mid \varvec{\beta })))}{\partial {\varvec{\beta }}} \right| {\mathbf {y}}, {\varvec{\beta }}\right] \right) ^\top . \end{aligned} \end{aligned}$$\end{document}

The computation of (17) is a straightforward by-product of the proposed algorithm. To approximate $\hat{M}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{M}}$$\end{document} , we only need to add the following update in each iteration

(18) $\begin{array}{c}{M}^{\left(t\right)}=M^{(t-1)}+{\gamma }_t(m(y,{\xi }^{\left(t\right)}\mid {\beta }^{\left(t\right)})-M^{(t-1)}),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} M^{(t)} = M^{(t-1)} + \gamma _t\big (m({\mathbf {y}},\varvec{\xi }^{(t)}\mid {\varvec{\beta }}^{(t)})-M^{(t-1)}\big ), \end{aligned}$$\end{document}

for $t\ge 2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t \ge 2$$\end{document} , where $M^{\left(1\right)}=m(y,{\xi }^{\left(1\right)}\mid {\beta }^{\left(1\right)})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M^{(1)} = m({\mathbf {y}},\varvec{\xi }^{(1)}\mid {\varvec{\beta }}^{(1)})$$\end{document} . We approximate $\hat{M}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{M}}$$\end{document} by the Polyak–Ruppert averaging ${\overline{M}}_n=\left({\sum }_{t=\varpi +1}^n{M}^{\left(t\right)}\right)/(n-\varpi )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{M}}_n = (\sum _{t=\varpi + 1}^n M^{(t)})/(n-\varpi )$$\end{document} . When the sequence ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} converges to $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} (see Theorem 1 for the convergence analysis), under mild conditions, Theorem 3.17 of Ben-veniste et al. (1990) suggests the convergence of $M^{\left(n\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M^{(n)}$$\end{document} to $\hat{M}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{M}}$$\end{document} with probability 1, which further implies the convergence of ${\overline{M}}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{M}}_n$$\end{document} to $\hat{M}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{M}}$$\end{document} . Note that we use the averaged estimator ${\overline{M}}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{M}}_n$$\end{document} as it tends to converge faster than the pre-average sequence $M^{\left(n\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M^{(n)}$$\end{document} . We point out that the updating rule for the diagonal matrix $D^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}^{(t)}$$\end{document} in Algorithm 1 makes use of such an averaged estimator.

Remark 9

(Burn-in size) Like MCMC algorithms, the proposed method also has a burn-in period, where parameter updates from that period are not used in the Polyak–Ruppert averaging. The choice of the burn-in size will not affect the asymptotic property of the method, but does affect the empirical performance. This is because, the parameter updates may be far away from the solution due to the effect of the starting point. Including them in the Polyak–Ruppert averaging may introduce a high bias. In our numerical analysis, the burn-in size $\varpi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varpi $$\end{document} is fixed to be sufficiently large in each of our examples. Adaptive choice of the burn-in size is possible; see S. Zhang et al. (2020b).

2.2. Example I: Item Factor Analysis

We now explain the details of using the proposed method to solve (5) for exploratory IFA. The computation is similar when replacing the $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} regularization by the elastic net regularization. For confirmatory IFA, the stochastic step is the same as that of exploratory IFA and the proximal update step is straightforward as no penalty is involved. Therefore, the details for the computation of confirmatory IFA are omitted here.

We first consider the stochastic step for solving (5). Note that ${\xi }_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_1$$\end{document} ,..., ${\xi }_N$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_N$$\end{document} are conditionally independent given data, and thus can be sampled separately. For each ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} , we sample its entries by Gibbs sampling. More precisely, each entry is sampled by adaptive rejection sampling (Gilks & Wild, 1992; S. Zhang et al., 2020b), as the conditional distribution of ${\xi }_{\mathrm{ik}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\xi _{ik}$$\end{document} given data and the other entries of ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} is log-concave. We refer the readers to S. Zhang et al. (2020b) for more explanations of this sampling procedure. If a normal ogive IFA is considered instead of the logistic model above, then we can sample ${\xi }_i^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }^{(t)}_i$$\end{document} by a similar Gibbs method with a data augmentation trick; see Chen & Zhang (2020a) for a review.

We now discuss the computation for the proximal step. Recall that $\beta =\{B,d_j,a_j,j=1,\dots ,J\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}= \{{\mathbf {B}}, d_j, {\mathbf {a}}_j, j =1,\ldots , J\}$$\end{document} . We denote

$\begin{array}{c}{\tilde{\beta }}^{\left(t\right)}={\beta }^{(t-1)}-{\gamma }_t{\left(D^{\left(t\right)}\right)}^{-1}{G}^{\left(t\right)}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\tilde{{\varvec{\beta }}}}^{(t)} = {\varvec{\beta }}^{(t-1)} - \gamma _t({\mathbf {D}}^{(t)})^{-1}{\mathbf {G}}^{(t)} \end{aligned}$$\end{document}

as the input of the scaled proximal operator. The parameter update is given by

$\begin{array}{c}{\beta }^{\left(t\right)}={\text{Prox}}_{\gamma ,g}^D\left({\tilde{\beta }}^{\left(t\right)}\right)=\underset{\beta }{\text{arg min}}\left(g,\left(\beta \right),+,\frac{1}{2{\gamma }_t},\sum _{i=1}^p,{\delta }_i^{\left(t\right)},{({\beta }_i-{\tilde{\beta }}_i)}^2\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\varvec{\beta }}^{(t)} = \text {Prox}_{\gamma , g}^{{\mathbf {D}}}({\tilde{{\varvec{\beta }}}}^{(t)}) = \mathop {\hbox {arg min}}\limits _{{\varvec{\beta }}} \left\{ g({\varvec{\beta }}) + \frac{1}{2\gamma _t} \sum _{i=1}^p \delta _i^{(t)}(\beta _i - {\tilde{\beta }}_i )^2\right\} , \end{aligned}$$\end{document}

where the parameter space

$\begin{array}{c}B=\{\beta \in R^p:b_{k{k}^{\prime }}=0,1\le k<k^{\prime }\le K,\sum _{k^{\prime }=1}^K{b}_{k{k}^{\prime }}^2=1,k=1,\dots ,K\},\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathcal {B}} = \{{\varvec{\beta }}\in {\mathbb {R}}^{p}: b_{kk'} = 0, 1\le k < k' \le K, \sum _{k' = 1}^K b_{kk'}^2 = 1, k = 1,\ldots , K\}, \end{aligned}$$\end{document}

and $g\left(\beta \right)=\lambda {\sum }_{j=1}^J{\sum }_{k=1}^K\left|a_{\mathrm{jk}}\right|+I_B\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }}) = \lambda \sum _{j = 1}^J\sum _{k=1}^K \vert a_{jk}\vert + I_{{\mathcal {B}}}({\varvec{\beta }})$$\end{document} only involves loading parameters $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} and parameters $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} for the covariance matrix.

We first look at the update for $d_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j$$\end{document} s. As the g function does not involve $d_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j$$\end{document} , its update is simply $d_j^{\left(t\right)}={\tilde{d}}_j^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j^{(t)} = {\tilde{d}}_j^{(t)}$$\end{document} , where ${\tilde{d}}_j^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{d}}_j^{(t)}$$\end{document} is the corresponding component in ${\tilde{\beta }}^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{{\varvec{\beta }}}}^{(t)}$$\end{document} . We then look at the update for the loading parameters $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} . Suppose that $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} corresponds to the $i_{a_{\mathrm{jk}}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i_{a_{jk}}$$\end{document} th component of $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} . Then the update is given by solving the optimization

$\begin{array}{c}{a}_{\mathrm{jk}}^{\left(t\right)}=\underset{a_{\mathrm{jk}}}{\text{arg min}}\lambda \left|a_{\mathrm{jk}}\right|+\frac{1}{2{\gamma }_t}{\delta }_{i_{a_{\mathrm{jk}}}}^{\left(t\right)}{(a_{\mathrm{jk}}-{\tilde{a}}_{\mathrm{jk}}^{\left(t\right)})}^2.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} a_{jk}^{(t)} = \mathop {\hbox {arg min}}\limits _{a_{jk}}~ \lambda |a_{jk}| + \frac{1}{2\gamma _t} \delta _{i_{a_{jk}}}^{(t)} (a_{jk} - {\tilde{a}}_{jk}^{(t)})^2. \end{aligned}$$\end{document}

As discussed in Remark 2, this optimization has a closed-form solution via soft-thresholding. We finally look at the update for $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} . Suppose that $b_{\mathrm{kl}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$b_{kl}$$\end{document} corresponds to the $i_{b_{\mathrm{kl}}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i_{b_{kl}}$$\end{document} th component of $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} . Then the update of $b_k$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {b}}_k$$\end{document} , the kth row of $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} , is given by solving the following optimization problem:

$\begin{array}{c}{b}_k^{\left(t\right)}=\underset{b_k:\Vert b_k\Vert =1,b_{k{k}^{\prime }}=0,k^{\prime }>k}{\text{arg min}}\left(\sum _{l=1}^K,{\delta }_{i_{b_{\mathrm{kl}}}}^{\left(t\right)},{(b_{\mathrm{kl}}-{\tilde{b}}_{\mathrm{kl}}^{\left(t\right)})}^2\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbf {b}}_k^{(t)} = \mathop {\hbox {arg min}}\limits _{{\mathbf {b}}_k: \Vert {\mathbf {b}}_k\Vert =1, b_{kk'} = 0, k' > k} \left\{ \sum _{l=1}^K \delta _{i_{b_{kl}}}^{(t)} (b_{kl} - {\tilde{b}}_{kl}^{(t)})^2 \right\} , \end{aligned}$$\end{document}

which can be easily solved by the method of Lagrangian multiplier.

2.3. Example II: Restricted LCA

We now provide a brief discussion on the computation for the restricted LCA model. First, the stochastic step is straightforward, as the posterior distribution for each ${\xi }_i$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }_i$$\end{document} is still a categorical distribution which can be sampled exactly. Second, the proximal step requires to solve a quadratic programming problem. Again, we denote

$\begin{array}{c}{\tilde{\beta }}^{\left(t\right)}={\beta }^{(t-1)}-{\gamma }_t{\left(D^{\left(t\right)}\right)}^{-1}{G}^{\left(t\right)}.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\tilde{{\varvec{\beta }}}}^{(t)} = {\varvec{\beta }}^{(t-1)} - \gamma _t({\mathbf {D}}^{(t)})^{-1}{\mathbf {G}}^{(t)}. \end{aligned}$$\end{document}

The proximal step requires to solve the following quadratic programming problem

(19) $\begin{array}{c}\begin{array}{cc}& \underset{\beta }{min}{(\beta -{\tilde{\beta }}^{\left(t\right)})}^{\top }{D}^{\left(t\right)}(\beta -{\tilde{\beta }}^{\left(t\right)}),\\ s.t.& \underset{\alpha ⪰q_j}{max}{\theta }_{j,\alpha }=\underset{\alpha ⪰q_j}{min}{\theta }_{j,\alpha }\ge {\theta }_{j,{\alpha }^{\prime }}\ge {\theta }_{j,0},\text{for}\text{all}{\alpha }^{\prime }⪰̸q_j,\\ & \text{and}{\nu }_0=0.\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned}&\min _{{\varvec{\beta }}} ({\varvec{\beta }}- {\tilde{{\varvec{\beta }}}}^{(t)})^\top {\mathbf {D}}^{(t)} ({\varvec{\beta }}- {\tilde{{\varvec{\beta }}}}^{(t)}), \\ s.t.~&\max _{\varvec{\alpha }\succeq {\mathbf {q}}_{j}} \theta _{j, \varvec{\alpha }} = \min _{\varvec{\alpha }\succeq {\mathbf {q}}_{j}} \theta _{j, \varvec{\alpha }} \ge \theta _{j, \varvec{\alpha }'}\ge \theta _{j,{\varvec{0}}}, ~\text{ for } \text{ all }~ \varvec{\alpha }'\nsucceq {\varvec{q}}_j,\\&~\text{ and }~ \nu _{{\mathbf {0}}} = 0. \end{aligned} \end{aligned}$$\end{document}

Quadratic programming is the most studied nonlinear convex optimization problem (Chapter 4, Boyd et al., 2004) and many efficient solvers exist. In our simulation study in Sect. 4.3, we use the dual method of Goldfarb & Idnani (1983) implemented in the R package quadprog (Turlach et al., 2019).

2.4. Comparison with Related Algorithms

We compare Algorithm 1 with several related algorithms in more details.

2.4.1. Robbins-Monro SA and Variants

The proposed method is closely related to the stochastic approximation approach first proposed in Robbins & Monro (1951), and its variants given in Gu & Kong (1998) and Cai (2010a) that are specially designed for latent variable model estimation. Note that the Robbins–Monro method is the first SGD method with convergence guarantee. Both the methods of Gu & Kong (1998) and Cai (2010a) approximate the original Robbins–Monro method by using MCMC sampling to generate an approximate stochastic gradient in each iteration, when an unbiased stochastic gradient is difficult to obtain. All these methods do not handle complex constraints or non-smooth objective functions.

When there is no constraint or penalty on parameters (i.e., $g\left(\beta \right)\equiv 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }}) \equiv 0$$\end{document} ), the proximal operator degenerates to an identity map. In this case, the proposed method is essentially the same as Gu & Kong (1998) and Cai (2010a), except for the sampling method in the stochastic step, the way the Hessian matrix is approximated, the specific choices of step size, and the averaging in the last step of the proposed method. Among these differences, the step size and the trajectory averaging are key to the advantage of the proposed method.

As pointed out in Remark 4, the Robbins–Monro procedure has the same general requirement on the step size as the proposed method. Specially, the Robbins–Monro procedure, as well as its MCMC variants (Gu & Kong 1998; Cai 2010a), typically let the step size ${\gamma }_t$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t$$\end{document} decay in the order 1/t as suggested by asymptotic theory (Fabian, 1968). However, this step is often too short at the early stage of the algorithm, resulting in poor performance in practice (Sect. 4.5.3., Spall, 2003). On the other hand, the proposed method adopts a longer step size. By further adopting Polyak–Ruppert averaging (Ruppert, 1988; Polyak & Juditsky, 1992), we show in Sect. 3 that the proposed method almost achieves the optimal convergence speed.

2.4.2. Perturbed Proximal Gradient Algorithm

Proximal gradient descent algorithm (Parikh & Boyd, 2014) is a non-stochastic algorithm for solving nonsmooth and/or constrained optimization algorithms. For example, the widely used gradient projection algorithm for oblique rotation in factor analysis (Jennrich, 2002) is a special case. The vanilla proximal gradient descent algorithm does not use the second-order information of the objective function and thus sometimes converges slowly. To improve convergence speed, proximal Newton-type methods have been proposed in Lee et al. (2014) that utilize the second-order information of the smooth part of the objective function.

The perturbed proximal gradient algorithm (Atchadé et al., 2017) solves a similar optimization problem as in (2) by combining the methods of stochastic approximation, proximal gradient decent, and Polyak–Ruppert averaging. The proposed method extends Atchadé et al., (2017) by adopting a Newton-type proximal update suggested in (2014). The method of Atchadé et al., (2017) can be viewed as a special case of the proposed one with $c_1=c_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_{1} = c_{2}$$\end{document} . As shown by simulation study in the sequel, thanks to the second-order information, the proposed method converges much faster than that of Atchadé et al., (2017). We also point out that the theoretical analysis of Atchadé et al., (2017) focuses on convex optimization, while in Sect. 3 we consider a more general setting of non-convex optimization that includes a wide range of latent variable model estimation problems as special cases.

2.4.3. Stochastic EM Algorithm

The proposed method is also closely related to the stochastic-EM algorithm (Celeux, 1985; Ip, 2002; Nielsen, 2000; S. Zhang et al., 2020b). The stochastic-EM algorithm is a similar iterative algorithm, consisting of a stochastic step and a maximization step in each iteration, where the stochastic step is the same as that in the proposed algorithm. The maximization step plays a similar role as the proximal step in the proposed algorithm. More precisely, when there is no constraint or penalty, the maximization step of the stochastic-EM algorithm obtains parameter update ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} by minimizing the negative complete data log-likelihood function $-logf(y,{\xi }^{\left(t\right)}\mid \beta )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$-\log f({\mathbf {y}}, \varvec{\xi }^{(t)}\mid {\varvec{\beta }})$$\end{document} , instead of a stochastic gradient update. It is also recommended to perform a trajectory averaging in the stochastic-EM algorithm (Nielsen, 2000; S. Zhang et al., 2020b), like the last step of the proposed algorithm. As pointed out in S. Zhang et al. (2020b), the stochastic EM algorithm can potentially handle constraints and non-smooth penalties on parameters by incorporating them into the maximization step.

The stochastic-EM algorithm is typically not as fast as the proposed method, which is revealed by simulation studies below. This is because, it requires to solve an optimization problem completely in each iteration, which is time consuming, especially when constraints and non-smooth penalties are involved. On the other hand, the proximal step of the proposed algorithm can often be efficiently performed.

3. Theoretical Properties

In what follows, we establish the asymptotic properties of the proposed algorithm, under suitable technical conditions. For readers who are not interested in the asymptotic theory, this section can be skipped without affecting the reading of the rest of the paper. Note that in this section, we view data as fixed and the randomness comes from sampling of the latent variables. The following expectation is taken with respect to latent variable $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} given data $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{y}}$$\end{document} and parameters $\beta ,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }},$$\end{document} denoted by $E(·\mid \beta )=\int ·{\pi }_{\beta }\left(\xi \right)d\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {E}}(\cdot \mid {\varvec{\beta }})=\int \cdot \pi _{{\varvec{\beta }}}(\varvec{\xi })d\varvec{\xi }$$\end{document} , where ${\pi }_{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{\varvec{\beta }}$$\end{document} is the posterior distribution for $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} given $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{y}}$$\end{document} and $\beta .$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}.$$\end{document} Let $\Vert ·\Vert$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert \cdot \Vert $$\end{document} denote the vector $l_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$l_2$$\end{document} -norm. Following the typical convergence analysis of non-convex optimization (e.g., Chapter 3, Floudas, 1995), we will first discuss the convergence of the sequence ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} to a stationary point of the objective function $h\left(\beta \right)+g\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }}) + g({\varvec{\beta }})$$\end{document} in Theorem 1, which follows the theoretical development in Duchi & Ruan (2018). Then with some additional assumptions on the local geometry of the objective function at the stationary point being converged to, we will show the convergence rate of the Polyak–Ruppert averaged sequence ${\overline{\beta }}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n$$\end{document} in Theorem 2 which extends the results of Atchadé et al. (2017) to the setting of non-convex optimization.

For a function $f:R^d\mapsto R\cup +\infty ,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f:{\mathbb {R}}^d\mapsto {\mathbb {R}}\cup +\infty ,$$\end{document} denote the Fréchet subdifferential (Chapter 8.B Rockafellar & Wets, 1998) of f at the point $x$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{x}}$$\end{document} by $\partial f\left(x\right),$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\partial f({\varvec{x}}),$$\end{document}

$\begin{array}{c}\partial f\left(x\right)=\left(z,\in ,R^p,:,f,\left(y\right),\ge ,f,\left(x\right),+,z^{\top },(y-x),+,o(\Vert y,-,x\Vert )\text{as}y,\to ,x\right).\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \partial f({\varvec{x}})=\left\{ {\varvec{z}} \in {\mathbb {R}}^{p}: f({\varvec{y}}) \ge f({\varvec{x}})+ {\varvec{z}}^\top ({\varvec{y}}-{\varvec{x}})+o(\Vert {\varvec{y}}-{\varvec{x}}\Vert ) \text{ as } {\varvec{y}} \rightarrow {\varvec{x}}\right\} . \end{aligned}$$\end{document}

Define the set of stationary points of the objective function as

$\begin{array}{c}{B}^{\ast }=\{\beta \in B:\exists x\in \partial h\left(\beta \right)+\partial g\left(\beta \right)\text{with}{x}^{\top }(y-\beta )\ge 0,\text{for}\text{all}y\in B\}.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathcal {B}}^* = \{{\varvec{\beta }}\in {\mathcal {B}}: \exists \ {\mathbf {x}} \in \partial h({\varvec{\beta }}) + \partial g({\varvec{\beta }}) ~\text{ with }~ {\mathbf {x}}^\top ({\mathbf {y}}-{\varvec{\beta }}) \ge 0, ~\text{ for } \text{ all }~ {\mathbf {y}}\in {\mathcal {B}}\}. \end{aligned}$$\end{document}

Note that the global minimum $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} is a stationary point, i.e., $\hat{\beta }\in B^{\ast }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}} \in {\mathcal {B}}^*$$\end{document} . In addition, when the objective function is smooth, i.e., $g\left(\beta \right)\equiv 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g({\varvec{\beta }}) \equiv 0$$\end{document} , then $B^{\ast }=\{\beta \in B:\nabla h\left(\beta \right)=0\},$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}^* = \{{\varvec{\beta }}\in {\mathcal {B}}: \nabla h({\varvec{\beta }}) = 0 \},$$\end{document} which is the standard definition of stationary points set for a smooth function.

The following assumptions are assumed for our objective function.

1. H1. $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}$$\end{document} is compact and contains finite stationary points. For stationary points $\beta ,{\beta }^{\prime }\in B^{\ast },$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }},{\varvec{\beta }}^\prime \in {\mathcal {B}}^*,$$\end{document} $h\left(\beta \right)+g\left(\beta \right)=h\left({\beta }^{\prime }\right)+g\left({\beta }^{\prime }\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})+ g({\varvec{\beta }})= h({\varvec{\beta }}^\prime )+g({\varvec{\beta }}^\prime )$$\end{document} if and only if $\beta ={\beta }^{\prime }.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}= {\varvec{\beta }}^\prime .$$\end{document}

2. H2. $H(\xi ,\beta )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$H(\varvec{\xi },{\varvec{\beta }})$$\end{document} is a differentiable function with respect to $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} for given $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} and let $G_{\beta }\left(\xi \right)=\partial H(\xi ,\beta )/\partial \beta .$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}_{{\varvec{\beta }}}(\varvec{\xi })= \partial H(\varvec{\xi },{\varvec{\beta }})/\partial {\varvec{\beta }}.$$\end{document} Define function $M_{ϵ}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$M_\epsilon $$\end{document} : $\Theta \times \Xi \to R_{+}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Theta \times \Xi \rightarrow {\mathbb {R}}_+$$\end{document} as

   $\begin{array}{c}{M}_{ϵ}(\beta ;\xi )=\underset{{\beta }^{\prime }\in B,\Vert {\beta }^{\prime }-\beta \Vert <ϵ}{sup}\Vert G_{{\beta }^{\prime }}\left(\xi \right)\Vert .\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} M_\epsilon ({\varvec{\beta }};\varvec{\xi }) = \sup _{{\varvec{\beta }}'\in {\mathcal {B}},\Vert {\varvec{\beta }}'-{\varvec{\beta }}\Vert <\epsilon }\Vert {\mathbf {G}}_{{\varvec{\beta }}'}(\varvec{\xi })\Vert . \end{aligned}$$\end{document}
   There exists ${ϵ}_0>0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\epsilon _0>0$$\end{document} such that for all $0<ϵ<{ϵ}_0,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$0<\epsilon <\epsilon _0,$$\end{document}
   $\begin{array}{c}E[M_{ϵ}{(\beta ;\xi )}^2\mid \beta ]<\infty \text{for all}\beta \in B.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbb {E}}[M_\epsilon ({\varvec{\beta }};\varvec{\xi })^2\mid {\varvec{\beta }}]<\infty \text { for all }{\varvec{\beta }}\in {\mathcal {B}}. \end{aligned}$$\end{document}

3. H3. There exists ${ϵ}_0>0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\epsilon _0>0$$\end{document} such that for all ${\beta }^{\prime }\in B,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^\prime \in {\mathcal {B}},$$\end{document} there exists $\lambda (\xi ,{\beta }^{\prime })\ge 0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda (\varvec{\xi },{\varvec{\beta }}^\prime )\ge 0$$\end{document} such that

   $\begin{array}{c}\beta \mapsto H(\xi ,\beta )+\frac{\lambda (\xi ,{\beta }^{\prime })}{2}{\Vert \beta -{\beta }_0\Vert }^2\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\varvec{\beta }}\mapsto H(\varvec{\xi },{\varvec{\beta }}) + \frac{\lambda (\varvec{\xi },{\varvec{\beta }}^\prime )}{2}\Vert {\varvec{\beta }}-{\varvec{\beta }}_0\Vert ^2 \end{aligned}$$\end{document}
   is convex on the set $\{\beta :\Vert \beta -{\beta }^{\prime }\Vert \le {ϵ}_0\}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{\beta }}: \Vert {\varvec{\beta }}-{\varvec{\beta }}^\prime \Vert \le \epsilon _0\}$$\end{document} for any ${\beta }_0,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}_0,$$\end{document} and $E[\lambda (\xi ,{\beta }^{\prime })\mid \beta ]<\infty .$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {E}}[\lambda (\varvec{\xi },{\varvec{\beta }}^\prime )\mid {\varvec{\beta }}]<\infty .$$\end{document}

4. H4. The stochastic gradient $G_{{\beta }^{(t-1)}}\left({\xi }^{\left(t\right)}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}_{{\varvec{\beta }}^{(t-1)}}(\varvec{\xi }^{(t)})$$\end{document} is a Monte Carlo approximation of $\nabla h\left({\beta }^{(t-1)}\right).$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nabla h({\varvec{\beta }}^{(t-1)}).$$\end{document} That is, if computationally feasible, we take ${\xi }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }^{(t)}$$\end{document} as an exact sample from ${\pi }_{{\beta }^{(t-1)}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{{\varvec{\beta }}^{(t-1)}}$$\end{document} , where, as defined earlier, ${\pi }_{{\beta }^{(t-1)}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{{\varvec{\beta }}^{(t-1)}}$$\end{document} is the posterior distribution of $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} given $y$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{y}}$$\end{document} and ${\beta }^{(t-1)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t-1)}$$\end{document} . If not, we sample ${\xi }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }^{(t)}$$\end{document} from a Markov kernel $P_{{\beta }^{(t-1)}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$P_{{\varvec{\beta }}^{(t-1)}}$$\end{document} with invariant distribution ${\pi }_{{\beta }^{(t-1)}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{{\varvec{\beta }}^{(t-1)}}$$\end{document} .

5. H5. Define

   (20) $\begin{array}{c}\begin{array}{cc}& {\beta }_{\gamma }^{+}\left(\xi \right)=\underset{x\in B}{\text{argmin}}\left({\left[G_{\beta }\left(\xi \right)\right]}^{\top },(x-\beta ),+,g,\left(x\right),+,\frac{1}{2\gamma },{\left(x,-,\beta \right)}_D^2\right),\\ & U_{\gamma }(\xi ;\beta )=\frac{1}{\gamma }(\beta -{\beta }_{\gamma }^{+}\left(\xi \right)),\end{array}\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \begin{aligned}&{\varvec{\beta }}_{\gamma }^+(\varvec{\xi }) = \underset{{\varvec{x}} \in {\mathcal {B}}}{{\text {argmin}}}\left\{ [{\varvec{G}}_{{\varvec{\beta }}}(\varvec{\xi })]^\top ({\varvec{x}}-{\varvec{\beta }})+{\mathbf {g}}({\varvec{x}})+\frac{1}{2 \gamma }\left\| {\varvec{x}}-{\varvec{\beta }}\right\| ^{2}_{{\mathbf {D}}}\right\} ,\\&{\varvec{U}}_{\gamma }(\varvec{\xi };{\varvec{\beta }}) = \frac{1}{\gamma }({\varvec{\beta }}- {\varvec{\beta }}_{\gamma }^+(\varvec{\xi })), \end{aligned} \end{aligned}$$\end{document}
   where step size satisfy ${\sum }_{t=1}^{\infty }{\gamma }_t=\infty ,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{t=1}^\infty \gamma _t = \infty ,$$\end{document} ${\sum }_{t=1}^{\infty }{\gamma }_t^2<\infty .$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{t=1}^\infty \gamma _t^2<\infty .$$\end{document} Then with probability 1,
   $\begin{array}{c}\underset{n\to \infty }{lim}\sum _{t=1}^n{\gamma }_t\left(U_{{\gamma }_t},({\xi }^{\left(t\right)};{\beta }^{(t-1)}),-,E,[U_{{\gamma }_t}({\xi }^{\left(t\right)};{\beta }^{(t-1)})\mid {\beta }^{(t-1)}]\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \lim _{n\rightarrow \infty }\sum _{t=1}^n\gamma _t\left( {\varvec{U}}_{\gamma _t} (\varvec{\xi }^{(t)};{\varvec{\beta }}^{(t-1)}) - {\mathbb {E}}[{\varvec{U}}_{\gamma _t}(\varvec{\xi }^{(t)};{\varvec{\beta }}^{(t-1)})\mid {\varvec{\beta }}^{(t-1)}]\right) \end{aligned}$$\end{document}
   exists and is finite.

We remark that conditions H1 through H5 are quite mild. Condition H1 imposes mild requirements on the compactness of the parameter space and the properties of the stationary points of the objective function. Specifically, the compactness of the parameter space is often assumed when analyzing stochastic optimization problems without assuming convexity, see e.g., Gu & Kong (1998), Nielsen (2000), Cai (2010b), and Duchi & Ruan (2018). It also requires that the objective function has different values at different stationary points. Conditions H2 and H3 require the complete-data log-likelihood function $H(\xi ,·)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$H(\varvec{\xi },\cdot )$$\end{document} is locally Lipschitzian and weakly convex, respectively. These conditions hold when the complete-data log-likelihood function $H(\xi ,·)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$H(\varvec{\xi },\cdot )$$\end{document} is Lipschitzian and convex on the entire parameter space. Requiring locally Lipschitzian and weakly convex enables our theory to be applicable to a wider range of problems. Similar conditions are imposed in Duchi & Ruan (2018). For the examples that we consider in Sects. 1.2 and 1.3, these two conditions are satisfied because $H(\xi ,\beta )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$H(\varvec{\xi },{\varvec{\beta }})$$\end{document} is smooth and convex in $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} . Condition H4 is automatically satisfied according to the way the latent variables are sampled in Algorithm 1. Finally, H5 is a key condition for the convergence of the sequence ${\beta }^{\left(t\right)}.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}.$$\end{document} When exact samples from the posterior distribution are used, Lemma 1 guarantees that H5 is satisfied. With approximate samples from an MCMC algorithm, H5 may still hold when the bias from the MCMC samples is small.

Lemma 1

Define the filtration of $\sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma $$\end{document} -algebra $F_{t-1}=\sigma \left({\beta }^{\left(0\right)},,,{\xi }^{\left(k\right)},,,0,\le ,k,\le ,t,-,1\right).$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {F}}_{t-1} = \sigma \left( {\varvec{\beta }}^{(0)}, \varvec{\xi }^{(k)}, 0 \le k \le t-1\right) .$$\end{document} $\xi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }$$\end{document} is a sample from ${\pi }_{\beta }.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\pi _{\varvec{\beta }}.$$\end{document} Let

$\begin{array}{c}{ϵ}_{\gamma }(\xi ;\beta )=U_{\gamma }(\xi ;\beta )-E[U_{\gamma }(\xi ;\beta )\mid \beta ],\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \varvec{\epsilon }_{\gamma }(\varvec{\xi };{\varvec{\beta }}) = {\varvec{U}}_{\gamma }(\varvec{\xi };{\varvec{\beta }}) - {\mathbb {E}}[{\varvec{U}}_{\gamma }(\varvec{\xi };{\varvec{\beta }})\mid {\varvec{\beta }}], \end{aligned}$$\end{document}

then ${\gamma }_t{ϵ}_{{\gamma }_t}({\xi }^{\left(t\right)},{\beta }^{(t-1)})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t\varvec{\epsilon }_{\gamma _t}(\varvec{\xi }^{(t)},{\varvec{\beta }}^{(t-1)})$$\end{document} is a square-integrable martingale difference sequence adapted to $F_{t-1},$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {F}}_{t-1},$$\end{document} and with probability 1, ${lim}_n{\sum }_{t=1}^n{\gamma }_t{ϵ}_{{\gamma }_t}({\xi }^{\left(t\right)},{\beta }^{(t-1)})$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lim _n \sum _{t=1}^n \gamma _t \varvec{\epsilon }_{\gamma _t}(\varvec{\xi }^{(t)},{\varvec{\beta }}^{(t-1)})$$\end{document} exists and is finite.

Theorem 1

Apply Algorithm 1 to optimization problem (11) with step size ${\gamma }_t=t^{-\frac{1}{2}-ϵ},ϵ\in (0,\frac{1}{2}]$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t=t^{-\frac{1}{2}-\epsilon },\epsilon \in (0,\frac{1}{2}]$$\end{document} , for which conditions H1-H5 hold. Then with probability 1, the sequence ${\beta }^{\left(n\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(n)}$$\end{document} converges to a stationary point in $B^{\ast }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}^*$$\end{document} .

We remark that the convergence of the proposed method is similar to that of the EM algorithm. In fact, for marginal maximum likelihood estimation that is non-convex, the EM algorithm also only guarantees the convergence to a stationary point (Wu, 1983). Moreover, when the objective function has a single stationary point (e.g., when the objective function is strictly convex), then Theorem 1 guarantees global convergence.

The convergence of ${\beta }^{\left(n\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(n)}$$\end{document} guarantees the convergence of the Polyak–Ruppert averaging sequence ${\overline{\beta }}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n$$\end{document} . However, Theorem 1 does not provide information on the convergence speed. In what follows, we establish the convergence speed of ${\overline{\beta }}_n$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_n$$\end{document} . Without loss of generality, by Theorem 1, we assume that ${\beta }^{\left(n\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(n)}$$\end{document} converges to ${\beta }_{\ast }\in B^{\ast }.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}_* \in {\mathcal {B}}^*.$$\end{document}

1. H6. There exists $\delta >0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta >0$$\end{document} , such that $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} is strongly convex in $B_1=\{\beta \in B:\Vert \beta -{\beta }_{\ast }\Vert \le \delta \}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}_1 = \{{\varvec{\beta }}\in {\mathcal {B}}: \Vert {\varvec{\beta }}- {\varvec{\beta }}_*\Vert \le \delta \}$$\end{document} and $\nabla h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nabla h({\varvec{\beta }})$$\end{document} is Lipschitz in $B_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}_1$$\end{document} with Lipschitz constant L.

2. H7. For $\beta ,{\beta }^{\prime }\in B_1,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }},{\varvec{\beta }}' \in {\mathcal {B}}_1,$$\end{document} any $\gamma >0,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma >0,$$\end{document} and diagonal matrix $D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {D}}$$\end{document} with diagonal entries ${\delta }_i\in [c_1,c_2],$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\delta _i\in [c_1,c_2],$$\end{document} the following conditions hold.

   1. (i) $g\left({\text{Prox}}_{\gamma ,g}^D,\left(\beta \right)\right)-g\left({\beta }^{\prime }\right)\le -\frac{1}{\gamma }{\left({\text{Prox}}_{\gamma ,g}^D,\left(\beta \right),-,{\beta }^{\prime },,,{\text{Prox}}_{\gamma ,g}^D,\left(\beta \right),-,\beta \right)}_D$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g\left( {\text {Prox}}^{{\mathbf {D}}}_{\gamma , g}({\varvec{\beta }})\right) -g\left( {\varvec{\beta }}^{\prime }\right) \le -\frac{1}{\gamma }\left\langle {\text {Prox}}^{{\mathbf {D}}}_{\gamma , g}({\varvec{\beta }})-{\varvec{\beta }}^{\prime }, {\text {Prox}}^{{\mathbf {D}}}_{\gamma , g}({\varvec{\beta }})-{\varvec{\beta }}\right\rangle _{{\mathbf {D}}}$$\end{document} .

   2. (ii) ${\left({\text{Prox}}_{\gamma ,g}^D,\left(\beta \right),-,{\text{Prox}}_{\gamma ,g}^D,\left({\beta }^{\prime }\right)\right)}_D^2+{\left(\left({\text{Prox}}_{\gamma ,g}^D,\left(\beta \right),-,\beta \right),,-,,\left({\text{Prox}}_{\gamma ,g}^D,\left({\beta }^{\prime }\right),-,{\beta }^{\prime }\right)\right)}_D^2\le {\left(\beta ,-,{\beta }^{\prime }\right)}_D^2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\left\| {\text {Prox}}^{{\mathbf {D}}}_{\gamma , g}({\varvec{\beta }})-{\text {Prox}}^{{\mathbf {D}}}_{\gamma ,{} g}\left( {\varvec{\beta }}^{\prime }\right) \right\| _{{\mathbf {D}}}^{2} +\left\| \left( {\text {Prox}}^{{\mathbf {D}}}_{\gamma , g}({\varvec{\beta }})-{\varvec{\beta }}\right) \!-\!\left( {\text {Prox}}^{{\mathbf {D}}}_{\gamma , g}\left( {\varvec{\beta }}^{\prime }\right) -{\varvec{\beta }}^{\prime }\right) \right\| _{{\mathbf {D}}}^{2} \le \left\| {\varvec{\beta }}-{\varvec{\beta }}^{\prime }\right\| _{{\mathbf {D}}}^{2}$$\end{document} .

   3. (iii) ${sup}_{\gamma \in (0,c_1/L]}{sup}_{\beta \in B_1}{\gamma }^{-1}\left({\text{Prox}}_{\gamma ,g}^D,\left(\beta \right),-,\beta \right)<\infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sup _{\gamma \in (0,c_1 / L]} \sup _{{\varvec{\beta }}\in {\mathcal {B}}_1} \gamma ^{-1}\left\| {\text {Prox}}_{\gamma , g}^{{\mathbf {D}}}({\varvec{\beta }})-{\varvec{\beta }}\right\| <\infty $$\end{document} .

3. H8. For a measurable function $V:\Xi \to [1,+\infty ),$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V: \Xi \rightarrow [1,+\infty ),$$\end{document} a signed measure $\mu$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mu $$\end{document} on the $\sigma$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma $$\end{document} -field of $\Xi ,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Xi ,$$\end{document} and a function $f:\Xi \to R,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f: \Xi \rightarrow {\mathbb {R}},$$\end{document} define

   $\begin{array}{c}{\left|f\right|}_V\stackrel{\text{def}}{=}\underset{\xi \in \Xi }{sup}\frac{\left|f\right(\xi \left)\right|}{V\left(\xi \right)},{\Vert \mu \Vert }_V\stackrel{\text{def}}{=}\underset{{f,\left|f\right|}_V\le 1}{sup}\left(\int ,f,d,\mu \right).\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} |f|_{V} {\mathop {=}\limits ^{ \text{ def } }} \sup _{\varvec{\xi } \in \Xi } \frac{|f(\varvec{\xi })|}{V(\varvec{\xi })}, \quad \Vert \mu \Vert _{V} {\mathop {=}\limits ^{ \text{ def } }} \sup _{f,|f|_{V} \le 1}\left| \int f \mathrm {d} \mu \right| . \end{aligned}$$\end{document}
   There exist $\lambda \in (0,1),b<\infty ,m\ge 4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda \in (0,1), b<\infty , m\ge 4$$\end{document} and a measurable function $W:\Xi \to [1,+\infty )$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$W: \Xi \rightarrow [1,+\infty )$$\end{document} such that
   $\begin{array}{c}\underset{\beta \in B_1}{sup}{\left(G_{\beta }\right)}_W<\infty ,\underset{\beta \in B_1}{sup}{P}_{\beta }{W}^m\le \lambda W^m+b,\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \sup _{{\varvec{\beta }}\in {\mathcal {B}}_1}\left| {\mathbf {G}}_{{\varvec{\beta }}}\right| _{W}<\infty , \quad \sup _{{\varvec{\beta }}\in {\mathcal {B}}_1} P_{{\varvec{\beta }}} W^{m} \le \lambda W^{m}+b, \end{aligned}$$\end{document}
   where $G_{\beta }\left(\xi \right)=\partial H(\xi ,\beta )/\partial \beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {G}}_{{\varvec{\beta }}}(\varvec{\xi })= \partial H(\varvec{\xi },{\varvec{\beta }})/\partial {\varvec{\beta }}$$\end{document} and $P_{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$P_{{\varvec{\beta }}}$$\end{document} is the Markov kernel defined in condition H4. In addition, for any $\ell \in (0,m],$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\ell \in (0,m], $$\end{document} there exists $C<\infty ,\rho \in (0,1)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C<\infty , \rho \in (0,1)$$\end{document} such that for any $\xi \in \Xi ,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }\in \Xi ,$$\end{document}
   $\begin{array}{c}\underset{\beta \in B_1}{sup}{\left(P_{\beta }^n,(\xi ,·),-,{\pi }_{\beta }\right)}_{W^{\ell }}\le C{\rho }^n{W}^{\ell }\left(\xi \right).\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \sup _{{\varvec{\beta }}\in {\mathcal {B}}_1}\left\| P_{{\varvec{\beta }}}^{n}(\varvec{\xi }, \cdot )-\pi _{{\varvec{\beta }}}\right\| _{W^{\ell }} \le C \rho ^{n} W^{\ell }(\varvec{\xi }). \end{aligned}$$\end{document}

4. H9. There exists a constant C such that for any $\beta ,{\beta }^{\prime }\in B_1,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }},{\varvec{\beta }}^{\prime } \in {\mathcal {B}}_1, $$\end{document}

   $\begin{array}{c}{\left(G_{\beta },-,G_{{\beta }^{\prime }}\right)}_W+\underset{\xi \in \Xi }{sup}\frac{{\left(P_{\beta },(\xi ,·),-,P_{{\beta }^{\prime }},(\xi ,·)\right)}_W}{W\left(\xi \right)}+{\left({\pi }_{\beta },-,{\pi }_{{\beta }^{\prime }}\right)}_W\le C\left(\beta ,-,{\beta }^{\prime }\right).\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \left| {\mathbf {G}}_{{\varvec{\beta }}}-{\mathbf {G}}_{{\varvec{\beta }}^{\prime }}\right| _{W}+\sup _{\varvec{\xi }\in \Xi } \frac{\left\| P_{{\varvec{\beta }}}(\varvec{\xi }, \cdot )-P_{{\varvec{\beta }}^{\prime }}(\varvec{\xi }, \cdot )\right\| _{W}}{W(\varvec{\xi })}+\left\| \pi _{{\varvec{\beta }}}-\pi _{{\varvec{\beta }}^{\prime }}\right\| _{W} \le C\left\| {\varvec{\beta }}-{\varvec{\beta }}^{\prime }\right\| . \end{aligned}$$\end{document}

We provide a few remarks on conditions H6-H9, which are needed for establishing the convergence speed in addition to conditions H1–H5. Condition H6 requires that the smooth part of the objective function is strongly convex and its derivative is Lipschitz continuous in a small neighborhood of ${\beta }_{\ast }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}_*$$\end{document} . Specifically, $h\left(\beta \right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$h({\varvec{\beta }})$$\end{document} being strongly convex in $B_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {B}}_1$$\end{document} means that there exists a positive constant C, such that ${(\nabla h\left(\beta \right)-\nabla h\left({\beta }^{\prime }\right))}^{\top }(\beta -{\beta }^{\prime })\ge C{\Vert \beta -{\beta }^{\prime }\Vert }^2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$(\nabla h({\varvec{\beta }})-\nabla h({\varvec{\beta }}'))^\top ({\varvec{\beta }}- {\varvec{\beta }}') \ge C\Vert {\varvec{\beta }}- {\varvec{\beta }}'\Vert ^2$$\end{document} , for any $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} and ${\beta }^{\prime }\in B_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}' \in {\mathcal {B}}_1$$\end{document} . Condition H7 imposes some requirements on the non-smooth part of the objective function, with regard to the proximal operator. As verified in Lemma C.1, H7 holds when g is a generalized function that indicates constraints or when g is locally Lipschitz continuous and convex that holds when g is a $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} regularization function. Thus, H7 holds for the examples we consider in Sects. 2.2 and 2.3. Conditions H8 and H9 impose mild regularity conditions on the stochastic gradient in a local neighborhood of ${\beta }_{\ast }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}_*$$\end{document} , especially when the stochastic gradients are generated by a Markov kernel. These conditions are used to control the bias caused by MCMC sampling. H8 is essentially a uniform-in- $\beta$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}$$\end{document} ergodic condition and H9 is a local Lipschitzian condition on the Markov kernel. These regularity conditions are commonly adopted in the stochastic approximation literature (Benveniste et al., 1990; Andrieu et al., 2005; Fort et al., 2016) and have been shown to hold for general families of MCMC kernels including Metropolis–Hastings and Gibbs samplers (Andrieu & Moulines, 2006; Fort et al., 2011; Schmidt et al., 2011).

Theorem 2

Suppose that H1–H9 hold. Then there exists a constant C, such that for the Polyak–Ruppert averaging sequence ${\overline{\beta }}_n=\frac{1}{n}{\sum }_{t=1}^n{\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\bar{{\varvec{\beta }}}}_{n} = \frac{1}{n}\sum _{t=1}^n {\varvec{\beta }}^{(t)}$$\end{document} from Algorithm 1,

(21) $\begin{array}{c}E\Vert {\overline{\beta }}_n-{\beta }_{\star }{\Vert }^2\le C{n}^{-\frac{1}{2}+\epsilon }.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\mathbb {E}}\Vert {\bar{{\varvec{\beta }}}}_n - {\varvec{\beta }}_\star \Vert ^2 \le C n^{-\frac{1}{2} + \varepsilon }. \end{aligned}$$\end{document}

Note that the expectation is taken with respect to ${\xi }^{\left(1\right)},\dots ,{\xi }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }^{(1)},\ldots ,\varvec{\xi }^{(t)}$$\end{document} given ${\beta }^{\left(0\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(0)}$$\end{document} and ${\xi }^{\left(0\right)}.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\xi }^{(0)}.$$\end{document}

We now provide a few remarks regarding the convergence speed (21). First, the small positive constant $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} comes from the requirement on step size that ${\sum }_{t=1}^{\infty }{\gamma }_t^2<\infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{t=1}^\infty \gamma _t^2 < \infty $$\end{document} in H5. Since ${\sum }_{t=1}^{\infty }{\gamma }_t^2<\infty$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sum _{t=1}^\infty \gamma _t^2 < \infty $$\end{document} is satisfied when ${\gamma }_t=\mu t^{-\frac{1}{2}-\epsilon }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t = \mu t^{-\frac{1}{2} - \varepsilon }$$\end{document} , for any $\mu ,\epsilon >0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mu , \varepsilon > 0$$\end{document} , the convergence speed of $E\Vert {\overline{\beta }}_n-{\beta }_{\star }{\Vert }^2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {E}}\Vert {\bar{{\varvec{\beta }}}}_n - {\varvec{\beta }}_\star \Vert ^2$$\end{document} can be arbitrarily close to $O\left(n^{-\frac{1}{2}}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$O(n^{-\frac{1}{2}})$$\end{document} by choosing an arbitrarily small $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} . Second, this $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} might be an artifact due to our proof strategy to overcome the non-convexity of the problem. In fact, if the objective function is convex, similar to Atchadé et al. (2017), we can choose $ϵ=0$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\epsilon = 0$$\end{document} and then prove under similar conditions that $E\Vert {\overline{\beta }}_n-{\beta }_{\star }{\Vert }^2\le C{n}^{-\frac{1}{2}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {E}}\Vert {\bar{{\varvec{\beta }}}}_n - {\varvec{\beta }}_\star \Vert ^2 \le C n^{-\frac{1}{2}}$$\end{document} . Lastly, it is well-known that for non-smooth convex optimization, the minimax optimal convergence rate is $O\left(n^{-\frac{1}{2}}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$O(n^{-\frac{1}{2}})$$\end{document} ; see Chapter 3, Nesterov (2004). In this sense, our algorithm is almost minimax optimal, when $\epsilon$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varepsilon $$\end{document} is very close to zero. It is well-known that Polyak–Ruppert averaging typically improves the convergence speed of a slowly convergent sequence (Ruppert, 1988, Bonnabel, 2013).

4. Simulation Study

4.1. Study I: Confirmatory IFA

Table 1. Comparison of five stochastic algorithms.

In the first study, we compare the performance of four variants of the proposed method and the stochastic EM (StEM) algorithm. The five methods, including their abbreviations, are given in Table 1. For a fair comparison, the same Gibbs sampling method is used. We further explain the differences below.

1. 1. USP is the method that we recommend. It has a step size ${\gamma }_t$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t$$\end{document} close to $t^{-1/2}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$t^{-1/2}$$\end{document} , applies Polyak–Ruppert averaging, and uses a quasi-Newton update in the proximal step.

2. 2. The USP-PPG method is the perturbed proximal gradient method that is implemented the same as the USP method except that $c_1=c_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_1 = c_2$$\end{document} so that it does not involve a quasi-Newton update. $c_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_1$$\end{document} is set to be 1 without tuning in this study.

3. 3. The USP-RM1 method is implemented the same as the USP method, except that ${\beta }^{\left(n\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(n)}$$\end{document} from the last iteration is taken as the estimator instead of applying Polyak–Ruppert averaging. This method is very similar to a Robbins–Monro algorithm, except for the update of parameters $B$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {B}}$$\end{document} for the covariance matrix where constraints involve.

4. 4. The USP-RM2 method is the same as USP-RM1, except that we set the step size ${\gamma }_t=1/t$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t = 1/t$$\end{document} which is the asymptotic optimal step size for the Robbins–Monro algorithm (Fabian, 1968).

5. 5. The implementation of the StEM algorithm is the same as USP, except for the proximal step. Instead of making stochastic gradient update, StEM obtains ${\beta }^{\left(t\right)}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{\beta }}^{(t)}$$\end{document} by completely solving an optimization problem

   $\begin{array}{c}{\beta }^{\left(t\right)}=\underset{\beta \in B}{\text{arg max}}H({\xi }^{\left(t\right)},\beta )+g\left(\beta \right).\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} {\varvec{\beta }}^{(t)} = \mathop {\hbox {arg max}}\limits _{{\varvec{\beta }}\in {\mathcal {B}}} H(\varvec{\xi }^{(t)}, {\varvec{\beta }}) + g({\varvec{\beta }}). \end{aligned}$$\end{document}
   In our implementation, this optimization problem is solved by making the quasi-Newton proximal update (14) iteratively until convergence.

We consider a confirmatory IFA setting with only two factors (i.e., $K=2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=2$$\end{document} ), so that an EM algorithm with sufficient numbers of quadrature points and EM steps can be used to obtain a more accurate approximation of $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} that will be used as the standard when comparing the five methods. We emphasize that it is important to compare the convergence speed of difference algorithms based on $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} rather than the true model parameters. This is because, under suitable conditions, these algorithms converge to $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} rather than the true model parameters. If we compare the algorithms based on the true model parameters, the difference in the convergence speed cannot be observed clearly, as the statistical error (i.e., the difference between $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} and the true model parameters) tends to dominate the computational errors (i.e., the difference between $\hat{\beta }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{{\varvec{\beta }}}}$$\end{document} and the results given by the stochastic algorithms).

More precisely, we consider sample size $N=1000$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N = 1000$$\end{document} and the number of items $J=20$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$J = 20$$\end{document} . The design matrix $Q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}}$$\end{document} is specified by the assumptions that items 1 through 5 only measure the first factor, items 6 through 10 only measure the second factor, and items 11 through 20 measure both. The intercept parameters $d_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j$$\end{document} are drawn i.i.d. from the standard normal distribution, and the non-zero loading parameters are drawn i.i.d. from a uniform distribution over the interval (0.5, 1.5). The variances of the two factors are set to be 1 and the covariance is set to be 0.4. Under these parameters, 100 independent datasets are generated, based on which the five methods are compared. To ensure a fair comparison, the true parameters are used as the starting point for all the methods. In addition, 1000 iterations are run (i.e., $n=1000$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n = 1000$$\end{document} ) for each method, instead of using an adaptive stopping criterion. For USP, USP-PPG, and StEM, the burn-in size $\varpi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varpi $$\end{document} is chosen to be 500. All algorithms are implemented in C++ and run on the same platform¹ using a single core.

Figure 1. The boxplot of mean squared errors for estimated parameters from the five methods.

Figure 2. The boxplot of mean squared errors for estimated parameters from ‘USP,’ ‘USP-RM1,’ and ‘StEM’ method.

The results regarding the accuracy of the proposed methods are given in Figs. 1 and 2 that are based on the following performance metrics. Specifically, for the intercept parameters $d_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j$$\end{document} , the following mean squared error (MSE) is calculated for each simulated dataset and each method,

$\begin{array}{c}\frac{1}{J}\sum _{j=1}^J{\left({\tilde{d}}_j,-,{\hat{d}}_j\right)}^2,\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \frac{1}{J}\sum _{j=1}^{J} \left( {\tilde{d}}_j - {\hat{d}}_j\right) ^2, \end{aligned}$$\end{document}

where ${\hat{d}}_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{d}}_j$$\end{document} , which is treated as the global optimum, is obtained by an EM algorithm with 31 Gaussian–Hermite quadrature points per dimension, and ${\tilde{d}}_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{d}}_j$$\end{document} is given by one of the five stochastic methods after 1000 iterations. Similarly, the MSEs for the loading parameters and for the correlation ${\sigma }_{12}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma _{12}$$\end{document} between the factors are calculated, where the MSE for the loading parameters is calculated for the unrestricted ones, i.e.,

$\begin{array}{c}\frac{{\sum }_{j,k}{1}_{\{q_{\mathrm{jk}}\ne 0\}}{({\tilde{a}}_{\mathrm{jk}}-{\hat{a}}_{\mathrm{jk}})}^2}{{\sum }_{j,k}{1}_{\{q_{\mathrm{jk}}\ne 0\}}}.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \frac{\sum _{j,k} 1_{\{q_{jk}\ne 0\}} ({\tilde{a}}_{jk} - {\hat{a}}_{jk})^2 }{\sum _{j,k} 1_{\{q_{jk}\ne 0\}} }. \end{aligned}$$\end{document}

Again, ${\hat{a}}_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\hat{a}}_{jk}$$\end{document} is given by the EM algorithm, and ${\tilde{a}}_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\tilde{a}}_{jk}$$\end{document} is given by one of the five methods.

Figure 1 compares the accuracy of all the five methods. As we can see, the USP, USP-RM1, and StEM methods have much smaller MSEs than the USP-PPG and USP-USP-RM2 methods. Since the USP-PPG method only differs from the USP method by whether using a quasi-Newton update, the inferior performance of USP-PPG implies the importance of the second-order information in the stochastic proximal gradient update. As the USP-RM2 method only differs from USP-RM1 by their step sizes, the inferior performance of USP-RM2 is mainly due to the use of short step size.

In Fig. 2, we zoom in to further compare the USP, USP-RM1, and StEM methods. First, we see that the USP method performs the best among the three, for all the parameters. As the USP-RM1 method is the same as the USP method except for not applying Polyak–Ruppert averaging, this result suggests that averaging does improve accuracy. Moreover, the USP method and the StEM method only differ by the way the parameters are updated, where the USP method takes a quasi-Newton proximal update, while the StEM method completely solves an optimization problem. It is likely that the way parameters are updated in the USP method yields more smoothing (i.e., averaging) than the StEM, which leads to the outperformance of the USP method.

Table 2. The elapsed time (seconds) for the five methods in confirmatory IFA.

On the computational efficiency, we show in Table 2 the elapsed time for the five methods. ‘USP’, ‘USP-RM1’, ‘USP-PPG’, and ‘USP-RM2’ share similar computation time since their floating point operations per iteration are at the same level. ‘StEM’ is most time consuming because an inner loop of optimization is involved in each iteration. In summary, the proposed USP algorithm is computationally the most efficient among the five algorithms, in the sense that it achieves the highest accuracy (see Figs. 1 and 2), within a similar or smaller amount of time (see Table 2).

4.2. Study II: Exploratory IFA by Regularization

In the second study, we apply the proposed method to regularized estimation for exploratory IFA as discussed in Sect. 1.2. We consider increasing sample size $N=1000,2000,4000,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N=1000, 2000, 4000,$$\end{document} eighty items and five correlated latent factors (i.e., $J=80,K=5$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$J=80,K=5$$\end{document} ). The true loading matrix is sparse, where the items each factor loads on are given in Table 3. Similar to Study I, the intercept parameters $d_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$d_j$$\end{document} are drawn i.i.d. from the standard normal distribution, and the non-zero loading parameters $a_{\mathrm{jk}}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$a_{jk}$$\end{document} are drawn i.i.d. from a uniform distribution over the interval (0.5, 1.5). The elements of covariance matrix $\Sigma ={\left({\sigma }_{k,k}\right)}_{5\times 5}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Sigma =(\sigma _{k,k})_{5\times 5}$$\end{document} are set to be ${\sigma }_{k,k^{\prime }}=1,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma _{k,k'}=1,$$\end{document} for $k=k^{\prime }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k=k'$$\end{document} and ${\sigma }_{k,k^{\prime }}=0.4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\sigma _{k,k'}=0.4$$\end{document} for $k\ne k^{\prime }.$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k\ne k'.$$\end{document}

Table 3. The sparse loading structure in the data generation IFA model.

For each sample size, 50 independent datasets are generated. In the proposed algorithm, we adopt a burn-in size $\varpi =50$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varpi =50$$\end{document} and stop based on the criterion discussed in Sect. 2, where the stopping threshold is set to be ${10}^{-3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$10^{-3}$$\end{document} . A decreasing penalty parameter ${\lambda }_N=\sqrt{logJ/N}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\lambda _N = \sqrt{\log J /N}$$\end{document} is used to ensure estimation consistency (Chapter 6, Bühlmann & van de Geer, 2011). Other implementation details can be found in Sect. 2.2. The algorithm in this example is implemented in C++ and is run on the same platform as in Study I. Although a regularized EM algorithm (Sun et al., 2016) can also solve this problem, it suffers from a very high computational cost. Due to the five-dimensional numerical integrals involved, it takes a few hours to fit one dataset. We thus do not consider it here.

We focus on the accuracy in the estimation of the loading matrix $A={\left(a_{\mathrm{jk}}\right)}_{J\times K}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {A}} = (a_{jk})_{J\times K}$$\end{document} . Note that although the rotational indeterminacy issue is resolved in this regularized estimator, the loading matrix can still only be identified up to column swapping. That is, two estimates of the loading matrix have the same objective function value, if one can be obtained by swapping the columns of the other. The following mean-squared-error measure is used that takes into account column swapping of the loading matrix

(22) $\begin{array}{c}\underset{A^{\prime }\in P\left(\tilde{A}\right)}{min}\left(\frac{1}{\mathrm{JK}},{\Vert A^{\prime }-A\Vert }_F^2\right),\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \min _{{\mathbf {A}}'\in {\mathcal {P}}(\tilde{{\mathbf {A}}})}\left\{ \frac{1}{JK}\Vert {\mathbf {A}}' - {\mathbf {A}}\Vert _F^2\right\} , \end{aligned}$$\end{document}

where ${\Vert ·\Vert }_F$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\Vert \cdot \Vert _F$$\end{document} is the Frobenius norm, $A$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {A}}$$\end{document} is the true loading matrix, $\tilde{A}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tilde{{\mathbf {A}}}$$\end{document} is the output of Algorithm 1, and $P\left(\tilde{A}\right)$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {P}}(\tilde{{\mathbf {A}}})$$\end{document} denotes the set of $J\times K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$J\times K$$\end{document} matrices that can be obtained by swapping the columns of $\tilde{A}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tilde{{\mathbf {A}}}$$\end{document} .

Results are given in Tables 4 and 5. In Table 4, we see that the MSE for the loading matrix is quite small and decreases as the sample size grows, suggesting that consistency of the regularized estimator. In Table 5, the quantiles of time consumption under different sample sizes are given, which suggests the computational efficiency of the proposed method.

Table 4. The mean squared errors for estimated loading parameters in exploratory IFA with $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} regularization.

Table 5. The elapsed time (seconds) for exploratory IFA with $L_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$L_1$$\end{document} regularization.

4.3. Study III: Restricted LCA

In this study, we apply the proposed method to the estimation of a restricted latent class model as discussed in Sect. 1.3, where the optimization involves complex inequality constraints. Specifically, data are from a Deterministic Input, Noisy ‘And’ gate (DINA) model (Junker & Sijtsma, 2001) that is a special restricted latent class model. Note that the DINA assumptions are only used in the data generation. We solve optimization (10) which is based on a general restricted latent class model considered in Xu (2017) instead of the DINA model, mimicking the practical situation when the parametric form is unknown.

We consider a test consisting of twenty items (i.e., $J=20$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$J = 20$$\end{document} ) that measure four binary attributes (i.e., $K=4$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$K=4$$\end{document} ). Three sample sizes are considered, including $N=1000,2000$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N=1000, 2000$$\end{document} , and 4000. The design matrix $Q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}}$$\end{document} is given in Table 6. In addition, the guessing and slipping parameters $s_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$s_j$$\end{document} and $g_j$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$g_j$$\end{document} of the DINA model are drawn i.i.d. from a uniform distribution over the interval (0.05, 0.2), which gives the values of ${\theta }_{j,\alpha }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{j, \varvec{\alpha }}$$\end{document} . That is,

$\begin{array}{c}{\theta }_{j,\alpha }=\left(\begin{array}{cc}log((1-s_j)/s_j),& \text{if}\alpha ⪰q_j,\\ log(g_j/(1-g_j)),& \text{otherwise.}\end{array}\right)\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \theta _{j, \varvec{\alpha }} =\left\{ \begin{array}{cl} \log ( (1-s_j)/s_j), &{} ~\text{ if }~ \varvec{\alpha } \succeq {\mathbf {q}}_{j}, \\ \log ( g_j/(1-g_j)), &{} ~\text{ otherwise. } \end{array}\right. \end{aligned}$$\end{document}

Finally, we let ${\nu }_{\alpha }=0,$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu _{\varvec{\alpha }} = 0,$$\end{document} for all $\alpha \in {\{0,1\}}^K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\alpha } \in \{0, 1\}^K$$\end{document} , so that $P(\xi =\alpha )=1/2^K$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {P}}(\varvec{\xi } = \varvec{\alpha }) = 1/2^K$$\end{document} . According to the results in Xu (2017), the model parameters are identifiable, given the $Q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}}$$\end{document} -matrix in Table 6.

Table 6. The design matrix $Q$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {Q}}$$\end{document} for the restricted LCA model.

For each sample size, 50 independent datasets are generated. The proposed algorithm adopts a burn-in size $\varpi =50$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varpi =50$$\end{document} and stops based on the criterion discussed in Sect. 2, where the stopping threshold is set to be ${10}^{-3}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$10^{-3}$$\end{document} . Other implementation details can be found in Sect. 2.3. The following metrics are used to evaluate the estimation accuracy. For item parameters ${\theta }_{j,\alpha }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{j, \varvec{\alpha }}$$\end{document} , the MSE is calculated as

$\begin{array}{c}\frac{1}{J\times 2^K}\sum _{j=1}^J\sum _{\alpha \in {\{0,1\}}^K}{\left({\tilde{\theta }}_{j,\alpha },-,{\theta }_{j,\alpha }\right)}^2.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \frac{1}{J\times 2^K}\sum _{j=1}^J\sum _{\varvec{\alpha }\in \{0,1\}^K}\left( {\tilde{\theta }}_{j,\varvec{\alpha }}-\theta _{j,\varvec{\alpha }}\right) ^2. \end{aligned}$$\end{document}

For structural parameters ${\nu }_{\alpha }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu _{\varvec{\alpha }}$$\end{document} , the MSE is calculated as

$\begin{array}{c}\frac{1}{2^K-1}\sum _{\alpha \in {\{0,1\}}^K,\alpha \ne 0}{\left({\tilde{\nu }}_{\alpha },-,{\nu }_{\alpha }\right)}^2.\end{array}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\begin{aligned} \frac{1}{2^K-1}\sum _{\varvec{\alpha }\in \{0,1\}^K,\ \varvec{\alpha }\ne {\varvec{0}}}\left( {\tilde{\nu }}_{\varvec{\alpha }} - \nu _{\varvec{\alpha }}\right) ^2. \end{aligned}$$\end{document}

Our results are given in Tables 7 and 8. As we can see, the estimation becomes more accurate as the sample size increases for both sets of parameters. It confirms that the current model is identifiable as suggested by Xu (2017) and thus can be consistently estimated.

Table 7. The MSE for item parameters ${\theta }_{j,\alpha }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\theta _{j, \varvec{\alpha }}$$\end{document} in the restricted latent class model.

Table 8. The MSE for structural parameters ${\nu }_{\alpha }$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\nu _{\varvec{\alpha }}$$\end{document} in the restricted latent class model.

5. Concluding Remarks

In this paper, a unified stochastic proximal optimization framework is proposed for the computation of latent variable model estimation. This framework is very general that applies to a wide range of estimators for almost all commonly used latent variable models. Comparing with existing stochastic optimization methods, the proposed method not only solves a wider range of problems including regularized and constrained estimators, but also is computationally more efficient. Theoretical properties of the proposed method are established. These results suggest that the convergence speed of the proposed method is almost optimal in the minimax sense.

The power of the proposed method is shown via three examples, including confirmatory IFA, exploratory IFA by regularized estimation, and restricted latent class analysis. Specifically, the proposed method is compared with several stochastic optimization algorithms, including a stochastic-EM algorithm and a Robbin–Monro algorithm with MCMC sampling, in the simulation study of confirmatory IFA, where there is no complex constraint or penalty. Using the same starting point and the same number of iterations, the proposed one is always more accurate than its competitors. The simulation studies on exploratory IFA and restricted latent class analysis further show the power of the proposed method for handling optimization problems with non-smooth penalties and complex inequality constraints.

The implementation of the proposed algorithm involves several tuning parameters. First, we need to choose a step size ${\gamma }_t$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t$$\end{document} . Our theoretical results suggest that ${\gamma }_t=t^{-0.5-ϵ}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t =t^{-0.5 - \epsilon }$$\end{document} for any $ϵ\in (0,0.5]$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\epsilon \in (0,0.5]$$\end{document} , and a smaller $ϵ$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\epsilon $$\end{document} leads to faster convergence. In practice, we suggest to set ${\gamma }_t=t^{-0.51}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t =t^{-0.51}$$\end{document} that performs well in all our simulations. This choice of step size is very different from the choice of ${\gamma }_t=t^{-1}$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\gamma _t =t^{-1}$$\end{document} in the MCMC stochastic approximation algorithms. Second, a burn-in size $\varpi$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varpi $$\end{document} is needed. The burn-in in the proposed algorithm is similar to the burn-in in MCMC algorithms. It does not affect the asymptotic convergence of the algorithm but improves the finite sample performance. In practice, the burn-in size can be decided similarly as in MCMC algorithms by monitoring the parameter updates using trace plots. Third, two positive constraints $c_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_1$$\end{document} and $c_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_2$$\end{document} are needed to regularize the second-order matrix in the scaled proximal update. Depending on the scale of each particular problem, we suggest to choose $c_1$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_1$$\end{document} to be sufficiently small and $c_2$ \documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$c_2$$\end{document} to be sufficiently large. It is found that the performance of our algorithm is not sensitive to their choices. Finally, a stopping criterion is needed. We suggest to stop the iterative update by monitoring a window of successive differences in parameter updates.

The proposed framework may be improved from several aspects that are left for future investigation. First, the sampling strategy in the stochastic step needs further investigation. Although in theory any reasonable MCMC sampler can yield the convergence of the algorithm, a good sampler will lead to superior finite sample performance. More sophisticated MCMC algorithms need to be investigated regarding their performance under the proposed framework. Second, methods for parallel and distributed computing need to be developed. As we can see, many steps of Algorithm 1 can be performed independently. This enables us to design parallel and/or distributed computing systems for solving large-scale and/or distributed versions of latent variable model estimation problems (e.g., fitting models for assessment data from online learning platforms and large-scale mental health records). Finally, the performance of the proposed method under other latent variable models needs to be investigated. For example, the proposed method can also be applied to latent stochastic process models (e.g., Chow et al., 2016; Chen & Zhang, 2020) that are useful for analyzing intensive longitudinal data. These models bring additional challenges, as stochastic processes need to be sampled in the stochastic step of our algorithm.

In summary, the proposed method is computationally efficient, theoretically solid, and applicable to a broad range of latent variable model inference problems. Like the EM algorithm as the standard tool for low-dimensional latent variable models, we believe that the proposed method may potentially serve as the standard approach to the estimation of high-dimensional latent variable models.