贝叶斯机器学习:最大熵及高斯分布

高斯分布[1],也被称为正态分布,广泛应用于连续型随机变量分布的模型中。对于一元变量x的情形。高斯分布可以写成下列的形式:

N(xμ,σ2)=1(2πσ2)12exp{12σ2(xμ)2}

其中μ是均值,σ2是方差。对于D维向量x,多元高斯分布的形式为:

(1)N(xμ,Σ)=1(2π)D21(detΣ)12exp{12(xμ)TΣ1(xμ)}

其中,μ是一个D维均值向量,Σ是一个D×D协方差矩阵(covariance matrix)detΣΣ的行列式。

1 高斯分布的物理意义

1.1 做为最大熵分布的高斯分布

高斯分布可以从多个不同的角度来理解。例如,对于一个一元实值向量,使得熵取得最大值的是高斯分布。这个性质对于多元高斯分布也成立。

熵(entropy) 的概念最早起源于物理学,是在热力学平衡的背景中介绍的。后来,熵称为描述统计力学中的无序程度的度量。在统计力学中,玻尔兹曼熵(Boltzmann's entropy)[2][3]的定义为:

SklnΩ

这里Ω为系统宏观态的重数,k为玻尔兹曼常量。

关于统计力学的术语我们用一个例子[4]来简要做一下介绍。考虑将3个不同的球放进3个不同的箱子中,全部的33=27种可能放法如下图所示:

这27种不同结果中的每一种都称为微观态(microstate)。通常,在统计力学中,为了知道系统的微观态,我们必须清楚每个粒子的状态,在我们这个例子中是每个球划分到箱子中去的状态。如果更一般地指定状态——例如仅仅说各个箱子里有多少个球,我们称它为宏观态(macrostate)。当然,如果知道一个系统的微观态(比如{ab   c  }),那么我们肯定也能知道它的宏观态(比如箱子中的球个数分别为2、1、0)。但反过来却不行:知道箱子中的球个数分别为2、1、0并没有告诉我们每个球的状态,因为有3个微观态都对应于这个宏观态。对应于给定宏观态的微观态数量称为该宏观态的重数(multiplicity),在这种情况下为3。
更一般地,考虑将N个球放进M个箱子里,各个箱子里有n1,n2,,nM个球时宏观态的重数为我们在博客《概率论沉思录:初等抽样论》中提到的多项式系数:

Ω(n1,n2,,nM)=N!iMni!

其中i=1Mni=N

对于将N个不同的球放进M个不同的箱子所形成的这样一个系统,各个箱子里有n1,n2,,nM个球时系统的熵为

S=klnN!iMni!=klnN!ki=1Mlnni!

现在我们在此基础上乘以一个缩放参数1N,忽略掉常数k,并考虑N1,根据斯特林近似(Stirlings approximation) 我们有lnN!NlnNN,于是

1NS1N(NlnNN)1Ni=1M[nilnnini]=lnNi=1MniNlnni=i=1M(niN)ln(niN)

斯特林近似为:

N!2πN(Ne)N

N1时,这个公式十分准确。该公式可以采用下列的直观方式进行理解。N!是从1NN个因子的乘积,一个十分粗略的估计是把每个因子都替换成N,这就是N!NN。这是一个过高的估计,因为几乎所有的因子都比N小,平均下来每个因子都大了大约e倍,也就是说

N!(Ne)N

这仍比N!差了一个大数字的因子,大约是2πN。但若N是一个大数字,那么N!会是一个非常大的数字,所以这个(大数字)因子就可以被省略了。如果我们只关心N!的对数,通常上式就已经足够准确:

lnN!NlnNN

取极限N,并保持比值niN固定,我们可以得到

limNi=1M(niN)ln(niN)=i=1Mpilnpi

这里pi=limN(niN)为给定系统宏观态时,一个球被分配到第i个箱子的概率。对于一个球,它被分配到的箱子可以表述成一个离散随机变量XX一共有M个状态),且p(X=xi)=pi (i=1,,M)。这样,随机变量X的熵就为

H[p]=i=1Mp(xi)lnp(xi)

