IEICE TRANS. FUNDAMENTALS, VOL.E105–A, NO.6 JUNE 2022

906

PAPER

Supervised Audio Source Separation Based on Nonnegative Matrix
Factorization with Cosine Similarity Penalty
Yuta IWASE† , Nonmember and Daichi KITAMURA†a) , Member

SUMMARY In this study, we aim to improve the performance of audio
source separation for monaural mixture signals. For monaural audio source
separation, semisupervised nonnegative matrix factorization (SNMF) canachieve higher separation performance by employing small supervised signals. In particular, penalized SNMF (PSNMF) with orthogonality penalty
is an effective method. PSNMF forces two basis matrices for target andnontarget sources to be orthogonal to each other and improves the separation accuracy. However, the conventional orthogonality penalty is basedon an inner product and does not affect the estimation of the basis matrix
properly because of the scale indeterminacy between the basis and activation matrices in NMF. To cope with this problem, a new PSNMF withcosine similarity between the basis matrices is proposed. The experimental
comparison shows the efficacy of the proposed cosine similarity penalty insupervised audio source separation.
key words: audio source separation, nonnegative matrix factorization,orthogonality, cosine similarity

1.

Fig. 1 Matrix decomposition by NMF. Amplitude spectrogram of audio
signal is input as nonnegative matrix.

Introduction

Audio source separation is a technique of separating or extracting individual audio sources from the observed mixture signal. In particular, nonnegative matrix factorization(NMF) [1], [2] has been utilized in many situations [3]–[7],depending on the recording conditions and applications. As
shown in Fig. 1, NMF is an algorithm that decomposes anobserved nonnegative matrix X into two nonnegative matrices, a basis matrix F and an activation matrix Q. The basis matrix includes frequently appearing spectral patterns in
the observed matrix as basis vectors, and their time-varying
gains are included in the activation matrix. Audio source
separation can be achieved by clustering these estimatedcomponents into each source.
In recent years, deep learning methods have attracted
wide attention and achieved high performance even in thesource separation field. However, when the training dataof target sources are excessively limited, the performance ofdeep-learning-based approaches degrades markedly. In sucha situation, supervised or semisupervised NMF [8], [9] isstill a promising algorithm, in which spectral patterns of thetarget sources are pretrained using only one sequence of theirsample sounds, e.g., octave notes of the target instruments.
In this paper, we only focus on a problem of music source
separation based on SNMF, where only few sample notes of
the target instruments are available as a supervision.
Manuscript received November 12, 2021.
Manuscript publicized December 8, 2021.
† The authors are with the National Institute of Technology,
Kagawa College, Takamatsu-shi, 761-8058 Japan.
a) E-mail: kitamura-d@t.kagawa-nct.ac.jp
DOI: 10.1587/transfun.2021EAP1149

Fig. 2 SNMF algorithm. In training stage, sample sound of target source
is decomposed by simple NMF, and supervised basis matrix F is obtained,
which represents spectral patterns of target source. In separation stage,
mixture sound is decomposed while F is fixed.

Figure 2 shows the SNMF-based source separation algorithm that consists of training and separation stages. In
the training stage, an amplitude spectrogram of the sample
sound of the target source, X, is decomposed by simple
NMF, and the supervised basis matrix F is obtained. In theseparation stage, an amplitude spectrogram of the observed
mixture Y is decomposed into the target spectrogram F G
and the other (nontarget) spectrogram HU using the fixed
supervised basis matrix F.
When the target and nontarget sources in the observed
mixture contain similar spectra, the separation performanceof SNMF degrades. This is due to the fact that similar spectra can be represented by either the supervised basis matrix
F or the nontarget basis matrix H. In this case, either of thefollowing problems occur: (a) some components of the target
source are incorrectly included in HU, or (b) some components of the nontarget sources are inappropriately captured
by F G. Since SNMF utilizes the sample sound of only the

Copyright © 2022 The Institute of Electronics, Information and Communication Engineers

IWASE and KITAMURA: SUPERVISED AUDIO SOURCE SEPARATION BASED ON NMF WITH COSINE SIMILARITY PENALTY

907

