机器学习 [白板推导]

开个文章补一下白版推导吧,现在回过头来看这些都还是挺有用的,不把数学底层逻辑理清总有一种学了又好像没学的感觉。


1. 绪论

1.1. 机器学习的两大流派——频率派和贝叶斯派

1.1.1. 频率派:统计机器学习

  • 建模方法:概率派将模型参数 \theta 视为一个常量,样本 X_{N\times p} 是一个随机变量,使用极大似然估计(Maximum Likelihood Estimate,MLE)方法建模求得 \theta 的解,即 

\theta _{MLE}=argmax_{\left \{ \theta \right \}}\log[P(X|\theta )]

1.1.2. 贝叶斯派:概率图模型

  • 建模方法:认为模型参数 \theta 是一个随机变量,其遵从的分布 p(\theta )被称为先验分布(没见到样本就有的概率分布),根据贝叶斯公式计算后验概率 

 p(\theta |X)=\frac{p(X|\theta )\cdot p(\theta )}{p(X)}=\frac{p(X|\theta )\cdot p(\theta )}{\int_{\theta }p(X|\theta )\cdot p(\theta )d\theta } 

        ,则寻找最大后验估计(Maximum a Posteriori estimation,MAP),即

\begin{aligned} \theta _{MAP} &=argmax_{\left \{ \theta \right \}}p(\theta |X)\\ &=argmax_{\left \{ \theta \right \}}p(X|\theta )\cdot p(\theta )\\ & =argmax_{\left \{ \theta \right \}}[\log p(X|\theta ) + \log p(\theta )] \end{aligned}

        ,和MLP相比仅多了一项 \log p(\theta )

  • 求解后验概率时,分母的积分计算是复杂的,因此通常采用蒙特卡洛采样法来代替;
  • 当先验分布遵从高斯分布时,

p(\theta )=a\cdot e^{\frac{-\theta ^2}{2\sigma ^2}}

        ,a为常数,因此 

\log p(\theta )=\log a+(-\frac{\theta ^2}{2\sigma ^2})

        ,相当于在MLP优化函数中增加一项L2正则化

  • 贝叶斯预测:已有数据为 X ,预测新数据 \widetilde{x} 时,需要计算新的后验概率

p(\widetilde{x}|X)=\int_\theta p(\widetilde{x}|\theta )\cdot p(\theta |X)d\theta

1.2. 数学基础

1.2.1. 矩阵迹的性质

tr(AB)=tr(BA)

tr(ABC)=tr(CAB)=tr(BCA)

1.2.2.  克罗内克积(Kronecker Product)

 A_{m\times n}\bigotimes B_{p\times q} = \begin{bmatrix} a_{11}B & \cdots & a_{1n}B\\ \vdots &\ddots & \vdots\\ a_{m1}B & \cdots & a_{mn}B \end{bmatrix}_{mp\times nq}

A的元素保持相对位置分别乘以B

1.2.3. 奇异值分解(Singular Value Decomposition,SVD)

  • 类似于特征值分解,奇异值分解也是将矩阵分解为需要的矩阵,而特征值分解只能处理方针,奇异值分解则是针对非方针提出;
  • 定义:对于一个矩阵A_{m\times n},存在两个酉矩阵U_{m\times m}V_{n\times n}U^TU=I_mV^TV=I_{n},使得A=U\Sigma V^T,其中\Sigma_{m\times n}主对角线上元素为奇异值,其他位置元素为0;
  • 求解:对A^TA进行特征值分解得到的特征向量经过正交单位化后拼接为V(假设A=U\Sigma V^T成立,则A^TA=V\Sigma^T U^TU\Sigma V^T=V\left (\Sigma^T\Sigma \right ) V^T,下面的证明同理),再对AA^T进行特征值分解得到特征向量经正交单位化后拼接为U,而A^TAAA^Tn个特征值的根号即为奇异值,也就是\Sigma主对角线元素;
  • 性质
    • 对于U的第i\vec{u}_ii\leqslant n)及V对应列\vec{v}_i,有A\vec{v}_i=\sigma_i\vec{u}_i,或\vec{u}_i^TA\vec{v}_i=\sigma_i,其中\sigma_i\Sigma对应位置奇异值,而当i> ni\neq j时,\vec{u}_i^TA\vec{v}_j=0
    • 奇异值通常从大到小排列,但递减得特别快,很多情况下前10%甚至1%的奇异值之和就占了全部奇异值之和的99%以上,因此有时可以用一部分来描述奇异值分解:A_{m\times n}\approx U_{m\times k}\Sigma_{k\times k} V_{k\times n}^T,其中kn小得多,这个性质也是PCA的重要基础;

1.2.4. 矩阵求导

  • \frac{d\left | A \right |}{dA} = \left | A \right|(A^{-1})^T
  • \frac{d(A^{-1})}{dA} = -X^{-1}dXX^{-1}
  •  \frac{\partial tr(AB) }{\partial A}=B^{T}
  •  \frac{\partial \left | A \right |}{\partial A}=\left | A \right |\cdot A^{-1}
  • \frac{\partial ln\left | A \right |}{\partial A}=A^{-1}

因此\frac{\partial (\vec{x}^TA\vec{x})}{\partial \vec{x}}=2A\vec{x}

  • \frac{\partial \vec{x} \vec{x}^T}{\partial \vec{x}}= \vec{x}\bigotimes I_p+I_p\bigotimes \vec{x}

1.2.5. Jensen不等式

  • 对任意凸函数,其上两点(a,f(a))(b,f(b)),有f(\frac{a+b}{2})\geqslant \frac{f(a)+f(b)}{2}
  • 推广对任意t\in [0,1],有f(a+t(b-a))\geqslant f(a)+t(f(b)-f(a))
  • 推广对任意随机变量 x ,有f(EX) \geqslant E[f(x)]

2. 高斯分布

2.1. 高维高斯分布的理解

2.1.1. 马氏距离

  • 对于两个向量\vec{x_1}\vec{x_2},其马氏距离为(\vec{x_1}-\vec{x_2})^T\Sigma^{-1}(\vec{x_1}-\vec{x_2})\Sigma为每个特征对应距离的权重,是一个正定对称矩阵,当\Sigma=I时,马氏距离退化为欧氏距离(p维空间两个坐标点的空间距离)(\vec{x_1}-\vec{x_2})^T(\vec{x_1}-\vec{x_2})=\sqrt{\sum _{i=1}^{p}(x_{1i}-x_{2i})^2}

2.1.2. 推导

  • 高斯分布的概率密度函数为:

p(\vec{x})=\frac{1}{2\pi^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}\exp(-\frac{1}{2}(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu})) 

        令\Delta =(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu}),因为\Sigma是一个正定对称矩阵,因此可以对其进行矩阵分解,存在正交矩阵U和对角阵\Lambda ,使得

\begin{matrix} \begin{aligned}\Sigma &=U\Lambda U^T\\ &=\left ( \vec{u}_1, \cdots , \vec{u}_p \right )\begin{pmatrix} \lambda _1 & \cdots &0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & \lambda _p \end{pmatrix} \begin{pmatrix} \vec{u}_1^T\\ \vdots \\ \vec{u}_p^T \end{pmatrix} \\ &=\sum _{i=1}^{p}\vec{u}_i\lambda _i\vec{u}_i^T \end{aligned} & &\begin{aligned}\Sigma^{-1} &=\sum _{i=1}^{p}\vec{u}_i\lambda^{-1} _i\vec{u}_i^T \end{aligned} \end{matrix}

        因此可得

\begin{aligned}\Delta &=(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu})\\ & =(\vec{x}-\vec{\mu})^T\left [\sum _{i=1}^{p}\vec{u}_i\lambda^{-1} _i\vec{u}_i^T \right ](\vec{x}-\vec{\mu})\\ & =\sum _{i=1}^{p}(\vec{x}-\vec{\mu})^T\vec{u}_i\lambda^{-1} _i\vec{u}_i^T(\vec{x}-\vec{\mu})\\ &\overset{y_i=(\vec{x}-\vec{\mu})^T\vec{u}_i}{=}\sum _{i=1}^{p}y_i\lambda^{-1} _iy_i= \sum _{i=1}^{p}\frac{y_i^2}{\lambda}\end{aligned}

        最终得到的结果在等于定值的情况下(例如\frac{y_1^2}{\lambda_1}+\frac{y_2^2}{\lambda_2}=1)为一个高维空间的椭球体,其中经过的投影y_i=(\vec{x}-\vec{\mu})^T\vec{u}_i即为将该球体球心移至坐标原心,再对每个 \vec{u}_i 坐标轴分别进行投影的过程(下图中长短轴分别为 \lambda _1^{\frac{1}{2}} 和 \lambda _2^{\frac{1}{2}} );

2.2. 极大似然估计

2.2.1. 问题定义:

  • 对于一组数据Data:X_{N\times p}=(\vec{x}_1, \vec{x}_2, \cdots ,\vec{x}_N)^T,\vec{x}_i \overset{iid}{\sim}N(\vec{\mu} ,\Sigma ),其未知参数为\Theta =(\vec{\mu },\Sigma ),对其进行极大似然估计;

2.2.2. 过程推导

  • 似然函数

 \begin{aligned}L(\Theta )&=\log[\prod_{i=1}^N p(\vec{x}_i|\Theta) ]\\ &=\sum _{i=1}^N \log [p(\vec{x}_i|\Theta)] \end{aligned}

        将概率密度函数带入得:

\begin{aligned}L(\Theta ) =\sum_{i=1}^N \left [ -\frac{p}{2}\log(2\pi)-\frac{1}{2}\log(\left | \Sigma \right |)-\frac{1}{2}(\vec{x}_i-\vec{\mu})^T\Sigma ^{-1}(\vec{x}_i-\vec{\mu}) \right ]\end{aligned}

        对\vec{\mu}求导,并令其为0,则有:

\begin{aligned} \frac{dL(\Theta )}{d\vec{\mu}} &=\sum _{i=1}^N (\vec{x}_i-\vec{\mu})^T(\Sigma ^{-1}+{\Sigma^{-1}}^T)\\ &=\sum _{i=1}^N 2(\vec{x}_i-\vec{\mu})^T\Sigma ^{-1}=0\end{aligned}\Rightarrow\hat{\vec{\mu}} = \frac{\sum _{i=1}^N\vec{x}_i}{N}

        再对\Sigma求导:

\begin{aligned} {dL(\Theta )} &=tr(dL(\Theta))\\ &=-\frac{N}{2}\Sigma ^{-1}d\Sigma- \frac{1}{2}\sum _{i=1}^Ntr((\vec{x}_i-\vec{\mu})^T(-\Sigma^{-1}d\Sigma\Sigma^{-1})(\vec{x}_i-\vec{\mu}))\\ &=-\frac{N}{2}\Sigma ^{-1}d\Sigma+ \frac{1}{2}\sum _{i=1}^Ntr((\vec{x}_i-\vec{\mu})^T\Sigma^{-1})(\Sigma^{-1}(\vec{x}_i-\vec{\mu}))d\Sigma\\ &\overset{tr(a^Tb)=tr(ba^T)}{=}-\frac{N}{2}\Sigma ^{-1}d\Sigma+ \frac{1}{2}\sum _{i=1}^Ntr((\Sigma^{-1}(\vec{x}_i-\vec{\mu}))(\vec{x}_i-\vec{\mu})^T\Sigma^{-1})d\Sigma\\ \end{aligned}

        又根据df=tr((\frac{df}{dx})^Tdx),可得:

\frac{dL(\Theta)}{d\Sigma} =-\frac{N}{2}\Sigma ^{-1}+\frac{1}{2}\sum_{i=1}^N\Sigma^{-1}(\vec{x}_i-\vec{\mu})(\vec{x}_i-\vec{\mu})^T\Sigma^{-1}

        令其为0,可得:

\hat{\Sigma}=\frac{1}{N}\sum_{i=1}^N(\vec{x}_i-\vec{\mu})(\vec{x}_i-\vec{\mu})^T

  • 对高斯分布进行极大似然估计是一种有偏估计,均值无偏而方差有偏,估计得到的\hat{\Sigma}=\frac{1}{N}\sum_{i=1}^N(\vec{x}_i-\vec{\mu})(\vec{x}_i-\vec{\mu})^TE(\hat{\Sigma})=\frac{N-1}{N}\Sigma,而真正的无偏估计为\hat{\Sigma}=\frac{1}{N-1}\sum_{i=1}^N(\vec{x}_i-\vec{\mu})(\vec{x}_i-\vec{\mu})^T

2.3. 高斯分布的局限性

  • 协方差矩阵 \Sigma 尺寸为p\times p,虽然其为实对称,但参数量依然有\frac{p^2+p}{2}个,对于机器学习类问题,特征维度往往很大,需要学习的参数量也就很大,常见的解决方法有:
    • 将 \Sigma 矩阵近似为对角阵;
    • 利用 x 的线性变换将 \Sigma 矩阵变为单位阵(或单位阵的k倍);
  • 复杂性不够,很多数据的分布难以用一个分布拟合;

2.4. 高维高斯分布的边缘概率分布

2.4.1. 问题定义

  • 对于一组数据Data:X_{N\times p}=(\vec{x}_1, \vec{x}_2, \cdots ,\vec{x}_N)^T,\vec{x}_i \overset{iid}{\sim}N(\vec{\mu} ,\Sigma ),可将其分为两组\vec{x}_i=\begin{pmatrix} \vec{x}_{ia}\\ \vec{x}_{ib} \end{pmatrix},其中 \vec{x}_{ia} 为 m 维向量, \vec{x}_{ib} 为 n 维向量,m+n=p,因此则可以将其他参数划分为\vec{\mu}=\begin{pmatrix} \vec{\mu}_{a}\\ \vec{\mu}_{b} \end{pmatrix}\Sigma =\begin{pmatrix} \Sigma_{aa} & \Sigma_{ab}\\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix}
  • 此时概率密度函数 p(\vec{x}) 可以认为是联合概率分布 p(\vec{x}_a,\vec{x}_b)
  • 需要求的是边缘概率密度 p(\vec{x}_a),p(\vec{x}_b)

2.4.2. 定理基础

  • 高斯分布的线性定理:若 \vec{x}_i \overset{iid}{\sim}N(\vec{\mu} ,\Sigma ) ,则对于 \vec{y}=A\vec{x}+\vec{b} ,有 \vec{y}_i \overset{iid}{\sim}N(A\vec{\mu}+\vec{b} ,A\Sigma A^T )

2.4.3. 过程推导

  • 首先对于 \vec{x}_{ia},\vec{x}_{ib} ,可以写成矩阵形式:\vec{x}_{ia}=(I_m,0)\cdot\vec{x}_i,\vec{x}_{ib}=(0,I_n)\cdot\vec{x}_i
  • 则可以令A_1 = (I_m,0),A_2=(0,I_n),可得 E(\vec{x}_a)=A_1\vec{\mu}=\vec{\mu}_a,D(\vec{x}_a)=A_1\Sigma A_1^T=\Sigma_{aa} ;
  • 因此 \vec{x}_{ia} \overset{iid}{\sim}N(\vec{\mu}_a ,\Sigma_{aa} ) ,同理可得 \vec{x}_{ib} \overset{iid}{\sim}N(\vec{\mu}_b ,\Sigma_{bb} ); 

2.5. 高维高斯分布的条件概率分布

2.5.1. 问题定义

  • 同2.4.,\vec{x}_i=\begin{pmatrix} \vec{x}_{ia}\\ \vec{x}_{ib} \end{pmatrix}\vec{\mu}=\begin{pmatrix} \vec{\mu}_{a}\\ \vec{\mu}_{b} \end{pmatrix}\Sigma =\begin{pmatrix} \Sigma_{aa} & \Sigma_{ab}\\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix}
  • 需要求的是条件概率密度 p(\vec{x}_a|\vec{x}_b),p(\vec{x}_b|\vec{x}_a) ;

2.5.2. 过程推导

  • 根据schur complementary矩阵*,可设定一个中间量 \vec{x}_{b\cdot a}=\vec{x}_b-\Sigma_{ba}\Sigma_{aa}^{-1}\vec{x}_a ;
  • 对其求期望和方差可得:E(\vec{x}_{b\cdot a})=E(\vec{x}_b-\Sigma_{ba}\Sigma_{aa}^{-1}\vec{x}_a)=\vec{\mu}_b-\Sigma_{ba}\Sigma_{aa}^{-1}\vec{\mu}_aD(\vec{x}_{b\cdot a})=D(\vec{x}_b-\Sigma_{ba}\Sigma_{aa}^{-1}\vec{x}_a)=\Sigma_{bb}-\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab}
  • 记 \vec{\mu}_{b\cdot a}=\vec{\mu}_b-\Sigma_{ba}\Sigma_{aa}^{-1}\vec{\mu}_a ,\Sigma_{bb\cdot a}=\Sigma_{bb}-\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab},则\vec{x}_{b\cdot a} \overset{iid}{\sim}N(\vec{\mu}_{b\cdot a} ,\Sigma_{bb\cdot a} )
  • 因此\vec{x}_b=\vec{x}_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}\vec{x}_a,则E(\vec{x}_b|\vec{x}_a)=\vec{\mu}_{b\cdot a}+\Sigma_{ba}\Sigma_{aa}^{-1}\vec{x}_a=\vec{\mu}_{b}+\Sigma_{ba}\Sigma_{aa}^{-1}(\vec{x}_a-\vec{\mu}_a)D(\vec{x}_b|\vec{x}_a)=\Sigma_{bb\cdot a},则\vec{x}_{b}|\vec{x}_{a} \overset{iid}{\sim}N(\vec{\mu}_{b}+\Sigma_{ba}\Sigma_{aa}^{-1}(\vec{x}_a-\vec{\mu}_a) ,\Sigma_{bb}-\Sigma_{ba}\Sigma_{aa}^{-1}\Sigma_{ab} ),同理可得\vec{x}_{a}|\vec{x}_{b} \overset{iid}{\sim}N(\vec{\mu}_{a}+\Sigma_{ab}\Sigma_{bb}^{-1}(\vec{x}_b-\vec{\mu}_b) ,\Sigma_{aa}-\Sigma_{ab}\Sigma_{bb}^{-1}\Sigma_{ba})

2.6. 线性高斯模型对称求解

2.6.1. 问题定义

  • 已知 \vec{x}\sim N(\vec{x}|\vec{\mu},\Lambda ^{-1}),\vec{y}|\vec{x}\sim N(\vec{y}|A\vec{x}+\vec{b},L ^{-1}) ,求 p(\vec{y}),p(\vec{x}|\vec{y}) ;
  • \vec{y}=A\vec{x}+\vec{b}+\vec{\varepsilon },\vec{\varepsilon }\sim N(0, L^{-1})\vec{\varepsilon }\vec{x}相互独立;

2.6.2. 过程推导

  • 首先求 p(\vec{y}) ,E(\vec{y})=A\vec{\mu}+\vec{b}D(\vec{y})=A\Lambda ^{-1} A^T+L^{-1},可得\vec{y}\sim N(A\vec{\mu}+\vec{b} ,A\Lambda ^{-1} A^T+L^{-1} )
  • 再求解p(\vec{x}|\vec{y}),直接求解比较复杂,因此我们先求\vec{y}\vec{x}的联合概率分布,再通过2.5的结论套用即可得到结果;
  • 易知\vec{z}=\begin{pmatrix} \vec{x}\\ \vec{y} \end{pmatrix}\sim N(\begin{bmatrix} \vec{\mu}\\ A\vec{\mu}+\vec{b} \end{bmatrix},\begin{bmatrix} \Lambda ^{-1} & \Delta \\ \Delta^T &A\Lambda ^{-1} A^T+L^{-1} \end{bmatrix}),其中\Delta是需要求解的;

\begin{aligned}\Delta &=Cov(\vec{x},\vec{y})\\ &=E[(\vec{x}-E(\vec{x}))(\vec{y}-E(\vec{y}))^T]\\ &=E[(\vec{x}-\vec{\mu})(A\vec{x}+\vec{b+\vec{\varepsilon }}-A\vec{\mu}-\vec{b})^T]\\ &=E[(\vec{x}-\vec{\mu})(\vec{x}-\vec{\mu})^TA^T]+E[(\vec{x}-\vec{\mu})]E[\vec{\varepsilon}]\\ &=D(\vec{x})A^T=\Lambda ^{-1}A^T \end{aligned}

  • \vec{x}=\vec{z},\vec{x}_a=\vec{x},\vec{x}_b=\vec{y},即可套用2.5的计算p(\vec{x}_a|\vec{x}_b)的结论计算p(\vec{x}|\vec{y})

3. 线性回归

3.1. 最小二乘法

3.1.1. 问题定义

  • Data={(\vec{x}_1,y_1),(\vec{x}_2,y_2),\cdots,(\vec{x}_n,y_n)}X=(\begin{matrix} \vec{x}_1 & \vec{x}_2 & \cdots & \vec{x}_n \end{matrix})^T{}_{N\times p}Y=(\begin{matrix} y_1 & y_2 & \cdots & y_n \end{matrix}){}_{N\times1}
  • 对于增广矩阵\hat{X}=(\begin{matrix} \hat{\vec{x}}_1 & \hat{\vec{x}}_2 & \cdots &\hat{ \vec{x}}_n \end{matrix})^T=(\begin{matrix} \begin{bmatrix} 1\\ \vec{x}_1 \end{bmatrix} & \begin{bmatrix} 1\\ \vec{x}_2 \end{bmatrix} & \cdots & \begin{bmatrix} 1\\ \vec{x}_n \end{bmatrix} \end{matrix})^T{}_{N\times (p+1)},需要找到一个参数矩阵W{}_{(p+1)\times1},使得W=\underset{W}{argmin}L(W)=\left \|\hat{X}W-Y \right \|=\sum _{i=1}^N(W^T\hat{\vec{x}}_i-y_i )^2

3.1.2. 求出解析解

  • L(W) = (W^T\hat{X}^T-Y^T)(\hat{X}W-Y)=W^T\hat{X}^T\hat{X}W-2W^T\hat{X}^TY+Y^TY
  • 将这个损失函数对参数矩阵求导可得:\frac{\partial L(W)}{\partial W}=2\hat{X}^T\hat{X}W-2\hat{X}^TY
  • 令其为0,可得\hat{W}=(\hat{X}^T\hat{X})^{-1}X^TY,其中X^\dagger =(\hat{X}^T\hat{X})^{-1}X^T被称为\hat{X}的伪逆;

3.2. 概率角度

3.2.1. 问题定义

  • 同样是Data={(\vec{x}_1,y_1),(\vec{x}_2,y_2),\cdots,(\vec{x}_n,y_n)}X=(\begin{matrix} \vec{x}_1 & \vec{x}_2 & \cdots & \vec{x}_n \end{matrix})^T{}_{N\times p},增广矩阵\hat{X}=(\begin{matrix} \hat{\vec{x}}_1 & \hat{\vec{x}}_2 & \cdots &\hat{ \vec{x}}_n \end{matrix})^T=(\begin{matrix} \begin{bmatrix} 1\\ \vec{x}_1 \end{bmatrix} & \begin{bmatrix} 1\\ \vec{x}_2 \end{bmatrix} & \cdots & \begin{bmatrix} 1\\ \vec{x}_n \end{bmatrix} \end{matrix})^T{}_{N\times (p+1)},参数矩阵W{}_{(p+1)\times1}Y=(\begin{matrix} y_1 & y_2 & \cdots & y_n \end{matrix}){}_{N\times1}
  • 可以知道,y_i 和 \vec{x}_i 有着一定的线性关系,但真实数据存在随机误差,因此可以定义随机误差\varepsilon \sim N(0, \sigma^2),使得y=\hat{\vec{x}}^TW+\varepsilon,则可得y|\hat{\vec{x}}W\sim N(\hat{\vec{x}}W,\sigma^2)

3.2.2. 求出解析解

  • 使用最大似然估计法求解参数,由y|\hat{\vec{x}}W\sim N(\hat{\vec{x}}W,\sigma^2),可得似然函数L(W)=\sum_{i=1}^N\log[p(y_i|\hat{\vec{x}};W)]=N\log(\frac{1}{\sqrt{2\pi}\sigma})-\frac{1}{2\sigma^2}(y_i-\hat{\vec{x}}_i^TW)^2,因此优化目标变成了W=\underset{W}{argmin}L(W)=\left \|\hat{X}W-Y \right \|=\sum _{i=1}^N(W^T\hat{\vec{x}}_i-y_i )^2,与最小二乘法相同;

3.3. 正则化

3.3.1. 为什么要正则化

  • 从前面的推论可以知道参数矩阵的解析解\hat{W}=(\hat{X}^T\hat{X})^{-1}X^TY,一般来说N\gg p,但特殊情况样本量较小时,可能导致\hat{X}^T\hat{X}不可逆,无法求出解析解;
  • 从拟合情况看,如果不能满足N\gg p,维数较高而样本不足,则容易导致过拟合问题,对于过拟合情况可以用三种方法:
    • 采集更多数据;
    • 特征选择/特征提取(PCA)
    • 正则化

3.3.2. 范数

  • 通俗来讲,一个向量的L_p范数可以定义为\left \| \vec{x} \right \|_p=\sqrt[p]{\sum_{i=1}^m\left | x_i \right |^p}
    • L1范数:Lasso,即为向量各元素的绝对值之和;
    • L2范数:Ridge,岭回归,\left \| \vec{x} \right \|_2=\sqrt[2]{\sum_{i=1}^m\left | x_i \right |^2}=\vec{x}^T\vec{x}

3.3.3. L2正则化

  • 正则化的框架:J(W)=L(W)+\lambda P(W),即在原有损失函数上加上一个正则化惩罚;
  • 对于L2正则化,则对应的优化函数为J(W)=L(W)+\lambda W^TW,这个优化函数也被称为岭回归;
  • 对其进行化简和求导(同3.1.2.)可以得到新的解析解\hat{W}=(\hat{X}^T\hat{X}+\lambda I)^{-1}X^TY,解决了\hat{X}^T\hat{X}不可逆的问题

