DETERMINED BLIND SOURCE SEPARATION VIA PROXIMAL SPLITTING ALGORITHM
Kohei Yatabe† and Daichi Kitamura‡∗†‡

Department of Intermedia Art and Science, Waseda University, Tokyo, Japan
Department of Information Physics and Computing, The University of Tokyo, Tokyo, Japan
ABSTRACTThe state-of-the-art algorithms of determined blind source separation (BSS) methods based on the independent component analysis(ICA) have gained computational efficiency by the majorizationminimization (MM) principle with a price of losing flexibility. Thatis, replacing and comparing different source models are not easy insuch MM-based framework because it requires efforts to derive a
new algorithm each time when one changes the model. In this paper,a general framework for obtaining an ICA-based BSS algorithm is
proposed so that a source model can easily be replaced because onlya single line of the algorithm must be modified. A sparsity-based extension of the independent vector analysis and a low-rankness-based
BSS model using the nuclear norm are also proposed to demonstrate
the simplicity and easiness of the proposed framework.
Index Terms— Independence-based separation, frequency domain independent component analysis (FDICA), independent vector
analysis (IVA), primal-dual splitting algorithm, proximity operator.
1. INTRODUCTION
Blind source separation (BSS) is methodology for recovering source
signals from multiple mixtures without any knowledge about the
mixing system. Let a convolutive mixing process be approximatedin time-frequency domain asx[t, f ] ≈ A[f ]s[t, f ],(1)

where x = [x1 , x2 , . . . xM ]T is an observation obtained by M microphones, s = [s1 , s2 , . . . sN ]T is a source signal to be recovered,A[f ] is an M ×N mixing matrix, and t and f are indices of time andfrequency, respectively. Then, the aim of BSS is to recover N source
signals s from the mixtures x. In a determined or overdeterminedsituation (M ≥ N ), many of the BSS problems are formulated as
an estimation problem of finding an N × M demixing matrix W [f ]which is a left inverse of A[f ] (i.e., W [f ]A[f ] = I), and the source
signals are recovered by simple multiplication:
W [f ]x[t, f ] ≈ W [f ]A[f ]s[t, f ] = s[t, f ].

(2)

For the sake of simplicity, only a determined situation (M = N ) is
considered in this paper.
For estimating a demixing matrix W [f ], statistical independence between source signals is often assumed that leads to a
family of independence-based BSS algorithms. Arguably, independent component analysis (ICA) [1] applied in frequency domain(FDICA) [2–6] is one of the most famous methods among them.
However, FDICA suffers from the so-called permutation problem [7–10], and thus some recent developments on BSS aim to∗ This work was partly supported by JSPS Grant-in-Aid for Research Ac-

tivity Start-up (17H06572,17H07191).

avoid it by considering more sophisticated models of source signals.
For instance, independent vector analysis (IVA) [11–13] assumes
co-occurrence among the frequency components in each source, andindependent low-rank matrix analysis (ILRMA) [14–16] assumes
low-rankness on spectrogram of each source. The key to successof these methods is to incorporate prior knowledge of source signals into their formulations. That is, improvement brought by thesemethods relies on the preciseness of their source models. Therefore,seeking a better model is the important process for developing anovel and effective BSS method.
However, recent algorithms [16–18] cannot be applied to adifferent source model directly because they are specialized toeach method. These state-of-the-art algorithms are based on the
majorization-minimization (MM) principle which requires speciallydesigned upper-bound of the objective functions. That is, one has toderive a new algorithm each time as the source model is modified.
Therefore, it may take a lot of time to examine a new source model
especially when the model is a complicated one. If a single algorithm can handle a large number of source models without effort,
discovering a better source model should become much easier thatpossibly boosts the development of BSS.
In this paper, a flexible framework for independence-based BSS
is proposed based on a proximal splitting algorithm [19–22]. The
usefulness of the splitting algorithm comes from its capability of
splitting an optimization problem into several easier subproblems
which are handled by the proximity operators. We take advantage ofthis feature to split the ICA-based BSS problem into two parts so thatdifferent source models can easily be combined by modifying onlya single line of the algorithm. The proposed framework is tested
by introducing three new source models (a sparsity-based extension
of IVA, a low-rankness-based model using the nuclear norm, and asparsity-based extension of the low-rank model), and potentiality ofthe proposed method is indicated by the result.
2. INDEPENDENCE-BASED BSS
As introduced in the previous section, independence-based BSS
methods aim to estimate M × M demixing matrices {W [f ]}F
f =1
which approximately recover the source signals from the observations as W [f ]x[t, f ] ≈ s[t, f ]. Many of them fall into a minimization problem of the following form:
Minimize P(W [f ]x[t, f ]) −{W [f ]}F
f =1

