机器学习（十）—

http://antkillerfarm.github.io/

对Σ\Sigma的限制（续）

这实际上也就是方法一中对角线元素的均值，反映到二维高斯分布图上就是椭圆变成圆。

当我们要估计出完整的Σ\Sigma时，我们需要m≥n+1m\ge n+1才能保证在最大似然估计下得出的Σ\Sigma是非奇异的。然而在上面的任何一种假设限定条件下，只要m≥2m\ge 2就可以估计出限定的Σ\Sigma。

这样做的缺点也是显而易见的，我们认为特征间相互独立，这个假设太强。接下来，我们给出一种称为因子分析（factor analysis）的方法，使用更多的参数来分析特征间的关系，并且不需要计算一个完整的Σ\Sigma。

利用多元高斯分布密度函数计算积分的技巧

I(A,b,c)=∫xexp(−12(xTAx+xTb+c))dx

I(A,b,c)=\int_x\exp\left(-\frac{1}{2}(x^TAx+x^Tb+c)\right)\mathrm{d}x

其中A为对称正定矩阵，b为向量。对于上面这样的积分，可以使用“完全配方法”（completion-of-squares）的数学技巧求解。

因为

xTAx+xTb+c=(x−h)TA(x−h)+k

x^TAx+x^Tb+c=(x-h)^TA(x-h)+k

其中h=−A−1b2,k=c−bTA−1b4h=-\frac{A^{-1}b}{2},k=c-\frac{b^TA^{-1}b}{4}。

所以

I(A,b,c)=∫xexp(−12((x−h)TA(x−h)+k))dx=∫xexp(−12(x−h)TA(x−h)−k/2)dx=exp(−k/2)⋅∫xexp(−12(x−h)TA(x−h))dx

\begin{align} I(A,b,c)&=\int_x\exp\left(-\frac{1}{2}((x - h)^TA(x - h)+k)\right)\mathrm{d}x \\&=\int_x\exp\left(-\frac{1}{2}(x - h)^TA(x - h)-k/2\right)\mathrm{d}x \\&=\exp(-k/2)\cdot\int_x\exp\left(-\frac{1}{2}(x - h)^TA(x - h)\right)\mathrm{d}x \end{align}

令μ=h,Σ=A−1\mu=h,\Sigma=A^{-1}，则：

I(A,b,c)=(2π)n/2∣Σ∣1/2exp(k/2)⋅∫x1(2π)n/2∣Σ∣1/2exp(−12(x−μ)TΣ−1(x−μ))dx

I(A,b,c)=\frac{(2\pi)^{n/2}\lvert\Sigma\rvert^{1/2}}{\exp(k/2)}\cdot\int_x\frac{1}{(2\pi)^{n/2}\lvert\Sigma\rvert^{1/2}}\exp\left(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\right)\mathrm{d}x

公式右侧的被积分函数，正好是多元高斯分布密度函数，因此该积分值为1。于是：

I(A,b,c)=(2π)n/2∣Σ∣1/2exp(k/2)

I(A,b,c)=\frac{(2\pi)^{n/2}\lvert\Sigma\rvert^{1/2}}{\exp(k/2)}

注：原始讲义里，Chuong B. Do写的《Gaussian processes》的附录A.1和本节内容类似，但推导过程有问题，疑似笔误，特更换为维基百科中的例子。（矩阵的完全配方那块的变换，我能推导出维基百科的结果，但推导不出Chuong B. Do的结果。）如有错误，望读者指出。

边缘和条件高斯分布

假设x由两个随机向量组成（可以看作是将之前的x(i)x^{(i)}分成了两部分）。

x=[x1x2]

x=\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}

其中x1∈Rr,x1∈Rsx_1\in R^r,x_1\in R^s，则x实际上是r+sr+s维向量。

假设x∼N(μ,Σ)x\sim N(\mu,\Sigma)，其中：

μ=[μ1μ2],Σ=[Σ11Σ21Σ12Σ22]

\mu=\begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix},\Sigma=\begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}

因为协方差矩阵是对称矩阵，因此Σ12=ΣT21\Sigma_{12}=\Sigma_{21}^T。

Cov(x)=Σ=[Σ11Σ21Σ12Σ22]=E[(x−μ)(x−μ)T]=E[(x1−μ1x2−μ2)(x1−μ1x2−μ2)]=E[(x1−μ1)(x1−μ1)T(x2−μ2)(x1−μ1)T(x1−μ1)(x2−μ2)T(x2−μ2)(x2−μ2)T]

