Category Archives: Computer Science

强化学习的入门方法

强化学习是人工智能领域一个非常重要的分支,在游戏 AI 领域的应用是无处不在的,包括大家耳熟能详的 AlphaGo 和王者荣耀,其背后都使用了大量的强化学习方法。随着人工智能技术在互联网公司的日渐普及,人们对人工智能的期待也越来越高,也期待着人工智能未来能够做更多的事情。除了各种门禁系统,人脸识别系统在各个公司的应用之外,在游戏领域,强化学习也扮演着举足轻重的作用。

当前的强化学习技术除了使用在游戏 AI 领域之外,也在机器人控制,推荐系统等领域有着很多应用,也取得了不少商业价值。但强化学习的学习还是存在一定的门槛,除了一些在线的 tutorial 和零零散散的资料之外,其实还是需要一本相对完整和优秀的教材来辅助高校的老师和学生,只有这样,强化学习才能够成为一门课进入高校课堂。

上海交大的强化学习团队教师张伟楠,沈键,俞勇合力撰写了一本关于强化学习的,理论与实践相结合的书籍《动手学强化学习》。本书的每个章节除了必要的理论介绍,公式推导之外,还包含了相关的 Python Notebook。Notebook 中包含了算法过程和可运行代码,方便读者阅读相关内容并立刻进行实践。

《动手学强化学习》

《动手学强化学习》是一本相对完整的书籍,读者并不需要知道太多的预备知识就可以直接上手这门课。本书包括三个部分:

  • 第一部分:强化学习基础;
  • 第二部分:强化学习进阶;
  • 第三部分:强化学习前沿。

在第一部分,作者们介绍了强化学习的基本原理,多臂老虎机,马尔可夫决策过程,动态规划算法等内容;在第二部分,作者介绍了 DQN 算法及其改进算法,策略梯度算法,Actor-Critic 算法等诸多算法,让作者能够在第一部分的基础上更进一步的理解强化学习;在第三部分,作者们撰写了前沿的强化学习算法,包括离线强化学习,多智能体强化学习等内容,让读者在之前的基础上迅速进入到强化学习的前沿领域。

本书的另外一个亮点就是 Python Notebook 的代码丰富完善,读者可以一边阅读强化学习的基础知识,一边进行代码的开发工作。基于理论和实战相结合的办法,读者可以尽快掌握强化学习的实战方法。学习这本书需要读者有一定的 Python 开发基础,并且这本书有相应的强化学习视频,方便读者尽快掌握这本书。

Advertisement

如何入门 PyTorch?

PyTorch 是学术界和工业界用来研究和开发深度学习模型的常用工具之一,自 2016 年 9 月 PyTorch 在 GitHub 上开放了 pytorch alpha-1 版本以来,至今已有整整 6 年时间。从一开始只有 PyTorch 的 Tutorial 至今,很多网站和博客都介绍过 PyTorch 的安装和使用方法。

最近,国内的人民邮电出版社推出了一本关于 PyTorch 的教材《PyTorch 深度学习实战》,作者是:Eli Stevens,Luca Antiga,Thomas Viehmann,译者是:牟大恩。本书适用于学校的学生和想了解 PyTorch 的公司开发人员,通过学习这本书,大家都能够掌握使用 PyTorch 开发深度学习模型。

PyTorch 深度学习实战

本书的章节分成三个部分:

  • 第一部分:PyTorch 核心,介绍了深度学习的基础概念和一些简单的神经网络,张量等知识;
  • 第二部分:使用现实数据来学习,包括肺癌的早期检测,训练分类器,增强学习,语义分割等知识;
  • 第三部分:模型的部署,使用 Flask 服务来部署模型,导出模型等知识。

通常情况下,很多 Tutorial 和网络上的博客只会介绍前两部分,毕竟学会了前两部分就可以正常的开发模型,清洗数据。在学术界的绝大部分工作一般也只需要掌握前两部分,在工业界的一些 AI Labs 前两部分也能够满足绝大部分的需求。但是对于业务部门而言,只有前两部分的知识可能就不太够了,还需要有在线上部署模型,提供模型使用的能力。

在移动设备被大家大量使用的今天,很多情况下需要将模型部署在移动设备上,那么仅仅学会离线训练模型,进行数据的预测是远远不够的。有的情况下需要进行必要的模型压缩。除了 Python 之外,可能还需要掌握一些 C++ 的使用方法,同时掌握多个技能,才能够在工业界的业务部门完成各种各样的工作。

除了 Linux 和 macOS 之外,PyTorch 在 2018 年开始逐渐支持 Windows 系统。对于学生党而言,如果只有 Windows 系统,也能够基于现有的电脑进行深度学习的模型开发。

本书会教会读者如何使用 PyTorch 创建神经网络和深度学习系统,它帮助读者快速从零开始构建一个真实示例:肿瘤图像分类器。在构建模型过程中,它涵盖了整个深度学习管道的关键实践,包括 PyTorch 张量 API,用 Python 加载数据,监控训练以及将结果进行可视化展示。

本书主要内容会包括以下几个部分:

  • 训练深层神经网络;
  • 实现模块和损失函数;
  • 使用 PyTorch Hub 预先训练的模型;
  • 探索在 Jupyter Notebooks 中编写示例代码。

本书适用于对深度学习感兴趣的 Python 程序员,了解深度学习的基础知识对阅读本书有一定的帮助,但读者无须具有使用 PyTorch 或其他深度学习框架的经验。本书也有配套的 Python 源代码,方便读者进行下载和学习。

社交网络之间的帐号映射

帐号映射的整体介绍

在现实生活中,用户通常会同时使用多个社交网络,例如国外的 Twitter,Instagram,Facebook,也可能使用微信,QQ,微博等国内的产品。基于这些产品的不同定位,用户自身的社交网络会有很大的差异,那么如何通过机器学习算法找到一个人的社交网络帐号就成为了一个有趣的问题。在学术界,有学者对开源的 Facebook,Twitter 等社交网络数据进行了研究,设计了一套帐号映射的技术方案。

帐号映射这个课题有很多的别名,例如:

  • Social identity linkage;
  • User identity linkage;
  • User identity resolution;
  • Social network reconciliation;
  • User account linkage inference;
  • Profile linkage;
  • Anchor link prediction;
  • Detecting me edges;

帐号映射的目的就是将社交网络上这些看似不同的帐号映射到自然人:帐号(user accounts)-> 真实的自然人(real natural person)。令社交网络 \mathcal{G}=\mathcal{G}(\mathcal{U},\mathcal{E}) 是一个图,顶点是帐号 \mathcal{U}=\{u_{1},\cdots,u_{n}\}, 边是由帐号之间的连线 \mathcal{E}\subseteq \mathcal{U}\times\mathcal{U} 所构成的。

帐号映射(User Identity Linkage)的定义是:给定两个社交网络 \mathcal{G}^{s}=(\mathcal{U}^{s},\mathcal{E}^{s})\mathcal{G}^{t}=(\mathcal{U}^{t},\mathcal{E}^{t}), 其目标是找到一个函数 \mathcal{F}:\mathcal{U}^{s}\times\mathcal{U}^{t}\rightarrow\{0,1\} 使得,
\mathcal{F}(u^{s},u^{t})=\begin{cases}1,\text{ if } u^{s} \text{ and } u^{t} \text{ belong to same person,}\\ 0, \text{ otherwise.} \end{cases}
其中 u^{s}\in\mathcal{U}^{s}, u^{t}\in\mathcal{U}^{t}.

上述函数 \mathcal{F} 就是模型需要学习的目标函数,进一步地,对于 u^{s}\in\mathcal{U}^{s}, u^{t}\in\mathcal{U}^{t}, 学习得到的预测函数 \mathcal{\hat{F}}(u^{s},u^{t})=p\in [0,1] 表示两个帐号属于同一个自然人的概率值。

一个帐号在社交网络的属性包括很多方面,例如:

画像属性(Profile Features)内容属性(Content Features)社交网络属性(Network Features)
ID(社交网络的唯一标识)时间戳(timestamp)关注,被关注,好友关系(Friendship)
身份证 ID(identity card)语音(speech)点赞(like)
手机号(phone)视频(video)评论(comment)
昵称(username)图片(image)@(at)
头像(head image)文本(text)收藏(collect)
性别(gender)设备信息(device)消息(message)
年龄(age)wifi 信息(wifi)回复(reply)
邮箱(email)地理位置(gps)
个人网页(url)
职业(occupation)
常见的社交网络数据

一般情况下,

  1. 画像属性:绝大多数帐号都会有基础的画像信息;
  2. 内容属性:不活跃的帐号较难获取;
  3. 社交网络属性:线上的社交网络关系并不代表线下的社交网络关系,存在一定的噪声数据。

帐号映射的技术框架可以基于特征工程来做,然后使用有监督算法,无监督算法,或者半监督算法来进行帐号对之间的训练和预测。

帐号映射的技术框架

帐号映射的特征工程

帐号的画像特征

对于社交网络 \mathcal{G}=\mathcal{G}(\mathcal{U},\mathcal{E}) 中的一个帐号 u\in\mathcal{U} 而言,用 \overrightarrow{p_{u}}=(p_{u}^{1},\cdots,p_{u}^{m}) 来表示画像特征向量,其中 m\geq 1 表示画像属性特征的个数。对于两个社交网络 \mathcal{G}^{s},\mathcal{G}^{t} 的帐号 u^{s},u^{t} 而言,可以得到相应的画像特征向量 \overrightarrow{p_{u^{s}}},\overrightarrow{p_{u^{t}}}, 然后可以用基于距离(distance-based)或者基于频率(frequence-based)的方法来获得向量的距离或者相似性。换句话说,就是通过加权平均算法来获得结果:

sim(\overrightarrow{p_{u^{s}}},\overrightarrow{p_{u^{t}}})=\sum_{1\leq i\leq m}w_{i}\cdot sim(p_{u^{s}}^{i},p_{u^{t}}^{i}), 或者

dis(\overrightarrow{p_{u^{s}}},\overrightarrow{p_{u^{t}}})=\sum_{1\leq i\leq m}w_{i}\cdot dis(p_{u^{s}}^{i},p_{u^{t}}^{i}),

其中 sim 表示相似度,dis 表示距离。

基于距离的方法(distance-based)的方法很多,例如:

  1. 文本(text field)之间的距离可以考虑用 Jaro-Winkler distance,Jaccard similarity,Levenshtein distance 等方法;
  2. 图像(visual field)之间的距离可以考虑用 mean square error,dot product,angular distance,peak signal-to-noise ratio,Levenshtein distance 等方法;

基于频率的方法(frequence-based)可以考虑 bag-of-word model,TF-IDF model,Markov-chain model 等方法;

帐号的内容特征

帐号所产生的内容数据包括三个部分:

  1. 时间上的数据(temporal):时间戳的数据;
  2. 空间上的数据(spatial):帐号的设备数据,IP 数据,WIFI 数据,地理位置数据等等;
  3. 内容上的数据(post):帐号所产生的内容数据,包括但不限于视频,文本,语音,图片等。

用数学公式来描述就是 \overrightarrow{c_{u}}=\{(t_{1},s_{1},p_{1}), \cdots,(t_{m},s_{m},p_{m})\}, 其中 t_{i},s_{i},p_{i}(1\leq i\leq m) 分别表示时间戳,空间数据,内容数据。

在某个时间段内,帐号所产生的内容特征可以提炼出用户在社交网络上的行为数据和内容数据,形成一个行为序列。通过这个行为序列,可以得到用户的内容特征。

  1. 基于兴趣的特征(Interest-based):可以基于内容数据判断帐号对哪些内容更感兴趣;
  2. 基于风格的特征(Style-based):基于内容数据得到帐号的写作风格,例如常用词语等;
  3. 基于轨迹的特征(Trajectory-based):基于帐号的行为轨迹数据,包括设备,IP,WIFI,地理位置以及相应的时间戳,得到帐号的足迹(footprint)。

帐号的社交网络特征

社交网络包括两种:

  1. 局部社交网络(local network):查看帐号的邻居(关注,被关注,好友关系)等诸多数据;
  2. 全局社交网络(global network):查看帐号在全局数据中的位置情况;

对于帐号的社交网络特征,包括以下两种常见形式:

  1. 基于邻居的特征(Neighborhood-based):共同好友数,共同邻居个数,Jaccard Coefficient,Overlap Coefficient,Dice Coefficient,Adamic/Adar score;
  2. 基于嵌入的特征(Embedding-based):通过计算帐号在相应的社交网络的嵌入特征,然后计算特征之间的距离或者相似性。

帐号映射的建模思路

机器学习的常见算法包括有监督算法(supervised model),无监督算法(unsupervised model)和半监督算法(semi-supervised model)。

基于上面的特征工程,加上合适的权重之和可以得到一个分数(score),也就是:

\mathcal{F}(u^{s},u^{t})=\alpha \cdot sim_{p}(\overrightarrow{p_{u^{s}}}, \overrightarrow{p_{u^{t}}})+\beta\cdot sim_{c}(\overrightarrow{p_{u^{s}}}, \overrightarrow{p_{u^{t}}})+\gamma\cdot sim_{n}(\overrightarrow{p_{u^{s}}}, \overrightarrow{p_{u^{t}}}),

其中 sim_{p}, sim_{c}, sim_{n} 分别表示画像(profile),内容(content),社交网络(network)之间的相似度。

如果有样本的话,全体样本是 \mathcal{Q}=\{(u^{s},u^{t}),\forall u^{s}\in\mathcal{U}^{s},\forall u^{t}\in\mathcal{U}^{t}\}, 那么正样本是 \mathcal{M}=\{(u^{s},u^{t}),u^{s}\in\mathcal{U}^{s},u^{t}\in\mathcal{U}^{t}, u^{s} \text{ and } u^{t} \text{ belong to the same person}\}, 负样本 \mathcal{N}=\mathcal{Q}-\mathcal{M}. 在实际使用的时候,要注意采样的比例和负样本的选择方法。

帐号映射的评价指标

对于 u^{s}\in\mathcal{U}^{s}, u^{t}\in\mathcal{U}^{t}, 假设
\mathcal{Q}=\{(u^{s},u^{t}),\forall u^{s}\in\mathcal{U}^{s},\forall u^{t}\in\mathcal{U}^{t}\} 是所有的帐号对,
\mathcal{M}=\{(u^{s},u^{t}),u^{s}\in\mathcal{U}^{s},u^{t}\in\mathcal{U}^{t}, u^{s} \text{ and } u^{t} \text{ belong to the same person}\} 是所有属于相同自然人的帐号对,
\mathcal{N}=\mathcal{Q}-\mathcal{M} 是所有属于不同自然人的帐号对。
\mathcal{A}= {被算法映射成相同自然人的帐号对},\mathcal{B}= {被算法映射成不同自然人的帐号对};用 TP, TN, FN, FP 来描述就是:

  • True Positive(TP):|\mathcal{A} \cap \mathcal{M}|;
  • True Negative(TN):|\mathcal{B}\cap\mathcal{N}|;
  • False Negative(FN):|\mathcal{B}\cap\mathcal{M}|;
  • False Positive(FP):|\mathcal{A}\cap\mathcal{N}|;

那么,Precision = TP / (TP + FP)=|\mathcal{A}\cap\mathcal{M}|/|\mathcal{A}|, Recall = TP / (TP + FN)= |\mathcal{A}\cap\mathcal{M}|/|\mathcal{M}|.

另外,F1=2\cdot Precision \cdot Recall / (Precision + Recall), Accuracy=(TP+TN)/(TP+TN+FP+FN).

部分论文细节

  1. Liu, Jing, et al. “What’s in a name? An unsupervised approach to link users across communities.” Proceedings of the sixth ACM international conference on Web search and data mining. 2013. 本篇文章主要是基于用户的名字来识别跨网络的用户的,提取用户的特征之后,使用 SVM 分类器来进行识别;
  2. Riederer, Christopher, et al. “Linking users across domains with location data: Theory and validation.” Proceedings of the 25th International Conference on World Wide Web. 2016. 本篇文章主要是基于用户的内容特征来进行建模;
  3. Labitzke, Sebastian, Irina Taranu, and Hannes Hartenstein. “What your friends tell others about you: Low cost linkability of social network profiles.” Proc. 5th International ACM Workshop on Social Network Mining and Analysis, San Diego, CA, USA. 2011. 本篇论文是根据社交网络中的用户邻居数据,来判断用户之间相似性的。

参考文献

  1. Shu, Kai, et al. “User identity linkage across online social networks: A review.” Acm Sigkdd Explorations Newsletter 18.2 (2017): 5-17.
  2. Liu, Jing, et al. “What’s in a name? An unsupervised approach to link users across communities.” Proceedings of the sixth ACM international conference on Web search and data mining.
  3. Riederer, Christopher, et al. “Linking users across domains with location data: Theory and validation.” Proceedings of the 25th International Conference on World Wide Web. 2016.
  4. Labitzke, Sebastian, Irina Taranu, and Hannes Hartenstein. “What your friends tell others about you: Low cost linkability of social network profiles.” Proc. 5th International ACM Workshop on Social Network Mining and Analysis, San Diego, CA, USA. 2011.

近似最近邻搜索算法 ANNOY(Approximate Nearest Neighbors Oh Yeah)

在搜索的业务场景下,基于一个现有的数据候选集(dataset),需要对新来的一个或者多个数据进行查询(query),返回在数据候选集中与该查询最相似的 Top K 数据。

Google:the two-tower neural network model

最朴素的想法就是,每次来了一个新的查询数据(query),都遍历一遍数据候选集(dataset)里面的所有数据,计算出 query 与 dataset 中所有元素的相似度或者距离,然后精准地返回 Top K 相似的数据即可。

但是当数据候选集特别大的时候,遍历一遍数据候选集里面的所有元素就会耗费过多的时间,其时间复杂度是 O(n), 因此,计算机科学家们开发了各种各样的近似最近邻搜索方法(approximate nearest neighbors)来加快其搜索速度,在精确率和召回率上面就会做出一定的牺牲,但是其搜索速度相对暴力搜索有很大地提高。

在这个场景下,通常都是欧式空间里面的数据,形如 \bold{x}=(x_{1},\cdots,x_{n})\in \mathbb{R}^{n}, 其中 n 是欧氏空间的维度。常用的距离公式包括:

  • Manhattan 距离:L1 范数;
  • Euclidean 距离:L2 范数;
  • Cosine 距离:1 – Cosine 相似度;
  • 角距离:用两个向量之间的夹角来衡量两者之间的距离;
  • Hamming 距离:一种针对 64 维的二进制数的 Manhattan 距离,相当于 \mathbb{R}^{64} 中的 L1 范数;
  • Dot Product 距离\bold{x}\cdot \bold{y}=\sum_{i=1}^{n}x_{i}\cdot y_{i}.
欧氏空间的数据点聚类

在近似最近邻搜索(ANN)领域,有很多开源的算法可以使用,包括但不限于:

  • Annoy(Approximate Nearest Neighbors Oh Yeah);
  • ScaNN(Scalable Nearest Neighbors);
  • Faiss(Billion-scale similarity search with GPUs);
  • Hnswlib(fast approximate nearest neighbor search);

本文将会重点介绍 Annoy 算法及其使用案例;

ANN 的 benchmark

Annoy 的算法思想

本文以 \mathbb{R}^{2} 中的点集来作为案例,介绍 annoy 算法的基本思想和算法原理。

二维欧氏空间中的点集

用 n 表示现有的文档个数,如果采用暴力搜索的方式,那么每次查询的耗时是 O(n), 采用合适的数据结构可以有效地减少查询的耗时,在 annoy 算法中,作者采用了二叉树这个数据结构来提升查询的效率,目标是把查询的耗时减少至 O(\ln(n)).

用中垂线来作为切分平面

刚开始的时候,在数据集中随机选择两个点,然后用它们的中垂线来切分整个数据集,于是数据集就被分成了蓝绿两个部分。然后再随机两个平面中各选出一个顶点,再用中垂线进行切分,于是,整个平面就被切成了四份。

继续切分

用一颗二叉树来表示这个被切分的平面就是:

用二叉树来表示切分平面

后续继续采用同样的方式进行切分,直到每一个平面区域最多拥有 K 个点为止。当 K = 10 时,其相应的切分平面和二叉树如下图所示。

持续随机切分
K = 10 的时候,每个子平面最多 10 个点
相应的二叉树

下面,新来的一个点(用红色的叉表示),通过对二叉树的查找,我们可以找到所在的子平面,然后里面最多有 K = 10 个点。从二叉树的叶子节点来看,该区域只有 7 个点。

红叉表示需要查询的点
搜索路径

在 ANN 领域,最常见的两个问题是:

  1. 如果我们想要 Top K 的点,但是该区域的点集数量不足 K,该怎么办?
  2. 如果真实的 Top K 中部分点不在这个区域,该怎么办?

作者用了两个技巧来解决这个问题:

  1. 使用优先队列(priority queue):将多棵树放入优先队列,逐一处理;并且通过阈值设定的方式,如果查询的点与二叉树中某个节点比较相似,那么就同时走两个分支,而不是只走一个分支;
  2. 使用森林(forest of trees):构建多棵树,采用多个树同时搜索的方式,得到候选集 Top M(M > K),然后对这 M 个候选集计算其相似度或者距离,最终进行排序就可以得到近似 Top K 的结果。

同时走两个分支的的示意图:

随机生成多棵树,构建森林的示意图:

随机生成多颗树

Top K 的查询方法:

选择候选区域
获得区域中的所有点
区域中 Top K 点排序

Annoy 算法原理:

构建索引:建立多颗二叉树,每颗二叉树都是随机切分的;

查询方法
1. 将每一颗树的根节点插入优先队列;
2. 搜索优先队列中的每一颗二叉树,每一颗二叉树都可以得到最多 Top K 的候选集;
3. 删除重复的候选集;
4. 计算候选集与查询点的相似度或者距离;
5. 返回 Top K 的集合。

Annoy 的编程实践

Annoy 的安装:

pip install annoy

Annoy 的 Python 接口函数

常用的 Annoy Python 接口函数包括以下内容:

  • a = AnnoyIndex(f, metric):f 指的是向量的维度,metric 表示度量公式。在这里,Annoy 支持的度量公式包括:”angular”, “euclidean”, “manhattan”, “hamming”, “dot”;
  • a.add_item(i, v):i 是一个非负数,表示 v 是第 i 个向量;
  • a.build(n_trees, n_jobs=-1):n_trees 表示树的棵数,n_jobs 表示线程个数,n_jobs=-1 表示使用所有的 CPU 核;
  • a.save(fn, prefault=False):表示将索引存储成文件,文件名是 fn;
  • a.load(fn, prefault=False):表示将索引从文件 fn 中读取出来;
  • a.unload():表示不再加载索引;
  • a.get_nns_by_item(i, n, search_k=-1, include_distances=False):返回在索引中的第 i 个向量 Top n 最相似的向量;如果不提供 search_k 值的话,search_k 默认为 n_tree * n,该指标用来平衡精确度和速度;includ_distances=True 表示返回的时候是一个包含两个元素的 tuple,第一个是索引向量的 index,第二个就是相应的距离;
  • a.get_nns_by_vector(v, n, search_k=-1, include_distances=False):返回与向量 v Top n 最相似的向量;如果不提供 search_k 值的话,search_k 默认为 n_tree * n,该指标用来平衡精确度和速度;includ_distances=True 表示返回的时候是一个包含两个元素的 tuple,第一个是索引向量的 index,第二个就是相应的距离;
  • a.get_item_vector(i):返回添加索引的时候的第 i 个向量;
  • a.get_distance(i, j):返回第 i 个向量与第 j 个向量的距离;
  • a.get_n_items():返回索引中的向量个数;
  • a.get_n_trees():返回索引中的树的棵数;
  • a.on_disk_build(fn):在一个文件 fn 中构建索引,并不在内存中构建索引;
  • a.set_seed(seed):用给定的种子初始化随机数生成器,需要在建树,添加向量构建索引之前使用该函数;

影响 annoy 算法效率和精度的重要参数:

  • n_trees:表示树的棵数,会影响构建索引的时间。值越大表示最终的精度越高,但是会有更多的索引;
  • search_k:值越大表示搜索耗时越长,搜索的精度越高;如果需要返回 Top n 最相似的向量,则 search_k 的默认值是 n_trees * n;

Annoy 的使用案例

from annoy import AnnoyIndex
import random

# f 表示向量的维度
f = 40
# 'angular' 是 annoy 支持的一种度量;
t = AnnoyIndex(f, 'angular')  # Length of item vector that will be indexed
# 插入数据
for i in range(1000):
    v = [random.gauss(0, 1) for z in range(f)]
    # i 是一个非负数,v 表示向量
    t.add_item(i, v)

# 树的数量
t.build(10) # 10 trees

