# Source Identification for Mixtures of Product Distributions

Spencer L. Gordon Engineering and Applied Science, California Institute of
Technology<EMAIL_ADDRESS>Bijan Mazaheri Engineering and Applied
Science, California Institute of Technology<EMAIL_ADDRESS>Yuval Rabani
The Rachel and Selim Benin School of Computer Science and Engineering, The
Hebrew University of Jerusalem, Jerusalem 9190416, Israel,
<EMAIL_ADDRESS>Research supported in part by NSFC-ISF grant 2553-17
and by NSF-BSF grant 2018687. Part of this work was done while visiting
Caltech. Leonard J. Schulman Engineering and Applied Science, California
Institute of Technology<EMAIL_ADDRESS>Research supported in part by
NSF grants CCF-1618795, 1909972.

###### Abstract

We give an algorithm for source identification of a mixture of $k$ product
distributions on $n$ bits. This is a fundamental problem in machine learning
with many applications. Our algorithm identifies the source parameters of an
identifiable mixture, given, as input, approximate values of multilinear
moments (derived, for instance, from a sufficiently large sample), using
$2^{O(k^{2})}n^{O(k)}$ arithmetic operations. Our result is the first explicit
bound on the computational complexity of source identification of such
mixtures. The running time improves previous results by Feldman, O’Donnell,
and Servedio (FOCS 2005) and Chen and Moitra (STOC 2019) that guaranteed only
learning the mixture (without parametric identification of the source). Our
analysis gives a quantitative version of a qualitative characterization of
identifiable sources that is due to Tahmasebi, Motahari, and Maddah-Ali (ISIT
2018).

##  1 Introduction

The main result of this paper is an algorithm for source identification of
finite mixtures of binary product distributions. Our algorithm identifies the
source parameters of a mixture of $k$ product distributions on $n\geq 3k-3$
observable bits that satisfy a sufficient condition for identifiability, which
is that at least $3k-3$ of the observable bits are $\zeta$-separated
(definition below). The algorithm uses, as inputs, empirical frequencies that
are assumed to deviate from their true value by an error of at most
$\varepsilon<\zeta^{O(k^{2}\log k)}$. (This assumption will be satisfied, with
high probability, for a sample of size roughly $1/\varepsilon^{2}$.) The
algorithm identifies the parameters of the mixture model to an accuracy of
$\zeta^{-O(k^{2}\log k)}\varepsilon$. This is also roughly (up to a different
constant hidden by the big-Oh notation) the statistical distance between the
empirical distribution and the output (i.e., the learned) distribution. The
runtime of our algorithm is $2^{O(k^{2})}n^{O(k)}$ arithmetic operations.111We
note that the runtime relates to the post-sampling computation, after
aggregating the empirically observed frequencies. There are good reasons for
making this distinction. Collecting the sample and computing the frequencies
is computationally trivial. It can often be done under a streaming model or in
parallel. Or the frequencies might be available from an external source. 222In
this introduction we suppress dependence on mixture weights. Under a stronger
assumption that all observable bits are $\zeta$-separated, the runtime
improves to $2^{O(k^{2})}n$.

To formulate the problem, consider observable random variables
$X_{1},\ldots,X_{n}$ that are distributed on a common range $R$. In a finite
mixture model, the joint distribution on these random variables is governed by
a hidden or latent random variable $H$ supported on $\\{1,\ldots,k\\}$, such
that $X_{1},\ldots,X_{n}$ are statistically independent conditional on $H$. We
consider the case of a finite range $R$. The hardest and most fundamental case
is when the range is binary (i.e., $R=\\{0,1\\}$) and there are no further
constraints relating the distributions at the observables; the case of larger
$R$ reduces to this case (see [FOS08], for example).

Finite mixture models were pioneered in the late 1800s in [New86, Pea94] in
the context of applications in astronomy and the mathematical theory of
evolution. It is difficult to do justice to the vast literature in statistics
on mixture models; see, e.g., the surveys [EH81, TSM85, Lin95, MLR19]. The
computational complexity of learning mixture models was studied starting with
the seminal papers [KMR+94, CGG01, Das99, FM99]. The machine learning
community has shown recent interest in learning mixtures of product
distributions for pattern recognition. Motivating applications abound in
population genetics, bioinformatics, image recognition, text classification,
and other areas, e.g., [PSD00, JWL+05, JV04, JV02]. The algorithms in this
context are primarily based on the method of iterative Expectation
Maximization (EM) clustering (e.g., [JGHV04, LWPA16, PKM16, CPnR00]). As
detailed further below, in this genre provable guarantees of source
identification are provided in [NMR20], but the runtime depends exponentially
on the sample size, and a very large $n$ is required. In the theory of
computing literature, special cases and variants of learning mixtures of
product distributions were considered in [CGG01, FM99, CR08, AGM12, AFH+12,
AHK12, RSS14, LRSS15, KKM+19, CM19, GMRS20]. A more detailed account on the
previous results [FOS08, CM19] on mixture learning is given below.

The question of source identification dates back in the statistics literature
at least to 1950 [KR50, Koo50, Tei63, Bli64, YS68]. This is still a thriving
area of research (see, e.g., [CPnR00, AMR09, VS15, TMMA18, RVS20, ADXR20]).
From a computational perspective, the “learning” vs. “source identification”
contrast was raised in [FM99], in the context that interests us here of
learning a mixture of several (in that work, two) binary product
distributions. “Learning” means computing any hypothetical model that
generates observable statistics close to the empirical ones. “Source
identification” means computing a model that is close in parameter space to
the true model underlying the empirical statistics. Clearly, source
identification is a more challenging goal that implies also learning. Source
identification is desirable for a variety of reasons. A key reason is the
danger of overfitting a model to the empirical data, thereby ruining its
predictive guarantees (see, for instance, [KF09] ch. 16). Another reason is
that source identification is necessary for the quantification of causal
relations between the hidden variable and the observable variables (see
[Pea09, SGS00] for an introduction to graphical causal models)—here
identification, but not learning, provides explanatory value along with the
possibility of effective intervention. Finally, whether in the causal context
or otherwise, identification provides the statistician with semantics—an
actual characterization of the process generating the data—which may be far
more useful than the mere ability to artificially generate samples from the
same distribution.

The distinction between learning and source identification has two main
aspects. First, a distribution on observable variables might be inherently
explainable by two (or more) far-apart models, ruling out source
identification. Thus, we must restrict our attention to a class of
identifiable models, for which the mapping of model to distribution on
observables is one-to-one. In this paper, we consider the case that there are
at least $2k-1$ observables that are $\zeta$-separated. An observable bit $X$
is separated if the $k$ conditional probabilities $\Pr[X=1\mid H=j]$,
$j=1,\dots,k$ are mutually distinct, and $\zeta$-separated if every two values
differ by at least $\zeta>0$.333Notice that $\zeta$ cannot exceed
$\frac{1}{k-1}$. Having $2k-1$ observables that are separated is a sufficient,
and in general necessary condition for identification [TMMA18]. (Depending on
the parameters, fewer sources may suffice. But for example $2k-1$ are always
necessary when the $X_{i}$ are iid conditional on $H$, which we call the
“power distribution” case [RSS14], since terms of the form
$\prod_{i=1}^{n}\Pr[X_{i}=1\mid H=j]$ are replaced by $(\Pr[X_{1}=1\mid
H=j])^{n}$.)

Second, even if a model is identifiable in the limit of perfect statistics
(infinite sample size), the available empirical statistics might be
insufficiently accurate, yet still allow for the learning objective. It is
obvious, though, that a learning algorithm that runs on data produced by an
identifiable source and is required to achieve sufficiently high accuracy also
implicitly identifies the source.444This is because the space of source
parameters is compact and the mapping to observable statistics is continuous.
However, a learning algorithm might learn a mixture model without this
implying that the empirical data enables identification. In fact, the ground-
breaking work of [FOS08] gives an $(nk/\varepsilon)^{O(k^{3})}$ time algorithm
for learning a model of $k$-mixture of $n$ binary product distributions that
generates a distribution on the observables within statistical distance
$\varepsilon$ of the empirical statistics. A faster
$k^{O(k^{3})}(n/\varepsilon)^{O(k^{2})}$ time learning algorithm is given in
[CM19].555Notice, in particular, footnote 2 in that paper. The algorithms in
these papers are not guaranteed, under any assumptions, to identify the source
to any particular accuracy $\varepsilon$; in fact the algorithms succeed even
in cases where the source is not identifiable.

Our contribution is twofold. First, we prove a quantitative version of the
identifiability criterion: namely, for any given $\zeta>0$, we establish
(roughly; see later for exact statements) that any two models which differ by
$\eta$ in parameter space, differ in their statistics by at least
$\eta\zeta^{O(k^{2}\log k)}$. 666To simplify the informal discussion, we state
most of the bounds in the rest of this paragraph just for
$\eta\geq\zeta^{O(k^{2}\log k)}$. This, of course, implies that the algorithms
of [FOS08, CM19] can be used for source identification, assuming
$\zeta$-separation and sufficiently small target accuracy
$\varepsilon=\zeta^{O(k^{2}\log k)}$. Second, we improve substantially over
the runtime of these two algorithms, so that (even under the conditions under
which our result imples that those algorithms can perform identification), our
source identification algorithm is more efficient. Specifically, if we have at
least $3k-3$ observables that are $\zeta$-separated, our algorithm requires
empirical statistics accuracy $\zeta^{O(k^{2}\log k)}$ (i.e., sample size
$\zeta^{-O(k^{2}\log k)}$) and has a runtime of $2^{O(k^{2})}n^{O(k)}$
arithmetic operations. Our algorithm can also identify the source using the
minimum of $2k-1$ $\zeta$-separated observables; but then we require input
accuracy $\zeta^{O(k^{3})}$ (but the same $2^{O(k^{2})}n^{O(k)}$ runtime). If
all observables are $\zeta$-separated, then the runtime improves to
$2^{O(k^{2})}n$ 777Here we suppress a $O(\log\log(\varepsilon^{-1}))$ term.
(with the same input accuracies, according to the number of variables). Our
contributions establish quantitative bounds on the qualitative sufficiency of
$\zeta$-separation of [TMMA18]; the results in that paper are entirely non-
algorithmic. The only explicit algorithmic result on source identification
that we are aware of is [NMR20]. Their algorithm, which runs under a somewhat
more general assumption than $\zeta$-separation, requires complete enumeration
over the choice of mixture constituent that generated each sample point, for a
sufficiently large sample. Thus, it is prohibitively expensive, requiring at
least $n=\exp(k^{2})$ observable random variables (as compared with $3k-3$),
and runtime that is doubly exponential, namely $k^{\exp(k^{2})}$.

The main idea underlying our algorithm is the following. Given sufficiently
many $\zeta$-separated observables, for which we have sufficiently accurate
empirical multilinear moments,888We call these empirical moments as we expect
them to obtained by sampling; but our theorems depend only on their being
sufficiently accurate. we show how to construct “synthetic bits” for which we
can compute highly accurate power moments, i.e., the moments that occur in the
far more restricted problem of power distributions. The higher moments of
these synthetic bits are created out of linear combinations of multilinear
moments of the original bits. This mechanism in its idealized form (i.e., for
perfect statistics) suffices to re-prove the theorem of [TMMA18] that $2k-1$
separated observables suffice for source identification. The next challenge we
face is to bound the coefficients of the multilinear monomials in these linear
combinations, as this affects the required accuracy of the empirical
statistics (thus, the required sample size). The synthetic bits method reduces
the problem to the special case of identification of a mixture of $k$ power
distributions, i.e., when the $X_{i}$-s are iid conditional on $H$. This
special case is effectively an extension of the theory of orthogonal
polynomials on the reals and the classical moment problem [Sch17, Sim15], and
methods such as Prony’s method or the Matrix Pencil method were shown to solve
it [RSS14, LRSS15, KKM+19, GMRS20]. Despite the power distributions case being
so highly constrained, it has useful applications, e.g., in reconstructing
population histories and in learning topic models (see the above references).
We can use these existing algorithms to recover the mixture of power
distributions on synthetic bits. That in turn enables recovery of the mixture
of product distributions on the observable bits. The best algorithm for power
distributions to date [GMRS20] requires estimates of the first $2k$ moments of
the synthetic bits to within accuracy $\zeta^{O(k)}$, and has runtime
$O(k^{2+o(1)})$. Thus, this component of our runtime is cheap; the runtime
bound of our algorithm is dominated by the exhaustive search for $3k-3$
observables that are $\zeta$-separated (unless all observables are known to be
$\zeta$-separated), and by the construction of the $2k$ synthetic bits.

##  2 Preliminaries and main theorem

##### Notation

There is a hidden variable $H$ ranging in $[k]$, and $n$ binary observable
variables $X_{i}$; we write ${\mathbf{m}}_{ij}:=\Pr[X_{i}=1\mid H=j]$.

