Detection of multiple changes in mean by sparse parameter estimation

. The contribution is focused on detection of multiple changes in the mean in a one-dimensional stochastic process by sparse parameter estimation from an overparametrized model. The authors’ approach to change point detection differs entirely from standard statistical techniques. A stochastic process residing in a bounded interval with changes in the mean is estimated using dictionary (a family of functions, the so-called atoms, which are overcomplete in the sense of being nearly linearly dependent) and consisting of Heaviside functions. Among all possible representations of the process we want to ﬁnd a sparse one utilizing a signiﬁcantly reduced number of atoms. This problem can be solved by (cid:96) 1 -minimization. The basis pursuit algorithm is used to get sparse parameter estimates. In this contribution the authors calculate empirical probability of successful change point detection as a function depending on the number of change points and the level of standard deviation of additive white noise of the stochastic process. The empirical probability was computed by simulations where locations of change points were chosen randomly from uniform distribution. The authors’ approach is compared with LASSO algorithm, (cid:96) 1 trend ﬁltering and selected statistical methods. Such probability decreases with increasing number of change points and/or standard deviation of white noise. The proposed method was applied on the time series of nuclear magnetic response during the drilling of a well.


Introduction
Detection of changes in the mean helps us to discover time points at which some intrinsic properties of a stochastic process change.This covers a broad range of real-world problems and has been actively discussed.The standard statistical approach to this problem is described, for example, in [1].The article is focused on detection of multiple changes in the mean in a one-dimensional stochastic process using sparse parameter estimation from an overparametrized model.Sparse modeling is a quickly developing area on the intersection of statistics, machine-learning and signal processing.It can be effectively used in the case of large data sets fitting when standard numerical methods are unstable.Several methods of large-scale data analysis can be found, for example, in [2].Detection of changes has originally arisen in the context of quality control.Nowadays, we can find a wide range of fields where change point problem is applied, such as epidemiology, medicine (rhythm analysis), ecology, signal processing etc.When analyzing measurements coming from real problems often some events or anomalies are manifested by abrupt changes which may not be quite apparent in their early stage, detection of hearing loss in audiological frequencies [3] being an example out of many illustrating such a situation.Detection of changes in processes need not involve only changes in the mean, but also changes in the variance or other characteristics.For example, the article [4] deals with detection of changed segment in a binomial sequence.
The authors' approach to change point detection is quite different from the standard statistical techniques.The statistical approaches usually use for change point estimation the maximum likelihood method or non-parametric methods, such as methods based on M and R estimates (see [1,5]).In the case of uncorrelated white noise in the model and one change in the mean, the authors Neubauer and Veselý [6] showed that there is no significant difference among these methods.The proposed method of the change point detection is based on entirely different methods, built on a completely different approach.It is a very versatile approach to signal modeling, assuming that the signal is sparse; therefore, it can be described by a small number of suitable functions.A stochastic process residing in a bounded interval with changes in the mean is estimated using dictionary (a family of functions, the so-called atoms) and consisting of Heaviside functions.Among all the possible representations of the process we want to find a sparse one utilizing a significantly reduced number of atoms.This problem can be solved by 1 -minimization (see [8,9]) or by some of the greedy algorithms which are compared to the 1 -minimization in [10].The basis pursuit algorithm is used to get sparse parameter estimates.Basic properties of this approach to one change point detection were studied in [6] by simulations and included a performance comparative study with standard statistical procedures as well.The impact of serial correlation of noise on change point detection is described in [7].The contribution [11] deals primarily with detection of two change points.According to these results, the basis pursuit approach proposes an alternative technique of the change point detection.In this article the authors calculate empirical probability of successful change point detection as a function depending on the number of change points and the level of standard deviation of additive white noise of the stochastic process.The empirical probability was computed by simulations where locations of change points were chosen randomly from uniform distribution.As an alternative to the basis pursuit algorithm, LASSO method [12], 1 trend filtering [13] and selected statistical method were taken to judge the proposed method.The last part of the article demonstrates how the aforementioned methods can automate the change point detection in the time series of nuclear magnetic response during the drilling of a well.

