7、切换线性动态系统中的期望传播与广义期望传播方法

关注

原创最新推荐文章于 2026-06-11 14:24:44 发布 45 阅读 AI 写同款· GEO 优化 ›

切换线性动态系统中的期望传播与广义期望传播方法

1. 引言

许多现实世界的问题可以用扩展的经典线性高斯动态系统模型来描述，这些模型引入了（未观察到的）离散状态指示符。在这类扩展模型中，离散指示符决定了过程在特定时间遵循的转移和观测模型。例如，机动目标跟踪与估计、多目标问题、非高斯干扰、未知模型参数、传感器故障以及不同趋势等问题，都可以在条件高斯状态空间模型框架下进行表述。

然而，尽管扩展模型具有强大的表达能力，但精确估计后验分布却是一个难题。一般来说，精确的滤波、平滑或预测后验分布的复杂度会随着观测数量呈指数增长。即使只需要指示变量的边缘分布，该问题仍然是NP难的。

为了解决这个问题，我们引入了一种确定性近似方案，特别适合于寻找平滑的单时间片和双时间片后验分布。该方案可以看作是对先前提出的假设密度滤波方法的对称反向传播和迭代方案。

2. 符号与问题描述

在切换线性动态系统（SLDS）中，假设观测到的 $T$ 个 $d$ 维观测序列 $y_{1:T}$ 是由一阶马尔可夫过程产生的噪声观测。潜在空间由一个 $q$ 维的连续状态 $x_t$ 和一个可以取 $M$ 个值的离散状态 $s_t$ 组成。

条件分布如下：
- 状态转移：$p(x_t|x_{t - 1}, s_{t - 1}, s_t, \theta) \sim N(x_t; A_{s_{t - 1}, s_t}x_{t - 1}, Q_{s_{t - 1}, s_t})$
- 观测模型：$p(y_t|x_t, s_{t - 1}, s_t, \theta) \sim N(y_t; C_{s_{t - 1}, s_t}x_t, R_{s_{t - 1}, s_t})$
- 离散状态转移：$p(s_t|s_{t - 1}, \theta) = \Pi_{s_{t - 1} \to s_t}$

在 $t = 1$ 时，$p(s_1|\theta) = \pi_{s_1}$，$p(x_1|s_1, \theta)$ 是一个具有已知参数的高斯分布。

我们的目标是估计滤波和平滑后的单时间片和双时间片后验分布。然而，由于状态变化历史未知，需要考虑所有可能的指示变量序列 $s_{1:T}$。条件高斯（CG）分布族在边缘化操作下不封闭，因此对所有可能序列求和的结果比简单的CG分布更复杂，条件于 $s_t$ 的后验分布是一个具有 $M^{T - 1}$ 个分量的高斯混合分布。

3. 假设密度滤波

3.1 局部近似

精确滤波后验分布的分量数量会随着每次新观测增加 $M$ 倍。为了避免这种复杂度的指数增长，我们采用一种近似方法：先结合观测证据 $y_t$，然后用一个“最佳拟合”的条件高斯分布来近似得到的关于 $z_t \equiv {s_t, x_t}$ 的后验分布。这里的“最佳拟合”定义为最小化原始分布和近似分布之间的Kullback - Leibler（KL）散度的条件高斯分布。

KL散度定义为：
$KL(\hat{p}(z_t)||q(z_t)) \equiv \sum_{z_t} \hat{p}(z_t) \log \frac{\hat{p}(z_t)}{q(z_t)}$

我们引入符号 $\hat{q}(z_t) = Collapse(\hat{p}(z_t)) \equiv \arg\min_{q \in CG} KL(\hat{p}(z_t)||q(z_t))$ 来表示这个近似过程。

3.2 和积算法

通过在每个递归步骤中进行局部投影，近似滤波器的计算复杂度从指数级降低到了线性级。这种近似前向传播过程被称为假设密度滤波，也称为广义伪贝叶斯2（GPB 2）。

假设密度滤波算法如下：

