Nonnegative Matrix Factorization
with Markov-Chained Basesfor Modeling Time-Varying Patternsin Music Spectrograms
Masahiro Nakano1 , Jonathan Le Roux2 , Hirokazu Kameoka2 , Yu Kitano1 ,Nobutaka Ono1 , and Shigeki Sagayama1
1

Graduate School of Information Science and Technology, The University of Tokyo
2
NTT Communication Science Laboratories, NTT Corporation

Abstract. This paper presents a new sparse representation for polyphonic music signals. The goal is to learn the time-varying spectral patterns of musical instruments, such as attack of the piano or vibrato ofthe violin in polyphonic music signals without any prior information. We
model the spectrogram of music signals under the assumption that they
are composed of a limited number of components which are composedof Markov-chained spectral patterns. The proposed model is an extension of nonnegative matrix factorization (NMF). An eﬃcient algorithm
is derived based on the auxiliary function method.
Keywords: Nonnegative matrix factorization, Sparse signal representation, Source separation, Markov chain, Auxiliary function.

1Introduction

The use of sparse representation in acoustic signal processing has been a veryactive area of research in recent years, with very eﬀective algorithms based on
nonnegative matrix factorization (NMF) [1] and sparse coding. These are typically based on a simple linear model.
NMF, in particular, has been applied extensively with considerable successto various problems including automatic music transcription, monaural soundsource separation [2]. NMF is able to project all signals that have the same
spectral shape on a single basis, allowing one to represent a variety of phenomenaeﬃciently using a very compact set of spectrum bases. However, because NMF isalso fundamentally a dimension reduction technique, a lot of information on theoriginal signal is lost. This is in particular what happens when assuming thatthe spectrum of the note of a musical instrument can be represented through asingle spectral basis whose amplitude is modulated in time, while its variations
in time are actually much richer. For example, a piano note would be moreaccurately characterized by a succession of several spectral patterns such as“attack”, “decay”, “sustain” and “release”. As another example, singing voicesand stringed instruments feature a particular musical eﬀect, vibrato, which canV. Vigneron et al. (Eds.): LVA/ICA 2010, LNCS 6365, pp. 149–156, 2010.
c Springer-Verlag Berlin Heidelberg 2010
150

M. Nakano et al.

be characterized by its ”depth”, (the amount of pitch variation), and its ”speed”,(the speed at which the pitch varies). Learning such time-varying spectra withstandard NMF would require to use a large number of bases, and some postprocessing to group the bases into single events.
In this paper, we propose a new sparse representation, “NMF with Markovchained bases” for modeling time-varying patterns in music spectrograms. Ourmodel represents a single event as a succession of spectral patterns. The proposedmodel is presented in Section 2, together with the derivation of an eﬃcientalgorithm to optimize its parameters. We present basic experimental results inSection 3.

2

NMF with Markov-Chained Bases

2.1

Presentation of the Model

Most algorithms for unsupervised sound source separation are based on a signal
model where the magnitude or power spectrogram Y = (Yω,t )Ω×T ∈ R≥0,Ω×T ,where ω = 1, · · ·, Ω is a frequency bin index, and t = 1, · · ·, T is a time frameindex, is factorized into nonnegative parameters, H = (Hω,d )Ω×D ∈ R≥0,Ω×Dand U = (Ud,t )D×T ∈ R≥0,D×T . This can be written asDYω,t =Hω,d Ud,t ,(1)d=1where D is the number of bases hd = [H1,d , · · ·, HΩ,d ]. The term component isused to refer to one basis hd and its time-varying gain Ud,t . The bases can beconsidered as spectral patterns which are frequently observed.
Hopefully, one component should represent a single event. However, the spectrum of instrument sounds is actually in general nonstationary. Each source willthus tend to be modeled as a sum several components, leading to the diﬃcult
problem of determining which source each component belongs to. Recently, anextension of NMF where temporal activations become time/frequency activations based on a source/ﬁlter model [3] have been proposed to overcome thisproblem. However, it has not been clariﬁed whether the source/ﬁlter model isﬁt for auditory stream composed of various origin, such as the piano. For example, “attack” caused by a keystroke has energies on wide-band spectrum, while“sustain” has a harmonic structure.
In contrast with the above approaches, we focus on the hierarchical structureof the sounds produced by musical instruments, and model the spectrogram ofmusic signals under the assumption that they are composed of spectral patterns
which have themselves a limited number of Markov-chained states. Concretely,
we assume that each basis hd has Φ states, the transitions between those states
being constrained and only one state being activated at each time t.
We attempt to model the spectrogram again based on H = (Hω,φ,d )Ω×Φ×D ∈R≥0,Ω×Φ×D Here, P = (Pφ,t,d )Φ×T ×D ∈ R≥0,Φ×T ×D is binary to show which(φ)
basis is activated at time t, i.e., Pφ,t,d = 1 if hd is activated at time t, andPφ,t,d = 0 otherwise. Then, the proposed model can be written as

