六、概率PCA

发布于 2023-07-17 23:38:26 字数 15768 浏览 0 评论 0 收藏 0

定义隐变量 $ MathJax-Element-417 $ ，它属于低维空间（也称作隐空间，即隐变量所在的空间）。假设 $ MathJax-Element-607 $ 的先验分布为高斯分布： $ MathJax-Element-548 $ ，其均值为 $ MathJax-Element-691 $ ，协方差矩阵为 $ MathJax-Element-751 $ 。
定义观测变量 $ MathJax-Element-422 $ ，它属于高维空间。假设条件概率分布 $ MathJax-Element-439 $ 也是高斯分布： $ MathJax-Element-424 $ 。其中：均值是的 $ MathJax-Element-607 $ 线性函数， $ MathJax-Element-426 $ 为权重， $ MathJax-Element-427 $ 为偏置；协方差矩阵为 $ MathJax-Element-435 $ 。
则PPCA 模型生成观测样本的步骤为：
- 首先以概率 $ MathJax-Element-516 $ 生成隐变量 $ MathJax-Element-607 $ 。
- 然后观测样本 $ MathJax-Element-753 $ 由如下规则生成： $ MathJax-Element-432 $ 。
  其中 $ MathJax-Element-445 $ 是一个 $ MathJax-Element-897 $ 维的均值为零、协方差矩阵为 $ MathJax-Element-435 $ 的高斯分布的噪声： $ MathJax-Element-436 $ 。

6.1 参数求解

6.1.1 解析解

可以利用最大似然准则来确定参数 $ MathJax-Element-437 $ 的解析解。
根据边缘概率分布的定义有：
$ p(\mathbf{\vec x})=\int p(\mathbf{\vec x}\mid \mathbf{\vec z}) d\mathbf{\vec z} $
由于 $ MathJax-Element-516 $ 、 $ MathJax-Element-439 $ 均为高斯分布，因此 $ MathJax-Element-440 $ 也是高斯分布。假 $ MathJax-Element-441 $ 的其均值为 $ MathJax-Element-442 $ ，协方差为 $ MathJax-Element-483 $ 。则：
$ \vec\mu^\prime = \mathbb E[\mathbf{\vec x}] =\mathbb E[\mathbf W\mathbf{\vec z}+\vec\mu +\vec\epsilon] = \vec\mu\\ \mathbf C =cov[\mathbf{\vec x}] = \mathbb E[(\mathbf W\mathbf{\vec z}+\vec\mu +\vec\epsilon)(\mathbf W\mathbf{\vec z}+\vec\mu +\vec\epsilon)^T]\\ =\mathbb E[\mathbf W\mathbf{\vec z}\mathbf{\vec z}^T\mathbf W]+\mathbb E[\vec\epsilon \vec\epsilon^T]+\vec\mu\vec\mu^T=\mathbf W\mathbf W^T+\sigma^2\mathbf I+\vec\mu\vec\mu^T $
推导过程中假设 $ MathJax-Element-607 $ 和 $ MathJax-Element-445 $ 是相互独立的随机变量。
因此 $ MathJax-Element-446 $ 。
给定数据集 $ MathJax-Element-778 $ ，则对数似然函数为：
$ \mathcal L = \log p(\mathbb D;\mathbf W,\vec\mu,\sigma^2)= \sum_{i=1}^N \log p(\mathbf{\vec x}_i;\mathbf W,\vec\mu,\sigma^2)\\ =-\frac {Nn}{2} \log(2\pi) -\frac{N}{2}\log |\mathbf C|-\frac 12 \sum_{i=1}^N(\mathbf{\vec x}_i-\vec\mu)^T\mathbf C^{-1}(\mathbf{\vec x}_i-\vec\mu) $
其中 $ MathJax-Element-705 $ 这里表示行列式的值。
求解 $ MathJax-Element-449 $ ，解得 $ MathJax-Element-450 $ 。
对数据集 $ MathJax-Element-778 $ 进行零均值化，即： $ MathJax-Element-545 $ 。则有：
- $ MathJax-Element-453 $ ，因此 $ MathJax-Element-454 $ 。
  其中 $ MathJax-Element-455 $ 。
