A MAJORIZATION-MINIMIZATION ALGORITHM WITH PROJECTED GRADIENTUPDATES FOR TIME-DOMAIN SPECTROGRAM FACTORIZATION
Hideaki Kagami† , Hirokazu Kameoka‡ , Masahiro Yukawa†† Dept. Electronics and Electrical Engineering, Keio University, Japan.
‡ NTT Communication Science Laboratories, NTT Corporation, Japan.
sequence of observed spectra can be approximated reasonably wellwhen each basis spectrum expresses the spectrum of an underlying
audio event that occurs frequently over the entire observed range.
Thus, with music signals, each basis spectrum usually becomes
the spectrum of a frequently used pitch in the music piece. In a
supervised/semi-supervised setting, NMF is first applied to train thebasis spectra of each sound source using the individually recorded
audio samples. At test time, NMF is applied to the spectrogram of atest mixture signal, where the subsets of the basis spectra are fixedat the pretrained spectra. In this way, the signals of the underlying
source components can be separated out using the Wiener filter
obtained with the estimated spectrograms of the individual sources.
Although the NMF approach has been shown to be successful,one limitation is that it assumes the additivity of magnitude (orpower) spectra, which holds only approximately, and does not takeaccount of phase information. To overcome this limitation, we havepreviously proposed a framework called complex NMF [7], where
the complex spectrum observed at each time frame is modeled as
the sum of components, each of which is described by the multiplication of a static basis spectrum, a time-varying amplitude anda time-varying phase spectrum. With a similar motivation, Parryet al. [8] and Févotte et al. [9] have independently proposed a
generative model of the complex spectrogram obtained with theshort-time Fourier transform (STFT) of a mixture signal, where
the power spectrogram of each component is modeled as a rank-1
matrix and the phase spectrogram is treated as uniformly distributed
latent variables. It can be shown that when each element of the complex spectrogram is assumed to independently follow a zero-mean
complex normal distribution, the maximum likelihood estimationof the model parameters amounts to fitting the NMF model to an
observed power spectrogram using the Itakura-Saito (IS) divergence
as a goodness-of-fit criterion. This approach is called IS-NMF. A
similar kind of generative model using a complex Cauchy distribution instead of a complex normal distribution has recently been
proposed [10].
Although these phase-aware NMF variants treat each elementof the phase spectrogram as an independent parameter (or latentvariable), in fact the phases of time-frequency components are constrained and dependent on each other. This is because the spectrograms obtained with typical time-frequency transforms (such as the
STFT and the wavelet transform) are redundant representations ofa signal. For example, the STFT spectrogram is computed by concatenating the Fourier transforms of overlapping short-time framesof the signal. Hence, all the elements of the STFT spectrogram mustsatisfy a certain condition to ensure that the waveforms within theoverlapping segment of consecutive frames are consistent [11]. The
shortcomings of the complex NMF and IS-NMF frameworks are thatthey fail to take account of this kind of redundancy. Moreover, thetime domain signal converted from the estimated complex spectrogram is in general not “optimal” unless the redundancy is taken into

ABSTRACT
We previously introduced a framework called time-domain spectrogram factorization (TSF), which realizes nonnegative matrix
factorization (NMF)-like source separation in the time domain.
This framework is particularly noteworthy in that, while maintaining the ability of NMF to obtain a parts-based representation ofmagnitude spectra, it allows us to (i) circumvent the commonly
made assumption with the NMF approach that the magnitude spectra of source components are additive and (ii) take account of theinterdependence of the phase/amplitude components at different
time-frequency points. In particular, the second factor has been
overlooked despite its potential importance. Our previous study
revealed that the conventional TSF algorithm was relatively slow
due to large matrix inversions, and the early stopping of the algorithm often resulted in poor separation accuracy. To overcome thisproblem, this paper presents an iterative TSF solver using projectedgradient updates. Simulation results show that the proposed TSF
approach yields higher source separation performance than NMF
and the other variants including the original TSF.
Index Terms— Audio source separation, non-negative matrix
factorization(NMF), projected gradient method
1. INTRODUCTION
Many sound recordings are mixtures of multiple sound sources. Audio source separation, i.e., the inverse operation of the mixing process, has long been a formidable challenge in the field of audio signal
processing [1].
For speech enhancement tasks, i.e., audio source separation
tasks where the mixture is given by speech and noise, deep neural network-based approaches have recently proved powerful [2].
A recently proposed deep-network-based approach [3] has further
made it possible to deal with “cocktail party” scenarios where theinterference is also speech. Although these methods have beenshown to work well when a large number of training samples areavailable, the supervised/semi-supervised non-negative matrix factorization (NMF) approach [4–6] still remains attractive for audio
source separation tasks particularly when only a limited number oftraining data are available. With the NMF approach, the magnitude(or power) spectrogram of a mixture signal, represented as a nonnegative matrix Y , is factorized into a product of two non-negative
matrices H and U . This can be interpreted as approximating theobserved spectra at each time frame as a linear sum of basis spectra
scaled by time-varying amplitudes, and amounts to approximating
the observed spectrogram as the sum of rank-1 spectrograms. The
This work was supported by the Support Center for Advanced
Telecommunications Technology Research (SCAT) and JSPS Grants-in-Aid(15K06081, 15K13986, 15H02757, 26730100).

