text:
\section*{Acknowledgements}
Solt Kov\'acs and Peter B\"uhlmann have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 786461 CausalStats - ERC-2017-ADG). Axel Munk and Housen Li are funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - EXC 2067/1- 390729940. The authors would like to thank Alexandre M\"{o}sching for careful comments.

\bibliographystyle{apalike}%

\section{Introduction}
\label{Introduction}
Change point detection, also called break point detection, tackles the problem of estimating the location of abrupt structural changes for ordered data, such as ordering by time or space (e.g.~the position on the genome). One can distinguish between online (sequential) and offline (retrospective) detection problems. We will primarily focus on the latter setup where all ordered observations are available beforehand and only point to online detection in connection with our methods and results for the detection of a single change point. The goal is usually twofold: to estimate the number and to estimate the location of change points. Applications include detecting changes in copy number variation \citep{CBS,ZhSi07}, ion channels \citep{Hotz_ionchannel}, financial time series \citep{BaPe98,Kim_finance,DaHo12}, climate data \citep{Reeves_climatology}, environmental monitoring systems \citep{Londschien}, among many others.
One can further differentiate with respect to parametric versus non-parametric and univariate versus multivariate methods, or whether there is dependency across consecutive observations, see for example recent reviews of \cite{review_Truong} and \citet{review_Niu}.

We focus on computational and algorithmic aspects of change point detection. Two common algorithmic approaches are optimal partitions via dynamic programming (see \citealp{jackson2005algorithm} and \citealp{FKLW08}) and greedy procedures (e.g.~binary segmentation (BS, \citealp{Vostrikova}) and its variants). The former approaches are mainly investigated for univariate data. For references for some exceptions, e.g. involving high-dimensional data, see references further below. In the univariate scenario, typical approaches include  e.g.~$\ell_0$ penalization methods (\citealp{BKLMW09}; PELT by \citealp{Killick_etal} and FPOP by \citealp{Maidstone}) and multiscale methods (SMUCE by \citealp{Frick_etal} and FDRSeg by \citealp{LMS16}), which are known to be statistically minimax optimal. However, finding the optimal partition requires in the worst case at least quadratic run time. In contrast, BS is typically faster and easier to adapt to more general scenarios, but worse in terms of estimation performance than methods finding the optimal partitioning. \cite{Fryzlewicz_WBS} proposed wild binary segmentation (WBS) and \cite{Baranowski} proposed a similar narrowest over threshold (NOT) method. Both of these BS-like procedures improve on estimation performance of BS, but lose some of its computational efficiency. Recently, \cite{SeedBS} proposed seeded binary segmentation (SeedBS) combining the best of both worlds, i.e.~improved estimation performance of BS, while keeping computational efficiency.

Already for simple univariate cases, increasingly larger data sets with long time series recently led to the development of more efficient (univariate) approaches (\citealp{Maidstone, intelligent_sampling, Fryzlewicz_WBS2, SeedBS}). Computational issues are even much more pronounced for multivariate problems. This is in particular true for recently emerging high-dimensional change point detection approaches (e.g.~\citealp{LeonBuhl, RoyMichailidis, GibberdNelson, Gibberd_Roy, Bybee_Atchade_JMLR, Avanesov_theory, Wang_highdim_mean_change, Wang_Willet_regression, Londschien, Wang_Yu_Rinaldo_theory}) or in some recent multivariate non-parametric change point detection methods (e.g.~\citealp{Padilla_Yu_Wang_Rinaldo_multivariate_nonparametric}).
Many of these approaches rely on regularized estimators such as the lasso \citep{lasso_original} or the graphical lasso \citep{glasso}, where even a single fit is costly. Performing a full grid search in order to find a single split point requires as many fits as there are observations. Even with warm-starts, neighbouring fits with one additional observation are not straightforward to update (unlike e.g.~means in univariate cases) and the number of fits is the main driver of computational cost. Thus, for a few hundred or thousand observations, full grid search based methods (including BS and its variants and even more costly dynamic programming based approaches) can be very slow, beyond what is acceptable. This is a main motivation for our work here, namely to avoid too many fits with piecewise constant (one or more change points) model structure.

\subsection{Our contribution and related work}
The key idea is to replace the exhaustive search by an \emph{adaptive search} that dynamically determines the next search location given the previous ones. We will show that this leads to only a logarithmic  number of evaluations, while still being statistically efficient in certain scenarios. Such an adaptive search is possible mainly because of a rather general observation. In numerous change point detection problems, population gain functions (describing the expected gain when splitting into two parts at a given split point) have a specific piecewise quasiconvex structure such that local maxima correspond to change points, and the presence of a single change point implies only one global maximum.

In the single change point case, we propose for the first time, to our best knowledge, an algorithm with $O(\log T)$ function evaluations leading to asymptotically (nearly) minimax optimal change point detection as opposed to full grid search based methods requiring at least $O(T)$ evaluations for univariate Gaussian changing mean problems. This result is surprising and fundamental, as essentially one loses nothing in terms of detectable signals and in terms of localization error (in an asymptotic sense) compared to full grid search. While a direct and practical application could be in developing novel online change point detection methods,  we rather focus here on extending the methodology to asymptotically minimax optimal detection with sublinear number of evaluations also for multiple change points. In univariate Gaussian changing mean problems, our methodology is significantly faster only if cumulative sums have been pre-computed, a reasonable assumption e.g.~in online change point detection. This also suggests to store data in cumulative sums format for offline analysis in this particular case. Otherwise we are limited by the $O(T)$ cost of calculating the cumulative sums, which then is only a marginal improvement compared to currently available $O(T\log T)$ time methods. However, the situation is markedly different in more complex models where the number of model fits involved in each evaluation of the gain function are the main driver of computational cost. A motivating example is given in~Section~\ref{motivating_example}.

Some problem specific solutions for gaining computational speedup in certain high-dimensional setups arose, e.g. for changing (Gaussian) graphical models (\citealp{Hallac_network, Bybee_Atchade_JMLR}) or changing linear regression coefficients (\citealp{Kaul_2step, Kaul_multiple}), but these may not be easy to adapt to other scenarios. The following two proposals may be seen as most related to our approach. Our idea of avoiding a high number of model fits is vaguely related to the procedure of \citet{Kaul_2step} for change point detection in high-dimensional linear regression. For an initial split point they fit appropriate models which are then kept for evaluating (an approximation of) the gain function for candidate split points on the full grid, but without updating the costly high-dimensional fits.
However, in scenarios where the evaluation of the gain on the full grid has a comparable cost as the model fit itself (e.g.~calculating and updating means in the univariate Gaussian change in mean setup), their procedure would not lead to any speedup. \citet{intelligent_sampling} on the other hand propose a procedure specific for the univariate case. They thin out the number of evaluations by searching on a small subsample to obtain preliminary estimates, and then use a dense sample in neighbourhoods for the final change point estimates. A rough subsample is not well suited for high-dimensional problems, unlike our approach of keeping all samples but avoiding evaluations on the full grid.

\subsection{A motivating example}
\label{motivating_example}
Instead of problem specific solutions, we first propose a computationally attractive new and general methodology for searching for single change points that we call \emph{optimistic search} (including naive, advanced and combined versions). Figure~\ref{Fig:glasso_single_cpt} illustrates its computational efficiency in an example of change point detection for a $200$-dimensional Gaussian graphical model (based on an estimator discussed in Section~\ref{highdim_simulations} and Appendix~\ref{Appendix_highdim}). The aim is to find the maximum of the gain function (black curve). Evaluating the gain function at a single split point $s \in\{ 1, \ldots, T=2000\}$ requires two graphical lasso \citep{glasso} fits: one for the segment $(1, s]$ and one for $(s, T]$. For finding the maximum, the full grid search, i.e.~considering all possible split points, took roughly $100$ times longer than the shown approximation using naive optimistic search. The latter evaluated the gain only at two initial points (marked by two zeros) and subsequently further $14$ split points (marked by colored numbers). Hence, the number of costly graphical lasso fits was massively reduced compared to the full gain curve (in black). The best guess for the change point is the overall maximum at observation $402$, very close to the actual true change point at observation $400$ (in green). Optimistic search does not know the shape of the black curve in advance and yet after two initial starting points it is capable of adaptively picking only $14$ further evaluations to return a remarkably good candidate at observation $423$ (in blue) with only a slightly worse gain than for the global maximum. When searching for multiple change points, optimistic search can be combined flexibly with existing algorithms, e.g.~BS or its variants including SeedBS and WBS. The respective optimistic algorithms offer considerable speedup compared to their traditional full grid search based counterparts, similar to the single change point example shown in Figure~\ref{Fig:glasso_single_cpt}. As an advantage of combining with BS-type procedures, the methodology remains flexible and general to be used for various change point detection scenarios. As opposed to problem specific solutions, our methods are computationally efficient in the range from other high-dimensional problems or complex time series models with costly fits all the way to classical univariate change point detection.

\begin{figure}
\centering
\includegraphics[width=0.82\textwidth]{figures/fig1.pdf} \\
\caption{Finding a single change point with full grid search (black) and naive optimistic search (red) in a $200\times 200$-dimensional covariance change example with underlying graphical lasso fits. Optimistic search starts with two initial evaluations marked by the two zeros and then evaluates further $14$ split points adaptively, in the order marked by the respective colored numbers shown in the zoomed in part~b). The true underlying change point at observation~$400$ is marked in green and the final candidate returned by optimistic search at observation $423$ in blue. The overall maximum of the black gain curve found by full grid search is at observation $402$.}
\label{Fig:glasso_single_cpt}
\end{figure}

\subsection{Outline and announcing our results}
The basis of our methodology is the optimistic search which we introduce in Section~\ref{Methodology_single}. It is capable of finding a local maximum of the population gain function with $O(\log T)$ evaluations for an interval with $T$ observations. We consider the classical univariate Gaussian change in mean setup, detailed in Section~\ref{s:gauss}. In the scenario when only a single change point is present, the naive version of optimistic search (Theorem~\ref{th:noss}) is able to detect the only change point if
$$
\text{minimal jump size}\, \times \,\text {minimal segment length} \,  \gtrsim \,\sqrt{\log (T)/ T}
$$
which is suboptimal. In contrast, the advanced and the combined optimistic searches (Theorem~\ref{th:aoss}) are able to detect the change point under the weakest condition (up to a possibly removable log factor)
\begin{equation}\label{eq:opt}
\text{minimal jump size} \,\times\, \sqrt{\text{minimal segment length}}\, \gtrsim \,\sqrt{\log(T)/T}
\end{equation}
thus being optimal in minimax sense. Moreover, all three (naive, advanced and combined) optimistic searches achieve the localization error at the sampling rate up to a possible log factor.
In Section~\ref{Methodolgy_multiple} we extend the methodology for multiple change points in combination with BS and SeedBS with differing advantages and disadvantages in terms of estimation performance guarantees and the achievable computational speedups. In particular, the combination of SeedBS and the advanced or the combined optimistic search leads to a minimax optimal performance, namely under the same condition as \eqref{eq:opt} the number of change points is identified correctly, and meanwhile the location of each change point is estimated at the best rate that is available in literature, see Theorem~\ref{th:osbs}.

The number of evaluations can be even sublinear in terms of sample size~$T$, if some prior knowledge on the minimal segment length is available. For instance, if the minimal segment length is known to be of order $T^{-\omega}$ with $\omega \in [0,1]$, then we need $O(\min\{T^{\omega}\log T,\, T\})$ evaluations. In particular, SeedBS combined with optimistic search has a linear run time in the worst case, even without any prior knowledge on the minimal spacing between change points. In Section~\ref{Simulations} we show simulations both in single and multiple change point cases for the univariate Gaussian setup demonstrating that estimation performance remains remarkably close to the respective full grid search variants. We present simulations from an example with high-dimensional covariance matrices in Section~\ref{highdim_simulations}, as we believe our methodology is most useful in such scenarios with costly fits. We conclude with a discussion in Section~\ref{Discussion}. Additional material as well as proofs are given in the Appendix.

\paragraph{Notation} For a real number $r$, we define downward rounding as $\lfloor r \rfloor = \max\{n\in \mathbb{Z} \,:\, n \le r\}$ and upward rounding as $\lceil r \rceil = \min\{n \in \mathbb{Z}\,:\,n \ge r\}$. For two sequences of positive real numbers $\{a_T\}_{T=1}^\infty$ and $\{b_T\}_{T=1}^\infty$, we write $a_T \lesssim b_T$, or equivalently $a_T = O(b_T)$, if $\limsup_{T \to \infty} a_T / b_T < \infty$, and write $a_T \asymp b_T$ if both $a_T \lesssim b_T$ and $b_T \lesssim a_T$\,.

\section{Gaussian mean shifts with constant variance}\label{s:gauss}
Throughout this paper, as a toy example we will consider Model~\ref{Gaussian_setup}. Recent proposals for this classical and well-studied problem can be found in \cite{Frick_etal, Fryzlewicz_WBS} and \cite{SeedBS} with references therein providing a good overview. We believe that the understanding for this setup can be easily generalized to other scenarios.

\begin{Model}
\label{Gaussian_setup}
Assume that observations $X_1, \ldots, X_T$ are independent and
\begin{equation*}
\begin{split}
 X_{\cpt_{0}T+1}(=X_1), \ldots, X_{\cpt_1 T} & \sim \mathcal{N}(\mu_0, \sigma^2)\,, \\
& \vdots	\\
 X_{\cpt_{\ncp}T+1}, \ldots, X_{\cpt_{\ncp+1}T}(=X_{T}) 	& \sim \mathcal{N}(\mu_{\ncp}, \sigma^2)\,,
\end{split}
\end{equation*}
where $\{\cpt_i \,:\, i = 1, \ldots, \ncp\}$ gives the location of change points satisfying
$$
0 = \cpt_0 < \tau_1 < \cdots < \cpt_{\ncp+1} = 1\quad \text{and}\quad \cpt_i T \in \N\,,
$$
means $\mu_i \neq \mu_{i-1}$ for $i=1,\ldots,\ncp$ give the levels on segments, and the common standard deviation $\sigma > 0$ is known.
Assume w.l.o.g.~$\sigma = 1$.
Moreover, define the minimal segment length $\lambda$ as
$$
\lambda \equiv \lambda_T = \min_{i = 0,\ldots,\ncp} (\cpt_{i+1} - \cpt_i)\,,
$$
and the minimal jump size $\delta$ as
$$
\delta \equiv \delta_T = \min_{i = 1,\ldots, \ncp} \delta_i \qquad \text{with} \quad \delta_i = \abs{\mu_i - \mu_{i-1}}\,.
$$
\end{Model}

The goal of change point detection is to estimate the number $\ncp$ and the location $\cpt_i$ of the true underlying change points from realizations $X_1, \ldots, X_{T}\,$. A common criterion for determining the best split point is the CUSUM statistics \citep{Page54}, defined for an interval $(l,r]$ and a split point $s$ as
\begin{equation}
\label{CUSUM}
    \csm_{(l,r]}(s) = \sqrt{\frac{r-s}{n(s-l)}}\sum_{t=l+1}^s X_t -
        \sqrt{\frac{s-l}{n(r-s)}}\sum_{t=s+1}^r X_t\,,
\end{equation}
with integers $0 \leq l < s < r \leq T$ and $n = r-l\,$. The CUSUM statistic is the generalized likelihood ratio test for a single change point at location $s$ in the interval $(l,r]$ against a constant signal. The population counterpart of $|\csm_{(l,r]}(\cdot)|$, i.e.
\begin{equation}\label{eq:pgf}
|\csm_{(l,r]}^*(\cdot)|\, : \,s \mapsto \left|\sqrt{\frac{r-s}{n(s-l)}}\sum_{t=l+1}^s \E{X_t} -
        \sqrt{\frac{s-l}{n(r-s)}}\sum_{t=s+1}^r \E{X_t}\right|
\end{equation}
has its maximum at one of the underlying change points. In noisy cases, the best split point candidate when dividing the segment $(l,r]$ into two parts is the location of the maximal absolute CUSUM statistics
$$    \hat s_{(l,r]} = \argmax_{s\in{\{l+1, \dots, r-1\}}} |\csm_{(l,r]}(s)|\,.
$$
We refer to the function $|\csm_{(l,r]}(\cdot)|$ as a gain function, denoted by $G(\cdot)$, because the square of it describes gains, namely the reductions in squared errors when fitting separate means on the left and right segments for split points in the segment $(\ell,r]$. In general scenarios, gain functions $G(\cdot)$, similar to $|\csm_{(l,r]}(\cdot)|$, measure the benefits of introducing a split point, see Section~\ref{highdim_simulations}. Note that the gain functions $G(\cdot)$ are initially defined on a discrete grid of split points, but for easier interpretability we will implicitly look at their continuously embedded versions (via e.g.~linear interpolation).

\section{Optimistic search for finding a single change point}
\label{Methodology_single}

In this section, we focus on the scenario with a single change point, and introduce three versions of optimistic search in such a scenario, which are going to be at the core of our methodology for multiple change points.

\subsection{Naive optimistic search}

\begin{algorithm}[t]
\caption{naive optimistic search}\label{Alg:OS}
\begin{algorithmic}[1]
\Require $R - L > 2, ~L,R \in \N$; and step size $\nu \in (0, 1)$ with $ 1/2$ by default \\
\textbf{initialize:} $l\leftarrow L$, $r\leftarrow R$ and $s \leftarrow  \floor{(L+\nu R)/(1+\nu)}$
\Function{\textit{nOS}}{$l, s, r \mid \nu, L, R$}
\If {$r-l \leq 5$} \Comment{Stopping condition for recursion}
\State $\hat s_{(L,R]} \leftarrow \argmax\limits_{s\in\{l+1, \dots, r-1\}} G_{(L,R]}(s)$
	\Comment{Search over all points if less than $5$ remain}
\State \textbf{return} $\hat s_{(L,R]}$
\EndIf
\If {$r-s>s-l$} \Comment{Pick a new probe point in larger segment}
\State $w \leftarrow \lceil r - (r - s)\nu \rceil$
\If {$G_{(L,R]}(w) \geq G_{(L,R]}(s)$}
\State \Call{\textit{nOS}}{$s, w, r \mid \nu, L, R$}
\Else
\State \Call{\textit{nOS}}{$l, s, w \mid \nu, L, R$}
\EndIf
\Else
\State $w \leftarrow \lfloor l + (s-l)\nu \rfloor$
\If {$G_{(L,R]}(w) \geq G_{(L,R]}(s)$}
\State \Call{\textit{nOS}}{$l, w, s \mid \nu, L, R$}
\Else
\State \Call{\textit{nOS}}{$w, s, r \mid \nu, L, R$}
\EndIf
\EndIf
\EndFunction
\end{algorithmic}
\end{algorithm}
We describe first the \emph{naive} version of searching within a segment $(L,R]$ in Algorithm~\ref{Alg:OS}. The procedure is similar to the golden section search (\citealp{Kiefer1953, Avriel1966, Avriel1968}), typically used to find the global extremum of unimodal functions.  Algorithm~\ref{Alg:OS} splits an interval into three segments recursively and discards one of the outer segments in each iteration. For unimodal functions with one peak this search recovers the global maximum. If there is only a single change point contained in $(L,R]$, then the (continuously embedded) population gain function in \eqref{eq:pgf} is unimodal with the single peak at the true underlying change point. What happens if there are multiple change points contained in $(L,R]$ is described in Appendix~\ref{Appendix_OBS}. Optimism is required in noisy scenarios, and hence the naming of the method, as noisy counterparts are rather ``wiggly'' functions following the shape given by the underlying population gain function only approximately.

When initializing by calling $\nOS(L, \floor{(L+\nu R)/(1+\nu)}, R \mid \nu, L, R)$ for an interval $(L,R]$, the search first probes the point $s = L + \nu/(1+\nu) \cdot (R-L)$ vs. the point $w = R - \nu/(1+\nu) \cdot (R-L)$ (up to rounding), i.e.~the two first probe points are equally distant from $L$ and $R$ respectively. Depending on the gains at the probe points $s$ and $w$ either $(L,s]$ or $(w,R]$ is then discarded. The possible decisions for the general case when the search is already narrowed down to the sub-interval $(l,r]$ is depicted in Figure~\ref{Fig:optim-search}. Note that in general, the lengths of the two candidate intervals for discarding are not necessarily equal. In case a) we have $G_{(L,R]}(s) < G_{(L,R]}(w)$ and hence the less promising blue area is discarded. In case b) we have $G_{(L,R]}(s) > G_{(L,R]}(w)$ and the red area is discarded. Note that one of the previous probe points will be one of the new boundary points while the other probe point is going to be one of the new probe points in the middle with gain that is thus at least as high as for the new boundary. This leads to a ``triangular structure'' that one probe point in the middle has a higher gain than both boundary points, throughout the search. While giving an intuition why the algorithm works, this property is essential when considering multiple change points (see Appendix~\ref{Appendix_OBS_theory}).

\begin{figure}
\centering
\includegraphics[width=0.85\textwidth]
{figures/fig2.pdf} \\
\caption{Naive optimistic search step $nOS(l, s, r \mid \nu, L, R)$ within the current leftover segment $(l,r]\subseteq(L,R]$ given the previous evaluation at $s$ and step size $\nu$. As $r-s>s-l$, the new probe point $w$ is taken within $(s,r]$ as $w = \lceil r - (r - s)\nu \rceil$. Depending on the gain $G_{(L,R]}(w)$ vs. $G_{(L,R]}(s)$, one either continues with $nOS(s, w, r \mid \nu, L, R)$ (discarding the blue part) or $nOS(l, s, w \mid \nu, L, R)$ (discarding the red part).}
\label{Fig:optim-search}
\end{figure}

\begin{lemma} \label{Lem:nOS_speed}
The naive optimistic search (Algorithm~\ref{Alg:OS}) with a fixed step size $\nu=1/2$ will always terminate within $O(\log(R-L))$ and hence at most $O(\log T)$ steps (i.e., number of evaluations of the gain function).
\end{lemma}