- 对数似然函数（忽略常数项 $ MathJax-Element-456 $ ）：
  $ \mathcal L= \log p(\mathbb D;\mathbf W,\vec\mu,\sigma^2)=-\frac N2\log|\mathbf C|-\frac 12 \sum_{i=1}^N \mathbf{\vec x}_i^T\mathbf C^{-1}\mathbf{\vec x}_i =-\frac N2\left[\log|\mathbf C|+tr(\mathbf C^{-1}\mathbf S)\right] $
  其中 $ MathJax-Element-457 $ 为协方差矩阵。记：
  $ \mathbf X=(\mathbf{\vec x}_1^T,\cdots,\mathbf{\vec x}_N^T)^T=\begin{bmatrix} \mathbf{\vec x}_1^T\\ \vdots\\ \mathbf{\vec x}_N^T \end{bmatrix}= \begin{bmatrix} x_{1,1}&x_{1,2}&\cdots&x_{1,n}\\ x_{2,1}&x_{2,2}&\cdots&x_{2,n}\\ \vdots&\vdots&\ddots&\vdots\\ x_{N,1}&x_{N,2}&\cdots&x_{N,n}\\ \end{bmatrix} $
  则 $ MathJax-Element-458 $ 。
Tipping and Bishop(1999b) 证明：
- $ MathJax-Element-844 $ 的所有驻点都可以写做： $ MathJax-Element-478 $ 。其中：
  - $ MathJax-Element-461 $ 的列由协方差矩阵 $ MathJax-Element-666 $ 的任意 $ MathJax-Element-586 $ 个特征向量组成。
  - $ MathJax-Element-464 $ 是对角矩阵，其元素是协方差矩阵 $ MathJax-Element-666 $ 对应的 $ MathJax-Element-586 $ 个特征值 $ MathJax-Element-515 $ 。
  - $ MathJax-Element-468 $ 是任意一个正交矩阵。
- 当 $ MathJax-Element-586 $ 个特征向量被选择为前 $ MathJax-Element-586 $ 个最大的特征值对应的特征向量时， $ MathJax-Element-844 $ 取得最大值。其它的所有解都是鞍点。
假定协方差矩阵 $ MathJax-Element-666 $ 的特征值从大到小排列 $ MathJax-Element-473 $ ，对应的 $ MathJax-Element-897 $ 个特征向量为 $ MathJax-Element-475 $ 。
则最大似然准则得到的解析解为：
- $ MathJax-Element-476 $ ，它由前 $ MathJax-Element-586 $ 个特征向量组成。 $ MathJax-Element-478 $ 。
- $ MathJax-Element-479 $ ，它就是与丢弃的维度相关连的平均方差。
$ MathJax-Element-693 $ 是正交矩阵，因此它可以视作 $ MathJax-Element-586 $ 维隐空间的一个旋转矩阵。
根据 $ MathJax-Element-482 $ ，则 $ MathJax-Element-483 $ 与 $ MathJax-Element-693 $ 无关。这表明： $ MathJax-Element-503 $ 在隐空间中具有旋转不变性，因此 $ MathJax-Element-693 $ 可以选任意一个正交矩阵。
- 这代表某种形式的统计不可区分性，因此有多个 $ MathJax-Element-775 $ 都会产生同样的密度函数 $ MathJax-Element-503 $ 。
- 当 $ MathJax-Element-497 $ 时， $ MathJax-Element-490 $ 。此时 $ MathJax-Element-491 $ ，它表示 $ MathJax-Element-775 $ 的列是对 $ MathJax-Element-522 $ 进行缩放，缩放比例为 $ MathJax-Element-494 $ 。
  由于 $ MathJax-Element-495 $ 是正交的，因此 $ MathJax-Element-775 $ 的列也是正交的。
- 当通过解析解直接求解时，可以直接令 $ MathJax-Element-497 $ 。但是当通过共轭梯度法或者EM 算法求解时， $ MathJax-Element-693 $ 可能是任意的，此时 $ MathJax-Element-775 $ 的列不再是正交的。
  如果需要 $ MathJax-Element-775 $ 是正交的，则需要恰当的后处理。
