Bayesian Inference for Nonnegative Matrix Factorisation Models
Ali Taylan Cemgil
CUED/F-INFENG/TR.609
July 2008

University of Cambridge
Department of Engineering
Trumpington Street
Cambridge CB2 1PZ
United Kingdom
Email: atc27@cam.ac.uk
This research is funded by the Engineering and Physical Sciences Research Council (EPSRC) under the grant EP/D03261X/1
entitled “Probabilistic Modelling of Musical Audio for Machine Listening”).

Bayesian Inference for Nonnegative Matrix Factorisation Models∗Ali Taylan Cemgil
Signal Processing and Comms. Lab, University of Cambridge
Department of Engineering, Trumpington Street, Cambridge, CB2 1PZ, UK
atc27@eng.cam.ac.uk
July 2008
Abstract
We describe non-negative matrix factorisation (NMF) with a Kullback-Leibler error measure in a statistical framework, with a hierarchical generative model consisting of an observation and a prior component. Omitting the prior leads
to standard NMF algorithms as special cases, where maximum likelihood parameter estimation is carried out via the
Expectation-Maximisation (EM) algorithm. Starting from this view, we develop Bayesian extensions that facilitate more
powerful modelling and allow full Bayesian inference via variational Bayes or Monte Carlo. Our construction retains
conjugacy and enables us to develop models that fit better to real data while retaining attractive features of standard
NMF such as fast convergence and easy implementation. We illustrate our approach on model order selection and image
reconstruction.

1 Introduction
Nonnegative Matrix Factorisation (NMF) was introduced by Lee and Seung [16] as an alternative to k-means clusteringand principal component analysis (PCA) for data analysis and compression. In NMF, given a W × K nonnegative matrix
X = {xν,τ } where ν = 1 :W , τ = 1 :K, we seek positive matrices T and V such thatX
xν,τ ≈ [T V ]ν,τ =tν,i vi,τ
iwhere i = 1 :I. We will refer to the W ×I matrix T as the template matrix, and I ×K matrix V the excitation matrix. The
key property of NMF is that T and V are constrained to be positive matrices. This is in contrast with PCA, where there
are no positivity constraints or k-means clustering where each column of V is constrained to be a unit vector. Subject tothe positivity constraints, we seek a solution to the following minimisation problem(T, V )∗=

arg min D(X||T V )
T,V

(1)

Here, the function D is a suitably chosen error function. One particular choice for D, on which we will focus here, is theinformation (Kullback-Leibler) divergence, which we write as
¶
Xµ
λν,τ
D(X||Λ) = −(2)xν,τ log− λν,τ + xν,τxν,τν,τUsing Jensen’s inequality [2] and concavity of log x, it can be shown , that D(·) is nonnegative and D(X||Λ) = 0 if andonly if X = Λ. The objective in (1) could be minimised by any suitable optimisation algorithm. Lee and Seung [16] haveproposed a very efficient variational bound minimisation algorithm that has attractive convergence properties and which
has been successfully applied in various applications in signal analysis and source separation, e.g. [23, 13, 1].
The interpretation of NMF, like SVD (singular value decomposition), as a low rank matrix approximation is sufficientfor the derivation of a useful inference algorithm, yet this view arguably does not provide the complete picture about
assumptions underlying the statistical properties of X. Therefore, we describe NMF from a statistical perspective asa hierarchical model. In our framework, the original nonnegative multiplicative update equations of NMF appear as
an expectation-maximisation (EM) algorithm for maximum likelihood estimation of a conditionally Poisson model via
data augmentation. Starting from this view, we develop Bayesian extensions that facilitate more powerful modelling
and allow more sophisticated inference, such as Bayesian model selection. Inference in the resulting models can be
carried out easily using variational (structured mean field) or Markov Chain Monte Carlo (Gibbs sampler). The resulting
algorithms outperform existing NMF strategies and open up the way for a full Bayesian treatment for model selection
via computation of the marginal likelihoods (the evidence), such as estimating the dimensions of the template matrix or
regularising overcomplete representations via automatic relevance determination.
∗ This research is funded by the Engineering and Physical Sciences Research Council (EPSRC) under the grant EP/D03261X/1).

1

2 The statistical perspective
The interpretation of NMF as a low rank matrix approximation is sufficient for the derivation of an inference algorithm, yet
this view arguably does not provide the complete picture. In this section, we describe NMF from a statistical perspective.
This view will pave the way for developing extensions that facilitate more realistic and flexible modelling as well as moresophisticated inference, such as Bayesian model selection.
Our first step is the derivation of the information divergence error measure from a maximum likelihood principle. We
consider the following hierarchical model:
T ∼ p(T |Θt )V ∼ p(V |Θv )Xxν,τ =sν,i,τ

sν,i,τ ∼ PO(sν,i,τ ; tν,i vi,τ )

(3)
(4)

i

Here, PO(s; λ) denotes the Poisson distribution of the random variable s ∈ N0 with nonnegative intensity parameter λwherePO(s; λ) = exp(s log λ − λ − log Γ(s + 1))
and Γ(s + 1) = s! is the gamma function. The priors p(T |·) and p(V |·) will be specified later. We call the variables
Si = {sν,i,τ } latent sources. We can analytically marginalise out the latent sources S = {S1 . . . SI } to obtain themarginal likelihoodXYX
log p(X|T, V ) = log
p(X|S)p(S|T, V ) = log
PO(xν,τ ;tν,i , vi,τ )
ν,τ

S

=

XX
ν

i

(xν,τ log[T V ]ν,τ − [T V ]ν,τ − log Γ(xν,τ + 1))(5)τThis result follows from the well known superposition property of Poisson random variablesP[14], namely when si ∼PO(si ; λi ) and x = s1 +s2 +· · ·+sI then the marginal probability is given by p(x) = PO(x; i λi ). The maximisationof this objective in T and V is equivalent to the minimisation of the information divergence in (2). In the derivation oforiginal NMF in [17], this objective is stated first; the S variables are introduced implicitly later during the optimisation
on T and V . In the sequel, we show that this algorithm is actually equivalent to EM, ignoring the priors p(T |·) and p(V |·).

2.1 Maximum Likelihood and the EM algorithm
The loglikelihood of the observed data X can be written as
LX (T, V ) ≡log

Xp(X|S)p(S|T, V ) ≥S

X
Sq(S) logp(X, S|T, V )≡ BEM [q]
q(S)

(6)

where q(S) is an instrumental distribution, that is arbitrary provided that the sum on the right exists; q can only vanish at aparticular S only when p does so. Note that this defines
Pa lower bound to the loglikelihood. It can be shown via functional
derivatives and imposing the normalisation condition S q(S) = 1 via Lagrange multipliers that the lower bound is tightfor the exact posterior of the latent sources, i.e.
arg max BEM [q]

=

q(S)

p(S|X, T, V )Hence the loglikelihood can be maximised iteratively
E Step
M Stepq(S)(n) = p(S|X, T (n−1) , V (n−1) )(T (n) , V (n) ) = arg max hlog p(S, X|T, V )iq(S)(n)T,V

R

Here, hf (x)ip(x) = p(x)f (x)dx, the expectation of some function f (x) with respect to p(x).In the E step, we computethe posterior distribution of S. This defines a lower bound on the likelihood
B (n) (T, V |T (n−1) , V (n−1) )= hlog p(S, X|T, V )iq(S)(n)For many models in the exponential family, which includes (4), the expectation on the right depends on the sufficient
statistics of q(S)(n) and is readily available; in fact “calculating q(S)” should be literally taken as calculating the sufficient
statistics of q(S). The lower bound is readily obtained as a function of these sufficient statistics and maximisation in the
M Step yields a fixed point equation.