# 存储索引成为文件
t.save('test.ann')

# 读取存储好的索引文件
u = AnnoyIndex(f, 'angular')
u.load('test.ann') # super fast, will just mmap the file

# 返回与第 0 个向量最相似的 Top 100 向量;
print(u.get_nns_by_item(0, 1000)) # will find the 1000 nearest neighbors

# 返回与该向量最相似的 Top 100 向量;
print(u.get_nns_by_vector([random.gauss(0, 1) for z in range(f)], 1000))

# 返回第 i 个向量与第 j 个向量的距离;
# 第 0 个向量与第 0 个向量的距离
print(u.get_distance(0, 0))
# 第 0 个向量与第 1 个向量的距离
print(u.get_distance(0, 1))

# 返回索引中的向量个数;
print(u.get_n_items())
# 返回索引中的树的棵数;
print(u.get_n_trees())

# 不再加载索引
print(u.unload())

基于 hamming 距离的 annoy 使用案例:

from annoy import AnnoyIndex

# Mentioned on the annoy-user list
bitstrings = [
 '0000000000011000001110000011111000101110111110000100000100000000',
    '0000000000011000001110000011111000101110111110000100000100000001',
    '0000000000011000001110000011111000101110111110000100000100000010',
    '0010010100011001001000010001100101011110000000110000011110001100',
    '1001011010000110100101101001111010001110100001101000111000001110',
    '0111100101111001011110010010001100010111000111100001101100011111',
    '0011000010011101000011010010111000101110100101111000011101001011',
    '0011000010011100000011010010111000101110100101111000011101001011',
    '1001100000111010001010000010110000111100100101001001010000000111',
    '0000000000111101010100010001000101101001000000011000001101000000',
    '1000101001010001011100010111001100110011001100110011001111001100',
    '1110011001001111100110010001100100001011000011010010111100100111',
]

# 将其转换成二维数组
vectors = [[int(bit) for bit in bitstring] for bitstring in bitstrings]

# 64 维度
f = 64
idx = AnnoyIndex(f, 'hamming')
for i, v in enumerate(vectors):
    idx.add_item(i, v)

# 构建索引
idx.build(10)
idx.save('idx.ann')
idx = AnnoyIndex(f, 'hamming')
idx.load('idx.ann')
js, ds = idx.get_nns_by_item(0, 5, include_distances=True)

# 输出索引和 hamming 距离
print(js, ds)

基于欧几里得距离的 annoy 使用案例:

from annoy import AnnoyIndex
import random

# f 表示向量的维度
f = 2
# 'euclidean' 是 annoy 支持的一种度量;
t = AnnoyIndex(f, 'euclidean')  # Length of item vector that will be indexed
# 插入数据
t.add_item(0, [0, 0])
t.add_item(1, [1, 0])
t.add_item(2, [1, 1])
t.add_item(3, [0, 1])

# 树的数量
t.build(n_trees=10) # 10 trees
t.save('test2.ann')

u = AnnoyIndex(f, 'euclidean')
u.load('test2.ann') # super fast, will just mmap the file

print(u.get_nns_by_item(1, 3)) # will find the 1000 nearest neighbors

print(u.get_nns_by_vector([0.1, 0], 3))

参考资料

  1. GitHub 的 Annoy 开源代码:https://github.com/spotify/annoy
  2. Nearest neighbors and vector models – part 2 – algorithms and data structures:https://erikbern.com/2015/10/01/nearest-neighbors-and-vector-models-part-2-how-to-search-in-high-dimensional-spaces.html
  3. ann-benchmark 的效果测试:https://github.com/erikbern/ann-benchmarks

字符串相似度的数学原理和开源工具

在 DNA 测序,蛋白质测序,计算语言学等研究领域,其研究对象可以是一个字符串,也可以是一个短文本,甚至一篇完整的文章。例如:

  1. 在蛋白质测序领域,其案例形如:Cys-Gly-Leu-Ser-Trp;
  2. 在 DNA 测序领域,其案例形如:AGCTTCGA;
  3. 在计算语言学领域,其案例形如:it is rainy today;

在以上的科学领域中,如何对字符串进行研究就显得尤其重要,一种方式是直接对字符和字符串本身进行研究与分析,另一种方式是对字符,单词,句子进行嵌入式的操作,将其映射成欧式空间中的向量,然后再对其进行数据分析和机器学习建模。本篇文章将会从字符与字符串本身出发,分析字符串与字符串之间的性质。

如果需要研究字符串之间的相似性与距离,大致有以下两种方案:

  1. 基于序列的方法(sequence-based);
  2. 基于集合的方法(set-based);

顾名思义,可以将字符串看成序列,然后使用序列的相似性或者距离来计算两者之间的性质;也可以将字符串看成集合(或者多重集合),然后使用集合(或者多重集合)的相似性来计算两者之间的性质。

基于集合的方法

n 元语法

n 元语法(n-gram)是自然语言处理中最常见的模型之一,它指的是文本中连续出现的 n 个词语。例如 X="\text{it is rainy today.}" 其 n 元语法分别是:

  1. 1 元语法(unigram):it, is, rainy, today
  2. 2 元语法(bigram):it is, is rainy, rainy today
  3. 3 元语法(trigram):it is rainy, is rainy today

对于一个长度为 m 的字符串 X=X[0,\cdots,m-1]=[X[0],\cdots,X[m-1]] 而言,其 unigram,bigram,trigram,multiset(多重集合)分别是:

  1. 1 元语法集合:unigram(X)=\{X[i], 0\leq i\leq m-1\};
  2. 2 元语法集合:bigram(X)=\{X[i]X[i+1], 0\leq i\leq m-2\};
  3. 3 元语法集合:trigram(X)=\{X[i]X[i+1]X[i+2], 0\leq i\leq m-3\};
  4. 1 元语法多重集合:multiset(X)=\{X[i]:d[i], 0\leq i\leq m-1\},

其中 d[i] 表示 X[i] 出现的次数。通过 n 元语法,我们可以将一个字符串转换成一个集合,然后通过计算集合之间的相似性来评估字符串的相似性。对于字符串 X,Y 而言,其相似度可以转换为:

sim(X,Y)=\begin{cases}fun\_sim(unigram(X),unigram(Y)),\\ fun\_sim(bigram(X),bigram(Y)),\\ fun\_sim(trigram(X),trigram(Y)).\end{cases}

其中 fun\_sim 表示集合的相似度计算函数。

dis(X,Y)=\begin{cases}fun\_dis(unigram(X),unigram(Y)),\\ fun\_dis(bigram(X),bigram(Y)),\\ fun\_dis(trigram(X),trigram(Y)).\end{cases}

其中 fun\_dis 表示集合的距离计算函数。

集合的相似性或者距离有很多方法可以计算。对于集合 A,B 而言,其 fum\_sim 的选型就包括但不限于以下几种:

\text{Jaccard Coefficient}(A,B)=|A\cap B|/|A\cup B|;

\text{Overlap Coefficient}(A,B)=|A\cap B|/\min(|A|,|B|);

\text{Dice Coefficient}(A,B)=2\cdot|A\cap B|/(|A|+|B|);

\text{Cosine Coefficent}(A,B)=|A\cap B|/\sqrt{|A|\cdot|B|};

\text{Tversky Index}(A,B)=|A\cap B|/(|A\cap B|+\alpha\cdot|A-B|+\beta\cdot|B-A|);

其中,\alpha, \beta \geq 0. 同时,不仅集合之间可以进行交集,并集的计算,多重集合之间同样可以进行类似的操作,于是上述方法同样可以应用在多重集合上。

而在字符串距离中,有一个特殊的方法叫做 bag distance,也就是将字符串转换成 unigram 的多重集合,然后计算其距离。令 bag(X),bag(Y) 表示字符串 X,Y 的 unigram 多重集合,然后 bag distance 就可以定义为 \max(|bag(X)-bag(Y)|,|bag(Y)-bag(X)|). 其中 |\cdot| 表示多重集合的 SIZE,减法 - 表示多重集合的差集。

基于序列的方法

最长公共子序列与最长公共子串

在数据结构这门课中,在讲解动态规划的部分,一般都会提到最长公共子序列(Longest Common Subsequence)和最长公共子串(Longest Common Substring)。而子序列和子串其实定义是不一样的。对于序列 [x_{0},x_{1},x_{2},\cdots,x_{n}] 而言,其子序列(subsequence) [x_{n_{1}},\cdots,x_{n_{k}}] 指的是从原始的序列中通过去除某些元素但不破坏余下元素的相对位置(在前或者在后)而形成的新序列。子串(substring)是相对于一个字符串而言,它是其原始字符串中的完整一段。例如:对于“苹果手机”而言,“苹手”是其子序列,但“苹手”并不是子串。

在这里,如果需要计算两个字符串 XY 的最长公共子序列的长度和最长公共子串的长度,就需要使用动态规划方面的知识,构建其边界条件和动态转移方程。

首先,我们来计算两个字符串 XY 的最长公共子序列。假设 mn 分别是字符串 XY 的长度,i.e. X=X[0,\cdots,m-1], Y=Y[0,\cdots,n-1].

L(m,n)=L(X[0,\cdots,m-1], Y[0,\cdots,n-1]) 表示字符串 X[0,\cdots,m-1]Y[0,\cdots,n-1] 的最长公共子序列的长度。则可以得到其状态转移方程如下:

L(m,n)=\begin{cases} L(m-1,n-1)+1, \text{ if } X[m-1] == Y[n-1] \\ \max(L(m-1,n), L(m, n-1)), \text{ else }. \end{cases}

其边界条件是 L(0,j)=L(i,0)=0, \forall 0\leq i\leq m, 0\leq j\leq n. 返回 L(m,n) 即可表示最长公共子序列的长度。

例如,X="abcd", Y="ced", 则有

L(4,3)=L("abcd", "ced")=L("abc","ce")+1

=\max(L("abc", "c"),L("ab","ce"))+1=\max(1,0)+1=2.

其次,我们来计算两个字符串 XY 的最长公共子串。假设 mn 分别是字符串 XY 的长度,i.e. X=X[0,\cdots,m-1], Y=Y[0,\cdots,n-1].

L(m,n)=L(X[0,\cdots,m-1], Y[0,\cdots,n-1]) 表示字符串X[0,\cdots,m-1]Y[0,\cdots,n-1] 的最长公共子串的长度。则可以得到其状态转移方程如下:

L(m,n)=\begin{cases} L(m-1,n-1)+1, \text{ if } X[m-1]==Y[n-1] \\ 0, \text{ else }. \end{cases}

其边界条件是 L(0,j)=L(i,0)=0, \forall 0\leq i\leq m, 0\leq j\leq n. 返回:\max_{0\leq i\leq m,0\leq j\leq n}L(i,j) 即可表示最长公共子串的长度。

例如:X = "abc", Y="bcd", 则有最长公共子串的长度是 L(3,3)=2.

Jaro 和 Jaro-Winkler 相似度

首先,我们来定义两个字符串之间的 Jaro 相似度。对于两个字符串 X, Y 而言,如果两个字符 X[i], Y[j] 满足以下两个条件:

  1. X[i]=Y[j];
  2. |i-j|\leq [\max(|X|,|Y|)/2]-1;

