Laetitia Andrieu

Guy Cohen
email: guy.cohen@mail.enpc.fr

Felisa J Vázquez-Abad

Stochastic Programming with Probability Constraints

Keywords: Probability constraints, stochastic programming, stochastic gradient algorithm, stochastic approximation

In this work we study optimization problems subject to a failure constraint. This constraint is expressed in terms of a condition that causes failure, representing a physical or technical breakdown. We formulate the problem in terms of a probability constraint, where the level of "confidence" is a modelling parameter and has the interpretation that the probability of failure should not exceed that level. Application of the stochastic Arrow-Hurwicz algorithm poses two difficulties: one is structural and arises from the lack of convexity of the probability constraint, and the other is the estimation of the gradient of the probability constraint. We develop two gradient estimators with decreasing bias via a convolution method and a finite difference technique, respectively, and we provide a full analysis of convergence of the algorithms. Convergence results are used to tune the parameters of the numerical algorithms in order to achieve best convergence rates, and numerical results are included via an example of application in finance.

Introduction 1.Constrained Optimization in a Stochastic Setting

Optimization Theory provides a convenient approach to formulate and solve problems involving conflicting objectives, which is generally the challenge present in decision making situations. The main idea is to aggregate as many objectives as possible into a single objective function, which may be straightforward when those objectives are amenable to an expression into a common unit, say, a currency unit as dollar or euro. In this objective aggregation, weights are allocated to each term in order to reflect preferences or priorities. However, there might be other objectives that can hardly be expressed in a unit commensurable with the previous ones (examples to come hereafter). In such a case, it is better to introduce those other objectives through constraints, that is, each such objective should not exceed a prescribed level. The constraint levels are set a priori, as are the weights for the different terms in the cost function.

Duality Theory provides the tools to evaluate the sensitivity of the optimal solution (cost) to those prescribed constraint levels. In mathematical terms, let u be the decision variable in a Hilbert space U, J : U → R the cost function, and Θ : U → R d the constraint function. We consider problems of the type:

min u∈U ad J(u) s.t. Θ(u) ≤ α , (1)
where U ad is an "admissible" or "feasible" closed convex subset of U and inequalities in the constraints involving Θ are understood componentwise. Introduce the multiplier λ (in R d + ) and the Lagrangian L(u, λ) = J(u) + λ, Θ(u) -α ,

where •, • denotes the scalar product. Kuhn-Tucker optimality conditions characterize an optimal multiplier λ ♯ which can be interpreted as the sensitivity of the optimal cost function J(u ♯ ) (where u ♯ denotes the solution of problem [START_REF] Arrow | Studies in nonlinear programming[END_REF]) with respect to α (up to a change of sign). When random factors affect the outcomes of a decision, a classical approach is to assume that the probability distribution of those factors is known and to appeal to stochastic optimization. Call ξ the corresponding random variable, then the objective function is usually expressed in terms of an expectation of some cost function of the form J(u) = E j(u, ξ) .

In the stochastic situation, modelling choices for aggregation of objectives, weights and constraints are similar to the deterministic case. However a new question also arises regarding constraints, namely, constraints can be formulated in various ways: "almost surely", "in expectation", "in probability", etc.

The first possibility ("almost sure" constraints) means that certain quantities θ(u, ξ) depending on decision variables and affected by random factors should satisfy equality or inequality for "almost all" values of those random factors (according to their probability distributions). This is in particular the case of constraints which express "laws of nature" which are part of the mathematical model of the system under consideration. However, regarding objectives or "wishes", such strict constraints are generally inappropriate from the economic or simply realistic point of view. Suppose for example that a pressure should not exceed a certain level beyond which death will almost surely happen. First of all, observe that it is hard to aggregate such an objective (actually, that to stay alive) with other more economic objectives which aim at saving money. Second, under the constraint that the pressure "almost never" exceeds the dangerous level, the operation can be extremely costly if not simply impossible. That is, some risk must be accepted for the operation to be economically viable.

The second possibility (constraints "in expectation") means that, given a decision, the expected value Θ(u) = E θ(u, ξ) of a critical quantity (a pressure in our example) should not exceed a certain level. Such a formulation is generally mathematically attractive, but it is difficult to understand how much risk is involved in choosing such or such prescribed level. Indeed, given a decision u, the pressure (to keep on with our example) becomes a random variable θ(u, ξ) with a certain distribution (which is affected by the chosen decision), and the only thing one ask is that the first moment (the expectation) of this random variable stay below a prescribed level, but with no direct control on how much of the probability mass will lie beyond that prescribed level.

The third possibility advocated to (constraints "in probability" or "probabilistic constraints") means that one accepts that the critical quantity (the pressure, say) remains under the prescribed level not "almost always" as earlier, but with a certain probability whose value must be chosen. In mathematical terms, one now considers the problem min

u∈U ad E j(u, ξ) s.t. P θ(u, ξ) ≤ α ≥ π . (3)
This chosen probability value π exactly reflects the risk one is ready to assume (in contrast with the previous approach of constraints in expectation). As discussed earlier, duality should then help in evaluating the sensitivity of the optimal cost function with respect to this accepted, but arbitrarily fixed, level of risk.

Quantitative vs. Qualitative Risk Measures

We now motivate the interest of probability constraints in contrast with other measures of risk. Before choosing a risk measure, it is very important to know which type of failure we are interested in: quantitative failure or qualitative failure. For example, a power supply company would minimize its cost under the constraint of supplying the demand. If that demand cannot be fully supplied, it matters to know which percentage of it will not be covered and during which amount of time. This is what we mean by "quantitative failure": introducing a penalty for the total amount of demand not supplied directly into the cost function, or choosing to constrain a quantity which accounts for the amount of supply failure is appropriate in that situation. On the contrary, when simply going beyond a certain threshold causes death, it not does matter to know by which amount that threshold has been exceeded -this is what we mean by "qualitative failure" -but it does matter to know the likelihood of going beyond that critical threshold. Probability constraints are particularly adapted to this latter situation.

In fact, because of the mathematical difficulties raised by probability constraints, these constraints must be exclusively used in the case of qualitative failure problems. For quantitative failure problems, there are other risk measures with better mathematical properties (e.g. convexity), like Conditional Value-at-Risk (CVaR) for instance. Introduced in [START_REF] Rockafellar | Conditional Value-at-Risk for General Loss Distributions[END_REF], CVaR is one of the most popular risk measure in finance. CVaR is the average of a random variable for the worst scenarios. Denote by α u (π) the quantile function of the distribution of θ(u, •) with confidence level π (also called Value-at-Risk). Then, CVaR, denoted by φ π (u), is defined by

φ π (u) = E θ(u, ξ) | θ(u, ξ) ≥ α u (π) .
The risk constraint will be then φ π (u) ≤ α, where α represents the accepted level of risk, and π is fixed a priori.

Notice that the critical threshold α in the probability constraint is generally provided by technical considerations, whereas π characterizes the level of risk one is ready to accept. That is, the decision maker may bargain about the constraint level π but not on that threshold α which is a technical data. With the CVaR approach, this α disappears from the formulation and we believe that this is a weakness of this approach. Moreover, in the case of "qualitative failure", there is no meaning in averaging values of θ beyond a threshold which is supposed to be fatal.

About this Paper

Problem (3) is the class of problems considered in this paper. Its advantage is again the fact that the meaning of constraints in terms of risk assumed is of immediate perception. Its drawback is its mathematical difficulty.

