Gilles Pagès
email: gilles.pages@upmc.fr.

applications for numerics

Keywords: Optimal vector quantization, greedy quantization, quantization tree, Lloyd's I algorithm, Competitive Learning Vector Quantization, stochastic gradient descent, learning algorithms, Zador's Theorem, Feynman-Kac's formula, variational inequality, optimal stopping, quasi-Monte

may come

Introduction to vector quantization and its applications

for numerics

Introduction to optimal vector quantization and its 1 Introduction to vector quantization 1.1 Signal transmission, information

The history of optimal vector quantization theory goes back to the 1950's in the Bell laboratories where researches were carried out to optimize signal transmission by appropriate discretization procedures. Two kinds of "stationary" signal can be naturally considered: either a deterministic -more or less periodic -signal (denoted by (x t ) t≥0 ) or a stochastic signal (denoted by (X t ) t≥0 ), considered under its stationary regime, supposed to be ergodic. In both cases, these signals share an averaging property as will be seen further on. Vector quantization can be briefly introduced as follows.

Let Γ = {x 1 , . . . , x N } be a subset of R d (d ≥ 1) of size (at most) N ≥ 1, called a quantization grid or simply a quantizer at level N if Γ has exactly cardinality N i.e. if the elementary quantizers x i are pairwise distinct. When d = 1 the numbering of the elementary quantizers x i is a priori made consistent with the natural order on the real line so that i → x i is non-decreasing.

In what follows, except specific mention, | . | will denote the canonical Euclidean norm on R d (although many of the stated results remain true or admit variants for more general norms).

A Γ-valued quantization function (also called quantizer) is simply any Borel function q : R d → Γ. A naive idea is to transmit at time t the stochastic signal q(X t ) instead of X t itself inducing a resulting pointwise error |X tq(X t )| One proceeds likewise for a deterministic signal with a resulting error |x tq(x t )|.

⊲ Deterministic signal : Let p ∈ (0, +∞). Assume that 1 t = ξq(ξ) L p (µ) < +∞.

⊲ Stationary ergodic stochastic signal : We consider again p ∈ (0, +∞). Assume the process (X t ) t≥0 is stationary. Then, X t has the same marginal distribution, say µ, for every t ∈ R + . Moreover, if E |X t | p = R d |ξ| p µ(dξ) < +∞, then X tq(X t ) L p (P) = X 0q(X 0 ) L p (P) = ξq(ξ) L p (µ) < +∞.

Moreover, if the process (X t ) t≥0 is ergodic, ergodic pointwise Birkhoff's Theorem ensures that P-a.s. lim t→+∞

1 t t 0 |X s -q(X s )| p ds 1 p = ξ -q(ξ) L p (µ) < +∞.
At this stage, several questions arise to optimize the transmission. Based on what precedes, we will mainly adopt from now on the static point of view of an R d -valued random vector X, defined on a probability space (Ω, A, P), with distribution µ. It corresponds to the value of X t at any time t or to the asymptotic behavior of the signal (x(t)) t≥0 . More general situations of quantization or coding can be investigated in Information Theory which take into account the dynamics of the (ergodic) process leading to the most general Shannon's source coding theorem. For these deeper aspects from Information Theory, we refer to the general distortion theory as analyzed by large deviation methods in [START_REF] Dembo | Large deviations techniques and applications[END_REF] and the references therein.

Question 1 How to optimally choose the Γ-valued quantization function q (Geometric optimization)?

It is clear that, whatever the quantization function q : R d → Γ is, one has

|ξ -q(ξ)| ≥ dist(ξ, Γ)
where dist(ξ, A) = inf a∈A |ξ -a| denotes the distance of ξ to the set A ⊂ R d (with respect to the current norm). One easily checks that equality holds in the above inequality if and only if q is a Borel 2 nearest neighbour projection i.e. q = π Γ defined by

= π Γ (ξ) = N i=1 x i 1 C i (Γ) (ξ)
where the N -tuple of subsets C i (Γ) 1≤i≤N is a Borel partition of (R d , Bor(R d )) satisfying

∀ i = 1, . . . , N, C i (Γ) ⊂ ξ ∈ R d : |ξ -x i | = min 1≤j≤N |ξ -x j | .
Such a partition of R d is called a Voronoi partition (or sometimes tessellation) induced by Γ. When the norm | . | is Euclidean the closures C i (Γ) of the cells are non-empty polyhedral closed convex sets (intersection of finitely many half-spaces defined by median hyperplanes of the couples of points (x i , x j ), i = j). One easily shows that

ξ ∈ R d : |ξ -x i | < min 1≤j≤N, j =i |ξ -x j | ⊂ • Ci (Γ) ⊂ C i (Γ) ⊂ ξ ∈ R d : |ξ -x i | = min 1≤j≤N |ξ -x j | .
The inclusions at both ends of the inclusion chain can be replaced by equalities in an Euclidean framework.

Then, for a given (static) random vector having values in R d , one defines a Voronoi Γ-quantization of X by Γ as

X Γ = π Γ (X).
Remark. For more developments on the non-Euclidean framework, like e.g. the ℓ r -norms defined by

|ξ| r = |ξ 1 | r + • • • + |ξ d | r 1 r , r ∈ [1, +∞), or |ξ| ∞ = max 1≤j≤d |ξ j |, ξ = (ξ 1 , . . . , ξ d ) ∈ R d
, we refer to [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF], Chapter 1. This leads us to define for p ∈ (0, +∞), the L p -mean quantization error induced by a grid Γ by e p (Γ, X) = Xπ Γ (X) L p (P) = dist(X, Γ) L p (P) = min Note that, from a computational point of view, the computation of π Γ (ξ) is very demanding when the size N is large since it amounts to a nearest neighbour search. We will come back to that point further on in Section 3 devoted to numerical aspects of (optimal) quantization grid computation.

It is clear that, if we define the optimal L p -mean quantization problem by e p,N (X) = inf Γ, card(Γ)≤N e p (Γ, X) (1.3) where card(Γ) denotes the cardinality of the grid Γ ⊂ R d , then e p,N (X) = inf

(x 1 ,...,x N )∈(R d ) N
G p,N (x 1 , . . . , x N )

1 p .
Note that, in fact, e p,N (X) only depends on the distribution µ of X. So we will occasionally write e p,N (µ) instead of e p,N (X). This follows from the easy remark that a grid Γ with less than N elements can be represented by always an N -tuple in which each element of the grid appears as a component at least once.

Proposition 1.1. Let p ∈ (0, +∞). Assume that X ∈ L p R d (P) i.e. The proof of this proposition is postponed to Section 2.1. The grid Γ (N,p) , the corresponding Ntuples x (N,p) (there are N ! N -tuples obtained by permutations of the components if the grid has full size N ) as well as the (Borel) nearest neighbour projections q Γ (N,p) are all called L p -optimal quantizers.

Of course a crucial question in view of possible applications is to compute such L p -optimal quantizers at level N , especially in higher dimension.

When d = 1 and µ = U ([0, 1]), then, for any p ∈ (0, +∞), the mid-point grid Γ (N,p) = 2i-1 2N , i = 1, . . . , N is the unique optimal L p -quantizer at level N . The attached weights are all equal to w (p,N ) i = 1

N , i = 1, . . . , N and the resulting optimal L p -quantization error is given for every N ≥ 1 by e p,N U ([0, 1]) = 1 2(1 + p) 1/p N .

(1.4)

More generally the question of the rate of decay of e p,N (X) is the central question of optimal vector quantization theory. It will be investigated further on in Section 2.3.

Application to signal transmission

As mentioned in the introduction, this application of (optimal) quantization goes back to the very origin of quantization theory in the 1950's. Imagine one has access to an L p -optimal quantization grid, say for p = 2 (quadratic case in an Euclidean setting). For convenience, we assume that Γ = {x 1 , . . . , x N } is a grid (possibly optimal) such that P X ∈ 1≤i≤N ∂C i (Γ) = µ 1≤i≤N ∂C i (Γ) = 0 e.g. because µ assigns no mass to hyperplanes.

What is the information "contained" in X Γ = π Γ (X)? Or equivalently, in probabilistic terms, what are the characteristics of the distribution of X Γ ? 1. Its state space Γ = {x 1 , . . . , x N }, 2. Its "companion" weights w i = w (Γ) i = P( X Γ = x i ) = P(X ∈ C i (Γ)) = µ(C i (Γ)), i = 1, . . . , N .

If X is a random vector with a known simulatable distribution µ, one can pre-compute these weights w i with an arbitrary accuracy by a large scale Monte Carlo simulation since, owing to the Strong Law of large Numbers,

w i = P-a.s. lim M →+∞ 1 ≤ m ≤ M : |X m -x i | < min j =i |X m -x j | M , i = 1, . . . , M,
where (X m ) m≥1 is a sequence of i.i.d. random vectors with distribution µ. In case of a not too large dataset (which is more common for signal transmission) one compute the weights w i by assigning its Voronoi class to each data (when the dataset is too large a statistical approach is adopted).

Coding the (quantized) signal Let Γ = {x 1 . . . , x N } ⊂ R d be a grid (possibly non-optimal at this stage) of size N and let P(Γ) be the set of distributions whose support is exactly Γ. In order to transmit a Γ-valued signal from a sender A to a receiver B, A will transmit a codeword C i = C(x i ) representative of x i instead of (an accurate enough approximation of) x i itself. For simplicity we will assume that the coding function C maps Γ into the set {0, 1} (N) of finite {0, 1}-valued sequences (this means that we adopt a dyadic coding procedures). The set {0, 1} is called a 2-alphabet (1 ). Our first request on the function C is identifiability i.e. that B can always recover x i from C i or equivalently that C is injective. To design the codewords (C i ) 1≤i≤N , one aims at minimizing the mean transmission cost κ, also known as the mean length of the message. This is in fact very old problem which goes back to the origins of Information Theory introduced by Claude Shannon in [START_REF] Shannon | Discrepancy and integration of continuous functions[END_REF].

Let us focus for a while on this coding problem. The mean transmission cost κ(N ) for a grid of size N is clearly defined by

κ(N ) = N i=1 w i × length(C i ).
A first (not so naive) idea is to re-index the points x i by a permutation σ so that i → w σ(i) is nonincreasing. Without loss of generality, we may assume from now on that σ is identity (though, for one-dimensional distributions, it is not consistent in general with the natural order of the points x i on the real line). Then, it is intuitive (but in fact not mandatory) to devise the coding function C so that i → length(C i ) is non-decreasing since, doing so, the more often a code is transmitted, the shorter it will be. In case of equality (like for the uniform distribution over Γ conventions have to be made).

A naive approach is to simply code x i through the regular dyadic expression ī2 of i which needs 1 + ⌊log 2 i⌋ digits (where ⌊ξ⌋ denotes the lower integer part of ξ ∈ R). This yields

κ(N ) = N i=1 w i 1 + ⌊log 2 i⌋ = 1 + N i=1 w i ⌊log 2 i⌋ ≤ 1 + ⌊log 2 N ⌋.
The transmission relies on the fact that both A and B share the codebook i.e. a one-to-one correspondence

x i ←→ ī2 .

(1.5)

⊲ A toy example. Imagine that, to transmit a uniformly distributed signal over the unit interval [0, 1], we first optimally quantize it using the mid-point grid Γ (N ) = 2i -1 2N , i = 1, . . . , N . This is equivalent to transmit a uniformly distributed signal over {1, . . . , N } thanks to the codebook so that, as far as transmission is concerned, the grid Γ (N ) itself plays no role. The resulting mean transmission cost κ(N ) is equal to

κ(N ) = 1 + 1 N N i=1 ⌊log 2 i⌋ ∼ log 2 N/e as N → +∞.
To be more precise, once noted that the dyadic entropy H 2 μUnif N of the uniform distribution μUnif N over {1, . . . , N } (or equivalently on Γ (N ) ) is equal to log 2 N , we can show that

c -= lim inf N κ(N ) -H 2 μUnif N ≤ lim sup N κ(N ) -H 2 μUnif N ) = c +
where c -≈ -2, 8792 and c + ≈ -0.9139.

⊲ Instantaneous coding. However, this approach is definitely too naive. In practice, A does not send one isolated codeword but a sequence of codewords. Such a coding is not satisfactory, mainly because it is not self-punctuated. To be decodable, an extra symbol (space, comma, etc) is needed to isolate the codewords. Doing so amounts to adding one symbol to the alphabet (with a special status since it cannot be repeated, like the large space in Morse coding). But this lowers the global performance of the coding system since it induces de facto switching from a 2-alphabet to a 3-alphabet coding function C, the third symbol having moreover a lower status of "under-symbol". To overcome this problem, the idea, again due to Shannon in his seminal 1948 paper [START_REF] Shannon | Discrepancy and integration of continuous functions[END_REF], is to devise self-punctuated codes. This relies on two conditions. First we ask the coding process to be uniquely decodable in the sense that the concatenation of codewords

C(x 1 ) • • • C(x N ) uniquely characterizes the concatenation x 1 • • • x N .
The additional condition which defines an instantaneous coding system is that a codeword can never be the prefix of another or, equivalently, no codeword can be obtained as the concatenation of another codeword and further symbols of the alphabet (here 0 and 1 digits). One easily checks that an instantaneous coding procedure is always self-punctuated.

Unfortunately, it is also straightforward to check that the naive dyadic coding (1.5) formerly mentioned which consists in writing in base 2 every index i is not an instantaneous coding system since, e.g., 22 = 10 and 52 = 101. Let us illustrate on a simple example how an instantaneous coding procedure look. We consider the following coding procedure of the set of indices {1, 2, 3, 4}:

C(1) = 0, C(2) = 10, C(3) = 110, C(4) = 111.
Such a code is uniquely decodable (e.g. 0110111100110 can be uniquely decoded as the string 134213). Furthermore it is clearly instantaneous (thus 010111110010 can be parsed on line as 0, 10, 111, 110, 10 i.e. the string 12432).

If we consider the uniform distribution μUnif 4 over {1, 2, 3, 4}, the resulting mean cost transmission (or mean length) is equal to κ μUnif

4 := 1 4 (1 + 2 + 3 + 3) = 9 4
whereas the naive dyadic coding of the indices seemingly yields 8 4 = 2. However, the implementable version of this naive dyadic coding (1.5), i.e. including an extra symbol like ",", has a mean length equal to 3 > 9 4 . This can be up to 30% more symbol consuming than the above instantaneous code! Now, let us consider a general distribution μN (exactly) supported by {1, . . . , N } (or equivalently by a grid Γ N of size N ) and a priori not uniform. Assume we have access to the weights w i = μN {i} (i.e. to μ itself). We define the dyadic entropy H 2 (μ) of μ by

H 2 (μ) = - N i=1 w i log 2 w i .
Then, the following classical theorem from Information Theory holds (see [START_REF] Cover | Elements of Information Theory[END_REF], Chapter 5, Theorem 5.3.1 and Section 5.4).

Theorem 1.1. For any instantaneous dyadic coding procedure C : {1, . . . , N } → {0, 1} (N) of the distribution μ, its mean transmission cost κ μ(N ) satisfies

κ(μ N ) ≥ H 2 (μ N ).
(1.6)

Furthermore, there exists (at least) one instantaneous coding procedure such that

κ μ(N ) < H 2 (μ N ) + 1. (1.7)
For a proof of this result based on Kraft's inequality, which is too far from the scope of this paper, we refer to [START_REF] Cover | Elements of Information Theory[END_REF]. Furthermore when a sequence (Y n ) n≥0 of {1, . . . , N }-valued signals to be transmitted is stationary with marginal (invariant) distribution μN and ergodic, it is possible by aggregating n of them to show (with obvious notations, see again [START_REF] Cover | Elements of Information Theory[END_REF]) that

κ(Y 1 , . . . , Y n ) → H 2 (μ) as n → +∞.
(1.8)

Examples: (a) The Huffman code: It was the first optimal instantaneous code -devised in Huffman's PhD thesis (see also [START_REF] Huffman | A method for the construction of minimum redundancy codes[END_REF]). Its length sequence (ℓ * i ) 1≤i≤N can be obtained as the solution to the integer optimization problem (ℓ i denotes the length of a code C i ):

ℓ * = argmin 2 -ℓ i ≤1 w i ℓ i so that H 2 (μ N ) ≤ κ Huf (μ N ) = w i ℓ * i ≤ H 2 (μ N )+1.
For an explicit construction of the Huffman code (and not only of its length sequence!), we refer again to [START_REF] Cover | Elements of Information Theory[END_REF], Sections 5.6 & 5. Let us simply mention that the codes are obtained by the concatenation of labels given to the edges (1 for "right" edges, 0 for "left edges" starting from the root) of successive trees built from the (non-increasing) monotony of the weights w i . The se successive trees are obtained by summing up the lower probabilities (starting from wN-1 := w N +w N -1 ), with appropriate conventions in case of equality like with uniform distributions).

(b) The Shannon coding (see exercise 5.28 in [START_REF] Cover | Elements of Information Theory[END_REF]): Still assume that the weights of the distribution μN satisfy 0 < w N ≤ • • • ≤ w 1 < 1. Let F μN denote the strict-cumulative distribution function of μN defined by N, where ⌈ξ⌉ denotes the upper integer part of the real number ξ. Elementary computations show that Shannon's code is instantaneous and that its mean transmission cost κ Shan S(μ N ) also satisfies

F μN i = j<i w j . Set ℓ i = ⌈-log 2 w i ⌉ and C i = ⌊2 ℓ i F μ i ⌋, i = 1, . . . ,
H 2 (μ N ) ≤ κ Shan S(μ N ) < H 2 (μ N ) + 1.
Global error induced by the transmission of a quantized signal Let us bring back quantization into the game by considering a continuous signal which needs to be quantized in order to reduce its transmission cost. Let us briefly compare from a quantitative viewpoint two modes of transmission for a signal.

⊲ Direct transmission. Let (X t ) t≥0 be a stochastic stationary signal with marginal distribution µ defined on a probability space (Ω, A, P) and Γ = {x 1 , . . . , x N }. To transmit the Γ-quantization X Γ of the random signal X = X t0 at time t 0 , the resulting quadratic mean quantization error is equal to X -X Γ L 2 (P) + 2 -r = e 2 (Γ, µ) + 2 -r where 2 -r is the dyadic transmission accuracy of any of the elementary quantizers x i . In fact this corresponds to a fixed transmission cost κ = r + 1 i.e. number of dyadic digits used to transmit these values. Common values for r lie between 10 and 20 (having in mind that 2 -10 = 1 1024 ≈ 10 -3 ). ⊲ Signal transmission using the codebook. If the receiver B uses the codebook (C i ←→ x i ) 1≤i≤N for the decoding phase (2 ), the resulting mean quadratic transmission error will be equal to X -X Γ L 2 (P) = e 2 (Γ, µ) whereas the mean unitary transmission cost is κ μ(N ) where μ is the distribution of the quantized signal X Γ . In this second case, there is a connection between the transmission error and the transmission cost that will be made more precise in Section 2.3 when the grid Γ is L 2 -optimal at level N for µ.