X[i],Y[j] 被称为匹配(matching),其中 |\cdot| 表示 string 的长度,[\cdot] 表示高斯取整函数。在此定义下计算出 X,Y 的匹配字符为 X',Y'. 从定义可以得到 X',Y' 的长度是一样的,记为 m.t=[\#\{0\leq i\leq m-1:X'[i]\neq Y'[i]\}/2] 称为 transposition。于是,Jaro 相似度就可以定义为:

\text{Jaro Similarity}(X,Y)=\begin{cases}0,\text{ if } m=0;\\ \frac{1}{3}\bigg(\frac{m}{|X|}+\frac{m}{|Y|}+\frac{m-t}{m}\bigg).\end{cases}

例如,

  1. X="martha",Y="marhta", 可以得到 X'="martha",Y'="marhta", 于是 m=6,t=1, 因此,Jaro 相似度是
    Jaro("martha","marhta")=\frac{1}{3}\bigg(\frac{6}{6}+\frac{6}{6}+\frac{6-1}{6}\bigg)=0.94.
  2. X="DWAYNE",Y="DUANE", 可以得到 X'="DANE",Y'="DANE", 于是 m=4,t=0, 因此,Jaro 相似度是
    Jaro("DWAYNE","DUANE")=\frac{1}{3}\bigg(\frac{4}{6}+\frac{4}{5}+\frac{4-0}{4}\bigg)=0.82.
  3. X="DIXON",Y="DICKSONX", 可以得到 X'="DION",Y'="DION", 于是 m=4,t=0, 因此,Jaro 相似度是
    Jaro("DIXON","DICKSONX")=\frac{1}{3}\bigg(\frac{4}{5}+\frac{4}{8}+\frac{4-0}{4}\bigg)=0.77.
  4. X="arnab",Y="aranb", 可以得到 X'="arnab",Y'="aranb", 于是 m=5,t=1, 因此,Jaro 相似度是
    Jaro("arnab","aranb")=\frac{1}{3}\bigg(\frac{5}{5}+\frac{5}{5}+\frac{5-1}{5}\bigg)=0.93.

其次,我们来定义两个字符串之间的 Jaro-Winkler 相似度。对于字符串 X,Y 而言,其 Jaro-Winlker 相似度定义为:

\text{Jaro-Winkler Similarity}(X,Y)=\text{Jaro Similarity}(X,Y)

+ \ell\cdot p\cdot(1-\text{Jaro-Winkler Similarity}(X,Y)),

其中 p 表示系数(默认是 0.1,可以调整),\ell 表示 |X|, |Y| 的最长前缀子串的 SIZE,并且不超过 4。

例如:X="arnab",Y="aranb", 可以得到 \ell=2. 因此,Jaro-Winkler 相似度为 0.9333333+0.1\cdot 2\cdot(1-0.9333333)=0.946667.

Levenshtein 距离

Levenshtein 距离是编辑距离的一种形式,所谓编辑距离,指的是在两个字符串之间,由一个转成另外一个所需要的最少编辑次数。在这里的编辑操作包括:

  1. 更换:将一个字符更换为另一个字符;
  2. 删除:删除一个字符;
  3. 插入:插入一个字符;

对于字符串 X,Y 而言,|X|,|Y| 分别表示字符串 X,Y 的长度。则可以用动态规划的思想来计算 Levenshtein 距离。用 Lev(i,j) 表示 X 的前 i 个字符和 Y 的前 j 个字符的 Levenshtein 距离,L(|X|,|Y|) 表示两个字符串的 Levenshtein 距离。

Lev(i,j)=\begin{cases} j, \text{ if } i=0,\\ i, \text{ else if } j=0,\\ \min\begin{cases}Lev(i-1,j)+1,\\ Lev(i,j-1)+1,\\ Lev(i-1,j-1)+1_{X[i-1]\neq Y[j-1]} \end{cases}\text{ otherwise }.\end{cases}

其中 1 表示指示函数,1\leq i\leq |X|, 1\leq j\leq |Y|.

例如:kitten -> sitting 的 Levenshtein 距离就是 3。因为

  1. kitten -> sitten:替换 k -> s;
  2. sitten -> sittin:替换 e -> i;
  3. sittin -> sitting:增加 g.

距离与相似度的转换

对于相似度 sim 和距离 dist 而言,只要找到某个递减函数就可以将其互相转换。

  1. f:[0,+\infty)\rightarrow [-1,1] 是将距离转换成相似度的函数,则可以表示为 sim = f\circ dist, 其中,f 是严格递减函数,值域属于 [-1,1], f(0)=1.
  2. g:[-1,1]\rightarrow [0,+\infty) 是将相似度转换为距离的函数,则可以表示为 dist=g\circ sim, 其中,g 是严格递减函数,值域属于 [0,+\infty), g(1)=0.

例如:

  1. Levenshtein 相似度就可以用函数 1-lev(X,Y)/\max(|X|,|Y|) 来转换;
  2. 基于最长公共子串/最长公共子序列的相似度既可以使用函数 1-lcs(X,Y)/\max(|X|,|Y|), 也可以使用函数 1-lcs(X,Y)/\min(|X|,|Y|) 来计算。

开源工具 py_stringmatching

如果需要计算两个字符串的相似性,可以直接使用开源工具 py_stringmatching。通过 py_stringmatching 的官网可以看到,其主要功能分成两个部分,第一部分是 Tokenizers 的介绍,第二部分是各种相似度的计算。

Tokenizer

在第一部分 Tokenizers 中,本质上就是将一个字符串切分成一个序列。其主要的切分函数分成五种,分别是:

  1. Alphabetic Tokenizer:返回最长连续的英文序列;
  2. Alphanumeric Tokenizer:返回最长连续的英文/数字序列;
  3. Delimiter Tokenizer:根据某个指定的字符串来进行切分;
  4. Qgram Tokenizer:基于 Q 元语法的切分;
  5. Whitespace Tokenizer:基于空格的切分。

从以上的切分方式来看,py_stringmatching 更适用于英文,因为中文需要使用专门的切词工具。这里的 return_set=True 指的是返回 set(去除重复的元素)。

首先来看 Alphabetic Tokenize 的案例:

from py_stringmatching import AlphabeticTokenizer

al_tok = AlphabeticTokenizer()
print(al_tok.tokenize('algebra88analysis, geometry#geometry.'))

# 输出:['algebra', 'analysis', 'geometry', 'geometry']

al_tok = AlphabeticTokenizer(return_set=True)
print(al_tok.tokenize('algebra88analysis, geometry#geometry.'))

# 输出:['algebra', 'analysis', 'geometry']

其次来看 Alphanumeric 的案例:

from py_stringmatching import AlphanumericTokenizer
alnum_tok = AlphanumericTokenizer()
print(alnum_tok.tokenize('algebra9,(analysis), geometry#.(geometry).88'))
# 输出:['algebra9', 'analysis', 'geometry', 'geometry', '88']

alnum_tok = AlphanumericTokenizer(return_set=True)
print(alnum_tok.tokenize('algebra9,(analysis), geometry#.(geometry).88'))
# 输出:['algebra9', 'analysis', 'geometry', '88']

然后看 Delimiter Tokenizer 的案例:

from py_stringmatching import DelimiterTokenizer
delim_tok = DelimiterTokenizer()
print(delim_tok.tokenize('algebra analysis geometry  geometry'))
# 输出:['algebra', 'analysis', 'geometry', 'geometry']

delim_tok = DelimiterTokenizer(delim_set={'$ #$'})
print(delim_tok.tokenize('algebra$ #$analysis'))
# 输出:['algebra', 'analysis']

delim_tok = DelimiterTokenizer(delim_set={',', '.'})
print(delim_tok.tokenize('algebra,analysis,geometry.geometry'))
# 输出:['algebra', 'analysis', 'geometry', 'geometry']

delim_tok = DelimiterTokenizer(delim_set={',', '.'}, return_set=True)
print(delim_tok.tokenize('algebra,analysis,geometry.geometry'))
# 输出:['algebra', 'analysis', 'geometry']

再次来看 Q 元语法的案例:QgramTokenize 的参数包括:

  • qval = 2:q 元数组;
  • padding = True:是否需要加上前后缀;
  • prefix_pad =‘#’:前缀;
  • suffix_pad = ‘$’:后缀;
  • return_set = True:是否去重
from py_stringmatching import QgramTokenizer
qgram_tok = QgramTokenizer()
print(qgram_tok.tokenize('algebra'))
# 输出:['#a', 'al', 'lg', 'ge', 'eb', 'br', 'ra', 'a$']

qgram_tok = QgramTokenizer(qval=3)
print(qgram_tok.tokenize('algebra'))
# 输出:['##a', '#al', 'alg', 'lge', 'geb', 'ebr', 'bra', 'ra$', 'a$$']

qgram_tok = QgramTokenizer(padding=False)
print(qgram_tok.tokenize('algebra'))
# 输出:['al', 'lg', 'ge', 'eb', 'br', 'ra']

qgram_tok = QgramTokenizer(prefix_pad='^', suffix_pad='!')
print(qgram_tok.tokenize('algebra'))
# 输出:['^a', 'al', 'lg', 'ge', 'eb', 'br', 'ra', 'a!']

最后来看 Whitespace Tokenize 的案例:

from py_stringmatching import WhitespaceTokenizer
ws_tok = WhitespaceTokenizer()
print(ws_tok.tokenize('algebra analysis geometry geometry topology'))
# 输出:['algebra', 'analysis', 'geometry', 'geometry', 'topology']

print(ws_tok.tokenize('algebra analysis geometry  geometry topology'))
# 输出:['algebra', 'analysis', 'geometry', 'geometry', 'topology']

ws_tok = WhitespaceTokenizer(return_set=True)
print(ws_tok.tokenize('algebra analysis geometry  geometry topology'))
# 输出:['algebra', 'analysis', 'geometry', 'topology']

Similarity Measures

下面来介绍开源工具 py_stringmatching 的常见相似度函数。

Bag Distance,Levenshtein 函数的使用:

from py_stringmatching import BagDistance
bd = BagDistance()
print(bd.get_raw_score('algebra', 'algebraic'))
# 输出:2

print(bd.get_sim_score('algebra', 'algebraic'))
# 输出:0.7778

from py_stringmatching import Levenshtein
lev = Levenshtein()
print(lev.get_raw_score('algebra', 'algebraic'))
# 输出:2

print(lev.get_sim_score('algebra', 'algebraic'))
# 输出:0.7778

其中 bag,Levenshtein 相似度的定义是 1 - raw\_score/\max(|X|,|Y|), 其中 |\cdot| 表示字符串的长度, X,Y 表示字符串。

Cosine,Dice,Jaccard,OverlapCoefficient,TverskyIndex 函数的使用:

  • 这些函数的输入都是 set 或者 list;
  • 如果是需要比较字符串,则可以使用 Tokenizer 函数将其切分或者转换,例如 1 gram;
  • 这些函数的 get_raw_score 与 get_sim_score 是一样的,因为都是相似度函数。
from py_stringmatching import Cosine, Dice, Jaccard, OverlapCoefficient
cos = Cosine()
print(cos.get_raw_score(['algebra'], ['algebra', 'analysis']))
# 输出:0.7071

dice = Dice()
print(dice.get_raw_score(['algebra'], ['algebra', 'analysis']))
# 输出:0.6667

jaccard = Jaccard()
print(jaccard.get_raw_score(['algebra'], ['algebra', 'analysis']))
# 输出:0.5

overlap = OverlapCoefficient()
print(overlap.get_raw_score(['algebra'], ['algebra', 'analysis']))
# 输出:1.0

tversky = TverskyIndex()
print(tversky.get_raw_score(['algebra'], ['algebra', 'analysis']))
# 输出:0.6667

Jaro,Jaro Winkler 函数的使用:

  • 这两个函数的输入都是 string;
  • 这两个函数的 get_raw_score 与 get_sim_score 是一样的;
from py_stringmatching import Jaro
jaro = Jaro()
print(jaro.get_raw_score('algebra', 'algebraic'))
print(jaro.get_sim_score('algebra', 'algebraic'))
# 输出:0.9259

from py_stringmatching import JaroWinkler
jw = JaroWinkler()
print(jw.get_raw_score('algebra', 'algebraic'))
print(jw.get_sim_score('algebra', 'algebraic'))
# 输出:0.9556

Hamming 距离的使用:

  • 只能适用于长度一致的字符串;
  • 相似度的计算是通过 1-raw\_score/length 得到的;
from py_stringmatching import HammingDistance
hm = HammingDistance()
print(hm.get_raw_score('algebra', 'algebri'))
# 输出:1
print(hm.get_sim_score('algebra', 'algebri'))
# 输出:0.8571

开源工具 FuzzyWuzzy

FuzzyWuzzy 是一个计算 STRING 相似度的开源工具库,其值域是 [0,100]. 如果需要计算相似度,直接除以 100 即可。

from fuzzywuzzy import fuzz
print(fuzz.partial_ratio('algebra is interesting', 'algebraic is good'))
# 输出:65
print(fuzz.partial_token_sort_ratio('algebra is interesting', 'algebraic is good'))
# 输出:59
print(fuzz.partial_token_set_ratio('algebra is interesting', 'algebraic is good'))
# 输出:100
print(fuzz.token_sort_ratio('algebra is interesting', 'algebraic is good'))
# 输出:62
print(fuzz.token_set_ratio('algebra is interesting', 'algebraic is good'))
# 输出:62

开源工具 python-Levenshtein

Levenshtein 距离也可以使用开源工具 python-Levenshtein 来计算,形如:

from Levenshtein import distance
print(distance('algebra', 'algebraic'))
# 输出 2

参考资料:

  1. 最长公共子串:https://www.geeksforgeeks.org/longest-common-substring-dp-29/
  2. 最长公共子序列:https://www.geeksforgeeks.org/longest-common-subsequence-dp-4/
  3. Jaro 与 Jaro-Winkler 相似度:https://www.geeksforgeeks.org/jaro-and-jaro-winkler-similarity/
  4. Levenshtein 距离:https://en.wikipedia.org/wiki/Levenshtein_distance
  5. n 元语法:https://zh.wikipedia.org/wiki/N%E5%85%83%E8%AF%AD%E6%B3%95
  6. 开源工具 py_stringmatching:https://anhaidgroup.github.io/py_stringmatching/v0.4.2/index.html
  7. 开源工具 FuzzyWuzzy:https://github.com/seatgeek/fuzzywuzzy
  8. 开源工具 python-Levenshtein:https://github.com/ztane/python-Levenshtein/
  9. Cui, Yi, et al. “Finding email correspondents in online social networks.” World Wide Web 16.2 (2013): 195-218.
  10. Rieck, Konrad, and Christian Wressnegger. “Harry: A tool for measuring string similarity.” The Journal of Machine Learning Research 17.1 (2016): 258-262.

复杂网络中的节点相似性

在机器学习领域,很多时候需要衡量两个对象的相似性,特别是在信息检索,模式匹配等方向上。一般情况下,相似性与距离是都是为了描述两个对象之间的某种性质。但在实际使用的时候,需要根据具体的情况来选择合适的相似度或者距离函数。

相似性与距离

首先,我们来看一下相似性函数的含义。对于两个对象 x,y \in X, 相似性函数 s:X\times X\rightarrow \mathbb{R} 是将 X\times X 映射到实数域 \mathbb{R} 的有界函数,i.e. 存在上下界使得 s_{min}\leq s\leq s_{max}, 它具有以下两个性质:

  1. 自反性:s(x,x)=s_{max} 对于所有的 x\in X 都成立;
  2. 对称性:s(x,y)=s(y,x) 对于所有的 x,y\in X 都成立;

一般情况下,不要求相似度函数具有三角不等式的性质。相似度越大,表示两个元素越相似;相似度越小,表示两个元素越不相似。

相似度

其次,我们来看一下距离函数的含义。对于两个对象 x,y\in X, 距离函数 d:X\times X\rightarrow \mathbb{R}^{+}\cup\{0\} 是将 X\times X 映射到非负实数域的函数,它只存在下界 0, 并不存在上界,它具有以下三个性质:

  1. 自反性:d(x,x)=0 对于所有的 x\in X 都成立;
  2. 对称性:d(x,y)=d(y,x) 对于所有的 x,y\in X 都成立;
  3. 三角不等式:d(x,y)+d(y,z)\geq d(x,z) 对于所有的 x,y,z\in X 都成立。

距离越小,表示两个元素越近;距离越大,表示两个元素越远。

距离

相似度(Similarity)

对于欧式空间 \mathbb{R}^{n} 中的两个点 A=(a_{1},a_{2},\cdots,a_{n})B=(b_{1},b_{2},\cdots,b_{n}) 而言,可以多种方法来描述它们之间的相似性。

余弦相似度(Cosine Similarity)

\text{Cosine Similarity}(A,B)=\frac{A\cdot B}{||A||_{2}\cdot ||B||_{2}}=\frac{\sum_{i=1}^{n}a_{i}b_{i}}{\sqrt{\sum_{i=1}^{n}a_{i}^{2}}\cdot\sqrt{\sum_{i=1}^{n}b_{i}^{2}}}.

根据 Cauchy 不等式可以得到 Cosine Similarity 的取值范围是 [-1,1].

Pearson 相似度(Pearson Similarity)

\text{Pearson Similarity}(A,B)=\frac{cov(A,B)}{\sigma_{A}\cdot\sigma_{B}}=\frac{\sum_{i=1}(a_{i}-\overline{A})\cdot(b_{i}-\overline{B})}{\sqrt{\sum_{i=1}^{n}(a_{i}-\overline{A})^{2}}\cdot\sqrt{\sum_{i=1}^{n}(b_{i}-\overline{B})^{2}}}.

其中 \overline{A}=\sum_{i=1}^{n}a_{i}/n, \overline{B}=\sum_{i=1}^{n}b_{i}/n. 同样根据 Cauchy 不等式可以得到 Pearson Similarity 的取值范围是 [-1,1].

Pearson 系数

Dice 相似度(Dice Similarity)

\text{Dice Similarity}(A,B)=\frac{2\sum_{i=1}^{n}a_{i}b_{i}}{\sum_{i=1}^{n}(a_{i}^{2}+b_{i}^{2})},

其中 AB 不能同时是零点,并且由均值不等式可以得到 Dice Similarity 的范围也是 [-1,1].

除了欧式空间的点之外,在有的情况下需要对两个集合 AB 来做相似度的判断。特别地,欧式空间 \mathbb{R}^{n} 里面的点可以看成 n 个点所组成的集合。因此,下面的集合相似度判断方法同样适用于欧式空间的两个点。

Jaccard 相似度(Jaccard Similarity)

对于集合 AB 而言,

\text{Jaccard Similarity}=\frac{|A\cap B|}{|A\cup B|} = \frac{|A\cap B|}{|A|+|B|-|A\cap B|},

其中,|\cdot| 表示集合的势,并且 Jaccard 相似度的取值范围是 [0,1]. 越靠近 1 表示两个集合越相似,越靠近 0 表示两个集合越不相似。

两个集合 A 和 B

重叠相似度(Overlap Similarity)

对于集合 AB 而言,

\text{Overlap Similarity}=\frac{|A\cap B|}{\min\{|A|, |B|\}}

= \max\bigg\{\frac{|A\cap B|}{|A|}, \frac{|A\cap B|}{|B|}\bigg\}

= \max\{P(B|A), P(A|B)\},

其中 P(B|A), P(A|B) 指的是条件概率,意思分别是 A 发生的时候 B 同时发生的概率,B 发生的时候 A 同时发生的概率。重叠相似度的另外一个名称是 Hub Promoted(HP),它主要用于计算两个集合的重叠程度。

Overlap Similarity

类似的,可以将重叠相似度中的 min 函数换成 max 函数,那就是所谓的 Hub Degressed(HD),用公式来描述就是

\text{HD}(A,B)=\frac{|A\cap B|}{\max\{|A|,|B|\}},

它可以用于描述两个集合不重叠的程度。

距离(Distance)

欧氏距离(Euclidean Distance)

d_{2}(A,B)=\sqrt{\sum_{i=1}^{n}(a_{i}-b_{i})^{2}}.

另外,如果将 2 进行推广,则可以引导出 L^{p}(1\leq p\leq  +\infty) 距离如下:

d_{p}(A,B)=\bigg(\sum_{i=1}^{n}|a_{i}-b_{i}|^{p}\bigg)^{\frac{1}{p}}, 其中 p\geq 1.

d_{\infty}(A,B)=\max_{1\leq i\leq n}|a_{i}-b_{i}|.

欧氏距离

复杂网络中的节点相似性

在复杂网络 G=(V,E) 中,G 表示顶点集合,E 表示边的集合。为了简单起见,这里暂时是考虑无向图的场景。对于顶点 x \in V 而言,N(x) 表示其邻居的集合。在复杂网络中,同样需要描述两个顶点 x,y\in V 的相似性,于是可以考虑以下指标。

无标度网络

共同邻居相似度(Common Neighbours Similarity)

对于两个顶点 x,y\in V 而言,如果它们的共同邻居数越多,表示它们的相似度越高,反之,相似度越低。

CN(x,y)=|N(x)\cap N(y)|=\sum_{u\in N(x)\cap N(y)}1.

共同邻居

所有邻居相似度(Total Neighbours Similarity)

类似地,将顶点 xy 的邻居求并集,也可以得到一个指标,TN(x,y)=|N(x)\cup N(y)|.

Preferential Attachment

PA(x,y)=|N(x)|\cdot |N(y)|, 它将 xy 的邻居数乘起来,获得一个指标。

Jaccard 相似度(Jaccard Similarity)

如果将两个节点 xy 的邻居分别作为两个集合 N(x), N(y), J(x,y)=CN(x,y)/TN(x,y) 就可以作为顶点 xy 的 Jaccard 相似度指标,其相似度是通过邻居来衡量的。

Sorensen-Dice 相似度(Sorensen-Dice Similarity)

SI(x,y)=\frac{2|N(x)\cap N(y)|}{|N(x)|+|N(y)|},

该相似度与 Jaccard 相似度有恒等变换,J(x,y)=\frac{SI(x,y)}{2-SI(x,y)}SI(x,y)=\frac{2\cdot J(x,y)}{1+J(x,y)}.

Hub Promoted 相似度

该相似度描述了顶点 xy 的重叠程度,

HP(x,y) = \frac{|N(x)\cap N(y)|}{\min\{|N(x)|,|N(y)|\}}.

Hub Depressed 相似度

HD(x,y)=\frac{|N(x)\cap N(y)|}{\max\{|N(x)|,|N(y)|\}}.

复杂网络的社区

好友度量(Friend Measure)

\text{Friend-measure}(x,y)=\sum_{u\in N(x)}\sum_{v\in N(y)}\delta(u,v),

其中 \delta 用于判断 u,v 之间是否有边相连接。如果相连接,则取值为 1, 否则取值为 0.

Adamic Adar 相似度(Adamic Adar Similarity)

A(x,y)=\sum_{u\in N(x)\cap N(y)}\frac{1}{\ln |N(u)|},

因此,0\leq A(x,y)\leq \frac{CN(x,y)}{\ln(2)}. 事实上,当 u\in N(x)\cap N(y) 时,|N(u)|\geq 2. A(x,y) 越大,表示顶点 xy 的相似度就越高;反之,如果 A(x,y) 越小,表示顶点xy 的相似度就越低。Adamic Adar Algorithm 相当于在共同邻居的计算上增加了权重,如果 x,y 的共同邻居 u 拥有较多的邻居,则降低权重,否则增加权重。

Resource Allocation 相似度(Resource Allocation Similarity)

RA(x,y)=\sum_{u\in N(x)\cap N(y)}\frac{1}{|N(u)|},

该相似度函数与 Adamic Adar 相似度类似,只是分母上没有增加对数函数而已。

参考文献:

  1. Silva, Thiago Christiano, and Liang Zhao. Machine learning in complex networks. Vol. 2016. Switzerland: Springer, 2016.
  2. Barabási, Albert-László. Network science. Cambridge university press, 2016.
  3. Wang, Peng, et al. “Link prediction in social networks: the state-of-the-art.” Science China Information Sciences 58.1 (2015): 1-38.

随机图模型

数学家是一种把咖啡变成定理的机器。
Alfred Renyi

A mathematician is a machine for turning coffee into theorems.
Alfred Renyi

随机图的历史

在 1959 和 1968 年期间,数学家 Paul Erdos 和 Alfred Renyi 发表了关于随机图(Random Graph)的一系列论文,在图论的研究中融入了组合数学和概率论,建立了一个全新的数学领域分支—随机图论。

随机图的案例 p=0.01

随机图的定义

本文只关注无向图的场景。顾名思义,随机图(Random Graph)就是将一堆顶点随机的连接上边。好比在地上撒了一堆豆子,而豆子之间是否用线来相连是根据某个概率值来确定的。通常来说,对于随机图而言有两种定义方式

  1. 【定义一】给定 NM, G_{1}(N,M) 的定义是随机从 N 个顶点和 M 条边所生成的所有图集合中选择一个。其中,这样的图集合的势是 C(N(N-1)/2, M), 因此获得其中某一个图的概率是 1/C(N(N-1)/2, M).
  2. 【定义二】给定 NpG_{2}(N,p) 的定义是有 N 个顶点,并且两个顶点之间以概率 p\in[0,1] 来决定是否连边。

事实上,这两个定义是等价的,N 个顶点的图最多拥有的边数是 N(N-1)/2,G_{1}(N,M) 恰好有 M 条边,并且它们分配的概率是均等的,因此两个顶点之间是否存在边的概率就是 p = M/(N(N-1)/2), 这里的 C 指的是组合数。i.e.

G_{1}(N,M) = G_{2}(N, \frac{M}{N(N-1)/2}).

另一方面,对于 G_{2}(N,p) 而言,顶点两两之间是否存在边的概率是 p,N 个顶点的图最多拥有 N(N-1)/2 条边,于是边数为 pN(N-1)/2. i.e.

G_{2}(N,p)=G_{1}(N,pN(N-1)/2).

进一步地,通过以上两个公式可以得到:

G_{1}(N,M)=G_{2}(N,\frac{M}{N(N-1)/2}) = G_{1}(N,M).

在定义一中,可以直接算出所有顶点的平均度是 \langle k\rangle = 2 M /N. 但如果要计算图的其余指标,用第二种定义 G_{2}(N,p) 反而更加容易,因此后续将会重点关注第二种定义,为方便起见,记号简化为 G(N,p) = G_{2}(N,p).

随机图的度

图的度(degree)指的是对于某个顶点而言,与它相关联的边的条数。对于随机图 G(N,p) 而言,它的边数大约是 pN(N-1)/2, 最多与该节点相连接的顶点数为 N-1, 整个图的顶点平均度是(边数 * 2) / 顶点数,用记号 \langle k\rangle 来表示,意味着顶点平均度是 \langle k\rangle = p(N-1) \sim pN,N 充分大的时候成立。换言之,

p \sim \langle k\rangle / N.

顶点上的值就是该顶点的度

对于随机图 G(N,p) 中的一个顶点 i 而言,我们想计算它恰好有 d 条边的概率值。事实上,对于除了 i 之外的 N-1 个点而言,有 d 个顶点与 i 相连,N-1-d 个顶点与 i 不相连,其概率是 p^{d}(1-p)^{N-1-d}, 同时需要从这 N-1 个点中选择 d 个点,因此,顶点 i 的度恰好是 d 的概率是

p_{d}=C(N-1, d)\cdot p^{d}\cdot (1-p)^{N-1-d}.

特别地,当 d\ll N 时,上述概率近似于泊松分布(Possion Distribution)。事实上,p=\langle k\rangle / (N-1) 并且

C(N-1,d) = (N-1)(N-2)\cdots(N-d+1)/ d! \sim (N-1)^{d} / d!,

(1-p)^{N-1-d} \sim (1-\langle k\rangle /(N-1))^{N-1-d} \sim e^{-\langle k\rangle},

因此,在 d\ll N 时,p_{d} 近似于泊松分布,

p_{d} \sim \langle k\rangle^{d}e^{-\langle k\rangle}/d!.

随机图的连通分支

对于随机图 G(N, p) 而言,它的连通分支个数是与顶点的平均度 \langle k\rangle 息息相关的。特别地,当 \langle k\rangle=0 时,每个顶点都是孤立的,连通分支个数为 N;\langle k\rangle=N-1 时,任意两个顶点都有边相连接,整个图是完全图,连通分支的个数是 1. 顶点的平均度从 0N-1 的过程中,连通分支的个数从 N 演变到 1, 最大连通分支顶点数从 1 演变到 N, 那么在这个变化的过程中,最大连通分支的顶点数究竟是怎样变化的呢?是否存在一些临界点呢?数学家 Erdos 和 Renyi 在 1959 年的论文中给出了答案:

对于随机图 G(N,p) 而言,用 N_{G} 表示最大连通分支的顶点个数,那么对于图的平均度 \langle k\rangle 而言,

  1. \langle k\rangle = Np < 1, 那么 N_{G} = O(\ln(N));
  2. \langle k\rangle = Np = 1, 那么 N_{G} = O(N^{2/3});
  3. \langle k\rangle = Np \in (1, \ln(N)), 那么巨连通分支(Giant Component)存在,同时存在很多小的连通分支,在临界点 1 的附近时, N_{G} \sim (p-p_{c})N, 这里 p_{c}=1/N;
  4. \langle k\rangle = Np \in (\ln(N),+\infty), 那么图 G 是全连通图,i.e. N_{G}=N.

在这个定理中,对于顶点的平均度 \langle k\rangle 而言,存在两个临界点,分别是 1\ln(N).\langle k\rangle < 1 时,巨连通分支不存在,所有连通分支的量级都在 O(\ln(N)) 以下;当 \langle k\rangle = 1 时,巨连通分支开始出现,量级大约是 O(N^{2/3});1<\langle k\rangle <\ln(N) 时,随机图存在一个巨连通分支和很多小的连通分支;当 \langle k\rangle > \ln(N) 时,图是连通图。

图的平均度的临界点

整个定理的证明有点复杂,但本文将会介绍两个临界点的计算。先来考虑第一个临界点 \langle k\rangle = 1 的情况:

N_{G} 来表示随机图 G 中的最大连通分支的顶点个数,u 表示图 G 中不在最大连通分支的顶点比例,i.e.

u=(N-N_{G})/N = 1 - N_{G}/N= 图的顶点不在最大连通分支的概率。

对于不在最大连通分支的顶点 i 而言,其余的 N-1 个顶点分成两种情况,Case(1):要么 i 与之不相连,此时概率是 1-p; Case(2):要么 i 与之相连,但此时的顶点不能在最大连通分支中,那就只能在剩下的 uN 个顶点中,其概率是 pu. 于是,对于所有顶点而言,它不在最大连通分支的概率是 (1-p+pu)^{N-1}. 于是,

u=(1-p+pu)^{N-1}=(1-p(1-u))^{N-1}.

根据 p\sim\langle k\rangle /N\lim_{N\rightarrow +\infty}(1+x/N)^{N}=e^{x} 可以得到当 N 充分大时,有

u = (1-p(1-u))^{N-1} = (1-(1-u)\langle k\rangle /N)^{N-1} \sim e^{-(1-u)\langle k\rangle}.

s= 1-u = N_{G}/N, 它表示最大连通分支的顶点个数在所有顶点个数的占比,从而可以得到近似方程

1-s=e^{-\langle k\rangle s}.

g(s) = 1 - s - e^{-\langle k\rangle s},g(0) = 0, g(1) = -e^{\langle k\rangle}<0. 它的导数是 g'(s) = - 1 + \langle k\rangle e^{-\langle k\rangle s}, 通过计算可以得到:

  1. \langle k\rangle \leq 1 时,g'(s)<0(0,1) 上成立,i.e. g(s) = 0[0,1] 上的唯一解是 s=0, 换言之,N_{G}/N = s \rightarrow 0;
  2. \langle k\rangle > 1 时,g'(s)>0(0,\ln\langle k\rangle/\langle k\rangle) 成立,g'(s)<0(\ln\langle k\rangle /\langle k\rangle,1) 成立。换言之,g(s)=0[0,1] 上除了零之外还有解 s_{0}\in(0,1). 此时会存在巨连通分支,N_{G}/N = s_{0}\in (0,1) 是解。

因此,最大连通分支的顶点数在这个点会出现突变,1 是该方程的第一个临界点,并且是出现巨连通分支的临界点。

再来考虑第二个临界点 \langle k\rangle = \ln(N) 的情况。对于极限状况而言,假设仅有一个顶点不在最大连通分支中,那么 s = N_{G}/N = (N-1)/N, 此刻,

1/N=1-s=e^{-\langle k\rangle s}=e^{-\langle k\rangle (N-1)/N},

两边求对数可以得到 \langle k\rangle = \ln(N), 因此,\ln(N) 也是一个临界点,并且是出现全连通图的临界点。

随机图的六度分离

六度分离又称为小世界现象,它的含义是在地球上任意选择两个人,他们之间最多相隔 6 个相识关系。换言之,来自世界上任何地方的两个人都可以通过不超过 6 个相识关系所连接起来。

The Six Degrees of Larry Stone

图中两个顶点的距离定义为两个顶点之间的最短路径长度,图的直径就是图中任意两点的距离的最大值。对于随机图 G(N,p) 而言,如果 \langle k\rangle \leq 1 则是不连通的,因此通常只需要考虑 \langle k\rangle>1 的情况,甚至只考虑 \langle k\rangle >\ln(N) 的全连通图。任取一个顶点 i,则有

  1. \langle k\rangle 个距离为 1 的顶点;
  2. \langle k\rangle^{2} 个距离为 2 的顶点;
  3. \langle k\rangle^{3} 个距离为 3 的顶点;
  4. \langle k\rangle^{d_{max}} 个距离为 d_{max} 的顶点;

同时,G(N,p) 而言,顶点的个数为 N, 这意味着 \langle k\rangle + \langle k\rangle^{2}+\cdots+\langle k\rangle^{d_{max}}\leq N. 通过等比级数的公式可以得到 \langle k\rangle^{d_{max}} \leq N, 因此,

d_{max} = O(\ln(N)/\ln(\langle k\rangle)).

而随机图的直径的量级是与 d_{max} 成正比的,因此,随机图的直径量级同样是 O(\ln(N)/\ln(\langle k\rangle)). 如果 N = 10^9 并且每个人认识 \langle k\rangle = 200 个人,于是随机图的直径量级是 \ln(6*10^9) / \ln(200) = 4.25 < 6.

参考文献

  1. Erdos Renyi Model:https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model
  2. Giant Component:https://en.wikipedia.org/wiki/Giant_component
  3. Erdős P, Rényi A. On the evolution of random graphs[J]. Publ. Math. Inst. Hung. Acad. Sci, 1960, 5(1): 17-60.
  4. Albert R, Barabási A L. Statistical mechanics of complex networks[J]. Reviews of modern physics, 2002, 74(1): 47.
  5. 《巴拉巴西网络科学》,艾伯特-拉斯洛·巴拉巴西(Albert-LászlóBarabási),2020.

主动学习(active Learning)

主动学习背景介绍

机器学习的研究领域包括有监督学习(Supervised Learning)无监督学习(Unsupervised Learning)半监督学习(Semi-supervised Learning)强化学习(Reinforcement Learning)等诸多内容。针对有监督学习和半监督学习,都需要一定数量的标注数据,也就是说在训练模型的时候,全部或者部分数据需要带上相应的标签才能进行模型的训练。但是在实际的业务场景或者生产环境中,工作人员获得样本的成本其实是不低的,甚至在某些时候是相对较高的,那么如何通过较少成本来获得较大价值的标注数据,进一步地提升算法的效果就是值得思考的问题了。

机器学习

在工业界的图像标注领域,虽然有 ImageNet 这个学术界和工业界都在使用的图像数据库,但是在很多特殊的业务场景上,从业人员依旧需要想尽办法去获取业务标注数据。在安全风控领域,黑产用户相对于正常用户是偏少的,因此,如何通过极少的黑产用户来建立模型则是值得思考的问题之一。在业务运维领域,服务器,app 的故障时间相对于正常运行的时间也是偏少的,必然会出现样本不均衡的情况。因此,在这些业务领域,要想获得样本和构建模型,就必须要通过人力的参与。那么如何通过一些机器学习算法来降低人工标注的成本就是从业者需要关注的问题了。毕竟需要标注 100 个样本和需要标注成千上万的样本所需要的人力物力是截然不同的。

在学术界,同样有学者在关注这方面的问题,学者们通过一些技术手段或者数学方法来降低人们标注的成本,学者们把这个方向称之为主动学习(Active Learning)。在整个机器学习建模的过程中有人工参与的部分和环节,并且通过机器学习方法筛选出合适的候选集给人工标注的过程。主动学习(Active Learning)的大致思路就是:通过机器学习的方法获取到那些比较“难”分类的样本数据,让人工再次确认和审核,然后将人工标注得到的数据再次使用有监督学习模型或者半监督学习模型进行训练,逐步提升模型的效果,将人工经验融入机器学习的模型中。

在没有使用主动学习(Active Learning)的时候,通常来说系统会从样本中随机选择或者使用一些人工规则的方法来提供待标记的样本供人工进行标记。这样虽然也能够带来一定的效果提升,但是其标注成本总是相对大的。

用一个例子来比喻,一个高中生通过做高考的模拟试题以希望提升自己的考试成绩,那么在做题的过程中就有几种选择。一种是随机地从历年高考和模拟试卷中随机选择一批题目来做,以此来提升考试成绩。但是这样做的话所需要的时间也比较长,针对性也不够强;另一种方法是每个学生建立自己的错题本,用来记录自己容易做错的习题,反复地巩固自己做错的题目,通过多次复习自己做错的题目来巩固自己的易错知识点,逐步提升自己的考试成绩。其主动学习的思路就是选择一批容易被错分的样本数据,让人工进行标注,再让机器学习模型训练的过程。

那么主动学习(Active Learning)的整体思路究竟是怎样的呢?在机器学习的建模过程中,通常包括样本选择,模型训练,模型预测,模型更新这几个步骤。在主动学习这个领域则需要把标注候选集提取和人工标注这两个步骤加入整体流程,也就是:

  1. 机器学习模型:包括机器学习模型的训练和预测两部分;
  2. 待标注的数据候选集提取:依赖主动学习中的查询函数(Query Function);
  3. 人工标注:专家经验或者业务经验的提炼;
  4. 获得候选集的标注数据:获得更有价值的样本数据;
  5. 机器学习模型的更新:通过增量学习或者重新学习的方式更新模型,从而将人工标注的数据融入机器学习模型中,提升模型效果。
主动学习的流程

通过这种循环往复的方法,就可以达到人工调优模型的结果。其应用的领域包括:

  1. 个性化的垃圾邮件,短信,内容分类:包括营销短信,订阅邮件,垃圾短信和邮件等等;
  2. 异常检测:包括但不限于安全数据异常检测,黑产账户识别,时间序列异常检测等等。

主动学习的模型分类包括两种,第一种是流式的主动学习(Sequential Active Learning),第二种是离线批量的主动学习(Pool-based Active Learning)。在不同的场景下,业务人员可以选择不同的方案来执行。

主动学习的三种场景

而查询策略(Query Strategy Frameworks)就是主动学习的核心之处,通常可以选择以下几种查询策略:

  1. 不确定性采样的查询(Uncertainty Sampling);
  2. 基于委员会的查询(Query-By-Committee);
  3. 基于模型变化期望的查询(Expected Model Change);
  4. 基于误差减少的查询(Expected Error Reduction);
  5. 基于方差减少的查询(Variance Reduction);
  6. 基于密度权重的查询(Density-Weighted Methods)。

不确定性采样(Uncertainty Sampling)

顾名思义,不确定性采样的查询方法就是将模型中难以区分的样本数据提取出来,提供给业务专家或者标注人员进行标注,从而达到以较快速度提升算法效果的能力。而不确定性采样方法的关键就是如何描述样本或者数据的不确定性,通常有以下几种思路:

  1. 置信度最低(Least Confident);
  2. 边缘采样(Margin Sampling);
  3. 熵方法(Entropy);

Least Confident

对于二分类或者多分类的模型,通常它们都能够对每一个数据进行打分,判断它究竟更像哪一类。例如,在二分类的场景下,有两个数据分别被某一个分类器预测,其对两个类别的预测概率分别是:(0.9,0.1) 和 (0.51, 0.49)。在此情况下,第一个数据被判定为第一类的概率是 0.9,第二个数据被判定为第一类的概率是 0.51,于是第二个数据明显更“难”被区分,因此更有被继续标注的价值。所谓 Least Confident 方法就是选择那些最大概率最小的样本进行标注,用数学公式描述就是:

x_{LC}^{*}=argmax_{x}(1-P_{\theta}(\hat{y}|x))=argmin_{x}P_{\theta}(\hat{y}|x),

其中 \hat{y}=argmax_{y}P_{\theta}(y|x),这里的 \theta 表示一个已经训练好的机器学习模型参数集合。\hat{y} 对于 x 而言是模型预测概率最大的类别。Least Confident 方法考虑那些模型预测概率最大但是可信度较低的样本数据。

Margin Sampling

边缘采样(margin sampling)指的是选择那些极容易被判定成两类的样本数据,或者说这些数据被判定成两类的概率相差不大。边缘采样就是选择模型预测最大和第二大的概率差值最小的样本,用数学公式来描述就是:

x_{M}^{*}=argmin_{x}(P_{\theta}(\hat{y}_{1}|x)-P_{\theta}(\hat{y}_{2}|x)),

其中 \hat{y}_{1}\hat{y}_{2} 分别表示对于 x 而言,模型预测为最大可能类和第二大可能类。

特别地,如果针对二分类问题,least confident 和 margin sampling 其实是等价的。

Entropy

在数学中,可以使用熵(Entropy)来衡量一个系统的不确定性,熵越大表示系统的不确定性越大,熵越小表示系统的不确定性越小。因此,在二分类或者多分类的场景下,可以选择那些熵比较大的样本数据作为待定标注数据。用数学公式表示就是:

x_{H}^{*}=argmax_{x}-\sum_{i}P_{\theta}(y_{i}|x)\cdot \ln P_{\theta}(y_{i}|x),

相较于 least confident 和 margin sample 而言,entropy 的方法考虑了该模型对某个 x 的所有类别判定结果。而 least confident 只考虑了最大的概率,margin sample 考虑了最大的和次大的两个概率。

不确定性采样的差异性

基于委员会的查询(Query-By-Committee)

除了考虑单个模型的不确定性采样方法之外,还可以考虑多个模型的场景,这就是类似集成学习的方法。通过多个模型投票的模式,来选择出那些较“难”区分的样本数据。在 QBC(Query-By-Committee)的技术方案中,可以假设有 C 个模型,其参数分别是 \{\theta^{(1)},\cdots,\theta^{(C)}\},并且这些模型都是通过数据集 \mathcal{L} 的训练得到的。

如果不需要考虑每个模型的检测效果,其实可以考虑类似不确定性采样中的 least confident 和 margin sampling 方法。可以选择某一个分类器难以区分的样本数据,也可以选择其中两三个分类器难以区分的数据。但是如果要考虑所有模型的分类效果的时候,则还是需要熵(Entropy)或者 KL 散度等指标。因此,QBC 通常也包括两种方法:

  1. 投票熵(Vote Entropy):选择这些模型都无法区分的样本数据;
  2. 平均 KL 散度(Average Kullback-Leibler Divergence):选择 KL 散度较大的样本数据。

投票熵(Vote Entropy)

对于这种多模型 \{\theta^{(1)},\cdots,\theta^{(C)}\} 的场景而言,可以用熵来衡量样本数据被这些分类器区分的难易程度,如果这些分类器都把样本数据划分到某一类,则容易区分;如果分类器把样本数据划分到多类,则表示难以区分,需要重点关注。用数学公式表达就是:

x_{VE}^{*}=argmax_{x}-\sum_{i}\frac{V(y_{i})}{C}\cdot\ln\frac{V(y_{i})}{C},

其中 y_{i} 表示第 i 类,求和符号表示将所有的类别 i 相加,V(y_{i}) 表示投票给 y_{i} 的分类器个数,C 表示分类器的总数,并且 \sum_{i}V(y_{i})=C

平均 KL 散度(Average KL Divergence)

KL 散度可以衡量两个概率之间的“距离”,因此可以用 KL 散度计算出那些偏差较大的数据样本。用数学公式来描述就是:

x_{KL}^{*}=argmax_{x}\frac{1}{C}\sum_{c=1}^{C}D(P_{\theta^{(c)}}||P_{\mathcal{C}}),

其中 P_{\mathcal{C}}(y_{i}|x)=\frac{1}{C}\sum_{c=1}^{C}P_{\theta^{(c)}}(y_{i}|x) 也是概率分布,D(P_{\theta^{(c)}}||P_{\mathcal{C}}) 表示两个概率的 KL 散度。

期望模型变化(Expected Model Change)

模型变化最大其实可以选择那些使得梯度变化最大的样本数据。

期望误差减少(Expected Error Reduction)

可以选择那些通过增加一个样本就使得 loss 函数减少最多的样本数据。

方差减少(Variance Reduction)

选择那些方差减少最多的样本数据。

基于密度权重的选择方法(Density-Weighted Methods)

有的时候,某个数据点可能是异常点或者与大多数数据偏差较大,不太适合做样本选择或者区分,某些时候考虑那些稠密的,难以区分的数据反而价值更大。于是,可以在使用不确定性采样或者 QBC 方法的时候,将样本数据的稠密性考虑进去。用数学公式表示就是:

x_{ID}^{*}=argmax_{x}\phi_{A}(x)\cdot\bigg(\frac{1}{U}\sum_{u=1}^{U}sim(x,x^{(u)})\bigg)^{\beta},

在这里,\phi_{A} 表示某个不确定性采样方法或者 QBC 方法,\beta 表示指数参数,x^{(u)} 表示第 u 类的代表元,U 表示类别的个数。加上权重表示会选择那些与代表元相似度较高的元素作为标注候选集。

B 附近的点信息量会大于 A 附近的点

总结

在主动学习(Active Learning)领域,其关键在于如何选择出合适的标注候选集给人工进行标注,而选择的方法就是所谓的查询策略(Query Stategy)。查询策略基本上可以基于单个机器学习模型,也可以基于多个机器学习模型,在实际使用的时候可以根据情况来决定。整体来看,主动学习都是为了降低标注成本,迅速提升模型效果而存在的。主动学习的应用场景广泛,包括图像识别,自然语言处理,安全风控,时间序列异常检测等诸多领域。后续笔者将会持续关注这一领域的发展并撰写相关文档。

参考资料

  1. Settles, Burr. Active learning literature survey. University of Wisconsin-Madison Department of Computer Sciences, 2009.
  2. Aggarwal, Charu C., et al. “Active learning: A survey.” Data Classification: Algorithms and Applications. CRC Press, 2014. 571-605.

大数据领域的近似分析方法(一)

基数估算问题

基数估算(Cardinality Estimation),也称为 count-distinct problem,一直是大数据领域的重要问题之一。顾名思义,基数估算就是为了估算在一批数据中,它的不重复元素有多少个。

这个问题的应用场景十分广泛。例如:对于 Google 主页面而言,同一个账户可能会访问 Google 主页面多次。于是,在诸多的访问流水中,如何计算出 Google 主页面每天被多少个不同的账户访问过就是一个重要的问题。那么对于 Google 这种访问量巨大的网页而言,其实统计出有十亿 的访问量或者十亿零十万的访问量其实是没有太多的区别的,因此,在这种业务场景下,为了节省成本,其实可以只计算出一个大概的值,而没有必要计算出精准的值。

从数学上来说,基数估计这个问题的详细描述是:对于一个数据流 x_{1},x_{2}, \cdots, x_{s} 而言,它可能存在重复的元素,用 n 来表示这个数据流的不同元素的个数,i.e. n=|\{x_{1},\cdots,x_{s}\}|, 并且这个集合可以表示为 \{e_{1},\cdots,e_{n}\}. 目标是:使用 m 这个量级的存储单位,可以得到 n 的估计值 \hat{n}, 其中 m\ll n, 并且估计值 \hat{n} 和实际值 n 的误差是可以控制的。

如果是想得到精确的基数,可以使用字典(dictionary)这一个数据结构。对于新来的元素,可以查看它是否属于这个字典;如果属于这个字典,则整体计数保持不变;如果不属于这个字典,则先把这个元素添加进字典,然后把整体计数增加一。当遍历了这个数据流之后,得到的整体计数就是这个数据流的基数了。

cardinality_estimation_naive_solution
Naive Solution

这种算法虽然精准度很高,但是使用的空间复杂度却很高。那么是否存在一些近似的方法,可以估算出数据流的基数呢?其实,在近几十年,不少的学者都提出了很多基数估算的方法,包括 LogLog,HyperLogLog,MinCount 等等。下面将会简要的介绍一下这些方法。

cardinality_estimation_survey_table_1
基数估计的部分算法

HyperLogLog 的理论介绍

HyperLogLog 是大数据基数统计中的常见方法,无论是 Redis,Spark 还是 Flink 都提供了这个功能,其目的就是在一定的误差范围内,用最小的空间复杂度来估算一个数据流的基数。

Spark
Spark 的 Logo

HyperLogLog 算法简要思路是通过一个 hash 函数把数据流 \mathcal{D} 映射到 \{0,1\}^{\infty}, 也就是说用二进制来表示数据流中的元素。每一个数据流中的元素 x 都对应着一个 0,1 序列。

在介绍 HyperLogLog 之前,我们可以考虑这个实际的场景。在一个抛硬币的场景下,假设硬币的正面对应着 1, 硬币的反面对应着 0; 依次扔出 0,0,0,1 的概率是多少?通过概率计算可以得到是这个概率是 1/2^{4}=1/16. 那么相当于平均需要扔 16 次,才会获得 0001 这个序列。反之,如果出现了 0001 这个序列,说明起码抛了 16 次硬币。

考虑这样一个 0,1 序列,w=w_{1}w_{2}\cdots, w_{i}\in\{0,1\}, i\geq 1,k 表示第一个 1 出现的位置。也就是说 w_{1}=w_{2}=\cdots=w_{k-1}=0. 那么在扔硬币的场景下,出现这样的序列平均至少需要扔 2^{k} 次。对于一批大量的随机的 0,1 序列, 可以根据第一个 1 出现的位置来估算这批 0,1 序列的个数。也就是说:

  • 出现序列 1XXXXX 意味着不重复的元素估计有 2^1=2 个;
  • 出现序列 01XXXX 意味着不重复的元素估计有 2^2=4 个;
  • 出现序列 001XXX 意味着不重复的元素估计有 2^3=8 个;
  • 出现序列 0001XX 意味着不重复的元素估计有 2^4=16 个。

于是,对于随机的 0,1 序列,可以定义函数 \rho(w_{1}w_{2}\cdots) 来表示 1 出现的第一个位置。i.e. \rho(1XXXXX)=1, \rho(01XXXX)=2, \rho(001XXX)=3, \rho(0001XX)=4.

简单来看,其实 HyperLogLog 的基数统计就使用了这样的思想,通过二进制中 1 出现的第一个位置来估算整体的数量。首先把这批元素通过 hash 函数处理成 0,1 序列,然后把这批 0,1 序列都放入 1 个桶,然后通过计算这个桶里面所有 0,1 序列的 \rho(w) 的最大值,就可以预估出整体的数量。i.e. M=\max_{w}\rho(w), 整体的数量预估是 2^{M}=2^{\max_{w}\rho(w)}.

  • 1 个桶:计算出 M=\max_{w}\rho(w), 预估不重复的元素个数是 2^{M}.

那么如果只有 1 个桶,其实是会存在一定的偏差的。为了解决这个问题,一种想法就是重复以上操作,从 hash 函数开始处理成 0,1 序列,每次都把这批 0,1 序列放入 1 个桶,每次获得一个 M 值。总共操作 m 次,第 j 次操作得到的值记为 M_{j}; 于是就可以对 \{M_{1},\cdots,M_{m}\} 进行均值处理,可以使用以下方法:

  • 算术平均数:M=\sum_{j=1}^{m}M_{j};
  • 几何平均数:M=\sqrt[m]{M_{1}\cdots M_{m}};
  • 调和平均数:M=m/\sum_{j=1}^{m}M_{j}^{-1};
  • 中位数:M = median\{M_{1},\cdots,M_{m}\}.

从而可以预估整体的数量为 2^{M}.

如果按照以上的步骤进行操作,就是需要重复进行多次操作,在足够多的情况下,其实是没有必要那么操作的。HyperLogLog 也是用了多个桶,但是用了一个截断的技巧。对于一个 0,1 序列 x=\cdots x_{b+2}x_{b+1}x_{b}\cdots x_{1}, HyperLogLog 从某个位置 b 开始,低位 x_{b}\cdots x_{1} 用于决定桶的序号,也就是第几个桶。桶的个数就是 m=2^{b}, 高位 \cdots x_{b+2}x_{b+1} 用于估算放在桶里面的元素个数。

每次都可以获得一个值,也就是桶里面第一次出现

  • 第 1 个桶:计算出 M_{1}=\max_{w}\rho(w), 预估元素个数 2^{M_{1}};
  • 第 2 个桶:计算出 M_{2}=\max_{w}\rho(w), 预估元素个数 2^{M_{2}};
  • ….
  • m 个桶:计算出 M_{m}=\max_{w}\rho(w), 预估元素个数 2^{M_{m}};

均值的计算,HyperLogLog 使用了调和平均数 m/\sum_{j=1}^{m}2^{-M_{j}} 来估算桶里面的元素个数,那么在有 m 个桶的情况下,整体的元素个数就可以估算为 E=m^{2}\cdot\bigg(\sum_{j=1}^{m}2^{-M_{j}}\bigg)^{-1}.

hyperloglog_algorithm_simple
原始的 HyperLogLog 算法

其中的 \alpha_{m}=\bigg(\int_{0}^{+\infty}\bigg(\log_{2}\bigg(\frac{2+u}{1+u}\bigg)\bigg)^{m}du\bigg)^{-1}.m=1 的时候,\int_{0}^{+\infty}\log_{2}\bigg(\frac{2+u}{1+u}\bigg)du 是发散的;当 m\geq 2 的时候,\int_{0}^{+\infty}\bigg(\log_{2}\bigg(\frac{2+u}{1+u}\bigg)\bigg)^{m}du 是收敛的。因此,在使用这个算法的时候最好放入 m\geq 2 个桶。

HyperLogLog 分成两块,第一块就是 add 模块,用于分桶和统计;

for v in M:
    set x := h(v);
    set j = 1 + <x(b),...,x(2),x(1)>;
    set w := x(b+1)x(b+2)...; 
    set M[j] := max(M[j], \rho(w));

在 HyperLogLog 算法中,对于集合 \mathcal{M} 中的每一个元素 v\in\mathcal{M}, 可以通过 hash 函数转换成一个 0,1 序列 h(v)=x=<\cdots x_{b+2}x_{b+1}x_{b}\cdots x_{1}>, 其中 x_{1} 表示二进制中的最低位,x_{2} 表示次低位。

然后可以通过 <x_{b}\cdots x_{1}>_{2} 来计算放在第 j 个桶,这里 j=1+<x_{b}\cdots x_{1}>_{2}. 同时将 x 的高位拿出来,也就是 x_{b+1}x_{b+2}\cdots, 计算这批序列的 \rho 函数的最大值,然后记为 M_{j}; 这一步也可以称为 merge 模块,也就是进行更新合并。

compute Z := (\sum_{j=1}^{m}2^{-M[j]})^{-1};

上一步就是 HyperLogLog 的另外一步,count 模块,于是,进一步估算出 E=\alpha_{m}m^{2}Z.

HyperLogLog 的空间复杂度特别低,大约是 O(m\log_{2}\log_{2}n) 这个量级的,其中 m 是桶的个数,n 是基数。HyperLogLog 的时间复杂度则是 O(n), 只需要遍历一遍所有元素即可得到最终结果。

  • 假设基数为 2^{k}, 二进制就是 k 位,1 最晚就会出现在第 k 个位置上;而 k 只需要 log_{2}k 个 bit 就能够存储;
  • 假设基数为 2^{64}, 二进制就是 64 位,1 最晚会出现在第 64 个位置上;而 64 需要 6 个 bit 就可以存储。

hyperloglog_algorithm_compare
算法对比

在论文中,论文 “Hyperloglog: the analysis of a near-optimal cardinality estimation algorithm” 的作者们针对各种算法进行了对比,其实 HyperLogLog 的空间复杂度是非常小的,并且误差也在可控的范围内。

hyperloglog_algorithm_theorem
HyperLogLog 的定理证明

在论文 “Hyperloglog: the analysis of a near-optimal cardinality estimation algorithm” 作者们得到上述定理,精准的给出了 HyperLogLog 算法的误差估计。因此,HyperLogLog 算法其实是有数学定理证明的。

以上只是获得了理论上的 HyperLogLog 算法,但是在实战中,其实是需要进行微调的。主要的微调部分是根据理论中的 E 值来进行调整。将 E 值进行调整的话,情况可以分成三种:

  • 小范围;
  • 中等范围;
  • 大范围;

hyperloglog_algorithm_practical
实战中的 HyperLogLog

Case(1):小范围

在小范围的情况下,E\leq 5m/2, 此时的基数相对于桶的数量而言不算太多,因此可能存在多个空桶,需要进行调整。

可以思考这样一个问题:假设有 m 个桶,同时有 n 个球,把这 n 个球随机往这 m 个桶里面扔,每个球只能够进入一个桶,那么空桶个数的期望是多少个?

Answer:假设 A_{1},\cdots,A_{m}m 个桶,

P(A_{j}=\emptyset)=\bigg(1-\frac{1}{m}\bigg)^{n} 表示桶 A_{j} 空的概率;

P(A_{j}=\emptyset \cap A_{k}=\emptyset)=\bigg(1-\frac{2}{m}\bigg)^{n} 表示桶 A_{j}, A_{k} 同时为空的概率(j\neq k);

那么空桶个数的期望就是 m\cdot\bigg(1-\frac{1}{m}\bigg)^{n},m,n 充分大的时候,约为 me^{-n/m} 个。

因此,在小范围的情况下,如果空桶的个数 V\neq 0, 那么可以更新为 m\ln(m/V). 事实上,可以通过 V=me^{-n/m} 解出 n=m\ln(m/V).

Case(2):中等范围

E 值不作调整。

Case(3):大范围

E>2^{32}/30, 那么更新为 E^{*}=-2^{32}\ln(1-E/2^{32}), 其中 E^{*}>E.

通过这样的方法,E^{*} 的误差大约在 \pm 1.04/\sqrt{m} 左右。

除此之外,\alpha_{m} 其实也可以用近似值来代替,毕竟如下公式的计算是有一定的成本的。

\alpha_{m}=\bigg(\int_{0}^{+\infty}\bigg(\log_{2}\bigg(\frac{2+u}{1+u}\bigg)\bigg)^{m}du\bigg)^{-1}.

近似的值为 \alpha_{16}=0.673, \alpha_{32}=0.697, \alpha_{64}=0.709, \alpha_{m}=0.7213/(1+1.079/m)m\geq 128.

HyperLogLog 的案例分析

有一个关于 HyperLogLog 的 demo 网站可以看到 HyperLogLog 的算法过程,其链接是 http://content.research.neustar.biz/blog/hll.html

在这个 demo 中,作者对比了 LogLog 和 HyperLogLog 的区别和运行过程,有助于大家理解整个过程。其中 LogLog 与 HyperLogLog 的区别就在与它们平均值的处理方式不一样,前者是使用算术平均值,后者是使用调和平均值。

  • LogLog\alpha_{m}\cdot m\cdot 2^{\sum_{j=1}^{m}M_{j}/m};
  • HyperLogLog\alpha_{m}\cdot m^{2}\cdot\bigg(\sum_{j=1}^{m}2^{-M[j]}\bigg)^{-1};

hyperloglog_demo_1
HyperLogLog Demo:初始化

hyperloglog_demo_2
第一个 hash 值:3852172429

3852172429 的二进制是:11100101100110110111110010001101,可以划分为100 110110111110010 001101。最后的六位是 001101,十进制就是 13,那么这个数字就会被放入第 13 个桶;而 110110111110010(从低位到高位看),\rho 函数的值就是 2;于是在第 13 个桶就会把 0 更新成 2。

hyperloglog_demo_3
第二个 hash 值:2545698499

2545698499 的二进制是 10010111101111000100011011000011,用同样的分析可得结论。

hyperloglog_demo_4
第三个 hash 值:2577699815

2577699815 的二进制是 10011001101001001001001111100111。

hyperloglog_demo_5
第四个 hash 值:775376803

775376803 的二进制是 101110001101110100111110100011。

hyperloglog_demo_6
运行结束

从以上的 Demo 运行过程可以看出,整个 HyperLogLog 的算法逻辑还是相对清晰的,其整个算法的亮点应该在于借助了抛硬币的场景,用抛硬币的结果来估算抛硬币的次数。

参考文献:

  1. Count-distinct Problem 的维基百科:https://en.wikipedia.org/wiki/Count-distinct_problem
  2. Heule, Stefan, Marc Nunkesser, and Alexander Hall. “HyperLogLog in practice: algorithmic engineering of a state of the art cardinality estimation algorithm.” Proceedings of the 16th International Conference on Extending Database Technology. 2013.
  3. Flajolet, Philippe, et al. “Hyperloglog: the analysis of a near-optimal cardinality estimation algorithm.” 2007.
  4. HyperLogLog 的 demo 网站:http://content.research.neustar.biz/blog/hll.html

时间序列异常检测—节假日效应的应对之道

在时间序列异常检测中,通常有一个较为常见的场景就是“节假日效应”。所谓节假日效应,指的就是在节假日的时候,其时间序列的走势跟日常有着明显的差异性,但是又属于正常的情况。从国内 2020 年的节假日安排可以看出,一年中有好几个关键的假日:

  1. 元旦:1 天;
  2. 春节:7 天;
  3. 清明节:3 天;
  4. 五一劳动节:5 天;
  5. 端午节:3 天;
  6. 国庆节:8 天。

在这些节假日的时候,为了调休,自然也会带来工作日上的调整。例如:在 2020 年 1 月 19 日,2020 年 2 月 1 日是需要上班的(虽然今年受疫情影响最终也没上班)。因此,在这些节假日进行调整和变化的时候,各种各样的业务指标(时间序列)通常也会发生变化,变得跟以往的走势不太一致。因此,如何解决节假日效应的时间序列异常检测就是业务上所面临的问题之一。

£¨Í¼±í£©[Éç»á]2020Äê½Ú¼ÙÈշżٰ²ÅŹ«²¼
2020 年的放假安排
清华大学的 Netman 实验室在 2019 年发表了一篇论文,专门用于解决时间序列异常检测中的节假日效应问题,论文的标题是《Automatic and Generic Periodic Adaptation for KPI Anomaly Detection》。在本文中,所用的时间序列是关于各种各样的业务指标的,包括搜索引擎,网上的应用商店,社交网络数据等等。作者们针对 KPI(Key Performance Indicator)做了时间序列异常检测,并且发明了一种方法来避免节假日效应的问题。论文针对时间序列的工作日(work days),休息日(off days),节假日(festival)做了必要的区分,然后将时间序列的不同时间段进行合理地拆分和组装,再进行时间序列异常检测,从而在一定的程度上解决节假日效应问题。

节假日效应_Fig1
实际案例(1)

在实际的案例中,我们可以看到,同一条时间序列的走势在工作日(work day),休息日(off day),春节(Spring Festival)明显是不一样的。因此,根据工作日的时间序列走势来预测春节的走势明显是不太合理的;同理,根据春节的走势来预测休息日的走势也会带来一定的偏差的。那么如何解决节假日效应的问题就成为了本篇论文的关键之一。

节假日效应_Fig3
实际案例(2)

在上图中,我们可以看到论文中使用的数据都具有某种周期性(Periodicity)。KPI A,B,C 都是具有明显具有工作日和周末特点的,在工作日和周末分别有着不同的形状;KPI D 则是关于网上应用商店周五促销的,因此在周五周六的时候,其实时间序列会出现一个尖峰(peak);KPI E 的话则是每隔 7 天,会有两个尖刺,然后并且迅速恢复;KPI F 的话则是可以看出时间序列在十一的走势跟其余的时间点明显有区别。除此之外,对于一些做旅游,电商等行业的公司,其节假日效应会更加突出一点,而且不同的业务在节假日的表现其实也是不一样的。有的时间序列在节假日当天可能会上涨(电商销售额),有的时间序列在节假日当天反而会下降(订车票,飞机票的订单量)。因此,在对这些时间序列做异常检测的同时,如何避免其节假日效应就是一个关键的问题了。

而在实际处理的时候,通常也会遇到几个常见的问题;

  1. 周期性的多样性:通过实际案例可以看出,对于不同的时间序列,其周期是完全不一样的,而且在不同的周期上也有着完全不同的表现;
  2. KPI 数量巨大:这个通常来说都是智能运维领域中的常见问题;
  3. 周期的漂移:一般来说,通过时间序列的走势我们只能够看出一个大致的变化,但是具体到细节的话,周期是存在一定的波动的。例如不一定恰好是 7 天,有可能是 7 天加减 5 分钟之类的周期。这个跟业务的具体场景有关系,也跟当时的实际情况有关。

于是,基于这些挑战,作者们希望提出一个健壮的机器学习算法来解决这个问题,本文的系统被作者们称之为 Period,正好也象征着解决节假日效应这个寓意。

Period_Fig7
Period 的整体架构

从论文中可以看出 Period 的整体架构如上图所示,包括两个部分:

  1. 离线周期性检测(offline periodicity detection);
  2. 在线适应性异常检测(online anomaly detection adaptation)。

在第一部分,每一条时间序列都会被按天切分成很多子序列(subsequence),然后将其聚集起来,把相似的时间序列放在一类,不相似的放在另外一类;在第二部分,新来的时间序列会根据其具体的日期,分入相应的聚类,然后用该类的时间序列异常检测方法来进行异常检测。

Period_Fig6
Period 的核心思路

从上图可以看到 Period 的核心思路(core idea)。在本文使用的数据中,时间序列的长度较长,一般来说都是好几个月到半年不等,甚至更长的时间。对于一条时间序列(a given KPI),可以将它的历史数据(historical data)进行按天切分,获得多个子序列(sub KPIs)。对于这多个子序列,需要进行聚类以得到不同类别。或者按照日历直接把时间序列的工作日(work day),休息日(off day),春节(spring festival)序列进行切分,将工作日放在一起,休息日放在一起,春节放在一起。把这些子序列进行拼接就可以得到三条时间序列数据,分别是原时间序列的工作日序列(work day subsequence),休息日序列(off day subsequence),春节序列(spring festival subsequence)。然后分别对着三条时间序列训练一个异常检测的模型(例如 Holt-Winters 算法,简写为 HW)。对于新来的时间序列,可以根据当日具体的日期(工作日,休息日或者春节)放入相应的模型进行异常检测,从而进一步地得到最终的结果。

在离线周期性检测的技术方案里面,是需要对时间序列进行周期性检测(Periodicity Detection)。而周期性检测有多个方案可以选择。第一种就是周期图方法(Periodogram),另外一种就是自相关函数(Auto-correlation function)。但是在这个场景下,用这些方法就不太合适了。作者们提出了别的解决方案。

在本文中,作者们提出了一种 Shape-based distance(SBD)的方法,针对两条时间序列 X=(x_{1},x_{2},\cdots,x_{m})Y=(y_{1},y_{2},\cdots,y_{m}),提出了相似性的计算方法。

X_{(s)}=\begin{cases}(0,\cdots,0,x_{1},\cdots,x_{m-s}), &\text{ if } s\geq 0 \\ (x_{1-s},x_{1-s+1},\cdots,x_{m},0,\cdots,0), &\text{ else } s<0.\end{cases}

其中 0 的个数都是 |s|. 进一步可以定义,当 s\in[-w,w]\cap\mathbb{Z} 时,

CC_{s}(X,Y)=\begin{cases}\sum_{i=1}^{m-s}x_{i}\cdot y_{s+i}, &\text{ if } s\geq 0 \\ \sum_{i=1}^{m+s}x_{i-s}\cdot y_{i}, &\text{ else } s<0.\end{cases}

于是,选择令 CC_{s}(X,Y) 归一化之后的最大值作为 X,Y 的相似度,i.e.

NCC(X,Y)=\max_{s\in[-w,w]\cap\mathbb{Z}}\frac{CC_{s}(X,Y)}{\|x\|_{2}\cdot\|y\|}.

而基于 SBD 的距离公式则可以定义为:

SBD(X,Y) = 1-NCC(X,Y).

Period_Fig8
Periodicity Drift

那么为什么需要考虑一个漂移量 s 呢,因为在一些实际的情况下,时间序列是会存在漂移的,例如上图所示。该时间序列在 10 月 30 日,31 日,11 月 1 日 都出现了一个凸起,但是如果考虑它的同比图,其实是可以清楚地看出该时间序列就存在了漂移,也就是说并不是在一个固定的时间戳就会出现同样的凸起,而是间隔了一段时间。这就是为什么需要考虑 s 的由来。

Period_Alg1
聚类的命名算法

通过相似性和距离的衡量工具,我们可以将时间序列进行聚类,然后通过上述算法也可以对每一个聚类的结果进行命名。

Period_Table1
实验数据

Period_Table2
实验数据的聚类结果

在本文中,针对以上六条时间序列,作者们做了详细的分析,也对其余的 50 条时间序列进行了实验。其使用的方法包括 HW,TSD,Diff,MA,EWMA,Donut。在 HW 中,针对不同的日期使用了不同的方法,例如 HW-day,HW-week,HW-period;其余的方法也是针对不同的日期来做的。

Period_Table3
实验方法

Period_Table12
实验效果

从实验效果来看,Period 方法的话相对于其他方法有一定的优势。

结论:Period 方法包括两个部分,第一部分是离线周期性检测,第二部分是在线适应性异常检测。通过这样的方法,可以有效地减缓时间序列异常检测受节假日效应的影响。除此之外,想必未来也会有其余学者提出相应的问题和解决方案,敬请期待。

FluxRank: 如何快速地进行机器故障定位

在运维领域,服务侧的异常会由多方面的原因造成,有的时候是因为网络的抖动,有的时候是因为机器的故障,有的时候甚至是因为人为的变更。本篇博客会介绍一种机器异常定位的方法,论文是来自于清华 Netman 实验室的《FluxRank:A Widely-Deployable Framework to Automatically Localizting Root Cause Machines for Software Service Failure Mitigation》。本篇论文主要介绍了如何从服务的故障定位到局部异常的机器,也就是说在发现服务故障的同时,进一步推断出是由哪些机器出现问题而导致的。

通常来说,在服务异常(例如服务的耗时长,失败数上涨)的时候,需要运维人员通过历史上的经验迅速定位到是哪个业务,哪个模块,甚至哪台服务器出现了故障。而人工定位的速度总是会出现瓶颈的,无论对模块的判断,还是机器的判断,都依赖于人工所积累的经验。而每个人的经验却各不相同,并且经验的传承也需要一定的时间成本。那么如何基于人工运维的经验来构建模型,进一步地提升异常定位的速度就是智能运维的关键之处之一。

FluxRank_fig_1
从告警到故障恢复

对于一条业务指标(时间序列)而言,大多数情况下是处于正常的状态(normal)。但是如果出现了错误的变更,发布了错误的程序,或者服务器突然出现了故障,都会导致业务指标出现变化,就从正常(normal)变成异常(abnormal)。这个时候就会出现一个故障的开始时间,也就是 failure start time T_{f},这个时间戳是运维领域非常重要的时间戳,它由异常检测(anomaly detection)产生,无论在告警收敛(alarm convergence)还是根因分析(root cause analysis)都非常依赖这个时间戳。而另外一个时间戳虽然没有故障开始时间那么重要,但是也有着其实用价值,那就是缓和开始时间(mitigation start time),它表示故障虽然还没有恢复,但是出于稍微平稳的走势,并没有持续恶化。在出现了故障之后,通常都会发送相应的告警给运维人员,那么在发送告警的时候,如果将异常定位的结果随之带出,则会大大减少运维人员排障的时间。在故障缓和的时间内,运维人员通常需要进行必要的操作来排查故障,例如切换流量(switch Traffic),回滚版本(Rollback Version),重启实例(Restart Instances),下线机器等操作。除此之外,为了定位问题(Root Cause Analysis),运维人员需要分析源码(Code Analysis),查看日志(Log Analysis)等一系列操作。如果能够将这一系列操作融入相应的机器学习模块中,将会节省运维人员大量的排障时间。

贝叶斯网络

通常来说,故障定位也称为根因分析或者根源分析(Root Cause Analysis),都是为了排查产生这次故障的原因。在机器学习领域,为了进行因果分析(Causal Analysis),则需要使用相应的模型来进行建模。其中较为经典的统计分析方法则是贝叶斯分析法,其中的贝叶斯网络(Bayesian Network)则是经典模型之一。下面来看一个简单的例子。

假设降雨(Rain)的概率是 0.2,不降雨的概率是 0.8;而洒水器(Sprinkler)是否开启会受到降雨的影响,其条件概率与下图所示。而降雨或者洒水器都会导致草湿润(Grass Wet),其概率分布如下图所示。那么可以问如下问题:

  1. 如果草已经湿润,求降雨的概率是多少?
  2. 如果草已经湿润,求没有降雨且洒水器开启的概率是多少?

BayesianNetwork_2
贝叶斯网络的经典案例

而这一类的问题可以通过贝叶斯公式来进行解答。从表格来看:

从 Rain 的表格可得:P(R=T)=0.2, P(R=F)=0.8

从 Rain 和 Sprinkler 的表格可得:P(S=T|R=F)=0.4, P(S=F|R=F)=0.6P(S=T|R=T)=0.01, P(S=F|R=T)=0.99

从 Grass Wet 和 Sprinkler,Rain 的表格可得:P(W=T|S=F, R=F)=0.0, P(W=F|S=F,R=F)=1.0P(W=T|S=F,R=T)=0.8, P(W=F|S=F,R=T)=0.2P(W=T|S=T,R=F)=0.9, P(W=F|S=T,R=F)=0.1P(W=T|S=T,R=T)=0.99, P(W=F|S=T,R=T)=0.01.

针对问题 1,需要计算条件概率 P(R=T|W=T)。从 Bayes 公式可以得到:P(R=T|W=T) = P(R=T,W=T)/P(W=T)。分别计算分子分母即可:

P(R=T,W=T)=P(R=T,S=T,W=T)+P(R=T,S=F,W=T)

= P(W=T|R=T,S=T)P(S=T|R=T)P(R=T) + P(W=T|R=T,S=F)P(S=F|R=T)P(R=T)

= 0.99*0.01*0.2+0.8*0.99*0.2=0.16038

P(W=T)=P(W=T,S=T,R=T)+P(W=T,S=F,R=T)+P(W=T,S=T,R=F)+P(W=T,S=F,R=F)

=P(W=T|S=T,R=T)P(S=T|R=T)P(R=T) + P(W=T|S=F,R=T)P(S=F|R=T)P(R=T) + P(W=T|S=T,R=F)P(S=T|R=F)P(R=F)+P(W=T|S=F,R=F)P(S=F|R=F)P(R=F)

= 0.99*0.01*0.2+0.8*0.99*0.2+0.9*0.4*0.8+0.0*0.6*0.8=0.44838,

那么如果草已经湿润,求降雨的概率是 P(R=T|W=T)=P(R=T,W=T)/P(W=T)=0.16038/0.44838=0.3577.

另外一个题目可以用类似的方法进行求解,在此不再赘述。

虽然贝叶斯算法能够计算出条件概率,例如本次故障是由哪些原因导致的,但是这个需要长期收集数据,需要对历史数据进行积累,才能通过人工或者统计的方法得到以上表格的条件概率。但是在实际的环境中是较难获取这些数据的,需要大数据平台的支持,因此需要探索其他的解决方案。

FluxRank

在本论文中,为了克服贝叶斯网络模型中的一些问题,针对子机异常定位的场景,设计了一套技术方案,作者们称之为 FluxRank。

FluxRank_fig_2
FluxRank 的整体框架

FluxRank 这一模块的触发需要服务指标(Service KPI)的异常,因此需要对服务指标(Service KPI)进行异常检测。这里的服务指标通常指的是业务指标,包括某块 APP 的在线人数,某个接口的成功率,某个视频网站的卡顿数等指标。当服务指标出现了异常的时候,就启动 FluxRank 模块进行异常机器定位。

如果按照人工处理的流程来看,分成几个步骤:

  1. 异常检测部分:通过设定阈值或者某个简单的规则来进行异常检测,包括服务的 KPI(Service KPI)和机器的 KPI(machine KPIs);
  2. 手工检查异常的时间段,并且查看在异常的时间段内发生了什么情况;
  3. 运维人员根据自身的业务经验来对机器的故障程度做人工排序;
  4. 运维人员根据自身的业务经验来对故障进行处理,并且人工给出处理方案。

那么 FluxRank 所面临的挑战就有以下几点:

  1. 如何衡量海量 KPIs 的变化程度?在这里不仅有服务的 KPIs,还有机器的 KPIs。而机器的 KPIs 包括内存,硬盘,IO,CPU等诸多固定的指标,那么如何对这些海量的 KPI 曲线进行变化程度的衡量,为后续的指标排序做准备就成为了一个难点;
  2. 如何对 KPIs 进行异常性或者重要性的聚类,让运维人员能够一眼看出每个聚簇的差异或者异常程度?
  3. 如何对 KPIs 聚类的结果进行排序?

为了解决以上的问题,FluxRank 的框架有以下几个贡献点:

  1. 基于 Kenel Density Estimation 用于衡量海量 KPIs 在某一个时间段的变化程度和异常程度;
  2. 基于上一步生成的异常程度,对诸多机器所形成的特征使用距离公式或者相似度公式,然后使用 DBSCAN 聚类算法来对机器进行聚类;
  3. 在排序部分,对上一步的机器聚类结果进行排序;

Change Quantification

首先,来看一下 Change Quantification 是怎么样做出来的。这里的 Change Quantification 使用与衡量机器 KPIs 的变化程度,称之为 change degree。Change degree 可以用于 CPU,内存,IO 等诸多机器指标。为了达到衡量变化程度,需要一个非常重要的信息,那就是变化的开始时间,change start time,也就是说在哪个时刻时间序列开始出现了变化。于是在 Change Quantification 部分,就分成两部分:(1)用 absolute derivative 或者 CUSUM 算法获得变化开始时间(change start time);(2)用 Kernel Density Estimation(KDE)来计算变化程度(change degree)。

FluxRank_fig_1
重要的时间戳

正如上图所示,针对服务 KPIs(ervice KPIs),存在两个关键的时间点,那就是失败开始时间(Failure Start Time)T_{f} 和缓和开始时间(Mitigation Start Time)T_{m}。在失败开始时间 T_{f} 之前,可能有的机器已经出现了故障,因此变化开始时间(Change Start Time)T_{c} 小于或者等于 T_{f}。通常情况下,一个或者多个机器故障会在半小时(30 mins)甚至更短的时间内引发服务故障,因此,只需要假设 w_{1}=30 即可。关键时间点的排序为 T_{f}-w_{1}<T_{c}\leq T_{f}<T_{m}

对于服务 KPIs 的异常检测,FluxRank 中提到了两种方法:分别是 absolute derivative 和 CUSUM 方法。

  1. absolute derivative 方法:个人理解就是对时间序列进行一阶差分操作,然后对一阶差分来做时间序列异常检测,例如 3-sigma 等方法,一旦有明显的变化,就说明当前的时间点出现了突增或者突降;与该方法比较类似的一种方法是:MAD(Median Absolute Deviation)。对于一条时间序列 X=[x_{1},\cdots,x_{n}] 而言,MAD 定义为 MAD = median_{1\leq i\leq n}(|x_{i}-median(X)|),而每个点的异常程度可以定义为:s_{i}=(x_{i}-median(X))/MAD = (x_{i}-median(X))/median_{1\leq i\leq n}(|x_{i}-median(X)|).s_{i} 较大或者较小的时候,表示上涨或者下降的异常程度。通过设置相应的阈值,同样可以获得时间序列的异常开始时间。
  2. CUSUM 算法也是用于时间序列异常检测的。对于一条时间序列 X=[x_{1},x_{2},\cdots,x_{n}],可以预估它的目标值(target value)\mu_{0},通常可以用均值来估计,也需要计算出这条时间序列的标准差 \sigma。通常设定 \mu_{1}=\mu_{0}+\delta\sigmaK=\delta\sigma/2=|\mu_{1}-\mu_{0}|/2。而 Tabular CUSUM 指的是迭代公式 C_{i}^{+}=\max[0,x_{i}-(\mu_{0}+K)+C_{i-1}^{+}]C_{i}^{-}=\max[0,(\mu_{0}-K)-x_{i}+C_{i-1}^{-}],初始值是 C_{0}^{+}=C_{0}^{-}=0。当累计偏差 C_{i}^{+} 或者 C_{i}^{-} 大于 H=5\sigma 的时候,表示 x_{i} 出现了异常,也就是 out of control。通过这个值,可以获得时间序列开始异常的时间。

从论文的描述来看,作者是使用 absolute derivative 来做异常检测的,并且定位其异常开始时间的准确率较高。

Change Degree

其次,我们来看一下变化程度(Change Degree)是怎么计算出来的,通过之前的计算,我们已经可以获得一些关键的时间戳,例如 T_{f}, T_{c}, T_{m} 等时间戳。根据变化开始时间(change start time)T_{c},同样需要设置一个窗口值 w_{2},例如 60 分钟(1 小时)。可以从两个时间段获取数据,正常时间段 [T_{c}-w_{2},T_{c}),异常时间段 [T_{c},T_{m}],分别获取到数据 \{x_{i}\}\{x_{j}\},前者是在变化开始时间之前的数据点,后者是在变化开始之后的数据点。于是,作者们通过概率值来计算变化程度 P(\{x_{j}\}|\{x_{i}\}),意思就是计算一个条件概率,在观察到 \{x_{i}\} 之后,得到 \{x_{j}\} 的概率值。

为了计算以上概率值,需要简化模型,因此这里需要假设 \{x_{j}\} 是独立同分布(iid)的,于是 P(\{x_{j}\}|\{x_{i}\})=\prod_{j=1}^{\ell}P(x_{j}|\{x_{i}\}),在这里 \ell 表示集合 \{x_{j}\} 的元素个数。 为了分别得到其上涨和下降到概率,则需要计算:

P_{o}(\{x_{j}\}|\{x_{i}\}) = \prod_{j=1}^{\ell}P(X\geq x_{j}|\{x_{i}\}),

P_{u}(\{x_{j}\}|\{x_{i}\}) = \prod_{j=1}^{\ell}P(X\leq x_{j}|\{x_{i}\}),

其中 P_{o}(\{x_{j}\}|\{x_{i}\}) 表示上涨的程度,P_{u}(\{x_{j}\}|\{x_{i}\}) 表示下降的程度。如果不想处理连乘的话,则需要处理连加:

o=-\frac{1}{\ell}\sum_{j=1}^{\ell}\ln P(X\geq x_{j}|\{x_{i}\}),

u =-\frac{1}{\ell}\sum_{j=1}^{\ell}\ln P(X\leq x_{j}|\{x_{i}\}).

在这里,作者们使用了三种概率分布函数,分别是 Beta 分布(Beta distribution),泊松分布(Poisson distribution),高斯分布(Gaussian distribution)。

Beta 分布的概率密度函数(probabilisty density function)是 f(x;\alpha,\beta) = x^{\alpha-1}(1-x)^{\beta-1}/B(\alpha,\beta),其中 B(\alpha,\beta)=\Gamma(\alpha)\Gamma(\beta)/\Gamma(\alpha+\beta)。在机器 KPIs 中,CPU 等指标可以用 Beta 分布;

泊松分布的概率密度函数是 f(x;\lambda)=\lambda^{x}e^{-\lambda}/x!,在机器 KPIs 中,SYS_OOM 用于衡量超出内存的频率,可以用泊松分布来做。

高斯分布的概率密度函数 f(x;\mu,\sigma) = e^{-(x-\mu)^{2}/2\sigma^{2}}/(\sqrt{2\pi}\sigma)

根据论文中的陈述,机器 KPIs 分别适用于以下概率分布:

FluxRank_table_1
机器指标遵循的概率分布

通过以上公式,可以计算出每一个机器的每一个指标的 ou 两个值。

Digest Distillation

再来看一下 Digest Distillation 部分,在此部分需要对机器的 KPIs 进行聚类操作;那么就需要构造特征向量和距离函数,再加上聚类算法即可获得结果。

每一个机器的特征向量是由之前计算的 Change Degree 形成的,由于每台机器的 KPIs 都是一样的,因此可以对它们的 KPIs 的 change degree 进行排列。假设每台机器有 k 个 KPIs,那么这台机器所对应的向量就是 (o_{0},u_{0},\cdots,o_{k},u_{k})

在描述向量的相似性方面,可以使用相关性的系数,包括 Pearson 系数,Kendall tau 系数,Spearman 系数。对于两条时间序列而言,X=[x_{1},\cdots,x_{n}]Y=[y_{1},\cdots,y_{n}]

Pearson 系数指的是:\rho_{X,Y}=\sum_{i=1}^{n}(x_{i}-\overline{x})\cdot(y_{i}-\overline{y})/\sqrt{\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}\cdot\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}}, 其中 \overline{x}=\sum_{i=1}^{n}x_{i}/n\overline{y}=\sum_{i=1}^{n}y_{i}/n

