Block Coordinate Descent Algorithms for Auxiliary-Function-Based Independent Vector Extraction

Recently, maximum likelihood approaches have been proposed for BSE in which the background noise components are modeled as stationary Gaussians. These methods include independent vector extraction (IVE [21, 22, 23]) and overdetermined IVA (OverIVA [24, 25, 26, 27]), which will be collectively referred to as IVE in this paper. When $K$ non-Gaussian sources are all super-Gaussian (as defined in Assumption 2), IVE can use a majorization-minimization (MM [28]) optimization algorithm developed for an auxiliary-function-based ICA (AuxICA [29]) and an auxiliary-function-based IVA (AuxIVA [30, 31, 32]). In this paper, we only deal with an IVE based on the MM algorithm [24, 25, 26, 27, 23] and always assume that the super-Gaussian distributions are defined by Assumption 2 in Section III-A.

In the MM-based IVE, each separation matrix $W\in\mathbb{C}^{M\times M}$ is optimized by iteratively minimizing a surrogate function of the following form:

Here, $\bm{w}_{1},\ldots,\bm{w}_{K}\in\mathbb{C}^{M}$ are filters that extract the target-source signals and $W_{\operatorname*{z}}$ is a filter that extracts the $M-K$ noise components. $V_{1},\ldots,V_{K},V_{\operatorname*{z}}\in\mathbb{C}^{M\times M}$ are the Hermitian positive definite matrices that are updated in each iteration of the MM algorithm. Interestingly, when $K=M$ or when viewing the second term of $g(W)$ as $\operatorname*{tr}\left(W_{\operatorname*{z}}^{h}V_{\operatorname*{z}}W_{\operatorname*{z}}\right)=\sum_{i=K+1}^{M}\bm{w}_{i}^{h}V_{\operatorname*{z}}\bm{w}_{i}$, the problem of minimizing $g(W)$ has been discussed in ICA/IVA literature [33, 34, 35, 36, 37]. Among the algorithms developed so far, block coordinate descent (BCD [38]) methods with simple analytical solutions have attracted much attention in the field of audio source separation because they have been experimentally shown to give stable and fast convergence behaviors. A family of these BCD algorithms [29, 30, 31, 32], summarized in Table I, is currently called an iterative projection (IP) method. The IP method has been widely applied not only to ICA and IVA but also to audio source separation methods using more advanced source models (see [39, 40, 41, 42, 43] for the details of such methods).

When we consider directly applying the IP methods developed for AuxIVA to the BSE problem of minimizing $g(W)$ with respect to $W$, AuxIVA-IP1 [30] in Table I, for instance, will cyclically update $\bm{w}_{1}\rightarrow\cdots\rightarrow\bm{w}_{K}\rightarrow\cdots\rightarrow\bm{w}_{M}$ one by one. However, in the BSE problem, the signals of interests are $K$ non-Gaussian sources, and most of the optimization of $W_{\operatorname*{z}}$ should be skipped.

Therefore, in a previous work [24], an algorithm (IVE-OC in Table I) was proposed that cyclically updates $\bm{w}_{1}\rightarrow W_{\operatorname*{z}}\rightarrow\cdots\rightarrow\bm{w}_{K}\rightarrow W_{\operatorname*{z}}$ one by one with a computationally cheap formula for $W_{\operatorname*{z}}$. In this work [24], the updating equation for $W_{\operatorname*{z}}$ was derived solely from the (weak) orthogonal constraint (OC [22, 21, 23]) where the sample correlation between the separated $K$ non-Gaussian sources and the $M-K$ background noises equals zero. (Note that the $K$ non-Gaussian sources are not constrained to be orthogonal in the OC.) Although the algorithm has successfully reduced the computational cost of IVA by nearly a factor of $K/M$, the validity of imposing the heuristic OC remains unclear.

After that, IVE-OC has been extended by removing the OC from the model. One such extension is a direct method that can obtain a global minimizer of $g(W)$ when $K=1$ [25, 27]. The other extension is a BCD algorithm for $K\geq 2$ that cyclically updates the pairs $(\bm{w}_{1},W_{\operatorname*{z}})\rightarrow\cdots\rightarrow(\bm{w}_{K},W_{\operatorname*{z}})$ one by one [26], but the computational cost is not so cheap due to the full update of $W_{\operatorname*{z}}$ in each iteration. These algorithms are called IVE-IP2 in this paper (see Table I).

