Population Priors for Matrix Factorization
Authors: Sohrab Salehi
Achille Nazaret
Sohrab P Shah
David Blei
Venue: Accepted by TMLR ()

Abstract: We develop an empirical Bayes prior for probabilistic matrix factorization.  
Matrix factorization models each cell of a matrix with two latent variables, one
associated with the cell's row and one associated with the cell's column.  
How to set the priors of these two latent variables?  
Drawing from empirical Bayes principles, we consider estimating the priors from data, to
find those that best match the populations of row and column latent vectors.  
Thus we develop the twin population prior.
We develop a variational inference algorithm to simultaneously learn the empirical priors
and approximate the corresponding posterior. We evaluate this approach with both
synthetic and real-world data on diverse applications: movie ratings, book ratings, single-cell gene expression data, and musical preferences. Without needing to tune Bayesian hyperparameters, we
find that the twin population prior leads to high-quality predictions, outperforming manually tuned priors.

Published in Transactions on Machine Learning Research (12/2024)
Population Priors for Matrix Factorization
Sohrab Salchi
edu
Computational Oncology, Department of Epidemiology and Biostatistics, Memorial Sloon Kettering Cancer Center, New York, NY, USA Irving Instttute for Cancer Dynomics, Columbia University, New York, NY, USA
Achille Nazaret
aon2108@cohumbie edu
Department Computer Science, Columbia University, New York NY, USA
Sohrab P Shah
shahs3 ory
Computational Oncology, Department of Epidemiology and Biostatistics, Memorial Sloan Kettering Cancer Center, New York, NY, USA
David M Blei
david. bleobblbbaall edu
Department of Stalistics, Columbia University, New York, NY, USA
Duta Science Institute, Colombia University, New York, NY, USA
Department of Computer Science, Cohumbia University, New York NY, USA
Reviewed on OpenReview: https: openreview. net/ formm AT9G5s1ppp
Abstract
We develop empirical Bayes prior for probabilistic matrix factorization. Matrix factorization models each cell of a matrix with two latent variables, one associated with the cell's row and one associated with the cell's column. How to set the priors of these two latent variables? Drawing from empirical Bayes principles, we consider estimating the priors from data, find those that best match the populations of row and column latent vectors. Thus we develop the twin population prior. We develop variational inference algorithm to simultaneously learn the empirical priors and approximate the corresponding posterior. We evaluate this approach with both synthetic and real-world data diverse applications: movic ratings, book ratings, single-cell gene expression data, and musical preferences. Without needing to tune Bayesian hyperparameters, we find that the twin population prior leads to high-quality predictions, outperforming mamually tuned priors.
1 Introduction
This paper is about empirical Bayes methods for sctting the priors in Bayesian matrix factorization (Mnih & Salakhutdinov, 2007; Gopalan et al., 2015). Matrix lactorization models each cell of a mairix with two latent variables, one associated with its row and one associated with its column. Matrix factorization has found broad applications across many fields, inchuding studying consumer behavior, understanding legislative patterns, asscssing pharmacoutical impacts, and exploring social networks (Gopalan et al., 2015; Koren et al., 2009; Gerrish & Blei, 2011; Jamali Ester, 2010).

---

Published in Transactions on Machine Learning Research (12/2024)
Suppose X1 is the observed entry in row and cell j, such user i's rating of movie j. As hierarchical model, a matrix factorization generates the data from the following process:
U² ~ (Ui))U))
(1)
Vj ~ Toal
(2)
X1, ~ P(Xes U, V)).
Here, U2 and V, are per row and per column specific latent vectors, and Trow and Tool are the priors distributions.
This formulation encompasssss many factorization models. In Gaussian matrix factorization (Mnih & Salakhutdinov, 2007), the priors are Gaussians and X1, is drawn from a Gaussian with mean U, Vj and variance o². In Poisson matrix factorization (Canny, 2004; Dunson & Herring, 2005; Gopalan et al., 2015), the priors are over the positive reals and X,, is drawn from a Poisson with rate U.! Vy.
An observed matrix of data X then defines posterior distribution P(U, V X) over the row variables and the column variables. The posterior can provide interpretations of the data and an avenue to form predictions about missing entries, for example, for a recommendation system.
The prior distributions on the row variables and column variables significantly impact the quality of the model posterior. How should we set them? Practitioners typically assume a simple parametric family for the priors, such as a Gaussian or a Gamma, and then find the prior hyperparameters best suited to the data, e.g.; with cross-validation (Salakhutdinov & Mnih, 2008; Schmidt et al., 2009). This approach can be effective, but it is expensive and only allows for priors from simple parametric family.
In this paper, we develop an empirical Bayes (EB) methodology for setting the priors (Robbins, 1992; Efron, 2012), learning them from data. The EB idea is to set the priors using the data, for example by finding the one that maximizes the marginal log- likeliiiooodd of the data. The EB idea is applicable to any prior distribution from which are repeatedly drawn multiple independent variables. For example, all the U, are independently drawn from the same Trow and the V are independently drawn from the same Tool- EB has found applications in applied sciences in diverse fields such as astronomy (Bovy et al., 2011), actuarial sciences (Bühlmann & Gisler, 2005), genomics (Smyth, 2005; Love al., 2014), economics (Frost & Savarino, 1986; Angrist et al., 2017), and survey sampling (Rao & Molina, 2015). EB priors have been successfully employed in simple hierarchical models, such as in variational autoencoders (Kingma & Welling, 2013; Tomezak & Welling, 2018; Kim & Mnih, 2018).
Matrix factorization, however, provides a different type of application of EB. In matrix factorization, there are two priors to set, one for the row variables and one for the column variables, and the same data contains information about them, namely the observed data X,, contains information about both U, and V,. Thus, we will find empirical Bayes priors for both row variables and column variables. The result is the twin empirical Bayes prior (TwinEB), a practical EB method for matrix factorization. Other methods for EB on matrix factorization include Wang & Stephens (2021) and da Silva et al. (2023); we discuss these in section 2.
Specifically, we model the two priors with mixture distributions, one for each prior. We use mixtures since they are flexible family of distributions that can approximate wide range of distributions (Titterington et al., 1985; Nguyen al., 2020). We use variational inference (Blei et al., 2017; Wainwright & Jordan, 2008) to simultaneously estimate the priors and approximate the corresponding posterior. We verify the efficiency and the robustness of this approach with real-world data about recommendation systems and computational biology, and for both Poisson matrix factorization and Gaussian matrix factorization.
We summarize the contributions as follows:
We develop the twin population prior, an EB prior for Bayesian matrix factorization.
2. We derive a variational inference algorithm to approximate the matrix factors' posteriors and learn the twin population prior simultaneously.
2

---

Published in Transactions on Machine Learning Research (12/2024)
3. We study the twin population prior both synthetic and real data and with two types of factoriza- tion. The automatically learned EB prior performs as well as the best prior chosen redrospectively.
2 Related work
One approach to setting the hyperparameters of priors involves using hicrarchical Bayesian models (Gelman, 2006). In this line of work, the prior's hyperparameters are treated as unknown and assigned prior with hyperpriors to be determined. Gopalan et al. (2015) employs this method for introducing hierarchical Poisson factorization, and Levitin et al. (2019) applies it to gene signature discovery from scRNAser data. Our work involves learnable priors (mixtures), which leads to a more flexible class models while keeping a simpler model as no extra variables for the hyperparameters are introduced.
A fruitful direction of research for setting priors for matrix factorization has also been the use of empirical Bayes, methodology by which one tries to find the prior that matches, in some sense, the population distribution of the data. In Rukat et al. (2017), the authors set single and fixed values for the prior hyperparameter basod on the expected value of the observed matrix. In Wang & Stephens (2021), the authors formulate prior elicitation for matrix factorization (MF) as steps in a variational expectation maximization algorithm, and in that context, find that it is equivalent to solving the EB normal means (EBNM) problem (Jiang& & Zhang, 2009). Our proposed approach is closely related to Wang & Stephens (2021), but we propose to update both priors simultancously which bypasses potentially costly numerical integration required in the EBNM step. A related line of research makes the connection to EBNM via random-matrix theory, focusing specifically on denoising principal components (Zhong et al., 2022). Closely related to TwinEB, da Silva et al. (2023) optimizes hyperparameters to match the prior predictive distribution to statistics computed from the data, and derives closed-form solutions for PMF and its hierarchical extensions. In Section 4, we compare our method to Wang & Stephens (2021) and da Silva et al. (2023). In Appendix E, we give a more detailed review these methods.
Wang et al. (2024) uses auxiliary information improve inference in matrix factorization. They assume the existence of an auxiliary matrix G, comprising per user extra features. Their method is closest to that of Wang & Stephens (2021), except that they define the per row latent variables a function of the auxiliary matrix G. They use an ensemble of regression trees to model the functional dependence of the row-wise latents on the auxiliary variables. They use a Gaussian likelihood, with a mean as a linear function of the user and item representations.
Van Linh et al. (2020) incorporate auxiliary information about the items (columns). They assume descrip- tions of each item exists. They use pre-computed word-embeddings to train a neural network to learn user reprosentations. They use Poisson likelihood with rate computed as the linear combination of the user and item representations.
Other flexible distribution learning methods can be used to specify expressive priors. For instance, Zhou et al. (2020) extend variational autoencoders (VAEs) and use normalizing flows (Papamakarios et al., 2021) to learn flexible priors on the latent embeddings of users and items. They combine user and item auxiliary information with these latent emboddings to learn user and item representation that are then passod through multi-layered perceptrons to generate the parameter of a Bernoulli distribution.
Empirical Bayes (EB) priors have been explored in the context of variational autoencoders (VAEs) (Kingma & Welling, 2013) under the name aggregated posterior, average encoding distribution (Hoffman & Johnson, 2016) or VampPrior (Tomczak & Welling, 2018). The VampPrior learns an amorlized posterior and a prior over the latent variables using a shared neural network. It models a prior on the rows only and was used to address postcrior collapse (Tomezak & Welling, 2018) or to learn disentangled representations (Kim & Mnih, 2018). Our method focuses on matrix factorization, and we derive EB priors for both row and column latent variables.
3

---