If $\nu=1/2$, the lengths of the three intervals within the current search interval $l$ to $r$ have relative sizes 1:1:1, 1:1:2 or 2:1:1 (up to rounding). Hence, at each step we discard at least a quarter of the current search interval, which then leads to the logarithmic number of evaluations as claimed in Lemma~\ref{Lem:nOS_speed}. Of course, one can easily generalize this statement for $\nu\neq1/2$. We take $\nu=1/2$ as a default, but in general $\nu$ can be interpreted as (relative) step size, expected to reflect some kind of trade off between computational performance and estimation accuracy. The choice of evaluating the last $5$ points remaining (in line $3$) is somewhat arbitrary and could be also set to e.g.~$10$ or~$20$, but empirically this only gives slightly better localization of the split points. In rare cases, when $r-s = s-l$ or $G_{(l,r]}(w) = G_{(l,r]}(s)$ one could also take a (pseudo) random choice or incorporate additional information (e.g.~variance in the segments for Model~\ref{Gaussian_setup}) to decide on the new probe point and the segment to discard.

While population gain functions are smooth in between change points, noisy variants are rather ``wiggly''. We investigate the resulting localization error as follows.

\begin{theorem}[Naive optimistic search]\label{th:noss}
Under Model~\ref{Gaussian_setup} with a single change point, i.e.~$\ncp =1$, we assume that the minimal segment length $\lambda$ and the minimal jump size $\delta$ satisfy
\begin{equation}\label{eq:area}
\delta \lambda \sqrt{T} \ge C_0 \sqrt{\log T}
\end{equation}
for some large enough constant $C_0$. Then:
\begin{enumerate}[i.]
    \item The naive optimistic search (Algorithm~\ref{Alg:OS}) on $(0,T]$ provides an estimated change point $\hat\cpt = \hat s_{(0,T]}/T$ such that
    $$
    \lim_{T \to \infty}\Prob{\abs{\hat\cpt - \cpt} \le C_1 \frac{\log T}{\delta^2 T}} = 1 \qquad \text{with some constant } C_1\,.
    $$
    \item Let $\hat\cpt_i$, $i = 1,\ldots,\hat\ncp$, be the estimated change points (scaled by $1/T$) by the combination of naive optimistic search and BS, resulting in OBS (Algorithm~\ref{Alg:OBS}), with threshold $\gamma$. If $\gamma$ is chosen as
    $$
    C_2 \sqrt{\log T} \le \gamma \le C_3 \delta \lambda \sqrt{T}\qquad \text{with constants } C_2, \, C_3\,,
    $$
    then there is a constant $C_4$ such that
    $$
    \lim_{T \to \infty} \Prob{\hat \ncp = \ncp = 1 \text{ and } \abs{\hat\tau_1 - \tau} \le C_4\frac{\log T}{\delta^2 T}} = 1\,.
    $$
\end{enumerate}
\end{theorem}

In Theorem~\ref{th:noss}, part i.~states that the naive optimistic search detects the only change point with a localization error of order $\log T / (\delta^2 T)$, which is minimax optimal up to a possible log-factor (see e.g.~Lemma~2 in \citealp{WYR20}, and also Section~\ref{ss:theory}); Part ii.~states that treating the problem as multiple change points without knowing $\ncp = 1$, with a proper threshold, e.g.~$\gamma \asymp \sqrt{\log T}$, the OBS method, which will be discussed in details in Section~\ref{s:obs}, is able to estimate the number of change points correctly, and localize the only change point with the same error rate.
Condition~\eqref{eq:area} can be intuitively interpreted as a condition on the \emph{area} of the bump in the signal. In comparison with the weakest condition to ensure consistency of change point estimation, which is $\delta \sqrt{\lambda T} \gtrsim \sqrt{\log\log T}$ (see \citealp{Liu_Gao_Samworth_dyadic}), the naive optimistic search is sub-optimal.
However, in the particular case that $\lambda \asymp 1$, $\delta \sqrt{T}\gtrsim \sqrt{\log T}$, so the naive optimistic search is then nearly optimal in the sense that we lose only a log factor. This is the situation where the length of the left segment is comparable to that of the right segment, which we thus refer to as a ``balanced'' scenario. On the other hand, ``unbalanced'' scenarios are ones where the lengths of shorter and longer segments are very different. Theorem~\ref{th:noss} implies that naive optimistic search will not work well for very unbalanced scenarios, i.e., when the change point is very close to one of the boundaries. Further, we stress that the suboptimal condition \eqref{eq:area} cannot be improved for the naive optimistic search, see the next example. In this sense, the suboptimality is intrinsic to the naive optimistic search, rather than an artifact in our theoretical analysis. However, some improved versions of optimistic search lead to optimality as discussed in Theorem~\ref{th:aoss}.

\begin{example}
We consider a specific example of Model~\ref{Gaussian_setup} with $\ncp = 1$ and $\cpt_1 = \lambda \le 1/3$. For notational simplicity, let the step size be $\nu = 1/2$ in the naive optimistic search and $T/3 \in \N$. In the first step of naive optimistic search, we check the gain function at $T/3$ and $2T/3$. In order to avoid wrongly discarding $(0, T/3]$, we have to ensure $$\abs{\csm_{(0,T]}\Bigl(\frac{T}{3}\Bigr)} \ge \abs{\csm_{(0,T]}\Bigl(\frac{2T}{3}\Bigr)}\,.$$
In fact, we have
\begin{multline*}
    \Prob{\abs{\csm_{(0,T]}\Bigl(\frac{T}{3}\Bigr)} <\csm_{(0,T]}\Bigl(\frac{2T}{3}\Bigr)
    \le\Prob{\abs{\csm_{(0,T]}\Bigl(\frac{T}{3}\Bigr)} < \abs{\csm_{(0,T]}\Bigl(\frac{2T}{3}\Bigr)}} \\
    \le 2\,\Prob{\abs{\csm_{(0,T]}\Bigl(\frac{T}{3}\Bigr)} <\csm_{(0,T]}\Bigl(\frac{2T}{3}\Bigr)}\,.
\end{multline*}
Elementary calculation relying on properties of Gaussian distribution reveals
$$
\Prob{\abs{\csm_{(0,T]}\Bigl(\frac{T}{3}\Bigr)} <\csm_{(0,T]}\Bigl(\frac{2T}{3}\Bigr)} = \Phi\Bigl(-\delta \lambda\sqrt{\frac{T}{2}}\Bigr)\Phi\Bigl(\delta \lambda\sqrt{\frac{3T}{2}}\Bigr)\,,
$$
with $\Phi$ the distribution function of a standard Gaussian random variable. Thus, if and only if $\delta \lambda \sqrt{T} \to \infty$ it holds that
$$
\Prob{\abs{\csm_{(0,T]}\Bigl(\frac{T}{3}\Bigr)} \ge \abs{\csm_{(0,T]}\Bigl(\frac{2T}{3}\Bigr)}} \to 1 \qquad \text{as}\quad T \to \infty\,.
$$
Note that $\delta \lambda \sqrt{T} \to \infty$ is, up to a log factor, equivalent to the condition~\eqref{eq:area}, which guarantees that the probability of making a mistake in the first step of naive optimistic search vanishes eventually.
\end{example}

\subsection{Advanced and combined optimistic search}
In Algorithm~\ref{Alg:aOS} we propose \emph{advanced optimistic search} that improves on the naive version to tackle also unbalanced cases with the change point very close to the boundary. The main idea is to check a preliminary set of dyadic locations (up to rounding) to localize the change point approximately and then apply naive optimistic search in a suitable (balanced) neighborhood around the preliminary estimate in order to achieve a better localization. The preliminary estimate and the two locations marking its neighborhood are chosen from the dyadic locations, namely, as the location of the biggest gain as well as the closest dyadic neighbors thereof (to the left and to the right). This leads to a ``triangular structure'' for the gains of the three selected locations, where the probe point in the middle is guaranteed to have a gain at least as high as the boundaries. This essential property is afterwards maintained in the consecutive steps of naive optimistic search. Intuitively, as the dyadic points cover more points on the boundaries, the advanced search is suitable even in very unbalanced scenarios where the naive version fails. From a theoretical perspective, this modification even leads to minimax optimality as follows.

\begin{theorem}[Advanced/combined optimistic search]\label{th:aoss}
Under Model~\ref{Gaussian_setup} with a single change point, i.e.~$\ncp =1$, we assume that the minimal segment length $\lambda$ and the minimal jump size $\delta$ satisfy
\begin{equation}\label{eq:best}
\delta  \sqrt{\lambda T} \ge C_0 \sqrt{\log T}
\end{equation}
for some large enough constant $C_0$. Then:
\begin{enumerate}[i.]
    \item The advanced optimistic search (Algorithm~\ref{Alg:aOS}) and the combined optimistic search (Algorithm~\ref{Alg:combinedOS}) on $(0,T]$ provides an estimated change point $\hat\cpt = \hat s_{(0,T]}/T$ such that
    $$
    \lim_{T \to \infty}\Prob{\abs{\hat\cpt - \cpt} \le C_1 \frac{\log T}{\delta^2 T}} = 1 \qquad \text{with some constant } C_1\,.
    $$
    \item Let $\hat\cpt_i$, $i = 1,\ldots,\hat\ncp$, be the estimated change points (scaled by $1/T$) by the combination of advanced/combined optimistic search and BS, resulting in OBS (Algorithm~\ref{Alg:OBS}), with threshold $\gamma$. If $\gamma$ is chosen as
    $$
    C_2 \sqrt{\log T} \le \gamma \le C_3 \delta \sqrt{\lambda T}\qquad \text{with constants } C_2, \, C_3\,,
    $$
    then there is a constant $C_4$ such that
    $$
    \lim_{T \to \infty} \Prob{\hat \ncp = \ncp = 1 \text{ and } \abs{\hat\tau_1 - \tau} \le C_4\frac{\log T}{\delta^2 T}} = 1\,.
    $$
\end{enumerate}
\end{theorem}

Similar to the case of naive optimistic search (Theorem~\ref{th:noss}), it is shown in Theorem~\ref{th:aoss} that the advanced optimistic search (or the combined optimistic search introduced later) is able to localize the only change point at the best possible rate up to a possible log-factor, but now under a much weaker condition~\eqref{eq:best} instead. In comparison with the weakest condition $\delta \sqrt{\lambda T} \gtrsim \sqrt{\log \log T}$ by~\citet{Liu_Gao_Samworth_dyadic}, we lose here only a log factor, which seems to be possibly removable, see Section~\ref{ss:theory}. Therefore, the advanced optimistic search (or the combined optimistic search) possesses the (nearly) statistical minimax optimality like the full grid search, which checks every possible split point in $\{1, \ldots, T\}$. Note that optimistic searches only require $O(\log T)$ computations (which will be  explained later in detail), in sharp contrast to $O(T)$ required by the full grid search. It is a surprising fact that computational speed-ups come at almost no cost of statistical performance at all. That is,
``free lunch'' is possible!

In Algorithm~\ref{Alg:aOS}, checking the set of dyadic points $S$ of order $O(\log T)$ was inspired by \cite{Liu_Gao_Samworth_dyadic} and \cite{SeedBS}. The former authors used dyadic points in change point testing problems, and the seeded interval construction of the latter authors also relies on a similar idea. While seeded intervals are formulated more broadly for covering multiple change point cases, they can be adapted easily to the scenario of only a single change point, keeping $O(\log T)$ intervals. We also remark that the idea of a preliminary check of dyadic locations could be applied to other methods as well in order to obtain improvements. For example, the best split among dyadic locations could serve as initial split point candidate for the method of \cite{Kaul_2step}.

\begin{algorithm}[H]
\caption{advanced optimistic search}\label{Alg:aOS}
\begin{algorithmic}[1]
\Require $R-L > 2; ~L,R \in \mathbb{N}$ and step size $\nu \in (0,1)$ with $1/2$ by default
\Function{\textit{aOS}}{$\nu, L, R$}
\State $n \leftarrow \lfloor\log_2((R-L)/2)\rfloor$
\State $S \leftarrow \bigl\{ \floor{L+2^{-k}(R-L)},\,\ceil{R-2^{-k}(R-L)}\,:\, k =1, \ldots, n\bigr\}$ \Comment{Dyadic locations}
\State $s_* \leftarrow \argmax\limits_{s \in S} G_{(L,R]}(s)$ \Comment{Find the best split point on the ``dyadic'' grid}
\If{$s_* \le (R+L)/2$}
\State $l \gets \floor{s_* -(s_*-L)/2}$  and  $r \leftarrow \ceil{s_*+(s_*- L)}$ \Comment{``dyadic neighbours'' of $s_*$}
\Else
\State $l \leftarrow \floor{s_* - (R-s_*)}$  and $r \leftarrow \ceil{s_*+(R-s_*)/2}$ \Comment{``dyadic neighbours'' of $s_*$}
\EndIf
\State $\hat s_{(L,R]} \leftarrow$ \Call{\textit{nOS}}{$l, s_*, r \mid \nu, L, R$} \Comment{Naive optimistic search on $(l,r]$ containing $s_*$}
\State \textbf{return} $\hat s_{(L,R]}$
\EndFunction
\end{algorithmic}
\end{algorithm}

In practice, slightly altered versions of the advanced optimistic search might be equally viable e.g.~in certain high-dimensional settings and we present one in Appendix~\ref{Appendix_aos} that we also used for simulations. Furthermore, combining the results of naive and advanced optimistic search (Algorithm~\ref{Alg:combinedOS}), thus referred to as the \emph{combined optimistic search}, leads to slightly better empirical performance than the individual searches, but at a slightly higher computational cost. Also from a theory point of view, the combined optimistic search enjoys the same statistical minimax optimality as the advanced version.

\begin{algorithm}[H]
\caption{combined optimistic search}\label{Alg:combinedOS}
\begin{algorithmic}[1]
\Require $R-L > 2; ~L,R \in \mathbb{N}$ and step size $\nu \in (0,1)$ with $1/2$ by default
\Function{\textit{cOS}}{$\nu, L, R$}
\State $\hat s_a \leftarrow \aOS(\nu, L, R)$ \Comment{Advanced optimistic search}
\State $\hat s_n \leftarrow \nOS(L, \floor{(L + \nu R)/(1+\nu)}, R\mid\nu, L, R)$ \Comment{Naive optimistic search}
\If{$G_{(L,R]} (\hat s_a) \ge G_{(L, R]}(\hat s_n)$}
\State $\hat s_{(L,R]} \leftarrow \hat s_a$
\Else
\State $\hat s_{(L,R]} \leftarrow \hat s_n$
\EndIf
\State \textbf{return} $\hat s_{(L,R]}$
\EndFunction
\end{algorithmic}
\end{algorithm}

\begin{lemma} \label{Lem:cOS_speed}
Advanced and combined optimistic search (Algorithm~\ref{Alg:aOS} and Algorithm~\ref{Alg:combinedOS}) with a fixed step size $\nu=1/2$ terminate within $O(\log(R-L))$ and thus at most $O(\log T)$ steps (i.e., number of gain function evaluations).
\end{lemma}
The proof of Lemma~\ref{Lem:cOS_speed} is trivial, as the search over dyadic points in the advanced optimistic search (line 4 in Algorithm~\ref{Alg:aOS}) needs $O(\log (R-L))$ evaluations of the gain function (due to the logarithmic cardinality of the set~$S$), similar to the subsequent naive optimistic search (line 9 in Algorithm~\ref{Alg:aOS}, cf.~Lemma~\ref{Lem:nOS_speed}). Thus, overall, in the worst case one has $O(\log (R-L)) \le O(\log T)$ evaluations for the advanced optimistic search, and thus also for the combined optimistic search (which is the same as for the naive version).

We were primarily concerned with the number of evaluations of the gain function as that is the main driver of computational cost in recently emerging more complex, e.g.~high-dimensional and/or non-parametric problems mentioned in Section~\ref{Introduction}, or even scenarios where the underlying time series models are fitted via numerical optimization. Reducing the number of necessary fits using our optimistic search techniques could lead to methods that are computationally much more efficient in such applications. An example is given in Section~\ref{highdim_simulations}. Costly estimators such as the lasso, graphical lasso, random forests or time series methods used in complex scenarios often cannot be efficiently updated when adding one more observation. Hence, going through all possible split points in a full grid search requires basically new fits for each considered split point on the grid. This leads to severe computational issues for such estimators that we avoid using our optimistic search strategies. Note that the lack of efficient updates for such costly estimators is in contrast to many low-dimensional problems where efficient updates are typically possible. Computational gains of optimistic search compared to full grid search will hence depend on whether fits from neighboring segments can be re-used for efficient updates in the traditional full grid search based approach. Moreover, the final computational complexity with optimistic search depends primarily on the cost of single fits.

For the univariate Gaussian setting, the overall computational cost is only $O(\log T)$ if cumulative sums have been pre-computed and are freely available, as in that case each evaluation is possible in $O(1)$ time. Otherwise the $O(T)$ cost of calculating the cumulative sums becomes dominant. We remark that availability of cumulative sums (or similar ``sufficient statistics'' for the evaluation of the gain function) should be interpreted as a practical guidance when storing data for offline change point analysis, and is also a reasonable assumption in online change point detection. Given that online change point detection is essentially concerned with the detection of a single change point and advanced/combined optimistic search is capable of detecting change points even in very unbalanced scenarios, our technique could lead to fast online change point detection methodologies, see also Section~\ref{Discussion}.

\subsection{Scenarios beyond the univariate Gaussian setting}
Some technical modifications may be necessary for optimistic searches for scenarios beyond the univariate Gaussian setting. For instance, in high-dimensional models requiring a minimal number of observations, we have to exclude split points that are too close to the boundaries to ensure a sufficient number of observations for the fits (see e.g.~the estimator described in Appendix~\ref{Appendix_highdim} that we use in simulations in Section~\ref{highdim_simulations}). Advanced optimistic search is particularly easy to adapt: from the preliminary set of dyadic locations one can simply exclude those that are within a certain distance to the boundaries (see the version in Appendix~\ref{Appendix_aos}).

\section{Methodology and theory for multiple change points}
\label{Methodolgy_multiple}

We consider now the setup of multiple change points, and investigate how our methodology can be extended to such a more ambitious setup in order to still have a sublinear number of evaluations of the gain function and yet with theoretical guarantees for the estimation performance. Even binary segmentation (BS, \citealp{Vostrikova}), which is considered as a fast approximation to the optimal segmentation, performs a full grid search to find the best single split point over each of the considered intervals, resulting typically in $O(T\log T)$ evaluations. Another recently proposed fast method, Seeded Binary Segmentation \citep{SeedBS} needs $O(T\log T)$ evaluations in the worst case, again due to the full grid search. We aim for far fewer evaluations in the following.

\subsection{Optimistic Binary Segmentation (OBS)}\label{s:obs}
A very natural avenue is to speed up BS by performing fewer model fits in each step when searching for the new split point candidate using our optimistic search techniques. We call this approach \emph{optimistic binary segmentation} (OBS, see Algorithm~\ref{Alg:OBS}).

\begin{algorithm}[H]
\caption{Optimistic Binary Segmentation}\label{Alg:OBS}
\begin{algorithmic}[1]

\Require $0<\nu \leq 0.5; ~\gamma\in\mathbb{R}; ~m\in\mathbb{N}^+; ~L<R; ~L,R \in \mathbb{N}^+$
\Function{OBS}{$L, R \mid \nu, \gamma, m$}
\If {$R-L < m$} 	\Comment{Insufficient observations remaining}
\State STOP
\Else
\State $\hat s_{(L,R]}\gets$ split point returned by optimistic search on $(L,R]$ with step size~$\nu$.
\If {$g_{(L,R]}(\hat s_{(L,R]}) \ge \gamma$}														\Comment{Sufficient gain necessary}
\State Add $\hat s_{(L,R]}$ to the set of estimated change points.
\State \Call{OBS}{$L, \hat s_{(L,R]} \mid \nu, \gamma, m$}		\Comment{Further search on left segment}
\State \Call{OBS}{$\hat s_{(L,R]}, R \mid \nu, \gamma, m$}		\Comment{Further search on right segment}
\Else
\State STOP
\EndIf
\EndIf
\EndFunction
\end{algorithmic}
\end{algorithm}
OBS is started by the call OBS($1, T|\ldots$). Algorithm \ref{Alg:OBS} resembles BS, the difference is the way the ``best'' split point is approximated via the optimistic search (in line $5$) instead of the full grid search that is used in BS. Again, the choice of naive, advanced or combined optimistic search in line $5$ leads to the respective versions for OBS. The parameter $\gamma$~is a threshold defining how big the gains need to be in order to continue with further splitting, $m$~is a requirement for the minimal number of observations, which might be necessary in high-dimensional scenarios, and $\nu$~determines the step size for the chosen optimistic search variant (naive, advanced or combined). If the tuning parameters $\gamma$ and $m$ are set correctly, OBS needs $O(k_0  \log T)$ model fits ($k_0$ being the number of change points and typically much smaller than~$T$) compared to $O(T\log T)$ for BS. As the computational cost for OBS is adaptive to the number of change points, it is typically even faster than OSeedBS (introduced later in Section~\ref{Methodology_OSeedBS}). However, there are also some drawbacks. In some special cases with signal cancellation effects, optimistic search might fail in the presence of multiple change points even in the population case (see Example~\ref{Example_failure_OBS}) and general consistency results are thus not possible without further assumptions or slight modifications of the method (see Appendix~\ref{App:OBS_failure} for further details).

\begin{example}
\label{Example_failure_OBS}
Consider an example of Model~\ref{Gaussian_setup} with $\ncp=3$ change points at $\cpt_1 = 1/8,\, \cpt_2 =3/16,\, \cpt_3 = 1/4$ and the levels of segments $\mu_0 = 0,\, \mu_1 = 1,\, \mu_2 = -1,\, \mu_3 = 0$ (see case 3) in the right panel of Figure~\ref{f_star_plot} in the Appendix for a visualization). To ease notation, we set step size $\nu = 1/2$ for optimistic search and sample size $T = 2^n$ with $n\in\mathbb{N}$ and $n\ge3$. Note first that the population version of CUSUM statistics in~\eqref{eq:pgf} satisfies
$$
\abs{\csm_{(0,T]}^*(\lfloor T/3\rfloor)} = \abs{\csm_{(0,T]}^*(\lceil 2T/3\rceil)} = 0\,.
$$
This implies that the first step of naive optimistic search will make a mistake with probability $50\%$. Note further
$$
\abs{\csm_{(0,T]}^*(2^{-k}T)} = \abs{\csm_{(0,T]}^*(T-2^{-k}T)} = 0\quad \text{for }
k = 1, \ldots, n\,.
$$
Then the dyadic search in advanced optimistic search will end up at every dyadic location $\bigl\{2^{-k}T,(1-2^{-k})T \,:\, k = 1, \ldots,n\bigr\}$ with equal probability. Therefore, this example illustrates a scenario where all three optimistic searches might get stuck in a flat part of the population gain function and thus will fail with non-negligible probability.
\end{example}

