# Efficient Algorithms for Sparse Moment Problems without Separation

Zhiyuan Fan<EMAIL_ADDRESS>IIIS, Tsinghua University  Jian Li
<EMAIL_ADDRESS>IIIS, Tsinghua University

###### Abstract

We consider the sparse moment problem of learning a $k$-spike mixture in high
dimensional space from its noisy moment information in any dimension. We
measure the accuracy of the learned mixtures using transportation distance.
Previous algorithms either assume certain separation assumptions, use more
recovery moments, or run in (super) exponential time. Our algorithm for the
1-dimension problem (also called the sparse Hausdorff moment problem) is a
robust version of the classic Prony’s method, and our contribution mainly lies
in the analysis. We adopt a global and much tighter analysis than previous
work (which analyzes the perturbation of the intermediate results of Prony’s
method). A useful technical ingredient is a connection between the linear
system defined by the Vandermonde matrix and the Schur polynomial, which
allows us to provide tight perturbation bound independent of the separation
and may be useful in other contexts. To tackle the high dimensional problem,
we first solve the 2-dimensional problem by extending the 1-dimension
algorithm and analysis to complex numbers. Our algorithm for the high
dimensional case determines the coordinates of each spike by aligning a 1-d
projection of the mixture to a random vector and a set of 2d-projections of
the mixture. Our results have applications to learning topic models and
Gaussian mixtures, implying improved sample complexity results or running time
over prior work.

##  1 Introduction

###  1.1 Background

We study the problem of learning a statistical mixture of $k$ discrete
distributions over a common discrete domain $[d]=\\{1,2,\cdots,d\\}$ from
noisy measurements of its first few moments. In particular, we consider
mixture $\boldsymbol{\vartheta}$ which is a $k$-spike distribution supported
on the $d-1$ dimensional simplex
$\Delta_{d-1}=\\{\boldsymbol{x}=(x_{1},\cdots,x_{d})\in\mathbb{R}^{d}\mid\sum_{i=1}^{d}x_{i}=1,x_{i}\geq
0\,\,\forall i\in[d]\\}$. Each point in $\Delta_{d-1}$ represents a discrete
distribution over $[d]$.

The sparse Hausdorff moment problem: We first define the sparse moment problem
for $d=2$, where the mixture is a $k$-spike distribution supported on
$\Delta_{1}\cong[0,1]$. We call the model the $k$-coin model (i.e., a mixture
of $k$ Bernoulli distributions). Denote the underlying mixture as
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$. Here,
$\boldsymbol{\alpha}=\\{\alpha_{1},\alpha_{2},\cdots,\alpha_{k}\\}$ and
$\boldsymbol{w}=\\{w_{1},w_{2},\cdots,w_{k}\\}\in\Delta_{k-1}$, where
$\alpha_{i}\in[0,1]$ specifies the $i$th Bernoulli distribution and
$w_{i}\in[0,1]$ is the probability of mixture component $\alpha_{i}$. For
mixture $\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$, the
$t$-th moment is defined as

$M_{t}(\boldsymbol{\vartheta})=\int_{[0,1]}\alpha^{t}\boldsymbol{\vartheta}(\mathrm{d}\alpha)=\sum_{i=1}^{k}w_{i}\alpha_{i}^{t}.$

Our goal is to learn the unknown parameters
$(\boldsymbol{\alpha},\boldsymbol{w})$ of the mixture $\boldsymbol{\vartheta}$
given the first $K$ noisy moment values $M^{\prime}_{t}$ with
$|M^{\prime}_{t}-M_{t}(\boldsymbol{\vartheta})|_{\infty}\leq\xi$ for $1\leq
t\leq K$. The problem is also known as the sparse Hausdorff moment problem or
the $k$-coin model [57, 45, 31]. We call $K$ the moment number and $\xi$ the
moment precision/accuracy, both of them are important parameters of the
problem. We measure the quality of our estimation
$\widetilde{\boldsymbol{\vartheta}}$ in term of the transportation distance
between probability distributions, i.e.,
$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(\epsilon)$. Using transportation distance as the metric is advantageous for
several reasons: (1) if the desired accuracy $\epsilon$ is much smaller than
the minimum separation $\zeta=\min\|\alpha_{i}-\alpha_{j}\|$, the estimation
$\widetilde{\boldsymbol{\vartheta}}$ must be a per-spike recovery since it
must contain a spike that is sufficiently close to an original spike in
$\boldsymbol{\vartheta}$ (given the weight of the spike is lower bounded); (2)
if the desired accuracy $\epsilon$ is larger than $\zeta$ or the minimum
weight $w_{\min}$, we are allowed to be confused by two spikes that are very
close, or miss a spike with very small weight, thus potentially avoiding the
inverse dependency on $\zeta$ and $w_{\min}$, which are otherwise unavoidable
if we must recover every spike. Generally, for this problem, we are interested
in understanding the relations and trade-offs among the moment number $K$, the
moment accuracy/precision $\xi$, and the accuracy $\epsilon$ (i.e.,
transportation distance at most $\epsilon$) of our estimation
$\widetilde{\boldsymbol{\vartheta}}$, and designing efficient algorithms for
recovering $\boldsymbol{\vartheta}$ to the desired accuracy.

The sparse moment problem in higher dimensions: For the higher dimensional
case, suppose the underlying mixture is
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$ where
$\boldsymbol{\alpha}=\\{\boldsymbol{\alpha}_{1},\ldots,\boldsymbol{\alpha}_{k}\\}$
($\boldsymbol{\alpha}_{i}\in\Delta_{d-1}$ which are $d$-dimensional points in
a bounded domain) and
$\boldsymbol{w}=\\{w_{1},\ldots,w_{k}\\}\in\Delta_{k-1}$. For
$\boldsymbol{\alpha}_{i}=(\alpha_{i,1},\ldots,\alpha_{i,d})$ and
$\mathbf{t}=(t_{1},\cdots,t_{d})\in\mathbb{Z}_{+}^{d}$, we denote
$\boldsymbol{\alpha}_{i}^{\mathbf{t}}=\alpha_{i,1}^{t_{1}}\cdots\alpha_{i,d}^{t_{d}}$,
and the $\mathbf{t}$-moment of $\boldsymbol{\vartheta}$ is

$M_{\mathbf{t}}(\boldsymbol{\vartheta})=\int_{\Delta_{d-1}}\boldsymbol{\alpha}^{\mathbf{t}}\boldsymbol{\vartheta}(\mathrm{d}\boldsymbol{\alpha})=\sum_{i=1}^{k}w_{i}\boldsymbol{\alpha}_{i}^{\mathbf{t}}.$

We are given noisy access of $\mathbf{t}$-moment
$M_{\mathbf{t}}(\boldsymbol{\vartheta})$ for $\|\mathbf{t}\|_{1}\leq K$. Here
we also call $K$ the moment number. Since there are $K^{O(d)}$ such moments,
when $d$ is large, we also consider the case where we have noisy access to the
moments of the projection of $\boldsymbol{\vartheta}$ onto a lower dimension
subspace. Such noisy moments can often be obtained by an affine transformation
of the samples in concrete applications such as topic models. Also motivated
by topic models, we measure the accuracy of our estimation
$\widetilde{\boldsymbol{\vartheta}}$ in term of the $L_{1}$-transportation
distance between probability distributions (see the formal definition in
Section 2), which is a more stringent metric than $L_{2}$-transportation. For
example, the distance between two discrete distributions $(1/d,\ldots,1/d)$
and $(2/d,\ldots,2/d,0,\ldots,0)$ is $1$ in $L_{1}$ metric of $\mathbb{R}^{d}$
but only $\sqrt{1/d}$ in $L_{2}$ metric. These two distributions should be
regarded as very different distributions in topic model and thus $L_{1}$
should be a more appropriate distance measure for our setting.

The problems we study are noisy versions of the classical moment problems [57,
44] and have applications to a variety of unsupervised learning scenarios,
including learning topic models [35, 12, 54], learning Gaussian mixtures [61],
collaborative filtering [36, 43], learning mixture of product distributions
[28, 32] and causal inference [33]. For example, one important application of
the above model is learning topic models. We are given a corpus of documents.
We adopt the popular “bag of words” model and take each document as an
unordered multiset of words. The assumption is that there is a small number
$k$ of “pure” topics, where each topic is a distribution over the underlying
vocabulary of $d$ words. A $K$-word document (i.e., a string in $[d]^{K}$) is
generated by first selecting a topic $p\in\Delta_{d-1}$ from the mixture
$\boldsymbol{\vartheta}$, and then sampling $K$ words from this topic. $K$ is
called the snapshot number and we call such a sample a $K$-snapshot of $p$.
For all $\|\mathbf{t}\|_{1}\leq K$, we can obtain noisy estimates of
$M_{\mathbf{t}}(\boldsymbol{\vartheta})$ or
$M_{\mathbf{t}}(\Pi_{B}(\boldsymbol{\vartheta}))$ (the projection of
$\boldsymbol{\vartheta}$ onto a lower dimensional subspace $B$) using the
$K$-snapshots (i.e., the documents in the corpus) . Our goal here is to
recover the mixture $\boldsymbol{\vartheta}$. We will discuss the application
to topic models in Section 5 and also to Gaussian mixture in Section 6. The
mixture model is indeed a special case of the mixture models of $k$ product
distributions [28, 32], which is further a special case of the more general
problem of the mixture models of Bayesian networks [33].

###  1.2 Previous Work and Connections

We first discuss prior work on the $k$-coin model ($\boldsymbol{\vartheta}$ is
supported on $[0,1]$). Rabani et al. [55] showed that recovering a $k$-spike
mixture requires at least the first $K=2k-1$ moments in the worst case (even
without any noise). Moreover, they showed that for the topic learning problem,
if one use $K=c(2k-1)$ snapshot samples (for any constant $c\geq 1)$ and wish
to achieve an accuracy $O(1/k)$ in terms of transportation distance,
$e^{\Omega(K)}$ samples are required (or equivalently, the moment accuracy
should be at most $e^{-\Omega(K)}$). On the positive side, they solved the
problem using sample complexity
$\max\\{(1/\zeta)^{O(k)},(k/\epsilon)^{O(k^{2})}\\}$ (or moment accuracy
$\min\\{\zeta^{O(k)},(\epsilon/k)^{O(k^{2})}\\}$), where $\zeta$ is a lower
bound of the minimum separation ). In fact, their algorithm is a variant of
the classic Prony method and the (post-sampling) running time is
$\operatorname{poly}(k)$ where their bottleneck is solving a special convex
quadratic programming instance. Li et al. [45] provided a linear programming
based algorithm that requires the moment accuracy $(\epsilon/k)^{O(k)}$ (or
sample complexity $(k/\epsilon)^{O(k)}$), matching the lower bound in [55].
However, after computing the noisy moments, the running time of their
algorithm is super-exponential $k^{O(k^{2})}$. Motivated by a problem in
population genetics that reduces to the k-coin mixture problem, Kim et al.
[42] analyzed the Matrix Pencil Method, which requires moment accuracy
$\zeta^{O(k)}\cdot w_{\min}^{2}\cdot\epsilon$. The matrix pencil method
requires solving a generalized eigenvalue problem, which needs $O(k^{3})$
time. Wu and Yang [61] studied the same problem in the context of learning
Gaussian mixtures. In fact, they showed that learning 1-d Gaussian mixture
models with same variance can be reduced to learning the $k$-coin model. Their
algorithm achieves the optimal dependency on the moment accuracy without
separation assumption. Their recovery algorithm is based on SDP and runs in
time $O(k^{6.5})$. Recently, Gordon et al. [31] showed that under the
separation assumption, it is possible to recover all $k$ spikes up to accuracy
$\epsilon$, using the first $2k$-moments with accuracy $\zeta^{O(k)}\cdot
w_{\min}\cdot\epsilon$. Their algorithm is also a variant of Prony’s method
and it runs in time $O(k^{2+o(1)})$. We summarize the prior results in Table
1.

For the high-dimensional case, especially when $d$ is large, it is costly to
obtain all moments $M_{\mathbf{t}}(\boldsymbol{\vartheta})$ for all
$\mathbf{t}$ with $\|\mathbf{t}\|_{1}\leq K$. Indeed, most previous work on
topic modeling solved the problem using noisy moment information of the
projection of $\boldsymbol{\vartheta}$ onto some lower dimensional subspaces.
As different previous work used different projection methods, it is difficult
to compare the accuracy of their moment estimations, as in the 1-dimensional
case. Hence, we directly state the sample complexity result in the context of
topic modeling. Rabani et al. [55] showed that it is possible to produce an
estimate $\widetilde{\boldsymbol{\vartheta}}$ of the original mixture
$\boldsymbol{\vartheta}$ such that
$\operatorname{Tran}(\boldsymbol{\vartheta},\widetilde{\boldsymbol{\vartheta}})\leq\epsilon$,
using $\operatorname{poly}(n,k,1/\epsilon)$ 1- and 2-snapshot samples, and
$(k/\epsilon)^{O(k^{2})}\cdot(w_{\min}\zeta)^{-O(k^{2})}$ $K$-snapshot
samples, under minimal separation assumptions. Li et al. [45] provided an LP-
based learning algorithm which uses almost the same number of samples (with a
slightly worse polynomial for 2-snapshot samples), but without requiring any
separation assumptions. The number of the $(2k-1)$-snapshot samples used in
both [55, 45] is $(k/\epsilon)^{O(k^{2})}$ while that of the lower bound (for
1-d) is only $e^{\Omega(k)}$. Recently, Gordon et al. [31] showed that under
the minimal separation assumption, the sample complexity can be reduced to
$(k/\epsilon)^{O(k)}\cdot w_{\min}^{-O(1)}\cdot\zeta^{-O(k)}$. We summarize
the prior results in Table 2.

Reference | Moment Number $(K)$ | Moment Accuracy $(\xi)$ | Running Time | Separation
---|---|---|---|---
[55] | $2k-1$ | $\min\\{\zeta^{O(k)},(\epsilon/k)^{O(k^{2})}\\}$ | $\textrm{poly}(k)$ | Required
[45] | $2k-1$ | $(\epsilon/k)^{O(k)}$ | $(k/\epsilon)^{O(k^{2})}$ | No Need
[42] | $2k-1$ | $\zeta^{O(k)}\cdot w_{\min}^{2}\cdot\epsilon$ | $O(k^{3})$ | Required
[61] | $2k-1$ | $(\epsilon/k)^{O(k)}$ | $O(k^{6.5})$ | No Need
[31] | $2k$ | $\zeta^{O(k)}\cdot w_{\min}\cdot\epsilon$ | $O(k^{2+o(1)})$ | Required
Theorem 1.1 | $2k-1$ | $(\epsilon/k)^{O(k)}$ | $O(k^{2})$ | No Need
Table 1: Algorithms for the $k$-coin problem. The last column indicates
whether the algorithm needs the separation assumption. In particular, the
algorithms that require the separation assumption need to know $\zeta$ and
$w_{\min}$ where $\zeta$ is the minimum separation, i.e.,
$\zeta\leq\min_{i\neq j}|\alpha_{i}-\alpha_{j}|\leq 1/k$, and
$w_{\min}\leq\min_{i}w_{i}$.

###  1.3 Our Contributions

The $k$-Coin model. We first state our result for the 1-dimensional problem.

######  Theorem 1.1.

(The $k$-Coin model) Let $\boldsymbol{\vartheta}$ be an arbitrary $k$-spike
distribution over $[0,1]$. Suppose we have noisy moments $M^{\prime}_{t}$ such
that $|M_{t}(\boldsymbol{\vartheta})-M^{\prime}_{t}|\leq(\epsilon/k)^{O(k)}$
for $0\leq t\leq 2k-1$. We can obtain a mixture
$\widetilde{\boldsymbol{\vartheta}}$ such that
$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(\epsilon)$ and our algorithm only uses $O(k^{2})$ arithmetic operations.

Our techniques: Our algorithm for the $k$-coin model is also based on the
classic Prony’s algorithm [24]. Suppose the true mixture is
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$ where
$\boldsymbol{\alpha}=[\alpha_{1},\cdots,\alpha_{k}]^{\top}$ and
$\boldsymbol{w}=[w_{1},\cdots,w_{k}]^{\top}$. In Prony’s algorithm, if we know
the moment vector $M$ exactly, all $\alpha_{i}$s can be recovered from the
roots of a polynomial whose coefficients are the entries of the eigenvector
$\boldsymbol{c}$ (corresponding to eigenvalue 0) of a Hankel matrix (see e.g.,
the matrix ${\mathcal{H}}_{k+1}$ in [31]). However, if we only know the noisy
moment vector $M^{\prime}$, we only get a perturbation of the original Hankel
matrix. Recent analyses of Prony’s method [55, 31] aim to bound the error of
recovery from the perturbed Hankel matrix. If $\alpha_{i}$ and $\alpha_{j}$
($i\neq j$) are very close, or any $w_{i}$ is very small, the 2nd smallest
eigenvalue of the Hankel matrix is also very close to 0, hence the the
eigenvector $\boldsymbol{c}$ corresponding to eigenvalue 0 is extremely
sensitive to small perturbation (see e.g.,[59]). Hence, recent analysis of
Prony’s method [55, 31] requires that $\alpha_{i}$ and $\alpha_{j}$ are
separated by at least $\zeta$ and the minimum $w_{i}$ is lower bounded by
$w_{\min}$ to ensure certain stability of the eigenvector (hence the
coefficients of the polynomial), and the corresponding sample complexity
becomes unbounded when $\zeta$ or $w_{\min}$ approaches to 0. Our algorithm
(Algorithm 1) can be seen as a robust version of Prony’s algorithm. Instead of
computing the smallest eigenvector of the perturbed Hankel matrix, we solve a
ridge regression to obtain a vector $\widehat{\boldsymbol{c}}$, which plays a
similar role as $\boldsymbol{c}$ (Line 2 in Algorithm 1). However, we do not
show $\widehat{\boldsymbol{c}}$ is close to $\boldsymbol{c}$ (in fact, they
can be very different). Instead, we adopt a more global analysis to show that
the moment vector of the estimated mixture is close to the true moment vector
$M$, and by a moment-transportation inequality (see Section B), we can
guarantee the quantity of our solution in terms of transportation distance.

Another technical challenge lies in bounding the error in recovering the
weight vector $\boldsymbol{w}$. In the noiseless case, the weight is simply
the solution of $V_{\boldsymbol{\alpha}}\boldsymbol{w}=M$ where
$V_{\boldsymbol{\alpha}}$ is a Vandemonde matrix (see Equation (4)). It is
known that Vandermonde matrices tend to be badly ill-conditioned [30, 53, 49]
(with a large dependency on the inverse of the minimal separation). Hence,
using standard condition number-based analysis, slight perturbations of
$V_{\boldsymbol{\alpha}}$ and $M$ may result in an unbounded error. 111 Note
that for Vandermonde matrix $V=(a^{j}_{i})_{0\leq i<k,1\leq j\leq k}$, its
determinant is $\det(V)=(-1)^{k(k-1)/2}\prod_{p<q}(a_{p}-a_{q})$. Hence,
$\det(V)^{-1}$ depends inversely on the separation.  However, we show
interestingly, in our case, that we can bound the error independent of the
separation (Lemma 3.9) via the connection between Vandermonde linear systems
and Schur polynomial (Lemma 3.8). Our technical are useful in analyzing the
perturbation of Vandermonde linear systems in other contexts as well.