target source, there exists a trade-off between the problems(a) and (b).
Penalized SNMF (PSNMF) [9] is a technique that addresses the above-mentioned trade-off problem. PSNMF
forces the nontarget basis matrix H to be orthogonal to the supervised basis matrix F by adding an orthogonality penalty
in the optimization of the separation stage. This is direct
solution for the problem (a). We can avoid lack of the target
source components in the estimated matrix F G, although
the problem (b) might be encouraged. Since semisupervised
techniques aim to accurately extract the target source froma mixture, ensuring sound quality of the estimated target
source is an important objective for many applications, e.g.,music editing software. In addition, by tuning intensity ofthe orthogonality penalty, PSNMF can control the trade-off
between the problems (a) and (b), which enables us to apply
audio source separation to various objectives.
2,
The penalty term proposed in [9] is defined as kF T H kFr
a sum of the inner products between all the bases in F andH (k · kFr is the Frobenius norm). Since the spectral patternsin H become dissimilar from those in F, the separation performance is greatly improved. However, the formulation of
conventional PSNMF is incorrect because the orthogonality
2 does not affect the estimation of H proppenalty kF T H kFr
erly owing to the scale indeterminacy between the basis andactivation matrices in NMF.
In this paper, to solve the above-mentioned problem,
two new orthogonality penalties based on cosine similarity are introduced to PSNMF. Since cosine similarity does
not depend on scales of vectors, the proposed methods canappropriately force the nontarget spectral bases in H to be orthogonal to those in F even if the scale indeterminacy exists
between H and U. A convergence-guaranteed multiplicative
update rule of the proposed methods is derived on the basis of the majorization-minimization (MM) algorithm [10].
The validity of the proposed methods is confirmed by musicsource separation experiments.
2.

Conventional Methods

2.1

NMF and SNMF

The optimization problem in NMF [1], [2] is formulated as
min D(X kF Q) s.t. fik , qk j 0 ≥ 0 ∀i, j 0, k,
F,Q

(1)

NMF-based audio source separation [9], [11]. The KL divergence is defined as

Õ
mi j
mi j log
D(M k N ) =− mi j + ni j ,
(2)
ni j
i, j
where mi j and ni j are the elements of matrices M and N ,respectively, and j = 1, 2, . . . , J. Thus, F and Q can be
estimated by solving the minimization problem of Eq. (1).
In SNMF [8], the supervised basis matrix F is pretrained by applying simple NMF to the sample signal X ofthe target audio source in the training stage. The basis matrix
F contains the frequently appearing spectral patterns of thetarget source as K column vectors (bases). In the separation stage, the amplitude spectrogram of the mixture signal
Y ∈ RI×J≥0 is decomposed using the supervised basis matrix
F as follows:
min D(Y kF G + HU)

G,H ,Us.t. gk j , hil , ul j ≥ 0 ∀i, j, k, l,(3)

where G ∈ RK×J
≥0 is the activation matrix for F. More×L and U ∈ R L×J are the basis and activation
over, H ∈ RI≥0≥0matrices for representing nontarget sources, and gk j , hil ,and ul j are the elements of G, H, and U, respectively, andl = 1, 2, . . . , L represents the index of bases in H. Ideally,
the components of the target source in Y are extracted as
F G, and the components of the other nontarget sources aremodeled by HU, resulting in the separation of target andnotarget sources. However, when the target and nontarget
sources contain similar spectra, the components can be represented by either F or H. This ambiguity will cause thefollowing problem: part of the target source is captured by
HU or part of the nontarget sources is mixed into F G, which
degrades the accuracy of source separation.
2.2

Conventional PSNMF Based on Inner ProductTo solve the above-mentioned problem, PSNMF [9] wasproposed. In PSNMF, the nontarget bases in H are forced
to be as orthogonal as possible to the supervised bases in F
by imposing an orthogonality penalty to the cost function of
the separation stage as follows:
min D(Y kF G + HU) + µPinner (F, H)

G,H ,U
0where X ∈ RI×J
≥0 is a nonnegative observed matrix and is anamplitude spectrogram in this paper. In addition, F ∈ RI×K
≥00
and Q ∈ RK×J
are the basis and activation matrices, fik≥0
and qk j 0 are the elements of F and Q, and i = 1, 2, . . . , I,
j 0 = 1, 2, . . . , J 0, and k = 1, 2, . . . , K represent the indices offrequency bins, time frames, and basis vectors, respectively.
Moreover, D(M kN ) is a divergence function between two
I×J
input matrices M ∈ RI×J
≥0 and N ∈ R ≥0 . In this paper, weonly consider the generalized Kullback–Leibler (KL) divergence because it is experimentally confirmed that the KLdivergence-based NMF provides the highest performance fors.t. gk j , hil , ul j ≥ 0 ∀i, j, k, l,(4)

where µ > 0 is the weight coefficient and Pinner (F, H) isdefined as
2
Pinner (F, H) = kF T H kFr

!2
=

Õ Õ
k,l

fik hil.

(5)i

The penalty term Pinner (F, H) corresponds to the sum ofsquared inner products between two bases in F and H.

