# Quantization of normal distribution

For $$n\in\mathbb{N}$$, denote by $$\mathcal{Q}_n$$ the set of all probability measures on $$\mathbb{R}$$ that are supported on at most $$n$$ points.

Question: Is it known which element in $$\mathcal{Q}_n$$ is closest to the standard normal distribution with respect to the $$p$$-Wasserstein distance (for some $$p \geq 1$$)? I.e., if $$\mathcal{N}$$ is the standard normal distribution, can we calculate $$\arg\min_{\mu \in \mathcal{Q}_n} W_p(\mu, \mathcal{N})$$?

Remarks: In particular, I am interested in specific optimizers for small $$n$$ (around $$n=10$$). If there are provably reliable numerical methods that can calculate the above quantizers, that would also answer my question.

I found some papers in this direction, but they all appeared to look for quantizers of the form $$\mu = \frac{1}{n}\sum_{i=1}^n \delta_{x_i}$$, while I am really interested in the general case $$\mu = \sum_{i=1}^n w_i \delta_{x_i}$$ for $$0 \leq w_i \leq 1$$, $$\sum_{i=1}^n w_i = 1$$.

$$\newcommand{\N}{\mathcal N}\newcommand{\vpi}{\varphi}$$For $$\mu=\sum_{i=1}^n w_i\delta_{x_i}$$ with $$0\le w_i\le 1$$ and $$\sum_{i=1}^n w_i=1$$, we have the following (see e.g. page 4): \begin{align*} W_p(\mu,\N)^p&=\int_0^1 du\,|F^{-1}(u)-G^{-1}(u)|^p \\ &=S(x,u):=\sum_{j=1}^n\int_{u_j}^{u_{j+1}} du\,|x_j-g(u)|^p, \end{align*} where $$F$$ is the cdf of $$\mu$$, $$G$$ is the cdf of $$\N$$, $$F^{-1}$$ is the generalized inverse of $$F$$, $$g:=G^{-1}$$ is the inverse of $$G$$, $$x:=(x_1,\dots,x_n)$$, $$-\infty, $$u:=(u_1,\dots,u_n)$$, and $$\begin{equation} u_j:=\sum_{i=1}^{j-1} w_i, \tag{0} \end{equation}$$ so that $$\begin{equation} u_1=0\le u_2\le\cdots\le u_n\le u_{n+1}=1. \tag{1} \end{equation}$$

We want to minimize $$S(x,u)$$ in $$x,u$$. By continuity and compactness, for each $$x$$, the minimum of $$S(x,u)$$ in $$u$$ is attained. Let $$u(x)$$ denote a corresponding minimizer. Then without loss of generality the inequalities in (1) are strict (otherwise, one of the corresponding $$w_i$$'s is $$0$$, and so, the cardinality of the support set of $$\mu=\sum_{i=1}^n w_i\delta_{x_i}$$ can be reduced).

For each $$j\in\{2,\dots,n-1\}$$, the partial derivative of $$S(x,u)$$ in $$u_j$$ is $$|x_{j-1}-g(u_j)|^p-|x_j-g(u_j)|^p=0$$ at $$u=u(x)$$, whence $$g(u_j)=(x_{j-1}+x_j)/2$$ and $$\begin{equation} u_j=G\big((x_{j-1}+x_j)/2\big). \tag{2} \end{equation}$$ So, with the substitution $$t=g(u)\iff u=G(t)$$, we get \begin{align*} S_*(x):=\min_u S(x,u)=S(x,u(x))\\ =\sum_{j=1}^n\int_{(x_{j-1}+x_j)/2}^{(x_{j+1}+x_j)/2} dt\,\vpi(t)|x_j-t|^p, \tag{3} \end{align*} where $$\vpi$$ is the pdf of $$\N$$, $$(x_{j-1}+x_j)/2:=-\infty$$ for $$j=1$$ and $$(x_{j+1}+x_j)/2:=\infty$$ for $$j=n$$. If $$p$$ is a natural number, then the integrals in (3) can be expressed in terms of $$G$$ and $$\vpi$$.

So, the minimum of $$W_p(\mu,\N)^p$$ is $$\begin{equation*} \min_{x,u} S(x,u)=\min_x S_*(x), \end{equation*}$$ and the latter minimum can be found for $$n\le10$$ by such provably reliable numerical methods as the interval arithmetic method.

For $$p=1$$ and $$n=10$$, the best values of $$x_1,\dots,x_n$$ found by Mathematica are $$-1.94975, -1.31679, -0.87372, -0.503282, -0.164679, \\ 0.164678, 0.503281, 0.873719, 1.31679, 1.94975.$$ The corresponding optimal $$w_i$$'s can then be found by (0) and (2). Here are the graphs $$\{(t,F(t))\colon|t|<2.5\}$$ (gold) and $$\{(t,G(t))\colon|t|<2.5\}$$ (blue): • Thanks for the answer! I don't quite understand how the interval arithmetic method is used to solve problem (3). As I understand it, additional to the interval arithmetic method there has to be some optimization procedure or optimality criteria that is applied? I don't yet see what the important structure of $x \mapsto S_*(x)$ is which is utilized to solve (3). (For instance, would any other density $\varphi$ be valid as well?). Could you perhaps elaborate on that? Apr 12, 2021 at 6:52
• @Steve : I think the main benefit of this answer is the reduction of the number of variables from $2n-1$ to $n$. As for the structure of $S_*(x)$, note that its partial derivatives in $x_j$ depend only on at most three variables: $x_ j,x_{j-1},x_{j+1}$. So, these partial derivatives can be bounded comparatively easily using the interval arithmetic. This should help quite a bit. I have not done such specific calculations -- which is of course a lot of work, but which I think can be done. Apr 12, 2021 at 17:26
• Previous comment continued: If, instead of $\varphi$, you have some other density, then the calculations could be more difficult, depending of course on the density. Apr 12, 2021 at 17:26
• Thanks for the explanation. I definitely see the benefit of the reformulation. Do I understand correctly that to obtain the concrete values for $n=10$ you posted (so, what mathematica is doing), one calculates numerically where the derivative of $S_*$ is zero? I am wondering whether one can show that $S_*$ really has a unique global minimum. $S_*$ is not convex as far as I can tell? Apr 12, 2021 at 17:54
• @Steve : No, I did not mean finding numerically zeroes of derivatives of $S_*$. What I meant is this: According to the interval arithmetic method, partition $\mathbb R$ into, say, $k$ subintervals and then accordingly partition $\mathbb R^n$ into $k^n$ $n$-dimensional boxes. Using the partial derivatives (possibly of higher orders), bound $S_*$ from below on each box. The lower bounds for most of the boxes will be greater then the quasi-minimum of $S_*$ found by Mathematica (say). Work similarly with each of the remaining boxes, etc. Apr 12, 2021 at 19:09