arXiv:1209.5654v2 [math.PR] 30 Sep 2016

Bernoulli 23(1), 2017, 670–709
DOI: 10.3150/14-BEJ680

Nonasymptotic analysis of adaptive and
annealed Feynman–Kac particle models
FRANÇOIS GIRAUD1,2,* and PIERRE DEL MORAL2,**
1

CEA-CESTA, 33114 Le Barp, France.
INRIA Bordeaux Sud-Ouest, team ALEA, Domaine Universitaire, 351, cours de la Libération,
33405 Talence Cedex, France. E-mail: * francois.giraud@ens-cachan.org
3
School of Maths and Stats. University of New South Wales, Sydney, Australia.
E-mail: ** pierre.del-moral@inria.fr
2

Sequential and quantum Monte Carlo methods, as well as genetic type search algorithms can
be interpreted as a mean field and interacting particle approximations of Feynman–Kac models
in distribution spaces. The performance of these population Monte Carlo algorithms is strongly
related to the stability properties of nonlinear Feynman–Kac semigroups. In this paper, we analyze these models in terms of Dobrushin ergodic coefficients of the reference Markov transitions
and the oscillations of the potential functions. Sufficient conditions for uniform concentration
inequalities w.r.t. time are expressed explicitly in terms of these two quantities. We provide
an original perturbation analysis that applies to annealed and adaptive Feynman–Kac models,
yielding what seems to be the first results of this kind for these types of models. Special attention is devoted to the particular case of Boltzmann–Gibbs measures’ sampling. In this context,
we design an explicit way of tuning the number of Markov chain Monte Carlo iterations with
temperature schedule. We also design an alternative interacting particle method based on an
adaptive strategy to define the temperature increments. The theoretical analysis of the performance of this adaptive model is much more involved as both the potential functions and the
reference Markov transitions now depend on the random evolution on the particle model. The
nonasymptotic analysis of these complex adaptive models is an open research problem. We initiate this study with the concentration analysis of a simplified adaptive models based on reference
Markov transitions that coincide with the limiting quantities, as the number of particles tends
to infinity.
Keywords: adaptive models; Feynman–Kac models; interacting particle systems; sequential
Monte Carlo methods

Introduction
Feynman–Kac (abbreviate FK ) particle methods, also called sequential, quantum or diffusion Monte Carlo methods, are stochastic algorithms to sample from a sequence of

This is an electronic reprint of the original article published by the ISI/BS in Bernoulli,
2017, Vol. 23, No. 1, 670–709. This reprint differs from the original in pagination and
typographic detail.
1350-7265

c

2017 ISI/BS

2

F. Giraud and P. Del Moral

complex high-dimensional probability distributions. These stochastic simulation techniques are of current use in numerical physics [1, 5, 33] to compute ground state energies
in molecular systems. They are also used in statistics, signal processing and information
sciences [6, 15, 17, 20] to compute posterior distributions of a partially observed signal
or unknown parameters. In the evolutionary computing literature, these Monte Carlo
methods are used as natural population search algorithms for solving optimization problems. From the pure mathematical viewpoint, these advanced Monte Carlo methods are
an interacting particle system (abbreviate IPS ) interpretation of FK models. For a more
thorough discussion on these models, we refer the reader to the monograph [16], and the
references therein. The principle (see also [17] and the references therein) is to approximate a sequence of target probability distributions (ηn )n by a large cloud of random
samples termed particles or walkers. The algorithm starts with N independent samples
from η0 and then alternates two types of steps: an acceptance-rejection scheme equipped
with a selection type recycling mechanism, and a sequence of free exploration of the state
space.
In the recycling stage, the current cloud of particles is transformed by randomly duplicating and eliminating particles in a suitable way, similarly to a selection step in models
of population genetics. In the Markov evolution step, particles move independently one
each other (mutation step).
This method is often used for solving sequential problems, such as filtering (see, e.g.,
[6, 15, 21, 28]). In other interesting problems, these algorithms also turn out to be efficient
to sample from a single target measure η. In this context, the central idea is to find a judicious interpolating sequence of measures (ηk )0≤k≤n with increasing sampling complexity,
starting from some initial distribution η0 , up to the terminal one ηn = η. Consecutive
measures ηk and ηk+1 are sufficiently similar to allow for efficient importance sampling
and/or acceptance-rejection sampling. The sequential aspect of the approach is then an
“artificial way” to introduce the difficulty of sampling gradually. In this vein, important
examples are provided by annealed models. More generally, a crucial point is that large
population sizes allow to cover several modes simultaneously. This is an advantage compared to standard MCMC methods that are more likely to be trapped in local modes.
These sequential samplers have been used with success in several application domains,
including rare events simulation (see [9]), stochastic optimization and more generally
Boltzmann–Gibbs measures sampling [17].
Up to now, IPS algorithms have been mostly analyzed using asymptotic (i.e., when
number of particles N tends to infinity) techniques, notably through fluctuation theorems
and large deviation principles (see, e.g., [6, 11, 14, 15, 19, 20, 23, 25, 35] and [16] for an
overview).
Some nonasymptotic theorems have been recently developed [10, 18], but unfortunately none of them apply to analyze annealed and adaptive FK particle models. On
the other hand, these type of nonhomogeneous IPS algorithms are of current use for
solving concrete problems arising in numerical physics and engineering sciences (see,
e.g., [4, 13, 27, 30, 34, 37–39]). By the lack of nonasymptotic estimates, these particle
algorithms are used as natural heuristics.

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

3

The main contribution of this article is to analyze these two classes of time nonhomogeneous IPS models. Our approach is based on semigroup techniques and on an original
perturbation analysis to derive several uniform estimates w.r.t. the time parameter.
More precisely, in the case of annealed type models, we estimate explicitly the stability properties of FK semigroup in terms of the Dobrushin ergodic coefficient of the
reference Markov chain and the oscillations of the potential functions. We combine these
techniques with nonasymptotic theorems on Lp -mean error bounds [25] and some useful
concentration inequalities [22, 26]. Then we provide parameter tuning strategies that allow to deduce some useful uniform concentration inequalities w.r.t. the time parameter.
These results apply to nonhomogeneous FK models associated with cooling temperature
parameters. In this situation, the sequence of measures ηn is associated with a nonincreasing temperature parameter. We mention that other independent approaches, such
as Whiteley’s [41] or Schweizer’s [40], are based on, for example, drift conditions, hyperboundedness, spectral gaps, or nonasymptotic biais and variance decompositions. These
approaches lead to convergence results that may also apply to noncompact state spaces.
To our knowledge, these techniques are restricted to nonasymptotic variance theorems
and they cannot be used to derive uniform and exponential concentration inequalities.
It seems also difficult to extend these approaches to analyze the adaptive IPS model
discussed in the present article. To solve these questions, we develop a perturbation technique of stochastic FK semigroups. In contrast to traditional FK semigroup, the adaptive
particle scheme is now based on random potential functions that depend on a cooling
schedule adapted to the variability and the adaptation of the random populations.
The rest of the article is organized as follows. In a preliminary section, we recall a few
essential notions related to Dobrushin coefficients or FK semigroups. We also provide
some important nonasymptotic results we use in the further development of the article.
Section 2 is concerned with the semigroup stability analysis of these models. We also
provide a couple of uniform Lp -deviations and concentration estimates. In Section 3, we
apply these results to Boltzmann–Gibbs models associated with a decreasing temperature
schedule. In this context, IPS algorithm can be interpreted as a sequence of interacting
simulated annealing algorithms (abbreviate ISA). We design an explicit way of tuning the
number of Markov chain Monte Carlo iterations with the temperature schedule. Finally,
in Section 4, we propose an alternative ISA method based on an original adaptive strategy
to design on the flow the temperature decrements. We provide a nonasymptotic study for
a simplified model, based on a perturbation analysis. We end the article with Lp -deviation
estimates as well as a couple of concentration inequalities.

Statement of some results
Feynman–Kac particle algorithms consist in evolving a particle system ζn = (ζn1 , . . . , ζnN )
of size N , on a given state space E. Their interacting evolution is decomposed into two
genetic type transitions: a selection step, associated with some positive potential function
Gn ; and a mutation step, where the selected particles evolve randomly according to a
given Markov transition Mn (a more detailed description of these IPS algorithms is

4

F. Giraud and P. Del Moral

P
provided in Section 1.4). In this context, the occupation measures ηnN := N1 1≤i≤N δζni
are N -approximations of a sequence of measures ηn defined by the FK recursive formulae:
ηn (f ) =

ηn−1 (Gn × Mn .f )
,
ηn−1 (Gn )

for all bounded measurable function f on E (a more detailed discussion on these evolution
equations is provided in Section 1.3.1).
To describe with some precision the main results of the article, we consider the pair of
parameters (gn , bn ) defined below.
Gn (x)
x,y∈E Gn (y)

gn := sup

and bn = β(Mn ) := sup |Mn (x, A) − Mn (y, A)|.
x,y∈E
A⊂E

The quantity β(Mn ) is called the Dobrushin ergodic coefficient of the Markov transition
Mn . One of our first main results can be basically stated as follows.
Theorem 1. We assume that
sup gp ≤ M

and

p≥1

sup bp ≤
p≥1

a
M (1 + a)

for some finite constant M < ∞ and some a ∈ (0, 1). In this situation, for any n ≥ 0,
N ≥ 1, y ≥ 0 and f ∈ O1 (E), the probability of the event
ηnN (f ) − ηn (f ) ≤

√
r1⋆ (1 + h0 (y)) + r2⋆ N y
N

is greater than 1 − e−y , where r1⋆ and r2⋆ are some constants that are explicitly defined in
√
terms of (a, M ), and h0 (y) = 2y + 2 y.
In Section 2.2, under the same assumptions of Theorem 1, we also prove uniform Lp mean error bounds as well as new concentration inequalities for unnormalized particle
models. We also extend the analysis to the situation where gn −→ 1.
n→+∞

We already mentioned that the regularity conditions on bn may appear difficult to check
since the Markov kernels are often dictated by the application under study. However,
we can deal with this problem as soon as we can simulate a Markov kernel Kn such
that ηn .Kn = ηn . Indeed, to stabilize the system, the designer can “add” several MCMC
evolution steps next to each Mn -mutation step. From a more formal viewpoint, the target
sequence (ηn )n is clearly also solution of the FK measure-valued equations associated
with the Markov kernels Mn′ = Mn .Knmn , where iteration numbers mn are to be chosen
loosely. This system is more stable since the corresponding b′n satisfy
b′n = β(Mn′ ) ≤ bn .β(Knmn ) ≤ bn .β(Kn )mn .

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

5

In such cases, Theorem 1 and its extension provide sufficient conditions on the iteration
numbers mn to ensure the convergence and the stability properties of the algorithm.
These results apply to stochastic optimization problems. Let V : E → R be a bounded
potential function, βn a sequence which tends to infinity, and m a reference measure on
E. It is well known that the sequence of Boltzmann–Gibbs measures
ηn (dx) ∝ e−βn .V (x) m(dx)
concentrates on V ’s global minima (in the sense of m-essinf(V )). In the above display, ∝ stands for the proportional sign. One central observation is that these measures can be interpreted as a FK flow of measures associated with potential functions
Gn = e−(βn −βn−1 ).V and Markov kernels Mn = Kβmnn k0 where Kβn is a simulating annealing kernel (see Section 3.2) and mn and k0 are given iteration parameters. In the further
development of this section, we let K be the proposal transition of the simulated annealing transition Kβ . In this context, the IPS methods can be used to minimize V . The
conditions on bn and gn can be turned into conditions on the temperature schedule βn
and the number of MCMC iterations mn . Moreover, combining our results with standard
concentration properties of Boltzmann–Gibbs measures, we derive some convergence results in terms of optimization performance. In this notation, our second main result is
basically stated as follows.
Theorem 2. Let us fix a ∈ (0, 1). We assume that for any x ∈ E, K k0 (x, ·) ≥ δν(·) for
some measure ν on E, some δ > 0 and some k0 ≥ 1. We also assume that the temperature
increments ∆p := βp −βp−1 and the iteration numbers mp satisfy the following conditions:
sup ∆p ≤ ∆
p≥1

and

mp ≥