\begin{align} Cov(x)&=\Sigma=\begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix} \\&=E\begin{bmatrix} (x-\mu)(x-\mu)^T \end{bmatrix}=E\begin{bmatrix} \begin{pmatrix} x_1-\mu_1 \\ x_2-\mu_2 \end{pmatrix} & \begin{pmatrix} x_1-\mu_1 \\ x_2-\mu_2 \end{pmatrix} \end{bmatrix} \\&=E\begin{bmatrix} (x_1-\mu_1)(x_1-\mu_1)^T & (x_1-\mu_1)(x_2-\mu_2)^T \\ (x_2-\mu_2)(x_1-\mu_1)^T & (x_2-\mu_2)(x_2-\mu_2)^T \end{bmatrix} \end{align}

因此，E[x1]=μ1,Cov(x1)=E[(x1−μ1)(x1−μ1)T]=Σ11E[x_1]=\mu_1,Cov(x_1)=E[(x_1-\mu_1)(x_1-\mu_1)^T]=\Sigma_{11}。可见，多元高斯分布的边缘分布仍然是多元高斯分布。

下面讨论一下条件高斯分布。

p(x1|x2)=p(x1,x2)p(x2)=1(2π)n/2∣Σ∣1/2exp(−12(x−μ)TΣ−1(x−μ))∫x1p(x1,x2;μ,Σ)dx1=1Z1exp⎧⎩⎨−12([x1x2]−[μ1μ2])T[V11V21V12V22]([x1x2]−[μ1μ2])⎫⎭⎬

\begin{align} p(x_1\vert x_2)&=\frac{p(x_1,x_2)}{p(x_2)}=\frac{\frac{1}{(2\pi)^{n/2}\lvert\Sigma\rvert^{1/2}}\exp\left(-\frac{1}{2}(x-\mu)^T\Sigma^{-1}(x-\mu)\right)}{\int_{x_1}p(x_1,x_2;\mu,\Sigma)\mathrm{d}x_1} \\&=\frac{1}{Z_1}\exp\left\{-\frac{1}{2}\left(\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}-\begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}\right)^T\begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \end{bmatrix}\left(\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}-\begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}\right)\right\} \end{align}

其中的Z1Z_1是和x1x_1无关的部分，可看作常数，下面的ZiZ_i也是同理。

Σ−1=V=[V11V21V12V22]

\Sigma^{-1}=V=\begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \end{bmatrix}

因为：

([x1x2]−[μ1μ2])T[V11V21V12V22]([x1x2]−[μ1μ2])=(x1−μ1)TV11(x1−μ1)+(x1−μ1)TV12(x2−μ2)+(x2−μ2)TV21(x1−μ1)+(x2−μ2)TV22(x2−μ2)

\begin{align} &\left(\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}-\begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}\right)^T\begin{bmatrix} V_{11} & V_{12} \\ V_{21} & V_{22} \end{bmatrix}\left(\begin{bmatrix} x_1 \\ x_2 \end{bmatrix}-\begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}\right) \\&=(x_1-\mu_1)^TV_{11}(x_1-\mu_1)+(x_1-\mu_1)^TV_{12}(x_2-\mu_2) \\&\qquad+(x_2-\mu_2)^TV_{21}(x_1-\mu_1)+(x_2-\mu_2)^TV_{22}(x_2-\mu_2) \end{align}

保留上式中与x1x_1有关的部分，可得：

p(x1|x2)=1Z2exp(−12(xT1V11x1−2xT1V11μ1+2xT1V12(x2−μ2)))

p(x_1\vert x_2)=\frac{1}{Z_2}\exp\left(-\frac{1}{2}\left(x_1^TV_{11}x_1-2x_1^TV_{11}\mu_1+2x_1^TV_{12}(x_2-\mu_2)\right)\right)

使用上一节中的完全配方技巧，可得：

p(x1|x2)=1Z3exp(−12(x1−μ1|2)TV11(x1−μ1|2))

p(x_1\vert x_2)=\frac{1}{Z_3}\exp\left(-\frac{1}{2}(x_1-\mu_{1\vert 2})^TV_{11}(x_1-\mu_{1\vert 2})\right)

其中：

μ1|2=μ1−V−111V12(x2−μ2)(1)

