\section{Introduction}
A Bayesian Network is a directed acyclic graph $G=(\VV,E)$, on a set of $n$ random variables (identified with the vertices); a Bayesian Network Distribution (BND) is
a probability distribution on the rv's that is Markovian on the graph. That is to say, the joint distribution on the variables can be factored as $\prod_{i=1}^n \mathcal{P}[V_i=v_i\mid \pa(V_i)]$ where $\pa(V_i)$ is the assignment to the parents of $V_i$.
A $k$-MixBND on $G$ is a convex combination of $k$ BNDs. Such a distribution arises naturally if there is an additional, ``confounding'' random variable $U$ with range $\{1,\ldots,k\}$. Then the distribution on the variables in $G$ (the ``observables'') is the BND of a larger graph in which $U$ is an extra vertex having directed edges to every vertex of $G$. \emph{The main complexity parameter of the problem is $k$.}

Models of this type (and even more general) are fundamental to research in Causal Inference, where the chief difficulty has to do with performing statistical inference for collections of random variables which, besides some known potential direct effects on each other (the directed edges), are also affected by unseen confounders.

One extremely special case has been of longstanding interest in the theory literature: where $G$ is empty. Such a distribution is a mixture of $k$ product distributions or $k$-MixProd. See Fig.~\ref{small-examples}.
\begin{figure}[h]
    \[ \entrymodifiers={+[o][F-]} \xymatrix{
        *{} & \dottedsmallc{U} \ar@{..>}[ld] \ar@{..>}[d] \ar@{..>}[rd] \ar@{..>}[rrd] &*{} & *{} & *{}
        &
        *{} & \dottedsmallc{U} \ar@{..>}[ld] \ar@{..>}[d] \ar@{..>}[rd] \ar@{..>}[rrd] &*{} & *{}
        \\
        V_1 \ar[r] \ar@/_1.5pc/[rr] & V_2 \ar[r] & V_3 \ar[r] & V_4 & *{}
        &
        V_1  & V_2  & V_3 & V_4
        } \]
    \caption{(a) A small Bayesian Network, with latent variable $U$. (b) In the empty graph, a $k$-MixBND is a $k$-MixProd (mixture of product distributions).}  \label{small-examples}  \end{figure}

In this paper we study the \emph{Identification} problem for $k$-MixBNDs. Specifically, given the graph $G$, and given a joint distribution $\mathcal{P}$ on the variables (vertices), recover up to small statistical error (a) the mixture weights, up to a permutation of
the constituents, and (b) for every mixture constituent and for every vertex $V$, its conditional distribution
given each possible setting to its parents. We assume that each individual random variable (vertex) ranges in a finite set of size $\leq d$, which we show w.l.o.g.\ (and at mild cost in complexity), can be taken as $2$.
Clearly, $k$-MixBND models are not always identifiable.
For instance, if $k>1$ and $G$ consists of a single vertex; more generally if two constituents of the mixture produce the same BND, or if one BND is the convex combination of two others. Thus,
another contribution of our paper is to establish mild conditions on the size of $G$ and on the constituent
BND's, that suffice to guarantee identifiability.

The following theorem gives our main results in somewhat simplified form (see Theorems~\ref{thm: product distributions},\ref{thm: main k-BND} for more detail). Here every observable is binary; $\zeta$-separated and $\zeta$-informative are conditions (see later) enabling \emph{quantitative} guarantees on identifiability; the mixture weights are $\pi_i$; and ${\varepsilon}$ represents the distance of the output model from the true model.

\begin{theorem}\label{thm: main}
The algorithms given later w.h.p.\ identify a $k$-MixBND, with the following runtime (including sample) complexity.
\begin{enumerate} 
\item Suppose $\Delta=0$ (i.e., $k$-MixProd) and the $n\geq 3k-3$ observables are 
$\zeta$-separated. \\ The runtime is $(1/\zeta)^{O(k \log k)} (n \log n) (1/\min \pi_i)^{O(\log k)} (1/\varepsilon)^{2}$.
\item Suppose $\Delta\geq 1$, $n \geq \Omega(k)$, and the $k$-MixBND is $\zeta$-informative. \\ The runtime is 
$ (1/\zeta)^{O(k(\Delta^2+ \log k))}(n \log n) (1/\min \pi_i)^{O(\log k)} (1/\varepsilon)^{2}.$
\end{enumerate}
\end{theorem}

In the $k$-MixProd result it is actually sufficient that there \emph{exist} a subset of $3k-3$ \, $\zeta$-separated random variables. Under this weaker assumption, the runtime above is multiplied by $n^{O(k)}$. (The sample size is unaffected.)

\subsection{Background}
To the best of our knowledge, the only prior work on $k$-MixBND estimation problem in the theory literature deals with the
special case of an empty graph $G$ (the joint distribution conditional on the confounder $U$ is
a product distribution), see Fig.~\ref{small-examples}(b). This is
better known as the problem of identifying a mixture of product distributions, $k$-MixProd. It has been studied for nearly 30
years~\cite{KMRRSS94,CGG01,FM99,FOS08,CR08,tahmasebi2018identifiability,ChenMoitra19,gordon2021source} (and we do not even try to list the extensive analogous literature for parametrized distributions over $\mathbb R$). There are actually two versions of the problem: (1)
\emph{Learning} the model, namely, producing any model consistent with (or close to) the observations; (2) \emph{Identifying} the model, namely, producing the true model (or one close to it).
These are incomparable problems, as the second goal is stronger but not always possible. Our paper is concerned entirely with Identification, because the underlying motivation is to be able to predict the effect of interventions in a system (represented by the graph). Such predictions will not be valid when the model is not identifiable.

In~\cite{FOS08} a seminal algorithm for $k$-MixProd was given. Its running
time, for mixtures on $n$ binary variables ($n$ suff.\ large), is $n^{O(k^3)}$.
This was improved in~\cite{ChenMoitra19} to $k^{O(k^3)} n^{O(k^2)}$. The most recent algorithm~\cite{gordon2021source}
 identifies a mixture of $k$ product distributions on at least $3k-3$ variables in time
$2^{O(k^2)} n^{O(k)}$, under a mild ``separation'' condition that excludes unidentifiable instances (see 
below). Under somewhat stricter separation, the time complexity improves to $2^{O(k^2)} n$. In this
paper we improve the latter results, giving upper bounds with exponent $O(k\log k)$ rather than $O(k^2)$. This is nearly optimal, as it is
known that even in the special case that the observables are known to be iid conditional on $U$---call this the 
$k$-MixIID problem---exponent $\Omega(k)$ is unavoidable~\cite{RSS14}. Moreover $n\geq 2k-1$ observable variables are generally necessary.
(There are interesting special cases in which much better efficiency is possible.
\cite{ChenMoitra19} study ``mixtures of sub-cubes,'' $k$-MixProd with all bit probabilities
in $\{0,\frac 1 2, 1\}$, and achieve complexity $n^{O(\log k)}$.)

It turns out that the $k$-MixIID problem plays a special role at the base of a tower of reductions. On the one hand it can be solved using a relatively simple, two-century-old method of Prony~\cite{Prony1795}, which connects it to the classical Hausdorff moment problem; it can also be solved by the Matrix Pencil Method~\cite{RSS14,LRSS15,KKMMR18,gordon2020sparse}. Both analyses are quite recent however.
In turn, the $k$-MixProd problem was solved by a nontrivial reduction to $k$-MixIID~\cite{gordon2021source}. And, in the present paper, we show that, by a rather complex reduction involving subjecting the model to many conditionings, the $k$-MixBND problem can in turn be reduced to $k$-MixProd.

Thus, one contribution of our paper is to add to the corpus of evidence of the amazing power of the $k$-MixIID problem. (Apart from its many direct applications---op.\ cit.) It turns out that if one wants to identify mixtures of Bayesian network distributions, then
solving this extremely restricted case is (possibly up to a logarithmic factor
in the exponent) the complexity bottleneck.

As already noted, a $k$-MixBND can be unidentifiable (i.e., there may be two distinct models mapping to the same distribution on the observables). There are two fundamental ways this can occur. One is that the particular distribution is degenerate. We will rule out such degeneracy by assumptions we call $\zeta$-separated, and later, $\zeta$-informative. The second is that the graph $G$ is too ``impoverished''---the simplest way for this to occur is if $n$ is too small. (As noted, even in the empty graph, $n\geq 2k-1$ is needed information-theoretically.) But it can also happen if $G$ has too many edges. A sufficient condition to rule out such obstacles is that $n\geq \Omega(k (\Delta+1)^4)$.

The parameter $k$ is vital. In a Bayesian network where there is an
unobserved $U$ which can directly affect every observed variable, causal identification is not possible in
general. Each value $u\in U$ could determine completely the values of
all observable variables. There is no way to distinguish this model from models that include direct causal
relationships among the observables. However, the additional assumption of bounded range for $U$
changes the situation, and in this paper we show that there indeed are many graphs for which causal
identification is possible.

\subsection{Bayesian networks and causal identification}
The primary application of Bayesian networks (a term coined in~\cite{pearl1985bayesian}) is to model complex
causal relations. Concrete applications abound: mining data on diverse populations or multi-purpose
usage, bioinformatics, etc. In a wide range of complex natural and artificial systems causal relations are probabilistic,
and moreover inquiry using controlled experiments is impossible or prohibitively expensive. In such settings
researchers often have to resort to collecting samples of an assortment of measurements that are related in
a pattern of probabilistic causal relations that govern the joint distribution. The conditional dependencies among
the components of the joint distribution are traditionally modeled as a Bayesian network. Often, practitioners
have domain-specific knowledge that leads to the construction of a Bayesian network representing the
direct causal relationships among the variables under examination. However, if there are multiple populations
for which the strength of the causal relationships differ, the observed statistics are a convex combination
of the statistics of the different populations. Multiplicity of populations, or equivalently an unobserved discrete
confounding variable $U$, will create joint distributions on the observed variables that obscure the true causal
behavior of the system. On the other hand, while such deviations of the observables from the assumed causal
structure pose a challenge, they also contain valuable information about the presence and behavior of the hidden
variable $U$ that created them. For some critical algorithms in this literature see~\cite{shpitserP06,huangV09}, and for a broader survey~\cite{Pea09,SpGlS00,PetersJS17,KF09}.