算法7.1 假设密度滤波
1. 初始化：\(\hat{q}_1(z_1) = \alpha_1(z_1) \propto \psi_1(z_1)\)
2. 对于 \(t = 2, 3, \ldots, T\)：
    - 构建双时间片联合分布：\(\hat{p}_t(z_{t - 1}, z_t) \propto \alpha_{t - 1}(z_{t - 1})\psi_t(z_{t - 1}, z_t)\)
    - 边缘化得到单时间片滤波边缘分布：\(\hat{p}_t(z_t) = \sum_{z_{t - 1}} \hat{p}_t(z_{t - 1}, z_t)\)
    - 用KL散度意义下最接近的CG分布近似 \(\hat{p}_t(z_t)\)：\(\hat{q}_t(z_t) = Collapse(\hat{p}_t(z_t))\)
    - 更新前向消息：\(\alpha_t(z_t) = \hat{q}_t(z_t)\)

其中，$\psi_t$ 是由模型方程定义的势函数：
- $\psi_t(z_{t - 1}, z_t) \equiv p(y_t, z_t|z_{t - 1}, \theta) = p(y_t|x_t, s_{t - 1}, s_t, \theta)p(x_t|x_{t - 1}, s_{t - 1}, s_t, \theta)p(s_t|s_{t - 1}, \theta)$
- $\psi_1(z_0, z_1) \equiv p(y_1, z_1|\theta) = p(y_1|x_1, s_1, \theta)p(x_1|s_1, \theta)p(s_1|\theta)$

4. 期望传播

4.1 反向传播

在建立了前向传播的确定性近似之后，自然会考虑类似的反向传播（平滑）过程。之前的一些尝试都在投影到条件高斯族之外引入了额外的近似。

平滑过程依赖于反向消息 $\beta_t(z_t)$，在精确情况下，$\beta_t(z_t) \propto p(y_{t + 1:T}|z_t, \theta)$，使得 $\alpha_t(z_t)\beta_t(z_t) \propto p(z_t|y_{1:T}, \theta)$。

然而，反向消息 $\beta_t$ 通常不能被归一化，因此不能直接通过投影反向消息来定义反向传播。我们提出一种基于期望传播的近似方法，先构建信念，然后对其进行近似，最后从近似后的信念中推导出新的消息。

期望传播算法如下：

算法7.2 切换线性动态系统的期望传播
1. 初始化：\(\alpha_0 = \beta_T = 1\)
2. 前向传播：对于 \(t = 1, 2, \ldots, T\)，\(t' \equiv t\)
    - 构建双时间片信念：\(\hat{p}_t(z_{t - 1}, z_t) \propto \alpha_{t - 1}(z_{t - 1})\psi_t(z_{t - 1}, z_t)\beta_t(z_t)\)
    - 边缘化得到单时间片边缘分布：\(\hat{p}_t(z_{t'}) = \sum_{z_{t''}} \hat{p}_t(z_{t - 1}, z_t)\)，其中 \(t'' \equiv \{t - 1, t\} \setminus t'\)
    - 找到在Kullback - Leibler意义下最接近 \(\hat{p}_{t'}(z_{t'})\) 的近似分布：\(\hat{q}_{t'}(z_{t'}) = Collapse(\hat{p}_{t'}(z_{t'}))\)
    - 通过除法推断新消息：\(\alpha_t(z_t) = \frac{\hat{q}_t(z_t)}{\beta_t(z_t)}\)
3. 反向传播：对于 \(t = T - 1, T - 2, \ldots, 1\)，\(t' \equiv t - 1\)
    - 构建双时间片信念：\(\hat{p}_t(z_{t - 1}, z_t) \propto \alpha_{t - 1}(z_{t - 1})\psi_t(z_{t - 1}, z_t)\beta_t(z_t)\)
    - 边缘化得到单时间片边缘分布：\(\hat{p}_t(z_{t'}) = \sum_{z_{t''}} \hat{p}_t(z_{t - 1}, z_t)\)，其中 \(t'' \equiv \{t - 1, t\} \setminus t'\)
    - 找到在Kullback - Leibler意义下最接近 \(\hat{p}_{t'}(z_{t'})\) 的近似分布：\(\hat{q}_{t'}(z_{t'}) = Collapse(\hat{p}_{t'}(z_{t'}))\)
    - 通过除法推断新消息：\(\beta_{t - 1}(z_{t - 1}) = \frac{\hat{q}_{t - 1}(z_{t - 1})}{\alpha_{t - 1}(z_{t - 1})}\)
