# Change-in-Slope Optimal Partitioning Algorithm
in a Finite-Size Parameter Space

Vincent Runge
Université Paris-Saclay, CNRS, Univ Evry,
Laboratoire de Mathématiques et Modélisation d’Evry,
91037, Evry-Courcouronnes, France111E-mail<EMAIL_ADDRESS>Marco
Pascucci
Université Paris-Saclay, CNRS, Univ Evry,
Laboratoire de Mathématiques et Modélisation d’Evry,
91037, Evry-Courcouronnes, France  Nicolas Deschamps de Boishebert
Université Paris-Saclay, CNRS, Univ Evry,
Laboratoire de Mathématiques et Modélisation d’Evry,
91037, Evry-Courcouronnes, France

###### Abstract

We consider the problem of detecting change-points in univariate time series
by fitting a continuous piecewise linear signal using the residual sum of
squares. Values of the inferred signal at slope breaks are restricted to a
finite set of size $m$. Using this finite parameter space, we build a dynamic
programming algorithm with a controlled time complexity of $O(m^{2}n^{2})$ for
$n$ data points. Some accelerating strategies can be used to reduce the
constant before $n^{2}$. The adapted classic inequality-based pruning is
outperformed by a simpler “channel” method on simulations. Besides, our finite
parameter space setting allows an easy introduction of constraints on the
inferred signal. For example, imposing a minimal angle between consecutive
segment slopes provides robustness to model misspecification and outliers. We
test our algorithm with an isotonic constraint on an antibiogram image
analysis problem for which algorithmic efficiency is a cornerstone in the
emerging context of mobile-health and embedded medical devices. For this
application, a finite state approach can be a valid compromise.

Keywords: Multiple change-point detection, change in slope, pruned dynamic
programming, isotonic constraint, unimodal constraint, robust inference

##  1 Introduction

Detecting change-points in time series is a long-standing problem in
statistics that have been tackled in many different ways since the 1950s.
Originally developed for quality control in manufacturing [25], it has spread
to the most modern sciences such as genomics [14, 5], neuroscience [22, 7] or
climate change [3, 35] among many others. Due to the wide range of possible
modeling assumptions and problem settings, change-point detection remains
today an active scientific field, especially challenging for large data [24].

###  1.1 Change in Mean

Last decades, researchers have been principally focused on the change-in-mean
problem, which consists in inferring a piecewise constant signal from the
data. This multiple change-point problem deals with finding the location and
the number of changes in a time series $y_{1:n}$ of length $n$. Among the
$2^{n-1}$ possible change-point vectors, the problem requires to select one
vector for which a criterion is minimized or close to the minimum. Algorithmic
efficiency is then a central challenge in order to rapidly infer a good
change-point vector candidate.

Binary segmentation (BS) is one of the best known algorithm [28, 29] for
solving this problem. It returns an approximate solution in time
$O(n\log(n))$. Many other approaches have been proposed: see [33] for a review
and [9] for some comparisons. A well-known alternative method is to use a
dynamic programming approach [4] called optimal partitioning (OP) [18], which
returns the best change-point vector minimizing the criterion222The problem
can have many best change-point vectors with the same minimal value given a
criterion. However, this multiple solution is an event of zero probability
with real-valued time series.. Its principle consists in finding the best last
segment in consecutive truncated time series $y_{1:t}$ for $t$ from $1$ to
$n$. Time complexity of this algorithm is of order $O(n^{2})$ as for each
$y_{1:t}$ we test the $t$ available positions (candidates) for potential best
last change-point. In the rest of this paper we will focus on approaches using
dynamic programming algorithms.

In recent years, pruning ideas were proposed to reduce the set of candidates
to consider for each truncated time series $y_{1:t}$, so that the complexity
of OP became comparable with BS in some simulation studies. Two main classes
of pruning are available: inequality-based pruning implemented in the PELT
algorithm [21] and functional pruning [26, 23] implemented in the FPOP
algorithm. PELT is based on a recursion on the position of the last change and
remains close to OP (only a conditional expression is added). PELT is
efficient only when the number of change-points is proportional to data
length. FPOP is based on a recursion on the value of the last segment
parameter and is more difficult to code: we need to update a functional cost
at any step, that is a continuous piecewise quadratic function (in Gaussian
model). However, this strategy achieves the best possible pruning and the log-
linear time complexity on many data-sets [23].

###  1.2 Change in Slope

All these algorithms were developed to solve the change-in-mean problem. The
change-in-slope problem is much less studied, even though many applications
need such models [35, 19, 30, 32]. This problem consists in fitting a
continuous piecewise linear signal based on minimizing the residual sum of
squares of each segment. A few approaches directly tackle this problem [2, 8].
The latter is based on an FPOP-like algorithm and introduces a functional cost
parametrized by the continuous variable $\phi$ being the inferred value of the
last considered data. Using this cost, they derive a one-parametric functional
update and eventually get a time complexity between $O(n^{2})$ and $O(n^{3})$
(see simulations in Subsection 5.2).

The guideline of our work was the construction of a simplified efficient
algorithm for change-in-slope problem, capable of returning a close-to-optimal
change-point vector in a controlled time complexity of order $O(n^{2})$.
Furthermore, we were also interested in the possibility to easily constrain
the inference as for change-in-mean algorithm FPOP [15]. Our algorithm can
restrict the inference to an increasing signal but also force a minimal angle
value between consecutive segment slopes. This latter option provides
robustness to model misspecification and outliers, which is a desired features
for many applications.

Our algorithm, called slopeOP, is an original PELT-like OP algorithm for
change-in-slope problem, based on a finite set (of size $m$) of allowed values
for the inferred signal at slope breaks, which we called states. Given a time
series $y_{1:n}$ of size $n$, like OP, our algorithm searches for the best
last change-point in $y_{1:t}$ for $t$ from $1$ to $n$. But, differently from
OP, it looks for it in a finite set of couples (state,time) of size $m\times
t$, where the state is in the finite states set. Summing up for all states and
for $t$ from $1$ to $n$, we get an overall running time of order
$O(m^{2}n^{2})$.

The use of states is not only a bypass for reducing running time and enforce
constraints. This restricted inference is also the right level of complexity
for some applications as for example in online beat detection, for which the
musician only needs integer beat-per-minute information. Other applications in
embedded or resource-limited systems require time efficency. In this case,
reducing the number of states is a possible trade-off between running time and
precision, as we demonstrate with our application to antibiogram image
analysis.

###  1.3 Parameter Estimation

The criterion to optimize contains a penalty that is linear in the number of
change-points. It is classically set to $2k\hat{\sigma}^{2}\log(n)$ for
change-in-mean problem [34, 31], with $\hat{\sigma}^{2}$ the estimated noise
variance and $k$ the number of changes. For change-in-slope problem, theorems
for a similar penalty have been given in asymptotic regime [36]. As variance
estimation is needed, we proposed a robust estimator of the variance based on
an adaptation to the difference estimator of Hall et al. [12]. We provide a
wide simulation study, showing the importance of developing algorithms
dedicated to changes in slope with continuity constraint. We highlight in
particular the difficulty for setting the penalty value for small data size as
this penalty also depends on the signal shape. The minimum angle constraint
helps inferring a signal closer to the true one, even if we use a sub-optimal
penalty value. The minimum-angle constraint can be safely used to remove
outliers in most cases, if adapted to the expected signal smoothness.

###  1.4 Outline

