# A log-linear time algorithm for constrained changepoint detection

Toby Dylan Hocking<EMAIL_ADDRESS>
Guillem Rigaill<EMAIL_ADDRESS>
Paul Fearnhead<EMAIL_ADDRESS>
Guillaume Bourque<EMAIL_ADDRESS>
###### Abstract

Changepoint detection is a central problem in time series and genomic data.
For some applications, it is natural to impose constraints on the directions
of changes. One example is ChIP-seq data, for which adding an up-down
constraint improves peak detection accuracy, but makes the optimization
problem more complicated. We show how a recently proposed functional pruning
technique can be adapted to solve such constrained changepoint detection
problems. This leads to a new algorithm which can solve problems with
arbitrary affine constraints on adjacent segment means, and which has
empirical time complexity that is log-linear in the amount of data. This
algorithm achieves state-of-the-art accuracy in a benchmark of several genomic
data sets, and is orders of magnitude faster than existing algorithms that
have similar accuracy. Our implementation is available as the PeakSegPDPA
function in the coseg R package, https://github.com/tdhock/coseg

###### Contents

  1. 1 Introduction
    1. 1.1 Contributions and organization
  2. 2 Related work
  3. 3 Isotonic regression and changepoint models
    1. 3.1 Classical isotonic regression
    2. 3.2 Segment neighborhood changepoint model
    3. 3.3 Reduced isotonic regression
  4. 4 Functional pruning algorithms for constrained changepoint models
    1. 4.1 Equivalent optimization space
    2. 4.2 Dynamic programming update rules
    3. 4.3 Example and comparison with unconstrained case
    4. 4.4 The PeakSeg up-down constraint
    5. 4.5 General affine inequality constraints between adjacent segment means
  5. 5 Results on peak detection in ChIP-seq data
    1. 5.1 Empirical time complexity in ChIP-seq data
    2. 5.2 Test accuracy in ChIP-seq data
  6. 6 Discussion and conclusions
  7. 7 Reproducible Research Statement
  8. 8 Acknowledgements
  9. A Proof of optimality of dynamic programming algorithm
  10. B Algorithm pseudocode
    1. B.1 GPDPA for reduced isotonic regression
    2. B.2 MinLess algorithm
    3. B.3 Implementation details
    4. B.4 Penalized version of reduced isotonic regression
    5. B.5 Generalized Functional Pruning Optimal Partitioning Solvers

##  1 Introduction

Changepoint detection is a central problem in fields such as finance or
genomics, where $n$ data are gathered in a sequence over time or space. Many
models define the optimal changepoints using maximum likelihood, resulting in
a discrete optimization problem. Multiple changepoint detection models seek
the optimal $K$ segments ($K-1$ changes), which amounts to optimizing
likelihood parameters over a space that contains $O(n^{K-1})$ discrete
arrangements of changepoints. In general this problem can be solved in
$O(Kn^{2})$ time using the original dynamic programming algorithm of Auger and
Lawrence (1989). Recently proposed pruning techniques reduce the number of
changepoints considered by the algorithm, thus reducing time complexity to
$O(Kn\log n)$ while maintaining optimality (Rigaill, 2010; Johnson, 2013;
Maidstone et al., 2016).

In “unconstrained” changepoint models, there are no contraints between model
parameters on separate segments. To regularize and obtain a more interpretable
model, it is often desirable to introduce constraints between model parameters
before and after changepoints. For example, the main problem that motivates
this paper is peak detection in ChIP-seq data, which provide noisy
measurements of protein binding or modification throughout a genome (Bailey et
al., 2013). An up-down constrained changepoint detection model has been shown
to achieve state-of-the-art peak detection accuracy in ChIP-seq data (Hocking
et al., 2015). The constraints of this model force an up change in the segment
mean parameter after each down change, and vice versa. The fastest existing
solver for this problem is the Constrained Dynamic Programming Algorithm
(CDPA), which has two issues. First, it is a heuristic algorithm that is not
guaranteed to recover the optimal solution. Second, its $O(Kn^{2})$ quadratic
time complexity is too slow for use on large data sets. In this paper we
propose a new algorithm that fixes both of these issues.

###  1.1 Contributions and organization

We begin by discussing previous research into pruning techniques for solving
unconstrained changepoint detection problems (Section 2), then state the
constrained optimization problems (Section 3). Our main contribution is
Section 4, which generalizes the functional pruning technique of Rigaill
(2010), thus providing a new Generalized Pruned Dynamic Progamming Algorithm
(GPDPA) for solving a class of constrained changepoint detection problems. We
show that the GPDPA achieves state-of-the-art speed and accuracy in genomic
data with several different labeled patterns (Section 5), then conclude by
discussing the significance of our contributions (Section 6).

| No pruning | Functional pruning
---|---|---
Unconstrained | Dynamic Programming Algorithm (DPA) | Pruned DPA (PDPA)
| Optimal solution, $O(Kn^{2})$ time | Optimal solution, $O(Kn\log n)$ time
| Auger and Lawrence (1989) | Rigaill (2010); Johnson (2013)
Up-down constrained | Constrained DPA (CDPA) | Generalized Pruned DPA (GPDPA)
| Sub-optimal solution, $O(Kn^{2})$ time | Optimal solution, $O(Kn\log n)$ time
| Hocking et al. (2015) | This paper
Table 1: Our contribution is the Generalized Pruned Dynamic Programming
Algorithm (GPDPA), which uses a functional pruning technique to compute the
constrained optimal $K-1$ changepoints in a sequence of $n$ data, in $O(Kn\log
n)$ time on average.

##  2 Related work

There are many efficient algorithms available for computing the optimal $K-1$
changepoints in $n$ data points. Auger and Lawrence (1989) proposed an
$O(Kn^{2})$ algorithm for computing the sequence of models with $1,\dots,K$
segments. Jackson et al. (2005) consider a related approach, which introduces
a penalty for each changepoint, rather than fixing the number of changepoints.
Their $O(n^{2})$ algorithm computes the single model for a given penalty
constant $\lambda$. Both of these algorithms recover the optimal solution, and
follow from using dynamic programming updates (Bellman, 1961) to recursively
compute the maximum likelihood from 1 to $n$ data points. Alternatively there
are methods which are computationally faster but are not guaranteed to find
the optimal segmentation. The most popular of these is the binary segmentation
algorithm which has $O(Kn)$ worst-case time complexity (Scott and Knott,
1974). An L1 relaxation of this problem is known as the fused lasso signal
approximator, for which efficient solvers also exist (Hoefling, 2010).

Several pruning methods have been recently proposed in order to reduce time
complexity, while maintaining optimality. Rigaill (2010) and Johnson (2011)
independently discovered a functional pruning technique, which results in
algorithms with $O(n\log n)$ average time complexity. Killick et al. (2011)
proposed an inequality pruning technique, which results in an algorithm with
average time complexity from $O(n)$ to $O(n^{2})$, depending on the number of
changes. Maidstone et al. (2016) provides a clear discussion on the
differences between the two pruning techniques.

All algorithms discussed thus far are for solving problems with no constraints
between adjacent segment mean parameters, but there are many examples of
constrained changepoint detection models. Rather than searching all possible
changepoints and likelihood parameters, the idea is to use a constraint in
order to search a smaller, more interpretable model space. For example,
Haiminen et al. (2008) propose an $O(Kn^{2})$ algorithm for unimodal
regression, which enforces no up changes after the first down change. Hocking
et al. (2015) proposed an $O(Kn^{2})$ algorithm for peak detection, which
enforces a down change after each up change, and vice versa.

Isotonic regression is another example of a constrained changepoint detection
model. There is no limit on the number of segments $K$, but the segment means
are constrained to be non-decreasing. This problem can be solved in $O(n)$
time using the pool-adjacent-violators algorithm (Mair et al., 2009), or in
$O(n\log n)$ time using a dynamic programming algorithm (Rote, unpublished).
An L1 relaxation of this problem is known as nearly-isotonic regression
(Tibshirani et al., 2011). A problem known as reduced isotonic regression
occurs by imposing an additional constraint of $K$ segments (Schell and Singh,
1997). The techniques for solving this problem lead to sub-quadratic time
algorithms (Hardwick and Stout, 2014), but do not generalize to other kinds of
constraints (such as unimodal regression or peak detection).

Our contribution in this paper is proving that the functional pruning
technique can be generalized to constrained changepoint models (Table 1). Our
resulting Generalized Pruned Dynamic Programming Algorithm (GPDPA) enjoys
$O(Kn\log n)$ time complexity, and works for any changepoint model with affine
constraints between adjacent segment means (including isotonic regression,
unimodal regression, and peak detection).

##  3 Isotonic regression and changepoint models

Although our proposed algorithm can solve many constrained changepoint
detection problems (Section 4.5), we will simplify our discussion by
emphasizing the isotonic regression model.

###  3.1 Classical isotonic regression

The classical isotonic regression model is defined as the most likely sequence
of non-decreasing segment means. More precisely, assume that the data
$\mathbf{y}\in\mathbb{R}^{n}$ are a realization of a probability distribution
with mean parameter $\mathbf{m}\in\mathbb{R}^{n}$. For example, assuming
$y_{t}\sim\mathcal{N}(m_{t},\sigma^{2})$ and performing maximum likelihood
inference results in a convex minimization problem with affine constraints,

$\displaystyle\operatorname*{minimize}_{\mathbf{m}\in\mathbb{R}^{n}}$
$\displaystyle\ \ \sum_{t=1}^{n}\ell(y_{t},m_{t})$ (1) subject to
$\displaystyle\ \ m_{t}\leq m_{t+1},\,\forall t<n.$

The convex loss function
$\ell:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ in the case of the
Gaussian likelihood is the square loss $\ell(y,m)=(y-m)^{2}$. This
optimization problem (1) is referred to as isotonic regression, and can be
efficiently solved in $O(n)$ time using the Pool-Adjacent-Violators Algorithm
(PAVA) (Best and Chakravarti, 1990).

Since isotonic regression imposes no limit on the number of changepoints
($m_{t}<m_{t+1}$), it tends to overfit. For example, consider the toy data set
$\mathbf{y}=\left[\begin{array}[]{cccccc}2&5&30&34&600&621\end{array}\right]\in\mathbb{R}^{6}$.
Because these data are strictly increasing, the isotonic regression (1)
solution is the trivial model $m_{t}=y_{t}$. However, these data contain only
two large changes. To recover these changes, we could instead use the segment
neighborhood model, which we discuss in the next section.