3.3.4. 从贝叶斯角度理解线性回归和L2正则化

  • 从贝叶斯派的角度求解线性回归问题,则是已知XY的分布p(\vec{x})p(y),对于参数矩阵假设一个先验分布p(W),采用最大后验估计法求解W=\underset{W}{argmax}p(W|y)
  • 假设参数矩阵服从正态分布,即W\sim N(\vec{0}, \Sigma_0),则根据贝叶斯公式:p(W|y)=\frac{p(y|W)\cdot p(W)}{p(y)},因为p(y)已知,所以\max(p(W|y))=\max\left \{\frac{p(y|W)\cdot p(W)}{p(y)} \right \}=\max\left \{p(y|W)\cdot p(W)\right \}
  • 根据3.2.1.可知y|\hat{\vec{x}}W\sim N(\hat{\vec{x}}W,\sigma^2I),因为p(\vec{x})已知,所以也可以视为y|W\sim N(\hat{\vec{x}}W,\sigma^2I),因此p(y|W)\cdot p(W)=\frac{1}{\sqrt{2\pi}\sigma^2}\frac{1}{\sqrt{2\pi}\sigma_0^2}\cdot \exp{\left\{-\frac{(y-\hat{\vec{x}}^TW)^T(y-\hat{\vec{x}}^TW)}{2\sigma^2}-\frac{W^T\Sigma_0^{-1}W}{2}\right\}}
  • \Sigma_0=\sigma_0^2I(即方差各向同性)时,\frac{W^T\Sigma_0^{-1}W}{2}=\frac{W^TW}{2\sigma_0^2}=\frac{\left \| W \right \|^2}{2\sigma_0^2},因此有

\begin{aligned}W&=\underset{W}{argmax}\left \{p(W|y) \right \}\\ &=\underset{W}{argmax}\left \{ p(y|W)\cdot p(W)\right \} \\ &=\underset{W}{argmax}\left \{\log( p(y|W)\cdot p(W))\right \} \\ &=\underset{W}{argmin}\left \{\frac{(y-\hat{\vec{x}}^TW)^T(y-\hat{\vec{x}}^TW)}{2\sigma^2}+\frac{\left \| W \right \|^2}{2\sigma_0^2} \right \} \\&=\underset{W}{argmin}\left \{(y-\hat{\vec{x}}^TW)^T(y-\hat{\vec{x}}^TW)+\frac{\sigma^2}{\sigma_0^2} \left \| W \right \|^2\right \} \\\end{aligned}

        此时的优化函数和L2正则化是相同的,只是将 \lambda 替换为了\frac{\sigma^2}{\sigma_0^2}

  • 使用全量数据求解参数矩阵W(求解线性回归问题):
    • \begin{aligned} P(W|X,Y)&\propto \exp\left \{ -\frac{1}{2\sigma^2}(Y-XW)^T(Y-XW)-\frac{1}{2}W^T\Sigma_0^{-1}W \right \}\\ &= \exp\left \{ -\frac{1}{2\sigma^2}(Y^TY-2Y^TXW+W^TX^TXW)-\frac{1}{2}W^T\Sigma_0^{-1}W \right \}\\ &\propto \exp\left \{ \underset{\text{Quadratic term}}{\underbrace{-\frac{1}{2}\left [ W^T(\sigma^{-2}X^TX+\Sigma_0^{-1})W \right ]}}-\underset{\text{Linear term}}{\underbrace{\sigma^{-2}Y^TXW}} \right \}\\ \end{aligned},其中Quadratic term表示变量的二次项,Linear term表示变量的一次项;
    • 根据高斯分布标准形式:\begin{aligned} P(\vec{x})&=\frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}\exp\{-\frac{1}{2}\left[(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu}) \right ]\}\\ &=\frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}\exp\{-\frac{1}{2}\left[ \vec{x}^T\Sigma^{-1}\vec{x}-2\vec{\mu}^T\Sigma^{-1}\vec{x}+\vec{\mu}^T\Sigma^{-1}\vec{\mu} \right ]\}\\ \end{aligned},根据二次项对应系数矩阵可以求得方差矩阵\Sigma,进而根据一次项对应系数矩阵可以求得均值向量\vec{\mu}
    • 因此将上面两个式子对应联立求解,可得:
      • \Sigma_w=(\sigma^{-2}X^TX+\Sigma_0^{-1})^{-1}
      • \mu_w=\sigma^{-4}(X^TX+\sigma^{2}\Sigma_0^{-1})^{-1}X^TY
      • P(W|X,Y)\sim N(\mu_w,\Sigma_w),当\Sigma_0=\sigma_0^2I时,取均值结果与3.3.3. L2正则化的解析解形式相同,只是将 \lambda 替换为了\frac{\sigma^2}{\sigma_0^2}

3.4. 线性回归的特性

4. 线性分类

4.1. 线性分类的典型模型

  • 硬分类:输出结果只有0或1这种离散结果;
    • 线性判别分析 Fisher
    • 感知机
  • 软分类:会输出0-1之间的值作为各个类别的概率;
    • 概率生成模型:高斯判别分析GDA、朴素贝叶斯,主要建模的是p(x,y)
    • 概率判别模型:逻辑回归,主要建模的是p(y|x)

4.2. 感知机

4.2.1. 基本模型

  • 模型:f(\vec{x})=sign(\hat{\vec{x}}^TW)
  • 思想:错误驱动,先随机初始化一个W,研究错误分类的样本来调整W
  • 策略:使用损失函数L(W)=\sum_{\hat{\vec{x}}\in D}-y_i\cdot \hat{\vec{x}}^TW,其中D表示错分类的样本;

4.3. 线性判别分析 Fisher

4.3.1. 问题定义

  • 对于一个二分类问题,将样本分为X_{c_1}=\left \{ \hat{\vec{x}}_i|y_i=+1 \right \}X_{c_1}=\left \{ \hat{\vec{x}}_i|y_i=-1 \right \},这两组的样本数分别为N_1N_2N_1+N_2=N
  • 寻找一个投影超平面W,使所有样本点在这个平面的投影可以做到类内间距小,类间间距大

4.3.2. 过程推导

  • 样本\hat{\vec{x}}在超平面W上的投影可以表示为z=\hat{\vec{x}}^T\cdot W\bar{z}=\frac{1}{N}\sum_{i=1}^Nz_i=\frac{1}{N}\sum_{i=1}^N\hat{\vec{x}}_i^T\cdot W
  • 对投影求方差S_z=\frac{1}{N}\sum_{i=1}^N(z_i-\bar{z})^2=\frac{1}{N}\sum_{i=1}^N(\hat{\vec{x}}_i^T\cdot W-\bar{z})^2
    • \bar{z}_1=\frac{1}{N_1}\sum_{i=1}^{N_1}\hat{\vec{x}}_i^T\cdot WS_{z_1}=\frac{1}{N_1}\sum_{i=1}^{N_1}(\hat{\vec{x}}_i^T\cdot W-\bar{z}_1)^2
    • \bar{z}_2=\frac{1}{N_2}\sum_{i=1}^{N_2}\hat{\vec{x}}_i^T\cdot WS_{z_2}=\frac{1}{N_2}\sum_{i=1}^{N_2}(\hat{\vec{x}}_i^T\cdot W-\bar{z}_2)^2
  • 类内间距分别为S_{z_1}S_{z_2},类间间距为(\bar{z}_1-\bar{z}_2)^2
  • 为了尽可能类内间距小,类间间距大,将目标函数定义为J(W)=\frac{(\bar{z}_1-\bar{z}_2)^2}{S_{z_1}+S_{z_2}},则W=\underset{W}{argmax}J(W)
  • 对目标函数进行化简:
    • \begin{aligned} (\bar{z}_1-\bar{z}_2)^2&=(\frac{1}{N_1}\sum_{i=1}^{N_1}\hat{\vec{x}}_i^T\cdot W-\frac{1}{N_2}\sum_{i=1}^{N_2}\hat{\vec{x}}_i^T\cdot W)^2\\ &=[(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})^TW]^2\\ &=W^T(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})^TW\end{aligned}
    • \begin{aligned}S_{z_1}&=\frac{1}{N_1}\sum_{i=1}^{N_1}(\hat{\vec{x}}_i^T\cdot W-\bar{z}_1)^2\\&=\frac{1}{N_1}\sum_{i=1}^{N_1}W^T(\vec{x}_i-\bar{\vec{x}}_{C_1})(\vec{x}_i-\bar{\vec{x}}_{C_1})^TW\\ & =W^T\left [\frac{1}{N_1}\sum_{i=1}^{N_1}(\vec{x}_i-\bar{\vec{x}}_{C_1})(\vec{x}_i-\bar{\vec{x}}_{C_1})^T \right ]W\\ &=W^TS_{C_1}W\end{aligned}
    • 同理可得S_{z_2}=W^TS_{C_2}W
    • 所以目标函数化为J(W)=\frac{(\bar{z}_1-\bar{z}_2)^2}{S_{z_1}+S_{z_2}}=\frac{W^T(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})^TW}{W^T(S_{C_1}+S_{C_2})W}
  • 定义类内方差S_w=S_{C_1}+S_{C_2},类间方差S_b=(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})^T,因此目标函数被表示为J(W)=\frac{W^TS_bW}{W^TS_wW},对其求导得\frac{\partial J(W)}{\partial W}=2S_bW(W^TS_wW)^{-1}+W^TS_bW\cdot(-1)\cdot (W^TS_wW)^{-2}\cdot2S_wW,令其为0可得

\begin{aligned}W&=\frac{W^TS_wW}{W^TS_bW}S_w^{-1}S_bW\\&=\frac{W^TS_wW}{W^TS_bW}S_w^{-1}(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})^TW\\&=\frac{W^TS_wW\cdot (\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})^TW}{W^TS_bW}S_w^{-1}(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2}) \end{aligned}

        因为W是一个单位向量,所以我们只关心其方向而不关心其长度,所以最终得到W\propto S_w^{-1}(\bar{\vec{x}}_{C_1}-\bar{\vec{x}}_{C_2})

4.4. 逻辑回归

4.4.1. 基本思想

  • 在线性回归中引入非线性激活函数,使其可以将回归结果映射为概率值,作为类别的概率;
  • 将这个模型看做一个条件概率分布的建模,输入 \hat{\vec{x}} ,通过model=p(y|\hat{\vec{x}}),输出 y 的离散取值;

4.4.2. Sigmoid 激活函数

  • 基本公式:\sigma (z)=\frac{1}{1+e^{-z}}
  • 特殊取值:
    • z\rightarrow -\infty时,\lim \sigma (z)=0
    • z=0时,\sigma (z)=\frac{1}{2}
    • z\rightarrow \infty时,\lim \sigma (z)=1
  • 函数图像:

4.4.3. 模型推导

  • 根据条件概率建模的思想:

\begin{aligned}p_1&=p(y=1|\hat{\vec{x}})=\sigma(\hat{\vec{x}}^TW)= \frac{1}{1+e^{-\hat{\vec{x}}^TW}}\\ p_0&=p(y=0|\hat{\vec{x}})=1-\sigma(\hat{\vec{x}}^TW)= \frac{e^{-\hat{\vec{x}}^TW}}{1+e^{-\hat{\vec{x}}^TW}}\\ \end{aligned}

  • 因此将整个模型写作p(y|x)=p_1^y\cdot p_0^{1-y},当y=0时,p=p_0,当y=1时,p=p_1
  • 用极大似然估计法求解模型:

\begin{aligned} \hat{W}&=\underset{W}{argmax }\log p(Y|X)\\ &=\underset{W}{argmax }\sum_{i=1}^N\left [ y_i\cdot \log\sigma(\hat{\vec{x}}^TW)+(1-y_i)\cdot \log(1-\sigma(\hat{\vec{x}}^TW)) \right ]\\\end{aligned}

  • 求梯度得

\begin{aligned} \bigtriangledown grad_W &=\sum_{i=1}^N \left [y_i\cdot (1-\sigma(\hat{\vec{x}}^TW))\cdot\hat{\vec{x}} - (1-y_i)\cdot \sigma(\hat{\vec{x}}^TW)\cdot\hat{\vec{x}} \right ]\\ &=\sum_{i=1}^N \left [y_i-\sigma(\hat{\vec{x}}^TW) \right ]\cdot\hat{\vec{x}}\\ \end{aligned}

4.5. 高斯判别分析

4.5.1. 概率判别式模型与概率生成式模型的区别

  • 概率判别式模型主要计算条件概率密度p(y|\vec{x}),取令该概率最大的y为分类结果;
  • 概率生成式模型并不需要计算具体的p(y|\vec{x})值,而是直接思考p(y=1|\vec{x})p(y=0|\vec{x})的结果谁更大,根据贝叶斯公式p(y|\vec{x}) = \frac{p(\vec{x}|y)\cdot p(y)}{p(\vec{x})},将目标函数变为\hat{y}=\underset{y}{argmax} p(y|\vec{x})=\underset{y}{argmax} p(\vec{x}|y)\cdot p(y),其中p(y=1)=\phip(y=0)=1-\phi,可以合并为p(y)=\phi ^y\cdot (1-\phi)^{1-y}

4.5.2. 高斯概率假设

  • 在高斯判别模型中,假设条件分布\vec{x}|y是遵从高斯概率分布的,即:\vec{x}|y=0\sim N(\mu_1, \Sigma),\vec{x}|y=1\sim N(\mu_2, \Sigma)
  • 使用对数似然求解目标函数,可得

\begin{aligned} L(\theta)&=\log\prod_{i=1}^N p(\vec{x}_i,y_i)\\ &=\sum _{i=1}^N\left [\log p(\vec{x}_i|y_i)+\log p(y_i) \right ] \\ &=\sum _{i=1}^N\left \{ y_i\cdot \left [-\frac{1}{2}\left (\vec{x}_i-\vec{\mu}_1 \right )^T\Sigma^{-1} \left(\vec{x}_i-\vec{\mu}_1 \right )-\frac{p}{2}\log2\pi--\frac{1}{2}\log|\Sigma| \right ]+\left (1-y_i \right )\cdot \left [-\frac{1}{2}\left (\vec{x}_i-\vec{\mu}_2 \right )^T\Sigma^{-1} \left(\vec{x}_i-\vec{\mu}_2 \right )-\frac{p}{2}\log2\pi--\frac{1}{2}\log|\Sigma| \right ] +y_i \cdot\log\phi + (1-y_i)\cdot\log(1-\phi)\right \} \end{aligned}

4.5.3. 求解模型

  • 先求解参数\phi,先求偏导\frac{\partial L(\theta)}{\partial \phi}=\sum_{i=1}^N\left (\frac{y_i}{\phi}-\frac{1-y_i}{1-\phi} \right ),令其为0,求得\hat{\phi}=\frac{1}{N}\sum_{i=1}^{N}y_i,也就是样本中各个标签出现的频率即为最优概率值;
  • 再求解参数\vec{\mu}_1,由于其他和\vec{\mu}_1无关的部分求偏导后都得0,所以从目标函数中单独取出和\vec{\mu}_1相关的部分,即\hat{\vec{\mu}}_1=\underset{\vec{\mu}_1}{argmax}\sum_{i=1}^Ny_i\cdot \left [-\frac{1}{2}\left (\vec{x}_i-\vec{\mu}_1 \right )^T\Sigma^{-1} \left(\vec{x}_i-\vec{\mu}_1 \right ) \right ],即\hat{\vec{\mu}}_1=\underset{\vec{\mu}_1}{argmax}\sum_{\vec{x_i}\in C_1}\left [-\frac{1}{2}\left (\vec{x}_i-\vec{\mu}_1 \right )^T\Sigma^{-1} \left(\vec{x}_i-\vec{\mu}_1 \right ) \right ],其中C_1表示\left \{ \vec{x}_i |y_i=1\right \},因此求解方法等同于高维高斯分布的极大似然估计(2.2.2.),结果为\hat{\vec{\mu}}_1=\frac{\sum_{i=1}^Ny_i\cdot \vec{x}_i}{N_1}=\frac{\sum_{\vec{x}_i\in C_1}\vec{x}_i}{N_1},同理\hat{\vec{\mu}}_2=\frac{\sum_{i=1}^N\left (1-y_i \right )\cdot \vec{x}_i}{N_2}=\frac{\sum_{\vec{x}_i\in C_2}\vec{x}_i}{N_2}
  • 最后求\Sigma,同样取出目标函数中和\Sigma相关的部分,得:\begin{aligned} \hat{\Sigma}=\underset{\Sigma}{argmax}\sum_{\vec{x}_i\in C_1}\left [-\frac{1}{2}\left (\vec{x}_i-\vec{\mu}_1 \right )^T\Sigma^{-1} \left(\vec{x}_i-\vec{\mu}_1 \right ) \right ]+\sum_{\vec{x}_i\in C_2}\left [-\frac{1}{2}\left (\vec{x}_i-\vec{\mu}_2 \right )^T\Sigma^{-1} \left(\vec{x}_i-\vec{\mu}_2 \right ) \right ] - \frac{N}{2}\log |\Sigma| \end{aligned},对其求偏导得\begin{aligned} \frac{dL(\theta)}{d\Sigma} =-\frac{N}{2}\Sigma ^{-1}+\frac{1}{2}\sum_{\vec{x}_i\in C_1}\Sigma^{-1}(\vec{x}_i-\vec{\mu}_1)(\vec{x}_i-\vec{\mu}_1)^T\Sigma^{-1}+\frac{1}{2}\sum_{\vec{x}_i\in C_2}\Sigma^{-1}(\vec{x}_i-\vec{\mu}_2)(\vec{x}_i-\vec{\mu}_2)^T\Sigma^{-1} \end{aligned},令其为0得\begin{aligned} \hat{\Sigma} =\frac{1}{N}\left [\sum_{\vec{x}_i\in C_1}(\vec{x}_i-\vec{\mu}_1)(\vec{x}_i-\vec{\mu}_1)^T+\sum_{\vec{x}_i\in C_2}(\vec{x}_i-\vec{\mu}_2)(\vec{x}_i-\vec{\mu}_2)^T \right ]=\frac{N_1\cdot S_1+N_2\cdot S_2}{N} \end{aligned},其中S_1=\frac{1}{N_1}\sum_{\vec{x}_i\in C_1}(\vec{x}_i-\vec{\mu}_1)(\vec{x}_i-\vec{\mu}_1)^TS_2=\frac{1}{N_2}\sum_{\vec{x}_i\in C_2}(\vec{x}_i-\vec{\mu}_2)(\vec{x}_i-\vec{\mu}_2)^T,分别为两个类内样本方差;

4.6. 朴素贝叶斯 (Naive Bayes Classifier)

4.6.1. 基本思想

  • 所有朴素贝叶斯家族的算法都是基于朴素贝叶斯假设,又叫条件随机场假设,即假设各个特征之间相互独立;
  • 朴素贝叶斯模型是最简单的概率图模型;
  • 模型方法和高斯判别分析较为接近,这里不做重复;

5. 降维

5.1. 降维的意义

5.1.1. 过拟合问题

  • 通常模型通过训练集数据进行训练,若其测试集的效果明显低于训练集,这是一种不理想的效果,称为过拟合;
  • 过拟合常采用三种方法应对:
    • 采集更多数据,提高模型泛化能力;
    • 正则化,抑制模型拟合能力;
    • 降维,提取数据中的有效信息;
      • 直接降维:特征选择
      • 线性降维:PCA
      • 非线性降维:流形

5.1.2. 维度灾难

  • 从几何角度,若一组数据的特征维度是2,假设两个特征的值都是有上界的,则其构成的样本空间近似一个正方(长)形,数据在这个样本空间中分布;
  • 现在要对这组数据进行分类任务,模型映射到样本空间中变成了一条分隔线,将两个类别的样本分隔到了线的两边;
  • 为了便于从几何角度直观理解维数灾难,我们假设现在的模型分割线是一个圆,则这个分类模型的示意图如下:

  • 将两个特征归一化,样本空间变为了一个边长为1的正方形,面积为1,模型变成了一个半径小于等于\frac{1}{2}的圆,其包含区域面积为\pi r^2,占总空间比例小于\leqslant \pi/4
  • 而当模型维度变高时,例如维度为3,此时归一化的样本空间变成了一个三维立方体,体积为1,模型变成了一个三维球面,其包含区域体积为\frac{4\pi}{3}r^3,占总空间比例小于\leqslant \pi/6
  • 当维度继续升高时,模型表示的这个高维超球面包含区域所占总样本空间的超立方体比例将会越来越小,当维度非常大(趋于无穷)时,超球面包含区域将接近0,;
  • 分类边界是一个圆/球面,这是模型内部算法所决定的,而无论这个圆/球面的半径有多大,维度的升高都会导致分类边界仅仅只是样本空间的一个小点,这样的小点也必然无法有好的性能;
  • 现在将分类边界也泛化为任意几何体,模型越复杂,通常这个几何体越复杂,但当模型算法确定后,必然存在某个适宜维度使得该几何体可以在样本空间中很好地完成分类任务,而当维度远远大于该适宜维度时,模型都会不再适用;

5.2. 预备统计知识

  • 现有数据

Data:X=(\vec{x}_1, \vec{x}_2, \cdots, \vec{x}_n)^T_{N\times p}=\begin{pmatrix} \vec{x}_1^T\\ \vec{x}_2^T\\ \vdots \\ \vec{x}_N^T \end{pmatrix}

  • 均值为\bar{\vec{x}}=\frac{1}{N}\sum_{i=1}^N\vec{x}_i=\frac{1}{N}X^T\cdot 1_N,其中1_N=(1, 1, \cdots,1)^T是一个全一列向量;
  • 样本方差为S=\frac{1}{N}\sum_{i=1}^N(\vec{x}_i-\bar{\vec{x}})(\vec{x}_i-\bar{\vec{x}})^T,并经过如下推导:

\begin{aligned} S&=\frac{1}{N}\sum_{i=1}^N(\vec{x}_i-\bar{\vec{x}})(\vec{x}_i-\bar{\vec{x}})^T\\ &=\frac{1}{N}(\vec{x}_1-\bar{\vec{x}}, \vec{x}_2-\bar{\vec{x}}, \cdots, \vec{x}_N-\bar{\vec{x}})\begin{pmatrix} \vec{x}_1^T-\bar{\vec{x}}^T\\ \vec{x}_2^T-\bar{\vec{x}}^T\\ \vdots\\ \vec{x}_N^T-\bar{\vec{x}}^T \end{pmatrix}\\ &=\frac{1}{N}\left (X^T-\bar{\vec{x}}\cdot1_N^T \right )\left (X^T-\bar{\vec{x}}\cdot1_N^T \right )^T\\ &=\frac{1}{N}\left ( X^T-\frac{1}{N}X^T1_N1_N^T \right )\left ( X^T-\frac{1}{N}X^T1_N1_N^T \right )^T\\ &=\frac{1}{N}X^T\left ( I_N-\frac{1}{N}1_N1_N^T \right )\left ( I_N-\frac{1}{N}1_N1_N^T \right )^TX \end{aligned}\begin{aligned} \end{aligned}

  •  令H=I_N-\frac{1}{N}1_N1_N^T称为中心矩阵,则有S=\frac{1}{N}X^THH^TX,易知H^T=H,因此逆推可得H^TX=HX=(I_N-\frac{1}{N}1_N1_N^T)X=X-1_N\bar{\vec{x}}^T=\begin{pmatrix} \vec{x}_1^T-\bar{\vec{x}}^T\\ \vec{x}_2^T-\bar{\vec{x}}^T\\ \vdots \\ \vec{x}_N^T-\bar{\vec{x}}^T \end{pmatrix},因此对于任意数据X_{N\times p}HX即为对其进行中心化(每个样本减去均值);
  • H^2=I_N-\frac{2}{N}1_N1_N^T+\frac{1}{N^2}1_N1_N^T1_N1_N^T=I_N-\frac{2}{N}1_N1_N^T+\frac{1}{N}1_N1_N^T=H,因此可进一步化简S=\frac{1}{N}X^THX

5.3. 主成分分析(Principal Component Analysis,PCA)

5.3.1. 核心思想

  • 一个中心:原始特征空间的重构(特征之间从相关到无关),即找到一组线性无关的基,可以将原始特征空间的信息尽可能保留地投影到新空间,新空间特征维数低于原空间,这些基被称为主成分;
  • 两个基本点:
    • 最大投影方差;
    • 最小重构距离;

5.3.2. 从最大投影方差看PCA

  • 设一个投影方向为\vec{\mu}_1(单位向量),则某个样本\vec{x}_i\vec{\mu}_1方向的投影值为\vec{x}_i^T\cdot \vec{\mu}_1,因此可得,只需找到q个基(相互线性无关)U_{p\times q}=(\vec{\mu}_1, \vec{\mu}_2, \cdots, \vec{\mu}_q)q<p,则将样本\vec{x}_i投影到新特征空间的坐标为\vec{z}_i^T=\vec{x}_i^TU=\vec{x}_i^T\cdot (\vec{\mu}_1, \vec{\mu}_2, \cdots, \vec{\mu}_q)=(\vec{x}_i^T\vec{\mu}_1, \vec{x}_i^T\vec{\mu}_2, \cdots, \vec{x}_i^T\vec{\mu}_q),因此特征变换可以记作Z_{N\times q}=X_{N\times p}U_{p\times q},通过这个线性映射即可完成降维;
  • 单看一个投影方向\vec{\mu}_1,其投影结果为\vec{h}_1=X\vec{\mu}_1,投影均值为\bar{h}_1=\frac{1}{N}\sum_{i=1}^Nh_{1i}=\frac{1}{N}\sum_{i=1}^N\vec{x}_i^T\vec{\mu}_1=\bar{\vec{x}}^T\vec{\mu}_1,投影方差为

