On Perturbed Proximal Gradient Algorithms
Yves F. Atchadé University of Michigan, 1085 South University, Ann Arbor, 48109, MI,
United States,

yvesa@umich.edu

arXiv:1402.2365v4 [] 19 Nov 2016

Gersende Fort LTCI, CNRS, Telecom ParisTech, Université Paris-Saclay. 46, rue Barrault 75013 Paris, France,
gersende.fort@telecom-paristech.fr
Eric Moulines LTCI, CNRS, Telecom ParisTech, Université Paris-Saclay. 46, rue Barrault 75013 Paris, France,

eric.moulines@polytechnique.edu

Editors: Léon Bottou

Abstract
We study a version of the proximal gradient algorithm for which the gradient is
intractable and is approximated by Monte Carlo methods (and in particular Markov
Chain Monte Carlo). We derive conditions on the step size and the Monte Carlo
batch size under which convergence is guaranteed: both increasing batch size and
constant batch size are considered. We also derive non-asymptotic bounds for an
averaged version. Our results cover both the cases of biased and unbiased Monte
Carlo approximation. To support our findings, we discuss the inference of a sparse
generalized linear model with random effect and the problem of learning the edge
structure and parameters of sparse undirected graphical models.
Keywords: Proximal Gradient Methods; Stochastic Optimization; Monte Carlo
approximations; Perturbed Majorization-Minimization algorithms.

1. Introduction
This paper deals with statistical optimization problems of the form:
(P)

min F (θ)

with F = f + g .

θ∈Rd

This problem occurs in a variety of statistical and machine learning problems, where f
is a measure of fit depending implicitly on some observed data and g is a regularization
term that imposes structure to the solution. Typically, f is a differentiable function
with a Lipschitz gradient, whereas g might be non-smooth (typical examples include
sparsity inducing penalty).
H 1 The function g : Rd → [0, +∞] is convex, not identically +∞, and lower semicontinuous. The function f : Rd → R is convex, continuously differentiable on Rd
and there exists a finite non-negative constant L such that, for all θ, θ0 ∈ Rd ,
k∇f (θ) − ∇f (θ0 )k ≤ Lkθ − θ0 k ,
where ∇f denotes the gradient of f .
def

We denote by Θ the domain of g: Θ = {θ ∈ Rd : g(θ) < ∞}.
1

H 2 The set argminθ∈Θ F (θ) is a non empty subset of Θ.
In this paper, we focus on the case where f + g and ∇f are both intractable. This
setting has not been widely considered despite the considerable importance of such
models in statistics and machine learning. Intractable likelihood problems naturally
occur for example in inference for bayesian networks (e.g. learning the edge structure and the parameters in an undirected graphical models), regression with latent
variables or random effets, missing data, etc... In such applications, f is the negated
log-likelihood of a conditional Gibbs measure πθ known only up to a normalization
constant and the gradient of ∇f (θ) is typically expressedRas a very high-dimensional
integral w.r.t. the associated Gibbs measure ∇f (θ) = Hθ (x)πθ (dx). Of course,
this integral cannot be computed in closed form and should be approximated. Most
often, some forms of Monte Carlo integration (such as Markov Chain Monte Carlo,
or MCMC) is the only option.
To cope with problems where f +g is intractable and possibly non-smooth, various
methods have been proposed. Some of these works focused on stochastic sub-gradient
and mirror descent algorithms; see Nemirovski et al. (2008); Duchi et al. (2011); Cotter
et al. (2011); Lan (2012); Juditsky and Nemirovski (2012a,b). Other authors have
proposed algorithms based on proximal operators to better exploit the smoothness
of f and the properties of g (see e.g. Combettes and Wajs (2005); Hu et al. (2009);
Xiao (2010); Juditsky and Nemirovski (2012a,b)).
The current paper focuses on the proximal gradient algorithm (see e.g. Beck and
Teboulle (2010); Combettes and Pesquet (2011); Parikh and Boyd (2013) for literature
review and further references). The proximal map (Moreau (1962)) associated to g is
defined for γ > 0 and θ ∈ Rd by:


1
def
2
Proxγ,g (θ) = argminϑ∈Θ g(ϑ) +
kϑ − θk .
(1)
2γ
Note that under H1, there exists an unique point ϑ minimizing the RHS of (1) for
any θ ∈ Rd and γ > 0. The proximal gradient algorithm is an iterative algorithm
which, given an initial value θ0 ∈ Θ and a sequence of positive step sizes {γn , n ∈ N},
produces a sequence of parameters {θn , n ∈ N} as follows:
Algorithm 1 (Proximal gradient algorithm) Given θn , compute
θn+1 = Proxγn+1 ,g (θn − γn+1 ∇f (θn )) .

(2)

When γn = γ for any n, it is known that the iterates of the proximal gradient
algorithm {θn , n ∈ N} (Algorithm 1) converges to θ∞ , this point is a fixed point of
the proximal-gradient map
def

Tγ (θ) = Proxγ,g (θ − γ∇f (θ)) .

(3)

Under H1 and H2, when γn (0, 2/L] and inf n γn > 0, it is indeed known that the
iterates of the proximal gradient algorithm {θn , n ∈ N} defined in (2) converges to a
2

point in the set L of the solutions of (P ) which coincides with the fixed points of the
mapping Tγ for any γ ∈ (0, 2/L)
def

L = argminθ∈Θ F (θ) = {θ ∈ Θ : θ = Tγ (θ)} .

(4)

(see e.g. (Combettes and Wajs, 2005, Theorem 3.4. and Proposition 3.1.(iii))).
Since ∇f (θ) is intractable, the gradient ∇f (θn ) at n-th iteration is replaced by
an approximation Hn+1 :
Algorithm 2 (Perturbed Proximal Gradient algorithm) Let θ0 ∈ Θ be the initial solution and {γn , n ∈ N} be a sequence of positive step sizes. For n ≥ 1, given
(θ0 , . . . , θn ) construct an approximation Hn+1 of ∇f (θn ) and compute
θn+1 = Proxγn+1 ,g (θn − γn+1 Hn+1 ) .

(5)

We provide in Theorem 2 sufficient conditions on the perturbation ηn+1 = Hn+1 −
∇f (θn ) to obtain the convergence of the perturbed proximal gradient sequence given
by (5). We then consider an averaging scheme of the perturbed proximal gradient algorithm: given non-negative weights
Theorem 3 provides non-asymptotic
P {an , n ∈ N},P
bound of the deviation between nk=1 ak F (θk )/ nk=1 ak and the minimum of F . Our
results complement and extend Rosasco et al. (2014); Nitanda (2014); Xiao and Zhang
(2014).
R
We then consider the case where the gradient ∇f (θ) = X Hθ (x)πθ (dx) is defined as an expectation (see H3 in section 3). In this case, at each iteration ∇f (θn )
Pmn+1
(j)
is approximated by a Monte Carlo average Hn+1 = m−1
n+1
j=1 Hθn (Xn+1 ) where
(j)

mn+1 is the size of the Monte Carlo batch and {Xn+1 , 1 ≤ j ≤ mn+1 } is the Monte
Carlo batch. Two different settings are covered. In the first setting, the samples
(j)
{Xn+1 , 1 ≤ j ≤ mn+1 } are conditionally independent and identically distributed
(i.i.d.) with distribution πθn . In such case, the conditional expectation of Hn+1 given
all the past iterations, denoted by E [ Hn+1 | Fn ] (see section 3), is equal to ∇f (θn ).
(j)
In the second setting, the Monte Carlo batch {Xn+1 , 1 ≤ j ≤ mn+1 } is produced by
(j)

running a MCMC algorithm. In such case, the conditional distribution of Xn+1 given
the past is no longer exactly equal to πθn which implies that E [ Hn+1 | Fn ] 6= ∇f (θn ).
Theorem 4 (resp. Theorem 6) establish the convergence of the sequence {θn , n ∈
N} when the batch size mn is either fixed or increases with the number of iterations n.
(j)
When the Monte Carlo batch {Xn+1 , 1 ≤ j ≤ mn+1 } is i.i.d. conditionally to the past
the two theorems essentially say that with probability
one, {θn , nP
∈ N} converges to an
P
2 /m
element of the set of minimizer L as soon as n γn = +∞ and n γn+1
n+1 < ∞.
Hence, one can choose either a fixed step size γn = γ and a batch size {mn , n ∈ N}
increasing at least linearly (up to a logarithmic factor); or a decreasing step size and
(j)
a fixed batch size mn = m. When {Xn+1 , 1 ≤ j ≤ mn+1 } is produced by a MCMC
algorithm (under appropriate
essentially say that the same
Passumptions) our
P theorems
2
convergence result holds
if
γ
=
∞
and
γ
<
∞
when mn = m is constant
n
n
n n+1
P
across iterations or n γn+1 /mn+1 < ∞ if the batch size is increased.
Theorem
4 and Theorem
6 also provide non asymptotic bounds for the difference
P
P
∆n = nk=1 ak F (θk )/ nk=1 ak − min F in Lq -norm for q ≥ 1. When the batch size
3

sequence mn+1 increases linearly at each iteration while the step size γn+1 is held
constant, ∆n = O(ln n/n). We recover (up to a logarithmic factor) the rate of the
proximal gradient algorithm. If we now compare the complexity of the algorithms
in terms of the number of simulations N needed (and not the number of iterations),
the error bound decreases like O(N −1/2 ). The same error bound can be achieved by
√
choosing a fixed batch size and a decreasing step size γn = O(1/ n).
In section 4, these results are illustrated with the problem of estimating a highdimensional discrete graphical models. In section 5, we consider high-dimensional
random effect logistic regression model. All the proofs are postponed to section 6.

2. Perturbed proximal gradient algorithms
The key property to study the behavior of the sequence the perturbed proximal
gradient algorithm is the following elementary lemma which might be seen as a deterministic version of the Robbins-Siegmund lemma (see e.g. (Polyak, 1987, Lemma 11,
Chapter 2)). It replaces in our analysis (Combettes, 2001, Lemma 3.1) for quasi-Fejer
sequences and modified Fejer monotone sequences (see Lin et al. (2015)). Compared
to the Robbins-Siegmund Lemma, the sequence (ξn )n is not assumed to be nonnegative. When applied in the stochastic context as in Section 3, the fact that the result
is purely deterministic and deals with signed perturbations ξn allows more flexibility
in the study of the dynamics.
Lemma 1 Let {v
Pn , n ∈ N} and {χn , n ∈ N} be non-negative sequences and {ξn , n ∈
N} be such that n ξn exists. If for any n ≥ 0,
vn+1 ≤ vn − χn+1 + ξn+1
then

P

n χn

< ∞ and limn vn exists.

Proof See Section 6.2.1
Applied with vn = kθn − θ? k for some θ? ∈ L, this lemma is the key result for the
proof of the following theorem, which provides sufficient conditions on the stepsize
sequence {γn , n ∈ N} and on the approximation error :
def