Kendall tau 系数指的是:如果 (x_{i}>x_{j}y_{i}>y_{j}) 或者 (x_{i}<x_{j}y_{i}<y_{j}),那么称之为 concordant;如果 (x_{i}<x_{j}y_{i}>y_{j}) 或者 (x_{i}>x_{j}y_{i}<y_{j}),称之为 discordant;如果 x_{i}=x_{j} 或者 y_{i}=y_{j},则既不是 concordant,也不是 discordant。那么 Kendall tau 定义为 [\text{(number of concordant pairs)}-\text{(number of disordant paris)}] / [n(n-1)/2]

Spearman 系数指的是:通过原始序列变成秩次变量(rank)(从大到小降序排列即可),x_{i} 将会对应到 x_{i}',后者表示 x_{i} 在从大到小排序之后的序列 \{x_{i}\}_{1\leq i\leq n} 的位置,称之为秩次(rank),得到序列 X'=[x_{1}',\cdots,x_{n}']。对原始序列 Y=[y_{1},\cdots,y_{n}] 作同样的操作,得到 Y'=[y_{1}',\cdots,y_{n}']。一个相同的值在一列数据中必须有相同的秩次,那么在计算中采用的秩次就是数值在按从大到小排列时所在位置的平均值。如果没有相同的 rank,那么使用公式 r_{s} = 1-6\sum_{i=1}^{n}d_{i}^{2}/(n(n^{2}-1)) 进行计算,其中 d_{i}=x_{i}'-y_{i}';如果存在相同的秩次,则对 X'=[x_{1}',\cdots,x_{n}']Y'=[y_{1}',\cdots,y_{n}'] 来做 Pearson 系数即可,也就是 \rho_{X',Y'}