As opposed to SeedBS, BS is not a minimax optimal method, hence, OBS would not be minimax optimal either, even if the issue with special cases similar to the one presented in Example~\ref{Example_failure_OBS} is solved. Moreover, model selection in BS and maybe even more so for OBS is more challenging compared to SeedBS. Hence, we advocate OBS as a very fast heuristic approach that performs reasonably in practice, in particular given its speed (see Section~\ref{highdim_simulations} for simulations). We present more details on OBS including explanations for the mentioned special cases and comparisons to BS in Appendix~\ref{Appendix_OBS}.

\subsection{Optimistic Seeded Binary
\label{Methodology_OSeedBS}
Segmentation (OSeedBS)} The main idea of seeded binary segmentation (SeedBS, \citealp{SeedBS}) is to apply searching for a single change point in various intervals with the hope that some of these intervals contain only a single change point, where the detection is ``easy''. While the best split point in each interval is a candidate, the decision which candidates to declare finally as a change points depends on a subsequent selection step. The intervals are called seeded intervals and they are constructed deterministically (see Definition~\ref{def:seeded_intervals}). SeedBS is thus very similar to wild binary segmentation (WBS, \citealp{Fryzlewicz_WBS}) and the narrowest over threshold method (NOT, \citealp{Baranowski}). The latter two procedures use random intervals instead of the deterministic ones as intervals which in general leads to total length and number of considered intervals to be larger and thus computationally more expensive.

\begin{definition}[Seeded intervals of \citealp{SeedBS}]
\label{def:seeded_intervals}
Let $a \in [1/2,1)$ denote a given decay parameter. Let $I_1=(0,T]$. For $k = 2,\ldots, \lceil\log_{1/a}(T)\rceil$ (i.e.~logarithm with base $1/a$) define the $k$-th layer $\mathcal I_k$ as the collection of $n_k$ intervals of initial length $l_k$ that are evenly shifted by the deterministic shift $s_k$:
\begin{equation*}
    \mathcal I_k = \bigcup\limits_{i=1}^{n_k} \{ (\lfloor (i-1)s_k \rfloor,
        \lceil(i-1)s_k + l_k \rceil ] \},
\end{equation*}
where
$n_k = 2\lceil (1/a)^{k-1} \rceil - 1$,
$l_k  = Ta^{k-1}$ and
$s_k  = (T-l_k)/(n_k-1)$.
The overall collection of seeded intervals is
\begin{equation*}
    \mathcal I = \bigcup\limits_{k=1}^{\lceil\log_{1/a}(T)\rceil} \mathcal I_k.
\end{equation*}
\end{definition}

Note that $\mathcal I$ contains $O(T)$ intervals, as it is constructed to guarantee appropriate coverage even in frequent change point scenarios with up to $O(T)$ change points. On the other hand, when assuming that there is only a single change point, all intervals that do not have a starting point at $1$ or do not have an end point at $T$ can be discarded, reducing the number of intervals to $O(\log T)$. In the case of multiple change points, assuming a certain minimal spacing between change points also allows to discard intervals that are too short. Note that $\abs{\mathcal I_k}$ (the number of intervals on the $k$-th layer) increases exponentially and thus skipping the lowest layers with a length below a pre-specified minimal segment length $m$ reduces the number of intervals drastically.

\begin{algorithm}[H]
\caption{Optimistic Seeded Binary Segmentation}\label{Alg:OSeedBS}
\begin{algorithmic}[1]
\Require a decay parameter $a\in[1/2,1)$, a minimal segment length $m \ge 2$ as well as tuning parameters for the selected optimistic search (advanced or combined)
\Function{OSeedBS}{}
\State $\mathcal I \gets$  seeded intervals (Definition \ref{def:seeded_intervals}) with decay $a$ and at least $m$ observations.
\For{$(L,R] \in \mathcal I$}
\State \qquad $\hat s_{(L,R]}\gets$ the split point returned by optimistic search on $(L,R]$.
\EndFor
\State Apply some selection method to $\bigl(\hat s_{(L,R]},G_{(L,R]}(\hat s_{(L,R]})\bigr), (L,R] \in \mathcal I$, to output the final change point estimates.
\EndFunction
\end{algorithmic}
\end{algorithm}

The difference between SeedBS and its optimistic counterpart OSeedBS is essentially in line $4$ of Algorithm~\ref{Alg:OSeedBS}, where we perform optimistic search rather than full grid search. The choice of naive, advanced or combined optimistic search in line $4$ leads to the respective versions for OSeedBS. The selection method in line $5$ can be for example greedy or narrowest over threshold (NOT) based selection, see \cite{SeedBS} for details. The computational times of OSeedBS depend of course on the chosen optimistic search type, but even more critical is the choice of minimal segment length $m$. If $m=O(T)$, only a handful of intervals are considered, with $O(\log T)$ evaluations each, and thus $O(\log T)$ evaluations overall. For the other extreme, when $m$ is very small, many intervals need to be generated and thus the main driver of the number of evaluations (and hence, the computational cost) is the number of considered intervals. For $m=2$, the number of intervals and also the total number of evaluations is $O(T)$. The estimation performance of course also depends on the choice of $m$. If chosen too small, estimation performance will be bad as change points within short segments cannot be detected. Thus, intuitively, $m$ offers some kind of trade-off between estimation performance and computational efforts. Note that such trade-offs are inherent also in other methods, see e.g.~the number of random intervals chosen in WBS.

\begin{theorem}\label{th:osbs}
Under Model~\ref{Gaussian_setup}, we assume that the minimal segment length $\lambda$ and the minimal jump size $\delta$ satisfy the same condition~\eqref{eq:best} as in Theorem~\ref{th:aoss}. Assume that there is an a priori known lower bound $\lambda_*$ of all segment lengths, i.e.
\begin{equation}\label{eq:mlen}
\lambda \ge \lambda_* \asymp T^{-\omega} \qquad \text{for some constant } \omega \in [0,\,1]\,.
\end{equation}
By $\hat\ncp$ and $\hat\cpt_1 < \cdots < \hat\cpt_{\hat \ncp}$ denote respectively the number and the locations of estimated change points by OSeedBS (Algorithm~\ref{Alg:OSeedBS}), where the NOT is chosen as the selection method, and the seeded intervals are constructed with decay $a \in [1/2,1)$ and minimal length $m = \lfloor \lambda_* T/3\rfloor$. Then:
\begin{enumerate}[i.]
\item
There are generic constants $C_1$, $C_2$, independent of $T$, $\omega$, $a$ and $m$, such that, given the threshold for the selection method
$\gamma = C_1 \sqrt{\log T}$,
$$
\lim_{T \to \infty}\Prob{\hat \ncp = \ncp,\,\, \max_{i = 1, \ldots, \ncp} \delta_i^2\abs{\hat \cpt_i - \cpt_i} \le C_2 \frac{\log T}{T}} = 1\,.
$$
\item
The number of evaluations is $O(\min\{T^\omega \log T,\, T\})$.
\end{enumerate}
\end{theorem}

We emphasize that assumption~\eqref{eq:mlen} is mainly needed for computational purpose, which ensures a sublinear number of evaluations as specified in part ii.~of Theorem~\ref{th:osbs}. If cumulative sums have been pre-computed and are freely available, then the overall computational cost itself is also $O(\min\{T^\omega \log T,\, T\})$, i.e., it equals the number of evaluations. However, if cumulative sums are not available, then the $O(T)$ cost of calculating cumulative sums becomes dominant and the overall computational cost is $O(T)$, see Section~\ref{ss:mcpt} for further details.
For the statistical guarantee in part i., an assumption on the minimal spacing, i.e.~\eqref{eq:mlen}, becomes obvious if we choose $\omega = 1$ and $m =2$, since it is often pointless to work on a higher resolution than the sampling rate $1/T$ without further model assumption, and thus in this sense it imposes no restriction at all. In case of multiple change points, the signal strength condition \eqref{eq:best} is the weakest one that still allows for detection. This coincides with the best known results (e.g.~\citealp{Frick_etal}; \citealp{Baranowski}; \citealp{Chan_Chen2017} and \citealp{WYR20}, to name a few) with only difference in multiplying constants.
Following the proof in Appendix~\ref{ss:mcpt}, we can easily replace it by
$$
\min_{i=1,\ldots,\ncp} \Bigl(\min\{\cpt_{i+1}-\cpt_i,\,\cpt_i - \cpt_{i-1}\}\delta_i^2 T\Bigr)\,\, \gtrsim \,\,\log T\qquad \text{as} \quad T \to \infty.
$$
This is slightly more general, as it allows for frequent large jumps and meanwhile small jumps over long segments (cf.~\citealp{ChKi19}). However, we prefer the current version as in Theorem~\ref{th:osbs}, for notational simplicity. Note, moreover, that the localization rate reported in part i.~of Theorem~\ref{th:osbs} is minimax optimal up to a possible log factor.

WBS, and also the similar NOT method, can be speeded up as well using optimistic search. However, in the worst case with very short spacings between change points, e.g.~in frequent change point scenarios with up to $O(T)$ change points, these two methods need to draw up to $O(T^2)$ random intervals, which prohibits sublinear number of evaluations overall. Nonetheless, we expect substantial computational gains using optimistic search in connection with many other multiple change point detection techniques compared to the respective full grid search based counterparts.

\section{Simulations}
\label{Simulations}
In order to illustrate the search methods and to examine the corresponding theory, we consider simulation settings first in the univariate change in mean case for a single and then for multiple change points. Furthermore, we illustrate speedups in a high-dimensional setup.

\subsection{Single change point in the mean of univariate Gaussian observations}
\begin{example} \label{setup:single_mean_shift}
Let $X_{1}, \ldots, X_{100}\sim\mathcal{N}(0,\sigma^2)$ and $X_{101}, \ldots, X_{100+n}\sim\mathcal{N}(0.5,\sigma^2)$ be independent observations with a single change in the mean value at observation $100$.
\end{example}
The top part of Table~\ref{tab:simulation_results} shows the localization error of the change point estimates found by the naive, advanced and combined optimistic search, as well as the full grid search for various choices of $n$ (from 100 to 5,000) for three different noise levels ($\sigma = 0.5, 1, 1.5)$. In addition, Figure~\ref{Fig:single_cpt} in Appendix~\ref{Appendix_Simulations} shows found change points using various search methods in 1,000 simulations for a ``balanced'' ($n=200$) and an ``unbalanced'' ($n=\text{5,000}$) scenario. Recall that we call scenarios where the length of the shorter versus the longer segment is not too different as balanced, while scenarios where the change point is relatively close to the boundaries as unbalanced. While the naive optimistic search clearly struggles when the length of the two segments is very unbalanced, i.e., when $n$ is large, in particular when the noise level is high, the advanced optimistic search has a stable performance even in such unbalanced scenarios. However, for the more balanced scenarios (up to $n=400$ for $\sigma=1.5$ for example), the naive version can be better than the advanced one (and even better than the full grid search for e.g.~$n=200$). Given that the naive optimistic search has a preference towards points in the middle, and in the balanced scenarios the true change point lies in there, this should not come as a surprise. The combined optimistic search, as expected, has a slightly improved performance compared to the advanced search (in particular for the rather balanced scenarios). The full grid search has the best performance for the rather challenging scenarios that are very unbalanced and/or have a high noise level, but the advanced and combined optimistic search strategies come remarkably close. We note that increasing absolute errors in change point location for higher values of $n$ despite having more available observations is actually reasonable, as there are meanwhile more potential candidates for change points on the grid. This is also compatible with the theoretical bound $\asymp \log T /\delta^2$.  Overall, the simulation results in Table~\ref{tab:simulation_results} confirm our theoretical results that the naive optimistic search is not consistent for unbalanced signals, while the advanced and combined versions are. Moreover, the price to pay in terms of localization error for all noise levels seems to be fairly low. In terms of computation, the bottom part of Table~\ref{tab:simulation_results} shows that the number of evaluations for optimistic search variants can be orders of magnitude smaller compared to full grid search, in particular if $n$ is large.

\begin{table}
\caption{\label{tab:simulation_results} Simulation results for Example~\ref{setup:single_mean_shift} for various choices of noise level $\sigma$ and number of observations $n$ from the second segment. Reported are average absolute differences between the true change point at location $100$ and the location of the best single split point found by the respective search method (top three blocks) as well as the average number of evaluations (bottom block). Values are averaged over $10000$ simulations (rounded to two digits) and in parentheses the corresponding standard deviations (rounded to integers). The average number of evaluations for noise levels $\sigma=0.5$ and $\sigma=1.5$ are not reported as they are very similar to the case $\sigma=1$.}
\centering
\fbox{\begin{tabular}{llcccc}
&   &   &   &   & \\
&  & \multicolumn{4}{c}{average absolute estimation error for search methods}\\
\hline
\em noise level& $n$& \em naive        & \em advanced     & \em combined       & \em full search\\
\hline
&               100  &   3.38 (7)   & 2.77 (4)     & 2.88 (5)       &  3.24 (5) \\
&               200  &   2.72 (4)   & 4.22 (7)     & 2.95 (5)       &  3.17 (5) \\
&               300  &   3.43 (7)   & 4.45 (8)     & 3.21 (5)       &  3.16 (5) \\
$\sigma=0.5$&   400  &   4.68 (10)  & 3.95 (6)     & 3.37 (5)       &  3.16 (5) \\
&               500  &   6.55 (27)  & 4.24 (8)     & 3.09 (5)       &  3.08 (5) \\
&               1000 &  13.75 (74)  & 3.84 (6)     & 3.35 (5)       &  3.08 (5) \\
&               2000 & 171.74 (387) & 3.92 (7)     & 3.26 (6)       &  3.01 (4) \\
&               5000 &1021.12 (1338)& 3.92 (7)     & 3.52 (6)       &  3.05 (5) \\
\hline
&               100  &  15.86   (20) &  15.26  (23) &   15.07  (21)  &     16.79  (22) \\
&               200  &  12.37   (18) &  28.93  (43) &   15.78  (26)  &     17.44  (28) \\
&               300  &  19.50   (34) &  26.91  (45) &   19.30  (35)  &     17.73  (33) \\
$\sigma=1$&     400  &  30.58   (56) &  26.02  (54) &   20.14  (42)  &     17.85  (37) \\
&               500  &  50.09   (87) &  26.97  (59) &   21.06  (49)  &     18.80  (44) \\
&               1000 & 136.75  (240) &  29.70  (94) &   24.59  (81)  &     21.24  (72) \\
&               2000 & 544.70  (547) &  35.73 (160) &   34.16 (156)  &     24.21 (116) \\
&               5000 &1948.79 (1328) &  48.08 (341) &   51.94 (354)  &     38.34 (298) \\
\hline
&               100  &  25.24   (25)  &   33.95   (35) &  31.70    (32)   &   34.19    (33) \\
&               200  &  23.77   (29)  &   60.82   (62) &  39.03    (50)   &   42.05    (52) \\
&               300  &  41.23   (54)  &   65.17   (82) &  50.79    (72)   &   48.55    (72) \\
$\sigma=1.5$&   400  &  62.98   (85)  &   70.69  (107) &  58.85    (95)   &   56.11    (93) \\
&               500  &  96.54  (114)  &   82.27  (134) &  70.03   (121)   &   62.41   (115) \\
&               1000 & 253.11  (291)  &  121.14  (256) &  114.73  (243)   &   98.52   (226) \\
&               2000 & 739.92  (534)  &  202.01  (504) &  203.74  (493)   &   156.51  (434) \\
&               5000 &2171.28 (1211)  &  436.96 (1269) &  455.99 (1260)   &   355.35 (1123) \\
\hline
\hline
&   &   &   &   & \\
&  & \multicolumn{4}{c}{average number of evaluations for search methods}\\
\hline
\em noise level& $n$& \em naive        & \em advanced     & \em combined       & \em full search\\
\hline
&            100 & 16.18 (1)  &  25.10 (1)  &  41.28 (2)    &     199 (0) \\
&            200 & 17.31 (1)  &  25.92 (2)  &  43.24 (2)    &     299 (0) \\
$\sigma=1$&  500 & 19.08 (1)  &  29.34 (2)  &  48.43 (2)    &     599 (0) \\
&           1000 & 19.36 (1)  &  30.95 (1)  &  50.31 (2)    &    1099 (0) \\
&           2000 & 21.37 (1)  &  33.00 (1)  &  54.36 (2)    &    2099 (0) \\
&           5000 & 23.69 (1)  &  35.02 (1)  &  58.71 (2)    &    5099 (0) \\
\end{tabular}}
\end{table}

\subsection{Multiple change points in the mean of univariate Gaussian observations}
We illustrate the performance in a univariate change in mean scenario, but with multiple change points. Example~\ref{setup:blocks_signal} describes the blocks signal of \citet{DoJo94} with the noise level as used by \cite{Fryzlewicz_WBS} in his simulation study.

\begin{example} \label{setup:blocks_signal}
Consider a total of $2048$ observations with $11$ change points at locations $205, 267$, $308, 472$, $512, 820$, $902, 1332$, $1557, 1598$ and $1659$ as well as mean values $0, 14.64$, $-3.66, 7.32$, $-7.32, 10.98$, $-4.39, 3.29$, $19.03$, $7.68, 15.37$ and $0$ between the change points to which independent Gaussian noise with a standard deviation of $\sigma = 10$ is added.
\end{example}

Most of the points in the left panel of Figure~\ref{Fig:blocks_signal} line up on the diagonal indicating that the segmentation found by full grid search in the seeded intervals (SeedBS) and combined optimistic search in the seeded intervals (OSeedBS) have exactly the same localization error in terms of Hausdorff distance. This holds for the whole range of considered minimal segment lengths of $2$ to $128$ observations. In the right panel of Figure~\ref{Fig:blocks_signal} the points still line up around the diagonal, but somewhat more spread out. Hence, with the naive optimistic search the segmentation found overall, differs more often. Nonetheless, the average Hausdorff distances across the $100$ simulations and various minimal segment length constraints (as indicated by the vertical, respectively horizontal dashed lines) are still very close for the full grid search and the naive optimistic search. Hence, the average performance of both naive and combined optimistic search is very close to that of the full search. Instead of the search method, what seems to matter more for the overall localization error in terms of Hausdorff distance is the minimal segment length constraint. The best minimal segment length constraints from the investigated ones seem to be $32$ and~$64$ as indicated by the averages. This is not entirely surprising as the five shortest distances between two consecutive change points in the blocks signal are $40, 41, 41, 61$ and~$62$ observations. When the minimal segment length constraint is smaller, e.g.~only~2 observations, then the chance of including a false positive coming from a very short segment increases and hence, Hausdorff distances can be marginally worse. On the other hand, when the minimal segment length constraint for the seeded intervals is $128$, one might miss some of the true change points, resulting in a clearly worse average performance. Overall, it turns out that the optimistic variants of OSeedBS have a competitive average performance compared to the full grid search based SeedBS. Moreover, as long as the minimal segment length constraints are short enough to guarantee coverage of each single change points, both SeedBS and variants of OSeedBS perform well.

\begin{figure}
\centering
\includegraphics[width=0.8\textwidth]
{figures/fig3.pdf} \\
\caption{Pairwise plots of Hausdorff distances of the location of the best $11$ change point candidates (with greedy selection) compared to the true ones in $100$ simulations for SeedBS (decay $a=1/\sqrt{2}$) with various minimal segment length constraints and full grid search in each seeded interval (horizontal-axis) versus combined optimistic search (vertical-axis, left) and naive optimistic search (vertical-axis, right). The vertical and horizontal dashed lines indicate the average Hausdorff distances over the $100$ simulations for the respective minimal segment length constraint and search method within the seeded intervals. Note the logarithmic scales on both axes.}
\label{Fig:blocks_signal}
\end{figure}

\subsection{A high-dimensional example}
\label{highdim_simulations}
First, we introduce a changing covariance setup in Model~\ref{GGM_setup}. Specific instances used in simulations follow in Example~\ref{setup:high-dim}. Lemma~\ref{Lem:GGM_convexity} below reveals that the population gain function for this model shares the piecewise (in between change points) convex structure similar to the univariate Gaussian change in mean setup (see Appendix~\ref{Appendix_OBS}). The piecewise (quasi)convex structure is at the core of optimistic search based change point detection and thus motivates the application of optimistic search also for Model~\ref{GGM_setup}.

\begin{Model}
\label{GGM_setup}
Assume that observations $X_1, \ldots, X_T$ are independent and
\begin{equation*}
\begin{split}
 X_{\cpt_{0}T+1}(=X_1), \ldots, X_{\cpt_1 T} & \sim \mathcal{N}(0, \Sigma_0)\,, \\
& \vdots	\\
 X_{\cpt_{\ncp}T+1}, \ldots, X_{\cpt_{\ncp+1}T}(=X_{T}) 	& \sim \mathcal{N}(0, \Sigma_\ncp)\,,
