Convergence rate of a simulated annealing algorithm with noisy
observations
Clément Bouttier

∗1,2,3

and Ioana Gavra

†3

arXiv:1703.00329v1 [stat.ML] 1 Mar 2017

1

Airbus Operations SAS, 316 route de Bayonne 31060 Toulouse Cedex 9, France
2
ENAC, 7 avenue Édouard Belin 31055 Toulouse Cedex 4, France
3
Institut de Mathématiques de Toulouse, Université Toulouse III 118 route de Narbonne 31062
Toulouse Cedex 9, France
December, 2016

Abstract
In this paper we propose a modified version of the simulated annealing algorithm for solving a
stochastic global optimization problem. More precisely, we address the problem of finding a global
minimizer of a function with noisy evaluations. We provide a rate of convergence and its optimized
parametrization to ensure a minimal number of evaluations for a given accuracy and a confidence level
close to 1. This work is completed with a set of numerical experimentations and assesses the practical
performance both on benchmark test cases and on real world examples.

1

Introduction

We are interested in an algorithm that solves the stochastic global optimization problem:
Find x? = arg min Eω (U (x, ω)),

(*)

x∈E

where x is a decision variable belonging to some large space E, ω is a random variable and U is the cost,
a positive and bounded real valued function. We do not make any assumption on the regularity of U . We
only expect it to be rapidly evaluable: typically, U is the result of some short numerical simulation. We
do not make any distribution assumption for the random inputs ω themselves but only on the outputs U .
We assume the code has some robustness property in the sense that, at some point x, it is either infinite
for all ω or bounded uniformly in ω.
This problem is twofold: we must both estimate and minimize the expectation of the cost. A simple and
general approach consists in the minimization of a sample average of Monte Carlo estimators:
N
X
b ω,N (U (x, ω)) := 1
U (x, ωi ),
E
N
i=1

for any given i.i.d. sample (ωi )1≤i≤N of size N distributed according to the distribution PΩ (x) of ω and
for any x ∈ E. Such an estimator consistently estimates Eω (U (x, ω)) for any given x. Nevertheless,
its accuracy is directly linked to N and thus to the computational effort. One can thus wonder if a
computationally efficient procedure using that estimator can return a solution to the initial problem
given a certain level of accuracy .
∗
†

clement.bouttier@airbus.com / clement.bouttier@math.univ-toulouse.fr
ioana.gavra@math.univ-toulouse.fr

1

1.1

Previous works: different types of algorithms

There were many attempts to solve this stochastic problem across several research communities. We give
a brief survey of them in what follows.
In the case E is finite and U takes its values in [0, 1], problem (*) is labelled as a ”simple regret bandit
optimization problem” by the bandit community. Indeed it can be seen as the problem of choosing, among
a small finite set of slot machines providing random rewards, the one with the best expected reward by
playing a minimal number of times. This is not the classical setting of bandit optimization which usually
seeks for the ”cumulative regret”. As algorithm proposals for the simple regret context often extends cumulative regret concepts [3], we focus on them first. The Upper Confidence Bound (UCB, [4]) algorithm
aimed at building sequences of confidence bounds around the estimated expected cost of each element of
the search space. If the space is too large this can be prohibitive. These were several attempts to bypass
this issue by adding some assumptions on the regularity of the cost function around its optimum. We can
mention HOO [9] that produced guarantees about the cumulative regret for a continuous Lipschitz cost
function with known Lipschitz constant. In the same framework StoSOO [23] relaxed this last assumption. Both algorithms were however not very efficient in practice if the search space is multidimensional.
Indeed they still required some uniform exploration of the state space in the first phase. This could
lead to numerical difficulties when the dimension was growing. The Adaptive-treed bandit algorithm [11]
partially solved this issue by adapting the exploration step using a Lipschitz constant per dimension.
These algorithms could all be viewed as stochastic variations around the classical branch and bound algorithm [22], which was extensively studied by the optimization community. We can mention the very
popular DIRECT algorithm [20], from which the StoSOO procedure was inspired.
Finally, let us mention the computer experiment community that introduced another popular global optimization method for dealing with the stochastic case, the so called Efficient Global Optimization (EGO,
[21]) based on expected improvement. The convergence rate of this method was already investigated in
[10] in a noise free context. This algorithm focused on minimizing the number of cost evaluations because
it considered a setting where the cost evaluations were very time consuming. As a result, in order to
select each evaluation point, it required a higher computational effort and memory storage per iteration
than other optimization methods. Such a method could therefore turn out to under-perform in a setting
where the computational cost ratio between selection and evaluation was inverted.
A typical algorithm that was known to perform well in the case of time-cheap cost evaluations was the
simulated annealing (SA) as mentioned by Locatelli in [19]: ”The latter algorithms (mainly EGO) often
outperform SA algorithms from the point of view of the number of function evaluations to reach a given
relative accuracy on the standard test functions from [12], but usually require a higher computational
effort per iteration. Typical advantages of SA algorithms are their very mild memory requirements and
the small computational effort per iteration. If the cost of a function evaluation is very high, then even
a considerable computational effort per iteration may be negligible with respect to the cost of a function
evaluation, and algorithms which require few function evaluations are preferable; otherwise, also the computational effort per iteration should be taken into account, and from this point of view SA algorithms are
often better than other algorithms.”
However SA algorithms have been designed and extensively studied in a context where the exact cost
could be observed. We recall below some basic facts in the noiseless case (Section 1.2) and then present
the noisy case which is the setting addressed in this paper (Section 1.3).

1.2

Simulated Annealing without noise

Let E be some finite search space and J : E → R+ a function that we want to minimize, called cost
thereafter.
Simulated Annealing is a classical global optimization method. It aims at building a sequence of
elements from E whose last element is drawn from a uniform probability law on the subset of global
2

minima of J. In other words it aims at sampling from the following distribution
µ? =

1Sopt
,
|Sopt |

where Sopt = {x, J(x) = miny∈E J(y)} and |.| denotes the cardinality of a set. Such a sampling is of course
not straightforward but one can notice that this distribution can be rewritten in the following form:
?

∀x ∈ E, µ (x) = lim P
T →0

e

−J(x)
T

y∈E e

−J(y)
T

,

P −J
−J
and it is well-known that the Gibbs distributions of the form µT = e T / e T are efficiently sampled
for reasonably low temperatures T ∈ R+ using the Metropolis-Hastings algorithm [1]. A quite natural attempt is therefore to build a sequence of sequences obtained using Metropolis-Hastings algorithm
for a set of decreasing temperatures. In particular, at a very low temperature, the Metropolis-Hastings
algorithm generates exploratory moves that are accepted with very low probabilities, which makes it a
very bad sampler. Therefore it is necessary to first encourage exploration by using a sampling at higher
temperatures. A lower bound on the temperature at each step ensuring a convergence in probability of
the algorithm has been provided by [16]. At the same time another proof of convergence using modern
semi-group representation of Markov processes has been obtained by [17]. The obtained bounds are less
explicit but contain information about the convergence rate and the proof scheme is much more general.
We set our work in the continuity of this last work and use similar notations.

1.3

Simulated Annealing with noisy evaluations

As mentioned previously, our main interest is to extend such a method of simulated annealing to the
stochastic case:
Find x? = arg min Eω (U (x, ω))
x∈E

where ω is a random input of a bounded cost function U whose expectation can only be numerically estimated through Monte Carlo simulations. In other words, we consider J(.) = Eω (U (., ω)). This question
is not novel and several attempts were made to address this problem theoretically in the 90’s. [14] were
probably the first ones to introduce the notion of simulated annealing with noisy measurements. They
assumed an additive Gaussian noise independent of the evaluation point and gave a sufficient condition
for the decrease of the variance σk2 of this noise, to ensure convergence of the algorithm to the optimal
set. [15] extended the results to distributions that are more peaked around zero than the Gaussian distribution. Their convergence result can be stated roughly as follows:
Theorem 1 ([15]). Let (Xk )k∈N denote the sequence of states in E visited by the simulated annealing
algorithm with Monte Carlo sampling of the noisy measurements. If:
(i) the convergence conditions from [16] are satisfied
(ii) ∃ > 0 such that the standard error of the noise at step k of the algorithm σk2 = O(k −(2+) )
then ∀x ∈ E,

lim P(Xk = x) = µ? (x), where µ? is the uniform distribution on the global minima of

k→+∞

the expected cost.



This result provided a first answer to our question about the convergence of the algorithm in the
stochastic case. However the convergence statement above did not give any information about the convergence rate of the algorithm. Following the noise-free proof of [1], [18] provided an extension of this
3

statement to the noisy case with bounded variance and introducing a state dependent noise. He obtained
the same constraint on the decrease of the variance and the same convergence statement. He also highlighted the need for an extended result concerning the rate of convergence and for numerical experiments.
Indeed on this second point we can mention the works of [13] and [8] that addressed this issue. [13]
made a very interesting proposition in the framework of Gaussian noise. He proposed to use the noise of
measurement to drive the simulated annealing, i.e., accept a move if the estimated cost of the proposed
solution is lower than the one of the current solution. Using an analogy with the Glauber acceptance
mechanism, which is a symmetric alternative to the Metropolis-Hasting mechanism [1], he proposed a far
more efficient criteria for the variance decrease, i.e., σk = O(log(k)−2 ). Unfortunately he only provided
a few numerical examples to validate his statement and a theoretical proof is still missing.

1.4

Main contributions

In this paper we consider a simulated annealing algorithm based on mini-batches of increasing size. More
precisely, at each iteration, the expected cost is estimated by Monte Carlo sampling
of increasing sizes.
b ω (U (xk , ω)) = 1/Nk PNk U (xk , ωi ), where Nk is
The estimated cost at step k of the algorithm is thus E
i=1
an increasing sequence and ωi are i.i.d. random variables having the same law as ω. The cost can be
b ω (U (xk , ω)) = Eω (U (xk , ω)) + ζω (xk ), where ζω (xk ) is some bounded random variable.
written also as E
2
We denote σk := Var(ζω (xk )), the variance of post-sampling noise. As it is directly linked to the number
of measurements made during the mini-batch, it can be tuned by the user.

Rate of convergence for all variances of polynomial decay. In the sequel we first show that
theoretical guarantees of Theorem 1 can be extended to sub-Gaussian random variables (e.g., bounded
noise distributions) with stronger convergence results for this algorithm. Indeed we show that convergence
can be ensured if the number of measurements is chosen such that σk = O(k −(α/2) ) with α > 0, which
corresponds to Nk of the order k α . One can observe that, as opposed to [15], the convergence still holds
for α ≤ 2. This is summarized in Theorem 2.
We derive the rate of convergence of the procedure (Theorem 4) and optimize it (Corollary 1) with
respect to the noisy simulated annealing algorithm parameters in order to provide a minimal total number
of measurements at given accuracy and confidence requirements. This leads to the optimal value α = 2
for which the number of cost evaluations increases fast enough to ensure almost the same convergence rate
as in the noise-free case. This shows that the convergence rate is limited by the concentration speed of the
Gibbs measure around its modes. According to our concentration result, increasing the estimation effort
cannot increase the performance of the algorithm above this limit. On the other hand the convergence
still holds for a decreased estimation effort (α < 2) as soon as the cooling schedule is slowed consequently.

Computational cost in the general case. Finally, we derive an upper bound on the computational time-complexity of our simulated annealing algorithm (with noisy measurements). This quantity
is roughly of the order of:
e

C1 log 1
δ


,

