前言

最近在学习拐点检测的相关问题，发现 C.Truong 的论文对拐点检测的整个流程和目前主流的一些算法介绍的比较清楚，所以在这里进行了一些记录以及总结，并且对 Truong 发布的 ruptures 库做了一些简单的介绍。如果想要进行更深入的研究，请参考原论文和 ruptures。

文章目录

前言
概览
- 问题定义
- 符号定义
- 研究方法
损失函数
- Piecewise i.i.d. signals
- - 例一：Mean-shifts for normal data
  - 例二：Mean-shifts and scale-shifts for normal data
  - 例三：Count data
- Linear model
- - 例四：Autoregressive model
- Kernel change point detection
- Mahalanobis-type metric
- 损失函数总结
应用于固定 K 的查找算法
- 最优检测 Optimal Detection
- - 简介
  - 伪代码
- 滑动窗口 Window sliding
- - 简介
  - 伪代码
- 二分分割 Binary segmentation
- - 简介
  - 伪代码
- 自底向上分割 Bottom-up segmentation
- - 简介
  - 伪代码
应用于可变 K 值的查找算法
- l 0 l_{0} l0 惩罚项（ l 0 l_{0} l0 penalty）
- - 简介
  - Pelt 算法
  - 伪代码
- l 1 l_{1} l1 惩罚项（ l 1 l_{1} l1 penalty）
- 其他惩罚项
贝叶斯方法
ruptures 库简介
总结

概览

问题定义

拐点检测名为 change point detection，对于一条不平缓的时间序列曲线，认为存在一些时间点 ( t 1 , t 2 , . . . , t k ) (t_{1},t_{2},...,t_{k}) (t1,t2,...,tk)，使得曲线在这些点对应的位置发生突变，这些时间点对应的曲线点称为拐点，在连续的两个拐点之间，曲线是平稳的。

拐点检测算法的质量通过算法输出拐点与实际观测到的拐点的差值绝对值除以样本数来评估。
∀ k , ∣ t ^ k − t k ⋆ ∣ / T , \forall k, \quad\left|\hat{t}_{k}-t_{k}^{\star}\right| / T, ∀k,∣∣t^k−tk⋆∣∣/T,
理想情况下，当样本数T无穷大时，误差应该减少到0，这种性质称为满足渐近一致性 (asymptotic consistency.)
max ⁡ k = 1 , … , K ∗ ∣ t ^ k − t k ⋆ ∣ / T ⟶ p 0 \max _{k=1, \ldots, K^{*}}\left|\hat{t}_{k}-t_{k}^{\star}\right| / T \stackrel{p}{\longrightarrow} 0 k=1,…,K∗max∣∣t^k−tk⋆∣∣/T⟶p0

符号定义

y a . . b y_{a..b} ya..b 表示时间点 a 和 b 之间的时间序列，因此完整信号为 y 0.. T y_{0..T} y0..T。

对于给定的拐点索引 t t t，它的关联分数 associate fraction 称为拐点分数 change point fractions ，公式为：
τ : = t / T ∈ ( 0 , 1 ] . \tau :=t / T \in(0,1]. τ:=t/T∈(0,1].

拐点分数的集合 τ = { τ 1 , τ 2 , … } \boldsymbol{\tau}=\left\{\tau_{1}, \tau_{2}, \ldots\right\} τ={τ1,τ2,…} 写作 ∣ τ |\boldsymbol{\tau} ∣τ|

研究方法

一般思路是构造一个对照函数 contrast function，目标是将对照函数的值最小化。
V ( t , y ) : = ∑ k = 0 K c ( y t k … t k + 1 ) ， V(\mathbf{t}, y) :=\sum_{k=0}^{K} c\left(y_{t_{k} \ldots t_{k+1}}\right)， V(t,y):=k=0∑Kc(ytk…tk+1)，
其中 c ( ⋅ ) c(\cdot) c(⋅) 表示用来测量拟合度 goodness-of-fit 的损失函数 cost function, 损失函数的值在均匀的子序列上较低，在不均匀的子序列上较高。

基于离散优化问题 discrete optimization problem，拐点的总数量记为 K K K ：

如果 K K K 是固定值，估算的拐点值为：
min ⁡ ∣ t ∣ = K V ( t ) , \operatorname{min}_{|\mathbf{t}|=K} V(\mathbf{t}), min∣t∣=KV(t),
如果 K K K 不是固定值，估算的拐点值为：
min ⁡ t V ( t ) + pen ⁡ ( t ) . \min _{t} V(\mathbf{t})+\operatorname{pen}(\mathbf{t}). tminV(t)+pen(t).
其中 p e n ( t ) pen(t) pen(t)为对 t t t 的惩罚项项

在这种方法论下，拐点检测的算法包含以下三个元素：

选择合适的损失函数来测算子序列的均匀程度 homogeneity，这与要检测的变化类型有关
解决离散优化问题
合理约束拐点的数量，确定使用固定的 K K K 还是用 p e n ( ) pen() pen() 来惩罚 penalizing 不固定的数量

损失函数

损失函数与被检测拐点的变化类型有关，下文介绍的损失函数使用在子序列 y a . . b y_{a..b} ya..b 上，其中 1 < = a < b < = T 1 <= a < b <= T 1<=a<b<=T

Piecewise i.i.d. signals

首先介绍一种使用 i . i . d . i.i.d. i.i.d. 变量对信号 y y y 进行建模的一般参数模型。

i . i . d . i.i.d. i.i.d. Independent and identically distributed，独立同分布，指一组随机变量中每个变量的概率分布都相同，且这些随机变量互相独立。

假设 y 1 , . . . , y T y_{1},...,y_{T} y1,...,yT 是独立随机变量，概率分布 probability distribution f ( y t ∣ θ ) f\left(y_{t} | \theta\right) f(yt∣θ) 依赖于向量值参数 θ \theta θ vector-value parameter， θ \theta θ 表示预测到发生突然变化的兴趣点。假设存在一系列真实的拐点 t ⋆ = { t 1 ⋆ , … } \mathrm{t}^{\star}=\left\{t_{1}^{\star}, \ldots\right\} t⋆={t1⋆,…}，并且满足：
y t ∼ ∑ i = 0 K f ( ⋅ ∣ θ k ) , ( t k ⋆ < t ≤ t k + 1 ⋆ ) . y_{t} \sim \sum_{i=0}^{K} f\left(\cdot | \theta_{k}\right) , \left(t_{k}^{\star}<t \leq t_{k+1}^{\star}\right). yt∼i=0∑Kf(⋅∣θk),(tk⋆<t≤tk+1⋆).
在这种场景下，拐点检测通过最大似然估计 maximum likelihood estimation 实施，相关的损失函数是 negative log-likelihood：
c i . i . d ( y a . . b ) : = − sup ⁡ θ ∑ t = a + 1 b log ⁡ f ( y t ∣ θ ) . c_{\mathrm{i.i.d}}\left(y_{a . . b}\right) :=-\sup _{\theta} \sum_{t=a+1}^{b} \log f\left(y_{t} | \theta\right). ci.i.d(ya..b):=−θsupt=a+1∑blogf(yt∣θ).
分布的选择取决于对数据的先验知识。历史上，高斯分布通常用于模拟均值和极差（Scale = 最大值 - 最小值）变化。后来的大部分文献选用了其他参数分布，最突出的是关于一般指数族的分布。

例一：Mean-shifts for normal data

均值漂移模型是拐点检测文献研究最多的模型之一， y 1 , … , y T y_{1}, \dots, y_{T} y1,…,yT 为具有分段常数均值和常数方差的一系列独立正态随机变量。这种情况下，损失函数 c i . i . d . ( ⋅ ) c_{i.i.d.}(\cdot) ci.i.d.(⋅) 为：
c L 2 ( y a . . b ) : = ∑ t = a + 1 b ∥ y t − y ‾ a . . b ∥ 2 2 c_{L_{2}}\left(y_{a . . b}\right) :=\sum_{t=a+1}^{b}\left\|y_{t}-\overline{y}_{a . . b}\right\|_{2}^{2} cL2(ya..b):=t=a+1∑b∥yt−ya..b∥22
其中 y ‾ a . . b \overline{y}_{a..b} ya..b 是子信号 y a . . b y_{a..b} ya..b 的经验均值 empirical mean，这种损失函数一般也称为平方误差损失 quadratic error loss。一般应用于DNA序列和地理信号处理。

例二：Mean-shifts and scale-shifts for normal data

这是一种对均值漂移模型的自然拓展，这种模型也允许方差发生突然变化， y 1 , … , y T y_{1}, \dots, y_{T} y1,…,yT 为具有分段常数均值和常数方差的一系列独立正态随机变量。这种情况下，损失函数 c i . i . d . ( ⋅ ) c_{i.i.d.}(\cdot) ci.i.d.(⋅) 为：
c Σ ( y a . . b ) : = log ⁡ det ⁡ Σ ^ a . . b c_{\Sigma}\left(y_{a .. b}\right) :=\log \operatorname{det} \widehat{\Sigma}_{a .. b} cΣ(ya..b):=logdetΣ a..b
其中 Σ ^ a . . b \widehat{\Sigma}_{a .. b} Σ a..b 是子信号 y a . . b y_{a..b} ya..b 的经验协方差矩阵 empirical covariance matrix，这种损失函数可用于检测随机（不要求高斯）变量的前两个时刻的变化，即使它是加入 c i . i . d . ( ⋅ ) c_{i.i.d.}(\cdot) ci.i.d.(⋅) 的高斯似然。一般应用于股票市场时间序列和医药信息处理。

例三：Count data

y 1 , … , y T y_{1}, \dots, y_{T} y1,…,yT 为具有分段常数速率参数的一系列独立泊松分布 Poisson distributed 随机变量，这种模型的损失函数 c poisson ( ⋅ ) c_{\text {poisson}}(\cdot) cpoisson(⋅) 为：
c Poisson ( y a . . b ) : = − ( b − a ) y ‾ a . . b log ⁡ y ‾ a . . b c_{\text { Poisson }}\left(y_{a . . b}\right) :=-(b-a) \overline{y}_{a .. b} \log \overline{y}_{a .. b} c Poisson (ya..b):=−(b−a)ya..blogya..b
这种模型常应用于建模统计数据。

从理论上讲，变化点检测将一致性属性视为最大似然估计方法。实际上，已经表明，一般情况下，随着样本数量增长为无穷大，估测的最优拐点在概率上收敛于真实的拐点。拐点检测在 c i . i . d . ( ⋅ ) c_{i.i.d.}(\cdot) ci.i.d.(⋅) 下满足渐进一致性。总体而言，参数化的损失函数是有用的并且理论上满足渐进一致性，然而，如果数据不能很好的近似于有效的参数化分布，拐点检测的效果也会受到很大的影响，这迫使我们使用其他的模型。

Linear model

当发生突然变化的变量之间存在线性关系时，可以使用线性模型，这种变化在相关文献中通常被称为结构化变化。这个公式引出了几个著名的模型比如自回归模型 autoregressive (AR) model 和多元回归模型 multiple regressions 等。

在数学上，信号 y 被视为单变量响应变量，两种协变量的信号 x = { x t } t = 1 T x=\left\{x_{t}\right\}_{t=1}^{T} x={xt}t=1T 和 z = { z t } t = 1 T z=\left\{z_{t}\right\}_{t=1}^{T} z={zt}t=1T 分别记为 R p − v a l u e d \mathbb{R}^{p} -valued Rp−valued 和 R q − v a l u e d \mathbb{R}^{q}-valued Rq−valued，假设有一系列真实的拐点 t ⋆ = { t 1 ⋆ , … , t K ⋆ } \mathbf{t}^{\star}=\left\{t_{1}^{\star}, \ldots, t_{K}^{\star}\right\} t⋆={t1⋆,…,tK⋆}，一般线性模型为：
∀ t , t k ⋆ < t ≤ t k + 1 ⋆ , y t = x t ′ u k + z t ′ v + ε t ( k = 0 , … , K ) . \forall t, t_{k}^{\star}<t \leq t_{k+1}^{\star}, \quad y_{t}=x_{t}^{\prime} u_{k}+z_{t}^{\prime} v+\varepsilon_{t} \quad(k=0, \ldots, K). ∀t,tk⋆<t≤tk+1⋆,yt=xt′uk+zt′v+εt(k=0,…,K).
其中 v ∈ R q v \in \mathbb{R}^{q} v∈Rq 和 u k ∈ R p u_{k} \in \mathbb{R}^{p} uk∈Rp 为未知回归参数， ε t \varepsilon_{t} εt 为噪声。此模型也被称为部分结构化变化模型 partial structural change model，y 和 x 的线性关系发生突变，而 y 和 z 的关系是一个常量，如果删除项 z t ′ v z^{\prime}_{t} v zt′v 则为纯结构化变化模型 pure structural change model。相关的损失函数基于最小二乘残差，公式为：
c linear ( y a . . b ) : = min ⁡ u ∈ R p , v ∈ R q ∑ t = a + 1 b ( y t − x t ′ u − z t ′ v ) 2 . c_{\text { linear }}\left(y_{a . . b}\right) :=\min _{u \in \mathbb{R}^{p}, v \in \mathbb{R}^{q}} \sum_{t=a+1}^{b}\left(y_{t}-x_{t}^{\prime} u-z_{t}^{\prime} v\right)^{2}. c linear (ya..b):=u∈Rp,v∈Rqmint=a+1∑b(yt−xt′u−zt′v)2.
存在一个封闭形式的公式，因为它是一个简单的最小二乘回归。

理论上，对于协变量分布 ( x t , z t ) (x_{t},z_{t}) (xt,zt)，噪声分布 ε t \varepsilon_{t} εt，拐点间的距离为 t k + 1 ⋆ − t k ⋆ t_{k+1}^{\star}-t_{k}^{\star} tk+1⋆−tk⋆，估测的拐点在温和假设下概率上渐进回归于真实的拐点。还有其他的渐近性质，例如拐点位置的渐近分布和 V ( t ^ K ) V(\hat{t}_{K}) V(t^K) 的限制分布（可用于设计统计检验）。

