机器学习之潜在狄利克雷分配(LDA)变分EM算法及python实现

最新推荐文章于 2024-03-27 12:14:36 发布

董蝈蝈

最新推荐文章于 2024-03-27 12:14:36 发布

阅读量2.3k

点赞数 10

分类专栏：机器学习 NLP 文章标签： python 机器学习算法人工智能

本文链接：https://blog.csdn.net/weixin_41566471/article/details/106649159

版权

机器学习同时被 2 个专栏收录

17 篇文章 3 订阅

订阅专栏

NLP

8 篇文章 2 订阅

订阅专栏

LDA变分EM算法

变分推理(variational inference)是贝叶斯学习中常用的，含有隐变量模型的学习方法。变分推理和上一章节的吉布斯采样不同，吉布斯采样的算法通过随机抽样的方法近似计算模型的后验概率，而变分推理通过解析的方法计算模型的后验概率的近似值。

本章节内容参考李航博士的《统计学习方法》
书中是简化版本的变分EM算法推导，本章节进行完整的推导。
上一章节介绍了LDA的吉布斯采样方法.

1. 变分推理

变分推理的思想如下：
假设模型是联合概率分布 $p (x, z)$ ，其中 $x$ 是观测变量， $z$ 是隐变量，包括参数。目标是学习模型的后验概率 $p (z ∣ x)$ 。但是这个分布复杂无法直接求解，所以考虑用概率分布 $q (z)$ 来近似条件概率分布 $p (z ∣ x)$ ,之后用KL散度 $K L (q (z) ∣ ∣ p (z ∣ x))$ 计算两者之间的相似度， $q (z)$ 称为变分分布。

KL散度可以写成以下形式：

$\begin{aligned} KL(q(z) || p(z|x)) &= \sum_z q(z) \log \frac {q(z)} {p(z|x)} \\ &= \sum_z q(z) \log q(z) - \sum_z q(z) \log (p(z|x)) \\ & = \sum_z q(z) \log q(z) - \sum_z q(z) \log \frac {p(z, x)}{p(x)} \\ & = \sum_z q(z) \log q(z) - \sum_z q(z) \log p(z, x) + \sum_z q(z) \log p(x) \\ & = \sum_z q(z) \log q(z) - \sum_z q(z) \log p(z, x) + \log p(x) \\ & = \log p(x) -\{ E_q[\log p(x,z)] - E_q[\log q(z)] \} \end{aligned}$

因为KL散度是大于等于0，但且仅当两个分布一致时位0，所以可以得到：

$\log p(x) \ge E_q[\log ps(x,z)] - E_q[\log q(z)]$

不等式右断是左端的下界，左端称为证据(evidence)，有端称为证据下界(evidence lower bound, ELBO)，记做：

$L(q) = E_q[\log p(x,z)] - E_q[\log q(z)]$

因为 $x$ 是观测量，所以 $p (x)$ 是定值，所以通过求解证据下界的最大化来求取 $q (z)$ ，因为证据下界越大，越靠近 $p (x)$ ，BL散度越小，说明 $q (z)$ 和 $p (z ∣ x)$ 的分布越接近。

其实证据下界也可以由Jensen不等式得到，如下推导：

$\begin{aligned} \log p(x) & = \log \int_z p(x, z) d z \\ & = \log \int_z p(x, z) \frac {q(z)} {q(z)} d z \\ & = \log \int_z q(z) \frac {p(z,x)} {q(z)} d z \\ & = \log [E_q \frac{p(x,z)}{q(z)}] \\ & \ge E_q[\log \frac {p(x,z)} {q(z)}] \\ &= E_q[\log p(x,z)] - E_q[\log q(z)] \end{aligned}$

通常假设 $q (z)$ 对 $z$ 的所有分量都是相互独立的（实际是条件独立于参数），即满足

$q(z) = q(z_1)q(z_2)...q(z_n)$

这时的变分分布称为平均场(mean filed).

2. LDA的EM算法推导