\end{split}
\end{equation*}
where $\{\cpt_i \,:\, i = 1, \ldots, \ncp\}$ gives the location of change points satisfying
$$
0 = \cpt_0 < \tau_1 < \cdots < \cpt_{\ncp+1} = 1\quad \text{and}\quad \cpt_i T \in \N\,.
$$
The means are $0\in\mathbb{R}^p$ while $\Sigma_i \neq \Sigma_{i-1}\in\mathbb{R}^{p\times p}$ for $i=1,\ldots,\ncp$ give the covariances on the segments.
\end{Model}
We first introduce $\Sigma_{(l,r]}$ as the convex combination of the covariance matrices within the segment $(l,r]\subseteq(0,1]$ with the weights given by the relative segment lengths within $(l,r]$, that is,
$$
\Sigma_{(l,r]} = \sum_{i = 0}^\ncp \frac{\abs{(\cpt_i,\,\cpt_{i+1}]\cap (l,\,r]}}{\abs{(l,\,r]}}\, \Sigma_i\,,
$$
where $\abs{(a,\,b]} = b-a$ is the length of an interval $(a,\,b]$. We consider the following function
\begin{equation}\label{e:gainHD}
    G^*_{(l,r]}(s)
    =
    \frac{r-l}{T}\log(|\Sigma_{(l,r]}|)
    -\frac{s-l}{T}\log(|\Sigma_{(l,s]}|)
    -\frac{r-s}{T}\log(|\Sigma_{(s,r]}|)\,,
\end{equation}
where $\abs{A}$ denotes the determinant of a square matrix $A$\,. Note that $G^*_{(l,r]}(\cdot)$ in \eqref{e:gainHD} can be thought of as the continuously embedded population gain function corresponding to multivariate Gaussian log-likelihood based segmentation. That is, we replace  the empirical covariance matrices by the corresponding true ones in the generalized likelihood ratio test for a single change point at location $s$, $l < s < r$,  against a constant signal.

\begin{lemma} \label{Lem:GGM_convexity}
The function $G^{*}_{(l,r]}(\cdot)$ defined in \eqref{e:gainHD} is piecewise (i.e.~in between change points) convex, and up to some special cases as detailed in the proof (see Appendix~\ref{App:GGM_convexity_proof}) even strictly convex.
\end{lemma}
Lemma~\ref{Lem:GGM_convexity} in particular implies that in the presence of a single change point in $(l,r]$, $G^*_{(l,r]}$ is unimodal and in case $(l,r]$ contains multiple change points, then each strict local maximum corresponds to a change point. In case that $G^{*}$ is even piecewise strictly convex, then Lemma~\ref{Lem:OS} from the Appendix implies that optimistic search would find a strict local maximum and hence a change point. In practice, we do not have access to the population gain function $G^{*}$, but only a noise variant thereof. Given this motivation, in the following we focus on how well optimistic search based change point detection methods perform in the changing covariance setup of Model~\ref{GGM_setup}.

\cite{SeedBS} considered a simulation setup involving change point detection in (high-dimensional) Gaussian graphical models as an example to demonstrate that seeded intervals can be computationally much more attractive compared to the random intervals utilized in WBS in high-dimensional settings with computationally expensive fits. We will show that further speedups in such computationally challenging change point detection problems for many available algorithms can be easily obtained utilizing our optimistic search strategies. In Example~\ref{setup:high-dim} we take up on this setup by \cite{SeedBS} and add a similar, but somewhat simpler scenario.

\begin{example}
\label{setup:high-dim}
Let $\Sigma_{ij} = \exp{(-\frac{1}{2} \abs{s_i-s_j})}$ with $s_i - s_{i-1} = 0.75, i = 2,\ldots,20$ be a chain network model (see e.g.~Example~4.1 in \citealp{fan2009_chain_network}) with $p=20$ variables. A modified version $\tilde \Sigma$ is obtained by replacing the top left $5\times5$ block of $\Sigma$ by a $5$-dimensional identity matrix.
\begin{itemize}
    \item[a)] In the setup from \cite{SeedBS} we set in Model~\ref{GGM_setup} $\Sigma_0=\Sigma, \Sigma_1=\tilde\Sigma, \Sigma_2=\Sigma, \Sigma_3=\tilde\Sigma$, etc.~and draw $100$ observations for each segment until obtaining a total of $T=\text{2,000}$ observations. Hence, there are $20$ segments of length $100$ each, with a total of $\ncp=19$ change points.
    \item[b)] In Model~\ref{GGM_setup} we set $\Sigma_0=\Sigma_2=\Sigma_4=\Sigma, \Sigma_1=\Sigma_3=\Sigma_5=\tilde\Sigma$, and draw $550, 300, 700, 250, 100$ and $100$ observations for the respective segments, obtaining again a total of $T=\text{2,000}$ observations, but this time with $6$ segments, i.e., a total of $\ncp=5$ change points.
\end{itemize}
\end{example}

We consider change point detection where the gain function is based on the multivariate Gaussian log-likelihood, where the underlying precision matrices are obtained using the graphical lasso precision matrix estimator of \citet{glasso}. More details are in Appendix~\ref{Appendix_highdim}, see equation~\eqref{glasso_change_est}. The graphical lasso is rather costly especially when repeatedly fitting at each possible split point $s$ on a grid. Note that the essential problem is that the estimator $\hat\Omega_{(u,s]}^{\mathrm{glasso}}$ for a segment $(u, s]$ cannot be efficiently updated (not even using warm starts) to obtain $\hat\Omega_{(u,s+1]}^{\mathrm{glasso}}$ for the segment $(u,s+1]$. Hence, the overall number of graphical lasso fits is the main driver of computational time. Note that Example~\ref{setup:high-dim} is actually not a truly high-dimensional scenario, but for higher dimensions~$p$ the computational times went beyond reasonable for the full grid search based approaches that we wanted to include as baselines.

Our goal is to compare the estimation performance and computational times of three full grid search based approaches (BS, SeedBS and WBS) and the three optimistic search based counterparts (OBS, OSeedBS and OWBS) when using the computationally intense graphical lasso based estimator described in Appendix~\ref{Appendix_highdim}. In the examples below we use advanced optimistic search for OBS, OSeedBS and OWBS.

We simulated $10$ different data sets from Example~\ref{setup:high-dim}. Setup~a) is shown on the top of Figure~\ref{Fig:glasso_multiple_changes} while setup~b) is shown in the middle. For SeedBS and OSeedBS five different decays for the seeded intervals, and for WBS as well as OWBS five values $M$ for the number of random intervals were investigated, leading to $5$ respective point clouds for each of these four methods in the top and middle plots of Figure~\ref{Fig:glasso_multiple_changes}. In order to eliminate the effect of model selection, for all algorithms we selected greedily as many change points as the true underlying number (with some exceptions for WBS with small~$M$).

\begin{figure}
\centering
\includegraphics[width=0.9\textwidth]
{figures/gl1.pdf} \\
\includegraphics[width=0.9\textwidth]
{figures/gl2.pdf} \\
\vspace{-0.1cm}
\includegraphics[width=0.9\textwidth]
{figures/gl3.pdf} \\
\vspace{-0.2cm}
\caption{Estimation performance (in terms of Hausdorff distance) and computational times on the changing Gaussian graphical model from Example~\ref{setup:high-dim} with setup~a) on the top and setup~b) in the middle for $10$ simulated data sets each. The symbols differentiate between the basic algorithms and the colors indicate whether full grid search or optimistic search was used (``Optimistic'' in the legend indicating the use of the optimistic variant). The five point clouds for SeedBS and OSeedBS correspond to decay parameters $a=1/2, 1/2^{1/2}, 1/2^{1/4}, 1/2^{1/8}, 1/2^{1/16}$ for the underlying seeded intervals, while the five point clouds for WBS and OWBS correspond to $M = 100, 200, 400, 800, 1600$ random intervals. For all $3\times2$ algorithms we picked the segmentation corresponding to the true number of change points from the solution path with some exceptions for WBS with lower values of $M$ for setup~a) where the solution path contained less than $19$ candidates. Here we took all proposed change points from the solution path. The plot on the bottom shows Hausdorff distances (on a logarithmic scale) of the solution paths for OBS (left) and BS (right) for the $10$ different simulated data sets from setup~a) shown in the top panel. The vertical dashed line at $19$ indicates the true number of change points in setup~a). The colors on the bottom plot are only used as a visual aid to easier distinguish the paths, and are unrelated to the color legend of the upper panels. The colored dots indicate the segmentation where the best Hausdorff distance from the respective solution path was first attained.}
\label{Fig:glasso_multiple_changes}
\end{figure}

The top and middle panels of  Figure~\ref{Fig:glasso_multiple_changes} show roughly speedups of factor~$30$ for OBS compared to BS, factor~$35$ for OWBS versus WBS and factor~10--14 for OSeedBS versus SeedBS in both of the considered setups. The computational times for SeedBS, OSeedBS, WBS and OWBS across setups~a) and~b) of Example~\ref{setup:high-dim} are very similar due to the fact that both setups have $T=\text{2,000}$ observations and we used the same type of seeded respective random interval settings for the two examples. BS and OBS on the other hand are adaptive to the number of change points (if the stopping criterion is set appropriately) and thus considerably faster for setup~b) which only has $5$ change points compared to the $19$ for setup~a). We provide a further discussion on the speedups of different approaches and potential benefits of combining them in Appendix~\ref{Appendix_computational_gains_highdim}.

The disadvantage of WBS (and OWBS) is that a large number of random intervals might be needed to obtain good estimation performance. Compared to that, for SeedBS (and OSeedBS), the decay parameter has much less influence on the estimation performance and the deterministic seeded interval construction seems computationally much more attractive than the random ones. BS and OBS do very well on setup~b), only slightly worse than e.g.~SeedBS, but on the more challenging setup~a), BS and OBS are clearly worse than results obtained with SeedBS (and OSeedBS). Note also that the difficulty of model selection does not appear in the estimation performances shown here, as we took as many estimated change points as the true number. In general, model selection, i.e. setting appropriate thresholds for the detection is much more difficult for BS (or OBS) than for SeedBS (or OSeedBS) and thus, in challenging scenarios BS (or OBS) is expected to perform worse than SeedBS (or OSeedBS) when incorporating the model selection part as well. For easier scenarios, e.g.~setup~b) in the middle of Figure~\ref{Fig:glasso_multiple_changes}, BS (or OBS) is not much worse than SeedBS (or OSeedBS).

As expected, we see slightly worse estimation performance for the optimistic variants compared to the full grid search. Bigger differences only occur in setup~a) for OBS compared to BS (top panel of Figure~\ref{Fig:glasso_multiple_changes}). The segmentation found by OBS in most of the cases has a Hausdorff distance of around either $100$ or $200$ (indicating that 1--2 change points were missed) compared to around $100$ or below for BS. The colored dots in bottom plot of Figure~\ref{Fig:glasso_multiple_changes} indicating the segmentation with best Hausdorff distance from the solution path (i.e. for varying number of estimated change points) should ideally lie on the vertical dashed lines. However, both for BS and OBS it seems like occasionally false positives are picked prior to the candidates corresponding to true change points (leading to over-segmentation in a sense) as one often has to include a few more change point candidates in order to reach the best possible Hausdorff distance attainable with the given solution path. This problem is more pronounced for OBS, explaining the worse performance of OBS in the top plot of Figure~\ref{Fig:glasso_multiple_changes} showing the segmentations according to the true number of change points. While the quality of estimated change point is not bad, OBS increasingly faces the difficulty to select the candidates corresponding to true change points rather than picking false positives. In other words, the difficult model selection problem for BS in challenging scenarios can be even further amplified for OBS.

Overall, we demonstrated that optimistic search is a very generic technique leading to massive computational speedups in complex high-dimensional change point detection problems for many existing change point detection algorithms, while not sacrificing much in terms of estimation performance.

\section{Discussion}
\label{Discussion}
We introduced optimistic search strategies that avoid the full grid search and thus lead to computationally fast change point detection methods applicable in various change point detection problems. For the classical univariate Gaussian changing mean setup we proved that the advanced and combined optimistic search are asymptotically minimax optimal for detecting a single change point with only a logarithmic number of evaluations of the gain function. For multiple change point problems we combined optimistic search with seeded binary segmentation, leading to asymptotically minimax optimal detection while having superior runtime compared to existing approaches. Our methodology is most relevant for complex change point detection problems with computationally expensive model fits, as demonstrated by the massive computational gains in examples involving high-dimensional Gaussian graphical models.

\subsection{Further applications of optimistic search}
We see various options to extend and refine our methodology. First and foremost, online change point detection is essentially concerned with the detection of a single change point with only a handful of observations after the change happened. In such unbalanced scenarios, advanced and combined optimistic search could be particularly helpful. Besides the advantage of only a logarithmic number of evaluations, a logarithmic runtime of the entire corresponding algorithm requires the availability of cumulative sums or similar ``sufficient statistics''. The latter assumption is realistic in online detection, as one can continuously update the necessary statistic (e.g.~cumulative sums). Applications could be e.g. in volatility changes for financial time series \citep{Spo09}. We refer to \citet{Basseville_Nikiforov_book} and \citet{Tartakovsky_etal_book} for an overview and many additional references on online change point detection.

We mostly focused on SeedBS as a tool to extend our optimistic search methods to multiple change point detection. However, one could speed up any other algorithm which relies on grid search, including WBS, NOT, or the Circular Binary Segmentation (CBS) algorithm of \citealp{CBS} which is popular in genomic applications involving copy number variations. If evaluation of the gain function (without the model fitting) is considerably cheaper than computing the model fits (which is typical in high-dimensional scenarios), then optimistic search can also be used for finding a candidate for the initial split point in the two step procedure of \cite{Kaul_2step} to obtain further refinements and possible improvements of their estimate. These authors perform a full grid search for the comparably cheap evaluation of the gain function relying on a good initial split point and corresponding fixed fit for that specific split point. Our optimistic search would propose initial split points that are close to the true change point and thus considerably better than the initial split $T/2$ used in \cite{Kaul_2step} in the absence of prior information on the location of the change point, or similar non-informative initial split points obtained by random guessing.

Another prime example for grid search are tree algorithms \citep{CART}: here, one can replace the underlying grid search by our optimistic search, at least for continuous covariates. Instead of the usual reduction in mean squared error criterion, our search might allow for more complex split criteria for which full grid search would be too expensive. Trees have direct implications for Random Forests \citep{Breiman2001}, possibly opening another range of applications both for our procedures for single and multiple change point detection. Admittedly, these scenarios could be farther away from the idealized piecewise convex setting of our gain functions as the underlying distributional changes might be smooth rather than piecewise constant (and the tree algorithm enforces a piecewise constant approximation). Hence, here and in general, it will be interesting to see how well optimistic search would perform in such ``misspecified'' scenarios.

At the core of our \emph{optimistic search} strategy for change point detection is the optimization of a piecewise quasi-convex function. As such, there might be connections to a variety of different applications such as finding the maximum of a response surface for e.g.~most effective drug dose \citep{Facer_Mueller}, finding density ridges \citep{Chen_Genovese_Wasserman}, or speeding up optimization in problems where the evaluation of the gradient is too costly and thus one resorts to (noisy) function evaluations \citep{Wang_Balakrishnan_Singh}. Further related examples for the latter could be also hyperparameter selection in some expensive non-parametric procedures.

\subsection{Computational considerations for multiple change points}
In the presence of a single change point, advanced and combined optimistic search only require $O(\log T)$ evaluations of the gain function and yet they are optimal statistically. Can this be generalized for multiple change points? Assuming there are $\ncp$ change points, is there a procedure that localizes change points optimally with $O(\ncp \log T)$ (or maybe $O(\ncp(\log T)^r)$ for some $r\in\N$) evaluations and yet has provable statistical guarantees? OBS, in principle, is a candidate for such a procedure provided that the special cases with signal cancellation effects could be eliminated (see Appendix~\ref{Appendix_OBS}). Note that OSeedBS is optimal statistically (Theorem~\ref{th:osbs}), however, its computational cost is not adapted to the number of change point $\ncp$ but to the minimal spacing between them (which may not be known). The reason is that all seeded intervals above a certain length are considered. Here, an alternative would be to adaptively choose only a subset of all possible seeded intervals instead of considering all of them, thus trying to reduce the computational cost even further. Overall, there seem to be multiple paths to pursue and the precise characterization whether in multiple change point scenarios there is a trade-off between computational cost (number of evaluations) and statistical performance (in terms of localization error or permitted signal classes) remains an open question. Another challenge is, whether going below a logarithmic number of evaluations is possible. Last, but not least, we would like to emphasize that for multiple change points, even without knowing cumulative sums and the minimal segment length between change points, OSeedBS with NOT selection has $O(T)$ computational time in the worst case. We are not aware of other methods that have such a favorable worst case complexity and are yet statistically optimal. We find it surprising that for $T$ observations $O(T)$ worst case complexity is possible.

An open question remains what good ways are for storing data to allow fast change point detection. In the univariate Gaussian change in mean setup storing cumulative sums enables sublinear computational times for change point detection. In other situations, other proper ``sufficient statistics'' might be suitable.  For instance, for detection of changes in quantiles (cf.~\citealp{vanegas2020multiscale}), one might consider to store the data in a tree structure, e.g.~AVL tree \citep{AVL62}.

\subsection{Comments on the theory}\label{ss:theory}
In case of multiple change points, the weakest condition available in literature which allows to detect change points is of the form
\begin{equation}\label{eq:optc}
\delta\sqrt{\lambda T} \ge \sqrt{2 r\log \frac{1}{\lambda}} + b_T
\end{equation}
for any sequence $b_T \to \infty$ arbitrarily slow. In case of only a single bump, it is known that $r = 1$ in \eqref{eq:optc} is optimal (cf.~\citealp{Chan_Walther2013}), while in case of general multiple change points the optimal $r$ remains still unknown, and the best known result is $r = 4 + \epsilon$ with arbitrarily small $\epsilon > 0$ by \citet{Chan_Chen2017} and $r=1$ under certain additional sparsity condition by \citet{JCL10}, which improves the result that $r= 8$ for the bounded number of change points and $r=72$ for the unbounded number of change points, reported earlier in \citet{Frick_etal}. Without further model assumption, it often makes no sense to pursue changes faster than the sampling rate, so we usually have $\lambda \gtrsim 1/T$\,. In the particular case of $\lambda \asymp T^{-\omega}$ with $\omega \in (0,\,1]$, condition~\eqref{eq:optc} is equivalent to
$$
\delta\sqrt{\lambda T} \ge \sqrt{2r\omega \log T}\,.
$$
This is exactly what we require for optimistic search based methods (e.g.~OSeedBS), up to a multiplying constant. Thus, it remains open whether the difference in constants is the price we have to pay, or whether it is simply due to our proof techniques.
Meanwhile, the best known result for multiple change points in the literature concerning the localization error is  $\abs{\hat \cpt_i - \cpt_i}\lesssim \log T / (\delta_i^2T)$, and it is not clear what the best multiplying constant is. Such a localization error rate coincides exactly with our result, see Theorem~\ref{th:osbs}, and again we are not sure whether we have to pay a price in terms of constant or not.

In the particular case of a single change point, the weakest detection condition was recently shown to be $\delta^2\lambda T\gtrsim \log \log T$ by exploring the monotonicity structure, see \citet{Liu_Gao_Samworth_dyadic}.
Moreover, \citet{GHZ20} derived minimax rates in $\ell^2$ risk for isotonic piecewise constant signals. Based on $\ell_2$ risk, we can obtain a better localization error rate for the single change point scenario, which is $\abs{\hat \cpt_1 - \cpt_1}\lesssim \log\log T / (\delta^2 T)$, via the same argument as in \citet{lin2016approximate}. Thus, in terms of both detection condition and localization errors, our result on the advanced and combined optimistic searches has an additional factor of order $\log T/\log\log T$ compared to the optimal ones, see Theorem~\ref{th:aoss}. However, removing such a factor seems to be possible, because the $\log T$ factor in the current result is due to the control of noise part, and there are in fact at most $O(\log T)$ evaluations (which are de facto signal dependent but not noise dependent). Thus by the union bound and the sub-Gaussian tail, it might be possible to control the noise part by $\log \log T$ instead, which offers an interesting avenue for future research in this direction.

\section{Another advanced optimistic search variant}
\label{Appendix_aos}

\begin{algorithm}[H]
\caption{advanced optimistic search, version~2}\label{Alg:aOS2}
\begin{algorithmic}[1]
\Require $0<\nu \leq 0.5; ~\delta<(R-L)/4; ~L,R, \delta \in \mathbb{N}^+$
\Function{\textit{aOS}}{$\nu, \delta, L, R$}
\State $i \leftarrow \lfloor\log_2((R-L)/2)\rfloor$
\State $S \leftarrow \{ L+2, L+4,\ldots,L + 2^i\} \bigcup \{R - 2^i, \ldots, R-4, R-2\}$ \Comment{Preliminary set of dyadic locations that we fine tune in lines $4$ to $8$}
\State Remove points $s\in S$ s.t. $\abs{s-L}<\delta$ or $\abs{s-R}<\delta$ \Comment{Update $S$ by removing points within distance $\delta$ from the boundaries}
\If {$(L + \lfloor(R-L)/2\rfloor)  - (L + 2^i) > 2^{i-1}$}
\State Add $L + \lfloor(R-L)/2\rfloor$ to $S$ \Comment{Add the middle point if the gap is too large}
\EndIf
\If {$(R - 2^i)  - (L + 2^i) < 2^{i-1}$}
\State Replace $L + 2^i$ and $R - 2^i$ in $S$ by $L + \lfloor(R-L)/2\rfloor$  \Comment{Modify the middle points if the gap is too small}
\EndIf
\State $t \leftarrow \argmax\limits_{s \in S} G_{(L,R]}(s)$ \Comment{Localize the best split point on the ``dyadic'' grid $S$}
\State Denote $l$ the largest $s\in S$ s.t. $s<t$ or if there is no such $s$, then set $l=L+(t-L)/2$
\State Denote $r$ the smallest $s\in S$ s.t. $t<s$ or if there is no such $s$, then set $r=R+(t-R)/2$
\State $\hat s_{(L,R]} \leftarrow$ \Call{\textit{nOS}}{$l, t, r \mid \nu, L, R$} \Comment{Run naive optimistic search on the neigbourhood $(l,r]$ of $t$}
\State \textbf{return} $\hat s_{(L,R]}$
\EndFunction
\end{algorithmic}
\end{algorithm}

