CONSTRAINED AND REGULARIZED VARIANTS OF NON-NEGATIVE MATRIXFACTORIZATION INCORPORATING MUSIC-SPECIFIC CONSTRAINTS
Hirokazu Kameoka1,2, Masahiro Nakano2, Kazuki Ochiai1, Yutaka Imoto1, Kunio Kashino2, Shigeki Sagayama1
1

Graduate School of Information Science and Technology, The University of Tokyo
2
NTT Communication Science Laboratories, NTT Corporation
ABSTRACTMusic spectrograms typically have many structural regularities thatcan be exploited to help solve the problem of decomposing a given
spectrogram into distinct musically meaningful components. In thispaper, we introduce new variants of the non-negative matrix factorization concept that incorporate music-speciﬁc constraints.
Index Terms— Non-negative Matrix Factorization, Music transcription, Source separation
1. INTRODUCTION
Non-negative Matrix Factorization (NMF) [1] is a relatively recent
technique, which is used to decompose a data matrix Y into two
factors H and U with non-negative entries. Owing to its simplicity and intuitive decomposition, it has attracted a lot of attention in
many scientiﬁc and engineering areas in recent years. In the area ofmusic signal processing, one successful approach involves applying
NMF to a magnitude spectrogram (time-frequency representation)interpreted as a non-negative matrix [2], where the spectrogram Y
is factorized into the product of a basis matrix H consisting of spectrum atoms and an activation matrix U consisting of time-varying
amplitudes associated with these atoms. An important feature ofthis approach is that it is capable of ﬁnding a ﬁnite set of spectrum
atoms that are considered to be the dominant elements constituting
the observed spectrogram, in an unsupervised manner. This ability
is proven to be very powerful and has enabled the NMF approach tobe applied with notable success to many tasks including monauralsource separation, noise reduction, music transcription, bandwidth
expansion, and missing data imputation.
NMF usually considers decompositions that are approximativein nature. That is, the decomposition is performed so that the product X = HU should approximate the original data Y as well as
possible. Although the approximation error can be minimized as desired by using a larger number of spectrum atoms, the obtained basis
vectors will be less likely to represent the spectra of meaningful audio events as the number of bases becomes larger. As an extreme
example, if we use the same number of bases as the number of thefrequency bins, we obtain an exact reconstruction Y = X = HU
under a trivial solution H = I (where I is an identity matrix) andU = Y , which is no longer a meaningful decomposition. On theother hand, if we use a relatively small number of bases, the chanceof obtaining a meaningful decomposition will increase, but the resulting decomposition will become less accurate. This illustrates thefact that, for classical NMF, there is a trade-off between the accuracyof the reconstruction and the meaningfulness of the decomposition.
To achieve an accurate approximation and thus explain the data with
the NMF model as well as possible, we must use many bases. To obtain a meaningful decomposition while using many bases, we must

978-1-4673-0046-9/12/$26.00 ©2012 IEEE

incorporate reasonable assumptions other than non-negativity to further constrain the model HU .
Another issue concerning NMF is the local optimum problem.
Typically most cost functions used in NMF are difﬁcult to optimize
analytically with respect to H and U . In addition, they are notjointly convex in both of the arguments H and U . Thus, many existing algorithms developed for NMF are guaranteed to converge to oneof the stationary points but not necessarily to a global optimum. In
particular, when applying NMF to music signals, we would want toavoid local optimum solutions that are “musically” unacceptable. A
good strategy to guarantee that we obtain a musically likely solution
would be to incorporate music-speciﬁc constraints into the model
HU or into the optimization problem.
Using these considerations as a basis, we have been concernedwith developing improved variants of NMF, tailored speciﬁcally formusic signals, by utilizing reasonable assumptions that we can makeabout music spectrograms. This paper introduces our ongoing work
along with some new ideas on NMF variants incorporating musicspeciﬁc constraints.
2. NMF WITH TIME-VARYING BASIS SPECTRA
2.1. Motivation
When applying the classical NMF to music spectrograms, we may
expect the spectra of a single note produced by a musical instrument to be represented using a single basis spectrum scaled by timevarying amplitudes. However, its variations in time are actuallymuch richer. A piano note would be more accurately characterized
by a succession of several basis spectra corresponding to, for example, “attack”, “decay”, “sustain” and “release” segments. As anotherexample, singing voices and string instruments feature a particular
musical effect, vibrato, which can be characterized by its “depth”(the range of pitch variation), and its “speed” (the rate at which thepitch varies). Learning such time-varying spectra with the classical NMF would require the use of a large number of bases, andsome postprocessing to group the bases into single events. However, as indicated earlier, blindly increasing the number of bases will
not necessarily give meaningful decompositions at the NMF stage.
If we want to increase the number of bases while maintaining themeaningfulness of the decomposition, we will need some additional
constraints and/or regularization to correspondingly reduce the degree of freedom of the model. In this section, we brieﬂy review themodel described in [3, 4], which is designed to undertake the decomposition and group the basis spectra simultaneously based on
the concept of time-varying basis spectra.
2.2. Model
Let us begin by dividing the bases into groups each of which we
expect to correspond to a single note produced by a particular musical instrument. The NMF model X = (Xω,t )Ω×T can then be

