引言
在现实世界中,我们常常遇到多种信号混合在一起的情况。盲源分离(Blind Source Separation,BSS)旨在从混合信号中恢复出原始信号,而不依赖于对信号源或混合过程的先验知识。
问题的数学表述
基本模型
假设我们有 n n n 个源信号 s ( t ) = [ s 1 ( t ) , s 2 ( t ) , … , s n ( t ) ] T ∈ R n \mathbf{s}(t) = [s_1(t), s_2(t), \ldots, s_n(t)]^T \in \mathbb{R}^n s(t)=[s1(t),s2(t),…,sn(t)]T∈Rn,它们经过某种未知的混合矩阵 A ∈ R m × n \mathbf{A} \in \mathbb{R}^{m \times n} A∈Rm×n 混合后,得到 m m m 个观测信号 x ( t ) = [ x 1 ( t ) , x 2 ( t ) , … , x m ( t ) ] T ∈ R m \mathbf{x}(t) = [x_1(t), x_2(t), \ldots, x_m(t)]^T \in \mathbb{R}^m x(t)=[x1(t),x2(t),…,xm(t)]T∈Rm。考虑加性噪声 n ( t ) ∈ R m \mathbf{n}(t) \in \mathbb{R}^m n(t)∈Rm,混合模型可以表示为:
x ( t ) = A s ( t ) + n ( t ) \mathbf{x}(t) = \mathbf{A}\mathbf{s}(t) + \mathbf{n}(t) x(t)=As(t)+n(t)
在频域中,混合模型可以表示为:
X ( ω ) = A ( ω ) S ( ω ) + N ( ω ) \mathbf{X}(\omega) = \mathbf{A}(\omega)\mathbf{S}(\omega) + \mathbf{N}(\omega) X(ω)=A(ω)S(ω)+N(ω)
其中 X ( ω ) \mathbf{X}(\omega) X(ω)、 S ( ω ) \mathbf{S}(\omega) S(ω) 和 N ( ω ) \mathbf{N}(\omega) N(ω) 分别是 x ( t ) \mathbf{x}(t) x(t)、 s ( t ) \mathbf{s}(t) s(t) 和 n ( t ) \mathbf{n}(t) n(t) 的傅里叶变换。
盲源分离的目标是找到一个分离矩阵 W ∈ R n × m \mathbf{W} \in \mathbb{R}^{n \times m} W∈Rn×m,使得:
y ( t ) = W x ( t ) ≈ P D s ( t ) \mathbf{y}(t) = \mathbf{W}\mathbf{x}(t) \approx \mathbf{P}\mathbf{D}\mathbf{s}(t) y(t)=Wx(t)≈PDs(t)
其中 y ( t ) ∈ R n \mathbf{y}(t) \in \mathbb{R}^n y(t)∈Rn 是对源信号的估计, P ∈ R n × n \mathbf{P} \in \mathbb{R}^{n \times n} P∈Rn×n 是置换矩阵, D ∈ R n × n \mathbf{D} \in \mathbb{R}^{n \times n} D∈Rn×n 是对角矩阵。
混合模型的数学分类
1. 瞬时线性混合
瞬时线性混合是最简单的混合模型,假设混合过程是线性的且没有时间延迟:
x i ( t ) = ∑ j = 1 n a i j s j ( t ) , i = 1 , 2 , … , m x_i(t) = \sum_{j=1}^{n} a_{ij}s_j(t), \quad i = 1,2,\ldots,m xi(t)=j=1∑naijsj(t),i=1,2,…,m
这可以表示为矩阵形式: x ( t ) = A s ( t ) \mathbf{x}(t) = \mathbf{A}\mathbf{s}(t) x(t)=As(t),其中 A ∈ R m × n \mathbf{A} \in \mathbb{R}^{m \times n} A∈Rm×n 是常数矩阵。
2. 卷积混合
卷积混合考虑了信号在传播过程中的时间延迟和多径效应:
x i ( t ) = ∑ j = 1 n ∑ k = 0 L − 1 a i j ( k ) s j ( t − k ) , i = 1 , 2 , … , m x_i(t) = \sum_{j=1}^{n} \sum_{k=0}^{L-1} a_{ij}(k)s_j(t-k), \quad i = 1,2,\ldots,m xi(t)=j=1∑nk=0∑L−1aij(k)sj(t−k),i=1,2,…,m
其中 a i j ( k ) a_{ij}(k) aij(k) 是从源 j j j 到传感器 i i i 的有限冲激响应(FIR)滤波器的第 k k k 个系数, L L L 是滤波器长度。
用Z变换表示为:
X ( z ) = A ( z ) S ( z ) \mathbf{X}(z) = \mathbf{A}(z)\mathbf{S}(z) X(z)=A(z)S(z)
其中 A ( z ) ∈ C m × n \mathbf{A}(z) \in \mathbb{C}^{m \times n} A(z)∈Cm×n 是混合矩阵在Z域中的表示:
A ( z ) = ∑ k = 0 L − 1 A k z − k \mathbf{A}(z) = \sum_{k=0}^{L-1} \mathbf{A}_k z^{-k} A(z)=k=0∑L−1Akz−k
对应的分离系统为:
Y ( z ) = W ( z ) X ( z ) \mathbf{Y}(z) = \mathbf{W}(z)\mathbf{X}(z) Y(z)=W(z)X(z)
其中 W ( z ) ∈ C n × m \mathbf{W}(z) \in \mathbb{C}^{n \times m} W(z)∈Cn×m 是分离矩阵在Z域中的表示。
3. 非线性混合
非线性混合假设信号之间的混合是非线性的:
x ( t ) = f ( s ( t ) ) \mathbf{x}(t) = \mathbf{f}(\mathbf{s}(t)) x(t)=f(s(t))
其中 f : R n → R m \mathbf{f}: \mathbb{R}^n \rightarrow \mathbb{R}^m f:Rn→Rm 是一个非线性函数。更具体地,可以表示为:
x i ( t ) = f i ( s 1 ( t ) , s 2 ( t ) , … , s n ( t ) ) + n i ( t ) , i = 1 , 2 , … , m x_i(t) = f_i(s_1(t), s_2(t), \ldots, s_n(t)) + n_i(t), \quad i = 1,2,\ldots,m xi(t)=fi(s1(t),s2(t),…,sn(t))+ni(t),i=1,2,…,m
针对非线性混合,我们通常使用泛函级数展开来近似非线性函数:
f i ( s ) = ∑ j = 1 n a i j ( 1 ) s j + ∑ j 1 = 1 n ∑ j 2 = 1 n a i j 1 j 2 ( 2 ) s j 1 s j 2 + ∑ j 1 = 1 n ∑ j 2 = 1 n ∑ j 3 = 1 n a i j 1 j 2 j 3 ( 3 ) s j 1 s j 2 s j 3 + ⋯ f_i(\mathbf{s}) = \sum_{j=1}^{n} a_{ij}^{(1)}s_j + \sum_{j_1=1}^{n}\sum_{j_2=1}^{n} a_{ij_1j_2}^{(2)}s_{j_1}s_{j_2} + \sum_{j_1=1}^{n}\sum_{j_2=1}^{n}\sum_{j_3=1}^{n} a_{ij_1j_2j_3}^{(3)}s_{j_1}s_{j_2}s_{j_3} + \cdots fi(s)=j=1∑naij(1)sj+j1=1∑nj2=1∑naij1j2(2)sj1sj2+j1=1∑nj2=1∑nj3=1∑naij1j2j3(3)sj1sj2sj3+⋯
其中 a i j 1 j 2 … j r ( r ) a_{ij_1j_2\ldots j_r}^{(r)} aij1j2…jr(r) 是 r r r 阶Volterra核的系数。
主要分离理论与算法
独立成分分析(ICA)
独立成分分析是最广泛使用的盲源分离方法,基于以下关键假设:
- 源信号之间是统计独立的
- 最多只有一个源信号服从高斯分布
- 混合矩阵是满秩的
数学原理
ICA的目标是找到一个变换 W \mathbf{W} W,使得变换后的信号 y = W x \mathbf{y} = \mathbf{W}\mathbf{x} y=Wx 中的各个分量尽可能独立。
两个随机变量 y 1 y_1 y1 和 y 2 y_2 y2 的独立性可以通过它们的联合概率密度函数与边缘概率密度函数的关系来表示:
p ( y 1 , y 2 ) = p ( y 1 ) p ( y 2 ) p(y_1, y_2) = p(y_1)p(y_2) p(y1,y2)=p(y1)p(y2)
对于多个随机变量,其互信息(mutual information)定义为:
I ( y 1 , y 2 , … , y n ) = ∑ i = 1 n H ( y i ) − H ( y 1 , y 2 , … , y n ) I(y_1, y_2, \ldots, y_n) = \sum_{i=1}^{n} H(y_i) - H(y_1, y_2, \ldots, y_n) I(y1,y2,…,yn)=i=1∑nH(yi)−H(y1,y2,…,yn)
其中 H ( y ) H(y) H(y) 表示微分熵:
H ( y ) = − ∫ p ( y ) log p ( y ) d y H(y) = -\int p(y) \log p(y) dy H(y)=−∫p(y)logp(y)dy
互信息等于零当且仅当所有随机变量相互独立。因此,ICA的目标可以表示为最小化互信息:
min W I ( y 1 , y 2 , … , y n ) \min_{\mathbf{W}} I(y_1, y_2, \ldots, y_n) WminI(y1,y2,…,yn)
由于直接计算互信息比较困难,实际算法常使用其他度量,如非高斯性、负熵等。
非高斯性与负熵
根据中心极限定理,多个独立随机变量的线性组合会比原始变量更接近高斯分布。因此,最大化非高斯性可以用来分离独立成分。
峰度(kurtosis)是衡量非高斯性的一种方式,定义为:
kurt ( y ) = E { y 4 } − 3 ( E { y 2 } ) 2 \operatorname{kurt}(y) = E\{y^4\} - 3(E\{y^2\})^2 kurt(y)=E{y4}−3(E{y2})2
超高斯分布(如拉普拉斯分布)的峰度为正,而次高斯分布(如均匀分布)的峰度为负。
负熵(negentropy)是另一种更稳健的非高斯性度量,定义为:
J ( y ) = H ( y gauss ) − H ( y ) J(y) = H(y_{\text{gauss}}) - H(y) J(y)=H(ygauss)−H(y)
其中 y gauss y_{\text{gauss}} ygauss 是与 y y y 具有相同方差的高斯随机变量。
负熵的近似计算公式为:
J ( y ) ≈ ∑ i = 1 p k i [ E { G i ( y ) } − E { G i ( ν ) } ] 2 J(y) \approx \sum_{i=1}^{p} k_i[E\{G_i(y)\} - E\{G_i(\nu)\}]^2 J(y)≈i=1∑pki[E{Gi(y)}−E{Gi(ν)}]2
其中 G i G_i Gi 是非二次非线性函数,常用的选择有:
G 1 ( u ) = 1 a 1 log cosh ( a 1 u ) , G 2 ( u ) = − exp ( − u 2 / 2 ) G_1(u) = \frac{1}{a_1} \log \cosh(a_1 u), \quad G_2(u) = -\exp(-u^2/2) G1(u)=a11logcosh(a1u),G2(u)=−exp(−u2/2)
其中 1 ≤ a 1 ≤ 2 1 \leq a_1 \leq 2 1≤a1≤2 是一个常数, ν \nu ν 是均值为0、方差为1的高斯随机变量, k i k_i ki 是正常数。
FastICA算法
FastICA是一种快速高效的ICA实现算法,基于固定点迭代方法。其目标函数是最大化每个分离成分的非高斯性,同时保持各个分离成分之间的正交性。
预处理步骤:
- 中心化数据: x ← x − E [ x ] \mathbf{x} \leftarrow \mathbf{x} - E[\mathbf{x}] x←x−E[x]
- 白化:首先计算协方差矩阵 C = E [ x x T ] \mathbf{C} = E[\mathbf{xx}^T] C=E[xxT],进行特征值分解得到 C = E D E T \mathbf{C} = \mathbf{E}\mathbf{D}\mathbf{E}^T C=EDET,然后计算白化矩阵 V = D − 1 / 2 E T \mathbf{V} = \mathbf{D}^{-1/2}\mathbf{E}^T V=D−1/2ET,白化后的数据为 z = V x \mathbf{z} = \mathbf{V}\mathbf{x} z=Vx,满足 E [ z z T ] = I E[\mathbf{zz}^T] = \mathbf{I} E[zzT]=I
对每个成分 w i \mathbf{w}_i wi 的迭代方法:
- 初始化 w i \mathbf{w}_i wi 为单位范数的随机向量
- 更新: w i + ← E [ z g ( w i T z ) ] − E [ g ′ ( w i T z ) ] w i \mathbf{w}_i^+ \leftarrow E[\mathbf{z}g(\mathbf{w}_i^T\mathbf{z})] - E[g'(\mathbf{w}_i^T\mathbf{z})]\mathbf{w}_i wi+←E[zg(wiTz)]−E[g′(wiTz)]wi
- 单位化: w i ← w i + / ∥ w i + ∥ \mathbf{w}_i \leftarrow \mathbf{w}_i^+ / \|\mathbf{w}_i^+\| wi←wi+/∥wi+∥
- 如果 w i \mathbf{w}_i wi 未收敛,返回步骤2
- 对于 j < i j < i j<i 的所有 j j j,进行偏移正交化: w i ← w i − ∑ j = 1 i − 1 ( w i T w j ) w j \mathbf{w}_i \leftarrow \mathbf{w}_i - \sum_{j=1}^{i-1}(\mathbf{w}_i^T\mathbf{w}_j)\mathbf{w}_j wi←wi−∑j=1i−1(wiTwj)wj
- 归一化: w i ← w i / ∥ w i ∥ \mathbf{w}_i \leftarrow \mathbf{w}_i / \|\mathbf{w}_i\| wi←wi/∥wi∥
其中 g ( u ) g(u) g(u) 是 G ( u ) G(u) G(u) 的导数, g ′ ( u ) g'(u) g′(u) 是 g ( u ) g(u) g(u) 的导数。常用的 g ( u ) g(u) g(u) 函数有:
g 1 ( u ) = tanh ( a 1 u ) , g 2 ( u ) = u exp ( − u 2 / 2 ) g_1(u) = \tanh(a_1 u), \quad g_2(u) = u \exp(-u^2/2) g1(u)=tanh(a1u),g2(u)=uexp(−u2/2)
自然梯度算法
自然梯度算法是另一种常用的ICA算法,它考虑了参数空间的黎曼几何结构。自然梯度的更新规则为:
Δ W ∝ η [ ( I − E { g ( y ) y T } ) W ] \Delta \mathbf{W} \propto \eta [(\mathbf{I} - E\{\mathbf{g}(\mathbf{y})\mathbf{y}^T\})\mathbf{W}] ΔW∝η[(I−E{g(y)yT})W]
其中 g ( y ) = [ g ( y 1 ) , g ( y 2 ) , … , g ( y n ) ] T \mathbf{g}(\mathbf{y}) = [g(y_1), g(y_2), \ldots, g(y_n)]^T g(y)=[g(y1),g(y2),…,g(yn)]T, η \eta η 是学习率。
主成分分析(PCA)
主成分分析是一种广泛使用的降维方法,也常作为ICA的预处理步骤。PCA通过寻找数据的最大方差方向来降维。
数学原理
给定数据矩阵 X ∈ R m × N \mathbf{X} \in \mathbb{R}^{m \times N} X∈Rm×N,其中 m m m 是特征维度, N N N 是样本数量。首先将数据中心化: X ← X − μ 1 T \mathbf{X} \leftarrow \mathbf{X} - \mathbf{\mu}\mathbf{1}^T X←X−μ1T,其中 μ = 1 N X 1 \mathbf{\mu} = \frac{1}{N}\mathbf{X}\mathbf{1} μ=N1X1 是均值向量, 1 \mathbf{1} 1 是全1向量。
PCA的目标是找到一组正交单位向量 { u 1 , u 2 , … , u m } \{\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_m\} {u1,u2,…,um},使得投影后的方差最大:
u i = arg max ∥ u ∥ = 1 , u ⊥ u j , j < i 1 N ∑ j = 1 N ( u T x j ) 2 = arg max ∥ u ∥ = 1 , u ⊥ u j , j < i u T C u \mathbf{u}_i = \arg\max_{\|\mathbf{u}\|=1, \mathbf{u} \perp \mathbf{u}_j, j<i} \frac{1}{N}\sum_{j=1}^{N}(\mathbf{u}^T\mathbf{x}_j)^2 = \arg\max_{\|\mathbf{u}\|=1, \mathbf{u} \perp \mathbf{u}_j, j<i} \mathbf{u}^T\mathbf{C}\mathbf{u} ui=arg∥u∥=1,u⊥uj,j<imaxN1j=1∑N(uTxj)2=arg∥u∥=1,u⊥uj,j<imaxuTCu
其中 C = 1 N X X T \mathbf{C} = \frac{1}{N}\mathbf{X}\mathbf{X}^T C=N1XXT 是数据的协方差矩阵。
这个问题的解可以通过对协方差矩阵 C \mathbf{C} C 进行特征值分解得到:
C = U Λ U T \mathbf{C} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^T C=UΛUT
其中 U = [ u 1 , u 2 , … , u m ] \mathbf{U} = [\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_m] U=[u1,u2,…,um] 是特征向量矩阵, Λ = diag ( λ 1 , λ 2 , … , λ m ) \mathbf{\Lambda} = \text{diag}(\lambda_1, \lambda_2, \ldots, \lambda_m) Λ=diag(λ1,λ2,…,λm) 是特征值对角矩阵,且 λ 1 ≥ λ 2 ≥ … ≥ λ m \lambda_1 \geq \lambda_2 \geq \ldots \geq \lambda_m λ1≥λ2≥…≥λm。
PCA降维到 r r r 维可以通过保留前 r r r 个主成分实现:
Y = U r T X \mathbf{Y} = \mathbf{U}_r^T\mathbf{X} Y=UrTX
其中 U r = [ u 1 , u 2 , … , u r ] \mathbf{U}_r = [\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_r] Ur=[u1,u2,…,ur], Y ∈ R r × N \mathbf{Y} \in \mathbb{R}^{r \times N} Y∈Rr×N 是降维后的数据。
PCA与白化
白化是ICA的重要预处理步骤,其目的是将数据转换为零均值、单位协方差矩阵的形式。白化变换可以通过PCA实现:
Z = Λ − 1 / 2 U T X \mathbf{Z} = \mathbf{\Lambda}^{-1/2}\mathbf{U}^T\mathbf{X} Z=Λ−1/2UTX
其中 Z ∈ R m × N \mathbf{Z} \in \mathbb{R}^{m \times N} Z∈Rm×N 是白化后的数据,满足 1 N Z Z T = I \frac{1}{N}\mathbf{Z}\mathbf{Z}^T = \mathbf{I} N1ZZT=I。
白化的优点在于:
- 降低数据的复杂性,简化后续ICA算法
- 消除数据之间的二阶相关性
- 减少估计参数的数量,从 m n mn mn 减少到 n ( n − 1 ) / 2 n(n-1)/2 n(n−1)/2
非负矩阵分解(NMF)
非负矩阵分解适用于数据为非负的情况,如图像、文本词频、音频功率谱等。NMF的目标是将非负矩阵 V ∈ R + m × n \mathbf{V} \in \mathbb{R}_{+}^{m \times n} V∈R+m×n 分解为两个非负矩阵的乘积:
V ≈ W H \mathbf{V} \approx \mathbf{W}\mathbf{H} V≈WH
其中 W ∈ R + m × r \mathbf{W} \in \mathbb{R}_{+}^{m \times r} W∈R+m×r, H ∈ R + r × n \mathbf{H} \in \mathbb{R}_{+}^{r \times n} H∈R+r×n,且 r ≤ min ( m , n ) r \leq \min(m, n) r≤min(m,n)。
目标函数
NMF常用的目标函数有:
- Frobenius范数(对应欧氏距离):
D F ( V ∣ ∣ W H ) = 1 2 ∥ V − W H ∥ F 2 = 1 2 ∑ i = 1 m ∑ j = 1 n ( V i j − ( W H ) i j ) 2 D_F(\mathbf{V}||\mathbf{W}\mathbf{H}) = \frac{1}{2}\|\mathbf{V} - \mathbf{W}\mathbf{H}\|_F^2 = \frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{n}(V_{ij} - (\mathbf{W}\mathbf{H})_{ij})^2 DF(V∣∣WH)=21∥V−WH∥F2=21i=1∑mj=1∑n(Vij−(WH)ij)2
- Kullback-Leibler散度(当数据可以被解释为概率分布时):
D K L ( V ∣ ∣ W H ) = ∑ i = 1 m ∑ j = 1 n ( V i j log V i j ( W H ) i j − V i j + ( W H ) i j ) D_{KL}(\mathbf{V}||\mathbf{W}\mathbf{H}) = \sum_{i=1}^{m}\sum_{j=1}^{n}\left(V_{ij}\log\frac{V_{ij}}{(\mathbf{W}\mathbf{H})_{ij}} - V_{ij} + (\mathbf{W}\mathbf{H})_{ij}\right) DKL(V∣∣WH)=i=1∑mj=1∑n(Vijlog(WH)ijVij−Vij+(WH)ij)
- Itakura-Saito散度(适用于音频功率谱):
D I S ( V ∣ ∣ W H ) = ∑ i = 1 m ∑ j = 1 n ( V i j ( W H ) i j − log V i j ( W H ) i j − 1 ) D_{IS}(\mathbf{V}||\mathbf{W}\mathbf{H}) = \sum_{i=1}^{m}\sum_{j=1}^{n}\left(\frac{V_{ij}}{(\mathbf{W}\mathbf{H})_{ij}} - \log\frac{V_{ij}}{(\mathbf{W}\mathbf{H})_{ij}} - 1\right) DIS(V∣∣WH)=i=1∑mj=1∑n((WH)ijVij−log(WH)ijVij−1)
- β-散度(统一了上述三种散度):
D β ( V ∣ ∣ W H ) = 1 β ( β − 1 ) ∑ i = 1 m ∑ j = 1 n ( V i j β − ( W H ) i j β − β V i j β − 1 ( V i j − ( W H ) i j ) ) D_{\beta}(\mathbf{V}||\mathbf{W}\mathbf{H}) = \frac{1}{\beta(\beta-1)}\sum_{i=1}^{m}\sum_{j=1}^{n}\left(V_{ij}^{\beta} - (\mathbf{W}\mathbf{H})_{ij}^{\beta} - \beta V_{ij}^{\beta-1}(V_{ij} - (\mathbf{W}\mathbf{H})_{ij})\right) Dβ(V∣∣WH)=β(β−1)1i=1∑mj=1∑n(Vijβ−(WH)ijβ−βVijβ−1(Vij−(WH)ij))
其中 β = 2 \beta = 2 β=2 对应Frobenius范数, β → 1 \beta \rightarrow 1 β→1 对应KL散度, β → 0 \beta \rightarrow 0 β→0 对应IS散度。
优化算法
NMF的优化算法主要有两类:乘性更新规则和交替非负最小二乘法。
乘性更新规则(以Frobenius范数为例):
W i j ← W i j ( V H T ) i j ( W H H T ) i j W_{ij} \leftarrow W_{ij} \frac{(\mathbf{V}\mathbf{H}^T)_{ij}}{(\mathbf{W}\mathbf{H}\mathbf{H}^T)_{ij}} Wij←Wij(WHHT)ij(VHT)ij
H i j ← H i j ( W T V ) i j ( W T W H ) i j H_{ij} \leftarrow H_{ij} \frac{(\mathbf{W}^T\mathbf{V})_{ij}}{(\mathbf{W}^T\mathbf{W}\mathbf{H})_{ij}} Hij←Hij(WTWH)ij(WTV)ij
对于KL散度,更新规则为:
W i j ← W i j ∑ k = 1 n V i k H j k / ( W H ) i k ∑ k = 1 n H j k W_{ij} \leftarrow W_{ij} \frac{\sum_{k=1}^{n} V_{ik}H_{jk}/(\mathbf{W}\mathbf{H})_{ik}}{\sum_{k=1}^{n} H_{jk}} Wij←Wij∑k=1nHjk∑k=1nVikHjk/(WH)ik
H i j ← H i j ∑ k = 1 m W k i V k j / ( W H ) k j ∑ k = 1 m W k i H_{ij} \leftarrow H_{ij} \frac{\sum_{k=1}^{m} W_{ki}V_{kj}/(\mathbf{W}\mathbf{H})_{kj}}{\sum_{k=1}^{m} W_{ki}} Hij←Hij∑k=1mWki∑k=1mWkiVkj/(WH)kj
稀疏编码(Sparse Coding)
稀疏编码假设源信号可以用少量基函数的线性组合表示,即信号在某个基或字典下是稀疏的。数学表示为:
x = D s + n \mathbf{x} = \mathbf{D}\mathbf{s} + \mathbf{n} x=Ds+n
其中 D ∈ R m × p \mathbf{D} \in \mathbb{R}^{m \times p} D∈Rm×p 是字典或基矩阵, s ∈ R p \mathbf{s} \in \mathbb{R}^p s∈Rp 是稀疏系数向量(绝大多数元素为零), n \mathbf{n} n 是噪声。
数学优化问题
稀疏编码通常表述为以下优化问题:
- 基于l0范数(计数稀疏度):
min s ∥ x − D s ∥ 2 2 s.t. ∥ s ∥ 0 ≤ K \min_{\mathbf{s}} \|\mathbf{x} - \mathbf{D}\mathbf{s}\|_2^2 \quad \text{s.t.} \quad \|\mathbf{s}\|_0 \leq K smin∥x−Ds∥22s.t.∥s∥0≤K
其中 ∥ s ∥ 0 \|\mathbf{s}\|_0 ∥s∥0 表示向量 s \mathbf{s} s 中非零元素的个数, K K K 是稀疏度参数。
- 基于l1范数(凸松弛):
min s ∥ x − D s ∥ 2 2 + λ ∥ s ∥ 1 \min_{\mathbf{s}} \|\mathbf{x} - \mathbf{D}\mathbf{s}\|_2^2 + \lambda \|\mathbf{s}\|_1 smin∥x−Ds∥22+λ∥s∥1
其中 ∥ s ∥ 1 = ∑ i = 1 p ∣ s i ∣ \|\mathbf{s}\|_1 = \sum_{i=1}^{p} |s_i| ∥s∥1=∑i=1p∣si∣ 是l1范数, λ \lambda λ 是正则化参数。
- 字典学习问题:
min D , S ∥ X − D S ∥ F 2 + λ ∥ S ∥ 1 s.t. ∥ d j ∥ 2 = 1 , ∀ j \min_{\mathbf{D}, \mathbf{S}} \|\mathbf{X} - \mathbf{DS}\|_F^2 + \lambda \|\mathbf{S}\|_1 \quad \text{s.t.} \quad \|\mathbf{d}_j\|_2 = 1, \forall j D,Smin∥X−DS∥F2+λ∥S∥1s.t.∥dj∥2=1,∀j
其中 X = [ x 1 , x 2 , … , x n ] ∈ R m × n \mathbf{X} = [\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_n] \in \mathbb{R}^{m \times n} X=[x1,x2,…,xn]∈Rm×n 是数据矩阵, S = [ s 1 , s 2 , … , s n ] ∈ R p × n \mathbf{S} = [\mathbf{s}_1, \mathbf{s}_2, \ldots, \mathbf{s}_n] \in \mathbb{R}^{p \times n} S=[s1,s2,…,sn]∈Rp×n 是稀疏系数矩阵, d j \mathbf{d}_j dj 是字典 D \mathbf{D} D 的第 j j j 列。
算法实现
稀疏编码的常用算法包括:
-
匹配追踪(Matching Pursuit,MP):
MP是一种贪婪算法,它在每一步选择与残差最相关的字典原子:
j ∗ = arg max j ∣ ⟨ r ( i ) , d j ⟩ ∣ j^* = \arg\max_j |\langle \mathbf{r}^{(i)}, \mathbf{d}_j \rangle| j∗=argjmax∣⟨r(i),dj⟩∣
其中 r ( i ) \mathbf{r}^{(i)} r(i) 是第 i i i 步的残差, d j \mathbf{d}_j dj 是字典的第 j j j 列。
然后更新系数和残差:
s j ∗ ← s j ∗ + ⟨ r ( i ) , d j ∗ ⟩ s_{j^*} \leftarrow s_{j^*} + \langle \mathbf{r}^{(i)}, \mathbf{d}_{j^*} \rangle sj∗←sj∗+⟨r(i),dj∗⟩
r ( i + 1 ) ← r ( i ) − ⟨ r ( i ) , d j ∗ ⟩ d j ∗ \mathbf{r}^{(i+1)} \leftarrow \mathbf{r}^{(i)} - \langle \mathbf{r}^{(i)}, \mathbf{d}_{j^*} \rangle \mathbf{d}_{j^*} r(i+1)←r(i)−⟨r(i),dj∗⟩dj∗
-
正交匹配追踪(Orthogonal Matching Pursuit,OMP):
OMP在每一步选择与残差最相关的字典原子,但与MP不同,它在每一步都重新计算所有已选择原子的系数:
s Λ ( i ) = arg min s Λ ( i ) ∥ x − D Λ ( i ) s Λ ( i ) ∥ 2 2 \mathbf{s}_{\Lambda^{(i)}} = \arg\min_{\mathbf{s}_{\Lambda^{(i)}}} \|\mathbf{x} - \mathbf{D}_{\Lambda^{(i)}}\mathbf{s}_{\Lambda^{(i)}}\|_2^2 sΛ(i)=argsΛ(i)min∥x−DΛ(i)sΛ(i)∥22
其中 Λ ( i ) \Lambda^{(i)} Λ(i) 是截至第 i i i 步已选择的原子指标集合, D Λ ( i ) \mathbf{D}_{\Lambda^{(i)}} DΛ(i) 是由这些原子组成的子字典。
-
LASSO(Least Absolute Shrinkage and Selection Operator):
LASSO通过求解l1正则化问题来获得稀疏解:
min s 1 2 ∥ x − D s ∥ 2 2 + λ ∥ s ∥ 1 \min_{\mathbf{s}} \frac{1}{2}\|\mathbf{x} - \mathbf{D}\mathbf{s}\|_2^2 + \lambda \|\mathbf{s}\|_1 smin21∥x−Ds∥22+λ∥s∥1
常用的算法有坐标下降法、LARS(Least Angle Regression)、近似消息传递(AMP)等。
盲源分离的挑战与局限性
数学不确定性问题
盲源分离存在固有的数学不确定性,主要包括:
1. 排列不确定性
分离出的信号顺序与原始信号可能不同,可以表示为一个置换矩阵 P \mathbf{P} P:
y = P s \mathbf{y} = \mathbf{P}\mathbf{s} y=Ps
其中 P \mathbf{P} P 是 n × n n \times n n×n 的置换矩阵,仅在每行和每列有一个元素为1,其余为0。
2. 尺度不确定性
分离出的信号幅度可能与原始信号不同,可以表示为一个对角矩阵 D \mathbf{D} D:
y = D s \mathbf{y} = \mathbf{D}\mathbf{s} y=Ds
其中 D = diag ( d 1 , d 2 , … , d n ) \mathbf{D} = \text{diag}(d_1, d_2, \ldots, d_n) D=diag(d1,d2,…,dn) 是对角矩阵, d i ≠ 0 d_i \neq 0 di=0。
综合考虑这两种不确定性,分离结果可以表示为:
y = P D s \mathbf{y} = \mathbf{P}\mathbf{D}\mathbf{s} y=PDs
在盲源分离问题中,我们只能确保 W A = P D \mathbf{W}\mathbf{A} = \mathbf{P}\mathbf{D} WA=PD 而非单位矩阵 I \mathbf{I} I。
欠定与过定问题
从系统可解性的角度,盲源分离问题可分为三种情况:
1. 正定系统( m = n m = n m=n)
当观测信号数量等于源信号数量时,系统是正定的。理论上,如果混合矩阵 A \mathbf{A} A 是可逆的,那么存在唯一的分离矩阵 W = A − 1 \mathbf{W} = \mathbf{A}^{-1} W=A−1(忽略排列和尺度不确定性)。
2. 欠定系统( m < n m < n m<n)
当观测信号数量少于源信号数量时,系统是欠定的。在这种情况下,线性代数理论告诉我们,方程组 x = A s \mathbf{x} = \mathbf{A}\mathbf{s} x=As 有无穷多解。
为了使问题可解,需要引入额外的约束或先验信息,如信号的稀疏性。此时,可以将问题表述为:
min s ∥ s ∥ 0 s.t. x = A s \min_{\mathbf{s}} \|\mathbf{s}\|_0 \quad \text{s.t.} \quad \mathbf{x} = \mathbf{A}\mathbf{s} smin∥s∥0s.t.x=As
或者使用l1范数的凸松弛:
min s ∥ s ∥ 1 s.t. x = A s \min_{\mathbf{s}} \|\mathbf{s}\|_1 \quad \text{s.t.} \quad \mathbf{x} = \mathbf{A}\mathbf{s} smin∥s∥1s.t.x=As
3. 过定系统( m > n m > n m>n)
当观测信号数量多于源信号数量时,系统是过定的。这种情况下,可以利用额外的观测信号来提高分离的鲁棒性。
过定系统中,由于噪声的存在,方程 x = A s \mathbf{x} = \mathbf{A}\mathbf{s} x=As 通常无法精确满足。此时,可以采用最小二乘估计:
s ^ = ( A T A ) − 1 A T x \hat{\mathbf{s}} = (\mathbf{A}^T\mathbf{A})^{-1}\mathbf{A}^T\mathbf{x} s^=(ATA)−1ATx
或者正则化最小二乘:
s ^ = ( A T A + λ I ) − 1 A T x \hat{\mathbf{s}} = (\mathbf{A}^T\mathbf{A} + \lambda\mathbf{I})^{-1}\mathbf{A}^T\mathbf{x} s^=(ATA+λI)−1ATx
其中 λ > 0 \lambda > 0 λ>0 是正则化参数。
非线性混合的复杂性
现实中的混合过程通常是非线性的,这给盲源分离带来了巨大挑战。处理非线性混合的方法主要有:
1. 核方法(Kernel Methods)
核方法通过非线性映射 Φ : R m → F \Phi: \mathbb{R}^m \rightarrow \mathcal{F} Φ:Rm→F 将数据映射到高维特征空间 F \mathcal{F} F,然后在特征空间中寻求线性分离。核函数 k ( x , y ) = ⟨ Φ ( x ) , Φ ( y ) ⟩ k(\mathbf{x}, \mathbf{y}) = \langle \Phi(\mathbf{x}), \Phi(\mathbf{y}) \rangle k(x,y)=⟨Φ(x),Φ(y)⟩ 允许我们在不显式计算映射的情况下进行操作。
核PCA的目标函数为:
max ∥ α ∥ = 1 α T K α \max_{\|\mathbf{\alpha}\|=1} \mathbf{\alpha}^T\mathbf{K}\mathbf{\alpha} ∥α∥=1maxαTKα
其中 K \mathbf{K} K 是核矩阵, K i j = k ( x i , x j ) K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j) Kij=k(xi,xj)。
核ICA可以通过先进行核PCA降维,然后在新空间中应用线性ICA来实现。
2. 局部线性近似(Local Linear Approximation)
局部线性近似方法假设非线性映射在局部区域可以用线性变换近似。这类方法首先对数据进行聚类,然后在每个簇内应用线性ICA。
数学上,可以表示为在每个区域 R i \mathcal{R}_i Ri 中:
x ≈ A i s + b i , ∀ x ∈ R i \mathbf{x} \approx \mathbf{A}_i\mathbf{s} + \mathbf{b}_i, \quad \forall \mathbf{x} \in \mathcal{R}_i x≈Ais+bi,∀x∈Ri
其中 A i \mathbf{A}_i Ai 是局部线性变换矩阵, b i \mathbf{b}_i bi 是偏置向量。
3. 流形学习(Manifold Learning)
流形学习方法假设数据分布在低维流形上,通过保持局部几何结构来进行非线性降维。
常用的算法有等距映射(Isomap)、局部线性嵌入(LLE)、拉普拉斯特征映射(Laplacian Eigenmaps)等。
例如,LLE的优化问题为:
min Y ∑ i = 1 N ∥ y i − ∑ j ∈ N i w i j y j ∥ 2 2 s.t. Y Y T = I \min_{\mathbf{Y}} \sum_{i=1}^{N} \|\mathbf{y}_i - \sum_{j \in \mathcal{N}_i} w_{ij}\mathbf{y}_j\|_2^2 \quad \text{s.t.} \quad \mathbf{Y}\mathbf{Y}^T = \mathbf{I} Ymini=1∑N∥yi−j∈Ni∑wijyj∥22s.t.YYT=I
其中 Y = [ y 1 , y 2 , … , y N ] ∈ R d × N \mathbf{Y} = [\mathbf{y}_1, \mathbf{y}_2, \ldots, \mathbf{y}_N] \in \mathbb{R}^{d \times N} Y=[y1,y2,…,yN]∈Rd×N 是低维表示, N i \mathcal{N}_i Ni 是点 i i i 的近邻集合, w i j w_{ij} wij 是重构权重。