log(e∆. osc(V ) (1 + a)/a)eosc(V ).βp
δ

i
for some constant ∆. For all ε > 0, let pN
n (ε) be the proportion of particles (ζn ) so that
V (ζni ) ≥ Vmin + ε. Then, for any n ≥ 0, N ≥ 1, y ≥ 0 and for all ε′ < ε, the probability of
the event
√
′
e−βn (ε−ε ) r1⋆ (1 + h0 (y)) + r2⋆ N y
pN
(ε)
≤
+
n
mε′
N
√
′
−y
is greater than 1 − e , with mε′ = m(V ≤ Vmin + ε ), h0 (y) = 2y + 2 y, and the same
constants (r1⋆ , r2⋆ ) as the ones stated in Theorem 1 (with M = e∆. osc(V ) ).

It is instructive to compare the estimates in the above theorem with the performance
analysis of the traditional simulated annealing model (abbreviate SA). Firstly, most of the
literature on SA models is concerned with the weak convergence of the law of the random
states of the algorithm. When the initial temperature of the scheme is greater than some
critical value, using a logarithmic cooling schedule, it is well known that the probability
for the random state to be in the global extrema levels tends to 1, as the time parameter
tends to ∞. The cooling schedule presented in Theorem 2 is again a logarithmic one.

6

F. Giraud and P. Del Moral

In contrast to the SA model (see, e.g., Theorem 1 in [32], and Theorem 4.3.16 in [3]),
Theorem 2 allows to quantify the performance analysis of the ISA model in terms of
uniform concentration inequalities that does not depend on a critical parameter.
In practice, choosing the sequence of increments ∆n = (βn − βn−1 ) in advance can cause
computational problems. To solve this problem, adaptive strategies, where increment ∆n
depends on the current set of particles ζn−1 , are of common use in the engineering
community (see, e.g., [13, 27, 34, 37, 39]). In this context, we propose to study the case
where the increment ∆N
n is chosen so that
N

N
ηn−1
(e−∆n ·V ) = ε,

where ε > 0 is a given constant (see Section 4.2 for a detailed description of the algorithm). Computationally speaking, ε is the expectation of the proportion of particles
which are not concerned with the recycling mechanism in the selection step. We interpret this particle process as a perturbation of a theoretical FK sequence ηn associated
with a theoretical temperature schedule βn . As we mentioned in the abstract and the
Introduction, the theoretical analysis of this class of adaptive particle model is much
more involved as the potential functions and the mutations transitions now depend on
the evolution of the particle populations. We consider a simplified adaptive model associated with the limiting reference Markov transitions. A precise description of these
adaptive particle algorithms, and their reduced simplified versions are provided in Sections 4.1, 4.2 and 4.3.
Our main result is the following Lp -mean error estimate, where (ηnN ) denotes the
sequence of empirical measures associated with the particle system described in (4.3)
page 30.
Theorem 3. For any p ≥ 1, n ≥ 0, N ≥ 1 and any bounded by 1 function f , we have
p 1/p

E(|ηnN (f ) − ηn (f )| )
with cn =

Vmax e∆n Vmax
ε·ηn−1 (V )

n
n
Bp X Y
≤√
(bi gi (1 + ci )),
N k=0 i=k+1

, ∆n = βn − βn−1 and Bp defined below.
2p
B2p
=

(2p)!
;
2p .p!

2p+1
B2p+1
=

(2p + 1)!
√
.
2p + 1

2p .p!

(0.1)

Under appropriate regularity conditions on the parameters bn , gn , cn , we mention that
these Lp -mean error bounds also provide uniform concentration inequalities.
The proofs of Theorems 1, 2, 3 and related uniform exponential estimates are detailed,
respectively, in Sections 2.2, 3.2 and 4.4.

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

7

1. Some preliminaries
1.1. Basic notation
Let (E, r) be a complete, separable metric space and let E be the σ-algebra of Borel
subsets of E. Denote by P(E) the space of probability measures on E. Let B(E) be the
space of bounded, measurable, real-valued functions on E.
If µ ∈ P(E),
f ∈ B(E) and K, K1 , K2 are Markov kernels on E, then µ(f ) denotes the
R
quantity E f (x)µ(dx), K1 .K2 denotes the Markov kernel defined by
(K1 .K2 )(x, A) =

Z

K1 (x, dy)K2 (y, A),

E

K.f denotes the function defined by
K.f (x) =

Z

K(x, dy)f (y)

E

and µ.K denotes the probability measure defined by
Z
µ.K(A) =
K(x, A)µ(dx).
E

If G is a positive, bounded function on E, then ψG : P(E) → P(E) denotes the
Boltzmann–Gibbs transformation associated with G, defined by
∀µ ∈ P(E), ∀f ∈ B(E)

ψG (µ)(f ) =

µ(G × f )
.
µ(G)

For any f ∈ B(E), let kf k∞ = supx∈E |f (x)| and osc(f ) = (fmax − fmin ). Let B1 (E) ⊂
B(E) be the subset of functions f so that kf k∞ ≤ 1, and O1 (E) ⊂ B(E) be the subset
of functions f so that osc(f ) ≤ 1. For any random variable X : Ω → R defined on some
probability space (Ω, F , P), and any p ≥ 1, kXkp stands for the Lp norm E(|X|p )1/p . Let
PΩ (E) be the set of random probability measures on E. For all p ≥ 1, we denote by dp
the distance on PΩ (E) defined for all random measures µ̂, ν̂ by
dp (µ̂, ν̂) =

sup kµ̂(f ) − ν̂(f )kp .

f ∈O1 (E)

Finally, for any x ∈ E, δx stands for the Dirac measure centered on x.

1.2. Dobrushin Ergodic coefficient
Let us recall here the definitions as well as some simple properties that will be useful in
the following.

8

F. Giraud and P. Del Moral

Definition 4. Let µ, ν ∈ P(E). The total variation distance between µ and ν is defined
by
kµ − νktv = sup{|µ(A) − ν(A)|; A ∈ E},
or in an equivalent way
kµ − νktv = sup{|µ(f ) − ν(f )|; f ∈ B1 (E), f ≥ 0}.
Definition 5. To each Markov kernel K on E, is associated its Dobrushin ergodic coefficient β(K) ∈ [0, 1] defined by
β(K) = sup{K(x, A) − K(y, A); x, y ∈ E, A ∈ E},
or in an equivalent way


kµ.K − ν.Kktv
; µ, ν ∈ P(E), µ 6= ν .
β(K) = sup
kµ − νktv
The parameter β(K) characterizes mixing properties of the Markov kernel K. Note
that function β is an operator norm, in the sense that β(K1 .K2 ) ≤ β(K1 ).β(K2 ), for any
couple of Markov kernels K1 , K2 . By definition, for any measures µ, ν ∈ P(E) and any
Markov kernel K, we have kµ.K − ν.Kktv ≤ β(K).kµ − νktv . Otherwise, for any function
f ∈ B(E),
osc(K.f ) ≤ β(K) · osc(f ).

(1.1)

Further details on these ergodic coefficients can be found in the monograph [16]. Let us
give here a lemma that we will need hereinafter (a proof is given in the Appendix, page
37).
Lemma 6. Let µ, ν ∈ P(E) and G a positive, bounded function on E satisfying
supx,y∈E G(x)
G(y) ≤ g, for some finite constant g ≥ 0. In this situation, we have
kΨG (µ) − ΨG (ν)ktv ≤ g.kµ − νktv .

1.3. Feynman–Kac models
We recall here some standard tools related to FK models. They provide useful theoretical
background and notation to formalize and analyze IPS methods (see, e.g., [20, 22, 25]
for further details).
1.3.1. Evolution equations
Consider a sequence of probability measures (ηn )n , defined by an initial measure η0 and
recursive relations
ηn−1 (Gn × Mn .f )
(1.2)
∀f ∈ B(E)
ηn (f ) =
ηn−1 (Gn )

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

9

for positive functions Gn ∈ B(E) and Markov kernels Mn with Mn (x, ·) ∈ P(E) and
Mn (·, A) ∈ B1 (E). This is the sequence of measures we mainly wish to approximate with
the IPS algorithm. In an equivalent way, (ηn )n can be defined by the relation
ηn = φn (ηn−1 ),
where φn : P(E) → P(E) is the FK transformation associated with potential function Gn
and Markov kernel Mn and defined by
φn (ηn−1 ) = ψGn (ηn−1 ).Mn
with
ψGn (ηn−1 )(dx) :=

1
Gn (x)ηn−1 (dx).
ηn−1 (Gn )

The next formula provides an interpretation of the Boltzmann–Gibbs transformation in
terms of a nonlinear Markov transport equation
Z
ΨGn (ηn−1 )(dy) = (ηn−1 Sn,ηn−1 )(dy) := ηn−1 (dx)Sn,ηn−1 (x, dy)
with the Markov transition Sn,ηn defined below
Sn,ηn−1 (x, dy) = εn .Gn (x)δx (dy) + (1 − εn .Gn (x))ΨGn (ηn−1 )(dy)
(for any constant εn > 0 so that εn .Gn ≤ 1). This implies
ηn = ηn−1 Kn,ηn−1

with Kn,ηn−1 = Sn,ηn−1 Mn .

(1.3)

Therefore, ηn can be interpreted as the distributions of the random states X n of a Markov
chain whose Markov transitions
P(X n+1 ∈ dy|X n = x) := Kn+1,ηn (x, dy)

(1.4)

depend on the current distribution ηn = Law(X n ).
We finally recall that the measures ηn admit the following functional representations:
ηn (f ) =

γn (f )
γn (1)

(1 stands for the unit function) with the unnormalized FK measures γn defined by the
formulae
γ0 = η0 ;

γn (f ) = γn−1 (Gn × Mn .f ).

(1.5)

Comparing this definition with (1.2), it is clear that the normalizing constant γn (1)
satisfies
n
Y
γn (1) =
ηp−1 (Gp ).
(1.6)
p=1

10

F. Giraud and P. Del Moral

The special interest given to this quantity will be motivated in Section 3.1.
1.3.2. Feynman–Kac semigroup
An important point is that the semigroup transformations
φp,n := φn ◦ φn−1 ◦ · · · ◦ φp+1
admit a comparable structure as each of the φk . To be more precise, for each integer p,
let us define the unnormalized integral operator Qp
∀f ∈ B(E)

Qp .f = Gp .Mp .f

(1.7)

and the composition operators Qp,n defined by the backward recursion
Qp,n = Qp+1 .(Qp+2 · · · Qn ) = Qp+1 .Qp+1,n .

(1.8)

We use the convention Qn,n = Id for p = n. Comparing these definitions with (1.5), it is
clear that γn = γn−1 .Qn and more generally
γn = γp .Qp,n .
for any p ≤ n. The semigroup φp,n can be expressed in terms of Qp,n with the following
formulae:
φp,n (µ)(f ) =

µ(Qp,n .f )
µ(Qp,n .1)

for any f ∈ B(E) and µ ∈ P(E). Finally, if we set
Pp,n .f =

Qp,n .f
Qp,n .1

and Gp,n = Qp,n .1

then we find that
φp,n (µ)(f ) =

µ(Gp,n .Pp,n .f )
,
µ(Gp,n )

or in other words: φp,n (µ) = ψGp,n (µ).Pp,n .

1.4. The interacting particle system model
The central idea is to approximate the measures ηn by simulating an interacting particle
system (ζn )n = (ζn1 , . . . , ζnN )n of size N so that
ηnN =

1 X
δζni −→ ηn .
N ↑∞
N
1≤i≤N

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

11

Of course, the main issue is to make precise and to quantify this convergence. The particle
model is defined as follows.
We start with N independent samples ζ0 = (ζ01 , . . . , ζ0N ) from η0 . Then the particle
dynamics alternates two genetic type transitions.
During the first step, every particle ζni evolves to a new particle ζbni randomly chosen
with the distribution
Sn+1,ηnN (ζni , dx) := εn+1 .Gn+1 (ζni )δζni (dx) + (1 − εn+1 .Gn+1 (ζni ))ΨGn+1 (ηnN )(dx)

with the updated measures
ΨGn+1 (ηnN ) =

N
X
j=1