In this paper, we discuss an approach relying upon Lagrangian duality and stochastic gradient to solve [START_REF] Cohen | Auxiliary problem principle and decomposition of optimization problems[END_REF]. The use of stochastic gradient is based on the reformulation of constraints in probability as constraints in expectation, using an indicator function. As usual with stochastic gradient, we assume that the functions involved in the problem (here, j and θ) are known explicitly but that the probability law governing the "noise" ξ is not, or that the computation of expectations of the variables involved is out of reach or too costly. It is rather assumed that an external mechanism delivers samples of ξ which are used in the iterative algorithm.

Writing the probability as an expectation opens the possibility of using stochastic gradient algorithms, but it also raises the difficulty of handling a discontinuous function, namely the indicator function. We will discuss various ways of overcoming that difficulty.

The rest of the paper is organized as follows. In §2, we present the analysis of the problem, and our resolution strategy, a stochastic Arrow-Hurwicz algorithm. In §3, we describe two structural difficulties of stochastic programming under probability constraint. To implement a stochastic Arrow-Hurwicz algorithm, we need to handle the probability function gradient. In §4, the question we are interested in is therefore: how to compute stochastic estimates of the probability function gradient? In order to answer this question, we propose two methods that allow to obtain biased stochastic gradient estimates, namely Approximation by Convolution (AC) ad Finite Differences (FD). We defer to a forthcoming paper to propose techniques based on integration by parts ideas and providing unbiased (or consistent) estimates, and to compare them with the biased estimates studied hereafter.

We consider a very basic portfolio optimization problem under a probability constraint and use this example throughout the rest of the paper to illustrate and compare the AC and FD techniques. Section 5 is devoted to the convergence analysis of the proposed methods. Finally, §6 reports numerical experiments with the Arrow-Hurwicz algorithm.

2 Analysis of the Problem

Review of Main Difficulties

Probability constraints provide a straightforward risk formulation with an immediate intuitive interpretation. But at the same time, it is well known that such constraints raise important mathematical difficulties, such as the lack of convexity or connectedness of the feasible subset. Indeed, even if θ is a convex function of u for almost all values ξ, the constraint in (3) may not define a convex feasible subset in U (which can even be not connected, if not empty). Those convexity or connectedness (or emptiness) properties depend of course on the properties of θ as a function of its two arguments u and ξ, on the probability distribution of the random variable ξ, on the level α of constraint required and on the level π of probability required. One may refer to [START_REF] Kall | Stochastic Programming[END_REF] for a discussion on those convexity properties, and to [START_REF] Henrion | On the Connectedness of Probabilistic Constraint Sets[END_REF] for connectedness properties.

In [START_REF] Kall | Stochastic Programming[END_REF], the authors prove that if θ(•, •) is jointly convex in (u, ξ) and the probability measure is quasi-concave, then the feasible subset of (3) is convex. But those assumptions seem to us to be rather strong in practice, notably the joint convexity property. Indeed, there are numerous situations in which the decision variable multiplies the random variable, as in the portfolio problem presented in §3, or in a quite other domain, when one wants to model the breakdown of an actuator, in which case the random variable must be able to kill the action the decision variable. In all those situations, the joint convexity property is not realistic.

Mathematical Approach for Programming under Probability Constraint

Before explaining our resolution strategy, we review some basic results on the stochastic Arrow-Hurwciz algorithm [START_REF] Arrow | Studies in nonlinear programming[END_REF][START_REF] Culioli | Decomposition-coordination algorithms in stochastic optimization[END_REF]. First of all, starting with the deterministic constrained optimization problem (1) and assuming that there exists a saddle point of the Lagrangian (2) over U ad × R d + , the (deterministic) Arrow-Hurwicz algorithm consists in performing successive minimization and maximization steps to search for this saddle point:

u k+1 = Π U ad u k -ε k ∇ u J(u k ) + ∇ u Θ(u k )λ k , (4a)
λ k+1 = Π + λ k + ρ k Θ(u k+1 ) -α , (4b)
where Π U ad is the projection onto U ad and Π + is the projection onto the cone R d + .

Stochastic Arrow-Hurwicz Algorithm

The stochastic Arrow-Hurwicz algorithm is typically used to solve a stochastic optimization problem with constraint in expectation:

min u∈U ad E j(u, ξ) s.t. E θ(u, ξ) ≤ α , (5)
where the calculation of expectations is basically difficult if not impossible. The stochastic algorithm overcomes this difficulty by simultaneously approximating the saddle point and the expectations by a Monte-Carlo like technique. It is in fact a combination of the idea of the Monte-Carlo method with the iterative procedure of gradient methods in optimization.

We do assume that a saddle point (u ♯ , λ ♯ ) over U ad × R d + exists for the Lagrangian associated with problem (5) (hence u ♯ is a solution of ( 5)). Observe that this Lagrangian L (2) is equal to the expectation of ℓ(u, λ, •) = j(u, •)+ λ, θ(u, •)-α . We use unbiased estimates of the gradients of L in u and λ obtained with the corresponding gradients of ℓ evaluated at independent drawings ξ k of ξ supposed to follow the probability law of ξ. More specifically, at stage k of the algorithm, u k and λ k being the current estimates of the solution, 1. we draw an independent sample (according to the law P of ξ), or we observe a new independent sample ξ k+1 , 2. we compute the stochastic gradients ∇ u j(u k , ξ k+1 ) and ∇ u θ(u k , ξ k+1 ), 3. we update u k+1 and λ k+1 as follows:

u k+1 = Π U ad u k -ε k ∇ u j(u k , ξ k+1 ) + ∇ u θ(u k , ξ k+1 ) λ k , (6a)
λ k+1 = Π + λ k + ρ k θ(u k+1 , ξ k+1 ) -α . (6b)
Under essentially measurability and convexity assumptions, assuming that the Lagrangian of the problem admits a saddle point, and with

k∈N ε k = +∞ , k∈N (ε k ) 2 < +∞ (and the same for ρ k ) ,
it is shown in [START_REF] Culioli | Optimisation stochastique sous contraintes en espérance[END_REF] that this algorithm converges in the sense that primal {u k } k∈N and dual {λ k } k∈N sequences are bounded a.s. and that {u k } k∈N a.s. weakly converges to some solution u ♯ of (5).

Mathematical Approach: Strategy and Difficulties

From now on, we assume that the critical or risky event is defined by a single (scalar) inequality, that is, θ is R-valued. Let I R + denote the indicator function of the positive half-line. The principle of our resolution strategy is first to replace the probability constraint by a constraint in expectation

-P (u) ≤ -π , (7)
where P (u) = P θ(u, ξ) ≤ α and this probability is evaluated as an expectation:

P θ(u, ξ) ≤ α = E I R + α -θ(u, ξ) , (8)
then resort to duality, and lastly resort to the stochastic Arrow-Hurwicz algorithm. Observe that w.r.t. the general formulation (1), Θ is now -P and the constraint level α is now -π.

There are major difficulties with probability constraints.

• First of all, as we recalled in §2.1, convexity is not preserved. Therefore, existence of a saddle point of the Lagrangian is not granted; in this case, we should resort to augmented Lagrangian techniques to increase the chance that a saddle point does exist. However, this raises new problems because the nonlinearities involved in the augmented Lagrangian formula cannot be combined straightforwardly with expectation to yield obvious stochastic gradient algorithms. This issue of using augmented instead of ordinary Lagrangians goes beyond the scope of this paper and is not considered here.

• We rather address here another difficulty: to replace a probability constraint by a constraint in expectation, we need to handle the indicator function (see ( 8)); but this indicator function involves a discontinuity which may, nevertheless, be smoothed by the expectation operation; however, the stochastic Arrow-Hurwicz algorithm is based on the consideration of a unique sample drawn at each iteration; obtaining a stochastic gradient is therefore not trivial. As it will be shown later on, we propose two ways of overcoming this difficulty: Approximation by Convolution (AC) and Finite Differences (FD). Both approaches will lead us to consider algorithms such as [START_REF] Culioli | Optimisation stochastique sous contraintes en espérance[END_REF] in which either a smooth approximation θ of function θ will be used in both equations (6a) and (6b), or an approximation of its gradient will be used in (6a), leading to a stochastic Arrow-Hurwicz algorithm with biased stochastic estimates of the Lagrangian gradients.