However, in the very simple case of the uniform distribution U ([0, 1]) over the unit interval, we can establish a direct relation between quadratic mean transition error and mean transmission cost κ when both the quantization and the instantaneous coding are optimal. The optimal quadratic quantization of U ([0, 1]) is the uniform distribution μUnif N over the N -mid-point whose dyadic entropy is exactly

H 2 (μ U nif N ) = log 2 N .
Plugging this equality in (1.7) yields κ μN ≤ log 2 (N ). In turn, plugging this inequality in the quantization error bound (1.4) yields that the lowest achievable mean transmission error, for a prescribed mean transmission cost κ, approximately satisfies

2 -(κ+1) √ 3 ≤ L 2 -Mean transmission error(κ) ≤ 2 -κ √ 3 .
A less sharp (reverse) formulation is

-log 2 Transmission error(κ) ∼ κ as κ → +∞.
This result appears as the most elementary version of a seminal result from Information Theory: Shannon's source coding theorem (in one dimension). Its extension to more general distributions µ on R d will be possible, once established the sharp convergence rate of the L 2 -optimal mean quantization error of general distributions on R d stated in Section 2.3 (Zador's Theorem). We enhanced in the above lines a static random signal presentation but the adaptation to a stationary process or a quasi-periodic (as defined above in terms of weak convergence of its time empirical measure) signal is straightforward (in particular for stationary ergodic signal one may take advantage of the improvement provide by (1.8), using n-aggregates of the signal, to reduce the range of the two-sided inequality (1.6)-(1.7) in Theorem 1.1. where δ a denotes the Dirac mass at a ∈ R d . In such a framework, n is usually large, say 10 6 or more, and optimal quantization can be viewed as a model for clustering i.e. the design of a set of N prototypes of the dataset, with N ≪ n, obtained as a solution to the mean quadratic (or more generally L p -) optimal quantization at level N ≥ 1 of the distribution µ (p ∈ (0, +∞) being fixed). This reads as the (L p -)minimization problem min

What else is quantization for?

(x 1 ,...,x N )∈(R d ) N 1 n n k=1 min 1≤i≤N |ξ k -x i | p .
The existence of such an optimal N -quantization grid Γ (N,p) of prototypes follows from the above Proposition 1.1. Such a distribution does assign mass to hyperplanes and in particular to the boundaries of polyhedral Voronoi cells. However, owing to Theorems 4.1 and 4.2 in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF] (p.38), we know that the boundaries of the Voronoi cells induced by an optimal grid Γ (N,p) are always µ-negligible.

Once an optimized grid of N prototypes has been computed (see Section 3 devoted to the algorithmic aspects), it can be used to produce an automatic classification of the dataset by making up "clusters" of points of the dataset following the nearest neighbour rule among the prototypes. Formulated equivalently, one defines the N clusters as the "trace" of the dataset on the Voronoi cells

C i (Γ (N ) ), i = 1, . . . , N .
From a mathematical point of view, investigations on this topic are carried out be replacing the deterministic dataset (ξ k ) 1≤k≤n by a sequence of i.i.d. random vectors (X k ) k≥0 defined on a probability space (Ω, A, P) with distribution µ. The quantities of interest become, in short, the sequence of optimization problems induced by the random empirical measures µ n (ω, dξ) = 1 n n k=1 δ X k (ω) (dξ), ω ∈ Ω. This has given rise to a huge literature in Statistics and has known a kind of renewal with the emergence of clustering methods in the "Big Data" world, see [START_REF] Biau | On the performance of clustering in Hilbert spaces[END_REF]. We consider, for every ω ∈ Ω, the optimization problem min

(R d ) N 1 n n k=1 min 1≤i≤N |X k (ω) -x i | p = R d min 1≤i≤N |ξ -x i | p µ n (ω, dξ) .
(1.9)

The main connection with optimal quantization is the following: Assume that µ(B(0; 1)) = 1. For every ω ∈ Ω, there exists (at least) an optimal N -tuple x (N ) (ω, n) for the above problem which satisfies

E e 2 x (N ) (ω, n), µ -e 2,N (µ) ≤ C min   N d n , d N 1-2 d log n n  
where C > 0 is a positive universal real constant.

1.3.2 From Numerical integration (I) . . .

Another way to take advantage of optimal quantization emerged in the 1990's (see [START_REF] Pagès | Voronoi Tessellation, space quantization algorithms and numerical integration[END_REF]). As we know that for a sequence (Γ (N,p) ) N ≥1 of L p -optimal grids of size N with N → +∞ we have X -X Γ N,p L p (P) = e p,N (X) → 0 i.e. X Γ N,p → X as N → +∞ in L p (hence in distribution). It can be shown (see [START_REF] Delattre | Quantization of probability distributions under norm-based distortion measures[END_REF]) that, in fact, this convergence also holds in an a.s. sense although we will make little use of this feature in what follows.

In particular, if a function F : R d → R is bounded and continuous, then E F ( X Γ N,p ) → E F (X) as N → +∞. On the other hand, using the characteristics (x

(N ) i , w (N ) i
) 1≤i≤N of the distribution of X Γ N,p , we derive a very simple weighted cubature formula

E F ( X Γ N,p ) = N i=1 w (N ) i F x (N ) i .
(1.10)

When F has more regularity (and possibly not bounded), we can easily establish precise error bound for this quantization based cubature formula.

First order quantization based cubature formula. More generally, we consider a function F : R d → R which is locally α-Hölder continuous in the sense that there exists α ∈ (0, 1], β ≥ 0, and a real constant

[F ] α,β such that ∀ x, y ∈ R d , |F (x) -F (y)| ≤ [F ] α,β |x -y| α 1 + |x| β + |y| β .
Then, for every conjugate Hölder exponents (p, q) ∈ [1, +∞],

E F (X) -E F ( X Γ ) ≤ [F ] α,β E |X -X Γ | α 1 + |X| β + | X Γ | β ) ≤ [F ] α,β,loc X -X Γ α L αp (P) 1 + X β L βq (P) + X Γ β L βq (P) .
In particular, if p = 1 α , one gets

E F (X) -E F ( X Γ ) ≤ [F ] α,β X -X Γ α 1 1 + X β L β 1-α (P) + X Γ β L β 1-α (P) (1.11)
(with the convention . 0

L β 1-1 (P) = 1). If F is α-Hölder continuous with Lipschitz coefficient [F ] α = 1 3 [F ] 1,0 , then E F (X) -E F ( X Γ ) ≤ F ] α X -X Γ α L α (P) .
(1.12)

From the cubature formula (1.12) and using that bounded Hölder functions characterize the weak convergence of probability measures, we derive the following corollary about L p -optimal quantizers (by considering α = p ∧ 1).

Corollary 1.1. Let (Γ (N ) ) N ≥1 be a sequence of (L p , µ)-quantizers with Γ (N ) of size N and such that e p (Γ (N ) , µ) → 0 as N → +∞. Let μN denote the distribution of X Γ (N ) . Then

μN = N i=1 µ C i (Γ (N ) ) δ x (N ) i (w) -→ µ as N → +∞.
(1.13) where (w)

-→ denotes the weak convergence of distributions.

. . . to Numerical Probability (conditional expectation)

One of the main problem investigated in the past twenty years in Numerical Probability has been the numerical computation of conditional expectations, mostly motivated by problems arising in finance for the pricing of derivative products of American style or more generally known as "callable". It is also a challenging problem for the implementation of numerical schemes for Backward Stochastic Differential Equations (see [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF][START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF]), Stochastic PDEs (see [START_REF] Gobet | Discretization and simulation of the Zakai Equation[END_REF]), for non-linear filtering [START_REF] Pagès | Optimal quantization methods for nonlinear filtering with discrete-time observations[END_REF][START_REF] Pham | Approximation by quantization of the filter process and applications to optimal stopping problems under partial observation[END_REF] or Stochastic Control Problems (see [START_REF] Pagès | An Optimal Markovian Quantization Algorithm for Multidimensional Stochastic Control Problems[END_REF][START_REF] Brandejsky | Optimal stopping for partially observed piecewise-deterministic Markov processes[END_REF][START_REF] Brandejsky | Numerical method for impulse control of piecewise deterministic Markov processes[END_REF]). Further references are valuable in the survey paper [START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF] devoted to applications of optimal vector quantization to Numerical Probability. The specificity of these problems in the probabilistic world is that, whatever the selected method is to suffer in some way or another, from the curse of dimensionality. Optimal quantization trees (introduced in [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF]) is one of the numerical methods designed to cope with this problem (with regression and Monte Carlo-Malliavin method, see [START_REF] Longstaff | Valuing American options by simulation: a simple least-squares approach[END_REF], [START_REF] Fournié | Some applications of Malliavin calculus to Monte Carlo methods in finance[END_REF]). The precise connection between vector quantization and conditional expectation computation can be summed up in the proposition below. We consider a couple of random vectors (X, Y ) : (Ω, A, P) → R d × R q and the regular version Q of the conditional distribution operator of X given Y , defined on every bounded or non-negative Borel function f : R d → R, by

Qf (y) = E f (X) | Y = y .
Then, Qf is a Borel function on R d . We define the Lipschitz ratio of a function f :

R d → R by [f ] Lip = sup x =y |f (y)-f (x)| |x-y|
≤ +∞. We make the following Lipschitz continuity propagation assumption on Q: there exists

[Q] Lip ∈ R + such that ∀ f : R d → R, Borel function, [Qf ] Lip ≤ [Q] Lip [f ] Lip .
(1.14) Proposition 1.2. Assume that the conditional distribution operator Q of X given Y satisfies the above Lipschitz continuity propagation property (1.14). Let Γ X ⊂ R d and Γ Y ⊂ R q be two quantization grids of X and Y respectively.

(a) Quadratic case. Assume X, Y ∈ L 2 (P). Let f : R d → R be a Lipschitz continuous function and let g : R d → R be a Borel function with linear growth. Then

E f (X) | Y -E g( X Γ X ) | Y Γ Y 2 L 2 (P) ≤ [Qf ] 2 Lip Y -Y Γ Y 2 L 2 R q (P) + f (X) -g( X Γ X ) 2 L 2 (P) so that if g = f , E f (X) | Y -E f ( X Γ X ) | Y Γ Y 2 L 2 R q (P) ≤ [Qf ] 2 Lip Y -Y Γ Y 2 L 2 (P) + [f ] 2 Lip X -X Γ X 2 L 2 (P) .
(b) L p -case. Assume X, Y ∈ L p (P), p ∈ [1, +∞) and let f and g be like in (a). Then

E f (X) | Y -E g( X Γ X ) | Y Γ Y L p R q (P) ≤ (2 -δ p,2 )[Qf ] Lip Y -Y Γ Y L p (P) + f (X) -g( X Γ X ) L p (P)
where δ p,p ′ denotes the Kronecker symbol. In particular, if g = f , one has

E f (X) | Y -E f ( X Γ X ) | Y Γ Y L p (P) ≤ (2 -δ p,2 )[Qf ] Lip Y -Y Γ Y L p R q (P) + [f ] Lip X -X Γ X L p (P) . Proof. (a) We decompose E f (X) | Y -E f ( X Γ X ) | Y Γ Y into two (L 2 (P)-orthogonal) terms E f (X) | Y -E f ( X Γ X ) | Y Γ Y = E f (X) | Y -E E(f (X) | Y | Y Γ Y (1) + E E(f (X) | Y | Y Γ Y -E g( X Γ X ) | Y Γ Y (2)
.

To check the announced L ( P)-orthogonality, we note that (2) is σ( Y Γ Y )-measurable; hence, the characterization of conditional expectation given Y Γ Y impliesE (1) × (2) = 0. On the other hand, the very definition of conditional expectation given Y Γ Y as the best approximation in L 2 R q (P) by a square integrable σ( Y Γ Y )-measurable random vector implies in turn

E (1) 2 = E Qf (Y ) -E(Qf (Y ) | Y Γ Y ) 2 ≤ E Qf (Y ) -Qf ( Y Γ Y ) 2 ≤ [Qf ] 2 Lip Y -Y Γ Y 2 L 2 (P) .
On the other hand, using that Y Γ Y is σ(Y )-measurable, we first derive from the chain rule for conditional expectation that

(2) = E f (X) | Y Γ Y -E g( X Γ X ) | Y Γ Y = E f (X) -g( X Γ X ) | Y Γ Y .
Using now that conditional expectation is an L 2 -contraction, we deduce that

E (2) 2 ≤ f (X) -g( X Γ X ) 2 L 2 (P) ≤ f (X) -g( X Γ X ) 2 L 2 R q (P) . When g = f , the conclusion is straightforward. (b) We start from the classical Minkowski Inequality E f (X) | Y -E g( X Γ X ) | Y Γ Y L p (P) ≤ Qf (Y ) -E(Qf (Y ) | Y Γ Y ) L p (P) + E f (X) | Y Γ Y -E g( X Γ X ) | Y Γ Y L p (P)
where we used like in (a

) that E Qf (Y ) | Y Γ Y = E f (Y ) | Y Γ Y . Now, still owing to Minkowski's Inequality Qf (Y ) -E(Qf (Y ) | Y Γ Y ) L p (P) ≤ Qf (Y ) -Qf ( Y Γ Y ) L p (P) + E Qf ( Y Γ Y ) -Qf (Y ) | Y Γ Y L p (P) so that Qf (Y ) -E(Qf (Y ) | Y Γ Y ) L p (P) ≤ 2 Qf (Y ) -Qf ( Y Γ Y ) L p (P) ≤ 2[Qf ] Lip Y -Y Γ Y L p
R q (P) . Note that when p = 2 the above coefficient 2 can be cancelled using again, like in (a), that conditional expectation given Y Γ Y is the best approximator in L 2 (P) by σ( Y Γ Y )-measurable square integrable random vectors. On the other hand,

E f (X) | Y Γ Y -E g( X Γ X ) | Y Γ Y L p (P) ≤ f (X) -g( X Γ X ) L p (P) .
The case g = f follows immediately. This completes the proof.

To conclude this section, we make the connection between these cubature formulas and the L 1 -Wasserstein distance W 1 defined by Proof. (a) The inequality sup

W 1 (µ, ν) = inf E P |X -Y |, X, Y : (Ω, A, P) → R d , X d = µ, Y
[F ]p≤1
EF (X) -EF X Γ ≤ X -X Γ p L p (P) is straightforward (set e.g. α = p and β = 0 in (1.11) and note [F ] p = 1 3 [F ] p,0 ). The equality follows by noting that the function

F p defined for every ξ ∈ R d by F p (ξ) = min 1≤i≤N |ξ -x i | p is p-Hölder with [F ] p = 1. (b) Let X : (Ω, A, P) → R d with distribution P X = µ. It is clear, as already seen, that if Y : (Ω, A, P) → R d is such that Γ Y = Y (Ω) has at most N values then |X -Y | ≥ dist(X, Γ Y ) = |X -X Γ Y | so that X -X Γ Y 1 ≤ E|X -Y |.
As a consequence e 1,N (µ) ≤ W 1 (µ, P N ). Conversely, it follows from the definition of e 1,N (µ) in (1.3) that e 1,N (µ) ≥ W 1 (µ, P N ) since it is defined as an infimum over less random vectors (only those of the form Y = q(X) of X where q : R d → R takes at most N values) that the L 1 -Wasserstein distance. This completes the proof.

Application to Numerical Analysis (Automatic meshing)

1.4.1 Representation of the solution of parabolic equation, Feynman-Kac's formula Let b : [0, T ] × R d → R d and a : [0, T ] × R d → S + (d, R
) be two continuous functions with at most linear and quadratic growth in x, uniformly with respect to t ∈ [0, T ], respectively (S + (d, R) denotes the set of d × d symmetric non-negative matrices). Let f : R d → R be a Borel function with polynomial growth. We want to solve numerically the following parabolic partial differential equation (P DE), either by a Monte Carlo simulation or by a quadrature formula

∂u ∂t + Lu = 0, u(T, .) = f (1.15)
where, denoting by (.|.) the canonical inner product on R d ,

Lu = (b|∇u) + 1 2
Tr(a∇ 2 u).

(1.16)

⊲

Step 1 (The Feynman-Kac representation formula). This fundamental connection between diffusion process and (parabolic) PDEs is summed up in the following theorem. (a) Then the function u admits the following representation as an expectation:

∀ x ∈ R d , ∀ t ∈ [0, T ], u(t, x) = E f (X t,x T )
where (X x,t s ) s∈[t,T ] denotes the unique solution to the Stochastic Differential Equation (SDE)

dX t,x s = b(s, X t,x s )ds + σ(s, X t,x s )dW s , X t,x t = x, s ∈ [t, T ], (1.17)
(starting from x ∈ R d at time t ∈ [0, T ] and defined on [t, T ]) where W is a q-dimensional standard Brownian motion defined on a probability space (Ω, A, P).

Owing to the Markov property, an alternative formulation is given by

∀ t ∈ [0, T ], E f (X T ) | X t = u(t, X t ) a.s.
for any solution (X t ) t∈[0,T ] of the above SDE defined over the whole interval [0, T ] starting at a finite random vector X 0 independent of W . In particular

u(t, x) = E f (X T ) | X t =
x (in the sense that it is a regular version of the conditional expectation as x varies).

(b) Time homogeneous diffusion coefficients: If b(t, x) = b(x) and σ(t, x) = σ(x) (no dependence of b and σ in t), then the representation can be written

∀ x ∈ R d , ∀ t ∈ [0, T ], u(t, x) = E f (X 0,x T -t ). (1.18)
Proof. (a) Itô's formula applied to the function u and the process (s, X t,x s ) s∈[t,T ] between t and T yields

u(T, X t,x T ) = u(t, x) + T t ∂u ∂t + Lu (s, X t,x s ) =0 ds + T t ∇ x u(s, X t,x s )|σ(s, X t,x s )dW s .
The integral in "ds" is zero since u satisfies the parabolic PDE (1.15) and one easily establishes that the local martingale null at 0 defined by the Brownian stochastic integral is a true martingale, null at 0, owing to the growth control assumption made on ∇ x u. Then, one gets E u(T, X t,x T ) = u(t, x).

(b) One writes Itô's formula between 0 and Tt to u(Tt, X 0,x t ) and proceeds as above.

Remark. In the time homogeneous case, one can proceed by verification. Under smoothness assumption on b and σ, say C 2 with bounded existing derivatives and Hölder second order partial derivatives, one shows, using the tangent process of the diffusion, that the function u(t, x) defined by (1.18) is C 1,2 in (t, x). Then, the above claim (b) shows the existence of a solution to the parabolic P DE (1.15).

⊲

Step 2a (Monte Carlo simulation). Assume for the sake of simplicity that we want to compute a numerical approximation of u(0, x) = E f (X 0,x T ) i.e. that t = 0. At this stage, the idea is to replace the diffusion by its Euler scheme with step T n (n ≥ 1) starting at x: let t n k = kT n , k = 0, . . . , n be the uniform mesh of [0, T ] with step T n . It is recursively defined as follows (to alleviate notations, we drop the dependance in 0,x of the Euler scheme):

Xn t n k+1 = Xn t n k + T n b(t n k , Xn t n k ) + T n σ(t n k , Xn t n k )U (n) k+1 , k = 0, . . . , n Xn t n 0 = Xn 0 = x where (U (n) k ) k=1,.
..,n is an i.i.d. sequence of N (0; I q )-distributed random vectors representative of the Brownian increments i.e.

W t n k -W t n k-1 = T n U (n) k , k = 1, . . . , n.
The quantity Ef ( Xn T ) (have in mind that T = t n n ) is the counterpart of E f (X T ) for the Euler scheme. Assume that b and σ are Lipschitz continuous in (t, x). Then, on the one hand,

sup n≥1 max 0≤k≤n | Xn t n k | L p (P) + max t∈[0,T ] |X t | L p (P) ≤ κ p,b,σ,T 1 + |x| (1.19)
and, on the other hand, the discrete time Euler scheme strongly converges to X for the sup norm in every L p (P) at rate 1 n in the following sense max k=0,...,n

| Xn t n k -X t n k L p (P) ≤ C p,b,σ,T T n 1 + |x| . As a consequence, E f ( Xn T ) → E f (X T ) with a O 1
n -rate as the step T n goes to 0 (i.e. n → +∞) if f is Lipschitz continuous (or even locally Lipschitz continuous with polynomial growth). The convergence still holds, without rate, if f is continuous with polynomial growth.

In fact if b, σ and f are smooth enough then the so-called weak error E f ( Xn T ) -E f (X T ) can be investigated directly by more analytic methods. As a result, a (faster) O 1 n -ate can be established (see [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF]). This rate can be extended to bounded Borel functions f provided σ satisfies a uniform ellipticity property (or even a hypo-ellipticity assumption "à la Hörmander" for a modified Euler scheme) as proved in a celebrated Bally-Talay's paper (see [START_REF] Bally | The distribution of the Euler scheme for stochastic differential equations: I. Convergence rate of the distribution function[END_REF]). This yields

u(0, x) = E f (X 0,x T ) = E f ( Xn T ) + O 1 n .
The point of interest at this stage is of course that the expectation E f ( Xn T ) can be computed by simulation since the Euler scheme can be straightforwardly simulated as soon as b and σ are computable functions (and X 0 itself can be simulated). So, we can implement a Monte Carlo simulation to compute Ef ( Xn T ) i.e. simulate M i.i.d. copies ( Xn T ) m m=1,...,M of the above Euler scheme at time T = t n n and approximate E f ( Xn T ) by the strong Law of Large Numbers

E f ( Xn T ) ≈ 1 M M m=1 f ( Xn T ) m
(simply because a.s. convergence holds as M → +∞). This second error (known as the Monte Carlo or the statistical error) is of order O( 1 √ M ) owing to the Central Limit Theorem which provides (asymptotic) confidence intervals for an a priori prescribed given confidence level involving the asymptotic variance Var f ( Xn

T ) = E f ( Xn T ) -E f ( Xn T ) 2 = Ef ( Xn T ) 2 -E f ( Xn T ) 2 .
In turn, this quantity can be expressed by expectations of functions of Xn T , consequently it can be computed on line as a companion parameter of the Monte Carlo simulation. (By the way note that one often has Var f ( Xn T ) ≈ Var f (X T )), either because f is continuous or because the diffusion is "elliptic enough", see above). For more details on these elementary aspects of the Monte Carlo method, we refer to classical textbooks devoted Monte Carlo simulation and Numerical Probability (see [START_REF] Lapeyre | Introduction to Monte-Carlo Methods for Transport and Diffusion Equations[END_REF] for a more PDE oriented introduction to Monte Carlo method or [START_REF] Glasserman | Monte Carlo Methods in Financial Engineering[END_REF][START_REF] Pagès | Introduction to Numerical Probability and Applications to Finance[END_REF] for more connections with Finance, among many others).