\mu_{1\vert 2}=\mu_1-V_{11}^{-1}V_{12}(x_2-\mu_2)\tag{1}

即：

x1|x2∼N(μ1−V−111V12(x2−μ2),V−111)

x_1\vert x_2\sim N(\mu_1-V_{11}^{-1}V_{12}(x_2-\mu_2),V_{11}^{-1})

另，根据分块矩阵的求逆法则，可得：

Σ−1=[Σ11Σ21Σ12Σ22]−1=[(Σ11−Σ12Σ−122Σ21)−1−Σ−122Σ21(Σ11−Σ12Σ−122Σ21)−1−(Σ11−Σ12Σ−122Σ21)−1Σ12Σ−122(Σ22−Σ21Σ−111Σ12)−1]

\Sigma^{-1}=\begin{bmatrix} \Sigma_{11} & \Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}^{-1}=\begin{bmatrix} (\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21})^{-1} & -(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21})^{-1}\Sigma_{12}\Sigma_{22}^{-1} \\ -\Sigma_{22}^{-1}\Sigma_{21}(\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21})^{-1} & (\Sigma_{22}-\Sigma_{21}\Sigma_{11}^{-1}\Sigma_{12})^{-1} \end{bmatrix}

因此：

Σ1|2=V−111=Σ11−Σ12Σ−122Σ21(2)

\Sigma_{1\vert 2}=V_{11}^{-1}=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}\tag{2}

因子分析的例子

下面通过一个简单例子，来引出因子分析背后的思想。

假设我们有m=5个2维的样本点xix^{i}，如下：

按照因子分析模型，样本点的生成过程如下：

1.我们首先认为在1维空间（这里k=1），存在着按正态分布生成的m个点z(i)z^{(i)}，即：

z(i)∼N(0,I)

z^{(i)}\sim N(0,I)

这里的I是单位矩阵。

2.使用变换矩阵Λ∈Rn×k\Lambda\in R^{n\times k}，将z(i)z^{(i)}映射到n维空间中，即Λz(i)\Lambda z^{(i)}。

3.使用n维向量μ\mu，将Λz(i)\Lambda z^{(i)}移动到样本的中心点μ\mu，即μ+Λz(i)\mu+\Lambda z^{(i)}

4.样本点不可能这么规则，在模型上会有一定偏差，因此我们需要将上步生成的点做一些扰动（误差）。这里添加一个n维的扰动向量ϵ∼N(0,Ψ)\epsilon \sim N(0,\Psi)。

综上可得:

x(i)=μ+Λz(i)+ϵ

x^{(i)}=\mu+\Lambda z^{(i)}+\epsilon

x|z∼N(μ+Λz,Ψ)

x \vert z\sim N(\mu+\Lambda z,\Psi)

由以上的直观分析，我们知道了因子分析其实就是认为：高维样本点实际上是由低维样本点经过高斯分布、线性变换、误差扰动生成的，因此高维数据可以使用低维来表示。

线性回归的概率模型

在进一步讨论因子分析模型之前，我们首先讨论一下，和它类似的线性回归的概率模型。

从概率的角度看，线性回归中的y(i)y^{(i)}可以看作是预测函数hθ(x)h_\theta(x)加上扰动后的结果。即：

y(i)=θTx(i)+ϵ(i),ϵ(i)∼N(0,σ2)

y^{(i)}=\theta^Tx^{(i)}+\epsilon^{(i)},\epsilon^{(i)}\sim N(0,\sigma^2)

p(ϵ(i))=12π−−√σexp(−(ϵ(i))22σ2)

p(\epsilon^{(i)})=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(\epsilon^{(i)})^2}{2\sigma^2}\right)

p(y(i)|x(i);θ)=12π−−√σexp(−(y(i)−θTx(i))22σ2)

p(y^{(i)}\vert x^{(i)};\theta)=\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2}\right)

ℓ(θ)=log∏i=1m12π−−√σexp(−(y(i)−θTx(i))22σ2)=∑i=1mlog12π−−√σexp(−(y(i)−θTx(i))22σ2)=mlog12π−−√σ−1σ2⋅12∑i=1m(y(i)−θTx(i))2