Structural Difficulties of Programming under Probability Constraint

In this section, we focus on two structural difficulties of optimization problem with probability constraints. The first one is related to the non convexity of probability constraint: we show what are the consequences of this non convexity on the stochastic Arrow-Hurwicz algorithm.

The second one concerns the behavior of the probability constraint multiplier in some particular cases.

The Non Convexity of Probability Constraint

Consider the following optimization problem

min u∈R 1 2 (u -1) 2 s.t. P(u ≤ ξ) ≥ π , (9)
where ξ is a normal random variable with mean value -2 and standard deviation 0.1.

In order to point out the first difficulty, we study the qualitative behavior of the underlying deterministic problem, namely that in which the probability constraint is expressed with help of the cumulative distribution function F of ξ: indeed, P(u ≤ ξ) = 1 -F (u) and therefore, the constraint in ( 9) can be replaced by

1 -F (u) ≥ π (10)
without of course altering the corresponding Kuhn-Tucker multiplier.

The Lagrangian of problem [START_REF] Kall | Stochastic Programming[END_REF], with constraint written as in [START_REF] Kushner | Stochastic approximation and recursive algorithms[END_REF], is

L(u, λ) = 1 2 (u -1) 2 + λ F (u) -1 + π ;
and, the Kuhn-Tucker necessary conditions of optimality allow for the calculation of the solution which is, for example with π = 0.7,

u ♯ = -2.05244 and λ ♯ = 0.877913 .
As expected, u ♯ takes the maximal possible value to satisfy the constraint, that is, the (1 -π)-th percentile of the distribution: F (-2.05244) = 0.3, so the constraint is active, and saturated. Figure 1 represents the Lagrangian surface in the (u, λ) domain. For λ = 0, we recognize We insist on the following two points. First, our approach in this paper is based on stochastic estimates of the gradients of the Lagrangian, not on their exact computation, which we assume impossible. Naturally, we cannot expect the stochastic algorithm to behave better than its underlying average driving vector field, which we will study directly. Second, with some probability distributions there is a way to manipulate the constraints in order to preserve convexity. In particular with normal distributions, the map ln(1 -F (•)) is concave [START_REF] Prekopa | Stochastic Programming[END_REF], which leads to a convex formulation of the constraint. If we seem to overlook this remark in the following treatment, this is because the difficulty we try to point out in this very simple case is a fortiori likely to occur in more general situations when the above clever manipulations are no longer possible: recall that we do not assume knowledge of the distribution of θ(u, ξ).

Let us now consider the ODE associated with the Arrow-Hurwicz algorithm,

u = -L ′ u (u, λ) = -J ′ (u) -λ F ′ (u) , (11a) λ = L ′ λ (u, λ) = F (u) -1 + π . ( 11b
)
At u ♯ = 1, the unconstrained optimal solution, one has that J ′ (u ♯ ) = 0 and F ′ (u ♯ ) = 1.47 × 10 -195 , because u ♯ happens to be in the tail of the distribution. Therefore, even for very large values of λ, L ′ u (u ♯ , λ) remains very close to 0; in other words, if the (continuous) algorithm ( 11) is started at (or close to) (u ♯ , λ), for practically any λ, u will stay at u ♯ ! At the same time, if u ♯ doesn't satisfy the constraint, one has that F (u ♯ ) > 1 -π. It follows that L ′ λ (u ♯ , λ) > 0: λ increases almost indefinitely! This is better illustrated by the vector field of the ODE, shown in Figure 2. The white zone corresponds to the basin of attraction of the optimal solution. In the grey zone, the algorithm is driven more or less indefinitely towards large values of λ in a valley corresponding to the unconstrained solution u ♯ = 1.

This example shows that even a deterministic algorithm may, if started on the "wrong" side, wander away from the actual optimal solution. Stochastic versions of the algorithm are expected to behave erratically, and even if the current values of (u k , λ k ) are in the basin of attraction of the Kuhn-Tucker point, random observations may take the algorithm to other regions away from the optimal solution.

Degeneracy of the Probability Constraint Multiplier

Let us now consider the following portfolio optimization problem. This very simple problem allows us to point up another structural difficulty of probability constraint. This example will also be used in the remainder of this paper to illustrate our various approaches.

We borrow a capital which we have to pay off at the end of the period with an interest rate l. We can invest a proportion u of this capital at the fixed rate b, invest a proportion v at the random rate ξ, and finally consume the available remainder, which brings a satisfaction measured by a concave nondecreasing function f . We assume of course that E(ξ) > l, in other words, risk is rewarding. We try to maximize the sum of the satisfaction provided by consumption and by the expected final capital. We also want to be in a position to pay off the capital and the interests at the end of the period, with a probability of a least p. In this case, the optimization problem can be stated as follows:

max u,v E f (1 -u -v) + (1 + b)u + (1 + ξ)v s.t. u ≥ 0 , v ≥ 0 , u + v ≤ 1 , P (1 + b)u + (1 + ξ)v ≥ 1 + l ≥ π . Let l = 0.15, b = 0.2, f (x) = -x 2 /2 + 2x , F (ξ) =        0 if ξ < ξ -σ , 1 16 3 ξ- ξ σ 5 -10 ξ- ξ σ 3 + 15 ξ- ξ σ + 8 if ξ < ξ + σ , 1 otherwise , ( 12
)
where F is the distribution function. For numerical experiments, we set ξ = 0.4 and σ = 3. To identify this problem with (3), consider the equivalent minimization problem with cost function

j(u, v, ξ) = -f (1 -u -v) -(1 + b) u -(1 + ξ) v .
Let also

P (u, v) = P (1 + b) u + (1 + ξ) v ≥ 1 + l . ( 13
)
This problem is now formulated as min u≥0,v≥0 this function is not convex. In fact, it is convex for probability levels below 0.57. For probability levels close to 0.57, the risk of not being in a position to pay off the capital and the interests is important; the investment in the risky asset v decreases to zero, whereas simultaneously, that in the secure asset u increases. The optimal cost, which was until then a convex function of the required probability level, becomes concave. Above 0.65, v is zero, u is equal to (1 + l)/(1 + b) = 0.95833 in order to satisfy the probability constraint, and the optimal cost becomes constant. This example shows another structural difficulty of optimization under a probability constraint, namely the degeneracy of the probability constraint multiplier. Indeed, for π small enough, the secure asset, u is zero at optimum. For π large enough, the risky asset, v, is zero at optimum. In the latter case, the event (1 + b) u + (1 + ξ)v ≥ 1 + l can only have probability 0 or 1. That is to say, at v = 0 and u = (1 + l)/(1 + b) = 0.95833, the function (13) exhibits a discontinuity (see Figure 4). When π is large and v ♯ = 0, the probability (13) can only take values 0 or 1 (depending on the value of u), that is, this probability is strictly larger than the required level π when the constraint is met. Clearly the constraint is not "saturated", because there is no equality, and consequently the corresponding multiplier is zero; small changes in π will not affect the solution. However, the constraint is "active", that is, the solution (u ♯ , v ♯ ) of the problem with the probability constraint is different from the solution (u * , v * ) without it.

