Stochastic Proximal Gradient Algorithms for Penalized Mixed Models

TL;DR

This paper introduces SAPG, a stochastic proximal gradient algorithm for penalized likelihood maximization in complex models where gradients are approximated via Monte Carlo methods, with proven convergence.

Contribution

The paper presents the first convergence analysis of SAPG, a novel stochastic approximation proximal gradient algorithm for non-linear mixed effect models.

Findings

SAPG converges under biased and unbiased Monte Carlo approximations.

Compared to existing methods, SAPG effectively handles high-dimensional non-linear mixed models.

Application to pharmacokinetic data demonstrates practical utility.

Abstract

Motivated by penalized likelihood maximization in complex models, we study optimization problems where neither the function to optimize nor its gradient have an explicit expression, but its gradient can be approximated by a Monte Carlo technique. We propose a new algorithm based on a stochastic approximation of the Proximal-Gradient (PG) algorithm. This new algorithm, named Stochastic Approximation PG (SAPG) is the combination of a stochastic gradient descent step which - roughly speaking - computes a smoothed approximation of the past gradient along the iterations, and a proximal step. The choice of the step size and the Monte Carlo batch size for the stochastic gradient descent step in SAPG are discussed. Our convergence results cover the cases of biased and unbiased Monte Carlo approximations. While the convergence analysis of the Monte Carlo-PG is already addressed in the literature (see Atchad\'e et al. [2016]), the convergence analysis of SAPG is new. The two algorithms are compared on a linear mixed effect model as a toy example. A more challenging application is proposed on non-linear mixed effect models in high dimension with a pharmacokinetic data set including genomic covariates. To our best knowledge, our work provides the first convergence result of a numerical method designed to solve penalized Maximum Likelihood in a non-linear mixed effect model.

Figures

Equations

$${\mathrm{argmax}}_{\theta \in {\mathbb{R}}^d}F(\theta ),\text{where}F:=\ell -g,$$

$$\begin{array}{rl}\nabla \ell (\theta )& =\nabla \varphi (\theta )+\Psi (\theta )\overline{S}(\theta ),\\ & \text{with}\overline{S}(\theta ):={\int }_{\mathsf{Z}}S(z){\pi }_{\theta }(z)\nu (\mathrm{d}z);\end{array}$$

$$\Vert \nabla \ell (\theta )-\nabla \ell ({\theta }^{\prime })\Vert \le L\Vert \theta -{\theta }^{\prime }\Vert ;$$

$$\theta \mapsto \ell (\theta ):=\log {\int }_{\mathsf{Z}}\exp \left(\varphi (\theta )+⟨S(z),\psi (\theta )⟩\right)\nu (\mathrm{d}z),$$

$$\begin{array}{rl}\nabla \ell (\theta )& =\nabla \varphi (\theta )+J\psi (\theta ){\int }_{\mathsf{Z}}S(z){\pi }_{\theta }(z)\nu (\mathrm{d}z),\\ & \text{with}{\pi }_{\theta }(z):=\frac{\exp (⟨S(z),\psi (\theta )⟩)}{{\int }_{\mathsf{Z}}\exp (⟨S(u),\psi (\theta )⟩)\nu (\mathrm{d}u)},\end{array}$$

$$\begin{array}{rl}\theta \mapsto & \ell (\theta ):=\\ & \sum _{k=1}^N⟨S({\mathsf{Y}}_k),\theta ⟩-N\log {\int }_{\mathsf{Z}}\exp \left(⟨S(z),\theta ⟩\right)\nu (\mathrm{d}z).\end{array}$$

$$\begin{array}{rl}\nabla \ell (\theta )& =\sum _{k=1}^{N}S({\mathsf{Y}}_k)-N{\int }_{\mathsf{Z}}S(z){\pi }_{\theta }(z)\nu (\mathrm{d}z),\\ & \text{with}{\pi }_{\theta }(z):=\frac{\exp \left(⟨S(z),\theta ⟩\right)}{\int \exp \left(⟨S(u),\theta ⟩\right)\nu (\mathrm{d}u)}.\end{array}$$

$${\mathsf{Y}}_k=f(t_k,Z^{(k)})+{\epsilon }_k,$$

$${\mathrm{Prox}}_{\gamma ,g}(\theta ):={\mathrm{argmin}}_{\tau \in \Theta }\left\{g(\tau )+\frac{1}{2\gamma }\Vert \theta -\tau {\Vert }^2\right\}.$$

$$g(\theta):=\left\{\begin{array}[]{ll}0&\text{if $\theta\in\Theta$}\\ +\infty&\text{otherwise},\end{array}\right.$$

$$\left(\operatorname{Prox}_{\gamma,g_{\lambda,\alpha}}(\theta)\right)_{i}\\ =\frac{1}{1+\gamma\lambda(1-\alpha)}\left\{\begin{array}[]{ll}0&\text{if $|\theta_{i}|\leq\gamma\lambda\alpha$,}\\ \theta_{i}-\gamma\lambda\alpha&\text{if $\theta_{i}\geq\gamma\lambda\alpha$,}\\ \theta_{i}+\gamma\lambda\alpha&\text{if $\theta_{i}\leq-\gamma\lambda\alpha$.}\end{array}\right.$$

$$\begin{array}{rl}{\theta }_{n+1}& ={\mathrm{Prox}}_{{\gamma }_{n+1},g}\left({\theta }_n+{\gamma }_{n+1}\nabla \ell ({\theta }_n)\right)\\ & ={\mathrm{Prox}}_{{\gamma }_{n+1},g}\left({\theta }_n+{\gamma }_{n+1}\{\nabla \varphi ({\theta }_n)+\Psi ({\theta }_n)\overline{S}({\theta }_n)\}\right).\end{array}$$

$$\begin{array}{rl}& F(\theta )=\ell (\theta )-g(\theta )\\ & \ge \ell ({\theta }_n)-⟨\nabla \ell ({\theta }_n),\theta -{\theta }_n⟩-\frac{1}{2{\gamma }_{n+1}}\Vert \theta -{\theta }_n{\Vert }^2-g(\theta ).\end{array}$$

$$S_{n+1}^{\mathrm{mc}}:=\frac{1}{m_{n+1}}\sum _{j=1}^{m_{n+1}}S(Z_{j,n}).$$

$$S_{n+1}^{\mathrm{sa}}:=(1-{\delta }_{n+1})S_n^{\mathrm{sa}}+\frac{{\delta }_{n+1}}{m_{n+1}}\sum _{j=1}^{m_{n+1}}S(Z_{j,n}),$$

$$\mathcal{Q}(\theta |{\theta }^{\prime }):=\varphi (\theta )+⟨\overline{S}({\theta }^{\prime }),\psi (\theta )⟩.$$

$${\theta }_{n+1}={\mathrm{argmax}}_{\theta \in \Theta }\left\{\varphi (\theta )+⟨\overline{S}({\theta }_n),\psi (\theta )⟩-g(\theta )\right\}.$$

$$\begin{array}{rl}{\theta }_{n+1}\text{s.t.}& \varphi ({\theta }_{n+1})+⟨\overline{S}({\theta }_n),\psi ({\theta }_{n+1})⟩-g({\theta }_{n+1})\\ & \ge \varphi ({\theta }_n)+⟨\overline{S}({\theta }_n),\psi ({\theta }_n)⟩-g({\theta }_n).\end{array}$$

$$\begin{gathered} \Vert \nabla \varphi (\theta )-\nabla \varphi ({\theta }^{\prime })+\left(J\psi (\theta )-J\psi ({\theta }^{\prime })\right)s\Vert \\ \le L\Vert \theta -{\theta }^{\prime }\Vert . \end{gathered}$$

$${\theta }_{n+1}={\mathrm{Prox}}_{{\gamma }_{n+1},g}\left({\theta }_n+{\gamma }_{n+1}\left\{\nabla \varphi ({\theta }_n)+\Psi ({\theta }_n)S_{n+1}\right\}\right),$$

$$\sum _n{\gamma }_{n+1}⟨\Psi ({\theta }_n)\left(S_{n+1}-\overline{S}({\theta }_n)\right),T_{{\gamma }_{n+1}}({\theta }_n)⟩,$$

$$\sum _n{\gamma }_{n+1}\Psi ({\theta }_n)\left(S_{n+1}-\overline{S}({\theta }_n)\right),$$

$$\sum _n{\gamma }_{n+1}^2\Vert \Psi ({\theta }_n)\left(S_{n+1}-\overline{S}({\theta }_n)\right){\Vert }^2,$$

$$T_{\gamma }(\theta ):={\mathrm{Prox}}_{\gamma ,g}(\theta +\gamma \left\{\nabla \varphi (\theta )+\Psi (\theta )\overline{S}(\theta )\right\}),$$

$${\mathcal{F}}_n:=\sigma \left(Z_{i,j},i\le m_{j+1},j\le n-1\right),$$

$$\mathbb{E}\left[S_{n+1}|{\mathcal{F}}_n\right]=\overline{S}({\theta }_n).$$

$$\begin{gathered} \Vert \nabla \varphi (\theta )-\nabla \varphi ({\theta }^{\prime })\Vert +\Vert \Psi (\theta )-\Psi ({\theta }^{\prime })\Vert +\Vert \overline{S}(\theta )-\overline{S}({\theta }^{\prime })\Vert \\ \le L\Vert \theta -{\theta }^{\prime }\Vert , \end{gathered}$$

$$P^0(x,A):=\text{mathbbm}{1}_A(x),$$

$$P^k(x,A):=\int P^{k-1}(x,\mathrm{d}z)P(z,A),k\ge 1.$$

$$\xi P(A):=\int \xi (\mathrm{d}z)P(z,A),A\in \mathcal{Z};$$

$$Pf(z):=\int f(y)P(z,\mathrm{d}y).$$

$${\left|f\right|}_W:=\underset{\mathsf{Z}}{\sup }\frac{\Vert f\Vert }{W},{\Vert \nu \Vert }_W:=\underset{f:{|f|}_W\le 1}{\sup }\left|\int f\mathrm{d}\nu \right|;$$

Full text

Keywords: Proximal-Gradient algorithm; Stochastic Gradient; Stochastic EM algorithm; Stochastic Approximation; Non-linear mixed effect models.

1 Introduction

Many problems in computational statistics reduce to the maximization of a criterion

[TABLE]

and the functions $\ell,g$ satisfy

H 1**.**

the function $g:\;\mathbb{R}^{d}\to\left[0,+\infty\right]$ is convex, not identically $+\infty$, and lower semi-continuous.

H 2**.**

the function $\ell:\mathbb{R}^{d}\to\mathbb{R}\cup\{-\infty\}$ is continuously differentiable on $\Theta:=\{\theta\in\mathbb{R}^{d}:g(\theta)+|\ell(\theta)|<\infty\}$ and its gradient is of the form

[TABLE]

$\nabla$* denotes the gradient operator and $\pi_{\theta}\mathrm{d}\nu$ is a probability distribution on a measurable subset $(\mathsf{Z},\mathcal{Z})$ of $\,\mathbb{R}^{p}$. The measurable functions $\nabla\phi:\mathbb{R}^{d}\to\mathbb{R}^{d}$ and $\Psi:\mathbb{R}^{d}\to\mathbb{R}^{d\times q}$ are known but the expectation $\overline{S}$ of the function $S:\mathsf{Z}\to\mathbb{R}^{q}$ with respect to $\pi_{\theta}\mathrm{d}\nu$ may be intractable. Furthermore, there exists a finite non-negative constant $L$ such that for all $\theta,\theta^{\prime}\in\Theta$,*

[TABLE]

$\|\cdot\|$* is the Euclidean norm.*

Examples of functions $\ell$ satisfying Eq. (2) are given below. We are interested in numerical methods for solving Eq. (1), robust to the case when neither $\ell$ nor its gradient have an explicit expression.

Such an optimization problem occurs for example when computing a penalized maximum likelihood estimator in some parametric model indexed by $\theta\in\mathbb{R}^{d}$: $\ell$ denotes the log-likelihood of the observations $\mathsf{Y}$ (the dependence upon $\mathsf{Y}$ is omitted) and $g$ is the penalty term.

The optimization problem Eq. (1) covers the computation of the maximum when the parameter $\theta$ is restricted to a closed convex subset $\Theta$ of $\mathbb{R}^{d}$; in that case, $g$ is the characteristic function of $\Theta$ i.e. $g(\theta)=0$ for any $\theta\in\Theta$ and $g(\theta)=+\infty$ otherwise. It also covers the case when $g$ is the ridge, the lasso or the elastic net penalty; and more generally, the case when $g$ is the sum of lower semi-continuous non-negative convex functions.

A first example of such a function $\ell$ is given by the log-likelihood in a latent variable model with complete likelihood from the $q$-parameter exponential family (see e.g. Bickel and Doksum (2015) and Bartholomew et al. (2011) and the references therein). In that case, $\ell$ is of the form

[TABLE]

where $\left<a,b\right>$ denotes the scalar of two vectors $a,b\in\mathbb{R}^{l}$, $\phi:\mathbb{R}^{d}\to\mathbb{R}$, $\psi:\mathbb{R}^{d}\to\mathbb{R}^{q}$ and $S:\mathsf{Z}\to\mathbb{R}^{q}$ are measurable functions, and $\nu$ is a $\sigma$-finite positive measure on $(\mathsf{Z},\mathcal{Z})$. The quantity $\theta\mapsto\phi(\theta)+\left<S(Z),\psi(\theta)\right>$ is known as the complete log-likelihood, and $Z$ is the latent data vector. Under regularity conditions, we have

[TABLE]

where $\operatorname{J}\psi(\theta)$ denotes the transpose of the jacobian matrix of the function $\psi$ at $\theta$.

A second example is given by the log-likelihood of $N$ independent observations $(\mathsf{Y}_{1},\cdots,\mathsf{Y}_{N})$ from a log-linear model for Markov random fields. In this model, $\ell$ is given by

[TABLE]

The function $\theta\mapsto\int_{\mathsf{Z}}\exp\left(\left<S(z),\theta\right>\right)\nu(\mathrm{d}z)$ is known as the partition function. Under regularity conditions, we have

[TABLE]

In these two examples, the integrals in Eqs. (4) to (7) are intractable except for toy examples: neither the function $\ell$ nor its gradient are available. Nevertheless, all the integrals in Eqs.  (4)-(7) can be approximated by a Monte Carlo sum (see e.g. Robert and Casella, 2004). In the first example, this Monte Carlo approximation consists in imputing the missing variables $z$; it is known that such an imputation is far more efficient when the Monte Carlo samples are drawn under $\pi_{\theta}\mathrm{d}\nu$, i.e. the a posteriori distribution of the missing variables given the observations (see Eq. (5)) than when they are drawn under the a priori distribution. This remark is the essence of the Expectation Maximization (EM) algorithm (introduced in Dempster et al., 1977), a popular iterative procedure for maximizing the log-likelihood $\ell$ in latent variable models.

In this paper, we are interested in first order optimization methods to solve Eq. (1), that is methods based on the gradient. In Section 2.1, we describe two stochastic first-order descent methods, which are stochastic perturbations of the Proximal-Gradient (PG) algorithm (introduced in  Combettes and Pesquet (2011); see also  Beck and Teboulle (2009); Parikh and Boyd (2013) for literature reviews on Proximal-Gradient algorithms). The two algorithms are the Monte Carlo Proximal-Gradient algorithm (MCPG) and the Stochastic Approximation Proximal-Gradient algorithm (SAPG), which differ in the approximation of the gradient $\nabla\ell$ and more precisely, of the intractable integral $\overline{S}(\theta)$ (see Eq. (2)). In MCPG, at each iteration $n$ of the algorithm, this expectation evaluated at the current point $\theta_{n}$ is approximated by a Monte Carlo sum computed from samples $\{Z_{1,n},\cdots,Z_{m_{n+1},n}\}$ approximating $\pi_{\theta_{n}}\mathrm{d}\nu$. In SAPG, the approximation is computed as a Monte Carlo sum based on all the points drawn during all the previous iterations of the algorithm $\{Z_{i,j},i\leq m_{j+1},j\leq n\}$.

When $\ell$ is the log-likelihood of a latent variable model, we prove in Section 2.2 that our algorithms are Generalized EM algorithms (see e.g. McLachlan and Krishnan, 2008; Ng et al., 2012) combined with a stochastic E-step: in MCPG and SAPG, the stochastic E-step mimics respectively the E-step of the Monte Carlo EM (Wei and Tanner, 1990; Levine and Fan, 2004) and the E-step of the Stochastic Approximation EM (see e.g. Delyon et al., 1999).

Section 3 is devoted to the convergence analysis of MCPG and SAPG. These algorithms can be seen as perturbed Proximal-Gradient algorithms when the perturbation comes from replacing the exact quantity $\overline{S}(\theta_{n})$ by a Monte Carlo approximation $S_{n+1}$ at each iteration of the algorithm. Our convergence analysis covers the case when the points $\{Z_{1,n},\cdots,Z_{m_{n+1},n}\}$ are sampled from a Markov chain Monte Carlo sampler (MCMC) with target distribution $\pi_{\theta_{n}}\mathrm{d}\nu$ - and therefore, it also covers the case of i.i.d. draws. This implies that the estimator $S_{n+1}$ of $\overline{S}(\theta_{n})$ may be biased. There exist many contributions in the literature on the convergence of perturbed Proximal-Gradient algorithms when $\ell$ is concave, but except in the works by Atchadé et al. (2017) and Combettes and Pesquet (2015), most of them assume that the error $S_{n+1}-\overline{S}(\theta_{n})$ is unbiased and gets small when $n\to\infty$ (see e.g. Rosasco et al., 2014; Combettes and Pesquet, 2016; Rosasco et al., 2016; Lin et al., 2015). In this paper, we provide sufficient conditions for the almost-sure convergence of MCPG and SAPG under the assumption that $\ell$ is concave and with no assumptions on the bias of $S_{n+1}-\overline{S}(\theta_{n})$. The convergence analysis of MCPG is a special case of (Atchadé et al., 2017, Section 4); to our best knowledge, the convergence of SAPG is a new result.

Practical implementation is discussed in Section 4. Some guidelines are given in Section 4.2 to choose the sequences involved in the stochastic approximation procedures. Then, MCPG and SAPG are compared through a toy example in Section 4.3. A more challenging application to penalized inference in a mixed effect model is detailed in Section 5. Mixed models are applied to analyze repeated data in a population of subjects. The $N$ independent vectors of observations $(\mathsf{Y}_{k},k=1,\ldots,N)$ of the $N$ subjects are modeled by

[TABLE]

with individual latent variable $Z^{(k)}$ independent of the measurement error vector $\varepsilon_{k}$ and $f$ the regression function that depends on the vector of observation times $t_{k}$. Mixed models thus enter the class of models given by Eq. (4) with latent variables $Z=(Z^{(1)},\ldots,Z^{(N)})$. When a covariate model is introduced, the number of covariates can be large, but with only a few of them being influential. This is a sparse estimation problem and the selection problem can be treated through the optimization of a penalized version of the log-likelihood Eq. (4). In non-linear mixed models, the optimization problem is not explicit and stochastic penalized versions of EM (Bertrand and Balding, 2013; Ollier et al., 2016; Chen et al., 2017) have been proposed. To our best knowledge, stochastic Proximal-Gradient algorithms have not been proposed for mixed models.

2 Stochastic Proximal-Gradient based algorithms

In this section, we describe first-order based algorithms for solving Eq. (1) under the assumptions H1 and H2, when the expectation $\overline{S}(\theta)$ in Eq. (2) is intractable.

2.1 The MCPG and SAPG algorithms

Both MCPG and SAPG are iterative algorithms, each update relies on the combination of a gradient step and a proximal operator. The proximal map (Moreau (1962), see also Bauschke and Combettes (2011); Parikh and Boyd (2013)) associated to a convex function $g$ is defined for any $\gamma>0$ and $\theta\in\mathbb{R}^{d}$ by

[TABLE]

Note that under H1, for any $\gamma>0$ and $\theta\in\mathbb{R}^{d}$, there exists an unique point $\tau$ minimizing the RHS of Eq. (9). This proximal operator may have an explicit expression. When $g$ is the characteristic function

[TABLE]

for some closed convex set $\Theta\subseteq\mathbb{R}^{d}$, then $g(\theta)$ is the projection of $\theta$ on $\Theta$. This projection is explicit for example when $\Theta$ is an hyper-rectangle. Another example of explicit proximal operator is the case associated to the so-called elastic net penalty i.e. $g_{\lambda,\alpha}(\theta):=\lambda\left(\frac{1-\alpha}{2}\sum_{i=1}^{d}\theta_{i}^{2}\right.$ $\left.+\alpha\sum_{i=1}^{d}|\theta_{i}|\right)$ with $\theta=(\theta_{1},\cdots,\theta_{d})$, $\lambda>0$ and $\alpha\in\left(0,1\right]$, then for any component $i\in\{1,\cdots,d\}$,

[TABLE]

The Proximal-Gradient algorithm for solving the optimization problem Eq. (1) produces a sequence $\{\theta_{n},n\geq 0\}$ as follows: given a $\left(0,1/L\right]$-valued sequence $\{\gamma_{n},n\geq 0\}$,

[TABLE]

This update scheme can be explained as follows: by H2, we have for any $L\leq\gamma_{n+1}^{-1}$,

[TABLE]

This minorizing function is equal to $F(\theta_{n})$ at the point $\theta_{n}$; the maximization (w.r.t. $\theta$) of the RHS yields $\theta_{n+1}$ given by Eq.  (10). The Proximal-Gradient algorithm is therefore a Minorize - Majorization (MM) algorithm and the ascent property holds: $F(\theta_{n+1})\geq F(\theta_{n})$ for all $n$. Sufficient conditions for the convergence of the Proximal-Gradient algorithm Eq. (10) can be derived from the results by Combettes and Wajs (2005); Parikh and Boyd (2013) or from convergence analysis of MM algorithms (see e.g. Zangwill, 1969; Meyer, 1976).

In the case $\overline{S}(\theta)$ can not be computed, we describe two strategies for a Monte Carlo approximation. At iteration $n+1$, given the current value of the parameter $\theta_{n}$, $m_{n+1}$ points $\{Z_{1,n},\cdots,Z_{m_{n+1},n}\}$ from the path of a Markov chain with target distribution $\pi_{\theta_{n}}\mathrm{d}\nu$ are sampled. A first strategy consists in replacing $\overline{S}(\theta_{n})$ by a Monte Carlo mean:

[TABLE]

A second strategy, inspired by stochastic approximation methods (see e.g. Benveniste et al., 1990; Kushner and Yin, 2003) consists in replacing $\overline{S}(\theta_{n})$ by a stochastic approximation

[TABLE]

where $\{\delta_{n},n\geq 0\}$ is a deterministic $\left[0,1\right]$-valued sequence. These two strategies yield respectively the Monte Carlo Proximal-Gradient (MCPG) algorithm (see Algorithm 1) and the Stochastic Approximation Proximal- Gradient (SAPG) algorithm (see Algorithm 2).

In Section 3, we prove the convergence of MCPG to the maximum points of $F$ when $\ell$ is concave, for different choices of the sequences $\{\gamma_{n},m_{n},n\geq 0\}$ including decreasing or constant step sizes $\{\gamma_{n},n\geq 0\}$ and respectively, constant or increasing batch size $\{m_{n},n\geq 0\}$. We also establish the convergence of SAPG to the maximum points (in the concave case); only the case of a constant batch size $\{m_{n},n\geq 0\}$ and a decreasing step size $\{\gamma_{n},n\geq 0\}$ is studied, since this framework corresponds to the Stochastic Approximation one from which the update rule Eq. (12) is inherited (see details in Delyon et al., 1999). From a numerical point of view, the choice of the sequences $\{\gamma_{n},n\geq 0\}$, $\{\delta_{n},n\geq 0\}$ and $\{m_{n},n\geq 0\}$ is discussed in Section 4: guidelines are given in Section 4.2 and the behavior of the algorithm is illustrated through a toy example in Section 4.3.

2.2 Case of latent variable models from the exponential family

In this section, we consider the case when $\ell$ is given by Eq. (4). A classical approach to solve penalized maximum likelihood problems in latent variables models with complete likelihood from the exponential family is the Expectation-Maximization (EM) algorithm or a generalization called the Generalized EM (GEM) algorithm (Dempster et al., 1977; McLachlan and Krishnan, 2008; Ng et al., 2012). Our goal here, is to show that MCPG and SAPG are stochastic perturbations of a GEM algorithm.

The EM algorithm is an iterative algorithm: at each iteration, given the current parameter $\theta_{n}$, the quantity $\mathcal{Q}(\theta|\theta_{n})$, defined as the conditional expectation of the complete log-likelihood under the a posteriori distribution for the current fit of the parameters, is computed:

[TABLE]

The EM sequence $\{\theta_{n},n\geq 0\}$ for the maximization of the penalized log-likelihood $\ell-g$ is given by (see McLachlan and Krishnan, 2008, Section 1.6.1.)

[TABLE]

When $\overline{S}(\theta)$ is intractable, it was proposed to replace $\overline{S}(\theta_{n})$ in this EM-penalized algorithm by an approximation $S_{n+1}$

• see Algorithm 3. When $S_{n+1}=S_{n+1}^{mc}$ (see Eq. (11)), this yields the so-called Monte Carlo-EM penalized algorithm (MCEM-pen), trivially adapted from MCEM proposed by Wei and Tanner (1990); Levine and Fan (2004). Another popular strategy is to replace $\overline{S}(\theta_{n})$ by $S_{n+1}^{\mathrm{sa}}$ (see Eq. (12)) yielding to the so-called Stochastic Approximation-EM penalized algorithm (SAEM-pen) - (see Delyon et al. (1999) for the unpenalized version).

When the maximization of Eq. (13) is not explicit, the update of the parameter is modified as follows, yielding the Generalized EM-penalized algorithm (GEM-pen):

[TABLE]

This update rule still produces a sequence $\{\theta_{n},n\geq 0\}$ satisfying the ascent property $F(\theta_{n+1})\geq F(\theta_{n})$ which is the key property for the convergence of EM (see e.g. Wu, 1983). Here again, the approximations defined in Eq. (11) and Eq. (12) can be plugged in the GEM-pen update Eq. (15) when $\overline{S}$ is not explicit.

We show in the following proposition that the sequence $\{\theta_{n},n\geq 0\}$ produced by the Proximal-Gradient algorithm Eq. (10) is a GEM-pen sequence since it satisfies the inequality Eq. (15). As a consequence, MCPG and SAPG are stochastic GEM-pen algorithms.

Proposition 1**.**

Let $g$ satisfying H1 and $\ell$ be of the form Eq. (4) with continuously differentiable functions $\phi:\mathbb{R}^{d}\to\mathbb{R}$, $\psi:\mathbb{R}^{d}\to\mathbb{R}^{q}$ and $S:\mathsf{Z}\to\mathbb{R}^{q}$. Set $\Theta:=\{g+|\ell|<\infty\}$. Define $\overline{S}:\Theta\to\mathbb{R}^{q}$ by $\overline{S}(\theta):=\int_{\mathsf{Z}}S(z)\,\pi_{\theta}(z)\,\nu(\mathrm{d}z)$ where $\pi_{\theta}$ is given by Eq. (5). Assume that there exists a constant $L>0$ such that for any $s\in\overline{S}(\Theta)$, and any $\theta,\theta^{\prime}\in\Theta$,

[TABLE]

Let $\{\gamma_{n},n\geq 0\}$ be a (deterministic) positive sequence such that $\gamma_{n}\in\left(0,1/L\right]$ for all $n\geq 0$.

Then the Proximal-Gradient algorithm Eq. (10) is a GEM-pen algorithm for the maximization of $\ell-g$.

The proof is postponed in Appendix A. The assumption Eq. (16) holds when $\Theta$ is compact and $\overline{S}$ (resp. $\phi$ and $\psi$) are continuous (resp. twice continuously differentiable). Note also that for any $\theta$, $\theta^{{}^{\prime}}\in\Theta$ and $s\in\overline{S}(\Theta)$, we have $\left(\operatorname{J}\psi(\theta)-\operatorname{J}\psi(\theta^{\prime})\right)s=0$ if $\overline{S}(\theta)\in\mathrm{Ker}(\operatorname{J}\psi(\theta^{\prime}))$ for any $\theta$, $\theta^{{}^{\prime}}\in\Theta$.

3 Convergence of MCPG and SAPG

The convergence of MCPG and SAPG is established by applying recent results from Atchadé et al. (2017) on the convergence of perturbed Proximal-Gradient algorithms. (Atchadé et al., 2017, Theorem 2) applied to the case $\nabla\ell(\theta)$ is of the form $\nabla\phi(\theta)+\Psi(\theta)\overline{S}(\theta)$, where $\overline{S}(\theta)$ is an intractable expectation and $\nabla\phi,\Psi$ are explicit, yields

Theorem 2**.**

Assume H1, H2, $\theta\mapsto\ell(\theta)$ is concave, and the set $\mathcal{L}:=\operatorname{argmax}_{\theta\in\Theta}F(\theta)$ is a non empty subset of $\Theta$. Let $\{\theta_{n},n\geq 0\}$ be given by

[TABLE]

with a $\left(0,1/L\right]$-valued stepsize sequence $\{\gamma_{n},n\geq 0\}$ satisfying $\sum_{n}\gamma_{n}=+\infty$. If the series

[TABLE]

converge, where

[TABLE]

then there exists $\theta_{\infty}\in\mathcal{L}$ such that $\lim_{n}\theta_{n}=\theta_{\infty}$.

We check the conditions of Theorem 2 in the case $S_{n+1}$ is resp. given by Eq. (11) for the proof of MCPG and by Eq. (12) for the proof of SAPG. Our convergence analysis is restricted to the case $\ell$ is concave; to our best knowledge, the convergence of the perturbed Proximal-Gradient algorithms when $\ell$ is not concave is an open question.

The novelty in this section is Proposition 5 and Theorem 6 which provide resp. a control of the $L_{2}$-norm of the error $S_{n+1}^{\mathrm{sa}}-\overline{S}(\theta_{n})$ and the convergence of SAPG. These results rely on a rewriting of $\left(S_{n+1}^{\mathrm{sa}}-\overline{S}(\theta_{n})\right)$ taking into account that $S_{n+1}^{\mathrm{sa}}$ is a weighted sum of the function $S$ evaluated at all the samples $\{Z_{i,j},i\leq m_{j+1},j\leq n\}$ drawn from the initialization of the algorithm. This approximation differs from a more classical Monte Carlo approximation (see Theorems 3 and 4 for the convergence of MCPG, which are special cases of the results in Atchadé et al. (2017)).

We allow the simulation step of MCPG and SAPG to rely on a Markov chain Monte Carlo sampling: at iteration $(n+1)$, the conditional distribution of $Z_{j+1,n}$ given the past is $P_{\theta_{n}}(Z_{j,n},\cdot)$ where $P_{\theta}$ is a Markov transition kernel having $\pi_{\theta}\mathrm{d}\nu$ as its unique invariant distribution. The control of the quantities $S_{n+1}-\overline{S}(\theta_{n})$ requires some ergodic properties on the kernels $\{P_{\theta_{n}},n\geq 0\}$ along the path $\{\theta_{n},n\geq 0\}$ produced by the algorithm. These properties have to be uniform in $\theta$, a property often called the “containment condition” (see e.g. the literature on the convergence of adaptive MCMC samplers, for example Andrieu and Moulines (2006); Roberts and Rosenthal (2007); Fort et al. (2011b)). There are therefore three main strategies to prove the containment condition. In the first strategy, $\Theta$ is assumed to be bounded, and a uniform ergodic assumption on the kernels $\{P_{\theta},\theta\in\Theta\}$ is assumed. In the second one, there is no boundedness assumption on $\Theta$ but the property $\mathbb{P}(\limsup_{n}\|\theta_{n}\|<\infty)=1$ has to be established prior the proof of convergence; a kind of local boundedness condition on the sequence $\{\theta_{n},n\geq 0\}$ is then applied - see e.g. Andrieu and Moulines (2006); Fort et al. (2011b). The last strategy consists in showing that $\mathbb{P}(\sup_{n}\rho_{n}\|\theta_{n}\|<\infty)=1$ for some deterministic sequence $\{\rho_{n},n\geq 0\}$ vanishing to zero when $n\to\infty$ at a rate compatible with the decaying ergodicity rate - see e.g. Saksman and Vihola (2010). The last two strategies are really technical and require from the reader a strong background on controlled Markov chain theory; for pedagogical purposes, we therefore decided to state our results in the first context: we will assume that $\Theta$ is bounded.

By allowing MCMC approximations, we propose a theory which covers the case of a biased approximation, called below the biased case: conditionally to the past

[TABLE]

the expectation of $S_{n+1}$ is not $\overline{S}(\theta_{n})$: $\mathbb{E}\left[S_{n+1}|\mathcal{F}_{n}\right]\neq\overline{S}(\theta_{n})$. As soon as the samplers $\{P_{\theta},\theta\in\Theta\}$ are ergodic enough (for example, under H4a) and H4b)), the bias vanishes when the number of Monte Carlo points $m_{n}$ tends to infinity. Therefore, the proof for the biased case when the sequence $\{m_{n},n\geq 0\}$ is constant is the most technical situation since the bias does not decay. It relies on a specific decomposition of the error $S_{n+1}-\overline{S}(\theta_{n})$ into a martingale increment with bounded $L^{2}$-moments, and a remainder term which vanishes when $n\to\infty$ even when the batch size $m_{n}$ is constant. Such a behavior of the remainder term is a consequence of regularity properties on the functions $\nabla\phi$, $\Psi$, $\overline{S}$ (see H3c)), on the proximity operator (see H3d)) and on the kernels $\{P_{\theta},\theta\in\Theta\}$ (see H4c)).