Vectors are row vectors unless otherwise indicated. For $S\subseteq[n]$ define
the random variable $X_{S}=\prod_{i\in S}X_{i}$. We make extensive use of
Hadamard product for vectors $u=(v_{1},\ldots,v_{k})$,
$v=(v_{1},\ldots,v_{k})$:

$\displaystyle\odot:\mathbb R^{k}\times\mathbb R^{k}$ $\displaystyle\to\mathbb
R^{k}$ $\displaystyle u\odot v$ $\displaystyle=(u_{1}v_{1},\ldots,u_{k}v_{k})$

The identity for this product is the all-ones vector ${\mathbb{1}}$. We
associate with vector $u$ the linear operator
$u_{\odot}=\operatorname{diag}(u)$, a $k\times k$ diagonal matrix, so that

$v\cdot u_{\odot}=u\odot v.$

Let ${\mathbf{m}}_{i}$ be the row vector
$({\mathbf{m}}_{i1},\ldots,{\mathbf{m}}_{ik})$. Let
${\mathbf{M}}\in\mathbb{R}^{2^{[n]}\times k}$ be the matrix with rows indexed
by subsets $S=\left\\{i_{1},\dotsc,i_{s}\right\\}\subseteq[n]$, with rows

${\mathbf{M}}_{S}={\mathbf{m}}_{i_{1}}\odot{\mathbf{m}}_{i_{2}}\odot\dotsb\odot{\mathbf{m}}_{i_{s}}.$

In particular, ${\mathbf{M}}_{\emptyset}={\mathbb{1}}$ and
${\mathbf{M}}_{\left\\{i\right\\}}={\mathbf{m}}_{i}$ for all $i\in[m]$.

Observe that source identification is not possible if ${\mathbf{M}}$ has less
than full column rank, i.e., $\operatorname{rank}{\mathbf{M}}<k$, as then the
mixing weights cannot be unique.

For a collection of subsets $\mathcal{S}\subseteq 2^{[n]}$, let
${\mathbf{M}}[\mathcal{S}]$ denote the restriction of ${\mathbf{M}}$ to the
rows ${\mathbf{M}}_{S},S\in\mathcal{S}$. E.g.,
${\mathbf{M}}={\mathbf{M}}[2^{[n]}]$.

##### The empirical multi-linear moments

For a finite sample drawn from the model, we let
$\operatorname{{\tilde{g}}}(S)$ be the empirical estimate of
$\operatorname{\mathbb{E}}[X_{S}]$, i.e., the fraction of samples for which
$\prod_{i\in S}X_{i}=1$. These $\operatorname{{\tilde{g}}}(S)$ for
$S\subseteq[n]$ are the complete list of “observables” of the model. Each
converges, in the infinite-sample limit, to the value
$\operatorname{g}(S):=\operatorname{\mathbb{E}}[X_{S}]$,

$\operatorname{g}(S)={\mathbf{M}}_{S}\pi^{\mathsf{T}}={\mathbf{M}}_{S}\pi_{\odot}{\mathbb{1}}^{\mathsf{T}}.$

##### Comparison to the iid case

The fact that we have access only to multi-linear moments is the key
constraint of the problem. Compare with the “power” case, i.e., when we know
in advance that all rows ${\mathbf{m}}_{i}$ are identical. (That is, the
$X_{i}$ are iid conditional on $H$.) Then observing bit $X_{2}$ alone tells us
nothing about $X_{1}$ that we do not learn from observing bit $X_{1}$. But
observing the product $X_{1}X_{2}$ does tell us something new: in fact it is
distributed as would be a binary observable whose row had entries
$({\mathbf{m}}_{11}^{2},\ldots,{\mathbf{m}}_{1k}^{2})$. With $2k$ rows, we
obtain the first $2k$ moments of the $k$-sparse distribution corresponding to
${\mathbf{m}}_{1}$ (i.e., the real-valued, $k$-sparse distribution which
places atomic probability $\pi_{j}$ at ${\mathbf{m}}_{1j}\in\mathbb R$, also
called a “$k$-spike” distribution). From this point on, one may apply the
time-honored method of Prony to identify the source. For details, a runtime
analysis, and further references, see [GMRS20].

Much of the interest of the present problem, by contrast, is due precisely to
the fact that we cannot read out higher moments of the distributions
corresponding to any of the bits $X_{i}$, because the various rows do not have
any assumed relationship. How to nonetheless obtain higher moments of
individual rows, is the challenge our algorithm tackles.

##### Main theorem

In what follows $\zeta$ is an assumed separation parameter, and $\pi_{\min}$
is an assumed lower bound on mixture weights.

######  Theorem 1.

($i$) Given access to the joint statistics of $n$ observable bits among which
at least $3k-3$ which are $\zeta$-separated, with all statistics available to
additive accuracy $\varepsilon$ for $\varepsilon\leq(\pi_{\min})^{O(\log
k)}\zeta^{O(k^{2}\log k)}$, our algorithm runs in time $2^{O(k^{2})}n^{O(k)}$
and computes the model parameters ($\pi$ and all row values
${\mathbf{m}}_{ij}$) to within accuracy $\varepsilon\zeta^{-O(k^{2}\log
k)}(\pi_{\min})^{-O(\log k)}$.

($ii$) If all rows are $\zeta$-separated, then a slightly simpler version of
our algorithm identifies the source (to within the same accuracy, given the
same input accuracy, as in ($i$)), in runtime $2^{O(k^{2})}n$ (i.e., linear in
$n$).

($iii$) If only $2k-1$ $\zeta$-separated rows are available, another simpler
version of our algorithm also identifies the source, but the loss factor on
the accuracy is $(\pi_{\min})^{-k}\zeta^{-O(k^{3})}$ and consequently one must
start with ${\varepsilon}\leq(\pi_{\min})^{k}\zeta^{O(k^{3})}$ which requires
sample complexity comparable to prior work, but achieving the same improved
runtime as in ($i$).

##  3 Algorithm

##### Further definitions

Let $V_{(i)}$ be the subspace of $\mathbb R^{k}$ spanned by the $i$’th
standard basis vector, and $P_{(i)}$ the projection onto $V_{(i)}$ w.r.t.
usual inner product (i.e., as a matrix, all $0$’s except a $1$ in entry
$(i,i)$). For a separated vector $v=(v_{1},\ldots,v_{k})$ define the Lagrange
interpolation polynomials $p_{v,i}$ by $p_{v,i}(v_{j})=\delta_{ij}$. We then
have the matrix equation

$p_{v,i}(v_{\odot})=P_{(i)}.$

Write $p_{v,i}(x)=\sum_{j=0}^{k-1}p_{v,i,j}x^{j}$. Define the “coefficient
norm” of a polynomial by $\|p_{v,i}\|_{c}=\sum_{0}^{k-1}j|p_{v,i,j}|$.

It will be very useful to put our observables in matrix form. Let
$S,T\subseteq[n]$ be disjoint sets and take any $\mathcal{A}\subseteq
2^{S},\mathcal{B}\subseteq 2^{T}$. Then the matrix
$\mathbf{C}_{\mathcal{B}\mathcal{A}}$ is observable (meaning every entry of it
is a function of the joint statistics of the observable random variables
$X_{1},\ldots,X_{n}$), where

$\mathbf{C}_{\mathcal{B}\mathcal{A}}\coloneqq{\mathbf{M}}[\mathcal{B}]\pi_{\odot}{\mathbf{M}}[\mathcal{A}]^{\top}.$
(1)

Let
${\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\coloneqq\left[\operatorname{{\tilde{g}}}(B\cup
A)\right]_{B\in\mathcal{B},A\in\mathcal{A}}$ be the corresponding matrix of
empirical moments.

We’ll use $\left\|\cdot\right\|$ to denote the operator norm for matrices
(with respect to the Euclidean norm in the domain and range).

##### Strategy

If ${\mathbf{m}}$ possesses $2k-1$ or $3k-3$ (depending on the version of the
algorithm; this affects only runtime) $\zeta$-separated rows of
${\mathbf{m}}$, and given empirical statistics within the required distance
from perfect statistics, our algorithm will identify the model (to within
similar accuracy). This includes even rows which are not themselves
$\zeta$-separated; all we need is that some $2k-1$ or $3k-3$ rows be
$\zeta$-separated. The algorithm has the following structure: range over all
$n^{O(k)}$ subsets of the rows; run the identification algorithm using that
set. If the set does not suffice for identification (which can happen only if
the set includes some non-$\zeta$-separated rows), this will be flagged by the
algorithm. Any two such runs which do terminate successfully, must result in
very close parameter reconstructions.

###  3.1 Synthetic bits and bootstrapping

We start with the $2k-1$-observables version of the algorithm. Our final
algorithm in section 3.1.2, which uses $3k-3$ observables and achieves better
runtime, will use a slightly more complicated construction, but the main ideas
are present in the simpler variant here.

####  3.1.1 Constructing higher moments of a row, using $2k-2$ other rows

In what follows we show how to compute moments of arbitrary degree of the
$k$-spike distribution associated with any bit $X_{i}$, given access to any
additional $2k-2$ $\zeta$-separated observable bits. For concreteness, let
$X_{1}$ (corresponding to the row ${\mathbf{m}}_{1}$) be the variable for
which we want to find higher moments. (We don’t require that
${\mathbf{m}}_{1}$ be $\zeta$-separated for the moment computation, although
that will be needed in a subsequent step of the algorithm.) Let
$S=\\{2,\ldots,k\\}$, $T=\\{k+1,\ldots,2k-1\\}$ be the indices of the $2k-2$
other $\zeta$-separated rows, partitioned into two sets of $k-1$ rows each.

The only thing that our statistics tell us about row ${\mathbf{m}}_{1}$,
alone, is its first moment:
$\operatorname{\mathbb{E}}(X_{1})={\mathbf{m}}_{1}\pi_{\odot}{\mathbb{1}}^{\top}$.
Equivalently this quantity is also the expectation the $k$-spike distribution
associated with $X_{1}$. It will be critical to obtain higher moments of this
distribution. The second moment is equal to
$({\mathbf{m}}_{1}\odot{\mathbf{m}}_{1})\pi_{\odot}{\mathbb{1}}^{\top}$, and
more generally (with ${\mathbf{m}}_{1}^{\odot r}$ denoting the $r$-fold
Hadamard product of ${\mathbf{m}}_{1}$ with itself), the $r$’th moment is
${\mathbf{m}}_{1}^{\odot r}\pi_{\odot}{\mathbb{1}}^{\top}$ for any $r$. We
will need to have the $1,\dotsc,2k$’th moments in order to solve for
${\mathbf{m}}_{1}$ and $\pi$.

We show in Section 5 that there exist subsets
$\mathcal{A}=\left\\{A_{1},\dotsc,A_{k}\right\\}\subseteq 2^{S}$,
$\mathcal{B}=\left\\{B_{1},\dotsc,B_{k}\right\\}\subseteq 2^{T}$ of size $k$
each such that
${\mathbf{A}}\coloneqq{\mathbf{M}}[\mathcal{A}]\in\mathbb{R}^{k\times k}$ and
${\mathbf{B}}\coloneqq{\mathbf{M}}[\mathcal{B}]\in\mathbb{R}^{k\times k}$ are
invertible. Moreover, we’ll have $A_{1}=\emptyset=B_{1}$ so that the first row
of ${\mathbf{M}}[A]$ and ${\mathbf{M}}[B]$ will be the vector ${\mathbb{1}}$.
Then the matrix
$\mathbf{C}_{\mathcal{B}\mathcal{A}}={\mathbf{B}}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$
(as defined in (1)) is an invertible, observable matrix.

Now consider the vector

$v_{1}\coloneqq{\mathbf{m}}_{1}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}=(\operatorname{\mathbb{E}}[X_{1}X_{A_{1}}],\dotsc,\operatorname{\mathbb{E}}[X_{1}X_{A_{k}}]).$

Each coordinate
$\operatorname{\mathbb{E}}[X_{1}X_{A_{i}}]=\operatorname{\mathbb{E}}[X_{A_{i}\cup\left\\{1\right\\}}]$
is a multi-linear moment and is therefore observable. In the algorithm, our
starting point will be approximations to multi-linear moments, and we will
address that shortly, but for now we write down the algorithm as it would run
with perfect statistics.

While we’d like to find ${\mathbf{m}}_{1}^{\odot 2}$ directly, we won’t be
able to do so, and instead we’ll settle for computing the vector $u_{1}$
satisfying

$u_{1}{\mathbf{B}}\coloneqq{\mathbf{m}}_{1}.$

This $u_{1}$ is unique since ${\mathbf{B}}$ is invertible, and it gives the
coefficients needed to express ${\mathbf{m}}_{1}$ as a linear combination of
the rows of ${\mathbf{B}}$. We can in fact compute $u_{1}$ from observables,
because

$u_{1}=v_{1}\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}.$ (2)

We can verify (2) algebraically:

$v_{1}(\mathbf{C}_{\mathcal{B}\mathcal{A}})^{-1}{\mathbf{B}}=v_{1}(A^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}{\mathbf{B}}^{-1}{\mathbf{B}}=v_{1}({\mathbf{A}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}={\mathbf{m}}_{1}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}({\mathbf{A}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}={\mathbf{m}}_{1}.$
(3)

Conceptually what (3) means is that $u_{1}{\mathbf{B}}$, a linear combination
of rows of ${\mathbf{B}}$, defines a random variable that has been synthesized
out of the $X_{B_{\ell}}$ random variables, and shares the same expectations
as $X_{1}$ conditional on any setting of $H$. We call this the method of
synthetic bits. Our algorithm consists of repeatedly bootstrapping synthetic
bits.

To be more explicit, use $u_{1}=(u_{11},\dotsc,u_{1k})$ to define a new random
variable $Y\coloneqq\sum_{\ell=1}^{k}u_{1\ell}X_{B_{\ell}}$.

$Y$ has been synthesized out of $X_{B_{1}},\dotsc,X_{B_{k}}$. $Y$ has the same
expectation as $X_{1}$, conditioned on the value of the hidden variable $H$,
since

$\operatorname{\mathbb{E}}(X_{1}\mid
H=j)=({\mathbf{m}}_{1})_{j}=\sum_{\ell=1}^{k}u_{1\ell}\operatorname{\mathbb{E}}(X_{B_{\ell}}\mid
H=j)=\operatorname{\mathbb{E}}(Y\mid H=j).$

Moreover, given the value of $H$, $X_{1}$ and $Y$ are independent, because
$\\{1\\}\cap B_{\ell}=\emptyset$.

##### Bootstrapping to obtain the second moment.

As a consequence of (3), we can perform the bootstrapping step which is at the
heart of our algorithm. We are interested in obtaining the vector

$v_{2}\coloneqq{\mathbf{m}}_{1}^{\odot 2}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$

which is a linear image (under the mapping
$\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$) of a random variable distributed as
the product of two random variables which are independent conditional on $H$;
each distributed as $X_{1}$ conditional on $H$. This we get by setting

$v_{2}\coloneqq(\operatorname{\mathbb{E}}(X_{1}YX_{A_{1}}),\dotsc,\operatorname{\mathbb{E}}(X_{1}YX_{A_{k}})),$

Since $A_{1}=\emptyset$, the entry
$(v_{2})_{1}=\operatorname{\mathbb{E}}(X_{1}Y)={\mathbf{m}}_{1}^{\odot
2}\pi_{\odot}{\mathbb{1}}^{\mathsf{T}}$ is exactly our desired second moment.
To get access to $v_{2}$, we observe that even though we don’t necessarily
have two independent copies of ${\mathbf{m}}_{1}$ among our rows, our
synthetic bit $Y$ provides the needed independence from $X_{1}$. As a matter
of notation, for our collection $\mathcal{B}\subseteq 2^{T}$, we define
$\mathcal{B}+\\{1\\}$ to consist of the sets $B_{\ell}\cup\\{1\\}$ for each
$B_{\ell}\in\mathcal{B}$.

Now we can write

$v_{2}={\mathbf{m}}_{1}^{\odot
2}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}=({\mathbf{m}}_{1}\odot(u_{1}{\mathbf{B}}))\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}=u_{1}\mathbf{C}_{\mathcal{B}+\\{1\\},\mathcal{A}}$
(4)

Since we already have $u_{1}$, and since
$\mathbf{C}_{\mathcal{B}+\\{1\\},\mathcal{A}}$ is observable, expression (4)
can be used to compute $v_{2}$.

##### Bootstrapping to all moments.

The generalization of the second moment computation is this. Given vector
$u_{r-1}$ defined by
$u_{r-1}\coloneqq{\mathbf{m}}_{1}^{\odot(r-1)}{\mathbf{B}}^{-1}$, we define
$v_{r}$ and $u_{r}$ and show how to compute them:

Definition Computation $\displaystyle v_{r}$
$\displaystyle\coloneqq{\mathbf{m}}_{1}^{\odot
r}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$
$\displaystyle\quad=({\mathbf{m}}_{1}\odot(u_{r-1}{\mathbf{B}}))\pi_{\odot}{\mathbf{A}}^{\top}$
$\displaystyle=u_{r-1}\mathbf{C}_{\mathcal{B}+\\{1\\},\mathcal{A}}$ (5)
$\displaystyle u_{r}$ $\displaystyle\coloneqq{\mathbf{m}}_{1}^{\odot
r}{\mathbf{B}}^{-1}$ $\displaystyle\quad={\mathbf{m}}_{1}^{\odot
r}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}{{\mathbf{A}}^{\mathsf{T}}}^{-1}\pi_{\odot}^{-1}{\mathbf{B}}^{-1}$
$\displaystyle=v_{r}\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}.$ (6)

In the actual algorithm, we’ll be working with empirical approximations of
$v_{r},u_{r}$, and $\mathbf{C}_{\mathcal{B}\mathcal{A}}$, which we’ll denote
$\tilde{v}_{r}$, $\tilde{u}_{r}$, and
${\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}$, respectively. A key part of
the technical work will be in bounding error amplification. In order to
compute $\tilde{v}_{r}$ and $\tilde{u}_{r}$, we’ll use (following (5),(6)) the
following assignments, where all quantities are empirical estimates:

$\displaystyle\tilde{v}_{r}$
$\displaystyle\coloneqq\tilde{u}_{r-1}{\tilde{\mathbf{C}}}_{\mathcal{B}+\\{1\\},\mathcal{A}}$
(7) $\displaystyle\tilde{u}_{r}$
$\displaystyle:=\tilde{v}_{r}({\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}})^{-1}$
(8)

It is important to note that the bootstrapping can be performed only because
$\tilde{u}_{r-1}$ places weight solely upon rows in $\mathcal{B}\subseteq
2^{T}$. Then the bootstrapping yields an expression for $\tilde{v}_{r}$ that
contains moments involving $X_{1}$ as well as the entries in $T$. But once we
compute $\tilde{u}_{r}$, we regain a coefficient vector for a synthetic bit,
that again places weight solely upon rows in $\mathcal{B}$.

The fact that both the computation of $\tilde{v}_{r}$ (in (7)) and of
$\tilde{u}_{r}$ (in (8)) work stably and not only in the perfect-statistics
limit, relies upon the following:

######  Corollary 2.

If all empirical multilinear moments are within
$\varepsilon<\zeta^{\Omega(k^{2})}$ of their true values, then
$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}-\mathbf{C}_{\mathcal{B}\mathcal{A}}\right\|\leq
k\varepsilon$ and
$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}-{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}\right\|\leq\zeta^{-O(k^{2})}\pi_{\min}^{-2}\varepsilon$.

###### Proof.

This will be an immediate consequence of Lemma 10. ∎

From this, we can bound the increase in error due to each subsequent
application of (7) and (8).

######  Lemma 3.

If all empirical multilinear moments are within
$\varepsilon<\zeta^{\Omega(k^{2})}$ of their true values, then for all $i$,
$\left\|\tilde{u}_{i}-u_{i}\right\|_{\infty}<\zeta^{-O(k^{2})}\pi_{\min}^{-2}\left\|\tilde{u}_{i-1}-u_{i-1}\right\|_{\infty}$
and
$\left\|\tilde{v}_{r}-v_{r}\right\|\leq\zeta^{-O(k^{2})}\pi_{\min}^{-2}\left\|\tilde{v}_{i-1}-v_{i-1}\right\|_{\infty}$.

###### Proof.

This follows from Lemma 11, which is itself a consequence of Lemma 10. ∎

Nevertheless, if we compute each $\tilde{v}_{r}$ and $\tilde{u}_{r}$ as
described above, we would need to start with accuracies
${\varepsilon}<\zeta^{\Omega(k^{3})}$ in order to retain accuracy after
bootstrapping $O(k)$ times, as is required by the above algorithm. To avoid
this, we’ll need to reduce the number of iterations by using $3k-3$, rather
than $2k-1$, $\zeta$-separated rows.

####  3.1.2 Improved error control: constructing higher moments using $3k-3$
$\zeta$-separated rows

In order to avoid performing $k$ iterations to compute $\tilde{v}_{k}$ and
$\tilde{u}_{k}$, we’ll use another set of rows,
$\mathcal{B}^{\prime}=\left\\{B^{\prime}_{1},\dotsc,B^{\prime}_{k}\right\\}\subseteq
2^{T^{\prime}}$ of size $k$, where $T^{\prime}$ is disjoint from $S$ and $T$,
and $B^{\prime}_{1}=\emptyset$. Now we’ll introduce
${\mathbf{B}}^{\prime}\coloneqq{\mathbf{M}}[{\mathbf{B}}^{\prime}]$ and
$\mathbf{C}_{\mathcal{B}^{\prime}\mathcal{A}}\coloneqq{\mathbf{B}}^{\prime}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$,
both of which will be invertible as before.

Previously, $u_{i}$ was a linear combination of the rows of ${\mathbf{B}}$
such that $u_{i}{\mathbf{B}}={\mathbf{m}}_{1}^{\odot i}$. We’ll introduce a
new, but similar, sequence of vectors $u^{\prime}_{i}$ where
$u^{\prime}_{i}{\mathbf{B}}^{\prime}={\mathbf{m}}_{1}^{\odot i}$ and
$u^{\prime}_{i}$ is obtained from $v_{i}$ by
$u^{\prime}_{i}=v_{i}\mathbf{C}_{\mathcal{B}^{\prime}\mathcal{A}}^{-1}$. In
the algorithm, we’ll only have access to the approximations
$\tilde{u}^{\prime}_{i}$ and
${\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}}$ and we’ll compute
each successive $\tilde{u}^{\prime}_{i}$ by
$\tilde{u}^{\prime}_{i}=\tilde{v}_{i}{\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}}^{-1}$.

The advantage we obtain over the more straightforward process of the previous
section, results from our ability to compute $v_{2^{i}}$ using $u_{2^{i-1}}$
and $u^{\prime}_{2^{i-1}}$; in that way, we are able to get away with
performing only $1+\lg k$ iterations to compute any of $v_{1},\dotsc,v_{2k}$.

To describe the computation, we first define the sum of two collections of
subsets $\mathcal{U},\mathcal{V}\subseteq 2^{[n]}$ by

$\mathcal{U}+\mathcal{V}\coloneqq\left\\{U\cup
V:U\in\mathcal{U},V\in\mathcal{V}\right\\}.$

Now let $x$ and $y$ be vectors indexed by the the subsets in $\mathcal{B}$ and
$\mathcal{B}^{\prime}$, respectively. Recall that $\mathcal{B}\subseteq
2^{T},\mathcal{B}^{\prime}\subseteq 2^{T^{\prime}}$, where $T\cap
T^{\prime}=\emptyset$. Then
$\left|\mathcal{B}+\mathcal{B}^{\prime}\right|=k^{2}$ and each subset in
$\mathcal{B}+\mathcal{B}^{\prime}$ can be uniquely written as $B_{\ell}\cup
B^{\prime}_{j}$ for $\ell,j\in[k]$. We define the Kronecker product $(x\otimes
y)\in\mathbb{R}^{\mathcal{B}+\mathcal{B}^{\prime}}$ to be the vector indexed
by subsets in $\mathcal{B}+\mathcal{B}^{\prime}$ given by

$(x\otimes y)_{B_{\ell}\cup B^{\prime}_{j}}\coloneqq
x_{B_{\ell}}y_{B^{\prime}_{j}}$

for any $\ell,j\in[k]$. To access $v_{2^{i}}$, we write

$v_{2^{i}}={\mathbf{m}}_{1}^{\odot
2^{i}}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}=((u_{2^{i-1}}{\mathbf{B}})\odot(u^{\prime}_{2^{i-1}}{\mathbf{B}}))\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}=(u_{2^{i-1}}\otimes
u^{\prime}_{2^{i-1}})\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}},$

expressing the row $v_{2^{i}}$ as the linear of combination of the $k^{2}$
rows corresponding to the subsets in $\mathcal{B}+\mathcal{B}^{\prime}$.

Of course, we’ll need to compute $v_{\ell}$ for $\ell$ not a power of $2$. We
can do this using a slight modification of the recursive procedure where for
$i=1,\dotsc,\lg k+1$ and $j=1,\dotsc,2^{i}$, we compute

$\tilde{v}_{\ell}\coloneqq(\tilde{u}_{j}\otimes\tilde{u}^{\prime}_{2^{i}})\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}$
(9)

where $\ell=j+2^{i}$, and we’ve computed $\tilde{u}_{j}$ in a prior iteration
for all $j\leq 2^{i}$.

Under this modification, each $\tilde{v}_{i}$ is produced in only $\log k$
iterations each of which involves a matrix multiplication by
${\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}$ or a convolution followed
by multiplication by
$\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}$, each step of
which can increase the error in total by $\zeta^{-O(k^{2}\log
k)}\pi_{\min}^{-2}$. By starting with empirical moments accurate to within
$\zeta^{\Omega(k^{2}\log k)}\pi_{\min}^{2k}$, we can ensure that the resulting
vectors $\tilde{v}_{i}$ and $\tilde{u}_{i}$ are sufficiently close to start
solving for ${\mathbf{m}}_{1}$ and $\pi$.