F
X

log |det(W [f ])|,

(3)

f =1

where P is a real-valued penalty function corresponding to thesource model. For example, with some constant C,P(y[t, f ]) = C ky[t, f ]k1 = C

M X
T X
F
X
m=1 t=1 f =1

|ym [t, f ]|(4)

recovers the traditional FDICA with the Laplace distribution as thesource model (which is separable for each f ), andP(y[t, f ]) = C ky[t, f ]k2,1 = C

M X
T X
F
X
m=1 t=1|ym [t, f ]|2

1
2

(5)

f =1

obtains IVA whose source model is the spherical Laplace distribution. ILRMA can also be interpreted as Eq. (3) withP(y[t, f ]) = CMXDR (ym [t, f ]),

(6)

m=1

where DR (ym [t, f ]) is a measure of low-rankness based on the
Itakura–Saito non-negative matrix factorization (IS-NMF) [23]:DR (ym [t, f ]) =
(7)


T
F
R
2
XX
X [m] [m]|ym [t, f ]|
min
+ log
ϕf,r ψr,t .
PR[m] [m][m][m]
ϕf,r ≥0,ψr,t ≥0 t=1 f =1
r=1
r=1 ϕf,r ψr,t
From this perspective, it is clear that the difference of performance
among these methods is owing to goodness of the penalty function
P. Therefore, performance of a BSS method can be improved byfinding a better source model and the corresponding penalty function P. For seeking a better model, it is convenient to have a single
algorithm that can handle a large number of source models without
spending time for its derivation.

which accepts not only smooth functions but also non-differentiable
functions, including sparsity inducing norms in Eqs. (4) and (5), andnon-finite functions such as − log in Eq. (3). That is, difficulty associated with properties of the functions are eliminated within theproximity operator thanks to the proximity term. Therefore, difficulty of an optimization problem of the form Eq. (8) only depends
on computational complexity of the corresponding proximity operators. Fortunately, proximity operators of several functions related to
ICA-based BSS can be computed quite efficiently as follows [19,20]:py + y 2 + 4µ
prox−µ log [ y ] =
,
(11)
2



µ
proxµk·k1 [y] m [t, f ] = 1 −ym [t, f ],(12)|ym [t, f ]| +
proxµk·k2,1 [y] m [t, f ] =µym [t, f ],(13)
1 − PF1( f =1 |ym [t, f ]|2 )2 +where (·)+ = max{0, ·}.
3.2. Reformulating optimization problem of ICA-based BSS
To apply the primal-dual splitting algorithm, Eq. (3) is reformulated
into the form of Eq. (8). Firstly, for considering the proximity operator, the second term is modified. Since determinant of a matrix
can be expressed in terms of the singular values as |det(W [f ])| =QM
2m=1 σm (W [f ]), Eq. (3) can be rewritten as follows :3. PROPOSED FRAMEWORKIn this section, a general algorithm for solving the BSS problems in
the form of Eq. (3) is proposed by applying a primal-dual splitting
algorithm to the reformulated version of Eq. (3).
3.1. Primal-dual splitting algorithm
Let us consider a minimization problem of the form:
Minimize g(w) + h(Lw),w(8)

where g and h are proper lower-semicontinuous convex functions,and L is a matrix. A primal-dual splitting algorithm [21] solves thisproblem by iterating the following procedure1 :
 [k]


H [k]
 w
,
 e = proxµ1 g w − µ1 µ2 L y

[k]
[k]
 z = y + L(2we − w ),