NMF with Markov-Chained Bases for Modeling Time-Varying Patterns

151

Frequency

state transitionTime1 23state #1 2

3activation

state #

Markov-chainedbasis

21

3

1state32activation
12

3

1 2

3stateFig. 1. Diagram of NMF with Markov-chained bases

Yω,t =

Hω,φ,d Pφ,t,d Ud,t .

(2)

d,φ

Note that P does not show the probabilities that state φ is activated at time t
and the class S of all possible P is deﬁned by the topology of the Markov chain,for example left-to-right model or ergodic model.
2.2

Problem Setting

Given an observed spectrogram, we would like to ﬁnd the optimal estimatesof H, U and P. The standard NMF algorithm proposed by Lee and Seung [1]performs the decomposition by minimizing the reconstruction error between theobservation and the model while constraining the matrices to be entry-wise nonnegative. Our estimation can also be written as optimization problem, similarly
to standard NMF. Various measures for standard NMF have been proposed. We
choose here the following β-divergence, which has been widely used, βy + (β − 1)xβ − βyxβ−1.
(3)
Dβ (y|x) =β(β − 1)Note that the above deﬁnition can be properly extended by continuity to β = 0and β = 1 [5]. Let θ denote the set of all parameters {(Hω,φ,d )Ω×Φ×D , (Ud,t )D×T ,(Pφ,t,d )Φ×T ×D }. We then need to solve the following optimization problem:minimize J (θ) = ω,t Dβ Yω,t | d,φ Hω,φ,d Pφ,t,d Ud,t(4)subject to ∀ω, φ, d, Hω,φ,d ≥ 0, ∀d, t, Ud,t ≥ 0, P ∈ S .
In this paper, we assume that the transition probabilities of the Markov-chainedbases are uniform. Thus, the cost of each of the paths is regarded as constantand ignored in Eq. (4). In the following, we will refer to the algorithm to solve
NMF with Markov-chained bases as “MNMF”.

152

M. Nakano et al.

Note that our model can also be expressed as Factorial Hidden Markov Modelfor speciﬁc βs. The case of β = 0 reduces to Factorial scaled Hidden Markov
Model which reduces to NMF with Itakura-Saito (IS) divergence and Gaussian
Scaled Mixture Model [4]. These probabilistic model may achieve the extensionof NMF with the Euclidean distance (β = 2), the generalized Kullback-Leibler
divergence (β = 1) and IS divergence (β = 0) based on a statistical approach to
NMF [5,6]. However, it is not clariﬁed whether statistical approach can apply to
NMF with various measures, such as β-divergence.
2.3

Iterative Algorithm