We provide pseudocode for this “$3k-3$ rows” version of the algorithm, see
Fig. 1.

##### Flagging a failure condition

If the chosen rows $S\cup T\cup T^{\prime}$ fail to all be $\zeta$-separated,
the algorithm might fail. However, we will detect such failure. The conditions
that we actually need so that the algorithm should work, are these: (a)
${\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}$ and
${\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime},\mathcal{A}}$ should have a large
least singular value. (It does not actually matter whether all rows we use are
$\zeta$-separated, that was merely a sufficient condition for this well-
conditioning.) We compute this singular value explicitly and simply dismiss
the triple $S,T,T^{\prime}$ if this condition fails. (b) The first row of
$\mathcal{B}$, namely $\\{1\\}$ in the numbering used in the pseudocode for
Algorithm 1, should be $\zeta$-separated. If condition (a) holds but this
condition fails, we will detect the failure in line 13 of the algorithm,
because the Hankel matrix will have insufficient eigenvalue gap (see Cor. 12
of [GMRS20]).

###  3.2 Solving the power distribution problem

Once we’ve computed $\tilde{v}_{1},\dotsc,\tilde{v}_{2k}$, we have access to
all of the moments of the $k$-spike distribution corresponding to observable
$X_{1}$. (The $i$th moment is the first entry of the corresponding vector
$\tilde{v}_{i}$). Recall, these are the moments of a mixture of $k$ Bernoulli
random variables, where the $r$’th moment corresponds to drawing a mixture
component $j$ with probability $\pi_{j}$, then setting the Bernoulli random
variable to $1$ with probability ${\mathbf{m}}_{1j}^{r}$. The problem of
recovering the parameters (i.e., the vectors ${\mathbf{m}}_{1}$ and $\pi$)
from approximate moments of this form has been extensively studied, and many
algorithms have been provided. We use the algorithm
$\operatorname{\textsc{LearnPowerDistribution}}$ from [GMRS20], which on
inputs accurate to within $\varepsilon$, outputs parameters
$\tilde{{\mathbf{m}}}_{1}$ and $\tilde{\pi}$ to within accuracy
$\frac{1}{\pi_{\min}}\zeta^{-O(k)}\varepsilon$ while running in time
(arithmetic operations) $k^{2+o(1)}$.

###  3.3 Recovering the remaining parameters

Once we have estimates for ${\mathbf{m}}_{1}$ and $\pi$, we can simply solve
for ${\mathbf{A}}$, ${\mathbf{B}}$, and ${\mathbf{B}}^{\prime}$ using the fact
that
${\mathbf{A}}^{\mathsf{T}}=\pi_{\odot}^{-1}(\operatorname{Vdm}({\mathbf{m}}_{1}))^{-1}V$
where

$\tilde{V}=(v_{0};\dotso;v_{k-1})$

is the matrix with rows $v_{i}$ and $\operatorname{Vdm}({\mathbf{m}}_{1})$ is
the Vandermonde matrix with rows ${\mathbf{m}}_{1}^{\odot i}$ for
$i=0,\dotsc,k-1$. Note that here we finally do require that ${\mathbf{m}}_{1}$
be $\zeta$-separated. We can thus solve for ${\mathbf{A}}^{\mathsf{T}}$.

To solve for ${\mathbf{B}}$ we use
${\mathbf{B}}=\mathbf{C}_{\mathcal{B}\mathcal{A}}({\mathbf{A}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}$.
Likewise for ${\mathbf{B}}^{\prime}$ we use
${\mathbf{B}}^{\prime}=\mathbf{C}_{\mathcal{B}^{\prime}\mathcal{A}}({\mathbf{A}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}$.

Now for any row $i$ not already computed, we need only pick any other set
$\mathcal{S}=\left\\{S_{1},\dotsc,S_{k}\right\\}$ of $k$ linearly independent
rows supported on a set not containing $i$, and we can solve for
${\mathbf{m}}_{i}$ by writing

${\mathbf{m}}_{i}=(\operatorname{\mathbb{E}}[X_{i}X_{S_{1}}],\dotsc,\operatorname{\mathbb{E}}[X_{i}X_{S_{k}}]){{\mathbf{M}}[\mathcal{S}]^{\mathsf{T}}}^{-1}\pi_{\odot}^{-1}$

In particular, by setting $\mathcal{S}=\mathcal{A}$ we can solve for all rows
in $[n]\setminus S$ and by setting $\mathcal{S}=\mathcal{B}$ we can solve for
all rows in $[n]\setminus T$. Together, this suffices to solve for all rows.

###  3.4 Runtime

The algorithm contains three main parts:

  1. 1.

Find disjoint $S,T,T^{\prime}\subset[n]$ and $\mathcal{A}\subset 2^{S}$,
$\mathcal{B}\subset 2^{T}$, and $\mathcal{B}^{\prime}\subset 2^{T^{\prime}}$,
which is complexity $n^{O(k)}2^{O(k^{2})}$. First we require $n^{O(k)}$
iterations to check all possible disjoint $S,T,T^{\prime}$. Then,
$2^{O(k^{2})}$ operations are required in each iteration to check all size $k$
subsets of the $2^{k-1}$ rows of ${\mathbf{M}}[2^{S}],{\mathbf{M}}[2^{T}]$ and
${\mathbf{M}}[2^{T^{\prime}}]$.

  2. 2.

Nested loops to compute higher order moments $\tilde{v}$, and the
corresponding $\tilde{u}$. This step takes time $O(\text{poly}(k))$.

  3. 3.

Applying the power distribution result. This can be done in time
$O(k^{2+o(1)}+k(\log^{2}k)\log\log(\varepsilon^{-1}))$ (see Corollary 15).

This gives runtime complexity of
$n^{O(k)}2^{O(k^{2})}+O(k^{2+o(1)}+k(\log^{2}k)\log\log(\varepsilon^{-1}))$.
If all sources are $\zeta$-separated, we do not need to iterate over choices
of $S,T,T^{\prime}$, so the runtime improves to
$2^{O(k^{2})}+O(k^{2+o(1)}+k\log^{2}(k)\log\log(\varepsilon^{-1}))$.

Algorithm 1 Identifies a mixture of product distributions given $3k-3$
$\zeta$-separated observable bits

1:Let $\mathcal{B}\subseteq
2^{\left\\{1,\dotsc,k-1\right\\}},\mathcal{B}^{\prime}\subseteq
2^{\left\\{k,\dotsc,2k-2\right\\}},\mathcal{A}\subseteq
2^{\left\\{2k-1,\dotsc,3k-3\right\\}}$, with
$\left|\mathcal{A}\right|=\left|\mathcal{B}\right|=\left|\mathcal{B}^{\prime}\right|=k$
maximizing
$\min\\{\sigma_{k}({\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}),\sigma_{k}({\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}})\\}$.
If this min is below $\pi_{\min}\zeta^{O(k^{2})}$, terminate. Denote
$\mathcal{B}=\left\\{B_{1},\dotsc,B_{k}\right\\}$,
$\mathcal{B}^{\prime}=\left\\{B^{\prime}_{1},\dotsc,B^{\prime}_{k}\right\\}$,
and $\mathcal{A}=\left\\{A_{1},\dotsc,A_{k}\right\\}$. Without loss of
generality $B_{1}=\\{1\\}$.

2:$\tilde{v}_{0}\leftarrow(\operatorname{{\tilde{g}}}(A_{1});\dotso;\operatorname{{\tilde{g}}}(A_{k}))$.

3:$\tilde{v}_{1}\leftarrow(\operatorname{{\tilde{g}}}(A_{1}\cup\left\\{1\right\\});\dotso;\operatorname{{\tilde{g}}}(A_{k}\cup\left\\{1\right\\}))$.

4:$\tilde{u}_{1}\leftarrow\tilde{v}_{1}({\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}})^{-1}$.

5:$\tilde{u}^{\prime}_{1}\leftarrow\tilde{v}_{1}({\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}})^{-1}$.

6:for $i=2,\dotsc,\log k+1$ do

7: for $j=1,\dotsc,2^{i-1}$ do

8:
$\tilde{v}_{2^{i-1}+j}\leftarrow(\tilde{u}_{j}\otimes\tilde{u}^{\prime}_{2^{i-1}}){\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}$

9:
$\tilde{u}_{2^{i-1}+j}\leftarrow\tilde{v}_{2^{i-1}+j}({\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}})^{-1}$.

10: end for

11:
$\tilde{u}^{\prime}_{2^{i}}\leftarrow\tilde{v}_{2^{i}}({\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}})^{-1}$.

12:end for

13:Let $H_{k+1}$ be the $(k+1)\times(k+1)$ Hankel matrix with entries given by
$[H_{k+1}]_{i,j=0}^{k}=(\tilde{v}_{i+j})_{1}$. If the second-smallest
eigenvalue of $H_{k+1}$ is below $\frac{\pi_{\min}}{2}(\zeta/16)^{2k-2}$,
terminate.

14:$\tilde{{\mathbf{m}}}_{1},\tilde{\pi}\leftarrow\operatorname{\textsc{LearnPowerDistribution}}(H_{k+1})$.

15:$\tilde{V}\leftarrow(\tilde{v}_{0};\dotso;\tilde{v}_{k-1})$.

16:$\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1})\leftarrow(\tilde{{\mathbf{m}}}_{1}^{\odot
0};\dotso;\tilde{{\mathbf{m}}}_{1}^{\odot(k-1)})$.

17:$\tilde{{\mathbf{A}}}^{\mathsf{T}}\leftarrow\tilde{\pi}_{\odot}^{-1}(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}\tilde{V}$.