2.1 证据下界的定义

在这里插入图片描述
如上图所示，实线圈是LDA的模型参数，而虚线圈代表了变分参数。

和上一章节的参数一致，随机变量 $\theta, \varphi$ 的联合概率分布是：

其中 $w$ 是可观测变量， $\theta, \varphi, z$ 是隐变量， $\alpha, \beta$ 是参数。

定义基于平均场的变分分布

$q(\theta, \varphi, z| \gamma, \rho, \eta) = q(\theta| \gamma) q(\varphi|\rho) q(z| \eta)$

其中 $\gamma, \rho$ 分别是 $\theta, \varphi$ 的狄利克雷分布参数， $\eta$ 是 $z$ 的多项分布参数。

由此，可以得到文本的证据下界

$L(\gamma,\eta,\rho,\alpha,\beta) = E_q[\log p(w, z,\theta, \varphi | \alpha, \beta)] - E_q[\log p(\theta, \varphi, z| \gamma, \rho, \eta)]$

将上式展开

变分参数增加了各个隐变量之间的独立性

下面就是依据变分参数和模型参数将上面每一项展开

第一项：

$\begin{aligned} E_q[\log p(\theta|\alpha)] &= E_q[\log \prod_{m=1}^M p(\theta_m|\alpha)] \\ & = E_q[\log \prod_{m=1}^M \frac {\Gamma(\sum_{k=1}^K \alpha_k)} {\prod_{k=1}^K \Gamma(\alpha_k)} \prod_{k=1}^K \theta_{mk}^{\alpha_k-1} ] \\ & = \sum_{m=1}^M E_q[\log \Gamma(\sum_{k=1}^K \alpha_k) - \sum_{k=1}^K \log \Gamma(\alpha_k) + \sum_{k=1}^K (\alpha_k-1) \log \theta_{mk} ] \\ & = \sum_{m=1}^M \{\log \Gamma(\sum_{k=1}^K \alpha_k) - \sum_{k=1}^K \log \Gamma(\alpha_k) + \sum_{k=1}^K (\alpha_k-1)E_q[ \log \theta_{mk}] \} \\ &= \sum_{m=1}^M \{\log \Gamma(\sum_{k=1}^K \alpha_k) - \sum_{k=1}^K \log \Gamma(\alpha_k) + \sum_{k=1}^K (\alpha_k-1)[\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml})]\} \\ \end{aligned}$

其中 $\Psi$ 是digamma函数，是对数伽马函数的一阶导数。

第二项：

$\begin{aligned} E_q[\log p(\varphi|\beta)] &= E_q[\log \prod_{k=1}^K p(\varphi_k|\beta)] \\ & =E_q[\log \prod_{k=1}^K \frac {\Gamma (\sum_{v=1}^V \beta_v)} {\prod_{v=1}^V \Gamma(\beta_v)} \prod_{v=1}^V \varphi_{kv}^{\beta_v-1}] \\ &= \sum_{k=1}^K E_q[\log \Gamma (\sum_{v=1}^V \beta_v) - \sum_{v=1}^V \log \Gamma(\beta_v) + \sum_{v=1}^V (\beta_v -1) \log \varphi_{kv}] \\ &= \sum_{k=1}^K \{ \log \Gamma (\sum_{v=1}^V \beta_v) - \sum_{v=1}^V \log \Gamma(\beta_v) + \sum_{v=1}^V (\beta_v -1) E_q[ \log \varphi_{kv}] \}\\ &= \sum_{k=1}^K \{ \log \Gamma (\sum_{v=1}^V \beta_v) - \sum_{v=1}^V \log \Gamma(\beta_v) + \sum_{v=1}^V (\beta_v -1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \}\\ \end{aligned}$

第三项：

$\begin{aligned} E_q[\log p(z|\theta)] &= E_q[\log \prod_{m=1}^M \prod_{n=1}^{N_m} \prod_{k=1}^K \theta_{mk}^{I(z_{mn} =k)}] \\ & = \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1} E_q[I(z_{mn} =k)] E_q[\log \theta_{mk}] \\ & = \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1} \eta_{mnk} E_q[\log \theta_{mk}] \\ & = \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1} \eta_{mnk} \{\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \} \end{aligned}$

