基于HBIC准则的混合PPCA的有效模型选择
注意:本博客来源于Efficient Model Selection for Mixtures of Probabilistic PCA via Hierarchical BIC,仅用于学习,转载请注明出处!
1 引言
混合MPCA模型提供了一种结合局部PCAs的主要的方法。在MPCA中,聚类和降维是同时实现的,并且参数估计可以很容易通过流行的EM算法被极大似然估计和最大后验估计方法得到。此外,在高维数据的密度建模中,MPCA提供了比高斯混合模型更显著的优势,因为它能够在过拟合的full GMM和欠拟合的diagonal/spheral GMM中提供一种合理的权衡。
尽管MPCA有大的吸引力,但它仍然有一个挑战性的问题,那就是混合成分的数量 m m m的确定和每个成分的子空间维度 k j k_j kj的确定。 为此,通常采用两阶段的方法来确定:stage1主要对一系列候选模型进行参数估计,stage2基于一种模型选择标准来选择最佳的模型。 在此,通常都假定共同- k k k-MPCA,也就是说假定每一个成分的子空间维度 k j = k k_j=k kj=k,然后在所有可能的集合 { m , k } \{m,k\} {m,k}所对应的一系列模型中进行穷举,然后基于模型选择标准选择最佳模型。如果没有共同- k k k的假定,在一个更大的集合 { m , k j } \{m,k_j\} {m,kj}中搜索会更加耗时。
直观上来看,不同的成分有不同的局部子空间维度 k j k_j kj是非常自然的,事实上,在图像处理中,通过仔细分配 k j k_j kj,图像压缩的质量将会显著提高。在手写字符识别中,带有不同的 k j k_j kj的成分的MPCA与共同- k k k模型相比将会获得更高的检验似然和更低的分类错误率。也有很多尝试来解决不同的局部子空间维度 k j k_j kj的问题。例如:“Resolution-based complexity control for gaussian mixture”一文中,限制所有分量的噪声方差 σ j 2 = σ 2 \sigma_j^2=\sigma^2 σj2=σ2,而每个分量的子空间维度 k j k_j kj是由大于 σ 2 \sigma^2 σ2的协方差矩阵的特征值的个数所决定的。除了这个限制,最优的 σ 2 \sigma^2 σ2是通过验证数据集来决定的。在PCA中的经典的方法是保证方差比大于一个阈值,这也在MPCA中被采用。缺点是 k j k_j kj的决定和参数估计是通过不同的损失函数来实现的。
为了有效地学习GMM,Figueiredo and Jain 提出了一种one-stage的方法,成为FJ算法。不像两阶段方法,此算法把把参数估计和模型选择集合到一个算法中。此外,在初始化上,此算法比标准的EM算法有更少的敏感性。可是,对于MPCA的FJ算法的应用有两个局限:其一,只有在局部子空间维度 k j k_j kj给定时才能应用到MPCA模型上;其二,在我们的实验中,分量消除的步骤使用的约束对一些基础的真实数据集来说都太强了。
变分贝叶斯方法被建议用来拟合MPCA模型。不像ML/MAP方法,VB保持模型参数的分布并且它的下界自然地作为模型选择的标准,在MPCA的VB处理中,不同的子空间的 k j k_j kj可以通过自动相关确定(automatic relemance determination ,ARD)自动的确定。为了进一步自动地确定混合成分的数量 m m m,一种birth/death的操作被融入VBMPCA方法中或者是采用一种非参数的处理方法。这样的方法导致参数估计和模型选择是在同一个算法中,可是,局部子空间 k j k_j kj的决定在这些基于VB的方法中都纯粹依赖于ARD方法,但它的结果对于预先设定的最大子空间维度的值是敏感的。
为了有效地学习MPCA,我们采用一种有效的模型选择标准称为HBIC,理论上,HBIC是VB下界的大样本极限,并且BIC是HBIC的进一步近似,为了有效地基于HBIC来学习MPCA,提出了两种算法:一种两阶段算法的变种和一种一阶段算法的变种。这两种算法都消除了common- k k k的限制,并且同时进行参数估计和模型选择,此外,one-stage算法也克服了FJ算法的两种局限。
2 Mixtures of probabilistic PCA(MPCA)
在MPCA模型的框架下,每一个的 d d d维的数据向量 x n \mathbf{x}_{n} xn都是i.i.d的样本 X = { x n } n = 1 N \mathbf{X}=\left\{\mathbf{x}_{n}\right\}_{n=1}^{N} X={xn}n=1N,它的生成需要两步。第一,基于在约束 ∑ j = 1 m π j = 1 \sum_{j=1}^{m} \pi_{j}=1 ∑j=1mπj=1下的分布 p ( j ) = π j , j = 1 , ⋯ , m p(j)=\pi_j,j = 1,\cdots,m p(j)=πj,j=1,⋯,m生成一个自然数 j j j。第二,给定自然数 j j j, x n \mathbf{x}_{n} xn由下文中限制性因子分析模型生成:
x n ∣ j = A j y n j + μ j + ϵ n j y n j ∼ N ( 0 , I ) , ϵ n j ∼ N ( 0 , σ j 2 I ) \begin{array}{l} \mathbf{x}_{n} \mid j=\mathbf{A}_{j} \mathbf{y}_{n j}+\boldsymbol{\mu}_{j}+\boldsymbol{\epsilon}_{n j} \\ \mathbf{y}_{n j} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}), \quad \boldsymbol{\epsilon}_{n j} \sim \mathcal{N}\left(\mathbf{0}, \sigma_{j}^{2} \mathbf{I}\right) \end{array} xn∣j=Ajynj+μj+ϵnjynj∼N(0,I),ϵnj∼N(0,σj2I)
I \mathbf{I} I表示单位阵。 μ j \boldsymbol{\mu}_{j} μj表示 d d d维的均值向量, A j \mathbf{A}_{j} Aj是一个 d × k j d \times k_j d×kj的因子载荷阵, y n j \mathbf{y}_{n j} ynj是一个独立的 k j k_j kj维的潜在因子向量, σ j 2 \sigma_{j}^{2} σj2是第 j j j个分量的噪声方差。很明显,这是 m m m个PPCA子模型的混合,其中混合比例为 π j ′ s \pi_j's πj′s。不想传统的因子分析,假定 ϵ n j \boldsymbol{\epsilon}_{n j} ϵnj是对角协方差,MPCA假定每一个分量有一个标量协方差。
定义: θ ≡ { π j , θ j ; j = 1 , … , m } , θ j = ( A j , μ j , σ j 2 ) \boldsymbol{\theta} \equiv\left\{\pi_{j}, \boldsymbol{\theta}_{j} ; j=1, \ldots, m\right\}, \boldsymbol{\theta}_{j}=\left(\mathbf{A}_{j}, \boldsymbol{\mu}_{j}, \sigma_{j}^{2}\right) θ≡{πj,θj;j=1,…,m},θj=(Aj,μj,σj2)
Σ j = A j A j ′ + σ j 2 I \Sigma_{j}=\mathbf{A}_{j} \mathbf{A}_{j}^{\prime}+\sigma_{j}^{2} \mathbf{I} Σj=AjAj′+σj2I
在MPCA模型下,观测数据的对数似然为:
L
(
X
∣
θ
)
=
∑
n
=
1
N
log
[
∑
j
=
1
m
π
j
p
(
x
n
∣
θ
j
)
]
\mathcal{L}(\mathbf{X} \mid \boldsymbol{\theta})=\sum_{n=1}^{N} \log \left[\sum_{j=1}^{m} \pi_{j} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}_{j}\right)\right]
L(X∣θ)=n=1∑Nlog[j=1∑mπjp(xn∣θj)]
其中,
p
(
x
n
∣
θ
j
)
=
(
2
π
)
−
d
/
2
∣
Σ
j
∣
−
1
/
2
⋅
exp
{
−
1
2
(
x
n
−
μ
j
)
T
Σ
j
−
1
(
x
n
−
μ
j
)
}
\begin{aligned} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}_{j}\right)=&(2 \pi)^{-d / 2}\left|\Sigma_{j}\right|^{-1 / 2} \\ & \cdot \exp \left\{-\frac{1}{2}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{j}\right)^{T} \Sigma_{j}^{-1}\left(\mathbf{x}_{n}-\mu_{j}\right)\right\} \end{aligned}
p(xn∣θj)=(2π)−d/2∣Σj∣−1/2⋅exp{−21(xn−μj)TΣj−1(xn−μj)}
ML估计
θ
^
\hat{\boldsymbol{\theta} }
θ^定义为:
θ
^
=
arg
max
θ
{
L
(
X
∣
θ
)
}
.
\hat{\boldsymbol{\theta}}=\underset{\boldsymbol{\theta}}{\arg \max }\{\mathcal{L}(\mathbf{X} \mid \boldsymbol{\theta})\} .
θ^=θargmax{L(X∣θ)}.
如果能获得参数
θ
\boldsymbol{\theta}
θ的先验分布
p
(
θ
)
p(\boldsymbol{\theta})
p(θ),然后MAP估计
θ
^
\hat{\boldsymbol{\theta} }
θ^定义为:
θ ^ = arg max θ { L ( X ∣ θ ) + log p ( θ ) } . \hat{\boldsymbol{\theta}}=\underset{\boldsymbol{\theta}}{\arg \max }\{\mathcal{L}(\mathbf{X} \mid \boldsymbol{\theta})+\log p(\boldsymbol{\theta})\} . θ^=θargmax{L(X∣θ)+logp(θ)}.
2.1 EM 算法
给定分量的数量 m m m和子空间维度 k = ( k 1 , k 2 , ⋯ , k m ) \mathbf{k}=\left(k_{1}, k_{2}, \cdots, k_{m}\right) k=(k1,k2,⋯,km),MPCA中的参数 θ \boldsymbol{\theta} θ可以通过EM算法求解。
考虑完全数据 ( X , Z ) = { x n , z n } n = 1 N (\mathbf{X}, \mathbf{Z})=\left\{\mathbf{x}_{n}, \mathbf{z}_{n}\right\}_{n=1}^{N} (X,Z)={xn,zn}n=1N,其中 z n = ( z n 1 , ⋯ , z n M ) \mathbf{z}_{n} = (z_{n1},\cdots,z_{nM}) zn=(zn1,⋯,znM), z n j z_{nj} znj是一个示性变量如果 x n \mathbf{x}_{n} xn来自分量 j j j,则 z n j z_{nj} znj为1,其他为0.
完全数据似然函数可表达为:
L c ( X , Z ∣ θ ) = ∑ n = 1 N ∑ j = 1 m z n j log ( π j p ( x n ∣ θ j ) ) \mathcal{L}_{c}(\mathbf{X}, \mathbf{Z} \mid \boldsymbol{\theta})=\sum_{n=1}^{N} \sum_{j=1}^{m} z_{n j} \log \left(\pi_{j} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}_{j}\right)\right) Lc(X,Z∣θ)=n=1∑Nj=1∑mznjlog(πjp(xn∣θj))
E-step:在给定
X
\mathbf{X}
X 和
θ
(
t
)
\boldsymbol{\theta}^{(t)}
θ(t)时计算
L
c
\mathcal{L}_{c}
Lc的期望:
Q
(
θ
∣
θ
(
t
)
)
=
∑
j
=
1
m
I
I
j
(
θ
j
∣
θ
(
t
)
)
Q\left(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)}\right)=\sum_{j=1}^{m} \mathrm{II}_{j}\left(\boldsymbol{\theta}_{j} \mid \boldsymbol{\theta}^{(t)}\right)
Q(θ∣θ(t))=j=1∑mIIj(θj∣θ(t))
其中,
I
I
j
(
θ
j
∣
θ
(
t
)
)
=
∑
n
−
1
N
E
(
z
n
j
)
log
(
π
j
p
(
x
n
∣
θ
j
)
)
\mathrm{II}_{j}\left(\boldsymbol{\theta}_{j} \mid \boldsymbol{\theta}^{(t)}\right)=\sum_{n-1}^{N} \mathbb{E}\left(z_{n j}\right) \log \left(\pi_{j} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}_{j}\right)\right)
IIj(θj∣θ(t))=n−1∑NE(znj)log(πjp(xn∣θj))
I
I
j
\mathrm{II}_{j}
IIj仅仅依赖第
j
j
j个分量的
(
π
j
,
θ
j
)
\left(\pi_{j}, \boldsymbol{\theta}_{j}\right)
(πj,θj),定义期望
R
n
j
(
θ
(
t
)
)
≜
E
(
z
n
j
)
R_{n j}\left(\boldsymbol{\theta}^{(t)}\right) \triangleq \mathbb{E}\left(z_{n j}\right)
Rnj(θ(t))≜E(znj)是数据点
x
n
\mathbf{x}_{n}
xn属于第
j
j
j个分量的后验概率:
R
n
j
(
θ
(
t
)
)
=
P
(
z
n
j
=
1
∣
θ
(
t
)
)
=
π
j
(
t
)
p
(
x
n
∣
θ
j
(
t
)
)
∑
k
=
1
m
π
k
(
t
)
p
(
x
n
∣
θ
k
(
t
)
)
R_{n j}\left(\boldsymbol{\theta}^{(t)}\right)=P\left(z_{n j}=1 \mid \boldsymbol{\theta}^{(t)}\right)=\frac{\pi_{j}^{(t)} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}_{j}^{(t)}\right)}{\sum_{k=1}^{m} \pi_{k}^{(t)} p\left(\mathbf{x}_{n} \mid \boldsymbol{\theta}_{k}^{(t)}\right)}
Rnj(θ(t))=P(znj=1∣θ(t))=∑k=1mπk(t)p(xn∣θk(t))πj(t)p(xn∣θj(t))
此外,我们定义一个局部样本协方差矩阵:
S
j
(
t
+
1
)
=
1
N
π
j
(
t
+
1
)
∑
n
=
1
N
R
n
j
(
θ
(
t
)
)
(
x
n
−
μ
j
(
t
+
1
)
)
(
x
n
−
μ
j
(
t
+
1
)
)
′
\mathbf{S}_{j}^{(t+1)}=\frac{1}{N \pi_{j}^{(t+1)}} \sum_{n=1}^{N} R_{n j}\left(\boldsymbol{\theta}^{(t)}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{j}^{(t+1)}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{j}^{(t+1)}\right)^{\prime}
Sj(t+1)=Nπj(t+1)1n=1∑NRnj(θ(t))(xn−μj(t+1))(xn−μj(t+1))′
M-step: 对于ML估计,在约束 ∑ j = 1 m π j = 1 \sum_{j=1}^{m} \pi_{j}=1 ∑j=1mπj=1下关于 θ \boldsymbol{\theta} θ最大化Q函数产生如下更新方程:
π j ( t + 1 ) = 1 N ∑ n = 1 N R n j ( θ ( t ) ) μ j ( t + 1 ) = 1 N π j ( t + 1 ) ∑ n = 1 N R n j ( θ ( t ) ) x n . σ j 2 ( t + 1 ) = ( d − k j ) − 1 ∑ ℓ = k j + 1 d λ j ℓ A j ( t + 1 ) = U j ( Λ j − σ j 2 ( t + 1 ) I ) 1 / 2 \begin{array}{c} \pi_{j}^{(t+1)}=\frac{1}{N} \sum_{n=1}^{N} R_{n j}\left(\theta^{(t)}\right) \\ \mu_{j}^{(t+1)}=\frac{1}{N \pi_{j}^{(t+1)}} \sum_{n=1}^{N} R_{n j}\left(\boldsymbol{\theta}^{(t)}\right) \mathbf{x}_{n} . \\ \sigma_{j}^{2(t+1)}=\left(d-k_{j}\right)^{-1} \sum_{\ell=k_{j}+1}^{d} \lambda_{j \ell} \\ \mathbf{A}_{j}^{(t+1)}=\mathbf{U}_{j}\left(\Lambda_{j}-\sigma_{j}^{2(t+1)} \mathbf{I}\right)^{1 / 2} \end{array} πj(t+1)=N1∑n=1NRnj(θ(t))μj(t+1)=Nπj(t+1)1∑n=1NRnj(θ(t))xn.σj2(t+1)=(d−kj)−1∑ℓ=kj+1dλjℓAj(t+1)=Uj(Λj−σj2(t+1)I)1/2
其中, U j = ( u j 1 , … , u j k j ) , Λ j = diag ( λ j 1 , … , λ j k j ) \mathbf{U}_{j}=\left(\mathbf{u}_{j 1}, \ldots, \mathbf{u}_{j k_{j}}\right), \quad \Lambda_{j}=\operatorname{diag}\left(\lambda_{j 1}, \ldots, \lambda_{j k_{j}}\right) Uj=(uj1,…,ujkj),Λj=diag(λj1,…,λjkj), { λ j ℓ } ℓ = 1 d , { u j ℓ } ℓ = 1 d \left\{\lambda_{j \ell}\right\}_{\ell=1}^{d},\left\{\mathbf{u}_{j \ell}\right\}_{\ell=1}^{d} {λjℓ}ℓ=1d,{ujℓ}ℓ=1d是 S j ( t + 1 ) \mathbf{S}_{j}^{(t+1)} Sj(t+1)降序特征值和特征向量。
对于MAP估计, θ ( t + 1 ) \boldsymbol{\theta}^{(t+1)} θ(t+1)可以通过如下函数获得:
θ
(
t
+
1
)
=
arg
max
θ
{
Q
(
θ
∣
θ
(
t
)
)
+
log
p
(
θ
)
}
.
\boldsymbol{\theta}^{(t+1)}=\underset{\boldsymbol{\theta}}{\arg \max }\left\{Q\left(\boldsymbol{\theta} \mid \boldsymbol{\theta}^{(t)}\right)+\log p(\boldsymbol{\theta})\right\} .
θ(t+1)=θargmax{Q(θ∣θ(t))+logp(θ)}.
其中,
θ
(
t
+
1
)
\boldsymbol{\theta}^{(t+1)}
θ(t+1)的更新公式依赖于先验分布
p
(
θ
)
p(\boldsymbol{\theta})
p(θ)的选择。
3 模型选择方法
模型选择方法从一个计算的角度,可以划分为两类:two-stage和one-stage方法。
3.1 two-stage methods
two-stage 方法在一个集合 ( m , k ) (m, \mathbf{k}) (m,k),其中 m m m的范围 [ m m i n , m m a x ] [m_{min},m_{max}] [mmin,mmax],每一个 k j k_j kj的范围是 [ k m i n , k m a x ] [k_{min},k_{max}] [kmin,kmax],假定最优的模型被包含在里面,two-stage方法首先获得每一个模型的ML/MAP估计 θ ^ ( m , k ) \hat{\boldsymbol{\theta}}(m, \mathbf{k}) θ^(m,k),然后利用模型选择标准 L ∗ ( m , k , θ ^ ( m , k ) ) \mathcal{L}^{*}(m, \mathbf{k}, \hat{\boldsymbol{\theta}}(m, \mathbf{k})) L∗(m,k,θ^(m,k))来选择最佳的模型:
(
m
^
,
k
^
)
=
arg
max
(
m
,
k
)
{
L
∗
(
m
,
k
,
θ
^
(
m
,
k
)
)
}
.
(\hat{m}, \hat{\mathbf{k}})=\underset{(m, \mathbf{k})}{\arg \max }\left\{\mathcal{L}^{*}(m, \mathbf{k}, \hat{\boldsymbol{\theta}}(m, \mathbf{k}))\right\} .
(m^,k^)=(m,k)argmax{L∗(m,k,θ^(m,k))}.
已经提出了很多模型选择标准,例如AIC,BIC,LEC,ICL,MDL等,但是所有的标准都可以被写成如下形式:
L ∗ ( m , k , θ ^ ( m , k ) ) = L ( X ∣ θ ^ ( m , k ) ) − P ( θ ^ ( m , k ) ) \mathcal{L}^{*}(m, \mathbf{k}, \hat{\boldsymbol{\theta}}(m, \mathbf{k}))=\mathcal{L}(\mathbf{X} \mid \hat{\boldsymbol{\theta}}(m, \mathbf{k}))-\mathcal{P}(\hat{\boldsymbol{\theta}}(m, \mathbf{k})) L∗(m,k,θ^(m,k))=L(X∣θ^(m,k))−P(θ^(m,k))
其中, P ( θ ^ ( m , k ) ) \mathcal{P}(\hat{\boldsymbol{\theta}}(m, \mathbf{k})) P(θ^(m,k))是惩罚项,用于惩罚更高的 ( m , k ) (m, \mathbf{k}) (m,k)。
由于BIC理论上的相合性和满意的实际应用性能,BIC对于混合模型来说更流行。MDL和BIC有相同的形式,但是它来自于信息论或者编码论。
对于 m m m个分量的MPCA,BIC的惩罚项如下:
P
b
i
c
(
θ
^
)
=
1
2
[
∑
j
=
1
m
(
D
j
+
1
)
−
1
]
log
N
\mathcal{P}_{b i c}(\hat{\boldsymbol{\theta}})=\frac{1}{2}\left[\sum_{j=1}^{m}\left(\mathcal{D}_{j}+1\right)-1\right] \log N
Pbic(θ^)=21[j=1∑m(Dj+1)−1]logN
其中,
D
j
=
d
(
k
j
+
1
)
−
k
j
(
k
j
−
1
)
/
2
+
1
\mathcal{D}_{j}=d\left(k_{j}+1\right)-k_{j}\left(k_{j}-1\right) / 2+1
Dj=d(kj+1)−kj(kj−1)/2+1是
θ
j
\boldsymbol{\theta}_{j}
θj的自有参数的个数。可是BIC使用整个样本量
N
N
N来决定第
j
j
j个分析器
D
j
\mathcal{D}_{j}
Dj。
为了使用two-stage方法,通常假定common- k k k模型,否则是太耗时了。
3.2 one-stage 方法
不像两阶段方法,one-stage方法将参数估计和混合模型的分量数 m m m的确定集成到一个单独的算法中,这在计算上更有效。一个有代表性的例子就是FJ算法。作者提出了一种MML标准:
θ ^ = arg max θ { L ( X ∣ θ ) − P m m l ( θ ) } \hat{\boldsymbol{\theta}}=\underset{\theta}{\arg \max }\left\{\mathcal{L}(\mathbf{X} \mid \boldsymbol{\theta})-\mathcal{P}_{\boldsymbol{m m} l}(\boldsymbol{\theta})\right\} θ^=θargmax{L(X∣θ)−Pmml(θ)}
其中,
P
m
m
l
(
θ
)
=
∑
j
=
1
m
D
j
2
log
(
N
π
j
12
)
+
m
2
log
N
12
+
∑
j
=
1
m
D
j
+
1
2
\mathcal{P}_{m m l}(\boldsymbol{\theta})=\sum_{j=1}^{m} \frac{\mathcal{D}_{j}}{2} \log \left(\frac{N \pi_{j}}{12}\right)+\frac{m}{2} \log \frac{N}{12}+\sum_{j=1}^{m} \frac{\mathcal{D}_{j}+1}{2}
Pmml(θ)=j=1∑m2Djlog(12Nπj)+2mlog12N+j=1∑m2Dj+1
这个标准的显著特征是:之前的AIC和BIC准则和参数都没有关系,MML准则却不同。1) 惩罚项包含模型参数
π
j
′
s
\pi_j's
πj′s和似然函数一起被联合优化;2)新的标准关于
θ
\boldsymbol{\theta}
θ和
D
j
\mathcal{D}_{j}
Dj同时最大化。为了做到这些,作者提出了一种one-stage算法(FJ算法),实现只需要在原来的EM算法中做简单的修改就可以了。
FJ算法有两个局限:第一,只有在局部子空间维度 k j k_j kj给定时才能应用到MPCA模型上;第二,在FJ算法中,关于参数 θ \boldsymbol{\theta} θ最大化目标函数的更新方程和上述EM算法的更新方程是一样的,除了混合比例 π j ′ s \pi_j's πj′s的更新方程改变为:
π j ( t + 1 ) = max { 0 , ∑ n = 1 N R n j ( θ ( t ) ) − D j 2 } ∑ j = 1 m max { 0 , ∑ n = 1 N R n j ( θ ( t ) ) − D j 2 } \pi_{j}^{(t+1)}=\frac{\max \left\{0, \sum_{n=1}^{N} R_{n j}\left(\boldsymbol{\theta}^{(t)}\right)-\frac{\mathcal{D}_{j}}{2}\right\}}{\sum_{j=1}^{m} \max \left\{0, \sum_{n=1}^{N} R_{n j}\left(\boldsymbol{\theta}^{(t)}\right)-\frac{\mathcal{D}_{j}}{2}\right\}} πj(t+1)=∑j=1mmax{0,∑n=1NRnj(θ(t))−2Dj}max{0,∑n=1NRnj(θ(t))−2Dj}
上述这个公式自动地消除了数据支持不是很好的第 j j j个分量,即 ∑ n = 1 N R n j < = D j / 2 \sum_{n=1}^{N} R_{n j}<=\mathcal{D}_{j}/2 ∑n=1NRnj<=Dj/2。尽管上述公式具有吸引人的特征,但我们也发现,这个约束太强了而不能在用在很多真实的数据集。
4 Hierarchical BIC(H-BIC)
H-BIC准则如下:
P
h
b
i
c
(
θ
^
)
=
∑
j
=
1
m
D
j
2
log
(
N
π
^
j
)
+
m
−
1
2
log
N
\mathcal{P}_{h b i c}(\hat{\boldsymbol{\theta}})=\sum_{j=1}^{m} \frac{\mathcal{D}_{j}}{2} \log \left(N \hat{\pi}_{j}\right)+\frac{m-1}{2} \log N
Phbic(θ^)=j=1∑m2Djlog(Nπ^j)+2m−1logN
其中
D
j
\mathcal{D}_{j}
Dj是
θ
j
\boldsymbol{\theta}_{j}
θj自由参数的个数,
π
^
j
\hat{\pi}_j
π^j是ML/MAP估计。
HBIC是VB下界的大样本极限,证明请详细看相关论文。
4.1 BIC和MML的区别和联系
- Hbic和bic的联系:其一,对于单分量混合模型 m = 1 m=1 m=1时,HBIC退化成BIC;其二,当 m > 1 m>1 m>1时,BIC比HBIC惩罚的更重。其三,在大样本情形下,BIC是HBIC的进一步近似;其四,对于HBIC,是通过分层形式的BIC来实现的,因为1)****BIC for π j ′ s \pi_j's πj′s: 因为约束条件 ∑ j = 1 m π j = 1 \sum_{j=1}^{m} \pi_{j}=1 ∑j=1mπj=1和N个数据点都致力于估计 π j ′ s \pi_j's πj′s,对应的 π j ′ s \pi_j's πj′s的BIC惩罚项为 m − 1 2 log N \frac{m-1}{2} \log N 2m−1logN. 2) 每一个分量的局部BIC: N π ^ j N \hat{\pi}_{j} Nπ^j可以看做有效的样本量。
回顾在MPCA模型下的数据生成的两步过程,首先基于 π j \pi_j πj生成分量标签 j j j,然后基于 θ j \boldsymbol{\theta}_{j} θj生成数据 x n \mathbf{x}_n xn,因此,上述标准可以看做是分层BIC准则,每一层对应数据点的生成步骤。
- HBIC和MML的区别和联系
1)可以发现,hbic中包含的是 π ^ j \hat{\pi}_j π^j,而MML包含的是需要和似然函数一起优化的 π j \pi_j πj,因此MML的解并不是ML估计和MAP估计。此外, 12 / e = 4.4146 12/e=4.4146 12/e=4.4146,因此mml的惩罚项可重写为:
P m m l ( θ ) ≈ ∑ j = 1 m D j 2 log ( N π j 4.4146 ) + m 2 log N 4.4146 \mathcal{P}_{m m l}(\boldsymbol{\theta}) \approx \sum_{j=1}^{m} \frac{\mathcal{D}_{j}}{2} \log \left(\frac{N \pi_{j}}{4.4146}\right)+\frac{m}{2} \log \frac{N}{4.4146} Pmml(θ)≈j=1∑m2Djlog(4.4146Nπj)+2mlog4.4146N
可以看到:mml的惩罚比HBIC的更轻。
由于FJ算法遭遇了很强的约束条件,下文中基于HBIC提出相关的算法。其次,mml选择的模型通常带有太多的谬误的分量(仅仅包含很少的数据点),主要是因为mml的惩罚项对太小的 N π j N \pi_j Nπj不敏感。特别地,当 N π j < 4.4 N \pi_j<4.4 Nπj<4.4,事实上,该分量的惩罚是负的,因此可能会遇到不显著的分量。
5 两个算法
显然,在MPCA中,经典的使用BIC的two-stage方法也可以用来实现HBIC。可是,即使限制common- k k k,这仍然是耗时的。本节主要提出两个有效的算法:两阶段变种和一阶段变种算法,且这两种算法都是基于修改的EM算法。
5.1 修改的EM算法
在这里,我们将维度 k j ′ s k_j's kj′s看做参数,并引入一个额外的M-step关于 k j ′ s k_j's kj′s最大化HBIC。通过这样,参数估计和关于 k j ′ s k_j's kj′s的模型选择可以被同时解决。
当给定分量数
m
m
m时,目标函数为:
L
∗
(
θ
,
k
)
=
L
(
θ
)
−
∑
j
=
1
m
D
j
2
log
(
N
π
j
)
−
m
−
1
2
log
N
\mathcal{L}^{*}(\boldsymbol{\theta}, \mathbf{k})=\mathcal{L}(\boldsymbol{\theta})-\sum_{j=1}^{m} \frac{\mathcal{D}_{j}}{2} \log \left(N \pi_{j}\right)-\frac{m-1}{2} \log N
L∗(θ,k)=L(θ)−j=1∑m2Djlog(Nπj)−2m−1logN
目标是估计
(
k
^
,
θ
^
(
k
^
)
)
(\hat{\mathbf{k}}, \hat{\boldsymbol{\theta}}(\hat{\mathbf{k}}))
(k^,θ^(k^))。
k
^
\hat{\mathbf{k}}
k^是具有最大HBIC的模型。但是,注意到,最大化上述的函数类似于最大化mml目标函数,这并不能ML/MAP估计
θ
\boldsymbol{\theta}
θ。因此,我们需要使用包含
θ
^
\hat{\boldsymbol{\theta}}
θ^的EM 算法。
当给定 ( π j ( t + 1 ) , μ j ( t + 1 ) ) ′ s \left(\pi_{j}^{(t+1)}, \mu_{j}^{(t+1)}\right)^{\prime} \mathrm{s} (πj(t+1),μj(t+1))′s,新的对应的惩罚函数为:
Q
∗
(
θ
∣
θ
(
t
)
,
π
(
t
+
1
)
)
=
∑
j
=
1
m
I
j
∗
(
θ
∣
θ
(
t
)
,
π
(
t
+
1
)
)
−
m
−
1
2
log
N
Q^{*}\left(\theta \mid \theta^{(t)}, \pi^{(t+1)}\right)=\sum_{j=1}^{m} \mathrm{I}_{j}^{*}\left(\theta \mid \theta^{(t)}, \pi^{(t+1)}\right)-\frac{m-1}{2} \log N
Q∗(θ∣θ(t),π(t+1))=j=1∑mIj∗(θ∣θ(t),π(t+1))−2m−1logN
其中,
I
I
j
∗
(
θ
∣
θ
(
t
)
,
π
(
t
+
1
)
)
=
I
I
j
(
θ
∣
θ
(
t
)
)
−
D
j
2
log
(
N
π
j
(
t
+
1
)
)
\mathrm{II}_{j}^{*}\left(\theta \mid \theta^{(t)}, \pi^{(t+1)}\right)=\mathrm{II}_{j}\left(\theta \mid \theta^{(t)}\right)-\frac{\mathcal{D}_{j}}{2} \log \left(N \pi_{j}^{(t+1)}\right)
IIj∗(θ∣θ(t),π(t+1))=IIj(θ∣θ(t))−2Djlog(Nπj(t+1))
修改的EM算法按如下步骤接待:
step1: 最大化Q关于
θ
\theta
θ获得
θ
(
t
+
1
)
\theta^{(t+1)}
θ(t+1)
step2:最大化
Q
∗
Q^*
Q∗关于
k
\mathbf{k}
k获得
k
(
t
+
1
)
\mathbf{k}^{(t+1)}
k(t+1)
k
j
k_j
kj的闭式表达按如下获得:
k
j
(
t
+
1
)
=
arg
min
k
{
N
π
j
(
t
+
1
)
(
∑
l
=
1
k
log
λ
j
l
+
(
d
−
k
)
⋅
log
∑
l
=
k
+
1
d
λ
j
l
d
−
k
)
+
D
j
log
(
N
π
j
(
t
+
1
)
)
}
\begin{array}{l} k_{j}^{(t+1)}=\underset{k}{\arg \min }\left\{N \pi_{j}^{(t+1)}\left(\sum_{l=1}^{k} \log \lambda_{j l}+(d-k) \cdot \log \frac{\sum_{l=k+1}^{d} \lambda_{j l}}{d-k}\right)\right. \\ \left.\quad+\mathcal{D}_{j} \log \left(N \pi_{j}^{(t+1)}\right)\right\} \end{array}
kj(t+1)=kargmin{Nπj(t+1)(∑l=1klogλjl+(d−k)⋅logd−k∑l=k+1dλjl)+Djlog(Nπj(t+1))}
其中,
tr
(
Σ
j
−
1
S
j
(
t
+
1
)
)
=
d
\operatorname{tr}\left(\Sigma_{j}^{-1} \mathbf{S}_{j}^{(t+1)}\right)=d
tr(Σj−1Sj(t+1))=d,丢掉与
k
j
k_j
kj不相关的项就得到上式。
5.2 最大后验估计
ML估计可能会接触到参数空间的边界,因此会遭受奇异值问题。在这种情况下,HBIC的值不能被计算。当分量数
m
m
mh和子空间维度
k
j
k_j
kj比最优的大的时候,这种情况将会频繁发生。为了处理这种问题,使用MAP估计,因为他比ML估计更稳定,HBIC总是可以计算。
使用如下先验:
p ( θ ) = p ( π ) ∏ j p ( μ j ) p ( Σ j ) ∝ ∏ j ∣ Σ j ∣ − g 2 exp { − 1 2 tr ( Σ j − 1 B ) } \begin{aligned} p(\theta)=& p(\pi) \prod_{j} p\left(\mu_{j}\right) p\left(\Sigma_{j}\right) \\ & \propto \prod_{j}\left|\Sigma_{j}\right|^{-\frac{g}{2}} \exp \left\{-\frac{1}{2} \operatorname{tr}\left(\Sigma_{j}^{-1} \mathbf{B}\right)\right\} \end{aligned} p(θ)=p(π)j∏p(μj)p(Σj)∝j∏∣Σj∣−2gexp{−21tr(Σj−1B)}
其中B是正定阵,上述先意味着: Σ j \Sigma_j Σj使用带有 g g g个自由度的逆wishart先验, π j \pi_j πj和 μ j \mu_j μj使用无信息先验。
对于固定的 k j ′ s k_j's kj′s,我们只需要替代 S j ( t + 1 ) \mathbf{S}_{j}^{(t+1)} Sj(t+1)为:
S
~
j
(
t
+
1
)
=
1
N
π
j
(
t
+
1
)
+
g
(
N
π
j
(
t
+
1
)
S
j
(
t
+
1
)
+
B
)
=
(
1
−
γ
)
S
j
(
t
+
1
)
+
γ
B
g
\begin{aligned} \tilde{\mathbf{S}}_{j}^{(t+1)} &=\frac{1}{N \pi_{j}^{(t+1)}+g}\left(N \pi_{j}^{(t+1)} \mathbf{S}_{j}^{(t+1)}+\mathbf{B}\right) \\ &=(1-\gamma) \mathbf{S}_{j}^{(t+1)}+\gamma \frac{\mathbf{B}}{g} \end{aligned}
S~j(t+1)=Nπj(t+1)+g1(Nπj(t+1)Sj(t+1)+B)=(1−γ)Sj(t+1)+γgB
其中,
γ
=
g
/
(
N
π
j
(
t
+
1
)
+
g
)
.
\gamma=g /\left(N \pi_{j}^{(t+1)}+g\right) .
γ=g/(Nπj(t+1)+g).
如果我们将B设为标量协方差矩阵,则 S ~ j ( t + 1 ) \tilde{\mathbf{S}}_{j}^{(t+1)} S~j(t+1)实现一个与RDA相似的正则化
5.3 消除分量
当 N π j ≤ n 0 N \pi_{j} \leq n_{0} Nπj≤n0时消除分量, n 0 n_0 n0是人为设定的。