# Robust Bayesian Tensor Factorization with Zero-Inflated Poisson Model and
Consensus Aggregation

###### Abstract.

Tensor factorizations (TF) are powerful tools for the efficient representation
and analysis of multidimensional data. However, classic TF methods based on
maximum likelihood estimation underperform when applied to zero-inflated count
data, such as single-cell RNA sequencing (scRNA-seq) data. Additionally, the
stochasticity inherent in TFs results in factors that vary across repeated
runs, making interpretation and reproducibility of the results challenging. In
this paper, we introduce Zero Inflated Poisson Tensor Factorization (ZIPTF), a
novel approach for the factorization of high-dimensional count data with
excess zeros. To address the challenge of stochasticity, we introduce
Consensus Zero Inflated Poisson Tensor Factorization (C-ZIPTF), which combines
ZIPTF with a consensus-based meta-analysis. We evaluate our proposed ZIPTF and
C-ZIPTF on synthetic zero-inflated count data and synthetic and real scRNA-seq
data. ZIPTF consistently outperforms baseline matrix and tensor factorization
methods in terms of reconstruction accuracy for zero-inflated data. When the
probability of excess zeros is high, ZIPTF achieves up to $2.4\times$ better
accuracy. Additionally, C-ZIPTF significantly improves the consistency and
accuracy of the factorization. When tested on both synthetic and real scRNA-
seq data, ZIPTF and C-ZIPTF consistently recover known and biologically
meaningful gene expression programs. All our data and code are available at:
https://github.com/klarman-cell-observatory/scBTF and
https://github.com/klarman-cell-observatory/scbtf_experiments.

.

###### 2010 Mathematics Subject Classification:

Primary: 11E76,11P05,12D15,14N10

**footnotetext: These authors contributed equally to this work

Daniel Chafamo 1,2∗, Vignesh Shanmugam 1,3,4∗, Neriman Tokcan1,5∗

1 Broad Institute of MIT and Harvard, Cambridge, MA, USA

2Klarman Cell Observatory, Broad Institute of MIT and Harvard, Cambridge, MA,
USA

3Department of Pathology, Brigham and Women’s Hospital, Boston, MA

4Harvard Medical School, Boston, MA, USA

5Department of Mathematics, University of Massachusetts Boston, Boston, MA

##  1\. Introduction

Tensors are multi-way arrays that extend matrices to higher dimensions and
provide a natural way to represent multidimensional data. Traditional matrix
methods _matricize_ tensors, limiting their ability to exploit the intrinsic
multi-way structure of the data [24]. Tensor factorization extends matrix
factorization to higher dimensions while preserving the said intrinsic
structure and enabling the discovery of complex interactions within the data.
Several variants of tensor factorization methods exist, among which
Candecomp/Parafac (CP) and Tucker are the most widely used [20]. Tensor
factorizations have found applications in fields such as computer vision,
neuroscience, genomics, recommender systems, and social network analysis [24,
32, 42, 14, 43, 20, 21].

Classic tensor factorization methods using maximum likelihood estimation (MLE)
can be unreliable when applied to sparse count data [10]. Bayesian Poisson
Tensor Factorization (BPTF)—a higher-order extension of Poisson matrix
factorization—is used to overcome the limitations of the MLE approach when
dealing with high-dimensional count data. BPTF provides advantages such as the
ability to incorporate prior knowledge, perform model selection, and quantify
uncertainty in parameter estimates [17, 37, 19]. Highly-dispersed count data
with excess number of zeros is common in various fields such as healthcare
(e.g., hospital readmissions), genomics (e.g., gene expression levels), social
sciences (e.g., user behaviors), and insurance claims [11, 39]. The Zero-
Inflated Poisson (ZIP) distribution is a better model for such data compared
to the Poisson distribution [27, 11, 15], and has been successfully used in
recommendation systems and other applications [39].

In addition to modeling the distribution of data and noise appropriately,
another issue to be addressed is the inherent randomness of tensor
factorization algorithms. This leads to varying results for multiple runs and
negatively impacts the interpretability and reproducibility [24, 41]. In this
paper, we propose a novel approach for stable tensor factorization which is
robust for high-dimensional sparse count data with excess zeros (Section 3).
We claim three main contributions:

  * •

We propose a novel factorization approach for high-dimensional sparse count
data with excess zeros, namely Zero Inflated Poisson Tensor Factorization
(ZIPTF), which utilizes the Bayesian ZIP model (Section 3).

  * •

To address the discussed randomness issue, we develop a meta-analysis method
that generalizes consensus matrix factorization [25] and incorporates novel
techniques to improve the stability and interpretability of the factorization
results (Section 3.4, Figure 1). We specifically focus on its integration with
ZIPTF, namely Consensus-ZIPTF (C-ZIPTF). Nonetheless, our method is
generalizable to other factorization approaches.

  * •

We provide an extensive evaluation on three different datasets: (1) synthetic
zero-inflated count tensors with increasing probability $\Phi$ of excess zeros
(Section 4.1); (2) synthetic multi-sample single-cell RNA sequencing (scRNA-
seq) data (Section 4.2); (3) real scRNA-seq dataset of immune cells stimulated
with interferon beta (Section 4.3). We compare ZIPTF and C-ZIPTF against
baseline matrix and tensor factorization methods. Our results indicate that
ZIPTF outperforms the baselines in terms of reconstruction accuracy for zero-
inflated data. Specifically, for $\Phi=0.8$, ZIPTF achieves an average
explained variance of 0.92, compared to a maximum of 0.38 achieved by the
baseline models. Additionally, C-ZIPTF significantly improves the consistency
and accuracy of the factorization results. Finally, both ZIPTF and C-ZIPTF
successfully capture biologically meaningful gene expression programs (GEPs)
and result in factors with higher Pearson correlations to known GEPs.

##  2\. Tensor preliminaries

This section presents the foundational concepts and notations for tensors,
with most of the notation borrowed from [24]. We denote the
$(i_{1},i_{2},\ldots,i_{N})$-th entry of an $N$\- way tensor
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ as
$\mathcal{X}_{i_{1}i_{2}\ldots i_{N}}$. The Frobenius norm of a tensor is
similar to the matrix Frobenius norm:

(2.1)
$\|\mathcal{X}\|_{F}=\sqrt{\sum_{i_{1}=1}^{I_{1}}\sum_{i_{2}=1}^{I_{2}}\ldots\sum_{i_{N}=1}^{I_{N}}\mathcal{X}_{i_{1}i_{2}\ldots
i_{N}}^{2}}.$

An $N$-way tensor $\mathcal{Y}$ is called a rank-1 tensor if it can be written
as outer product of $N$ vectors, i.e, $\mathcal{Y}=u^{(1)}\otimes
u^{(2)}\otimes\ldots\otimes u^{(N)}$ with $\mathcal{Y}{i_{1}i_{2}\ldots
i_{N}}=u{i_{1}}^{(1)}u_{i_{2}}^{(2)}\ldots u_{i_{N}}^{(N)}$. A rank $R\geq 1$
approximation to the tensor $\mathcal{X}\in\mathbb{R}^{I_{1}\times
I_{2}\ldots\times I_{N}}$ can given as:

