# Identification of Mixtures of Discrete Product Distributions in Near-Optimal
Sample and Time Complexity

Spencer L. Gordon  Engineering and Applied Science, California Institute of
Technology, Pasadena CA 91125, USA<EMAIL_ADDRESS>,  Erik Jahn
Engineering and Applied Science, California Institute of Technology, Pasadena
CA 91125, USA<EMAIL_ADDRESS>,  Bijan Mazaheri  Engineering and Applied
Science, California Institute of Technology, Pasadena CA 91125, USA.
<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>and  Leonard J. Schulman
Engineering and Applied Science, California Institute of Technology, Pasadena
CA 91125, USA<EMAIL_ADDRESS>
###### Abstract.

We consider the problem of _identifying,_ from statistics, a distribution of
discrete random variables $X_{1},\ldots,X_{n}$ that is a mixture of $k$
product distributions. The best previous sample complexity for $n\in O(k)$ was
$(1/\zeta)^{O(k^{2}\log k)}$ (under a mild separation assumption parameterized
by $\zeta$). The best known lower bound was $\exp(\Omega(k))$.

It is known that $n\geq 2k-1$ is necessary and sufficient for identification.
We show, for any $n\geq 2k-1$, how to achieve sample complexity and run-time
complexity $(1/\zeta)^{O(k)}$. We also extend the known lower bound of
$e^{\Omega(k)}$ to match our upper bound across a broad range of $\zeta$.

Our results are obtained by combining (a) a classic method for robust tensor
decomposition, (b) a novel way of bounding the condition number of key
matrices called Hadamard extensions, by studying their action only on
flattened rank-1 tensors.

Research supported by NSF CCF-1909972 and CCF-2321079, by ISF grants 3565-21
and 389-22, and by BSF grant 2023607.

##  1\. Introduction

###  1.1. The problem and our results.

This paper resolves the sample and runtime complexity of identification of
mixtures of product distributions, a problem introduced almost thirty years
ago in [22], and further studied in [6, 12, 19, 11, 4, 31, 5, 15]. In this
problem, an observer collects samples from the distribution over $n$ binary
(or otherwise drawn from a small finite set) random variables
$X_{1},X_{2},\dots,X_{n}$. The samples are collected from a mixture of $k$
distinct sub-populations. Conditional on the sub-population
$U\in\\{1,2,\dots,k\\}$ from which a sample is drawn, the random variables are
independent. The sub-population $U$ from which a sample is collected is a
random variable chosen, for each sample of $n$ bits, independently according
to the frequencies of the sub-populations in the entire population. However,
the observer does not know these frequencies, and does not get an indication
of the sub-population from which a sample was drawn; the observer only sees
the values of the $n$ observable random bits. The goal is to reconstruct the
probabilistic model that generated the collected samples, namely, the
frequencies of the sub-populations and the conditional product distributions
on $(X_{1},X_{2},\dots,X_{n})\in\\{0,1\\}^{n}$.

Most of the above literature discusses the problem of learning the model. That
is, the goal is to produce some model with similar statistics (as measured,
for instance, by KL-divergence) on the observables as the model that generated
the samples. This does not necessarily guarantee similarity in the parameter
space of the models. In this paper, we focus on the stricter goal of
identifying the model. That is, our goal is to produce a model whose
parameters are sufficiently close to the true underlying model to generate
also similar statistics on the observables. It is known that identification is
not always possible, as there exist some distributions on the observables that
can be generated by more than one model. However, a mild condition of
separability, namely that on each observable the distribution is different
among the sub-populations, guarantees identifiability information
theoretically [31]. We shall use $\zeta$ to denote the minimum difference
between sub-population distributions on an observable. This will be defined
precisely later. Note that identification also implies learning.

Our main result is an algorithm that identifies the parameters of any
$\zeta$-separated mixture of product distributions using
$(1/\zeta)^{O(k)}(1/{\pi_{\min}})^{O(1)}(1/\varepsilon)^{2}$ samples and
runtime, up to additive error $\varepsilon$. Here ${\pi_{\min}}$ is the
minimum frequency of a sub-population. The result holds if there are at least
$2k-1$ $\zeta$-separated observables, which is known to be a necessary
condition [32, 3]. This result greatly improves upon the best previously known
complexity for identification (and learning), that required
$(1/\zeta)^{O(k^{2}\log k)}(1/{\pi_{\min}})^{O(\log k)}(1/\varepsilon)^{2}$
[15] samples from $3k-3$ $\zeta$-separated observables. Furthermore, we show
that the sample complexity of identification (for constant $\varepsilon$) is
at least $(1/(k\zeta))^{\Omega(k)}$ (note that $\zeta\leq\frac{1}{k}$ always).
This generalizes the previously known lower bound of $\exp(-\Omega(k))$ that
held only for $\zeta=\Theta(1/k)$ [29]. Hence, our results are essentially
optimal, both in terms of the number of $\zeta$-separated observables needed,
and in terms of the sample and runtime complexity (excluding the case of
$\zeta=\frac{1}{k^{1+o(1)}}$, where a small gap remains).

For large $n$, $n=\omega(k)$, if a subset of $2k-1$ $\zeta$-separated
observables is known, then the runtime bounds pick up an additional factor of
$n$ (to identify all the remaining observables). Otherwise, if the required
subset exists but is not known, then the runtime picks up an additional factor
of $n^{O(k)}$ (to enumerate over all possibilities). In both cases, the sample
complexity increases by a factor of $\log n$.

###  1.2. Related work and motivation.

The seminal work of [11] solves the learning problem for general $k$ in sample
and runtime complexity $n^{O(k^{3})}$. This was improved in [5] to
$k^{O(k^{3})}n^{O(k^{2})}$. (That paper also studied the problem of learning a
“mixture of subcubes” of the hypercube, which is the special case where each
random bit $X_{i}$ is either fixed or uniformly distributed; for this case,
they showed sample and runtime complexity of $n^{O(\log k)}$.) The
identification problem was first solved in better complexity in [15], where it
was shown how to identify a mixture of $k$ product distributions on $3k-3$
$\zeta$-separated variables with sample and runtime complexity of
$(1/\zeta)^{O(k^{2}\log k)}$. That paper reduced the problem to a special case
of identifying a mixture of $k$ distributions on independent and identically
distributed bits. The latter problem is solved using an elegant two-century-
old method of Prony [7] coupled with a robustness analysis given in [14].

As pointed out in [11] and elsewhere, the case of observables taking values in
a finite set reduces to the case of binary observables. Moreover, it is known
[10] that with polynomial overhead in the size of the range of the
observables, the problem reduces to the case of identifying the conditional
expectations of real-valued observables that are independent conditional on
the sub-population. Thus, we shall focus in this paper on this real-valued
version of the problem.

A concrete context in which identification of mixture distributions comes up
is causal inference. When data is drawn from multiple sources or sub-
populations it is said to contain a latent class [8], which mirrors an
unidentified mixture source. Standard procedures for identifying causal
effects require considering the distributions within each latent class
separately to control for potential confounding effects [28]. Such an approach
is generally impossible unless these within-source probability distributions
can be identified.

More broadly, the theory of causal inference relies at its core upon _Bayesian
networks_ of random variables [27, 30, 28]; such a network imposes conditional
independencies among random variables of the system. An important scenario is
that several latent classes are subject to the same “system mechanics” (i.e.,
Bayesian network), but have different statistics. In this case, the problem of
identifying the model is a far-reaching generalization of the problem of
identifying mixtures of product distributions. There is a recent algorithm for
the more general problem [16]; and that algorithm uses as an essential (and
complexity-bottleneck) subroutine, any algorithm for identifying mixtures of
product distributions. Thus the improvements of the present paper, carry over
directly to that application.

###  1.3. Our methods.

We study the so-called Hadamard extensions that were also used to derive the
bound in [15]. We give a new and much more powerful bound on the condition
number of the Hadamard extensions. This bound alone would improve the sample
complexity of the algorithm in [15] to $(1/\zeta)^{O(k\log k)}$. We gain
further improvement as follows. Instead of reducing the problem to identifying
mixtures of iid (synthetic) bits, we reduce the problem to a tensor
decomposition problem, where the tensor components are guaranteed to be well-
conditioned (which makes the decomposition unique). An algorithm for tensor
decomposition in this setting was given thirty years ago in [25] and later
analyzed for robustness in [18, 2]. (This algorithm can also be seen as a
generalization of the matrix pencil method [21], applied to the iid case in
[23].) Adapting the tensor decomposition algorithm to our setting and
analyzing it using our new condition number bound, then yields the sample
complexity of $(1/\zeta)^{O(k)}$.

###  1.4. Comparison with the parametric case.

The literature on mixture models for parametric families (exponential
distributions, Gaussians in $\mathbb{R}$ or $\mathbb{R}^{d}$, etc.) is even
more extensive and older than for discrete mixture models. It is essential to
realize a fundamental difference between the types of problems. In general,
data is generated by (unseen) selection of a sub–population $j$ $(1\leq j\leq
k)$, followed by (seen) sampling of $n$ independent samples _from the $j$-th
distribution_. In almost every parametric scenario (think e.g., of a mixture
of $k$ Gaussians or exponential distributions on the line), _$n=1$ is
sufficient_ in order to (in the limit of many repetitions) identify the model.
This is fundamentally untrue in the non-parametric case; we have already
mentioned that a lower bound of $n\geq 2k-1$ was shown in [29]; this threshold
for $n$ is called there the “aperture” of the problem. To see, for starters,
why the aperture must be larger than $n=1$, consider a single binary variable
with $k=2$ equiprobable sources (i.e., $\Pr(U=0)=\Pr(U=1)=1/2$), one of which
has $\Pr(X=1\mid U=0)=\frac{3}{4}$ and the other $\Pr(X=1\mid
U=1)=\frac{1}{4}$. If we see after each selection of a source only a single
sample of $X$, it is impossible to distinguish between the above mixture and a
mixture in which $\Pr(X=1\mid U=0)=1$ and $\Pr(X=1\mid U=1)=0$. With access to
multiple independent samples from the _same_ source, however, we get empirical
estimates of higher moments of the distribution, and at the critical aperture
can identify the model.

###  1.5. Organization.

Section 2 formally states the identification problem for mixtures of product
distributions and sets up the key mathematical objects needed for our work.
Section 3 describes our algorithm and states our upper bounds on its sample
complexity. Along with the algorithm pointers are provided to the main steps
of the analysis in Sections 4, 5. In section 6 we prove lower bounds on the
sample complexity of the identification problem. Finally, in section 7 we
discuss potential further directions of research on our topic.

##  2\. Results and preliminaries

###  2.1. The k-MixProd problem

Consider $n$ real, compactly supported random variables $X_{1},\dots,X_{n}$
that are independent conditional on a latent random variable $U$ with range
$[k]=\\{1,\dots,k\\}$. Given iid samples of the joint distribution of
$(X_{1},\dots,X_{n})$, we want to identify the distribution of $U$, given by
$\pi_{j}:=\Pr(U=j)\;(j\in[k])$, and the conditional expectations
${\mathbf{m}}_{ij}=\operatorname{\mathbb{E}}(X_{i}\mid
U=j)\;(i\in[n],j\in[k])$. Hence, the model parameters for our problem are
given by a vector $(\pi,{\mathbf{m}})\in\Delta^{k-1}\times\mathbb{R}^{n\times
k}$, where $\Delta^{k-1}$ denotes the $(k-1)$-simplex.

Set $X_{S}=\prod_{i\in S}X_{i}$, so $\operatorname{\mathbb{E}}(X_{S}\mid
U=j)=\prod_{i\in S}{\mathbf{m}}_{ij}$. The mapping of the model to the
statistics is then given by:

(1) $\displaystyle\gamma_{n}:\Delta^{k-1}\times[0,1]^{n\times
k}\to\mathbb{R}^{2^{[n]}}$ (2)
$\displaystyle\gamma_{n}(\pi,{\mathbf{m}})(S)=\operatorname{\mathbb{E}}(X_{S})=\sum_{j=1}^{k}\pi_{j}\operatorname{\mathbb{E}}(X_{S}\mid
U=j)=\sum_{j=1}^{k}\pi_{j}\prod_{i\in S}{\mathbf{m}}_{ij}$