Published in Transactions on Machine Learning Research (12/2024)
3 Empirical Bayes priors for probabilistic matrix factorization
Our goal is to develop empirical Bayes (EB) priors for Bayesian matrix factorization models. We will focus here on Poisson matrix factorization (PMF). In Appendix B, we derive EB priors for Gaussian matrix factorization (GMF).
With matrix factorization, the presence of repeated and identically distributed latent variables for each row and each column provides the opportunily to learn their prior distribution from data. This is form of empirical Bayes (Robbins, 1992; Efron, 2012) that prescribes population prior (see Section 3.1).
This population prior aims to align the model's marginal distribution of observations with the observed population distribution. In the special case of mairix factorization, there are two distinct populations: the populations row vectors, and the population of column vectors. With TwinEB, we learn one prior for each population. This form of hicrarchical modeling without introducing an extra layer of latent variables.
Notations. We employ boldface symbols to distinguish vectors from scalars. We the notations:
N the number of rows (e.g,, individuals), D the number of columns (e.g,, features).
X e NNXD; the observed data matrix; X is the entry in row i and column j.
Prow (.) is the unknown distribution of rows vectors, and Pa (X.) the one of columns.
U E RL and Vj E R are model latent variables for row and column j: of dimension L.
Trow and Tool are the prior distribution of the latent variables.
3.1 Background: Population priors for simple hierarchical models
In this paper, we assume to have data X that is drawn from unkown distribution P*, also called population distribution, that would like to model with probabilistic model consisting of latent variables Z, prior over the latent variables n, and likelihood function P(X Z).
A crucial step in Bayesian statistics is the choice of the prior distribution; done arbitrarily, it can lead to suboptimal posterior inforence (Wang et al., 2021). We choose to follow an empirical Bayes principle that prescribes population prior (Hoffman & Johnson, 2016; Tomczak & Welling, 2018). This prior, by design, aligns the model's marginal distribution of observations with the population distribution P*(X).
We first focus family of latent variable models called simple hicrarchical models (Agrawal & Domke, 2021). The joint distribution factorizes as follows:
II (2.))((; 21);
(4)
where π is the prior distribution of z, To simplify notations, we then focus on the marginal likelihood of a single observation, æ, and its corresponding local latent variable z.
An empirical Bayes criterion that the marginal distribution of observations under the model, denoted P. (x), should match with their true population distribution p*(x) (Ignatiadis & Wager, 2022), that is:
(w) = Px(x)
= f (() z) dz.
Our goal is to set such that Equation (5) holds. The expression for the prior म that satisfies this conditions
is:
T(2) / Pz(z / px))d
(6)
Ep-(w) [Px x)),
4

---

Published in Transactions on Machine Learning Research (12/2024)
where (z(z x) is the (local) posterior distribution of the latent variable given the observation under the model. The definition presents two issues: the unknown true population distribution P*(x), and the fact that the target prior π is on both sides of Equation (6), explicitly on the left and implicitly via the posterior on the right. The research literature has approximated Equation (6) with Monte Carlo estimates of P*(x) and vonaational inference of Px(2 a)) (Hoffman & Johnson, 2016; Tomczak Welling, 2018).
3.2 Population priors for probabilistic matrix factorization
Our goal is to develop population priors for Bayesian matrix factorization models. The challenge is that unlike simple hierarchical models, there is no distinction between local and global latent variables, rather latent variables denote row- and column-specift random variables.
3.2.1 Twin population priors
We establish population priors for two latent variables, one for the latent variables of the rows Trow (U:) and one for the latent variables of the columns (V)). These priors will match two different populations, one of the row vectors and one of the column vectors:
Prow (X..) = Population distribution of row vectors
Poll (x.) Population distribution of column vectors
We recall that the populations distributions unkown distributions, from which we sample the rows and the columns.
We begin with population priors for row latent variables. As in section 3.1, we specify the prior based on an empirical Bayes principle such that the true marginal distribution the rows Pr (Xi.)) is aligned with the distribution ion of the rows under the model, that
f
Prow (X.:)= = Trow (U.) II P(Xis UU, Vj)dU,
(7)
for fixed set of column variables V. For Equation (7) to hold, a population prior should be used:
Trow (Ug V) =Ep, Prom (X
s..) [P ( (U, X,,,, V)],
(8)
where explicted the dependence of the EB prior the column latent variables V. Similarly for columns. for fixed set of row latent variables U, the empirical Bayes criterion is:
f Tool Vi)
Pio (X.s)
II P(X U. Vj)dV, 1=1
(9)
The prior that satisfies this criterion the column population prior:
Tool UV; Epp, (x I (V) Xx,,U)
(10)
Since there are two populations in need of prior specification, we call Equations (8 and 10) the twin population priors.
We have established the form of population priors in probabilistic matrix factorization. Next, we focus on how to estimate the twin population priors and how to approximate posterior inference under them.
In the remainder of the paper, we focus on the Poisson matrix factorization (PMF). We derive the priors for Gaussian matrix factorization (GMF) in Appendix B.
5

---