4. 可能的迭代：重复前向和反向传播步骤，直到收敛

4.2 迭代

如果在算法7.2的步骤3中不进行近似，即进行精确的边缘化，那么一次前向传播和一次反向传播就足够了。但在当前设置下，由于“collapse”操作不是线性的，前向和反向消息会相互影响。因此，我们可以迭代算法7.2的步骤1到4，以找到尽可能一致的局部近似。

4.3 支持性

算法7.2中需要考虑支持性问题，即新消息的构建是否会导致相关信念仍然可归一化。在第一次前向传播中，这个条件自动满足。但在一般情况下，由于步骤4中的除法操作，更新消息后相邻的双时间片势函数可能不可归一化。如果检测到归一化问题，可以通过阻尼更新消息来解决。

下面是阻尼更新的公式：
$\alpha_{new_t} = \epsilon\alpha_{* t} + (1 - \epsilon)\alpha {old_t}$

5. 自由能最小化

迭代算法7.2的步骤可以看作是对一个目标函数的最小化过程。在精确情况下，寻找平滑的单时间片和双时间片后验分布的任务可以表述为一个最小化问题。

引入变分分布 $p(z_{1:T})$ 得到一个目标函数 $F$，有时被称为自由能：
$- \log Z = \min_p F(p) \equiv \min_p KL(p||p^*) - \log Z$

其中，$p^ $ 是精确后验分布 $p(z_{1:T}|y_{1:T}, \theta)$，最小化是在与 $p^ $ 相同域上的所有适当分布上进行的。

在势函数 $\psi_t$ 的形式下，精确后验分布可以表示为：
$p^*(z_{1:T}|y_{1:T}, \theta) = \frac{1}{Z} \prod_{t = 1}^T \psi_t(z_{t - 1}, z_t)$

将其代入 $F$ 中得到：
$F(p) = E(p) - S(p)$
$E(p) \equiv - \sum_{t = 1}^T \sum_{z_{t - 1}, t} p(z_{t - 1}, z_t) \log \psi_t(z_{t - 1}, z_t)$
$S(p) \equiv - \sum_{z_{1:T}} p(z_{1:T}) \log p(z_{1:T})$

由于精确联合后验分布在链上的形式，我们可以将最小化问题限制在由单时间片和双时间片边缘分布参数化的所有可能链上。

为了简化计算，我们将单时间片边缘分布 $q_t(z_t)$ 限制为条件高斯分布，得到近似自由能：
$F_{EP}(\hat{p}, \hat{q}) \equiv \sum_{t = 1}^T \sum_{z_{t - 1}, t} \hat{p}(z_{t - 1}, z_t) \log \frac{\hat{p}(z_{t - 1}, z_t)}{\psi_t(z_{t - 1}, z_t)} - \sum_{t = 2}^{T - 1} \sum_{z_t} \hat{q}(z_t) \log \hat{q}(z_t)$

同时，一致性要求被弱化，重叠的双时间片边缘分布只需要在期望上一致。

迭代前向和反向传播可以看作是在约束条件下寻找 $F_{EP}$ 最小值的启发式方法。以下定理描述了这种关系：
定理7.1：信念集合 $\hat{p} t(z {t - 1}, t)$ 和 $\hat{q} t(z_t)$ 构成算法7.2的不动点，当且仅当它们是 $F {EP}$ 在适当约束下的零梯度点。

下面是该定理的证明流程：

graph TD
    A[构建拉格朗日函数] --> B[求偏导数]
    B --> C[令偏导数为0]
    C --> D[得到不动点条件]
    D --> E[与算法7.2的不动点条件对比]
    E --> F[得出结论]