如果分布p(xi)在几个值周围有尖锐的峰值,熵就会相对较低。如果分布p(xi)相对平衡地跨过许多值,那么熵就会相对较高。下图展示了两个概率分布在30个箱子上的直方图:

可以看到,熵值越大,H越宽(最大的熵值产生于均匀分布,此时的熵值为H=ln(130)3.40)。

由于0pi1,因此熵是非负的(由于当x0时,xlnx0,我们约定0ln0=0)。当pi=1且所有其他的pji=0时,熵取得最小值0。在ip(xi)=1(概率归一化)的约束下,使用拉格朗日乘数法可以找到熵的最大值。因此,我们要最大化

H~=i=1Mp(xi)lnp(xi)+λ(i=1Mp(xi)1)

d(H~)dp(xi)=0,解得p(xi)=eλ1,然后代回约束条件ip(xi),得到λ=1lnM,于是最优解为p(xi)=1M(此时所有的p(xi)都相等),此时熵取得最大值lnM

我们可以把熵的定义扩展到连续随机变量X的概率分布p(x),方法如下。首先把连续随机变量X的取值范围切分成宽度为Δ的小段区间。根据均值定理(mean value theorem)[5],对于每个这样的小段区间Δ,一定在区间中存在一个值xi使得

iΔ(i+1)Δp(x)dx=p(xi)Δ

现在我们可以这样量化连续随机变量X:只要X落在第i个区间中,我们就把X赋值为xi。因此观察到值xi的概率为p(xi)Δ。这样就变成了离散的分布,这种情况下熵的形式为:

HΔ=i=1(p(xi)Δ)ln(p(xi)Δ)=i=1p(xi)Δlnp(xi)lnΔ

推导时我们使用了ip(xi)Δ,这是因为ip(xi)Δ=iiΔ(i+1)Δp(x)dx=1。我们现在省略上式右侧的第二项lnΔ,然后考虑极限Δ0,此时有:

limΔ0{i=1p(xi)Δlnp(xi)}=p(x)lnp(x)dx

其中,右侧的量被称为微分熵(differential entropy)。我们看到,熵的离散形式与连续形式的差是lnΔ,这在极限Δ0的情形下发散。这反映出一个事实:具体化一个连续随机变量需要大量的比特位。和离散情形类似,我们将微分熵记为关于概率密度函数p(x)的函数:

H[p]=p(x)lnp(x)dx

在离散分布的情况下,我们看到最大熵对应于变量所有可能状态的均匀分布。现在让我们考虑连续随机变量的最大熵。为了让这个最大值有一个合理的定义,除了保留归一化的约束之外,还有必要限制p(x)的均值(一阶矩)和方差(二阶矩)。之所以需要限制其方差,是因为当方差增大时,熵也会无限制地增加,因此除非我们给定固定的方差σ2,否则寻找哪一个分布有最大熵这个问题是没有意义的。之所以需要限制其均值,是因为在不改变熵的条件下一个分布可以被随意地改变,因此我们需要再加一个均值为μ的约束以获得一个唯一的解[6]。综上,我们最大化微分熵的时候要施加下面三个约束:

p(x)dx=1xp(x)dx=μ(xμ)2p(x)dx=σ2

那么这个问题的拉格朗日泛函如下:

J(p)=H[p]+λ1(p(x)dx1)+λ2(xp(x)dxμ)+λ3((xμ)2p(x)dxσ2)

根据变分法,可以得到泛函导数为:

δJδp(x)=(lnp(x)+1)+λ1+λ2x+λ3(xμ)2

δJδp(x)=0(对x),解得p(x)=exp{1+λ1+λ2x+λ3(xμ)2}

函数f的函数被称为泛函(functional)J(f)。正如许多情况下对一个函数求关于以向量的各元素为变量的偏导数一样,我们可以使用泛函导数(functional derivative),即在任意特定的x值,对一个泛函J(f)求关于函数f(x)的导数,这也被称为变分导数(variational derivative)。泛函J的关于函数f在点x处的泛函导数被记作δJδf(x)

对于可微分函数f(x)以及带有连续导数的可微分函数g(y,x)(其中y=f(x)),有下列结论:

δδf(x)g(f(x),x)dx=yg(f(x),x)

为了使上述等式更加直观,我们可以把f(x)看做是一个有着无穷不可数多个元素的向量。在这里,这种关系式中描述的对于特定的x,关于f(x)的泛函导数可以类比为对于特定的下标i,关于向量θRn的第i个元素的偏导数:

θijg(θj,j)=θig(θi,i)

为了关于一个向量优化某个函数,我们可以求出这个函数关于这个向量的梯度,然后找到使这个梯度中每一个元素都为0的点。类似地,我们可以通过寻找一个函数使得泛函导数在每个点上都等于0,从而来优化一个泛函。

为了满足所有的约束,我们可以令λ1=1lnσ2πλ2=0λ3=12σ2,从而得到

p(x)=1(2πσ2)12exp{(xμ)22σ2}

因此最大化微分熵的分布是高斯分布N(xμ,σ2)。这也是当我们不知道真实的分布时,往往使用高斯分布的一个原因。因为高斯分布拥有最大的熵,我们通过这个假定来保证了最小可能量的结构。

X[a,b],无其它约束条件,则此时最大熵分布就是该区间上的均匀分布,这里可以和X为离散随机变量的情形联系起来。

如果我们求高斯分布的微分熵,我们会得到:

H(p)=p(x)lnp(x)dx=1(2πσ2)12exp{(xμ)22σ2}[ln1(2πσ2)12(xμ)22σ2]=ln1(2πσ2)12+12σ21(2πσ2)12exp{(xμ)22σ2}(xμ)2dxσ2=12{ln(2πσ2)+1}

我们可以看到熵随着分布宽度(即σ2)的增加而增加。这个结果也表明,与离散熵不同,微分熵可以为负,因为对于上式而言,当σ2<12πe时,H(p)<0

1.2 做为随机变量和分布的高斯分布

当我们考虑多个随机变量之和的时候,也会产生高斯分布。根据林德伯格所证明的中心极限定理(central limit theorem)(为我们在博客《概率论沉思录:初等假设检验》中提到过的伯努利试验的棣莫弗-拉普拉斯极限定理的推广),设{Xk}k=1N是相互独立且具有共同分布的随机变量序列,假定μ=E[Xk]σ2=Var[Xk]都存在,并令SN=X1++XN,则当N时,有SNN(Nμ,σN)(依分布)。

例如,考虑N个随机变量X1,,XN,每一个都是区间[0,1]上的均匀分布,然后考虑N个随机变量的均值SNN=1N(X1++XN)的分布。对于大的N,这个分布趋向于高斯分布,如下图所示:

在实际应用中,随着N的增加,分布会很迅速收敛为高斯分布。

2 高斯分布的性质

2.1 高斯分布的几何性质

观察式(1)中的多元高斯分布的形式,考虑其中在指数位置上出现的二次型

(2)(xμ)TΣ1(xμ)

由于协方差矩阵Σ是对称矩阵,那么Σ1也是对称矩阵(注意Σ1的第j,i个元素为(1)i+j[detΣij/detΣ],其中Σij是从Σ中除去第i行和第j列后得到的矩阵,而对于对称矩阵ΣdetΣij=detΣji)。我们假定Σ是正定的,那么Σ1也是正定的(后面我们会证明)。于是,式(2)xμ马⽒距离(Mahalanobis distance)Δ的平方。当Σ是单位阵时,就变成了欧氏距离。

给定D×D的对称正定矩阵A,则从点x到原点的马氏距离Δ由正定二次型0<Δ2=xTAxx0)确定;反过来,一个正定二次型xTAx>0x0)可以解释为点x到原点的马氏距离Δ的平方[7]。此外,从点x到任意一固定点μ的马氏距离Δ的平方由公式(xμ)TA(xμ)给出。

马氏距离可以理解为在新的基底(也即A的规范正交特征向量组成的基底)下两点之间的距离。例如,设D=2,则到原点的马氏距离为常数c的点x=(x1,x2)T满足xTAx=c2。对A进行谱分解得:

A=λ1e1e1T+λ2e2e2T

所以