第四项：

$\begin{aligned} E_q[\log p(w|z,\varphi)] &= E_q[\log \prod_{m=1}^M \prod_{n=1}^{N_m} \prod_{v=1}^V \prod_{k=1}^K \varphi_{kv} ^{I(z_{mn}=k, w_{mn}=v)}]\\ & = E_q[ \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{v=1}^V \sum_{k=1}^K I(z_{mn}=k, w_{mn}=v) \log \varphi_{kv}]\\ & = \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{v=1}^V \sum_{k=1}^K I(z_{mn}=k,) I(w_{mn}=v) E_q[\log \varphi_{kv}] \\ &= \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{v=1}^V \sum_{k=1}^K \eta_{nmk} w_{mn}^v [\Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks})] \\ \end{aligned}$

第五项：

$\begin{aligned} E_q[\log q(\theta| \gamma)] &= E_q[\log \prod_{m=1}^M q(\theta_m|\gamma_m)] \\ & = E_q[\log \prod_{m=1}^M \frac {\Gamma(\sum_{k=1}^K \gamma_{mk})} {\prod_{k=1}^K \Gamma(\gamma_{mk})} \prod_{k=1}^K \theta_{mk}^{\gamma_{mk} - 1}] \\ &=\sum_{m=1}^M E_q[\log \Gamma(\sum_{k=1}^K \gamma_{mk}) - \sum_{k=1}^K \log \Gamma (\gamma_{mk}) + \sum_{k=1}^K (\gamma_{mk} - 1) \log \theta_{mk}] \\ &= \sum_{m=1}^M \{ \log \Gamma(\sum_{k=1}^K \gamma_{mk}) - \sum_{k=1}^K \log \Gamma (\gamma_{mk}) + \sum_{k=1}^K (\gamma_{mk} - 1) E_q[\log \theta_{mk}] \}\\ &= \sum_{m=1}^M \{ \log \Gamma(\sum_{k=1}^K \gamma_{mk}) - \sum_{k=1}^K \log \Gamma (\gamma_{mk}) + \sum_{k=1}^K (\gamma_{mk} - 1) [\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml})] \}\\ \end{aligned}$

第六项：

$\begin{aligned} E_q[\log q(\varphi|\rho)] &= E_q[\log \prod_{k=1}^K q(\varphi_k)| \rho_k] \\ & = E_q[\log \prod_{k=1}^K \frac {\Gamma (\sum_{v=1}^V \rho_{kv})} {\prod_{v=1}^V \Gamma(\rho_{kv})} \prod_{v=1}^V \varphi_{kv}^{\rho_{kv}-1}] \\ & = \sum_{k=1}^K E_q[\log \Gamma (\sum_{v=1}^V \rho_{kv}) - \sum_{v=1}^V \log \Gamma(\rho_{kv}) + \sum_{v=1}^V (\rho_{kv}-1) \log \varphi_{kv}] \\ &= \sum_{k=1}^K \{\log \Gamma (\sum_{v=1}^V \rho_{kv}) - \sum_{v=1}^V \log \Gamma(\rho_{kv}) + \sum_{v=1}^V (\rho_{kv}-1) E_q[\log \varphi_{kv}] \} \\ &= \sum_{k=1}^K \{\log \Gamma (\sum_{v=1}^V \rho_{kv}) - \sum_{v=1}^V \log \Gamma(\rho_{kv}) + \sum_{v=1}^V (\rho_{kv}-1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \} \\ \end{aligned}$

第七项：

$\begin{aligned} E_q[\log q(z| \eta)] &= E_q[\log \prod_{m=1}^M \prod_{n=1}^{N_m} q(z_{mn}|\eta_{mn})]\\ & = E_q[\log \prod_{m=1}^M \prod_{n=1}^{N_m} \prod_{k=1}^K \eta_{mnk}^{I(z_{mn}=k)}]\\ & = \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1}^K E_q[I(z_{mn}=k)] E_q[\log \eta_{mnk}] \\ & = \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1}^K \eta_{mnk} \log \eta_{mnk} \end{aligned}$