Reference | Snapshot Number $(K)$ | Sample Complexity for $K$-Snapshots | Separation
---|---|---|---
[55] | $2k-1$ | $(k/\epsilon)^{O(k^{2})}\cdot(w_{\min}\zeta)^{-O(k^{2})}$ | Required
[45] | $2k-1$ | $(k/\epsilon)^{O(k^{2})}$ | No Need
[31] | $2k$ | $(k/\epsilon)^{O(k)}\cdot w_{\min}^{-O(1)}\cdot\zeta^{-O(k)}$ | Required
Theorem 5.3 | $2k-1$ | $(k/\epsilon)^{O(k)}$ | No Need
Table 2: Algorithms for recovering the mixture over $\Delta_{d-1}$. All
algorithms requires $\operatorname{poly}(1/\epsilon,d,k)$ many 1- and
2-snapshots. The last column indicates whether the algorithm needs the minimum
separation assumption.

Higher Dimensions. Now state our result for the high dimensional case. Here,
we assume we have noisy access to the moments of the linear projections
$\Pi_{R}(\boldsymbol{\vartheta})$ of $\boldsymbol{\vartheta}$ onto lines or
2-d planes (one may refer to the formal definition of
$\Pi_{R}(\boldsymbol{\vartheta})$ in Section 4.2).

######  Theorem 1.2.

Let $\boldsymbol{\vartheta}$ be an arbitrary $k$-spike mixture supported in
$\Delta_{d-1}$. Suppose we can access noisy $\mathbf{t}$-moments of
$\Pi_{R}(\boldsymbol{\vartheta})$, with precision $(\epsilon/(dk))^{O(k)}$,
where $\Pi_{R}(\boldsymbol{\vartheta})$ is the projected measure obtained by
applying linear transformation $R$ (with $\|R\|_{\infty}\leq 1$) of
$\boldsymbol{\vartheta}$ onto any $h$-dimensional subspace we choose, for all
$\|\mathbf{t}\|\leq K=2k-1$ and $h=1,2$ (i.e., lines and 2-d planes). We can
construct a $k$-spike mixture $\widetilde{\boldsymbol{\vartheta}}$ such that
$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(\epsilon)$ using only $O(dk^{3})$ arithmetic operations with high
probability.

Applications to Learning Topic Models and Gaussian Mixtures. Our result for
the high dimensional case can be easily translated into an algorithm for topic
modeling. We first apply the dimension reduction developed in [45] to find a
special subspace of dimension $O(k)$, then apply our Algorithm 2 to the
projection of $\boldsymbol{\vartheta}$ on this subspace. See Section 5 and
Theorem 5.3. We obtain for the first time the worst case optimal sample
complexity for $K$-snapshots for the high dimensional case, improving previous
work [55, 45, 31] and matching the lower bound even for the 1-dimensional case
[55]. A comparison between our results and prior results can be found in Table
2.

We also study the problem of learning the parameters of Gaussian mixtures. We
assume all Gaussian components share a variance parameter, following the
setting studied in [61]. For the 1-dimensional setting, we can leverage our
algorithm for the $k$-coin model. Our algorithm achieves the same sample
complexity as in [61], but an $O(k^{2})$ post-sampling running time, improving
over the $O(k^{6.5})$ time SDP-based algorithm developed in [61]. See Theorem
6.2.

For the high dimensional setting, we can use the dimension reduction technique
in [45] or [26] to reduce the dimension $d$ to $O(k)$. The dimension reduction
part is not a bottleneck and we assume $d=O(k)$ for the following discussion.
We show in Section 6.2 that we can utilize our algorithm for the high
dimensional sparse moment problem and obtain an algorithm without any
separation assumption with sample complexity $(k/\epsilon)^{O(k)}$. Note that
the algorithm in [61] requires a sample size of $(k/\epsilon)^{O(k)}\cdot
w_{\min}^{-O(1)}\cdot\zeta^{-O(k)}$, which depends on the separation parameter
$\zeta$ between Gaussian distributions. Recently, Doss et al. [26] removed the
separation assumption and achieved the optimal sample complexity
$(k/\epsilon)^{O(k)}$. Comparing with [26], the sample complexity of our
algorithm is the same, but our running time is substantially better: during
the sampling phase, the algorithm in [26] requires
$O(n^{5/4}\operatorname{poly}(k))$ time where $n=(k/\epsilon)^{O(k)}$ is the
number of samples (for each sample, they need to update
$n^{1/4}\operatorname{poly}(k)$ numbers, since their algorithm requires
$n^{1/4}$ 1-d projections) while our algorithm only needs
$O(n\operatorname{poly}(k))$ time (we only need one 1-d projection and
$d=O(k)$ 2-d projections). The post-sampling running time of the algorithm in
[26] is exponential $(k/\epsilon)^{O(k)}$ while our algorithms runs in
polynomial time $\operatorname{poly}(k)$. For more details, see Section 6.

Finally, we argue that although the sample complexity is exponential (for the
above problems), improving the post-sampling running time is also very
important. During the sampling phase, we only need to keep track of the first
few moments (e.g., using basic operations such as counting or adding numbers),
hence the sampling phase can be easily distributed, streamed, or implemented
in inexpensive computing devices. Our mixture recovery algorithm only requires
the moment information (without storing the samples) and runs in time
polynomial in $k$ (not the sample size). Moreover, one may well have other
means to measure the moments to achieve the desired moment accuracy (e.g., via
longer documents, prior knowledge, existing samples etc.) Exploiting other
settings and applications is an interesting further direction.

Our techniques: We first solve the 2-dimensional problem. This is done by
simply extending the previous 1-dimensional algorithm and its analysis to
complex numbers. The real and imaginary parts of a complex location can be
used to represent the two dimensions, respectively. While most ideas are
similar to the 1-dimensional case, the analysis requires extra care in various
places. In particular, it also requires us to extend the moment-transportation
inequality to complex numbers (Appendix B).

Our algorithm for the high dimensional case uses the algorithms for 1- and
2-dimensional cases as subroutines. We first pick a random unit vector
$\boldsymbol{r}$ and learn the projection $\widetilde{\boldsymbol{\phi}}$ of
$\boldsymbol{\vartheta}$ onto $\boldsymbol{r}$ using the 1-dimensional
algorithm. It is not difficult to show that for two spikes that are far apart,
with reasonable probability, their projections on $\boldsymbol{r}$ are also
far apart. Hence, the remaining task is to recover the coordinates
($\in\mathbb{R}^{d}$) of each spike in $\widetilde{\boldsymbol{\phi}}$. Then,
for each dimension $i\in[d]$, we learn the 2d-linear map
$\widetilde{\boldsymbol{\vartheta}}_{i}$ of $\boldsymbol{\vartheta}$ onto the
2-dimensional subspace spanned by $\boldsymbol{r}$ and the $i$th dimensional
axis. By assembling the projected measure $\widetilde{\boldsymbol{\vartheta}}$
and $\widetilde{\boldsymbol{\vartheta}}_{i}$, we show that we can recover the
$i$th coordinates of all spikes.

We remark that the improvement of our algorithm over prior work can be
attributed to the usage of the 2-dimensional projections. Note that all prior
works [55, 45, 61, 31, 26] recover the high dimensional distribution by
somehow “assembling” many 1-dimensional projections. In particular, if one
projects $\boldsymbol{\vartheta}$ to $O(k)$ directions such as in [55, 61,
31], all existing techniques require the separation assumption for the
recovery. 222 It is possible to convert an algorithm with separation
assumption to an algorithm without separation assumption by merging close-by
spikes and removing spikes with small weight. However, the resulting sample
complexity is not optimal. For example, sample complexity
$(k/\epsilon)^{O(k)}\cdot w_{\min}^{-O(1)}\cdot\zeta^{-O(k)}$ in [31] can be
converted to one without separation with sample complexity
$(k/\epsilon)^{O(k^{2})}$. We omit the details.  Other works without
separation condition such as [45, 26] requires to project
$\boldsymbol{\vartheta}$ to a net of exponentially many directions. We hope
this idea of using 2-d projections is useful in solving other high dimensional
recovery problems.

###  1.4 Other Related Work

Learning statistical mixture models has been studied extensively for the past
two decades. A central problem in this area was the problem of learning a
mixture of high-dimensional Gaussians, even robustly [21, 22, 7, 60, 40, 1,
29, 13, 39, 11, 50, 47, 8, 61, 26]. Many other structured mixture models have
also been studied (see e.g., [41, 19, 10, 51, 20, 28, 40, 16, 23, 46]). Our
problem is closely related to topic models which have also been studied
extensively recently [6, 3, 4, 5]. Most work all make certain assumptions on
the structure of the mixture, such as the pure topics are separated, Dirichlet
prior, the existence of anchor words, or certain rank conditions. Some
assumptions (such as [6, 3]) allow one to use documents of constant length
(independent of the number of pure topics $k$ and the desired accuracy
$\epsilon$).

The problem we study can be seen as a sparse version of the classic moment
problem in which our goal is to invert the mapping that takes a measure to the
sequences of moments [57, 44]. When the measure is supported on a finite
interval, the problem is known as the Hausdorff moment problem.

Our problem is also related to the super-resolution problem in which each
measurement takes the form $v_{\ell}=\int e^{2\pi\ell
t}\mathrm{d}\boldsymbol{\vartheta}(t)+\xi_{\ell}$ where $\xi_{\ell}$ is the
noise of the measurement. We can observe the first few $v_{\ell}$ for
$|\ell|\leq K$ where $K$ is called cutoff frequency, and the goal is to
recover the original signal $\boldsymbol{\vartheta}$. There is a long history
of the estimation problem. The noiseless can also be solved by Prony’s method
[24], ESPRIT algorithm [56] or matrix pencil method [37], if $K\geq k$ (i.e.,
we have $2k+1$ measurements). The noisy case has also been studied extensively
and a central goal is to understand the relations between the cutoff
frequency, the size of measure noises, and minimum separation [25, 15, 14, 49,
38, 17]. Various properties, such as the condition number of the Vandermonde
matrix, also play essential roles in this line of study [49]. The relation
between the Vandermonde matrix and Schur polynomial is also exploited in [17].

##  2 Preliminaries

We are given an statistical mixture $\boldsymbol{\vartheta}$ of $k$ discrete
distributions over $[d]=\\{1,2,\cdots,d\\}$. Each discrete distribution
$\boldsymbol{\alpha}_{i}$ can be regarded as a point in the $(d-1)$-simplex
$\Delta_{d-1}=\\{\boldsymbol{x}=(x_{1},\cdots,x_{d})\in\mathbb{R}^{d}\mid\sum_{i=1}^{d}x_{i}=1,x_{i}\geq
0\,\,\forall i\in[d]\\}$. We use
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$ to represent the
mixture where
$\boldsymbol{\alpha}=\\{\boldsymbol{\alpha}_{1},\boldsymbol{\alpha}_{2},\cdots,\boldsymbol{\alpha}_{k}\\}\subset\Delta_{d-1}$
are the locations of spikes and
$\boldsymbol{w}=\\{w_{1},w_{2},\cdots,w_{k}\\}\in\Delta_{k-1}$ in which
$w_{i}$ is the probability of $\boldsymbol{\alpha}_{i}$. Since the dimension
of $\Delta_{d-1}$ is $d-1$, we say the dimension of the problem is $d-1$.
Moreover, we use $\mathrm{Spike}(\Delta_{d-1},\Delta_{k-1})$ to denote the set
of all such mixtures where $\Delta_{d-1}$ indicates the domain of
$\boldsymbol{\alpha}_{i}$ and $\Delta_{k-1}$ indicates the domain of
$\boldsymbol{w}$. Our algorithms may produce negative or complex weights as
intermediate results, so we further denote
$\Sigma^{\mathbb{R}}_{d-1}=\\{\boldsymbol{x}=(x_{1},\cdots,x_{d})\in\mathbb{R}^{d}\mid\sum_{i=1}^{d}x_{i}=1\\}$
and
$\Sigma^{\mathbb{C}}_{d-1}=\\{\boldsymbol{x}=(x_{1},\cdots,x_{d})\in\mathbb{C}^{d}\mid\sum_{i=1}^{d}x_{i}=1\\}$.

For the 1-dimensional case (which is called the $k$-coin problem or sparse
Hausdorff moment problem), we have
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})\in\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$,
where $\boldsymbol{\alpha}=\\{\alpha_{1},\alpha_{2},\cdots,\alpha_{k}\\}$ is a
set of $k$ discrete points in $[0,1]$. In this scenario, for each
$t\in\mathbb{N}$, we denote the $t$th moment as

$M_{t}(\boldsymbol{\vartheta})=\int_{[0,1]}\alpha^{t}\boldsymbol{\vartheta}(\mathrm{d}\alpha)=\sum_{i=1}^{k}w_{i}\alpha_{i}^{t}.$

For the higher dimensional case, we use
$t=(t_{1},t_{2},\cdots,t_{d})\in\mathbb{Z}_{+}^{d}$ such that
$\|t\|_{1}=K(=2k-1)$ to denote a multi-index. In addition, we denote
$\boldsymbol{\alpha}_{i}^{\mathbf{t}}=\alpha_{i,1}^{t_{1}}\alpha_{i,2}^{t_{2}}\cdots\alpha_{i,d}^{t_{d}}$
for every discrete distribution
$\boldsymbol{\alpha}_{i}=(\alpha_{i,1},\alpha_{i,2},\cdots,\alpha_{i,d})\in\mathbb{R}^{d}$.
The moment vector of $\boldsymbol{\vartheta}$ is then defined as

$\displaystyle
M(\boldsymbol{\vartheta})=\left\\{M_{\mathbf{t}}(\boldsymbol{\vartheta})\right\\}_{\|\mathbf{t}\|_{1}\leq
K}\text{ where
}M_{\mathbf{t}}(\boldsymbol{\vartheta})=\int_{\Delta_{d-1}}\boldsymbol{\alpha}^{\mathbf{t}}\boldsymbol{\vartheta}(\mathrm{d}\boldsymbol{\alpha})=\sum_{i=1}^{k}w_{i}\boldsymbol{\alpha}_{i}^{\mathbf{t}}.$
(1)

We define the moment distance between two mixtures $\boldsymbol{\vartheta}$
and $\boldsymbol{\vartheta}^{\prime}$ as the $L_{\infty}$ norm of the
difference between corresponding moment vectors,

$\operatorname{Mom}_{K}(\boldsymbol{\vartheta},\boldsymbol{\vartheta}^{\prime}):=\max_{|\mathbf{t}|\leq
K}|M_{\mathbf{t}}(\boldsymbol{\vartheta})-M_{\mathbf{t}}(\boldsymbol{\vartheta}^{\prime})|.$

Transportation Distance: Recall that for any two probability measures $P$ and
$Q$ defined over $\mathbb{R}^{d}$, the $L_{1}$-transportation distance
$\operatorname{Tran}(P,Q)$ is defined as

$\displaystyle\operatorname{Tran}(P,Q):=\inf\left\\{\int\|x-y\|_{1}\,\mathrm{d}\mu(x,y):\mu\in
M(P,Q)\right\\},$ (2)

where $M(P,Q)$ is the set of all joint distributions (also called coupling) on
$\mathbb{R}^{d}\times\mathbb{R}^{d}$ with marginals $P$ and $Q$.
Transportation distance is also called Rubinstein distance, Wasserstein
distance or earth mover distance in the literature.

Let $1\text{-}\mathsf{Lip}$ be the set of 1-Lipschitz functions on
$\mathbb{R}^{d}$, i.e.,
$1\text{-}\mathsf{Lip}:=\\{f:\mathbb{R}^{d}\rightarrow\mathbb{R}\mid|f(x)-f(y)|\leq\|x-y\|_{1}\text{
for any }x,y\in\mathbb{R}^{d}\\}$. We need the following important theorem by
Kantorovich and Rubinstein (see e.g., [27]), which states the dual formulation
of transportation distance:

$\displaystyle\operatorname{Tran}(P,Q)=\sup\left\\{\int f\mathrm{d}(P-Q):f\in
1\text{-}\mathsf{Lip}\right\\}.$ (3)

##  3 An Efficient Algorithm for the $k$-Coin Problem

In this section, we study the $k$-coin model, i.e., $\boldsymbol{\vartheta}$
is a $k$-spike distribution over $[0,1]$. We present an efficient algorithm
reconstructing the mixture from the first $K=2k-1$ noisy moments. We denote
$\boldsymbol{\vartheta}:=(\boldsymbol{\alpha},\boldsymbol{w})\in\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$
as the ground truth mixture where
$\boldsymbol{\alpha}=[\alpha_{1},\cdots,\alpha_{k}]^{\top}\in\Delta_{1}^{k}$
and $\boldsymbol{w}=[w_{1},\cdots,w_{k}]^{\top}\in\Delta_{k}$. Let
$M(\boldsymbol{\vartheta})=[M_{0}(\boldsymbol{\vartheta}),\cdots,M_{2k-1}(\boldsymbol{\vartheta})]^{\top}$
be the ground truth moment vector containing moments of degree at most
$K=2k-1$. We denote the noisy moment vector as
$M^{\prime}=[M_{0}^{\prime},\cdots,M_{2k-1}^{\prime}]$ such that the error is
bounded by $\|M^{\prime}-M(\boldsymbol{\vartheta})\|_{\infty}\leq\xi$. We
further assume $M^{\prime}_{0}=M_{0}(\boldsymbol{\vartheta})=1$ (because it is
a distribution) and the noise satisfies $\xi\leq 2^{-\Omega(k)}$ (this is
necessary in light of the lower bound in [55] or Lemma B.1).

For $M=[M_{0},\cdots,M_{2k-1}]^{\top}$ and
$\boldsymbol{\alpha}=[\alpha_{1},\cdots,\alpha_{k}]^{\top}$, we denote

$\displaystyle A_{M}:=\begin{bmatrix}M_{0}&M_{1}&\cdots&M_{k-1}\\\
M_{1}&M_{2}&\cdots&M_{k}\\\ \vdots&\vdots&\ddots&\vdots\\\
M_{k-1}&M_{k}&\cdots&M_{2k-2}\end{bmatrix},\quad
b_{M}:=\begin{bmatrix}M_{k}\\\ M_{k+1}\\\ \vdots\\\
M_{2k-1}\end{bmatrix},\quad
V_{\boldsymbol{\alpha}}:=\begin{bmatrix}\alpha_{1}^{0}&\alpha_{2}^{0}&\cdots&\alpha_{k}^{0}\\\
\alpha_{1}^{1}&\alpha_{2}^{1}&\cdots&\alpha_{k}^{1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\alpha_{1}^{2k-1}&\alpha_{2}^{2k-1}&\cdots&\alpha_{k}^{2k-1}\end{bmatrix}.$
(4)

Note that $A_{M}$ is a Hankel matrix and $V_{\boldsymbol{\alpha}}$ is a
Vandemonde matrix.

Algorithm 1 Reconstruction Algorithm in 1-Dimension

1:function OneDimension($k,M^{\prime},\xi$)

2:
$\widehat{\boldsymbol{c}}\leftarrow\arg\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|A_{M^{\prime}}\boldsymbol{x}+b_{M^{\prime}}\|_{2}^{2}+\xi^{2}\|\boldsymbol{x}\|_{2}^{2}$
$\triangleright$
$\widehat{\boldsymbol{c}}=[\widehat{c}_{0},\cdots,\widehat{c}_{k-1}]^{\top}\in\mathbb{R}^{k}$

3:
$\widehat{\boldsymbol{\alpha}}\leftarrow\textrm{roots}(\sum_{i=0}^{k-1}\widehat{c}_{i}x^{i}+x^{k})$
$\triangleright$
$\widehat{\boldsymbol{\alpha}}=[\widehat{\alpha}_{1},\cdots,\widehat{\alpha}_{k}]^{\top}\in\mathbb{C}^{k}$

4:
$\overline{\boldsymbol{\alpha}}\leftarrow\textrm{project}_{\Delta_{1}}(\widehat{\boldsymbol{\alpha}})$
$\triangleright$
$\overline{\boldsymbol{\alpha}}=[\overline{\alpha}_{1},\cdots,\overline{\alpha}_{k}]^{\top}\in\Delta_{1}^{k}$

5:
$\widetilde{\boldsymbol{\alpha}}\leftarrow\overline{\boldsymbol{\alpha}}+\textrm{Noise}(\xi)$
$\triangleright$
$\widetilde{\boldsymbol{\alpha}}=[\widetilde{\alpha}_{1},\cdots,\widetilde{\alpha}_{k}]^{\top}\in\Delta_{1}^{k}$

6:
$\widehat{\boldsymbol{w}}\xleftarrow{O(1)-\textrm{approx}}\arg\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M^{\prime}\|_{2}^{2}$
$\triangleright$
$\widehat{\boldsymbol{w}}=[\widehat{w}_{1},\cdots,\widehat{w}_{k}]^{\top}\in\mathbb{R}^{k}$

7:
$\widetilde{\boldsymbol{w}}\leftarrow\widehat{\boldsymbol{w}}/(\sum_{i=1}^{k}\widehat{w}_{i})$
$\triangleright$
$\widetilde{\boldsymbol{w}}=[\widetilde{w}_{1},\cdots,\widetilde{w}_{k}]^{\top}\in\Sigma_{k-1}$

8:
$\widetilde{\boldsymbol{\vartheta}}\leftarrow(\widetilde{\boldsymbol{\alpha}},\widetilde{\boldsymbol{w}})$
$\triangleright$
$\widetilde{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\Delta_{1},\Sigma_{k-1})$

9:
$\widecheck{\boldsymbol{w}}\leftarrow\arg\min_{\boldsymbol{x}\in\Delta_{k-1}}\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},(\widetilde{\boldsymbol{\alpha}},\boldsymbol{x}))$