The main asset of this approach is that it is dimension free, in the sense that its complexity grows more or less linearly with the dimension d of the random vector of interest and with little influence of the ellipticity of the function a, at least when the function f is regular as we saw before.

⊲ Step 2b (Quantization based cubature formula). If one has many computations to carry out with the same operator L, i.e. for various functions f , it may be interesting to replace the Monte Carlo simulation by a cubature formula based on an optimal quantization of Xn T . To perform this quantization, as it will be seen further on in Section 3, one can rely on a stochastic optimization procedure which can be viewed as a kind of compressed Monte Carlo simulation. In that perspective, one faces now the following chain of approximations

u(0, x) = E f (X 0,x T ) ≈ E f ( Xn T ) ≈ E f Xn T Γ (N )
.

where Γ (N ) is an optimal (quadratic) quantization grid for the random vector Xn T . ⊲ Toward automatic meshing An alternative to the direct quantization procedure is to consider the grid Γ (N ) as a starting point to produce an optimized mesh for the numerical solving of the original P DE by deterministic schemes like finite element or finite volumes methods, etc. In such an approach, an optimal grid needs to be produced at each discretization time t n k . This approach has been widely investigated by Gunzberger's group in Florida (USA) (see e.g. [START_REF] Du | Centroidal Voronoi Tessellations: Applications and Algorithms[END_REF] and the references therein). More recently, a new concept of quantization (dual quantization, see [START_REF] Pagès | Intrinsic stationarity for vector quantization: Foundation of dual quantization[END_REF]) has refined this point for view by switching from Voronoi diagram to a direct approach of Delaunay triangulation which provides grids better fitted to deterministic numerical analysis methods in medium dimensions.

From optimal stopping theory to variational inequalities

Discrete time optimal stopping theory in a Markov framework. We consider a standard discrete time Markovian framework: let (X k ) 0≤k≤n be an R d -valued (F k ) 0≤k≤n -Markov chain defined on the filtered probability space (Ω, A, (F k ) 0≤k≤n , P) with transitions

P k (x, dy) = P X k+1 ∈ dy | X k = x , k = 0, . . . , n -1,
so that for every bounded or non-negative Borel function f :

R d → R, P k f (x) = R d f (y)P k (x, dy) and E f (X k+1 ) | F k = E f (X k+1 ) | X k = P k (f )(X k ).
From now on, we denote by F the filtration (F k ) 0≤k≤n . Intuitively, F k is a σ-field of A which represents the observable (or available) information at time k. Let Z = (Z k ) 0≤k≤n be an F-adapted obstacle/payoff sequence of non-negative integrable random variables of the form

0 ≤ Z k = f k (X k ) ∈ L 1 (Ω, F k , P), k = 0, . . . , n.
In term of modeling, this can be understood as follows: an agent plays a stochastic game. Each round of the game takes place at time k ∈ {0, . . . , n}. The random variable Z k represents the reward when leaving the game at time k. The question: "Is there an optimal way to quit the game in order to maximize the gain?"

By "quitting the game" we mean leaving possibly at a random time τ : Ω → {0, . . . , n} but always honestly i.e. in such a way that, for every ℓ ∈ {0, . . . , n}, the event

{τ = ℓ} = ω ∈ Ω | τ (ω) = ℓ ∈ F ℓ .
This means that if the agent adopts this strategy τ , all the information that leads him/her to leave the game at a time ℓ is available at time ℓ = τ (ω) after having observed the history of the game up to time ℓ i.e. (X k (ω)) 0≤k≤τ (ω) . Such a strategy is reasonable provided that the payoff sequence (Z ℓ ) ℓ=0,...,n has been supposed adapted to the filtration F as we did. Such a random variable is called an F-stopping time.

Imagine that the agent enters the game at time k ∈ {0, . . . , n}. The aim of the agent is to attain the optimal possible mean gain given the available information at time k, namely

U k = P-esssup E f τ (X τ ) | F k , τ : (Ω, A) → {k, . . . , n}, F-stopping time (1.20)
with an optimal mean gain given by E U k . The next question is to know whether there is an optimal stopping time, when starting the game at time k, i.e. a {k, . . . , n}-valued F-stopping time τ k satisfying

U k = E Z τ k | F k .
For more details on this topic we refer to [START_REF] Neveu | Martingales à temps discret[END_REF] or [START_REF] Lamberton | Introduction to stochastic calculus applied to finance[END_REF] (Chapter 2) or, more recently, [START_REF] Lamberton | Optimal stopping and American options[END_REF].

The sequence U = (U k ) 0≤k≤n is known as the (P, F)-Snell envelope of the sequence (Z k ) 0≤k≤n .

From a numerical point of view, we want to compute, or at least approximate, this Snell envelope, especially at time 0 and the related optimal stopping time τ 0 (if any).

The first important result of discrete time optimal stopping theory is the following Backward Dynamic Programming Principle (BDP P ).

Proposition 1.4. (a) The (P, F)-Snell envelope (U k ) 0≤k≤n satisfies the following BDDP :

U n = Z n and U k = max Z k , E U k+1 | F k , k = 0, . . . , n -1, (1.21)
and τ k = min ℓ ∈ {k, . . . , n} | U ℓ = Z ℓ is an optimal stopping time at time k i.e.

U k = E Z τ k | F k .
(b) Furthermore, for every k ∈ {0, . . . , n}, there exists a Borel function u k : R d → R such that

U k = u k (X k ), k = 0, . . . , n,
where the sequence (u k ) 0≤k≤n satisfies

u n = f n and u k = max f k , P k u k+1 , k = 0, . . . , n -1.
Proof. (a) We prove this claim by a backward induction on k. The fact that U n = Z n is obvious since τ n = n is the only {n}-valued stopping time (hence optimal at time n).

Now let k ∈ {0, . . . , n -1}. Assume that τ k+1 = min ℓ ∈ {k + 1, . . . , n} | U ℓ = Z ℓ is an optimal F-stopping time at time k + 1 i.e. U k+1 = E Z τ k+1 | F k+1 . It follows that E U k+1 | F k = E E Z τ k+1 | F k+1 | F k = E Z τ k+1 | F k ≤ U k since τ k+1 ≥ k + 1 is
U k ≥ max Z k , E U k+1 | F k .
To prove the reverse inequality and establish the BDP P at time k, we consider a {k, . . . , n}-valued F-stopping time τ . Then, using that {τ ≥ k

+ 1} = c {τ ≤ k} ∈ F k , E Z τ | F k = Z k 1 {τ =k} + E Z τ ∨(k+1) | F k 1 {τ ≥k+1} P-a.s. Now, using that τ ∨ (k + 1) is a {k + 1, . . . , n}-valued F-stopping time, E Z τ ∨(k+1) | F k = E E Z τ ∨(k+1) | F k+1 | F k P-a.s. ≤ E U k+1 | F k P-a.s.
by the definition (1.20) of U k+1 . As a consequence,

E Z τ | F k ≤ Z k 1 {τ =k} + E U k+1 | F k 1 {τ ≥k+1} ≤ max Z k , E U k+1 | F k .
Taking the P-esssup over all such {k, . . . , n}-valued F-stopping times, we get

U k ≤ max Z k , E U k+1 | F k P-a.s.
which in turn implies that the Snell envelope satisfies (1.21) at time k. Let us deal now with the optimal stopping time. One checks from its definition that

τ k = k1 {U k =Z k } + τ k+1 1 {U k =Z k } .
Using that both events

{U k = Z k } and {U k = Z k } lie in F k and that {U k = Z k } ⊂ E(U k+1 |F k ) = U k ,
we get the following string of equalities

E Z τ k | F k = Z k 1 {U k =Z k } + E Z τ k+1 | F k 1 {U k =Z k } P-a.s. = U k 1 {U k =Z k } + E E Z τ k+1 | F k+1 | F k 1 {U k =Z k } P-a.s. = U k 1 {U k =Z k } + E U k+1 | F k 1 {U k =Z k } P-a.s. = U k 1 {U k =Z k } + U k 1 {U k =Z k } P-a.s. = U k .
(c) This straightforwardly follows from the fact that, owing to the definition of the Markov transitions, if

U k+1 = u k+1 (X k+1 ) for k ∈ {0, . . . , n -1}, E U k+1 | F k = E U k+1 | X k = P k u k+1 (X k ) P-a.s. so that U k = max f k (X k ), E U k+1 | F k = max f k (X k ), P k u k+1 (X k ) = u k (X k ).
Remarks. • Note that claims (a) and (b) do not make use of the Markov property so that they remain true for any F-adapted integrable sequence Z = (Z k ) 0≤k≤n of non-negative random variables.

• The above optimal stopping times may be not unique but it can be shown that τ k as defined above is always the lowest stopping time for the game starting at time k (in the sense that another optimal stopping time τk for the game starting at k satisfies τk ≥ τ k a.s.). Moreover it follows from the above proof that the sequence of optimal stopping times (τ k ) 0≤k≤n satisfies the dual backward dynamic programing principle