E j(u, v, ξ) s.t. u + v ≤ 1, -P (u, v) ≤ -π with Lagrangian L(u, v; λ 1 , λ 2 ) = E j(u, v, ξ) + λ 1 (u + v -1) + λ 2 (π -P (u, v)) .
In the remainder, in order to guarantee existence of a saddle point of the Lagrangian, we consider only probability levels below 0.57 (otherwise, one should resort to augmented Lagrangian techniques, but, as it has been said earlier, this issue goes beyond the scope of this paper). For example, for a probability level of 0.24, the primal-dual optimal solution is

u ♯ = 0 , v ♯ = 0.50407 , λ ♯ 1 = 0 , λ ♯ 2 = 0.08815 . ( 14
)

Stochastic Estimates of Probability Function Gradient

As it was mentioned at §2.2.2, in order to use a stochastic Arrow-Hurwicz algorithm, we need to handle the probability function gradient, that is, to obtain a stochastic estimate of the gradient of (see [START_REF] Henrion | On the Connectedness of Probabilistic Constraint Sets[END_REF])

P (u) = E I R + α -θ(u, ξ) . (15)
It is well known that this gradient is difficult to compute; we may refer to [START_REF] Uryasev | Introduction to the Theory of Probabilistic Functions and Percentiles (Value-at-Risk)[END_REF] for a discussion of this topic. Recall that in our case, replacing the probability constraint by a constraint in expectation raises the difficulty of handling an indicator function, which is a discontinuous function. One way of dealing with this problem is to appeal to a technique based on convolution to derive a smooth approximation of this discontinuous function. Alternatively, we can obtain a stochastic estimate of the gradient of this function, based on a single sample drawing of ξ, by appealing to a finite difference technique, and we rely upon the multiplication of such drawings along the iterative algorithm to smooth out this crude estimate.

Approximation by Convolution Method (AC)

General Theory

The basic principle of this approach is to smooth out the indicator function appearing in (15) so that differentiation underneath expectation becomes possible. Consider a smooth function h : R → R with the following properties : h as a unique maximum at x = 0,

∀x, h(x) ≥ 0; h(x) = h(-x); +∞ -∞ h(x) dx = 1 . (16)
We will give a few examples of such functions later on and will consider only functions with finite support although this is not absolutely necessary. With any other function φ : R → R, and r a small positive number, the convolution

φ r (x) = 1 r ∞ -∞ φ(y)h x -y r dy ,
can be viewed as an approximation of φ since h(•/r)/r approximates the Dirac function (in the sense of convergence of distributions) at 0 when r tends to zero. The function φ r is differentiable with

φ ′ r (x) = 1 r 2 ∞ -∞ φ(y) h ′ x -y r dy .
This technique is widely known as the "mollifier" technique [START_REF] Ermoliev | The Minimization of Semicontinuous Functions, Mollifier Subgradients[END_REF]. We now apply it to I R + : recall (15) and define

P r (u) = 1 r E +∞ -∞ I R + (y) h α -θ(u, ξ) -y r dy = 1 r E +∞ 0 h y -α + θ(u, ξ) r dy
(here, we have used the fact that h is an even function)

= E p r (u, ξ) (17)
with

p r (u, ξ) = 1 r +∞ 0 h y -α + θ(u, ξ) r dy . (18)
Then

(p r ) ′ u (u, ξ) = 1 r 2 θ ′ u (u, ξ) +∞ 0 h ′ y -α + θ(u, ξ) r dy = - 1 r h θ(u, ξ) -α r θ ′ u (u, ξ) , (19)
and clearly

P ′ r (u) = E (p r ) ′ u (u, ξ) . (20)
Therefore, for any sample ξ, (19) can be considered as a stochastic estimate of P ′ (u), albeit a biased one; however, this bias vanishes when r approaches 0. In what follows, we evaluate the bias and the variance of this estimate as a function of r.

Remark 1. In the same way, according to (17), ( 18) can be considered a biased estimate of

P (u) whereas p(u, ξ) = I R + α -θ(u, ξ) (21)
is an unbiased one. In Equation (6b) of the iterative algorithm, we may either use the unbiased estimate or the biased one, consistently with that used in (6a) for Θ ′ . The latter option has the advantage of preserving the specific geometry of vector fields of Arrow-Hurwicz algorithms (with some symmetry, or skew-symmetry, properties, according to the point of view). The former option may seem preferable as long as it avoids seemingly unnecessary bias or approximation. Both options will be tested later on in §6. Therefore, the next theorem deals with both the estimates (18) and (19) in order to cover all variants.

Theorem 2. The random variable (or vector) ξ is supposed to admit a density q(ξ). For the random variable X u (•) = θ(u, •) depending on the parameter u, we assume that the induced probability law also admits a density denoted q Xu (x) and that this density is at least twice continuously differentiable with L 1 first and second order derivatives. Then, for any sample drawing ξ following the probability density q, the expression (18) provides a biased estimate of P (u) with a bias in O(r 2 ) and a variance in O(1).

For the pair of random variables (or vectors)

X u (•), Y u (•) = θ(u, •), θ ′ u (u,
•) depending on the parameter u, we assume that the induced joint probability law admits a density denoted q XuYu (x, y) and that this density is at least twice continuously differentiable in x with integrable L 1 first and second order derivatives. Then, for any sample drawing ξ following the probability density q, the expression (19) provides a biased estimate of P ′ (u) with a bias in O(r 2 ) and a variance in O(1/r).

Proof. Consider first (17)-(18)

. With the induced probability law for the random variable X u , one has that

P r (u) = 1 r +∞ -∞ +∞ 0 h y -α + x r q Xu (x) dy dx .
Using Fubini theorem and the change of variable z = (y -α + x)/r in the integral in x yields

P r (u) = +∞ 0 +∞ -∞ h(z) q Xu (rz -y + α) dz dy .
With the smoothness assumptions on q Xu , the Taylor expansion of this term for r near 0 yields

P r (u) = +∞ 0 +∞ -∞ h(z) q Xu (α -y) + rz q ′ Xu (α -y) + r 2 z 2 2 q ′′ Xu (α -y) + O(r 3 )z 3 dz dy = +∞ 0 q Xu (α -y) dy + r 2 2 σ 2 h +∞ 0 q ′′ Xu (α -y) dy + O(r 3 )
by using (16) on the one hand and by introducing

σ 2 h = +∞ -∞ z 2 h(z) dz (22)
on the other hand. The term of order 0 in r can be written as

α -∞
q Xu (t) dt and, as such, is recognized to be equal to P(X u ≤ α), that is, P (u). Therefore, P r (u) differs from P (u) by an O(r 2 ) term (proportional to σ 2 h ). As for the variance of the estimate (18), it is equal to the second order moment E p r (u, ξ) 2 from which the square of E p r (u, ξ) must be subtracted. The latter is close to P (u) 2 up to a term of order O(r 2 ). Therefore we concentrate on the second order moment which can be written, according to (18),

E p r (u, ξ) 2 = 1 r 2 E +∞ 0 h y -α + θ(u, ξ) r dy 2 = 1 r 2 +∞ -∞ +∞ 0 h y -α + x r dy 2 q Xu (x) dx = +∞ -∞ +∞ x-α r h(z) dz 2 q Xu (x) dx
using the change of variable z = (y -α + x)/r in the integral in y,

≤ +∞ -∞ +∞ -∞ h(z) dz 2 q Xu (x) dx since h(•) ≥ 0, = 1
according to (16) (last equality) and the fact that q Xu is a probability density. The proof regarding the bias of P ′ r (u) w.r.t. P ′ (u) may follow one of the following two paths: either a similar result is proved for the derivative of a function whenever the function itself is approximated by another function up to an O(r 2 ) term; or, with (19)-(20), we perform similar calculations to those we have just performed with (17)-(18). Let us sketch this second path. Considering ( 19)-(20) and the pair X u (•), Y u (•) , we have that