where C1 is some constant depending on the cost function itself as detailed in Corollary 1. The provided
bound exhibit an exponential dependency in 1/ and log 1/δ. This is comprehensive regarding the generality of the considered problem.

4

Computational cost in the absence of local minimum If the function has no local minimum apart
from the global minimum (e.g., a convex function evaluated on a finite set) the temperature schedule can
be adapted and the computational cost becomes of the order of:
C2 log 1δ


!3
,

where C2 is a constant detailed in Corollary 2. This second bound increases in a polynomial way with
respect to 1/ and log 1/δ. This is a very positive result as it shows that the noisy simulated annealing
algorithm recovers the state-of-the-art convergence guaranties if stronger hypotheses on the cost are
considered.
Numerical experiments. We provide numerical evidence indicating that the numerically observed
requirements in [13], i.e., σk = O(log(k)), do not hold for a Metropolis-Hastings Acceptance criteria. We
apply the noisy simulated annealing on classical non convex optimization test cases with different level of
noise, but also perform a test on a real-world example, i.e., an aircraft trajectory optimization problem
using a black-box aircraft performance model.

1.5

Aircraft trajectory optimization

As a leading example for this setting, we consider the problem of optimizing commercial aircraft trajectories with respect to a combination of fuel consumption and flight duration:.
Z

?

tf

Find u = arg min g(x(tf ), tf ) +
u

t0

−ṁ(x(s), u(s))ds

s.t. ∀t > t0 ẋ(t) = f (x(t), u(t))
x(t0 ) = x0

d(tf ) = df ,
where x is the state of the aircraft, m its mass, ṁ its instantaneous fuel consumption, d the ground
distance it has flown over, u the path control, f the instantaneous dynamic and g the terminal cost
function. The path control u is the combination of the thrust rating δT and the lift coefficient Cl .

u=

Cl
δT



Estimates of the cost of trajectories are usually obtained through numerical integration of the flight
dynamic equations, f :



ẋ = 



V̇
γ̇
ḣ
d˙
ṁ





 
 
=
 
 

1
(T (h, V, δT ) − D(h, V, CL )) m
− g sin γ
1
(L(h, V, CL ) − mg cos γ) mV
V sin(γ)
V cos(γ)
ηT (h, V )








where T is the thrust, D the drag, L the lift, η the specific fuel consumption, γ the path angle, V the speed
of the aircraft and h its height. These equations involve some terms like the aerodynamic drag coefficient
(CD ) or maximal propulsion effort (Tmax ) who are estimated using interpolation of experimental local
5

measurements.
T (h, V, δT ) = δT Tmax (h, V )
1
L(h, V, CL ) = ρ(h, V )SV 2 CL
2
1
D(h, V, CL ) = ρ(h, V )SV 2 CD (CL , V )
2
No analytic solution is therefore available nor conceivable. Moreover the relation between cost and
trajectory control parameters cannot reasonably be assumed to be convex. At last, the cost estimation
relies on some predicted flight conditions including atmospheric ones. Hence, real-flight costs can thus
deviate substantially from their predictions and some uncertainty propagation method must be applied to
obtain an accurate estimate of the expected flight costs. In other words the function we want to minimize
can only be evaluated with a certain random error, which corresponds exactly to the setting of this paper.
Finally, the computational efficiency is a key ingredient as it must be performed only a few hours before
the planned flight. For more information about aircraft trajectory optimization we refer to [6].
This example completely fits our requirement as the computation of the cost of one single complete
trajectory is quite fast, i.e., less than a second. Therefore, the EGO algorithm [21] would not be suited
for this application. On the other hand formulations based on the DIRECT algorithm [20] would suffer
strongly from the dimension of the problem. An additional element that motivates the use of simulated
annealing is the fact that in the case of trajectory optimization the set of admissible controls is not known
in advance as it is path dependent. We can only ensure that this set is connected. This implies in
particular that no projection on the constraints can be performed and excludes the projected stochastic
gradient descent for example. In the case of simulated annealing, a very simple step can bypass this
issue. By setting the value of the cost to infinity when the trajectory evaluator returns an error we ensure
staying in the admissible domain. Consequently, a feasible solution and a conservative approximation of
the admissible domain are the only requirements to initiate the algorithm in this setting.

1.6

Outline of the paper

Our paper is organized as follows. In Section 2 we present the noisy simulated algorithm and our main
theoretical result. In Sections 3, 4 and 5 we provide the proof of this statement. More precisely, in
Section 3 we compute the infinitesimal generator of the noisy simulated annealing algorithm. In Section
4 we compare it to the one of the noise-free simulated annealing algorithm from [17]. This enables us
to derive a differential inequality for a L2 distance between the distributions of the two previously mentioned processes. Integrating by applying Grönwall’s Lemma Section 5, we obtain obtain our convergence
result. In the same section, we show how to tune the parameters of the algorithm in order to optimize
the performance bound and give the corresponding computational cost. In Section 6 we propose some
numerical insight on synthetic and real data experiments.

Acknowledgments.
We thank Sébastien Gadat for introducing this topic to us, making this collaboration possible and for
fruitful discussions and helpful insights, and Sébastien Gerchinovitz for all his constructive advice and
ideas.

2

Noisy Simulated Annealing algorithm: statement and convergence
result

We first present our extended version of the simulated annealing to the stochastic case, whose pseudo-code
can be found in Algorithm 1.
6

2.1

Noisy Simulated Annealing algorithm (NSA)

