Sparse encoding for more-interpretable feature-selecting representations in probabilistic matrix factorization
Authors: Joshua C Chang
Patrick Fletcher
Jungmin Han
Ted L Chang
Shashaank Vattikuti
Bart Desmet
Ayah Zirikly
Carson C Chow
Venue: ICLR 2021 Poster (2021)
Keywords: poisson matrix factorization
generalized additive model
probabilistic matrix factorization
bayesian
sparse coding
interpretability
factor analysis

Abstract: Dimensionality reduction methods for count data are critical to a wide range of applications in medical informatics and other fields where model interpretability is paramount. For such data, hierarchical Poisson matrix factorization (HPF) and other sparse probabilistic non-negative matrix factorization (NMF) methods are considered to be interpretable generative models. They consist of sparse transformations for decoding their learned representations into predictions. However, sparsity in representation decoding does not necessarily imply sparsity in the encoding of representations from the original data features.  HPF is often incorrectly interpreted in the literature as if it possesses encoder sparsity. The distinction between decoder sparsity and encoder sparsity is subtle but important. Due to the lack of encoder sparsity, HPF does not possess the column-clustering property of classical NMF -- the factor loading matrix does not sufficiently define how each factor is formed from the original features. We address this deficiency by self-consistently enforcing encoder sparsity, using a generalized additive model  (GAM), thereby allowing one to relate each representation coordinate to a subset of the original data features. In doing so, the method also gains the ability to perform feature selection. We demonstrate our method on simulated data and give an example of how encoder sparsity is of practical use in a concrete application of representing inpatient comorbidities in Medicare patients.

Published as a conference paper at ICLR 2021
SPARSE ENCODING FOR MORE-INTERPRETABLE FEATURE-SELECTING REPRESENTATIONS IN PROBA- BILISTIC MATRIX FACTORIZATION
Joshua C. Chang, Patrick A. Fletcher, Jungmin Han
Ted L. Chang meds rrata ted @ @mederrata.com
National Institutes of Health & mede mata (josh patrick jungmin| @ @mederrata.com
Shashaank Vattikuti
Bart Desmet, Ayah Zirikly, Carson Chow
Walter Reed Army Institute of Research & medé
National Institutes of Health
bart. ayah.zirikly carson.che chhoo/ @ @nih.gov
shashaank @ mederrata.com
ABSTRACT
Dimensionality reduction methods for count data are critical to a wide range of applications in medical informatics and other fields where model interpretability is paramount. For such data, hierarchical Poisson matrix factorization (HPF) and other sparse probabilistic non-negative matrix factorization (NMF) methods are considered be interpretable generative models. They consist of sparse transforma- tions for decoding their learned representations into predictions. However, sparsity representation decoding does not necessarily imply sparsity in the encoding of representations from the original data features. HPF is often incorrectly inter- preted in the literature as if possesses encoder sparsity. The distinction between decoder sparsity and encoder sparsity is subtle but important. Due to the lack of encoder sparsity, HPF does not possess the column-cluttrring property of classical NMF - the factor loading matrix does not sufficiently define how each factor is formed from the original features. We address this deficiency by self-consistently enforcing encoder sparsity, using generalized additive model (GAM), thereby allowing one to relate each representation coordinate to a subset of the original data features. In doing so, the method also gains the ability to perform feature selection. We demonstrate our method on simulated data and give an example of how encoder sparsity of practical use concrete application of representing inpatient comorbidities in Medicare patients.
INTRODUCTION
For many inverse problems featuring high-dimensional count matrices. such as those found in healthcare, model interpretability is paramount. Building interpretable nigh-performing solutions is technically challenging and requires flexible frameworks. A general approach to these problems to structure solutions into pipelines; if each step interpretable, one can achieve interpretability of the overall larger model. A common first step in modeling high- sigh--dmensional data sets is to use dimensionality reduction to find tractable data representations (also called factors or embeddings), that are then fed into downstream analyses. Our goal is to develop a dimension reduction scheme for count matrices such that the reduced representation has an innate interpretation in terms of the original data features.
Interpretability versus explainability. We seek latent data representations that are not only post-hoc explainable (Laugel et al.) 2019 Caruana et al. 2020), but also intrinsically interpretable (Rudin 2019). Our definition of intrinsic interpretability requires clarity in the relationship between predictors and prediction, and meaningfulness of interactions and latent variables.
Post-hoc explanations are based on subjective examination of a solution through the lens of subject- matter expertise. For black-box models that lack intrinsic interpretability, these explanations

---