P ′ r (u) = - 1 r h
x -α r y q XuYu (x, y) dx dy (remember y may be a vector of the same dimension as u and dy should be understood adequately). From here, we proceed as previously with the change of variable z = (x -α)/r which yields P ′ r (u) = -h(z) y q XuYu (rz + α, y) dz dy .

Then, a Taylor expansion of q XuYu w.r.t. its first argument for r near 0 yields, for the same reasons as previously,

P ′ r (u) = -y q XuYu (α, y) dy + r 2 2 σ 2 h ∂ 2 q XuYu (α, y) ∂x 2 y dy + O(r 3 ) . ( 23
)
Assuming that the first term in the right-hand side above is equal to P ′ (u) (see Claim 3 hereafter), we obtain again that P ′ r (u) differs by an O(r 2 ) term. The variance is equal to the second order moment E p ′ r (u, ξ)

2 from which we must subtract E p ′ r (u, ξ) 2
. The latter is close to P (u) 2 up to O(r 2 ). As for the former, we have that

E p ′ r (u, ξ) 2 = 1 r 2 h 2 θ(u, ξ) -α r θ ′ u (u, ξ) 2 q(ξ) dξ = 1 r 2 h 2 x -α r y 2 q XuYu (x, y) dx dy = 1 r h 2 (z) y 2 q XuYu (rz + α, y) dz dy .
From here, we proceed as earlier with a Taylor expansion for r close to 0, and it should be clear that the above expression is of order 1/r with a coefficient which can be bounded by a term proportional to the square of the L 2 norm of h. The same consideration is still valid for the variance itself.

Claim 3. The first term in the right-hand side of ( 23) is equal to P ′ (u). We sketch the proof of this fact here. For any smooth function f : R → R, consider

F (u) = E f θ(u, ξ) = f θ(u, ξ) q(ξ) dξ = f (x) q Xu (x) dx .
Then,

F ′ (u) = f ′ θ(u, ξ) θ ′ u (u, ξ) q(ξ) dξ = f ′ (x) y q XuYu (x, y) dx dy .
Integrating by parts in the integral in x, one gets

F ′ (u) = - f (x) y ∂q XuYu ∂x (x, y) dx dy .
If f is not smooth enough for this calculation to be immediately justified, one can consider a sequence of smooth approximations converging to f in order to establish this formula. Let us now use it for f (•) = I R + (α -•). Then F (u) = P (u) and

P ′ (u) = -y I R + (α -x) ∂q XuYu ∂x (x, y) dx dy = -y α -∞
∂q XuYu ∂x (x, y) dx dy = -y q XuYu (α, y) dy , which is exactly the expected result.

Remark 4. As a side remark, observe that q Xu (α) = q XuYu (α, y) dy , and that q XuYu (α, y)/q Xu (α) = q Yu (y

| X u = α) ,
that is, the conditional density of Y u knowing that X u = α. Therefore, the first term in the right-hand side of ( 23) can be written as -E Y u | X u = α × q Xu (α). We conclude that

P ′ (u) = -q θ(u,•) (α) × E θ ′ u (u, •) | θ(u, •) = α .
Remark 5. Observe that, although we started with the idea of a smooth function h, the expression (19) of the estimate and the analysis in the proof of Theorem 2 does not involve more than the function h itself (not its derivatives), so that we may as well consider non smooth functions (and even discontinuous functions at the ends of its support).

In conclusion, the variance of the stochastic estimate (19) blows up like A/r as r goes to 0 (where A can be bounded from above by something proportional to the square of the L 2 norm of h), that of (18) remains of order O(1), whereas the square of the bias of both estimates goes to 0 as B 2 r 4 (where B is proportional to σ 2 h -see ( 22)). If the estimate of P ′ (u) is rather based on the average of N expressions as (19) for N independent drawings of ξ, the variance will blow up as A/(N r) whereas the square of the bias will still behave as B 2 r 4 . Therefore, the best trade-off between variance and bias is realized by that r which minimizes the mean square error (MQE; this is the sum of the variance and of the square of the bias) equal to A/(N r) + B 2 r 4 : the "best" r is thus A/(4B 2 N ) 1/5 . This yields a MQE estimated to (5A 4/5 B 2/5 )/(4N ) 4/5 . Therefore, in the choice of function h, it is meaningful to pay attention to the quantity σ

4/5 h h 8/5 L 2 . Remember B is proportional to σ 2
h and A is proportional to h 2 L 2 . The bias of the AC estimate goes to 0 with r: this parameter r allows for a trade-off between mean and variance which should be adapted to the number of samples available (as just discussed) or visited in one run in the context of an iterative algorithm, as discussed later on in §5.3.

Practical Aspects and Application to Example of §3.2

Define

I(x) = 1 if -1 ≤ x ≤ 1, 0 otherwise.
Table 1 proposes 6 functions with bounded supports that can play the role of function h (see ( 16)) and compares them from the point of view of their constants σ

4/5 h h 8/5
L 2 (last column), the relevance of which has just been explained. The column h(0) is provided to help identifying the functions with their graphs displayed in Figure 5.

We observe that the fourth function, namely h(x) = 3(1 -x 2 )I(x)/4 is the one to retain because it offers the smallest value in the last column of the table. We now apply the technique to the example of §3.2 again. The estimates for P ′ u (u, v) and P ′ v (u, v) based on this technique and on a given sample ξ read as follows: This MQE will be compared with that obtained by the next approach, namely finite differences. We performed some exact computations of bias and variance with the help of Mathematica for those estimates evaluated at the optimal solution (see [START_REF] Uryasev | Introduction to the Theory of Probabilistic Functions and Percentiles (Value-at-Risk)[END_REF]) and with h equal to the fourth function in Table 1. We have found:

(p r ) ′ u (u, v, ξ) = 1 + b r h (1 + α) -(1 + b)u -(1 + ξ)v r ; ( 24a
) (p r ) ′ v (u, v, ξ) = 1 + ξ r h (1 + α) -(1 + b)u -(1 + ξ)v r . ( 24b
) h h(0) σ 2 h h 2 L 2 σ 4/
E(p r ) ′ u (u ♯ , v ♯ , •) = 0.62 -0.096r 2 + 0.012r 4 , Var(p r ) ′ u (u ♯ , v ♯ , •) = 0.45/r -0.39 -0.05r + O(r 2 ) , E(p r ) ′ v (u ♯ , v ♯ , •) = 1.18 -0.36r 2 + 0.06r 4 , Var(p r ) ′ v (u ♯ , v ♯ , •) = 1.62/r -1.39 -0.35r + O(r 2 ) .
If the estimates are based on the average over N samples, the MQE of the AC estimates are obtained by considering Var(r)/N + (E(r)) 2 -(E(0)) 2 . In those expressions, we consider the terms in 1/N r and r 4 only in order to tune r as a function of N . This computation is done for the sum of the MQE's related to the two components of p ′ r (that is, for the mean square norm of the vector estimate error -we denote it MQE(r, N )). This yields r = 1.30 N -1/5 . Finally, we plug this value of r into MQE(r, N ) to get the following function of N (again, calculations are exact, up to Mathematica accuracy, even if results are displayed in a truncated form):

1.98 N 4/5 - 1.78 N + O(N -6/5 ) . (25)
those two events being of course disjoint. Using the mean value theorem (or Taylor representation) for the function θ(•, ξ) ∈ C 1 , we have that θ(u + c1 j , ξ) = θ(u, ξ) + cθ ′ u j v + (ξ), ξ , and similarly θ(u

-c1 j , ξ) = θ(u, ξ) - cθ ′ u j v -(ξ), ξ . Therefore, E ∇ c u j p(u, ξ) 2 = 1 4c 2 P {α + cθ ′ u j v -(ξ), ξ < θ(u, ξ) ≤ α -cθ ′ u j v + (ξ), ξ } + P {α -cθ ′ u j v + (ξ), ξ < θ(u, ξ) ≤ α + cθ ′ u j v -(ξ), ξ } , (27)
Thanks to (H1), we can bound each of the above two probabilities by