(9)
 
 e
 y = z − proxh/µ2 z ,e ye ) + (1 − α)(w[k] , y[k] ),(w[k+1] , y[k+1] ) = α(w,where µ1 > 0 and µ2 > 0 are step sizes, and 2 > α > 0 is a relaxation factor which can adjust the speed of convergence (α = 1 isthe standard speed, and α > 1 accelerates and α < 1 slows down thealgorithm). The important feature of this algorithm is that each function in the problem is handled through the proximity operator [20],ih1ky − zk22 ,(10)proxµg [y] = arg min g(z) +z
2µ1 For the sake of simplicity, details of the algorithm are omitted in thispaper. The reader may refer to [21] and references therein for detailed explanations and the condition for convergence. Note that, although several
primal-dual algorithms can handle more complicated problem than Eq. (8),this paper only focus on that form because it is sufficient for the proposal.

Minimize P(W [f ]x[t, f ]) −{W [f ]}F
f =1

F X
M
X

log σm (W [f ]),(14)

f =1 m=1

where σm (W [f ]) is the mth singular value of W [f ].
Then, the optimization variables {W [f ]}F
f =1 are vectorized toform a single vector. Let w be an M 2 F -dimensional vector corresponding to the demixing filters {W [f ]}F
f =1 ,
T
T
T T
w = [w[1] , w[2] , . . . , w[F ] ] ,w[f ] = V(W [f ]) (15)where V is a linear operator converting a matrix into a vector,
V(W [f ]) = [W1,1 [f ], . . . , W1,M [f ], W2,1 [f ], . . . , WM,M [f ]]T, (16)and let M be a linear operator converting it back into the matrix,
M(w)[f ] = W [f ],(17)which will be used as w[f ] = V(M(w[f ])) = V(M(w)[f ]). With
these notations, Eq. (14) can be expressed as follows:
Minimize P(Xw) −
w

F X
M
X

log σm (M(w)[f ]),(18)

f =1 m=1

where X is a matrix constructed from the observed data x[t, f ] as
X = blkdiag(χ[1], χ[2], . . . , χ[F ]),χ[f ] = blkdiag(χ[f ], χ[f ], . . . , χ[f ]), (M times)
χ[f ] = [τ1 [f ], τ2 [f ], . . . , τM [f ]],

(19)
(20)
(21)

τm [f ] = [xm [1, f ], xm [2, f ], . . . , xm [T, f ]]T ,(22)

2 Note that, while Eq. (3) is only defined for square matrices, the extended
formulation in Eq. (14) allows rectangular demixing matrices (i.e., N 6= M ).
Therefore, the proposed method does not require the pre-processing method
using principle component analysis (PCA) in an over-determined situation,
which might improve the performance as discussed in [15].

blkdiag(·) is an operator constructing a block-diagonal matrix byconcatenating inputted matrices diagonally, τm [f ] is T × 1, χ[f ] is
T ×M , χ[f ] is M T ×M 2 , and X is F M T ×F M 2 .
Let the second term in Eq. (18) be shortly denoted by I:
I(w) = −F X
MXlog σm (M(w)[f ]).

(23)

Algorithm 1 PDS-BSS
1: Input: X, w[1] , y[1] , µ1 , µ2 , α
2: Output: w[K+1]
3: for k = 1, . . . , K do
e = proxµ1 I [ w[k] − µ1 µ2 X H y[k] ]
4:
w
5:
6:

f =1 m=1

Then, Eq. (18) can be rewritten as follows:
Minimize I(w) + P(Xw),w(24)

which is in the same form of Eq. (8).

e − w[k] )z = y[k] + X(2we = z − prox 1 P [ z ]yµ2

7:y[k+1] = αey + (1 − α)y[k][k+1]
e + (1 − α)w[k]
8:
w
= αw
9: end for

Algorithm 2 PDS-BSS-multiPenalty
[1]

[1]

1: Input: X, w[1] , y1 , . . . , yQ , µ1 , µ2 , α

3.3. Proposed algorithm
Since the BSS problem is reformulated into the form of Eq. (8), theprimal-dual splitting algorithm in Eq. (9) can be applied to it if theproximity operator of each function is efficiently computable.
It is known that a proximity operator of an orthogonally invariant function can be evaluated by applying the corresponding proximity operator to the singular values of the matrix [20]. By regarding− log σm in Eq. (23) as − log |σm |, the proximity operator of I(w)is obtained as the similar form to the singular value thresholding:e(proxµI [w])[f ] = V( U Σ(M(w)[f]) V H ),(25)

2: Output: w[K+1]
3: for k = 1, . . . , K do
PQ[k] e = proxµ1 I [ w[k] − µ1 µ2 X H
4:w]
q=1 yq
5:
for q = 1, . . . , Q do[k]e − w[k] )6:zq = yq + X(2w7:eq = zq − prox 1 Pq [ zq ]
y
[k+1]

µ2

[k]

= αe
yq + (1 − α)yq
8:
yq
9:
end fore + (1 − α)w[k]10:w[k+1] = αw
11: end forwhere W = U ΣV H is the singular value decomposition of W ,eΣ(W) = diag(prox−µlog [σ1 (W )], . . . , prox−µlog [σM (W )]),
(26)
prox−µlog [·] is in Eq. (11), and diag(·) is the operator constructing
a diagonal matrix from inputted scalars. In other words, applying
the proximity operator of −µ log to each singular value of W [f ], foreach frequency independently, gives proxµI [·]. It is worth mentioning that this operation is stable because it does not magnify kwk2
much [see Eq. (11)] in contrast to MM algorithms which involve
inversions of matrices that sometimes lead to instability.
By using Eq. (25) as the main building block, the primal-dual
splitting algorithm for BSS problems is obtained as in Algorithm 1.
In the 6th line, the proximity operator of P is required. For FDICA
and IVA in the form of Eqs. (4) and (5), the proximity operators are
given in Eqs. (12) and (13), respectively. Any other penalty function
P can be incorporated by only changing proxP/µ2 [·] in that line,
and therefore this algorithm can be used to test performance of several source models without efforts on modifying the code. Note thatan iterative algorithm can be used to evaluate proxP/µ2 [·] since it
is defined through the optimization problem in Eq. (10). That is, it
is still possible to apply the proposed algorithm even when the corresponding proxP/µ2 [·] does not admit a closed form solution. For
example, the penalty function of ILRMA in Eq. (6) yields the optimization problem similar to IS-NMF which could be solved by theexisting algorithms. Therefore, by using the proposed algorithm, incorporating existing source models into ICA-based BSS should beeasy if the computational complexity is not a concern.
3.4. Proposed algorithm with multiple penalty functions
A source model consists of two or more penalty functions,
Q

Minimize I(w) +
w

X
q=1

Pq (Xw),(27)

can also be handled by the proposed algorithm. Algorithm 1 is easily extended to deal with this problem by vertically concatenating thematrix X as L = [X T , . . . , X T ]T (Q times). Then, the primal-dual
splitting algorithm for Eq. (27) can be derived by simply applying
Eq. (9) to its vertically concatenated version, which is summarizedin Algorithm 2. As presented in the algorithm, each penalty function
Pq is independently handled by the corresponding proximity operator. Therefore, this algorithm is applicable to a complicated source
model consisting of several simple functions {Pq }Q
q=1 . Such models may include sparse+low-rank model of robust PCA [24, 25] andmulti-group sparsity of harmonic/percussive source separation [26].
3.5. Normalization of data matrix
As introduced in Section 2, an ICA-based BSS method usually has a
specific constant C which is derived from the corresponding statistical model. Although this constant is important, it may not be easy
to derive C for some non-standard models. To circumvent this complication, a normalization method convenient for the proposed algorithm is considered here. Since the step sizes µ1 and µ2 in Eq. (9)should be chosen to satisfy µ1 µ2 kLk2s ≤ 1 [21]3 , they can be set to
one if kLks = 1, where k · ks denotes the spectral norm. That is, if
the data matrix X is normalized asp
e = X/( Q kXk ),(28)
Xsthen the step sizes can be set to µ1 = 1 and µ2 = 1. Remind that
α can be arbitrarily chosen from (0, 2) [21] or chosen as 1 whichbypasses the last line of Eq. (9). Thus, with this normalization, there
is nothing to worry about in regards to the choice of the parameters.
3 Because of the − log σm term, the problem in Eq. (24) is non-convex,
which requires an additional analysis for convergence. Such analysis will beconsidered in the future work, and thus we omitted theoretical details in thispaper. Nevertheless, experimental results in the next section show that thischoice of step sizes, derived for convex problems, seems work properly.

5.5

3.5

Proposed

5

4.5
0

100

200

Mixture A

300

4

400

Number of iterations

500

Proposed

3

2.5
0

100

200

Mixture B

300

400

Number of iterations

500

sparse & low-ranksparse IVA

6

5.5
0

6

SDR improvement [dB]6

6.5

AuxIVASDR improvement [dB]4.5

AuxIVASDR improvement [dB]SDR improvement [dB]

6.5

5.5
5

4.5

low-rank
200

400

Mixture A

600

800

Number of iterations

1000

sparse & low-ranksparse IVA low-rank

4

3.5
3
0

200

400

Mixture B

600

800

Number of iterations

1000

Fig. 1. Comparison between AuxIVA [18] and the proposed algorithm. The black lines with circles denotes AuxIVA, while the bluelines with cross marks denotes the proposed algorithm.

Fig. 2. Performance of several models solved by the proposed algorithm. Each curve represents each model proposed in Section 4.2,
while black and blue lines are the ones from Fig. 1 for comparison.

4. EXPERIMENTS

to sparsify the spectrograms of separated signals by inducing group
sparsity through k · k2,1 , while I prevents W [f ] from being rank deficient. Any other penalty function possibly improves the accuracy
if it promotes a property of source signals appropriately [29].
One of such properties is the low-rankness of a spectrogram [14–16]. Then, it may be natural to consider the nuclear norm k · k∗because it is well-known as a function inducing low-rankness of amatrix. Therefore, we propose the following model (low-rank):

The proposed algorithm was tested by applying it to speech mixtures. The database used in the experiment was a part of SiSEC2011[27] called dev1 of the UND task4 . Live recording (liverec)of four female speech sources recorded by two microphones (5 cm
spacing) was chosen as the test data. For making the problem determined, two pairs of sources were considered: Mixture A consistsof two sources arrived from −50◦ and 45◦ and Mixture B consistsof two sources arrived from −10◦ and 15◦ , where 0◦ corresponds
to the normal direction to the microphone array. The reverberation
time was 130 ms, and 128-ms-long Hann window with 64-ms shift
was used. The initial value of demixing matrices w[1] was set to theidentity matrices (W [f ] = I for all f ). The parameters in the proposed algorithm were µ1 = 1, µ2 = 1, and α = 1.75, and all initialvalues of y in Algorithms 1 and 2 were set to the zero vector.
4.1. Comparison with MM algorithm
In order to compare with the state-of-the-art MM algorithm in [18],the proposed algorithm was applied to IVA whose source model isthe spherical Laplace distribution:e
Minimize I(w) + kXwk
2,1 ,
w

w

(30)

e into the correwhere M∗ is the operator converting the vector Xw
sponding matrices (spectrograms of the separated signals), and thenuclear norm is evaluated for each spectrogram. The proximity operator of k · k∗ is the famous singular-value thresholding [20] which
also easily adapts the proposed algorithm to this problem.
Further, it can be presumed that an additional sparsity-inducing
term improves separation performance. Thus, to demonstrate thebenefit of flexibility of the proposed algorithm, sparse IVA,ee
Minimize I(w) + kXwk
2,1 + λ kXwk1 ,w(31)

and a sparse & low-rank model,(29)

where k · k2,1 is in Eq. (5). By simply incorporating the proximity
operator in Eq. (13) to Algorithm 1, the proposed method is adapted
to this problem. The performance at each iteration count was measured by the signal-to-distortion ratio (SDR) [28]. In the MM algorithm [18] which is denoted by AuxIVA, the statistical constant C
in Eq. (5) was correctly considered.
Figure 1 shows the results of IVA solved by AuxIVA and the
proposed algorithm. AuxIVA reached rapidly to specific values of
SDR, where 50 or 100 iteration seems enough for these cases. Although the proposed algorithm required more iterations, it reached tothe same values with several hundreds iterations. Running time periteration of each method was 80.8 ms for AuxIVA and 47.4 ms forthe proposed algorithm (with Core i5-7200U processor and MATLAB 2017a). This result indicates that the proposed algorithm with
the normalization rule in Section 3.5 properly works as the state-ofthe-art algorithm with comparable running time.
4.2. Application to several source models
To demonstrate simplicity and easiness of the proposed algorithm,
three additional formulations with different source models are considered here. From Eq. (29), it can be clearly seen that IVA attempts
4 Available at http://sisec2011.wiki.irisa.freMinimize I(w) + kM∗ (Xw)k∗,ee
Minimize I(w) + kM∗ (Xw)k∗ + λ kXwk1 ,w(32)

are also proposed, where λ > 0 is a weighting parameter balancing
the influence of the norms.
Figure 2 shows the results of these three models solved by the
proposed algorithm5 (λ = 0.002) together with IVA in Fig. 1 forcomparison. For both mixtures, sparse IVA and sparse & low-rank
models resulted in higher SDR than IVA. Although the modificationof the proposed algorithm was able to be done within a few minutes, we were able to test the performance of these four BSS models.
This fact illustrates the advantage of the proposal in this paper which
could contribute to discovery of a new BSS model.
5. CONCLUSIONSFor solving ICA-based BSS problems, this paper proposed a generalalgorithm which admits a complicated source model consisting ofmultiple penalty terms. It can handle a large number of models witha slight modification that allows quick investigation of the model’s
performance. For illustrative examples, three models based on sparsity and low-rankness of the spectrograms were also proposed.
√
5 For this experiment,
2 in Eq. (28) was omitted for Eqs. (31) and (32)in order to set the trajectories of the four models to the similar ones for easiercomparison. It should not be omitted in practice for the sake of reliability.

6. REFERENCES[1] P. Comon, “Independent component analysis, A new concept?,” Signal Process., vol. 36, no. 3, pp. 287–314, 1994.
[2] P. Smaragdis, “Blind separation of convolved mixtures in thefrequency domain,” Neurocomputing, vol. 22, no. 1, pp. 21–34,
1998.
[3] S. Araki, R. Mukai, S. Makino, T. Nishikawa, and
H. Saruwatari, “The fundamental limitation of frequency
domain blind source separation for convolutive mixtures ofspeech,” IEEE Trans. Speech Audio Process., vol. 11, no. 2,
pp. 109–116, Mar. 2003.
[4] H. Sawada, R. Mukai, S. Araki, and S. Makino, “Convolutive blind source separation for more than two sources in thefrequency domain,” in Proc. IEEE Int. Conf. Acoust., Speech,
Signal Process., May 2004, vol. 3, pp. 885–888.
[5] H. Buchner, R. Aichner, and W. Kellermann, “A generalizationof blind source separation algorithms for convolutive mixtures
based on second-order statistics,” IEEE Trans. Speech Audio
Process., vol. 13, no. 1, pp. 120–134, Jan. 2005.
[6] H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee, andK. Shikano,“Blind source separation based on a fastconvergence algorithm combining ICA and beamforming,”
IEEE Trans. Audio, Speech, Lang. Process., vol. 14, no. 2, pp.
666–678, Mar. 2006.
[7] S. Kurita, H. Saruwatari, S. Kajita, K. Takeda, and F. Itakura,“Evaluation of blind signal separation method using directivity
pattern under reverberant conditions,” in Proc. IEEE Int. Conf.
Acoust., Speech, Signal Process., 2000, vol. 5, pp. 3140–3143.
[8] N. Murata, S. Ikeda, and A. Ziehe, “An approach to blindsource separation based on temporal structure of speech signals,” Neurocomputing, vol. 41, no. 1, pp. 1–24, 2001.
[9] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust and precise method for solving the permutation problem
of frequency-domain blind source separation,” IEEE Trans.
Speech Audio Process., vol. 12, no. 5, pp. 530–538, Sep. 2004.

[15] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Relaxation of rank-1 spatial constraint in
overdetermined blind source separation,” in Proc. Eur. Signal
Process. Conf., Aug. 2015, pp. 1261–1265.
[16] 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, Lang. Process., vol.
24, no. 9, pp. 1626–1641, Sep. 2016.
[17] N. Ono and S. Miyabe, “Auxiliary-function-based independent
component analysis for super-gaussian sources,” in Proc. Int.
Conf. Latent Variable Anal. Signal Separation, 2010, pp. 165–
172.
[18] N. Ono, “Stable and fast update rules for independent vector
analysis based on auxiliary function technique,” in Proc. IEEE
Workshop Appl. Signal Process. Audio Acoust., Oct. 2011, pp.
189–192.
[19] P. L. Combettes and J.-C. Pesquet, Proximal Splitting Methods
in Signal Processing, pp. 185–212, Springer, New York, 2011.
[20] N. Parikh and S. Boyd, “Proximal algorithms,” Found. Trends
Optim., vol. 1, no. 3, pp. 127–239, 2014.
[21] N. Komodakis and J. C. Pesquet, “Playing with duality: An
overview of recent primal-dual approaches for solving largescale optimization problems,” IEEE Signal Process. Mag., vol.
32, no. 6, pp. 31–54, Nov. 2015.
[22] M. Burger, A. Sawatzky, and G. Steidl, First Order Algorithms
in Variational Image Processing, pp. 345–407, Springer,
Cham, 2016.
[23] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix
factorization with the Itakura–Saito divergence. With application to music analysis,” Neural Comput., vol. 21, no. 3, pp.
793–830, 2009.
[24] P. S. Huang, S. D. Chen, P. Smaragdis, and M. HasegawaJohnson, “Singing-voice separation from monaural recordings
using robust principal component analysis,” in Proc. Int. Conf.
Acoust., Speech, Signal Process., 2012, pp. 57–60.

[10] H. Sawada, S. Araki, and S. Makino, “Measuring dependence of bin-wise separated signals for permutation alignment
in frequency-domain BSS,” in Proc. IEEE Int. Symp. Circuits
Syst., May 2007, pp. 3247–3250.

[25] Y. Ikemiya, K. Itoyama, and K. Yoshii, “Singing voice separation and vocal F0 estimation based on mutual combination ofrobust principal component analysis and subharmonic summation,” IEEE/ACM Trans. Audio, Speech, Lang. Process., vol.
24, no. 11, pp. 2084–2095, Nov. 2016.

[11] A. Hiroe, “Solution of permutation problem in frequency domain ICA, using multivariate probability density functions,” inProc. Int. Conf. Independent Compon. Anal. Blind Source Separation, 2006, pp. 601–608.

[26] N. Ono, K. Miyamoto, J. Le Roux, H. Kameoka, and
S. Sagayama, “Separation of a monaural audio signal into harmonic/percussive components by complementary diffusion on
spectrogram,” in Proc. Eur. Signal Process. Conf., 2008.

[12] T. Kim, T. Eltoft, and T.-W. Lee, “Independent vector analysis: An extension of ica to multivariate components,” in Proc.
Int. Conf. Independent Compon. Anal. Blind Source Separation, 2006, pp. 165–172.

[27] S. Araki, F. Nesta, E. Vincent, Z. Koldovský, G. Nolte,
A. Ziehe, and A. Benichoux, “The 2011 signal separation evaluation campaign (SiSEC2011): - audio source separation -,” inProc. Int. Conf. Latent Variable Anal. Signal Separation, 2012,
pp. 414–422.

[13] T. Kim, H. T. Attias, S. Y. Lee, and T. W. Lee, “Blind source
separation exploiting higher-order frequency dependencies,”
IEEE Trans. Audio, Speech, Lang. Process., vol. 15, no. 1, pp.
70–79, Jan. 2007.
[14] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Efficient multichannel nonnegative matrix factorization exploiting rank-1 spatial model,” in Proc. IEEE Int.
Conf. Acoust., Speech Signal Process., Apr. 2015, pp. 276–280.

[28] E. Vincent, R. Gribonval, and C. Fevotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio, Speech, Lang. Process., vol. 14, no. 4, pp. 1462–1469, Jul.
2006.
[29] K. Yatabe and Y. Oikawa, “Phase corrected total variation foraudio signals,” in Proc. IEEE Int. Conf. Acoust., Speech Signal
Process., Apr. 2018.