单指标时间序列异常检测——基于重构概率的变分自编码（VAE）代码实现（详细解释）

admin管理员组

文章数量:1630012

1. 编写目的

不少论文都是基于VAE完成的异常检测，比如 Donut 、Bagel。尽管 Donut 实现的模型很容易通过继承于重写父类方法的方式实现一个 VAE-baseline，并且 Bagel 中自带了一个 VAE-baselina（感兴趣的小伙伴可以前去查看一下源码），但为了简化过程，详细解释 VAE 用于单指标时间序列异常检测的方法，我重新实现了一个简单的 VAE-baselina，并进行了详细的解释，希望可以帮助到需要的小伙伴们。

可以的话，点个赞并发表评论吧 ~ 一定回复 ~~

2. 参考资料

名称	链接
vae + 重构概率 => 异常检测	https://blog.csdn/smileyan9/article/details/109255466
donut 进行了一定的改进的 vae	https://blog.csdn/smileyan9/article/details/112307506
Bagel 条件变分自编码	https://smileyan.blog.csdn/article/details/113463339
tensorflow官网实现->卷积 vae	https://tensorflow.google/tutorials/generative/cvae
龙书作者 vae 源码	https://github/dragen1860/TensorFlow-2.x-Tutorials/blob/master/12-VAE/main.py

版权声明

未经本人 (smile-yan) 允许，不得转发与转载。

3. 源代码

整个项目的源码地址为：https://github/smile-yan/vae-anomaly-detection-for-timeseries，觉得可以的话，顺手点个星星吧，感谢~

3.1 网络结构

图片源地址为龙书作者开源地址 https://github/dragen1860/TensorFlow-2.x-Tutorials/tree/master/12-VAE

我们简化整个过程，写到 VAE 的初始化函数中如下所示：

class VAE(tf.keras.Model):
    def __init__(self, latent_size=4):
        super(VAE, self).__init__()
        # 与输入数据对接
        self.fc1 = tf.keras.layers.Dense(100)
        # fc1 => μ and log σ^2
        self.fc2 = tf.keras.layers.Dense(latent_size)
        self.fc3 = tf.keras.layers.Dense(latent_size)

        # decode  
        self.fc4 = tf.keras.layers.Dense(100)
        # 试图还原原始数据
        self.fc5 = tf.keras.layers.Dense(120)

所以这个部分可以概述为 fc1 与 fc2 构成了 encode 的过程，此过程完成以后，我们可以得到由 μ \mu μ 和 log ⁡ σ 2 \log \sigma^2 logσ2 组成的隐变量。

而 decode 过程则是分两步还原原始数据。

3.2 encode 过程

初始化过程只是初始化我们要用的变量，真正的 编码过程从这个函数开始，这个过程是非常简单的，可以看作盲盒降维操作。

    def encode(self, x):
        """encode过程，返回 μ 和 log σ^2
        :param x: 单窗口数据
        :return:  μ 和 log σ^2
        """
        h = tf.nn.relu(self.fc1(x))
        # mu, log_variance
        return self.fc2(h), self.fc3(h)

3.3 decode 过程

为了重用方便，将 decode 按照是否使用 sigmoid 函数分为两个过程，这个与 https://tensorflow.google/tutorials/generative/cvae 的 decode 函数添加一个参数的效果是一样的。

    def decode_logits(self, z):
        h = tf.nn.relu(self.fc4(z))
        return self.fc5(h)

    def decode(self, z):
        return tf.nn.sigmoid(self.decode_logits(z))

3.4 概率密度方法（Probability Density Function）

decode 过程可以分为两个步骤，如下面右图所示，可以理解为 g ( z ) g(z) g(z) 对应的是一个数据分布，而 g ( z ) g(z) g(z) 以后的就是重构数据，也就是说重构数据可以理解为从这个分布中采样而得到的，当然，深度神经网络采样过程是模糊的，所以这样的对应关系可以理解为采样得到的，至于怎么采样就是神经网络的参数调整过程了。

所以训练目标，就是使得得到的分布 g ( z ) g(z) g(z) 中采样得到的 x ′ x' x′ 尽可能地接近于观测数据 x x x。接下来我们计算重构概率就是基于这个理论：

使用正常数据训练模型后，对于测试数据 x t x_t xt​：

• 如果它是正常数据，那么它在 g ( z ) g(z) g(z) 对应的分布中，分布密度较大；
• 如果它是异常数据，那么它在 g ( z ) g(z) g(z) 对应的分布中，分布密度较小。