We drop the subscript $n$ and write $\gamma$ when $n$ is implied. The
$k$-MixProd identification problem is to invert $\gamma_{n}$, i.e., to recover
$(\pi_{j})_{j\in[k]}$ and $({\mathbf{m}}_{ij})_{i\in[n],j\in[k]}$ (up to
permuting the set $[k]$). This task is interesting in two versions, exact
identification of $(\pi,{\mathbf{m}})$ from $\gamma_{n}(\pi,{\mathbf{m}})$
(i.e., from perfect statistics), and approximate identification of
$(\pi,{\mathbf{m}})$ from noise-perturbed statistics
$\operatorname{{\tilde{g}}}$, i.e., from
$\operatorname{{\tilde{g}}}\in\mathbb{R}^{\\{0,1\\}^{n}}$ that is close to
$\gamma_{n}(\pi,{\mathbf{m}})$. To make the latter goal precise we need to
specify metrics on the domain and range of $\gamma_{n}$. These are
$L_{\infty}$ metrics, up to relabelings of the latent variable. ($S_{k}$
denotes the symmetric group on $k$ letters.)

(3) $\displaystyle
d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\pi^{\prime},{\mathbf{m}}^{\prime}))$
$\displaystyle:=\min_{\rho\in
S_{k}}\max\\{\max_{j}|\pi_{j}-\pi^{\prime}_{\rho(j)}|,\max_{i,j}|{\mathbf{m}}_{i,j}-{\mathbf{m}}^{\prime}_{i,\rho(j)}|\\}$
(4) $\displaystyle
d_{\mathrm{stat}}(\operatorname{g},\operatorname{g}^{\prime})$
$\displaystyle:=\max_{S\subseteq[n]}|\operatorname{g}(S)-\operatorname{g}^{\prime}(S)|.$

The mapping $\gamma_{n}$ is not everywhere injective, so the $k$-MixProd model
identification problem is not always feasible. To guarantee identifiability we
need to make the following assumptions:

  1. (a)

($\zeta$-separation) each variable $X_{i}$ is $\zeta$-_separated_ , i.e.
$|{\mathbf{m}}_{ij}-{\mathbf{m}}_{ij^{\prime}}|\geq\zeta$ for all $j\neq
j^{\prime}\in[k]$;

  2. (b)

(non-degenerate prior) for each $j\in[k]$, we have
$\pi_{j}\geq{\pi_{\min}}>0$;

  3. (c)

(sufficiently many observables) there are at least $n\geq 2k-1$ variables
$X_{i}$.

Let $\mathcal{D}_{n,\zeta,{\pi_{\min}}}$ denote the space of the $k$-MixProd
models with $n$ variables satisfying assumptions (a) and (b). Formally:

(5)
$\displaystyle\mathcal{D}_{n,\zeta,{\pi_{\min}}}=\\{(\pi,{\mathbf{m}})\in\Delta^{k-1}\times[0,1]^{n\times
k}\mid\min_{j}\pi_{j}\geq{\pi_{\min}},\;\forall i\;\min_{j\neq
j^{\prime}}|{\mathbf{m}}_{ij}-{\mathbf{m}}_{ij^{\prime}}|\geq\zeta\\},$

Theorem 11 shows (quite apart from its algorithmic content) that for $n\geq
2k-1$, if $(\pi,{\mathbf{m}})$ is a model in
$\mathcal{D}_{n,\zeta,{\pi_{\min}}}$, then any model whose statistics are
close (in $d_{\mathrm{stat}}$) to those of $(\pi,{\mathbf{m}})$, must also be
close to $(\pi,{\mathbf{m}})$ in $d_{\mathrm{model}}$.

We now introduce some mathematical concepts that will be needed for our work.

###  2.2. Hadamard extensions and related definitions

Subsets of $[n]$ will typically be denoted by a capitalized variable in
ordinary font: $S\subseteq[n]$.

######  Definition 1.

Given a matrix ${\mathbf{A}}$ of any dimensions, let ${\mathbf{A}}_{i*}$
denote the $i$’th row of ${\mathbf{A}}$, and ${\mathbf{A}}_{*j}$ the $j$’th
column of ${\mathbf{A}}$. Where clear from context we write ${\mathbf{A}}_{i}$
instead of ${\mathbf{A}}_{i*}$. For $S$ a set of rows, ${\mathbf{A}}[S]$
denotes the submatrix of ${\mathbf{A}}$ consisting of the rows in $S$.

######  Definition 2 (Hadamard product).

The _Hadamard product_ is the mapping
$\odot:\mathbb{R}^{[k]}\times\mathbb{R}^{[k]}\to\mathbb{R}^{[k]}$ which, for
row vectors $u=\left(u_{1},\dotsc,u_{k}\right)$ and
$v=\left(v_{1},\dotsc,v_{k}\right)$, is given by $u\odot
v\coloneqq(u_{1}v_{1},\ldots,u_{k}v_{k})$. Equivalently, using the notation
$v_{\odot}=\operatorname{diag}(v)$, the Hadamard product is $u\odot v=u\cdot
v_{\odot}$. The identity element for the Hadamard product is the all-ones row
vector ${\mathbb{1}}$.

######  Definition 3 (Hadamard extension).

For ${\mathbf{n}}\in\mathbb{R}^{n\times p}$, the _Hadamard extension_ of
${\mathbf{n}}$, written $\mathbb{H}({\mathbf{n}})$, is the $2^{n}\times p$
matrix with rows $\mathbb{H}({\mathbf{n}})_{S}$ for all $S\subseteq[n]$,
where, for $S=\\{i_{1},\ldots,i_{\ell}\\}$,
$\mathbb{H}({\mathbf{n}})_{S}={\mathbf{n}}_{i_{1}}\odot\cdots\odot{\mathbf{n}}_{i_{\ell}}$;
equivalently $\mathbb{H}({\mathbf{n}})_{S,j}=\prod_{i\in S}{\mathbf{n}}_{ij}$.
In particular $\mathbb{H}({\mathbf{n}})_{\emptyset}={\mathbb{1}}$, and for all
$i\in[n]$, $\mathbb{H}({\mathbf{n}})_{\\{i\\}}={\mathbf{n}}_{i}$.

This construction first to our knowledge appeared (not under this name) in
[5].

######  Definition 4.

The singular values of a real matrix ${\mathbf{A}}$ are denoted
$\sigma_{1}({\mathbf{A}})\geq\sigma_{2}({\mathbf{A}})\geq\dotsm$. The
$L_{2}\to L_{2}$ operator norm is denoted $\left\|{\mathbf{A}}\right\|$. The
condition number of ${\mathbf{A}}$ is denoted
$\kappa({\mathbf{A}})=\left\|{\mathbf{A}}\right\|\cdot\left\|{\mathbf{A}}^{-1}\right\|$.

######  Definition 5 (Vandermonde matrix).

The Vandermonde matrix $\operatorname{Vdm}(m)\in\mathbb{R}^{k\times k}$
associated with a row vector $m\in\mathbb{R}^{k}$ has entries
$\operatorname{Vdm}(m)_{ij}=(m_{j})^{i}$ for
$i\in\left\\{0,1,\dotsc,k-1\right\\}$ and $j\in\left\\{1,2,\dotsc,k\right\\}$.
We also write $\operatorname{Vdm}(m,r)$ for the $r\times k$ matrix with
entries $\operatorname{Vdm}(m,r)_{ij}=(m_{j})^{i}$.

###  2.3. Multilinear moments

The data we obtain from our samples will be estimates of
$\operatorname{\mathbb{E}}[X_{S}]=\operatorname{\mathbb{E}}[\prod_{i\in
S}X_{i}]$ for all subsets $S\subseteq[n]$. We call these the _multilinear
moments_ of the distribution, since they are multilinear in the rows
${\mathbf{m}}_{i}$. Observe that
$\operatorname{\mathbb{E}}[X_{S}]=\sum_{j}\pi_{j}\prod_{i\in
S}{\mathbf{m}}_{ij}=(\mathbb{H}({\mathbf{m}}))_{S}\cdot\pi$, or equivalently,
$\operatorname{\mathbb{E}}[X_{S}]=({\mathbf{m}}_{i_{1}}\odot{\mathbf{m}}_{i_{2}}\odot\dotsb\odot{\mathbf{m}}_{i_{s}})\pi$
where $S=\\{i_{1},i_{2},\dotsc,i_{s}\\}$. Hence, the vector of statistics for
a model $(\pi,{\mathbf{m}})$ is given by
$\gamma(\pi,{\mathbf{m}})=\mathbb{H}({\mathbf{m}})\pi$. Observe that source
identification is not possible if $\mathbb{H}({\mathbf{m}})$ has less than
full column rank, i.e., if $\operatorname{rank}(\mathbb{H}({\mathbf{m}}))<k$,
as then the mixing weights cannot be unique.

######  Definition 6.

Given disjoint sets $S,T\subseteq\\{2,\dots,n\\}$, define

$\displaystyle\mathbf{C}_{ST}$
$\displaystyle=\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}},$
$\displaystyle\mathbf{C}_{ST,1}$
$\displaystyle=\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot{\mathbf{m}}_{1\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}$

Note, for $A\subseteq S,B\subseteq T$,

$\displaystyle(\mathbf{C}_{ST})_{A,B}$
$\displaystyle=\gamma_{n}(\pi,{\mathbf{m}})(A\cup B),$
$\displaystyle(\mathbf{C}_{ST,1})_{A,B}$
$\displaystyle=\gamma_{n}(\pi,{\mathbf{m}})(A\cup B\cup\\{1\\})$

Consequently $\mathbf{C}_{ST}$ and $\mathbf{C}_{ST,1}$ are observable, that
is, every one of their entries is a statistic which the algorithm receives (a
noisy version of) as input.

###  2.4. Tensor decomposition

A matrix ${\mathbf{A}}$ is rank $1$ if and only if it can be written as
${\mathbf{A}}=\mathbf{uv}^{\mathsf{T}}$ for some vectors
$\mathbf{u},\mathbf{v}$. This concept can be generalized for tensors:

######  Definition 7.

A $3$-way tensor ${\mathcal{T}}\in\mathbb{R}^{d_{1}\times d_{2}\times d_{3}}$
is said to be of rank $1$ if there exist vectors
$u\in\mathbb{R}^{d_{1}},v\in\mathbb{R}^{d_{2}},z\in\mathbb{R}^{d_{3}}$ such
that for all $i,j,k$:

${\mathcal{T}}_{ijk}=u_{i}\cdot v_{j}\cdot z_{k}.$

Now, the rank of any tensor ${\mathcal{T}}$ can be defined as the minimum
number of rank-1-tensors that sum up to ${\mathcal{T}}$. Equivalently, a
tensor of rank $r$ has the following decomposition:

######  Definition 8.

A $3$-way tensor ${\mathcal{T}}\in\mathbb{R}^{d_{1}\times d_{2}\times d_{3}}$
has a rank-$r$-decomposition if there exist matrices
$\mathbf{U}\in\mathbb{R}^{d_{1}\times r},\mathbf{V}\in\mathbb{R}^{d_{2}\times
r},\mathbf{Z}\in\mathbb{R}^{d_{3}\times r}$ such that for all $i,j,k$:

${\mathcal{T}}_{ijk}=\sum_{\ell=1}^{r}\mathbf{U}_{i\ell}\mathbf{V}_{j\ell}\mathbf{Z}_{k_{\ell}}.$

We write ${\mathcal{T}}=[\mathbf{U},\mathbf{V},\mathbf{Z}]$ and call
$\mathbf{U},\mathbf{V},\mathbf{Z}$ the factor matrices or tensor components of
${\mathcal{T}}$.

In general, the rank-$r$-decomposition of a tensor ${\mathcal{T}}$ need not be
unique, but a classical result of Kruskal [24] gives sufficient conditions for
uniqueness.

######  Definition 9.

The Kruskal rank of a matrix ${\mathbf{A}}$ is the largest number $r$ such
that any $r$ columns of ${\mathbf{A}}$ are linearly independent.

######  Theorem 10 (Kruskal [24]).

The rank-$r$-decomposition of a three-way tensor
${\mathcal{T}}=[\mathbf{U},\mathbf{V},\mathbf{Z}]$ is unique up to scaling and
permuting the columns of the factor matrices if

$k_{\mathbf{U}}+k_{\mathbf{V}}+k_{\mathbf{Z}}\geq 2r+2,$

where $k_{\mathbf{U}},k_{\mathbf{V}},k_{\mathbf{Z}}$ denote the Kruskal rank
of the matrices $\mathbf{U},\mathbf{V},\mathbf{Z}$ respectively.

