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$.


1 Answer 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<x_1\le\dots\le 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):

enter image description here

  • $\begingroup$ 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? $\endgroup$
    – Steve
    Commented Apr 12, 2021 at 6:52
  • $\begingroup$ @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. $\endgroup$ Commented Apr 12, 2021 at 17:26
  • $\begingroup$ 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. $\endgroup$ Commented Apr 12, 2021 at 17:26
  • $\begingroup$ 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? $\endgroup$
    – Steve
    Commented Apr 12, 2021 at 17:54
  • $\begingroup$ @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. $\endgroup$ Commented Apr 12, 2021 at 19:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.