孤立森林:大数据背景下的最佳异常检测算法之一

孤立森林或“iForest”是一个非常漂亮和优雅简单的算法，可以用很少的参数来识别异常。原始的论文对广大的读者来说是容易理解的，并且包含了很少的数学知识。在这篇文章中，我将解释为什么iForest是目前最好的大数据异常检测算法，提供算法的总结，算法的历史，并分享一个代码实现。

为什么iForest是目前最好的大数据异常检测算法

iForest有着基于ROC性能和精度的一流的综合性能。iForest在各种数据集上的性能均优于大多数其他异常值检测（OD）算法。我从Python离群值检测包（PyOD）的作者那里获取了基准数据，并在Excel中应用了行向绿-红渐变条件格式。深绿色表示数据集的最佳算法，深红色表示性能最差的算法：

我们看到，iForest在大多数数据集中均处于领先地位，如我所计算的均值，中位数和标准差行的颜色所示。iForest的相同优异结果也适用于N次精度：

可扩展性。iForest是性能最快的算法。预期在所有数据集中，基于PCA和基于直方图的离群值（HBOS）都更快。k最近邻（KNN）慢得多，并且随着更多的观测值N而扩展得非常厉害。

我已经成功建立了孤立森林，其中包含在集群环境中以分钟为单位的包含100M个观测值和36列的数据集。这样的数据如果使用sk-learn的KNN（）速度上简直无法忍受。

算法要点总结

一下可以认为是10页原始论文的总结，如果不想深入研究，看一下要点就可以了。

1. 大多数其他离群值检测(OD)算法试图构建“正常”实例的配置文件，然后标记不符合该配置文件的实例。iForest通过利用异常的固有特性明确地孤立异常记录:它们的协变量集合具有不寻常的值。
2. 由于计算量大，现有方法仅限于低维数据和小数据大小。举例：尝试对大数据使用sklearn.neighbor.KNeighborsClassifier吗？
3. 另外，iForest具有低开销的特点。细节：外部节点的数量为n，因为每个观测值n都是独立的。内部节点的总数显然为n-1，而节点的总数为2n-1。因此，我们了解了为什么内存需求是有界的并且随n线性增长。
4. 孤立树节点定义：T是无子外部节点或具有一个测试且恰好有两个子节点（Tₗ，Tᵣ）的内部节点。要构建iTree，我们通过随机选择属性q和拆分值p递归地将X划分为：（i）树达到高度限制，（ii）所有观测值都孤立在其自己的外部节点上，或者（iii） 所有数据的所有属性值都相同。
5. 路径长度。观测值x的路径长度h（x）通过从根节点横穿iTree直到横向终止于外部节点的边数x来度量。E（h（x））是来自孤立树集合的h（x）的平均值。可以从平均路径长度E（h（x））得出异常分数s（x，n）：s（x，n）= 2 ^ [-E（h（x））/ c（n ）]。基本上，s和E（h（x））之间存在单调关系（有关详细信息，请参见附录末尾）。我不会涉及术语c（n），所以我可以保持简短，但是对于任何给定的静态数据集来说，它都是常数。
6. 仅要求用户设置两个变量：要构建的树数和子采样大小。作者利用生成的高斯分布数据进行了实验，这些实验表明如何在很少的树和较小的子样本的情况下相对快速地实现平均路径长度的收敛。
7. 小的次抽样(样本的样本)解决了沼泽化和掩蔽问题。对于异常检测而言，输入数据太大而造成了沼泽化和掩蔽。沼泽化是指将“正常”观测结果误认为“异常”观测结果，因为它被异常所包围，而掩蔽则相反。换句话说，当为一棵树提供包含大部分异常的样本时，一个正常的数据点可能看起来异常。作者用x光检查的数据提供了这种现象的例子。
8. 小的子样本允许每个孤立树被特殊化，因为每个子样本包含一组不同的异常或甚至没有异常
9. iForest不依赖于任何距离或基于密度的测量来识别异常，所以它速度快，计算成本低，这就引出了下一个问题
10. 线性时间复杂度，O(n)通俗地说，这意味着运行时间随着输入的大小线性增加。

算法的历史