##  3\. The algorithm

###  3.1. Reducing $k$-MixProd to tensor decomposition

To motivate our algorithm, we first discuss a way of solving the $k$-MixProd
identification problem for $n\geq 2k-1$ given perfect statistics. Consider
three disjoint sets $S,T,U\subseteq[n]$ of $\zeta$-separated observables, such
that $|S|=|T|=k-1$ and $|U|=1$. For convenience, we index rows so that
$U=\\{1\\}$. Consider the vector of perfect statistics
$\operatorname{g}=\gamma_{2k-1}({\mathbf{m}}[S\cup T\cup\\{1\\}])$ that
corresponds to the observables in $S\cup T\cup\\{1\\}$. We can naturally view
$\operatorname{g}$ as a three-way tensor
${\mathcal{T}}\in\mathbb{R}^{2^{S}\times 2^{T}\times 2^{\\{1\\}}}$ whose
entries are given by ${\mathcal{T}}_{A,B,C}=\operatorname{g}_{A\cup B\cup C}$
for $A\subseteq S,B\subseteq T,C\subseteq\\{1\\}$. Since we have

$\displaystyle{\mathcal{T}}_{A,B,C}=\sum_{j=1}^{k}\pi_{j}\prod_{i\in A\cup
B\cup
C}{\mathbf{m}}_{ij}=\sum_{j=1}^{k}\mathbb{H}({\mathbf{m}}[S])_{A,j}\cdot\mathbb{H}({\mathbf{m}}[T])_{B,j}\cdot(\mathbb{H}({\mathbf{m}}_{1})\cdot\pi_{\odot})_{C,j},$

the tensor ${\mathcal{T}}$ can be decomposed as
${\mathcal{T}}=\left[\mathbb{H}({\mathbf{m}}[S]),\mathbb{H}({\mathbf{m}}[T]),\mathbb{H}({\mathbf{m}}_{1})\pi_{\odot}\right]$.
It will follow as a “qualitative” corollary of Theorem 14 that
$\mathbb{H}({\mathbf{m}}[S])$ and $\mathbb{H}({\mathbf{m}}[T])$ have full
column rank, which implies both matrices have Kruskal rank $k$. Moreover,
$\zeta$-separation of $X_{1}$ implies that the matrix
$\mathbb{H}({\mathbf{m}}_{1})\cdot\pi_{\odot}$ has Kruskal rank $2$. Hence, by
Theorem 10, the decomposition of ${\mathcal{T}}$ is unique. We deduce that
identifying the model parameters ${\mathbf{m}}$ and $\pi$ from perfect
statistics (provided in ${\mathcal{T}}$) is equivalent to computing the unique
tensor decomposition of ${\mathcal{T}}$. An efficient algorithm for computing
tensor decomposition in this setting has first been given by [25] and later
analyzed for stability by [18], [2]. The idea is to first project the
components of ${\mathcal{T}}$ down to their image, i.e. find matrices
$\mathbf{U},\mathbf{V}\in\mathbb{R}^{2^{k-1}\times k}$ such that
$\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])$ and
$\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])$ are invertible. Then,
compute
${\hat{\mathbf{C}}}_{ST}=\mathbf{U}^{\mathsf{T}}\mathbf{C}_{ST}\mathbf{V}$ and
${\hat{\mathbf{C}}}_{ST,1}=\mathbf{U}^{\mathsf{T}}\mathbf{C}_{ST,1}\mathbf{V}$
(see Definition 6) from the given statistics. Now, the key observation is that
the tensor components $\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])$ and
$\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])$ can be found as the
eigenvectors of ${\hat{\mathbf{C}}}_{ST,1}{\hat{\mathbf{C}}}_{ST}^{-1}$ and
${\hat{\mathbf{C}}}_{ST,1}^{\mathsf{T}}({\hat{\mathbf{C}}}_{ST}^{\mathsf{T}})^{-1}$
respectively. This is because

$\displaystyle{\hat{\mathbf{C}}}_{ST,1}{\hat{\mathbf{C}}}_{ST}^{-1}$
$\displaystyle=\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot{\mathbf{m}}_{1\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\mathbf{V}\cdot\left(\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\mathbf{V}\right)^{-1}$
$\displaystyle=\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot{\mathbf{m}}_{1\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\mathbf{V}\cdot\left(\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\mathbf{V}\right)^{-1}\cdot\pi_{\odot}^{-1}\cdot\left(\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\right)^{-1}$
(6)
$\displaystyle=\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\cdot{\mathbf{m}}_{1\odot}\cdot\left(\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\right)^{-1},$

and

$\displaystyle{\hat{\mathbf{C}}}_{ST,1}^{\mathsf{T}}({\hat{\mathbf{C}}}_{ST}^{\mathsf{T}})^{-1}$
$\displaystyle=\left(\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot{\mathbf{m}}_{1\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\mathbf{V}\right)^{\mathsf{T}}\cdot\left(\left(\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\cdot\pi_{\odot}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\mathbf{V}\right)^{\mathsf{T}}\right)^{-1}$
$\displaystyle=\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])\cdot\pi_{\odot}\cdot{\mathbf{m}}_{1\odot}\cdot\mathbb{H}({\mathbf{m}}[S])^{\mathsf{T}}\mathbf{U}\cdot\left(\mathbb{H}({\mathbf{m}}[S])^{\mathsf{T}}\mathbf{U}\right)^{-1}\cdot\pi_{\odot}^{-1}\cdot\left(\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])\right)^{-1}$
(7)
$\displaystyle=\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])\cdot{\mathbf{m}}_{1\odot}\cdot\left(\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])\right)^{-1}.$

In both cases, the eigenvalues of the matrices above are given by the entries
of ${\mathbf{m}}_{1}$. Crucially, $\zeta$-separation of ${\mathbf{m}}_{1}$
allows us to match up the columns of
$\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])$ and
$\mathbf{V}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])$ (and guarantees numerical
stability). The original entries of $\mathbb{H}({\mathbf{m}})$ and $\pi$ can
then be found from linear systems. In fact, notice that

(8)
$\displaystyle\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\pi=\mathbf{U}^{\mathsf{T}}\operatorname{g}[2^{S}],$

and for any row ${\mathbf{m}}_{i}$ with $i\notin S$, we have

(9)
$\displaystyle\mathbf{U}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\pi_{\odot}{\mathbf{m}}_{i}^{\mathsf{T}}=\mathbf{U}^{\mathsf{T}}\operatorname{g}(R\cup\\{i\\})_{R\subseteq
S}.$

At this point, $\pi$ and ${\mathbf{m}}_{i}$ are the only unknowns in the
equations above. We now give Algorithm 1, which performs the identification
procedure described above.

Algorithm 1 (adapted from [25], [2]) Identifies a mixture of product
distributions on $2k-1$ binary variables given the joint distribution.

⬇

1 Input: Two disjoint sets $S,T$ of $\zeta$-separated, binary observables
$X_{i}$, each of cardinality $k-1$; a single $\zeta$-separated, binary
observable that is not part of $S,T$, wlog $X_{1}$; vector
$\operatorname{{\tilde{g}}}\in\mathbb{R}^{2^{S\cup T\cup\\{1\\}}}$
($\operatorname{{\tilde{g}}}(R)$ is the empirical approximation to
$E(X_{R})=\gamma(\pi,{\mathbf{m}})(R)$).

2 Construct ${\tilde{\mathbf{C}}}_{ST}$ and ${\tilde{\mathbf{C}}}_{ST,1}$ by
$({\tilde{\mathbf{C}}}_{ST})_{AB}=\operatorname{{\tilde{g}}}(A\cup B)$ and
$({\tilde{\mathbf{C}}}_{ST,1})_{AB}=\operatorname{{\tilde{g}}}(A\cup
B\cup\\{1\\})$ for $A\subseteq S,B\subseteq T$.

3 Set $\hat{\mathbf{U}}\in\mathbb{R}^{2^{|T|}\times k}$ to be the top $k$ left
singular vectors of ${\tilde{\mathbf{C}}}_{ST}$ and
$\hat{\mathbf{V}}\in\mathbb{R}^{2^{|S|}\times k}$ to be the top $k$ right
singular vectors of ${\tilde{\mathbf{C}}}_{ST}$.

4
${\hat{\mathbf{C}}}_{ST}\leftarrow\hat{\mathbf{U}}^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\hat{\mathbf{V}}$,
${\hat{\mathbf{C}}}_{ST,1}\leftarrow\hat{\mathbf{U}}^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST,1}\hat{\mathbf{V}}$

5 Set $\hat{\mathbf{S}}$ to be the eigenvectors of
${\hat{\mathbf{C}}}_{ST,1}({\hat{\mathbf{C}}}_{ST})^{-1}$ (sorted from highest
eigenvalue to lowest).

6 Set $\hat{\mathbf{T}}$ to be the eigenvectors of
${\hat{\mathbf{C}}}_{ST,1}^{\mathsf{T}}({\hat{\mathbf{C}}}_{ST}^{\mathsf{T}})^{-1}$
(sorted from highest eigenvalue to lowest).

7
$\tilde{\pi}\leftarrow\hat{\mathbf{S}}^{-1}\cdot\hat{\mathbf{U}}^{\mathsf{T}}\left(\operatorname{{\tilde{g}}}(R)_{R\subseteq
S}\right)$

8 for every $i\in T\cup\\{1\\}$,
$\tilde{{\mathbf{m}}}_{i}\leftarrow\left((\operatorname{{\tilde{g}}}(R\cup\left\\{i\right\\}))_{R\subseteq
S}\right)^{\mathsf{T}}\cdot\hat{\mathbf{U}}\cdot{\left(\hat{\mathbf{S}}^{\mathsf{T}}\right)}^{-1}\cdot\tilde{\pi}_{\odot}^{-1}$.

9 for every $i\in S$,
$\tilde{{\mathbf{m}}_{i}}\leftarrow\left((\operatorname{{\tilde{g}}}(R\cup\left\\{i\right\\}))_{R\subseteq
T}\right)^{\mathsf{T}}\cdot\hat{\mathbf{V}}\cdot{\left(\hat{\mathbf{T}}^{\mathsf{T}}\right)}^{-1}\cdot\tilde{\pi}_{\odot}^{-1}$.

_Some comments on the algorithm:_ first, it is not necessary that the sets
$S,T$ are of size $k-1$; typically, $\lceil\lg k\rceil$ will suffice. It is
not even necessary that observables are $\zeta$-separated. These assumptions
guarantee that $\sigma_{k}(\mathbf{C}_{ST})$ is large, but the latter, along
with the argument that this is w.h.p. reproduced for
$\sigma_{k}({\tilde{\mathbf{C}}}_{ST})$, is sufficient for the success of the
algorithm.

Second, we do not need to start with knowledge of $S,T$. Given $n$ variables
of which an unknown subset of $2k-1$ variables are $\zeta$-separated, we can
simply perform the algorithm for all possible choices of subsets $S,T$ of size
$k-1$ (and an additional single row). Then, we choose the computed model whose
statistics are closest to the observed statistics as the final output. This
exhaustive search can increase the runtime by a factor of about $n^{2k}$; but
actually all these complexities are only exponential in the _actual_ needed
size of $S\cup T$, which as noted, for generic ${\mathbf{m}}$ will be as small
as $\lceil\lg k\rceil$, which makes the algorithm far more attractive in
practice.

Our main result essentially states that Algorithm 1 is robust to noise, i.e.
still performs identification accurately when only presented with approximate
statistics.

######  Theorem 11.

Let $n=2k-1$ and fix any ${\varepsilon}\in(0,\zeta/2)$. Let
$(\pi,{\mathbf{m}})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}}$. If Algorithm 1 is
given approximate statistics $\operatorname{{\tilde{g}}}$ on
$(X_{1},\ldots,X_{n})$ as input, satisfying
$d_{\mathrm{stat}}(\gamma(\pi,{\mathbf{m}}),\operatorname{{\tilde{g}}})<{\pi_{\min}}^{O(1)}\zeta^{O(k)}{\varepsilon}$,
then in runtime $\exp(O(k))$ the algorithm outputs
$(\tilde{\pi},\tilde{{\mathbf{m}}})$ s.t.

(10)
$d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\tilde{\pi},\tilde{{\mathbf{m}}}))<{\varepsilon}.$