IEICE TRANS. FUNDAMENTALS, VOL.E105–A, NO.6 JUNE 2022

908

Therefore, the nontarget basis matrix H is estimated by taking the following perspectives into account: (i) the divergence between Y and F G + HU should be minimized and(ii) the bases in H should be as orthogonal as possible to thebases in F. As a result, the source separation performanceof PSNMF is greatly improved from that of SNMF [9].
The update rules for the variables G, H, and U basedon the MM algorithm [10] are respectively derived as follows [9]:
Í
yi j
Í
Í
i k 0 fi k 0 gk 0 j + l 0 hil 0 ul 0 j fikÍ
gk j ← gk j ·
,
(6)
i fik
Í
yi j
Í
Í
j k 0 fi k 0 gk 0 j + l 0 hil 0 ul 0 j ul jÍhil ← hil · Í,(7)
j ul j + µ fik i0 fi0 k hi0 l
Í
yi j
Í
Í
i k 0 fi k 0 gk 0 j + l 0 hil 0 ul 0 j hilÍ,(8)
ul j ← ul j ·i hilwhere yi j is the element of Y .
3.

Proposed Method

3.1

Motivation and Strategy

In conventional PSNMF Eq. (4), the penalty term
Pinner (F, H) is minimized to increase the dissimilarity between F and H. However, in practice, the penalty term
Pinner (F, H) can be minimized by simply multiplying H by
α > 0 as2
lim Pinner (F, αH) = lim αkF T H kFr

α→+0

(9)

In this case, the activation matrix U can be updated as U ←α−1 U so that the value of divergence does not change, i.e.,D(Y kF G + αH[α U]) = D(Y kF G + HU).
−1

3.2

Proposed Method 1: Log-Cos PSNMF

3.2.1

Cost Function

In Log-Cos PSNMF, we consider the following optimization
problem:

α→+0

= 0.

L1 norm of the activation matrix is imposed to the cost function of simple NMF. For the sparse NMF, in [13], a norm
constraint of bases is newly imposed to the NMF optimization to avoid the scale indeterminacy problem. Although
PSNMF may be solved by introducing the same constraint
as in [13], the theoretical convergence in this approach cannot be ensured. This is because the method [13] utilizes
fraction-based update rules with positive and negative terms
in the gradient. This heuristic update rule has often appearedin the history of NMF-based methods (e.g., [14], [15]), but
its theoretical convergence with the norm constraint has notbeen proven.
In this paper, to solve the problem in conventional
PSNMF, we propose to utilize cosine similarity as an orthogonality penalty term. Since cosine similarity depends
on only the angle between two input vectors, we can measure
the orthogonality regardless of the lengths of the vectors andthe scale indeterminacy in NMF does not affect the penalty
term. We propose two types of cosine similarity penalty:
logarithmic cosine similarity (hereafter referred to as LogCos PSNMF) and simple cosine similarity (hereafter referred
to as Cos PSNMF) penalties. In addition, we derive MMalgorithm-based (convergence-guaranteed) update rules forboth Log-Cos PSNMF and Cos PSNMF and obtain a tuningfree and easy-to-use optimization algorithm.

(10)

For this reason, the penalty term Pinner (F, H) does not affectthe optimization of H, and the minimization problem Eq. (4)is equivalent to Eq. (3). Therefore, the orthogonalization
between the bases in F and H is not properly performed inconventional PSNMF Eq. (4).
To cope with this problem, in [9], PSNMF is implemented with a basis normalization process hl ← hl /khl kand ul ← khl kul after every update of H so that the matrix HU does not change, where hl = [h1l , h2l , · · · , hIl ]T ,ul = [ul1, ul2, · · · , ulJ ]T , and k · k is an arbitrary norm. Although the source separation performance is improved bythis heuristic implementation technique, this iterative normalization process is not a fundamental (or mathematical)
solution. In addition, this basis normalization can increase
the value of cost function Eq. (4), hence theoretical convergence of the MM algorithm (the update rules Eqs. (6)–(8))
is lost.
A similar problem was pointed out in the context ofsparse NMF [12], where sparse regularization based on the

min J1 s.t. gk j , hil , ul j ≥ 0 ∀i, j, k, l,

G,H ,U(11)

where
J1 = D(Y kF G + HU) + µPlogcos (F, H),(12)ÍÕi fik hil
Plogcos (F, H) =
log 
(13)
 12 Í
 12