\begin{aligned} S_{\vec{h}_1}&=\frac{1}{N}\sum_{i=1}^N(h_{1i}-\bar{h}_1)^2\\ &=\frac{1}{N}\sum_{i=1}^N(\vec{x}_i^T\vec{\mu}_1-\bar{\vec{x}}^T\vec{\mu}_1)^2\\ &=\frac{1}{N}\sum_{i=1}^N\vec{\mu}_1^T(\vec{x}_i-\bar{\vec{x}})(\vec{x}_i-\bar{\vec{x}})^T\vec{\mu}_1\\ &=\vec{\mu}_1^TS\vec{\mu}_1 \end{aligned}

  • 为了最大化投影方差,将寻找\vec{\mu}_1的问题定义为一个条件优化问题,即\hat{\vec{\mu}}_1=\underset{\vec{\mu}_1}{argmax}(\vec{\mu}_1^TS\vec{\mu}_1)s.t.\vec{\mu}_1^T\vec{\mu}_1=1,利用拉格朗日算子法得\pounds (\vec{\mu}_1,\lambda)=\vec{\mu}_1^TS\vec{\mu}_1+\lambda(1-\vec{\mu}_1^T\vec{\mu} _1),求导得\frac{\partial \pounds }{\partial \vec{\mu}_1}=2S\vec{\mu}_1-\lambda\cdot 2\vec{\mu}_1,令其为0可得S\vec{\mu}_1=\lambda\vec{\mu}_1,因此\vec{\mu}_1是原样本方差S的特征向量时,投影方差可以被最大化;
  • 因此可以求出投影方差S_{\vec{h}_1}=\vec{\mu}_1^TS\vec{\mu}_1=\lambda_1\vec{\mu}_1^T\vec{\mu}_1=\lambda_1,也就是\vec{\mu}_1对应的特征值即为投影方差,因此PCA就是对原样本方差矩阵S求特征值和特征向量,由于实对称矩阵必可相似对角化,所以必然可以找到p个线性无关的特征向量,将其按照特征值大小排序,取前q个特征值最大的特征向量作为投影方向,每个样本通过这q个投影方向算得新的坐标值即为新的特征,构成了新的q维特征空间,实现了降维;

5.3.3. 从最小重构距离看PCA

  • 先看看最小化投影距离:从几何角度,设某个样本\vec{x}_i和投影方向\vec{\mu}_1之间的夹角为\theta,则投影值\vec{x}_i^T\cdot \vec{\mu}_1=\left \| \vec{x}_i \right \|\cdot \cos \theta,而投影距离为\left \| \vec{x}_i \right \|\cdot \sin \theta=\sqrt{\left \| \vec{x}_i \right \|^2-(\vec{x}_i^T\cdot \vec{\mu}_1)^2},因为不同投影方向的坐标值大小不同,因此为了将不同投影方向的情况进行公平对比,可以进行零均值化,此时投影距离的平方为(\vec{x}_i-\bar{\vec{x}})^T (\vec{x}_i-\bar{\vec{x}})-\vec{\mu}_1^T (\vec{x}_i-\bar{\vec{x}}) (\vec{x}_i-\bar{\vec{x}})^T\vec{\mu}_1,而优化目标为

\begin{aligned} \hat{\vec{\mu}}_1&=\underset{\vec{\mu}_1}{argmin}\sum_{i=1}^N\left [(\vec{x}_i-\bar{\vec{x}})^T (\vec{x}_i-\bar{\vec{x}})-\vec{\mu}_1^T (\vec{x}_i-\bar{\vec{x}}) (\vec{x}_i-\bar{\vec{x}})^T\vec{\mu}_1 \right ]\\ &=\underset{\vec{\mu}_1}{argmin}(-\vec{\mu}_1^T S\vec{\mu}_1)\underset{\vec{\mu}_1}{argmax}(\vec{\mu}_1^T S\vec{\mu}_1) \end{aligned}

        同样约束为s.t.\vec{\mu}^T\vec{\mu}=1,与最大投影方差方法完全一致;

  • 从重构角度,一个向量视为该向量空间内基向量的线性组合,例如\vec{x}=(1,2)=1\cdot(1,0)+2\cdot(0,1)=(1,2)\begin{pmatrix} (1,0)\\ (0,1) \end{pmatrix}= (h_1,h_2)\begin{pmatrix} \vec{v}_1\\ \vec{v}_2 \end{pmatrix},但在另一个向量空间中由不同的基向量会得到不同的坐标表示,例如\begin{aligned}\vec{x}=(1,2)= (h_1',h_2')\begin{pmatrix} \vec{v}_1'\\ \vec{v}_2' \end{pmatrix}= (\sqrt{2},1)\begin{pmatrix} (\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})\\ (0, 1) \end{pmatrix}=\sqrt{2}\cdot (\frac{1}{\sqrt{2}}, \frac{1}{\sqrt{2}})+1\cdot (0, 1) \end{aligned},因此在上文的投影中,通过找到q个投影方向可以将样本\vec{x}_i投影到新的q维向量空间,坐标为\vec{z}_i^T=(\vec{x}_i^T\vec{\mu}_1, \vec{x}_i^T\vec{\mu}_2, \cdots, \vec{x}_i^T\vec{\mu}_q),因此可以用新空间的坐标和基向量将\vec{x}_i表示为\hat{\vec{x}}_i=U\vec{z}_i,当q<p即降维时,基向量信息不完整,即\hat{\vec{x}}_i\neq \vec{x}_i,无法完全重构,当q=p时可以完全重构;(其实上文中\vec{z}_i^T=\vec{x}_i^TU\Rightarrow \vec{z}_i=U^T\vec{x}_i,其中U由实对称矩阵S的单位特征向量组成,所以当p=q时有U^{-1}=U^T,因此\vec{z}_i=U^T\vec{x}_i=U^{-1}\vec{x}_i\Rightarrow\vec{x}_i=U\vec{z}_i可推成立);
  • 因此也可以理解为,当U是按照特征值排序的前q\vec{\mu}_i时,可以表示为\begin{aligned} \hat{\vec{x}}_i=U\vec{z}_i=\sum_{k=1}^q\left (\vec{\mu}_k\cdot z_{ik} \right )=\sum_{k=1}^q\left [\vec{\mu}_k\cdot \left (\vec{\mu}_k^T\vec{x}_i \right )\right ]=\sum_{k=1}^q\left (\vec{\mu}_k \vec{\mu}_k^T \right )\vec{x}_i\end{aligned},同理\begin{aligned}\vec{x}_i=\sum_{k=1}^p\left [\vec{\mu}_k\cdot \left (\vec{\mu}_k^T\vec{x}_i \right )\right ]\end{aligned}\begin{aligned}\left (\vec{x}_i -\hat{\vec{x}_i} \right )=\sum_{k=q+1}^p\left [\vec{\mu}_k\cdot \left (\vec{\mu}_k^T\vec{x}_i \right )\right ]\end{aligned}
  • 因为降维任务会导致\hat{\vec{x}}_i\neq \vec{x}_i,因此在给定q的情况下,PCA的优化目标也可以理解为最小化重构代价,即

\begin{aligned} \hat{U}_{p\times q}&=\underset{U}{argmin}\frac{1}{N}\sum_{i=1}^N\left \| \hat{\vec{x}}_i - \vec{x}_i \right \|^2\\ &=\underset{U}{argmin}\frac{1}{N}\sum_{i=1}^N\left \| \sum_{k=q+1}^p\left [\vec{\mu}_k\cdot \left (\vec{\mu}_k^T\vec{x}_i \right )\right ] \right \|^2\\ &=\underset{U}{argmin}\frac{1}{N}\sum_{i=1}^N\left [\sum_{k=q+1}^p\left (\vec{\mu}_k^T\vec{x}_i \right )^2\left \| \vec{\mu}_k\cdot \right \|^2 \right ]\\ &=\underset{U}{argmin}\frac{1}{N}\sum_{i=1}^N\left [\sum_{k=q+1}^p\left (\vec{\mu}_k^T\vec{x}_i \right )^2 \right ] \end{aligned}

  • 因为不同投影方向的坐标值大小不同,因此为了将不同投影方向的情况进行公平对比,可以进行零均值化,即\begin{aligned} \hat{U}_{p\times q}&=\underset{U}{argmin}\frac{1}{N}\sum_{i=1}^N\left [\sum_{k=q+1}^p\left [\vec{\mu}_k^T(\vec{x}_i-\bar{\vec{x}}_i)(\vec{x}_i-\bar{\vec{x}}_i)^T\vec{\mu}_k \right ] \right ]=\underset{U}{argmin}\sum_{k=q+1}^p\left [\vec{\mu}_k^TS\vec{\mu}_k \right ]\\ \end{aligned}
  • 因此从最大化投影方差角度,PCA是在选择q个特征值最大的特征向量,从最小化重构距离的角度,PCA是在选择p-q个特征值最小的特征向量,将其舍弃,两个问题的等价的;

5.3.4. 从SVD角度看PCA

  • 从5.2.中可知,S=\frac{1}{N}X^THH^TX,若将H^TX作奇异值分解,即H^TX=V\Sigma U^T,则S=\frac{1}{N}U\Sigma^T V^TV\Sigma U^T=U(\frac{1}{N}\Sigma^T\Sigma) U^T,因此对S进行特征值分解可以得到投影方向UU也可以是H^TX进行奇异值分解得到的结果,而由奇异值分解的性质(1.2.3.)可知,最大的一小部分奇异值之和就占了所有奇异值之和的99%,因此取前几个奇异值最大的投影方向可以最大程度保留原数据的信息,也就是\hat{\vec{x}}_i\approx \vec{x}_i
  • 又由奇异值分解的第一个性质(1.2.3.)可得,H^TXU=V\Sigma=(\sigma_1\cdot\vec{v}_1, \sigma_2\cdot\vec{v}_2,\cdots,\sigma_p\cdot\vec{v}_p)_{N\times p},其中\sigma_i\vec{v}_i=H^TX\vec{u}_i即为所有样本从原始特征空间经方向\vec{u}_i投影的坐标值(等同于\vec{z}_i^T=\vec{x}_i^TU=\vec{x}_i^T\cdot (\vec{\mu}_1, \vec{\mu}_2, \cdots, \vec{\mu}_q)=(\vec{x}_i^T\vec{\mu}_1, \vec{x}_i^T\vec{\mu}_2, \cdots, \vec{x}_i^T\vec{\mu}_q),左乘中心矩阵H^T便得到相同结果),所以奇异值分解的矩阵V其实是坐标矩阵,因此可以再令T=\frac{1}{N}H^TXX^TH=\frac{1}{N}V\Sigma U^TU\Sigma^T V^T=V\left (\frac{1}{N}\Sigma \Sigma^T \right ) V^T,经过特征值分解可以得到坐标矩阵V,这个过程被称为主坐标分析(Principal Coordinate Analysis,PCoA);

6. 支持向量机(Support Vector Machine,SVM)

6.1. 相关概念

6.1.1. 分类

  • 硬间隔(Hard-margin)SVM,又叫最大间隔分类SVM;
  • 软间隔(Soft-margin)SVM;
  • 核(Kernel)SVM;

6.2. 硬间隔SVM

6.2.1. 问题定义:

  • 之前的分类算法中,模型致力于寻找一个超平面,可以将训练集中的两类样本分开,但这个超平面理论上存在无数个,它们在训练集上的分类效果可能是相同的,但模型的鲁棒性是不同的,为了寻找最鲁邦的模型,SVM便想要寻找一个所有样本总间隔最大的超平面;
  • 首先定义间隔:假设通过一个超平面f(\vec{x})=sign\left (\vec{w}^T\vec{x}+b \right ),可以完美拟合训练集所有样本,即对任意样本i,都有y_i(\vec{w}^T\vec{x}_i+b)>0,则所有样本中最小的间隔为margin(\vec{w})=\underset{\vec{w},\vec{x}_i}{\min}\frac{1}{\left \| \vec{w} \right \|}\left | \vec{w}\vec{x}_i+b \right |
  • 因此硬间隔SVM的优化目标即为最大化这个间隔,使得分类边界具有最高的鲁棒性,即\begin{aligned}\underset{\vec{w}}{\max} \ margin(\vec{w})=\underset{\vec{w}}{\max}\ \underset{\vec{x}_i}{\min}\frac{1}{\left \| \vec{w} \right \|}\left | \vec{w}^T\vec{x}_i+b \right |=\underset{\vec{w},b}{\max}\ \frac{1}{\left \| \vec{w} \right \|}\underset{\vec{x}_i}{\min}\left | \vec{w}^T\vec{x}_i +b\right | \end{aligned}

6.2.2. 问题化简:

  • 因为y_i\in \{-1,1\},所以可以写作\left | \vec{w}^T\vec{x}_i+b \right |=y_i(\vec{w}^T\vec{x}_i+b),因此优化目标可以写作\underset{\vec{w},b}{\max} \ margin(\vec{w})=\underset{\vec{w},b}{\max}\ \frac{1}{\left \| \vec{w} \right \|}\underset{\vec{x}_i}{\min}y_i(\vec{w}^T\vec{x}_i+b )
  • 又因为对超平面\vec{w},我们只关心其方向,对\left \| \vec{w} \right \|无要求,因此对\vec{w}^T\vec{x}_i+b我们只关心正负,无所谓值的大小,因此对y_i(\vec{w}^T\vec{x}_i+b)我们只关心是否大于0,无所谓值的大小,因此为了简化优化目标,可以令\underset{\vec{w},\vec{x}_i,b}{\min}\ y_i( \vec{w}^T\vec{x}_i+b )=1,则优化目标变为\underset{\vec{w}}{\max} \ margin(\vec{w})=\left\{\begin{matrix} \underset{\vec{w}}{\max}\ \frac{1}{\left \| \vec{w} \right \|}\cdot 1=\underset{\vec{w}}{\min}\left \| \vec{w} \right \|=\underset{\vec{w}}{\min}\ \vec{w}^T\vec{w}\\ s.t.\ \ y_i(\vec{w}^T\vec{x}_i +b)\geqslant 1 \end{matrix}\right.,变为一个凸优化问题;

6.2.3. 问题求解:

  • 优化问题首先考虑使用拉格朗日乘数法:
    • \pounds (\vec{w},\vec{\lambda},b)=\vec{w}^T\vec{w}+\sum_{i=1}^N\lambda_i\left [ 1-y_i(\vec{w}^T\vec{x}_i+b) \right ],将带约束的原问题变为无约束优化问题:\left\{\begin{matrix} \underset{\vec{w},b}{\min} \ \underset{\vec{\lambda}}{\max} \ \pounds (\vec{w},\vec{\lambda},b)\\ s.t. \lambda_i\geqslant 0 \end{matrix}\right.
    • 若原问题的约束y_i(\vec{w}^T\vec{x}_i+b )\geqslant 1不满足,即1-y_i(\vec{w}^T\vec{x}_i+b )> 0,则\underset{\vec{\lambda}}{\max} \ \pounds (\vec{w},\vec{\lambda},b)\rightarrow \infty,反之若都满足,则\underset{\vec{\lambda}}{\max} \ \pounds (\vec{w},\vec{\lambda},b)=\vec{w}^T\vec{w},因此优化目标其实是\left\{\begin{matrix} \underset{\vec{w}}{\min} (\infty ,\vec{w}^T\vec{w})\\ s.t. \lambda_i\geqslant 0 \end{matrix}\right.,因此这个优化目标函数为了取到\underset{\vec{w}}{\min} (\infty ,\vec{w}^T\vec{w})=\vec{w}^T\vec{w},会调整\vec{w}使得y_i(\vec{w}^T\vec{x}_i+b )\geqslant 1对所有\hat{\vec{x}}_i都满足,由此原问题的约束成立;
    • 而此时\pounds (\vec{w},\vec{\lambda},b)=\vec{w}^T\vec{w}+\sum_{i=1}^N\lambda_i\left [ 1-y_i(\vec{w}^T\vec{x}_i+b) \right ]=\vec{w}^T\vec{w},因此对任意i,有\lambda_i\left [ 1-y_i(\vec{w}^T\vec{x}_i+b) \right ]=0,则\lambda_i=01-y_i(\vec{w}^T\vec{x}_i+b) =0,由于绝大多数情况,只有两个点离超平面最近(存在更多点距离超平面也是0的理论可能,但实际数据中几乎不存在),即1-y_i(\vec{w}^T\vec{x}_i+b) =0,其他情况都是\lambda_i=0(但这里不能直接带入\lambda_i=0,因为这个是目标函数为了达到优化而隐含的约束,如果直接带入则变成了优化之间的已知约束,问题会变得无法求解);

    • 这两个最近的样本我们称其为\vec{x}_my_m=1)和\vec{x}_ny_n=-1),这两个向量被称为支持向量;
  • 变为对偶问题:
    • 弱对偶关系:\min \ \max F \geqslant \max \ \min F(宁做凤尾,不做鸡头),当中间取等号时变为强对偶关系,而线性约束的二次优化问题天然满足强对偶性(此处不做数学证明);
    • 因此这里可以将优化目标变为\left\{\begin{matrix} \ \underset{\vec{\lambda}}{\max} \ \underset{\vec{w},b}{\min} \ \pounds (\vec{w},\vec{\lambda},b)\\ s.t. \lambda_i\geqslant 0 \end{matrix}\right.
  • 求偏导得到结果:
    • \frac{\partial \pounds }{\partial b}=-\sum_{i=1}^N\lambda_iy_i=0
    • \frac{\partial \pounds }{\partial \vec{w}}=2\vec{w}-\sum_{i=1}^N\lambda_iy_i\vec{x}_i=0\Rightarrow\vec{w}^*= \frac{1}{2}\sum_{i=1}^N\lambda_iy_i\vec{x}_i
    • \begin{aligned} \pounds (\vec{w},\vec{\lambda},b)&=\frac{1}{4}\sum_{i=1}^N\sum_{j=1}^N\lambda_i\lambda_jy_iy_j(\vec{x}_i^T\vec{x}_i )-\frac{1}{2}\sum_{i=1}^N\sum_{j=1}^N\lambda_i\lambda_jy_iy_j(\vec{x}_i^T\vec{x}_i )+\sum_{i=1}^N\lambda_i\\ &=-\frac{1}{4}\sum_{i=1}^N\sum_{j=1}^N\lambda_i\lambda_jy_iy_j(\vec{x}_i^T\vec{x}_i )+\sum_{i=1}^N\lambda_i \\ \end{aligned}
  • 因为对\left\{\begin{matrix} \ \underset{\vec{\lambda}}{\max} \ \underset{\vec{w},b}{\min} \ \pounds (\vec{w},\vec{\lambda},b)\\ s.t. \lambda_i\geqslant 0 \end{matrix}\right.的优化必须使\vec{\lambda}\vec{w}满足:只有两个点\vec{x}_m\vec{x}_n离超平面最近,即1-y_i(\vec{w}^T\vec{x}_i+b) =0,其他情况都是\lambda_i=0,将这个条件带入\sum_{i=1}^N\lambda_iy_i=0可得:\lambda_m=\lambda_n,\ \ \lambda_{i\neq m,n}=0
  • 带入目标函数得\begin{aligned} \pounds (\vec{w},\vec{\lambda},b)&=-\frac{1}{4}\left ( \lambda_m^2(\vec{x}_m^T\vec{x}_m )+\lambda_n^2(\vec{x}_n^T\vec{x}_n )-2\lambda_m\lambda_n(\vec{x}_m^T\vec{x}_n ) \right )+\lambda_m+\lambda_n=-\frac{1}{4}\lambda_m^2\left \|\vec{x}_m-\vec{x}_n \right \|+2\lambda_m\end{aligned},对其求偏导\frac{\partial \pounds }{\partial \lambda_m}=-\frac{1}{2}\lambda_m\left \| \vec{x}_m-\vec{x}_n \right \|+2=0\Rightarrow \lambda_m=\frac{4}{\left \| \vec{x}_m-\vec{x}_n \right \|},带入得\begin{aligned} \pounds (\vec{w},\vec{\lambda},b)=-\frac{1}{4}\lambda_m^2\left \|\vec{x}_m-\vec{x}_n \right \|+2\lambda_m=\frac{4}{\left \| \vec{x}_m-\vec{x}_n \right \|}\end{aligned},因此\max\pounds (\vec{w},\vec{\lambda},b)=\min\left \| \vec{x}_m-\vec{x}_n \right \|,找到两个类中,间距最小的两个样本,计算其欧氏距离可以算得\vec{\lambda}
  • 再带入\vec{w}^*= \sum_{i=1}^N\lambda_iy_i\vec{x}_i=\lambda_m(\vec{x}_m-\vec{x}_n)b^*=y_m-\lambda_m\vec{x}_m^T(\vec{x}_m-\vec{x}_n)即可得到分类边界;

6.3. 软间隔SVM

6.3.1. 引入软间隔的思想

  • 硬间隔SVM假设所有样本都可以被完美分类,但实际情况很难做到,因此需要允许一些错误,但又要尽可能控制错误,便有了软间隔SVM的思想;
  • 相比于硬间隔,软间隔SVM取消了约束s.t.\ \ y_i(\vec{w}^T\vec{x}_i +b)\geqslant 1,而是采用一个损失函数来代替,即设定一组参数\vec{\xi}=\{ \xi_1, \xi_2, \cdots, \xi_N\}, \ \xi_i \geqslant 0,目标函数变为\left\{\begin{matrix} \underset{\vec{w}, b, \vec{\xi}}{\min}\ \vec{w}^T\vec{w}+C\sum_{i=1}^N\xi_i\\ s.t. \ y_i(\vec{w}^T\vec{x}_i+b)\geqslant 1-\xi_i, \ \xi_i\geqslant 0 \end{matrix}\right.
  • 通过这个目标函数,原本间隔外的点(y_i(\vec{w}^T\vec{x}_i+b)\geqslant 1)的参数\xi_i会被优化到0,而间隔内的点无论是否可以被分类正确,都可以尽可能使其向两边分散;
  • C是一个乘法因子,可以作为模型超参数,其作用有点像正则化,C的值越大,\vec{\xi}被惩罚得越厉害,间隔小的点越少;

6.3.2. 软间隔SVM的求解

  • 前面的部分和硬间隔相同,化为无约束目标函数和对偶优化问题,分别对\vec{w}b求偏导后,最终目标函数化简如下:\left\{\begin{matrix} \ \underset{\vec{\lambda}}{\min} \ \frac{1}{4}\sum_{i=1}^N\sum_{j=1}^N\lambda_i\lambda_jy_iy_j(\vec{x}_i^T\vec{x}_i )-\sum_{i=1}^N\lambda_i\\ s.t. \ 0\leqslant \lambda_i\leqslant C ,\ \sum_{i=1}^N\lambda_iy_i=0\end{matrix}\right.\vec{w}^*= \frac{1}{2}\sum_{i=1}^N\lambda_iy_i\vec{x}_i
  • 同样先求得\vec{\lambda}后(这里需要用到SMO迭代求解算法,过于复杂,后补),即可得到超平面;

7. 核方法

7.1. 背景

  • 对于严格线性可分问题,感知机、SVM等算法都可以提供很好的解决方案,但对于非线性分类问题,需要一个非线性转换来将其变为线性问题;
  • 高维空间往往比低维空间更容易线性可分,因此可以通过一些非线性变换将已知维度计算出新的维度,便有可能将非线性问题变为线性问题;
  • 然而在实际问题中,可以将非线性问题转化为线性问题的变换函数\phi (\vec{x})可能是非常复杂的,而根据概率派思想,将一个问题最终看做一个优化问题,则往往需要计算内积,如\left (\phi(\vec{x}_i) - \phi(\bar{\vec{x}}) \right )^T\left (\phi(\vec{x}_j) - \phi(\bar{\vec{x}}) \right ),这会给优化带来非常可怕的复杂度;
  • 事实上,我们所关心的是最终的内积的值,而不需要\phi (\vec{x})的结果,因此如果有一个函数可以一步得到内积,即K(\vec{x}_i, \vec{x}_j) \equiv \left (\phi(\vec{x}_i) - \phi(\bar{\vec{x}}) \right )^T\left (\phi(\vec{x}_j) - \phi(\bar{\vec{x}}) \right ),无需计算\phi (\vec{x})即可直接得到内积结果,就可以很大程度简化优化的难度,因此核方法诞生了,K(\vec{x}_i, \vec{x}_j)即为(正定)核函数;

7.2. 定义

  • 核(无需过多关注):若存在函数K使得K:X \times X\rightarrow \mathbb{R},其中X为特征空间,\mathbb{R}为全体实数域,则称K为核函数;
  • 正定核:对于一个核函数K,若存在非线性变换\phi (\vec{x})\phi (\vec{x})属于希尔伯特空间,希尔伯特空间指完备的、可能是无限维的、被赋予内积的线性空间),使得K(\vec{x}_i, \vec{x}_j) \equiv \phi(\vec{x}_i) ^T\phi(\vec{x}_j),换言之任意N个样本的Gram矩阵\begin{bmatrix} K(\vec{x}_1, \vec{x}_1) & K(\vec{x}_1, \vec{x}_2) & \cdots & K(\vec{x}_1, \vec{x}_N)\\ K(\vec{x}_2, \vec{x}_1) & K(\vec{x}_2, \vec{x}_2) & \cdots & K(\vec{x}_2, \vec{x}_N)\\ \vdots & \vdots & \ddots &\vdots \\ K(\vec{x}_N, \vec{x}_1) & K(\vec{x}_N, \vec{x}_2) & \cdots & K(\vec{x}_N, \vec{x}_N) \end{bmatrix}是半正定的,则称K为正定核函数;

8. 指数族分布

8.1. 背景

  • 指数族分布的定义式p(\vec{x}\ |\ \vec{\eta})=h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x})-A(\vec{\eta}) \right \},其中:
    • \vec{\eta} 为参数向量;
    • A(\vec{\eta}) 为对数配分函数(log partition function),配分函数表示包含了所有状态信息的函数,常用于归一化等用途,例如将某函数\hat{p}(\vec{x}\ |\ \vec{\theta})归一化为概率分布,可以设z=\int \hat{p}(\vec{x}\ | \ \vec{\theta}),则p(\vec{x} \ | \ \vec{\theta})=\frac{1}{z}\hat{p}(\vec{x}\ |\ \vec{\theta})可以使\int p(\vec{x} \ | \ \vec{\theta})=1,即被归一化为概率分布。在指数族分布中,p(\vec{x}\ |\ \vec{\eta})=h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x}) \right \}\cdot /\exp\left \{A(\vec{\eta}) \right \},则z = \exp\left \{ A(\vec{\eta})\right \},A(\vec{\eta})=\log z即为对数配方函数;
    • \phi (\vec{x})代表充分统计量,即对于一组数据\left \{ \vec{x}_i \right \}_{i\ \text{from 1 to}\ N}\phi (\vec{x})包含了其所有统计信息,例如\phi (\vec{x})=\left [ E(\vec{x}), D(\vec{x}) \right ]等等;
  • 指数族分布的特点(学完概率图模型和变分推断后补)