(2.2)
$\mathcal{X}=\mathcal{\tilde{X}}+\mathcal{E}~{}\text{where}~{}\mathcal{\tilde{X}}=\sum_{r=1}^{R}a^{(1)}_{r}\otimes
a^{(2)}_{r}\otimes\dots\otimes a^{(N)}_{r},$

$A^{(i)}=[a^{(i)}_{1}\ldots a^{(i)}_{R}]\in\mathbb{R}^{I_{i}\times R},1\leq
i\leq N$ is the factor matrix along the $i-$th mode, and
$\mathcal{E}\in\mathbb{R}^{I_{1}\times I_{2}\ldots\times I_{N}}$. The
factorization given in Eqn. (3.9) is often referred to as the CP (Candecomp /
Parafac) decomposition which is a special case of Tucker decomposition (see
[24] for details). The approximation can be concisely expressed as
$\tilde{\mathcal{X}}=[[A^{(1)},~{}A^{(2)},\ldots,~{}A^{(N)}]].$ In this paper,
we impose a non-negativity constraint on factors to improve their
interpretability. The primary method for solving Eqn. (3.9) involves using the
maximum likelihood estimation (MLE) approach, which entails minimizing the
following error:

(2.3)
$\min_{A^{(1)},A^{(2)},\ldots,A^{(N)}}||\mathcal{X}-\tilde{\mathcal{X}}||_{F}.$

Iterative algorithms such as multiplicative updates, alternating least, and
gradient descent are commonly utilized for Eqn. (2.3) [1, 24, 41]. The MLE
approach often assumes Gaussian noise [24, 41].

##  3\. Bayesian tensor factorization and consensus aggregation

###  3.1. Bayesian Poisson tensor factorization

Traditional tensor factorization methods using MLE are unstable when applied
to zero-inflated count data [10]. Bayesian Poisson Tensor Factorization (BPTF)
extends the Poisson Matrix Factorization method to higher dimensions and
utilizes Bayesian inference to obtain a point estimate and offers benefits
such as uncertainty quantification, realistic noise assumptions, and
principled inclusion of prior information [17, 19, 37, 44]. This section
presents a general framework for BPTF with Variational Inference (VI) for
high-dimensional count data. Let $\mathcal{X}\in\mathbb{R}^{I_{1}\times
I_{2}\ldots\times I_{N}}$ be the observed count data drawn from the Poisson
distribution and with the CP decomposition as given in Eqn. (3.9). Let
$I=i_{1}i_{2}\ldots i_{N}\in\overline{I}=\\{i_{1}i_{2}\ldots i_{N}~{}:~{}1\leq
i_{j}\leq I_{j},~{}1\leq j\leq N\\}$, then

(3.1) $X_{I}\approx
Poisson(\lambda_{I})~{}~{}~{}\text{where}~{}\mathcal{X}_{I}\approx\tilde{\mathcal{X}}_{I}=\sum_{r=1}^{R}a^{(1)}_{i_{1}r}a^{(2)}_{i_{2}r}\ldots
a^{(N)}_{i_{N}r}\approx\lambda_{I}.$

BPTF uses Gamma priors to regularize the estimation of the latent factors [8,
46, 5]. The Gamma distribution, which is characterized by a shape parameter
$\alpha>0$ and a rate parameter $\alpha\beta>0$, is employed as a sparsity-
inducing prior [37, 17, 8]. Then for each $a^{(k)}_{jr}$ in Eqn. (3.1), we
have:

(3.2) $a^{(k)}_{jr}\approx Gamma(\alpha,\alpha\beta^{(k)}),~{}1\leq k\leq N,$

with the expectation $E[a^{(k)}_{jr}]=\frac{1}{\beta^{(k)}}$ and
$Var[a^{(k)}_{jr}]=\frac{1}{\alpha{\beta^{(k)}}^{2}}$. The posterior
distribution given by
$P(A^{(1)},A^{(2)},\ldots,A^{(N)}|\mathcal{X},\mathcal{H})$ is intractable due
to the inability to compute the evidence, given a model hyperparameter set
$\mathcal{H}=\\{\alpha,B^{(1)},B^{(2)},\ldots,B^{(N)}\\}$ [5]. BPTF uses VI
and assumes a variational family of distributions
$\mathcal{Q}_{V}=\mathcal{Q}(A^{(1)},A^{(2)},\ldots,A^{(N)};V^{(1)},\ldots,V^{(N)})$
which is indexed by a set of variational parameters $V^{(k)},1\leq k\leq N$
[5, 4]. We employ a fully-factorized mean-field approximation assuming that
$\mathcal{Q_{V}}(A^{(1)},A^{(2)},\ldots,A^{(N)})=\prod_{k=1}^{N}\mathcal{Q_{V}}(A^{(k)};V^{(k)}),$
where
$\mathcal{Q}(a^{(k)}_{jr};V^{(k)}_{jr})=Gamma(a^{(k)}_{jr};\gamma^{(k)}_{jr},\delta^{(k)}_{jr}),1\leq
k\leq N.$ The variational family $Q$ used here is similar to the one employed
in Bayesian Poisson Matrix Factorization [8, 16, 33]. BPTF fits variational
parameters by minimizing the Kullback-Leibler (KL) divergence between the true
posterior distribution and $\mathcal{Q_{V}}$, which is equal to maximizing the
evidence lower bound (ELBO) [5, 4, 17]:

(3.3)
$ELBO(V)=E_{Q_{V}}[log(P(\mathcal{X},A^{(1)},A^{(2)},\ldots,A^{(N)}|\mathcal{H})]+H(Q_{V}).$

where $H(Q_{V})$ is the entropy for $Q_{V}.$ Coordinate ascent algorithms are
commonly used to maximize the ELBO by iteratively optimizing each variational
parameter while fixing the others until convergence, monitored by the relative
change in the ELBO [4, 5]. From Eqn. (3.1), we have the total
$n=\sum_{I\in\overline{I}}\mathcal{X}_{I}\approx~{}Poisson(\Lambda)$ where
$\Lambda=\sum_{I\in\overline{I}}\lambda_{I}.$ We can use the Poisson-
Multinomial connection to express $\mathcal{X}$ given $n$ as
$Multinomial(n,\pi)$ where $(\pi)_{I}=\frac{\lambda_{I}}{\Lambda}$, and update
variational parameters using this auxiliary distribution [5, 37, 8, 26]:

(3.4) $\displaystyle\gamma^{(k)}_{jr}$ $\displaystyle=$
$\displaystyle\alpha+\sum_{\begin{subarray}{c}i_{1}i_{2}\ldots
i_{N}\in\overline{I}\\\ i_{k}=j\end{subarray}}\mathcal{X}_{i_{1}i_{2}\ldots
i_{n}}\frac{\mathbb{G}_{Q_{V}}\big{[}\prod_{s=1}^{N}a^{(s)}_{i_{s}r}\big{]}}{\sum_{t=1}^{R}\mathbb{G}_{Q_{V}}\big{[}\prod_{s=1}^{N}a^{(s)}_{i_{s}t}\big{]}},$
(3.5) $\displaystyle\delta^{(k)}_{jr}$ $\displaystyle=$
$\displaystyle\alpha\beta^{(k)}+\sum_{i_{1}i_{2}\ldots
i_{N}\in\overline{I}}E_{Q_{V}}\big{[}\prod_{1\leq s\neq k\leq
N}a^{(s)}_{i_{s}r}~{}\big{]},$

where $E_{Q_{V}}[.]$ and $\mathbb{G}_{Q_{V}}=exp(E_{Q_{V}}[log(.)])$ denote
arithmetic and geometric expectations. Since $Q_{V}$ is fully factorized, the
expectations in Equations (3.4) and (3.5) can be expressed as a product of
individual expectations [5]. Specifically, for $a^{(s)}_{i_{s}r},$

(3.6)
$E_{Q_{V}}[a^{(s)}_{i_{s}r}]=\frac{\gamma^{(s)}_{i_{s}r}}{\delta^{(s)}_{i_{s}r}}~{}~{}\text{and}~{}~{}\mathbb{G}_{Q_{V}}[a^{(s)}_{i_{s}r}]=\frac{exp(\Psi(\gamma^{(s)}_{i_{s}r}))}{\delta^{(s)}_{i_{s}r}},$

where $\Psi$ is the digamma function (logarithmic derivative of the gamma
function). An empirical Bayes approach can be used to update the
hyperparameters $\beta^{(k)},1\leq k\leq N$, in conjunction with the
variational parameters [8, 5]:

(3.7)
$\beta^{(k)}=\big{(}\sum_{j=1}^{I_{j}}\sum_{r=1}^{R}\mathbb{E}_{Q_{V}}[a^{(k)}_{jr}]\big{)}^{-1}.$

The variational inference algorithm for BPTF is fully specified by the set of
update equations Equations (3.4), (3.5), and (3.7).

###  3.2. Zero-inflated Poisson tensor factorization (ZIPTF)

Poisson models may not always be sufficient to model count data with excess
zeros, and zero-inflated models can often provide a better fit [27, 39]. The
Zero-Inflated Poisson (ZIP) model assumes that the counts in the tensor
$\mathcal{X}$ can be modeled as a mixture of a point mass at zero and a
Poisson distribution with parameter $\lambda$. Let $\mathcal{X}$ be a count
data in $\mathbb{R}^{I_{1}\times I_{2}\ldots\times I_{N}}.$ We define the
index set $\overline{I}$ as the collection of all possible indices, i.e.,
$\overline{I}=\\{i_{1}i_{2}\ldots i_{N}~{}:~{}1\leq i_{j}\leq I_{j},~{}1\leq
j\leq N\\}.$ We say $\mathcal{X}$ has Zero-inflated Poisson (ZIP) distribution
if for every $I\in\overline{I}:$

(3.8)
$P(\mathcal{X}_{I}=x_{I})=p_{I}{\mathbb{1}_{x_{I}=0}}+(1-p_{I})\frac{e^{-\lambda}\lambda^{x_{I}}}{x_{I}!},$

where the outcome variable $x_{I}$ has non-negative integer values,
$\lambda_{I}$ is the expected Poisson count, and $p_{I}$ is the probability of
extra zeros [27]. As an abbreviation, we write it as $\mathcal{X}_{I}\sim
ZIP(\lambda_{I},p_{I}).$ The ZIP can be considered as the _product_ of a
Poisson random variable $\mathcal{Y}_{I}\sim Poisson(\lambda_{I})$ and an
independent Bernoulli variable $\Phi_{I}\sim Bernoulli(p_{I})$ [11]. The
Bernoulli variable $\Phi_{I}$ takes the value of $1$ when $\mathcal{X}_{I}$ is
equal to $0$, due to the Bernoulli component, and takes the value of $0$
otherwise.

We consider the low rank $R\geq 1$ decomposition of the zero-inflated count
tensor $\mathcal{X}:$

(3.9) $\mathcal{X}\approx\sum_{r=1}^{R}a^{(1)}_{r}\otimes
a^{(2)}_{r}\otimes\dots\otimes a^{(N)}_{r}.$

Hence, for $I=i_{1}i_{2}\ldots i_{N}$, the reconstruction
$\sum_{r=1}^{R}a^{(1)}_{i_{1}r}a^{(2)}_{i_{2}r}\ldots a^{(N)}_{i{N}r}$ can be
interpreted as the mean of the distribution from which the observed count
$\mathcal{X}_{I}$ is assumed to be sampled. Then we have:

(3.10) $\mathcal{X}_{I}\sim
ZIP(\lambda_{I}=\sum_{r=1}^{R}a^{(1)}_{i_{1}r}a^{(2)}_{i_{2}r}\ldots
a^{(N)}_{i{N}r},p_{I}).$

###  3.3. Variational Inference for ZIPTF

For given position $I=i_{1}i_{2}\ldots i_{N},$ we consider the rank $R$
decomposition in Eqn. (3.10). In Bayesian Poisson factorizations, the Gamma
distribution is utilized as a prior to induce sparsity, and it is assumed that
each latent factor matrix $A^{(k)}=[a^{(k)}_{1}\ldots
a^{(k)}_{R}]\in\mathbb{R}_{+}^{I_{k}\times R},~{}1\leq k\leq N$, follows a
Gamma distribution [8, 37]. Therefore, for each $a^{(k)}_{jr}$ in Eqn. (3.10),
we have:

(3.11) $a^{(k)}_{jr}\sim Gamma(\alpha^{(k)},\beta^{(k)}),\quad 1\leq k\leq N,$

where $\alpha^{(k)}>0$ and $\beta^{(k)}>0$ represent the shape and rate
parameters of the distribution, with the expectation
$E[a^{(k)}_{jr}]=\frac{\alpha^{(k)}}{\beta^{(k)}}$ and
$Var[a^{(k)}_{jr}]=\frac{\alpha^{(k)}}{{\beta^{(k)}}^{2}}$. Additionally, for
ZIP models a latent variable $\xi$ is introduced to capture the hidden state
of the probability of extra zeros which specify $\Phi\sim Bernoulli(p_{I})$
[39, 5]. Let $S(.)$ denote the logistic sigmoid function, given by
$S(x)=\frac{1}{1+e^{-x}}$, then:

(3.12) $\xi=S(\zeta)~{}~{}\text{where}~{}~{}\zeta\sim Normal(\mu,\sigma).$

Let $Z=\\{A^{(1)},A^{(2)},\ldots,A^{(N)},\Phi\\},$ consider the posterior
distribution $P(Z|\mathcal{X},\mathcal{H}),$ given a model hyperparameter set
$\mathcal{H}=\\{\alpha^{(1)},\beta^{(1)},\alpha^{(2)},\ldots,\beta^{(2)},\ldots,\alpha^{(N)},\beta^{(N)},\mu,\sigma\\}.$

Variational inference approximates the true posterior distribution using a
family of probability distributions $\mathcal{Q}$ over hidden variables [5].
This family of distributions is characterized by free parameters, and the key
assumption is that each latent variable is independently distributed given
these parameters. We assume a variational family of distributions
$\mathcal{Q}$ indexed by a set of variational parameters
$V=\\{\gamma^{(1)},\delta^{(1)},\gamma^{(2)},\delta^{(2)},\dots,\gamma^{(N)},\delta^{(N)},\overline{\mu},\overline{\sigma}\\}$
where $(\gamma^{(k)},\delta^{(k)})$ are variational shape and rate parameters
of the Gamma distribution for the latent factor along the $k-$th mode, and
$(\overline{\mu},\overline{\sigma})$ are the variational parameters for
$\zeta$. We use a fully factorized mean-field approximation [5] and the
variational distribution factors as the following:

(3.13)
$\mathcal{Q}(A^{(1)},A^{(2)},\ldots,A^{(N)},\Phi)=\mathcal{Q}(\Phi;\overline{\mu},\overline{\sigma})\prod_{k=1}^{N}\mathcal{Q}(A^{(k)};\gamma^{(k)},\delta^{(k)}).$

where $a^{(k)}_{jr}\sim Gamma(\gamma^{(k)}_{jr},\delta^{(k)}_{jr})$ and
$\Phi_{I}\sim Bernoulli(S(\zeta))~{}\text{for}~{}\zeta\sim
Normal(\overline{\mu},\overline{\sigma})$. The goal is to choose a member
$q^{*}$ of the variational family variational distributions which minimizes
the Kullback-Leibler (KL) divergence of the exact posterior from
$\mathcal{Q}$:

(3.14)
$q^{*}(Z)=\arg\min_{q(Z)\in\mathcal{Q}}D_{KL}\left(q(Z)\middle\|P(Z|\mathcal{X},\mathcal{H})\right).$

Upon examining the KL divergence, we encounter a significant challenge: it
involves the true posterior distribution $P(Z|\mathcal{X},\mathcal{H})$, which
is not known. Nevertheless, we can rewrite the KL divergence as follows:

(3.15) $\displaystyle
D_{KL}\left(q(Z)\middle\|P(Z|\mathcal{X},\mathcal{H})\right)$ $\displaystyle=$
$\displaystyle\int
q(Z)\log\left(\frac{q(Z)}{P(Z|\mathcal{X},\mathcal{H})}\right)dZ~{}~{}~{}~{}~{}~{}~{}~{}`$
(3.16) $\displaystyle=$ $\displaystyle\int
q(Z)\log\left(\frac{q(Z)P(\mathcal{X},\mathcal{H})}{P(Z,\mathcal{X},\mathcal{H})}\right)dZ$
(3.17) $\displaystyle=$
$\displaystyle\log\left(P(\mathcal{X},\mathcal{H})\right)\int q(Z)dZ-\int
q(Z)\log\left(\frac{P(Z,\mathcal{X},\mathcal{H})}{q(Z)}\right)dZ~{}~{}~{}~{}~{}~{}~{}~{}$
(3.18) $\displaystyle=$
$\displaystyle\log\left(P(\mathcal{X},\mathcal{H})\right)-\int
q(Z)\log\left(\frac{P(Z,\mathcal{X},\mathcal{H})}{q(Z)}\right)dZ.~{}~{}~{}$

The second term in Eqn. (3.18) is called Evidence Lower Bound (ELBO). We know
that the KL divergence is non-negative, therefore,
$\log\left(P(\mathcal{X},\mathcal{H})\right)\geq~{}\text{ELBO}(q(Z))=\int
q(Z)\log\left(\frac{P(Z,\mathcal{X},\mathcal{H}}{q(Z)}\right)dZ.$

(3.19) $\displaystyle\text{ELBO}(q(Z))$ $\displaystyle=$ $\displaystyle\int
q(Z)\log\left(P(Z,\mathcal{X},\mathcal{H})\right)dZ-\int
q(Z)\log\left(q(Z)\right)dZ$ (3.20) $\displaystyle=$ $\displaystyle
E_{q(Z)}[\log\left(P(\mathcal{X},Z,\mathcal{H}\right)]-E_{q(Z)}[\log q(Z)].$

The evidence lower bound serves as a transformative tool that converts
intractable inference problems into optimization problems that can be tackled
using gradient-based methods [5].

Coordinate ascent algorithms are frequently employed in maximizing the
evidence lower bound (ELBO)[5, 39]. However, these algorithms require tedious
gradient calculations and may not scale well for very large datasets [18, 34].
Closed-form coordinate-ascent updates are applicable to conditionally
conjugate exponential family models, but they necessitate analytic computation
of various expectations for each new model[18, 34].

Stochastic Variational Inference (SVI) [18] offers a more efficient algorithm
by incorporating stochastic optimization [35]. This technique involves
utilizing noisy estimates of the gradient of the objective function. To
maximize the evidence lower bound (ELBO), we employ a stochastic optimization
algorithm known as the Black Box Inference Algorithm [34]. This algorithm
operates by stochastically optimizing the variational objective using Monte
Carlo samples from the variational distribution to compute the noisy gradient
(see Section 2, [34] for details). By doing so, it effectively alleviates the
burden of analytic computations and provides a more efficient approach to ELBO
maximization.

Figure 1. Overview of the consensus meta-analysis approach discussed in
Section 3.4 for the 3-way tensor $\mathcal{X}$.

###  3.4. Generic consensus-based tensor factorization

Selecting the number of components in tensor factorization is challenging [24,
1]. The dependence on initial guesses for latent factors can lead to
substantially different factor sets across repeated runs, making it difficult
to interpret the results [24, 1, 41]. We typically select the minimum value of
$R$ in Eqn. (3.9) that maximizes the explained variance of the approximation,
defined as follows:

(3.21) $\text{explained
variance}=1-\frac{||\mathcal{X}-\tilde{\mathcal{X}}||_{F}}{||\mathcal{X}||_{F}}.$

Our goal is not solely to improve the explained variance, but also to ensure
the interpretability and stability of the factors. We generalize the consensus
meta-analysis approach, which has been previously used for matrix
factorization [25], and include novel techniques to enhance the stability. The
overview of the proposed pipeline is depicted in Figure 1. In the remainder of
this section, we will refer to the steps ①-⑤ given in the figure.

Running a generic rank $R$ factorization given in Eqn. (3.9) for
$\mathcal{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ with
$M$ different random seeds yields the sets of non-negative factor matrices
$\\{(A^{(1)})_{m},(A^{(2)})_{m},\ldots,(A^{(N)})_{m}\\},~{}1\leq m\leq M$
(Step ①). For a chosen modality $k~{}(1\leq k\leq N$), we can aggregate and
normalize the factor matrices from independent runs (Step ②):

(3.22)
$\overline{A^{(k)}}=\Big{[}\frac{(A^{(k)})_{1}}{||(A^{(k)})_{1}||_{F}}~{}\frac{(A^{(k)})_{2}}{||(A^{(k)})_{2}||_{F}}~{}\ldots\frac{(A^{(k)})_{M}}{||(A^{(k)})_{M}||_{F}}\Big{]}\in
I_{k}\times(R\times M).$

The cophenetic correlation coefficient, commonly used to select ranks for
matrix factorizations [7], assumes one-to-one mapping between features and
factors based on maximum loadings. However, this assumption may not be valid
when a feature contributes significantly to multiple factors.

Our method for selecting the rank and evaluating factorization stability
involves clustering column factors of aggregated matrices and fixing the
initial guess to ensure reliability. Initially, we perform K-means clustering
[31] on the columns of the aggregated factor matrix $\overline{A^{(k)}}$ with
$K=R$ (Step ③). The resulting cluster sets are given as
$C^{(k)}_{i}=\\{\text{columns of}~{}\overline{A^{(k)}}~{}\text{assigned to
cluster}~{}i\\},1\leq i\leq R.$ The Local Outlier Factor algorithm [6] is used
to remove outliers by considering the local density deviation of a data point
compared to its neighbors. We evaluate the _goodness_ of the clustering with
the silhouette coefficient [36], computed as (b-a)/max(a,b), where a is the
average intra-cluster distance, and b is the average inter-cluster distance.
The silhouette coefficient ranges from -1 to 1, with higher values indicating
more coherence. After clustering, we obtain the consensus factors
$f^{(k)}_{C_{i}}$, where $1\leq i\leq R$, by computing the median value of the
factors in each cluster (Step ④) and form the consensus matrix:

(3.23) $\overline{A^{(k)}}_{C}=[f^{(k)}_{C_{1}}f^{(k)}_{C_{2}}\ldots
f^{(k)}_{C_{R}}]\in\mathbb{R}^{I_{k}\times R},1\leq k\leq N.$

We perform the decomposition using $\overline{A^{(k)}}_{C}$ as the fixed
initial guess for the $k$-th modality to obtain the final factor matrices
(Step ⑤).

Notice that if ZIPTF is employed as the factorization method in Step ①
described above, we refer to the resulting factorization as C-ZIPTF.

##  4\. Experiments

Here we present results showing the superior performance of C-ZIPTF across
multiple evaluation metrics. We implemented C-ZIPTF in Python, using the
probabilistic programming language Pyro [3]. Our presentation in Section 3.1
focused on Bayesian tensor factorization frameworks that utilize Poisson and
Zero-Inflated Poisson based models. However, our implementation is designed to
be more versatile and can accommodate different types of noise models. We
conduct three different evaluations to assess the performance of our proposed
method. First, we compare ZIPTF with alternative factorization methods on
simulated tensors with known factors and ZIP noise and evaluate the benefits
of using a ZIP model and the inclusion of the consensus approach (Section
4.1). Second, we evaluate the performance of our method on simulated single-
cell RNA sequencing data and compare it with other matrix and tensor
factorization methods at the task of recovering identity and activity gene
expression programs (GEPs), Section 4.2. Finally, we demonstrate the ability
of our method to capture biologically meaningful gene expression programs by
applying it to a real single-cell RNA sequencing dataset of immune cells
stimulated with interferon beta (Section 4.3).

Figure 2. ZIPTF compared to alternative factorization methods on a synthetic
tensor with known factors and ZIP noise, and stability comparison between
ZIPTF and C-ZIPTF. (a) We calculated the explained variance of ZIPTF and
alternative methods for different levels of extra zeros. (b) Cosine similarity
between factors obtained on repeat runs for ZIPTF and C-ZIPTF. (c) Cosine
similarity between inferred factors and original factors for ZIPT and C-ZIPT.

###  4.1. Synthetic tensor experiment

To evaluate the performance of C-ZIPTF on synthetic data, we generate tensors
using known factors and Poisson noise with varying degrees of zero inflation
and measure the accuracy of different methods at recovering the original
factors. To generate a tensor $\mathcal{T^{\prime}}\in\mathbb{R}^{I\times
J\times K}$, we first create three factor matrices
$A\in\mathbb{R}_{+}^{I\times R}$, $B\in\mathbb{R}_{+}^{J\times R}$, and
$C\in\mathbb{R}_{+}^{K\times R}$, with elements drawn from a Gamma
distribution with shape $\alpha=3$ and rate $\beta=0.3$, where $R$ is the
desired true rank. We then construct a tensor $\mathcal{T}$ by taking the sum
of the outer product of the corresponding columns of the matrices, i.e.,
$\mathcal{T}=[[A,B,C]]$. Finally, we generate the tensors
$\mathcal{T}^{\prime}$ by sampling from a ZIP distribution with mean
$\mathcal{T}$ and a given probability of extra zeros, denoted by $\Phi$ in
Section 3.2.

Zero-inflated Poisson model results in higher explained variance For the first
experiment, we ran ZIPTF without consensus aggregation to evaluate the
advantages of using the ZIP model alone. For comparison, we considered three
alternative tensor factorization methods: Non-Negative CP decomposition via
Alternating-Least Squares (NNCP-ALS) [24], Bayesian Tensor Factorization with
Truncated Gaussian model (TGTF)[17], and Bayesian Tensor Factorization with
Gamma Poisson model (GPTF)[37]. We conducted 20 trials, generating a new
simulated tensor $\mathcal{T}^{\prime}$ of shape $10\times 20\times 300$ and
rank 9 each time and running each factorization method on the tensor for a
fixed maximum number of iterations ($max\\_iter=1000$). We evaluated the
performance of the methods using the explained variance (3.21) of the
approximation generated by each factorization. ZIPTF consistently outperformed
both NNCP-ALS and the Bayesian Tensor Factorization methods without the ZIP
model, as shown in Figure 2(a). At a zero probability of excess zeros, all
four methods showed similar and nearly perfect explained variance. However, as
the excess zero level increased, the performance of the other methods
deteriorated rapidly. At the highest probability of excess zeros simulated,
$\Phi=0.8$, the average explained variance of the ZIPTF approximation was
0.974, with a 95% confidence interval (CI) [0.962, 0.987], about $2.4\times$
better than the second highest explained variance of 0.338, 95% CI [0.334,
0.342] achieved by the Gamma Poisson model. We also note that the difference
in explained variance between NNCP-ALS and the Bayesian methods other than
ZIPTF is minimal compared to the difference to ZIPTF. This indicates that the
superiority of ZIPTF arises from using the appropriate noise model.

Consensus aggregation leads to more consistent factorization After
demonstrating ZIPTF’s superior performance in modeling zero-inflated count
data, we examine the benefits of consensus aggregation. We generate tensors of
shape $40\times 20\times 2000$ and rank 9 with known factors and Zero-Inflated
Poisson noise as described above and evaluate the recovered factors by running
ZIPTF with and without consensus aggregation. For this experiment, we fix the
probability of excess zeros $\Phi=0.6$. We compare the internal consistency of
factors obtained from multiple runs of the decompositions. For our simulated
tensor $\mathcal{T}^{\prime},$ assume that we have two rank R approximations
$[[A,B,C]]$ and $[[D,E,F]]$ corresponding to different randomly initialized
runs. To measure the similarity between factorizations, we calculate:

(4.1) $\textit{cosine
score}([[A,B,C]],[[D,E,F]])=\frac{1}{R}\sum_{i=1}^{R}\max_{1\leq j\leq
R}~{}\cos(a_{i},d_{j})\cos(b_{i},e_{j})\cos(c_{i},f_{j}).\ $

We evaluate the similarity of factors recovered from 20 randomly initialized
runs of both ZIPTF and C-ZIPTF using the cosine score given in Eqn. 4.1. We
observe that the factors recovered from C-ZIPTF are more consistent with one
another compared to those recovered from ZIPTF, as shown in Figure 2(b). The
consensus approach makes C-ZIPTF more robust, reducing the impact of the
inherent stochasticity of the factorization process and resulting in a more
stable set of factors.

Consensus aggregation leads to more accurate recovery of original factors We
assess the accuracy of both ZIPTF and C-ZIPTF in recovering the original
factors used to create the tensor with $\Phi=0.6$. We perform 20 randomly
initialized runs of each method and compare the recovered factors to the
original factors using the cosine score. Figure 2(c) demonstrates that C-ZIPTF
outperforms ZIPTF in recovering the original factors.

###  4.2. Synthetic single-cell RNA-Seq data analysis

We test the performance of C-ZIPTF on single-cell RNA sequencing (scRNA-seq)
data, which is prone to zero inflation due to technical limitations that
result in dropout events [28]. To evaluate the effectiveness of C-ZIPTF, we
compared its performance with other matrix and tensor factorization methods
using a synthetic scRNA-seq dataset. We used the Splatter simulation framework
[45] which was adapted to Python in a previous study [25] to generate the
synthetic scRNA-seq dataset.

Figure 3. The average Pearson correlation between the true gene expression
programs (GEPs) used in the simulation and the inferred GEPs obtained from
various factorization methods. The results are presented for three different
signal intensity levels (0.25, 0.5, 0.75), which are indicated by the mean
log2 fold change (log2FC) of simulated differentially expressed genes.

The simulation framework utilizes a Gamma-Poisson hierarchical model with
hyper-parameters estimated from real data. Technical dropouts are simulated by
randomly replacing some of the simulated counts with zeros using a Bernoulli
distribution. The complete details of the simulation framework and parameters
used are provided in Section 6. The synthetic dataset consists of 15,000 cells
and 5,000 genes from six donors, with five gene expression programs defining
cell type identities and three gene expression programs defining donor-
specific activity. We evaluated the performance of various factorization
methods, including ZIPTF, C-ZIPTF, Non-negative Matrix Factorization (NMF)
[29], Consensus NMF (CNMF) [25], and NNCP-ALS, with the goal of recovering the
eight gene expression programs embedded in the synthetic scRNA-seq dataset.
For NMF and CNMF, the rank 8 decomposition is performed with a maximum of 1000
iterations for convergence after normalizing the cell-by-gene count matrix to
counts per million (CPM). For the tensor-based approaches, we construct the
observed tensor by pseudobulking the cell-by-gene counts matrix. We cluster
the cells to obtain tentative _cell type_ groupings and generate pseudobulk
counts by summing all the counts for each donor, cell type, and gene. This
creates a tensor of shape $D\times C\times G$, where D, C, and G represent the
number of donors, cell types, and genes, respectively. We then normalize the
pseudobulk tensor to CPM and apply tensor factorization methods with rank 8
and perform 1000 iterations for each method.

To evaluate the performance of the factorization methods in recovering the
eight true gene expression programs (GEPs), we computed the Pearson
correlation [2] between each of the eight latent factors in the gene mode
obtained via factorization and the original GEPs. This correlation was used to
establish a one-to-one alignment between the factors and the GEPs. We
calculated the average Pearson Correlation between each factor and its
corresponding GEP as the accuracy score of the method. The results of our
analysis are presented in Figure 3 for three different levels of simulated
intensity for the activity GEPs (mean log2 fold change of differentially
expressed genes (log2FC) $\in\\{0.25,0.5,0.75\\}$). We observed that when the
signal is strong enough (log2FC $=0.75$), CNMF and C-ZIPTF perform comparably.
However, when the signal intensities are lower (log2FC $\in\\{0.25,0.5\\}$),
C-ZIPTF clearly outperforms all other methods. Additionally, we found that
ZIPTF without consensus aggregation also performs better than the other
factorization methods, indicating that both the use of the ZIP model and the
consensus aggregation improve the accuracy of the method in recovering GEPs.

###  4.3. Real single-cell RNA-Seq data analysis

We applied C-ZIPTF to a real-world single-cell RNA sequencing dataset of
peripheral blood mononuclear cells (PBMCs) from patients with Lupus, reported
in [22]. We obtained the single-cell RNA-Seq data from GEO using accession
number GSE96583 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96583.
As described in [22], the dataset contains 29,065 cells from eight patients,
which are divided into stimulated and control groups, with the former being
treated with interferon beta (IFN-$\beta$), a cytokine that modulates the
transcriptional profile of immune cells. As part of the preprocessing step, we
filter out multiplets and cells without a cell type assignment. Additionally,
we remove samples and cell types that constitute less than 2 percent of cells.
After these filtering steps, the dataset contained 14 samples, 7 control and 7
stimulated, and 6 cell types: CD4 T-cells, CD14+ Monocytes, B-cells, CD8
T-cells, NK- cells, and FCGR3A+ Monocytes. In order to facilitate biological
interpretability of factors and reduce noise in the tensor formed we removed
genes that are either not provided with HGNC symbols [38], or had a total
count of less than 50 across all cells. Finally, we create a pseudobulk tensor
by summing up the raw counts for each cell type, sample, and gene. The
resulting pseudobulk data tensor has dimensions S $\times$ C $\times$ G (14
$\times$ 6 $\times$ 9,276), where S, C and G denote the number of samples,
cell types and genes respectively. We normalize the tensor such that each
sample-cell type pair has a total of $10^{6}$ counts. We first determined the
optimal rank for the data by running C-ZIPTF with a range of ranks from 2 to
14 and just 5 restarts. Some of the metrics we considered in deciding the
optimal rank including explained variance and silhouette score are shown in
Figure 4. We select rank 8, which exhibited a high explained variance of
0.969.

Figure 4. Metrics used in selecting the optimal rank for running C-ZIPTF on
real single-cell RNA sequencing dataset [22] of immune cells stimulated with
interferon beta (IFN-$\beta$).

As depicted in Figure 5, C-ZIPTF successfully identifies both cell type
identity and condition-specific gene expression programs (GEPs). Notably,
factor 4 represents an identity GEP that remains active in all B-cells,
irrespective of the condition. The genes exhibiting the highest loadings for
this factor are well-established B-cell markers, such as MS4A1, CD79A, and
BANK1 [12]. Furthermore, we performed gene set enrichment analysis [40] of
these factors using GSEApy [13] in Python. This analysis revealed enrichment
pathways consistent with B-cell characteristics, including B-cell activation
and the B-cell receptor signaling pathway.

Figure 5. Full set of factors recovered by running C-ZIPTF on real single-cell
RNA sequencing dataset [22] of immune cells stimulated with interferon beta
(IFN-$\beta$). Each row represents a factor, and the first three columns
display the three modes: sample, cell type, and gene. The $y-$axis in the
sample and cell type modes represent the loading of the sample or cell type on
that factor. The gene mode exhibits the top 20 genes associated with the
factor. The last column provides the top 3 enriched terms obtained from a gene
set enrichment analysis.

Conversely, factor 1 and factor 6 capture distinct gene expression programs
that are specifically activated in IFN-$\beta$ stimulated samples. Factor 1
captures a cross-cell-type response to IFN-$\beta$ stimulation, whereas factor
6 represents a Monocyte-specific response. These findings align with previous
studies that have reported a Monocyte-specific response to IFN-$\beta$
stimulation [30]. Furthermore, gene set enrichment analysis revealed
enrichment in pathways such as the cellular response to type I interferon and
inflammatory response, among others. For a comprehensive list of factors
identified by C-ZIPTF and their associated gene expression programs, please
refer to Figure 5.

##  5\. Conclusion

Zero-inflated count data is a common phenomenon in a wide range of fields,
including genomics, finance, risk management, healthcare, and social sciences.
However, traditional tensor factorization methods have limited effectiveness
when dealing with zero-inflated data, often yielding inaccurate and unstable
results across runs with different initializations. To overcome these
challenges, we propose ZIPTF, a Bayesian tensor factorization model that is
specifically tailored to zero-inflated count data. Additionally, we introduce
a generic meta-analysis framework for consensus-driven tensor factorization.
By combining these two approaches, we develop a novel method called C-ZIPTF
that achieves both high accuracy and stability, and outperforms state-of-the-
art baselines on synthetic and real data. Our proposed method provides a
useful tool for researchers in various fields to gain deeper insights into
their data.

The Bayesian approach for tensor factorizations offers several other
advantages over maximum likelihood estimation-based methods for tensor
factorization. These include the ability to incorporate prior knowledge,
perform model selection, and quantify uncertainty in the parameter estimates.
However, it is important to note that Bayesian methods can be computationally
expensive and require careful specification of prior distributions, which may
require expert knowledge. Moreover, the tensor methods discussed in this paper
rely on a multilinear factorization form and may be inadequate for capturing
more complex, nonlinear relations in the data. To overcome this limitation,
one possible solution is to integrate a kernelized approach into the
factorization. In future work, we plan to focus on the careful design of
kernel functions that would enable us to effectively capture nonlinear
patterns in the data.

##  6\. Supplementary Material

###  6.1. Implementation of ZIPTF and C-ZIPTF

We present a Python implementation of a versatile Bayesian Tensor
Factorization method using Variational Inference. Our implementation leverages
Pyro [3], a probabilistic programming framework built on PyTorch. The
BayesianCP class inherits from torch.nn.Module and offers functionalities for
model fitting and summarizing the posterior distribution of factor matrices.
During model fitting, Stochastic Variational Inference (SVI) is employed with
an Adam optimizer [23, 18]. The current implementation supports three models:
Zero Inflated Poisson model (ZIPTF), a Gamma Poisson model (GPTF) [37], and a
Truncated Gaussian model (TGTF) [17].

###  6.2. Implementation of baseline methods

As mentioned in Section 6.1, we utilize the same implementation for the other
Bayesian tensor factorization approaches (Gamma Poisson Bayesian Tensor
Factorization and Truncated Gaussian Bayesian Tensor Factorization) as the
ZIPTF method and the code is provided at https://github.com/klarman-cell-
observatory/scBTF and https://github.com/klarman-cell-
observatory/scbtf_experiments. For the remaining baselines used in our
comparisons we use the following implementations:

  * •

Non-negative Matrix Factorization (NMF): We use the Python implementation
provided in the scikit-learn package. https://scikit-
learn.org/stable/modules/generated/sklearn.decomposition.non_negative_factorization.html.

  * •

Consensus Non-negative Matrix Factorization (cNMF): We use the Python
implementation described in [25], and provided on GitHub
https://github.com/dylkot/cNMF/tree/master

  * •

Non-negative CP via Alternating-Least Squares (NNCP-ALS): We use the Python
implementation provided in the Tensorly package.
http://tensorly.org/stable/modules/generated/tensorly.decomposition.non_negative_parafac_hals.html

###  6.3. Simulation details

We use a Python adaptation of the Splatter [45] statistical framework given in
[25] to simulate single-cell RNA-Seq data. The core of the simulation is a
Gamma-Poisson distribution used to generate a cell-by-gene count matrix. While
the original Splatter framework supports the simulation of both expression
outlier genes and technical dropout (random knockout of counts), the Python
adaptation in [25] only keeps outlier expression simulation. Since our method
is specifically adapted to handle dropout noise in single-cell data, we add
back the modeling of dropout to the Python adaptation. Specifically, after
sampling counts from a Poisson distribution, we simulate dropout noise by
calculating the probability of a zero for each gene from its mean expression
and using that to randomly replace some of the simulated counts with zeros
employing a Bernoulli distribution as described in [45].

The distribution of expression values prior to incorporating differential
expression was determined based on parameters estimated from a random sample
of 8000 cells from an organoid dataset as described in [25]. Specifically, the
library size of a cell is sampled from a Lognormal distribution derived from a
Normal distribution with a mean of 7.64 and a standard deviation of 0.78. The
mean expression of a gene is sampled from a Gamma distribution with a mean of
7.68 and a shape of 0.34. With the probability of 0.00286, a gene will be an
outlier from this Gamma distribution and will instead be sampled from a
Lognormal distribution derived from a Normal distribution with a mean of 6.15
and standard deviation of 0.49. Additionally, we set a 5% doublet rate.
Doublets are formed by randomly sampling a pair of cells, combining their gene
counts, and downsampling such that the total count equals the larger of the
two.

## References

  * [1] E. Acar, T. Kolda, D. M. Dunlavy, All-at-once optimization for coupled matrix and tensor factorizations, arXiv:1105.3422 (2011).
  * [2] J. Benesty, J. Chen, Y. Huang, and I. Cohen, Pearson correlation coefficient, Noise reduction in speech processing (2009), pp. 37–40, Springer.
  * [3] E. Bingham, J.P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. Szerlip, P. Horsfall, and N.D. Goodman, N.D, Pyro: Deep Universal Probabilistic Programming, Journal of Machine Learning Research (2018), 19(108) , 1-6.
  * [4] C. Bishop, Pattern Recognition and Machine Learning, New York:Springer (2006).
  * [5] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, Variational Inference: A Review for Statisticians, Journal of the American Statistical Association (2017), 112:518, 859-877, DOI: 10.1080/01621459.2017.1285773
  * [6] M.M. Breunig, H.-P.Kriegel, R.T. Ng, and J. Sander, LOF: Identifying density-based local outliers, Proceedings of the 2000 ACM SIGMOD international conference on Management of data—SIGMOD ‘00, Dallas, TX, USA.
  * [7] J.P. Brunet, P. Tamayo, T.R. Golub, J.P. Mesirov, and E. S. Lander, Metagenes and Molecular Pattern Discovery Using Matrix Factorization, Proceedings of the National Academy of Sciences of the United States of America (2004), 101(12), 4164–4169, http://www.jstor.org/stable/3371580.
  * [8] A. Cemgil, Bayesian inference for nonnegative matrix factorisation models, Computational Intelligence and Neuroscience (2009).
  * [9] J.C. Chen, The nonnegative rank factorizations of nonnegative matrices, Linear Algebra Appl.(1984), vol. 62, pp. 207–217.
  * [10] E.C. Chi and T.G. Kolda, On tensors, sparsity, and nonnegative factorizations, SIAM Journal on Matrix Analysis and Applications (2012), 33 1272–1299.
  * [11] M. Chigona and C. Gaetan, Semiparametric zero-inflated Poisson models with applications to animal abundance studies, Environmetrics (2007), 18(3), 303-314.
  * [12] C. Domínguez Conde, C. Xu, L.B. Jarvis, and others, Cross-tissue immune cell analysis reveals tissue-specific features in humans, Science(2022), 376(6594), eabl5197.
  * [13] Z. Fang, X. Liu, and G. Peltz, GSEApy: a comprehensive package for performing gene set enrichment analysis in Python, Bioinformatics, Volume 39, Issue 1, January (2023), btac757, https://doi.org/10.1093/bioinformatics/btac757
  * [14] E. Frolovand I. Oseledets, Tensor methods and recommender systems, WIREs Data Mining Knowl Discov (2017), 7: e1201, https://doi.org/10.1002/widm.1201.
  * [15] S.K. Ghosh, P Mukhopadhyay, and J.C. Lu, Bayesian analysis of zero-inflated regression models, Journal of Statistical Planning and Inference (2006), 136(4), 1360-1375.
  * [16] P. Gopalan, J. Hofman, and D. Blei, Scalable recommendation with Poisson factorization, Proceedings of the Thirty-First Conference on Uncertainty in Artificial Intelligence (2015).
  * [17] J. L. Hinrich, K. H. Madsen, and M. Mørup, The probabilistic tensor decomposition toolbox, Machine Learning: Science and Technology (2020), available online at https://github.com/JesperLH/prob-tensor-toolbox.
  * [18] M.D. Hoffman, D.M. Blei, C. Wang, and J. Paisley,Stochastic variational inference(2013), Journal of Machine Learning Research.
  * [19] C. Hu, P. Rai, C. Chen, M. Harding, and L Carin, Scalable Bayesian Non-negative Tensor Factorization for Massive Count Data, Machine Learning and Knowledge Discovery in Databases (2015), ECML PKDD.
  * [20] Y. Ji, Q. Wang, X. Li, and J. Liu, A Survey on Tensor Techniques and Applications in Machine Learning, IEEE Access(2019), vol. 7, pp. 162950-162990, doi: 10.1109/ACCESS.2019.2949814.
  * [21] I. Jung, M. Kim, S. Rhee, S. Lim, and S. Kim, MONTI: A Multi-Omics Non-negative Tensor Decomposition Framework for Gene-Level Integrative Analysis, Frontiers in genetics (2021), 12, 682841. https://doi.org/10.3389/fgene.2021.682841.
  * [22] H.M Kang, M. Subramaniam, and others, Multiplexed droplet single-cell RNA-sequencing using natural genetic variation, Nature biotechnology (2018), 36(1), 89–94, https://doi.org/10.1038/nbt.4042.
  * [23] D.P. Kingma and J. Ba, Adam: A method for stochastic optimization (2014), arXiv preprint arXiv:1412.6980.
  * [24] T. G. Kolda and B. W. Bader, Tensor decompositions and applications, SIAM review 51 (2009), no. 3, 455–500.
  * [25] D. Kotliar, A. Veres, M. A. Nagy, S. Tabrizi, E. Hodis, D.A. Melton, and P.C. Sabeti, Identifying gene expression programs of cell-type identity and cellular activity with single-cell RNA-Seq, eLife (2019), 8, e43803, https://doi.org/10.7554/eLife.43803
  * [26] B. Kurt, Bayesian poison tensor factorization for dummies (2018), https://github.com/bariskurt/bptf.
  * [27] D. Lambert, Zero-Inflated Poisson Regression, with an Application to Defects in Manufacturing. Technometrics (1992), 34(1):1-14.
  * [28] D. Lähnemann, J. Köster, E. Szczurek, D. J. McCarthy, S. C. Hicks, M. D. Robinson, C. A. Vallejos, K. R. Campbell, and others, Eleven grand challenges in single-cell data science, Genome biology (2020), 21, 1–35.
  * [29] D. D. Lee and H. S. Seung, Learning the parts of objects by non-negative matrix factorization, Nature (1999), vol. 401, no. 6755, pp. 788–791.
  * [30] M. Lotfollahi, S. Rybakov, and others, Biologically informed deep learning to query gene programs in single-cell atlases, Nat Cell Biol (2023), 25, 337–350,https://doi.org/10.1038/s41556-022-01072-x
  * [31] J. MacQueen, Some methods for classification and analysis of multivariate observations (1967).
  * [32] A. Özdemir, M. A. Iwen, and S. Aviyente, A multiscale approach for tensor denoising, IEEE Statistical Signal Processing Workshop (2016), pages 1–5.
  * [33] J. Paisley, D. Blei, and M. Jordan, Bayesian nonnegative matrix factorization with stochastic variational inference, Handbook of Mixed Membership Models and Their Applications (2014). Chapman and Hall/CRC.
  * [34] R. Ranganath, S. Gerrish, and D. Blei, Black box variational inference (2014), Artificial intelligence and statistics, PMLR (pp. 814-822)
  * [35] H. Robbins and S. Monro,A stochastic approximation method(1951), The annals of mathematical statistics, 400-407.
  * [36] P.J. Rousseeuw, Silhouettes: a graphical aid to the interpretation and validation of cluster analysis, Journal of computational and applied mathematics (1987), 20, 53-65.
  * [37] A. Schein, J. Paisley, D. Blei, David, and H. Wallach, Bayesian Poisson Tensor Factorization for Inferring Multilateral Relations from Sparse Dyadic Event Counts, Proceedings of the 21st ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (2015), https://doi.org/10.1145/2783258.278341410.1145/2783258.2783414.
  * [38] R.L. Seal, B. Braschi, and others, Genenames.org: the HGNC resources in 2023, Nucleic Acids Res. PMID: 36243972, DOI:10.1093/nar/gkac888.
  * [39] M. Simchowitz, Zero-inflated Poisson factorization for recommendation systems, Junior Independent Work (advised by D. Blei), Princeton University, Department of Mathematics (2013).
  * [40] A. Subramanian, P. Mootha, and others, Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles, Proceedings of the National Academy of Sciences (2005), 102(43), 15545-15550.
  * [41] G. Tomasi and R. Bro, A comparison of algorithms for fitting the parafac model, Comput. Stat. Data Anal. (2006),50 1700–34.
  * [42] M. Wang, J. Fischer, and Y. S. Song, Three-way clustering of multi-tissue multi-individual gene expression data using semi non-negative tensor decomposition, The annals of applied statistics (2019), 13(2), 1103–1127, https://doi.org/10.1214/18-aoas1228
  * [43] H. Yang, G. He, and Y. Dong, Nonnegative Tensor Decomposition and It’s Applications in Image Processing, Proceedings of the International Conference on Computer Science and Information Technology, ICCSIT (2008), 212 - 217. 10.1109/ICCSIT.2008.99.
  * [44] Y.K. Yılmaz and A.T. Cemgil, Probabilistic Latent Tensor Factorization in Latent Variable Analysis and Signal Separation, Lecture Notes in Computer Science (2010), vol 6365. Springer, Berlin, Heidelberg, https://doi.org/10.1007/978-3-642-15995-4_43.
  * [45] L. Zappia, B. Phipson, and A. Oshlack, Splatter: simulation of single-cell RNA sequencing data, Genome Biol (2017),18 174, https://doi.org/10.1186/s13059-017-1305-0.
  * [46] M. Zhou, L. Hannah, D. Dunson, and L. Carin, Beta-negative binomial process and Poisson factor analysis. arXiv:1112.3605 (2011).