机器学习中的高斯分布
写在最前：
本系列计划介绍机器学习领域一些基本问题的原理，由于本人是贝叶斯学派的拥趸，因此这些问题都会尽力从贝叶斯理论的视角给出解释，包括首先介绍概率论的基础知识，接着用最直接的最大似然方法解释回归和分类问题，进而上升到贝叶斯解释。用高斯过程求解回归和分类问题是一个非常有意思的视角，但受限于实际问题的规模，往往不会这么做。但高斯过程方法可以让我对核方法有更深入的理解，进而思考如何通过稀疏化降低问题的规模，从而满足实际应用需要，支持向量机（SVM）即是这类例子。概率图模型（PGM）是一种结构化的表示问题的方法，更简洁、更直观，我们实际遇到的很多机器学习模型都可以转换为概率图模型，且概率图模型已经发展至拥有完整的理论基础。以下系列文章将先从传统角度介绍隐马尔科夫模型（HMM）和卡尔曼滤波（KF）方法，进而展示二者如何作为概率图模型的子集，从概率图模型的理论推导出同样的结果。本系列文章大部分内容学自《Pattern Recognition and Machine Learning》一书，该书号称机器学习领域的圣经，如果你是机器学习领域的爱好者或从业者，建议深入研读本书。虽然是英文图书，但你只需要最简单的英文能力，以及最基础的数学知识即可理解本书。
以下为本文内容：
由于高斯分布具有良好的数学特性，以及现实生活中遇到的很多连续随机变量都是高斯的，因此，高斯分布成为一种最常见的概率分布。本文介绍联合高斯分布的边缘、条件分布，以及高斯分布的贝叶斯定理，这是线性高斯模型的基础。本文介绍的内容作为工具性结论会在以下系列文章中反复用到。
多维高斯分布一个重要的特性是，如果两组随机变量的联合分布是高斯的，那么在已知一组随机变量的条件下，另一组随机变量也是高斯分布；同样的，任何一组随机变量的边缘分布也是高斯的。用数学符号表达就是，假设\mathbf{x}、\mathbf{y}是两组随机变量，如果p(\mathbf{x}, \mathbf{y})服从联合高斯分布，则p(\mathbf{x} | \mathbf{y})或p(\mathbf{y} | \mathbf{x})也服从高斯分布；且p(\mathbf{x})或p(\mathbf{y})也服从高斯分布。
一个服从高斯分布的D维随机变量\mathbf{x}，它的概率密度函数可以表示为：
\begin{equation} p(\mathbf{x})\sim N(\mathbf{x} | \mathbf{\mu}, \Sigma)=\frac{1}{(2\pi)^{D/2} |\Sigma|^{1/2}}exp\left\{ -\frac{1}{2}(\mathbf{x} - \mathbf{\mu})^T\Sigma^{-1}(\mathbf{x} - \mathbf{\mu}) \right\} \end{equation}
其中，\mathbf{\mu}为均值向量，\Sigma为方差矩阵，D为\mathbf{x}的维数。可以看到，高斯分布概率密度函数的指数项是\mathbf{x}的二次函数，进一步展开为-\frac{1}{2}\left( \mathbf{x}^T\Sigma^{-1} \mathbf{x} - 2 \mathbf{x}^T\Sigma^{-1}\mathbf{\mu} + const \right)，其中const是不必关心的常数。这就告诉我们，如果已知一个随机变量服从高斯分布，那么只需要找到指数项的二次部分和一次部分，就可以知道这个分布的均值和方差。这一技巧将在下面反复用到。
一、联合高斯分布的条件分布
已知一个D维随机变量\mathbf{x}服从高斯分布p(\mathbf{x})\sim N(\mathbf{x} | \mathbf{\mu}, \Sigma)。将\mathbf{x}拆分成不相交的两部分\mathbf{x_a}和\mathbf{x_b}，不妨假设\mathbf{x_a}是\mathbf{x}的前M维变量，\mathbf{x_b}是剩下的D - M维变量，则有
\begin{equation} \mathbf{x} =\begin{pmatrix} \mathbf{x_a} \\ \mathbf{x_b}\end{pmatrix} \end{equation} ，\mathbf{\mu} = \begin{pmatrix} \mathbf{\mu_a} \\ \mathbf{\mu_b}\end{pmatrix}，\Sigma = \begin{pmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix}，其中\Sigma_{ab}^T = \Sigma_{ba}。
“联合高斯分布的条件分布”的问题是，如果已知\mathbf{x_b}，求\mathbf{x_a}的分布，即p(\mathbf{x_a} | \mathbf{x_b})。由于已经知道这是一个高斯分布，也就是要求条件均值\mathbf{\mu}_{a|b}和条件方差\Sigma_{a|b}。
为了保持记号的简洁，很多时候也可以用精度矩阵表示方差矩阵。精度矩阵定义为方差矩阵的逆，即
\Lambda \equiv \Sigma ^{-1} = \begin{pmatrix} \Lambda_{aa} & \Lambda_{ab} \\ \Lambda_{ba} & \Lambda_{bb} \end{pmatrix}，其中\Lambda_{ab}^T = \Lambda_{ba}，
从而，可以将指数项（忽略系数-1 / 2）分块展开得到三部分：
\begin{align} & (\mathbf{x} - \mathbf{\mu})^T\Sigma^{-1}(\mathbf{x} - \mathbf{\mu}) \\ = & \begin{pmatrix} \mathbf{x}_a - \mathbf{\mu}_a \\ \mathbf{x}_b - \mathbf{\mu}_b\end{pmatrix} ^ T \begin{bmatrix} \Lambda_{aa} & \Lambda_{ab} \\ \Lambda_{ba} & \Lambda_{bb} \end{bmatrix} \begin{pmatrix} \mathbf{x}_a - \mathbf{\mu}_a \\ \mathbf{x}_b - \mathbf{\mu}_b\end{pmatrix} \\ = & (\mathbf{x}_a - \mathbf{\mu}_a)^T \Lambda_{aa} (\mathbf{x}_a - \mathbf{\mu}_a) + 2(\mathbf{x}_a - \mathbf{\mu}_a)^T \Lambda_{ab} (\mathbf{x}_b - \mathbf{\mu}_b) + (\mathbf{x}_b - \mathbf{\mu}_b)^T \Lambda_{bb} (\mathbf{x}_b - \mathbf{\mu}_b) \end{align}
把\mathbf{x_a}以外的变量都当做已知量，利用前面的技巧，首先找出\mathbf{x_a}的二次项为\mathbf{x_a}^T \Lambda_{aa} \mathbf{x_a}，因此\Lambda_{a|b} = \Lambda_{aa} = \Sigma_{a|b}^{-1}；再找出\mathbf{x_a}的一次项为-2 \mathbf{x_a}^T \left[ \Lambda_{aa} \mathbf{\mu}_{a} - \Lambda_{ab}(\mathbf{x}_b - \mathbf{\mu}_b) \right]，所以，\mathbf{\mu}_{a|b} = \mathbf{\mu}_a - \Lambda_{aa}^{-1} \Lambda_{ab} (\mathbf{x}_b - \mathbf{\mu}_b)。剩下的都是常数项，可以直接合入分布函数的归一化常数，因此不需要特别关注。
二、联合高斯分布的边缘分布
该问题要求\mathbf{x}_a的边缘分布，即p(\mathbf{x}_a)。由于已经知道这是一个高斯分布，也就是要求均值E\left[ \mathbf{x}_a \right]和方差cov\left[ \mathbf{x}_a \right]。
由贝叶斯定理可知，
p(\mathbf{x}_a) = \frac{p(\mathbf{x}_a, \mathbf{x}_b)}{p(\mathbf{x}_b | \mathbf{x}_a)}，即
\begin{align} \ln{p(\mathbf{x}_a)} & = \ln{p(\mathbf{x}_a, \mathbf{x}_b)} - \ln{p(\mathbf{x}_b | \mathbf{x}_a)} \\ &= -\frac{1}{2} \left[ (\mathbf{x} - \mathbf{\mu})^T\Lambda(\mathbf{x} - \mathbf{\mu}) - (\mathbf{x}_{b} - \mathbf{\mu}_{b|a})^T\Lambda_{b|a}(\mathbf{x}_{b} - \mathbf{\mu}_{b|a}) \right] + const \end{align}
由上一节的结论可知，p(\mathbf{x}_b | \mathbf{x}_a) \sim N(\mathbf{x}_b | \mathbf{\mu}_{b|a}, \Lambda_{b|a})， 其中，\mathbf{\mu}_{b|a} = \mathbf{\mu}_b - \Lambda_{bb}^{-1} \Lambda_{ba} (\mathbf{x}_a - \mathbf{\mu}_a)，\Lambda_{b|a} = \Lambda_{bb} ，因此
\begin{align} & (\mathbf{x}_{b} - \mathbf{\mu}_{b|a})^T\Lambda_{b|a}(\mathbf{x}_{b} - \mathbf{\mu}_{b|a}) \\ =& \left[ \mathbf{x}_{b} - \mathbf{\mu}_b + \Lambda_{bb}^{-1} \Lambda_{ba} (\mathbf{x}_a - \mathbf{\mu}_a) \right] ^T\Lambda_{bb} \left[ \mathbf{x}_{b} - \mathbf{\mu}_b + \Lambda_{bb}^{-1} \Lambda_{ba} (\mathbf{x}_a - \mathbf{\mu}_a) \right] \\ =& (\mathbf{x}_b - \mathbf{\mu}_b)^T \Lambda_{bb} (\mathbf{x}_b - \mathbf{\mu}_b) + 2(\mathbf{x}_b - \mathbf{\mu}_b)^T \Lambda_{ba} (\mathbf{x}_a - \mathbf{\mu}_a) + (\mathbf{x}_a - \mathbf{\mu}_a)^T \Lambda_{ab} \Lambda_{bb}^{-1} \Lambda_{ba} (\mathbf{x}_a - \mathbf{\mu}_a) \end{align}
该展开式的前两项与(\mathbf{x} - \mathbf{\mu})^T\Lambda(\mathbf{x} - \mathbf{\mu})的展开式的后两项相等，因此
\begin{align} lnp(\mathbf{x}_a) & = -\frac{1}{2} (\mathbf{x}_a - \mathbf{\mu}_a)^T (\Lambda_{aa} - \Lambda_{ab} \Lambda_{bb}^{-1} \Lambda_{ba}) (\mathbf{x}_a - \mathbf{\mu}_a) + const \end{align}，
可见，
\begin{align} E\left[ \mathbf{x}_a \right] &= \mathbf{\mu}_a \\ cov\left[ \mathbf{x}_a \right] &= (\Lambda_{aa} - \Lambda_{ab} \Lambda_{bb}^{-1} \Lambda_{ba})^{-1} \end{align}
根据线性代数基本知识，分块矩阵的逆矩阵可表示为：
\begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{pmatrix} ^{-1} = \begin{pmatrix} \mathbf{M} & -\mathbf{M}\mathbf{B}\mathbf{D}^{-1} \\ -\mathbf{D}^{-1} \mathbf{C} \mathbf{M} & \mathbf{D}^{-1} + \mathbf{D}^{-1}\mathbf{C}\mathbf{M}\mathbf{B}\mathbf{D}^{-1}\end{pmatrix}
其中\mathbf{M} = (\mathbf{A} - \mathbf{B} \mathbf{D}^{-1}\mathbf{C}) ^ {-1}，(\mathbf{A} - \mathbf{B} \mathbf{D}^{-1}\mathbf{C})称为矩阵\begin{pmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{pmatrix}的舒尔补（Schur complement）。因此，由\begin{pmatrix} \Lambda_{aa} & \Lambda_{ab} \\ \Lambda_{ba} & \Lambda_{bb} \end{pmatrix} ^{-1} = \begin{pmatrix} \Sigma_{aa} & \Sigma_{ab} \\ \Sigma_{ba} & \Sigma_{bb} \end{pmatrix}可知，
cov\left[ \mathbf{x}_a \right] = (\Lambda_{aa} - \Lambda_{ab} \Lambda_{bb}^{-1} \Lambda_{ba})^{-1} = \Sigma_{aa}。
该问题也可以用另一种方法求解。边缘分布可以通过对联合分布积分求得
p(\mathbf{x}_a) = \int{p(\mathbf{x}_a, \mathbf{x}_b)}d\mathbf{x}_b
将被积项中关于\mathbf{x}_b的部分提出来，由于具有高斯分布的形式，因此可以直接积掉，将剩下关于\mathbf{x}_a的部分写成标准高斯分布函数的形式，从而得到与上述相同的结论。感兴趣的可以用这种方式推导一下。
三、高斯分布的贝叶斯定理
该问题已知一个多维随机变量\mathbf{x}服从高斯分布：
p(\mathbf{x}) \sim N(\mathbf{x} | \mu, \Lambda^{-1})，
另一个多维随机变量\mathbf{y}在已知\mathbf{x}的条件下也服从高斯分布，且其均值为\mathbf{x}的线性变换，即
p(\mathbf{y} | \mathbf{x}) \sim N(\mathbf{y} | A \mathbf{x} + \mathbf{b}, L^{-1})，
其中，A和\mathbf{b}是线性变换参数，L是精度矩阵。此式是线性高斯模型的一个例子，在介绍卡尔曼滤波一节会重点讲到。
该问题要求\mathbf{y}的边缘分布p(\mathbf{y})，以及\mathbf{x}的后验分布p(\mathbf{x} | \mathbf{y})。
为求p(\mathbf{y})和p(\mathbf{x} | \mathbf{y})，可以先求联合分布p(\mathbf{x} , \mathbf{y})，然后就可以利用前两节的结论直接得到结果。为此，定义\mathbf{z} = \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix}，则
\begin{align} \ln{p(\mathbf{z})} &= \ln{p(\mathbf{x})} + \ln{p(\mathbf{y} | \mathbf{x})} \\ &= -\frac{1}{2} \left[ (\mathbf{x} - \mathbf{\mu})^T\Lambda(\mathbf{x} - \mathbf{\mu}) + (\mathbf{y} - A \mathbf{x} - \mathbf{b})^T L (\mathbf{y} - A \mathbf{x} - \mathbf{b}) \right] + const \end{align}
由于是求\mathbf{z}的分布，那么要从上式的方括号中找到关于\mathbf{z}的形如\begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix} ^ T \begin{pmatrix} \mathbf{?}_{xx} & \mathbf{?}_{xy} \\ \mathbf{?}_{yx} & \mathbf{?}_{yy} \end{pmatrix} \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix} 的二次项，以及形如-2 \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix} ^ T \mathbf{?}的一次项。通过观察可以直接写出，关于\mathbf{x}的二次项为\mathbf{x}^T(\Lambda + A^TLA)\mathbf{x}，关于\mathbf{x}、\mathbf{y}的二次项为-\mathbf{x}^T A^TL \mathbf{y}，关于\mathbf{y}、\mathbf{x}的二次项为-\mathbf{y}^T LA \mathbf{x}，关于\mathbf{y}的二次项为\mathbf{y}^T L \mathbf{y}。因此，关于\mathbf{z}的高斯分布的精度矩阵为
R = \begin{pmatrix} \mathbf{?}_{xx} & \mathbf{?}_{xy} \\ \mathbf{?}_{yx} & \mathbf{?}_{yy} \end{pmatrix} = \begin{pmatrix} \Lambda + A^TLA & -A^TL \\ -LA & L\end{pmatrix}，
利用分块矩阵的逆矩阵公式可以写出方差矩阵为
cov[\mathbf{z}] = R^{-1} = \begin{pmatrix} \Lambda^{-1} & \Lambda^{-1} A^T \\ A \Lambda^{-1} & L^{-1} + A\Lambda^{-1}A^T\end{pmatrix}。
上式方括号中，关于\mathbf{x}、\mathbf{y}的一次项为-2(\mathbf{x}^T \Lambda \mu - \mathbf{x}^T A^T L \mathbf{b} + \mathbf{y} ^T L \mathbf{b} ) = -2 \begin{pmatrix} \mathbf{x} \\ \mathbf{y} \end{pmatrix} ^ T \begin{pmatrix} \Lambda \mu - A^T L\mathbf{b} \\ L \mathbf{b} \end{pmatrix}，
从而，关于\mathbf{z}的高斯分布的均值向量为
E\left[ \mathbf{z} \right] = R^{-1} \begin{pmatrix} \Lambda \mu - A^T L\mathbf{b} \\ L \mathbf{b} \end{pmatrix} = \begin{pmatrix} \mathbf{\mu} \\ A\mathbf{\mu} + \mathbf{b}\end{pmatrix}
至此，\mathbf{x}、\mathbf{y}的联合分布已求出。利用本文第二节的结论，\mathbf{y}的均值向量和方差矩阵分别为
\begin{align} E\left[ \mathbf{y} \right] & = A\mathbf{\mu} + \mathbf{b} \\ cov\left[ \mathbf{y} \right] &= L^{-1} + A\Lambda^{-1}A^T \end{align}
利用本文第一节的结论，\mathbf{x}的后验分布的均值向量和方差矩阵分别为
\begin{align} E\left[ \mathbf{x} | \mathbf{y} \right] & = (\Lambda + A^TLA)^{-1} \left[ A^T L (\mathbf{y} - \mathbf{b}) + \Lambda \mathbf{\mu} \right] \\ cov\left[ \mathbf{x} | \mathbf{y} \right] &= (\Lambda + A^TLA)^{-1} \end{align}