8.2. 高斯分布的指数族分布推导

  • 对于一组符合一维高斯分布的数据,设参数为\vec{\theta}=[\mu , \sigma^2],其概率密度函数为:\begin{aligned} p(x|\theta)&=\frac{1}{\sqrt{2\pi}\cdot\sigma}\exp\left \{ -\frac{(x-\mu)^2}{2\sigma^2} \right \}\\ &=\exp\left \{ -\frac{1}{2}\log(2\pi\cdot\sigma^2) \right \}\cdot \exp\left \{ -\frac{1}{2\sigma^2}(x^2-2\mu \cdot x) - \frac{\mu^2}{2\sigma^2} \right \}\\ &=\exp\left \{ \left [\begin{matrix} \frac{\mu}{\sigma^2} & -\frac{1}{2\sigma^2} \end{matrix} \right ]\cdot\left [ \begin{matrix} x\\ x^2 \end{matrix} \right ] - \left (\frac{\mu^2}{2\sigma^2} +\frac{1}{2}\log(2\pi\cdot\sigma^2) \right )\right \}\\ \end{aligned}
  • 因此可得:\vec{\eta}=\left [\begin{matrix} \frac{\mu}{\sigma^2} & -\frac{1}{2\sigma^2} \end{matrix} \right ]^T,\eta_1=\frac{\mu}{\sigma^2},\eta_2=-\frac{1}{2\sigma^2},同时\phi(x)=\left [ \begin{matrix} x & x^2 \end{matrix} \right ] ^TA(\vec{\eta})=\frac{\mu^2}{2\sigma^2} +\frac{1}{2}\log(2\pi\cdot\sigma^2) =-\frac{\eta_1^2}{4\eta_2}+\frac{1}{2}\log(-\frac{\pi}{\eta_2}),即得到了一维高斯分布的指数族分布形式;

8.3. 关于对数配分函数的一些结论

  • 由于p(\vec{x}\ |\ \vec{\eta})=h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x}) \right \}/\exp\left \{ A(\vec{\eta})\right \},可得\exp\left \{ A(\vec{\eta})\right \}\cdot p(\vec{x}\ |\ \vec{\eta})=h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x}) \right \}
  • 再两边积分可得:\int \exp\left \{ A(\vec{\eta})\right \}\cdot p(\vec{x}\ |\ \vec{\eta})d\vec{x}=\exp\left \{ A(\vec{\eta})\right \}=\int h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x}) \right \}d\vec{x}
  • 两边对\vec{\eta}求导为:\exp\left \{ A(\vec{\eta})\right \}\cdot A'(\vec{\eta})=\int h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x}) \right \}\cdot \phi(\vec{x}) d\vec{x}
  • 因此\begin{aligned} A'(\vec{\eta})&=\int \left [h(\vec{x})\cdot \exp\left \{\vec{\eta}^T\cdot \phi(\vec{x}) \right \}/\exp\left \{ A(\vec{\eta})\right \} \right ] \cdot \phi(\vec{x}) d\vec{x}\\&=\int p(\vec{x}\ |\ \vec{\eta}) \cdot \phi(\vec{x}) d\vec{x}=E_{p(\vec{x}\ |\ \vec{\eta})}[\phi(\vec{x}) ] \end{aligned}
  • 进一步也可推导得A''(\vec{\eta})=D_{p(\vec{x}\ |\ \vec{\eta})}[\phi(\vec{x}) ],推导过程略;

8.4. 最大熵角度解读指数族分布

  • 概率的信息量:-\log p
  • 熵的定义:信息量的期望,E_{p}[-\log p]=\int -p(x)\cdot \log p(x)dx \ \ \ \text{or} \ \ \ -\sum_xp(x)\log p(x)
  • 最大熵:令熵值最大的概率分布情况,通常分布越趋向于等可能/均匀分布,熵值越大,因此最大熵也是等可能的量化分析;
  • 对于一组离散变量的数据Data=\{ x_1, x_2, \cdots, x_N \},其分布未知,但可以用经验分布代替,即\hat{p}(x)=\frac{count(x)}{N},同时设一向量函数\vec{f}(x)=[\begin{matrix} f_1(x) & f_2(x) & \cdots & f_Q(x) \end{matrix}],则可以对其求期望E_{\hat{p}}[\vec{f}(x)]=\vec{\Delta },因此可以将最大熵问题设为一个优化问题,即\left\{\begin{matrix} min \sum_xp(x)\log p(x)\\ s.t.\ \sum_xp(x)=1,\ \ E_p[\vec{f}(x)]=E_{\hat{p}}[\vec{f}(x)]=\vec{\Delta} \end{matrix}\right.
  • 使用拉格朗日乘数法求解,\pounds (p,\lambda_0, \vec{\lambda})=\sum_xp(x)\log p(x)+\lambda_0[1-\sum_xp(x)]+\vec{\lambda}^T(\vec{\Delta}-E_p[\vec{f}(x)]),对p(x)求导得\frac{\partial \pounds }{\partial p}=\sum_x(\log p+1)-\sum_x\lambda_0-\sum_x\vec{\lambda}^T\vec{f}(x),令其为0,则解得p(x)=\exp\{ \vec{\lambda}^T\vec{f}(x)+\lambda_0-1 \},因此当\vec{f}(x)=\phi(x)时,可以令\vec{\eta}=\vec{\lambda}h(x)\cdot \exp\{ -A(\vec{\eta}) \}=\exp\{ \lambda_0-1 \},此时p(x)为指数族分布;

9. 概率图模型

9.1. 背景

  • 概率模型的四大法则:
    • 加法法则p(x_1)=\int p(x_1, x_2)dx_2
    • 乘法法则p(x_1, x_2)=p(x_1)\cdot p(x_2|x_1)=p(x_2)\cdot p(x_1|x_2)
    • 链式法则p(x_1, x_2, \cdots, x_p)=\prod _{i=1}^Np(x_i|x_1, \cdots, x_{i-1})
    • 贝叶斯法则p(x_2|x_1)=\frac{p(x_1, x_2)}{p(x_1)}=\frac{p(x_1, x_2)}{\int p(x_1, x_2)dx_2}=\frac{p(x_2)\cdot p(x_1| x_2)}{\int p(x_1, x_2)dx_2}
  • 概率模型的困境:当样本特征维度较高时,计算非常复杂;
  • 简化方法:
    • 朴素贝叶斯假设:假设特征之间相互独立,则p(x_1, x_2, \cdots, x_p)=\prod_{i=1}^N p(x_i),但该假设过强,很多时候无法保证成立;
    • 马尔可夫假设:假设每个特征变量只与前面的m个特征有关,与其他特征都独立,即p(x_i|x_1, x_2, \cdots, x_{i-1}, x_{i+1}, \cdots, x_p)=p(x_i| x_{i-1}, \cdots , x_{i-m}),称为m阶马尔可夫假设,但该假设依然过强,且依赖关系的建设较为随意,不符合真实分布;
    • 条件独立性假设:假设特征之间存在一组三个集合,\{ x_i \}_A, \{ x_i \}_B, \{ x_i \}_C,三个集合互不相交且无顺序性,使得\{ x_i \}_A\perp \{ x_i \}_B |\{ x_i \}_C,即在给定集合\{ x_i \}_C作为条件的情况下,\{ x_i \}_A\{ x_i \}_B内变量相互独立,条件独立性假设为了解决朴素贝叶斯假设和马尔可夫假设存在的问题而提出,是概率图模型的理论基础;
  • 概率图模型可以由以下分类:
    • 表示(Representation):
      • 有向图:贝叶斯网络;
      • 无向图:马尔可夫网络;
      • 高斯图(连续):有高斯贝叶斯网络和高斯马尔可夫网络;
    • 推断(Inference):
      • 精确推断;
      • 近似推断:确定性近似推断——变分推断;随机近似推断(MCMC);
    • 学习(Learning):
      • 参数学习:完备数据或隐变量学习(EM);
      • 结构学习;
  • 概率图模型,又叫贝叶斯网络,常见模型有(后续章节都会介绍):
    • 单一变量:朴素贝叶斯;
    • 混合变量:高斯混合模型(GMM);
    • 时间:
      • 马尔可夫链;
      • 高斯过程(无限维高斯分布);
      • 将时间贝叶斯网络与混合变量结合,即为动态模型:
        • 隐马尔可夫模型,离散时间;
        • 线性动态系统(LDS):卡曼滤波等(Kalman Filter),连续时间;
        • Particle Filter,非高斯,非线性;
    • 空间:高斯贝叶斯网络; 

9.2. 条件独立性的图结构表示

  • 下图是最简单的概率图模型,其概率模型可以写作p(a,b,c)=p(a)\cdot p(b|a)\cdot p(c|a),这就是条件独立性用图结构表示的方法,条件变量在父节点,独立变量是同层的无边节点,该过程也被称为因子分解;

  • 条件独立性的三种图结构:
    • T2T(tail to tail):如下图,概率模型已经写过是p(a,b,c)=p(a)\cdot p(b|a)\cdot p(c|a),可以理解为a,b,c三个随机变量相互之间可能有联系(如b=2\cdot a+1,\ c=a+5,此时b,c也有联系),也可能只有a,b之间、a,c之间有联系(如b=2\cdot x+a,\ c=y+5\cdot a,此时b,c没有联系,但一旦a被观测(值已定,不再是随机变量),且b,c依然是随机变量(刚刚说的第二种情况),则b,c独立(如果是第一种情况,那b,c是定值,也可以按照独立处理);

    • H2T(head to tail):如下图,概率模型为p(a,b,c)=p(b)\cdot p(a|b)\cdot p(c|b),若b被观测,则a,c独立

    • H2H(head to head):如下图,当c未被观测时,概率模型为p(a,b,c)=p(a,b)\cdot p(c|a,b)=p(a)\cdot p(b)\cdot p(c|a,b),此时a,b独立,但如果c被观测,此时c=f(a)=g(b),a=f^{-1}(g(b))a,b不再独立。值得注意的是,若c后继子节点被观测,同样会导致a,b连通,二者不再独立;

9.3. D-Seperation介绍 

  • 利用条件独立性的图结构来划分整个概率图模型的集合,即为D划分(D-Seperation)方法,具体来说,图的部分节点分别属于\{ x_i \}_A, \{ x_i \}_B, \{ x_i \}_C三个集合,三者互不相交且无顺序性,当满足下面两个条件时,使得\{ x_i \}_A\perp \{ x_i \}_B |\{ x_i \}_C,即在给定集合\{ x_i \}_C作为条件的情况下,\{ x_i \}_A\{ x_i \}_B内变量相互独立:
    • \{ x_i \}_A中节点a\{ x_i \}_C中节点cT2T或H2T结构存在时,其中间的连通节点b一定在\{ x_i \}_B中;
    • \{ x_i \}_A中节点a\{ x_i \}_C中节点cH2H结构存在时,其中间的连通节点b一定不在\{ x_i \}_B中;
  • 同时可以对某一个节点求其相关节点时,可以用全局节点作为条件求概率密度(设x_{-i}=x_1, x_2, \cdots, x_{i-1}, x_{i+1}, \cdots, x_p,即为除了x_i以外所有节点),即:\begin{aligned} p(x_i|x_{-i})&=\frac{p(x_i, x_{-i})}{p(x_{-i})}=\frac{p(\vec{x})}{\int p(\vec{x}) dx_i}=\frac{\prod _{j=1}^pp(x_j|x_{pa(j)})}{\int \prod _{j=1}^pp(x_j|x_{pa(j)}) dx_i} \\&= \frac{\prod _{x_j\text{ relate }x_i}p(x_j|x_{pa(j)}) \cdot \prod _{x_k\text{ not relate }x_i}p(x_k|x_{pa(k)})}{\int \prod _{x_j\text{ relate }x_i}p(x_j|x_{pa(j)})dx_i \cdot \int \prod _{x_k\text{ not relate }x_i}p(x_k|x_{pa(k)})dx_i}\\&= \frac{\prod _{x_j\text{ relate }x_i}p(x_j|x_{pa(j)}) \cdot \prod _{x_k\text{ not relate }x_i}p(x_k|x_{pa(k)})}{\int \prod _{x_j\text{ relate }x_i}p(x_j|x_{pa(j)})dx_i \cdot \prod _{x_k\text{ not relate }x_i}p(x_k|x_{pa(k)})}\\&= \frac{\prod _{x_j\text{ relate }x_i}p(x_j|x_{pa(j)})}{\int \prod _{x_j\text{ relate }x_i}p(x_j|x_{pa(j)})dx_i} \end{aligned}只留下了与x_i相关的节点;
  • x_i相关的节点如下图所示,只有x_i的父节点、子节点及其子节点的其他父节点与x_i有关,这种关系图被称为马尔可夫毯(Markov Blanket);

9.4. 马尔可夫随机场

  • 马尔可夫随机场(Markov random field),又叫概率无向图模型(probabilistic undirected graphical model),每个节点表示一个随机变量,边表示变量之间的相关性,可以用于表示联合概率分布;
  • 马尔可夫随机场有三种马尔可夫性质:
    • 成对马尔可夫性(Pair Markov Property):设\{x\}_{-i-j}为除了x_ix_j的其他所有节点,则当x_ix_j没有边连接时,则x_i\perp x_j \ | \ \{x\}_{-i-j}
    • 局部马尔可夫性(Local Markov Property):设\{x\}_{connect(i)}为与x_i连接的所有节点,则x_i\perp \left \{ x \right \}_{-i-connect(i)} \ | \ \{x\}_{connect(i)}
    • 全局马尔可夫性(Global Markov Property):设\{x\}_{A}为一个节点的集合,若存在三个集合则\{x\}_{A}\{x\}_{B}\{x\}_{C},其中\{x\}_{A}中任意节点与\{x\}_{B}中节点想要连接,都需要经过\{x\}_{C}中节点,则\left \{ x \right \}_A\perp \left \{ x \right \}_{B} \ | \ \{x\}_{C}
  • (clique):对于一组节点,如果任意两两之间都存在边连接,则称为一个团,若此时这个团无法再加入新的节点,加入任意一个节点都不能满足两两连通,则称为最大团(类比最大全连通子图);
  • Hammersley-Clifford定理:可以将马尔可夫随机场的联合概率分布表示为:p(\vec{x})=\frac{1}{Z}\prod _{C=1}^K \Psi _C(x_C),其中:
    • Z=\sum _{x_1}\sum _{x_2} \cdots \sum _{x_p}\prod _{C=1}^K\Psi (x_C)为一个权,用于将等号右侧化为概率值:
    • \Psi _C(x_C)是图模型中某个最大团{x}_C的势函数,势函数的概念由统计物理得来,可以表示为\Psi _C(x_C)=\exp\{ -E(x_C) \}E(x_C)为能函数,如此可以将概率密度化简为:p(\vec{x})=\frac{1}{Z}\prod _{C=1}^K \exp\{ -E(x_C) \}=\frac{1}{Z} \exp\{ -\sum _{C=1}^KE(x_C) \},该形式满足之前章节的指数族分布的形式,即可以说,马尔可夫随机场等价于指数族分布,等价于最大熵原理

9.5. 推断(Inference) 

  • 推断的本质工作就是求概率,或求分布,例如边缘概率、条件概率、后验概率等;
  • 推断有两种类型:
    • 精确推断:
      • Variable Elimination (VE),变量消除法,具体应用见17.3. CRF的边缘概率密度——Smoothing;
      • Belief Propagation (BP),又叫Sum-Product算法,由VE改进而来,是精确推断的重要算法,常用于树结构的推断;
      • Junction Tree Algorithm,是BP在普遍图结构上的拓展算法;
    • 近似推断:
      • Loop Belief Propagation,常用于有环图;
      • 蒙特卡罗推断(Mente Carlo Inference):例如Importance Sampling,MCMC等;
      • variational Inference

10. EM期望最大算法(expectation maximization)

10.1 算法简介

  • 概率模型求解时,很多时候数据只有观测变量,因为缺少隐变量而无法求得解析解,也很难通过极大似然估计法对完整的模型参数进行估计;
  • EM算法即为含有隐变量的极大似然估计法(或含有隐变量的极大后验概率估计法);
  • 举例:有并不均匀的 A,B,C 三个硬币,正面朝上的概率分别为 o,p,q ,首先抛掷 A ,若正面朝上就抛掷 B ,若反面朝上就抛掷 C ,记录最终的抛掷结果为 x_i=1\ or\ 0 ,但硬币 A 的抛掷结果作为隐变量未能观测,使得模型参数 o,p,q 不能精准求解;
  • 若使用极大似然估计法,则\begin{aligned}\theta=\underset{\theta}{\arg \max} \ \log P(X|\theta) \end{aligned},忽略了隐变量,但由于\begin{aligned} P(X|\theta)=\int_Z P(X,Z|\theta) dZ=\int_Z P(X|Z,\theta)\cdot P(Z|\theta) dZ=E_{Z|\theta}[P(X|Z,\theta)] \end{aligned},因此可以使用生成式概率模型的思想,认为X是由Z生成的,来对模型参数进行求解
  • EM算法分为两步,E步是对 P(X|\theta) 求期望,M步是迭代模型参数 \theta 使得期望最大;

10.2 算法收敛性证明

  • 首先直接看公式\begin{aligned}\theta^{(t+1)}=\underset{\theta}{\arg \max} \ E_{z|x,\theta^{(t)}}\left [\log P(x, z|\theta)\right ]=\underset{\theta}{\arg \max} \ \int_z\left [\log P(x, z|\theta)\cdot P(z|x,\theta^{(t)})\right ]dz \end{aligned},其中:
    • x 是样本数据,z是隐变量, \theta 是模型参数;
    • \log P(x, z|\theta) 被称为完备数据概率,因为其包含了 x,z 的联合概率分布;
    • p(z|x,\theta^{(t)}) 是后验概率;
  • 这个公式是通过调整参数 \theta ,使得 x,z 的完备数据的联合概率密度期望最大,因此被称为期望最大算法,核心思路类似于最大似然估计;
  • 具体来看,参数是逐步更新的,也就是每次给定当前的参数 \theta^{(t)},调整 \theta 使得完备数据概率的期望最大,并将其赋值给 \theta^{(t+1)} ,即完成了参数更新,直至收敛;
  • 接下来看收敛性证明:
    • 根据极大似然估计思想,若保证 \log P(x|\theta^{(t)})\leqslant \log P(x|\theta^{(t+1)}),则可以保证EM算法的有效性和收敛性;
    • 已知\log P(x|\theta)=\log \frac{P(x,z|\theta)}{P(z|x,\theta)}=\log P(x,z|\theta)-\log P(z|x,\theta),将两边对 z 积分:
      • 等号左边因为所有变量都与 z 无关,所以积分值为\log P(x|\theta)\cdot 1=\log P(x|\theta)
      • 等号右边为\int_z P(z|x,\theta^{(t)})\cdot \log P(x,z|\theta)dz-\int_z P(z|x,\theta^{(t)})\cdot \log P(z|x,\theta)dz,分别记作Q(\theta,\theta^{(t)})H(\theta,\theta^{(t)}),因此等号右边为Q(\theta,\theta^{(t)})-H(\theta,\theta^{(t)})
      • 若能保证Q(\theta^{(t+1)},\theta^{(t)})\geqslant Q(\theta,\theta^{(t)})H(\theta^{(t+1)},\theta^{(t)})\leqslant H(\theta,\theta^{(t)}),则可以保证\log P(x|\theta^{(t)})\leqslant \log P(x|\theta^{(t+1)}),从而保证期望上升;
        • Q(\theta,\theta^{(t)})即为定义(第一行的公式),因此Q(\theta^{(t+1)},\theta^{(t)})即为Q(\theta,\theta^{(t)})的最大值,因此Q(\theta^{(t+1)},\theta^{(t)})\geqslant Q(\theta,\theta^{(t)})恒成立;
        • \begin{aligned} &H(\theta^{(t+1)},\theta^{(t)}) - H(\theta^{(t)},\theta^{(t)})\\ =& \int_z P(z|x,\theta^{(t)})\cdot \log P(z|x,\theta^{(t+1)})dz-\int_z P(z|x,\theta^{(t)})\cdot \log P(z|x,\theta^{(t)})dz\\ =&\int_z P(z|x,\theta^{(t)})\cdot \log \frac{P(z|x,\theta^{(t+1)})dz}{P(z|x,\theta^{(t)})dz}=-KL\left (P(z|x,\theta^{(t)})\ \|\ P(z|x,\theta^{(t)}) \right )\leqslant 0 \end{aligned}(或者不从KL散度的角度,也可以用E[\log x ] \leqslant \log E[x] 证明),因此H(\theta^{(t+1)},\theta^{(t)})\leqslant H(\theta^{(t)},\theta^{(t)})得证;

10.3 ELBO+KL散度导出公式

  • 已知\log P(x|\theta)=\log \frac{P(x,z|\theta)}{P(z|x,\theta)}=\log P(x,z|\theta)-\log P(z|x,\theta),设隐变量 z 的先验概率分布为 q(z) ,则:

\begin{aligned} \log P(x|\theta)&=\log P(x,z|\theta)-\log P(z|x,\theta)\\&=\log P(x,z|\theta)-\log q(z)-\left [\log P(z|x,\theta) -\log q(z) \right ] \\&=\log \frac{P(x,z|\theta)}{q(z)}-\log \frac{P(z|x,\theta) }{q(z)} \end{aligned}

  • 将上式两边对 z 积分,左边=\int _z \log P(x|\theta) \cdot q(z) dz=\log P(x|\theta)
  • 右边\begin{aligned} =\int_z \log \frac{P(x,z|\theta)}{q(z)} \cdot q(z) dz-\int_z \log \frac{P(z|x,\theta)}{q(z)} \cdot q(z) dz=ELBO-KL\left [\ q(z)\ ||\ P(z|x,\theta) \ \right ] \end{aligned},ELBO为证据下界evidence lower bound;
  • 当且仅当 q(z)==P(z|x,\theta) 时,KL\left [\ q(z)\ ||\ P(z|x,\theta) \ \right ]=0\log P(x|\theta)=ELBO,此时优化 ELBO 即为优化 \log P(x|\theta),但这种方法并不保证可以找到 \log P(x|\theta)的全局最大值;

10.4 Jensen不等式导出公式

  • 根据Jensen不等式,f(E(x))\geqslant E[f(x)]\begin{aligned} \log P(x|\theta)&=\log \left [\int _z P(x,z|\theta) dz \right ] =\log \left [\int _z \frac{P(x,z|\theta) }{q(z)}\cdot q(z) dz \right ] \\ &=\log E_{q(z)}\left [ \frac{P(x,z|\theta)}{q(z)} \right ]\geqslant E_{q(z)}\left [ \log\frac{P(x,z|\theta)}{q(z)} \right ] \end{aligned}当且仅当\begin{aligned} \frac{P(x,z|\theta)}{q(z)}=C \end{aligned}时等号成立;
  • \begin{aligned}\frac{P(x,z|\theta)}{C}=q(z) \end{aligned},两边对 z 积分,得\frac{1}{C}\int_z P(x,z|\theta) dz=\frac{1}{C}P(x|\theta)=\int_z q(z) dz=1,即P(x|\theta)=C,因此\begin{aligned} q(z)=\frac{P(x,z|\theta)}{P(x|\theta)}=P(z|x,\theta) \end{aligned},与ELBO+KL散度得到相同结论; 

10.5 广义EM

  • 根据上面的推导,EM算法的目标函数\log P(X|\theta)等效为ELBO:E_{q(z)}[\log \frac{P(X,Z|\theta)}{q(z)}]=E_{q(z)}[\log P(X,Z|\theta)]-\int_z\log q(z)dz,当且仅当q(z)==P(z|x,\theta)时,ELBO达到最大值ELBO=\log P(x|\theta)
  • 但在某些场景中,后验概率P(z|x,\theta)无法直接求解,因此可以对q(z)==P(z|x,\theta)进行同步迭代,即为广义EM算法:
    • E步:q^{(t+1)}(z)=\underset{q}{argmax}\ ELBO
    • M步:\theta^{(t+1)}=\underset{\theta}{argmax}\ ELBO

11. 高斯混合模型(Gaussian Mixture Module)

11.1. 模型介绍

  • 高斯混合模型即为将某个变量的分布视为多个高斯分布叠加而成;

  •  每个样本由各个高斯分布生成的隐变量加权得来,即对每个样本 x_i ,都有一组参数(c_1, c_2, \cdots , c_K),使得p(x_i)=c_1\cdot N(\mu_1,\sigma_1^2)+c_2\cdot N(\mu_2,\sigma_2^2)+\cdots c_K\cdot N(\mu_K,\sigma_K^2)
  • 注:分布来自于各个高斯分布,但其实某个样本真正的来源只有一个,举例:先抛硬币A,若正面朝上则抛B,反之抛C,记录最终抛的结果,每个样本得到一个记录值,可能来自于B也可能来自于C,B的概率分布为p_B=N(\mu_1,\sigma_1^2),C的概率分布为p_C=N(\mu_2,\sigma_2^2),A朝上的概率为c_1,朝下的概率为c_2,此时最终记录值的概率分布为p(x_i=1)=c_1\cdot p_B(x_i=1)+c_2\cdot p_C(x_i=1),是B和C分别的概率加权叠加,但记录值(或者叫观测值)的真正来源只有1个,B或者C;

11.2. 极大似然估计

  • 首先明白,极大似然估计无法求解GMM,这里介绍其中的困难,引出EM算法求解;
  • \begin{aligned} \hat{\theta}_{MLE}&=\underset{\theta}{argmax} \log P(X) \\ &=\underset{\theta}{argmax} \log \prod _{i=1}^N p(x_i)\\ &=\underset{\theta}{argmax} \sum_{i=1}^N\log p(x_i) \\ &=\underset{\theta}{argmax} \sum_{i=1}^N\log \sum_{j=1}^K[p_j\cdot p(x_i|\mu_j,\sigma_j^2)] \\ \end{aligned}
  • 因为log中有连加,且隐变量的概率分布未知,所以这种目标函数无法用极大似然估计求出解析解,需要用迭代逼近的方法计算,例如基于梯度的方法和EM算法;

11.3. EM求解

  • E步:
    • 首先设隐变量\gamma_{j,k}=1 \ or \ 0,表示某个观测样本x_j来自某个高斯模型N(\mu_k,\sigma_k^2),每个样本的该变量可以表示为一个K维独热向量;
    • 列出似然函数\begin{aligned} P(Y,\Gamma|\theta)&=\prod _{j=1}^N \prod_{k=1}^K \left [ p_k \cdot p_N(x_j|\mu_k, \sigma_k) \right ]^{\gamma_{j,k}}\\ &=\prod _{k=1}^K \prod_{j=1}^N \left [ p_k \cdot\frac{1}{\sqrt{2\pi} \sigma_k}\exp\{-\frac{(x_j-\mu_k)^2}{2\sigma_k^2}\} \right ]^{\gamma_{j,k}}\\ \end{aligned}
    • 对数似然函数为\begin{aligned} \log P(Y,\Gamma|\theta)=\sum_{k=1}^K \sum_{j=1}^N \gamma_{j,k}\left [\log p_k -\frac{1}{2}\log 2\pi-\log\sigma_k-\frac{1}{2\sigma_k^2}(x_j-\mu_k)^2 \right ] \\ \end{aligned}
    • 对其求期望得目标函数:\begin{aligned} E_{\Gamma|X,\theta}\left [\log P(X,\Gamma|\theta) \right ]=\sum_{k=1}^K \sum_{j=1}^N E_{\gamma_{j,k}|X,\theta}\left [\gamma_{j,k} \right ]\left [\log p_k -\frac{1}{2}\log 2\pi-\log\sigma_k-\frac{1}{2\sigma_k^2}(x_j-\mu_k)^2 \right ] \\ \end{aligned}
    • \begin{aligned} E_{\gamma_{j,k}|X,\theta}\left [ \gamma_{j,k} \right ]=\frac{p_k\cdot p_N(x_j|\mu_k,\sigma_k)}{\sum_{m=1}^Kp_m\cdot p_N(x_j|\mu_m,\sigma_m)} \end{aligned},设n_k=\sum_{j=1}^N\gamma_{j,k}=\sum_{j=1}^NE _{\gamma_{j,k}|X,\theta}\left [ \gamma_{j,k} \right ]\sum_{k=1}^Kn_k=1,则目标函数的期望\begin{aligned} E_{\Gamma|X,\theta}\left [\log P(X,\Gamma|\theta) \right ]=\sum_{k=1}^K \left \{n_k\cdot \log p_k+\sum_{j=1}^N E_{\gamma_{j,k}|X,\theta}\left [\gamma_{j,k} \right ]\left [ -\frac{1}{2}\log 2\pi-\log\sigma_k-\frac{1}{2\sigma_k^2}(x_j-\mu_k)^2 \right ] \right \}\\ \end{aligned}
  • M步:
    • 对期望函数的各个参数求偏导并使其为0,得到结果(期望简化为E):
      • \begin{aligned} \hat{\mu}_k=\frac{\sum_{j=1}^NE\left[ \gamma_{j,k} \right ]x_j}{\sum_{j=1}^NE\left[ \gamma_{j,k}\right ]} \end{aligned}
      • \begin{aligned} \hat{\sigma}_k^2=\frac{\sum_{j=1}^NE\left[ \gamma_{j,k} \right ]\left (x_j -\mu_k \right )^2}{\sum_{j=1}^N E\left[ \gamma_{j,k}\right ]} \end{aligned}
      • \begin{aligned} \hat{p}_k=\frac{n_k}{N} \end{aligned}

12. 变分推断(Variational Inference)

12.1. 背景回顾

12.1.1. 两种思路,两种问题

  • 概率派:建立概率模型,定义目标函数/损失函数,将问题视为优化问题,得到解析解/数值解(无法得到解析解时,使用梯度下降等方法求得数值解),例如:线性回归、SVM、EM等算法;
  • 贝叶斯派:建立概率图模型,求出后验概率p(\Theta|X),并根据后验概率进行决策p(\tilde{x}|X)=\int _{\Theta}p(\tilde{x}|\Theta)p(\Theta|X)d\Theta,将问题视为积分问题,后验概率有时可以精确推断,有时只能进行近似推断,而近似推断的确定性推断即为变分推断,近似推断还有一种随机推断则是MCMC方法;

12.1.2. 推断与生成

  • 在概率派思想中,遇到无法观测的隐变量时,模型参数无法通过极大似然估计等方法求出解析解,此时通常使用EM等迭代逼近的方法求出数值解;
  • 在贝叶斯派思想中,模型参数 \theta 遵循某种概率分布(先验),与隐变量 Z,观测变量 X 共同构成概率图模型,相互之间存在相关性和条件独立性,而这个模型中有两种问题:推断与生成:
    • 推断:由观测变量的样本集合估计隐变量的概率分布,与EM算法的思路相同,通过最大化 \log P(X) 来得到 Z 的分布,即求解P(Z\ | \ X)
    • 生成:估计隐变量与观测变量之间的联系,并利用大量随机的隐变量生成新的观测变量(VAE),即求解P(X\ | \ Z)

12.2. 公式推导

  • 根据贝叶斯派思想,模型参数 \theta 也满足某种概率分布,而变分推断是要对参数和隐变量都进行近似推断,因此设 X 为观测数据样本Z 为待推断变量,包含隐变量与模型参数(与EM略有不同);
  • 同EM中的推导(这里再写一遍,\log P(X)=\log P(X,Z)-\log P(Z|X)=\log\frac{P(X,Z)}{q(Z)}-\log \frac{P(Z|X)}{q(Z)},左右对Z求积分得\begin{aligned} \log P(X)=\int_Z q(Z)\cdot \log \frac{P(X,Z)}{q(Z)}\ dZ-\int_Zq(Z)\cdot \log \frac{P(Z|X)}{q(Z)}\ dZ=ELBO+KL[q(Z)\ ||\ P(Z|X)]\end{aligned}),ELBO即为变分下界,\begin{aligned}ELBO=\int_Zq(Z)\log P(X,Z) dZ-\int_Z q(Z) \log q(Z)dZ=J(q)+H(Z) \end{aligned}
  • 引入平均场理论(mean field theory),设 Z 中存在M个部分,每两个部分两两独立,则有q(Z)=\prod_{i=1}^M q^{(i)}(Z^{(i)}),则\begin{aligned} J(Z)=\int_Z \prod_{i=1}^M q^{(i)}(Z^{(i)})\log P(X,Z)\ dZ^{(1)}dZ^{(2)}\cdots dZ^{(M)} \end{aligned},此时单独求解某个部分的分布q_j(Z_j),则可以固定其他部分,即\begin{aligned} J(Z)&=\int_{Z^{(j)}} q^{(j)}(Z^{(j)}) \cdot \left [ \int_{Z^{(i\neq j)}}\prod_{i\neq j}^M q^{(i)}(Z^{(i)})\log P(X,Z)\ dZ^{(1)}dZ^{(2)}\cdots dZ^{(M)} \right ] dZ^{(j)}\\ &=\int_{Z^{(j)}} q^{(j)}(Z^{(j)}) \cdot E_{\prod_{i\neq j}^Mq^{(j)}(Z^{(j)})}\left [ \log P(X,Z) \right ] dZ^{(j)} \end{aligned}
  • 同时\begin{aligned} H(Z)&=-\int_Zq(Z)\log q(Z) dZ=-\int_{Z^{(1)}\cdot Z^{(2)}\cdots Z^{(M)}}\prod _{i=1}^M q^{(i)}(Z^{(i)})\cdot \left [\sum_{i=1}^M \log q^{(i)}(Z^{(i)}) \right ] \ dZ^{(1)}dZ^{(2)}\cdots dZ^{(M)} \\ &=\sum_{k=1}^M\left [-\int_{Z^{(1)}\cdot Z^{(2)}\cdots Z^{(M)}}\prod _{i=1}^M q^{(i)}(Z^{(i)})\cdot \log q^{(k)}(Z^{(k)})\ dZ^{(1)}dZ^{(2)}\cdots dZ^{(M)} \right ]\\ &=\sum_{k=1}^M\left \{-\int_{Z^{(k)}}q^{(k)}(Z^{(k)})\cdot \log q^{(k)}(Z^{(k)})\ dZ^{(k)}\left [\prod _{i\neq k}^M\int_{Z^{(i)}} q^{(i)}(Z^{(i)})dZ^{(i)} \right ] \right \}\\ &=-\sum_{k=1}^M\int_{Z^{(k)}}q^{(k)}(Z^{(k)})\cdot \log q^{(k)}(Z^{(k)})\ dZ^{(k)} \end{aligned},当其他部分固定,单独求解q^{(j)}(Z^{(j)})时,\begin{aligned} H(Z)&=-\int_{Z^{(j)}}q^{(j)}(Z^{(j)})\cdot \log q^{(j)}(Z^{(j)})\ dZ^{(j)} +C\end{aligned}
  • 因为\begin{aligned} E_{\prod_{i\neq j}^Mq^{(i)}(Z^{(i)})}\left [ \log P(X,Z) \right ] \end{aligned}积分的结果是一个关于(X,Z_j)的对数概率密度函数,因此可以设\begin{aligned} \log \hat{P}(X,Z^{(j)})= E_{\prod_{i\neq j}^Mq^{(i)}(Z^{(i)})}\left [ \log P(X,Z) \right ] \end{aligned},则变分下界为\begin{aligned} ELBO=L(Z)+H(Z)&=\int_{Z^{(j)}}q^{(j)}(Z^{(j)}) \log \hat{P}(X,Z^{(j)})dZ^{(j)}-\int_{Z^{(j)}}q^{(j)}(Z^{(j)}) \log q^{(j)}(Z^{(j)}) dZ^{(j)} +C\\ &=\int_{Z^{(j)}}q^{(j)}(Z^{(j)})\frac{\log \hat{P}(X,Z^{(j)})}{q^{(j)}(Z^{(j)})} dZ^{(j)}+C\\ &=-KL\left [ q^{(j)}\ ||\ \hat{P}(X,Z^{(j)}) \right ]+C \end{aligned},当且仅当\begin{aligned} q_j == \hat{P}(X,Z^{(j)}) \end{aligned}时,ELBO达到最大;
  • 如此反复迭代求解,以坐标上升法求出所有参数和隐变量;

12.3. 随机梯度上升变分推断SGVI

12.3.1. 公式推导

  • 用另一种思路求解变分推断的隐变量,设 \theta 是模型参数,隐变量 Z 遵从某种确定的分布,例如指数族分布,则求解 Z 即为求解分布中的参数 \phi,即\hat{\phi}=\underset{\phi}{\arg\max}\ ELBO
  • 目标函数无解析解,只能使用迭代逼近的方法求数值解,梯度上升法是一个合适的算法,对目标函数求梯度,(积分和求导符号可以互换位置)\begin{aligned} \bigtriangledown_{\phi} ELBO&=\bigtriangledown_{\phi}\ E_{q_{\phi}}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] \\ &=\bigtriangledown_{\phi}\ \int q_{\phi}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ]dZ \\ &=\int\ \bigtriangledown_{\phi}\left\{ q_{\phi}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] \right \} dZ \\ &=\int\ \bigtriangledown_{\phi} q_{\phi}\cdot\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] dZ+\int\ q_{\phi} \cdot \bigtriangledown_{\phi}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] dZ \\ \end{aligned}
  • 右边第二项中\begin{aligned}\log P_{\theta}(X,Z) \end{aligned}\phi 无关,梯度为0,因此\begin{aligned} \int\ q_{\phi} \cdot \bigtriangledown_{\phi}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] dZ &= \int\ q_{\phi} \cdot \bigtriangledown_{\phi}\left[-\log q_{\phi} \right ] dZ \\ &=- \int\ \bigtriangledown_{\phi}q_{\phi} dZ\\ &= -\bigtriangledown_{\phi}\int\ q_{\phi} dZ\\ &= -\bigtriangledown_{\phi}\ 1=0 \end{aligned}
  • 因此\begin{aligned} \bigtriangledown_{\phi} ELBO&=\int\ \bigtriangledown_{\phi} q_{\phi}\cdot\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] dZ \\ &=\int\ q_{\phi} \cdot \bigtriangledown_{\phi} \log q_{\phi} \cdot\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] dZ \\ &=E_{q_{\phi}} \left \{\bigtriangledown_{\phi} \log q_{\phi} \cdot\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] \right \}\\ \end{aligned},(其中\begin{aligned} q_{\phi} \cdot \bigtriangledown_{\phi}\log q_{\phi} =\bigtriangledown_{\phi}q_{\phi} \\ \end{aligned}由上一段推出),由此将目标函数梯度的求解变为一个求期望问题,所以可以通过蒙特卡洛方法求出期望的近似解,从而进行梯度上升迭代计算;

