NMF及语音分离

这一部分在之前的笔记中曾经有过讨论，在这里重新整理，补充缺失的内容。

非负矩阵分解 NMF基本原理

数据通常表示为矩阵形式：

字典学习，低秩近似，因子分析，潜在语义分析等方法，通常会将/将矩阵分解为字典矩阵和行为矩阵的乘积，以实现降维，拆分（源分离），插值等工作。

非负矩阵分解要求：

数据 $V$ 以及因子 $W,H$ 具有非负的值。

因子 $W$ 的非负性确保了字典的可解释性，因为模式 $w_k$ 和样本 $v_n$ 属于相同的空间。
因子 $H$ 的非负性倾向于生成基于局部特征的表示，因为禁止使用减法组合。

具体来说，给定一个$F\times N$的非负矩阵$V$，寻求分解：

$V \simeq WH$

其中，$W \in \mathbb{R}{+}^{F \times K}, H \in \mathbb{R}{+}^{K \times N}$。选择$FK+KN\ll FN$，以实现降维的目的。

NMF的优化

$W$和$H$通常由优化得来，具体来说，我们选择$WH$，以最小化：

$\min _{\mathbf{W}, \mathbf{H} \geq \mathbf{0}} D(\mathbf{V} \mid \mathbf{W H})=\sum_{f n} d\left([\mathbf{V}]_{f n} \mid[\mathbf{W H}]_{f n}\right),$

其中，$d$是一个标量损失函数。

一种常用的选择是$\beta$散度：