将上面的公式带入证据下界

$\begin{aligned} L(\gamma,\eta,\rho,\alpha,\beta) &= E_q[\log p(\theta|\alpha)] + E_q[\log p(\varphi|\beta)] + E_q[\log p(z|\theta)] + E_q[\log p(w|z,\varphi)] \\ &- E_q[\log q(\theta| \gamma)] - E_q[\log q(\varphi|\rho)] - E_q[\log q(z| \eta)] \\ & = \sum_{m=1}^M \{\log \Gamma(\sum_{k=1}^K \alpha_k) - \sum_{k=1}^K \log \Gamma(\alpha_k) + \sum_{k=1}^K (\alpha_k-1)[\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml})]\} \\ &+ \sum_{k=1}^K \{ \log \Gamma (\sum_{v=1}^V \beta_v) - \sum_{v=1}^V \log \Gamma(\beta_v) + \sum_{v=1}^V (\beta_v -1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \}\\ & + \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1}^K \eta_{mnk} \{\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \} \\ & + \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{v=1}^V \sum_{k=1}^K \eta_{nmk} w_{mn}^v [\Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks})] \\ & - \sum_{m=1}^M \{ \log \Gamma(\sum_{k=1}^K \gamma_{mk}) - \sum_{k=1}^K \log \Gamma (\gamma_{mk}) + \sum_{k=1}^K (\gamma_{mk} - 1) [\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml})] \}\\ & - \sum_{k=1}^K \{\log \Gamma (\sum_{v=1}^V \rho_{kv}) - \sum_{v=1}^V \log \Gamma(\rho_{kv}) + \sum_{v=1}^V (\rho_{kv}-1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \} \\ & - \sum_{m=1}^M \sum_{n=1}^{N_m} \sum_{k=1}^K \eta_{mnk} \log \eta_{mnk} \end{aligned}$

2.2 变分参数的估计

求 $\gamma_{mk}$ :
先抽取证据下界中含有 $\gamma_{mk}$ 的项

$\begin{aligned} L(\gamma_{mk}) = & (\alpha_k-1)[\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \\ & + \sum_{n=1}^{N_m} \eta_{mnk} \{\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \\ & - \log \Gamma(\sum_{k=1}^K \gamma_{mk}) + \log \Gamma (\gamma_{mk}) - (\gamma_{mk} - 1) [\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml})] \end{aligned}$

求导：

$\begin{aligned} \frac {\partial L(\gamma_{mk})} {\partial \gamma_{mk}} = & (\alpha_k-1)[\Psi'(\gamma_{mk}) - \Psi'(\sum_{l=1}^K \gamma_{ml})] \\ & + \sum_{n=1}^{N_m} \eta_{mnk} \{\Psi'(\gamma_{mk}) - \Psi'(\sum_{l=1}^K \gamma_{ml}) \\ & - \Psi(\sum_{l=1}^K \gamma_{ml})+ \Psi(\gamma_{mk}) \\ &- (\gamma_{mk} - 1) [\Psi'(\gamma_{mk}) - \Psi'(\sum_{l=1}^K \gamma_{ml})] \\ & - \Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \\ & = [(\alpha_k-1) + \sum_{n=1}^{N_m} \eta_{mnk}- (\gamma_{mk} - 1) ] [\Psi'(\gamma_{mk}) - \Psi'(\sum_{l=1}^K \gamma_{ml})] \end{aligned}$

令上式为零，又有 $\Psi'(\gamma_{mk}) - \Psi'(\sum_{l=1}^K \gamma_{ml}) \ne 0$ ，故

$(\alpha_k-1) + \sum_{n=1}^{N_m}\eta_{mnk} - (\gamma_{mk} - 1) = 0$

解得：

$\gamma_{mk} = \alpha_k + \sum_{n=1}^{N_m} \eta_{mnk}$
求 $\rho_{kv}$ :
先抽取证据下界中含有 $\rho_{kv}$ 的项

$\begin{aligned} L(\rho_{kv}) &= (\beta_v -1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \\ & + \sum_{m=1}^M \sum_{n=1}^{N_m} \eta_{nmk} w_{mn}^v [\Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks})] \\ & - \log \Gamma (\sum_{v=1}^V \rho_{kv}) + \log \Gamma(\rho_{kv}) - (\rho_{kv}-1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \\ \end{aligned}$

求导：

$\begin{aligned} \frac {\partial L(\rho_{kv})} {\partial \rho_{kv}} &= (\beta_v -1) [\Psi'(\rho_{kv} )- \Psi'( \sum_{s=1}^V \rho_{ks} ) ] \\ & + \sum_{m=1}^M \sum_{n=1}^{N_m} \eta_{nmk} w_{mn}^v [\Psi'(\rho_{kv}) - \Psi'(\sum_{s=1}^V \rho_{ks})] \\ & - \Psi( \sum_{s=1}^V \rho_{ks} ) + \Psi(\rho_{kv} ) - (\rho_{kv}-1) [\Psi'(\rho_{kv} )- \Psi'( \sum_{s=1}^V \rho_{ks} ) ] \\ & - \Psi(\rho_{kv} )+ \Psi( \sum_{s=1}^V \rho_{ks} ) \\ & = [(\beta_v -1) + \sum_{m=1}^M \sum_{n=1}^{N_m} \eta_{nmk} w_{mn}^v - (\rho_{kv}-1)] [\Psi'(\rho_{kv}) - \Psi'(\sum_{s=1}^V \rho_{ks})] \\ & = 0 \end{aligned}$

因为 $\Psi'(\rho_{kv}) - \Psi'(\sum_{s=1}^V \rho_{ks}) \ne 0$ ,所以得到：

$\rho_{kv} = \beta_v + \sum_{m=1}^M \sum_{n=1}^{N_m} \eta_{nmk} w_{mn}^v$
求 $\eta_{mnk}$
先抽取证据下界中含有 $\eta_{mnk}$ 的项，又因为 $\sum_{k=1}^K \eta_{mnk} = 1$ , 所以有约束的最优化问题拉格朗日函数为

$\begin{aligned} L(\eta_{mnk}) & = \eta_{mnk} \{\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \} \\ & + \sum_{v=1}^V \eta_{nmk} w_{mn}^v [\Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks})] \\ & - \eta_{mnk} \log \eta_{mnk} \\ & + \lambda(\eta_{mnk} - 1) \end{aligned}$