Basis pursuit and change point detection
In this paragraph we briefly describe the method based on basis pursuit algorithm (BPA) for the detection of the change point in the sample path {y t } in one-dimensional stochastic process {Y t }.We assume a deterministic functional model on a bounded interval I described by the dictionary G = {G j } j∈J with atoms G j ∈ L 2 (I) and with additive white noise e on a suitable finite discrete mesh T ⊂ I: where x ∈ sp({G j } j∈J ), e t ∼ WN (0, σ 2 ), t ∈ T , σ > 0, and J is a big finite indexing set.Smoothed function x = j∈J ξj G j =: G ξ minimizes on T 1 -penalized optimality measure (1/2) y − Gξ 2 2 as follows: where λ = σ 2 ln (cardJ) is a smoothing parameter chosen according to the softthresholding rule commonly used in wavelet theory.This choice is natural because one can prove that with any orthonormal basis G = {G j } j∈J the shrinkage via soft-thresholding produces the same smoothing result x (see [9]).Such approaches are also known as basis pursuit denoising (BPDN).Solution of this minimization problem with λ close to zero may not be sparse enough: we are searching small F ⊂ J such that x ≈ j∈F ξj G j is a good approximation.The procedure of BPDN is described in [6].
We build our dictionary from heaviside-shaped atoms on L 2 (R) derived from a fixed "mother function" via shifting and scaling following the analogy with the construction of wavelet bases.
We construct an oversized shift-scale dictionary G = {G a,b } a∈A,b∈B derived from the "mother function" by varying the shift parameter a and the scale (width) parameter b between values from big finite sets A ⊂ R and B ⊂ R + , respectively (J = A × B), on a bounded interval I ⊂ R spanning the space H = sp({G a,b }) a∈A,b∈B , where In the simulations below I = [0, 1], T = {t/T } (typically with mesh size T = 100), A = {t/T } T −t0 t=t0 (t 0 is a boundary trimming, t 0 = 5 was used in the simulations) and scale b fixed to zero (B = {0}).Clearly the atoms of such Heaviside dictionary are normalized on I, i.e.G a,0 2 = 1.Some examples of Heaviside functions are displayed in Fig. 1.Neubauer and Veselý [6] proposed the method of change point detection if there is just one change point in a one-dimensional stochastic process (or in its sample path).We briefly describe this method.We would like to find a change point in a stochastic process where µ, δ = 0, t 0 c < T − t 0 are unknown parameters and t are independent identically distributed random variables with zero mean and variance σ 2 .The parameter c indicates the change point in the process.Using the basis pursuit algorithm we obtain some significant atoms, we calculate the correlation between the significant atoms and the analyzed process.As a BPDN estimator of the change point c we take the shift parameter of the significant atom matching one of three alternative rules: the atom is the most significant one, it exhibits the highest correlation with the data or it produces the highest value of geometric mean of the correlation coefficient and its contribution to the total 1 norm.In the case of one change point, the estimator based on the highest correlation performed slightly better.Now let us assume the model with p change points where µ, δ 1 , . . ., δ p = 0, t 0 c 1 < • • • < c p < T − t 0 are unknown parameters and t are independent identically distributed random variables with zero mean and variance σ 2 .
We use the method of change point estimation described above for detection of p change points c 1 , . . ., c p in model (2).Instead of finding only one significant atom we identify either p most significant atoms, or p significant atoms exhibiting the highest correlation, or p significant atoms having the highest value of geometric mean of the correlation coefficient and the contribution to the total 1 norm.The shift parameters of such atoms determine direct BPDN estimators for the change points c 1 , . . ., c p .Another possibility is to apply the procedure of one change point detection p-times in sequence (sequential BPDN estimator).In the first step we identify one change point in the process Y t , then we subtract the contribution of such a significant atom from the process (by linear regression) and we apply the same method to the new process Y t .This sequence is repeated until we get p change point estimates.The shift parameters of selected atoms are again identifiers of the change points c 1 , . . ., c p .
Regarding the theoretical background, we consider three well-known criteria.
The so-called spark of a matrix (dictionary) G is the minimum number of columns of G which are linearly dependent.It can be shown (see [14]) that if a solution x of our inverse problem is found which satisfies x 0 < spark(G)/2 then x is the sparsest and unique solution to this problem.However, in general, establishing spark is an NP-hard problem so in practice this is intractable.Nevertheless, due to a clear structure of our dictionary G, we know that its spark is full, i.e. 101 in our case.Thus, if a solution x would be found with x 0 45 then we have the true solution.This actually motivates to exploit the sparsity-seeking approach.However, one must be aware of two major facts: First, this result only holds in case without noise and the noisy case is not solved in the literature yet.Second, the statement does not tell anything about how the solution should be obtained!Another measure of goodness of the dictionary is the mutual coherence (see [8,14]) which is defined pairwise as µ(G) = max i =j g T i g j for 2 -normalized columns.This value can be determined computationally -in our case it is 0.99 which means that the atoms are very close to each other.The relationship between µ and spark is known and thanks to that there exists a theorem saying that if a solution x is found satisfying x 0 < (1/2)(1 + 1/µ(G)) then it is the unique one.But in our case the value of the right side only 1.0051 which means that this would be true only if the sparsity of the solution would be 1 (without noise, again).In our case of multiple change-points this is contradictory, so the coherence does not help much here.
The third property of a dictionary is the RIP (restricted isometry property, see [15]) which is again not possible to compute in polynomial time, but can be used to specify a condition under which the solution can be found via 1 -minimization (which we use).We found numerically that the RIP constant δ k which could be found in a combinatorial way even for such low values as k = 2 and k = 3 is in the order of 10 1 ; RIP requires δ k < 1. However this is not a surprise, since the only matrices which are known to obey the RIP are certain classes of random matrices, see [8,15].
To conclude, no one of the state-of-the art quantities measuring the goodness of the dictionary really helps, so we rely on the simulation results (which is in fact not uncommon in the sparse-land publications area), see [16,17].

The LASSO
The LASSO (least absolute selection and shrinkage operator) is a shrinkage and selection method for linear regression.It minimizes the sum of squared errors, with a bound on the sum of the absolute values of the coefficients.Using the usual linear regression notation the LASSO estimate is defined in [12] as follows: Rewriting this in terms of operators and language of atoms and dictionaries, we obtain a denoising approach similar to basis pursuit.

1 trend filtering
In [13] the authors propose a variation on Hodrick-Prescott filtering, which they call 1 trend filtering.The trend is estimated as the minimizer of the objective function which can be written in the matrix form as where χ 0 is a smoothing parameter, x = (x 1 , . . ., x T ) ∈ R T , y = (y 1 , . . ., y n ) ∈ R T and D ∈ R (T −2)×T is the matrix The 1 trend estimate is piecewise linear in t.The points where the slope of the estimated trend is changed can be interpreted as abrupt changes (change points) in the process.Model (2) describes only changes in the mean.We adapt the proposed method for detection change point of this kind.The argument appearing in the second term of (3), It is zero when and only when three points x t−1 , 2x t and x t+1 are on the line.In the case of the piecewise constant trend, we replace the second difference with the first difference x t − x t−1 .We will minimize the function in the matrix form (1/2) y − x 2 2 + χ Dx 1 , where D ∈ R (T −1)×T is the matrix .
The estimated trend is piecewise constant.The points where the estimated trend is changed can be taken as estimates of change points.