In this work, we propose BCD algorithms for IVE, which are summarized in Table I with comparisons to the previous BCDs (IPs). The followings are the contributions of this paper.

(i) We speed up the previous IVE-IP2 for $K\geq 2$ by showing that $W_{\operatorname*{z}}$’s update can be omitted in each iteration of BCD without changing the behaviors of the algorithm (Section IV-F). This is attained by carefully exploiting the stationary Gaussian assumption for the noise components. In an experiment of speaker source enhancement, we confirmed that the computational cost of the conventional IVE-IP2 is consistently reduced.

(ii) We provide a comprehensive explanation of IVE-IP2 for $K=1$ (Sections IV-D and IV-H). Interestingly, it turns out to be a method that iteratively updates the maximum signal-to-noise-ratio (MaxSNR [44, 45]) beamformer and the power spectrum for the unique ($K=1$) target-source signal. The experimental result shows that this algorithm has much faster convergence than conventional algorithms. Note that IVE-IP2 for $K=1$ was developed independently and simultaneously in our conference paper [27] and by Scheibler–Ono [25].

(iii) We reveal that IVE-OC [24], which was developed with the help of the heuristic orthogonal constraint (OC), can also be obtained as a BCD algorithm for our proposed IVE without the OC (Section IV-G). (This interpretation of IVE-OC as a BCD algorithm was described in our conference paper [27], but it is not mathematically rigorous. We here provide a rigorous proof of this interpretation.)

(iv) We further extend the proposed IVE-IP2 for the semiblind case where the transfer functions for $L$ ($1\leq L\leq K$) sources of interests can be given a priori (Section V). We call the proposed semiblind method Semi-IVE. In Semi-IVE, $L$ separation filters, e.g., $\bm{w}_{1},\ldots,\bm{w}_{L}$, which correspond to the known transfer functions are optimized by iteratively solving the linear constrained minimum variance (LCMV [44]) beamforming algorithm, and the remaining $\bm{w}_{L+1},\ldots,\bm{w}_{K}$ (and $W_{\operatorname*{z}}$) are optimized by the full-blind IVE-IP2 algorithm. We experimentally show that when $L\geq K-1$ transfer functions are given Semi-IVE yields surprisingly fast convergence.

Organization: The rest of this paper is organized as follows. The BSE and semi-BSE problems are defined in Section II. The probabilistic model for the proposed IVE is compared with related methods in Section III. The optimization algorithms are developed separately for the BSE and semi-BSE problems in Sections IV and V. The computational time complexity for these methods is discussed in Section VI. In Sections VII and VIII, experimental results and concluding remarks are described.

Let $\mathcal{S}_{+}^{d}$ or $\mathcal{S}_{++}^{d}$ denote the set of all Hermitian positive semidefinite (PSD) or positive definite (PD) matrices of size $d\times d$. Let $\operatorname*{GL}(d)$ denote the set of all (complex) nonsingular matrices of size $d\times d$. Also, let $\bm{0}_{d}\in\mathbb{C}^{d}$ be the zero vector, let $O_{i,j}\in\mathbb{C}^{i\times j}$ be the zero matrix, let $I_{d}\in\mathbb{C}^{d\times d}$ be the identity matrix, and let $\bm{e}_{k}$ be a vector whose $k$-th element equals one and the others are zero. For a vector or matrix $A$, $A^{\top}$ and $A^{h}$ represent the transpose and the conjugate transpose of $A$. The element-wise conjugate is denoted as $A^{\ast}=(A^{\top})^{h}$. For a matrix $A\in\mathbb{C}^{i\times j}$, $\operatorname*{Im}A$ is defined as the subspace $\{A\bm{x}\mid\bm{x}\in\mathbb{C}^{j}\}\subset\mathbb{C}^{i}$. For a vector $\bm{v}$, $\|\bm{v}\|=\sqrt{\bm{v}^{h}\bm{v}}$ denotes the Euclidean norm.

Throughout this paper, we discuss IVA and IVE using the terminology from audio source separation in the short-term Fourier transform (STFT) domain.