22.1.1 The E Step
To derive the posterior of the latent sources, we observe thatp(S|X, T, V )For the model in (4), we havelog p(S, X|T, V )== p(S, X|T, V )/p(X|T, V )

Ã
XX X
ν

τ

(7)

(−tν,i vi,τ + sν,i,τ log(tν,i vi,τ )i− log Γ(sν,i,τ + 1)) + log δ(xν,τ −X!
sν,i,τ )(8)i

It follows from (4), (7), (8) and (5)log p(S|X, T, V )

Ã
Ã
XX X

=

ν

τ

sν,i,τ log(tν,i vi,τ /X!
tν,i0 vi0 ,τ ) − log Γ(sν,i,τ + 1)i0i

+ log Γ(xν,τ + 1) + log δ(xν,τ −XX=
where pν,i,τ ≡ tν,i vi,τ /

ν

P

X

!
sν,i,τ )ilog M(sν,1,τ , . . . , sν,I,τ ; xν,τ , pν,1,τ , . . . , pν,I,τ )(9)τ

i0 tν,i vi ,τ are the cell probabilities. Here, M denotes a multinomial distribution defined by00µM(s; x, p)

=

¶
I
X
X
Y
psi i
x
si ) = δ(x −si )x!
ps11 ps22 · · · psII δ(x −s1 s2 . . . sI
s!
i
i
i=1 i

where s = {s1 , s2 , . . . , sI } and p = {p1 , p2 , . . . , pI } and p1 + p2 + · · · + pI = 1. Here, pi , i = 1 . . . I are the cell
probabilities and x is the index parameter where s1 + s2 + · · · + sI = x. The Kronecker delta function is defined by
δ(x) = 1 when x = 0, and δ(x) = 0 otherwise. It is a standard result that the marginal mean is
hsi i

= xpi

i.e., the expected value each source si is a fraction of the observation, where the fraction is given by the corresponding
cell probability.
2.1.2 The M Step
It is indeed a good news that the posterior has an analytic form since now the M step can be calculated easilyÃXX X
hlog p(S, X|T, V )ip(S|X,T,V ) =(−tν,i vi,τ + hsν,i,τ i log(tν,i vi,τ )
ν

τ

i

*

− hlog Γ(sν,i,τ + 1)i) +log δ(xν,τ −X+!
sν,i,τ )i

Fortunately, for maximisation w.r.t. T and V , the last two difficult terms are merely constant and we need only to maximise
the simpler objective
Ã
!
´
XX X³
(n)
−tν,i vi,τ + hsν,i,τ i log(tν,i vi,τ )Q(T, V ) =ντiwhere we only need the expected value of the sources given by the previous values of the templates and excitations(n) (n)(n)hsν,i,τ i=tν,i vi,τ
xν,τ P (n) (n)i0 tν,i0 vi0 ,τ(n)Maximisation of the objective Q and substituting hsν,i,τ i gives the following fixed point equations:
X (n) X∂Q
(n)vi,τ +
hsν,i,τ i /tν,i
= −
∂tν,i
τ
τ
P
(n) P
(n) (n)
X (n)
X
(n+1)
(n)
τ xν,τ vi,τ /i0 tν,i0 vi0 ,τ
(n)vi,τ = tν,itν,ihsν,i,τ i /
=
P (n)
τ
τ
τ vi,τ
3(10)

V
i = 1...I

Θv

vi,1

Θt

tν,i

···

vi,τ

···

vi,K

···

sν,i,τ

···sν,i,Kvi,τ
T
sν,i,τ

tν,i

S

xν,τ
sν,i,1

i
ν

xν,1

X

xν,τ

xν,K
ν = 1...WτFigure 1: (Left) A schematic description of the NMF model with data augmentation. (Right) Graphical model withhyperparameters.
Each source element sν,i,τ is PoissonP
P distributed with intensity tν,i vi,τ . The observations are given byxν,τ = i sν,i,τ . In matrix notation, we write X =Si . We can analytically integrate out over S. Due to superposition
property of Poisson distribution, intensities add up and we obtain hXi = T V . Given X, the NMF algorithm is shownto seek the maximum likelihood estimates of the templates T and excitations V . In our Bayesian treatment, we further
assume elements of T and V are Gamma distributed with hyperparameters Θ.

∂Q
∂vi,τ
(n+1)
vi,τ

=

−

X (n) X(n)tν,i +
hsν,i,τ i /vi,τ
ν

=

X

ν

hsν,i,τ i

(n)

X (n)
(n)
/
tν,i = vi,τνν

P (n) (n)(n)
i0 tν,i0 vi0 ,τν tν,i xν,τ /P (n)ν tν,iP(11)

The equations (10) and (11) are identical to the multiplicative update rules of [17]. However, our derivation via data
augmentation obtains the same result as an EM algorithm. It is interesting to note that in the literature, NMF is often
described as EM-like; here, we show that it is actually just an EM algorithm. We see that the efficiency of NMF is dueto the fact that the W × I × K object hSi need not be explicitly calculated as we only need its marginal statistics (sumsacross τ or ν).
We note that our model is valid when X is integer valued. See [15] for a detailed discussion about consequences ofthis issue. Here, we assume that for nonnegative real valued X̃ we only consider the integer part, i.e. we let X̃ = X + E
where E is a noise matrix with entries uniformly drawn in [0, 1). In practice, this is not an obstacle when the entries of X
are large.
The interpretation of NMF as a maximum likelihood method in a Poisson model is mentioned in the original NMF
paper [16] and discussed in more detail by [13, 22]. Kameoka in [13] focuses on the optimisation and gives an equivalentdescription using auxiliary function maximisation. In contrast, the auxiliary variables can be viewed as model variables(the sources s) that are analytically integrated out [22]. However, none of these approaches provided a full Bayesiantreatment.

2.2 Hierarchical Prior Structure
Given the probabilistic interpretation, it is possible to propose various hierarchical prior structures to fit the requirementsof an application. Here we will describe a simple choice where we have a conjugate priortν,i∼G(tν,i ; atν,i , btν,i /atν,i )vi,τ ∼ G(vi,τ ; avi,τ , bvi,τ /avi,τ )(12)Here, G denotes the density of a gamma random variable x ∈ R+ with shape a ∈ R+ and scale b ∈ R+ defined byG(x; a, b)=exp((a − 1) log x − x/b − log Γ(a) − a log b)

The primary motivation for choosing a Gamma distribution is computational convenience: Gamma distribution is the
conjugate prior to Poisson intensity. The indexing highlights the most general case where there are individual parameters
for each element tν,i and vi,τ . Typically, we don’t allow many free hyperparameters but tie them depending upon therequirements of an application. See figure 1 for an example. As an example, consider a model where we tie the hyperparameters such as atν,i = at , btν,i = bt , avi,τ = av and bvi,τ = bv for i = 1 . . . I, ν = 1 . . . W and τ = 1 . . . K. This model issimple to interpret, where each component of the templates and the excitations is drawn independently from the Gammafamily shown in Figure 2. Qualitatively, the shape parameter√ a controls the sparsity of the representation. Rememberthat G(x; a, b/a) has the mean b and standard deviation b/ a. Hence, for large a, all coefficients will have more or lessthe same magnitude b and typical representations will be full. In contrast, for small a, most of the coefficients will be4a = 0.05

Gamma(x; a, 1/a)
5
a=100

4

a=1

3
2
a = 50a=10
1
0

a=0.01
0

vi

a=1
a=0.1

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

i

2

Figure 2: (Left) The family of densities p(v; a, b = 1) = G(v; a, b/a) with the same mean hvi = b = 1. (Right) Typicalindependent draws vi ∼ p(v; a, 1) for a = [0.05, 1, 50]. Small values of a enforce (sparser representations) and large
values tie all values to be close to a nonzero mean (nonsparse representation).

very close to zero and only very few will be dominating, hence favoring a sparse representation. The scale parameter b isadapted to give the expected magnitude of each component.
To model missing data, i.e., when some of the xν,τ are not observed, we define a mask matrix M = {mν,τ }, the same
size as X where, mν,τ = 0, if xν,τ is missing and 1 otherwise (see appendix A.4 for details). Using the mask variables,
the observation model with missing data can be written asYmp(X|S)p(S|T, V ) =(p(xν,τ |sν,1:I,τ )p(sν,1:I,τ |tν,1:I , v1:I,τ )) ν,τν,τ

The hierarchical model in (12) is potentially more flexible than the basic model of (4), in that it allows a lot of freedom
for more realistic modelling. First of all, the hyperparameters can be estimated from examples of a certain class of source
to capture the invariant features. Another possibility is Bayesian model selection, where we can compare alternativemodels in terms of their marginal likelihood. This enables one to estimate the model order, for example, the optimum
number of templates to represent a source.

3 Full Bayesian Inference
Below, we describe various interesting problems that can be cast to Bayesian inference problems. In signal analysis
and feature extraction with NMF we may wish to calculate the posterior of templates and excitations given data andhyperparameters Θ ≡ (Θt , Θv ). Another important quantity is the marginal likelihood (also known as the evidence)
where
Z
X
p(X|Θ) =dT dVp(X|S)p(S|T, V )p(T, V |Θ)S

The marginal likelihood can be used to estimate the hyperparameters, given examples of a source classΘ∗= arg max p(X|Θ)Θ

or to compare two given models via Bayes factorsl(Θ1 , Θ2 )=p(X|Θ1 )p(X|Θ2 )This latter quantity is particularly useful for comparing different classes of models. Unfortunately, the integrations
required can not be computed in closed form. In the sequel, we will describe the Gibbs sampler and variational Bayes asapproximate inference strategies.

3.1 Variational Bayes
We sketch here the Variational Bayes (VB) [9, 2] method to bound the marginal loglikelihood as
XZp(X, S, T, V |Θ)d(T, V )q log
LX (Θ) ≡ log p(X|Θ) ≥= hlog p(X, S, V, T |Θ)iq + H[q] ≡ BV B [q]
q
S

5

(13)

where, q = q(S, T, V ) is an instrumental distribution and H[q] is its entropy. The bound is tight for the exact posteriorq(S, T, V ) = p(S, T, V |X, Θ), but as this distribution is complex we assume a factorised form for the instrumental
distribution by ignoring some of the couplings present in the exact posterior


Ã
!
Y
Y
Y
Y
q(S, T, V ) = q(S)q(T )q(V ) =q(sν,1:I,τ )  q(tν,i )  q(vi,τ ) ≡
qα
ν,τ

ν,i

i,τ

α∈C

where α ∈ C = {{S}, {T }, {V }} denotes set of disjoint clusters. Hence, we are no longer guaranteed to attain theexact marginal likelihood LX (Θ). Yet, the bound property is preserved and the strategy of VB is to optimise the bound.
Although the best q distribution respecting the factorisation is not available in closed form, it turns out that a local optimum
can be attained by the following fixed point iteration:
³
´
(n+1)
(14)
qα
∝ exp hlog p(X, S, T, V |Θ)iq(n)¬αwhere q¬α = q/qα . This iteration monotonically improves the individual factors of the q distribution, i.e. B[q (n) ] ≤B[q (n+1) ] for n = 1, 2, . . . given an initialisation q (0) . The order is not important for convergence, one could visit blocks
in arbitrary order. However, in general, the attained fixed point depends upon the order of the updates as well as thestarting point q (0) (·). We choose the following update order in our derivations,³´q(S)(n+1) ∝ exp hlog p(X, S, T, V |Θ)iq(T )(n) q(V )(n)(15)³´
q(T )(n+1) ∝ exp hlog p(X, S, T, V |Θ)iq(S)(n+1) q(V )(n)(16)³´
q(V )(n+1) ∝ exp hlog p(X, S, T, V |Θ)iq(S)(n+1) q(T )(n+1)(17)

3.2 Variational update equations and sufficient statistics
The expectations of hlog p(X, S, T, V |Θ)i are functions of the sufficient statistics of q (see the expression in the appendix
A.2). The update equation for the latent sources (15) leads to the following:Ã!
XX
q(sν,1:I,τ ) ∝ exp(sν,i,τ (hlog tν,i i + hlog vi,τ i) − log Γ(sν,i,τ + 1)) δ(xν,τ −sν,i,τ )ipν,i,τi∝ M(sν,1,τ , . . . , sν,i,τ , . . . , sν,I,τ ; xν,τ , pν,1,τ , . . . , pν,i,τ , . . . , pν,I,τ )X= exp(hlog tν,i i + hlog vi,τ i)/exp(hlog tν,i i + hlog vi,τ i)

(18)

i

hsν,i,τ i= xν,τ pν,i,τ

(19)

These equations are analogous to the multinomial posterior of EM given in 9; only the computation of cell probabilities
is different. The excitation and template distributions and their sufficient statistics follow from properties of the gamma
distribution.
Ã
Ã
!
!
t
X
X
a
ν,i
q(tν,i ) ∝ exp (atν,i +
hsν,i,τ i − 1) log(tν,i ) −+hvi,τ i tν,i
bν,i
τ
τ
¡
¢
t
t
∝ G tν,i ; αν,i, βν,itαν,i≡ atν,i +

X

Ã
t
βν,i
≡

hsν,i,τ iτtt
exp(hlog tν,i i) = exp(Ψ(αν,i))βν,iÃ
q(vi,τ )∝ exp (avi,τ +
∝

X

t
t
htν,i i = αν,iβν,iÃhsν,i,τ i − 1) log vi,τ −

ν
¡
¢
v
v
G vi,τ ; αi,τ , βi,τv
αi,τ≡ avi,τ +X!−1atν,i Xhvi,τ i
+
btν,i
τ

!
!
avi,τ Xhtν,i i vi,τ
+
bvi,τ
ν

Ã
v
βi,τ
≡

hsν,i,τ iνvv
exp(hlog vi,τ i) = exp(Ψ(αi,τ
))βi,τ

!−1avi,τ Xhtν,i i
+
bvi,τ
ν

v
v
hvi,τ i = αi,τβi,τ6

3.3 Efficient implementation
One of the attractive features of NMF is easy and efficient implementation. In this section, we derive that the update
equations of section 3.2 in compact matrix notation to illustrate that these attractive properties are retained for the full
Bayesian treatment. A subtle but key point in the efficiency of the algorithm is that we can avoid explicitly storing andcomputing the W × I × K object hSi, as we only need the marginal statistics during optimisation. Consider Eq. 18 and
19. We can writeXX
hsν,i,τ i =xν,τ pν,i,τττ

=

exp(hlog tν,i i)XÃxν,τ /

Ã
X

τ

Σt

!!
exp(hlog tν,i0 i) exp(hlog vi0 ,τ i)exp(hlog vi,τ i)i0

Lt .∗((X ./(Lt Lv ))L>v)=

Here, the denominator has to be nonzero. In the last line, we have represented the expression in compact notation where
we define the following matrices:Xtt
Et = {htν,i i}, Lt = {exp(hlog tν,i i)}, Σt = {hsν,i,τ i}, At = {atν,i }, Bt = {btν,i }, α t = {αν,i}, β t = {βν,i}τXEv = {hvi,τ i}, Lv = {exp(hlog vi,τ i)}, Σv = {vv
hsν,i,τ i}, Av = {avi,τ }, Bv = {bvi,τ }, αv = {αi,τ}, β v = {βi,τ}ν
×IThe matrices subscripted with t are in RW
and with v are in RI×K
. For notational convenience, we define .∗ and ./++as element wise matrix multiplication and division respectively and 1W as a W × 1 vector of ones. After straightforward
substitutions, we obtain the variational nonnegative matrix factorisation algorithm, that can compactly be expressed as inpanel Algorithm 1).