10: report
$\widecheck{\boldsymbol{\vartheta}}\leftarrow(\widetilde{\boldsymbol{\alpha}},\widecheck{\boldsymbol{w}})$
$\triangleright$
$\widecheck{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$

11:end function

Our algorithm is a variant of Prony’s algorithm. The pseudocode can be found
in Algorithm 1. The algorithm takes the number of spikes $k$, the noisy moment
vector $M^{\prime}$ which
$\|M^{\prime}-M(\boldsymbol{\vartheta})\|_{\infty}\leq\xi$ and the moment
accuracy $\xi$ as the input. We describe our algorithm as follows. Let
$\boldsymbol{c}=[c_{0},\cdots,c_{k-1}]^{\top}\in\mathbb{R}^{k}$ such that
$\prod_{i=1}^{k}(x-\alpha_{i})=\sum_{i=0}^{k-1}c_{i}x^{i}+x^{k}$ be the
characteristic vector of locations $\alpha_{i}$. We first perform a ridge
regression to obtain $\widehat{\boldsymbol{c}}$ in Line 2. Note that
$A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}=\mathbf{0}.$
(see Lemma 3.1). Hence, $\widehat{\boldsymbol{c}}$ serves a similar role as
$\boldsymbol{c}$ (note that $\widehat{\boldsymbol{c}}$ is not necessarily
close to $\boldsymbol{c}$ without the separation assumption). From Line 3 to
Line 5, we aim to obtain estimations of the positions of the spikes, i.e.,
$\alpha_{i}$s. We first solve the roots of polynomial
$\sum_{i=0}^{k-1}\widehat{c}_{i}x^{i}+x^{k}$. For polynomial root findings,
note that some roots may be complex without any separation assumption.
333Gordon et al. [31] prove that all roots are real and separated under the
minimal separation assumption.  Hence, we need to project the solutions back
to $\Delta_{1}$ and inject small noise, ensuring that all values are distinct
and still in $\Delta_{1}$. We note that any noise of size at most $\xi$
suffices. After recovering the positions of the spikes, we aim to recover the
corresponding weights $\widetilde{\boldsymbol{w}}$ by the linear regression
defined by the Vandemonde matrix $V_{\boldsymbol{\alpha}}$ (Line 6). We
normalized the weight vector in Line 7. We note that
$\widetilde{\boldsymbol{w}}$ may still have some negative components. In Line
9, we find the closest $k$-spike distribution in
$\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$, which is our final output. The
details for implementing the above steps in $O(k^{2})$ time can be found in
Appendix C.1.

###  3.1 Error Analysis

Now, we start to bound the reconstruction error of the algorithm. The
following lemma is well known (see e.g., [18, 31]). We provide a proof in the
appendix for completeness.

######  Lemma 3.1.

Let
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})\in\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$
where $\boldsymbol{\alpha}=[\alpha_{1},\cdots,\alpha_{k}]^{\top}$ and
$\boldsymbol{w}=[w_{1},\cdots,w_{k}]^{\top}$. Let
$\boldsymbol{c}=[c_{0},\cdots,c_{k-1}]^{\top}\in\mathbb{R}^{k}$ such that
$\prod_{i=1}^{k}(x-\alpha_{i})=\sum_{i=0}^{k-1}c_{i}x^{i}+x^{k}$. For all
$i\geq 0$, the following equation holds
$\sum_{j=0}^{k-1}M_{i+j}(\boldsymbol{\vartheta})c_{j}+M_{i+k}(\boldsymbol{\vartheta})=0.$
In matrix form, the equation can be written as follows:
$A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}=\mathbf{0}.$

Next, the following lemma shows that the intermediate result
$\widehat{\boldsymbol{c}}$ a good estimation for the solution of
$A_{M(\boldsymbol{\vartheta})}\boldsymbol{x}+b_{M(\boldsymbol{\vartheta})}=\boldsymbol{0}$
with a small norm:

######  Lemma 3.2.

Let $\boldsymbol{c}=[c_{0},\cdots,c_{k-1}]^{\top}\in\mathbb{R}^{k}$ such that
$\prod_{i=1}^{k}(x-\alpha_{i})=\sum_{i=0}^{k-1}c_{i}x^{i}+x^{k}$. Suppose
$\|M^{\prime}-M(\boldsymbol{\vartheta})\|_{\infty}\leq\xi$. Let
$\widehat{\boldsymbol{c}}=[\widehat{c}_{0},\cdots,\widehat{c}_{k-1}]^{\top}\in\mathbb{R}^{k}$
be the intermediate result (Line 2) in Algorithm 1. Then,
$\|\boldsymbol{c}\|_{1}\leq 2^{k}$, $\|\widehat{\boldsymbol{c}}\|_{1}\leq
2^{O(k)}$ and
$\|A_{M(\boldsymbol{\vartheta})}\widehat{\boldsymbol{c}}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}\leq
2^{O(k)}\cdot\xi$.

###### Proof.

From Vieta’s formulas (see e.g., [9]), we have
$c_{i}=\sum_{S\in\binom{[k]}{k-i}}\prod_{j\in S}(-\alpha_{j}).$ Thus,

$\|\boldsymbol{c}\|_{1}=\sum_{i=0}^{k-1}|c_{i}|\leq\sum_{S\in
2^{[k]}}\prod_{j\in S}|\alpha_{j}|=\prod_{i=1}^{k}(1+|\alpha_{i}|)\leq 2^{k}.$

where the last inequality holds because $|\alpha_{i}|\leq 1$ for all $i$.

From Lemma 3.1, we can see that
$\|A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}=0$.
Therefore,

$\displaystyle\|A_{M^{\prime}}\boldsymbol{c}+b_{M^{\prime}}\|_{\infty}$
$\displaystyle\leq\|A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}+\|A_{M^{\prime}}\boldsymbol{c}-A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}\|_{\infty}+\|b_{M^{\prime}}-b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq\|A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}+\|(A_{M^{\prime}}-A_{M(\boldsymbol{\vartheta})})\boldsymbol{c}\|_{\infty}+\|b_{M^{\prime}}-b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq\|A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}+\|A_{M^{\prime}}-A_{M(\boldsymbol{\vartheta})}\|_{\infty}\|\boldsymbol{c}\|_{1}+\|b_{M^{\prime}}-b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq\|A_{M(\boldsymbol{\vartheta})}\boldsymbol{c}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}+\|M^{\prime}-M(\boldsymbol{\vartheta})\|_{\infty}(\|\boldsymbol{c}\|_{1}+1)$
$\displaystyle\leq 2^{O(k)}\cdot\xi.$

The fourth inequality holds since
$(A_{M^{\prime}}-A_{M(\boldsymbol{\vartheta})})_{i,j}=(M^{\prime}-M(\boldsymbol{\vartheta}))_{i+j}$
and
$(b_{M^{\prime}}-b_{M(\boldsymbol{\vartheta})})_{i}=(M^{\prime}-M(\boldsymbol{\vartheta}))_{i+k}$.
From the definition of $\widehat{\boldsymbol{c}}$, we can see that

$\displaystyle\|A_{M^{\prime}}\widehat{\boldsymbol{c}}+b_{M^{\prime}}\|_{2}^{2}+\xi^{2}\|\widehat{\boldsymbol{c}}\|_{2}^{2}$
$\displaystyle\leq\|A_{M^{\prime}}\boldsymbol{c}+b_{M^{\prime}}\|_{2}^{2}+\xi^{2}\|\boldsymbol{c}\|_{2}^{2}$
$\displaystyle\leq
k\|A_{M^{\prime}}\boldsymbol{c}+b_{M^{\prime}}\|_{\infty}^{2}+\xi^{2}\|\boldsymbol{c}\|_{1}^{2}$
$\displaystyle\leq 2^{O(k)}\cdot\xi^{2}.$

The second inequality holds since
$\|\boldsymbol{x}\|_{2}\leq\|\boldsymbol{x}\|_{1}$ and
$\|\boldsymbol{x}\|_{2}\leq\sqrt{k}\|\boldsymbol{x}\|_{\infty}$ holds for any
vector $\boldsymbol{x}\in\mathbb{R}^{k}$.

Now, we can directly see that

$\displaystyle\|A_{M^{\prime}}\widehat{\boldsymbol{c}}+b_{M^{\prime}}\|_{\infty}$
$\displaystyle\leq\|A_{M^{\prime}}\widehat{\boldsymbol{c}}+b_{M^{\prime}}\|_{2}\leq
2^{O(k)}\cdot\xi,$ $\displaystyle\|\widehat{\boldsymbol{c}}\|_{1}$
$\displaystyle\leq\sqrt{k}\|\widehat{\boldsymbol{c}}\|_{2}\leq 2^{O(k)}.$

Finally, we can bound
$\|A_{M(\boldsymbol{\vartheta})}\widehat{\boldsymbol{c}}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
as follows:

$\displaystyle\|A_{M(\boldsymbol{\vartheta})}\widehat{\boldsymbol{c}}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq\|A_{M^{\prime}}\widehat{\boldsymbol{c}}+b_{M^{\prime}}\|_{\infty}+\|A_{M^{\prime}}\widehat{\boldsymbol{c}}-A_{M(\boldsymbol{\vartheta})}\widehat{\boldsymbol{c}}\|_{\infty}+\|b_{M^{\prime}}-b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq\|A_{M^{\prime}}\widehat{\boldsymbol{c}}+b_{M^{\prime}}\|_{\infty}+\|A_{M^{\prime}}-A_{M(\boldsymbol{\vartheta})}\|_{\infty}\|\widehat{\boldsymbol{c}}\|_{1}+\|b_{M^{\prime}}-b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq\|A_{M^{\prime}}\widehat{\boldsymbol{c}}+b_{M^{\prime}}\|_{\infty}+\|M^{\prime}-M(\boldsymbol{\vartheta})\|_{\infty}(\|\widehat{\boldsymbol{c}}\|_{1}+1)$
$\displaystyle\leq 2^{O(k)}\cdot\xi.$

∎

Using this result, we can show that $\widetilde{\boldsymbol{\alpha}}$ in Line
5 is a good estimation for the ground truth spikes, $\boldsymbol{\alpha}$. In
particular, the following lemma shows that every ground truth spike
$a_{i}\in\boldsymbol{\alpha}$ has a close spike
$\widetilde{\alpha}_{j}\in\widetilde{\boldsymbol{\alpha}}$ where the closeness
depends on the weight of the spike $w_{i}$. We note that this property
tolerates permutation and small spikes weights, enabling us to analyze
mixtures without separations.

######  Lemma 3.3.

Let
$\widehat{\boldsymbol{\alpha}}=[\widehat{\alpha}_{1},\cdots,\widehat{\alpha}_{k}]^{\top}\in\mathbb{C}^{k}$
be the intermediate result (Line 3) in Algorithm 1. Then, the following
inequality holds:

$\displaystyle\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widehat{\alpha}_{j}|^{2}\leq
2^{O(k)}\cdot\xi.$

###### Proof.

Since $\widehat{\boldsymbol{c}}$ and $\boldsymbol{\alpha}$ are real, from the
definition of $\widehat{\boldsymbol{c}}$,

$\displaystyle\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widehat{\alpha}_{j}|^{2}=\sum_{i=1}^{k}w_{i}\left|\prod_{j=1}^{k}(\alpha_{i}-\widehat{\alpha}_{j})\right|^{2}=\sum_{i=1}^{k}w_{i}\left|\sum_{j=0}^{k-1}\widehat{c}_{j}\alpha_{i}^{j}+\alpha_{i}^{k}\right|^{2}=\sum_{i=1}^{k}w_{i}\left(\sum_{j=0}^{k-1}\widehat{c}_{j}\alpha_{i}^{j}+\alpha_{i}^{k}\right)^{2}.$

Hence, we can see that

$\displaystyle\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widehat{\alpha}_{j}|^{2}$
$\displaystyle=\sum_{i=1}^{k}w_{i}\left(\sum_{p=0}^{k-1}\sum_{q=0}^{k-1}\widehat{c}_{p}\widehat{c}_{q}\alpha_{i}^{p+q}+2\sum_{p=0}^{k-1}\widehat{c}_{p}\alpha_{i}^{p+k}+\alpha_{i}^{2k}\right)$
$\displaystyle=\sum_{p=0}^{k-1}\sum_{q=0}^{k-1}\widehat{c}_{p}\widehat{c}_{q}\sum_{i=1}^{k}w_{i}\alpha_{i}^{p+q}+2\sum_{p=0}^{k-1}\widehat{c}_{p}\sum_{i=1}^{k}w_{i}\alpha_{i}^{p+k}+\sum_{i=1}^{k}w_{i}\alpha_{i}^{2k}$
$\displaystyle=\sum_{p=0}^{k-1}\widehat{c}_{p}\left(\sum_{q=0}^{k-1}M_{p+q}(\boldsymbol{\vartheta})\widehat{c}_{q}+M_{p+k}(\boldsymbol{\vartheta})\right)+\sum_{p=0}^{k-1}\widehat{c}_{p}M_{p+k}(\boldsymbol{\vartheta})+M_{2k}(\boldsymbol{\vartheta})$
$\displaystyle=\sum_{p=0}^{k-1}\widehat{c}_{p}\Bigl{(}A_{M(\boldsymbol{\vartheta})}\widehat{c}+b_{M(\boldsymbol{\vartheta})}\Bigr{)}_{p}+\sum_{p=0}^{k-1}\widehat{c}_{p}M_{p+k}(\boldsymbol{\vartheta})+M_{2k}(\boldsymbol{\vartheta}).$

According to Lemma 3.1,
$M_{p+k}(\boldsymbol{\vartheta})=-\sum_{q=0}^{k-1}M_{p+q}(\boldsymbol{\vartheta})c_{q}$,
so

$\displaystyle\sum_{p=0}^{k-1}\widehat{c}_{p}M_{p+k}(\boldsymbol{\vartheta})+M_{2k}(\boldsymbol{\vartheta})$
$\displaystyle=\sum_{p=0}^{k-1}\widehat{c}_{p}\left(-\sum_{q=0}^{k-1}M_{p+q}(\boldsymbol{\vartheta})c_{q}\right)+\left(-\sum_{q=0}^{k-1}M_{k+q}(\boldsymbol{\vartheta})c_{q}\right)$
$\displaystyle=-\sum_{q=0}^{k-1}c_{q}\left(\sum_{p=0}^{k-1}M_{p+q}(\boldsymbol{\vartheta})\widehat{c}_{p}+M_{k+q}(\boldsymbol{\vartheta})\right)$
$\displaystyle=-\sum_{q=0}^{k-1}c_{q}\Bigl{(}A_{M(\boldsymbol{\vartheta})}\widehat{c}+b_{M(\boldsymbol{\vartheta})}\Bigr{)}_{q}.$

Combining the above equations, we finally get

$\displaystyle\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widehat{\alpha}_{j}|^{2}$
$\displaystyle=\sum_{i=0}^{k-1}(\widehat{c}_{i}-c_{i})\Bigl{(}A_{M(\boldsymbol{\vartheta})}\widehat{c}+b_{M(\boldsymbol{\vartheta})}\Bigr{)}_{i}$
$\displaystyle\leq(\|\widehat{c}\|_{1}+\|c\|_{1})\|A_{M(\boldsymbol{\vartheta})}\widehat{c}+b_{M(\boldsymbol{\vartheta})}\|_{\infty}$
$\displaystyle\leq(2^{k}+2^{O(k)})\cdot 2^{O(k)}\cdot\xi$
$\displaystyle(\textrm{Lemma}~{}\ref{lm:dim1-step1})$ $\displaystyle\leq
2^{O(k)}\cdot\xi.$

∎

We then give the following simple property that allows us to upper bound the
impact of injected noise. The proof can be found in the appendix.

######  Lemma 3.4.

Let $b>0$ be some constant. For $a_{1},\cdots,a_{k}\in\mathbb{R}$ and
$a_{1}^{\prime},\cdots,a_{k}^{\prime}\in\mathbb{R}$ such that $a_{i}\in[0,b]$
and $|a_{i}^{\prime}-a_{i}|\leq\xi$ holds for all $i$. In case that $\xi\leq
n^{-1}$, then,

$\displaystyle\prod_{i=1}^{k}a_{i}^{\prime}\leq\prod_{i=1}^{k}a_{i}+O(b^{k}\cdot
k\xi).$

With the lemma, we are able to show that the projection (Line 4) and the
injected noise (Line 5) do not introduce much extra error for our estimation
on the positions of ground truth spikes.

######  Lemma 3.5.

Let
$\widetilde{\boldsymbol{\alpha}}=[\widetilde{\alpha}_{1},\cdots,\widetilde{\alpha}_{k}]^{\top}\in\mathbb{R}^{k}$
be the intermediate result (Line 5) in Algorithm 1. Then,

$\displaystyle\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widetilde{\alpha}_{j}|^{2}\leq
2^{O(k)}\cdot\xi.$

###### Proof.

Let
$\overline{\boldsymbol{\alpha}}=[\overline{\alpha}_{1},\cdots,\overline{\alpha}_{k}]^{\top}\in\mathbb{R}^{k}$
be the set of projections (Line 4). Since the projected domain $\Delta_{1}$ is
convex,
$|\alpha_{i}-\overline{\alpha}_{j}|\leq|\alpha_{i}-\widehat{\alpha}_{j}|$
holds. Thus,

$\displaystyle\prod_{j=1}^{k}|\alpha_{i}-\overline{\alpha}_{j}|^{2}\leq\prod_{j=1}^{k}|\alpha_{i}-\widehat{\alpha}_{j}|^{2}.$

Recall that $\widetilde{\alpha}_{j}$ is obtained from $\overline{\alpha}_{j}$
by adding an independent noise of size at most $\xi$. We have
$|\alpha_{i}-\widetilde{\alpha}_{j}|\leq|\alpha_{i}-\overline{\alpha}_{j}|+\xi$.
Applying Lemma 3.4 by regarding $|\alpha_{i}-\widetilde{\alpha}_{j}|$ as
$a_{j}$ and $|\alpha_{i}-\widehat{\alpha}_{j}|$ as $a_{j}^{\prime}$, from
$|\alpha_{i}-\widetilde{\alpha}_{j}|\leq|\alpha_{i}|+|\widetilde{\alpha}_{j}|\leq
2$, we can conclude that

