Repeated Update of Demixing Vectors in Independent Low-rank Matrix Analysis for Better Separation Open Access

Blind source separation (BSS) aims to separate a mixture of multiple signals into individual signals [2]. It has been applied in a vast number of applications, such as multichannel audio processing [11], music information retrieval [1], wireless communication systems [10], soundscape information retrieval [9], and brain imaging [29].

Independent vector analysis (IVA) considers frequency-wise correlations and avoids permutation ambiguity by assuming a multivariate distribution of sources [3, 5]. The demixing matrices of IVA are estimated by updating with the gradient descent method, which is sensitive to a step-size parameter and may be unstable. Auxiliary-function-based IVA (AuxIVA) has been proposed to improve convergence and performance [21]. AuxIVA estimates demixing matrices faster and more stably than the conventional IVA with no tuning parameters using the majorization-minimization algorithm [8]. Some extensions of AuxIVA for online/real-time processing [27] and its application to hearing aids [26] have been proposed. To further improve the source model, AuxIVA was extended to independent low-rank matrix analysis (ILRMA) [7] by exploiting a low-rank source model of nonnegative matrix factorization. ILRMA models the spectrograms of separated sources more accurately and thus achieves a higher separation performance than AuxIVA. Many extensions to other source models of ILRMA were proposed [6, 14, 15]. Update rules of demixing matrices and source models have been extended from IVA to ILRMA.

Update rules of demixing matrices of AuxIVA have been extensively studied in recent years. Iterative projection (IP) was initially proposed for AuxIVA and updates each row vector of a demixing matrix (demixing vector) per iteration [21]. The pairwise update rule called iterative projection 2 (IP2) was also proposed and updates two demixing vectors simultaneously in each iteration [20]. We refer to the original IP as IP1 to distinguish it from IP2 in this paper. As an alternative approach, iterative source steering (ISS), which updates the entire demixing matrix using elementary row operations, has been proposed [24]. ISS can eliminate matrix inversions and thus updates more rapidly while keeping the same separation performance as IP1. All these update methods are derived by minimizing the likelihood function of the input mixture signal, leading to a system of quadratic equations called hybrid exact-approximate joint diagonalization (HEAD) [30, 31] (also known as “sequentially drilled” joint congruence (SeDJoCo) [28]). When the number of sources is three or more, HEAD is still an open problem. IP2 can yield a global optimal solution of HEAD only for two sources by solving a generalized eigenvalue problem. IP1, IP2, and ISS are available for ILRMA because ILRMA includes AuxIVA as a special case [7]. The original ILRMA employs IP1 to update demixing vectors. ILRMA with IP2 or ISS has also been proposed [17, 18]. Interestingly, it has been reported that IVA with IP2 [19] and ILRMA with IP2 [18] showed a higher performance than those with IP1, even using the same source model. These results imply that how to update the parameters also contributes the better separation performance. Similarly, it has been shown that a slower updating of source model parameters tends to improve the separation performance [13].

Motivated by these studies, we here focus on the parameter updates in ILRMA. The contribution of this study is summarized as follows:

1. We propose a new algorithm of ILRMA named ILRMA-Rep. It repeatedly updates the demixing vectors by fixing the source model at each iteration. It works as approximately updating all the demixing vectors at once.

2. We derive faster update rules of ILRMA-Rep based on the matrix inversion lemma (MIL). Especially the update rule based on IP2 with MIL, which we call IP2-MIL, is new and applicable to other BSS methods, such as online AuxIVA.

3. We conducted simulated separation experiments with 144 music signals to validate the efficacy of our proposed method. The results show that the proposed methods improve the separation performance.

The rest of this paper is organized as follows. We briefly summarize the background of the conventional ILRMA and its update rules in Section 2. Section 3 describes ILRMA with repeated updates of demixing vectors (ILRMA-Rep), the proposed ILRMA updates with repeated updates of demixing vectors. We evaluate the performance of ILRMA-Rep and compare it with that of the conventional ILRMA with IP1, IP2, and ISS in Section 4. Finally, Section 5 concludes this paper.