Í
2
2
k,l
i fik
i hil
"
Õ
Õ
Õ
1
=
log
fik hil − log
fik2
2
i
i
k,l
#
Õ
1
hil2 .
(14)
− log
2
i
The penalty term Plogcos (F, H) corresponds to the sum of
logarithmic cosine similarities for all combinations of supervised bases in F and nontarget bases in H. The logarithm
function in Eq. (13) is employed to decompose the cosine
similarity as Eq. (14), resulting in a simpler form of theMM-algorithm-based update rules compared with those of
Cos PSNMF.

IWASE and KITAMURA: SUPERVISED AUDIO SOURCE SEPARATION BASED ON NMF WITH COSINE SIMILARITY PENALTY909

and only if
3.2.2

Derivation of MM-Algorithm-Based Update Rules

Since it is difficult to directly minimize the cost function
Eq. (12), we use the MM algorithm [10] as in conventional
PSNMF. Note that since the penalty term Plogcos (F, H) in
Eq. (12) does not depend on the activation matrices G andU, the update rules for G and U that minimize Eq. (12) areequivalent to Eqs. (6) and (8), respectively.
The cost function Eq. (12) can be rewritten as follows:
"
!
Õ
Õ
Õ
yi j log yi j − yi j log
J1 =fik gk j +hil ul j
i, j

k

l

#
Õ

− yi j +fik gk j +Õkhil ul j

l

"
+µ

Õ

log

Õ

fik hil −

i

k,l

Õ
1
log
fik2
2
i

#
Õ
1
2
− log
hil .
2
i

(15)

The second, sixth, and eighth terms of Eq. (15) containthe sum of the variables (gk j , hil , and ul j ) in a logarithm
function, and it is difficult to calculate the stationary point
w.r.t. the variables. In the MM algorithm, we design an upperbound function to indirectly optimize the cost function whose
stationary points are difficult to calculate. Since the second
term, i.e., the negative logarithmic function (− log(·)), is aconvex function, by applying Jensen’s inequality, we candesign the upper bound function as follows:!
ÕÕ
hil ul j
fik gk j +
− log
l

k

= − logÕ αi jk fik gk j
αi jk

k

+

Õ βi jl hil ul j
βi jl

l

fik gk j≤−−
αi jk log
αi jk
k
Õ

!

Õ
l

hil2 = − log

Õ

i

γil

i

≤−

Õ

γil log

i

= −2

Õ
i

The sixth term of Eq. (15) is a positive logarithm(+ log(·)) function that includes the sum of the variables(hil ). Since the positive logarithmic function is a concave
function, by applying the tangent-line inequality, the upperbound function can be designed as
!
Õ
1 Õ
log
fik gk j ≤fik gk j − δkli + log δkli ,δkli ii(21)
where δkli > 0 is an auxiliary variable. The equality in
Eq. (21) holds if and only ifÕδkli =
fik hil .
(22)
iFrom Eqs. (16), (17), and (21), the upper bound function
J1+ of J1 can be designed as
J1 ≤ J1+Õ=yi j log yi ji, jfik gk j Õ
hil ul j
βi jl log
αi jk log+
− yi j
αi jk
βi jl
l
k
#
Õ
Õ
− yi j +fik gk j +hil ul j

!

Õ

l

!
1
fik hil − δkli + log δkli+µ
δkli i
k,l
Õ
Õ
1
− log
fik2 −
γil log hil
2
i
i
#
1Õ
+
γil log γil .
2 i
Õ

(23)

The update rule for hil can be derived by solving∂J1+ /∂hil = 0 w.r.t. hil as

γil
γil

(19)
(20)

k

hil2
hil2

h2
γil = Í il 2 .
i 0 hi 0 lÕwhere Íαi jk > 0Íand βi jl > 0 are auxiliary variables thatsatisfy k αi jk + l βi jl = 1. Similarly to Eq. (16), we design
the upper bound function of the eighth term of Eq. (15) as
Õ

(18)

ik 0

"

hil ul j
βi jl log
,
βi jl
(16)

− log

fik gk jÍ,
k 0 f gk 0 j + l hil ul jhil ul jÍ
βi jl = Í,k fik gk j + l0 hil0 ul0 j

αi jk = Í,γil log hil +ÍÕγil log γil , (17)iÍwhere γil > 0 is an auxiliary variable that satisfies i γil =1. The equality in Eqs. (16) and (17) respectively holds if

j yi j βi jl + µ

hil = Íj ul j + µ

Í

k γil
.
fi k
k δk l iÍ(24)

By substituting the equality conditions Eqs. (19), (20), and(22) for Eq. (24), we obtain the multiplicative update rule forhil asIEICE TRANS. FUNDAMENTALS, VOL.E105–A, NO.6 JUNE 2022

910

Í

j

Í

hil ← hil ·yi jÍ

k 0 fi k 0 gk 0 j +Í

