5
$\begingroup$

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

$\endgroup$

1 Answer 1

7
$\begingroup$

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

$\endgroup$
6
  • $\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

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .