CACGMM

今天谈一下多通道分离算法里面一个比较成熟的系列，spatial clustering。在声源分布比较分散的情况下，这个算法有着稳定的表现，mask估计也比较准确。NTT在这块做了很多工作，在CHiME3到CHiME5的历次测评中，此算法的合理使用均拿到了十分不错的结果。因为聚类算法都是无监督的，所以相比用网络去估计mask，它的一个好处是不存在数据层面的mismatch，加上复杂度上不算高，online实现也比较简单，因此我觉得它是一套具有优势和潜力的算法。

我很早之前（应该是17年的时候）就接触过此类算法，是在CHiME4数据上做CGMM的mask估计^[1]。当时文章看完之后，照猫画虎也能很快的把程序写出来，跑出一个结果。代码上问题主要有两个，一是慢，二是数值上不稳定。程序写的慢是我最难以忍受的两个地方之一（其二是写的丑），加上当时要处理的数据有点多，所以在代码层面就是想办法减少循环（python和matlab道理一样），算法上就是化简运算，最后多任务并行进行提交。关于数值稳定性，因为涉及到埃尔米特矩阵的求逆，遇到奇异矩阵，EM过程就直接退出了，所以代码中一般需要对这些情况进行预防和处理。还有一点就是聚类过程在每个sub-band上是完全独立的，所以聚类的结果在不同的频率上是不一致的，存在一个置换问题，一般是需要一个后处理重新对齐之后才能使用的。CGMM在处理增强问题中（也就是估计speech&noise两类mask），通过合理的初始化，可以免去做这类后处理的需要，因此我后续一般都只用CGMM做多通道的语音增强，而没有在分离任务中使用。

Paderborn那边有几个人一直在用CACGMM^[2]做分离任务，比如在CHiME5里面的GSS和network的unsupervised训练，最近集中看了一下他们的code^[3]和文章^[4,5,6]，并把setk中的cluster部分重写了一下，也顺便写几点小结在这里分享一下。

CG vs CACG

这两者的EM过程在化简之后是完全一致的。CGMM的建模对象是直接的观测向量$\mathbf{y}_{tf}$，假定分布服从0均值的CGMM：

$$\mathbf{y}_{tf} \sim \sum_k \alpha_f^k \mathcal{N}(\mathbf{y}_{tf}; 0, \phi_{tf}^k\mathbf{R}_f^k) \tag{1}$$

其中

$$\begin{align} \mathcal{N}(\mathbf{y}_{tf}; 0, \phi_{tf}^k \cdot \mathbf{R}_f^k) &= \frac{1}{\det(\phi_{tf}^k \cdot \mathbf{R}_f)} e^{-\mathbf{y}_{tf}^H (\phi_{tf}^k\mathbf{R}_f^k)^{-1}\mathbf{y}_{tf}} \notag \\ & = \frac{1}{\det(\mathbf{R}_f) \cdot (\phi_{tf}^k)^M} e^{-\mathcal{K}_{tf}^k/ \phi_{tf}^k} \tag{2} \end{align}$$

注：在写分布表达式的时候，我把常数项（比如$\pi$之类的）都去掉了，不影响$\gamma_{tf}^k$的计算和$Q$函数值的相对大小，$\mathcal{K}_{tf}^k = \mathbf{y}_{tf}^H (\mathbf{R}_f^k)^{-1}\mathbf{y}_{tf}$

我们最终要的TF-mask就是E步中得出的后验概率$\gamma_{tf}^k$

$$\gamma_{tf}^k = \frac{\alpha_f^k \mathcal{N}(\mathbf{y}_{tf}; 0, \phi_{tf}^k\mathbf{R}_f^k)}{\sum_c \alpha_f^k \mathcal{N}(\mathbf{y}_{tf}; 0, \phi_{tf}^k\mathbf{R}_f^k) } \tag{3}$$

CGMM在M步的参数更新为：