We briefly outline (for binary observables) how we solve $k$-MixBND using $k$-MixProd as a subroutine.
The algorithm requires a
set of $3k-3$ nodes of $G$ with disjoint \emph{Markov Boundaries.} See Figure~\ref{fig:disjointmb}
for an illustration. In particular, if the maximum degree in the undirected skeleton of $G$ is $\Delta\ge 2$ and the
number of observables is $n \geq (\Delta+1)^4 (3k-3)$, then this condition is satisfied. (If the skeleton of $G$ happens to be a path,
then we only need a milder condition that $n\ge 2(3k-3)$.)
\begin{figure}[h]
    \centering
    \scalebox{.60}{
        \begin {tikzpicture}[-latex ,auto ,node distance =2 cm and 2 cm ,on grid , semithick, state/.style ={ circle, draw, minimum width =1 cm}]
        \filldraw[color=blue, fill=blue!4, ultra thick](-.8,2.8) rectangle (4.8, -.8);
        \filldraw[color=red, fill=red!4, ultra thick](5.2,2.8) rectangle (14.8, 1.2);
        \filldraw[ultra thick, color = green!60!black, fill=green!4](11.2,-.8) -- (11.2, .8) -- (17.2, .8) -- (17.2, 2.8) -- (18.8, 2.8) -- (18.8, .8) -- (20.8, .8)-- (20.8, - .8) -- cycle;
        \filldraw[color=orange, fill=orange!4, ultra thick](23.2,2.8) rectangle (24.8, -.8);
        \node[state, fill=blue!50] (X1) at (0,0) {$V_1$};
        \node[state] (V1) at (2,2) {$V_2$};
        \node[state] (V2) at (4,0) {$V_3$};
        \node[state] (V3) at (6,2) {$V_4$};
        \node[state] (V4) at (8,0) {$V_5$};
        \node[state, fill=red!50] (X2) at (10,2) {$V_6$};
        \node[state] (V5) at (12,0) {$V_7$};
        \node[state] (V6) at (14,2) {$V_8$};
        \node[state, fill=green!50] (X3) at (16,0) {$V_9$};
        \node[state] (V7) at (18,2) {$V_{10}$};
        \node[state] (V8) at (20,0) {$V_{11}$};
        \node[state] (V9) at (24,0) {$V_{12}$};
        \node[state, fill=orange!50] (X4) at (24,2) {$V_{13}$};
        \node[state, dashed] (U) at (12,  6) {$U$};

        \path[ultra thick] (X1) edge (V2) (V1) edge (V2) (V1) edge (V3) (V2) edge (V4) (V3) edge (X2) (X2) edge (V6) 
(V4) edge (V5) (V5) edge (X3) (X3) edge (V7) (X3) edge (V8) (V8) edge (V9) (V9) edge (X4);
        \path[ultra thick, dashed] (U) edge[bend right = 35] (X1) (U) edge[bend right = 25] (V1) (U) edge[bend right = 35] (V2) (U) edge[bend right = 35] (V2) (U) edge[bend right = 20] (V3) (U) edge[bend right = 25] (V4) (U) edge[bend right = 15] (X2) (U) edge[] (V5) (U) edge[bend left = 20] (V6) (U) edge[bend left = 25] (X3) (U) edge[bend left = 20] (V7) (U) edge[bend left = 35] (V8) (U) edge[bend left = 35] (V9) (U) edge[bend left = 35] (X4);
        \end{tikzpicture}}
    \caption
        A Bayesian network. Four vertices with disjoint Markov boundaries are indicated.} \label{fig:disjointmb}
\end{figure}
We then execute a sequence
of ``runs" of the $k$-MixProd algorithm on the selected set of variables, for all the possible restrictions (assignments of
values) to their Markov boundaries. The main challenge is then to align the outcomes of all the runs, and then to
propagate the conditional probabilities through the entire network, in order to recover the model, a procedure we call
``unzipping.'' Thus, the reduction to $k$-MixProd relies on three contributions of this paper: the construction of a good 
collection of runs, the alignment algorithm, and the Bayesian unzipping algorithm.

Each individual run of $k$-MixProd employs (a variation on) the method of \emph{bootstrapping synthetic bits,} introduced in~\cite{gordon2021source}. Our contribution here
which leads to the near-optimal bound of $\exp(k \log k)$ 
 (rather than $\exp(O(k^2))$)
is a method to bound the singular values of the matrix of multilinear moments of the
distribution on observables.

\subsection{Related work}
Some loosely related work includes learning hidden Markov models~\cite{hsuKZ12,anandkumar2012method,SharanKLV17},
an incomparable line of work to our question, but with somewhat similar motivation. In the same vein, some papers study learning
mixtures of Markov chains from observations of random paths through the state space~\cite{BGK04,GKV16}. These models,
too, differ substantially from the models addressed in this paper, and pose very different challenges. 
Literature on causal structure learning~\cite{spirtes2000constructing, causationreview}
answers the question of identifying the \emph{presence} of hidden confounders.  Fast Causal Inference (FCI) harnesses
observed conditional independence to learn causal structure, which can detect the presence of unobserved variables when
the known variables are insufficient to explain the observed behavior. This literature includes the MDAG problem in which the
DAG structure may depend upon the hidden variable, see~\cite{thiesson1998} for heuristic approaches to this problem. Other
related works study causal inference in the presence of visible ``proxy'' variables which are influenced by a latent
confounder~\cite{miao2018identifying, kuroki2014measurement}.  This has more recently given rise to multiple causal inference,
which assumes multiple causes which are independent when conditioned on a
confounder~\cite{heckerman2018accounting,ranganath2018multiple,wang2018blessings}. These settings are given with
respect to a known causal inference task and significantly less general assumptions on the behavior of the proxy variables.
Finite mixture models have been the focus of intense research for well over a century, since pioneering work in the late
1800s~\cite{Newcomb86,Pearson94}, and doing justice to the vast literature that emanated from this endeavor is impossible
within the scope of this paper. See, e.g., the surveys~\cite{Everitt1981,TSM85,lindsay1995mixture,McLachlanLR19}.

\section{Mixtures of Product Distributions} \label{sec:empty}
In this section, we describe and analyze an algorithm for identifying a mixture of product distributions
(a.k.a.\ $k$-MixBNDs on the empty graph). It improves substantially the previously
best sample size and time complexity in~\cite{gordon2021source}, and it nearly matches the known
lower bound that holds even in the simpler case that the observables are iid~\cite{RSS14}. The chief novelty of this section is not the algorithm, which is similar to that in~\cite{gordon2021source} although one significant change has been necessary, but the analysis, which is entirely different and stronger.

For this reason we arrange the section in a perhaps unusual order: after some necessary definitions, first we give the crucial new theorem which is essential to the sample complexity bound. Only afterward we give the algorithm and the theorem which analyzes its performance.

\subsection{Problem statement. Hadamard extensions and their condition number}

\paragraph{Problem statement and notation.}
We have a hidden variable $U$ ranging in $[k]$, with distribution $\pi_j:=\Pr(U=j)$.
There are $n$ binary observable variables $V_1,V_2,\dots,V_n$; their distributions are given by the rows ${\mathbf{m}}_i$, with entries
${\mathbf{m}}_{ij}:= \Pr[V_i=1\mid U=j]$, of an $n\times k$ matrix ${\mathbf{m}}$.
Thus, the identification problem is to recover $(\pi_j)_{j \in [k]}$
and $({\mathbf{m}}_{ij})_{i\in [n],j\in[k]}$ (up to permuting $[k]$) to within high accuracy.

Identification is not always possible, as there can be disparate models with the same observable distributions. To set the conditions for identifiability, we need the following definition.
\begin{definition}
    An observable $V_i$ is \emph{$\zeta$-separated} if $\min_{j\ne j'} |{\mathbf{m}}_{ij} - {\mathbf{m}}_{ij'}| >\zeta$.
\end{definition}
It is known that a mixture of product distributions (with $\pi_{\text{min}}>0$) is identifiable (from perfect statistics) if it has at least $2k-1$ \, $0$-separated observables~\cite{tahmasebi2018identifiability}; in general, $2k-1$ are  necessary~\cite{RSS14}. Clearly, getting a sample size bound requires also a quantitative assumption, hence the $\zeta$-separation assumption. Our running time and sample complexity guarantee require a slightly larger number of observables: $3k-3$. It should be noted that the post-sampling runtime, although exponential in $k$, is hardly affected by $\zeta$. The bottleneck resource is sample size.

(A similar algorithm works with just $2k-1$ \, $\zeta$-separated observables,
and its post-sampling runtime is about the same, but its sample size requirement is larger.)

As in earlier works~\cite{ChenMoitra19,gordon2021source}, the
algorithm and its analysis make extensive use of the Hadamard extension of a matrix ${\mathbf{m}}$.
\begin{definition}
    The {\em Hadamard product} for row vectors $u=\Paren{u_1,\dotsc,u_k}$, $v=\Paren{v_1,\dotsc,v_k}$ is
    the mapping $\odot: \mathbb R^k \times \mathbb R^k  \to \mathbb R^k$ given by $u \odot v \coloneqq (u_1v_1,\ldots,u_kv_k)$.
    Equivalently, using the  linear operator $v_\odot = \diag(v)$,
    the Hadamard product is $u \odot v = u \cdot v_\odot$. The identity element for the Hadamard product is the all-ones vector ${\mathbb{1}}$. \end{definition}

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

\begin{definition} $\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$.) \end{definition}

In order to prevent some confusions, in the rest of this section, for a matrix $A$ we
denote by $A_{i*}$ its $i$-th row $i$ and by $A_{*j}$ its $j$-th column. Thus e.g., row ${\mathbf{m}}_i$ becomes ${\mathbf{m}}_{i*}$.

\begin{theorem}\label{thm: singular}
    Let ${\mathbf{m}}$ have $k-1$ rows, every row $\zeta$-separated, and let ${\mathbf{M}} = \mathbb H({\mathbf{m}})$. List the singular values of ${\mathbf{M}}\pi_\odot$ in decreasing order: $\sigma_1\geq \sigma_2 \geq \ldots$, and note $\sigma_{k+1}=\ldots=\sigma_{2^{k-1}}=0$.
    Then
    $$ \sigma_k \geq 
        \frac{(\min \pi_i) (\zeta/2\sqrt{5})^{k-1}}{\sqrt{k(k^2+1)^{k-1}}}.
    $$