Algorithm 1 Variational Nonnegative Matrix Factorisation
1: Initialise :
(0)

Lt

(0)

=

Et

∼ G(·; At , Bt ./ At )

(0)
L(0)
∼ G(·; Av , Bv ./ Av )v = Ev(n−1)

2: for n = 1 . . . MAXITER do
3:
Source sufficient statistics

4:

5:
6:

Σt

(n)

:=

Lt

(n−1)

Σ(n)
v

:=

L(n−1).∗(Ltv.∗(((X .∗ M ) ./(Lt(n−1) >>

L(n−1)))Lv(n−1) )v(n−1)((X .∗ M )./(LtL(n−1)
)))
v

Means
Et

(n)

:=

(n)
(n)
α t .∗ β t(n)(n)
α t = At + Σ t

Ev(n)

:=

(n)
α (n)
v .∗ β v

α (n)= Av + Σ(n)
v
v

´
³
>
(n)
β t = 1./ At ./Bt + M Ev(n−1)
³
´
(n) >
β (n)
= 1./ Av ./Bv + Et
Mv

Optional: Compute Bound (See appendix, (28))
Means of Logs
(n)

Lt

=

(n)
α(n)
exp(Ψ(αt )) .∗ β t

(n)
α(n)
L(n)
= exp(Ψ(αvv )) .∗ β v