12.3.2. SGVI简化

  • 对于SGVI中的梯度求解,常常用到蒙特卡洛方法求期望,但采样过程中有一定概率得到一个接近0的q_{\phi}值,则其中的\log q_{\phi}会变成一个很小很小的负数,导致采样结果的方差很大,本身梯度上升法就是一种迭代计算的数值解,过程中存在一定偏差,如果采样过程的误差过大则会使得整个算法收敛缓慢,甚至无法得到最优解;
  • 为了解决上述问题,也为了简化求梯度的过程,引入一个重参数化机制
    • 原本模型中 X 是观测变量, Z 是隐变量, Z 的分布有参数\phi,与 X 无关;
    • 重新定义 Z ,设一个新的随机变量 \varepsilon ,遵循某种已知分布p(\varepsilon ),则 z=g_{\phi}(\varepsilon ,X),此时q_{\phi}(Z|X)dZ=p(\varepsilon )d\varepsilonE_{q_{\phi}}(*)=E_{p_{\varepsilon }}(*)
  • 则重新推导梯度式:\begin{aligned} \bigtriangledown_{\phi} ELBO&=\bigtriangledown_{\phi}\ E_{q_{\phi}}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] \\ &=\bigtriangledown_{\phi}\ E_{p_{\varepsilon }}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ] \\ &=E_{p_{\varepsilon }}\left\{\bigtriangledown_{\phi}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ]\right\}\\ &=E_{p_{\varepsilon }}\left\{\bigtriangledown_{Z}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ]\cdot \bigtriangledown_{\phi}Z \right\}\\ &=E_{p_{\varepsilon }}\left\{\bigtriangledown_{Z}\left[\log P_{\theta}(X,Z)-\log q_{\phi} \right ]\cdot \bigtriangledown_{\phi}g_{\phi}(\varepsilon ,X) \right\}\\ \end{aligned},此时再采用蒙特卡洛方法则不会出现误差;

13. 马尔可夫链蒙特卡洛法(Markov Chain & Monte Carlo,MCMC)

13.1. 蒙特卡洛方法简介

  • 在推断中的近似推断问题中,常常会使用随机近似推断方法,例如最经典的MCMC,核心思想是蒙特卡洛方法,即基于采样的随机近似方法;
  • 例如,要计算推断中的期望E_{z}\left[f(z) \right ],假设随机变量z是连续/无限取值的离散,则对其求期望常常是困难的(连续-积分/无限离散-无穷级数),此时可以随机采样N个样本\left \{ z_1,z_2,\cdots,z_N \right \},则期望的估计值为E_{z}\left[f(z) \right ]=mean\left \{ f(z_1),f(z_2),\cdots,f(z_N) \right \}
  • 蒙特卡洛方法中常常使用以下采样方法:
    • 概率分布采样:利用变量的概率分布p(z)进行随机数生成,但这个方法仅限概率分布简单的情况;
    • Adjection Sampling (接收-拒绝采样法):
      • 当随机变量的概率分布p(z)较为复杂时,可以设一个简单分布q(z),以及寻找一个常数M,使得定义域上的任意\forall \ z_i,\ \ M\cdot q(z)\geqslant p(z)
      • 引入一个新的参数\alpha(z)=\frac{p(z)}{M\cdot q(z)},称为接收率,其值在0-1之间,且越接近1表示p(z)M\cdot q(z)越接近;
      • 此时可以从简单分布q(z)中进行采样,每采样一个z_i的同时也从均匀分布U(0, 1)中采样一个u_i,当u_i\leqslant \alpha(z_i)时,接受采样结果,否则拒绝,需要重新采样;
      • 这种采样方法又叫接受-拒绝法,其优点是容易实现,缺点是效率低,如果p(z)所含面积占M\cdot q(z)所含面积比例过小,则采样很容易被拒绝,需要重新采样,类似的,如果待采样变量维度过高,也很容易采样被拒绝;
    • 重要度采样:同样设一个简单分布q(z),则原本难算的期望E_z\left[f(x) \right ]=\int_zp(z)f(z)dz=\int_z\frac{p(z)}{q(z)}f(z)\cdot q(z)dz=E_{q(z)}\left[\frac{p(z)}{q(z)}f(z) \right ]可以化为另外一个期望,使得原本从p(z)中采样的困难转为从简单分布q(z)中采样,而\frac{p(z)}{q(z)}称为重要度;

13.2. 马尔科夫链简介

  • 马尔可夫链p(x_t)=p(x|x_1, x_2, \cdots, x_{t-1}),其中x_t表示序列上某个步/时间的样本,若一阶马尔可夫,则p(x_t)=p(x|x_{t-1}),序列上所有随机变量的取值集合相同,被称为状态空间 Sp(x_t)=p(x|x_{t-1})被称为 t 时刻的状态转移函数
  • 这里只讨论一阶马尔可夫链,因为更高阶的马尔可夫链可以推导转换为一阶马尔科夫链,同时若任意的 s,有P(x_{t+s}|x_{t+s-1})=P(x_t|x_{t-1}),则称为时间齐次的马尔可夫链,这里也只讨论时间齐次的马尔科夫链;
  • 对于状态空间 S 全部是离散值的马尔可夫链,其所有的状态转移函数可以写作一个矩阵,即\left \{ \ p_{ij}=P(X_t=i|X_{t-1}=j) \ \right \},称为状态转移矩阵,满足\sum_{i}p_{ij}=1
  • 基于状态转移函数,马尔可夫链可以表示为,对 t 时刻,取值为 i 的概率为P(x_t=i)=\sum_{j}p_{ij}\cdot P(x_{t-1}=j)=\sum_{j}p_{ij}\cdot \sum_kp_{jk} P(x_{t-2}=k)=\cdots,为了方便表示,设矩阵\Pi_{S\times T}=\{ \pi_i(t) \}=P(x_t=i),即为将每一步的计算值存储下来便于后面的计算,其中\vec{\pi}(0) =\left [ \pi_{S_0}(0), \pi_{S_1}(0), \cdots \right ]^T被称为初始状态分布,常常是个独热向量,表示马尔可夫链从某个固定状态开始;
  • 如果马尔可夫链的状态空间上存在某个分布\vec{\pi}=\left [ \pi_{S_0}, \pi_{S_1}, \cdots \right ]^T,使得\pi=P\pi,则称该分布为平稳分布,如果马尔可夫链进入平稳分布,则其接下来的状态分布较为平稳,变化很小,同时一个马尔可夫链可能存在唯一的平稳分布,无穷多个平稳分步,或不存在平稳分布;
  • 马尔可夫链的性质:
    • 可约性:若对于任意初始状态 x_0=i 和任意状态 j ,都至少存在一个时刻 t>0 ,使得P(x_t=j|x_0=i)>0,则说明任意初始状态经过一定时间可以到达任意状态,这种情况称马尔可夫链式不可约的,否则,说明从某个状态出发会永远到不了其他某些状态,这个马尔可夫链称为可约的;
    • 周期性:对任意初始状态 x_0=i,若其经过一定时间后可以回到状态 x_t=i,且所有最短途径时长\left \{ t_i \right \}的最大公约数不为1,则称该马尔科夫链具有周期性,否则是非周期的;
    • 正常返:定义 p_{ij}^t 表示P(x_t=i,x_{t-m}\neq i|x_0=j),即从状态 j 出发,经过 t 时间首次来到状态 i 的概率,若对任意 j 和 i ,都有\lim_{t \to \infty }p_{ij}^t>0,则称为正常返的马尔科夫链,否则称为不是正常返;
    • 可逆性:对于一个马尔可夫链,若存在一种状态分布\vec{\pi}=[\pi_0, \pi_1, \cdots],使得任意两个状态i,j都有P(x_t=i|x_{t-1}=j)\cdot P(x_{t-1}=j)=P(x_t=j|x_{t-1}=i)\cdot P(x_{t-1}=i),即p_{ij}\pi_j=p_{ji}\pi_i(该式被称作细致平衡方程,detailed balance),则称该马尔可夫链是可逆的;
  • 马尔可夫链性质的定理:
    • 不可约且非周期的有限状态马尔可夫链一定有唯一平稳分布存在;
    • 不可约且非周期,正常返的马尔可夫链一定有唯一平稳分布存在;
    • 遍历定理:对于一个不可约,非周期,正常返的马尔可夫链,设其唯一的平稳分布为\vec{\pi},则对于任意初始状态x_0=j,有\lim_{t \to \infty } P(x_t=i|x_0=j)=\pi_i,也就是无论初始状态,未来经过无限时刻,状态分布都会趋近于平稳分布,同时对于一个定义在状态空间上的函数f(x),x\in S,则lim_{t \to \infty}\frac{1}{t}\sum_{k=1}^tf(x_k) \to E_{\vec{\pi}}\ [f(x)]=\sum _if(S_i)\pi_i
    • 满足细致平稳方程的分布就是该马尔可夫链的平稳分布;

13.3. MCMC算法

  • MCMC适用于随机变量高维且相互不独立,密度函数复杂的情况,在这样的分布中求某个函数的期望E_{p(x)}\left [f(x) \right ],主要思想是在随机变量的状态空间(所有取值)上定义一个满足便利定理的马尔可夫链,则其平稳分布下函数值的加权平均即为函数的数学期望;
  • 因此MCMC计算期望的方法为\hat{E}[f]=\frac{1}{n-m}\sum_{i=m+1}^{n}f(x_i),其中 m 是一个较大的数,使得t>m时的 x_t 分布趋于平稳分布,在m之前的遍历称为燃烧期;
  • 因为马尔科夫链满足遍历定理,所以随机游走的起点不影响最终结果,因为都会收敛到平稳分布,这个方法比接收-拒绝采样法效率高;