这个定理表明，迭代前向和反向传播可以解释为拉格朗日乘子空间中的不动点迭代。

6. 广义期望传播

目标函数 $F_{EP}$ 与Bethe自由能密切相关。Bethe自由能具有相同的形式，但不限制伪边缘分布的参数族，并保持强一致性约束。

引入期望传播作为Bethe自由能的近似，激发了一种广义期望传播算法，类似于Kikuchi对Bethe自由能的扩展。对于切换线性动态系统，这允许我们在近似边缘分布中保留更多分量，并通过反向传播和迭代方案改进滤波器中的局部贪婪近似。

在基本的EP自由能中，最小化是关于双时间片边缘分布（外簇）和它们的重叠（单时间片边缘分布）的信念。在Kikuchi扩展中，外簇更大，最小化是关于外簇、它们的直接重叠、重叠的重叠等的信念。

我们定义对称扩展的外簇，其大小由 $0 \leq \kappa \leq \lfloor \frac{T - 2}{2} \rfloor$ 表示：
$z_{outer}(i) = {s_{i:i + 2(\kappa + 1) - 1}, x_{i + \kappa, i + \kappa + 1}}$
$z_{over}(i) = z_{outer}(i) \cap z_{outer}(i + 1)$

广义期望传播算法如下：

算法7.3 切换线性动态系统的广义期望传播
1. 初始化：\(\alpha_0 = \beta_N = 1\)，消息初始化为条件高斯势函数，使得 \(\tilde{q}(z_{over}(i)) = \alpha_i(z_{over}(i))\beta_i(z_{over}(i))\) 被归一化
2. 前向传播：对于 \(i = 1, 2, \ldots, N - 1\)，\(i' \equiv i\)
    - 构建外簇信念：\(\tilde{p}_i(z_{outer}(i)) = \frac{\alpha_{i - 1}(z_{over}(i - 1))\Psi(i)(z_{outer}(i))\beta_i(z_{over}(i))}{Z_i}\)，其中 \(Z_i = \sum_{z_{outer}(i)} \alpha_{i - 1}(z_{over}(i - 1))\Psi(i)(z_{outer}(i))\beta_i(z_{over}(i))\)
    - 计算单时间片边缘分布：\(\tilde{p}_i(z_{over}(i')) = \sum_{z_{outer}(i) \setminus z_{over}(i')} \tilde{p}_i(z_{outer}(i))\)
    - 找到在Kullback - Leibler意义下最接近 \(\tilde{p}_{i'}(z_{over}(i'))\) 的近似分布：\(\tilde{q}_{i'}(z_{over}(i')) = Collapse(\tilde{p}_{i'}(z_{over}(i')))\)
    - 通过除法推断新消息：\(\alpha_i(z_{over}(i)) = \frac{\tilde{q}_i(z_{over}(i))}{\beta_i(z_{over}(i))}\)
3. 反向传播：对于 \(i = N, N - 1, \ldots, 2\)，\(i' \equiv i - 1\)
    - 构建外簇信念：\(\tilde{p}_i(z_{outer}(i)) = \frac{\alpha_{i - 1}(z_{over}(i - 1))\Psi(i)(z_{outer}(i))\beta_i(z_{over}(i))}{Z_i}\)
    - 计算单时间片边缘分布：\(\tilde{p}_i(z_{over}(i')) = \sum_{z_{outer}(i) \setminus z_{over}(i')} \tilde{p}_i(z_{outer}(i))\)
    - 找到在Kullback - Leibler意义下最接近 \(\tilde{p}_{i'}(z_{over}(i'))\) 的近似分布：\(\tilde{q}_{i'}(z_{over}(i')) = Collapse(\tilde{p}_{i'}(z_{over}(i')))\)
    - 通过除法推断新消息：\(\beta_{i - 1}(z_{over}(i - 1)) = \frac{\tilde{q}_{t - 1}(z_{over}(i - 1))}{\alpha_{i - 1}(z_{over}(i - 1))}\)