τ k = k1 {U k =Z k } + τ k+1 1 {U k =Z k } . (1.22)
This second backward dynamic programming principle -sometimes called dual -is often used in regression methods to compute the Snell envelope (see e.g. Longstaff-Schwarz's paper [START_REF] Longstaff | Valuing American options by simulation: a simple least-squares approach[END_REF]).

Approximation of the Snell envelope by a quantization tree. The starting idea of the quantization tree method originally introduced in [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF] is to approximate the whole Markovian dynamics of the chain X = (X k ) 0≤k≤n using a sequence of quantizations (X k ) 0≤k≤n to produce a skeleton of the distribution of X supported by the tree made up by the quantization grids Γ k = {x k 1 , . . . , x k N k } of X k , k = 0, . . . , n and the transitions weights w k ij between states x k i and x k+1 j defined for every k ∈ {0, . . . , n -1} by

w k ij = P X k+1 = x k+1 j | X k = x k i , 1 ≤ i ≤ N k , 1 ≤ j ≤ N k+1 .
Although we will rely on these transitions below, it is important to have in mind that the sequence of quantizations ( X k ) 0≤k≤n is not a Markov chain in general.

At this stage, the idea is to mimic the BDP P (1.21) satisfied by the Snell envelope (U k ) 0≤k≤n by replacing X k by its Γ k -quantization X k . For every k ∈ {0, . . . , n}, we replace X k by a function

X k = π k (X k ) of X k taking values in a grid Γ k , (where π k : R d → Γ k is a Borel function).
We will see in the next section an optimal way to specify the function π k (including its value "grid"

Γ k = π k (R d )) by minimizing the induced L p -mean quadratic error X k -X k L p (P) .
The key point, keeping in mind as it has already been emphasized that the sequence ( X k ) 0≤k≤n has no reason to share a Markov property, is to force this Markov property in the BDP P . Doing so, we introduce a Quantized Backward (pseudo-)Dynamic Programming Principle for the sequence of quantized payoff/obstacle (f k ( X k )) 0≤k≤n reading

(QBDP P ) ≡ U n = f n ( X n ), U k = max f k ( X k ), E U k+1 | X k . (1.23)
The forcing of the Markov property is obtained by directly conditioning by the single random vector X k rather than by the σ-field

F k := σ( X ℓ , 0 ≤ ℓ ≤ k).
Then it is straightforward by a backward induction that, for every k ∈ {0, . . . , n},

U k = u k ( X k ), u k : R d → R + , Borel function.
Then, from a computational point of view, (1.23) reads

u n (x n i ) = f n (x n i ), 1 ≤ i ≤ N n , u k (x k i ) = max f k (x k i ), N k+1 j=1 w k ij u k+1 (x k+1 j ) , 1 ≤ i ≤ N k , 1 ≤ j ≤ N k+1 , 0 ≤ k ≤ n -1, (1.24)
where

Γ k = {x k 1 , . . . , x k N k }, k = 0, . . . , n. See subsection 2.
3.4 for details on the practical implementation, including the computation by Monte Carlo simulation of the transition weights w k ij .

Error bounds. The following theorem establishes the control on the approximation of the true Snell envelope (U k ) 0≤k≤n by the quantized pseudo-Snell envelope ( U k ) 0≤k≤n using the L p -mean approximation errors X k -X k L p (P) .

Theorem 1.3 (see [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF] (2001), [START_REF] Pagès | Intrinsic stationarity for vector quantization: Foundation of dual quantization[END_REF] (2011)). Assume that all functions f k : R d → R + are Lipschitz continuous and that the transitions P k (x, dy) = P(X k+1 ∈ dy|X k = x) propagate Lipschitz continuous i.e.

[P k ] Lip = sup [g] Lip ≤1
[P k g] Lip < +∞, k = 0, . . . , n.

Set [P ] Lip = max 0≤k≤n-1 [P k ] Lip and [f ] Lip = max 0≤k≤n [f k ] Lip . Let p ∈ [1, +∞). We assume that n k=1 X k L p (P) + X k L p (P) < +∞.
(a) For every k ∈ {0, . . . , n},

U k -U k L p (P) ≤ 2[f ] Lip n ℓ=k [P ] Lip ∨ 1 n-ℓ X ℓ -X ℓ L p (P) . (b) If p = 2, for every k ∈ {0, . . . , n}, U k -U k L 2 (P) ≤ √ 2[f ] Lip n ℓ=k [P ] Lip ∨ 1 2(n-ℓ) X ℓ -X ℓ 2 L 2 (P) 1 2
.

(1.25)

Remark. In fact the error bounds established in this theorem remain mathematically true as soon as the "quantizations" X k are simply σ(X k )-measurable i.e. when X k = π k (X k ) where π k : R d → R d is a Borel functions of X k for every k ∈ {0, . . . , n}, provided the integrability assumptions on X k which give sense to the QBDP P as defined in (1.23) are satisfied. If the functions π k take countably many values, even (1.24) still makes sense; otherwise sums should be replaced by integrals with respect to conditional measures and the computational tractability is usually lost.

Proof.

Step 1. First, we control the Lipschitz constants of the functions u k . It follows form the elementary inequality

| sup i∈I a i -sup i∈I b i | ≤ sup i∈I |a i -b i |, a i , b i ∈ R, i ∈ I, that [u k ] Lip ≤ max [f k ] Lip , [P k u k+1 ] Lip ≤ max [f ] Lip , [P k ] Lip [u k+1 ] Lip with the convention [u n+1 ] Lip = 0. A straightforward backward induction yields [u k ] Lip ≤ [f ] Lip [P ] Lip ∨ 1 n-k . (1.26)
Step 2. We focus on claim (b) (quadratic case p = 2). First, we derive from Proposition 1.2(a) applied to X = X k+1 and

Y = X k , Q = P k and f = u k+1 , g = u k+1 and h = u k that E U k+1 |X k -E U k+1 | X k 2 L 2 (P) ≤ [P k u k+1 ] 2 X k -X k 2 L 2 (P) + u k+1 (X k+1 ) -u k+1 ( X k+1 ) 2 L 2 (P) . (1.27)
Now, it follows from both original and quantized dynamic programming formulas that

|U k -U k | ≤ max |f k (X k ) -f k ( X k )|, E U k+1 |X k -E U k+1 | X k so that |U k -U k | 2 ≤ |f k (X k ) -f k ( X k )| 2 + E U k+1 |X k -E U k+1 | X k 2 .
Taking expectation and plugging (1.27) in the above inequality yields for every k ∈ {0, . . . , n -1},

U k -U k 2 L 2 (P) ≤ [f ] 2 Lip + [P ] 2 Lip [u k+1 ] 2 Lip X k -X k 2 L 2 (P) + U k+1 -U k+1 2 L 2 (P)
still with the convention [u n+1 ] Lip = 0. Now, using (1.26), we obtain

[f ] 2 Lip + [P ] 2 Lip [u k+1 ] 2 Lip ≤ [f ] 2 Lip + [P ] 2 Lip 1 ∨ [P ] Lip 2(n-(k+1)) ≤ 2[f ] 2 Lip 1 ∨ [P ] Lip 2(n-k) . Consequently U k -U k 2 L 2 (P) ≤ 2 n-1 ℓ=k [f ] 2 Lip 1 ∨ [P ] Lip 2(n-ℓ) X ℓ -X ℓ 2 L 2 (P) + [f ] 2 Lip X n -X n 2 L 2 (P) ≤ 2[f ] 2 Lip n ℓ=k 1 ∨ [P ] Lip 2(n-ℓ) X ℓ -X ℓ 2 L 2 (P)
which completes the proof. Claim (a) can be established following the above lines and Claim (b) of Proposition 1.2.

Example of application: the Euler scheme. Let ( Xn

t n k
) 0≤k≤n be the Euler scheme with step T n of the d-dimensional diffusion (X 0,x t ) t∈[0,T ] , solution to the SDE (1.17). It defines a homogeneous Markov chain with transition

P n k g(x) = E g x + T n b(t n k , Xn t n k ) + σ(t n k , Xn t n k ) T n Z , Z L ∼ N (0, I q ). If f is Lipschitz continuous, P n k g(x) -P n k g(x ′ ) 2 ≤ [g] 2 Lip E x -x ′ + T n b(t n k , x) -b(t n k , x ′ ) + T n σ(t n k , x) -σ(t n k , x ′ ) Z 2 ≤ [g] 2 Lip x -x ′ + T n b(t n k , x) -b(t n k , x ′ ) 2 + σ(t n k , x) -σ(t n k , x ′ ) 2 T n ≤ [g] 2 Lip |x -x ′ | 2 1 + T n [σ] 2 Lip + 2T n [b] Lip + T 2 n 2 [b] 2 Lip
where A = Tr(AA * ), Tr stands for the trace of a square matrix, A * stands for the transpose of the d × q-matrix A. The coefficient [σ] Lip should be understood as the Lipschitz coefficient of σ : (R d , | . |) → (M(d, q, R), . ). As a consequence

[ P n k g] Lip ≤ 1 + C b,σ,T T n [g] Lip , k = 0, . . . , n -1,
where

C b,σ,T = [b] Lip + 1 2 [b] 2 Lip T + [σ] 2 Lip i.e. [ P n ] Lip ≤ 1 + C b,σ,T T n .
Let ( Ūk ) 0≤k≤n denote the (F W

U k -U k L 2 (P) ≤ √ 2[f ] Lip n ℓ=k 1 + C b,σ,T T n 2(n-ℓ) X ℓ -X ℓ 2 L 2 (P) 1 2 ≤ √ 2[f ] Lip n ℓ=k e 2C b,σ,T (T -t n ℓ ) X ℓ -X ℓ 2 L 2 (P) 1 2
(1.28)

≤ √ 2[f ] Lip e C b,σ,T T n ℓ=k X ℓ -X ℓ 2 L 2 (P) 1 2
.

(1.29)

The fact of interest concerning the upper bounds(1.29) and (1.29) that should be emphasized is that the real constants involved in them only depend on b, σ and T but not on n (except for the sum itself of course). We will see further on that the penultimate inequality can be used to "calibrate" the sizes of the quantization grids associated to the X k (see section 2.3.4). Finally note that the above computations hold more generally for the Euler scheme of a diffusion driven by a Lévy processes Z with Lévy measure ν satisfying ν(z 2 ) < +∞ (i.e. Z t ∈ L 2 for every t ∈ [0, T ]).

Connection with parabolic variational inequalities. We consider the parabolic variational inequality on [0, T ] × R d defined by

max f -u, ∂u ∂t + Lu (t, x) = 0, (t, x) ∈ [0, T ) × R d , u(T, . ) = f (T, . ) (1.30)
where L is the operator introduced in (1.16) (which is for the probabilist the generator of the diffusion (1.17)). We assume that f : [0, T ] × R d → R + is (at least) continuous with polynomial growth in the space variable x, uniformly in t ∈ [0, T ].

From now on we will switch to a completely heuristic reasoning in order to highlight in a simpler way the connection between the above variational inequality and optimal stopping theory in continuous time. This connection holds through a probabilistic representation formula involving the diffusion process (1.17) in the same spirit as that which holds for parabolic P DEs through the Feynman-Kac formula. A probabilistic representation of this variational inequality (1.30) is provided, under appropriate conditions that we will not detail here (see [START_REF] Bensoussan | Applications des inéquations variationnelles en contrôle stochastique[END_REF][START_REF] Bensoussan | On the theory of option pricing[END_REF] or, more recently, [START_REF] Lamberton | Optimal stopping and American options[END_REF] and the references therein for a rigorous presentation in various settings), by the continuous time optimal stopping problem, related to the diffusion process (X x,t ) s∈[t,T ] solution to Equation (1.17) and the obstacle process

Z t := f (t, X x,0 t ), t ∈ [0, T ]
. This obstacle process is F W -adapted and continuous (hence predictable) where F W denotes the augmented filtration of the Brownian motion W . The function f having polynomial growth in x uniformly in t ∈ [0, T ], then sup t∈[0,T ] Z t ∈ L 1 (P). We define the P-Snell envelope (U t ) t∈[0,T ] by

U t = P-esssup E f (τ, Z τ ) | F t , τ ∈ T W t,T (1.31)
where T W t,T denotes the set of F W -stopping times τ : (Ω, A, P) → [0, T ] i.e. [t, T ]-valued random times satisfying ∀ s ∈ [t, T ], {τ ≤ s} ∈ F W s . (This definition implies that for every s ∈ [t, T ], {τ = s} ∈ F W s but the converse -which is required for technical reasons -is usually not true since [0, T ] is not countable.) One shows (see e.g. [START_REF] Shiryaev | Optimal Stopping rules, 2 nd edition[END_REF]) that under these conditions, there exists a function u :

[0, T ] × R d → R + such that U t = u(t, X x,0 t ), t ∈ [0, T ],
and, in terms of réduite,

u(t, x) = sup E f (τ, X x,t τ ), τ ∈ T W t,T , t ∈ [0, T ], x ∈ R d .
Unfortunately, even in simple frameworks, this function u is not smooth enough, say C

= P-esssup E f (τ, Z τ ) | F t n k , τ ∈ T n k,n .
The sequence (U n t n k ) 0≤k≤n is nothing but the P-Snell envelope of (Z t n k ) 0≤k≤n viewed as a discrete time optimal stopping problem associated to the Markov chain (X x,0 t n k ) 0≤k≤n . Its transitions are given by P k (ξ, dy) = P(X t n k+1 ∈ dy | X t n k = ξ), k = 0, . . . , n -1. They are usually not explicit and, more important in practice, even not simulatable: more generally, exact simulation of time samples of a d-dimensional diffusion process is impossible (at a reasonable cost) as soon as d ≥ 2.

Then, keeping in mind that U n

t n n = U n T , the sequence (U n t n k ) 0≤k≤n satisfies the following BDP P U n t n n = f (T, X x,0 T ), U n t n k = max f (X 0,x t n k ), E U n t n k+1 | X 0,x t n k , 0 ≤ k ≤ n -1.
Consequently, one may consider that U t n k is close enough to U n t n k so that the sequence (U t n k ) 0≤k≤n approximately satisfies

U T = f (T, X x,0 T ), max f (X x,0 t n k ) -U t n k , E U t n k+1 | X x,0 t n k -U t n k ≈ 0 or equivalently u(T, X x,0 T ) = f (T, X x,0 T ), 0 = max f (X 0,x t n k ) -u(t n k , X x,0 t n k ), E u(t n k+1 , X x,0 t n k+1 ) | X 0,x t n k -u(t n k , X x,0 t n k ) .
If, as a second step, we assume that, in fact, the function u has the requested regularity C 1,2 to apply Itô's Lemma between t n k and t n k+1 , and, moreover, that u has polynomial growth in x, uniformly in t ∈ [0, T ]. Then one obtains u(t n k+1 , X x,0

t n k+1 ) = u(t n k , X x,0 t n k ) + t n k+1 t n k ∂u ∂t + Lu (s, X x,0 s )ds + t n k+1 t n k ∇ x u | σ (s, X x,0 s dW s
where the stochastic integral with respect to W is a true martingale under our assumptions on b, σ (linear growth) and f (polynomial growth). As a consequence,

E u(t n k+1 , X x,0 t n k+1 ) | X t n k -u(t n k , X x,0 t n k ) = E t n k+1 t n k ∂u ∂t + Lu (s, X x,0 s )ds P-a.s. ≈ ∂u ∂t + Lu (t n k , X x,0 t n k ) T n P-a.s.
Finally, this leads to

u(T, X x,0 T ) = f (T, X x,0 T ), max f (X t n k ) -u(t n k , X x,0 t n k ), ∂u ∂t + Lu (t n k , X x,0 t n k ) ≈ 0 P-a.s..
Letting n → +∞, one may reasonably guess, switching back to continuous time that, P-a.s.,

u(T, X x,0 T ) = f (T, X x,0 T ), max f (X t ) -u(t, X x,0 t ), ∂u ∂t + Lu (t, X x,0 t ) = 0.
If the support of X x,0 t is the whole space R d for every t ∈ (0, T ), this suggests that u is a solution to the above variational inequality (1.30).

The simplest setting in which the above heuristic reasoning is not trivial can be made rigorous is when X x,0 t = x + W t , relying on distribution theory and extended versions of Itô's lemma (see [START_REF] Lamberton | Optimal stopping and American options[END_REF] for a rigorous variational approach using distribution theory).

In practice, this approach is usually used in a backward way as a verification theorem: starting from a solution u to (1.30) (possibly in a viscosity sense) with enough regularity, one checks that (u(t, X x,0 t )) t∈[0,T ] is the P-Snell envelope of the payoff/obstacle process Z t = f (t, X x,0 t ) (see again [START_REF] Lamberton | Optimal stopping and American options[END_REF] and [START_REF] Bensoussan | Applications des inéquations variationnelles en contrôle stochastique[END_REF]).

From a numerical point of view, this connection suggests numerical methods to solve (1.30), at least in higher dimension when deterministic schemes become inefficient due to the curse of dimensionality: any probabilistic method devised to compute the réduite of a discrete time optimal stopping problem with a Markov structure process can be viewed as a way to solve (1.30). Having in mind that, in these probabilistic methods, the discrete time Markov process needs to be simulated at a moment or another, so we cannot rely on a sample (X t n k ) 0≤k≤n of the diffusion process itself at least when d ≥ 2 ( 4), with of course, the noticeable exception of diffusion processes of the form X x,0 t = ϕ(t, x, W t ) like e.g. the Black-Scholes model.

Consequently, for numerical purposes, an additional step is necessary : one has to introduce the discrete time Euler scheme Xn,x,0 = ( Xn,x,0

t n k
) 0≤k≤n with step T n (with obvious notations) of the diffusion (X x,0 t ) t∈[0,T ] and the Snell envelope Ū n = ( Ū n t n k ) 0≤k≤n related to the payoff/obstacle process Zn

t n k = f (t n k , Xn,x,0 t n k
), k = 0, . . . , n. All these approximations of the "true" original Snell envelope have been investigated in full details in [START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF], depending on the regularity of the payoff function and of the coefficients b and σ (or a = σσ * ).

Note that for the sake of simplicity, we did not take into account a zero order term in (1.30). This is done in the cited literature and corresponds to a discounting factor appearing in the payoff/obstacle process.

The resulting convergence rates are of the form O(n -1 2 ) between ( Ū n t n k ) 0≤k≤n and (U n t n k ) 0≤k≤n and, under a quasi-convexity assumption (in the space variable) on f , O(n -1 ) between (U n t n k ) 0≤k≤n and (U t n k ) 0≤k≤n This holds for various type of convergence (on compacts set with respect to x or in L p (P)).

Conversely, we must of course emphasize that in low dimension (say d ≤ 3), the numerical strategy is the exact converse of what precedes: one relates the optimal stopping theorem to the parabolic variational inequality and implement a deterministic solver based on finite difference or finite elements/volume methods. So is the case for the seminal paper on the computation of American Put option in a Black-Scholes model, see [START_REF] Lamberton | Variational inequalities and the pricing of American options[END_REF] (see also [START_REF] Villeneuve | Parabolic A.D.I. methods for pricing american option on two stocks[END_REF] for exchange options). For an overview of P DE methods for option pricing, we refer to [START_REF] Achdou | Computational Methods for Option Pricing[END_REF].

Optimal vector quantization

As emphasized in what precedes, whatever the application fields are, we need to have "good" and, if possible, "optimal", quantization grids at various levels N . This section is devoted to the existence of optimal quantizers (via the proof of Proposition 1.1(b)) and their properties as well as the numerical methods to compute them.

Existence of optimal quantization grids (proof of Proposition 1.1(a))

First assume that p ∈ (0, +∞) and X ∈ L p R d (P). We will proceed by induction on the level N . First we note that, for every level N ≥ 1, the L p -mean quantization error function defined by

(x 1 , . . . , x N ) -→ min 1≤i≤N |X -x i | L p (P)
is clearly 1-Lipschitz continuous with respect to the ℓ ∞ -norm on (R d ) N defined by |(x 1 , . . . , x N )| ℓ ∞ := max 1≤i≤N |x i |. This is a straightforward consequence of Minkowski's inequality and the more elementary inequality | min 1≤i≤N a imin 1≤i≤N b i | ≤ max 1≤i≤N |a ib i |. As a consequence, the L p -distortion function defined for every x 1 , . .

. , x N ) ∈ (R d ) N by G p,N (x 1 , . . . , x N ) = min 1≤i≤N |X -x i | p L p (P)
is continuous.

If p ∈ (0, 1], on directly shows that G p,N is p-Hölder continuous.

⊲ N = 1. The non-negative continuous function G p,1 clearly goes to +∞ as |x 1 | → +∞. Hence, G p,1 attains a minimum at a so-called L p -median x (1) (which is clearly unique when p > 1 by a strict convexity argument). So {x

(1)
1 } is an optimal quantization grid at level 1.

⊲ N ⇒ N + 1. Assume there exists x (N ) ∈ (R d ) N such that G p,N (x (N ) ) = min (R d ) N G p,N . Set Γ (N ) = {x (N )
i , i = 1, . . . , N }. Then, either supp(µ) \ Γ (N ) = ∅ and any N + 1-tuple of (R d ) N which "exhausts" the grid Γ (N ) makes the function G p,N +1 equal to 0 which is its lowest possible value, or there exists ξ N +1 ∈ supp(µ) \ {x

λ N +1 = E |X -X Γ * | p = E |X -ξ N +1 | p ∧ min 1≤i≤N |X -x i | p < E min 1≤i≤N |X -x i | p = e p,N (X) p .
Hence, the set

K N +1 = x ∈ (R d ) N : G p,N (x) ≤ λ N +1
is non-empty since it contains all the N + 1-tuples which "exhaust" the elements of Γ * and closed since G p,N is continuous. In fact, we will show that it is also a bounded subset of (R d ) N . Let

x(k) = (x(k) 1 , . . . , x(k) N +1 ), k ∈ N, be a K N +1 -valued sequence of N + 1-
G p,N +1 (x(k)) ≥ lim inf k min i∈I |X -x i (∞)| p L p (P) ≥ e p,|I| (X) p .
The fact that the sequence x(k) k∈N is K N +1 -valued implies that e p,|I| (X) ≤ λ N +1 < e N (X). In turn, this implies that |I| = N + 1 i.e. the sequence of N + 1-tuples x(k) k≥0 is bounded. As a consequence, the set K N +1 is compact and the function G p,N +1 attains a minimum over K N +1 which is obviously its absolute minimum.

Remarks. • The proof of Claim (b) and (c) are as follows: if N ≥ 2 and P P X ∈ C i (Γ (N ) = 0 for some index i, then e p,N (X) = e p,N -1 (X) which is impossible by the above proof if supp(µ) has at least N points. Moreover, if (z N ) N ≥1 is an everywhere dense sequence in R d , then 0 ≤ e p,N (X) ≤ e p {z 1 , . . . , z N } → 0 as N → +∞ by the Lebesgue dominated convergence theorem.

• This existence result admits many extensions, in particular in infinite dimension when R d is replaced by a separable Hilbert space or, more generally, a reflexive Banach space. It holds true for L 1 -spaces as well (see [START_REF] Graf | Optimal quantizers for Radon random vectors in a Banach space[END_REF] or Section 2.3.2 for a brief introduction to functional quantization).

• When d ≥ 2, argmin G p,N is never reduced to a single N -tuple, simply because argmin G p,N is left stable under the action of the permutations of {1, . . . , N }. Even from a geometrical viewpoint, as soon as d ≥ 2, uniqueness of the corresponding quantization grid usually fails, e.g. because of the transformations leaving the distribution µ of X invariant. So is the case for the normal distribution N (0; I d ) which is invariant by all orthogonal transforms. But there are also examples (see [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]) for which optimal grids at level N do not even make up a "connected" set.

However, in 1-dimension, it has been proved (see e.g. [START_REF] Kieffer | Exponential rate of convergence for Lloyd's method I[END_REF]) that, as soon as µ is absolutely continuous with a log-concave density, there exists exactly one optimal quantizer at level N (characterized by its stationarity, see the next section). Such distributions recalled strictly unimodal.

Stationary quantization grid

Throughout this section, we still consider a random vector X : (Ω, A, P) → R d with distribution µ on (R d , Bor(R d )). (ii) Self-consistency of the centroids : ∀ i ∈ {1, . . . , N },

x i = C i (Γ) ξµ(dξ) µ(C i (Γ)) .
Note that item (ii) is equivalent to

X Γ = E X | X Γ . (2.32)
We will see further on in Section 3.1 (see Corollary 3.1 and the comments that follow) that L 2optimal grids are always stationary since they satisfy condition (i) (see Theorem 4.3, p.38, in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]), and are consequently critical points of the quadratic distortion function G 2,N which implies condition (ii).

Note that almost all numerical algorithms devised to search for (at least locally) optimal quantizers at a given level N are based on this stationary property (see Section 3).

Remark. This property is important for applications but is unfortunately satisfied by few quantizers. A new notion of quantization, called dual quantization has been recently developed (see [START_REF] Pagès | Intrinsic stationarity for vector quantization: Foundation of dual quantization[END_REF] for a theoretical introduction and [START_REF] Pagès | Dual Quantization for random walks with application to credit derivatives[END_REF][START_REF] Pagès | Optimal Delaunay et Voronoi quantization methods for pricing American options, Numerical methods in Finance[END_REF] for applications to Numerical Probability) in which a reverse stationarity property is satisfied by all dual quantizers. Typically for dual quantization, one has

X = E X Γ,dual | X .
This dual quantization requires the introduction of an exogenous noise and relies on the Delaunay triangulation rather than on the Voronoi diagram (hence its name "dual").

2.2.1 Application to quantization based cubature formula (II): second order Proposition 2.1. Let Γ be a µ-stationary quantizer and let X : (Ω, A, P) → R d be a µ-distributed random vector such that E |X| 1+α < +∞ with distribution µ with moment 1 + α. Let F : R d → R be a C 1 function with an α-Hölder gradient ∇F . We denote by [∇F ] α its α-Hölder ratio. Then

E F (X) | X Γ -F ( X Γ ) 1 ≤ [∇F ] α X -X Γ 1+α L 1+α (P)
so that, by taking expectation,

E F (X) -E F ( X Γ ) ≤ [∇F ] α X -X Γ 1+α L 1+α (P) (2.33) In particular, if ∇F is Lipschitz continuous E F (X) -E F ( X Γ ) ≤ [∇F ] Lip X -X Γ 2 L 2 (P) .
Remark. Since we know from Claim (c) in Proposition 1.1 that e 2,N (X) ↓ 0 as N → +∞, it is clear that if we consider optimal quadratic (hence stationary) grids Γ

(N ) satisfying e 2,N (X) = X -X Γ (N ) L 2 (P) which implies that X -X Γ (N ) 1+α L 2 (P) = o X -X Γ (N )
L 2 (P) as N → +∞. So, the error bounds in the above cubature formulas become infinitely smaller than those obtained in Section 1.3.2. Moreover, a sharp convergence rate of convergence for e 2,N (X) is proved in Section 2.3 below (Theorem 2.1, Zador's Theorem).

Proof. It follows from a second order Taylor expansion that

F (X) -F ( X Γ ) -∇F ( X Γ )|X -X Γ ≤ [∇F ] α X -X Γ 1+α .
The stationarity property reads X Γ = E X | X Γ ) and, as ∇F ( X Γ ) is σ( X Γ )-measurable, we obtain

E ∇F ( X Γ )|X -X Γ | X Γ = ∇F ( X Γ )| E X -X Γ =0 | X Γ ) = 0.

Application to convex functions

If F : R d → R is convex and X, F (X) ∈ L 1 (P), then for every stationary quantization grid, one has

E F ( X Γ ) = E F E X | X Γ ≤ E E F (X) | X Γ = E F (X).
This shows that quantization based cubature formulas provide lower bounds for the true expectation when F is convex. It should be noticed, that, one has the reverse inequality by the dual quantization based cubature formula E F (X) ≥ EF X Γ,dual

(see [START_REF] Pagès | Intrinsic stationarity for vector quantization: Foundation of dual quantization[END_REF]) which can be implemented with the same stationary grid Γ.

Optimal quantization rate

Rate results

This rate is ruled by two results: the first one is asymptotic and is known as Zador's Theorem (see [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]), the second one is non-asymptotic and, in some way, universal.

Theorem 2.1. (a) Zador's Theorem (asymptotic rate, see [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]). Let X ∈ L p+η R d (P) for some η > 0, with distribution P X = µ. Let | . | denote (here) any norm on R d . We define ϕ = dµ dλ d (where λ d denotes the Lebesgue measure on (R d , Bor(R d ))) the density of the non-singular part of µ. Then 3 and ( 5)

J p,d ∼ d 2πe 1 2
as d → +∞.

(2.34) (b) Extended Pierce's Lemma (non-asymptotic rate, see [START_REF] Luschgy | Functional quantization rate and mean regularity of processes with an application to Lévy processes[END_REF]). Let p, η > 0. There exists a real constant κ d,p,η ∈ (0, +∞) such that, for every R d -valued random vector X : (Ω,

A, P) → R d , ∀ N ≥ 1, e p,N (X) ≤ κ d,p,η σ p+η (X)N -1 d (2.35)
where, for every r ∈ (0, +∞), σ r (X) = inf a∈R d Xa L r (P) denotes the L p -standard deviation of X.

Claim (a), known as Zador's Theorem, was first established in [START_REF] Zador | Asymptotic quantization error of continuous signals and the quantization dimension[END_REF] for the uniform distribution on the unit hypercube [0, 1] d , using a self-similarity argument "à la Hammersley". The extension to more general distributions is due to Bucklew and Wise in [START_REF] Bucklew | Multidimensional Asymptotic Quantization Theory with r th Power Distortion Measures[END_REF], with a gap in the proof for non-compactly supported distributions. This gap was filled in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF], Chapter II, Theorem 6.2, p.78, where the most general case ("L p+ "-integrable random vectors) has been finally rigorously established.

When µ is purely singular, the above result remains true but the above asymptotics is not the right one since the limit in this normalization is 0. Other asymptotics can emerge e.g. for (uniform) distributions on fractal sets, leading to the definition of a quantization dimension (see again [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF], Chapter III, for an introduction to quantization on fractal sets).

Comments and first applications

⊲ Comments. (a) A natural question is to compare optimal quantization methods with "standard" square grid methods on [0, 1] d when N = m d , m ∈ N * (such values of N are the only ones for which it is possible to implement the square grid cubature formulas). Let

Γ d (m) = d ℓ=1 2i ℓ -1 2m
1≤i ℓ ≤m be a hyper-cubic "product quantization" grid. This grid induces, when p = 2, owing to the Pythagorus Theorem and the formula (1.4) for the uniform distribution U ([0, 1]) on the unit interval, a mean quadratic quantization error given by

e 2,N Γ d (m), U ([0, 1] d ) 2 = d 1 2 √ 3 m 2 = d 12m 2 so that e 2,N Γ d (m), U ([0, 1] d ) = d 12 N -1 d .
If we consider that d is large enough so that, by (2.34), Jp,d ≈ d 2πe 1 2 , one derives by the above Zador Theorem,

e 2,N U ([0, 1] d ) ≈ d 2πe N -1 d ≈ d 17.0795 N -1 d .
This shows that, as d grows to infinity, the ratio between true optimal and optimal product quantization remains (asymptotically) bounded and, in fact, very close to 1 since 2πe 12 ≈ 1.1930 . . . ! But the main difference between these two points of view is that optimal quantization grids do exist for every size/level N whereas product grids become increasingly sparse in higher dimensions.

(b) If Γ N = {x (N ) 1 , . . . , x (N )
N , N ≥ 1, denotes a sequence of (L p , µ)-optimal N -quantizers, we saw in Equation (1.12) that the weighted distribution μN of X Γ (N ) weak converges to the original distribution µ of X. But what about the regular empirical measures

µ N = 1 N N i=1 δ x (N ) i , N ≥ 1 ?
The answer is provided by the empirical measure theorem established by Graf & Luschgy in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF] (see Theorem 7.5, p.96) and [START_REF] Delattre | Quantization of probability distributions under norm-based distortion measures[END_REF].

Theorem 2.2. If X ∈ L p+η (P) for some η > 0, then, with the notations of Zador's Theorem, if ϕ ≡ 0,

µ N (w) -→ µ = ϕ d d+p R d ϕ d d+p dλ d .λ d as N → +∞.
where

-→ denotes weal convergence of probability measures.

Remark. In fact the theorem holds true for (sub-)sequences of asymptotically optimal quantizers in the sense that they satisfy Zador sharp rate).