xTAx=xT(λ1e1e1T+λ2e2e2T)x=λ1(e1Tx)2+λ2(e2Tx)2

y1=e1Tx,y2=e2Tx,则y=(y1,y2)T可视为点x=(x1,x2)T在新的基底e1,e2下的坐标,基变换关系为(y1y2)=(e1e2)T(x1x2)=(e1Te2T)(x1x2)。此时,y=(y1,y2)T在以原点为中心的椭圆上,其方程为:

λ1y12+λ2y22=c2

该椭圆的轴即为做为新基底的规范正交特征向量e1,e2。我们可以发现将x=cλ11/2e1代入椭圆方程得到xTAx=λ1(cλ11/2e1Te1)2=c2,可见椭圆沿e1方向的半轴长为cλ11/2;同理,x=cλ21/2e2也给出了椭圆沿e2方向的半轴长cλ21/2。到原点的马氏距离为常数c的点的位置如下图所示(D=2,λ1<λ2):

如果D>2,到原点的马氏距离为常数c=xTAx的点x=(x1,,xD)在椭球面λ1(e1Tx)2++λD(eDTxT)2=c2上,其轴由A的规范正交特征向量给出。沿ei方向的半轴长为cλi,i=1,2,,p,其中λ1,,λDA的特征值。

下图直观地说明了马氏距离相比欧氏距离具有优越性的情况:

该图描述了重心(样本均值)在点Q的一组点。PQ的欧氏距离大于OQ的欧氏距离。然而,P点却比原点O更像是属于这一组点内的点。如果我们用马氏距离来度量距离,则P距离Q就会比O距离Q要近了。

至于马氏距离中的矩阵A如何确定则和样本各维度的标准差有关。马氏距离相当于将样本各维度除以样本标准差之后(也即“标准化”之后),再采用欧氏距离的公式进行计算,感兴趣的读者可以阅读《Applied Multivariate Statistical Analysis》1.5节。

如果令多元高斯分布中的马氏距离的平方Δ2=(xμ)TΣ1(xμ)为常数c2,那么多元高斯分布的概率密度也为常数,这也就意味着到均值μ的马氏距离相等的点拥有相同的赋概。我们前面提到过,这些点在一个椭球面上,我们将其称为轮廓线(contours)

常数概率密度轮廓线={x(xμ)TΣ1(xμ)=c2}=中心在μ

为了确定椭球面的轴方向和半轴长,我们对协方差矩阵的逆矩阵Σ1进行谱分解:

(3)Σ1=i=1D1λieieiT

其中λ1,,λDΣ的特征值,e1,,eD为与之相伴的规范正交特征向量。

关于Σ1,我们有下列结论:
结论Σ是正定的,其逆为Σ1,则Σe=λeΣ1e=(1λ)e

所以Σ的特征值-特征向量对(λ,e)对应于Σ1的特征值-特征向量对(1λ,e)。且Σ1也是正定的。

证明 对于正定的Σ以及一个特征向量e0,我们有0<eTΣe=eT(Σe)=eT(λe)=λeTe=λ。而且e=Σ1(Σe)=Σ1(λe)=λΣ1e,用λ>0除,得到Σ1e=(1λ)e。于是,(1λ,e)Σ1的一对特征值-特征向量。因此,我们对任意D维向量x0

xTΣ1x=xT(i=1D1λieieiT)x=i=1D(1λi)(eiTx)2>0

由此得出Σ1是正定的。

将式(3)Σ1的谱分解结果代入式(2)所示的二次型中,有

c2=(xμ)TΣ1(xμ)=(xμ)T(i=1D1λieieiT)(xμ)=i=1D(eiT(xμ))2λi

yi=eiT(xμ),则y可视为点xμ在新的基底{ei}下的坐标,基变换关系为(y1yD)=(e1eD)T(xμ)=(e1TeDT)(xμ)。此时,y=(y1,,yD)T在以μ为中心的椭球面上,其方程为:

i=1Dyi2λi=c2

该椭球面的轴即为做为新基底的规范正交特征向量e1,,eD。我们可以发现将xi=μ+cλi1/2ei代入椭圆方程得到(xμ)TΣ1(xμ)=(cλi1/2eiTei)2λi=c2,可见椭圆沿ei方向的半轴长为cλi1/2i=1,,D)。在D=2的情况下,该椭球面为中心为μ的椭圆,如下图所示(λ2<λ1):

