\section{Introduction}
Many learning and inference tasks in high-dimensional statistics boil down to estimating / approximating the joint probability of a set of random variables (RVs).
Having the joint probability allows us to construct predictors that are optimal in a certain probabilistic sense,
e.g., the \textit{maximum a posteriori} (MAP) estimator or the \textit{minimum mean squared error} (MMSE) predictor. However, in the high-dimensional regime, directly estimating the joint probability via `structure-free' methods such as sample averaging is considered infeasible \cite{wainwright2019high}. Suppose that there are $N$ RVs, where each has an $I$-dimensional alphabet. To estimate the joint probability reliably, we generally need $S \gg \Omega(I^{N})$ `diverse' $N$-dimensional examples (without any entries missing). This is because the probability of encountering most $N$-tuples is very low. Therefore, only a small portion of the empirical distribution will be non-zero given a reasonable amount of data samples---this makes the approach considerably inaccurate in most cases.

Many workarounds have been proposed to circumvent this challenge in the literature. For example, linear approximations have been widely used, e.g., the linear MMSE (LMMSE) estimator \cite{kay1993fundamentals}.
Logistic regression, kernel methods, and neural networks can be understood as nonlinear function approximation-based counterparts \cite{bishop2006pattern,murphy2012machine}. These are effective surrogates, but do not directly address the fundamental challenge in estimating high-dimensional joint probability from limited samples. Another commonly adopted route in statistical learning is to make use of problem-specific structural assumptions such as graphical models and prior distributions to help reduce the model complexity \cite{jordan1998learning,tipping2001sparse,blei2003latent}---but making these assumptions restricts the methods to specific applications.

Very recently, Kargas {\it et al.} proposed a new framework for {\it blindly} estimating the \textit{joint probability mass function} (PMF) of $N$ discrete finite-alphabet RVs \cite{kargas2017tensors}. Building upon a link between high-dimensional joint PMFs and low-rank tensors, the work in \cite{kargas2017tensors} shows that any $N$-dimensional joint PMF admits a naive Bayes {\it representation}. In addition, if the RVs are `reasonably dependent', the joint PMF can be recovered via jointly decomposing the three-dimensional marginal PMFs (which are third-order tensors). The approach in \cite{kargas2017tensors} does not involve any structural information among the RVs or function approximation and provably guarantee the recoverability of the joint PMF---if the three-dimensional marginal PMFs are accessible.

The work in \cite{kargas2017tensors} has shown promising results. However, many challenges remain. First, accurately estimating three-dimensional marginal PMFs is not a trivial task, which requires enumerating the co-occurrences of three RVs. This is particularly hard if the alphabets of the RVs are large and the data is sparse.
Second, tensor decomposition is a challenging optimization problem \cite{hillar2013most}. Factoring a large number of latent factor-coupled tensors as in \cite{kargas2017tensors} also gives rise to scalability issues. A natural question is: can we use pairwise marginals, i.e., joint PMFs of {\it two} RVs, to achieve the same goal of joint PMF recovery? Pairwise marginals are much easier to estimate relative to triple ones. In addition, the pairwise marginals give rise to probability mass matrices other than tensors---which may circumvent tensor decomposition in algorithm design and thus lead to more lightweight solutions.

Notably, a recent work in \cite{yeredor2019maximum} offers an {\it expectation maximization} (EM) algorithm that directly estimates the latent factors of the $N$th-order probabilistic tensor from the `raw data', instead of working with a large number of three-dimensional marginals. The EM algorithm is well-motivated, since it is associated with the {\it maximum likelihood estimator} (MLE) under the naive Bayes model in \cite{kargas2017tensors}. In addition, it admits simple and economical updates, and thus is much more scalable relative to the coupled tensor decomposition approach. However, unlike the tensor-based approach in \cite{kargas2017tensors}, it is unclear if the recoverability of the joint PMF can be guaranteed using the EM and the MLE. In addition, since the ML estimation problem is nonconvex, convergence properties of the EM algorithm are unclear.

\smallskip

\noindent
{\bf Contributions.}
In this work, we propose a new framework that utilizes pairwise marginals to recover the joint PMF of an arbitrary number of discrete finite-alphabet RVs.
We address both the theoretical aspects (e.g., recoverability) and the algorithmic aspects. Our contributions are as follows:

\noindent
$\bullet$ {\bf A Pairwise Marginal-based Framework.} On the theory side, we offer in-depth analyses for the recoverability of the joint PMF from pairwise marginals.
We show that the joint PMF can be recovered via jointly factoring the pairwise marginals under nonnegativity constraints in a latent factor-coupled manner.
Unlike \cite{kargas2017tensors} which leverages the uniqueness of tensor decomposition to establish recoverability, our analysis utilizes the model identifiability of \textit{nonnegative matrix factorization} (NMF)---leading to recoverability conditions that have a different flavor compared to those in \cite{kargas2017tensors}.
More importantly, accurately estimating pariwise marginals is much easier relative to the three-dimensional counterparts, in terms of the amount of data samples needed.

\noindent
$\bullet$ {\bf Provable Algorithms.} On the algorithm side, we propose to employ a block coordinate descent (BCD)-based algorithm for the formulated coupled NMF problem.
More importantly, we propose a pragmatic and easy-to-implement initialization approach, which is based on performing a Gram--Schmidt (GS)-like structured algorithm on a carefully constructed `virtual NMF' model.
We show that this approach works provably well even when there is modeling noise---e.g., noise induced by finite-sample estimation of the pairwise marginals. This is in contrast to the method in \cite{kargas2017tensors}, whose recoverability guarantees are based on the assumption that the third-order marginals are perfectly accessible.
We also show that the EM algorithm from \cite{yeredor2019maximum} can provably improve upon proper initialization (e.g., that given by the GS-like algorithm) towards the target latent model parameters that are needed to reconstruct the joint PMF.

\noindent
$\bullet$ {\bf Extensive Validation.} We test the proposed approach on a large variety of synthetic and real-world datasets. In particular, we validate our theory via performing classification on four different UCI datasets and three movie recommendation tasks.

\smallskip

Part of the paper was submitted to the Asilomar signal processing conference 2020 \cite{ibrahim2020recover}. The journal version additionally includes new theoretical results, e.g., the recoverability analysis for the coupled NMF formulation, and a BCD-based CNMF algorithm, the detailed proofs for the Gram--Schmidt-like procedure, and the optimality analysis for the EM algorithm. More real data experiments are also presented.

\smallskip

\noindent
{\bf Notation.} We use $x,\bm x,\bm X,\underline{\boldsymbol{X}}$ to denote a scalar, vector, matrix and  tensor, respectively. $\mathbb{D}_{\rm KL}$ denotes the KL divergence. $\kappa(\bm X)$ denotes the condition number of the matrix $\bm X$ and is given by $\kappa(\bm X) = \frac{\sigma_{\max}(\bm X)}{\sigma_{\min}(\bm X)}$ where $\sigma_{\max}$ and $\sigma_{\min}$ are the largest and smallest singular values of $\bm X$, respectively. $\|\bm x\|_0$ denotes the `zero norm' of the vector $\bm x$, i.e., the number of nonzero elements in $\bm x$. $\bm X \ge \bm 0$ implies that all the elements of $\bm X$ are nonnegative. ${\rm cone}(\bm X)$ represents the conic hull formed by the columns of $\bm X$ and ${\rm conv}(\bm X)$ denotes the convex hull formed by the columns of $\bm X$. $\emptyset$ denotes the empty set. $\|\bm X\|_{2}$ represents the 2-norm of matrix $\bm X$ and $\|\bm X\|_{2} = \sigma_{\max}(\bm X)$.  $\|\bm X\|_{\rm F}$ denotes the Frobenius norm of $\bm X$. $\|\bm x\|_2$ and $\|\bm x\|_1$ denote $\ell_2$ and $\ell_1$ norm of vector $\bm x$, respectively. ${\rm vec}(\bm X)$ denotes the vectorization operator that concatenates the columns of $\bm X$. $^{\!\top\!}$ and $^\dag$ denote transpose and pseudo-inverse, respectively. $|{\cal C}|$ denotes the cardinality of set ${\cal C}$. ${\rm Diag}(\bm x)$ is a diagonal matrix that holds the entries of $\bm x$ as its diagonal elements. `$\odot$' denotes the Khatri-Rao product.

\section{Background}
Consider a set of discrete and finite-alphabet RVs, i.e., $Z_1,\ldots,Z_N$. We will use ${\sf Pr}(i_1,\ldots,i_N)$ as the shorthand notation to represent $${\sf Pr}\left(Z_1=z_1^{(i_1)},\ldots,Z_N=z_N^{(i_N)}\right)$$ in the sequel, where $\{z_n^{(1)},\ldots,z_n^{(I_n)}\}$ denotes the alphabet of $Z_n$.
Suppose that we do not have access to the joint PMF ${\sf Pr}(i_1,\ldots,i_N)$.
Instead, we have access to the lower-dimensional marginal distributions, e.g., ${\sf Pr}(i_j,i_k,i_\ell)$ and ${\sf Pr}(i_j,i_k)$ for different $j,k,\ell$. Kargas \textit{et al.} asked the following research question \cite{kargas2017tensors}:

\vspace{.15cm}
\noindent
\boxed{
	\begin{minipage}{.98\linewidth}
		\noindent
		Can the joint distribution of $Z_1,\ldots,Z_N$ be recovered from the low-dimensional marginals without knowing any structural information about the underlying probabilistic model?
	\end{minipage}
}
\vspace{.15cm}

Note that estimating a high-dimensional joint distribution from lower-dimensional marginals is challenging, but not entirely impossible---if some structural information between the RVs (e.g., tree or Markov chain) is known; see discussions in \cite{kargas2017tensors}.
However, learning the structural information from data {\it per se} is often a hard problem \cite{jordan1998learning}.