$\displaystyle\prod_{j=1}^{k}|\alpha_{i}-\widetilde{\alpha}_{j}|^{2}\leq\prod_{j=1}^{k}|\alpha_{i}-\overline{\alpha}_{j}|^{2}+O(2^{k}\cdot
k\xi).$

Combining two inequalities, we get the desired inequality:

$\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widetilde{\alpha}_{j}|^{2}\leq\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widehat{\alpha}_{j}|+O(2^{k}\cdot
k\xi)\leq 2^{O(k)}\cdot\xi.$

∎

Now, we start to bound the error in Line 6 which is a linear regression
defined by the Vandermonde matrix. We first introduce Schur polynomial which
has a strong connection to Vandermonde matrix.

######  Definition 3.6 (Schur polynomial).

Given a partition $\lambda=(\lambda_{1},\cdots,\lambda_{k})$ where
$\lambda_{1}\geq\cdots\geq\lambda_{k}\geq 0$. The Schur polynomial
$s_{\lambda}(x_{1},x_{2},\cdots,x_{k})$ is defined by the ratio

$s_{\lambda}(x_{1},x_{2},\cdots,x_{k})=\left|\begin{matrix}x_{1}^{\lambda_{1}+k-1}&x_{2}^{\lambda_{1}+k-1}&\cdots&x_{k}^{\lambda_{1}+k-1}\\\
x_{1}^{\lambda_{2}+k-2}&x_{2}^{\lambda_{2}+k-2}&\cdots&x_{k}^{\lambda_{2}+k-2}\\\
\vdots&\vdots&\ddots&\vdots\\\
x_{1}^{\lambda_{k}}&x_{2}^{\lambda_{k}}&\cdots&x_{k}^{\lambda_{k}}\end{matrix}\right|\cdot\left|\begin{matrix}1&1&\cdots&1\\\
x_{1}^{1}&x_{2}^{1}&\cdots&x_{k}^{1}\\\ \vdots&\vdots&\ddots&\vdots\\\
x_{1}^{k-1}&x_{2}^{k-1}&\cdots&x_{k}^{k-1}\end{matrix}\right|^{-1}.$

A Schur polynomial is a symmetric function because the numerator and
denominator are both alternating (i.e., it changes sign if we swap two
variables), and a polynomial since all alternating polynomials are divisible
by the Vandermonde determinant $\prod_{i<j}(x_{i}-x_{j})$. 444It is easy to
see that if a polynomial $p(x_{1},\ldots,x_{n})$ is alternating, it must
contain $(x_{i}-x_{j})$ as a factor . The following classical result (see,
e.g., [58]) shows an alternative way to calculate Schur polynomial. It plays a
crucial role in bounding the error of a linear system defined by the
Vandermonde matrix (Lemma 3.9).

######  Theorem 3.7.

(e.g. [58]) Let $\textrm{SSYT}(\lambda)$ be the set consists all semistandard
Young tableau of shape $\lambda=(\lambda_{1},\cdots,\lambda_{k})$. Then, Schur
polynomial $s_{\lambda}(x_{1},x_{2},\cdots,x_{k})$ can also be computed as
follows:

$\displaystyle
s_{\lambda}(x_{1},x_{2},\cdots,x_{k})=\sum_{T\in\textrm{SSYT}(\lambda)}\prod_{i=1}^{k}x_{i}^{T_{i}}.$

where $T_{i}$ counts the occurrences of the number $i$ in $T$.

In a semistandard Young tableau, we allow the same number to appear more than
once (or not at all), and we require the entries weakly increase along each
row and strictly increase down each column As an example, for
$\lambda=(2,1,0)$, the list of semistandard Young tableau are:

1 & 1
2 1 & 1
3 1 & 2
2 1 & 2
3 1 & 3
2 1 & 3
3 2 & 2
3 2 & 3
3

The corresponding Schur polynomial is

$\displaystyle
s_{(2,1,0)}(x_{1},x_{2},x_{3})=x_{1}^{2}x_{2}+x_{1}^{2}x_{3}+x_{1}x_{2}^{2}+x_{1}x_{2}x_{3}+x_{1}x_{3}x_{2}+x_{1}x_{3}^{2}+x_{2}^{2}x_{3}+x_{2}x_{3}^{2}.$

Considering the special case where $\lambda=(j-k+1,1,1,\cdots,1)$, this
theorem implies the following result. For completeness, we also present a
self-contained proof of the lemma in the appendix.

######  Lemma 3.8.

For all $j\geq k$ and distinct numbers $a_{1},\cdots,a_{k}\in\mathbb{R}$ (or
$a_{1},\cdots,a_{k}\in\mathbb{C}$),

$\left|\begin{matrix}a_{1}^{j}&a_{2}^{j}&\cdots&a_{k}^{j}\\\
a_{1}^{1}&a_{2}^{1}&\cdots&a_{k}^{1}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{1}^{k-1}&a_{2}^{k-1}&\cdots&a_{k}^{k-1}\end{matrix}\right|\cdot\left|\begin{matrix}1&1&\cdots&1\\\
a_{1}^{1}&a_{2}^{1}&\cdots&a_{k}^{1}\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{1}^{k-1}&a_{2}^{k-1}&\cdots&a_{k}^{k-1}\end{matrix}\right|^{-1}=\prod_{i=1}^{k}a_{i}\cdot\sum_{s\in(k)^{j-k}}\prod_{i=1}^{k}a_{i}^{s_{i}}$

where $(a)^{b}\subset\\{0,\cdots,b\\}^{a}$ is the set that contains all vector
$s$ in which $\sum_{i=1}^{a}s_{i}=b$.

Using this result, we can show the reconstructed moment vector
$V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}$ is close to the
ground true moment vector $M(\boldsymbol{\vartheta})$. The next lemma shows
the discrete mixture with support $\widetilde{\boldsymbol{\alpha}}$ can well
approximate a single spike distribution at $\alpha$ in terms of the moments,
if any of the spikes in $\widetilde{\boldsymbol{\alpha}}$ is close to
$\alpha$.

######  Lemma 3.9.

Let $M(\alpha)$ be the vector
$[\alpha^{0},\alpha^{1},\cdots,\alpha^{2k-1}]^{\top}$. Suppose $F$ is either
$\\{\mathbb{R}$ or $\mathbb{C}\\}$. For $\alpha\in F$ with $|\alpha|\leq 1$
and
$\widetilde{\boldsymbol{\alpha}}=[\widetilde{\alpha}_{1},\cdots,\widetilde{\alpha}_{k}]^{\top}\in
F^{k}$ with all $\widetilde{\alpha}$s distinct and
$\|\widetilde{\boldsymbol{\alpha}}\|_{\infty}\leq 1$. We have,

$\min_{\boldsymbol{x}\in
F^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M(\alpha)\|_{2}\leq
2^{O(k)}\prod_{j=1}^{k}|\alpha-\widetilde{\alpha}_{j}|.$

###### Proof.

Let $\Delta\widetilde{\alpha}_{j}=\widetilde{\alpha}_{j}-\alpha$. Consider
$\boldsymbol{x}^{*}=[x_{1}^{*},\cdots,x_{k}^{*}]^{\top}$ such that

$\displaystyle\begin{bmatrix}1&1&\cdots&1\\\
\Delta\widetilde{\alpha}_{1}^{1}&\Delta\widetilde{\alpha}_{2}^{1}&\cdots&\Delta\widetilde{\alpha}_{k}^{1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\Delta\widetilde{\alpha}_{1}^{k-1}&\Delta\widetilde{\alpha}_{2}^{k-1}&\cdots&\Delta\widetilde{\alpha}_{k}^{k-1}\end{bmatrix}\begin{bmatrix}x_{1}^{*}\\\
x_{2}^{*}\\\ \vdots\\\ x_{k}^{*}\\\ \end{bmatrix}=\begin{bmatrix}1\\\ 0\\\
\vdots\\\ 0\\\ \end{bmatrix}.$

We see that for each $j$,

$\displaystyle(V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}^{*}-M(\alpha))_{j}$
$\displaystyle=\sum_{t=1}^{k}x^{*}_{t}\widetilde{\alpha}_{t}^{j}-\alpha^{j}=\sum_{t=1}^{k}x^{*}_{t}\sum_{p=0}^{j}\binom{j}{p}\Delta\widetilde{\alpha}_{t}^{p}\alpha^{j-p}-\alpha^{j}$
$\displaystyle=\sum_{p=1}^{j}\binom{j}{p}\alpha^{j-p}\sum_{t=1}^{k}x_{t}^{*}\Delta\widetilde{\alpha}_{t}^{p}.$

When $1\leq j<k$, clear we have
$\sum_{t=1}^{k}x_{t}^{*}\Delta\widetilde{\alpha}_{t}^{j}=0.$

Now, we bound $\sum_{t=1}^{k}x_{t}^{*}\Delta\widetilde{\alpha}_{t}^{j}$ for
$j\geq k$. Considering the relation between the inverse matrix the
corresponding adjoint matrix, we have for example

$\displaystyle
x_{1}^{*}=\left|\begin{matrix}\Delta\widetilde{\alpha}_{2}^{1}&\cdots&\Delta\widetilde{\alpha}_{k}^{1}\\\
\vdots&\ddots&\vdots\\\
\Delta\widetilde{\alpha}_{2}^{k-1}&\cdots&\Delta\widetilde{\alpha}_{k}^{k-1}\end{matrix}\right|\cdot\left|\begin{matrix}1&1&\cdots&1\\\
\Delta\widetilde{\alpha}_{1}^{1}&\Delta\widetilde{\alpha}_{2}^{1}&\cdots&\Delta\widetilde{\alpha}_{k}^{1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\Delta\widetilde{\alpha}_{1}^{k-1}&\Delta\widetilde{\alpha}_{2}^{k-1}&\cdots&\Delta\widetilde{\alpha}_{k}^{k-1}\end{matrix}\right|^{-1}.$

Thus, for $j\geq k$, we have that

$\displaystyle\sum_{t=1}^{k}x_{t}^{*}\Delta\widetilde{\alpha}_{t}^{j}$
$\displaystyle=\left|\begin{matrix}\Delta\widetilde{\alpha}_{1}^{j}&\Delta\widetilde{\alpha}_{2}^{j}&\cdots&\Delta\widetilde{\alpha}_{k}^{j}\\\
\Delta\widetilde{\alpha}_{1}^{1}&\Delta\widetilde{\alpha}_{2}^{1}&\cdots&\Delta\widetilde{\alpha}_{k}^{1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\Delta\widetilde{\alpha}_{1}^{k-1}&\Delta\widetilde{\alpha}_{2}^{k-1}&\cdots&\Delta\widetilde{\alpha}_{k}^{k-1}\end{matrix}\right|\cdot\left|\begin{matrix}1&1&\cdots&1\\\
\Delta\widetilde{\alpha}_{1}^{1}&\Delta\widetilde{\alpha}_{2}^{1}&\cdots&\Delta\widetilde{\alpha}_{k}^{1}\\\
\vdots&\vdots&\ddots&\vdots\\\
\Delta\widetilde{\alpha}_{1}^{k-1}&\Delta\widetilde{\alpha}_{2}^{k-1}&\cdots&\Delta\widetilde{\alpha}_{k}^{k-1}\end{matrix}\right|^{-1}$
$\displaystyle=\prod_{t=1}^{k}\Delta\widetilde{\alpha}_{t}\cdot\sum_{s\subseteq(k)^{j-k}}\prod_{t=1}^{k}\Delta\widetilde{\alpha}_{t}^{s_{t}}$

where the last equation holds due to Lemma 3.8. According to
$|\Delta\widetilde{\alpha}_{t}|\leq|\alpha|+|\widetilde{\alpha}_{t}|\leq 2$
and $(k)^{j-k}$ has $\binom{j}{k}$ terms, we have

$\displaystyle\left|\sum_{t=1}^{k}x_{t}^{*}\Delta\widetilde{\alpha}_{t}^{j}\right|\leq
2^{k}\cdot 2^{j-k}\prod_{t=1}^{k}|\Delta\widetilde{\alpha}_{t}|.$

Plugging in this result, we can see that

$\displaystyle|(V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}^{*}-M(\alpha))_{j}|$
$\displaystyle=\left|\sum_{p=1}^{j}\binom{j}{p}\alpha^{j-p}\sum_{t=1}^{k}x_{t}^{*}\Delta\widetilde{\alpha}_{t}^{p}\right|\leq\sum_{p=k}^{j}\binom{j}{p}\cdot
2^{k}\cdot 2^{p-k}\prod_{t=1}^{k}|\Delta\widetilde{\alpha}_{t}|$
$\displaystyle\leq 2^{k}\cdot
2^{2j}\prod_{t=1}^{k}|\Delta\widetilde{\alpha}_{t}|.$

Therefore, we have that

$\displaystyle\min_{\boldsymbol{x}\in
F^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M(\alpha)\|_{2}$
$\displaystyle\leq\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}^{*}-M(\alpha)\|_{2}\leq\sqrt{\sum_{j=k}^{2k-1}\left(2^{k}\cdot
2^{2j}\prod_{t=1}^{k}|\Delta\widetilde{\alpha}_{t}|\right)^{2}}\leq
2^{O(k)}\prod_{j=1}^{k}|\alpha-\widetilde{\alpha}_{j}|.$

∎

Lemma 3.5 shows that every discrete distribution in ground truth
$\boldsymbol{\vartheta}$ has a close spike in recovered positions
$\widetilde{\boldsymbol{\alpha}}$ and thus can be well approximated according
to the above lemma. Therefore, the mixture of the discrete distributions can
also be well approximated with support $\widetilde{\boldsymbol{\alpha}}$. As
the following lemma shows, solving linear regression in Line 6 finds a weight
vector $\widetilde{w}$ which the corresponding moment vector
$V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}$ is close to the
ground truth $M(\boldsymbol{\vartheta})$.

######  Lemma 3.10.

Let
$\widehat{\boldsymbol{w}}=[\widehat{w}_{1},\cdots,\widehat{w}_{k}]^{\top}\in\mathbb{R}^{k}$
be the intermediate result (Line 6) in Algorithm 1. Then,

$\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{\infty}\leq
2^{O(k)}\cdot\xi^{\frac{1}{2}}.$

###### Proof.

Firstly, by triangle inequality, we can see

$\displaystyle\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{2}$
$\displaystyle\leq\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M^{\prime}\|_{2}+\|M(\boldsymbol{\vartheta})-M^{\prime}\|_{2}$
$\displaystyle\leq
O(1)\cdot\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M^{\prime}\|_{2}+\|M(\boldsymbol{\vartheta})-M^{\prime}\|_{2}$
$\displaystyle\leq
O(1)\cdot\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M(\boldsymbol{\vartheta})\|_{2}+O(1)\cdot\|M(\boldsymbol{\vartheta})-M^{\prime}\|_{2}$

where the first and third inequalities hold due to triangle inequality, and
the second inequality holds since
$\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M^{\prime}\|_{2}^{2}\leq
O(1)\cdot\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M^{\prime}\|_{2}^{2}$
(we only find an $O(1)$-approximation in this step).

For the first term, we can see that

$\displaystyle\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M(\boldsymbol{\vartheta})\|_{2}$
$\displaystyle=\min_{\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{k}\in\mathbb{R}^{k}}\left\|\sum_{i=1}^{k}w_{i}V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}_{i}-\sum_{i=1}^{k}w_{i}M(\alpha_{i})\right\|_{2}$
$\displaystyle\leq\sum_{i=1}^{k}w_{i}\min_{\boldsymbol{x}\in\mathbb{R}^{k}}\|V_{\widetilde{\boldsymbol{\alpha}}}\boldsymbol{x}-M(\alpha_{i})\|_{2}$
$\displaystyle\leq
2^{O(k)}\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widetilde{\alpha}_{j}|$
$\displaystyle(\text{Lemma~{}\ref{lm:step3}})$ $\displaystyle\leq
2^{O(k)}\sqrt{\sum_{i=1}^{k}w_{i}\prod_{j=1}^{k}|\alpha_{i}-\widetilde{\alpha}_{j}|^{2}}$
$\displaystyle(\text{AM-QM Inequality},\boldsymbol{w}\in\Delta_{k-1})$
$\displaystyle\leq 2^{O(k)}\cdot\sqrt{\xi}.$
$\displaystyle(\text{Lemma~{}\ref{lm:dim1-step2}}).$

For the second term,
$\|M(\boldsymbol{\vartheta})-M^{\prime}\|_{2}\leq\sqrt{k}\|M(\boldsymbol{\vartheta})-M^{\prime}\|_{\infty}\leq\sqrt{k}\cdot\xi$.
As a result,

$\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{\infty}\leq\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{2}\leq
O(1)\cdot 2^{O(k)}\cdot\sqrt{\xi}+O(1)\cdot\sqrt{k}\cdot\xi\leq
2^{O(k)}\cdot\zeta^{\frac{1}{2}}.$

∎

The moment-transportation inequality (Lemma B.1) requires the input to be
signed measure. So, we normalize the recovered spikes in Line 7. The following
lemma shows the moment vector would not change too much after normalization.

######  Lemma 3.11.

Let
$\widetilde{\boldsymbol{\vartheta}}=(\widetilde{\boldsymbol{\alpha}},\widetilde{\boldsymbol{w}})$
be the intermediate result (Line 8) in Algorithm 1. Then, we have

$\operatorname{Mom}_{K}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
2^{O(k)}\cdot\xi^{\frac{1}{2}}.$

###### Proof.

According to the definition of $V_{\widetilde{\boldsymbol{\alpha}}}$ and
$M(\boldsymbol{\vartheta})$, we have
$(V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta}))_{1}=\sum_{i=1}^{k}\widehat{w}_{i}-1$.
Therefore,
$|\sum_{i=1}^{k}\widehat{w}_{i}-1|\leq|(V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta}))_{1}|\leq\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{\infty}\leq
2^{O(k)}\cdot\xi^{\frac{1}{2}}.$

Note that
$\|\widetilde{\boldsymbol{w}}-\widehat{\boldsymbol{w}}\|_{1}=\|\widetilde{\boldsymbol{w}}\|_{1}-\|\widehat{\boldsymbol{w}}\|_{1}=|(\sum_{i=1}^{k}\widehat{w}_{i})^{-1}-1|$.
For $2^{O(k)}\cdot\xi^{\frac{1}{2}}\leq 1/2$, we can conclude that
$\|\widetilde{\boldsymbol{w}}-\widehat{\boldsymbol{w}}\|_{1}=|(\sum_{i=1}^{k}\widehat{w}_{i})^{-1}-1|\leq
2|\sum_{i=1}^{k}\widehat{w}_{i}-1|\leq 2^{O(k)}\cdot\xi^{\frac{1}{2}}.$

Moreover,

$\displaystyle\operatorname{Mom}_{K}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})$
$\displaystyle=\|M(\widetilde{\boldsymbol{\vartheta}})-M(\boldsymbol{\vartheta})\|_{\infty}=\|V_{\widetilde{\boldsymbol{\alpha}}}\widetilde{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{\infty}$
$\displaystyle\leq\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{\infty}+\|V_{\widetilde{\boldsymbol{\alpha}}}\widetilde{\boldsymbol{w}}-V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}})\|_{\infty}$
$\displaystyle\leq\|V_{\widetilde{\boldsymbol{\alpha}}}\widehat{\boldsymbol{w}}-M(\boldsymbol{\vartheta})\|_{\infty}+\|V_{\widetilde{\boldsymbol{\alpha}}}\|_{\infty}\|\widetilde{\boldsymbol{w}}-\widehat{\boldsymbol{w}}\|_{1}$
$\displaystyle\leq 2^{O(k)}\cdot\xi^{\frac{1}{2}}.$

where the first inequality holds due to triangle inequality and third holds
since $\|V_{\widetilde{\boldsymbol{\alpha}}}\|_{\infty}\leq 1$. ∎