l 0 hil 0 ul 0 j

j ul j + µÍkul j + µK Í h0 ilh2
i

fi k
Í

i0 l

.

i 0 fi 0 k hi 0 l(25)Note that Eq. (25) coincides with Eq. (7) when µ = 0.
By iterating the update rule Eq. (25), we can estimate
H that tends to be orthogonal to F. However, in Log-Cos
PSNMF, the penalty term Plogcos (F, H) becomes −∞ whenhl is perfectly orthogonal to any bases in F or hl = 0 for any
l. In such a case, the update rule Eq. (25) is undefined, andthe solution of minimization problem Eq. (11), J1 → −∞, ismeaningless. To avoid this inherent problem, the following
flooring process [16] is applied in each iteration:hil ← max(hil , ε),(26)

where ε is the machine epsilon and max(·) is a function that
returns the largest value of the input arguments.
3.3

Cost Function

Log-Cos PSNMF utilizes a logarithm function in the penalty
Plogcos (F, H) to decompose the cosine similarity function.
This logarithm function makes the cost function J1 unbounded, namely, the minimum value of J1 becomes −∞.
We apply the flooring Eq. (26) to avoid the meaningless solution, but such heuristic treatment is not essential as well as thebasis normalization described in Sect. 3.1. To cope with thisproblem, we propose Cos PSNMF, which utilizes a simple
cosine similarity penalty. Although derivation and update
rules become more complicated compared with those in LogCos PSNMF, we can still obtain a convergence-guaranteed
optimization algorithm in this formulation.
We consider the following optimization problem:
min J2 s.t. gk j , hil , ul j ≥ 0 ∀i, j, k, l,

G,H ,U(27)

where
J2 = D(Y kF G + HU) + µPcos (F, H),ÍÕi fik hilPcos (F, H) =Í
 21 Í 12 .
2
2
k,l
i fik
i hil

i

where il > 0 is an auxiliary variable that satisfies
The equality in Eq. (30) holds if and only ifi il = 1.

Íh2
il = Í il 2 .
i0 hi0 l(31)

From Eqs. (16) and (30), the upper bound function J2+
of J2 can be designed as

Proposed Method 2: Cos PSNMF

3.3.1

this function hinders the direct calculation of a stationary
point w.r.t. hil . Since this function is convex, we can design
the upper bound function by applying Jensen’s inequality as
follows:! − 12! 1
ÕÕ h2 − 2
il
2
=
il
hil
il
i
i
!− 1
Õ
hil2 2
≤
il
il
i
Õ 3
=
il2 hil−1,
(30)

(28)
(29)

In Cos PSNMF, the logarithm function in Plogcos (F, H) isomitted, and the minimum value of the cost function Eq. (28)is bounded by zero.

J2 ≤ J2+Õ=yi j log yi ji, jhil ul j
fik gk j Õ
βi jl log+
− yi jαi jk log
αi jk
βi jl
l
k
#
Õ
Õ
− yi j +fik gk j +hil ul j

!

Õ

k

l

!
! − 21ÕÕ  Õ
3 h 0

i
l
2
2
.

fi0 k i00l
+µ
fi0 k

hi00l i0 ,i00k,l  i0

(32)

The update rule for hil can be derived by solving∂J2+ /∂hil = 0 w.r.t. hil . The derivative can be obtainedas! − 12Õ  −yi j βi jlÕ  Õ
2
+ ul j + µ
fi0 k

hil
j
k  i0!#Õ 3 1
1 23 Õ
− 2 ilfi0 k hi0l + fik
i20l
= 0.
(33)
hi0l
hil i0 ,i0i ,iRearrangement of Eq. (33) gives the quadratic equation
ail hil2 + bil hil + cil = 0,

(34)

where
3.3.2

Derivation of MM-Algorithm-Based Update Rules

Similarly to Log-Cos PSNMF, direct minimization of thecost function Eq. (27) is difficult, and we apply the MM
algorithm. The update rules of G and U are equivalent toEqs. (6) and (8), respectively.
ÍThe penalty term Pcos (F, H) includes ( i hil2 )−1/2 , and! − 21ail =Õ

ul j + µ

j

Õ
k

fik

Õ
i0

fi20 k

!
1

,
hi0l
i0 ,i

Õ

3
2
i0 l

(35)
bil = −

Õj

yi j βi jl ,(36)

IWASE and KITAMURA: SUPERVISED AUDIO SOURCE SEPARATION BASED ON NMF WITH COSINE SIMILARITY PENALTY

911

3
2

cil = −µil! − 12Õ Õ!
Õfi20 ki0k

fi0 k hi0l .

(37)