Published as a conference paper at ICLR 2021
produced using inexact simpler approximating models (typically local linear regressions). These explanations can be misleading (Laugel et al. 2019).
Disentangled autoencoders. Disentangled variational autoencoders (Higgins et al., 2016 Tomczak & Welling 2017 Deng et al. 2017) are deep learning models that are inherently mindful of post-hoc model explainability. Like other autoencoders, these models encoder-decode structured (see Definitions 2nd where the encoder generates dimensionally reduced representations.
Definition 1. The encoder transformation maps input data features to latent representations
Definition 2. The decoder transformation maps latent representations to predictions
Disentangled autoencoders use a combination of penalties (Higgins et al., 2016 Hoffman et al. 2017) and structural constraints (Ainsworth al.) 2018) encourage statistical independence in representations, facilitating explanation. These methods arose in computer vision and have demonstrated empirical utility in producing nonlinear factor models where the factors are conceptually sensible. Yet, due to the black-box nature of deep learning, explanations for how the factors are generated from the data, using local saliency maps for instance, are unreliable or imprecise (Laugel et al.] 2019 [Slack et al., 2020, Arun et al. (020). In imaging applications, where the features are raw pixels, this type of interpretability is unnecessary. However, when modeling structured data problems, one often wishes to learn the effects of the individual data features.
Probabilistic matrix factorization. Probabilistic matrix factorization methods are related to au- toencoders (Mnih & Salakhutdinov] 2008). These methods are often presented in the context of recommender systems. these cases, rows of the input matrix are attributed to users, and columns (features) are attributed to items. Probabilistic matrix factorization methods are bi-linear in item- and user-specific effects, de-convolving them in manner similar to item response theory (Chang et al.) 2019). In applications with non-negative data, non-negative sparse matrix factorization methods further improve on interpretability by computing predictions using only additive terms (Lee & Seung 1999).
For count matrices, Gopalan et al. (2014) introduced hierarchical Poisson matrix factorization (HPF). Suppose Y (yos) U matrix of non-negative integers, where each row corresponds to user and each column corresponds to an item (feature). Adopting their notation, Gopalan et al. (2014) formulated their model as
OuklKn, a ~ Gamma (a, E.. BEs Bealmic c Gamma (c,r)).
Musl O B Poisson
(I)
where θ = (0wx) ) is a x K matrix, and B (3xx is the representation decoder matrix. Additional priors 7); Gamma (, /d' and Eu Gamma(a', a' /b') model item and user-specific variability in the dataset, and a', V. d' E R are hyper-pant pareerrarrmtters... The row vector o,, = (fw1,..., iia....,,uuu Oak) constitutes K- K-dimensional representation of the user, and the matrix B (Bks) decodes the representation into predictions on the user's counts.
In HPF, the gamma priors the decoder matrix B (Bks) enforce non-negativity. Because the gamma distribution can have density at zero, these priors also allow for sparsity where only a few of the entries are far from zero. Sparsity, non-negativity, and the simple bi-linear structure of the likelihood in HPF combine to yield simple interpretation of the model: in HPF, predictive density for each matrix element is formed using linear combination of subset of representation elements, where the elements of B determine the relative additive contributions of each of the elements (Fig. [T). However, the composition of each latent factor in terms the original items is not explicitly determined but arises from Bayesian inference (Fig TP)..
Limitations of HPF. Classical non-negative matrix factorization (NMF) is often touted for having a column-clusterie property (Ding et al. (2005), where data features are grouped into coherent factors. The standard HPF of Eq. lacks this property. In HPF, while each prediction a linear combination of a subset of factors, each factor is not necessarily a linear combination of subset of features (depicted in Fig. [b).
The transformation matrix B defines a decoder (Def. 2 like classic autoencoder. A corresponding encoding transformation (Def. [ does not explicitly appear in the formulation of Eq. Determining the composition of factors is not simply matter of reading the decoding matrix row-wise. Mathemat-
2

---

Published as a conference paper at ICLR 2021
(a) standard HPF
(c) Interpreting factors
Aus
Ru2
Standard HPF
Auz
Ouz
0g a, (=(0, B....
(Y)
X
Incurrett:
Par
Yn?
encoding
A00
representation
decoding
prediction
(b) sparsely enneeelyeeeeded HPF
Spassly-ensooee encoded HPF
inw.
S (P
Avi
T
Out
3_.,
of
Also currect:
X
Mas
Vuto
Au3o
cuccoding
rxpresentation
decoding
prodiction
Figure 1: Interpreting hierarchical sparse probabilistic matrix factorization (HPF). (a) Stan- dard HPF: Rates Aui for Doisson-distribute predictions sparse linear combinations of the learned representation as defined by the decoding matrix; this matrix does not define how representations derived from the input data. (b) Sparsely-encodel HPF (proposed method); the mapping
from input data to representation is given explicitly by a sparse encoding matrix. (c) Interpreting representations: It is tempting but misleading read the decoding matrices row-wise in determining the feature subsets that contribute to forming representation coordinate. Representations 0,, are computed by inferring the statistics of an associated joint posterior distribution T(... Y) the sets of non-sparse entries in rows of the decoding matrices do not necessarily correspond to feature sets that determine the representation components. However, for sparsely-encode HPF, the representations are explicit functions of subsets of features.
ically, sparsity in decoding does not imply sparsity in encoding analogous how pseudo-inverses of sparse matrices are not necessarily sparse.
HPF is also unable to perform feature selection. By Eq. predictions formed by weighting representations using values from columns of the decoding matrix a feature's corresponding terms in the decoding matrix will be near zero if and only if that feature's mean is zero. Exclusion of feature column from the decoding matrix yields no information on whether that feature plays part in generating representations. Also, this deficiency can cause HPF to erroneously imply structure when none is present, as demonstrated in Fig. 3 on factorization of pure Poisson noise.
Our contributions. We propose a method to consistssitt constrain HPF so that its corresponding encoding transformation is explicit. In doing so, we improve interpretability of HPF, and give it the ability to perform feature selection. Constraining HPF in this manner also makes it more suitable to training with large datasets because the representation matrix does not need to be stored in memory. Using medical claims case study, we demonstrate how our method facilitates reparameterization of decision rules in representation space into corresponding rules on the original data features.
2
METHODS
In this section we describe our extension to HPF that resolves its limitations. Our method yields representations that have explicit sparse dependence on relevant features the input data.
IMPROVING HPF BY CONSTRAINING IT 2.1
Our augmented HPF model takes the form
Oux PullA,yuse = En 2e
Vui B.B.p Poisson (A. (
Eoneees +pps
(2)
Bx Normal (0, 1/ax),
where prior distributions for the model parameters defined later in this section.
3

---

Published as a conference paper at ICLR 2021
The key point is that the encoder function (that computes Byk)) is an explicit function of the input data, formulated using a generalized additive model (GAM) (Rigby & Stasinopoulos 2005; Hastie & Tibshirani 1987 Klein et al. 2015). The encoding matrix A (ast) controls how features map into the representation.
To allow the model to perform automatic feature selection, also incorporate non-negative item-specific gain term Pi as a background Poisson rate for item é that is intended to be independent of the factor model. We also slightly generalize the likelihood of Eq. by by giving each feature an associated link function fi, which models nonlinearities without sacrificing interpretability.
The distributions of the parameters of the encoder are learned self-consistently with other model parameters. In the process, one is training not only the generative model, but also the subsequent Bayesian inference of mapping data representation by learning the statisties of the posterior distribution,
O Yu, Y ~ /f T B, ply YdBdpp
(3)
where the generative process has been marginalized. In short, the model of Eq. 2 usee the marginal posterior distribution of the encoding matrix A to reparameterize this Bayesian inference. Doing so amortizes this inference. making it trivial to apply the model to new data in order to compute new representations. It also allows us impose desirable constraints on the representations themselves.
In the original HPF, the parameters &. are used to account for variability in user activity (row sums). Similarly, the na parameters account for variability in item popularity (in column sums). To simplify the method, we pre-set these parameters (based on some training data) to fu = 1 and 7): = 0 Σ. Vui. encoder parameters Oss so we can generalize weakly-informative and other scale-dependent priors within the model to disparate datasets, as is common in preprocessing for Bayesian statistical inference problems (Gelman al.) (017). One may also model over-dispersed data using Eu = U Σ. Ba//., ... 2n Bucs to account for document-size variability.
where 7/s absorbed into the function f.. Doing so de-scales
We use sparsity to achieve feature selection. We encourage the elements ag and y to be mutually exclusive by using the decomposition
[s] $7 Horseshoe + (1,
Orok
Mik
+ S
[uis u2a uik Horseshoe + (1, 1/V/T
(4)
S
PilT, Wi, 8,
w; Normal* (0, 10).
where the non-negative version of the Horseshoe* + prior Carvalho et al. 2009 2010 Polson Scott 2011) is the hierarchical Bayesian model
ry A3, r ~ Normal (0, AçT)
x Horseshoe" + (Ao: To)
Ag Cauchy (0, Ao)
(5)
T ~ Cauchy (0, To).
This concentrates marginal distributions of vector components near zero. Additionally, it minimally shrinks large components, resulting lower bias compared to lasso and other alternatives (Bhadra) et al. 2015a b) 2019 Piironen & Vehtari 2017b). The horseshoe has previously been applied in other factorization methods, including autoencoders (Ghosh & Doshi-Veler, 2017b) and item response theory (Chang et al. 2019), but not to probabilistic matrix factorization.
Applied to the parameters s sparsity discourages variables that load into the factor model from leaking into the corresponding background rate term pp, Conversely, variables that load into i are discouraged from appearing in agk-
Finally, is often done in variational autoencoders, regularize the representation by placing unit half normal priors on its components
0u Y, - Σ gi (Juc)oua Normat (0, 1).
The choices of the encoding function f. and decoding functions 9; are application-specine and may be learned (Rigby & Stasinopoulos 2005). So not to distract from our focus improving
4

---

Published as a conference paper at ICLR 2021
the interpretabilityy of pre-existing matrix factorization approaches, fix these functions here. In standard Poisson matrix factorization approaches, f((z) = 9:(x) = I, Vi. Equivalently, we choose to rescale the inputs that
fa(z) = 7/5 and 9,(x) = f,- (z) = -/m
(7)
Another choice for these functions can be motivated by Poisson regression with logarithmic log (*/n. + 1). For maximum interpretability, restricting fi and 9: to monotonically increasing functions where 9:(0) fa( 0
link function, where fi(x) e - 1, and (z) = (x)
results in order-preserving representations that are zero when the corresponding feature counts are zero.
2.2 INTERPRETING REPRESENTATION AND DERIVED QUANTITIES
In constraining the encoder mapping using the generalized additive model of Eq. we regain the column-clustering property of classical non-negative matrix factorization methods: each representa- tion component is explicitly determined from well-defined subset (cluster) of the data features. In Fig. [T), we demonstrate how the encoding matrix can be read determine the composition of the factors. Consequently, decision rules over the representation can be easily expressed as decision rules over the original features,
Ouk E (a, - Σ gy(Un))agk E (4/Eu,//e.)
JEZt
where S2 is the subset of features that determines factor k. As we will demonstrate in our main case study, Eq. [ iis useful for inverting clustering rules defined over the representations.
2.3 INFERENCE
The model of Eq. 2 is a generalized linear factor model that we have mathematically related to probabilistic autoencoder. When augmenting HPF with explicit encoder inference, as we have done, one obtains a probabilistic autoencoder. This suggests that previous work can serve as a guide for training, especially work done on using the horseshoe prior Bayesian neural networks (Ghosh &
Velez 2017a Ghosh et al. 2018 Louizos et al. 2017).
In particular, Ghosh et al. (2018) investigated structured variational approximations of inference of Bayesian neural networks that use the horseshoe prior and found them to have similar predictive power as mean-field variational approximations. The disadvantage is the extra computational cost of inferring covariance matrices. For these reasons, we focus on mean-fiell black-box box variational inference, using Ghosh et al. (2018) as a guide, noting consistency of their scheme with other works that have investigated variational inference problems using the horseshoe prior (Wand et al. 2011 Louizos et al. 2017).
structured approximations
As in Ghosh & Doshi-Veler! (2017a): Ghosh al. (0018); Chang al. (2019), for numerical stability, we reparameterize the Cauchy distributions in terms of the auxiliary inverse Gamma representation (Makalic & Schmidt 2016).
x² Inverse-Gamma (1/2, 1/x) Severse-Gamms (1/, 1/o9).
x ~ Cauchy+ (0,
We perform approximate Bayesian inference using fully-factorized mean-field Automatic Differenti- ation Variational Inference (ADVI) (Kucukelbir et al. 2017). For all matrix elements, we utilized Gaussians, and coupled these to inverse- Gamma distributions for the scale parameters, as investigated in Wand al. (20IT). For our use cases, we implemented a minibatch training regimen common to machine learning, with stepping given by the Adam optimizer (Kingmal Ba 2017) combined with the Lookahead algorithm (Zhang al. 2019) for stabilization.
Bayesian sparsity methods concentrate marginal distributions for parameters near zero. To further refine the Bayesian parameter densities so that some parameters are identically zero in distribution, one may use projection-based sparsification (Pironen & Vehtari 2017b). Finally, one can assess predictive power without model refitting using approximate leave-one-out cross validation (LOO) using the Widely Applicable Information Criterion (WAIC) (Watanabe) 2010;; Gelman al. 2014
Pironen & Vehtari 2017a Vehtari et al., 2017 Chang 2019) or Pareto smoothed importance
sampling 1.00 (PSTS-L.OO) VVettari et al. 2017)
5

---

Published as a conference paper at ICLR 2021
3 EXPERIMENTS
We implemented the method in censorfoo-prrbbablii (Dillon et al. 2017), modifying the inference routines for implementing our variational approximation. Our implementation be found at github: : mederr spmf along with notebooks reproducing our simulation results. We present here simulation results and application to medical claims data.
3.1 SIMULATIONS
To demonstrate the properties of our method, we factorized synthetic datasets of: a) completely random noise with underlying structure, b) mixture of random noise and linear structure, and a mixture of random noise and nonlinear structure.
( bote
Leneatly
() Nonirearly
buttes
beckarinnd
1/0
0.1 56 03
12
14
06
Figure 2: Factorization of simulated datasets. The (mean) effective encoding matrix A = (assk) ) for each factor process, placed common color scale, and the posterior distribution of the background process rate i by item for a) Poisson(1) noise, where there is no relationship between the features, b) linear factor model where every third variable is generated from dense factor model and the other variables are Poisson(1) noise, c) nonlinear factor model where every third variable is generated from dense nonlinear factor model and the other variables are Poisson(1) noise. See Fig. 3 for standard HPF on these datasets for comparison.
For a), we sampled a 50, 000 x 30 Poisson(1) random matrix. Figure 2h) shows the inferred mean encoder matrix A along with the posterior distributions for each of the background components pa, We see that all features are excluded from the encoding matrix, showing up instead as background noise.
Next, we created test system where there is underlying linear structure mixed with noise. For this system, we again used I = 30 features and put every third feature into dense system by generating random 10 x 10 decoding matrix B, sampling representations from a non non-negative truncated normal distribution, and sampling counts according to the generative process of Eq. [ For the remaining features, we used Poisson(1) noise. After simulating 50, 000 records by this process, we performed factorization again into K = 3 dimensions. The results of this factorization are shown in Figure ( where it is clear that every third feature falls into the overall factor model and the remaining features show up as background noise.
As an example of factorization under model mismatch, we generated random data with derlying nonlinear structure. Here again, we used every third feature, simulating and θ before. However, simulated counts for these features using the model Yur Poisson (Ex OutBxs/2) (xp(-EX E₁ (u40e/22 + (Ex (=r8me/2))) Factorization of this dataset is shown in Figure Again, it is clear that every third feature falls into the overall factor model and the remaining features are load into the background process with rates near 1, indicating that even when the model is mis-specified, can successfully separate structure from noise. In Supplemental Fig. SI present the same factorizations while using the logarithmic link function described in demonstrating robustness to mis-specification of the link functions as measured using the
Sec.
WAIC.

---

Published as a conference paper at ICLR 2021
a) Poiesoni L) noise
b) Lineslly strnectreed
30
25
06 09
Nonlinoarly stractared
Figure 3: Decoder matrices B (Bes) for standard HPF factorization of the synthetic datasets of Fig. using the python package hpfrec. Shown are posterior means.
Comparison to standard HPE. In standard HPF (Gopalan et al. 2014) only decoder matrices are inferred, and encoders synthetic datasets using standard hierarchical Poisson matrix factorization (Gopalan et al.] 2014) found in the hpf rrec package. In all three examples, the standard HPF fails remove independent noise items from the factor model. In contrast, our method excludes all irrelevant features from the factor model (Fig. 2). In this case, the encoder not explicitly solved. decoder matrix B row-wise, to say for instance that the first factor in Fig, 33 is determined from items (2, 3, 7, ... 29). However, results from standard HPF are often erroneously interpreted in this manner, suggesting that there is structure to the dataset even when none is present. Additionally, for this reason, the generative process fails to adequately fit the data in that it correlates sources of independent noise.
explicitly reconstructed. Fig. demonstrates factorization of the same
incorrect read the
3.2 COMORBIDITIES FROM BILLING CODES
As real-world case study, we used a 5% sample of the Medicare Limited Data Set (LDS) over the years 2009 201 1 to discover representation of inpatient comorbidity during a hospital visit. The LDS consists of de-identified medical claims data for Medicare and Medicaid recipients across the United States. Pursuant to a data use agreement, the Center for Medicare and Medicaid Services (CMS) provides a 5% sample of this dataset for research purposes. single hospital visit consists of multiple claims across providers and types of services. No standard method for grouping claims into hospital visits within claims data exists. heuristic algorithm (often called a grouper algorithm) used reconstruct medical and billing events during visit type of service from claims. Our grouper algorithm collapsed the claims into = 1, 949, 788 presumptive inpatient visits. Diagnostic codes were then made bourser-grained by backing off from the original 13,000 ICD-9-CM to 136 clinically-releannt categories using the top two levels of the CCS multilevel classification. Within each visit, counted the number of codes that fell into each of the CCS categories. Fig. 4 presents the encoding matrix A = (ask) for factorization of comorbidities into four dimensions, the transpose of the decoding matrix B (Bkz). and the vector of background process rates 4P (vx) for the same model.
The values in the encoding matrix provide coefficients that are used in Eq. 6lto produce weighted sum of billing code counts, which is then used to formulate a representation. One may read the encoding order to determine the feature composition of each of the representation factors.
matrix column-wiss Being able to do so facilitates interpretation of the factor model, by allowing one to understand single factor at a time by focusing on subsets of the original features. For example, conceptually, easy to what factor represents. It is computed by tallying up various billing codes that pertain to lung ailments - lung cancer (CCS 2.3), several broad respiratory disorders (CCS 8.x), and heart disease (CCS 7.2). The relative weights of these codes are depicted by the color of the shading.
The interpretation of matrix factorization is not provided by the decoding matrix, which the sole output of standard HPF. The decoding matrix provides only incomplete picture of the structure of the data. Recall from Fig. T that the decoding matrix not to be read row-wise (or column- wise in the transpose depiction of Fig. Looking solely the decoding, and reading it in this incorrect way,
7

---

Published as a conference paper at ICLR 2021
Medicare patient comorbidity
Encoding (A)
Decoding
Infectious parashic diseases
Nesplasms
Endoctine, nutiitional metabellc diseases &
Bload blood-forming or gans
Mental
Nervous system disorders
Cire datory cisorders
Respiratory disorders
Digestive
Genitsuriarry pregnancy;
childbith: & the
subostaneous bissue disorder
Musculaskkeetal
& connective bssile disorder
Congerital anomalies
Cortain conditons
orio hating perinatal peroo
Injury and paiscring
Symptoms, signs, & illidefnedd conditions inflaencing Unchasstied
Figure 4: Medicare comorbidity factorization for inpatient visits based on medical claims from a 5% sample of the Medicare Limited Dataset (LDS), in four factor dimensions. Prior to factorization, we mapped each raw ICD diagnostic code into the second tier of the Clinical Classification Software (CCS), counting the number of codes present within each broad category. Shown are posterior means. Left: encoding A = (axc), middle: decoding B = (Bes)" right: background y = (yps)
one might erroneously conclude that lower respiratory disorders (CCS 8.8) are the main determinant of factor 2. However, this conclusion is incorrect diagnoses of heart disease (CCS 7.2) are the main determinant.
The decoding matrix also provides misleading insights on feature selection. For example, lung cancer (CCS 2.3) appears only very faintly in the decoding. For this reason, one might come to the erroneous conclusion that it is not important will only faintly register in the decoder matrix. The rate of lung cancer diagnoses is low, yet lung cancer diagnoses are predictive of and coincide with other respiratory issues. Hence, lung cancer (CCS 2.3) appears strongly in the encoder matrix.
feature. However, recall that by Eq. [ relatively rare features
After computing representations for the entire dataset, one may cluster patients into diagnostic groups. One way of doing so is through stratification, for instance into low, medium, and high groups for each of the four factors in Fig. 4 Doing so based on quantiles yields thresholds between the groups, defined over representations. Using Eq. 8 one can easily convert these thresholds into decision rules on the counts. For example, patient presenting with one or more lung cancer-related diagnoses would generally be placed in the medium or high strata for factor 2, depending on the number of other respiratory billing codes in that visit. Recall that the input data is sparse so in general the low strata for each representation would encompass people who had no or very few of the associated diagnoses. As a first step a modeling pipeline. one could use the strata to segment a larger overall model, so that the model has both local and global behavior, while making it easy to interpret how individual or collective billing codes contribute to an overall prediction.
DISCUSSION AND CONCLUSION
We introduced a constrained HPF where an encoder transformation is learned self-consistently with the matrix factorization. By imposing sparsity on the encoder, rather than on the decoder, improve
8

---

Published as a conference paper at ICLR 2021
interpretability of HPF. We demonstrated the approach simulated data, showing that the method can successfully separate structure from noise, even when the model is mis-specified. We also presented comorbidity factorization case study.
Although we focused on Poisson factorization, our central argument holds for other sparse matrix factorization methods. Sparse decoding matrices (loading matrices) inferred using these methods are generally not orthogonal. Unlike in classical or orthogonally PCA, the transpose of these decoding matrices does not correspond to their pseudo-inverse. Hence, decoding matrices should never be interpreted row-wise (Fig. Tc).
LIMITATIONS AND EXTENSIONS
Our method relies on the horseshoe prior for sparsification. The horseshoe prior relies on scaling hyperparameters, which control the effective sparsity of the method. In order to make these priors scale asymptotically with data size (Piironen & Vehtari 20176), we chose to scale this prior using 3/VUI. Empirically, this choice. along with the scaling of the priors on 8k4, led to desirable behavior in simulations like those in Fig. 2 under several combinations of U and I. In effect, we have taken the liberty of formulating our method based on these considerations so that it is usable without needing to manually choose hyperparameters. One may wish to rescale the horseshoe prior in order control sparsity. Further guidance to the regularization scale will require analysis outside of the scope of this manuscript.
We note that one could also formulate our method using the sparsifying priors found in Gopalan et al. (2014). The chief advantage of the original HPF formulation is in how yields explicit variational inference updates. However, our method yields well to ADVI, achieving convergence with learning rate of 0.05 in approximately 100 epochs in all included examples.
A limitation of our method, shared by standard HPF, is that a generalized linear model does not have the expressivity of nonlinear paradigms such as deep learning. For some applications, with sufficient data, nonlinear models may be more performant. We note that one could place non- non-linearity in either Sa or 9i. without compromising interpretability of the representation. These functions may be learned using Gaussian processes (Chang et al.) 2014) splines, even neural networks, making the method more like other probabilistic autoencoders. So long as the conditions of monotonicity and a fixed point at y 0 are maintained, the overall method remains interpretable. However, the simplicity of HPF offers statistical advantages that help it generalize better than deep learning except when there is enough data to learn any true nonlinearities in the true generating process. Additionally, while we do not explore this, a strength of Bayesian modeling is that it provides a principled approach to incorporating prior information. One could encourage or discourage features from co-factoring by setting suitable priors on the encoding and decoding matrices.
ACKNOWLEDGEMENTS
We thank the Innovation Center of the Center for Medicare and Medicaid services for providing access to the CMS Limited Dataset. We also thank the mede team (particularly Joe Maisog) for their support in various data pre-procd pre-processing tasks. JCC, AZ, and BD were supported in-part by the US Social Security Administration. PF, JH. CCC, and SV were supported by the Intramural Research Program of the NIH, NIDDK. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) (Towns et al., 2014), which is supported by National Science Foundation grant number ACI- ACI-1548562 through allocation TG-DMS 19002. We also thank Amazon Web Services for providing computational resources and Boost Labs LLC for helping with visualization.
REFERENCES
Samuel Ainsworth, Nicholas Foti, Adrian KC Lee, and Emily Fox. Interpretable VAEs for nonlinear group factor analysis. 007.006802 [cs, stat], February 2018. URL arXiv: 1802.06765.
abs/18022 06765
Nishanth Arun, Nathan Gaw, Praveer Singh, Ken Chang, Mehak Aggarwal, Bryan Chen, Katha- rina Hoebel, Sharut Gupta, Jay Patel, Mishka Gidwani. Julius Adebayo, Matthew D. Li, and Jayashree Kalpathy-Crames Assessing the (Un)Trustwrtthines of Saliency Maps for Lo-

---

Published as a conference paper at ICLR 2021
calizing Abnormalities in Medical Imaging. arXiv:2008.0276: [cs], August 2020. URL
http://arxiyy httpxiiv.ogg/bbb/000533 .270008..7666
arXiv: 2008.02766.
Anindya Bhadra, Jyotishka Datta, Nicholas G. Polson, and Brandon Willard. The Horseshoe+
Estimator of Uptra-Sparse Signals. arXiv:/502.000 50500000000000 [math, stat), February 201 URLhttp:// arXiv:
arxiv. org/ abs / 1502. 00560 1502.00560.
tex.ids:
Anindya Bhadra, Jyotishka Datta, Nicholas G. Polson, and Brandon T. Willard. Default Bayesian analysis with global-local shrinkage priors. arXiv: 15 10.03516 [stat], October 2015b. URL
http:///arruuu.rrr/////5130... 03516
arXiv: 1510.03516.
Anindya Bhadra, Jyotishka Datta, Nicholas G. Polson, and Brandon Willard. Lasso Meets Horse- shoe: A Survey. Statist. Sci., 34(3):405-44) 427, August 2019. ISSN 0883-4237. 2168-8745, doi: 10.1214/19-STS70.. URLhttps:: http///ootteecici..ccccc wttp///cei...aooee..o.. org/
Publisher: Institute of Mathematical Statistics.
Rich Caruana, Scott Lundberg, Marco Tulio Ribeiro, Harsha Nori, and Samuel Jenkins. Intelligible and Explainable Machine Learning: Best Practices and Practical Challenges. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD '20, 35 3-11-2 New York, NY, USA, August 2020, Association for Computing Machinery. ISBN 978-1-4503-7998-4 doi: 10.1145533944486.466655 URL (tRttttt::://oot http:tts///....///0000..
3394486. 3406707
Carlos M. Carvalho, Nicholas G. Polson, and James G. Scott. Handling Sparsity via the Horseshoe. In Artificial Intelligence and Statistics, pp. 73-80, April 2009, URL/sttp://procee Section: Machine Learning.
mlr. wpwss////carvaaaaaeeaa
html ISSN: 1938-7228
Carlos M. Carvalho, Nicholas G. Polson, and James G. Scottt. The horseshoe estimator for sparse
signals. Biometrika, 97(2)4655-45) 2010. ISSN 0006-3444. /tLhhtts::/ /etttps:/www... www.wsttor..
orrg stable/2573409
Publisher: [Oxford University Press, Biometrika Trust].
Joshua C. Chang. Predictive Bayesian selection of multistep Markov chains, applied to the
detection of the hot hand and other statistical dependencies in free throws. Royal Soci- doi: 10.10988/sos..821174 182174
ery Open Science, 6(3):182174, March 2019.
URL https
/
Joshua C. Chang, Van M. Savage, and Tom Chou. A Path- Integral Approach to Bayesian Inference for Inverse Problems Using the Semiclassical Approximation. Stat Phys. 15733:582-602, November 2014. ISSN
doi: 10. 10.1007/s10955-014- 1007/31095-014-1053
URL httpt:::/tpp.ooroeeeeoo
Joshua C. Chang, Shashaank Vattikuti, and Carson C. Chow. Probabilistially-attteecoo horseshoe- disentangled multidomain item-response theory models. arXiv: 1912.0231 [cs, stat], December 2019. URL http:
www.oorg//aa.... /abs////a/////// 1912. C2351 arXiv: 1912.02351.
Zhiwei Deng, Rajitha Navarathna, Peter Carr, Stephan Mandt, Yisong Yue, Iain Matthews, and Greg Mori. Factorized Variational Autoencoders for Modeling Audience Reactions to Movies. In 20/7 IEEE Conference on Computer Vision and Patlern Recognition (CVPR), pp. 6014-6023, July 2017. ISSN: 1063-6919.
doi: 10.1 D..1099/VVP.2207..633
Joshua V. Dillon, Ian Langmore, Dustin Tran, Eugene Brevdo, Srinivas Vasudevan, Dave Moore, Brian Patton, Alex Alemi, Matt Hoffman, and Rif A. Saurous. TensorFlow Distributions.
arXiv:1711.10604 [cs, stat], November 2017. URL http:://axiiii wwa.xi.oor//ab//1111 10604 arXiv: 1711.10604.
Chris Ding, Xiaofeng He, and Horst D. Simon. On the equivalence of nonnegative matrix factorization and spectral clustering. In in SIAM International Conference on Data Mining, 2005.
Andrew Gelman, Jessica Hwang, and Aki Vehtari. Understanding predictive information criteria for Bayesian models. Stat Comput, 24(6):997-0016 November 2014. ISSN 1573-1375, doi: 10.10 s1 112220033-941622 /ttp:////////////ooo////////00/0
10

---

Published as a conference paper at ICLR 2021
Andrew Gelman, Daniel Simpson, and Michael Betancourt. The Prior Can Often Only Be Understood in the Context the Likelihood. Entropy, 19(10))555, October 2017. doi: 10.3390/e191003 URL / wwww mdpp. com. 1099- 4309994400/1110/555353 Number: 10 Publisher: Multi-
disciplinary Digital Publishing Institute.
Soumya Ghosh and Finale Doshi-Velez. Model Selection in Bayesian Neural Networks via Horseshoe Priors. DrXiv:/705.10388 [stat]. May 2017a. UKLhttoe//axx....r///////1100...
arXiv: 1705.10388.
Soumya Ghosh and Finale Doshi-Veled. Model Selection in Bayesian Neural Networks via Horseshoe Priors. May 2017b.
Soumya Ghosh, Jiayu Yao, and Finale Doshi- Structured Variational Learning of Bayesian Neural Networks with Horseshoe Priors. arXiv://800.55975 [cs, stat], June 2018. URL http
7/arxiv. wttp:/..rrrabbs/110050 05975 1806.05975.
tex.ids: ghashhstnntturrddaaaattiaaa...n earning201 arXiv:
Prem Gopalan, Jake M. Hofman, and David M. Blei. Scalable Recommendation with Poisson Factorization. 1rXivv131.11000 [cs, stat), May 2014. URL /tL/ttt:://axxxv.crr///oa 1704 arXiv: 1311.1704.
Trevor Hastie and Robert Tibshirani. Generalized Additive Models: Some Applications. Jour- nal of the American Statistical Association, 82(398):37-336,
June 1987. ISSN 0162-
1459. doi: 021.0//00001119010.112 129872249599777.. 10478440 Publisher: Taylor Francis _eprint:
1171:0011111410199071177878941 URL https: /wwwww www.wttaddoonine...mm
1ttp://11.1000001111111111 http://www.aaaaiii.cccmmooocooio 10474444
Irina Higgins, Loie Matthey, Arka Pal, Christopher Burgess, Xavier Glorot, Matthew Botvinick,
Shakir Mohamed, and Alexander Lerchner. Det: Learning Basic Visual Concepts with
Constrained Variational Framework. November 2016. URL http:::://ppeeeee.....ee
Eorum?1d=Sy2f2U9gg]
Matt Hoffmman, Carlos Riquelme, and Matthew Johnson. The Beta VAE's Implicit Prior. 2017. URL
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization. 41Xiv:1412.6980 arXiv: 1412.6980.
[cs], January 2017. URL http::tt::// //axxiv. org/ abtpaaass./1./eeeea//ee/aaaaaata 6980
Nadja Klein, Thomas Kneib, and Stefan Lang. Bayesian Generalized Additive Models for Location, Scale, and Shape for Zero-Inflated and Overdispersed Count Data. Journal of the American Statistical Association, 110(509):405 419, January 2015. ISSN 0162-1459 doi: 10.1080/
21621459.2014.912955 URL httts: //coi orq / 10. http:///........0000000001000000
2014. 912955
Publisher: Taylor & Francis _eprint:
Alp Kucukelbir, Dustin Tran, Rajesh Ranganath, Andrew Gelman, and David M. Blei. Automatic differentiation variational inference. J. Mach. Learn. Res., 18(1):430-74, January 2017. ISSN 1532-4435.
Thibault Laugel, Marie-Jeanne Lesot, Christophe Marsala, Xavier Renard, and Marcin De- tyniecki. The Dangers of Post-hoc Interpretability: Unjustified Counterfactual Explanations. (00ivv9907.0922 [cs, stat], July 2019. URL http://aaxxiv.orrr/bb//// /992000099/ 09294 arXiv: 1907.09294.
Daniel D. Lee and H. Sebastian Seung. Learning the parts of objects by non-negative matrix factorization. Nature, 401((655) 791, October 1999. ISSN 1476-4687. doi: 10. 1038/44565. URL :ttp:::://wwww.. Publishing Group.
com art cotmm/arttileee/455555 Number: 6755 Publisher: Nature
Christos Louizos, Karen Ulirich, and Max Welling. Bayesian Compression for Deep Leam- ing. 1rXiv:1705.088 (cs, stat], November 2017. URL (tLLkktt:///axxex http://axxxv..../ababaaaae/a 08665 arXiv: 1705.08665.

---

Published as a conference paper at ICLR 2021
Enes Makalic and Daniel F. Schmidt. A simple sampler for the horseshoe estimator. IEEE Signal January Process. Lett. 23(1):179- Januaryy 2016. ISSN 1070-9908, 1558-2361. doi: 10.1109/SSP.
2015.2503725. URL http: http:://axxi..orrtetet httpgg/bss/1000000 03884
arXiv: 1508.03884.
Andriy Mnih and Russ R Salakhutdinov. Probabilistic Matrix Factorization. In J. C. Platt, D. Koller, Y. Singer, and S. T. Roweis (eds.), Advances Neural Information Processing Systems 20,
pp. 1257-1264. Curran Associates, Inc., 2008. URL http: /ttp:r///pppp....pe nips httpp.ccc/ppppr/ 3208 probabilistic orrooba:lsscccaratri--fefrorroo.
Juho Piironen and Aki Vehtari. Comparison of Bayesian predictive methods for model selection. Stat Comput, 27(3) May 2017a. ISSN 1573-1375, doi: 10. 11577/111222211---11
URL
https://doo org.orrgggo 10.
Juho Piironen and Aki Vehtari. Sparsity information and regularization in the horseshoe and other shrinkage priors. Electron. J. Statiss., (10)00001100111110100555 2017b. ISSN 1935-7524. doi: 10.1214/
17-E/S1337SI, / http://pp.pbbba/oooobob.
arXiv: 1707.01694.
Nicholas G. Polson and James G. Scott. On the half-Cauchy prior for global scale parameter. 4937 arXiv:
arXiv: 11t4:491044477 [stat], April 2011. URL/tttp///xxxx..r/////3.
1104.4937.
R. A. Rigby and D. M. Stasinopoulos. Generalized additive models for location, scale and shape. Journal of the Royal Statistical Sociery: Series (Applied Statisties), 54
(3):507-554, 2005. ISSN 1467-9876.
doi: 10.1111/.1467-97777 001.6620000000000..0 URL 1111/ 1467 9876.
2005. 2005.001001000001111112211111111 .eprint:
TOP SERIES C RESEARCH 10..000000000.
Cynthia Rudin. Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nature Machine Intelligence, ((5)206-2155 May 2019. ISSN 2522- 5839, doi: 10. 12.008//42256-019-00446 s42256 019-0048--001
URL https: / www. naturatuueee.. wtt..oommtttttee//co Number: 5 Publisher: Nature Publishing Group.
Dylan Slack, Sophie Hilgard, Emily Jia, Sameer Singh, and Himabindu Lakkaraju. Fooling LIME and SHAP: Adversarial Attacks Post hoc Explanation Methods. February 2020. URL http::/axxiv..t/////////// http:/xwwwwxxxx../////oaa///////
[cs, stat],
02508
arXiv: 1911.02508.
Jakub M. Tomezak and Max Welling. VAE with a VampPrior. orXiv:1705.0711 [cs, stat], May 2017. URL TtR.htt::2 //trxv/rvvvvvivi org /bs 1705. 00000010.0000
07120
arXiv: 1705.07120.
J. Towns, T. Cockerill, M. Dahan, I. Foster, K. Gaither, A. Grimshaw, V. Hazlewood, S. Lathrop, D. Lifka, G. D. Peterson, R. Roskies, R. Scott, and N. Wilkins-Diehr, XSEDE: Accelerating scientific discovery. Computing in Science & Engineering, 16(5):62-74, October 2014. ISSN 1521-9615. doi: 20.1009/ URL ool. 2059/MCSE.2004.
Aki Vehtari, Andrew Gelman, and Jonah Gabry. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat Comput, 27(5):1413-1432 September 2017.
ISSN 1573-1375 doi: 10. 10.1007/s112220 122-0111122--1-999966 511222-016-9633
Matthew P. Wand, John T. Ormerod, Simone A. Padoan, and Rudolf Frühwirth. Mean Field Variational Bayes for Elaborate Distributions. Bayesian Anal., 6(4):847-900, December 1. ISSN 1936- 0975, 1931-6690. doi: 10. 10.1214/11-BA311 URL
ba/ 16/115 1339616546]
Sumio Watanabe. Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory. Journal of Machine Learning Research, 11 (Dec):3571-3594, 2010. ISSN ISSN 1533-7928, URL http:/ /wwp///www////////////////.. imlr. wwg.pppers/
v1 w11/watd http///ttaabeel0...at html
Michael R. Zhang, James Lucas, Geoffrey Hinton, and Jimmy Ba. Lookahead Optimizer: k steps forward, 1 step back.
[cs, stat], December 2019. URL http:// ///r//v arXiv: 1907.08611
org//bb 900.. 08610
12

---

Published as a conference paper at ICLR 2021
(iisonn)) pobe
Lineary struduned
Nonincaaly structuned
belgread Intenriate
Ta lagigaaa Sutten
buckyruual budaes
223
076
1/0
10.3
Figure S1: Factorization based on the logarithmic link function of Section 2.1 the synthetic dataset of Fig. 2
f(( x/n
f log(z/n) (3.54 0.02)
Poisson(1) noise
(3.54 0.02) x 102
x 10
Linearly structured
(4.45 0.03) 10a
(4.43 - 0.03)
10
Nonlinearly structured
(4.13 0.03)
x 10a
(4.13 + 0.03) x 10
Table 1: Model comparison using WAIC (#: standard error) for factorizatons of the synthetic data. Lower better.
Here we provide additional experiments. We note that the code for reproducing these experiments and those in the main manuscript can be found at github:mederratt/spnnf These supplemental examples can be found at Google Collaboratory Please refer to the notebooks therein, where the details behind hyperparameters and optimization. In general we did little tuning of the method beyond tuning the learning rate for stable inference.
also find
S1 COMPARING CHOICES OF f.g
In method, without refitting by using the WAIC. Here we provide an example of factorizations under different functions f. and compare the models. In Fig. ST. we performed factorization of the synthetic datasets of Fig. 2 using logarithmic link function of Section 2.1 Although the model is mis-specified, the key structure of the data is still exposed and irrelevant features are removed.
free choose functions Zi We evaluate models for predictive power
We then used WAIC to compare the use of the log link function versus the identity function. On the basis of predictive accuracy, the two models are similar as shown in Table [ so the method is not sensitive to this choice.
S2 SAMPLE SIZES
For systematic exploration of how sample size affects results, we used the nonlinearly generated synthetic dataset of Fig. and examined factorization as we varied N. Fig. [22 presents examples of
these factorizations using the standard HPF link functions of Eq. 7 and Fig. S3 presents factorizations
using the logarithmic link of Section 2.T
For U sufficiently large, the factorizations successfully remove the irrelevant background features. However, the structure of the factors is inconsistent as changes. Examining the correlation matrix of this dataset (Fig, S4] sheds light this behavior. Since the true generating data for this example is dense in every third feature, these features are highly correlated. Hence, without a spare substructure select, the factorization settles of the many sparse approximations the truly dense process.
SI

---

Published as a conference paper at ICLR 2021
k
1-
2-
U 500
2
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
5000
1
24
50000
2
13
14
18
19
24
Figure S2: Factorization under different data set sizes of the nonlinearly generated data of Fig.
50
k
1
2-
12
13
14
18
19
20
24
25
29
U = 500
0-
K
2
3
12
13
14
18
19
20
23
24
25
28
29
U=
k
1
2-
50000
2
10
12
13
14
15
17
18
19
20
23
24
25
26
28
29
Figure S3: Factorization under different data set sizes of the nonlinearly generated data of Fig. using the logarithmic link function.
S2

---

Published as a conference paper at ICLR 2021
0.25
0.20
0.15
22 22
27
Figure $4: Correlation of the nonlinear synthetic dataset features
S3


Reviews:
[{"review_id": "tz4cG9Mcmtg", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "** Description\n\nThis paper provides a new approach for finding a sparse encoding of count data matrices and hence automatically achieve feature selection.\n\n** Pros\n\nThe proposed technique is clearly efficient and practical.  It identifies a failing in traditional hierarchical Poisson matrix\nfactorisation (HPF) and proposes a solution.  This approach is tested on real world datasets where its usefulness is demonstrated.\n\n** Cons\n\nFor me the weakness of the approach is that the solution proposed appears very ad hoc with no probabilistic basis.  Particularly in\nmatrix factorisation problems where there is often no easy way to establish ground truth: this often leads to setting hyperparameters arbitrarily.  To be fair the authors point to ways of assessing predictive power which might help in this domain.  A less ad hoc model might allow a more principled approach to choosing hyperparameters.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "kZH2Siy_r6-", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This manuscript revisits the hierarchical Poisson matrix factorization (HPF) promoting sparsity in both the representation and the decoding function with a Horseshoe+ prior. Sparsity on both sides is put forward for interpretation purposes, and to create a column-clustering property. In addition, the proposed approach caters for non-linear decoding using a GAM model.\n\nThe main practical benefit compared to sparsifying priors (as in Gopalan et al. (2014)) is that the contributed method can be solved with automatic differentiation variational inference and hence stochastic solvers, which makes it in theory easier to scale, though this improvement is not demonstrated empirically in the manuscript.\n\nThe manuscript is well written: precise and well articulated. It develops well the theoretical point that sparsity in encoding and decoding is important. However, the practical value of the contribution is not strongly demonstrated. In the real-life application, on comorbidity data, the sparsity is demonstrated as expected, but the benefit compared to other approaches is difficult to gauge, whether it is to a non-sparse approach, or an approach based on sparsifying priors. The benefit of the GAM decoding is not clear.\n\nWhile there are good theoretical arguments for the model, they only partly convince in practical terms: a full analytics pipeline has made aspects to it, and the arguments might not be as important as they seem in practice. It could help to perform more empirical comparison, and to study more the contribution in the context of full analyses. The empirical demonstrations show that the model exhibit the properties that it was designed for: sparsity in the decoding. What are practical consequences of these properties in real-life applications?\n\nIn the bigger picture, it is unclear to see how this contribution positions itself in terms of practical benefits in the vast literature on latent factors with distangling approaches (distangling autoencoders, various matrix factorizations including NMF with different losses).\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "PCM4HTQSqUC", "reviewer": "AnonReviewer5", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "Summary:\nThe paper proposes a PMF variant that replaces the free-form representation for latent factors from explicit mapping from inputs, to improve interpretability. A reasonable probabilistic formulation with carefully selected priors is provided, but inference is carried out using generic tools, and the method is illustrated on artificial data and a simple comorbidity application.\n\nReasons for score:\nThe proposed model is interesting and well motivated, and the technical details regarding the choice of priors for encouraging sparsity are good and match current recommendations. The model itself is a bit counter-intuitive, explaining a generally desirable property (latent variables following a reasonably chosen prior distribution, free from additional computational constraints) as a limitation and proceeds to replace it with a simplified mapping from inputs, but as the mapping itself is well justified in terms of sparsity the overall construct still makes sense. The presentation angle is somewhat narrow and some connections are missed (e.g. the authors do not explain this as amortizing the inference for the latent variables, but explain the model in terms of encoders/decoders), but this is not a major issue.\n\nMy biggest issue with the paper concerns inference. The description is limited, only referring to a specific algorithm (ADVI) without specifying all details (mean-field vs full-rank approximation), and there is no discussion or analysis on how well it works. The authors do say that the algorithm converged fast, but do not show this in any experiment. More importantly, the convergence does not yet guarantee the approximation is good and there are well-known cases for which ADVI does not really work that well. Expanding both the discussion and empirical demonstration of this would be critical. Now you say \"one may use...\" and \"one can access...\" with references to specific techniques for evaluating the quality, which gives the impression you have not actually done that.  As I presume you implemented the model in Stan (no point in using ADVI if not; there are better stochastic VI methods around that are also easier to implement), would you be able to compare the inference results against HMC at least in some small-scale problem?\n\nPros:\n1. The idea is insightful and matches well the application needs.\n2. The priors match current literature on suggestions for sparsity-inducing priors.\n\nCons:\n1. Very limited coverage of inference, which is an important aspect even if carried out by an external software. Both theoretical and empirical evidence is missing, even though methods for evaluating the approximation quality are referred to.\n2. Figures 3 and 4 are pretty, but the spherical representation does not seem to add anything here and only results in the plots taking too much space while being slightly more difficult to read. More generally, the comorbidity example is a bit superficial and could have been developed a bit further. \n\nQuestions:\n1. How was ADVI applied?\n2. How did you check the approximation is good? Did you apply WAIC/PSIS-LOO or just say that it could be done?\n\nModifications after discussion:\nIncreased score by one since the revised paper clarifies the missing details on inference and also improves the motivation.", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors present a hierarchical factorization of the Poisson matrix and explain why sparcity in the encoder is important for interpretability. The reviewers appreciated the contribution of the paper and highlighted the advantage of such an approach for users. The authors have improved their initial version by adding more detail on inferences and experiments.  The decision is to accept the paper.
Decision: Accept (Poster)