Moreover, this output is essentially unique, in the sense that: any (not
necessarily $\zeta$-separated) model $(\pi^{\prime},{\mathbf{m}}^{\prime})$
with
$d_{\mathrm{stat}}(\gamma(\pi,{\mathbf{m}}),\gamma(\pi^{\prime},{\mathbf{m}}^{\prime}))<{\pi_{\min}}^{O(1)}\zeta^{O(k)}{\varepsilon}$
also satisfies (10).

######  Corollary 12.

For $n=2k-1$ random variables $X_{i}$, sample complexity

$(1/\zeta)^{O(k)}(1/{\pi_{\min}})^{O(1)}(1/\varepsilon)^{2}$

suffices to compute a model that w.h.p. satisfies (10).

######  Corollary 13.

Let the number of observables be $n\geq 2k-1$, and of these let some $2k-1$ be
$\zeta$-separated. If this subset is known, then sample complexity

$\log n\cdot(1/\zeta)^{O(k)}(1/{\pi_{\min}})^{O(1)}(1/\varepsilon)^{2}$

and post-sampling runtime $n\cdot\exp(O(k))$ suffices to compute a model that
w.h.p. satisfies (10). If the subset is not known, then the same sample
complexity and post-sampling runtime $n^{2k}\cdot\exp(O(k))$ suffices to
compute a model that w.h.p. satisfies (10).

Corollary 12 follows from Theorem 11 by standard Chernoff bounds and a union
bound. The proofs of Theorem 11 and Corollary 13 are in Section 5.

###  3.2. Annotated steps of the algorithm

This section describes what each non-trivial line of the algorithm (Alg. 1)
accomplishes, and in each case points to the part of the analysis necessary to
justify it. We first set some definitions for the analysis. Let $n=2k-1$, let
$\\{2,\dots,n\\}$ be the disjoint union of sets $S,T$ each of size $k-1$. Let
$\operatorname{g}=\gamma_{n}(\pi,{\mathbf{m}})$ and let
$\operatorname{{\tilde{g}}}\in\mathbb{R}^{\\{0,1\\}^{n}}$ be the empirical
statistics. The analysis relies on assuming that
$(\pi,{\mathbf{m}})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}}$. In the sequel let

(11)
${\mathbf{d}}:=d_{\mathrm{stat}}(\operatorname{{\tilde{g}}},\operatorname{g}).$

Now for the line-by-line:

  * •

Line 3: The SVD of ${\tilde{\mathbf{C}}}_{ST}$ can be computed with high
numerical accuracy in time $\exp(O(k))$ using, for instance, the Golub-Kahan-
Reinsch algorithm (see Lemma 18 and [13], chapter 8). The span of the top $k$
left and right singular vectors of ${\tilde{\mathbf{C}}}_{ST}$ approximate the
images of $\mathbb{H}({\mathbf{m}}[S])$ and $\mathbb{H}({\mathbf{m}}[T])$
respectively.

  * •

Line 4: We project on the top $k$ singular vectors of
${\tilde{\mathbf{C}}}_{ST}$ to get an invertible matrix
${\hat{\mathbf{C}}}_{ST}$. Lemma 19 shows that ${\hat{\mathbf{C}}}_{ST}$ is
well-conditioned, given that ${\mathbf{d}}$ is small enough. Note that the two
matrices ${\hat{\mathbf{C}}}_{ST}$ and ${\hat{\mathbf{C}}}_{ST,1}$ can be
arranged to a $k\times k\times 2$-tensor $\hat{{\mathcal{T}}}$ that is close
to the tensor
${\mathcal{T}}=\left[\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]),\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T]),\mathbb{H}({\mathbf{m}}_{1})\pi_{\odot}\right]$
by Lemma 20.

  * •

Lines 5,6: These two lines implement the core of the tensor decomposition
algorithm from [25], as explained in the previous section. Hence, the matrices
$\hat{\mathbf{S}}$ and $\hat{\mathbf{T}}$ approximate the tensor components
$\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])$ and
$\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])$ of ${\mathcal{T}}$.
The diagonalization steps can be performed with high numerical accuracy in
time $\exp(O(k))$ using, for instance, the algorithm from [1]. For the error
bounds, we rely on Theorem 21 from [2].

  * •

Lines 7-9: Here, we solve for the model parameters $\pi$ and ${\mathbf{m}}$,
using the equations (8) and (9). For error control, we apply Lemma 22.

##  4\. The condition number bound

The key to the sample-complexity and runtime bounds for our algorithm lies in
the following condition number bound for the Hadamard Extension.

######  Theorem 14.

  1. (1)

Let ${\mathbf{m}}$ consist of $k-1$ $\zeta$-separated rows in
$\mathbb{R}^{k}$, and observe that the singular values of
$\mathbb{H}({\mathbf{m}})=\mathbb{H}({\mathbf{m}})$ satisfy
$\sigma_{1}(\mathbb{H}({\mathbf{m}}))\geq\ldots\geq\sigma_{k}(\mathbb{H}({\mathbf{m}}))\geq
0=\sigma_{k+1}(\mathbb{H}({\mathbf{m}}))=\ldots=\sigma_{2^{k-1}}(\mathbb{H}({\mathbf{m}}))$.
Then

(12)
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))>\frac{1}{\sqrt{k}}\left(\frac{\zeta}{2\sqrt{5}}\right)^{k-1}=:{\bm{\sigma}}.$

  2. (2)

Let $(\pi,{\mathbf{m}})\in\mathcal{D}_{2k-1,\zeta,{\pi_{\min}}}$ and let
$\mathbf{C}_{ST}$ be as in Defn. 6. Then

$\sigma_{k}(\mathbf{C}_{ST})>{\pi_{\min}}{\bm{\sigma}}^{2}=\frac{{\pi_{\min}}}{k}\left(\frac{\zeta}{2\sqrt{5}}\right)^{2k-2}.$

######  Definition 15.

$\mathbb{H}_{p}=\\{\mathbb{H}({\mathbf{n}}):{\mathbf{n}}\in\mathbb{R}^{[k-1]\times[p]}\\}$.
(So $\mathbb{H}_{1}$ consists of rank-1 tensors of order $k-1$.)

The proof of Theorem 14 relies on the following insight. Since
$\mathbb{H}({\mathbf{m}})$ has dimensions $2^{k-1}\times k$,
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))$ characterizes the least norm of
$\mathbb{H}({\mathbf{m}})\cdot v$ ranging over any unit vector $v$ (all norms
in $L_{2}$), but it does _not_ characterize the least norm of vectors of the
form $h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})$, which is $0$ as the left-
kernel of $\mathbb{H}({\mathbf{m}})$ is of course very large. The insight is
that it _does_ become possible to bound $\sigma_{k}(\mathbb{H}({\mathbf{m}}))$
in terms of such vectors $h$, _provided $h$ is restricted to rank $1$
tensors._ With this in mind we define:

$\tau({\mathbf{m}})=\min_{0\neq
h\in\mathbb{H}_{1}}\left\|h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})\right\|/\left\|h\right\|.$

Proof of Theorem 14. To show Part 2 from Part 1: The SVD implies there is a
$k$-dimensional space $V$ s.t. $\forall v\in V$,
$\|v^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}}[S])\|\geq{\bm{\sigma}}\|v\|$.
Further, for all $w\in\mathbb{R}^{k}$,
$\|w^{\mathsf{T}}\cdot\pi_{\odot}\|\geq{\pi_{\min}}\|w\|$. And for all
$w\in\mathbb{R}^{k}$,
$\|w^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\|\geq{\bm{\sigma}}\|w\|$.
So $\forall v\in V$,
$\|v\cdot\mathbf{C}_{TS}\|\geq{\pi_{\min}}{\bm{\sigma}}^{2}\|v\|$.

In order to establish Part 1 we prove the following two lemmas.

######  Lemma 16.

$\sigma_{k}(\mathbb{H}({\mathbf{m}}))\geq\tau({\mathbf{m}})/\sqrt{k}$.

######  Lemma 17.

$\tau({\mathbf{m}})>(\zeta/2\sqrt{5})^{k-1}$.

Proof of Lemma 16. Consider $v\in\mathbb{R}^{k}$, $\|v\|=1$, achieving
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))$, i.e.,
$r\coloneqq\mathbb{H}({\mathbf{m}})\cdot v$ satisfies
$\left\|r\right\|=\sigma_{k}(\mathbb{H}({\mathbf{m}}))$. W.l.o.g. the order of
coordinates is such that $|v_{k}|\geq 1/\sqrt{k}$. The last column of
$\mathbb{H}({\mathbf{m}})$ is then:

(13)
$\mathbb{H}({\mathbf{m}})_{*k}=\frac{1}{v_{k}}\left(r-\sum_{j=1}^{k-1}v_{j}\mathbb{H}({\mathbf{m}})_{*j}\right).$

Now we carefully choose $h\in\mathbb{H}_{1}$ based on ${\mathbf{m}}$. Define
the column vector ${\mathbf{n}}\in\mathbb{R}^{[k-1]}$ by
${\mathbf{n}}_{i}\coloneqq-1/{\mathbf{m}}_{ii}$ ($1\leq i\leq k-1$); and let

(14)
$h\coloneqq\left(\prod_{1}^{k-1}{\mathbf{m}}_{ii}\right)\mathbb{H}({\mathbf{n}})=({\mathbf{m}}_{11};-1)\odot\ldots\odot({\mathbf{m}}_{k-1,k-1};-1).$

For $j\neq k$ we have:

$\displaystyle h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})_{*j}$
$\displaystyle=$ $\displaystyle\sum_{S}(-1)^{|S|}\left(\prod_{i\notin
S}{\mathbf{m}}_{ii}\right)\left(\prod_{i\in
S}{\mathbf{m}}_{ij}\right)=\prod_{i=1}^{k-1}({\mathbf{m}}_{ii}-{\mathbf{m}}_{ij})=0.$

So, $h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})_{*j}=0$ for $j=1,\dotsc,k-1$.
For $j=k$, we apply (13) to evaluate
$h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})_{*k}$:

$\left(h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})\right)_{k}=\frac{1}{v_{k}}h^{\mathsf{T}}\cdot\left(r-\sum_{j=1}^{k-1}v_{j}\mathbb{H}({\mathbf{m}})_{*j}\right)=\frac{1}{v_{k}}\left(h^{\mathsf{T}}\cdot
r-\sum_{j=1}^{k-1}v_{j}h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})_{*j}\right)=\frac{1}{v_{k}}h^{\mathsf{T}}\cdot
r.$

The norm of $h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})$ is then upper
bounded by

$\displaystyle\left\|h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}})\right\|$
$\displaystyle=\left|(h^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}}))_{k}\right|$
$\displaystyle=\frac{1}{v_{k}}h^{\mathsf{T}}\cdot r$
$\displaystyle\leq\sqrt{k}\left\|h\right\|\left\|r\right\|$
$\displaystyle=\sqrt{k}\left\|h\right\|\sigma_{k}(\mathbb{H}({\mathbf{m}})).$

∎

Proof of Lemma 17. Consider any ${\mathbf{G}}\in\mathbb{H}_{1}$, say
${\mathbf{G}}=\mathbb{H}({\mathbf{g}})$, ${\mathbf{g}}\in\mathbb{R}^{[k-1]}$.
Then
$({\mathbf{G}}^{\mathsf{T}}\cdot\mathbb{H}({\mathbf{m}}))_{j}=\sum_{S}{\mathbf{G}}_{S}\mathbb{H}({\mathbf{m}})_{S,j}=\prod_{1}^{k-1}(1+{\mathbf{g}}_{i}{\mathbf{m}}_{i,j})$.
We also note that
$\|{\mathbf{G}}\|=\sqrt{\prod_{1}^{k-1}(1+{\mathbf{g}}_{i}^{2})}$.

We now show that there is some $j$ such that
$\prod_{i=1}^{k-1}\left|\frac{1+{\mathbf{g}}_{i}{\mathbf{m}}_{ij}}{\sqrt{1+{\mathbf{g}}_{i}^{2}}}\right|$
is large. First, for any $i$ for which ${\mathbf{g}}_{i}\geq\frac{1}{2}$,
there is at most one $j$ s.t. ${\mathbf{m}}_{ij}\leq\zeta$; exclude these
$j$’s. Next, for each $i$ for which ${\mathbf{g}}_{i}<\frac{1}{2}$, there is
at most one $j$ s.t.
$\left|\frac{1}{{\mathbf{g}}_{i}}+{\mathbf{m}}_{ij}\right|\leq\zeta/2$;
exclude these $j$’s. For the remainder of the argument fix any $j$ which has
not been excluded. Since ${\mathbf{m}}$ has $k$ columns while
${\mathbf{g}}\in\mathbb{R}^{k-1}$, such a $j$ exists. We now lower bound
$\left|\frac{1+{\mathbf{g}}_{i}{\mathbf{m}}_{ij}}{\sqrt{1+{\mathbf{g}}_{i}^{2}}}\right|$
for each $i$; there are three cases.

  1. (1)