978-1-5090-4117-6/17/$31.00 ©2017 IEEE

561

ICASSP 2017

importance of this weight parameter will be shown later. The firstterm of J (θ) becomes 0 when the magnitude spectrograms of all
the S have exactly rank-1 structures. It is important to note thatas with complex NMF, this model allows the components to cancel
each other out, and therefore some constraint is needed to induce
the sparsity of U . The second term of J (θ) is introduced for thispurpose, which we define as the ℓp norm, where λ > 0 weighs
the importance of the sparsity cost relative to the fitting cost. For
0 < p < 2, it promotes sparsity if the norm of U is bounded. To
bound U , we assume∑ 2
Hk,l = 1.
(3)

Fig. 1. Illustration of time-frequency basis functions.

kaccount, since the complex spectrogram of the converted signal is nolonger the same as the original complex spectrogram.
To exploit the intrinsically redundant structure of spectrograms explicitly, we previously introduced a framework called
time-domain spectrogram factorization (TSF) [12], which realizes
NMF-like source separation in the time domain. This framework
is particularly noteworthy in that, while maintaining the ability toobtain a parts-based representation of magnitude spectra, it directly
decomposes a time-domain signal into each underlying component.
Postprocessing to obtain a time domain signal from the estimatedspectrogram as typically required in the conventional frameworks
is implicitly involved in the optimization process. However, ourprevious study revealed that the conventional TSF algorithm was
relatively slow due to large matrix inversions, and the early stoppingof the algorithm often resulted in poor separation accuracy. This
paper proposes an iterative algorithm leveraging the projected gradient method to solve the TSF optimization problem in a tractablemanner.

2.2. Majorization-minimization algorithm
Let F (θ) be a cost function that we wish to minimize with respect
to θ. F + (θ, θ̄) is defined as an auxiliary function for F (θ) if itsatisfies
F (θ) = min F + (θ, θ̄).

(4)

θ̄

Here, we call θ̄ an auxiliary variable.
Theorem 1 ([13]). F (θ) is non-increasing under the following updates:θ̄ ← argmin F + (θ, θ̄),

(5)

θ̄

θ ← argmin F + (θ, θ̄).

(6)

θNote that the convergence of the algorithm is still guaranteed aslong as F + is decreased with respect to θ at each iteration. Thus,
F + does not need to be exactly minimized.

2. PROBLEM FORMULATION
2.1. Objective function

2.3. Auxiliary function of J

We denote an observed signal at time n by y[n], and the signal of theentire period by y = [y[1], . . . , y[N ]]T ∈ RN , where (·)T stands forthe vector (matrix) transpose. While the NMF approach decomposes
an observed magnitude spectrogram into the sum of rank-1 spectrograms, TSF decomposes y into the sum of L signal components:
y=

L
∑

sl ,

One difficulty as regards the current optimization problem comespfrom the nonsmoothness of |ψ H
k,m sl | and |Ul,m | for 0 < p < 1
in J . We can design an auxiliary function for J with a convenient
form
∑
1
2
ψ H sl − Hk,l Ul,m cl,k,m
J + (θ, θ̄) =βl,k,m k,m
l,k,m
∑[
]
2
+λ
p|Vl,m |p−2 Ul,m
+ (2 − p)|Vl,m |p , (7)

(1)

l=1