Our derivation is based on a principle called the auxiliary function method,
similar to [1]. Let G(θ) denote an objective function to be minimized w.r.t. aparameter θ. A function G+ (θ, θ̂) which satisﬁes G(θ) = minθ̂ G+ (θ, θ̂) is thencalled an auxiliary function for G(θ), and θ̂ an auxiliary variable. The function
G(θ) is easily shown to be non-increasing through the following iterative update
rules: θ̂(s+1) ← argminθ̂ G+ (θ(s) , θ̂) and θ(s+1) ← argminθ G+ (θ, θ̂(s+1) ), whereθ̂(s+1) and θ(s+1) denote the updated values of θ̂ and θ after the s-th step.
An auxiliary function for standard NMF with β-divergence has been proposed [8]. This strategy can apply to our problem.
θ̂ denotes auxiliary variables{(λω,t,k )Ω×T ×K , (Zω,t )Ω×T } (∀k, λω,t,k ≥ 0, k λω,t,k = 1, Zω,t ∈ R) for convenience. We obtain the following auxiliary function:⎧ (β)⎪ ω,t − Yω,t Q(β−1)(β < 1)ω,t ⎨R(β) Yω,t(β−1)++J (θ, θ̂) =(5)(1 ≤ β ≤ 2) ,Q − Yω,t Qω,t
⎪ ω,tβ(β − 1)(β−1)ω,tω,t ⎩ (β)
Qω,t − Yω,t Rω,t(β > 2)(β)(β)βwhere Qω,t = (1/β) d,φ λω,t,d Pφ,t,d (Hω,φ,d Ud,t /λω,t,d )β and Rω,t = (1/β)Zω,t
β−1+Zω,t( d,φ Hω,φ,d Pφ,t,d Ud,t − Zω,t ). J + (θ, θ̂) is minimized w.r.t. θ̂ whenφ Hω,φ,d Pφ,t,d Ud,tλω,t,d = , Zω,t =Hω,φ,d Pφ,t,d Ud,t .
(6)
φ,d Hω,φ,d Pφ,t,d Ud,tφ,dMinimizing J + (θ, θ̂) w.r.t. P ∈ S is a search problem
P ← argmin J + (θ, θ̂)P∈S,(7)

which can be straightforwardly solved using the Viterbi algorithm. Diﬀerentiating J + (θ, θ̂) partially w.r.t. Hω,φ,t and Ud,t , and setting to zero, we obtain
update rules for Hω,φ,d and Ud,t :  β−2ϕ(β)
Max{2−β, 0}Pφ,t,d Ud,tt αω,t Yω,t γω,t,dHω,φ,d ←,(8) β−1 Min{1−β, 0}Pφ,t,d Ud,tt αω,t γω,t,d  β−2ϕ(β)ω αω,t Yω,t γω,t,dUd,t ← Ud,t,(9) β−1γω,t,dω αω,t 
where αω,t = d,φ Hω,φ,d Pφ,t,d Ud,t , γω,t,d = φ Hω,φ,d Pφ,t,d and ϕ(β) = 1/(2 −β) (β < 1), 1 (1 ≤ β ≤ 2), 1/(β − 1) (β > 2).