ηn+1 = Hn+1 − ∇f (θn ) ,

(6)

for the sequence {θn , n ∈ N} to converge to a point θ∞ in the set L of the minimizers
of F . Denote by h·, ·i the usual inner product on Rd associated to the norm k · k.
Theorem 2 Assume H1 and H2. Let {θn , n ∈ N}
P be given by Algorithm 2 with step
sizes satisfying γn ∈ (0, 1/L] for any n ≥ 1 and n γn = +∞. If the following series
converge
X
X
X
2
γn+1 Tγn+1 (θn ), ηn+1 ,
γn+1 ηn+1 ,
γn+1
kηn+1 k2 ,
(7)
n≥0

n≥0

then there exists θ∞ ∈ L such that limn θn = θ∞ .
4

n≥0

Proof See Section 6.2.2.
Theorem 2 applied with ηn+1 = 0 provides sufficient conditions for the convergence
of
P
Algorithm 1 to L: the algorithm converges as soon as γn ∈ (0, 1/L] and n γn = +∞.
Sufficient conditions for the convergence of {θn , n ∈ N} are also provided in Combettes and Wajs (2005). P
When applied to our settings, Theorem 3.4. in Combettes
and Wajs (2005) require n kηn+1 k < ∞ and inf n γn > 0, which for instance cannot
accommodate the fixed Monte Carlo batch size stochastic algorithms considered in
this paper. The same limitation applies to the analysis of the stochastic quasi-Fejer
iterations
(see Combettes and Pesquet (2015a)) which in our particular case requires
P
n γn+1 kηn+1 k < ∞. These conditions are weakened in Theorem 2. However in all
fairness we should mention that unlike the present work, Combettes and Wajs (2005)
and Combettes and Pesquet (2015a) deal with infinite-dimensional problems which
raises additional technical difficulties, and study algorithms that include a relaxation
parameter. Furthermore, in the case where ηn ≡ 0, larger values of the stepsize γn
are allowed (γn ∈ (0, 2/L].)
Let {a0 , · · · , anP
} be non-negative real numbers. Theorem 3 provides a control of
the weighted sum nk=1 ak (F (θk ) − min F ).
Theorem 3 Assume H1 and H2. Let {θn , n ∈ N} be given by Algorithm 2 with
γn ∈ (0, 1/L] for any n ≥ 1. For any non-negative weights {a0 , · · · , an }, any θ? ∈ L
and any n ≥ 1,
n
X
ak {F (θk ) − min F } ≤ Un (θ? )
k=1

where Tγ and ηn are given by (3) and (6) respectively and
n

1X
Un (θ? ) =
2
def

k=1



ak
ak−1
−
γk
γk−1



kθk−1 − θ? k2 +
−

n
X

a0
kθ0 − θ? k2
2γ0

ak hTγk (θk−1 ) − θ? , ηk i +

k=1

n
X

ak γk kηk k2 . (8)

k=1

Proof See Section 6.2.3.
When
with ηn = 0, Theorem 3 gives
Papplied
Pnan explicit bound of the difference ∆n =
n
A−1
a
F
(θ
)
−
min
F
where
A
=
j
n
n
j=1 j
k=1 ak for the (exact) proximal gradient
sequence {θn , n ∈ N} given by Algorithm 1. When the sequence {an /γn , n ≥ 1} is
−1
non decreasing, (8) shows that ∆n = O(an A−1
n γn ).
Taking ak = 1 for any k ≥ 0 provides a bound for the cumulative regret. When
ak = 1, γk = 1/L for any k ≥ 0, (Schmidt etPal., 2011, Proposition 1) provides a
bound of order O(1) under the assumption that n kηn+1 k < ∞. Using the inequality
| T1/L (θk ) − θ? , ηk+1 | ≤ kθk − θ? kkηk+1 k (see Lemma 9), the upper bound Un (θ? )
in (8) is also O(1).
When an = γn for any n ≥ 0, then supn Un (θ? ) < ∞ under the assumptions that
the series
X
X
γn hTγn (θn−1 ) − θ? , ηn i ,
γn2 kηn k2 ,
n

n

5

converge. In this case, we have

!−1 
 Pn

n
X
γ F (θk )
k=1
.
Pn k
− min F = O 
γk
γ
k
k=1
k=1

Consider the weighted averaged sequence {θ̄n , n ∈ N} defined by
def

θ̄n =

n
1 X
ak θk .
An

(9)

k=1


Pn
Under H1 and H2, F is convex so that F θ̄n ≤ A−1
n
k=1 ak F (θk ). Therefore,
Theorem 3 also provides convergence rates for F (θ̄n ) − min F .

3. Stochastic Proximal Gradient algorithm
In this section, it is assumed that Hn+1 is a Monte Carlo approximation of ∇f (θn ),
where ∇f (θ) satisfies the following assumption:
H 3 for all θ ∈ Θ,
Z
∇f (θ) =

Hθ (x)πθ (dx) ,

(10)

X

for some probability measure πθ on a measurable space (X, X ) and an integrable function (θ, x) 7→ Hθ (x) from Θ × X to Θ.
Note that X is not necessarily a topological space, even if, in many applications,
X ⊆ Rd .
Assumption H3 holds in many problems (see section 4 and section 5). To approximate ∇f (θ), several options are available. Of course, when the dimension of the state
space X is small to moderate, it is always possible to perform a numerical integration
using either Gaussian quadratures or low-discrepancy sequences. Another possibility
is to approximate these integrals: nested Laplace approximations have been considered recently for example in Schelldorfer et al. (2014) and further developed in Ogden
(2015). Such approximations necessarily introduce some bias, which might be difficult
to control. In addition, these techniques are not applicable when the dimension of the
state space X becomes large. In this paper, we rather consider some form of Monte
Carlo approximation.
When sampling πθ is doable, then an obvious choice is to use a naive Monte Carlo
(j)
estimator which amounts to sample a batch {Xn+1 , 1 ≤ j ≤ mn+1 } independently of
the past values of the parameters {θj , j ≤ n} and of the past draws i.e. independently
of the σ-algebra
def
(j)
Fn = σ(θ0 , Xk , 0 ≤ k ≤ n, 0 ≤ j ≤ mk ) .
(11)
We then form

mn+1

Hn+1 =

m−1
n+1

X
j=1

6

(j)

Hθn (Xn+1 ) .

Conditionally to Fn , Hn+1 is an unbiased estimator of ∇f (θn ). The batch size mn+1
can either be chosen to be fixed across iterations or to increase with n at a certain
rate. In the first case, Hn+1 is not converging. In the second case, the approximation
error is vanishing. The fixed batch-size case is closely related to Robbins-Monro
stochastic approximation (the mitigation of the error is performed by letting the
stepsize γn → 0); the increasing batch-size case is related to Monte Carlo assisted
optimisation; see for example Geyer (1994).
The situation that we are facing in section 4 and section 5 is more complicated
because direct sampling from πθ is not an option. Nevertheless, it is fairly easy to
construct a Markov kernel Pθ with invariant distribution πθ . Monte Carlo Markov
Chains (MCMC) provide a set of principled tools to sample from complex distributions
(j)
over large dimensional spaces. In such case, conditional to the past, {Xn+1 , 1 ≤ j ≤
mn+1 } is a realisation of a Markov chain with transition kernel Pθn and started from
(m )
Xn n (the last sample draws in the previous minibatch).
Recall that a Markov kernel P is an application on X × X , taking values in [0, 1]
such that for any x ∈ X, P (x, ·) is a probability measure on X ; and for any A ∈ X ,
x 7→ P (x, A) is measurable. Furthermore, if P is a Markov kernel on X, we denote by
k the k-th iterate of P defined recursively as P 0 (x, A) = 1 (x), and P k (x, A) =
P
A
R k−1
P
(x, dz)P (z, A), k ≥ 1. Finally, the kernel P acts on probability measure: for
any probability measure µ on X , µP is a probability measure defined by
Z
def
µP (A) =
µ(dx)P (x, A),
A ∈ X;
def

def

and P acts on positive measurable functions: for a measurable function f : X → R+ ,
P f is a function defined by
Z
def
P f (x) =
f (y) P (x, dy).
We refer the reader to Meyn and Tweedie (2009) for the definitions and basic properties of Markov chains.
In this Markovian setting, it is possible to consider the fixed batch case and
the increasing batch case. From a mathematical standpoint, the fixed batch case is
trickier, because Hn+1 is no longer an unbiased estimator of ∇f (θn ), i.e. the bias Bn
defined by
i
Pmn+1 h
def
(j)
E
H
(X
)
F
− ∇f (θn )
Bn = E [ Hn+1 | Fn ] − ∇f (θn ) = m−1
n
θ
n
n+1
n+1
j=1
P
(0)
mn+1 j
= m−1
(12)
n+1
j=1 Pθn Hθn (Xn+1 ) − ∇f (θn ) ,
does not vanish. When mn = m is small, the bias can even be pretty large, and the
way the bias is mitigated in the algorithm requires substantial mathematical developments, which are not covered by the results currently available in the literature (see
e.g. Combettes and Pesquet (2015a); Rosasco et al. (2014); Combettes and Pesquet
(2015b); Rosasco et al. (2015); Lin et al. (2015)).
To capture in a common unifying framework these two different situations we
assume that
7

H 4 Hn+1 is a Monte Carlo approximation of the expectation ∇f (θn ) :
mn+1

Hn+1 =

m−1
n+1

X

(j)

Hθn (Xn+1 ) ;

j=1
(j)

for all n ≥ 0, conditionally to the past, {Xn+1 , 1 ≤ j ≤ mn+1 } is a Markov chain
(m )

(m0 )

started from Xn n and with transition kernel Pθn (we set X0
θ ∈ Θ, Pθ is a Markov kernel with invariant distribution πθ .

= x? ∈ X). For all

For a measurable function V : X → [1, ∞), a signed measure µ on the σ-field of X,
and a function f : X → R, define
Z
|f (x)|
def
def
|f |V = sup
,
kµkV = sup
f dµ .
x∈X V (x)
f,|f |V ≤1
H 5 There exist λ ∈ (0, 1), b < ∞, p ≥ 2 and a measurable function W : X → [1, +∞)
such that
sup |Hθ |W < ∞ ,
sup Pθ W p ≤ λW p + b .
θ∈Θ

θ∈Θ

