Covariate-moderated Empirical Bayes Matrix Factorization
Authors: William R.P. Denault
Karl Tayeb
Peter Carbonetto
Jason Willwerscheid
Matthew Stephens
Venue: NeurIPS 2025 poster ()
Keywords: matrix factorization
non-negative matrix factorization
variational inference
empirical Bayes
spatial transcriptomics
Primary area: machine_learning_for_sciences
TL;DR: We introduce a modular matrix factorization framework called "covariate-moderated empirical Bayes matrix factorization" (cEBMF) that can leverage side information to improve the factorization through the use of adaptive priors.

Abstract: Matrix factorization is a fundamental method in statistics and machine learning for inferring and summarizing structure in multivariate data. Modern data sets often come with "side information" of various forms (images, text, graphs) that can be leveraged to improve estimation of the underlying structure. However, existing methods that leverage side information are limited in the types of data they can incorporate, and they assume specific parametric models. Here, we introduce a novel method for this problem, *covariate-moderated empirical Bayes matrix factorization* (cEBMF). cEBMF is a modular framework that accepts any type of side information that is processable by a probabilistic model or a neural network. The cEBMF framework can accommodate different assumptions and constraints on the factors through the use of different priors, and it adapts these priors to the data. We demonstrate the benefits of cEBMF in simulations and in analyses of spatial transcriptomics and collaborative filtering data. A PyTorch-based implementation of cEBMF with flexible priors is available at https://github.com/william-denault/cebmf_torch.

Covariate-moderated Empirical Bayes Matrix Factorization
William R. P. Denault, Karl Tayeb, Peter Carbonetto & Matthew Stephens
Departments of Statistics and Human Genetics
University of Chicago
Chicago, II. 60637, USA
Jason Willwerscheid
Mathematics and Computer Science
Providence College Providence, RI 02918, USA jvillver@providence. edu
Abstract
Matrix factorization is a fundamental method in statistics and machine learning for inferring and summarizing structure in multivariate data. Modern data sets often come with "side information" of various forms (images, text, graphs) that can be leveraged to improve estimation of the underlying structure. However, existing methods that leverage side information are limited in the types of data they can incorporate, and they assume specific parametric models. Here, we introduce novel method for this problem, covariate mederiate- empirical Bayes matrix factorization (cEBMF). cEBMF a modular framework that accepts any type of side information that is processable by probabilistic model or neural network. The cEBMF framework can accommodate different assumptions and constraints the factors through the use of different priors, and it adapts these priors to the data. We demonstrate the benefits of cEBMF simulations and analyses of spatial transcriptomies and collaborative filtering data. A Py TyTochh- based implementation of cEBMF with flexible priors is available at https:// gi thub. com/wt 1llliimm www.dddenanauuuueeeee torch.
Introduction
Matrix factorization methods, which include principal component analysis (PCA), factor analysis, and non-negative matrix factorization (NMF) [1-3], very widely used methods for inferring latent structure from data, performing exploratory data analyses, and visualizing large data sets (e.g. [4-6)). Matrix factorization methods also instrumental in other statistical analyses such adjusting for unobserved confounding [7]. When factorizing a matrix, say Z, the matrix may be accompanied by additional factorization algorithm toward a more accurate or interpretable factorization. A recent prominent example of this in genomics research is spatial transcriptomics data [8], which is expression profiled in many genes at many spatial locations ("pixels") [9]. For a variety of reasons, one typically seeks factorize Z. the matrix of gene expression profiles. But the 2-d coordinates of the pixels also provide important information about the biological context of the cells; for example. we might expect nearby pixels to belong to the same cell type or tissue region. Therefore, "spatially aware" matrix factorization methods have recently been proposed for spatial transcriptomics data [10-12].
column data-"side information" --that may be able "guide" the matrix
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

Ground-Iruth labeling on 2-d co-ordinate system
G cEBMF prior
B PCA
C NMF
k=1 1
1.00
200-
1.04
0.75
100
POP
PC2
0-
0.50
0.25
-100-
1.0
k= 2
-200-
1.0
0.00
0.25
0.50
0.75
1.00
-5
0
-200
100
D. PC1
100 200 300
0.5
D EBMF
E Spatial PCA
F cEBMF
0.0
0.02
0.0
0.5
1.0
1
0.4
k=
PC2
pCP
b-
PC2
1.0
0.5
-0.4
-0.04
0.0
-0.08
-0.04
0
0.5
1.0
PC1
PC1
PC1
Figure 1: Toy simulation illustrating cEBMF for learning matrix factorization, LF". In this example. Z is a 1,000 x 200 matrix. Each of the = 1,000 data points is assigned to one of three clusters (orange, light blue, dark blue). Points near each other tend be assigned the same cluster, except near the boundaries (A). Without the side information (the 2-d coordinates in A), PCA, NMF and EBMF with K 3 factors cluster some points accurately, but many other points are not clustered accurately (B-D). By contrast, Spatial PCA [11] and our new method, cEBMP, which both incorporate the side information into the prior, more accurately cluster the points (E, F). (For consistency of visualization, the top 2 PCs of the L matrices from NMF, EBMF and cEBMF are shown.) Spatial PCA assumes the data points are spatial, whereas cEBMF does make this assumption; instead, it has flexible prior that is adapted to the data. This learned prior is shown in G: the color of the points shows the prior probability that row 1, column & of L is nonzero (blue = low prior probability, red high prior probability). See Sections 3 and 4 for definitions and additional details.
In this paper, we describe a novel matrix factorization framework that allows high-dimensional row and column data to guide the matrix factorizations without having to make specific assumptions about how these data inform the factorization. For example, although our framework can be applied to data that exhibit spatial properties, it does not assume or require that the data be spatial. Our framework is also flexible in that includes many existing approaches as special cases, including unconstrained matrix factorization [13, 14], non-negative matrix factorization [15], semi-non-negative matrix factorization [16], and more recent methods that incorporate side information [17]. These features are achieved by taking an empirical Bayes approach, building on the recent empirical Bayes matrix factorization (EBMF) framework [13, 14]. In particular, we extend the EBMF approach of [13] with adaptive priors that are modified by the side information. We call this approach "covariate- moderated empirical Bayes matrix factorization," " or "cEBMF" for short. See Fig. I for toy example that illustrates the key features of cEBMF.
Related work
The literature on matrix factorization methods that incorporate side information is quite extensive. The different methods make different modeling assumptions, and are typically motivated by certain types of data. Although it is not possible to review all relevant literature here, we discuss a few of the most important or related methods.
Several variants of the topic model- can be viewed as matrix factorizations with "sum-to- one" constraints on L and F [18] incopoaceee side information in different ways; for example, the correlated topie model [19] and the structural topic model [20] incorporate document-level side information into the priors L. Collective matrix factorization (CMF) [21-23] has gained

---