${\mathbf{g}}_{i}\geq 1/2,{\mathbf{m}}_{ij}>\zeta$. Then
$\left|\frac{1+{\mathbf{g}}_{i}{\mathbf{m}}_{ij}}{\sqrt{1+{\mathbf{g}}_{i}^{2}}}\right|\geq{\mathbf{m}}_{ij}>\zeta$.

  2. (2)

$-1/2<{\mathbf{g}}_{i}<1/2$. Then
$\left|\frac{1+{\mathbf{g}}_{i}{\mathbf{m}}_{ij}}{\sqrt{1+{\mathbf{g}}_{i}^{2}}}\right|>\sqrt{\frac{({\mathbf{m}}_{ij}-2)^{2}}{5}}\geq
1/\sqrt{5}$.

  3. (3)

${\mathbf{g}}_{i}\leq-1/2,|\frac{1}{{\mathbf{g}}_{i}}+{\mathbf{m}}_{ij}|>\zeta/2$.
Then
$\left|\frac{1+{\mathbf{g}}_{i}{\mathbf{m}}_{ij}}{\sqrt{1+{\mathbf{g}}_{i}^{2}}}\right|=\left|\frac{{\mathbf{g}}_{i}(\frac{1}{{\mathbf{g}}_{i}}+{\mathbf{m}}_{ij})}{{\mathbf{g}}_{i}\sqrt{1+1/{\mathbf{g}}_{i}^{2}}}\right|>\frac{\zeta}{2\sqrt{5}}$.

We therefore have $\tau({\mathbf{m}})>(\zeta/2\sqrt{5})^{k-1}$. ∎

Part 1 is an immediate consequence of the two lemmas. ∎

##  5\. Analysis of the algorithm

To begin with, we will need the following result on the numerical accuracy of
computing the SVD:

######  Lemma 18 ([13], section 5.4.1).

Given a matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ and precision
$\varepsilon>0$, the Golub-Kahan-Reinsch algorithm computes an approximate SVD
given by $\mathbf{U}\in\mathbb{R}^{m\times
m},\mathbf{\Sigma}\in\mathbb{R}^{m\times n},\mathbf{V}\in\mathbb{R}^{n\times
n}$ such that

  * •

$\mathbf{\Sigma}$ is a diagonal matrix;

  * •

$\mathbf{U}=\mathbf{W}+\mathbf{\Delta U}$ and
$\mathbf{V}=\mathbf{Z}+\mathbf{\Delta V}$, where $\mathbf{W},\mathbf{Z}$ are
unitary and $||\mathbf{\Delta U}||,||\mathbf{\Delta V}||<\varepsilon$;

  * •

$\mathbf{W}\mathbf{\Sigma}\mathbf{Z}^{\mathsf{T}}=\mathbf{A}+\mathbf{\Delta
A}$ with $||\mathbf{\Delta A}||<\varepsilon||\mathbf{A}||$.

To analyze Algorithm 1, let
${\mathbf{d}}=d_{\mathrm{stat}}(\operatorname{{\tilde{g}}},\operatorname{g})=\max\\{||{\tilde{\mathbf{C}}}_{ST}-\mathbf{C}_{ST}||_{\infty},||{\tilde{\mathbf{C}}}_{ST,1}-\mathbf{C}_{ST,1}||_{\infty}\\}$
as in Definition (11). Furthermore, assume that the SVD in line 3 of the
algorithm is calculated with precision $\varepsilon={\mathbf{d}}$ according to
the statement of Lemma 18. Then, we get the following bounds:

######  Lemma 19.

  1. (a)

$\sigma_{k}({\tilde{\mathbf{C}}}_{ST})\geq{\pi_{\min}}{\bm{\sigma}}^{2}-2^{k}{\mathbf{d}}$

  2. (b)

$\sigma_{k}({\hat{\mathbf{C}}}_{ST})\geq{\pi_{\min}}{\bm{\sigma}}^{2}-2^{k+2}{\mathbf{d}}$

###### Proof.

(a) By a classical perturbation bound of Weyl [33], we have
$|\sigma_{k}({\tilde{\mathbf{C}}}_{ST})-\sigma_{k}(\mathbf{C}_{ST})|\leq||{\tilde{\mathbf{C}}}_{ST}-\mathbf{C}_{ST}||$.
Moreover, we know $||{\tilde{\mathbf{C}}}_{ST}-\mathbf{C}_{ST}||\leq
2^{k}||{\tilde{\mathbf{C}}}_{ST}-\mathbf{C}_{ST}||_{\infty}\leq
2^{k}{\mathbf{d}}$. Combining this with Part 2 of Theorem 14, we get the
desired result.

(b) Let $\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\mathsf{T}}$ be an approximate
SVD for ${\tilde{\mathbf{C}}}_{ST}$ satisfying the conditions of Lemma 18 with
precision $\varepsilon$. In particular, $\mathbf{U}=\mathbf{W}+\mathbf{\Delta
U}$, $\mathbf{V}=\mathbf{Z}+\mathbf{\Delta V}$ and
$\mathbf{W}\mathbf{\Sigma}\mathbf{Z}^{\mathsf{T}}={\tilde{\mathbf{C}}}_{ST}+\mathbf{\Delta}{\tilde{\mathbf{C}}}_{ST}$.
We assume that the columns of $\mathbf{\Sigma}$ are ordered by magnitude of
their diagonal entries. Let
$\hat{\mathbf{U}},\hat{\mathbf{V}},\hat{\mathbf{W}},\hat{\mathbf{Z}},\mathbf{\Delta}\hat{\mathbf{U}},\mathbf{\Delta}\hat{\mathbf{V}}$
denote the first $k$ columns of the respective matrices. We get

$\displaystyle\sigma_{k}({\hat{\mathbf{C}}}_{ST})$
$\displaystyle=\sigma_{k}(\hat{\mathbf{U}}^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\hat{\mathbf{V}})=\sigma_{k}((\hat{\mathbf{W}}+\mathbf{\Delta}\hat{\mathbf{U}})^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}(\hat{\mathbf{Z}}+\mathbf{\Delta}\hat{\mathbf{V}}))$
$\displaystyle\geq\sigma_{k}(\hat{\mathbf{W}}^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\hat{\mathbf{Z}})-||(\mathbf{\Delta}\hat{\mathbf{U}})^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\hat{\mathbf{Z}}||-||\hat{\mathbf{W}}^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\mathbf{\Delta}\hat{\mathbf{V}}||-||(\mathbf{\Delta}\hat{\mathbf{U}})^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\mathbf{\Delta}\hat{\mathbf{V}}||$
$\displaystyle\geq\sigma_{k}(\hat{\mathbf{W}}^{\mathsf{T}}{\tilde{\mathbf{C}}}_{ST}\hat{\mathbf{Z}})-3\varepsilon||{\tilde{\mathbf{C}}}_{ST}||$
$\displaystyle\geq\sigma_{k}(\hat{\mathbf{W}}^{\mathsf{T}}({\tilde{\mathbf{C}}}_{ST}+\mathbf{\Delta}{\tilde{\mathbf{C}}}_{ST})\hat{\mathbf{Z}})-||\hat{\mathbf{W}}^{\mathsf{T}}\mathbf{\Delta}{\tilde{\mathbf{C}}}_{ST}\hat{\mathbf{V}}||-3\varepsilon||{\tilde{\mathbf{C}}}_{ST}||$
$\displaystyle\geq\sigma_{k}(({\tilde{\mathbf{C}}}_{ST}+\mathbf{\Delta}{\tilde{\mathbf{C}}}_{ST}))-4\varepsilon||{\tilde{\mathbf{C}}}_{ST}||$
$\displaystyle\geq\sigma_{k}({\tilde{\mathbf{C}}}_{ST})-5\varepsilon||{\tilde{\mathbf{C}}}_{ST}||.$

Since the entries of ${\tilde{\mathbf{C}}}_{ST}$ are all at most $1$, we have
$||{\tilde{\mathbf{C}}}_{ST}||\leq 2^{k-1}$ and with
$\varepsilon={\mathbf{d}}$, the result follows. ∎

At this point, we can view ${\hat{\mathbf{C}}}_{ST}$ and
${\hat{\mathbf{C}}}_{ST,1}$ as the two slices of a $k\times k\times 2$-tensor
$\hat{{\mathcal{T}}}$ that is close to the tensor
${\mathcal{T}}=[\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]),\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T]),\mathbb{H}({\mathbf{m}}_{1})\cdot\pi_{\odot}]$
(with the two slices
${\mathcal{T}}(:,:,0)=\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST}\hat{\mathbf{V}}$
and
${\mathcal{T}}(:,:,1)=\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST,1}\hat{\mathbf{V}}$).
The following lemma bounds the distance between $\hat{{\mathcal{T}}}$ and
${\mathcal{T}}$, and shows that the first two components of ${\mathcal{T}}$
are full rank and well-conditioned.

######  Lemma 20.

Suppose that ${\mathbf{d}}\leq{\pi_{\min}}{\bm{\sigma}}^{2}/(k2^{2k+2})$, then

  1. (a)

$||{\hat{\mathbf{C}}}_{ST}-\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST}\hat{\mathbf{V}}||_{\infty},||{\hat{\mathbf{C}}}_{ST,1}-\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST,1}\hat{\mathbf{V}}||_{\infty}<2^{2k}{\mathbf{d}}$

  2. (b)

$\kappa(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])),\kappa(\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T]))\leq
k^{2}2^{2k+1}/({\pi_{\min}}{\bm{\sigma}}^{2}).$

###### Proof.

(a) Remember that
$\hat{\mathbf{U}}=\hat{\mathbf{W}}+\mathbf{\Delta}\hat{\mathbf{U}},\hat{\mathbf{V}}=\hat{\mathbf{Z}}+\mathbf{\Delta}\hat{\mathbf{V}}$,
where
$||\mathbf{\Delta}\hat{\mathbf{U}}||,||\mathbf{\Delta}\hat{\mathbf{V}}||<{\mathbf{d}}$
and the columns of $\hat{\mathbf{W}},\hat{\mathbf{Z}}$ are orthonormal. In
particular, this implies that entries of $\hat{\mathbf{U}}$ and
$\hat{\mathbf{V}}$ are bounded by $1+k{\mathbf{d}}\leq 2$. Hence, we have

$\displaystyle||{\hat{\mathbf{C}}}_{ST}-\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST}\hat{\mathbf{V}}||_{\infty}=||\hat{\mathbf{U}}^{\mathsf{T}}({\tilde{\mathbf{C}}}_{ST}-\mathbf{C}_{ST})\hat{\mathbf{V}}||_{\infty}$
$\displaystyle\leq$
$\displaystyle||\hat{\mathbf{U}}^{\mathsf{T}}||_{\infty}\cdot
2^{k-1}\cdot||{\tilde{\mathbf{C}}}_{ST}-\mathbf{C}_{ST}||_{\infty}\cdot||\hat{\mathbf{V}}||_{\infty}\cdot
2^{k-1}\leq 2^{2k}{\mathbf{d}}.$

The result follows analogously for ${\hat{\mathbf{C}}}_{ST,1}$.

(b) First, since the entries of $\mathbb{H}({\mathbf{m}}[S])$ are bounded by
$1$, we have
$\sigma_{1}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]))\leq
k\cdot||\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])||_{\infty}\leq
k2^{k}$, and similarly,
$\sigma_{1}(\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\hat{\mathbf{V}})\leq
k2^{k}$. Using part (a), Lemma 19, and the assumption on ${\mathbf{d}}$, we
get