Our theory also covers the unbiased case i.e. when

[TABLE]

We therefore establish the convergence of MCPG and SAPG by strengthening the conditions H1 and H2 with

H 3**.**

1. a)

$\ell$* is concave and the set $\mathcal{L}:=\operatorname{argmax}_{\Theta}F$ is a non-empty subset of $\Theta$.* 2. b)

$\Theta$* is bounded.* 3. c)

There exists a constant $L$ such that for any $\theta,\theta^{\prime}\in\Theta$,

[TABLE]

where for a matrix $A$, $\|A\|$ denotes the operator norm associated with the Euclidean vector norm. 4. d)

$\sup_{\gamma\in\left(0,1/L\right]}\sup_{\theta\in\Theta}\gamma^{-1}\|\operatorname{Prox}_{\gamma,g}(\theta)-\theta\|<\infty$.

Note that the assumptions H3b)-H3c) imply Eq. (3) and $\sup_{\theta\in\Theta}\left(\|\nabla\phi(\theta)\|+\|\Psi(\theta)\|+\|\overline{S}(\theta)\|\right)<\infty$. When $\Theta$ is a compact convex set, then H3d) holds for the elastic net penalty, the Lasso or the fused Lasso penalty. (Atchadé et al., 2017, Proposition 11) gives general conditions for H3d) to hold.