7:Optional: Update Hyperparameters (See appendix, section A.5)
8: end for

Similarly, an iterative conditional modes (ICM) algorithm can be derived to compute the maximum a-posteriori (MAP)
solution (see section A.4)
V
T:= (Av + V .∗(T > ((M .∗ X) ./(T V )))) ./(Av ./ Bv + T > M ):= (At + T .∗(((M .∗ X) ./(T V ))V > )) ./(At ./ Bt + M V > )

(20)
(21)

Note that when the shape parameters go to zero, i.e. At , Av → 0, we obtain the maximum likelihood NMF algorithm.

3.4 Markov Chain Monte Carlo, the Gibbs sampler
Monte Carlo methods [10, 18] are powerful computational techniques to estimate expectations of formN

E = hf (x)ip(x) ≈1 X
f (x(i) ) = ẼN
N n=1
7

(22)

where x(i) are independent samples drawn from p(x). Under mild conditions on f , the estimate ẼN converges to thetrue expectation for N → ∞. The difficulty here is obtaining independent samples {x(i) }i=1...N from complicateddistributions.
The Markov Chain Monte Carlo (MCMC) techniques generate subsequent samples from a Markov chain defined by atransition kernel T , that is one generates x(i+1) conditioned on x(i)x(i+1)

T (x|x(i) )∼Note that the transition kernel T is not needed explicitly in practice; all is needed is a procedure to sample a new configuration given the previous one. Perhaps surprisingly, even though subsequent samples are correlated, provided that
T satisfies certain ergodicity conditions, (22) remains still valid and estimated expectations converge to their true values
when number of samples N , goes to
R infinity [10]. To design a transition kernel T such that the desired distribution is thestationary distribution, i.e. p(x) = dx0 T (x|x0 )p(x0 ), many alternative strategies can be employed [18]. One particularly
convenient and simple procedure is the Gibbs sampler where one samples each block of variables from full conditional
distributions. For the NMF model, a possible Gibbs sampler is
S (n+1)T (n+1)V (n+1)

p(S|T (n) , V (n) , X, Θ)

∼

(n)

∼
∼

(n+1)

p(T |V , S, X, Θ)p(V |S (n+1) , T (n+1) , X, Θ)

(23)
(24)
(25)

Note that this procedure implicitly defines a transition kernel T (·|·). It can be shown [10] that the stationary distribution of
T is the exact posterior p(S, T, V |X, Θ). Eventually, the Gibbs sampler converges regardless of the order that the blocks
are visited, provided that each block is visited infinitely often in the limit n → ∞. However the rate of convergence isvery difficult to assess as it depends upon the order of the updates as well as the starting configuration (T (0) , V (0) , S (0) ).
It is instructive to contrast above equations (23)-(25) with the variational update equations(15)-(17) : algorithmically thetwo approaches are quite similar. The pseudo-code is given in Algorithm 2.
3.4.1 Marginal Likelihood estimation with Chib’s method
The marginal likelihood can be estimated from the samples generated by the Gibbs sampler using a method proposed by
Chib [6]. Suppose we have run the block Gibbs sampler until convergence and have N samples{T (n) }n=1:N{V (n) }n=1:N{S (n) }n=1:N

The marginal likelihood is (omitting hyperparameters Θ)p(X)=p(V, T, S, X)
p(V, T, S|X)

(26)

This equation holds for all points (V, T, S). We choose a point in the configuration space; provided that the distribution
is unimodal, a good candidate is a configuration near the mode (T̃ , Ṽ , S̃). The numerator in (26) is easy to evaluate. The
denominator isp(Ṽ , T̃ , S̃|X) = p(Ṽ |T̃ , S̃, X)p(T̃ |S̃, X)p(S̃|X)= p(Ṽ |T̃ , S̃)p(T̃ |S̃)p(S̃|X)

The first term is the full conditional so it is available for the Gibbs sampler. The third term isZ
p(S̃|X) =N

dV dT p(S̃|V, T, X)p(V, T |X) ≈The second term is trickier

1 X
p(S̃|V (n) , T (n) , X)
N n=1

(27)

Z
p(T̃ |S̃) =dV p(T̃ |V, S̃)p(V |S̃)The first term here is the full conditional. However, the original Gibbs run gives us only samples from p(V |X), not(m)(m)
p(V |S̃). The idea is to run the Gibbs sampler for M further iterations where we sample from (VS̃ , TS̃ ) ∼ p(V, T |S =S̃), i.e. with S clamped at S̃. The resulting estimate isMp(T̃ |S̃)

≈

1 X
(m)
p(T̃ |VS̃ , S̃)
M m=18

Chib’s method estimates the marginal likelihood as follows:
log p(X|Θ)= log p(Ṽ , T̃ , S̃, X|Θ) − log p(Ṽ , T̃ , S̃|X, Θ)≈log p(Ṽ , T̃ , S̃, X|Θ) − log p(Ṽ |T̃ , S̃, Θ)
− log

M
X
m=1

(m)

p(T̃ |VS̃, S̃, Θ) − log

NX

p(S̃|V (n) , T (n) , X, Θ) + log(M N )n=1

Algorithm 2 Gibbs sampler for Nonnegative Matrix Factorisation
1: Initialize :
T (0)

=

V (0) ∼ G(·; Av , Bv )∼ G(·; At , Bt )2: for n = 1 . . . MAXITER do
3:
Sample Sources
4:
for τ = 1 . . . K, ν = 1 . . . W do(n)5:pν,1 :I,τ = T (n−1) (ν, 1 :I) .∗ V (n−1) (1 :I, τ )> ./(T (n−1) (ν, 1 :I)V (n−1) (1 :I, τ ))
6:
7:

(n)

S (n) (ν, 1 :I, τ ) ∼ M(sν,1:I,τ ; xν,τ , pν,1:I,τ )end for
X (n)
(n)
Σt
=
Sν,i,τ

=
Σ(n)
v

τ

8:

9:

X

(n)

Sν,i,τ

νSample Templates
³
´
β (n)
= 1./ At ./Bt + 1W (V (n−1) 1K )>
t

α(n)
t

=

At + Σt

(n)

T (n)

∼

G(T ; α t , β t )

α(n)
v

=

Av + Σ(n)
v

V (n)

∼