一个伟大的新想法和更广泛的采纳之间可能有几十年的滞后性。例如，logistic 函数在1845年被发现，在1922年被重新发现，现在被现代数据科学家用于logistic 回归。近几十年来，一个新想法和它被广泛采用之间的滞后时间已经缩短了，但这仍然是一个有争议的很长的时间。iForest于2008年首次共享，直到2018年底才发布具有商业可行性的应用程序!时间表如下:

12/2008 - iForest发布的原始论文

07/2009 - iForest作者最后一次修改他们的代码实现代码

10/2018- h2o团队为R和Python用户提供iForest代码

01/2019 - PyOD发布面向Python用户的离群点检测(OD)工具包代码

08/2019 - LinkedIn工程团队发布Spark/Scala实现iForest代码

代码的实现

由于本文是关于大数据的，所以假设是AWS集群环境。代码脚手架(QA对象的代码，调试等等)大多被省略了。

我发现iForest可以轻松快速地处理750万行和36个特性，按分钟的顺序完成计算。

Python (h2o):

import h2o # h2o automated data cleaning well for my dataset
import pkg_resources
##################################################################
# print packages + versions for debugging/future reproducibility #
##################################################################dists = [d for d in pkg_resources.working_set] # Filter out distributions you don't care about and use.
dists.reverse()
dists
##################################################################
# initialize h2o cluster and load data
##################################################################
h2o.init() # http://docs.h2o.ai/h2o/latest-stable/h2o-r/docs/reference/h2o.init.htmlimport pyarrow.parquet as pq # allow loading of parquet files
import s3fs                  # for working in AWS s3
s3 = s3fs.S3FileSystem()
df = pq.ParquetDataset('s3a://datascience-us-east-1/anyoung/2_processedData/stack_parquetFiles', filesystem=s3).read_pandas().to_pandas()# check input data loaded correctly; pretty print .shape
print('(' + '; '.join(map('{:,.0f}'.format, 
                          df.shape)) + ')')# if you need to sample data
df_samp_5M = df.sample(n=5000000, frac=None, replace=False, weights=None, random_state=123, axis=None)# convert Pandas DataFrame object to h2o DataFrame object
hf = h2o.H2OFrame(df)# drop primary key column
hf = hf.drop('referenceID', axis = 1) # referenceID causes errors in subsequent code# you can omit rows with nas for a first pass
hf_clean = hf.na_omit()# pretty print .shape with thousands comma separator
print('(' + '; '.join(map('{:,.0f}'.format, 
                          hf.shape)) + ')')from h2o.estimators import H2OIsolationForestEstimator
from h2o.estimators 
import H2OIsolationForestEstimator
fullX = ['v1', 
         'v2', 
         'v3'
        ]# split h2o DataFrame into 80/20 train/test
train_hf, valid_hf = hf.split_frame(ratios=[.8], seed=123)# specify iForest estimator models
isolation_model_fullX = H2OIsolationForestEstimator(model_id = "isolation_forest_fullX.hex", seed = 123)
isolation_model_fullX_cv = H2OIsolationForestEstimator(model_id = "isolation_forest_fullX_cv.hex", seed = 123)# train iForest models
isolation_model_fullX.train(training_frame = hf, x = fullX)
isolation_model_fullX_cv.train(training_frame = train_hf, x = fullX)# save models (haven't figured out how to load from s3 w/o permission issues yet)
modelfile = isolation_model_fullX.download_mojo(path="~/", get_genmodel_jar=True)
print("Model saved to " + modelfile)# predict models
predictions_fullX = isolation_model_fullX.predict(hf)# visualize results
predictions_fullX["mean_length"].hist()

如果你的数据具有想要用iForest验证的标签，那么您可以比较正常实例集与异常实例集的分布，并与原始数据集进行进一步的推断。例如，你可以通过原始数据集中不同的特征组合来查看计数，如下所示:

N = df.count()
df[['v1', 'v2', 'id']].groupby(['v1', 'v2']).count() / N
df[['v1', 'v3', 'id']].groupby(['v1', 'v3']).count() / N
...

并与iForest确定的正常/异常实例集进行比较，如下图所示:

##################################################################
# column bind predictions from iForest to the original h2o DataFrame
##################################################################
hf_X_y_fullX = hf.cbind(predictions_fullX)
##################################################################
# Slice using a boolean mask. The output dataset will include rows 
# with column value meeting condition
##################################################################
mask = hf_X_y_fullX["label"] == 0
hf_X_y_fullX_0 = hf_X_y_fullX[mask,:]
mask = hf_X_y_fullX["label"] == 1
hf_X_y_fullX_1 = hf_X_y_fullX[mask,:]
##################################################################
# Filter to only include records that are clearly normal
##################################################################
hf_X_y_fullX_ml7 = hf_X_y_fullX[hf_X_y_fullX['mean_length'] >= 7]
hf_X_y_fullX_0_ml7 = hf_X_y_fullX_1[hf_X_y_fullX_0['mean_length'] >= 7]
hf_X_y_fullX_1_ml7 = hf_X_y_fullX_3[hf_X_y_fullX_1['mean_length'] >= 7]
##################################################################
# Convert to Pandas DataFrame for easier counting/familiarity
##################################################################
hf_X_y_fullX_ml7_df = h2o.as_list(hf_X_y_fullX_ml7, use_pandas = True)
hf_X_y_fullX_0_ml7_df = h2o.as_list(hf_X_y_fullX_0_ml7, use_pandas = True)
hf_X_y_fullX_1_ml7_df = h2o.as_list(hf_X_y_fullX_1_ml7, use_pandas = True)
##################################################################
# Look at counts by combinations of variable levels for inference
##################################################################
hf_X_y_fullX_ml7_df[['v1', 'v2', 'id']].groupby(['v1', 'v2']).count()
hf_X_y_fullX_0_ml7_df = h2o.as_list(hf_X_y_fullX_0_ml7, use_pandas = True)

# Repeat above for anomalous records:##################################################################
# Filter to only include records that are clearly anomalous
##################################################################
hf_X_y_fullX_ml3 = hf_X_y_fullX[hf_X_y_fullX['mean_length'] < 3]
hf_X_y_fullX_0_ml3 = hf_X_y_fullX_1[hf_X_y_fullX_0['mean_length'] < 3]
hf_X_y_fullX_1_ml3 = hf_X_y_fullX_3[hf_X_y_fullX_1['mean_length'] < 3]
##################################################################
# Convert to Pandas DataFrame for easier counting/familiarity
##################################################################
hf_X_y_fullX_ml3_df = h2o.as_list(hf_X_y_fullX_ml3, use_pandas = True)
hf_X_y_fullX_0_ml3_df = h2o.as_list(hf_X_y_fullX_0_ml3, use_pandas = True)
hf_X_y_fullX_1_ml3_df = h2o.as_list(hf_X_y_fullX_1_ml3, use_pandas = True)

完成所有这些，并数据导出到Excel中，以快速生成一些累积分布函数，像这样:

参考文献

F. T. Liu, K. M. Ting, and Z.-H. Zhou. Isolation forest. In: Proceedings of the 8th IEEE International Conference on Data Mining (ICDM’08), Pisa, Italy, 2008, pp.413–422. [code] This paper won the Theoretical/Algorithms Runner-Up Best Paper Award at IEEE ICDM’08
Zhao, Y., Nasrullah, Z. and Li, Z., 2019. PyOD: A Python Toolbox for Scalable Outlier Detection. Journal of machine learning research (JMLR), 20(96), pp.1–7.

附录

算法要点/小结，第5点具体内容:

E(h(x))→c(n)， s→0.5(观测不明显异常)

当E(h(x))→0,s→1(异常)

当E(h(x))→n−1,s→0(不是异常)

所以用这篇论文的话来总结一下上面的内容:

(a)如果实例返回s非常接近1，那么它们肯定是异常的，

(b)如果实例的s远远小于0.5，那么它们被视为正常实例是相当安全的，并且

(c)如果所有实例返回s≈0.5，那么整个样本实际上没有任何明显的异常。

有助于说明异常得分、s和平均路径长度E(h(x))之间关系的图表

作者：Andrew Young

deephub翻译组