Suppose that $K$ super-Gaussian target-source signals and a stationary Gaussian noise signal of dimension $M-K$ are transmitted and observed by $M$ sensors. In this paper, we only consider the case where $1\leq K<M$. The observed signal in the STFT domain is modeled at each frequency bin $f=1,\ldots,F$ and time-frame $t=1,\ldots,T$ as

where $s_{i}(f,t)\in\mathbb{C}$ and $\bm{z}(f,t)\in\mathbb{C}^{M-K}$ are the STFT coefficients of target source $i\in\{1,\ldots,K\}$ and the noise signal, respectively. $\bm{a}_{i}(f)\in\mathbb{C}^{M}$ is the (time-independent) acoustic transfer function (or steering vector) of source $i$ to the sensors, and $A_{\operatorname*{z}}(f)\in\mathbb{C}^{M\times(M-K)}$ is that of the noise. It is assumed that the source signals are statistically mutually independent.

In the blind source extraction (BSE) problem, we are given an observed mixture $\bm{x}=\{\bm{x}(f,t)\}_{f,t}$ and the number of target sources $K$. From these inputs, we seek to estimate the spatial images $\bm{x}_{1},\ldots,\bm{x}_{K}$ for the target sources, which are defined as $\bm{x}_{i}(f,t)\coloneqq\bm{a}_{i}(f)s_{i}(f,t)\in\mathbb{C}^{M}$, $i=1,\ldots,K$.

In the semiblind source extraction (Semi-BSE) problem, in addition to the BSE inputs, we are given $L$ transfer functions $\bm{a}_{1},\ldots,\bm{a}_{L}$ for $L$ super-Gaussian sources, where $1\leq L\leq K$. From these inputs, we estimate all the target-source spatial images $\bm{x}_{1},\ldots,\bm{x}_{K}$. If $L=K$, then Semi-BSE is known as a beamforming problem.

Motivation to address Semi-BSE: In some applications of audio source extraction, such as meeting diarization [46], the locations of some (but not necessarily all) point sources are available or can be estimated accurately, and their acoustic transfer functions can be obtained from, e.g., the sound propagation model [47]. For instance, in a conference situation, the attendees may be sitting in chairs with fixed, known locations. On the other hand, the locations of moderators, panel speakers, or audience may change from time to time and cannot be determined in advance. In such a case, by using these partial prior knowledge of transfer functions, Semi-BSE methods can improve the computational efficiency and separation performance of BSE methods. In addition, since there are many effective methods for estimating the transfer function of (at least) a dominant source [48, 49, 50], there is a wide range of applications where Semi-BSE can be used to improve the performance of BSE.

In the mixing model (1), the mixing matrix

is square and generally invertible, and hence the problem can be converted into one that estimates a separation matrix $W(f)\in\operatorname*{GL}(M)$ satisfying $W(f)^{h}A(f)=I_{M}$, or equivalently, satisfying

where we define

Denote by $\bm{s}_{i}(t)=[\,s_{i}(1,t),\ldots,s_{i}(F,t)\,]^{\top}\in\mathbb{C}^{F}$ the vector of all the frequency components for source $i$ and time-frame $t$. The proposed IVE exploits the following three assumptions. Note that Assumption 2 was introduced for developing AuxICA [29] and AuxIVA [30, 31, 32].

Assumption 3 plays a central role for deriving several efficient algorithms for IVE. Despite this assumption, as we experimentally confirm in Section VII, the proposed IVE can extract speech signals even in a diffuse noise environment where the noise signal is considered super-Gaussian or nonstationary and has an arbitrary large spatial rank.

With the model defined by (6)–(13), the negative loglikelihood, $g_{0}(W)\coloneqq-\frac{1}{T}\log p(\bm{x}\mid W)$, can be computed as

where $W\coloneqq\{W(f)\}_{f=1}^{F}$ are the variables to be optimized.

If we assume that the $M-K$ noise components also independently follow super-Gaussian distributions, then the IVE model coincides with that of IVA [5] or AuxIVA [30, 31, 32]. The following are the two advantages of assuming the stationary Gaussian model (13).

(i) As we confirm experimentally in Section VII, when we optimize the IVE model, separation filters $\bm{w}_{1},\ldots,\bm{w}_{K}$ extract the top $K$ highly super-Gaussian (or nonstationary) signals such as speech signals from the observed mixture while $W_{\operatorname*{z}}$ extracts only the background noises that are more stationary and approximately follow Gaussian distributions. On the other hand, in IVA, which assumes super-Gaussian noise models, $K$ ($<M$) target-source signals need to be chosen from the $M$ separated signals after optimizing the model.

