USAD 复现

这篇论文的复现工作中，我的工作量和贡献量都较低，这里只简单总结我参与的复现全过程。

论文：USAD: Unsupervised Anomaly Detection on Multivariate Time Series

Github 仓库地址：USAD

背景

老师已经基于 USAD 方法的官方代码实现了在我们数据集上的实验，但是代码的性能和论文声明的性能差距很大。因为官方代码只有 ROC 曲线的绘制代码，而没有计算 F1 的代码，没有给出异常检测的阈值，而论文的参数附录页也没有明确给出，考虑是阈值设置的问题。

按照相关工作的通行办法，阈值一般是训练集的自身的 loss 最大值或 99% 分位数。但是老师验证过这些常见的设置，效果都不理想。

复现过程

接手任务前，老师告诉我他已经写了一段简单的 brute force 代码来检查最优的阈值，但是这段代码运行的时间太长，没有跑出结果来。

best_result_info = None
best_threshold = None
for index in range(0,y_pred.shape[0]):
    print(index)
    current_threshold = y_pred[index]
    current_y_pred_labels = []
    for index in range (0,y_pred.shape[0]):
        if y_pred[index] >= current_threshold:
            current_y_pred_labels.append(1)
        else:
            current_y_pred_labels.append(0)

    current_result = calc_point2point(np.array(current_y_pred_labels),np.array(y_test))
    if best_result_info == None:
        best_result_info = current_result
        best_threshold = current_threshold
    else:
        if current_result[0] > best_result_info[0]:
            best_result_info = current_result
            best_threshold = current_threshold

简单看一下，这段代码的复杂度是 \(O(n^2)\)，枚举每个元素的阈值，然后计算对应的 F1 值。在已经确定每个样本的检测 loss 的情况下，这种方法一定能获得最优解。

这段代码运行时间过长，首先考虑是数据范围过大的问题，打印一下数据集样本总数，发现大约是 \(4 \times 10^5\) 的样本。按照算法竞赛的经验，一般 CPU 在执行 C++ 代码时大约能达到每秒 \(10^8\) 次运算，Python 代码只能达到 \(10^6\) 运算。同时，考虑到这段代码的复杂度是 \(O(n^2)\)，所以这段代码的运行时间大约是 \(4 \times 10^5 \times 4 \times 10^5 \div 10^6 = 1600\) 秒，大约半小时。

按照老师说法，这段代码至少运行了一天，这不太符合实际的复杂度计算，因此考虑代码中的 calc_point2point 函数复杂度过高，超过了 \(O(n)\)，或者整个代码被 print(index) 这个耗时的 I/O 操作拖慢了运行时间。

calc_point2point 是实验室自己写的方法，源代码所在的库我无法查看，因此先用 sklearn 的 f1_score 函数替换这个函数，同时，把代码中的 range 步长设置为 100，取消 I/O 的 print 操作。

优化效果很显著，总体运行时间马上降低到了 10 分钟左右，得到的最优 F1 是 0.735，论文声明的结果是 0.79，差距较小。

在等待代码运行结果的同时，我去确认论文中关于 F1 检测阈值的表述，论文中是这样描述的：

As not all of the anomaly detection methods used for comparison provide a mechanism to select anomaly thresholds, we tested pos- sible anomaly thresholds for every model and report the results linked to the highest F1 score.

由于并非所有用于比较的异常检测方法都提供了选择异常阈值的机制，我们测试了每个模型的可能异常阈值，并报告与最高 F1 分数相关的结果。

考虑到 USAD 的官方开源代码中也没有涉及到异常阈值的选择，我认为这个阈值就是通过类似的 brute force 方法得到的……

最终方案

在确认了论文的阈值选择方法大概率是 brute force 枚举最优阈值的情况下，我做了一个简单的优化，来保证代码能够在比较合理的时间内运行完毕：

best_result_info = None
best_threshold = None

sorted_y_pred = np.sort(y_pred)
current_threshold = sorted_y_pred[(y_pred.shape[0] * 3) // 4]

cnt = 0
while current_threshold <= maxn:
    print(current_threshold, best_result_info, best_threshold)
    current_y_pred_labels = []
    for i in range(0, y_pred.shape[0]):
        if y_pred[i] >= current_threshold:
            current_y_pred_labels.append(1)
        else:
            current_y_pred_labels.append(0)

    current_result = f1_score(np.array(current_y_pred_labels), np.array(y_test))
    if best_result_info is None:
        best_result_info = current_result
        best_threshold = current_threshold
    else:
        if current_result > best_result_info:
            best_result_info = current_result
            best_threshold = current_threshold
    current_threshold += gap

    cnt = cnt + 1
    if cnt == 100:
        cnt = 0
        print("Now:", current_threshold, "Best F1: ", best_result_info, "best threshold: ", best_threshold)

print("Best F1: ", best_f1)
print("Best threshold: ", best_threshold)

这段代码选择最终 loss 函数排序后的 3/4 位置值作为初始阈值，然后每次增加一个 gap，直到最大值。这样的代码可以在合理的时间内运行完毕。这个 gap 经过实际测试，在 SWaT 上设置为 0.0001，WADI 上设置为 0.1，可以实现性能和开销的 trade-off。

总结

这个复现工作的难度较低，只需要一些简单的复杂度分析和代码优化，就可以获得和论文相似的结果。这个复现工作的价值在于，我意识到并不是已经发表的论文都是优秀的或者说严谨的。

采用 brute force 方法来选择最优阈值明显是一种投机取巧的做法，在实际工程中很难让人接受，而这个工作的效能非常依赖这个阈值选取的效果，所以这个工作实际上只能作为一个小规模数据集上运行的 demo，可能并不适合实际工程应用。