18:$\tilde{{\mathbf{B}}}\leftarrow{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}(\tilde{{\mathbf{A}}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}$.

19:$\tilde{{\mathbf{B}}^{\prime}}\leftarrow{\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}}(\tilde{{\mathbf{A}}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}$.

20:For every $i\in[n]\setminus[k]$,
$\tilde{{\mathbf{m}}}_{i}\leftarrow\left(\operatorname{{\tilde{g}}}(A_{1}\cup\left\\{i\right\\}),\dotsc,\operatorname{{\tilde{g}}}(A_{k}\cup\left\\{i\right\\})\right)^{\mathsf{T}}{\left(\tilde{{\mathbf{A}}}^{\mathsf{T}}\right)}^{-1}\tilde{\pi}_{\odot}^{-1}$.

21:For every $i\in\left\\{2,\dotsc,k\right\\}$,
$\tilde{{\mathbf{m}}_{i}}\leftarrow\left(\operatorname{{\tilde{g}}}(B_{1}\cup\left\\{i\right\\}),\dotsc,\operatorname{{\tilde{g}}}(B_{k}\cup\left\\{i\right\\})\right)^{\mathsf{T}}{\left(\tilde{{\mathbf{B}}}^{\mathsf{T}}\right)}^{-1}\tilde{\pi}_{\odot}^{-1}$.

##  4 Analyzing the Algorithm

As we have seen the algorithm consists of bootstrapping steps which, as
indicated in Eqns. (7), 8, lift us “forward” from $\tilde{u}_{r-1}$ to
$\tilde{v}_{r}$ and then “back” to $\tilde{u}_{r}$. We must now control the
loss in accuracy of the statistics, in each of these steps. It turns out that
the first of these is easier and less expensive in accuracy; while the second,
in which we “invert” from approximate statistics to obtain a linear
combination of sources, is harder and also more expensive. In this section we
show how to control these steps. We rely for this control on a condition
number bound (which is not in itself algorithmic and is due entirely to the
$\zeta$-separation), which will be given in Section 5. Throughout the
analysis, we’ll assume that every multi-linear moment we use is known with
additive error bounded by

$\varepsilon\coloneqq\zeta^{C_{1}k^{2}\log k}\pi_{\min}^{C_{2}\log k}$

for constants $C_{1},C_{2}$. Choosing $C_{1}\coloneqq 60,C_{2}\coloneqq 8$ is
sufficient to give us final error $\zeta^{\Omega(k^{2}\log k)}$.

###  4.1 Bounding $\left\|\tilde{u}_{j}-u_{j}\right\|$ for $j\leq 2k$

######  Definition 4.

We define $\beta\coloneqq\frac{(\zeta/2)^{k-1}}{3k^{3}}$. We’ll frequently use
the bound $\beta\geq\zeta^{3k}$.

The following Lemma is a consequence of Theorem 22 to be proven in the next
section:

######  Lemma 5.

When the input mixture contains $3k-3$ $\zeta$-separated rows can find
disjoint sets $S,T,T^{\prime}\subseteq[n]$ of size $k-1$ each and subsets
$\mathcal{A}\subseteq 2^{S},\mathcal{B}\subseteq
2^{T},\mathcal{B}^{\prime}\subseteq 2^{T^{\prime}}$ with
$\left|\mathcal{A}\right|,\left|\mathcal{B}\right|=k$ such that the matrices
${\mathbf{A}}\coloneqq{\mathbf{M}}[\mathcal{A}]$,
${\mathbf{B}}\coloneqq{\mathbf{M}}[\mathcal{B}]$, and
${\mathbf{B}}^{\prime}\coloneqq{\mathbf{M}}[\mathcal{B}^{\prime}]$ satisfy

  1. 1.

The first row of ${\mathbf{A}}$, ${\mathbf{B}}$, and ${\mathbf{B}}^{\prime}$
is the all-ones vector, ${\mathbb{1}}$.

  2. 2.

$\sigma_{k}({\mathbf{A}}),\sigma_{k}({\mathbf{B}}),\sigma_{k}({\mathbf{B}}^{\prime})\geq\beta^{k}2^{-3k/2}k^{-3/2}$,
and

  3. 3.

$\sigma_{\max}({\mathbf{M}}[\mathcal{A}]),\sigma_{\max}({\mathbf{M}}[\mathcal{B}]),\sigma_{\max}({\mathbf{M}}[\mathcal{B}])\leq
k$.

Moreover, the first row of each of these matrices is ${\mathbb{1}}$. Finally,
the derived matrices
$\mathbf{C}_{\mathcal{B}\mathcal{A}}={\mathbf{B}}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$
and
$\mathbf{C}_{\mathcal{B}^{\prime}\mathcal{A}}={\mathbf{B}}^{\prime}\pi_{\odot}{\mathbf{A}}^{\mathsf{T}}$
satisfy

  1. 1.

$\sigma_{\max}(\mathbf{C}_{\mathcal{B}\mathcal{A}}),\sigma_{\max}(\mathbf{C}_{\mathcal{B}^{\prime}\mathcal{A}})\leq
k^{2}$, and

  2. 2.

$\sigma_{k}(\mathbf{C}_{\mathcal{B}\mathcal{A}}),\sigma_{k}(\mathbf{C}_{\mathcal{B}^{\prime}\mathcal{A}})\geq\beta^{2k}2^{-3k}k^{-3}\pi_{\min}$.

###### Proof.

This follows immediately from Theorem 22, the definition
$\mathbf{C}_{\mathcal{B}\mathcal{A}}={\mathbf{M}}[\mathcal{B}]\pi_{\odot}{\mathbf{M}}[\mathcal{A}]^{\mathsf{T}}$,
and the min-max characterization of the first and last singular values. ∎

######  Corollary 6.

$\left\|(\mathbf{C}_{\mathcal{B}\mathcal{A}})^{-1}\right\|\leq\zeta^{-10k^{2}}\pi_{\min}^{-1}$.

###### Proof.

$\left\|(\mathbf{C}_{\mathcal{B}\mathcal{A}})^{-1}\right\|\leq(\zeta^{-3k})^{2k}2^{3k}k^{3}\pi_{\min}^{-1}\leq\zeta^{-10k^{2}}\pi_{\min}^{-1}$
∎

######  Corollary 7.

$\left\|{\mathbf{A}}^{-1}\right\|,\left\|{\mathbf{B}}^{-1}\right\|,\left\|({\mathbf{B}}^{\prime})^{-1}\right\|\leq\zeta^{-6k^{2}}$.

######  Lemma 8.

$\left\|u_{i}\right\|\leq\zeta^{-6k^{2}}$, and
$\left\|v_{i}\right\|\leq\zeta^{-1}$.

###### Proof.

We observe that $\left\|u_{i}\right\|=\left\|{\mathbf{m}}_{1}^{\odot
i}{\mathbf{B}}^{-1}\right\|\leq
k\left\|{\mathbf{B}}^{-1}\right\|\leq\beta^{-k}2^{3k/2}k^{3/2}\leq\zeta^{-6k^{2}}$.
On the other hand, $\left\|v_{i}\right\|\leq k\leq\zeta^{-1}$, since $v_{i}$
is a vector of moments of products of Bernoulli random variables. ∎

######  Observation 9.

$\left\|\mathbf{C}_{\mathcal{B}\mathcal{A}}\right\|,\left\|\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}\right\|\leq
k^{3}$. If all moments are within $\varepsilon$ of their true values,
$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\right\|,\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}\right\|\leq
2k^{3}$.

######  Lemma 10.

If all multilinear moments are within $\varepsilon$ of their true values, then

$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}-\mathbf{C}_{\mathcal{B}\mathcal{A}}\right\|_{2}\leq\zeta^{-1}\varepsilon,\qquad\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}-\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}\right\|\leq\zeta^{-2}\varepsilon,$

and

$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}-\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}\right\|\leq\zeta^{-26k^{2}}\pi_{\min}^{-2}\varepsilon.$

###### Proof.

The first two inequalities just use
$\left\|\cdot\right\|_{2}\leq\left\|\cdot\right\|_{F}$. For the final
inequality we use Lemma 33:

$\displaystyle\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}-\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}\right\|\leq
2\left\|\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}\right\|^{2}\left\|\mathbf{C}_{\mathcal{B}\mathcal{A}}-{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\right\|\leq\beta^{-4k}2^{6k+1}k^{7}\pi_{\min}^{-2}\varepsilon$

∎

######  Lemma 11.

When the assumptions of Lemma 10 are satisfied, we have for any $i\in[2k]$ and
$j=\left\lceil\log i\right\rceil$,

$\left\|\tilde{v}_{i}-v_{i}\right\|,\left\|\tilde{u}_{j}-u_{j}\right\|,\left\|\tilde{u}^{\prime}_{j}-u^{\prime}_{j}\right\|\leq\zeta^{-42ik^{2}}\pi_{\min}^{-2i}\varepsilon.$

###### Proof.

Recall that we initialize the algorithm with

$\tilde{v}_{1}\leftarrow(\operatorname{{\tilde{g}}}(A_{1}\cup\left\\{1\right\\}),\
\dotsc,\operatorname{{\tilde{g}}}(A_{k}\cup\left\\{1\right\\})),\qquad\tilde{u}_{1}\leftarrow\tilde{v}_{1}({\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}})^{-1},\qquad\tilde{u}^{\prime}_{1}\leftarrow\tilde{v}_{1}({\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}})^{-1}.$

First, we observe that $\left\|\tilde{v}_{1}-v_{1}\right\|\leq\varepsilon$ by
assumption. Since $\tilde{u}_{1},\tilde{u}^{\prime}_{1}$ are computed in the
same manner here as in the loop, we’ll bound that error in the induction. Now
assume that the claim holds up to $i-1$. Recall that in each iteration of the
outer loop we compute

$\tilde{v}_{2^{i}}\leftarrow(\tilde{u}_{2^{i-1}}\otimes\tilde{u}^{\prime}_{2^{i-1}}){\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}},\qquad\tilde{u}_{2^{i}}\leftarrow\tilde{v}_{2^{i}}({\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}})^{-1},\qquad\tilde{u}^{\prime}_{2^{i}}\leftarrow\tilde{v}_{2^{i}}{\tilde{\mathbf{C}}}_{\mathcal{B}^{\prime}\mathcal{A}}.$

We’ll first focus on bounding
$\left\|\tilde{v}_{2^{i}}-v_{2^{i}}\right\|_{\infty}$. To do this we write

$\displaystyle\tilde{v}_{2^{i}}-v_{2^{i}}$
$\displaystyle=(\tilde{u}_{2^{i-1}}\otimes\tilde{u}^{\prime}_{2^{i-1}}){\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}-(u_{2^{i-1}}\otimes
u^{\prime}_{2^{i-1}})\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}$

and letting $w=\tilde{u}_{2^{i-1}}-u_{2^{i-1}}$,
$w^{\prime}=\tilde{u}^{\prime}_{2^{i-1}}-u_{2^{i-1}}$, and
$E={\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}-\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}$
we can bound the norm of the difference as follows, using the bilinearity of
the Kronecker product, Lemma 8, and the induction hypothesis:

$\displaystyle\left\|\tilde{v}_{2^{i}}-v_{2^{i}}\right\|$
$\displaystyle=\left\|(\tilde{u}_{2^{i-1}}\otimes\tilde{u}^{\prime}_{2^{i-1}}){\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}-(u_{2^{i-1}}\otimes
u^{\prime}_{2^{i-1}})\mathbf{C}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}\right\|$
$\displaystyle\leq\left\|(w\otimes\tilde{u}^{\prime}_{2^{i-1}}){\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}\right\|+\left\|(\tilde{u}_{2^{i-1}}\otimes
w^{\prime}){\tilde{\mathbf{C}}}_{\mathcal{B}+\mathcal{B}^{\prime},\mathcal{A}}\right\|+\left\|(\tilde{u}_{2^{i-1}}\otimes\tilde{u}^{\prime}_{2^{i-1}})E\right\|$
$\displaystyle\leq
2\zeta^{-42(i-1)k^{2}}\pi_{\min}^{-2(i-1)}\zeta^{-6k^{2}}k^{3}\varepsilon+\zeta^{-12k^{2}}\zeta^{-2}\varepsilon$
$\displaystyle\leq\zeta^{-(42(i-1)+16)k^{2}}\pi_{\min}^{-2(i-1)}\varepsilon$

Now we can bound $\left\|\tilde{u}_{2^{i}}-u_{2^{i}}\right\|_{\infty}$ by
observing that

$\displaystyle\tilde{u}_{2^{i}}-u_{2^{i}}$
$\displaystyle=\tilde{v}_{2^{i}}{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}-v_{2^{i}}\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}.$

Let $z=\tilde{v}_{2^{i}}-v_{2^{i}}$ and let
$D={\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}^{-1}-\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}$.
The above equation becomes

$\displaystyle\tilde{u}_{2^{i}}-u_{2^{i}}=(v_{2^{i}}+z)(\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}+D)-v_{2^{i}}\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}=v_{2^{i}}\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}+z\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}+zD$

and after taking norms and using the triangle inequality we obtain

$\displaystyle\left\|\tilde{u}_{2^{i}}-u_{2^{i}}\right\|$
$\displaystyle\leq\left\|v_{2^{i}}D\right\|+\left\|z\mathbf{C}_{\mathcal{B}\mathcal{A}}^{-1}\right\|+\left\|zD\right\|$

By Corollary 6, Lemma 8 and the induction hypothesis, we get

$\displaystyle\left\|\tilde{u}_{2^{i}}-u_{2^{i}}\right\|$
$\displaystyle\leq\zeta^{-1}\zeta^{-26k^{2}}\pi_{\min}^{-2}\varepsilon+\zeta^{-(42(i-1)+16)k^{2}}\pi_{\min}^{-2(i-1)}\varepsilon\zeta^{-16k^{2}}\pi_{\min}^{-1}+\zeta^{-(42(i-1)+16)k^{2}}\pi_{\min}^{-2(i-1)}\zeta^{-26k^{2}}\pi_{\min}^{-2}\varepsilon$
$\displaystyle\leq\zeta^{-26k^{2}-1}\pi_{\min}^{-2}\varepsilon+\zeta^{-(42(i-1)+20)k^{2}}\pi_{\min}^{-2(i-1)-1}\varepsilon+\zeta^{-(42(i-1)+42)k^{2}}\pi_{\min}^{-2(i-1)}\pi_{\min}^{-2}\varepsilon$
$\displaystyle\leq\zeta^{-42ik^{2}}\pi_{\min}^{-2i}\varepsilon$

For $j$ not a power of $2$, we can do the same analysis, and since the error
bound is increasing in $j$, the result will follow. ∎

######  Corollary 12.

Algorithm 1 will produce vectors $\tilde{v}_{i}$ for $i\leq 2k$ satisfying

$\left\|\tilde{v}_{i}-v_{i}\right\|\leq\zeta^{-42k^{2}(\log
k+1)}\pi_{\min}^{-2(\log k+1)}\varepsilon.$

###  4.2 Applying the power distribution result

######  Definition 13.

Given a mixture $\mathcal{M}$ of $k$ Bernoulli random variables with
probabilities $m_{1},\dotsc,m_{k}$ and mixing probabilities
$\pi_{1},\dotsc,\pi_{k}$, respectively, let
$[\mathcal{H}_{k+1}]_{i,j=0}^{k}=\mu_{i+j}$ be the matrix of moments of the
distribution.

######  Theorem 14 (Theorem 17 from [GMRS20]).

Given a mixture $\mathcal{M}=(m,\pi)$ as above where $m$ is $\zeta$-separated,
there is an algorithm, $\operatorname{\textsc{LearnPowerDistribution}}$, that
takes a Hankel matrix
$[\tilde{\mathcal{H}}_{k+1}]_{i,j=0}^{k}=\tilde{\mu}_{i+j}$ of approximate
moments of $\mathcal{M}$ satisfying
$\left\|\tilde{\mathcal{H}}_{k+1}-\mathcal{H}_{k+1}\right\|_{2}\leq\pi_{\min}2^{-\gamma}\zeta^{16k}$
(for some $\gamma\geq 1$) and outputs a model
$\tilde{\mathcal{M}}=(\tilde{m},\tilde{\pi})$ satisfying