$\displaystyle\sigma_{k}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST}\hat{\mathbf{V}})$
$\displaystyle\geq\sigma_{k}({\hat{\mathbf{C}}}_{ST})-||{\hat{\mathbf{C}}}_{ST}-\hat{\mathbf{U}}^{\mathsf{T}}\mathbf{C}_{ST}\hat{\mathbf{V}}||\geq\sigma_{k}({\hat{\mathbf{C}}}_{ST})-k2^{2k}{\mathbf{d}}$
$\displaystyle\geq{\pi_{\min}}{\bm{\sigma}}^{2}-2^{k+2}{\mathbf{d}}-k2^{2k}{\mathbf{d}}\geq{\pi_{\min}}{\bm{\sigma}}^{2}/2$

Hence, we deduce

$\displaystyle{\pi_{\min}}{\bm{\sigma}}^{2}/2$
$\displaystyle\leq\sigma_{k}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])\pi_{\odot}\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\hat{\mathbf{V}})$
$\displaystyle\leq\sigma_{k}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]))\cdot\sigma_{1}(\pi_{\odot})\cdot\sigma_{1}(\mathbb{H}({\mathbf{m}}[T])^{\mathsf{T}}\hat{\mathbf{V}})$
$\displaystyle\leq\sigma_{k}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]))\cdot^{\cdot}k2^{k}.$

We conclude that
$\kappa(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]))=\sigma_{1}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]))/\sigma_{k}(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]))\leq
2\cdot(k2^{k})^{2}/({\pi_{\min}}{\bm{\sigma}}^{2})$, and the result follows
analogously for
$\kappa(\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T]))$. ∎

The following result from [2] provides us with error bounds for the core step
of the tensor decomposition algorithm we are using.

######  Theorem 21 ([2], Theorem 2.3).

Let $\varepsilon>0$ and ${\mathcal{T}},\hat{{\mathcal{T}}}$ be two $k\times
k\times 2$-tensors, such that

  * •

${\mathcal{T}}=[\mathbf{X},\mathbf{Y},\mathbf{Z}]$ with
$\mathbf{X},\mathbf{Y}\in\mathbb{R}^{k\times
k},\mathbf{Z}\in\mathbb{R}^{2\times k}$;

  * •

$\kappa(\mathbf{X}),\kappa(\mathbf{Y})\leq\kappa$;

  * •

the entries of $(\mathbf{Z}_{1k}\mathbf{Z}_{2k}^{-1})_{k}$ are
$\zeta$-separated;

  * •

$||\mathbf{X}_{i}||_{2},||\mathbf{Y}_{i}||_{2},||\mathbf{Z}_{i}||_{2}$ are
bounded by a constant;

  * •

$||\hat{{\mathcal{T}}}-{\mathcal{T}}||_{\infty}<\varepsilon\cdot\mbox{poly}(1/\kappa,1/k,\zeta)$;

then the eigenvectors of ${\mathcal{T}}(:,:,1){\mathcal{T}}(:,:,0)^{-1}$ and
${\mathcal{T}}(:,:,1)^{\mathsf{T}}({\mathcal{T}}(:,:,0)^{\mathsf{T}})^{-1}$
approximate the columns $\mathbf{X}_{i}$ and $\mathbf{Y}_{i}$ respectively, up
to permutation and additive error $\varepsilon$.

(_Comment:_ This is slightly more specific than the theorem statement in the
reference, but it easily follows from the general statement.)

Finally, we will make use of the following classical result on perturbations
of linear systems:

######  Lemma 22 ([20], Section 7.1).

Let $\mathbf{A}x=b$ and $({\mathbf{A}}+\mathbf{\Delta}{\mathbf{A}})y=b+\Delta
b$, where $||\mathbf{\Delta}{\mathbf{A}}||\leq\gamma||{\mathbf{A}}||$ and
$||\Delta b||\leq\gamma||b||$, and assume that
$\gamma\cdot\kappa({\mathbf{A}})<1$. Then,

$\frac{||x-y||_{2}}{||x||_{2}}\leq\frac{2\gamma\cdot\kappa({\mathbf{A}})}{1-\gamma\cdot\kappa({\mathbf{A}})}.$

Now, we have all the tools to prove Theorem 11.

Proof of Theorem 11. Fix $\varepsilon\in(0,\zeta/2)$ and a model
$(\pi,{\mathbf{m}})\in\mathcal{D}_{2k-1,\zeta,{\pi_{\min}}}$ with statistics
$\operatorname{g}$, and suppose that Algorithm 1 is given approximate
statistics $\operatorname{{\tilde{g}}}$ with
$||\operatorname{g}-\operatorname{{\tilde{g}}}||_{\infty}={\mathbf{d}}\leq\varepsilon\cdot({\pi_{\min}}\zeta^{k})^{C}$
for some large constant $C$. Consider the $k\times k\times 2$-tensor
$\hat{{\mathcal{T}}}$ that consists of the slices ${\hat{\mathbf{C}}}_{ST}$
and ${\hat{\mathbf{C}}}_{ST,1}$ and
${\mathcal{T}}=\left[\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S]),\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T]),\mathbb{H}({\mathbf{m}}_{1})\pi_{\odot}\right]$.
Define
$\kappa=\max\\{\kappa(\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])),\kappa(\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T]))\\}$.
By Lemma 20, we have
$\kappa\leq\frac{1}{{\pi_{\min}}}\cdot\left(\frac{1}{\zeta}\right)^{O(k)}$ and
$||\hat{{\mathcal{T}}}-{\mathcal{T}}||_{\infty}\leq\exp(O(k)){\mathbf{d}}$.
Hence, by Theorem 21, we get
$||\hat{\mathbf{S}}-\hat{\mathbf{U}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[S])||_{\infty},||\hat{\mathbf{T}}-\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])||_{\infty}\leq\exp(O(k)){\mathbf{d}}\cdot\mbox{poly}(\kappa,k,1/\zeta)=\left(\frac{1}{{\pi_{\min}}}\right)^{O(1)}\left(\frac{1}{\zeta}\right)^{O(k)}{\mathbf{d}}$,
after possibly permuting the columns of $\hat{\mathbf{S}}$ and
$\hat{\mathbf{T}}$. Now, $\tilde{\pi}$ is defined via a linear system that is
a perturbation of equation (8). If ${\mathbf{d}}$ is small enough, then the
conditions of Lemma 22 are satisfied with $\gamma=\varepsilon/(4\kappa)$, and
we get $||\pi-\tilde{\pi}||_{2}\leq\varepsilon$. Similarly, we get
$\tilde{{\mathbf{m}}_{i}}$ from a perturbed version of equation (9) and we can
use Lemma 22 with $\gamma=\varepsilon/(4\kappa\cdot{\pi_{\min}})$ to deduce
$||{\mathbf{m}}_{i}-\tilde{{\mathbf{m}}_{i}}||_{2}\leq\varepsilon$. Changing
the norms to $||.||_{\infty}$ incurs at most another factor of $k$, by which
we can decrease ${\mathbf{d}}$, so then the output of Algorithm 1 satisfies
$d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\tilde{\pi},\tilde{{\mathbf{m}}}))<\varepsilon$.
To prove the second part of the theorem, suppose the model
$(\tilde{\pi},\tilde{{\mathbf{m}}})$ has statistics
$\operatorname{{\tilde{g}}}$, and again
$||\operatorname{g}-\operatorname{{\tilde{g}}}||_{\infty}\leq\varepsilon\cdot({\pi_{\min}}\zeta^{k})^{C}$.
By Lemma 19, ${\tilde{\mathbf{C}}}_{ST}$ has rank $k$. Hence, by equation (6),
$\tilde{{\mathbf{m}}}_{1}$ is the vector of eigenvalues of
${\hat{\mathbf{C}}}_{ST,1}{\hat{\mathbf{C}}}_{ST}^{-1}$. At the same time, by
the first part of the theorem, these eigenvalues give an
$\varepsilon$-approximation to ${\mathbf{m}}$ (after permuting them), and
$\varepsilon<\zeta/2$ implies that they must be separated. We essentially
proved in section 3.1 through equations (6) - (9) that
${\tilde{\mathbf{C}}}_{ST}$ having rank $k$ and separation of
$\tilde{{\mathbf{m}}}_{1}$ are sufficient conditions for Algorithm 1 to
perfectly recover $(\tilde{\pi},\tilde{{\mathbf{m}}})$ given perfect
statistics
$\operatorname{{\tilde{g}}}=\gamma(\tilde{\pi},\tilde{{\mathbf{m}}})$. But
then, the first part of the theorem implies
$d_{\mathrm{model}}((\tilde{\pi},\tilde{{\mathbf{m}}}),(\pi,{\mathbf{m}}))\leq\varepsilon$,
as desired. ∎

Proof of Corollary 13. Suppose we have $n\geq 2k-1$ variables and we know a
subset $A\subseteq[n]$ of $2k-1$ $\zeta$-separated variables. Then, Algorithm
1 can be used to identify ${\mathbf{m}}[A]$ and $\pi$ up to error
$\varepsilon$ in runtime $\exp(O(k))$ using
$(1/\zeta)^{O(k)}(1/{\pi_{\min}})^{O(1)}(1/\varepsilon)^{2}$ many samples. For
each $i\notin A$, we can then compute

$\displaystyle\tilde{{\mathbf{m}}_{i}}=\left((\operatorname{{\tilde{g}}}(R\cup\left\\{i\right\\}))_{R\subseteq
T}\right)^{\mathsf{T}}\cdot\hat{\mathbf{V}}\cdot{\left(\hat{\mathbf{T}}^{\mathsf{T}}\right)}^{-1}\cdot\tilde{\pi}_{\odot}^{-1}$

as in line 9 of the algorithm (where $T\subseteq A$ is a set of size $k-1$ and
$\hat{\mathbf{T}}$ approximates
$\hat{\mathbf{V}}^{\mathsf{T}}\mathbb{H}({\mathbf{m}}[T])$). This takes
runtime at most $n\cdot\exp(O(k))$. Given that
$||(\operatorname{{\tilde{g}}}(R\cup\left\\{i\right\\}))_{R\subseteq
T}-(\operatorname{g}(R\cup\left\\{i\right\\}))_{R\subseteq
T}||_{\infty}\leq\varepsilon({\pi_{\min}}\zeta^{k})^{C}$ for some large enough
$C$, the same analysis as for Theorem 11 asserts that
$\tilde{{\mathbf{m}}_{i}}$ is an $\varepsilon$-approximation to
${\mathbf{m}}_{i}$ (here, $\operatorname{g}$ is the vector of perfect
statistics). Hence, computing $\varepsilon$-approximations to
${\mathbf{m}}_{i}$ for all $i\notin A$ requires obtaining
$(\varepsilon({\pi_{\min}}\zeta^{k})^{C})$-approximations to $(n-|A|)\cdot
2^{k-1}$ entries of the observable moment vector $\operatorname{g}$. Standard
Chernoff bounds and a union bound show that this is possible with $\log
n\cdot(1/\zeta)^{O(k)}(1/{\pi_{\min}})^{O(1)}(1/\varepsilon)^{2}$ many
samples.
If the subset of $2k-1$ $\zeta$-separated variables is not known (but
guaranteed to exist), we can simply run Algorithm 1 to identify
${\mathbf{m}}[A]$ for all $\binom{n}{2k-1}$ possible guesses $A$ of this
subset (stopping the algorithm when encountering any issues that might occur
when ${\hat{\mathbf{C}}}_{ST}$ is not invertible). Then, we compute the vector
of statistics
$\gamma(\tilde{\pi},\tilde{{\mathbf{m}}}[A])=\mathbb{H}(\tilde{{\mathbf{m}}}[A])\tilde{\pi}$
for each valid output $(\tilde{\pi},\tilde{{\mathbf{m}}}[A])$ and choose the
model $(\tilde{\pi},\tilde{{\mathbf{m}}}[A])$ that minimizes the distance
${\mathbf{d}}(A)=||\gamma(\tilde{\pi},\tilde{{\mathbf{m}}}[A])-\operatorname{{\tilde{g}}}[2^{A}]||_{\infty}$.
All this takes runtime at most $n^{2k}\exp(O(k))$. Let $A^{*}$ be the correct
guess and suppose we have
$||\operatorname{{\tilde{g}}}[2^{A^{*}}]-\operatorname{g}[2^{A^{*}}]||_{\infty}\leq\delta$.
Then, by Theorem 11, we get
$d_{\mathrm{model}}((\tilde{\pi},\tilde{{\mathbf{m}}}[A^{*}]),(\pi,{\mathbf{m}}[A^{*}]))\leq\delta\cdot({\pi_{\min}}\zeta^{k})^{-C}$
for some large, positive constant $C$. Hence, we have (after possibly
permuting $\tilde{\pi},\tilde{{\mathbf{m}}}[A^{*}]$) that
$||\pi-\tilde{\pi}||_{\infty}\leq\delta\cdot({\pi_{\min}}\zeta^{k})^{-C}$ and
$||\mathbb{H}(\tilde{{\mathbf{m}}}[A^{*}])-\mathbb{H}({\mathbf{m}}[A^{*}])||_{\infty}\leq
2^{2k}\delta\cdot({\pi_{\min}}\zeta^{k})^{-C}$. This implies