Moreover, the weight can still be negative after normalization. To find a
reconstruction in original space, we want to find a mixture of discrete
distributions in $\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$ that is close to
$\widetilde{\boldsymbol{\vartheta}}$. However, this step can cast a huge
impact on the moments. Thus, we directly estimate the influence in terms of
the transportation distance instead of the moment distance. We note the
transportation distance in the next lemma is defined possibly for non-
probability measures (see Equation (6) in Appendix A for details).

######  Lemma 3.12.

Let $\mathcal{C}$ be a compact convex set in $\mathbb{R}^{d}$ and let
$\widetilde{\boldsymbol{\vartheta}}=(\widetilde{\boldsymbol{\alpha}},\widetilde{\boldsymbol{w}})\in\mathrm{Spike}(\mathcal{C},\Sigma_{k-1})$
be a $k$-spike distribution over support $\mathcal{C}$. Then, we have

$\min_{\widecheck{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\widetilde{\boldsymbol{\alpha}},\Delta_{k-1})}\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\widecheck{\boldsymbol{\vartheta}})\leq
2\min_{\overline{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\mathcal{C},\Delta_{k-1})}\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\overline{\boldsymbol{\vartheta}}).$

Note that the minimization on the left is over support
$\widetilde{\boldsymbol{\alpha}}$.

The proof of the above lemma can be found in Appendix C.2. Finally, we are
ready to bound the transportation distance error of the reconstructed
$k$-spike distribution:

######  Lemma 3.13.

Let
$\widecheck{\boldsymbol{\vartheta}}=(\widetilde{\boldsymbol{\alpha}},\widecheck{\boldsymbol{w}})$
be the final result (Line 10) in Algorithm 1. Then, we have

$\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(k\xi^{\frac{1}{4k-2}}).$

###### Proof.

Combining Lemma B.1 and Lemma 3.11, we can see

$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(k\cdot(2^{O(k)}\cdot\xi^{\frac{1}{2}})^{\frac{1}{2k-1}})\leq
O(k\xi^{\frac{1}{4k-2}}).$

Moreover, from Lemma 3.12, we have

$\displaystyle\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\widecheck{\boldsymbol{\vartheta}})\leq
2\min_{\overline{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\Delta_{1},\Delta_{k-1})}\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\overline{\boldsymbol{\vartheta}})\leq
2\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})=O(k\xi^{\frac{1}{4k-2}}).$

Finally, by triangle inequality,
$\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\widetilde{\boldsymbol{\vartheta}})+\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(k\xi^{\frac{1}{4k-2}}).$ ∎

By choosing $\xi=(\epsilon/k)^{O(k)}$, the previous lemma directly implies
Theorem 1.1.

##  4 Efficient Algorithms for Higher Dimensions

In this section, we solve the problem in higher dimensions. We first present
the algorithm for 2-dimension. The algorithm for higher dimensions use the
algorithms for both 1- and 2-dimension as subroutines.

###  4.1 An Efficient Algorithm for the 2-Dimensional Problem

We can generalize the 1-dimension algorithm described in Section 3 to
2-dimension. Let
$\boldsymbol{\vartheta}:=(\boldsymbol{\alpha},\boldsymbol{w})\in\mathrm{Spike}(\Delta_{2},\Delta_{k-1})$
be the underlying mixture where
$\boldsymbol{\alpha}=[\boldsymbol{\alpha}_{1},\cdots,\boldsymbol{\alpha}_{k}]^{\top}$
for $\boldsymbol{\alpha}_{i}=(\alpha_{i,1},\alpha_{i,2})$ and
$\boldsymbol{w}=[w_{1},\cdots,w_{k}]^{\top}\in\Delta_{k-1}$. The true moments
can be computed according to
$M_{i,j}(\boldsymbol{\vartheta})=\sum_{t=1}^{k}w_{t}\alpha_{t,1}^{i}\alpha_{t,2}^{j}$.
The input is the noisy moments $M^{\prime}_{i,j}$ in $0\leq i,j,i+j\leq 2k-1$
such that $|M^{\prime}_{i,j}-M_{i,j}(\boldsymbol{\vartheta})|\leq\xi$. We
further assume that $M^{\prime}_{0,0}=M_{0,0}(\boldsymbol{\vartheta})=1$ and
$\xi\leq 2^{-\Omega(k)}$.

The key idea is very simple: a distribution supported in $\mathbb{R}^{2}$ can
be mapped to a distribution supported in the complex plane $\mathbb{C}$. In
particular, we define the complex set
$\Delta_{\mathbb{C}}=\\{a+b{\mathfrak{i}}\mid(a,b)\in\Delta_{2}\\}$. Moreover,
we denote
$\boldsymbol{\beta}=[\beta_{1},\cdots,\beta_{k}]^{\top}:=[\alpha_{1,1}+\alpha_{1,2}{\mathfrak{i}},\cdots,\alpha_{k,1}+\alpha_{k,2}{\mathfrak{i}}]^{\top}\in\Delta_{\mathbb{C}}^{k}$,
and define
$\boldsymbol{\phi}:=(\boldsymbol{\beta},\boldsymbol{w})\in\mathrm{Spike}(\Delta_{\mathbb{C}}^{k},\Delta_{k-1})$
to be the complex mixture corresponding to $\boldsymbol{\vartheta}$. The
corresponding moments of $\boldsymbol{\phi}$ can thus be defined as
$G_{i,j}(\boldsymbol{\phi})=\sum_{t=1}^{k}w_{t}(\beta_{t}^{\dagger})^{i}\beta_{t}^{j}$,
where $\beta_{t}^{\dagger}$ is the complex conjugate of $\beta_{t}$. In fact
$G_{i,j}$ can be computed from $M_{i,j}$ (see Lemma D.1). Therefore, our task
reduces to the recovery the complex mixture $\boldsymbol{\phi}$ given noisy
moments $G^{\prime}_{i,j}$. The algorithm and the proof of the following
theorem are very similar to that in Section 3. Handling complex numbers
requires several minor changes and we present the full details in Appendix D.
For the analysis, we need to generalize the moment-transportation inequality
to complex domain as well (see Appendix B.2).

######  Theorem 4.1.

(The 2-dimensional Problem) Let
$\boldsymbol{\vartheta}\in\mathrm{Spike}(\Delta_{2},\Delta_{k-1})$. Suppose we
have noisy moments $M_{i,j}$ ($0\leq i,j,i+j\leq 2k-1$) up to precision
$(\epsilon/k)^{O(k)}$. We can obtain a mixture
$\widetilde{\boldsymbol{\vartheta}}$ such that
$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(\epsilon)$ and our algorithm runs in time $O(k^{3})$.

###  4.2 An Efficient Algorithm for $d>3$

Now, we present our algorithm for higher dimensions. Let
$\boldsymbol{\vartheta}:=(\boldsymbol{\alpha},\boldsymbol{w})\in\mathrm{Spike}(\Delta_{d-1},\Delta_{k-1})$
be the underlying mixture where
$\boldsymbol{\alpha}=[\boldsymbol{\alpha}_{1},\cdots,\boldsymbol{\alpha}_{k}]^{\top}\in\mathbb{R}^{k\times
d}$ such that
$\boldsymbol{\alpha}_{i}=[\alpha_{i,1},\cdots,\alpha_{i,d}]^{\top}$ and
$\boldsymbol{w}=[w_{1},\cdots,w_{k}]^{\top}$. We denote
$\boldsymbol{\alpha}_{:,i}=[\alpha_{1,i},\cdots,\alpha_{k,i}]$ as the vector
of the $i$th coordinates of the spikes. We write
$\widetilde{\boldsymbol{\alpha}}=[\boldsymbol{\alpha}_{:,1},\cdots,\boldsymbol{\alpha}_{:,k}]^{\top}$.

Since the number of different moments is exponential in $d$, we consider the
setting in which one can access the noisy moments of some linear maps of
$\boldsymbol{\vartheta}$ onto lower dimensional subspaces, as in previous work
[55, 45] (such noisy moments can be easily obtained in applications such as
topic modeling, see Section 5, or learning Gaussian mixtures, see Section 6).
For some
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})\in\mathrm{Spike}(\mathbb{R}^{d},\Delta_{k-1})$
and a vector $\boldsymbol{r}\in\mathbb{R}^{d}$, we denote
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}):=(\boldsymbol{\alpha}\boldsymbol{r},\boldsymbol{w})\in\mathrm{Spike}(\mathbb{R},\Delta_{k-1})$
as the projected distribution of $\boldsymbol{\vartheta}$ along vector
$\boldsymbol{r}$. For a matrix
$R=[\boldsymbol{r}_{1},\cdots,\boldsymbol{r}_{p}]\in\mathbb{R}^{d\times p}$,
we also denote
$\Pi_{R}(\boldsymbol{\vartheta}):=(\boldsymbol{\alpha}R,\boldsymbol{w})\in\mathrm{Spike}(\mathbb{R}^{p},\Delta_{k-1})$
as the projected distribution of $\boldsymbol{\vartheta}$ along multiple
dimensions $\boldsymbol{r}_{1},\cdots,\boldsymbol{r}_{p}$. In this section, we
consider mappings to 1-dimensional and 2-dimensional subspaces, which suffice
for our purposes.

We assume the moment information is revealed by a noisy moment oracle. A
moment oracle corresponding to $k$-spike distribution $\boldsymbol{\vartheta}$
is defined to be a function $M^{\prime}(\cdot)$ that takes a matrix (or a
vector) $R$ with $\|R\|_{\infty}\leq 1$ and generates a noisy vector
$M^{\prime}(R)$ such that
$\|M^{\prime}(R)-M(\Pi_{R}(\boldsymbol{\vartheta}))\|_{\infty}\leq\xi$ and
$M^{\prime}_{\mathbf{0}}(R)=M_{\mathbf{0}}(\Pi_{R}(\boldsymbol{\vartheta}))=1$
where $\xi\leq 2^{-\Omega(k)}$. Recall that the moment vector $M$ is defined
in (1).

Let $e_{t}\in\mathbb{R}^{d}$ be the unit vector that takes value $1$ only
$t$-th coordinates and $0$ in the remaining coordinates. Denote
$\mathbf{1_{d}}=(1,\cdots,1)\in\mathbb{R}^{d}$. Let $\mathbb{S}_{d-1}$ be the
unit sphere in $\mathbb{R}^{d}$.

Algorithm 2 Reconstruction Algorithm in High Dimension

1:function HighDimension($k,M^{\prime}(\cdot),\xi$)

2: $\boldsymbol{r}\leftarrow\mathbb{S}_{d-1}$

3:
$\widetilde{\boldsymbol{\phi}}^{\prime}\leftarrow\textrm{OneDimension}(k,M^{\prime}(\frac{\boldsymbol{r}+\mathbf{1}}{4}),\xi)$
$\triangleright$
$\widetilde{\boldsymbol{\phi}}^{\prime}=(\widetilde{\boldsymbol{y}}^{\prime},\widetilde{\boldsymbol{w}}^{\prime})\in\mathrm{Spike}(\Delta_{1},\Delta_{k-1})$

4:
$\widetilde{\boldsymbol{\phi}}\leftarrow([4\widetilde{\boldsymbol{y}}^{\prime}-\mathbf{1}_{d}],\widetilde{\boldsymbol{w}}^{\prime})$
$\triangleright$
$\widetilde{\boldsymbol{\phi}}=(\widetilde{\boldsymbol{y}},\widetilde{\boldsymbol{w}})$

5: for $t\in[1,d]$ do

6:
$\widehat{\boldsymbol{\phi}}_{t}^{\prime}\leftarrow\textrm{TwoDimension}(k,M^{\prime}([\frac{\boldsymbol{r}+\mathbf{1}}{4},\frac{e_{t}}{2}]),\xi)$
$\triangleright$
$\widehat{\boldsymbol{\phi}}_{t}^{\prime}=([\widehat{\boldsymbol{y}}_{:,t}^{\prime},\widehat{\boldsymbol{\alpha}}_{:,t}^{\prime}],\widehat{\boldsymbol{w}}_{:,t}^{\prime})\in\mathrm{Spike}(\Delta_{2},\Delta_{k-1})$

7:
$\widehat{\boldsymbol{\phi}}_{t}\leftarrow([4\boldsymbol{y}_{:,t}^{\prime}-\mathbf{1}_{d},2\widehat{\boldsymbol{\alpha}}_{:,t}^{\prime}],\widehat{\boldsymbol{w}}_{:,t}^{\prime})$
$\triangleright$
$\widehat{\boldsymbol{\phi}}_{t}=([\widehat{\boldsymbol{y}}_{:,t},\widehat{\boldsymbol{\alpha}}_{:,t}],\widehat{\boldsymbol{w}}_{:,t})$

8: for $j\in[1,k]$ do

9:
$s^{*}_{j,t}\leftarrow\arg\min_{s:\widehat{w}_{s,t}\geq\frac{\sqrt{\epsilon}}{k}}|\widetilde{y}_{j}-\widehat{y}_{s,t}|$

10: $\widetilde{\alpha}_{j,t}\leftarrow\widehat{\alpha}_{s_{j,t}^{*},t}$

11: end for

12:
$\widetilde{\boldsymbol{\alpha}}_{:,t}\leftarrow[\widetilde{\alpha}_{1,t},\cdots,\widetilde{\alpha}_{k,t}]^{\top}$
$\triangleright$ $\widetilde{\boldsymbol{\alpha}}\in\mathbb{R}^{k}$

13: end for

14:
$\widetilde{\boldsymbol{\alpha}}\leftarrow[\widetilde{\boldsymbol{\alpha}}_{:,1},\cdots,\widetilde{\boldsymbol{\alpha}}_{:,d}]^{\top}$
$\triangleright$
$\widetilde{\boldsymbol{\alpha}}=[\widetilde{\boldsymbol{\alpha}}_{1},\cdots,\widetilde{\boldsymbol{\alpha}}_{k}]\in\mathbb{R}^{d\times
k}$

15:
$\widetilde{\boldsymbol{\vartheta}}\leftarrow(\widetilde{\boldsymbol{\alpha}},\widetilde{\boldsymbol{w}})$
$\triangleright$
$\widetilde{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\mathbb{R}^{d},\Delta_{k-1})$

16:
$\widecheck{\boldsymbol{\alpha}}\leftarrow\textrm{project}_{\Delta_{d-1}}(\widetilde{\boldsymbol{\alpha}})$
$\triangleright$
$\widecheck{\boldsymbol{\alpha}}=[\widecheck{\boldsymbol{\alpha}}_{1},\cdots,\widecheck{\boldsymbol{\alpha}}_{k}]\in\Delta_{d-1}^{k}$

17: report
$\widecheck{\boldsymbol{\vartheta}}\leftarrow(\widecheck{\boldsymbol{\alpha}},\widetilde{\boldsymbol{w}})$
$\triangleright$
$\widecheck{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\Delta_{d-1},\Delta_{k-1})$

18:end function

The pseudocode can be found in Algorithm 2. The input of the algorithm
consists of the number of spikes $k$, the moment error bound $\xi$ and the
noisy moment oracle $M^{\prime}(\cdot)$ mentioned above. Now, we describe the
implementation details of our algorithm.

We first generate a random vector $\boldsymbol{r}$ in Line 2 by sampling from
unit sphere
$\mathbb{S}_{d-1}=\\{\boldsymbol{r}\in\mathbb{R}^{d}\mid\|\boldsymbol{r}\|_{2}=1\\}$.
Note that by a probabilistic argument, the distance between spikes are roughly
kept after the projection along $\boldsymbol{r}$ (see Lemma 4.2). Then, we aim
to use the one dimension algorithm to recover
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$. However, the support of
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$ is contained in $[-1,1]$ but
not $[0,1]$. In Line 3, we apply the 1-dimension algorithm (Algorithm 1) to
the noisy moments of a shifted and scaled map
$\Pi_{(\boldsymbol{r}+\mathbf{1}_{d})/4}(\boldsymbol{\vartheta})$ whose
support is in $[0,1]$. Then we scale and shift back to obtain the result
$\widetilde{\boldsymbol{\phi}}$ in Line 4. Intuitively,
$\widetilde{\boldsymbol{\phi}}$ is close to
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$.

Now, we try to recover the coordinates of the spikes. In particular, we would
like to find, for each spike of $\widetilde{\boldsymbol{\phi}}$, the
coordinates of the corresponding location in $\Delta_{d-1}$. This is done by
considering each dimension separately. For each coordinate $t\in[d]$, we run
the 2-dimension (Algorithm 3) to recovery the projection
$\Pi_{[\boldsymbol{r},e_{t}]}(\boldsymbol{\vartheta})$, that is a linear map
of $\boldsymbol{\vartheta}$ onto the subspace spanned by vector
$\boldsymbol{r}$ and axis $e_{t}$. However,
$[\boldsymbol{r},e_{t}]\boldsymbol{\alpha}_{i}$ is in $[-1,1]\times[0,1]$.
Again, in Line 6, we apply the algorithm to the noisy moments of a shifted and
scaled projection $\Pi_{B_{t}}(\boldsymbol{\vartheta})$ where
$B_{t}=[\frac{\boldsymbol{r}+\mathbf{1}_{d}}{4},\frac{e_{t}}{2}]$, and scale
the result back in Line 7. The result is denoted as
$\widehat{\boldsymbol{\phi}}_{t}$.

Next, we assemble the 1-d projection $\widetilde{\boldsymbol{\phi}}$ and 2-d
projections $\widehat{\boldsymbol{\phi}}_{t},t\in[d]$. Due to the noise,
$\widetilde{\boldsymbol{\phi}}_{t}$ and $\widetilde{\boldsymbol{\phi}}$ may
have different supports along $\boldsymbol{r}$. So we assign the coordinates
of $\widetilde{\boldsymbol{\phi}}$ according to the closest spike of
$\widehat{\boldsymbol{\phi}}_{t}$ with a relative large weights in Line 8 to
Line 12. Due to the error, the support of the reconstructed distribution
$\widetilde{\boldsymbol{\vartheta}}$ may not locate in $\Delta_{d-1}$. Thus
the final step in Line 16 is to project each individual spike back to
$\Delta_{d-1}$.

We note the bottleneck of this algorithm is Line 6 where we have to solve the
2-dimensional problem $d$ times. Hence, the total running time of the
algorithm is $O(dk^{3})$.

###  4.3 Error Analysis

First, we show the random vector $\boldsymbol{r}$ has a critical property that
two distant spikes are still distant after the projection with high
probability. This result is standard and similar statements are known in the
literature (e.g., Lemma 37 in [61]).

######  Lemma 4.2.

Let
$\mathrm{Supp}=(\boldsymbol{\alpha}_{1},\cdots,\boldsymbol{\alpha}_{k})\subset\Delta_{d-1}$
such that $d>3$. For every constant $0<\eta<1$ and a random vector
$\boldsymbol{r}\leftarrow\mathbb{S}_{d-1}$, with probability at least
$1-\eta$, for any
$\boldsymbol{\alpha}_{i},\boldsymbol{\alpha}_{j}\in\mathrm{Supp}$, we have
that

$|\boldsymbol{r}^{\top}(\boldsymbol{\alpha}_{i}-\boldsymbol{\alpha}_{j})|\geq\frac{\eta}{k^{2}d}\|\boldsymbol{\alpha}_{i}-\boldsymbol{\alpha}_{j}\|_{1}.$

###### Proof.

Let $0<\epsilon<1$ be some constant depending on $k$ and $d$ which we fix
later. Assume $\boldsymbol{r}=[r_{1},\cdots,r_{d}]^{\top}$. For a fixed pair
$(\boldsymbol{\alpha}_{i},\boldsymbol{\alpha}_{j})$, we have