$\left\|\tilde{m}-m\right\|_{\infty},\left\|\tilde{\pi}-\pi\right\|_{\infty}\leq
2^{-\gamma}$

using $O(k^{2}\log
k+k\log^{2}k\cdot\log(\log\zeta^{-1}+\log\pi_{\min}^{-1}+\gamma))$ arithmetic
operations.

######  Corollary 15.

The output $(\tilde{{\mathbf{m}}}_{1},\tilde{\pi})$ of
$\operatorname{\textsc{LearnPowerDistribution}}$ in line 14 of Algorithm 1
will satisfy

$\left\|\tilde{{\mathbf{m}}}_{1}-{\mathbf{m}}_{1}\right\|,\left\|\tilde{\pi}-\pi\right\|\leq\zeta^{-42k^{2}(\log
k+1)-16k-1}\pi_{\min}^{-2\log k-3}\varepsilon.$

This step will use $O(k^{2}\log k+k\log^{2}k\cdot\log\log(\varepsilon^{-1}))$
arithmetic operations.

###### Proof.

Every entry $(\tilde{v}_{i})_{1}$ satisfies
$\left\|(\tilde{v}_{i})_{1}-(v_{i})_{1}\right\|\leq\zeta^{-42k^{2}(\log
k+1)}\pi_{\min}^{-2(\log k+1)}\varepsilon$ so

$\left\|\tilde{\mathcal{H}}_{k+1}-\mathcal{H}_{k+1}\right\|\leq\zeta^{-42k^{2}(\log
k+1)}\pi_{\min}^{-2(\log k+1)}$

which implies that

$\left\|\tilde{{\mathbf{m}}}_{1}-{\mathbf{m}}_{1}\right\|_{\infty},\left\|\tilde{\pi}-\pi\right\|_{\infty}\leq\zeta^{-42k^{2}(\log
k+1)-16k}\pi_{\min}^{-2(\log k+1)-1}\varepsilon.$

Finally, we add a factor of $\zeta^{-1}$ to convert back to the Euclidean norm
to get the stated bound. ∎

###  4.3 Solving for the rest of the model

Once we’ve computed $\tilde{{\mathbf{m}}}_{1}$ and $\tilde{\pi}$, we’ll use
them to compute the remaining model parameters. In this section we bound the
additional error introduced by these computations.

######  Observation 16.

$\left\|\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1})-\operatorname{Vdm}({\mathbf{m}}_{1})\right\|\leq\zeta^{-1}\left\|\tilde{{\mathbf{m}}}_{1}-{\mathbf{m}}_{1}\right\|$.

######  Observation 17 (Claim 26 in [GMRS20]).

$\left\|\operatorname{Vdm}({\mathbf{m}}_{1})^{-1}\right\|\leq
2^{k}/\zeta^{k-1}\leq\zeta^{-2k}$ when ${\mathbf{m}}_{1}$ is
$\zeta$-separated.

######  Lemma 18.

The computed $\tilde{{\mathbf{A}}}$ produced by Algorithm 1 will satisfy
$\left\|\tilde{{\mathbf{A}}}-{\mathbf{A}}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-20k-6}\pi_{\min}^{-2\log k-5}\varepsilon$.

###### Proof.

First, we observe that
$\left\|\tilde{V}\right\|\leq\left\|V\right\|\leq\zeta^{-2}$ and
$\left\|\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1})^{-1}\right\|\leq\zeta^{-3k}$
by Lemma 33 and Observation 17. Now
$\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-16k-2}\pi_{\min}^{-2\log k-5}\varepsilon$ by Lemma 33. Thus,

$\displaystyle\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}\right\|\left\|\tilde{V}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-18k-5}\pi_{\min}^{-2\log k-5}\varepsilon.$

Now
$\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}-(\operatorname{Vdm}({\mathbf{m}}_{1}))^{-1}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-20k-2}\pi_{\min}^{-2\log k-3}\varepsilon$ by Lemma 33, so

$\left\|\tilde{\pi}_{\odot}^{-1}\right\|\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}-(\operatorname{Vdm}({\mathbf{m}}_{1}))^{-1}\right\|\left\|\tilde{V}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-20k-4}\pi_{\min}^{-2\log k-4}\varepsilon.$

Finally, $\left\|\tilde{V}-V\right\|\leq\zeta^{-42k^{2}(\log
k+1)-1}\pi_{\min}^{-2\log k+1}\varepsilon$ so that

$\left\|\pi_{\odot}^{-1}\right\|\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}\right\|\left\|\tilde{V}-V\right\|\leq\zeta^{-42k^{2}(\log
k+1)-3k-1}\pi_{\min}^{-2\log k-3}\varepsilon.$

Putting these together, we easily obtain

$\displaystyle\left\|\tilde{{\mathbf{A}}}-{\mathbf{A}}\right\|$
$\displaystyle=\left\|\tilde{\pi}_{\odot}^{-1}(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}\tilde{V}-\pi_{\odot}^{-1}(\operatorname{Vdm}({\mathbf{m}}_{1}))^{-1}V\right\|$
$\displaystyle\leq\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}\right\|\left\|\tilde{V}\right\|+\left\|\tilde{\pi}_{\odot}^{-1}\right\|\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}-(\operatorname{Vdm}({\mathbf{m}}_{1}))^{-1}\right\|\left\|\tilde{V}\right\|$
$\displaystyle\qquad+\left\|\pi_{\odot}^{-1}\right\|\left\|(\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1}))^{-1}\right\|\left\|\tilde{V}-V\right\|(\left\|\operatorname{Vdm}({\mathbf{m}}_{1})^{-1}E_{2}\right\|+\left\|E_{1}V\right\|_{\infty})+\left\|w\right\|\left\|\operatorname{Vdm}(\tilde{{\mathbf{m}}}_{1})^{-1}\right\|\left\|\tilde{V}\right\|$
$\displaystyle\leq\zeta^{-42k^{2}(\log k+1)-20k-6}\pi_{\min}^{-2\log
k-5}\varepsilon.\qed$

######  Lemma 19.

The matrices $\tilde{{\mathbf{A}}},\tilde{{\mathbf{B}}},$ and
$\tilde{{\mathbf{B}}}^{\prime}$ satisfy

$\left\|\tilde{{\mathbf{A}}}^{-1}\right\|,\left\|\tilde{{\mathbf{B}}}^{-1}\right\|,\left\|(\tilde{{\mathbf{B}}}^{\prime})^{-1}\right\|\leq\zeta^{-7k^{2}}.$

###### Proof.

By Lemma 33 and Corollary 7. ∎

######  Lemma 20.

The matrices $\tilde{{\mathbf{B}}}^{\prime}$ and $\tilde{{\mathbf{B}}}$
produced by Algorithm 1 will satisfy

$\left\|\tilde{{\mathbf{B}}}^{\prime}-{\mathbf{B}}^{\prime}\right\|,\left\|\tilde{{\mathbf{B}}}-{\mathbf{B}}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-14k^{2}-20k-16}\pi_{\min}^{-2\log k-6}\varepsilon.$

###### Proof.