(TGSWGPE[=M*\?

(TGSWGPE[=M*\?

NMF with Markov-Chained Bases for Modeling Time-Varying Patterns










6KOG=U?










6KOG=U?

(a) Original spectrogram(TGSWGPE[=M*\?

153(b) MNMF-ltr, K = 1, Φ = 4








6KOG=U?

(c) Markov-chained bases, K = 1, Φ = 4(d) Activation, K = 1, Φ = 4Fig. 2. Original spectrogram of the extract of the piano (MIDI) (a), reconstructed
spectrograms (b), the bases (c) and the activation (d) learned by MNMF. MNMF
was able to decompose the evolution of the spectrum as a succession of several part,“attack”, “sustain”, “decay” and “release”.

2.4

Update Scheduling

The algorithm described above converges quickly. However, it often falls into
unexpected stationary points. H and U at the early stage of iterations are notuseful for estimating P. P ﬁxed using unrealistic H and U induces in turn H
and U in unexpected directions. We thus improve the update rule for P by(s)introducing an updating schedule. Here, let the scheduling parameter, kφ , at

(s+1)
(s)
(S)
(s+1)
s-th iteration satisfy ∀φ, kφ ≥ 0,
≥ k1 , k1 = 1, kφ≤
φ kφ = 1, k1
(s)

(S)

kφ , kφ= 0 (φ = 2, · · ·, Φ). We replace the update rule, Eq. (7) byΦ(s+1)
kn(s+1) P̂φ−(n−1),t,d ,Pφ,t,d ←(10)n=1where ∀n, P̂−(n−1),t,d = P̂Φ−(n−1),t,d and (P̂φ,t,d )Φ×T ×D = argminP∈S J + (θ, θ̂).
As a result, at each iteration the auxiliary function is not minimized anymore.
However, convergence is guaranteed.

3Simulation Results

In this section, some results on the application of our algorithm to audio signals.
All data were downmixed to mono and downsampled to 16kHz. The magnitude
spectrogram was computed using the short time Fourier transform with 32 ms
long Hanning window and with 16 ms overlap. The state transitions of the basesin MNMF were modeled using left-to-right (-ltr) and ergodic (-erg) models.
At ﬁrst, we tested whether the algorithm was able to learn in an unsupervised
way the time-varying spectral patterns from notes with a unique pitch. The
proposed method was applied to a piano note (C3) synthesized from MIDI, apiano note (C3) recorded from RWC-MDB-I-2001 No.1 [7] and a violin note

(TGSWGPE[=M*\?

M. Nakano et al.

(TGSWGPE[=M*\?

154




6KOG=U?





6KOG=U?

(c) MNMF-ltr, K = 1, Φ = 4(a) Original spectrogram(TGSWGPE[=M*\?





6KOG=U?

(d) Markov-chained bases, K = 1, Φ = 4(e) Activation, K = 1, Φ = 4

(TGSWGPE[=M*\?

(TGSWGPE[=M*\?

Fig. 3. Original spectrogram of the extract of the piano (RWC database) (a), reconstructed spectrograms (b) and (d), the bases (c) and the activation (d) learned by
MNMF. Time-varying spectral patterns are also learned.





6KOG=U?

(TGSWGPE[=M*\?

(TGSWGPE[=M*\?













6KOG=U?



5VCVG(d) MNMF-erg, K = 1, Φ = 4
5VCVG



6KOG=U?

(c) MNMF-ltr, K = 1, Φ = 4








(b) Standard NMF, K = 1(a) Original spectrogram
6KOG=U?










6KOG=U?

6KOG=U?

(e) State transition of MNMF-ltr(f) State transition of MNMF-erg

Fig. 4. Original spectrogram of the extract of the violin (RWC database) (a), reconstructed spectrograms (b), (c) and (d), and the activation (e) and (f) learned by
MNMF. The topology of the Markov chain aﬀects the state transitions.

(A ) recorded from RWC-MDB-I-2001 No.15. We used MNMF with β = 1 (thegeneralized Kullback-Leibler divergence). As shown in Fig. 2, 3 and 4, timevarying spectral patterns are learned in an unsupervised way.
Next, we applied our model to a mixture of vocal signals taken from RWCMDB-I-2001 No.45. The sequence is composed of 3 notes (D CFCA ): ﬁrst, eachnote is played alone in turn, then all combination of two notes are played andﬁnally all notes are played simultaneously. The result is shown in Fig. 5.

NMF with Markov-Chained Bases for Modeling Time-Varying Patterns

155(a) Mixed spectrogram(b) Partial spectrogram played by all notes. (c) Portion of estimated component: A(d) Portion of estimated component: F(e) Portion of estimated component: DFig. 5. Mixed spectrogram (a) and (b), and one portion of estimated spectrogram ofeach component (c), (d) and (e)Table 1. Source Separation Performance

Algorithms
NMF(D = 3)
NMF(D = 6)
NMF(D = 9)
NMF(D = 12)
MNMF-erg(Φ = 4)
MNMF-erg(Φ = 5)

RWC signal
Number of
H and U SNR(dB)
2919
3.55
5838
7.11
8757
10.63
11676
13.51
7536
9.62
9075
11.02

Real-world audio signal
Number of
H and U SNR(dB)
2871
6.49
5742
9.61
8613
10.50
11484
11.70
7488
9.57
9027
9.84

Finally, our model was applied to sound source separation. The tested signals are RWC signal (as shown in Fig. 5) and audio data recorded in realworld conditions from male vocal in the room size of 5.5m × 3.5m × 3mby IC recorder, whose sequence is composed similar to RWC signal (A, Dand E instead of D CF and A ). We use the signal-to-noise ratio (SNR) between each component and source as the measure for assigning components to(m)sources. The measure was calculated between the magnitude spectrograms Iω,t(n)and Iˆω,t of the mth reference and the nth separated component, respectively,

(m)
(m)
(n)
SNR = 10 log10 ( ω,t (Iω,t )2 / ω,t (Iω,t − Iˆω,t )2 ) . The SNR was averaged overall the sources to get the separation performance and each algorithm was run
10 times with random initializations. We set the number of bases to D = 3

156

M. Nakano et al.

for MNMF. This means that one component with Φ spectral patterns may beexpected to represent one source having the perceveid pitch of a note with vibrato. Standard NMF was used as the baseline. NMF with D = 3 was expected
to separate each source into a single component, while for NMF with D > 3
each source was expected to be split into the sum of several components. Acomponent n is assigned to the source m which leads to the highest SNR. As
reported by [2], using NMF with a large number of components and clustering
them using the original signals as references may produce unrealistically goodresults. One should thus keep in mind when comparing the results of our methodwith those of NMF that the experimental conditions were strongly biased in favor of NMF. As shown in Table 1, we can see that MNMF performs as well as
standard NMF when the total number of parameters H and H is similar (Pis excluded), although again for NMF the original signals need to be used tocluster the extracted components.

4Concluding Remarks

We developed a new framework for the sparse representation of audio signals.
The proposed model is extension of NMF, in which the bases consist of statetransition. We derived an eﬃcient algorithm for the optimization of the model
based on the auxiliary function method. Future work will includes the extensionof our model to automatic estimation of the number of bases and states.

Acknowledgements
This research was partially supported by CrestMuse Project under JST from Mext,
Japan, and Grant-in-Aid for Scientiﬁc Research (KAKENHI) (A) 20240017.

References
1. Lee, D.D., Seung, H.S.: Algorithms for non-negative matrix factorization. In Proc.
NIPS, vol. 13, pp. 556–562 (December 2001)
2. Virtanen, T.: Monaural sound source separation by nonnegative matrix factorization with temporal continuity and sparseness criteria. IEEE Transactions on Audio,Speech, and Language Processing 15, 1066–1074 (2007)
3. Hennequin, R., Badeau, R., David, B.: NMF with time-frequency activations tomodel non stationary audio events. In: Proc. ICASSP, pp. 445–448 (March 2010)
4. Ozerov, A., Fevotte, C., Charbit, M.: Factorial scaled hidden Markov model forpolyphonic audio representation and source separation. In Proc. WASPAA (2009)
5. Févotte, C., Bertin, N., Durrieu, J.-L.: Nonnegative matrix factorization with theitakura-saito divergence with application to music analysis. Neural Computation 21,793–830 (2009)
6. Févotte, C., Cemgil, A.T.: Nonnegative matrix factorizations as probabilistic inference in composite models. In: Proc. EUSIPCO, vol. 47, pp. 1913–1917 (2009)
7. Goto, M., Hashiguchi, H., Nishimura, T., Oka, R.: RWC music database: Popular,classical, and jazz music database. In: Proc. ISMIR, pp. 287–288 (2002)
8. Nakano, M., Kameoka, H., Le Roux, J., Kitano, Y., Ono, N., Sagayama, S.:
Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with β-divergence. In: Proc. MLSP (2010)