i0 ,i

Thus, the stationary point of J2+ w.r.t. hil is given by thequadratic formula:q−bil ± b2il − 4ail cil
hil =
,
(38)
2ail
where ± in Eq. (38) should be determined so that the r.h.s. of
Eq. (38) becomes nonnegative.
Similarly to Eq. (25), a unified update rule of hil canbe obtained by substituting the equality conditions Eqs. (19),(20), and (31) for Eq. (38). Since the unified update rule
becomes complicated, we instead show the update rules ofail , bil , and cil :
Õ
ail ←
ul j
j

! − 12 ! − 23 Õ 2 Õ 2


,
fik
hi0l
hi00l
+µ


 i0 ,i

i0
i00
k


(39)
Õ
hil ul jÍbil ← −
yi j Í
,(40)
k fik gk j + l0 hil0 ul0 jjÕÕcil ← −µ Íhil2
2i0 hi0 l

fi20 k! 32! − 12Õ Õ
ki0

fi20 kfi0 k hi0l .

i0 ,iThe up-to-date variable hil is obtained by calculating Eq. (38)after the update of Eqs. (39)–(41). Note that Eqs. (38)–(41)coincide with Eq. (7) when µ = 0.
Experiment4.1

Conditions

1
hil ,
(42)
dl
ul j ← dl ul j ,(43)Íwhere dl = i hil . This normalization was performed after
the update of all the parameters in each iteration. Also,
for Log-Cos PSNMF, we performed both Eq. (26) (after the
update of hil ) and Eqs. (42) and (43) to avoid numericalinstability.
hil ←

4.2
4.2.1

Results
Hyperparameter Tuning Using Development Dataset

!
Õ

(41)

4.

X for obtaining the supervised basis matrix F.
The source-to-distortion ratio (SDR) [18] was used as
the evaluation score, which includes both the quality of theseparated target sound (absence of artificial noise) and thedegree of separation (absence of nontarget source components). The initial values of each matrix were set to uniformly distributed random values in (0, 1). The window andshift lengths in short-time Fourier transform were set to 92.9
ms and 46.4 ms, respectively. The numbers of bases for thetarget and nontarget sources were K = 27 (24 notes and threecommon spectra) and L = 50, respectively.
For Inner-Prod. PSNMF, as described in Sect. 3.1, wemust apply the following basis normalization to validate thepenalty term Pinner (F, H):

To confirm the validity of the proposed methods, we compare
the performances of simple SNMF, conventional PSNMF [9](Inner-Prod. PSNMF), Log-Cos PSNMF, and Cos PSNMF.
In the experiments, we used the artificial audio sources
produced by a Yamaha MU-1000 synthesizer in songKitamura [6], [17] as the development and test datasets. Two
instrumental melodies were selected from the eleven instruments, namely, oboe, trumpet, horn, flute, violin, clarinet,
piano, harpsichord, trombone, bassoon, and cello, and wemixed them with the same power to produce a two-source
mixture monaural signal Y , resulting in 90 mixture signals. Then, we randomly split 90 mixture signals into 45
development dataset and 45 test dataset. The development
dataset was used to find the optimal hyperparameter µ ineach method, and the test dataset was used for performancecomparison. For the sample sound in the training stage, twooctave ascending notes of the target instrument were used as

Figure 3 shows average SDR behaviors for the development dataset in each method. We can confirm that two
proposed methods, Log-Cos PSNMF and Cos PSNMF, outperform Inner-Prod. PSNMF at the optimal hyperparametersetting. This is because the penalty terms Plogcos (F, H) andPcos (F, H) appropriately force the orthogonality between
the supervised and nontarget bases in the optimization of H,
whereas the conventional penalty term Pinner (F, H) does notdirectly affect the optimization because of the scale indeterminacy in NMF; the scale of H simply decreases when we
set a large µ value. However, since the proposed methods

Fig. 3 Average SDR behaviors of development dataset. Weight coefficient µ that provides highest average SDR in development dataset is used
for performance comparison with test dataset.

IEICE TRANS. FUNDAMENTALS, VOL.E105–A, NO.6 JUNE 2022

912
Table 2p values obtained by pairwise one-sided Welch’s t test.

Method of group A
SNMF
SNMF
SNMF
Inner-Prod. PSNMF
Inner-Prod. PSNMF
Log-Cos PSNMF
Table 3

Table 1

Average and median values for test dataset.

MethodSNMF
Inner-Prod. PSNMF
Log-Cos PSNMF
Cos PSNMF

Average [dB]
6.07
7.01
7.75
7.82

Median [dB]
5.95
6.81
7.69
8.23