$$\begin{align} \alpha_f^k & = \frac{1}{T} \sum_t \gamma_{tf}^k \notag \\ \mathbf{R}_f^k & = \frac{\sum_t \gamma_{tf}^k \cdot \mathbf{y}_{tf} \mathbf{y}_{tf}^H \cdot(\phi_{tf}^k)^{-1}} {\sum_t \gamma_{tf}^k} \tag{4} \\ \phi_{tf}^k & = \frac{1}{M} \mathbf{y}_{tf}^H (\mathbf{R}_f^k)^{-1}\mathbf{y}_{tf} = \frac{\mathcal{K}_{tf}^k}{M} \notag \end{align}$$

可以发现$(4)$中的$\phi_{tf}^k$带入$(2)$，可得

$$\mathcal{N}(\mathbf{y}_{tf}; 0, \phi_{tf}^k\mathbf{R}_f^k) = \frac{e^{-M}}{\det(\mathbf{R}_f) \cdot (\mathcal{K}_{tf}^k)^M} \tag{5}$$

由此，$\gamma_{tf}^k$和$\mathbf{R}_f^k$的表达式也可以继续简化为：

$$\gamma_{tf}^k = \frac{\alpha_f^k}{\det(\mathbf{R}_f^k)(\mathcal{K}_{tf}^k)^M} \cdot \left(\sum_k \frac{\alpha_f^k}{\det(\mathbf{R}_f^k)(\mathcal{K}_{tf}^k)^M} \right)^{-1} \tag{6}$$

和

$$\mathbf{R}_f^k = M \cdot \frac{\sum_t \gamma_{tf}^k \cdot \mathbf{y}_{tf} \mathbf{y}_{tf}^H \cdot(\mathcal{K}_{tf}^k)^{-1}} {\sum_t \gamma_{tf}^k} \tag{7}$$

所以从式子$(6,7)$来看，已经没有必要维护$\phi_{tf}^k$这个中间参数了，而它们正是CACGMM的M步，根据$(5)$式，可知道新的分布可以表示为

$$\mathcal{A}(\mathbf{z}_{tf}; \mathbf{B}_f) = \frac{1}{\det(\mathbf{B}_f) \cdot \left(-\mathbf{z}_{tf}^H \mathbf{B}_f^{-1}\mathbf{z}_{tf} \right)^M} \tag{8}$$

$\mathbf{B}_{f}$主要用于区分$\mathbf{R}_f$。从参数量的角度来说，$\mathcal{A}(\cdot)$相比$\mathcal{N}(\cdot)$少维护了一个$\phi_{tf}^k$。此外还有一点是CACGMM的建模对象是完全的空间信息 $\mathbf{z}_{tf} = \mathbf{y}_{tf} / |\mathbf{y}_{tf}|$，可以证明在$\mathbf{y}_{tf} \sim \mathcal{N}(\cdot)$时，$\mathbf{z}_{tf} \sim \mathcal{A}(\cdot)$。

Implementation

这类算法的训练都是EM的过程，每个类别和每个sub-band上的更新都是相对独立的，所以比较直观的写法一定有两个循环（外层频率$F$，内层混合模型的类别$K$）。循环带来的问题就是效率很低，而逻辑上本身是允许在$F\times K$个类别层面上并行的，这里推荐使用十分强大的np.einsum函数进行多维矩阵和向量的运算。关于它网上有一些样例和介绍，我在这里简单举两个例子说明一下如何使用它计算上述几个公式中的变量值。

先提前声明几个变量和对应的shape

• obs $\mathbf{y} \in \mathbf{C}^{F \times M \times T}$
• gamma $\gamma \in \mathbf{R}^{K \times F \times T}$
• kernel $\mathcal{K} \in \mathbf{R}^{K \times F \times T}$

1. 计算式$(7)$中的$\mathbf{R}_{f}^k$

   denominator = np.sum(gamma, -1, keepdims=True)
   # (K x F x T) @ (F x M x T) @ (F x M x T) => (K x F x M x M)
   R = M * np.einsum("...t,...xt,...yt->...xy", gamma / kernel, obs, obs.conj())
   R = R / denominator[..., None]