In addition, for any ` ∈ (0, p], there exist C < ∞ and ρ ∈ (0, 1) such that for any
x ∈ X,
(13)
sup kPθn (x, ·) − πθ kW ` ≤ Cρn W ` (x) .
θ∈Θ

Sufficient conditions for the uniform-in-θ ergodic behavior (13) are given e.g. in of
Fort et al. (2011) Lemma 2.3., in terms of aperiodicity, irreducibility and minorization
conditions on the kernels {Pθ , θ ∈ Θ}. Examples of MCMC kernels Pθ satisfying this
assumption can be found in (Andrieu and Moulines, 2006, Proposition 12), (Saksman
and Vihola, 2010, Proposition 15), (Fort et al., 2011, Proposition 3.1.), (Schreck et al.,
2013, Proposition 3.2.), (Allassonnière and Kuhn, 2015, Proposition 1), (Fort et al.,
2015, Proposition 3.1.).
The proof of the results below consists in verifying the conditions of Theorem 2
Pmn+1
(j)
with the error term defined by ηn+1 = m−1
n+1
j=1 Hθn (Xn+1 ) − ∇f (θn ). If the
approximation is unbiased in the sense that E [ ηn+1 | Fn ] = 0, then {ηn , n ∈ N} is a
martingale increment sequence. In all the other cases, we decompose ηn+1 as the sum
of a martingale increment term and a remainder term. When the batch size {mn , n ∈
N} is increasing, the martingale increment sequence can be set to ηn+1 − E [ ηn+1 | Fn ]
and the remainder term E [ ηn+1 | Fn ] will be shown to be vanishingly small. When
the batch size {mn , n ∈ N} is constant, then E [ ηn+1 | Fn ] does not vanish. A more
subtle definition of the martingale increment has to be done, introducing the Poisson
equation for Markov chain (see Proposition 19 in section 6).
3.1 Monte Carlo approximation with fixed batch-size
We first study the case when mn = m for any n ∈ N. Theorem 4 provides sufficient conditions for the convergence towards the limiting set L and for a bound for
P
n
k=1 ak F (θk ) − min F . Consider the following assumption
8

H6

(i) there exists a constant C such that for any θ, θ0 ∈ Θ
|Hθ − Hθ0 |W + sup
x

kPθ (x, ·) − Pθ0 (x, ·)kW
+ kπθ − πθ0 kW ≤ C kθ − θ0 k .
W (x)

(ii) supγ∈(0,1/L] supθ∈Θ γ −1 kProxγ,g (θ) − θk < ∞.
P
(iii)
n |γn+1 − γn | < ∞.
Assumption H6-(i) requires a Lipschitz-regularity in the parameter θ of the Markov
kernel Pθ which, for MCMC algorithms, is inherited under mild additional conditions
from the Lipschitz regularity in W -norm of the target distribution. Such conditions have been worked out for general families of MCMC kernels including HastingsMetropolis dynamics, Gibbs samplers, and hybrid MCMC algorithm; see for example
Proposition 12 in Andrieu and Moulines (2006), the proof of Theorem 3.4. in Fort
et al. (2011), Lemmas 4.6. and 4.7. in Fort et al. (2015) and the references therein. It
is a classical assumption when studying Stochastic Approximation with conditionally
Markovian dynamic (see e.g. Benveniste et al. (1990), Andrieu et al. (2005), Fort et al.
(2014)).
We prove in Proposition 11 that when g is proper, convex, Lipschitz on Θ, then
H6-(ii) is satisfied. In particular, if Θ is a closed convex set, H6-(ii) is satisfied with
the Lasso or fused Lasso penalty. If Θ is a compact convex set, then H6-(ii) is satisfied
by the elastic-net penalty.
For a random variable Y , denote by kY kLq = (E[|Y |q ])1/q .
Theorem 4 Assume Θ is bounded. Let {θn , n ≥ 0} be given by Algorithm 2 with
γn ∈ (0, 1/L] for any n ≥ 0. Assume H1–H5, mn = m ≥ 1 and, if the Monte Carlo
approximation is biased, assume also H6.
P
P
(i) Assume that n γn = ∞ and n γn2 < ∞. With probability one, there exists
θ∞ ∈ L such that limn→∞ θn = θ∞ .
(ii) For any q ∈ (1, p/2] there exists a constant C such that for any non-negative
numbers {a0 , · · · , an }
n
X

ak {F (θk ) − min F }

k=1



Lq
n

a0 X ak
ak−1
≤C +
−
+
γ0
γk
γk−1
k=1

n
X
k=1

!1/2
a2k

+

n
X
k=1

ak γk + υ

n
X


|ak − ak−1 |

k=1

and
n
X

ak {E[F (θk )] − min F }

k=1

≤C

!
n
n
n
X
X
a0 X ak
ak−1
+
−
+
ak γk + υ
|ak − ak−1 |
γ0
γk
γk−1
k=1

k=1

k=1

where υ = 0 if the Monte Carlo approximation is unbiased and υ = 1 otherwise.
9

Proof The proof is postponed to Section 6.3.
When an = 1 and γn = (n + 1)−1/2 , Theorem 4 shows that when n → ∞,


n
X
1
−1
.
n
F (θk ) − min F
=O √
n
q
k=1

L

√
An upper bound O(ln n/ n) can be obtained from Theorem 4 by choosing an = γn =
(n + 1)−1/2 .
3.2 Monte Carlo approximation with increasing batch size
The key property to discuss the asymptotic behavior of the algorithm is the following
result
Proposition 5 Assume H3, H4 and H5. There exists a constant C such that w.p. 1
for any n ≥ 0,
−p/2

(mn )
),
kE [ ηn+1 | Fn ] k ≤ C m−1
n+1 W (Xn

E [kηn+1 kp |Fn ] ≤ C mn+1 W p (Xn(mn ) ) .

Proof The first inequality follows from (12) and (13). The second one is established
in (Fort and Moulines, 2003, Proposition 12).

Theorem 6 Assume Θ is bounded. Let {θn , n ≥ 0} be given by Algorithm 2 with
γn ∈ (0, 1/L] for any n ≥ 0. Assume H1–H5.
P 2
P
−1
γ = +∞,
(i) Assume
n γn+1 mn+1 < ∞ and, if the approximation is biP n n −1
ased, n γn+1 mn+1 < ∞. With probability one, there exists θ∞ ∈ L such that
limn→∞ θn = θ∞ .
(ii) For any q ∈ (1, p/2], there exists a constant C such that for any non-negative
numbers {a0 , · · · , an }
n
X

ak {F (θk ) − min F }

k=1

Lq



n
a0 X ak
ak−1

≤C
+
−
+
γ0
γk
γk−1
k=1

n
X

!1/2
a2k m−1
k

+

k=1

n
X

ak γk m−1
k +υ

k=1

n
X



ak m−1
k

k=1

and
n
X

ak {E[F (θk )] − min F }

k=1
n

≤C

n

n

X
X
a0 X ak
ak−1
+
−
+
ak γk m−1
+υ
ak m−1
k
k
γ0
γk
γk−1
k=1

k=1

!
,

k=1

where υ = 0 if the Monte-Carlo approximation is unbiased and υ = 1 otherwise.
10

Proof See Section 6.4.
Theorem 6 shows that when n → ∞,
n
X
k=1

!−1
ak

n
X


ak F (θk ) − min F

k=1

=O

ln n
n



Lq

by choosing a fixed stepsize γn = γ, a linearly increasing batch-size mn ∼ n and a
uniform weight an = 1. Note thatPthis is the rate after n iterations of the Stochastic
Proximal Gradient algorithm but nk=1 mk = O(n2 ) Monte Carlo samples. Therefore,
√
the rate of convergence expressed in terms of complexity is O(ln n/ n).

4. Application to network structure estimation
To illustrate the algorithm we consider the problem of fitting discrete graphical models
in a setting where the number of nodes in the graph is large compared to the sample
size. Let X be a nonempty finite set, and p ≥ 1 an integer. We consider a graphical
model on Xp with joint probability mass function


p


X
X
1
fθ (x1 , . . . , xp ) =
θkk B0 (xk ) +
θkj B(xk , xj ) ,
(14)
exp


Zθ
k=1

1≤j<k≤p

for a non-zero function B0 : X → R and a symmetric non-zero function B : X × X →
R. The term Zθ is the normalizing constant of the distribution (the partition function), which cannot (in general) be computed explicitly. The real-valued symmetric
matrix θ defines the graph structure and is the parameter of interest. It has the same
interpretation as the precision matrix in a multivariate Gaussian distribution.
We consider the problem of estimating θ from N realizations {x(i) , 1 ≤ i ≤ N }
(i)
(i)
from (14) where x(i) = (x1 , . . . , xp ) ∈ Xp , and where the true value of θ is assumed sparse. This problem is relevant for instance in biology (Ekeberg et al. (2013);
Kamisetty et al. (2013)), and has been considered by many authors in statistics and
machine learning (Banerjee et al. (2008); Höfling and Tibshirani (2009); Ravikumar
et al. (2010); Guo et al. (2010); Xue et al. (2012)).
The main difficulty in dealing with this model is the fact that the log-partition
function log Zθ is intractable in general. As a result, most of the existing works
estimate θ by using the sub-optimal approach of replacing the likelihood function by
a pseudo-likelihood function. One notable exception that tackles the log-likelihood
function is Höfling and Tibshirani (2009), using an active set strategy (to preserve
sparsity), and the junction tree algorithm for computing the partial derivatives of the
log-partition function. However, the success of this strategy depends crucially on the
sparsity of the solution1 . We will see that Algorithm 2 implemented with a MCMC
1. Indeed the implementation of their algorithm in the BMN package is very sensitive to the sparsity
of the solution, and their solver typically fails to converge if the regularization parameter is not
large enough to produce a sufficiently sparse solution. In our numerical experiments, we were not
able to obtain a successful run from their package for p = 100.

11

approximation of the gradient gives a simple and effective approach for computing
the penalized maximum likelihood estimate of θ.
Let Mp denote the space of p×p symmetric matrices equipped with the (modified)
Frobenius inner product
def

X

hθ, ϑi =

def

θjk ϑjk , with norm kθk =

p
hθ, θi.

1≤k≤j≤p

Equipped with this norm, Mp is the same space as the Euclidean space Rd where
d = p(p + 1)/2. Using a `1 -penalty on θ, we see that the computation of the penalized
maximum likelihood estimate of θ is a problem of the form (P) with F = −` + g
where
N
E
X
1 XD
`(θ) =
θ, B̄(x(i) ) − log Zθ and g(θ) = λ
|θjk | ;
N
i=1

1≤k≤j≤p

the matrix-valued function B̄ : Xp → Rp×p is defined by
B̄kk (x) = B0 (xk )

B̄kj (x) = B(xk , xj ) , k 6= j .