4. 迭代：重复前向和反向传播步骤，直到收敛

选择不同的 $\kappa$ 值可以在计算复杂度和近似自由度之间进行权衡。当 $\kappa = 0$ 时，得到与第4节中EP/Bethe自由能等价的结果；当 $\kappa = \lfloor \frac{T - 2}{2} \rfloor$ 时，得到一个强连接树，找到的后验分布是精确的。

下面是不同方法的对比表格：
| 方法 | 复杂度 | 额外近似 |
| ---- | ---- | ---- |
| EP1 | (O(M^2T)) | 无 |
| EP | (O(M^2TI)) | 无 |
| ABP | (O(M^2T)) | 有 |
| Gibbs | (O(KMT)) | 无 |
| RBPS - M | (O(KMT)) | 无 |
| RBPS - 10M | (O(KMT)) | 无 |

综上所述，期望传播和广义期望传播方法为切换线性动态系统的推理提供了一种有效的近似方案。通过迭代和扩展外簇，可以在计算复杂度和近似精度之间取得平衡。但这些方法仍然存在一些问题，如迭代方案的收敛性、数值稳定性以及近似质量的保证等，需要在未来的研究中进一步解决。

切换线性动态系统中的期望传播与广义期望传播方法

7. 替代反向传播方法

7.1 近似反向消息

在算法的平滑过程中，存在两种先前的方法值得探讨，它们与算法7.2有所不同。

第一种方法是替代反向传播（ABP），其前向传播与假设密度滤波相同，但反向传播不同。ABP基于传统的卡尔曼平滑器形式，以对平滑后验分布 $p(x_t|s_t, y_{1:T})$ 和 $p(s_t|y_{1:T})$ 的近似作为递归基础。该方法分别处理离散和连续的潜在状态，这导致了额外的近似。

对于连续潜在状态 $x_t$ ，近似如下：
$p(x_t, x_{t + 1}|s_t, s_{t + 1}, y_{1:T}) \approx p(x_t, x_{t + 1}|s_t, s_{t + 1}, y_{1:t})q(x_{t + 1}|s_{t + 1}, y_{1:T})$

对于离散潜在状态 $s_t$ ，近似如下：
$p(s_t, s_{t + 1}|y_{1:T}) \approx q(s_{t + 1}|y_{1:T})p(s_t|s_{t + 1}, y_{1:t})$

ABP的额外近似源于分别处理离散和连续潜在状态，未充分利用马尔可夫链的性质，而期望传播（EP）则利用了该性质，不需要额外近似。

7.2 部分平滑

第二种方法是文献中的滤波器，它将滤波和平滑估计近似为单峰分布，其前向传播也称为广义伪贝叶斯1（GPB 1）。在平滑过程中，该方法不平滑离散状态 $s_t$ 的概率，仅对连续状态 $x_t$ 使用标准卡尔曼平滑器递归，就好像没有状态切换一样。

下面是这两种替代方法与算法7.2的对比：
| 方法 | 处理方式 | 额外近似 | 平滑离散状态 |
| ---- | ---- | ---- | ---- |
| ABP | 分别处理离散和连续状态 | 有 | 是 |
| 部分平滑 | 近似为单峰分布 | 有 | 否 |
| 算法7.2 | 统一处理 | 无 | 是 |

8. 实验

为了验证所提出方法的有效性，我们进行了一系列实验，比较了期望传播（EP）、替代反向传播（ABP）和最先进的采样技术的性能。

8.1 与精确后验的比较

我们通过从共轭先验中抽取参数生成模型，并生成长度为8的序列。对于每个任务，我们使用以下方法计算近似的单时间片平滑后验分布：
- EP1 ：算法7.2进行一次前向 - 反向传播，计算复杂度为 $O(M^2T)$。
- EP ：算法7.2迭代直到收敛或最多20次迭代，计算复杂度为 $O(M^2TI)$。
- ABP ：使用上述的替代反向传播方法，复杂度为 $O(M^2T)$。
- Gibbs ：使用高效的Gibbs采样器生成1000个样本，计算成本为 $O(KMT)$。
- RBPS - M ：使用Rao - Blackwellised粒子平滑器，粒子数等于状态数，复杂度为 $O(KMT)$。
- RBPS - 10M ：使用Rao - Blackwellised粒子平滑器，粒子数为状态数的10倍。