Before stating the ergodicity conditions on the kernels $\{P_{\theta},\theta\in\Theta\}$, let us recall some basic properties on Markov kernels. A Markov kernel $P$ on the measurable set $(\mathsf{Z},\mathcal{Z})$ is an application on $\mathsf{Z}\times\mathcal{Z}$, taking values in $\left[0,1\right]$ such that for any $x\in\mathsf{Z}$, $P(x,\cdot)$ is a probability measure on $\mathcal{Z}$; and for any $A\in\mathcal{Z}$, $x\mapsto P(x,A)$ is measurable. Furthermore, if $P$ is a Markov kernel, $P^{k}$ denotes the $k$-th iterate of $P$ defined by induction as

[TABLE]

Finally, the kernel $P$ acts on the probability measures: for any probability measure $\xi$ on $\mathcal{Z}$, $\xi P$ is a probability measure defined by

[TABLE]

and $P$ acts on the positive measurable functions: for a measurable function $f:\mathsf{Z}\to\mathbb{R}_{+}$, $Pf$ is a measurable function defined by

[TABLE]

We refer the reader to Meyn and Tweedie (2009) for the definitions and basic properties on Markov chains. Given a measurable function $W:\mathsf{Z}\to\left[1,+\infty\right)$, define the $W$-norm of a signed measure $\nu$ on $\mathcal{Z}$ and the $W$-norm of a function $f:\mathsf{Z}\to\mathbb{R}^{d}$:

[TABLE]

these norms generalize resp. the supremum norm of a function and the total variation norm of a measure.

Our results are derived under the following conditions on the kernels:

H 4**.**

1. a)

There exist $\lambda\in\left(0,1\right]$, $b<\infty$ and a measurable function $W:\mathsf{Z}\to\left[1,+\infty\right)$ such that

[TABLE] 2. b)

There exist constants $C<\infty$ and $\rho\in\left(0,1\right)$ such that for any $z\in\mathsf{Z}$ and $n\geq 0$,

[TABLE] 3. c)

There exists a constant $C$ such that for any $\theta,\theta^{\prime}\in\Theta$,

[TABLE]

Sufficient conditions for the uniform-in-$\theta$ ergodic behavior H4b) are given e.g. in  (Fort et al., 2011a, Lemma 2.3.): this lemma shows how to deduce such a control from a minorization condition and a drift inequality on the Markov kernels. Examples of MCMC kernels $P_{\theta}$ satisfying these assumptions can be found in (Andrieu and Moulines, 2006, Proposition 12) and (Saksman and Vihola, 2010, Proposition 15) for the adaptive Hastings-Metropolis algorithm, in (Fort et al., 2011a, Proposition 3.1.) for an interactive tempering sampler, in (Schreck et al., 2013, Proposition 3.2.) for the equi-energy sampler, and in (Fort et al., 2015, Proposition 3.1.) for a Wang-Landau type sampler.

Theorem 3 establishes the convergence of MCPG when the number of points in the Monte Carlo sum $S_{n+1}^{\mathrm{mc}}$ is constant over iterations and the step size sequence $\{\gamma_{n},n\geq 0\}$ vanishes at a convenient rate. It is proved in (Atchadé et al., 2017, Theorem 4).

Theorem 3**.**

Assume H1, H2, H3a-c) and H4a-b). Let $\{\theta_{n},n\geq 0\}$ be the sequence given by Algorithm 1 with a $\left(0,1/L\right]$-valued sequence $\{\gamma_{n},n\geq 0\}$ such that $\sum_{n}\gamma_{n}=+\infty$ and $\sum_{n}\gamma_{n}^{2}<\infty$, and with a constant sequence $\{m_{n},n\geq 0\}$.

In the biased case, assume also H3d) and H4c) and $\sum_{n}|\gamma_{n+1}-\gamma_{n}|<\infty$.

Then, with probability one, there exists $\theta_{\infty}\in\mathcal{L}$ such that $\lim_{n}\theta_{n}=\theta_{\infty}$.

Theorem 4 establishes the convergence of MCPG when the number of points in the Monte Carlo sum $S_{n+1}^{\mathrm{mc}}$ is increasing; it allows a constant stepsize sequence $\{\gamma_{n},n\geq 0\}$. It is proved in (Atchadé et al., 2017, Theorem 6).

Theorem 4**.**

Assume H1, H2, H3a-c) and H4a-b). Let $\{\theta_{n},n\geq 0\}$ be the sequence given by Algorithm 1 with a $\left(0,1/L\right]$-valued sequence $\{\gamma_{n},n\geq 0\}$ and an integer valued sequence $\{m_{n},n\geq 0\}$ such that $\sum_{n}\gamma_{n}=+\infty$ and $\sum_{n}\gamma_{n}^{2}/m_{n}<\infty$.

In the biased case, assume also $\sum_{n}\gamma_{n}/m_{n}<\infty$.

Then, with probability one, there exists $\theta_{\infty}\in\mathcal{L}$ such that $\lim_{n}\theta_{n}=\theta_{\infty}$.

MCPG and SAPG differ in their approximation of $\overline{S}(\theta_{n})$ at each iteration. We provide below a control of this error for a constant or a polynomially increasing batch size $\{m_{n},n\geq 0\}$, and polynomially decreasing stepsize sequences $\{\gamma_{n},n\geq 0\}$ and $\{\delta_{n},n\geq 0\}$.

Proposition 5**.**

Let $\gamma_{\star},\delta_{\star},m_{\star}$ be positive constants and $\beta\in\left[0,1\right)$, $\alpha\geq\beta$, $c\geq 0$. Set $\gamma_{n}=\gamma_{\star}n^{-\alpha}$, $\delta_{n}=\delta_{\star}n^{-\beta}$ and $m_{n}=m_{\star}n^{c}$. Assume H1 to H4. Then

[TABLE]

The proof is given in Appendix C. This proposition shows that when applying MCPG with a constant batch size $(c=0)$, the error $S_{n+1}^{\mathrm{mc}}-\overline{S}(\theta_{n})$ does not vanish; this is not the case for SAPG, since even when $c=0$, the error $S_{n+1}^{\mathrm{sa}}-\overline{S}(\theta_{n})$ vanishes as soon as $\alpha>\beta>0$. Since the case "constant batch size" is the usual choice of the practitioners in order to reduce the computational cost of the algorithm, this proposition supports the use of SAPG instead of MCPG.

We finally study the convergence of SAPG without assuming that the batch size sequence $\{m_{n},n\geq 0\}$ is constant, which implies the following assumption on the sequences $\{\gamma_{n},\delta_{n},m_{n},n\geq 0\}$.

H 5**.**

The step size sequences $\{\gamma_{n},n\geq 0\}$, $\{\delta_{n},n\geq 0\}$ and the batch size sequence $\{m_{n},n\geq 0\}$ satisfy

1. a)

$\gamma_{n}\in\left(0,1/L\right]$, $\delta_{n}\in\left(0,1\right)$, $m_{n}\in\mathbb{N}$, $\sum_{n}\gamma_{n}=+\infty$, $\sum_{n}\gamma_{n}^{2}<\infty$,

[TABLE]

where $\mathsf{D}_{n}:=\sum_{k\geq n}\left(\prod_{j=n}^{k}(1-\delta_{j})\right)$. 2. b)

Furthermore,

[TABLE]

Let us comment this assumption in the case the batch size sequence $\{m_{n},n\geq 0\}$ is constant. This situation corresponds to the "stochastic approximation regime" where the number of draws at each iteration is $m_{n}=1$ (or say, $m_{n}=m$ for any $n$), and it also corresponds to what is usually done by practitioners in order to reduce the computational cost. When $\delta_{n}=\delta_{\star}\in\left(0,1\right)$ for any $n\geq 0$, then $\mathsf{D}_{n}=\delta_{\star}^{-1}$ for any $n\geq 0$. This implies that the condition H5 is satisfied with polynomially decreasing sequences $\gamma_{n}\sim\gamma_{\star}/n^{\alpha}$ with $\alpha\in\left(1/2,1\right]$ (and $m_{n}=m$ for any $n$).

When $\delta_{n}\sim\delta_{\star}\,n^{-\beta}$ for $\beta\in\left(0,1\right)$, then $\mathsf{D}_{n}=O(n^{\beta})$ (see Lemma 9). Hence, using Lemma 9, H5a) and H5b) are satisfied with $\gamma_{n}\sim\gamma_{\star}n^{-\alpha}$ where $\beta<(1+\beta)/2<\alpha\leq 1$, and $m_{n}=m$ for any $n$.

We can not have $\delta_{n}=\delta_{\star}n^{-1}$ since it implies $\mathsf{D}_{n}=+\infty$ for any $n\geq 0$.

Theorem 6**.**

Assume H1, H2, H3 and H4a-b).