所以计算一个测试数据的重构概率，实质就是计算这个 g   ( z ) g\ (z) g (z) 分布中的分布密度。

考虑到没有基础的小伙伴们，现在我们快速需要介绍一下分布密度的计算：

这里是参考维基百科的内容 https://en.wikipedia/wiki/Probability_density_function :

更一般的情况，我们会计算分布密度的对数值，这里做一个简单的公式推导：

log_normal_pdf  ( x ; μ , σ 2 ) = log ⁡ 1 σ 2 π   e − 1 2 ( x − μ σ ) 2 = − log ⁡   ( σ 2 π ) − 1 2 ( x − μ σ ) 2 = − log ⁡ σ − 1 2 log ⁡ 2 π − 1 2 ( x − μ σ ) 2 = − 1 2 ( log ⁡ σ 2 + log ⁡ 2 π + ( x − μ ) 2 σ 2 ) (1) \text{log\_normal\_pdf}\ (x; \mu, \sigma^2) = \log \frac{1}{\sigma\sqrt{2\pi}} \ e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \\ = - \log\ (\sigma \sqrt{2\pi}) - {\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \\ = - \log \sigma - \frac{1}{2} \log 2\pi - {\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \\ = -\frac{1}{2}\bigl(\log \sigma^2 + \log 2\pi + \frac{(x-\mu)^2}{\sigma^2}\bigr) \tag{1} log_normal_pdf (x;μ,σ2)=logσ2π ​1​ e−21​(σx−μ​)2=−log (σ2π ​)−21​(σx−μ​)2=−logσ−21​log2π−21​(σx−μ​)2=−21​(logσ2+log2π+σ2(x−μ)2​)(1)

所以我们可以得到一个 log_normal_pdf 的 python 代码，如下所示：

def log_normal_pdf(sample, mean, log_var, axis=1):
    log2pi = tf.math.log(2. * np.pi)
    return tf.reduce_sum(-.5 * ((sample - mean) ** 2. * tf.exp(-log_var) + log_var + log2pi), axis=axis)

在后面的内容中我们会多次用到这个计算方法。

3.5 重参数化

重参数化过程就是建立 z z z 与 μ \mu μ 、 σ \sigma σ 之间的关系，为了简化这个关系一般都会直接使用最简单的线性关系（同时也要考虑到梯度下降的需要），所以会令

z = μ + ε ⋅ σ (2) z = \mu + \varepsilon \cdot \sigma \tag{2} z=μ+ε⋅σ(2)

这个对应的python 代码实现非常简单：

def reparameterize(mu, log_var):
    """重参数化，计算隐变量 z = μ + ε ⋅ σ
    :param mu:  均值
    :param log_var: 方差的 log 值
    :return: 隐变量 z
    """
    std = tf.exp(log_var * 0.5)
    eps = tf.random.normal(std.shape)

    return mu + eps * std

3.6 VAE 的损失函数 ELBO 的计算

在 《VAE 模型基本原理简单介绍》 中我们已经比较了解了VAE的面貌，VAE 的损失函数 ELBO 的计算方法如下公式：

log ⁡ p ( x ) ≥ ELBO = E q ( z ∣ x ) [ log ⁡ p ( x , z ) q ( z ∣ x ) ] (3) \log p(x) \ge \text{ELBO} = \mathbb{E}_{q(z|x)}\left[\log \frac{p(x, z)}{q(z|x)}\right] \tag{3} logp(x)≥ELBO=Eq(z∣x)​[logq(z∣x)p(x,z)​](3)

为了方便我们利用对数函数的性质拆解成三个式子的和：

ELBO = E q ( z ∣ x ) [ log ⁡ p ( x ∣ z ) + log ⁡ p ( z ) − log ⁡ q ( z ∣ x ) ] (4) \text{ELBO} = \mathbb{E}_{q(z|x)}\left[ \log p(x| z) + \log p(z) - \log q(z|x) \right] \tag{4} ELBO=Eq(z∣x)​[logp(x∣z)+logp(z)−logq(z∣x)](4)

像这种条件概率计算一般可以使用蒙特卡洛方法进行求解，也就转换成 从特定分布 q ( z ∣ x ) q(z|x) q(z∣x) 中采样得到 z z z ，然后计算一下式子的值：
log ⁡ p ( x ∣ z ) + log ⁡ p ( z ) − log ⁡ q ( z ∣ x ) (5) \log p(x| z) + \log p(z) - \log q(z|x) \tag{5} logp(x∣z)+logp(z)−logq(z∣x)(5)

现在开始写代码求解这个式子：

def compute_loss(model, x):
    # 使用模型的encode方法对输入数据x进行编码，得到均值向量mean和对数方差向量log_var
    mean, log_var = model.encode(x)
    # 应用变分自编码器(VAE)中的重新参数化技巧(reparameterization trick)，根据mean和log_var生成隐变量z
    z = reparameterize(mean, log_var)
    # 使用模型的decode_logits方法对生成的隐变量z进行解码，得到重构数据的对数似然概率x_logit（logits形式）
    x_logit = model.decode_logits(z)
    # 计算原始输入数据x与重构数据x_logit之间的交叉熵损失(cross entropy loss)
    # 使用TensorFlow内置的sigmoid_cross_entropy_with_logits函数，同时考虑sigmoid激活
    cross_ent = tf.nn.sigmoid_cross_entropy_with_logits(logits=x_logit, labels=x)
    # 对交叉熵损失求和，得到负的对数似然概率(log likelihood) log_p_x_z
    log_p_x_z = -tf.reduce_sum(cross_ent)
    # 计算隐变量z关于标准正态分布（均值为0，方差为1）的对数概率密度函数(log probability density function, log pdf)
    log_p_z = log_normal_pdf(z, 0., 0.)
    # 计算隐变量z关于编码分布（由mean和log_var确定）的对数概率密度函数(log pdf)
    log_q_z_x = log_normal_pdf(z, mean, log_var)
    # 计算总损失，包括重建损失(-log_p_x_z)、KL散度项(log_p_z - log_q_z_x)，
    # 并取平均值（tf.reduce_mean），作为整个变分自编码器的负对数似然损失
    return -tf.reduce_mean(log_p_x_z + log_p_z - log_q_z_x)

因为最终我们需要的是一个double类型的数值，所以最后一步是 tf.reduce_mean。

3.7 训练过程

@tf.function
def train_step(model, x, optimizer):
    with tf.GradientTape() as tape:
        loss = compute_loss(model, x)
    gradients = tape.gradient(loss, model.trainable_variables)
    optimizer.apply_gradients(zip(gradients, model.trainable_variables))

3.8 重构概率计算函数

这个地方一定要参考原论文 《Variational Autoencoder based Anomaly Detection using Reconstruction Probability》，关于重构概率的计算，原文中的描述为：
The reconstruction probability that is calculated here is the Monte Carlo estimate of E q ϕ ( z ∣ x ) [ log ⁡ p θ ( x ∣ z ) ] E_{q_\phi(z|x)}[\log p_\theta (x|z)] Eqϕ​(z∣x)​[logpθ​(x∣z)], the second term of the right hand side of equation (7).

摘录其中的算法截图如下：

为了观看更加方便，我们把 x ( i ) x^{(i)} x(i) 看作一次测试数据，这里的对应是一个窗口的数据，所以整个过程，我们写一下：

1. μ \mu μ, σ \sigma σ = encoder ( x ) \text{encoder}(x) encoder(x)；
2. 进行 L L L 次重参数化，得到 L L L 个 z z z，也可以理解为 L L L 个正态分布 z ∼ N ( μ z ( i ) , σ z ( i ) ) z\sim \mathcal{N}\bigl(\mu_{z^{(i)}},\sigma_{z^{(i)}}\bigr) z∼N(μz(i)​,σz(i)​)；
3. 对于每一个 z ∼ N ( μ z ( i ) , σ z ( i ) ) z\sim \mathcal{N}\bigl(\mu_{z^{(i)}},\sigma_{z^{(i)}}\bigr) z∼N(μz(i)​,σz(i)​)，进行 decode 操作，并计算 每一次重构数据与原始数据（观测数据）之间的差异。换句话说，decode 之后，在添加 sigmoid 之前，计算 log_p_x_z 作为本次数据的重构概率。

def reconstruction_prob(self, x, L=10):
    # 使用模型的encode方法对输入数据x进行编码，得到均值向量mean和对数方差向量log_var
    mean, log_var = self.encode(x)

    # 初始化一个空列表samples_z，用于存储L次采样的隐变量z
    samples_z = []

    # 对于L次采样
    for i in range(L):
        # 应用变分自编码器(VAE)中的重新参数化技巧(reparameterization trick)，根据mean和log_var生成隐变量z
        z = reparameterize(mean, log_var)
        # 将生成的隐变量z添加到samples_z列表中
        samples_z.append(z)

    # 初始化重构概率reconstruction_prob为0.
    reconstruction_prob = 0.

    # 对于采样得到的所有隐变量z
    for z in samples_z:
        # 使用模型的decode_logits方法对生成的隐变量z进行解码，得到重构数据的对数似然概率x_logit（logits形式）
        x_logit = self.decode_logits(z)
        # 计算原始输入数据x与重构数据x_logit之间的交叉熵损失(cross entropy loss)
        # 使用TensorFlow内置的sigmoid_cross_entropy_with_logits函数，同时考虑sigmoid激活
        cross_ent = tf.nn.sigmoid_cross_entropy_with_logits(logits=x_logit, labels=x)
        # 对交叉熵损失求和，得到负的对数似然概率(log likelihood) log_p_x_z
        log_p_x_z = -tf.reduce_sum(cross_ent)
        # 将当前隐变量z对应的log_p_x_z累加到总重构概率reconstruction_prob中
        reconstruction_prob += log_p_x_z

    # 取平均值，返回L次采样得到的重构概率的平均值作为总重构概率
    return reconstruction_prob / L

3.9 数据窗口化处理

给定一个时序数据，不断切片，处理成多个窗口数据。

3.10 寻找最合适的重构概率阈值

异常检测是一个二分类任务，因此需要寻找一个合适的重构概率阈值，与每一个点的重构概率进行比较大小，从而确定是否为异常。

关于如何找到最合适的阈值可以参考我以前的一篇博客 快速求解 best F1-score 以及对应的阈值。

4. 直接入手 demo

4.1 clone 到本地

$ git clone git@github:smile-yan/vae-anomaly-detection-for-timeseries.git

4.2 安装依赖

依赖比较少，主要就是 tensorflow 2.x

$ cd vae-anomaly-detection-for-timeseries
$ pip install -r requirement.txt

4.3 运行 demo

$ python main.py

输出效果如下：

The size of train_value_windows is 12179.
The size of test_value_windows is 5151.
Epoch: 1/2, test set ELBO: -38.2476, train time elapse : 77.35 s, test time elapse : 10.38 s
Epoch: 2/2, test set ELBO: -37.8864, train time elapse : 76.15 s, test time elapse : 10.94 s
Epoch: 1/2, test set ELBO: -38.2476, train time elapse : 77.35 s, test time elapse : 10.38 s
Epoch: 2/2, test set ELBO: -37.8864, train time elapse : 76.15 s, test time elapse : 10.94 s
The best threshold: 0.3978
The best f1-score: 0.4375

4.4 任何疑问欢迎留言

如果查看源码遇到不理解之处，结合本文依然不能理解的话，请在下面评论，一定想办法解决问题。感谢支持~~

如果觉得不错，请务必点个赞吧~~您的支持是我更新的最大动力。

4.5 答疑摘录

【答疑 1】

感谢 thunderdu 的评论支持 ~ 您的评论很有意义，对改善博客质量非常有帮助 ~ 感谢 ~

问题描述：关于重构概率的计算，可能与论文原意有所差异。具体内容为：

重构概率计算公式为：

reconstruction probability(i) = 1 L ∑ l = 1 L p θ ( x ( i ) ∣ μ x ^ ( i , l ) , σ x ^ ( i , l ) ) \text{reconstruction\ probability(i)} =\frac{1}{L} \sum^L_{l=1}{p_\theta (x^{(i)}|\mu_{\hat x^{(i,l)}},\sigma_{\hat x^{(i,l)}}}) reconstruction probability(i)=L1​l=1∑L​pθ​(x(i)∣μx^(i,l)​,σx^(i,l)​)

在 博客1 计算过程为

在 本博客 的计算过程为

关键问题： 本博客重构概率的计算是错误的，或者说是不准确的。

问题解释：本博客中，关于重构概率的计算，是通过 “先重构，再计算差异” 的方法实现的。重构是通过已经训练好的模型的 decode 过程完成，计算差异是通过计算交叉熵而得到。本博客的计算方法是有问题的，不严谨的。 但从 “重构误差” 的字面意思上看，也只是能讲得通的。这里主要是当初参考 TensorFlow-2.x-Tutorials 实现的。

优化方案：推荐使用 博客1 的方法求解重构概率。当然也可以在本代码基础上 fork 后修改，改完以后提一个 merge 请求，我们一起完善这份代码，非常感谢 ~

5. 总结

VAE 用于异常检测的 demo 拖了好久，向各位道歉 ~ 主要是写博客确确实实不能当饭吃，官方也不会打赏啥的，而我只求个点赞，感谢各位的理解与支持。

提前祝各位中秋快乐~

Smileyan
2022.9.6 14:50

本文标签： 序列概率重构异常指标