The essential part of Algorithm~\ref{Alg:aOS2} is line $3$ and lines $9$ to $13$, while lines $4$ to $8$ are of secondary importance. Changing the default $\delta=1$ (in line~$4$) to some higher value is mostly useful/necessary e.g.~in high-dimensional models requiring a minimal number of observations and lines $5$ to $8$ are a slight improvement ensuring that the dyadic points starting from the left boundary $L$ and the dyadic points from the right boundary $R$ can be merged in the middle of the segment around $L+(R-L)/2$ appropriately.

\section{Motivation and special cases for Optimistic Binary Segmentation}
\label{Appendix_OBS}

\subsection{Intuition for OBS}
Figure \ref{Fig:bs_vis} shows an illustrative example for BS (left side) and OBS (right side) in the population case for the univariate Gaussian model \ref{Gaussian_setup} with piecewise constant mean in segments $S_1, S_2, S_3$ and $S_4$. The plotted function (solid line) describes the gain in terms of squared errors when splitting the underlying data at the point given by the $x$-axis. The colored dashes show the boundaries for the underlying segments in iterations II-V. Note that each time the solid line has a local maximum, the coefficients change and a true change point occurs. Moreover, each dot along the gain function represents a model fit in the respective search step of the algorithm within a segment. Whenever the algorithms perform a split, as marked by the magenta crosses, a split point is added in the next iteration, indicated by the vertical colored dashes. If a segment does not contain any more change points, there will be one more search step before concluding that no more splitting is possible.

\begin{figure
\centering
\includegraphics[width = \textwidth]{figures/figApp1.pdf} \\
\caption[Illustrative example for BS compared to OBS]{Illustrative example (in the population case) for BS (left) and OBS (right). The plotted function (solid line) describes the gain depending on the split location and the underlying data. The dots represent the model fits. The dashed vertical lines mark the true underlying change points, whereas the colored dashes mark the already found change points in each iteration. Finally, the magenta crosses mark the selected point for splitting in each iteration and thus the boundary for the underlying segments in the consecutive iteration.}
\label{Fig:bs_vis}
\end{figure}

For BS each split occurs from top to bottom at the respective global maximum in each segment, leading to the final segmentation at the bottom. Optimistic search, however, does not guarantee to find the global maximum in each segment. A key insight regarding OBS, as depicted on the right hand side is, that the sequence of the splits does not matter for obtaining the same final segmentation result at the bottom, as long as each split occurs at one of the true change points. In each search step one is satisfied with finding a local maximum instead of the global maximum within the considered segment, as each local maximum is also one of the true change points. For OBS, the order of splits from top to bottom on the right hand side is different than for BS. In the first split, OBS chooses the local maximum at the boundary of $S_1$ and $S_2$ instead of the global maximum at the boundary of $S_2$ and $S_3$. Choosing a different split point in the first iteration also leads to different subsequent splits. However, the final segmentation results on the bottom are the same for both BS and OBS. Also, OBS needs significantly fewer model fits, as indicated by the black dots along the gain function, as it suffices to look at a very small subset of points in order to find a local maximum in each search~step.

\subsection{Theory for OBS in the population case}
\label{Appendix_OBS_theory}
If there is more than one change point, multiple local maxima may exist even in the population case. In many change point detection problems local maxima of the gain function also correspond to change points. Hence, as long as the search recovers one of the local maxima, we do not include any false positive and after we split our data at the found split point, we can later on still detect the change point that was the cause for the global maximum in the previous iteration. This is not entirely obvious, as the underlying gain functions typically change depending on the boundaries $L$ and $R$. However, the basic structure remains that local maxima still correspond to change points. In the following Lemma~\ref{Lem:OS} we present sufficient conditions for our previous claim,~i.e.~when the naive optimistic search returns at least a local maximum.

\begin{definition}
A function $f:X \mapsto \mathbb{R}$ on a vector space $X$ is quasiconvex if for any $x,y \in X$ and $\lambda \in (0,1)$, one has
\begin{equation*}
f((1-\lambda)x +\lambda y) \leq max\{f(x),f(y)\}
\end{equation*}
and strictly quasiconvex if
\begin{equation*}
f((1-\lambda)x +\lambda y) < max\{f(x),f(y)\}.
\end{equation*}
\end{definition}

\begin{lemma} \label{Lem:OS}
If the (continuously embedded version of the) gain function $G_{(L,R]}(\cdot)$ is piecewise (between the change points) strictly quasiconvex with
$G_{(L,R]}(s) \geq G_{(L,R]}(L) = G_{(L,R]}(R)$ $\forall s\in (L,R)$,
then the naive optimistic search (Algorithm \ref{Alg:OS}) with step size $\nu=1/2$ finds a strict local maximum of $G_{(L,R]}$ after at most $O(\log(R-L))$, and hence after at most $O(\log T)$, steps.
\end{lemma}

One might think that arbitrary search techniques that find the global maximum for unimodal functions in $O(\log (R-L))$ steps (e.g.~ternary search) would work just as well and find local maxima, but this is not the case. Essential in our proposed search is that we keep the information about the previous step in the sense that one of the previous probe points will be one of the new boundary points while the other probe point is going to be one of the new probe points with gain that is thus at least as high as for the new boundary. As we assume that $G_{(L,R]}(s) \geq G_{(L,R]}(L) = G_{(L,R]}(R)$ $\forall s\in (L,R)$, even for the first probe point it is true that the gain is at least as big as on the boundaries and for the above described reason, each step keeps this ``triangular structure'' where the probe point in the middle is guaranteed to be at least as high (and typically even strictly higher) than the boundaries. This triangular structure is also kept by the pre-screening through dyadic locations in advanced optimistic search and hence the result and proof would be analogous to the one presented for the naive version. Gain functions usually fulfill $G_{(L,R]}(L) = 0 = G_{(L,R]}(R)$ and a split somewhere within $(L,R]$ typically leads to a positive gain as one is able to fit two models instead of a single one.

\begin{proof}[Proof of Lemma \ref{Lem:OS}]
We need to show that there exists no iteration after which no more strict local maxima are left in the search domain (Claim~1) and that the number of iterations is of order $\mathcal{O}(\log (R-L))$ (Claim~2). Assume Algorithm~\ref{Alg:OS} is initialized on the interval $(L,R]$ containing at least one strict local maximum and w.l.o.g. assume that the search has been narrowed down to the interval $(l,r]$ that is now split into $(l,s], (s,w]$ and $(w,r]$ by the two current probe points $s$ and $w$. We can distinguish two cases for Claim~1:
\begin{enumerate}[i.]
    \item  There are strict local maxima in at least two of the three segments.
    \item Only one of the segments contains strict local maxima.
\end{enumerate}
Since we only discard one segment in each iteration, in case~i.~the next search interval will still contain at least one local maximum. Regarding case~ii., if the local maxima are in $(s,w]$ the next iteration will also contain the local maxima since $(s,w]$ is contained in the new iteration no matter whether $(l,s]$ or $(w,r]$ is discarded. Hence, the only way that the next iteration would not contain a local maximum is, if in case (ii) the maxima are all located in $(l,s]$ but $G_{(L,R]}(s) \leq G_{(L,R]}(w)$ such that $(l,s]$ is discarded. (The symmetric counterpart when all local maxima are located in $(w,r]$ is analogous).

Assume now that there is no strict local maximum located in $(s,r]$ and that $G_{(L,R]}(s) \leq G_{(L,R]}(w)$ such that we make a mistake at this iteration. From the two probe points one of them was already evaluated in the previous iteration. If we assume that $s$ was evaluated previously, then $G_{(L,R]}(s) \geq G_{(L,R]}(r)$ as the probe point's gain is bigger than that of the boundaries (due to the ``triangular structure'' discussed earlier, below Lemma \ref{Lem:OS}). Thus, $G_{(L,R]}(r) \leq G_{(L,R]}(s) \leq G_{(L,R]}(w)$. In this case, however, $G$ cannot be strictly quasiconvex on $(s,r]$ as the gain at $w\in(s,r]$ is not strictly smaller than both of the gains on the boundaries $s$ and $r$ (and hence also their maximum). Hence, there needs to be at least one more strict local maximum in $(s,r]$, leading to a contradiction. If we assume that $w$ was evaluated previously, then $G_{(L,R]}(w) \geq G_{(L,R]}(r)$ due to the ``triangular structure'', as well as $G_{(L,R]}(s) \leq G_{(L,R]}(w)$ by assumption, which overall again leads to a contradiction with strict quasiconvexity. Claim 2 has been proven in Lemma~\ref{Lem:nOS_speed}.
\end{proof}

In the population case, if each local maximum corresponds to a change point (which is the case in many change point detection problems beyond the univariate Gaussian change in mean problem, see Lemma~\ref{Lem:GGM_convexity} for an example) and if the gain functions are \emph{strictly} piecewise quasiconvex, one would recover a true change point with optimistic search, and thus, with iterative application of it as in OBS (Algorithm \ref{Alg:OBS}), eventually all of the change points would be found (if $m$ for the minimal segment length and $\gamma$ for the minimal required segment length are chosen appropriately). However, as discussed next, there are unfortunately some special cases where these conditions are not fulfilled and OBS might fail.

\subsection{Special cases where OBS might fail}
\label{App:OBS_failure}
Strict quasiconvexity is essential for Lemma~\ref{Lem:OS}. Otherwise, in some regions, the gain function could be flat and the search could make a mistake when deciding which side of the current search interval to discard. Unfortunately, this can occur in some multiple change point detection problems. Luckily, cases where this would happen are fairly special (see e.g.~the middle and right panels of Figure~\ref{f_star_plot}).
Figure~\ref{f_star_plot} shows population gain functions in terms of squared errors (top) for the signals on the bottom. In the left panel mean values in the segments were drawn independently $\mathcal{N}(0,1)$. The middle and right panels show special cases where flat parts in the gain curve occur. In the Gaussian change in mean problem such flat parts in the population gain (denoted $G^{*}$ here) within a segment $(l,r]$ occur only if $\mu_{(l,\tau_{j}]}^{*} = \mu_{j+1} = \mu_{(\tau_{j+1}, r]}^{*}$ for some $j \in \{1,\ldots,\ncp\}$, i.e.~if the means before the segment and after the segment are equal to the mean within the segment.
To see this, one can look at the second derivative of the population gains with respect to the split point $s$ (similar to Lemma~\ref{Lem:GGM_convexity} on covariance matrices).
Case 2) on the middle of Figure~\ref{f_star_plot} shows such a scenario with the same location of the change points as in case~1), but with up/down jumps that exactly cancel each other resulting in flat parts for $G^{*}$. Case 3) on the right panel shows Example~\ref{Example_failure_OBS}. Overall, we can conclude that the population gain functions are piecewise convex and in most scenarios even piecewise strictly convex, unless some very special signal cancellation is happening \citep[see also][Lemma 2.2]{ven92}.

\begin{figure}
\includegraphics[width=0.95\textwidth]{figures/figApp2.pdf}
\caption{The (population) gains in terms of squared errors (top) together with the underlying mean values in the segments (bottom). Case 3) on the right panel depicts Example~\ref{Example_failure_OBS}.}
\label{f_star_plot}
\end{figure}

Flat parts of the gain function for the Gaussian change in mean example can only happen if the mean value for a segment is exactly the same as the mean of all the data before this segment and the mean of all the data after this segment. In this cases there is no gain in splitting. In a very adverse scenario (see Example~\ref{Example_failure_OBS} plotted on the right panel of Figure~\ref{f_star_plot}) it could happen both for the naive and advanced optimistic search that it gets stuck in a flat region which means that one stops searching for further change points even though there would be some left. However, in general, even if signal cancellation occurs for some parts of a signal, optimistic search may not fail. If at least one of the probe points that are evaluated during the search has a positive gain, one already starts to climb up to one of the local maxima. For example, in the middle panel of Figure~\ref{f_star_plot} optimistic search with $\nu=1/2$ would even find the global maximum in the first iteration even though some flat regions occur. Also note that splitting at the found maximum changes the shape of the gain function in the resulting two segments such that the signal canceling phenomenon may not even happen in the searches later on.

However, as long as special cases exist even in the population case where OBS could fail, the same failures will occur also in the noisy case. Some ideas that might lead to (theoretical) guarantees could be the following:
\begin{enumerate}[i.]
    \item Restricting to signals where signal cancellation does not happen, for example if all jumps are upwards e.g.\ a stair-like signal.
    \item The essential problem of optimistic search is that in some cases it can get stuck in a flat part which results in missing some change points. One could try to show that in this case at least a subset of change points will be found (and no false positives are picked on the way).
    \item One could choose a smaller step size or adapt the step size to the current segment length. As more split points are evaluated, it is more likely that at least one is evaluated with a positive gain and thus one starts to climb up to a local maximum. The hope is that the step sizes can be translated into guarantees of the length of segments that can be found. However, clearly, this way one loses some of the computational efficiency.
    \item When the gains of both probe points are equal, the choice whether to throw away $(l,s]$ or $(w,r]$ is taken randomly. However, if one looks at the sum of squared errors of the segment $(l,s]$ vs. $(w,r]$ (in the Gaussian mean change case), then the segment containing a change point will have larger value as it contains additionally to the usual variance at least one pair of up and down jumps. Hence, keeping the segment with the larger sum of squares is advisable. This idea could also help to speed up OSeedBS in the sense that one could prioritize which seeded intervals are the most promising based on their variance. However, it is somewhat unclear how far this variance check idea generalizes beyond the Gaussian change in mean problem.
\end{enumerate}

\section{Additional material on the univariate Gaussian simulations}
\label{Appendix_Simulations}

Figure~\ref{Fig:single_cpt} shows found change points using various search methods in each 1,000 simulations for a balanced ($n=200$) and an unbalanced ($n=\text{5,000}$) scenario. The failure of the naive optimistic search in most cases for the unbalanced scenario is again clearly visible, while for the advanced, and in particular the combined optimistic search, the found change points very often lie exactly on the diagonal when compared to the full grid search and hence exactly the candidate proposed by the full grid search were found. Together with the simulation results in Table~\ref{tab:simulation_results}, this confirm our theoretical results that the naive optimistic search is not consistent for very unbalanced signals, while the advanced and combined versions are.

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{figures/figApp3.png} \\
\caption{Pairwise plots of found change points using different optimistic search methods (horizontal-axis) versus the ones returned by the full grid search (vertical-axis) for a noise level $\sigma = 1.5$ and $n=200$ (top) as well as $n=\text{5,000}$ (bottom) in $1,000$ simulations from Example~\ref{setup:single_mean_shift}. The vertical and horizontal lines indicate the location of the true change point at observation~$100$.}
\label{Fig:single_cpt}
\end{figure}

\section{Change point detection for high-dimensional Gaussian graphical models}
\label{Appendix_highdim}

In the following, we briefly describe an estimator introduced by \cite{Londschien} for change point detection in high-dimensional Gaussian graphical models, as this is the basis of all change point detection algorithms (BS, SeedBS, WBS and their optimistic variants) that we investigate in Section~\ref{highdim_simulations}. For a segment $(u,v]$ with $0 \leq u < v \leq T$ let $S_{(u,v]}$ denote the empirical covariance matrix within that segment. Let $0 < \delta < 1/2$ be some required minimal relative segment length and $\lambda>0$ be a regularization parameter. According to a proposal of \cite{Londschien}, for a segment $(u,v]$ with $v-u>2\delta T$ we define the split point candidate as
\begin{multline}
    \label{glasso_change_est}
    \hat\eta_{(u,v]}
= \argmax\limits_{s \in \{ u+\delta T,\ldots, v-\delta T\} }
L_T(\hat\Omega_{(u, v]}^{\mathrm{glasso}}; S_{(u, v]}) \\
- \left(L_T(\hat\Omega_{(u, s]}^{\mathrm{glasso}}; S_{(u,s]}) +
L_T(\hat\Omega_{(s, v]}^{\mathrm{glasso}}; S_{(s, v]}) \right),
\end{multline}
where
\begin{equation*}
L_T(\Omega; S_{(u,v]}) =
\frac{v-u}{T}\left(\Tr(\Omega S_{(u,v]}) - \log(|\Omega|) \right)
\end{equation*}
is a multivariate Gaussian log-likelihood based loss in the considered segment $(u,v]$ (scaled according to its length), and $\hat\Omega_{(u,v]}^{\mathrm{glasso}}$ is the graphical lasso precision matrix estimator \citep{glasso} with a scaled regularization parameter $\sqrt{T/(v-u)}\lambda$, i.e.,
\begin{equation*}
\label{eq:glasso_estimator}
\hat\Omega_{(u,v]}^{\mathrm{glasso}} = \argmin\limits_{\mathbb{R}^{p \times p} \ni \Omega \succ 0} \Tr(\Omega S_{(u,v]}) - \log(|\Omega|) + \sqrt{(T / (v- u))}\lambda \|\Omega\|_1.
\end{equation*}

Each split point $s$ in \eqref{glasso_change_est} requires fitting two graphical lasso estimators. While various algorithms for computing exact or approximate solutions for the graphical lasso estimator exist, scalings of $O(p^3)$ (or worse) are common (assuming that the input covariance matrices have been pre-computed and that no special structures such as the block diagonal screening proposed by \cite{Witten2011_connected_components} as well as \cite{Mazumder2012_connected_components} can be exploited). Hence, the graphical lasso is rather costly especially when repeatedly fitting at each possible split point $s$ on the grid $u+\delta T,\ldots, v-\delta T$. Note that the essential problem is that it is not easy to re-use the estimator $\hat\Omega_{(u,s]}^{\mathrm{glasso}}$ for the segment $(u, s]$ to obtain $\hat\Omega_{(u,s+1]}^{\mathrm{glasso}}$ for the segment $(u,s+1]$. One could use $\hat\Omega_{(u,v]}^{\mathrm{glasso}}$ as a warm start, but not all algorithms that have been developed to compute graphical lasso fits are guaranteed to converge with warm starts (see \citealp{Mazumder2012_dpglasso}) and even the ones that do converge would not save orders of magnitude in terms of runtime. Note, that this lack of efficient updates is common for more complex (e.g.~high-dimensional) scenarios and it is in sharp contrast with e.g.~change point detection for univariate Gaussian variables. There, one needs to calculate means, but the mean for the segment $(u,s+1]$ can be updated in $O(1)$ cost if the mean for the segment $(u,s]$ is already available, and hence the computational cost is typically proportional to the total length of considered segments. In contrast, in the estimator in~\eqref{glasso_change_est}, the number of graphical lasso fits, as given by the number of considered split points, is the main driver of computational cost. Our optimistic search techniques rely on evaluating far fewer split points $s$ than the full grid search and thus provide an option for massive computational speedups. Of course, the price to pay is having no guarantee to obtain exactly the optimal split point, but the ``optimistic'' approximation to $\hat\eta_{(u,v]}$, i.e., the one obtained via the optimistic search, is still fairly good, see the simulations in Section~\ref{highdim_simulations}.

In the simulations, we used the \textbf{glasso} \textsf{R} package, available on CRAN, for the graphical lasso fits. For all six methods, we set $\delta=0.01$, i.e.\ skipping $20$ observations on the boundaries of each considered interval and overall, no change points were searched in intervals containing less than $60$ observations. We set $\lambda = 0.007$. Regularization in these examples is not essential in the sense that we do not have truly high-dimensional scenarios, but for split points close to the boundaries of the search interval and in short intervals, where the number of observations is close to $p$, regularization can be still helpful. We could have increased $p$ in order to cover truly high-dimensional setups in our simulations, but given the scaling $O(p^3)$ of the graphical lasso, this very quickly goes beyond reasonable computational times for the full grid search based approaches that we want to include as references in terms of achievable estimation error.

\subsection{Proof of Lemma~\ref{Lem:GGM_convexity}}
\label{App:GGM_convexity_proof}
First note that in the following we only consider the case $l<\cpt_i < s < \cpt_{i+1}< r$, but similar arguments can be used also in the presence of a single change point in $(l,r]$ or when considering a split point $s$ in the segment from $l$ to the first change point within $(l,r]$.
Recall that $\Sigma_{(l,r]}$ denotes the convex combination of the covariance matrices within the segment $(l,r]\subseteq(0,1]$ with the weights given by the relative segment lengths within $(l,r]$. In particular, for $l<\cpt_i < s < \cpt_{i+1}< r$,
\begin{equation*}
\Sigma_{(l,s]}
=\frac{1}{s-l}
\left(
(\cpt_i - l)\Sigma_{(l,\cpt_i]}
+
(s-\cpt_i)\Sigma_{i+1}
\right)
\end{equation*}
and
\begin{equation*}
\Sigma_{(s,r]}
=\frac{1}{r-s}
\left(
(\cpt_{i+1}-s)\Sigma_{i+1}
+
(r - \cpt_{i+1})\Sigma_{(\cpt_{i+1}, r]}
\right),
\end{equation*}
where $\Sigma_{i+1} = \Sigma_{(\cpt_{i}, \cpt_{i+1}]}$ is the covariance matrix in the $i+1$-st segment $(\cpt_{i}, \cpt_{i+1}]$. We seek to find the first and second derivatives of $G^*_{(l,\,r]}(s)$. First note that,
\begin{equation*}
\frac{\partial}{\partial s} \Sigma_{(l,s]}
=
\frac{\cpt_{i} - l}{(s-l)^2}
\left(\Sigma_{i+1} - \Sigma_{(l,\cpt_i]} \right)
\end{equation*}
and
\begin{equation*}
\frac{\partial^2}{\partial s^2} \Sigma_{(l,s]}
=
-\frac{2(\cpt_{i} - l)}{(s-l)^3}
\left(\Sigma_{i+1} - \Sigma_{(l,\cpt_i]} \right)
=
-\frac{2}{s-l}\left(\frac{\partial}{\partial s} \Sigma_{(l,s]}\right).
\end{equation*}
We need the following expressions (see e.g.~\citealp{matrix_cookbook}) for derivatives of an invertible matrix $A(s)$ depending on $s$,
\begin{equation*}
\frac{\partial}{\partial s} \log(|A(s)|) = \Tr\left( A(s)^{-1} \cdot \frac{\partial}{\partial s} A(s)\right);
\end{equation*}
\begin{equation*}
\frac{\partial}{\partial s} \Tr(A(s)) = \Tr\left( \frac{\partial}{\partial s} A(s)\right);
\end{equation*}
\begin{equation*}
\frac{\partial}{\partial s} A(s)^{-1} = -A(s)^{-1} \cdot \frac{\partial}{\partial s} A(s)  \cdot A(s)^{-1}.
\end{equation*}
We next compute the first and second derivatives of $G^{*}_{(l,r]}(s)$. Recall from equation~\eqref{e:gainHD} that
\begin{equation*}
    G^*_{(l,r]}(s)
    =
    \frac{r-l}{T}\log(|\Sigma_{(l,r]}|)
    -\frac{s-l}{T}\log(|\Sigma_{(l,s]}|)
    -\frac{r-s}{T}\log(|\Sigma_{(s,r]}|)\,.