对于任意一个方向向量 $ MathJax-Element-510 $ ，其中 $ MathJax-Element-502 $ ，分布 $ MathJax-Element-503 $ 在该方向上的方差为 $ MathJax-Element-504 $ 。
- 如果 $ MathJax-Element-510 $ 与 $ MathJax-Element-511 $ 正交，即 $ MathJax-Element-510 $ 是被丢弃的特征向量的某个线性组合，则有 $ MathJax-Element-508 $ 。
  因此有 $ MathJax-Element-509 $ 。
- 如果 $ MathJax-Element-510 $ 就是 $ MathJax-Element-511 $ 其中之一，即 $ MathJax-Element-512 $ ，则有 $ MathJax-Element-513 $ 。
  可以看到：沿着特征向量 $ MathJax-Element-518 $ 方向上的方差 $ MathJax-Element-515 $ 由两部分组成：
  - 单位方差的分布 $ MathJax-Element-516 $ 通过线性映射 $ MathJax-Element-775 $ 之后，在 $ MathJax-Element-518 $ 方向上的投影贡献了： $ MathJax-Element-519 $ 。
  - 噪声模型在所有方向上的各向同性的方差贡献了： $ MathJax-Element-523 $ 。
因此：PPCA 正确的描述了数据集 $ MathJax-Element-841 $ 沿着主轴方向（即 $ MathJax-Element-522 $ 方向）的方差，并且用一个单一的均值 $ MathJax-Element-523 $ 近似了所有剩余方向上的方差。
当 $ MathJax-Element-524 $ 时，不存在降维的过程。此时有 $ MathJax-Element-525 $ ， $ MathJax-Element-526 $ 。根据正交矩阵的性质： $ MathJax-Element-527 $ ，以及 $ MathJax-Element-528 $ ，则有： $ MathJax-Element-529 $ 。
由于计算时需要用到 $ MathJax-Element-530 $ ，这涉及到一个 $ MathJax-Element-619 $ 矩阵的求逆。
可以考虑简化为： $ MathJax-Element-532 $ ，其中 $ MathJax-Element-533 $ 。计算复杂度从 $ MathJax-Element-556 $ 降低到了 $ MathJax-Element-535 $ 。

6.1.2 EM算法解

在PPCA 模型中，数据集 $ MathJax-Element-841 $ 中的每个数据点 $ MathJax-Element-917 $ 都对应一个隐变量 $ MathJax-Element-937 $ ，于是可以使用EM 算法来求解模型参数。
实际上PPCA 模型参数已经可以得到精确的解析解，看起来没必要使用EM 算法。但是EM 算法具有下列优势：
- 在高维空间中，使用EM 算法而不是直接计算样本的协方差矩阵可能具有计算上的优势。
- EM 算法的求解步骤也可以推广到因子分析模型中，那里不存在解析解。
- EM 算法可以为缺失值处理提供理论支撑。
观测变量为 $ MathJax-Element-753 $ ，隐变量为 $ MathJax-Element-607 $ ，则完全数据为 $ MathJax-Element-541 $ ，其中 $ MathJax-Element-778 $ 、 $ MathJax-Element-543 $ 。其中对数据集 $ MathJax-Element-778 $ 进行零均值化，即： $ MathJax-Element-545 $ 。
根据后验概率的定义以及高斯分布的性质，后验概率 $ MathJax-Element-571 $ 。
完全数据的对数似然函数为：
$ \log p(\mathbb D,\mathbb Z;\mathbf W,\sigma^2)=\sum_{i=1}^N\log p(\mathbf{\vec x}_i,\mathbf{\vec z}_i)=\sum_{i=1}^N[\log p(\mathbf{\vec x}_i\mid \mathbf{\vec z}_i)+\log p(\mathbf{\vec z}_i)] $
其中： $ MathJax-Element-547 $ 、 $ MathJax-Element-548 $
- E 步：计算期望：
  $ \mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})} \log p(\mathbb D,\mathbb Z;\mathbf W,\sigma^2)=-\sum_{i=1}^N\left[\frac n2\log(2\pi\sigma^2)+\frac 12 tr(\mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i\mathbf{\vec z}_i^T])\\ + \frac{1}{2\sigma^2}||\mathbf{\vec z}_i||_2^2-\frac{1}{\sigma^2}\mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i]^T\mathbf W^T\mathbf{\vec x}_i \\ +\frac{1}{2\sigma^2}tr(\mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i\mathbf{\vec z}_i^T]\mathbf W^T\mathbf W)+\frac d2\log(2\pi) \right] $
  其中： $ MathJax-Element-549 $ 表示计算期望的概率分布为后验概率分布 $ MathJax-Element-550 $ 。假设此时迭代的参数为 $ MathJax-Element-551 $ 、 $ MathJax-Element-552 $ ，则有：
  $ \mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i]=\mathbf M^{-1}_{old}\mathbf W^{T}_{old}\mathbf{\vec x}_i\\ \mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i\mathbf{\vec z}^T_i]=cov_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i]+\mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i]\mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})}[\mathbf{\vec z}_i]^T\\ =\sigma^{2}_{old}\mathbf M^{-1}_{old}+(\mathbf M^{-1}_{old}\mathbf W^{T}_{old}\mathbf{\vec x}_i)(\mathbf M^{-1}_{old}\mathbf W^{T}_{old}\mathbf{\vec x}_i)^T $