$\displaystyle{\mathbf{d}}(A^{*})$
$\displaystyle=||\gamma(\tilde{\pi},\tilde{{\mathbf{m}}}[A^{*}])-\operatorname{{\tilde{g}}}[2^{A^{*}}]||_{\infty}$
$\displaystyle\leq||\mathbb{H}(\tilde{{\mathbf{m}}}[A^{*}])\tilde{\pi}-\mathbb{H}({\mathbf{m}}[A^{*}])\pi||_{\infty}+||\operatorname{g}[2^{A^{*}}]-\operatorname{{\tilde{g}}}[2^{A^{*}}]||_{\infty}$
$\displaystyle\leq||\mathbb{H}(\tilde{{\mathbf{m}}}[A^{*}])\tilde{\pi}-\mathbb{H}({\mathbf{m}}[A^{*}])\tilde{\pi}||_{2}+||\mathbb{H}({\mathbf{m}}[A^{*}])\tilde{\pi}-\mathbb{H}({\mathbf{m}}[A^{*}])\pi||_{2}+\delta$
$\displaystyle\leq||\mathbb{H}(\tilde{{\mathbf{m}}}[A^{*}])-\mathbb{H}({\mathbf{m}}[A^{*}])||_{2}\cdot||\tilde{\pi}||_{2}+||\mathbb{H}({\mathbf{m}}[A^{*}])||_{2}\cdot||\tilde{\pi}-\pi||_{2}+\delta$
$\displaystyle\leq
2^{k}||\mathbb{H}(\tilde{{\mathbf{m}}}[A^{*}])-\mathbb{H}({\mathbf{m}}[A^{*}])||_{\infty}\cdot\sqrt{k}+2^{k}\cdot\sqrt{k}||\tilde{\pi}-\pi||_{\infty}+\delta$
$\displaystyle\leq 2^{4k}\delta\cdot({\pi_{\min}}\zeta^{k})^{-C}.$

If $\delta<2^{-4k}({\pi_{\min}}\zeta^{k})^{2C}\varepsilon$, then the minimal
distance ${\mathbf{d}}(A)$ is at most
$\varepsilon({\pi_{\min}}\zeta^{k})^{C}$. Hence, by Theorem 11, the model
$(\tilde{\pi},\tilde{{\mathbf{m}}}[A])$ is an $\varepsilon$-approximation of
$(\pi,{\mathbf{m}}[A])$ and we can proceed as in the first part of the proof
to get a full $\varepsilon$-approximation to $(\pi,{\mathbf{m}})$. Ensuring
that
$||\operatorname{{\tilde{g}}}[2^{A}]-\operatorname{g}[2^{A}]||_{\infty}\leq
2^{-4k}({\pi_{\min}}\zeta^{k})^{2C}\varepsilon$ for all subsets of
$A\subseteq[n]$ of size $2k-1$ requires $\log(n^{2k})\cdot
2^{8k}({\pi_{\min}}\zeta^{k})^{-4C}\varepsilon^{-2}=\log
n\cdot(1/{\pi_{\min}})^{O(1)}(1/\zeta)^{O(k)}(1/\varepsilon)^{2}$ many
samples. ∎

##  6\. Lower bounds

The following theorem shows that our algorithmic results are optimal when
$\zeta$ is small enough.

######  Theorem 23.

Let $n=2k-1,\zeta\leq\frac{1}{8k},{\pi_{\min}}\leq\frac{1}{4k},\varepsilon>0$
and $\varepsilon<\min\\{\frac{{\pi_{\min}}}{4\sqrt{k}},\zeta\\}$. Then, there
exist models $(\pi,{\mathbf{m}})$,
$(\pi^{\prime},{\mathbf{m}}^{\prime})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}}$
such that
$d_{\mathrm{stat}}(\gamma(\pi,{\mathbf{m}}),\gamma(\pi^{\prime},{\mathbf{m}}^{\prime}))\leq(k\zeta)^{\Omega(k)}\varepsilon$,
but
$d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\pi^{\prime},{\mathbf{m}}^{\prime}))>\varepsilon$.

######  Corollary 24.

For $n=2k-1$ random variables $X_{i}$, sample complexity

$(1/(k\zeta))^{\Omega(k)}(1/\varepsilon)$

is necessary to compute a model that w.h.p. satisfies (10).

Note that (assuming ${\pi_{\min}}$ is not smaller than $\zeta^{O(k)}$), the
upper and lower sample complexity bounds in Corollary 12 and Corollary 24
match in the case that $\zeta\leq k^{-1-\delta}$ for some arbitrary small
$\delta>0$. Only when $\zeta$ comes closer to its maximal possible value of
$\frac{1}{k-1}$, there is a gap between upper and lower bound. In the edge
case, when $\zeta=\Theta(\frac{1}{k})$, the upper bound evaluates to
$\exp{(O(k\log k))}$ and the lower bound evaluates to $\exp{(\Omega(k))}$. We
remark that the lower bound can be shown to be tight for $k$-MixIID over the
entire range of the parameter $\zeta$.

To prove Theorem 23, we start with the following lemma:

######  Lemma 25.

Let $n\geq 1,\varepsilon>0$ and
$\varepsilon<\min\\{\frac{{\pi_{\min}}}{4\sqrt{k}},\zeta\\}$. Suppose
$(\pi,{\mathbf{m}})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}}$ and
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))=\sigma<\frac{1}{2}$. Then, there exists
$\hat{\pi}$ with $\min_{j}\hat{\pi}_{j}\geq\frac{1}{4}{\pi_{\min}}$ and such
that
$d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\hat{\pi},{\mathbf{m}}))>\varepsilon$
but
$d_{\mathrm{stat}}(\gamma(\pi,{\mathbf{m}}),\gamma(\hat{\pi},{\mathbf{m}}))\leq
4k\sigma\cdot\varepsilon$.

###### Proof.

Let $(\pi,{\mathbf{m}})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}}$ and
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))=\sigma<\frac{1}{2}$. Let
$\tilde{\mathbb{H}({\mathbf{m}})}$ be the best rank-$(k-1)$-approximation of
$\mathbb{H}({\mathbf{m}})$. By the Eckart-Young Theorem [9], we have
$||\tilde{\mathbb{H}({\mathbf{m}})}-\mathbb{H}({\mathbf{m}})||=\sigma$. Let
$\alpha\in\mathbb{R}^{k}$ be a vector in the right kernel of
$\tilde{\mathbb{H}({\mathbf{m}})}$ with $||\alpha||_{2}=1$. Let $\mathbb{1}$
denote the all-ones vector in $\mathbb{R}^{k}$ and let $e_{1}$ denote the
vector whose first entry is 1 and all other entries are zero. We have

$\displaystyle|\mathbb{1}^{\mathsf{T}}\alpha|=|(\mathbb{H}({\mathbf{m}})\alpha)_{1}|=|((\mathbb{H}({\mathbf{m}})-\tilde{\mathbb{H}({\mathbf{m}})})\alpha)_{1}|\leq||\mathbb{H}({\mathbf{m}})-\tilde{\mathbb{H}({\mathbf{m}})}||\cdot||\alpha||_{2}\leq\sigma.$

Now, define
$\hat{\pi}=\pi+2\sqrt{k}\varepsilon\cdot(\alpha-(\mathbb{1}^{\mathsf{T}}\alpha)e_{1})$.
First, we check that $\hat{\pi}$ is a valid probability vector. By our
assumptions, we have

$\displaystyle||2\sqrt{k}\varepsilon\cdot(\alpha-(\mathbb{1}^{\mathsf{T}}\alpha)e_{1})||_{\infty}\leq
2\sqrt{k}\varepsilon\cdot(||\alpha||_{\infty}+|\mathbb{1}^{\mathsf{T}}\alpha|)<\frac{{\pi_{\min}}}{2}\cdot(1+\sigma)\leq\frac{3}{4}{\pi_{\min}},$

hence, all the entries of $\hat{\pi}$ are larger than
$\frac{1}{4}{\pi_{\min}}$. Moreover, we have

$\displaystyle\sum_{j=1}^{k}\hat{\pi}_{j}=\mathbb{1}^{\mathsf{T}}\pi+2\sqrt{k}\varepsilon\cdot(\mathbb{1}{{}^{\mathsf{T}}}\alpha-\mathbb{1}^{\mathsf{T}}\alpha)=0.$

Now, recall the definition of $d_{\mathrm{model}}$ as

$\displaystyle
d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\hat{\pi},{\mathbf{m}}))=\min_{\rho\in
S_{k}}\max\\{\max_{j}|\pi_{j}-\hat{\pi}_{\rho(j)}|,\max_{i,j}|{\mathbf{m}}_{i,j}-{\mathbf{m}}_{i,\rho(j)}|\\}.$

For any permutation $\rho$ that is not the identity, we have
$\max_{i,j}|{\mathbf{m}}_{i,j}-{\mathbf{m}}_{i,\rho(j)}|\geq\zeta>\varepsilon$
by $\zeta$-separation of ${\mathbf{m}}$. For the identity permutation, we have

$\displaystyle\max_{j}|\pi_{j}-\hat{\pi}_{j}|=||2\sqrt{k}\varepsilon\cdot(\alpha-(\mathbb{1}^{\mathsf{T}}\alpha)e_{1})||_{\infty}\geq
2\varepsilon\cdot(||\alpha||_{2}-||(\mathbb{1}^{\mathsf{T}}\alpha)e_{1}||_{2})\geq
2\varepsilon\cdot(1-\sigma)>\varepsilon,$

so we conclude
$d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\hat{\pi},{\mathbf{m}}))>\varepsilon$.
Moreover, we have

$\displaystyle
d_{\mathrm{stat}}(\gamma(\pi,{\mathbf{m}}),\gamma(\hat{\pi},{\mathbf{m}}))$
$\displaystyle=||\mathbb{H}({\mathbf{m}})\pi-\mathbb{H}({\mathbf{m}})\hat{\pi}||_{\infty}$
$\displaystyle=2\sqrt{k}\varepsilon\cdot||\mathbb{H}({\mathbf{m}})(\alpha-(\mathbb{1}^{\mathsf{T}}\alpha)e_{1})||_{\infty}$
$\displaystyle\leq
2\sqrt{k}\varepsilon\cdot\left(||\mathbb{H}({\mathbf{m}})\alpha||_{\infty}+|\mathbb{1}^{\mathsf{T}}\alpha|\cdot||\mathbb{H}({\mathbf{m}})||_{\infty}\right)$
$\displaystyle\leq
2\sqrt{k}\varepsilon\cdot\left(\sqrt{k}\cdot||(\mathbb{H}({\mathbf{m}})-\tilde{\mathbb{H}({\mathbf{m}})})\alpha||_{2}+\sigma\right)$
$\displaystyle\leq
2\sqrt{k}\varepsilon\cdot\left(\sqrt{k}\sigma+\sigma\right)\leq
4k\sigma\cdot\varepsilon.$

∎

Lemma 25 reduces the challenge of finding models with large model distance but
small statistical distance to finding a model with a small $k$’th singular
value of the Hadamard extension $\mathbb{H}({\mathbf{m}})$. To do this, we
will need the following known result on Vandermonde matrices (also see
Definition 5).

######  Lemma 26 ([20], section 22.1).

For any row vector $m\in\mathbb{R}^{k}$, we have

$\displaystyle||\operatorname{Vdm}(m)^{-1}||_{\infty}\geq\max_{i}\prod_{j\neq
i}\frac{\max\\{1,|m_{j}|\\}}{|m_{i}-m_{j}|}.$

######  Lemma 27.

Let $n\geq k-1$. For any $\zeta\leq\frac{1}{k}$, there exists a matrix
${\mathbf{m}}\in[0,1]^{n\times k}$ with $\zeta$-separated columns and
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))\leq n2^{n}\cdot(k\zeta)^{k}$.

###### Proof.