\end{equation*}
Consider first the middle part of $G^{*}_{(l,r]}(s)$, i.e.\
$$
L^*(s) := -\frac{s-l}{T}\log(\abs{\Sigma_{(l,s]}})\,.
$$
Then for the first derivative
$$
L^{*'}(s)
=
-\frac{1}{T}\log(|\Sigma_{(l,s]}|)
-\frac{s-l}{T}
\Tr\left( \Sigma_{(l,s]}^{-1} \cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right),
$$
and for the second derivative
\begin{align*}
L^{*''}(s) = &
-\frac{1}{T}\Tr\left( \Sigma_{(l,s]}^{-1} \cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right)
-\frac{1}{T}
\Tr\left( \Sigma_{(l,s]}^{-1} \cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right) \\
& \mbox{}\qquad{}\qquad{}\qquad{}
-\frac{s-l}{T}
\frac{\partial}{\partial s}
\left(
\Tr\left( \Sigma_{(l,s]}^{-1} \cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right)
\right)\\
= &
-\frac{2}{T}\Tr\left( \Sigma_{(l,s]}^{-1} \cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right)
-\frac{s-l}{T}
\Tr\left(
-\Sigma_{(l,s]}^{-1} \cdot \frac{\partial}{\partial s} \Sigma_{(l,s]}  \cdot \Sigma_{(l,s]}^{-1}
\cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right)\\
& \mbox{}\qquad{}\qquad{}\qquad{}
-\frac{s-l}{T}
\Tr\left( \Sigma_{(l,s]}^{-1} \cdot
\frac{\partial^2}{\partial s^2} \Sigma_{(l,s]}
\right)\\
= &
\frac{s-l}{T}
\Tr\left(
\Sigma_{(l,s]}^{-1} \cdot \frac{\partial}{\partial s} \Sigma_{(l,s]}  \cdot \Sigma_{(l,s]}^{-1}
\cdot
\frac{\partial}{\partial s} \Sigma_{(l,s]}
\right)
\\
= & \frac{s-l}{T}
\left\|{\Sigma_{(l,s]}^{-1/2} \cdot \frac{\partial}{\partial s} \Sigma_{(l,s]}  \cdot \Sigma_{(l,s]}^{-1/2}}\right\|_{F}^2 \ge 0.
\end{align*}
By symmetry, we can obtain similarly for the right part of $G^{*}_{(l,r]}(s)$,
$$
\frac{\partial^2}{\partial s^2}\left(-\frac{r-s}{T}\log(|\Sigma_{(s,r]}|)\right) = \frac{r-s}{T}
\left\|{\Sigma_{(s,r]}^{-1/2} \cdot \frac{\partial}{\partial s} \Sigma_{(s,r]}  \cdot \Sigma_{(s,r]}^{-1/2}}\right\|_{F}^2 \ge 0.
$$
As the left part of $G^{*}_{(l,r]}(s)$ is constant, we have
$G^{*''}_{(l,r]}(s)\geq 0$ for $\cpt_i < s < \cpt_{i+1}$ and hence, $G^{*}_{(l,r]}$ is convex in between change points $\cpt_i$ and $\cpt_{i+1}$. Moreover, we see that with the exception of the special cases when $\Sigma_{(l,\cpt_i]} = \Sigma_{i+1} = \Sigma_{(\cpt_{i+1},s]}$, $G^{*}_{(l,r]}$ is even strictly convex in the interval $(\cpt_{i}, \cpt_{i+1}]$. Note that for such special cases $G^{*}_{(l,r]}(s)=0$ for arbitrary $s \in (\cpt_i, \cpt_{i+1}]$, i.e.~the population gain function is flat in between change points $\cpt_i$ and $\cpt_{i+1}$. These special cases are thus very similar to those for the univariate Gaussian change in mean case discussed in Appendix~\ref{App:OBS_failure}  (see e.g.~Figure~\ref{f_star_plot}). Note that special cases can only occur in the presence of two or more change points within the considered segment. In particular, in case $\cpt_i$ is the single change point contained in $(l,r]$, $G^{*}_{(l,r]}$ is strictly convex in $(l,\cpt_i]$ and strictly convex in $(\cpt_i,r]$.

\section{Comments on computational gains for high-dimensional simulations}
\label{Appendix_computational_gains_highdim}

The achievable speedups using optimistic search in general are dependent on the cost of the model fit in each segment (how they depend on the number of observations $T$ and the dimensionality $p$), whether there are possibilities to update neighboring fits efficiently, but also on the length of the series, the number of change points, which basic algorithm (BS, SeedBS, WBS or yet another one) is used with which specific tuning parameters, etc. Nonetheless, we would like to further comment on some of the observed computational gains in the high-dimensional simulations presented in Section~\ref{highdim_simulations}.

The biggest computational gains for optimistic search occur when the underlying search intervals are long. Random intervals have expected length $O(T)$ and thus many of them are comparably long. For these long intervals we gain a lot by optimistic search. However, the lengths of the intervals in lower layers of seeded intervals are quite short (decaying exponentially) and what becomes dominant in that case is the number of very short intervals. For example, while there is only a single interval containing 2,000 observations (first layer), there were more than sixty intervals on the lowest layer we considered with the minimally required segment length of $m=60$ observations. This explains why the speedup for OSeedBS versus SeedBS is a factor 10--14, while for BS and WBS we could achieve factor $30$ or more. Skipping the last few layers of seeded intervals would have saved considerable computational time for OSeedBS, which is important to keep in mind when interpreting the results from Figure~\ref{Fig:glasso_multiple_changes}. From a practical perspective, when utilizing OSeedBS, one should thus limit the number of covered layers in order to consider fewer of the very short intervals that are a driver of computational cost. However, the minimal segment length in seeded intervals cannot be too large either, as in that case one is risking not covering each single change points sufficiently (similar to what happened in the shown examples for WBS and OWBS with a small number of random intervals $M$). The choice for the minimal segment length for seeded intervals might come fairly naturally in some applications, where segments below a certain size are uninteresting or when considering high-dimensional problems requiring a minimal number of observations for fitting reasonable models.

A pragmatic approach could be to combine the best of both worlds from OBS and OSeedBS. For example, find a first set of change points with fewer number of seeded intervals and then, to protect against the possibility that there could be even further change points that were not discovered due to having chosen a too large minimal segment length, in between the first found change points from the seeded intervals, one could perform a further OBS-like search that adapts better to the number of change points within these shorter search intervals. This way adaptively one could invest more computational effort if there is evidence for further change points beyond the ones found by the rough first set of seeded intervals, but without the need to go over each and every very short interval as would be the case with further layers of seeded intervals containing very short intervals. Thus, one could keep computational advantages from OBS and at the same time exploit the better expected estimation performance of OSeedBS.

\input{sections/proofs}

\section{Proofs of statistical guarantees}

Here we provide proofs of statistical guarantees for optimistic searches in terms of consistency and localization rates. For ease of reading, we rewrite Model~\ref{Gaussian_setup} as
$$
X_t = f_t + \xi_t \qquad \text{for}\quad t = 1, \ldots, T\,,
$$
where $f_t = \sum_{i = 0}^\ncp \mu_i \ind_{(\cpt_i, \cpt_{i+1}]}(t/T)$ and $\xi_t \stackrel{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0,1)$, i.e.~standard Gaussian distributed. Let $X = (X_1,\ldots, X_T)^{\top} $, $f = (f_1, \ldots, f_T)^{\top}$ and $\xi = (\xi_1, \ldots, \xi_T)^{\top}$ in $\R^T$. Then it holds that $\xi\sim \mathcal{N}(0, I_T)$ and $X \sim \mathcal{N}(f, I_T)$ with $I_T$ the identity matrix in $\R^{T \times T}$. By $\csm^{g}_{(l,r]}(s)$, we denote the CUSUM statistics in \eqref{CUSUM} being applied to $g$ instead of $X$, for $g = f,\, \xi$ or $X$. In particular, we have $\csm^{X}_{(l,r]}(s) = \csm_{(l,r]}(s)$ and $\abs{\csm^{f}_{(l,r]}(s)} = \abs{\csm_{(l,r]}^*(s)}$\,. In all the proofs, we try to give constants as explicitly as possible, but those constants might not be the best ones.

\subsection{CUSUM of noise}

We start with some properties of the CUSUM statistics $\csm^{\xi}_{(l,r]}(s)$. The first result is about the maximal size of $\csm^{\xi}_{(l,r]}(s)$, which is a known fact and can be easily proven by the union bound (i.e.~Boole's inequality). For completeness, we provide a proof here.

\begin{lemma}\label{le:max}
Let $\xi \in \R^T$ and $\xi \sim \mathcal{N}(0, I_T)$. Then as $T \to \infty$
$$
 \Prob{\max_{\stackrel{0 \le l < s < r \le T}{l,s,r\in\N}} \abs{\csm^{\xi}_{(l,r]}(s)} < \sqrt{6 \log T}} \to 1\,.
$$
\end{lemma}

\begin{proof}
Note that $\csm^{\xi}_{(l,r]}(s) \sim \mathcal{N}(0,1)$ for all possible $l$, $s$ and $r$. Then
\begin{align*}
    \Prob{\max_{\stackrel{0 \le l < s < r \le T}{l,s,r\in\N}} \abs{\csm^{\xi}_{(l,r]}(s)} \ge \sqrt{6 \log T}} & \le \frac{T^3}{6} 2\Bigl(1-\Phi\bigl(\sqrt{6\log T}\bigr)\Bigr) \\
    & \le \frac{T^3}{3\sqrt{2 \pi \cdot 6 \log T}}\exp\left(-\frac{6\log T}{2}\right)\\
    & = \frac{1}{6\sqrt{3 \pi \log T}} \to 0\qquad\text{as}\quad T \to \infty\,,
\end{align*}
where the second inequality is due to the Mills ratio and $\Phi$ is the cumulative distribution function of a standard Gaussian random variable.
\end{proof}

The next result examines correlations between  $\csm^{\xi}_{(l,r]}(s)$ and $\csm^{\xi}_{(l,r]}(t)$, which is especially useful when $s$ and $t$ are close to each other.

\begin{lemma}\label{le:cor}
Let $\xi \in \R^T$ and $\xi \sim \mathcal{N}(0, I_T)$ and define
$$
\mathcal{I}_T = \left\{(l,s,t,r)\,:\, 0 \le l < s <t < r \le T,\,l,s,t,r \in \N\right\}\,.
$$
Then as $T \to \infty$
$$
\Prob{\max_{(l,s,t,r) \in \mathcal{I}_T}  \abs{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)}\sqrt{\frac{\min\{s-l,r-t\}}{t-s}} \le 6\sqrt{\log T}} \to 1\,.
$$
\end{lemma}

\begin{proof}
We split $\mathcal{I}_T$ into two subsets
$$
\mathcal{I}_{T}^0 = \left\{(l,s,t,r)\in \mathcal{I_T}\,: \,\frac{t-s}{\min\{s-l,r-t\}} \ge 1\right\}\,,
$$
and
$$
\mathcal{I}_{T}^1 = \left\{(l,s,t,r)\in \mathcal{I_T}\,: \,\frac{t-s}{\min\{s-l,r-t\}} \le 1\right\}\,.
$$
For $\mathcal{I}_T^0$, by Lemma~\ref{le:max} we obtain
\begin{align*}
& \Prob{\max_{(l,s,t,r) \in \mathcal{I}_T^0}  \abs{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)}\sqrt{\frac{\min\{s-l,r-t\}}{t-s}} \le 6\sqrt{\log T}} \\
\ge \,& \Prob{\max_{(l,s,t,r) \in \mathcal{I}_T^0}  \abs{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)} \le 6\sqrt{\log T}} \\
\ge\, & \Prob{\max_{\stackrel{0 \le l < s < r \le T}{l,s,r\in\N}} \abs{\csm^{\xi}_{(l,r]}(s)} \le 3\sqrt{ \log T}} \to 1 \quad\text{as}\quad T\to \infty\,.
\end{align*}

We next consider $\mathcal{I}_T^1$\,. Note that $ \csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)$ is Gaussian distributed with mean zero. This implies that we only need to control its variance. Let $n=r-l$. Then
\begin{multline}\label{eq:vcs}
\var{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)} \\
= \sum_{i = l+1}^s\left(\sqrt{\frac{r-s}{n(s-l)}}-\sqrt{\frac{r-t}{n(t-l)}}\right)^2 +
\sum_{i = s+1}^t \left(\sqrt{\frac{s-l}{n(r-s)}}+\sqrt{\frac{r-t}{n(t-l)}}\right)^2 \\+
\sum_{i = t+1}^r \left(\sqrt{\frac{s-l}{n(r-s)}}-\sqrt{\frac{t-l}{n(r-t)}}\right)^2\,.
\end{multline}
For the first term on the right hand side of~\eqref{eq:vcs}, we apply the mean value theorem to the function $x\mapsto\sqrt{(r-x)/(x-l)}$ with $x \in [s,\,t]$, and obtain
\begin{align*}
    &\sum_{i = l+1}^s\left(\sqrt{\frac{r-s}{n(s-l)}}-\sqrt{\frac{r-t}{n(t-l)}}\right)^2 &&\\
    = \,\,& \frac{n(s-l)(t-s)^2}{4(r-\theta)(\theta- l)^3} && \text{for some } \theta \in [s,t]\,, \\
   \le\,\,& \frac{n(t-s)^2}{4(r-\theta)(\theta- l)^2} &&\text{since } \theta - l \ge s -l \,,\\
   \le \,\, &\frac{(t-s)^2}{2  \min\{r-\theta, \theta -l\}^2}&& \text{since }\max\{r-\theta, \theta -l\} \ge \frac{n}{2} \,,\\
   \le \,\, & \frac{(t-s)^2}{2 \min\{r-t, s -l\}^2}&&\text{since } r-\theta \ge r -t,\, \theta - l \ge s - l\,,\\
   \le \,\,&\frac{t-s}{2 \min\{r-t, s -l\}}  && \text{since } (l,s,t,r) \in \mathcal{I}_T^1\,.
\end{align*}
For the second term on the right hand size of~\eqref{eq:vcs}, we have \begin{align*}
&\sum_{i = s+1}^t \left(\sqrt{\frac{s-l}{n(r-s)}}+\sqrt{\frac{r-t}{n(t-l)}}\right)^2 && \\
\le \,\, &2(t-s) \left(\frac{s-l}{n(r-s)} + \frac{r-t}{n(t-l)}\right) && \text{since }(a + b)^2\le 2a^2 + 2b^2 \text{ for } a, b \in \R\,,\\
\le \,\, & 2(t-s)\left(\frac{1}{r-s} + \frac{1}{t-l}\right) && \text{since } s-l < n,\, r-t < n\,,\\
\le \,\, & \frac{4(t-s)}{\min\{r-t, s-l\}}&&\text{since }r-s > r-t,\, t-l > s-l\,.
\end{align*}
For the third term on the right hand size of~\eqref{eq:vcs}, in a similar way as for the first term we obtain
$$
\sum_{i = t+1}^r \left(\sqrt{\frac{s-l}{n(r-s)}}-\sqrt{\frac{t-l}{n(r-t)}}\right)^2 \le \frac{t-s}{2 \min\{r-t, s -l\}}\,.
$$
Combining the three terms, we have
$$
\var{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)} \le
\frac{5(t-s)}{\min\{r-t, s-l\}}\,.
$$
In such a case, it holds
$$
\Prob{\abs{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)} \ge x\sqrt{\frac{5(t-s)}{\min\{r-t, s-l\}}}} \le 2\bigl(1-\Phi(x)\bigr)\quad \text{for }x \ge 0\,.
$$
Then, by the union bound (cf.~the proof of Lemma~\ref{le:max}), we can prove as $T\to\infty$\,,
$$
\Prob{\max_{(l,s,t,r) \in \mathcal{I}_T^1}  \abs{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)}\sqrt{\frac{\min\{s-l,r-t\}}{t-s}} \le 6\sqrt{\log T}} \to 1\,.
$$
Since $\mathcal{I}_T = \mathcal{I}_T^0 \cup \mathcal{I}_T^1$, we have shown the assertion of the lemma.
\end{proof}

As a remark, we stress that for Lemmata~\ref{le:max} and~\ref{le:cor} one can choose the constant in front of $\sqrt{\log T}$ sufficiently large to ensure the convergence of probability towards one at any given polynomial rate $O(T^{-k})$ with $k > 0$.
\subsection{Single change point}

The following proofs rely on the observation that the localization error is no larger than the minimal length of search intervals that contain the only true change point. Moreover, assume that the search interval in a step still contains the change point, i.e.~no mistake has been made yet. Then excluding the segment containing the true change point, and thus making a mistake, can only happen when both probe points lie to the left of the change point or when both probe points lie to the right of the change point. In such cases, in order to avoid wrongly excluding the segment containing the true change point, we have to ensure that the difference of $\csm_{(0,T]}^f(\cdot)$ at two investigated probe points is larger than the oscillation ($\asymp \sqrt{\log T}$) caused by the noise.

\begin{proof}[Proof of Theorem~\ref{th:noss}]
Note that part~ii.~is a strengthened version of part~i., as it requires further to identify the correct number of change points, i.e.~$\hat\ncp = \ncp$. Thus, we only need to prove part~ii., and then part~i.~follows immediately. Since the CUSUM statistic in~\eqref{CUSUM} remains unchanged if we use instead $X + C$ for every constant $C\in \R$, we assume w.l.o.g.
$$
f(x) = \delta\ind_{(0, \lambda]}(x) \qquad x \in (0,1]
$$
which in particular implies $\cpt = \lambda$. For notational simplicity, we further assume there is no rounding in defining probe points $s$ and $w$ (cf.~lines 1, 7 and 13 of Algorithm~\ref{Alg:OS}) and $\nu= 1/2$, which implies that the three parts within each search window have relative lengths $1:1:1$, $1:1:2$ or $2:1:1$.

We introduce the event
\begin{multline}\label{eq:good}
    \mathcal{A} = \left\{\max_{\stackrel{0 \le l < s < r \le T}{l,s,r\in\N}} \abs{\csm^{\xi}_{(l,r]}(s)} < \sqrt{6 \log T}\right\}\, \bigcap\\
    \left\{\max_{(l,s,t,r) \in \mathcal{I}_T} \abs{\csm_{(l,r]}^\xi(s)-\csm_{(l,r]}^\xi(t)} \sqrt{\frac{\min\{s-l,r-t\}}{t-s}} \le 6\sqrt{\log T}\right\}
\end{multline}
with $\mathcal{I}_T$ defined in Lemma~\ref{le:cor}. By Lemmata~\ref{le:max} and~\ref{le:cor}, it holds that $\Prob{\mathcal{A}} \to 1$ as $T \to \infty$. Then, it is sufficient to prove $\abs{\hat\cpt -\cpt} \lesssim \log T / (\delta^2 T)$ and $\hat\ncp = \ncp $ on the event $\mathcal{A}$. Thus, in what follows we always assume that the event $\mathcal{A}$ occurs. Let the constant $C_0$ in \eqref{eq:area} satisfy
\begin{equation}\label{eq:nc0}
C_0 \ge 292\sqrt{6}\,.
\end{equation}

Consider now the naive optimistic search on $(0,T]$. We number the search intervals naturally as $I_1, I_2, \ldots$ such that $I_1 \supsetneqq I_2 \supsetneqq \cdots$, i.e.~$I_1 = (0,\,T]$, and either $I_2 = (0,\,2T/3]$ or $I_2 = (T/3,\, T]$,\textellipsis  and define
$$
k_* = \min\left\{k \,\mid\, \min\{s, w\} \le \cpt T  \text{ with $s,w$ the probe points of $I_k$}\right\},
$$
i.e., when some probe point drops in $[0, \cpt T]$ for the first time. We claim
\begin{equation}\label{eq:clm}
  \cpt T \equiv \lambda T \in I_k \qquad \text{for all } k \le k_*\,.
\end{equation}