5365

ICASSP 2012

expressed as
Xω,t =

Ik
K X
X

(i)

(i)

Hω,k Uk,t ,

(a)

(b)

(c)

(d)

(1)

k=1 i=1

where k is the group index, i is the basis index in each group, Ik is
the number of basis spectra in group k, and ω and t are frequency
and time indices, respectively.
P Note that so far this is the same as theclassical NMF model with k Ik bases. We shall now introduce thefollowing assumptions to constrain the model:
1. For each time t, exactly one basis spectrum in each group isallowed to be activated.
2. For each group k, the order in which the basis spectra appear
is governed by a Markov chain.
By incorporating assumption 1, the model can be rewritten as
Xω,t =K
X(Z)

Hω,kk,t Uk,t ,

(2)

k=1

where Zk,t ∈ {1, . . . , Ik } denotes a hidden state variable indicating
which basis spectrum is supposed to be activated at time t. Noticethat the superscript i is dropped from U in (2) as it is no longer nec(i)essary since we are assuming Uk,t = 0 for i = Zk,t . From assumption 2, the path of the state variables Zk,1 , . . . , Zk,T is governed bya state transition probability p(Zk,t = a|Zk,t−1 = b) = πk,a,b .
Now, let D(x, y) be a discrepancy measure between x and ysuch that D(x, y) ≥ 0 and D(x, y) = 0 only if x = y. We can thendeﬁne a goodness-of-ﬁt measure between the observed spectrogram
Y and the current NMF model X
X
D(Yω,t , Xω,t ).
(3)
J (Y , X) =ω,t

If we deﬁne D(x, y) as the I-divergence [5]D(x, y) = y logy− (y − x),x(4)then minimizing J (Y , X) with respect to X is known to amountto maximizing the Poisson likelihoodp(Y |X) =Y

Poisson(Yω,t ; Xω,t ),

(5)

ω,t

where Poisson(y; x) = xy e−x /y!. From assumption 2, p(Z) is
deﬁned byp(Z) =Ykp(Zk,1 |πk0 )TY−1p(Zk,t |Zk,t−1 , ı k ),

(6)

t=2

where ı k = {πk,a,b }1≤a≤Ik ,1≤b≤Ik .
Readers are referred to [3, 4] for detailed derivations of the parameter inference algorithms under this setting: [3] describes an algorithm for ﬁnding the Maximum A Posteriori (MAP) estimates of
H, U and Z. [4] describes a nonparametric Bayesian formulationof the proposed NMF variant, making it possible to infer the numberof groups, K, and the number of basis spectra in each group, Ik ,along with the posterior distributions of H, U and Z.

5366