13.3.1. Metropolis-Hastings算法

  • 在MCMC算法中,第一个难点是定义一个马尔可夫链用于随机游走,使其满足遍历定理,无论从任何起点出发都将收敛到平稳分布,对于离散分布,要定义状态转移矩阵\left \{ p_{ij}=P(x_t=i|x_{t-1}=j) \right \},对于连续分布,要定义状态转移核函数p(i,j)=P(x_t=i|x_{t-1}=j)
  • 在MH算法中,可以定义一个建议分布 q(i,j) 较为简单易于操作,以及一个接收分布 \alpha (i,j)=\min \{ 1,\frac{p(j)q(i,j)}{p(i)q(j,i)} \} ,则转移核p(i,j)=q(i,j)\alpha (i,j)=\left\{\begin{matrix} q(i,j) & p(j)q(i,j)\geqslant p(i)q(j,i)\\ q(j,i)\cdot \frac{p(i)}{p(j)} & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.
  • 通俗理解,MH算法的流程为,定义一个简易分布q(i,j),作为转移函数进行随机游走,每次从状态 x_{t-1}=j 转移到状态 i 时,计算接受分布\alpha (i,j)=\min \{ 1,\frac{p(j)q(i,j)}{p(i)q(j,i)} \},并在0-1之间采样随机数u,若u<\alpha (i,j)则状态转移成功,x_{t}=i,否则停留在原状态x_{t}=j
  • \begin{aligned} p(j)p(i,j)&=p(j)\left\{\begin{matrix} q(i,j) & p(j)q(i,j)\geqslant p(i)q(j,i)\\ q(j,i)\cdot \frac{p(i)}{p(j)} & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.\\ &=\left\{\begin{matrix} p(j)q(i,j) & p(j)q(i,j)\geqslant p(i)q(j,i)\\ p(i) q(j,i) & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.\\ &=p(i)\left\{\begin{matrix} q(i,j)\cdot \frac{p(j)}{p(i)} & p(j)q(i,j)\geqslant p(i)q(j,i)\\ q(j,i) & p(j)q(i,j)< p(i)q(j,i) \end{matrix}\right.\\ &=p(i)p(j,i) \end{aligned},满足细致平衡方程,因此上面所定义的转移核函数是可逆的,即存在某个分布p(x)使得p(i)p(j,i)=p(j)p(i,j),同时该分布为平稳分布;
  • Metropolis选择:将建议分布设为可逆分布,即q(i,j)=q(j,i),此时接收分布简化为\alpha(i,j)=\min\{ 1,\frac{p(j)}{p(i)} \},这是MH算法最初采用的方法,将q(i,j)定义为一个以 j 为均值,方差为常数矩阵的正态分布,或是其他在 j 附近概率较高的简单分布;
  • 单分量Metropolis-Hastings算法:当多元变量分布抽样困难时,可以逐步对变量的每一个维度的分量进行抽样,此时建议分布变为q(x_t^{(m)}|x_t^{(-m)},x_{t-1}^{(m)}),其中x_t^{(m)}为t时刻状态的第m个分量,x_t^{(-m)}=\{ x_t^{(1)}, x_t^{(2)}, \cdots, x_t^{(m-1)}, x_{t-1}^{(m+1)}, \cdots, x_{t-1}^{(K)} \},类似于每次固定K-1个维度,更新一个维度的过程,同时接受分布\alpha (x_t^{(m)}|x_{t-1}^{(m)}, x_t^{(-m)})=\min\{ 1,\frac{p(x_{t}^{(m)}|x_{t}^{(-m)})q(x_{t-1}^{(m)}|x_{t}^{(m)}, x_{t}^{(-m)})}{p(x_{t-1}^{(m)}|x_{t}^{(-m)})q(x_t^{(m)}|x_{t-1}^{(m)}, x_{t}^{(-m)})} \}

13.3.2. 吉布斯抽样

  • 满条件分布:对于一个多元概率联合分布p(x),x=[x^{(1)}, x^{(2)}, \cdots, x^{(K)}]^T,定义一个条件概率分布p(x^I|x^{-I}),其中I是变量维度的一个子集,即x^I=\{ x^{(m)},m\in I \}, x^{-I}=\{ x^{(m)},m\notin I \},也就是所有维度都在括号里出现,无论在 | 之前还是之后,这个条件概率分布叫做满条件分布,满条件分布满足性质:\frac{p(x^I_{t-1}|(x^{-I}_{t-1})}{p(x^I_{t}|(x^{-I}_{t})}=\frac{p(x_{t-1})}{p(x_{t})},左边计算大大方便,因此可以用满条件分布简化一些计算;
  • 吉布斯抽样的方法类似于单分量MH算法,不同的是吉布斯抽样使用 x_{t-1} 的满条件分布代替建议分布,由于满条件分布比例于原分布的性质,使用满条件分布可以完全模拟原分布p(x),因此无需MH算法中的接受分布,或者说接受率为1;

14. 隐马尔可夫模型(Hidden Markov Model,HMM)

14.1. 背景

  • 概率图模型是贝叶斯派思想的成果,在概率图模型中加入时间的维度,就成为了动态系统(Dynamic Model),其中每个样本是一个时间序列,序列上每个点有观测变量和隐变量,按照隐变量取值类型分为:变量取值离散——HMM;取值连续线性——卡曼滤波(Kalman Filter);取值连续非线性——例子滤波(Particle Filter);
  • 贝叶斯派思想的主要问题有 Learning 和 Inference,这些问题也同样是HMM中所关注的,其中Learning是通过样本数据(通常是观测变量X)求模型参数,而Inference又会分为不同形式:
    • Decoding:P(Z|X),通过所有观测变量序列推断所有隐变量序列;
    • Probability of Evidence: P(X|\theta),通过模型参数推断所有观测变量序列;
    • Filtering:P(z_t|x_1,\cdots,x_t),通过历史观测信息推断当前时刻隐状态(实时推断,常见于Online计算);
    • Smoothing:P(z_t|x_1,\cdots,x_T),通过所有观测信息推断当前时刻隐状态(复盘式推断,常见于Offline计算);
    • Prediction:P(z_{t+1}|x_1,\cdots,x_t),通过历史观测信息推断后续时刻隐状态,也可以预测观测变量;

14.2. 模型详解

14.2.1. 模型定义

  • 先定义观测变量和状态变量的符号:HMM中隐变量取值是离散的,因此可以设观测变量序列为 O=\{ o_1,o_2,\cdots,o_T \},o_t\in V=\{ v_1,v_2,\cdots,v_M \},隐变量序列为   I=\{ i_1,i_2,\cdots,i_T \},i_t\in Q=\{ q_1,q_2,\cdots,q_N \}
  • HMM常用的定义方式为 \lambda=(\pi, A, B),其中 \pi 为初始状态,A=\{ a_{ij} \} 为状态转移矩阵,矩阵中元素a_{ij}=P(i_{t+1}=q_j|i_t=q_i)B=\{ b_{j}(o_t) \} 为发射矩阵,矩阵中元素b_{j}(o_t)=P(o_t=v_k|i_t=q_j)
  • HMM满足两个基本假设:
    • n阶齐次马尔可夫性:某个时刻隐状态只与前n个时刻有关,即P(i_{t+1}|i_1,i_2,\cdots,i_t,o_1,o_2,\cdots,o_t)=P(i_{t+1}|i_t)
    • 观察独立性:某个时刻观测变量只与当前时刻隐变量有关,即P(o_{t}|i_1,i_2,\cdots,i_t,o_1,o_2,\cdots,o_{t-1})=P(o_{t}|i_{t})
  • HMM的三个基本问题:
    • Evaluation:求 P(O|\lambda)
    • Learning:求 \lambda=\underset{\lambda}{argmax}P(O|\lambda)
    • Decoding:求 \hat{I}=\underset{I}{argmax}P(I|O,\lambda)

14.2.2. Evaluation

  • 前向算法:
    • 正常求解P(O|\lambda),可以对其进行分解:\begin{aligned} P(O|\lambda)&=\sum_I P(I,O|\lambda)=\sum_I P(O|I,\lambda)\cdot P(I|\lambda)\\ \end{aligned},其中\begin{aligned} P(I|\lambda)&=P(i_1,i_2,\cdots,i_T|\lambda)=P(i_T|i_1,i_2,\cdots,i_{T-1},\lambda)\cdot P(i_1,i_2,\cdots,i_{T-1}|\lambda)\\ &=P(i_T|i_{T-1},\lambda)\cdot P(i_1,i_2,\cdots,i_{T-1}|\lambda)=a_{i_{T-1},i_T}\cdot P(i_1,i_2,\cdots,i_{T-1}|\lambda)\\ &=\cdots=\pi (i_1) \cdot \prod _{t=2}^Ta_{i_{t-1},i_t} \end{aligned}\begin{aligned} P(O|I,\lambda)=\prod_{t=1}^Tb_{i_t}(o_t) \end{aligned},因此带入原式可得:\begin{aligned} P(O|\lambda)&=\sum_I P(O|I,\lambda)\cdot P(I|\lambda)\\ &=\sum_I\left [\pi (i_1) \cdot \prod _{t=2}^Ta_{i_{t-1},i_t}\cdot \prod_{t=1}^Tb_{i_t}(o_t) \right ]\end{aligned},但其中\begin{aligned} \sum_I=\sum_{t_1}\sum_{t_2}\cdots \sum_{t_T} \end{aligned}的遍历计算量较大,因此需要考虑使用递推的思想简化,即为前向传播算法;
    • \begin{aligned} \alpha_t(q_i)=P(o_1,\cdots,o_t,i_t=q_i|\lambda),\alpha_T(q_i)=P(O,i_T=q_i|\lambda) \end{aligned},只要能求出\begin{aligned} \alpha_T(q_i) \end{aligned},即有\begin{aligned} P(O|\lambda)=\sum_{i=1}^NP(O,i_t=q_i|\lambda)=\sum_{i=1}^N\alpha_T(q_i) \end{aligned}
    • \begin{aligned} \alpha_T(q_i) \end{aligned}的引入需要简化计算,现在推导 \begin{aligned} \alpha_{t+1}(q_i) \end{aligned} 与 \begin{aligned} \alpha_{t}(q_i) \end{aligned} 的递推式:\begin{aligned} \alpha_{t+1}(q_j)&=P(o_1,\cdots,o_{t+1},i_{t+1}=q_j|\lambda)\\ &=\sum_{i=1}^NP(o_1,\cdots,o_{t+1},i_{t+1}=q_j,i_t=q_i|\lambda)\\ &=\sum_{i=1}^NP(o_{t+1}|o_1,\cdots,o_t,i_{t+1}=q_j,i_t=q_i|\lambda)\cdot P(o_1,\cdots,o_t,i_{t+1}=q_j,i_t=q_i|\lambda)\\ &=\sum_{i=1}^NP(o_{t+1}|i_{t+1}=q_j,\lambda)\cdot P(i_{t+1}=q_j|o_1,\cdots,o_t,i_t=q_i,\lambda)\cdot P(o_1,\cdots,o_t,i_t=q_i|\lambda)\\ &=\sum_{i=1}^NP(o_{t+1}|i_{t+1}=q_j,\lambda)\cdot P(i_{t+1}=q_j|i_t=q_i,\lambda)\cdot P(o_1,\cdots,o_t,i_t=q_i|\lambda)\\ &=\sum_{i=1}^N b_j(o_{t+1})\cdot a_{ij}\cdot\alpha_t(q_i) \end{aligned},由此即可从 \alpha_1(i_1) 一路计算到 \alpha_T(i_T) ,从而完成前向算法;
  • 后向算法:
    • 同样是为了求解P(O|\lambda)而提出,假设已有信息是T时刻的变量,从后往前递推的一种算法;
    • \beta_t(q_j)=P(o_{t+1},\cdots, o_{T}|i_t=q_j,\lambda),\beta_1(q_j)=P(o_{2},\cdots, o_{T}|i_1=q_j,\lambda),则有\begin{aligned} P(O|\lambda)&=P(o_1,\cdots,o_T|\lambda)\\ &=\sum_{j=1}^NP(o_1,\cdots,o_T,i_1=q_j)\\ &=\sum_{j=1}^NP(o_1,\cdots,o_T|i_1=q_j)\cdot\underset{\pi_j}{ \underbrace{P(i_1=q_j)}}\\ &=\sum_{j=1}^N\underset{=P(o_1|i_1=q_j)=b_j(o_1)}{\underbrace{P(o_1|o_2,\cdots,o_T,i_1=q_j)}}\cdot P(o_2,\cdots,o_T|i_1=q_j)\cdot\pi_j\\ &=\sum_{j=1}^Nb_j(o_1)\cdot \beta_1(q_j)\cdot\pi_j \end{aligned}
    • 同样,需要推导递推式:\begin{aligned} \beta_t(q_j)&=P(o_{t+1},\cdots, o_{T}|i_t=q_j,\lambda) \\ &=\sum_{k=1}^NP(o_{t+1},\cdots,o_T,i_{t+1}=q_k|i_t=q_j,\lambda)\\ &=\sum_{k=1}^N\underset{P(o_{t+1},\cdots,o_T|i_{t+1}=q_k,\lambda)}{\underbrace{P(o_{t+1},\cdots,o_T|i_{t+1}=q_k,i_t=q_j,\lambda)}}\cdot \underset{a_{jk}}{\underbrace{P(i_{t+1}=q_k|i_t=q_j,\lambda)}}\\ &=\sum_{k=1}^N\underset{=P(o_{t+1}|i_{t+1}=q_k,\lambda)=b_k(o_{t+1})}{\underbrace{P(o_{t+1}|o_{t+2},\cdots,o_T,i_{t+1}=q_k,\lambda)}}\cdot \underset{\beta_{t+1}(q_k)}{\underbrace{P(o_{t+2},\cdots,o_T|i_{t+1}=q_k,\lambda)}}\cdot a_{jk}\\ &=\sum_{k=1}^Nb_k(o_{t+1})\cdot \beta_{t+1}(q_k) \cdot a_{jk}\end{aligned}

14.2.3. Learning

  • Learning问题即为求解模型参数的过程,从EM思想入手,EM算法的公式为\begin{aligned}\theta^{(t+1)}=\underset{\theta}{\arg \max} \ E_{z|x,\theta^{(t)}}\left [\log P(x, z|\theta)\right ]=\underset{\theta}{\arg \max} \ \int_z\left [\log P(x, z|\theta)\cdot P(z|x,\theta^{(t)})\right ]dz \end{aligned},带入HMM的各个变量得到公式:\lambda^{(s+1)}=\underset{\lambda}{\arg\max}\sum_{I}\log P(O,I|\lambda)\cdot P(I|O,\lambda^{(s)})(因为HMM有时序性,为了避免和时间变量t重复,这里改用s,表示迭代的某个step),其中\lambda^{(s)}=(\pi^{(s)},A^{(s)},B^{(s)}),且P(I|O,\lambda^{(s)})=\frac{P(I,O|\lambda^{(s)})}{P(O|\lambda^{(s)})},分母的值与 \lambda 无关,因此式子还可以写作\lambda^{(s+1)}=\underset{\lambda}{\arg\max}\sum_{I}\log P(O,I|\lambda)\cdot P(O,I|\lambda^{(s)})=\underset{\lambda}{\arg\max}Q(\lambda,\lambda^{(s)})
  • 易知P(O,I|\lambda)=\pi_{i_1}\cdot \prod_{t=2}^Ta_{i_{t-1},i_t}\cdot \prod_{t=1}^Tb_{i_t}(o_t),带入公式得\begin{aligned} Q(\lambda,\lambda^{(s)})&=\sum_{I}\log P(O,I|\lambda)\cdot P(O,I|\lambda^{(s)})\\ &=\sum_I\left [ \left( \log\pi_{i_1}+\sum_{t=2}^T\log a_{i_{t-1},i_t}+\sum_{t=1}^T \log b_{i_t}(o_t) \right) \cdot P(O,I|\lambda^{(s)}) \right] \end{aligned},再分别使用拉格朗日乘子法等方法求解各个变量(推导过程略,简单繁杂,可以自己尝试):
    • \pi_{i}^{(s+1)}=\frac{P(O,i_1=q_i|\lambda^{(s)})}{P(O|\lambda^{(s)})}
    • a_{ij}^{(s+1)}=\frac{\sum_{t=1}^{T-1}P(O,i_{t}=q_i,i_{t+1}=q_j|\lambda^{(s)})}{\sum_{t=1}^{T-1}P(O,i_t=q_i|\lambda^{(s)})}
    • b_{j}^{(s+1)}(k)=\frac{\sum_{t=1}^{T}P(O,i_{t}=q_j|\lambda^{(s)})\cdot I(o_t=v_k)}{\sum_{t=1}^{T}P(O,i_{t}=q_j|\lambda^{(s)})}

14.2.4. Decoding

  • Decoding过程核心任务为通过观测变量的信息计算最有可能的隐状态序列,即 P(i_{t+1}|o_1,\cdots,o_t,\lambda) (也和预测问题一起讨论,例如可以预测P(o_{t+1}|o_1,\cdots,o_t,\lambda)); 
  • 近似算法:
    • 近似算法使用贪心算法的思想,在每一步取概率最大的隐状态,得到一个隐状态序列;
    • 给定模型 \lambda 和观测序列 O ,则时间 t 处在某个状态的概率为 \gamma_t(q_j)=\frac{\alpha_t(q_j)\cdot \beta_t(q_j)}{P(O|\lambda)}=\frac{\alpha_t(q_j)\cdot \beta_t(q_j)}{\sum_{k=1}^N \alpha_t(q_k)\cdot \beta_t(q_k)},而最有可能的状态即为i_t^*=\underset{1\leqslant j\leqslant N}{\arg\max}\ \gamma_t(j)
  • 维特比算法(Viterbi):
    • 维特比算法是采用动态规划的方法求解HMM的解码问题;
    • 首先定义目标概率值:\delta _t(j)=\underset{i_1,i_2,\cdots,i_{t-1}}{\max}P(i_1,\cdots,i_{t-1},i_t=q_j,o_1,\cdots,o_t|\lambda),也就是到达 t 时刻的目标隐状态所经过概率最大的路径的概率值(类似于最短路径等问题,这类问题非常适合使用动态规划方法)(此时求出的 \delta_t(j) 只是到达状态 q_j 的各个路径的概率中,最大的一个概率值,并不能得到具体的路径);
    • 则最大概率路径的递推式为:\delta_{t+1}(k)=\left[\underset{1\leqslant j\leqslant N}{\max} \delta_t(j)\cdot a_{jk} \right]\cdot b_k(o_{t+1})
    • 通过递推以及\delta_1(j)=\pi_j\cdot b_j(o_1),可以算得最终时刻各个状态的最大概率,则最终时刻的状态为i_T^*=\underset{1\leqslant j\leqslant N}{\arg\max}\ \delta_T(q_j)
    • 此时还需要另一个递推方法来求出完整的状态序列(路径),因此定义逆推概率\Psi _t(j)=\underset{1\leqslant j\leqslant N}{\arg\max}\ \delta_{t-1}(q_k)\cdot a_{kj},表示当 t 时刻隐状态是 q_j 时,上一个节点 i_{t-1} 的最可能的状态,即i_{t-1}^*=\Psi_{t}(i_{t}^*),依次逆推直至 i_1

14.2.5. 其他推断任务

  • Filtering:公式(简化)为 P(i_t|o_{1:t}),可以化简为:\begin{aligned} P(i_t|o_{1:t})=\frac{P(i_t,o_{1:t})}{P(o_{1:t})}\propto P(i_t,o_{1:t})=\alpha_t(i_t) \end{aligned}
  • Smoothing:公式为 P(i_t|o_{1:T}),同样可以化简为\begin{aligned} P(i_t|o_{1:T})=\frac{P(i_t,o_{1:T})}{P(o_{1:T})}&\propto P(i_t,o_{1:T})\\ &=P(i_t,o_{1:t},o_{t+1:T})\\ &=\underset{\underset{\beta_t(i_t)}{\underbrace{P(o_{t+1:T}|i_t)}}}{\underbrace{P(o_{t+1:T}|i_t,o_{1:t})}}\cdot \underset{\alpha_t(i_t)}{\underbrace{P(i_t,o_{1:t})}}\\ &=\beta_t(i_t)\cdot\alpha_t(i_t) \end{aligned},又被称为前向后向算法;
  • Prediction:公式为P(i_{t+1}|o_{1:t}),可以化简为\begin{aligned} P(i_{t+1}|o_{1:t})&=\sum_{i_t}P(i_{t+1},i_t|o_{1:t})\\ &=\sum_{i_t}\underset{\underset{a_{i_{t},i{t+1}}}{\underbrace{P(i_{t+1}|i_t)}}}{\underbrace{P(i_{t+1}|i_t,o_{1:t}})}\cdot \underset{\text{Filtering:}\ \propto \alpha_t(i_t)}{\underbrace{P(i_t|o_{1:t})}}\\ &\propto \sum_{i_t}a_{i_{t},i{t+1}}\cdot \alpha_t(i_t) \end{aligned}

15. 线性动态系统(卡曼滤波,Kalman Filter)

15.1. 概述

15.1.1. 背景介绍 

  • 变量随时间变化的系统叫做动态系统,其中隐变量取值离散的是隐马尔可夫模型(HMM),而隐变量取值连续的分为线性动态系统和粒子滤波(Particular Filter),前者各变量和变量各时刻满足线性关系和高斯分布,后者则不满足;
  • 线性关系:时刻之间:z_t=A\cdot z_{t-1}+B+\varepsilon;变量之间:x_t=C\cdot z_{t}+D+\delta;其中 \varepsilon,\delta 是噪声,在线性动态系统中服从高斯分布;

15.1.2. 模型定义

  • 设噪声 \varepsilon \sim N(0,Q),\delta \sim N(0,R),同时设初始状态z_1\sim N(\mu_1,\Sigma_1),则可以通过高斯分布依次计算 z_t|z_{t-1}\sim N(A\cdot z_{t-1}+B,Q),x_t|z_{t}\sim N(C\cdot z_{t}+D,R) 得到所有变量的分布,因此模型参数为 \theta =(A,B,C,D,Q,R,\mu_1,\Sigma_1)
  • 线性动态系统和HMM一样,也服从一阶齐次马尔可夫假设和观测独立假设(见14.2.1.);

15.2. Filtering

  • 首先,filtering公式为P(z_t|x_{1:t}),是一种Online算法,通过 t 时刻及之前的观测变量来推测当前状态,这既是模型求解(Learning)的必要信息,同时也是对predictionP(z_t|x_{1:t-1})结果的修正(Correction);
  • filtering也可以递推求解,推导为:\begin{aligned} P(z_t|x_{1:t})&=\frac{P(z_t,x_{1:t})}{P(x_{1:t})}=\frac{\overset{P(x_t|z_t)}{\overbrace{P(x_t|z_t,x_{1:t-1})}}\cdot P(z_t,x_{1:t-1})}{P(x_{1:t})}\\ &=\frac{P(x_t|z_t)\cdot P(z_t|x_{1:t-1})\cdot P(x_{1:t-1})}{P(x_{1:t})}\\ &=\underset{\text{Emission Probability}}{\underbrace{P(x_t|z_t)}}\cdot \underset{\text{Prediction of t-1}}{\underbrace{P(z_t|x_{1:t-1})}}\cdot \underset{\text{Constant}}{\underbrace{\frac{P(x_{1:t-1})}{P(x_{1:t})}}}\\ \end{aligned},其中Prediction可以求解为:\begin{aligned} P(z_t|x_{1:t-1})&=\int_{z_{t-1}}P(z_t,z_{t-1}|x_{1:t-1})dz_{t-1}\\ &=\int_{z_{t-1}}\underset{\underset{\text{State Transition Probability}}{\underbrace{P(z_t|z_{t-1})}}}{\underbrace{P(z_t|z_{t-1},x_{1:t-1})}}\cdot \underset{\text{Filtering of t-1}}{\underbrace{P(z_{t-1}|x_{1:t-1})}}dz_{t-1}\\ \end{aligned},因此可以通过t-1时刻的filtering计算t时刻的prediction,进而计算t时刻的filtering,循环迭代进行求解;
  • 总结来说,Filtering的过程可以写作以下步骤:
    • t=1,\begin{aligned} P(z_1|x_1)=\frac{P(z_1,x_1)}{P(x_1)}=\frac{P(x_1|z_1)\cdot P(z_1)}{\int_{z_1}P(x_1|z_1)\cdot P(z_1)dz_1} \end{aligned},并用上面的公式计算Prediction\begin{aligned} P(z_2|x_1) \end{aligned}
    • t=2,根据上面的递推关系计算\begin{aligned} P(z_2|x_1,x_2) \end{aligned}(作为filtering的更新,update,以及prediction的修正,correction)和\begin{aligned} P(z_3|x_1,x_2) \end{aligned}(新的prediction);
    • ……
    • t=T,计算\begin{aligned} P(z_T|x_{1:T-1}) \end{aligned}
  • 在整个计算过程中,Filtering和Prediction的计算结果均为高斯分布,其参数(均值和标准差)可以通过之前时刻的计算结果以及模型参数算出,推导过程参考2. 高斯分布中高斯分布的性质,这里不再详细介绍;

16. 粒子滤波(Particular Filtering)

16.1. 背景介绍

  • 在15. 线性动态系统(卡曼滤波)中,可以通过filtering和prediction交替迭代的方法求得隐状态序列,但应对非线性非高斯模型,无法求出解析解,只能通过采样模拟的方法求出近似解,如此便诞生了粒子滤波;
  • 粒子滤波应用范围非常之广,且已经产生了许多研究成果和变体,这里只介绍最基本的粒子滤波算法;

16.2. 算法介绍

16.2.1. 重要度采样及针对filtering的改进 —— 序贯重要度采样(Sequential Importance Sampling,SIS)

  • 首先回顾重要度采样,给定一个复杂的分布p(x),需要计算期望E[f(x)]时,直接从p(x)中采样较为困难,此时可以设置一个简单的建议分布q(x),此时E[f(x)]=\int_z \left[f(x)p(x)\right]dz=\int_z \left[f(x)\cdot\frac{p(x)}{q(x)}\cdot q(x)\right]dz,采样变得容易;
  • 将重要度采样应用于动态系统的filtering问题时,对每个时刻 t 和每个样本x_t^{(i)},都需要计算重要度\begin{aligned} w_t^{(i)}=\frac{p(z_t^{(i)}|x_{1:t}^{(i)})}{q(z_t^{(i)}|x^{(i)}_{1:t})} \end{aligned},而p(z_t^{(i)}|x_{1:t}^{(i)})的计算是复杂的,还要计算N次(N个样本);
  • 将filtering问题的计算目标进行调整,由 p(z_t^{(i)}|x_{1:t}^{(i)}) 改为 p(z_{1:t}^{(i)}|x_{1:t}^{(i)}) (同样可以达到目标),此时重要度为\begin{aligned} w_t^{(i)}\propto \frac{p(z_{1:t}^{(i)}|x_{1:t}^{(i)})}{q(z_{1:t}^{(i)}|x^{(i)}_{1:t})} \end{aligned},分别推导分子和分母的递推式,如果可以建立递推关系,则可以简化运算:
    • 对于分子,可以推导如下\begin{aligned} p(z_{1:t}^{(i)}|x_{1:t}^{(i)})&=\frac{p(z_{1:t}^{(i)},x_{1:t}^{(i)})}{p(x_{1:t}^{(i)})}\\ &=\frac{\overset{\overset{\text{Emission Probability}}{\overbrace{p(x_t^{(i)}|z_t^{(i)})}}}{\overbrace{p(x_t^{(i)}|z_{1:t}^{(i)},x_{1:t-1}^{(i)})}}\cdot p(z_{1:t}^{(i)},x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\\ &=\frac{p(x_t^{(i)}|z_t^{(i)})\cdot \overset{\overset{\text{State Transition Probability}}{\overbrace{p(z_t^{(i)}|z_{t-1}^{(i)})}}}{\overbrace{p(z_t^{(i)}|z_{1:t-1}^{(i)},x_{1:t-1}^{(i)})}}\cdot p(z_{1:t-1}^{(i)},x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\\ &=\frac{p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)})\cdot p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})\cdot p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\\ &=\left [\frac{p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\cdot p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)}) \right ]\cdot p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})\\ \end{aligned}
    • 对于分母,可以设 q(z_{1:t}^{(i)}|x_{1:t}^{(i)}) 满足 q(z_{1:t}^{(i)}|x_{1:t}^{(i)})=q(z_{t}^{(i)}|z_{1:t-1}^{(i)},x_{1:t}^{(i)})\cdot q(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})
    • 将分子分母带入得\begin{aligned} w_t^{(i)}&\propto \frac{p(z_{1:t}^{(i)}|x_{1:t}^{(i)})}{q(z_{1:t}^{(i)}|x^{(i)}_{1:t})} \\ &=\frac{\left [\frac{p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}\cdot p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)}) \right ]\cdot p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}{q(z_{t}^{(i)}|z_{1:t-1}^{(i)},x_{1:t}^{(i)})\cdot q(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}\\ &=\underset{\text{Constant}}{\underbrace{\frac{p(x_{1:t-1}^{(i)})}{p(x_{1:t}^{(i)})}}}\cdot \left [\frac{p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)}) }{\underset{q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})}{\underbrace{q(z_{t}^{(i)}|z_{1:t-1}^{(i)},x_{1:t}^{(i)})}}}\right ]\cdot \underset{w_{t-1}^{(i)}}{\underbrace{\frac{p(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}{q(z_{1:t-1}^{(i)}|x_{1:t-1}^{(i)})}}} \\ &\propto \left [\frac{p(x_t^{(i)}|z_t^{(i)})\cdot p(z_t^{(i)}|z_{t-1}^{(i)}) }{q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})}\right ]\cdot w_{t-1}^{(i)} \end{aligned},如此便建立了递推计算关系,可以通过蒙特卡洛的重要度采样方法计算filtering问题;
  • 为了采样能收敛于真正的期望,需要对重要度进行归一化,即保证\sum_{i=1}^N w_t^{(i)}=1,这也就是为什么递推关系中可以化简常数,只要得到正比的关系即可;
  • 上述算法即为基础的粒子滤波算法,其中对状态值的每一次抽样结果z_t^{(i)}\sim q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})被称为一个粒子;