⊲ A brief look back at Information theory. Let X : (Ω, A, P) → R d be a (static) square integrable random vector/signal X with a non-singular distribution µ. We consider an optimal quadratic quantization grid Γ (N ) of size N for µ. Let μN denote the distribution of its optimal quadratic quantization X Γ (N ) We saw in Section 1.2 that the mean transmission cost κ μN , when coded by an appropriate instantaneous code, satisfies the two-sided inequality (1.6)-(1.7) involving the dyadic entropy of H 2 (μ N ) of μN . ⊲ A universal asymptotic bound. Let w

(N ) i = μN x (N ) i
, i = 1, . . . , N , be the weights attached to the distribution μN . We know, by classical optimization arguments, that

H 2 μN = - N i=1 w (N ) i log 2 w (N ) i ≤ log 2 N so that we derives that 2 κ μN -1 ≤ N ≤ 2
κ μN . Now, let κ > 0 be a prescribed mean transmission cost. Plugging this universal bound in (2.1) with N = ⌈2 κ-1 ⌉, we derive that the quadratic mean transmission error X -X Γ N L 2 (P) satisfies

L 2 -mean transmission error(κ) J 2,d ϕ 1 2 L d d+2 (λ d ) 2 -κ-1 d as κ → +∞.
In particular, expressed in terms of quantization level, we get this coarser equivalent

-log 2 L 2 -mean transmission error(N ) ∼ κ(N ) d as N → +∞.
This provides the flavor of Shannon's Source Coding Theorem (for a general formulation and a proof in a more general framework, we refer to [START_REF] Dembo | Large deviations techniques and applications[END_REF]).

⊲A sharper, but heuristic, bound. It is commonly shared in the Information Theory community (see [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF]) that an absolutely distribution µ with a continuous bounded density ϕ > 0 satisfies (ϕ log ϕ is bounded and)

w (N ) i ∼ c N ϕ d d+2 (x (N ) i ) N as N → +∞
uniformly over the elementary quantizers x

(N ) i
lying in a fixed compact set, where c N is a normalizing constant. Partial results in that direction have been recently established in [START_REF] Graf | The local quantization behaviour of absolutely continuous probabilities[END_REF] for a wide class of absolutely continuous distributions (including among many others normal distributions, gamma distributions, hyper-exponential distributions).

It follows from the empirical measure theorem (2.2) that lim

N c N = R d ϕ d d+2 (ξ)dξ -1
.

If we assume a sharper result, namely w

(N ) i -c N ϕ d d+2 (x (N ) i ) N
→ 0 as N → +∞, uniformly on compact sets, elementary computations based on repeated applications of Theorem 2.2 -the empirical mean theorem -yield that

C µ = -lim N H 2 μN -log 2 N = log R d ϕ 2 d+2 dλ d + d d + 2 R d ϕ log 2 ϕ dλ d .
Noting that

R d ϕ 2 d+2 dλ d = R d ϕ -d d+2
dµ and applying Jensen's Inequality to the distribution µ and the concave function log imply that C µ > 0 (in accordance with the fact that the uniform distribution µ U nif N over {1, . . . , N } or any set of size N has the highest possible entropy equal to log 2 N ). This finally yields the slightly sharper bound for the quadratic mean signal transmission error

L 2 -mean transmission error(κ) J 2,d ϕ 1 2 L d d+2 (λ d ) 2 1-Cµ d 2 -κ d as κ → +∞.
⊲ From vector to functional quantization. Most questions raised in the former sections can be formulated in an infinite dimensional setting where R d is replaced by a separable Hilbert space or more generally a (reflexive or not) separable Banach space (E, | . | E ). From a more probabilistic point of view E is often a functional space like L 2 R ([0, T ], λ) or C([0, T ], R). This corresponds to the path space of a bi-measurable stochastic process (X t ) t∈[0,T ] defined on a probability space (Ω, A, P) so that for every ω ∈ Ω, X(ω) = (t → X t (ω)) ∈ E: viewing things that way makes X an E-valued (infinitely dimensional) random vector. For this reason, the quantization problem in infinite dimension is often known as functional quantization. Thus, existence of optimal quantizers can be established in the case of reflexive Banach spaces like L p R ([0, T ], dt), 1 < p < +∞, following the proof of Proposition 1.1, using this time a weak continuity lower semi-continuity argument (see e.g. [20, 45, ?]). it also holds true for the L 1 R (dt) space, but not in full flu generality for the space C([0, T ], R). However, no such general rate result as Zador's Theorem is available, but optimal quantization rates (sometimes sharp) can be established for various classes of stochastic processes, including the (fractional) Brownian motion and the Brownian bridge, scalar Brownian diffusions processes, etc, viewed as random variables having values on their path spaces. We refer to [START_REF] Luschgy | Functional quantization of stochastic processes[END_REF][START_REF] Luschgy | Sharp asymptotics of the functional quantization problem for Gaussian processes[END_REF][START_REF] Dereich | The coding complexity of diffusion processes under L p ([0, 1])-norm distortion, Stochastic Process[END_REF][START_REF] Luschgy | Functional quantization rate and mean regularity of processes with an application to Lévy processes[END_REF][START_REF] Graf | Optimal quantizers for Radon random vectors in a Banach space[END_REF] and the references therein for results in that direction. It turns out that in many situations these convergence rates for optimal functional quantization error "lives" in a (log N ) -r scale (r = 1 2 for the standard Brownian motion and bridge, diffusion processes, r = H for the fractional Brownian motion with Hurst constant H, r = m + 1 2 for the m-fold integrated Brownian motion, etc). More generally, when dealing with Gaussian processes this quantization rate is closely related to the rate of decay of the eigenvalues of its covariance operator (Karhunen-Loève eigen-system). For more general stochastic processes, a connection can be made between the mean L p ([0T ), dt)-functional quantization rate and the regularity of the mapping t → X t from [0, T ] to L p ([0, T ], dt).

Numerical applications have also been devised, mostly based on the cubature formulas (see below) which can be straightforwardly extended to Hilbert or Banach frameworks (see e.g. [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF] for applications to the pricing and hedging of path-dependent options). But we will not go further in that direction in this paper.

Numerical integration (III): rates

Assume X ∈ L 2 R d (P) (still with a distribution µ). The above results show that, provided one has access to quadratic optimal quantization grids Γ N for various levels/sizes, it is possible to produce cubature formulas with respect to µ which have the following properties.

⊲ α-Hölder functions, α ∈ (0, 1]. An α-Hölder function F : R d → R can be integrated using the Npoint cubature formula (1.10) with an accuracy (upper-bounded by) X -X Γ N L 1 (P) owing to (1.12).

Let Γ (N ) be an optimal quadratic quantization grid at level N for the distribution µ and that X ∈ L 2+η R d (P) for some η > 0. We use that X - ⊲ Locally α-Hölder functions. Let us come back to the general error bound (1.11) for α-Höder functions with α + β-polynomial growth at infinity. If α + β ≥ 1 then β 1-α ≥ 1 so that, combining stationarity and Jensen's inequality, we get

X Γ N L 1 (P) ≤ X -X Γ (N ) L 2 (P) = e 2,N (µ) ≤ C d,µ N -1
X Γ (N ) L β 1-α (P) = E X | X Γ (N ) L β 1-α (P) ≤ X L β 1-α (P)
which makes the right hand side of (1.11) only depending on (norms of) X, namely

E F (X) -E F ( X Γ (N ) ) ≤ [F ] α,β C d,µ 1 + 2 X β L β 1-α (P) N -α d .
⊲ Functions with α-Hölder gradient. If F is continuously differentiable on R d with an α-Hölder gradient α ∈ (0, 1] and X ∈ L 2+η R d (P), note that in the cubature formula (2.33), one has

X -X Γ (N ) L 1+α (P) ≤ X -X Γ (N ) L 2 (P) = e 2,N (µ) = O N -1 d )
since the grid Γ (N ) is stationary (i.e. satisfies (2.32)) so that

E F (X) -E F ( X Γ N ) ≤ C d,µ [∇F ] α N -1+α d . (2.36)
This last result somewhat extends the one established in [START_REF] Mirebeau | Greedy bisection generates optimally adapted triangulations[END_REF] for compactly supported absolutely continuous distributions µ on R d . ⊲ Numerical tests. For various numerical tests carried out on these cubature formulas, including the spatial Richardson-Romberg extrapolation not developed here, we refer to [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF] and the survey [START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF].

Application of (quadratic) optimal quantization to optimal stopping problems

Assume that the characteristics of the quantization tree (Γ k , w k ) 0≤k≤n related to the Markov chain (X k ) 0≤k≤n is made up with grids such that Γ k = Γ

(N k ) k of size N k is L 2 -optimal for the marginal distribution µ k of X k at every time k ∈ {0, . . . , n}. Let N = N 0 + • • • + N n
denote the total number of points (or nodes) used in the quantization tree.

If one sets N k = N := N n+1 , (or N -1 n if X 0 is deterministic), it follows from the upper-bound (1.25) that the (optimally) quantized approximation of the Snell envelope related to the optimal stopping problem (1.31) and their related réduites satisfy

E u 0 ( X Γ 0 0 ) -E u 0 (X 0 ) ≤ U 0 -U 0 L 2 (P) ≤ C X √ n N -1 d .
In the special case where the Markov chain is the Euler scheme ( Xn

t n k
) 0≤k≤n with step T n of a diffusion process with drift b and diffusion coefficient σ, then the constant C X = C b,σ,T does not depend on the time discretization parameter n. In particular, if X 0 = x 0 ∈ R d (with the obvious notations u n and u n ) we have

u n 0 (x 0 ) -u n 0 (x 0 ) ≤ C b,σ,T (N -1) -1 d n 1 d + 1 2 .
Although we will not discuss this point here in details, we must have in mind that the time discretization error bounds established in [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF] behave under Lipschitz assumption on b, σ and f (in

x uniformly in t ∈ [0, T ]) is O(n -α ) with α = 1
2 when the discrete time approximating Markov chain is the Euler scheme. When the discrete time approximating Markov chain is the sampled process itself (X t n k ) 0≤k≤n and the function f is quasi-convex with linear growth in x uniformly in t ∈ [0, T ], the rate holds with α = 1 (this class of functions includes convex functions f with linear growth in x and functions having a Lipschitz continuous gradient ∇ x f (in x) uniformly in t in both cases).

In fact, especially in a diffusion framework, it is important for practical purposes to optimize the quantization tree, especially the sizes N k of the grids, subject to the above global constraint 0≤k≤n N k = N by taking advantage of the non asymptotic bound (2.35) for the mean quadratic quantization error i.e.

min N 0 +•••+Nn=N n k=0 N -2 d k Xn,x t n k 2 L 2+η (P)
(with the additional constraint that all the N k s are non-zero positive integers). In practice, Xn

t n k L 2+η (P)
is not known explicitly and one may replace this quantity by its classical upper bound, atypically Xn,x

t n k 2+η ≤ e C b,σ,T t n k (1 + X 0 L 2+η (P)
) which holds true as soon as both functions b and σ satisfy a linear growth assumption (in x, uniformly in t ∈ [0, T ]): this follows from (1.19) after integrating the starting value with respect to the distribution of X 0 (see [START_REF] Bouleau | Numerical methods for stochastic processes[END_REF][START_REF] Pagès | Introduction to Numerical Probability and Applications to Finance[END_REF] among others for details). This (approximately) leads to