上述公式（13）有时也会替换为最小绝对偏差之和，理论上也满足渐近一次性：
c linear , L 1 ( y a … b ) : = min ⁡ u ∈ R p , v ∈ R q ∑ t = a + 1 b ∣ y t − x t ′ u − z t ′ v ∣ . c_{\text {linear}, L_{1}}\left(y_{a \ldots b}\right) :=\min _{u \in \mathbb{R}^{p}, v \in \mathbb{R}^{q}} \sum_{t=a+1}^{b}\left|y_{t}-x_{t}^{\prime} u-z_{t}^{\prime} v\right|. clinear,L1(ya…b):=u∈Rp,v∈Rqmint=a+1∑b∣yt−xt′u−zt′v∣.
线性模型一般应用于金融数据处理。

例四：Autoregressive model

分段自回归模型 Piecewise autoregressive models 也属于线性模型，通过将公式（12）中的参数 x t x_{t} xt 设定为 x t = [ y t − 1 , y t − 2 , … , y t − p ] x_{t}=\left[y_{t-1}, y_{t-2}, \dots, y_{t-p}\right] xt=[yt−1,yt−2,…,yt−p] 并删除项 z t ′ v z^{\prime}_{t} v zt′v，信号 y 实际上是对有序 p 的分段自回归。由此产生的损失函数 c A R ( ⋅ ) c_{AR}(\cdot) cAR(⋅) 能够检测过程中的自回归系数变化。
c A R ( y a . . b ) : = min ⁡ u ∈ R P ∑ t = a + 1 b ( y t − x t ′ u ) 2 . c_{AR}(y_{a..b}) := \min _{u \in \mathbb{R} P} \sum_{t=a+1}^{b}\left(y_{t}-x_{t}^{\prime} u\right)^{2}. cAR(ya..b):=u∈RPmint=a+1∑b(yt−xt′u)2.
分段自回归过程是一种特殊的随时间变化 time-varying 的 ARMA 过程。理论上，即使引入了噪声 ε t \varepsilon_{t} εt（例如动态平均过程），也能保持一致性。因此，在实践中，拐点检测可以应用于完整的ARMA过程。

分段回归模型一般应用于分段弱平稳信号，比如 EEG/ECG 时间序列，fMRI 时间序列（functional magnetic resonance imaging）和语音识别任务。

Kernel change point detection

可以将原始信号通过核函数 kernel function 映射到高维空间中执行拐点检测，这种技术在机器学习的很多技术中都得到了应用，如 SVM（support vector machine）和聚类 clustering。

原始信号 y 可以通过一个核函数 k ( ⋅ , ⋅ ) : R d × R d → R k(\cdot, \cdot) : \mathbb{R}^{d} \times \mathbb{R}^{d} \rightarrow \mathbb{R} k(⋅,⋅):Rd×Rd→R 映射到一个 RKHS（reproducing Hilbert space）空间 H \mathcal{H} H，相关的映射函数 ϕ : R → H \phi : \mathbb{R} \rightarrow \mathcal{H} ϕ:R→H 对任何样本 y s , y t ∈ R d y_{s}, y_{t} \in \mathbb{R}^{d} ys,yt∈Rd 隐式定义为：
ϕ ( y t ) = k ( y t , ⋅ ) ∈ H and ⟨ ϕ ( y s ) ∣ ϕ ( y t ) ⟩ H = k ( y s , y t ) . \phi\left(y_{t}\right)=k\left(y_{t}, \cdot\right) \in \mathcal{H} \quad \text { and } \quad\left\langle\phi\left(y_{s}\right) | \phi\left(y_{t}\right)\right\rangle_{\mathcal{H}}=k\left(y_{s}, y_{t}\right). ϕ(yt)=k(yt,⋅)∈H and ⟨ϕ(ys)∣ϕ(yt)⟩H=k(ys,yt).
RKHS norm ∥ ⋅ ∥ H \|\cdot\|_{\mathcal{H}} ∥⋅∥H 隐式定义为：
∥ ϕ ( y t ) ∥ H 2 = k ( y t , y t ) . \left\|\phi\left(y_{t}\right)\right\|_{\mathcal{H}}^{2}=k\left(y_{t}, y_{t}\right). ∥ϕ(yt)∥H2=k(yt,yt).
直觉上，信号被映射到高位空间 H \mathcal{H} H 上会变成分段常数，于是目标变为检测嵌入信号中的均值漂移，相关的衡量“平均分散” average scatter 的损失函数为：
c kernel ( y a … b ) : = ∑ t = a + 1 b ∥ ϕ ( y t ) − μ ‾ a . . b ∥ H 2 c_{\text { kernel }}\left(y_{a \ldots b}\right) :=\sum_{t=a+1}^{b}\left\|\phi\left(y_{t}\right)-\overline{\mu}_{a . . b}\right\|_{\mathcal{H}}^{2} c kernel (ya…b):=t=a+1∑b∥ϕ(yt)−μa..b∥H2
其中 μ ‾ a . . b ∈ H \overline{\mu}_{a . . b} \in \mathcal{H} μa..b∈H 是嵌入信号的经验均值 { ϕ ( y t ) } t = a + 1 b \left\{\phi\left(y_{t}\right)\right\}_{t=a+1}^{b} {ϕ(yt)}t=a+1b，由于众所周知的“核技巧” kernel trick,，不需要显式计算映射数据样本。实际上，在简单的代数操作之后，kernel 损失函数可以重写为：
c kernel ( y a … b ) : = ∑ t = a + 1 b k ( y t , y t ) − 1 b − a ∑ s , t = a + 1 b k ( y s , y t ) . c_{\text { kernel }}\left(y_{a \ldots b}\right):=\sum_{t=a+1}^{b} k\left(y_{t}, y_{t}\right)-\frac{1}{b-a} \sum_{s, t=a+1}^{b} k\left(y_{s}, y_{t}\right). c kernel (ya…b):=t=a+1∑bk(yt,yt)−b−a1s,t=a+1∑bk(ys,yt).
任何 kernel 都可以插入到这个公式中，但最常用的 kernel 为高斯核 Gaussian kernel（也称为径向基函数 radial basis function）：
c t b f ( y a . . b ) : = ( b − a ) − 1 b − a ∑ s , t = a + 1 b exp ⁡ ( − γ ∥ y s − y t ∥ 2 ) . c_{\mathrm{tbf}}\left(y_{a . . b}\right) :=(b-a)-\frac{1}{b-a} \sum_{s, t=a+1}^{b} \exp \left(-\gamma\left\|y_{s}-y_{t}\right\|^{2}\right). ctbf(ya..b):=(b−a)−b−a1s,t=a+1∑bexp(−γ∥ys−yt∥2).
其中 γ > 0 \gamma>0 γ>0，被称为带宽参数 bandwidth parameter，在该情况下使用线性 kernel 在形式上等同于对分段常数信号使用 c L 2 ( ⋅ ) c_{L_{2}}(\cdot) cL2(⋅)。

总体而言，kernel 拐点检测是一种非参数化的免模型 model-free 方法，可以用于没有先验条件，不满足基础分部形式的分段 i.i.d. 信号。常被应用于生理信号（如脑电波）以及音频时间序列分割任务。

Mahalanobis-type metric

在拐点检测场景中，Mahalanobis-type metric 可以与度量学习算法一起使用来度量信号，也常见于聚类方法。