(ii) As we reveal in Section IV-G, in the proposed IVE, it suffices to optimize $\operatorname*{Im}W_{\operatorname*{z}}=\{W_{\operatorname*{z}}\bm{v}\mid\bm{v}\in\mathbb{C}^{M-K}\}\subset\mathbb{C}^{M}$, the subspace spanned by $W_{\operatorname*{z}}\in\mathbb{C}^{M\times(M-K)}$, instead of $W_{\operatorname*{z}}$. Because we can optimize $\operatorname*{Im}W_{\operatorname*{z}}$ very efficiently (Section IV-G), IVE can reduce the computation time of IVA, especially when $K\ll M$.

The proposed IVE is inspired by OverIVA with an orthogonal constraint (OC) [24]. This conventional method will be called IVE-OC in this paper.

IVE-OC [24] was proposed as an acceleration of IVA [30] for the case where $K<M$, while maintaining its separation performance. The IVE-OC model is defined as the proposed IVE by replacing the noise model from (13) to (16) and introducing two additional constraints:

The first constraint (18) may be applicable because there is no need to extract the noise components (see [24, 26] for details). The second constraint (19), called an orthogonal constraint (OC [21]), was introduced to help the model distinguish between the target-source and noise signals.

OC, which forces the sample correlation between the separated target-source and noise signals to be zero, can equivalently be expressed as

which together with (18) imply

where we define

From (18) and (23), it turns out that $W_{\operatorname*{z}}$ is uniquely determined by $W_{\operatorname*{s}}$ in the IVE-OC model. Hence, in a paper on IVE-OC [24], an algorithm was proposed in which $W_{\operatorname*{z}}$ is updated based on (18) and (23) immediately after updating any other variables $\bm{w}_{1},\ldots,\bm{w}_{K}$ to always impose OC on the model. Although the algorithm was experimentally shown to work well, its validity from a theoretical point of view is unclear because the update rule for $W_{\operatorname*{z}}$ is derived solely from the constraints (18)–(19) and does not reflect an objective of the optimization problem for parameter estimation, such as minimizing the negative loglikelihood.

In this paper, we develop BCD algorithms for the maximum likelihood estimation of the proposed IVE that does not rely on OC, and identify one such algorithm (IVE-IP1 developed in Section IV-G) that exactly coincides with the conventional algorithm for IVE-OC. This means that OC is not essential for developing fast algorithms in IVE-OC. Due to removing OC from the model, we can provide other more computationally efficient algorithms for IVE, which is the main contribution of this paper.

We briefly describe how to apply the conventional MM technique developed for AuxICA [29] to the proposed IVE.

In IVE as well as AuxICA/AuxIVA, a surrogate function $g$ of $g_{0}$ is designed with an auxiliary variable $r$ that satisfies

Then, variables $W$ and $r$ are alternately updated by iteratively solving

for $l=1,2,\ldots$ until convergence. In the same way as in AuxIVA [30], by applying Proposition 5 in Appendix A to the first term of $g_{0}(W)$, problem (28) in IVE comes down to solving the following $F$ subproblems:

where $r^{(l)}\coloneqq\{r_{i}^{(l)}(t)\in\mathbb{R}_{\geq 0}\}_{i,t}$ and

Here, the computation of (33)–(34) corresponds to the optimization of (27). Recall that $G^{\prime}$ in (33) is the first derivative of $G\colon\mathbb{R}_{\geq 0}\to\mathbb{R}$ (see Assumption 2). To efficiently solve (29), we propose BCD algorithms in the following subsections. From the derivation, it is guaranteed that the objective function is monotonically nonincreasing at each iteration in the MM algorithm.¹¹1 Due to space limitations, showing the convergence rate and other convergence properties of the proposed algorithms will be left as a future work.

No algorithms have been found that obtain a global optimal solution for problem (29) for general $K,M\in\mathbb{N}$. Thus, iterative algorithms have been proposed to find a local optimal solution. Among them, a family of BCD algorithms has been attracting much attention (e.g., [29, 30, 31, 32, 24, 26, 25, 27, 39, 40, 41, 42]) because they have been experimentally shown to work faster and more robustly than other algorithms such as the natural gradient method [52]. The family of these BCD algorithms (specialized to solve (29)) is currently called an iterative projection (IP) method.