eliminate the scale indeterminacy problem, the optimal µ
becomes relatively peaky compared with that of the conventional PSNMF.
4.2.2

SDR Comparison Using Test Dataset

The SDR values of each method for the test dataset werecompared, where the hyperparameter µ in each PSNMF was
set to the optimal value that provides the highest scores inthe development dataset (Fig. 3). The violin plot of 45 SDR
results is shown in Fig. 4. Also, average and median valuesof the results are summarized in Table 1. We can confirm
that both Log-Cos PSNMF and Cos PSNMF outperformthe Inner-Prod. PSNMF in terms of average and medianvalues. In particular, Cos PSNMF provides the best medianvalue, whose improvement from that of Inner-Prod. PSNMF
is 1.42 dB. This result shows that the optimal settings of thehyperparameter µ in Log-Cos PSNMF and Cos PSNMF do
not strongly depend on the type of instruments.
4.2.3

Statistical Testing

To evaluate the validity of SDR differences in Table 1, twotypes of statistical testing were applied, namely, the pairwise one-sided Welch’s t test [19] and the pairwise Brunner–
Munzel (BM) test [20]. The null hypotheses of the one-sided
Welch’s t test and the BM test are respectively as follows:• The true means µA and µB of normally distributed samples in two groups A and B satisfy µA ≥ µB , where
variances of them are unequal.
• When we randomly select samples s A and sB from eachof two groups A and B, the probabilities of s A ≥ sB
and s A < sB are equal (stochastic equality).

p value [%]0.33890.00010.00003.94342.061442.8965

p values obtained by pairwise BM test.

Method of group A
SNMF
SNMF
SNMF
Inner-Prod. PSNMF
Inner-Prod. PSNMF
Log-Cos PSNMF

Fig. 4 Violin plot of 45 SDR results in test dataset. In each method,white circle indicates median value, gray vertical line shows range of 25–75
percentiles, and violin-shape curve is estimated distributions.

Method of group B
Inner-Prod. PSNMF
Log-Cos PSNMF
Cos PSNMF
Log-Cos PSNMF
Cos PSNMF
Cos PSNMF

Method of group B
Inner-Prod. PSNMF
Log-Cos PSNMF
Cos PSNMF
Log-Cos PSNMF
Cos PSNMF
Cos PSNMF

p value [%]3.15760.00560.00005.90940.981159.7709

Note that BM test does not assume the normal distribution
for observed data.
The p values obtained by each test are shown in Tables 2
and 3, respectively. When we compare the conventional andproposed methods, the null hypotheses in each test can berejected with satisfactory confidence. For example, the pvalues for comparing Inner-Prod. PSNMF and Cos PSNMF
are 2.06% and 0.98% in one-sided Welch’s t test and BM
test, respectively, which are satisfactory low. These results
show the efficacy of the proposed methods.
5.

Conclusion

In this paper, we proposed new algorithms for SNMF using the cosine similarity penalty, where the scale indeterminacy in NMF does not affect the optimization of basis
and activation matrices. From the results of experiments
using mixtures of instrumental sources, we confirmed thatthe proposed PSNMF using the cosine-similarity-based orthogonality penalty can improve the separation performancecompared with the conventional PSNMF using the innerproduct-based orthogonality penalty. Since the parametertuning for the proposed PSNMF becomes peaky, the prediction of the optimal µ is our important future work.
Although we only focus on SNMF in this paper, it isworth mentioning that the cosine similarity penalty can beapplied even in the simple (unsupervised) NMF decomposition Eq. (1) for all the bases in F. Since the orthogonality
between all the bases can be maximized, NMF with cosine
similarity can be used for a discriminative basis learning [21]and as an alternative approach to sparse NMF.
Acknowledgments
This work was partially supported by JSPS Grants-in-Aid
for Scientific Research 19K20306 and 19H01116.
References[1] D.D. Lee and H.S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, vol.401, no.6755, pp.788–
791, 1999.

IWASE and KITAMURA: SUPERVISED AUDIO SOURCE SEPARATION BASED ON NMF WITH COSINE SIMILARITY PENALTY