16.2.2. 重采样

  • 在应用上述粒子滤波算法的过程中,很多时候存在粒子退化问题,即只有少数状态拥有真正有效的重要度 w_t^{(i)} ,其他大多数状态的重要度接近0,但算法每次都需要对所有 z_t^{(i)} 以 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)}) 的概率抽样,而 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)}) 与重要度 w_t^{(i)} 之间没有关联,这意味着可能有很多情况 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)}) 较大,w_t^{(i)} 较小,抽样结果充斥着不重要的信息,或是反过来,有效信息在抽样结果中占比很低,这些情况都会使计算资源浪费,模型难以收敛,效果变差等等;
  • 重采样的思想用于解决上述问题,具体来说,当使用概率 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)}) 完成抽样后,计算每个粒子的重要度,并对重要度计算CDF,再在0-1的均匀分布上抽取N个分裂因子\mu_t^{(i)},每个 \mu_t^{(i)} 落在CDF的哪个粒子区间内,就把哪个粒子复制一次(有时也叫粒子分裂);
  • 例如,抽样得到的结果为z1,z2,z3(从小到大排列),三者的重要度分别为0.1,0.1,0.8,计算CDF为0-0.1、0.1-0.2、0.8-1.0三段,则需要抽取三个分裂因子,假设为0.15,0.38,0.54,第一个落入z2区间(0.1-0.2),则z2复制一次,另外两个落入z3区间(0.2-0.8),则z3复制两次,通过这个方法尽可能让每次抽样的粒子拥有接近1/N的重要度(或者说重要度高的粒子分裂的次数也多,更加能被有效利用,而重要度低的粒子可以适当舍弃);
  • 当样本的维度过高时粒子退化问题将更加可怕,此时重采样的必要性更高;

16.2.3. Sampling-Importance-Resampling(SIR)

  •  除了重采样外,还有一种方法可以解决16.2.2.中提到的粒子退化问题,就是设置一个更加合理的 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)}) ,有一个方法是令 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})=p(z_t^{(i)}|z_{t-1}^{(i)}),此时重要度的递推关系为w_t^{(i)}=p(x_t^{(i)}|z_t^{(i)})\cdot w_{t-1}^{(i)},这种改进版的粒子滤波再加上重采样,可以进一步提升算法性能,叫做SIR;
  • 为什么要令 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})=p(z_t^{(i)}|z_{t-1}^{(i)}),其实是逻辑上的考虑,filtering的目标为计算p(z_t^{(i)}|x_{1:t}^{(i)})(或p(z_{1:t}^{(i)}|x_{1:t}^{(i)})),而这个问题中粒子滤波的方案是要对 z_t^{(i)} 进行采样来计算均值,一个很简单的思想是,在已知 x_t^{(i)} 的前提下,p(x_t^{(i)}|z_t^{(i)}) 越大的 z_t^{(i)} 很多时候越有可能成为真正的状态,因此应该赋予其越高的权重,这就是要设置一个特殊的 q(z_{t}^{(i)}|z_{t-1}^{(i)},x_{1:t}^{(i)})=p(z_t^{(i)}|z_{t-1}^{(i)}) 的原因,就是为了凑出w_t^{(i)}=p(x_t^{(i)}|z_t^{(i)})\cdot w_{t-1}^{(i)},然而这是一个朴素的思想,后人对这个建议分布的取法进行了很多研究,这里只介绍最基本的粒子滤波,因此不做赘述;

17. 条件随机场(Conditional Random Field,CRF)

17.1. 背景

  • 机器学习分类模型中,有硬分类和软分类两种主流思想,其中硬分类模型有支持向量机SVM(最大化几何间隔)、感知机PLA(误差驱动优化)、线性判别分析LDA(类间大,类内小)等,通过严格的决策边界来区分多类样本;而软分类模型有逻辑回归(及Softmax回归等)和朴素贝叶斯模型等;
  • 同时机器学习分类模型也有概率生成模型和概率判别模型两种思想:
    • 软分类中的逻辑回归是概率判别模型,优化思想是最大熵(见8.4.,给定均值和方差而无其他信息时,熵值最大的分布是高斯分布,同时最大熵思想是使模型贴近已有样本,而样本无法提供的信息使其尽量服从高斯分布);
    • 朴素贝叶斯是概率生成模型,并从中衍生出了隐马尔可夫模型(可以看做时序版本的朴素贝叶斯,同时认为每个观测变量每个时刻都有一个隐状态,为了求解这个隐状态,引入了齐次马尔可夫假设和观测独立假设);
  • 从软分类的概率生成模型和概率判别模型的研究成果中,结合了最大熵思想和隐马尔可夫模型的优点,诞生了最大熵马尔可夫模型(Maximum Entropy Markov Model,MEMM),这是一个概率判别模型,即求解的问题是P(Y|X),同时是一种新的马尔可夫模型,其引入了一个全局变量,使得观测独立假设被打破,这在实际问题中更加合理,样本信息被更多地利用;

  • 上述图中阴影部分为MEMM求解过程的的一个步骤,即P(y_2|x_2,y_1),但这样的条件概率求出以后需要进行归一化才能使其符合概率分布特性(求和为1),这样的操作为局部归一化,但每次归一化时没有考虑全局信息,会导致标签偏差问题(Label Bias Problem),而将MEMM的概率图模型变为无向图即可解决这个问题,因此诞生了条件随机场CRF;
  • 标签偏差问题是在Learning问题有效收敛的前提下,decoding任务中出现了明显偏差的现象,因为MEMM和CRF中的隐状态常常是标签(用y_t表示,例如词性标注中具体词语是观测变量,词性如名词、动词、形容词为隐变量),因此叫做标签偏差问题;
  • 举例来说:
    • 在经典的John Lafferty论文中的例子是:
      • 假如样本数据有15个 r→i→d,5个 r→o→d,忽视发射矩阵,仅考虑状态转移,那么在Learning阶段,理想的概率图模型应该如下图所示,此时在decoding阶段,给定观测变量 r→o→d,得到 i2→i4→i5 的隐变量序列;

      • 但如果在o2位置有更多的隐变量取值空间,则Learning阶段可能学习到下图所示的,则在decoding阶段给定观测变量 r→o→d时,会在第一步走向i1的路径,从而得到 i1→i3→i5 的隐变量序列,这显然是不正确的;

    •  西湖大学NLP课程中给出例子为:如下图,根据频率统计,d4出现了2次,比d6频率更高,但根据维特比算法求出二者的概率,却是d6比d4概率高,因为l3有两种转移的可能,即便两种可能出现的频率更高(l2→l1有2次,l3→l3有2次),但经过归一化后都为0.5的概率值,而l2只有一种转移的可能,即便出现频率低(只有l2→l2有1次),但经过归一化后概率值依然为1,因此d6算得的概率更高,这种现象本质上是局部归一化导致的;
    • Toutanova在论文中也讨论了标签偏差现象的原因,提出有一种情况是发射概率过高,例如下图同样是词性标注问题,Will to fight是句子,是观测变量,NN MODAL TO VB是词性(NN指名词,MODAL指情态动词,TO指to专属词性,VB指动词),而MEMM在decoding任务中需要计算的是P(y_t=\text{[TO]}\ |\ y_{t-1},x_t=\text{[to]}),而因为TO是to的专属词性,因此只要见到t时刻观测变量为to即可直接标注词性为TO,无需纳入前一时刻Will的词性作为考虑,这样的模型Learning阶段会无法学习到NN→TO和MODAL→TO之间的状态转移关系,而事实上情态动词后面是基本不搭配to的,因此对Will的词性出现了标签偏差问题,这样的问题本质上也是对to→TO的发射概率进行了局部归一化所导致;
  • 因此CRF将MEMM改为无向图连接,使得局部归一化变为了全局归一化,解决了标签偏差问题,同时打破了齐次马尔可夫假设,使得模型更加有效地利用信息;

17.2. CRF的概率密度函数——Decoding

17.2.1. 无向图模型的概率密度函数

  • 团:对于一个无向图,其中某个子图满足各个节点相互连通,则作为原无向图的一个团,而对于一个团,若其无法添加一个新节点使其变为一个更大的团(也不能等于原无向图),则称为最大团;
  • 概率因子分解法:
    • 对于一个无向图形式的概率图模型,要求解其概率密度函数时,常常用到概率因子分解法,具体来说是找到图模型的所有最大团,对每个最大团求解势函数,相乘并归一化得到概率密度函数;
    • 公式:设概率图模型中所有变量集合为X=\{ x_1,x_2,\cdots,x_m\},每个 x_i 为一个节点,对图划分得到K个最大团C=\{ c_1,\cdots,c_K \},每个团的势函数为\Psi _i(x_j\in c_i),则联合概率密度为P(X)=\frac{1}{Z}\prod_{i=1}^K\Psi_i(x_j\in c_i),其中归一化因子Z=\int_X \prod_{i=1}^K\Psi_i(x_j\in c_i)\ dX
    • 举例:
  • 可以进一步简化概率因子分解的公式:P(X)=\frac{1}{Z}\prod_{i=1}^K\Psi_i(x_j\in c_i)=\frac{1}{Z}\prod_{i=1}^K\exp\{-E_i(x_j\in c_i)\},其中E(x)为能量函数,再设F(x)=-E(x),则P(X)=\frac{1}{Z}\prod_{i=1}^K\exp\{F_i(x_j\in c_i)\}=\frac{1}{Z}\exp\left [\sum_{i=1}^K\{F_i(x_j\in c_i)\} \right ]

17.2.2. CRF的概率密度函数(参数形式)

  • 下图是一个典型的CRF概率图,其中每个最大团是c_t=\{ y_{t-1},y_t,x_{1:T} \},三个节点组成;

  • 由概率因子分解可得,P(Y|X)=\frac{1}{Z}\exp\left [\sum_{t=1}^T\{F_t(y_{t-1},y_t,x_{1:T}\} \right ],因为decoding问题中观测变量x_{1:T}给定,则可以写作F_t(y_{t-1},y_t,x_{1:T})=F(y_{t-1},y_t)=F_a(y_{t-1})+F_b(y_{t})+F_c(y_{t-1},y_t)
  • 由于在计算 F_t 时,y_{t-1} 已知,所以只关注 F_b(y_{t}) 和 F_c(y_{t-1},y_t) ;其中 F_b(y_{t}) 所表达的信息比较类似HMM中的发射矩阵,而 F_c(y_{t-1},y_t) 类似状态转移矩阵,可以通过这个思路来设计合理的函数使得模型更加合理和有效,例如设 F_b(y_t)=\sum_{i=1}^K\eta_i\cdot g_i(y_t)F_a(y_{t-1})+F_c(y_{t-1},y_t)=\sum_{i=1}^L\lambda_i\cdot f_i(y_{t-1},y_t),其中根据应用场景设计f_i(\cdot)g_i(\cdot),而将 \lambda_i 和 \eta_i 作为可学习参数,这就完成参数形式的概率密度函数;
  • 整理参数形式的概率密度函数为:\begin{aligned} P(Y|X)&=\frac{1}{Z}\exp\left \{\sum_{t=1}^T\left [\sum_{i=1}^K\eta_i\cdot g_i(y_t)+\sum_{i=1}^L\lambda_i\cdot f_i(y_{t-1},y_t) \right ] \right \}\\ &=\frac{1}{Z}\exp\left \{\sum_{t=1}^T\left [\vec{\eta}^T\cdot \vec{g}(y_t)+\vec{\lambda}^T\cdot \vec{f}(y_{t-1},y_t) \right ] \right \}\\ \end{aligned},其中\vec{\eta}=\{ \eta_1,\eta_2,\cdots,\eta_K \}^T\vec{g}(\cdot)=\{g_1(\cdot),g_2(\cdot),\cdots,g_K(\cdot) \}^T, \vec{\lambda} 和 \vec{f}(\cdot) 同理;

17.3. CRF的边缘概率密度——Smoothing

  • 在Decoding中求出概率密度函数P(Y|X),还需要求解边缘概率密度函数P(y_t|X),先直接写出公式\begin{aligned} P(y_t|X)=\int_{y_{1},\cdots,y_{t-1},y_{t+1},\cdots,y_T}P(Y|X)=\int_{y_{1},\cdots,y_{t-1},y_{t+1},\cdots,y_T}\frac{1}{Z}\prod_{i=1}^T\Psi_i(y_{i-1},y_i,X) \end{aligned},这个连乘与积分的计算是非常复杂的(指数级时间复杂度),需要使用动态规划方法递推求解;
  • 拆解\begin{aligned} P(y_t|X)=\int_{y_{1},\cdots,y_{t-1}}\int_{y_{t+1},\cdots,y_T}\frac{1}{Z}\prod_{i=1}^T\Psi_i(y_{i-1},y_i,X) \end{aligned},则可以令\begin{aligned} \alpha_t(y_t)=\int_{y_{1},\cdots,y_{t-1}}\prod_{i=1}^t\Psi_i(y_{i-1},y_i,X) \end{aligned}\begin{aligned} \beta_t(y_t)=\int_{y_{t+1},\cdots,y_{T}}\prod_{i=t+1}^T\Psi_i(y_{i-1},y_i,X) \end{aligned},因此\begin{aligned} P(y_t=i|X)=\alpha_t(i)\cdot \beta_t(i)\end{aligned}
  • 建立递推关系(也是一种前向后向算法):
    • \begin{aligned} \alpha_t(y_t)&=\int_{y_{1},\cdots,y_{t-1}}\prod_{i=1}^t\Psi_i(y_{i-1},y_i,X) \\ &=\int_{y_{t-1}}\Psi_t(y_{t-1},y_t,X)\left [ \int_{y_{t-2}}\Psi_{t-1}(y_{t-2},y_{t-1},X)\left [ \cdots\left [ \int_{y_{2}}\Psi_2(y_{1},y_2,X)\left [ \int_{y_{1}}\Psi_1(y_{0},y_1,X) \right ] \right ] \right ] \right ] \\ &=\int_{y_{t-1}}\Psi_t(y_{t-1},y_t,X)\cdot \alpha_{t-1}(y_{t-1}) \end{aligned}
    • \begin{aligned} \beta_t(y_t)&=\int_{y_{t+1},\cdots,y_{T}}\prod_{i=t+1}^T\Psi_i(y_{i-1},y_i,X) \\ &=\int_{y_{t+1}}\Psi_{t+1}(y_{t},y_{t+1},X)\left [ \int_{y_{t+2}}\Psi_{t+2}(y_{t+1},y_{t+2},X)\left [ \cdots\left [ \int_{y_{T-1}}\Psi_{T-1}(y_{T-2},y_{T-1},X)\left [ \int_{y_{T}}\Psi_T(y_{T-1},y_T,X) \right ] \right ] \right ] \right ] \\ &=\int_{y_{t+1}}\Psi_{t+1}(y_{t},y_{t+1},X)\cdot \beta_{t+1}(y_{t+1}) \end{aligned}

17.4. CRF的参数估计——Learning

  • CRF的Learning就是根据已有的样本数据来最大化似然函数,即\theta=<\vec{\lambda},\vec{\eta}>=\underset{\vec{\lambda},\vec{\eta}}{\arg\max}\ \prod_{i=1}^NP(Y^{(i)}|X^{(i)}),等价于求对数似然函数\begin{aligned} <\vec{\lambda},\vec{\eta}>&=\underset{\vec{\lambda},\vec{\eta}}{\arg\max}\ \log\prod_{i=1}^N\frac{1}{Z^{(i)}}\exp\left \{\sum_{t=1}^T\left [\vec{\eta}^T\cdot \vec{g}(y_t^{(i)})+\vec{\lambda}^T\cdot \vec{f}(y_{t-1}^{(i)},y_t^{(i)}) \right ] \right \}\\ &=\underset{\vec{\lambda},\vec{\eta}}{\arg\max}\ \sum_{i=1}^N\left \{-\log Z^{(i)}(X^{(i)},\vec{\lambda},\vec{\eta})+\sum_{t=1}^T\left [\vec{\eta}^T\cdot \vec{g}(y_t^{(i)})+\vec{\lambda}^T\cdot \vec{f}(y_{t-1}^{(i)},y_t^{(i)}) \right ] \right \}\\ \end{aligned}
  • 对上述最大化似然函数,可以使用梯度上升法逼近数值解,即\begin{aligned} \bigtriangledown _{\vec{\lambda}}=\sum_{i=1}^N- \frac{1}{Z^{(i)}}\frac{\partial Z^{(i)}}{\partial \vec{\lambda}}+\sum_{i=1}^N\sum_{t=1}^T\vec{f}(y_{t-1}^{(i)},y_t^{(i)}) \end{aligned},其中\begin{aligned} Z^{(i)}=\int_Y\exp\left\{\sum_{t=1}^T\left [\vec{\eta}^T\cdot \vec{g}(y_t^{(i)})+\vec{\lambda}^T\cdot \vec{f}(y_{t-1}^{(i)},y_t^{(i)}) \right ]\right \} dY \end{aligned},则\begin{aligned} \sum_{i=1}^N- \frac{1}{Z^{(i)}}\frac{\partial Z^{(i)}}{\partial \vec{\lambda}} &=\sum_{i=1}^N- \int_Y\underset{P(Y^{(i)}|X^{(i)})}{\underbrace{\frac{ \exp\left\{\sum_{t=1}^T\left [\vec{\eta}^T\cdot \vec{g}(y_t^{(i)})+\vec{\lambda}^T\cdot \vec{f}(y_{t-1}^{(i)},y_t^{(i)}) \right ]\right \}}{Z^{(i)}}}}\cdot \vec{f}(y_{t-1}^{(i)},y_t^{(i)}) dY\\ &=\sum_{i=1}^N- E\left [\vec{f}(y_{t-1}^{(i)},y_t^{(i)}) \right ] \end{aligned}
  • 因为\begin{aligned} \int_YP(Y|X)\cdot\vec{f}(y_{t-1},y_t)dY&=\int_{y_1}\cdots\int_{y_T}\left [P(y_1,\cdots,y_T|X)\cdot\vec{f}(y_{t-1},y_t) \right ]\ dy_1\cdots dy_T \\ &=\int_{y_{t-1}}\int_{y_t}\left [P(y_{t-1},y_t|X)\cdot\vec{f}(y_{t-1},y_t) \right ]\ dy_{t-1}dy_t\end{aligned},其中P(y_{t-1},y_t|X)可以用前向后向算法求得(见17.3. CRF的边缘概率密度——Smoothing),将该式带入上一步,可以求得\begin{aligned} \bigtriangledown _{\vec{\lambda}}\end{aligned}\begin{aligned} \bigtriangledown _{\vec{\eta}}\end{aligned}也可以同理计算;

18. 高斯网络(Gaussian Network,GN)

18.1. 背景

  • 高斯网络也属于概率图模型中的一种,变量离散的概率图模型按照有向和无向分为贝叶斯网络(Bayesian Network,BN)和马尔可夫网络(Markov Network,MN),变量连续的概率图模型是高斯网络,其中假设每个变量都遵循高斯分布,同样按照有向无向分为贝叶斯高斯网络(Gaussian Bayesian Network,GBN)和马尔可夫高斯网络(Gaussian Markov Network,GMN);
  • 设总随机变量 \vec{x} 是一个p维向量,即高斯网络中有p个节点,则存在总均值 \vec{\mu} 和协方差矩阵 \Sigma ,此时概率密度函数为P(\vec{x})=\frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}\cdot \exp\{-\frac{1}{2}(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu})\}
  • 对协方差矩阵\Sigma=\begin{pmatrix} \sigma_{11} & \cdots & \sigma_{1p} \\ \vdots & \ddots & \vdots \\ \sigma_{p1} & \cdots & \sigma_{pp} \end{pmatrix},有x_i\perp x_j\Leftrightarrow \sigma_{ij}=0,表示变量之间的(完全)局部独立性,同时存在精度矩阵\Lambda =\Sigma^{-1}=\begin{pmatrix} \lambda_{11} & \cdots & \lambda_{1p} \\ \vdots & \ddots & \vdots \\ \lambda_{p1} & \cdots & \lambda_{pp} \end{pmatrix},有x_i\perp x_j|\{ x_{k,k\neq i,k\neq j}\}\Leftrightarrow \lambda_{ij}=0,表示变量之间的条件独立性(原因在18.3. 高斯马尔可夫网络 中有详细介绍);

18.2. 高斯贝叶斯网络(GBN)

  •  在贝叶斯网络中,变量的因子分解为P(\vec{x})=\prod_{i=1}^pP(x_i|\vec{x}_{parent(i)}),其中\vec{x}_{parent(i)}是父节点变量,而高斯贝叶斯网络中的每个节点x_i|\vec{x}_{parent(i)}\sim N(\mu_i+W_i^T\vec{x}_{parent(i)},\sigma_i^2),即每个节点为高斯分布,且与其父节点呈线性关系;
  • 用向量表示为\vec{x}-\vec{\mu}=W(\vec{x}-\vec{\mu})+S\cdot \vec{\epsilon},其中W=\left [ w_{ij} \right ]_{p\times p}w_{ij}表示变量x_i与其父节点x_j的线性关系系数(若不是父节点则为0),S=\text{diag}(\sigma_i)_{p\times p}为标准差矩阵,\vec{\epsilon}\sim N(0, 1)为一组随机数;
  • 通过化简可得\vec{x}-\vec{\mu}=(I-W)^{-1}S\cdot \vec{\epsilon}(假设(I-W))可逆,则cov(\vec{x}-\vec{\mu})=\left [(I-W)^{-1}S \right ]^T\left [ (I-W)^{-1}S \right ]cov(\vec{\epsilon})\\=\left [(I-W)^{-1}S \right ]^T\left [ (I-W)^{-1}S \right ]

18.3. 高斯马尔可夫网络(GMN)

  • 高斯马尔可夫网络是变量连续的高斯网络中的无向图模型,其建模的思路主要是变量之间的条件独立性,概率密度函数可以统一为P(\vec{x})=\frac{1}{Z}\prod_{i=1}^p\left [\Psi_i(x_i)\cdot \prod_{j}\Psi_{i, j}(x_i,x_j) \right ],其中\Psi_i(x_i)表示了每个变量(节点)的势函数,而\Psi_{i, j}(x_i,x_j)表示与其相关的变量之间相关性(边)的势函数;
  • 将多维高斯分布的概率密度函数写成高斯马尔可夫网络的形式,\begin{aligned} P(\vec{x})&=\frac{1}{(2\pi)^{\frac{p}{2}}|\Sigma|^{\frac{1}{2}}}\cdot\exp\{-\frac{1}{2}(\vec{x}-\vec{\mu})^T\Sigma^{-1}(\vec{x}-\vec{\mu})\} \\ &\propto \exp\{-\frac{1}{2}(\vec{x}-\vec{\mu})^T\Lambda (\vec{x}-\vec{\mu})\} \\ &= \exp\{-\frac{1}{2}\left (\vec{x}^T\Lambda \vec{x} - 2\vec{\mu}^T\Lambda \vec{x} +\vec{\mu}^T\Lambda \vec{\mu} \right ) \} \\ &\propto \exp\{-\frac{1}{2}\vec{x}^T\Lambda \vec{x} + \left (\Lambda\vec{\mu} \right )^T \vec{x} \} \\ \end{aligned}
    • 其中对每个变量x_i,其节点势函数为 A(-\frac{1}{2}\lambda_{ii}x_i^2+h_ix_i) ,其中 A 为常数, \lambda_{ii} 为精度矩阵中的元素, h_i 为向量 \Lambda\vec{\mu} 中的元素;
    • 边的势函数为B(\lambda_{ij}x_ix_j),表示了节点之间的相关性,因此若 \lambda_{ij}=0,则两个变量满足条件独立性;
  • 利用精度矩阵 \Lambda ,可以将无向图和多维高斯分布有效联系,建立一个可以建模条件独立性的模型,其中每个变量的边缘概率密度为x_i|x_{j,j\neq i}\sim N(\sum_{j\neq i}\frac{\lambda_{ij}}{\lambda_{ii}}x_j,\lambda_{ii}^{-1}),这个结论可以通过2.5. 高斯分布的条件概率分布 中的计算得到;