N k =     a 2d d+2 k n ℓ=0 a 2d d+2 ℓ N     ∨ 1,
where a k is a known upper-bound of Xn,x t n k L 2+η (P) for every k = 0, . . . , n. The resulting bound is of the form

u n 0 (x 0 ) -u n 0 (x 0 ) ≤ C b,σ,T N -1 d n k=0 a 2d d+2 k 1 2 + 1 d .
⊲ Numerical tests: Extensive numerical experiments on quantization based schemes for pricing American options, but also non-linear filtering, stochastic control have been carried out, including the introduction and analysis of a Richardson-Romberg extrapolation method which dramatically improves the convergence rate. For more details we refer to [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF][START_REF] Bally | Error analysis of the quantization algorithm for obstacle problems[END_REF][START_REF] Bally | A quantization method for pricing and hedging multidimensional American style options[END_REF][START_REF] Pagès | Optimal Delaunay et Voronoi quantization methods for pricing American options, Numerical methods in Finance[END_REF] (the last reference is devoted to both Voronoi and dual vector quantization applied to the the pricing of American style derivatives) and [START_REF] Pagès | Optimal quantization methods and applications to numerical problems in finance[END_REF][START_REF] Pham | Approximation by quantization of the filter process and applications to optimal stopping problems under partial observation[END_REF][START_REF] Corsi | Numerical approximation by quantization of control problems in finance under partial observation, chapter from Mathematical Modeling and Numerical Methods in Finance[END_REF] (for non-linear filtering, stochastic control applied to Finance) and the references therein.

3 How to get optimal (quadratic) quantization?

The foundation of almost every numerical method devoted to the computation of optimal quantizers is to establish the differentiability of the L p -distortion function G p,N (at level N ). In what follows

∂A = Ā\ • A denotes the boundary of A ⊂ R d .
3.1 Differentiability of the (quadratic) distortion function

Proposition 3.1. Let X ∈ L 2 R d (P) with distribution P X = µ. If x = (x 1 , . . . , x N ) ∈ (R d
) N has pairwise distinct components (i.e. x i = x j , i = j) and P X ∈ 1≤i≤N ∂C i (x) = 0 (with an obvious abuse of notation concerning the Voronoi cells), the L 2 -distortion function G 2,N is (finite and) differentiable at x and

∇G 2,N (x) = 2 C i (x) (x i -ξ)µ(dξ) 1≤i≤N = 2 E1 {X∈C i (x)} (x i -X) 1≤i≤N . (3.37)
Note that the µ-stationarity property is also known in the (statistical) literature as the selfconsistency property and such quantizers as self-consistent prototypes (see [START_REF] Tarpey | Principal Points and Self-Consistent Points of Elliptical Distributions[END_REF]). Moreover, this also gave rise to a new terminology (borrowed from statistics) to emphasize that each x i is the µcentre of mass of its Voronoi cell: one speaks of the points x i as the centroids of their Voronoi cells.

Proof. First note that, as the N -tuple x has pairwise distinct components, all the interiors • Ci (x), i = 1, . . . , N , of the Voronoi cells induced by x are non-empty. For every ξ / ∈ 1≤i≤N ∂C i (x), i.e µ(dξ)-a.s., one has for every i ∈ {1, . . . , N },

∂ ∂x i min 1≤j≤N |x j -ξ| 2 = 1 {ξ∈ • Ci(x)} ∂|x i -ξ| 2 ∂x i = 2 1 {ξ∈ • Ci(x)} (x i -ξ).
On the other hand, for every x,

x ′ ∈ (R d ) N , the function G 2,N is locally Lipschitz continuous since |G 2,N (x ′ ) -G 2,N (x)| ≤ R d min 1≤j≤N |x j -ξ| -min 1≤j≤N |x ′ j -ξ| min 1≤j≤N |x j -ξ| + min 1≤j≤N |x ′ j -ξ| µ(dξ) ≤ max 1≤j≤N |x j -x ′ j | R d max 1≤j≤N |x j | + |ξ| + min 1≤j≤N |x j | + |ξ| µ(dξ) ≤ C µ |x -x ′ | ∞ 1 + |x| ∞ + |x ′ | ∞ .
As a consequence, G 2,N is differentiable at x by the local interchange Lebesgue differentiation Theorem.

Remarks. • In fact, when p > 1, the L p -distortion function (with respect to an Euclidean norm) is also differentiable at N -tuple having pairwise distinct components with gradient

∇G p,N (x) = p C i (x) x i -ξ |x i -ξ| |x i -ξ| p-1 µ(dξ) 1≤i≤N = p E 1 {X∈C i (x)} x i -X |x i -X| |x i -X| p-1
1≤i≤N .

An extension to the case p ∈ (0, 1] does exist under appropriate continuity and integrability assumptions on the distribution µ so that µ({a}) = 0 for every a and the function a → R d |ξ -a| p-1 µ(dξ) remains bounded on compact sets of R d . More generally, a differentiation result exist for strictly convex smooth norms (see Lemma 2.5, p.28 in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]).

This leads to the following corollary which is the starting point of all numerical methods to compute optimal (or at least locally optimal) quantizers.

Corollary 3.1. Let X ∈ L 2
R d (P) with distribution µ. Any grid Γ (N ) attached to an N -tuple x (N ) which minimizes the quadratic distortion function G 2,N is a stationary quantizer (at level card(Γ (N ) ) ≤ N ).

By "attached" we mean that

Γ (N ) = {x (N ) i , 1 ≤ i ≤ N }. (b) Any grid Γ N attached to a critical point x (N ) of G 2,N
and whose induced Voronoi partition has a µ-negligible boundary is a stationary quantizer at level (at most) N .