such that the magnitude spectrogram of each component is as close
to a rank-1 structure as possible. Here, we use ψ k,m ∈ CN todenote a basis function for the time-frequency transform, where k
and m are the frequency and time indices, respectively. By using
ψ k,m , the magnitude spectrogram of sl can be written as |ψ H
k,m sl |,where (·)H stands for the Hermitian transpose of a vector. Hence, we
formulate TSF as the optimization problem of minimizing
∑
∑
1
2
(|ψ H
|Ul,m |p ,J (θ):=
k,m sl |−Hk,l Ul,m ) +2λβl,k,ml,ml,k,m∑subject to y =
sl ,
(2)

l,m

by using the inequalities (A.30) and (A.31) given in the appendix,
where θ̄ = {c, V } with c = {cl,k,m } and V = {Vl,m }. A minimizer of the auxiliary function with respect to C and V has thefollowing closed-form expression:Hcl,k,m = ψ Hk,m sl /|ψ k,m sl |,(8)

Vl,m = Ul,m .

(9)3. OPTIMIZATION WITH PROJECTED GRADIENTlwith respect to θ = {H, U , S, β} where H = {Hk,l }, U ={Ul,m∑}, S = {sl } and β = {βl,k,m } with βl,k,m > 0 satisfying l βl,k,m = 1. This problem can be seen as a weighted leastsquares problem where βl,k,m is the reciprocal of the weight. The

By using the method of Lagrange multipliers, we obtain the update
rule for sl as
sl = Ψ−1l (dl − µ),

562

(10)

∑

where
Ψl =

∑ 2Re[ψ k,m ψ Hk,m ]βl,k,m

k,m

dl =

,

