Coupled Compound Poisson Factorization

0.4 MB·

by  Mehmet E. Basbug

Coupled Compound Poisson Factorization
MehmetE.Basbug [email protected]
PrincetonUniversity,35OldenSt.,Princeton,NJ08540USA
BarbaraE.Engelhardt [email protected]
PrincetonUniversity,35OldenSt.,Princeton,NJ08540USA
7
1 Abstract thatthedistributionofM dependsonlyonYobs.Themech-
0 anism is called not missing at random (NMAR) when the
We present a general framework, the cou-
2 distributionofM dependsonY .
pled compound Poisson factorization (CCPF), mis
n
to capture the missing-data mechanism in ex- When data are MCAR or MAR, the maximum likelihood
a
J tremely sparse data sets by coupling a hierar- estimatesforthedatageneratingmodelparametersdonot
9 chical Poisson factorization with an arbitrary change when M is taken into account (Little & Rubin,
data-generating model. We derive a stochas- 2014). In these cases, missing-data mechanism is said

] tic variational inference algorithm for the re- to be ignorable. When data are NMAR, however, it be-
G
sulting model and, as examples of our frame- comesmoreeffectiveforparameterestimationtoconsider
L
work, implement three different data-generating the joint likelihood of the data generating model and the
.
s models—amixturemodel,linearregression,and missingness-encoding model. In a Bayesian framework,
c
factor analysis—to robustly model non-random the missingness-encoding model may be represented as a
[
missingdatainthecontextofclustering,predic- probabilistic model; this model may be coupled with a

1 tion,andmatrixfactorization. Inallthreecases, arbitrary data generating model. When the missing-data
v
wetestourframeworkagainstmodelsthatignore mechanismisnon-ignorable,wehypothesizethatidentify-
8
themissing-datamechanismonlargescalestud- ingthemechanismcorrectlywillimproveinferenceinthe
5
0 ieswithnon-randommissingdata,andweshow jointdatageneratingmodel.
2 that explicitly modeling the missing-data mech-
Inthetheoreticalmachinelearningliterature,themissing-
0 anism substantially improves the quality of the
. dataproblemisoftendiscussedwithinthelimitedattribute
1 results,asmeasuredusingdataloglikelihoodon
observability framework (Birkendorf et al., 1998; Cesa-
0 aheld-outtestset.
Bianchietal.,2011;Hazan&Koren,2012;Kukliansky&
7
1 Shamir,2015)whereitisassumedthatthelearnercontrols
: which attributes it may observe. These models address a
v 1.Introduction
i variety of data generating models. Chechik et al. (2008)
X
The statistical theory of missing data developed by Little considered the situation where the learner does not con-
r & Rubin (2014) starts with an important distinction be- trol attribute observability in the linear regression setting.
a
tween a missing-data pattern and a missing-data mecha- Similarly,Hazanetal.(2015)proposedanon-probabilistic
nism. A missing-data pattern is an indicator matrix, M, algorithm for the classification problem under the low-
thatdescribeswhichmatrixvaluesaremissing. Amissing- rankassumptionandextremesparsity.Probabilisticmodels
data mechanism, on the other hand, captures the relation- for low-rank approximation of extremely sparse matrices
ship between the missing-data pattern and the data gener- areabundant,andincludeprobabilisticmatrixfactorization
atingmodel. WhenthedistributionofM,themissingness- (PMF) (Salakhutdinov & Mnih, 2011), non-negative ma-
encoding model, does not depend on the observed data, trix factorization (NMF) (Lee & Seung, 1999), and their
Y , or the missing data, Y , the missing-data mech- variants.
obs mis
anism is characterized as missing completely at random
On the applied side, a motivating missing-data problem
(MCAR). In contrast, missing at random (MAR) indicates
with extreme sparsity is collaborative filtering—creating
a predictive ranking of items for each user given observa-
tions of users’ preferences—where capturing the missing-
Preliminarywork datapatterniscrucialforaccurateranking(Huetal.,2008).
CoupledCompoundPoissonFactorization
Table1.Sixcommonadditiveexponentialdispersionmodels. Gaussian,gamma,inverseGaussian,Poisson,binomial,andnegative
binomialdistributionsinadditiveEDMform.Forthegammadistribution,aandbrefertotheshapeandrateparameters,respectively.
DISTRIBUTION θ κ Ψ(θ) h(x,κ)
GAUSSIAN N(µ,σ2) µ/σ2 σ2 θ2/2 −x2 − 1log(2πκ)
2κ 2
GAMMA Ga(a,b) −b a −log(−θ) (κ−1)logx−logΓ(κ)
√ √
INV. GAUSSIAN IG(µ,λ) − λ λ − −2θ −κ2 +log(κ)− 1log(2πx3)
2µ2 2x 2
POISSON Po(λ) logλ 1 eθ xlogκ−logΓ(x+1)
BINOMIAL Bi(r,p) log( p ) r log(1+eθ) log(cid:0)κ(cid:1)
1−p x
NEG. BINOMIAL NB(r,p) logp r −log(1−eθ) log(cid:0)x+κ−1(cid:1)
x
State-of-the-artprobabilisticcollaborativefilteringmodels havelargelydriventhedevelopmentofcollaborativefilter-
arebasedoneitherPMForPoissonfactorization(PF),the ingmodelsthatincludeNMARdata;however,theproblem
probabilistic counterpart of NMF (Cemgil, 2009). For in- ofNMARdataexistinabroadrangeofanalytictasksthat
stance, weighted matrix factorization (WMF) (Hu et al., have been underdeveloped up to now. Motivated by these
2008) is a PMF model fit to binarized data with a het- two observations, we develop an explicit NMAR missing
eroscedastic variance term. The exposure matrix factor- data mechanistic model coupled with an arbitrary genera-
ization model (ExpoMF) uses a WMF model conditioned tive model framework. We show the benefits of including
on a Bernouilli exposure matrix to capture the missing- an explicit missing-data mechanistic model on three spe-
data pattern (Liang et al., 2015). Another successful im- cificgenerativemodeltasks:mixturemodelsforclustering,
plicitfeedbackmodelishierarchicalPoissonfactorization linearregressionforprediction,andlatentfactormodelsfor
(HPF)(Gopalanetal.,2013). BuildingontheHPF,hierar- matrixfactorization.
chicalcompoundPoissonfactorization(HCPF)usesPFto
We start with data generating models with an additive ex-
encode the missing-data pattern and extends the PF struc-
ponential dispersion model output. This large collection
ture with a general additive exponential dispersion model
of models includes Gaussian, gamma, binomial mixture
togeneratedata(Basbug&Engelhardt,2016). HCPFisa
models, PMF, HPF, and linear regression models among
flexiblemodelthatcanbeusedtofactorizecontinuousreal-
others. Weexploretherelationshipbetweenthedatagen-
valuedornon-negativedataaswellasnon-negativediscrete
erating model and the missingness-encoding model, and
data.
we identify the need for a missing-data mechanism that
The increasing popularity of PF-based methods in collab- can capture a heteroscedastic relationship between obser-
orative filtering can be attributed to two factors. First, the vationsandthemissing-datapattern. Morespecifically,we
Gamma-PoissondistributionsofPFareconjugatewhereas empirically show that the variance of an observation y
ij
the Gaussian-Bernoulli distributions of PMF models are is a function of probability of missingness of that obser-
not. Thus, probabilistic inference in PF-based models for vation, Pr(M ), in three large data sets. To address this
ij
large scale problems is more straightforward and compu- issue, we propose the coupled compound Poisson factor-
tationally tractable (Gopalan et al., 2013; Basbug & En- ization (CCPF) framework as a missing-data mechanism
gelhardt, 2016). Second, in HPF and HCPF, hierarchical for NMAR data with extreme sparsity. We prove that the
structure models the user activity and the item popularity CCPF model reduces to the data generating model when
inanaturalway. Inparticular,Gopalanetal.(2013)show themissingness-encodingmodelisignorable. Thefirstim-
thattheheavytailGammapriorsinHPFaccuratelycapture plication of this result is that it is sufficient to update the
the user behavior and item popularity by using posterior parameters of the data-generating model using the non-
predictivechecks. missing entries only. The second implication is that the
statistician is able to describe the data generating model
Anequallyimportantbutmorecomplexaspectofimplicit
ignoring the possible impact of the missingness-encoding
feedback models is the missing-data mechanism. Marlin
model on heteroscedasticity. When heteroscedasticity ex-
& Zemel (2009) give ample empirical evidence motivat-
istsinthedata,themissingness-encodingmodelwithinthe
ingtheneedfor,andshowingtheperformancebenefitsof,
CCPF framework will accurately capture and control for
explicitly modeling NMAR missing-data mechanisms in
thiscomplexityintheanalytictask.
collaborative filtering models. Recommendation systems
CoupledCompoundPoissonFactorization
3.50 3.00
amazon
movielens
3.00 2.50
yelp
2.50 2.00
e
Varianc 2.00 Scaling 1.50
d 1.50 d 1.00
e e
aliz ect
m p
or 1.00 Ex 0.50
N
0.50 0.00
0.00 -0.50
-0.50 -1.00
-0.20 0.00 0.20 0.40 0.60 0.80 1.00 1.20 0.00 0.20 0.40 0.60 0.80 1.00
Probabilityofnon-missingness Probabilityofnon-missingness
(a) (b)
Figure1.Exploratoryanalysisforheteroscedasticity.a)Scatterplotsofthenormalizedvarianceoftheresidualsoftestdatasetsafter
fittingaPMFmodeltonon-missingentrieswithrespecttoqunatilesoftheprobabilityofnon-missingnesscalculatedbyfittinganHPF
modeltothebinarizedfullmatrixforamazon(blue)movielens(green)andyelp(red)datasets,b)expectedscalingversusprobability
ofnon-missingnessforctakingvalues0.5(solid), 0.125(dashed), −0.125(dash-dot), −0.5(dotted)orangeforexponentiallinkage
functionandpurpleforlinearlinkagefunction.
2.AdditiveExponentialDispersionModels thenX =(cid:80) X ∼p (x;θ,(cid:80) κ ).
+ j j Ψ j j
We start by describing the family of additive exponential
We will exploit this property to provide theoretical jus-
dispersion models and its defining characteristic, additiv-
tification of compounding in our coupling framework.
ity. Exponential dispersion models (EDM) have been pri-
Moreimportantly,inferenceincompoundPoissonadditive
marily used as error distributions for generalized linear
EDMs is relatively straightforward due to the availability
models (?)for a comprehensive treatment of the theory of
of the conditional and marginal densities (Basbug & En-
EDMs,see][]jorgensen1997theory.AdditiveEDMs,asub-
gelhardt,2016).
familyofEDMs,includesGaussian,gamma,inverseGaus-
sian, Poisson, binomial, and negative binomial distribu-
3.CouplingFramework
tions, among others (see Table 1). Following Basbug &
Engelhardt(2016),wedefineadditiveEDMsasfollows.
Lety ∼p (θ ,κ)Pr(θ |ϑ)beanarbitrarydatagener-
ij Ψ ij ij
(cid:8)Dpe(fiΨn,θit,iκo)n|θ1∈. ΘA=dfaommi(lyΨ)o⊆f Rd,iκstr∈ibRut+io+n(cid:9)s iFsΨcalle=d adtiisnpgermsioodneplasruacmhettheartκp,Ψa(θnidj,Pκr)(iθsijan|adϑd)itiivseaEcDoMnjuwgaitthe
anadditiveexponentialdispersionmodelif prior hierarchical model with parameters ϑ. Examples of
p (x)=exp(xθ−κΨ(θ)+h(x,κ)) such data generating models include Gaussian, Poisson,
(Ψ,θ,κ)
gamma,andbinomialmixturemodels,PMF,HPF,andlin-
whereθisthenaturalparameter,κisthedispersionparam-
earregressionmodels.
eter,Ψ(θ)isthebaselog-partitionfunction,andh(x,κ)is
thebasemeasure. First,weinvestigatetherelationshipbetweenthemissing-
datapatternandthedata-generatingmodel. WefitaPMF
Withthedefinitionabove, weseethatthesumofadditive model to three different data sets with MCAR assump-
EDMs sharing Ψ and θ has the same base log-partition tion, i.e., only the non-missing entries are sampled dur-
function and natural parameter. The following theorem ing training. We also fit an HPF model to the binarized
makesthisstatementconcrete: full matrix for each of the data sets. For a test set of
Theorem 1. (Jorgensen, 1997) Let X ...X be a se- non-missingentries,wecalculatetheresidualsunderPMF
1 M
quence of additive EDMs such that X ∼ p (x;θ,κ ), (y − yˆ ) and the probability of non-missingness under
j Ψ j ij ij
CoupledCompoundPoissonFactorization
Table2.DatasetstatisticsNumberofsamples,attributes,non-missingentries,andtheratioofmissingentriestothetotalnumberof
entries(sparsity),themaximumandminimumobservationvalues.
DATASET #SAMPLES #ATTRIBUTES SPARSITY #NON-MISSING MIN MAX
AMAZON 256,059 74,258 0.999970 568,454 1 5
MOVIELENS 247,753 9,732 0.990723 22,369,101 1 10
NETFLIX 480,189 17,770 0.988224 100,483,024 1 5
YELP 552,339 77,079 0.999948 2,225,213 1 5
WORDPRESS 86,661 78,754 0.999915 581,508 1 1013
GEUVADIS 9,358 462 0.462121 2,325,461 -19.3 12.6
HPF, Pr(M = 0). We discretize data into 100 equal- • Foreachrowi=1,...,C
ij I
sizedbinsbasedonsamplequantilesofprobabilityofnon-
1. Sampleactivityr ∼Ga(ρ,(cid:37))
missingness, Pr(M = 0). For each bin, we calculate i
ij
2. For each component k, sample factor weight
the variance of the residuals under PMF, Var(y −yˆ ).
ij ij
u ∼Ga(η,r )
Fig. 1(a) shows the scatter plot of variance within each ik i
bin and the quantiles of probability of non-missingness
• Foreachcolumnj =1,...,C
J
for amazon, movielens and yelp data sets. We clearly
see that there is a linear/exponential relationship between 1. Samplepopularitywj ∼Ga(υ,ν)
Var(yij−yˆij)andPr(Mij =0)inallthreedatasets. We 2. For each component k, sample factor weight
alsotestforheteroscedasticitybyfittingadoublegeneral- v ∼Ga(ζ,w )
jk j
ized linear model (DGLM) (Dunn & Smyth, 2012) to the
residuals using the probability of non-missingness as the • Foreachiandj
regressor. We conclude that there is strong heteroscedas-
1. Sampleinteractionvariable
ticityinallthreedatasets(maxp≤6×10−17). (cid:80)
n ∼Po(Λ = u v )
ij ij k ik jk
To address the heteroscedastic structure in data analysis 2. ifn is0,theny ismissing
ij ij
with missing observations, we present a framework for 3. else
NMAR missing-data mechanism that can handle large-
(a) sample natural parameter θ from the data
ij
scale, extremely sparse data with two additional desirable
generatingmodelPr(θ |ϑ)
ij
properties: i) compatibility with the low-rank assumption
(b) sampleobservation
forthemissing-datapattern,M,ii)convergencetotheho-
y ∼p (θ ,φ(n )κ).
moscedastic model y ∼ p (θ ,κ) when y is almost ij Ψ ij ij
ij Ψ ij ij
surelymissing(Pr(M =0)→0).
ij
First, we summarize the model assumptions behind HPF
Let the missingness-encoding model be an HPF (Gopalan
when used to capture the missing-data pattern. HPF is
et al., 2013; Basbug & Engelhardt, 2016). In HPF, an en-
anextensionofPoissonfactorization(PF)whichisessen-
try, y , is missing when a draw from the model, n , is
ij ij tially a non-negative matrix factorization model. Imagine
zero; otherwise, y = n . HCPF has a more flexible ob-
ij ij we have a movie ratings data set, for which an entry is
servation model, i.e., y ∼ p (θ,n κ), when n (cid:54)= 0.
ij Ψ ij ij missingwhenauserhasnotratedamovie. InPFsetting,
However, HCPF assumes that the natural parameter θ is
ij the missing-data pattern is approximated with an interac-
thesameforallobservations;therefore,thedata-generating
tion matrix where each entry is a latent Poisson random
model is a fixed distribution and not a complex structure variable (n for the ith user and jth movie). The inter-
ij
such as a mixture model or a regression model. In our action variable, n , is the sum of k interaction contribu-
ij
framework,weassumeanarbitrarydatageneratingmodel tions, n , which are themselves are Poisson distributed
ijk
yij ∼ pΨ(θij,κ)Pr(θij | ϑ) whose output is an additive (n ∼ Po(u v )). Thelowrankassumptionk << C
ijk ik jk I
EDM.
impliesthatthereexistlatentgroupsofusersandthefactor
We couple the missingness-encoding and the data gener- weightuik modelsthemembershipofitokth group. Dif-
atingmodelsbyscalingthedispersionparameterκwitha ferentthanPF,HPFhasanotherrandomvariablerimodel-
linkage function φ(n ). The full generative model of the inghowactivetheuseriis(i.e.howmanymoviesuserihas
ij
coupledcompoundPoissonfactorization(CCPF)isasfol- rated). Ifauserisactive,thenthefactorweightuikandthe
lows: interactioncontributionnijk isadjustedaccordingly. Sim-
ilarly,thereexistlatentgroupofmovies(k <<C )where
J
thefactorweightv modelsthemembershipofmoviejto
jk
kth genre. Thepopularityvariablew controlshowmany
j
CoupledCompoundPoissonFactorization
Algorithm1SVIforCCPF tion on movie ratings where a user group may dislike a
Initialize: Hyperparametersc,η,ζ,ρ,(cid:37),υ,ν movie genre and down-vote. Furthermore, the latent fac-
Setthelearningratest =t =t andτ =τ =t−ξ tors of movies or users in terms interaction can be vastly
i j 0 i j 0
repeat different than the latent factors in terms preference. This
Sampleanobservationy uniformly bringsustothedatageneratingmodel. InCCPF,wehave
ij
CalculateΛ =(cid:80) αuαv /(βuβv ) the ability to choose an arbitrary data-generating model
ij k ik jk ik jk
ifyij ismissingthen Pr(θij | ϑ)forthemovieratings,perhapsanothermatrix
q(n )=δ factorization model. We model the relationship between
ij 0
else the missingness-encoding model and the data-generating
CalculatethesufficientstatisticE[Ψ(θij)]fromthe modelviathelinkagefunctionφ(nij).
data-generatingmodel
Second, we note that the probability distribution of the
Update the local variational parameters of the
missing-datapattern,M,isgivenby
missingness-encodingmodel
Pr(M =1)=Pr(y ismissing)=e−Λij.
ij ij
q(n =n)∝exp{−κφ(n)E[Ψ(θ )]
ij ij
Λn In the generative model, we have the linkage function
+h(y ,φ(n)κ)} ij
ij n! φ(nij)scalingthedispersionparameter,κ.Theconditional
ϕ ∝exp{Ψ(αu)−logβu meanandvarianceofyij aregivenby
ijk ik ik
+Ψ(αv )−logβv (cid:9) E[y |n ]=φ(n )κΨ(cid:48)(θ )
jk jk ij ij ij ij
Var(y |n )=φ(n )κΨ(cid:48)(cid:48)(θ ).
Update the variational parameters of the data gen- ij ij ij ij
eratingmodelusingE[φ(n )]
ij
endif ThedatageneratingmodelwithGaussianobservationsisa
Updateglobalvariationalparameters specialcaseofthisframeworkbecauseθij isalsoadditive.
Thus, we can use separate scaling functions φ (n ) and
µ ij
αr =(1−τ )αr+τ (ρ+Kη) φ (n )thatresultin
i i i i σ ij
αw =(1−τ )αw+τ (υ+Kζ)
j j j j E[y |n ]=φ (n )µ
ij ij µ ij ij
αu =(1−τ )αu +τ (η+C E[n ]ϕ )
ik i ik i I ij ijk Var(y |n )=φ (n )σ2.
ij ij σ ij
αv =(1−τ )αv +τ (ζ+C E[n ]ϕ )
jk j jk j U ij ijk
(cid:32) (cid:33)
βr =(1−τ )βr+τ (cid:37)+(cid:88)αiuk Ignorable missing-data mechanism: When φ(nij) =
i i i i βu φI(nij) = 1,thedatageneratingmodelisdecoupledfrom
k ik the missingness-encoding model. This is equivalent to
(cid:32) (cid:33)
βw =(1−τ )βw+τ ν+(cid:88)αjvk trainingthetwomodelsseparately.
j j j j βv
k jk
(cid:32) (cid:33) Linear missing-data mechanism: To obtain a linear
βu =(1−τ )βu +τ αir +C αjvk relationship between the probability of missingness and
ik i ik i βr Iβv the dispersion parameter, we set the linkage function as
i jk
(cid:32)αw αu (cid:33) φ(nij) = φL(nij) = 1 − c + c(−1)nij+1. The expec-
βjvk =(1−τj)βjvk+τj βwj +CUβiuk tation of φL(nij) under the zero truncated Poisson (ZTP)
j ik distributionnij ∼ZTP(Λij)isgivenby
Updateuserlearningrateti =ti+1andτi =t−i ξ E(cid:2)φL(nij)=1−c+c(−1)nij+1(cid:3)=1−c+ce−Λij.
Updateitemlearningratet =t +1andτ =t−ξ
j j j j In Fig. 1(b), the purple lines indicatethe expected scaling
(Optional)Updatehyperparameterscandκ
as the probability of non-missingness changes between 0
untilvalidationsetloglikelihoodconverges
and1.Asseeninthefigure,whenc>0,thereisaninverse
linearrelationshipbetweenthetwo. Inotherwords,asthe
probabilityofnon-missingnessincreasestheexpecteddis-
times the jth movie has been rated. The key point is that persiondecreaseslinearly. Thisiscompatiblewiththeem-
theinteractionsbetweengroupsofusersandthegroupsof piricalfindingsinFig.1(a).
movies are non-negative. A group of users can only ‘not
interact’ with a movie genre, hence n cannot be neg- Exponential missing-data mechanism: Another possi-
ijk
ative. This is different than performing matrix factoriza- ble linkage function is exponential. Setting φ(n ) =
ij
CoupledCompoundPoissonFactorization
φE(n ) = 1 − c + cn with c > 0 implies that, as Pr(θ |ϑ)canbeestimatedusingthenon-missingentries
ij ij ij
the probability of non-missingness increases, we expect only,aswouldbedoneforhomoscedasticmodel.
a greater dispersion. Fixing c = 1, we get a standard
For computationally tractable learning, we use stochastic
compoundPoissonadditiveEDMmodel(Basbug&Engel-
variational inference (SVI) (Hoffman et al., 2013), which
hardt,2016). Whenc < 0,thereisaninverserelationship
minimizes the lower bound on the expected posterior log
betweentheprobabilityofnon-missingnessandthedisper-
likelihoodunderavariationaldistribution. Themeanfield
sion. The expectation of φE(n ) under the ZTP is given
ij variationaldistributionforCCPFisgivenby
by
E(cid:2)φE(nij)=1+c(nij −1)(cid:3)=1−c+ceeΛΛiijjΛ−ij1. q(ri |αir,βir)q(uik |αiuk,βiuk)q(wj |αjw,βjw)
q(v |αv ,βv )q(nˆ |ϕ )q(n )q(θ )q(ϑ),
In Fig. 1(b), the orange lines show the expected scaling jk jk jk ij ij ij ij
withrespecttotheprobabilityofnon-missingnessfordif-
ferent values of c. When c > 0, the expected scaling in-
where the distributions of the variational approximation
creasesexponentiallyastheprobabilityofnon-missingness
are as in the generative distribution. For gamma distri-
increases. Expectedscalingdecreasesexponentially,when
butions, α is the shape and β is the rate parameter. We
c < 0. Since dispersion is always positive, one needs to
makeuseofthemultinomialrepresentationofPoissonfac-
be careful in choosing c. Exponential decay with a small
tors as in (Cemgil, 2009; Gopalan et al., 2013; Basbug
c might be a better choice than the linear relationship for
& Engelhardt, 2016). To update the variational parame-
certaindatasetsasseeninFig.1(a).
ters of the data generating model, we need the sufficient
AnothertakeawayfromFig.1(b)isthattheexpectedscal- statisticE[φ(nij)]. Similarly,weneedthesufficientstatis-
ingconvergesto1.0astheprobabilityofnon-missingness tics E[Ψ(θij)] to update the variational parameters of the
becomes 0. With the following theorem, we provide a missingness-encoding model. The SVI algorithm for the
strongerresult. generalCCPFframeworkissummarizedinAlgorithm1.
Theorem 2. Let N be a zero-truncated Poisson random
variablewithparameter Λ, andlet Y ...Y bei.i.d. ad- 4.Experiments
1 N
ditive EDM random variables such that Y ∼ p (θ,cκ),
j Ψ
WeapplytheCCPFframeworktothreeclassesofdatagen-
where 0 < c < 1. Let Y ∼ p (θ,(1 − c)κ) be an-
0 Ψ
eratingmodels: mixturemodels, matrix factorization, and
otherrandomvariablewiththesamelog-partitionfunction
linearregression. Foreachmodel,wecomparetoCCPFto
Ψ and the natural parameter θ. Define the linear linkage
. data-generatingmodelsthatignorethemissing-datamech-
functionasφL(N)=1−c+c(−1)N+1andtheexponen-
. anism. In each case, the CCPF approach outperforms the
tial linkage function as φE(N) = 1−c+cN. Then the
probabilitydensitiesofYL =. Y + (cid:80)N (−1)n+1Y comparisonsintermsoftheloglikelihoodofheld-outnon-
andYE =. (cid:80)N Y ar+egiven0by n=1 n missingentries.
+ n=0 n
Mixturemodels,linearregression,andmatrixfactorization
YL ∼p (θ,φL(N)κ) arethreeusefuldatageneratingmodelsfortheanalysisof
+ Ψ
YE ∼p (θ,φE(N)κ). high-dimensional data. Mixture models represent the par-
+ Ψ titioning of observations into subgroups, and are used for
Furthermore,bothYL andYE convergetoY ∼ p (θ,κ) exploratorydataanalysis.Matrixfactorizationdecomposes
+ + Ψ
indistributionasΛgoestozero;thatis a matrix of observations into two lower dimensional ma-
trices, allowing each observation from the original matrix
YL →D Y as Λ→0 to be represented as the weighted linear combination of
+
a lower dimensional space. Linear regression models the
Y+E →D Y as Λ→0. relationship between regressors and observations; in our
case,theobservationsaremultivariatewithalargenumber
Theorem 2 implies that the distribution of an observation ofdimensions.1 Allthreeofthesemodelsmaybeapplied
y thatismissingalmostsurely(Λ → 0)convergestothe tothesamedata,dependingontheobjectivesoftheanaly-
ij
homoscedastic observation model, i.e., p (θ ,κ). From sis. ThegoalofCCPFistoaccountformissingdataunder
Ψ ij
the model building aspect, this is a useful property, since variousmodelingparadigms.
we can design the data-generating model independently
1NotethatCCPFisonlyapplicabletoalinearregressionset-
andsimplypluginthemissingness-encodingmodeltocap-
tingswherewehavefullaccesstoregressors,butmultivariateob-
tureanysuspectedheteroscedasticity. Anothercorollaryof servationsaremissing. Itisnotapplicabletodatawithmissing
thistheoremisthattheparametersofdatageneratingmodel covariates.
CoupledCompoundPoissonFactorization
Table3.Non-missingtestsetloglikelihoodformixturemodels. Loglikelihoodpernon-missingtestentryforclusteringperformed
acrosssamplesandattributes.GMMandPMMstandforGaussianmixturemodelandPoissonmixturemodel,respectively.CCPF-GMM
andCCPF-PMMarethevariantsofthesetwomodelsinourframework.
AXIS SAMPLING MOVIELENS AMAZON NETFLIX YELP GEUVADIS
GMM SAMPLE NONZERO -2.082 -1.549 -1.427 -1.709 -2.139
GMM ATTRIBUTE NONZERO -2.064 -1.667 -1.491 -1.637 -2.173
CCPF-GMM SAMPLE FULL -2.121 -1.487 -1.419 -1.708 -2.012
CCPF-GMM ATTRIBUTE FULL -2.102 -1.669 -1.480 -1.615 -2.083
PMM SAMPLE NONZERO -3.554 -1.793 -3.682 -1.912 N/A
PMM ATTRIBUTE NONZERO -3.674 -1.886 -1.918 -1.793 N/A
CCPF-PMM SAMPLE FULL -3.713 -1.789 -3.315 -1.889 N/A
CCPF-PMM ATTRIBUTE FULL -3.655 -1.913 -2.001 -1.801 N/A
4.1.Datasets 2016). WesetN usingtheexpectedrangeofλ andy
tr ij ij
aswellasthemaximumlikelihoodestimatesofE[Ψ(θ )]
ij
We analyze four user behavior data sets and one gene ex-
and κ under the assumption that θ is fixed. We have an
ij
pression data set; Table 2 outlines the characteristics of
inverseGammaprioronκwithshapeparameter1.01and
each. Theuserbehaviordataincludesmultipleratingsdata
scaleparameter1.0. WenotethattheinverseGammaisthe
sets: amazon contains ratings of fine food (McAuley &
conjugatepriortoκfortheexponentialdispersionmodels
Leskovec, 2013), netflix consists of movie ratings (Bell &
under the saddle-point approximation (Jorgensen, 1997).
Koren,2007),andyelpcomprisesvenueratings;ratingsfor
Additionally,wehaveaGaussianpriorwithzeromeanand
allthreeofthesedatasetsrangefrom1to5. Wealsocon-
0.1standarddeviationonthelinkageparameterc. Wethen
sider movielens, another movie data set with ratings rang-
calculate the MAP estimates of κ and c using stochastic
ingfrom1to10. Inthisdata,eachmoviehas1129corre-
gradientdescentwithsmoothedgradientsbytakingtheav-
spondingexogenousvariablesrepresentingtheassociation
erageof1000gradients.
with a predetermined tag. The tag weights are calculated
fromthecrowdsourcedtag-movieassociationdata(Harper
4.3.MixtureModelResults
&Konstan,2015). Behaviordataisnotalwaysbounded—
wordpressisasocialmediainteractiondatasetofusersand We first consider Gaussian mixture models (GMM) with
blogswheretheresponseisthenumberoflikesauserhad spherical Gaussian priors and Poisson mixture models
for a given blog. The geuvadis data shows that CCPF is (PMM) with Gamma priors as our data-generating mod-
applicablebeyonddiscreteuserbehaviordata; itisagene els. We derive SVI algorithms modified to accommodate
expressiondatasetof9,358genesforeachof462individ- coupling(seeAppendix).
uals.
We fit GMM, PMM, CCPF-GMM and CCPF-PMM to
movielens,amazon,netflixandyelpdatasets. Weperform
4.2.ExperimentalDetails
clusteringacrossusers(samples)anditems(attributes)sep-
We held out 20% and 1% of the non-missing entries for arately. In the continuous gene expression data set (geu-
testing (Ytest) and validation, respectively. Test log like- vadis),wefitGMMandCCPF-GMMtoidentifyclusterof
obs
lihood of the non-missing entries (L ) under CCPF is individuals (samples) and genes (attributes). We first take
NM
calculatedas thelogofthegeneexpressionlevels. Wecomparemodels
performanceintermsoftestsetloglikelihood.
(cid:88) (cid:88)Ntr
LNM = log pΨ(yitjest;θij,φ(n)κ)ZTP(n|Λij). In Table 3, we observe that CCPF-GMM outperforms
Ytest n=1 GMM in amazon, netflix, yelp and geuvadis for cluster-
obs
ingbothsamplesandattributes. Onlyinmovielens,wesee
To set the hyper parameters of the missingness-encoding
GMM beating CCPF-GMM. In PMM comparisons, there
modelwefollowedthemethodpresentedforHCPF(Bas-
isnoconsistentpattern. Thiscanbeexplainedbyourpre-
bug & Engelhardt, 2016). We fix K = 160, ξ = 0.7 and
viousobservationthattheimpactofmissingness-encoding
t = 10,000 after an empirical study on small data sets.
0 model on the data-generating model is significant for the
We estimate the sparsity level from the number of non-
dispersionparameterbutnotsignificantforthemeanterm.
missing entries and from the sparsity level we calculated
In the case of GMM, the scaling of dispersion parameter
E[n ]. ForheavytailGammapriors,wesetν = (cid:37) = 0.1
ij only affects the variance term. In PMM, φ(n ) scales
(cid:112) ij
andυ = ρ = 0.01. Wethensetη = ρ E[nij]/K/(cid:37)and λ whichisthemeanandthevarianceoftheobservation;
(cid:112) ij
ζ = υ E[n ]/K/ν as done in (Basbug & Engelhardt,
ij
CoupledCompoundPoissonFactorization
Table4.Non-missing test set log likelihood for factorization models. PMF stands for the probabilistic matrix factorization, HPF
standsforhierarchicalPoissonfactorization.CCPF-PMFandCCPF-HPFarethevariantsofthesetwomodelsinourframework.
SAMPLING MOVIELENS AMAZON NETFLIX YELP WORDPRESS
PMF NONZERO -2.089 -1.703 -1.421 -1.952 -3.379
HPF NONZERO -2.104 -1.868 -1.634 -1.881 -2.146
CCPF-PMF FULL -2.009 -1.630 -1.383 -1.734 -3.046
CCPF-HPF FULL -2.013 -1.835 -1.693 -1.878 -1.922
mean ratings across users. SVI for the resulting coupled
Table5.Comparison for regression models. Log likelihood,
modelisgivenintheAppendix.
RMSE and R2 values for hierarchical linear regression models
onmovielensdataset.CCPF-Regr.isthevariantoftheregression We fit the regression model as well as its coupled version
modelinourframework. to movielens data set where the attribute vector is the tag
SAMPLING TLL RMSE R2 weightsdiscussedearlier. Thegoalistocaptureuserpref-
REGRESSION NONZERO -1.921 1.731 0.333 erencesoverarbitrarytagsandusethisinformationwithin
CCPF-Regr. FULL -1.907 1.694 0.360 the collaborative filtering setting. In Table 5, we observe
that coupling improves the regression performance. We
therefore,couplingismoreeffectiveinGMM. alsonotethattheregressionapproachachievesthehighest
testloglikelihoodamongmixturemodelandmatrixfactor-
4.4.MatrixFactorizationResults izationapproaches.
Next, we considered two major matrix factorization mod-
5.Conclusion
els. PMF is the probabilistic counterpart of regularized
SVDwherethepenaltytermsrelatetothesphericalGaus- In this work, we present the coupled compound Pois-
sianpriorsonfactors(Salakhutdinov&Mnih,2011). HPF son factorization (CCPF) that models the missing-data
is another probabilistic matrix factorization model where mechanism in extremely sparse data sets by coupling
the factor contributions are non-negative (Gopalan et al., a missingness-encoding model with an arbitrary data-
2013). Previously, we utilized HPF for the missingness- generating model. We derive stochastic variational infer-
encoding model; however, HPF can also be trained only ence algorithm for our CCPF models and, as examples of
on the non-missing entries. We derive SVI algorithm our framework, implement instances of a mixture model,
for PMF and HPF accommodating the coupling from the linear regression, and matrix factorization. We compare
missingness-encodingmodel(seeAppendix). ourmodelwiththedatageneratingmodelsthatignoresthe
We fit PMF, HPF, CCPF-PMF and CCPF-HPF to movie- missing-data mechanism on large scale studies and show
lens,amazon,netflix,yelpandwordpressdatasets. Table4 that explicitly modeling the missing-data mechanism sub-
stantiallyimprovestestloglikelihoodandothermetricsrel-
summarizes the comparison of models in terms of test set
evanttotheanalysisofinterest.
loglikelihood. Inratingsdatasets,CCPF-PMFisthebest
performing algorithm with a clear margin in most cases.
Similartothemixturemodelanalysis, weobserveamore
substantialimprovementinPMFtoCCPF-PMFtransition
than HPF to CCPF-HPF transition. Social media activ-
itydataset,wordpress,exhibitsadifferentcharacteristic—
PoissonmodelsoutperformPMFandCCPF-PMF.Thiscan
beattributedtothehighlydispersedPoisson-likeresponse
distribution.
4.5.LinearRegressionResults
Finally, we considered a hierarchical linear regression as
the data-generating model. Let x be the attribute vec-
j
tor for item j. We describe the data-generating model by
y ∼ N(b Tx ,σ2)whereeachuser, i, hasauniqueco-
ij i j
efficientvectorb . WehaveaGaussianprioronthecoeffi-
i
cientvectorsthatissettothemaximumlikelihoodestimate
of the coefficient vector when the response is fixed to the
CoupledCompoundPoissonFactorization
References Jorgensen, Bent. The theory of dispersion models. CRC
Press,1997.
Basbug, Mehmet E and Engelhardt, Barbara E. Hierar-
chical compound Poisson factorization. Proceedings of Kukliansky, Doron and Shamir, Ohad. Attribute efficient
theInternationalConferenceonMachineLearning, pp. linear regression with distribution-dependent sampling.
17951803, July 2016. URL http://arxiv.org/ InProceedingsofthe32ndInternationalConferenceon
abs/1604.03853. MachineLearning(ICML-15),pp.153–161,2015.
Bell, Robert M and Koren, Yehuda. Lessons from the
Lee,DanielDandSeung,HSebastian. Learningtheparts
netflix prize challenge. ACM SIGKDD Explorations ofobjectsbynon-negativematrixfactorization. Nature,
Newsletter,9(2):75–79,2007.
401(6755):788–791,1999.
Birkendorf, Andreas, Dichterman, Eli, Jackson, Jef-
Liang, Dawen, Charlin, Laurent, McInerney, James, and
frey, Klasner, Norbert, and Simon, Hans Ulrich.
Blei, David M. Modeling user exposure in recommen-
On restricted-focus-of-attention learnability of boolean
dation. arXivpreprintarXiv:1510.07025,2015.
functions. MachineLearning,30(1):89–123,1998.
Little,RoderickJAandRubin,DonaldB. Statisticalanal-
Cemgil, Ali Taylan. Bayesian inference for nonnegative
ysiswithmissingdata. JohnWiley&Sons,2014.
matrixfactorisationmodels. ComputationalIntelligence
andNeuroscience,2009,2009. Marlin,BenjaminMandZemel,RichardS. Collaborative
prediction and ranking with non-random missing data.
Cesa-Bianchi, Nicolo, Shalev-Shwartz, Shai, and Shamir,
InProceedingsofthethirdACMconferenceonRecom-
Ohad. Efficient learning with partially observed at-
mendersystems,pp.5–12.ACM,2009.
tributes. The Journal of Machine Learning Research,
12:2857–2878,2011. McAuley,JulianJohnandLeskovec,Jure. Fromamateurs
to connoisseurs: modeling the evolution of user exper-
Chechik, Gal, Heitz, Geremy, Elidan, Gal, Abbeel, Pieter,
tisethroughonlinereviews. InProceedingsofthe22nd
and Koller, Daphne. Max-margin classification of data
international conference on World Wide Web, pp. 897–
withabsentfeatures. TheJournalofMachineLearning
908,2013.
Research,9:1–21,2008.
Salakhutdinov,RuslanandMnih,Andriy.Probabilisticma-
Dunn,PeterKandSmyth,GordonK. dglm: Doublegen-
eralizedlinearmodels. Rpackageversion,1(2),2012. trixfactorization. Citeseer,2011.
Gopalan,Prem,Hofman,JakeM,andBlei,DavidM.Scal-
able recommendation with poisson factorization. arXiv
preprintarXiv:1311.1704,2013.
Harper, F Maxwell and Konstan, Joseph A. The movie-
lens datasets: History and context. ACM Transactions
onInteractiveIntelligentSystems(TiiS),5(4):19,2015.
Hazan, Elad and Koren, Tomer. Linear regression with
limited observation. In Proceedings of the 29th Inter-
national Conference on Machine Learning (ICML-12),
pp.807–814,2012.
Hazan, Elad, Livni, Roi, and Mansour, Yishay. Classifi-
cation with low rank and missing data. arXiv preprint
arXiv:1501.03273,2015.
Hoffman, Matthew D, Blei, David M, Wang, Chong, and
Paisley,John.Stochasticvariationalinference.TheJour-
nal of Machine Learning Research, 14(1):1303–1347,
2013.
Hu,Yifan,Koren,Yehuda,andVolinsky,Chris. Collabora-
tivefilteringforimplicitfeedbackdatasets.InDataMin-
ing,2008.ICDM’08.EighthIEEEInternationalConfer-
enceon,pp.263–272.Ieee,2008.
CoupledCompoundPoissonFactorization
A.MixtureModel
A.1.GaussianMixtureModel
Generativemodelcanbedescribedas
• Foreachu=1,...,C andforeachcomponentk =1,...,K
U
1. Samples ∼N(η,ρ)
uk
• Foreachuandi
1. drawn fromthemissingness-encodingmodel
ui
2. sampley ∼N(φ (n )(cid:80) s ,φ (n )σ2).
ui 1 ui k uk 2 ui
LetthevariationaldistributionbeQ(s )=N(as ,bs ),thenstochasticvariationalupdatesare
uk uk uk
(cid:32)ησ2+E[φ /φ ]ρC y −E[φ2/φ ]ρC (cid:80) as (cid:33)
as ←(1−t−ξ)as +t−ξ 1 2 I ui 1 2 I j(cid:54)=k uj (1)
uk uk σ2+E[φ /φ ]ρC
1 2 I
(cid:18) ρσ2 (cid:19)
bs ←(1−t−ξ)bs +t−ξ . (2)
uk uk σ2+ρC E[φ /φ ]
I 1 2
A.2.PoissonMixtureModel
Generativemodelcanbedescribedas
• Foreachu=1,...,C andforeachcomponentk =1,...,K
U
1. Samples ∼Ga(η,ρ)
uk
• Foreachuandi
1. drawn fromthemissingness-encodingmodel
ui
(cid:80)
2. sampley ∼Po(φ(n ) s ).
ui ui k uk
LetthevariationaldistributionbeQ(s )=Ga(as ,bs ),thenstochasticvariationalupdatesare
uk uk uk
ϕ ∝exp{Ψ(as )−logbs } (3)
uk uk uk
as ←(1−t−ξ)as +t−ξ(η+C ϕ y ) (4)
uk uk I uk ui
bs ←(1−t−ξ)bs +t−ξ(ρ+C E[φ]). (5)
uk uk I