(n)
G(V ; α (n)v ,βv )

(n)(n)

Sample Excitations³´(n−1) > >β (n)= 1./ Av ./Bv + (1>
) 1K
v
WT

10: end for

4 Simulations
Our goal is to illustrate our approach in a model selection context. We first illustrate that the variational approximation
to the marginal likelihood is close to the one obtained from the Gibbs sampler via Chib’s method. Then, we compare thequality of solutions we obtain via Variational NMF and compare them to the original NMF on a prediction task. Finally,
we focus on reconstruction quality in the overcomplete case where the standard NMF is subject to overfitting.
Model Order Determination: To test our approach, we generate synthetic data from the hierarchical model in(12) with W = 16, K = 10 and the number of sources Itrue = 5. The inference task is to find the correct number ofsources, given X. The hyperparameters of the true model are set to atν,i = at = 10, btν,i = bt = 1, avi,τ = av = 1,bvi,τ = bv = 100. In the first experiment the hyperparameters are assumed to be known and in the second are jointly
estimated from data, using hyperparameter adaptation. We evaluate the marginal likelihood for models with the numberof templates I = 1 . . . 10, with the Gibbs sampler using Chib’s method and variational lower bound B via variational
Bayes. We run the Gibbs sampler for MAXITER = 10000 steps following a burn-in period of 5000 steps; then we clamp
the sources S and continue the simulation for a further 10000 steps to estimate quantities required by Chib’s method. We
run the variational algorithm until convergence of the bound or 10000 iterations, whichever occurs first. In Figure 3-top,
we show a comparison of the variational estimate with the average of 5 independent runs obtained via Chib’s method.
We observe, that both methods give consistent results. In Figure 4, we show the lower bound as a function of modelorder I, where for each I, the bound is optimised independently by jointly optimising hyperparameters at , bt , av and bv
using the equations derived in the Appendix. We observe, that the correct model order can be inferred even when thehyperparameters are unknown a-priori. This is potentially useful for estimation of model order from real data.
As real data, we use a version of the Olivetti face image database ( K = 400 images of 64 × 64 pixels available at
http://www.cs.toronto.edu/˜roweis/data/olivettifaces.mat). We further downsampled to 16×16
or 32 × 32 pixels, hence our data matrix X is 162 × 400 or 322 × 400. We use a model with tied hyperparameters as
atν,i = at , btν,i = bt , avi,τ = av and bvi,τ = bv , where all hyperparameters are jointly estimated. In Figure 4, bottom,
we show results of model order determination for this dataset with joint hyperparameter adaptation. Here, we run thevariational algorithm for each model order I = 1 . . . 100 independently and evaluate the lower bound after optimising the9hyperparameters. The Gibbs sampler is not found practical and is omitted here. The lower bound behaves as is expectedfrom marginal likelihood, reflecting the tradeoff between too many and too few templates. Higher resolution implies more
templates, consistent with our intuition that detail requires more templates for accurate representation.

Log Evidence (Lower Bound)

−1500
−2000
−2500
−3000
Variational Bound
Chib’s

−3500
Itrue
−4000

1

2

3

4

5
6
I (Model order)

7

8

9

10

Log Evidence (Lower Bound)

−1500
−2000
−2500
−3000
−3500
−4000

1

2

3

4

5
6
I (Model order)

7

8

9

10

Log Evidence (Lower Bound)

−1000

−1500

−2000

−2500

1

2

3

4

5
6
I (Model order)

7

8

9

10

Figure 3: Model selection results. (top) Comparison of model selection by variational bound (squares) and marginal
likelihood estimated by Chib’s (circles) method. The hyperparameters are assumed to be known. (Middle) Box-plot ofmarginal likelihood estimates by Chib’s method using 5000, 10000 and 10000 iterations for burn-in, free and clamped
sampling. The boxes show the lower quartile, median, and upper quartile values. (Bottom) Model selection by variationalbound when hyperparameters are unknown and jointly estimated.

6

5

−4.5

−1.8

x 10

x 10

B (Lower Bound)

B (Lower Bound)

−1.9

−5

−5.5

−2
−2.1
−2.2

−6
−2.3

−6.5

10

20

30

40

50
60
I (Model Order)

70

80

90

−2.4

100

10

20

30

40

50
60
I (Model Order)

70

80

90

100

Figure 4: Model selection using variational bound with adapted hyperparameters on face data 16 × 16 with I ∗ = 27 (left)and 32 × 32 with I ∗ = 42 (right).
We also investigate the nature of the representations (see Figure 5). Here, for each independent run, we fix the valuesof shape parameters to (at , av ) = [(10, 10), (0.1, 0.1), (10, 0.2), (10, 0.5)] and only estimate bt and bv . This corresponds
to enforcing sparse or non-sparse t and v. Each columns shows I = 36 templates estimated from the dataset conditionedon hyperparameters. The middle image is the same template image above weighted with the excitations corresponding
to the reconstruction (the expected value of the predictive distribution) below. Here, we clearly see the effect of thehyperparameters. In the first condition (at , av ) = (10, 10), the prior does not enforces sparsity to the templates andexcitations. Hence, for the representation of a given image, there are many active templates. In the second condition, we
try to force both matrices to be sparse with (at , av ) = (0.1, 0.1) . Here, the result is not satisfactory as isolated componentsof the templates are zeroed, giving a representation that looks like one contaminated by “salt-and-pepper” noise. The third
condition ((at , av ) = (10, 0.2)) forces only the excitations to be sparse. Here, we observe that the templates correspond to10

at = 10 av = 10

at = 0.1 av = 0.1

at = 10 av = 0.2

at = 0.5 av = 10

B = −2051618.1672

B = −2519222.8848

B = −2081517.199

B = −1965293.3874

Figure 5: Templates, excitations for a particular example and the reconstructions obtained for different hyperparametersettings. B is the lower bound for the whole dataset.

some average face images. Qualitatively, each image is reconstructed using a superposition of a few of these templates. In
the final representation, we enforce sparsity in the templates but not in the excitations. Here, our estimate finds templates
that correspond to parts of individual face images (eyebrows, lips etc.). This solution, intuitively corresponding to aparsimonious representation, also is the best in terms of the marginal likelihood. With proper initialisation, our variational
procedure is able to find such solutions.
Prediction: We now compare variational Bayesian NMF with the maximum likelihood NMF on a missing data
prediction task. To illustrate the self regularisation effect, we set up an experiment in which we select a subset of the facedata consisting of 50 images. From half of the images, we remove a the same patch (Fig. 6) and try to predict missing
pixels. This is a rather small dataset for this task, as we have only 10 images for each of the 5 different persons, and half ofthese images have missing data at the same spot. We measure the quality of the prediction in terms of signal-to-noise ratio(SNR). The missing values are reconstructed using the mean of the predictive distribution Xpred ≡ hXiPO(X;T ∗ V ∗ ) =T ∗ V ∗ where T ∗ and V ∗ are point estimates of the template and excitation matrix. We compare our variational algorithmwith the classical NMF. For each algorithm, we test two different versions. The variational algorithms differ in how
we estimate T ∗ and V ∗ . In the first variational algorithm, we use a crude estimate of T ∗ and V ∗ as the mean of theapproximating q distribution. In the second condition, after convergence of hyperparameters via VB, we reinitialise T and
V randomly and switch to an ICM algorithm (see Eq.21). This strategy finds a local mode (T ∗ , V ∗ ) of the exact posterior
distribution. In NMF, we test two initialisation strategies: in the first condition, we initialise the templates randomly. In
the second we set them equal to the images in the dataset with random perturbations.
In Fig. 6, we show the reconstruction results for a typical run, for a model with I = 100 templates. Note that thisan overcomplete model, with twice as many templates as there are images. To characterise the nature of the estimatedtemplateand
we use thePexcitation matrices,P√√ sparseness criteria [11] of an m × n matrix X, defined as Sparseness(X) =2 1/2( mn − ( i,j |Xi,j |)/( i,j Xi,j) )/( mn − 1). This measure is 1 when the matrix X has only a single non-zero
entry and 0 when all entries are equal. We see that the variational algorithms are superior in this case in terms of SNR aswell as the visual quality of the reconstruction. This is perhaps not surprising, since with maximum likelihood estimation;
if the model order is not carefully chosen, generalisation performance is poor: the “noise” in the observed data is fitted
but the prediction quality drops on new data. An interesting observation is that highly sparse solutions (either in templates
or excitations) do not give the best result, the solution that balances both seems to be the best in this setting. This example
illustrates that sparseness in itself may not be necessarily a good criteria to optimise; for model selection, the marginal
likelihood should be used as the natural quantity.

11

Sparseness (Templates)

1

0.6

Data (25% Missing)

Variational Bayes NMF

ML NMF1
7.1 dB
VB+ICM
13.3 dB

0.4

ML NMF2
0.2
0Ground truth

VB
12.6 dB

0.8

10.4 dB

00.5
Sparseness (Excitations)

Variational Bayes+ICM NMF

ML NMF1

1
ML NMF2

Figure 6: Results of a typical run. (Top) Example images from the data set. Left to right. The ground truth, data withmissing pixels. The reconstructions of VB, VB+ICM, and ML NMF with two initialisation strategies (1=random, 2=to
image). (Bottom) Comparison of the reconstruction accuracy of different methods in terms of SNR (in dB), organised
according to the sparseness of the solution.

4.1 Discussion and Conclusions
In this paper, we have investigated NMF from a statistical perspective. We have shown that KL minimisation formulation
the original algorithm can be derived from a probabilistic model where the observations are superposition of I independent
Poisson distributed latent sources. Here, the template and excitation matrices turn out to be latent intensity parameters.
The interpretation of NMF as a maximum likelihood method in a Poisson model is mentioned in the original NMF paper[16] and discussed in more detail by [13, 22]. [13] focuses on the optimisation and gives an equivalent description
using auxiliary function maximisation. In contrast, [22] illustrates that the auxiliary variables can be viewed as model
variables (the sources s) that are analytically integrated out. The novel observation in the current article is the exact
characterisation of the approximating distribution q(S) or full conditionals p(S|T, V, X) as a product of multinomial
distributions, leading to a richer approximation distribution than a naive mean field or single site Gibbs (which would
freeze due to deterministic p(X|S)). This conjugate form leads to significant simplifications in full Bayesian integration.
Apart from the conditionally Gaussian case, NMF with KL objective seems to be unique in this respect. For several otherdistance metrics D(.||.), we find full Bayesian inference not as practical, as p(S|T, V, X) is not standard.
We have also shown that the standard NMF algorithm with multiplicative update rules is in fact an EM algorithm withdata augmentation. Extending upon this observation, we have developed an hierarchical model with conjugate Gamma
priors. We have developed a variational Bayes algorithm and a Gibbs sampler for inference in this hierarchical model. We
have also developed methods for estimating the marginal likelihood for model selection.
Our simulations suggest that the variational bound seems to be a reasonable approximation to the marginal likelihoodand can guide model selection for NMF. The computational requirements are comparable to the ML NMF. A potentiallytime consuming step in the implementation of the variational algorithm is the evaluation of the Ψ function, but this step
can also be replaced by a simple piecewise polynomial approximation since exp(Ψ(x)) ≈ x − 0.5 for x > 5.
We first compare the variational inference with a Gibbs sampler. In our simulations, we observe that both algorithms
give qualitatively very similar results, both for inference of templates and excitations as well as model order selection.
We find the variational approach somewhat more practical as it can be expressed as simple matrix operations, where boththe fixed point equations as well as the bound can be compactly and efficiently implemented using matrix computation
software. In contrast, our Gibbs sampler is computationally more demanding and the calculation of marginal likelihood
is somewhat more tricky. With our implementation of both algorithms the variational method is faster by a factor ofaround 13. Reference implementations of both algorithms in matlab are available from the following url: http://
www-sigproc.eng.cam.ac.uk/˜atc27/bnmf/.
In terms of computational requirements, the variational procedure has several advantages. First, we circumvent sampling from multinomial variables, which is the main computational bottleneck with the Gibbs sampler. Whilst efficientalgorithms are developed for multinomial sampling [7], the procedure is time consuming when the number of latentsources I is large. In contrast, the variational method estimates the expected sufficient statistics directly by elementary12matrix operations. Another advantage is hyperparameter estimation. In principle, it is possible to maximise the marginal
likelihood via a Monte Carlo EM procedure [21, 19], yet this potentially requires many more iterations. In contrast,
the evaluation of the derivatives of the lower bound is straightforward and can be implemented without much additional
computational cost.
The efficiency of the Gibbs sampler
improved by working out the distribution of the sufficient statistics of
P could beP
sources directly (namely quantities τ sν,i,τ or ν sν,i,τ ) to circumvent multinomial sampling. Unfortunately, for thesum of binomial random variables with different cell probability parameters, the sum does not have a simple form but
various approximations are possible [4].
From a modelling perspective, our hierarchical model provides some attractive properties. It is easy to incorporate
prior knowledge about individual latent sources via hyper parameters and one can easily capture variability in the templates
and excitations that is potentially useful for developing robust techniques. The prior structure here is qualitatively similar
to an entropic prior [3, 20] and we find qualitatively similar representations to the ones found by NMF reported earlier
by [16, 11]. However, none of the above mentioned methods provide an estimate of the marginal likelihood, which isuseful for model selection. Our generative model formulation can be extended in various ways to suit the specific needsof particular applications. For example, one can enforce more structured prior models such as chains or fields [22]. As asecond possibility, the Poisson observation model can be replaced with other models such as clipped Gaussian, Gamma
or Gaussians which lead to alternative source separation algorithms. For example, the case of Gaussian sources where theexcitations and templates correspond to the variances is discussed in [5].
Our main contribution here is the development of a principled and practical way to estimate both the optimal sparsity
criteria and model order, in terms of marginal likelihood. By maximising the bound on marginal likelihood, we havea method where all the hyperparameters can be estimated from data, and the appropriate sparseness criteria is foundautomatically. We believe that our approach provides a practical improvement to the highly popular NMF algorithm
without incurring much additional computational cost.

Acknowledgements
We like to thank to Nick Whiteley, Tuomas Virtanen and Paul Peeling for fruitful discussion and for their comments on
earlier drafts of this paper.

A Standard Distributions in Exponential Form, their sufficient statistics andentropies• GammaG(λ; a, b) ≡hλiG =
HG [λ] ≡µ¶1
exp +(a − 1) log λ − λ − log Γ(a) − a log bbabhlog λiG = Ψ(a) + log(b)− hlog GiG = −(a − 1)Ψ(a) + log b + a + log Γ(a)Here, Ψ denotes the digamma function defined as Ψ(a) ≡ d log Γ(a)/da.
• Poisson
PO(s; λ) = exp(s log λ − λ − log Γ(s + 1))
hsiPO = λ• Multinomial
M(s; x, p) = δ(x −XÃsi ) exp log Γ(x + 1) +i

hsi iMIX!
(si log pi − log Γ(si + 1))

i=1

= xpi

Here, s = {s1 , s2 , . . . , sI } and p = {p1 , p2 , . . . , pI } and p1 + p2 + · · · + pI = 1. Here, pi , i = 1 . . . I are the cell
probabilities and x is the index parameter where s1 + s2 + · · · + sI = x. The entropy is given as:
HM [sν,1:I,τ ] =− log Γ(xν,τ + 1) −IX

hsν,i,τ i log pν,i,τ

i=1

+

I
X

hlog Γ(sν,i,τ + 1)i −i=1*
log δ(xν,τ −

X
i

13

+
sν,i,τ )A closed form expression for the entropy is not known due to hlog Γ(s + 1)i terms, but asymptotic expansions exist[8, 12]. Computationally efficient sampling from a multinomial distribution is not trivial, see [7] for a comparison ofvarious methods and detailed discussion of tradeoffs.

A.1 Summary of the Generative Model
Indicesi = 1 . . . I Source index
ν = 1 . . . W Row (frequency bin) index
τ = 1 . . . K Column (time frame) indextν,i Template variable at ν’th row of the i’th source∼ G(tν,i ; atν,i , btν,i /atν,i )tν,ivi,τ Excitation variable of the i’th source at τ ’th columnvi,τ

G(vi,τ ; avi,τ , bvi,τ /avi,τ )∼sν,i,τ Source variable of i’th source at ν’th row (frequency bin) and τ ’th column (time frame)sν,i,τ

∼PO(sν,i,τ ; tν,i vi,τ )xν,τ Observation at ν’th row (frequency bin) and τ ’th column (time frame)Xxν,τ ∼sν,i,τ
i

A.2 Expression of the full joint distribution
Here, φ ≡ p(X, S, T, V |Θ) = p(X|S)p(S|T, V )p(V |Θv )p(V |Θt )XXXlog φ =(−tν,i vi,τ + sν,i,τ log(tν,i vi,τ ) − log Γ(sν,i,τ + 1))
ν

+

ν

+

X

sν,i,τ )i+(atν,i − 1) log tν,i −