It is easy to see
Pthat in this example, Problem (P) admits at least one solution θ?
that satisfies λ 1≤k≤j≤p |θjk | ≤ p log |X|, where |X| denotes the size of X. To see this,
P
(i)
note that since fθ (x) isPa probability, −`(θ) = −N −1 N
i=1 log fθ (x ) ≥ 0. Hence
F (θ) ≥ g(θ) → ∞, as 1≤k≤j≤p |θjk | → ∞ and since F is continuous, we conclude
that it admits at least one minimizer θ? that satisfies F (θ? ) ≤ F (0) = log Z0 =
p log |X|. As a result, and without any loss
P of generality, we consider Problem (P)
with the penalty g replaced by g(θ) = λ 1≤k≤j≤p |θjk | + 1(θ), where 1(θ) = 0 if
maxij |θij | ≤ (p/λ) log |X|, and 1(θ) = +∞ otherwise. Hence in this problem, the
domain of g is Θ = {θ ∈ Mp : maxij |θij | ≤ (p/λ) log |X|}.
Upon noting that (14) is a canonical exponential model, (Shao, 2003, Section
4.4.2) shows that θ 7→ −`(θ) is convex and
∇`(θ) =

Z
N
1 X
B̄(x(i) ) −
B̄(z)fθ (z)µ(dz) ,
N
Xp

(15)

i=1

where µ is the counting measure on Xp . In addition, (see Appendix B)

k∇`(θ) − ∇`(ϑ)k ≤ p (p − 1)osc2 (B) + osc2 (B0 ) kθ − ϑk,

(16)

where for a function B̃ : X × X → R, osc(B̃) = supx,y,u,v∈X |B̃(x, y) − B̃(u, v)|.
Therefore, in this example, the assumption H1 and H2 are satisfied.
The representation of the P
gradient in (15) shows that H3 holds, with πθ (dz) =
N
−1
(i)
fθ (z)µ(dz), and Hθ (z) = N
i=1 B̄(x ) − B̄(z). Direct simulation from the distribution fθ is rarely feasible, so we turn to MCMC. These Markov kernels are easy to
construct, and can be constructed in many ways. For instance if the set X is not too
large, then a Gibbs sampler (see e.g. Robert and Casella (2005)) that samples from
the full conditional distributions of fθ can be easily implemented. In the case of the
12

Gibbs sampler, since Xp is a finite set, Θ is compact, fθ (x) > 0 for all (x, θ) ∈ Xp × Θ,
and, θ 7→ fθ (x) is continuously differentiable, the assumptions H4, H5 and H6(i)-(ii)
automatically hold with W ≡ 1. We should point out that the Gibbs sampler is
a generic algorithm that in some cases is known to mix poorly. Whenever possible
we recommend the use of specialized problem-specific MCMC algorithms with better
mixing properties.
Illustrative example We consider the particular case where X = {1, . . . , M },
B0 (x) = 0, and B(x, y) = 1{x=y} , which corresponds to the well known Potts model.
We report in this section some simulation results showing the performances of the
stochastic proximal gradient algorithm. We use M = 20, B0 (x) = x, N = 250 and for
p ∈ {50, 100, 200}. We generate the “true” matrix θtrue such that it has on average p
non-zero elements below the diagonal which are simulated from a uniform distribution
on (−4, −1) ∪ (1, 4). All the diagonal elements are set to 0.
p
By trial-and-error we set the regularization parameter to λ = 2.5 log(p)/n for all
the simulations. We implement Algorithm 2, drawing samples from a Gibbs sampler
to approximate the gradient. We compare the following two versions of Algorithm 2:
1. Solver 1: A version with a fixed Monte Carlo batch size mn = 500, and decreas1
ing step size γn = 25
p n0.7 .
2. Solver 2: A version with increasing Monte Carlo batch size mn = 500 + n1.2 ,
√1
and fixed step size γn = 25
p 50 .
We run Solver 2 for Niter = 5p iterations, where p ∈ {50, 100, 200} is as above. And
we set the number of iterations of Solver 1 so that both solvers draw approximately
the same number of Monte Carlo samples. For stability in the results, we repeat
the solvers 30 times and average the sample paths. We evaluate the convergence of
each solver by computing the relative error kθn − θ∞ k/kθ∞ k, along the iterations,
where θ∞ denotes the value returned by the solver on its last iteration. Note that we
compare the optimizer output to θ∞ , not θtrue . Ideally, we would like to compare the
iterates to the solution of the optimization problem. However in the present setting a
solution is not available in closed form (and there could be more than one solution).
Furthermore, whether the solution of the optimization problem approaches θ? is a
complicated statistical problem2 that is beyond the scope of this work. The relative
errors are presented on Figure 1 and suggest that, when measured as function of
resource used, Solver 1 and Solver 2 have roughly the same convergence rate.
def

2Senn Precn
We also compute the statistic Fn = Sen
which measures the recovery of the
n +Precn
sparsity structure of θ∞ along the iteration. In this definition Senn is the sensitivity,
and Precn is the precision defined as

P
Senn =

1{|θ |>0} 1{|θ∞,ij |>0}
P n,ij
, and Precn =
j<i 1{|θ∞,ij |>0}

j<i

P

j<i

1{|θn,ij |>0} 1{|θ∞,ij |>0}
.
j<i 1{|θ∞,ij |>0}

P

2. this depends heavily on n, p, the actual true matrix θtrue , and depends also heavily the choice of
the regularization parameter λ

13

0

20

40

60

80

0.2

0.4

0.6

0.8

1.0

Solver 1
Solver 2

0.0

0.2

0.4

0.6

0.8

1.0

Solver 1
Solver 2

0.0

0.0

0.2

0.4

0.6

0.8

1.0

Solver 1
Solver 2

0

50

100

150

0

500

1000

1500

2000

2500

3000

0

20

40

60

80

1.10
0.85

0.90

0.95

1.00

1.05

Solver 1
Solver 2

0.80

0.85

0.90

0.95

1.00

1.05

Solver 1
Solver 2

0.80

0.80

0.85

0.90

0.95

1.00

1.05

Solver 1
Solver 2

1.10

1.10

Figure 1: Relative errors plotted as function of computing time for Solver 1 and Solver
2.

0

50

100

150

0

500

1000

1500

2000

2500

3000

Figure 2: Statistic Fn plotted as function of computing time for Solver 1 and Solver
2.

The values of Fn are presented on Figure 2 as function of computing time. It shows
that for both solvers, the sparsity structure of θn converges very quickly towards that
of θ∞ . We note also that Figure 2 seems to suggest that Solver 2 tends to produce
solutions with slightly more stable sparsity structure than Solver 1 (less variance on
the red curves). Whether such subtle differences exist between the two algorithms (a
diminishing step-size and fixed Monte Carlo size versus a fixed step-size and increasing
Monte Carlo size) is an interest question. Our analysis does not deal with the sparsity
structure of the solutions, hence cannot offer any explanation.

5. A non convex example: High-dimensional logistic regression with
random effects
We numerically investigate the extension of our results to a situation where the assumptions H2 and H3 hold but H1 is not in general satisfied and the domain Θ is
not bounded. The numerical study below shows that the conclusions reached in sec14

tion 2 and section 3 provide useful information to tune the design parameters of the
algorithms.
5.1 The model
We model binary responses {Yi }N
i=1 ∈ {0, 1} as N conditionally independent realizations of a random effect logistic regression model,

ind.
Yi |U ∼ Ber s(x0i β + σzi0 U) , 1 ≤ i ≤ N ,

(17)

where xi ∈ Rp is the vector of covariates, zi ∈ Rq are (known) loading vector, Ber(α)
denotes the Bernoulli distribution with parameter α ∈ [0, 1], s(x) = ex /(1 + ex ) is the
cumulative distribution function of the standard logistic distribution. The random
effect U is assumed to be standard Gaussian U ∼ Nq (0, I).
The log-likelihood of the observations at θ = (β, σ) ∈ Rp × (0, ∞) is given by
`(θ) = log

Z Y
N

1−Yi
φ(u)du ,
s(x0i β + σzi0 u)Yi 1 − s(x0i β + σzi0 u)

(18)

i=1

where φ is the density of a Rq -valued standard Gaussian random vector. The number
of covariates p is possibly larger than N , but only a very small number of these
covariates are relevant which suggests to use the elastic-net penalty


1−α
2
λ
kβk2 + αkβk1 ,
(19)
2
P
where λ > 0 is the regularization parameter, kβkr = ( pi=1 |βi |r )1/r and α ∈ [0, 1]
controls the trade-off between the `1 and the `2 penalties. In this example,


1−α
2
g(θ) = λ
kβk2 + αkβk1 + 1(0,+∞) (σ) ,
(20)
2
where 1A (x) = +∞ is x ∈
/ A and 0 otherwise. Define the conditional log-likelihood
of Y = (Y1 , . . . , YN ) given U (the dependence upon Y is omitted) by
`c (θ|u) =

N
X




Yi x0i β + σzi0 u − ln 1 + exp x0i β + σzi0 u

,

i=1

and the conditional distribution of the random effect U given the observations Y and
the parameter θ
πθ (u) = exp (`c (θ|u) − `(θ)) φ(u) .
(21)
The Fisher identity implies that the gradient of the log-likelihood (18) is given by
Z
∇`(θ) =

 )