求导

$\begin{aligned} \frac {\partial L(\eta_{mnk})} {\partial \eta_{mnk}} & = \Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml}) \\ & + \sum_{v=1}^V w_{mn}^v \Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks}) \\ & - 1 - \log \eta_{mnk} + \lambda \\ & = 0 \end{aligned}$

则：

$\log \eta_{mnk} = [\Psi(\gamma_{mk}) + \sum_{v=1}^V w_{mn}^v \Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks})) ] + [\lambda -1 - \Psi(\sum_{l=1}^K \gamma_{ml})]$

当 $k$ 变动时，后面一项不影响 $\eta_{mnk}$ ，所以可以写作

$\eta_{mnk} \ltimes \exp(\Psi(\gamma_{mk}) + \sum_{v=1}^V w_{mn}^v \Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks}))$

在具体实现过程中要对 $\eta_{mnk}$ 归一化处理。

2.3 模型参数的估计

求 $\alpha_k$
先抽取证据下界中含有 $\alpha_k$ 的项

$L(\alpha_k) = \sum_{m=1}^M \{\log \Gamma(\sum_{k=1}^K \alpha_k) - \log \Gamma(\alpha_k) + (\alpha_k-1)[\Psi(\gamma_{mk}) - \Psi(\sum_{l=1}^K \gamma_{ml})]\}$