19. 高斯过程

19.1. 高斯过程简介

  • 对于一个无限维的随机变量 \vec{\xi},其连续域 T 是连续的,可以是时间域(动态系统)也可以是空间域(无限维特征),以时间域为例,在时间上的不同点观测得到观测序列,时间点记作 t_1,t_2,\cdots,t_n \in T,则该变量的一个样本为一条观测序列\left \{\xi _{t_1},\xi _{t_2},\cdots,\xi _{t_n} \right \}
  • 对于高斯过程中的动态随机变量,其每个时刻的观测值服从当前时刻的一个高斯分布,即\xi_{t}\sim N(\mu_t,\sigma_t^2),而整个观测序列服从多维高斯分布 \vec{\xi}\sim N(\vec{\mu},\Sigma)
  • 该高斯过程通常记作GP(m(t),K(s,t)),其中:
    • m(t) 为均值函数,m(t)=\mu_t=E[\xi_t]
    • K(s,t) 为协方差函数,又叫核函数,K(s,t)=Cov(\xi_s,\xi_t)=E[(\xi_s-m(s))\cdot (\xi_t-m(t))]

19.2. 高斯过程回归

19.2.1. 从贝叶斯线性回归到高斯过程回归

  • 先回顾贝叶斯线性回归:
    • 根据样本数据X=\left[\hat{\vec{x}}_1,\hat{\vec{x}}_2,\cdots,\hat{\vec{x}}_N\right]^T,Y=\left[ y_1,y_2,\cdots,y_N \right]^T,设存在映射系数矩阵,使得y_n=\hat{\vec{x}}_n^TW+\epsilon,\epsilon\in N(0,\sigma^2),则可以解得W|X,Y\sim N(\vec{\mu}_w,\Sigma_w),其中\vec{\mu}_w=\sigma^{-2}A^{-1}X^TY,\Sigma_w=A^{-1},A=\sigma^{-2}X^TX+\Sigma_p^{-1}具体推导过程见3.3.4.);
    • 进一步地,给定一个新样本 \hat{\vec{x^*}},计算其预测结果 f(\hat{\vec{x^*}})=W^T\hat{\vec{x^*}}\sim N({\hat{\vec{x^*}}}^T\vec{\mu}_w,\ {\hat{\vec{x^*}}}^T\Sigma_w{\hat{\vec{x^*}}})y^*=W^T\hat{\vec{x^*}}+\epsilon\sim N({\hat{\vec{x^*}}}^T\vec{\mu}_w,\ {\hat{\vec{x^*}}}^T\Sigma_w{\hat{\vec{x^*}}}+\sigma^2),然而贝叶斯线性回归只能解决线性问题,对于现实世界中更多的非线性问题无能为力;
  • 为了应对非线性问题,常用的方法是使用一个映射函数,将低维特征映射到高维空间,变成一个高维空间的线性可分问题,我们定义一个非线性转换 \Phi =\left [ \phi(\hat{\vec{x}}_1), \phi(\hat{\vec{x}}_2), \cdot, \phi(\hat{\vec{x}}_N) \right ]^T,则预测函数变为f(\hat{\vec{x^*}})=W^T\phi(\hat{\vec{x^*}})\sim N(\sigma^{-2}{\phi(\hat{\vec{x^*}}})^T(A^{-1}\Phi^TY),\ \phi({\hat{\vec{x^*}}})^TA^{-1}\phi({\hat{\vec{x^*}}})),其中A=\sigma^{-2}\Phi^T\Phi+\Sigma_p^{-1}
  • 然而这里又会遇到一个问题,A^{-1}计算非常困难,需要用到Woodbury Formula公式(公式非常复杂不作介绍,这里只讲简化思路),简化思路如下:
    • \begin{aligned} A&=\sigma^{-2}\Phi^T\Phi+\Sigma_p^{-1}\\ A\Sigma_p&=\sigma^{-2}\Phi^T\Phi\Sigma_p+I\\ A\Sigma_p\Phi^T&=\sigma^{-2}\Phi^T\Phi\Sigma_p\Phi^T+\Phi^T=\sigma^{-2}\Phi^T(\underset{\bar{K}}{\underbrace{\Phi\Sigma_p\Phi^T}}+\sigma^{2}I) \\ \Sigma_p\Phi^T&=\sigma^{-2}A^{-1}\Phi^T(\bar{K}+\sigma^{2}I) \\ \Sigma_p\Phi^T(\bar{K}+\sigma^{2}I)^{-1}&=\sigma^{-2}A^{-1}\Phi^T\\ \left [\phi({\hat{\vec{x^*}}})^T\Sigma_p\Phi^T \right ](\bar{K}+\sigma^{2}I)^{-1}Y&=\sigma^{-2}\phi({\hat{\vec{x^*}}})^TA^{-1}\Phi^TY \end{aligned},等式右侧正好是f(\hat{\vec{x^*}})均值,也就是期望;
    • 用相同方法可以得到f(\hat{\vec{x^*}})方差\left [\phi({\hat{\vec{x^*}}})^T\Sigma_p\phi({\hat{\vec{x^*}}}) \right ]-\left [\phi({\hat{\vec{x^*}}})^T\Sigma_p\Phi^T \right ](\bar{K}+\sigma^2I)^{-1}\left [\Phi\Sigma_p\phi({\hat{\vec{x^*}}}) \right ]
  • 可见均值和方差中出现的\bar{K}\left [ \cdot \right ]都为\phi(\cdot)\Sigma_p\phi(\cdot)的形式,则设核函数:K(x,x')=\phi(x)\Sigma_p\phi(x'),则可以带入均值和方差:
    • 均值:\vec{\mu}_{\phi}=\left [ K(\hat{\vec{x^*}},\hat{\vec{x}}_1), K(\hat{\vec{x^*}},\hat{\vec{x}}_2),\cdots, K(\hat{\vec{x^*}},\hat{\vec{x}}_N) \right ](\bar{K}+\sigma^{2}I)^{-1}Y,其中\bar{K}=\begin{bmatrix} K(\hat{\vec{x}}_1,\hat{\vec{x}}_1) & K(\hat{\vec{x}}_1,\hat{\vec{x}}_2) & \cdots & K(\hat{\vec{x}}_1,\hat{\vec{x}}_N)\\ K(\hat{\vec{x}}_2,\hat{\vec{x}}_1) & K(\hat{\vec{x}}_2,\hat{\vec{x}}_2) & \cdots & K(\hat{\vec{x}}_2,\hat{\vec{x}}_N)\\ \vdots & \vdots & \ddots & \vdots \\ K(\hat{\vec{x}}_N,\hat{\vec{x}}_1) & K(\hat{\vec{x}}_N,\hat{\vec{x}}_2) & \cdots & K(\hat{\vec{x}}_N,\hat{\vec{x}}_N)\end{bmatrix}
    • 方差:\Sigma_{\phi}=K(\hat{\vec{x^*}},\hat{\vec{x^*}})-\left [ K(\hat{\vec{x^*}},\hat{\vec{x}}_1), K(\hat{\vec{x^*}},\hat{\vec{x}}_2),\cdots, K(\hat{\vec{x^*}},\hat{\vec{x}}_N) \right ](\bar{K}+\sigma^{2}I)^{-1}\begin{bmatrix}K(\hat{\vec{x^*}},\hat{\vec{x}}_1) \\ K(\hat{\vec{x^*}},\hat{\vec{x}}_2) \\ \vdots \\ K(\hat{\vec{x^*}},\hat{\vec{x}}_N) \end{bmatrix}
    • 核函数K(x,x')=\phi(x)\Sigma_p\phi(x')=\phi(x)\Sigma_p^{\frac{1}{2}}\Sigma_p^{\frac{1}{2}}\phi(x')=\left (\Sigma_p^{\frac{1}{2}}\phi(x) \right )^T\left (\Sigma_p^{\frac{1}{2}}\phi(x') \right ),化简为了内积的形式,因此简化了上面的计算;
  • 上述过程实际上就是高斯过程回归的推导,可以理解为(贝叶斯线性回归+非线性转换+核技巧),而这是高斯过程回归的一个视角,叫权重空间视角,还有另外一个函数空间视角将在后面阐述,两个视角会得到相同结果;

19.2.2. 从权重空间到函数空间

  •  根据上述推导结果,f(\hat{\vec{x^*}})\sim N(\vec{\mu}_{\phi},\Sigma_{\phi}),在已知样本数据XY的前提下,\vec{\mu}_{\phi},\Sigma_{\phi}都是\hat{\vec{x^*}}的函数,这个很像前面提到的高斯过程 \xi_t\sim GP(m(t),K(s,t)) 的形式,因此可以从高斯过程来理解;
  • 将待预测(待回归)变量设为X^*=\left[ \hat{\vec{x^*}}_1, \hat{\vec{x^*}}_2, \cdots,\hat{\vec{x^*}}_N \right]^T,则可以将f(\vec{x})f(\hat{\vec{x^*}})看做高斯过程,则 y\sim N(m(\hat{\vec{x}}),K(\hat{\vec{x}},\hat{\vec{x}'})+\sigma^2I)y^*\sim N(m(\hat{\vec{x^*}}),K(\hat{\vec{x^*}},\hat{\vec{x^*}'})+\sigma^2I),因此得到联合概率密度\begin{bmatrix} Y\\ f(\hat{\vec{x^*}}) \end{bmatrix}\sim N(\begin{bmatrix} \mu(\hat{\vec{x}})\\ \mu(\hat{\vec{x^*}}) \end{bmatrix},\begin{bmatrix} K(X,X)+\sigma^2I & K(X,X^*)\\ K(X^*,X) & K(X^*,X^*)\end{bmatrix}),再根据多维高斯分布的条件概率密度公式(2.5),求得\vec{\mu}^*=E[f(\hat{\vec{x^*}})]=K(X^*,X)\left [K(X,X)+\sigma^2I \right ]^{-1}\left [Y-\mu(X) \right ]+\mu(X^*)\begin{aligned} \Sigma^*=D[f(\hat{\vec{x^*}})]=K(X^*,X^*)-K(X^*,X)\left [K(X,X)+\sigma^2I \right ]^{-1}K(X,X^*) \end{aligned},与上面推导的权重空间结果相同,但更加简单;

20. 受限玻尔兹曼机(Restricted Bottzmann Machine,RBM)

20.1. 背景

  • 在马尔可夫随机场(无向图的概率图模型)中,基于对图划分的最大图进行因子分解,即P(x)=\frac{1}{Z}\prod_{i=1}^K\Psi_i(x_{c_i}),其中势函数 \Psi_i(x_{c_i}) 的值需要严格大于0,一种有效的方法为设势函数为指数族分布,即 \Psi_i(x_{c_i})=\exp\{-E(x_{c_i})\},则P(x)=\frac{1}{Z}\prod_{i=1}^K\Psi_i(x_{c_i})=\frac{1}{Z}\exp \{ -\sum_{i=1}^KE(x_{c_i}) \}=\frac{1}{Z}\exp \{ -E(x) \},这样的概率分布被称为玻尔兹曼分布(吉布斯分布);
  • 对于一个变量符合玻尔兹曼分布的无向图模型,将其变量分为观测变量和隐变量,即\vec{x}=\left [x_1,\cdots,x_p \right ]^T=\begin{bmatrix} \vec{h}\\ \vec{v} \end{bmatrix}=\left [h_1,\cdots,h_m,v_1,\cdots,v_n \right ]^T,便是玻尔兹曼机(BM),但这样的模型在Inference问题中无法精确推断,利用MCMC进行近似推断又计算量太大,因此需要进行简化,便有了受限玻尔兹曼机(RBM),简化思路为:h 和 v 之间有连接,但各自内部无连接;

20.2. 模型介绍

20.2.1. 模型表示

  • 从图的结构角度来表示模型,可以将能量函数写为以下形式:E(\vec{x})=E(\vec{h}, \vec{v})=-(\vec{h}^TW\vec{v}+\vec{\alpha}^T\vec{v}+\vec{\beta}^T\vec{h}),因此模型的概率密度函数为:P(\vec{x})=\frac{1}{Z}\exp\{\vec{h}^TW\vec{v}\}\cdot \exp\{\vec{\alpha}^T\vec{v}\}\cdot \exp\{\vec{\beta}^T\vec{h}\}

20.2.2. 模型推断

  • 在推断中需要求解的是三个问题:P(\vec{h}|\vec{v})P(\vec{v}|\vec{h}) 和 P(\vec{v}),前两个计算方法类似,因此只推导一个即可;
  • 设 h_{-i}=\{h_j\}_{j\neq i} 为 \vec{h} 中除了 h_i 以外的变量,因为受限玻尔兹曼机的假设,内部节点之间满足条件独立性,即 h_i\perp h_{-i} |\vec{v},因此P(\vec{h}|\vec{v})=\prod_{i=1}^mP(h_i|\vec{v})=\prod_{i=1}^mP(h_i|h_{-i},\vec{v})
  • 由于实际问题中RBM隐变量节点多数是离散变量,这里讨论最简单的0-1分布,其他离散分布形式可以同理计算,因此\begin{aligned} P(h_i=1|h_{-i},\vec{v})=\frac{P(h_i=1,h_{-i},\vec{v})}{P(h_{-i},\vec{v})}= \frac{P(h_i=1,h_{-i},\vec{v})}{P(h_i=1,h_{-i},\vec{v})+P(h_i=0,h_{-i},\vec{v})}\end{aligned}
  • 进一步地,将概率密度函数分解为 h_i 和 h_{-i} 相关的项:\begin{aligned} E(\vec{h}, \vec{v})&=-(\vec{h}^TW\vec{v}+\vec{\alpha}^T\vec{v}+\vec{\beta}^T\vec{h})\\ &=-\left(\sum_{j=1}^m\sum_{k=1}^nh_jw_{jk}v_k+\sum_{k=1}^n\alpha_kv_k+\sum_{j=1}^m\beta_jh_j \right )\\ &=-\left [\left(\sum_{j\neq i}\sum_{k=1}^nh_jw_{jk}v_k+\sum_{k=1}^n\alpha_kv_k+\sum_{j\neq i}\beta_jh_j \right )+\left (\sum_{k=1}^nh_iw_{ik}v_k+\beta_ih_i \right ) \right ]\\ &=\bar{H}(\vec{h},\vec{v})+H(h_i,\vec{v}) \end{aligned},因此\begin{aligned} P(h_i=1,h_{-i},\vec{v})&=\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})+H(h_i=1,\vec{v}) \right\}\\ &=\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})\right\}\cdot\exp\left\{\sum_{k=1}^nw_{ik}v_k+\beta_i \right\}\end{aligned},同理\begin{aligned} P(h_i=0,h_{-i},\vec{v})&=\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})+H(h_i=0,\vec{v}) \right\}\\ &=\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})\right\}\end{aligned},则\begin{aligned} P(h_i=1|h_{-i},\vec{v})&= \frac{\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})\right\}\cdot\exp\left\{\sum_{k=1}^nw_{ik}v_k+\beta_i \right\}}{\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})\right\}\cdot\exp\left\{\sum_{k=1}^nw_{ik}v_k+\beta_i \right\}+\frac{1}{Z}\exp\left\{\bar{H}(\vec{h},\vec{v})\right\}}\\&=\frac{1}{1+\exp\left\{-\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right ) \right\}}\end{aligned}
  • 上述推导得到的概率密度函数其实就是Sigmoid函数 \sigma(x)=\frac{1}{1+e^{-x}},则\begin{aligned} P(h_i=1|\vec{v})=\sigma\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right ) \end{aligned},同理可得\begin{aligned} P(h_i=0|\vec{v})=1-\sigma\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right )=\frac{\exp\left\{-\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right ) \right\}}{1+\exp\left\{-\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right ) \right\}} \end{aligned}进而得\begin{aligned} P(\vec{h}|\vec{v})=\prod_{i=1}^mP(h_i|\vec{v})=\left [ \sigma\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right ) \right ]^q\left [ 1-\sigma\left (\sum_{k=1}^nw_{ik}v_k+\beta_i \right ) \right ]^{(m-q)}\end{aligned},其中 q 为 h_i= 1 的个数;
  • 无论是看RBM的概率图模型,还是概率密度函数与Sigmoid激活函数的相关性,RBM似乎与神经网络都有着很多相近之处,事实上神经网络就是从RBM发展而来,在后面章节中会有所介绍;
  • 为了方便推导P(\vec{v}),我们将参数 W 进行分割,即 W=\begin{bmatrix} \vec{w}_1\\ \vdots\\ \vec{w}_m \end{bmatrix},则总的概率密度函数中的\vec{h}^TW\vec{v}=\begin{bmatrix} h_1, & \cdots, & h_m \end{bmatrix}\begin{bmatrix} \vec{w}_1\\ \vdots\\ \vec{w}_m \end{bmatrix}\vec{v}=\sum_{i=1}^mh_i\vec{w}_i\vec{v}
  • \begin{aligned} P(\vec{v})&=\int_{\vec{h}}P(\vec{v},\vec{h})\\ &=\int_{\vec{h}}\left [\frac{1}{Z}\exp\{\vec{h}^TW\vec{v}\}\cdot \exp\{\vec{\alpha}^T\vec{v}\}\cdot \exp\{\vec{\beta}^T\vec{h}\} \right ]\\ &=\frac{1}{Z}\exp\{\vec{\alpha}^T\vec{v}\}\left [\int_{h_1}\cdots\int_{h_m}\exp\left\{\sum_{i=1}^m\left (h_i\vec{w}_i\vec{v}+\beta_ih_i \right )\right\} \right ] \\ &=\frac{1}{Z}\exp\{\vec{\alpha}^T\vec{v}\}\left [\int_{h_1}\exp\left\{h_1\vec{w}_1\vec{v}+\beta_1h_1 \right\}\cdots\int_{h_m}\exp\left\{h_m\vec{w}_m\vec{v}+\beta_mh_m \right\} \right ] \\ \end{aligned},继续以h_i为0-1分布的离散分布为例,\begin{aligned} \int_{h_i}\exp\left\{h_i\vec{w}_i\vec{v}+\beta_ih_i \right\} =1+\exp\left\{\vec{w}_i\vec{v}+\beta_i \right\}\end{aligned},则\begin{aligned} P(\vec{v})&=\frac{1}{Z}\exp\{\vec{\alpha}^T\vec{v}\}\prod_{i=1}^m\left (1+\exp\left\{\vec{w}_i\vec{v}+\beta_i \right\} \right ) \\ \end{aligned},为了便于化简,令1+\exp\left\{\vec{w}_i\vec{v}+\beta_i \right\}=\exp\{\log(1+\exp\left\{\vec{w}_i\vec{v}+\beta_i \right\})\},则\begin{aligned} P(\vec{v})&=\frac{1}{Z}\exp\{\vec{\alpha}^T\vec{v}\}\prod_{i=1}^m\exp\left\{\log\left (1+\exp\left\{\vec{w}_i\vec{v}+\beta_i \right\} \right )\right\} \\ &=\frac{1}{Z}\exp\left\{\vec{\alpha}^T\vec{v}+\sum_{i=1}^m\log\left (1+\exp\left\{\vec{w}_i\vec{v}+\beta_i \right\} \right )\right\} \\ \end{aligned},其中的累加项非常类似Softplus函数 softplus(x)=\log(1+e^x),该函数图像接近ReLU激活函数(可以看做ReLU的平滑版,相比ReLU没有“死点”,缓解了其梯度消失问题,但计算略微复杂)是最接近脑神经元的激活模型,且一阶导函数为Softmax函数,因此可得\begin{aligned} P(\vec{v})&=\frac{1}{Z}\exp\left\{\vec{\alpha}^T\vec{v}+\sum_{i=1}^mSoftplus\left (\vec{w}_i\vec{v}+\beta_i \right ) \right\} \\ \end{aligned}

25. 近似推断(Approximate Inference)

  • 推断的动机
    • 推断——即为对原因的追溯,具体表现为对隐变量的推断p(Z|X)
    • 同时推断也是模型学习的需要;
  • 推断的困难
    • 在现实很多问题中,精确推断要么无法做到,要么复杂度指数级增高;
      • 对有向图来说,还存在explain away的困难(概率图模型中head2head结构使得本来独立的变量变得相关,其相关关系求解起来很困难);
      • 对无向图来说,存在节点之间的mutual interaction的困难;

32. 变分自编码器(VAE)

32.1. 模型表示

  • VAE是将概率图模型的思路(变分推断)和神经网络相结合的模型;
  • 在VAE之前的生成模型中,主流思想是一种autoregression的方法,即一种补全的思想,例如捕捉图像中像素点之间的相关性,使用一个随机的像素生成一整张图,本质上是用观测变量的一个或多个维度生成全部维度;
  • VAE则是一种LVM(latent variable models),即引入高阶隐变量,基于高阶隐变量的建模p(x|z),使得模型直接从隐变量 z 中采样即可得到生成的数据,更好地学习到显示的数据分布,十分优雅;
  • 要使模型学习到p(x|z),就需要计算后验概率p(z|x)=\frac{p(x|z)p(z)}{p(x)},使其最大化的过程中学习p(x|z),但因为隐变量z的维度很高,因此直接计算p(x)=\int _zp(x,z)dz\int _zp(x|z)p(z)dz非常困难,因此就需要用到变分推断的方法;
  • VAE模型可以认为是K=\infty的高斯混合模型(GMM),而隐变量z设为服从高斯分布z\sim N(\mu_z,\Sigma_z),实际问题中为了简化,可以设z\sim N(0,I)
  • 对于隐变量到观测变量的生成过程,可以将模型参数化,即z|x\sim N(\mu_{\theta}(x), \Sigma_{\theta}(x)),即每个样本 x^{(i)} 都有对应一组参数 \theta^{(i)} 对应一个高斯分布,而该样本的隐变量服从这个唯一的高斯分布,每个高斯分布的参数都可以用神经网络拟合;
  • 因此,VAE模型引入两个假设
    • 预设定隐变量z\sim N(0, I)是个高维标准高斯分布;
    • 对于一组模型参数\thetaz|x\sim N(\mu_{\theta}(x), \Sigma_{\theta}(x))
  • VAE模型结构也由编码器和解码器组成:
    • X\xrightarrow[Encoder]{​{\theta_1}}Z,学习z|x\sim N(\mu_{\theta_1}(x), \Sigma_{\theta_1}(x))
    • Z\xrightarrow[Decoder]{​{\theta_2}}X,学习从随机抽取的隐变量z映射到新生成数据的参数;

32.2. 模型训练

  • 根据EM算法和变分推断中的推导,对模型参数 \theta 进行推断,就是最大化 \log p_{\theta}(X) 的迭代过程,同时也是最大化变分下界ELBO的过程,因此模型目标函数即为ELBO=E_{q_{\phi}(Z|X)}\left[ \log P_{\theta}(X|Z) \right ]-KL(q_{\phi}(Z|X)\ ||\ P_{\theta}(Z))
  • 对ELBO的第二项,基于两个假设,可得P_{\theta}(Z)=N(0,I),没有需要学习的参数,q_{\phi}(z^{(i)}|x^{(i)})= N(\mu_{\theta}(x^{(i)}), \Sigma_{\theta}(x^{(i)})),分别用两个神经网络拟合均值和方差的对数(方差恒大于0,其对数的值域是全体实数,更符合神经网络的训练,无需再用激活函数什么的限定),完成编码器的训练;
  • 对ELBO的第一项,因为一个样本 x^{(i)} 对应一个隐变量 z^{(i)} ,因此z^{(i)}|x^{(i)}\sim N(\mu_{\theta}(x^{(i)}), \Sigma_{\theta}(x^{(i)}))应当是一个尖锐的单峰分布,所以E_{q_{\phi}(Z|X)}\left[ \log P_{\theta}(X|Z) \right ]\approx \log P_{\theta}(X|Z),无需对所有z求期望,简化了无限维z无法求积分的问题,而最大化\log P_{\theta}(X|Z)的过程类似于分类/回归任务,输入为Z输出为X,若X服从离散分布则视为分类任务,使用交叉熵作为损失函数,若X服从连续分布则视为回归任务,使用MSE作为损失函数,完成解码器训练;
  • 但VAE存在一个问题,当通过编码器得到隐变量 z 的分布后,从该分布中随机采样得到一个生成的隐变量样本,送入解码器生成数据,但此时编码器的信息会丢失,即采样得到的这个新的 z 与已知的 X 解耦,失去了联系,不知道 z 是从哪个分布中得到的,同时解码器和编码器各自训练,反向传播断开,使得生成效果会很差;
  • 为了解决上述问题,使用重参数化技巧,即通过编码器得到某个分布N(\mu_{\theta}(x^{(i)}), \Sigma_{\theta}(x^{(i)}))后,放弃原本的随机采样策略,引入一个误差项 \varepsilon,而新生成的隐变量 z=\mu+\varepsilon \cdot \Sigma,如此反向传播被打通,编码器和解码器可以同步训练,生成效果提升; 

33. 流模型

  • 0
    点赞
  • 7
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值