P {θ(u, ξ) ∈ (α -Kc, α + Kc]} , ( 28
)
where K is the uniform bound on θ ′ u j . Our goal is now to evaluate the behavior of this probability when c is approaching 0. Let m be the dimension of ξ; g is supposed to be R-valued. Consider any solution ξ of

θ(u, ξ) = α . (29)
The case when no such solution exists for some u will be discussed later on at Remark 7. If θ ′ ξ (u, ξ) = 0 as assumed in (H3), then the manifold of solutions of ( 29) is locally of dimension less than or equal to m -1. The set of ξ's involved in the event in (28) is locally a set with a "backbone" given by this manifold around ξ, and a "thickness" which is proved to be of order O(c). Indeed, with a Taylor expansion of θ(u, •) around ξ, we get

θ(u, ξ + y) = α + θ ′ ξ (u, ξ), y + O( y 2 ) .
In this expression, we need only consider y's which are (asymptotically as c goes to 0) parallel to the gradient θ ′ ξ (u, ξ) (that is, the component in the kernel of the linear form defined by this gradient is useless). It should now be obvious that to match variations of g around α which are of order c, we need only consider y's which are also of order c in norm. If this holds true for any ξ in the manifold of solutions of (29), then the probability (28) is of order O(c) and the second order moment (27) of our estimate (and consequently the variance too) is bounded by an O(c -1 ).

If (H3) does not hold but (H4) does, then the same reasoning can be repeated (for a Taylor expansion of the next order) with y's which are now orthogonal to the kernel of the Hessian θ ′′ ξ 2 (u, ξ) (this component is non zero thanks to (H4)) and it should be clear that to compensate for variations of g of order c, we now need y's which are of order O(c 1/2 ) in norm. This also gives the order of the probability (28) and then, the bound on the variance is in O(c -3/2 ).

We could continue like that by removing assumption (H4) but introducing an assumption (H5), and so on and so forth. Ultimately, with no particular assumptions, (28) is of order O(1) and the variance is of order O(c -2 ).

Remark 7. Suppose that for some u, there exists no solutions to (29). Then, since g is assumed to be at least continuous in ξ, this means that for all ξ, θ(u, ξ) is always either strictly less or strictly greater than α, in which cases P (u) (see (15)) assumes either the value 1 or 0 (which are extreme values for P ).

If θ(u, •) can be bounded away from α, then the probability (28) will be 0 for c small enough. This is the good case for the variance of the estimate. But θ(u, •) may also approach α asymptotically, and, with heavy tails for the density q of ξ, it is not possible to give a better bound for (28) than O(1). Here is an example. Let θ(u, ξ) = u -e -ξ (u and ξ are both scalar). Consider the probability P(θ(0, ξ)) ∈ [-c, c] for c small, that is, P(ξ ≥ -ln c). Assume the density q(ξ) has a positive support and that it is equal to a(1 + ξ) -(1+a) I R + (ξ) with a an arbitrary small positive number. Then P(ξ ≥ -ln c) = (1 -ln c) -a . For c positive and below e, (1 -ln c) -1 ≥ c, hence this probability is larger than c a . Since a is positive and arbitrarily small, we cannot clearly make this case enter the case of better bounds obtained with assumptions (H3) or (H4).

Application to the Example and Comparison with the AC Method

We have used (26) (for the two components of the gradient, that in u and that in v) to our example and evaluated, once again with the help of Mathematica, the mean and variance of those estimates at the optimal solution [START_REF] Uryasev | Introduction to the Theory of Probabilistic Functions and Percentiles (Value-at-Risk)[END_REF]. The results are as follows:

E ∇ c u p(u ♯ , v ♯ , •) = 0.62 -0.23c 2 + 0.06c 4 + O(c 7 ) , Var ∇ c u p(u ♯ , v ♯ , •) = 0.31 c -0.39 -0.12c + O(c 2 ) , E ∇ c v p(u ♯ , v ♯ , •) = 1.18 -1.49c 2 -42.25c 4 -199.41c 6 + O(c 7 ) , Var ∇ c v p(u ♯ , v ♯ , •) = 0.59 c -0.39 -0.74c + O(c 2 ) .
Following the same procedure as for the AC estimate, the MQE for the gradient vector estimate based on N independent samples is obtained by Var(c)/N + (E(c)) 2 -(E(0)) 2 ; in this expression, the dominant terms in 1/N c and in c 4 only are retained to tune c as a function of N . This yields c = 0.63N -1/5 and an optimal MQE equal to:

1.79

N 4/5 - 1.78 N + O(N -6/5 ) . (30)
Compared with (25) which was obtained with the AC estimate, this is asymptotically slightly better. However, a more careful inspection with complete expressions of the MQEs shows that this conclusion becomes true only for N above about 11000. Hence one may say that the AC and the FD methods yield approximately the same performances.

5 Convergence Analysis

Stochastic Algorithms

Consider algorithm [START_REF] Culioli | Optimisation stochastique sous contraintes en espérance[END_REF]. With Θ(u) = E θ(u, ξ) and J(u) = Ej(u, ξ), an equilibrium point (u ♯ , λ ♯ ) of this algorithm solves the system of Kuhn-Tucker optimality conditions of problem (1): for all positive ε and ρ,

u ♯ = Π U ad u ♯ -ε ∇ u J(u ♯ ) + ∇ u Θ(u ♯ ) λ ♯ , (31a)
λ ♯ = Π + λ ♯ + ρ Θ(u ♯ ) -α . (31b)
We will write algorithm (6) (with ρ k proportional to ε k ) more compactly: we set x = (u, λ) and write

x k+1 = Π(x k -ε k ψ k ) , (32)
Lemma 8. Consider the iteration (32) and assume that

k ε k = +∞ , (37a)
k ε k B k < ∞ a.s., (37b) k
(ε k ) 2 V k < ∞ . (37c)
Then, a.s., x k has the same asymptotic behavior as the solution of (36).

This result follows from [START_REF] Kushner | Stochastic approximation and recursive algorithms[END_REF]Chap. 5].

Proposition 9. Consider the case when the estimate (19) (and possibly (18) too) is (are) used in the stochastic algorithm (6) (with ρ k proportional to ε k ) with the following choices of the stepsize ε k and of the "mollifier" parameter r k :

ε k = k -γ , r k = k -β/2 , ( 38
)
for β and γ positive. Then the conditions of Lemma 8 are satisfied if

γ ≤ 1, β + γ > 1, 2γ -β/2 > 1 . (39)
Proof. The first condition (39) is required by (37a). Theorem 2 states that the bias β+γ) ; therefore (37b) is satisfied under the second condition (39). As for the variance -2γ ; thus (37c) is satisfied under the third condition (39).