Let $\{\theta_{n},n\geq 0\}$ be the sequence given by Algorithm 2 and applied with sequences $\{\gamma_{n},\delta_{n},m_{n},n\geq 0\}$ verifying H5a).

In the biased case, assume also H4c) and H5b).

Then with probability one, there exists $\theta_{\infty}\in\mathcal{L}$ such that $\lim_{n}\theta_{n}=\theta_{\infty}$.

Proof.

The proof is in Section D. ∎

4 Numerical illustration in the convex case

In this section, we illustrate the behavior of the algorithms MCPG and SAPG on a toy example. We first introduce the example and then give some guidelines for a specific choice of the sequences $\{\delta_{n},n\geq 0\}$, $\{\gamma_{n},n\geq 0\}$. Finally, the algorithms are compared more systematically on repeated simulations.

4.1 A toy example

The example is a mixed model, where the regression function is linear in the latent variable $Z$. More precisely, we observe data $(\mathsf{Y}_{1},\cdots,\mathsf{Y}_{N})$ from $N$ subjects, each individual data being a vector of size $J$: $\mathsf{Y}_{k}:=(\mathsf{Y}_{k1},\ldots,\mathsf{Y}_{kJ})$. For the subject $k$, $k=1,\cdots,N$, $\mathsf{Y}_{kj}$ is the $j$-th measurement at time $t_{kj}$, $j=1,\cdots,J$. It is assumed that $\{\mathsf{Y}_{k},k=1,\cdots,N\}$ are independent and for all $k=1,\cdots,N$,

[TABLE]

that is, a linear regression model with individual random intercept and slope, the $\mathbb{R}^{2}$-valued vector being denoted by $Z^{(k)}$. The latent variable is $Z=(Z^{(1)},\ldots,Z^{(N)})$. Furthermore,

[TABLE]

here, $\theta\in\mathbb{R}^{2(D+1)}$ is an unknown parameter and the design matrix $X_{k}\in\mathbb{R}^{2\times 2(D+1)}$ is known

[TABLE]

The optimization problem of the form Eq. (1) that we consider is the log-likelihood $\ell(\theta)$ penalized by a lasso penalty: the objective is the selection of the influential covariates

[TABLE]

on the two components of $Z^{(k)}$. We thus penalize all the elements except $\theta_{1}$ and $\theta_{D+2}$ which correspond to the two intercepts; hence, we set

[TABLE]

The above model is a latent variable model with complete log-likelihood equal to - up to an additive constant

[TABLE]

It is of the form $\phi(\theta)+\left<S(z),\psi(\theta)\right>$ by setting (with $(\cdot)^{\prime}$ denoting the transpose of a matrix)

[TABLE]

The a posteriori distribution $\pi_{\theta}$ is a Gaussian distribution on $\mathbb{R}^{2N}$, equal to the product of $N$ Gaussian distributions on $\mathbb{R}^{2}$:

[TABLE]

Hence, $\overline{S}(\theta)$ is explicit and given by

[TABLE]

with

[TABLE]

Finally, note that in this example, the function $\ell$ is explicit and given by (up to an additive constant)

[TABLE]

Thus $\ell$ is a concave function. Furthermore, in this toy example, $\theta\mapsto\nabla\ell(\theta)$ is linear so that the Lipschitz constant $L$ is explicit and equal to

[TABLE]

where for a matrix A, $\|A\|_{2}$ denotes the spectral norm. Finally, we assumed that $\Theta=\{\theta\in\mathbbm{R}^{2(D+1)}|\|\theta\|<10^{4}\}$ to fulfill the theoretical boundedness assumption. The MCMC algorithm includes a projection step on $\Theta$ if necessary. But in practice, it never happens.

A data set is simulated using this model with $N=40$, $J=8$, $D=300$ and $t_{kj}\in\{0.25,4,6,8,10,12,14,16\}$, $\forall k\in\{1,...,N\}$. The design components $(X_{k1},\ldots,X_{kD})$ (see Eq. (20)) are drawn from a centered Gaussian distribution with covariance matrix $\Gamma$ defined by $\Gamma_{rr^{\prime}}=0.5^{|r-r^{\prime}|}$ ($r,r^{\prime}=1,...,300$). To sample the observations, we use a parameter vector $\theta^{\star}$ defined as follows: $\theta_{1}^{\star}=\theta_{D+2}^{\star}=1$; the other components are set to zero, except $12$ components randomly selected ($6$ among the components $\{2,\cdots,D+1\}$ and $6$ among the components $\{D+3,\cdots,2D+2\}$) and chosen uniformly in $\left[0.5,1.5\right]$ - see the last row on Figure 7.

4.2 Guidelines for the implementation

In this section, we give some guidelines on the choice of the sequences $\{\delta_{n},n\geq 0\}$ and $\{\gamma_{n},n\geq 0\}$. We illustrate the results on single runs of each algorithm. We use the same random draws for all the algorithms to avoid potential differences due to the randomness of the simulations. Similar results have been observed when simulations are replicated. We refer to Section 4.3 for replicated simulations.

Classical sequences $\{\delta_{n},n\geq 0\}$ and $\{\gamma_{n},n\geq 0\}$ are of the form:

[TABLE]

[TABLE]

Impact of $\gamma_{\star}$ and $\delta_{\star}$ on the transient phase: the theoretical study on the asymptotic behavior of SAPG and MCPG is derived under the assumption that $\gamma_{n}\leq 1/L$: when $\alpha>0$, this property holds for any $n$ large enough. In this section, we illustrate the role of $\gamma_{n},\delta_{n}$ for small values of $n$ that is, in the transient phase of the algorithm. In Figure 1, we display the behavior of MCPG and SAPG for two different values of the initial point $\theta_{n=0}$: on the left, it corresponds to a standard initialization ($\theta_{n=0}=(0,\cdots,0)$) while on the right, it corresponds to a poor initialization - which mimics what may happen in practice for challenging numerical applications.

On both plots, we indicate by a vertical line the smallest $n$ such that $\gamma_{n}\leq 1/L$ - remember that in this example, $L$ is explicit (see Eq. (24)). The plots show the estimation of component #245, as a function of the number of iterations $n$. In all cases, $n_{\alpha}=n_{\beta}=0$, $\alpha=0.75$, $m_{n}=60$, and for SAPG, $\beta=0.5$. The dotted blue curve displays a run of SAPG when $(\gamma_{\star},\delta_{\star})=(0.009,0.2)$; the dashed-dotted yellow curve displays a run of SAPG when $(\gamma_{\star},\delta_{\star})=(0.009,0.5)$; the dashed red curve displays a run of SAPG when $(\gamma_{\star},\delta_{\star})=(0.009,0.8)$; the green solid curve displays a run of MCPG when $\gamma_{\star}=0.009$.

The stability of MCPG during the transient phase depends crucially on the first values of the sequence $\{\gamma_{n},n\geq 0\}$. Then when $n$ is large enough so that $\gamma_{n}\leq 1/L$ (after the vertical line), MCPG is more stable and gets smoother. For SAPG, a small value of $\delta_{\star}$ implies an important impact of the initial point $\theta_{n=0}$. When this initial point is poorly chosen, a small value of $\delta_{\star}$ delays the convergence of SAPG. A value of $\delta_{\star}$ around $0.5$ is a good compromise.

Role of ${\alpha}$ and ${\beta}$: Figure 2 displays the behavior of SAPG for different values of $\alpha$ and $\beta$ with $(\gamma_{\star},\delta_{\star})=(0.015,0.5)$, $n_{\alpha}=n_{\beta}=0$ and $m_{n}=60$. The plots show that the larger the parameter $\alpha$ is, the longer the transient phase is. We then recommend to set $\alpha$ close to $0.6$. The parameter $\beta$ seems to have an impact only when $\alpha$ is close to $1$. Therefore, we recommend to set $\delta_{n}$ constant during the transient phase ($n_{\beta}>0$) and then to decrease it rapidly in the convergence phase.

Random stepsize sequence $\{\gamma_{n},n\geq 0\}$: The convergence of the SAPG algorithm can suffer from the scale difference of the parameters, when run with the same stepsize sequence $\{\gamma_{n},n\geq 0\}$ applied to each component of $\theta_{n}$.

Ideally each component of $\theta_{n}$ should have a specific $\gamma_{n}$ value adapted to its scale. But it can be time-consuming to find, by hand-tuning, a sequence that ensures a fast and stable convergence of the algorithm. As an alternative, we suggest to use a matrix-valued random sequence $\{\Gamma^{n},n\geq 0\}$ and replace the update rule of SAPG by

[TABLE]

We propose to define the matrix $\Gamma_{n+1}$ as a diagonal matrix with entries $\Gamma^{n+1}_{ii}$ depending on $H^{n}_{ii}$, where $H^{n}$ is an approximation of the hessian of the likelihood $\ell(\theta)$ (we give an example of such an approximation in Section 5). Through numerical experiments, we observed that asymptotically, $H^{n}$ converges. Hence, to ensure a stepsize sequence decaying like $O(n^{-\alpha})$ asymptotically, we propose the following definition of the random sequence:

[TABLE]

4.3 Long-time behavior of the algorithm

In this section, we illustrate numerically the theoretical results on the long term convergence of the algorithms MCPG, SAPG and SAEM-pen (i.e. Algorithm 3 applied with $S_{n+1}=S_{n+1}^{\mathrm{sa}}$) and EM-pen on the toy model. In this example, the exact algorithm EM-pen (see Eq. (14)) applies: the quantity $\overline{S}(\theta)$ is an explicit expectation under a Gaussian distribution $\pi_{\theta}$. Therefore, we use this example (i) to illustrate the convergence of the three stochastic methods to the same limit point as EM-pen, (ii) to compare the two approximations $S_{n+1}^{\mathrm{mc}}$ and $S_{n+1}^{\mathrm{sa}}$ of $\overline{S}(\theta_{n})$ in a GEM-pen approach, and (iii) to study the effect of relaxing the M-step by comparing the GEM-pen and EM-pen approaches namely SAPG and SAEM-pen.

The sequences $\{\gamma_{n},n\geq 0\}$ and $\{\delta_{n},n\geq 0\}$ are defined as follows: $(\gamma_{\star},\delta_{\star})=(0.004,0.5)$, and $n_{\alpha}=n_{\beta}=0$; three different pairs $(\alpha,\beta)$ are considered: $(\alpha,\beta)=(0.9,0.4)$, $(\alpha,\beta)=(0.6,0.1)$, and $(\alpha,\beta)=(0.5,0.5)$. The algorithms are implemented with a fixed batch size $m_{n}=60$. $100$ independent runs of each algorithm are performed. For the penalty term, we set $\lambda=50$. In MCPG, SAPG and SAEM-pen, the simulation step at iteration $(n+1)$ relies on exact sampling from $\pi_{\theta_{n}}$ - see Eq. (21); therefore, in this toy example, the Monte Carlo approximation of $\overline{S}(\theta_{n})$ is unbiased.