We prove this claim by induction. Note that $\cpt T \in I_1 = (0,\, T]$\,. If $k < k_*$ and $\cpt T \in I_k$, then the probe points $s,w$ of $I_k$ lie to the right of $\cpt T$, and the left end point of $I_k$ is $0$\,. We apply the mean value theorem to the function $x \mapsto \sqrt{(T-x)/x}$ with $x \in [\min\{s,w\},\,\max\{s,w\}]$ and obtain
\begin{align*}
&\csm_{(0,T]}^f(\min\{s,w\}) - \csm_{(0,T]}^f(\max\{s,w\}) \\
=\, & \lambda \delta \sqrt{T}\left(\sqrt{\frac{T- \min\{s,w\}}{\min\{s,w\}}}- \sqrt{\frac{T - \max\{s,w\}}{\max\{s,w\}}}\right) \\
= \, & \frac{\lambda \delta T^{3/2}}{2(T - \theta)^{1/2}\theta^{3/2}}\abs{s-w}\qquad \text{for some } \theta \in [\min\{s,w\},\, \max\{s,w\}]\,.
\end{align*}
Since $(T-\theta)^{1/2}\theta^{1/2} \le T/2$, and $\abs{s-w}\ge \max\{s,w\}/3$ (which is due to $\nu = 1/2$ and the left end point of $I_k$ being $0$), we further obtain
$$
\csm_{(0,T]}^f(\min\{s,w\}) - \csm_{(0,T]}^f(\max\{s,w\})
\ge  \frac{\lambda\delta \sqrt{T}}{\max\{s,w\}}\abs{s-w} \ge \frac{\lambda\delta\sqrt{T}}{3}\,.
$$
By the assumption~\eqref{eq:area} with $C_0$ in~\eqref{eq:nc0}, we have
\begin{align*}
    &\csm_{(0,T]}^X(\min\{s,w\}) - \csm_{(0,T]}^X(\max\{s,w\}) \\
    =\, & \csm_{(0,T]}^f(\min\{s,w\}) - \csm_{(0,T]}^f(\max\{s,w\}) \\
    & {}\qquad\qquad{} + \csm_{(0,T]}^{\xi}(\min\{s,w\}) - \csm_{(0,T]}^{\xi}(\max\{s,w\}) \\
    \ge \, &
    \frac{\lambda\delta\sqrt{T}}{3} - 2\sqrt{6 \log T} \\
    >\, & 95\sqrt{6\log T}
    >  0\,.
\end{align*}
As the probe point $\min\{s,w\}$ has the higher test statistics, we keep it and discard the outer right segment in Algorithm~\ref{Alg:OS}. Thus, $\max\{s,w\}$ becomes the new right end point of the next search interval, i.e., $I_{k+1} = (0, \max\{s,w\}]$ and $\cpt \in I_{k+1}$. Thus, by induction, we have proved the claim in~\eqref{eq:clm}.

We deal next with $k > k_*$. Then the probe points $s,w$ of $I_k$ satisfy that $s,w \le 3 \lambda T$. Note that the true change point $\cpt T$ can be wrongly excluded from consecutive search intervals only when both probe points $s,w$ lie on the same side of $\cpt T$.

If $s,w \le \cpt T$ (i.e.~left side of $\cpt T$), then again by the mean value theorem
\begin{align*}
    &\csm_{(0,T]}^f(\max\{s,w\}) - \csm_{(0,T]}^f(\min\{s,w\})\\
=\, & (1-\lambda) \delta \sqrt{T}\left(\sqrt{\frac{\max\{s,w\}}{T - \max\{s,w\}}} - \sqrt{\frac{\min\{s,w\}}{T- \min\{s,w\}}}\right) \\
= \, & \frac{(1-\lambda) \delta T^{3/2}}{2(T - \theta)^{3/2}\theta^{1/2}}\abs{s-w}\qquad \text{for some } \theta \in [\min\{s,w\}, \max\{s,w\}] \,,\\
\ge \, & \frac{\delta}{4\max\{s,w\}^{1/2}}\abs{s-w}\qquad \text{since $\lambda\le 1/2$ and $0 \le \theta\le \max\{s,w\}$}\,,\\
\ge \, & \frac{\delta}{4(3\lambda T)^{1/2}}\abs{s-w} \qquad\qquad \text{since $s,w \le 3 \lambda T$}\,.
\end{align*}

If $s,w \ge \cpt T$ (i.e.~right side of $\cpt T$), then we obtain similarly
$$
\csm_{(0,T]}^f(\min\{s,w\}) - \csm_{(0,T]}^f(\max\{s,w\}) \ge \frac{\delta}{6(3\lambda T)^{1/2}}\abs{s-w}\,.
$$
Combining both cases, we can show, when $\abs{s-w} > 36\sqrt{2\lambda T \log T}/\delta$,
\begin{multline*}
\qquad\abs{\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)} \ge \frac{\delta}{6(3\lambda T)^{1/2}}\abs{s-w} \\
> 2 \sqrt{6 \log T} \ge \abs{ \csm_{(0,T]}^\xi(s) - \csm_{(0,T]}^\xi(w)}\qquad
\end{multline*}
and such a difference in $\csm_{(0,T]}^f$ has the same sign as that in $\csm_{(0,T]}^X$. This implies that the gain value at the probe point that is closer to $\cpt T$ is higher than that at the other probe point. Thus, in Algorithm~\ref{Alg:OS}, we do not make a mistake in which part to discard and $\cpt T \in I_{k+1}$\,.

We consider even larger $k$ such that $\cpt T \in I_k$ and $\abs{s-w} \le 36\sqrt{2\lambda T \log T}/\delta$ with $s,w$ the probe points of $I_k$\,. For probe points being this close, the CUSUM values of the noise become highly correlated (cf. Lemma~\ref{le:cor}), and this allows for a more subtle difference in the CUSUM values of the signal to be seen from that corresponding to noisy observations. This idea is made precise in the following. By the assumption~\eqref{eq:area} and \eqref{eq:nc0}, we have
$$
\min\{s,w\} \ge \lambda T - 3\abs{s- w} \ge \lambda T - 108\delta^{-1}\sqrt{2\lambda T \log T} \ge 2^{-1/2}\lambda T\,,
$$
and
$$
T-\max\{s,w\} \ge T-(\lambda T + 3\abs{s-w}) \ge \bigl(1-\frac{5}{4}\lambda\bigr)T \ge 2^{-1/2}\lambda T\,,
$$
with the last inequality due to $\lambda \le 1/2$\,.
Then, the event $\mathcal{A}$ guarantees
$$
\abs{\csm_{(0,T]}^{\xi}(s) - \csm_{(0,T]}^{\xi}(w)} \le 6 \cdot \sqrt{\frac{2\abs{s-w} \log T}{\lambda T}}\,.
$$
Recall that when $s,w$ lie on the same side of $\cpt T$ (which are the only cases where we might wrongly drop one part of the search interval), we have
$$
\abs{\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)}\ge \frac{\delta}{6(3\lambda T)^{1/2}}\abs{s-w}\,.
$$
Thus, when $\abs{s-w} >  6^5\log T/\delta^2$, we have
\begin{multline*}
\abs{\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)}\ge \frac{\delta}{6(3\lambda T)^{1/2}}\abs{s-w} \\
> 6 \cdot \sqrt{\frac{2\abs{s-w} \log T}{\lambda T}} \ge \abs{\csm_{(0,T]}^{\xi}(s) - \csm_{(0,T]}^{\xi}(w)}\,.
\end{multline*}
Then, values of ${\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)}$ and  ${\csm_{(0,T]}^X(s) - \csm_{(0,T]}^X(w)}$ have the same sign. Thus, as argued earlier, we make no mistake in  Algorithm~\ref{Alg:OS} and thus $\cpt T\in I_{k+1}$\,.

Recall that the localization error is no larger than the minimal length of search intervals that still contain the only true change point. Therefore, by the assumption~\eqref{eq:area} with $C_0$ in~\eqref{eq:nc0}, we obtain the localization error
$$
\abs{\hat\cpt - \cpt} \le \varepsilon_T \equiv 4 \cdot 6^5\cdot \frac{\log T}{\delta^2T} \le \frac{\lambda}{2}\,,
$$
i.e.~$C_1=C_4 = 4\cdot 6^5$ in the statement of this Theorem. We further note
\begin{equation}\label{eq:max}
    \begin{aligned}
    \abs{\csm_{(0,T]}^X(\hat\cpt T)} \ge\, & \abs{\csm_{(0,T]}^f(\hat\cpt T)} - \abs{\csm_{(0,T]}^\xi(\hat\cpt T)} \\
    \ge\, & \abs{\csm_{(0,T]}^f\bigl((\cpt - \varepsilon_T) T\bigr)} - \sqrt{6 \log T}\\
    = \,&(1-\lambda)\delta\sqrt{\frac{(\lambda - \varepsilon_T)T}{1-\lambda + \varepsilon_T}} - \sqrt{6 \log T} \\
    \ge \,& \frac{\delta}{2} \sqrt{\frac{\lambda T}{2- \lambda}} -\sqrt{6 \log T} \qquad \text{as } \varepsilon_T \le \frac{\lambda}{2} \text{ and }\lambda \le \frac{1}{2}\,,\\
    \ge \, & \frac{\delta \sqrt{2\lambda}}{2} \cdot \sqrt{\frac{\lambda T}{2}} -\sqrt{6 \log T} \qquad \text{as } 0 \le \lambda \le \frac{1}{2}\,,\\
    \ge \, & \frac{145}{292}\delta \lambda \sqrt{T}\qquad \quad \qquad\, \text{by assumption \eqref{eq:area} and \eqref{eq:nc0}}.
\end{aligned}
\end{equation}

As in \citet[Lemma 2.8]{ven92}, we obtain
\begin{multline*}
   \max\left\{
\max_{t \in (0, \hat\cpt T]\cap \mathbb{N}} \abs{\csm_{(0,\hat\cpt T]}^X(t)}, \max_{t \in (\hat\cpt T, T]\cap \mathbb{N}} \abs{\csm_{(\hat\cpt T, T]}^X(t)}\right\} \\
\le 2 \delta \sqrt{\varepsilon_T T} + \sqrt{6 \log T} = 145\sqrt{6\log T}\,.
\end{multline*}
Thus, if the threshold $\gamma$ is chosen such that
$$
145\sqrt{6\log T} \le \gamma \le \frac{145}{292} \delta \lambda \sqrt{T}\,,
$$
i.e., $C_2 = 145\sqrt{6}$ and $C_3 = 145/292$ in the statement of this Theorem, then
$$
\hat \ncp = \ncp = 1\quad \text{ and }\quad \abs{\hat\cpt - \cpt} \le \varepsilon_T \qquad \text{on the event $\mathcal{A}$\,.}
$$
The assertion of part~ii.~follows since $\lim_{T \to \infty}\Prob{\mathcal{A}} = 1$.
\end{proof}

\begin{proof}[Proof of Theorem~\ref{th:aoss}]
As argued in the proof of Theorem~\ref{th:noss}, we only need to prove part ii., and can w.l.o.g.~assume
$$
f(x) = \delta\ind_{(0, \lambda]}(x) \qquad x \in (0,1]\,.
$$
For ease of notation, we further assume $\nu = 1/2$ and that there is no rounding in determining the dyadic search locations and the probe points in all search intervals. We treat advanced and combined optimistic searches separately.

\emph{Advanced optimistic search.}
We introduce the same event $\mathcal{A}$ as in~\eqref{eq:good}.
Since $\lim_{T \to \infty} \Prob{\mathcal{A}} = 1$ by Lemmata~\ref{le:max} and~\ref{le:cor}, we only need to show the localization error and the correct estimation of the number of change points on the event $\mathcal{A}$. To this end, we always assume $\mathcal{A}$ to hold below. Let constant $C_0$ in  assumption~\eqref{eq:best} satisfy
\begin{equation}\label{eq:ac0}
C_0 \ge 2^6\cdot 3^2\cdot\sqrt{6} + 2^3\cdot \sqrt{3}\,.
\end{equation}

Let us start with the dyadic search on $(0,T]$, i.e.~lines 2--4 in Algorithm~\ref{Alg:aOS}. Since $\cpt\equiv \lambda \le 1/2$ by assumption, there exists an integer $k_0 \ge 1$ such that $\cpt \in (2^{-k_0-1}, 2^{-k_0}]$\,. We claim
$$
s_* \in S \cap [2^{-k_0-1}T, 2^{-k_0}T] \equiv \{2^{-k_0-1}T, 2^{-k_0}T\}\,,
$$
where $S$  and $s_*$ are defined in lines 3 and 4 of Algorithm~\ref{Alg:aOS}, respectively.

In fact, for $s < 2^{-k_0-1} T$ and $s \in S$, it holds
\begin{align*}
& \csm_{(0,T]}^f(2^{-k_0-1}T) - \csm_{(0,T]}^f(s) \\
\ge \, & \csm_{(0,T]}^f(2^{-k_0-1}T) - \csm_{(0,T]}^f(2^{-k_0-2}T) \\
=\, & (1- \lambda) \delta \sqrt{T} \left(\sqrt{\frac{2^{-k_0-1}}{1-2^{-k_0-1}}} - \sqrt{\frac{2^{-k_0-2}}{1-2^{-k_0-2}}}\right) \\
= \, & \frac{(1-\lambda)\delta\sqrt{T}}{2^{k_0+3}\theta^{1/2}(1-\theta)^{3/2}}\qquad \text{for some } \theta \in [2^{-k_0-2}, 2^{-k_0-1}]\,.
\end{align*}
The last equality above is due to the mean value theorem. Since $2^{k_0} \le \lambda^{-1} $, $1-\lambda \ge 2^{-1}$, $1-\theta \le 1$ and $\theta \le \lambda$, we further obtain
$$
\csm_{(0,T]}^f(2^{-k_0-1}T) - \csm_{(0,T]}^f(s) \ge \frac{2^{-1} \delta \sqrt{T}}{\lambda^{-1}\cdot 2^3  \cdot\lambda^{1/2} \cdot 1^{3/2}} =  2^{-4}\delta\sqrt{\lambda T} \,.
$$

For $s > 2^{-k_0} T$ and $s \in S$, by similar calculations, we obtain
\begin{align*}
& \csm_{(0,T]}^f(2^{-k_0}T) - \csm_{(0,T]}^f(s) \\
\ge \, & \csm_{(0,T]}^f(2^{-k_0}T) - \csm_{(0,T]}^f\bigl((2^{-k_0}+2^{-k_0-1})T\bigr) \\
= \, & \frac{\lambda\delta\sqrt{T}}{2^{k_0+2}\theta^{3/2}(1-\theta)^{1/2}}\qquad \text{for some }\theta \in [2^{-k_0},\,3\cdot2^{-k_0-1}]\,, \\
\ge \, & \frac{\lambda^2 \delta \sqrt{T}}{4\theta^{3/2}} \qquad \qquad\qquad \quad\text{ since } 2^{k_0} \le \lambda^{-1} \text{ and } \theta \ge 0\,,\\
\ge \,& \frac{\delta\sqrt{\lambda T}}{12\sqrt{3}} \qquad \qquad\qquad \qquad\text{since } \theta \le 3 \cdot 2^{-k_0-1} \le 3 \lambda\,.
\end{align*}

Thus, combining both cases, we can show that on the event $\mathcal{A}$,
\begin{align*}
& \max_{s \in \{2^{-k_0-1}T,2^{-k_0}T\}}\abs{\csm_{(0,T]}^X(s)}\,\, - \max_{s \in S\, \setminus\, \{2^{-k_0-1}T,2^{-k_0}T\}}\abs{\csm_{(0,T]}^X(s)} \\
\ge \, & \max_{s \in \{2^{-k_0-1}T,2^{-k_0}T\}}\csm_{(0,T]}^f(s) \,\,- \max_{s \in S\, \setminus\, \{2^{-k_0-1}T,2^{-k_0}T\}}\csm_{(0,T]}^f(s) - 2\sqrt{6 \log T} \\
\ge \, & \min\Bigl\{\csm_{(0,T]}^f(2^{-k_0}T) \,\, - \max_{s\in S,\, s> 2^{-k_0}T}\csm_{(0,T]}^f(s),\\
& \qquad{}\qquad{}\qquad{}
\csm_{(0,T]}^f(2^{-k_0-1}T)\,\,
- \max_{s\in S,\, s < 2^{-k_0-1}T}\csm_{(0,T]}^f(s)\Bigr\} - 2\sqrt{6 \log T}\\
\ge \, & \frac{\delta \sqrt{\lambda T}}{12\sqrt{3}} - 2\sqrt{6 \log T}\\
\ge \, & 43\sqrt{\log T} > 0
\end{align*}
because of \eqref{eq:best} and \eqref{eq:ac0}. This proves our claim.

Since $s_* \le 2^{-k_0}T \le T/2$, the naive optimistic search (line 9 in Algorithm~\ref{Alg:aOS}) starts with the search interval $I_1 \equiv (s_*/2,\, 2s_*]$\,. We number naturally the later search intervals as $I_2, I_3, \ldots$ such that $I_1 \supsetneqq I_2 \supsetneqq I_3 \supsetneqq \cdots$\,. Note that $\cpt T \in (2^{-k_0-1}T,\,2^{-k_0}T] \subseteq (s_*/2,\, 2s_*] \equiv I_1$\,. In order to avoid wrongly excluding the true change point in later steps, we have to ensure that the difference of $\csm_{(0,T]}^f(\cdot)$ at two probe points is larger than the oscillation ($\asymp \sqrt{\log T}$) caused by the noise when both probe points lie on the same side of the change point.

Assume now that $\cpt T \in I_k$ for some $k \ge 1$, and let $s$ and $w$ be the two probe points of $I_k$. Note that $s,w \ge s_*/2 \ge \lambda T/4\,$. By the mean value theorem and similar calculations as above we obtain, in case of $s,w \le \cpt T$,
\begin{align*}
& \abs{\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)} \\
= \, & \frac{(1-\lambda)\delta T^{3/2}}{2\theta^{1/2}(T-\theta)^{3/2}}\abs{s- w}\qquad \text{for some $\theta \in [\min\{s,w\}, \max\{s,w\}]$}\,, \\
\ge \, & \frac{\delta \abs{s-w}}{4\sqrt{\lambda T}}\,,
\end{align*}
and in case of $s,w \ge \cpt T$,
\begin{align*}
& \abs{\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)} \\
= \, & \frac{\delta \lambda T^{3/2}}{2\theta^{3/2}(T - \theta)^{1/2}}\abs{s-w}\qquad \text{for some $\theta \in [\min\{s,w\}, \max\{s,w\}]$}\,, \\
\ge \, & \frac{\delta \abs{s-w}}{6 \sqrt{3\lambda T}}\,.
\end{align*}
Note further that $T - \max\{s,w\} \ge T - 3s_*/2 \ge s_*/2 \ge \lambda T/4\,$. Then, on the event $\mathcal{A}$, it holds
$$
\abs{\csm_{(0,T]}^\xi(s) - \csm_{(0,T]}^\xi(w)} \le 6\cdot\sqrt{\frac{\abs{s-w}\log T }{\lambda T /4}} = 12 \cdot \sqrt{\frac{\abs{s-w}\log T}{\lambda T}}\,.
$$
Therefore, values of $\csm_{(0,T]}^f(s) - \csm_{(0,T]}^f(w)$ and $\csm_{(0,T]}^X(s) - \csm_{(0,T]}^X(w)$ have the same sign (and thus $\cpt T \in I_{k+1}$) provided that
$$
\frac{\delta\abs{s-w}}{6\sqrt{3\lambda T}} > 12 \cdot \sqrt{\frac{\abs{s-w}\log T}{\lambda T}}\,.
$$
This condition is equivalent to
$$
\abs{s- w} >  3^5\cdot2^6\cdot \frac{\log T}{\delta^2}\,.
$$
Recall that the localization error is no larger than the minimal length of search intervals that still contain the only true change point. Thus, on the event $\mathcal{A}$, by~\eqref{eq:best} and~\eqref{eq:ac0}, we have the localization error
$$
\abs{\hat \cpt - \cpt} \le \varepsilon_T \equiv 3^5\cdot 2^8\cdot \frac{\log T}{\delta^2 T} \le \frac{\lambda}{2} \,,
$$
i.e.~$C_1 = C_4 = 3^5\cdot 2^8$ in the statement of this Theorem.
Further, by similar calculations of \eqref{eq:max} in the proof of Theorem~\ref{th:noss}, we obtain on the event $\mathcal{A}$
\begin{align*}
\abs{\csm_{(0,T]}^X(\hat\cpt T)} \ge\, & \abs{\csm_{(0,T]}^f(\hat\cpt T)} - \abs{\csm_{(0,T]}^\xi(\hat\cpt T)} \\
    \ge\, & (1-\lambda)\delta\sqrt{\frac{(\lambda - \varepsilon_T)T}{1-\lambda + \varepsilon_T}} - \sqrt{6 \log T} \\
    \ge \,& \frac{\delta\sqrt{\lambda T}}{2\sqrt{2}}- \sqrt{6 \log T}\qquad \text{as } \varepsilon_T \le \frac{\lambda}{2} \text{ and } 0 \le \lambda \le \frac{1}{2}\,,\\
    \ge \, & \frac{3^2\cdot 2^4 \cdot\sqrt{2} + 1}{4\sqrt{2} + 2^6\cdot 3^2}\delta  \sqrt{\lambda T}\qquad \quad \qquad\, \text{by assumption \eqref{eq:best} and \eqref{eq:ac0}}.
\end{align*}
Similar as in \citet[Lemma 2.8]{ven92}, we can show, on the event $\mathcal{A}$,
\begin{multline*}
\max\left\{
\max_{t \in (0, \hat\cpt T]\cap \mathbb{N}} \abs{\csm_{(0,\hat\cpt T]}^X(t)}, \max_{t \in (\hat\cpt T, T]\cap \mathbb{N}} \abs{\csm_{(\hat\cpt T, T]}^X(t)}\right\} \\
\le 2 \delta \sqrt{\varepsilon_T T} + \sqrt{6 \log T} = (3^2\cdot 2^5\cdot\sqrt{3}+\sqrt{6})\sqrt{\log T}\,.
\end{multline*}
Therefore, on the event $\mathcal{A}$,
$$
\hat \ncp = \ncp = 1\quad \text{ and }\quad \abs{\hat\cpt - \cpt} \le \varepsilon_T
$$
if the threshold $\gamma$ is chosen such that
$$
(3^2\cdot2^5\cdot\sqrt{3}+\sqrt{6})\sqrt{\log T} \le \gamma \le \frac{3^2\cdot2^4\cdot\sqrt{2}+1}{4\sqrt{2}+2^6\cdot 3^2} \delta \sqrt{\lambda T}\,,
$$
i.e.~$C_2 = 3^2\cdot2^5\cdot\sqrt{3} + \sqrt{6}$ and $C_3 = (3^2\cdot2^4\cdot\sqrt{2} + 1)/(4\sqrt{2}+2^6\cdot 3^2)$ in the statement of this Theorem.
Recall that  $\lim_{T \to \infty}\Prob{\mathcal{A}} = 1$, so the assertion of advanced optimistic search is proven.