The outline of this paper is as follows. Section 2 introduces the model and
describes the optimal partitioning algorithm for changes in slope with a
simple description of the continuity constraint. We also give the three
proposed constrained inference modes: isotonic, unimodal and minimal angle.
Pruning-like approaches for speeding up the algorithm are exposed in Section
3. Section 4 is dedicated to variance estimation of time series with slopes.
In the simulation study in Section 5 we show the benefit of using the
continuity hypothesis by measuring the mean squared error (MSE) between the
true signal and the inferred one. We also search for the range of penalties
minimizing the MSE for 4 different simulation scenarios. We then compare the
efficiency of the two proposed pruning methods and also compare the time
complexity with the FPOP-like challenging algorithm CPOP [8]. Eventually, we
show how using the minimal-angle constraint enhances the robustness of the
algorithm. In last Section 6 we apply our algorithm to an antibiogram image
analysis problem to find a unique change-point. In this application where non-
decreasing signals are expected, the use of a monotonicity constraint improves
detection precision. In this example, although the true signal is not expected
to have a finite number of states, the finite-state trade-off allows a
sensible gain in time and still a good analysis precision. In the context of
mobile-health applications and hardware with limited computing capacity, where
efficiency is important, such a trade-off could be fundamental.

Our change-in-slope algorithms are available into an R package on CRAN called
slopeOP333https://CRAN.R-project.org/package=slopeOP and on
github444https://github.com/vrunge/slopeOP.

##  2 Change-in-slope Optimal Partitioning

###  2.1 Model and Cost Function

We define the set of states $\mathcal{S}$ as a finite set of accessible real
values for beginning and ending values in inferred segments. When we write
$s=s_{min},\dots,s_{max}$, the variable $s$ goes through all the values of
$\mathcal{S}$ from the smallest one to the biggest one. For computational
efficiency we recommend to have $\\#\mathcal{S}=m<<n$ but this is not
mandatory.

Data are generated by the model

$Y_{t}=s_{i}+\frac{s_{i+1}-s_{i}}{\tau_{i+1}-\tau_{i}}(t-\tau_{i})+\epsilon_{t}\,,\quad
t=\tau_{i}+1,\dots,\tau_{i+1}\,,\quad i=0,\dots,k\,,$

with $0=\tau_{0}<\tau_{1}<\dots<\tau_{k}<\tau_{k+1}=n$,
$s_{0},\dots,s_{k+1}\in\mathcal{S}$ and
$\epsilon_{t}\sim\mathcal{N}(0,\sigma^{2})$ identically and independently
distributed (iid).

This modeling imposes to define the cost for fitting data $y_{(\tau+1):t}$,
$\tau<t$, with linear interpolation from value $s_{1}\in\mathcal{S}$ to value
$s_{2}\in\mathcal{S}$ as a residual sum of squares:

$\mathcal{C}(y_{(\tau+1):t},s_{1},s_{2})=\sum_{i=\tau+1}^{t}\bigg{(}y_{i}-\Big{(}s_{1}+(s_{2}-s_{1})\frac{i-\tau}{t-\tau}\Big{)}\bigg{)}^{2}\,.$
(1)

The value $s_{1}$ is “unseen” as the cost $(y_{\tau}-s_{1})^{2}$ obtained at
index $\tau$ is not present in the summation. Note that this cost function
doesn’t use a minimisation as in PELT: we cannot guarantee that splitting a
segment into two segments would result in a smaller overall cost. This
property is at the core of the PELT pruning rule. However, this choice of cost
function will give us a simple update rule including the continuity condition
for consecutive segments.

###  2.2 Optimization Problem

The slopeOP problem consists in finding the optimal partitioning of a time
series $y_{1:n}=(y_{1},\dots,y_{n})\in\mathbb{R}^{n}$ of size $n$ into $k+1$
consecutive segments $(y_{\tau_{i}+1},\dots,y_{\tau_{i+1}})$ optimizing a
penalized risk

$Q_{n}^{slope}=\min_{\begin{subarray}{c}\tau=(\tau_{1},\dots,\tau_{k})\in{\mathbb{N}}^{k}\\\
\tau_{0}=0\,,\,\tau_{k+1}=n\\\
(s_{0},\dots,s_{k+1})\in\mathcal{S}^{k+2}\end{subarray}}\sum_{i=0}^{k}\big{\\{}\mathcal{C}(y_{(\tau_{i}+1):\tau_{i+1}},s_{i},s_{i+1})+\beta\big{\\}}-\beta\,,$
(2)

where the states defined inside the cost function provide the continuity
constraint between successive segments and $\beta$ is a positive penalty. This
penalty controls the amount of evidence we need to add a change: the greater
this quantity, the less is $k$. We emphasise that $k$ is an unknown quantity.
The penalty value is often set to $2\hat{\sigma}^{2}\log(n)$ [34, 31, 36] with
$\hat{\sigma}^{2}$ an estimation of the variance $\sigma^{2}$.

###  2.3 Dynamic Programming Algorithm

To address the continuity constraint in an efficient algorithm we introduce
the function