$\Pr_{\boldsymbol{r}\sim\mathbb{S}_{d-1}}[|\boldsymbol{r}^{\top}(\boldsymbol{\alpha}_{i}-\boldsymbol{\alpha}_{j})|<\epsilon\|\boldsymbol{\alpha}_{i}-\boldsymbol{\alpha}_{j}\|_{2}]=\Pr_{\boldsymbol{r}\sim\mathbb{S}_{d-1}}[|r_{1}|<\epsilon]$

Suppose $\boldsymbol{r}$ is obtained by first sampling a vector from unit ball
$\mathbb{B}_{d}$ and then normalizing it to the unit sphere
$\mathbb{S}_{d-1}$. We can see that

$\displaystyle\Pr_{\boldsymbol{r}\sim\mathbb{S}_{d-1}}[|r_{1}|<\epsilon]\leq\Pr_{\boldsymbol{r}\sim\mathbb{B}_{d}}[|r_{1}|<\epsilon]$

where the right hand side is equal to the probability that $\boldsymbol{r}$
lies in the slice $[-\epsilon,\epsilon]$.

Let $V_{d}$ be the volume of $d$-dimensional ball. We have

$\displaystyle\Pr_{\boldsymbol{r}\sim\mathbb{B}_{d}}[|r_{1}|<\epsilon]$
$\displaystyle=\frac{\int_{-\epsilon}^{\epsilon}V_{d-1}(1-x^{2})^{\frac{d-1}{2}}\mathrm{d}x}{V_{d}}\leq\epsilon\sqrt{d}$

where the last inequality holds since $V_{d}/V_{d-1}\geq 2/\sqrt{d}$.

By union bound over all pairs
$(\boldsymbol{\alpha}_{i},\boldsymbol{\alpha}_{j})$, we can see that the
failure probability can be bounded by

$\Pr_{\boldsymbol{r}\sim\mathbb{S}_{d-1}}[\exists\boldsymbol{\alpha}_{i},\boldsymbol{\alpha}_{j}\in\mathrm{Supp}:|\boldsymbol{r}^{\top}(A_{i}-A_{j})|<\epsilon\|A_{i}-A_{j}\|_{2}]<k^{2}\epsilon\sqrt{d}.$

Notice that
$\|\boldsymbol{\alpha}\|_{1}\leq\sqrt{d}\|\boldsymbol{\alpha}\|_{2}$ for any
$\boldsymbol{\alpha}\in\mathbb{R}^{d}$, and let
$\epsilon=\frac{\eta}{k^{2}\sqrt{d}}$. We conclude that

$\Pr_{\boldsymbol{r}\sim\mathbb{S}_{d-1}}[\forall\boldsymbol{\alpha}_{i},\boldsymbol{\alpha}_{j}\in\mathrm{Supp}:|\boldsymbol{r}^{\top}(A_{i}-A_{j})|\geq\frac{\eta}{k^{2}d}\|A_{i}-A_{j}\|_{1}]>1-\eta.$

∎

Next, we show that we can recover reconstruct the higher dimension mixture
$\boldsymbol{\vartheta}$, using its 1-d projection
$\widetilde{\boldsymbol{\phi}}\sim\Pi_{[\boldsymbol{r}]}(\boldsymbol{\vartheta})$,
by assigning each spike $(\widetilde{y}_{i},\widetilde{w}_{i})$ in the one
dimension distribution a location in $\mathbb{R}^{d}$. In a high level,
consider a clustering of the spikes in
$\Pi_{[\boldsymbol{r}]}(\boldsymbol{\vartheta})$. The spikes in one cluster
still form a cluster in the original $\boldsymbol{\vartheta}$ by Lemma 4.2
since the projection approximately keeps the distance. Thus assigning the
clusters a location in $\mathbb{R}^{d}$ can produce to a good estimation for
$\boldsymbol{\vartheta}$. The following lemma formally proves that we can
reconstruct the original mixture using 2-dimension projections.

######  Lemma 4.3.

Let
$\widetilde{\boldsymbol{\vartheta}}\in\mathrm{Spike}(\mathbb{R}^{d},\Delta_{k-1})$
be the intermediate result (Line 15) in Algorithm 2. Let $\epsilon$ be the
smallest number such that
$\operatorname{Tran}(\widetilde{\boldsymbol{\phi}},\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))\leq\epsilon$
and
$\operatorname{Tran}(\widehat{\boldsymbol{\phi}}_{t},\Pi_{[\boldsymbol{r},e_{t}]}(\boldsymbol{\vartheta}))\leq\epsilon$
for all $t$. If Lemma 4.2 holds for $\boldsymbol{\vartheta}$, we have that

$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(k^{3}d^{2}\sqrt{\epsilon}).$

###### Proof.

Partition the spikes of $\boldsymbol{\vartheta}$ into clusters
$C_{1},\cdots,C_{m}$ such that

  * •

$|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{j}|\leq
k\sqrt{\epsilon}$ if $\boldsymbol{\alpha}_{i}$ and $\boldsymbol{\alpha}_{j}$
are in the same cluster $C_{p}$.

  * •

$|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{j}|>\sqrt{\epsilon}$
if $\boldsymbol{\alpha}_{i}$ and $\boldsymbol{\alpha}_{j}$ are in different
clusters.

We define the distance between spike $\widetilde{y}_{j}$ and cluster $C_{p}$
to be $\min_{i\in
C_{p}}|\widetilde{y}_{j}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}|$.
Clearly, there exists at most one cluster $C_{p}$ such that the distance
between the spike and the cluster is less than $\frac{\sqrt{\epsilon}}{2}$. In
this case,
$|\widetilde{y}_{j}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}|>\frac{\sqrt{\epsilon}}{2}$
holds for all $j\notin C_{p}$. Let $\widetilde{C}_{p}$ be the set of all index
$j$ such that the distance between spike $\widetilde{y}_{j}$ and cluster
$C_{p}$ is less than $\frac{\sqrt{\epsilon}}{2}$. From this definition, the
distance between the spikes in $\widetilde{C}_{p}$ and the spikes not in
$C_{p}$ is at least $\frac{\sqrt{\epsilon}}{2}$. By the performance guarantee
of our 1-dimensional algorithm, the total transportation distance between
$\widetilde{\boldsymbol{\phi}}$ and
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$ is no more than $\epsilon$.
Hence, if we consider couple $C_{p}$ and $\widetilde{C}_{p}$ together, we can
see that

$|\sum_{i\in
C_{p}}w_{i}-\sum_{j\in\widetilde{C}_{p}}\widetilde{w}_{j}|<2\sqrt{\epsilon}.$

Recall the transportation distance
$\operatorname{Tran}(\boldsymbol{\vartheta},\widetilde{\boldsymbol{\vartheta}})$
can be also defined as the solution of the following program over some joint
distribution $\mu$ (i.e., coupling).

min
$\displaystyle\sum_{i=1}^{k}\sum_{j=1}^{k}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}$
s.t. $\displaystyle
w_{i}=\sum_{j=1}^{k}\mu_{i,j},\widetilde{w}_{j}=\sum_{i=1}^{k}\mu_{i,j},\mu\geq\mathbf{0}.$

Intuitively, a near-optimal way for generating such distribution is to first
couple spikes in the same cluster with as much weight as possible, then couple
spikes between different clusters. Formally, consider a two-stage procedure
that construct $\mu$. In the first stage, increase $\mu_{i,j}$ as much as
possible without exceeding the boundary for $(i,j)$ that $i\in C_{p}$ and
$j\in\widetilde{C}_{p}$ for some common $p$. In the second stage, fill
$\mu_{i,j}$ to fulfill all constraints. Since the restriction set is convex,
$C_{p}$ are disjoint, and $\widetilde{C}_{p}$ are disjoint, this two-stage
process ensures

$\displaystyle\sum_{i\in
C_{p}}\sum_{j\in\widetilde{C}_{p}}\mu_{i,j}=\min\\{\sum_{i\in
C_{p}}w_{i},\sum_{j\in\widetilde{C}_{p}}\widetilde{w}_{j}\\}.$

As a result, the coupling distribution satisfies

$\displaystyle\sum_{i\in
C_{p}}\sum_{j\notin\widetilde{C}_{p}}\mu_{i,j}\leq|\sum_{i\in
C_{p}}w_{i}-\sum_{j\in\widetilde{C}_{p}}\widetilde{w}_{j}|\leq
2\sqrt{\epsilon}.$

Now we begin to show $\mu$ is a good joint distribution for
$\operatorname{Tran}(\boldsymbol{\vartheta},\widetilde{\boldsymbol{\vartheta}})$.
Since $\\{C_{p}\\}$ is a partition of $[1,k]$, we have

$\displaystyle\operatorname{Tran}(\boldsymbol{\vartheta},\widetilde{\boldsymbol{\vartheta}})$
$\displaystyle\leq\sum_{i=1}^{k}\sum_{j=1}^{k}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}=\sum_{p}\sum_{i\in
C_{p}}\sum_{j=1}^{k}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}.$

Next, we bound the RHS. If $\sum_{i\in C_{p}}\sum_{j=1}^{k}\mu_{i,j}\leq
2k\sqrt{\epsilon}$ for some specific $p$, we directly get

$\sum_{i\in
C_{p}}\sum_{j=1}^{k}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}\leq
2k\sqrt{\epsilon}$

since $\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}\leq
1$ always holds. Otherwise we have $\sum_{i\in C_{p}}w_{i}=\sum_{i\in
C_{p}}\sum_{j=1}^{k}\mu_{i,j}>2k\sqrt{\epsilon}$. In this case, we decompose
the last summation by

$\displaystyle\sum_{i\in
C_{p}}\sum_{j=1}^{k}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}$
$\displaystyle=\sum_{i\in
C_{p}}\sum_{j\notin\widetilde{C}_{p}}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}+\sum_{i\in
C_{p}}\sum_{j\in\widetilde{C}_{p}}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}.$

For the first term, we have $\sum_{i\in
C_{p}}\sum_{j\notin\widetilde{C}_{p}}\mu_{i,j}\leq 2\sqrt{\epsilon}$ as we
stated above and
$\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}\leq 1$
according to their definitions. Thus the first term can be bounded by
$2\sqrt{\epsilon}$. For the second term, we have $\sum_{i\in
C_{p}}\sum_{j\in\widetilde{C}_{p}}\mu_{i,j}\leq 1$ according to the property
of coupling distribution and the following lemma shows
$\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}$ is
small. We delay the proof for readability.

######  Lemma 4.4.

If $\sum_{i\in C_{p}}w_{i}>2k\sqrt{\epsilon}$, then
$\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}\leq
O(k^{3}d^{2}\sqrt{\epsilon})$ for all $i\in C_{p}$ and
$j\in\widetilde{C}_{p}$.

As a result,

$\displaystyle\sum_{i\in
C_{p}}\sum_{j=1}^{k}\mu_{i,j}\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}\leq
O(k^{3}d^{2}\sqrt{\epsilon}).$

Hence, we can conclude that
$\operatorname{Tran}(\boldsymbol{\vartheta},\widetilde{\boldsymbol{\vartheta}})\leq
O(k^{3}d^{2}\sqrt{\epsilon}).$ ∎

###### Proof of Lemma 4.4.

We prove the statement by showing $\widetilde{\boldsymbol{\alpha}}_{j}$ and
$\boldsymbol{\alpha}_{i}$ are close for each dimension $t$.

Firstly, for each $t$, we claim that there must exist some spike
$(\widehat{y}_{s,t},\widehat{\alpha}_{s,t})$ of weight $\widehat{w}_{s,t}$ in
$\widehat{\boldsymbol{\phi}}_{t}$ such that
$\widehat{w}_{s,t}\geq\sqrt{\epsilon}$ and
$|\widehat{y}_{s,t}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime}}|+|\widehat{\alpha}_{s,t}-\alpha_{i^{\prime},t}|\leq\sqrt{\epsilon}$
for some $i^{\prime}\in C_{p}$. If it is not the case, that is, each spike in
$\widehat{\boldsymbol{\phi}}_{t}$ is either of weight less than
$\sqrt{\epsilon}$ or it is of distance $\sqrt{\epsilon}$ from spikes in
$C_{p}$, then, each unit of the weights in $C_{p}$ suffer a transportation
cost of at least $\sqrt{\epsilon}$ after first $k\times\sqrt{\epsilon}$
weights (Since spikes with low weight is not sufficient for covering $C_{p}$
while other spikes are far away). Recall the total weight in $C_{p}$ is at
least $\sqrt{\epsilon}$. So in this case, the transportation distance between
$\widehat{\boldsymbol{\phi}}$ and
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$ is at least
$\operatorname{Tran}(\widehat{\boldsymbol{\phi}}_{t},\Pi_{[\boldsymbol{r},e_{t}]}(\boldsymbol{\vartheta}))\geq(2k\sqrt{\epsilon}-\sqrt{\epsilon}\times
k)\times\sqrt{\epsilon}>k\epsilon$ which contradicts to our assignment for
$\epsilon$. Since we have
$|\widetilde{y}_{j}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime}}|\leq\frac{\sqrt{\epsilon}}{2}$
for some $i^{\prime\prime}\in C_{p}$ according to our construction of
$\widetilde{C}_{p}$ and
$|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime}}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime}}|\leq
k\sqrt{\epsilon}$ according to the clustering, with triangular inequality, we
have

$|\widetilde{y}_{j}-\widehat{y}_{s,t}|\leq|\widehat{y}_{s,t}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime}}|+|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime}}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime}}|+|\widetilde{y}_{j}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime}}|\leq(1.5+k)\sqrt{\epsilon}.$

Since $s$ satisfies the condition in Line 9 of Algorithm 2, the minimization
ensures
$|\widetilde{y}_{j}-\widehat{y}_{s^{*}_{j,t},t}|\leq(1.5+k)\sqrt{\epsilon}.$

Moreover, we also claim that there must exist some spike
$(\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime}},\alpha_{i^{\prime\prime\prime},t})$
in $\Pi_{[\boldsymbol{r},e_{t}]}(\boldsymbol{\vartheta})$ such that
$|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime}}-\widehat{y}_{s^{*}_{j,t},t}|+|\alpha_{i^{\prime\prime\prime},t}-\widehat{\alpha}_{s^{*}_{j,t},t}|\leq
2\sqrt{\epsilon}$. If it is not the case, then the spike
$(\widehat{y}_{s^{*}_{j,t},t},\widehat{\alpha}_{s^{*}_{j,t},t})$, which is of
weight at least $\sqrt{\epsilon}$, is distant from spikes in
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$. This makes the transportation
between $\widehat{\boldsymbol{\phi}}$ and
$\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$ exceed
$\operatorname{Tran}(\widehat{\boldsymbol{\phi}}_{t},\Pi_{[\boldsymbol{r},e_{t}]}(\boldsymbol{\vartheta}))\geq\sqrt{\epsilon}\times
2\sqrt{\epsilon}>2\epsilon$, which contradicts to our assignment for
$\epsilon$. Since we have
$|\widetilde{y}_{j}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime\prime}}|\leq\frac{\sqrt{\epsilon}}{2}$
for some $i^{\prime\prime\prime\prime}\in C_{p}$ according to our construction
of $\widetilde{C}_{p}$ and
$|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime\prime}}|\leq
k\sqrt{\epsilon}$ for each specific $i\in C_{p}$. Recall we already proved
$|\widetilde{y}_{j}-\widehat{y}_{s^{*}_{j,t},t}|\leq(1.5+k)\sqrt{\epsilon}.$
By triangle inequality, we have

$\displaystyle\quad\
|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime}}|+|\alpha_{i^{\prime\prime\prime},t}-\widehat{\alpha}_{s^{*}_{j,t},t}|$
$\displaystyle\leq|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime\prime}}|+|\widetilde{y}_{j}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime\prime}}|+|\widetilde{y}_{j}-\widehat{y}_{s^{*}_{j,t},t}|+|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime}}-\widehat{y}_{s^{*}_{j,t},t}|+|\alpha_{i^{\prime\prime\prime},t}-\widehat{\alpha}_{s^{*}_{j,t},t}|$
$\displaystyle\leq(4+2k)\sqrt{\epsilon}.$

According to Lemma 4.2, we have
$|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime}}|\geq\frac{\eta}{k^{2}d}\|\boldsymbol{\alpha}_{i}-\boldsymbol{\alpha}_{i^{\prime\prime\prime}}\|_{1}\geq\frac{\eta}{k^{2}d}|\alpha_{i,t}-\alpha_{i^{\prime\prime\prime},t}|$.
Thus,

$\displaystyle|\widetilde{\alpha}_{j,t}-\alpha_{i,t}|$
$\displaystyle=|\widehat{\alpha}_{s^{*}_{j,t},t}-\alpha_{i,t}|$
$\displaystyle\leq|\alpha_{i,t}-\alpha_{i^{\prime\prime\prime},t}|+|\alpha_{i^{\prime\prime\prime}}-\widehat{\alpha}_{s^{*}_{j,t},t}|$
$\displaystyle\leq|\alpha_{i,t}-\alpha_{i^{\prime\prime\prime},t}|+\frac{k^{2}d}{\eta}|\alpha_{i^{\prime\prime\prime}}-\widehat{\alpha}_{s^{*}_{j,t},t}|$
$\displaystyle\leq\frac{k^{2}d}{\eta}(|\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i}-\boldsymbol{r}^{\top}\boldsymbol{\alpha}_{i^{\prime\prime\prime}}|+|\alpha_{i^{\prime\prime\prime},t}-\widehat{\alpha}_{s^{*}_{j,t},t}|)$
$\displaystyle\leq\frac{k^{2}d}{\eta}\cdot(4+2k)\sqrt{\epsilon}\leq
O(k^{3}d\sqrt{\epsilon})$

where the equality dues to the assignment of
$\widetilde{\alpha}_{j,t}\leftarrow\widehat{\alpha}_{s_{j,t}^{*},t}$ in Line
10 of Algorithm 2 and the second inequality dues to $\eta<1$. Taking summation
over $t$, we conclude that

$\displaystyle\|\widetilde{\boldsymbol{\alpha}}_{j}-\boldsymbol{\alpha}_{i}\|_{1}\leq\sum_{t=1}^{d}|\widetilde{\alpha}_{j,t}-\alpha_{i,t}|\leq
O(k^{3}d^{2}\sqrt{\epsilon}).$

∎

Note that we can bound $\epsilon$ using Lemma 3.13 and Lemma D.8. As a result,
we can provide performance guarantee of our algorithm for the high dimensional
case.

######  Theorem 4.5.

Let
$\widecheck{\boldsymbol{\vartheta}}=(\widecheck{\boldsymbol{\alpha}},\widetilde{\boldsymbol{w}})$
be the final result (Line 17) in Algorithm 2. Then with probability at least
$1-\eta$ for any fixed constant $0<\eta<1$,

$\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq(kd\xi^{\frac{1}{k}})^{O(1)}.$

###### Proof.

Since $\widetilde{\boldsymbol{\phi}}^{\prime}$ is generated by one dimension
algorithm, according to Lemma 3.13, we have

$\operatorname{Tran}(\widetilde{\boldsymbol{\phi}}^{\prime},\Pi_{\frac{\boldsymbol{r}+\mathbf{1}}{4}}(\boldsymbol{\vartheta}))\leq
O(k\xi^{\frac{1}{4k-2}}).$

Note that $\Pi_{\frac{\boldsymbol{r}+\mathbf{1}}{4}}(\boldsymbol{\vartheta})$
can be transformed to $\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$ by
transforming the coordinates by $\alpha\rightarrow 4\alpha-1$. We have

$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))\leq
4\operatorname{Tran}(\widetilde{\boldsymbol{\phi}}^{\prime},\Pi_{\frac{\boldsymbol{r}+\mathbf{1}}{4}}(\boldsymbol{\vartheta}))\leq
O(k\xi^{\frac{1}{4k-2}}).$