∑ 2Re[cl,k,m ψ k,m ]Hk,l Ul,m,βl,k,m)( ∑ )−1 ( ∑
µ=
Ψl
Ψ−1
l dl − y .

{}V := s̃ ∈ RLN | Gs̃ = y .

(11)

(20)

Then, the projection onto the affine subspace V is given by [14]PV (s̃) = s̃ − GT (GGT )−1 (Gs̃ − y).

(12)

k,m

l

l sl can then be written equivalently as y = Gs̃. Also define

(21)

T

By considering GG = LI N and
I N I N ... I N
 I N I N ... I N T ∈ RLN ×LN ,
G G=
..




.

(13)

l

The conventional TSF algorithm uses (10) as the update rule for sl .
As can be seen from the above, one problem with the conventionalalgorithm is that it involves large matrix inversions, resulting in avery slow algorithm. In [12], we showed that in a particular casewhere β1,k,m = · · · = βL,k,m (i.e., ∀l, βl,k,m = 1/L), Ψl becomes an identity matrix, which means that the update rule of sl canbe computed without matrix inversions. However, it transpired that
this β setting results in a poor separation accuracy. We show in thissection that using the projected gradient method to update sl allows
us to sidestep these computations and leads to a reasonably efficient
optimization algorithm.

(22)

I N I N ... I N
we obtain
1[PV (s̃)]l = sl −
L

( L
∑

)
sl′ − y

.

(23)

l′ =1

Accordingly, we seek to minimize J + with respect to sl by alternately performing (15) and (23). In fact, ρ gives the upper boundof κ, and this ensures γ ∈ (0, 2/κ) so that the projected gradientalgorithm converge to a minimizer of J + under the linear constraint[15].

3.1. Update rule for sl
3.2. Summary of proposed algorithm

In this subsection, we solely consider the update for sl and fix all
the other variables included in θ and θ̄. We therefore regard J + asa function of sl , and denote it by J + (sl ). The partial derivative ofJ + with respect to sl is given as∇sl J + (sl ) = Ψl sl − dl .

By setting the partial derivative of J + with respect to each elementof U at zero, the update rule for U is obtained as∑
Hk,l |ψ H
k,m sl |/βl,k,m
Ul,m = ∑ k2
.
(24)
p−2
k Hk,l /βl,k,m + λp|Vl,m |(14)Each term of (14) can be computed efficiently. Indeed, when ψ k,mis defined as the STFT basis function, ψ H
the (k, m)
k,m sl represents∑
element of the STFT of sl . Since the operator k,m Re[ψ k,m ·]in (11) corresponds to the inverse STFT process, Ψl sl can beobtained by computing the inverse STFT of a vector containing2ψ H
k,m sl /βl,k,m . Analogously, dl in (12) can be obtained by computing the inverse STFT of a vector containing 2Hk,l Ul,m cl,k,m/βl,k,m . We can thus efficiently update sl using the gradient-based
update rulesl ← sl − γ∇sl J + (sl ),Similarly, we obtain the update rules for H and β by using themethod of Lagrange multipliers as follows:∑Hm Ul,m |ψ k,m sl |/βl,k,m,(25)
Hk,l = √∑ ∑H
2m Ul,m |ψ k,m sl |/βl,k,m )k(
βl,k,m = ∑

l

(15)

|ψ H
k,m sl | − Hk,l Ul,m

.

(26)

Overall, the proposed iterative algorithm is summarized as follows.

where γ ∈ (0, 2/κ) is step size. Here, κ > 0 is the Lipschitz constant of the gradient operator ∇sl J + , and is given by the largest
eigenvalue of the matrix Ψl . For the sake of computational efficiency, we use the following step size:
2(16)
γ := ∈ (0, 2/κ),ρ∑∑ ∑ 2{w(n − αm )}2
ρ :=tr (Ψl ) =,(17)βl,k,ml|ψ Hk,m sl | − Hk,l Ul,m

1. Initialize H, U and S.
2. Update c and V using (8), (9).
3. Update S and γ using (14), (16)–(17), (23).
for t = 1 : iteration
(
)
(t+1)
(t)
(t)
sl
= PV sl − γ∇sl J + (sl )end
4. Update U , H and β using (24)–(26).

l,n k,m

4. RELATION TO COMPLEX NMF

where w(n) is the window function and αm is the hop size of the
STFT.
Define


s1
 . 
LN

s̃ := (18)
 ..  ∈ R ,

It is interesting to note that the auxiliary function J + has a similar
form to that of complex NMF [7]. The aim of complex NMF isto approximate an observed complex spectrogram Yk,m using thefollowing model:∑
Yk,m ≈
Hk,l Ul,m ejϕl,k,m ,(27)sL
G := [I N · · · I N ] ∈ RN ×LN ,

l

(19)

where ϕl,k,m denotes the phase spectrogram of the l-th signal component. In [7], the cost function that must be minimized with respect

where I N is the N × N identity matrix. The linear constraint y =563

test segments, and take the average SNR improvement obtained withthe three repeated rounds. With all these methods, the basis spectra
were pretrained using the individual tracks of the training data, andthen source separation was performed on the test data. All the audio
samples were monaural and sampled at 22.05kHz. The STFT was
computed using a square-root Hanning window that was 32 ms longwith a 16 ms overlap. With both methods, 6 basis spectra wereassigned to each track. Thus, for 5-track recordings, a total of 30
basis spectra were used for the separation. Fig. 2 shows the SNR
improvements after the separations with the six methods. These
results show that proposed TSF performed better than NMF, CNMF
and original TSF.

Fig. 2. SNR improvements.

to Hk,l , Ul,m and ϕl,k,m is defined as∑∑∑
2|Ul,m |p ,Hk,l Ul,m ejϕl,k,m + 2λI(θ) =
Yk,m −k,ml

6. CONCLUSION
This paper presented an efficient TSF algorithm that performs source
separation in the time domain. A noteworthy advantage of TSF is
that, while preserving the desirable property of NMF, it exploits(i) the additive nature of time-domain signals and (ii) the interdependence of the phase/amplitude components at different timefrequency points. The proposed iterative algorithm was built basedon a majorization-minimization scheme leveraging the projectedgradient method. Simulation results showed that the proposed TSF
algorithm yielded better source separation performance than NMF.

l,m

(28)
with θ = {H, U , ϕ} where H = {Hk,l }, U = {Ul,m } and ϕ ={ϕl,k,m }. It can be shown that∑

1
2
Xl,k,m − Hk,l Ul,m ejϕl,k,m
βl,k,m
l,k,m
∑{
}
2
+λ
p|Vl,m |p−2 Ul,m
+ (2 − p)|Vl,m |p ,I + (θ, θ̄) =(29)
Appendix A. INEQUALITIESl,mis an auxiliary function for I, where θ̄ = {β, X}, β = {βl,k,m }and X = {Xl,m } are the auxiliary∑ variables. Here, βl,k,m can be
any positive number∑ satisfying l βl,k,m = 1, and Xl,k,m mustsatisfy Yk,m = l Xl,k,m . Xl,k,m can be viewed as an estimate of
the complex spectrogram of the l-th signal component.
By comparing (7) and (29), we see that Xl,k,m and ejϕl,k,m in
the auxiliary function of the complex NMF objective are analogous
to sl and cl,k,m in that of the present TSF objective, respectively.
One drawback with complex NMF is that Xl,k,m is not guaranteed tosatisfy the explicit condition that complex spectrograms must satisfy.
This implies that the complex NMF algorithm is designed to searchfor optimal parameters in an unnecessarily large solution space. By
contrast, the present algorithm always ensures that the estimate of
the complex spectrogram of the l-th latent component ψ Hk,m sl is
associated with a time-domain signal sl , keeping the search within
an appropriate solution space.

We summarize the inequalities that are used to design the auxiliary
function.
Lemma 1. For any complex number z and any complex number csatisfying |c| = 1, it holds that
− |z| ≤ −Re[c∗ z].

(A.30)

Equality holds when c = z/ |z|.
Lemma 2. When 0 < p < 2, for any real or complex number x, it
holds that
2 |x|2 ≤ p |v|p−2 |x|2 + 2 − p |v|p .

(A.31)

Equality holds when v = x.
7. REFERENCES[1] S. Roweis, “One microphone source separation,” in NIPS,
2000, vol. 13, pp. 793–799.

5. EXPERIMENTAL RESULTS[2] Y. Wang and D. Wang, “Towards scaling up classificationbased speech separation,” IEEE Trans. Audio, Speech, Language Process., vol. 21, no. 7, pp. 1381–1390, 2013.

We quantitatively compared the source separation performance ofthe proposed time-domain spectrogram factorization (TSF) by conducting supervised source separation experiments, original TSF[12], complex NMF (CNMF) [7], and NMF using the Euclidean
distance (EUC), Kullback-Leibler (KL) divergence and ItakuraSaito (IS) divergence by conducting supervised source separation
experiments. The original TSF minimizes (2) in a particular casewhere β1,k,m = · · · = βL,k,m (i.e., ∀l, βl,k,m = 1/L). We usedprofessionally produced music recordings from the SiSEC 2013
database, available at https://sisec.wiki.irisa.fr/, as the experimental
data. Each recording was a mixture of multiple tracks, each of which
was produced by a single instrument or singer. The tracks were alsoavailable separately. We performed 3-fold cross validation. We
partitioned each recording into three segments, used one segment
as the test data and the other two segments as the training data, repeated signal-to-noise (SNR) evaluations three times with different[3] J. R. Hershey, Z. Chen, J. Le Roux, and S. Watanabe, “Deepclustering: Discriminative embeddings for segmentation andseparation,” in Proc. ICASSP 2016, 2016, pp. 31–35.
[4] P. Smaragdis, B. Raj, and M. Shashanka,“Supervisedand semi-supervised separation of sounds from single-channel
mixtures,” in Proc. LVA/ICA 2010, 2010, pp. 140–148.
[5] Y. Morikawa and M. Yukawa, “A sparse optimization approachto supervised NMF based on convex analytic method,” inProc. ICASSP, May 2013, pp. 6078–6082.
[6] H. Kagami and M. Yukawa, “Supervised nonnegative matrix
factorization with Dual-Itakura-Saito and Kullback-Leibler divergences for music transcription,” in Proc. EUSIPCO, Hungary, Aug. 2016, pp. 1138–1142.

564

[7] H. Kameoka, N. Ono, K. Kashino, and S. Sagayama, “Complex NMF: A new sparse representation for acoustic signals,”in Proc. IEEE ICASSP, April 2009, pp. 3437–3440.
[8] R.M. Parry and I. Essa, “Phase-aware non-negative spectrogram factorization,” Independent Component Analysis andSignal Separation, pp. 536–543, 2007.
[9] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix
factorization with the itakura-saito divergence: With application to music analysis,” Neural computation, vol. 21, no. 3, pp.
793–830, 2009.
[10] A. Liutkus, D. Fitzgerald, and R. Badeau, “Cauchy nonnegative matrix factorization,” in Proc. WASPAA 2015, 2015.
[11] J. Le Roux, H. Kameoka, N. Ono, and S. Sagayama, “Fast signal reconstruction from magnitude STFT spectrogram basedon spectrogram consistency,” in Proc. DAFx-10, 2010, pp.
397–403.
[12] H. Kameoka, “Multi-resolution signal decomposition with
time-domain spectrogram factorization,” in Proc. ICASSP
2015, 2015, pp. 86–90.
[13] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in Advances in Neural Information Processing Systems 13, pp. 556–562. 2001.
[14] D. G. Luenberger, Optimization by Vector Space Methods,
New York: Wiley, 1969.
[15] A. A. Goldstein, “Convex programming in Hilbert space,”
Bull. Amer. Math. Soc. 70, pp. 709–710, 1964.

565