求偏导得

$\frac {\partial L(\alpha_k)} {\alpha_k} = M[\Psi(\sum_{l=1}^K \alpha_l) - \Psi(\alpha_k)] + \sum_{m=1}^M [\Psi (\gamma_{mk} ) - \Psi(\sum_{l=1}^K \gamma_{ml})]$

这个结果是 $L$ 对 $\alpha$ 得梯度 $g(\alpha)$

在对 $\alpha_l$ 求偏导得

$\frac {\partial^2 L} {\partial \alpha_k \partial \alpha_l} = M[\Psi'(\sum_{l=1}^M \alpha_l) - \delta(k, l) \Psi'(\alpha_k)]$

这个结果是Hessian矩阵 $H(\alpha)$ 。
用以下公式迭代，得到参数 $\alpha$ 的估计值

$\alpha_{new} = \alpha_{old} - H(\alpha_{old})^{-1} g(\alpha_{old})$
求 $\beta_v$
先抽取证据下界中含有 $\beta_v$ 的项

$L(\beta_v) = \sum_{k=1}^K \{ \log \Gamma (\sum_{s=1}^V \beta_s) - \log \Gamma(\beta_v) + (\beta_v -1) [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ] \}$

同样的

$g(\beta) = \frac {\partial L} {\partial \beta_k} = K[\Psi(\sum_{s=1}^V \beta_s) - \Psi(\beta_v)] + \sum_{k=1}^K [\Psi(\rho_{kv} )- \Psi( \sum_{s=1}^V \rho_{ks} ) ]$

$H(\beta) = \frac {\partial L} {\partial \beta_v \partial \beta_s} = K[\Psi'(\sum_{s=1}^V \beta_s) - \delta(v,s) \Psi'(\beta_v)]$

用以下公式迭代：

$\beta_{new} = \beta_{old} - H(\beta_{old})^{-1}g(\beta_{old})$

变分EM算法
E步：求解变分参数

$\gamma_{mk} = \alpha_k + \sum_{n=1}^{N_m} \eta_{mnk}$

$\rho_{kv} = \beta_v + \sum_{m=1}^M \sum_{n=1}^{N_m} \eta_{nmk} w_{mn}^v$

$\eta_{mnk} \ltimes \exp(\Psi(\gamma_{mk}) +\sum_{v=1}^V w_{mn}^v \Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks}))$

M步：求模型参数

$\alpha_{new} = \alpha_{old} - H(\alpha_{old})^{-1} g(\alpha_{old})$

$\beta_{new} = \beta_{old} - H(\beta_{old})^{-1}g(\beta_{old})$

3. LDA的python实现

代码实现变量名称与理论推导参数的对应关系
:param gamma: $\gamma$
:param rho: $\rho$
:param eta: $\eta$
:param alpha: $\alpha$
:param beta: $\beta$

import numpy as np
from scipy.special import digamma, polygamma, gammaln
from copy import deepcopy