$d_\beta(x \mid y) \stackrel{\text { def }}{=}\left\{\begin{array}{cl}\frac{1}{\beta(\beta-1)}\left(x^\beta+(\beta-1) y^\beta-\beta x y^{\beta-1}\right) & \beta \in \mathbb{R} \backslash\{0,1\} \\x \log \frac{x}{y}+(y-x) & \beta=1 \\\frac{x}{y}-\log \frac{x}{y}-1 & \beta=0\end{array}\right.$

散度在$\beta = 0$时是IS散度，$\beta = 1$时是$KL$散度，$\beta = 2$时是欧拉散度。

$\beta$的不同取值对应了不同的散度形式

平方欧几里得距离 / 二次损失（$β=2$）
广义 Kullback-Leibler（KL）散度（$β=1$）
Itakura-Saito（IS）散度（$β=0$）

性质

齐次性：$d\beta(\lambda x | \lambda y)=\lambda^\beta d\beta(x | y)$
$d_{\beta}(x||y)$ 是 $y$ 的凸函数，当 $1≤β≤2$ 时成立

IS散度，是Itakura和Saito基于短时语音谱最大似然估计推导出的，用于度量两个频谱之间的拟合程度的散度，具有良好的听觉效果。IS散度具有以下特点：

尺度不变性（scale invariance）：高低能量的谱分量具有相同的重要性；
特别适合动态范围大的数据，如短时音频谱；
具备明确的统计解释：对应于叠加信号中 W 和 H 的最大似然估计，在 Gamma 噪声下亦可解释；
可采用乘法更新规则（如 Dhillon 和 Sra 的方法），但其理论收敛性仍未完全证明。

乘法更新规则：

Lee和Seung证明，在$\beta = 1,2$时，如以下规则的更新可以使得代价函数下降：
$\begin{aligned}& H \leftarrow H \cdot \frac{W^{\top}\left((W H)^{\beta-2} \odot V\right)}{W^{\top}(W H)^{\beta-1}} \\& W \leftarrow W \cdot \frac{\left((W H)^{\beta-2} \odot V\right) H^{\top}}{(W H)^{\beta-1} H^{\top}}\end{aligned}$
而后对于Bregman散度家族，即通式为：$d_\phi(x \mid y)=\phi(x)-\phi(y)-\nabla \phi(y)(x-y)$的散度家族，Dhillon和Sra提出了以下更新规则：
$\begin{aligned} & H \leftarrow H \cdot \frac{W^{\top}\left(\nabla^2 \phi(W H) \odot V\right)}{W^{\top}\left(\nabla^2 \phi(W H) \odot W H\right)} \\ & W \leftarrow W \cdot \frac{\left(\nabla^2 \phi(W H) \odot V\right) H^{\top}}{\left(\nabla^2 \phi(W H) \odot W H\right) H^{\top}} \end{aligned}$
在实践中，IS散度应用以上规则表现出收敛趋势，但无收敛性证明。

尺度不变性

$\beta$散度满足$d\beta(\gamma x \mid \gamma y)=\gamma^\beta d\beta(x \mid y)$，因而在$\beta = 0$时，散度是尺度不变的，矩阵的高值和低值具有相同的权重，拟合误差不会因为分量大小而有精度上的区别。

IS-NMF的语音建模

在进一步讨论之前，我们来看一个简单的正弦信号的分离问题。待分离信号如下

首先我们计算时频图，如图可见时频图高度冗余。对时频图进行NMF计算，得到的W和H两个矩阵。

从结果来看分离后的$W$和$H$分别包含了两个主模式的频率和时间信息。分别称其为频谱基和时间激活。将$W$和$H$的对应行列再度相乘，得到两个重构的时频图，并所ISTFT可以得到两个分离的信号。

在基础NMF中，接收到的信号被认为是一系列独立随机向量（声源）的线性组合。

这个建模方法过于简单，多数声源的信号不能被简单的表示为一个随机变量，一个随机变量中包含的信息也常不足以重构一个信号。

可以看到上图随着使用的基的数量的增加，重构信号逐渐接近于原始信号。而实际上原始信号仅是双信号的叠加。

接下来关注IS-NMF的语音建模。约定原矩阵为$V$，其中的元素记作$v_{fn}$。分解后，重构矩阵$\hat V = WH$。对于列向量，使用小写加粗字体表示，如$W = [\mathbf w_1, \mathbf w_2,…,\mathbf w_K]$表示；对于行向量，使用普通小写字母表示，如$H = [h_1^T,h_2^T,…,h_K^T]^T$。

IS-NMF 作为高斯分量求和模型的最大似然估计

考虑生成模型：

$\mathbf{x}_n=\sum_{k=1}^K \mathbf{c}_{k, n},$

其中$\mathbf{x}n=\sum{k=1}^K c{k, n},$ 且$c{k, n} \sim \mathcal{N}c\left(0, h{k n} \operatorname{diag}\left(\mathbf{w}k\right)\right),$ 满足多元高斯分布。假设$c{kn} $c_{fn}$ x{fn}$条件独立，且$x{fn}\sim \mathcal{N}c\left(0, h{k n} \operatorname{diag}\left(\mathbf{w}k\right)\right)$。记$v{fn} = |x_{fn}|^2$。则从矩阵$X$最大似然估计得到的$W$和$H$等价于对$V$做以IS散度作为代价函数的非负矩阵分解。

$\begin{aligned}C_{M L, 1}(\mathbf{W}, \mathbf{H}) & =-\sum_{n=1}^N \sum_{f=1}^F \log \mathcal{N}_c\left(x_{f n} \mid 0, \sum_k w_{f k} h_{k n}\right) \\& =N F \log \pi+\sum_{n=1}^N \sum_{f=1}^F \log \left(\sum_k w_{f k} h_{k n}\right)+\frac{\left|x_{f n}\right|^2}{\left(\sum_k w_{f k} h_{k n}\right)} \\& \stackrel{\mathrm{c}}{=} \sum_{n=1}^N \sum_{f=1}^F d_{I S}\left(\left|x_{f n}\right|^2 \mid \sum_k w_{f k} h_{k n}\right)\end{aligned}$

复高斯分布的一维对数密度：
$\log \mathcal{N}_c(x \mid \mu, \lambda)=-\log (\pi \lambda)-\frac{|x-\mu|^2}{\lambda}$

在音频信号的情况下

对于每一个音源（这里假设有两路），各音源的STFT可以表示为：

$s_{1, n}=\sum_{k=1}^{K_1} c_{k, n}, \quad s_{2, n}=\sum_{k=K_1+1}^{K_1+K_2} c_{k, n}, \quad \text { avec } K_1+K_2=K .$

而接受的音频信号的STFT可以被表示为$\mathbf xn =\mathbf s{1,n}+\mathbf s{2,n}$。在$\mathbf x_n = [x{1,n}, x{2,n},…,x{Fn}]$中，帧率列$n = 1,…,N$，频率索引$f = 1,…，F$。与之前的简单建模不同，每一音源的STFT被奸魔为若干元素性分量。每个分量具有功率谱密度（PSD）$\mathbf wk$，并由随帧变化的激活系数$h{kn}$在时间上调制。

传统的方法可能会现在训练数据上学习PSD，在将混合信号的谱图$|X|^2$进行分解到已知的字典：
$W=\left[\mathbf w_1, \ldots, \mathbf w_{K_1}, \mathbf w_{K_1+1}, \ldots, \mathbf w_{K_1+K_2}\right]$
利用最大似然估计和IS-NMF的相关关系，可以不预先进行估计。

定性的，矩阵外积$\mathbf w_k\mathbf h_k^T$在二维坐标中呈现行向量×列向量的结构，而对于时频图：

行向量只随频率变化，因而$\mathbf w_k$将仅能提取频率变化信息，因而能够提取频率模板；

列向量只随时间变化，因而$\mathbf h_k$将仅能提取时间信息，因而能提取时间包络信息

额外的，对于给定$W$和$H$，可以通过维纳滤波得到最小均方误差估计，从而实现对信号的分解：

$\hat{c}_{k, f n}=\frac{w_{f k} h_{k n}}{\sum_{l=1}^K w_{f l} h_{l n}} x_{f n}$

对于固定的$(f,n)$，维纳增益之和为1，分解能量守恒。

考虑${c{k,fn}}$与$x{fn}$的联合分布。两者都复高斯分布，因而联合分布仍未多元复高斯分布：
$\begin{aligned} &\left[\begin{array}{c} c_{f n} \\ x_{f n} \end{array}\right] \sim \mathcal{N}_c(0, \underbrace{\left[\begin{array}{cc} \operatorname{diag}\left(\sigma_{f n}^2\right) & \sigma_{f n}^2 \\ \sigma_{f n}^{22} & \lambda_{f n} \end{array}\right]}_{\Sigma_{f n}}),\\ &\text { avec } \sigma_{f n}^2=\left[\sigma_{1, f n}^2, \ldots, \sigma_{K, f n}^2\right]^{\top}, \lambda_{f n}=\sum_k \sigma_{k, f n}^2 \text {. } \end{aligned}$

对于多元复高斯分布，其条件分布$z|y$仍为复高斯分布，其均值和方差分别为：
$\mathbb{E}[\mathbf{z} \mid y]=\boldsymbol{\Sigma}_{z y} \boldsymbol{\Sigma}_{y y}^{-1} y, \quad \operatorname{Cov}[\mathbf{z} \mid y]=\boldsymbol{\Sigma}_{z z}-\boldsymbol{\Sigma}_{z y} \boldsymbol{\Sigma}_{y y}^{-1} \boldsymbol{\Sigma}_{y z} .$
因而考虑后验分布$c{k,fn}|x{fn}$，有后验均值和后验方差：
$\begin{aligned}&\mathbb E[c_{k,fn}\mid x_{fn} ]=\frac{\sigma_{k, f n}^2}{\lambda_{f n}} x_{f n}=\frac{w_{f k} h_{k n}}{\sum_l w_{f l} h_{l n}} x_{f n}\\& \operatorname{Var}\left[c_{k, f n} \mid x_{f n}\right]=\sigma_{k, f n}^2-\frac{\left(\sigma_{k, f n}^2\right)^2}{\lambda_{f n}}=w_{f k} h_{k n}-\frac{\left(w_{f k} h_{k n}\right)^2}{\sum_l w_{f l} h_{l n}} \end{aligned}$
Wiener滤波方法期望在最小均方误差的准则下，从观测信号中恢复出原始信号。即期望找到一个函数，使得在给定观测$x$的条件下，使得条件均方误差：
$\operatorname{MSE}(x)=\mathbb{E}\left[|c-\hat{c}(x)|^2 \mid x\right]$
最小。使用后验均值$m(x) = \mathbb E[c\mid x]$（条件期望）表示分解MSE。
$E\left[(c-\hat{c}(x))^2 \mid x\right]=E\left[(c-m(x))^2 \mid x\right]+(m(x)-\hat{c}(x))^2,$
第一项是定值，而第二项仅在$\hat c = m$时为零。因此MMSE最小时应有$\mathbb E[c{k,fn}\mid x{fn} ]=\hat{c}_{k, f n}$。

因而得到上式。

由此式实现音源信号的重构。

乘性噪音

对于生成模型，$V = (WH) \odot E$。设$E$时均值为$1$的独立同分布乘性Gamma噪音。则对$W,H$做最大似然估计，等价于将$V$分解为$V = WH$，代价函数采用IS散度。

上定理在此不再证明。若实际噪声确实近似乘性且服从 Gamma（或更一般的、均值为 1 的乘性分布），则 IS-NMF 的代价函数就是该噪声模型下的负对数似然。

多通道非负矩阵分解 MNMF

IS-NMF解决了基础NMF中的语音建模问题，而只能处理单通道数据的问题则再MNMF中得到了解决。。为此，一种解决方法是多通道非负矩阵分解（Multichannel Nonnegative Matrix Factorization，MNMF）。

MNMF的基本原理

考虑第$i\in{1,\ldots,I}$个观测通道的采样信号为：

$\tilde{x}_i(t)=\sum_{j=1}^J \sum_{\tau=0}^{L-1} \tilde{a}_{i j}(\tau) \tilde{s}_j(t-\tau)+\tilde{b}_i(t)$

其中，$\tilde{a}_{i j}(\tau)$是因果的脉冲响应，$\tilde{b}_i(t)$是加性噪音，$\tilde s_j$是第$j$个源的时域信号。即每个麦克风包含的信号不仅包含声音的直接路径，还包含回波信号。在短时傅里叶变换域中表示为：

$x_{i, f n}=\sum_{j=1}^J a_{i j, f} s_{j, f n}+b_{i, f n}\Leftrightarrow \mathbf{X}_{f n}=\mathbf{A}_f \mathbf{S}_{f n}+\mathbf{b}_{f n}$

研究期望将源的功率谱图分解为两个非负矩阵的乘积$|S{fn}|^2\simeq W{fn}H{fn}$，由此系统转化为在给定观察$X{fn}$的情况下，联合估计源谱因子${W{fn},H{fn}}$和混合矩阵$A_{fn}$。

MNMF的音频建模

与IS-NMF相同，MNMF认为一个音源$s{jfn}$由多个表现为复高斯分布的成分$c{kfn}$组成。假设$K\ge J$，且$\left{\mathcal{K}j\right}{j=1}^J$是$\mathcal K = {1,\ldots,K}$的一个非平凡划分。认为第$j$的音源是$# \mathcal K_j$个潜在成分的总和。

$s_{j f n}=\sum_{k \in \mathcal{K}_j} c_{k f n} \quad \text { avec } c_{k f n} \sim \mathcal{N}_{\mathbb{C}}\left(0, w_{f k} h_{k n}\right)$

成分互相独立，因此有：

$s_{j f n} \sim \mathcal{N}_{\mathbb{C}}\left(0, \sum_{k \in \mathcal{K}_j} w_{f k} h_{k n}\right)$

由IS-NMF的证明，给定源$Sj = [s{jfn}]_{fn}$，最大似然估计$W_j$和$H_j$等价于将$|S_j|^2$分解为$W_jH_j$。

$-\log p\left(\mathbf{S}_j \mid \mathbf{W}_j, \mathbf{H}_j\right) \cong \sum_{f n} d_{I S}\left(\left|s_{j f n}\right|^2 \|\left[\mathbf{W}_j \mathbf{H}_j\right]_{f n}\right)$

由于对成分的假设，成分是“有限”且“可重复的”，即不同音源可能包含相同的成分，而总的成分数量小于$K$。因此，可以将混合模型重新写作：

$\mathbf{x}_{f n}=\mathbf{A}_f \mathbf{c}_{f n}+\mathbf{b}_{f n}$

其中，$\mathbf{c}{f n}=\left[c{1, f n}, \ldots, c{K, f n}\right]^{\top} \in \mathbb{C}^{K \times 1}$为有限种成分。因此，模型被简化作一个有$I$个通道和$K$个高思远的线性混合模型。记$\boldsymbol{\Sigma}{c, f n}=\operatorname{diag}\left(\left[w{f k} h{k n}\right]\right)$。

噪音被假设为平稳且空间上不相关的：

$b_{i f n} \sim \mathcal{N}_{\mathbb{C}}\left(0, \sigma_{i, f}^2\right)$

噪音成分除了模拟常规意义上的噪音之外，同样有利于解释模型不匹配和量化带来的“噪音”。因此引入此噪音对避免数据不稳定性和加速收敛是有必要。记$\boldsymbol{\Sigma}{b, f}=\operatorname{diag}\left(\left[\sigma{i, f}^2\right]_i\right)$。

NMNF的优化方法

在当前建模下，所有参数包括$\boldsymbol{\theta}=\left{\mathbf{A}, \mathbf{W}, \mathbf{H}, \boldsymbol{\Sigma}b\right}$。根据之前的假设， $x{fn}$是一个零均值的复高斯变量，协方差为：

$\boldsymbol{\Sigma}_{x, f n}(\boldsymbol{\theta})=\mathbf{A}_f \boldsymbol{\Sigma}_{s, f n} \mathbf{A}_f^H+\boldsymbol{\Sigma}_{b, f}$

$\Sigma_{s,fn}$是源的协方差矩阵。此时，最大似然估计等价于最小化：

$\mathcal{C}_1(\boldsymbol{\theta})=\sum_{f n} \operatorname{trace}\left(\mathbf{x}_{f n} \mathbf{x}_{f n}^H \boldsymbol{\Sigma}_{x, f n}^{-1}\right)+\log \operatorname{det} \boldsymbol{\Sigma}_{x, f n}$

$x_{fn}$条件概率密度函数可以写作：
$p\left(x_{f n} \mid \theta\right)=\frac{1}{\pi^I \operatorname{det}\left(\Sigma_{x, f n}\right)} \exp \left(-x_{f n}^{\mathrm{H}} \Sigma_{x, f n}^{-1} x_{f n}\right) .$
负对数似然即为：
$\mathcal{L}(\theta)=\sum_{f, n}\left[x_{f n}^{\mathrm{H}} \Sigma_{x, f n}^{-1} x_{f n}+\log \operatorname{det}\left(\Sigma_{x, f n}\right)+I \log \pi\right]$
最后一项为常数，可以忽略。将二次型写作迹形式，即可得到上式。

在优化中，通常施加$\sumi\left|a{i j, f}\right|^2=1$，且$\sumj w{f k}=1$的约束以减少不确定性。这些约束使得$A_f$的列对通道归一化，$W$的列对频率归一化，幅值信息保留在$H$中随时间变化。