Given $\zeta$, define ${\mathbf{m}}$ as the matrix with $n$ identical rows of
the form ${\mathbf{m}}_{i}=(0,\zeta,2\zeta,\dots,(k-1)\cdot\zeta)$. According
to Lemma 26, we have

$\displaystyle\sigma_{k}(\operatorname{Vdm}({\mathbf{m}}_{1}))=||\operatorname{Vdm}({\mathbf{m}}_{1})^{-1}||^{-1}\leq
k||\operatorname{Vdm}({\mathbf{m}}_{1})^{-1}||_{\infty}^{-1}\leq
k\cdot(k-1)!\zeta^{k-1}=k!\zeta^{k-1}.$

Now, consider the matrix
$\operatorname{Vdm}({\mathbf{m}}_{1},n)=\binom{\operatorname{Vdm}({\mathbf{m}}_{1})}{\mathbf{R}}$,
where $\mathbf{R}\in\mathbb{R}^{(n-k)\times k}$ denotes the last $n-k$ columns
of $\operatorname{Vdm}({\mathbf{m}}_{1},n)$. Note that $||\mathbf{R}||\leq
n\cdot||\mathbf{R}||_{\infty}\leq n\cdot((k-1)\cdot\zeta)^{k}$. Let
$q\in\mathbb{R}^{k},||q||_{2}=1$ such that
$||\operatorname{Vdm}({\mathbf{m}}_{1})q||_{2}$ is minimized. We have

$\displaystyle||\operatorname{Vdm}({\mathbf{m}}_{1},n)q||_{2}\leq||\operatorname{Vdm}({\mathbf{m}}_{1})q||_{2}+||\mathbf{R}q||_{2}\leq\sigma_{k}(\operatorname{Vdm}({\mathbf{m}}_{1}))+||\mathbf{R}||\leq
n\cdot k^{k}\cdot\zeta^{k},$

hence, $\sigma_{k}(\operatorname{Vdm}({\mathbf{m}}_{1},n))\leq n\cdot
k^{k}\cdot\zeta^{k}$. Finally, note that $\mathbb{H}({\mathbf{m}})$ is a
matrix with $2^{n}$ rows that are all duplicates of some row in
$\operatorname{Vdm}({\mathbf{m}}_{1},n)$. Hence, we have
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))\leq n2^{n}\cdot(k\zeta)^{k}$. ∎

Proof of Theorem 23. Fix $n,\zeta,{\pi_{\min}}$, and $\varepsilon$ as in the
statement of Theorem 23. Let $\pi:=(1/k,\dots,1/k)^{\mathsf{T}}$. According to
Lemma 27, there exists ${\mathbf{m}}$ such that
$(\pi,{\mathbf{m}})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}}$ and
$\sigma_{k}(\mathbb{H}({\mathbf{m}}))\leq
n2^{n}\cdot(k\zeta)^{k}=(2k-1)2^{2k-1}\cdot(k\zeta)^{k}<\frac{1}{2}$. Now, by
Lemma 25, there exists $\hat{\pi}$ such that
$(\hat{\pi},{\mathbf{m}})\in\mathcal{D}_{n,\zeta,{\pi_{\min}}},d_{\mathrm{model}}((\pi,{\mathbf{m}}),(\hat{\pi},{\mathbf{m}}))>\varepsilon$
and
$d_{\mathrm{stat}}(\gamma(\pi,{\mathbf{m}}),\gamma(\hat{\pi},{\mathbf{m}}))\leq
4k\sigma_{k}(\mathbb{H}({\mathbf{m}}))\cdot\varepsilon\leq(k\zeta)^{\Omega(k)}\cdot\varepsilon$.
∎

##  7\. Discussion

Two larger questions remain to be addressed in this area. First, it would be
of great interest to achieve a similar sample complexity as in Theorem 11 for
the more general “learning” task. Second, it is open to characterize the set
of models, for which identification (even with perfect statistics) is
possible. There exist identifiable models $(\pi,{\mathbf{m}})$ where none of
the rows of ${\mathbf{m}}$ has fully-separated entries. Recent work [17] gives
a sufficient condition for identification that is less restrictive (though
more complicated) than $\zeta$-separation. However, there is no known way of
obtaining in that less restrictive framework the quantitative bounds on noise-
stability, which are essential to this paper. (Quantification in that
framework would likely be misguided, anyway, given that it is a complex yet
not tight characterization; for example it excludes mixtures of subcubes [5].)
Note that some kind of separation assumption is unavoidable if we insist on
$\mathcal{L}_{\infty}$-reconstruction of the model parameters: as one example,
if there are $j,j^{\prime}$ s.t.
${\mathbf{m}}_{ij}={\mathbf{m}}_{ij^{\prime}}$ for all $i$, then it is
impossible to determine $\pi_{j}$ and $\pi_{j^{\prime}}$. However, it might be
possible to completely eliminate the separation assumption in favor of
settling for reconstruction in transportation (Wasserstein) distance. This was
achieved for $k$-MixIID (the $k$-MixProd problem where all observables are
conditionally identically distributed) in [26], and improved in [10]. It is an
open question whether these ideas can be extended to $k$-MixProd, with the
goal being transportation-cost reconstruction of each of the rows of
${\mathbf{m}}$.

We remark that even though Theorem 11 and its corollaries are given for
$\zeta$-separated models, our algorithm itself also works under weaker
conditions. In fact, all it needs is just one $\zeta$-separated observable and
two more disjoint sets of observables, each of which have a Hadamard extension
(see below) with good condition number. These requirements can, for models in
“general position,” be met with as few as $2\lg k+1$ observables. The sample
complexity and runtime of the algorithm scale singly-exponentially in the
number of observables actually used. Hence, it is possible that, except for a
small set of “adversarial” models, model identification typically be achieved
with much lower sample complexity and runtime than can be guaranteed for the
worst case. In fact, this line of thought has already been pursued for the
tensor decomposition problem to which we reduce $k$-MixProd. It was shown that
the time complexity of the tensor decomposition problem significantly improves
when each tensor component of the input tensor is perturbed by small random
noise [2]. Note that one can never solve the $k$-MixProd Identification
problem with fewer than $\lg k$ observables, because then the Hadamard
extensions cannot be full rank, and then the mapping from model to statistics,
cannot be injective.

## References

  * [1] J. Banks, J. Garza-Vargas, A. Kulkarni, and N. Srivastava.  Pseudospectral shattering, the sign function, and diagonalization in nearly matrix multiplication time.  In 2020 IEEE 61st Annual Symposium on Foundations of Computer Science (FOCS), pages 529–540, 2020.  doi:10.1109/FOCS46700.2020.00056.
  * [2] A. Bhaskara, M. Charikar, A. Moitra, and A. Vijayaraghavan.  Smoothed analysis of tensor decompositions.  In Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing, STOC ’14, page 594–603, New York, NY, USA, 2014. Association for Computing Machinery.  doi:10.1145/2591796.2591881.
  * [3] W. R. Blischke.  Estimating the parameters of mixtures of binomial distributions.  Journal of the American Statistical Association, 59(306):510–528, 1964.  doi:10.1080/01621459.1964.10482176.
  * [4] K. Chaudhuri and S. Rao.  Learning mixtures of product distributions using correlations and independence.  In Proc. 21st Ann. Conf. on Learning Theory - COLT, pages 9–20. Omnipress, 2008.  URL: http://colt2008.cs.helsinki.fi/papers/7-Chaudhuri.pdf.
  * [5] 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.  doi:10.1145/3313276.3316375.
  * [6] 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.  doi:10.1137/S0097539798342496.
  * [7] R. de Prony.  Essai expérimentale et analytique.  J. Écol. Polytech., 1(2):24–76, 1795.
  * [8] J. A. Rhodes E. S. Allman, C. Matias.  Identifiability of parameters in latent structure models with many observed variables.  Ann. Statist., 37(6A):3099–3132, 2009.  doi:10.1214/09-AOS689.
  * [9] C. Eckart and G. Young.  The approximation of one matrix by another of lower rank.  Psychometrika, 1:211–218, 1936.  doi:10.1007/BF02288367.
  * [10] Z. Fan and J. Li.  Efficient algorithms for sparse moment problems without separation, 2022\.  arXiv:2207.13008.
  * [11] 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.  doi:10.1137/060670705.
  * [12] Y. Freund and Y. Mansour.  Estimating a mixture of two product distributions.  In Proc. 12th Ann. Conf. on Computational Learning Theory, pages 53–62, July 1999.  doi:10.1145/307400.307412.
  * [13] G. H. Golub and C. F. Van Loan.  Matrix Computations.  The Johns Hopkins University Press, 4th edition, 2013.
  * [14] S. Gordon, B. Mazaheri, L. J. Schulman, and Y. Rabani.  The sparse Hausdorff moment problem, with application to topic models, 2020.  arXiv:2007.08101.
  * [15] S. L. Gordon, B. Mazaheri, Y. Rabani, and L. J. Schulman.  Source identification for mixtures of product distributions.  In Proc. 34th Ann. Conf. on Learning Theory - COLT, volume 134 of Proc. Machine Learning Research, pages 2193–2216. PMLR, 2021.  URL: http://proceedings.mlr.press/v134/gordon21a.html.
  * [16] S. L. Gordon, B. Mazaheri, Y. Rabani, and L. J. Schulman.  Causal inference despite limited global confounding via mixture models.  In Proc. CLeaR, 2023.  URL: www.cclear.cc/2023/AcceptedPapers.
  * [17] S. L. Gordon and L. J. Schulman.  Hadamard extensions and the identification of mixtures of product distributions.  IEEE Transactions on Information Theory, 68(6):4085–4089, 2022\.  doi:10.1109/TIT.2022.3146630.
  * [18] N. Goyal, S. Vempala, and Y. Xiao.  Fourier PCA and robust tensor decomposition.  In Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing, STOC ’14, page 584–593, New York, NY, USA, 2014. Association for Computing Machinery.  doi:10.1145/2591796.2591875.
  * [19] P. Hall and X.-H. Zhou.  Nonparametric estimation of component distributions in a multivariate mixture.  The Annals of Statistics, 31(1):201 – 224, 2003.  doi:10.1214/aos/1046294462.
  * [20] N. J. Higham.  Accuracy and Stability of Numerical Algorithms.  Society for Industrial and Applied Mathematics, second edition, 2002.  doi:10.1137/1.9780898718027.
  * [21] Y. Hua and T. K. Sarkar.  Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise.  IEEE Transactions on Acoustics, Speech, and Signal Processing, 38(5):814–824, 1990.  doi:10.1109/29.56027.
  * [22] 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.  doi:10.1145/195058.195155.
  * [23] 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.  doi:10.1007/978-3-030-17083-7_9.
  * [24] J. B. Kruskal.  Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics.  Linear Algebra and its Applications, 18(2):95–138, 1977.  doi:10.1016/0024-3795(77)90069-6.
  * [25] S. E. Leurgans, R. T. Ross, and R. B. Abel.  A decomposition for three-way arrays.  SIAM Journal on Matrix Analysis and Applications, 14(4):1064–1083, 1993.  doi:10.1137/0614071.
  * [26] 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.  doi:10.1145/2746539.2746584.
  * [27] J. Pearl.  Bayesian networks: A model of self-activated memory for evidential reasoning.  Technical Report CSD-850021, R-43, UCLA Computer Science Department, June 1985.
  * [28] J. Pearl.  Causality.  Cambridge, 2nd edition, 2009.
  * [29] 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.  doi:10.1145/2554797.2554818.
  * [30] P. Spirtes, C. Glymour, and R. Scheines.  Causation, Prediction and Search.  MIT Press, second edition, 2000.
  * [31] B. Tahmasebi, S. A. Motahari, and M. A. Maddah-Ali.  On the identifiability of finite mixtures of finite product measures.  (Also in “On the identifiability of parameters in the population stratification problem: A worst-case analysis,” Proc. ISIT’18 pp. 1051-1055.), 2018.  URL: https://arxiv.org/abs/1807.05444.
  * [32] H. Teicher.  Identifiability of mixtures.  Annals of Mathematical Statistics, 32:244–248, 1961.
  * [33] H. Weyl.  Das asymptotische verteilungsgesetz der eigenwerte linearer partieller differentialgleichungen (mit einer anwendung auf die theorie der hohlraumstrahlung.  Mathematische Annalen, 71(4):441–479, 1912.  doi:10.1007/BF01456804.