On Figure 3, for the three algorithms MCPG, SAPG and SAEM-pen, the evolution of an approximation of $\|S_{n+1}-\bar{S}(\theta_{n})\|_{2}$ with iterations $n$ is plotted, where, for a random variable $U$, $\|U\|_{2}:=\sqrt{\mathbb{E}\left[\|U\|^{2}\right]}$. This $L_{2}$-norm is approximated by a Monte Carlo sum computed from $100$ independent realizations of $S_{n+1}$; here, $\overline{S}(\theta_{n})$ is explicit (see Eq. (22)). SAEM-pen and SAPG behave similarly; the $L_{2}$-norm converges to [math], and the convergence is slower when $(\alpha,\beta)=(0.6,0.1)$ - this plot illustrates the result stated in Proposition 5, Section 3. This convergence does not hold for MCPG because the size $m_{n}$ of the Monte Carlo approximation is kept fixed.

We compared the limiting vectors $\lim_{n}\theta_{n}$ obtained by each algorithm, over the $100$ independent runs. They are all equal, and the limiting vector is also the limiting value $\theta_{\infty}$ of the EM-pen algorithm. In order to discuss the rate of convergence, we show the behavior of the algorithms when estimating the component $\#245$ of the regression coefficients; this component was chosen among the non-null component of $\theta_{\infty}$. Figure 4 shows the boxplot of $100$ estimations of the component $\#245$ of the vector $\theta_{n}$, when $n=5,25,50,500$, for the algorithms MCPG, SAPG and SAEM-pen with $(\alpha,\beta)=(0.9,0.4)$. Here, SAPG and MCPG behave similarly, with a smaller variability among the $100$ runs than SAEM-pen. SAEM-pen converges faster than SAPG and MCPG which was expected since they correspond respectively to stochastic perturbations of EM-pen and GEM-pen algorithms. Figure 5 shows the boxplot of $100$ estimations by MCPG, SAPG and SAEM-pen of the component $\#245$ after $n=500$ iterations with different values for the parameters $\alpha$ and $\beta$. We observe that the three algorithms give similar final estimates for the three conditions on parameters $\alpha$ and $\beta$. This is due to the fact that with $n_{\alpha}=n_{\beta}=200$, the algorithms have already attained the convergence phase when $n=200$. This allows the algorithms to quickly converge toward the limit points when $n>200$.

Figure 6 shows the convergence of a Monte Carlo approximation of $n\mapsto\mathbb{E}\left[F(\theta_{n})\right]$ based on $100$ independent estimations $\theta_{n}$ obtained by three different algorithms: EM-pen, MCPG, SAPG and SAEM-pen run with $(\alpha,\beta)=(0.9,0.4)$ and $m_{n}=60$. Here again, all the algorithms converge to the same value and EM-pen and SAEM-pen converge faster than MCPG and SAPG. We observe that the path of SAPG is far more smooth than the path of MCPG.

Finally, Figure 7 shows the support of the vector $\lim_{n}\theta_{n}$ (where the component $\theta_{1}$ and $\theta_{302}$ are removed) estimated by MCPG, SAPG, SAEM-pen and EM-pen (the estimated support is the same for the four algorithms). The frequency, among $100$ independent runs, for each component to be in the support of the limit value $\lim_{n}\theta_{n}$, is displayed. Algorithms are implemented with $(\alpha,\beta)=(0.9,0.4)$ and $m_{n}=60$. For all algorithms, we observe that most of the non-null components of $\lim_{n}\theta_{n}$ are non-null components of $\theta^{\star}$. Note also that the stochastic algorithms MCPG, SAPG and SAEM-pen converge to the same vector as EM-pen.

5 Inference in non-linear mixed models for pharmacokinetic data

In this section, SAPG is applied to solve a more challenging problem. The objective is to illustrate the algorithm in cases that are not covered by the theory. The application is in pharmacokinetic analysis, with non-linear mixed effect models (NLMEM); in this application, the penalized maximum-likelihood inference is usually solved by the SAEM-pen algorithm, possibly combined with an approximation of the M-step when it is non explicit. This section also provides a numerical comparison of SAPG and SAEM-pen. Both algorithms have a simulation step; in this more challenging application, it will rely on a Markov chain Monte Carlo (MCMC) sampler - see Section 5.1. Therefore, for both algorithms, $\overline{S}(\theta)$ is approximated by a biased Monte Carlo sum.

We start with a presentation of the statistical analysis and its translation into an optimization problem; we then propose a modification of the SAPG by allowing a random choice of the stepsize sequence $\{\gamma_{n},n\geq 0\}$, to improve the numerical properties of the algorithm. We conclude the section by a comparison of the methods on a pharmacokinetic real data set.

5.1 The non-linear mixed effect model

Pharmacokinetic data are observed along time for $N$ patients. Let $\mathsf{Y}_{k}$ be the vector of the $J$ drug concentrations observed at time $t_{kj}$ ($j\in\{1,\ldots,J\}$) for the $k$-th patient ($k\in\{1,\ldots,N\}$). The kinetic of the drug concentration is described by a non-linear pharmacokinetic regression model $f$, which is a function of time $t$ and unobserved pharmacokinetic parameters $Z^{(k)}$. These parameters are typically the rates of absorption or elimination of the drug by the body. An example is detailed below. The variability among patients is modeled by the randomness of the hidden variables $Z^{(k)}$. These pharmacokinetic parameters may be influenced by covariates, such as age, gender but also genomic variables. Among these high dimension factors, only few of them are correlated to $Z^{(k)}$. Their selection can thus be performed by optimizing the likelihood with a sparsity inducing penalty, an optimization problem that enters problem Eq. (1). However, the likelihood is generally not concave, that is, through this example, we explore beyond the framework in which we are able to prove the convergence of MCPG and SAPG (see Section 3).

Let us now detail the model and the optimization problem. The mixed model is defined as

[TABLE]

where the measurement errors $\epsilon_{kj}$ are centered, independent and identically normally distributed with variance $\sigma^{2}$. Individual parameters $Z^{(k)}$ for the $k$-th subject is a $R$-dimensional random vector, independent of $\epsilon_{kj}$. In a high dimension context, the $Z^{(k)}$’s depend on covariates (typically genomics variables) gathered in a matrix design $X_{k}\in\mathbbm{R}^{R\times(D+1)R}$. The distribution of $Z^{(k)}$ is usually assumed to be normal with independent components

[TABLE]

where $\mu\in\mathbbm{R}^{(D+1)R}$ is the mean parameter vector and $\Omega$ is the covariance matrix of the random parameters $Z^{(k)}$, assumed to be diagonal. The unknown parameters are $\theta=\left(\mu,\Omega_{11},\cdots,\Omega_{RR},\sigma^{2}\right)\in\mathbb{R}^{R(D+1)}\times\left(0,+\infty\right)^{R+1}$.

A typical function $f$ is the two-compartmental pharmacokinetic model with first order absorption, describing the distribution of a drug administered orally. The drug is absorbed from the gut and reaches the blood circulation where it can spread in peripheral tissues. This model corresponds to $f=\frac{A_{c}}{V_{c}}$ with $A_{c}$ defined as

[TABLE]

with $A_{d}(0)=Dose$, $A_{c}(0)=0$, $A_{p}(0)=0$ and where $A_{d}$, $A_{c}$, $A_{p}$ are the amount of drug in the depot, central and peripheral compartments, respectively; $V_{c}$, $V_{p}$ are the volume of the central compartment and the peripheral compartment, respectively; $Q$ and $Cl$ are the inter compartment and global elimination clearances, respectively. To assure positiveness of the parameters, the hidden vector is

[TABLE]

It is easy to show that the model described by Eqs. (28)-(29) belongs to the curved exponential family (see Eq. (4)) with minimal sufficient statistics:

[TABLE]

and $S(z):=\mathrm{Vect}\left(S_{11}(z),\cdots,S_{1N}(z),S_{2}(z),S_{3}(z)\right)$, $\psi:=\mathrm{Vect}\left(\psi_{11},\cdots,\psi_{1N},\psi_{2},\psi_{3}\right)$. The function $\phi$ is given by $\phi(\theta)=-JN\log(\sigma)-\frac{N}{2}\log(|\Omega|)-\frac{1}{2}\sum_{k}(X_{k}\mu)^{\prime}\Omega^{-1}(X_{k}\mu)$.

The selection of genomic variables that influence all coordinates of $Z^{(k)}$ could be obtained by optimizing the log-likelihood penalized by the function $g(\theta)=\lambda\|\mu\|_{1}$, the $L_{1}$ norm of $\mu$ with $\lambda$ a regularization parameter.

However, this estimator is not invariant under a scaling transformation (ie $\tilde{Z}^{(k)}=bZ^{(k)}\mbox{, }\tilde{\mu}=b\mu\mbox{ and }\tilde{\Omega}_{rr}^{1/2}=b\Omega_{rr}^{1/2}$) (see e.g. (Lehmann and Casella, 2006)). In our high dimension experiments, the scale of the hidden variables has a non negligible influence on the selection of the support. To be more precise, let us denote, for $r\in\{1,\ldots,R\}$,

[TABLE]

the coordinates corresponding to the $r$-th pharmacokinetic parameter of function $f$. When the variance $\Omega_{rr}$ of the random parameters $Z_{r}^{(k)}$ is low, the algorithms tend to select too many covariates. This phenomenon is strengthened with a small number of subjects as random effect variances are more difficult to estimate. A solution is to consider the following penalty

[TABLE]

that makes the estimator invariant under scaling transformation. It was initially proposed by Städler et al. (2010) to estimate the regression coefficients and the residual error’s variance in a mixture of penalized regression models. However, the resulting optimization problem is difficult to solve directly because the variance of the random effect $\Omega_{rr}$ appears in the penalty term. Therefore, we propose a new parameterization

[TABLE]

and $\tilde{\theta}:=\{\tilde{\mu},\Sigma_{11},\cdots,\Sigma_{RR},\sigma^{2}\}\in\mathbb{R}^{R(D+1)}\times\left(0,+\infty\right)^{R+1}$. Then, the optimization problem is the following:

[TABLE]

This problem can be solved using MCPG, SAPG or SAEM-pen algorithms. Indeed, the complete log-likelihood is now - up to an additive constant -

[TABLE]

It is again a complete likelihood from the exponential family, with the statistic $S$ unchanged and the functions $\phi$ and $\psi$ given by - up to an additive constant -

[TABLE]

With these definitions of $\phi,\psi$ and $g$, the M-step of SAEM-pen amounts to compute the optimum of a convex function, which is solved numerically by a call to a cyclical coordinate descent implemented in the R package glmnet (Friedman et al., 2010).

MCMC sampler.

In the context of non-linear mixed models, simulation from $\pi_{\theta_{n}}\mathrm{d}\nu$ can not be performed directly like in the toy example. We then use a MCMC sampler based on a Metropolis Hastings algorithm to perform the simulation step. Two proposal kernels are successively used during the iterations of the Metropolis Hastings algorithm. The first kernel corresponds to the prior distribution of $\Sigma Z^{(k)}$ that is the Gaussian distribution $\mathcal{N}(X_{k}\tilde{\mu}_{n},I)$. The second kernel corresponds to a succession of $R$ uni-dimensional random walk in order to update successively each component of $Z^{(k)}$. The variance of each random walk is automatically tuned to reach a target acceptance ratio following the principle of an adaptive MCMC algorithm (Andrieu and Thoms, 2008).