In the rest of this paper, we use lower- and uppercase bold symbols for vectors and matrices, respectively. Lowercase normal symbols denote vector entries or scalars. For example, x_kft denotes the kth entry of a vector x_ft. The transpose and the conjugate transpose of a vector and a matrix are respectively denoted as (·)^⊤ and $(⋅)H$⁠. We denote ℂ and ℝ₊ as the sets of all complex numbers and non-negative real numbers, respectively. Unless specified otherwise, indices f, t, and m always take the ranges from 1 to F, T, and M, respectively. Also, we omit the bounds of sums or products over these indices when they span the ranges. For example, Σ_k is the sum over all k in the range of k = 1,…,K, and Σ_ft is the double-sum over f = 1,…, F and t = 1,…, T. Similarly, {W_f}_f denotes the set of W_f for all f ; {x_kft}_kft denotes the set of x_kft for all k, f, and t, for example.

We formulate the short-time Fourier transform (STFT) domain BSS as

where x_ft∈ ℂ^M denotes the observed mixture signals recorded by M microphones and s_ft∈ ℂ^K denotes the source signals of K sources, respectively, at frequency f and time t. Here, A_f = [α_1f… α_Kf] ∈ ℂ^M×K is the mixing matrix whose kth column vector α_Kf corresponds to the steering vector from source k to each microphne. The goal of BSS is to estimate the frequency-wise demixing matrices,

such that the estimated source is

given only mixture signals x_ft, where y_ft∈ ℂ^K denotes the estimated signals. The kth row vector of the demixing matrix w_kf is called the demixing vector for source k. We henceforth consider the determined situation, i.e., the number of sources equals that of microphones M = K.

We briefly summarize ILRMA as a maximum likelihood estimation problem. The likelihood function is derived with the following assumptions.

1. All the sources are statistically independent; their joint distribution is the product of their individual distributions.

2. The estimated signals follow this zero-mean multivariate complex Gaussian distribution,

   $py(ykft)=1πrkftexp(−|ykft|2rkft),$

   (4)

   where r_kft ∈ ℝ₊ is called the variance of the distribution corresponding to the kth source in a time-frequency point f and t.

3. The variance r_kft is decomposed into L parts as

   $rkft=∑l=1Lbkflcklt,$

   (5)

   where b_kft and c_klt ∈ ℝ₊ are called the basis and activity coefficients, respectively. Henceforth, we call these coefficients the source model parameters.

The kth estimated signal is denoted as $ykft=wkfHxft$ from (3). By using these assumptions, we can determine the likelihood of the observation as

where p_x is the probability density function of the observed signals and 𝒲 is the set defined as 𝒲 = {W_f}_f. The determinant term is the Jacobian that comes from the change of variable. The goal of ILRMA is to estimate W_f by minimizing the following function:

where ℬ and 𝒞 are the sets defined as {b_kfl}_kfl and {c_klt}_klt, respectively. The aim of ILRMA is to estimate the demixing vectors w_kf and the source model parameters b_kfl,_cklt that minimize (9) with only the observed mixture x_ft.

We can derive the following multiplicative update rules of the source model parameters by applying the auxiliary-function method to (9) [7]:

These update rules guarantee the convergence of (9).

We summarized the well-used update rules of demixing vectors for AuxIVA and ILRMA, which is based on our proposed method in the next section.

We calculate the derivative of (9) with respect to w_kf (k = 1, …, K) to derive the update rules of demixing vectors. Then, we obtain the following system of quadratic equations:

where U_kf ∈ ℂ^K×K is called the weighted covariance matrix defined as

This system of quadratic equations is called hybrid exact-approximate joint diagonalization (HEAD) [28, 30, 31]. When K = 2, the closed-form solution of (12) can be derived by solving a generalized eigenvalue problem [19, 25]. However, for K ≥ 3, the closed-form solution of the HEAD problem has not yet been found. Instead, we minimize the objective function (9) with respect to only one demixing vector w_kf while keeping the others fixed. The resulting update rule is given by

It is referred to as IP1 [21]. Algorithm 1 summarizes the conventional ILRMA with IP1.

The closed-form solution of HEAD is available for two sources by solving the following generalized eigenvalue problem:

where λ_1f and λ_2f (λ_1f ≥ λ_2f) are the corresponding eigenvalues of w_1f and w_2f. Therefore, the two demixing vectors w_1f, w_2f can be updated simultaneously [20, 32].

For K ≥ 3, the joint update method for two demixing vectors in the case of three or more sources has been proposed [19]. This method was originally proposed for AuxIVA and achieved a higher separation performance with fewer iteration steps. We obtain the following system of 2K quadratic equations by calculating ∂ℒ+ /∂w_mf = 0 and ∂ℒ+ / ∂w_nf = 0 for all k ≠ m, n and m ≠ n:

The pairwise update rule of w_𝓁 (𝓁 = m, n) is derived by solving (16)-(18) as follows [19];

where z_𝓁f is the eigenvector of Z_𝓁f (see [20] for details). Note that the choice of the two indices m and n is arbitrary as long as m ≠ n.

Instead of estimating the demixing vector w_kf in IP1, ISS [24] updates the entire demixing matrix by estimating a new vector v_kf = [v_1kf… v_Kkf]^⊤,

where the update rule that minimizes (9) with respect to v_mkf is given by

Furthermore, from the demixing model $ykft=wkfHxft$⁠, the following are the inverse-free update rules of v_kf and the output estimated signal y_ft:

IVA with IP2 [19] and ILRMA with IP2 [18] showed a higher performance than those with IP1. These results imply that the performance will be further improved if more demixing vectors, ideally all of them, are updated simultaneously such as

However, it leads to the HEAD problem and no closed-form solution has yet been found for K ≥ 3 as mentioned in the previous section.

Instead of solving the HEAD problem in a closed-form manner, we propose an “approximately” simultaneous update of demixing vectors by simply repeating the demixing vector updates with IP1, IP2, or ISS several times while keeping the source model parameters fixed. Since each update reduces the objective function of (9), it is expected that the demixing matrix converges the solution of the HEAD problem by repeating the demixing matrix update. We can see this as follows. From (12), we rewrite the HEAD problem as

where Ε_κ is the K × K identity matrix. Therefore, we can check how close the demixing matrix is to the solution of HEAD by visualizing the left side of (27). Figure 1 shows one example of convergence to the solution of HEAD by iteration steps. Each colormap shows the element-wise absolute value of the left-hand side of (27). The initial demixing matrices were set to the identity matrices. The weighted covariance matrices were set to the random Hermitian matrices. As shown in Figure 1, the repeated updates of demixing vectors considerably improve the estimated solution of HEAD. This example implies that a set of repeated updates of demixing vectors can work as a simultaneous update of demixing vectors.

The outline of the proposed procedures is summarized in Algorithm 2. We call this algorithm ILRMA-Rep. Figure 2 illustrates how the demixing vectors and the source model parameters are updated.

Updates of demixing matrices with IP1 or IP2 include a matrix inversion. It could cause a large computational complexity for repeating the update of demixing vectors. However, the matrix inversion lemma (MIL) is available to reduce complexities because the covariance matrices are fixed within each “repeat.” In this section, we derive an efficient algorithm for the repeated update of demixing vectors.

MIL states that:

where X is an invertible square matrix of shape K × K, Ψ and Ω are matrices of shape K × J and J × K (J ≤ K) respectively, and E_J is the J × J identity matrix.

Let us consider avoiding the matrix inversion of IP1 included in (13). Let $Af=Wf−1$ and a_kf= A_fe_k, if we update A_f together with W_f, (13) can be easily calculated as $w^kf=Ukf−1akf$ where $w^kf$ denote the left side of (13). Then, the point is how to update A_f without the matrix inversion. We rewrite (13) as

Then, by using the matrix inversion lemma with Ψ ← e_k and $Ω←(w^kf−wkf)H$ in (28), we obtain the update rule of A_f without matrix inversion. Finally, the update rules without the matrix inversion are summarized as follows.

These update rules were first derived in the context of online BSS [27]. Henceforth, we call this version IP1-MIL.

Likewise, by setting Ψ ← [e_me_n] and $Ω←[(w^mf−wmf)(w^nf−wnf)]H$ in (28), the efficient update rule of IP2 is given by

where E₂ is the 2 × 2 identity matrix. Note that G is a matrix of shape 2 × 2 and thus the cost of its inversion is much less than those of W_f and U_kf. To the best of our knowledge, this is the first derivation of IP2 with the MIL. Henceforth, we call this version IP2-MIL.