\end{theorem}
\begin{proof} The proof relies on the following insight. Since ${\mathbf{M}}$ has dimensions $2^{k-1} \times k$, $\sigma_k$ characterizes the least norm of ${\mathbf{M}} v$ ranging over any unit vector $v$, but it does \emph{not} characterize the least norm of vectors of the form $h {\mathbf{M}}$; the left-kernel of ${\mathbf{M}}$ is of course very large. The insight is that it \emph{does} become possible to bound $\sigma_k$ in terms of such vectors $h$, \emph{provided $h$ is restricted to rank $1$ tensors.}

    Let $\tau({\mathbf{m}},\pi)=\min_{0\neq h\in \mathbb H_1} \Norm{h^\mathsf{T} {\mathbf{M}} \pi_\odot}/\Norm{h}$.
    The proof relies on two lemmas. 
    \begin{lemma}
        $\sigma_k({\mathbf{M}}\pi_{\odot}) \ge \tau({\mathbf{m}}, \pi) / \sqrt{k \prod_1^{k-1}(1+\|{\mathbf{m}}_{i*}\|^2)}$. \end{lemma}
    \begin{proof}
        Consider $v\in \mathbb R^k$, $\|v\|=1$, achieving $\sigma_k$, i.e.,
        $r\coloneqq \mathbb H({\mathbf{m}}) \pi_\odot v$ satisfies $\Norm{r}=\sigma_k$. W.l.o.g.\ $\|v_k\|\geq 1/\sqrt{k}$. Now
        $\mathbb H({\mathbf{m}})_{*k}=\frac{1}{\pi_k v_k} \Paren{r-\sum_{j=1}^{k-1} \pi_j v_j \mathbb H({\mathbf{m}})_{*j}}$.

        Now we carefully choose $h \in \mathbb H_1$ based on ${\mathbf{m}}$. Define ${\mathbf{n}} \in \mathbb R^{[k-1]}$ by ${\mathbf{n}}_i\coloneqq-1/{\mathbf{m}}_{ii}$; and let $h \coloneqq \Paren{\prod_1^{k-1} {\mathbf{m}}_{ii}} \mathbb H({\mathbf{n}})$. For $j \neq k$ we have
        \begin{eqnarray*}
            h^\mathsf{T} \mathbb H({\mathbf{m}})_{*j}
            & = & \Paren{\prod_1^{k-1} {\mathbf{m}}_{ii}} \sum_S {\mathbf{n}}_S {\mathbf{m}}_{S,j}
            =  \Paren{\prod_1^{k-1} {\mathbf{m}}_{ii}} \sum_S (-1)^{|S|} \Paren{\prod_{i \in S} 1/{\mathbf{m}}_{ii}}\Paren{\prod_{i \in S} {\mathbf{m}}_{ij}} \\
            & = & \sum_S (-1)^{|S|} \Paren{\prod_{i \notin S} {\mathbf{m}}_{ii}} \Paren{\prod_{i \in S} {\mathbf{m}}_{ij}}
            =  \prod_{i=1}^{k-1} ({\mathbf{m}}_{ii}-{\mathbf{m}}_{ij}) = 0.
        \end{eqnarray*}
        So, $h^\mathsf{T} \mathbb H({\mathbf{m}})_{*j} \pi_j = 0$ for $j = 1,\dotsc,k-1$.
        For $j=k$ we can also show that $ h^\mathsf{T} \mathbb H({\mathbf{m}})_{*j} \pi_j$ is small:
        \[
            \Paren{h^\mathsf{T} \mathbb H({\mathbf{m}})\pi_\odot}_k
             =  \frac{1}{v_k} h^\mathsf{T} \Paren{r-\sum_{j=1}^{k-1} \pi_j v_j \mathbb H({\mathbf{m}})_{*j}}
            =  \frac{1}{v_k} h^\mathsf{T} r - \sum_{j=1}^{k-1}\pi_j v_j \Paren{h^\mathsf{T} \mathbb H({\mathbf{m}})_{*j}} 
             =  \frac{1}{v_k} h^\mathsf{T} r.
        \]
        The norm of $h^\mathsf{T} \mathbb H({\mathbf{m}})\pi_{\odot}$ is then upper bounded by
        \begin{eqnarray*}
            \Norm{h^\mathsf{T} \mathbb H({\mathbf{m}})\pi_\odot}
            & = & \Abs{(h^\mathsf{T} \mathbb H({\mathbf{m}})\pi_{\odot})_k} = \frac{1}{v_k}h^\mathsf{T} r
            \leq  \sqrt{k} \Norm{h} \Norm{r}                                                          \\
            & = & \sqrt{k}\Norm{h}\sigma_k = \sigma_k \sqrt{k\prod_{1}^{k-1}\Paren{1+{\mathbf{m}}_{ii}^2}} \quad \text{due to $h$ being a rank $1$ tensor} \qedhere
        \end{eqnarray*}
    \end{proof}

    Next we lower bound $\tau({\mathbf{m}},\pi)$.
    \begin{lemma}  $\tau({\mathbf{m}},\pi) \geq (\min \pi_i) (\zeta/2\sqrt{5})^{k-1}$. \end{lemma} \begin{proof} Consider
        any ${\mathbf{G}} \in \mathbb H_1$, say ${\mathbf{G}} = \mathbb H(\gg)$, $\gg\in\mathbb R^{[k-1]}$. Then
        $({\mathbf{G}}^\mathsf{T} \mathbb H({\mathbf{m}}) \pi_\odot)_j = \pi_j \sum_S \gg_S {\mathbf{m}}_{S,j} = \pi_j \prod_1^{k-1} (1+\gg_i {\mathbf{m}}_{i,j})$.
        We also note that $\|{\mathbf{G}}\|=\sqrt{\prod_1^{k-1} (1+\gg_i^2)}$.

        We now show that there is some $j$ such that $\prod_{i=1}^{k-1} \Abs{ \frac{1+\gg_i {\mathbf{m}}_{ij}}{\sqrt{1+\gg_i^2}} }$ is large.
        First, for any $i$ for which $\gg_i\geq \frac 1 2$, there is at most one $j$ s.t.\ ${\mathbf{m}}_{ij}< \zeta$; exclude these $j$'s. Next,
        for each $i$ for which $\gg_i<\frac 1 2$, there is at most one $j$ s.t.\ $\Abs{\frac{1}{\gg_i}+{\mathbf{m}}_{ij}}<\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
        $\gg \in \mathbb R^{k-1}$, such a $j$ exists. We consider three cases for $i$. In each case we lower bound
        $\Abs{ \frac{1+\gg_i {\mathbf{m}}_{ij}}{\sqrt{1+\gg_i^2}}} $.
        \begin{enumerate}
            \item $\gg_i\geq 1/2, {\mathbf{m}}_{ij} \geq \zeta$. Then $\Abs{\frac{1+\gg_i {\mathbf{m}}_{ij}}{\sqrt{1+\gg_i^2}}} \geq {\mathbf{m}}_{ij} \geq \zeta$.
            \item $-1/2 \leq \gg_i< 1/2$. Then $\Abs{\frac{1+\gg_i {\mathbf{m}}_{ij}}{\sqrt{1+\gg_i^2}}} \geq \sqrt{\frac{({\mathbf{m}}_{ij}-2)^2}{5}} \geq 1/\sqrt{5}$.
                  (This does not depend on $j$.)
            \item $\gg_i< -1/2, |\frac1{\gg_i}+{\mathbf{m}}_{ij}|\geq \zeta/2$. Then $\left| \frac{1+\gg_i {\mathbf{m}}_{ij}}{\sqrt{1+\gg_i^2}}\right| =\left| \frac{\gg_i (\frac1{\gg_i}+{\mathbf{m}}_{ij})}{\gg_i \sqrt{1+1/\gg_i^2}}\right|
                      \geq \frac{\zeta}{2\sqrt{5}}$.
        \end{enumerate}
        We therefore have $\tau({\mathbf{m}},\pi) \geq (\min \pi_i) (\zeta/2\sqrt{5})^{k-1}$.\end{proof}
    Finally for the theorem combine the above lemmas and the fact that $\|{\mathbf{m}}_{i*}\|\leq k$ for any $i$.
\end{proof}

\subsection{Algorithm for mixtures of product distributions: overview and theorem}
We use the synthetic bits method introduced
in~\cite{gordon2021source}.
Let $V_S $ be the random variable $ \prod_{i\in S} V_i$.
Our entire access to information about the source is that for every
$S\subseteq [n]$ we have an empirical estimate $\emom(S)$ of $\operatorname{\mathbb{E}}[V_S]$, defined as the fraction of samples
in which $\prod_{i\in S} V_i = 1$. Notice that $\operatorname{\mathbb{E}}[V_S] = {\mathbf{M}}_S \pi_\odot {\mathbb{1}}$
(where as previously, ${\mathbf{M}}=\mathbb H({\mathbf{m}})$). This is a multilinear
moment of the ${\mathbf{m}}_{ij}$.
For the algorithm outlined here to produce ${\mathbf{m}}$ and $\pi$ within additive error $\varepsilon$ (in each
coordinate), we shall need the $\emom(S)$ that we use to be accurate to within
$\varepsilon (\min \pi_i)^{O(\log k)} \zeta^{O(k \log k)}$. This can be achieved with high probability using a
sample of size $(1/\varepsilon)^2 (1/\min \pi_i)^{O(\log k)} (1/\zeta)^{O(k \log k)}$. To avoid introducing more
notation here, the algorithm is described with respect to the exact statistics. The complete analysis in
Appendix~\ref{app: product} includes the analysis of the effects of sampling error.

${\mathbf{M}}$ has $k$ columns.
Thus, if $n$ is not too small, it will have linear dependencies
among its rows.
We will create multiple ``synthetic'' copies of a single row of ${\mathbf{M}}$ by taking linear
combinations of other rows. If each copy uses a distinct set of rows of ${\mathbf{m}}$ in the combinations that are added up, the corresponding
synthetic variables are statistically independent (conditional on $U$). This enables computing higher moments of a random variable (in particular of the $V_i$)
and ultimately reduces the mixtures of products problem to identification of mixtures of products of iid variables. This is an easier problem which can be solved
within the required guarantees using either the matrix pencil method~\cite{KKMMR18} or Prony's
method~\cite{gordon2020sparse} (the differences between the two methods are inconsequential for the application
here). Finally, after recovering a good estimate for one row and $\pi$, it is a simple task to solve for the rest of the
model.

The implementation of the above outline as follows. We need three disjoint sets $S,T,T'\subset [n]$, each of size
$k-1$, such that each of these $3k-3$ variables are $\zeta$-separated.
These can be found by enumeration over at
most $n^{O(k)}$ choices, and the rest of the algorithm needs to be repeated for each choice, until a solution is
found.

\begin{definition} For a matrix $A$ and a set
    of row indices $S$, denote by $A[S]$ the sub-matrix derived from $A$ by taking the rows in $S$.
    The {\em sum of two collections} of subsets $\mathcal{X},\mathcal{Y}\subseteq 2^{[n]}$ is defined as
    $\mathcal{X} + \mathcal{Y} \coloneqq \Set{X\cup Y: X\in \mathcal{X}, Y\in \mathcal{Y}}$.
    Let $S,T\subset [n]$ be two disjoint sets. Let $x$ and $y$, respectively, be vectors indexed by the
    subsets in $\mathcal{X} = 2^S$ and $\mathcal{Y} = 2^T$, respectively. The {\em Kronecker product}
    $x \otimes y$ is the vector indexed by the subsets in $\mathcal{X} + \mathcal{Y}$ given by
    $(x \otimes y)_{X\cup Y} \coloneqq x_{X}y_{Y}$, for all $X\in \mathcal{X}$ and $Y\in \mathcal{Y}$.  Note that $\mathcal{X} + \mathcal{Y} = 2^{S \cup T}$.
    (Notice that every set $Z\in \mathcal{X} + \mathcal{Y}$ can be written in this case
    uniquely as $Z = X\cup Y$ for $X\in \mathcal{X}$ and $Y\in \mathcal{Y}$.) 
\end{definition}