\subsection{Tensor Algebra}
The work in \cite{kargas2017tensors} utilizes a connection between joint PMFs and low-rank tensors under the \textit{canonical polyadic decomposition} (CPD) model \cite{harshman1970foundations,sidiropoulos2017tensor}.
If an $N$-way tensor $\underline{{\bm X}} \in \mathbb{R}^{I_1 \times I_2 \times \cdots \times I_N}$ has CP rank $F$,
it can be written as:
\begin{equation}\label{eq:individ}
\underline{{\bm X}}(i_1,i_2,\ldots,i_N) = \sum_{f=1}^F \boldsymbol{\lambda}(f)  \prod_{n=1}^N {\bm A}_n(i_n,f),
\end{equation}
where  ${\bm A}_n \in \mathbb{R}^{I_n \times F}$ is called the mode-$n$ latent factor.
In the above, ${\bm \lambda}=[ \bm \lambda(1),\ldots,\bm \lambda(F)]^{\!\top\!}$ with $\|{\bm \lambda}\|_0=F$ is employed to `absorb' the norms of columns (so that the columns of $\boldsymbol{A}_n$ can be normalized w.r.t. a certain $\ell_q$ norm for $q\geq 1$).
We use  the shorthand notation $$\underline{{\bm X}} = [\![ {\bm \lambda}, {\bm A}_1,\ldots,{\bm A}_N ]\!]$$ to denote the CPD in \eqref{eq:individ}.
Note that $F\gg I_n$ for $n=1,\ldots,N$ could happen. In fact, the tensor rank can reach $O(I^{N-1})$ if $I_1=\ldots=I_N=I$ \cite{sidiropoulos2017tensor}.

\subsection{Joint PMF Recovery: A Tensor Perspective}
A key takeaway in \cite{kargas2017tensors} is that {\em any} joint PMF admits a naive Bayes (NB) model {\em representation}; i.e., any joint PMF can be generated from a latent variable model with just one hidden variable. It follows that the joint PMF of $\{Z_n\}_{n=1}^N$ can always be decomposed as
\begin{equation}\label{eq:pmf_latent_var}
\begin{aligned}
{\sf Pr}(i_1,i_2,\ldots,i_N) = \sum_{f=1}^F {\sf Pr}(f) \prod_{n=1}^N {\sf Pr}(i_n | f)
\end{aligned},
\end{equation}
where ${\sf Pr}(f) := {\sf Pr}(H = f) $ is the prior distribution of a latent variable $H$ and ${\sf Pr}(i_n | f) := {\sf Pr}( Z_n = z_n^{(i_n)} | H = f) $ are the conditional distributions.
Remarkably, such representation is {\em universal} if one allows the hidden variable to have a sufficiently rich finite alphabet:
\begin{theorem}\label{prop:F_upperbound}\cite{kargas2017tensors}
	The maximum $F$ needed in order to describe an arbitrary PMF using a naive Bayes model is bounded by
	$F\leq \underset{k}{\min}\prod_{\substack{n=1 \\ n \neq k}}^N I_n.$
\end{theorem}
To see this, one can represent any joint PMF as an $N$th-order tensor by letting $\underline{\bm X}(i_1,\ldots,i_N)={\sf Pr}(i_1,\ldots,i_N)$ and  ${\bm A}_n(i_n,f) = {\sf Pr}(i_n | f),~\boldsymbol{\lambda}(f) = {\sf Pr}(f).$
Since any nonnegative tensor admits a nonnegative polyadic decomposition with nonnegative latent factors, any joint PMF admits a naive Bayes representation with finite latent alphabet.

A relevant question is: what does `low tensor rank' mean for the joint PMF tensors? In \cite{kargas2017tensors}, the authors argue that reasonably dependent RVs lead to a joint PMF tensor with relatively small $F$---which is normally the case in statistical learning and inference. For example, classification is concerned with inferring $Z_N$ (label) from $Z_1,\ldots,Z_{N-1}$ (features), and the dependence amongs the RVs are leveraged for building up effective classifiers.