###  3.2 Segment neighborhood changepoint model

The segment neighborhood model of Auger and Lawrence (1989) uses the same cost
function as isotonic regression, but a different constraint set. There is no
constraint on the direction of changes, but there must be exactly $K\leq n$
distinct segments ($K-1$ changes).

$\displaystyle\operatorname*{minimize}_{\mathbf{m}\in\mathbb{R}^{n}}$
$\displaystyle\ \ \sum_{t=1}^{n}\ell(y_{t},m_{t})$ (2) subject to
$\displaystyle\ \ \sum_{t=1}^{n-1}I(m_{t}\neq m_{t+1})=K-1.$

This optimization problem is non-convex since the model complexity is the
number of changepoints, measured via the non-convex indicator function $I$.
Nonetheless, the optimal solution can be computed in $O(Kn^{2})$ time using
the standard dynamic programming algorithm (Auger and Lawrence, 1989). By
exploiting the structure of the convex loss function $\ell$, the pruned
dynamic programming algorithm of Rigaill (2010) computes the same optimal
solution in faster $O(Kn\log n)$ time.

Unlike isotonic regression, the segment neighborhood model does not constrain
the direction of the changes. Thus, for some data sets $\mathbf{y}$, the
segment neighborhood model may recover a change down ($m_{t}>m_{t+1}$). For
applications where isotonic regression is used, it would be desirable to
compute a model with $K$ non-decreasing segment means. This results in the
reduced isotonic regression problem, which we introduce in the next section.

###  3.3 Reduced isotonic regression

The idea of fitting a non-decreasing function with a limited number of
changepoints has been previously described as reduced isotonic regression
(Schell and Singh, 1997). Combining the constraints of the isotonic regression
(1) and segment neighborhood (2) problems gives

$\displaystyle\operatorname*{minimize}_{\mathbf{m}\in\mathbb{R}^{n}}$
$\displaystyle\ \ \sum_{t=1}^{n}\ell(y_{t},m_{t})$ (3) subject to
$\displaystyle\ \ \sum_{t=1}^{n-1}I(m_{t}\neq m_{t+1})=K-1,$ $\displaystyle\ \
m_{t}\leq m_{t+1},\,\forall t<n.$

In the next section, we explain how functional pruning can be used for solving
this and related changepoint problems.

##  4 Functional pruning algorithms for constrained changepoint models

We begin by discussing an algorithm for solving the reduced isotonic
regression problem, then explain how the algorithm generalizes to other
constrained changepoint problems.

###  4.1 Equivalent optimization space

The reduced isotonic regression problem (3) has $n$ segment mean variables
$m_{t}$, one for each data point $t$. To derive our algorithm, we re-write the
problem in terms of the mean $u_{k}\in\mathbb{R}$ and endpoint
$t_{k}\in\\{1,\dots,n\\}$ for each segment $k\in\\{1,\dots,K\\}$.

######  Definition 1 (Reduced isotonic regression optimization space).

Let $(\mathbf{u},\mathbf{t})\in{\mathcal{I}}^{n}_{K}$ be the set of non-
decreasing segment means $u_{1}\leq\cdots\leq u_{K}$ and increasing
changepoint indices $0=t_{0}<t_{1}<\cdots<t_{K-1}<t_{K}=n$.

Each segment mean $u_{k}$ is assigned to data points
$\tau\in(t_{k-1},t_{k}]\subset\\{1,\dots,n\\}$, resulting in the following
cost for each segment $k\in\\{1,\dots,K\\}$,

$h_{t_{k-1},t_{k}}(u_{k})=\sum_{\tau=t_{k-1}+1}^{t_{k}}\ell(y_{\tau},u_{k}).$
(4)

The reduced isotonic regression problem can be equivalently written as

$\operatorname*{minimize}_{(\mathbf{u},\mathbf{t})\in{\mathcal{I}}^{n}_{K}}\sum_{k=1}^{K}h_{t_{k-1},t_{k}}(u_{k})$
(5)

Rather than explicitly summing over data points $i$ as in problem (3), this
problem uses the equivalent sum over segments $k$.

$C_{1,1}(\mu)=(\mu-2)^{2}$$\min_{\mu}C_{1,1}(\mu)$$C^{\leq}_{1,1}(\mu)=$$\min_{x\leq\mu}C_{1,1}(x)$0123401234segment
mean $\mu$costMin-less computation for data point 1

constrained$C_{2,2}(\mu)=$unconstrained$(\mu-1)^{2}$$C^{\leq}_{1,1}(\mu)+(\mu-1)^{2}$0.02.55.07.501234segment
mean $\mu$costCost of 2 segments up to data point 2

Figure 1: Comparison of previous unconstrained algorithm (grey) with new
algorithm that constrains segment means to be non-decreasing (red), for the
toy data set $\mathbf{y}=[2,1,0,4]\in\mathbb{R}^{4}$ and the square loss.
Left: rather than computing the unconstrained minimum (constant grey
function), the new algorithm computes the min-less operator (red), resulting
in a larger cost when the segment mean is less than the first data point
($\mu<2$). Right: adding the cost of the second data point $(\mu-1)^{2}$ and
minimizing yields equal means $u_{1}=u_{2}=1.5$ for the constrained model and
decreasing means $u_{1}=2,\,u_{2}=1$ for the unconstrained model.

###  4.2 Dynamic programming update rules

Optimization problem (5) has $K$ segment mean variables $u_{k}$ and $K-1$
changepoint index variables $t_{k}$. Minimizing over all variables except the
last segment mean $u_{K}$ results in the following definition of the optimal
cost.

######  Definition 2 (Optimal cost with last segment mean $\mu$).

Let $C_{K,n}(\mu)$ be the optimal cost of the segmentation with $K$ segments,
up to data point $n$, with last segment mean $\mu$:

$C_{K,n}(\mu)=\min_{(\mathbf{u},\mathbf{t})\in{\mathcal{I}}^{n}_{K}\ |\
u_{K}=\mu}\ \left\\{\sum_{k=1}^{K}h_{t_{k-1},t_{k}}(u_{k})\right\\}.$ (6)

As in the PDPA of Rigaill (2010), our proposed dynamic programming algorithm
uses an exact representation of the $C_{k,t}:\mathbb{R}\rightarrow\mathbb{R}$
cost functions. Each $C_{k,t}(\mu)$ is represented as a piecewise function on
intervals of $\mu$. This is implemented as a linked list of FunctionPiece
objects in C++ (for details see Section B). Each element of the linked list
represents a convex function piece, and implementation details depend on the
choice of the loss function $\ell$ (for an example using the square loss see
Section 4.3).

In the original unconstrained PDPA, computing the $C_{k,t}(u_{k})$ function
requires taking the minimum of $C_{k,t-1}(u_{k})$ (a function of the last
segment mean $u_{k}$) and
$\hat{C}_{k-1,t-1}=\min_{u_{k-1}}C_{k-1,t-1}(u_{k-1})$ (the constant loss
resulting from an unconstrained minimization with respect to the previous
segment mean $u_{k-1}$). The main novelty of our paper is the discovery that
this update can also be computed efficiently for constrained problems. For
example in reduced isotonic regression the second term is no longer a
constant, but instead a function of $u_{k}$,
$C_{k-1,t-1}^{\leq}(u_{k})=\min_{u_{k-1}\leq u_{k}}C_{k-1,t-1}(u_{k-1})$,
which we refer to as the min-less operator (Figure 1, left).

######  Definition 3 (Min-less operator).

Given any real-valued function $f:\mathbb{R}\rightarrow\mathbb{R}$, we define
the min-less operator of that function as $f^{\leq}(\mu)=\min_{x\leq\mu}f(x)$.

The min-less operator is used in the following Theorem, which states the
update rules used in our proposed algorithm.

######  Theorem 1 (Generalized Pruned Dynamic Programming Algorithm for
reduced isotonic regression).

The optimal cost functions $C_{k,t}$ can be recursively computed using the
following update rules.

  1. 1.

For $k=1$ we have $C_{1,1}(\mu)=\ell(y_{1},\mu)$, and for the other data
points $t>1$ we have

$C_{1,t}(\mu)=C_{1,t-1}(\mu)+\ell(y_{t},\mu)$ (7)

  2. 2.

For $k>1$ and $t=k$ we have

$C_{k,k}(\mu)=\ell(y_{k},\mu)+C_{k-1,k-1}^{\leq}(\mu)$ (8)

  3. 3.

In all other cases we have

$C_{k,t}(\mu)=\ell(y_{t},\mu)+\min\\{C_{k-1,t-1}^{\leq}(\mu),\,C_{k,t-1}(\mu)\\}.$
(9)

###### Proof.

Case 1 and 2 follow from Definition 2, and there is a proof for case 3 in
Section A.

∎

The dynamic programming algorithm requires computing $O(Kn)$ cost functions
$C_{k,t}$. As in the original pruned dynamic programming algorithm, the time
complexity of the algorithm is $O(KnI)$ where $I$ is the number of intervals
(convex function pieces; candidate changepoints) that are used to represent
the cost functions. The theoretical maximum number of intervals is $I=O(n)$,
implying a time complexity of $O(Kn^{2})$ (Rigaill, 2015). However, this
maximum is only achieved in pathological synthetic data sets, such as a
monotonic increasing data sequence. The average number of intervals in real
data sets is empirically $I=O(\log n)$, as we will show in Section 5.1. Thus
the average time complexity of the algorithm is $O(Kn\log n)$.

###  4.3 Example and comparison with unconstrained case

To clarify the discussion, consider the toy data set
$\mathbf{y}=\left[\begin{array}[]{cccccc}2&1&0&4\end{array}\right]\in\mathbb{R}^{4}$
and the square loss $\ell(y,\mu)=(y-\mu)^{2}$. The first step of the algorithm
is to compute the minimum and the maximum of the data (0,4) in order to bound
the possible values of the segment mean $\mu$. Then the algorithm computes the
optimal cost in $k=1$ segment up to data point $t=1$:

$C_{1,1}(\mu)=(2-\mu)^{2}=4-4\mu+\mu^{2}\text{ (for $\mu\in[0,4]$)}$ (10)