Let $\mathcal{A} = 2^S$, $\mathcal{B} = 2^T$, and $\mathcal{B}' = 2^{T'}$. Obviously, Theorem~\ref{thm: singular}
applies to all three matrices ${\mathbf{M}}[\mathcal{A}]$, ${\mathbf{M}}[\mathcal{B}]$, and ${\mathbf{M}}[\mathcal{B}']$.
W.l.o.g., $\{1\}\in\mathcal{B}$. We generate the first $2k$ power moments
of $V_1$ as follows.
Define $\CC_{\calB\calA} = {\mathbf{M}}[\mathcal{B}]\pi_\odot {\mathbf{M}}[\mathcal{A}]^{\mathsf{T}}$, and similarly
$\CC_{\calB'\calA} = {\mathbf{M}}[\mathcal{B}']\pi_\odot {\mathbf{M}}[\mathcal{A}]^{\mathsf{T}}$ and
$\CC_{\calB+\calB', \calA} = {\mathbf{M}}[\mathcal{B} + \mathcal{B}']\pi_\odot {\mathbf{M}}[\mathcal{A}]^{\mathsf{T}}$. Notice that
the entries of all of these matrices are multilinear moments, and hence we
have good empirical estimates of these matrices with high probability.
We will compute vectors $v^{(1)},v^{(2)},\dots,v^{(2k)}\in\R^{\mathcal{A}}$
iteratively as follows. Initialize $v^{(1)} = (\mom(S\cup\{1\}))_{S\in\mathcal{A}}$.
Notice that $v^{(1)}_{\{1\}} = \mom(\{1\}) = {\mathbf{m}}_1\pi_\odot{\mathbb{1}}$.
Then, for $i = 2,3,\dots,2k$, writing $i = 2^r+ j$ for $j\in\{1,2,\dots,2^r\}$,
put
$$
    v^{(i)} = \left(\left(v^{(j)} \CC_{\calB\calA}^{-1}\right) \otimes \left(v^{(2^r)} \CC_{\calB'\calA}^{-1}\right)\right) \CC_{\calB+\calB', \calA}.
$$
A simple induction shows that $v^{(i)}_{\{1\}} = {\mathbf{m}}_1^{\odot i}\pi_\odot{\mathbb{1}}$, the
$i$-th power moment of $V_1$. The depth of the recurrence for computing the $v^{(i)}$-s
does not exceed $O(\log k)$, and this, along with Theorem~\ref{thm: singular}, controls
the accumulation of the error from the empirical values of the multilinear moments. Given
the $2k$ moments of $V_1$, we can recover $\pi$ and ${\mathbf{m}}_1$ using, for instance, the
above-mentioned result on Prony's method. The rest of ${\mathbf{m}}$ can be recovered easily
using this information and the empirical estimates of bilinear moments.

A detailed and rigorous description and analysis of the algorithm for mixtures of product
distributions appears in Appendix~\ref{app: product}. In particular, it includes a proof of
the following theorem, summarizing the result.
\begin{theorem}\label{thm: product distributions}
    There is an algorithm for identifying mixtures of $k$ product distributions on $n$ bits with the following
    guarantees. Let $\varepsilon > 0$ and suppose that at least $3k-3$ observable bits are $\zeta$-separated.
    The algorithm runs in time
    $(1/\varepsilon)^2 (1/\min \pi_i)^{O(\log k)} (1/\zeta)^{O(k \log k)} n^{O(k)}$.\footnote{This includes the
        sample size. The number of samples that is required is
        $(\log n/\varepsilon^2) (1/\min \pi_i)^{O(\log k)} (1/\zeta)^{O(k \log k)}$; each sample has $n$ bits.}
    With high probability (over the distribution of the samples), the algorithm computes $\tilde{\pi}$ and
    $\tilde{{\mathbf{m}}}$ such that $\max_i |\tilde{\pi}_i - \pi_i|\le \varepsilon$ and $\max_{ij} |\tilde{{\mathbf{m}}}_i^j - {\mathbf{m}}_i^j|\le \varepsilon$.
\end{theorem}

\emph{Comment:}
In Section~\ref{sec:genl}, the general $k$-MixBND inference problem will rely on applications of this algorithm as a black box.
The improvement given in this section is required in order to achieve the stated sample size and time complexity of solving the
general case.

\section{General Bayesian Networks} \label{sec:genl}
\subsection{Describing Bayesian networks}

A Bayesian network consists of a directed acyclic graph (DAG) $G=(\VV,E)$, where $\VV$ is a set of random variables
whose values are distributed according to a joint probability distribution $\mathcal{P}$. The distribution $\mathcal{P}$ must maintain
the following consistency requirement with respect to $G$. For a vertex $V\in \VV$, let $\Pa(V) = \{U\in\VV:\ (U,V)\in E\}$
be the set of the parents of $V$. (Similarly, let $\Ch(V)$ denote the set of children of $V$.) Let $\VV = \{V_1,V_2,\dots,V_n\}$
be an enumeration of the nodes of $G$. Then, $\mathcal{P}$ factors into
$$
    \mathcal{P}[V_1=v_1,V_2=v_2,\dots,V_n=v_n] = \prod_{i=1}^n \mathcal{P}[V_i=v_i\mid \pa(V_i)],
$$
where $\pa(V_i)$ denotes the assignment $V_j = v_j$ for all  $V_j\in \Pa(V_i)$. In particular, for every $V\in\VV$, conditional
on the values of $\Pa(V)$, $V$ is independent of all other nodes which are non-descendants of $V$ in $G$.

To induce independence among vertices, we will use the concept of the \emph{Markov boundary}~\cite{pearl2014probabilistic}.
\begin{definition}[Markov Boundary] \label{def:mb}
    For a vertex $Y$ in a DAG $G=(\VV,E)$, the \emph{Markov boundary} of $Y$, denoted $\Mb(Y)$, is defined by
    \[ \Mb(Y) \coloneqq \Pa(Y) \cup \Ch(Y) \cup \Pa(\Ch(Y)) \setminus \{Y\}. \]
\end{definition}

\begin{definition}  We'll introduce a parameter $\gamma(G)$ which will appear in the complexity of the identification procedure, which is defined by $\gamma(G) \coloneqq \max_{V\in \VV} \Card{\Mb(V)}$. \label{def:gamma}
\end{definition}

\begin{observation} For any $X,Y \in \VV$, $X\in \Mb(Y) \iff Y \in \Mb(X)$ \label{obs:mb mb}
\end{observation}
Conditioning on the Markov boundary of a vertex allows us to ensure that it is $d$-separated,  and thus independent, from the rest of the graph.
\begin{lemma}[See \cite{pearl2014probabilistic}]\label{lem:mbindependent} For any vertex $V\in \VV$ and subset $S\subseteq V\setminus (\Mb(V)\cup \Set{V})$, $\mathcal{P}(V\mid \Mb(V),S) = \mathcal{P}(V\mid \Mb(V))$.
\end{lemma}

The final definition will be useful when considering vertices in a topological ordering.
\begin{definition}[Top] \label{def:top}
    We will use $\Top(V)$ to denote $\Mb(V) \setminus \Ch(V)$.
\end{definition}

\begin{observation} We can decompose \begin{align*}P(\mb(Y), y) & = \ChProbA{y} P(y\mid \ttop(Y))\mathcal{P}_u(\ttop^a(Y)) \\&= \ChProbA{y} P(y\mid \pa(Y)) \mathcal{P}_u(\ttop^a(Y))
    \end{align*} since $Y$ is $d$-separated from $\ttop(Y)\setminus \pa(Y)$ given $\pa(Y)$. \label{obs:top decomposition}
\end{observation}

\paragraph{Uppercase/lowercase conventions}
Following notation in causal inference literature,  we will use lowercase letters to denote assignments.  For example
$\mathcal{P}(v \mid u) = \mathcal{P}(V=v \mid U=u)$.  Following this convention,  we will write $\pa(V)$, $\ch(V)$, and $\mb(X)$
to denote assignment to the parents $\Pa(V)$, children $\Ch(V)$, and Markov boundary $\Mb(V)$.

\subsection{Applying a $k$-MixProd algorithm}
Our algorithm will induce instances of independence through post-selected conditioning. We will then run a black box $k$-MixProd algorithm on this post-selected subset of the data.  We will now formally define this concept.
\subsubsection{Describing runs}
\begin{definition}[Run]
    A \emph{run} over a graph $G = (\VV,E)$ is a tuple $a=(\mathcal{I}^a, f^a)$ where $\mathcal{I}^a \subseteq \VV$ are variables that we will focus on by conditioning on assignments to the set
    \[ \MCond^a \coloneqq \Paren{\cup_{X \in \mathcal{I}^a:d(X)< d^a_{\max}} \Mb(X)} \bigcup \Paren{\cup_{X\in \mathcal{I}^a:d(X)=d^a_{\max}} \Pa(X)} \] where $d^a_{\max} \coloneqq \max_{X\in \mathcal{I}^a} d(X)$. The value of the assignment is given by $f^a: \MCond^a \to \Set{0,1}$.  We will restrict our attention through the entire paper to \emph{well-formed runs},  i.e.  runs for which
    \[ \mathcal{I}^a \cap \MCond^a = \emptyset. \] We'll call $\mathcal{I}^a$ the \emph{independent set} for $a$, and $\MCond^a$ the \emph{conditioning set}.
\end{definition}

\begin{definition}[Depth of a vertex]
    Given a DAG $G=(\VV, E)$ and any vertex $V\in \VV$, let $d_G(V)$ be the \emph{depth} of $V$ in $G$, i.e. the length of the shortest path from a source vertex to $V$ in $G$. When $G$ is clear from context, we'll omit the subscript.
\end{definition}

\begin{definition}[Bottom vertex in a run] \label{def:bottom vertex} Given a run $a = (\mathcal{I}^a, f^a)$, we'll define the set of bottom vertices of the run to be the subset $\mathcal{J}^a \subseteq \mathcal{I}^a$ of vertices with maximal depth among the vertices in $\mathcal{I}^a$.
\end{definition}

\paragraph{Superscript notation}
We'll write $\mb^a(V)$,$\pa^a(V)$, $\ch^a(V)$ to refer to the assignment to the Markov boundary of $V$, parents of $V$, and children of $V$ as set by run $a$.\footnote{Any quantities parameterized by a run will take the parameter as a superscript.} In a similar spirit, we'll occasionally write $v^0$ to denote the assignment  $V=0$.

\begin{definition}[Distribution induced by a run]
    For any run well-formed run $a$, there is a distribution on the variables in $\mathcal{I}^a$ induced by the run, which we'll denote
    $\mathcal{P}^a(\cdot)$. The distribution is given by
    \[ \mathcal{P}^a(\cdot) = \mathcal{P}(\cdot \mid \mcond^a), \] where $\mcond^a$ is the assignment to $\MCond^a$ in keeping with our conventions.
\end{definition}

The outputs of applying a $k$-MixProd algorithm to  $\mathcal{P}^a(\mathcal{I}^a)$ are a matrix $m^a \in [0,1]^{\Card{\mathcal{I}^a}\times k}$ and a mixture vector $\pi^a \in [0,1]^k$ (satisfying $\sum_u \pi^a_u = 1$) defined by
\begin{align}
    m^a_{iu} & \coloneqq \mathcal{P}^a(X_i=1\mid U_a=u) = \mathcal{P}(X_i=1\mid U_a = u, \MCond^a), & \quad \forall i\in \mathcal{I}^a, u\in [k] \\\pi^a_{u} &\coloneqq \mathcal{P}^a(U_a = u) = \mathcal{P}(U_a=u\mid \MCond^a),&\quad \forall j\in [k],
\end{align} where $U_a$ is a source over $[k]$ with a distribution that can be permuted to equal the distribution of $U\mid \MCond^a$.  Note that because a mixture of products is invariant to a permutation of the labels of the mixture components, a mixture of products algorithm has total freedom to output any permutation of the labels for the source variable it's considering.  In other words, the labels of $U_a$ correspond to a unique unknown label of $U$. We will see how to handle this this in Section~\ref{sec: aligning source labels}.

\subsubsection{Assumptions} \label{sec: assumptions} In order to apply a $k$-MixProd algorithm on the distribution $\mathcal{P}^a(\mathcal{I}^a)$ induced by a run, we need to ensure that the requirements of $k$-MixProd identifiabilty are satisfied.

Moreover, in order for us to be able to stitch together the information gleaned from multiple runs into a coherent mixture distribution over all of $\VV$, we'll need to be able to uniquely determine a mapping from mixture constituents to probability assignments at each variable.  Both of these requirements will be met by distributions satisfying the following criterion:

\begin{definition}[$\zeta$-informative mixture] Let $\mathcal{P}$ be distribution over a set $\Set{0,1}^\VV$ that is a mixture of $k$ distributions that are individually Markov with respect to a DAG $G=(\VV,E)$, i.e., $\mathcal{P}_u(\VV)$ is Markov with respect to $G$ for a variable $U \in [k]$ with some unknown distribution. We'll say that $\mathcal{P}$ is $\zeta$-\emph{informative} if for every variable $V\in \VV$,
    \begin{enumerate}
        \item for every assignment $\mb^a(V)$ to the Markov boundary of $V$ we have that $V$ is $\zeta$-separated in the distribution $\mathcal{P}(V\mid \mb^a(V))$ and that $\mathcal{P}(V\mid \mb^a(V)) \in (\zeta, 1-\zeta)$ for all $u\in [k]$, and
        \item for every assignment $\pa^a(V)$ to the parents of $V$ we have that $V$ is $\zeta$-separated in the distribution $\mathcal{P}(V\mid \pa^a(V))$
    \end{enumerate} \label{def:f-zeta-im}
\end{definition}

\begin{lemma} If $\mathcal{P}$ is $\zeta$-informative, then for any assignment $s$ to a subset of vertices $S\subseteq \VV$, $\mathcal{P}(s) \geq \zeta^{\Card{S}}$. \label{lem:probabilities not too small}
\end{lemma}
\begin{proof} Fix a topological ordering $\Seq{X_1,\dotsc,X_n}$ on the vertices in $\VV$. We'll use induction on the latest variable appearing in $S$, $X_i$. Assume the claim holds for all subsets $R \subseteq \Set{X_1,\dotsc, X_{i-1}}$, and write
    \[ \mathcal{P}(s) = \mathcal{P}(s(X_i)\mid s_{(i-1)})\mathcal{P}(s_{(i-1)}). \]

    Now \begin{align*}\mathcal{P}(s(X_i)\mid s_{(i-1)}) &= \sum_{\pa(X_i)\in \Set{0,1}^{\Pa(X_i)}}\mathcal{P}(s(X_i)\mid \pa(X_i))\mathcal{P}(\pa(X_i)\mid s_{(i-1)}) \\&\geq \zeta \sum_{\pa(X_i)\in \Set{0,1}^{\Pa(X_i)}}\mathcal{P}(\pa(X_i)\mid s_{(i-1)}) = \zeta.\end{align*}  Applying the induction hypothesis, we get that
    \[ \mathcal{P}(s) \geq \zeta \cdot\zeta^{\Card{S} - 1}. \]
\end{proof}

\begin{definition}
    We'll say that a distribution over $\VV$ is \emph{positive} if the probability of any event is $> 0$. \label{def:pos}
\end{definition}

\begin{observation}
    If $\mathcal{P}$ is $\zeta$-informative, then $\mathcal{P}$ is positive.
\end{observation}

We will assume henceforth that our distribution $\mathcal{P}$ is a $\zeta$-informative mixture. This guarantees that the variables in $\mathcal{I}^a$ are separated in $\mathcal{P}^a$ for any run $a$ with $\Card{\mathcal{I}^a} \geq 3k-3$, allowing us to run a $k$-MixProd algorithm.

Our algorithm will begin by executing a collection of runs of the $k$-MixProd algorithm $\mathcal{A}$.  In order to succeed,  we need to ensure that each of these runs is on at least $3k-3$ independent variables, motivating the following definition:
\begin{definition}
    An individual run $a=(\mathcal{I}^a, f^a)$ is \emph{$N$-independent} if $\Card{\mathcal{I}^a} \geq N$.
\end{definition}

\subsection{Using different runs to learn source-specific parameters}
In contrast to the $k$-MixProd setting, in the $k$-MixBND setting we are not always fortunate enough to have a sufficient number of independent variables. Instead, we will have to condition on variables in order to induce an independent set using $d$-separation rules. Specifically, our algorithm will consist of many calls to a black box $k$-MixProd algorithm, where each call looks at the distribution on a subset $\mathcal{I} \subseteq \VV$ conditioned on an assigment to the Markov boundaries of the variables in $\mathcal{I}$ such that we $d$-separate each pair of variables.

\subsection{Aligning source labels across different runs} \label{sec: aligning source labels}
Each run of the $k$-MixProd algorithm will return $\mathcal{P}^a(V\mid U_a = u)$ for some arbitrary permutation $U_a$ of the variable. We need to align all of the outputs so that different outputs share the same permutation of $U$. This can be done if we know ahead of time that a variable in one run's independent set should have the same probability distribution as it does in the other run.
\begin{definition}[Aligned mixture distributions]\label{def:two_mixtures_aligned}
    A pair of $k$-component mixture distributions $\mathcal{P}^{(1)}, \mathcal{P}^{(2)}$ over tuples of random variables $U^{(1)}, (X^{(1)}_1,\dotsc, X^{(1)}_m)$ and $U^{(2)}, (X^{(2)}_1,\dotsc, X^{(2)}_m)$, respectively, is alignable if there exist $i_1, i_2 \in [m]$ such that $\mathcal{P}_u(X^{(1)}_{i_1}) = \mathcal{P}_u(X^{(2)}_{i_2})$ for all $u\in [k]$ and $X^{(1)}$ is separated with respect to $U$. We'll call the random variable $X^{(1)}_{i_1}$ (which is the same as $X^{(2)}_{i_2}$) \emph{the alignment variable}.
\end{definition}

\begin{definition}[Alignment spanning tree]
    We say a set of $\ell$ mixture distributions $\mathcal{P}^{(1)}, \ldots, \mathcal{P}^{(\ell)}$, is alignable if there exists an undirected spanning tree with edges $E \subseteq \Set{\mathcal{P}^{(1)}, \ldots, \mathcal{P}^{(\ell)}} \times \Set{\mathcal{P}^{(1)}, \ldots, \mathcal{P}^{(\ell)}}$ where each edge corresponds to $\AV( \mathcal{P}^{(i)},  \mathcal{P}^{(j)})$,  variables at which $\mathcal{P}^{(i)},  \mathcal{P}^{(j)}$ are aligned.
\end{definition}

The alignment step will take the output from a number of alignable runs and permute the mixture labels,  assigning parameters to a global set of sources.  Pseudocode for this step is given in Algorithm~\ref{alg:align}.
\begin{algorithm}[H]
    \caption{Alignment}\label{alg:align}
    \begin{algorithmic}[1]
        \Require A set of distributions $\mathcal{P}^{(1)}, \ldots,  \mathcal{P}^{(\ell)}$ with $\mathcal{P}^{(i)}$ a distribution with parameters $\mathcal{P}^{(i)}(x^{(i)} \mid u^{(i)})$ for $X^{(i)} \in \mathcal{X}^{(i)}$ and $\mathcal{P}(u^{(i)})$.  In addition, we have a spanning tree with edges $E$ and alignment variables $\AV(\mathcal{P}^{(i)}, \mathcal{P}^{(j)}) \in \mathcal{X}^{(i)} \times \mathcal{X}^{(j)}$.
        \Ensure $\mathcal{P}_u(X_i^{(j)})$ for $X_i^{(j)} \in \mathcal{X}^{(j)}$ for $j \in [\ell]$.  Also $\mathcal{P}(u)$ for each $u \in [k]$.
        \State Let $T$ be an oriented tree on the undirected spanning tree with edges $E$ with all vertices having a directed path to an arbitrary fixed vertex $\mathcal{P}^{(t)}$.
        \For{each edge $\mathcal{P}^{(i)} \to \mathcal{P}^{(j)}$ in $T$ and corresponding alignment variables $\AV(\mathcal{P}^{(i)}, \mathcal{P}^{(j)}) = (X_{\AV}^{(i)}, X_{\AV}^{(j)})$}
        \State Let $\sigma^{i \to j}$ be the permutation on the sources that minimizes $\Norm{\mathcal{P}^{(i)}(X^{(i)}_{\AV} \mid \sigma^{i \to j}(u^{(i)})) - \mathcal{P}^{(j)}(X^{(j)}_{\AV} \mid u^{(j)})}_{\infty}$.
        \EndFor
        \For{each vertex $\mathcal{P}^{(s)}$}
        \State Let $p = (s,  p_1, \ldots,  p_{r} ,t)$ be the directed path in $T$ from $\mathcal{P}^{(s)}$ to $\mathcal{P}^{(t)}$.
        \State Let $\sigma^{s\to t} = \sigma^{s \to p_1} \cdot \sigma^{p_1 \to p_2} \dotsb \cdot \sigma^{p_r \to t}$ be the composition of permutations along the path from $\mathcal{P}^{(s)}$ to $\mathcal{P}^{(t)}$.
        \State Set $\mathcal{P}_u(X_i^{(s)})) \gets \mathcal{P}^{(s)}(X_i^{(s)} \mid \sigma^{s \to t}(u^{(s)}))$ for all $i \in [n]$.
        \State Set $\mathcal{P}(u) \gets \mathcal{P}^{(s)}(\sigma^{s \to t}(u^{(s)}))$
        \EndFor
    \end{algorithmic}
\end{algorithm}

\subsection{Bayesian unzipping: recovering parameters per source}
Recall that our algorithm uses Markov boundary conditioning to induce independent variables.  Hence,  after aligning the sources in runs of the $k$-MixProd algorithm we will have access to $\mathcal{P}_u(Y \mid \Mb(Y))$ for each $Y \in V$.  Our goal is to obtain $\mathcal{P}_u(Y\mid \Pa(Y))$.  Note that
\begin{equation}\label{eq:mbdecomp}
    \mathcal{P}_u(y^1 \mid \mb(Y)) = \frac{\mathcal{P}_u(y^1 , \mb(Y))}{\mathcal{P}_u(y^1 , \mb(Y)) + \mathcal{P}_u(y^0 , \mb(Y))}
\end{equation}
The terms in this fraction are all of the same form and can be factored according to the DAG into
\begin{equation*}
    \mathcal{P}_u(y, \mb^a(Y)) = \mathcal{P}_u(\ttop(Y)) \mathcal{P}_u(y \mid \pa^a(Y)) \underbrace{\prod_{V \in \Ch(Y)}\mathcal{P}_u(v^a \mid f^{a}(\Pa(V) \setminus \{Y\}),  y)}_{\mathcal{P}_u(\ch^a(Y) \mid \ttop^a(Y), y)}.
\end{equation*} See Figure~\ref{fig:MBdecomp} for a concrete example of this decomposition.
After substitution into \eqref{eq:mbdecomp} we see that $\mathcal{P}_u(\ttop(Y))$ appears in both the numerator and denominator because it is independent of the assignment to $Y$. We can cancel these terms out leaving only the following terms:
\begin{enumerate}
    \item $\mathcal{P}_u(y^0 \mid \pa^a(Y))$ and $\mathcal{P}_u(y^1 \mid \pa^a(Y))$, which must sum to 1.
    \item  $\mathcal{P}_u(\ch^a(Y) \mid \ttop^a(Y), y^0)$ and  $\mathcal{P}_u(\ch^a(Y) \mid \ttop^a(Y), y^1)$ which are both the product of the desired parameters of variables later in the topological ordering.  We can ensure we have access to these terms by solving for the parameters of $V \in \VV$ in a reverse-topological ordering.
\end{enumerate}
We can substitute $1-\mathcal{P}_u(y^1\mid \pa^a(Y))$ for $\mathcal{P}_u(y^0\mid \pa^a(Y))$ in the expanded version of \eqref{eq:mbdecomp} to obtain a single equation with only $\mathcal{P}_u(y^1\mid \pa^a(Y))$ as an unknown, which we can then solve.  The pseudocode for this process is given in Algorithm~\ref{alg:bayesian-unzipping}.
\begin{figure}[h]
    \setlength{\belowcaptionskip}{-10pt}
    \centering
    \scalebox{.4}{
        \begin {tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}, -latex ,auto ,node distance =2 cm and 2 cm ,on grid , semithick, state/.style ={ circle, draw, minimum width =1 cm}]
        \filldraw[color = red, fill=red!4](-2.8, 2.8) -- (4.8, 2.8) -- (4.8, -.8) -- (3.2, -.8) -- (3.2, 1.2) -- (-2.8, 1.2) -- cycle;
        \node[state] (v) {$Y$};
        \node[state] (p1)[above left = of v] {$V_1$};
        \node[state] (p2)[above right = of v] {$V_2$};
        \node[state] (c1)[below left = of v] {$V_5$};
        \node[state] (c2)[below right = of v] {$V_4$};
        \node[state] (x)[above right = of c2] {$V_3$};
        \path (p1) edge (v) (p2) edge (v) (v) edge (c1) (v) edge (c2) (x) edge (c2) (c2) edge (c1);
        \end{tikzpicture}}
    \scalebox{.4}{
        \begin {tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}, -latex ,auto ,node distance =2 cm and 2 cm ,on grid , semithick, state/.style ={ circle, draw, minimum width =1 cm}]
        \filldraw[color = green, fill=green!4](-2.8, 2.8) rectangle (2.8, -.8);
        \node[state] (v) {$Y$};
        \node[state, line width=3pt] (p1)[above left = of v] {$V_1$};
        \node[state, line width=3pt] (p2)[above right = of v] {$V_2$};
        \node[state] (c1)[below left = of v] {$V_5$};
        \node[state] (c2)[below right = of v] {$V_4$};
        \node[state, line width=3pt] (x)[above right = of c2] {$V_3$};
        \path (p1) edge (v) (p2) edge (v) (v) edge (c1) (v) edge (c2) (x) edge (c2) (c2) edge (c1);
        \end{tikzpicture}}
    \scalebox{.4}{
        \begin {tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}, -latex ,auto ,node distance =2 cm and 2 cm ,on grid , semithick, state/.style ={ circle, draw, minimum width =1 cm}]
        \filldraw[color = blue, fill=blue!4](-.8, .8) rectangle (4.8, -2.8);
        \node[state, line width=3pt] (v) {$Y$};
        \node[state, line width=3pt] (p1)[above left = of v] {$V_1$};
        \node[state, line width=3pt] (p2)[above right = of v] {$V_2$};
        \node[state] (c1)[below left = of v] {$V_5$};
        \node[state] (c2)[below right = of v] {$V_4$};
        \node[state, line width=3pt] (x)[above right = of c2] {$V_3$};
        \path (p1) edge (v) (p2) edge (v) (v) edge (c1) (v) edge (c2) (x) edge (c2) (c2) edge (c1);
        \end{tikzpicture}}
    \scalebox{.4}{
        \begin {tikzpicture}[baseline={([yshift=-.5ex]current bounding box.center)}, -latex ,auto ,node distance =2 cm and 2 cm ,on grid , semithick, state/.style ={ circle, draw, minimum width =1 cm}]
        \filldraw[color = orange, fill=orange!4](-2.8, .8) rectangle (2.8, -2.8);
        \node[state, line width=3pt] (v) {$Y$};
        \node[state, line width=3pt] (p1)[above left = of v] {$V_1$};
        \node[state, line width=3pt] (p2)[above right = of v] {$V_2$};
        \node[state] (c1)[below left = of v] {$V_5$};
        \node[state, line width=3pt] (c2)[below right = of v] {$V_4$};
        \node[state, line width=3pt] (x)[above right = of c2] {$V_3$};
        \path (p1) edge (v) (p2) edge (v) (v) edge (c1) (v) edge (c2) (x) edge (c2) (c2) edge (c1);
        \end{tikzpicture}}
    \caption{We can decompose $\mathcal{P}_u(v_1,v_2,  v_3, y,  v_4, v_5) = \color{red}\mathcal{P}_u(v_1, v_2, v_3)\color{black} \color{green!60!black}\mathcal{P}_u(y \mid v_1, v_2)\color{black}\color{blue}\mathcal{P}_u(v_4 \mid y, v_3)\color{black}  \color{orange!70!black} \mathcal{P}_u(v_5 \mid y, v_4) \color{black}$.  This figure is restricted to $\Mb(Y)\cup \Set{Y}$. $U$ and any other variables in the graph are omitted for clarity.
    }\label{fig:MBdecomp}