Gn+1 (ζnj )
δ j.
PN
k ζn
k=1 Gn+1 (ζn )

This transition can be interpreted as an acceptance-rejection scheme with a recycling
mechanism. In the second step, the selected particles ζbni evolve randomly according to
the Markov transitions Mn+1 . In other words, for any 1 ≤ i ≤ N , we sample a random
i
state ζn+1
with distribution Mn+1 (ζbni , dx).
In view of (1.4), if we replace ηnN by ηn , then ζn coincide with N independent copies of
the Markov chains X n defined in (1.3). On the other hand, by the law of large numbers,
we have η0N ≃ η0 so that
η1N ≃ η0N .K1,η0N ≃ η0 .K1,η0 = η1 .
Iterating this approximation procedure, the empirical measure ηnN is expected to approximate ηn at any time n ≥ 0. As for the unnormalized measures γn , we define
γnN (1) =

n
Y

N
ηp−1
(Gp )

p=1

Qn
N
(Gp ). Let us
(mimicking formula (1.6)) and more generally γnN (f ) = ηnN (f ) × p=1 ηp−1
mention (see, e.g., [15]) that these particle models provide an unbiased estimate of the
unnormalized measures; that is we have that
∀f ∈ B(E)

E(γnN (f )) = γn (f ).

In addition to the analysis of ηnN ’s convergence, the concentration properties of the unbiased estimators γnN (1) around their limiting values γn (1) will also be considered thereafter.

1.5. Some nonasymptotic results
To quantify the FK semigroup stability properties, it is convenient to introduce the
following parameters.

12

F. Giraud and P. Del Moral

Definition 7. For any integers p < n, we set
bn := β(Mn )

and

Gn (x)
x,y∈E Gn (y)

gn := sup

bp,n := β(Pp,n ),
and

Gp,n (x)
.
x,y∈E Gp,n (y)

gp,n := sup

The quantities gp,n , and, respectively, bp,n , reflect the oscillations of the potential
functions Gp,n , and, respectively, the mixing properties of the Markov transition Pp,n
associated with the FK semigroup φp,n described in Section 1.3.2. Several contraction
inequalities of φp,n w.r.t. the total variation norm or different types of relative entropies
can be derived in terms of these two quantities (see, for instance, [16]).
Let
√
h0 := x 7→ 2(x + x) and
(1.9)
x √
h1 := x 7→ + 2x.
3
The performance analysis developed in Sections 2.1 and 3.2 is partly based on the three
nonasymptotic inequalities presented below.
Firstly, the following Lp -mean error bound for all f ∈ B1 (E) is inspired by [25] (see
Section 2.2.2):
p 1/p
E(|ηnN (f ) − ηn (f )| )

n
Bp X
gk,n bk,n ,
≤√
N k=0

(1.10)

where Bp are the constants introduced in (0.1), page 6. A detailed proof of this estimate
is given in the Appendix, page 38.
Secondly, the following concentration inequality is derived in [26] (see Theorem 1.2
and Section A3). For all f ∈ O1 (E) and any y ≥ 0, we have
r 

rn
2y
N
≤ e−y ,
(1.11)
P ηn (f ) − ηn (f ) ≥ (1 + h0 (y)) + β n
N
N
where rn and β n are constants so that

n
X

3

gp,n
bp,n ,
r
≤
4

n


p=0

n
X

2

2


gp,n
b2p,n .
 βn ≤ 4
p=0

To be more precise with the derivation of these estimates, the upper bound on rn is
based on the quadratic remainder estimates provided in the end of page 14 and on page
28 (line 17) in [26]. The function h0 coincides with the function ǫ−1
0 used in [26], and we

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

13

2

have used the estimate provided on page 16 (line 13). The upper bound on β n is a direct
consequence of the estimates provided on page 9 (line 2) and on page 28 (line 13) in [26].
Thirdly, the following concentration inequality for unnormalized particle models γnN is
provided in [22] (see Theorem 6.13).
∀ǫ ∈ {+1, −1} and ∀y ≥ 0:




 N
r̄(n)
y
ǫ
γ (1)
P
≥
≤ e−y ,
(1.12)
log n
h0 (y) + τn⋆ σ̄n2 h1
n
γn (1)
N
N.σ̄n2
where quantities τn⋆ , σ̄n2 and r̄(n) can be estimated this way:
• τn⋆ = sup0≤q≤n τq,n , where τq,n satisfy
n−1

τq,n ≤

4X
gq,p .gp+1 .bq,p ;
n p=q

Pn−1
τ
2
where σq satisfy σq ≤ 1;
• σ̄n2 = q=0 σq2 ( τq,n
⋆ )
n
• r̄(n) satisfy
8 X
3
gp+1 .gq,p
.bq,p .
r̄(n) ≤
n

(1.13)

(1.14)

0≤q≤p<n

2. Nonasymptotic theorems
The formulae (1.10), (1.11) and (1.12) provide explicit nonasymptotic estimates in terms
of the quantities gp,n and bp,n . Written this way, they hardly apply to any IPS parameters
tuning decision, since the only known or calculable objects are generally the reference
Markov chain Mp and the elementary potential functions Gp . We thus have to estimate
gp,n and bp,n with some precision in terms of the gp and bp . This task is performed
in Section 2.1. In the second section, Section 2.2, we combine these estimates with the
concentration results presented in Section 1.5 to derive some useful uniform estimates
w.r.t. the time parameter.

2.1. Semigroup estimates
We start with a series of technical lemmas.
Lemma 8. Let K be a Markov kernel and G a positive function on E satisfying
supx,y∈E G(x)
G(y) ≤ g, for some finite constant g. In this situation, we have that
K.G(x)
≤ 1 + β(K)(g − 1).
x,y∈E K.G(y)
sup

14

F. Giraud and P. Del Moral