实验结果表明，算法7.2几乎总是优于传统方法。尽管粒子平滑器和Gibbs采样器在 $K$ 趋于无穷时能得到精确结果，但在有限时间内，算法7.2能以更短的计算时间获得更准确的结果。ABP的额外近似对结果有负面影响。

以下是实验结果的mermaid流程图：

graph LR
    A[生成模型和序列] --> B[选择方法计算近似后验]
    B --> C{方法}
    C -->|EP1| D1[一次前向 - 反向传播]
    C -->|EP| D2[迭代直到收敛]
    C -->|ABP| D3[替代反向传播]
    C -->|Gibbs| D4[Gibbs采样]
    C -->|RBPS - M| D5[粒子数等于状态数的粒子平滑]
    C -->|RBPS - 10M| D6[粒子数为10倍状态数的粒子平滑]
    D1 --> E[比较结果]
    D2 --> E
    D3 --> E
    D4 --> E
    D5 --> E
    D6 --> E
    E --> F[得出结论]

8.2 与Gibbs采样的比较

在长观测序列（长度为100）的实验中，我们将确定性近似（EP）与5次独立的Gibbs采样运行进行比较。结果显示，所有近似结果相对接近，这表明所有方法都能给出合理的结果。

通过分析误差关系，我们得到EP均方误差的有效界：
$E_{EP} \leq 2 (E_{Gibbs} + D_{Gibbs, EP})$

其中，$E_{EP}$ 是EP的均方误差，$E_{Gibbs}$ 是Gibbs采样的均方误差，$D_{Gibbs, EP}$ 是Gibbs采样和EP之间的均方距离。实验表明，确定性近似与Gibbs采样结果的差距与单个Gibbs采样运行之间的差距处于同一数量级。

8.3 更大外簇的影响

为了探究算法7.3中 $\kappa$ 的影响，我们生成长度为10的序列，并使用不同的 $\kappa$ 值（$0, 1, 2, 3, 4$）计算近似的单节点后验分布。结果显示，随着 $\kappa$ 的增加，后验均值的精度系统地提高。当 $\kappa = 4$ 时，理论上能得到精确结果，但由于实现中的舍入误差，实际结果存在一些差异。

9. 讨论

本文讨论的确定性近似方案为条件高斯状态空间模型的推理问题提供了一种有效的解决方法。与精确推理的指数复杂度相比，新的近似方法的计算时间与观测数量呈线性关系。

该方法基于期望传播算法，是对先前假设密度滤波方法的对称反向传播。与其他替代反向传播方法相比，其优势在于前向和反向传播的基本原理相同，且在反向传播中不需要额外的近似，这显著提高了精度。

通过迭代前向和反向传播步骤，可以找到更优的近似。这一过程与一种变体的Bethe自由能相关，迭代的不动点对应于该能量的极值点。

广义期望传播扩展了基本算法，通过保留更多的分量和更新局部近似，提高了近似的精度。实验表明，增加外簇的大小可以系统地提高近似的准确性。

然而，这些方法仍存在一些问题。首先，算法7.2和7.3的迭代方案不能保证收敛，虽然有更复杂的双循环技术，但速度较慢，需要对消息传递算法进行适当调整。其次，算法中涉及条件高斯势的乘法和边缘化操作，可能导致数值问题，目前缺乏一种既保证数值稳定性又避免反向传播中额外近似的实现。最后，对于近似质量的保证仍然是一个挑战，尽管对于一般的切换线性动态系统，常数因子近似仍然是NP难的，但对于特定子类模型或特定观测序列，可能有望给出保证。

综上所述，期望传播和广义期望传播方法为切换线性动态系统的推理提供了有前景的基础，但仍有许多问题需要在未来的研究中解决。