\end{figure}

\begin{algorithm}[H]
    \caption{Bayesian Unzipping}\label{alg:bayesian-unzipping}
    \begin{algorithmic}[1]
        \Require A collection of runs $\mathcal{A}$ of size at most $2^{O(\Delta^2)}$ and their aligned output.  For each $V_i$ and assignment to its parents $\pa(V_i)$ there must be some run with $V_i$ in its independent set with parents conditioned to $\pa(V_i)$ .
        \Ensure $\widetilde{\mathcal{P}}_u(Y \mid \Pa(Y))$
   
        \State Fix a topological ordering on the vertices in $\VV$, $\Seq{X_1,X_2,\dotsc,X_n}$.
        \For{$i=n,n-1,\dotsc,1$}
        \For{each assignment $\pa(X_i)$ to $\Pa(X_i)$}
        \State Let $a$ be a run in $\mathcal{A}$ with $\pa^a(X_i) = \pa(X_i)$.
        \For{$u=1,\dotsc,k$}
        \For{$b=0,1$}

        \If{$X_i$ is a bottom vertex}
        \State Set $\widetilde{\mathcal{P}}_u(x_i^b\mid \pa^a(X_i))\gets \widetilde{\mathcal{P}}_u^a(x_i^b)$.
        \Else
        \State Set $\widetilde{\mathcal{P}}_u(x_i^b\mid \pa(X_i)) \gets \frac{\widetilde{\mathcal{P}}_u^a(x_i^b)\widetilde{\mathcal{P}}_u(\ch^a(X_i)\mid \ttop^a(X_i),x_i^{1-b})}{\widetilde{\mathcal{P}}_u^a(x_i^b)\widetilde{\mathcal{P}}_u(\ch^a(X_i)\mid \ttop^a(X_i),x_i^{1-b})+\widetilde{\mathcal{P}}_u^a(x_i^{1-b})\widetilde{\mathcal{P}}_u(\ch^a(X_i)\mid \ttop^a(X_i),x_i^b)}$.
        \EndIf
        \EndFor
        \EndFor
        \EndFor
        \EndFor
    \end{algorithmic}