现在我们考虑多元高斯分布是否是归一化的。我们之前使用了变量替换y=T(x)=(e1TeDT)(xμ),于是在新的基底下的积分公式可以经由变量替换[8][9]表示为

N(y0,Σ)dy=N(T(x)0,Σ)|detJ|dx

这里J为仿射变换T的Jacobian矩阵:

J=(T1(x)x1T1(x)xDTD(x)x1TD(x)xD)=(e11e1DeD1eDD)=(e1TeDT)

U=(e1eD)(此时J=(e1TeDT)=UT),由于新的基底e1,,eD是规范正交的,因此U是正交矩阵,此时我们有:

|detJ|=|detUT|=(detUT)2=detUTdetU=det(UTU)=detI=1

而我们发现又多元高斯分布N(y0,Σ)可以分解成D个独立一元高斯分布N(yi0,λi)的乘积:

N(yμ,Σ)=1(2π)D21(detΣ)12exp{12i=1Dyi2λi}=i=1D1(2πλi)12exp{yi22λi}=i=1DN(yi0,λi)

(其中我们用到了结论detΣ=iλi)因此

N(y0,Σ)dy=i=1DN(yi0,λi)dyi=1

于是我们有

N(xμ,Σ)dx=N(T(x)0,Σ)dx=N(y0,Σ)dy=1

至此我们证明了多元高斯分布是归一化的。

2.2 高斯分布的一阶矩和二阶矩

现在我们考察多元高斯分布的一阶矩和二阶矩,这可以提供参数μΣ的描述。多元高斯分布下x的期望为:

E[x]=1(2π)D21(detΣ)12exp{12(xμ)TΣ1(xμ)}xdx=1(2π)D21(detΣ)12exp{12zTΣ1z}(z+μ)dz=1(2π)D21(detΣ)12exp{12zTΣ1z}zdz+μ

其中在第二个等式中我们使用了z=xμ进行了变量替换,在第三个等式中我们利用了N(z0,Σ)dz=1。现在我们来考虑第一个积分项,我们注意到指数位置是z的偶函数,而z是奇函数,因此被积函数是奇函数[10],又由于积分区域关于原点对称,因此第一个积分项为0。于是我们有

E[x]=μ

因此我们把μ称为多元高斯分布的均值。

现在我们考虑多元高斯分布的二阶矩。在一元变量的情形下,二阶矩由E[x2]给出。对于多元高斯分布,有D2个由E[xixj]给出的二阶矩,可以聚集在一起组成矩阵E[xxT]。这个矩阵可以表示为:

E[xxT]=1(2π)D21(detΣ)12exp{12(xμ)TΣ1(xμ)}xxTdx=1(2π)D21(detΣ)12exp{12zTΣ1z}(z+μ)(z+μ)Tdz

其中在第二个等式中我们再次使用了z=xμ来进行变量替换,涉及到zμTμzT的交叉项将再次变为0,而μμT也可以拿出。于是我们有

E[xxT]=1(2π)D21(detΣ)12exp{12zTΣ1z}zzTdz+μμT

接下来考虑第一个积分项,我们采用和之前类似的做法,对Σ1进行谱分解得Σ1=i1λieieiT,并使用变量替换令yi=eiTz(基变换关系为(y1yD)=(e1TeDT)z,因此z=(e1eD)y=iyiei),于是有

1(2π)D21(detΣ)12exp{12zTΣ1z}zzTdz=1(2π)D21(detΣ)12i=1Dj=1DeiejTexp{12k=1Dyk2λk}yiyjdy=1(2π)D21(detΣ)12i=1DeieiTexp{12k=1Dyk2λk}yi2dy=i=1DeieiTλi=Σ

其中第二个等式是由于当ij时,被积函数是奇函数,导致积分项为0;第三个等式是由于E[yi2]=N(yi0,λi)yi2dyi=λi+02=λi;最后一个等式是根据Σ的谱分解得到。

这样,我们就得到了