considerable interest, but CMF is based on ideas that are quite different from cEBMF; like cEBMF, CMF assumes that the side information is the form of a matrix, but unlike cEBMF, CMF assumes that the side information factorizes in a similar way to Z. Clearly, this assumption will not make sense for some types of data. Another prominent theme in matrix factorization with side information incorporating group-level or categorical information, including ontological data. Among the methods in this area are CTPF [24] and the method of [25]. Another important class of methods related to cEBMF from the deep learning literature are variational autoencoders (VAE) [26], conditional variational autoencoders (cVAE) [27] and neural collaborative filtering (NCF) [28]. These methods generalize the concept of matrix factorization to nonlinear embeddings.
The method that is most closely related to cEBMF MFAI [17] (see also [29] for related ideas). MFAI is in fact special case of cEBMF in which the priors on F are normal and the prior means are informed by the covariates. Like cEBMF, MFAI allows these priors to adapted separately for each dimension k. However, MFAI is not nearly as general as cEBMP; it implements only a single model, a single prior family with a specific parametric form, a specific procedure for fitting these priors (using gradient boosted tree methods [30]), and it only accommodates row-wise side information.
Several matrix factorization methods have been developed specifically for spatial transcriptomies data. Spatial PCA [1I] models the spatial similarity among rows of L using Gaussian process prior. (Spatial PCA is similar to GP-LVM [31]. See also [32].) An NMF version of this approach generates "parts-based representations" guided by the spatial context of the data points [12]. More recently. IRIS [33] regularizes the matrix factors through a penalty function that encodes the spatial information in graph (see also [34]).
3 Covariate-moderated empirical Bayes matrix factorization
3.1 Background: empirical Bayes matrix factorization
Empirical Bayes matrix factorization (EBMF) [13] is a flexible framework for matrix factorization: approximates matrix Z E Raxp as the product of two low-rank matrices,
Z = LF".
(1)
where E Ra F E RoxK and K 2 1. (In our applications, K < n,,p. P.) EBMF assumes a normal model of the data,
Z L= EE, Ei NOO TT),
in which N(,, o denotes the normal distribution with mean fr and variance o², and the residual variances Tiy may vary by row (i) or by column (j) or both. (EBMF, and by extension cEBMF, also allows to contaín missing values [13], which important in many applications of matrix factorization, including collaborative filtering: see Sec. 4.2.) EBMF assumes prior distributions for elements of and F, which are themselves estimated among pre-specified prior families Gr and Gyku, respectively:
Lu 9k
9k E Ge
k == I,.., K
Jik 9k
(
g
(
k= K.
The flexibility of EBMF comes from the wide range of possible prior families (including non- parametric families) [35]. Different choices of prior family correspond to different existing matrix factorization methods. For example, all families Gr and Gy are the family of zero-mean normal priors, then LFT is similar truncated singular value decomposition (SVD) [36, 37). When the prior families are all point-normal (mixture of a point mass at zero and a zero-centered normal), one obtains empirical Bayes versions of sparse SVD sparse factor analysis (38-40) The prior families can also constrain L and F; for example, families that only contain distributions with non-negative support result in empirical Bayes versions of NMF. In summary, EBMF (2-3) is a highly flexible modeling framework for matrix factorization that includes important previous methods as special cases, but also many new combinations (e.g., [41]).
3.2 The cEBMF modeling framework
In covariate-moderated EBMF (cEBMF), we assume that we have some "side information" (covari- ates) for rows and/or columns of Z [42, 43]. Let æ; denote the available information for the ith row
3

---

of Z, and let y, denote the available information for the jth column of Z. In principle, 2; and y can be any information processable by a neural net (text, graph, image, other structured data), but for simplicity we assume that this information is stored as a matrix. Therefore, let X E R be a matrix containing information on the rows of Z. with æ; corresponding to the ith row of X (e.g.. 2; might contain the 2-d coordinate of cell i). Similarly, let Y E T2xx contain information on the columns of Z. with y/, corresponding to the jth row of Y. In cEBMF, we incorporate this side information into the model through parameterized priors.
la 9k (ar.);
(6)
9k
(0) (xx.) E ge,ks
k 1......, K k k= I........AA K,
fzx 9k (gn),
(1)
(1) (2:) E Gymss
9k
where (x.)) is a probability distribution within the family Ge,s. parameterized by Hú and ge (f) (y,)
is a probability distribution within the family GI parameterized by y:-
A key limitation of many existing approaches is that they integrate the side information using restrictive parametric models that may or may not be appropriate for the particular application. Another limitation that the priors chosen for these methods may make strong or perhaps unrealistic assumptions about the structure underlying the data; for example, Gaussian process priors, which have been used in matrix factorization (e.g.. [11, 31, 44),, typically assume that the factors vary smoothly in space, which makes it difficult to capture sharp changes at boundaries [45]. Existing methods also typically rely on hyperparameters that need to be tuned or selected (c.g., using cross-validation).
To address these issues, we propose cEBME, method that:
1. Can leverage a large variety of models (e.g., multinomial regression, multilayer perceptron, graphical neural nets, convolutional neural nets) to integrate the side information into the prior.
2. Can use families of priors that are flexible in form and thus do not make strong assumptions.
3. Allows automatic selection of the hyperparameters in (4) via an empirical Bayes approach.
More formally, we fit a prior for each column & of L, which maps each vector of covariates æ; to a given element 9% (0) (ax) E Ge,ks, and similarly for each column & of F. In Sec. 3.3, we describe simple yet general algorithm that simultaneously learns the factors L, F and the priors 9k (æ;)). 9k. (y,)) A PyTorch-based [46] implementation of cEBMF with several different parameterized prior families is available at https: / //tpt/gtttt.tt////////// cwm/wimmw........ 11 1iam-maaai denault/cebmt ttrch
3.2.1 An illustration: cEBMF with side information on factor sparsity
Here we illustrate the implementation of the cEBMF framework using a simple yet broadly applicable prior family. This prior family assumes that the covariates X, Y only inform the pattern of sparsity- that is, the placement of zeros-in L and F. This type of prior is of particular interest for matrix factorization because matrix factorizations are typically invariant to rescaling, and therefore priors that inform the magnitudes of la and fjs are difficult to design. (By "invariant rescaling," we mean that
the likelihood or objective does not change if we replace LFT by LF", where L LD. = FD= and D is an invertible diagonal matrix.) We define this prior family as
Ges = {g: g(u) (1 π(x ()) + T 0)x,0)g1 (w:w),
(5)
in which (o(n) denotes the point-mass at zero, 9n (w; w) denotes the density of some probability
distribution Ø1 (w) on E R. and Rm denotes the covariate. For example, when Φi is the normal distribution and w specifies the mean and variance, (5) is family of parameterized "spike-and-slab" priors [47], and cEBMF with Ge Gs. Gy, h Gss implements version of sparse factor analysis [38-40] in which the sparsity of the factors is informed by the covariates. (Note that the "ss" in G. is short for "spike-and-slab.") Alternatively, if 91 is a distribution with support only on non-negative numbers, such as an exponential distribution, then cEBMF implements version of sparse NMF. The free parameters are 0, which control the weight on the "spike", ò0, and w, which control the shape of the "slab", 91. One simple parameterization of π (xx 0) uses a logistic regression model,
न (x, 0) = (0u + Σm x,O2).
(6)
where (bxx) := 1/(1+e=> denotes the sigmoid function, and 0 E Rm+11 Most of the parameterized prior families used in this paper and in the cEBMF software are variants elaborations on Ges.

---

3.3
The cEBMF learning algorithm
key feature of the cEBMF modeling framework that the algorithm for fitting the priors and estimating the factorization is simple to describe and often straightforward to implement. In brief, the cEBMF learning algorithm reduces complex model-fitting task series of simpler subproblems. Fach of these subproblems involves fitting a covariate-moderated variant of an empirical Bayes normal means (EBNM) model [35]. This also has the advantage of making the cEBMF framework and software modular so that method that solves covariate-moderatee EBNM problem can be "plugged in" the generic cEBMF algorithm.
3.3.1 Background: empirical Bayes normal means
Given observations 3, E R with known standard deviations 8; 0, 1,..., n, the normal means model [48-50] is
B ind N ((, s
(7)
in which the "true" means & E R are unknown. It is further assumed that the unknown means are
B iid.
g E G.
(8)
where G is some pre-specified family probability distributions.
The empirical Bayes (EB) approach to fitting the normal means model (7-8) exploits the fact that the noisy observations 3, contain information not only about the underlying means 3; but also how the means are collectively distributed, EB "borrows information" across observations estimate g: typically this is done by maximizing the marginal likeeeili of (7-8). The unknown means 3; are typically estimated by their posterior means (given the estimate of g).
To adapt the EBNM model (7-8) to cEBMF, we allow the prior for the ith unknown mean to depend on additional data d, and parameters 0,
ind. ~ g(d, 0) E G,
B
so that each combination of 0 and d; maps to element of G. We refer to this as "covariate- EBNM" (cEBNM). Solving the cEBNM problem therefore involves two key computations:
1. Estimate the model parameters.
Compute
= argmax C(0), 0 Rm
(10)
where ((0) denotes the marginal likelihood.
((0) == p(, 8, S, = II N(3 g(;; 9(( 0) d3,,
(11)
in which / denotes the density of (Bi, 8t) B. and g di, 0) denotes the density of (d, 0) Bi.
(B,....,
(SI...., S.)), D is matrix storing di,..., dn. N (8:; Bi, sf)
2. Compute posterior summaries. estimated prior,
Compute summaries from the posterior distributions given the
9 8f) de Q) p(B: | [,, 84. 0,D) J N3.;1, g( d. 0)
N(A. B,
(12)
For many classical prior families, such as the spike-and-sla and-slab family in Sec. 3.2.1, the integrals in (11) and (12) can be computed analytically. More generally, standard numerical techniques such as Gauss-Hermite quadrature may provide reasonably fast and accurate solutions for prior families that do not result in closed-form integrals since the integrals in (11) and (12) are one-dimensional. As result, Ô can often be obtained efficiently using off-the-shelf optimization algorithms even when the chosen priors do not result in analytical integrals.
In summary, solving the cEBNM problem consists of finding a mapping from the known quantities (3, s, D) to tuple (0, Q). where each (d, 0) maps an element g(d,, 0) G, and q is the posterior
5

---

distribution of the unknown means, ((33 : IT. p(6; Br, st, 0, D) To facilitate the description of the cEBMF algorithm below, we denote this mapping as
cEBNM( D) (6,
(13)
Note that in practice the full posterior ((3) is not needed; the first and second posterior moments sufficient (see Sec. 3.3.2). Any prior family is admissible under the cEBMF framework so long as the mapping (13) is computable (either numerically or analytically).
3.3.2
Algorithm
Given a method for solving the cEBNM problem (Sec. 3. the cEBMP model can be fitted using a simple coordinate ascent algorithm. In brief, the cEBMF algorithm maximizes an objective function-the evidence lower bound (ELBO) [51] under a variational approximation with conditional independence assumptions on L and F (see the Appendix) iterating over the following updates for each factor k = 1..... K until some stopping criterion is met:
1. Disregard the kth factor in R. the n x p matrix expected residuals, R = R R laf..
2. For each i ,..... ,t, compute the least-squares estimates of lik, denoted lax, and the standard deviations o/ of these estimates,
L = Ty $/1 = 1
(14)
(15)
where S5x and fike denote, respectively, the first and second posterior moments of fzk-
3. Update 9) Ge by solving (10), in which we make the following substitutions in (10):
B, +- lik, 3; + s/ I,... n, D + X,9 Ge,k-
4. Making the same substitutions in (12), update the posterior means la = (Li Ent
and posterior second moments L (L )
5. Perform updates similar to those in Steps 2-4 update fx. f2 and gk (f)
-5
6. Update the matrix of expected residuals, R R -+af%.
These steps are iterated until some stopping criterion is met. The algorithm must be initialized with initial estimates of L, F. The expected residuals are then initialized as R LF". Note that to simplify presentation we have omitted some details such as how update the residual variances Tij These and other details are provided in the Appendix.
4 Experiments
4.1 Simulations
To assess the benefits of cEBMF, we compared cEBMF with other matrix factorization methods in simulated data sets. We compared with several methods that do not use side information, including EBMF (flashier R package [13, 35]), penalized matrix decomposition ("PMD"; PMA R package (39]), and a variational autoencoder (VAE) [26] implemented in PyTorch [46]. We also compared with other methods that use side information, including MFAI (mfair R package [17].. Spatial PCA [11], conditional VAE (AV) [27], and neural collaborative filtering (NCF) [28]. cVAE and NCF were also implemented in PyTorch. Note that Spatial PCA accepts only specific type of side information, the 2-d coordinates of the data points, so was not included in all the simulations.
We compared the methods in four simulation scenarios designed to capture a range of settings where one might perform a matrix factorization analysis, with or without side information: (1) "sparsity- -riven covariate" setting in which the covariates only informed the sparsity pattern of L and F: (2) an "uninformative covariate" setting which the covariates provided no information about the true matrix factorization; (3) "tiled-clustering" setting in which L depended on the 2-d location of the data points; and (4) a "shifted tiled-clustering" setting which the cEBMF priors were unable recover the true data generating process. The latter scenario was used to assess cEBMF under model

---

misspecification. We simulated 100 data sets in each setting, and we assessed the ability of each method recover the true matrix factorization as measured by root mean squared error (RMSE) between the true matrix factorization LF" and estimated matrix factorization LFT, More detailed descriptions of the simulations and the methods compared are given in the Appendix.
The results are summarized in Fig. 2. cEBMF was gener- ally more accurate than the other methods, particularly when the covariates were informative; cEBMF achieved the greatest gains over EBMF the tiled- clustering setting where the co- variates were also the most infor- mative. Reassuringly, cEBMF did not perform worse than EBMF in settings with an unin- formative covariate or a prior that was misspecified ("shifted tiled- clustering"). The deep learning approaches (VAE, VAE,, NCF) generally performed worse than the other methods. cVAE some- times outperformed VAE when the side information was highly informative, such as in the tiled- clustering scenario, but did not provide improvements over VAE the more challenging sparsity- driven scenario. We also ran Spa- tial PCA on the cilst-clustering and shifted tilee-clustering data sets where the factors were partly driven by the 2-d locations of the data points. Despite the fact that Spatial PCA can exploit the side information, it had worse accu- racy than EBMF which did not use the side information. This may be because Spatial PCA makes assumptions (c.g.. orthogonal factors) that cEBMF except in the shifted tilee-clustering setting: MFAI is a much less flexible model than cEBMF and therefore its performance models were misspecified in the shifted tiled-clustering setting, but perhaps MFAI was the least misspecified.) Additional results including comparisons with other methods (PCA/SVD, Sparse SVD [38), CMF (21]) are in the Appendix.
Unintormative covariate
Sparsity-driven covariate
0.5
RMSE
0.3
RMSE
0.5
0.3
GRBMP
EBMP
MFAI
PMO
VAE
oYAE
NCP
GERMP
EBMP
MFAI
PMD
oVAE
NDF
Tiled-dustering
Shifted fled-clustering
2.0
30
RMSE
RISE
0.5
LEBMF
EBMP
MFAI
PMO
SpatiaPCA
VAE
oVAE
NCF
CEBMF
EBMP
MFAI
PMD
SpallalPCA
VAE
oVAB
NOR
Figure 2: Performance of the different matrix factorization meth- ods in the simulated data sets. Each boxplot summarizes the root mean squared errors (RMSEs) across 100 simulations in that scenario. (Lower RMSES are better.) See Figures 5-8 in the Appendix for additional results from the simulations. Note Fig. 1 shows results from one of the tiled-clustering simulations in detail.
not met by simulations. MFAI generally performed worse than EBMF and
sensitive to the appropriateness of its modeling assumptions. (All
4.2 Collaborative filtering
To provide quantitative assessment of the matrix factorization methods in real data, we ran the same methods on the MovieLens 100K data [52], a standard collaborative filtering benchmark in which the goal is to predict the unobserved elements of the matrix. Here, Z was 1,682 x 943 matrix containing integer-valued movie ratings, with rows corresponding to movies and columns corresponding to users. Since most (93%) of the movie ratings were missing, this example highlights the ability of these methods to handle missing data (unlike most NMF methods). The side information X was a 1,682 19 binary matrix containing information about the movie's genre (comedy, adventure, etc). We held out some of the moving ratings at random, and used these held-out ratings as a test set.
We ran EBMF and cEBMF so as to produce non-negative matrix factorizations, which is common in collaborative filtering (e.g., [21]). Therefore, the results we labeled these methods as "EBNMF" and "<EBNMF.. (Note that MFAI cannot produce non-negative matrix factorizations.) To enforce non-
7

---

training proportion 50%
training proportion 70%
training proportion 90%
1.1
1.0
1.0
1.0
RMSE
RMSE
RMSE
0.9
0.9
0.9
0.8
dEBNMF
EBNMP
MFAI
VAE
oVAE
NCF
CEBNMM
EBNMF
MFAI
VAE
cVAE
NCF
CEBNMF
EBNMF
MFAI
VAE
cVAE
NCF
Figure 3: Prediction performance of different matrix factorization and deep learning methods in the Mov 100K data [52]. Training proportion X% means that X% of the movie ratings were used in training, and the remaining (100- X)% were used to evaluate accuracy (measured using RMSE). The results at each training proportion are from 10 random training-test splits.
negativity in and F, we used mixture ex-eooenniill priors. For cEBNMF, the side information was incorporated into the priors on L using a multi-layer perceptron. (See the Appendix for further details.) Since we didn't know the true number of factors, K was chosen adaptively in EBMF, cEBMF and MFAL (We set an upper limit of 7 for EBMF and cEBMF, and 12 for MFAL)
The results are summarized in Fig. 3. Both MFAI and CEBNMF were able to the side information (the movie genres) to improve over EBNMF, and all three matrix factorization methods were more accurate than the deep learning methods. We conjecture that the deep learning methods would have performed better with more data (such as the more recently released MovieLens data sets that are much larger). On the MovieLens 100K data, cEBNMF yielded overall the best prediction accuracy across the different training-ser set splits.
4.3 Spatial transcriptomics
Although cEBMF was not specifically designed for spatial data, here we show that cEBMF also yields compelling results from spatial transcriptomics data [8] by exploiting the side information, the spatial locations of the data points. We illustrate this using a data set [53, 54] that has been annotated by domain experts and has been used in several papers to benchmark methods for spatial transcriptomics (e.g.. [11, 55-57J). The data were collected from 12 slices of the human dorsolateral prefrontal cortex (DLPFC). After data preprocessing, each slice contained about 4,000 pixels and expression measured in about 5,000 genes (n x 4000, p ~ 5,000).
In this application, our aim was to generate a srarts-based representation of the data, with the hopes that the "parts" would resolve biologically interpretable units (e.g., cell types, tissue regions, gene programs) [12, 58, 59]. This is a fundamentally different aim from the previous examples: in the previous examples, the goal was generate accurate low-dimensional representations, but did not ask whether the individual dimensions were accurate or interpretable. With this aim in mind, we to produce non-negative matrix factorizations ("CEBNMF") using the same priors that were used for the MovieLens data. We compared to two other non-negative matrix factorizations that did leverage the side information: NMF implemented in the R package NNLM [15], and EBMF with point-exponential priors ("EBNMF"). (The point-exponential prior is a simplification of the ofotuure-o-exppoenntaal prior with a single exponential component in which the rate parameter is also learned.). We also compared to several of the methods that were considered in the previous experiments, including methods such as Spatial PCA and eVAE that make use of the side information, and others that do not.
cEBMF
Spatial PCA deserves special mention because it was specifically designed for spatial transcriptomics data [11]. Although Spatial PCA does not produce a parts-based decomposition, the Spatial PCA software automatically clusters the data points after projection onto the principal components (PCs), and this clustering can be compared to the non-negative matrix factorizations. Following [11], we
We used the data prepared by the authors of the SpatiaLLIBD package [53] which were made available for download https: / httt//er.eee..ceeeeeeee 1ibd. org/ hetp//pataaallB/../
8

---

slide
human labeling
EBNMF
CEBNME
MFAI clustering
Isbel
factor
factor
cluster
Spatial PCA
cVAE clustering
NCF
eluste
eluster
cluster
cluster
slide 10 human labeling
labe
GEBNMF
clustering
factor
clustrr
Spatial PCA + clustering
VAE clustering
cVAE clustering
NCF clustering
diutter
cluster
cluster
cluster
Figure 4: Results on slides 4 (top row) and 10 (bottom row) of the DLPFC spatial transcriptomics data [53, 54]. For the NMF, EBMF and cEBMF results, each data point ("pixel") shown as a pie chart using the relative values of ith row of L (after performing an "LDA-style" post- processing of L, F [12]). (Higher-resolution versions of these images are available online at https: //github.co cetistuub.oom/eeee wwlllllmmdenalltl//mmmm produce a non-negative matrix factorization, we clustered the low-dimensional embeddings using the same approach that was used in Spatial PCA [11]. CMF results for these two slices and additional results on all slices are given the Appendix.
_errpeeeee__mre_t_eeeeeee Since MFAI and other methods did not
computed the top 20 PCs, then we ran the walk-trap clustering algorithm [60] on the PCs. Additional details of the Spatial PCA analysis and the other methods are given in the Appendix.
Figure 4 shows results on two of the slices, with additional results on all 12 slices provided in the Appendix. The manual annotations the left-hand side should be viewed useful reference point, but not necessarily the "ground truth". (Consider that the data-driven annotations might identify previously unknown cellular structures.) EBNMF, <EBNMF and MFAI adapted the number of factors K to the data (with upper limits of 50, 20 and 9, respectively). For the other methods, the number of clusters set match the manual annotation. Qualitatively, some of the factors from NMF and EBMF seem to correspond to the expert-labeled regions, but several other factors appear to be capturing other substructures that have no obvious spatial quality. Comparatively, the cEBNMF results in slices 4 and 9 capture the expert labeling much more closely, with most factors showing clear spatial quality. The clusters from Spatial PCA, MFAI and VAE also capture spatial structure and expert labeling well, but with some notable exceptions, e.g., Spatial PCA cluster 7 in slices 4 and 9, (The Spatial PCA software performed additional post-procsssing on the clusters which why these clusters look less "noisy" than the others.) e eVAE, despite using the side information, did not seem to

---

rows due number of rows (n)
method
software
10a
10*
10s
10
ЕВMF
flashier 35]
0.8
2.5
36.9
163.1
сEBMF
"orch"
5.2
35.5
416.3
3,403.6
MFAI
mfair [17]
45.4
251.3
11,293.2
Spatial PCA
SpatialPCA [11]
234.8
8,213.7
Table 1: Running times of matrix factorization methods data sets in which varied n. the number of rows in X and Z. The numbers in the table are the average running times (in seconds) from 10 simulated data sets, *Available at https: //tpt/ggtbb//////// com/w
homm/wii wiwlill waiiimmmmmallt//aaaeaaeaa http://mmmmoorr///
improve over VAE. The CMF results were comparatively poor (Fig. 9 in the Appendix), reflecting the imappropriateness of the CMF assumptions in this setting. Note that the NMF methods can capture continuous variation in expression within and across cell types or regions as well as the expectation that some pixels might reflect combinations of cellular structures whereas the clustering cannot.
4.4 Scalability benchmark
cEBMF can also handle much larger data sets than the MovieLens and DLPFC data sets considered above. (One reason we did not use larger data sets was to allow for comparison with methods that do not scale well to large data sets.) To illustrate this, we ran EBMF and CEBMF on "tiled-clustering" data sets (using the same priors described Sec. 4.1), which the data sets were simulated with different numbers of rows, We compared the running times with two other matrix factorization methods that make use of the side information, MFAI and Spatial PCA (Table 1). While cEBMF had considerably higher running times than EBMF on the same data, cEBMF completed in much less time on average than MFAI and Spatial PCA. Further, while cEBMF was able to handle data sets with million rows, Spatial PCA struggled to analyze data sets with 100,000
to its high memory usage; for example. Spatial PCA needed approximately 300 GB of memory for 100,000. MFAI crashed frequently in data sets with n. 2 100,000 rows (only 2 of 10 of the completed at n 100,000). Note this benchmark was performed on a computer with 32 GB memory, NVIDIA GeForce RTX 4070 GPU and AMD Ryzen TM 9 7940HS CPU (8 cores, 16 threads). The EBMF and cEBMF algorithms were run for at most 20 iterations. See also the Appendix where we describe some of the computational properties of the cEBMF algorithmic framework.
Conclusions
We have introduced cEBMF, general and flexible framework for matrix factorization that (i) incorporates side information through flexible eevariate-deeeddee priors and (ii) learns these priors from the data using empirical Bayes ideas. Considerable effort has gone into optimizing the software implementation building on our previous work on this topic [13, 35]. As result, cEBMF scales well to large data sets with, say, hundreds of thousands millions of experiments highlight the importance of using matrix factorization models that make appropriate assumptions about the data performed competitively against other matrix factorization methods and deep learning approaches that make use of the side information. Because the priors in cEBMF can take the form of virtually any probabilistic model optimized via equations (10-11) I), our framework opens the door to incorporating other types of side information, including images and graphs.
and/or columns. Our
sufficiently flexible to adapt to the data. In experiments, cEBMF
Note: R and Python code implementing the experiments is available at bttps://githut com/ CBBMe_eepeeeeeieeeeeeeeee www.eeeppprmnt.e and Py Torch-bhaaas implementation of cEBMF is avail- able at at https: //gi http//gttttutttt////////tt/... thub. com/wil whs//wllllll iam denault t/ cebmf torch.
willi am- dww.ult/aeeeeeeeeeeeee CEBM/
Acknowledgments
We thank the staff the Research Computing Center at the University of Chicago for providing the high-performance computing resources used to implement the numerical experiments. We also thank
10

---

Deaglan Bartlett and Augustin Marignier for their help in developing the PyTorch implementation. This work was supported in part by NIH grant R01HG002585 (to M.S.) and Eric and Wendy Schmidt Al in Science Postdoctoral Fellowship. Schmidt Sciences, LLC program (to W.R.P.D.). Additional support came from the University of Chicago Data Science Institute through the 2024 AI+Science Research Initiative. We also thank the reviewers for their feedback.
References
[1] D. Lee and H. S. Seung. Learning the parts of objects by non-negative matrix factorization. Nature, (01(6555)778-79)
1999.
[2] D. D. Lee and H. S. Seung. Algorithms for non-negative matrix factorization. In Advances in Neural Information Processing Systems, volume 13, pages 556-562, 2001.
[3] N. Gillis. Nonnegative matrix factoriiaii...n.... Society for Industrial and Applied Mathematics, Philadelphia, PA, 2021.
[4] Sainburg, M. Thielk, and T. Q. Gentner. Finding, visualizing, and quantifying latent structure across diverse animal vocal repertoires. PLoS Computational Biology, 16(10)=e10022))) 1008228, 2020.
[5] T. A. Alexander, R. A. Irizarry, and H. C. Bravo. Capturing discrete latent structures: choose LDs over PCs. Biostatistics, 24(1):1-16, 2023.
[6] Novembre and M. Stephens. Interpreting principal component analyses of spatial population genetic variation. Nature Genetics, 40(5):646 649, 2008
[7] Leek and J. D. Storey. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genetics, 169), 2007.
[8] V. Marx. Method of the year: spatially resolved transcriptomics, Nature Methods, 18:9-14,
2021.
[9] K. Vandereyken, A. Sifrim, B. Thienpont, and T. Voet. Methods and applications for single-cell and spatial multi-omics. Nature Reviews Genetics, 24(8):494-515, 2023.
[10] B. Velten, M. Braunger, R. Argelaguet, D. Arnol, J. Wirbel, D. Bredikhin, Zeller, and
O. Stegle. Identifying temporal and spatial patterns of variation from multimodal data using MEFISTO. Nature Methods, 19:179-186, 2022.
[11] L. Shang and X. Zhou. Spatially aware dimension reduction for spatial transcriptomics. Nature Communications, 13:7203, 2022.
[12] W. Townes and B. E. Engelhardt. Nonnegative spatial factorization applied spatial genomics. Nature Methods, 20:229- 238, 2023.
[13] W. Wang and M. Stephens. Empirical Bayes matrix factorization. Journal of Machine Learning Research, 22(120): 104 2021.
[14] X. Zhong, C. Su, and Z. Fan. Empirical Bayes PCA in high dimensions. Journal of the Royal Statistical Society, Series B, 84(3):853887,, 2022.
[15] X. Lin and P. C. Boutros. Optimization and expansion of non-negative matrix factorization. BMC Bioinformatics, 21:7, 2020.
[16] C. H. Ding, T. Li, and M. Jordan. Convex and semi-nonnegative matrix factorizations. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(1):45-55, 2010.
[17] Z. Wang, F. Zhang, C. Zheng, X. Hu, M. Cai, and C. Yang. MFAI: scalable Bayesian matrix factorization approach leveraging auxiliary information. Journal of Computational and Graphical Statistics, (5(4):1339-1344 2024.
[18] P. Carbonetto, A. Sarkar, Z. Wang, and M. Stephens. Non-negative matrix factorization algorithms generally improve topic model fits. arXiv, 2105.13440, 2025.
[19] J. Lafferty and D. Blei. Correlated topic models. In Advances in Neural Information Processing Systems, volume 18, pages 147-154, 2005.
[20] M. Roberts, B. M. Stewart, and M. Airoldi, A model of text for experimentation in the social sciences. Journal of the American Statistical Association, 11 1116515):888-000 2016.

---

[21] A. P. Singh and G. J. Gordon. Relational learning via collective matrix factorization. In Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 650-658, 2008
[22] H. Lee and S. Choi. Group nonnegative matrix factorization for EEG classification. In Proceedings of the 12th International Conference on Artificial Intelligence and Statistics, pages 320-327, 2009.
[23] G. Bouchard, D. Yin, and S. Guo. Convex collective matrix factorization. In Proceedings of the 16th International Conference Artificial Intelligence and Statistics, pages 144-152, 2013.
[24] K. Gopalan, L. Charlin, and D. Blei. Content-bact based recommendations with Poisson factoriza- tion. In Advances in Neural Information Processing Systems, volume 27, 2014.
[25] C. Hu, P. Rai, and L. Carin. Non-negative matrix factorization for discrete data with hierarchical side-information. In Proceedings of the 19th International Conference on Artificial Intelligence
and Statistics, pages 1124-1132, 2016.
[26] D. P. Kingma and M. Welling. Auto-encoding variational Bayes. In 2nd International Confer- ence Learning Representations, 2014.
[27] K. Sohn, H. Lee, and X. Yan. Learning structured output representation using deep conditional generative models. In Advances in Neural Information Processing Systems, volume 28, 2015.
[28] X. Hc, L. Liao, H. Zhang, L. Nie, X. Hu, and T.-S. Chua. Neural collaborative filtering. In Proceedings of the 26th International Conference on World Wide Web, pages 173-182. International World Wide Web Conferences Steering Committee, 2017.
[29] Porteous, A. Asuncion, and M. Welling. Bayesian matrix factorization with side informa- tion and Dirichlet process mixtures. Proceedings of the 24th AAAI Conference on Artificial Intelligence, pages 563-568, 2010.
[30] J. H. Friedman. Greedy function approximation: a gradient boosting machine. Annals of Statistics, 29(5): 189- 2001.
[31] N. D. Lawrence and A. J. Moore. Hicrarchical Gaussian process latent variable models. In Proceedings of the 24th International Conference on Machine Learning, pages 481-488, 2007.
[32] Zhou, Shan, A. Banerjee, and G. Sapiro. Kernelized probabilistic matrix factorization: ex- ploiting graphs and side information. In Proceedings of the 2012 SIAM International Conference on Data Mining, pages 403 414, 2012.
[33] Y. Ma and X. Zhou. Accurate and efficient integrative referennce-ifformed spatial domain detection for spatial transcriptomics. Nature Methods, 21:1231-1244, 2024.
[34] D. Cai, X. He, J. Han, and T. S. Huang. Graph regularized nonnegative matrix factorization for data representation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(8): 1548-1560, 2011.
[35] J. Willwerscheid, P. Carbonetto, and M. Stephens. ebnm: an R package for solving the empirical Bayes normal means problem using a variety of prior families. Journal of Statistical Software 1 14(3))1-32, 2025.
[36] S. Nakajima and M. Sugiyama. Theoretical analysis of Bayesian matrix factorization. Journal of Machine Learning Research, 12(79):2583-2648 2011.
[37] S. Nakajima, M. Sugiyama, S. D. Babacan, and R. Tomioka. Global analytic solution of fully- observed variational Bayesian matrix factorization. Journal of Machine Learning Research, 14: 1-37, 2013.
[38] D. Yang, Z. Ma, and A. Buja. A sparse singular value decomposition method for high- dimensional data. Journal of Computational and Graphical Statistics, 23(4):923-944 2014.
[39] D. M. Witten, R. Tibshirani, and T. Hastie. A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. Biostatistics, 1033)511-533, 2009.
[40] B. Engelhardt and M. Stephens. Analysis of population structure: a unifying framework and novel methods based on sparse factor analysis. PLoS Geneties, 609)x1001111. 2010.
12

---

[41] Y. Liu, P. Carbonetto, J. Willwerscheid. S. A. Oakes, K. F. Macleod, and M. Stephens. Dissect- ing tumor transcriptional heterogeneity from single-cell RNA-seq data by generalized binary covariance decomposition. Nature Genetics, 37:263-273, 2025.
[42] I. Virshup, S. Rybakov, F. J. Theis, P. Angerer, and F. A. Wolf. anndata: access and store annotated data matrices. Journal of Open Source Software, 9(101):437,, 2024.
[43] R. P. Adams, G. E. Dahl, and I. Murray, Incorporating side information in probabilistic matrix factorization with Gaussian processes. In Proceedings of the 26th Conference on Uncertainty in Artificial Intelligence, 2010.
[44] N. Lawrence. Probabilistic non non-linear principal component analysis with Gaussian process latent variable models. Journal of Machine Learning Research, 6(60):1783-1811 2005.
[45] H.-M. Kim, B. K. Mallick, and C. C. Holmes. Analyzing nonstationary spatial data using piecewise Gaussian processes. Journal of the American Statistical Association, 100(470): 653-668, 2005.
[46] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Te- jani, S. Chilamkurthy, B. Steiner, L. Fang, Bai, and S. Chintala. PyTorch: an imperative style, high perfforerfe deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, H'Alché-Buc, E. Fox. and R. Gameet, editors, Advances in Neural Information Processing Systems, volume 32, 2019.
[47] Mitchell and J. Beauchamp. Bayesian variable selection in linear regression. Journal of the American Statistical Association, 13(404):1023-1033 1988.
[48] H. Robbins. Asymptotically subminimax solutions of compound statistical decision problems. In Proceedings of the Second Berkeley Symposium on Mathematical Statistics and Probability 1951, vol. II, pages 131-149. University of California Press, Berkeley and Los Angeles, CA, 1951.
[49] B. Efron and C. Morris. Limiting the risk of Bayes and empirical Bayes estimators Part II: the empirical Bayes case. Journal of the American Statistical Association, 67(337):130-139)) 1972.
[50] M. Stephens. False discovery rates: a new deal. Biostatistics, 18(2):275-22) 2244, 2017.
[51] D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: review for statisticians. Journal of the American Statistical Association, 112(518):859-8777 2017.
[52] F. M. Harper and J. A. Konstan. The MovieLens datasets: history and context. ACM Transactions on Interactive Intelligent Systems, 5(4):19, 2015.
[53] B. Pardo, A. Spangler, L. M. Weber, S. C. Page, S. C. Hicks, A. E. Jaffe, K. Martinowich, K. R. Maynard, and L. Tooladd-Toores.. spatialLIBD: an R/Bioconductor package to visualize spatially-resolved transcriptomics data. BMC Genomics, 23:434, 2022.
[54] K. R. Maynard, L. Collado-Torre L. M. Weber, C. Uytingco. B. K. Barry, S. R. Williams, J. L. Catallini, M. N. Tran, Z. Besich, M. Tippani, J. Chew, Y. Yin, J. E. Kleinman, T. M. Hyde, N. Rao, S. C. Hicks, K. Martinowich, and A. E. Jaffe. Transcriptome-scall spatial gene expression in the human dorsolateral prefrontal cortex. Nature Neuroscience, 24:425 436, 2021.
[55] M. Varrone, D. Tavernari, A. Santamaria-Martiner L. A. Walsh, and G. Ciriello. CellCharter reveals spatial cell niches associated with tissue remodeling and cell plasticity. Nature Genetics, 56:74-84, 2024.
[56] E. Zhao, M. R. Stone, X. Ren, J. Guenthoer, K. S. Smythe, T. Pulliam, S. R. Williams, C. R. Uytingco, S. E. B. Taylor, Nghiem, J. H. Bielas, and R. Gottardo. Spatial transcriptomies at subspot resolution with BayesSpace. Nature Biotechnology, 09:1375-1384 2021.
[57] J. Zhu, L. Shang, and X. Zhou. SRTsim: spatial pattern preserving simulations for spatially resolved transcriptomics. Genome Biology, 24:39, 2023.
[58] P. Carbonetto, K. Luo, A. Sarkar, A. Hung. K. Tayeb, S. Pott, and M. Stephens. GoM DE: interpreting structure sequence count data with differential expression analysis allowing for grades of membership. Genome Biology, 24:236, 2023.
13

---

[59] K. K. Dey, C. Hsiao, and M. Stephens. Visualizing the structure of RNA-seq expression data using grade of membership models. PLoS Geneties, 13()):e000659, 2017.
[60] Pons and M. Latapy. Computing communities in large networks using random walks. In Proceedings of the 20th International Symposium on Computer and Information Sciences, pages 2005.
[61] D. M. Blei. A. Y. Ng, and M. I. Jordan. Latent Dirichlet Allocation. Journal of Machine Learning Research, 3:993-1022, 2003.
[62] L.. K. Saul and M. Jordan. Exploiting tractable substructures in intractable netw orks. In D. S. Touretzky, M. C. Mozer, and M. E. Hasselmo, editors, Advances in Neural Information Processing Systems 8, pages 486-492. MIT Press, 1996.
[63] Z. Ghahramani and G. E. Hinton. Variational learning for switching state-space models. Neural Computation, 12(4):831-864 2000.
[64] M. A. van de Wiel, D. Te Beest, and M. M. Münch. Learning from a lot: empirical Bayes for high-dimensional model-based prediction. Scandinavian Journal of Statistics, 46(1):2-25, 2019.
[65] Y. Kim, Wang. Carbonetto, and M. Stephens. A flexible empirical Bayes approach to multiple linear regression and connections with penalized regression. Journal of Machine Learning Research, 25(185):1-59, 2024.
[66] Morgante, Carbonetto, G. Wang, Y. Zou, Sarkar, and M. Stephens. A flexible empirical Bayes approach to multivariate multiple regression, and its improved accuracy in predicting multi-tissue gene expression from genotypes. PLoS Genetics, 19(7)e1001039, 2023.
[67] S. Kullback and R. A. Leibler. On information and sufficiency. Annals of Mathematical Statistics, 22(1):79-866 1951.
[68] G. Wang. A. Sarkar, P. Carbonetto, and M. Stephens. A simple new approach to variable selection in regression, with application to genetic fine mapping. Journal of the Royal Statistical Society, Series B, 82(5):1233-1300, 2020.
[69] S. Wright. Coordinate descent algorithms. Mathematical Programming, 151:3-34 2015.
[70] Y. Kim, P. Carbonetto, M. Stephens, and M. Aniteseu. A fast algorithm for maximum likeli- hood estimation of mixture proportions using sequential quadratic programming. Journal of Computational and Graphical Statistics, 29(2):261-273, 2020.
[71] F. Chollet and others. Keras, 2015. URL
io.
[72] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, Moore, G. Murray, B. Steiner, P. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: system for large- scalescale machine learning. In Proceedings of the 12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16), pages 265-283, 2016.
14

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
The claims made in the abstract and introduction are supported by (i) a qualitative comparison of cHBMF to related work and (ii) empirical assessments in variety data sets.
Justification:
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We empirically assessed the performance cEBMF the situation where the prior misspecified.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "I "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and prool's
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
15

---

Answer: [Yes]
Justification:
The key theoretical results are found in the Appendix. We have provided
proof' of Proposition I and Lemma 1, and we have clearly stated the assumptions.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We have provided detailed descriptions of the experiments, including the software used. We have also provided the code that was used to generate all the results in the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g.. in the of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily
model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing verifying the results.
16

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
We have provided the R and Python code implementing the methods and experiments. The data sets were either included or instructions for accessing the data sets were given.
Justification:
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips cc/
for more details.
While we encourage the release of code and data, we understand that this might be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: / /nips.cc/ppll wwtppppwwccgg.eeeeeeeeecoeeepeee http://gggiddee/eeeuububiiiiiii www.eesuummsssoooooo..oooo
for more details.
The authors should provide instructions data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: Many of these details are given in the Appendix. Additional implementation details can
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: We summarize some results using boxplots with the conventional definitions for whiskers, box bounds, center line and outliers.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by emror bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
17

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out range (e.g. negative еπτог rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification: Details of the computing hardware used in the scability benchmark are given in Sec. 4.4.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: //neur:////uriss....
Answer: [Yes]
Justification: We have reviewed the NeurIPS Code of Ethics and work abides by it the best of our knowledge.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., s if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification:
We don't believe that our work has any obvious direct social impacts.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
18

---

If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), faimess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: This work did not release data or models that would be considered to have strong potential for misuse.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging. and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
All the original data sources are credited/cited.
Guidelines:
The answer NA means that the paper does existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
19

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We provide the source code for new methods, and it is accompanied with detailed documentation.
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/moddl as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
This work does not involve crowdsourcing or research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: this work does not involve research with human subjects.
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
20

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: LLMs were not were not used in this work.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https:/ /ttpar/////u.....//e /c/Coonferencee//22551/.. for what should or should be described.
21

---

Appendices
(1) A Derivations and additional method details
The cEBMF algorithm (Sec. 3.3.2) is variational empirical Bayes algorithm [61-64] that is formulated as solving the following optimization problem:
argmax ELBO(g.9,.))
(16)
9.9.7
where g is shorthand for the priors (6)
9K q is distribution (L, F), and ELBO 9, T) is the "Evidence Lower BOund" (ELBO) [51], lower bound to the "evidence", log p(Z 9, T):
9K 96 (1)
(
ELBO (g == E., Doggp( L, F.]] + Eg log pLL, ((L,
FU,,FXXYY X. XX
}
(17)
See [13, 65, 66] for other variational empirical Bayes algorithms derived in a similar way.
To achieve tractable update expressions for the model parameters, we approximate the posterior (L, F) so that factorizes over all elements and elements (sometimes called "mean field approximation):
(L, F) (L)
(L) II (la)
(18)
K
of
(F)
j=
With this factorization (or conditional independence) qconstraint on the right-hand part of the ELBO can be immediately decomposed into of expectations over the individual elements of L and F, so we have
ELBO T) = Logg((Z L, F..)
Ih (1) (1i:
}
Eg
log
Che
k=1
9k ( (fsi vi)
Eg
log
}
(19)
(E) where 9% (e; denotes the density of gb (xx) at l, and of (f: (f: denotes the density of gk (x:))
(
at f.
A.1 Updating the factors
The following proposition formally connects the updates of the individual factors h 1...., ., (Step
2-4 of the algorithm in Sec. 3.3.2) to learning a covariate-moderrate EBNM model (Sec. 3.3.1).
Proposition 1. Let & = "(kk)) denote the kth column of L, let fx (fiA,... Fpk ) denote the kth column of F, let Cs [x), fx [fk] L further define
E, [ and f? = Eylff]] and we
(ell)) = (E()
(20)
P II
(f(
of (fx)).
(21)
1=1
22

---

Let RRt denote the n x matrix of expected residuals (with elements ( that ignores the contribution of the kth factor,
F = Zij
Σ
Taw Flw..
(22)
Also define i t, w, T). fzz t, w, T), s(w, 7) and of (w, T) as vector-valued functions which the individual vector elements given by
Ei=1 Tajzyزي/
(Z L, T)
(23)
[s T)
Di=1 Tiz
(Z, w,
(24)
[8] (w,
(ww T) (E)=1 Tiy w 1122
(25)
of (WWW)) T) (E wi) 1/2
(26)
Then using the definition of the ELBO in (17) and the cEBNM mapping defined in (13), have that
argmax, 19
4..9 3. E
ELBO
CEBNM(E(R* fr, 8 se T), X, Ge,,)
(27)
argmax, ( EGG.
FILBO (g.9, +) CFBNM(f(R* La Ey 7), (Cg T) Y (y
(28)
Note that this identity requires slight change the definition of the cEBNM mapping (13) as returning the priors g(d,, Q) rather than the parameter estimates f themselves.
(6) Proof. Starting from (19), we expand on the parts of ELBO that involve of or 9 or both:
1
9k (0) ( x.)
}
ELBO(q g. r)
Σ Eg lautt 2b:4a] +
Σ E
log
+ const,
(29)
2 1=1
(lak)
2m1
lik
where "const" is placeholder for the terms in the ELBO that do not depend on the kth factor, and we define
aik =
(30)
==1
du Σ Toyyly S
(31)
The identity (27) then follows from Lemma similarly.
(given below). The other identity (28) is proved
A.2 Updating the residual variances
Focussing on the part of the ELBO depends on т, we have
ELBO(g. g) = 2 1=1
(log TTy + const,
(32)
in which "const" is a placeholder for the other terms in the ELBO that do not involve т, and F?, is the expected squared difference between the observation zi and the value predicted by the matrix factorization:
P2 := [,i -2,)
(33)
where
K
2ig
Lik fike
(34)
k=1
23

---

If one makes the modeling assumption that all the residual variances are the same, i.c, Tij = then from (32) the update for works out
11 X
p
(35)
2 i= Σ マリ=
If instead one makes the weaker modeling assumption that the residual variances are the same in each column, i.e., Tij = J = ,,...,P, then the updates work out
(36)
Similarly, for row-specific residual variances the updates are
P
Ti
(37)
For all these expressions, the squared differences Fij casily computed given the conditional independence assumptions of the fully-faaee fallo--fe approximation (18):
K
2
K
K
Zij Σ
(W) - Σ (Tandgh))
(38)
R=1
k=1
in which we have defined le == Byllet), JSx= == Byl. = = and == -.,
A.3
Covariate-moderateed EBNM
To complete the proof of Proposition 1, remains to show that the identity (27) satisfied at the objective function given by (29). (And similarly for the identity (28).) This connection is made in the following lemma.
Lemma 1. Consider the cEBNM mapping defined in (13). An equivalent definition of this mapping
(0.0) () argmaxe, F(O, 12 3, D),
(39)
where
F(O, q: D) = 1
g(B de 0)
E Ta. 2b,33] +
Eyg
log
(40)
(:(3s)
1=1
1=1
and g(3; d: 0) denotes the density of g(d; 0) at B, and we further define ((3) = I (:()))
(41)
0; = 1/s²
(42)
b. 8//87-
(43)
Proof. We begin with the ELBO for the cEBNM model (7, 9):
ELBO(6, q; 3, D) log L()) - Dxt.(g Ppos).
(44)
where (g p) denotes the Kullback Leibler (K-L) divergence from a distribution to distri- bution [67], L(0) is the marginal likelihood defined in (11), and Ppea) is the (exact) posterior
distribution, Ppost(()) := [T p(3; B, 8;, O, (see eq. 12). Since DKL(q p)) is always zero or greater, and is exactly zero when p q, we have that argmax, ELBO(0, B, s, D) = Ppest and max, ELBO((, g; B, s, D) log L(0). Next, basic identity of the ELBO (see for example Appendix B of [68]) is that the ELBO can be rewritten as
Egg log g(365 0)
H
mathrmEL.BOO q: 3, D) E.g (log B, s)] +
(45)
(1,)))
To complete the proof, expand terms in the likeliioldood in (45):
(B. 3,) + const,
log 1( [3,8) 2
(46)
where "const" a placeholder for terms that do not involve q (or g). Plugging this identity into (45), and with bit of additional algebraic manipulation, we recover (40).
24

---

Algorithm cEBMF algorithm
Require: n x p data matrix, Z; covariate "side information" matrices, X (n n and Y ny): K, the number of factors; the prior families Ge,k and Gf... k = l...., K; and initial estimates of the first and second moments of (n x K), F (p x K), which are denoted by L, F. L F2. Compute the expected residuals, R Z LF repeat
Update the residual variances using (35), (36) or (37).
for k = . K do
Remove the effect of the kth factor from the expected residuals, R$ R + L
Perform a single-factor update for factor l (Algorithm 2).
Update the expected residuals, R R fxJk ft end for until some convergence criterion is met return F, ..... 9K 91 (f)
9K
Algorithm 2 cEBMF single-factor factor update
Require: covariate or "side information" matrices, X (n x ng) and Y (p x ny); E (K)...,K), the dimension update; the prior families Gz, and Gy,s5 implementation of cEBNM(3,e 8, D) + (0, q) (eq. 13) for prior families G Gr and G Gysi the expected residuals, R*; estimates of the second moments, L (n K). F2 (p x K); and the residuals
variances, T.
1. B + [(R f2.
(f 2. 8 + 8
3. (gk ( CEBNM(3, s, X, 9c.k)
4. Compute posterior moments la == Eyllik] and e3 Eyl I....,n.
5. β + f(Rk, Ee,l,, T)
6. 8 + 8 (L
7. (
aaf) CEBNM( 8, Y.
8. Compute posterior moments fyx Ey and F3. == Eglffah. 1.....p.
9. return L, Eg ft. fhe 9k 9
(1) (f)
A.4 Detailed algorithms
In summary, the cEBMF algorithm is a block co-ordinate ascent algorithm [69] for finding local maximum of the ELBO (17), in which the "blocks" i.e., the subsets of parameters to be updated are the individual factors = I..... K (Sec. A. and the residual variances T (Sec. A.2). This co-ordinate ascent algorithm is described in Algorithm 1, and the single-factor update is described in Algorithm (And it is described informally in Sec. 3.3.2.) In practice, we run Algortithm until the increase in the ELBO across two successive iterations is smaller than some specified tolerance, until we have reached an upper bound on the number of iterations.
Two features of the empirical Bayes approach to matrix factorization discussed in [13] are worth highlighting here. First, there a simple stepwise procedure for obtaining good initial estimates of L and F by introducing the factors sequentially. This was called "greedy initialization" in [13]. Second, instead fixing the number factors, the EBMF approach can also select K automatically
(8) by adapting the priors 9x. . separately for each factor k. The idea is that factors that are not useful for explaining variation in the data should produce priors that are concentrated near zero (this feature of course requires that the chosen prior families Gr,k, Gy include distributions that are concentrated near zero). Therefore, can initially set to a large value, and the cEBMF algorithm will automatically determine appropriate number of factors by "shrinking" the unneeded factors.
(f)
A.4.1 Computational complexity
Since cEBMF a modeling and algorithmic framework, and not a specific method or algorithm, we cannot give the exact computational complexity of Algorithm 1. However, provide some rules
25

---

of thumb. Steps 1, 2, 5 and 6 in Algorithm 2 (also, Steps 1 and 2 in Sec. 3.3.2) involve preparing the inputs for the CEBNM solver (13). Since these steps do not depend on the prior families Gr.k J we can give their computational complexity: when Z is "dense" (non-sparse) matrix, the time complexity for updating single factor is O(np); when Z is sparse, the complexity is O(S), where S is the number of nonzero entries in Z. (Note this requires careful implementation that avoids directly storing R). Steps 3 and 7 in Algorithm (or Steps 3 and 4 in Sec. 3.3.2) will depend on the details of the cEBNM solver and the type of side information. However, when the priors on L, F are simple and involve low-dimensional covariates, the other steps are expected to dominate, in which case the complexity of Algorithm 1 is expected to be (пpp)) or O(SK).
B Details of the experiments
B.1 Simulations
We simulated data sets from different cEBMF models. In all cases, the data were generated with homoskedastic noise, Taj T.
Sparsity-driven covariate simulations. This simulation was intended to illustrate the behaviour of cEBMF when provided with simple row and column-covariates that inform only the sparsity of and F (and not the magnitudes of their elements). The side information was stored in 1,000 x 10 and 200 x 10 matrices and Y, and the 1,000 x 200 matrix Z was simulated using a simple cEBMF model with K 2 and with anike-and-slab priors chosen to ensure that 90% of the elements of LF" were zero. Specifically, we used the following priors:
lak ~ Tu + (1 - T.s. (0,,
fox= au + (1 (0,1)
(47)
TiA: (0 xx;)
Ojk == N(w yı),
in which 0g and w/ were chosen to achieve 90% zeros in LFT.
Uninformative covariate simulations. To verify that cEBMF was robust to situations in which the side information was not helpful, we considered an "uninformative covariate" setting in which the covariates were just noise. The data sets were simulated in the same way as in the sparsity-driven covariate simulations except that the true factors were simulated as lek Tōg + (1 T) N), 1). fix adn + (1 (O, 1), with n, chosen to achieve a target sparsity of 90% zeros in LF"..
Tiled-clustering simulations. In this setting, we simulated rank-3 matrix factorizations in which L but not F depended the 2-d locations of the data points. (One of these simulations is shown in Fig. 1.) This was accomplished as follows. First, we generated a periodic tiling of [0, 1] x [o, 1]. randomly labeling each tile 1, 2 3. For each data point
I,..., 1,000, we sampled its 2-d location uniformly from [0, 1] x [o, 1.]. then we set Esp- = 1 if the data point fell in the tile with label k. otherwise L 0. The 200 3 matrix F was simulated from a scale mixture of zero-centered normals that did not depend on tile membership.
Shifted tiled-clustering simulations. To assess robustness to model misspecification, we simulated data using prior that could not be recovered by the prior family used in cEBMF. These simulations were the same as the tiled-clusterin clustering simulations except that we generated the ith row of L as follows: (1, 2, 3) if data point and (2, 3, 1) if data point was in the tile with label 3.
in the tile with label 1; (3, 1,2) if data point was in the tile with label 2,
B.2 Additional details on the methods compared
We first describe how the methods were run on the simulated data sets. Modifications to the methods for the MovieLens and spatial transcriptomics data are given in the main text, with additional technical details below. For all methods, when possible to do so, we set the rank, K, to match the rank of the simulated matrix factorization.
26

---

For EBMF and cEBMF, we assumed homoskedastic noise (Ti = =), and the prior families were chosen to align with how the data were simulated (except for the shifted-tiled clustering simulations, which were intended to illustrate the methods" behaviour when the priors were misspecified). For EBMF, the prior families were all elaborations of the "spike and slab" priors (Sec. 3.2.1). For cEBMF, the priors were of the same form as EBMF in which the mixture weights were parameterized using either a multinomial regression (i.e., a single-layer neural network with a softmax link function) or a multilayer perceptron.
In the pparsity-drives covariate and uninformative covariate simulations, the priors for L and F in EBMF were all scale mixtures of normals with fixed grid of scales [50, 70J. For cEBMF, we used priors of the same form, except that side information was incorporated into the prior mixture weights as follows using priors of the following form:
g(d,, 0) To 0)0o +
Σ Tm (d,, (0, om)) m==
(48)
The mixture weights To..... TM were implemented using a standard multinomial regression model with the softmax link function.
In the tiled-clustering and shifted tiled-clustering simulations, the true L was always non-negative. Therefore, we chose the prior families in EBMF and CEBMF to produce semi-non-negaaive matrix factorizations [16] with non-negative L. Specifically, we assigned mixture-of cexxture-of-epoonee of-exponenssals priors to L, similar to the only [35]. And we assigned the seale-mixtre-oo-norre priors, same as above, to F. In cEBMF, side information was incorporated into the mixture weights in the prior in a manner similar above:
m-w-tmre-o--ooommaeeee priors, except with support for non-negative numbers
M
g(d,, 0) = to (d:, 0700
+
Σ Tm O) exp(Am), m=1
(49)
in which exp(A) denotes the exponential distribution with scale parameter X, and Am-1 < Am-
m = .,, M. As before, the mixture weights Tu. TM were implemented using standard multinomial regression model with the softmax link function.
The deep learning methods (VAE, VAE, NCF) were all implemented in PyTorch. All the models were trained for 50 epochs using the Adam optimizer with learning rate 0.001 and batch size 64. VAE had three hidden layers (of width 128, 64 and 30) in both the encoder and decoder (20 hidden dimensions). ReI.U activations were used throughout. We use the ELBO from [26] to train the model. We proceeded similarly for the cVAE, conditioning both the encoder and decoder the available covariate data X, Y. For c cVAE, we used the training objective from [27]. NCF models Z using two separate multilayer perceptrons for the row and column covariates [28]. The multilayer perceptrons were implemented in similar way to the VAE encoders and decoders; that is, three hidden layers (of width 128, 64 and 30) with RELU activations. The penalty parameters in PMD were tuned via cross-ve valss-valiottion as recommended by the authors. SSVD (R package "ssvd") was run with its default values.
For the spatial transcriptomies data, we fit cEBMF and EBMF using gene-specific residual variances, oîj =o We used mixture-of-ex -fx-o-nentxppottialllll priors for F, and the parameterized cixture-of-expooentti priors (49) for in which the mixture weights were learned using a multilayer perceptron instead of multinomial regression. The multilayer perceptions were defined as sequential models with dense layer with 64 units and ReLU activations. We use two subsequent dense layers, each with 64 units, and ReLU activations using an L2 regularization coefficient of 0.001 to prevent overfitting. These regularized layers were followed by dropout layer (with a dropout rate of 0.5). The subsequent layers were four dense layers each with 64 units, ReLU activations and L2 regularization coefficient of 0.001. The final layer was a dense layer with a softmax activation. These models were trained during each single-factar factor update using 300 epochs and a batch size of 1,500.
In the simulations, cEBMF was implemented in R, which learning the parameterized priors was performed using the Keras R interface [71] to TensorFlow [72]. For the Movielens and spatial transcriptomics data sets, we used the PyTorch-based implementation of cEBMF which we have made available as a Python package on GitHub.
27

---

C
Additional results
Unkiaratt 09.00 +01
Figure 5: Simulation results from the "uninformative covariate" setting in which the data were simulated under different noise levels, T. Note that for improved visualization the RMSE shown on the log-scale.
0.0
m0
Figure 6: Simulation results from the "sparsity-driven covariate" setting in which the data were simulated under different noise levels, Note that for improved visualization the RMSE shown on the scale.
28

---

83
30 248
289
1050
1800 OFF 080 1RY
W/ 12 20PEN3200 AV.
10:003
335033
31.
Figure 7: Simulation results from the "tiled-clustering" setting in which the data were simulated under different noise levels, Note that for improved visualization the RMSE is shown the log-scale. (spaPCA = Spatial PCA)
Shilled tied-clustering model
noise sd
noise sd
2
3.0
30-
1.0
0.9-
RMSE
CEBMF
EBMF
MFN
PMD
8810
BVD
6BBMF
EBMF
MFW
FMD
SBVD
8VD
Figure Simulation results from the "shifted tiled-clustering" setting in which the data were simulated under different noise levels, Note that for improved visualization the RMSE is shown on the log-scale.
29

---

numan labe ed
GEBNMF
CMF
MFAI clustering
Sperial PCA clustering
13
14
A5
18
numan Isbered
GEBNMF
CMF
MFAI clustering
Spotial PCA clustering
Figure 9: Additional results on slides 4 (top) and 10 (bottom) of the DLPFC spatial transcriptomics data. See Fig. for additional information about these results.
30

---

EBNMF
human labeled
GEBNMF
ESNMF
humen labeled
Spetial PCA
CEBNMF
EBNMF
NMF
K2
K3
K1
12
64
K5
K7
16
Spatial
EBNMF
human Iabeled
Spatial PCA
cEBNMF
EBNMF
NMF
Figure 10: Selected results on slices (top row) through 6 (bottom row) of the DLPFC spatial transcriptomics data.
31

---

human labeled
Spatial PCA
CEBNMF
EBNMF
NMF
human Isbeled
Spatial PCA
GEBNMF
EBNMF
NMF
GEBNMF
EBNMF
labeled
GEBNMF
ESNMF
human lsbeled
Spatial PCA
cEBNMF
EBNMF
Figure 11: Selected results on slices (top row) through 12 (bottom row) of the DLPFC spatial transcriptomics data.
32


Reviews:
[{"review_id": "wk1exoJGIf", "reviewer": "Reviewer_NLUp", "summary": "The paper proposes cEBMF, a modular empirical‑Bayes matrix‑factorization framework that can ingest side information (text, images, graphs, coordinates) via covariate‑dependent prior families.\nThe framework unifies and generalizes prior methods (EBMF, NMF, Spatial PCA, MFAI, CMF) as special cases through flexible prior families and automatic hyper‑parameter learning.\nThe authors derive a simple coordinate‑ascent algorithm that reduces the learning problem to repeated covariate‑moderated empirical‑Bayes normal‑means sub‑problems optimizing a variational approximation, enabling plug‑and‑play priors and scalability.\nIn experiments, state‑of‑the‑art accuracy is shown across simulations, collaborative filtering (MovieLens 100K) and spatial transcriptomics (human DLPFC), while outperforming EBMF, MFAI, Spatial PCA, PMD, VAE/cVAE, NCF.\nAn open‑source R/Python implementation is provided including runtime benchmarks where cEBMF scales markedly better than covariate-dependent alternatives MFAI and Spatial PCA.", "questions": "How sensitive is the performance to the choice of covariate model?\nWhat guarantees (if any) do you have about convergence or posterior consistency in this setting?\nIs your method applicable in settings where covariates are high-dimensional or partially missing?", "limitations": "No Limitations Section is provided. Arguably, some limitations are implicitly addressed, by describing the limited scope of the methods.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "(+) generalizes many existing matrix factorization methods as special cases\n(+) a computationally efficient implementation is provided\n(+) The paper is written with generality in mind\n(+) The experimental evaluation is extensive\n(-) due to generality it was somewhat work for me as a reader to figure out, which priors and models have actually been used. Sometimes instantiating special cases helps clarity.\n(-) compared to EBMF, the covariate-dependent cEBMF seems like an obvious extension of known pieces. One might even argue that it is a special case of EBMF where the prior family is parametrized using covariates.\n(-) no theoretical results are provided beyond derivations of the algorithm and methods.", "was_revised": true, "final_justification": "Thank you for the response and clarifications. In the light of the comments, I deem my original assessment appropriate.", "rebuttal": "We thank the reviewer for the time and the detailed review of our work.\n\n**How sensitive is the performance to the choice of covariate model**? \n\nRegarding your question regarding the choice of the covariate model, we perform a simulation study in which the covariate model cannot recover the data-generative model (see shifted tiled clustering, section 4). We show that the results are very similar to those of the EBMF, and thus using a completely misspecified covariate model has a rather limited impact on the cEBMF performance. See Figure 2 on the right panel.\n\n**What guarantees (if any) do you have about convergence or posterior consistency in this setting**?\n\nSimilar to most fully factorized variational approximations, we can only guarantee that the fitting procedure will converge to a fixed point that is a local maximum of the ELBO, but we don't have posterior consistency results in this setting. Future work may focus on getting more precise on posterior consistency.\n\n**Is your method applicable in settings where covariates are high-dimensional or partially missing**?\n\nOur approach is applicable when the covariates are high-dimensional. See our sparsity-inducing simulations, where we have 1000 covariates for factors of length 200, and among the 1000 covariates, only three are actually affecting the factor sparsity. We show that cEBMF outperforms the other methods in this setting.\n\nRegarding missing values, this depends completely on the covariate model selected by the user; some have \"built-in\" missing value routines. In our paper, we use neural nets that modulate the prior weight of different mixture components.  However, it is possible to use XGBoost instead of a neural network, and XGBoost has some missing data routines already built-in. In this sense, cEBMF can be viewed as a \"wrapper\" that can benefit from the covariate model perks.\n\nWe hope that we have answered your questions and that you may consider increasing our rating."}, {"review_id": "wF2F3doJK0", "reviewer": "Reviewer_XAok", "summary": "This paper introduces cEBMF, a framework that integrates side-information into Empirical Bayes Matrix Factorization by learning both spike-and-slab mixing weights and slab parameters as differentiable functions of covariates.  The authors demonstrate improved reconstruction accuracy over EBMF, MFAI, Spatial PCA, and VAE-based baselines on simulated data, MovieLens 100K, and spatial transcriptomics, and show scalability to up to one million rows.", "questions": "* How about the robustness to prior mis-specification? I suggest to include an ablation showing performance when the true slab distribution differs from the chosen parametric family.\n\n* The paper instantiates the covariate-modulated prior using a specific parametric family. Can the authors characterize the representational capacity of that family—i.e. what classes of true slab distributions it can approximate well? Moreover, what would it take to replace this with a fully nonparametric density estimator (e.g. a normalizing flow or kernel density network)? I recognize that this adds optimization and identifiability challenges, but such an analysis could clarify the upper bound on prior flexibility and guide future nonparametric extensions.", "limitations": "The author does not describe the limitation in the main part of the paper but mentioned several in the checklist.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths\n\n* New Prior Modeling: developed a new covariate→prior-parameter mappings in an EBMF framework;\n\n* Modular Learning  Structure: Decomposes the overall inference into repeated calls to a cEBNM solver, enabling plug-and-play of different solver backends.\n\n* Empirical Validation across Domains: Provides thorough experiments on synthetic scenarios, recommendation data, and spatial omics, demonstrating consistent gains when side information is informative, and no degradation when it is not.\n\n## Weakness \n\n* By design, cEBMF cannot capture nonlinear interactions or higher-order effects; no comparison is made to deep nonlinear baselines (e.g. self-attention recommenders).\n\n* Omits comparisons to modern graph-based or attention-based methods (e.g. GC-MC, LightGCN, SASRec, BERT4Rec), which also integrate side information.\n\n* Requires pre-specifying a family for g_1 (Gaussian, Exponential, mixtures). When the chosen family mismatches data, performance may suffer; robustness to mis-specification is not studied.\n\n* No formal analysis of identifiability, convergence rates, or generalization bounds—particularly important when covariates weakly correlate with latent factors.", "was_revised": true, "final_justification": "The authors’ responses have addressed some of my concerns. However, the lack of comparison with methods from other families remains unresolved. Therefore, I will keep my original rating of 4", "rebuttal": "We thank the reviewer for the time and the detailed review of our work.\n\n\n**How the  method handles prior mis-specification?**\n\nThis is an important point raised by the reviewer. We would like to highlight that we actually did what the reviewer is wondering. We performed a simulation study in which the covariate model cannot recover the data-generative model (shifted tiled clustering). We show that the results are very similar to the EBMF ones and so using a completely misspecified covariate model has rather limited impact on the cEBMF performance. See Figure 2 most right panel\nRegarding your questions about the choice of prior family and the use of mixture priors, the \"spike-and-slab\" prior presented in Sec. 3.2.1 is intended only to illustrate the key ideas. (Although, to be clear, the \"spike-and-slab\" prior is an example of a mixture prior—it is a mixture prior with two mixture components.)\n\n**Other questions**\n\n For the simulations and spatial transcriptomics data application, we used slightly more flexible mixture priors; the details are given in Appendix C of the paper. Certainly, other priors could be used and may be more appropriate for different data sets or for different types of side information. We have not considered the use of hierarchical Dirichlet models nor multi-kernel Gaussian processes, and it is hard to comment on the benefits of such approaches without being more specific about the type of data or side information (categorical, spatial, temporal, etc) being considered. It is straightforward to use cEBMF with hierarchical Dirichlet models or multi-kernel Gaussian processes prior, given that you can solve the cEBNM mapping for these priors. So, more generally, cEBMF can be used with any prior model in which the cEBNM mapping can be computed or at least numerically approximated (see Equation 13).\n\n**Limitations**\n\nRegarding the limitations acknowledgment, while we agree that we don't have a specific section about it, we spent a great deal of time describing the limited scope of the method in the first section of our manuscript.\n\n\nWe hope that we have answered your questions and that you may consider increasing our rating."}, {"review_id": "MOrEUgN8st", "reviewer": "Reviewer_cne9", "summary": "This paper introduces a covariate-moderated empirical Bayes matrix factorization framework (cEBMF) for incorporating side information into matrix factorization. The authors first extend the EBMF model by allowing the priors over latent factor entries to depend on covariates (side information), enabling the model to adaptively learn sparsity or structure guided by auxiliary data. Then, a modular variational algorithm is proposed, where matrix factorization is learned by solving covariate-modulated empirical Bayes normal means (cEBNM) subproblems iteratively. The framework supports various prior families, including spike-and-slab for sparsity or non-negativity. Experiments on simulated data, MovieLens 100K data, and spatial transcriptomics data show that cEBMF achieves superior performance and scalability compared to existing methods like EBMF, MFAI, Spatial PCA, and deep learning baselines.", "questions": "i) While the authors included comparisons between cEBMF and EBMF, would it be possible to include an ablation study exploring different prior families? Additional visualizations showing performance variation under different priors would help.\nii) How robust is cEBMF to noisy or misspecified covariates? A sensitivity analysis or simulation study on covariate quality would be valuable.\niii) Have the authors considered methods for automatically inferring the number of latent components (K)? For example, a Bayesian nonparametric approach using a Dirichlet Process Mixture or similar mechanism might offer a flexible alternative.", "limitations": "Yes.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths\ni) The paper presents a well-founded theoretical framework with thorough derivations and proofs. The proposed algorithms are clearly described and easy to follow.\nii) Visualizations are informative, with multiple figures that effectively illustrate the experimental results. Applications to MovieLens 100K and spatial transcriptomics are particularly compelling.\niii) The incorporation of covariates (side information) into the factor model is innovative, and the extension of EBMF is both well-motivated and systematically developed.\n\nWeaknesses\ni) The paper lacks an ablation study evaluating the impact of covariates or the choice of different prior families on model performance.\ni)) The distinction between “EBMF” and “EBNMF” is somewhat unclear and may be confusing for readers unfamiliar with the naming convention.", "was_revised": true, "final_justification": "I am satisfied about authors' response during the rebuttal period and appreciate their effort in trying to address my concerns in a short period of time, especially to my last point. \nOriginally I already leaned to accept their paper so I would keep my rating.", "rebuttal": "We thank the reviewer for the time and the detailed review of our work.\n\n**The reviewer asks if it would be possible to include an ablation study exploring different prior families**? \n\nRegarding your question on the choice of the prior model, we perform a simulation study in which the prior model cannot recover the data-generative model (see shifted tiled clustering, section 4). We show that the results are very similar to the EBMF ones, and so using a completely misspecified covariate model has rather limited impact on the cEBMF performance.  See Figure 2 most right panel\n \nFurther, the reviewer suggests that **\"A sensitivity analysis or simulation study on covariate quality would be valuable.\"**\n\nWe included a simulation study in which we use uninformative covariates (see Figure 2 most left panel). When the covariates are uninformative, cEBMF has similar performance to EBMF, which is a sign of the good behaviour of the fitting procedure, as the EBMF model corresponds to a cEBMF model without covariates.\n\nLastly, the reviewer asks: **\"Have the authors considered methods for automatically inferring the number of latent components (K)? For example, a Bayesian nonparametric approach using a Dirichlet Process Mixture or similar mechanism might offer a flexible alternative.\"**\n \nThis is a very good and important point. Due to the page limit, we discuss the choice of K in the appendix (A.4). We do not consider the Dirichlet Process Mixture but rather a more classical and simpler approach from the Empirical Bayes literature.\nWe write in appendix (A.4)\n\n\"A noteworthy aspect of empirical Bayes methods for matrix factorization,\nhighlighted in Bishop 1999, is their inherent\nability to automatically determine $K$ . This\nfeature comes from the fact that the maximum-likelihood estimates for\n$g_k^{(\\ell)}$ and $g_k^{(f)}$ can converge to a point-mass at zero\n(assuming that the families $G_{\\ell, k}$ and\n$G_{f, k}$ include a point-mass at zero). Thus, if $K$ is\ninitially set to a large value, certain row or column factors will\nconverge to zero, which effectively \"zeroes out\" the corresponding\ncomponent $k$.\n\nSo as long as the prior model family contains  the point mass at 0, then cEBMF automatically estimates K This means that if K is set sufficiently large, then some loading/factor combinations will be optimized to be exactly 0.\n\nWe hope that we have answered your questions and that you may consider increasing our rating."}, {"review_id": "55kJdiNvOK", "reviewer": "Reviewer_6qut", "summary": "The paper introduces cEBMF, a flexible and modular framework for matrix factorization. The method integrates various types of \"side information\" by incorporating them into adaptive, covariate-dependent prior distributions. cEBMF learns these priors directly from the data, enhancing the inference and summarization of multivariate data structures. The benefits of cEBMF are demonstrated through simulations and applications to spatial transcriptomics and MovieLens datasets.", "questions": "### Choice of $K$\n\nThe cEBMF framework demonstrates impressive flexibility in incorporating various prior assumptions and types of side information. However, the rank $K$ is treated as a fixed hyperparameter in the current experimental setup.\n  \nGiven the modular nature of cEBMF and its stated ability to accommodate different prior families, have the authors considered extending the framework with non-parametric Bayesian approaches, such as Dirichlet Process priors, to allow for automatic inference of the optimal rank $K$ ? What are the theoretical or practical challenges in implementing such an extension, and are there plans for future work in this direction?\n\n### Limited sample scenario with large $K$\n\nA common challenge in matrix factorization models is dealing with scenarios where the number of latent factors, $K$, is relatively large compared to the available data size.\n  \nWhat is the expected behavior and potential limitations of cEBMF in such \"small data, large $K$\" regimes? Specifically, since training the function $\\pi(\\cdot, \\cdot)$ (which determines the covariate-dependent prior probabilities) in a \"small data size, large K\" regime seems to present significant challenges despite its flexibility, how might the empirical Bayes nature and the covariate moderation mechanisms of cEBMF help to mitigate issues like overfitting or unstable parameter estimation in these challenging settings? \n\n### Optimization\n\nThe cEBMF learning algorithm relies on an iterative coordinate ascent approach to maximize the ELBO. The performance and convergence of such optimization procedures can often be sensitive to the initial values of the latent factors and prior parameters.\n  \nCould the authors provide more details on the initialization strategies employed for the latent matrices ($L$ and $F$) and the parameters of the covariate-moderated priors ($\\theta$​)? How do these choices impact the final results and convergence speed, particularly in the challenging scenarios presented in the experiments? Does the framework incorporate any specific robust initialization techniques or strategies to mitigate potential issues arising from poor initialization?", "limitations": "no issue", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n\n- cEMBF's modular and flexible design is attractive and the novelty of incorporating any type of side information is practical. In addition, the proposed algorithm tries to learn the covariate-dependent prior distribution representing the side information is inspiring.\n\n- The cEMBF framework generalizes different matrix factorization algorithms with different regularizations under the view of covariate-dependent priors. It is simple yet inspiring to see well-established algorithms like sparse NMF can also be viewed under this framework. The scalability of the framework is also promising. \n\n### Weaknesses\n\n- The performance of cEMBF on small dataset, large $K$ regime is unclear, yet it is a very common scenario in real application\n- The choice of the latent dimension is $K$ for all experimental setup, while given the flexibility claimed in cEMBF, one should be able to incorporate Dirichlet process prior so $K$ is learnable.\n- While the performance of many ELBO optimization procedures can indeed be sensitive to the choice of initialization, the provided paper content does not explicitly detail the specific initialization strategies used for the latent matrices (L and F) or the prior parameters $\\theta$.", "was_revised": true, "final_justification": "As per the system’s instructions, paste my extended, prior final justification:\n\nI am satisfied with the authors’ responses provided during the rebuttal period. I have also reviewed the questions, comments, and exchanges between the other reviewers and the authors. Several noteworthy points were raised, and in my view, the majority of the authors’ responses were appropriate and well-reasoned. I therefore maintain that my original assessment remains valid.", "rebuttal": "We thank the reviewer for the time and the detailed review of our work.\n\n**Regarding the choice of K**\nThis is a very good and important point. Due to the page limit we discuss the choice of K in the appendix (A.4). We do not consider the Dirichlet Process Mixture but rather a more classical and simpler approach from the Empirical Bayes literature. \n\nWe write in the appendix (A.4)  \"A noteworthy aspect of empirical Bayes methods for matrix factorization, highlighted in Bishop 1999, is their inherent ability to automatically determine $K$ . This feature comes from the fact that the maximum-likelihood estimates for $g_k^{(\\ell)}$ and $g_k^{(f)}$ can converge to a point-mass at zero (assuming that the families $G_{\\ell, k}$ and $G_{f, k}$ include a point-mass at zero). Thus, if $K$ is initially set to a large value, certain row or column factors will converge to zero, which effectively \"zeroes out\" the corresponding\ncomponent $k$.\n\nSo, as long as the prior model family contains the point mass at 0, then cEBMF automatically estimates K. This means that if K is set sufficiently large, then some loading/factor combinations will be optimized to be exactly 0.\n\n**Regarding the choice of prior**\n\nCertainly, other priors could be used and may be more appropriate for different data sets or for different types of side information. We have not considered the use of hierarchical Dirichlet models nor multi-kernel Gaussian processes, and it is hard to comment on the benefits of such approaches without being more specific about the type of data or side information (categorical, spatial, temporal, etc) being considered. It is straightforward to use cEBMF with hierarchical Dirichlet models or multi-kernel Gaussian processes prior, given that you can solve the cEBNM mapping for these priors. So, more generally, cEBMF can be used with any prior model in which the cEBNM mapping can be computed or at least numerically approximated (see equation 13).\n\n**Regarding the initialization**\n\nThis is good point raised by the reviewer. For unconstrained factorization, we simply start by performing a truncated SVD of the data matrix and using the U and V  matrices as initialization for L and F. For non-negative matrix factorization and semi-non-negative factorization, we use a fast NMF routine (in R or in Python) to provide the initialization for L and F. We found that when using cEBMF for non-negative matrix factorization, using fast NMF routines for the initial starting point for L and F helps compared to using SVD initialization.\n\nFor the covariate-moderated priors model, we have not experimented with specific initialization, but for future work, we aim at using a pretrained model for specific applications (e.g., a model for cell recognition), and we hope to obtain some speedup as well as performance improvement from those.\n\nWe will add those details in the appendix in the final version.\n\nThe reviewer wrote: **\"A common challenge in matrix factorization models is dealing with scenarios where the number of latent factors is relatively large compared to the available data size\".**\n\nWe think that, in general, the main goal is to find a lower space representation of a relatively large data matrix. So, usually k is small compared with the data size. For example, in single-cell applications, the number of genes is >10000, the number of cells is >10000, but k is of the order 100.\n\nWe hope that we have answered your questions and that you may consider increasing our rating."}, {"review_id": "eiRbB1Fq2q", "reviewer": "Reviewer_KUPC", "summary": "This paper proposes Covariate-moderated Empirical Bayes Matrix Factorization (cEBMF), a general framework that extends EBMF by allowing covariate-dependent priors on latent factors. The key idea is to learn a mapping from covariates to priors via flexible models (e.g., neural networks), enabling the integration of high-dimensional side information. Empirical results demonstrate strong performance across synthetic and real-world datasets, with clear gains over competing methods and favorable scalability.", "questions": "1. Add experimental validation to demonstrate that the proposed method can leverage more complex side information, such as incorporating semantic and visual information of items in the task of collaborative filtering.\n2. Add discussion on why an end-to-end neural network is not directly used to model all types of information and accomplish the target task.\n3. Add some new baselines about deep learning.\n4. Add the discussion of the tasks that Matrix Factorization (MF) are suitable for.", "limitations": "Even with the introduction of side information, matrix factorization (MF) remains a linear model, which seems to limit its applicability.", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths\n1. This paper introduces cEBMF, a general framework integrating any side information (text, images, graphs) via covariate-dependent priors. This surpasses existing methods limited to specific data types (e.g., spatial coordinates).\n2. This paper allows plug-and-play priors, enabling adaptation to diverse problems.\n3. The proposed method is supported by solid theoretical foundations.\n\nWeaknesses\n1. This method can be regarded as a two-stage approach: it first uses a neural network to process side information, treating it as a prior, and then performs matrix factorization. I think the paper lacks sufficient discussion on why an end-to-end neural network is not directly used to model all types of information and accomplish the target task.\n2. The deep learning baselines (e.g., VAE, NCF) are relatively simple and may not fully represent the state of the art. The paper would benefit from a more nuanced positioning.\n3. Lacks analysis of the experiment results — why the proposed linear model outperforms nonlinear models in recommendation tasks.", "was_revised": true, "final_justification": "My concerns still remain in spite of the author reponse. So I would keep my original score of borderline reject.", "rebuttal": "We thank the reviewer for the time and the detailed review of our work.\n\nWe agree on the fact that the semantic and visual information of items was not used in the collaborative filtering task. Collaborative filtering is not our main field. Could the reviewer suggest some good benchmarks that have this information? \n\nThe reviewer wrote, **\"This method can be regarded as a two-stage approach: it first uses a neural network to process side information, treating it as a prior, and then performs matrix factorization. I think the paper lacks sufficient discussion on why an end-to-end neural network is not directly used to model all types of information and accomplish the target task. \"**\n\nThe reviewer might be unfamiliar with Empirical Bayes procedures, but this two-step procedure is the very essence of Empirical Bayes (EB) (we recommend Efron's book on Large-scale inference for an introduction). The nice aspect of EB is that it doesn't require the user to specify some hyperparameters prior to the analysis. Further, the reviewer may have missed that the prior model (here neural nets) is retrained whenever we are fitting a new factor, and the overall procedure is not quite as described by the reviewer.\n\nThe reviewer wrote **\"Lacks analysis of the experiment results — why the proposed linear model outperforms nonlinear models in recommendation tasks.\"**\n\nWe agree that it may at first glance look surprising, but there are several instances of real data examples in which linear models are ahead of non-linear models. For example, in genetics, polygenic risk score methods (which are linear models) are still much better than the nonlinear/ML-based methods despite having millions of data points.\n\nWe hope that we have answered your questions and that you may consider increasing our rating."}]

Author rebuttal:

Meta-review: This paper introduces cEBMF, an empirical Bayes matrix factorization framework that flexibly incorporates side information into priors. Reviewers agreed that the framework is novel, modular, and well-presented, with empirical results across simulations, MovieLens, and spatial transcriptomics, as well as robustness to prior or covariate misspecification.

The main point of contention comes from Reviewer KUPC, who questions the positioning of cEBMF relative to modern end-to-end deep models and notes that the current evaluation (e.g., MovieLens 100K with limited covariates) may not fully substantiate broad claims. Other reviewers raised technical questions about robustness and initialization, which were addressed in rebuttal, and overall leaned to acceptance.

On balance, I find the contribution substantive and broadly useful. It offers a flexible extension of empirical Bayes factorization that accommodates diverse side information. While additional comparisons to state-of-the-art deep baselines would further strengthen the work, the majority of reviewers support acceptance, and I concur.
Decision: Accept (poster)