913[2] D.D. Lee and H.S. Seung, “Algorithms for non-negative matrix factorization,” Proc. Neural Information Processing Systems, pp.556–
562, 2000.
[3] T. Virtanen, “Monaural sound source separation by nonnegative
matrix factorization with temporal continuity and sparseness criteria,” IEEE Trans. Audio, Speech, Language Process., vol.15, no.3,
pp.1066–1074, 2007.
[4] A. Ozerov, C. Févotte, and M. Charbit, “Factorial scaled hidden
Markov model for polyphonic audio representation and source separation,” Proc. WASPAA, pp.121–124, 2009.
[5] H. Kameoka, M. Nakano, K. Ochiai, Y. Imoto, K. Kashino, and
S. Sagayama, “Constrained and regularized variants of non-negative
matrix factorization incorporating music-specific constraints,” Proc.
ICASSP, pp.5365–5368, 2012.
[6] D. Kitamura, H. Saruwatari, H. Kameoka, Y. Takahashi, K. Kondo,
and S. Nakamura, “Multichannel signal separation combining directional clustering and nonnegative matrix factorization with spectrogram restoration,” IEEE/ACM Trans. Audio, Speech, Language
Process., vol.23, no.4, pp.654–669, 2015.
[7] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari,“Determined blind source separation unifying independent vector
analysis and nonnegative matrix factorization,” IEEE/ACM Trans.
Audio, Speech, Language Process., vol.24, no.9, pp.1626–1641,
2016.
[8] P. Smaragdis, B. Raj, and M. Shashanka, “Supervised and semisupervised separation of sounds from single-channel mixtures,” Proc.
LVA/ICA, pp.414–421, 2007.
[9] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano, Y. Takahashi, andK. Kondo, “Music signal separation based on supervised nonnegative
matrix factorization with orthogonality and maximum-divergencepenalties,” IEICE Trans. Fundamentals, vol.E97-A, no.5, pp.1113–
1118, May 2014.
[10] D.R. Hunter and K. Lange, “Quantile regression via an MM algorithm,” J. Comput. Graph. Stat., vol.9, no.1, pp.60–77, 2000.
[11] D. FitzGerald, M. Cranitch, and E. Coyle, “On the use of the beta
divergence for musical source separation,” Proc. Irish Signals Syst.
Conf., 2009.
[12] P.D. O’Grady and B.A. Pearlmutter, “Discovering speech phones
using convolutive non-negative matrix factorisation with a sparsenessconstraint,” Neurocomputing, vol.72, no.1, pp.88–101, 2008.
[13] J. Le Roux, F. Weninger, and J.R. Hershey, “Sparse NMF - halfbaked or well done?,” Mitsubishi Electric Research Lab. Technical
Report, TR2015-023, 2015.
[14] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence: With application to music
analysis,” Neural Computat., vol.21, no.3, pp.793–830, 2009.
[15] A. Liutkus, D. Fitzgerald, and R. Badeau, “Cauchy nonnegative
matrix factorization,” Proc. WASPAA, 2015.
[16] N. Takahashi and M. Seki, “Multiplicative update for a class ofconstrained optimization problems related to NMF and its global
convergence,” Proc. EUSIPCO, pp.438–442, 2016.
[17] D. Kitamura, “Open dataset: songKitamura” http://d-kitamura.net/
dataset_en.html, accessed Oct. 8. 2020.
[18] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio, Speech,
Language Process., vol.14, no.4, pp.1462–1469, 2006.
[19] B.L. Welch, “The generalization of ‘Student’s’ problem when several different population variances are involved,” Biometrika, vol.34,
no.1/2, pp.28–35, 1947.
[20] E. Brunner and U. Munzel, “The nonparametric Behrens-Fisher
problem: Asymptotic theory and a small-sample approximation,”
Biometrical J., vol.42, no.1, pp.17–25, 2000.
[21] H. Nakajima, D. Kitamura, N. Takamune, H. Saruwatari, and N. Ono,
“Bilevel optimization using stationary point of lower-level objectivefunction for discriminative basis learning in nonnegative matrix factorization,” IEEE Signal Process. Lett., vol.26, no.6, pp.818–822,
2019.

Yuta Iwase
is a bachelor student of Advanced Course in Electronics, Information andCommunication Engineering, National Instituteof Technology, Kagawa College. His researchinterests include nonnegative matrix factorization and audio signal processing.

Daichi Kitamura
received the Ph.D. degree from SOKENDAI, Hayama, Japan. He
joined The University of Tokyo in 2017 as a
Research Associate, and he moved to National
Institute of Technology, Kagawa Collage as an
Assistant Professor in 2018. His research interests include audio source separation, statisticalsignal processing, and machine learning. He
was the recipient of the Awaya Prize Young Researcher Award from The Acoustical Society ofJapan (ASJ) in 2015, Ikushi Prize from Japan
Society for the Promotion of Science in 2017, Best Paper Award from IEEE
Signal Processing Society Japan in 2017, Itakura Prize Innovative Young
Researcher Award from ASJ in 2018, and IEEE Signal Processing Society(SPS) 2019 Young Author Best Paper Award from IEEE SPS in 2020.