Benchmarks: Runtime, False Alarms and Coverage

Youzhi Yu
University of Chicago

Fifty methods share one interface in this package, which makes it easy to swap one for another and easy to forget that they do not cost the same thing. This article answers three questions a reader choosing a method has, from the measurements that ship with the package:

1. How long does it take, and at what series length does it stop being usable at all?
2. How often does it fire when there is nothing there?
3. Do the intervals cover?

How these numbers were produced, and what they are not

The tables come from the cpt_runtimes, cpt_null_sizes and cpt_calibration data sets, generated by the scripts in data-raw/ for this release on one Linux x86-64 machine (AMD EPYC) with R 4.4.1 and the engine versions current then, and read here when the page is built. The runtimes are one fit per cell: they answer “will it finish?”, and are not a basis for saying one engine is some factor faster than another. Re-measure on your own hardware before quoting a time.

1. Runtime

Simulated Gaussian series with four mean changes of two noise standard deviations (three coordinates for the multivariate engines), each fit in its own R process after a warm-up call, with a two-minute cap (five at a million). A length is attempted only when the previous one finished in under ten seconds.

rt <- cpt_runtimes
rt$cell <- ifelse(rt$status == "ok", format(round(rt$seconds, 2), nsmall = 2),
                  rt$status)
wide <- reshape(as.data.frame(rt[, c("method", "n", "cell")]),
                idvar = "method", timevar = "n", direction = "wide")
names(wide) <- c("method", "1e3", "1e4", "1e5", "1e6")[seq_along(wide)]
wide[is.na(wide)] <- ""
fastest <- rt[rt$n == 1e6 & rt$status == "ok", ]
# Fastest at a million first, then the rest alphabetically.
rank_1e6 <- match(wide$method, fastest$method[order(fastest$seconds)])
wide <- wide[order(is.na(rank_1e6), rank_1e6, wide$method), ]
knitr::kable(wide, row.names = FALSE,
             caption = "Seconds per fit (blank: not attempted)")

Seconds per fit (blank: not attempted)

method	1e3	1e4	1e5	1e6
mosum	0.01	0.01	0.11	0.54
amoc	0.01	0.01	0.04	0.59
binsegrcpp	0.01	0.01	0.05	0.62
fpop	0.01	0.01	0.06	0.83
pelt	0.01	0.01	0.07	1.01
pilliat	0.01	0.02	0.12	1.56
fastcpd	0.01	0.06	0.44	4.01
decafs	0.09	0.15	0.63	7.04
esac	0.01	0.08	1.28	12.71
segmented	0.09	0.21	2.29	13.95
not	0.09	0.24	1.83	17.30
geomcp	0.11	0.23	1.82	18.36
wbs	0.06	0.30	2.55	25.84
tguh	0.20	0.47	3.43	33.95
idetect	0.02	0.32	5.42	60.98
bcp	0.64	5.30	40.27
beast	0.17	1.00	15.35
bfast	0.96	56.48
binseg	0.01	0.01	0.07	error
bocpd	4.53	timeout
buishand	19.91
cpm	0.01	0.13	11.76
cpop	3.74	timeout
ecp	20.43
envcpt	1.08	41.00
fcov	timeout
fmean	31.58
hdcov	13.38
hdreg	error
hsmuce	0.05	0.11	timeout
inspect	0.14	0.54	error
kcp	58.37
kwc	0.16	timeout
network	73.22
np	0.04	1.90	106.07
npmojo	2.07	timeout
nsp	8.00	108.72
ocd	timeout
pettitt	0.04	0.45	35.63
segneigh	0.16	10.95
smuce	0.90	0.10	timeout
sn	2.19	96.05
snht	33.73
strucchange	3.18	timeout
taylor	13.42
var	38.27
wbs2	1.10	12.22
wbsts	0.27	6.50	timeout

15 methods fit a million points inside the cap. The pruned optimal-partitioning family and the single-change test are fast enough that length is not a consideration; the random-interval searches cost more because they fit many sub-intervals by construction, which is also why they find short changes the single-pass methods miss. The methods that stop at 1,000 or 10,000 are quadratic or worse in the length (ecp, bocpd, segneigh’s cost matrix), solve a harder problem (cpop’s continuous piecewise-linear fit), or were written for the lengths their fields see (strucchange, taylor). None is broken; they are methods for short series, and choosing one for a long series is a modelling decision rather than a performance problem to tune around. cpt_methods() carries each method’s measured cost and max_n, and cpt_recommend(n = ) leaves out the ones that did not finish at the length you give it.

2. False alarms on pure noise

Independent Gaussian noise with no changepoint anywhere, for the univariate engines the runtime table shows finishing each length quickly: 50 replicates at 1,000, 10 at 10,000 and 6 at 100,000. Every changepoint reported is a false alarm.

ns <- cpt_null_sizes
ns$cell <- paste0(format(ns$mean_fp, nsmall = 2), " (", ns$reps, ")")
wide_ns <- reshape(as.data.frame(ns[, c("method", "n", "cell")]),
                   idvar = "method", timevar = "n", direction = "wide")
names(wide_ns) <- c("method", "1e3", "1e4", "1e5")[seq_along(wide_ns)]
wide_ns[is.na(wide_ns)] <- ""
knitr::kable(wide_ns[order(wide_ns$method), ], row.names = FALSE,
             caption = "Mean false changepoints per series (replicates)")

Mean false changepoints per series (replicates)

method	1e3	1e4	1e5
amoc	0.00 (50)	0.00 (10)	0.00 (6)
bcp	0.02 (50)	0.00 (10)
beast	0.00 (50)	0.00 (10)
bfast	0.00 (50)
binseg	0.00 (50)	0.00 (10)	0.00 (6)
binsegrcpp	0.00 (50)	0.00 (10)	0.00 (6)
bocpd	0.00 (50)
buishand	0.08 (50)
cpm	0.34 (50)	0.40 (10)
cpop	0.00 (50)
decafs	0.00 (50)	0.00 (10)	0.00 (6)
ecp	0.16 (50)
envcpt	0.00 (50)
fastcpd	0.00 (50)	0.00 (10)	0.00 (6)
fpop	0.00 (50)	0.00 (10)	0.00 (6)
hsmuce	0.28 (50)	0.40 (10)
idetect	0.06 (50)	0.00 (10)
mosum	0.04 (50)	0.20 (10)	0.00 (6)
not	0.04 (50)	0.00 (10)	0.00 (6)
np	1.94 (50)	5.10 (10)
nsp	0.02 (50)
pelt	0.00 (50)	0.00 (10)	0.00 (6)
pettitt	0.06 (50)	0.00 (10)
segmented	1.00 (50)	1.00 (10)	1.00 (6)
segneigh	0.02 (50)	0.00 (10)
smuce	0.28 (50)	0.40 (10)
sn	0.06 (50)
snht	0.04 (50)
strucchange	0.00 (50)
taylor	0.32 (50)
tguh	0.00 (50)	0.00 (10)	0.00 (6)
wbs	0.04 (50)	0.00 (10)	0.00 (6)
wbs2	0.58 (50)
wbsts	0.02 (50)	0.10 (10)

Most engines raise none at their defaults. The rows that do are worth reading correctly:

• segmented always reports the breakpoint it is asked for (npsi = 1), change or no change; test it with cpt_test() before reporting it.
• smuce and hsmuce default to alpha = 0.5, the level that is optimal for estimating the number of changes, and which by design overestimates it in up to half of all series. Pass a smaller alpha for a test.
• cpm is a sequential monitor used in batch mode: its threshold is a false-alarm rate per observation, so it reports about n / arl0 false changepoints on any series. Since 0.6.0 the default arl0 scales with the length (five times n, up to the 50,000 cpm ships), and the result records the expected count.
• np over-detects on Gaussian noise, and more so as the series grows.

Zero false alarms in a handful of replicates is consistent with a true rate anywhere below a few tenths per series, so a zero here means “not obviously anti-conservative” rather than “exact”. And it is measured on independent noise: under autocorrelation almost every one of those zeros becomes a positive count (cpt_noise_benchmark, and the pitfalls article).

3. Interval coverage

The calibration article has the measured coverage of each provenance of cpt_confint():

subset(cpt_calibration, kind == "coverage")[, c("call", "nominal",
                                                 "realised", "mcse")]
#> # A tibble: 3 × 4
#>   call                                                 nominal realised   mcse
#>   <chr>                                                  <dbl>    <dbl>  <dbl>
#> 1 "cpt_detect(x, method = \"strucchange\")"               0.95    0.953 0.0047
#> 2 "cpt_detect(x, method = \"smuce\", alpha = 0.05)"       0.95    1     0.0005
#> 3 "cpt_confint(cpt_detect(x), method = \"bootstrap\")"    0.95    0.968 0.0056

Two provenances are not in it yet. NSP’s regions hold a global guarantee across all intervals and survive selection, so they are conservative by design. Posterior intervals from bcp can be close to vacuous at 95%: its posterior probability of a change is sharply peaked but has diffuse tails, so pushing the credible level out sweeps in most of the series. cpt_confint() warns when that happens, and its source column says which provenance each interval has, because a posterior credible interval and a confidence interval are not interchangeable.