$Q_{t}:\left\\{\begin{array}[]{lll}\mathcal{S}&\to&\mathbb{R}\\\
u&\mapsto&Q_{t}(u)\,,\end{array}\right.$

which is the optimal penalized cost up to position $t$ with the inferred value
at this position equal to $u$. Our goal is then to update the set

$\mathcal{Q}_{t}=\\{Q_{t}(u),u=s_{min},\dots,s_{max}\\}\,,$

at any time step $t\in\\{1,\dots,n\\}$, where $s_{min}$ and $s_{max}$ are
bounds that can be determined in a pre-processing step. We find the objective
value $\displaystyle Q_{n}^{slope}$ with the minimization $\displaystyle
Q_{n}^{slope}=\min_{s_{min}\leq u\leq s_{max}}\big{\\{}Q_{n}(u)\big{\\}}$.

######  Proposition 2.1.

The update rule with continuity constraint takes the form

$Q_{t}(v)=\min_{0\leq t^{\prime}<t}\bigg{(}\min_{s_{min}\leq u\leq
s_{max}}\big{\\{}Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\big{\\}}\bigg{)}\,,$
(3)

where the state $u$ in $Q_{t^{\prime}}$ and the cost function including states
realizes the continuity constraint. At the initial step we have
$Q_{0}(v)=-\beta$ for all $v\in\mathcal{S}$.

###### Proof.

The update rule is obtained by direct computations (we write “vectors $\tau$
and $s$” to mean the constraints in (2)):

$Q_{t}(v)=\min_{\begin{subarray}{c}\text{vectors }\tau\text{ and }s\\\
\tau_{k+1}=t\,,\,s_{k+1}=v\end{subarray}}\sum_{i=0}^{k}\big{\\{}\mathcal{C}(y_{(\tau_{i}+1):\tau_{i+1}},s_{i},s_{i+1})+\beta\big{\\}}-\beta$
$=\min_{\begin{subarray}{c}\text{vectors }\tau\text{ and
}s\end{subarray}}\sum_{i=0}^{k-1}\big{\\{}\mathcal{C}(y_{(\tau_{i}+1):\tau_{i+1}},s_{i},s_{i+1})+\beta\big{\\}}-\beta+\mathcal{C}(y_{(\tau_{k}+1):t},s_{k},v)+\beta$
$=\\!\\!\\!\\!\\!\\!\min_{\begin{subarray}{c}0\leq t^{\prime}<t\\\ s_{min}\leq
u\leq s_{max}\end{subarray}}\\!\\!\\!\min_{\begin{subarray}{c}\text{vectors
}\tau\text{ and }s\\\
\tau_{k}=t^{\prime}\,,\,s_{k}=u\end{subarray}}\sum_{i=0}^{k-1}\big{\\{}\mathcal{C}(y_{(\tau_{i}+1):\tau_{i+1}},s_{i},s_{i+1})+\beta\big{\\}}-\beta+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta$
$=\min_{\begin{subarray}{c}0\leq t^{\prime}<t\\\ s_{min}\leq u\leq
s_{max}\end{subarray}}\left(Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\right)\,.$

∎

The corresponding optimal partitioning algorithm based on (3) is given in
Algorithm 1.

Input : data $y_{1:n}$, set of states
$\mathcal{S}=\\{s_{min},\dots,s_{max}\\}$ and penalty $\beta>0$

1 $Q$ = matrix of size $(n+1)\times m$

2 $cp$ = matrix of size $n\times m$

3 $U$ = matrix of size $n\times m$

4 for  _$v=s_{min}$ to $s_{max}$_ do

5 $Q(0,v)\leftarrow-\beta$

6 end for

7for  _$t=1$ to $n$_ do

8 for  _$v=s_{min}$ to $s_{max}$_ do

9 $\displaystyle Q(t,v)\leftarrow\min_{0\leq
t^{\prime}<t}\left(\min_{s_{min}\leq u\leq
s_{max}}\big{\\{}Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\big{\\}}\right)$

10 $\displaystyle cp(t,v)\leftarrow\operatorname*{argmin}_{0\leq
t^{\prime}<t}\left(\min_{s_{min}\leq u\leq
s_{max}}\big{\\{}Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\big{\\}}\right)$

11 $\displaystyle U(t,v)\leftarrow\operatorname*{argmin}_{s_{min}\leq u\leq
s_{max}}\left(\min_{0\leq
t^{\prime}<t}\big{\\{}Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\big{\\}}\right)$

12 end for

13

14 end for

return $cp$, $U$ and $Q(n,\cdot)$

Algorithm 1 Change-in-slope Optimal Partitioning (slopeOP)

The backtracking step is a little different from the standard optimal
partitioning algorithm as we need to take into account the states.

Input : Output of Algorithm 1 : $cp$, $U$ and
$\\{Q(n,v)\,,\,v\in\mathcal{S}\\}$

1 $chpts\leftarrow()$

2 $\displaystyle states\leftarrow(\operatorname*{argmin}_{s_{min}\leq v\leq
s_{max}}\\{Q(n,v)\\})$

3 $t\leftarrow n$

4

5while  _$t >0$_ do

6 $chpts\leftarrow(t,chpts)$

7 $states\leftarrow(U(t,states(0)),states)$

8

9 $t\leftarrow cp(t,states(1))$

10

11 end while

return $chpts$ and $states$

Algorithm 2 Backtracking Algorithm

###  2.4 Inference with Constraints

Besides its unambiguous computational time efficiency in $O(m^{2}n^{2})$, one
of the benefit of this finite-state approach is the possibility to simply
enforce some constraints in the inference. The quantity to optimize is then

$Q_{n}^{ctt}=\min_{\begin{subarray}{c}\text{vectors }\tau\text{ and }s\\\
(\tau,s)\in
U_{k+1}\end{subarray}}\sum_{i=0}^{k}\big{\\{}\mathcal{C}(y_{(\tau_{i}+1):\tau_{i+1}},s_{i},s_{i+1})+\beta\big{\\}}-\beta\,,$
(4)

where $U_{k+1}$ is the set of constraints containing the positions (couples)
that can be used in the minimization. We use notation $\nu_{i}$ for a couple
$(\tau_{i},s_{i})$ position-state,
$(\tau,s)=((\tau^{\prime},n),(s^{\prime},v))$ is a sequence of couples with
the last one being $(n,v)$ and $U_{k+1}^{(n,v)}$ is the set $U_{k+1}$ of
sequences of couples with a last couple equal to $(n,v)$. A dynamic
programming algorithm solving exactly (4) can be build if in $U_{k+1}^{(n,v)}$
the information over the path $(\nu_{1},...,\nu_{k})$ can be conveyed
alongside this path to the current $\nu_{k}$ position. That is:

$U_{k+1}^{(n,v)}=\\{(\tau,s)\,|\,(\tau^{\prime},s^{\prime})\in
U_{k}^{\nu_{k}}\,,\,f(\nu_{1},...,\nu_{k},(n,v))=1\\}$
$=\\{(\tau,s)\,|\,(\tau^{\prime},s^{\prime})\in
U_{k}^{\nu_{k}}\,,\,g(M(\nu_{k}),\nu_{k},(n,v))=1\\}\,,$

with $g:\mathbb{R}\times(\\{0,\dots,n\\}\times\mathcal{S})^{2}\to\\{0,1\\}$ a
validity test associated to the constraint. Function $f$ is reduced to $g$.
The “memory” function $M:\\{0,\dots,n\\}\times\mathcal{S}\to\mathbb{R}$
summarizes the information of the path $(\nu_{1},...,\nu_{k})$ and is
associated in practice to the cost $Q_{t^{\prime}}(u)$ of the update rule. If
$g$ equals to $0$, the couple $\nu_{k}=(t^{\prime},u)$ can not be considered
in the minimization.

With this definition for $U_{k+1}$ we get the update rule:

$Q_{t}^{ctt}(v)=\min_{\begin{subarray}{c}\tilde{\nu}=(t^{\prime},u)\in\\{0,\dots,t-1\\}\times\mathcal{S}\\\
\nu=(t,v)\,,\,g(M(\tilde{\nu}),\tilde{\nu},\nu)=1\end{subarray}}\bigg{(}Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\bigg{)}\,.$
(5)

In this rule, the positions $\tilde{\nu}$ of the minimization operator are
taken into account according to the $g$ function value depending only on the
current couple $\nu=(t,v)$ and past events memorized by function $M(\cdot)$.
We illustrate this constraint approach by the examples of isotonic, unimodal
and robust inferences.

Isotonic case. The update rule for isotonic constraint is:

$Q_{t}^{iso}(v)=\min_{0\leq t^{\prime}<t}\bigg{(}\min_{s_{min}\leq u\leq
v}\big{\\{}Q_{t^{\prime}}(u)+\mathcal{C}(y_{(t^{\prime}+1):t},u,v)+\beta\big{\\}}\bigg{)}\,,$
(6)

which is close to (3) but with a constrained minimization for variable $u$.
Here $M(\tilde{\nu})=0$ for all $\tilde{\nu}$ and $g(0,\tilde{\nu},\nu)=1$ if
and only if $u\leq v$ for $\tilde{\nu}=(t^{\prime},u)$ and $\nu=(t,v)$. We
illustrate the use of the isotonic constraint in Section 6 with an antibiogram
image analysis problem.

Unimodal case. $M(\tilde{\nu})$ is equal to $1$ if at position
$\tilde{\nu}=(t^{\prime},u)$ no decreasing segment has been yet inferred.
Otherwise $M(\tilde{\nu})=0$. $g(1,\cdot,\cdot)=1$ and

$g(0,\tilde{\nu}=(t^{\prime},u),\nu=(t,v))=\left\\{\begin{array}[]{ll}0&\hbox{if}\,\,u<v\,,\\\
1&\hbox{if}\,\,u\geq v\,.\\\ \end{array}\right.$

Value $M(\nu)$ is determined just after setting $Q(\nu)=Q_{t}(v)$ using the
slope sign information on the segment $(\tilde{\nu},\nu)$.

Minimal angle case. $M(\tilde{\nu})\in\mathbb{R}$ is the slope value of the
last inferred segment at position $\tilde{\nu}$. Here the value
$g(M(\tilde{\nu}),\tilde{\nu},\nu)$ computes the angle deviation in degree
between slope $M(\tilde{\nu})$ ans slope formed by $(\tilde{\nu},\nu)$ If the
obtained value is less that a threshold we have
$g(M(\tilde{\nu}),\tilde{\nu},\nu)=0$ and
$g(M(\tilde{\nu}),\tilde{\nu},\nu)=1$ otherwise. With this constraint, we tend
to be robust to outliers with the right level of smoothness (the threshold).
We also hope to get some kind of robustness to model misspecification (see
simulations in Subsection 5.3).

In the R package slopeOP, we also provide a method returning the best change-
point vector with a given fixed number of segments.

##  3 Accelerating Strategies

As revealed by Algorithm 1 the double loop for variables $t$ and $s$ (lines 7
and 8) is time consuming. We first give a formula for a constant-time
computation of cost (1). Then we propose two accelerating strategies aiming at
reducing the set of values to consider for the search for minimum in matrix
$Q$. The first one is a PELT-like method based on inequalities of type
“$Q_{t^{\prime}}(u)+\mathcal{C}(y_{{(t^{\prime}+1)}:t},u,v)>Q_{t}(v)$”. The
second one considers direct comparisons between elements of the minimization
in (3). We test the efficiency of those two accelerating rules on simulations
in Subsection 5.2. As the constraints can force the choice of non-minimal
elements in (5) these methods can only be implemented in non-constraint
setting (at the exception of the isotonic constraint).

###  3.1 Efficient Cost Computation

######  Proposition 3.1.

The cost function (1) can be computed in constant time with the formula

$\mathcal{C}(y_{(t^{\prime}+1):t},u,v)=S_{t}^{2}-S_{t^{\prime}}^{2}-\frac{2}{t-t^{\prime}}\Big{(}\big{(}ut-
vt^{\prime}\big{)}\big{(}S_{t}^{1}-S_{t^{\prime}}^{1}\big{)}+(v-u)\big{(}S_{t}^{+}-S_{t^{\prime}}^{+}\big{)}\Big{)}$
$\quad\quad\quad\quad\quad\quad+\frac{v^{2}-u^{2}}{2}+\frac{u^{2}+uv+v^{2}}{3}(t-t^{\prime})+\frac{(v-u)^{2}}{6(t-t^{\prime})}\,,$

where $t^{\prime}<t\in\\{0,\dots,n\\}$, $(u,v)\in\mathcal{S}^{2}$ and

$S^{1}_{t}=\sum_{i=1}^{t}y_{i}\quad,\quad
S^{2}_{t}=\sum_{i=1}^{t}y_{i}^{2}\quad\hbox{and}\quad
S^{+}_{t}=\sum_{i=1}^{t}iy_{i}\quad\hbox{for all}\,\,t\in\\{1,\dots,n\\}\,.$

The proof is straightforward by direct (tedious) calculations. The vectors
$(S_{t}^{1})_{t=1,\dots,n}$, $(S_{t}^{2})_{t=1,\dots,n}$ and
$(S_{t}^{+})_{t=1,\dots,n}$ are computed in a pre-processing step in $O(n)$
time and saved. This is a marginal cost compared to the cost for Algorithm 1.
Notice that the $O(m^{2}n^{2})$ time complexity is obtained using the result
of Proposition 3.1.

###  3.2 Inequality-based Pruning

With the following Proposition 3.2 we build a standard inequality-based
pruning that has the specificity to take into account future data. At each
position $v\in\mathcal{S}$, we need to update a set of positions:

$\mathcal{U}(t,v)\subset\\{(t^{\prime},u)\,,\,t^{\prime}\in\\{0,\dots,t-1\\}\,,\,u\in\\{s_{min},\dots,s_{max}\\}\\}\,,$

to get $\mathcal{U}(t+1,v)$. There are $m$ such sets. Therefore, we actually
chose to transfer information at constant state value.

######  Proposition 3.2.

At time $t$ and position $v$, suppose that there exists $t^{\prime}<t$ and
$u\neq v$ such that

$Q_{t^{\prime}}(u)+\mathcal{C}(y_{{(t^{\prime}+1)}:t},u,v)>Q_{t}(v)\,.$ (7)

We define $S_{t}^{t^{\prime}}=\sum_{i=t^{\prime}+1}^{t}(i-t^{\prime})y_{i}$
and coefficients
$(\alpha_{t}^{-},\alpha_{t}^{+},\gamma^{-}_{t},\gamma^{+}_{t})$ such that

$\alpha_{t}^{+}T+\gamma^{+}_{t}\leq\sum_{i=t+1}^{T-1}y_{i}\frac{T-i}{T-t}\leq\alpha_{t}^{-}T+\gamma^{-}_{t}\,,\quad
T=t+1,\dots,n\,.$

We also define for all integers $T>t$ the affine in $T$ functions:

$\left\\{\begin{array}[]{l}\displaystyle
f^{+}(T)=\left(\alpha^{+}_{t}-\frac{u+2v}{6}\right)T+t^{\prime}\frac{u+2v}{6}-\frac{v-u}{12(t-t^{\prime})}+\frac{S_{t}^{t^{\prime}}}{t-t^{\prime}}+\gamma^{+}_{t}\,,\\\
\displaystyle
f^{-}(T)=\left(\alpha^{-}_{t}-\frac{u+2v}{6}\right)T+t^{\prime}\frac{u+2v}{6}-\frac{v-u}{12(t-t^{\prime})}+\frac{S_{t}^{t^{\prime}}}{t-t^{\prime}}+\gamma^{-}_{t}\,.\end{array}\right.$

If $f^{+}(t+1)\geq 0$ and $f^{+}(n)\geq 0$ (case $v-u>0$) or if
$f^{-}(t+1)\leq 0$ and $f^{-}(n)\leq 0$ (case $v-u<0$) then the position
$(t^{\prime},u)$ doesn’t have to be considered for further iterations $T>t$
for computing $Q_{T}(v)$ in Algorithm 1. We remove $(t^{\prime},u)$ from
$\mathcal{U}(T,v)$ for all $T>t$.

###### Proof.

In order to prune we need to force the inequality

$\mathcal{C}(y_{{(t^{\prime}+1)}:T},u,v)\geq\mathcal{C}(y_{{(t^{\prime}+1)}:t},u,v)+\mathcal{C}(y_{{(t+1)}:T},v,v)\,,$

to be true for all $T\geq t+1$ and fixed $t^{\prime},t,u,v$. Using the cost
definition (1), the inequality is equivalent to

$(v-u)\Bigg{(}\sum_{i=t^{\prime}+1}^{t}y_{i}\frac{(i-t^{\prime})(T-t)}{(t-t^{\prime})(T-t^{\prime})}+\sum_{i=t+1}^{T-1}y_{i}\frac{T-i}{T-t^{\prime}}\Bigg{)}$
$\quad\quad\quad\quad\quad\quad\quad\geq(v-u)(T-t)\Bigg{(}\frac{u+2v}{6}+\frac{v-u}{12(t-t^{\prime})(T-t^{\prime})}\Bigg{)}$

and then

$(v-u)g_{t}(T)\geq(v-u)\Bigg{(}(T-t^{\prime})\Big{(}\frac{u+2v}{6}\Big{)}+\frac{v-u}{12(t-t^{\prime})}-\frac{S_{t}^{t^{\prime}}}{t-t^{\prime}}\Bigg{)}$
(8)

with

$g_{t}(T)=\left\\{\begin{array}[]{ll}\displaystyle\sum_{i=t+1}^{T-1}y_{i}\frac{T-i}{T-t}&\hbox{
for }T=t+2,\dots,n\,,\\\ 0&\hbox{ for }T=t+1\,.\end{array}\right.$

As the right hand side of inequality (8) is linear in $T$ we choose the
following strategy: find coefficients
$(\alpha_{t}^{-},\alpha_{t}^{+},\gamma^{-}_{t},\gamma^{+}_{t})$ of an upper
and lower linear approximation such that

$\alpha_{t}^{+}T+\gamma^{+}_{t}\leq
g_{t}(T)\leq\alpha_{t}^{-}T+\gamma^{-}_{t}\,,\quad T=t+1,\dots,n\,.$

We then introduce the linear-in-$T$ functions $f^{-}$ and $f^{+}$ defined in
the Proposition. To prove that the inequality $f^{+}(T)\geq 0$ holds for all
$T=t+1$ to $n$ we only need to have $f^{+}(t+1)\geq 0$ and $f^{+}(n)\geq 0$
due to the linearity of $f^{+}$. With the same argument for $f^{-}$ the result
is proven. ∎

The use of additional parameters (the states) in cost function leads to a more
complex pruning rule than the one for standard PELT [21] as the overall cost
is not guarantee to decrease as the number of change-point increases. However,
in case $u=v$ we have a simpler result:

######  Proposition 3.3.

At time $t$, if there exists $t^{\prime}\in\\{0,\dots,t-1\\}$ such that

$Q_{t^{\prime}}(v)+\mathcal{C}(y_{{(t^{\prime}+1)}:t},v,v)>Q_{t}(v)\,,$ (9)

then the position $t^{\prime}$ doesn’t have to be considered for further
iterations $T>t$ for computing $Q_{T}(v)$ in Algorithm 1.

###### Proof.

We use the fact that we have

$\mathcal{C}(y_{{(t^{\prime}+1)}:t},v,v)+\mathcal{C}(y_{(t+1):T},v,v)=\mathcal{C}(y_{{(t^{\prime}+1)}:T},v,v)\,,$

for any $t^{\prime}<t<T$. Adding the quantity $\mathcal{C}(y_{(t+1):T},v,v)$
to (9) and using this equality we get

$Q_{t^{\prime}}(v)+\mathcal{C}(y_{{(t^{\prime}+1)}:T},v,v)+\beta>Q_{t}(v)+\mathcal{C}(y_{(t+1):T},v,v)+\beta\,,$

which means that $(t^{\prime},v)$ is a choice less optimal than $(t,v)$ for
the minimization in $Q_{T}(v)$. ∎

We have a second pruning rule for v-constant positions which is original in
the sense that a recent position can be pruned by an older one.

######  Proposition 3.4.

At time $t$, if there exists $t^{\prime}\in\\{0,\dots,t-1\\}$ such that

$Q_{t^{\prime}}(v)+\mathcal{C}(y_{{(t^{\prime}+1)}:t},v,v)\leq Q_{t}(v)\,,$
(10)

then the position $t$ doesn’t have to be considered for further iterations
$T>t$ for computing $Q_{T}(v)$ in Algorithm 1.

The proof is the same as for Proposition 3.3.

######  Remark 3.1.

One of the two inequalities (9) and (10) is always satisfied so that the
subset $\mathcal{U}(t,v)\cap\\{(t^{\prime},v)\,,\,t^{\prime}=1,\dots,t-1\\}$
contains at most one element.

###  3.3 Channel Pruning

We propose a new rule for speeding-up the algorithm, called the “channel
method”. It does not prune the position $(t^{\prime},u)$ in the classical
sense, as a non-considered value can be further reintegrated in the
minimization (in a set $\mathcal{U}(t,v)$). However its simplicity could help
to reduce time complexity.

We expect the cost $Q$ to verify some structural properties. In particular,

$q_{t}=\left\\{\begin{array}[]{lll}\mathcal{S}&\to&\mathbb{R}\\\
v&\mapsto&Q_{t}(v)\end{array}\right.$

should have a minimum for $v$ “near” the data, which means that $q_{t}$ is,
for most $t$, decreasing until a state $v_{l}$ and increasing from a state
$v_{u}\geq v_{l}$. The function

$C_{(t^{\prime}+1):t}^{\tilde{v}}=\left\\{\begin{array}[]{lll}\mathcal{S}&\to&\mathbb{R}\\\
v&\mapsto&\mathcal{C}(y_{(t^{\prime}+1):t},v,\tilde{v})\end{array}\right.$

has a global minimum that can be easily found (it’s a quadratic function). The
idea is to study the variations of
$q_{t^{\prime}}+C_{(t^{\prime}+1):t}^{\tilde{v}}$ in order to leave out the
state values for which we know they are away from the argminimum.

######  Proposition 3.5.

We consider the function $q_{t^{\prime}}+C^{\tilde{v}}_{(t^{\prime}+1):t}$
with $t^{\prime}+1<t$. The minimum of this function is inside the interval

$I_{t^{\prime}}^{\tilde{v}}=\mathcal{S}\setminus(S_{min}\cup S_{max})\,,$

with

$\left\\{\begin{array}[]{l}S_{min}=\\{s_{min},\dots,\min(v_{l}^{-1},[v^{*}]^{-1})\\}\,,\\\
S_{max}=\\{\max(v_{u}^{+1},[v^{*}]^{+1}),\dots,s_{max}\\}\,,\end{array}\right.$

where $[\cdot]$ is the nearest “value in set $\mathcal{S}$” operator and
$(\cdot)^{+1}$ (resp. $(\cdot)^{-1}$) denotes the following (resp. preceding)
value in the ordered set $\mathcal{S}$. Function $q_{t^{\prime}}$ is
decreasing on $\\{s_{min},\dots,v_{l}\\}$ and increasing on
$\\{v_{u},\dots,s_{max}\\}$. We also know that the argminimum of
$C^{\tilde{v}}_{(t^{\prime}+1):t}$ is given by formula

$v^{*}=\frac{6}{(t-t^{\prime}-1)(2(t-t^{\prime})-1)}\sum_{i=t^{\prime}+1}^{t}\Big{[}(t-i)y_{i}\Big{]}-\tilde{v}\frac{t-t^{\prime}+1}{2(t-t^{\prime})-1}\,.$

###### Proof.

$v^{*}$ is easy to compute as the argminimum of the quadratic function
$C^{\tilde{v}}_{(t^{\prime}+1):t}$. The set $S_{min}$ (resp. $S_{max}$)
corresponds to the set on which we know that function
$q_{t^{\prime}}+C^{\tilde{v}}_{(t^{\prime}+1):t}$ is decreasing (resp.
increasing) with the next (resp. previous) value in $\mathcal{S}$ giving a
smaller output. ∎

Notice that the closest state search $[v^{*}]$ is a $O(1)$ operation if
$\mathcal{S}$ is made of consecutive integers or a $O(\log(n))$ operation in a
general ordered list. We define the matrix $\mathcal{Q}\in\mathbb{R}^{m\times
n}$ with $\mathcal{Q}_{ij}=Q_{j}(s_{i})$. The name “channel” comes from the
fact that the matrix $\mathcal{Q}$ is restricted by a channel with fixed
$v_{u}$ and $v_{l}$ states for each column of this matrix. For all
$\mathcal{U}(t,v)$, this channel is updated to $I_{t^{\prime}}^{\tilde{v}}$
after computing the value $v^{*}$ and using Proposition 3.5.

##  4 Variance Estimation

With real data-sets the noise level $\sigma^{2}$ in time series is an unknown
quantity that has to be estimated. In change-in-mean problems there exist many
estimators for the variance of a time-series as for example the mean absolute
deviation (MAD) or HALL estimators [12]. They have the property to be slightly
sensitive to change-points when the number of changes is small compared to
data length. In slope problems we face the additional difficulty to remove the
slope effect. We suggest to adapt the HALL estimator to our change-in-slope
problem. In classic mean problem, we have for HALL of order $3$:

$\hat{\sigma}^{2}_{mean}=\frac{1}{n-3}\sum_{j=1}^{n-3}\Big{(}\sum_{k=0}^{3}d_{k}Y_{j+k}\Big{)}^{2}\,,$

with $d_{0}=0.1942$, $d_{1}=0.2809$, $d_{2}=0.3832$ and $d_{3}=-0.8582$. We
apply Hall on the successive differences (“HALL Diff”):

$\hat{\sigma}^{2}_{slope}=\frac{1}{(n-4)\Delta}\sum_{j=1}^{n-4}\Big{(}\sum_{k=0}^{3}d_{k}(Y_{j+k+1}-Y_{j+k})\Big{)}^{2}\,,$

with
$\Delta=d_{0}^{2}+(d_{1}-d_{0})^{2}+(d_{2}-d_{1})^{2}+(d_{3}-d_{2})^{2}+d_{3}^{2}=1.527507$
a normalization coefficient.

To evaluate the quality of the HALL Diff estimator, we simulate $100$ time
series of length $100$ with $2$ signals, one piecewise linear and the other
sinusoidal as in Figure 1.

Figure 1: Examples for the linear and sinusoidal time series with noise
$\sigma=1$ (left) and $\sigma=5$ (right).

In Table 1, we present results for estimating the standard deviation $\sigma$
for $100$ simulated time series for each level of noise ($\sigma=1,\dots,5$).
We easily conclude that the HALL Diff method gives very accurate results, even
with a non-linear (sinusoidal) model. We also have a better precision of the
estimated sigma for HALL Diff than for MAD.

method | signal | $\sigma$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---|---|---
MAD | slope | mean($\hat{\sigma}$) | 1.21 | 2.12 | 3.08 | 4.04 | 5.10
|  | sd.($\hat{\sigma}$) | 0.13 | 0.23 | 0.37 | 0.45 | 0.68
| sinus | mean($\hat{\sigma}$) | 2.16 | 2.58 | 3.37 | 4.25 | 5.24
|  | sd.($\hat{\sigma}$) | 0.15 | 0.27 | 0.40 | 0.50 | 0.62
HALL | slope | mean($\hat{\sigma}$) | 1.84 | 2.51 | 3.37 | 4.31 | 5.28
|  | sd.($\hat{\sigma}$) | 0.047 | 0.12 | 0.22 | 0.26 | 0.37
| sinus | mean($\hat{\sigma}$) | 3.50 | 3.88 | 4.50 | 5.21 | 5.99
|  | sd.($\hat{\sigma}$) | 0.03 | 0.08 | 0.14 | 0.26 | 0.30
HALL Diff | slope | mean($\hat{\sigma}$) | 1.01 | 1.99 | 2.99 | 4.04 | 5.04
|  | sd.($\hat{\sigma}$) | 0.092 | 0.19 | 0.28 | 0.34 | 0.47
| sinus | mean($\hat{\sigma}$) | 1.02 | 1.98 | 3.00 | 3.97 | 4.98
|  | sd.($\hat{\sigma}$) | 0.09 | 0.17 | 0.24 | 0.35 | 0.41
Table 1: Variance estimation with MAD, HALL and HALL Diff estimators

##  5 Simulation Study

Our simulation study is split into three parts. We first compare the mean
squared error (MSE) between the true signal and the inferred one with many
different algorithms and study the impact of choosing different $\beta$
penalty values. In the second part, we select our main competitor CPOP to
challenge it in terms of computational efficiency. Finally, we consider
misspecified time series with an heavy-tailed noise and explore the capacity
of slopeOP with a minimal angle between segments to infer a good model on a
range of penalties.

###  5.1 MSE Competition with Other Algorithms

We consider the 5 following methods:

  * •

OP-2D : we fit the data with the PELT algorithm but the cost function on a
segment is the residual sum of squares between data and a linear regression
(2-dimensional fit). In that case, we don’t have any continuity constraint;

  * •

FPOP : The standard efficient gaussian FPOP algorithm [23] used on differenced
data $z_{t}=y_{t+1}-y_{t}$, $t=1,\dots,n-1$. The continuity is obtained by
construction;

  * •

RFPOP : Same algorithm as the previous FPOP but with robust biweight loss
(with robust parameter $K=3\sigma$) [10];

  * •

CPOP : A FPOP-like algorithm for changes in slope with continuity constraint
[8]. The set of states is here infinite-dimensional ($\mathbb{R}$) and not
finite as for slopeOP;

  * •

slopeOP : Our finite-state change-in-slope OP algorithm.

Figure 2: The first row presents 3 time-series of size $50$, $300$ and $500$
that are obtained by the black hat-shaped signal with an additive Gaussian
noise with $\sigma=12$. For a range of values for
$b=\frac{\beta}{\sigma^{2}\log(n)}$ we plot the average log(MSE) curve for
algorithms FPOP, RFPOP and slopeOP for $10$ simulated time series. The first
two algorithms give performances far away from slopeOP. In the third row we
show an example of an inferred signal by FPOP with a $b$ value corresponding
to the minimum in log(MSE) curves. In last row we plot the segments obtained
by the biggest possible $b$ value before having a unique segment with the same
data. We observe that we fail to obtain only $2$ segments and that we
obviously obtain an unrealistic result.

The first idea to detect change in slope consists in looking at changes in
data $z_{t}=y_{t+1}-y_{t}$ with any well-known change-in-mean algorithm. We
use here the efficient FPOP and robust FPOP algorithms for the inference. We
highlight in Figure 2 the fact that with a simple hat-shaped model and a low
level of noise, the FPOP algorithms give bad performances as the number of
datapoint increases. Therefore, we left out theses two methods and study only
the behavior of three algorithms: OP-2D, CPOP and slopeOP.

Figure 3: The plain black line is the true signal. For each scenario we
simulate time series with an additive Gaussian noise of increasing standard
deviation ($\sigma=3,\sigma=12,\sigma=24$). Each plot shows a realisation of
this procedure with $500$ data points. The number of segments for scenarios
$1$ to $4$ is respectively $2$, $7$, $6$ and $8$.

We now consider 4 scenarios of increasing complexity with 3 levels of noise as
shown in Figure 3. The beginning and ending values of all segments are
integers in $\\{0,\dots,60\\}$ and we choose a set of states
$\mathcal{S}=\\{-10,\dots,70\\}$ to reduce side effects in limit points ($0$
and $60$).
We simulate $30$ time series of length $500$ for each of the couple (scenario,
noise) and draw the graph of the MSE (between the true signal and the inferred
signal) with respect to $b=\frac{\beta}{\sigma^{2}\log(n)}$, where $\sigma$ is
the chosen known level of noise. For each of the $30$ time series, the $3$
algorithms are run for $b=0.1,0.2,\dots,4.9,5$ and each point of the MSE curve
is the mean over the $30$ independent simulations. For the smallest MSE in
each of the $12$ MSE curves, we present the $30$ corresponding results in the
boxplots. The last column shows the obtained segments for the best (smallest
MSE) of the $30$ results in boxplots (for OP-2D and slopeOP).

We discuss the results for the intermediate level of noise ($\sigma=12$)
exposed in Figure 4, other results are exposed in A. For all scenarios CPOP
and slopeOP give very close and similar good results compared to OP-2D. The
lack of continuity constraint explains the higher MSE value for OP-2D in
scenarios $1$, $3$ and $4$. For all scenarios and all noise levels, the
optimal penalty in our simulations with $n=500$ is close to $b=2$ for slopeOP
and CPOP as expected in asymptotic regime [36].

Figure 4: Results for scenarios $1$ to $4$ with noise $\sigma=12$. Each point
of the MSE curve is the mean over $30$ independent simulations. For the
smallest MSE in each of the $12$ MSE curves, we present the $30$ corresponding
results in the boxplots. The last column shows the obtained segments for the
best of the $30$ results in boxplots for OP-2D and slopeOP.

###  5.2 Computational Efficiency

In this section we consider hat-shaped time series with a signal from $10$ to
$50$ and $61$ integers states ($\mathcal{S}=\\{0,\dots,60\\}$) for inference.

We first compare time efficiency for the two accelerating methods: channel-
based rule and inequality-based pruning. In Figure 5 we easily understand that
the channel rule is much more efficient for all level of noise. It may be
possible to get similar pruning efficiency in some very particular situations
but we conclude with these simulations that the default rule should always be
the channel method.

Figure 5: Without pruning the dynamic programming algorithm needs to scan
$m^{2}n(n-1)$ elements for finding mimima at each step. We compare the
proportion of elements that the algorithm scans for channel-based and
inequality-based prunings with a hat-shaped signal with different noise
levels. For all noise levels, the channel method is much more efficient.

We run the slopeOP and CPOP algorithms over $10$ simulations for each $n$ from
$n=100$ to $n=1500$ by a step of size $10$. In Figure 6, we plot the mean time
for each algorithm and two levels of noise in log-log scale. The results
confirms the quadratic complexity for slopeOP and a complexity for CPOP
between $n^{2}$ and $n^{3}$ (closer to $n^{3}$ with a higher level of noise).
For example, with $n=1500$ and $\sigma=24$ we have a mean time of $43s$ for
CPOP and $21s$ for slopeOP.

Figure 6: Computational time of CPOP versus slopeOP for hat-shaped data with
noise level $\sigma=3$ and $\sigma=24$ in log-log scale. The coefficient $q$
in complexity $O(n^{q})$ is equal to $2.68$ and $2.97$ for CPOP for $\sigma=3$
and $\sigma=24$, respectively. For slopeOP we get $1.93$ and $2.01$ with the
channel pruning option.

###  5.3 A Minimum Angle Between Consecutive Segments

We generate time series in scenario $1$ of the hat-shaped signal with an
heavy-tailed noise (a Student distribution with a degree of freedom equal to
$3$) with $\sigma=24$. Simple computations give an angle of about $153^{o}$
between the two segments. We choose to fix a minimal angle between segments to
$130^{o}$. We explore the value of the MSE returned by the standard slopeOP
algorithm compared with the same algorithm with smoothing option for a range
of penalty values. We also return the number of inferred segments. Each data-
point is the mean over $10$ simulations.

Results show that the MSE as well as the number of segments is lower, whatever
the penalty value. With a misspecified model or in presence of outliers, the
minimum angle option (when we expect no small angle) improves the result of
slopeOP. The penalty value has a reduced impact on the algorithm, which can be
important for applications when the time series to analysis does not have a
Gaussian noise structure.

Figure 7: For a range of penalties from $\beta=0.5\sigma^{2}\log(n)$ to
$\beta=2.5\sigma^{2}\log(n)$, $10$ hat-shaped time series have been simulated
with $n=500$ for each penalty value and we ran slopeOP with a smoothing option
(minAngle = 130) and without constraint (std). On the left panel A, we compare
the MSE between the inferred signal and the true one for the two methods; on
the right panel B, we compare the number of inferred segments for the two
methods. Results highlight the stability with respect to the penalty of the
algorithm with minimal angle option.

##  6 Application to Antibiogram Image Analysis

Antibiotic Susceptibility Testing (AST) is a microbiology test used to guide
antibiotic prescription in bacterial diseases by determining the
susceptibility of bacteria to different antibiotics.

In disk diffusion AST, or antibiogram, cellulose disks impregnated with
specific antibiotics are placed on the surface of Petri-dishes previously
inoculated with a microorganism. Agar plates are incubated so that bacteria
can grow everywhere, except around the cellulose disks that contain
antibiotics to which the bacteria are susceptible. Then the diameter of the
zone of non-growth (inhibition) surrounding each antibiotic disk is measured
and compared to known minimum diameters to determine the susceptibility (the
comparison is done with a sensibility of one millimeter)[20]. Diameters can be
read by eye with a ruler, but several automatic reading systems exist which
process pictures of incubated plates [27, 17].

The inhibition zone boundary is usually clear and easy to measure, but it can
sometimes be hazy and fuzzy. In most cases, zone edges should be read at the
point of complete inhibition[1], but determining this point can be
challenging.

Several image processing procedures have been proposed for measuring
inhibition diameters[13, 11, 6], but none of them focuses specifically on the
issue of fuzzy borders. We propose the use of slopeOP to measure inhibition
diameters in these cases. To our knowledge, a MSE-minimizing optimal
segmentation algorithm has never been used to measure AST’s inhibition
dimeters.

We designed a processing pipeline that uses slopeOP to read inhibition
diameters from the picture of an AST and tested it on one hundred inhibition
zones presenting fuzzy borders. Fifty inhibition zones were taken from
standard Mueller-Hinton growth medium antibiograms, the other fifty form blood
agar antibiograms, which are darker and less contrasted.

Figure 8: Picture of a part of an antibiogram showing an antibiotic disk and
its inhibition zone (left). On the right a plot of the intensity radial
profile and the segmentation with slopeOP with isotonic constraint.

After normalizing the intensity of the AST picture, a sub-region of each
inhibition zone was selected and centered on the the antibiotic disk. Then the
radial intensity profile $I(r)$ was extracted: we recorded $I$, the average
intensity value of the 10 most intense pixels (the scale of the images is
$\simeq 10\frac{pixel}{mm}$) lying at fixed distance $r$ from the sub-region
center (see figure 8). This method is a simple way to detect bacteria even if
the inhibition zone in not a perfect disk. The average profile size in the
data-set is 230 point.

slopeOP was used to segment the radial profile and determine the inhibition
zone radius. For the segmentation with slopeOP we take the signal starting at
3.5mm, i.e. just after the plateau corresponding to the pellet disk (from this
point, the ideal signal is supposed to be isotonic).

slopeOP is used with $\mathcal{S}=\\{I(r)_{min},...,I(r)_{max}\\}$ and
$\beta=255$ (the max gray value of a 8-bit image). We tested both with and
without the isotonic constraint (as the observed signal is supposed to be
isotonic).

After the segmentation, we measure the inhibition radius $r_{inib}$ as the
distance corresponding to the first change-point. The inhibition diameter is
obtained as $d_{inhib}=2\times r_{inhib}$.

As a comparison, we measured the diameters with the method suggested by
Gavoille et al.[11] which uses a student t-test thus considering both the
intensity and the texture of the pixels around the antibiotic disk.

The so measured diameters were compared with manual (by eye) measurements, by
calculating the absolute difference.

|  | diameter diff quantiles
---|---|---
AST type | algorithm | 25% | 50% | 75%
M-H | slopeOP | 0.29 | 0.91 | 2.92
| slopeOP-iso | 0.24 | 0.69 | 1.71
| t-test | 0.50 | 1.50 | 3.25
blood | slopeOP | 0.72 | 1.72 | 3.00
| slopeOP-iso | 0.47 | 1.28 | 2.38
| t-test | 2.00 | 9.50 | 16.75
Table 2: Agreement of auto and manual reading. We observe the distribution of
$\Delta{d}=|d_{a}-d_{c}|$ the absolute difference between the automatically
measured diameter $d_{a}$ and the control value $d_{c}$.

The results (Table 2) show a good consistency with the manual measurements. In
the case of Mueller-Hinton antibiorgams, half of the diameter differences
between our procedure and control are below the test sensibility (1mm) and 75%
are smaller than 1.5mm, which is comparable with the acknowledged inter-
operator variability due to the subjectivity of the measurement ($\pm$ 1mm)
[16]. The results on the blood agar ASTs are slightly worse because of the
lower contrast in the images and a consequent decreased signal-noise ratio. In
both cases the isotonic constraint yields improved results.

|  | diameter diff quantiles | speed | mean number of changepoints
---|---|---|---|---
|  | 25% | 50% | 75% |  |
agar | state density |  |  |  |  |
M-H | 1 | 0.24 | 0.69 | 1.71 | 1.00 | 4.25
| 1/2 | 0.19 | 0.64 | 1.84 | 2.97 | 4.33
| 1/4 | 0.21 | 0.77 | 1.55 | 7.74 | 4.24
| 1/8 | 0.38 | 0.83 | 2.57 | 17.92 | 4.25
blood | 1 | 0.47 | 1.28 | 2.38 | 1.00 | 3.92
| 1/2 | 0.43 | 1.04 | 1.84 | 2.76 | 3.90
| 1/4 | 0.70 | 1.52 | 2.34 | 6.69 | 3.94
| 1/8 | 1.07 | 2.05 | 3.83 | 14.53 | 3.96
Table 3: Performance of slopeOP-isotonic at various states density. The
calculation speed is reported as the inverse of the average execution time
(normalized to density=1). Changepoints=average number of changepoints found.
NOTE: Mean execution time for slopeOP-iso1 is 5ms @ 2,3 GHz Intel Core i5.

Successively we tested the performance of slopeOP with isotonic constraint
when reducing the density of states by 2,4 and 8. We defined the states space
as $\mathcal{S}=\\{I_{min},I_{min}+s,I_{min}+2s,...,I_{max}\\}$ and repeated
the measurements on the data-set at each value of $s=2,4,8$. The results are
reported in Table 3 and show a neat improvement in the execution speed by
reducing the states density. Although the precision of the measurement
decreases with decreasing states density, the measured diameter are still in
reasonable accord with the control values.

## Acknowledgement

We thank Guillem Rigaill (LaMME Evry, IPS2 Paris-Saclay) for joint mentoring
of our intern (co-author Nicolas Deschamps de Boishebert), for the carefull
reading of this work and the many usefull feedback.

##  Appendix A Simulation Results for the $4$ Scenarios

Figure 9: Results for scenarios $1$ to $4$ with noise $\sigma=3$ Figure 10:
Results for scenarios $1$ to $4$ with noise $\sigma=24$

## References

  * [1] Reading guide - eucast disk diffusion method for antimicrobial susceptibility testing.  page 3. EUCAST, 2020.
  * [2] R. Baranowski, Y. Chen, and P. Fryzlewicz.  Narrowest-over-threshold detection of multiple change points and change-point-like features.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 81(3):649–672, 2019.
  * [3] C. Beaulieu, J. Chen, and J. L. Sarmiento.  Change-point analysis as a tool to detect abrupt climate variations.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1962):1228–1249, 2012.
  * [4] R. E. Bellman and S. Dreyfus.  Applied dynamic programming, princeton unlv.  Press, Princeton, NJ, 1962.
  * [5] H. Cao and W. Biao Wu.  Changepoint estimation: another look at multiple testing problems.  Biometrika, 102(4):974–980, 2015.
  * [6] L. F. Costa, E. S. da Silva, V. T. Noronha, I. Vaz-Moreira, O. C. Nunes, and M. M. de Andrade.  Development of an automatic identification algorithm for antibiogram analysis.  Computers in biology and medicine, 67:104–115, 2015.
  * [7] I. Cribben and Y. Yu.  Estimating whole-brain dynamics by using spectral clustering.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 66(3):607–627, 2017.
  * [8] P. Fearnhead, R. Maidstone, and A. Letchford.  Detecting changes in slope with an l0 penalty.  Journal of Computational and Graphical Statistics, pages 1–11, 2018\.
  * [9] P. Fearnhead and G. Rigaill.  Relating and comparing methods for detecting changes in mean.  Stat, page e291.
  * [10] P. Fearnhead and G. Rigaill.  Changepoint detection in the presence of outliers.  Journal of the American Statistical Association, 9 2017.
  * [11] A. Gavoille, B. Bardy, and A. Andremont.  Measurement of inhibition zone diameter in disk susceptibility tests by computerized image analysis.  Computers in biology and medicine, 24(3):179–188, 1994.
  * [12] P. Hall, J. Kay, and D. Titterinton.  Asymptotically optimal difference-based estimation of variance in nonparametric regression.  Biometrika, 77(3):521–528, 1990.
  * [13] G. Hejblum, V. Jarlier, J. Grosset, and A. Aurengo.  Automated interpretation of disk diffusion antibiotic susceptibility tests with the radial profile analysis algorithm.  Journal of clinical microbiology, 31(9):2396–2401, 1993.
  * [14] T. Hocking, G. Rigaill, and G. Bourque.  Peakseg: constrained optimal segmentation and supervised penalty learning for peak detection in count data.  In International Conference on Machine Learning, pages 324–332, 2015.
  * [15] T. D. Hocking, G. Rigaill, P. Fearnhead, and G. Bourque.  A log-linear time algorithm for constrained changepoint detection.  arXiv preprint arXiv:1703.03352, 2017.
  * [16] M. Hombach, F. P. Maurer, T. Pfiffner, E. C. Böttger, and R. Furrer.  Standardization of operator-dependent variables affecting precision and accuracy of the disk diffusion method for antibiotic susceptibility testing.  Journal of clinical microbiology, 53(12):3864–3869, 2015.
  * [17] M. Hombach, R. Zbinden, and E. C. Böttger.  Standardisation of disk diffusion results for antibiotic susceptibility testing using the sirscan automated zone reader.  BMC microbiology, 13(1):225, 2013.
  * [18] B. Jackson, J. D. Scargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, and T. T. Tsai.  An algorithm for optimal partitioning of data on an interval.  IEEE Signal Processing Letters, 12(2):105–108, 2005.
  * [19] S. Jamali, P. Jönsson, L. Eklundh, J. Ardö, and J. Seaquist.  Detecting changes in vegetation trends using time series segmentation.  Remote Sensing of Environment, 156:182–195, 2015.
  * [20] J. H. Jorgensen and J. D. Turnidge.  Susceptibility test methods: Dilution and disk diffusion methods*.  In Manual of Clinical Microbiology, Eleventh Edition, pages 1253–1273. American Society of Microbiology, 2015.
  * [21] R. Killick, P. Fearnhead, and I. A. Eckley.  Optimal detection of changepoints with a linear computational cost.  Journal of the American Statistical Association, 107(500):1590–1598, 2012.
  * [22] L. Koepcke, G. Ashida, and J. Kretzberg.  Single and multiple change point detection in spike trains: comparison of different cusum methods.  Frontiers in systems neuroscience, 10:51, 2016.
  * [23] R. Maidstone, T. Hocking, G. Rigaill, and P. Fearnhead.  On optimal multiple changepoint algorithms for large data.  Statistics and Computing, 27(2):519–533, Mar 2017.
  * [24] National Research Council.  Frontiers in massive data analysis, 2013.
  * [25] E. S. Page.  Continuous inspection schemes.  Biometrika, 41(1/2):100–115, 1954.
  * [26] G. Rigaill.  A pruned dynamic programming algorithm to recover the best segmentations with 1 to k_max change-points.  Journal de la Société Française de Statistique, 156(4):180–205, 2015.
  * [27] M. Sanchez, B. S. Del Saz, E. Loza, F. Baquero, and R. Canton.  Evaluation of the osiris video reader system for disk diffusion susceptibility test reading.  Clinical microbiology and infection, 7(7):352–357, 2001.
  * [28] A. J. Scott and M. Knott.  A cluster analysis method for grouping means in the analysis of variance.  Biometrics, 30(3):507–512, 1974.
  * [29] A. Sen and M. S. Srivastava.  On tests for detecting change in mean.  The Annals of statistics, pages 98–108, 1975.
  * [30] Y. Sowa, A. D. Rowe, M. C. Leake, T. Yakushi, M. Homma, A. Ishijima, and R. M. Berry.  Direct observation of steps in rotation of the bacterial flagellar motor.  Nature, 437(7060):916, 2005.
  * [31] S. Tickle, I. Eckley, P. Fearnhead, and K. Haynes.  Parallelization of a common changepoint detection method.  Journal of Computational and Graphical Statistics, 29(1):149–161, 2020.
  * [32] R. Todd Ogden and G. L. Collier.  On detecting and modeling deterministic drift in long run sequences of tapping data.  Communications in Statistics-Theory and Methods, 28(3-4):977–987, 1999.
  * [33] C. Truong, L. Oudre, and N. Vayatis.  Selective review of offline change point detection methods.  Signal Processing, 167:107299, 2020.
  * [34] Y.-C. Yao.  Estimating the number of change-points via schwarz’criterion.  Statistics & Probability Letters, 6(3):181–189, 1988.
  * [35] M. Yu and E. Ruggieri.  Change point analysis of global temperature records.  International Journal of Climatology, 39(8):3679–3688, 2019.
  * [36] C. Zheng, I. A. Eckley, and P. Fearnhead.  Consistency of a range of penalised cost approaches for detecting multiple changepoints.  arXiv preprint arXiv:1911.01716, 2019.