Z (X
N
x
∇θ `c (θ|u) πθ (u) du =
(Yi − s(x0i β + σzi0 u)) 0 i
πθ (u) du .
zi u
i=1

15

The Hessian of the log-likelihood ` is given by (see e.g.(McLachlan and Krishnan,
2008, Chapter 3))


∇2 `(θ) = Eπθ ∇2θ `c (θ|U) + Covπθ (∇θ `c (θ|U))
where Eπθ and Covπθ denotes the expectation and the covariance with respect to the
distribution πθ , respectively. Since
∇2θ `c (θ|u)

=−

N
X

s(x0i β

σzi0 u)

+

1−

s(x0i β

+

   0
 xi
xi
,
zi0 u zi0 u

σzi0 u)

i=1

R
and supθ∈Θ kuk2 πθ (u) du < ∞ (see Appendix A), ∇2 `(θ) is bounded on Θ. Hence,
∇`(θ) satisfies the Lipschitz condition showing that H1 is satisfied.
5.2 Numerical application
The assumption H3 is satisfied with πθ given by (21) and
Hθ (u) = −

N
X

(Yi −

F (x0i β

+

σzi0 u))

i=1




xi
.
zi0 u

(22)

The distribution πθ is sampled using the MCMC sampler
R proposed in Polson et al.
(2013) based on data-augmentation. We write −∇`(θ) = Rq ×RN Hθ (u)π̃θ (u, w) dudw
where π̃θ (u, w) is defined for u ∈ Rq and w = (w1 , · · · , wN ) ∈ RN by
!
N
Y

0
0
π̃θ (u, w) =
π̄PG wi ; xi β + σzi u πθ (u) ;
i=1

in this expression, π̄PG (·; c) is the density of the Polya-Gamma distribution on the
positive real line with parameter c given by

π̄PG (w; c) = cosh(c/2) exp −wc2 /2 ρ(w)1R+ (w) ,
P
where ρ(w) ∝ k≥0 (−1)k (2k +1) exp(−(2k +1)2 /(8w))w−3/2 (see (Biane et al., 2001,
Section 3.1)). Thus, we have
π̃θ (u, w) = Cθ φ(u)

N
Y


exp σ(Yi − 1/2)zi0 u − wi (x0i β + σzi0 u)2 /2 ρ(wi )1R+ (wi ) ,

i=1

P
0
where ln Cθ = −N ln 2 − `(θ) + N
i=1 (Yi − 1/2)xi β. This target distribution can be
sampled using a Gibbs algorithm: given the current value (ut , wt ) of the chain, the
next point is obtained by sampling ut+1 under the conditional distribution of u given
wt , and wt+1 under the conditional distribution of w given ut+1 . In the present case,
these conditional distributions are given respectively by
π̃θ (u|w) ≡ Nq (µθ (w); Γθ (w))

π̃θ (w|u) =

N
Y
i=1

16

π̄PG (wi ; |x0i β + σzi0 u|)

with
Γθ (w) =

I +σ

2

N
X

!−1
wi zi zi0

,

µθ (w) = σΓθ (w)

i=1

N
X


(Yi − 1/2) − wi x0i β zi .

i=1

(23)
Exact samples of these conditional distributions can be obtained (see (Polson et al.,
2013, Algorithm 1) for sampling under a Polya-Gamma distribution). It has been
shown by Choi and Hobert (2013) that the Polya-Gamma Gibbs sampler is uniformly
ergodic. Hence H5 is satisfied with W ≡ 1. Checking H6 is also straightforward.
We test the algorithms with N = 500, p = 1, 000 and q = 5. We generate the N ×p
covariates matrix X columnwise, by sampling a stationary
RN -valued autoregressive
p
model with parameter ρ = 0.8 and Gaussian noise 1 − ρ2 NN (0, I). We generate
the vector of regressors βtrue from the uniform distribution on [1, 5] and randomly set
98% of the coefficients to zero. The variance of the random effect is set to σ 2 = 0.1.
We consider a repeated measurement setting so that zi = ediq/N e where {ej , j ≤ q} is
the canonical basis of Rq and d·e denotes the upper integer part. With such a simple
expression for the random effect, we will be able to approximate the value F (θ) in
order to illustrate the theoretical results obtained in this paper. We use the Lasso
penalty (α = 1 in (19)) with λ = 30.
We first illustrate the ability of Monte Carlo Proximal Gradient algorithms to find
a minimizer of F . We compare the Monte Carlo proximal gradient algorithm
√
(i) with fixed batch size: γn = 0.01/ n and mn = 275 (Algo 1); γn = 0.5/n and
mn = 275 (Algo 2).
(ii) with increasing batch size: γn = γ = 0.005, mn = 200 + n (Algo 3); γn =
√
√
γ = 0.001, mn = 200 + n (Algo 4); and γn = 0.05/ n and mn = 270 + d ne
(Algo 5).
Each algorithm is run for 150 iterations. The batch sizes {mn , n ≥ 0} are chosen
so that after 150 iterations, each algorithm used approximately the same number
of Monte Carlo samples. We denote by β∞ the value obtained at iteration 150. A
path of the relative error kβn − β∞ k/kβ∞ k is displayed on Figure 3[right] for each
algorithm; a path of the sensitivity Senn and of the precision Precn (see section 4 for
the definition) are displayed on Figure 4. All these sequences are plotted versus the
total number of Monte Carlo samples up to iteration n. These plots show that with
a fixed batch-size (Algo 1 or Algo 2), the best convergence is obtained with a step
√
size decreasing as O(1/ n); and for an increasing batch size (Algo 3 to Algo 5), it
is better to choose a fixed step size. These findings are consistent with the results
in section 3. On Figure 3[left], we report on the bottom row the indices j such that
βtrue,j is non null and on the rows above, the indices j such that β∞,j given by Algo 1
to Algo 5 is non null.
We now study the convergence of {F (θn ), n ∈ N} where θn is obtained by one of
the algorithms described above. We repeat 50 independent runs for each algorithm
and estimate E [F (θn )] by the empirical mean over these runs. On Figure 5[left],
n 7→ F (θn ) is displayed for several runs of Algo 1 and Algo 3. The figure shows
17

4

10
Algo 1

Algo 1
Algo 2

2

10

0

10
Algo 2

−2

10

−4

Algo 3

10

0

0.5

1

1.5

2

2.5

3

3.5

4
4

x 10
4

10

Algo 4

Algo 3
Algo 4
Algo5

2

10
Algo 5

0

10

−2

10

Beta True

−4

0

100

200

300

400

500

600

700

800

900

1000

10

0

0.5

1

1.5

2

2.5

3

3.5

4
4

x 10

Figure 3: [left] The support of the sparse vector β∞ obtained by Algo 1 to Algo 5;
for comparison, the support of βtrue is on the bottom row. [right] Relative
error along one path of each algorithm as a function of the total number of
Monte Carlo samples.

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2
0

0.2

Algo 1
Algo 2
0

0.5

1

1.5

2

2.5

3

3.5

0

4

Algo 1
Algo 2
0

0.5

1

1.5

2

2.5

3

3.5

4

x 10

1

1

0.8

0.8

0.6

0.6

0.4

0.4
Algo 3
Algo 4
Algo5

0.2
0

4
4

x 10

0

0.5

1

1.5

2

2.5

3

3.5

Algo 3
Algo 4
Algo5

0.2
0

4

0

4

0.5

1

1.5

2

2.5

3

3.5

4
4

x 10

x 10

Figure 4: The sensitivity Senn [left] and the precision Precn [right] along a path,
versus the total number of Monte Carlo samples up to time n

18

5

4

10

10

Algo 1
Algo 3
Algo 4

Algo 1
Algo 3
4

10

3

10

2

10

0

50

100

150

5

10

3

Algo 1
Algo 3

10

4

10

3

10

2

10

0

0.5

1

1.5

2

2.5

3

3.5

4

0

0.5

1

1.5

4

2

2.5

3

3.5

4
4

x 10

x 10

Figure 5: [left] n 7→ F (θn ) for several independent runs. [right] E [F (θn )] versus the
total number of Monte Carlo samples up to iteration n

that all the paths have the same limiting value, which is approximately F? = 311;
we observed the same behavior on the 50 runs of each algorithm. On Figure 5[right],
we report the Monte Carlo estimation of E[F (θn )] versus the total number of Monte
Carlo samples used up to iteration n for the best strategies in the fixed batch size
case (Algo 1) and in the increasing batch size case (Algo 3 and Algo 4).

6. Proofs
6.1 Preliminary lemmas
Lemma 7 Assume that g is lower semi-continuous and convex. For θ, θ0 ∈ Θ and
γ>0


1
g Proxγ,g (θ) − g(θ0 ) ≤ − Proxγ,g (θ) − θ0 , Proxγ,g (θ) − θ .
(24)
γ
For any γ > 0 and for any θ, θ0 ∈ Θ,


k Proxγ,g (θ)−Proxγ,g (θ0 )k2 +k Proxγ,g (θ)−θ − Proxγ,g (θ0 )−θ0 k2 ≤ kθ−θ0 k2 . (25)
Proof See (Bauschke and Combettes, 2011, Propositions 4.2., 12.26 and 12.27).

Lemma 8 Assume H1 and let γ ∈ (0, 1/L]. Then for all θ, θ0 ∈ Θ,


−2γ F (Proxγ,g (θ))−F (θ0 ) ≥ k Proxγ,g (θ)−θ0 k2 +2 Proxγ,g (θ) − θ0 , θ0 − γ∇f (θ0 ) − θ .
(26)
0
If in addition f is convex, then for all θ, θ , ξ ∈ Θ,



− 2γ F Proxγ,g (θ) − F (θ0 ) ≥ k Proxγ,g (θ) − θ0 k2
+ 2 Proxγ,g (θ) − θ0 , ξ − γ∇f (ξ) − θ − kθ0 − ξk2 . (27)
19

Proof Since ∇f is Lipschitz, the descent lemma implies that for any γ −1 ≥ L
f (p) − f (θ0 ) ≤ ∇f (θ0 ), p − θ0 +

1
kp − θ0 k2 .
2γ

(28)

This inequality applied with p = Proxγ,g (θ) combined with (24) yields (26). When
f is convex, f (ξ) + h∇f (ξ), θ0 − ξi − f (θ0 ) ≤ 0 which, combined again with (24) and
(28) applied with (p, θ0 ) ← (Proxγ,g (θ), ξ) yields the result.

Lemma 9 Assume H1. Then for any γ > 0, θ, θ0 ∈ Θ,
kθ − γ∇f (θ) − θ0 + γ∇f (θ0 )k ≤ (1 + γL)kθ − θ0 k ,
0

0

kTγ (θ) − Tγ (θ )k ≤ (1 + γL)kθ − θ k .

(29)
(30)

If in addition f is convex then for any γ ∈ (0, 2/L],
kθ − γ∇f (θ) − θ0 + γ∇f (θ0 )k ≤ kθ − θ0 k ,
0

(31)

0

kTγ (θ) − Tγ (θ )k ≤ kθ − θ k .

(32)

Proof (30) and (32) follows from (29) and (31) respectively by the Lipschitz property
of the proximal map Proxγ,g (see Lemma 7). (29) follows directly from the Lipschitz
property of f . It remains to prove (31). Since f is a convex function with Lipschitzcontinuous gradients, (Nesterov, 2004, Theorem 2.1.5) shows that, for all θ, θ0 ∈ Θ,
L h∇f (θ) − ∇f (θ0 ), θ − θ0 i ≥ k∇f (θ) − ∇f (θ0 )k2 . The result follows.

def

def

Lemma 10 Assume H1. Set Sγ (θ) = Proxγ,g (θ − γH) and η = H − ∇f (θ). For
any θ ∈ Θ and γ > 0,
kTγ (θ) − Sγ (θ)k ≤ γkηk .
(33)
Proof We have kTγ (θ) − Sγ (θ)k = k Proxγ,g (θ − γ∇f (θ)) − Proxγ,g (θ − γH)k and
(33) follows from Lemma 7.

6.2 Proof of section 2
6.2.1 Proof of Lemma 1
P
P
def
Set wn = vn + k≥n+1 ξk + M with M = − inf n k≥n ξk so that inf n wn ≥ 0. Then
0 ≤ wn+1 ≤ vn − χn+1 + ξn+1 +

X

ξk + M ≤ wn − χn+1 .

k≥n+2

{wn ,Pn ∈ N} is non-negative
Pand non increasing; therefore it converges. Furthermore,
n
0 ≤ k=0 χk ≤ w0 so that n χn < ∞. The convergence of {wn , n ∈ N} also implies
the convergence of {vn , n ∈ N}. This concludes the proof.
20

6.2.2 Proof of Theorem 2
Let θ? ∈ L, which is not empty by H2; note that F (θ? ) = min F . We have by (27)
applied with θ ← θn − γn+1 Hn+1 , ξ ← θn , θ0 ← θ? , γ ← γn+1
kθn+1 − θ? k2 ≤ kθn − θ? k2 − 2γn+1 (F (θn+1 ) − min F ) − 2γn+1 hθn+1 − θ? , ηn+1 i .
We write θn+1 − θ? = θn+1 − Tγn+1 (θn ) + Tγn+1 (θn ) − θ? . By Lemma 10, kθn+1 −
Tγn+1 (θn )k ≤ γn+1 kηn+1 k so that,
− hθn+1 − θ? , ηn+1 i ≤ γn+1 kηn+1 k2 − Tγn+1 (θn ) − θ? , ηn+1 .
Hence,
kθn+1 − θ? k2 ≤ kθn − θ? k2 − 2γn+1 (F (θn+1 ) − min F )
2
+ 2γn+1
kηn+1 k2 − 2γn+1 Tγn+1 (θn ) − θ? , ηn+1 . (34)
P
− min F ) < ∞ and limn kθn −
Under (7) and (34), Lemma 1 shows that n γn (F (θn ) P
θ? k exists. This implies that supn kθn k < ∞. Since n γn = +∞, there exists a
subsequence {θφn , n ∈ N} such that limn F (θφn ) = min F . The sequence {θφn , n ≥ 0}
being bounded, we can assume without loss of generality that there exists θ∞ ∈ Rd
such that limn θφn = θ∞ .
Let us prove that θ∞ ∈ L. Since g is lower semi-continuous on Θ, lim inf n g(θφn ) ≥
g(θ∞ ) so that θ∞ ∈ Θ. Since F is lower semi-continuous on Θ, we have

min F = lim inf F (θφn ) ≥ F (θ∞ ) ≥ min F ,
n→∞

showing that F (θ∞ ) = min F .
By (34), for any m and n ≥ φm
2

2

kθn+1 − θ∞ k ≤ kθφm − θ∞ k − 2

n
X

γk+1 { Tγk+1 (θk ) − θ∞ , ηk+1 + γk+1 kηk+1 k2 } .

k=φm

For any  > 0, there exists m such that the RHS is upper bounded by . Hence, for
any n ≥ φm , kθn+1 − θ∞ k2 ≤ , which proves the convergence of {θn , n ∈ N} to θ∞ .
6.2.3 Proof of Theorem 3
Let θ? ∈ L; note that F (θ? ) = min F . We first apply (27) with θ ← θj − γj+1 Hj+1 ,
ξ ← θj , θ0 ← θ? , γ ← γj+1 :

F (θj+1 ) − min F ≤ (2γj+1 )−1 kθj − θ? k2 − kθj+1 − θ? k2 − hθj+1 − θ? , ηj+1 i .
Multiplying both sides by aj+1 gives:


aj+1 F (θj+1 ) − min F



1
≤
2




aj+1 aj
aj
−
kθj − θ? k2
kθj − θ? k2 +
γj+1 γj
2γj
aj+1
−
kθj+1 − θ? k2 − aj+1 hθj+1 − θ? , ηj+1 i .
2γj+1
21

Summing from j = 0 to n − 1 gives
n

n

X
an
1X
kθn − θ? k2 +
aj {F (θj ) − min F } ≤
2γn
2
j=1

j=1

−

n
X



aj−1
aj
−
γj
γj−1



aj hθj − θ? , ηj i +

j=1

kθj−1 − θ? k2

a0
kθ0 − θ? k2 . (35)
2γ0

We decompose hθj − θ? , ηj i as follows:
hθj − θ? , ηj i = θj − Tγj (θj−1 ), ηj + Tγj (θj−1 ) − θ? , ηj
By Lemma 10, we get

.

≤ γj kηj k2 which concludes the proof.

θj − Tγj (θj−1 ), ηj

6.3 Proof of Section 3.1
The proof of Theorem 4 is given in the case m = 1; we simply denote by Xn the
(1)
sample Xn . The proof for the case m > 1 can be adapted from the proof below, by
substituting the functions Hθ (x) and W (x) by
m

H θ (x1 , · · · , xm ) =

m

1 X
Hθ (xk )
m

W (x1 , · · · , xm ) =

k=1

1 X
W (xk ) ;
m
k=1

the kernel Pθ and its invariant measure πθ by
Z
P θ (x1 , · · · , xm ; B) =

Z
···

Pθ (xm , dy1 )

m
Y

Pθ (yk−1 , dyk )1B (y1 , . . . , ym ) ,

k=2

Z
π θ (B) =

Z
···

πθ (dy1 )

m
Y

Pθ (yk−1 , dyk )1B (y1 , . . . , ym ) ,

k=2

for any (x1 , . . . , xm ) ∈ Xn and B ∈ X ×n .
6.3.1 Preliminary results
Proposition 11 Assume that g is proper convex and Lipschitz on Θ with Lipschitz
constant K. Then, for all θ ∈ Θ,
kProxγ,g (θ) − θk ≤ Kγ .
Proof For all θ ∈ Θ, we get by Lemma 7
0 ≤ γ −1 kθ − Proxγ,g (θ)k2 ≤ g(θ) − g(Proxγ,g (θ)) ≤ K kθ − Proxγ,g (θ)k .

22

(36)

Proposition 12 Assume H1, H2 and Θ is bounded. Then
sup kTγ (θ)k < ∞ .

sup

γ∈(0,1/L] θ∈Θ

If in addition H6-(ii) holds, then there exists a constant C such that for any θ, θ̄ ∈ Θ,
γ, γ̄ ∈ (0, 1/L]

Tγ (θ) − Tγ̄ (θ̄) ≤ C γ + γ̄ + kθ − θ̄k .
Proof Let θ? such that for any γ > 0, θ? = Tγ (θ? ) (such a point exists by H2 and
(4)). We write Tγ (θ) = (Tγ (θ) − θ? )+θ? . By Lemma 9, there exists a constant C such
that for any θ ∈ Θ and any γ ∈ (0, 1/L], kTγ (θ) − θ? k ≤ 2 kθ − θ? k ≤ 2 kθk + 2 kθ? k.
This concludes the proof of the first statement. We write Tγ (θ) − Tγ̄ (θ̄) = Tγ (θ) −
Tγ̄ (θ) + Tγ̄ (θ) − Tγ̄ (θ̄). By Lemma 7
Tγ̄ (θ) − Tγ̄ (θ̄) ≤ θ − θ̄ − γ̄∇f (θ) + γ̄∇f (θ̄) ≤ kθ − θ̄k + γ̄ sup k∇f (θ)k .
θ∈Θ

By H1 and since Θ is bounded, supθ∈Θ k∇f (θ)k < ∞. In addition, using again
Lemma 7,
kTγ (θ) − Tγ̄ (θ)k ≤ (γ + γ̄) sup k∇f (θ)k + kProxγ,g (θ) − Proxγ̄,g (θ)k .
θ∈Θ

We conclude by using
kProxγ̄,g (θ) − Proxγ,g (θ)k ≤ kProxγ̄,g (θ) − θk + kθ − Proxγ,g (θ)k
≤ (γ + γ̄)

sup

sup γ −1 kProxγ,g (θ) − θk .

γ∈(0,1/L] θ∈Θ

Lemma 13 Assume H5 and H6-(i).
b θ (x) such that supθ∈Θ H
bθ
(i) There exists a measurable function (θ, x) 7→ H
and for any (θ, x) ∈ Θ × X,

W

<∞

Z
b θ (x) − Pθ H
b θ (x) = Hθ (x) −
H

Hθ (y)πθ (dy) .

(ii) There exists a constant C such that for any θ, θ0 ∈ Θ,
b θ − Pθ0 H
b θ0
Pθ H
Proof See (Fort et al., 2011, Lemma 4.2).

23

W

≤C

θ − θ0

.

(37)

Lemma 14 Assume H4 and H5. Then, supn E [W p (Xn )] < ∞.
Proof Conditionnally to the past Fj−1 , the conditional distribution of Xj is Pθj−1 (Xj−1 , ·).
Therefore, we write


E [W p (Xn )] = E [E [ W p (Xn ) | Fn−1 ]] = E Pθn−1 W p (Xn−1 ) .
We then use the drift inequality to obtain E [W p (Xn )] ≤ λE [W p (Xn−1 )] + b. The
proof then follows from a trivial induction.

Lemma 15 Assume H1, H6-(ii) and Θ is bounded. There exists a constant C such
that w.p.1, for all n ≥ 0,
kθn+1 − θn k ≤ Cγn+1 (1 + kηn+1 k) .
Proof We write
θn+1 − θn = θn+1 − Proxγn+1 ,g (θn ) + Proxγn+1 ,g (θn ) − θn .
Since by Lemma 7, θ 7→ Proxγ,g (θ) is Lipschitz for any γ > 0, we get
θn+1 − Proxγn+1 ,g (θn ) = Proxγn+1 ,g (θn − γn+1 ηn+1 − γn+1 ∇f (θn )) − Proxγn+1 ,g (θn )


≤ γn+1 kηn+1 + ∇f (θn )k ≤ γn+1 kηn+1 k + sup k∇f (θ)k .
θ∈Θ

By H1, w.p.1. supθ∈Θ k∇f (θ)k < ∞; hence, there exists C1 such that w.p.1. for all
n ≥ 0, θn+1 − Proxγn+1 ,g (θn ) ≤ C1 γn+1 (1 + kηn+1 k). Finally, under H6-(ii), there
exists a constant C2 such that, w.p.1.,
−1
sup γn+1
Proxγn+1 ,g (θn ) − θn ≤
n

sup

sup γ −1 kProxγ,g (θ) − θk ≤ C2 .

γ∈(0,1/L] θ∈Θ

This concludes the proof.

Lemma 16 Assume H1, H4, H5 and Θ is bounded. There exists a constant C such
that w.p.1, for all n ≥ 0, kηn+1 k ≤ CW (Xn+1 ).
Proof By H4 and H5, kηn+1 k ≤ (supθ∈Θ |Hθ |W ) W (Xn+1 ) + supθ∈Θ k∇f (θ)k. The
result follows since ∇f is Lipschitz by H1, and since W ≥ 1.

6.3.2 Proof of Theorem 4
The proof of the almost-sure convergence consists in verifying
assumptions of
P the
2 kη
2
Theorem 2. Let us start with the proof that almost-surely, n γn+1
n+1 k < ∞.
24

This property is a consequence of Lemma 17 applied with an ← γn2 . It remains to
prove that almost-surely
X
X
γn ηn < ∞,
γn+1 Tγn+1 (θn ), ηn+1 < ∞;
n

n

P

note that they are both of the form n γn+1 Aγn+1 (θn )ηn+1 with, respectively, Aγ (θ)
equal to the identity matrix, and Aγ (θ) = Tγ (θ). In the case the Monte Carlo is unbiased, we apply Proposition 18 with an ← γn and
P Aγ (θ) equal to the identity matrix
and we obtain the almost-sure convergence of n γn ηn ; we then apply Proposition 18
with
P an ← γn and Aγ (θ) = Tγ (θ), and we obtain the almost-sure convergence of
n γn+1 Tγn+1 (θn ), ηn+1 - note that by Proposition 12, Tγ (θ) satisfies the assumptions on Aγ (θ). In the case the Monte Carlo is biased, the steps are the same except
we use Proposition 19 instead of Proposition 18.
For the control of the moments, we use Theorem 3 and again Lemma 17 and
Proposition 18 for the unbiased case (or Proposition 19 for the biased case).
Lemma 17 Assume H1, H4, H5 and Θ is bounded.
P
P
2
(i) If ak ≥ 0 and ∞
n≥1 an kηn k < ∞.
k=1 ak < ∞ then with probability one,
(ii) for any q ∈ [1, p/2], there exists a constant C such that for any non-negative
numbers {a1 , · · · , an },
n
X

ak kηk k

2

k=1

≤C
Lq

n
X

ak .

k=1

Proof We write

E


X


 X
an+1 kηn+1 k2  ≤ sup E kηn+1 k2
an+1 .
n

n≥0

n≥0

By Lemma 14 and Lemma 16, supn kηn+1 kL2 < ∞ so the RHS is finite. By the
Minkovski inequality, we write since ak > 0,
n
X

≤ sup kηn k2L2q

ak+1 kηk+1 k2

k=0

Lq

n

n+1
X

ak .

k=1

The supremum is finite by Lemma 14 and Lemma 16.

Proposition 18 Assume H1, H3, H4, H5, Θ is bounded and the Monte Carlo approximation is unbiased. Let {an , n ∈ N} be a deterministic positive sequence and
{Aγ (θ), γ ∈ (0, 1/L] , θ ∈ Θ} be deterministic matrices such that
sup

sup kAγ (θ)k < ∞ .

γ∈(0,1/L] θ∈Θ

25

(38)

(i) If

2
n≥0 an

P

< ∞, then the series

P

n≥0 an+1 Aγn+1 (θn )ηn+1

converges P-a.s.

(ii) For any q ∈ (1, p/2], there exists a constant C such that
n
X

≤C

ak+1 Aγk+1 (θk )ηk+1

k=0

n
X

!1/2
a2k+1

.

k=0

Lq



Proof P
Since θn ∈ Fn , we have E an+1 Aγn+1 (θn ) ηn+1 |Fn = 0, thus showing that
{Mn = nk=0 ak+1 Aγk+1
P (θk )ηk+1 , n ∈ N} is a martingale. This martingale converges
almost-surely if S = n≥0 a2n+1 kAγn+1 (θn )k2 kηn+1 k2 < ∞ P-a.s. (see e.g. (Hall and
Heyde, 1980, Theorem 2.17)). Using (38) and Lemma 17, S < ∞ P-a.s.
Consider now the Lq -moment of Mn . We apply (Hall and Heyde, 1980, Theorem
2.10): for any q ∈ (1, p/2], there exists a constant C such that for any n ≥ 0,
n
X

≤C

ak+1 Aγk+1 (θk )ηk+1

k=0

n
X

!1/2
2
ak+1 Aγk+1 (θk )ηk+1 Lq

.

k=0

Lq

Lemma 14 and Lemma 16 imply that supn kηn+1 kLq < ∞; we then conclude with
(38).

Proposition 19 Assume H1, H3–H6 and Θ is bounded. Let {an , n ≥ 0} be a positive
sequence and {Aγ (θ), γ ∈ (0, 1/L] , θ ∈ Θ} be (deterministic) function-valued matrices
such that there exists CA and for any γ, γ̄ ∈ (0, 1/L] and θ, θ̄ ∈ Θ

Aγ (θ) − Aγ̄ (θ̄) ≤ CA γ + γ̄ + θ − θ̄ . (39)
sup sup kAγ (θ)k < ∞ ,
γ∈(0,1/L] θ∈Θ

P 2
P
< ∞,
n an < ∞ and
n |an+1 − an | < ∞ then the series
a
A
(θ
)η
converges
P-a.s.
n≥0 n+1 γn+1 n n+1

(i) P
If

P

n an γn

(ii) For any q ∈ (1, p/2], there exists a constant C such that
n
X
k=0

ak+1 Aγk+1 (θk )ηk+1
Lq



≤C 1+


n
X
k=0

!1/2
a2k+1

+

n
X

|ak+1 − ak | +

k=1

Proof
(i) By H4 and Lemma 13-(i), we write
b θ (Xn+1 ) − Pθ H
b θ (Xn+1 )
ηn+1 = H
n
n
n

 

b θ (Xn ) + Pθ H
b θ (Xn ) − Pθ H
b θ (Xn+1 )
b θ (Xn+1 ) − Pθ H
= H
n
n
n
n
n
n+1
n+1


b θ (Xn+1 ) − Pθ H
b θ (Xn+1 ) .
+ Pθn+1 H
n
n
n+1
26

n
X
k=1

ak γk





.

We prove successively that w.p.1,

b θ (Xn+1 ) − Pθ H
b θ (Xn ) < ∞ ,
H
n
n
n

(40)


b θ (Xn+1 ) < ∞ ,
b θ (Xn ) − Pθ H
Pθn H
n
n+1
n+1

(41)



b θ (Xn+1 ) − Pθ H
b θ (Xn+1 ) < ∞ .
an+1 Aγn+1 (θn ) Pθn+1 H
n
n
n+1

(42)

X

an+1 Aγn+1 (θn )



n

X

an+1 Aγn+1 (θn )



n≥0

X
n≥0

b θ (Xn+1 ) − Pθ H
b θ (Xn ), n ∈ N} is a martingale increment
Proof of (40) By H4, {H
n
n
n
w.r.t. the filtration {Fn , n ≥ 0}. The proof is along the same lines as the proof of
Proposition 18 upon noting that by Lemma 13 and H5, there exists C such that w.p.1
for all n ≥ 0,
b θ (Xn )k ≤ C (W (Xn+1 ) + W (Xn )) .
b θ (Xn+1 ) − Pθ H
kH
n
n
n
P
b θ (Xn ) with ∆n+1 = an+1 Aγ (θn )−
Proof of (41) The sum is equal to n≥0 ∆n+1 Pθn H
n+1
n
an Aγn (θn−1 ). On one hand, by Lemma 13 and H5, there exists C such that w.p.1 for
all n ≥ 0,
b θ (Xn ) ≤ C W (Xn ) .
Pθn H
n
On the other hand, by (39), Lemma 15 and Lemma 16, there exists C such that a.s.


for all n ≥ 0,
k∆n+1 k ≤ C |an+1 − an | + an (γn + γn+1 ) W (Xn ) .


By Lemma 14, supn E W 2 h(Xn ) < ∞. Therefore,
i by (39) and the assumptions on
P
b
{an , n ≥ 0}, we have n E ∆n+1 Pθn Hθn (Xn ) < ∞; which concludes the proof.
Proof of (42) By (39) and Lemma 13, there exists a constant C such that w.p.1 for
any n


b θ (Xn+1 ) ≤ C kθn+1 − θn k W (Xn+1 ) .
b θ (Xn+1 ) − Pθ H
Aγn+1 (θn ) Pθn+1 H
n
n
n+1
By Lemma 15 and Lemma 16, there exists a constant C such that w.p.1,
kθn+1 − θn k W (Xn+1 ) ≤ Cγn+1 W 2 (Xn+1 ) .


P
From Lemma 14 and the assumptions on {an , n ≥ 0}, n an+1 γn+1 E W 2 (Xn+1 ) <
∞ from which (42) follows.
for all n ≥ 0,

(ii) We start from the same decomposition of ηn+1 in three terms. The first one
is a martingale, and following the same lines as in the proof of Proposition 18, we
obtain
!1/2
n
n


X
X
b θ (Xn+1 ) − Pθ H
b θ (Xn )
ak+1 Aγ (θn ) H
≤C
a2
.
n+1

n

n

k=0

k+1

n

Lq

27

k=0

For the second term, we write
n
X

ak+1 Aγk+1 (θk )




b θ (Xk+1 )
b θ (Xk ) − Pθ H
Pθk H
k
k+1
k+1

k=0

b θ (Xn+1 )
b θ (X0 ) − an+1 Aγ (θn )Pθ H
≤ a1 Aγ1 (θ0 )Pθ0 H
n+1
n+1
n+1
0
n
X
b θ (Xk ) .
+
∆k+1 Pθk H
k
k=1

By the Minkovski inequality, it is easily seen that there exists a constant C such that
n
X

ak+1 Aγk+1 (θk )




b θ (Xk ) − Pθ H
b θ (Xk+1 )
Pθk H
k
k+1
k+1

k=0

Lq

≤

1 + an+1 +

n 
X

|ak+1 − ak | + ak (γk + γk+1 )



!
.

k=1

Finally, for the last term, following the same computations as above, we have by the
Minkovski inequality
n
X

ak+1 Aγk+1 (θk )




b θ (Xk+1 ) − Pθ H
b θ (Xk+1 )
Pθk+1 H
k+1
k
k

k=0

≤C
Lq

n
X

ak+1 γk+1 .

k=0

6.4 Proof of Theorem 6
We write ηn+1 = Bn + (ηn+1 − Bn ) where Bn is given by (12). Observe that {ηn+1 −
Bn , n ∈ N} is a martingale-increment sequence. Sufficient conditions for the almostsure convergence of a martingale and the control of Lq -moments can be found in (Hall
and Heyde, 1980, Theorems 2.10 and 2.17). Then the proof follows from Proposition 5
and Lemma 14.

Appendix A. section 4
By using the Cauchy-Schwartz inequality, it holds
1/2

Z

Z
exp(`c (θ|u)) φ(u)du ≥
Z

exp(0.5`c (θ|u)) φ(u)du

2 Z
 Z

4
exp(`c (θ|u)) kuk φ(u) du ≤
exp(0.5`c (θ|u))φ(u) du
exp(3`c (θ|u)/2)kuk φ(u)du
2

28

which implies that
Z

exp(`c (θ|u))kuk2 φ(u)du
R
exp(`c (θ|v))φ(v)dv
Z
1/2
4
≤
exp (3`c (θ|u)/2) kuk φ(u) du
R

2

kuk πθ (u)du =

Since exp(`c (θ|u)) ≤ 1 (it is the likelihood of i.i.d. Bernoulli variables) and
q(2 + q), we have
Z
p
sup kuk2 πθ (u)du ≤ q(2 + q) .

R

kuk4 φ(u)du =

θ∈Θ

Appendix B. section 5
For θ, ϑ ∈ Θ, the (i, j)-th entry of the matrix ∇`(θ) − ∇`(ϑ) is given by
Z
(∇`(θ) − ∇`(ϑ))ij =

Z
B̄ij (x)πϑ (dx) −

B̄ij (x)πθ (dx).

Xp

Xp

For t ∈ [0, 1] let
def

πt (dz) = exp

B̄(z), tϑ + (1 − t)θ



Z
/

exp



B̄(x), tϑ + (1 − t)θ

µ(dx),

defines a probability measure on Xp . It is straightforward to check that
Z
(∇`(θ) − ∇`(ϑ))ij =
and that t 7→
d
dt

R

Z
B̄ij (x)π1 (dx) −

B̄ij (x)π0 (dx),

B̄ij (x)πt (dx) is differentiable with derivative


Z

Z
B̄ij (x)πt (dx) =

Z



B̄ij (x) B̄(x) −

= Covπt

B̄(z)πt (dz), ϑ − θ πt (dx),

B̄ij (X), B̄(X), ϑ − θ ,

where the covariance is taken assuming that X ∼ πt . Hence
Z

1

(∇`(θ) − ∇`(ϑ))ij =

dt Covt B̄ij (X), B̄(X), ϑ − θ
0
sX
≤ osc(B̄ij )
osc2 (B̄kl )kθ − ϑk2 .



k≤l

This implies the inequality (16).
Acknowledgments: We are grateful to George Michailidis for very helpful discussions. This work is partly supported by NSF grant DMS-1228164.
29

References
S. Allassonnière and E. Kuhn. Convergent Stochastic Expectation Maximization
algorithm with efficient sampling in high dimension. Application to deformable
template model estimation. Comput. Stat. Data An., 91:4–19, 2015.
C. Andrieu and E. Moulines. On the ergodicity properties of some adaptive MCMC
algorithms. Ann. Appl. Probab., 16(3):1462–1505, 2006.
C. Andrieu, E. Moulines, and P. Priouret. Stability of stochastic approximation under
verifiable conditions. SIAM J. Control Optim., 44(1):283–312, 2005.
O. Banerjee, L. El Ghaoui, and A. d’Aspremont. Model selection through sparse
maximum likelihood estimation for multivariate Gaussian or binary data. J. Mach.
Learn. Res., 9:485–516, 2008.
H. Bauschke and P.L. Combettes. Convex analysis and monotone operator theory
in Hilbert spaces. CMS Books in Mathematics/Ouvrages de Mathématiques de la
SMC. Springer, New York, 2011. ISBN 978-1-4419-9466-0. With a foreword by
Hédy Attouch.
A. Beck and M. Teboulle. Gradient-based algorithms with applications to signalrecovery problems. In Convex optimization in signal processing and communications, pages 42–88. Cambridge Univ. Press, Cambridge, 2010.
A. Benveniste, M. Métivier, and P. Priouret. Adaptive algorithms and stochastic
approximations, volume 22 of Applications of Mathematics (New York). SpringerVerlag, Berlin, 1990.
P. Biane, J. Pitman, and M. Yor. Probability laws related to the Jacobi theta and
Riemann zeta functions, and Brownian excursions. Bull. Amer. Math. Soc. (N.S.),
38(4):435–465 (electronic), 2001. ISSN 0273-0979.
H.M. Choi and J. P. Hobert. The polya-gamma gibbs sampler for bayesian logistic
regression is uniformly ergodic. Electronic Journal of Statistics, 7:2054–2064, 2013.
P.L. Combettes. Inherently parallel Algorithms in Feasibility and Optimization and
their Applications, chapter Quasi-Fejerian analysis of some optimization algorithms,
pages 115–152. Elsevier Science, 2001.
P.L. Combettes and J.C. Pesquet. Proximal splitting methods in signal processing. In
Fixed-point algorithms for inverse problems in science and engineering, volume 49
of Springer Optim. Appl., pages 185–212. Springer, New York, 2011.
P.L. Combettes and J.C. Pesquet. Stochastic Quasi-Fejer block-coordinate fixed point
iterations with random sweeping. SIAM J. Optim., 25(2):1221–1248, 2015a.
P.L. Combettes and J.C. Pesquet. Stochastic Approximations and Perturbations in Forward-Backward Splitting for Monotone Operators. Technical report,
arXiv:1507.07095v1, 2015b.
30

P.L. Combettes and V. Wajs. Signal recovery by proximal forward-backward splitting.
Multiscale Modeling and Simulation, 4(4):1168–1200, 2005.
A. Cotter, O. Shamir, N. Srebro, and K. Sridharan. Better mini-batch algorithms via
accelerated gradient methods. In J. Shawe-taylor, R.s. Zemel, P. Bartlett, F.c.n.
Pereira, and K.q. Weinberger, editors, Advances in Neural Information Processing
Systems 24, pages 1647–1655. 2011.
J. Duchi, E. Hazan, and Y. Singer. Adaptive subgradient methods for online learning
and stochastic optimization. J. Mach. Learn. Res., 12:2121–2159, 2011. ISSN
1532-4435.
M. Ekeberg, C. Lövkvist, Y. Lan, M. Weigt, and E. Aurell. Improved contact prediction in proteins: Using pseudolikelihoods to infer potts models. Phys. Rev. E, 87:
012707, 2013.
G. Fort and E. Moulines. Convergence of the Monte Carlo expectation maximization
for curved exponential families. Ann. Statist., 31(4):1220–1259, 2003. ISSN 00905364.
G. Fort, E. Moulines, and P. Priouret. Convergence of adaptive and interacting
Markov chain Monte Carlo algorithms. Ann. Statist., 39(6):3262–3289, 2011. ISSN
0090-5364.
G. Fort, E. Moulines, M. Vihola, and A. Schreck. Convergence of Markovian Stochastic Approximation with discontinuous dynamics. Technical report, arXiv 
1403.6803, 2014.
G. Fort, B. Jourdain, E. Kuhn, T. Lelièvre, and G. Stoltz. Convergence of the WangLandau algorithm. Mathematics of Computation, 84:2297–2327, 2015.
C.J. Geyer. On the convergence of Monte Carlo maximum likelihood calculations. J.
Roy. Statist. Soc. Ser. B, 56(1):261–274, 1994.
J. Guo, E. Levina, G. Michailidis, and J. Zhu. Joint structure estimation for categorical Markov networks. Technical report, Univ. of Michigan, 2010.
P. Hall and C.C. Heyde. Martingale Limit Theory and its Application. Academic
Press, 1980.
H. Höfling and R. Tibshirani. Estimation of sparse binary pairwise Markov networks
using pseudo-likelihoods. J. Mach. Learn. Res., 10:883–906, 2009.
C. Hu, W. Pan, and J.T. Kwok. Accelerated gradient methods for stochastic optimization and online learning. In Y. Bengio, D. Schuurmans, J. Lafferty, C. K. I
Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems, pages 781–789, 2009.
A. Juditsky and A. Nemirovski. First-order methods for nonsmooth convex large-scale
optimization, i: General purpose methods. In S. Sra, S. Nowozin, and S. Wright,
editors, Oxford Handbook of Innovation, pages 121–146. MIT Press, Boston, 2012a.
31

A. Juditsky and A. Nemirovski. First-order methods for nonsmooth convex largescale optimization, ii: Utilizing problem’s structure. In S. Sra, S. Nowozin, and
S. Wright, editors, Oxford Handbook of Innovation, pages 149–181. MIT Press,
Boston, 2012b.
H. Kamisetty, S. Ovchinnikov, and D. Baker. Assessing the utility of coevolutionbased residue-residue contact predictions in a sequence-and structure-rich era. Proceedings of the National Academy of Sciences, 2013. doi: 10.1073/pnas.1314045110.
G. Lan. An optimal method for stochastic composite optimization. Math. Program.,
133(1-2, Ser. A):365–397, 2012. ISSN 0025-5610.
J. Lin, L. Rosasco, S. Villa, and D.X. Zhou. Modified Fejer Sequences and Applications. Technical report, arXiv:1510:04641v1 math.OC, 2015.
G.J. McLachlan and T. Krishnan. The EM algorithms and Extensions. WileyInterscience; 2 edition, 2008.
S. Meyn and R.L. Tweedie. Markov chains and stochastic stability. Cambridge University Press, Cambridge, second edition, 2009. ISBN 978-0-521-73182-9. With a
prologue by Peter W. Glynn.
J. Moreau. Fonctions convexes duales et points proximaux dans un espace hilbertien.
CR Acad. Sci. Paris Sér. A Math, 255:2897–2899, 1962.
A. Nemirovski, A. Juditsky, G. Lan, and A. Shapiro. Robust stochastic approximation
approach to stochastic programming. SIAM J. Optim., 19(4):1574–1609, 2008.
ISSN 1052-6234.
Y.E. Nesterov. Introductory Lectures on Convex Optimization, A basic course. Kluwer
Academic Publishers, 2004.
A. Nitanda. Stochastic proximal gradient descent with acceleration techniques. NIPS,
2014.
H.E. Ogden. A sequential reduction method for inference in generalized linear mixed
models. Electron. J. Statist., 9(1):135–152, 2015.
N. Parikh and S. Boyd. Proximal algorithms. Foundations and Trends in Optimization, 1(3):123–231, 2013.
N. G. Polson, J. G. Scott, and J. Windle. Bayesian inference for logistic models using
Polya-Gamma latent variables. J. Am. Stat. Assoc., 108(504):1339–1349, 2013.
B.T. Polyak. Introduction to Optimization. xx, 1987.
P. Ravikumar, M.J. Wainwright, and J.D. Lafferty. High-dimensional Ising model
selection using `1 -regularized logistic regression. Ann. Statist., 38(3):1287–1319,
2010.
32

C.P. Robert and G. Casella. Monte Carlo Statistical Methods. Springer Texts in
Statistics. Springer; 2nd edition, 2005.
L. Rosasco, S. Villa, and B.C. Vu. Convergence of a Stochastic Proximal Gradient
Algorithm. Technical report, arXiv:1403.5075v3, 2014.
L. Rosasco, S. Villa, and B.C. Vu. A Stochastic Inertial Forward-Backward
Splitting Algorithm for multi-variate monotone inclusions. Technical report,
arXiv:1507.00848v1, 2015.
E. Saksman and M. Vihola. On the ergodicity of the adaptive Metropolis algorithm
on unbounded domains. Ann. Appl. Probab., 20(6):2178–2203, 2010.
J. Schelldorfer, L. Meier, and P. Bühlmann. GLMMLasso: an algorithm for highdimensional generalized linear mixed models using `1 -penalization. J. Comput.
Graph. Statist., 23(2):460–477, 2014.
M.W. Schmidt, N. Le Roux, and F. Bach. Convergence rates of inexact proximalgradient methods for convex optimization. In NIPS, pages 1458–1466, 2011. see
also the technical report INRIA-00618152.
A. Schreck, G. Fort, and E. Moulines. Adaptive Equi-energy sampler : convergence and illustration. ACM Transactions on Modeling and Computer Simulation
(TOMACS), 23(1):Art 5., 2013.
J. Shao. Mathematical Statistics. Springer texts in Statistics, 2003.
L. Xiao. Dual averaging methods for regularized stochastic learning and online optimization. J. Mach. Learn. Res., 11:2543–2596, 2010. ISSN 1532-4435.
L. Xiao and T. Zhang. A Proximal Stochastic Gradient Method with Progressive
Variance Reduction. SIAM J. Optim., 24:2057–2075, 2014.
L. Xue, H. Zou, and T. Cai. Non-concave penalized composite likelihood estimation
of sparse ising models. Ann. Statist., 40:1403–1429, 2012.

33