Another key observation in \cite{kargas2017tensors} is that the marginal distribution of any subset of RVs is an induced CPD model. From the law of total probability,
if one marginalize the joint PMF down to ${\sf Pr}(i_j,i_k,i_\ell)$ $\forall j,k,\ell \in [N],$ the expression is
$${\sf Pr}(i_j,i_k,i_\ell) = \sum_{f=1}^{F}   {\sf Pr}(f) {\sf Pr}(i_j | f) {\sf Pr}(i_k | f) {\sf Pr}(i_\ell | f).$$
Let $\underline{\bm X}_{jk\ell}(i_j,i_k,i_\ell)={\sf Pr}(i_j,i_k,i_\ell)$. Then, we have $\underline{{\bm X}}_{jk\ell} = [\![ \boldsymbol{\lambda},{\bm A}_j,{\bm A}_k,{\bm A}_\ell ]\!]$,
where $\{{\bm A}_n\}_{n=1}^N$ and ${\bm \lambda}$ are defined as before. Another important observation here is that
{\it the marginal PMFs $\underline{\bm X}_{jk\ell}$'s and the joint PMF $\underline{\bm X}$ share the same $\bm A_{j}$, $\boldsymbol{A}_{k}$, $\bm A_{\ell}$ and $\bm \lambda$}.
It is readily seen that if the $\underline{\bm X}_{jk\ell}$'s admit essentially unique CPD, then $\bm A_n$'s and $\bm \lambda$ can be identified from the marginals---and thus ${\sf Pr}(i_1,\ldots,i_N)$ can be recovered via identifying the latent factors.

The tensor-based approach in \cite{kargas2017tensors} has shown that recovering high-dimensional joint PMF from low-dimensional marginals is possible. However, a couple of major hurdles exist for applying it in practice. First, estimating three-dimensional marginals ${\sf Pr}(i_j,i_k,i_\ell)$ is still not easy, since one needs many co-realizations of three RVs. For sparse datasets, this is particularly hard. Second, tensor decomposition is a hard computation problem \cite{hillar2013most}. The coupled tensor decomposition (CTD) algorithm in \cite{kargas2017tensors} that involves many tensors is even more challenging in practice.

\subsection{Maximum Likelihood Estimation and EM} \label{sec:mle_em}

A recent work in \cite{yeredor2019maximum} offers an alternative for estimating ${\sf Pr}(f)$ and ${\sf Pr}(i_n|f)$.
Instead of working with third-order statistics, the model estimation problem is formulated as a maximum likelihood (ML) estimation problem in \cite{yeredor2019maximum}, and an EM algorithm is employed to handle the ML estimator (MLE).
The upshot of this approach is that the EM updates are much more scalable relative to the CTD updates in \cite{kargas2017tensors}.
However, there are two challenges of applying EM. First, the {\it identifiability} of the latent variables in \eqref{eq:pmf_latent_var} is unclear under the ML framework, since ML estimation does not imply identifiability; see discussions in the literature, e.g., \cite{huang2014putting}. Hence, the recoverability of the joint PMF may not be guaranteed.
Second, as the authors of \cite{yeredor2019maximum} have noticed, the EM algorithm sometimes does not converge to a good solution if randomly initialized.
This is perhaps due to the nonconvex nature of the ML estimation problem.

\section{Proposed Approach}
To advance the task of joint PMF recovery from marginal distributions, we propose a pairwise marginal-based approach.
Note that the pairwise marginal is given by
\begin{align*}
& {\sf Pr}(i_j,i_k)=\sum_{f=1}^F{\sf Pr}(f){\sf Pr}(i_j|f){\sf Pr}(i_k|f) \\&  \Longleftrightarrow
\boldsymbol{X}_{jk} = \boldsymbol{A}_j{\bm D}(\bm \lambda)\boldsymbol{A}_k^{\!\top\!},~\boldsymbol{X}_{jk}(i_j,i_k)={\sf Pr}(i_j,i_k),
\end{align*}
where $\bm D(\bm \lambda) = {\rm Diag}(\bm \lambda)$.

In practice, the pairwise marginals ${\bm X}_{jk}$'s are estimated from realizations of the joint PMF. Consider a set of realizations (data samples) of ${\sf Pr}(Z_1,\ldots,Z_N)$, denoted as $\{\bm d_{s} \in\mathbb{R}^{N} \}_{s=1}^S$.
Assuming that there are no missing observations, the following sample averaging estimator can be employed:
\begin{align*}
\widehat{\bm X}_{jk}(i_j,i_k) &= \frac{1}{S}\sum_{s=1}^S \mathbb{I}\left[\bm d_s(j)=z_j^{(i_j)},\bm d_s(k)=z_k^{(i_k)}\right],
\end{align*}
where $\bm d_s(n)$ denotes the realization of $Z_n$ in the $s$-th sample and $\mathbb{I}[E]=1$ if the event $E$ happens and $\mathbb{I}[E]=0$ otherwise.

{
Notably, with the same amount of data (i.e., with a fixed $S$), the second-order statistics can be estimated to a much higher accuracy, compared to the third-order ones \cite{han2015minimax}.
This is particularly articulated when the dataset is sparse.
To see this, assume that every RV has an alphabet size $I$ and is observed with probability $p$ (which means that every entry of $\bm d_s$ is observed with probability $p$---and $\bm d_s$ is very sparse when $p$ is small). Then, on average $Sp^2$ and $Sp^3$ co-realizations of every pair and triple of RVs, respectively, can be observed in the samples.
Under such cases, Theorem~1 in \cite{han2015minimax} shows that the average estimation error (over all possible joint distributions) for $\underline{\bm X}_{jk\ell}$ in terms of $\ell_1$ distance is $O(\sqrt{\nicefrac{I^3}{Sp^3}})$,
while that for $\boldsymbol{X}_{jk}$ is $O( \sqrt{\nicefrac{I^2}{Sp^2}})$. One can see that the difference is particularly significant when $p$ is small---i.e., when the dataset is considerably sparse.}

\subsection{Challenges}
Despite the sample complexity appeal of using pairwise marginals, the {\it recoverability} of the joint PMF under such settings is nontrivial to establish. The main analytical tool for proving recoverability in \cite{kargas2017tensors} is the essential uniqueness of tensor decomposition---i.e., the latent factors of low-rank tensors are identifiable up to resolvable ambiguities under mild conditions.
However, pairwise distributions such as $\bm X_{jk}= \boldsymbol{A}_j{\bm D}(\bm \lambda)\boldsymbol{A}_k^{\!\top\!}$ are matrices, and low-rank matrix decomposition is in general \textit{nonunique}.
A natural thought to handle the identifiability problem would be employing certain NMF tools \cite{fu2018nonnegative,gillis2014and}, since the latent factors are all nonnegative, per their physical interpretations.
However, the identifiability of NMF models holds only if $F\leq \min\{I_j,I_k\}$ (and preferably $F\ll\min\{I_j,I_k\}$).
The pairs $\bm X_{jk}= \boldsymbol{A}_j{\bm D}(\bm \lambda)\boldsymbol{A}_k^{\!\top\!}\in\mathbb{R}^{I_j\times I_k}$ inherit the inner dimension $F$ (i.e., the column dimension of $\bm A_j$) from the joint PMF of all the variables, which is the tensor rank of an $N$th-order tensor. As mentioned, the tensor rank $F$ could be much larger than the $I_j$'s.
Hence, one may not directly use the available NMF uniqueness results on individual $\boldsymbol{X}_{jk}$'s to argue for joint PMF recoverability.
Nonetheless, as we will show shortly, NMF identifiability can be applied onto carefully constructed {\it coupled} NMF models to establish recoverability.

\subsection{Preliminaries: NMF and Model Identifiablility}
To see how we approach these challenges, let us first briefly introduce some pertinent preliminaries of NMF.
Consider a nonnegative data matrix $\boldsymbol{X}\in\mathbb{R}^{L\times K}$. Assume that the matrix is generated by $\boldsymbol{X}=\boldsymbol{W}\H^{\!\top\!}$, where $\boldsymbol{W}\in\mathbb{R}^{L\times F}$ and $\H \in \mathbb{R}^{K\times F}$. Also assume that $\boldsymbol{W}\geq \bm 0$ and $\H\geq \bm 0$ and $F\leq \min\{L,K\}$. NMF aims to factor $\boldsymbol{X}$ for identifying $\boldsymbol{W}$ and $\H$ up to certain trivial ambiguities.

A number of conditions on $\boldsymbol{W}$ and/or $\H$ allow us to establish identifiability of the latent factors. The first one is the so-called separability condition:
\begin{definition} (Separability) \label{def:sep} \cite{donoho2003does}
	If $\bm H\geq \bm 0$, and  $\bm \varLambda =\{ l_1,\ldots,l_F \}$ such that $\bm H(\bm \varLambda,:)=\bm \Sigma$ holds, where $\bm \Sigma ={\rm Diag}(\alpha_1,\ldots,\alpha_F)$ and $\alpha_f>0$, then, $\H$ satisfies the {\it separability condition}. When $\bm \varLambda =\{ l_1,\ldots,l_F \}$ satisfies $\|\H(l_f,:)-\bm e_f\|_2\leq \varepsilon$ for $f=1,\ldots,F$, $\bm H$ is called $\varepsilon$-separable.
\end{definition}
If one of $\boldsymbol{W}$ and $\H$ satisfies the separability condition and the other has full column rank, there exist polynomial-time algorithms that can provably identify $\boldsymbol{W}$ and $\H$ up to scaling and permutation ambiguities.
Among these algorithms, a number of Gram--Schmidt-like lightweight algorithms exist; see, e.g., \cite{Gillis2012,fu2014self,arora2012practical,VMAX,MC01}.
A more relaxed condition is as follows:
\begin{definition}\label{def:suff}  (Sufficiently Scattered) \cite{huang2014non,fu2015blind,fu2016robust}
	Assume that $\bm H\geq \bm 0$ and
	${\cal C}\subseteq {\rm cone}\{\bm H^{\!\top\!}\}$
	where 	$\mathcal{C} = \{\bm{x} \in\mathbb{R}^F~|~ \bm{x}^{\!\top\!}\mathbf{1} \geq \sqrt{F-1}\|\bm{x}\|_2\}$.
	In addition, assume that ${\sf cone}\{\H^{\!\top\!}\} \not\subseteq {\sf cone}\{\bm Q\} $ for any orthonormal $\boldsymbol{Q}\in\mathbb{R}^{K\times K}$ except for the permutation matrices.
	Then, $\bm H$ is called {\it sufficiently scattered}.
\end{definition}
The sufficiently scattered condition subsumes the separability condition as a special case. An important result is that if $\boldsymbol{W}$ and $\H$ are both sufficiently scattered, then the model $\boldsymbol{X}=\boldsymbol{W}\H^{\!\top\!}$ is unique up to scaling and permutation ambiguities \cite{huang2014non}. The two conditions are illustrated in Fig.~\ref{fig:cone}.
More detailed discussions and illustrations can be seen in \cite{fu2018nonnegative,fu2018identifiability}.

In principle, the `taller' $\bm H$ is, the chance that $\bm H$ (approximately) attains the separability or sufficiently scattered condition is higher, if $\H$ is generated following a certain continuous distribution. This intuition was formalized in \cite{ibrahim2019crowdsourcing}, under the assumption that $\H(l,:)$'s are drawn from the probability simplex uniformly at random and then positively scaled. The work in \cite{ibrahim2019crowdsourcing} also shows that the sufficiently scattered condition is easier to satisfy relative to separability, under the same generative model; see Lemma~\ref{lem:Lm} in Appendix~\ref{app:lemmata} in the supplementary materials.

\begin{figure}[t!]
	\centering
	\includegraphics[width= 0.30\linewidth]{figures/simplex3}
	\includegraphics[width= 0.25\linewidth]{figures/simplex6}
	\put(-35,180){$\mathcal{C}$}
	\caption{Left: an $\bm H$ that satisfies the separability condition. Right: an $\H$ that satisfies the sufficiently scattered condition. The triangle corresponds to probability simplex $\Delta_3$; and the dots correspond to rows of $\H$ rescaled to have unit $\ell_1$-norm. The inner circle on $\Delta_3$ corresponds to the intersection of the second-order cone ${\cal C}$ and $\Delta_3$; and the shaded region is the intersection of $\cone{\bm H^{\!\top\!}}$ and $\Delta_3$.}
	\label{fig:cone}
	\vspace{-.35cm}
\end{figure}

\subsection{Recoverability Analysis}
Our goal is to identify $\boldsymbol{A}_n$ and $\bm \lambda$ from the available pairwise marginals $\boldsymbol{X}_{jk}=\boldsymbol{A}_j\bm D(\bm \lambda)\boldsymbol{A}_k^{\!\top\!}$'s, so that the joint PMF can be reconstructed following \eqref{eq:pmf_latent_var}.
Note that applying existing NMF techniques onto $\boldsymbol{X}_{jk}$ is unlikely to work, since when $F> \min \{I_j,I_k\}$, `fat' $\boldsymbol{A}_j\in\mathbb{R}^{I_j\times F}$'s can never be separable or sufficiently scattered. To circumvent this, our idea is as follows.
Consider a splitting of the indices of the $N$ variables, i.e., ${\cal S}_1=\{\ell_1,\ldots,\ell_M\}$ and ${\cal S}_2=\{\ell_{M+1},\ldots,\ell_{N}\}$ such that
${\cal S}_1\cup{\cal S}_2=\{1,\ldots,N\}$ and ${\cal S}_1\cap {\cal S}_2=\emptyset.$
Then, we construct the following matrix:
\begin{equation}\label{eq:constructedX}
\begin{aligned}
\widetilde{\bm X}=& \begin{bmatrix} \bm X_{\ell_1\ell_{M+1}}& \ldots & {\bm X}_{\ell_1\ell_N} \\ \vdots & \vdots &\vdots \\ {\bm X}_{\ell_M\ell_{M+1}} & \ldots & {\bm X}_{\ell_M\ell_N} \end{bmatrix} \\
= &\underbrace{\begin{bmatrix}
	\bm A_{\ell_1}\\ \vdots\\ \bm A_{\ell_M}
	\end{bmatrix}}_{\bm W}\underbrace{{\bm D}({\bm \lambda })[{\bm A}_{\ell_{M+1}}^{\!\top\!},\ldots,{\bm A}_{\ell_{N}}^{\!\top\!}]}_{{\bm H}^{\!\top\!}}.
\end{aligned}
\end{equation}
Note that in this `virtual NMF' model, $\boldsymbol{W}$ and $\H$ matrices have sizes of $MI\times F$ and $(N-M)I\times F$, respectively, if $I_1=\ldots=I_N=I$.
The idea is to construct $\widetilde{\bm X}$ such that $F\leq\min\{MI, (N-M)I \}$ so that $\boldsymbol{W}$ and $\H$ may satisfy the separability or sufficiently scattered condition---and thus $\boldsymbol{A}_1,\ldots,\boldsymbol{A}_N$
can be identified via factoring $\widetilde{\boldsymbol{X}}$.
This is a straightforward application of the NMF identifiability results introduced in the previous subsection.

The challenge lies in finding a suitable splitting of ${\cal S}_1,{\cal S}_2$ such that $\boldsymbol{W}$ and $\H$ are sufficiently scattered and is highly nontrivial.
Note that even if $\boldsymbol{W}$ and $\H$ are given, checking if these matrices are sufficiently scattered is hard \cite{huang2014non}.
It is also impractical to try every possible index splitting.
To address this challenge, we consider the following factorization problem:
\begin{subequations}\label{eq:coupled_mat}
	\begin{align}
	{\rm find}&~\boldsymbol{A}_n,~n=1,\ldots,N,~\bm \lambda\\
	{\rm subject~to}&~\boldsymbol{X}_{jk}=\boldsymbol{A}_j\bm D(\bm \lambda)\boldsymbol{A}_k^{\!\top\!},~\forall  j,k \in{\bm \Omega} \label{eq:fic_feas}\\
	&~\bm 1^{\!\top\!}\bm A_j=\bm 1^{\!\top\!},~\bm A_j\geq \bm 0\\
	&~\bm 1^{\!\top\!}\bm \lambda =1,~\bm \lambda\geq{\bm 0},
	\end{align}
\end{subequations}
where $\bm \Omega$ contains the index set of $(j,k)$'s such that $j<k$ and the joint PMF ${\sf Pr}(i_j,i_k)$ is accessible.
Note that the above can be understood as a latent factor-coupled NMF problem.
Regarding the identifiability of the conditional PMFs and the latent prior, we have the following result:
\begin{theorem}[Recoverability]\label{prop:prop2}
	Assume that that ${\sf Pr}(i_j,i_k)$'s for $j,k\in{\bm \Omega}$ are available and that ${\sf Pr}(f)\neq 0$ for $f=1,\ldots,F$.
	Suppose that there exists ${\cal S}_1=\{\ell_1,\ldots,\ell_M\}$ and ${\cal S}_2=\{\ell_{M+1},\ldots,\ell_Q\}$ such that $Q\leq N$
	and $ {\cal S}_1\cup{\cal S}_2\subseteq \{1,\ldots,N\},~{\cal S}_1\cap {\cal S}_2=\emptyset. $
	Also assume the following conditions hold:\\
	\textbf{(i)} the matrices $[\boldsymbol{A}_{\ell_1}^{\!\top\!},\ldots,\boldsymbol{A}_{\ell_M}^{\!\top\!}]^{\!\top\!}$ and $[\boldsymbol{A}_{\ell_{M+1}}^{\!\top\!},\ldots,\boldsymbol{A}_{\ell_{Q}}^{\!\top\!}]^{\!\top\!}$ are \textit{sufficiently scattered};\\
	\textbf{(ii) }all pairwise marginal distributions ${\sf Pr}(i_j,i_k)$'s for $j\in{\cal S}_1$ and $k\in{\cal S}_2$ are available; \\
	\textbf{(iii)} every $T$-concatenation of $\boldsymbol{A}_n$'s, i.e., $[\boldsymbol{A}_{n_1}^{\!\top\!},\ldots,\boldsymbol{A}_{n_T}^{\!\top\!}]^{\!\top\!}$, is a full column rank matrix, if $I_{n_1}+\ldots+I_{n_T}\geq F$;\\
	\textbf{(iv) }for every $j\notin{\cal S}_1\cup{\cal S}_2$ there exists a set of $r_t\in{\cal S}_1\cup{\cal S}_2$ for $t=1,\ldots,T$ such that ${\sf Pr}(i_j,i_{r_t})$ or ${\sf Pr}(i_{r_t},i_j)$ are available.

	Then, solving Problem~\eqref{eq:coupled_mat} recovers ${\sf Pr}(i_j|f)$ and ${\sf Pr}(f)$ for $j=1,\ldots,N, ~f=1,\dots, F$, thereby the joint PMF ${\sf Pr}(i_1,\ldots,i_N)$.
\end{theorem}

The proof is relegated to Appendix~\ref{app:prop1}, which is reminiscent of identifiability of a similar coupled NMF problem that arises in crowdsourced data labeling \cite{ibrahim2019crowdsourcing}, with proper modifications to handle the cases where $F>I_n$ (since the crowdsourcing model always has $F=I_n$, which is a simpler case).
In particular, conditions (iii) and (iv) are employed to accommodate the more challenging cases.
The criterion spares one the effort for first finding ${\cal S}_1$ and ${\cal S}_2$ and then constructing the matrix $\widetilde{\boldsymbol{X}}$. Instead, under Theorem~\ref{prop:prop2}, the goal is simply finding an approximate solution of \eqref{eq:coupled_mat}---which is much more approachable in practice.

\begin{remark}
{
Under our pairwise marginal based approach, a key enabler for recovering the joint PMF is the {\it identifiability} of $\boldsymbol{A}_n$'s up to an {\it identical} column permutation ambiguity.
At first glance, the identifiability of the latent factors seems to be irrelevant since what we truly care is the ambient joint PMF tensor $\underline{\bm X}\in\mathbb{R}^{I_1\times\ldots\times I_N}$.
However, since the tensor is heavily marginalized (i.e.,     one only observes pairwise marginals), the identifiability of $\bm A_n$'s is leveraged on to assemble the joint PMF following the model in \eqref{eq:pmf_latent_var} under our framework.
}
\end{remark}

\begin{comment}

Theorem~\ref{prop:prop2} does not impose any restrictions on $F$, and thus can be generally applicable.
Nevertheless, to establish recoverability, one has to ensure that two `virtual matrices'  (i.e., $\boldsymbol{W}$ and $\H$) comprising different conditional PMF matrices to be sufficiently scattered---which might be hard to guarantee in practice.

Nevertheless, we find that a stronger result can be derived if $F$ is below a certain threshold.
To simplify the statement, let us assume $I_1=\ldots=I_N=I$. We have:
\begin{theorem}[Enhanced Recoverability]
	\label{prop:completion}
	Assume that ${\sf Pr}(f)\neq 0$ for $f=1,\ldots,F$, and that ${\sf Pr}(i_j,i_k)$'s for all $j,k$ are available and ${\sf Pr}(i_k,i_j)={\sf Pr}(i_j,i_k)$\footnote{This means that the RVs are pairwise exchangeable, which normally holds if the realizations are i.i.d.}.
	If $\boldsymbol{Z}=[\boldsymbol{A}_1^{\!\top\!},\ldots,\boldsymbol{A}_N^{\!\top\!}]^{\!\top\!}\in\mathbb{R}^{NI\times F}$ is separable or sufficiently scattered and $F\leq (N-1)I - 1$, then, solving the problem in \eqref{eq:coupled_mat} recovers ${\sf Pr}(i_j|f)$ and ${\sf Pr}(f)$ for $j=1,\ldots,N, ~f=1,\dots, F$, thereby the joint PMF ${\sf Pr}(i_1,\ldots,i_N)$.
\end{theorem}
The proof is relegated to Appendix~\ref{app:prop2}. Theorem~\ref{prop:completion} is stronger in the following sense.
In Theorem~\ref{prop:prop2}, the recoverability of the joint PMF depends on if $\boldsymbol{W}=[\boldsymbol{A}_{\ell_1}^{\!\top\!},\ldots,\boldsymbol{A}_{\ell_M}^{\!\top\!}]^{\!\top\!}$ and $\H=[\boldsymbol{A}_{\ell_{M+1}}^{\!\top\!},\ldots,\boldsymbol{A}_{\ell_N}^{\!\top\!}]^{\!\top\!}$ are sufficiently scattered. This requires that there exists at least one splitting strategy for $\boldsymbol{A}_1,\ldots,\boldsymbol{A}_N$ under which two sufficiently scattered/separable matrices can be constructed. However, under Theorem~\ref{prop:completion}, one only needs that all the conditional PMFs can together construct a sufficiently scattered $\boldsymbol{Z}$. Note that $\boldsymbol{Z}$ could thus have many more rows than those of $\boldsymbol{W}$ and $\H$.
In fact, it was shown in \cite{ibrahim2019crowdsourcing} that if the rows of $\boldsymbol{Z}$ are randomly drawn from the nonnegative orthant, then the chances for $\boldsymbol{Z}$ being separable or sufficiently scattered are much higher when the number of rows of $\boldsymbol{Z}$ grows.
Therefore, Theorem~\ref{prop:completion} is a substantially improved recoverability result relative to Theorem~\ref{prop:prop2}.

\end{comment}

\section{Algorithms and Performance Analyses}
In this section, we develop algorithms under the coupled NMF framework and analyze performance under realistic conditions, e.g., cases with finite samples and missing values.
\subsection{BCD for Problem~\eqref{eq:coupled_mat}} \label{sec:alg1}
To handle Problem \eqref{eq:coupled_mat}, we employ the standard procedure based on \textit{block coordinate descent} (BCD) for handling coupled tensor/matrix decomposition algorithms in the literature \cite{kargas2017tensors,kargas2019learning,ibrahim2019crowdsourcing}.
To be specific, we cyclically minimize the constrained optimization problem w.r.t. $\boldsymbol{A}_k$, when fixing $\boldsymbol{A}_j$ for all $j\neq k$ and $\bm \lambda$. Consequently, the subproblem is as follows:
\begin{subequations}\label{eq:sub}
	\begin{align}
	\minimize_{\boldsymbol{A}_k}&~\sum_{j\in{\bm \Omega}_k}{\sf dist}\left( \boldsymbol{X}_{jk},\boldsymbol{A}_j\bm D(\bm \lambda)\boldsymbol{A}_k^{\!\top\!}\right)\\
	{\rm subject~to}&~\bm 1^{\!\top\!}\bm A_k=\bm 1^{\!\top\!},~\bm A_k\geq \bm 0,
	\end{align}
\end{subequations}
where ${\bm \Omega}_k$ is the index set of $j$ such that ${\sf Pr}(i_j,i_k)$ is available and ${\rm dist}(\cdot,\cdot)$ is a certain distance measure. Under a number of distance measures, e.g., the Euclidean distance or the Kullback--Leibler (KL) distance, the subproblems are convex and thus can be solved via off-the-shelf algorithms.
In this work, we employ the KL divergence as the distance measure since it is natural for measuring distance between PMFs.
Consequently, the subproblems are
\begin{subequations}\label{eq:subkl}
	\begin{align}
	\boldsymbol{A}_k &\leftarrow\arg\min_{ \begin{subarray}{c} \bm 1^{\!\top\!}\bm A_k=\bm 1^{\!\top\!}\\\bm A_k\geq \bm 0\end{subarray}}~\sum_{j\in{\bm \Omega}_k}\mathbb{D}_{\sf KL}\left( \boldsymbol{X}_{jk},\boldsymbol{A}_j\bm D(\bm \lambda)\boldsymbol{A}_k^{\!\top\!}\right), \label{eq:Aupdate}\\
	\bm \lambda &\leftarrow\arg\min_{ \begin{subarray}{c} \bm 1^{\!\top\!}\bm \lambda=1\\ \bm \lambda\geq \bm 0\end{subarray}}~\sum_{j,k\in{\bm \Omega}}\mathbb{D}_{\sf KL}\left( \boldsymbol{X}_{jk},\boldsymbol{A}_j\bm D(\bm \lambda)\boldsymbol{A}_k^{\!\top\!}\right). \label{eq:lambda_update}
	\end{align}
\end{subequations}
We employ the {\it mirror descent} algorithm to solve the two subproblems. Mirror descent admits closed-form updates for handling KL-divergence based linear model fitting under simplex constraints (which is also known as exponential gradient in the literature); see, e.g., \cite{arora2012practical}. The algorithm {CNMF-OPT} is described in \eqref{algo:CNMF-OPT}.
\begin{algorithm}[h!]
	\SetKwInOut{Input}{input}
	\SetKwInOut{Output}{output}
	\SetKwRepeat{Repeat}{repeat}{until}

	\Input{data samples $\{\bm d_s\}_{s=1}^S$}

	 estimate second order statistics $\widehat{\boldsymbol{X}}_{jk }$;

	 get initial estimates $\{\widehat{\boldsymbol{A}}_n\}_{n=1}^N$ and $\widehat{\bm \lambda}$ (e.g., using {CNMF-SPA});

		\Repeat{some stopping criterion is reached}{

		\For{$n=1$ {\bf to} $N$}{

		  update $\widehat{\boldsymbol{A}}_n \leftarrow \eqref{eq:Aupdate}$ using mirror descent;
		}

		 update $\widehat{\bm \lambda} \leftarrow\eqref{eq:lambda_update}$ using mirror descent;
		}
		\Output{estimates $\{\widehat{\boldsymbol{A}}_n\}_{n=1}^N$, $\widehat{\bm \lambda}$}
			\caption{{CNMF-OPT}}\label{algo:CNMF-OPT}
\end{algorithm}

Notably, the work in \cite{kargas2017tensors} also has some experiments (without recoverability analysis) using pairwise marginals to recover the joint PMF using a similar formulation (using Euclidean distance).
Our rationale for reaching the formulation in \eqref{eq:coupled_mat} is very different, which is motivated by the NMF theory.
Also, the result using pairwise marginals in \cite{kargas2017tensors} seems not promising.
In this work, our analysis in Theorems~\ref{prop:prop2} has suggested that the reason why pairwise marginals do not work in \cite{kargas2017tensors} is likely due to optimization pitfalls, rather than lack of recoverability.
In the next subsection, we show that with a carefully designed initialization scheme, accurately recovering joint PMFs from pairwise marginals is viable.

\subsection{Key Step: Gram--Schmidt-like Initialization}	\label{sec:alg2}
Directly applying the BCD algorithm to the joint PMF recovery problem has a number of challenges. First, the CNMF problem is nonconvex, and thus global optimality is not ensured.
Second, the computational complexity grows quickly with $N$ and $I$, giving rise to potential scalability issues.
In this subsection, we offer a fast algorithm that can provably identify ${\sf Pr}(i_n|f)$ and ${\sf Pr}(f)$ up to a certain accuracy---under more stringent conditions relative to those in Theorem~\ref{prop:prop2}. The algorithm is reminiscent of a Gram--Schmidt-like algorithm for NMF with careful construction of a `virtual NMF model'.

Let us assume that there is a splitting ${\cal S}_1=\{\ell_1,\ldots,\ell_M\}$ and ${\cal S}_2=\{\ell_{M+1},\ldots,\ell_{N}\}$ such that the $\widetilde{\boldsymbol{X}}=\boldsymbol{W}\H^{\!\top\!}$ in \eqref{eq:constructedX} admits a full rank $\boldsymbol{W}$ and a separable $\H$.
Under the separability condition, we have ${\bm H}(\bm \varLambda,:) =\bm \Sigma= {\rm Diag}(\alpha_1,\ldots,\alpha_F)$ and ${\bm W}\bm \Sigma = \widetilde{\boldsymbol{X}}(\bm \varLambda,:).$ Hence, the coupled NMF task boils down to identifying the index set $\bm \varLambda$.  The so-called \textit{successive projection algorithm} (SPA) from the NMF literature \cite{Gillis2012,MC01,fu2014self,arora2012practical} can be employed for this purpose. A remark is that SPA is a Gram--Schmidt-like algorithm, which only consists of norm comparison and orthogonal projection.

Once ${\bm W}$ is identified, one can recover $\bm A_{\ell_n}\in \mathbb{R}^{I_{\ell_n} \times F}$ for $\ell_n\in{\cal S}_1$ up to identical column permutations, by extracting the corresponding rows of $\bm W$ (cf. line 6 in Algorithm~\ref{algo:CNMF-SPA}).
Unlike general NMF models, since every column of $\boldsymbol{A}_n$ is a conditional PMF, there is no scaling ambiguity.
The $\H$ matrix can be estimated using (constrained) least squares, and $\boldsymbol{A}_{\ell_n}$ for $\ell_n\in {\cal S}_2$ can then be extracted in a similar way.
Denote \eqref{eq:constructedX} as $\widetilde{\bm X} = {\bm W} \bm D(\bm \lambda) \widetilde{\bm H}^{\top}$, where $\widetilde{\bm H} = \begin{bmatrix} \bm A_{\ell_{M+1}}^{\top} \ldots \bm A_{\ell_N}^{\top} \end{bmatrix}^{\!\top\!}$.
Then, the PMF of the latent variable can be estimated via
$
{ \bm \lambda} = ( \widetilde{\bm H} \odot \boldsymbol{W})^{\dagger}\text{vec}(\widetilde{\bm X}),$
where we have used the fact that the Khatri-Rao product $\widetilde{\bm H} \odot \boldsymbol{W}$ has full column rank since both ${\boldsymbol{W}}$ and $\widetilde\H$ have full column rank.

A remark is that the SPA-estimated $\widehat{\bm A}_{\ell_n}$ is at best the column-permuted version of $\boldsymbol{A}_{\ell_n}$ \cite{Gillis2012}.
However, since the permutation ambiguity across all the $\bm A_n$'s and the $\bm \lambda$ are identical per the above procedure, the existence of column permutations does not affect the ``assembling'' of ${\sf Pr}(i_n|f)$ and ${\sf Pr}(f)$ to recover ${\sf Pr}(i_1,\ldots,i_N)$.
We refer to this procedure as {\it coupled NMF via SPA} ({CNMF-SPA}); see Algorithm \ref{algo:CNMF-SPA}.

From Algorithm~\ref{algo:CNMF-SPA}, one can see that the procedure is rather simple since SPA is a Gram--Schmidt-like procedure and the other steps are either convex quadratic programming or least squares.
Hence, {CNMF-SPA} can be employed for initializing computationally more intensive algorithms, e.g., {CNMF-OPT}.

	\begin{algorithm}[h]
	\SetKwInOut{Input}{input}
	\SetKwInOut{Output}{output}
	\SetKwRepeat{Repeat}{for}{end}

	\Input{data samples $\{\bm d_s\}_{s=1}^S$ and $M$}
	 estimate second order statistics $\widehat{\boldsymbol{X}}_{jk }$;

    split $\{1,\ldots,N\}$ into ${\cal S}_1=\{1,\ldots,M\}$ and ${\cal S}_2=\{M+1,\ldots,N\}$;

	 Construct $\widetilde{\bm X}$

	 Estimate $\widehat{\bm W}$ using the {\bf SPA algorithm} \cite{Gillis2012} to select $\bm \varLambda$;

	\For{$n=1$ {\bf to} $M$}{
	 $\widehat{\bm A}_n \leftarrow \widehat{\bm W} ((n-1)I_n+1:nI_n,:)$;

	 normalize columns of $\widehat{\bm A}_n$ with respect to $\ell_1$ norm;}

 $\widehat{\bm H} \leftarrow \underset{{\bm H} \ge 0}{\text{arg~min}}~\|\widetilde{\bm X}- \widehat{\bm W} {\bm H}^{{\!\top\!}}\|_{\rm F}^2$

	\For{$n=M+1$ {\bf to} $N$}{
	 $\widehat{\bm A}_n \leftarrow \widehat{\bm H} ((n-1)I_n+1:nI_n,:)$;

		 normalize columns of $\widehat{\bm A}_n$ with respect to $\ell_1$ norm;}

	 $\widetilde{\bm W}^{{\!\top\!}} \leftarrow \begin{bmatrix} \widehat{\bm A}_1^{{\!\top\!}}, \hdots, \widehat{\bm A}_M^{{\!\top\!}} \end{bmatrix}$ ;

	$\widetilde{\bm H}^{{\!\top\!}} \leftarrow  \begin{bmatrix} \widehat{\bm A}_{M+1}^{\top} \hdots \widehat{\bm A}_{N}^{\top} \end{bmatrix}^{\top} $;

 $\widehat{\bm \lambda} \leftarrow ( \widetilde{\bm H} \odot \widetilde{\bm W})^{\dagger}\text{vec}(\widetilde{\bm X})$;

	 \Output{estimates $\{\widehat{\boldsymbol{A}}_n\}_{n=1}^N$, $\widehat{\bm \lambda}$.}
	 \caption{{CNMF-SPA}}\label{algo:CNMF-SPA}

\end{algorithm}
\vspace{-.25cm}

\subsection{Performance Analysis of CNMF-SPA} \label{sec:analysisSPA}
In principle, to make {CNMF-SPA} work, one needs to identify a splitting scheme ${\cal S}_1$ and ${\cal S}_2$ such that $\bm H$ in \eqref{eq:constructedX} satisfies the separability condition.
It is impossible to know such a splitting {\it a priori}.
Testing all combinations of ${\cal S}_1$ and ${\cal S}_2$ gives rise to a hard combinatorial problem.
In addition, the pairwise marginals $\boldsymbol{X}_{jk}$ are not perfectly estimated in practice.
In this subsection, we offer performance analysis of {CNMF-SPA} by taking the above aspects into consideration.

From our extensive experiments, we observed that using the `naive' construction ${\cal S}_1=\{1,\ldots,M\}$ and ${\cal S}_2=\{M+1,\ldots,N\}$ seems to work reasonably well for various tasks.
To establish theoretical understanding to such effectiveness, recall that we use $S$ and $p$ to denote the number of available realizations of ${\sf Pr}(Z_1,\ldots,Z_N)$ and the probability of observing each entry of data sample $\bm d_s$, respectively.
For simplicity, we assume that $I_n=I$ for all $n$.
We impose the following generative model for $\boldsymbol{A}_m$:

\begin{Assumption}\label{ass:A}
	Assume that every row of $\boldsymbol{A}_n\in\mathbb{R}^{I_m\times F}$ for all $n$ is generated from the $(F-1)$-probability simplex uniformly at random and then positively scaled by a scalar, so that $\bm 1^{\!\top\!}\bm A_n=\bm 1^{\!\top\!}$ is respected\footnote{Such scaling exists under mild conditions; see \cite[Proposition 1]{yang2019learning}.}.
\end{Assumption}
\color{black}
Under the above settings, we show that the following holds:
\begin{theorem} \label{thm:spabound}
	Assume that $\|\widehat{\bm X}_{ij}(:,q)\|_1 \ge \eta>0$ for any $q,i,j$.
	Also, assume that
	$M \ge  F/I $, $p \geq \left( \frac{8}{S}\log(4/\delta) \right)^{1/2}$,
	\begin{align*}
	S & = \Omega\left(\frac{M^2I{\rm log}(1/\delta)}{\sigma^2_{\rm max}({\bm{W}})\eta^2\varepsilon^2p^2}\right),\\
	N &= M+ \Omega\left(\frac{\varepsilon^{-2(F-1)}}{IF}{\rm log}\left(\frac{F}{\delta}\right)\right),
	\end{align*}
	where $\varepsilon = \frac{M\text{\rm min}\left(\frac{1}{2\sqrt{F-1}},\frac{1}{4}\right)}{2\kappa(\bm{W})(1+80\kappa^2(\bm{W}))}$.
	Then,   under the defined $S,p$ and Assumption~\ref{ass:A}, {CNMF-SPA} outputs $\widehat{\boldsymbol{A}}_m$'s such that
	\begin{align} \label{eq:thmspa}
	\min_{\bm \Pi:~\text{permuation}}\|\widehat\boldsymbol{A}_m\bm \Pi - \boldsymbol{A}_m\|_2 = O\left(\kappa^2(\bm W)\sqrt{F} \zeta\right)
	\end{align}
	for $m\in{\cal S}_1$ with a probability greater than or equal to $1-\delta$, where $\zeta = \max(\sigma_{\max}(\bm W)\varepsilon,\nicefrac{M\sqrt{I\log(1/\delta)}}{\eta p \sqrt{S}})$.
\end{theorem}
The proof is relegated to Appendix~\ref{app:spa}.
Theorem~\ref{thm:spabound} asserts that using {CNMF-SPA} on the constructed $\widetilde{\boldsymbol{X}}$ in \eqref{eq:constructedX} approximately recovers the joint PMF---if the $N,I,S$ are large enough.
Note that if $\boldsymbol{A}_m$ for all $m \in{\cal S}_1$ can be accurately estimated, the estimation accuracy of
$\bm A_n$ for all $n \in{\cal S}_2$ and $\bm \lambda$ can also be guaranteed and quantified, following the standard sensitivity analyses of least squares; see, e.g., \cite{wedin1973perturbation}.

\begin{remark} \label{rem:theorem_spa}
Theorem~\ref{thm:spabound} is not entirely surprising.
The insight behind the proof is to recast the constructed $\widetilde{\boldsymbol{X}}$ as an equivalent noisy NMF model with the separability condition holding exactly on the right latent factor.
This way, the noise robustness analysis for SPA \cite{Gillis2012} can be applied.
The key is to quantify the noise bound in this equivalent model. Such `virtual noise' is contributed by the finite sample-induced error for estimating the pairwise marginals and the violation to the separability condition.
In particular,
classical concentration theorems are leveraged to quantify the former, and Assumption~\ref{ass:A} for the latter.
Note that Assumption~\ref{ass:A} is more of a working assumption, rather than an exact model for capturing reality---it helps formalize the fact that $\H$ in \eqref{eq:constructedX} more likely attains the $\varepsilon$-separability condition when $|{\cal S}_2|$ grows \cite{ibrahim2019crowdsourcing}. In principle, if $\boldsymbol{A}_m$ is drawn from any joint absolutely continuous distribution, a similar conclusion can be reached---but using the uniform distribution simplifies the analysis.
\end{remark}

\color{black}

\subsection{More Discussions: EM Meets CNMF}
As mentioned, a recent work in \cite{yeredor2019maximum} proposed an EM algorithm for solving the maximum likelihood estimation problem under \eqref{eq:pmf_latent_var}.
To see the idea, let $f_s \in \{1,\dots,F\}$ be the realization of the `latent variable' $H$ in the $s$th realization of ${\sf Pr}(Z_1,\ldots,Z_N)$. The joint log-likelihood of the observed data $\{\bm d_s\}_{s=1}^S$ and the corresponding latent variable $\{f_s\}_{s=1}^S$ as a function of the latent factor $\bm \theta =\begin{bmatrix}{\rm vec}(\bm A_1)^{\!\top\!},\dots,{\rm vec}(\bm A_N)^{\!\top\!}, \bm \lambda^{{\!\top\!}}\end{bmatrix}^{\!\top\!}$ can be expressed as
\begin{align}\label{eq:loglike}
{\cal L}(\{\bm d_s,f_s\}_{s=1}^S;\bm \theta)
&=\log\left(\prod_{s=1}^S \bm \lambda(f_s) \prod_{n=1}^N\prod_{i=1}^{I_n} \bm A_n(i,f_s)^{\mathbb{I}(\bm d_s(n)=z_{n}^{(i)})}\right).
\end{align}
The EM algorithm estimates $\bm \theta$ via attempting to maximize the log-likelihood. The resulting updates are simple and economical. However, the optimality of the algorithm is unclear. It was also observed in \cite{yeredor2019maximum} that the algorithm sometimes converges to some undesired solutions if not carefully initialized.

Naturally, our idea is using {CNMF-SPA} to initialize the EM algorithm in \cite{yeredor2019maximum}. As one will see, this combination oftentimes offers competitive performance
This combination also features an appealing runtime performance, since both algorithms consist of lightweight updates. We refer to this procedure as CNMF-SPA-initialized EM ({CNMF-SPA-EM}).

In this work, we also present performance characterizations for the EM algorithm.
To proceed, we make the following assumption:

\begin{Assumption}\label{ass:D}
	Define
	$$\overline{D}_1 = \min_{f \neq f'}\frac{1}{N}\sum_{n=1}^N p\mathbb{D}_{\rm KL}(\bm A_n(:,f),\bm A_n(:,f')),\quad
	\overline{D}_2 = \frac{2}{N}\min_{f\neq f'} \log(\bm \lambda(f)/\bm \lambda(f'))$$
	and $\overline{D}=\nicefrac{(\overline{D}_1 + \overline{D}_2)}{2}$.
	Assume that  ${\bm A}_n(i,f) \ge \rho_1$ and ${\bm \lambda}(f) \ge \rho_2$ for all $n,i,f$, and that $\bm A_n, \bm \lambda$ and the initial estimates $	\widehat{\bm A}^0_n, \widehat{\bm \lambda}^0$ satisfy
	\begin{align*}
	    |\widehat{\bm A}^0_n(i,f)-{\bm A}_n(i,f)| &\le \delta_1:=\frac{4}{\rho_1(4+ \overline{D})}\\
	    |\widehat{\bm \lambda^0}(f)-{\bm \lambda}(f)|& \le \delta_2:= \frac{4}{\rho_2(4+ N\overline{D})}.
	\end{align*}