This function can be stored for all values of $\mu$ via the three real-valued
coefficients ($\text{constant}=4$, $\text{linear}=-4$, $\text{quadratic}=1$).
To compute the optimal cost in $K=2$ segments, we first compute the min-less
operator (red curve on left of Figure 1),

$C_{1,1}^{\leq}(\mu)=\begin{cases}4-4\mu+\mu^{2}&\text{ if
}\mu\in[0,2],\,\mu^{\prime}=\mu,\\\ 0+0\mu+0\mu^{2}&\text{ if
}\mu\in[2,4],\,\mu^{\prime}=2.\end{cases}$ (11)

This function can be stored as a list of two intervals of $\mu$ values, each
with associated real-valued coefficients. In addition, to facilitate recovery
of the optimal parameters, we store the previous segment mean $\mu^{\prime}$
and endpoint (not shown). Note that $\mu^{\prime}=\mu$ means that the equality
constraint is active ($u_{1}=u_{2}$).

By adding the first min-less function $C_{1,1}^{\leq}(\mu)$ to the cost of the
second data point $(\mu-2)^{2}$ we obtain the optimal cost in $K=2$ segments
up to data point $t=2$,

$C_{2,2}(\mu)=\begin{cases}5-6\mu+2\mu^{2}&\text{ if
}\mu\in[0,2],\,\mu^{\prime}=\mu,\\\ 1-2\mu+1\mu^{2}&\text{ if
}\mu\in[2,4],\,\mu^{\prime}=2.\end{cases}$ (12)

Note that the minimum of this function is achieved at $\mu=1.5$ which occurs
in the first of the two function pieces (red curve on right of Figure 1), with
an equality constraint active. This implies the optimal model up to data point
$t=2$ with $k=2$ non-decreasing segment means actually has no change
($u_{1}=u_{2}=1.5$). In contrast, the minimum of the cost computed by the
unconstrained algorithm is at $u_{2}=1$ (grey curve on right of Figure 1),
resulting in a change down from $u_{1}=2$.

###  4.4 The PeakSeg up-down constraint

The PeakSeg model described by Hocking et al. (2015) is the most likely
segmentation where the first change is up, all up changes are followed by down
changes, and all down changes are followed by up changes. More precisely, the
constrained optimization problem can be stated as

$\displaystyle\operatorname*{minimize}_{\begin{subarray}{c}\mathbf{u}\in\mathbb{R}^{K}\\\
0=t_{0}<t_{1}<\cdots<t_{K-1}<t_{K}=n\end{subarray}}$ $\displaystyle\ \
\sum_{k=1}^{K}h_{t_{k-1},t_{k}}(u_{k})$ (13) subject to $\displaystyle\ \
u_{k-1}\leq u_{k}\ \forall k\in\\{2,4,\dots\\},$ $\displaystyle\ \ u_{k-1}\geq
u_{k}\ \forall k\in\\{3,5,\dots\\}.$

pruning at data $t=34$cost up to data $t=35$pruning at data
$t=35$$C^{\geq}_{2,34}$$C_{3,34}$$M_{3,34}=\min\\{C_{3,34},C^{\geq}_{2,34}\\}$$C_{3,35}=\ell_{35}+$$M_{3,34}$$\ell_{35}$$M_{3,34}$$C_{3,35}$$C^{\geq}_{2,35}$$M_{3,35}=\min\\{C_{3,35},C^{\geq}_{2,35}\\}$0.00.10.20.30.40.000.250.500.750.000.250.500.750.000.250.500.75mean
$u_{3}$ of segment 3cost value $C(u_{3})$cost
typeaaaa$\ell_{t}(u_{3})=\ell(y_{t},u_{3})$cost of data
$t$$C^{\geq}_{2,t}=$cost ofnon-increasingchange after $t$min
$M_{3,t}$$C_{3,t}=$costup to data $t$ Figure 2:  Demonstration of GPDPA for
the PeakSeg model (13) with $k=3$ segments. Cost functions are stored as
piecewise functions on intervals (black dots show limits between function
pieces). Left: the min $M_{3,34}$ is the minimum of two functions:
$C^{\geq}_{2,34}$ is the cost if the second segment ends at data point $t=34$
(the min-more operator forces a non-increasing change after), and $C_{3,34}$
is the cost if the second segment ends before that. Middle: the cost
$C_{3,35}$ is the sum of the min $M_{3,34}$ and the cost of the next data
point $\ell_{35}$. Right: in the next step, all previously considered
changepoints are pruned (cost $C_{3,35}$), since the model with a the second
segment ending at data point $t=35$ is always less costly ($C^{\geq}_{2,35}$).

Our proposed Generalized Pruned Dynamic Programming Algorithm (GPDPA) can be
used to solve the PeakSeg problem. The initialization $k=1$ is the same as in
the reduced isotonic regression solver (Section 4.2). The dynamic programming
updates for even $k\in\\{2,4,\dots\\}$ are also the same. However, to
constrain non-increasing changes, the updates for odd $k\in\\{3,5,\dots\\}$
are

$C_{k,t}(\mu)=\ell(y_{t},\mu)+\min\\{C_{k-1,t-1}^{\geq}(\mu),\,C_{k,t-1}(\mu)\\},$
(14)