Its main consequence is the following: if card(supp(µ)) ≥ N , any optimal grid has full size N by Proposition 1.1(b) and can be reduced, up to a permutation, to an N -tuple x (N ) . Moreover, owing to Theorems 4.1 and 4.2 in [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF] (p. [START_REF] Lamberton | Variational inequalities and the pricing of American options[END_REF][START_REF] Kieffer | Exponential rate of convergence for Lloyd's method I[END_REF], we know that µ 1≤i≤N ∂C i (x) = 0 and µ C i (x (N ) ) > 0, i = 1, . . . , N , i.e. the boundaries of its Voronoi partition is µ-negligible. Then, one easily checks that x (N ) is in fact a both global and local minimum of G 2,N and G 2,N is differentiable at x (N ) and ∇G 2,N x (N ) = 0 i.e.

∀ i ∈ {1, . . . , N },

C i (x (N )) (x i -ξ)µ(dξ) = 0 or equivalently x (N ) i = C i (x (N )) ξµ(dξ) µ(C i (x (N ) )) = E X | X ∈ C i (x (N ) )
which in turn can be rewritten

X Γ (N ) = E X | X Γ (N ) where Γ (N ) = x (N ) i , i = 1, . . . , N . (3.38)

Competitive Learning Vector Quantization

The Competitive Learning Vector Quantization (CLV Q) is a stochastic gradient descent deriving from the quadratic distortion G 2,N , viewed as a potential function to be minimized. In a deterministic approach, the associated gradient descent is a zero search of its gradient ∇G 2,N . It formally reads

x(k + 1) = x(k) -γ k+1 ∇G 2,N x(k) , x(0) ∈ Hull supp(µ) N , (3.39)
where Hull supp(µ) denotes the closed convex hull of the support of the distribution µ. The sequence (γ k ) k≥1 is a sequence of positive step parameters satisfying the so-called decreasing step assumptions k≥1 γ k = +∞ and k≥1 γ 2 k < +∞. When (or if) ∇G 2,N is itself differentiable (enough) one can replace this step γ k+1 by the inverse of the Hessian of G 2,N , leading to the classical Newton-Raphson procedure (opening the door to many of its variants like Levenberg-Markardt algorithm, etc).

The procedure is (mathematically) well-defined if µ assigns no mass to hyperplanes and x(0) has pairwise distinct components. Indeed, when these conditions are fulfilled, one shows that x i (k + 1) lies in the interior of the Voronoi cell of x i (k) which prevents the components to get stuck as a result of an iteration. Thus, the existence of ∇G 2,N for the next iteration is preserved.

Unfortunately lim inf max i |x i |→+∞ G 2,N (x) < +∞ (one only has lim inf

min i |x i |→+∞ G 2,N (x) = +∞), hence,
the function G 2,N is not a standard potential function for optimization purposes and the classical convergence results existing for this class of zero search procedures (and their variant involving the Hessian like Newton-Raphson) do not apply here: the boundedness of the sequence (x(k)) k≥0 cannot be established, except when the distribution µ has a compact support (so that its closed convex hull is compact too).

In fact, there is an upstream problem: beyond dimension d ≥ 2 or 3, an accurate computation of the components of G 2,N (x) for a given x ∈ (R d ) N becomes too costly to be used in such a recursive procedure as soon as µ has an infinite support (especially if it is atomless/continuous). When computing simple integrals on R d , it is well-known that; as the dimension d increases, one has to switch from cubature formulas to Monte Carlo simulations, provided µ can be simulated. For the same reason and, in the same spirit, we will switch from the above (deterministic) gradient descent to its stochastic counterpart, taking advantage of its representations (3.37) as an expectation, provided the distribution µ can be simulated at reasonable low cost.

Such a stochastic gradient descent formally reads in our setting (the factor 2 is dropped on purpose)

x(k + 1) = x(k) -γ k+1 1 {X k+1 ∈C i (x)} (x i (k) -X k+1 ) 1≤i≤N , x(0) ∈ Hull supp(µ) N ,
x(0) having pairwise distinct components. We simply replaced mutatis mutandis the expectation operator E by a sequence of independent copies (X k ) k≥1 . The iterates (x(k)) k≥0 are now random vectors (though we will still denote them by script letters). When µ assigns no mass to hyperplanes, one shows that the algorithm is well defined i.e. no components of x(k + 1) get stuck as a result of an iteration. As could be expected from the deterministic framework, and for the same reason induced by the behavior of G 2,N at infinity, it is again -we should say all the more! -hopeless to apply the classical counterpart convergence theorems for stochastic gradient like those established e.g. in [START_REF] Duflo | Algorithmes stochastiques, coll[END_REF][START_REF] Kushner | Stochastic approximation and recursive algorithms and applications, 2 nd edition[END_REF][START_REF] Benveniste | Algorithmes adaptatifs and approximations stochastiques[END_REF]. Of course, its main asset is that it can be implemented very easily when µ is simulatable. In fact, under various names (k-means, Competitive Learning Vector Quantization algorithm, nuŕes dynamiques, etc), it has been widely implemented for years in the communities of Artificial Neural Networks, Data Mining and, more recently, Machine Learning, etc, as a clustering procedure producing prototypes and an automatic classifier (see below). Note however that when µ is the empirical measure of a dataset (ξ k ) 1≤k≤n , the standard/original "deterministic" gradient descent (3.39) at a given level N ≪ n can always be implemented since each computation of ∇G 2,N will require finitely many operations involving a known function (the one appearing inside the expectation operator in (3.37)). In this form, the procedure is also known as Forgy's algorithm or batch k-means procedure. However, in practice, it can be too much time consuming when the size n of the dataset is too large. As an alternative, practitioners often come back to the above stochastic gradient by sampling at random uniformly in the dataset the input datum at each iteration.

When implemented at level N ≥ 1 for optimal classification purposes, the resulting procedure is designed to -hopefully -converge toward a(n at least local) minimum of G 2,N or equivalently the quadratic mean quantization error G N ,∈ . If so, it produces N prototypes of the dataset (the resulting N -quantizer). The classification is finally obtained by assigning to each datum of the dataset the label of its Voronoi cell among the prototypes.

At this stage, let us describe more precisely what this CLV Q algorithm really does from a geometric point of view at each iteration.

Stochastic Gradient Descent formula (Quantizer updating). Still starting from a starting N -tuple x(0) = x 1 (0), . . . , x N (0) in Hull supp(µ) N , we update the procedure from k to k + 1 in a two fold procedure which can be analyzed and interpreted as follows: let x(k) := x 1 (k), . . . , x N (k) ∈ (R d ) N be the value of the state vector of the algorithm at the k th iteration (the running vector of "prototypes").

⊲ Competition phase: This phase is also known as the winner selection stage; it amounts to solving the following nearest neighbour search

i win (k + 1) ∈ argmin i∈{1,...,N } x i (k) -X k+1 .
When the input X k+1 falls in a general position (i.e. not on a median hyperplane), i win (k + 1) is uniquely defined. In case of conflict, when X k+1 falls on a median hyperplane and has subsequently (at least) two nearest neighbours, a rule to determine the winner has to be defined. Picking up i win (k + 1) at random among all the nearest neighbours seems the most natural rule. Note that when µ assigns no mass to hyperplanes (e.g. because it is absolutely continuous) this situation P-a.s. never occurs.

⊲ Learning phase: During this phase, the current set of prototypes x(k) = (x i (k)) i=1,...,N is updated to incorporate the information provided by the input X k+1 . As a result the nearest neighbour x i win (k+1) (k) is moved closer to k + 1 by a dilatation centered at X k+1 with a ratio lower than 1. To be more precise   

x i win (k+1) (k + 1) = Dilat(X k+1 , 1γ k+1 )(x i win (k+1) (k))

x i (k + 1) = x i (k) for every i = i win (k + 1).

where Dilat(ξ, ρ) denotes a dilation centered at ξ with ratio ρ ∈ [0, 1] defined by

Dilat(ξ, ρ)(u) = ξ + ρ(u -ξ), u ∈ R d .
First note that, if for every k ∈ N * , the step γ k ∈ (0, 1) then x i win (k + 1) is a convex combination of x i win (k) and X k+1 . Consequently, as the (pairwise distinct) components of x(0) lie in Hull supp(µ) , a straightforward induction shows that so will be the case of all the iterations x(k), k ≥ 0.

As a consequence all these clustering procedures (batch or stochastic) are always implemented under the step assumption ∀ k ∈ N * , γ k ∈ (0, 1).

One deduces that, if the N -tuple x(k) has pairwise distinct components, this feature is preserved by the learning phase. So that the above procedure is well-defined, up to the convention to be made in case of conflict between several components x j (k) in the competitive phase.

The name of the procedure -Competitive Learning Vector Quantization algorithm -is a synthesis of these two basic phases.

The heuristics, rigorously proved in few situations (e.g. when d = 1 and µ has a log-concave density on a bounded interval) is that the procedure a.s. converges, toward an at least local minimum of G 2,N . But one must have in mind that, as soon as the structural dimension d increases, the "landscape" of the potential function G 2,N -the quadratic distortion -has many local minima and various kinds of saddle points. For (partial) theoretical results on the convergence of the CLV Q algorithm, we may refer to [START_REF] Pagès | Voronoi Tessellation, space quantization algorithms and numerical integration[END_REF] (in the special case where X has a compactly supported distribution) and to [START_REF] Pollard | A central limit theorem for k-means clustering[END_REF] for a weak rate of convergence (a Central Limit Theorem under the assumption that the algorithm converges toward the global minimum).

Several L p -variants have also been investigated more recently, in particular the k-median case corresponding to p = 1 (see [START_REF] Cardot | A fast and recursive algorithm for clustering large dataset with k-medians[END_REF]).

Adaptive on line computation of the "companion parameters". Assume that the above CLV Q procedure a.s. converges toward an N -tuple x (N, * ) = x (N, * ) 1 , . . . , x (N, * ) N having values in Hull supp(µ) . One can compute on line the associated weights and the resulting mean quadratic quantization error:

• Weights w (N, * ) i = P( X x * = x (N, * ) i
), i = 1, . . . , N : for every i = 1, . . . , N ,

w i [k + 1] := (1 -γ k+1 )w i [k] + γ k+1 1 {i win =i} a.s.
-→ w (N, * ) i on the event x(k) → x (N, * ) .

• Quadratic Distortion error G 2,N (x) = min 1≤i≤N |X -x i | 2 : G N [k+1] := (1-γ k+1 )G N [k]+γ k+1 |x i win (k+1) (k)-X k+1 | 2 a.s.
-→ G 2,N (x (N, * ) ) on the event x(k) → x (N, * ) .

Note that, since the "ingredients" involved in the above companion procedures are those used in the competition learning phase, they add (almost) no extra C.P.U. time cost, especially if one has in mind (see below) that the costly part of the algorithm (as well as that of the Lloyd I procedure described below) lies in the nearest neighbour search of the "competition phase". For the convergence of the online adaptive version of the companion procedures, we refer to [START_REF] Bally | A quantization algorithm for solving discrete time multidimensional optimal stopping problems[END_REF].

In some way the CLVQ algorithm can be seen as a Non Linear Monte Carlo Simulation devised to design an optimal skeleton of the distribution µ of X, this skeleton being the distribution of X x (N, * ) (with an obvious abuse of notation), that is the N -tuple x (N, * ) itself and its companion weight vector w (N, * ) i 1≤i≤N .

A fixed point algorithm: Lloyd I procedures

⊲ The original Lloyd I procedure. This time, we start directly from the stationary Equation (3.38). The idea is to use this identity to devise a fixed point procedure. Let Γ(0) ⊂ Hull supp(µ) (with N pairwise distinct elements). Then for every integer k ≥ 0, one updates the current grid Γ(k) as follows:

       (i)Centroid updating: X Γ(k+1) = E X | X Γ(k) = E X1 {X∈C i (Γ(k)} P(X∈C i (Γ(k))) ,
(ii)Voronoi cell updating: X Γ(k+1) ← X Γ(k+1) .

(3.40)

This procedure clearly lives in Hull supp(µ) by a convexity argument. Note that (i) can re-written (with obvious notations) as

x i (k + 1) = E X | X Γ(k) = x i (k) , i = 1, . . . , N.
As set, this procedure is a pseudo-algorithm since: it cannot be implemented as long as we do not know explicitly the quantities E X1 {X∈C i (Γ(k))} and P X ∈ C i (Γ(k)) at each step k.

⊲ Practical implementation: the randomized Lloyd I procedure.

-When µ is the empirical measure of a dataset, the above procedure can be implemented: one computes the expectation and the probability appearing in the ratio (3.40).

-If the distribution µ has an infinite support (e.g. because µ has a density) the situation becomes more involved. When d = 1, it is often possible to compute the two quantities of interest which appear in the conditional expectation of (3.40)(i) since both are terms of the ratio are 1-dimensional integrals. This is sometimes still possible in low dimension, say d = 2, still using appropriate cubature formulas for numerical integration (see e.g. the algorithms available on the website www.qhull.org and the references therein). Nevertheless, as d increases, the above procedure becomes rapidly intractable: to compute both expectations in (3.40)(i) there is no alternative to the Monte Carlo method (or Quasi-Monte Carlo). One replaces the distribution of X by the empirical measure of a large M -sample of (independent copies) (X m ) 1≤m≤M of X defined a probability space (Ω, A, P), namely

µ = P X ← 1 M M m=1 δ X m (ω),
considering that, for large values of M , both distributions are P(dω)-a.s. get close enough. Then phase (i) (centroid updating) becomes Γ(k + 1) = x i (k + 1, M ), i = 1, . . . , N where

(i) M ≡ x i (k + 1, M ) = M m=1 X m 1 {X m ∈C i (Γ(k))} card {1 ≤ m ≤ M, X m ∈ C i (Γ(k))} , i = 1, . . . , N. (3.41)
This algorithm is often called randomized Lloyd's I algorithm. Like for the CLV Q algorithm, the convergence results for the Lloyd I procedure are still partial, even in its original (deterministic) form (3.40) (see e.g. [START_REF] Pagès | Pointwise convergence of the Lloyd algorithm in higher dimension[END_REF]).

Two examples of (nearly) optimal quadratic quantizations are reproduced in Figures 2 and3. They have been obtained by a hybrid stochastic optimization procedure mixing the randomized Lloyd I algorithm and the CLV Q algorithm. Furthermore, it relies on a "splitting approach" that is a levelby-level computation of the grids (see [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF] as concerns the normal distribution).

For more details about these numerical stochastic optimization procedures used to produce optimal quantizers at level N , we refer e.g. to [START_REF] Benveniste | Algorithmes adaptatifs and approximations stochastiques[END_REF][START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF] for CLV Q, [START_REF] Kieffer | Exponential rate of convergence for Lloyd's method I[END_REF][START_REF] Du | Centroidal Voronoi Tessellations: Applications and Algorithms[END_REF][START_REF] Pagès | Pointwise convergence of the Lloyd algorithm in higher dimension[END_REF] for (randomized) Lloyd's I procedure or more applied textbooks like [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF].

3.4 Nearest neighbor search: how to speed it up. . .

In both above described procedures -CLV Q and randomized Lloyd's I procedures -above, the most time consuming phase is by far the nearest neighbour search which determines the Voronoi cell in which the new input X k+1 falls at step k + 1. This nearest neighbour search is well-known to have an exploding complexity when d and N are large. Its reduction or at least its control is known as a highly challenging problem in computer science.

There are several methods to speed it up, at least in medium dimension. We consider in what follows, either an M -sample (X m (ω)) m=1,...,M obtained by simulation of the distribution µ or a true dataset denoted in both cases by (ξ m ) m=1,...,M to alleviate notations.

⊲ Partial Distance Search principle (P DS, Chen). This idea is quite simple and very efficient: a nearest neighbour search amounts to check, in an Euclidean framework, whether a squared norm 2 is lower than a record value, say δ 2 rec . for N = 1 up to 10 4 and d = 1, . . . , 10. The precise procedure is described in a companion which can be downloaded. On this website are also made available functional quantization grids of the standard Brownian motion over the interval [0, 1], the Brownian bridge, as well as a detailed procedure to compute grids for the (normalized) Ornstein-Uhlenbeck process and its bridge (6 ).

|ξ| 2 = |(ξ 1 , . . . , ξ d )| 2 = (ξ 1 ) 2 + • • • + (ξ d )
4 Greedy quantization 4.1 What is greedy quantization?

We still consider p ∈ (0, +∞) and X : (Ω, A, P) → R d an L p -integrable random vector. The starting idea of optimal greedy quantization is to determine a sequence (a N ) N ≥1 of points of R d which is recursively optimal step by step or level by level with respect to the L p -mean quantization criterion. We mean that, if we denote a (N ) = {a 1 , . . . , a N }, N ≥ 1, then the points a N are recursively defined by a (0) = ∅ and a N +1 ∈ argmin ξ∈R d e p a (N ) ∪ {ξ}, X , N ≥ 0. (4.42)

Note that this section is no longer amide survey but is made of new material, developed in an extended version in [START_REF] Luschgy | Greedy vector quantization, pré-pub[END_REF]. A priori the sequence is not uniquely defined since the above function may attain its minimum at several points of R d . Note that if N = 0, a 1 is simply the L p -median of the distribution µ of X and is subsequently unique whenever p > 1 by a strict convexity argument. This idea to design not only optimal N -tuples but an optimal sequence which, hopefully, will produce N -tuples with a rate optimal behaviour as N → +∞ is very natural and can be compared to sequences with low discrepancy in Quasi-Monte Carlo methods. In fact such sequences have already been investigated but only in an L 1 setting when X has a compact support as a model of short term planning (by contrast with long term planning, see [START_REF] Brancolini | Long-term planning versus short term planning in the asymptotical location problem[END_REF]). Our aim in this section, which is made of new results, is to solve this greedy optimization problem for general distributions µ = P X in any L p -space in two directions: first establish the existence of such L p -optimal greedy sequences and then evaluate the rate of decay of e p (a (N ) , X) to 0 as the quantization level N goes to infinity.

For convenience, we will introduce for every subset A ⊂ R d and every ξ ∈ R d , the notation dist(ξ, A) = inf a∈A |ξ -a| for the distance between ξ and A, where | . | denotes the canonical Euclidean norm in what follows (however, except for algorithmic aspects the results that follow are true with any norm on R d ). B(ξ, r) will denote the closed ball centered at ξ ∈ R d with radius r > 0.

Proposition 4.1. (a) If X ∈ L p R d (P), then the sequence of optimization problems (4.42) admits at least one solution (a N ) N ≥1 where a 1 is the L p -median of the distribution µ. Moreover, the sequence e p (a (n) , X) 1≤n≤N is (strictly) decreasing as long as

N ≤ |supp(µ)| (so that a n / ∈ a (n-1) , 1 ≤ n ≤ N ). (b) Any solution (a N ) N ≥1 to (4.42) satisfies lim N e p (a (N ) , X) = 0 i.e. lim N →+∞ R d min 1≤i≤N |ξ -a i | p µ(dξ) = 0.
Such a solution is called an L p -optimal greedy quantization sequence.

Proof. (a) We proceed by induction. When N = 1, the existence of a 1 obviously follows from the fact that a → E(|X -a| p ) is continuous and goes to infinity as |a| → +∞. Assume there exists a 1 , . . . a N such that e p (a (k) , X) = min a∈R d e p (a (k-1) ∪ {a}, X) for k = 2, . . . , N . If supp(µ) ⊂ {a 1 , . . . , a N } then for every a ∈ R d , e p (a (N ) ∪ {a}, X) = e p (a (N ) , X). Otherwise, let ξ * ∈ supp(µ) \ {a 1 , . . . , a N }. It is clear that |ξξ * | < dist(ξ, a (N ) ) on the ball B ξ * , 1 4 dist(ξ * , a (N ) )

Greedy quantization is rate optimal

In this section we answer the question whether a greedy sequence may be rate optimal i.e. that e p (a (N ) , X) goes to zero at the same rate as a sequence of optimal quantizers at level N when N → +∞. Following [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF], we define for every b ∈ (0, +∞) the b-maximal function associated to a greedy sequence (a N ) N ≥1 by

∀ ξ ∈ R d , Ψ b (ξ) = sup N ≥1 λ d B(ξ, b dist(ξ, a (N ) )) µ B(ξ, b dist(ξ, a (N ) )) ∈ [0, +∞].
It is clear that Ψ b (ξ) > 0 for every ξ = a 1 (L p -median).

Note that this notion of maximal function (originally introduced in [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF]) can be naturally defined with respect to a sequence of grids (Γ N ) N ≥1 where Γ N has size N . The theorem below yields a criterion based on the integrability of the maximal function Ψ b which implies that an L p -optimal greedy quantization sequence is L p -rate optimal (in the sense of Zador's theorem). More practical and easy-to-check criterions are given further on. Proof. First, if µ is a Dirac mass δ a for some a ∈ R d , then a 1 = a and e p (a (N ) , X) = 0 for every integer N ≥ 1. Otherwise, we rely on the following micro-macro inequality established inWe rely on the following micro-macro inequality established in [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF] (see Equation (3.4) in the proof of Theorem 2).

∀ ξ ∈ R d , dist(ξ, a (N ) ) p ≤ C p,b µ B(y, bd(y, a (N ) ))
e p (a (N ) , X) pe p (a (N ) ∪ {ξ}, X) p where b ∈ (0, 1 2 ) and C p,b is a positive real constant depending on p and b. Then, it follows that

e p (a (N ) ∪ {ξ}, X) p ≤ e p (a (N ) , X) p - 1 C p,b µ B(ξ, b dist(ξ, a (N ) )) λ d B(ξ, b dist(ξ, a (N ) )) b d dist(ξ, a (N ) ) p+d V d
where V d denotes the hyper-volume of the unit ball with respect to the current norm on ∈ (0, +∞) is a normalizing real constant). Then, integrating the above inequality with respect to ν yields

R d i.e. V d = λ d B | . | (0; 1) . This implies that e p (a (N ) ∪ {ξ}, X) p ≤ e p (a (N ) , X) p - 1 C p,b,d 1 Ψ b (ξ) dist(ξ, a (N ) ) p+d where C p,b,d = C p,b /(b d V d ) ∈ (0, +∞). Note that µ({a 1 }) < 1 since µ is not a Dirac mass, so that
R d e p a (N ) ∪ {ξ}, X p ν(dξ) ≤ e p (a (N ) , X) p - 1 C p,b,d R d dist(ξ, a (N ) ) p+d ν(dξ) Ψ b (ξ) .
The case of non-compactly supported distributions with radial density functions Lemma 4.2 (see Corollary 3 in [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF]). If X ∈ L p+η R d (P) for some η > 0 with an essentially radial distribution µ = ϕ.λ d in the sense that

ϕ = h(| . | 0 ) on B | . | 0 (0, R) c with h : (R, +∞) → R + , non-increasing and | . | 0 a norm on R d . (4.43) Let (Γ N ) N ≥1 be a sequence of N -quantizers such that R d dist(ξ, Γ N ) p µ(dξ) → 0. If there exists a real constant c > 1 such that R d ϕ(c ξ) -p p+d µ(dξ) = R d ϕ(c ξ) -p p+d ϕ(ξ)dξ < +∞ (4.44) then Ψ b ∈ L p p+d (µ).
In fact Corollary 3 in [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF] is stated to be used only with L p -optimal quantizers so the above formulation includes minor modifications. Combining this lemma with Theorem 4.1 yields the following proposition.

exists R 0 , R ′ 0 ∈ R, R ′ 0 ≥ R 0 such that supp(µ) ⊂ [R 0 , +∞) and f |[R ′ 0 ,+∞) is non-increasing. (4.45)
This criterion is e.g. satisfied by the gamma distributions on R + (including the exponential distributions).

A criterion for non-radial density functions can also be derived from Corollary 4 in [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF], see [START_REF] Luschgy | Greedy vector quantization, pré-pub[END_REF] for details.

From now on in this section, we denote by a N,p N ≥1 an (L p , µ)-optimal greedy quantization sequence f and a (N ),p N ≥1 the resulting sequence of N -quantizers. Further comments and provisional remarks. • If E|X| p+η < +∞ for some η > 0 µ = P X has a non-zero absolutely continuous component ϕ.λ d with respect to the Lebesgue measure, Zador's Theorem implies lim inf

N N 1 d e p (a (N ) , X) ≥ lim N N 1 d e p (X) = J p,d ϕ 1 p L p p+d (λ d ) > 0.
In that setting, we derive that e p (a (N ) , X) ≍ N -1 d as N → +∞. • About the sharpness of optimal greedy sequences. It is pointed out in [START_REF] Brancolini | Long-term planning versus short term planning in the asymptotical location problem[END_REF] (see Theorem 4.10) that, by contrast with what is observed with regular (i.e. global) optimization at level N (Zador's Theorem), N 1 d e p (a (N ) , X) usually does not converge when (a N ) N ≥1 is a greedy L p -optimal sequence. The counter-example is exhibited in the one-dimensional framework for the uniform distribution U ([0, 1]) where an L 1 -optimal greedy sequence (a N,1 ) N ≥ is analyzed. It is shown that for this sequence N e p a (N ),1 , U ([0, 1]) does not converge as N goes to infinity and 1 < lim inf N e 1 a (N ),1 , U ([0,1])

e 1,N U ([0,1]) < lim sup N e 1 a (N ),1 ,U ([0,1]) e 1,N U ([0,1]) < +∞.
Other numerical experiments reproduced below in the quadratic case p = 2, still with U ([0, 1]), on the one hand and, on the other hand, in [START_REF] Luschgy | Greedy vector quantization, pré-pub[END_REF], with the scalar and bivariate normal distributions N (0; 1) and N (0; I 2 ) (which have an unbounded support), strongly suggest that, for more general absolute continuous distributions µ on (R d , Bor(R d )), (L p , µ)-optimal greedy sequence(s) (a N,p ) N ≥1 satisfy lim inf N e p a (N ),p , µ e p,N µ > 1.

However, no proof of this fact is known to us so far. As a consequence, owing to Theorem 2.2 (and the remark that follows), this would prove that no subsequence extracted from an (L p , µ)-greedy optimal sequence (a N,p ) N ≥1 can produce a (subsequence) a (N ′ ),p ) N ≥1 of asymptotically L p -optimal N ′ -quantizers for µ.

• Rate optimality of non-greedy sequences. Another natural question arises at this stage: "Are there rate optimal sequences for the L p -mean quantization error which are not solution to the greedy problem?" To answer -positively -to this question, let us consider the celebrated dyadic Van der Corput (VdC ) sequence, viewed as a quantization sequence rather than a sequence with low discrepancy. Let us recall that the VdC sequence is defined by where J 1,1 = lim N e 1,N U ([0, 1]) . This lim inf is achieved by the subsequence N ′ = 2 n-1 , n ≥ 1, and the lim sup with the subsequence N = 3 2 .2 n = 3.2 n-1 . So we can answer to the first question: it does exist rate optimal sequences for the L 1 -mean quantization error which are not solutions to the greedy problem (4.42). In fact this example shows that it even exists rate optimal sequences (ξ N ) N ≥1 containing subsequence of quantizers (ξ (N ′ ) ) N ≥1 which are asymptotically L 1 -rate optimal quantizers: so is the case of the Van der Corput sequence with the above subsequence N ′ = 2 n-1 .

∀ N ≥ 1, ξ N =
On the other hand he sequence a (N ),1 N ≥1 has been investigated in in [START_REF] Brancolini | Long-term planning versus short term planning in the asymptotical location problem[END_REF], numerical computations carried out in this paper still suggest that lim inf As for the lim inf, one checks again that no subsequence of a (N ),2 N ≥1 can be asymptotically L 2 -optimal.

As for the lim sup, one checks again that the quadratic optimal greedy sequence a N,2 N ≥1 outperforms the VdC sequence from the lim sup criterion since 1.13401 < 1.67706 = 3 √ 5 4 . Both these results and numerical experiments naturally lead to two open theoretical questions, left for further investigations:

(1) Does it exist distributions µ (with a non-countable support) for which an L p -optimal greedy sequence (a N,p ) N ≥1 produces a sequence of asymptotically optimal N -quantizers(a (N ),p ) N ≥1 i.e. such that ep(a (N ),p ,µ) e p,N (µ)

→ 1 as N goes to +∞; or at least subsequences (a (N ′ ),p ) N ≥1 having this property ?

(2) Does an (L p , µ)-optimal greedy sequence for a distribution µ on R d produce the lowest possible value for lim sup N N 1 d e p,N a (N ),p , µ among all R d -valued sequences (ξ N ) N ≥1 ?

The answer to the first question is probably negative (based on our numerical experiments) whereas the second one remains more uncertain.

Algorithmic aspects

We adopt notations of Sections 3. Practical computation of an optimal greedy sequence of quantizers relies on obvious variants of the historical algorithms (CLV Q and Lloyd I) implemented recursively: to switch from level N to N + 1, one first adds a (N + 1) th point (sampled from the support of the distribution µ) to the N -tuple (a 1 , . . . , a N ) computed during the first N stages of the optimization procedure. This makes up the starting (N + 1)-tuple for the modified CLV Q or Lloyd I procedure. Then, one implements one of these two optimization procedures with the following restriction: all formerly computed components a i , 1 ≤ i ≤ N , are frozen, and only the new point is moved, following the standard rules. Thus, if implementing a greedy CLV Q algorithm, when the (N + 1) th component is the "winner" in the competition phase (i.e. this (N + 1) th component at the m th is the nearest neighbour to the new input stimulus, say X m+1 ). For Lloyd's I procedure, this (N + 1) th component is the only one to be updated following (3.41), the other N components remaining frozen as well.

For more details about these greedy variants, especially Lloyd's I, we refer to [START_REF] Luschgy | Greedy vector quantization, pré-pub[END_REF] where convergence proofs are provide (compete in 1-diemsnion, partial in hgher dimension).

We reproduce in Figure 4 the graph of the N → N e 2,N a (N ) , U ([0, 1]) where (a N ) N ≥1 is an L 2 -optimal greedy quantization sequence for the uniform distribution U ([0, 1]).

Greedy quantization versus Quasi-Monte Carlo?

Of course, for every integer N ≥ 1, the weights induced by the µ-mass of the Voronoi cells associated to a (N ) define canonically a sequence of N -tuples which usually cannot be "arranged" into a sequence, even up to a re-scaling. When considering the unit hypercube [0, 1] d as a state spaces in d dimension, it is easy natural to compare an optimal greedy sequence with respect to the uniform distribution U ([0, 1] d ) and the so-calleduniformly distributed sequences usually implemented in the Quasi-Monte Carlo method.

Let us recall that a sequence (ξ N ) N ≥1 is uniformly distributed over [0, 1] d if the empirical measures f (u)du. This means that the weights associated to a uniformly distributed sequence are by definition equal to 1 N (i.e. are equal to 1 up to the normalization factor 1/N ). We will see that the cost induced by considering these uniform weights 1 N is essentially log N . This follows from Proinov's Theorem (see [START_REF] Proinov | Discrepancy and integration of continuous functions[END_REF]) recalled below, which evaluates precisely the rate of convergence of these sequences on Lipschitz continuous functions.

In the Quasi-Monte Carlo (QM C) method, the performance of an N -tuple (ξ 1 , . . . where the bases p i , i = 1, . . . , d are the first d prime numbers), the Faure sequences, the Sobol' sequences and the unifying framework developed by Niederreiter [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]. For definitions of these sequences and numerical tests on various numerical integration problems we refer to [START_REF] Bouleau | Numerical methods for stochastic processes[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF][START_REF] Pagès | Sequences with low discrepancy and pseudo-random numbers : theoretical results and numerical tests[END_REF][START_REF] Pagès | Introduction to Numerical Probability and Applications to Finance[END_REF]. Although such a rate has never been proved to be the lowest possible, its optimality is a commonly shared opinion in the QM C community (however see [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] for a review of existing lower bounds).

The striking fact with these sequences satisfying (4.48), called sequences with low discrepancy, is that, when implemented on the class of functions with finite variation on [0, 1] d , the Koksma-Hlawka inequality implies that, for every such function f :

[0, 1] d → R [0,1] d f dλ d - 1 N N i=1 f (ξ i ) ≤ V (f )D * N (ξ 1 . . . , ξ N )
where V (f ) denotes the variation of the function f . So it induces for this specific class of functions a rate of numerical integration of order O (log N ) d

N

. In one dimension (d = 1), the above notion of finite variation coincides with the standard definition of finite variation in real analysis. When d ≥ 2, several definitions can be given, the most popular being the finite variation in the Hardy & Krause sense (as described e.g. in [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]). Another -slightly less general but more elementary -being the finite variation in the signed measure sense developed in [START_REF] Bouleau | Numerical methods for stochastic processes[END_REF] (see also [START_REF] Pagès | Introduction to Numerical Probability and Applications to Finance[END_REF]). Unfortunately, as the dimension d increases, the set of functions with finite variation (in any of the above senses) becomes somewhat "sparse" among the set of all real-valued Borel functions defined on [0, 1] d . So, one may have doubts about this striking performance when dealing with practical simulation problems. This is confirmed by their behaviour on the more natural space of Lipschitz continuous functions which is ruled by the following Proinov Theorem.

Figure 1 :

 1 Figure 1: Quantization of a scalar (periodic) signal (B. Wilbertz)

t 0 δ

 0 x(s) ds weakly converges as t → +∞ toward a distribution µ on (R d , Bor(R d )) such that R d |ξ| p µ(dξ) < +∞. If the quantization function q is µ-a.s. continuous and, e.g., lim sup t→+∞ 1 t t 0 |x(s)| r ds < +∞ for some r > p, )q(x(s))| p ds 1 p = R d |ξq(ξ)| p µ(dξ) 1 p

  1≤i≤N |Xx i | L p (P) (1.1) = min 1≤i≤N |ξx i | L p (µ) = R d min 1≤i≤N |ξx i | p µ(dξ)

  R d |ξ| p µ(dξ) < +∞ so that the distortion function G p,N is finite everywhere on (R d ) N . (a) The distortion function G p,N attains a minimum at an N -tuple x (N,p) = (x (N,p) 1 , . . . , x (N,p) N ). (b) If card supp(µ) ≥ N , then the corresponding grid Γ (N,p) = x (N,p) 1 , . . . , x (N,p) N has full size N and for every Voronoi partition C i (Γ (N ) ) 1≤i≤N of R d induced by Γ (N ) , P(X ∈ C i (Γ (N ) ) > 0. (c) The sequence N → e p,N (X) (strictly) decreases as long as N ≤ |supp(µ)| and lim N e p,N (X) = 0.

1. 3 . 1

 31 Data mining, clustering, automatic classificationLet (ξ k ) 1≤k≤n be an R d -valued dataset and µ the uniform distribution over this dataset -the empirical measure of the dataset -defined by

=

  denotes the identity in distribution.Proposition 1.3. Let X ∈ L p R d (P), p ∈ (0, 1], with distribution µ and let Γ = {x 1 , . . . , x N }. (a) For every p ∈ (0, 1], X -X Γ p L p (P) = sup [F ]p≤1 |EF (X) -EF X Γwhere [F ] p = sup x =y |F (x)-F (y)| |x-y| p denotes the p-Hölder coefficient of the function F : R d → R. (b) If P N denotes the set of probability measures with a support having at most N points in R d , then W 1 (µ, P N ) = e 1,N (µ).

Theorem 1 . 2 (

 12 Feynman-Kac's representation formula). Assume (for simplicity) that the functions b and a are such that the above PDE (1.16) has a unique C 1,2 ([0, T ] × R d ) solution u whose gradient ∇ x u has polynomial growth in x, uniformly in t ∈ [0, T ]. Let σ : R d → M(d, q) ( 3 ) such that a = σσ * (where * stands for matrix transposition). Assume that b and σ are continuous on [0, T ] × R d and, at least, Lipschitz continuous in x, uniformly in t ∈ [0, T ].

  in particular a {k, . . . , n}-valued F-stopping time. It implies the inequality in the last line owing to the definition (1.20) of the Snell envelope. Since U k ≥ Z k , by considering the deterministic stopping time τ = k, we finally get

t n k )

 k 0≤k≤n -Snell envelope of the payoff process f k ( Xn t n k ) 0≤k≤n and ( U k ) 0≤k≤n the quantized pseudo-Backward Dynamic Principle associated to a quantized version of this payoff process ( Xn t n k ) 0≤k≤n as defined by (1.23). Applying the control established in claim (b) of the above theorem yields with obvious notations

i

  , i = 1, . . . , N }. In the second case, let Γ * = Γ (N ) ∪{ξ N +1 } and let (C * i (Γ * )) 1≤i≤N +1 be an induced Voronoi partition where• CN+1 (Γ * ) is the Voronoi cell of ξ N +1 . As ξ N +1 / ∈ Γ (N ) , it is clear that • CN+1 (Γ * ) = ∅ and that |Xξ N +1 | < min 1≤i≤N |Xx (N ) i |on the interior of this cell. Furthermore, P(X ∈ C(Γ * ) N +1 ) = µ • CN+1 (Γ * ) > 0 since ξ N +1 ∈ • CN+1 (Γ * ) and ξ N +1 lies in the support of µ. Note that everywhere on (R d ) N , one has |Xξ N +1 | ∧ min 1≤i≤N |Xx i | ≤ min 1≤i≤N |Xx i |, so that, combining both inequalities yields

Definition 2 . 1 .

 21 A (quantization) grid Γ = {x 1 , . . . , x N } ⊂ R d of size N is µ-stationary if the following two conditions are satisfied  µ-negligibility of the boundary of the Voronoi diagram : µ 1≤i≤N ∂C i (Γ) = 0.

lim n N 1 d 1 d

 11 e p,N (X) = J p,d ϕ 1 p L p p+d (λ d ) where J p,d = inf N ≥1 N e p,N U ([0, 1] d ) ∈ (0, +∞) (depends on the current norm). When this norm is the canonical Euclidean norm on R d , one has J 2,2 = 5 18 √

d

  where C d,µ is a real constant only depending on the dimension d and the ((2 + η)-pseudo standard deviation of) distribution µ given by (2.35) in Theorem 2.1 (b).

Figure 2 :

 2 Figure 2: Optimal N -quantizer (N = 500) of the bi-variate Normal distribution depicted with its Voronoi tessellation (with J. Printems).

Figure 3 :

 3 Figure 3: Optimal N -quantization (N = 500) of W 1 , sup t∈[0,1] W t depicted with its Voronoi tessellation, W standard Brownian motion (with B. Wilbertz).

Theorem 4 . 1 . 1 d

 411 Let p ∈ (0, +∞) and let µ = P X be such thatR d |ξ| p µ(dξ) < +∞. Assume that there exists b ∈ (0, 1 2 ) such that Ψ b ∈ L p p+d (µ). Then lim supN N e p (a (N ) , X) < +∞.

  Ψ b ∈ L p p+d (µ), we can define the probability distribution ν = κ b,p,d Ψ

Proposition 4 . 3 ( 1 d

 431 Non-compact support with radial density). If X ∈ L p+η R d (P) for some η > 0 with an essentially radial distribution in the sense of (4.43). If, furthermore, ϕ satisfies (4.44), then any greedy L p -optimal sequence (a N ) N ≥1 is rate optimal i.e. satisfies lim sup N N e p (X, a (N ) ) < +∞. This case includes e.g. all the (centered) hyper-exponential distributions of the form µ = ϕ.λ d with ϕ(ξ) = κ a,b,c |ξ| c 0 e -a|ξ| b 0 , ξ ∈ R d , a, b > 0, c > -d and | . | 0 denotes any norm on R d and subsequently all hyper-exponential distributions distributions since L p -mean-quantization errors is invariant by translation of the random vector X. In particular, this includes all normal and Laplace distributions Remark. In one dimension, (4.43) can be replaced mutatis mutandis by a one-sided variant: if there

r k=0 n k 2 4 = J 1 , 1

 2411 k+1 where N = n r 2 r + • • • + n 0 , n i ∈ {0, 1}, i = 1, . . . , r. (4.46) ⊲ L 1 -mean quantization problem. Elementary computations carried out with the L 1 -mean quantization error modulus, not reproduced here, show that lim inf N N e 1 (ξ 1 , . . . , ξ N , [0, 1]) = 1 and lim sup N N e 1 (ξ 1 , . . . , ξ N , [0

N N e 1 √ 3 = J 2 , 1 √ 5 4 × J 2 , 1 where J 2 , 1

 132142121 (a (N ),1 , U ([0, 1])) ≈ 1, 02 × J 1,1 and lim supN N e 1 (a (N ),1 , U ([0, 1])) ≈ 1, 09 × J 1,1The second inequality the sequence a N,1 N ≥1 outperforms the VdC sequence from the lim sup criterion since 1.09 < 1.125 = 9/8. ⊲ L 2 -mean quantization problem. The same phenomenons are confirmed in the quadratic case since lim inf N N e 2 (ξ 1 , . . . , ξ N , [0, 1]) = 1 2 and lim sup N N e 2 (ξ 1 , . . . , ξ N , [0, 1]) = 3 = lim N N e 2,N (U ([0, 1])) = inf N N e 2,N (U ([0, 1])). Using the greedy Lloyd I procedure briefly described in the next Section 4.3, we observe numerically that lim inf N N e 2 (a (N ),2 , U ([0, 1])) ≈ 1.02732 × J 2,1 > J 2,1 and lim sup N N e 2 (a (N ),2 , U ([0, 1])) ≈ 1.13401 × J 2,1 .

Figure 4 :

 4 Figure 4: N → N e 2 a (N ),2 , U ([0, 1]) , N = 1, . . . , 10 000, (a N,2 ) N ≥1 L 2 -optimal greedy sequence; flat solid line (---) ≡ J 2,1 =

δ

  ξ i , N ≥ 1, weakly converges toward the Lebesgue measure λ d on [0, 1] d . In particular this means that, for every bounded λ d -a.s. continuous function f : [0, 1] d → R,

1

  , ξ N ) ∈ ([0, 1] d ) N is measured by the Kolmogorov-Smirnov distance between the d-dimensional cumulative distribution function of its empirical measure ν N and the uniform distribution U ([0, 1] d ), namely the so-called star discrepancy defined byD * N (ξ 1 . . . , ξ N ) = sup {ξ i ∈[[0,u]]}λ d [[0, u]] (4.47)where[[0, u]] = d ℓ=1 [0, u ℓ ], u = (u 1 , . . . , u d ) ∈ [0, 1] d .Several sequences ξ = (ξ N ) N ≥1 have been exhibited (see[START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]) whose star discrepancy at the origin satisfy for a real constant C(ξ) ∈ (0, +∞),∀ N ≥ 1, D * N (ξ 1 . . . , ξ N ) ≤ C(ξ) (1 + log N ) d N . (4.48) Among them, in 1-dimension, the dyadic VdC sequence as defined by (4.46) (and its p-adic generalizations VdC (p), p ≥ 2) and, in d-dimension, the Halton sequences (made up with VdC (p i ) sequences

Theorem 4 . 2 .

 42 (Proinov,[START_REF] Proinov | Discrepancy and integration of continuous functions[END_REF]) Assume R d is equipped with the ℓ ∞ -norm |(ξ 1 , . . . , ξ d )| ∞ = max 1≤i≤d |ξ i |. For every continuous function f : [0, 1] d → R, let w(f, δ) := sup ξ, ξ ′ ∈[0,1] d , |ξ-ξ ′ |∞≤δ |f (ξ)f (ξ ′ )|, δ ∈ (0, 1).denote the related uniform continuity modulus of f (with range δ). (a) Let (ξ 1 , . . . , ξ N ) ∈ ([0, 1] d ) N . For every continuous function f : [0, 1] d → R and every integer N ≥ 1, [0,1] d f dλ d -1 N N i=1 f (ξ i ) ≤ C d w f, D * N (ξ 1 , . . . , ξ N )

More generally, if C is X (N) -valued, then X is called a card(X )-alphabet.

The sender A only needs a codebook to discriminate the elementary quantizers xi i.e. a codebook where all xi are known with a fixed length ℓ1 (dyadic) bits in its dyadic representation. The receiver B may need arbitrary accurate values for the elementary quantizers xi in his/her codebook in order to reconstruct a posteriori the signal.

matrices with d rows and q columns.

When d = 1 an exact (and efficient) simulation method has been devised for diffusion processes in[START_REF] Beskos | Exact simulation of diffusions[END_REF]; unfortunately it deeply depends on the scalar feature of the diffusion.

where a k ∼ b k means that a k = u k b k with lim k u k = 1.

All downloads for scientific and non-commercial purposes are free of charges.

d

Acknowledgement:

The author thanks B. Jourdain and the referee for their careful reading of the manuscript. All errors are mine.

To test the inequality |ξ| 2 < δ 2 rec , one proceeds by induction on the coordinates as described in Algorithm 1 below.

Algorithm 1 (Chen's Partial Distance Search protocol) record = true; aux = 0; ℓ = 1 while aux + (ξ ℓ ) 2 < δ 2 rec do ℓ := ℓ + 1; end while if ℓ = d + 1 then record = f alse; end if

The testing procedure is stopped as soon as the partial computation of the squared norm is higher than the record value. This procedure is known as Chen's Partial Distance Search protocol. Adaptation to other common norms is straightforward. It can still be improved in medium dimensions, say up to d = 10, by storing the data in a tree structure as described below.

⊲ Space partitioning tree (see [START_REF] Friedman | An Algorithm for Finding Best Matches in Logarithmic Expected Time[END_REF]). As a first step, one divides the dataset into two subsets of size M 2 according to the median of the first coordinate. Then, one divides again each sub-sample following the second coordinate and so on. At each step, the procedure is applied either to a new coordinate or, when the dimension is exhausted, to the one with the highest index. Each leaf of the resulting tree is labelled. The data are coded in the tree by their leaf label and a local index inside its assigned leaf. After this pre-processing phase, whose complexity is O(M log M ), the search cost of a single datum in the tree is O(log M ). It is commonly shared by users that the depth of the tree should be fixed to 7 to get the best compromise. In many situations of interest, this choice seems more or less independent of the size of the dataset and its structure.

⊲ Principal Axis Tree (McNames, see [START_REF] Mcnames | A Fast Nearest-Neighbor Algorithm Based on a Principal Axis Search Tree[END_REF]). The additional improvement brought by this approach lies in a preliminary Principal Component Analysis (P CA) inducing a change of coordinates which make the search more efficient, especially in view of applying the P DS principle: in the new basis induced by the P CA the absolute value of coordinates tends to decrease when the coordinate index increases implying that Chen's P DS procedure stops sooner in average. ⊲ Rough Self-Quantization (Corlay, see [START_REF] Corlay | Some aspects of optimal quantization and applications to finance[END_REF], chapter 1). The underlying idea is completely different and competes with the P AT approach when designing the search tree. It consists of a cascade of rough pre-quantizations of the dataset with N 0 (N 0 ≈ 7) prototypes (or classes) at each level. This again speeds up Chen's P DS principle of course.

A different approach, not detailed here, is to devise approximate nearest neighbor search procedures, often based on some anchor points (see [START_REF] Gersho | Vector Quantization and Signal Compression[END_REF] for examples and the references therein).

Where to download optimal quantization grids?

Optimized quadratic quantization grids of the d-dimensional normal distributions N (0; I d ), computed on the occasions of various numerical experiments (pricing of multi-asset American style option, portfolio management, nonlinear filtering, swing option pricing, etc) by the hybrid "splitting method" described above can be downloaded from the website www.quantize.maths-fi.com which satisfies µ B ξ * , 1 4 dist(ξ * , a (N ) ) > 0. Consequently, e p (a (N ) ∪ {ξ * }, X) < e p (a (N ) , X). Now let

This yields a contradiction which in turn implies that K 0 N +1 is a compact set. On the other hand, ξ → e p (a (N ) ∪ {ξ}, X) is clearly Lipschitz continuous on R d , hence it attains its minimum on K 0 N +1 which is clearly its absolute minimum.

(b) It is clear that, for every ξ ∈ R d , min 1≤i≤N |ξ-a i | is non-increasing and converges toward inf

) .

Let N 0 be a positive integer such that,

We consider the (N 0 + 1)-quantizer a (N 0 ) ∪ {ξ 0 }. On the one hand,

and, on the other hand,

which yields a contradiction. Hence ℓ ∞ = 0 which completes the proof.

Jensen's Inequality applied to the convex function u → u

On the other hand, it is clear that

p so that, finally, if we set A N = e p (a (N ) , X) p , N ≥ 1, this sequence satisfies for every integer N ≥ 1, the recursive inequality

where κ b,p,d = κ Remark. When µ has no absolutely continuous part with respect to the Lebesgue measure, it is likely that, like for standard optimal vector quantization in Zador's Theorem, this rate is not optimal. The natural conjecture should be that greedy quantization sequence(s) go to 0 at the same rate as that obtained for sequences of optimal quantizers which is not N -1 d when the distribution µ is singular (see e.g. [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF]).

We produce below easy-to-check criterions that can be straightforwardly borrowed (and slightly adapted) from [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF] which ensure that Ψ b has the required integrability property so that any greedy sequence is rate optimal.

The case of compactly supported distributions. Lemma 4.1 (see Lemma 1 in [START_REF] Graf | Distortion mismatch in the quantization of probability measures[END_REF]). If X ∈ L p R d (P) has a distribution µ having an absolutely continuous component and if (Γ N ) N ≥1 is a sequence of N -quantizers such that R d dist(ξ, Γ N ) p µ(dξ) → 0 then the maximal functions Ψ b lie in all L r loc (µ) for every r ∈ (0, 1) i.e.

∀ r ∈ (0, 1), ∀ b ∈ (0, +∞), ∀R ∈ (0, +∞),

By combining this result (with r = p p+d ) with Proposition 4.1, we derive the following proposition which generalizes the one stated in [START_REF] Brancolini | Long-term planning versus short term planning in the asymptotical location problem[END_REF] for absolutely continuous distributions µ having a convex support on R d . where C d ∈ (0, +∞) is a universal optimal real constant only depending on d. In particular, if the function f :

(b) In particular if (ξ N ) N ≥1 is a sequence with low discrepancy in the above sense, then for every integer N ≥ 1, The above claim (b) emphasizes the fact that considering uniform weights 1 N induces the loss of a log N factor compared to an (L 1 , U ([0, 1]))-optimal greedy (or simply rate optimal) sequence since, for such a greedy sequence (a N ) N ≥1 , one has, for every N ≥ 1, ). However, by contrast, one checks that the cubature formulas based on optimal quantization turn out to be efficient for much lower values of N than sequences with low discrepancy (see e.g. the numerical experiments carried out in [START_REF] Pagès | Optimal quadratic quantization for numerics: the Gaussian case[END_REF] for the pricing of European derivatives).