2. 计算$\mathcal{K}_{tf}^k = \mathbf{y}_{tf}^H(\mathbf{R}_f^k)^{-1}\mathbf{y}_{tf}$

   # (F x M x T) @ (F x K x M x M) @ (F x M x T) => (K x F x T)
   K = np.einsum("...xt,...xy,...yt->...t", obs.conj(), R_inv, obs)

3. 求$(\mathbf{R}_f^k)^{-1}$和$\det(\mathbf{R}_f^k)$

   解决上面这些变量的计算之后，还有一个求逆的问题，即$(\mathbf{R}_f^k)^{-1}$。$(\mathbf{R}_f^k)^{-1}$是一个Hermitian矩阵，但 Hermitian矩阵不一定可逆，举一个极端的例子，$T = 1$，那么$\mathbf{R}_f^k$的秩为1，则不可逆（矩阵奇异）。

   实际运算中可以通过保存特征值和特征向量的方法存储协方差矩阵，这样方便矩阵求逆和行列式的计算。对于$\mathbf{R}_{tf}^k$做特征值分解得：

   $$\mathbf{R}_{tf}^k = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^H$$

   $\mathbf{\Lambda}$为特征值对角矩阵，如果全正，则$\mathbf{R}_{tf}^k$为正定矩阵，一定可逆，且特征向量正交，满足$\mathbf{U}\mathbf{U}^H = \mathbf{I}$，由此可知：

   $$(\mathbf{R}_{tf}^k)^{-1} = \mathbf{U}(\mathbf{\Lambda})^{-1}\mathbf{U}^H$$

   根据“矩阵的秩等于特征值的乘积”：

   $$\log \det(\mathbf{R}_{tf}^k) = \log \sum \text{diag}(\mathbf{\Lambda})$$

   在$\mathbf{R}_{tf}^k$不可逆时，$\mathbf{\Lambda}$不能保证非负，因此在EM过程中可以强制加一个最小非负约束。此外通常的解决方法一般都是用diag loading，就是对矩阵加上一个对角量$\epsilon \mathbf{I}$，这种方式不能完全避免矩阵的奇异，而且$\epsilon$的选取也比较关键。

从实现上基本就是这几个地方需要注意一下，整体下来训练速度上明显提升不少，结果也比较正常。CGMM里面如果解决增强问题，初始化可以用：

$$\begin{align} \mathbf{R}_{tf}^n &= \mathbf{I} \notag \\ \mathbf{R}_{tf}^s &= \sum_t \mathbf{y}_{tf}^H\mathbf{y}_{tf} \cdot 1/T \notag \end{align}$$

但是也不能完全避免permutation问题的出现，可以在一些样例里面看到一些bad case。CACGMM的初始化完全随机，EM输出的$\gamma$需要在过一遍align算法，这部分算法的原理尚未细看，但是对最终结果的影响是很大的（一旦align处理的不好，完全看不出估计的mask中谱的形状，下面放了一个比较正常的样例，cluster2是noise成分），感兴趣的可以参考一下这里。最终的代码也已经经过测试，merge到setk中去了。

Reference

[1] Higuchi T, Ito N, Araki S, et al. Online MVDR beamformer based on complex Gaussian mixture model with spatial prior for noise robust ASR[J]. IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2017, 25(4): 780-793.
[2] Ito N, Araki S, Nakatani T. Complex angular central Gaussian mixture model for directional statistics in mask-based microphone array signal processing[C]//2016 24th European Signal Processing Conference (EUSIPCO). IEEE, 2016: 1153-1157.
[3] https://github.com/fgnt/pb_bss
[4] Drude L, Hasenklever D, Haeb-Umbach R. Unsupervised training of a deep clustering model for multichannel blind source separation[C]//ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019: 695-699.
[5] Drude L, Heymann J, Haeb-Umbach R. Unsupervised training of neural mask-based beamforming[J]. arXiv preprint arXiv:1904.01578, 2019.
[6] Kanda N, Boeddeker C, Heitkaemper J, et al. Guided Source Separation Meets a Strong ASR Backend: Hitachi/Paderborn University Joint Investigation for Dinner Party ASR[J]. arXiv preprint arXiv:1905.12230, 2019.