where the min-more operator is defined for any function
$f:\mathbb{R}\rightarrow\mathbb{R}$ as $f^{\geq}(\mu)=\min_{x\geq\mu}f(x)$.
Figure 2 shows the geometric interpretation of the min-more operator, along
with an example of how the $\min\\{\\}$ operation performs pruning. We
implemented this algorithm using the Poisson loss $\ell(y,\mu)=\mu-y\log\mu$,
since our application in Section 5 is on count data
$y\in\mathbb{Z}_{+}=\\{0,1,2,\dots\\}$. We implemented this algorithm in C++,
and our free/open-source code is available as the PeakSegPDPA function in the
coseg R package for constrained optimal segmentation
(https://github.com/tdhock/coseg). Implementation details can be found in
Section B.

###  4.5 General affine inequality constraints between adjacent segment means

In this section we briefly discuss how our proposed Generalized Pruned Dynamic
Programming Algorithm (GPDPA) can be used to solve any optimization problem
with affine inequality constraints between adjacent segment means. For each
change $k\in\\{1,\dots,K-1\\}$, let $a_{k},b_{k},c_{k}\in\mathbb{R}$ be
arbitrary coefficients that define affine functions
$g_{k}(u_{k},u_{k+1})=a_{k}u_{k}+b_{k}u_{k+1}+c_{k}$. The changepoint
detection problem with general affine constraints is

$\displaystyle\operatorname*{minimize}_{\begin{subarray}{c}\mathbf{u}\in\mathbb{R}^{K}\\\
0=t_{0}<t_{1}<\cdots<t_{K-1}<t_{K}=n\end{subarray}}$ $\displaystyle\ \
\sum_{k=1}^{K}h_{t_{k-1},t_{k}}(u_{k})$ (15) subject to $\displaystyle\ \
\forall k\in\\{1,\dots,K-1\\},$ $\displaystyle\ \ g_{k}(u_{k},u_{k+1})\leq 0.$

Some examples of models that are special cases:

  1. 1.

If we take all $a_{k},b_{k},c_{k}=0$ then the constraints are trivially
satisfied, we recover the unconstrained segment neighborhood problem (2).

  2. 2.

If we take all $a_{k}=1$, $b_{k}=-1$ and $c_{k}=0$ we recover the reduced
isotonic regression problem (5).

  3. 3.

For the PeakSeg problem (13), we take all $c_{k}=0$. For odd
$k\in\\{1,3,\dots\\}$ we take $a_{k}=1$, $b_{k}=-1$ and for even
$k\in\\{2,4,\dots\\}$ we take $a_{k}=-1$, $b_{k}=1$.

To solve these problems, we need to compute the analog of the min-less/more
operator, which we call the constrained minimization operator. For any cost
function $f:\mathbb{R}\rightarrow\mathbb{R}$ and constraint function
$g:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$, we define the constrained
minimization operator $f^{g}:\mathbb{R}\rightarrow\mathbb{R}$ as

$f^{g}(u_{k})=\min_{u_{k-1}:g(u_{k-1},u_{k})\leq 0}f(u_{k-1}).$ (16)

When $g$ is affine, the constrained minimization operator is either non-
decreasing or non-increasing. In this case it can be computed using a simple
algorithm that scans the piecewise function $f$ either from left to right or
right to left. When a local minimum is found, its value is recorded, and a
constant function piece is added (for details see pseudocode for MinLess
algorithm in Section B.2). The constrained minimization operator is used in
the following general dynamic programming update rule which can be used to
compute the solution to (15)

$C_{k,t}(\mu)=\ell(y_{t},\mu)+\min\\{C_{k,t-1}(\mu),\,C_{k-1,t-1}^{g_{k-1}}(\mu)\\}.$
(17)

We note that this update rule is valid for constraint functions $g$ more
general than affine functions. However, the closed-form computation of the
constrained minimization operator (16) would possibly be much more difficult
for these more general constraint functions (e.g. quadratic constraint
functions).

##  5 Results on peak detection in ChIP-seq data

The real data analysis problem that motivates this work is the detection of
peaks in ChIP-seq data (Bailey et al., 2013), which are typically represented
as a vector of non-negative counts $\mathbf{y}\in\mathbb{Z}_{+}^{n}$ of
aligned sequence reads for $n$ continguous bases in a genome. Data sizes are
between $n=10^{5}$ (maximum of the benchmark we consider) and $n=10^{8}$
(largest region with no gaps in the human genome hg19). A peak detector can be
represented as a function $c(\mathbf{y})\in\\{0,1\\}^{n}$ for binary
classification at every base position. The positive class is peaks (genomic
regions with large values, representing protein binding or modification) and
the negative class is background noise (small values).

In the supervised learning framework of Hocking et al. (2016), a data set
consists of $m$ count data vectors $\mathbf{y}_{1},\dots,\mathbf{y}_{m}$ along
with labels $L_{1},\dots,L_{m}$ that identify regions with and without peaks.
Briefly, the number of errors $E[c(\mathbf{y}_{i}),L_{i}]$ is the total of
false positives (negative labels with a predicted peak) plus false negatives
(positive labels with no predicted peak). The benchmark consists of seven
histone ChIP-seq data sets, each with a different peak pattern (experiment
type, labeler, cell types). The goal in each data set is to learn the pattern
encoded in the labels, and find a classifier $c$ that minimizes the total
number of incorrectly predicted labels in a held-out test set:

$\operatorname*{minimize}_{c}\sum_{i=1}^{m}E\left[c(\mathbf{y}_{i}),L_{i}\right].$
(18)

Hocking et al. (2015) proposed a constrained dynamic programming algorithm
(CDPA) to approximately compute the optimal changepoints, subject to the
PeakSeg up-down constraint (Section 4.4). The CDPA has been shown to achieve
state-of-the-art peak detection accuracy, by classifying even-numbered
segments $k$ as peaks, and odd-numbered segments $k$ as background noise.
However, its quadratic $O(Kn^{2})$ time complexity makes it too slow to run on
large ChIP-seq data sets.

maxmedian, inter-quartile range1010025310001000087263169$n$ = data points to
segment (log scale)$I$ = intervals stored(log scale)

1 second1 minute1 hourGPDPA$O(n\log n)$PDPA$O(n\log
n)$CDPA$O(n^{2})$11001000010001000087263169$n$ = data points to segment (log
scale)seconds (log scale)

Figure 3: Empirical speed analysis on 2752 count data vectors from the histone
mark ChIP-seq benchmark. For each vector we ran the GPDPA with the up-down
constraint and a max of $K=19$ segments. The expected time complexity is
$O(KnI)$ where $I$ is the average number of intervals (function pieces;
candidate changepoints) stored in the $C_{k,t}$ cost functions. Left: number
of intervals stored is $I=O(\log n)$ (median, inter-quartile range, and
maximum over all data points $t$ and segments $k$). Right: time complexity of
the GPDPA is $O(n\log n)$ (median line and min/max band).

In this section, we show that our proposed GPDPA can be used to overcome this
speed drawback, while maintaining state-of-the-art accuracy. To show the
importance of enforcing the up-down constraint, we consider the unconstrained
Pruned Dynamic Programming Algorithm (PDPA) of Rigaill (2010) as a baseline
(Table 1). We also compare against two popular heuristics from the
bioinformatics literature, in order to demonstrate that constrained
optimization algorithms such as the CDPA and GPDPA are more accurate.

###  5.1 Empirical time complexity in ChIP-seq data

The ChIP-seq benchmark consists of seven labeled histone data sets. Overall
there are 2752 count data vectors $\mathbf{y}_{i}$ to segment, varying in size
from $n=87$ to $n=263169$ data. For each count data vector $\mathbf{y}_{i}$,
we ran each algorithm (CDPA, PDPA, GDPDA) with a maximum of $K=19$ segments.
This implies a maximum of 9 peaks (one for each even-numbered segment), which
is more than enough in these relatively small data sets. To analyze the
empirical time complexity, we recorded the number of intervals stored in the
$C_{k,t}$ cost functions (Section 4), as well as the computation time in
seconds.

As in the PDPA, the time complexity of our proposed GPDPA is $O(KnI)$, which
depends on the number of intervals $I$ (candidate changepoints) stored in the
$C_{k,t}$ cost functions (Rigaill, 2015). We observed that the number of
intervals stored by the GPDPA increases as a sub-linear function of the number
of data points $n$ (left of Figure 3). For the largest data set ($n=263169$),
the algorithm only stored median=16 and maximum=43 intervals. The most
intervals stored was 253 for one data set with $n=7776$. These results suggest
that our proposed GPDPA only stores on average $O(\log n)$ intervals (possible
changepoints), as in the original PDPA. The overall empirical time complexity
is thus $O(Kn\log n)$ for $K$ segments and $n$ data points.

We recorded the timings of each algorithm for computing models with up to
$K=19$ segments (a total of 10 peak models $k\in\\{1,3,\dots,19\\}$, from 0 to
9 peaks). Since $K$ is constant, the expected time complexity was $O(n^{2})$
for the CDPA and $O(n\log n)$ for the PDPA and GPDPA. In agreement with these
expectations, our proposed GPDPA shows $O(n\log n)$ asymptotic timings similar
to the PDPA (right of Figure 3).

It is clear that the $O(n^{2})$ CDPA algorithm is slower than the other two
algorithms, especially for larger data sets. For the largest count data vector
($n=263169$), the CDPA took over two hours, but the GPDPA took only about two
minutes. Our proposed GPDPA is nearly as fast as MACS (Zhang et al., 2008), a
heuristic from the bioinformatics literature which took about 1 minute to
compute 10 peak models for this data set.

The total computation time to process all 2752 count data vectors was 156
hours for the CDPA, and only 6 hours for the GPDPA (26 times faster). Overall,
these results suggest that our proposed GPDPA enjoys $O(n\log n)$ time
complexity in ChIP-seq data, which makes it possible to use for very large
data sets.

###  5.2 Test accuracy in ChIP-seq data

Figure 4: Four-fold cross-validation was used to estimate peak detection
accuracy. Each panel shows one of seven ChIP-seq data sets, labeled by
experiment (Broad H3K36me3), labeler (AM), and cell types (immune). Each black
circle shows test AUC in one of four cross-validation folds, the shaded grey
circle is the mean, and the vertical line is the maximum mean in each data
set. It is clear that the proposed GPDPA is just as accurate as the previous
state-of-the-art CDPA, and both are more accurate than the other baseline
methods.

For the optimal changepoint detection algorithms (CDPA, PDPA, GPDPA), the
prediction problem simplifies to selecting the number of segments
$K_{i}\in\\{1,3,\dots,19\\}$ for each data vector $i$, resulting in a
predicted peak vector $c^{K_{i}}(\mathbf{y}_{i})\in\\{0,1\\}^{n}$. We select
the number of segments using an oracle penalty
$K_{i}^{\lambda}=\operatorname*{arg\,min}_{k}l_{ik}+\lambda o_{ik}$ (Cleynen
and Lebarbier, 2014), where $l_{ik}$ is the Poisson loss and $o_{ik}$ is the
oracle model complexity for the model with $k$ segments for data vector $i$.
The problem thus simplifies to learning a scalar penalty constant $\lambda$,

$\operatorname*{minimize}_{\lambda}\sum_{i=1}^{m}E\left[c^{K_{i}^{\lambda}}(\mathbf{y}_{i}),L_{i}\right].$
(19)

To demonstrate that changepoint detection algorithms are more accurate than
typical heuristics from the bioinformatics literature, we also compared with
the MACS and HMCanBroad methods (Zhang et al., 2008; Ashoor et al., 2013).
MACS is a popular heuristic for data with a sharp peak pattern such as
H3K4me3, and HMCanBroad is a popular heuristic for data with a broad peak
pattern such as H3K36me3. Although they are not designed for supervised
learning, we trained them by performing grid search over a single significance
threshold parameter (qvalue for MACS and finalThreshold for HMCanBroad).

In each of the seven data sets in the histone benchmark, we performed four-
fold cross-validation and computed test AUC (area under the Receiver Operating
Characteristic curve) to estimate the accuracy of each algorithm. The previous
algorithm with state-of-the-art accuracy on this benchmark was the CDPA, which
enforces the up-down constraint on segment means. We expected our proposed
GPDPA to perform just as well, since it also enforces that constraint. In
agreement with our expectation, we observed that the CDPA and GPDPA yield
comparable test AUC in all seven data sets (Figure 4). In contrast, the
unconstrained PDPA had much lower test AUC in several data sets, because of
lower true positive rates. These results provide convincing evidence that the
constraint is necessary for optimal peak detection accuracy.

Since the baseline HMCanBroad algorithm was designed for data with a broad
peak pattern, we expected it to perform well in the H3K36me3 data. In
agreement with this expectation, HMCanBroad showed state-of-the-art test AUC
in two H3K36me3 data sets (broad peak pattern), but was very inaccurate in
four H3K4me3 data sets (sharp peak pattern). We expected the baseline MACS
algorithm to perform well in the H3K4me3 data sets, since it was designed for
data with a sharp peak pattern. In contrast to this expectation, MACS had test
AUC values much lower than the optimization-based algorithms in all seven data
sets (Figure 4). These results suggest that constrained optimal changepoint
detection algorithms are more accurate than the heuristics from the
bioinformatics literature.

##  6 Discussion and conclusions

Algorithms for changepoint detection can be classified in terms of time
complexity, optimality, constraints, and pruning techniques (Table 1). In this
paper, we investigated generalizing the functional pruning technique
originally discovered by Rigaill (2010) and Johnson (2011). We showed that the
functional pruning technique can be used to compute optimal changepoints
subject to affine constraints on adjacent segment mean parameters.

We showed that our proposed Generalized Pruned Dynamic Programming Algorithm
(GPDPA) enjoys the same log-linear $O(Kn\log n)$ time complexity as the
original unconstrained PDPA, when applied to peak detection in ChIP-seq data
sets (Figure 3). However, we observed that the up-down constrained GPDPA is
much more accurate than the unconstrained PDPA (Figure 4). These results
suggest that the up-down constraint is necessary for computing a changepoint
model with optimal peak detection accuracy. Indeed, we observed that the GPDPA
enjoys the same state-of-the-art accuracy as the previous best, the relatively
slow quadratic $O(Kn^{2})$ time CDPA.

We observed that the heuristic algorithms which are popular in the
bioinformatics literature (MACS, HMCanBroad) are much less accurate than the
optimal changepoint detection algorithms (CDPA, PDPA, GPDPA). In the past
these sub-optimal heuristics have been preferred because of their speed. For
example, the CDPA took 2 hours to compute 10 peak models in the largest data
set in the ChIP-seq benchmark, whereas the GPDPA took 2 minutes, and the MACS
heuristic took 1 minute. Using our proposed GPDPA, it is now possible to
compute highly accurate models in an amount of time that is comparable to
heuristic algorithms. Our proposed GPDPA can now be used as an optimal
alternative to heuristic algorithms, even for large data sets.

For future work we will be interested in exploring pruning techniques for
other constrained changepoint models. When the number of expected changepoints
grows with the number of data points, then $K=O(n)$ and our proposed GPDPA has
$O(n^{2}\log n)$ average time complexity (since it computes all models with
$1,\dots,K$ segments). We have already started modifying the GPDPA for optimal
partitioning (Jackson et al., 2005), which results in the Generalized
Functional Prunining Optimal Partitioning (GFPOP) algorithm (Section B.5). It
computes the $K$-segment model for a single penalty constant $\lambda$
(without computing models with $1,\dots,K-1$ segments) in $O(n\log n)$ time.

##  7 Reproducible Research Statement

The source code and data used to create this manuscript (including all
figures) is available at https://github.com/tdhock/PeakSegFPOP-paper

##  8 Acknowledgements

This work was supported by a Discovery Frontiers project grant, “The Cancer
Genome Collaboratory,” jointly sponsored by the Natural Sciences and
Engineering Research Council (NSERC), Genome Canada (GC), the Canadian
Institutes of Health Research (CIHR) and the Canada Foundation for Innovation
(CFI).

The supplementary materials begin on this page.

##  Appendix A Proof of optimality of dynamic programming algorithm

In this section we give a proof of Theorem 1.

###### Proof.

Case 1 and 2 follow from the definition of $C_{K,t}(u)$.

We now focus on case 3. First notice that by definition of $C_{K,t+1}(u)$
(i.e. the optimal segmentation) we must have $C_{K,t+1}(u)\leq
C_{K,t}(u)+\ell(y_{t},u)$ and also $C_{K,t+1}(u)\leq
C_{K-1,t}(u)+\ell(y_{t},u)$. Thus we have
$C_{K,t+1}(u)\leq\min\\{C_{K,t}(u),C_{K-1,t}(u)\\}+\ell(y_{t+1},u)$.

Now let us assume,

$C_{K,t+1}(u)<\min\\{C_{K,t}(u),C_{K-1,t}(u)\\}+\ell(y_{t+1},u).$

We will show that this lead to a contradiction.

We consider the optimal segmentation
$(\mathbf{u},\mathbf{t})\in{\mathcal{I}}_{t+1}^{K}$ which achieves the optimum
of $C_{K,t+1}(u)$. We consider two possible cases:

Scenario 1: $t_{K}<t$.

Define $\mathbf{t}^{\prime}$ such that for all $i<K$, we have
$t^{\prime}_{i}=t_{i}$ and $t^{\prime}_{K}=t$. We have
$(\mathbf{u},\mathbf{t}^{\prime})\in{\mathcal{I}}_{t}^{K}$. We can thus
decompose $C_{K,t+1}(u)$ as

$C_{K,t+1}(u)=\sum_{k=1}^{K}h_{t^{\prime}_{k-1},t^{\prime}_{k}}(u_{k})+\ell(y_{t+1},u).$

By assumption we would recover
$\sum_{k=1}^{K}h_{t^{\prime}_{k-1},t^{\prime}_{k}}(u_{k})<C_{K,t}(u)$ which is
a contradiction by definition of $C_{K,t}(u)$.

Scenario 2: $t_{K}=t$.

Define $\mathbf{t}^{\prime}$ such that for all $i<K-1$, we have
$t^{\prime}_{i}=t_{i}$ and $t^{\prime}_{K-1}=t$. Also define
$\mathbf{u}^{\prime}$ such that for all $k\leq K-1$, we have
$u^{\prime}_{k}=u_{k}$. Thus
$(\mathbf{u}^{\prime},\mathbf{t}^{\prime})\in{\mathcal{I}}_{t}^{K-1}$, and can
then decompose $C_{K,t+1}(u)$ as

$C_{K,t+1}(u)=\sum_{k=1}^{K}h_{t^{\prime}_{k-1},t^{\prime}_{k}}(u^{\prime}_{k})+\ell(y_{t+1},u).$

By assumption we would recover
$\sum_{k=1}^{K-1}h_{t^{\prime}_{k-1},t^{\prime}_{k}}(u^{\prime}_{k})<C_{K-1,t}(u)$
which is a contradiction by definition of $C_{K-1,t}(u)$.

∎

We have thus proved that the dynamic programming update rules can be used for
computing the optimal cost functions $C_{k,t}$.

##  Appendix B Algorithm pseudocode

In this section we give pseudocode for our proposed Generalized Pruned Dynamic
Programming Algorithm (GPDPA), and related algorithms.

###  B.1 GPDPA for reduced isotonic regression

We begin by providing a pseudocode solver for the simplest case, the reduced
isotonic regression problem. We propose the following data structures and sub-
routines for the computation:

  * •

FunctionPiece: a data structure which represents one piece of a $C_{k,t}(u)$
cost function (for one interval of mean values $u$). It has coefficients which
depend on the convex loss function $\ell$ (for the square loss it has three
real-valued coefficients $a,b,c$ which define a function $au^{2}+bu+c$). It
also has two real-valued elements for min/max mean values
$[\underline{u},\overline{u}]$ of this interval, meaning the function
$C_{k,t}(u)=au^{2}+bu+c$ for all $u\in[\underline{u},\overline{u}]$. Finally
it stores a previous segment endpoint $t^{\prime}$ (integer) and mean
$u^{\prime}$ (real).

  * •

FunctionPieceList: an ordered list of FunctionPiece objects, which exactly
stores a cost function $C_{k,t}(u)$ for all values of last segment mean $u$.

  * •

$\text{OnePiece}(y,\underline{u},\overline{u})$: a sub-routine that
initializes a FunctionPieceList with just one FunctionPiece $\ell(y,u)$
defined on $[\underline{u},\overline{u}]$.

  * •

$\text{MinLess}(t,f)$: an algorithm that inputs a changepoint and a
FunctionPieceList, and outputs the corresponding min-less operator $f^{\leq}$
(another FunctionPieceList), with the previous changepoint set to
$t^{\prime}=t$ for each of its pieces. This algorithm also needs to store the
previous mean value $u^{\prime}$ for each of the function pieces (see
pseudocode below).

  * •

$\text{MinOfTwo}(f_{1},f_{2})$: an algorithm that inputs two FunctionPieceList
objects, and outputs another FunctionPieceList object which is their minimum.

  * •

$\text{ArgMin}(f)$: an algorithm that inputs a FunctionPieceList and outputs
three values: the optimal mean $u^{*}=\operatorname*{arg\,min}_{u}f(u)$, the
previous segment end $t^{\prime}$ and mean $u^{\prime}$.

  * •

$\text{FindMean}(u,f)$ an algorithm that inputs a mean value and a
FunctionPieceList. It finds the FunctionPiece in $f$ with mean
$u\in[\underline{u},\overline{u}]$ contained in its interval, then outputs the
previous segment end $t^{\prime}$ and mean $u^{\prime}$ stored in that
FunctionPiece.

The above data structures and sub-routines are used in the following
pseudocode, which describes the GPDPA for solving the reduced isotonic
regression problem.

1: Input: data set $\mathbf{y}\in\mathbb{R}^{n}$, maximum number of segments
$K\in\\{2,\dots,n\\}$.

2: Output: matrices of optimal segment means $U\in\mathbb{R}^{K\times K}$ and
ends $T\in\\{1,\dots,n\\}^{K\times K}$

3: Compute min $\underline{y}$ and max $\overline{y}$ of $\mathbf{y}$.

4: $C_{1,1}\leftarrow\text{OnePiece}(y_{1},\underline{y},\overline{y})$

5: for data points $t$ from 2 to $n$:

6:
$C_{1,t}\leftarrow\text{OnePiece}(y_{t},\underline{y},\overline{y})+C_{1,t-1}$

7: for segments $k$ from 2 to $K$: for data points $t$ from $k$ to $n$: //
dynamic programming

8: $\text{min\\_prev}\leftarrow\text{MinLess}(t-1,C_{k-1,t-1})$ // this is
$C_{k-1,t-1}^{\leq}$

9: $\text{min\\_new}\leftarrow\text{min\\_prev}$ if $t=k$, else
$\text{MinOfTwo}(\text{min\\_prev},C_{k,t-1})$

10:
$C_{k,t}\leftarrow\text{min\\_new}+\text{OnePiece}(y_{t},\underline{y},\overline{y})$

11: for segments $k$ from 1 to $K$: // decoding for every model size $k$

12: $u^{*},t^{\prime},u^{\prime}\leftarrow\text{ArgMin}(C_{k,n})$

13: $U_{k,k}\leftarrow u^{*};\,T_{k,k}\leftarrow t^{\prime}$ // store mean of
segment $k$ and end of segment $k-1$

14: for segment $s$ from $k-1$ to $1$: // decoding for every segment $s<k$

15: if $u^{\prime}<\infty$: $u^{*}\leftarrow u^{\prime}$ // equality
constraint active, $u_{s}=u_{s+1}$

16: $t^{\prime},u^{\prime}\leftarrow\text{FindMean}(u^{*},C_{s,t^{\prime}})$

17: $U_{k,s}\leftarrow u^{*};\,T_{k,s}\leftarrow t^{\prime}$ // store mean of
segment $s$ and end of segment $s-1$

Algorithm 1 Generalized Pruned Dynamic Programming Algorithm (GPDPA) for
solving the reduced isotonic regression problem.

Algorithm 1 begins by computing the min/max on line 3. The main storage of the
algorithm is $C_{k,t}$, which should be initialized as a $K\times n$ array of
empty FunctionPieceList objects. The computation of $C_{1,t}$ for all $t$
occurs on lines 4–6.

The dynamic programming updates occur in the for loops on lines 7–10. Line 8
uses the MinLess sub-routine to compute the temporary FunctionPieceList
min_prev (which represents the function $C_{k-1,t-1}^{\leq}$). Line 9 sets the
temporary FunctionPieceList min_new to the cost of the only possible
changepoint if $t=k$; otherwise, it uses the MinOfTwo sub-routine to compute
the cost of the best changepoint for every possible mean value. Line 10 adds
the cost of data point $t$, and stores the resulting FunctionPieceList in
$C_{k,t}$.

The decoding of the optimal segment mean $U$ (a $K\times K$ array of real
numbers) and end $T$ (a $K\times K$ array of integers) variables occurs in the
for loops on lines 11–17. For a given model size $k$, the decoding begins on
line 12 by using the ArgMin sub-routine to solve
$u^{*}=\operatorname*{arg\,min}_{u}C_{k,n}(u)$ (the optimal values for the
previous segment end $t^{\prime}$ and mean $u^{\prime}$ are also returned).
Now we know that $u^{*}$ is the optimal mean of the last ($k$-th) segment,
which occurs from data point $t^{\prime}+1$ to $n$. These values are stored in
$U_{k,k}$ and $T_{k,k}$ (line 13). And we already know that the optimal mean
of segment $k-1$ is $u^{\prime}$. Note that the $u^{\prime}=\infty$ flag means
that the equality constraint is active (line 15). The decoding of the other
segments $s<k$ proceeds using the FindMean sub-routine (line 16). It takes the
cost $C_{s,t^{\prime}}$ of the best model in $s$ segments up to data point
$t^{\prime}$, finds the FunctionPiece that stores the cost of $u^{*}$, and
returns the new optimal values of the previous segment end $t^{\prime}$ and
mean $u^{\prime}$. The mean of segment $s$ is stored in $U_{k,s}$ and the end
of segment $s-1$ is stored in $T_{k,s}$ (line 17).

The time complexity of Algorithm 1 is $O(KnI)$ where $I$ is the complexity of
the MinLess and MinOfTwo sub-routines, which is linear in the number of
intervals (FunctionPiece objects) that are used to represent the cost
functions. There are pathological synthetic data sets for which the number of
intervals $I=O(n)$, implying a time complexity of $O(Kn^{2})$. However, the
average number of intervals in real data sets is empirically $I=O(\log n)$, so
the average time complexity of Algorithm 1 is $O(Kn\log n)$.

###  B.2 MinLess algorithm

The MinLess algorithm implements the min-less operator $f^{\leq}$ (Definition
3), which is an essential sub-routine of the GPDPA. The following sub-routines
are used to implement the MinLess algorithm.

  * •

$\text{GetCost}(p,u)$: an algorithm that takes a FunctionPiece object $p$, and
a mean value $u$, and computes the cost at $u$. For a square loss
FunctionPiece $p$ with coefficients $a,b,c\in\mathbb{R}$, we have
$\text{GetCost}(p,u)=au^{2}+bu+c$.

  * •

$\text{OptimalMean}(p)$: an algorithm that takes one FunctionPiece object, and
computes the optimal mean value. For a square loss FunctionPiece $p$ we have
$\text{OptimalMean}(p)=-b/(2a)$.

  * •

$\text{ComputeRoots}(p,d)$: an algorithm that takes one FunctionPiece object,
and computes the solutions to $p(u)=d$. For the square loss we propose to use
the quadratic formula. For other convex losses that do not have closed form
expressions for their roots, we propose to use Newton’s root finding method.
Note that for some constants $d$ there are no roots, and the algorithm needs
to report that.

  * •

$f.\text{push\\_piece}(\underline{u},\overline{u},p,u^{\prime})$: push a new
FunctionPiece at the end of FunctionPieceList $f$, with coefficients defined
by FunctionPiece $p$, on interval $[\underline{u},\overline{u}]$, with
previous segment mean set to $u^{\prime}$.

  * •

$\text{ConstPiece}(c)$: sub-routine that initializes a FunctionPiece $p$ with
constant cost $c$ (for the square loss it sets $a=b=0$ in $au^{2}+bu+c$).

1: Input: The previous segment end $t_{\text{prev}}$ (an integer), and
$f_{\text{in}}$ (a FunctionPieceList).

2: Output: FunctionPieceList $f_{\text{out}}$, initialized as an empty list.

3: $\text{prev\\_cost}\leftarrow\infty$

4: $\text{new\\_lower\\_limit}\leftarrow\text{LowerLimit}(f_{\text{in}}[0])$.

5: $i\leftarrow 0$; // start at FunctionPiece on the left

6: while $i<$ Length($f_{\text{in}}$): // continue until FunctionPiece on the
right

7: FunctionPiece $p\leftarrow f_{\text{in}}[i]$

8: if prev_cost = $\infty$: // look for min in this interval.

9: $\text{candidate\\_mean}\leftarrow\text{OptimalMean}(p)$

10: if $\text{LowerLimit}(p)<\text{candidate\\_mean}<\text{UpperLimit}(p)$:

11: $\text{new\\_upper\\_limit}\leftarrow\text{candidate\\_mean}$ // Minimum
found in this interval.

12: $\text{prev\\_cost}\leftarrow\text{GetCost}(p,\text{candidate\\_mean})$

13: $\text{prev\\_mean}\leftarrow\text{candidate\\_mean}$

14: else: // No minimum in this interval.

15: $\text{new\\_upper\\_limit}\leftarrow\text{UpperLimit}(p)$

16:
$f_{\text{out}}\text{.push\\_piece}(\text{new\\_lower\\_limit},\text{new\\_upper\\_limit},p,\infty)$

17: $\text{new\\_lower\\_limit}\leftarrow\text{new\\_upper\\_limit}$

18: $i\leftarrow i+1$

19: else: // look for equality of $p$ and prev_cost

20:
$(\text{small\\_root},\text{large\\_root})\leftarrow\text{ComputeRoots}(p,\text{prev\\_cost})$

21: if $\text{LowerLimit}(p)<\text{small\\_root}<\text{UpperLimit}(p)$:

22:
$f_{\text{out}}\text{.push\\_piece}(\text{new\\_lower\\_limit},\text{small\\_root},\text{ConstPiece}(\text{prev\\_cost}),\text{prev\\_mean})$

23: $\text{new\\_lower\\_limit}\leftarrow\text{small\\_root}$

24: $\text{prev\\_cost}\leftarrow\infty$

25: else: // no equality in this interval

26: $i\leftarrow i+1$ // continue to next FunctionPiece

27: if $\text{prev\\_cost}<\infty$: // ending on constant piece

28:
$f_{\text{out}}\text{.push\\_piece}(\text{new\\_lower\\_limit},\text{UpperLimit}(p),\text{ConstPiece}(\text{prev\\_cost}),\text{prev\\_mean})$

29: Set all previous segment end $t^{\prime}=t_{\text{prev}}$ for all
FunctionPieces in $f_{\text{out}}$

Algorithm 2 MinLess algorithm.

Consider Algorithm 2 which contains pseudocode for the computation of the min-
less operator. The algorithm initializes prev_cost (line 3), which is a state
variable that is used on line 8 to decide whether the algorithm should look
for a local minimum or an intersection with a finite cost. Since prev_cost is
initially set to $\infty$, the algorithm begins by following the convex
function pieces from left to right until finding a local minimum. If no
minimum is found in a given convex FunctionPiece (line 15), it is simply
pushed on to the end of the new FunctionPieceList (line 16). If a minimum
occurs within an interval (line 10), the cost and mean are stored (lines
11–12), and a new convex FunctionPiece is created with upper limit ending at
that mean value (line 16). Then the algorithm starts looking for another
FunctionPiece with the same cost, by computing the smaller root of the convex
loss function (line 20). When a FunctionPiece is found with a root in the
interval (line 21), a new constant FunctionPiece is pushed (line 22), and the
algorithm resumes searching for a minimum. At the end of the algorithm, a
constant FunctionPiece is pushed if necessary (line 28). The complexity of
this algorithm is $O(I)$ where $I$ is the number of FunctionPiece objects in
$f_{\text{in}}$.

The algorithm which implements the min-more operator is analogous. Rather than
searching from left to right, it searches from right to left. Rather than
using the small root (line 21), it uses the large root.

###  B.3 Implementation details

Some implementation details that we found to be important:

Weights

for data sequences that contain repeats it is computationally advantageous to
use a run-length encoding of the data, and a corresponding loss function. For
example if the data sequence 5,1,1,1,0,0,5,5 is encoded as $n=4$ counts
$y_{t}$ 5,1,0,5 with corresponding weights $w_{t}$ 1,3,2,2 then the Poisson
loss function for mean $\mu$ is $\ell(y_{t},w_{t},\mu)=w_{t}(\mu-
y_{t}\log\mu)$.

Mean cost

The text defines $C_{k,t}$ functions as the total cost. However for very large
data sets the cost values will be very large, resulting in numerical
instability. To overcome this issue we instead implemented update rules using
the mean cost. For weights $W_{t}=\sum_{i=1}^{t}w_{i}$, the update rule to
compute the mean cost is

$C_{k,t}(\mu)=\frac{1}{W_{t}}\left[\ell(y_{t},\mu)+W_{t-1}\min\\{C_{k,t-1}(\mu),\,C_{k-1,t-1}^{\leq}(\mu)\\}\right]$

Intervals in log(mean) space

For the Poisson model of non-negative count data $y_{t}\in\\{0,1,2,\dots\\}$
there is no possible mean $\mu$ value less than 0. We thus used $\log(\mu)$
values to implement intervals in FunctionPiece objects. For example rather
than storing $\mu\in[0,1]$ we store $\log\mu\in[-\infty,0]$.

Root finding

For the ComputeRoots sub-routine for the Poisson loss, we used Newton root
finding. For the larger root we solve $a\log\mu+b\mu+c=0$ (linear as
$\mu\rightarrow\infty$) and for the smaller root we solve $ax+be^{x}+c=0$
($x=\log\mu$, linear as $x\rightarrow-\infty$ and $\mu\rightarrow 0$). We stop
the root finding when the cost is near zero (absolute cost value less than
$10^{-12}$).

Storage

Since the dynamic programming update rule for $C_{k,t}$ only depends on
$C_{k-1,t-1}^{\geq}$ and $C_{k,t-1}$, these are the only functions that need
to be in memory, and the rest of the cost functions can be stored on disk
(until the decoding step). We used the Berkeley DB Standard Template Library
to store all the $C_{k,t}$ as a vector of FunctionPieceList objects.

###  B.4 Penalized version of reduced isotonic regression

Maidstone et al. (2016) proposed the Functional Pruning Optimal Partitioning
(FPOP) algorithm to solve the “penalized” or “optimal partitioning” version of
the segment neighborhood problem, where the constraint of $K$ segments is
replaced by a non-negative penalty $\lambda\in\mathbb{R}_{+}$ on the number of
changes in the objective function. Rather than computing all models from 1 to
$K$ segments (as in the PDPA), the FPOP algorithm computes the single model
with $K$ segments (without computing models from 1 to $K-1$ segments). The
same penalization idea can be applied to models with affine constraints
between adjacent segment means. The penalized version of the reduced isotonic
regression problem (3) can be stated as

$\displaystyle\operatorname*{minimize}_{\begin{subarray}{c}\mathbf{m}\in\mathbb{R}^{n}\\\
\mathbf{c}\in\\{0,1\\}^{n-1}\end{subarray}}$ $\displaystyle\ \
\sum_{t=1}^{n}\ell(y_{t},m_{t})+\lambda\sum_{t=1}^{n-1}I(c_{t}\neq 0)$ (20)
subject to $\displaystyle\ \ c_{t}=0\Rightarrow m_{t}=m_{t+1}$ $\displaystyle\
\ c_{t}=1\Rightarrow m_{t}\leq m_{t+1}.$

Note that the $c_{t}$ variable is a changepoint indicator. The same functional
pruning techniques used for the GPDPA can be exploited to create a solver for
this problem. This results in the Generalized Functional Pruning Optimal
Partitioning Algorithm (GFPOP, see Table 2).

| Segment Neighborhood | Optimal Partitioning
---|---|---
unconstrained | PDPA | FPOP
constrained | GPDPA | GFPOP
Table 2:  Algorithms for solving constrained and unconstrained versions of the
Segment Neighborhood and Optimal Partitioning problems. PDPA = Pruned Dynamic
Programming Algorithm, FPOP = Functional Pruning Optimal Partitioning, G =
Generalized (can handle affine constraints on adjacent segment means).

Let $\overline{C}_{\lambda,t}(u)$ be the penalized cost of the most likely
segmentation up to data point $t$, with last segment mean $u$. The
initialization for the first data point is
$\overline{C}_{\lambda,1}(u)=\ell(y_{1},u)$. The dynamic programming update
rule for all data points $t>1$ is

$\overline{C}_{\lambda,t}(u)=\ell(y_{t},u)+\min\\{\overline{C}_{\lambda,t-1}^{\leq}(u)+\lambda,\,\overline{C}_{\lambda,t-1}(u)\\}.$
(21)

The same sub-routines described in Section B.2 can be used to implement the
algorithm below, which solves the penalized reduced isotonic regression
problem (20).

1: Input: data set $\mathbf{y}\in\mathbb{R}^{n}$, penalty constant
$\lambda\geq 0$.

2: Output: vectors of optimal segment means $U\in\mathbb{R}^{n}$ and ends
$T\in\\{1,\dots,n\\}^{n}$

3: Compute min $\underline{y}$ and max $\overline{y}$ of $\mathbf{y}$.

4:
$\overline{C}_{\lambda,1}\leftarrow\text{OnePiece}(y_{1},\underline{y},\overline{y})$

5: for data points $t$ from 2 to $n$: // dynamic programming

6:
$\text{min\\_prev}\leftarrow\lambda+\text{MinLess}(t-1,\overline{C}_{\lambda,t-1})$

7:
$\text{min\\_new}\leftarrow\text{MinOfTwo}(\text{min\\_prev},\overline{C}_{\lambda,t-1})$

8:
$\overline{C}_{\lambda,t}\leftarrow\text{min\\_new}+\text{OnePiece}(y_{t},\underline{y},\overline{y})$

9:
$u^{*},t^{\prime},u^{\prime}\leftarrow\text{ArgMin}(\overline{C}_{\lambda,n})$
// begin decoding

10: $i\leftarrow 1;\,U_{i}\leftarrow u^{*};\,T_{i}\leftarrow t^{\prime}$

11: while $t^{\prime}>0$:

12: if $u^{\prime}<\infty$: $u^{*}\leftarrow u^{\prime}$

13:
$t^{\prime},u^{\prime}\leftarrow\text{FindMean}(u^{*},\overline{C}_{\lambda,t^{\prime}})$

14: $i\leftarrow i+1;\,U_{i}\leftarrow u^{*};\,T_{i}\leftarrow t^{\prime}$

Algorithm 3 Generalized Functional Pruning Optimal Partitioning (GFPOP) for
penalized reduced isotonic regression

Algorithm 3 begins by computing the min/max (line 3). The main storage of the
algorithm is $\overline{C}_{\lambda,t}$, which should be initialized as an
array of $n$ empty FunctionPieceList objects.

The dynamic programming recursion in this algorithm has only one for loop over
data points $t$ (line 5). The penalty constant $\lambda$ is added to all of
the function pieces that result from MinLess (line 6), before computing
MinOfTwo (line 7). The last step of each dynamic programming update is to add
the cost of the new data point (line 8).

The decoding process on lines 9–14 is essentially the same as the GPDPA
(Algorithm 1). The last segment mean and second to last segment end are first
stored on line 10 in $(U_{1},T_{1})$. For each other segment $i$, the mean and
previous segment end are stored on line 14 in $(U_{i},T_{i})$. Note that there
should be space to store $(U_{i},T_{i})$ parameters for up to $n$ segments.
However, there are usually less than $n$ segments, and the algorithm should
return a special flag for unused parameters, for example
$(U_{i}=\infty,T_{i}=-1)$.

The time complexity of Algorithm 3 is $O(nI)$, where $I$ is the time
complexity of the MinLess and MinOfTwo sub-routines. As in the GPDPA, the time
complexity of these sub-routines is linear in the number of intervals
(FunctionPiece objects) that are used to represent the
$\overline{C}_{\lambda,t}$ cost functions. Since the number of intervals in
real data is typically $I=O(\log n)$ (see Section 5.1), the overall time
complexity of Algorithm 3 is on average $O(n\log n)$.

###  B.5 Generalized Functional Pruning Optimal Partitioning Solvers

The GFPOP algorithm can solve problems with more general constraints than
reduced isotonic regression. Let $G=(V,E)$ be a directed graph that represents
the model constraints (for examples see Figure 5). The vertices
$V=\\{1,\dots,|V|\\}$ can be represented as integers, one for every distinct
state. The edges $E=\\{1,\dots,|E|\\}$ is another set of integers, each of
which represents one of the possible changes between states. Each edge/change
$c\in E$ has corresponding data
$(\underline{v}_{c},\overline{v}_{c},\lambda_{c},g_{c})$ which specifies a
transition from state $\underline{v}_{c}$ to state $\overline{v}_{c}$, with a
penalty of $\lambda_{c}\in\mathbb{R}_{+}$, and a constraint function
$g_{c}:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$.

Unconstrained
1$g_{0},\lambda$

Reduced
isotonic
regression
1$g_{\uparrow},\lambda$

Peak detection
peakbackground$g_{\uparrow},\lambda_{\uparrow}$$g_{\downarrow},\lambda_{\downarrow}$

Unimodal regression
up/downupdown$g_{\uparrow},\lambda_{1}$$g_{\uparrow},\lambda_{2}$$g_{\downarrow},\lambda_{3}$$g_{\downarrow},\lambda_{4}$

Figure 5: Examples of state graphs for four models. Nodes represent states and
edges represent changes. Each change has a corresponding penalty $\lambda$,
and a function $g$ that determines what types of changes are possible ($g_{0}$
any change, $g_{\uparrow}$ non-decreasing, $g_{\downarrow}$ non-increasing).
Even if there is no edge from a node to itself, it is still possible to stay
in the same state without introducing a changepoint and penalty. Note that the
Unimodal regression state X should be interpreted as “can change X” e.g.
up/down means “can change up/down.”

Reduced Isotonic Regression
$C_{1,t}$$C_{1,t+1}$$C_{1,t}$$C_{1,t}^{g_{\uparrow}}+\lambda$

Peak detection
$C_{\text{peak},t}$$C_{\text{bkg},t}$$C_{\text{peak},t+1}$$C_{\text{bkg},t+1}$$C_{\text{peak},t}$$C_{\text{peak},t}^{g_{\downarrow}}+\lambda_{\downarrow}$$C_{\text{bkg},t}$$C_{\text{bkg},t}^{g_{\uparrow}}+\lambda_{\uparrow}$

Unimodal regression
$C_{\text{up},t}$$C_{\text{up/down},t}$$C_{\text{down},t}$$C_{\text{up},t+1}$$C_{\text{up/down},t+1}$$C_{\text{down},t+1}$$C_{\text{up},t}$$C_{\text{up/down},t}$$C_{\text{down},t}$$C_{\text{up},t}^{g_{\uparrow}}+\lambda_{1}$$C_{\text{up/down},t}^{g_{\uparrow}}+\lambda_{2}\hskip
22.76228pt$$C_{\text{up/down},t}^{g_{\downarrow}}+\lambda_{3}$$C_{\text{down},t}^{g_{\downarrow}}+\lambda_{4}$

Figure 6: Computation graphs which represent the dynamic programming updates
(23) for the state graph models in Figure 5. Nodes represent cost functions
$C_{s,t}$ in each state $s$ at data $t$ and $t+1$. Edges represent inputs to
the $\min\\{\\}$ operation (solid for a changepoint, dotted for no change).
For example in reduced isotonic regression,
$C_{1,t+1}(u)=\ell(y_{t+1},u)+\min\\{C_{1,t}(u),C_{1,t}^{g_{\uparrow}}(u)+\lambda\\}$.

In the optimization problem below we also allow $c=0$, which implies no
penalty $\lambda_{0}=0$, and means no change:

$\displaystyle\operatorname*{minimize}_{\begin{subarray}{c}\mathbf{m}\in\mathbb{R}^{n},\
\mathbf{s}\in V^{n}\\\ \mathbf{c}\in\\{0,1,\dots,|E|\\}^{n-1}\end{subarray}}$
$\displaystyle\ \
\sum_{t=1}^{n}\ell(y_{t},m_{t})+\sum_{t=1}^{n-1}\lambda_{c_{t}}$ (22) subject
to $\displaystyle\ \ c_{t}=0\Rightarrow m_{t}=m_{t+1}\text{ and
}s_{t}=s_{t+1}$ $\displaystyle\ \ c_{t}\neq 0\Rightarrow
g_{c_{t}}(m_{t},m_{t+1})\leq 0\text{ and
}(s_{t},s_{t+1})=(\underline{v}_{c_{t}},\overline{v}_{c_{t}}).$

If some states are desired at the start or end, then those constraints
$s_{1}\in\underline{S},s_{n}\in\overline{S}$ can also be enforced. To compute
the solution to this optimization problem, we propose the following dynamic
programming algorithm.

Let $C_{s,t}(u)$ be the optimal cost with mean $u$ and state $s$ at data point
$t$. This quantity can be recursively computed using dynamic programming. The
initialization for the first data point is $C_{s,1}(u)=\ell(y_{1},u)$ for all
states $s$. The dynamic programming update rule for all data points $t>1$ is

$C_{s,t}(u)=\ell(y_{t},u)+\min\\{M_{s,t-1}(u),\,C_{s,t-1}(u)\\},$ (23)

where the minimum cost of all possible changes to state $s$ from time point
$t-1$ is

$M_{s,t-1}(u)=\min_{c\in
E_{s}}C^{g_{c}}_{\underline{v}_{c},t-1}(u)+\lambda_{c},$ (24)

and the set of all changes going to state $s$ is

$E_{s}=\\{c\in E\mid\overline{v}_{c}=s\\}.$ (25)

The computations required for the dynamic programming updates (23) can be
visualized using a computation graph (Figure 6).

The pseudocode for the algorithm which implements the dynamic programming
updates (23) is stated below.

1: Input: data $\mathbf{y}$ and weights $\mathbf{w}$ (both size $n$), number
of vertices/states $|V|$, starts $\underline{S}\subseteq V$, edges/transitions
$E$.

2: Allocate $|V|\times n$ array of optimal cost functions $C_{s,t}$, each
initialized to NULL.

3: for $t$ from $1$ to $n$:

4: if $t==1$:

5: for $s$ in $\underline{S}$: // initialize cost for all possible starting
states

6: $C_{s,t}\leftarrow\text{InitialCost}(y_{t},w_{t})$

7: else:

8: for $s$ from $1$ to $|V|$:

9: if $C_{s,t-1}$ is NOT NULL: // previous cost in this state has been
computed

10: $C_{s,t}\leftarrow C_{s,t-1}$ // cost of staying in this state (no change)

11: for ($\underline{v}$, $\overline{v}$, $\lambda$, ConstrainedCost) in $E$:

12: if $C_{\underline{v},t-1}$ is NOT NULL: // previous cost has been computed

13:
$\text{cost\\_of\\_change}\leftarrow\text{ConstrainedCost}(C_{\underline{v},t-1})$

14: $\text{cost\\_of\\_change.set}(\underline{v},t-1)$

15: $\text{cost\\_of\\_change.addPenalty}(\text{$\lambda$})$

16: if $C_{\overline{v},t}$ is NULL:

17: $C_{\overline{v},t}\leftarrow\text{cost\\_of\\_change}$

18: else:

19:
$C_{\overline{v},t}\leftarrow\text{MinOfTwo}(C_{\overline{v},t},\text{cost\\_of\\_change})$

20: for $s$ from $1$ to $|V|$:

21: if $C_{s,t}$ is NOT NULL:

22: $C_{s,t}\text{.addDataPoint}(y_{t},w_{t})$

23: Output: $|V|\times n$ array of optimal cost functions $C_{s,t}$.

Algorithm 4 Generalized Functional Pruning Optimal Partitioning Algorithm,
Dynamic Programming (GFPOP-DP)

The algorithm above performs several checks if $C_{s,t}$ is NULL or not (lines
9, 12, 16, 21). All costs are initialized as NULL (line 2). After having
performed the cost update for data $t$, a NULL cost $C_{s,t}$ means that state
$s$ is not feasible at data $t$. For each constraint function $g$ there is a
corresponding ConstrainedCost sub-routine that is mentioned on lines 11 and 13
(e.g. no constraint $g_{0}$ MinUnconstrained, non-decreasing change
$g_{\uparrow}$ MinLess, non-increasing change $g_{\downarrow}$ MinMore).

The average time and space complexity of Algorithm 4 is $O(|V|nI)$ where $|V|$
is the number of states and $I$ is the average number of of intervals stored
in the $|V|\times n$ array of $C_{s,t}$ cost functions. We observed that
$I=\log n$ in the empirical tests of the peak detection model on ChIP-seq data
(Section 5.1), so we expect that the average time complexity of Algorithm 4 is
$O(|V|n\log n)$.

Note that the algorithm above only performs the dynamic programming. The
decoding of optimal model parameters is achieved using the algorithm below.

1: Output: $|V|\times n$ array of optimal cost functions $C_{s,t}$, ends
$\overline{S}\subseteq V$.

2: Allocate $\mathbf{m}\in\mathbb{R}^{n}$ (mean),
$\mathbf{s}\in\mathbb{Z}^{n}$ (state), $\mathbf{t}\in\mathbb{Z}^{n}$ (segment
end).

3:
$u^{*},s^{*},t^{\prime},s^{\prime}\leftarrow\text{ArgMin}(C_{\cdot,n},\overline{S})$
// begin decoding

4: $i\leftarrow 1;\,m_{i}\leftarrow u^{*};\,s_{i}\leftarrow
s^{*};\,t_{i}\leftarrow n$

5: while $t^{\prime}>0$:

6: $i\leftarrow i+1;\,t_{i}\leftarrow t^{\prime}$

7:
$u^{*},s^{*},t^{\prime},s^{\prime}\leftarrow\text{ArgMin}(C_{s^{\prime},t^{\prime}})$

8: $m_{i}\leftarrow u^{*};\,s_{i}\leftarrow s^{*}$

9: Output: $\mathbf{m}$, $\mathbf{s}$, $\mathbf{t}$.

Algorithm 5 Generalized Functional Pruning Optimal Partitioning Algorithm,
decoding (GFPOP-decode)

## References

  * Ashoor et al. [2013] H. Ashoor, A. Hérault, A. Kamoun, F. Radvanyi, V. Bajic, E. Barillot, and V. Boeva.  HMCan: a method for detecting chromatin modifications in cancer samples using ChIP-seq data.  _Bioinformatics_ , 29(23):2979–2986, 2013.
  * Auger and Lawrence [1989] I. Auger and C. Lawrence.  Algorithms for the optimal identification of segment neighborhoods.  _Bull Math Biol_ , 51:39–54, 1989.
  * Bailey et al. [2013] T. Bailey, P. Krajewski, I. Ladunga, C. Lefebvre, Q. Li, T. Liu, P. Madrigal, C. Taslim, and J. Zhang.  Practical guidelines for the comprehensive analysis of ChIP-seq data.  _PLoS computational biology_ , 9(11), 2013.
  * Bellman [1961] R. Bellman.  On the approximation of curves by line segments using dynamic programming.  _Commun. ACM_ , 4(6):284–, June 1961.
  * Best and Chakravarti [1990] M. Best and N. Chakravarti.  Active set algorithms for isotonic regression; a unifying framework.  _Mathematical Programming_ , 47(1):425–439, 1990\.
  * Haiminen et al. [2008] N. Haiminen, A. Gionis, and K. Laasonen.  Algorithms for unimodal segmentation with applications to unimodality detection.  _Knowledge and Information Systems_ , 14(1):39–57, 2008.
  * Hardwick and Stout [2014] J. Hardwick and Q. Stout.  Optimal reduced isotonic regression.  _arXiv preprint arXiv:1412.2844_ , 2014.
  * Hocking et al. [2015] T. Hocking, G. Rigaill, and G. Bourque.  PeakSeg: constrained optimal segmentation and supervised penalty learning for peak detection in count data.  In _Proc. 32nd ICML_ , pages 324–332, 2015.
  * Hocking et al. [2016] T. Hocking, P. Goerner-Potvin, A. Morin, X. Shao, T. Pastinen, and G. Bourque.  Optimizing chip-seq peak detectors using visual labels and supervised machine learning.  _Bioinformatics_ , 2016.
  * Hoefling [2010] H. Hoefling.  A path algorithm for the fused lasso signal approximator.  _Journal of Computational and Graphical Statistics_ , 19(4):984–1006, 2010.
  * Jackson et al. [2005] B. Jackson, J. Scargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, and T. Tsai.  An algorithm for optimal partitioning of data on an interval.  _IEEE Signal Process Lett_ , 12:105–108, 2005.
  * Johnson [2011] N. Johnson.  _Efficient models and algorithms for problems in genomics_.  PhD thesis, Stanford, 2011.  https://purl.stanford.edu/jq411pj0455.
  * Johnson [2013] N. Johnson.  A Dynamic Programming Algorithm for the Fused Lasso and L0-Segmentation.  _Journal of Computational and Graphical Statistics_ , 22(2):246–260, 2013.
  * Killick et al. [2011] R. Killick, P. Fearnhead, and I. Eckley.  Optimal detection of changepoints with a linear computational cost.  _arXiv:1101.1438_ , Jan. 2011.
  * Maidstone et al. [2016] R. Maidstone, T. Hocking, G. Rigaill, and P. Fearnhead.  On optimal multiple changepoint algorithms for large data.  _Statistics and Computing_ , pages 1–15, 2016.  ISSN 1573-1375.
  * Mair et al. [2009] P. Mair, K. Hornik, and J. de Leeuw.  Isotone optimization in R: pool-adjacent-violators algorithm (PAVA) and active set methods.  _Journal of statistical software_ , 32(5):1–24, 2009.
  * Cleynen and Lebarbier [2014] A. Cleynen and E. Lebarbier.  Segmentation of the Poisson and negative binomial rate models: a penalized estimator.  _ESAIM: PS_ , 18:750–769, 2014.
  * Rigaill [2010] G. Rigaill.  Pruned dynamic programming for optimal multiple change-point detection.  arXiv:1004.0887, 2010.
  * Rigaill [2015] G. Rigaill.  A pruned dynamic programming algorithm to recover the best segmentations with 1 to kmax change-points.  _Journal de la Société Française de la Statistique_ , 156(4), 2015.
  * Rote [unpublished] G. Rote.  Isotonic regression by dynamic programming.  http://www.inf.fu-berlin.de/lehre/WS12/HA/isotonicregression.pdf, unpublished.
  * Schell and Singh [1997] M. Schell and B. Singh.  The reduced monotonic regression method.  _Journal of the American Statistical Association_ , 92(437):128–135, 1997.
  * Scott and Knott [1974] A. Scott and M. Knott.  A cluster analysis method for grouping means in the analysis of variance.  _Biometrics_ , 30:507–512, 1974.
  * Tibshirani et al. [2011] R. Tibshirani, H. Hoefling, and R. Tibshirani.  Nearly-isotonic regression.  _Technometrics_ , 53(1):54–61, 2011.
  * Zhang et al. [2008] Y. Zhang, T. Liu, C. Meyer, J. Eeckhoute, D. Johnson, B. Bernstein, C. Nusbaum, R. Myers, M. Brown, W. Li, et al.  Model-based analysis of ChIP-Seq (MACS).  _Genome Biol_ , 9(9):R137, 2008.