\end{Assumption}
Here, $\overline{D}_1$ can be understood as a measure for the `conditioning' of $\boldsymbol{A}_n$ under the KL divergence---a larger $\overline{D}_1$ implies more diverse columns of $\boldsymbol{A}_n$, and thus a better `condition number'.
In addition, $|\overline{D}_2|$ measures how far $\bm \lambda$ is away from the uniform distribution.
Using Assumption~\ref{ass:D}, we show that the EM algorithm can provably improve upon reasonable initial guesses (e.g., those given by {CNMF-SPA}), under some conditions.
To be specific, we have the following theorem:
\begin{theorem} \label{thm:em}
	Let $\delta_{\min} = \min(\delta_1,\delta_2)$.
	Assume that the following hold:
	\begin{align*}
	N &\ge \max\left(\frac{33\log(3SF/\mu)}{\rho_1\overline{D}_1},\frac{4\log(4SF^2/(3p\rho_2\mu))}{\overline{D}}\right),\\
	S &\ge \frac{192F^2\log(12NFI/\mu)}{p^2\rho_2^2\delta^2_{\min}},~\overline{D} \ge \max\left\{\frac{8-4\rho_1^2}{\rho_1^2},\frac{8-4\rho_2^2}{N\rho_2^2}\right\}.
	\end{align*}
	Then, under Assumption~\ref{ass:D},
	the EM algorithm in \cite{yeredor2019maximum} outputs $\widehat{\bm A}_n(i,f),\widehat{\bm \lambda}(f)$ that satisfy the following with a probability greater than or equal to $1-\mu$:
	\begin{align*}
	&	|\widehat{\bm A}_n(i,f)-\bm A_n(i,f)|^2 \le  \frac{48\log(12NFI/\mu)}{Sp\bm \lambda(f)}\leq \delta_1^2 ,\\
	&|\widehat{\bm \lambda}(f)-\bm \lambda(f)|^2 \le \frac{192F^2\bm \lambda(f)\log(12NFI/\mu)}{S}\leq \delta_2^2.
	\end{align*}