\begin{align} \ell(\theta)&=\log\prod_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2}\right)=\sum_{i=1}^m\log\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(y^{(i)}-\theta^Tx^{(i)})^2}{2\sigma^2}\right) \\&=m\log\frac{1}{\sqrt{2\pi}\sigma}-\frac{1}{\sigma^2}\cdot\frac{1}{2}\sum_{i=1}^m\left(y^{(i)}-\theta^Tx^{(i)}\right)^2 \end{align}

从上式可以看出采用极大似然估计和采用代价函数J(θ)J(\theta)的效果是一样的。其中：

J(θ)=12∑i=1m(y(i)−θTx(i))2

J(\theta)=\frac{1}{2}\sum_{i=1}^m\left(y^{(i)}-\theta^Tx^{(i)}\right)^2

因子分析模型

假设z和x的联合分布为：

[zx]∼N(μzx,Σ)

\begin{bmatrix} z \\ x \end{bmatrix}\sim N(\mu_{zx},\Sigma)

我们的任务就是求出μzx\mu_{zx}和Σ\Sigma。

因为：

E[x]=E[μ+Λz+ϵ]=μ+ΛE[z]+E[ϵ]=μ

E[x]=E[\mu+\Lambda z+\epsilon]=\mu+\Lambda E[z]+E[\epsilon]=\mu

所以：

μzx=[0⃗ μ]

\mu_{zx}=\begin{bmatrix} \vec{0} \\ \mu \end{bmatrix}

因为：

Σ=[ΣzzΣxzΣzxΣxx]

\Sigma=\begin{bmatrix} \Sigma_{zz} & \Sigma_{zx} \\ \Sigma_{xz} & \Sigma_{xx} \end{bmatrix}

所以我们只要分别计算这四个值即可。

因为z∼N(0,I)z\sim N(0,I)，所以Σzz=I\Sigma_{zz}=I。

Σzx=E[(z−E[z])(x−E[x])T]=E[(z−0)(μ+Λz+ϵ−μ)T]=E[z(Λz+ϵ)T]=E[z(Λz)T+zϵT]=E[zzTΛT+zϵT]=E[zzT]ΛT+E[zϵT]

\begin{align} \Sigma_{zx}&=E[(z-E[z])(x-E[x])^T]=E[(z-0)(\mu+\Lambda z+\epsilon-\mu)^T]=E[z(\Lambda z+\epsilon)^T] \\&=E[z(\Lambda z)^T+z\epsilon^T]=E[zz^T\Lambda^T+z\epsilon^T]=E[zz^T]\Lambda^T+E[z\epsilon^T] \end{align}

因为z和ϵ\epsilon是相互独立的随机变量，因此E[zϵT]=E[z]E[ϵT]=0E[z\epsilon^T]=E[z]E[\epsilon^T]=0。

又因为E[zzT]=Cov(z)=IE[zz^T]=Cov(z)=I，所以Σzx=ΛT\Sigma_{zx}=\Lambda^T。

Σxx=E[(x−E[x])(x−E[x])T]=E[(μ+Λz+ϵ−μ)(μ+Λz+ϵ−μ)T]=E[(Λz+ϵ)(ΛzT+ϵT)]=E[Λz(Λz)T+ϵ(Λz)T+ΛzϵT+ϵϵT]=E[ΛzzTΛT+ϵzTΛT+ΛzϵT+ϵϵT]=ΛE[zzT]ΛT+E[ϵzT]ΛT+ΛE[zϵT]+E[ϵϵT]=ΛIΛT+0+0+Ψ=ΛΛT+Ψ

\begin{align} \Sigma_{xx}&=E[(x-E[x])(x-E[x])^T]=E[(\mu+\Lambda z+\epsilon-\mu)(\mu+\Lambda z+\epsilon-\mu)^T] \\&=E[(\Lambda z+\epsilon)(\Lambda z^T+\epsilon^T)]=E[\Lambda z(\Lambda z)^T+\epsilon(\Lambda z)^T+\Lambda z\epsilon^T+\epsilon\epsilon^T] \\&=E[\Lambda zz^T\Lambda^T+\epsilon z^T\Lambda^T+\Lambda z\epsilon^T+\epsilon\epsilon^T] \\&=\Lambda E[zz^T]\Lambda^T+E[\epsilon z^T]\Lambda^T+\Lambda E[z\epsilon^T]+E[\epsilon\epsilon^T] \\&=\Lambda I\Lambda^T+0+0+\Psi=\Lambda \Lambda^T+\Psi \end{align}

把这些结果合在一起，可得：