B k of AC estimates is in O (r k ) 2 = O k -β , hence ε k B k = O k -(
V k , it is in O (r k ) -1 = O k β/2 which yields (ε k ) 2 V k = O k β/2
Proposition 10. Consider the case when the estimate (26) is used in (6a) (with ρ k in (6b) proportional to ε k ) with the following choices of the stepsize ε k and of the FD parameter c k :

ε k = k -γ , c k = k -β/2 , (40)
for β and γ positive. Then the conditions of Lemma 8 are satisfied if, in addition of assumptions (H1) and (H2) of Theorem 6, one has that

γ ≤ 1, β + γ > 1,      2γ -β/2 > 1 if (H3) is satisfied in Theorem 6, 2γ -3β/4 > 1 if (H4) is satisfied in Theorem 6, 2γ -β > 1 otherwise. (41)
The proof follows the same pattern as previously using the evaluations of Theorem 6, the only changes concerning V k .

Convergence Rate

Let β and δ be the integers such that:

B k = O(k -β ), V k = O(k -δ ) . (42)
Reference [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF] provides a comprehensive analysis of the convergence rates of algorithms of type (32) under the following assumption: close to its unique equilibrium point x ♯ (supposed to lie in towards a normal distribution of mean 0 and covariance Σ where:

B = lim k→∞ k β B k , H b = A -βI , H = A -(1 + δ)/2 I , R = lim k→∞ k δ E(∆M k (∆M k ) ⊤ | F k ) , ΣH + H ⊤ Σ = R ,
where ⊤ denotes transposition.

Remark 12. From the definition of H b above and the appearance of A -(1 + δ)/2 I in the definition of Σ, it is apparent that the strong stability condition (44) (case γ = 1) ensures that both these matrices are positive definite, so that Σ is well defined. Indeed, with our choices, H b and H are equal.

The Case of Arrow-Hurwicz Algorithms

We know discuss the properties of matrix A in the situation of Arrow-Hurwicz This matrix has been introduced in (43) in general, and the operator Ψ is defined by (33) in our case. Thus, A is the linearized version of that Ψ at the equilibrium point x ♯ , that is,

A = ∂ 2 L(u ♯ ,λ ♯ ) ∂u 2 ∂ 2 L(u ♯ ,λ ♯ ) ∂u ∂λ -∂ 2 L(u ♯ ,λ ♯ ) ∂λ ∂u -∂ 2 L(u ♯ ,λ ♯ ) ∂λ 2 = J ′′ (u ♯ ) + λ ♯ ⊤ Θ ′′ (u ♯ ) Θ ′ (u ♯ ) ⊤ -Θ ′ (u ♯ ) 0 (46)
However, among the constraints Θ, only those saturated (that is, satisfied with equality) at the equilibrium point should be taken into account together with their corresponding multipliers (that is, the non saturated constraints are virtually absent asymptotically). Under the assumptions that the gradients of saturated constraints are linearly independent (or, otherwise stated, the operator in the upper right-hand corner of the matrix is injective), and that the Hessian of the Lagrangian (that is, the operator in the upper left-hand corner) is positive definite, it can easily be proved that the real part of the eigenvalues of A are positive (see [START_REF] Bertsekas | Nonlinear Programming. 2nd Edition[END_REF]proof of Proposition 4.4.2]). This is condition (44) in the case γ < 1. When γ = 1, condition (44) is stronger and will be discussed shortly in the case of our example. Observe that if we assume that the only saturated dualized constraint is the probability constraint (which is the case in our example), then we should assume that the gradient of this probability function at the equilibrium is not zero.

Going back to example of §3. is not satisfied (with β = (1 + δ)/2 = 2/5). However, in the same way as we ignored multipliers corresponding to non saturated constraints because they are stuck to 0 asymptotically, we may consider that the part u of primal variables is "out of the game" ultimately because u is stuck to 0 (the constraint u ≥ 0 is saturated) at the end of the transient part of the algorithm (remember that the ODE (36a) is to be replaced by the more complex dynamics (36b) when following boundaries of the admissible domain). Therefore, we consider a reduced matrix A by keeping only the 2 × 2 lower right-hand block (corresponding to the pair (v, λ 2 )). The eigenvalues of this reduced matrix are 0.605 ± 1.014 i and now condition (44) is satisfied even for the case γ = 1.

Numerical Results

Algorithm [START_REF] Culioli | Optimisation stochastique sous contraintes en espérance[END_REF] has been used to solve the example of §3.2 with the AC and FD estimates.

u k+1 = Π U ad u k -ε k ∇ u j(u k , ξ k+1 ) -∇ u P (u k , ξ k+1 ) λ k , (47a)
λ k+1 = Π + λ k + ρ k π -P (u k+1 , ξ k+1 ) . (47b)
More precisely, for the AC method, ∇ u P (u, ξ) should be interpreted as the gradient estimate (19), applied to the example (see ( 24)); P (u, ξ) is either p(u, ξ) as in (21) or the biased estimate given by (18) (see Remark 1). We tested both versions numerically and there was no significant difference. The estimate (18) was retained for the rest of experiments. Of course, parameter r k is adjusted according to the rule r k = ak -1/5 where a is a positive constant to be tuned.

For the FD method, ∇ u P (u, ξ) is given by (26), applied to the example. Again, parameter c k is adjusted as c k = bk -1/5 where b is a positive constant to be tuned. For P (u, ξ), we used (21).

Numerical experiments are performed according to the following protocol:

• all runs of the algorithms start from the same initial conditions:

u 0 = 0.2, v 0 = 0.8, λ 0 1 = 0.5, λ 0 2 = 0.3 .
Recall that the solution is given by ( 14) and all results will be expressed in terms of differences with those optimal values (hence the equilibrium point for all variables is at 0).

• For AC and FD, 100 runs of the algorithms are performed using the same 100 sequences of pseudo-random numbers to generate Monte Carlo samples of ξ according to the distribution of this variable.

• 5000 thousands iterations are performed for each run.

• For AC and FD, averages of the differences x k -x ♯ are computed over the 100 runs together with their standard deviations. What will be shown on the plots are the trajectories of the "average ± standard deviation" of those quantities as functions of the iteration index k.

• The parameters a, b, d, e, f, g appearing in the following rules:

r k = a k 1/5 , c k = b k 1/5 , ε k = d e + k , ρ k = f g + k ,
are tuned by some trials to try to obtain the "best" results for both methods.

Figure 6 shows the plots for the four variables and for the AC (continuous line) and DF (dotted line) methods. Again what is displayed is the "average ± standard deviation" over 100 runs. Results obtained on this example are very close (with maybe a slight advantage to FD in the earliest iterations) with both methods. This confirm the estimation of variance and bias made with Mathematica around the optimum for the estimates obtained with the two methods.

Conclusions

This paper discussed the problem of stochastic optimization under probability constraints and in particular methods for solving them numerically. Although there exist other ways of taking care of risk considerations in decision problems under uncertainty, we discussed the fact Unfortunately, as shown by the discussion and examples in §3, probability constraints may be the source of several pathologies, and the loss of convexity is the most frequent one. Nevertheless, one must address the problem of numerical resolution with approaches which may fail in the worst cases but which may also succeed to solve nontrivial problems. Our strategy is based on duality and stochastic gradient algorithms. Duality, and the use of stochastic Arrow-Hurwicz algorithms, require the existence of a saddle point of the Lagrangian, which is not granted for the reasons advocated above. The use of augmented Lagrangians would certainly increase the chance of existence of saddle points but, in combination with stochastic algorithms, it raises new difficulties (namely, the operator of mathematical expectation would appear inside a nonlinear function). This new topic will be addressed in a forthcoming paper.

Apart from this problem of saddle point existence, the search of this saddle point by stochastic gradient algorithms is made possible by expressing the probability constraint as an expectation involving a discontinuous function. In this paper, we proposed two ways to overcome this difficulty, and we studied the convergence and convergence rate of the resulting algorithms. The two methods provide biased stochastic estimates of the constraint gradient. Although their implementation on a simple example showed a similar behavior, the theoretical results reveal that in more general situations, the "mollifier" (or "Approximation by Convolution" -AC) method should be of more general use and robustness then the "Finite Difference" (FD) method. We defer to a forthcoming paper to propose other estimation techniques providing unbiased estimates and based on techniques of integration by parts.

Still, the surface of this difficult field of numerical resolution of probability constrained stochastic optimization problems has been just scratched here, and several directions remain open for future investigations. For example, we have considered here only events (whose probability is constrained) which are described only by a scalar constraint and the case of events described by multidimensional constraints may raise new questions (although the techniques discussed in the present paper seem ready for an extension to this case).

Figure 1 :

 1 Figure 1: Lagrangian surface

Figure 2 :

 2 Figure 2: Vector field of the ODE

Figure 3 Figure 3 :

 33 Figure3represents the optimal cost as a function of probability level π. We observe that

Figure 4 :

 4 Figure 4: Graph of the probability function for (u, v) ∈ [0.9, 1] × [0, 0.1]

Figure 5 :

 5 Figure 5: Several possible h functions

Figure 6 :

 6 Figure 6: Average ± standard deviation for AC (solid line) and FD (dotted line) algorithms

Table 1 :

 1 Various h functions

	h	5	h	8/5 L 2

  2, matrix A (restricted to the variables (u, v, λ 2 )) is equal to 

		0.944 1.002 -0.621	
		1.002 1.211 -1.181	
		0.621 1.181	0
	with eigenvalues 0.974 ± 0.753 i and 0.207. As predicted, the real parts are positive but the
	smallest one is equal to 0.207 which is not greater than 2/5. Thus, condition (44) (case γ = 1)

Finite Differences (FD)

General Theory

The idea here is simply to evaluate the derivative w.r.t. each component u j of the expression inside expectation in (15) by the variation of this quantity, caused by, and divided by, the symmetric variation (u j + c) -(u j -c) = 2c for a sample ξ. We denote 1 j the vector of the same dimension as u with a 1 in the j-th component and 0 elsewhere. The FD stochastic estimate of P ′ u j is

It is wise to use the same sample ξ for the evaluation at u + c and u -c in order to reduce variance. A symmetric difference around u is also recommended. We notice that [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF] includes FD under a.s. continuity assumptions, which is not our case here, because the indicator functions are discontinuous.

The following theorem provides the analysis of bias and variance of this estimate w.r.t. the parameter c. Theorem 6. If P (see (15)) is three times continuously differentiable with bounded derivatives, the expression (26) provides a biased estimate of P ′ u j with a bias in O(c 2 ). If (H1) θ(•, ξ) is differentiable with derivatives bounded uniformly in ξ;

(H2) the probability measure of ξ has a density;

(H3) θ(u, •) is twice differentiable and, for all u, and for every solution ξ of θ(u, ξ) = α, we have that θ ′ ξ (u, ξ) = 0; then the variance of estimate (26) is in O(c -1 ).

If (H1) and (H2) still hold true but (H3) is replaced by (H4) θ(u, •) is three times differentiable and, whenever θ(u, ξ) = α for some ξ, and θ ′ ξ (u, ξ) = 0, we have that θ ′′ ξ 2 (u, ξ) = 0;

then the variance of estimate (26) is in O(c -3/2 ). Finally, under no particular assumptions on g, the best bound for the variance is in O(c -2 ).

Proof. With the smoothness assumption on P , one has that

which proves the claim on the bias.

To evaluate the variance of (26), we study its second order moment which differs from the variance by P ′ u j (u) 2 up to terms in O(c 2 ) as we have just seen. Consider

where Π stands for the projection operation on U ad × R d + and ψ k is driven by an underlying process of i.i.d. drawings ξ k+1 , independent of {x i } i≤k . Let F k be the filtration generated by {x k , {ξ i } i≤k } so that ψ k and x k+1 are F k+1 measurable.

With the stochastic estimates produced by the AC and FD techniques considered so far in this paper, we obtained biased estimates of ∇ u Θ (and the bias sometimes also affects the estimate of Θ itself), with a bias going to 0 as k → +∞. We will denote Ψ(x k ) the correct value of the vector field at x k , namely

that with which an equilibrium point satisfies (see (31)):

for all positive ε. Define the martingale difference ∆M k , the bias B k and the variance V k of {ψ k } by:

We will use references [START_REF] Kushner | Stochastic approximation and recursive algorithms[END_REF] and [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF] in which convergence results and convergence rates of algorithm (32) are provided. Essentially, if the nonlinear projection operation at the r.h.s. of (32) is missing, under conditions on the quantities (35) in connection with the step size ε k that we will recall below, the trajectory produced by (32) behaves a.s. as that of the deterministic ODE:

In the presence of the projection onto a closed convex set, the differential equation is more complex to write since it involves another process z in which z takes values in the orthogonal cone C(x) to the convex set at the current point x (hence this process effectively appears only at the border of the convex set). The ODE now reads

The role of z is to maintain x in the convex set, as it is the case for x k produced by (32). It is defined as the "minimum force" which achieves this goal.

Convergence

In this subsection, we recall the conditions which ensure that the stochastic Arrow-Hurwicz algorithm will behave as its ODE (36) and we refer to the previous subsection to deduce that primal iterates u k will converge, at least locally, towards the solution u ♯ (assumed unique) of the constrained optimization problem. We then apply those results to the case of biased gradient estimates provided by AC and FD methods to derive a policy on how to tune the parameters r (see ( 19)) and c (see ( 26)) as functions of the iteration index k in order to satisfy the convergence conditions.

the interior of the convex set onto which Π is the projection), function Ψ admits the following representation:

