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$.
 A: $\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):

