Experiments to support upcoming research papers about binary segmentation.
Finite Sample Complexity Analysis of Binary Segmentation
Toby Dylan Hocking
Binary segmentation is the classic greedy algorithm which recursively
splits a sequential data set by optimizing some loss or likelihood
function. Binary segmentation is widely used for change-point
detection in data sets measured over space or time, and as a
sub-routine for decision tree learning. In theory, using a simple
loss function like Gaussian or Poisson, its asymptotic time complexity
should be extremely fast: for
Paper: Comparing binsegRcpp with other implementations of binary segmentation for change-point detection
Abstract: Binary segmentation is a classic algorithm for detecting change-points in sequential data. In theory, using a simple loss function like Normal or Poisson, binary segmentation should be extremely fast for N data and S segments: asymptotically O(N S) time in the worst case, and O(N log S) time in the best case. In practice, existing implementations can be asymptotically slower, and can return incorrect results. We propose binsegRcpp, an R package which provides a correct C++ implementation, with the expected asymp- totic time complexity. We discuss several important C++ coding techniques, and include detailed comparisons with other implementations of binary segmentation: ruptures, fpop, changepoint, blockcpd, and wbs.
Replication materials: TGZ with data and R script for JSS.
- Figure 1 (models from 2 to 4 segments) and manuscript code/output: source article.Rnw, PDF.
- Figure 2 (DNA copy number example): R code, PNG.
- Figure 3 (cross-validation): R code, PDF.
- Figures 4-6 (list vs multiset, changepoint vs ruptures, throughput): R timing code, R figure code, Figure 4 PNG, Figure 4 tex, Figure 5 PNG, Figure 5 tex, Figure 6 PNG.
- Figure 7 (square loss): R timing code, R figure code, PNG.
- Figure 8 (normal change in mean and variance): R timing code, R figure code, PNG.
- Figure 9 (Poisson): R timing code, R figure code, PNG.
- Figure 10 (L1): R timing code, R figure code, PNG.
- Figure 11 (compare distributions): R figure code, PNG.
To reproduce the analyses and figures in the paper you can use the code below:
- Figure 1: synthetic data which require tie-breaking. R code, PNG figure.
- Figure 2: optimal binary trees. R code, PNG figure.
- Figure 3: number of candidate splits considered in real genomic ChIP-seq data from McGill benchmark. R code 1 makes CSV file, R code 2 makes PNG figure.
- Figure 4: number of candidate splits considered in real genomic copy number data from neuroblastoma benchmark. R code 1 makes CSV file, R code 2 makes PNG figure.
- Compare with Julia https://github.com/STOR-i/Changepoints.jl/blob/master/src/BS.jl maybe using https://cran.r-project.org/web/packages/JuliaCall/readme/README.html
- MATLAB https://www.mathworks.com/help/signal/ref/findchangepts.html#mw_748d5529-e05a-4b55-a3fe-2a12a5772d22
- NAG licensing and installation too complicated, https://www.nag.com/numeric/nl/nagdoc_latest/clhtml/g13/g13ndc.html, https://www.nag.com/content/calling-nag-fortran-routines-r, https://www.nag.com/content/nagfwrappers-r-package
- Try ruptures cumsum deepcharles/ruptures#245 (comment)
changepoint.bug.R makes figure for rkillick/changepoint#86
figure-compare-ruptures.R makes the figure below which plots asymptotic reference lines on top of the empirical median timing,
The figure above shows that
- first column: changepoint is clearly cubic,
- second/third columns: ruptures is clearly quadratic in worst case, and best case seems to be somewhere between quadratic and log-linear.
figure-timings-versions-data.R makes figure-timings-versions-data.csv
figure-timings-versions.R reads that and makes
figure-synthetic.R makes 
figure-neuroblastoma-iterations-data.R makes figure-neuroblastoma-iterations-data.csv and figure-neuroblastoma-iterations-bounds.csv
figure-neuroblastoma-iterations.R makes
figure-mcgill-iterations-data.R makes figure-mcgill-iterations-data.csv
figure-mcgill-iterations.R reads that and makes
Figure above shows that in real data achieves asymptotic best case complexity.
figure-best-heuristics.R makes figures below
figure-optimal-trees.R makes small pictures for slides like this
and big picture below
figure-compare-distributions.R makes
Figure above shows small asymptotic slowdown for Laplace median distribution.
figure-splits-loss.R makes new figure
which shows some empirical cumulative counts smaller than the best, which is possible since the “best” case time complexity bound refers to the total number of splits at the end, not at each iteration.
figure-splits-loss.R makes
Figure above shows that 0/1 seq can be used for worst case of l1,mean_norm,poisson loss, but for meanvar_norm we need 0/1/10/11 seq. Also linear data achieves best case for normal losses, and is very close to best case for l1 and poisson. TODO figure out a simple synthetic data sequence which achieves the best case for l1 and poisson.
figure-timings-laplace-data.R makes figure-timings-meanvar_norm-data.csv
figure-timings-laplace.R reads that file and makes
Figure above shows that ruptures looks asymptotically slower in best case.
figure-timings-meanvar_norm-data.R makes figure-timings-meanvar_norm-data.csv
figure-timings-meanvar_norm.R reads that file and makes
Figure above shows that
- blockcpd is about the same as binsegRcpp multiset.
- for worst case changepoint is faster up to very large model sizes, but asymptotically slower.
figure-timings-poisson-data.R makes figure-timings-poisson-data.csv
figure-timings-poisson.R reads that file and makes
Figure above shows that
- blockcpd about the same as binsegRcpp multiset.
- others consistent with other losses.
TODO compare both versions of blockcpd. Also compare with max.segs=n.data since that is what blockcpd does?
figure-neuroblastoma.R makes the figure below, which shows a real data set for which there are differences between binsegRcpp and ruptures/changepoint.
ruptures_bug.py and changepoint.bug.R used to report issues, deepcharles/ruptures#242 and rkillick/changepoint#69
figure-timings-data.R makes figure-timings-data.csv
figure-timings.R reads that and makes
Figure above was created using synthetic data which achieve the best/worst case time complexity of the binary segmentation algorithm. For each data set of a given size N in {2^2=4,8,16,32,…,2^20=1,048,576}, we run binary segmentation up to a max of N/2 segments (and not going to a larger N if the algo/case resulted in a time greater than 100 seconds). The timings suggest that changepoint R package uses a cubic algorithm (three nested for loops) whereas binsegRcpp uses an algorithm which is log-linear in the best case, and quadratic in the worst case. The ruptures python module seems to be asymptotically faster than changepoint but slower than binsegRcpp, maybe quadratic?
Figure above shows that loss for binsegRcpp is always less than loss for others, suggesting that there are bugs in the other implementations.
figure-select-segments-data.R computes simulations using a variety of model selection criteria, saving results to figure-select-segments-data.csv
figure-select-segments.R reads that result CSV file and makes