\end{theorem}
The proof is relegated to Appendix~\ref{app:thm_em} in the supplementary materials.
Our proof extends the analysis of a different EM algorithm proposed in \cite{zhang2014spectral} that is designed for learning the Dawid-Skene model in crowdsourcing. The EM algorithm there effectively learns a naive Bayes model when the latent variable admits a uniform distribution, i.e., $\bm \lambda(f)=1/F$ for all $f$.
The algorithm in \cite{zhang2014spectral} and the EM algorithm in \cite{yeredor2019maximum} are closely related, but the analysis in \cite{zhang2014spectral} does not cover the latter---since the EM in \cite{yeredor2019maximum} that does not restrict $\bm \lambda$ to be a uniform PMF.
Extending the analysis in \cite{zhang2014spectral} turns out to be nontrivial, which requires careful characterization for the `difficulty' introduced by $\bm \lambda$'s deviation from the uniform distribution.
Our definition of $\overline{D}_2$ and Theorem~\ref{thm:em} fill the gap.

\section{Experiments}
In this section, we use synthetic and real data experiments to showcase the effectiveness of the proposed algorithms.
\subsection{Baselines}\label{sec:baselines}
We use the CTD based method in \cite{kargas2017tensors} and the MLE-based method in \cite{yeredor2019maximum} as the major benchmarks, since they are general joint PMF estimation approaches as the proposed methods. Since the MLE based EM algorithm in \cite{yeredor2019maximum} may be sensitive to initialization, an iterative optimization algorithm designed for the same MLE formulation (cf. the alternating directions (AD) algorithm in \cite{yeredor2019maximum}) is used to initialize the EM in \cite{yeredor2019maximum}. This baseline is denoted as the {MLE-AD-EM} method in our experiments.

For the real data experiments, we test the proposed approaches on two core tasks in machine learning, namely, classification and recommender systems.
For the classification task, we use 4 different classifiers from the MATLAB
\textit{Statistics and Machine Learning Toolbox}, i.e., linear SVM, kernel SVM with radial basis function (RBF)), multinomial logistic regression and naive Bayes classifier.
For the user-item recommender systems, we use the celebrated \textit{biased matrix factorization} (BMF) method \cite{koren2009matrix} as a baseline. We also
compare with results obtained by \textit{global average} of the ratings, the
\textit{user average}, and the \textit{item average} for predicting
the missing entries in the user-item matrix.

All the algorithms are implemented in MATLAB 2018b and run on a desktop machine equipped with i7 3.40 GHZ CPU.
The iterative algorithms used in the experiments are stopped when the relative change in the objective function value is less than $10^{-6}$.

\begin{table}[t!]
	\centering
	\caption{MSE \& MRE Performance on Synthetic Data with $N=5,F=5,I=10,p=0.5,\varepsilon=0.1$.}

		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			\textbf{Algorithms} & \textbf{Metric} & $S=10^3$ & $S=10^4$ & $S=10^5$ & $S=10^6$ \\
			\hline
			CNMF-SPA [Proposed] & MSE & 0.0669 & 0.0270 & 0.0210 & \textbf{0.0204}  \\
			\hline
			CNMF-OPT [Proposed] & MSE & \textbf{0.0519} & \textbf{0.0242} & 0.0210 & \textbf{0.0204} \\
			\hline
			CNMF-SPA-EM [Proposed] & MSE & 0.0553 & 0.0247 & \textbf{0.0208} & \textbf{0.0204} \\
			\hline
			RAND-EM & MSE & 0.0734 & 0.0317 & 0.0289 & 0.0336 \\
			\hline
			CTD[Kargas \textit{et al.}]   & MSE &  0.1523 & 0.0256 & 0.0212 & 0.0205 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & MSE & 0.0909 & 0.0368 & 0.0266 & 0.0296 \\
			\hline
			\hline
			CNMF-SPA [Proposed] & MRE & 0.7965 & 0.3321 & 0.1052 & 0.0346  \\
			\hline
			CNMF-OPT [Proposed] & MRE & \textbf{0.6746} & \textbf{0.2327} & 0.0721 & 0.0239 \\
			\hline
			CNMF-SPA-EM [Proposed] & MRE & 0.6788 & 0.2341 & \textbf{0.0664} & \textbf{0.0219} \\
			\hline
			RAND-EM & MRE & 0.8078 & 0.2820 & 0.1845 & 0.2610 \\
			\hline
			CTD [Kargas \textit{et al.}] & MRE  & 0.9076 & 0.2920 & 0.0952 & 0.0329 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & MRE & 0.8363 & 0.3552 & 0.1384 & 0.1593 \\
			\hline
		\end{tabular}%
	\vspace{-.25cm}
	\label{tab:synthetic2}%
\end{table}%

\begin{comment}

\begin{table}[t!]
	\centering
	\caption{MSE \& MRE Performance on Synthetic Data with $N=5,F=5,I=10,p=0.5,\varepsilon=0.2$.}
	\resizebox{.9\linewidth}{!}{\Huge
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			\textbf{Algorithms} & \textbf{Metric} & $S=10^3$ & $S=10^4$ & $S=10^5$ & $S=10^6$ \\
			\hline
			CNMF-SPA [Proposed] & MSE & 0.0701 & 0.0261 & 0.0201 & 0.0200 \\
			\hline
			CNMF-OPT [Proposed] & MSE & \textbf{0.0532} & \textbf{0.0236} & 0.0207 & 0.0209 \\
			\hline
			CNMF-SPA-EM [Proposed] & MSE & 0.0612 & 0.0238 & 0.0208 & 0.0210 \\
			\hline
			RAND-EM & MSE & 0.0718 & 0.0264 & 0.0339 & 0.0272 \\
			\hline
			CTD[Kargas \textit{et al.}]   & MSE   & 0.1803 & 0.0250 & \textbf{0.0206} & \textbf{0.0207} \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & MSE & 0.0882 & 0.0332 & 0.0224 & 0.0292 \\
			\hline
			\hline
			CNMF-SPA [Proposed] & MRE & 0.8005 & 0.3796 & 0.2029 & 0.1724 \\
			\hline
			CNMF-OPT [Proposed] & MRE &  \textbf{0.7084} & \textbf{0.2530} & \textbf{0.1017} & 0.0541 \\
			\hline
			CNMF-SPA-EM [Proposed] & MRE & 0.7213 & 0.2545 & 0.1367 & 0.0911 \\
			\hline
			RAND-EM & MRE & 0.7975 & 0.2856 & 0.2883 & 0.2148 \\
			\hline
			CTD [Kargas \textit{et al.}] & MRE  & 0.9248 & 0.2992 & 0.1097 & \textbf{0.0320} \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & MRE & 0.8567 & 0.3168 & 0.1352 & 0.1895 \\
			\hline
		\end{tabular}%
	}
	\label{tab:synthetic_ep2}%
\end{table}%

\end{comment}

\begin{table}[t!]
	\centering
	\caption{MSE \& MRE Performance on Synthetic Data with $N=5,F=5,I=10,p=0.5,\varepsilon=0.3$.}
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			\textbf{Algorithms} & \textbf{Metric} & $S=10^3$ & $S=10^4$ & $S=10^5$ & $S=10^6$ \\
			\hline
			CNMF-SPA [Proposed] & MSE & 0.0880 & 0.0288 & 0.0238 & 0.0237 \\
			\hline
			CNMF-OPT [Proposed] & MSE & \textbf{0.0810} & \textbf{0.0247} & \textbf{0.0204} & 0.0210 \\
			\hline
			CNMF-SPA-EM [Proposed] & MSE & 0.0813 & 0.0249 & 0.0206 & 0.0217 \\
			\hline
			RAND-EM & MSE & 0.0883 & 0.0292 & 0.0216 & 0.0276 \\
			\hline
			CTD[Kargas \textit{et al.}]   & MSE   & 0.1534 & 0.0266 & \textbf{0.0204} & \textbf{0.0206} \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & MSE & 0.0911 & 0.0312 & 0.0313 & 0.0233 \\
			\hline
			\hline
			CNMF-SPA [Proposed] & MRE & 0.8933 & 0.4395 & 0.3201 & 0.3128 \\
			\hline
			CNMF-OPT [Proposed] & MRE & \textbf{0.8010} & 0.2675 & 0.1021 & 0.0784 \\
			\hline
			CNMF-SPA-EM [Proposed] & MRE & 0.8133 & \textbf{0.2658} & 0.1132 & 0.1151 \\
			\hline
			RAND-EM & MRE & 0.8214 & 0.3080 & 0.1390 & 0.2239 \\
			\hline
			CTD [Kargas \textit{et al.}] & MRE  & 0.9101 & 0.3298 & \textbf{0.1017} & \textbf{0.0346} \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & MRE & 0.8752 & 0.3593 & 0.2666 & 0.1413 \\
			\hline
		\end{tabular}%
	\vspace{-.25cm}
	\label{tab:synthetic_ep3}%
\end{table}%

\subsection{Synthetic-data Experiments} \label{sec:synthetic}
We consider $N=5$ RV's where each variable takes $I=10$ discrete values. The entries of the conditional PMF matrices (factor matrices) ${\bm A}_n\in \mathbb{R}^{I_n \times F}$ and the prior probability vector ${\bm \lambda} \in \mathbb{R}^{F}$ are randomly generated from a uniform distribution between 0 and 1 and the columns are rescaled to have unit $\ell_1$ norm with rank $F=5$.
We generate $S$ realizations of the joint PMF by randomly hiding each variable realization with observation probability $p=0.5$.
In particular, to test the performance of CNMF-SPA under various $\varepsilon$-separability conditions, we fix $M=3$ and manually ensure the $\varepsilon$-separability condition on $\bm H$ in \eqref{eq:constructedX} by fixing $\varepsilon \in\{0.1,0.3\}$.
To run the proposed algorithms, the low-order marginals are then estimated from the realizations via sample averaging. The {\it mean square error} (MSE) of the factors (see~\cite{huang2014non} for definition) and the {\it mean relative error} (MRE) of the recovered joint PMFs (see \cite{kargas2017tensors}) are evaluated.
We should mention that, ideally, MRE is more preferred for evaluation, but it is hard to compute (due to memory issues) for large $N$. MSE is an intermediate metric since it does not directly measure the recovery performance. Low MSEs imply good joint PMF recovery performance, but the converse is not necessarily true.
The results are averaged from 20 random trials.

Table \ref{tab:synthetic2}-\ref{tab:synthetic_ep3} present the MSE and MRE results for various values of $\varepsilon$ and $S$, where we manually control $\varepsilon$ when generating the $\boldsymbol{A}_n$'s. One can see that the proposed CNMF-SPA method works reasonably well for both $\varepsilon$'s under test.
It performs better when $\varepsilon$ is small, corroborating our analysis in Sec~\ref{sec:analysisSPA}. In addition, CNMF-SPA is quite effective for initializing the EM algorithm whereas randomly initialized EM (RAND-EM) and the MLE-AD-EM method from \cite{yeredor2019maximum} sometimes struggle to attain good performance. Notably, the proposed pairwise marginals based methods consistently outperform the three-dimensional marginals based method, i.e., CTD from \cite{kargas2017tensors}, especially when $S$ is small. This highlights the sample complexity benefits of the proposed methods. Another point is that when CNMF-SPA does not result in good accuracy, then the follow-up EM stage (i.e., CNMF-SPA-EM) does not perform well (e.g., when $S=10^3$). However, the proposed coupled factorization based approach, i.e., CNMF-OPT, still works well in this `sample-starved regime', showing the power of recoverability guaranteed problem formulation.