We’ll prove the claim for $\tilde{{\mathbf{B}}}$; the proof is identical for
$\tilde{{\mathbf{B}}}^{\prime}$. We can bound
$\tilde{{\mathbf{B}}}-{\mathbf{B}}$ using the same tools as in the previous
bounds. First, we bound
$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}-\mathbf{C}_{\mathcal{B}\mathcal{A}}\right\|\left\|\tilde{{\mathbf{A}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}\right\|\leq\zeta^{-7k^{2}-2}\pi_{\min}^{-1}\varepsilon$.
Now

$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\right\|\left\|\tilde{{\mathbf{A}}}^{-1}-{\mathbf{A}}^{-1}\right\|\left\|\tilde{\pi}^{-1}_{\odot}\right\|\leq\zeta^{-4}\zeta^{-42k^{2}(\log
k+1)-14k^{2}-20k-11}\pi_{\min}^{-2\log k-6}\varepsilon$

and

$\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\right\|\left\|\tilde{{\mathbf{A}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|\zeta^{-4}\zeta^{-7k^{2}}\zeta^{-42k^{2}(\log
k+1)-7k^{2}-16k-6}\pi_{\min}^{-2\log k-5}\varepsilon.$

The resulting bound is

$\displaystyle\left\|\tilde{{\mathbf{B}}}-{\mathbf{B}}\right\|$
$\displaystyle=\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}(\tilde{{\mathbf{A}}}^{\mathsf{T}})^{-1}\pi_{\odot}^{-1}-\mathbf{C}_{\mathcal{B}\mathcal{A}}({\mathbf{A}}^{\mathsf{T}})^{-1}\pi_{\odot}\right\|$
$\displaystyle\leq\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}-\mathbf{C}_{\mathcal{B}\mathcal{A}}\right\|\left\|\tilde{{\mathbf{A}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}\right\|+\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\right\|\left\|\tilde{{\mathbf{A}}}^{-1}-{\mathbf{A}}^{-1}\right\|\left\|\tilde{\pi}^{-1}_{\odot}\right\|+\left\|{\tilde{\mathbf{C}}}_{\mathcal{B}\mathcal{A}}\right\|\left\|\tilde{{\mathbf{A}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|$
$\displaystyle\leq\zeta^{-42k^{2}(\log k+1)-14k^{2}-20k-16}\pi_{\min}^{-2\log
k-6}\varepsilon.$

∎

######  Lemma 21.

Algorithm 1 will compute $\tilde{{\mathbf{m}}}_{i}$ satisfying
$\left\|\tilde{{\mathbf{m}}}_{i}-{\mathbf{m}}_{i}\right\|_{\infty}\leq\zeta^{-42k^{2}(\log
k+1)-14k^{2}-20k-19}\pi_{\min}^{-2\log k-8}\varepsilon$ for all $i$.

###### Proof.

We’ll compute the bound using the inversion of $\tilde{{\mathbf{B}}}$ since
this will give us the worst case. Let
$\tilde{y}=(\operatorname{{\tilde{g}}}(B_{1}\cup\left\\{i\right\\}),\dotsc,\operatorname{{\tilde{g}}}(B_{k}\cup\left\\{1\right\\})$
and let
$y=(\operatorname{g}(B_{1}\cup\left\\{i\right\\}),\dotsc,\operatorname{g}(B_{k}\cup\left\\{1\right\\})$.
We note that $\left\|\tilde{y}-y\right\|\leq\zeta^{-1}\zeta^{-1}\varepsilon$
by assumption, and $\left\|\tilde{y}\right\|\leq\zeta^{-2}$. Then

$\left\|\tilde{y}-y\right\|\left\|\tilde{{\mathbf{B}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}\right\|\leq\zeta^{-7k^{2}-2}\pi_{\min}^{-2}\varepsilon,$
$\left\|\tilde{y}\right\|\left\|\tilde{{\mathbf{B}}}^{-1}-{\mathbf{B}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-14k^{2}-20k-19}\pi_{\min}^{-2\log k-8}\varepsilon,$

and

$\left\|\tilde{y}\right\|\left\|\tilde{{\mathbf{B}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|\leq\zeta^{-42k^{2}(\log
k+1)-16k-7k^{2}-3}\pi_{\min}^{-2\log k-5}\varepsilon^{2}.$

Of the three terms, the middle one clearly dominates so that we get

$\displaystyle\left\|\tilde{{\mathbf{m}}}_{i}-{\mathbf{m}}_{i}\right\|$
$\displaystyle=\left\|\tilde{y}\tilde{{\mathbf{B}}}^{-1}\tilde{\pi}_{\odot}^{-1}-y{\mathbf{B}}^{-1}\pi_{\odot}^{-1}\right\|$
$\displaystyle\leq\left\|\tilde{y}-y\right\|\left\|\tilde{{\mathbf{B}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}\right\|+\left\|\tilde{y}\right\|\left\|\tilde{{\mathbf{B}}}^{-1}-{\mathbf{B}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}\right\|+\left\|\tilde{y}\right\|\left\|\tilde{{\mathbf{B}}}^{-1}\right\|\left\|\tilde{\pi}_{\odot}^{-1}-\pi_{\odot}^{-1}\right\|$
$\displaystyle\leq\zeta^{-42k^{2}(\log k+1)-14k^{2}-20k-19}\pi_{\min}^{-2\log
k-8}\varepsilon$

∎

##  5 The core stability bound

Recall we have set
$\beta\coloneqq(\zeta/2)^{k-1}/3k^{2}\geq\zeta^{2k}2^{-k}\geq\zeta^{3k}$.

######  Theorem 22.

Let $S$ be a set of $k-1$ $\zeta$-separated vectors
${\mathbf{m}}_{1},\ldots,{\mathbf{m}}_{k-1}$. Then there exists a set
$J\subseteq 2^{S}$, $|J|=k$, such that
$\sigma_{k}({\mathbf{M}}[J])\geq\beta^{k}2^{-3k/2}k^{-3/2}$ and
$\sigma_{\max}({\mathbf{M}}[J])\leq k$. The first row of ${\mathbf{M}}[J]$ is
${\mathbb{1}}$ (corresponding to $\emptyset\in J$).

We start with some definitions.

######  Definition 23.

For any subspace $U\subseteq\mathbb{R}^{n}$, let ${\mathcal{P}_{U}}$ denote
the orthogonal projection onto $U$ and let ${\mathcal{P}^{\perp}_{U}}$ denote
the orthogonal projection onto the orthogonal complement subspace. Then
$I={\mathcal{P}_{U}}+{\mathcal{P}^{\perp}_{U}}$.

In the following, $\left\|\cdot\right\|$ will always denote the $2\to 2$
operator norm of a matrix.

######  Definition 24 ($U$-operator norm).

Given a subspace $U\subseteq\mathbb{R}^{n}$, we can define the $U$-operator
norm of a matrix $C\in\mathbb{R}^{m\times n}$, denoted $\left\|C\right\|_{U}$,
as follows:

$\left\|C\right\|_{U}\coloneqq\max_{0\neq u\in
U}\left\|Cu\right\|/\left\|u\right\|.$

##### Bounding the condition number of ${\mathbf{M}}[2^{S}]$.

Most of our work will go into lower bounding the $k$’th singular value of
${\mathbf{M}}[2^{S}]$. This is where the $\zeta$-separation condition is
essential. (Even the non-quantitative result that ${\mathbf{M}}[2^{S}]$ has
full column rank is not trivial; indeed, the result of [TMMA18] that sources
can be identified from $2k-1$ separated bits, is implied by the non-
quantitative version of this section, omitted here, which employs the same
approach but is considerably shorter.)

######  Lemma 25.

Let $U$ be any subspace of $\mathbb{R}^{n}$ and let $C\in\mathbb{R}^{n\times
n}$. For all $n\geq 1$, $\left\|{\mathcal{P}^{\perp}_{U}}C^{n}\right\|_{U}\leq
n\|C\|^{n-1}\|{\mathcal{P}^{\perp}_{U}}C\|_{U}$.

###### Proof.

We need to show that $\forall n\geq 2\;\max_{0\neq u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}C^{n}u\right\|/\left\|u\right\|\leq
n\left\|C\right\|^{n-1}\max_{0\neq u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}Cu\right\|/\left\|u\right\|$. Using that
$\left\|{\mathcal{P}_{U}}\right\|,\left\|{\mathcal{P}^{\perp}_{U}}\right\|\leq
1$, we have

$\left\|{\mathcal{P}^{\perp}_{U}}C^{n}u\right\|=\left\|{\mathcal{P}^{\perp}_{U}}C({\mathcal{P}_{U}}+{\mathcal{P}^{\perp}_{U}})C^{n-1}u\right\|\leq\left\|{\mathcal{P}^{\perp}_{U}}C{\mathcal{P}_{U}}C^{n-1}u\right\|+\left\|{\mathcal{P}^{\perp}_{U}}C{\mathcal{P}^{\perp}_{U}}C^{n-1}u\right\|.$

So (using $\left\|C^{n-1}\right\|\leq\left\|C\right\|^{n-1}$ and
$\left\|{\mathcal{P}_{U}}\right\|\leq 1$ in the first term, and
$\left\|{\mathcal{P}^{\perp}_{U}}\right\|\leq 1$ in the second term)

$\max_{0\neq u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}C^{n}u\right\|/\left\|u\right\|\leq\left\|C\right\|^{n-1}\max_{0\neq
u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}Cu\right\|/\left\|u\right\|+\left\|C\right\|\max_{0\neq
u\in U}\left\|{\mathcal{P}^{\perp}_{U}}C^{n-1}u\right\|/\left\|u\right\|$

and applying induction this is

$\leq(\left\|C\right\|^{n-1}+(n-1)\left\|C\right\|^{n-1})\max_{0\neq u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}Cu\right\|/\left\|u\right\|.\qed$

Recall from Section 3 that the coefficient norms of the interpolation
polynomials are
$\left\|p_{v,i}\right\|_{c}=\sum_{0}^{k-1}j\left|p_{v,i,j}\right|$.

######  Lemma 26.

$\left\|{\mathcal{P}^{\perp}_{U}}P_{(i)}\right\|_{U}\leq\left\|p_{v,i}\right\|_{c}\cdot\left\|{\mathcal{P}^{\perp}_{U}}v_{\odot}\right\|_{U}$.

###### Proof.

We are to show that $\max_{0\neq u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}P_{(i)}u\right\|/\left\|u\right\|\leq\left\|p_{v,i}\right\|_{c}\max_{0\neq
u\in U}\left\|{\mathcal{P}^{\perp}_{U}}v_{\odot}u\right\|/\left\|u\right\|$.
We have

$\max_{u}\left\|{\mathcal{P}^{\perp}_{U}}P_{(i)}u\right\|/\left\|u\right\|=\max_{u}\left\|{\mathcal{P}^{\perp}_{U}}p_{v,i}(v_{\odot})u\right\|/\left\|u\right\|\leq\sum_{j}\max_{u}\left|p_{v,i,j}\right|\cdot\left\|{\mathcal{P}^{\perp}_{U}}(v_{\odot})^{j}u\right\|/\left\|u\right\|.$

Note that $\left\|v_{\odot}\right\|\leq 1$. So applying Lemma 25:

$\max_{0\neq u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}P_{(i)}u\right\|/\left\|u\right\|\leq\sum_{j}\left|p_{v,i,j}\right|j\max_{0\neq
u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}v_{\odot}u\right\|/\left\|u\right\|=\left\|p_{v,i}\right\|_{c}\max_{0\neq
u\in
U}\left\|{\mathcal{P}^{\perp}_{U}}v_{\odot}u\right\|/\left\|u\right\|.\qed$

######  Lemma 27.

Let ${\mathbb{1}}\in U\subsetneq\mathbb R^{k}$ (strict containment). Then
there is an $i$ s.t.
$\|{\mathcal{P}^{\perp}_{U}}P_{(i)}{\mathbb{1}}/\sqrt{k}\|>1/3k^{2}$.

###### Proof.

Let ${\varepsilon}=1/3k^{2}$ and suppose to the contrary that
$\|{\mathcal{P}^{\perp}_{U}}P_{(i)}{\mathbb{1}}/\sqrt{k}\|\leq{\varepsilon}$
for all $i$.

Fix some $u\in U$ and let $v_{i}=P_{(i)}u$ and
$u_{i}={\mathcal{P}_{U}}v_{i}={\mathcal{P}_{U}}P_{(i)}u$. Then since $\sum
P_{(i)}=I$, $\sum u_{i}={\mathcal{P}_{U}}\sum P_{(i)}u=u$.

Taking $u={\mathbb{1}}/\sqrt{k}\in U$, we find $\left\|v_{i}\right\|^{2}\geq
1/k$ for all $1\leq i\leq k$. (On the other hand of course
$\left\|v_{i}\right\|\leq 1$.) By assumption
${\varepsilon}\geq\left\|{\mathcal{P}^{\perp}_{U}}P_{(i)}{\mathbb{1}}/\sqrt{k}\right\|=\left\|{\mathcal{P}^{\perp}_{U}}v_{i}\right\|=\left\|{\mathcal{P}_{U}}v_{i}-v_{i}\right\|$.
So $\left\|{\mathcal{P}_{U}}v_{i}\right\|\geq 1/\sqrt{k}-{\varepsilon}$. For
$i\neq j$, $v_{i}^{*}v_{j}=0$, so
$\left|({\mathcal{P}_{U}}v_{i})^{*}({\mathcal{P}_{U}}v_{j})\right|=\left|(({\mathcal{P}_{U}}v_{i}-v_{i})+v_{i})^{*}(({\mathcal{P}_{U}}v_{j}-v_{j})+v_{j})\right|\leq\left|v_{i}^{*}v_{j}\right|+2{\varepsilon}+{\varepsilon}^{2}=2{\varepsilon}+{\varepsilon}^{2}$.
So we have $k$ vectors ${\mathcal{P}_{U}}v_{i}$ in $U$, each of length $\geq
1/\sqrt{k}-{\varepsilon}$, with inner products $\leq
2{\varepsilon}+{\varepsilon}^{2}$ in absolute value. Their Gram matrix has
diagonal entries $\geq(1/\sqrt{k}-{\varepsilon})^{2}\geq
1/k-2{\varepsilon}+{\varepsilon}^{2}$ and off-diagonal entries $\leq
2{\varepsilon}+{\varepsilon}^{2}$ in absolute value. Since
${\varepsilon}=1/3k^{2}$, this matrix is strictly diagonally dominant and
therefore nonsingular, contradicting the assumption that $U$ has dimension
$<k$. ∎

######  Lemma 28.

Let ${\mathbb{1}}\in U\subset\mathbb R^{k}$. Then
$\|{\mathcal{P}^{\perp}_{U}}v_{\odot}\|_{U}>\beta$.

###### Proof.

By Lemma 27, there is an $i$ s.t.
$\left\|{\mathcal{P}^{\perp}_{U}}P_{(i)}\right\|_{U}>1/3k^{2}$. Applying Lemma
26 to this $i$, we have that
$\left\|{\mathcal{P}^{\perp}_{U}}v_{\odot}\right\|_{U}\geq\frac{1}{3k^{2}\left\|p_{v,i}\right\|_{c}}$.
Now we need an upper bound on $\left\|p_{v,i}\right\|_{c}$. Recall
$p_{v,i}(x)=\frac{\prod_{j\neq i}(x-\lambda_{j})}{\prod_{j\neq
i}(\lambda_{i}-\lambda_{j})}$. Simply by upper-bounding all $\lambda_{j}$ by
$1$ and lower-bounding all separations by $\zeta$, we have the bound
$\left\|p_{v,i}\right\|_{c}\leq\zeta^{1-k}\sum_{\ell^{\prime}=0}^{k-1}\ell^{\prime}\binom{k-1}{\ell^{\prime}}\leq(k-1)(2/\zeta)^{k-1}$.
∎

######  Lemma 29.

Fix $\zeta$-separated vectors ${\mathbf{m}}_{1},\dotsc,{\mathbf{m}}_{k-1}$.
Then there exists a $k\times k$ matrix $V$, with rows $v_{1},\dotsc,v_{k}$,
such that:

  1. 1.

The row $v_{1}={\mathbb{1}}/\sqrt{k}$, and for $\ell\geq 1$,
$v_{\ell+1}\coloneqq{\mathbf{m}}_{\ell}\odot u$, where $u$ is a unit vector in
$U_{\ell}\coloneqq\operatorname{span}\left\\{v_{1},\dotsc,v_{\ell}\right\\}$.

  2. 2.

Any unit vector in $U_{\ell}$ can be formed as a linear combination of the
rows in $[{\mathbf{m}}_{R}]_{R\subseteq[\ell]}$ with coefficients bounded in
maximum magnitude by $\beta^{-\ell}$.

###### Proof.

We’ll use induction over $\ell$. For $\ell>1$, we’ll apply Lemma 28 to find a
unit vector $u^{*}\in U_{\ell-1}$ such that
$\left\|{\mathcal{P}^{\perp}_{U}}({\mathbf{m}}_{\ell}\odot
u^{*})\right\|>\beta$. We’ll set $v_{\ell}\coloneqq{\mathbf{m}}_{\ell}\odot
u^{*}$. Then
$U_{\ell}=\operatorname{span}(U_{\ell-1},{\mathcal{P}^{\perp}_{U}}({\mathbf{m}}_{\ell}\odot
u^{*}))$. Now $\left\|{\mathbf{m}}_{\ell}\odot u^{*}\right\|\leq 1$. Any
vector $u\in U_{\ell}$ can be written as
$c_{1}u^{(\ell-1)}+c_{2}({\mathbf{m}}_{\ell}\odot u^{*})$ for $u^{(\ell-1)}\in
U_{\ell-1}$ a unit vector and where $\left|c_{2}\right|\leq\beta^{-1}$ and
$\left|c_{1}\right|\leq\beta^{-1}$. (This is just the operator norm of the
inverse of $\begin{bmatrix}1&0\\\ \sqrt{1-\beta^{2}}&\beta\end{bmatrix}$.) By
the induction hypothesis, we can expand $u^{(\ell-1)}$ and $u^{*}$ in terms of
${{\mathbf{m}}_{S}}_{R\subseteq[\ell-1]}$ as follows:

$\displaystyle u$
$\displaystyle=c_{1}u^{(\ell-1)}+c_{2}({\mathbf{m}}_{\ell}\odot u^{*})$
$\displaystyle=c_{1}\sum_{R\subseteq[\ell-1]}\alpha_{R}{\mathbf{m}}_{R}+c_{2}{\mathbf{m}}_{\ell}\sum_{R\subseteq[\ell-1]}\alpha^{\prime}_{R}{\mathbf{m}}_{R}$
$\displaystyle=\sum_{R\subseteq[\ell-1]}\left(c_{1}\alpha_{R}{\mathbf{m}}_{R}+c_{2}\alpha^{\prime}_{R}{\mathbf{m}}_{R\cup\left\\{\ell\right\\}}\right)$

where
$\left\|\alpha\right\|_{\infty},\left\|\alpha^{\prime}\right\|_{\infty}\leq\beta^{\ell-1}$.
The claim follows immediately. ∎

######  Corollary 30.

For any unit column vector $z\in\mathbb{R}^{k}$, the matrix
${\mathbf{M}}[2^{S}]$ satisfies
$\left\|{\mathbf{M}}[2^{S}]z\right\|_{\infty}\geq\beta^{k}2^{-k}$.

###### Proof.

By Lemma 29, we know that we can write
$z^{\mathsf{T}}=\lambda^{\mathsf{T}}{\mathbf{M}}$ where
$\lambda\in\mathbb{R}^{2^{k}}$ and
$\left\|\lambda\right\|_{\infty}\leq\beta^{-k}$. Thus,
$1=\left\|z\right\|^{2}=\sum_{R\subseteq S}\lambda_{R}{\mathbf{m}}_{R}z$.
There must be some $R\subseteq S$ for which
$\left|\lambda_{R}{\mathbf{m}}_{R}z\right|\geq 1/2^{k}$. Since
$\left|\lambda_{R}\right|\leq\beta^{-k}$ we immediately get that
$\left\|{\mathbf{M}}z\right\|_{\infty}\geq\left|{\mathbf{m}}_{R}z\right|\geq\beta^{k}2^{-k}$.
∎

For a subset $T$ with $\left|T\right|>k-1$, the bound on the largest singular
value of ${\mathbf{M}}[2^{T}]$ increases to $k2^{\left|T\right|}$ while the
lower bound remains unchanged.

######  Corollary 31.

$\sigma_{\max}({\mathbf{M}}[2^{S}])\leq k2^{k-1}$ and
$\sigma_{k}({\mathbf{M}}[2^{S}])\geq\beta^{k}2^{-k}/k$.

###### Proof.

The largest singular value of ${\mathbf{M}}[2^{S}]$ is easily upper bounded by
$k2^{k-1}$. Lemma 29 gives the bound on $\sigma_{k}$. ∎

##### Bounding the condition number of a $k\times k$ submatrix of
${\mathbf{M}}[2^{S}]$.

We can now use the following result from [FOS08] to find a $k\times k$
submatrix that is similarly well-conditioned.

######  Lemma 32 (Corollary 6 in [FOS08]).

Let $A\in\mathbb{R}^{k\times n}$ with $k<n$, and let
$\sigma_{k}(A)\geq\varepsilon$. Then there exists a subset of the columns
$J\subseteq[n]$ with $\left|J\right|=k$ such that
$\sigma_{k}(A_{J})\geq\varepsilon/\sqrt{k(n-k)+1}$.

Proof of Theorem 22. The upper bound is trivial since all entries are in
$[0,1]$. The lower bound follows by applying Lemma 32 to Corollary 31. ∎

##  Appendix A Miscellaneous Proofs

######  Lemma 33.

For an invertible $n\times n$ matrix $M$ and a perturbed matrix $\tilde{M}$,
if $\left\|\tilde{M}-M\right\|=\varepsilon\leq\sigma_{n}(M)/2$, then

$\left\|\tilde{M}^{-1}-M^{-1}\right\|\leq
2\left\|M^{-1}\right\|^{2}\varepsilon,\quad\text{ and
}\left\|\tilde{M}^{-1}\right\|\leq 2\left\|M^{-1}\right\|.$

###### Proof.

First, we observe that

$\left\|\tilde{M}^{-1}\right\|=\frac{1}{\sigma_{n}(\tilde{M})}\leq\frac{1}{\sigma_{n}(M)-\sigma_{n}(M)/2}\leq
2\left\|M^{-1}\right\|.$

We use the identity
$\tilde{M}^{-1}-M^{-1}=\tilde{M}^{-1}\left(M-\tilde{M}\right)M^{-1}$.

$\left\|\tilde{M}^{-1}-M^{-1}\right\|=\left\|\tilde{M}^{-1}\left(M-\tilde{M}\right)M^{-1}\right\|\leq
2\left\|M^{-1}\right\|^{2}\left\|M-\tilde{M}\right\|.$

∎

## References

  * [ADXR20] B. Aragam, C. Dan, E. P. Xing, and P. Ravikumar.  Identifiability of nonparametric mixture models and Bayes optimal clustering.  Ann. Statist., 48(4):2277–2302, 2020.
  * [AFH+12] A. Anandkumar, D. P. Foster, D. J. Hsu, S. M. Kakade, and Y.-K. Liu.  A spectral algorithm for latent Dirichlet allocation.  In F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger, editors, Advances in Neural Information Processing Systems 25, pages 917–925. Curran Associates, Inc., 2012.
  * [AGM12] S. Arora, R. Ge, and A. Moitra.  Learning topic models — going beyond SVD.  In Proc. 53rd Ann. IEEE Symp. on Foundations of Computer Science, 2012.
  * [AHK12] A. Anandkumar, D. Hsu, and S. M. Kakade.  A method of moments for mixture models and hidden Markov models.  In Proc. 25th Ann. Conf. on Computational Learning Theory, pages 33.1–33.34, 2012.
  * [AMR09] E. Allman, C. Matias, and J. Rhodes.  Identifiability of parameters in latent structure models with many observed variables.  Annals of Statistics, 37:3099–3132, 2009.
  * [Bli64] W. R. Blischke.  Estimating the parameters of mixtures of binomial distributions.  Journal of the American Statistical Association, 59(306):510–528, 1964.
  * [CGG01] M. Cryan, L. Goldberg, and P. Goldberg.  Evolutionary trees can be learned in polynomial time in the two state general Markov model.  SIAM J. Comput., 31(2):375–397, 2001.  Prev. FOCS ’98.
  * [CM19] S. Chen and A. Moitra.  Beyond the low-degree algorithm: mixtures of subcubes and their applications.  In Proc. 51st Ann. ACM Symp. on Theory of Computing, pages 869–880, 2019.
  * [CPnR00] M. A. Carreira-Perpiñán and S. Renals.  Practical identifiability of finite mixtures of multivariate Bernoulli distributions.  Neural Computation, 12(1):141–152, 2000.
  * [CR08] K. Chaudhuri and S. Rao.  Learning mixtures of product distributions using correlations and independence.  In Proc. 21st Ann. Conf. on Computational Learning Theory, pages 9–20, 2008.
  * [Das99] S. Dasgupta.  Learning mixtures of Gaussians.  In Proc. 40th Ann. IEEE Symp. on Foundations of Computer Science, page 634–644, 1999.
  * [EH81] B. S. Everitt and D. J. Hand.  Mixtures of discrete distributions, pages 89–105.  Springer Netherlands, Dordrecht, 1981.
  * [FM99] Y. Freund and Y. Mansour.  Estimating a mixture of two product distributions.  In Proc. 12th Ann. Conf. on Computational Learning Theory, pages 183–192, July 1999.
  * [FOS08] J. Feldman, R. O’Donnell, and R. A. Servedio.  Learning mixtures of product distributions over discrete domains.  SIAM J. Comput., 37(5):1536–1564, 2008.
  * [GMRS20] S. L. Gordon, B. Mazaheri, Y. Rabani, and L. J. Schulman.  The sparse Hausdorff moment problem, with application to topic models.  ArXiv:2007.08101, 2020.
  * [JGHV04] A. Juan, J. García-Hernández, and E. Vidal.  EM initialisation for Bernoulli mixture learning.  In A. Fred, T. M. Caelli, R. P. W. Duin, A. C. Campilho, and D. de Ridder, editors, Structural, Syntactic, and Statistical Pattern Recognition, pages 635–643, Berlin, Heidelberg, 2004. Springer.
  * [JV02] A. Juan and E. Vidal.  On the use of Bernoulli mixture models for text classification.  Pattern Recognition, 35(12):2705–2710, 2002.
  * [JV04] A. Juan and E. Vidal.  Bernoulli mixture models for binary images.  In Proc. of the 17th International Conference on Pattern Recognition, volume 3, pages 367–370, 2004.
  * [JWL+05] Y. Ji, C. Wu, P. Liu, J. Wang, and K. R. Coombes.  Applications of beta-mixture models in bioinformatics.  Bioinformatics, 21(9):2118–2122, 2005.
  * [KF09] D. Koller and N. Friedman.  Probabilistic Graphical Models: Principles and Techniques.  MIT Press, 2009.
  * [KKM+19] Y. Kim, F. Koehler, A. Moitra, E. Mossel, and G. Ramnarayan.  How many subpopulations is too many? Exponential lower bounds for inferring population histories.  In L. Cowen, editor, Int’l Conf. on Research in Computational Molecular Biology, volume 11457 of Lecture Notes in Computer Science, pages 136–157. Springer, 2019.
  * [KMR+94] M. Kearns, Y. Mansour, D. Ron, R. Rubinfeld, R. Schapire, and L. Sellie.  On the learnability of discrete distributions.  In Proc. 26th Ann. ACM Symp. on Theory of Computing, pages 273–282, 1994.
  * [Koo50] T. C. Koopmans.  Statistical Inference in Dynamic Economic Models.  John Wiley and Sons, Inc., 1950.
  * [KR50] T. C. Koopmans and O. Reiersol.  The identification of structural characteristics.  Ann. Math. Statist., 21(2):165–181, 1950.
  * [Lin95] B. G. Lindsay.  Mixture models: theory, geometry and applications.  In NSF-CBMS regional conference series in probability and statistics, pages i–163. JSTOR, 1995.
  * [LRSS15] J. Li, Y. Rabani, L. J. Schulman, and C. Swamy.  Learning arbitrary statistical mixtures of discrete distributions.  In Proc. 47th Ann. ACM Symp. on Theory of Computing, pages 743–752, 2015.
  * [LWPA16] C. Li, B. Wang, V. Pavlu, and J. Aslam.  Conditional Bernoulli mixtures for multi-label classification.  In Proc. of the 33rd International Conference on Machine Learning, pages 2482–2491, 2016.
  * [MLR19] G. J. McLachlan, S. X. Lee, and S. I. Rathnayake.  Finite mixture models.  Annual Review of Statistics and Its Application, 6(1):355–378, 2019\.
  * [New86] S. Newcomb.  A generalized theory of the combination of observations so as to obtain the best result.  American Journal of Mathematics, 8(4):343–366, 1886.
  * [NMR20] A. Najafi, S. A. Motahari, and H. R. Rabiee.  Reliable clustering of Bernoulli mixture models.  Bernoulli, 26(2):1535–1559, 2020.
  * [Pea94] K. Pearson.  Contributions to the mathematical theory of evolution III.  Philosophical Transactions of the Royal Society of London (A.), 185:71–110, 1894.
  * [Pea09] J. Pearl.  Causality.  Cambridge, 2nd edition, 2009.
  * [PKM16] J. A. Palmer, K. Kreutz-Delgado, and S. Makeig.  A generalized multivariate logistic model and EM algorithm based on the normal variance mean mixture representation.  In IEEE Statistical Signal Processing Workshop, pages 1–5, 2016\.
  * [PSD00] J. K. Pritchard, M. Stephens, and P. Donnelly.  Inference of population structure using multilocus genotype data.  Genetics, 155(2):945–959, 2000.
  * [RSS14] Y. Rabani, L. J. Schulman, and C. Swamy.  Learning mixtures of arbitrary distributions over large discrete domains.  In Proc. 5th Conf. on Innovations in Theoretical Computer Science, pages 207–224, 2014.
  * [RVS20] A. Ritchie, R. A. Vandermeulen, and C. D. Scott.  Consistent estimation of identifiable nonparametric mixture models from grouped observations.  CoRR, abs/2006.07459, 2020.
  * [Sch17] K. Schmüdgen.  The Moment Problem, volume 277 of Graduate Texts in Mathematics.  Springer International Publishing, 2017.
  * [SGS00] P. Spirtes, C. Glymour, and R. Scheines.  Causation, Prediction and Search.  MIT Press, second edition, 2000.
  * [Sim15] B. Simon.  A comprehensive course in analysis.  American Mathematical Society, 2015.
  * [Tei63] H. Teicher.  Identifiability of finite mixtures.  Ann. Math. Statist., 34(4):1265–1269, 12 1963.
  * [TMMA18] B. Tahmasebi, S. A. Motahari, and M. A. Maddah-Ali.  On the identifiability of finite mixtures of finite product measures.  IEEE International Symposium on Information Theory (ISIT) 2018 and arXiv:1807.05444v1, 2018.
  * [TSM85] D. M. Titterington, A. F. M. Smith, and U. E. Makov.  Statistical Analysis of Finite Mixture Distributions.  John Wiley and Sons, Inc., 1985.
  * [VS15] R. A. Vandermeulen and C. Scott.  On the identifiability of mixture models from grouped samples.  ArXiv, abs/1502.06644, 2015.
  * [YS68] S. J. Yakowitz and J. D. Spragins.  On the identifiability of finite mixtures.  Ann. Math. Statist., 39(1):209–214, 02 1968.