This subsection outlines the time complexity of the demixing matrix update method for each iteration and frequency bin. For IP1 and IP2, the calculation of (W_fU_kf)⁻¹ in (13) and (19), respectively, takes O(K³) time at each source k. For ISS by (23), the calculation of v_kf takes O(K³) time. Therefore, the total time complexity of updating the demixing matrix W_f for all sources k = 1, …, K is O(K⁴) with these methods. In ILRMA-Rep, the covariance matrices U_kf are fixed, and the repeat of demixing matrix updates takes O(N_repK⁴) time, where N_rep is the number of repeats. The use of MIL allows for efficient updating of W⁻¹. More specifically, O(K⁴) can be computed for the first repeat only, and O(K³) for subsequent (N_rep − 1) repeats. Therefore, the total time complexity of IP1-MIL and IP2-MIL for all the sources takes the time complexity of O(N_repK³+ K⁴). Note that MIL cannot be applied to ISS because W⁻¹ is not calculated in ISS. In addition, we have an alternative for ISS using (24) instead of (23). As for ISS by (24), the calculation of v_kf takes O(TK) time, and the total time complexity of updating the demixing matrix W_f for all sources k = 1,…, K is O(TK²). From this, ILRMA-Rep using ISS by (24) takes O(N_repTK²) in total. Repeating the update by (24) should take more time than ISS by (23) because T ≫ K² in most cases for normal BSS. Therefore, we hereafter use (23) for ISS update in this paper. Table 1 summarizes the complexities of the update rules.

Since the runtime performance may differ from theoretical complexity, as described in Section 3.3, we first compare the runtime performance characteristics of the different methods. This simulation was run using Python language on a workstation powered by a 96-core AMD EPYC 7643 processor at 2.3 GHz. We generate sets of K random Hermitian matrices U_kf (k = 1, …, K) by the following procedure:

1. Generate random complex values

   $ykft,zkft∼N(0,1)(∀k,f,t),xkft:=ykft+jzkft(∀k,f,t).$

2. Generate random observed signals as x_ft := [x_1ft … x_Kft]^⊤ (∀f, t).

3. Calculate random covariance matrices as

   $Ukf←1T∑txftxftH∑f|xkft|2(∀k,f).$

The initial demixing matrices were set to the identity matrix. We set the number of time frame T to 100 and the number of frequency bins F to 257. For ISS, we update the demixing matrices by (23) in this experiment because T is large. We ran 20 iterations of IP1, IP1-MIL, IP2, IP2-MIL, and ISS by changing K = 2, 4, 6, 8, 10, 12, 14,16. The covariance matrices are fixed in each K. Table 2 shows the average runtimes for IP1, IP1-MIL, IP2, IP2-MIL, and ISS over 100 trials.

We can see that the fastest algorithm depends on the number of channels. For example, ISS is the best in K = 2, but IP1-MIL is the best in K = 6. The table also shows that IP1-MIL and IP2-MIL were almost consistently faster than IP1 and IP2, except for IP2 and IP2-MIL in K = 2. This exception might be because computing temporal variables in IP2-MIL take more time than the matrix inversion when K = 2. In the following evaluation of separation performance, IP1-MIL and IP2-MIL will respectively be referred to as IP1 and IP2 because the resulting outputs are identical to those of the original IP1 and IP2.

We evaluated the separation and convergence performance characteristics of our proposed ILRMA-Rep with convolutive music mixtures. The separation performance criterion is the scale-invariant source-to-distortion ratio (SI-SDR) [S.], defined as

where s ∈ ℝ^L is the reference signal and ŝ∈ ℝ^L is the estimated signal in the discrete time domain of length L. A higher SI-SDR indicates higher separation performance and quality. We also measured SI-SDR improvement (SI-SDRi) defined as the difference in SI-SDR between the estimated and mixture signals:

where x∈ ℝ^L is the mixture signal in the discrete time domain. In the following experiment, we use the individual reverberant signals at each microphone as reference signals.

We used the MUSDB18-HQ dataset [22] as a source signal, including 150 stereo-recorded songs. Each song in the dataset consists of four parts: bass, drums, other, and vocals. To perform the image-source method and reduce the computational load, we extracted the left channel of each signal and then downsampled each signal by 16 kHz from the original 44.1 kHz. We created simulated convolutive mixtures using the Python package pyroomacoustics [23]. As shown in Figure 3, the number of sound sources and microphones was four, and the microphone array was uniformly linear with a spacing of 2 cm. The reverberation time was approximately 200 ms.

We use the STFT with a 4096-points Hamming analysis window and half-overlap. We set the initial demixing matrices W_f to the identity matrices for all f. Also, we set the initial values of the source models b_kfi to 1 and c_klt to uniformly distributed values over ]0, 1] for all l, k, f, and t. We performed ILRMA-Rep with 100 iteration steps and several “repeates” N_rep in each iteration. After separation, the scale of the estimated signal was restored by back-projection onto the first microphone [16]. Then, we calculated the SI-SDRi for each music mixture and obtained its average along all the songs and channels.

We implemented all the algorithms in Python 3.9.6. For numerical stability, we performed the following regularization after the respective updates in each iteration:

where ε is a user-defined parameter. We set ε to 10⁻⁶ in this experiment.

Figure 4 shows the separation performance at the end of the 100 iterations for the number of bases L =1, 2, 5, 10, and 20. Note that the conventional ILRMA corresponds to N_rep = 1 in the brightest color bar in Figure 4, and the case of L = 1 in the leftmost plot in Figure 4 is nearly equivalent to that of AuxIVA.

For L =1, the SI-SDRi of all methods did not change markedly with increasing N_rep. For L = 2 or larger, SI-SDRi was steadily increased by N_rep in each L. These results support that performance can be improved by simply repeating the updates of demixing matrices.

One possible explanation for these results would be as follows. ILRMA is a difficult optimization problem in terms of two different types of parameters: demixing matrices and source model parameters (NMF parameters). It is a non-convex optimization, and then it should have local minima. For example, if the source model is optimized too fast, it may fit mixed (or insufficiently-separated) signals. Then, the demixing matrices are no longer updated, and a poor result might be obtained. Previous studies suggest that quickly updating the demixing matrices [18, 19] or slowly updating the source model parameters [13] leads to better separation performance. Our proposed method repeatedly updates the demixing matrices by fixing the source model parameters. It works by updating the demixing matrices faster than the source model parameters and thus improves the separation performance.

Figure 5 shows the relationship between SI-SDRi at the end of the iteration with N_rep = 1 and N_rep ≥ 2. Each point corresponds to the value before becoming the average in Figure 4. The dashed line in the figure is the line where the vertical and horizontal axes are equal; the points above the dashed line indicate that repeats of the demixing vector update improve SI-SDRi. For L =1, the points are clustered near the dashed line for any number of repetitions. By contrast, when L = 2 or larger, more points are clustered above the dashed line, and this tendency is stronger with a larger number of iterations. This result supports the observation that the separation performance improves with repeated updates of the demixing matrix.

Figure 6 shows the separation performance with the number of iteration steps for the number of bases L = 1, 2, 5, 10, and 20. For L =1, the convergence became faster as N_rep increased, but the final performance was almost the same, approximately 3dB. By contrast, for L = 2 or larger, both the convergence speed and the final performance became much higher than those for L = 1.

In this paper, we proposed a repeated update scheme of demixing vectors for independent low-rank matrix analysis. We also derived an efficient update rule for the repeated version of IP2 by applying the matrix inversion lemma. Music source separation experiments with more than 100 songs were performed to evaluate the separation performance. The experimental results implied that the proposed update scheme improved the separation performance. In particular, we experimentally found that when the number of bases is large, the separation performance is markedly improved by repeating parameter updates. In future work, we will further investigate the performance of repeated updates of demixing vectors on other source models using IP1, e.g., in independent deeply learned matrix analysis [12] or the multichannel variational autoencoder method [4].

Taishi Nakashima received his B.E. in Engineering from Osaka University, Osaka, Japan, in 2019 and his M.S. in Informatics from Tokyo Metropolitan University, Tokyo, Japan, in 2021. He is pursuing a Ph.D. at Tokyo Metropolitan University and is also a recipient of the JSPS Research Fellowship (DC1) from April 2021. He is an esteemed Student Member of the Acoustical Society of Japan (ASJ) and the IEEE Signal Processing Society (SPS). He received the 24th Best Student Presentation Award of ASJ and the 16th IEEE SPS Japan Student Conference Paper Award in 2022. His research interests primarily focus on blind source separation and acoustic signal processing.

Nobutaka Ono received his B.E., M.S., and Ph.D. degrees from the University of Tokyo, Japan, in 1996, 1998, and 2001, respectively. He became a research associate in 2001 and a lecturer in 2005 at the University of Tokyo. He moved to the National Institute of Informatics in 2011 as an associate professor and then to Tokyo Metropolitan University in 2017 as a full professor. His research interests include acoustic signal processing, especially microphone array processing, source localization and separation, machine learning, and optimization algorithms. He is a member of IEEE, EURASIP, APSIPA, IPSJ, IEICE, and ASJ. He was a member of IEEE Audio and Acoustic Signal Processing (AASP) Technical Committee from 2014 to 2019. He served as Associate Editor of IEEE Transactions on Audio, Speech, and Language Processing from 2012 to 2015. He received the best paper award at APSIPA ASC in 2018 and 2021 and Sadaoki Furui Prize Paper Award from APSIPA in 2021.