Fig. 1. Portions of the spectrograms of a mixture of 3 voices (a),estimated components corresponding to A (b), F (c), and D (d).
2.3. Related work
The model described above can be viewed as a factorial hidden
Markov model (HMM) [6]. Recently, several authors have independently proposed modeling spectrograms using factorial HMMs
[7, 8]. Let us consider a factorial HMM where the latent component generated from the k-th HMM is denoted by Ck,ω,t , andthe observed data Yω,t are assumed to be the sum of Ck,ω,t over k.
(Z)
Now, if we assume Ck,ω,t ∼ Poisson(Ck,ω,t ; Hω,kk,t Uk,t ), we canconﬁrm that the assumed factorial HMM amounts to the model de(Z)
scribed in 2.2. If we assume Ck,ω,t ∼ NC (Ck,ω,t ; 0, Hω,kk,t Uk,t ),2where NC (y; 0, x) ∝ (1/x)e−|y| /x , then it amounts to the model
presented in [7].
2.4. ExperimentFig. 1 shows the decomposition achieved by the present model. Here
we used a signal composed of 3 notes (DɼFɼA): ﬁrst, each noteis played alone in turn, then all two note combinations are played
and ﬁnally all the notes are played simultaneously. As can be seenfrom Fig. 1, the present model is capable of grouping together the
spectra originating from one voice even though there is a variationin pitch (vibrato).
3. NMF WITH BEAT STRUCTURE CONSTRAINT
3.1. Motivation
Music is highly structured in terms of the temporal regularity underlying the onset occurrences of notes. In general, the time between
consecutive onsets corresponds to multiples and fractions of the beatperiod, with small deviations in timing and tempo. As we wanted
each basis spectrum to correspond to a single note, the onsets ofthe basis spectra should have this rhythmic structure, which can beeffectively used to constrain the activation matrix U . In this section, we propose the introduction of a constrained variant of NMF,
in which the activation matrix model is parameterized by note onsets, beat locations and tempo. To our knowledge, this is the ﬁrst
NMF approach that explicitly incorporates constraints derived fromthe rhythmic structure of music. In the following, we describe only
the basic idea. For more details, please refer to [9].
3.2. Model
Let us consider a standard NMF modelX
Xω,t =
Hω,d Ud,t .

(7)dBased on the rhythmic structure of music, we make the following
assumptions as regards constraining the activity function Ud,t :First, from assumption 1, Ud,t should be written as
Ud,t =

Ld
X

Vd,l,t ,

C6
C5

Note

1. Each activity function consists of local activity patterns,
called “objects”, each of which we expect to correspond to asingle note activation.
2. Each object is characterized by a fast/slow rise at the onset
time followed by a 1 continuous contour.
3. The onset of each object is likely to be located on the multiples or fractions of the beat period.
4. The beat period varies gradually over time.

C4
C3
C2
3

(8)

6

9

12

Time [s]15

Fig. 2. Estimates of pitch, onsets/offsets, and beat locations of eachnote obtained with the present model applied to polyphonic data.

l=1

where Vd,l,t is the l-th object and Ld is the number of objects inthe d-th activity function. To incorporate assumption 2 into Vd,l,t ,we introduce a parametric model, which is expressed as the sum ofGaussians [11]Vd,l,t =MX
vd,l wd,l,m −(t−(m−1)φd,l −τd,l )2 /2φ2d,l√e,2πφd,lm=1(9)where vd,l , τd,l and φd,l are the total energy, onset time and parameter related to the duration of the l-th object, respectively.
wd,l,1 , . . . , wd,l,M are the weights associated with the M Gaussians that sum to unity, which determine the shape of the object.
Under this constrained model Xω,t , weQ assume a Poisson likelihood
as with 2.2, p(Y |H, v, w, ﬁ , ﬃ) = ω,t Poisson(Yω,t ; Xω,t ). To
avoid overﬁtting the shape of each object, we
Q place a Dirichlet priorover wd,l = {wd,l,m }1≤m≤M , p(w) = d,l Dirichlet(wd,l ; ¸),
Q xi −1and ¸ denotes the mostwhere Dirichlet(y; x) ∝i yi
expected shape of the object. To promote sparsity of the activity function, we place a generalized Gaussian prior over vd,l ,Qpp(v) = d,l GN (vd,l ; 0, λ, p), where GN (y; 0, λ, p) ∝ e−λ|y|and 0 < p < 2. Furthermore, to ensure that each basis spectrum
has a harmonic structure of a particularpitch, we shall also place a
Q
Gamma prior over H, p(H) = ω,d Gamma(Hω,d ; β H̄ω,d+1, β),where Gamma(x; a, b) ∝ xa−1 e−bx .
Next, to impose Assumption 3, we ﬁrst introduce a set of hyperparameters, = {ψj }1≤j≤J , where ψj corresponds to the timeinterval between the j-th and (j − 1)-th beat locations. With these
hyperparameters, we can design a Gaussian prior distribution over
the onset parameter τd,lY
p(ﬁ | ) =N (τd,l ; ρl , ν 2 ),
(10)
d,l
l/I

ρl =

X

ψj + (l/I − l/I ) ψl/I+1 ,

(11)

j=1
2

2

where N (y; x, σ 2 ) ∝ e−(y−x) /2σ , ρl denotes the most expected
location of the onset of the l-th object, ν 2 is the variance of the Gaussian indicating how much τd,l is allowed to deviate from ρl , I is thenumber of divisions per beat, and x denotes the largest integer not
greater than x. To impose Assumption 4, it is convenient to place a
Gaussian chain hyperprior overp( ) = p(ψ1 )JYp(ψj |ψj−1 ),