\end{algorithm}

\subsubsection{Recovering the distribution on sources} \label{sec:recoversourceprobability} Now consider some arbitrary run $a$ of the $k$-MixProd algorithm with conditioning $\mcond^a$. Since
$\mathcal{P}_u(\VV) = \prod_{V\in \VV} \mathcal{P}_u(V\mid \Pa(V))$, knowing $\mathcal{P}_u(V\mid \Pa(V))$ gives us access to all of $\mathcal{P}_u(\VV)$ after Bayesian unzipping. Thus, we can compute $\mathcal{P}_u(\mcond^a) = \mathcal{P}(\mcond^a\mid u)$. The $k$-MixProd algorithm will return $\mathcal{P}^a(U) = \mathcal{P}(U\mid \mcond^a)$ when run on $a$ (after source alignment).  Finally,  $\mathcal{P}(\mcond^a)$ is directly observable. Combining these terms in Bayes' rule lets us compute the distribution on $U$:
\begin{equation*}
    \mathcal{P}(u) = \frac{\mathcal{P}(u \mid \mcond^a) \mathcal{P}(\mcond^a)}{\mathcal{P}(\mcond^a \mid u)}
\end{equation*}
This equation is well defined under the assumption of positivity.

\subsection{A good collection of runs}
We must also be careful to ensure all of these runs can have their sources aligned and that all the desired parameters can be obtained via Bayesian unzipping.

\begin{definition}[Runs aligned at $X$/alignable pair of runs] A pair of runs $a,b$ are \emph{aligned at $X \in \VV$} if the distributions on $\mathcal{I}^a\mid \mcond^a$ and $\mathcal{I}^b \mid \mcond^b$ are aligned mixture distributions with alignment variable $X\mid \mcond^a$.
    We say that $a,b$ are \emph{alignable} if there is any variable aligning them. \label{def:alignment}
\end{definition}

\begin{observation}
    Two runs $a,b$ are aligned at $X\in \VV$ if and only if
    \begin{enumerate}
        \item $X \in \mathcal{I}^a \cap \mathcal{I}^b$,
        \item $\mb^a(X) = \mb^b(X)$, i.e, $f^a(\Mb(X)) = f^b(\Mb(X))$, and
        \item $X$ is separated given $\mb^a(X)$ (equivalently, given $\mb^b(X)$).
    \end{enumerate}
\end{observation}

\begin{definition} A collection of runs $\mathcal{A}$
    \begin{itemize}
        \item is \emph{alignable} if the \emph{undirected alignment graph} $H(\mathcal{A}) = (\mathcal{A}, E)$ with edges given by $E=\Set{(a,b):\text{$a$ and $b$ are alignable}}$ has a single connected component
        \item \emph{covers $X \in \VV$} if for every assignment $\pa(X)$ to $\Pa(X)$ there exists a run $a\in \mathcal{A}$ with $X \in \mathcal{I}^a$ and $\pa(X) = \pa^a(X)$.
    \end{itemize}
\end{definition}

Our algorithm will require \emph{good} collection of runs:
\begin{definition}[A good collection of runs]
    A collection of well-formed runs $\mathcal{A}$ is \emph{good} if it is \begin{itemize}
        \item alignable,
        \item every run is $(3k-3)$-independent,
        \item the collection covers every vertex in $\VV$,
        \item no vertex appears both as a bottom vertex and a non-bottom vertex, and
        \item every non-bottom vertex has depth at most $9k$.\footnote{The choice of $9k$ can be tightened to $9k-9$, 
        which we suppress.}
    \end{itemize}
\end{definition}

\begin{observation} In a good collection of runs, we can partition all vertices into bottom vertices and non-bottom vertices for the entire collection.
\end{observation}

\begin{definition}[Centers, Central Runs] A set of vertices $\XX = \Set{X_1,\dotsc, X_{3k-3}}\subseteq \VV$ will be called \emph{centers} if the Markov boundaries of the vertices in $\XX$ are disjoint. Given a set of centers $\XX$, a run $a$ is called a \emph{central run} if $\mathcal{I}^a = \XX$.
\end{definition}

The main result on good collections of runs that we'll utilize is the following:
\begin{lemma} Given a graph with max degree $\Delta$ satisfying $n \geq (3k-3)\times O(\Delta^4)$, we can find a set of centers $\XX = \Set{X_1,\dotsc, X_{(3k-3)}} \subseteq \VV$ of size $(3k-3)$ and depth at most $9k$, such that by running Algorithm~\ref{alg:good runs}, we obtain a good collection of runs $\mathcal{A}$ of size $O(2^{\Delta^2}n)$ assuming that the mixture $\mathcal{P}$ is $\zeta$-informative. \label{lem:main result on good collection}
\end{lemma}

The proof of this lemma is in Appendix~\ref{apx:good collections}, along with other sufficient conditions for the existence of good collections of runs.
\subsection{Outline of the algorithm}
Our algorithm will take in a good collection of runs $\mathcal{A}$ satisfying the size constraint obtained by Lemma~\ref{lem:main result on good collection}. While any good collection of runs will suffice for the correctness of our algorithm for sufficently large sample complexity, having this bound on the collection size ensures that our sample complexity bounds are sufficient.

Essentially our algorithm has four steps:

\begin{enumerate}
    \item Use a $k$-MixProd algorithm on $\mathcal{P}(\mathcal{I}^a\mid \MCond^a)$ for each run $a\in \mathcal{A}$ to compute $P(Y\mid \Mb(Y), U_a = u)$ for all variables $Y\in \VV$.
    \item Align the parameters obtained from the previous step to ensure that $U$ means the same thing across different runs, giving $\mathcal{P}_u(Y\mid \Mb(Y))$.
    \item Work backwards in topological order to recover $\mathcal{P}_u(Y\mid \Pa(Y))$ for each vertex $Y\in \VV$ via Bayesian unzipping.
    \item Compute $\mathcal{P}(U)$ by applying Bayes' law.
\end{enumerate}

The full procedure appears as Algorithm~\ref{alg:k-bnd}.