- M 步：求最大化：
  $ \mathbf W_{new},\sigma^2_{new}=\arg\max_{\mathbf W,\sigma^2} \mathbb E_{p(\mathbf{\vec z}\mid \mathbf{\vec x})} \log p(\mathbb D,\mathbb Z;\mathbf W,\sigma^2) $
  解得：
  - $ MathJax-Element-553 $ 。
  - $ MathJax-Element-554 $
EM 算法的物理意义：
- E 步涉及到数据点 $ MathJax-Element-917 $ 在隐空间上的正交投影。
- M 步涉及到隐空间的重新估计，使得满足最大似然，其中投影固定。
一个简单的类比：考虑二维空间中的一组小球，令一维隐空间用一个固定的杆表示。现在使用很多个遵守胡克定律（存储的能量正比于弹簧长度的平方）的弹簧将每个小球与杆相连。
- E 步：保持杆固定，让附着的小球沿着杆上下滑动，使得能量最小。这使得每个小球独立地到达对应于它在杆上的正交投影位置。
- M 步：令附着点固定，然后松开杆，让杆到达能量最小的位置。
重复E 步和M 步，直到满足一个收敛准则。
PPCA 的EM 算法的一个好处是大规模数据的计算效率。在高维空间中，EM 算法每次迭代所需的计算量都比传统的PCA 要小得多。
- PPCA 解析解算法中，对协方差矩阵进行特征分解的计算复杂度为 $ MathJax-Element-556 $ 。如果只需要计算前 $ MathJax-Element-586 $ 个特征向量和它们的特征值，则可以使用 $ MathJax-Element-558 $ 复杂度的算法。
  然后计算协方差矩阵本身需要 $ MathJax-Element-563 $ 的计算量，因此不适合大规模数据。
- PPCA的EM 算法没有显式建立协方差矩阵，其计算量最大的步骤是涉及到对数据集求和的操作，计算代价为 $ MathJax-Element-560 $ 。
  对于较大的 $ MathJax-Element-897 $ ，有 $ MathJax-Element-562 $ ，因此与 $ MathJax-Element-563 $ 相比，EM 算法的计算量大大降低。这可以抵消EM 算法需要多次迭代的缺陷。
PPCA 的EM 算法可以用一种在线的形式执行，其中 $ MathJax-Element-917 $ 被读入、处理。然后在处理下一个数据 $ MathJax-Element-856 $ 时丢弃 $ MathJax-Element-917 $
- E 步中需要计算的量（一个d 维向量和一个 $ MathJax-Element-618 $ 的矩阵 ) 可以分别对每个数据点单独计算。
- M 步中需要在数据点上累积求和，这个可以增量完成。
如果 $ MathJax-Element-657 $ 和 $ MathJax-Element-897 $ 都很大，则这种方式很有优势。