对于对于任何对称正半定矩阵 symmetric positive semi-definite matrix M ∈ R d × d M \in \mathbb{R}^{d \times d} M∈Rd×d，对于样本 y s , y t y_{s},y_{t} ys,yt，相关的伪度量 pseudo-metric ∥ ⋅ ∥ M \|\cdot\|_{M} ∥⋅∥M 为：
∥ y t − y s ∥ M 2 : = ( y t − y s ) ′ M ( y t − y s ) \left\|y_{t}-y_{s}\right\|_{M}^{2} :=\left(y_{t}-y_{s}\right)^{\prime} M\left(y_{t}-y_{s}\right) ∥yt−ys∥M2:=(yt−ys)′M(yt−ys)
设 y ‾ a . . b \overline{y}_{a . . b} ya..b 为子信号 y a . . b y_{a . . b} ya..b 的经验均值，对应的损失函数定义为：
c M ( y a . . b ) : = ∑ t = a + 1 b ∥ y t − y ‾ a . . b ∥ M 2 c_{M}\left(y_{a . . b}\right) :=\sum_{t=a+1}^{b}\left\|y_{t}-\overline{y}_{a . . b}\right\|_{M}^{2} cM(ya..b):=t=a+1∑b∥yt−ya..b∥M2
本质上，度量矩阵 M M M 为协方差矩阵的逆矩阵，记为 Mahalanobis metric：
M = Σ ^ − 1 M=\widehat{\Sigma}^{-1} M=Σ −1
其中 Σ ^ \widehat{\Sigma} Σ 为信号 y 的经验协方差矩阵。

如果把矩阵分解为 M = U ′ U M=U^{\prime} U M=U′U，可以得到：
∥ y t − y s ∥ M 2 = ∥ U y t − U y s ∥ 2 \left\|y_{t} -y_{s}\right\|_{M}^{2}=\left\|U y_{t}-U y_{s}\right\|^{2} ∥yt−ys∥M2=∥Uyt−Uys∥2
使用 Mahalanobis-type metric 仅涉及样本的线性处理，引入非线性的一种方法是将其与 kernel-based 转换相结合。对于给定的核函数 k ( ⋅ , ⋅ ) : R d × R d ↦ R k(\cdot, \cdot) : \mathbb{R}^{d} \times \mathbb{R}^{d} \mapsto \mathbb{R} k(⋅,⋅):Rd×Rd↦R 以及相关的映射 ϕ ( ⋅ ) : R d ↦ H \phi(\cdot) : \mathbb{R}^{d} \mapsto \mathcal{H} ϕ(⋅):Rd↦H，在特征空间 H \mathcal{H} H 上的 Mahalanobis-type metric ∥ H , M ∥ \|_{\mathcal{H}, M}\| ∥H,M∥ 为：
∥ ϕ ( y s ) − ϕ ( y t ) ∥ H , M 2 = ( ϕ ( y s ) − ϕ ( y t ) ) ′ M ( ϕ ( y s ) − ϕ ( y t ) ) . \left\|\phi\left(y_{s}\right)-\phi\left(y_{t}\right)\right\|_{\mathcal{H}, M}^{2}=\left(\phi\left(y_{s}\right)-\phi\left(y_{t}\right)\right)^{\prime} M\left(\phi\left(y_{s}\right)-\phi\left(y_{t}\right)\right). ∥ϕ(ys)−ϕ(yt)∥H,M2=(ϕ(ys)−ϕ(yt))′M(ϕ(ys)−ϕ(yt)).
其中 M 是一个定义在 H \mathcal{H} H 上的对称正半定矩阵，得到的成本函数为：
c H , M ( y a … b ) : = ∑ t = a + 1 b ∥ y t − y ‾ a . b ∥ H , M 2 c_{\mathcal{H}, M}\left(y_{a \ldots b}\right) :=\sum_{t=a+1}^{b}\left\|y_{t}-\overline{y}_{a . b}\right\|_{\mathcal{H}, M}^{2} cH,M(ya…b):=t=a+1∑b∥yt−ya.b∥H,M2

损失函数总结

应用于固定 K 的查找算法

下面展示当拐点数量K是固定值的时候，如何使用离散优化算法进行拐点检测。指定 t ^ K ( y ) \hat{\mathbf{t}}_{K}(y) t^K(y) 为最小化对照函数 V ( ⋅ ) V(\cdot) V(⋅) 的方法，写作：
t ^ K ( y ) : = arg ⁡ min ⁡ ∣ t ∣ = K V ( t , y ) . \hat{\mathbf{t}}_{K}(y) :=\underset{|\mathbf{t}|=K}{\arg \min } V(\mathbf{t}, y). t^K(y):=∣t∣=KargminV(t,y).

根据信号的上下文信息观测出的最优拐点集合记作 t ^ K \hat{\mathbf{t}}_{K} t^K。

最优检测 Optimal Detection

简介

应用于固定 K 的拐点查找算法是在索引集 { 1 , … , T } \{1, \ldots, T\} {1,…,T} 上的离散优化问题，所有可能的组合基数为 ( T − 1 K ) \left(\begin{array}{c}{T-1} \\ {K}\end{array}\right) (T−1K) ，数量太多使得无法进行全面的遍历所有可能组合，但是却可以有效的使用动态规划 dynamic programming 的方法，本文将此算法记作 Opt 。

实际上，动态规划依赖于目标函数 V ( ⋅ ) V(\cdot) V(⋅) 的附加性质。粗略地说，它包括递归地解决子问题，并满足于以下观察结果：
min ⁡ ∣ t ∣ = K V ( t , y = y 0.. T ) = min ⁡ 0 = t 0 < t 1 < ⋯ < t K < t K + 1 = T ∑ k = 0 K c ( y t k . . t k + 1 ) = min ⁡ t ≤ T − K [ c ( y 0.. t ) + min ⁡ t = t 0 < t 1 < ⋯ < t K − 1 < t K = T ∑ k = 0 K − 1 c ( y t k . . t k + 1 ) ] = min ⁡ t ≤ T − K [ c ( y 0.. t ) + min ⁡ ∣ t ∣ = K − 1 V ( t , y t . . T ) ] . \begin{aligned} \min _{|t|=K} V\left(\mathbf{t}, y=y_{0 .. T}\right) &=\min _{0=t_{0}<t_{1}<\cdots<t_{K}<t_{K+1}=T} \sum_{k=0}^{K} c\left(y_{t_{k} .. t_{k+1}}\right) \\ &=\min _{t \leq T-K}\left[c\left(y_{0 . . t}\right)+\min _{t=t_{0}<t_{1}<\cdots<t_{K-1}<t_{K}=T} \sum_{k=0}^{K-1} c\left(y_{t_{k}.. t_{k+1}}\right)\right] \\ &=\min _{t \leq T-K}\left[c\left(y_{0 . . t}\right)+\min _{|\mathbf{t}|=K-1} V\left(\mathbf{t}, y_{t .. T}\right)\right]. \end{aligned} ∣t∣=KminV(t,y=y0..T)=0=t0<t1<⋯<tK<tK+1=Tmink=0∑Kc(ytk..tk+1)=t≤T−Kmin[c(y0..t)+t=t0<t1<⋯<tK−1<tK=Tmink=0∑K−1c(ytk..tk+1)]=t≤T−Kmin[c(y0..t)+∣t∣=K−1minV(t,yt..T)].
直观来看，公式表明如果由 K-1 个元素组成的最优子信号 { y s } t T \left\{y_{s}\right\}_{t}^{T} {ys}tT 已知，那么计算出第一个拐点是容易的。接下来就可以通过递归的执行上述观察结果来得出完整的集合。结合对二次误差损失函数 quadratic error cost function 的计算，Opt 算法的复杂度为 O ( K T 2 ) \mathcal{O}\left(K T^{2}\right) O(KT2)。使用的损失函数越复杂，计算的复杂度就越高。由于 Opt 算法的复杂度是二次的，因此适合应用于比较短的信号，一般包含一百个样本点左右。

有兴趣的朋友可以阅读Rigaill和Yann在论文中提供的几种降低 Opt 算法计算负担的方法。

伪代码

滑动窗口 Window sliding

简介

滑动窗口算法（下文记为 Win）是一种快速近似的 Opt 的替代算法，该方案依赖于单个变化点检测程序并将其扩展以找到多个变化点。算法实施时，两个相邻的窗口沿着信号滑动。计算第一窗口和第二窗口之间的差异。对给定的损失函数，两个子信号间的差异为：
d ( y a . . t , y t . . b ) = c ( y a . . b ) − c ( y a . . t ) − c ( y t . . b ) ( 1 ≤ a < t < b ≤ T ) . d\left(y_{a . . t}, y_{t . . b}\right)=c\left(y_{a .. b}\right)-c\left(y_{a . . t}\right)-c\left(y_{t .. b}\right) \quad(1 \leq a<t<b \leq T). d(ya..t,yt..b)=c(ya..b)−c(ya..t)−c(yt..b)(1≤a<t<b≤T).
当这两个窗口包含不相似的片段时，计算得的差异值将会很大，可以检测到一个拐点。在离线设置中，计算完整的差异曲线并执行峰值搜索过程以找到拐点索引。

类似的，还有一种方法叫双样本检验（或同质性检验）two-sample test (or homogeneity test) ，这是一种统计假设检验程序，旨在评估两个样本群体的分布是否相同。

Win 的最大好处是复杂性低（与 c L 2 c_{L_{2}} cL2 搭配时为线性时间），易于实现，此外，任何单一拐点检测方法都可以加入到该方案中。有大量可用于差异的渐近分布的参考文献（用于几种成本函数），这对于在峰值搜索过程中找到合适的阈值很有用。然而，Win 通常是不稳定的，因为单个变化点检测是在小区域（窗口）上进行的，降低了它的统计功效。

伪代码

二分分割 Binary segmentation

简介

二分分割（下文记为 BinSeg）也是一种快速近似的 Opt 的替代算法，在拐点检测的应用场景下，BinSeg是最常用的方法之一，因为它概念简单，易于实现。

BinSeg 是一种顺序贪心算法 greedy sequential algorithm，在每次迭代中，执行单个变化点的检测并产生估计 t ^ ( k ) \hat{t}^{(k)} t^(k) 。在下文中，上标 ^k 指的是顺序算法的第 k 步。算法的示意图如下：

第一个估计的拐点 t ^ ( 1 ) \hat{t}^{(1)} t^(1) 为：
t ^ ( 1 ) : = arg ⁡ min ⁡ 1 ≤ t < T − 1 [ c ( y 0.. t ) + c ( y t . . T ) ] ⏟ V ( t = { t } ) . \hat{t}^{(1)} :=\underset{1 \leq t<T-1}{\arg \min } \underbrace{[c\left(y_{0 . . t}\right)+c\left(y_{t . . T}\right)]}_{V(\mathrm{t}=\{t\})}. t^(1):=1≤t<T−1argminV(t={t}) [c(y0..t)+c(yt..T)].
这个操作是贪心的，目的是找出使得总体损失最小的拐点，然后在 t ^ ( 1 ) \hat{t}^{(1)} t^(1) 的位置将信号分成了两部分，再对得到的子信号进行相同的操作，直到找不到更多的拐点为止。表达的正式一些，在 BinSeg 的 k 次迭代以后，下一个估算的拐点 t ^ ( k + 1 ) \hat{t}^{(k+1)} t^(k+1) 属于原始信号的 k + 1 个子片段之一，计算公式为：
t ^ ( k + 1 ) : = arg ⁡ min ⁡ t ∉ { t ^ ( 1 ) , … , t ^ ( k ) } V ( { t ^ ( 1 ) , … , t ^ ( k ) } ∪ { t } ) . \hat{t}^{(k+1)} :=\underset{t \notin\left\{\hat{t}^{(1)}, \ldots, \hat{t}^{(k)}\right\}}{\arg \min } V\left(\left\{\hat{t}^{(1)}, \ldots, \hat{t}^{(k)}\right\} \cup\{t\}\right). t^(k+1):=t∈/{t^(1),…,t^(k)}argminV({t^(1),…,t^(k)}∪{t}).
由于 K 的总数是已知的，因此 BinSeg 需要执行 K 步来估算所有的拐点，贪心检测的过程可以设计为一个双样本检验的过程。在使用二次损失误差 quadratic error loss C L 2 C_{L_{2}} CL2的情况下，BinSeg 的复杂度为 O ( T log ⁡ T ) \mathcal{O}(T \log T) O(TlogT)。

从理论上讲，BinSeg 被证明可以在均值模型下用二次误差 C L 2 C_{L_{2}} CL2 产生渐近一致的估计。更近期的结果表明，BinSeg 只有当任何两个相邻变化点之间的最小间距大于 T 3 / 4 T^{3 / 4} T3/4 时一致。一般来说，BinSeg的输出只是最优解的近似值，尤其是不能精确地检测到接近的拐点，原因在于估计的拐点 t ^ ( k ) \hat{t}^{(k)} t^(k) 不是从均匀分段估计的，并且每个估计都取决于先前的分析。

伪代码

目前有许多算法可以用于提高 BinSeg 的准确性，然而这些提高一般都会牺牲算法原有的简易性。感兴趣的朋友可以参考 Circular binary segmentation 和 wild binary segmentation algorithm。

自底向上分割 Bottom-up segmentation

简介

自底向上分割是 BinSeg 的自然对应，下文记作 BotUp，这种算法同样也是一种近似算法。跟 BinSeg 相反，BotUp 最开始就把源信号分割成很多子信号，然后序列化的把它们合并直到剩下 K 个拐点。在每个步骤中，所有潜在的变化点（分隔相邻子段的索引）按它们分开的段之间的差异进行排序。差异性最小的拐点被删除，由它分割的两个子段被合并。跟 BinSeg 的贪心算法相反，BotUp 被称为一种“慷慨”的算法。

算法示意图为：

对于给定的损失函数，两个子信号的差异性为：
d ( y a . . t , y t . . b ) = c ( y a . . b ) − c ( y a . . t ) − c ( y t . . b ) ( 1 ≤ a < t < b ≤ T ) . d\left(y_{a . . t}, y_{t . . b}\right)=c\left(y_{a .. b}\right)-c\left(y_{a . . t}\right)-c\left(y_{t .. b}\right) \quad(1 \leq a<t<b \leq T). d(ya..t,yt..b)=c(ya..b)−c(ya..t)−c(yt..b)(1≤a<t<b≤T).
BotUp 的优势在于线性复杂度（对T个样本使用二次损失函数 C L 2 C_{L_{2}} CL2）以及概念的简单性。但是它也存在一些问题：

如果实际的拐点不在最开始划分出来的点集中，那么 BotUp 也永远不会考虑它。
在第一次迭代中，融合过程可能是不可靠的，因为计算只是在小的区段上进行，统计意义比较小。
在现有文献中，BotUp 的研究不如 BinSeg，没有理论收敛性研究 theoretical convergence study 可用。

然而，An online algorithm for segmenting time series 的作者发现 BotUp 在十多个数据集上的表现优于 BinSeg。

伪代码

应用于可变 K 值的查找算法

接下来我们讨论拐点总数未知情况下的算法，在这种情况下，我们的主要注意力在于选择惩罚函数 penalty function pen() 以及惩罚等级 penalty level β \beta β。要解决的优化问题为：
min ⁡ t V ( t ) + pen ⁡ ( t ) \min _{t} V(t)+\operatorname{pen}(t) tminV(t)+pen(t)

l 0 l_{0} l0 惩罚项（ l 0 l_{0} l0 penalty）

简介

l 0 l_{0} l0 惩罚项也被称为线性惩罚项，被认为是最流行的惩罚项，一般而言， l 0 l_{0} l0 惩罚项记作 p e n l 0 \mathrm{pen}_{l_{0}} penl0，公式为：
pen ⁡ l 0 ( t ) : = β ∣ t ∣ . \operatorname{pen}_{l_{0}}(\mathrm{t}) :=\beta|\mathrm{t}|. penl0(t):=β∣t∣.
其中 β > 0 \beta > 0 β>0，负责权衡复杂性和拟合度 (controls the trade-off between complexity and goodness-of-fit)，直观而言，当出现的拐点越来越多时，它平衡了总成本 V ( ⋅ ) V(\cdot) V(⋅) 减少的幅度。有很多利用这个规则的例子，最著名的是 BIC (Bayesian Information Criterion) 和 AIC (Akaike Information Criterion)。

假设对信号 y 有一个分段 i.i.d. 模型，设定 c = c i . i . d c=c_{i . i . d} c=ci.i.d，对于给定的一组拐点 t，引出 BIC 的约束似然公式为：
2 V ( t ) + p log ⁡ T ∣ t ∣ . 2 V(\mathrm{t})+p \log T|\mathrm{t}|. 2V(t)+plogT∣t∣.
其中 p 为参数空间 parameter space 的维度（如 p = 1 p=1 p=1 表示平均值的单变量变化， p = 2 p = 2 p=2 表示平均值和极差的单变量变化等等），因此，BIC 等同于当 β = p / 2 log ⁡ T \beta=p / 2 \log T β=p/2logT 时的线性惩罚拐点检测：
pen ⁡ B I C ( t ) : = p 2 log ⁡ T ∣ t ∣ . \operatorname{pen}_{B I C}(\mathrm{t}) :=\frac{p}{2} \log T|\mathrm{t}|. penBIC(t):=2plogT∣t∣.
类似的，在固定方差的均值飘逸模型下 mean-shift model with fixed variance， c = c L 2 c=c_{L_{2}} c=cL2，BIC 惩罚项为：
pen ⁡ B I C , L 2 ( t ) : = σ 2 log ⁡ T ∣ t ∣ \operatorname{pen}_{B I C, L_{2}}(\mathrm{t}) :=\sigma^{2} \log T|\mathrm{t}| penBIC,L2(t):=σ2logT∣t∣
AIC 也是类似的，公式为：
pen ⁡ A I C , L 2 ( t ) : = σ 2 ∣ t ∣ \operatorname{pen}_{A I C, L_{2}}(\mathrm{t}) :=\sigma^{2}|\mathrm{t}| penAIC,L2(t):=σ2∣t∣
从理论上讲，拐点估计的一致性在某些情况下是可以证明的。研究最多的模型是均值漂移模型，当满足下列假设时，其估计的拐点数和拐点分数概率渐近收敛：

在无噪声设置中，如果惩罚值以适当的速率收敛到零，
如果噪声是高斯白噪声 Gaussian white noise 并且 p e n = p e n B I C , L 2 \mathrm{pen}=\mathrm{pen}_{\mathrm{BIC}, L_{2}} pen=penBIC,L2，
如果噪声是二阶静止过程（适当降低自相关功能）并且惩罚级别慢慢发散到无穷大（慢于T）

最近的研究表明对于 c = c k e r n e l c = c_{kernel} c=ckernel 的线性惩罚项的拐点检测估计的一致性会以 1 / T 1/T 1/T 的速率收敛到 0。对于其他损失函数（最值得注意的是 c i . i . d . c_{i.i.d.} ci.i.d.）的理论结果仅涉及成本总和到其真实值的收敛，没有关于一致性的结果。

在各种假设下渐近地收敛于概率：

Pelt 算法

实施线性惩罚项拐点检测的一个朴素的实施方案是对足够大的 K m a x K_{max} Kmax，对 K = 1 , . . . , K m a x K = 1, ... ,K_{max} K=1,...,Kmax 都进行执行 Opt，然后找出使得惩罚结果最小的部分，但是这种方案的平方复杂度是难以接受的，因此我们需要实施其他复杂度更低的方法。

Pelt 算法 Pruned Exact Linear Time 可以用于找出 p e n = p e n l 0 pen = pen_{l_{0}} pen=penl0 的准确解，该方法按顺序考虑每个样本，并基于明确的修剪规则决定是否从潜在的拐点集中丢弃它。假设片段长度是从均匀分布中随机选取的，Pelt 的复杂度是 O ( T ) \mathcal{O}(T) O(T)，尽管这种设定是不现实的(会产生短片段)，但 Pelt 依然比上述的朴素启发式快几个数量级。

使用 Pelt 时唯一需要调整的参数只有惩罚级别 β \beta β， β \beta β 的设定至少与样本数和信号维度有关， β \beta β 值越小，算法对变化的感知越敏感，可以选用较小的值来找到更多的拐点，也可以使用较大的值来仅关注显著变化。

伪代码

l 1 l_{1} l1 惩罚项（ l 1 l_{1} l1 penalty）

为了进一步降低具有线性惩罚的变化点检测的计算成本，已经提出了一种将 l 0 l_{0} l0 惩罚项转换为 l 1 l_{1} l1 惩罚项的替代公式。同样的理由也在机器学习的许多发展中起到了至关重要的作用，如稀疏回归，压缩感知，稀疏PCA等，在数值分析和图像去噪中，这种惩罚项也称为总变异正则项。