As we will see in the following subsections, all the IP algorithms summarized in Table I can be developed by exploiting the stationary condition, which is also called the first-order necessary optimality condition [53]. To simplify the notation, when we discuss (29), we abbreviate frequency bin index $f$ without mentioning it. For instance, $W(f)$ and $V_{i}(f)$ are simply denoted as $W$ and $V_{i}$ (without confusion).

To rigorously derive the algorithms, we always assume the following two technical but mild conditions (C1) and (C2) for problem (29):²²2 If (C1) is violated, problem (29) has no optimal solutions and algorithms should diverge to infinity (see [54, Proposition 1] for the proof). Conversely, if (C1) is satisfied, it is guaranteed that problem (29) has an optimal solution by Proposition 6 in Appendix A. In practice, the number of frames $T$ exceeds that of sensors $M$, and (C1) holds in general. Condition (C2) is satisfied automatically if we initialize $W$ as nonsingular. Intuitively, singular $W$ implies $-\log|\det W|=+\infty$, which will never occur during optimization.

(C1)

$V_{1},\ldots,V_{K},V_{\operatorname*{z}}\in\mathcal{S}_{+}^{M}$ are positive definite.

(C2)

Estimates of $W\in\mathbb{C}^{M\times M}$ are always nonsingular during optimization.

We focus on the case where $K=1$ and derive the proposed algorithm IVE-IP2³³3 The algorithm IVE-IP2 for $K=1$ was developed independently and simultaneously by Scheibler–Ono (called the fast IVE or FIVE [25]) and the authors [27] in the Proceedings of ICASSP2020. In this paper, we give a complete proof for the derivation of IVE-IP2 and add some remarks (see also Section IV-H for the projection back operation). that is expected to be more efficient than IVE-IP1. The algorithm is summarized in Algorithm 3.

When $(K,M)=(1,2)$ or $(2,2)$, problem (29) can be solved directly through a generalized eigenvalue problem [35, 31, 32]. We here extend this direct method to the case where $K=1$ and $M\geq 2$ in the following proposition.

By Proposition 2, under condition (C1), a global optimal solution for problem (29) can be obtained by updating $W=[\bm{w}_{1},W_{\operatorname*{z}}]$ using (41)–(44) with $Q_{1}=1$ and $Q_{\operatorname*{z}}=I_{M-1}$ and choosing the generalized eigenvalue $\lambda$ in (44) as the largest one. Moreover, this algorithm can be accelerated by omitting the computation of (42)–(43) and updating only $\bm{w}_{1}$ according to (41) and (44). It is applicable for extracting the unique target-source signal because the formulas for computing $\bm{w}_{1}$, $V_{1}$, and $V_{\operatorname*{z}}$ are independent of $W_{\operatorname*{z}}$. The obtained algorithm IVE-IP2 is shown in Algorithm 3.

Interestingly, because $V_{\operatorname*{z}}$ and $V_{1}$ can be viewed as the covariance matrices of the mixture and noise signals, the update rules (41) and (44) turn out to be a MaxSNR beamformer [44, 45]. Hence, IVE-IP2 can be understood as a method that iteratively updates the MaxSNR beamformer $\bm{w}_{1}$ and target-source signal $\bm{s}_{1}$.

Suppose $2\leq K<M$. The conventional IVE-IP2 developed in [26, Algorithm 2] is a BCD algorithm that cyclically updates the pairs $(\bm{w}_{1},W_{\operatorname*{z}}),\ldots,(\bm{w}_{K},W_{\operatorname*{z}})$ by solving the following subproblem for each $i=1,\ldots,K$ one by one:

It was shown in [26, Theorem 3] that a global optimal solution of problem (45) can be obtained through a generalized eigenvalue problem. In this paper, we simplify the result of [26, Theorem 3] in the next proposition, based on which we will speed up the conventional IVE-IP2 in Section IV-F.

In the conventional IVE-IP2, under condition (C2), problem (45) is solved by (46)–(50) where $B_{i,\operatorname*{z}}\coloneqq[\bm{b}_{i},B_{\operatorname*{z}}]$ is computed through the following generalized eigenvalue problem instead of using (51)–(52):