6.2 性质

概率PCA (probabilistic PCA:PPCA) 与传统的PCA 相比，有如下优势：
- 概率PCA 能够描述数据集的主要特征，如期望、方差等。
- 概率PCA 使用EM 算法求解。当只需要计算几个主要的特征向量的情况下，计算效率较高，它避免了计算协方差矩阵 $ MathJax-Element-754 $ 。
- 概率PCA 可以推广到混合概率分布，并且利用EM 算法训练。
- 概率PCA 采用似然函数，这使得它可以与其它的概率密度模型进行比较。
- 概率PCA 是一个生成式模型，可以用于生成新样本。
PPCA 考虑的是低维空间到高维空间的映射，这与传统的PCA 观点不同。传统PCA 观点是高维空间到低维空间的映射。
在PPCA 中，如果希望对数据进行降维，即从个高维空间映射到低维空间，则需要通过贝叶斯定理进行逆映射。
根据后验概率分布 $ MathJax-Element-571 $ ，任意一点 $ MathJax-Element-753 $ 在低维空间中的投影均值为 $ MathJax-Element-573 $ 。
- 如果取极限 $ MathJax-Element-574 $ ，则 $ MathJax-Element-575 $ ，这就是标准的PCA 模型。但此时后验协方差为0，概率密度函数变得奇异。
  但是如果使用EM 算法求解，仍然可以得到一个合法的结果。
- 对于 $ MathJax-Element-576 $ 的情况，低维空间中的投影会向着原点方向偏移。
目前在PCA 讨论中，假定低维空间的维度 $ MathJax-Element-586 $ 是给定的，实际上必须根据应用来选择一个合适的值。
- 如果用于数据可视化，则一般选择为 $ MathJax-Element-578 $ 。
- 如果特征值很自然的分成了两组：一组由很小的值组成，另一组由较大的值组成，两组之间有明显的区分。则 $ MathJax-Element-586 $ 选取为较大的特征值的数量。
  实际上这种明显的区分通常无法看到。
- 按照传统PCA 的方式：
  - 通过交叉验证法选取较好的 $ MathJax-Element-586 $ ，这种方式依赖于后续的模型。
  - 从算法原理的角度设置一个阈值，比如 $ MathJax-Element-581 $ ，然后选取使得下式成立的最小的 $ MathJax-Element-582 $ 的值：
    $ \frac{\sum_{i=1}^{d }\lambda_i}{\sum_{i=1}^{n}\lambda_i} \ge t $
    这种方式需要指定阈值 $ MathJax-Element-585 $ ，从而将 $ MathJax-Element-586 $ 的选择转移到 $ MathJax-Element-585 $ 的选择上。
- 基于PPCA 中的最大似然函数，使用交叉验证的方法，求解在验证集上对数似然函数最大的模型来确定维度的值。
  这种方法不依赖于后续的模型，但是缺点是计算量很大。
- 利用贝叶斯PCA 自动寻找合适的 $ MathJax-Element-586 $ 。
贝叶斯PCA ：在 $ MathJax-Element-775 $ 的每列上定义一个独立的高斯先验，每个这样的高斯分布都有一个独立的方差，由超参数 $ MathJax-Element-597 $ 控制。因此：
$ p(\mathbf W;\vec\alpha)=\prod_{i=1}^d\left(\frac{\alpha_i}{2\pi}\right)^{n/2}\exp\left(-\frac 12 \alpha_i\mathbf{\vec w}_i^T\mathbf{\vec w}_i\right) $
其中 $ MathJax-Element-770 $ 是 $ MathJax-Element-775 $ 的第 $ MathJax-Element-874 $ 列。
$ MathJax-Element-592 $ 可以通过最大化 $ MathJax-Element-593 $ 来迭代的求解。最终某个 $ MathJax-Element-597 $ 可能趋向于无穷大，对应的参数向量 $ MathJax-Element-770 $ 趋向于零，后验概率分布变成了原点处的 $ MathJax-Element-596 $ 函数。这就得到一个稀疏解。
这样低维空间的有效维度由有限的 $ MathJax-Element-597 $ 的值确定。通过这种方式，贝叶斯方法自动的在提升数据拟合程度（使用较多的向量 $ MathJax-Element-770 $ ）和减小模型复杂度（压制某些向量 $ MathJax-Element-770 $ ）之间折中。