E[xxT]=Σ+μμT

对于一元随机变量的方差,为了定义方差,我们在取二阶矩之前会减掉均值。类似地,对于多元变量的情形,把均值减掉同样很方便。这给出了随机变量x协方差(covariance) 的定义:

Cov[x]=E[(xE[x])(xE[x])T]

对于多元高斯分布这一特例,我们有

Cov[x]=E[xxTxμTμxT+μμT]=E[xxT]μμT=Σ

由于参数Σ控制了多元高斯分布下x的协方差,因此它被称为协方差矩阵。

虽然式(1)定义的高斯分布N(xμ,Σ)被广泛用作概率密度模型,但是它有着一些巨大的局限性。考虑分布中自由参数的数量。一个通常的对称协方差矩阵Σ1+2++D=D(D+1)2个独立参数,μ中有另外D个独立参数,因此总计有D(D+1)2+D=D(D+3)2个独立参数。对于大的D值,参数的总数随着D以平方的方式增长,导致对大的矩阵进行操作(如求逆)的计算变得不可行。解决这个问题的一种方式是使用协方差矩阵的限制形式。如果我们考虑对角的协方差矩阵,即Σ=diag(σi2),那么在概率密度模型中,我们就有总数2D个独立参数。此时常数概率密度轮廓线为与轴对齐的椭球。我们可以进一步地把协方差矩阵限制成正比于单位矩阵,也即Σ=σ2I,此时它被称为各向同性(isotropic)的协方差。这使得模型有D+1个独立的参数,并且常数概率密度轮廓线为球面。下图展示了通常的协方差矩阵、对角的协方差矩阵以及各向同性协方差矩阵的概率密度轮廓线(D=2):

高斯分布的另一个局限性是它本质上是单峰的(即只有一个最大值),因此不能很好地近似多峰分布。不过,相当多的多峰分布可以使用混合高斯分布来描述(参见博客《统计学习:EM算法及其在高斯混合模型(GMM)中的应用》)。

参考

  • [1] Bishop C M, Nasrabadi N M. Pattern recognition and machine learning[M]. New York: springer, 2006.
  • [2] Schroeder D V. An introduction to thermal physics[M]. Oxford University Press, 2020.
  • [3] 《维基百科:熵 (统计物理学)》
  • [4] Feller W. An introduction to probability theory and its applications, Volume 1[M]. John Wiley & Sons, 1991.
  • [5] Weisstein E W. CRC concise encyclopedia of mathematics[M]. Chapman and Hall/CRC, 2002.
  • [6] Bengio Y, Goodfellow I, Courville A. Deep learning[M]. Cambridge, MA, USA: MIT press, 2017.
  • [7] Johnson R A, Wichern D W. Applied multivariate statistical analysis[J]. 2002.
  • [8] Rudin W. Principles of mathematical analysis[M]. New York: McGraw-hill, 1964.
  • [9] Axler S. Linear algebra done right[M]. springer publication, 2015.
  • [10] 维基百科:《奇函数与偶函数》)

__EOF__

  • 本文作者: 猎户座
  • 本文链接: https://www.cnblogs.com/orion-orion/p/18688763
  • 关于博主: 研究生小菜一枚,机器学习半吊子,并行计算混子。
  • 版权声明: 欢迎您对我的文章进行转载,但请务必保留原始出处哦(*^▽^*)。
  • 声援博主: 如果您觉得文章对您有帮助,可以点击文章右下角推荐一下。
  • posted @   orion-orion  阅读(415)  评论(0编辑  收藏  举报
    相关博文:
    阅读排行:
    · DeepSeek-R1本地部署如何选择适合你的版本?看这里
    · 开源的 DeepSeek-R1「GitHub 热点速览」
    · 传国玉玺易主,ai.com竟然跳转到国产AI
    · 揭秘 Sdcb Chats 如何解析 DeepSeek-R1 思维链
    · 自己如何在本地电脑从零搭建DeepSeek!手把手教学,快来看看! (建议收藏)
    历史上的今天:
    2023-01-23 SICP:复数的直角和极坐标的表示(Python实现)
    点击右上角即可分享
    微信分享提示