Proof. Let x, y ∈ E be s.t. K.G(x) ≥ K.G(y). Let us write
β(K)(Gmax − Gmin )
K.G(x) K.G(x) − K.G(y)
+ 1.
=
+1≤
K.G(y)
K.G(y)
Gmin
We check the last inequality using the fact that
K.G(x) − K.G(y) ≤ osc(K.G) ≤ β(K). osc(G) = β(K).(Gmax − Gmin ).
R
On the other hand, we have K.G(y) = u G(u)K(y, du) ≥ Gmin . The desired result is
−Gmin )
now obtained taking the supremum over all (x, y) ∈ E 2 . Note that β(K)(GGmax
+1
min
is exactly equal to 1 + β(K)(g − 1).
This ends the proof of the lemma.

Lemma 9. Let M be a Markov kernel, Q a not necessarily normalized integral operator
satisfying supx,y∈E Q.1(x)
Q.1(y) ≤ g, for some finite constant g ≥ 1 and f a bounded, nonnegative function. In this situation, the Markov kernel P defined by
P.f (x) :=

M.Q.f (x)
M.Q.1(x)

satisfies the following property:
β(P ) ≤ g.β(M ).β(P ′ ).
In the above display formula, P ′ is the Markov transition defined by P ′ .f (x) :=

Q.f (x)
Q.1(x) .

Proof. Note that P.f (x) can be written in this way P.f (x) = ΨQ.1 (δx .M )(P ′ .f ). Thus,
for any x, y ∈ E, we have that
|P.f (x) − P.f (y)| = |(ΨQ.1 (δx .M ) − ΨQ.1 (δy .M ))(P ′ .f )|

≤ k(ΨQ.1 (δx .M ) − ΨQ.1 (δy .M ))ktv . osc(P ′ .f ).

Lemma 6 implies that
k(ΨQ.1 (δx .M ) − ΨQ.1 (δy .M ))ktv ≤ g.kδx .M − δy .M ktv ≤ g.β(M ).
Using (1.1), we have osc(P ′ .f ) ≤ β(P ′ ). osc(f ).
This ends the proof of the lemma.



Lemma 10. For any integers p ≤ n, we have
gp,n − 1 ≤
bp,n ≤

n
X

(gk − 1)

n
Y

bk .gk,n .

k=p+1

k=p+1

k−1
Y

(bi gi ),

(2.1)

i=p+1

(2.2)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

15

Proof. Let us prove (2.1). By definition, we have Gp,n = Qp,n .1. Combining (1.7) and
(1.8) applied to unit function, we have
Qp−1,n (1) = Qp .[(Qp+1 , . . . , Qn ).1] = Gp × Mp .(Qp,n .1).
This implies that the functions Gp,n satisfy the following “backward” relations:
Gn,n = 1;

Gp−1,n = Gp × Mp .Gp,n .

Then, for any x, y ∈ E, we deduce that
Gp−1,n (x) Gp (x) (Mp .Gp,n )(x)
×
.
=
Gp−1,n (y) Gp (y) (Mp .Gp,n )(y)
{z
}
| {z } |
E1

E2

Notice that E1 ≤ gp (by definition), and by Lemma 8, we have E2 ≤ 1 + β(Mp ).(gp,n −
1). This shows the following backward inequalities:
gn,n = 1;

gp−1,n ≤ gp (1 + bp (gp,n − 1)).

(2.3)

We end the proof of (2.1) by induction.
To prove (2.2), we use the formulae
Pp−1,n .f =
Q

Qp−1,n .f
Gp × Mp .Qp,n .f
Mp .Qp,n .f
=
=
.
Qp−1,n .1
Gp × Mp .Qp,n .1
Mp .Qp,n .1
.f

, we apply Lemma 9 to check that β(Pp−1,n ) ≤ β(Mp ).gp,n β(Pp,n ),
Recalling that Pp,n .f = Qp,n
p,n .1
from which we conclude that
bp−1,n ≤ bp .gp,n .bp,n .
We end the proof of (2.2) by induction.
This ends the proof of the lemma.



We end this section with a useful technical lemma to control the quantity gp,n bp,n .
Lemma 11. For any p ≤ n, we have
gp,n bp,n ≤
Proof. Using Lemma 10, we have
gp,n bp,n ≤ gp,n .

n
Y

k=p+1

bk gk,n

n
Y

k=p+1

(bk .gk−1,n ).

16

F. Giraud and P. Del Moral
= gp,n .(bp+1 gp+1,n ).(bp+2 gp+2,n ) · · · (bn−1 gn−1,n ).(bn gn,n )
|{z}
=1

= (gp,n bp+1 ).(gp+1,n bp+2 ) · · · (gn−1,n bn ) =

n
Y

bk gk−1,n .

k=p+1

This ends the proof of the lemma.



The term gp,n bp,n is central in the Lp -mean error bound (1.10). By Lemma 11, we have
n
n
X
Y

bk gk−1,n < +∞

=⇒

p=0 k=p+1

n
X

gp,n bp,n < +∞.

p=0

This gives a sufficient condition for a uniform Lp bound w.r.t. time n. gp,n bp,n is also
2
involved in the estimates of all the quantities defined in Section 1.5 such as rn , βn , and
others. In addition, by Lemma 6, we have the stability property
kφp,n (µ) − φp,n (ν)ktv ≤ gp,n bp,n kµ − νktv .
This shows that the term gp,n bp,n is central to quantify the stability properties of the
semigroup φp,n .

2.2. Uniform concentration theorems
To obtain uniform bounds w.r.t., the time horizon in (1.10), Lemma 11 naturally leads
to a sufficient condition of the following type:
sup bk .gk−1,n ≤ a

k≤n

for some a ∈ (0, 1).

In this situation, we prove that gp,n bp,n ≤ an−p and, therefore, using (1.10),
p 1/p

sup E(|ηnN (f ) − ηn (f )| )

f ∈B1 (E)

Bp 1
,
≤√
N 1−a

with the constants Bp introduced in (0.1). We then fix the parameter a ∈ (0, 1) and we
look for conditions on the bp so that bk gk−1,n ≤ a. This parameter a can be interpreted
as a performance degree of the N -approximation model. In order to explicit relevant
and applicable conditions, we study two typical classes of regularity conditions on the
potential functions Gp . The first one relates to bounded coefficients gp (Theorem 12). In
the second one, the parameters gp tend to 1 as p → ∞ (Theorem 13).
Theorem 12. We assume that
sup gp ≤ M
p≥1

and

sup bp ≤
p≥1

a
M (1 + a)

(2.4)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

17

for some finite M ≥ 1 and some a ∈ (0, 1). In this situation, we have the following uniform
estimates:
• The Lp -error bound:
sup dp (ηnN , ηn ) ≤

n≥0

Bp
√ .
2(1 − a) N

(2.5)

• For any n ≥ 0, N ≥ 1, y ≥ 0 and f ∈ O1 (E), the probability of the event
√
r1⋆ (1 + h0 (y)) + r2⋆ N y
N
ηn (f ) − ηn (f ) ≤
N
is greater than 1 − e−y , with the parameters
r1⋆ =

4M 2 (1 + a)2
1−a

and

(2.6)

√
2 2
.
r2⋆ = √
1 − a2

(2.7)

• For any n ≥ 0, N ≥ 1, ǫ ∈ {+1, −1}, and y ≥ 0, the probability of the event



 N
ǫ
r̃1
y
γ (1)
≤ h0 (y) + r̃2 .h1
log n
n
γn (1)
N
n.N
is greater than 1 − e−y , with the parameters r̃1 =
functions h0 and h1 are defined in (1.9), page 12.

8M 3 (1+a)2
1−a

and r̃2 =

4M
1−a .

Proof. Firstly, we prove the pair of inequalities

gp,n ≤ M (1 + a),
bp .gp−1,n ≤ a.
If gp ≤ M and bp ≤ b :=
gp,n − 1 ≤
We substitute b =

n
X

k=p+1

1 − (bM )n−p
M −1
≤
.
1 − bM
1 − bM

to find

M −1
+ 1 = (M − 1)(1 + a) + 1 ≤ M (1 + a).
1 − a/(1 + a)

For the second inequality, we have
bp .gp−1,n ≤

b(M − 1)
+b
1 − a/(1 + a)

= (1 + a)

The

(2.9)

then by Lemma 10

(M − 1)(bM )k−p−1 = (M − 1)

a
M(1+a)

gp,n ≤

a
M(1+a) ,

(2.8)

a
a(M − 1)
+
M (1 + a) M (1 + a)

18

F. Giraud and P. Del Moral
=



a
1
aM − a +
M
1+a

≤ a.
This ends the proof of (2.9).

• By Lemma 11 and (2.9), we have gp,n bp,n ≤ an−p . Combining this with (1.10), the
Lp -error bound (2.5) is clear. See the relation (A.1), page 38 for the factor 1/2.
• Let us prove (2.6), which is a consequence of the concentration inequality (1.11).
2
Combining the estimations of rn and βn given in Section 1.5 with (2.9) and
gp,n bp,n ≤ an−p , we deduce that
rn ≤

4M 2 (1 + a)2
1−a

2

and βn ≤

4
.
1 − a2

Equation (2.6) is obtained by making the suitable substitutions in (1.11).
• The last concentration inequality (2.8) is a consequence of (1.12) and (2.9). Let us
recall that Lemma 11 and (2.9) imply gp,n bp,n ≤ an−p . Then, from estimations (1.13)
and (1.14), we can easily show that the quantities τn⋆ and r̄(n) satisfy
τn⋆ ≤

4M
n(1 − a)

and r̄(n) ≤

8M 3 (1 + a)2
.
1−a

On the other hand, σ̄n2 is trivially bounded by n. Then we find that
r


2y(τn⋆ σ̄n )2
r̄(n)
r̄(n)
y
yτn⋆
⋆ 2
=
h0 (y) + τn σ̄n h1
h0 (y) +
+
.
2
N
N.σ̄n
N
3N
N
Finally, (2.8) is obtained from (1.12) by making the suitable substitutions.
This ends the proof of the theorem.



Let us now consider the case where gp decreases to 1 as p → ∞. The idea of the
forthcoming analysis is to find a condition on the bp so that the gp,n are uniformly
1+α
bounded w.r.t. n by gp+1
with
α=

a
>0
1−a


⇐⇒ a =


α
.
1+α

The concentration inequalities developed in Theorem 13 will be described in terms of
the parameters r3⋆ (n) and r̃3 (n), r̃4 , r̃5 (n) defined below.
r3⋆ (n) =

4u1 (n)
,
1−a

(2.10)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models
19


 r̃3 (n) = 16.u2(n) ,


1−a




4X
4.g1
r̃4 =
,
gn+1 .an ≤
(2.11)
3
3(1 − a)

n≥0


√




 r̃5 (n) = 4 2.u3 (n) .
1−a
The sequences u1 (n), u2 (n) and u3 (n) used in the above formulae are defined by

n
X

2(1+α)


u1 (n) = (1 − a)
gn−p+1 ap −→ 1,


n→∞


p=0



n


1 X 3+2α
u2 (n) =
g
−→ 1,
n→∞
n p=1 p




!1/2


n−1


1X 2


u (n) =
gp+1
−→ 1.

 3
n→∞
n
p=0

Notice that the sequence u1 (n) tends to 1 by dominated convergence. Sequences u2 (n)
and u3 (n) tend to 1 by Cesaro’s theorem.
Theorem 13.
We assume that gp ↓ 1 as p → ∞ and the sequence bp satisfies for any p ≥ 1
bp ≤

gpα − 1

gpα+1 − 1

( −→ a)
p→+∞

and

bp ≤

a
( −→ a).
gpα+1 p→+∞

(2.12)

In this situation, we have the following uniform estimates:
• The Lp -error bound
sup dp (ηnN , ηn ) ≤
n≥0

Bp
√
2(1 − a) N

with the constants Bp introduced in (0.1).
• For any n ≥ 0, N ≥ 1, y ≥ 0 and f ∈ O1 (E), the probability of the event
√
r3⋆ (n)(1 + h0 (y)) + r2⋆ N y
N
ηn (f ) − ηn (f ) ≤
N

(2.13)

(2.14)

is greater than 1 − e−y , with the bounded sequence r3⋆ (n) defined in (2.10), the parameter r2⋆ defined in (2.7) and the function h0 defined in (1.9), page 12.
• For any n ≥ 0, N ≥ 1, ǫ ∈ {+1, −1}, and y ≥ 0, the probability of the event
r




 N
√ 
y+ y
y
γ (1)
ǫ
y
≤ r̃3 (n)
+ r̃4
+ r̃5 (n)
log n
(2.15)
n
γn (1)
N
n.N
n.N

20

F. Giraud and P. Del Moral
is greater than 1 − e−y , with the parameters r̃3 (n), r̃4 , r̃5 (n) defined in (2.11).

Proof. Firstly, we prove that
(2.12)

=⇒



∀p ≤ n,

gp,n ≤ (gp+1 )1+α ,
gp−1,n .bp ≤ a.

(2.16)

The proof of the first inequality comes from a simple backward induction on p (with
fixed n), using formula gp−1,n ≤ gp (1 + bp (gp,n − 1)) (see (2.3)). For p = n, gp,n is clearly
1+α
smaller than gp+1
because gn,n = 1. The second assertion is now immediate.
1+α
Now we assume that gp,n ≤ gp+1
. In this case, gp−1,n ≤ gp1+α is met as soon as
bp ≤

gpα − 1

α+1
gp+1
−1

.

Notice that this estimate is met as soon as bp ≤
decreasing.

gpα −1
α+1
gp −1

, the sequence (gp )p being

• Now that we proved (2.16) (which implies gp,n bp,n ≤ an−p by Lemma 11), the Lp mean error bound (2.13) comes from a simple substitution in (1.10). See the relation
(A.1), page 38 for the factor 1/2.
• To prove (2.14), we focus on the quantities βn and rn arising in the concentration
inequality (1.11). With (2.16) and gp,n bp,n ≤ an−p , we readily verify that
2

βn ≤

4
.
1 − a2
2(1+α)

4
g1
, but another manipulation
The term rn can be roughly bounded by 1−a
provides a more precise estimate. Indeed, using the fact that bp,n .gp,n ≤ an−p and
1+α
gp,n ≤ gp+1
, we prove that

rn ≤ 4
≤

n
X
p=0

2
gp,n
an−p ≤ 4

n
X
p=0

2(1+α) n−p

gp+1

a

≤4

n
X

2(1+α)

gn−p+1 ap

p=0

4.u1 (n)
.
1−a

We prove (2.14) by making the suitable substitutions in (1.11).
• Let us prove the last concentration inequality (2.15). It is mainly a consequence of
the inequality (1.12). Starting from the following decomposition,
r


2r̄(n)
y
√
r̄(n)
2y(τn⋆ σ̄n )2
yτn⋆
⋆ 2
=
h0 (y) + τn σ̄n h1
(y + y) +
+
,
2
N
N.σ̄n
N
3N
N

(2.17)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

21

we need to find some refined estimates of the quantities τn⋆ , r̄n and (τn⋆ .σ̄n )2 . To
estimate τn⋆ , we notice that ∀q, gp+1 ≤ gp−q+1 , so that
τq,n

n−q−1
n−1
n−1
4X
4X
4 X
p−q
p−q
≤
gp+1 ap
gp+1 a
≤
gp−q+1 a
≤
n p=q
n p=q
n p=0
n−1

≤

4X
gp+1 ap .
n p=0

P
p
0
Finally, we have that τn⋆ ≤ 4.U
p≥0 gp+1 a ≤
n , where U0 =
We estimate r̄n , using the following inequalities:
r̄n ≤
≤

8
n

g1
1−a .

n−1 p

X

0≤q≤p<n

3
gp+1 .gq,p
.bq,p ≤

n
1 X 3+2α
8
.
·
g
1 − a n p=1 p
{z
}
|

8 X X 3+2α p−q
g
a
n p=0 q=0 p+1

=u2 (n)

Let us conduct a last useful estimation:
2
(τn⋆ .σ̄n )

≤

n−1
X

≤

n−1
X

≤

n−1
X

≤

q=0

q=0

q=0

2
τq,n

≤

n−1
X
q=0

n−1

4X
gp+1 .ap−q
n p=q

n−q+1
4 X
gp+q+1 .ap
n p=0

!2

!2

16 2
1
·g
·
n2 q+1 (1 − a)2

n−1
1X 2
16
×
g .
n.(1 − a)2 n q=0 q+1
| {z }
=(u3 (n))2

At last, we make the suitable substitutions in (2.17) and obtain the desired inequality
(2.15).
This ends the proof of the theorem.



22

F. Giraud and P. Del Moral

3. Interacting simulated annealing models
3.1. Some motivations
We consider the Boltzmann–Gibbs probability measure associated with an inverse “temperature” parameter β ≥ 0 and a given potential function V ∈ B(E) defined by
µβ (dx) =

1 −β.V (x)
e
m(dx),
Zβ

(3.1)

where m stands for some reference measure, and Zβ is a normalizing constant. We let βn
a strictly increasing sequence (which may tend to infinity as n → ∞). In this case, the
measures ηn = µβn can be interpreted as a FK flow of measures with potential functions
Gn = e−(βn −βn−1 )V and Markov transitions Mn chosen as being MCMC dynamics for
the current target distributions. Indeed, we have

 −βn−1 .V (x)
Zβ
e−βn .V (x)
e
m(dx) = n−1 |e−(βn −β{zn−1 )V (x)}
m(dx) .
µβn (dx) =
Zβn
Zβn
Zβn−1
=Gn (x)
|
{z
}
=µβn−1 (dx)

This shows µβn = ψGn (µβn−1 ). Let φn stand for the FK transformation associated with
potential function Gn and Markov transition Mn . We have
φn (µβn−1 ) = ψGn (µβn−1 ) = µβn .Mn = µβn .

Sampling from these distributions is a challenging problem in many application domains. The simplest one is to sample from a complex posterior distribution on some
Euclidian space E = Rd , for some d ≥ 1. Let x be a variable of interest, associated with
a prior density p(x) (easy to sample) with respect to Lebesgue measure dx on E, and
y a vector of observations, associated with a calculable likelihood model p(y|x). In this
context, we recall that p(y|x) is the density of the observations given the variable of
interest. The density p(x|y) of the posterior distribution η is given by Bayes’ formula
p(x|y) ∝ p(x).p(y|x).
In the case where η is highly multimodal, it is difficult to sample from it directly. As an
example, classic MCMC methods tend to get stuck in local modes for very long times.
As a result, they converge to their equilibrium η only on unpractical time-scales. To
overcome this problem, a common solution is to approximate the target distribution η
with a sequence of measures η0 , . . . , ηn with density
ηk (dx) ∝ p(x).p(y|x)βk dx,
where (βk )0≤k≤n is a sequence of number increasing from 0 to 1, so that η0 is the prior
distribution of density p(x), easy to sample, and the terminal measure ηn is the target

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

23

distribution η (see, for instance, [4, 30, 37, 38]). If we take V := x 7→ − log(p(y|x)) and
m(dx) := p(x) dx, then the ηk coincide with the Boltzmann–Gibbs measures µβk defined
in (3.1). In this context, IPS methods arise as being a relevant approach, especially if η
is multimodal, since the use of a large number of particles allows to cover several modes
simultaneously.
The normalizing constant γn (1) coincide with the marginal likelihood p(y). Computing
this constant is another central problem in model selection arising in hidden Markov chain
problems and Bayesian statistics.
Next, we present another important application in physics and chemistry, known as
free energy estimation. The problem starts with an unnormalized density of the form


H(ω, α)
q(ω|T, α) = exp −
,
k.T
where H(ω, α) is the energy function of state ω, k is Boltzmann’s constant, T is the
temperature and α is a vector of system characteristics. The free energy F of the system
is defined by the quantity
F (T, α) = −k.T. log(z(T, α)),
where z(T, α) is the normalizing constant of the system density. See, for instance, [7, 12,
29] for a further discussion on these ground state energy estimation problems.
Last, but not least, it is well known that Boltzmann–Gibbs measures’ sampling is
related to the problem of minimizing the potential function V . The central idea is that
µβ tends to concentrate on V ’s minimizers as the inverse temperature β tends to infinity.
To be more precise, we provide an exponential concentration inequality in Lemma 14.
In this context, the IPS algorithm can be interpreted as a sequence of interacting simulated annealing (abbreviate ISA) algorithms. As they involve a population of N individuals, evolving according to genetic type processes (selection, mutation), ISA methods
also belong to the rather huge class of evolutionary algorithms for global optimization.
These algorithms consist in exploring a state space with a population associated with
an evolution strategy, that is, an evolution based on selection, mutation and crossover.
See [31] or [36] and the references therein for an overview. As these algorithms involve
complex, possibly adaptive strategies, their analysis is essentially heuristic, or sometimes
asymptotic (see [8] for general convergence results on genetic algorithms). The reader will
also find in [24] a proof for a ISA method of the a.s. convergence to the global minimum
in the case of a finite state space, when the time n tends to infinity, and as soon as the
population size N is larger than a critical constant that depends on the oscillations of
the potential fitness functions and the mixing properties of the mutation transitions.
The results of the previous sections apply to the analysis of ISA optimization methods. Our approach is nonasymptotic since we estimate at each time n, and for a fixed
population size N the distance between the theoretical Boltzmann–Gibbs measure ηn
and its empirical approximation ηnN . In some sense, as βn → +∞, ηn tends to the Dirac
measure δx⋆ where x⋆ = argminx∈E V (x), as soon as there is a single global minimum. So

24

F. Giraud and P. Del Moral

intuitively, we guess that if this distance admits a uniform bound w.r.t. time n, then for
large time horizon we have ηnN ≃ δx⋆ .
In this section, we propose to turn the conditions of Theorems 12 and 13 into conditions
on the temperature schedule to use, and the number of MCMC steps that ensure a
given performance degree. Then we combine the concentration results of Section 2.2
with Lemma 14 to analyze the convergence of the IPS optimization algorithm.

3.2. An ISA optimization model
We fix an inverse temperature schedule βn and we set
• ηn (dx) = µβn (dx) =

1
−βn V (x)
m(dx);
Zβn e

• Gn (x) = e−∆n .V (x) ;
• and then gn = e∆n . osc(V ) ,

where ∆n are the increments of temperature ∆n = βn − βn−1 . We let Kβ the simulated
annealing Markov transition with invariant measure µβ and a proposition kernel K(x, dy)
reversible w.r.t. m(dx). We recall that Kβ (x, dy) is given by the following formulae:
Kβ (x, dy) = K(x, dy). min(1, e−β(V (y)−V (x)) )
∀y 6= x,
Z
Kβ (x, {x}) = 1 −
K(x, dy). min(1, e−β(V (y)−V (x)) ).
y6=x

Under the assumption K k0 (x, ·) ≥ δν(·) for any x with some integer k0 ≥ 1, some measure
ν and some δ > 0, one can show (see, e.g., [3], Lemma 4.3.11) that
β(Kβk0 ) ≤ (1 − δe−β∆V (k0 ) ).

(3.2)

In the above display formula, ∆V (k0 ) = supx∈E kV (x) − V (X(x))k∞ , where for all
x ∈ E, X(x) is a random variable of law δx .K k0 . In other words, ∆V (k0 ) is the maximum
potential gap one can obtain making k0 elementary moves with the Markov transition
K. It is clearly bounded by osc(V ). One way to control the mixing properties of the
k .m
ISA model is to consider the Markov transition Mp = Kβp0 p , the simulated annealing
kernel iterated k0 .mp times. In this case, the user has a choice to make on two tuning
parameters, namely the temperature schedule βp , and the iteration numbers mp . Note
that for all b ∈ (0, 1), condition bp ≤ b in turned into (1 − δe−βp ∆V (k0 ) )mp ≤ b, and this
last condition is satisfied if we have
mp ≥

log(1/b)e∆V (k0 ).βp
.
δ

(3.3)

We prove now a technical lemma that we will use in the following. It deals with
Boltzmann–Gibbs measures’ concentration properties.

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

25

Lemma 14. For any β > 0, and for all 0 < ε′ < ε, the Boltzmann–Gibbs measure µβ
satisfies
′

µβ (V ≥ Vmin + ε) ≤

e−β(ε−ε )
,
mε′

where mε′ = m(V ≤ Vmin + ε′ ) > 0.
The normalizing constant Zβ of the definition (3.1) is necessary equal to
RProof.
−βV
dm. Then we have
Ee
R
e−βV dm
V ≥Vmin +ε
R
µβ (V ≥ Vmin + ε) = R
−βV dm
−βV dm +
V <Vmin +ε e
V ≥Vmin +ε e
−1
Z
 Z
−βV
−βV
e
dm
≤
e
dm
.
V <Vmin +ε
V ≥Vmin +ε
|
{z
}
{z
}
|
A1

A−1
2

Firstly, it is clear that A1 ≤ e−β(Vmin +ε) . Secondly, ε′ < ε implies {V ≤ Vmin + ε′ } ⊂ {V <
Vmin + ε}, then we have
Z
′
e−βV dm ≥ m(V ≤ Vmin + ε′ )e−β(Vmin +ε ) .
A2 ≥
V ≤Vmin +ε′

We end the proof by making the appropriate substitutions.



Combining Lemma 14, the theorems of Section 2.2 (with indicator test function f =
1{V ≥Vmin +ε} ), and the Dobrushin ergodic coefficient estimate (3.2) (and the associated
remark (3.3)) we prove the following theorem, which is the complete version of Theorem 2.
Theorem 15. Let us consider the N-particles ISA algorithm associated with the potential
functions Gn and Markov transitions Mn defined in the beginning of this section. Let us
fix a ∈ (0, 1). For any ε > 0, and n ≥ 0, let pN
n (ε) denote the proportion of particles
(ζni ) s.t. V (ζni ) ≥ Vmin + ε. We assume that the inverse temperature schedule βp and the
iteration numbers mp satisfy one of these two conditions:
1. supp≥1 ∆p ≤ ∆ < ∞ (e.g., linear temperature schedule) and
mp ≥

log(e∆. osc(V ) (1 + a)/a)e∆V (k0 ).βp
.
δ

2. ∆p ↓ 0 (as p → ∞) and mp ≥ (osc(V ).∆p + log( a1 )) e

∆V (k0 ).βp

δ
′

.

In this situation, for any ε > 0, n ≥ 0, N ≥ 1, y ≥ 0, and ε ∈ (0, ε), the probability of
the event
√
′
e−βn (ε−ε ) ri⋆ (n)(1 + h0 (y)) + r2⋆ N y
N
+
pn (ε) ≤
mε′
N

26

F. Giraud and P. Del Moral

is greater than 1 − e−y . In the above display formula, the constant r2⋆ is defined in (2.7),
page 17, the function h0 in (1.9), page 12 and ri⋆ (n) is defined by:
– ri⋆ (n) = r1⋆ (see (2.7), page 17 with M = e∆. osc(V ) ) in the case of bounded ∆p ;
– ri⋆ (n) = r3⋆ (n) (see (2.10), page 18 with gn = e∆n osc(V ) ) in the second one.
−βn (ε−ε′ )

We distinguish two error terms. The first one, ( e m ′ ), is related to the concenε
tration of the Boltzmann–Gibbs
measure around the set of global minima of V . The
⋆
⋆√
r (n)(1+h0 (y))+r2 N y
second one, ( i
), is related to the concentration of the occupation
N
measure around the limiting Boltzmann–Gibbs measure. Besides the fact that Theorem 15 provides tuning strategies which ensure the performance of the ISA model, the
last concentration inequality is explicit in the relative importance of other parameters,
including the probabilistic precision y, the threshold t on the proportion of particles
possibly out of the area of interest, the final temperature βn and the population size
N . As long as one can have rough estimates of mε′ and osc(V ) (spatial scale and variation of the function V ), a simple
equation, deduced from this last theorem, such as
√
−βn (ε−ε′ )

r ⋆ (n)(1+h (y))+r ⋆

Ny

0
2
= 2t ) provides some orders of magnitude that can be
= i
(e m ′
N
ε
useful to roughly design an ISA model at first sight, which is generally a difficult task.
One natural way to choose ∆p is to look at n1 , the number of transitions of type Kβk0 we
need to proceed to move from βp to βn , and to compare it to n2 , the number of transitions
we need to proceed to move from βp to βq , and then from βq to βn , with βp < βq < βn .
Roughly speaking, we have seen that the convergence condition was bp ≤ gap , then, using
(3.3), we have

• n1 ≃ (osc(V ).∆p,n + log( a1 )) e
• n2 ≃ (osc(V

∆V (k0 ).βn

δ
∆V (k0 ).βq
).∆p,q + log( a1 )) e δ

,
+ (osc(V ).∆q,n + log( a1 )) e

∆V (k0 ).βn

δ

,

where ∆p,q := βp − βq for p > q. After some approximation technique we find that
n1 ≤ n2

⇐⇒

∆p,q ∆q,n ≥

log(1/a)
.
osc(V ).∆V (k0 )

This condition does not bring any relevantP
information in the case where ∆p −→ 0
n
N
except that the error decomposition ηnN − ηn = p=0 φp,n (ηpN )− φp,n φp (ηp−1
), underlying
our analysis, is not adapted to the case where ∆p −→ 0 (which can be compared to the
continuous time case). Nevertheless, this condition is interesting in the case of constant
inverse temperature steps. In this situation, the critical parameter ∆β is given by
s
log(1/a)
.
∆β =
osc(V ).∆V (k0 )
More precisely, above ∆β the algorithm needs to run too many MCMC steps to stabilize
the system. In the reverse angle, when the variation of temperature is too small, it is
difficult to reach the desired target measure.

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

27

4. An adaptive temperature schedule in ISA
As we already mentioned in the Introduction, the theoretical tuning strategies developed
in Section 3.2 are of the same order as the logarithmic cooling schedule of traditional simulated annealing [3, 8, 24]. In contrast to SA models, we emphasize that the performance
of the ISA models are not based on a critical initial temperature parameter. Another
advantage of the ISA algorithm is to provide at any time step an N -approximation of
the target measure with a given temperature. In other words, the population distribution
reflects the probability mass distribution of the Boltzmann–Gibbs measure at that time.
Computationally speaking, the change of temperature parameter ∆p plays an important
role. For instance, if ∆p is taken too large, the selection process is dominated by a minority of well fitted particles and the vast majority of the particles are killed. The particle
set’s diversity, which is one of the main advantage of the ISA method, is then lost. On
the contrary, if ∆p is taken too small, the algorithm does not proceed to an appropriate
selection. It wastes time by sampling from MCMC dynamics while the set of particles has
already reached its equilibrium. The crucial point is to find a relevant balance between
maintaining diversity and avoiding useless MCMC operations.
Designing such a balance in advance is almost as hard as knowing the function V in
advance. Therefore, it is natural to implement adaptive strategies that depend on the
variability and the adaptation of the population particles (see, e.g., [13, 27, 34, 37, 39] for
related applications). In the general field of evolutionary algorithms, elaborating adaptive
selection strategies is a crucial question (see, e.g., [2]) and a challenging problem to design
performant algorithms. In the case of ISA methods, the common ways to choose ∆p are
based on simple criteria such as the expected number of particle killed (see Section 4.2),
or the variance of the weights (effective sample size). All of these criteria are based on
the same intuitive idea; that is to achieve a reasonable selection. As a result, all of these
adaptive ISA models tend to perform similarly.
In [18], the reader will find a general formalization of adaptive IPS algorithms. The idea
is to define the adaptation as the choice of the times n at which the resampling occurs.
These times are chosen according to some adaptive criteria, depending on the current
particle set, or more generally on the past process. Under weak conditions on the criteria,
it is shown how the adaptive process asymptotically converges to a static process involving
deterministic interaction times when the population size tends to infinity. A functional
central limit theorem is then obtained for a large class of adaptive IPS algorithms.
The approach developed in the following section is radically different, with a special
focus on nonasymptotic convergence results for the ISA algorithm defined in Section 4.2.
The adaptation consists here in choosing the β increment ∆N
n+1 so that
N

ηnN (e−∆n+1 ·V ) = ε,
where ε > 0 is a given constant, at each iteration n. We show that the associated stochastic process can be interpretated as a perturbation of the limiting FK flow.

28

F. Giraud and P. Del Moral

4.1. Feynman–Kac representation
Let V ∈ B(E). To simplify the analysis, without any loss of generality, we assume Vmin =
0. Let us fix ε > 0. For any measure µ ∈ P(E), we let the function λµ defined by
[0, +∞) −→ (0, 1],

λµ = x 7→ µ(e−x·V ).

This function is clearly decreasing (λµ (0) = 1), convex, and differentiable infinitely.
Moreover, if µ({V = 0}) = 0, then it satisfies λµ (x) −→ 0 when x → +∞. Therefore, we
can define its inverse function κµ :
(0, 1] −→ [0, +∞),
κµ = ε 7→ x

so that µ(e−x·V ) = ε.

This function is again decreasing, convex, infinitely differentiable, takes value 0 for ε = 1,
and it satisfies κµ (ε) −→ +∞ when ε → 0+ .
Now, we let m be a reference measure on E s.t. m({V = 0}) = 0. We consider the sequence (βn )n and its associated Gibbs measures ηn = µβn ∝ e−βn V .m, defined recursively
by the equation
∆n+1 := (βn+1 − βn ) = κηn (ε).

(4.1)

In an equivalent way, we have
ληn (∆n+1 ) = ε

or ηn (e−∆n+1·V ) = ε.

The main objective of this section is to approximate these target measures. Formally
speaking, (ηn ) admits the FK structure described in Section 3.1, with potential functions
Gn (x) = e−∆n .V (x) , and some dedicated MCMC Markov kernels Mn . We let gn , bn ,
be the associated oscillations, Dobrushin ergodic coefficients and the corresponding FK
transformations φn .
The solving of the equation ηn (e−∆n+1 ·V ) = ε can be interpreted as a way to impose
some kind of theoretical regularity in the FK flow. Indeed, according to the formula (1.6)
(and definition (3.1)), it is equivalent to find ∆n+1 s.t.
γn+1 (1)
=ε
γn (1)



Zβn +∆n+1
.
=
Zβn

In other words, the sequence ∆n is defined so that the normalizing constants γn (1)
increase geometrically, with the ratio γn+1 (1)/γn (1) = ε. Notice that these increments
∆n are only theoretical, and the corresponding potential functions Gn are not explicitly
known.

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

29

4.2. An adaptive interacting particle model
As in the classic IPS algorithm, we approximate the measures ηn by simulating an interacting particle system (ζn )n = (ζn1 , . . . , ζnN )n of size N so that
ηnN =

1 X
δζni →N ↑∞ ηn .
N
1≤i≤N

We start with N independent samples from η0 and then alternate selection and mutation steps, as described in Section 1.4. As we mentioned above, in contrast to the
classic IPS model, the potential function Gn+1 arising in the selection is not known. The
selection step then starts by calculating the empirical increment ∆N
n+1 defined by
N (ε)
∆N
n+1 := κηn

P
i
−∆·V (ζn
)
N −∆N
n+1 ·V ) = ε. As the quantity η N (e−∆·V ) = 1
or ληnN (∆N
n+1 ) = ηn (e
n
1≤i≤N e
N
is easy to calculate for all ∆ ≥ 0, one can calculate ∆N
n+1 by, for example, performing a
dichotomy algorithm. If we consider the stochastic potential functions
N

−∆n+1 .V
GN
n+1 = e

then every particle ζni evolves to a new particle ζbni randomly chosen with the following
stochastic selection transition:
N
i
N
i
N
i
Sn+1,η
(ηnN )(dx).
N (ζn , dx) := Gn+1 (ζn )δζ i (dx) + (1 − Gn+1 (ζn ))ΨGN
n
n+1
n

In the above display formula, ΨGN
(ηnN ) stands for the updated measure defined by
n+1
Ψ GN
(ηnN ) =
n+1

N
X
j=1

j
GN
n+1 (ζn )
δ j.
PN
N
k ζn
k=1 Gn+1 (ζn )

Note that Vmin = 0 ensures 0 < GN
n+1 ≤ 1.
For the mutation step, we consider two models. The implementable one consists in
N
performing Markov transitions Mn+1
(ζbni , ·), defined as Mn+1 (ζbni , ·) by replacing βn+1 by
N
N
N
βn+1 = βn + ∆n+1 . These kernels must be defined to let the associated Boltzmann–Gibbs
measures stable, for instance, using some simulated annealing kernels (see Section 3.2,
page 24). Thus, conditionally to the previous particle set ζn , the new population of
particles ζn+1 is sampled from distribution
1
N
Law(ζn+1
, . . . , ζn+1
|ζn1 , . . . , ζnN )
N
N
N
N
= (δζn1 .Sn+1,η
N .Mn+1 ) ⊗ · · · ⊗ (δζ N .Sn+1,η N .Mn+1 ).
n
n
n

(4.2)

30

F. Giraud and P. Del Moral

We also consider a simplified model where the mutation transition is given by the
limiting transition Mn+1 (MCMC Markov kernel associated with theoretical temperature
βn+1 defined in Section 4.1):
1
N
Law(ζn+1
, . . . , ζn+1
|ζn1 , . . . , ζnN )
N
N
= (δζn1 .Sn+1,η
N .Mn+1 ) ⊗ · · · ⊗ (δζ N .Sn+1,η N .Mn+1 ).
n
n
n

(4.3)

The definition of ∆N
n+1 is to be interpreted as the natural approximation of the theoretical relation (4.1). On the other hand, it admits a purely algorithmic interpretation.
As a matter of fact, conditionally to the nth generation of particles (ζn1 , . . . , ζnN ), the
probability for any particle ζni to be accepted, that is, not affected by the recycling
i
−∆N
i
n+1 .V (ζn ) . Then the expectation of the number
mechanism, is given by GN
n+1 (ζn ) = e
P
i
N
of accepted particles is given by i e−∆n+1 .V (ζn ) . But it turns out that this quantity
N
is exactly N × ηnN (e−∆n+1 ·V ), which is equal to N.ε by definition of ∆N
n+1 . Therefore,
ε is an approximation of the proportion of particles which remain in place during the
selection step. In other words, at each generation n, the increment ∆N
n+1 is chosen so
that the selection step kills less than (1 − ε).N particles. This type of tuning parameter
is very important in practice to avoid degenerate behaviors.

4.3. A perturbation analysis
This section is mainly concerned with the convergence analysis of the simplified adaptive
model (4.3). The analysis of the adaptive model (4.2) is much more involved, and our
approach does not apply directly to study the convergence of this model.
Despite the adaptation, the sequence ηnN can be analyzed as a random perturbation
of the theoretical sequence ηn . Let us fix n and a population state ζn at time n. If φN
n+1
denotes the FK transformation associated with potential GN
n+1 and kernel Mn+1 , then,
N
N
by construction, the measure ηn+1
is close to φN
n+1 (ηn ). In particular, by the Khintchine’s
type inequalities presented in [16] (see Lemma 7.3.3 page 223), we have
∀f ∈ B1 (E)

p

1/p

N
N
E(|ηn+1
(f ) − φN
n+1 (ηn )(f )| |ζn )

Bp
≤√ ,
N

(4.4)

with the constants Bp introduced in (0.1). A simple, but important remark about the
Boltzmann–Gibbs transformations is that for any measure µ and any positive functions
G and G̃ we have
ψG̃ (µ) = ψG (ψG̃/G (µ)).
GN

N
Therefore, if we take Hn+1
:= Gn+1
, then the perturbed transformation φN
n+1 can be
n+1
written in terms of the theoretical one φn+1 by

φN
N
.
n+1 = φn+1 ◦ ψHn+1

(4.5)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

31

If we use an inductive approach, we face the following problem. Let η be a deterministic
measure (ηn in our analysis) and η̂ a random measure (ηnN in our analysis), close to η
under the dp distance (induction hypothesis). We also consider a Markov kernel M and
the potential functions
G = e−κη (ε).V ,

Ĝ = e−κη̂ (ε).V ,

Ĥ =

Ĝ
G

(4.6)

and we let φ (resp., φ̂) be the FK transformation associated with the potential function G
(resp., Ĝ). The question is now: how can we estimate dp (φ(η), φ̂(η̂)) in terms of dp (η, η̂)?
To answer to this question, we propose to achieve a two-step estimation. Firstly, we
estimate the distance between η̂ and ψĤ (η̂) (Lemma 16). Secondly, we analyze the stability properties of the transformation φ (Lemma 17). This strategy is summarized by
the following synthetic diagram:
η

/ φ(η)
φ

η̂
↓
ψĤ (η̂)

/ φ̂(η̂).
φ

Lemma 16. Let η ∈ P(E), η̂ ∈ PΩ (E) and let G, Ĝ, Ĥ be the positive functions on E
defined by the equations (4.6) for some ε > 0. If η({V = 0}) = η̂({V = 0}) = 0 (a.s.), then
we have
dp (ψĤ (η̂), η̂) ≤

Vmax · eκη (ε).Vmax
· dp (η̂, η).
ε · η(V )

Proof. We simplify the notation and we set x := κη (ε) and x̂ := κη̂ (ε).
We start with the following observation:
ψĤ (η̂)(f ) − η̂(f ) =

η̂(Ĥ.f )
η̂(Ĥ)

− η̂(f ) =

1
η̂[(Ĥ − η̂(Ĥ)).f ]
{z
}
η̂(Ĥ) |
| {z }
A2

(4.7)

A1

for any f ∈ B(E). We notice that Ĥ = Ĝ/G = e(x−x̂)·V , which leads to the lower bound
η̂(Ĥ) = η̂(e(x−x̂)·V ) ≥ η̂(e−x̂·V ) = ε. The last equality comes from the definition of x̂. We
just proved: |A1 | ≤ ε−1 . On the other hand, we have osc(Ĥ) = |e(x−x̂)·Vmax − 1|, so that
|A2 | ≤ η̂(|Ĥ − η̂(Ĥ)|) · kf k∞ ≤ osc(Ĥ) · kf k∞
≤ |e(x−x̂)·Vmax − 1| · kf k∞ .

(4.8)

32

F. Giraud and P. Del Moral

The quantity û := (e(x−x̂)·Vmax − 1) is intuitively small. Next, we provide an estimate
in terms of the functions λη and λη̂ . Given ω ∈ Ω, if x ≥ x̂, then we can write
Z x
−λ′η̂ (s) ds.
λη (x) −λη̂ (x) = λη̂ (x̂) − λη̂ (x) =
| {z }
x̂
=ε=λη̂ (x̂)

Furthermore, for all µ ∈ P(E) and s ≥ 0, we have

−λ′µ (s) = µ(V · e−sV ) ≥ µ(V · e−sVmax )
≥ µ(V ) · e−sVmax .
Then we have
λη (x) − λη̂ (x) ≥ η̂(V )

−1 −sVmax x
]x̂
[e
Vmax

= η̂(V )

e−xVmax (x−x̂)Vmax
− 1)
(e
Vmax

= η̂(V )

e−xVmax
û.
Vmax

By symmetry, we have
x ≤ x̂

=⇒

λη̂ (x) − λη (x) ≥ η̂(V )

e−xVmax
(−û).
Vmax

This yields the almost sure upper bound
|û| · η̂(V )

e−xVmax
≤ |λη (x) − λη̂ (x)|.
Vmax

Using the decomposition η̂(V ) = η(V ) + (η̂(V ) − η(V )), by simple manipulation, we
prove that
|û| ≤

Vmax exVmax
|η̂(V ) − η(V )|
|λη (x) − λη̂ (x)| +
· |û| .
|{z}
η(V )
η(V )
|
{z
}
|
{z
}
A5
A3
A4

p

Considering the L norm of the right-hand side of this inequality, one can check that
• A3 = (η − η̂)(e−x·V ) so, as osc(e−x·V ) ≤ 1, kA3 kp ≤ dp (η̂, η);
max
• as osc(V ) = Vmax , kA4 kp ≤ Vη(V
) · dp (η̂, η);
xVmax −x̂Vmax
−xVmax
• A5 = |û| = e
|e
−e
| ≤ exVmax .
Making the appropriate substitutions, we have
kûkp ≤

2Vmax ex·Vmax
· dp (η̂, η).
η(V )

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

33

Combining this result with (4.7) and (4.8), we check that
kψĤ (η̂)(f ) − η̂(f )kp ≤

2Vmax exVmax
· dp (η̂, η) · kf k∞.
ε · η(V )

)
by noticing that ψĤ (η̂)(f ) − η̂(f ) is equal to 0 for any
We finally go from kf k∞ to osc(f
2
min
constant function f , and by considering the above inequality taken for f˜ = f − fmax +f
,
2
osc(f )
˜
which satisfies kf k∞ = 2 .
This ends the proof of the lemma.


Lemma 17. Let η ∈ P(E), η̂ ∈ PΩ (E), and let φ be a FK transformation associated with
a positive function G and a Markov kernel M . If we set g := supx,y∈E G(x)/G(y) and
b := β(M ), then we have
dp (φ(η̂), φ(η)) ≤ g · b · dp (η̂, η).
Proof. Let us fix f ∈ B(E). We have
φ(η̂)(f ) − φ(η)(f ) =
=

η̂(G × M.f )
− φ(η)(f )
η̂(G)
η̂(G × [M.(f − φ(η)(f ))])
.
η̂(G)

Let f˜ = M.(f − φ(η)(f )). By property (1.1), f˜ satisfies osc(f˜) = osc(M.f ) ≤ b · osc(f ).
)
) = 0. So we obtain
Additionally we have η(G × f˜) = η(G × M.f ) − η(G × η(G×M.f
η(G)
φ(η̂)(f ) − φ(η)(f ) =

η̂(G × f˜) η(G × f˜)
1
=
−
(η̂ − η)(G × f˜)
η̂(G)
η̂(G)
η̂(G)
| {z }
=0



Gmax
G
˜
×f .
=
(η̂ − η)
η̂(G)
Gmax
max
Firstly, we notice that | G
η̂(G) | ≤
rewritten as

Gmax
Gmin

≤ g. On the other hand, we notice that f˜ can be

f˜ = M.(f − ψG (η)(M.f )) = (M.f ) − ψG (η)(M.f ).

G
× f˜) ≤ osc(f˜) ≤
It follows that f˜max ≥ 0 and f˜min ≤ 0. Under these conditions, osc( Gmax
b · osc(f ). We conclude that


 p 1/p
 p 1/p



Gmax
G
G
˜
˜
E
×f
×f
(η̂ − η)
≤ g · E (η̂ − η)
η̂(G)
Gmax
Gmax
≤ g · b · dp (η̂, η) · osc(f ).

34

F. Giraud and P. Del Moral

This ends the proof of the lemma.



Remark 18. Lemma 17 holds in the case where η is also a random measure (we now
note η̃ ∈ PΩ (E)) if one can find a σ-algebra F such that:
1. η̃ is F -measurable;
2. the (F -measurable) random variable
dF
p (η̂, η̃) :=

sup
f˜∈O1Ω,F (E)

p
1/p
E[|η̂(f˜) − η̃(f˜)| |F ]

is uniformly bounded on Ω. In the above definition, O1Ω,F (E) denotes the set of random,
F -measurable functions f˜: E → R satisfying osc(f˜) ≤ 1 a.s.
More precisely, under these conditions, we have
dp (φ(η̂), φ(η̃)) ≤ g · b · kdF
p (η̂, η̃)k∞ .
Proof. We fix f ∈ B(E) and use the same line of arguments as in the proof of Lemma 17.
Since the function f˜ = M.(f − φ(η)(f )) is F -measurable and satisfies


G
× f˜ ≤ osc(f˜) ≤ b · osc(f )
a.s.,
osc
Gmax
we have

 p 1/p

Gmax
G
˜
E
×f
(η̂ − η̃)
η̂(G)
Gmax
 

 p 1/p
G
≤ g · E E (η̂ − η̃)
× f˜
F
Gmax

 

G/Gmax × f˜
≤ g · b · osc(f ) · E E (η̂ − η̃)
b · osc(f )

p

F

1/p

≤ g · b · osc(f ) · kdF
p (η̂, η̃)k∞ .

This ends the proof.



4.4. Nonasymptotic convergence results
This section is mainly concerned with the proof of Theorem 3 stated on page 6. We also
deduce some concentration inequalities of the ISA adaptive model. In this section, (ηn )
denotes the sequence of theoretical measures defined in Section 4.1, and (ηnN ) denotes the
sequence of empirical measures associated with the particle system described in (4.3),
page 30. We start with the proof of Theorem 3.

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

35

P
Q
Proof of Theorem 3. We fix p ≥ 1 and we let ẽn = nk=0 ni=k+1 bi gi (1 + ci ). We
notice that this sequence can also be defined with the recurrence relation ẽn+1 = 1 +
gn+1 bn+1 (1 + cn+1 ) · ẽn starting at ẽ0 = 1. We also consider the following parameter:
√
2 N
· sup kη N (f ) − ηn (f )kp .
en :=
Bp f ∈O1 (E) n

We use an inductive proof to check that the proposition IH(n) = {en ≤ ẽn } is met at any
rank n. As η0N is obtained with N independent samples from η0 , IH(0) is given by the
Khintchine’s inequality (see the relation (A.1), page 38 for the factor 1/2). Now suppose
that IH(n) is satisfied.
According to the identity (4.5), we can write the following decomposition:
N
N
N
N
N
ηn+1
− ηn+1 = (ηn+1
− φN
n+1 (ηn )) + (φn+1 (ηn ) − ηn+1 )
N
N
(ηnN )) − φn+1 (ηn )) .
= (ηn+1
− φN
(η N )) + (φn+1 (ψHn+1
{zn+1 n } |
|
{z
}
A1

(4.9)

A2

Using (4.4), page 30, we know that for all function f ∈ O1 (E), we have
p

1/p

N
N
kA1 (f )kp = E(E[|ηn+1
(f ) − φN
n+1 (ηn )(f )| |ζn ])

Bp
≤ √
2 N
√

(see the relation (A.1), page 38 for the factor 1/2), so 2BpN kA1 (f )kp ≤ 1.
To estimate A2 , we start by decomposing (ψHn+1
N
(ηnN ) − ηn ) in this way:
ψHn+1
N
(ηnN ) − ηn = (ψHn+1
N
(ηnN ) − ηnN ) + (ηnN − ηn ) .
{z
} | {z }
|
Q1

By the induction hypothesis, we have
Lemma 16, we find that

√
2 N
Bp

Q2

· supf ∈O1 (E) kQ2 (f )kp ≤ ẽn . Therefore, by

√
2 N
· sup kQ1 (f )kp ≤ cn+1 · ẽn .
Bp f ∈O1 (E)
Thus, the measures ψHn+1
N
(ηnN ) and ηn satisfy
√
2 N
· sup kψ N (η N )(f ) − ηn (f )kp ≤ (1 + cn+1 ) · ẽn .
Bp f ∈O1 (E) Hn+1 n

36

F. Giraud and P. Del Moral

Applying Lemma 17, we also have
√
2 N
· sup kA2 (f )kp ≤ gn+1 bn+1 (1 + cn+1 ) · ẽn .
Bp f ∈O1 (E)
Back to (4.9), we conclude that en+1 ≤ 1 + gn+1 bn+1 (1 + cn+1 ) · ẽn = ẽn+1 .
This ends the proof of the theorem.



We are now in position to obtain a sufficient condition for uniform concentration
and Lp -mean error bounds w.r.t. time for the simplified adaptive model discussed in
Sections 4.2 and 4.3.
Corollary 19. If the condition bn gn (1 + cn ) ≤ a is satisfied for some a < 1 and any n,
then we have the uniform error bounds
dp (ηnN , ηn ) ≤

Bp
√
2(1 − a) N

(4.10)

for any p, with the constants Bp introduced in (0.1). In addition, for any f ∈ B1 (E), we
have the following concentration inequalities:
√
2
(4.11)
∀s ≥ 0
P(|ηnN (f ) − ηn (f )| ≥ s) ≤ r1 ( N .s)e−r2 N s
and
∀y ≥ 0



√
r(1 + 2y)
√
P |ηnN (f ) − ηn (f )| ≥
≤ e−y ,
N

(4.12)

with the parameters

r1 = e1/2 (1 − a),




1
r2 = (1 − a)2 ,
2



r = 1 .
1−a

Proof. The inequality (4.10) is a direct consequence of Theorem 3 (see the relation
(A.1), page 38 for the factor 1/2).
Let us fix n, f ∈ B1 (E) and set
X := |ηnN (f ) − ηn (f )| and ǫN :=

1
√ .
(1 − a) N

In this notation, we have kXkp ≤ Bp · ǫN for any p ≥ 1. Let us fix s ≥ 0. By Markov
inequality, for all t ≥ 0 we have
P(X ≥ s) = P(etX ≥ ets ) ≤ e−st E(etX ).

(4.13)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

37

Using the formulae (0.1), page 6, we estimate the Laplace transform E(etX ).
E(e

tX

X  tp .X p 
)=
E
p!
p≥0

≤

X t2p ǫ2p (2p)!
N

p≥0

(2p)!

2p .p!

+

2 2
ǫN /2

≤ (1 + tǫN )et

Taking the inequality (4.13) with t =
P(X ≥ s) ≤

X t2p+1 ǫ2p+1
N

p≥0

(2p + 1)!

(2p + 1)!
√
2p + 1

2p .p!

.

1
ǫN

( ǫsN − 1), we obtain

s −(1/2)[(s/ǫN )2 −1]
e
,
ǫN

which is equivalent to the first concentration inequality (4.11).
For the second one, let u = ǫsN . Since log(u) ≤ u − 1 we have
2

P(X ≥ s) = P(X ≥ ǫN · u) ≤ u · e−(u

/2−1/2)

2

= e−(u

/2−log u−1/2)

2

≤ e−(u

/2−u+1/2)

.

Let T the function defined by
[1, +∞) −→ [0, +∞),
T = v 7→

1
v2
−v+
2
2

T is bijective, and its inverse function is T −1 : y 7→ 1 +
∀y ≥ 0

√
2y. Thus, we have

P(X ≥ ǫN · T −1 (y)) ≤ e−y ,

which is equivalent to (4.12).
This ends the proof of the corollary.

Appendix
A.1. Proof of Lemma 6, page 8
For all η ∈ P(E), let Sη be the Markov kernel defined by
Sη (x, dy) =



G(x)
G(x)
ΨG (η)(dy).
δx (dy) + 1 −
Gmax
Gmax



38

F. Giraud and P. Del Moral

This kernel satisfies η.Sη = ΨG (η), since for all f ∈ B(E) we have


η(G)
η(G × f )
. ΨG (η)(f )
+ 1−
η.Sη (f ) =
| {z }
Gmax
Gmax

=η(G×f )/η(G)

=

η(G × f ) η(G × f ) η(G).η(G × f )
+
−
Gmax
η(G)
Gmax .η(G)

= ΨG (η)(f ).
We use the decomposition
ΨG (µ) − ΨG (ν) = µ.Sµ − ν.Sν = (µ − ν).Sµ + ν.(Sµ − Sν ).
On the one hand, k(µ − ν).Sµ ktv ≤ kµ − νktv .β(Sµ ) ≤ kµ − νktv .
On the other hand, for all nonnegative function f ∈ B1 (E),


G(x)
(Sµ − Sν )(f )(x) = 1 −
[ΨG (µ)(f ) − ΨG (ν)(f )]
{z
}
Gmax |
|
{z
} |−|≤kΨG (µ)−ΨG (ν)ktv
≤1−1/g

so |ν.(Sµ − Sν )(f )| ≤ (1 − g1 )kΨG (µ) − ΨG (ν)ktv .
Thus,


1
kΨG (µ) − ΨG (ν)ktv
kΨG (µ) − ΨG (ν)ktv ≤ kµ − νktv + 1 −
g
⇒

kΨG (µ) − ΨG (ν)ktv ≤ g.kµ − νktv .

This ends the proof of the lemma.

A.2. Proof of the Lp -mean error bound (1.10) (page 12)
We fix p ≥ 1. Let us start with this remark: for all µ̂, ν̂ ∈ PΩ (E), µ̂(f ) − ν̂(f ) = 0 for all
constant function f , so we have
sup kµ̂(f ) − ν̂(f )kp = 2

kf k∞ ≤1

sup

kµ̂(f ) − ν̂(f )kp

kf k∞ ≤1/2

= 2 sup kµ̂(f ) − ν̂(f )kp

(A.1)

osc(f )=1

= 2 · dp (µ̂, ν̂).
In the context of the interacting particle system described in Section 1.4, page 10, we
use the decomposition
ηnN − ηn =

n
X

k=0

N
φk,n (ηkN ) − φk,n (φk (ηk−1
))

(A.2)

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

39

N
(with the convention φ0 (η−1
) = η0 ). For all k ∈ {1, 2, . . ., n}, let Fk−1 be the σ-algebra
generated by the variable ζk−1 (we set F0 = {∅, E}). Remark 18, page 34, implies
N
N
k−1
dp (φk,n (ηkN ), φk,n (φk (ηk−1
))) ≤ gk,n · bk,n · kdF
(ηkN , φk (ηk−1
))k∞ ,
p

(A.3)

F

N
where dp k−1 (ηkN , φk (ηk−1
)) is the random variable defined by
N
dpFk−1 (ηkN , φk (ηk−1
)) :=

sup
Ω,F
f˜∈O1 k−1 (E)

p
1/p
N
E[|ηkN (f˜) − φk (ηk−1
)(f˜)| |Fk−1 ] .

Ω,F

In the above definition, O1 k−1 (E) denotes the set of random, Fk−1 -measurable functions f˜: E → R satisfying osc(f˜) ≤ 1 a.s.
Given the conditioning relation
N
1
, . . . , ζk−1
)
Law(ζk1 , . . . , ζkN |ζk−1
1
N
N
N
.Mk )
.Sk,ηk−1
= (δζk−1
.Sk,ηk−1
.Mk ) ⊗ · · · ⊗ (δζk−1

and the Khintchine’s type inequalities presented in [16], we have for all random, Fk−1 measurable function f˜ satisfying kf˜k∞ ≤ 1 a.s.:
Bp
p
1/p
N
E(|ηkN (f˜) − φk (ηk−1
)(f˜)| |Fk−1 ) ≤ √ ,
N
with the constants Bp introduced in (0.1), page 6. This is equivalent to
Bp
N
k−1
dF
(ηkN , φk (ηk−1
)) ≤ √ .
p
2 N
Finally, we combine this result with (A.3) and the decomposition (A.2) to obtain
n
Bp X
dp (ηnN , ηn ) ≤ √
gk,n bk,n .
2 N k=0

This ends the proof.

References
[1] Assaraf, R., Caffarel, M. and Khelif, A. (2000). Diffusion Monte Carlo with a fixed
number of walkers. Phys. Rev. E (3).
[2] Baker, J.E. (1985). Adaptive selection methods for genetic algorithms. ICGA1 101–111.
[3] Bartoli, N. and Del Moral, P. (2001). Simulation & Algorithmes Stochastiques.
Toulouse, France: Cépaduès.
[4] Bertrand, C., Hamada, Y. and Kado, H. (2001). MRI prior computation and parallel
tempering algorithm: A probabilistic resolution of the MEG/EEG inverse problem.
Brain Topography 14 57–68.

40

F. Giraud and P. Del Moral

[5] Caffarel, M. and Assaraf, R. (2000). A pedagogical introduction to quantum MonteCarlo. In Mathematical Models and Methods for Ab Initio Quantum Chemistry. Lecture
Notes in Chemistry 74 45–73. Berlin: Springer. MR1855574
[6] Cappé, O., Moulines, E. and Rydén, T. (2005). Inference in Hidden Markov Models.
Springer Series in Statistics. New York: Springer. MR2159833
[7] Ceperley, D.M. (1995). Path integrals in the theory of condensed helium. Rev. Modern
Phys. 67 279–355.
[8] Cerf, R. (1994). Une théorie asymptotique des algorithmes génétiques. Ph.D. thesis, Université de Montpellier 2, France.
[9] Cérou, F., Del Moral, P., Furon, T. and Guyader, A. (2012). Sequential Monte Carlo
for rare event estimation. Stat. Comput. 22 795–808. MR2909622
[10] Cérou, F., Del Moral, P. and Guyader, A. (2011). A nonasymptotic theorem for
unnormalized Feynman–Kac particle models. Ann. Inst. Henri Poincaré Probab. Stat.
47 629–649. MR2841068
[11] Chopin, N. (2004). Central limit theorem for sequential Monte Carlo methods and its
application to Bayesian inference. Ann. Statist. 32 2385–2411. MR2153989
[12] Ciccotti, G. and Hoover, W.G., eds. (1986). Molecular-Dynamics Simulation of
Statistical-Mechanical Systems. Proceedings of the International School of Physics “Enrico Fermi”, XCVII. Amsterdam: North-Holland. MR0960486
[13] Clapp, T. (2000). Statistical Methods in the Processing of Communications Data. Ph.D.
thesis, Cambridge Univ., Engineering Dept.
[14] Dawson, D.A. and Del Moral, P. (2005). Large deviations for interacting processes in
the strong topology. In Statistical Modeling and Analysis for Complex Data Problems
(P. Duchesne and B. Rémillard, eds.). GERAD 25th Anniv. Ser. 1 179–208. New
York: Springer. MR2189537
[15] Del Moral, P. (1996). Nonlinear filtering: Interacting particle solution. Markov Process.
Related Fields 2 555–579. MR1431187
[16] Del Moral, P. (2004). Feynman–Kac Formulae: Genealogical and Interacting Particle
Systems with Applications. Probability and Its Applications (New York). New York:
Springer. MR2044973
[17] Del Moral, P., Doucet, A. and Jasra, A. (2006). Sequential Monte Carlo samplers. J.
R. Stat. Soc. Ser. B. Stat. Methodol. 68 411–436. MR2278333
[18] Del Moral, P., Doucet, A. and Jasra, A. (2012). On adaptive resampling strategies
for sequential Monte Carlo methods. Bernoulli 18 252–278. MR2888706
[19] Del Moral, P. and Guionnet, A. (1998). Large deviations for interacting particle
systems: Applications to non-linear filtering. Stochastic Process. Appl. 78 69–95.
MR1653296
[20] Del Moral, P. and Guionnet, A. (1999). Central limit theorem for nonlinear filtering
and interacting particle systems. Ann. Appl. Probab. 9 275–297. MR1687359
[21] Del Moral, P. and Guionnet, A. (2001). On the stability of interacting processes with
applications to filtering and genetic algorithms. Ann. Inst. Henri Poincaré Probab.
Stat. 37 155–194. MR1819122
[22] Del Moral, P., Hu, P. and Wu, L. (2012). On the concentration properties of interacting
particle processes. Foundations and Trends in Machine Learning 3 225–389.
[23] Del Moral, P. and Ledoux, M. (2000). Convergence of empirical processes for interacting
particle systems with applications to nonlinear filtering. J. Theoret. Probab. 13 225–
257. MR1744985

Nonasymptotic analysis of adaptive and annealed Feynman–Kac particle models

41

[24] Del Moral, P. and Miclo, L. (1999). On the convergence and applications of generalized
simulated annealing. SIAM J. Control Optim. 37 1222–1250. MR1691939
[25] Del Moral, P. and Miclo, L. (2000). Branching and interacting particle systems
approximations of Feynman–Kac formulae with applications to non-linear filtering.
In Séminaire de Probabilités, XXXIV. Lecture Notes in Math. 1729 1–145. Berlin:
Springer. MR1768060
[26] Del Moral, P. and Rio, E. (2011). Concentration inequalities for mean field particle
models. Ann. Appl. Probab. 21 1017–1052. MR2830611
[27] Deutscher, J., Blake, A. and Reid, I. (2000). Articulated body motion capture by
annealed particle filtering. IEEE Conference on Computer Vision and Pattern Recognition 2 126–133.
[28] Doucet, A., de Freitas, N. and Gordon, N., eds. (2001). Sequential Monte Carlo
Methods in Practice. Statistics for Engineering and Information Science. New York:
Springer. MR1847783
[29] Frankel, D. and Smit, B. (1996). Understanding Molecular Simulation. New York: Academic Press.
[30] Giraud, F., Minvielle, P., Sancandi, M. and Del Moral, P. (2012) Rao-blackwellised
interacting Markov chain Monte Carlo for electromagnetic scattering inversion. Available at arXiv:stat.AP/1209.4006v2.
[31] Goldberg, D.E. (1989). Genetic Algorithms in Search, Optimization, and Machine Learning. Reading, MA: Addison-Wesley.
[32] Hajek, B. (1988). Cooling schedules for optimal annealing. Math. Oper. Res. 13 311–329.
MR0942621
[33] Hetherington, J.H. (1984). Observations on the statistical iteration of matrices. Phys.
Rev. A (3) 30 2713–2719. MR0763130
[34] Jasra, A., Stephens, D.A., Doucet, A. and Tsagaris, T. (2011). Inference for Lévydriven stochastic volatility models via adaptive sequential Monte Carlo. Scand. J. Stat.
38 1–22. MR2760137
[35] Künsch, H.R. (2005). Recursive Monte Carlo filters: Algorithms and theoretical analysis.
Ann. Statist. 33 1983–2021. MR2211077
[36] Le Riche, R., Schoenauer, M. and Sebag, M. (2007). Un état des lieux de
l’optimisation évolutionnaire et de ses implications en sciences pour l’ingénieur, chapter in modélisation numérique:. In Défis et Perspectives, Vol. 2, Traité Mécanique et
Ingénierie des Matériaux (P. Breitkopf and C. Knopf-Lenoir, eds.) 187–259. Hermes.
[37] Minvielle, P., Doucet, A., Marrs, A. and Maskell, S. (2010). A Bayesian approach
to joint tracking and identification of geometric shapes in video sequences. Image and
Vision Computing. 28 111–123.
[38] Neal, R.M. (2001). Annealed importance sampling. Stat. Comput. 11 125–139.
MR1837132
[39] Schäfer, C. and Chopin, N. (2013). Sequential Monte Carlo on large binary sampling
spaces. Stat. Comput. 23 163. DOI:10.1007/s11222-011-9299-z.
[40] Schweizer, N. (2012). Non-asymptotic error bounds for sequential MCMC and stability
of Feynman–Kac propagators. Working paper, University of Bonn.
[41] Whiteley, N. (2012). Sequential Monte Carlo samplers: Error bounds and insensitivity to
initial conditions. Stoch. Anal. Appl. 30 774–798. MR2966098
Received September 2012 and revised September 2014

