Abstract
Regularizing the leastsquares criterion with the total number of coefficient changes, it is possible to estimate timevarying (TV) autoregressive (AR) models with piecewiseconstant coefficients. Such models emerge in various applications including speech segmentation, biomedical signal processing, and geophysics. To cope with the inherent lack of continuity and the high computational burden when dealing with highdimensional data sets, this article introduces a convex regularization approach enabling efficient and continuous estimation of TVAR models. To this end, the problem is cast as a sparse regression one with grouped variables, and is solved by resorting to the group leastabsolute shrinkage and selection operator (Lasso). The fresh look advocated here permeates benefits from advances in variable selection and compressive sampling to signal segmentation. An efficient blockcoordinate descent algorithm is developed to implement the novel segmentation method. Issues regarding regularization and uniqueness of the solution are also discussed. Finally, an alternative segmentation technique is introduced to improve the detection of change instants. Numerical tests using synthetic and real data corroborate the merits of the developed segmentation techniques in identifying piecewiseconstant TVAR models.
1. Introduction
Autoregressive (AR) models have been the workhorse for parametric spectral estimation since they form a dense set in the class of continuous spectra and, in many cases, they approximate parsimoniously the spectrum of a given random process [1], Chap. 3]. These are among the main reasons why AR models have been widely adopted in various applications as diverse as speech modeling [25], electroencephalogram (EEG) signal analysis [6], and geophysics [7]. While AR modeling of stationary random processes is well appreciated, a number of signals encountered in real life are nonstationary. This justifies the growing interest toward nonstationary signal analysis and timevarying (TV) AR models, which arise naturally in speech analysis due to the changing shape of the vocal tract as well as in EEG signal analysis due to the changes in the electrical activity of neurons. If the TVAR coefficient trajectories can be well approximated by superimposing a small number of basis sequences, nonstationary modeling reduces to estimating via, e.g., leastsquares (LS), the basis expansion coefficients [7]. On the other hand, it has been welldocumented that piecewiseconstant AR systems excited by white Gaussian noise are capable of modeling realworld signals such as speech and EEG [68]. Piecewiseconstant AR models constitute a subset of TVAR models wherein AR coefficients change abruptly. In this case, basis expansion techniques fall short in estimating the change points [8].
Exploiting the piecewise constancy of TVAR models, several methods are available to detect the changing instants of the AR coefficients, and thus facilitate what is often referred to as signal segmentation. The literature on signal segmentation is large since the topic is of interest in signal processing, applied statistics, and several other branches of science and engineering. Recent advances can be mainly divided in two categories. The first class adopts regularized LS criteria in order to impose piecewiseconstant AR coefficients. To avoid "oversegmentation," the LS cost is typically regularized with the total number of changes [6]. The resulting estimator can be implemented via dynamic programming (DP), which incurs computational burden that scales quadratically with the signal dimension. For large data sets, such as those considered in speech processing, this burden refrains practitioners from applying DP to segmentation, and heuristics are pursued instead based on the generalized likelihood ratio test (GLRT), or, approximations of the maximum likelihood approach [2], [[7], p. 401], [9]. The second class of methods relies on Bayesian inference and Markov Chain Monte Carlo (MCMC) methods [35]. A distinct advantage of this class is that model order selection can be performed automatically, and a variable model order can be chosen per segment. However, Bayesian techniques are known to require large computational resources.
The algorithm for change detection of piecewiseconstant AR models developed in this article belongs to the first class of methods, and its first novelty consists in developing a new regularization function which encourages piecewiseconstant TVAR coefficients while being convex and continuous; hence, it can afford efficient convex optimization solvers. To this end, it is shown that the segmentation problem can be recast as a sparse regression problem. The regularization function in [6] is then relaxed with its tightest convex approximation. It turns out that the resultant change detector is a modification of the group leastabsolute shrinkage and selection operator (Lasso) [10].
With the emphasis placed on large data sets, a candidate algorithm for implementing the developed change detector is a blockcoordinate descent iteration, which is provably convergent to the group Lasso solution. Surprisingly, it turns out that each iteration of the blockcoordinate descent can be implemented at complexity that scales linearly with the signal dimension, thus encouraging its application to large data sets. Regularization tuning and uniqueness of the group Lasso solution are also discussed.
The second novelty of the present study is an alternative changepoint retrieval algorithm based on the smoothlyclipped absolute deviation (SCAD) regularization. The associated nonconvex problem is tackled by resorting to a local linear approximation (LLA), which yields iterated weighted group Lasso minimization problems that can be solved via blockcoordinate descent. Numerical tests using synthetic and real (speech and sound) data are performed to corroborate the capability of the developed algorithms to identify piecewiseconstant TVAR models.
The remainder of the article is structured as follows. Section 2 deals with piecewiseconstant TVAR model estimation preliminaries. In Section 3, the problem at hand is recast as a sparse linear regression, and the novel group Lasso approach is introduced. An efficient blockcoordinate descent algorithm is developed in Section 4, while tuning issues and uniqueness of the group Lasso solution are addressed in Section 5. Section 6 introduces a nonconvex segmentation method based on the SCAD regularization to enhance the sparsity of the solution, which translates to retrieving more precisely the change instants. Numerical tests are presented in Section 7, and concluding remarks are summarized in Section 8. The Appendix is devoted to technical proofs. Notation: Column vectors (matrices) are denoted using lowercase (uppercase) boldface letters; calligraphic letters are reserved for sets; (·)^{T }stands for transposition, denotes the Gaussian probability density function with mean μ and variance σ^{2}; ⊗ denotes the Kronecker product; 0_{L }is the Ldimensional column vector with all zeros, and I_{L }is the Ldimensional identity matrix. The ℓ_{p }norm of is defined as .
2. Preliminaries and problem statement
Let denotes the realization of an Lth order TVAR process obeying the discretetime inputoutput relationship
where v_{n }denotes the zeromean white input noise at time n with variance , and a_{ℓ,n }is the ℓth TVAR coefficient at time n. With and , (1) can be rewritten as
Suppose that abrupt changes in the spectrum of {y_{n}} occur, due to piecewiseconstant changes of a_{n}; that is,
for k = 0,1,..., K, where K denotes the number of abrupt changes in the TVAR spectrum, and n_{k }the time instant of the kth abrupt change. The interval [n_{k}, n_{k+1 } 1] is referred to as the kth segment. Without loss of generality, n_{0 }= 0 and n_{K+1 } 1 = N.
The goal is to estimate the instants where the given time series {y_{n}} is split into K + 1 (stationary) segments, and also the constant AR coefficients per segment, i.e., . The number of abrupt changes, namely K, is not necessarily known.
2.1. Optimum segmentation of TVAR processes
Modeling real world signals using AR processes is well motivated because, for a given continuous spectral density S(f), it is possible to find an AR process (of high enough order) whose spectral density is arbitrarily close to S(f) [11], p. 130]. On the other hand, depending on the underlying nonstationary phenomena, variations of the AR coefficients can be either slow or abrupt. The problem stated in Section 2 is often referred to as signal segmentation, and emerges in numerous applications ranging from speech processing [25] to EEG signal analysis [6]. Regularized LS has been the workhorse approach for analyzing this kind of nonstationary processes [1215]. Denoting with μ a positive tuning constant, a Schwarzlike regularization is typically adopted to estimate jointly the change points and the AR coefficients, i.e.,
The nonconvex regularization term not only captures the total number of changes, but also encourages piecewiseconstant . Clearly, the larger the μ, the smaller the total number of changes. The estimator in (4) is optimal in the maximum a posteriori (MAP) sense when the change occurrences are modeled as Bernoulli random variables, and [6]. In some problems, the total number of changes is known, and the following constrained version of (4) is adopted instead:
From a practical point of view, the minimization in (4) or (6) is challenging since an exhaustive search over all possible sets of change instants has to be performed. However, several techniques based on DP, simulated annealing and interactive conditional model algorithms have been developed to evaluate (4) [6,16]. Despite the fact that DP approaches solve (4) in polynomial time, the computational complexity is quadratic in N, which limits their applicability to signal segmentation in practice. In typical applications, N can be very large (up to several thousands), and even quadratic complexity cannot be afforded. On the other hand, when applied to real data, the performance of the estimator in (4) is not satisfactory [17].
To overcome these limitations of (4), heuristic approaches based on the GLRT are used in real world applications [[7], p. 401], [9,18,19]. However, GLRTbased change detectors are sensitive to modeling errors, and require fine tuning of the associated detection thresholds.
In what follows, a convex relaxation of the cost in (4) is advocated based on recent advances in sparse linear regression and compressive sampling. To this end, (4) is first reformulated to a sparse regression problem with nonconvex regularization that is successively relaxed through its tightest convex approximation. The consequent optimization rule will yield sparse vector estimators which result in surprisingly accurate retrieval of changepoints. Those are obtained by an efficient blockcoordinate descent iteration that incurs only linear computational burden and memory storage. Unlike (4), based on wellestablished results in statistics, it will be further argued that the resultant TVAR model estimates are a continuous function of the data.
3. Sparse linear regression and group lassoing
Let denotes the observation vector, , for n = 0,1,..., N, and , such that
Define the "difference" vector as
and . Observe that d_{n }= 0_{L }for n > 0 if and only if there is no change in the TVAR coefficients between time instants n  1 and n. Clearly, it is possible to recover from since
Let denotes a lower triangular matrix with all nonzero entries equal to one and , having the following structure:
Since a = (T ⊗ I_{L})d, an equivalent formulation of (4) in terms of can be given as [cf. (7)]
What makes the formulation in (11) attractive but also challenging is the nonconvex and discontinuous Schwarzlike regularization term. The latter "pushes" most of the vectors toward 0_{L}, while d_{0 }is not penalized. As a consequence, the vector is group sparse, and the nonzero group indexes correspond to the change instants of the TVAR coefficients. Recently, a convex model selector with grouped variables was put forth by [10], and successfully applied to biostatistics and compressive sampling [20]. It generalizes the (nongrouped) leastabsolute shrinkage and selection operator (Lasso) [21] to regression problems where the unknown vector exhibits sparsity in groups; hence, its name group Lasso. The crux of group Lasso is to relax the Schwarzlike regularization in (11) with its tightest convex approximation.
The group Lasso is advocated here for catching change points by estimating the difference vectors as
where λ is a positive tuning parameter. It is known that the group Lasso regularization encourages group sparsity; that is, for most n > 0 [10]. Again, the larger the λ, the sparser the .
The role of the regularization term in (12) is illustrated next through a simple example. Select L = 2 for simplicity, and let d:= [d_{1}, d_{2}]^{T}. Consider the family of penalties , where 0 < p ≤ 2. Figure 1 depicts the penalties for p = 2,1,0.5, and 0.1. Clearly, is convex for 1 ≤ p ≤ 2. On the other hand, is nonconvex for 0 < p < 1 but it comes close to (d) as p approaches 0. This demonstrates that ║d_{n}║_{2 }is the tightest convex approximation of (d_{n}). Furthermore, ║d_{n}║_{2 }is nondifferentiable at d_{n }= 0_{L}, which enables group Lasso to encourage group sparsity. Needless to say that convexity of the regularizing functions avoids the presence of local minima, and allows for solving the resulting optimization problem efficiently. To this end, an efficient blockcoordinate descent algorithm is developed next, with computational complexity per iteration that scales linearly with N. But first two remarks are in order.
Figure 1. Regularization family for L = 2 and (a)p = 2, (b)p = 1, (c)p = 0.5, and (d)p = 0.1.
Remark 1. Different from Schwarzlike regularization, Figure 1 illustrates that the group Lasso one grows unbounded. This makes the resultant estimator biased. Nevertheless, unlike the Schwarzlike regularization, the group Lasso one is continuous which renders the resulting estimator more stable when applied to real data; see also [10,22]. A continuous regularization function that reduces the bias of the group Lasso will be discussed in Section 6.
Remark 2. Convex relaxation for detecting changes in the mean of nonstationary processes was recently mentioned in [12], and analyzed in [17]. For the meanchange problem, the tightest convex approximation of the Schwarzlike regularized LS is provided by the Lasso, which can afford efficient solvers such as the leastangle regression (LARS) algorithm [23]. However, for the group Lasso cost proposed here to catch changes in TVAR models, an exact LARSlike solver is not available [10]; thus, the pursuit of efficient algorithms for solving (12) is well motivated. This is the theme of the ensuing section.
4. Blockcoordinate descent solver
The crux of blockcoordinate descent is to iterate minimizations of the function of interest over a group of variables, while keeping the rest fixed. Consider the objective function
and let denotes the provisional solution at iteration i1. The nth step of the ith blockcoordinate descent iteration entails minimization of J(d) only with respect to d_{n}, while retaining the provisional estimates at iteration i1, namely , and the newly updated blocks at iterations i, namely . Thus, blockcoordinate descent at the nth step of the ith iteration yields
for n = 0,1,..., N, and i > 0. Skipping constant terms, J(d) in (13) can be rewritten as
where R := X^{T}X, and r := X^{T}y. Upon defining and for n' ≥ n , it holds that
and
While for n = 0 (14) reduces to an LS problem, for n > 0, omitting again irrelevant terms, it can be rewritten as
with
The problem in (18) is a convex secondorder cone program (SOCP). Typically, L ≪ N and (18) can be solved with fast optimization solvers based on interior point methods [24], at worstcase complexity . Recently, it has been shown that the solution of (18) can be obtained as a function of the solution of the following scalar problem
whose solution is given by [25]
Finally, in (18) can be obtained from in (21) as
Notice that if , the solution of (18) is . Since it is expected that the solution of (12) is sparse, solving (18) is trivial most of the time. If can be obtained via interior point methods or by (numerically) solving the scalar equation in (21), which admits fast solvers via, e.g., NewtonRaphson iterations, as in [25].
Despite the fact that blockcoordinate descent is typically adopted for largesize sparse linear regression, what makes it particularly appealing for catching changepoints is the fact that the vector can be updated recursively in n due to the special structure of R in (16). Upon defining
it follows from (19) that
which shows that evaluating requires the vectors and Given from the (i  1)st iteration, and initializing and at n = 0 as and , it is possible to recursively evaluate and given , and from step n1 for n > 0 as
The blockcoordinate descent algorithm is summarized in Algorithm 1. Interestingly, matrix in (12) does not have to be stored since only and suffice to implement Algorithm 1. Thus, the memory storage and complexity to perform one blockcoordinate descent iteration grow linearly with N. This attribute renders the blockcoordinate descent appealing especially for largesize problems where DP approaches tend to be too expensive.
Regarding convergence, the ensuing assertion is a direct consequence of the results in [26].
Proposition 1. The iterates obtained by Algorithm 1 converge to the global minimum of (12); that is, .
Blockcoordinate descent will also be the basic building block for solving the nonconvex problem introduced in Section 6 to improve the retrieval of changepoints. But first, it is useful to consider two issues of the group Lasso change detector for TVAR models.
Initialize with for n = 1, ... , N
for i > 0 do
for n = 0,1,..., N do
if n = 0 then
else
else
Algorithm 1: Blockcoordinate descent algorithm
5. Regularization and uniqueness issues
Performance of model selection with grouped variables via group Lasso and related approaches has been analyzed in [10,20,27], while asymptotic analysis has been pursued in [28,29]. In this section, a couple of issues are investigated regarding the group Lasso cost function, and the uniqueness of its minimum.
5.1. Tuning the regularization parameter
Selection of λ is a critical issue since larger λ's promote sparser solutions, which translate to fewer changes in the TVAR spectrum. However, larger λ's increase the estimator bias as well. If the number of changes present are known a priori by other means, or, if a certain level of segmentation can be afforded, λ can be tuned accordingly by 'trial and error,' or by crossvalidation. But in general, analytic methods to automatically choose the "best" value of λ are not available. In essence, selecting the regularization parameters is more a matter of engineering art, rather than systematic science.
In this section, heuristic but useful guidelines are provided to choose λ based on rigorously established lower bounds of this parameter. Given in (10), define such that X = [X_{0}, X_{1},..., X_{N}]. To bound λ, we will rely on the following result; see Appendix 1 for the proof.
Proposition 2. If X_{0 }has full column rank, then with , if and only if
If λ exceeds a threshold, which is specified by the regression matrix and the observations, Proposition 2 asserts that and for n = 1,..., N. This, along with (9), implies that for n = 0, 1,..., N; that is, no change occurs in the coefficients of the TVAR process. To avoid this trivial (changefree) solution, the guideline provided by Proposition 2 is that λ must be chosen strictly less than λ*. Our extensive simulations suggest that setting λ equal to a small percentage of λ*, say 520%, results in satisfactory estimates.
5.2. Uniqueness of the sparse solution
Uniqueness of sparse linear regression with nongrouped variables has been studied in [3032]. Next, uniqueness issues in recovering sparse vectors with groupvariables are explored by exploiting the deterministic structure of the regression matrix in (12). The cost function in (12) is not strictly convex since X is a fat matrix, and the regularization term is not strictly convex; see also Figure 1. On the other hand, the blockcoordinate descent algorithm developed in Section 4 is guaranteed to converge to a global minimum. In the following, a criterion is introduced to check a posteriori whether the obtained solution is unique for a given groupsparsity level.
Traditionally, the support of a sparse vector is defined as the set of indexes corresponding to the nonzero entries. In the groupsparsity framework herein, a different definition of support is required. Indeed, the vector of interest here, namely comprises N + 1 groups of Ldimensional variables. Since the term d0 is not penalized in (12), almost surely. Define the group support (gsupport) of to be the set containing the indexes relative to the nonzero group of variables; that is, In the following, when denotes the gsupport of d, the set is assumed ordered; i.e., s_{j }< s_{k }for each j < k.
The following lemma establishes a property of the matrix X in (10); see Appendix 2 for the proof.
Lemma 1. If any L out of N + 1 vectors are linearly independent, for any gsupport such that , and s_{j } s_{k} ≥ L for each j ≠ k, the matrix has full column rank.
Lemma 1 asserts full column rank of the submatrix , if it is formed by the columns of X corresponding to the nonzero indexes of any sparse vector whose gsupport is sufficiently small, and the nonzero groups are sufficiently distant from each other.
Next, Lemma 1 is exploited to establish an interesting property for the solutions of (12); see Appendix 3 for the proof.
Proposition 3. If any L out of N+1 vectors are linearly independent, for any gsupport such that , and s_{j } s_{k} ≥ L for each j ≠ k, there exists at most one solution of (12) gsupported in .
Proposition 3 ensures that if is gsupported in , and is sufficiently sparse with nonzero groups sufficiently far apart, then is the only solution of (12) gsupported in .
Remark 3. Analysis of the group Lasso and its modifications has revealed that its performance can be close to the Schwarzregularized LS either when the regression matrix is sufficiently block incoherent, or, when the block restricted isometry property holds [10,27]. If the regression matrix can be chosen by the designer, and it is randomly drawn from selected distributions (e.g., Gaussian or Bernoulli), these analyzes provide useful connections between problems (11) and (12). In the problem at hand however, the regression matrix is fixed, and its blocks [X_{0}, X_{1},..., X_{N}] are highly correlated. In this case, the relationship between the solutions of (11) and (12) is much less understood, and constitutes an interesting future research direction.
6. Continuity, bias and the group SCAD
As already pointed out, convex relaxation of the Schwarzlike regularization was developed in [12,17] for the meanchange problem using the (nongrouped) Lasso. Numerical results in [17] reveal that the Lasso tends to detect a "cloud" of small change points around an actual change. Postprocessing via DP to select a few of the estimated change instants was proposed in [17]. Moreover, due to the bias introduced by the Lasso, once the change points are obtained, another step is required to reestimate the mean within a segment. In the following, a novel change detector is developed based on recent advances in model selection via nonconvex regularization. The resulting estimator reduces the bias of group Lasso and can afford a convergent optimization solver. The corresponding algorithm is based on iterative instantiations of weighted group Lasso, which is capable of enhancing the sparsity of the solution [33,34], and thus improving the precision of the detected change points.
Attributes of a "good" regularization function are delineated in [22], and three properties are identified to this end:
• Unbiasedness. The estimator has to be unbiased when the true unknown parameter has large amplitude.
• Sparsity. The estimator has to set smallamplitude coefficients to zero to reduce the number of variables.
• Continuity. The estimator has to be continuous in the data to avoid instability when estimating (non)zero variables.
To appreciate these properties, a simple setting is presented next. Consider estimating a scalar parameter, call it d, in additive noise v, based on the scalar observation y = d + υ. In its general form, the regularized LS approach yields
where p_{λ}(d) is the regularization function. The Lasso regularization is while the Schwarzlike regularization is
In both cases, the estimate can be found in closed form, and the dependence of on y is given as
The nonlinear estimation rules in (30) and (31) are depicted in Figure 2a,b, respectively, for λ = 2. Observe that both regularization functions effect sparsity, since coefficients with small amplitude are set to 0. The Schwarzlike regularization yields unbiased estimates, but the solution is not continuous with respect to y. Hence, small variations of y or λ may result in large variations of (this happens when one is uncertain whether to set the coefficient to 0 or not). On the other hand, the Lasso regularization possesses continuity but the estimates are biased, because in addition to small, largeamplitude coefficients are "shrunk" too.
Figure 2. Regularized LS estimation rules in the scalar case: (a) (group) Lasso, (b) Schwarztype regularization, and (c) SCAD.
To overcome the limitations of these regularization functions, the following smoothly clipped absolute deviation (SCAD) regularization can be used with a > 2 [22]
to obtain estimates given by
The data dependence of is depicted in Figure 2c for λ = 2 and a = 3.7. Observe that the SCAD enjoys the three aforementioned attributes of a desirable regularization function.
Motivated by this scalar example, we propose as an alternative to (12), the following group SCAD approach for catching changepoints:
The problem in (34) is nonconvex and its exact minimization is challenging due to the presence of local minima. Nevertheless, it is possible to generalize the iterated local linear approximation (LLA) of [34] in order to ensure converge to a stationary point of the cost in (34). Let denote the derivative of for d≥0; that is,
The idea behind the LLA is to approximate with its linear expansion around d_{o}; that is
Given a provisional estimate of at iteration j  1, namely , the iterated LLA of (34) is
for j = 1,..., J. Since are nonnegative constants, the cost in (37) is convex, and can be minimized using the blockcoordinate descent algorithm of Section 4. The role of the weights is to avoid penalizing terms that, most likely, are nonzero. Function is depicted in Figure 3 for λ = 2, a = 3.7, and d ≥ 0. It is clearly a decreasing function of its argument. More precisely, if 0 ≤ d ≤ λ, then hence, the regularization parameter is as large as the group Lasso. If λ ≤ d ≤ aλ, the regularization parameter is linearly decreasing until d ≥ aλ, where the regularization parameter is zero. Thus, the expression in (37) represents an iterated weighted group Lasso. Furthermore, if initialized with = 0_{(N+1)L}, the first iteration of (37) corresponds to the (unweighted) group Lasso, and few iterations of the type in (37) are required for convergence.
Figure 3. Derivative of the SCAD function.
Remark 4. In principle, one could also apply a blockcoordinate descent iteration to minimize (34) directly. The resulting iterates converge to a local minimum of (34) that depends on the starting points; see also [26]. In general, it is impossible to assess properties of this solution. Instead, the solution of the LLA in (37) has provable merits in estimating the true support of sparse signals [34].
Remark 5. Recently, greedy algorithms such as the matching pursuit and the orthogonal matching pursuit have been shown to approach the performance of (11) and (12) when the regression matrix is sufficiently block incoherent, or, when the block restricted isometry property holds [27]. In the problem at hand, wherein the regression matrix exhibits correlation among blocks, simulated tests have shown that greedy algorithms suffer from severe error propagation. In fact, a cloud of change points is typically declared around a true change point. For these reasons, greedy algorithms will not be considered hereafter.
Remark 6. As already mentioned in the Introduction, one advantage of Bayesian techniques is that automatic model selection can be performed on a persegment basis at the price of increasing computational cost. The convex relaxation developed in Section 3 can be modified to perform model selection too. Indeed, setting L to a prescribed upper bound on the model order, one may further regularize the cost in (12) to impose sparse AR coefficient vectors; that is,
The cost function in (38) is convex and effects a piecewiseconstant and sparse TVAR model. It is different from model order selection criteria in the sense that the selected nonzero AR coefficients do not necessarily have to be consecutive. A challenge associated with the optimization in (38) is that blockcoordinate descent algorithms do not converge, since the differentiable part is not separable groupwise [26]. The cost function in (38) resembles the fused Lasso developed in [35]. Efficient algorithms exploiting this link, and the structure of the problem at hand are currently under investigation.
7. Simulated tests
The merits of the novel approaches to catching changepoints in TVAR processes are assessed via numerical simulations using synthetic and real data.
7.1. Synthetic data
The signal of interest here is a realization of a TVAR process with N + 1 = 500, order L = 4, , and σ^{2 }= 10^{2}, exhibiting K = 2 abrupt changes in the spectrum at time n_{1 }= 100, and n_{2 }= 350. The AR model during the first segment (n ∈ [0,99]) has coefficients a_{0 }= [0.8000, 0.1500, 0.1940, 0.0280]^{T}, while in the second segment (n ∈ [100,349]) the AR coefficients are a_{1 }= [0.1200, 0.0245, 0.2787, 0.0693]^{T}. In the third segment (n ∈ [350, 499]), a_{0 }is in act. Figure 4 shows the true signal along with the group Lasso, group SCAD, and DPbased estimates of the TVAR coefficients obtained via (12), (37), and (6) with K = 2, respectively. The regularization parameter λ of the group Lasso and group SCAD was selected to return 2 change points. The blockcoordinate descent algorithm of Section 4 was used to solve (12) up to a maximum of 10^{3 }iterations or when The same stopping rule is used for each iteration in (37), and J = 5 outer iterations are run (the same stopping rules are adopted henceforth). Figure 5 depicts the variation of across time. Clearly, the instants where represents the change points retrieved. While the DP has detected changepoints at and the group Lasso and the group SCAD have detected changes at and confirming that the results returned by the DP programming, by the group Lasso, and by group SCAD are comparable. Since {y_{n}} was generated as in (1), any L out of N + 1 vectors are linearly independent almost surely, which means that Proposition 3 can be invoked to ensure uniqueness of the group Lasso TVAR model estimate. Notice from Figure 4 and Figure 5 that the estimated by the group Lasso are much smaller than the true ones due to the bias, which results in poor estimates of the AR coefficients. However, those estimated by the group SCAD are closer to the true one.
Figure 4. Synthetic data. From top to bottom: True signal, estimated of the TVAR coefficients by group Lasso, group SCAD, and DP, respectively.
Figure 5. Synthetic data. From top to bottom: True TVAR coefficients change, changes detected by group Lasso, group SCAD, and DP, respectively.
7.2. Real data: piano sound
Next, a piano sound of 0.5 s comprising three monochromatic notes is sampled at 8 kHz to obtain N + L + 1 = 4000 samples. The signal is depicted in Figure 6. A TVAR model with L = 1 is adopted, and λ is selected to be λ*/10. Figure 7 depicts the TVAR coefficients estimated by the group Lasso and the group SCAD. The estimated over time is displayed in Figure 8. Observe that the group Lasso captures the correct changes at time instants and along with a small false change at time instant n_{f }= 2770. Instead, the group SCAD retrieves two changes at time instants and Notice also that group SCAD exhibits a reduced bias in the estimated which results in better estimates of the TVAR coefficients.
Figure 6. Piano sound comprising three monochromatic notes.
Figure 7. Piano sound. Estimated of the TVAR coefficients by group Lasso (top), and group SCAD (bottom).
Figure 8. Piano sound. Changes detected by group Lasso (top), and group SCAD (bottom).
7.3. Real data: speech
Here, a speech signal of 0.5 s is sampled at 8 kHz, to obtain N + L + 1 = 4000 samples. The resultant time series depicted in Figure 9 comprises a descent diphthong /ai/ followed by an /o/, pronounced by another party. A TVAR model with L = 8 is adopted, and λ is selected to be λ*/10. The change of vocoid in the diphthong occurs approximately at time instant n_{1 }= 1500, while the /o/ occurs approximately at n_{2 }= 3000. Figure 10 shows the TVAR coefficients estimated by the group Lasso and group SCAD. In agreement with [17], the group Lasso tends to declare a cloud of change points in the proximity of an actual change, while the jumps of the group SCAD estimates are very sharp. Figure 11 depicts ║d_{n}║_{2 }over time. The group SCAD detects four segments with change points at , and . Clearly, the first segment corresponds to the /a/, the second, which is the shortest, to the transition of the diphthong, the third to the /i/, and the forth to the /o/. Observe that the group Lasso exhibits peaks around the actual change instants, and requires postprocessing via either DP as advocated in [17], or, by simply peak picking.
Figure 9. Speech signal: /ai//o/.
Figure 10. Speech signal. Estimated of the TVAR coefficients by group Lasso (top), and group SCAD (bottom).
Figure 11. Speech signal. Changes detected by group Lasso (top), and group SCAD (bottom).
Further tests are performed on a sampled speech that has been widely adopted for TVAR speech segmentation [35,7]. According to [7, p. 401], this signal belongs to a database designed by the French National Agency for Telecommunications (CNET) for testing and evaluating speech recognition algorithms. It consists of a noisy French speech recorded in a car at sampling frequency 8 kHz, prefiltered by a highpass filter with cutoff frequency equal to 150 Hz, and quantized with 16 bits per sample. The time series is shown in Figure 12 along with the changes caught by the group SCAD change detector, and the GLR algorithm of [18] whose results are reported in [35] for L = 2. The group SCAD is tuned to return the same number of changes as the GLR. The detected change instants are listed in Table 1.
Figure 12. French speech signal and changes detected by the group SCAD (dot) and GLR (dash) change TVAR detectors forL = 2.
Table 1. Change instants detected by the group SCAD and GLR algorithm
The first change detected by the GLR is at sample 445, while this change is not detected by the group SCAD. Interestingly, it is reported in [3] that this change is not relevant for segmentation purposes, and this fact is apparent by inspection of the true signal. Both algorithms have detected changes around samples 1750, 2100, 2800, and 3650. The group SCAD successfully removed the false change detected by the GLR at sample 3400. By inspecting the original signal, this change does not seem to be relevant. The group SCAD had detected a change instant around sample 1300 unlike the GLR. This change is also detected by advanced Bayesian techniques reported in [35]. The group SCAD has detected a change at sample 779, while the GLR at sample 645. Indeed, inspection of the original signal suggests that the detection of the GLR is preferable. Surprisingly, the group SCAD, unlike the GLR and the Bayesian techniques of [35], has detected a change at sample 2359. Observing the original signal around this point, there is a clear amplitude modulation that may cause a change in the TVAR coefficients, which existing algorithms have passed undetected.
A way to univocally compare the two algorithms is via the segmented prediction error (SPE). Assuming that K changes have been detected at instants , let denotes the LS estimates of the AR model in the kth segment, i.e., . The SPE is defined as , and represents the error in approximating the original signal {y_{n}} with a TVAR model exhibiting abrupt changes at instants . The GLR segmentation entails SPE_{glr} = 0.2638, while SPE_{g}_{scad }= 0.2578 for the group SCAD. Clearly, the group SCAD based segmentation seems preferable to that of the GLR algorithm.
Finally, since general analysis of the convergence rate for the blockcoordinate descent algorithm is not available, a simulated test assessing its converge speed is performed. Figure 13 depicts the normalized step size amplitude variations, namely , across the iteration index (i) for the first LLA of the group SCAD (which amounts to a group Lasso since for the speech signal in Figure 12. Observe that after a few hundred iterations, the normalized step amplitude size drops below 10^{6}. For practical purposes, the solution at this stage might be acceptable. Progressively, the speed of convergence slows down since the algorithm has to decide whether some components are truly zero or have small amplitudes, and large jumps in correspond to iteration indexes where the small components are set to zero. Once the correct vector support is determined, the speed of converge is fast. It is worth pointing out that only 10 min are required to perform 5,000 iterations of the blockcoordinate descent (implemented with Matlab) for the problem at hand, whereas advanced Bayesian techniques might require several hours (see e.g., [3]).
Figure 13. Speed of the blockcoordinate descent for the first LLA of the group SCAD.
8. Concluding remarks
Novel estimators were developed in this article for identification of piecewiseconstant TVAR models by exploiting recent advances in variable selection and compressive sampling. While traditional techniques consist in regularizing a LS criterion with the total number of coefficient changes, the novel estimator relies on a convex regularization function, which resembles the group Lasso and can afford efficient implementation using blockcoordinate descent iterations. The latter incurs computational burden that scales linearly with the number of data samples, thus being particularly attractive for largesize problems. Regularization tuning issues are discussed along with conditions for uniqueness of the estimated piecewiseconstant AR model. An alternative group smoothlyclipped absolute deviation regularization is also introduced, and an algorithm based on iterative weighted group Lasso minimizations is developed. Numerical tests using synthetic and real data confirm that the developed algorithms can effectively identify piecewiseconstant AR models of large size at manageable complexity, and outperform heuristic alternatives that are based on the GLRT.
Competing interests
The authors declare that they have no competing interests.
This paper was originally submitted to the IEEE Transactions on Signal Processing on April 6, 2010, and rejected on June 8, 2010  date on which we became aware of the independent work by H. Ohlsson, L. Ljung, and S. Boyd, "Segmentation of ARXmodels using sumofnorms regularization," Automatica, June 2010, which overlaps with the present contribution in the model and criterion employed, but it is distinct in the solution, application, and real data experimentation.
Appendix 1: proof of proposition 2
The necessary and sufficient firstorder optimality condition for to be the (unconstrained) minimum of (12), is that the subgradient of J(d) in (13) evaluated at contains the zero vector [[36], p. 92]; i.e.,
Defining as w := X^{T}(Xd  y) and , with for n = 0,1,..., N, the subgradient of J(d) evaluated at is given by
with such that ║s_{n}║_{2 }≤ 1. Using (40) and (41), Equation (39) translates to the following conditions:
(c1) w_{0 }= 0_{l}; and,
The changefree solution corresponds to having for n = 1,..., N. Thus, (c1) implies that , which is uniquely satisfied by , since X_{0 }has full column rank. Hence, and for n = 1,..., N hold if and only if (c2) is satisfied, which corresponds to ║w_{n}║_{2 }≤ λ for n = 1,..., N. Since , condition (c2) is satisfied if and only if .
Appendix 2: proof of lemma 1
Observe that
Notice that the first subsum in (42) comprises rank1 matrices, the last subsum comprises s_{1 }rank1 matrices, while the gth subsum comprises rank1 matrices for . Since is such that , and s_{j } s_{k} ≥ L for each j ≠ k, and any L out of N + 1 vectors are linearly independent, each of the summands in (42) has rank L. Thus, it is possible to find L linearly independent vectors such that the first subsum in (42) equals to with . Analogously, it is possible to find L linearly independent vectors such that the gth subsum in (42) can be written as with for . Finally, it is possible to find L linearly independent vectors such that the last subsum in (42) can be written as with . Thus, , and since are ( + 1) linearly independent vectors, has fullcolumn rank.
Appendix 3: proof of proposition 3
Suppose that and are two solutions of (12) with the same gsupport . Let denotes the matrix obtained by selecting the columns of X relative to the nonzero indexes dictated by the set . The minimization of (13) over the vector having gsupport in amounts to
From Lemma 1, has full column rank which implies that is strictly convex, and so is . Thus, (43) admits a unique solution.
Since both and are gsupported in by hypothesis, and their restrictions to are equal to û, it follows readily that .
Acknowledgements
This work was supported by MURI (AFOSR FA95501010567) grant.
References

P Stoica, RL Moses, Introduction to Spectral Analysis (PrenticeHall, New Jersey, 1997)

PM Djuric, A MAP solution to offline segmentation of signals. Proc of the International Conference on Acoustics, Speech, and Signal Processing (Adelaide, Australia, 1994) 4, pp. 505–508

N Dobigeon, JY Tourneret, M Davy, Joint segmentation of piecewise constant autoregressive processes by using a hierarchical model and a Bayesian sampling approach. IEEE Trans Signal Process 55(4), 1251–1263 (2007)

P Fearnhead, Exact Bayesian curve fitting and signal segmentation. IEEE Trans Signal Process 53(6), 2160–2166 (2005)

E Punskaya, C Andrieu, A Doucet, WJ Fitzgerald, Bayesian curve fitting using MCMC with applications to signal segmentation. IEEE Trans Signal Process 50(3), 747–758 (2002)

M Lavielle, Optimal segmentation of random processes. IEEE Trans Signal Process 46(5), 1365–1373 (1998)

M Basseville, IV Nikiforov, Detection of Abrupt Changes: Theory and Application (PrenticeHall, Englewood Cliffs, NJ, USA, 1993)

MG Hall, AV Oppenheim, AS Willsky, Timevarying parametric modeling of speech. Signal Process 5(3), 267–285 (1983)

D Rudoy, TF Quatieri, PJ Wolfe, Timevarying autoregressive tests for multiscale speech analysis. Proceedings of Interspeech (Brighton, UK, 2009) 1, pp. 2839–2842

M Yuan, Y Lin, Model selection and estimation in regression with grouped variables. J Royal Stat Soc, Ser B 68(1), 49–67 (2006)

PJ Brockwell, RA Davis, Time Series: Theory and Methods (SpringerVerlag, New York, NY, USA, 1990)

L Boysen, A Kempe, V Liebscher, A Munk, O Wittich, Consistencies and rates of convergence of jumppenalized leastsquares estimators. Annals Stat 37(1), 157–183 (2009)

M Lavielle, Using penalized contrasts for the changepoint problem. Signal Process 85(8), 1501–1510 (2005)

M Lavielle, E Moulines, Leastsquares estimation of an unknown number of shifts in a time series. J Time Series Anal 21(1), 33–59 (2000)

E Lebarbier, Detecting multiple changepoints in the mean of Gaussian process by model selection. Signal Process 85(4), 717–736 (2005)

G Winkler, V Liebscher, Smoothers for discontinuous signals, J. Nonparametric Stat 14(12), 203–222 (2002)

Z Harchaoui, C LevyLeduc, Catching changepoints with Lasso. Proceedings of the Advanced Neural Information Processes Systems (Vancouver, Canada, 2008) 20, pp. 161–168

U Appel, AV Brandt, Adaptive sequential segmentation of piecewise stationary time series. Inf Sci 29(1), 27–56 (1983)

A Willsky, H Jones, A generalized likelihood ratio approach to the detection and estimation of jumps in linear systems. IEEE Trans Autom Control 21(1), 108–112 (1976)

YC Eldar, M Mishali, Block sparsity and sampling over a union of subspaces. Proc of the 16th International Conference on Digital Signal Processing (Santorini, Greece, 2009) 1, pp. 1–8

R Tibshirani, Regression shrinkage and selection via the Lasso. J Royal Stat Soc Ser B 58(1), 267–288 (1996)

J Fan, R Li, Variable selection via nonconcave penalized likelihood and its oracle properties. J Am Stat Assoc 96, 1348–1360 (2001)

B Efron, T Hastie, I Johnstone, R Tibshirani, Least angle regression. Annals Stat 32, 407–499 (2004)

JF Sturm, Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones. Optim Methods Softw 1112, 625–653 (1999)

A Puig, A Wiesel, A Hero, A multidimensional shrinkagethresholding operator. Proceedings of the 15th Workshop on Statistical Signal Processing (Cardiff, UK, 2009) 18, pp. 363–366

P Tseng, Convergence of a block coordinate descent method for nondifferentiable minimization. J Optim Theory Appl 109(3), 475–494 (2001)

YC Eldar, P Kuppinger, H Bölcskei, Blocksparse signals: Uncertainty relations and efficient recovery. IEEE Trans Signal Process 58, 3042–3054 (2010)

FR Bach, Consistency of the group Lasso and multiple kernel learning. J Mach Learn Res 9, 1179–1225 (2008)

Y Nardi, A Rinaldo, On the asymptotic properties of the group Lasso estimator for linear models. Electron J Stat 2, 605–633 (2008)

JJ Fuchs, On sparse representations in arbitrary redundant bases. IEEE Trans. Inf Theory 50(6), 1341–1344 (2004)

I Gorodnitsky, B Rao, Sparse signal reconstruction from limited data using FOCUSS: a reweighted minimum norm algorithm. IEEE Trans Signal Process 45(3), 600–616 (1997)

J Tropp, Just relax: convex programming methods for identifying sparse signals in noise. IEEE Trans Inf Theory 52(3), 1030–1051 (2006)

EJ Candes, MB Wakin, S Boyd, Enhancing sparsity by reweighted ℓ_{1 }minimization. J Fourier Anal Appl 14(5), 877–905 (2008)

H Zou, R Li, Onestep sparse estimates in nonconcave penalized likelihood models. Annals Stat 36, 1509–1533 (2008)

R Tibshirani, M Saunders, S Rosset, J Zhu, k Knight, Sparsity and smoothness via the fused Lasso. J Royal Stat Soc Ser B 67(1), 91–108 (2005)

A Ruszczynski, Nonlinear Optimization (Princeton University Press, Princeton, NJ, USA, 2006)