Adaptive random stepsize sequences.

In the context of NLMEM, numerical experiments reveal that choosing a deterministic sequence $\{\gamma_{n},n\geq 0\}$ that achieve a fast convergence of SAPG algorithm could be difficult. Indeed, parameters to estimate are of different scales. For example, random effect and residual variances are constrained to be positive. Some of them are close to zero, some are not. As explained in Section 4.2, an alternative is to implement a matrix-valued random sequence $\{\Gamma^{n},n\geq 0\}$. The gradient and the hessian of the likelihood $\ell(\theta)$ can be approximated by stochastic approximation using the Louis principle (see McLachlan and Krishnan, 2008, Chapter 4). Let us denote $H_{n}$ the stochastic approximation of the hessian obtained at iteration $n$ as explained by Samson et al. (2007). Note that no supplementary random samples are required to obtain this approximation. Along the iterations, each diagonal entry of the matrix ${H}^{n}$ converges: this limiting value can be seen as a simple way to automatically tune a good $\gamma_{\star}$, that is parameter specific. The entries $\Gamma^{n+1}_{ii}$ are then defined by Eq. (27).

5.2 Simulated data set.

The convergence of the corresponding algorithms is illustrated on simulated data. Data are generated with the model defined by Eq. (5.1) and $N=40$, $J=12$, $D=300$. The design matrix $X_{k}$ is defined by Eq. (20), with components $(X_{k1},\ldots,X_{kD})$ drawn from $\mathcal{N}(0,\Gamma)$ with $\Gamma_{ii^{\prime}}=0.5^{|i-i^{\prime}|}$ ($i,i^{\prime}=1,...,300$). Parameter values are

[TABLE]

the other components are set to zero, except $\mu_{4}$ and $\mu_{912}$ that are set to $1$. The matrix $\Omega$ is diagonal with diagonal elements equal to $(0.16,0.16,0.16,0.04,0.04)$.

The penalty function is set to

[TABLE]

only the parameters corresponding to a covariate effect being penalized. The optimization problem Eq. (1) with regularization parameter $\lambda=190$ is solved on this dataset with SAEM-pen and SAPG; we run SAPG with the random sequence $\{\Gamma_{n},n\geq 0\}$ as described above (see (27)) with $n_{0}=9500$. For both algorithms, the stochastic approximation step size was set to:

[TABLE]

We set $\alpha=0.75$ and $\beta=0.499$. Figure 8 shows the convergence of SAEM-pen and three parameterizations of SAPG: i) a version with $\gamma^{\star}=0.005$ for all the components of $\theta$, ii) a version with $\gamma^{\star}=0.005$ for $\tilde{\mu}$, $\gamma^{\star}=0.0005$ for $\Sigma$ and $\gamma^{\star}=0.03$ for $\sigma$, and iii) a version with adaptive random step sizes. For the four algorithms, all the parameters corresponding to a covariate effect are estimated to zero except the two components $\mu_{4}$ and $\mu_{912}$. The version of SAPG with a same $\gamma^{\star}$ for all the component is the one that converge the most slowly. When the $\gamma^{\star}$ is tuned differently according the type of parameters, the convergence of SAPG is accelerated. Algorithms SAEM-pen and SAPG with adaptive random step sizes have a similar fast convergence profile.

Figure 9 presents the evolution of four entries of the matrix $\Gamma^{n}$ along the iterations of SAPG, corresponding to the components $\tilde{\mu}_{904}$, $\tilde{\mu}_{912}$, $\Sigma_{44}$ and $\sigma$. We can notice that they are not on the same scale. They vary during the first iterations and converge to limiting values before iteration $n_{0}=9500$. Then the step sizes decrease to 0, following the definition given in Eq. (27).

5.3 Application to real data

Algorithms SAEM-pen and SAPG with matrix-valued random sequence $\{\Gamma^{n},n\geq 0\}$ are applied to real data of the pharmacokinetic of dabigatran ($DE$) from two cross over clinical trials (Delavenne et al., 2013; Ollier et al., 2015). These $2$ trials studied the drug-drug interaction between $DE$ and different Pgp-inhibitors. From these $2$ trials, the pharmacokinetics of $DE$ are extracted from $15$ subjects with no concomitant treatment with Pgp-inhibitors. The concentration of dabigatran is measured at $9$ sampling times for each patient. Each subject is genotyped using the DMET$\textregistered$ microarray from Affymetrix. Single nucleotide polymorphisms (SNP) showing no variability between subjects are removed and $264$ SNP are included in the analysis.

Function $f$ of the non-linear mixed model is defined as the two compartment pharmacokinetic model with first order absorption previously described (see Eq. (5.1)) (Delavenne et al., 2013). The penalty function $g$ is defined by Eq. (32).

Because of the limited number of subjects, the influence of genetic covariates is only studied on $V_{c}$ and $Cl$ parameters, that characterize the elimination process and are the most likely to be influenced by the genetic. Finally, random effect variances of $Q$ and $V_{p}$ are set to $0.01$ in accordance with previously published population pharmacokinetic of dabigatran (Delavenne et al., 2013). The other variance parameters are estimated. The penalized likelihood problem (Eq. 31) is solved on the data with the SAEM-pen and SAPG algorithms, for $40$ different values of parameter $\lambda$. SAPG algorithm is run using the random sequence $\{\Gamma^{n},n\geq 0\}$ given in Eq. (27). The best regularization parameter $\lambda$ is chosen with a data-driven approach based on the EBIC criteria (Chen and Chen, 2008).

Figure 10 shows the results. The regularization paths of $Cl$ and $V_{c}$ parameters using both algorithms correspond to the evolution of covariate coefficient estimates as a function of the value of $\lambda$. They are reconstructed with low noise for both algorithms, are very similar for high values of $\lambda$ but less for lower values of $\lambda$.

Finally, the selected model has all covariates parameters set to zero. This means that none of the genetic covariates influence the distribution of the individual parameters. This result is not surprising given the low number of subjects and the fact that a large part of the inter individual variability is due to the dissolution process of the drug (Ollier et al., 2015) and is therefore not influenced by genetic covariates. This lack of relationship between dabigtran’s pharmacokinetic parameters and genetic covariates has already been highlighted in an other study (Gouin-Thibault et al., 2017).

6 Conclusion

In this work, we propose a new Stochastic Proximal-Gradient algorithm to solve penalized maximum likelihood problems when the likelihood is intractable: the the gradient is approximated through a stochastic approximation scheme. We provide a theoretical convergence analysis of this new algorithm and illustrate these results numerically on a simulated toy example in the case of a concave likelihood function. The robustness to the non concave case is explored through a more challenging application to population pharmacokinetic analysis relying on penalized inference in non-linear mixed effects models.

Appendix

Appendix A Proof of Proposition 1

Lemma 7**.**

Under the assumptions of Proposition 1, for any $\gamma\in\left(0,1/L\right]$, $s\in\overline{S}(\Theta)$ and any $\theta,\theta^{\prime}\in\Theta$,

[TABLE]

Proof.

Fix $\theta^{\prime}\in\Theta$ and $s\in\overline{S}(\Theta)$. The derivative of the function $\theta\mapsto\mathcal{L}(\theta):=\phi(\theta)+\left<s,\psi(\theta)\right>$ is $\nabla\phi(\theta)+\operatorname{J}\psi(\theta)s$ and this gradient is $L$-lipschitz. From a Taylor expansion to order $1$ at $\theta^{\prime}$ and since the gradient is Lipschitz, we have

[TABLE]

We then choose $s=\overline{S}(\theta^{\prime})$, use $L\leq 1/\gamma$ and conclude by the equality $2\left<a,b\right>-\|a\|^{2}=\|b\|^{2}-\|a-b\|^{2}$. ∎

Proof of Proposition 1

We prove that $\mathcal{Q}(\theta_{n+1}|\theta_{n})-g(\theta_{n+1})\geq\mathcal{Q}(\theta_{n}|\theta_{n})-g(\theta_{n})$ so that the sequence $\{\theta_{n},n\geq 0\}$ defined by Eq. (10) is a sequence satisfying Eq. (15).

By Lemma 7, it holds for any $\theta\in\Theta$ and any $\gamma\in\left(0,1/L\right]$

[TABLE]

Note that the RHS and the LHS are equal when $\theta=\theta_{n}$ so that for any point $\tau$ which maximizes the RHS, it holds $\mathcal{Q}(\tau|\theta_{n})-g(\tau)\geq\mathcal{Q}(\theta_{n}|\theta_{n})-g(\theta_{n})$. This concludes the proof upon noting that such a point $\tau$ is unique and equal to $\theta_{n+1}$ given by Eq. (10).

Appendix B Technical lemmas

Define

[TABLE]

Lemma 8**.**

For any $n\geq 2$, $\sum_{j=2}^{n}\Delta_{j+1:n}\,\delta_{j}=1-\Delta_{2:n}$.

Proof.

For any $j\leq n$, we have $\Delta_{j+1:n}-\Delta_{j:n}=\delta_{j}\Delta_{j+1:n}$ from which the result follows. ∎

Lemma 9**.**

Let $\beta\in\left(0,1\right)$ and $\delta_{\star}>0$. Set $\delta_{n}=\delta_{\star}n^{-\beta}$ for any $n\geq 1$. Then for any $k$ large enough,

[TABLE]

Furthermore, $\left|\delta_{n+1}\mathsf{D}_{n+2}-\delta_{n}\mathsf{D}_{n+1}\right|=O(1/n^{1+(1-\beta)\wedge\beta})$.

The proof of Lemma 9 relies on standard Taylor’s expansions with explicit formulas for the remainder. The proof is omitted.

Lemma 10**.**

Let $\beta\in\left[0,1\right)$ and $\delta_{\star}>0$. For any $r$, when $n\to\infty$,

[TABLE]

Proof.

We have

[TABLE]

Let $q_{\star}\geq 0$ such that for any $q\geq q_{\star}$, $q(1-\beta)+1-r>0$. For any constant $D>0$, there exist constants $C,C^{\prime}$ (whose value can change upon each appearance) such that

[TABLE]

This concludes the proof. ∎

Lemma 11**.**

Let $\{\mathsf{A}_{n},n\geq 0\}$ be a sequence of $d^{\prime}\times q$ matrices and $\{\sigma_{n},n\geq 0\}$ be a sequence of $q\times 1$ vectors. Let $\{S_{n}^{\mathrm{sa}},n\geq 0\}$ be given by Eq. (12). For any $n\geq 2$

[TABLE]

Proof.

By definition of $S_{n}^{\mathrm{sa}}$, it holds $\mathsf{A}_{n}\left(S_{n}^{\mathrm{sa}}-\sigma_{n-1}\right)=(1-\delta_{n})\mathsf{A}_{n-1}\left(S_{n-1}^{\mathrm{sa}}-\sigma_{n-2}\right)+B_{n}$ where

[TABLE]

By iterating, we have

[TABLE]