where A is a matrix with eigenvalues µ satisfying

Then, a direct application of Theorem 3.1 in [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF] gives the asymptotic mean square error (MSE) as a function of the algorithm parameters γ, β, δ and it states that:

Some comments are in order here regarding the application of this result to our situation. First, the authors of [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF] state than when x ♯ lies on the boundary of the feasible convex set, other techniques (e.g. large deviations) are required to establish convergence rates. In our case, we expect that the probability constraint is active at the optimum, hence the optimal dual variable should be strictly positive. In the example of §3.2, we also have positivity constraints on primal variables and that on u (the first primal component) is active at the optimum (see ( 14)). But it is felt that the projection is rather helpful in accelerating convergence for this component (see numerical results in the next section). We may consider that, asymptotically, u k is "frozen" at 0 and does not participate to the dynamics of the algorithm ultimately.

Second, condition (44) may not be satisfied. We will come back on this point in the next subsection. Nevertheless, we used the results of [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF] as guidelines for the choice of parameters β and γ to drive the primal solution to its equilibrium in the most efficient way.

That said, in order to achieve the fastest convergence rate, one should seek to maximize κ in (45) over the feasible set defined by (39) or (41) and the expression of δ as a function of β. For the case of AC estimates, δ = -β/2, the minimum of 2β and γ -β/2 is obtained when those two functions are equal, which yields β = 2γ/5 and a value of 4γ/5; because of the first condition (39), the maximal possible value is obtained with γ = 1, which yields β = 2/5 and κ = 4/5, and we check that this pair (β, γ) all conditions in (39). Observe that our heuristic reasoning at the end of §4.1.1 and §4.1.2 in order to tune the parameter r when N i.i.d. samples are available (here N is the iteration index k) yields the same results (see (25) in particular).

For the case of FD estimates, under assumption (H3) of Theorem 6, the calculations and conclusions are the same. Under assumption (H4), δ = -3β/4 and the optimal values are β = 4/11, γ = 1, κ = 8/11 which is of course worse than the previous case. Finally, in the worst case for FD, we get β = 1/3, γ = 1, κ = 2/3.

The following result is a direct application of [START_REF] L'ecuyer | Budget dependent convergence rate of stochastic approximation[END_REF]Th. 4.1 and 4.2]. This CLT gives additional information on the asymptotic behavior of the iterates of (32).