FluxRank_table_3
相似性函数的对比

通过作者们的实验,说明 Pearson 系数在这个数据集上效果最佳。在聚类算法的场景下,作者们同样对比了 KMeans,Gaussian Mixture,Hierarchical Clustering,DBSCAN 算法的效果,最后使用了 DBSCAN 的聚类算法。每一个聚类的结果,作者称之为一个 digest,也就是下图的 M1,M2 等聚类结果。

FluxRank_fig_6
聚类结果

Digest Ranking

最后,就是对聚类结果的排序工作。通过观察会发现:

  1. 变化开始时间(change start time)T_{c} 会在失败发生时间 T_{f} 之前;
  2. 不同的故障机器 KPIs 的 change start time 是非常接近的;
  3. 故障机器的一些 KPIs 的 change degree 是非常大的;
  4. 故障机器的占比是与故障原因相关的,故障机器越多说明故障越大;

在同一个模块下,如果出现故障机器的占比较大,那么故障将集中于这个模块下,可以通过 ratio 这个指标进行排序工作。

实验数据

在 FluxRank 论文中,作者们收集了 70 个真实的案例,然后根据实验效果获得了结果。

FluxRank_table_2
部分真实案例

在标记的时候,除了标记异常机器(Root Cause Machines,简称为 RCM)之外,也需要标记相关的指标(Relevant KPI,简称为 RK)。Root Cause Digest(简称为 RCD)把包括两个部分,不仅包括 RCM 的一个聚类结果,还包括聚类结果中的 top-five KPIs。