atν,itν,i − log Γ(atν,i ) − atν,i log(btν,i /atν,i )btν,i+(avi,τ − 1) log vi,τ −avi,τ
vi,τ − log Γ(avi,τ ) − avi,τ log(bvi,τ /avi,τ )
bvi,τ

i

XX
τ

log δ(xν,τ −

τ

XX
ν

+

τ

i

XX

i

A.3 The Variational Bound
The variational bound in (13) can be written as
LX (Θ)≡ log p(X|Θ) ≥ hlog φiq + H[q] = BV B

where the energy term is given by the expectation of the expression in section A.2 and H[q] denotes the entropy of thevariational approximation distribution q where the individual entropies are defined in section A.
XXXXXX
HG [vi,τ ]HM [sν,1:I,τ ] +HG [tν,i ] +H[q] = − hlog qi =
ν

ν

τ

i

i

τ

One potential problem is that this expression requires the entropy of a multinomial distribution for which there is no
known simple expression. This is due to terms of form hlog Γ(s + 1)i where only asymptotic expansions are known.

14Fortunately, the difficult terms in the energy term can be canceled by the corresponding terms in the entropy term and one
obtains the following expression that only depends on known sufficient statistics.
XXX
B = −
htν,i i hvi,τ i
ν

+

τ

ν

+

i

XX

Ã

hlog tν,i i atν,i − 1 +!
hsν,i,τ i

+

τ

i

XX
ν

X

−

i

XX
τ

Ã
hlog vi,τ i avi,τ − 1 +X!
hsν,i,τ i

ν

i

atν,i
htν,i i − log Γ(atν,i ) − atν,i log(bν,i /atν,i )
bν,i

XX

avi,τ
− v hvi,τ i − log Γ(avi,τ ) − avi,τ log(bvi,τ /avi,τ )
bi,τ
τ
i
Ã
!
I
XX
X
− log Γ(xν,τ + 1) −+hsν,i,τ i log pν,i,τ

+

ν

τ

ν

i

i

τ

i=1

XX¡
¢
t
t
t
t
t
+
−(αν,i− 1)Ψ(αν,i) + log βν,i+ αν,i+ log Γ(αν,i
)
XX¡
¢
v
v
v
v
v
+
−(αi,τ
− 1)Ψ(αi,τ) + log βi,τ+ αi,τ
+ log Γ(αi,τ)
After some careful manipulations, the following expression is obtained where log L denotes here element wise logarithmof matrix L.
XX
B =(−Et Ev − log Γ(X + 1))
ν

+
+
+
+
+

τ

XX
ν

τ

ν

i

ν

i

i

τ

i

τ

XX
XX
XX
XX

−X .∗ (((Lt .∗ log(Lt ))Lv + Lt (Lv .∗ log(Lv )))./(Lt Lv ) − log(Lt Lv ))−(At ./Bt ) .∗ Et − log Γ(At ) + At .∗ log(At ./Bt )αt )α t .∗(log β t + 1) + log Γ(α−(Av ./Bv ) .∗ Ev − log Γ(Av ) + Av .∗ log(Av ./Bv )αv )α v .∗(log β v + 1) + log Γ(α(28)

A.4 Handling Missing Data and MAP estimation
When there is missing data, i.e., when some of the xν,τ are not observed, computation is still straightforward in ourframework and can be accomplished by a simple modification to the original algorithm. We first define a mask matrix
M = {mν,τ }, same size as X where(
0 xν,τ is missingmν,τ =(29)
1 otherwise
Using the mask variables, the observation model with missing data can be written asYm(p(xν,τ |sν,1:I,τ )p(sν,1:I,τ |tν,1:I , v1:I,τ )) ν,τ
p(X|S)p(S|T, V ) =ν,τ

The prior is not affected. Hence, we merely replace the first two lines of the expression for the full joint distribution (givenin the appendix A.2) as
X
XX
mν,τ
(−tν,i vi,τ + sν,i,τ log(tν,i vi,τ ) − log Γ(sν,i,τ + 1))
log φ =
ν

+

τ

XX
ν

∝

mν,τ log δ(xν,τ −X

τConsequently, it is easy to see thatÃ
q(vi,τ ) ∝i

exp (avi,τ +sν,i,τ ) + . . .

iXÃmν,τ hsν,i,τ i − 1) log vi,τ −

ν
¡
¢
v
v
G vi,τ ; αi,τ, βi,τ

15!
!
avi,τ Xmν,τ htν,i i vi,τ
+
bvi,τ
ν

v
αi,τ
≡ avi,τ +

Ã

X

v
βi,τ
≡

mν,τ hsν,i,τ iν!−1avi,τ X+mν,τ htν,i ibvi,τ
νBy a derivation analogous to one detailed in section 3.3, we see that the excitation update equations in Algorithm 1, line
4, can be written using matrix notation as= Lv .∗(L>t ((M .∗ X)./(Lt Lv )))¡¢= Av + Σv
β v = 1./ Av ./Bv + Et> M

Σv
αv

The update rules for the templates are similar. Note that when there is no missing data, we have M = 1W 1>
K whichgives the original algorithm. The bound in (28) can also be easily modified for handling missing data. We merely replace
X ← M .∗ X and the first term Et Ev ← M .∗ Et Ev .
We conclude this subsection by noting that the standard NMF update equations, given in Equations (10) and (11) canbe also rewritten to handle missing data:
V (n+1)T (n+1)= V (n) .∗(T > ((M .∗ X) ./(T V ))) ./(T > M)
= T (n) .∗(((M .∗ X) ./(T V ))V > ) ./(MV > )Here, the denominator has to be nonzero. Similarly, an iterative conditional modes (ICM) algorithm can be derived tocompute the maximum a-posteriori (MAP) solution
V (n+1)T (n+1)(Av + V (n) .∗(T > ((M .∗ X) ./(T V )))) ./(Av ./ Bv + T > M)
(At + T (n) .∗(((M .∗ X) ./(T V ))V > )) ./(At ./ Bt + MV > )

=
=

(30)
(31)

Note that when the shape parameters go to zero, i.e. At , Av → 0, we obtain the maximum likelihood NMF algorithm.

A.5 Hyperparameter optimisation
The hyperparameters Θ = (Θt , Θv ) can be estimated by maximising the bound in 13. Below, we will derive the results
for the excitations; the results for templates are similar. The solution for shape parameters involves finding the zero of afunction f (a) − c where
f (a)
a∗

= log a − Ψ(a) + 1= f −1 (c)The solution can be found by Newton’s method by iteration of the following fixed point equationa(n+1)= a(n) −log(a(n) ) − Ψ(a(n) ) + 1 − cf (a(n) ) − c
(n)
=
a
−
= a(n) − ∆(n)f 0 (a(n) )1/a(n) − Ψ0 (a(n) )It is well known that Newton iterations can diverge if started away from the root. Occasionally, we observe that a canbecome negative. If this is the case, we set ∆(n) ← ∆(n) /2, and try again. The digamma function Ψ function and its
derivative Ψ0 are available in numeric computation libraries (e.g. in Matlab as psi(0,a) and psi(1,a), respectively).
The derivation of the hyperparameter update equations is straightforward:∂B∂avi,τ

=hlog vi,τ i −1hvi,τ i − Ψ(avi,τ ) − log bvi,τ + log avi,τ + 1 = 0
bvi,τ

ci,τ

= log avi,τ − Ψ(avi,τ ) + 1

ci,τ

≡

hvi,τ i− (hlog vi,τ i − log bvi,τ )bvi,τ
∂B∂bvi,τ
bvi,τ

=

avi,τ
1
hvi,τ i − avi,τ v = 0v
2(bi,τ )bi,τ= hvi,τ iTying parameters across τ as avi = avi,τ and bvi = bvi,τ yields
∂B
∂avi

=

ci

=

ci

=

X
τ

hlog vi,τ i −X 1
X
log bvi,τ + K = 0
hvi,τ i − KΨ(avi ) −
v
b
i,τ
τ
τ

log avi − Ψ(avi ) + 1
Ã
!
¢
1 X hvi,τ i ¡v− hlog vi,τ i − log bi,τ
K τ
bvi,τ
16

∂B
∂bvi
bvi

=

X avi,τ
τ
P

1 X v
ai,τ

hvi,τ i − v(bvi )2bi

τv
τ ai,τ hvi,τ i
P v
τ ai,τ

=Tying parameters across τ and i, av = avi,τ and bv = bvi,τ yields
c = log av − Ψ(av ) + 1
Ã
!
1 X X hvi,τ iv− (hlog vi,τ i − log bi,τ )c =KI τ i
bvi,τ
∂B
∂bv

=

bv

=

X X avi,ττ

1 XX vhvi−ai,τ
(bv )2bvi i τ i,τi
PP v
a hvi,τ ii
P τPi,τ vi
τ ai,τ

The derivation of the template parameters is exactly analogous. We can express the update equations once again incompact matrix notation
Z

C

A(n+1)v

Bv(n+1)(n)Ev(n) ./Bv(n) − log(L(n)
v ./Bv ))Z
Not tied
(Z1 )/K
Tie columns (over τ )K←>(1I Z)/I
Tie rows (over i) >(1I Z1K )/(KI) Tie all (over τ and i)←

SolveByNewton(A(n)
v , C)

(n)

Ev



(((A(n) .∗ E (n) )1 )./(A(n) 1 ))1>
v
v
v
K
K
K
←
(n)
(n)
> (n)
>

))
))./(1
A
.∗
E
1
((1
(A
v
v
v
I

I
I


1 ((1> (A(n) .∗ E (n) )1 )./(1> A(n) 1 ))1>
v
v
v
K
K
I
K
I
I

←

Not tiedTie columns (over τ )Tie rows (over i)
Tie all (over τ and i)

Here, we assume SolveByNewton(A0 , C) is a matrix valued function that finds root Ci,j = f (Ai,j ) for each elementof A, starting from the initial matrix A0 . If C is a scalar or vector, it is repeated over the missing index to implement
parameter tying. For example, if C is a I × 1 vector and A0 is I × K, we assume Ci = ci,τ for all τ = 1 . . . K
and the output is the same size as A0 . This is only a notational convenience, an actual implementation can be achieved
more efficiently. Again, the implementation of the template parameters is exactly analogous; merely replace above thesubscripts as v ← t, (i, τ ) ← (ν, i) and (I, K) ← (W, I).

References[1] P. Smaragdis W. Wang A. Cichocki, M. Mrup and R. Zdunek. Advances in nonnegative matrix and tensor factorization. Computational Intelligence and Neuroscience, 2008, 2008.
[2] Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
[3] W. Buntine. Variational extensions to EM and multinomial PCA. In ECML, LNAI 2430, pages 23–34. Springer, 2002.
[4] K. Butler and M. Stephens. The distribution of a sum of binomial random variables. Prepared for the Office of Naval Research
467, Stanford University, April 1993.
[5] A. T. Cemgil, P. Peeling, O. Dikmen, and S. J. Godsill. Prior structures for time-frequency energy distributions. In Proc. of IEEE
Workshop on Applications of Signal Processing to Audio and Acoustics, October 2007.
[6] S. Chib. Marginal likelihood from the gibbs output. JASA, 90(432):1313–1321, Dec. 1995.
[7] C. S. Davis. The computer generation of multinomial random variates. Computational Statistics and Data Analysis, 16:205–217,
1993.
[8] P. Flajolet. Singularity analysis and asymptotics of bernoulli sums. Theor. Comput. Sci., 275:371–387, 1999.
[9] Z. Ghahramani and M. Beal. Propagation algorithms for variational Bayesian learning. In Neural Information Processing Systems
13, 2000.

17

[10] W. R. Gilks, S. Richardson, and D. J Spiegelhalter, editors. Markov Chain Monte Carlo in Practice. CRC Press, London, 1996.
[11] P. O. Hoyer. Non-negative matrix factorization with sparseness constraints. J. Mach. Learn. Res., 5:1457–1469, 2004.
[12] P. Jacquet and W. Szpankowski. Entropy computations via analytic depoissonization. IEEE Transactions on Information theory,45(4), May 1999.
[13] H. Kameoka. Statistical Approach to Multipitch Analysis. PhD thesis, University of Tokyo, 2007.
[14] J. F. C. Kingman. Poisson Processes. Oxford Science Publications, 1993.
[15] J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama. On the relation between divergence-based minimization and maximumlikelihood estimation for the i-divergence. Submitted to IEEE Trans. Signal Processing, 2008.
[16] D. D. Lee and H. S. Seung. Learning the parts of objects with nonnegative matrix factorization. Nature, 401:788–791, 1999.
[17] D. D. Lee and H. S. Seung. Algorithms for non-negative matrix factorization. In NIPS, pages 556–562, 2000.
[18] J. S. Liu. Monte Carlo strategies in scientific computing. Springer, 2004.
[19] F. A. Quintana, J. S. Liu, and G. E. del Pino. Monte Carlo EM with importance reweighting and its applications in random effects
models. Computational Statistics and Data Analysis, 29:429–444, 1999.
[20] M.V. Shashanka, B. Raj, and P. Smaragdis. Sparse overcomplete latent variable decoposition of counts data. Neural Information
Processing Systems (NIPS), December 2007.
[21] M. A. Tanner. Tools for Statistical Inference: Methods for the Exploration of Posterior Distributions and Likelihood Functions.
Springer, New York, 3rd edition, 1996.
[22] T. O. Virtanen, A. T. Cemgil, and S. J. Godsill. Bayesian extensions to nonnegative matrix factorisation for audio signal modelling.
In Proc. of IEEE ICASSP 08, Las Vegas, 2008. IEEE.
[23] Tuomas Virtanen. Sound Source Separation in Monaural Music Signals. PhD thesis, Tampere University of Technology, November 2006.

18