Published in Transactions on Machine Learning Research (12/2024)
3.3 Twin EB prior for Poisson matrix factorization
In Poisson matrix factorization (PMF), the row and column latent variables U, and V; are non-negative L-vectors and the likelihood in Equation (3) is Poisson:
Xis U4 V3 Poisson V3)
(11)
for iE E jM]. [D].
The log- logllikelihodd of the data is:
log P(X V. log
P(X U0 Vj)
i=1==1
(12)
N,D
Σ
log Poisson (Xi, V,
Some methods place Gamma priors on V, and U, (Gopalan et al., 2015). Note that this is Bayesian formulation non-negative matrix factorization (Cemgil, 2009).
To compute the population prior for the rows, Frow (U,; V) = Epp.. (X...) a (U; X V)], we face two
problems. Namely, we do not know the true population distribution of the rows Prom (X,,), and the popu- lation prior (U;) appears both sides of the equality.
To find the population prior, first notice that Monte Carlo estimate of Equation (8) writes as:
N Σ Pa (U, Xe.: V).
Trow (U:V)
(13)
When the prior satisfies Equation (13), this property is called self-consistency (Laird, 1978).
The structure of Equations (13) suggests to use families of mixtures of parametric distributions to approxi- mate the row and column population priors (Tomezak & Welling, 2018). Mixtures can approximate complex distributions when their number of components increases while having the convenience of remaining para- metric (Titterington et al., 1985; Nguyen et al., 2020). We choose to model the priors by dropping their dependence on the other variable and express them as,
Kr
Trow (U:) = Peron (U,) =
Σ
Wp P ox (Ui),
(14)
k=1
Teal V1,) Pocas (V1)
Σ pa PuPun.m (V₂);
(15)
where K" and K are the number of components in the mixtures, f,ow {u, w) and Pcot (v. p}. The locations j E TR Kg,, xL and v E Rxx the scales a E RR.XL and n E &.. and the mixture weights w K, and p E A K. are the parameters the mixtures priors. As K, and Ke increase, the priors become more and more expressive (see Figure 4). Figure shows a graphical model representation matrix factorization with EB priors.
In classical empirical Bayes setup, the idea is to set the priors that maximize the marginal likelihood of the data:
0 = arg max log P(X; Drow, Ocol);
(16)
Drow 0
where
log P(X; Prow, Ocal) log / Po.. (U) (V)P VUaVV.
(17)
In Section 3.4, we find grow, Bcol the same time as we approximate the model posterior P(U, V X; Kroww

---

Published in Transactions on Machine Learning Research (12/2024)
H, 0,
D
U
X
Vj
N
V/ Th P
Figure 1: Twin population priors for Poisson matrix factorization model. Shaded nodes observed while other nodes represent latent random variables. The empty squares indicate that we will fit these priors to the data.
3.4 Posterior inference in PMF with twin EB priors
Given data X, goal is calculate the posterior P(U, V X; 0row Geal), which also depends choice of priors frow, Gca The challenges that this posterior is intractable (for any prior) and simultancously want to fit the priors to satisfy the EB criterion in Equation (16).
Our strategy will be as follows. We will use variational inference (V1) (Blei et al., 2017) to approximate the posterior, taking gradient steps in the variational objective with respect to the posterior approximation (the variational family).
the same time, however, the variational objective of VT is an approximation (lower bound) of the Oorgiaaginal from Equation (16). So we also take gradient steps with respect the EB priors to maximize it. The result is algorithm that simultancously approximates the posterior and learns the EB prior.
The variational posterior. Consider a parameterized mean-field variational family,
(A V x) = II 2
9A5, (U.,)
II 1xg, (Vz)),
(18)
This family has parameters for each row's latent vector and each column's latent vector, X; and Xद respec- lively. We further define A = [} and A == [Ag.]- The full set of variational parameters is A = A.).
From the perspective of posterior inference, our goal is to set aA to minimize the KL divergence to the exact posterior:
À arg KL(gN; P(U, V X; Ocal)).
(19)
A
In detail, the variational family is bank of Log-Normals:
x = (ab.,, bb,1),
(20)
NS. (azu, by.));
(21)
9x_, (U,,) CN(('' bl.))
(22)
2x2, (Viu) = CN((3,, by.)-
(23)
Each Log-Normal is parameterized by its natural parameters and b:
(-o log(a)
(log x²²
(N(r;,,)) & exp
2h
2b

---

Published in Transactions on Machine Learning Research (12/2024)
To minimize the KL divergence in Equation (19), VI optimizes the variational parameters A to, equivalently, maximize the evidence lower bound (EIBO) (Blei et. al., 2017):
C(X; A, Brow, Gcal) (U.V
x) loo P(X |U, V)]
E Dog P(U; Frow, Oenl) + log P(V; Orow, Pco1]
-Bx log qa (U, V [)]. (24)
Here, we use gradient ascent to maximize L( A) with respect to A (Ranganath et al., 2014). We further use stochastic reparameterization gradients to take such steps (Kingma Welling, 2013; Rezende et al., 2014).
Maximum marginal likelihood. At the same time, we would like to set the prior parameters to maximize the marginal likelihood of the data (Equation (16)). The variational objective in Equation (24) conveniently also provides a lower-bound on the marginal likelihood (Blei et al., 2017):
log P( Brow, Gcll 2 0))
(25)
So, we will also follow stochastic gradients of the ELBO with respect to the prior parameters frow, Gcol to maximize L(X; 0, A) with respect to Orow, Ocat- This stralegy has been used the context linear regression (Mukherjee et al., 2023).
Twin EB. Putting these two pieces together, our algorithm is stochastic gradient ascent of the ELBO with respect to two sets of parameters. In optimizing with respect to A, minimize the KL divergence between aA and the posterior; in optimizing with respect to Frow, Ocol we maximize the (approximate) marginal likelihood of the data.
We use the Adam algorithm for stochastic optimization (Kingma Ba, 2014) with a batch size of 128, and using ten particles obtain unbiased noisy estimates of the gradient of the ELBO via the reparameterivation trick (the particles are samples from QA to estimate the expectation Eg of the ELBO with Monte-Carlo).
The details of the algorithm are in Algorithm 1. Our implementation is available at https://gi com/ blei-lab/TwinEB.
Algorithm 1 Variational inference for Poisson matrix factorization with twin EB priors
Input: Data X, number of particles S, learning rate Ç, number of iterations T, number of components Kr, Ke; number of latent dimensions L.
Output: Variational posterior parameters A", prior parameters frow Ocol *
Initialize: A (0)) Frow (0) 0coll (0).
for /= to T do
for 1 to N, 1 to L do
for = to do (s)
Sample c ~N(0..). 1).
(t-1) ((t-l) (3) + G..)). Compute (,(0) exp(a i,l end for end for for j = to 1- to L do
(
for = to do Sample (14) (0, 1).
Compute ((s)
j.l = exp(a (t-1) (1-1))
+ 1/1 ((.))
end for end for Estimate L(X; Brow'
1-1 Acol (1-1) A (1-1)) using Monto-Carlo in Equation (24) with samples U(4), V(a)
in place of Eas.
A (t) Frow (1) Ocol (t) Adam( (V ((Promm Q).
end for return AT), frow (T) (T), 0cool
8

---

Published in Transactions on Machine Learning Research (12/2024)
MovieLens-1M
Fixed Col. prior scale
Fixed Row prior scale
HOLL
Test HOHH
-.355
-1.40
Row prior scale
Column prior sca
Ru1322b
Fixed Col. prior scale
Fixed prior scale
Test HOLL
Test HOLL
10"
10"
10"
10"
Row prior scale
Column prior scale
UserArtists
Fixed Co,, prior scale =
Fixed Row prior scale =
Test- HOLL
Test HOLL
10"
Row prior scale
Column prior scale
GoodBooks
Fixed Col. prior scale
Fixed Row prior scale
Teat HOLL
Test HOLL
-1.35
-1.40
10"
10'
10"
Row prior scale
Column prior scale
TwinEB (GMF)
TwinEB Single (GMF)
.... GMF
Figure 2: Twin population priors induce robustness to prior selection The held-out likelihood is sensitive to the choice of the prior hyper-parameters. GMF endowed with population priors both row and column latent variables, TwinEB (GMF), achieves comparable or better results than other methods. Each sub-panel displays held-out log-likelihood from adjusting the column prior variance with a fixed row prior, while the right sub-panel does the opposite, varying the row prior variance with a constant column prior. We demonstrate four datasets, from the top to bottom, MovieLens 1M, Ru1322b-scRNAeee Aseq, UserArtists, and GoodBooks. In all datasets, we set L 15. Similar results hold for other values of (see Appendix D.3).
4 Experiments
We study Algorithm in several real-world matrix factorization settings: book ratings, movie ratings, artist preferences, and single-cell RNA sequence gene counts. For all datasets, we studied both Gaussian and Poisson matrix factorization. In all datasets, (see Appendix D.3). We found that TwinEB performs well or better than mamually searching for parameterized prior, and performs as well better than setting simple parameterized prior by empirical Bayes.
set L 15. Similar results hold for other values of L
4.1 Datasets
In this section we first review four real-world datasots, ranging from pser-preeerenees genomics, and then explore the impact of twin population priors on the performance of Gaussian and Poisson matrix factorization. In our experiments, fix the number of row and column mixtures to K, 70 and K, 100 respectively.

---

Published in Transactions on Machine Learning Research (12/2024)
MovicLens 1M. This dataset comprises 1 million ratings from 6,000 users (rows) on 4,000 movies (columns) (Harper & Konstan, 2015). The ratings are on a scale of to 5. Sparsity of this dataset, defined as the number of nonzero elements divided by the total number of entries, is 0.04. We can use matrix factorization to capture different aspects of user preferences and movie characteristics. Specifically, U2a may signify user i's affinity for aspect (e.g., genre); while Vj, may represent the degree to which movie j exhibits aspect 1.
Ru1322b. We analyze single cell gene expression data from patient with small cell lung cancer (Chan et al., 2021). This dataset comprises 4,000 highly variable genes (columns) across 5,308 cells (rows). For the GMF family, we applied two-step transformation. First, we performed a log transformation on the counts after adding a pseudo-count of one. Then, we standardized the non-zero elements. Sparsity of this dataset is 0.16. Each entry of the matrix denotes the mumber of transcripts of gene j in cell i. We explain the gene expression matrix via I. gene-modules, with U. as the activity of module I in cell i, and Via as gene j's contribution to module Matrix factorization can be used as exploratory data analysis (finding gene modules associated with malignancy in cancer) or as component in more complex analysis (causal inference (Wang& & Blei, 2019)). See Appendix C for details on preprocessing of the sequencing data.
User Artists We use the data introduced in Cantador et al. (2011), comprising 92,834 user-listened artist relations, across 1,892 users and 17,632 artists, with a maximum value of 352, 698. Similar to Moviel Lns 1M, we use matrix factorization to explain different aspects user preferences and artist groupings. Sparsity of this dataset 0.003. We note that this dataset was analyzed in da Silva et al. (2023).
GoodBooks Contains 6 million ratings across 51, 51,288 users and 10,000 books". Ratings range from to 5. Sparsity of this dataset is 0.01.
4.2 Evaluation Metric and Baselines
We evaluate model performance using the likelihood of unseen interactions, that is, test holdout-likellhood (HOLL). To build intuition, we will use the example of user-item interactions where rows indicate users and columns items, and the entrics of the input matrix their obscrved interactions. One approach to testing performance is to hold-out a fraction of the entries of the input matrix, and to train the model on the held-in portion, and test performance on the held-out entries. In this approach, all users are observed during training. To test generalizability to cases when new users are observed, we can further hold-out a portion of users (rows) during training, and then test the performance on these hold-out users. Specifically, we will test the performance on masked entries of these hold-out users. This procedure measures strong generalization (Steck, 2019).
At test time, we need learn user-specific latents for the held-out rows, For this, we run the model in a special training step where the column latents are held fixed. Since these held-out users are nevertheless coming from the same pool training step. This may result in a more accurate measurement of out-of-disttrbbttion generalizability by breaking spurious correlations between items and users due to biases present in the current dataset. These biases may include the time-period of data collection, platforms used that may not cxist in a new datasct. Finally, we report the model performance on 30% of the unseen entries of the held-out users. See the supplement for the derivation of Equation (45) and more details on the experiments.
training users, use only 30% of the observed entries for this special
Baseline. We evaluate the performance of the PMF (GMF) model with TwinEB against multiple baselines, namely (i) TwinEB-Single, (ii) PMF (GMF). TwinEB-Sinn Single is a simple form of TwinEB where all latent dimensions have an identical prior, which we learn. PMF is prior of the same family as the TwinEB, but with fixed hyperparameters. We compare TwinEB against large choice of fixed parameters, akin to hyperparameter selection.
1 Accessed http:///tthhab.cmm//yygnmnnn
10

---

Published in Transactions on Machine Learning Research (12/2024)
4.3 Experimental Procedure and Results: Gaussian Matrix Factorization
We preprocess the data as follows. We standardize each column by subtracting the mean from non-zero entries and dividing the result by their standard deviation. We study two scenarios: (i) maintaining fixed prior row-wise variables while varying the prior on column-wise variables, and (ii) holding the prior on column-wise variables constant and adjusting the prior row-wise variables. We set the variational family as well as the mixture components for the population priors to be Gaussian. We set the fixed prior to NV 1) and vary the variance of the non-fixed one over (0.001, 0.01, 0.1, 1.0, 10). For example, in Figure 2, the left column of plots corresponds to fixing the prior on the wiise variables to P(V,) = =0 1) and varying the prior on the row-wise variables as P(U,) = -N(0, o") with o² E {0.001, 0.01, 1,
We compared the GMF model to TwinEB (GMF) and TwinEBBSingle (GMF). We treat zero entries missing values.
Figure 2 displays the outcomes for four real-world datasets. In GoodBooks and MovieLens- 1M datasets, TwinEB achieves the highest test HOLL. For the UserArtists dataset, TwinEB does better than the fixed prior. Similar results were obtained by varying L to other values and are reported in Supplementary Figures 5, in Appendix D.
We also compare to the method of Wang & Stephens (2021), we used the corresponding package flashr, that currently supports GMF. We designed an imputation experiment where we held-out 10% of entries of a standardized matrix, and compare reconstruction accuracy on non-zero entries. We applied flashr to the above four real-world datasets. Except in Ru1322b, the optimization objective in flashr encounters NaN values, which halt execution and terminate without results. For the Ru1322b dataset it achieves a mean absolute error of 0.64 vs 0.51 for the TwinEB. For a comparison on simulated data, please see Supplementary Table S1. We note that it is not straightforward to compare our method to that of Zhong et al. (2022); the software implementation does not immediately support missing data and imputation.
4.4 Experimental Procedure and Results:
Poisson Matrix Factorization
For the MovieLens 1M and GoodBooks datasets, we binarize the matrix, setting entries to 'one' if a user has rated a movie and to 'zero' otherwise. For the Ru1322b dataset, we normalize the rows such that the sum of all rows are equal; we then round each value to the nearest integer. This is to account for the effect of library size (Heumnos et al., 2023). We treat zeros as missing data. In all models, we set the variational family to be Log-Normal, and the prior (or the mixture components for the population priors) to be Gamma. Similar to GMF, we vary the row or column prior parameters, while keeping the other one fixed. We parameterize the Gamma prior by its mcan and variance Gamma(μ, 02), where μu = 2/3 and 0² = a/8². In cach scenario, we set the fixed prior to Gamma(1, 10). For the varying prior, we set its mean to μ = and change its variance along {0.01, 0.1, 0.25, 1.0, 10.0, 100.03. Figure 3 shows the results.
In datasets except for GoodBooks, TwinEB outperforms TwinEBSingle. In the UserArtists dataset, TwinEB out performs other methods. TwinFB, on average, takes about 1.5X the runtime of PMF (ranging form 1.1X to 2.1X). We point out that finding a good prior using grid-search incurs the sum of the cost of evaluations of the individual points in the grid, over 6.5X the running time of the PMF with TwinEB prior. We note that using a grid to find fixed priors (or alternatively, learning the scalar prior) does not guarantee reaching a best test HOLL Indeed the UserArtists datasot, TwinEB yields the best test HOLL.
We also compared TwinEB to the method of da Silva et al. (2023); given the data, it estimates values for the shape and rate parameters of the Gamma prior for both the row and column r.v.s. Table displays the estimated hyperparameters and the resulting test HOLL.
PMF equipped with fixed priors values calculated from this methods yields Test HOLL that are on par or slightly worse than the TwinEB method. This method yielded negative hyperparameter estimates for the M and GoodBooks datascts. These parameters are effectively not usable, since the parameters of Gamma must be positive. Hence, we emphasize that unlike ours, this methods cannot be used on all datasets. We explore additional baselines in Appendix D.2.
11

---

Published in Transactions on Machine Learning Research (12/2024)
MovieLens-1M
Fixed Col. prior var. =
10
Fixed Row prior var. 10
Test: HOLL
Test HOLL
-30
-35
16'
Row prior variance
Column variance
Ru1322b
Fixed Col. prior
10
Fixed Row prior 10
Test HOLL
-9.00
Test HOLL
-9.25
10"
10'
Row prior variance
Calumn prior variance
UserArtists
Col. prior var. 10
Fixed Row prior var. 10
Test HOLL
Test HOLL
1600
10'
10"
Row prior variance
Column prior variance
GoodBooks
Fixed Col. prior var. = 10
Fixed Row prior var. = 10
Test- HOLL
Test HOLL
-.00
-15
10"
10
10"
Row prior variance
Column prior variance
TwinEB (PMF)
TwinEB-Single (PMF)
PMF
Figure 3: Twin population priors induce robustness prior selection. As in Figure 2, but for PMF.
Table 1: Test Heldout Loglikelihood from running PMF with prior hyperparameters sets using the method of da Silva et al. (2023). For MovieLens- 1M and GoodBooks datasets, da Silva et al. (2023) yielded inadmissible hyperparameters.
test HOLL
Dalasel
da Silva et al. (2023)
TwinEB
Ru1322b
-9.53
-9.46
Artistsss
-1,583
-1,481
MovieLens-IM
NA
-3.04
GoodBooks
NA
-3.95
4.5 Simulation: Complexity of the Prior
Here, we examine the performance of the population priors as the number of mixture components are varied. To this end, we simulate a 1, 000 by 500 dataset, with dimensional row and column-wise r.v.s. We sample the and colurn-wise from a mixture of 15 and 20 Gamma distributions respectively, the
12

---

Published in Transactions on Machine Learning Research (12/2024)
Simulated Data Study of K
2.180
-2.885
Test HOLL
-2.190
-2.000
components:
components:
-2.
#Column components:
12.4 15
32
40
50
64
70
80
100
120
#Row components
Figure 4: Increasing the number of mixtures improves the model performance. Each line shows the average lest-HOLL over 10 seeds for fixed value for Key the number of column mixtures, and varying values of K,, the number of row mixtures. Here, set L = 64 to its true simulated value. The results similar for L 32 (see supplement).
and shape parameters of which are sampled from a Gamma(1, 1). The mixture weights are sampled from Dirichlet(eo..... .., eo) where the concentration parameter is co 10.
Performance Increase with Number of Mixture Components. The more mixture components, the more flexible the EB prior is. If learning EB prior is beneficial, then we expect the performance to increase with the number of mixture components. We run PMF with population priors with varying number row and column mixtures, setting the dimension of the latent variable to its true value, and, to avoid model misspecification, set the mixture prior family Gamma. We then report the average of the test HOLL across ten different seeds. Figure 4, the log-likelihood of test held-out data increases with the number of mixture components in the row prior.
5 Discussion
We introduced the twin population priors for probabilistic matrix factorization. We derived a method to estimate the corresponding posterior using Monte Carlo and variational inference. On real-world data, this method finds a prior as good as the best parametric prior chosen relrospectively.
Our method uses the widely applicable black box variational inference (BBVI); the availability and ease of use implementation brings flexible priors in the context of matrix factorization
of automatic differentiation, to a wider audience. For example, our proposed method works on multiple data types, both Gaussian and Poisson likelihoods, and it is easy to extend it to other reparameterizable likelihoods. In contrast, the methods of Wang et al. (2024) and da Silva et al. (2023) are tailored for only Gaussian and only count data only.
With TwinEB, the latent dimensions of the row-and column-wise variables is hyperparameter. In practice, we found that setting L at about 15 was generally good enough (used across all experiments) and ablation studies showed that the results were stable with the choice of L (GFM: Supplementary Figure 5 PMF: Supplementary Figure 6). However, an attractive property of the method of da Silva et al is that it can automatically set this value in closed form. Methods based on prior predictive statistics suggested by da Silva et al. (2023) could optimize the selection of L.
A limitation of our model compared others, is that it does not leverage auxiliary information. Auxiliary information information in addition the observed interaction matrix. For instance, for movie ratings, the additional information could be user-specific attributes including age, race, and occupation, while item-
13

---

Published in Transactions on Machine Learning Research (12/2024)
specific attributes could include description of the items. Methods like Wang al. (2024) and Van Linh et al. (2020) propose ways to integrate auxiliary information to improve the inference in matrix factorization.
We use a linear function to model the dependence of the observed interactions on the per-row per-column latent variables. This has the advantage of being interpretable and prevent overfitting (c.g., in the analysis of single cell gene expression data (Svensson et al., 2020)). However, the presence of non-linear interactions, this model would be misspecified. In that case, the link function could be changed, e-g., Zhou et al. (2020) use multi-layered perceptron (MLP) to model the dependence of the observed interactions on the per-row per-column latent variables.
Other flexible distribution learning methods can be used to specify expressive priors, such as normalizing flows used in Zhou et al. (2020), and one of our baseline methods, TwinNF (Appendix D.2). The model begins with a simple base distribution and iteratively transforms it into a complex distribution that can capture more intricate patterns the data.
One area of further work is extend this algorithm to lensor factorization (Kolda & Bader, 2009; Schein et al., 2015). While in matrix factorization, each entry of the observed matrix is explained via two latent variables, tensor factorization models will involve more. One detail to address is how to formally define the population distribution associated with each latent variable.
Investigating theoretical properties of the twin population prior a la da Silva et al. (2023) would offer insights into the behavior the priors before data observed that can a possible future research direction.
Broader Impact Statement
As a data-driven approach, our method is subject to potential pitsalls inherent in such methods, including tendency to inherit and amplify biases present in the data. If employed without proper supervision, this can lead to serious consequences including biased and/or discriminatory outcomes in real-world applications. These biases implicitly discriminate based on protected attributes (o.g., race, gender). One line of rescarch assumes specific protected features are known, and enforces fairness by decoupling these features from learnt latent variables (Togassh & Abe, 2022). If fairness criteria can be formulated in terms of regularization term, it can be reinterpreted as learning fairncss-aware priors and incorporated in our framework (Zhu et al., 2018).
6 Acknowledgements
This rescarch was funded in part through the NTH/NCI Cancer Center Support Grant P30 CA008748 (list all MSK authors); NIH grant 5K99CA277562-02 (S.S.), the Eric and Wendy Schmidt Center at the Broad Institute of MTT and Harvard (A.N.), the A ffiick Family Fund (A.N.), the National Science Foundation (NSF) grants 11S-2127869 and DMS-2311108 (D.B.), the Office of Naval Research grant N000142412244 (D.B). and the Simons Foundation (D.B.). SPS was partially supported by the Halvorsen Center for Computational Oncology and the MacMillan Center for the Non-Coding Cancer Genome.
References
Abhinav Agrawal and Justin Domke. Amortized Variational Inference for Simple Hicrarchical Models. NewrIPS, 34:21388 21399, 2021.
Joshua D Angrist, Peter D Hull, Parag A Pathak, and Christopher R Walters. Leveraging lotteries for school value-added: Testing and estimation. QJE, 132(2):871 919, 2017.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: A review for statisticians. JASA, 112(518):859-885 877, 2017.
Jo Bovy, David W. Hogg, and Sam T. Roweis. Extreme deconvolution: Inferring complete distribution functions from noisy, heterogeneous and incomplete observations. AOAS, 5(2B):1657 1677, 2011.
14

---

Published in Transactions on Machine Learning Research (12/2024)
Hans Btihlmann and Alois Gisler. A course in credibility theory and its applications, volume 317. Springer, 2005.
John Canny. GaP: a factor model for discrete data. SIGIR, pp. 122 222-129, 2004.
Iván Cantador, Peter Brusilovsky, and Tsvi Kullik. Second workshop information heterogeneity and fusion in recommender systems (hetrec2011). In ACM RecSys, pp. 387-388, 2011.
Ali Taylan Cemgil. Bayesian inference for nonnegative matrix factorisation models. COMPUT INTEL NEUROSC, 2009.
Joseph M Chan, Álvaro Quintanal-Villalonaa, Vianne Ran Gao, Yubin Xie, Viola Allaj, Ojasvi Chaudhary, Ignas Masilionis, Jacklynn Egger, Andrew Chow, Thomas Walle, et al. Signatures of plasticity, metastasis, and Immunosuppression in an atlas of human small cell lung cancer. Cancer Cell, 39(1)):447--4496 2021.
Eliezer de Souza da Silva, Tomasz KuL, Marcelo Hartmann, Arto Klami, et al. Prior specification for Bayesian matrix factorization via prior predictive matching. JMLR, 24(67):1-51, 2023.
David B Dunson and Amy H Herring. Bayesian latent variable models for mixed discrete outcomes. Bio- statistics, 6(1):11-25, 2005.
Bradley Efron. Large-Scale Inference: Empirical Bayes Methods for Estimation, Testing, and Prediction, volume Cambridge University Press, 2012.
Peter A Frost and James E Savarino. An empirical Bayes approach to efficient portfolio selection. JFQA, 21(3):293-305, 1986.
Andrew Gelman. Prior distributions for variance paramcters in hicrarchical models (comment on article by Browne and Draper). 2006.
Sean M Gerrish and David M Blei. Predicting legislative roll calls from text. In JCML, 2011.
Prem Gopalan, Jake M Hofman, and David M Blei. Scalable recommendation with hierarchical Poisson factorization. In UAI, pp. 226-335, 2015.
F Maxwell Harper and Joseph Konstan. The movielens datasets: History and context. THS, 5(4):1-19, 2015.
Lukas Heumos, Anna C. Schaar, et al. Best practices for single-cell analysis across modalities. Nat. Ren. Genet,, 24(8):550 572, 2023. ISSN 1471-0064.
Matthew D Hoffman and Matthew J Johnson. ELBO surgery: yet another way to carve up the variational evidence lower bound. 1(2), 2016.
Nikolaos Ignatiadis and Stefan Wager. Confidence intorvals for nonparametric empirical Bayes analysis. JASA, 117(539) 1166, 2022.
Tommi S Jaakkola and Michael Jordan. Bayesian parameter estimation via variational methods. Statistics and Compuling, 10:25-30 33, 2000.
Mohsen Jamali and Martin Ester. A matrix factorization technique with trust propagation for recommen-
dation in social networks. In RecSys, pp. 135 142, 2010.
Wenhua Jiang and Cun-Hui Zhang, General maximum likelihood empirical Bayes estimation of normal means. Ann. Stat., 2009.
Hyunjik Kim and Andriy Mnih. Disentangling by factorising. In ICML, pp. PMLR, 2018.
Diederik P Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. arXim: 12..:.12266.0000 12.0980, 2014.
Diederik P Kingma and Max Welling. Auto-Encoding variational Bayes. 1rXio:1312.6114, 2013.
15

---

Published in Transactions on Machine Learning Research (12/2024)
Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3):455-500 500, 2009.
Yehuda Koren, Robert Bell, and Chris Volinsky. Matrix factorization techniques for recommender systems. Computer, 42(8):30-337 37, 2009.
Nan Laird. Nonparametric maximum likelihood estimation of mixing distribution. JASA, 73(364):805 811, 1978.
Hanna Mendes Levitin, Jinzhou Yuan, Yim Ling Cheng, Francisco JR Ruiz, Erin C Bush, Jeftrey Bruce, Peter Canoll, Antonio lavarone, Anna Lasorella, David M Blei, al. De novo gene signature identification from single-cell RNA-seq with hierarchical Poisson factorization. Molecular systems biology, 15(2):08557, 2019.
Michael Love, Wolfgang Huber, and Simon Anders. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome biology, 15(12):1-21 2014.
Andriy Mnih and Russ Salakhutdinov. Probabilistic matrix factorization. NeurIPS, 20, 2007.
Sumit Mukherjoe, Bodhisattva Sen, and Subhabrata Sen. mean field approach to empirical Bayes estima- tion in high-dimensional linear regression. arX:::2309. 16843, 2023.
T Tin Nguyen, Hien Nguyen, Faicel Chamroukhi, and Geoffrey J McLachlan. Approximation by finite mixtures of contimuous density functions that vanish at infinity. Cogent Mathematics & Statistics, 7(1): 1750861, 2020.
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshmi- narayanan. Normalizing flows for probabilistic modeling and inference. JMLR, 22(57):1-64, 64, 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. NewrIPS, 32, 2019.
Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. In AISTATS, pp. 814 822. PMLR, 2014.
John NK Rao and Isabel Molina. Small area estimation. John Wiley & Sons, 2015.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate inference in deep generative models. In 1CML, pp. 1278 1286. PMLR, 2014.
Herbert E Robbins. An empirioal Bayes approach slatislics. Springer, 1992.
Tammo Rukat, Chris C Holmes, Michalis K Titsias, and Christopher Yau. Bayesian boolean matrix factori- sation. In ICML, pp. 2969-2978. PMLR, 2017.
Ruslan Salakhutdinov and Andriy Mnih. Bayesian probabilistic matrix factorization using Markov chain Monte Carlo. Tn ICML, pp. 880 880-887, 2008.
Aaron Schein, John Paisley, David M Blei, and Hanna Wallach. Bayesian poisson lensor factorization for inferring multilateral relations from sparse dyadic event counts. In SIGKDD, pp. 1054,,054, 2015.
Mikkel N Schmidt, Ole Winther, and Lars Kai Hansen. Bayesian non-negative matrix factorization. In ICA 2009, pp. 540- Springer, 2009.
Gordon K Smyth. Limma: linear models for microarray data. In Bioinformalies and compulational biology solutions using R and Bioconductor, pp. 397-420. Springer, 2005.
Harald Steck. Embarrassingly shallow autoencoders for sparse data. In WWW, pp. 3251 1751-227,, 2019.
Valentine Svensson, Adam Gayoso, Nir Yosef, and Lior Pachter. Interpretable factor models of singlo-cell RNA-seq via variational autoencoders. Bioinformatics, 36(11):3411 3444181-421,, 2020.
16

---

Published in Transactions on Machine Learning Research (12/2024)
David Michael Titterington, Adrian FM Smith, and Udi E Makov. Statistical analysis of finite mixture distributions. 1985.
Riku Togashi and Kenshi Abe. Fair matrix factorisation for large-scale recommender systems. arXiv, 2022. Jakub Tomezak and Max Welling, VAE with a VampPrior. In AISTATS, pp. 1214 1223. PMLR, 2018.
Ngo Van Linh, Due Anh Nguyen, Thai Binh Nguyen, and Khoat Than. Neural poisson factorization. IEEE Access, 8:106395-106407, 2020.
Martin J Wainwright and Michael I Jordan. Gruphical models, exponential families, and varialional inference. Now Publishers, Inc., 2008.
Wei Wang. Applications of adaptive Shrinkage in Multiple Statistical Problems. The University of Chicago, 2017.
Wei Wang and Matthew Stephens. Empirical Bayes matrix factorization. JMLR, (2(1):5332-5377 2021.
Yixin Wang and David M. Blei. The Blessings of multiple causes. JASA, 114(552):15744-58888 1596, 2019. Publisher: Taylor Francis.
Yixin Wang, David Blei, and John P Cunningham. Posterior collapse and latent variable non-identifaablity, NewrIPS, 34:5443 5455, 2021.
Zhiwei Wang, Fa. Zhang, Cong Zheng, Xianghong Hu, Mingxuan Cai, and Can Yang. MFAI: A Scalable Bayesian Matrix Factorization Approach to Leveraging Auxiliary Information. JCGS, pp. 11, 2024.
F Alexander Wolf, Philipp Angerer, and Fabian J Theis. SCANPY: large-scale single-cell gene expression data analysis. Genome biology, 19:1-5, 2018
Xinyi Zhong, Chang Su, and Zhou Fan. Empirical Bayes PCA in high dimensions. J. R. Stat. Soc., 84(3): 853-874 878, 2022.
Fan Zhou, Yuhua Mo, Goce Trajcevski, Kunpeng Zhang, Jin Wu, and Ting Zhong. Recommendation via collaborative autoregressive flows. Neural Networks, 126:52-66 64, 2020.
Ziwei Zhu, Xia Hu, and James Caverlee. Fairness-awai aware tensor- based recommendation. In ACM CIKM, pp. 1153-116,, 2018.
17

---

Published in Transactions on Machine Learning Research (12/2024)
A Introduction
This is the supplement for the manuscript titled PPopulation priors for matrix factorization'. The figures in this document supplement Figures 2, 3, and in the main text. We provide derivations for twin population priors for Gaussian matrix factorization, as well as computing the likelihood of the held-out data in section B. We then give more details about our experimental sctup, including the parameters used in training (c.g.. batch-size) in section C. Finally, we provide results for additional experiments in section D: for additional values of the latent dimension I. for the datasets studied in the main text.
Note that this document is accompanied by an archive file code. code.zip, that contains the source code, in- structions to install and run the code, and scripts to recreate the experiments and plot the figures in the manuscript. All scripts have been de-identified.
B Derivations
In this section, we give more details and derivations for the quantities defined in the main text. Specifically, section B.1 we derive the twin population priors for the Gaussian matrix factorization (GMF). In section B.2, derive the variational approximation to the twin population priors in GMF. Finally, in section B.3, we derive the expression for held-out likelihood in matrix factorization.
B.1 Twin population priors for Gaussian matrix factorization
In this section we derive population priors for Gaussian matrix factorization (GMF). In the classical GMF, the likelihood and priors on the latent variables are Gaussian (Mnih Salakhutdinov, 2007). The generative model as follows:
W₁, ~N(,0 oW),
1,1.=11.... I...,, Ze,1 NN(0, 02), 1-I...,,=11....
(26)
X1 U, ~ NO Z,,W,,,')), i=
where [Z,,] NXL and [W RDXL, with U, and V, and are L-vectors representing the row and column wise latent variables, and 0, Fw and oz are constant.
For TwinEB, our goal is to learn the prior for the row and column latent variables, so, in the main text, we construct mixture priors (Equations 14 and 15), where Pur.ox and Pe are Gaussian parameterized by their mean and scale parameters. That is,
Po-(U,) Σ
((U
(27)
km
Kr
PV.(V)) Σ
PvN(V;:Ps,n) mc),
(28)
(29)
where K" and K are the number of components in the mixtures, 0r = {u. a.w) and 0r (v, n, p). The locations E IR x. and e K.X. the scales e k, and m7 I2 K.XL and the mixture weights w E A Kr and p E A Ke are the parameters of the mixtures priors.
B.2 Posterior inference in GMF with twin population priors
Our goal to calculate the posterior P(U, V X.0",0" 8)), given data X, which also depends on our choice of priors 0r, 0 Our strategy is the same in the main text, namely, to simultaneously optimize for the
18

---

Published in Transactions on Machine Learning Research (12/2024)
parameters of the variational distribution and the TwinEB prior. We substitute the variational family in Equations (19) and (19) with Gaussian, as follows:
x = (4h,1- Hy,,),
(30)
Aj = (asa, bic);
(31)
1 (Ui,) Madaa 1/11;;
(32)
92 (V3,) = VV(a,, b51).
(33)
Each Gaussian is parameterized by its natural parameters and b:
N( a, b) x exp (az + bx²)
The log- ligllikeihhod of the data is:
N
D
log P(X Z, W)
log
(X
Ze, W,)
i=1==1
D
log.V
Xijs Σ 1=1
(34)
j=l
D
1
Xis -is-Ete
W
27
logov22 +
B.3 Held-out likelihood for matrix factorization
We use the log likelihood of held-out data as the score for each model. Let X = {Xir} denote Nout entries of the held-out rows that were maskod. We compute the likelihood of the masked entries via the posterior predictive distribution:
Nout
log P(Xout
X = log
II
P( Xin)
(35)
1=1
Nour log P(Xout X
(36)
For Bayesian matrix factorization we expand the summand in Equation (36) as follows:
// P(X.ut
log P(X:ut
X zs log
(UV)P(U, V X dU
(37)
//
se log
P(Xout
(.V)A V dd
(38)
log
E. U.V-gg(())
[P(xout V))
(39)
M
1
log
P(X."
(M)) (()),
(40)
M
m
where M is the number Monte Carlo samples from 4^ (U, V xin). In Equation (38) we approximate the true posterior P(U, V x with its variational counterpart gx. The log-likelihood score for the entire held-out data is then:
N
Nous
M Σ TTS
Σ
log P(Xout xia) =
Σ
log
P(Xout (T())) v(M)).
(41)
M
2
19

---

Published in Transactions on Machine Learning Research (12/2024)
C Experimental details
In this section we give more details on our experimental studies in the main manuscript. In section C.1 we describe preprocessing for the gene expression in the Ru1322-scRNAsse Aseq dataset. In section C.2, we specify the parameters used during training. In section C.4 we give a brief description the artifacts that accompany this supplementary material.
C.1 Preprocessing of the Ru1322-scRNAseg gene expression dataset
We use CollRanger version 6.0. to process the FASTQ and generate the unique molecular identifier (UMI) count matrices. We use scanpy to preprocess the data, and the seuratv3 algorithm select highly variable genes Wolf al. (2018).
C.2 Training details
We set a batch size of 128 in all our experiments. We ran Poisson and Gaussian matrix factorization experiments for maximum of 20, 000 iterations. By this step, all runs had convergod.
We initialized the learning rate for the and column variables, rlr and clr separately. We fix the initial 10 Monte Carlo
learning rate rlr E {0.01} and clr E {0.01}. In the experiments in the main text, samples to approximate the ELBO, while in the supplemental experiments, we use a single particle.
For the PMF experiments in the supplement, we subsample zeros as is standard (Gopalan et al., 2015). We uniformly randomly subsample the same number of zeros as non-zero values to estimate the likelihood. Concretely, let ((xout) denote the tog-likelihood of the masked entries of the held-out rows, Xout (Xtr). Then L(Xout can be decomposed the sum of non-zero and zero Xi.
L(Xout) L(Xout 20) + C(Xout
=0)=
(42)
where #2 and Xout = have cardinalities Nout nou-zero and Nout respectively.
We approximate Equation (42) by subsampling Nuub min(Nou zeros non-zoro) of the zeros. Let X Xeero denote multiset of zeros of cardinality Naub, then:
Nout Nout
Nout zero
((xout) x L(Xout
#0) +
Naub
X
Naub"
(43)
We ran experiment machine equipped with an NVIDIA A100 GPU with 80GB memory. We imple- mented all methods in pytorch (Paszke et al., 2019).
C.3
Evaluation
We evaluate each model based the likelihood of held-out data.
For given row-wise data split into held-in and held-out rows, let Xout = (Xout) denote the set of Nout masked entries of the held-out rows. We estimate the held-out log- likeliiiooddd as follows:
Nout
log P(X.uut
x
(44)
Nout
i=1
M Σ m=1
log
P(Xout
(mm) f(m)),
(45)
Nout
M
where M is the number of Monte Carlo samples from aA (U, V).
We randomly assign 20% of the rows as the test set, and the rest as training data. We then mask 20% of the entries random, and train the model this train set using ten random restarts. We use the 20% masked
20

---

Published in Transactions on Machine Learning Research (12/2024)
entries validation set. At test time, put aside 30% of the entries of the test at random these entries constitute the test set then we train the model on 40% of the rest of the entries. This procedure measures strong generalization (Steck, 2019).
To compute the held-out likelihood, after training the model on the held-in data, we fix its column parameters and then learn row parameters of the held-out rows. We then report the likelihood of the masked entries in the test set. For each model, we report the test HOLL of the random restart that achieved the best validation HOLL. In our experiments, we set M 500.
C.4 Implementation
Please refer to the code directory for the source code and instructions on how to run the model. The file README. md contains instructions for installing the software, and rumning it. Under the data directory, we included preprocessed data for the MovieLens-100f datasel, a smaller version of the MovieLens- 1M studied in the main text. In the notebooks directory, we put notebooks used for preprocessing the data, and plotting the figures in the manuscript. Finally, in the pipelines directory, we put nextflow scripts that recreate experiments that have run for the manuscript.
D Additional experiments
In this section we present additional experimental results. We show results for additional values of the latent dimension L. Concretely, section D.3 studies the elfects of the twin population priors Poisson and Gaussian matrix factorization, on three real world datasets, while section D.4 examines the sensitivity of the twin population priors to the choice of its hyper-parameters. We find that these results corroborate those that were presented in the main manuscript, that is, matrix factorization with traditional priors is sensitive to the choice of the hyper-parameters of the prior, and twin population priors is a robust. way set the prior in this family of models.
D.1 Simulated studies
We compare our method to that of Wang & Stephens (2021) on synthetic dataset, namely that used in Figure 4 in the main manuscript. Table S1 shows the mean absolute error of non-zero entries, averaged over 10 different repeats. The results are comparable, where TwinEB does slightly better.
Table S1: Comparison of methods based on Mean Absolute Error.
Method
Mean Absolute Error
Wang & Stephens (2021)
0.37
TwinEB (GMF)
0.35
D.2
Additional baselines
Here, we introduce two additional baseline methods. One based on normalizing flows (Papamakarios al., Gamma priors for our PMF experiments
2021), and the other is based on learning the hyper-parameters (close to Hierarchical Poisson Factorization of Gopalan et al. (2015)). We call both Twin, since are simultaneously learning the hyperparameters of the priors the row and column latents.
TWinHP In TwinHP, we infer the hyper-parameters of the two Gamma priors. This represents more exhaustive search in the hyper-parameter space of a model with uni-modal priors. Table $2 shows the results (a positive value indicates that TwinEB performs better). We note that in the Ru1322b dataset this method has the worst performance (as measured by Test Hold-Out likelihood), suggesting that the uni-modal prior is not well specified for that dataset. On the other hand, achieves the best HOLL in the MovieLens- 1M dataset, suggesting that uni-modal prior is more appropriate in this dataset.
21

---

Published in Transactions on Machine Learning Research (12/2024)
Table S2:
Comparison of TwinEB to TwinHP Test Heldout Loglikelihood (Delta-Tes) Test HOLL) across
various real-world datasets. Delta-Test HOLL Test HOLL Test HOLLTwinHP.
Datasct
Delta-Test HOLL
Ru1322b
0.54
UserArtish
211.74
MovicLens-1M.
-0.61
GoodBooks
0.01
TwinNF We implement Normalizing Flow (NF) (Papamakarios et al., 2021) prior by stacking 10 Planar flows (similar the method of Zhou al. (2020)) define a flexible prior distribution. Each Planar flow applies an affine transformation to the latent variables through the mapping:
f. = z u tanh(w"z + b),
where w are learnable parameters and b is a bias term. The determinant of the Jacobian is computed as:
det(Vf,) = + u² ((1 tanh"(w" + 6]) ww]
We do not claim that this implementation cannot be improved. Nevertheless, we find that TwinEB compares well with this NF prior, out-performing it in 3 out of the 4 real-world datascts. We show the result in the following table, where the second column shows the difference between Test-HOLL TwinEB and TwinNF (a positive value indicates that TWinEB performs better). We have used the same training procedure in the main text, with 10 seeds.
Table S3: As in Table S2, but for TwinNF, Delta-Test HOLL Test HOLLTwinEB - Test HOLL/WiONF.
Dataset
Delta-Test Test HOLL
Ru1322b
-0.12
Artistsssts
350.46
MovieLens- Le
0.18
GoodBooks
0.40
D.3
Additional values of latent dimensions
In the main text, we study the effects of twin population priors on the Poisson and Gaussian matrix factor- ization over four real datasets, with the dimension of the latent variables set to = 15. Here, we present results for additional values of L.
D.4 Performance Increases, then Plateaus with the Number of Mixture Components,
We show the performance of TwinEB on simulated data as we vary the number of row and column mixture components. Here, we add results for additional values of L,, and also study the MovieLens 100K dataset.. The result in Figures 7 suggests that as the number of mixture components increases, the performance of TwinEB improves. The gain is apparent when the number of row components K, increases. This is expected, our evaluation procedure, measures the generalizability for held-out rows, but not held-out columns.
In Figure we show the performance on real-world datasets as we vary the number of row and column mixtures. Again we find that, overall, as the performance first increases, until plateauing.
E Extended literature review
Here we review, in detail, the methods of Wang & Stephens (2021) and da Silva al. (2023). These methods closest to ours in hothodology.
22

---

Published in Transactions on Machine Learning Research (12/2024)
MovieLens-1Mr row prior scale 1.0
MovieLens- 1M: column prior scale 1.0
He dout log likellhood
He dout log likellhood
Latent dim.
-1.35
20
-1.35
30
40
-1.40
64
-1.40
Ru1322b: row prior scale = 1.0
Ru1322b: column prior scale 1.0
Haldout log likallhood
Heldout log likellhood
-1.2
-1.2
1.3
1.3
-1.4
-1.4
UserArtists: row prior scale 1.0
UserArtists: column prior scale 1.0
Heldout log likelihood
Heldout log likellhood
-140000
-140000
-160000
-160000
GoodBooks: row prior scale 1.0
GoodBooks: column prior scale 1.0
Heldout log likellhood
Heldout log likellhood
-1.4
-1.4
-1.5
-1.5
2
8
10
2
8
10
row prior scale
column prior scale
TwinEB (GMF)
GMF
Figure As in Figure 2, but with additional values latent dimension L.
E.1 Method of Wang & Stephens (2021)
Given the following model for observations X [Xt,]] E RND we can find the row and column factors by solving empirical Bayes normal means (EBNM) problem:
X Σ
Warai + Easj
1=1
U1 Jul Gu
(46)
Uj goy Go
Eig ~ NO0 1/Tis).
We set Ti,j = 1 for all i to simplify notation.
To solve for u and v, empirical Bayes (EB) prescribes that (1) we learn the Iu and go by maximizing the marginal log likelihood of the data, and (2) then return inference and v amounts to learning the induced postcrior by du and 9..
We mean-field variational FM, to solve this problem.
23

---

Published in Transactions on Machine Learning Research (12/2024)
MovieLens- 1M: row prior var 10
MovieLens- 1M: column prior var 10
dout log likelhhod
2.5
dout log likelihood
2.5
-3.0
Latent dim.
-3.0
30
64
Ru1322b: row prior var =
10
Ru1322b: column prior var = 10
Heldout log likellhood
Haldout log likelihood
-7.5
-7.5
-8.0
-8.0
-8.5
-8.5
9.0
99.
UserArtists: row prior var = 10
UserArtists: column prior var 10
Heldout log likelihood
Heldout log likellhood
-1400
-1400
-1600
-1600
-1800
GoodBooks: row prior var 10
GoodBooks: column prior var 10
Heldout log likellhood
Heldout log likelihood
-3.5
-3.5
-4.0
-4.0
-4.5
-4.5
-2
2
-2
row prior variance
column prior variance
TwinEB (PMF)
PMF
Figure As in Figure 3, but with additional values latent dimension L.
Define the ELBO, F(gus Sus Iv, Bv) our objective to be maximized as follows:
Flyu Iu: Bu.guv.ua Iv) = / q(u, v)
log log p(X,u, u. Ju.gv) gv. dudy q(u, v)
(47)
Assume L, the dimension of the row and column factors 1. Initialize Jus Dus u. and U, set = and until convergence do:
t++
(48)
(8) (t)
9v + arg max F F(gu, Tu: ol ((-1)) (4-1),g(0-1 (01-11)
(49)
qusQu
4 99 +- arg max F(g (t)
(8) (t)
99 qu: (v:
(50)
qo.ge
24

---

Published in Transactions on Machine Learning Research (12/2024)
L = 32
L = =4
Test HOLL
Test HOLL
-2.275
Test HOLL
-2.280
Test HOLL
2.19
-2.285
-2.20
-2.290
-2.21
15
32
32
#Row components (log)
#Row components (log)
#Column components: 1
#Column components: 2
#Column components: 4
#Column components: 8
#Column components: 16
#Column components: 20
#Column components: 32
Figure 7: As in Figure 4, but with additional column components as well as latent dimension 32.
E.2 Empirical Bayes normal means
For known = [4:], observations x = [xy]] are assumed to be distributed as follows:
g.geg
(51)
2; 0 id N (2,00,,0) 8:).
Empirical Bayes procodure prescribes two step process for inference over 0 [0,]:
1. Find a prior s(.) that maximizes the marginal likelihood of the data L(x):
9 arg max ((x) arg max II / g'(0.)p.: d.)d8.
(52)
2. Compute the induced posteriors by g:
P4(0
x,
II g(0:),,, (z: 00))
(53)
Solving this problem be seen mapping from the observations x and to the estimated prior ğ and induced posteriors (O x), that is:
(0(0),x(01 x)) + EBNM(x,s)
(54)
It turns out that in the variational EM algorithm, only the first two moments of the induced posterior are necessary, so we define them:
0, Eg 10
(55)
02, = Eg 103
x, g]
25

---

Published in Transactions on Machine Learning Research (12/2024)
E.3
Main proposition
arg max Fluu Iur Au: 10 166 1)(
(56)
where
Σ Yighy
ú(Y, w)
(57)
Σ3 Wj
-0.5
((w)
Σ W
(58)
Intuitively, the à, is normally distributed centered at uj, since for a fixed v we have:
x,: u,y Ei.: 7
0. Xi, (wo")-1
(59)
0 Σ₁ Xajuz
Σ oz
where we have computed independent regressions of the D-vector Yt, on the D-vector v. Given the normally distributed noise E, 0e ug ~ VWg 1).
E.4 Proof sketch
We show that maximizing the prior in EBNM is maximizing some FNM objective, We then show that by completing the square, we are solving the EBNM problem.
For the EBNM model of Equation (51), define the ELBO:
FNM (ge. 9) = = (o) log p(x, 0, l
(0) d0.
(60)
Note that since
FNM (g0.9) g) log Px (9) - -KL (( 9) Po x, x.9).
(61)
then
max N (qo.g) p(x g))
(62)
90
arg max FN (go.9) arg max log p(x g)-
(63)
go
Now we show that the complete log likelihood of the data in the EBNM model has squared form:
1
logps(x.00 0
2 (Iz 0,)² + (0) + const.
(64)
j=1
Therefore
1
+ Egg log 9(0)
PNN (go.9) 9) Egg
Σ
(Az03 -2B0,,)
(0) + const.
(65)
2
1=1
26

---

Published in Transactions on Machine Learning Research (12/2024)
Now to find the updates in Equation (50), we write down F() (the ELBO), then complete the square for (v) as needed to find the form of x A, and s(.) come from.
B, for the input to the EBNM problem. That where functions ú(.)
In their experiments, they set Gu = G. = 9, and try two settings for the family of the prior, namely, G E (SN, PN): the scale mixture normals, and a mixture of point mass at zero and a standard normal distribution. g = SN = En To,(0, σx) and G = PN = modo + (1 (0), For SN, Σb Th = 1 and ob are values lying on a user specified grid.
This paper is somewhat sparse in details. For instance, the scale minture of normals not explicitly defined anywhere in the paper. More details found in the first author's PhD thesis Wang (2017). The framework is implemented in the flashr R package.
Extensions to Binary, Binomial and Poisson data are explored in the first author's thesis Wang (2017) follows.
E.5 Extension to binary observations
The trick is to reparameterize the mean as log-odds, and then use a lower-bound trick for ELBO Jaakkola & Jordan (2000).
Assuming that Xis (-1,+1}}
Xij = 2 x Bern(pis)
= p(Xi, 2iz))
Pij 2ij
log - Pij
(66)
Zij = U₁Uj
4 ~ (.(
8/ (.).
We write the likelihood as:
1
P(Xg u. = P(Xis 2n)) = h(Yig,
=
(67)
1 + exp( eexppxX'pp (1)zz))
where 1(2) has (Light) lower bound with parameter E::
1(=) 2 M(&)) exp( E.
-
1
T(E=)
(h(Ex)
2
(68)
2
tanh(Ez/22
48.
Then replace P(X u, v) in F' with its lower bound H(u, X,E) and call this function Q(.) where H(). is defined as:
X Xijuay - Eij + - -61
H(u, v, X,g) = exp
(69)
2
where & = [Eey] E NDD are variational parameters. This new objective Q(.) has the following form:
Q( E, log (u, X,)) E, log p(u gu)p(v gv) E, log gu (u)
(70)
27

---

Published in Transactions on Machine Learning Research (12/2024)
We now maximize Q(1) with the following updates for the EBNM problem:
2j - 0.5 (X
(T(s))
(71)
8) 2 (S,,
E.6
Extension to Binomial observation
For Binomial data, they similar generative model to the Bernoulli case, Equation (66), but change the likelihood function to Binomial:
X;j 2 x Bin (niy.
(72)
They use an augmentation strategy with the Polya-Gas Gamma random variable. turns out this augmentation gives closed form for the expected values of the augmented r.v.s. and the updates identical to that of the Binary case.
E.7 Extension to Poisson observation
This is appendix B.2 the author's PhD thesis. The details are not given.
The generative model as follows:
Xi nij. Pij
Bin(nig: Pij)
Pij
log
Zij
(73)
Pij
R4 Poisson(Nig).
F Prior matching da Silva et al. (2023)
Summary. Given some target virtual statistics (say empirical moments of the observations), (da Silva et al., 2023) picks the hyper-parameters of the prior (A) such that the difference between the virtual statistios of the prior predictive distribution and the given virtual statistics are minimized.
Generic Bayesian matrix factorization takes the following form:
Row parameters Osr F(Jo,, o))
Row parameters Byk ~ F F( O()
Likelihood Yg ~ Fy Σ 0u
(74)
E[ Σ Og.3ge.
Poisson matrix factorization with Gamma priors takes the following form:
0sk iid F(us, oz) = Gamma b),
Byk F(Ug, o = Gamma(c, d),
lid
(75)
Yı Poisson
Σ
28

---

Published in Transactions on Machine Learning Research (12/2024)
where A4 b. og
1B 2. 0.03 -
b
A hierarchical Bayesian model p(Y, Z; X), has the following predictive distribution (PPD):
p(Y: / p( Z;)
(76)
Idea. Use Equation (76) to integrate out the parameters and search for hyper-parameters that match the data distribution well. Specifically, find a point estimate of À to used in posterior inference later.
Let (Y) denote the true population distribution of the observations.
(V) A(Y;
(77)
= f p(Y Z; A))(Z;))dZ
(78)
Note: in exponential family models with conjugate priors, an analytical form of PPD is available and the marginal likelihood of the data can be directly optimized to find solution for A.
Algorithm. Instead of using the true observation distribution pY(Y), use the observations virtual statistics of Y, fx.
Find such that Xx match target statistic T
In the paper, the authors use moments (mean and variance) of the observed data Y. First show closed form solution for PMF, Then provide more general case based gradient descent.
F.1 Analytical solution for Poisson matrix factorization
Given the PMF model in Equation (75), define:
Mean E[Yi;;
(79)
Variance V[Yiji
(80)
Covariance = plY.j. Yut:
(81)
hyper- parameearssssrss À No, oî uB, 03).
(82)
Proposition 1 For any entry of the virtual data matric Y [Yij} E NRM the mean and variance is given by:
B[Y A] K Kuaug
(83)
V[Yi; K [pana (1g0r)² (1o008)² (anas)²]
(84)
Proposition 2 For any pair of entries
Yig and Ya Y, their correlation is given by:
0 if i and j }
1
if i = t and j = iii-h t and j
PY Ynl
(85)
P1
P2
if i and j
where P1 = 3( [a[Yo]] and P2-
This yiclds the following for K:
TV EA
Ex
K
(86)
PIP2
[X[Y;]]
29

---

Published in Transactions on Machine Learning Research (12/2024)
for 1 - (p1 + p2).
For the Gamma hyper params, we'll have:
P2 V xY
TV Waly Ex(Yol
PI V
(87)
AWAl -By
bd EX ac
Vx(Yy)] PIP2
Given an observed data matrix Y, the methods can be used as follows:
1. Compute empirical mean Vx(Ya) and variance E [Yy] and covariances Åx Dy, Yul from Y.
Replace the virtual statistics Equations (86 and 87) with their empirical ones compute the Gamma hyperparameters a, c, and d.
Computing Px Vij, Yul from the single observed matrix Y not obvious, so we use the following estimator (Algorithm da Silva al. (2023)):
E Sx2 and B E Rx2
2. For s E S}:
3. Randomly pick row i, then randomly select two distinct indices in that row j, k:
As,1 Yijj Ax,2 + Ye,s
5. Randomly pick a column j, then randomly select two distinct indices in that column i, k:
6. B,, yesr Ba +- Kj
P1 < cor(A,11 A:2) and pa +- cor(Bs B,,, B.2).
8. Return and p
F.2 Gradient based optimization
When closed form solutions of the moments of the PPD are not available, the hyper-parameter optimization can be done as follows:
min d(T", Px),
(88)
X
for given (expert determined) T and f = Ex(g(Y)), and specified discrepancy measure d.
For instance, for target mean E" and variance V*:
d= (E" E (V*
s(Y) = (Y,Y²),
(89)
T(E,, Ez) (Ei, Es₂ E?)
The idea of optimization is: iteratively, draw samples from the PPD to estimate the virtual statistics Zx and solving for À by optimizing Equation (88).
To approximate the PPD moments, use Monte Carlo estimation and stochastic gradient descent, with repa- rameterization gradients (or REINFORCF), assuming that both d(.) and 7x are differentiable.
30

---

Published in Transactions on Machine Learning Research (12/2024)
For PMF, they have:
Ex[Y,]
~
Σ
Σ
B|y | ((e,,
SoSB
corppo ca
(90)
Ex[yy 31 = C
Σ
y.
yr-Poisson (0T B)
For given dataset, estimate E and V from the data, then optimize Equation (88) itcratively using GSD, where d(..-) is defined in Equation (89) and the moments in Equation (90).
31

---

Published in Transactions on Machine Learning Research (12/2024)
MovieLens- 1M
centered (70, 100)
#Column components
0.005
Test HOLL
10
0.000
30
50
-0.005
70
100
-0.010
150
1
10
30
50
70
100
150
#Row components
Ru1322b
centered at (70, 100)
0.2
#Column components
Test HOLL
0.1
70
0.0
100
150
10
70
150
#Row components
UserArtists
centered at (70, 100)
#Column components
Test HOLL
10
30
50
-50
70
-100
100
150
10
30
50
70
100
150
#Row components
GoodBooks
centered at (70, 100)
#Column components
Test HOLL
0.00
10
30
50
-0.05
70
100
-0.10
150
10
30
50
70
100
150
#Row components
Figure 8: TwinEB is robust to the number mixture components on real-world datasets. Test hold-out likelihood as a function of Kr, the number of row components (X-axis) and Ke, the number of column components (hue) for TwinEB with the Poisson likelihood. Each experiment from 10 seeds. The test hold-out likelihood is centered at the result of (K,, K2) (70, 100), highlighted by the gray dashed lines. Note that the range of the vertical axis is orders of magnitude smaller than the range of the real world experiments, suggesting that TwinEB is not sensitive to the choice of the number of mixture components.
32


Reviews:
[{"review_id": "MJMdOCKrog", "reviewer": "Reviewer_XJwX", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### **Strengths**\n\n1. **Novelty**: The introduction of twin population priors is a novel contribution that provides a principled solution to the problem of prior selection in Bayesian MF. The concept of empirically learning priors while leveraging the natural separation between rows and columns in a matrix is intuitive and effective.\n\n2. **Flexibility via Mixture Distributions**: By using mixture distributions for priors, the model can capture a wide range of data distributions, making it highly adaptable to various types of data and patterns.\n\n3. **Simultaneous Prior and Posterior Learning**: The variational inference approach allows for simultaneous learning of priors and posterior distributions, which reduces the need for manual tuning of hyperparameters and makes the method more robust.\n\n4. **Strong Empirical Results**: The authors provide convincing empirical evidence across different datasets, showing that their method outperforms traditional manually tuned priors and other empirical Bayes approaches. This is particularly relevant in scenarios with sparse data, such as movie recommendation systems, where the method demonstrates superior predictive accuracy.\n\n---\n\n### **Weaknesses**\n\n1. **Notational Inconsistencies**: The paper suffers from notational ambiguities, particularly in using $ \\theta_r $ and $ \\theta_c $, which denote distributions and latent variables. A clearer distinction between hyperparameters, distributions, and latent variables is necessary. For instance, using $ \\phi $ for hyperparameters and reserving $ \\theta $ for latent variables would improve readability. \n   - In equations like $ U_i \\sim P_{\\theta_r}(U) $ and  $ V_j \\sim P_{\\theta_c}(V) $, the parameters $ \\theta_r $ and $ \\theta_c $ are used both as the distributions' parameters and the latent variables' labels, which creates confusion. A more typical approach would be to separate hyperparameters and variables (e.g., using different symbols or subscripts).\n   - In other places, such as the factorization of the joint distribution (Equation 4), using $ \\theta $ as global and local latent variables is problematic. This contradicts typical Bayesian literature conventions, which recommend distinguishing hyperparameters (e.g., $ \\phi $) from random variables (e.g., $ z $).\n   - Potentially avoid using subscripts to introduce variables rather than using them to index variables' dimensions. The separation between random variables as conditional and hyperparameters could be achieved by using `|` and `;.` \n\n2. **Limited Comparison with Other Flexible Priors**: The paper does not sufficiently compare its method with other approaches incorporating flexible priors through side information, neural networks, or decision trees. Several methods, such as **Neural Poisson Factorization** (Ngo et al., 2020) and **MFAI: A Scalable Bayesian Matrix Factorization Approach** (Wang et al., 2024), use advanced techniques like neural networks and gradient-boosted trees to enhance prior specification. A discussion or comparison with these methods would provide a more comprehensive evaluation of the method’s strengths and weaknesses.\n\n3. **Absence of Comparison with Flow-Based Methods**: The paper does not explore flow-based methods, such as **Collaborative Autoregressive Flows (CAF)** (Zhou et al., 2020), which could offer a more flexible way of capturing complex posterior distributions. While the authors highlight the advantages of EB approach, a comparison with flow-based models that represent non-linear user-item interactions would provide a deeper understanding of the trade-offs involved.\n\n4. **Explanation of Latent Dimensionality Selection**: The rationale behind selecting the latent dimensionality $ L $ is not clearly explained. It would be useful for the authors to discuss how the dimensionality was chosen and whether methods like prior predictive statistics suggested by da Silva et al. (2023) could optimize the selection of $ L $. This would provide further rigor to the experimental setup.\n\n5. **No theoretical insight on the role of number mixture components**: some theoretical insight about the role the number of mixture components plays in properties of the posterior or prior predictive distribution could help understand the choice of such parameters.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces **twin population priors**, a novel approach for determining priors in **Bayesian Matrix Factorization (MF)** within an **Empirical Bayes (EB)** framework. Traditional methods for setting priors often rely on fixed parametric distributions (e.g., Gaussian, Gamma) or cross-validation, which can be computationally expensive and may not capture the data's complexity effectively. To address this, the authors propose data-driven priors for both row and column latent variables in the matrix, which are learned directly from the observed data. \n\nThe main contributions of the paper include:\n\n1. **Twin Population Priors**: A framework that separately learns empirical Bayes priors for both row and column latent variables, aligning the model’s marginal distributions with the true population distribution.\n   \n2. **Mixture Distributions for Flexibility**: Both row and column priors are modeled as mixtures of parametric distributions, offering a flexible prior that can capture more intricate patterns in the data.\n   \n3. **Joint Optimization via Variational Inference**: The authors develop a variational inference algorithm that jointly optimizes the twin population priors and the posterior distributions. This approach allows efficient learning of both the priors and the latent variables.\n\n4. **Comprehensive Empirical Evaluation**: The method is evaluated on multiple real-world datasets, including movie ratings and gene expression data, demonstrating improved performance over traditional parametric priors and comparable EB methods. Particularly, the twin population priors show robustness in sparse data scenarios.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. **Critical**: \n   - **Clarify Notation**: The notational ambiguities should be addressed to make the paper more accessible. The distinction between hyperparameters, random variables, and distributions must be clearly defined. Adopting more standard Bayesian notations would enhance the clarity and understanding of the model’s structure. A section in the background section or introduction clearly defining the notation used is important in this case (for example, observed variables, latent variables, hyperparameters, pdf/cdf, sampling distribution, empirical distributions). \n\n2. **Important**:\n   - **Expand the Discussion on Flexible Priors**: The literature review and discussion should incorporate a broader range of methods for flexible prior specification, including models that leverage side information, neural networks, and decision trees. This would give the reader a more holistic view of the field and place the contribution in the context of existing techniques.\n\n3. **Important**:\n   - **Comparison with Flow-Based Models**: The authors should include a comparison with flow-based methods, such as CAF, which offer flexible ways to approximate complex posterior distributions. This would provide valuable insights into the strengths and limitations of the variational inference approach used in the paper. Discussing the differences and limitations is relevant if such a comparison is impossible. \n\n4. **Important**:\n   - **Justification for Latent Dimensionality**: The paper should provide a more detailed explanation of how the latent dimensionality $ L $ was chosen. Incorporating prior predictive methods for optimizing $ L $ would add methodological rigor and strengthen the empirical results.\n\n5. **Important**\n- **Theoretical Discussion on Prior Predictive Properties**: A theoretical discussion on the prior predictive properties of the twin population priors would significantly strengthen the paper’s theoretical foundation. This could include calculations based on mixture distributions, as seen in prior work (e.g., Silva et al., 2023) combined with known moments for mixtures (https://math.stackexchange.com/questions/195911/calculation-of-the-covariance-of-gaussian-mixtures), and would also provide an opportunity to introduce a new baseline—**twin priors with moment matching**. Such a discussion would offer insights into the behavior of the priors before data is observed, which is crucial for understanding their impact on model performance.", "broader_impact_concerns": "The paper should mention potential ethical concerns in the context of recommender systems applications. Specifically, empirical priors derived from the data may inherit and amplify biases present in the data, which could lead to biased outcomes in real-world applications. Including a discussion of these potential biases and suggesting mitigation strategies (e.g., fairness-aware priors) would be a valuable addition to the broader impact section.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "DwBZRxVdot", "reviewer": "Reviewer_zHrL", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: This seems to be a novel approach, and the no-tuning feature is appealing.\n\nWeakness:  Lack of scope and clarity in numerical studies, so it is unclear how much advancement this method brings 'in performance' compares to existing approaches.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a new method for simultaneous probabilistic matrix factorization and prior estimation by maximizing ELBO with respect to both posterior distribution and prior parameters. The author claims that the algorithm performs better than manually tuning the hyper-parameters for the prior or selecting priors with Empirical Bayes.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "These changes are critical to securing the recommendation for acceptance. More details and numerical experiments are needed to support the advance of this method.\n1.  There should be more comparison with other methods in simulation and on real datasets, including Wang & Stephens (2021), da Silva et al. (2023) that demonstrates the strengths and weakness of each method. Now the comparison with Wang & Stephens (2021) is concluded with one sentence \"flashr crashes on all but the least sparse dataset, Ru1322b.\" More insights are needed.\n2. The selection of the real datasets needs to be reconsidered. All real datasets have less than 20% of non-zero elements. Why apply Gaussian models to them (Section 4.3)? \n3.  The methods under comparison have not been clearly described. Is 'row prior variance' the hyper-parameter for GMF? How is it used in GMF?\n4. \"We note that it is not straightforward to compare our method to that of Jiang & Zhang (2009);\" Is Jiang & Zhang (2009) a method for matrix factorization?\n\nThese changes will strengthen the work in my view:\n1. The description of \"population prior\" is too lengthy on pages 2-4.\n2. Abbreviations appear without definition, such as HOLL.", "broader_impact_concerns": "None.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "3NsEtxqBj8", "reviewer": "Reviewer_YEoN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "(+) Great to see the explicit comparison to da Silva et al and to Wang and Stevens.\n(+) The approach is elegant, re-using the ELBO for prior and posterior optimization.\n(-) Empirical results are rather weak. For some examples other methods perform better (Fig. 2). The authors downplay this as \"comparable\".\n(-) Some design choices could be better-justified (see requested changes)\n(-) Limited methodological novely", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper introduces an empirical Bayes prior over factor matrices for matrix factorization, specifically Poisson and Gaussian matrix factorization. These priors are modeled as mixture distributions. Inference is performed using a variational objective. Conveniently, the ELBO allows for optimization of both the variational parameters of the factor matrices' posterior, as well as the parameters of the priors' mixture distribution. The method is tested on held-out log-likelihood in comparison to similar methods on standard benchmark data. It performs on par, sometimes favorably, sometimes worse than existing methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The manuscript needs to be proofread for comma punctuations. E.g. \"Here, U and V are [...] latent vectors.\"\n\n- In Fig. 4, the held-out performance increases with the number of components. What is the limiting behavior when using even more components? Does it match your expectations?\n\n- How do you decide on the number of mixture components in your main experiments? Does it matter?\n\n- In the discussion, I would like to see more of a critical assessment where this method fits in or improves over existing work.\n\n- I may have missed it, but why can you not optimize row and column priors simultaneously (or rather iteratively in the same loop), instead of holding one fixed?\n\n\n## Detailed Comments\n\n\n### Introduction \nI find the formulation in the introduction somewhat confusing. E.g. you say \"Here \\mathbf{U}_i and \\mathbf{V}_j are per row and per column specific latent vectors\". I get what this means from the context, the sentence itself is unclear. Similarly, you say \"An observed matrix of data X then defined a posterior distribution p(U, V|X) over the row variables and the column variables.\" What are the \"row and column variables\"? Why not just say \"... over the latent factor matrices\" or so. Am I missing some subtlety here?\n\n\"It is applicable for priors over variables that have repeated draws from them\" -> Sentence is unclear.\n\n\"The data that informs us about them is overlapping\" -> What does this mean?\n\n- [FYI] Another example of empirically set matrix factorization population priors: https://proceedings.mlr.press/v70/rukat17a/rukat17a\n\n\n### 2 Related work\nYou say that mixture priors are more expressive than hierarchical priors. This is not generally true.\n\n### 3 Empirical Bayes\nWhat is P*? You're saying the \"true (unobserved) distribution of observations\". That seems nonsensical or, at least, requires some explanation.\n\nI don't understand the point of section 3.1. Could you explain why we need this explanation?\n\nYou say \"there is no distinction between local and global random variables\" - I think you want to say \"latent\" instead of \"random\".\n\nWhy do you choose a log normal? \n\n### 4 Experiments\n\"We studied Algorithm 1 ...\" -> You probably want to write this and what follows in the present tense.\n\nThe train/test split procedure is confusing. Can you please explain and motivate the choices more clearly?\n\nSections 4.2 and 4.3 are named \"Results\", but also explain the experimental procedure.\n\n### 5 Discussion\n\nThe last sentence is not complete.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: After the rebuttal, all reviewers backed acceptance of the paper.  There was clear agreement that the updates provided have addressed the key concerns originally raised.  My own assessments did not bring up any reason to go against this unanimous conclusion, and I therefore recommend acceptance of the paper.   As there are no major outstanding issues, I recommend that the paper is accepted "as is".
Decision: Accept (poster)