6.3 因子分析

因子分析：寻找变量之间的公共因子。
如：随着年龄的增长，儿童的身高、体重会发生变化，具有一定的相关性。假设存在一个生长因子同时支配这两个变量，那么因子分析就是从大量的身高、体重数据中寻找该生长因子。
因子分析 Factor Analysis:FA 是一个线性高斯隐变量模型，它与 PPCA 密切相关。
- 因子分析的定义与PPCA 唯一差别是：给定隐变量 $ MathJax-Element-607 $ 的条件下，观测变量 $ MathJax-Element-753 $ 的条件概率分布的协方差矩阵是一个对角矩阵，而不是一个各向同性的协方差矩阵。
  即： $ MathJax-Element-602 $ ，其中 $ MathJax-Element-615 $ 是一个 $ MathJax-Element-619 $ 的对角矩阵。因此也可以认为PPCA 是一种特殊情形的因子分析。
  如果对 $ MathJax-Element-753 $ 进行了零均值化，则 $ MathJax-Element-606 $ 。
- 与 PPCA 模型相同，因子分析模型假设在给定隐变量 $ MathJax-Element-607 $ 的条件下，观测变量 $ MathJax-Element-753 $ 的各分量 $ MathJax-Element-609 $ 是独立的。
在因子分析的文献中， $ MathJax-Element-775 $ 的列描述了观测变量之间的相关性关系，被称作因子载入factor loading 。
$ MathJax-Element-615 $ 的对角元素，表示每个变量的独立噪声方差，被称作唯一性uniqueness 。
观测变量的边缘概率分布为 $ MathJax-Element-612 $ ，其中 $ MathJax-Element-613 $ 。
与PPCA相同，FA模型对于隐空间的选择具有旋转不变性。
可以使用最大似然法来确定因子分析模型中的参数 $ MathJax-Element-775 $ 、 $ MathJax-Element-615 $ 的值。此时 $ MathJax-Element-616 $ 的最大似然解不再具有解析解，因此必须用梯度下降法或者EM 算法迭代求解。
EM 算法的迭代步骤为：
- E 步：用旧的参数求期望：
  $ \mathbb E[\mathbf{\vec z}_i]=\mathbf G\mathbf W^T\mathbf\Psi^{-1}\mathbf{\vec x}_i\\ \mathbb E[\mathbf{\vec z}_i\mathbf{\vec z}^T_i]=cov[\mathbf{\vec z}_i]+\mathbb E[\mathbf{\vec z}_i]\mathbb E[\mathbf{\vec z}_i]^T=\mathbf G+\mathbb E[\mathbf{\vec z}_i]\mathbb E[\mathbf{\vec z}_i]^T $
  其中 $ MathJax-Element-617 $ 。
  这里使用一个 $ MathJax-Element-618 $ 的矩阵求逆表达式，而不是 $ MathJax-Element-619 $ 的表达式。
- M 步：求最大化来获取新的参数。
  $ \mathbf W_{new}\leftarrow \left[\sum_{i=1}^N \mathbf{\vec x}_i\mathbb E[\mathbf{\vec z}_i]^T\right]\left[\sum_{i=1}^N\mathbb E[\mathbf{\vec z}_i\mathbf{\vec z}_i^T]\right]^{-1}\\ \mathbf \Psi_{new}\leftarrow \text{diag}\left[\mathbf S-\mathbf W_{new}\frac 1N\sum_{i=1}^N\mathbb E[\mathbf{\vec z}_i]\mathbf{\vec x}_i^T\right] $
  其中 $ MathJax-Element-620 $ 将所有非对角线上的元素全部设置为零。