(12)

j=2
1 This assumption is made according to the suggestion that continuity constraints can improve the performance of NMF-based source separation [10].

5367p(ψj |ψj−1 ) = N (ψj ; ψj−1 , σ 2 ).

(13)

Putting altogether, the posterior density of the unknown parameters is given byp(H, v, w, ﬁ , ﬃ, |Y )∝ p(Y |H, v, w, ﬁ , ﬃ)p(H)p(v)p(w)p(ﬁ | )p( ).

(14)

The MAP estimates of H, v, w, ﬁ , ﬃ, can be found through an
iterative algorithm consisting of parameter updates given in closed
form. See [9] for its derivation.
3.3. ExperimentFig. 2 shows an example of the application of the present model topolyphonic piano music. With the present model, we were able toestimate the pitch and onset of each note with a high detection rate,compared with the baseline NMF model [9].
4. NMF WITH TIMBRAL CLUSTERING CRITERION
4.1. Motivation
In general, each piece of music is typically played by only a handfulof musical instruments or sung by one or a few singers. Thus, theconstituent sounds contained in a single piece of music can probably
be grouped into a reasonably small number of clusters in a certain
feature space that best represents the timbral aspect of the instruments or human voice. Based on this expectation, we consider incorporating a timbral clustering criterion for H into the objectivefunction for NMF to effectively constrain the solution space of H.
4.2. Model
Let us begin by considering the standard NMF model given in (7).
We are concerned with the problem of ﬁnding a decomposition suchthat Y
HU in which the basis spectra are forced to be clusteredinto some number K of clusters in a certain feature space, where
we shall suppose that the value of K is given. Here we expect eachgroup k to consist of a set of basis spectra that are timbrally consistent, meaning that all the basis spectra assigned to the same cluster are expected to be mapped onto a single point in some timbral
space. Our preliminary experiments on an instrument identiﬁcation
task using monophonic data revealed that the mel-frequency cepstral coefﬁcient (MFCC) is a reasonably relevant feature for robustly
identifying the kinds of instruments, and so we choose here to assume the MFCC space as the feature space. The MFCC of the d-thbasis spectrum is, by deﬁnition, given by
X
X
cm,n log
fn,ω Hω,d ,(15)
Hm,d =nωwhere {fn,ω }1≤ω≤Ω denotes the n-th triangular mel-ﬁlter in themel-frequency ﬁlterbank and {cm,n }1≤n≤N denotes a set consisting of the discrete cosine transform coefﬁcients. It is worth notingthat this expression reduces to the mel-spectrum when cm,n = δm,n(where δ denotes Kronecker’s delta), and the log-power spectrum
when fn,ω = δn,ω and cm,n = δm,n . For each basis spectrumd, we introduce a corresponding set of binary indicator variablesrd,k ∈ {0, 1}, describing to which of the K clusters the d-th basis
spectrum is assigned, so that if the d-th basis spectrum is assigned
to cluster k then rd,k = 1, and rd,k = 0 for k = k. By using —k
to denote a prototype vector associated with the k-th cluster, we candeﬁne a cost function representing the sum of the distances of eachfeature vector to the prototype vector of the assigned cluster
XX‚‚2
R(H, r, —) =rd,k ‚Hd − —k ‚2
=

d

k

d

k

XX

X˛
˛
˛Hm,d − μm,k ˛2 ,
rd,k

(16)

m

where r = {rd,k }1≤d≤D,1≤k≤K and — = {—k }1≤k≤K . Hd :=(H1,d ,. . . , HM,d )T denotes the feature vector of the d-th basis spectrum. Our goal is to ﬁnd values for H, U , r, and — so as to minimize
J (Y , HU ) + λR(H, r, —),(17)

where λ > 0 is a regularization parameter. We can do this through an
iterative procedure in which each iteration involves four successive
steps corresponding to successive optimizations with respect to H,
U , r and —. Although the optimization with respect to H is mathematically intractable, we can derive a closed form update equationthat guarantees a certain decrease in the objective function by using
an auxiliary function approach. To do this, the ﬁrst step is to deﬁne
an upper bound function for the objective function (17). Owing tospace limitations, we only show an upper bound function R+ for the
regularization term R without proof:
R+ (H, r, —, q, j, , ‰, ﬃ) =
(18)
„X
«
X
X
X
ρn,ω,d 2rd,kAn,d+ h(ξn,d )
fn,ω Hω,dfn,ω Hω,d
n
ω
ω
H