[zx]∼N([0⃗ μ],[IΛΛTΛΛT+Ψ])(3)

\begin{bmatrix} z \\ x \end{bmatrix}\sim N\left(\begin{bmatrix} \vec{0} \\ \mu \end{bmatrix},\begin{bmatrix} I & \Lambda^T \\ \Lambda & \Lambda \Lambda^T+\Psi \end{bmatrix}\right)\tag{3}

从这个结论可以看出：x∼N(μ,ΛΛT+Ψ)x\sim N(\mu,\Lambda \Lambda^T+\Psi)

因此它的对数似然函数为：

ℓ(μ,Λ,Ψ)=log∏i=1m1(2π)n/2∣ΛΛT+Ψ∣1/2exp(−12(x(i)−μ)T(ΛΛT+Ψ)−1(x(i)−μ))

\ell(\mu,\Lambda,\Psi)=\log\prod_{i=1}^m\frac{1}{(2\pi)^{n/2}\lvert\Lambda \Lambda^T+\Psi\rvert^{1/2}}\exp\left(-\frac{1}{2}(x^{(i)}-\mu)^T(\Lambda \Lambda^T+\Psi)^{-1}(x^{(i)}-\mu)\right)

但这个函数是很难最大化的，需要使用EM算法解决之。

因子分析的EM估计

E-step比较简单。由公式1、2、3，可得：

μz(i)|x(i)=ΛT(ΛΛT+Ψ)−1(x(i)−μ)

\mu_{z^{(i)}\vert x^{(i)}}=\Lambda^T(\Lambda \Lambda^T+\Psi)^{-1}(x^{(i)}-\mu)

Σz(i)|x(i)=I−ΛT(ΛΛT+Ψ)−1Λ

\Sigma_{z^{(i)}\vert x^{(i)}}=I-\Lambda^T(\Lambda \Lambda^T+\Psi)^{-1}\Lambda

因此：

Qi(z(i))=1(2π)n/2∣Σz(i)|x(i)∣1/2exp(−12(x(i)−μz(i)|x(i))TΣ−1z(i)|x(i)(x(i)−μz(i)|x(i)))

Q_i(z^{(i)})=\frac{1}{(2\pi)^{n/2}\lvert\Sigma_{z^{(i)}\vert x^{(i)}}\rvert^{1/2}}\exp\left(-\frac{1}{2}(x^{(i)}-\mu_{z^{(i)}\vert x^{(i)}})^T\Sigma_{z^{(i)}\vert x^{(i)}}^{-1}(x^{(i)}-\mu_{z^{(i)}\vert x^{(i)}})\right)

M-step的最大化的目标是：

∑i=1m∫z(i)Qi(z(i))logp(x(i),z(i);μ,Λ,Ψ)Qi(z(i))dz(i)

\sum_{i=1}^m\int_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)},z^{(i)};\mu,\Lambda,\Psi)}{Q_i(z^{(i)})}\mathrm{d}z^{(i)}

下面我们重点求Λ\Lambda的估计公式。

首先将上式简化为:

∑i=1m∫z(i)Qi(z(i))logp(x(i)|z(i);μ,Λ,Ψ)p(z(i))Qi(z(i))dz(i)=∑i=1m∫z(i)Qi(z(i))[logp(x(i)|z(i);μ,Λ,Ψ)+logp(z(i))−logQi(z(i))]dz(i)=∑i=1mEz(i)∼Qi[logp(x(i)|z(i);μ,Λ,Ψ)+logp(z(i))−logQi(z(i))]

\begin{align} &\sum_{i=1}^m\int_{z^{(i)}}Q_i(z^{(i)})\log\frac{p(x^{(i)}\vert z^{(i)};\mu,\Lambda,\Psi)p(z^{(i)})}{Q_i(z^{(i)})}\mathrm{d}z^{(i)} \\&=\sum_{i=1}^m\int_{z^{(i)}}Q_i(z^{(i)})\left[\log p(x^{(i)}\vert z^{(i)};\mu,\Lambda,\Psi)+\log p(z^{(i)})-\log Q_i(z^{(i)})\right]\mathrm{d}z^{(i)} \\&=\sum_{i=1}^m E_{z^{(i)}\sim Q_i}\left[\log p(x^{(i)}\vert z^{(i)};\mu,\Lambda,\Psi)+\log p(z^{(i)})-\log Q_i(z^{(i)})\right] \end{align}