class LDAEM(object):
    def __init__(self, K=3):
        """
        变分EM算法实现LDA模型
        :param K: 主题个数，默认值为3
        :param M: 文档个数
        :param Nm: 第m个文档的单词个数
        :param N: 文本单词个数列表，(M,)
        :param V: 单词集合的个数
        :param tockens: 单词tockens
        :param gamma: 变分参数，对应文档主题，(M, K)
        :param rho: 变分参数，对应单词主题，(K, V)
        :param eta: 变分参数，对应文本中单词的主题(M, Nm, K)
        :param alpha: 模型参数文档主题theta的参数，(K,)
        :param beta: 模型参数单词主题varphi的参数，(V,)
        """
        self.K = K

        self.M = None
        self.V = None
        self.tockens = None
        self.N = None

        self.params = {
            'gamma': None,
            'rho': None,
            'eta': None,
            'alpha': None,
            'beta': None
        }

    def _init_params(self):
        """
        初始化参数
        :return:
        """
        # 初始化变分参数
        # 约束条件sum_k gamma_mk = 1, 即对于第m个文档主题之和为1
        gamma = np.random.dirichlet(100 * np.ones(self.K), self.M)

        # 约束条件sum_v rho_kv = 1，即对于第k个主题，所有单词之和为1
        rho = np.random.dirichlet(100 * np.ones(self.V), self.K)

        # 约束条件 sum_k eta_mnk=1，即对于第m个文档，第n个单词，文档主题概率之和为1
        eta = np.array([np.random.dirichlet(100 * np.ones(self.K), Nm) for Nm in self.N])

        # 初始化模型参数
        alpha = np.ones(self.K)
        beta = np.ones(self.V)

        self.params = {
            'gamma': gamma,
            'rho': rho,
            'eta': eta,
            'alpha': alpha,
            'beta': beta
        }

    def _update_gamma(self):
        """
         $$\gamma_{mk} = \alpha_k + \sum_{n=1}^{N_m} \eta_{mnk}$$
        :return:
        """
        alpha = deepcopy(self.params['alpha'])
        eta = deepcopy(self.params['eta'])

        # 预定义占位
        gamma = np.zeros((self.M, self.K))

        for m in range(self.M):
            gamma[m] = alpha + np.sum(eta[m], axis=0)

        # 归一化
        gamma = np.array([gamma[:, k] / np.sum(gamma, axis=1) for k in range(self.K)]).T
        return gamma

    def _update_rho(self):
        """
        $$\rho_{kv} = \beta_v + \sum_{m=1}^M \sum_{n=1}^{N_m}  \eta_{nmk} w_{mn}^v $$
        :return:
        """
        beta = deepcopy(self.params['beta'])
        eta = deepcopy(self.params['eta'])

        rho = np.zeros((self.K, self.V))

        for k in range(self.K):
            for v in range(self.V):
                sum_mn = 0
                for m in range(self.M):
                    for n in range(self.N[m]):
                        sum_mn += eta[m][n][k] * (self.texts[m][n] == self.tockens[v])
                rho[k][v] = sum_mn + beta[v]

        # 归一化
        rho = np.array([rho[k] / np.sum(rho[k]) for k in range(self.K)])
        return rho


    def _update_eta(self):
        """
        $$\eta_{mnk} \ltimes \exp(\Psi(\gamma_{mk})
        +\sum_{v=1}^V w_{mn}^v  \Psi(\rho_{kv}) - \Psi(\sum_{s=1}^V \rho_{ks}))$$
        计算中需要归一化
        :return:
        """
        gamma = deepcopy(self.params['gamma'])
        rho = deepcopy(self.params['rho'])

        eta = np.array([np.ones((n, self.K)) for n in self.N])

        for m in range(self.M):
            for n in range(self.N[m]):
                sum_k = 0
                eta_mn_k_list = np.zeros(self.K)
                for k in range(self.K):
                    sum_digamma_rho_k = 0
                    for v in range(self.V):
                        sum_digamma_rho_k += digamma(rho[k][v]) * (self.tockens[v] == self.texts[m][n])
                    a = digamma(gamma[m][k]) + sum_digamma_rho_k + digamma(np.sum(rho[k]))
                    # 为了防止指数增长中的内存溢出，限定一个阈值
                    if a > 20:
                        a = 20
                    sum_k += np.exp(a)
                    eta_mn_k_list[k] = np.exp(a)
                # 将k维度上的数据归一化后赋值
                eta[m][n] = eta_mn_k_list / sum_k

        return eta

    def _E_step(self):
        """
        更新变分参数，为了方便理解，将gamma,rho, eta分别计算
        :return:
        """
        # 更新gamma
        gamma = self._update_gamma()

        # 更新rho
        rho = self._update_rho()

        # 更新eta
        eta = self._update_eta()

        # 将结果写入
        self.params['gamma'] = gamma
        self.params['rho'] = rho
        self.params['eta'] = eta

    def _update_alpha(self, max_iter=1000, tol=0.1):
        """
        $$\alpha_{new} = \alpha_{old} - H(\alpha_{old})^{-1} g(\alpha_{old})$$
        其中g(alpha)是关于alpha的一阶导数
        H(alpha)是关于alpha的Henssian矩阵
        :return:
        """
        alpha = deepcopy(self.params['alpha'])
        gamma = deepcopy(self.params['gamma'])
        for _ in range(max_iter):
            alpha_old = alpha

            # 计算alpha的一阶导数
            # np.tile是将数据扩展
            g = self.M * (digamma(np.sum(alpha)) - digamma(alpha)) + \
                np.sum(
                    digamma(gamma) - np.tile(digamma(np.sum(gamma, axis=1)), (self.K, 1)).T,
                    axis=0
                )

            # 计算Hessen矩阵
            h = -1 * self.M * polygamma(1, alpha)
            z = self.M * polygamma(1, np.sum(alpha))
            c = np.sum(g / h) / (z ** (-1.0) + np.sum(h ** (-1.0)))

            # update alpha
            alpha = alpha - (g - c) / h

            # 设置终止条件
            if np.sqrt(np.mean(np.square(alpha - alpha_old))) < tol:
                break
        return alpha

    def _update_beta(self, max_iter=1000, tol=0.1):
        """
        $$\beta_{new} = \beta_{old} - H(\beta_{old})^{-1}g(\beta_{old})$$
        """
        beta = deepcopy(self.params['beta'])
        rho = deepcopy(self.params['rho'])

        for _ in range(max_iter):
            beta_old = beta

            g = self.K * (digamma(np.sum(beta)) - digamma(beta)) + \
                np.sum(
                    digamma(rho) - np.tile(digamma(np.sum(rho, axis=1)), (self.V, 1)).T,
                    axis=0
                )

            h = -1 * self.K * polygamma(1, beta)
            z = self.K * polygamma(1, np.sum(beta))
            c = np.sum(g / h) / (z ** (-1.0) + np.sum(h ** (-1.0)))

            beta = beta - (g - c) / h

            if np.sqrt(np.mean(np.square(beta - beta_old))) < tol:
                break
        return beta


    def _M_step(self):
        """
        更新模型参数，alpha, beta
        :return:
        """
        # 更新alpha
        alpha = self._update_alpha()

        # 更新beta
        beta = self._update_beta()

        self.params['alpha'] = alpha
        self.params['beta'] = beta

    def fit(self, texts, tokens, max_iter=10):
        """
        训练入口
        :param texts:
        :param tokens:
        :param max_iter:
        :return:
        """
        self.M = len(texts)
        self.tockens = tokens
        self.V = len(tokens)
        self.N = np.array([len(d) for d in texts])

        self.texts = texts

        self._init_params()

        for i in range(max_iter):
            print('iter: ', i + 1)

            self._E_step()
            self._M_step()

董蝈蝈

关注

10
点赞
踩
20

收藏

觉得还不错? 一键收藏
6
评论
机器学习之潜在狄利克雷分配(LDA)变分EM算法及python实现

LDA变分EM算法变分推理(variational inference)是贝叶斯学习中常用的，含有隐变量模型的学习方法。变分推理和上一章节的吉布斯采样不同，吉布斯采样的算法通过随机抽样的方法近似计算模型的后验概率，而变分推理通过解析的方法计算模型的后验概率的近似值。本章节内容参考李航博士的《统计学习方法》书中是简化版本的变分EM算法推导，本章节进行完整的推导。1. 变分推理变分推理的思想如下：假设模型是联合概率分布p(x,z)p(x,z)p(x,z)，其中xxx是观测变量，zzz是隐变量，
复制链接

扫一扫