d,k

X

|Bn,d,k | X
fn,ω Hω,d
φn,d
n
ω
d,k
X
X
X
−
rd,k
1[Bn,d,k < 0]|Bn,d,k |νn,ω,d log Hω,d ,
+

d,k

rd,k

X

1[Bn,d,k ≥ 0]

n

ω

H

where = denotes equality up to a term independent of H, 1[·]is the indicator function that takesPthe value 1 if its argument2
is truem cm,n /βm,n,d , Bn,d,k =P and 0 otherwise, An,d =−2 m αm,n,d μm,k cm,n /βm,n,d , and h(x) = 2 log(x)/x + 1/x2 .
αm,n,d ρn,ω,d , νn,ω,d , ξn,d and φn,d are auxiliary variables.
βm,n,dP
is an arbitrary parameter satisfying βm,n,d > 0 and n βm,n,d = 1.
An exact bound is achieved whenX
αm,n,d μm,k = cm,n log
fn,ω Hω,dω+ βm,n,d (Hm,d − μm,k ),fn,ω Hω,dρn,ω,d = νn,ω,d = X,
fn,ω Hω ,dξn,d = φn,d =

X

(19)
(20)

ω

fn,ω Hω,d .

(21)

ω

5368

Fig. 3. An example of the decomposition obtained with the proposed
regularized NMF.
4.3. ExperimentFig. 3 shows an example of the basis spectra/activations obtained
from polyphonic music played on a piano and a bass guitar (its spectrogram is shown lower right), where we assumed K = 2. The basis
spectra/activations corresponding to clusters 1 and 2 are different
colors. In this example, we were able to group together the spectra
originating from the same instrument automatically with the present
method, even though H and U were initialized randomly.
5. SUMMARY
This paper introduced our ongoing work along with some new ideas
on constrained variants of NMF that incorporate structural regularities underlying music. In the future, we plan to combine the ideas
introduced in this paper to construct a uniﬁed model. While we havefocused solely on the physical aspects of the regularities in music, weare also concerned with incorporating symbolic regularities, such as
tonality and harmony.
6. REFERENCES[1] D. D. Lee and H. S. Seung, “Learning the parts of objects with nonnegative matrix factorization,” Nature, vol. 401, pp. 788–791, 1999.
[2] P. Smaragdis and J. C. Brown, “Non-negative matrix factorization formusic transcription,” in Proc. WASPAA’03, 2003, pp. 177–180.
[3] M. Nakano, J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama, “Nonnegative matrix factorization with Markov-chained bases for modelingtime-varying patterns in music spectrograms,” in Proc. LVA/ICA’10,
2010, pp. 149–156.
[4] M. Nakano, J. Le Roux, H. Kameoka, T. Nakamura, N. Ono, and
S. Sagayama, “Bayesian nonparametric spectrogram modeling basedon inﬁnite factorial inﬁnite hidden Markov model,” in Proc. WASPAA’11, 2011, to appear.
[5] I. Csiszár, “I-divergence geometry of probability distributions and minimization problems,” The Annals of Probability, vol. 3, no. 1, pp. 146–
158, 1975.
[6] Z. Ghahramani and M. I. Jordan, “Factorial hidden Markov models,”
Machine Learning, vol. 29, pp. 245–273, 1997.
[7] A. Ozerov, C. Févotte, and M. Charbit, “Factorial scaled hidden
Markov model for polyphonic audio representation and source separation,” in Proc. WASPAA’09, 2009, pp. 121–124.
[8] G. J. Mysore, P. Smaragdis, and B. Raj, “Non-negative hidden Markov
modeling of audio with application to source separation,” in Proc.
LVA/ICA’10, 2010, pp. 140–148.
[9] K. Ochiai, H. Kameoka, and S. Sagayama, “Explicit beat structure
modeling for non-negative matrix faxtorization-based multipitch analysis,” in Proc. ICASSP’12, 2012, to appear.
[10] T. Virtanen, “Monaural sound source separation by nonnegative matrix
factorization with temporal continuity and sparseness criteria,” IEEE
Trans. ASLP, vol. 15, no. 3, pp. 1066–1074, 2007.
[11] H. Kameoka, T. Nishimoto, and S. Sagayama, “A multipitch analyzer
based on harmonic temporal structured clustering,” IEEE Trans. ASLP,
vol. 15, no. 3, pp. 982–994, 2007.