基于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} xnj=Ajynj+μj+ϵnjynjN(0,I),ϵnjN(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 πjs。不想传统的因子分析,假定 ϵ 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=1Nlog[j=1mπ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Σj1/2exp{21(xnμj)TΣj1(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=1Nj=1mznjlog(π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=1mIIj(θ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))=n1NE(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=1NRnj(θ(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)=N1n=1NRnj(θ(t))μj(t+1)=Nπj(t+1)1n=1NRnj(θ(t))xn.σj2(t+1)=(dkj)1=kj+1dλjAj(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=1m(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(kj1)/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=1m2Djlog(12Nπj)+2mlog12N+j=1m2Dj+1
这个标准的显著特征是:之前的AIC和BIC准则和参数都没有关系,MML准则却不同。1) 惩罚项包含模型参数 π j ′ s \pi_j's πjs和似然函数一起被联合优化;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 πjs的更新方程改变为:

π 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=1m2Djlog(Nπ^j)+2m1logN
其中 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 πjs: 因为约束条件 ∑ j = 1 m π j = 1 \sum_{j=1}^{m} \pi_{j}=1 j=1mπj=1和N个数据点都致力于估计 π j ′ s \pi_j's πjs,对应的 π j ′ s \pi_j's πjs的BIC惩罚项为 m − 1 2 log ⁡ N \frac{m-1}{2} \log N 2m1logN. 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=1m2Djlog(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 kjs看做参数,并引入一个额外的M-step关于 k j ′ s k_j's kjs最大化HBIC。通过这样,参数估计和关于 k j ′ s k_j's kjs的模型选择可以被同时解决。

当给定分量数 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=1m2Djlog(Nπj)2m1logN
目标是估计 ( 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=1mIj(θθ(t),π(t+1))2m1logN
其中,
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+(dk)logdkl=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(Σj1Sj(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(π)jp(μj)p(Σj)jΣj2gexp{21tr(Σj1B)}

其中B是正定阵,上述先意味着: Σ j \Sigma_j Σj使用带有 g g g个自由度的逆wishart先验, π j \pi_j πj μ j \mu_j μj使用无信息先验。

对于固定的 k j ′ s k_j's kjs,我们只需要替代 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πjn0时消除分量, n 0 n_0 n0是人为设定的。

5.4 两种算法

5.4.1 one-stage precedure

在这里插入图片描述

5.4.2 two-stage precedures

在这里插入图片描述

6 实验部分(详见相关论文)

  • 0
    点赞
  • 2
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

大浪淘沙_scc

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值