Algorithm 1 Noisy Simulated Annealing
procedure NSA(Inputs: Neighbourhoods structure (Sx )x∈S , Initial guess x0 , increasing function
β : R+ → R+ , Function t 7→ nt )
Initialize time t0 = 0
β0 = β(t0 )
for k from 0 to Maximal number of iterations do
Draw one solution candidate: x
etk ∈ Sxtk according to q0 (xk , ·)
Draw Ntk ∼ Poisson(ntk ) + 1
Draw 2Ntk simulation conditions independently:
⊗Ntk
k ) ∼ (P (e
k ) ∼ (P (x ))⊗Ntk and (e
ω1k , ..., ω
eN
(ω1k , ..., ωN
Ω xtk ))
Ω tk
tk
tk
b t ) and J(e
b xt ) using the Nt conditions:
Compute estimates J(x
k
P
N
t
k ),
k
b t )= 1
J(x
U
(x
,
ω
t
k
k
i
i=1
Ntk
1 PNtk
b
xtk , ω
eik )
J(e
xtk ) = Nt
i=1 U (e
k

Draw an exponential random variable ξk+1 of parameter 1
Update time tk+1 := tk + ξk+1
With probability e−βk bJ(extk )−J(xtk )c+ :
set xtk+1 := x
etk
Otherwise set xtk+1 := xtk
b

b

Increase the inverse of the temperature βk+1 := β(tk+1 )
end for
return xtk+1
end procedure
where bxc+ = 0 if x ≤ 0 and x if not.
As in the deterministic setting, the algorithm requires an initial feasible solution xt0 , a temperature
schedule Tt (we will mostly use its inverse βt = 1/Tt ), and a good neighbourhood structure. What we
mean by good will be specified in the definition of q0 . The algorithm explores the state space in the
following manner. After k iterations, at time tk , it selects a random neighbouring solution x
etk ∈ Sxtk (Sxt
b xt )
being the set of neighbours of xtk ) according to a proposition law. Then it compares the estimate J(e
k
b t ) of the current solution and then it decides to
of the cost of this new solution to the estimated cost J(x
k
substitute (or not) the new to the current:
b xt ) ≤ J(x
b t ), the move
• if the estimated cost of the new state is lower than the current one, i.e., J(e
k
k
is accepted, i.e., xtk+1 ← x
etk
b xt ) − J(x
b t )).
• if not, it is only accepted with a probability exp(−βtk (J(e
k
k
The time t is then updated using independent exponential random variables, enabling us to consider the
NSA as a continuous time Markov process.
7

2.2

General setting and notations

To state the convergence of Algorithm 1, we first need to describe formally the framework we are working
in. Notations introduced in this section are valid for the whole paper unless mentioned explicitly.
• Regarding the noise structure and the estimation procedure, we denote:
P
b the estimated cost: Jb : E → R+ , such that ∀ x ∈ E, J(x)
b
(J)
= N1 N
i=1 U (x, ωi ), where:
(ω1 , ..., ωN ) is a N i.i.d. vectors sequence drawn from distribution PΩ (x)
(ξk ) the time increments: (ξk )k∈N is a sequence of i.i.d. exponential random variables of parameter
1
P
(tk ) the jumping times: ∀k ∈ N, tk = ki=1 ξi .

(nt ) the samplig intensity: nt a continuous increasing function.

(Ntk ) the sample sizes: Nt1 , Nt2 , . . . Ntn are independent for all n ∈ N and all 0 < t1 < t2 . . . < tn
and
Ntk ∼ Poisson(ntk ) + 1,
We can make a few remarks about the different notations. The construction of Nt ensures that
its value is a strictly positive integer at all times. The reason why we choose to have a randomly
sized sample for the Monte Carlo estimation procedure is rather technical. It enables generating a
continuous transition probability as it can be noticed in Equation (3) and ease the formulation of
the infinitesimal generator (Equation (7)).
• About the state space, we denote:
(E) a finite state space.
(S) a neighbourhood structure such that E is connected with respect to it, i.e., S is a connected
graph containing all the points in E. For any x in E, we denote Sx the set of its direct
neighbours.
(µ0 ) the initial distribution, a probability measure that charges every point of a subset of interest
E 0 ⊂ E defined more precisely in (U ) ,
(q0 ) the proposition law, an irreducible and µ0 − reversible transition probability, i.e., ∀x, y ∈ E,
∞
P
(n)
q0 (x, y) = ∞ and µ0 (x)q0 (x, y) = µ0 (y)q0 (y, x). In addition we assume that for any x in
n=0

E, we have q0 (x, Sx ) = 1
Considering a finite search space E enables us to easily derive the spectral gap inequality in Theorem
3 and overcome differentiation-under-the-integral-sign issues in Equation (15). Nevertheless, it could
be replaced by coercivity assumptions on the function J, which could be more general but not really
well suited for the application we are looking for. It is our most restrictive assumption. Nevertheless
it is in line with previous works on noisy global optimization for example: [15], [18] or [13]. It
corresponds to a historical use of simulated annealing for problems with huge finite search space like
for the traveling salesman problem [2]. Mimicking [17], we might however relax this assumption of
finiteness. Nevertheless it requires more technicalities as in [1] and this is left for future work.
We assume that the algorithm can visit and start from every point in the solution space through
the connection assumption S and the definition of µ0 . The proposition law q0 defines the way a
new solution x
e is proposed to the NSA at each iteration. The irreducibility of q0 implies the fact
that one can go from any state x to any other state y using the neighbourhood structure S, in a
8

finite number of steps. The µ0 − reversibility is used to simplify the notations. A classical choice
1
[1] for q0 and µ0 is: ∀ x, y ∈ E, µ0 (x) = |E|
and q0 (x, y) = |S1x | , assuming every point in E to
have the same number of neighbors. However there are other possible choices for µ0 and q0 . This
last two assumptions are inherited from the classical Metropolis-Hasting sampling algorithm which
corresponds to the NSA algorithm with no cooling mechanism and no noise. They ensure that a run
in this setting, starting from any point of the search space, converges to a stationary distribution
which is the Gibbs measure associated to J.
• About the cost function, we consider:
(U ) the underlying cost: ∃M > 0 and ∃E 0 ⊂ E, such that U is bounded and non-negative on E 0 , i.e.
∀x ∈ E 0 , ∀ω, 0 ≤ U (x, ω) ≤ M and U is infinite on E\E 0 , i.e.,∀x ∈ E\E 0 , ∀ω, U (x, ω) = +∞.
The assumption about U being bounded is not restrictive. It reflects the practical setting where a
simulation code crashes out of the definition domains. We associate infinite costs to crashes and
thus (U ) is rather a consequence of (E).
• About the algorithm parametrization, we denote
(βt ) the inverse of the temperature: a positive increasing real function of t,
(α) the sampling size (expected number of simulations): ∃α ∈ R+ such that nt = (t + 1)α ,
bd
t
βt is usually chosen such that ∀t ∈ R+ , dβ
dt = 1+td for some b, d ∈ R+ , as it was shown by [16]
and [17] to be a necessary condition to ensure the convergence of the simulated annealing algorithm
for any cost function. There is no reason to expect that the noisy context would be more favorable
than the deterministic one. As suggested by the definition of α, we choose a polynomial growth of
the number of simulations for the cost estimation. We show later on in this work that this ensures
the convergence of the noisy simulated annealing for a good choice of α and b.

2.3

Tool for the analysis: the NSA process

We now present the mathematical formalization of the NSA algorithm’s underlying stochastic process.
First, for pedagogical purposes, we omit the temperature evolution and noisy measurements. The NSA
algorithm then becomes a simpler Markov chain exploring the state space E according to the Markovian
transition matrix whose elements are of the form:

q0 (x, y)e−βbJ(y)−J(x)c+ if y 6= x
(1)
P(x → y) = qβ (x, y) = 1 − P q (x, z) if y = x,
β

z∈E\x

This reflects the transition mechanism introduced at the beginning of this section. As the process is in
fact a continuous one, we must also consider the time component. NSA jumps happen at stochastic times
and the probability of acceptance depends on these times. Combining the law of the jumping times and
the previous mechanism, we can make their joint transition probability explicit:
• Let (e
χk , Tk )k∈N a E × R+ -valued Markov chain such that ∀k ∈ N, ∀y ∈ E, ∀u ∈ R+ :
P(e
χk+1 = y, Tk+1

+∞
Z
≥ u|e
χk , Tk ) =
qeβτ (e
χk , y)1[Tk ,+∞[ (τ )e−(τ −Tk ) dτ,

(2)

u

where
 β PN


t U (y,ω )−U (x,ω )c
− Nt b i=1
i
i +

t
q0 (x, y)ENt Eω1 ,...,ωNt e
qeβt (x, y) =
P

qeβt (x, z) if y = x
1 −
x6=z

9

if y 6= x

(3)

This is a similar construction to the one of the classical simulated annealing process [16]. The state
transition mechanism must also reflect the estimation procedure, therefore the form of Equation
(3) differs from Equation (1). As mentioned before the function t 7→ βt represents the inverse of
the temperature schedule and Nt is the random process described by Ntk . The jumping times, or
evaluation times of the process happen at times defined by the sequence Tk (cf. the definition of
tk ).
The chain (e
χk )k≥0 explores the state space E using a transition probability qeβt constructed in a
same way as the classical one, replacing the exact value of −βTk (J(y) − J(x))+ by its Monte Carlo
estimation. The expected value from the formula comes from the fact that, as mentioned in the
definition of Jb and in Algorithm 1, we use a random number of Monte Carlo shootings for the
estimations.
Finally, we obtain the NSA process by associating the two sub-processes as follows:
 
et = χ
et
ek if Tk ≤ t < Tk+1 . One
be the inhomogeneous Markov Process such that X
• Let X
t≥0

can see that this process is piecewise constant and jumps at exponential times from one candidate
et )t≥0 is just the continuous-time version of the noisy simulated
solution to another, in other words (X
annealing discrete time process,(e
χk )k≥0 .

Note that, if y ∈ E\E 0 then ∀x ∈ E 0 ,
et ∈ E 0 .
E 0 , then ∀t ≥ 0, X

2.4

e0 is chosen in
qeβt (x, y) = 0. Hence, if the initial solution X

Convergence result

We denote:
• m? the maximum depth of a well not containing a fixed global minimum of the function J. To be
more precise, we call a path from x to y any finite sequence x0 = x, x1 , . . . , xn = y such that for all
i, xi+1 ∈ Sxi . Let Pxy be the set of paths from x to y.
For a given path p ∈ Px,y , the elevation of the function J on p is max J(z). Minimizing this quantity
z∈p

over the set of possible paths Px,y , gives us the elevation of the cheapest path going from x to y.
Denote this elevation by:


Hx,y = min max J(z)
p∈Pxy

z∈p

Then
m? := max {Hx,y − max (J(y), J(x))}
x,y∈E

(4)

As represented on Figure 1, m? can also be understood as the highest energy barrier to climb to
go from one point to another in the search space in the easiest direction. As mentioned before, it
also represents the maximal depth of a well not containing a fixed global minimum. If x? is a global
minimum then:
m? = max {Hx? ,y − J(y)} .
y∈E

The definition provided here is equivalent to the classical one, i.e., the one provided in [16] and [17].
A proof of this statement can be found in Appendix B.
• γ(β) the spectral gap between 0 and the rest of the L2 (µβ ) spectrum of −Lβ , where Lβ is the
generator of the classical simulated annealing (for more details about Lβ see Section 3):
 Z

Z
Z
2
γ(β) := inf − φLβ φdµβ s.t.
|φ| dµβ = 1 and
φdµβ = 0
10

(5)

p ∈ Pyx?
J

min maxJ(z)

p∈Pyx? z∈p

m?

y

x? = argminJ(u)
u∈E

Figure 1: m? , maximal depth of local minima
Following [17], we know that given E, µ0 and U , there exists a constant c such that:
∀β ≥ 0,

γ(β) ≥ ce−m

?β

Remark that this lower bound is mainly informative for small values of β. In addition set:
∀x ∈ E, J(x) = Eω (U (x, ω)) and J ? = min Eω (U (x, ω)).
x∈E

We define χ the set of -optimal points in E, i.e.,
χ = {x : J(x) ≤ J ? + },

(6)

and denote c χ = E\χ , its complementary in E. We also write a ∧ b = min a, b.
Theorem 2. Consider the settings of Section 2.2,
if βt = b log(td + 1) and nt = (t + 1)α , with:
{m? b < 1 ∧ α/2} or {m? b = 1, α > 2 and d < 2cm? /M },
et ∈ c χ ) ≤ C(µβ (c χ ))1/2 .
P(X
t

then there exits C > 0 such that, ∀t ∈ R+ , ∀ > 0,



This theorem is a natural extension of the result provided by [17]. There are two main interesting
facts to point out. First, we obtain a balance between the expected number of Monte Carlo simulations
at each step of the algorithm and the inverse of the temperature, i.e., {m? b < 1 ∧ α/2} or {m? b = 1, α >
2 and d < 2cm? /M }. Reducing the growth rate α of the number of simulations below the quadratic rate
should be compensated by decreasing accordingly the temperature factor b. Second, the convergence is
stated in terms of a bound on the probability of not returning an optimal solution. Using the concentration speed of the Gibbs measure one can deduce a rate of convergence of the algorithm. Also the theorem
provides an insight on how the algorithm could be used in practice. A run of parallel noisy simulated
annealing would have a probability of returning a bad solution that would decrease in the power of the
number of runs. Nevertheless this benefit should be traded with an additional selection cost. Indeed, if
1
we obtain K solutions retrieved by K parallel NSA realization, we still face the problem of selecting the
best one. We only access estimates of the costs associated to each solution.

11

Sketch of the proof The proof of this theorem is divided into three parts. First, in Section 3, we compute the infinitesimal generator of the classical (Equation (11)) and noisy simulated annealing (Equation
(7)). Second, in Section 4, we compare them (Lemma 1) and third, in Section 5, we conclude about the
convergence using the Grönwall lemma (Equation (18)) and the convergence of the classical simulated
annealing (Equation (24)).
Convergence rate In the case m? b < 1 a finer bound can be deduced from Grönwall’s lemma and one
can obtain a more precise convergence rate for the algorithm (Theorem 4), which is roughly of the order
of:
et ∈ c χ ) ≤ Γ t((m
P(X

? −)b−min(1,α/2))/2

where Γ is some constant detailed in Theorem 4. In particular this implies that for fixed , δ > 0 we can
e ? ∈ c χ ) ≤ δ. This leads to a bound (Lemma 3) on the computational complexity,
find T ? such that P(X
T
T
∗
E Ncall , of the order of:
 2(α+1)/(min(1,α/2)−(m? −)b)

Γ
T∗
E Ncall ≤
.
δ

3

Proof, Part 1: Infinitesimal generator

In this section we use the semi-group characterization of the generator in order to prove that as soon
e β of the
as qeβt defined in Equation (3) is continuous with respect to t then the infinitesimal generator L
t
et can be written as:
Markov process X

X
e
(7)
f (y) − f (x) qeβt (x, y).
Lβt f (x) =
y∈E

We briefly recall the definition of the semi-group associated to a Markov process.
Definition 1. The semi-group (Pt,t+s )t≥0,s≥0 associated to the Markov process (Xt )t≥0 is a family of
probability kernels such that for all non-negative borelian functions:
∀t, s ∈ R+

Pt,t+s f (x) = E(f (Xt+s )|Xt = x)

Let (Pt,t+s )t≥0,s≥0 be the semi-group associated to the Markov process (Xt )t≥0 . The semi-group
characterization of its generator is given in the following definition:
Definition 2. The infinitesimal generator Lt of the Markov process (Xt )t≥0 is defined as the operator
such that for any bounded function f :
Lt f (x) = lim

s→0

Pt,t+s f (x) − Pt,t f (x)
s

We start by computing the infinitesimal generator Lβt of the process associated to the SA algorithm,
i.e., with no measurement noise, and then deduce the infinitesimal generator of the NSA algorithm. Using
similar notations to the ones of Section 2.3, we consider the noise free inhomogeneous Markov process,
(Xt )t≥0 constructed from the inhomogeneous Markov chain (χk )k∈N whose one step transition probability
is:

q0 (x, y)e−βTk (J(y)−J(x))+ if y 6= x
∀x, y ∈ E, qβTk (x, y) = 1 − P q (x, z) if y=x
βTk

z∈E\{x}

This is the natural extension of the simulated annealing process with discrete jumping times [16] to
the continuous time process. In this configuration, the jumping times are drawn from an i.i.d. sequence
12

of exponential random variables of parameter 1. In the homogeneous configuration, i.e., βt = β, the
infinitesimal generator has a classical form: Lβ = Qβ − Id where Qβ is the transition matrix associated
to qβ and Id denotes the identity. The extension to the generator of the non-homogeneous process is not
straightforward. Therefore we propose to detail the computations.
By definition, for any bounded function f :
Pt,t+s f (x) − Pt,t f (x)
s→0
s
P
f (y)P(Xt+s = y|Xt = x) − f (x)

Lβt f (x) = lim
= lim

y∈E

s→0

P
= lim

f (y)P(Xt+s

y∈E

s
= y, Ht+s − Ht ≥ 0|Xt = x) − f (x)
s

s→0

,

(8)

where Ht = max{k ∈ N : Tk < t} denotes the number of jumps before time t. Since Tk is a sum of
independent exponential variables of parameter 1, one can remark that Ht is in fact a Poisson process of
parameter 1.
In order to compute the above limit, we begin by calculating a more explicit form of the probabilities
above. We can divide these computations into three parts according to the number of jumps between t
and t + s:
P(Xt+s = y, Ht+s − Ht ≥ 0|Xt = x) = P(Xt+s = y, Ht+s − Ht = 0|Xt = x)

+ P(Xt+s = y, Ht+s − Ht = 1|Xt = x)

+ P(Xt+s = y, Ht+s − Ht ≥ 2|Xt = x).
The first case is straightforward, if there is no jump between t and t + s, the process will not change its
position and we thus have:
P(Xt+s = y, Ht+s − Ht = 0|Xt = x) = δx (y)e−s .
The second case is slightly more involved. Using the stationarity and the definition of Poisson processes,
the event that the algorithm goes from x to y, having only one jump between t and t + s, can be written
as:
P(Xt+s = y, Ht+s − Ht = 1|Xt = x) =P(Xt+s = y, ξ10 < s, s − ξ10 < ξ20 |Xt = x)
where ξ10 and ξ20 are two independent exponential random variables of parameter one.
Let ξ = (ξ10 , ξ20 ) and Ds = {(h1 , h2 ) ∈ R2 |h1 < s and h2 > s − h1 }. Also in what follows, for a random variable Y we denote fY its probability distribution. Using these notations and the fact that ξ is
independent of Xt , we can write:
Z
P(Xt+s = y, ξ ∈ Ds |Xt = x) =
f(Xt+s ,ξ)|Xt =x (y, h)dh
Ds
Z
=
fXt+s |ξ=h,Xt =x (y)fξ (h)dh
Ds
s Z +∞

Z
=

0

13

s−h1

qβt+h1 (x, y)e−h1 e−h2 dh1 dh2

The previous equality yields:
P(Xt+s = y, Ht+s − Ht = 1|Xt = x) = e

−s

Z
0

s

qβt+h1 (x, y)dh1 .

(9)

In the following we use the classical O(.) and o(.) notations: for all functions f and g defined on some
subset of R,
• f (x) = O(g(x)) as x → 0+ ⇐⇒ ∃σ, x0 > 0, |f (x)| ≤ σ|g(x)| for all 0 < x ≤ x0
• f (x) = o(g(x)) as x → 0+ ⇐⇒ lim

x→0+

f (x)
g(x)

= 0.

For the third term we can see that:
P(Xt+s = y, Ht+s − Ht ≥ 2|Xt = x) ≤ P(Ht+s − Ht ≥ 2) ≤ P(Hs ≥ 2)
Since Hs is a Poisson Process of parameter 1, one can check that for all s close to zero we have that
P(Hs ≥ 2) = 1 − P(Hs = 0) − P(Hs = 1) = 1 − e−s − se−s = O(s2 ).
This implies that when s is close to zero, the probability that the process goes from x to y between t
and t + s, with more than one jump is small in comparison to s:
P(Xt+s = y, Ht+s − Ht ≥ 2|Xt = x) = O(s2 )

(10)

Putting all the terms together and replacing them in Equation (2), we can rewrite the infinitesimal
generator as follows:
1
Lβt f (x) = lim
s→0 s

X


Z
−s
−s
f (y) δx (y)e + e
0

y∈E

s





qβ(t+τ ) (x, y)dτ − f (x)

1X
f (y)P(Xt+s = y, Ht+s − Ht ≥ 2|Xt = x)
s→0 s

+ lim

y∈E

Using the fact that f is bounded, E finite and the upper bound given by Equation (10), one can easily
check that the second term is zero. Hence we obtain:




Z s
X
1  −s 
Lβt f (x) = lim
e
f (x) +
f (y)
qβ(t+τ ) (x, y)dτ  − f (x)
s→0 s
0
y∈E


Z s
f (x)(e−s − 1)
e−s  X
= lim
+ lim
f (y)
qβ(t+τ ) (x, y)dτ 
s→0
s→0 s
s
0
y∈E

Noting the fact that qβt is continuous with respect to t and the following identity
e−s = 1 − s + O(s2 ),
we easily obtain the simplest form for the infinitesimal generator of the inhomogeneous Markov chain:

X
Lβt f (x) =
f (y) − f (x) qβt (x, y).
(11)
y∈E

We can remark that the explicit form of the transition probability qβt does not appear in the proof, hence
the result is completely general. The only necessary property of this transition probability is its continuity
with respect to t.
The fact that nt and βt are continuous functions ensures the continuity of transition probability qeβt ,
defined in Equation (3) . Therefore, following the same argument,one can deduce (7). Here we can see the
relevance of the randomness of Nt . An increasing deterministic sequence would generate a discontinuous
qeβt and would make difficult the use of derivations above.
14

4

Proof, part 2: Generators comparison

The fact that for a temperature schedule that decreases slowly enough, the process generated by the
classical Simulated Annealing converges to the set of global minima of J is well known. The Noisy
Simulated Annealing is a similar algorithm, built on the same principles except that the values of the
function J are replaced by an estimation each time its computation is needed. Therefore a tight relation
exists between both approaches. Furthermore, as we will show in this section, for a well chosen couple
(βt , nt ) the generators of the two algorithms will be ’close’ at large times. This is a key element of the
proof as it will imply a first condition for the ratio βt /nt .
Using the relations given by Equation (11) and Equation (7), the quantity of interest is:
X
e β f (x) = Lβ f (x) +
L
(f (y) − f (x))(qf
βt − qβt )(x, y).
t
t
y∈E

Hence quantifying the difference between the two generators can be reduced to bounding the difference
between the two probability transitions qβt and qeβt . Thus the main result of this section is the following
lemma.
√
+
Lemma 1. Let βt / nt → 0. There exist two functions −
t and t such that
∞

qβt − qβt )(x, y) ≤ +
−
t qβt (x, y)
t qβt (x, y) ≤ (e

∀ t ∈ R+ , ∀ x ∈ E 0 , ∀ y ∈ E,
and
lim −
t→+∞ t

= lim +
t = 0.
t→+∞

Before going into the proof of this lemma, we present some preliminaries. First it can be noticed that
for all x, y ∈ E, x 6= y we have:
(qf
βt − qβt )(x, y)


qf
βt
= qβt (x, y)
− 1 (x, y)
qβt



− Nβ b

PNt

U (y,ωi )−U (x,ωi )c+



e
= qβt (x, y) ENt Eω1 ,...,ωNt  −β bE(U (y,Ω))−E(U (x,Ω))c − 1 .
+
e t
t

i=1

Unless specified otherwise, in this section we always consider x 6= y . The case x = y is handled at the
end of the section. To simplify the notations we denote Xix,y := U (y, ωi )−U (x, ωi )−E(U (y, ωi )−U (x, ωi ))
and K x,y = E(U (y, ωi ) − U (x, ωi )). Hence,
(qf
βt − qβt )(x, y)





P Nt
x,y
x,y c +bK x,y c
−βt b N1
+
+
i=1 Xi +K
t
−1
.
= qβt (x, y) ENt Eω1 ,...,ωNt e
Noticing that,
∀ a, b ∈ R,

−|a| ≤ −ba + bc+ + bbc+ ≤ |a|,

we obtain the following bounds for (qf
βt − qβt )(x, y):

ENt Eω1 ,...,ωNt

β

| Nt

e

t

PNt

i=1

Xix,y |


−1

qf
β − qβt
≥( t
)(x, y) ≥ ENt Eω1 ,...,ωNt
qβt

 β PN

t X x,y |
−| Nt
i=1
i
e t
−1 .

(12)

In order to obtain a bound for the expectation of a function of Nt we need an estimation of the probability
that Nt takes values ’far’ from its expectation.
15

Lemma 2. There exist δ ∈ (0, 1) and a = |(1 − δ) (1 − log(1 − δ)) − 1| such that for all t > 0 we have:
P(Nt ≤ (1 − δ)nt ) ≤ e−ant .

Proof. We remind the reader that, as mentioned in the definition Ntk , at a fixed time t the process
Nt can be written as 1 + H, where H is a Poisson random variable of parameter nt .
Fix t ∈ R and δ ∈ (0, 1). We provide an upper bound for P(Nt ≤ (1 − δ)nt ) using the Cramer-Chernoff
method.
First we see that for any λ > 0, applying Markov’s inequality we have:
ENt [e−λNt ]
..
eλ(δ−1)nt
For t fixed Nt − 1 has the distribution of a Poisson random variable of parameter nt . Therefore by
direct computations we have:
X
nk
ENt [e−λNt ] =
(e−λ(k+1) e−nt t )
k!
P(Nt ≤ (1 − δ)nt ) = P(e−λNt > eλ(δ−1)nt ) ≤

k>0

= e−nt −λ

X (e−λ nt )k
k!

k>0
−λ n
t

= e−nt −λ ee
Putting all these elements together yields:

P(Nt ≤ (1 − δ)nt ) ≤ exp([−λ(δ − 1) − 1 + e−λ ]nt )

The idea is to choose λ and δ in order to obtain the smallest possible value for −λ(δ − 1) − 1 + e−λ . For
δ ∈ (0, 1), the minimum is reached at λ = − log(1 − δ) which is strictly positive. For such λ and δ, we
denote a = |(1 − δ) (1 − log(1 − δ)) − 1| and conclude the proof.

Now we have what we need in order to start the proof of Lemma 1.
Proof of Lemma 1.
First, recall that if y ∈ E\E 0 then ∀x ∈ E 0 , qeβt (x, y) = qβt (x, y) = 0. Hence Lemma Equation (1) is
trivially verified for x ∈ E and y ∈ E\E 0 .
Considering the inequalities given by Equation
sep- 
 β PN(12), the
 proof can be divided in two parts by
 studying
βt PNt
t X x,y |
| Nt
X x,y |
−|
arately the upper bound ENt Eω1 ,...,ωNt e t i=1 i
and the lower bound ENt Eω1 ,...,ωNt e Nt i=1 i
.

Upper bound

We have:

(qf
βt − qβt )(x, y) ≤ qβt (x, y) ENt Eω1 ,...,ωNt

 β PN

t X x,y |
| Nt
i=1
i
−1
.
e t

 β PN

t X x,y |
| t
First we will provide an estimate of Eω1 ,...,ωNt e Nt i=1 i |Nt for all Nt .
We start by rewriting this expectation as:
 β PN
 Z
t X x,y |
| t
Eω1 ,...,ωNt e Nt i=1 i |Nt =

 β PN

t X x,y |
| t
P e Nt i=1 i > u|Nt du

R+

Z
=

P
R+

Nt
X
i=1

Xix,y

log(u)Nt
>
|Nt
βt

16

!
du.

Since Xix,y = U (y, ωi ) − U (x, ωi ) − ENt Eωi (U (y, ωi ) − U (x, ωi )) is a centered random variable and U is
bounded on E 0 (see the definition of U ) there exists σ such that |Xix,y | ≤ σ (for example set σ = 2M ),
almost surely for all i. Therefore (Xix,y )1≤i≤Nt are sub-gaussian random variables ([7]) with variance
factor σ 2 , i.e.,
u2

∀u ≥ 0, max (P(Xix,y > u), P(−Xix,y > u)) ≤ e− 2σ2 .

Considering that (Xix,y )i≤Nt is a sequence of independent sub-Gaussian
 their sum is still a
P variables,
x,y
x,y
Nt
2
≤ Nt σ 2 and therefore:
sub-Gaussian variable. As Var(Xi ) ≤ σ for all i we have that Var
i=1 Xi
Nt
X

∀u ≥ 0, max P

i=1

!
Xix,y > u , P −

Nt
X

!!
Xix,y > u

i=1

≤e

−

u2
2σ 2 Nt

.

For more details about sub-gaussian variables we refer to [7]. We use this property of concentration
in order to get an estimate of the expectation:

Eω1 ,...,ωNt

 β PN

x,y
t
| Nt
i=1 Xi |
t
e
|Nt
Z
=

Nt
X

1

P
0

≤1+2

Xix,y

i=1

Z

+∞ −

e

1
2σ 2

log(u)Nt
>
|Nt
βt


√

log(u) Nt
βt

!

Z

+∞

du +

P

Nt
X

1

Xix,y

i=1

log(u)Nt
|Nt
>
βt

!
du

2

du.

1

Using a simple variable substitution s = log(u)
with λt = √σβNt , we get:
λt
t
 β PN

Z
+∞
t X x,y |
s2
| t
Eω1 ,...,ωNt e Nt i=1 i |Nt ≤ 1 + 2λt
eλt s e− 2 du
0
Z +∞
λ2
(s−λt )2
t
e− 2 ds
≤ 1 + 2λt e 2
0
Z +∞
λ2
u2
t
e− 2 du
≤ 1 + 2λt e 2
−λt

√

λ2
t

≤ 1 + 2 2πλt e 2 P(G > −λt )
√
λ2
t
≤ 1 + 2 2πλt e 2 (1 − P(G > λt )).
where G is a standard Gaussian. Thanks to the Taylor formula, we know that there exists a constant
0 < θ < 1, such that:
e

−(θλt )2
2

P(G > λt ) = P(G > 0) + λt √
=

2π

−(θλt )2
2

1
e
+ λt √
2
2π

−(λt )2

1
e 2
≥ + λt √
.
2
2π
This leads to:


β

| Nt

Eω1 ,...,ωNt e

t

P Nt

i=1

Xix,y |

√
λ2
t
|Nt ≤ 1 + ( 2πe 2 )λt .

17



Thus, replacing λt by its definition we see that we need an estimate of:
√



ENt 1 + ( 2πe

σ 2 βt2
2Nt


σβt
.
)√
Nt

In order to simplify the notations we denote: gt = σβt . Using the bound given by Lemma 2 we get:
"
ENt 1 +

√

g2

√t

2πe 2

= 1 + ENt
≤ 1 + ENt
≤1+

√

Nt

#

g
√t
Nt

g2

"
√

2πe 2

"
√

2πe

√t

Nt

g2



2

√t

Nt

gt2

#

g
√t
Nt

g
√ t (1[1,(1−δ)nt ] + 1[(1−δ)nt ,+∞) )
Nt


#

e 2(1−δ)nt
2
2πgt  p
+ egt −ant 
(1 − δ)nt

Finally we have obtained that under assumptions of Section 2.2:


 β PN

t X x,y |
| Nt
i=1
i
(qf
−1
βt − qβt )(x, y) ≤ qβt (x, y) ENt Eω1 ,...,ωNt e t
√



gt2
2(1−δ)nt



e
2
+ egt −ant 
≤ qβt (x, y) 2πgt  p
(1 − δ)nt
Hence it is natural to define +
t as:

+
t =

√


βt2 σ 2
2(1−δ)nt
e
2
2
2πβt σ  p
+ eβt σ −ant  .
(1 − δ)nt


√
Since βt / nt → 0 one can check that +
t goes to 0 when t goes to infinity. Here we can see once more the
∞
importance of the balance between the two parameters βt and nt .

Lower bound

Considering the left-hand side of Equation (12) we have:

 β PN

t X x,y |
−| Nt
i=1
i
t
(qf
−1
.
βt − qβt )(x, y) ≥ qβt (x, y) ENt Eω1 ,...,ωNt e


In a similar way we start by obtaining a lower bound for Eω1 ,...,ωNt
improve it using the probabilistic properties of Nt . First observe that:
18

β

−| Nt

e

t

P Nt

i=1

Xix,y |


|Nt

and after we



β

−| Nt

Eω1 ,...,ωNt e

t

P Nt

i=1

Xix,y |


|Nt =

Nt
X



Z

P Nt

β

−| Nt

P e

i=1

t

Xix,y |


> u|Nt du

R+

!
log(u)Nt
=
P −
>
|Nt du
βt
R+
i=1
!
Z 1
Nt
X
− log(u)Nt
x,y
P
Xi
<
=
|Nt du
βt
0
i=1
!
Z 1
Nt
X
− log(u)Nt
x,y
1−P
Xi
>
=
|Nt du
βt
0
i=1


Z 1
N
− log(u) 2
− t2
βt
2σ
e
≥1−2
du
Z

Xix,y

0

Again, this is due to the fact that the sum of Xix,y is sub-Gaussian with variance factor Nt σ 2 .
with λt = √σβNt , we get:
Using the same variable substitution as above: s = log(u)
λt
t
 β PN

Z
0
t X x,y |
s2
−| t
Eω1 ,...,ωNt e Nt i=1 i |Nt ≥ 1 − 2λt
eλt s e− 2 du
−∞

≥ 1 − 2λt e
≥ 1 − 2λt e

λ2
t
2

λ2
t
2

Z

0

e−

−∞
Z −λt

(s−λt )2
2

e−

u2
2

ds

du

−∞

√
λ2
t
≥ 1 − 2 2πλt e 2 P(G > λt )
where G is a standard N (0, 1) Gaussian. As seen before there exists some 0 < θ < 1, such that:
−(λt )2

e 2
1
.
P(G > λt ) ≥ + λt √
2
2π
This leads to


β

| Nt

Eω1 ,...,ωNt e

t

P Nt

i=1

Xix,y |

√
λ2
t
|Nt ≥ 1 − ( 2πe 2 )λt .


This expression has exactly the symmetric form to the one obtained in the upper bound part. Thus, the
lower bound is obtained the same way as the upper bound. We directly get:


gt2
√
gt
ENt 1 − ( 2πe 2Nt ) √
Nt


gt2
2(1−δ)nt
√
e
2
≥ 1 − 2πgt  p
+ egt −ant 
(1 − δ)nt
Now we define −
t :


βt2 σ 2
2(1−δ)nt



√
2 2
 e
−
+ eβt σ −ant  .
t = − 2πβt σ p
(1 − δ)nt

√
It is easy to see that −
0. This completes the
t goes to 0 when t goes to infinity as soon as βt / nt →
∞
proof of Lemma 1.


19

5

Proof, last part: rate of convergence in the general case

We first complete the proof of convergence as stated in Theorem 2 and then deduce the convergence rate
(Theorem 4) from it. This enables us to provide an upper bound on the minimal number of cost function
evaluations in Section 5.3.

5.1

Proof of Theorem 2

The proof of Theorem 2 follows the roadmap of Holley and Strook [17] and relies on the use of the
Grönwall lemma. We derive a differential inequality for the L2µβ -norm of the density measure of the NSA
t
process with respect to µβt and deduce an integrated version of it using the lemma. We then show that
bounding the L2µβ -norm of this density implies the convergence of the process to the optimal state space χ .
t

Proof. [Proof of Theorem 2] Our goal is to show that when t goes to infinity, the noisy simulated
annealing gets “close enough” to the classical simulated annealing. Therefore we denote by ft the Radonet with respect
Nikodym derivative of the probability density of the noisy simulated annealing process X
to the Gibbs measure µβt , i.e.:
dmt
ft =
(13)
dµβt
es )s≥0 at time t. A first remark is that ft (x) = 0 for all t ≥ 0 and all
where mt is the distribution of (X
0
x ∈ E \ E , since our process, by construction does not accept states out of E 0 .
e β is continuous and therefore the
Using the results obtained in Section 3 one can see that R+ 3 t → L
t
semi-group (Ps,t )0≤s≤t is smooth. Also by their definition the operators (Ps,t )0≤s≤t are linear and have
the following semi-group property: Ps,t+h = Ps,t ◦ Pt,t+h , for all 0 ≤ s < t and h > 0. Hence, for all
0 ≤ s ≤ t, we have:
d
Ps,t = Ps,t Lt .
(14)
dt
For details about the infinitesimal generator see Section 1.4. of [5].
As shown in Equation (23), bounding the L2 -norm of ft w.r.t. µβt , i.e., kft kµβt , ensures convergence
of the NSA algorithm. However it does not provide enough information about the convergence of mt to
µβt to deduce a fine convergence rate. This is why we study the evolution of kft − 1kµβt which controls
the distance between the two measures. If this quantity is bounded then we obtain the convergence of
the NSA algorithm. If moreover it converges to zero, it implies a stronger convergence rate. In order to
prove that, we deduce a differential inequality for kft − 1k2µβ . We start by computing its derivative:
t
X
∂t kft − 1k2µβ = ∂t kft k2µβ =∂t
ft2 (x)µβt (x)
t

t

x∈E

=2

X


ft (x)∂t

x∈E


X
mt
(x)µβt (x) +
ft (x)∂t µβt (x).
µβt
x∈E

Using the backward Kolmogorov equation given by Equation (14), for the first term we have:




X
X
X
mt
mt
(x)µβt (x) =
ft (x)∂t
(x)µβt (x) +
ft (x)∂t µβt (x)
ft (x)∂t
µβt
µβt
x∈E
x∈E
x∈E
i
Xh
X
mt
e β ft (x) mt (x) −
=
L
ft (x)
(x)∂t µβt (x).
t
µβt
x∈E
x∈E
R
Denote hJiµβt := Jdµβt the mean of J with respect to µβt . One can check that:
h
i
∂t µβt (x) = −βt0 J(x) − hJiµβt µβt (x).
20

(15)

Thus, we easily obtain the following equality:
X
X
e β ft )(x)µβ (x) + βt0
∂t kft − 1k2µβ = 2
ft (x)(L
(J(x) − hJiµβt )ft2 (x)µβt .
t
t

(16)

t

x∈E

x∈E

First, we focus on the first term of the right hand side of Equation (16). Since we try to control the
eβ
generator of the noisy simulated annealing by the generator of the classical one, it is natural to write L
t
e β − Lβ . This comparison leads to the computation:
as Lβt + L
t
t
X

e β ft )(x)µβ (x)
ft (x)(L
t
t

x∈E


=

X

ft (x)(Lβt ft )(x)µβt (x) +

x

X

ft (x) 

x∈E


X
y∈E


(ft (y) − ft (x))(qf
βt − qβt )(x, y) µβt (x).

We rewrite the last part of the second term using Lemma 1:
X
(ft (y) − ft (x))(qf
βt − qβt )(x, y)
y∈E

=

X
y∈E

≤ +
t
≤ +
t
≤

ft (y)(qf
βt − qβt )(x, y) −

X
y∈E

X

ft (y)qβt (x, y) − +
t

X
y∈E

X

ft (x))(qf
βt − qβt )(x, y)

ft (x)qβt (x, y) + +
t

X
y∈E

y∈E

ft (x)qβt (x, y) − −
t

X

ft (x)qβt (x, y)

y∈E

−
(ft (y) − ft (x))qβt (x, y) + (+
t − t )ft (x)

y∈E
+
t Lβt ft (x)

−
+ (+
t − t )ft (x).

Inserting this in the previous inequality, we get:
X
e β ft )(x)µβ (x)
ft (x)(L
t
t
x∈E

≤ (1 + +
t )

X

+
t )

X

≤ (1 +

x

−
ft (x)(Lβt ft )(x)µβt (x) + (+
t − t )

ft (x)(Lβt ft )(x)µβt (x) +

(+
t

x

−

−
t )

X

ft2 (x)µβt (x)

x∈E

X
x∈E

−
(ft2 (x) − 1)µβt (x) + (+
t − t ).

Therefore, using Equation (16), we obtain the following inequality:
d
kft − 1k2µβ
t
dt "

#

≤ 2 (1 + +
t )
+

βt0

X
x∈E

X
x

−
ft (x)(Lβt ft )(x)µβt (x) + (+
t − t )

(J(x) −

X
x∈E

−
(ft2 (x) − 1)µβt (x) + (+
t − t )

hJiµβt )ft2 (x)µβt

"

#

≤ 2 (1 + +
t )
+

βt0

X
x∈E

X
x

−
ft (x)(Lβt ft )(x)µβt (x) + (+
t − t )
2

(J(x) − hJiµβt )(ft (x) − 1) µβt (x) +
21

X
x∈E

2βt0

−
(ft2 (x) − 1)µβt (x) + (+
t − t )

X
x∈E

(J(x) − hJiµβt )(ft (x) − 1)µβt (x) .

In order to deal with the first sum we use an estimate of the spectral gap of Lβt . This is provided by
Theorem 2.1 of Holley and Strook[17].
Theorem 3 (Holley and Strook 88). Under assumptions of 2.2, there exist two positive constants 0 <
c ≤ C < +∞ such that ∀ β ∈ R+ ,
?

?

ce−βm ≤ γ(β) ≤ Ce−βm
R
R
where γ(β) = inf{− φLβ φ dµβ : kφkµβ = 1 and φdµβ = 0} and m? is the maximum depth of a well
containing a local minimum defined in Equation(4).

Remark The constant m? is always strictly positive as soon as the function has a strict local minimum
that is not global. This is generally the case in our setting. Also, we alwaysRhave m? ≤ M .
Following of the Proof. Using the definition of ft , one can see that ft dµβt = 1, hence applying
ft − 1
the theorem for φ =
gives:
kft − 1kµβt
−

X
x

?

φ(Lβt φ)(x)µβt (x) ≥ ce−βt m .

This and the definition of Lβt imply:
X
?
ft (x)(Lβt ft )(x)µβt (x) ≤ −ce−βt m kft − 1k2µβ .
t

x

J is a positive function bounded by M on E 0 . For all x ∈ E \ E 0 , the only points where J > M , we have
that J(x) = +∞ and therefore µβt (x) = 0. This implies that for all measurable functions g,
X
(J(x) − hJiµβt )g(x)µβt ≤ M kgkµβt .
x∈E

Putting all these terms together gives:


d
M 0
?
+
−
kft − 1k2µβ ≤ 2 −ce−βt m (1 + +
)
+
(
−

)
+
β
kft − 1k2µβ
t
t
t
t
t
dt
2 t
0
+ 2M βt kft − 1kµβt

(17)

−
+ 2(+
t − t ).

We denote ut = kft − 1k2µβ . Considering the fact that +
t is a positive function we have:
t

u0t



−βt m?

≤ 2 −ce

+

(+
t

√
+ 2M βt0 ut

−

−
t )


M 0
β ut
+
2 t

(18)

−
+ 2(+
t − t ).

Using that, ∀x ∈ R,

1 2
4x

+ 1 ≥ x, we get:


M
M 0
+
−
0
−βt m?
ut ≤ 2 −ce
+ (t − t ) + (
+
)β ut
2
4 t
−
+ 2M βt0 + 2(+
t − t ).

?

−
0
Let At = 2ce−βt m and Bt = 2(+
t − t ) + 2M βt .

22

(19)

Applying Grönwall’s Lemma for the previous relation gives:
Z t
Rt
Rt
−As +Bs ds
0
Bs e s −Ah +Bh dh ds.
ut ≤ u0 e
+

(20)

0

Under Assumptions 2.2, there exist b, d > 0 such that βt = b log(1 + td). This implies:

m? b
bd
1
0
−m? βt
and e
=
.
βt =
1 + td
1 + td
−
α
Using the definition of +
t , t and the fact that nt = (1 + t) one can check that:




1
1 log t
At = O m? b and Bt = O
∨
t tα/2
t

We can see that if Bt = o(At ) the second term of Equation (20) is bounded and gives us a finite upper
bound on ut . This happens as soon as:
m? b < 1 ∧ α/2
(21)

However, the condition given by Equation (21) is sufficient yet not necessary. For α > 2, Bt becomes
of the order O(1/t) and thus we can choose d in a way that preserves a finite upper bound of Equation
(20) even for m? b = 1. One can check by direct computation that this is true for any d < cm? /M .
Let βt and nt be chosen in order to comply to one of the two previously mentioned conditions. Then
there exists a constant K 0 such that
ut ≤ K 0 for all t ∈ R+
(22)
To complete the proof of Theorem 2 one can observe that for all t ∈ R+ , and all  > 0 :
et ∈ c χ ) = E(1[J ? +,+∞) (J(X
et ))).
P(X

Using the Cauchy-Schwarz inequality and the upper bound given by Equation (22) we obtain:
Z
e
E(1[J ? +,+∞ (J(Xt ))) = 1[J ? +,+∞ (J(x))ft dµβt (x)
R

≤

 1 Z
1
2
2
2
(ft ) dµβt (x)
1[J ? +,+∞) (J(x))dµβt (x)

Z

2

R

(23)

R

≤kft kL2µ (µβt (c χ ))1/2
βt

with K =

5.2

√

≤K(µβt (c χ ))1/2
K 0 + 1. This completes the proof of Theorem 2.



Convergence rate

A first rate of convergence can be deduced from Theorem 2 using the concentration speed of the Gibbs
measure on χ .
P

x∈c χ

c

µβt ( χ ) = P

x∈E

e−βt J(x)
e−βt J(x)
P

−βt J(x)
x∈c χ e
P
−βt J(x) +
−βt J(x)
x∈c χ e
x∈χ e
?
(|E| − |χ |)e−βt (J +)
0 + |χ |e−βt J ?

=P
≤


≤


|E|
− 1 (1 + td)−b
|χ |
23

(24)

As the dependency of K (Theorem 2) in b and α is not explicit, we can however not deduce an optimal
choice of (b, α) from this bound. This can be achieved if we assume that Equation (21) holds and
distinguish the two cases α ≤ 2 and α > 2. Indeed, we can then improve the bound on ut and derive
a more accurate convergence rate of the algorithm. This rate can then be optimized to obtain either an
upper bound of the probability of convergence to χ for a fixed computational budget or the minimal
computational budget at a fixed risk of convergence out of χ .
Theorem 4. Under assumptions of Section 2.2, suppose: βt = b log(td + 1), nt = (1 + td)α and m? b <
min(α/2, 1):
• if α ≥ 2, let b be such that m? b < 1 and let γ ∈ (0, α/2 − m? b),
Then, there exist Γγ , Γ2 > 0 such that for t large enough, for all  > 0,
? b−1−b)/2

et ∈ c χ ) ≤ Γγ Γ2 (1 + td)(m
P(X

+ Γ2 (1 + td)−b

• if α < 2, let b be such that m? b < α/2 and let γ ∈ (0, α/2 − m? b),
Then, there exist Γγ , Γ2 > 0 such that for t large enough, for all  > 0,
et ∈ c χ ) ≤ Γγ Γ2 (1 + td)(m
P(X

? b−α/2+γ−b)/2

+ Γ2 (1 + td)−b


Remark 1. γ is not a new parameter of the NSA algorithm. This is a technical element that enables the
tuning of the computational complexity bounds of Section 5.3. As shown in the Appendix A, Γγ is of the
order of 1/γ.
Remark 2. This two bounds display the trade off between the convergence rate of the Gibbs measure to
the uniform distribution over the global minima and the rate of convergence of the NSA process to the
Gibbs measure. For the first bound, considering α > 2 we recover the classical rate of convergence of
the simulated annealing in the noise free case. This corresponds to the result of [15]. The second bound
provides the rate of convergence for a choice of α < 2. It can be seen that b will have to be reduced to
ensure the convergence and thus this bound exhibits clearly the trade off between cooling and estimation.
Proof. Under assumptions of Theorem 4, the following bound on ut can be derived from Grönwall’s
Lemma (for details see Appendix A):
(
?
Γγ (1 + td)m b−1 if α ≥ 2
ut ≤
(25)
?
Γγ (1 + td)m b−α/2+γ if α < 2
et does not belong to the optimal set χ (cf.
Thus we can compute a new bound on the probability that X
6):
Z
et ∈ c χ ) = 1c χ (J(x))ft dµβ (x)
P(X
t
RZ
Z
c
=
1 χ (J(x))(ft − 1)dµβt (x) +
1c χ (J(x))dµβt (x)
R

≤

R+

Z
R+

2

(ft − 1) dµβt (x)

Z

q
≤ ut µβt (c χ ) + µβt (c χ )

1/2

12c χ (J(x))dµβt (x)

+ µβt (c χ )

R+

(26)

This means that if there exist (α, b) such that ut = O(µβt ) the convergence rate in the noisy case will be
of the same order as in the classical one, but for a smaller b.
24

Using the previous inequality, Equation (25) and the concentration rate of the Gibbs measure given
by Equation (24) we have:
c

e t ∈ χ ) ≤
P(X
where Γ2 =

5.3

|E|
|χ |

(
m? b−1−b
2
Γγ Γ2 (1 + td)
+ Γ2 (1 + td)−b
Γγ Γ2 (1 + td)

m? b−α/2+γ−b
2

+ Γ2 (1 +

if α ≥ 2

td)−b

if α < 2

− 1.



Computational complexity of NSA

Given the convergence rate of the algorithm, we can define T ? such that the confidence inequality constraint is satisfied at time T ? .
T be the number of cost function evaluations made by the NSA until time T . This is a random
Let Ncall
variable. We define the computational cost of the algorithm as the expectation of this random variable.
It can be written as:


X

T
= E
1Tk <T NTk  .
E Ncall
k≥1

Lemma 3. Let δ,  > 0, γ ∈ (0, α/2 − m? b) and
1
T =
d
?


max

2Γγ
δ

!
2/(min(1, α −γ)−m? b+b) 
 !
2
2Γ2 1/b
,
−1 .
δ

et ∈ c χ ) ≤ δ and the computational cost up to time T ? is bounded:
Then, for all t ≥ T ? , P(X

T∗

E Ncall

1
≤ max
d



2Γγ
δ

!
2(α+1)/(min(1, α −γ)−m? b+b) 

2
2Γ2 (α+1)/b
,
.
δ


Proof. In order to prove this statement, we use the inequalities from Theorem 4 treating each term
separately.
α
?
We consider T1 , T2 such that Γγ Γ2 (1 + dT1 )m b−min(1, 2 −γ)−b = δ/2 and Γ2 (1 + dT2 )−b = δ/2. This
implies:

1 + dT1 =

2Γγ
δ

2/(min(1, α −γ)−m? b+b)



2

and 1 + dT2 =

2Γ2
δ

1/b
.

(27)

Now we can define T ? , the time after which the current state of the NSA belongs to χ with probability
et ∈ χ ) ≥ 1 − δ:
at least 1 − δ, i.e., ∀t > T ? , P(X
T ? = max(T1 , T2 ).
?
We are interested in the computational cost up to time
PT , more precisely
 the expected number of
?
Monte Carlo simulations used up to T . This is given by E
k≥1 1Tk <T ? NTk . The value of this quantity

25

cannot be computed exactly, but it can easily be upper bounded.


 


X
X
E
1Tk <T ? NTk  =E E 
1Tk <T ? NTk  |(Tk )k=1···+∞ 
k≤1

k≥1


=E 


X

1Tk <T ? nTk 

k≥1


≤E 


X

1Tk <T ?  nT ?

k≥1

The last inequality is implied by the fact that nt is an increasing function. Since
variable of parameter T ? , using the definition of nt one can see that:


X
E
1Tk <T ? NTk  ≤T ? (1 + dT ? )α

P

k≥1 1Tk <T ?

is a Poisson

(28)

k≤1

1
≤ (1 + dT ? )α+1 .
d
We conclude using Equation (27).



The rate of growth of the total computation number is mainly driven by the exponent of 1δ in the cost
function. We are looking for the couple (α, b) that minimizes this quantity and fulfills the requirements
of Theorem 4. We can split the problem into two sub-problems:
Case 1:

α
2

−γ >1

min max
b,α

s.t.

α+1
2(α + 1)
,
?
1 − m b + b b


(29)

m? b < 1 and α − 2γ > 2
Case 2:

α
2

≤1

min max
b,α

s.t.
0<γ<

2(α + 1)
α+1
,
α/2 − γ − m? b + b b


(30)

α
− m? b and α − 2γ ≤ 2
2

The solution of Equation (29) is obvious, the minimal value for α and the maximal for b, i.e, α must
be as close to 2 as possible and b = m?1+ . As for Equation (30), we consider two sub-cases. First suppose
that:
2(α + 1)
α+1
≥
⇐⇒ α/2 − γ − m? b ≤ b.
(31)
?
α/2 − γ − m b + b
b

The function we want to minimize is strictly decreasing in α and strictly increasing in b, so its minimum
value is attained for the maximal value of α and the minimal value of b, under the domain constraints
given by Equation (30) and Equation (31), so the solution is:
α = 2(1 + γ) and b >
26

1
.
m? + 

(32)

In the second sub-case, supposing that the inequality Equation (31) is inverted, the problem can be
resumed at minimizing (α + 1)/b, a decreasing function with respect to b, for
b≤

α/2 − γ
and α ≤ 2.
(m? + )

Replacing b by its maximal value the objective function becomes a decreasing function in α, and therefore
we obtain the same solution as before, defined in Equation (32). This is a quite comprehensive result, as
it indicates that the lower the required accuracy in the solution space is, i.e.,  increases and thus the size
of χ does too, the faster the temperature can decrease to zero. We need to explore less the state space.
Corollary 1. For the optimal parameters choice defined in Equation (32), an -optimal solution is returned by NSA with probability 1 − δ at a computational cost at most :
1
d



2Γγ
δ

 m? + (3+2γ)


,

where Γγ is defined in Theorem 4.
This is rather costly but represents a general bound with few constraints on the function J. However,
if the function J has additional properties the bound can be significantly improved:
Corollary 2. Suppose that J has no well containing a local minimum, apart from the one containing the
global minimum, i.e. m? = 0, then an -optimal solution is returned by NSA with probability 1 − δ at a
computational cost at most :
!3
2 log 1δ
.
d
Remark 3. We recover the polynomial dependency in 1/ and log 1/δ of the state-of-the-art complexity
results (c.f. [25] and [24]) which are of the order of −2 log( 1δ ) for strongly convex cost functions. As
we relax this assumption and only consider cost functions with no local minimum, it seems coherent to
observe a slight degradation of the complexity.
Proof. In order to have an estimate of the computational cost in this setting we follow the same
method as before and highlight only the main steps of the proof. First remark that in this case, Theorem
3 states that there exist C, c > 0 such that ∀β ∈ R+ :
c ≤ γ(β) ≤ C

(33)

This changes the differential inequality obtained for ut = kft − 1k2µβ and thus Equation (19) becomes :
t



−
u0t ≤ 2 −c + (+
t − t ) + (



M
M 0
−
+
)β ut + 2M βt0 + 2(+
t − t ).
2
4 t

We can apply Grönwall’s Lemma and obtain the same type of inequality as before:
Z t
Rt
Rt
−A
+B
ds
s
s
+
ut ≤ u0 e 0
Bs e s −Ah +Bh dh ds.
0

−
0
where Bt has the same form as before, Bt = 2(+
t − t ) + 2M βt and At = 2c .
The convergence of ut towards 0 can be proved now for a larger class of functions nt , βt , since:

√
At = O(1) and Bt = O(βt / nt ∨ βt0 ).
27

(34)

We no longer need to impose βt = O(logt ). Let α, b, d > 0. Define
nt = (1 + t)α and βt = d(1 + t)b .
Using Equation (34) one can check that we have a finite upper bound on ut as soon as:
{b < α/2 ∧ 1, d > 0} or {b = 1, α ≥ 2, 0 < d < c}.
This in particular implies that the NSA algorithm converges a.s. to the set of global minimums of J.
Furthermore for the first set of conditions one can prove using the same technique as in Appendix A that
:
ut = O(tmax(b−1,b−α/2) ).
This means that there exits Γ0γ > 0 such that for t large enough ut ≤ Γ0γ t−γ , where γ = − max(b − 1, b −
α/2). Using this, Equation (26) and Equation (24), for t large enough, we get:
et ∈ c χ ) ≤ Γ0 t−γ/2 e−d(1+t)
P(X
γ
b /2

≤ e−d(1+t)

b /2

b

+ Γ2 e−d(1+t)

b /2

(Γ0γ t−γ/2 + Γ2 e−d(1+t)

)

−d(1+t)b /2

≤e

The last inequality is valid as soon as t > max



2 log(2Γ2 )
d

1
b

− 1, 2Γ0γ

2
γ


. This is not a restrictive

condition. Take for example the minimization of the k.k1 over the subset subset E = {x ∈ Zp , kxk∞ ≤ n}
for some n ∈ N. As Γ2 = kEk − 1 = (n + 1)p − 1, the time for which the first part of the condition is
fulfilled only grows linearly with the dimension of the search space. We show latter on that the optimal
choice for b is one and thus the second part of the condition can be omitted.
Let δ > 0 be a fixed. Using the previous inequality one can compute T ? such that the confidence
inequality constraint is satisfied:


2 log δ 1/b
?
T,δ = −
− 1.
d
Regarding the computational cost we remind the reader that Equation (28) implies:
T?
E(Ncall
)

?

≤ nT ? T ≤



2 log δ
−
d

 1+α
b

.

We can optimize this bound with respect to α and b in the same way as for Corollary 1. This leads to
α = 2 and b < 1 and thus to the desired results:


2 log δ 3
T?
E(Ncall ) ≤ −
.
d


6

Numerical experiments

In this section we first present some test cases, for which we use an additive Gaussian noise at each
evaluation. We recover the theoretical results introduced by [15]. In a second part we present some results
for the aircraft trajectory optimization problem. In this case the solution of the problem is unknown. We
can only observe the total cost improvement in comparison with a trajectory optimized for a similar but
deterministic setting.
28

U(x)

2

12

5

10

7

5

2

7

8

d*

3

1

13

11

10

4

0

6

15

12

9

14

S*

Figure 2: B. Hajek test case for the simulated annealing in a deterministic environment
Basic exemple
The first experimental setting we consider, was introduced in [16]. The cost function and the neighbourhood structure are represented on Figure 2. This is of particular interest as the function has two basins
from which it is hard to escape. B. Hajek has shown that the following holds:
Theorem 5 ([16]).
If βk = b log(k + 2), then

b ≤ d? ⇔ lim P(Xk ∈ S ? ) = 1, where d? is the maximum depth of a cup
k→∞

containing a local but not global minima. The depth of a cup is the maximal energy difference between
two of it states and (Xk )k∈N denotes the Markov chain generated by the classical simulated annealing. 
For a complete definition of d? , see [16].
We add Gaussian noises to the cost function of Figure 2 with different variance levels to highlight the
fact that if no sampling is performed the simulated annealing performance becomes rapidly very poor
as the variance increases. On the other hand it appears that the performance of the NSA for a linear
increase of the mean number of samples is as good a quadratic one. These results are summarized on
Figure 3.
Ackley test function
We introduce a second test case to further asses these observations. We consider the uniformly (2000
points) discretized version of the Ackley function in one dimension on [−100, 100]. This function has
many local minima as shown on Figure 4. Figures 5 displays the convergence results for different levels
of variance of the noise for each estimation schedule introduced in this paper. We observe that the only
case where the convergence is not impacted by the noise variance increase is the nt = t2 case.
These results highlight the fact that a logarithmic sampling schedule is not appropriate in general,
even in the Gaussian case. This invalidates partially the hypotheses introduced in [13]. A clear gap is
highlighted between the linear and the quadratic schedule.

Aircraft trajectory optimization
We use a black box trajectory evaluator for a long range commercial aircraft. We consider a direct
shooting method for optimizing the vertical part of the trajectory. As displayed on Figure 6, the vertical
path is made of a sequence of flight segments at constant altitude called steps. The transitions between
29

Figure 3: Convergence performance of NSA for the Hajek setting

Figure 4: Ackley 1D Test Function

30

Figure 5: Performance of the NSA algorithm as a function of the level of noise on the evaluation of the
cost function
those steps are called step climbs. This has been put in place by the international authorities to ease
the air traffic control. Aircraft can only fly at a finite set of altitudes. The steps climbs are transition
phases that must be very short. The Figure 6 is a conceptual. It does not reflect the real scale of the
different phases. Our optimization variables are the vectors of position of the steps and the vector of
steps’ altitude, denoted respectively x and h on Figure 6. The structure of this airspace strongly limits
the number of steps. We will only consider the problem with an a priori number of steps. There are two
main reasons why the aircraft might vary its altitude during a flight (optimizing fuel consumption and air
traffic control). Because of fuel consumption, the aircraft weight is decreasing during the flight. Analyzing
the laws of flight physics, it can be shown that there exists an altitude at which the fuel consumption per
flown distance unit is minimal. It can also be shown that this altitude increases as the weight decreases.
This last statement is however only true if there is no wind. It is easily understandable that for some

Figure 6: Aircraft trajectory, structure of the vertical path
31

Cost

Number of iterations

Figure 7: NSA descent: Aircraft trajectory optimization problem

Figure 8: Sampling of the cost function along the first step position for a 3 step vertical path

particular wind map configuration it might be preferable to target lower altitudes at lower weights.
The choice of the vertical path must be declared to the authorities before the flight to ensure traffic
manageability. Airlines operating aircraft have therefore a stochastic optimization problem to solve. This
is a stochastic problem for two main reasons. First, they only access predicted weather conditions that
suffer some uncertainty. Second, the airspace is not empty and sometimes air traffic controllers might
refuse some altitude changes because of the presence of other aircraft. As the weather, the traffic is not
known in advance.
We applied NSA to the problem of finding an optimal 3 steps configuration. An example of the current
solution cost evolution with respect to the number of iterations is displayed on Figure 7.
We observe a very quick convergence to a low cost trajectory. We do not claim it is a general behaviour.
It might be due to the structure of the cost function. Figure 8, shows how the cost evolves with respect to
the ground position of the first step. It is obviously not convex but has some regularity. We can observe
some flat parts. This explains why gradient based methods would fail solving this problem.
As for the previous problems we have observed that the increase sampling condition must be satisfied
to ensure a good behaviour of the algorithm.
32

A

Proof of bound Equation (25)

Proof.[Proof of Theorem 4]
Let βt = b log(td + 1) and nt = (1 + td)α .
Recall Equation (18):
 


M
M
0
−βt m?
+
−
βt0 ut + 2M βt0 + 2(+ − − )
ut ≤ 2 −ce
+ ( −  ) +
+
2
4
?

Let At = 2ce−βt m and Bt = 2(+ − − ) + 2M βt0 .
Applying Grönwall’s Lemma for the previous relation gives:
Z t
Rt
Rt
−A
+B
ds
s
s
Bs e s −Ah +Bh dh ds
ut ≤ u0 e 0
+
0

Under Assumptions 2.2, there exist b, d > 0 such that βt = b log(1 + td). This implies:
βt0

bd
?
=
and e−m βt =
1 + td



1
1 + td

m? b
.

Using the definition of − , + we have:


βt2 σ 2
2(1−δ)n
√
t
2 2
 e
+ − −
+ eβt σ −ant  .
t = 2 2πβt σ p
(1 − δ)nt
This implies that when t goes to infinity:




1
1 log t
At = O m? b and Bt = O
∨ α/2
t t
t
In order to highlight the mains ideas of the proof we will try to simplify
the

 notations as much as possible.
log t
1
First observe that for all α > 0 and all 0 < γ < α/2 , α/2 = o α/2−γ .
t
t
Hence we can assume there exist A, B > 0 and δ1 , δ2 such that:
At = A

d
d
and Bt ≤ B
δ
1
(1 + td)
(1 + td)δ2

where δ1 = m? b and δ2 = min(1, −γ + α/2). Since min(1, α/2) > m? b, and γ can be chosen arbitrarily
close to 0, we choose it such that δ1 < δ2 . This means that 0 < γ < α/2 − m? b.

Remark 4. The choice of γ influences the choice of B. If Bt = O
that Bt ≤ CB


log t
, there exists CB > 0 such
tα/2

log t
log t
, ∀t. The constant B is then such that ∀t, CB γ ≤ B. Hence we can choose:
α/2
t
t
B=

Let Tt1 = u0 e

Rt
0

−As +Bs ds

and Tt2 =

Rt

0 Bs e

Rt
s

CB
.
eγ

−Ah +Bh dh

33

(35)
ds.

The first term Tt1 is always easy to deal with and one can check that under the theorem’s assumptions
we always have Tt1 = o(1/tδ1 −δ2 ) when t goes to infinity. As for the second term, using a substitution
gives:
Z 1+td
B Rs1+td − Aδ + Bδ dh
2
h 1
h 2
ds
Tt ≤
e
s δ2
1
Z 1+td
1−δ2
1−δ1
A(1+td)1−δ1
B(1+td)1−δ2
B As1−δ
−
+
− Bs
1−δ
1−δ
1−δ2 ds
1
2
1
≤e
(36)
e
s δ2
1
For the last inequality we assume δ2 6= 1 which corresponds to the case α ≤ 2.
R 1+td B As1−δ1 − Bs1−δ2
1−δ1
Bs1−δ2
1−δ2 ds and f = As
Let It = 1
e 1−δ1
s
1−δ1 − 1−δ2 . Let T0 be such that for all s ≥ T0 ,
sδ 2
B
1/(δ2 −δ1 ) ). We can write I as follows:
sδ2 −δ1 ≥ B+1
t
A (for instance T0 = ( A + 1)
Z

T0

It =
1

= KT0
= KT0

Z 1+td
B fs
B fs
e
ds
+
e ds
δ2
sδ2
s
T0
Z 1+td
B
+
efs fs0 ds
δ2 (As−δ1 − Bs−δ2 )
s
T0
1+td Z 1+td

AB(δ2 − δ1 )sδ2 −δ1 −1 fs
B
fs
+
e
e ds
+
Asδ2 −δ1 − B
(Asδ2 −δ1 − B)2
T0
T0

A(δ2 − δ1 )s−δ1 −1
goes to 0 when s goes to infinity. Moreover one can check that for all
(Asδ2 −δ1 − B)2
s ≥ T0 this quantity is smaller than 1/2. Using this we get:
Since δ1 < δ2 ,


It ≤ KT0 +
and therefore for all t ≥ T0 :

This gives:
≤ 2e

−f1+td



fs

Asδ2 −δ1


It ≤ 2

Tt2

B
−B

B
Asδ2 −δ1 − B
B

e

fs

fs

1+td

e

T0

1
+ It
2

1+td
+ 2KT0
1

1+td

!

+ KT0
e
Asδ2 −δ1 − B
1




B
B
f1
≤2
+ 2 KT0 −
e
e−f1+td
A−B
A(1 + td)δ2 −δ1 − B

Regrouping the terms we obtain for all t ≥ T0 :



B
2B
f1
ut ≤ u0 + 2 KT0 −
e
e−f1+td +
A−B
A(1 + td)δ2 −δ1 − B
Since the first term is a O(e−ft ) and therefore a o(tδ1 −δ2 ) it is obvious that:


1
ut = O δ2 −δ1 when t → ∞
t
This means that for all γ ∈ (0, α/2 − m? b) there exists Γγ > 0 such that:
34

(37)

(38)

ut ≤ Γγ t−α/2+m

? b+γ

for all t ≥ T0 .

Since ft = O(t1−δ1 ), e−ft tδ1 −δ2 goes very fast to zero and therefore the size of Γγ is mainly driven by

Using Equation (35) one can see that:

Γγ '

(39)
2B
.
A

1
.
γ

If α > 2, δ2 = 1 and Equation (36) becomes of the form:
Z 1+td
1−δ1
A(1+td)1−δ1
B As1−δ
−
+B log(1+td)
2
1−δ
1
1 ds.
Tt ≤ e
e
s1+B
1
Remark 5. In this case one can choose γ = α/2 − 1, this way min(1, α2 − γ) = 1 and Γγ is minimal.


Definition of m?

B

In this section we prove that the definition of m? , i.e. Equation (4), is equivalent to the definition provided
by [17].




?
Lemma 4. Let, mHS := max min max J(z) − J(y) − J(x) + minu J(u) .
x,y∈E

p∈Pxy

z∈p

Then
m? = m?HS

Proof. Let x, y ∈ E and denote Hxy := min

p∈Pxy


max J(z) .



z∈p

First it can be noticed that if x is a global minimum of J then we have
Hx,y − J(y) − J(x) + min J(u) = Hx,y − J(y)
u

(40)

Thus m?HS ≥ Hx? ,y − J(y) for any y in E, where x? is a global minimum of J.
Recall that m? = max {Hxy − max (J(y), J(x))}. As the set of paths going from x to y containing a
x,y∈E

global minimum x? is a subset of the paths going from x to y, we have:
Hxy ≤ max (Hx? x , Hx? y )
Let x, y ∈ E such that m? = Hxy − max (J(y), J(x)),

m? ≤ max (Hx? x , Hx? y ) − max (J(y), J(x))
≤ max (Hx? x − J(x), Hx? y − J(y))
≤ m?HS

On the other hand, as ∀ x, y ∈ E, we have − min (J(y), J(x)) + minu J(u) ≤ 0, so
Hxy − J(y) − J(x) + min J(u) ≤ Hxy − max (J(y), J(x))
u

This implies m?HS ≤ m? , which completes the proof.
35



References
[1] E. Aarts and J. Korst. Simulated annealing and Boltzmann machines. New York, NY; John Wiley
and Sons Inc., 1988. 3, 4, 8, 9
[2] Emile HL Aarts and Jan HM Korst. Boltzmann machines for travelling salesman problems. European
Journal of Operational Research, 39(1):79–95, 1989. 8
[3] Jean-Yves Audibert and Sébastien Bubeck. Best arm identification in multi-armed bandits. In
COLT-23th Conference on Learning Theory-2010, pages 13–p, 2010. 2
[4] Peter Auer, Nicolò Cesa-Bianchi, and Paul Fischer. Finite-time analysis of the multiarmed bandit
problem. Machine learning, 47(2-3):235–256, 2002. 2
[5] Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and geometry of Markov diffusion
operators, volume 348. Springer Science & Business Media, 2013. 20
[6] John T. Betts. Survey of numerical methods for trajectory optimization. Journal of Guidance,Control,
and Dynamics, 21(2):193–207, March 1998. 6
Concentration inequalities:
[7] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.
nonasymptotic theory of independence. Oxford university press, 2013. 17

A

[8] J. Branke, S. Meisel, and C. Schmidt. Simulated annealing in the presence of noise. Journal of
Heuristics, 14(6):627–654, 2008. 4
[9] Sébastien Bubeck, Rémi Munos, Gilles Stoltz, and Csaba Szepesvári. X-armed bandits. Journal of
Machine Learning Research, 12:1587–1627, 2011. 2
[10] Adam D Bull. Convergence rates of efficient global optimization algorithms. The Journal of Machine
Learning Research, 12:2879–2904, 2011. 2
[11] Adam D Bull et al. Adaptive-treed bandits. Bernoulli, 21(4):2289–2307, 2015. 2
[12] Laurence Charles Ward Dixon and Giorgio Philip Szegö. Towards global optimisation 2. NorthHolland Amsterdam, 1978. 2
[13] T.MA. Fink. Inverse protein folding, hierarchical optimisation and tie knots. PhD thesis, University
of Cambridge, 1998. 4, 5, 8, 29
[14] S. B. Gelfand and S. K. Mitter. Simulated annealing with noisy or imprecise energy measurements.
J. Optim. Theory Appl., 62(1):49–62, 1989. 3
[15] W.J. Gutjahr and G.Ch. Pflug. Simulated annealing for noisy cost functions. Journal of Global
Optimization, 8(1):1–13, 1996. 3, 4, 8, 24, 28
[16] B. Hajek. Cooling schedules for optimal annealing. Mathematics of operations research, 13(2):311–
329, 1988. 3, 9, 10, 12, 29
[17] R. Holley and D. Stroock. Simulated annealing via Sobolev inequalities. Comm. Math. Phys.,
115(4):553–569, 1988. 3, 6, 8, 9, 10, 11, 20, 22, 35
[18] Tito Homem-de Mello. Variable-sample methods and simulated annealing for discrete stochastic
optimization. 2000. 3, 8
[19] Reiner Horst and Panos M Pardalos. Handbook of global optimization, volume 2. Springer Science
& Business Media, 2013. 2
36

[20] Donald R Jones, Cary D Perttunen, and Bruce E Stuckman. Lipschitzian optimization without the
lipschitz constant. Journal of Optimization Theory and Applications, 79(1):157–181, 1993. 2, 6
[21] Donald R Jones, Matthias Schonlau, and William J Welch. Efficient global optimization of expensive
black-box functions. Journal of Global optimization, 13(4):455–492, 1998. 2, 6
[22] John DC Little, Katta G Murty, Dura W Sweeney, and Caroline Karel. An algorithm for the traveling
salesman problem. Operations research, 11(6):972–989, 1963. 2
[23] Rémi Munos. From bandits to monte-carlo tree search: The optimistic principle applied to optimization and planning. 2014. 2
[24] Arkadi Nemirovski, D-B Yudin, and E-R Dawson. Problem complexity and method efficiency in
optimization. 1982. 27
[25] Michael Woodroofe. Normal approximation and large deviations for the robbins-monro process.
Probability Theory and Related Fields, 21(4):329–338, 1972. 27

37