Table~\ref{tab:synthetic3}-\ref{tab:synthetic4} show the MSE results of the algorithms tested on cases where $N=15,I=10$ and $F=10$ for various values of $p$ and $S$---recall that a small $p$ means the data has many missing observations.
The MRE is not presented since instantiating $15$th-order tensors is prohibitive in terms of memory.
The results for MLE-AD-EM is not presented since the AD stage consumes a large amount of time under this high-dimensional PMF setting.
We do not manually impose $\varepsilon$-separability here, with the hope that the $\bm H$ matrix is likely $\varepsilon$-separable under large $N$ and $I$, as argued in Theorem \ref{thm:spabound} and Remark~\ref{rem:theorem_spa}.
One can see that the proposed CNMF-based methods exhibit better MSE performance compared to other baselines.
Also, the advantage of CNMF-OPT is more articulated under challenging settings (e.g., small $p$ and small $S$), where one can see that CNMF-OPT almost always outperforms the other algorithms.

Table~\ref{tab:runtime} shows the average runtime of the algorithms under different problem sizes. One can see that for the smaller case (cf. the left column), all the methods exhibit good runtime performance except for MLE-AD-EM. The elongated runtime of MLE-AD-EM is due to the computationally expensive AD iterations used for initialization. For the high-dimensional case where $N=15,F=I=10$, CNMF-SPA and EM are much faster due to their computationally economical updates.
Considering the MSE and MRE results in Tables \ref{tab:synthetic2}-\ref{tab:synthetic4}, CNMF-SPA and CNMF-SPA-EM seem to provide a preferred good balance between speed and accuracy.
Furthermore, the proposed CNMF-OPT offers the best accuracy for sample-starved cases---at the expense of higher computational costs.