\emph{Combined optimistic search.} Let constant $C_0$ in  assumption~\eqref{eq:best} satisfy
\begin{equation}\label{eq:cc0}
C_0 \ge 2^6\cdot3^3\cdot\sqrt{2} + 2^3\cdot\sqrt{3}\,.
\end{equation}
We consider the same event $\mathcal{A}$, cf.~\eqref{eq:good}, as for the advanced optimistic search, and always assume $\mathcal{A}$ to occur. It has been  shown above that the estimated change point $\hat \cpt_{\mathrm{a}}$ by the advanced optimistic search satisfies
$$
\abs{\hat \cpt_{\mathrm{a}} - \cpt} \le \varepsilon_T \equiv 3^5\cdot 2^8\cdot \frac{\log T}{\delta^2 T}\qquad \text{on the event $\mathcal{A}$}\,.
$$
We claim now on the event $\mathcal{A}$,
\begin{equation}\label{eq:ball}
\min_{\abs{s/T - \cpt}\le\varepsilon_T}\abs{\csm_{(0,T]}^X(s)} \quad > \quad \max_{\abs{s/T - \cpt}\ge3\varepsilon_T}\abs{\csm_{(0,T]}^X(s)}\,.
\end{equation}
This claim~\eqref{eq:ball} will be proven in two steps.

Firstly, we prove an easier version of~\eqref{eq:ball}, which is on the event $\mathcal{A}$
\begin{equation*
\min_{\abs{s/T - \cpt}\le\varepsilon_T}\abs{\csm_{(0,T]}^X(s)} \quad > \quad \max_{\abs{s/T - \cpt}\ge\lambda/2}\abs{\csm_{(0,T]}^X(s)}\,.
\end{equation*}
Note that  $\varepsilon_T \le \lambda/6$ by assumption~\eqref{eq:best} with $C_0$ in \eqref{eq:cc0}.
Then
\begin{align*}
& \csm_{(0,T]}^f\Bigl((\lambda- \varepsilon_T) T \Bigr) - \csm_{(0,T]}^f\Bigl(\frac{3}{2}\lambda T\Bigr) \\
= \, & \delta (1- \lambda) \sqrt{\frac{(\lambda-\varepsilon_T)T}{1-\lambda + \varepsilon_T}} - \delta \lambda \sqrt{\frac{(1-3\lambda/2)T}{3\lambda/2}} \\
\ge \,& \delta (1- \lambda) \sqrt{\frac{5\lambda T/6}{1-5\lambda/6}} - \delta \lambda \sqrt{\frac{(1-3\lambda/2)T}{3\lambda/2}}\\
= \, & \frac{\delta (3- 2 \lambda) \sqrt{\lambda T}}{\sqrt{3(6-5\lambda)}\left((1-\lambda)\sqrt{15} + \sqrt{(2-3\lambda)(6-5\lambda)}\right)} \\
\ge \, & \frac{2\delta\sqrt{\lambda T}}{3\sqrt{30}+6\sqrt{6}} \,.
\end{align*}
Moreover, since the value of $\csm_{(0,T]}^f(\cdot)$ changes faster on $(0,\,\cpt T]$ than on $(\cpt T,\, T]$, we obtain, on the event $\mathcal{A}$,
\begin{align*}
&\min_{\abs{s/T - \cpt}\le\varepsilon_T}\abs{\csm_{(0,T]}^X(s)}  - \max_{\abs{s/T - \cpt}\ge\lambda/2}\abs{\csm_{(0,T]}^X(s)}\\
\ge \,& \min_{\abs{s/T - \cpt}\le\varepsilon_T}\abs{\csm_{(0,T]}^f(s)}  - \max_{\abs{s/T - \cpt}\ge\lambda/2}\abs{\csm_{(0,T]}^f(s)} - 2\sqrt{6 \log T} \\
=\,& \csm_{(0,T]}^f\Bigl((\lambda- \varepsilon_T) T \Bigr) - \csm_{(0,T]}^f\Bigl(\frac{3}{2}\lambda T\Bigr) - 2\sqrt{6 \log T} \\
\ge \,& \frac{2\delta\sqrt{\lambda T}}{3\sqrt{30} +6\sqrt{6}} - 2\sqrt{6 \log T} \\
\ge \, & 153 \sqrt{\log T} > 0,
\end{align*}
where the third inequality is due to assumption~\eqref{eq:best} with $C_0$ in~\eqref{eq:cc0}.

Secondly, we still need to show, on the event $\mathcal{A}$,
$$
\min_{\abs{s/T - \cpt}\le\varepsilon_T}\abs{\csm_{(0,T]}^X(s)} \quad > \quad \max_{3\varepsilon_T \le \abs{s/T - \cpt}\le\lambda/2}\abs{\csm_{(0,T]}^X(s)}\,.
$$
To that end, we consider arbitrarily $s_1 , s_2 \in \mathbb{N}$ such that $\abs{s_1/T - \cpt}\le\varepsilon_T$ and $3\varepsilon_T \le \abs{s_2/T - \cpt}\le\lambda/2$. Then, on the event $\mathcal{A}$,
\begin{equation}\label{eq:diffcsm}
\begin{aligned}
&\abs{\csm_{(0,T]}^X(s_1)} - \abs{\csm_{(0,T]}^X(s_2)} \\
= \, &\csm_{(0,T]}^X(s_1) - \csm_{(0,T]}^X(s_2)  \\
\ge\,& \csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2) - \abs{\csm_{(0,T]}^\xi(s_1) -\csm_{(0,T]}^\xi(s_2)} \\
\ge \, &\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2) - 6 \sqrt{\frac{2\abs{s_1-s_2}\log T}{\lambda T}}\,.
\end{aligned}
\end{equation}
In order to bound $\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2)$ from below, we deal with four cases separately. (Recall that $\lambda \equiv \cpt$\,.)
\begin{enumerate}[i.]
\item In case of $s_2 \le s_1 \le \lambda T$, we obtain by the mean value theorem
\begin{align*}
&\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2)\\
= \, & \delta(1-\lambda)\left(\sqrt{\frac{s_1 T}{T - s_1}}-\sqrt{\frac{s_2 T}{T - s_2}}\right)\\
= \,& \delta(1-\lambda)T^{3/2}\frac{s_1 - s_2}{2\theta^{1/2}(T - \theta)^{3/2}}\quad\qquad \text{for some }\theta \in [s_1, s_2]\,,\\
\ge \, &\frac{\delta \abs{s_1 - s_2}}{4\sqrt{\lambda T}}\,.
\end{align*}
\item In case of $\lambda T \le s_1 \le s_2$, similar calculation as in case i.~leads to
$$
\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2) \ge \frac{\sqrt{6}\delta \abs{s_1-s_2}}{9 \sqrt{\lambda T}}\,.
$$
\item For $s_2 \le \lambda T \le s_1$, we define $s_1' = 2\lambda T -s_1$ then $s_2 \le s_1' \le \lambda T$ and $\abs{s_1' - s_2} \ge \abs{s_1 - s_2}/2$. Thus,
\begin{align*}
&\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2)\\
\ge \,& \csm_{(0,T]}^f(s_1')-\csm_{(0,T]}^f(s_2) \\
\ge \, & \frac{\delta\abs{s_1' - s_2}}{4\sqrt{\lambda T}} \quad\qquad \qquad \qquad \text{by case i.,} \\
\ge \, & \frac{\delta\abs{s_1-s_2}}{8\sqrt{\lambda T}}\,.
\end{align*}

\item In case of $s_1 \le \lambda T \le s_2$, which is the most difficult one, we obtain
\begin{align*}
&\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2)\\
= \, & \delta (1-\lambda)\sqrt{\frac{s_1 T}{T - s_1}} - \delta \lambda \sqrt{\frac{(T - s_2)T}{s_2}} \\
= \, & \frac{\delta \sqrt{T}\left((1-2\lambda) s_1 s_2 - \lambda^2 T^2 + \lambda^2 (s_1 + s_2)T\right)}{\sqrt{s_2(T-s_1)}\left((1-\lambda)\sqrt{s_1s_2} + \lambda \sqrt{(T-s_1)(T-s_2)}\right)}\,.
\end{align*}
Since $1- \lambda \le 1$, $s_1 \le \lambda T$, $s_2 \le 3\lambda T/2$, $T-s_1 \le T$ and $T-s_2 \le T$, then
\begin{align*}
&\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2)\\
\ge \, & \frac{2\delta}{(3+\sqrt{6})(\lambda T)^{3/2}}\left((1-2\lambda) s_1 s_2 - \lambda^2 T^2 + \lambda^2 (s_1 + s_2)T\right) \\
= \, & \frac{2\delta}{(3+\sqrt{6})(\lambda T)^{3/2}} \bigl(\lambda(1-\lambda)(s_1+s_2-2\lambda T)T -(1-2\lambda)(\lambda T - s_1)(s_2- \lambda T)\bigr) \\
\ge \, & \frac{2\delta}{(3+\sqrt{6})(\lambda T)^{3/2}} \left(\lambda(1-\lambda)(s_1+s_2-2\lambda T)T-(1-\lambda)(\lambda T -s_1)\frac{\lambda T}{2}\right) \\
= \, & \frac{2\delta(1-\lambda)}{(3 + \sqrt{6})\sqrt{\lambda T}}\left((s_2 - \lambda T) - \frac{3}{2}(\lambda T -s_1)\right)\,.
\end{align*}
Because  further $s_2 -\lambda T \ge 3\varepsilon_T T$ and $\lambda T -s_1 \le \varepsilon_T T$ we have
$$
(s_2 - \lambda T) -\frac{3}{2}(\lambda T -s_1) \ge \frac{3\bigl((s_2 - \lambda T) + (\lambda T -s_1)\bigr)}{8} = \frac{3(s_2 - s_1)}{8}\,.
$$
Thus, by $1-\lambda \ge 1/2$, it further follows
$$
\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2) \ge \frac{3\delta \abs{s_1-s_2}}{8(3+\sqrt{6})\sqrt{\lambda T}}\,.
$$
\end{enumerate}
Combining all four cases, we obtain
\begin{align*}
&\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2) - 6 \sqrt{\frac{2\abs{s_1-s_2}\log T}{\lambda T}} \\
\ge\, & \frac{3\delta \abs{s_1-s_2}}{8(3+\sqrt{6})\sqrt{\lambda T}} - 6 \sqrt{\frac{2\abs{s_1-s_2}\log T}{\lambda T}} \\
= \, & \frac{3\delta\sqrt{\abs{s_1 - s_2}}}{8(3+\sqrt{6})\sqrt{\lambda T}}\left(\sqrt{\abs{s_1 - s_2}} - \frac{16(3\sqrt{2} + 2\sqrt{3})\sqrt{\log T}}{\delta}\right)\,.
\end{align*}
By $\abs{s_1 - s_2} \ge 2\varepsilon_T T$ and the definition of $\varepsilon_T$, we have, on the event $\mathcal{A}$,
\begin{equation*}
\csm_{(0,T]}^f(s_1)-\csm_{(0,T]}^f(s_2) - 6 \sqrt{\frac{2\abs{s_1-s_2}\log T}{\lambda T}}  \ge 22 \sqrt{\frac{\varepsilon_T \log T}{\lambda}} > 0\,.
\end{equation*}
Thus, due to~\eqref{eq:diffcsm} and the arbitrariness of $s_1$ and $s_2$, we have concluded the second step, and thus proven the claim~\eqref{eq:ball}.

By $\hat \cpt_{\mathrm{c}}$ we denote the estimated change point by the combined optimistic search. By construction, we have
$$
\abs{\csm_{(0,T]}^X(\hat \cpt_\mathrm{c} T)} \ge \abs{\csm_{(0,T]}^X(\hat \cpt_\mathrm{a} T)}\quad \text{and}\quad \abs{\hat \cpt_{\mathrm{a}} - \cpt}\le \varepsilon_T\,.
$$
Thus, by the claim~\eqref{eq:ball} it holds that $\abs{\hat \cpt_{\mathrm{c}} - \cpt} \le 3 \varepsilon_T$ on the event $\mathcal{A}$. That is, $C_1=C_4= 3^6\cdot 2^8$ in the statement of this Theorem. The rest of proof follows in exactly the same way as in that of advanced optimistic search. More precisely, we can show, if the threshold $\gamma$ is chosen as
$$
(2^5\cdot3^3 + \sqrt{6})\sqrt{\log T} \le \gamma \le \frac{2^4\cdot3^3\cdot\sqrt{2} +\sqrt{3}}{2^6\cdot3^3+4\sqrt{6}} \delta\sqrt{\lambda T}
$$
i.e.~$C_2 = 2^5\cdot3^3 + \sqrt{6}$ and $C_3 = (2^4\cdot3^3\cdot\sqrt{2} +\sqrt{3})/(2^6\cdot3^3+4\sqrt{6})$,  then the assertion of part~ii.~holds.
\end{proof}

\subsection{Multiple change points}\label{ss:mcpt}
The intuition why the OSeedBS in combination with the NOT selection performs optimally is that the selected seeded interval often contains only a single change point. On such a seeded interval, advanced and combined optimistic searches perform in a minimax optimal way. We make this precise in the following.

\begin{proof}[Proof of Theorem~\ref{th:osbs}]
We start with part i., and let the constant $C_0$ in assumption~\eqref{eq:best} satisfy
$$
C_0 \ge 2^5\cdot 3 + 2^3\cdot 3\cdot \sqrt{2} + 2^8\cdot 3^3\cdot \sqrt{6}\,.
$$
We set the threshold $\gamma = C_1 \sqrt{\log T}$ with an arbitrary constant
$$
C_1 \ge 2^6\cdot 3^3 \sqrt{2} + 2^3\cdot\sqrt{3}+\sqrt{6}\,.
$$
For notation simplicity, we assume there is no rounding, the decay $a = 1/2$ for seeded intervals, and the step size $\nu = 1/2$ for optimistic searches (otherwise, only the multiplying constants may be different). Define the event $\mathcal{A}$ as in~\eqref{eq:good}. We split the proof into several steps.

\emph{Step 1.}
Consider intervals $\bigl((\cpt_i - \lambda/2) T,\, (\cpt_i + \lambda/2)T\bigr]$ as potential backgrounds for change points $\cpt_i T$ for $i= 1,\ldots,\ncp$\,. By the construction of seeded intervals, we can find seeded intervals $(c_i -r_i, c_i + r_i]$ such that
\begin{align*}
(c_i - r_i,\, c_i + r_i] &\subseteq \bigl((\cpt_i -\frac{1}{2}\lambda)T,\, (\cpt_i + \frac{1}{2} \lambda)T] \\
r_i &\ge \frac{1}{6} \lambda T\\
\text{and}\qquad \abs{c_i - \cpt_i T} &\le \frac{1}{2}r_i\,.
\end{align*}
Note that $\abs{(c_i - r_i,\, c_i + r_i]} \ge \lambda T/3 \ge m$ and that $(c_i - r_i,\, c_i + r_i]$ contains only a single change point $\cpt_i T$ for every $i = 1,\ldots, \ncp$\,.

\emph{Step 2.}
By $\hat \cpt_i^0$ we denote the estimated change point (scaled by $1/T$) by either advanced or combined optimistic search on $(c_i - r_i,\,c_i + r_i]$.
Since there is a single change point in $(c_i - r_i,\,c_i + r_i]$, for every $i=1,\ldots,\ncp$,  from the proof of Theorem~\ref{th:aoss}, we can easily see, on the event $\mathcal{A}$,
$$
\abs{\hat \cpt_i^0 - \cpt_i} \le \varepsilon_{i,T}^0 \equiv 3^4 \cdot 2^{10}\cdot\frac{\log T}{\delta_i^2 T}\le \frac{\lambda}{24}
$$
and further
\begin{equation}\label{eq:value}
\abs{\csm^X_{(c_i - r_i,\,c_i + r_i]}(\hat \cpt_i^0 T)} \ge \left(\frac{C_0}{4\sqrt{3}} - \sqrt{6}\right) \sqrt{\log T} \ge \gamma \equiv C_1 \sqrt{\log T}\,.
\end{equation}

\emph{Step 3.}
Recall that the NOT selection rule selects the shortest seeded interval $I_*$ which has a CUSUM value above the threshold $\gamma$. Then, on the event $\mathcal{A}$, from~\eqref{eq:value}, it follows that the length of selected seeded interval $\abs{I_*} \le \lambda$, and thus $I_*$ contains at most one change point. If $I_*$  contains no change point, then on the event $\mathcal{A}$,
$$
\max_{s \in I_*} \abs{\csm_{I_*}^X(s)} \le \sqrt{6 \log T} < \gamma \equiv  C_1 \sqrt{\log T}\,.
$$
Thus, the selected seeded interval $I_*$ contains exactly one change point, conditioned on the event $\mathcal{A}$.

\emph{Step 4.}
We rewrite the selected seeded interval $I_*$ as $\bigl((\cpt_i - u) T, \,(\cpt_i + v) T \bigr]$ for some $i \in \{1,\ldots, \ncp\}$. Then, on the event $\mathcal{A}$
\begin{multline*}
\gamma \le \max_{s \in I_*} \abs{\csm_{I_*}^X(s)} \le \max_{s \in I_*} \abs{\csm_{I_*}^f(s)} + \max_{s \in I_*} \abs{\csm_{I_*}^\xi(s)} \\
= \abs{\csm_{I_*}^f(\tau_i T)} + \max_{s \in I_*} \abs{\csm_{I_*}^\xi(s)} \le \delta_i\sqrt{\frac{u v T}{u + v}} + \sqrt{6 \log T}\,,
\end{multline*}
which further implies
\begin{multline}
\label{eq:snr}
 \qquad\delta_i\sqrt{\min\{u,v\}T} \ge \delta_i\sqrt{\frac{u v T}{u + v}} \ge (C_1 - \sqrt{6})\sqrt{\log T} \\
 \ge \left(2^6\cdot 3^3 \sqrt{2} + 2^3\cdot\sqrt{3}\right) \sqrt{\log T}\,.\qquad
\end{multline}
Let $\hat \cpt_i$ be the estimated change point (scaled by $1/T$) by advanced or combined optimistic search applied to $I_*$\,. Then, from the proof of Theorem~\ref{th:aoss}, we can easily see that \eqref{eq:snr} implies
\begin{equation*}
\abs{\hat\cpt_i -\cpt_i} \,\,\le\,\, \varepsilon_{i,T} \,\,\equiv\,\,3^6\cdot 4^4\cdot\frac{\log T}{\delta_i^2 T}\qquad \text{on the event }\mathcal{A}\,.
\end{equation*}

\emph{Step 5.}
Steps 1-4 imply that all change points $\{\cpt_1, \ldots, \cpt_{\ncp}\}$ can be estimated by $\{\hat \cpt_1, \ldots, \hat \cpt_{\ncp}\}$ with errors $\varepsilon_{i,T}$ and in particular $\hat \ncp \ge \ncp$\,. It remains to show that $\hat \ncp = \ncp$\,. To this end, we note that after the detection of $\ncp$ change points, the remaining seeded intervals either contain no change point or contain one or two change points that are very close to the boundary. In the latter case, if a seeded interval contains say $\cpt_i$ then we must have the distance of $\cpt_i$ to one of the boundaries is no more than $\varepsilon_{i,T}$\,. Similar to \citet[Lemma 2.8]{ven92}, we can show that the maximum value of CUSUM statistics $\max_{s\in I}\abs{\csm_{I}^X(s)}$ on every remaining seeded interval $I$ is upper bounded by
$$
2 \max_{i = 1,\ldots,\ncp}\delta_i\sqrt{\varepsilon_{i,T}T} + \sqrt{6 \log T} \le (3^3 \cdot 2^5+\sqrt{6})\sqrt{\log T} < \gamma \equiv C_1\sqrt{\log T}\,,
$$
on the event $\mathcal{A}$\,. This further implies that  OSeedBS will stop after $\ncp$ steps, i.e.~$\hat\ncp = \ncp$\,. Then, part i.~is proven by noticing that $\Prob{\mathcal{A}} \to 1$ as $T \to \infty$\,.

\bigskip

Now we consider part ii.~of the Theorem. By $\mathcal{I}$ denote all seeded intervals of length no less than $m \asymp T^{1-w}$.  Then the total number of evaluations in advanced or combined optimistic search on all seeded intervals in $\mathcal{I}$ is
\begin{align*}
\sum_{I \in \mathcal{I}} \log \abs{I} \lesssim\, & \sum_{k = 1}^{\ceil{\log_{\frac{1}{a}}(T^\omega)}}a^{-(k-1)} \log(T a^{k-1})\\
\lesssim\, & \int_0^{\log_{\frac{1}{a}}(T^\omega)}a^{-x} \log(Ta^x) dx\\
\lesssim \,& \int_{0}^{T^\omega} \log\left(\frac{T}{x}\right) dx = T^\omega \bigl(1 + (1-\omega)\log T\bigr) \\
\lesssim\, & \min\{T^{\omega}\log T,\, T\}\,.
\end{align*}
The NOT selection only requires further $O(\abs{\mathcal{I}})=O(T^\omega)$ computations, because of the deterministic and nested structure of the seeded intervals $\mathcal{I}$. Thus, if cumulative sums are pre-computed, the final computation complexity is $O(\min\{T^{\omega}\log T,\, T\})$; otherwise, the computation of cumulative sums dominates, and leads to the final computation complexity of $O(T)$.
\end{proof}

redpajama_set_name: RedPajamaArXiv