通过对 FluxRank 进行实验,可以得到如下实验数据:

FluxRank_table_4
FluxRank 的实验结果

其中 Recall@K 指的是:Recall@K=\text{\# of cases whose top-k digests contain RCDs}/ \text{\# of all cases}, 或者 Recall@K=\text{\# of cases whose top-k machines contain RCMs}/\text{\# of all cases}.

参考资料

  1. FluxRank: A Widely-Deployable Framework to Automatically Localizing Root Cause Machines for Software Service Failure Mitigation,Ping Liu,Yu Chen,Xiaohui Nie,Jing Zhu,Shenglin Zhang,Kaixin Sui,Ming Zhang,Dan Pei,ISSRE 2019, Berlin, Germany, Oct 28-31, 2019。
  2. Introduction to Statistical Quality Control,6th edition,Douglas C.Montgomery。
  3. Bayesian Network:https://en.wikipedia.org/wiki/Bayesian_network

 

符号计算中的深度学习方法

符号计算

符号计算一直是计算数学的重要领域之一。在开源领域,Python 的 SymPy 就可以支持符号计算。在商业化领域,Maple,Matlab,Mathematica 都能够进行符号计算。它们不仅能够做简单的实数和复数加减乘除,还能够支持数学分析,线性代数,甚至各种各样的大学数学课程。

随着人工智能的进一步发展,深度学习不仅在图像识别,自然语言处理方向上发挥着自身的价值,还在各种各样的领域展示着自己的实用性。在 2019 年底,facebook 两位研究员在 arxiv 上挂出了一篇文章《Deep Learning for Symbolic Mathematics》,在符号计算方向上引入了深度学习的工具。

要想了解符号运算,就要先知道在计算机中,是怎么对数学公式进行表示的。较为常见的表达式形如:

  • 2 + 3 * (5 + 2)
  • 3x^{2}+\cos(2x)+1
  • \frac{\partial^{2}\psi}{\partial x^{2}} - \frac{1}{v^{2}}\frac{\partial^{2}\psi}{\partial t^{2}}

在这里,数学表达式通常都会被表示成树的结构,其中树的内部节点是由算子(operator),函数(function)组成的,叶子节点由数字,变量,函数等组成。例如:

dlsm_figure1.png
图 1

图 1 的三幅图分别对应着上面的三个数学表达式。

在 Python 的 SymPy 工具中,同样可以对数学公式进行展示。其表示方法就是用 sympy.srepr

>>> import sympy
>>> x, y = sympy.symbols("x y")
>>> expr = sympy.sin(x+y) + x**2 + 1/y - 10
>>> sympy.srepr(expr)
"Add(Pow(Symbol('x'), Integer(2)), sin(Add(Symbol('x'), Symbol('y'))), Integer(-10), Pow(Symbol('y'), Integer(-1)))"
>>> expr = sympy.sin(x*y)/2 - x**2 + 1/y
>>> sympy.srepr(expr)
"Add(Mul(Integer(-1), Pow(Symbol('x'), Integer(2))), Mul(Rational(1, 2), sin(Mul(Symbol('x'), Symbol('y')))), Pow(Symbol('y'), Integer(-1)))"

dlsm_figure2.png
图 2

SymPy 的 srepr 函数的输出用树状结构来表示就是形如图 2 这种格式。叶子节点要么是 x,y 这种变量,要么是 -1 和 2 这种整数。对于一元函数而言,例如 sin 函数,就是对应唯一的一个叶子。对于二元函数而言,例如 pow,mul,add,则是对应两个叶子节点。

论文方案

在 Deep Learning for Symbolic Mathematics 这篇论文中,作者们的大致思路是分成以下几步的:

  1. 生成数据;
  2. 训练模型;
  3. 预测结果;

第一步生成数据是为了让深度学习模型是大量的已知样本来做训练;第二步训练模型是用第一步的数据把端到端的深度学习模型进行训练;第三步预测结果是给一个函数或者一个微分方程,使用已经训练好的模型来预测结果,对预测出来的多个结果进行排序,选择最有可能的那个结果作为符号计算的值。

众所周知,深度学习的训练是依赖大量的样本数据的,那么要想用深度学习来解决符号计算的问题,就要解决样本少的问题。在这篇论文中,作者们把精力投入了三个领域,分别是:

  1. 函数的积分;
  2. 一阶常微分方程;
  3. 二阶常微分方程。

在生成数据之前,作者们对数据的范围进行了必要的限制:

  1. 数学表达式最多拥有 15 个内部节点;
  2. L = 11 表示叶子节点的值只有 11 个,分别是变量 x\{-5,-4,-3,-2,-1,1,2,3,4,5\}
  3. p_{1} = 15 表示一元计算只有 15 个,分别是 \exp, \log, \sqrt, \sin, \cos, \tan, \arcsin, \arccos, \arctan, sinh, cosh, tanh, arcsinh, arccosh, arctanh
  4. p_{2} = 4 表示二元计算只有四个,分别是 +, -, *, /;

意思就是在这个有限的范围内去生成深度学习所需要的数据集。

积分数据的生成

在微积分里面,积分指的是求导的逆运算,那么形如 (f',f) 的表达式就可以作为深度学习的积分训练数据。生成积分的话其实有多种方法:

第一种方法:前向生成(Forward Generation,简写为 FWD)。主要思路就是在以上的数据范围内随机生成各种各样的方程 f,然后使用 SymPy 或者 Mathematica 等工具来计算函数 f 的积分 F,那么 (f,F) 就可以作为一个训练集。当然,有的时候函数 f 的积分是无法计算出来的,那么这种计算表达式就需要进行放弃,就不能放入训练集。

第二种方法:反向生成(Backward Generation,简写为 BWD)。由于积分是求导的逆运算,可以在以上的数据范围内随机生成各种各样的方程 f,然后计算它们的导数 f',于是 (f',f) 就可以放入积分数据的训练集。

第三种方法:分部积分(Backward generation with integration by parts,简写为 IBP)。根据分部积分的公式,如果 F'=f, G'=g,那么 \int F(x)g(x)dx=F(x)G(x) - \int f(x)G(x)dx。对于两个随机生成的函数 F, G,可以计算出它们的导数 f, g。如果 fG 在训练集合里面,那么就把 Fg 的积分计算出来放入训练集合;反之,如果 Fg 在训练集合里面,那么就把 fG 的积分计算出来放入训练集合。如果 fGFg 都没有在训练集合,并且都无法积分出来,那么就放弃该条数据。

三种方法的比较可以参见表 1,从表 1 可以看出通过分部积分(Integration by parts)所获得的样本表达式相对于前向生成和反向生成都比较长。

dlsm_table1.png
表 1

一阶常微分方程的生成

一阶常微分方程只有该函数的一阶导数,因此在构造函数的时候,作者们用了一个技巧。在随机生成表达式的时候,叶子节点的元素只从 \{x, -5,-4,-3,-2,-1,1,2,3,4,5\} 选择一个,于是随机把其中的一个整数换成变量 c。例如:在 x\log(2/x) 中就把 2 换成 c,于是得到了一个二元函数 f(x,c) = x\log(c/x)。那么就执行以下步骤:

  1. 生成二元函数 f(x,c) = x\log(c/x)
  2. 求解 c,得到 c = xe^{f(x)/x}
  3. x 进行求导得到 0 = e^{f(x)/x}(1+f'(x)-f(x)/x) = 0
  4. 简化后得到 x+xf'(x) - f(x) =0,也就是 x+xy'-y=0

此时,一阶常微分方程的训练数据就是 (x+xy'-y=0, x\log(c/x))

二阶常微分方程不仅可能由该函数的一阶导数,还必须有二阶导数,那么此时就要求导两次:

  1. 生成三元函数 f(x,c_{1},c_{2}) = c_{1}e^{x}+c_{2}e^{-x}
  2. 求解 c_{2} 得到 c_{2} = f(x,c_{1},c_{2}) e^{x}- c_{1}e^{2x}
  3. x 进行求导得到 0 = e^{x}(\partial f(x,c_{1},c_{2})/\partial x + f(x,c_{1},c_{2})) - 2c_{1}e^{2x} = 0
  4. 求解 c_{1} 得到 c_{1} = e^{-x}(\partial f(x,c_{1},c_{2})/\partial x+f(x))/2
  5. x 进行求导得到 0 = e^{-x}(\partial^{2} f(x,c_{1},c_{2})/\partial x^{2} - f(x,c_{1},c_{2}))=0
  6. 简化后得到 \partial^{2} f(x,c_{1},c_{2})/\partial x^{2} - f(x,c_{1},c_{2})=0,也就是 y''-y=0

那么此时的二阶微分方程的训练数据就是 (y''-y=0, c_{1}e^{x}+c_{2}e^{-x})

需要注意的事情就是,在生成数据的时候,一旦无法求出 c_{1}, c_{2} 的表达式,该数据就需要放弃,重新生成新的数据。

数据处理

  1. 数学表达式的简化(expression simplification):例如 x+1+1 可以简化成 x+2\sin^{2}(x)+\cos^{2}(x) 可以简化成 1。
  2. 参数的简化(coefficient simplification):例如 \log(x^{2}) + c\log(x) 可以简化成 c\log(x)
  3. 无效表达式的过滤(invalid expression filter):例如 \sqrt{2}, \log(0) 等。

树状结构的表达式,是使用前缀表达式来写成一个句子的。例如 2+3 就写成 + 2 3,2 + x 就写成 + 2 x。

模型训练

在这里,作者们用了 Transformer 模型,8 attention heads,6 layers,512 个维度,并且发现用更复杂的模型并没有提升其效果。在预测的时候,使用不同的 Beam Size,其准确率是不一样的,在 Beam Size = 50 的时候,效果比较好。参见表 2。

dlsm_table2.png
表 2

在与其他数学软件的比较中, 作者限制了 Mathematica 的运行时间为 30s。

dlsm_table3.png

并且举出了几个现有模型能够计算出来,Mathematica 无法计算的例子。

dlsm_table4.png

在求解的等价性方面,可以根据 Score 逆序排列,然后这个例子的 Top10 都是等价的。

dlsm_table5.png

在使用 SymPy 的过程中,可以获得各种各样的积分表达式如下:

dlsm_table7.png

结论

符号计算在 1960 年代末就已经在研究了,有诸多的符号计算软件,例如 Matlab,Mathematica,Maple,PARI,SAGE 等。在这篇论文中,作者们使用标准的 seq2seq 模型来对生成的数据进行训练,然后获得积分,一阶常微分方程,二阶常微分方程的解。在传统符号计算的基础上,开拓了一种新的思路。

深度学习在时间序列分类中的应用

本篇博客将会分享几篇文章,其内容主要集中在深度学习算法在时间序列分类中的应用。

无论是图像分类,文本分类,还是推荐系统的物品分类,都是机器学习中的常见问题和应用场景。同样的,时间序列的分类问题也是研究时间序列领域的重要问题之一。近期,神经网络算法被用于物体识别,人脸识别,语音分类等方向中,于是有学者用深度学习来做时间序列的分类

假设

X=\{x_{1},\cdots,x_{n}\}

是一个长度为 n 的时间序列,高维时间序列

X = \{X^{1},\cdots,X^{M}\}

则是由 M 个不同的单维时间序列而组成的,对于每一个 1\leq i\leq M 而言,时间序列 X^{i} 的长度都是 n. 而时间序列的分类数据通常来说都是这种格式:数据集

D=\{(X_{1},Y_{1}),(X_{2},Y_{2}),\cdots,(X_{N},Y_{N})\}

表示时间序列与之相应的标签,而 Y_{i} 是 one hot 编码,长度为 K(表示有 K 个类别)。

整体来看,时间序列分类的深度学习方案大体是这个样子的:输入的是时间序列,通过某个神经网络算法进行端到端的训练,最后输出相应的分类概率。

ReviewFigure1.png

而做时间序列分类的深度学习算法分成生成式(Generative)判别式(Discriminative)两种方法。在生成式里面包括 Auto Encoder 和 Echo State Networks 等算法,在判别式里面,包括时间序列的特征工程和各种有监督算法,还有端到端的深度学习方法。在端到端的深度学习方法里面,包括前馈神经网络,卷积神经网络,或者其余混合模型等常见算法。

ReviewFigure5.png

深度学习算法在时间序列分类中的应用:Baseline

这一部分将会介绍用神经网络算法来做时间序列分类的 Baseline,其中包括三种算法,分别是多层感知机(MLP),FCN(Fully Convolutional Network)和 ResNet。其论文的全名是《Time Series Classification from Scratch with Deep Neural Networks: A Strong Baseline》。这篇论文中使用的神经网络框架如下图所示:

DNN_Baseline_结构1.png

多层感知机(MLP)模型使用了全连接层,每个隐藏层大约 500 个神经元,然后使用 ReLU 作为激活函数,同时使用 Dropout 来防止过拟合,最后一层是 Softmax 层。MLP 中一个基础的块包括:

\tilde{x} = f_{dropout, p}(x),

y = W\cdot \tilde{x} + b,

h = ReLU(y).

除了前馈神经网络之外,全卷积网络(FCN)同样可以作为时间序列的特征提取工具,一个卷积块包括:

y = W \otimes x + b,

s = BN(y),

h = ReLU(s),

在这里,\otimes 指的是卷积算子,BN 指的是 Batch Normalization,ReLU 则是激活函数。

Residual Network 是在 FCN 的基础上进行的改造。令 Block_{k} 来表示第 k 个卷积块,而 Residual 块就定义为:

h_{1} = Block_{k_{1}}(x),

h_{2} = Block_{k_{2}}(h_{1}),

h_{3} = Block_{k_{3}}(h_{2}),

y=h_{3} + x,

\hat{h}=ReLU(y).

其中,k_{1} = 64, k_{2} = 128, k_{3} = 128

ReviewFigure6_ResNet.png

评价指标

Mean Per Class Error (in Multi-class Classification only) is the average of the errors of each class in your multi-class data set. This metric speaks toward misclassification of the data across the classes. The lower this metric, the better.

模型的评价指标使用的是 Mean Per-Class Error,指的是在多分类场景下,每一类(Class)错误率的平均值。换句话说,一个数据集 D=\{d_{k}\}_{1\leq k\leq K} 是由 K 个类的元素构成的,每个类的标签是 C=\{c_{k}\}_{1\leq k\leq K},通过模型其实可以计算出模型对每一个类的错误率 e_{k},那么模型的 MPCE 就是:MPCE= \sum_{1\leq k\leq K} e_{k}/K.

其实验结论是:

DNN_Baseline_实验数据1.png

MSCNN

MSCNN 的全称是 Multi-Scale Convolutional Neural Networks,相应的论文是《Multi-Scale Convolutional Neural Networks for Time Series Classification》

在时间序列的分类算法里面,通常来说,可以分成以下几种:

  1. 基于距离的方法(distance-based methods):kNN,SVM(相似核),DTW;
  2. 基于特征的方法(feature-based methods):SVM,逻辑回归等;
  3. 基于神经网络的方法(neural network-based methods):CNN 等;

正如前文所提到的,一条时间序列通常可以写作 T=\{t_{1},\cdots,t_{n}\},其中 t_{i} 表示在时间戳 i 下的取值,并且时间序列 T 的长度是 n。在时间序列分类的场景下,每一条时间序列对应着唯一的一个标签(label),也就是说 D=\{(T_{i},y_{i})\}_{i=1}^{N}。其中 D 集合里面包含 N 条时间序列,每条时间序列 T_{i} 对应着一个标签 y_{i}y_{i} 表示分类值集合 \mathcal{C} = \{1,\cdots,C\} 中的元素,C\in \mathbb{Z}^{+}

MSCNN 的整体结构:

在 Multi-Scale Convolutional Neural Network(MSCNN)中,包括几个串行的阶段,

  1. 变换阶段(Transformation Stage):包括恒等变换,下采样,谱变换等变换方式,每一种方式都是一个分支,并且也是卷积神经网络的输入;
  2. 局部卷积(Local Convolution Stage):使用卷积层来对不同的输入提取特征,不同的输入分支之间是相互独立的,输出的时候都会经过一个最大值池化(max pooling)的过程;
  3. 整体卷积(Full Convolution Stage):把上一步提取到的特征进行拼接(concatenate),然后使用全连接层并且加上一个 softmax 层来做多分类。

如下图所示,MSCNN 是一个端到端的训练网络结构,所有参数都是通过后向传播算法得到的。

MSCNN结构1.png

首先来看神经网络的第一步,变换阶段(Transformation Stage),也就是神经网络的多尺度的输入。在不同的尺度下,神经网络能够提炼到不同类型的特征。长期的特征(long-term features)反映了时间序列的整体趋势,短期的特征(short-term features)反映了时间序列的局部的微妙变化。要想判断时间序列的形状,不仅要参考整体的特征,也要参考局部的特征,这两者对于判断时间序列的形状都具有一定的辅助作用。

在 Transformation Stage,identity map 指的是恒等变换,也就是说时间序列是原封不动的作为神经网络的输入数据。对于 Smoothing Transformation,指的就是对时间序列进行必要的平滑操作,将新的时间序列作为神经网络的输入数据。在这种情况下,我们可以对时间序列 T=\{t_{1},\cdots,t_{n}\} 进行移动平滑,i.e.

T^{\ell}=(x_{i}+x_{i+1}+\cdots+x_{i+\ell-1})/\ell, 0\leq i\leq n-\ell+1,

其中的 \ell\in \mathbb{Z}^{+} 表示窗口长度。对于不同的窗口长度 \ell,我们可以的到不同的时间序列平滑序列,但是它们的长度都是一样的,都是原始的时间序列长度 n

而下采样(down sampling)指的则是对时间序列的间隔进行抽样操作。假设时间序列 T=\{t_{1},\cdots,t_{n}\},下采样的比例是 k,也就是说我们每隔 k 个点保留时间序列的取值,i.e.

T^{k} = \{t_{1+k\cdot i}\}, 0\leq i \leq [(n-1)/k].

用这种方法,我们可以对 k=1,2,3,\cdots 来进行下采样的时间序列提取。在进行了恒等变换,平滑变换,下采样之后,时间序列就可以变成多种形式,作为神经网络的输入。

其次,在神经网络部分,本文使用了一维(1-D)的卷积层和最大值池化的方法来提取特征,并且在局部卷积阶段之后把提炼到的抽象特征进行拼接(concatenate)。拼接完了之后,持续使用卷积层和池化层进行特征的提取,然后使用全连接层(fully connected layers)和 softmax 层来进行时间序列类别的预测。

数据增强

在深度学习里面,由于是端到端的训练网络,因此是需要相对多的样本数据的,于是有的时候需要进行数据增强(data augmentation)。也就是在现有的基础上获得更多的训练数据。对于时间序列 T=\{t_{1},\cdots,t_{n}\},可以定义一个子序列:

S_{i:j} = \{t_{i}, t_{i+1},\cdots,t_{j}\}, 1\leq i,\leq j\leq n

对于正整数 s\in \mathbb{Z}^{+},可以生成 n-s+1 个子序列如下所示:

Slicing(T,s) = \{S_{1:s},S_{2:s+1},\cdots,S_{n-s+1:n}\},

这些子序列的标签与原始的时间序列 T 是一样的。

本文用到的数据集情况如下表所示:

MSCNN数据集1.png

实验数据如下图所示:

MSCNN实验1MSCNN数据集2

结论

本文使用了 MCNN 来对变换之后的时间序列进行特征提取,并且进行了端到端的模型训练。并且也讨论了卷积神经网络使用在 shapelet learning 上的一些逻辑和方法,然后解释了 MSCNN 在时间序列分类上能够有不错表现的原因。但是所有的 TSC 数据集都不算特别大,对端到端的训练模式有一定的限制。

GASF 和 GADF 方法

这篇文章《Imaging Time Series to Improve Classification and Imputation》介绍了如何把时间序列转换成图像,包括 GASF 方法和 GADF 方法。

假设时间序列是 X = \{x_{1},\cdots, x_{n}\},长度是 n,我们可以使用归一化方法把时间序列压缩到 [0,1] 或者 [-1,1]

\tilde{x}_{0}^{i} = (x_{i} - \min(X))/(\max(X) -\min(X)),

\tilde{x}_{-1}^{i} = ((x_{i} - \max(X)) + (x_{i}-\min(X)))/(\max(X) - \min(X)),

此时的 \tilde{x}_{0}^{i}\in[0,1], \forall 1\leq i\leq n\tilde{x}_{-1}^{i} \in [-1,1],\forall 1\leq i\leq n。于是可以使用三角函数来代替归一化之后的值。下面通用 \tilde{x}_{i} 来表示归一化之后的时间序列,令 \phi_{i} = \arccos(\tilde{x}_{i})\tilde{x}_{i} \in [-1,1]1\leq i\leq n。因此,\phi_{i}\in[0,\pi],于是,\sin(\phi_{i}) \geq 0

定义矩阵 GASF(Gramian Angular Summation Field)

GASF = [\cos(\phi_{i}+\phi_{j})]_{1\leq i,j\leq n}

于是,

GASF = [\cos(\phi_{i})\cdot \cos(\phi_{j}) - \sin(\phi_{i})\cdot \sin(\phi_{j})]_{n\times n}

\tilde{X}=(\cos(\phi_{1}),\cdots,\cos(\phi_{n}))^{T},可以得到

GASF = \tilde{X} \cdot \tilde{X}^{T} - \sqrt{I - \tilde{X}^{2}} \cdot \sqrt{I - \tilde{X^{T}}^{2}}

以上的都是 element 乘法和加法,I 表示单位矩阵。它的对角矩阵是

diag(GASF) = \{\cos(2\phi_{1}),\cdots, \cos(2\phi_{n})\}

= \{2\cos^{2}(\phi_{1})-1,\cdots,2\cos^{2}(\phi_{n})-1)\} = \{GASF_{ii}\}_{1\leq i\leq n}.

如果是使用 min-max normalization 的话,是可以从 diag(GASF) 反推出 \tilde{x}_{i} 的。因为,2 \tilde{x}_{i}^{2} - 1 = 2\cos^{2}(\phi_{i}) - 1 = GASF_{ii},可以得到 y_{i} = \sqrt{(GASF_{ii}+1)/2}

定义 GADF(Gramian Angular Difference Field)如下:

GADF = [\sin(\phi_{i}-\phi_{j})]_{1\leq i,j\leq n}

= [\sin(\phi_{i})\cdot cos(\phi_{j}) - \cos(\phi_{i})\cdot\sin(\phi_{j})]_{1\leq i,j\leq n}

= \sqrt{1-X^{2}}\cdot X^{T} - X \cdot \sqrt{1-(X^{T})^{2}}.

GASFGADF1.png

Markov Transition Field(MTF)

除了 GSAF 和 GSDF 之外,《Imaging Time Series to Improve Classification and Imputation》,《Encoding Time Series as Images for Visual Inspection and Classification Using Tiled Convolutional Neural Networks》,《Encoding Temporal Markov Dynamics in Graph for Time Series Visualization》也提到了把时间序列转换成矩阵 Image 的算法 MTF。在 pyts 开源工具库里面,也提到了 MTF 算法的源码。

假设时间序列是 X = \{x_{1},\cdots,x_{n}\},我们把它们的值域分成 Q 个桶,那么每一个 x_{i} 都可以被映射到一个相应的 q_{j} 上。于是我们可以建立一个 Q\times Q 的矩阵 Ww_{ij} 表示在桶 j 中的元素被在桶 i 中的元素跟随的概率,也就是说 w_{ij} = P(x_{t}\in q_{i}|x_{t-1}\in q_{j}),同时,它也满足 \sum_{j=1}^{Q}w_{ij} =1。于是,得到矩阵 W = (w_{ij})_{1\leq i,j\leq Q}

MTF1.png

除此之外,我们也能够计算一个迁移概率矩阵 M。其中 $m_{ij}$ 表示桶 i 中的元素迁移至桶 j 中的概率 P(q_{i}\rightarrow q_{j}),同样有 \sum_{1\leq j\leq Q} m_{ij} =1。因此,我们同样可以构造出一个 Q\times Q 的矩阵将时间序列可视化。

时间序列的降维方法有两种:

  1. 分段聚合(PAA):使用局部平均等方法,把时间序列进行降维;
  2. 核变换(Kernel):使用 Bivariate Gaussian 核或者均值核来把时间序列进行降维。

在把时间序列进行可视化之后,对于时间序列分类的场景,就可以使用 CNN 的技术方案来做了。如下图所示:

tildeCNN.png

其实验数据效果如下:

tildeCNN实验数据1.png

Time Le-Net

在本篇文章《Data Augmentation for Time Series Classification using Convolutional Neural Networks》中,主要用到了卷积神经网络来做时间序列的分类。

time_lenet_1.png

除此之外,也使用了不少数据增强(Data Augmentation)的技术。包括前面提到的 Window Slicing(WS)方法。也考虑了 Warping 的变换技巧,例如 Warping Ratio = 1/2 或者 2。这种时间扭曲指标比率可以通过交叉验证来选择。该方法叫做 Window Warping(WW)技术。

另外也有其余论文使用卷积神经网络做时间序列分类,例如《Convolutional neural networks for time series classification》,如下图所示:

CNN_figure1.png

Multi-Channels Deep Convolutional Neural Networks

在高维时间序列的分类中,有人提出用多通道的卷积神经网络来进行建模。

multichannelcnn_figure3.png

整体来看,分成四个部分。前三个部分作为特征提取工具,最后一层作为分类工具。

  1. Filter Layer:
  2. Activation Layer:
  3. Pooling Layer:
  4. Fully-Connected Layer:

实验对比数据如下:

multichannelcnn_table1.png

结论:

在本篇博客中,列举了一些深度学习算法在时间序列分类中的应用,也介绍了部分数据增强的方法和时间序列数据变换的方法。从以上各篇文章的介绍来看,深度学习在时间序列分类领域上应该是大有可为的。

参考资料:

  1. Wang Z, Yan W, Oates T. Time series classification from scratch with deep neural networks: A strong baseline[C]//2017 international joint conference on neural networks (IJCNN). IEEE, 2017: 1578-1585.
  2. Cui Z, Chen W, Chen Y. Multi-scale convolutional neural networks for time series classification[J]. arXiv preprint arXiv:1603.06995, 2016.
  3. Wang Z, Oates T. Imaging time-series to improve classification and imputation[C]//Twenty-Fourth International Joint Conference on Artificial Intelligence. 2015.
  4. Wang Z, Oates T. Encoding time series as images for visual inspection and classification using tiled convolutional neural networks[C]//Workshops at the Twenty-Ninth AAAI Conference on Artificial Intelligence. 2015.
  5. Liu L. Encoding Temporal Markov Dynamics in Graph for Time Series Visualization, Arxiv, 2016.
  6. Fawaz H I, Forestier G, Weber J, et al. Deep learning for time series classification: a review[J]. Data Mining and Knowledge Discovery, 2019, 33(4): 917-963.
  7. Zhao B, Lu H, Chen S, et al. Convolutional neural networks for time series classification[J]. Journal of Systems Engineering and Electronics, 2017, 28(1): 162-169.
  8. Le Guennec A, Malinowski S, Tavenard R. Data augmentation for time series classification using convolutional neural networks[C]. 2016.

时间序列的标签

本篇文章是为了介绍一种基于少量样本标记而获得更多样本的方法,论文的原文是《Label-Less: A Semi-Automatic Labeling Tool for KPI Anomalies》,是清华大学与多家公司(必示科技,中国建设银行等)的合作论文。

在时间序列异常检测中,因为标注的成本比较大,于是需要寻找一种较少而高效地标注时间序列异常点的方法。在该论文中,Alibaba,Tencent,Baidu,eBay,Sogou提供了上千条时间序列(每条时间序列大约是2-6个月的时间跨度),作者们进行了 30 条 KPIs 的标注工作。但是其标注成本依旧是很大的,于是作者们想到了一种异常相似搜索(anomaly similarity search)的算法,目标是对已经标注好的时间序列异常模式进行模版搜索。目的就是达到 label-less,也就是较少的标注而获得更多的标注数据。

在本篇论文中,在异常检测的过程中,作者们使用了时间序列的预测模型(time series prediction models)来获得时间序列的特征,使用了孤立森林(Isolation Forest) 来对时间序列的特征来做无监督的异常检测。并且其效果由于 one class svm 算法和 local outlier factor 算法。在搜索的部分,作者使用了加速版的 DTW 算法(accelerated dynamic time warping approach)来做相似度的搜索和模式的匹配。其中也尝试了各种技巧和方法,包括 constrained DTW,LB Keogh 方法,early stopping 算法等工具。

整个 Label-Less 的架构图如下表示:

Fig2.png

其中的 Operators 指的是业务运维人员,面对着无标记的多条时间序列曲线。系统首先会进行无监督的异常检测算法啊,包括时间序列的预处理(归一化等)操作,然后使用差分(Difference),移动平均算法(moving average),带权重的移动平均算法(weighted moving average),指数移动平均(ewma),holt winters,ARIMA 等算法来做特征的提取。此时,对于不同的时间序列预测工具,我们可以得到不同的预测值,然后把预测值减去实际值并且取绝对值,就得到时间序列的误差序列。i.e. |p_{i} - x_{i}| 就作为数据点 x_{i} 的特征。

在这种情况下,由于用了六个时间序列预测算法,因此原始的时间序列 X (n\times 1) 就可以变成特征矩阵 X' (n\times 1)。对于特征矩阵 X' 可以使用 isolation forest 来做无监督的异常检测并且做阈值的设定;如下图所示:

Fig3.png

而另外的一部分的异常相似搜索(anomaly similarity search)是在第一部分的基础上在做的,Unsupervised Anomaly Detection 会输出疑似异常或者候选异常,并且基于已知的异常模板(Anomaly Template)进行相似度的匹配,此时可以使用 accelerated DTW 算法,选择出最相似的 Top-K 异常,然后运维人员进行标注,得到更多的样本。

由于,对于两条长度分别是 mn 的时间序列,DTW 相似度算法的时间复杂度是 O(mn),因此在搜索的时候需要必要的加速工作。在这种地方,作者们使用了 LB-Kim,LB-Keogh,LB-Keogh-Reverse 算法来做搜索的加速工作。而这些的时间复杂度是 O(m+n)。整体的思路是,如果两条时间序列 qc 的 LB-Kim,LB-Keogh,LB-Keogh-Reverse 的下界大于某个阈值,则不计算它们之间的 DTW 距离。否则就开始计算 DTW。并且在计算 DTW 的时候,如果大于下界,则会提前终止(early stopping),不会继续计算下去。如果都没有大于阈值,则把这个候选曲线和 dist 距离放入列表,最后根据列表中的 dist 来做距离的逆序排列。

整体流程如下:

AnomalySimilaritySearch

其运行速度也比直接使用 DTW 快不少:

Table3.png

Label-Less 的交互页面如下所示:

图(a)表示使用无监督算法获得的疑似异常;

图(b)表示使用异常搜索算法获得的异常结果。

Fig10.png

下图则表示模板,m 表示模板的长度,c 表示相似的异常候选集个数;

Fig11

总结:

整体来看,本文提供了一种通过少量人工标注无监督算法相似度算法来获得更多样本的方法。在候选的时间序列条数足够多的时候,是可以进行时间序列的相似度匹配的。这给未来在运维领域提供海量的时间序列标注数据给予了一定的技术支持。

 

时间序列的联动分析

背景介绍

在互联网公司里面,通常都会监控成千上万的时间序列,用于保障整个系统或者平台的稳定性。在这种情况下,如果能够对多条时间序列之间判断其是否相关,则对于监控而言是非常有效的。基于以上的实际情况,清华大学与 Alibaba 集团在2019年一起合作了论文《CoFlux: Robustly Correlating KPIs by Fluctuations for Service Troubleshooting》,并且发表在 IWQos 2019 上。CoFlux 这个方法可以对多条时间序列来做分析,并且主要用途包括以下几点:

  1. 告警压缩和收敛;
  2. 推荐与已知告警相关的 Top N 的告警;
  3. 在已有的业务范围内(例如数据库的实例)构建异常波动传播链;

kpis.png

CoFlux 的整体介绍

从论文的介绍中来看,CoFlux 的输入和输出分别是:

输入:两条时间序列

输出:这两条时间序列的以下信息

  1. 波动相关性:两条时间序列是否存在波动相关性?
  2. 前后顺序:如果两条时间序列相关,那么它们的前后波动顺序是什么?是同时发生异常还是存在固定的前后顺序?
  3. 方向性:如果两条时间序列是波动相关的,那么它们的波动方向是什么?是一致还是相反?

Remark. CoFlux 的关键点就是并没有对时间序列做异常检测算法,而是直接从时间序列的历史数据(历史半个月或者一个月)出发,判断两条时间序列之间的波动相关性,并且进一步的分析先后顺序与波动方向。

从论文的介绍中来看,CoFlux 的流程图如下图所示:

coflux流程图1

如果两条时间序列 XY 存在波动相关性,则需要输出这两条时间序列的波动先后顺序和是否同向波动。如果两条时间序列 XY 并不存在波动相关性的话,则不需要判断波动先后顺序和是否同向波动。

coflux流程图2

CoFlux 的细节阐述

已知一个长度是 n 的时间序列 S=\{s_{1},\cdots,s_{n}\},对于任意一个 detector,可以得到一条关于 S 的预测值曲线 P=\{p_{1},\cdots,p_{n}\}。于是针对某个 detector 可以得到一个波动特征序列 E=\{\epsilon_{1},\cdots,\epsilon_{n}\},其中 \epsilon_{i} = s_{i} - p_{i}1\leq i\leq n。因此,一个detector 可以对应一个波动序列特征,也是一个时间序列。因此,对于 m 个 detector,可以对应 m 条波动特征序列,并且它们的长度都是 n

在 CoFlux 算法的内部,根据不同的参数使用了总共 86 个 detector,大致列举如下:

  • Difference:根据昨天,七天前的数据来做差分;
  • Holt-Winters:\{\alpha,\beta,\gamma\} \in \{0.2,0.4,0.6,0.8\}
  • 历史上的均值 & 历史上的中位数:1,2,3,4 周;
  • TSD & TSD 中位数:1,2,3,4 周;
  • Wavelet:1,3,5,7 天;
  • 移动平均算法:MA,WMA,EWMA。PS:根据作者们的说法,在这里,MA等方法并不适用。

detectors

根据直觉来看,

  • 对于任何一条时间序列 kpi,总有一个 detector 可以相对准确地提炼到其波动特征;
  • 如果两条时间序列 XY 波动相关,那么 X 的一个波动特征序列与 Y 的一个波动特征序列应该也是相关的;

Remark. 两条时间序列的波动特征可以对齐同一个 detector,也可以不做对齐工作。如果是前者的话,时间复杂度低;后者的话,时间复杂度高。

下图是从时间序列中提取波动特征曲线的案例:

fluxfeatures.png

提炼时间序列的波动曲线特征只是第一步,后续 CoFlux 还有几个关键的步骤:

  • 特征工程的扩大(amplify): 对波动序列特征进行放大,让某些波动序列特征更加明显;
  • Correlation Measurement:用于解决时间序列存在时间前后的漂移,两条时间序列之间存在 lag 的情况,因此需要对其中一条时间序列做平移操作;
  • CoFlux 考虑了历史数据(历史半个月或者一个月)作为参考,并且一个范围内的 kpi 数量不超过 60 条;

下面来一一讲解这些技术方案,对于每一条波动特征曲线(Flux-Features),按照以下几个步骤来进行操作:

Step 1:对波动特征曲线 E=\{\epsilon_{1},\cdots,\epsilon_{n}\} 做 z-score 的归一化,i.e.

\mu = \frac{\sum_{i=1}^{n}\epsilon_{i}}{n},
\delta = \sqrt{\frac{\sum_{i=1}^{n}(\epsilon_{i}-\mu)^{2}}{n}}.

Step 2:对归一化之后的波动特征曲线做特征放大(feature amplification):定义函数 f_{\alpha,\beta}(x) 如下:

f_{\alpha,\beta}(x)= \begin{cases} e^{\alpha\min(x,\beta)} - 1, \text{ when } x\geq 0,\\ -e^{\alpha\min(|x|,\beta)} + 1, \text{ when } x< 0. \end{cases}

E=\{\epsilon_{1},\cdots,\epsilon_{n}\} 放大之后的波动特征曲线(amplified flux feature)就是:\hat{E}=\{f(\epsilon_{1}),\cdots,f(\epsilon_{n})\}.

Step 3:对于两条放大之后的波动特征曲线(amplified flux features)G=\{g_{1},\cdots,g_{\ell}\}H=\{h_{1},\cdots,h_{\ell}\},可以计算它们之间的相关性,先后顺序,是否同向。

G_{s}= \begin{cases} \{0,\cdots,0,g_{1},\cdots, g_{\ell-s}\}, \text{ when } s\geq 0, \\ \{g_{1-s},\cdots,g_{\ell},0,\cdots,0\}, \text{ when } s< 0. \end{cases}

这里的 0 的个数是 |s| 个。其中,-\ell<s<\ell。特别地,当 s=0 时,G_{0}=\{g_{1},\cdots,g_{s}\}=G,那么我们可以定义 G_{s}H 的内积是:R(G_{s},H) = G_{s}\cdot H,

这里的 \cdot 指的是向量之间的内积(inner product)。同时可以定义相关性(Cross Correlation)为:CC(G_{s},H) = \frac{R(G_{s},H)}{\sqrt{R(G_{s},G_{s})\cdot R(H,H)}}.

由于波动有可能是反向的,那么在这里我们不仅要考虑相关性是大于零的情况,也需要考虑小于零的情况。于是,

minCC = \min_{-\ell<s<\ell}CC(G_{s},H),
maxCC = \max_{-\ell<s<\ell}CC(G_{s},H).

则最小值或者最大值的指标分别是

s_{1}=argmin_{-\ell<s<\ell}CC(G_{s},H),
s_{2}=argmax_{-\ell<s<\ell}CC(G_{s},H).


FCC(G,H) = \begin{cases} (minCC, s_{1}), \text{ when } |maxCC|<|minCC|, \\ (maxCC, s_{2}), \text{ when } |maxCC|\geq|minCC|. \end{cases}

从定义中可以看出,FCC(G,H) 是一个元组,里面蕴含着三个信息,分别是相关性,波动方向,前后顺序。FCC(G,H) \in [-1,1],越接近 1 或者 -1 就表示放大之后的波动特征曲线 GH 越相关。正值的 FCC(G,H) 表示 GH 的波动方向相同,是正相关;负值的 FCC(G,H) 表示 GH 的波动方向想法,是负相关。通过对 s<0 或者 s\geq 0 的分析就可以判断先后顺序。因此,CoFlux 方法的是通过对 FCC(G,H) 的分析来得到最终结果的。

在最后的相关性分析里面,其实伪代码正如论文中所示。先考虑是否存在相关性,再考虑基于相关性下的先后顺序和波动方向。

correlationmeasurement

 

CoFlux 的实战效果

从论文中看,CoFlux 的数据集基本上是小于 60 条时间序列曲线。其中包括 CPU,错误率,错误数,内存使用率,成功率等不同的指标。

datasets

datasets2.png

从运行时间上来看,对于一周的时间序列集合(< 60条)而言,CoFlux 基本上能够在 30 分钟内计算完毕,得到最终的运算结果。

executiontime.png

其效果的评价指标基本上就是机器学习中的常见评价指标了,准确率,召回率之类的。

评价指标

从 F1-Score 的评价指标来看,CoFlux 的效果优于其他算法。

experiments.png

告警压缩

如果对时间序列之间进行告警压缩的话,其实可以大量减少运维人员的工作量。在 CoFlux 里面,时间序列曲线被分成了三类,也就是三个颜色最深的模块。因此 21 条时间序列的告警量在实际中有可能只有三条告警。

alarmclustering

告警关联

在实际运维场景中,除了对告警进行压缩之外,也需要对告警进行关联性的分析。例如一条告警发生了,运维人员都希望知道与它相关的其他告警是什么,这样可以方便运维人员定位问题。

alarmcorrelation

构建告警关系链

在一些相对封闭的场景下,例如 mysql 数据库,通过对它里面的时间序列进行分析。不仅可以得到告警之间是否存在相关性,还可以对先后顺序,波动顺序进行分析。

mysql

 

结论

时间序列之间的联动分析是在运维领域场景下的常见技术,不仅可以做告警的压缩,也能够做告警的关联,还能够构建告警的关系链。在未来的工作中,作者们提到将会用深度学习的方法来进行关联和告警的分析,从而进一步加深对时间序列的研究。

计算机视觉中的注意力机制

引言

在机器翻译(Machine Translation)或者自然语言处理(Natural Language Processing)领域,以前都是使用数理统计的方法来进行分析和处理。近些年来,随着 AlphaGo 的兴起,除了在游戏AI领域,深度学习在计算机视觉领域,机器翻译和自然语言处理领域也有着巨大的用武之地。在 2016 年,随着深度学习的进一步发展,seq2seq 的训练模式和翻译模式已经开始进入人们的视野。除此之外,在端到端的训练方法中,除了需要海量的业务数据之外,在网络结构中加入一些重要的模块也是非常必要的。在此情形下,基于循环神经网咯(Recurrent Neural Network)的注意力机制(Attention Mechanism)进入了人们的视野。除了之前提到的机器翻译和自然语言处理领域之外,计算机视觉中的注意力机制也是十分有趣的,本文将会简要介绍一下计算机视觉领域中的注意力方法。在此事先声明一下,笔者并不是从事这几个领域的,可能在撰写文章的过程中会有些理解不到位的地方,请各位读者指出其中的不足。

LSTM_1

注意力机制

顾名思义,注意力机制是本质上是为了模仿人类观察物品的方式。通常来说,人们在看一张图片的时候,除了从整体把握一幅图片之外,也会更加关注图片的某个局部信息,例如局部桌子的位置,商品的种类等等。在翻译领域,每次人们翻译一段话的时候,通常都是从句子入手,但是在阅读整个句子的时候,肯定就需要关注词语本身的信息,以及词语前后关系的信息和上下文的信息。在自然语言处理方向,如果要进行情感分类的话,在某个句子里面,肯定会涉及到表达情感的词语,包括但不限于“高兴”,“沮丧”,“开心”等关键词。而这些句子里面的其他词语,则是上下文的关系,并不是它们没有用,而是它们所起的作用没有那些表达情感的关键词大。

在以上描述下,注意力机制其实包含两个部分

  1. 注意力机制需要决定整段输入的哪个部分需要更加关注;
  2. 从关键的部分进行特征提取,得到重要的信息。

通常来说,在机器翻译或者自然语言处理领域,人们阅读和理解一句话或者一段话其实是有着一定的先后顺序的,并且按照语言学的语法规则来进行阅读理解。在图片分类领域,人们看一幅图也是按照先整体再局部,或者先局部再整体来看的。再看局部的时候,尤其是手写的手机号,门牌号等信息,都是有先后顺序的。为了模拟人脑的思维方式和理解模式,循环神经网络(RNN)在处理这种具有明显先后顺序的问题上有着独特的优势,因此,Attention 机制通常都会应用在循环神经网络上面。

虽然,按照上面的描述,机器翻译,自然语言处理,计算机视觉领域的注意力机制差不多,但是其实仔细推敲起来,这三者的注意力机制是有明显区别的。

  1. 在机器翻译领域,翻译人员需要把已有的一句话翻译成另外一种语言的一句话。例如把一句话从英文翻译到中文,把中文翻译到法语。在这种情况下,输入语言和输出语言的词语之间的先后顺序其实是相对固定的,是具有一定的语法规则的;
  2. 在视频分类或者情感识别领域,视频的先后顺序是由时间戳和相应的片段组成的,输入的就是一段视频里面的关键片段,也就是一系列具有先后顺序的图片的组合。NLP 中的情感识别问题也是一样的,语言本身就具有先后顺序的特点;
  3. 图像识别,物体检测领域与前面两个有本质的不同。因为物体检测其实是在一幅图里面挖掘出必要的物体结构或者位置信息,在这种情况下,它的输入就是一幅图片,并没有非常明显的先后顺序,而且从人脑的角度来看,由于个体的差异性,很难找到一个通用的观察图片的方法。由于每个人都有着自己观察的先后顺序,因此很难统一成一个整体。

在这种情况下,机器翻译和自然语言处理领域使用基于 RNN 的 Attention 机制就变得相对自然,而计算机视觉领域领域则需要必要的改造才能够使用 Attention 机制。

LSTM_3

基于 RNN 的注意力机制

通常来说,RNN 等深度神经网络可以进行端到端的训练和预测,在机器翻译领域和或者文本识别领域有着独特的优势。对于端到端的 RNN 来说,有一个更简洁的名字叫做 sequence to sequence,简写就是 seq2seq。顾名思义,输入层是一句话,输出层是另外一句话,中间层包括编码和解码两个步骤。

而基于 RNN 的注意力机制指的是,对于 seq2seq 的诸多问题,在输入层和输出层之间,也就是词语(Items)与词语之间,存在着某种隐含的联系。例如:“中国” -> “China”,“Excellent” -> “优秀的”。在这种情况下,每次进行机器翻译的时候,模型需要了解当前更加关注某个词语或者某几个词语,只有这样才能够在整句话中进行必要的提炼。在这些初步的思考下,基于 RNN 的 Attention 机制就是:

  1. 建立一个编码(Encoder)和解码(Decoder)的非线性模型,神经网络的参数足够多,能够存储足够的信息;
  2. 除了关注句子的整体信息之外,每次翻译下一个词语的时候,需要对不同的词语赋予不同的权重,在这种情况下,再解码的时候,就可以同时考虑到整体的信息和局部的信息。

LSTM_4

注意力机制的种类

从初步的调研情况来看,注意力机制有两种方法,一种是基于强化学习(Reinforcement Learning)来做的,另外一种是基于梯度下降(Gradient Decent)来做的。强化学习的机制是通过收益函数(Reward)来激励,让模型更加关注到某个局部的细节。梯度下降法是通过目标函数以及相应的优化函数来做的。无论是 NLP 还是 CV 领域,都可以考虑这些方法来添加注意力机制。

LSTM_5

计算机视觉领域的 Attention 部分论文整理

下面将会简单的介绍几篇近期阅读的计算机视觉领域的关于注意力机制的文章。

Look Closer to See Better:Recurrent Attention Convolutional Neural Network for Fine-grained Image Recognition

在图像识别领域,通常都会遇到给图片中的鸟类进行分类,包括种类的识别,属性的识别等内容。为了区分不同的鸟,除了从整体来对图片把握之外,更加关注的是一个局部的信息,也就是鸟的样子,包括头部,身体,脚,颜色等内容。至于周边信息,例如花花草草之类的,则显得没有那么重要,它们只能作为一些参照物。因为不同的鸟类会停留在树木上,草地上,关注树木和草地的信息对鸟类的识别并不能够起到至关重要的作用。所以,在图像识别领域引入注意力机制就是一个非常关键的技术,让深度学习模型更加关注某个局部的信息。

RA_CNN_1

在这篇文章里面,作者们提出了一个基于 CNN 的注意力机制,叫做 recurrent attention convolutional neural network(RA-CNN),该模型递归地分析局部信息,从局部的信息中提取必要的特征。同时,在 RA-CNN 中的子网络(sub-network)中存在分类结构,也就是说从不同区域的图片里面,都能够得到一个对鸟类种类划分的概率。除此之外,还引入了 attention 机制,让整个网络结构不仅关注整体信息,还关注局部信息,也就是所谓的 Attention Proposal Sub-Network(APN)。这个 APN 结构是从整个图片(full-image)出发,迭代式地生成子区域,并且对这些子区域进行必要的预测,并将子区域所得到的预测结果进行必要的整合,从而得到整张图片的分类预测概率。

RA_CNN_2

RA-CNN 的特点是进行一个端到端的优化,并不需要提前标注 box,区域等信息就能够进行鸟类的识别和图像种类的划分。在数据集上面,该论文不仅在鸟类数据集(CUB Birds)上面进行了实验,也在狗类识别(Stanford Dogs)和车辆识别(Stanford Cars)上进行了实验,并且都取得了不错的效果。

RA_CNN_4

从深度学习的网络结构来看,RA-CNN 的输入时是整幅图片(Full Image),输出的时候就是分类的概率。而提取图片特征的方法通常来说都是使用卷积神经网络(CNN)的结构,然后把 Attention 机制加入到整个网络结构中。从下图来看,一开始,整幅图片从上方输入,然后判断出一个分类概率;然后中间层输出一个坐标值和尺寸大小,其中坐标值表示的是子图的中心点,尺寸大小表示子图的尺寸。在这种基础上,下一幅子图就是从坐标值和尺寸大小得到的图片,第二个网络就是在这种基础上构建的;再迭代持续放大图片,从而不停地聚焦在图片中的某些关键位置。不同尺寸的图片都能够输出不同的分类概率,再将其分类概率进行必要的融合,最终的到对整幅图片的鸟类识别概率。

因此,在整篇论文中,有几个关键点需要注意:

  1. 分类概率的计算,也就是最终的 loss 函数的设计;
  2. 从上一幅图片到下一幅图片的坐标值和尺寸大小。

只要获得了这些指标,就可以把整个 RA-CNN 网络搭建起来。

大体来说,第一步就是给定了一幅输入图片 X, 需要提取它的特征,可以记录为 W_{c}*X,这里的 * 指的是卷积等各种各样的操作。所以得到的概率分布情况其实就是 p(X) = f(W_{c}*X)f 指的是从 CNN 的特征层到全连接层的函数,外层使用了 Softmax 激活函数来计算鸟类的概率。

第二步就是计算下一个 box 的坐标 (t_{x}, t_{y}) 和尺寸大小 t_{\ell},其中 t_{x}, t_{y} 分别指的是横纵坐标,正方形的边长其实是 2*t_{\ell}。用数学公式来记录这个流程就是 [t_{x}, t_{y}, t_{\ell}] = g(W_{c}*X)。在坐标值的基础上,我们可以得到以下四个值,分别表示 x, y 两个坐标轴的上下界:

t_{x(t\ell)} = t_{x} - t_{\ell}, t_{x(br)} = t_{x} + t_{\ell},

t_{y(t\ell)} = t_{y} - t_{\ell}, t_{y(br)} = t_{y} + t_{\ell}.

局部注意力和放大策略(Attention Localization and Amplification)指的是:从上面的方法中拿到坐标值和尺寸,然后把图像进行必要的放大。为了提炼局部的信息,其实就需要在整张图片 X 的基础上加上一个面具(Mask)。所谓面具,指的是在原始图片的基础上进行点乘 0 或者 1 的操作,把一些数据丢失掉,把一些数据留下。在图片领域,就是把周边的信息丢掉,把鸟的信息留下。但是,有的时候,如果直接进行 0 或者 1 的硬编码,会显得网络结构不够连续或者光滑,因此就有其他的替代函数。

在激活函数里面,逻辑回归函数(Logistic Regression)是很常见的。其实通过逻辑回归函数,我们可以构造出近似的阶梯函数或者面具函数。

sigmoid_1

对于逻辑回归函数 \sigma(x) = 1/(1+e^{-kx}) 而言,当 k 足够大的时候,\sigma(x) \approx 1x \geq 0\sigma(x) \approx 0x<0。此时的逻辑回归函数近似于一个阶梯函数。如果假设 x_{0}<x_{1},那么 \sigma(x-x_{0}) - \sigma(x-x_{1}) 就是光滑一点的阶梯函数,\sigma(x-x_{0}) - \sigma(x-x_{1}) \approx 0x < x_{0} \text{ or } x > x_{1}\sigma(x-x_{0}) - \sigma(x-x_{1}) \approx 1x_{0}\leq x\leq x_{1}

因此,基于以上的分析和假设,我们可以构造如下的函数:X^{attr} = X \odot M(t_{x}, t_{y}, t_{\ell}), 其中,X^{attr} 表示图片需要关注的区域,M(\cdot) 函数就是 M(t_{x}, t_{y}, t_{\ell}) = [\sigma(x-t_{x(t\ell)}) - \sigma(x-t_{x(br)})]\cdot[\sigma(y-t_{y(t\ell)}) - \sigma(y-t_{y(br)})], 这里的 \sigma 函数对应了一个足够大的 k 值。

当然,从一张完整的图片到小图片,在实际操作的时候,需要把小图片继续放大,在放大的过程中,可以考虑使用双线性插值算法来扩大。也就是说:

X_{(i,j)}^{amp} = \sum_{\alpha,\beta=0}^{1}|1-\alpha-\{i/\lambda\}|\cdot|1-\beta-\{j/\lambda\}|\cdot X_{(m,n)}^{att},

其中 m = [i/\lambda] + \alpha, n = [j/\lambda] + \beta\lambda 表示上采样因子,[\cdot], \{\cdot\} 分别表示一个实数的正数部分和小数部分。

在分类(Classification)和排序(Ranking)部分,RA-CNN 也有着自己的方法论。在损失函数(Loss Function)里面有两个重要的部分,第一个部分就是三幅图片的 LOSS 函数相加,也就是所谓的 classification loss,Y^{(s)} 表示预测类别的概率,Y 表示真实的类别。除此之外,另外一个部分就是排序的部分,L_{rank}(p_{t}^{(s)}, p_{t}^{(s+1)}) = \max\{0,p_{t}^{(s)}-p_{t+1}^{(s+1)}+margin\}, 其中 p^{(s)} 表示在第 s 个尺寸下所得到的类别 t 的预测概率,并且最大值函数强制了该深度学习模型在训练中可以保证 p_{t}^{(s+1)} > p_{t}^{(s)} + margin,也就是说,局部预测的概率值应该高于整体的概率值。

L(X) = \sum_{s=1}^{3}\{L_{cls}(Y^{(s)},Y^{*})\} + \sum_{s=1}^{2}\{L_{rank}(p_{t}^{(s)},p_{t}^{(s+1)})\}.

RA_CNN_3

在这种 Attention 机制下,可以使用训练好的 conv5_4 或者 VGG-19 来进行特征的提取。在图像领域,location 的位置是需要通过训练而得到的,因为每张图片的鸟的位置都有所不同。进一步通过数学计算可以得到,t_{\ell} 会随着网络而变得越来越小,也就是一个层次递进的关系,越来越关注到局部信息的提取。简单来看,

\frac{\partial L_{rank}}{\partial t_{x}} \propto D_{top} \odot \frac{\partial M(t_{x},t_{y},t_{\ell})}{\partial t_{x}},

这里的 \odot 表示元素的点乘,D_{top} 表示之前的网络所得到的导数。

x\rightarrow t_{x(t\ell)}\frac{\partial M}{\partial t_{x}}<0;

x \rightarrow t_{x(br)}\frac{\partial M}{\partial t_{x}}>0;

其余情况,\frac{\partial M}{\partial t_{x}}=0.

y\rightarrow t_{y(t\ell)}\frac{\partial M}{\partial t_{y}}<0;

y \rightarrow t_{y(br)}\frac{\partial M}{\partial t_{y}}>0;

其余情况,\frac{\partial M}{\partial t_{y}}=0.

x \rightarrow t_{x(t\ell)}\text{ or } x \rightarrow t_{x(br)}\text{ or } y \rightarrow t_{y(t\ell)}\text{ or } y \rightarrow t_{y(br)}, \frac{\partial M}{\partial t_{\ell}}>0;

其余情况,\frac{\partial M}{\partial t_{\ell}}<0.

因此,t_{\ell} 在迭代的过程中会越来越小,也就是说关注的区域会越来越集中。

RA-CNN 的实验效果如下:

RA_CNN_5.png

 

Multiple Granularity Descriptors for Fine-grained Categorization

这篇文中同样做了鸟类的分类工作,与 RA-CNN 不同之处在于它使用了层次的结构,因为鸟类的区分是按照一定的层次关系来进行的,粗糙来看,有科 -> 属 -> 种三个层次结构。

MC_CNN_1

因此,在设计网络结构的过程中,需要有并行的网络结构,分别对应科,属,种三个层次。从前往后的顺序是检测网络(Detection Network),区域发现(Region Discovery),描述网络(Description Network)。并行的结构是 Family-grained CNN + Family-grained Descriptor,Genus-grained CNN + Genus-grained Descriptor,Species-grained CNN + Species-grained Descriptor。而在区域发现的地方,作者使用了 energy 的思想,让神经网络分别聚焦在图片中的不同部分,最终的到鸟类的预测结果。

MC_CNN_2MC_CNN_3

Recurrent Models of Visual Attention

在计算机视觉中引入注意力机制,DeepMind 的这篇文章 recurrent models of visual attention 发表于 2014 年。在这篇文章中,作者使用了基于强化学习方法的注意力机制,并且使用收益函数来进行模型的训练。从网络结构来看,不仅从整体来观察图片,也从局部来提取必要的信息。

DeepMind_1

DeepMind_2DeepMind_3

整体来看,其网络结构是 RNN,上一个阶段得到的信息和坐标会被传递到下一个阶段。这个网络只在最后一步进行分类的概率判断,这是与 RA-CNN 不同之处。这是为了模拟人类看物品的方式,人类并非会一直把注意力放在整张图片上,而是按照某种潜在的顺序对图像进行扫描。Recurrent Models of Visual Attention 本质上是把图片按照某种时间序列的形式进行输入,一次处理原始图片的一部分信息,并且在处理信息的过程中,需要根据过去的信息和任务选择下一个合适的位置进行处理。这样就可以不需要进行事先的位置标记和物品定位了。

DeepMind_4

正如上图所示,enc 指的是对图片进行编码,r_{i}^{(1)} 表示解码的过程,x_{i} 表示图片的一个子区域。而 y_{s} 表示对图片的预测概率或者预测标签。

 

Multiple Object Recognition with Visual Attention

这篇文章同样是 DeepMind 的论文,与 Recurrent Models of Visual Attention 不同之处在于,它是一个两层的 RNN 结构,并且在最上层把原始图片进行输入。其中 enc 是编码网络,r^{(1)}_{i} 是解码网络,r_{i}^{(2)} 是注意力网络,输出概率在解码网络的最后一个单元输出。

deep_recurrent_attention_model_1

在门牌识别里面,该网络是按照从左到右的顺序来进行图片扫描的,这与人类识别物品的方式极其相似。除了门牌识别之外,该论文也对手写字体进行了识别,同样取得了不错的效果。

deep_recurrent_attention_model_3deep_recurrent_attention_model_2

实验效果如下:

deep_recurrent_attention_model_4.png

总结

本篇 Blog 初步介绍了计算机视觉中的 Attention 机制,除了这些方法之外,应该还有一些更巧妙的方法,希望各位读者多多指教。

参考文献

  1. Look Closer to See Better:Recurrent Attention Convolutional Neural Network for Fine-grained Image Recognition,CVPR,2017.
  2. Recurrent Models of Visual Attention,NIPS,2014
  3. GitHub 代码:Recurrent-Attention-CNN,https://github.com/Jianlong-Fu/Recurrent-Attention-CNN
  4. Multiple Granularity Descriptors for Fine-grained Categorization,ICCV,2015
  5. Multiple Object Recognition with Visual Attention,ICRL,2015
  6. Understanding LSTM Networks,Colah’s Blog,2015,http://colah.github.io/posts/2015-08-Understanding-LSTMs/
  7. Survey on the attention based RNN model and its applications in computer vision,2016

 

时间序列的聚类

在机器学习领域,聚类问题一直是一个非常常见的问题。无论是在传统的机器学习(Machine Learning)领域,还是自然语言处理(Natural Language Processing)领域,都可以用聚类算法做很多的事情。例如在数据分析领域,我们可以把某个物品用特征来描述出来,例如该房子的面积,价格,朝向等内容,然后使用聚类算法来把相似的房子聚集到一起;在自然语言处理领域,通常都会寻找一些相似的新闻或者把相似的文本信息聚集到一起,在这种情况下,可以用 Word2Vec 把自然语言处理成向量特征,然后使用 KMeans 等机器学习算法来作聚类。除此之外,另外一种做法是使用 Jaccard 相似度来计算两个文本内容之间的相似性,然后使用层次聚类(Hierarchical Clustering)的方法来作聚类。

word2vec1

本文将会从常见的聚类算法出发,然后介绍时间序列聚类的常见算法。

机器学习的聚类算法

KMeans — 基于距离的机器学习聚类算法

KMeans 算法的目的是把欧氏空间 \mathbb{R}^{m} 中的 n 个节点,基于它们之间的距离公式,把它们划分成 K 个类别,其中类别 K 的个数是需要在执行算法之前人为设定的。

kmeans1

从数学语言上来说,假设已知的欧式空间点集为 \{x_{1},\cdots,x_{n}\},事先设定的类别个数是 K,当然 K\leq n 是必须要满足的,因为类别的数目不能够多于点集的元素个数。算法的目标是寻找到合适的集合 \{S_{i}\}_{1\leq i\leq K} 使得 argmin_{S_{i}}\sum_{x\in S_{i}}||x-\mu_{i}||^{2} 达到最小,其中 \mu_{i} 表示集合 S_{i} 中的所有点的均值。

上面的 ||\cdot|| 表示欧式空间的欧几里得距离,在这种情况下,除了使用 L^{2} 范数之外,还可以使用 L^{1} 范数和其余的 L^{p},p\geq 1 范数。只要该范数满足距离的三个性质即可,也就是非负数,对称,三角不等式。

层次聚类 — 基于相似性的机器学习聚类算法

层次聚类通常来说有两种方法,一种是凝聚,另外一种是分裂。

hierarchicalclustering1

所谓凝聚,其大体思想就是在一开始的时候,把点集集合中的每个元素都当做一类,然后计算每两个类之前的相似度,也就是元素与元素之间的距离;然后计算集合与集合之前的距离,把相似的集合放在一起,不相似的集合就不需要合并;不停地重复以上操作,直到达到某个限制条件或者不能够继续合并集合为止。

所谓分裂,正好与聚合方法相反。其大体思想就是在刚开始的时候把所有元素都放在一类里面,然后计算两个元素之间的相似性,把不相似元素或者集合进行划分,直到达到某个限制条件或者不能够继续分裂集合为止。

在层次聚类里面,相似度的计算函数就是关键所在。在这种情况下,可以设置两个元素之间的距离公式,例如欧氏空间中两个点的欧式距离。在这种情况下,距离越小表示两者之间越相似,距离越大则表示两者之间越不相似。除此之外,还可以设置两个元素之间的相似度。例如两个集合中的公共元素的个数就可以作为这两个集合之间的相似性。在文本里面,通常可以计算句子和句子的相似度,简单来看就是计算两个句子之间的公共词语的个数。

时间序列的聚类算法

通过以上的描述,如果要做时间序列的聚类,通常来说也有多种方法来做,可以使用基于距离的聚类算法 KMeans,也可以使用基于相似度计算的层次聚类算法。

时间序列的特征提取

之前写过很多时间序列特征提取的方法,无论是常见的时间序列特征,例如最大值,最小值,均值,中位数,方差,值域等内容之外。还可以计算时间序列的熵以及分桶的情况,其分桶的熵指的是把时间序列的值域进行切分,就像 Lebesgue 积分一样,查看落入那些等分桶的时间序列的概率分布情况,就可以进行时间序列的分类。除了 Binned Entropy 之外,还有 Sample Entropy 等各种各样的特征。除了时域特征之外,也可以对时间序列的频域做特征,例如小波分析,傅里叶分析等等。因此,在这种情况下,其实只要做好了时间序列的特征,使用 KMeans 算法就可以得到时间序列的聚类效果,也就是把相似的曲线放在一起。参考文章:时间序列的表示与信息提取。

在提取时间序列的特征之前,通常可以对时间序列进行基线的提取,把时间序列分成基线和误差项。而基线提取的最简单方法就是进行移动平均算法的拟合过程,在这种情况下,可以把原始的时间序列 \{x_{1},\cdots,x_{n}\} 分成两个部分 \{baseline_{1},\cdots,baseline_{n}\}\{residual_{1},\cdots,residual_{n}\}。i.e. x_{i} = baseline_{i} + residual_{i}。有的时候,提取完时间序列的基线之后,其实对时间序列的基线做特征,有的时候分类效果会优于对原始的时间序列做特征。参考文章:两篇关于时间序列的论文。

时间序列的相似度计算

如果要计算时间序列的相似度,通常来说除了欧几里得距离等 L^{p} 距离之外,还可以使用 DTW 等方法。在这种情况下,DTW 是基于动态规划算法来做的,基本想法是根据动态规划原理,来进行时间序列的“扭曲”,从而把时间序列进行必要的错位,计算出最合适的距离。一个简单的例子就是把 y=\sin(x)y=\cos(x) 进行必要的横坐标平移,计算出两条时间序列的最合适距离。但是,从 DTW 的算法描述来看,它的算法复杂度是相对高的,是 O(n^{2}) 量级的,其中 n 表示时间序列的长度。参考文章:时间序列的搜索。

dtw1

如果不考虑时间序列的“扭曲”的话,也可以直接使用欧氏距离,无论是 L^{1}, L^{2} 还是 L^{p} 都有它的用武之地。除了距离公式之外,也可以考虑两条时间序列之间的 Pearson 系数,如果两条时间序列相似的话,那么它们之间的 Pearson 系数接近于 1;如果它们之间是负相关的,那么它们之间的 Pearson 系数接近于 -1;如果它们之间没有相关性,Pearson 系数接近于0。除了 Pearson 系数之外,也可以考虑它们之间的线性相关性,毕竟线性相关性与 Pearson 系数是等价的。参考文章:时间序列的相似性。

除此之外,我们也可以用 Auto Encoder 等自编码器技术对时间序列进行特征的编码,也就是说该自编码器的输入层和输出层是恒等的,中间层的神经元个数少于输入层和输出层。在这种情况下,是可以做到对时间序列进行特征的压缩和构造的。除了 Auto Encoder 等无监督方法之外,如果使用其他有监督的神经网络结构的话,例如前馈神经网络,循环神经网络,卷积神经网络等网络结构,可以把归一化之后的时间序列当做输入层,输出层就是时间序列的各种标签,无论是该时间序列的形状种类还是时间序列的异常/正常标签。当该神经网络训练好了之后,中间层的输出都可以作为 Time Series To Vector 的一种模式。i.e. 也就是把时间序列压缩成一个更短一点的向量,然后基于 COSINE 相似度等方法来计算原始时间序列的相似度。参考文章:基于自编码器的时间序列异常检测算法基于前馈神经网络的时间序列异常检测算法。

总结

如果想对时间序列进行聚类,其方法是非常多的。无论是时间序列的特征构造,还是时间序列的相似度方法,都是需要基于一些人工经验来做的。如果使用深度学习的方法的话,要么就提供大量的标签数据;要么就只能够使用一些无监督的编码器的方法了。本文目前初步介绍了一些时间序列的聚类算法,后续将会基于笔者的学习情况来做进一步的撰写工作。

参考文献

  1. 聚类分析:https://en.wikipedia.org/wiki/Cluster_analysis
  2. Dynamic Time Warping:https://en.wikipedia.org/wiki/Dynamic_time_warping
  3. Pearson Coefficient:https://en.wikipedia.org/wiki/Pearson_correlation_coefficient
  4. Auto Encoder:https://en.wikipedia.org/wiki/Autoencoder
  5. Word2Vec:https://en.wikipedia.org/wiki/Word2vec,https://samyzaf.com/ML/nlp/nlp.html

时间序列的单调性

在时间序列的众多研究方向上,除了时间序列异常检测,时间序列的相似性,时间序列的趋势预测之外,无论是在量化交易领域还是其余领域,时间序列的单调性都是一个重要课题。本文将会对时间序列的单调性作简单的介绍。

连续函数的单调性

导数1

在微积分里面,通常都会研究可微函数的导数,因为导数是反映可微函数单调性的一个重要指标。假设 f(x) 是定义域 (a,b) 上的可导函数,那么某个点 x_{0}\in(a,b) 的导数则定义为:

f'(x_{0}) = \lim_{x\rightarrow x_{0}}\frac{f(x)-f(x_{0})}{x-x_{0}}.

对于区间 (a,b) 上的可导函数 f(x) 而言,假设 x_{0}\in (a,b)。如果 f'(x_{0})>0,那么在 x_{0} 的附近,f(x) 是严格单调递增函数;如果 f'(x_{0})<0,那么在 x_{0} 的附近,f(x) 是严格单调递减函数;如果 f'(x_{0})=0,则基于这个事实无法轻易的判断 f(x)x_{0} 附近的单调性。可以参考这两个例子:(1)f(x)=x^{2}x_{0}=0;(2)f(x) = x^{3}x_{0}=0。这两个例子在 x_{0}=0 的导数都是零,并且第一个例子在 x_{0}=0 附近没有单调性,x_{0}=0 就是最小值点;但是第二个例子在 x_{0}=0 处是严格递增的。

平方函数

立方函数

时间序列的单调性

通常来说,时间序列分成上涨下跌两种趋势。如果要严格来写的话,当 x_{n-i+1}<\cdots<x_{n} 时,表示时间序列在 [n-i+1,n] 这个区间内是严格单调递增的;当 x_{n-i+1}>\cdots>x_{n} 时,表示时间序列在 [n-i+1, n] 这个区间内是严格单调下跌的。但是,在现实环境中,较难找到这种严格递增或者严格递减的情况。在大部分情况下,只存在一个上涨或者下跌的趋势,一旦聚焦到某个时间戳附近时间序列是有可能存在抖动性的。所以我们需要给出一个定义,用来描述时间序列在一个区间内的趋势是上升还是下跌。

考虑时间序列 X_{N} = [x_{1},\cdots,x_{N}] 的一个子序列 [x_{i},x_{i+1},\cdots,x_{j}],其中 i<j。如果存在某个 k\in (i,j] 和一组非负实数 [w_{i}, w_{i+1},\cdots,w_{j}] 使得

\sum_{m=k}^{j}w_{m}x_{m} > \sum_{m=i}^{k-1} w_{m}x_{m}, 其中\sum_{m=k}^{j}w_{m} = \sum_{m=i}^{k-1}w_{m}.

就称时间序列 [x_{i},x_{i+1},\cdots,x_{j}]上涨的趋势。

如果存在某个 k\in (i,j] 和一组非负实数 [w_{i}, w_{i+1},\cdots,w_{j}] 使得

\sum_{m=k}^{j}w_{m}x_{m} < \sum_{m=i}^{k-1} w_{m}x_{m}, 其中 \sum_{m=k}^{j}w_{m} = \sum_{m=i}^{k-1}w_{m}.

就称时间序列 [x_{i},x_{i+1},\cdots,x_{j}]下跌的趋势。

时间序列的单调性 — 均线方法

虽然时间序列是离散的,但是却可以把连续函数的思想应用在上面。

假设现在有一个时间序列是 X = [x_{1},\cdots,x_{N}],可以考虑第 i 个点 x_{i} 附近的单调性,按照导数的思想来看就是:当 k\geq 1 时,

(x_{i+k}-x_{i})/((i+k)-i) = (x_{i+k}-x_{i})/k,
(x_{i} - x_{i-k})/(i-(i-k)) = (x_{i} -x_{i-k})/k.

考虑特殊的情形,假设 k=1,当第一个公式大于零时,表示 x_{i+1}>x_{i},i.e. 处于单调上升的趋势中。当第一个公式小于零时,表示 x_{i}<x_{i-1},i.e. 处于单调下降的趋势中。

但是,时间序列有可能有一定的波动性,也就是说时间序列有可能其实看上去是单调上升的,但是有一定的噪声或者毛刺。所以需要想办法处理掉一些噪声和毛刺。于是,就有人提出了以下几种方法。

双均线1

简单的移动平均算法

在时间序列领域,简单的移动平均算法 (Simple Moving Average) 是最常见的算法之一。假设原始的时间序列是 X=[x_{1},\cdots,x_{N}],如果考虑时间戳 n 的移动平均值,那就是考虑从时间戳 n 开始,历史上某个窗口上面的所有序列的平均值,用数学公式来描述就是:

M_{w}(n) = \frac{x_{n-w+1}+\cdots+x_{n}}{w} = \frac{\sum_{j=n-w+1}^{n}x_{j}}{w},

其中 w\geq 1 指的就是窗口的大小。

命题 1. 假设窗口值 \ell>s\geq 1M_{s}(n) - M_{\ell}(n) >0, 表示短线上穿长线,曲线有上涨的趋势;M_{s}(n) - M_{\ell}(n) <0, 表示短线下穿长线,曲线有下跌的趋势。

在这里,短线指的是窗口值 s 所对应的移动平均线,长线指的是窗口值 \ell 所对