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

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. Supplementary Information The online version contains supplementary material available at 10.1007/s11336-022-09863-9.


Introduction
Latent variable models have been playing a central role in psychometrics and related fields.
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 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 fullinformation 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 respective 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 simultane-ously 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; 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 The rest of the paper is organized as follows. In Section 2, we formulate latent variable model estimation as a general optimization problem which covers many commonly used estimators as special cases. In Section 3, a quasi-Newton stochastic proximal algorithm is proposed. Theoretical properties of the proposed algorithm are established in Section 4, 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 Section 5. We conclude with some discussions in Section 6. An R package has been developed and will be published online upon the acceptance of the current paper.

Problem Setup
We consider the estimation of a parametric latent variable model. We adopt a general setting, followed by concrete examples in Sections 2.2 and 2.3. Let Y be a random object representing observable data and let y be its realization. For example, in item factor analysis (IFA), Y represents (categorical) responses to all the items from all the respondents. A latent variable model specifies the distribution of Y by introducing a set of latent variables ξ P Ξ, where Ξ denotes the state space of the latent vector ξ. For example, in item factor analysis, ξ consists of the latent traits of all the respondents and Ξ is a Euclidian space. Let β " pβ 1 , ..., β p q J P B be a set of parameters in the model, where B denotes the parameter space. The goal is to estimate β given observed data y.
We consider an EB estimator which takes the form where f py, ξ | βq is a complete-data likelihood/pseudo-likelihood function that has an analytic form. We assume that the objective function lpβq is finite for any β P B and is also smooth in β.
The estimator is given by solving the following optimization problem where Rpβq 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.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.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.

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 ij P t0, 1u be a random variable denoting person i's response to item j and let y ij be its realization. Thus, we have Y " pY ij q NˆJ and y " py ij q NˆJ , where Y and y are the generic notations introduced in Section 2.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 ξ " pξ ik q NˆK , where ξ ik represents person i's kth factor. Recall that ξ is our generic notation for the latent variables in Section 2.1 and here the state space Ξ " R NˆK . Throughout this paper, we assume the number of factors K is known.
An IFA model makes the following assumptions: 1. ξ i " pξ i1 , ..., ξ iK q J , i " 1, ..., N , are independent and identically distributed (i.i.d.) random vectors, following a multivariate normal distribution N p0, Σq. The diagonal terms of Σ " pσ kk 1 q KˆK are set to one for model identification. As Σ is a positive semi-definite matrix, it is common to reparametrize Σ by Cholesky decomposition, where B " pb kk 1 q KˆK is a lower triangular matrix. Let b k be the kth row of B. Then }b k } " 1, k " 1, ..., K, since the diagonal terms of Σ are constrained to value 1.
2. Y ij given ξ i follows a Bernoulli distribution satisfying where d j and a j " pa j1 , ..., a jK q J are item-specific parameters. The parameters a jk are often known as the loading parameters.
3. Y i1 , ..., Y iJ are assumed to be conditionally independent given ξ i , which is known as the local independence assumption.
Note that we consider the most commonly used logistic model in (2.3). It is worth pointing out that the proposed algorithm also applies to the normal ogive (i.e. probit) model which assumes that PpY ij " 1 | ξ i q " Φpd j`a J j ξ i q. Under the current setting and using the reparametrization for Σ, our model parameters are β " tB, d j , a j , j " 1, ..., Ju.
The marginal likelihood function takes the form lpβq " where φpx | Bq is the density function for multivariate normal distribution N p0, BB J q. The K-dimensional integrals involved in (2.4) cause a high computational burden for a relatively large K (e.g., K ě 5).
IFA models are commonly used for both exploratory and confirmatory analyses. In exploratory IFA, an important problem is to learn a sparse loading matrix pa ij q JˆK from data, which facilitates the interpretation of the factors. One approach is by the L 1 -regularized estimator (Sun et al., 2016) which takes the form where the parameter space .., Ku, and the penalty term In Rpβq, λ ą 0 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 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 matrix has constraints, as reflected by the parameter space B. Here b kk 1 " 0 is due to that B is a lower triangle matrix and ř K k 1 "1 b 2 kk 1 " 1 is due to that the diagonal terms of Σ " BB J are all 1. We remark that it is possible to replace the L 1 penalty in Rpβq by other penalty functions for imposing sparsity, such as the elastic net penalty (Zou & Hastie, 2005) Rpβq " λ 1 where λ 1 , λ 2 ą 0 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 " pq jk q JˆK . If q jk " 0, then item j does not load on factor k and a jk is set to 0. Otherwise, a jk is freely estimated. These constraints lead to parameter space B " tβ : b kk 1 " 0, 1 ď k ă k 1 ď K; ř K k 1 "1 b 2 kk 1 " 1, a jk " 0 for q jk " 0, j " 1, ..., J, k " 1, ..., Ku. The MML estimator for confirmatory IFA is then given bŷ β " arg min βPB´l pβq. (2.8) 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β. The estimation of the asymptotic variance often requires to compute the Hessian matrix of lpβq atβ, which also involves intractable K-dimensional integrals. As we will see in Section 3.1, this Hessian matrix, as well as quantities taking a similar form, can be easily obtained as a by-product of the proposed algorithm.

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 " pY ij q NˆJ and y " py ij q NˆJ . In this context, Y ij " 1 means that item j is answered correctly and Y ij " 0 means an incorrect answer.
The restricted latent class model can be parameterized as follows.
3. Local independence is still assumed. That is, Y i1 , ..., Y iJ are conditionally independent given ξ i .
We consider a confirmatory setting where there exists a design matrix, similar to the Q-matrix in confirmatory IFA. With slight abuse of notation, we still denote Q " pq jk q JˆK , where q jk P t0, 1u. Here, q jk " 1 indicates that solving item j requires the kth skill and q jk " 0 otherwise. As will be explained below, this design matrix leads to equality and inequality constraints in model parameters.
Denote q j " pq j1 , ..., q jK q J as the design vector for item j. For α " pα 1 , ..., α K q J , we write α ľ q j , if α k ě q jk for all k P t1, ..., Ku, and write α ń q j , if there exists k such that α k ă q jk .
That is, α ľ q j if profile α has all the skills needed for solving item j and α ń q j if not.
The design information leads to the following constraints: who have mastered all the required skills have the same chance of answering the item correctly.
2. PpY ij " 1 | ξ i " αq ě PpY ij " 1 | ξ i " α 1 q if α ľ q j and α 1 ń q j . 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. PpY ij " 1 | ξ i " αq ě PpY ij " 1 | ξ i " 0q for all α. 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 bŷ When K is relatively large, the computation for solving (2.10) becomes challenging, due to both the summation over 2 K possible values of α in lpβq, and the large number of inequality constraints.

Stochastic Proximal Algorithm
In this section, we propose a quasi-Newton stochastic proximal algorithm for the computation of (2.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 Section 4. In what follows, we describe the algorithm in its general form in Section 3.1, followed by details for two specific models in Sections 3.2 and 3.3, and finally comparisons with related algorithms in Section 3.4.

General Algorithm
For ease of exposition, we introduce some new notations. We write the penalty function as the sum of two terms, Rpβq " R 1 pβq`R 2 pβq, where R 1 pβq is a smooth function and R 2 pβq is non-smooth. In the example of regularized estimation for exploratory IFA, R 1 pβq " 0 and R 2 pβq " λ ř J j"1 ř K k"1 |a jk |, when Rpβq is an L 1 penalty as in (2.6). When an elastic net penalty is used as in (2.7), R 1 pβq " λ 1 ř J j"1 ř K k"1 a 2 jk and R 2 pβq " λ 2 ř J j"1 ř K k"1 |a jk |. The optimization problem can be reexpressed as min β hpβq`gpβq, (3.1) where hpβq "´lpβq`R 1 pβq and g : R p Ñ R Y t`8u is a generalized function taking the form gpβq " R 2 pβq`I B pβq, where (3.2) Note that since both lpβq and R 1 pβq are smooth in β, hpβq is still smooth in β. The second term gpβq is non-smooth in β, unless it is degenerate (i.e., gpβq " 0). We further write Hpξ, βq "´log f py, ξ | βq`R 1 pβq, (3.3) which can be viewed as a complete-data version of hpβq 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 Prox D γ,g pβq " arg min where γ ą 0, D is a strictly positive definite matrix, and }¨} D is a norm defined by The choices of γ, D, 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 Section 3.4 for further comparison.
Update: At tth iteration where t ě 1, we perform the following two steps: 1. Stochastic step: Sample ξ from the conditional distribution of ξ given y, and obtain ξ ptq . The sampling can be either exact or approximated by MCMC.
2. Proximal step: Update model parameters by where T px; c 1 , c 2 q is a truncation function defined as (3.5) Iteratively perform these two steps until a stopping criterion is satisfied and let n be the last iteration number.
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 intuition 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 ξ is exact in the stochastic step, then G ptq is a stochastic gradient of the smooth part of our objective function, in the sense that EpG ptq | y, β pt´1q q " ∇hpβq| β"β pt´1q . If, in addition, there is no constraint or non-smooth penalty, i.e., gpβq " 0, then the proximal step degenerates to an SGD update β ptq " β pt´1q´γ t pD ptq q´1G ptq . 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 nonsmooth 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., Prox D γ,g pβq " β.
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β ptq " β pt´1q´γ t pD ptq q´1G ptq . Thenβ ptq is projected back to the feasible region B by the scaled proximal operator: which is a projection under the norm }¨} D . When D is an identity matrix as in the vanilla (i.e., non-scaled) proximal operator, then the projection is based on the Euclidian distance.
More generally, when the g function involves non-smooth penalties, then the proximal step can be viewed as minimizing the sum of gpβq and a quadratic approximation of hpβq at β pt´1q ; see Lee et al. (2014) for more explanations. We provide an example to facilitate the understanding. Suppose that Then Prox D γ,g pβ ptq q involves solving p optimization problems separately, each of which takes the formβ i " arg min It is well known that (3.6) has a closed-form solution given by soft-thresholding (see Chapter Remark 3 (Role of D ptq ). 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 ptq is used to approximate the Hessian matrix of hpβq at β ptq . When β ptq converges toβ, then δ see Remark 8 for more explanations.
In the proposed update, we choose D ptq to be a diagonal matrix for computational conve-nience. Specifically, as discussed in Remark 2, the proximal step is in a closed form when D ptq is a diagonal matrix. In addition, the proximal step requires to calculate the inverse of D ptq , whose complexity is much lower when D ptq 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 ptq to be non-diagonal. In fact, it is not difficult to generalize the BFGS updating formula for D ptq given in Lee et al. (2014) to a stochastic version.
Our choice of D ptq guarantees its eigenvalues to be constrained in the interval rc 1 , c 2 s. It rules out the singular situation when D ptq is not strictly positive definite. In the implementation, we set c 1 's to be a sufficiently small constant and set c 2 '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 Section 4, the convergence of the proposed method requires the step size to satisfy ř 8 t"1 γ t " 8 and ř 8 t"1 γ 2 t ă 8. This requirement is also needed in the Robbins-Monro algorithm. Here, we choose the step size γ t " µt´1 2´ε so that the above requirement is satisfied, where µ is a positive constant and ε is a small positive constant. As will be shown in Section 4, with sufficiently small ε, β n is almost optimal in terms of its convergence speed. We point out that ε is needed to prove the convergence ofβ n , under our non-convex setting. It is not needed, if the objective function (2.2) is convex; see Atchadé et al. (2017). The requirement of ε may be an artifact due to our proof strategy. Simulation results show that the algorithm converges well even if we set ε " 0. For the numerical analysis in this paper, we set ε " 10´2.
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 γ t " Op1{tq.
Remark 5 (Starting point). As the objective function (2.2) is typically non-convex for most latent variable models, the choice of the starting point β p0q 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 or spectral methods (H. . 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 ξ can be sampled exactly in the stochastic step, then G ptq is a stochastic gradient of hpβq. 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 ptq is only approximately unbiased.
As we show in Section 4, such G ptq may still yield convergence ofβ n .
Remark 7 (Stopping criterion). In the implementation of Algorithm 1, we stop the iterative update by monitoring a window of successive differences in β ptq . 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 Metroplis-Hasting Robins-Monro algorithm proposed by Cai (2010a).
Finally, as we explain in Remark 8, certain quantities, including the Hessian matrix of lpβq, 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 where mpy, ξ | βq is a given function with an analytic form and the conditional expectation E r¨| y, βs is with respect to the conditional distribution of ξ given y. The quantity (3.7) is intractable due to the high-dimensional integral with respect to ξ. One such example is the Hessian matrix of lpβq atβ as discussed in Section 2.2 that is a key quantity for the statistical inference ofβ. In fact, by Louis' formula (Louis, 1982), The computation of (3.7) is a straightforward by-product of the proposed algorithm. To approximateM , we only need to add the following update in each iteration We approximateM by the Polyak-Ruppert averagingM n " p ř n t" `1 M ptq q{pn´ q. When the sequence β ptq converges toβ (see Theorem 4.2 for the convergence analysis), under mild conditions, Theorem 3.17 of Benveniste et al. (1990) suggests the convergence of M pnq toM with probability 1, which further implies the convergence ofM n toM . Note that we use the averaged estimatorM n as it tends to converge faster than the pre-average sequence M pnq . We point out that the updating rule for the diagonal matrix D ptq in Algorithm 1 makes use of such an averaged estimator.
Remark 9 (Burn-in size). Like MCMC algorithms, the proposed method also has a burnin 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 is fixed to be sufficiently large in each of our examples. Adaptive choice of the burn-in size is possible; see S. .

Example I: Item Factor Analysis
We now explain the details of using the proposed method to solve (2.5) for exploratory IFA. The computation is similar when replacing the L 1 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 (2.5). Note that ξ 1 , ..., ξ N are conditionally independent given data, and thus can be sampled separately. For each ξ i , we sample its entries by Gibbs sampling. More precisely, each entry is sampled by adaptive rejection sampling (Gilks & Wild, 1992;, as the conditional distribution of ξ ik given data and the other entries of ξ i is log-concave. We refer the readers to S. Zhang et al. (2020) for more explanations of this sampling procedure. If a normal ogive IFA is considered instead of the logistic model above, then we can sample ξ ptq i 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 β " tB, d j , a j , j " 1, ..., Ju. We denoteβ ptq " β pt´1q´γ t pD ptq q´1G ptq as the input of the scaled proximal operator. The parameter update is given by where the parameter space and gpβq " λ ř J j"1 ř K k"1 |a jk |`I B pβq only involves loading parameters a jk and parameters B for the covariance matrix.
We first look at the update for d j s. As the g function does not involve d j , its update is simply d ptq j "d ptq j , whered ptq j is the corresponding component inβ ptq . We then look at the update for the loading parameters a jk . Suppose that a jk corresponds to the i a jk th component of β. Then the update is given by solving the optimization As discussed in Remark 2, this optimization has a closed-form solution via soft-thresholding.
We finally look at the update for B. Suppose that b kl corresponds to the i b kl th component of β. Then the update of b k , the kth row of B, is given by solving the following optimization which can be easily solved by the method of Lagrangian multiplier.

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 ξ i is still a categorical distribution which can be sampled exactly. Second, the proximal step requires to solve a quadratic programming problem. Again, we denotẽ The proximal step requires to solve the following quadratic programming problem min β pβ´β ptq q J D ptq pβ´β ptq q, and ν 0 " 0. (3.9) 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 Section 5.3, we use the dual method of Goldfarb & Idnani (1983) implemented in the R package quadprog (Turlach et al., 2019).

Comparison with Related Algorithms
We compare Algorithm 1 with several related algorithms in more details.
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., gpβq " 0), 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 γ t 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 (Section 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 Section 4 that the proposed method almost achieve the optimal convergence speed. c 1 " c 2 . 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 opti-mization, while in Section 4 we consider a more general setting of non-convex optimization that includes a wide range of latent variable model estimation problems as special cases.
Stochastic EM algorithm. The proposed method is also closely related to the stochastic-EM algorithm (Celeux, 1985;Ip, 2002;Nielsen, 2000;. 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 β ptq by minimizing the negative complete data log-likelihood function´log f py, ξ ptq | βq, instead of a stochastic gradient update.
It is also recommended to perform a trajectory averaging in the stochastic-EM algorithm (Nielsen, 2000;, like the last step of the proposed algorithm. As pointed out in S. , 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.

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 variable. The following expectation is taken with respect to latent variable ξ given data y and parameters β, denoted by Ep¨| βq " ş¨π β pξqdξ, where π β is the posterior distribution for ξ given y and β. Let }¨} denote the vector l 2 -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 β ptq to a stationary point of the objective function hpβq`gpβq in Theorem 4.2, 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β n in Theorem 4.3 which extends the results of Atchadé et al. (2017) to the setting of non-convex optimization.
For a function f : R d Þ Ñ R Y`8, denote the Fréchet subdifferential (Chapter 8.B Rockafellar & Wets, 1998) of f at the point x by Bf pxq, Define the set of stationary points of the objective function as B˚" tβ P B : D x P Bhpβq`Bgpβq with x J py´βq ě 0, for all y P Bu.
Note that the global minimumβ is a stationary point, i.e.,β P B˚. In addition, when the objective function is smooth, i.e. gpβq " 0, then B˚" tβ P B : ∇hpβq " 0u, which is the standard definition of stationary points set for a smooth function.
The following assumptions are assumed for our objective function.
H4. The stochastic gradient G β pt´1q pξ ptq q is a Monte Carlo approximation of ∇hpβ pt´1q q.
That is, if computationally feasible, we take ξ ptq as an exact sample from π β pt´1q , where, as defined earlier, π β pt´1q is the posterior distribution of ξ given y and β pt´1q . If not, we sample ξ ptq from a Markov kernel P β pt´1q with invariant distribution π β pt´1q .
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 Hpξ,¨q is locally Lipschitzian and weakly convex, respectively. These conditions hold when the completedata log-likelihood function Hpξ,¨q is Lipschitzian and convex on the entire parameter space.
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 4.2 guarantees global convergence.
The convergence of β pnq guarantees the convergence of the Polyak-Ruppert averaging sequenceβ n . However, Theorem 4.2 does not provide information on the convergence speed.
In what follows, we establish the convergence speed ofβ n . Without loss of generality, by Theorem 4.2, we assume that β pnq converges to β˚P B˚.
H7. For β, β 1 P B 1 , any γ ą 0, and diagonal matrix D with diagonal entries δ i P rc 1 , c 2 s, the following conditions hold.
H8. For a measurable function V : Ξ Ñ r1,`8q, a signed measure µ on the σ-field of Ξ, and a function f : Ξ Ñ R, define There exist λ P p0, 1q, b ă 8, m ě 4 and a measurable function W : Ξ Ñ r1,`8q such where G β pξq " BHpξ, βq{Bβ and P β is the Markov kernel defined in condition H4. In addition, for any P p0, ms, there exists C ă 8, ρ P p0, 1q such that for any ξ P Ξ, H9. There exists a constant C such that for any β, β 1 P B 1 , 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 β˚. Specifically, hpβq being strongly convex in B 1 means that there exists a positive constant C, such that p∇hpβq´∇hpβ 1 qq J pβ´β 1 q ě C}β´β 1 } 2 , for any β and β 1 P B 1 . 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 regularization function. Thus, H7 holds  Note that the expectation is taken with respect to ξ p1q , . . . , ξ ptq given β p0q and ξ p0q .
We now provide a few remarks regarding the convergence speed (4.2). First, the small positive constant ε comes from the requirement on step size that ř 8 t"1 γ 2 t ă 8 in H5. Since ř 8 t"1 γ 2 t ă 8 is satisfied when γ t " µt´1 2´ε , for any µ, ε ą 0, the convergence speed of E}β n´β‹ } 2 can be arbitrarily close to Opn´1 2 q by choosing an arbitrarily small ε. Second, this ε 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 and then prove under similar conditions that E}β n´β‹ } 2 ď Cn´1 2 . Lastly, it is well-known that for non-smooth convex optimization, the minimax optimal convergence rate is Opn´1 2 q; see Chapter 3, Nesterov (2004). In this sense, our algorithm is almost minimax optimal, when ε 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).

Study I: Confirmatory IFA
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. USP is the method that we recommend. It has a step size γ t close to t´1 {2 , applies Polyak-Ruppert averaging, and uses a quasi-Newton update in the proximal step.
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 so that it does not involve a quasi-Newton update. c 1 is set to be 1 without tuning in this study.
3. The USP-RM1 method is implemented the same as the USP method, except that β pnq 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 for the covariance matrix where constraints involve.
4. The USP-RM2 method is the same as USP-RM1, except that we set the step size γ t " 1{t which is the asymptotic optimal step size for the Robbins-Monro algorithm (Fabian, 1968).
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 β ptq by completely solving an optimization problem β ptq " arg max βPB Hpξ ptq , βq`gpβq.
In our implementation, this optimization problem is solved by making the quasi-Newton proximal update (3.4) iteratively until convergence.
We consider a confirmatory IFA setting with only two factors (i.e., K " 2), so that an EM algorithm with sufficient numbers of quadrature points and EM steps can be used to obtain a more accurate approximation ofβ 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β rather than the true model parameters. This is because, under suitable conditions, these algorithms converge toβ 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β and the true model parameters) tends to dominate the computational errors (i.e., the difference betweenβ and the results given by the stochastic algorithms).
More precisely, we consider sample size N " 1000 and the number of items J " 20. The design matrix Q 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 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 p0.5, 1.5q. 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) for each method, instead of using an adaptive stopping criterion. For USP, USP-PPG, and StEM, the burn-in size is chosen to be 500. All algorithms are implemented in C and run on the same platform 1 using a single core.
The results regarding the accuracy of the proposed methods are given in  .
Again,â jk is given by the EM algorithm, andã jk is given by one of the five methods. 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 Figure 5.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.

Study II: Exploratory IFA by Regularization
In the second study, we apply the proposed method to regularized estimation for exploratory IFA as discussed in Section 2.2. We consider increasing sample size N " 1000, 2000, 4000, eighty items and five correlated latent factors (i.e., J " 80, K " 5). 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 are drawn i.i.d. from the standard normal distribution, and the non-zero loading parameters a jk are drawn i.i.d. from a uniform distribution over the interval (0.5, 1.5). The elements of covariance matrix Σ " pσ k,k q 5ˆ5 are set to be σ k,k 1 " 1, for k " k 1 and σ k,k 1 " 0.4 for k ‰ k 1 .
For each sample size, 50 independent datasets are generated. In the proposed algorithm, we adopt a burn-in size " 50 and stop based on the criterion discussed in Section 3, where the stopping threshold is set to be 10´3. A decreasing penalty parameter λ N " a log J{N is used to ensure estimation consistency (Chapter 6, Bühlmann & van de Geer, 2011). Other implementation details can be found in Section 3.2. The algorithm in this MSE ofÃ N =1000 N =2000 N =4000 25% quantile 0.032 0.026 0.017 median 0.034 0.027 0.018 75% quantile 0.036 0.027 0.019 Table 4: The mean squared errors for estimated loading parameters in exploratory IFA with L 1 regularization.
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 " pa jk q JˆK . 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 min where }¨} F is the Frobenius norm, A is the true loading matrix,Ã is the output of Algorithm 1, and PpÃq denotes the set of JˆK matrices that can be obtained by swapping the columns ofÃ.
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.

Study III: Restricted LCA
In this study, we apply the proposed method to the estimation of a restricted latent class model as discussed in Section 2.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 (2.10) which is based on a general restricted latent class model considered in Xu (2017) instead of the DINA model, mimicing the practical situation when the parametric form is unknown.
We consider a test consisting of twenty items (i.e., J " 20) that measure four binary attributes (i.e., K " 4). Three sample sizes are considered, including N " 1000, 2000, and 4000. The design matrix Q is given in Table 6. In addition, the guessing and slipping parameters s j and g j 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 θ j,α . That is, Finally, we let ν α " 0, for all α P t0, 1u K , so that Ppξ " αq " 1{2 K . According to the results in Xu (2017), the model parameters are indentifiable, given the Q-matrix in Table 6.
For each sample size, 50 independent datasets are generated. The proposed algorithm adopts a burn-in size " 50 and stops based on the criterion discussed in Section 3, where the stopping threshold is set to be 10´3. Other implementation details can be found in   For structural parameters ν α , the MSE is calculated as 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.

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 γ t . Our theoretical results suggest that γ t " t´0 .5´ for any P p0, 0.5s, and a smaller leads to faster convergence. In practice, we suggest to set γ t " t´0 .51 that performs well in all our simulations. This choice of step size is very different from the choice of γ t " t´1 in the MCMC stochastic approximation algorithms. Second, a burn-in size 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 and c 2 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 to be sufficiently small and c 2 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, 2020b) 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.
Further since B is compact, there is a random variable B which is finite with probability 1, such that for t P N, }β ptq } ď B. Together with step size condition in H5, we have Thus γ t γt pξ ptq , β pt´1q q is a l 2 -summable martingale difference sequence adpated to F t´1 . By standard martingale convergence result (e.g., Thm. 5.3.33, Dembo, 2016), we have with probability 1, lim n ř n t"1 γ t γt pξ ptq , β pt´1q q exist and is finite.

B Proof of Theorem 1
In Theorem 1, we establish the convergence of β ptq to a stationary point β 8 P B˚using differential inclusion techniques in Duchi & Ruan (2018). The proposed method can be viewed as a special case of the general stochastic method discussed in Duchi & Ruan (2018) with a few differences.
With additional assumptions, a similar convergence result can be derived. In what follows, we first show the linear interpolation process of our stochastic updates is asymptotically equivalent to a differential inclusion, by verifying that conditions of Theorem 2 in Duchi & Ruan (2018) hold for our case. Then, cluster points of any trajectory of the limiting differ-ential inclusion are proved to be stationary points. Lastly, the convergence properties of our original sequence can be shown from the functional convergence.
First we define the linear interpolation of the iterates β pkq : βptq " β pkq`t´t k t k`1´tk pβ pk`1q´βpkq q and yptq " ErU γ k pξ pkq ; β pk´1q q | β pk´1q s for t P rt k , t k`1 q, and β t p¨q " βpt`¨q, t P R`be the time-shifted process.
In order to use Theorem 2 of Duchi & Ruan (2018), which is a general functional convergence theorem, conditions (i)-(iv) of Theorem 2 need to be verified for our case. Firstly, the boundness condition (i) holds as B is compact given H1; Non-summable but squaresummable steps size condition (ii) holds given H5; And we have verified (iii), which is the convergence of the summation of the weighted noise sequence, holds by Lemma 1; Lastly, condition (iv) holds similarly in our case for the close-value mapping´U pβq´N B pβq (see Lemma 10 in Duchi & Ruan (2018)), where U pβq " ∇hpβq`Bgpβq and N B pβq " tv P R p : xv, β 1´β y, for all β 1 P Bu is the normal cone for B at β.
Based on Theorem 2 and Theorem 3 of Duchi & Ruan (2018), for any sequences tτ k u 8 k"1 , the function sequence β τ k p¨q is relatively compact in CpR`, R p q and for any τ k Ñ 8, any limit point of tβ τ k p¨qu in CpR`, R p q satisfies βptq "βp0q`ż t 0 ypτ qdτ, where ypτ q P´U pβpτ qq´N B pβpτ qq.
So the sample path of our algorithm is asymptotically equivalent to the differential inclusion  (2018)).
Finally, according to Theorem 1 of Duchi & Ruan (2018), with probability 1, Consequently, given assumption H1, B is compact and B˚contains finite points, we have the objective value F pβ ptq q converges and all cluster points of the sequence tβ ptq u belong to B˚.
By further assumption that different stationary points in B˚have different objective values, we have β ptq converges to a stationary point in B˚, with probability 1.

C Proof of Theorem 4.3
Follow the proofs in Section 6 of Atchadé et al. (2017), we first prove several lemmas, then prove Theorem 4.3.
Lemma C.1. If g is convex and Lipschitz on B 1 with Lipschitz constant K, or g " I B p¨q. For β, β 1 P B 1 , any γ ą 0, and diagonal matrix D with diagonal entries δ i P rc 1 , c 2 s, c 2 ě c 1 ą 0, the following conditions hold.
Proof of Lemma C.1.
Given g is proper convex, Lipschitz on B 1 with Lipschitz constant is K and (i), we have Thus (iii) holds.
where the first inequality follows from Lemma C.1-(ii).
And by induction the proof is concluded.
Lemma C.6. Assume H1, H4, H5 and H8. If a t ě 0, for t ě 1, there exist a constant C such that Proof of Lemma C.6. By Minkowski inequality, as the supremum is finite based on Lemma C.4 and Lemma C.5.
The above inequality follows from assumption H7-(ii).
Proof of Theorem 4.3.

D Additional Simulation Results
We provide an additional simulation study to (1) assess the estimation of the asymptotic variances of parameter estimates and (2) assess the point estimation of the covariance between latent variables.
We consider a similar confirmatory IFA setting as in the simulation study I, with two factors, twenty items (i.e., K " 2, J " 20), and the same design matrix Q. The intercept parameters and non-zero loading parameters are drawn i.i.d. from the standard normal and a uniform distribution over the interval p0.5, 1.5q, respectively. The variances of two factors are set to be 1 and the covariance is set to be 0.4. For each of the three sample sizes N " 1000, 2000, 4000, 50 independent datasets are generated. We then apply the proposed USP method with 1000 burn-in size, and 4000 total iterations. Note that we use a larger burn-in size and a larger number of iterations here to ensure accurate computation of the asymptotic variances, because they tend to be more difficult to compute than the point estimates. The results from the UPS algorithm are compared with those from a standard EM algorithm that uses 31 quadrature points for each dimension.
We approximate the observed Fisher information matrix using the approach given in Remark 8.
Based on the approximated Fisher information matrix, we obtain the standard errors of parameter estimates. The obtained standard errors are compared with those given by the EM algorithm.
The results are given in Figure D.1. Each panel of Figure D.1 corresponds to a combination of a sample size and a type of parameters (loadings/intercepts/covariance). For each dataset and each parameter, we obtain the standard errors of the parameter estimate from the UPS and EM algorithms, respectively. These standard errors are shown as a point in the scatter plot, where the x-axis gives the standard error from the EM algorithm and the y-axis gives the standard error from the USP method. As we can see, all the points concentrate along the diagonal line, suggesting that the standard errors from the two algorithms are very close to each other.
We further assess the estimation of the covariance between the latent variables. The results are given in Figure D.2. For each sample size, we compute the squared difference between the estimate given by the USP algorithm and the true value (σ 12 " 0.4) and visualize the squared errors from the 50 datasets using a box plot. We see that all the squared errors are quite small and they decrease when the sample size increases.