\begin{algorithm}[H]
    \caption{The $k$-MixBND algorithm}\label{alg:k-bnd}
    \begin{algorithmic}[1]
        \Require A good collection of runs $\mathcal{A}$ of size at most $2^{O(\Delta^2)}$; a target accuracy $\varepsilon$.
        \Ensure $\widetilde{\mathcal{P}}_u(Y \mid \Pa(Y))$ and $\widetilde{\mathcal{P}}(U)$.
        \State Estimate $\widetilde{\mathcal{P}}^a(\mathcal{I}^a)$ for all runs $a\in \mathcal{A}$ using \[n \log n \cdot \varepsilon^{-2} (\Delta+1)^{O(k)}2^{O(\Delta^2 + \Delta k)}  (\min \pi_i)^{-O(\log k)} (1/\zeta)^{O(k \log k + \Delta^2 k)}\] samples from $\mathcal{P}$.
        \State Set $\varepsilon' \gets \varepsilon/(6(\Delta+1))^{18k}$.
        \For{each run $a\in \mathcal{A}$}
        \State Let $Q \coloneqq \mathcal{P}^a(\mathcal{I}^a)$.
        \State Set $\Paren{[m^a_{iu}]_{i\in \mathcal{I}^a,u\in [k]},(\pi^a_{u})_{u\in [k]}}\gets \textsc{LearnProductMixture}(Q, \varepsilon')$, where $m^a_{ij} = Q_u(x_i^1) = P_u(x_i^1\mid \mcond^a)$ and $\pi^a_{u} = P(U=u\mid \mcond^a)$
        \EndFor
        \State Run Algorithm~\ref{alg:align} to align the sources in the output of all the runs in $\mathcal{A}$.
       \State Run Algorithm~\ref{alg:bayesian-unzipping} to unzip the parameters.
        \State Fix any run $a\in \mathcal{A}$.
        \For{$j=1,2,\dotsc,k$}
        \State Set $\widetilde{\mathcal{P}}(u^j) \gets \frac{\widetilde{\mathcal{P}}(u^j \mid \mcond^a)\widetilde{\mathcal{P}}(\mcond^a)}{\widetilde{\mathcal{P}}(\mcond^a\mid u^j)}$.
        \EndFor
    \end{algorithmic}
\end{algorithm}

\subsection{Statement of main result for $k$-MixBND}
\begin{theorem}\label{thm: main k-BND}
If $\mathcal{P}(\VV) =\sum_{u}\mathcal{P}_u(\VV)\mathcal{P}(u)$ is a $\zeta$-informative mixture of $k$ distributions over a DAG $G=(\VV, E)$ and there exists a set of $3k-3$ centers, then there is an algorithm to compute estimates of all parameters to accuracy $\varepsilon$ with sample complexity (and runtime) bounded by  
\[n \log n \cdot \varepsilon^{-2} (\Delta+1)^{O(k)}2^{O(\Delta^2 + \Delta k)}  (\min \pi_i)^{-O(\log k)} (1/\zeta)^{O(k \log k + \Delta^2 k)}.\]
\end{theorem}

\section{Analyzing the $k$-MixBND algorithm}

\subsection{Applying the $k$-MixProd algorithm on each run}
We must collect the moments needed for each application of the $k$-MixProd algorithm and the accuracy with which the moments are estimated will show up in the output of the $k$-MixProd algorithm:

\begin{lemma} Suppose that we fix an $\varepsilon$ target accuracy and that $\mathcal{P}$ is $\zeta$-informative. Then with total sample complexity \[(n \log n/\varepsilon^2) (\Delta+1)^{O(k)}2^{O(\Delta^2 + \Delta k)}  (1/\min \pi_i)^{O(\log k)} (1/\zeta)^{O(k \log k + \Delta^2 k)},\] we can ensure that we have estimates of $\widetilde{\mathcal{P}}^a(\mathcal{I}^a)$ sufficiently accurate so that $\textsc{LearnProductMixture}(\mathcal{P}^a(\mathcal{I}^a), \varepsilon')$ (where $\varepsilon' = \varepsilon\zeta/(6(\Delta+1))^{9k}$) will return $\widetilde{\mathcal{P}}_u^a(\mathcal{I}^a)$ and $\widetilde{\mathcal{P}}^a(u)$ satisfying
    \[ \Abs{\widetilde{\mathcal{P}}_u^a(y^b \mid \pa^a(Y)) - \mathcal{P}_u^a(y^b\mid \pa^a(Y))} \leq \varepsilon\zeta/(6(\Delta+1))^{9k}\quad\text{and}\quad \Abs{\widetilde{\mathcal{P}}^a(u) - \mathcal{P}^a(u)} \leq \varepsilon\zeta/(6(\Delta+1))^{9k} \]
    for all $Y\in \mathcal{I}^a$, $b\in \Set{0,1}$, $u\in [k]$ in time $(6(\Delta+1))^{9k}\varepsilon^{-2}(1/\min \pi_i)^{O(k\log k)}(1/\zeta)^{O(k\log k)}n^{O(k)}$.
    Consequently, we will also have
    \[ \frac{\Abs{\widetilde{\mathcal{P}}_u^a(y^b \mid \pa^a(Y)) - \mathcal{P}_u^a(y^b\mid \pa^a(Y))}}{\mathcal{P}_u^a(y^b\mid \pa^a)} \leq \varepsilon/(6(\Delta+1))^{9k}\quad\text{and}\quad\frac{\Abs{\widetilde{\mathcal{P}}^a(u) - \mathcal{P}^a(u)}}{\mathcal{P}^a(u)} \leq \varepsilon/(6(\Delta+1))^{9k}.\]
\end{lemma}
\begin{proof} With the given sample complexity we have that the estimates of $\widetilde{\mathcal{P}}(\mathcal{I}^a)$ are within accuracy $\varepsilon\zeta/(6(\Delta+1))^{9k}$ for each run $a\in \mathcal{A}$ with high probability. The conclusion follows from the first part of Theorem~\ref{thm: main} on the behavior of the algorithm $\textsc{LearnProductMixture}$. Note that we require Lemma~\ref{lem:probabilities not too small} in order for the expected number of samples to observe one sample with $\cond^a = \Cond^a$ to be bounded above by $\zeta^{O(\Delta^2k)}$ for any run  $a\in \mathcal{A}$. 
\end{proof}

\subsection{Analysis of alignment}
\begin{lemma} Suppose $\mathcal{P}(V)$ is a $\zeta$-informative mixture of $k$ Bayesian networks.
    Given an alignable collection of runs $\mathcal{A}$ and the outputs of a $k$-MixProd algorithm $(\Paren{m^a_{iu}}_u, \Paren{\pi^a_u}_u)$ for each run $a\in \mathcal{A}$ satisfying
    \[ \Abs{\pi^a_u - \mathcal{P}^a(u)} < \zeta /2,\quad \Abs{m^a_{iu} - \mathcal{P}^a_u(x_i^1)} < \zeta / 2 \quad \text{for all }u\in [k], a\in \mathcal{A} \]  Algorithm~\ref{alg:align} returns $\mathcal{P}^a_u(X_i=1)$ for $i \in [n]$ and $ \mathcal{P}^a(u)$. \label{lem:alignment}
\end{lemma}

\begin{proof}
    There is a unique permutation $\sigma^{i\to j}$ such that
    \[\Norm{\mathcal{P}^{(i)}(X^{(i)}_{\AV} \mid \sigma^{i \to j}(u^{(i)})) - \mathcal{P}^{(j)}(X^{(j)}_{\AV} \mid u^{(j)})}_{\infty} < \frac{\zeta}{2}
    \]
    The algorithm chooses this permutation and aligns accordingly.
\end{proof}

\subsection{Recovering Parameters} \label{sec:recovering parameters}

\begin{lemma} Let $a$ be a run with $Y\in \mathcal{I}^a$ and let $\mathcal{P}$ be $\zeta$-informative for $\zeta > 0$. Then we can compute $\mathcal{P}_u(y^b\mid \pa^a(Y))$ for $b\in \Set{0,1}$ as follows:
    \begin{enumerate}
        \item If $Y\in \mathcal{J}_a$, then \[ \mathcal{P}_u(y^b\mid \pa^a(Y)) = \mathcal{P}_u^a(y^b). \]
        \item If $Y\notin \mathcal{J}_a$, then \begin{equation}
                  \mathcal{P}_u(y^b\mid \pa^a(y)) = \frac{\mathcal{P}_u^a(y^b)\mathcal{P}_u(\ch^a(y)\mid \ttop^a(y),y^{1-b})}{\mathcal{P}_u^a(y^b)\mathcal{P}_u(\ch^a(y)\mid \ttop^a(y),y^{1-b})+\mathcal{P}_u^a(y^{1-b})\mathcal{P}_u(\ch^a(y)\mid \ttop^a(y),y^b)}\label{eq:parameter recovery}\end{equation}
    \end{enumerate} \label{lem:recovery equations}
\end{lemma}
\begin{proof}
    In the following we'll fix $\ch^a \coloneqq \ch^a(Y)$, $\mb^a \coloneqq \mb^a(Y)$, $\pa^a \coloneqq \pa^a(Y)$, and $\ttop^a \coloneqq \ttop^a(Y)$.
    If $Y\in \mathcal{J}_a$, then $\mathcal{P}_u^a(y^b) = \mathcal{P}_u(y^b\mid \pa^a(Y))$, so the claim is trivially true.
    If $Y\notin \mathcal{J}_a$, then
    \begin{align*}
        \mathcal{P}_u^a(y^b)  = \mathcal{P}_u^a(y^b\mid \mb^a)
        = \frac{\mathcal{P}_u(y^b,\mb^a)}{\mathcal{P}_u(\mb^a)}
        = \frac{\mathcal{P}_u(y^b,\mb^a)}{\sum_{b' = 0,1}\widetilde{\mathcal{P}}_u(y^{b'}, \mb^a)}
    \end{align*} and we can expand $\mathcal{P}_u(y^b, \mb^a)$ using Observation~\ref{obs:top decomposition} as
    \[ \mathcal{P}_u(y^b, \mb^a) = \mathcal{P}_u(y^b\mid \pa^a)\mathcal{P}_u(\ch^a\mid \ttop^a, y^b) \mathcal{P}_u(\ttop^a). \] Substituting this back into the preceding equation, we obtain
    \begin{align*}
        \mathcal{P}_u^a(y^b)
         & =  \frac{\mathcal{P}_u(y^b\mid \pa^a)\mathcal{P}_u(\ch^a\mid \ttop^a,y^b)\mathcal{P}_u(\ttop^a)}{\sum_{b'\in \Set{0,1}}\mathcal{P}_u(y^{b'}\mid \pa^a(Y))\mathcal{P}_u(\ch^a\mid \ttop^a, y^{b'}) \mathcal{P}_u(\ttop^a)} \\
         & = \frac{\mathcal{P}_u(y^b\mid \pa^a)\mathcal{P}(\ch^a\mid \ttop^a,y^b)}{\sum_{b'\in \Set{0,1}}\mathcal{P}_u(y^{b'}\mid \pa^a)\mathcal{P}_u(\ch^a\mid \ttop^a y^{b'})}.
    \end{align*}

    We now multiply both sides by the denominator (which is non-zero since all terms are strictly positive by assumption) and then simplify to obtain
    \[ \mathcal{P}_u(y^b\mid \pa^a)\Paren{\mathcal{P}_u^a(y^{1-b})\mathcal{P}(\ch^a\mid \ttop^a,y^b)} + \mathcal{P}_u(y^{b}\mid \pa^a)\Paren{-\mathcal{P}_u^a(y^{1-b})\mathcal{P}_u(\ch^a\mid \ttop^a,y^{1-b})} = 0. \]
    When augmented with the equation \[ \mathcal{P}_u(y^b\mid \pa^a) + \mathcal{P}_u(y^{1-b}\mid \pa^a) = 1\] we have a system of two equations in $\mathcal{P}_u(y^b\mid \pa^a),\mathcal{P}_u(y^{1-b}\mid \pa^a)$ which is non-singular whenever
    \[ \Paren{\mathcal{P}_u^a(y^{1-b})\mathcal{P}(\ch^a\mid \ttop^a,y^b)} \neq \Paren{-\mathcal{P}_u^a(y^{1-b})\mathcal{P}_u(\ch^a\mid \ttop^a,y^{1-b})}.\] Since all probabilities occurring are strictly positive, this will always be the case and we can solve the system for $\mathcal{P}_u(y^b\mid \pa^a)$. The resulting equation is
    \[ \mathcal{P}_u(y^b\mid \pa^a) = \frac{\mathcal{P}_u^a(y^b)\mathcal{P}(\ch^a\mid \ttop^a,y^{1-b})}{\mathcal{P}_u^a(y^b)\mathcal{P}(\ch^a\mid \ttop^a,y^{1-b})+\mathcal{P}_u^a(y^{1-b})\mathcal{P}(\ch^a\mid \ttop^a,y^{b})}, \] proving the claim.
\end{proof}

The following standard inequalities will be useful throughout the analysis:
\begin{observation}
    \begin{align*} \frac{1+x}{1-x} & \leq 1+3x\quad\text{and}    & \quad 1-3x    & \leq \frac{1-x}{1+x} \quad & \text{for $x\in (0, 1/4)$;}                 \\
               1-2rx           & \leq (1-x)^r\quad\text{and} & \quad (1+x)^r & \leq 1+2rx          \quad  & \text{for $x\in (0,1),r\geq 1, rx \leq 1$.}
    \end{align*} \label{obs:inequalities}
\end{observation}

\begin{lemma}\label{lem:canrecover_inductive} Given access to $\Set{\widetilde{\mathcal{P}}_u^a(Y)}_{a\in \mathcal{A}, j\in [k]}$ for a good collection of runs $\mathcal{A}$ from a $\zeta$-informative distribution $\mathcal{P}$ satisfying $\Abs{\widetilde{\mathcal{P}}_u^a(y^b) - \mathcal{P}_u^a(y^b)}/\mathcal{P}_u^a(y^b) \leq \varepsilon$ for all $Y\in \VV$, $b\in \Set{0,1}$, $u\in [k]$ for $\varepsilon$ sufficiently small (and the estimated probabilities used as input to $k$-MixProd algorithms), the procedure in Algorithm~\ref{alg:bayesian-unzipping} will output $\widetilde{\mathcal{P}}_{u}(y^b\mid \Pa(Y))$ and $\widetilde{\mathcal{P}}(u)$ for all $Y\in \VV$, $b\in \Set{0,1}, u\in [k]$ satisfying
    \[ \frac{\Abs{\widetilde{\mathcal{P}}_u(y^b\mid \pa(Y)) - \mathcal{P}_u(y^b\mid \pa(Y))}}{\mathcal{P}_u(y^b\mid \pa(Y))} \leq (6(\Delta + 1))^{\ell}\varepsilon \] where $\ell$ is the distance from $Y$ to the nearest bottom vertex if $Y$ doesn't appear as a bottom vertex and is $0$ if $Y$ is a bottom vertex. \end{lemma}

\begin{proof}
    We fix a topological ordering on $\mathcal{V}$, let's say $X_1,X_2,\dotsc, X_n$. Now starting with the last vertex in topological order, $X_n$, and proceeding in decreasing order, we'll compute $\widetilde{\mathcal{P}}_u(X_i\mid \Pa(X_i))$. For bottom vertices, we set $\widetilde{\mathcal{P}}_u(x_i^b\mid \pa(X_i)) = \widetilde{\mathcal{P}}_u^a(x_i^b)$ for some run $a$ with $\pa(X_i) = \pa^a(X_i)$. It immediately follows that $P_u(X_i\mid \Pa(X_i))$ satisfies the desired relative error bound. Inductively, assume we've already computed $\widetilde{\mathcal{P}}_u(V_m\mid \Pa(V_m))$ for all $m > i$ satisfying the stated bounds, and that we also have access to $\widetilde{\mathcal{P}}^a_u(V_i) = \widetilde{\mathcal{P}}_u(V_i \mid \MCond^a)$ for a subset of the runs $a \in \mathcal{A}$ that cover $V_i$. To streamline notation, let $Y \coloneqq V_i$. Now fix a run $a\in \mathcal{A}$ with $Y\in \mathcal{I}^a$.

    Now we can bound each term above and below using the inductive hypothesis to get the desired result. In particular, we know that
    \begin{equation*}
        (1-\varepsilon)\mathcal{P}_u^a(y^b)                                        \leq \widetilde{\mathcal{P}}_u^a(y^b)                \leq (1+\varepsilon)\mathcal{P}_u^a(y^b)                                                                        \end{equation*} for $b\in \Set{0,1}$ and \begin{equation*}
        (1-(6(\Delta + 1))^{\ell-1}\varepsilon)\mathcal{P}_u(v^b\mid \pa^a(v^b)) \leq \widetilde{\mathcal{P}}_u(v^b\mid \pa^a(v^b))  \leq (1+(6(\Delta + 1))^{\ell-1}\varepsilon) \mathcal{P}_u(v^b\mid \pa^a(v^b))\end{equation*}
    for all $V\in \Ch(Y)$ and $b\in \Set{0,1}$ which implies that both the numerator and denominator of \eqref{eq:parameter recovery} are within a factor of $(1\pm \varepsilon)(1\pm (6(\Delta+1))^{\ell-1} \varepsilon)^\Delta$ of the correct value for those terms. It immediately follows that $\widetilde{\mathcal{P}}_u(y^b\mid \pa^a(Y))$ is bounded by
    \begin{align*}\mathcal{P}_u(y^b \mid \pa^a(Y)) \frac{1-\varepsilon}{1+\varepsilon} \Paren{\frac{1-(6(\Delta+1))^{\ell -1}\varepsilon}{1+(6(\Delta+1))^{\ell -1}\varepsilon}} &\leq \widetilde{\mathcal{P}}_u(y^1 \mid \pa^a(Y)) \\&\leq \mathcal{P}_u(y^1 \mid \pa^a(Y)) \frac{1+\varepsilon}{1-\varepsilon} \Paren{\frac{1+(6(\Delta+1))^{\ell -1}\varepsilon}{1-(6(\Delta+1))^{\ell -1}\varepsilon}}. \end{align*}

    We now use the inequalities from Observation~\ref{obs:inequalities} to simplify the bounds as follows:
    \begin{align*}\widetilde{\mathcal{P}}_u(y^b \mid \pa^a(Y))
         & \leq \mathcal{P}_u(y^b \mid \pa^a(Y)) \frac{1+\varepsilon}{1-\varepsilon} \Paren{\frac{1+(6(\Delta+1))^{\ell -1}\varepsilon}{1-(6(\Delta+1))^{\ell -1}\varepsilon}}^{\Delta} \\
         & \leq \mathcal{P}_u(y^b \mid \pa^a(Y))(1+3\varepsilon)(1+3(6(\Delta+1))^{\ell -1}\varepsilon)^{\Delta}                                                          \\
         & \leq \mathcal{P}_u(y^b \mid \pa^a(Y))(1+3(6(\Delta+1))^{\ell -1}\varepsilon)^{\Delta+1}                                                                 \\
         & \leq \mathcal{P}_u(y^b \mid \pa^a(Y))(1+2\cdot 3(\Delta + 1)(6(\Delta+1))^{\ell -1}\varepsilon)                                                         \\
         & \leq \mathcal{P}_u(y^b \mid \pa^a(Y))(1+(6(\Delta+1))^{\ell}\varepsilon).                                                                               \\
    \end{align*} The lower bound is analogous.
\end{proof}

\begin{lemma}\label{lem:k-bnd final error blowup} Given access to $\Set{\widetilde{\mathcal{P}}_u^a(Y)}_{a\in \mathcal{A}, j\in [k]}$ for a good collection of runs $\mathcal{A}$ from a $\zeta$-informative distribution $\mathcal{P}$ satisfying $\Abs{\widetilde{\mathcal{P}}_u^a(y^b) - \mathcal{P}_u^a(y^b)}/\mathcal{P}_u^a(y^b) \leq \varepsilon$ for all $Y\in \mathcal{I}^a$, $b\in \Set{0,1}$, $u\in [k]$ for $\varepsilon$ sufficiently small, the procedure in Algorithm~\ref{alg:bayesian-unzipping} will output $\widetilde{\mathcal{P}}_{u}(y^b\mid \Pa(Y))$ and $\widetilde{\mathcal{P}}(u)$ for all $Y\in \VV$, $b\in \Set{0,1}, u\in [k]$ satisfying
    \[ \frac{\Abs{\widetilde{\mathcal{P}}_u(y^b\mid \pa(Y)) - \mathcal{P}_u(y^b\mid \pa(Y))}}{\mathcal{P}_u(y^b\mid \pa(Y))} \leq (6(\Delta+1))^{9k}\varepsilon \] and \[ \frac{\Abs{\widetilde{\mathcal{P}}(u^j) - \mathcal{P}(u^j)}}{\mathcal{P}(u^j)} \leq 5\Delta^2\varepsilon. \]

    Since $\mathcal{P}_u(y^b \mid \pa(Y)), \mathcal{P}(u^j) \leq 1$, we also have that \[ \Abs{\widetilde{\mathcal{P}}_u(y^b\mid \pa(Y)) - \mathcal{P}_u(y^b\mid \pa(Y))} \leq (6(\Delta+1))^{9k}\varepsilon\quad \text{and}\quad\Abs{\widetilde{\mathcal{P}}(u^j) - \mathcal{P}(u^j)} \leq 5\Delta^2\varepsilon. \]\end{lemma}
\begin{proof}
    The error bound on $\widetilde{\mathcal{P}}_u(y^b\mid \pa(Y))$ follows from Lemma~\ref{lem:canrecover_inductive} above and the depth bound of $9k$ on non-bottom vertices.

    For the error bound on $\widetilde{\mathcal{P}}(u^j)$ we write
    \begin{align*} \widetilde{\mathcal{P}}(u^j) = \frac{\widetilde{\mathcal{P}}(u^j\mid \mcond^a)\widetilde{\mathcal{P}}(\mcond^a)}{\widetilde{\mathcal{P}}(\mcond^a\mid u^j)}
    \end{align*} and analyze each term separately. First, define $Z \coloneqq \An(\MCond^a)\setminus \MCond^a$ to be all ancestors of vertices conditioned upon in $a$ minus $\MCond^a$. Fix a topological order on the vertices in $Z \cup \MCond^a$, $\Seq{X_1,\dotsc, X_m}$. We can write
    \begin{align*}
        \widetilde{\mathcal{P}}(\mcond^a\mid u^j)
         & = \sum_{z\in \Set{0,1}^{Z}}\prod_{i=1}^m \widetilde{\mathcal{P}}(z(X_i)\mid z(X_1,\dotsc, X_{i-1}),\mcond^a)                                                                                 \\
         & = \sum_{z\in \Set{0,1}^{Z}}\prod_{X_i \in Z\cup \MCond^a} \widetilde{\mathcal{P}}(z(X_i)\mid z(\Pa(X_i)), \pa^a(X_i))                                                                        \\
         & = \sum_{z\in \Set{0,1}^{Z}}\Paren{\prod_{X_i \in Z} \widetilde{\mathcal{P}}(z(X_i)\mid z(\Pa(X_i)), \pa^a(X_i))}\Paren{\prod_{X_i \in \MCond^a} \widetilde{\mathcal{P}}(z(X_i)\mid z(\Pa(X_i)), \pa^a(X_i))}.
    \end{align*} Now let $C_z$ denote the value in the first grouped product for a given $z\in \Set{0,1}^Z$. Since $\sum_{z\in \Set{0,1}^Z} C_z = 1$, we can bound the error in the result by the error in the second grouped product:
    \[ (1-\varepsilon)^{\Card{\MCond^a}}\mathcal{P}(\mcond^a\mid u^j) \leq \widetilde{\mathcal{P}}(\mcond^a\mid u^j) \leq (1+\varepsilon)^{\Card{\MCond^a}}\mathcal{P}(\mcond^a\mid u^j). \] Finally using the bound $\Card{\Cond^a} < 2\Delta^2$ and the inequalities from Observation~\ref{obs:inequalities} we obtain
    \[ (1-2\Delta^2\varepsilon)\mathcal{P}(\mcond^a\mid u^j) \leq \widetilde{\mathcal{P}}(\mcond^a\mid u^j) \leq (1+2\Delta^2\varepsilon)\mathcal{P}(\mcond^a\mid u^j). \] By assumption $\widetilde{\mathcal{P}}(u^j\mid \cond^a) \in (1\pm \varepsilon)\mathcal{P}(u^j\mid \cond^a)$; $\widetilde{\mathcal{P}}(\cond^a)$ is also known up to multiplicative accuracy $\varepsilon$ since the sampling needed to estimate the distributions that are input to $k$-MixProd algorithms far exceed that needed to get an $\varepsilon$ estimate of this quantity.
    Thus, the resulting bound on $\widetilde{\mathcal{P}}(u^j)$ is
    \[ (1-2\Delta^2\varepsilon)(1-3\varepsilon)\mathcal{P}(u^j) \leq \widetilde{\mathcal{P}}(\mcond^a\mid u^j) \leq (1+2\Delta^2\varepsilon)(1+3\varepsilon)\mathcal{P}(\mcond^a\mid u^j) \] which can be simplified to
    \[ (1-2\Delta^2\varepsilon-3\varepsilon)\mathcal{P}(u^j) \leq \widetilde{\mathcal{P}}(\mcond^a\mid u^j) \leq (1+2\Delta^2\varepsilon+3\varepsilon)\mathcal{P}(\mcond^a\mid u^j) \] using $(1-x)(1-y) \geq (1-x-y)$ for $x,y\in [0,1)$ and $2\Delta^2\varepsilon + 3\varepsilon \leq 5\Delta^2\varepsilon$ for $\Delta \geq 1$.
\end{proof}

\newpage