Statistical methods
Consider model (1).Assuming σ 2 given, the unknown parameters c, µ and δ may be estimated by the least-squares method.The least-squares estimators ĉ, μ and δ of the parameters c, µ and δ are defined as solutions of the minimization problem In other words, the unknown parameters are estimated in such a way that the sum of squares of residuals is minimal.The estimates of the parameters µ and δ are (see [5] etc.) where ĉ is a solution of the maximization problem where we denote In the case of multiple change points, we can use the described statistical method of one change point detection via the following procedure.At first, find ĉ(1) by solving (5).Secondly, divide observations into two groups Y 1 , . . ., Y ĉ(1) and Y ĉ(1) +1 , . . .Y T and find the estimator in each group.The whole procedure is repeated until a "constant mean is obtained".This procedure is called "binary segmentation".
The statistical computing environment R, package "changepoint", offers several methods for multiple change point detection.For the purpose of comparison, we will employ two basic approaches: detection of multiple change points using binary segmentation for normally distributed data and the method based on cumulative sums with the use of binary segmentation (see [1,[18][19][20]).

Simulation study
For the purpose of performance study of the proposed method of multiple change point detection we use simulations of process (2).We put T = 100, µ = 0, the error terms are independent normally distributed with zero mean and the standard deviations σ = 0.1, 0.2, 0.5 and 1. Locations of change points were chosen at random uniformly in the interval [5,95], jumps δ i (i = 1, . . ., p) taking heights ±1 (+ or − at random with equal probability).The minimum distance between neighboring change points is 5.For each number of change points out of p = 1, . . ., 15 we calculated 500 simulations of process (2).A uniform distribution was used in the simulations to generate random position of change points.In terms of change point detection, the change points are distributed within intervals of equal length with the same probability.Thus, there is no preferred position of change points.If we choose, for example, a normal distribution with the mean in the middle of the interval, position change points are more likely close to this center.As was shown in [6], detection of change points located in the middle of the process is more successful than in the case that this jump is near the beginning or end of the process.From this perspective, the choice of the uniform distribution is the most appropriate.
We applied ten methods to detect change points: • the direct BPDN estimator based on p most significant atoms; • the direct BPDN estimator based on p significant atoms exhibiting the highest correlation with the data; • the direct BPDN estimator based on the highest geometric mean of the correlation coefficient and the contribution to the total 1 norm; • the sequential BPDN estimator based on the most significant atom; • the sequential BPDN estimator based on the significant atom with the highest correlation; • the sequential BPDN estimator based on the highest geometric mean of correlation coefficient and the contribution to the total 1 norm; • the LASSO algorithm; • the adapted 1 trend filtering; • the estimator for normally distributed data using the binary segmentation method; • the cumulative sums estimator using the binary segmentation method.
We calculate empirical probability of successful change point detection as empirical probability = number of successful detections number of simulations .
The detection was a success provided that all p change points could be correctly identified.The results are summarized by plots in Fig. 2. The empirical probability of the sequential methods is generally higher than the probability of the direct detection methods.We assumed that the difference between two neighbouring atoms is at least 5.In the case of sequential methods we exclude all atoms detected in preceding iterations.
Basis pursuit denoising estimates were computed in Matlab utilizing [21] and fourstep modification of BPDN (BPA4) from the toolbox [22].A detailed description of the BPA4 algorithm can be found in [23].The LASSO estimates were calculated in R using the package "lars".The algorithm of least angle regression is used to get LASSO estimates [24].For the purpose of 1 trend filtering, the Matlab code written by Koh et al. [25] was employed.The smoothing parameter χ was for σ = 0.1 and 0.2 set to 1, for σ = 0.5 to 2 and for σ = 1 to 5. The package "changepoint" in the statistical environment R was used to obtain the last two estimators.

Real data example
In this part of the article we will apply all of the aforementioned methods of change point detection to time series comprising measurements of nuclear magnetic response made during the drilling of a well (see [26,27]).The changes should correspond to the transition between the different strata of rock.We will analyze a subset of this dataset containing 1200 measurements (from the index 1551 to 2750).Table 1 contains the results obtained by 10 described methods, where 8 change points were estimated.For each method there is listed a sequence of change point estimates in the order they were detected, then these estimates are sorted in ascending order for clarity.Seven possible position of change points were detected (in the range 134-135, 316-317, 496-499, 919-920, 981-982 and 1041-1042).Only the method based on the LASSO algorithm did not find the change in 981-982.The estimates of the eighth change point whose position is not apparent differs for each method.Similarly to the simulation study, we assumed that the difference between two neighboring change point estimates was at least 5.The smoothing parameter χ in adapted 1 trend filtering was set to 10 5 .

Conclusion
According to the simulation results the basis pursuit approach proposes a reasonable detection method of change points in one-dimensional process.In the case of one change point, the method based on significant atoms with the highest correlation with the process yields better results than the method based on pure 1 minimization when the standard deviation of the white noise is large.This is not valid for more than one change point which can be drawn from the graphs of the empirical probability of successful change point detection.It is apparent from the graphs that the sequential methods performed better than the direct ones.The LASSO algorithm gives similar results for one or two change points compared to non-sequential approaches when the level of the noise is small.For more than two change points the LASSO algorithm offers lower empirical probability of successful change point detection.For larger standard deviation of the noise and small number of change points the LASSO gives a comparable result as the method based on the correlation coefficient.As we can deduce from the simulations, the adapted 1 filtering gives almost comparable results with the sequential methods (except the method based on the highest correlation).However, for the standard deviation 0.2 and high number of change points the empirical probability of successful detection of the two proposed methods is higher than for adapted 1 filtering.If one wants to compare the aforementioned statistical methods, it can be said according to the calculated empirical probability, the estimators using the cumulative sums perform worse than the other one.The sequential methods give, in most cases, better results than the statistical approaches.
In the case of large standard deviation and higher number of change points the differences among the proposed methods nearly vanish.
In general, with the increasing number of change points the empirical probability decreases as one can expect.The standard deviation of an additive white noise of the process also affects the probability of successful detection.Higher values of the standard deviation cause lower probability of successful detection.
As to the computational speed it should be noted that the LASSO algorithm, the adapted 1 trend filtering and the algorithms for statistical estimates are significantly faster than the algorithm BPA4 used for the basis pursuit denoising.In addition to the problem of change point detection, the described four-step procedure using basis pursuit algorithm can solve a number of problems of various types, such as signal processing [9], kernel smoothing [23], and ARMA models fitting [28].This 4-step version of basis pursuit algorithm usually gives better results than pure basis pursuit approach.The price for this refinement is the increase in computation time.The aim of the article was not to monitor the speed of algorithms but to focus on the ability of each method to correctly detect potential changes.The dictionary created from Heaviside functions (atoms) can be merged with another set of functions, i.e.Fourier dictionary (see [9]) or Gaussian functions.The resulting dictionary can be applied to smooth the signal (the process) containing abrupt changes (jumps).
The change point detection techniques may be useful in a broad range of real-world problems, for instance in the modeling of economical or environmental time series where jumps can occur.
In general, the scope of sparse modeling based on basis pursuit or other related algorithms [29] goes far beyond the change-point problem.
It is because two factors are crucial for successful modeling: choice of an adequate model and a numerically stable estimation procedure yielding reliable parameter estimates.Unfortunately the ideas about a choice of an appropriate model are often very vague and produce models where it is hard to balance the requirement on sufficient regularity of the model (as few parameters as possible to guarantee numerical stability) and feasible precision which forces the analyst to increase the number of model components typically leading to overparameterization accompanied with non-uniqueness and numerical instability of solutions.Then the obtained numerical solution (parameter estimates) may strongly depend on the algorithm applied and consequently many of them may accidentally produce instable solutions corrupted by round-off error propagation.
Thus there is a danger that statistical analysis on parameters is partially or fully inadequate because the results of numerical computations may be far from the theoretically exact ones (the random effects are buried in round-off errors and cannot be separated).
For example, the papers [28,30] may illustrate the above reasoning quite well.

Fig. 3 .
Fig. 3. Nuclear magnetic response during the drilling of a well.

Table 1 .
The results of change point detection for time series of nuclear magnetic response during the drilling of a well.