Similarly, according to Lemma D.8, we can see
$\operatorname{Tran}(\widehat{\boldsymbol{\phi}}_{t},\Pi_{[\boldsymbol{r},e_{t}]}(\boldsymbol{\vartheta}))\leq
O(k\xi^{\frac{1}{4k-2}}).$ Thus, from the definition of $\epsilon$, we have
$\epsilon\leq O(k\xi^{\frac{1}{4k-2}})$. With Lemma 4.3, we have

$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq
O(k^{3}d^{2}\sqrt{\epsilon})\leq(kd\xi^{\frac{1}{k}})^{O(1)}.$

Note that $\Delta_{d-1}$ is a convex set and
$\widecheck{\boldsymbol{\alpha}}_{j}$ are the projection of
$\widetilde{\boldsymbol{\alpha}}_{j}$ to the same set. Thus,
$\|\boldsymbol{\alpha}_{i}-\widecheck{\boldsymbol{\alpha}}_{j}\|_{1}\leq\|\boldsymbol{\alpha}_{i}-\widetilde{\boldsymbol{\alpha}}_{j}\|_{1}$
holds for all $A_{i}\in\Delta_{d-1}$ and $j$. Hence, the projection to
$\Delta_{d-1}$ does not increase the transportation distance. As a result,
$\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\widecheck{\boldsymbol{\vartheta}})\leq\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta}).$
Finally, by triangle inequality, we have

$\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq\operatorname{Tran}(\widetilde{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})+\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\widetilde{\boldsymbol{\vartheta}})\leq(kd\xi^{\frac{1}{k}})^{O(1)}.$

∎

By choosing $\xi=(\epsilon/(dk))^{O(k)}$, the previous theorem directly
implies Theorem 1.2.

##  5 Applications to Topic Models

In this section, we discuss the application of our results to topic models.
Here, the underlying vocabulary consists of $d$ words. We are given a corpus
of documents. We adopt the popular “bag of words” model and take each document
as an unordered multiset of words. The assumption is that there is a small
number $k$ of “pure” topics, where each topic is a distribution over $[d]$. A
$K$-word document (i.e., a string in $[d]^{K}$) is generated by first
selecting a topic $i\in\Delta_{d-1}$ from the mixture
$\boldsymbol{\vartheta}$, and then sampling $K$ words i.i.d. according to
$\boldsymbol{\alpha}_{i}$ from this topic. Here, $K$ is referred to as the
snapshot number and we call such a sample a $K$-snapshot of $p$. Our goal here
is to recover
$\boldsymbol{\vartheta}\in\mathrm{Spike}(\Delta_{d-1},\Delta_{k-1})$, which is
a discrete distribution over $k$ pure topics. Again, we measure the accuracy
in terms of $L_{1}$-transportation distance.

###  5.1 An Algorithm for $d=O(k)$

In this section, we first directly handle the high dimensional case using
Algorithm 2. We will perform dimension reduction when $d\gg k$, in Section
5.2.

######  Lemma 5.1.

There is an algorithm that can learn an arbitrary $k$-spike mixture of
discrete distributions supported in $\Delta_{d-1}$ for any $d$ within $L_{1}$
transportation distance $\epsilon$ with probability at least $0.99$ using
$(kd/\epsilon)^{O(k)}$ many $(2k-1)$-snapshots.

The following lemma is a high dimensional version for empirical error (see
e.g. [55], [45], [31]).

######  Lemma 5.2.

For matrix
$R=[\boldsymbol{r}_{1},\cdots,\boldsymbol{r}_{p}]\in\mathbb{R}^{d\times p}$
with $\|R\|_{\infty}\leq 1$, using $2^{O(kp)}\cdot\xi^{-2}\cdot\log(1/\eta)$
samples, we can obtain transformed empirical moments $M^{\prime}(R)$
satisfying

$\max_{|\mathbf{t}|\leq
K}|M^{\prime}_{\mathbf{t}}(R)-M_{\mathbf{t}}(\Pi_{R}(\boldsymbol{\vartheta}))|\leq\xi$

with probability at least $1-\eta$.

###### Proof.

Assume the words of a document are $x_{1},\cdots,x_{K}$. If the document is
from some topic $i$, we have $\mathbb{E}[x_{j}=c]=\alpha_{p,t}$ for all
$i,j,c$. For each word $x_{j}$, we generate a i.i.d. random vector
$\boldsymbol{\beta}_{j}=[\beta_{1,j},\cdots,\beta_{p,j}]^{\top}\in\\{0,1\\}^{p}$
such that $\mathbb{E}[\beta_{c,j}]=r_{c,x_{j}}$. For
$\mathbf{t}=(t_{1},\cdots,t_{p})$ with $|\mathbf{t}|\leq K$, we define
observable $h_{\mathbf{t}}=\mathbf{1}[\forall
c:\sum_{j=1}^{K}\beta_{c,j}=t_{c}]$.

Note that $h_{\mathbf{t}}$ is the normalized histogram for projected
distribution $\Pi_{R}(\boldsymbol{\vartheta})$. Hence, from the relation
between normalized histogram and standard moments, we can conclude that

$\displaystyle
M(\Pi_{R}(\boldsymbol{\vartheta}))=\left(\bigotimes_{t=1}^{p}\textrm{Pas}\right)h$

where $\textrm{Pas}\in\mathbb{R}^{(K+1)\times(K+1)}$ such that
$\textrm{Pas}_{i,j}=[j\geq i]\binom{j}{i}\binom{K+1}{i}$ and $\otimes$ is the
tensor product. Since $\|\textrm{Pas}\|\leq 2^{O(k)}$, we have
$\left(\bigotimes_{t=1}^{p}\textrm{Pas}\right)\leq 2^{O{(kp)}}$.

Let $\widehat{h}_{\mathbf{t}}$ be the empirical average of $h_{\mathbf{t}}$.
Since $h_{\mathbf{t}}$ is a Bernoulli variable, from Hoeffding’s inequality,
we have
$\Pr[|\widehat{h}_{\mathbf{t}}-\mathbb{E}[h_{\mathbf{t}}]|<\epsilon]\leq
2\exp(-2\epsilon^{2}s)$. By applying union bound on $\mathbf{t}$, we can see
that using $s>2^{O(kp)}\cdot\xi^{-2}\cdot\log(1/\eta)$ samples, we have
$|\widehat{h}_{\mathbf{t}}-\mathbb{E}[h_{\mathbf{t}}]|<\xi\cdot
2^{-\Omega(kp)}$ with probability at least $1-\eta$. Moreover, if we calculate
$M^{\prime}$ from $\widehat{h}$ using the relationship between $M$ and $h$,
this leads to $\max_{|\mathbf{t}|\leq
K}|M^{\prime}_{\mathbf{t}}(R)-M_{\mathbf{t}}(\Pi_{R}(\boldsymbol{\vartheta}))|\leq\xi$.
∎

###### Proof of Lemma 5.1.

Let $\xi=(\epsilon/kd)^{O(k)}$. From Theorem 4.5, the recovered $k$-spike
distribution $\widecheck{\boldsymbol{\vartheta}}$ satisfies
$\operatorname{Tran}(\widecheck{\boldsymbol{\vartheta}},\boldsymbol{\vartheta})\leq(kd\xi^{\frac{1}{k}})^{O(1)}=O(\epsilon)$.
According to Lemma 5.2, we can construct such a noisy moment oracle
$M^{\prime}(\cdot)$ using $(kd/\epsilon)^{O(k)}$ samples with high
probability. ∎

###  5.2 Dimension Reduction When $d\gg k^{\Omega(1)}$

In topic modeling, the number of words is typically much larger than the
number of pure topics (i.e., $d\gg k$). In this case, we face another
challenge to recover the mixture. First, there are
$\binom{d}{1}+\binom{d}{2}+\cdots+\binom{d}{k}=O(d^{k})$ many different
moments. To obtain all empirical moments accurate enough would require a huge
number of $2k-1$-snapshot samples. So if $d\gg k$, we reduce the dimension
from $d$ to $O(k)$. Now, we prove the following theorem. It improves on the
result in [55, 45], which uses more than $(k/\epsilon)^{O(k^{2})}$
$(2k-1)$-snapshots, and [31] which use
$(k/\epsilon)^{O(k)}\cdot(w_{\min}\zeta^{k})^{-O(1)}$ $2k$-snapshots and
requires the minimum separation assumption (See Table 2).

######  Theorem 5.3.

There is an algorithm that can learn an arbitrary $k$-spike mixture supported
in $\Delta_{d-1}$ for any $d$ within $L_{1}$ transportation distance
$O(\epsilon)$ with probability at least $0.99$ using
$\operatorname{poly}(d,k,\frac{1}{\epsilon})$,
$\operatorname{poly}(d,k,\frac{1}{\epsilon})$, many $1$-,$2$-,snapshots and
$(k/\epsilon)^{O(k)}$ many $(2k-1)$-snapshots.

Let $\boldsymbol{\vartheta}$ be the target $k$-spike mixture in $\Delta_{d-1}$
($d\gg k$). Suppose there is a learning algorithm ${\mathcal{A}}$ satisfying
the assumption in the theorem. We show how to apply ${\mathcal{A}}$ to the
projection of $\boldsymbol{\vartheta}$ to a subspace of dimension at most $k$
in $\mathbb{R}^{d}$. However, an arbitrary subspace of dimension at most $k$
is not enough, even the one spanning the $k$ spikes of
$\boldsymbol{\vartheta}$. This is because $L_{1}$ distance is not rotationally
invariant and not preserved under projection. For example, the $L_{1}$
distance between $(1/d,\ldots,1/d)$ and $(2/d,\ldots,2/d,0,\ldots,0)$ is 1 in
$\mathbb{R}^{d}$ but only $\sqrt{1/d}$ in the line spanned by the two points.
Hence, an accurate estimation of the projected measure
$\boldsymbol{\vartheta}_{B}=\Pi_{B}(\boldsymbol{\vartheta})$ may not translate
to an accurate estimation in the ambient space $\mathbb{R}^{d}$. Here, we can
use the dimension reduction method developed in [45], which shows that if
suffices to project the mixture to a special subspace $B$, such that a unit
$L_{1}$-ball in $B$ ($L_{1}$ measured in $\mathbb{R}^{d}$) is close to being
an $\widetilde{O}(d^{-1/2})$ $L_{2}$-ball ($\widetilde{O}$ hides factors
depending only on $k$ and $\epsilon$, but not $d$).

######  Lemma 5.4.

([45]) There is an algorithm that requires
$\operatorname{poly}(d,k,\frac{1}{\epsilon})$ many 1-,2- snapshots to
construct a subspace $B$ ($\dim(B)\leq k$) with the following useful
properties:

  1. 1.

Suppose $\\{b_{1},b_{2},\cdots,b_{m}\\}$ is an orthonormal basis of $B$ where
$m=\dim(B)\leq k$ (such a basis can be produced by the above algorithm as
well). For any $i\in[m]$, $\|b_{i}\|_{\infty}\leq
O(k^{3/2}\epsilon^{-2}d^{-1/2})$.

  2. 2.

Suppose we can learn an approximation
$\widetilde{\boldsymbol{\vartheta}}_{B}$, supported on
$\operatorname{Span}(B)$, of the projected measure
$\boldsymbol{\vartheta}_{B}=\Pi_{\operatorname{Span}(B)}(\boldsymbol{\vartheta})$
such that
$\operatorname{Tran}(\boldsymbol{\vartheta}_{B},\widetilde{\boldsymbol{\vartheta}}_{B})\leq\epsilon_{1}=\operatorname{poly}(1/k,\epsilon)$
(here $L_{1}$ is measured in $\mathbb{R}^{d}$, not in the subspace) using
$N_{1}(d)$, $N_{2}(d)$ and $N_{K}(d)$ 1-, 2-, and $K$-snapshot samples. Then
there is an algorithm for learning a mixture
$\widetilde{\boldsymbol{\vartheta}}$ such that
$\operatorname{Tran}(\boldsymbol{\vartheta},\widetilde{\boldsymbol{\vartheta}})\leq\epsilon$
using $O(N_{1}(d/\epsilon)+d\log d/\epsilon^{3})$,
$O(N_{2}(d/\epsilon)+O(k^{4}d^{3}\log n/\epsilon^{6}))$ and
$O(N_{K}(d/\epsilon))$ 1-, 2-, and $K$-snapshot samples respectively.

###### Proof of Theorem 5.3.

In light of Lemma 5.4, we only need to show how to learn the projection
$\boldsymbol{\vartheta}_{B}=\Pi_{\operatorname{Span}(B)}(\boldsymbol{\vartheta})$
using the algorithm developed in Section 5.1.

First, we show how to translate a $K$-snapshots from $\boldsymbol{\vartheta}$
to a $K$-snapshot in some new $k$-dimensional mixture
$\boldsymbol{\vartheta}_{\text{new}}$ (then we apply the algorithm in Section
5.1 for dimension $k$). Let $L=\sum_{i=1}^{m}\|b_{i}\|_{\infty}$. So $L\leq
O(k^{5/2}\epsilon^{-2}d^{-1/2})$. Let $f:[-L,L]\rightarrow[0,1]$ defined as
$f(x)=\frac{x}{2mL}+\frac{1}{2m}$. For each word in the original $K$-snapshot,
say $i\in[d]$, we translate it to $j$ with probability $q_{i,j}=f(b_{j,i})$,
for each $j\in[m]$, where $b_{j,i}$ denotes the $i$th coordinate of $b_{j}$,
and translate it into $m+1$ with probability
$q_{i,m+1}=1-\sum_{j=1}^{m}q_{i,j}$. Since we have
$\sum_{j=1}^{m}f(b_{j,i})=\frac{1}{2mL}\sum_{j=1}^{m}b_{j,i}+\frac{1}{2}$, we
know that $q_{i,1},q_{i,2},\cdots,q_{i,m+1}\in[0,1]$ and
$\sum_{j=1}^{k+1}q_{i,j}=1$. So this is a well defined $k$-spike mixture
$\boldsymbol{\vartheta}_{\text{new}}$ and its $K$ snapshots. Then we apply our
learning algorithm ${\mathcal{A}}$ for dimension $k$ to obtain some
$\boldsymbol{\vartheta}_{\text{new}}^{\prime}$ such that
$\operatorname{Tran}(\boldsymbol{\vartheta}_{\text{new}},\boldsymbol{\vartheta}_{\text{new}}^{\prime})\leq\epsilon_{1}=\operatorname{poly}(1/k,\epsilon)$
with $\operatorname{poly}(k,1/\epsilon_{1})=\operatorname{poly}(k,1/\epsilon)$
many $2k-1$-snapshots.

Next we can see $\boldsymbol{\vartheta}_{\text{new}}$ (in $\mathbb{R}^{k}$)
and $\boldsymbol{\vartheta}_{B}$ (in $\operatorname{Span}(B)$) are related by
an affine transform. For a spike
$\alpha_{i}\in\mathrm{Supp}(\boldsymbol{\vartheta})$, one can see that
$\Pi_{B}(\alpha_{i})=\sum_{j=1}^{m}a_{i,j}b_{j}$ , where
$a_{i,j}=\langle\alpha_{i},b_{j}\rangle$. $\alpha_{i}$ produces a new spike
$\beta_{i}$ in $\boldsymbol{\vartheta}_{\text{new}}$ and the mapping is as
follows: for each $j\in[k]$, we have

$\beta_{i,j}=\sum_{t=1}^{d}\alpha_{i,t}f(b_{j,t})=\sum_{t=1}^{d}\alpha_{i,t}(\frac{b_{j,t}}{2mL}+\frac{1}{2m})=\frac{a_{i,j}}{2mL}+\frac{1}{2m}$

($\beta_{i,j}$ denotes the $j$th coordinate of $\beta_{i}$). So
$a_{i,j}=2mL(\beta_{i,j}-\frac{1}{2m})$ and $g(x)=2mLx-L$ is the affine
transformation. Now we can translate
$\boldsymbol{\vartheta}_{\text{new}}^{\prime}$ into a mixture in $Span(B)$ by
applying $g(.)$ in each coordinate of each
$\beta_{i}\in\mathrm{Supp}(\boldsymbol{\vartheta}_{\text{new}}^{\prime})$, and
obtain an estimation of $\boldsymbol{\vartheta}_{B}$, say
$\boldsymbol{\vartheta}_{B}^{\prime}$.

It remains to show
$\operatorname{Tran}(\boldsymbol{\vartheta}_{B}^{\prime},\boldsymbol{\vartheta}_{B})\leq\epsilon_{1}$.
If we let $\operatorname{Tran}_{B}$ denote the transportation distance in
$\operatorname{Span}(B)$ (where if
$a=\sum_{j=1}^{k}a_{j}b_{j},c=\sum_{j=1}^{k}c_{j}b_{j}$ then
$\operatorname{Tran}_{B}(a,c)=\sum_{j=1}^{k}|a_{j}-c_{j}|$). Then
$\operatorname{Tran}_{B}(\boldsymbol{\vartheta}_{B}^{\prime},\boldsymbol{\vartheta}_{B})\leq
2mL\operatorname{Tran}(\boldsymbol{\vartheta}_{\text{new}}^{\prime},\boldsymbol{\vartheta}_{\text{new}})=\operatorname{poly}(\epsilon,\frac{1}{k})/\sqrt{d}$.
Hence, we have

$\operatorname{Tran}(\boldsymbol{\vartheta}_{B}^{\prime},\boldsymbol{\vartheta}_{B})\leq
d\max_{i}\|b_{i}\|_{\infty}\operatorname{Tran}_{B}(\boldsymbol{\vartheta}_{B}^{\prime},\boldsymbol{\vartheta}_{B})\leq\epsilon_{1}.$

where the first inequality holds since
$\operatorname{Tran}(a,c)=\|\sum_{j=1}^{k}(a_{j}-c_{j})b_{j}\|_{1}\leq
d\max_{i}\|b_{i}\|_{\infty}\sum_{j=1}^{k}(a_{j}-c_{j})$ for any two point
$a,c\in\operatorname{Span}(B)$. ∎

##  6 Applications to Gaussian Mixture Learning

In this section, we show how to leverage our results for sparse moment problem
to obtain improved algorithms for learning Gaussian mixtures. We consider the
following setting studied in [61, 26]. A $k$-Gaussian mixture in
$\mathbb{R}^{d}$ can be parameterised as
$\boldsymbol{\vartheta}_{N}=(\boldsymbol{\alpha},\boldsymbol{w},\Sigma)$.
Here,
$\boldsymbol{\alpha}=\\{\boldsymbol{\alpha}_{1},\boldsymbol{\alpha}_{2},\cdots,\boldsymbol{\alpha}_{k}\\}$
and $\boldsymbol{w}=\\{w_{1},w_{2},\cdots,w_{k}\\}\in\Delta_{k-1}$ where
$\boldsymbol{\alpha}_{i}\in\mathbb{R}^{d}$ specifies the mean of $i$th
component and $w_{i}\in[0,1]$ is the corresponding weight.
$\Sigma\in\mathbb{R}^{d\times d}$ is the common covariance matrix for all $k$
mixture components and we assume that $\Sigma$ is known in advance. We further
assume $\|\boldsymbol{\alpha}_{i}\|_{2}\leq 1$ and the maximal eigenvalue
$\|\Sigma\|_{2}$ is bounded by a constant. For $k$-Gaussian mixture
$\boldsymbol{\vartheta}_{N}=(\boldsymbol{\alpha},\boldsymbol{w},\Sigma)$, each
observation is distributed as

$\displaystyle\boldsymbol{\vartheta}_{N}\sim\sum_{i=1}^{k}w_{i}N(\boldsymbol{\alpha}_{i},\Sigma).$

We consider the parameter learning problem, that is, to learn the parameter
$\boldsymbol{\alpha}$ and $\boldsymbol{w}$ given known covariance matrix
$\Sigma$ and a set of i.i.d. samples from $\boldsymbol{\vartheta}_{N}$. The
model is also called Gaussian location mixture model [61, 26]. 555 Wu and Yang
[61] also studied the problem with unknown $\Sigma$. We leave it as a future
direction.