引入该策略的目的是检测带高斯噪声 Gaussian noise 的分段常量信号中的均值漂移现象，相关的损失函数是 C L 2 C_{L_{2}} CL2，在这种设定下，拐点检测优化问题被记作：
min ⁡ t V ( t ) + β ∑ k = 1 ∣ t ∣ ∥ y ‾ t k − 1 . . t k − y ‾ t k . . t k + 1 ∥ 1 \min _{\mathbf{t}} V(\mathbf{t})+\beta \sum_{k=1}^{|\mathbf{t}|}\left\|\overline{y}_{t_{k-1} .. t_{k}}-\overline{y}_{t_{k} .. t_{k+1}}\right\|_{1} tminV(t)+βk=1∑∣t∣∥∥ytk−1..tk−ytk..tk+1∥∥1
其中 y ‾ t k − 1 . . t k \overline{y}_{t_{k-1} .. t_{k}} ytk−1..tk 是子信号 y t k − 1 . . t k y_{t_{k-1} .. t_{k}} ytk−1..tk 的经验均值，因此， l 1 l_{1} l1 为：
pen ⁡ l 1 ( t ) : = β ∑ k = 1 ∣ t ∣ ∥ y ‾ t k − 1 . . t k − y ‾ t k . . t k + 1 ∥ 1 \operatorname{pen}_{l_{1}}(\mathrm{t}) :=\beta \sum_{k=1}^{|t|}\left\|\overline{y}_{t_{k-1} . . t_{k}}-\overline{y}_{t_{k} .. t_{k+1}}\right\|_{1} penl1(t):=βk=1∑∣t∣∥∥ytk−1..tk−ytk..tk+1∥∥1
有文献指出它实际上等同于凸优化问题 convex optimization problem，对于矩阵 Y : = [ y 1 , … , y T ] ′ ∈ R T × d Y:=\left[y_{1}, \ldots, y_{T}\right]^{\prime} \in \mathbb{R}^{T \times d} Y:=[y1,…,yT]′∈RT×d， S ∈ R T × T − 1 S \in \mathbb{R}^{T \times T-1} S∈RT×T−1 且：
S i j : = { 1 if i > j 0 if i ≤ j S_{i j} :=\left\{\begin{array}{l}{1 \text { if } i>j} \\ {0 \text { if } i \leq j}\end{array}\right. Sij:={1 if i>j0 if i≤j
Y ‾ \overline{Y} Y 和 S ‾ \overline{S} S 表示对矩阵 Y Y Y 和 S S S 进行 centering each column 所得的矩阵，centering each column 表示矩阵的每一列的每个元素减去该列均值后得到的矩阵。那么拐点检测优化问题的公式等价为：
min ⁡ Δ ∈ R ( T − 1 ) × d ∥ Y ‾ − S ‾ Δ ∥ 2 + β ∑ t = 1 T − 1 ∥ Δ t , ∙ ∥ 1 ⏟ ∥ Δ ∥ 1 , 1 \min _{\Delta \in \mathbb{R}^{(T-1) \times d}}\|\overline{Y}-\overline{S} \Delta\|^{2}+\beta \underbrace{\sum_{t=1}^{T-1}\left\|\Delta_{t, \bullet}\right\|_{1}}_{\|\Delta\|_{1,1}} Δ∈R(T−1)×dmin∥Y−SΔ∥2+β∥Δ∥1,1 t=1∑T−1∥Δt,∙∥1
估算的拐点为矩阵 Δ \Delta Δ 的有效行 support rows（不为0的行），矩阵 Δ \Delta Δ 为“跳跃矩阵” jump matrix，包含了信号 S ‾ Δ \overline{S} \Delta SΔ 均值漂移的位置和幅度，这个公式是著名的 Lasso Regression (“least absolute shrinkage and selection operator”)。可以直观看出 l 1 l_{1} l1 惩罚项将许多系数缩小为零，在这种情况下，它意味着它有利于 y 的分段常数近似 piecewise constant approximations。惩罚项 β \beta β 的值越高，允许得到的拐点越少。

在算法上，这种优化执行了最小角度回归 the least-angle regression (Lars)。当拐点总数的上限已知时，这种方法的复杂度为 O ( T log ⁡ T ) \mathcal{O}(T \log T) O(TlogT)。从理论角度来看，估算的拐点分数是渐近一致的。该结果证明了惩罚值 β \beta β 序列的适当收敛，拐点的数量也可以正确的预测。

其他惩罚项

改编的 BIC 标准是一种惩罚项，它不仅依赖于拐点数量，而且还依赖于拐点之间的距离，公式为：
pen ⁡ m B I C ( t ) : = 3 ∣ t ∣ log ⁡ T + ∑ k = 0 ∣ t ∣ + 1 log ⁡ ( t k + 1 − t k T ) . \operatorname{pen}_{\mathrm{mBIC}}(\mathrm{t}) :=3|\mathrm{t}| \log T+\sum_{k=0}^{|\mathrm{t}|+1} \log \left(\frac{t_{k+1}-t_{k}}{T}\right). penmBIC(t):=3∣t∣logT+k=0∑∣t∣+1log(Ttk+1−tk).
pen ⁡ m B I C ( t ) \operatorname{pen}_{\mathrm{mBIC}}(\mathrm{t}) penmBIC(t) 的第二项适用于均匀分布的拐点，在索引非常近的时候取得最大值，例如 t 1 = 1 , t 2 = 2 , … t_{1}=1, t_{2}=2, \ldots t1=1,t2=2,… 这种方式在实践中较难处理，因此一般用 p e n l 0 pen_{l_{0}} penl0 近似。详细内容可阅读文献。

还有一种惩罚项写作：
pen ⁡ L e b ( t ) : = ∣ t ∣ + 1 T σ 2 ( a 1 log ⁡ ∣ t ∣ + 1 T + a 2 ) . \operatorname{pen}_{\mathrm{Leb}}(\mathrm{t}) :=\frac{|\mathrm{t}|+1}{T} \sigma^{2}\left(a_{1} \log \frac{|\mathrm{t}|+1}{T}+a_{2}\right). penLeb(t):=T∣t∣+1σ2(a1logT∣t∣+1+a2).
详细信息可阅读文献。

贝叶斯方法

贝叶斯方法是拐点检测文献中的重要方法，这些方法可以适应关于拐点位置分布的先验知识，被广泛应用于语音识别，大脑图像，音频分割和生物信息学等特定学科。

贝叶斯方法通常假设一个未观察到的离散状态变量 s t ( t = 1 , … , T ) s_{t}(t=1, \ldots, T) st(t=1,…,T)，这个变量控制数据模型的产生。拐点为这个隐藏变量 S t S_{t} St 发生状态变化时的索引，假定信号 y 满足分段 i.i.d，记作：
y t ∼ f ( ⋅ ∣ θ s t ) y_{t} \sim f\left(\cdot | \theta_{s_{t}}\right) yt∼f(⋅∣θst)
其中 f ( ⋅ ∣ θ ) f(\cdot | \theta) f(⋅∣θ) 为对已观测到的参数 θ \theta θ 的概率分布， θ i \theta_{i} θi 是已观测状态 i 的激发参数 emission parameters，观测的联合概率为：
P ( y , [ s 1 , … , s T ] ) = ∏ t = 1 T f ( y t ∣ θ s t ) \mathbb{P}\left(y,\left[s_{1}, \ldots, s_{T}\right]\right)=\prod_{t=1}^{T} f\left(y_{t} | \theta_{s_{t}}\right) P(y,[s1,…,sT])=t=1∏Tf(yt∣θst)
贝叶斯算法和其他方法不同，它通过假设状态变量 S t S_{t} St 的分布并计算出对观测状态序列的最佳解释。一种很著名的贝叶斯模型为隐马尔科夫模型 Hidden Markov Model (HMM)，状态变量被假设为马尔科夫链 Markov chain，也就是说 S t S_{t} St 的值只跟 S t − 1 S_{t-1} St−1 有关。概率可以完全用转移矩阵 transition matrix A 描述：
A i j = P ( s t = j ∣ s t − 1 = i ) , A_{i j}=\mathbb{P}\left(s_{t}=j | s_{t-1}=i\right), Aij=P(st=j∣st−1=i),
HMM 由三个参数定义：转移矩阵 A，激发概率 f ( ⋅ ∣ θ ) f(\cdot | \theta) f(⋅∣θ) 以及初始状态概率 P ( s 1 = θ i ) \mathbb{P}\left(s_{1}=\theta_{i}\right) P(s1=θi)，这种情景下的的拐点检测相当于最大后验 posteriori（MAP）状态序列 s ^ 1 , … , s ^ T \hat{s}_{1}, \dots, \hat{s}_{T} s^1,…,s^T。如果提供参数的先验分布，则通常使用维特比算法 Viterbi algorithm 或通过随机采样方法（诸如蒙特卡罗马尔可夫链 Monte Carlo Markov Chain，MCMC）来执行计算。

HMM 的优势在于特殊场景下优秀的性能表现以及对任何分段 i.i.d. 信号的建模能力，是一种由很多实际应用的经过深思熟虑的方法。然而，HMM依赖于在许多情况下未经验证的若干假设，例如观察被认为是独立的、激发分布是参数化的、通常是高斯混合、不适用于长段等。也有几种相关算法被设计来克服这些限制。如：动态变换线性模型 Switching linear dynamic models 比独立激发模型具有更大的描述能力，它们适用于视觉跟踪的复杂任务。然而，由于需要调整许多参数，因此校准步骤显得更加麻烦。半马尔可夫模型 semi-Markov models 放宽了马尔可夫假设，半马尔可夫过程跟踪当前段的长度，更适合用几个拐点对状态变量建模。Dirichlet 隐马尔可夫模型 Dirichlet process hidden Markov models 非常适合执行拐点检测，隐式计算状态的数量并且从后验分布采样使得计算上变得容易。

ruptures 库简介

目前主要的用于拐点检测的库大部分都是基于 R 语言的，对于 python 而言，一个比较优秀的库是由 Charles Truong 等人开源的 ruptures 。本文简单介绍了 ruptures 的基本特性，更多详细介绍以及使用手册请阅读官方文档。

ruptures 的流程结构如图所示：

ruptures 的主要特性如下：

Search methods

ruptures 库包含了上文提到的主要查找算法，比如动态规划 dynamic programming，基于 l₀ 约束的检测 detection with a l₀ constraint，二元分割 binary segmentation，自下而上的分割 bottom-up segmentation 以及基于窗口的分割 window-based segmentation
Cost function

在 ruptures 库中，可以选用有参损失函数 parametric cost functions 来检测标准统计量 standard statistical quantities 的变化，如均值 mean，规模 scale，维度间的线性关系 linear relationship between dimensions，自回归系数 autoregressive coefficients 等。也可以使用无参损失函数 non-parametric cost functions (Kernel-based or Mahalanobis-type metric) 来检测分布式变化。
Constraints

无论拐点数量是否已知，上述方法都适用，需要注意的是，ruptures 在成本运算 cost budget 和线性惩罚项 linear penalty term 的基础上实现拐点检测。
Evaluation

评估指标可用于定量比较分割，以及用可视化的方法评估算法性能
Input

ruptures 能在可以转化为 Numpy array 的任何信号上实施，包括单因素信号和多因素信号 univariate or multivariate signal
Consistent interface and modularity

离散优化方法和损失函数是拐点检测的两个主要组成部分，这些在 ruptures 中都是具体的对象，因此本库的模块化程度很高，一旦有新的算法产生，即可无缝的集成到 ruptures 中。
Scalability

数据探索通常需要使用不同的参数集多次运行相同方法。为此，实现高速缓存以将中间结果保留在内存中，从而大大降低了在相同信号上多次运行相同算法的计算成本。

ruptures 还为具有速度限制的用户提供了信号采样和设置拐点最小距离的可能

总结

本文回顾了几种用于拐点检测的算法，描述了算法中的损失函数，搜索方法以及拐点数量约束。并简单介绍了基于 python 的拐点检测库 ruptures。拐点检测模块还有很多其他的优秀算法没有在文中展示，有兴趣更加深入的朋友可以参考 C.Truong 的论文第七章 Summary Table 中的文献整理或者的论文，比如 instance on Bayesian methods，in-depth theoretical survey 和 numerical comparisons in controlled settings。

拐点检测常用算法介绍相关推荐

拐点检测常用算法介绍（Ruptures）
拐点检测常用算法介绍最近在学习拐点检测的相关问题, 发现 C.Truong 的论文对拐点检测的整个流程和目前主流的一些算法介绍的比较清楚,所以在这里进行了一些记录以及总结,并且对 Truong 发 ...
拐点检测常用算法总结
目录概览问题定义符号定义研究方法损失函数概览问题定义拐点检测名为 change point detection,对于一条不平缓的时间序列曲线,认为存在一些时间点 ( t 1 , t 2 ...
根据大小分割大文本_场景文本检测—CTPN算法介绍
SIGAI特约作者:沪东三哥原创声明:本文为SIGAI 原创文章,仅供个人学习使用,未经允许,不得转载,不能用于商业目的. 其它机器学习.深度学习算法的全面系统讲解可以阅读<机器学习-原理.算 ...
数学建模竞赛常用算法介绍及对应国赛获奖论文分类整理分享
数学建模竞赛中应当掌握的算法: 数学建模国赛每年的题型都类似,除非是个人专业性很强,否则作者不太建议选择华为出的题,剩余的题型每年都类似,是有迹可循的,毕竟站在巨人的肩膀上看的更远.下面就介绍一些数模 ...
目标检测 - 主流算法介绍 - 从RCNN到DETR
目标检测是计算机视觉的一个非常重要的核心方向,它的主要任务目标定位和目标分类. 在深度学习介入该领域之前,传统的目标检测思路包括区域选择.手动特征提取.分类器分类.由于手动提取特征的方法往往很难满足目 ...
目标检测-EfficientDet算法介绍
文章与视频资源多平台更新微信公众号|知乎|B站|头条:AI研习图书馆深度学习.大数据.IT编程知识与资源分享,欢迎关注,共同进步~ 1. 引言文章:EfficientDet: Scalable ...
异常值检测常用算法及案例
异常值检测常用方法对历史数据进行异常值检测,对突发情况或者异常情况进行识别,避免因为异常值导致预测性能降低,并对其进行调整便于后续预测. 一.3-sigma原则异常值检测 3-Sigma原则又称为拉 ...
车道线检测相关算法介绍
车道线检测是计算机视觉领域的一个重要应用,常见的车道线检测算法包括以下几种: 1.基于边缘检测的算法该算法基于边缘检测原理,先对图像进行灰度化处理,然后使用Canny边缘检测算法提取边缘信息.最后, ...
超宽带定位中的TOA/TDOA两种最常用算法介绍
关注.星标公众号,直达精彩内容来源参考:华星智控整理:李肖遥 UWB最近一年都很火,自从苹果手机使用了UWB技术之后,可以说引领一波芯片热潮,而我参与到定位通信的项目中也有几年了,深深感受到做技术 ...

拐点检测常用算法介绍

前言