\begin{comment}
We have presented the simulation results with $\varepsilon =0.1$ in Table \ref{tab:synthetic2} in Sec.~\ref{sec:synthetic} (the parameter $\varepsilon$ which is associated with separability condition is defined in Def. \ref{def:sep}). Here, we present the results by further increasing $\varepsilon$  to show the effect of this parameter on the performance of the proposed methods. The results are presented in Tables~\ref{tab:synthetic_ep2} and \ref{tab:synthetic_ep3}. One can notice that as $\varepsilon$ increases, the estimation error for CNMF-SPA increases as expected. However the proposed algorithms CNMF-OPT and CNMF-SPA still outperform the baselines especially when the number of samples is not very large ($S \le 10^5)$.
\end{comment}

\begin{table}[t!]
	\centering
	\caption{MSE Performance on Synthetic Data with $N=15,F=10,I=10,p=0.5$.}
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			\textbf{Algorithms} & \textbf{Metric} & $S=10^3$ & $S=10^4$ & $S=10^5$ & $S=10^6$ \\
			\hline
			CNMF-SPA [Proposed] & MSE & 0.1183 & 0.1030 & 0.1063 & 0.1041 \\
			\hline
			CNMF-OPT [Proposed] & MSE & \textbf{0.0218} & \textbf{0.0042} & \textbf{0.0022} & 0.0020 \\
			\hline
			CNMF-SPA-EM [Proposed] & MSE & 0.0894 & 0.0110 & 0.0056 & \textbf{0.0018} \\
			\hline
			RAND-EM & MSE & 0.0376 & 0.0112 & 0.0149 & 0.0069 \\
			\hline
			CTD[Kargas \textit{et al.}]   & MSE &  0.0329 & 0.0359 & 0.0404 & 0.0355 \\
			\hline
		\end{tabular}%
	\vspace{-0.25cm}
	\label{tab:synthetic3}%
\end{table}%

\begin{table}[t!]
	\centering
	\caption{MSE Performance on Synthetic Data for $N=15,F=10,I=10,p=0.2$.}
		\begin{tabular}{|c|c|c|c|c|c|}
			\hline
			\textbf{Algorithms} & \textbf{Metric} & $S=10^3$ & $S=10^4$ & $S=10^5$ & $S=10^6$ \\
			\hline
			CNMF-SPA [Proposed] & MSE & 0.2884 & 0.1181 & 0.0987 & 0.0996 \\
			\hline
			CNMF-OPT [Proposed] & MSE & 0.2277 & \textbf{0.0423} & \textbf{0.0062} & \textbf{0.0021} \\
			\hline
			CNMF-SPA-EM [Proposed] & MSE & \textbf{0.2274} & 0.0672 & 0.0165 & 0.0070 \\
			\hline
			RAND-EM & MSE & 0.2333 & 0.0782 & 0.0349 & 0.0313 \\
			\hline
			CTD[Kargas \textit{et al.}]   & MSE &  0.4015 & 0.0862 & 0.0081 & 0.0260 \\
			\hline
		\end{tabular}%
	\label{tab:synthetic4}%
	\vspace{-0.25cm}
\end{table}%