###  6.1 Efficient Algorithm for $d=1$

In the 1-dimensional case, we denote the known variance by $\sigma$ which is
upper bounded by some constant. As shown in [61], the moments of Gaussian
mixture has a close connection with the moments of corresponding discrete
distributions $\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$.
For $x\sim N(\mu,1)$, we have ${\mathbb{E}}[H_{t}(x)]=\mu^{t}$ for Hermite
polynomial $H_{t}(x)$ which is defined as

$\displaystyle
H_{t}(x)=\sum_{i=0}^{t}h_{t,i}x^{i}=t!\sum_{j=0}^{\lfloor\frac{t}{2}\rfloor}\frac{(-1/2)^{j}}{j!(t-2j)!}x^{t-2j}.$
(5)

For 1-dimensional Gaussian mixture
$\boldsymbol{\vartheta}_{N}=(\boldsymbol{\alpha},\boldsymbol{w},\sigma)$ with
variance $\sigma$, the $t$th moment of the discrete distribution
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$ satisfies

$\displaystyle
M_{t}(\boldsymbol{\vartheta})={\mathbb{E}}_{x\sim\boldsymbol{\vartheta}_{N}}\left[\sum_{i=0}^{t}h_{t,i}\sigma^{t-i}x^{i}\right].$

The following lemma guarantees the performance of estimating the moment by
sampling.

######  Lemma 6.1.

(Lemma 5 of [61], restated) Let
$x_{1},\cdots,x_{n}\sim\boldsymbol{\vartheta}_{N}$ be $n$ independent samples.

$\widetilde{M}_{t}(\boldsymbol{\vartheta})=\frac{1}{n}\sum_{j=1}^{n}\left[\sum_{i=0}^{t}h_{t,i}\sigma^{t-i}x^{i}\right]$

is an unbiased estimator for $M_{t}(\boldsymbol{\vartheta})$, and the variance
of this estimator can be bounded by

$\operatorname{Var}[{\widetilde{M}_{t}(\boldsymbol{\vartheta})}]\leq\frac{1}{n}(\sigma
t)^{O(t)}.$

Thus we can compute the moment of original distribution and use our Algorithm
1 for recovering the parameter of $\boldsymbol{\vartheta}$. More concretely,
we can replace the last two lines (an SDP) in Algorithm 2 of [61] by our
Algorithm 1, and obtain the following theorem. The post-sampling time is
improved from from $O(k^{6.5})$ in [61] to $O(k^{2})$.

######  Theorem 6.2.

Let $\boldsymbol{\vartheta}_{N}$ be an arbitrary $k$-Gaussian mixture over
$\mathbb{R}$ with means $\alpha_{1},\cdots,\alpha_{n}$ and known variance
$\sigma$ bounded by some constant. There is an algorithm that can learn the
parameter $\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$ within
transportation distance $O(\epsilon)$, with probability at least $0.99$, using
$(k/\epsilon)^{O(k)}$ samples from the mixture $\boldsymbol{\vartheta}_{N}$.
Moreover, once we obtain the estimation of the moments of
$\boldsymbol{\vartheta}_{N}$, our algorithm only uses $O(k^{2})$ arithmetic
operations.

###### Proof.

Let $c$ be the constant to be determined. With $n=(k/\epsilon)^{ck}$ samples,
the variance of empirical moment can be bounded by

$\displaystyle\operatorname{Var}[M_{t}^{\prime}(\boldsymbol{\vartheta})]\leq\frac{1}{(k/\epsilon)^{ck}}\cdot
k^{O(k)}\leq(\epsilon/k)^{\Omega(k)}$

where the first inequality holds due to Lemma 6.1 and $\sigma$ is bounded by a
constant, and the second inequality holds by selecting a large enough constant
$c>0$. According to Chebyshev’s inequality, for each $t$, with probability at
least $1-0.01k^{-2}$,

$\displaystyle|M^{\prime}_{t}(\boldsymbol{\vartheta})-M_{t}(\boldsymbol{\vartheta})|\leq
0.1k\cdot(\epsilon/k)^{\Omega(k)}\leq(\epsilon/k)^{\Omega(k)}.$

By taking union bound, the probability that the above inequality holds for all
$0\leq t\leq 2k-1$ is greater than $0.99$. We can conclude the result by
applying Theorem 1.1 directly. ∎

###  6.2 Efficient Algorithm for $d>1$

For higher dimensional Gaussian mixture, we can reduce the problem to learning
the discrete mixture
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$, and leverage
Algorithm 2 to solve the problem. Assume the locations of all Gaussians are in
the unit ball. The error can be bounded easily according to Theorem 1.2.

Firstly, as shown in [26], one can use SVD to reduce a $d$-dimensions problem
to a $k$-dimension problem using $\operatorname{poly}(1/\epsilon,d,k)$
samples. Thus, we only need to consider settings with $d\leq k$. Similar to
the 1-d case, we can transform the problem of learning Gaussian mixture
$\boldsymbol{\vartheta}_{N}=(\boldsymbol{\alpha},\boldsymbol{w},\Sigma)$ to
the problem of learning the discrete mixture
$\boldsymbol{\vartheta}=(\boldsymbol{\alpha},\boldsymbol{w})$. In particular,
as we show below, we can estimate the moments of the projection of
$\boldsymbol{\vartheta}$ from the samples. After obtaining the noisy moment
information, we can apply Algorithm 2 to recover the parameter
$(\boldsymbol{\alpha},\boldsymbol{w})$. We mention that we need to modify
Algorithm 2 slightly: we change to domain from $\Delta_{d-1}$ to the unit
ball. We note this does not affect any parts of the proof since Theorem 1.2
only requires the projected space to be a convex domain.

The remaining task is to show how to estimate the projected moments for 1-d
and 2-d projections. To estimate the moments of 1-dimensional projection, we
can use the estimator in the last section. In particular, for the projected
measure $\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta})$ where $\boldsymbol{r}$
is an arbitrary unit vector, the $t$-th moment can be computed by

$\displaystyle
M_{t}(\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))={\mathbb{E}}_{\boldsymbol{x}\sim\boldsymbol{\vartheta}_{N}}\left[\sum_{i=0}^{t}h_{t,i}(\boldsymbol{r}^{\top}\Sigma\boldsymbol{r})^{t-i}(\boldsymbol{r}^{\top}\boldsymbol{x})^{i}\right].$

We can use sample average to estimate
$M_{t}(\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))$, and we denote the
estimation as
$\widetilde{M}_{t}(\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))$.

Now, we consider the problem of estimating the moments for 2-dimensional
projection along $\boldsymbol{r}_{1}$ and $\boldsymbol{r}_{2}$. First, we
compute the vector
$\boldsymbol{r}_{2}^{\prime}=\boldsymbol{r}_{2}-\frac{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}}{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}}\boldsymbol{r}_{1}$.
According to Gram–Schmidt process, $\boldsymbol{r}_{2}^{\prime}$ is
$\Sigma$-orthogonal to $\boldsymbol{r}_{1}$, i.e.,
$\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}^{\prime}=0$. Moreover,
$\boldsymbol{r}_{1}^{\top}\boldsymbol{x}$ and
$\boldsymbol{r}_{2}^{\prime\top}\boldsymbol{x}$ are independent Gaussian
distributions for variable $\boldsymbol{x}\sim N(0,\Sigma)$. Note that for
independent Gaussian variables $x_{1}\sim N(\mu_{1},1)$ and $x_{2}\sim
N(\mu_{2},1)$, we have
${\mathbb{E}}[H_{t_{1}}(x_{1})H_{t_{2}}(x_{2})]={\mathbb{E}}[H_{t_{1}}(x_{1})]{\mathbb{E}}[H_{t_{2}}(x_{2})]=\mu_{1}^{t_{1}}\mu_{2}^{t_{2}}$,
where $H$ is the Hermite polynomial defined in (5). As a result, the projected
moments along $\boldsymbol{r}_{1}$ and $\boldsymbol{r}_{2}^{\prime}$ can be
computed by

$\displaystyle
M_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}^{\prime}]}(\boldsymbol{\vartheta}))={\mathbb{E}}_{\boldsymbol{x}\sim\boldsymbol{\vartheta}_{N}}\left[\sum_{i_{1}=0}^{t_{1}}\sum_{i_{2}=0}^{t_{2}}h_{t_{1},i_{1}}h_{t_{2},i_{2}}(\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1})^{t_{1}-i_{1}}(\boldsymbol{r}_{2}^{\prime\top}\Sigma\boldsymbol{r}_{2}^{\prime})^{t_{2}-i_{2}}(\boldsymbol{r}_{1}^{\top}\boldsymbol{x})^{i_{1}}(\boldsymbol{r}_{2}^{\prime\top}\boldsymbol{x})^{i_{2}}\right].$

Moreover, the projected moments along $\boldsymbol{r}_{1}$ and
$\boldsymbol{r}_{2}$ can be computed by

$\displaystyle
M_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))=\sum_{i=0}^{t_{2}}\binom{t_{2}}{i}\left(\frac{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}}{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}}\right)^{i}M_{(t_{1}+i,t_{2}-i)}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}^{\prime}]}(\boldsymbol{\vartheta}))$

since
$\boldsymbol{r}_{2}^{\top}\boldsymbol{x}=\boldsymbol{r}_{2}^{\prime\top}\boldsymbol{x}+\frac{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}}{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}}\boldsymbol{r}_{1}^{\top}\boldsymbol{x}$
which implies
$(r_{2}^{\top}\boldsymbol{x})^{t_{2}}=\sum_{i=0}^{t_{2}}\binom{t_{2}}{i}\left(\frac{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}}{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}}\right)^{i}(r_{1}^{\top}\boldsymbol{x})^{i}(r_{2}^{\prime\top}\boldsymbol{x})^{t_{2}-i}$.
We use
$\widetilde{M}_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))$
to denote the sample average estimation of
$M_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))$.
We are ready to state the performance guarantee:

######  Theorem 6.3.

The algorithm can learn an arbitrary $d$-dimensional $k$-spike Gaussian
mixture within transportation distance $O(\epsilon)$ given the covariance
matrix $\Sigma$ using
$(k/\epsilon)^{O(k)}+\operatorname{poly}(1/\epsilon,d,k)$ samples.

###### Proof.

We first consider the problems of dimension $d\leq k$. The only thing we need
to show is that the above estimators
$\widetilde{M}_{t}(\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))$ and
$\widetilde{M}_{t}(\Pi_{\boldsymbol{r}_{1},\boldsymbol{r}_{2}}(\boldsymbol{\vartheta}))$
have sufficient accuracy with high probability using
$n=(k/\epsilon)^{\Theta(k)}$ samples. Clearly, both estimators are unbiased.
So it is sufficient to bound the variance of the estimators.

For 1-dimensional projected moments, we can use the same argument as Theorem
6.2. In particular, we can show that with probability at least $0.999$, for
all $0\leq t\leq 2k-1$,

$\displaystyle|\widetilde{M}_{t}(\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))-M_{t}(\Pi_{\boldsymbol{r}}(\boldsymbol{\vartheta}))|\leq(k/\epsilon)^{O(k)}.$

For 2-dimensional projected moments, we first bound the variance of estimator
along $r_{1}$ and $r_{2}^{\prime}$. With the same argument as Lemma 6.1, we
have

$\displaystyle\operatorname{Var}[\widetilde{M}_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}^{\prime}]}(\boldsymbol{\vartheta}))]\leq(\epsilon/k)^{\Omega(k)}.$

for all $t_{1}+t_{2}\leq 2k$ with probability at lest $0.999$. Moreover, since
we randomly chooses the same $r_{1}$ over a sphere of radius $\Theta(1)$ in
Line 6 of Algorithm 2, we have
$\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}\geq\Omega(\frac{\mu}{k^{2}}\|\boldsymbol{r}_{1}\|_{2}^{2}\|\Sigma\|_{2})$
with probability at least $1-\mu$. In this case, we have

$\frac{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}}{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}}\leq\frac{\|\boldsymbol{r}_{1}\|_{2}\|\Sigma\|_{2}\|\boldsymbol{r}_{2}\|}{\Omega(\frac{\mu}{k^{2}}\|\boldsymbol{r}_{1}\|_{2}^{2}\|\Sigma\|_{2})}\leq
O(k^{2}/\mu)$

where the last inequality holds because $\boldsymbol{r}_{1}$ and
$\boldsymbol{r}_{2}$ are of norm $\Theta(1)$. Since
$\sqrt{\operatorname{Var}[x+y]}\leq\sqrt{\operatorname{Var}[x]}+\sqrt{\operatorname{Var}[y]}$
for any random variables $x$ and $y$,

$\displaystyle\sqrt{\operatorname{Var}[\widetilde{M}_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))]}\leq\sum_{i=0}^{t_{2}}\binom{t_{2}}{i}\left(\frac{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{2}}{\boldsymbol{r}_{1}^{\top}\Sigma\boldsymbol{r}_{1}}\right)^{i}\sqrt{\operatorname{Var}[\widetilde{M}_{(t_{1}+i,t_{2}-i)}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}^{\prime}]}(\boldsymbol{\vartheta}))]}.$

As a result, by choosing $\mu=0.001$, we have

$\displaystyle\operatorname{Var}[\widetilde{M}_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))]\leq(\epsilon/k)^{\Omega(k)}$

holds for all $\boldsymbol{r}_{2}$ with probability at least $0.999$.
Conditioning on that this event holds, according Chebyshev’s inequality and a
union bound, with probability at least $1-0.001d$, for all $0\leq
t_{1},t_{2}\leq 2k-1$,

$\displaystyle|\widetilde{M}_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))]-M_{(t_{1},t_{2})}(\Pi_{[\boldsymbol{r}_{1},\boldsymbol{r}_{2}]}(\boldsymbol{\vartheta}))]|\leq(\epsilon/k)^{\Omega(k)}.$

Recall that the algorithm uses 1-d projection for one $\boldsymbol{r}$, 2-d
projection for one $\boldsymbol{r}_{1}$ and $d$ different $\boldsymbol{r}_{2}$
and the algorithm succeeds with high probability statement for random
$\boldsymbol{r}$. Hence, by union bound over all these events, the projected
moment estimations have sufficient accuracy with probability at least $0.995$.

If the dimension $d>k$, we can apply the dimension reduction in [26], which
requires an extra sample complexity of $\operatorname{poly}(1/\epsilon,d,k)$.
The total sample complexity is
$(k/\epsilon)^{O(k)}+\operatorname{poly}(1/\epsilon,d,k)$. ∎

Finally, we discuss the running time during the sampling phase and post-
sampling phase. The dimension reduction requires
$\operatorname{poly}(1/\epsilon,d,k)$ samples and $O(d^{3})$ time [26]. For
the recovery problem in dimension $k$, each $1$-d or $2$-d moment oracle can
be computed in $O(nk^{3})$ time where $n=(k/\epsilon)^{O(k)}$ is the number of
samples. Since we only requires $O(k)$ moment oracles, the sampling time can
be bounded by $O(nk^{4})$. This improves the
$O(n^{5/4}\operatorname{poly}(k))$ sampling time in [26] (their algorithm
requires 1-d projections to $O(n^{1/4}$ many directions). For post-sampling
running time, our Algorithm 2 runs in time $O(k^{4})$ (since $d\leq k$ by
dimension reduction). This improves the $O(n^{1/2}\operatorname{poly}(k))$
post-sampling time in [26].

##  7 Concluding Remarks

We provide efficient and easy-to-implement algorithms for learning the
$k$-mixture over discrete distributions from noisy moments. We measure the
accuracy in terms of transportation distance, and our analysis is independent
of the minimal separation. The techniques used in our analysis for the
1-dimensional case may be useful in the perturbative analysis of other
problems involving the Hankel matrix or the Vandermonde matrix (such as super-
resolution).

Our problem is a special case of learning mixtures of product distribution
[32], which is further a special case of the general problem of learning
mixtures of Bayesian networks [33], which has important applications in causal
inferences. In fact, the algorithms in [32, 33] used the algorithm for sparse
Hausdorff moment problem as a subroutine. It would be interesting to see if
our techniques can be used to improve the learning algorithms for these
problems.

##  8 Acknowledgement

We would like to thank Yuval Rabani and Leonard Schulman for their insightful
discussions. In particular, Leonard Schulman pointed out that Lemma 3.8 is
related to Schur polynomial. JL also would like to thank Xuan Wu for several
discussions around 2015 and a version of Lemma B.1 was proved back then,
independent from [61].

Appendix

##  Appendix A Transportation Distance for Non-probability Measures

The dual formulation (3) can be generalized to non-probability measures.
Specifically, for two signed measure $P$ and $Q$ defined over $\mathbb{R}^{d}$
with $\int P\mathrm{d}x=\int Q\mathrm{d}x=1$, we define its transportation
distance as (3). The equivalent primal formulation is as follows: Define
measure $(P-Q)^{+}$ and $(Q-P)^{+}$ in which
$(P-Q)^{+}(x)=\max\\{P(x)-Q(x),0\\}$. According to the dual form (3), one is
easy to check
$\operatorname{Tran}(P,Q)=\operatorname{Tran}((P-Q)^{+},(Q-P)^{+})$. Note that
both $(P-Q)^{+}$ and $(Q-P)^{+}$ are non-negative. As a result, we conclude
another definition for generalized transportation distance:

$\displaystyle\operatorname{Tran}(P,Q)=\inf\left\\{\int\|x-y\|_{1}\mathrm{d}\mu(x,y):\mu\in
M((P-Q)^{+},(Q-P)^{+})\right\\}.$ (6)

We also need to define the transportation distance over complex measures.
Denote $\mathbb{B}_{\mathbb{C}}=\\{\alpha\in\mathbb{C}:|\alpha|\leq 1\\}$. We
consider a complex measure $P$ such that $\int_{\mathbb{C}}\textrm{d}P(x)=1$.
Denote the real and imaginary components of a complex measure $P$ as
$P^{\textrm{r}}$ and $P^{\textrm{i}}$ respectively. In concrete,
$P(x)=P^{\textrm{r}}(x)+{\mathfrak{i}}P^{\textrm{i}}(x)$. We generalize
transportation distance to complex weights as

$\displaystyle\operatorname{Tran}(P,Q)=\operatorname{Tran}(P^{\textrm{r}},Q^{\textrm{r}})+\operatorname{Tran}(P^{\textrm{i}},Q^{\textrm{i}})$
(7)

By this definition, we can easily see that

$\displaystyle\operatorname{Tran}(P,Q)$ $\displaystyle=\sup\left\\{\int
f^{\textrm{r}}\mathrm{d}(P^{\textrm{r}}-Q^{\textrm{r}}):f^{\textrm{r}}\in
1\text{-}\mathsf{Lip}\right\\}+\sup\left\\{\int
f^{\textrm{i}}\mathrm{d}(P^{\textrm{i}}-Q^{\textrm{i}}):f^{\textrm{i}}\in
1\text{-}\mathsf{Lip}\right\\}$
$\displaystyle=\sup\left\\{\textrm{real}\left(\int(f^{\textrm{r}}-{\mathfrak{i}}f^{\textrm{i}})\mathrm{d}(P-Q)\right):f^{\textrm{r}},f^{\textrm{i}}\in
1\text{-}\mathsf{Lip}\right\\}$ $\displaystyle\leq 2\sup\left\\{\int
f\mathrm{d}(P-Q):f\in 1\text{-}\mathsf{Lip}_{\mathbb{C}}\right\\}$

where
$1\text{-}\mathsf{Lip}_{\mathbb{C}}:=\\{f:\mathbb{C}\rightarrow\mathbb{C}\mid|f(x)-f(y)|\leq|x-y|\text{