from which the lemma follows. ∎

Lemma 12**.**

Assume H4a). Let $\{S_{n}^{\mathrm{sa}},n\geq 0\}$ be given by Eq. (12). Then

[TABLE]

Proof.

By H4a), there exists a constant $C<\infty$ such that for any $n\geq 1$ and $1\leq j\leq m_{n}$, $\|S(Z_{j,n-1})\|^{2}\leq C\,W(Z_{j,n-1})$. In addition, by the drift assumption on the kernels $P_{\theta}$, we have

[TABLE]

Similarly, by using $Z_{0,n-1}=Z_{m_{n-1},n-2}$, we have

[TABLE]

A trivial induction shows that

[TABLE]

from which the first two results follow. For the third one: by Lemma 11 applied with $\mathsf{A}_{n}=I$ (the identity matrix) and $\sigma_{n}=0$, we have for any $n\geq 1$,

[TABLE]

By the Minkowsky inequality and the inequality $(a+b)^{2}\leq 2a^{2}+2b^{2}$, we have

[TABLE]

By definition, $\Delta_{2:n}\in\left[0,1\right]$ and by Lemma 8,

[TABLE]

Hence, $\sup_{n\geq 0}\mathbb{E}\left[\|S_{n}^{\mathrm{sa}}\|^{2}\right]<\infty$. ∎

Define the Proximal-Gradient operator

[TABLE]

Lemma 13**.**

Assume H1, H2 and H3. Let $\{S_{n}^{\mathrm{sa}},n\geq 0\}$ be given by (12). Then, for the sequence $\{\theta_{n},n\geq 0\}$ given by Algorithm 2,

1. (i)

There exists a constant $C$ such that almost-surely, for any $n\geq 0$,

[TABLE] 2. (ii)

There exists a constant $C^{\prime}$ such that almost-surely, for any $n\geq 0$,

[TABLE] 3. (iii)

There exists a constant $C"$ such that almost-surely, for any $n\geq 0$,

[TABLE]

Proof.

The proof of (i) is on the same lines as the proof of (Atchadé et al., 2017, Lemma 15), and is omitted. For (ii), we write by using H3b) and H3c),

[TABLE]

We then conclude by (i). The LHS in (iii) is upper bounded by

[TABLE]

Under H3, there exists a constant $C$ such that for all $\gamma,\gamma^{\prime}\in\left(0,1/L\right]$ and $\theta,\theta^{\prime}\in\Theta$ (see (Atchadé et al., 2017, Proposition 12))

[TABLE]

We then conclude by (i). ∎

Lemma 14**.**

Assume H4. For any $\theta\in\Theta$, there exists a function $\widehat{S}_{\theta}:\mathsf{Z}\to\mathbb{R}^{q}$ such that $S-\overline{S}(\theta)=\widehat{S}_{\theta}-P_{\theta}\widehat{S}_{\theta}$ and $\sup_{\theta\in\Theta}\left|\widehat{S}_{\theta}\right|_{\sqrt{W}}<\infty$. In addition, there exists a constant $C$ such that for any $\theta,\theta^{\prime}\in\Theta$,

[TABLE]

Proof.

Set $\widehat{S}_{\theta}(z):=\sum_{n\geq 0}\left(P_{\theta}^{n}S(z)-\overline{S}(\theta)\right)$. Observe that, when exists, this function satisfies $S-\overline{S}(\theta)=\widehat{S}_{\theta}-P_{\theta}\widehat{S}_{\theta}$. Note that under H4a)-H4b), there exist $C$ and $\rho\in\left(0,1\right)$ such that for any $\theta\in\Theta$,

[TABLE]

the RHS is finite, thus showing that $\widehat{S}_{\theta}$ exists. This inequality also proves that $\sup_{\theta}\left|\widehat{S}_{\theta}\right|_{\sqrt{W}}<\infty$. The Lipschitz property is established in  (Fort et al., 2011a, Lemma 4.2.) and its proof uses H4c). ∎

Appendix C Proof of Proposition 5

Throughout this section, set $\|U\|_{L_{2}}:=\mathbb{E}\left[\|U\|^{2}\right]^{1/2}$. By Lemma 11, $\|S_{n}^{\mathrm{sa}}-\overline{S}(\theta_{n-1})\|_{L_{2}}\leq\sum_{i=1}^{3}\mathcal{T}_{i,n}$ with

[TABLE]

Since $\Delta_{2:n}\leq\exp(-\delta_{\star}\sum_{j=2}^{n}j^{-\beta})$, then

[TABLE]

By H3c), Lemma 12 and Lemma 13, there exists a constant $C$ such that $\mathcal{T}_{2,n}\leq C\,\sum_{j=2}^{n}\Delta_{j:n}\gamma_{j-1}$. By Lemma 10, this yields $\mathcal{T}_{2,n}=O(n^{\beta-\alpha})$. For the last term, we use a martingale decomposition.

By Lemma 14, there exists a function $\widehat{S}_{\theta}$ such that

[TABLE]

and $\sup_{\theta\in\Theta}\left|\widehat{S}_{\theta}\right|_{\sqrt{W}}<\infty$. Hence, we write

[TABLE]

with

[TABLE]

we used that $Z_{0,j}=Z_{m_{j},j-1}$. Upon noting that $\partial M_{j}$ is a martingale-increment, and

[TABLE]

is a martingale-increment, we have by two successive applications of (Hall and Heyde, 1980, Theorem 2.10):

[TABLE]

By Lemma 10, this term is $O(n^{-(\beta+c)/2})$. For the second term, we write

[TABLE]

By Lemma 12 and Lemma 14, the RHS is $O(n^{-(\beta+c)}+n^{-(1+c)})$ so that this second term is $O(n^{-(\beta+c)})$. Finally, for the third term, by using Lemma 12, Lemma 13 and Lemma 14, we write

[TABLE]

Again by Lemma 10, this last term is $O(n^{-(\alpha+c)})$. Therefore, $\mathcal{T}_{3,n}=O(n^{-(\beta+c)/2})$.

Appendix D Proof of Theorem 6

Throughout the proof, we will write $S_{n+1}$ instead of $S_{n+1}^{\mathrm{sa}}$.

Proof of Theorem 6

We prove the almost-sure convergence of the three random sums given in Theorem 2. The third one is finite almost-surely since its expectation is finite (see Proposition 15). The first two ones are of the form $\sum_{n}\mathsf{A}_{n+1}\left(S_{n+1}-\overline{S}(\theta_{n})\right)$ where $\mathsf{A}_{n+1}$ is respectively

[TABLE]

Note that $\mathsf{A}_{n+1}\in\mathcal{F}_{n}$ (the filtration is defined by Eq. (17)). By Lemma 13 and H3b-c), for both cases, there exists a constant $C$ such that almost-surely, for any $n\geq 0$,

[TABLE]

We then conclude by Proposition 16.

Proposition 15**.**

Assume H4a) and

[TABLE]

Then there exists a constant $C$ such that

[TABLE]

Proof.

We write

[TABLE]

and conclude by Lemma 12. ∎

Proposition 16**.**

Let $\{\theta_{n},n\geq 0\}$ be given by Algorithm 2. Assume H1, H3, H4a-b) and H5a). In the biased case, assume also H4c) and H5b). Let $\{\mathsf{A}_{n},n\geq 0\}$ be a sequence of $d^{\prime}\times q$ random matrices such that for any $n\geq 0$, $\mathsf{A}_{n+1}\in\mathcal{F}_{n}$, and there exists a constant $C_{\star}$ such that almost-surely

[TABLE]

here $a_{n+1}:=\gamma_{n}\gamma_{n+1}+\gamma_{n}^{2}+|\gamma_{n+1}-\gamma_{n}|$. Then, almost-surely, the series $\sum_{n}\mathsf{A}_{n+1}\left(S_{n+1}-\overline{S}(\theta_{n})\right)$ converges.

By Lemma 11 applied with $\sigma_{n}=\overline{S}(\theta_{n})$, we decompose this sum into four terms:

[TABLE]

We have by using Eq. (35),

[TABLE]

By H5a), $\mathsf{D}_{2}<\infty$ so the RHS is finite thus implying that $\mathcal{T}_{1}$ is finite almost-surely.

Using Eq.(36), there exists a constant $C$ such that

[TABLE]

By H3b)-H3c), H5a) and Lemma 12, the RHS is finite thus implying that $\mathcal{T}_{2}$ is finite almost-surely.

Similarly, there exists a constant $C$ such that

[TABLE]

By H3c), the RHS is bounded (up to a multiplicative constant) by $\sum_{j\geq 2}\gamma_{j}\mathsf{D}_{j}\mathbb{E}\left[\|\theta_{j-1}-\theta_{j-2}\|\right]$; and by H5a) and Lemmas 12 and 13, this sum is finite. Hence $\mathcal{T}_{3}$ is finite almost-surely.

We give the proof of the convergence of the last term in the biased case: $\mathbb{E}\left[S(Z_{k,n})|\mathcal{F}_{n}\right]\neq\overline{S}(\theta_{n})$. The proof in the unbiased case corresponds to the following lines with $R_{j,1}=R_{j,2}=0$ and $\widehat{S}_{\theta}=S$. Set $\overline{\mathsf{D}}_{j}:=\delta_{j}(1+\mathsf{D}_{j+1})$. By Lemma 14, there exists $\widehat{S}_{\theta}$ such that

[TABLE]

and $\sup_{\theta\in\Theta}\left|\widehat{S}_{\theta}\right|_{\sqrt{W}}<\infty$. Hence, we have

[TABLE]

where

[TABLE]

Upon noting that $\mathbb{E}\left[\mathsf{A}_{j}\partial M_{j}|\mathcal{F}_{j-1}\right]=0$, the almost-sure convergence of the series $\sum_{j}\overline{\mathsf{D}}_{j}\mathsf{A}_{j}\partial M_{j}$ is proved by checking criteria for the almost-sure convergence of a martingale. By (35), there exists a constant $C$ such that

[TABLE]

By H5a), Lemma 12 and (Hall and Heyde, 1980, Theorem 2.10), the RHS is finite. (Hall and Heyde, 1980, Theorem 2.17) implies that $\sum_{j}\overline{\mathsf{D}}_{j}\mathsf{A}_{j}\,\partial M_{j}$ is finite almost-surely. For the second term, we write

[TABLE]

so that, by Lemmas 12 and 14, this series is finite almost-surely if $\sum_{j}\mathbb{E}\left[\|m_{j+1}^{-1}\overline{\mathsf{D}}_{j+1}\mathsf{A}_{j+1}-m_{j}^{-1}\overline{\mathsf{D}}_{j}\mathsf{A}_{j}\|\right]<\infty$. From Eq. (35) and Eq. (36), there exists a constant $C$ such that

[TABLE]

H5 and Lemma 12 imply that

[TABLE]

Finally, by (35), Lemmas 12 to 14, there exists a constant $C$ such that

[TABLE]

The RHS is finite by H5 thus implying that $\sum_{j}\overline{\mathsf{D}}_{j}\mathsf{A}_{j}R_{j,2}$ is finite almost-surely.