\begin{table}[t!]
\centering
   \begin{threeparttable}[b]

	\caption{Average Runtime (Sec.) of the Algorithms on Synthetic Data $S=10^5$ under Various $[N,F,I]$'s.}
		\begin{tabular}{|c|c|c|}
			\hline
			\textbf{Algorithms} & $[5,5,10]$ & $[15,10,10]$ \\
			\hline
			CNMF-SPA [Proposed] & \textbf{0.139} & \textbf{0.865}  \\
			\hline
			CNMF-OPT [Proposed] & 1.633 & 94.931 \\
			\hline
			CNMF-SPA-EM [Proposed] & 0.490 & 4.091 \\
			\hline
			RAND-EM & 0.353 & 3.217 \\
			\hline
			CTD [Kargas \textit{et al.}]   & 0.741 & 50.559 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] & 124.449 & $ 12.949\times 10^3~^\dagger$ \\
			\hline
		\end{tabular}%

	\label{tab:runtime

	      \begin{tablenotes}
     \item[$\dagger$] {\footnotesize This number is reported from a single trial due to the elongated runtime of the algorithm under this high-dimensional setting.}
   \end{tablenotes}
  \end{threeparttable}
	\vspace{-.25cm}
\end{table}%

\subsection{Real-data Experiments: Classification Tasks}
We conduct classification tasks on a number of UCI datasets (see \url{https://archive.ics.uci.edu}) using the joint PMF recovery algorithms and dedicated data classifiers as we detailed in Sec.~\ref{sec:baselines}.  The details of the datasets are given in Table~\ref{tab:uci}. Note that in the datasets, some of the features are continuous-valued. To make the joint PMF recovery methods applicable, we follow the setup in \cite{kargas2017tensors} to discretize those continuous features. We use $I_{\sf ave}=(1/N)\sum_{n=1}^N I_n$ to denote the averaged alphabet size of the discretized features for each dataset in Table~\ref{tab:uci} \footnote{ The detailed setup and source code are available at \url{https://github.com/shahanaibrahimosu/joint-probability-estimation}.}.
We split each dataset into training, validation and testing sets with a ratio of $50:20:30$.
For our approach, we estimate the joint PMF of the features and the label using the training set, and then predict the labels on the testing data by constructing an MAP predictor following \cite{kargas2017tensors}. The rank $F$ for the joint PMF recovery methods and the number of iterations needed for the methods involving EM are chosen using the validation set. We set $M=5$ for all datasets. We perform 20 trials for each dataset with randomly partitioned training/testing/validation sets, and report the mean and standard deviation of the classification accuracy.

Tables~\ref{tab:votes}-\ref{tab:mushroom} show the results for the UCI datasets `Votes', `Car', `Nursery' and `Mushroom', respectively.
 In all cases, one can see that the proposed approaches CNMF-SPA, CNMF-OPT and CNMF-SPA-EM are promising in terms of offering competitive performance for data classification.
 We should remark that our method is not designed for data classification, but estimating the joint PMF of the features and the label. The good performance on classification clearly support our theoretical analysis and the usefulness of the proposed algorithms.

 A key observation is that the CPD method in \cite{kargas2017tensors} does not perform as well compared to the proposed pairwise marginals based methods, especially when the number of data samples is small (e.g., the `Votes' dataset in Table~\ref{tab:votes}).  This echos our motivation for this work---using pairwise marginals is advantageous over using third-order ones when fewer samples are available.
 In terms of runtime, the proposed CNMF-SPA and CNMF-SPA-EM methods are faster than other methods by several orders of magnitude.

\begin{table}[t]
	\centering
	\caption{Details of UCI datasets.}

		\begin{tabular}{|l|r|r|r|r|}
			\hline
			\textbf{Dataset} & \textbf{$\#$ Samples} & \textbf{$\#$ Features} & \textbf{$I_{\text{avg}}$}  & \textbf{$\#$ Classes} \\
			\hline
			Votes & 435 & 17    & 2    & 2 \\
			\hline
			Car   & 1728  & 7     & 4     & 4 \\
			\hline

			Nursery & 12960 & 8    & 4    & 4 \\
			\hline
			Mushroom & 8124  & 22    & 6     & 2 \\
			\hline
		\end{tabular}%

	\label{tab:uci}%
		\vspace{-.25cm}
\end{table}%

\begin{table}[t!]
	\centering

	\caption{Classification Results on UCI Dataset `Votes'.}

			\begin{tabular}{|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{Accuracy (\%)} & \textbf{Time (sec.)} \\
			\hline
			CNMF-SPA [Proposed] & 90.07$\pm$1.30 & 0.005 \\
			\hline
			CNMF-OPT [Proposed] & \textbf{94.94$\pm$2.13} & 9.720 \\
			\hline
			CNMF-SPA-EM [Proposed] & 92.82$\pm$2.17 & 0.018 \\
			\hline
			CTD [Kargas \textit{et al.}] & 92.37$\pm$2.88 & 6.413 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] &90.99$\pm$4.18 & 10.556\\
			\hline
			SVM   & 94.34$\pm$1.68 & 0.057 \\
			\hline
			Logistic Regression & 94.05$\pm$1.24 & 0.301 \\
			\hline

			SVM-RBF & 92.21$\pm$2.66 & 0.012 \\
			\hline
			Naive Bayes & 90.46$\pm$1.98 & 0.032 \\
			\hline
	\end{tabular}

	\label{tab:votes}%
		\vspace{-.25cm}
\end{table}%

\begin{table}[t!]

	\centering
	\caption{Classification Results on UCI Dataset `Car'.}
				\begin{tabular}{|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{Accuracy (\%)} & \textbf{Time (sec.)} \\
			\hline
			CNMF-SPA [Proposed] & 70.31$\pm$2.17  & 0.007 \\
			\hline
			CNMF-OPT [Proposed] & 85.00$\pm$1.80 & 5.550 \\
			\hline
			CNMF-SPA-EM [Proposed] & \textbf{87.42$\pm$1.12} & 0.017\\
			\hline
			CTD [Kargas \textit{et al.}] & 84.21$\pm$1.27 & 3.601 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] &87.25$\pm$1.41 & 6.321\\
			\hline
			SVM   & 84.59$\pm$3.21 & 0.149 \\
			\hline
			Logistic Regression & 83.09$\pm$2.42 & 1.528 \\
			\hline

			SVM-RBF & 77.13$\pm$4.29 & 0.876 \\
			\hline
			Naive Bayes & 83.32$\pm$2.23 & 0.023 \\
			\hline
		\end{tabular}%
	\label{tab:uci_car}%
	\vspace{-.25cm}
\end{table}%
\begin{comment}
\begin{table}[t!]
	\centering
	\caption{Classification Results on UCI Dataset `Credit'}
	\resizebox{.8\linewidth}{!}{
		\begin{tabular}{|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{Accuracy (\%)} & \textbf{Time (s)} \\
			\hline
			CNMF-SPA [Proposed] & 84.47$\pm$2.12 & 0.008 \\
			\hline
			\multicolumn{1}{|l|}{CNMF-OPT [Proposed]} & \textbf{84.92$\pm$2.19} & 9.945 \\
			\hline
			CNMF-SPA-EM [Proposed] & 84.64$\pm$1.83 & 0.018 \\
			\hline
			CTD [Kargas \textit{et al.}] & 84.42$\pm$2.31 & 3.374 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] &84.52$\pm$1.90 & 23.334\\
			\hline
			SVM   & 84.18$\pm$2.04 & 0.026 \\
			\hline
			Linear Regression & 84.88$\pm$2.00 & 0.016 \\
			\hline
			Neural Net  & 84.20$\pm$3.20 & 0.107 \\
			\hline
			SVM-RBF & 81.11$\pm$2.25 & 0.009 \\
			\hline
			Naive Bayes & 84.17$\pm$2.07 & 0.016 \\
			\hline
		\end{tabular}%
	}
	\label{tab:uci_credit}%
\end{table}%
\end{comment}

\begin{table}[t!]

	\centering
	\caption{Classification Results on UCI Dataset `Nursery'.}
				\begin{tabular}{|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{Accuracy (\%)} & \textbf{Time (sec.)} \\
			\hline
			CNMF-SPA [Proposed] & 97.48$\pm$0.20 & 0.008 \\
			\hline
			CNMF-OPT [Proposed] & \textbf{98.16$\pm$0.22} & 6.719 \\
			\hline
			CNMF-SPA-EM [Proposed] & 98.04$\pm$0.35 & 0.041 \\
			\hline
			CTD [Kargas \textit{et al.}] & 97.70$\pm$0.22 & 4.318 \\
			\hline

			MLE-AD-EM [Yeredor \& Haardt] &98.04$\pm$0.33 & 8.920\\
			\hline
			SVM   & 97.42$\pm$0.79 & 0.706 \\
			\hline
			Logistic Regression & 98.02$\pm$0.15 & 5.019 \\
			\hline

			SVM-RBF & 70.92$\pm$4.70 & 0.738 \\
			\hline
			Naive Bayes & 98.03$\pm$0.20 & 0.046 \\
			\hline
		\end{tabular}
	\label{tab:nursery}%
	\vspace{-.25cm}
\end{table}%

\begin{table}[t!]
	\centering
	\caption{Classification Results on UCI Dataset `Mushroom'.}

		\begin{tabular}{|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{Accuracy (\%)} & \textbf{Time (sec.)} \\
			\hline
			CNMF-SPA [Proposed] & 91.86$\pm$6.33 & 0.025 \\
			\hline
			CNMF-OPT [Proposed] & 96.70$\pm$0.82 & 39.349 \\
			\hline
			CNMF-SPA-EM [Proposed] & \textbf{99.47$\pm$0.68} & 0.210 \\
			\hline
			CTD [Kargas \textit{et al.}] & 96.09$\pm$0.46 & 24.116 \\
			\hline
			SVM   & 97.60$\pm$0.28 & 34.443 \\
			\hline
			Logistic Regression & 96.75$\pm$0.62 & 2.875 \\
			\hline

			SVM-RBF & 97.16$\pm$0.41 & 1.523 \\
			\hline
			Naive Bayes & 94.60$\pm$0.56 & 0.063 \\
			\hline
		\end{tabular}%
	\label{tab:mushroom}%

		\vspace{-.25cm}
\end{table}%

\subsection{Real-data Experiments: Recommender Systems}
We test the approaches using the MovieLens 20M dataset \cite{harper2015movielens}. Again, following the evaluation strategy in \cite{kargas2017tensors}, we first round the ratings to the closest integers so that every movie's rating resides in $\{1,2,\ldots,5\}$. We choose different movie genres, namely, action, animation and romance, and select a subset of movies from each genre. Each subset contains 30 movies. Hence, for every subset, $N=30$, $Z_i$ represents the rating of movie $i$, and $Z_i$'s alphabet is $\{1,2,\ldots,5\}$. We predict the rating for a movie, say, movie $N$, by user $k$ via computing $\mathbb{E}[i_N|r_k(1),\ldots,r_k(N-1)]$, where $r_k(i)$ denotes the rating of movie $i$ by user $k$ (i.e., using the MMSE estimator). This can be done via estimating ${\sf Pr}(i_1,\ldots,i_N)$.
Note that we leave out the MLE-AD-EM algorithm for the MovieLens datasets due to its scalability challenges.

We create the validation and testing sets by randomly hiding  $20\%$ and $30\%$ of the dataset for each trial. The remaining $50\%$ is used for training (i.e., learning the joint PMF in our approach).
The results are taken from 20 random trials as before. In addition to the root MSE (RMSE) performance, we also report the \textit{mean absolute error} (MAE) of the predicted ratings, which is more meaningful in the presence of outlying trials; see the definitions of RMSE and MAE in \cite{kargas2017tensors}.

Table~\ref{tab:action}-\ref{tab:romance} present the results for three different subsets, respectively. One can see that the proposed methods are promising---their predictions are either comparable or better than the BMF approach. Note BMF is specialized for recommender systems, while the proposed approaches are for generic joint PMF recovery.
The fact that our methods perform better suggests that the underlying joint PMF is well captured by the proposed CNMF approach.

\begin{table}[t!]
	\centering
	\caption{Movie Recommendation Results on the Action set.  }
	\label{tab:action}%

		\begin{tabular}{|c|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{RMSE} & \textbf{MAE} & \textbf{Time (s)} \\
			\hline
			CNMF-SPA [Proposed] & 0.8497$\pm$0.0114 & 0.6663$\pm$0.0059 & 0.031 \\
			\hline
			CNMF-OPT [Proposed] & 0.8167$\pm$0.0035 & 0.6321$\pm$0.0040 & 70.018 \\
			\hline
			CNMF-SPA-EM [Proposed] & \textbf{0.7840$\pm$0.0025} & \textbf{0.5991$\pm$0.0031} & 2.424 \\
			\hline
			CTD [Kargas \textit{et al.}] & 0.8770$\pm$0.0088 & 0.6649$\pm$0.0076 & 52.253 \\
			\hline
			BMF   & 0.8011$\pm$0.0012 & 0.6260$\pm$0.0013 & 46.637 \\
			\hline
			Global Average  & 0.9468$\pm$0.0018 & 0.6956$\pm$0.0017 & -- \\
			\hline
			User Average & 0.8950$\pm$0.0010 & 0.6825$\pm$0.0010 & -- \\
			\hline
			Movie Average & 0.8847$\pm$0.0018 & 0.6982$\pm$0.0012 & -- \\
			\hline
		\end{tabular}%

	\vspace{-.25cm}
\end{table}%
\begin{table}[t!]
	\centering
	\caption{Movie Recommendation Results on the Animation set.}

		\begin{tabular}{|c|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{RMSE} & \textbf{MAE} & \textbf{Time (s)} \\
			\hline
			CNMF-SPA [Proposed] & 0.8705$\pm$0.0095 & 0.6798$\pm$0.0060 & 0.028 \\
			\hline
			CNMF-OPT [Proposed] & \textbf{0.8124$\pm$0.0031} & \textbf{0.6241$\pm$0.0041} & 61.018 \\
			\hline
			CNMF-SPA-EM [Proposed] & 0.8170$\pm$0.0075 & 0.6317$\pm$0.0086 & 2.424 \\
			\hline
			CTD [Kargas \textit{et al.}] & 0.8300$\pm$0.0053 & 0.6335$\pm$0.0029 & 48.253 \\
			\hline
			BMF   & 0.8408$\pm$0.0023 & 0.6553$\pm$0.0015 & 46.637 \\
			\hline
			Global Average  & 0.9371$\pm$0.0021 & 0.7042$\pm$0.0014 & -- \\
			\hline
			User Average & 0.8850$\pm$0.0009 & 0.6632$\pm$0.0011 & -- \\
			\hline
			Movie Average & 0.9027$\pm$0.0019 & 0.6900$\pm$0.0013 & -- \\
			\hline
		\end{tabular}%
	\label{tab:anim}%
	\vspace{-.25cm}
\end{table}%

\begin{table}[t!]
	\centering
	\caption{Movie Recommendation Results on the Romance set.}

		\begin{tabular}{|c|c|c|c|}
			\hline
			\textbf{Algorithm} & \textbf{RMSE} & \textbf{MAE} & \textbf{Time (s)} \\
			\hline
			CNMF-SPA [Proposed] & 0.9280$\pm$0.0066 & 0.7376$\pm$0.0076 & 0.032 \\
			\hline
			CNMF-OPT [Proposed] & 0.9076$\pm$0.0014 & 0.7123$\pm$0.0029 & 60.762 \\
			\hline
			CNMF-SPA-EM [Proposed] & \textbf{0.9057$\pm$0.0052} & \textbf{0.7106$\pm$0.0049} & 1.881 \\
			\hline
			CTD [Kargas \textit{et al.}] & 0.9498$\pm$0.0085 & 0.7416$\pm$0.0054 & 47.010 \\
			\hline
			BMF   & 0.9337$\pm$0.0007 & 0.7463$\pm$0.0009 & 31.823 \\
			\hline
			Global Average  & 1.0019$\pm$0.0007 & 0.8078$\pm$0.0008 & -- \\
			\hline
			User Average & 1.0195$\pm$0.0007 & 0.7862$\pm$0.0008 &  --\\
			\hline
			Movie Average & 0.9482$\pm$0.0007 & 0.7599$\pm$0.0007 &  --\\
			\hline
		\end{tabular}%
	\label{tab:romance}%
	\vspace{-.25cm}
\end{table}%

\section{Conclusion}

We proposed a new framework for recovering joint PMF of a finite number of discrete RVs from marginal distributions. Unlike a recent approach that relies on three-dimensional marginals, our approach only uses two-dimensional marginals, which naturally has reduced-sample complexity and a lighter computational burden.
We showed that under certain conditions, the recoverability of the joint PMF can be guaranteed---with only noisy pairwise marginals accessible. We proposed a coupled NMF formulation as the optimization surrogate for this task, and proposed to employ a Gram--Schmidt-like scalable algorithm as its initialization.
We showed that the initialization method is effective even under the finite-sample case (whereas the existing tensor-based method does not offer sample-complexity characterization). We also showed that the Gram--Schmidt-like algorithm can provably enhance the performance of an EM algorithm that is proposed from an ML estimation perspective for the joint PMF recovery problem.
We tested the proposed approach on classification tasks and recommender systems, and the results corroborate our analyses.

\color{black}

\bibliographystyle{IEEEtran}