SpringerOpen Newsletter

Receive periodic news and updates relating to SpringerOpen.

This article is part of the series Sparse Signal Processing.

Open Access Research

Group lassoing change-points in piecewise-constant AR processes

Daniele Angelosante1* and Georgios B Giannakis2

Author Affiliations

1 Asea Brown Boveri (ABB) Corporate Research Center, Baden, CH 5405, Switzerland

2 Department of Electrical and Computer Engineering, University of Minnesota, Minneapolis, MN 55455, USA

For all author emails, please log on.

EURASIP Journal on Advances in Signal Processing 2012, 2012:70  doi:10.1186/1687-6180-2012-70

The electronic version of this article is the complete one and can be found online at: http://asp.eurasipjournals.com/content/2012/1/70


Received:31 October 2011
Accepted:21 March 2012
Published:21 March 2012

© 2012 Angelosante and Giannakis; licensee Springer.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Regularizing the least-squares criterion with the total number of coefficient changes, it is possible to estimate time-varying (TV) autoregressive (AR) models with piecewise-constant 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 high-dimensional data sets, this article introduces a convex regularization approach enabling efficient and continuous estimation of TV-AR models. To this end, the problem is cast as a sparse regression one with grouped variables, and is solved by resorting to the group least-absolute 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 block-coordinate 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 piecewise-constant TV-AR 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 [2-5], 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 non-stationary. This justifies the growing interest toward non-stationary signal analysis and time-varying (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 TV-AR coefficient trajectories can be well approximated by superimposing a small number of basis sequences, non-stationary modeling reduces to estimating via, e.g., least-squares (LS), the basis expansion coefficients [7]. On the other hand, it has been well-documented that piecewise-constant AR systems excited by white Gaussian noise are capable of modeling real-world signals such as speech and EEG [6-8]. Piecewise-constant AR models constitute a subset of TV-AR 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 TV-AR 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 piecewise-constant 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 [3-5]. 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 piecewise-constant 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 piecewise-constant TV-AR 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 least-absolute 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 block-coordinate descent iteration, which is provably convergent to the group Lasso solution. Surprisingly, it turns out that each iteration of the block-coordinate 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 change-point retrieval algorithm based on the smoothly-clipped absolute deviation (SCAD) regularization. The associated non-convex problem is tackled by resorting to a local linear approximation (LLA), which yields iterated weighted group Lasso minimization problems that can be solved via block-coordinate descent. Numerical tests using synthetic and real (speech and sound) data are performed to corroborate the capability of the developed algorithms to identify piecewise-constant TV-AR models.

The remainder of the article is structured as follows. Section 2 deals with piecewise-constant TV-AR 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 block-coordinate 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 non-convex 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 lower-case (upper-case) boldface letters; calligraphic letters are reserved for sets; (·)T stands for transposition, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M162','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M162">View MathML</a> denotes the Gaussian probability density function with mean μ and variance σ2; ⊗ denotes the Kronecker product; 0L is the L-dimensional column vector with all zeros, and IL is the L-dimensional identity matrix. The p norm of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M163','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M163">View MathML</a> is defined as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M1">View MathML</a>.

2. Preliminaries and problem statement

Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M2">View MathML</a> denotes the realization of an Lth order TV-AR process obeying the discrete-time input-output relationship

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M3">View MathML</a>

(1)

where vn denotes the zero-mean white input noise at time n with variance <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M4">View MathML</a>, and aℓ,n is the th TV-AR coefficient at time n. With <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M164','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M164">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M165','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M165">View MathML</a>, (1) can be rewritten as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M5">View MathML</a>

(2)

Suppose that abrupt changes in the spectrum of {yn} occur, due to piecewise-constant changes of an; that is,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M6">View MathML</a>

(3)

for k = 0,1,..., K, where K denotes the number of abrupt changes in the TV-AR spectrum, and nk the time instant of the kth abrupt change. The interval [nk, nk+1 - 1] is referred to as the kth segment. Without loss of generality, n0 = 0 and nK+1 - 1 = N.

The goal is to estimate the instants <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M7">View MathML</a> where the given time series {yn} is split into K + 1 (stationary) segments, and also the constant AR coefficients per segment, i.e., <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M8">View MathML</a>. The number of abrupt changes, namely K, is not necessarily known.

2.1. Optimum segmentation of TV-AR 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 non-stationary 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 [2-5] to EEG signal analysis [6]. Regularized LS has been the workhorse approach for analyzing this kind of non-stationary processes [12-15]. Denoting with μ a positive tuning constant, a Schwarz-like regularization is typically adopted to estimate jointly the change points and the AR coefficients, i.e.,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M9">View MathML</a>

(4)

where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M10">View MathML</a> is defined as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M11">View MathML</a>

(5)

The non-convex regularization term <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M12">View MathML</a> not only captures the total number of changes, but also encourages piecewise-constant <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M13">View MathML</a>. 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M166','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M166">View MathML</a>[6]. In some problems, the total number of changes is known, and the following constrained version of (4) is adopted instead:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M14">View MathML</a>

(6)

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, GLRT-based 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 non-convex 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 change-points. Those are obtained by an efficient block-coordinate descent iteration that incurs only linear computational burden and memory storage. Unlike (4), based on well-established results in statistics, it will be further argued that the resultant TV-AR model estimates are a continuous function of the data.

3. Sparse linear regression and group lassoing

Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M167','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M167">View MathML</a> denotes the observation vector, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M15">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M16">View MathML</a> for n = 0,1,..., N, and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M168','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M168">View MathML</a>, such that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M17">View MathML</a>

(7)

Define the "difference" vector <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M169','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M169">View MathML</a> as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M18">View MathML</a>

(8)

and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M19">View MathML</a>. Observe that dn = 0L for n > 0 if and only if there is no change in the TV-AR coefficients between time instants n - 1 and n. Clearly, it is possible to recover <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M20">View MathML</a> from <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M21">View MathML</a> since

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M22">View MathML</a>

(9)

Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M170','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M170">View MathML</a> denotes a lower triangular matrix with all nonzero entries equal to one and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M171','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M171">View MathML</a>, having the following structure:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M23">View MathML</a>

(10)

Since a = (T IL)d, an equivalent formulation of (4) in terms of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M21">View MathML</a> can be given as [cf. (7)]

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M24">View MathML</a>

(11)

What makes the formulation in (11) attractive but also challenging is the non-convex and discontinuous Schwarz-like regularization term. The latter "pushes" most of the <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M25">View MathML</a> vectors toward 0L, while d0 is not penalized. As a consequence, the vector <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M26">View MathML</a> is group sparse, and the non-zero group indexes correspond to the change instants of the TV-AR 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 (non-grouped) least-absolute 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 Schwarz-like regularization in (11) with its tightest convex approximation.

The group Lasso is advocated here for catching change points by estimating the difference vectors as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M27">View MathML</a>

(12)

where λ is a positive tuning parameter. It is known that the group Lasso regularization encourages group sparsity; that is, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M28">View MathML</a> for most n > 0 [10]. Again, the larger the λ, the sparser the <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M29">View MathML</a>.

The role of the regularization term in (12) is illustrated next through a simple example. Select L = 2 for simplicity, and let d:= [d1, d2]T. Consider the family of penalties <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M30">View MathML</a>, where 0 < p ≤ 2. Figure 1 depicts the penalties <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M31">View MathML</a> for p = 2,1,0.5, and 0.1. Clearly, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M31">View MathML</a> is convex for 1 ≤ p ≤ 2. On the other hand, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M31">View MathML</a> is non-convex for 0 < p < 1 but it comes close to <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M32">View MathML</a>(d) as p approaches 0. This demonstrates that ║dn2 is the tightest convex approximation of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M33">View MathML</a> (dn). Furthermore, ║dn2 is non-differentiable at dn = 0L, 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 block-coordinate descent algorithm is developed next, with computational complexity per iteration that scales linearly with N. But first two remarks are in order.

thumbnailFigure 1. Regularization family <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M161','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M161">View MathML</a> for L = 2 and (a)p = 2, (b)p = 1, (c)p = 0.5, and (d)p = 0.1.

Remark 1. Different from Schwarz-like regularization, Figure 1 illustrates that the group Lasso one grows unbounded. This makes the resultant estimator biased. Nevertheless, unlike the Schwarz-like 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 non-stationary processes was recently mentioned in [12], and analyzed in [17]. For the mean-change problem, the tightest convex approximation of the Schwarz-like regularized LS is provided by the Lasso, which can afford efficient solvers such as the least-angle regression (LARS) algorithm [23]. However, for the group Lasso cost proposed here to catch changes in TV-AR models, an exact LARS-like 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. Block-coordinate descent solver

The crux of block-coordinate descent is to iterate minimizations of the function of interest over a group of variables, while keeping the rest fixed. Consider the objective function

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M34">View MathML</a>

(13)

and let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M35">View MathML</a> denotes the provisional solution at iteration i-1. The nth step of the ith block-coordinate descent iteration entails minimization of J(d) only with respect to dn, while retaining the provisional estimates at iteration i-1, namely <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M36">View MathML</a>, and the newly updated blocks at iterations i, namely <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M37">View MathML</a>. Thus, block-coordinate descent at the nth step of the ith iteration yields

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M38">View MathML</a>

(14)

for n = 0,1,..., N, and i > 0. Skipping constant terms, J(d) in (13) can be rewritten as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M39">View MathML</a>

(15)

where R := XTX, and r := XTy. Upon defining <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M40">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M41">View MathML</a> for n' n , it holds that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M42">View MathML</a>

(16)

and

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M43">View MathML</a>

(17)

While for n = 0 (14) reduces to an LS problem, for n > 0, omitting again irrelevant terms, it can be rewritten as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M44">View MathML</a>

(18)

with

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M45">View MathML</a>

(19)

The problem in (18) is a convex second-order cone program (SOCP). Typically, L N and (18) can be solved with fast optimization solvers based on interior point methods [24], at worst-case complexity <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M172','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M172">View MathML</a>. Recently, it has been shown that the solution of (18) can be obtained as a function of the solution of the following scalar problem

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M46">View MathML</a>

(20)

whose solution is given by [25]

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M47">View MathML</a>

(21)

Finally, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M48">View MathML</a> in (18) can be obtained from <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M49">View MathML</a> in (21) as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M50">View MathML</a>

(22)

Notice that if <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M51">View MathML</a>, the solution of (18) is <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M52">View MathML</a>. Since it is expected that the solution of (12) is sparse, solving (18) is trivial most of the time. If <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M53">View MathML</a> can be obtained via interior point methods or by (numerically) solving the scalar equation in (21), which admits fast solvers via, e.g., Newton-Raphson iterations, as in [25].

Despite the fact that block-coordinate descent is typically adopted for large-size sparse linear regression, what makes it particularly appealing for catching change-points is the fact that the vector <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M54">View MathML</a> can be updated recursively in n due to the special structure of R in (16). Upon defining

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M55">View MathML</a>

(23)

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M56">View MathML</a>

(24)

it follows from (19) that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M57">View MathML</a>

(25)

which shows that evaluating <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M54">View MathML</a> requires the vectors <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M58">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M59">View MathML</a> Given <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M60">View MathML</a> from the (i - 1)st iteration, and initializing <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M58">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M59">View MathML</a> at n = 0 as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M61">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M62">View MathML</a>, it is possible to recursively evaluate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M58">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M59">View MathML</a> given <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M63">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M64">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M65">View MathML</a> from step n-1 for n > 0 as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M66">View MathML</a>

(26)

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M67">View MathML</a>

(27)

The block-coordinate descent algorithm is summarized in Algorithm 1. Interestingly, matrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M173','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M173">View MathML</a> in (12) does not have to be stored since only <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M68">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M69">View MathML</a> suffice to implement Algorithm 1. Thus, the memory storage and complexity to perform one block-coordinate descent iteration grow linearly with N. This attribute renders the block-coordinate descent appealing especially for large-size 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M70">View MathML</a> obtained by Algorithm 1 converge to the global minimum of (12); that is, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M71">View MathML</a>.

Block-coordinate descent will also be the basic building block for solving the non-convex problem introduced in Section 6 to improve the retrieval of change-points. But first, it is useful to consider two issues of the group Lasso change detector for TV-AR models.

Given <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M72">View MathML</a>

Initialize with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M73','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M73">View MathML</a> for n = 1, ... , N

for i > 0 do

  for n = 0,1,..., N do

    if n = 0 then

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M61">View MathML</a>

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M74','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M74">View MathML</a>

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M75">View MathML</a>

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M76">View MathML</a>

    else

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M77">View MathML</a>

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M78">View MathML</a>

      <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M79">View MathML</a>

      if <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M80">View MathML</a>then

        <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M81">View MathML</a>

      else

        <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M82">View MathML</a>

         Algorithm 1: Block-coordinate 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 TV-AR 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 cross-validation. 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M174','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M174">View MathML</a> in (10), define <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M175','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M175">View MathML</a> such that X = [X0, X1,..., XN]. To bound λ, we will rely on the following result; see Appendix 1 for the proof.

Proposition 2. If X0 has full column rank, then <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M83">View MathML</a>with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M84">View MathML</a>, if and only if <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M85">View MathML</a>

If λ exceeds a threshold, which is specified by the regression matrix and the observations, Proposition 2 asserts that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M86">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M87">View MathML</a> for n = 1,..., N. This, along with (9), implies that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M88">View MathML</a> for n = 0, 1,..., N; that is, no change occurs in the coefficients of the TV-AR process. To avoid this trivial (change-free) 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 5-20%, results in satisfactory estimates.

5.2. Uniqueness of the sparse solution

Uniqueness of sparse linear regression with non-grouped variables has been studied in [30-32]. Next, uniqueness issues in recovering sparse vectors with group-variables 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 block-coordinate 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 group-sparsity level.

Traditionally, the support of a sparse vector is defined as the set of indexes corresponding to the non-zero entries. In the group-sparsity framework herein, a different definition of support is required. Indeed, the vector of interest here, namely <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M89">View MathML</a> comprises N + 1 groups of L-dimensional variables. Since the term d0 is not penalized in (12), <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M90">View MathML</a> almost surely. Define the group support (g-support) of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> to be the set containing the indexes relative to the nonzero group of variables; that is, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M92">View MathML</a> In the following, when <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M176','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M176">View MathML</a> denotes the g-support of d, the set is assumed ordered; i.e., sj < sk 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M93">View MathML</a>are linearly independent, for any g-support <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M188','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M188">View MathML</a>such that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M180','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M180">View MathML</a>, and |sj - sk| ≥ L for each j ≠ k, the matrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M94">View MathML</a>has full column rank.

Lemma 1 asserts full column rank of the submatrix <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178">View MathML</a>, if it is formed by the columns of X corresponding to the non-zero indexes of any sparse vector whose g-support is sufficiently small, and the non-zero 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M93">View MathML</a>are linearly independent, for any g-support <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M179','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M179">View MathML</a>such that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M180','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M180">View MathML</a>, and |sj - sk| L for each j ≠ k, there exists at most one solution of (12) g-supported in .

Proposition 3 ensures that if <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> is g-supported in , and is sufficiently sparse with non-zero groups sufficiently far apart, then <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> is the only solution of (12) g-supported in .

Remark 3. Analysis of the group Lasso and its modifications has revealed that its performance can be close to the Schwarz-regularized 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 [X0, X1,..., XN] 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 Schwarz-like regularization was developed in [12,17] for the mean-change problem using the (non-grouped) Lasso. Numerical results in [17] reveal that the Lasso tends to detect a "cloud" of small change points around an actual change. Post-processing 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 re-estimate the mean within a segment. In the following, a novel change detector is developed based on recent advances in model selection via non-convex 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 small-amplitude 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

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M95">View MathML</a>

(28)

where pλ(|d|) is the regularization function. The Lasso regularization is <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M96">View MathML</a> while the Schwarz-like regularization is

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M97">View MathML</a>

(29)

In both cases, the estimate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M190','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M190">View MathML</a> can be found in closed form, and the dependence of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M190','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M190">View MathML</a> on y is given as

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M98">View MathML</a>

(30)

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M99">View MathML</a>

(31)

The non-linear 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 Schwarz-like 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M190','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M190">View MathML</a> (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, large-amplitude coefficients are "shrunk" too.

thumbnailFigure 2. Regularized LS estimation rules in the scalar case: (a) (group) Lasso, (b) Schwarz-type 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]

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M100">View MathML</a>

(32)

to obtain estimates given by

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M101">View MathML</a>

(33)

The data dependence of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M102">View MathML</a> 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 change-points:

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M103">View MathML</a>

(34)

The problem in (34) is non-convex 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M104">View MathML</a> denote the derivative of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M105">View MathML</a>for d≥0; that is,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M106">View MathML</a>

(35)

The idea behind the LLA is to approximate <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M107">View MathML</a> with its linear expansion around do; that is

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M108">View MathML</a>

(36)

Given a provisional estimate of <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M109">View MathML</a> at iteration j - 1, namely <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M110">View MathML</a>, the iterated LLA of (34) is

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M111">View MathML</a>

(37)

for j = 1,..., J. Since <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M112">View MathML</a> are non-negative constants, the cost in (37) is convex, and can be minimized using the block-coordinate descent algorithm of Section 4. The role of the weights <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M112">View MathML</a> is to avoid penalizing terms that, most likely, are non-zero. Function <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M107">View MathML</a> 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M113">View MathML</a> 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M114">View MathML</a> = 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.

thumbnailFigure 3. Derivative of the SCAD function.

Remark 4. In principle, one could also apply a block-coordinate 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 per-segment 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,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M115">View MathML</a>

(38)

The cost function in (38) is convex and effects a piecewise-constant and sparse TV-AR model. It is different from model order selection criteria in the sense that the selected non-zero AR coefficients do not necessarily have to be consecutive. A challenge associated with the optimization in (38) is that block-coordinate descent algorithms do not converge, since the differentiable part is not separable group-wise [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 change-points in TV-AR 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 TV-AR process with N + 1 = 500, order L = 4, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M166','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M166">View MathML</a>, and σ2 = 10-2, exhibiting K = 2 abrupt changes in the spectrum at time n1 = 100, and n2 = 350. The AR model during the first segment (n ∈ [0,99]) has coefficients a0 = [-0.8000, -0.1500, 0.1940, -0.0280]T, while in the second segment (n ∈ [100,349]) the AR coefficients are a1 = [0.1200, 0.0245, -0.2787, -0.0693]T. In the third segment (n ∈ [350, 499]), a0 is in act. Figure 4 shows the true signal along with the group Lasso, group SCAD, and DP-based estimates of the TV-AR 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 block-coordinate descent algorithm of Section 4 was used to solve (12) up to a maximum of 103 iterations or when <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M116">View MathML</a> 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M117">View MathML</a> across time. Clearly, the instants where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M118">View MathML</a> represents the change points retrieved. While the DP has detected change-points at <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M119">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M120">View MathML</a> the group Lasso and the group SCAD have detected changes at <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M121">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M122">View MathML</a> confirming that the results returned by the DP programming, by the group Lasso, and by group SCAD are comparable. Since {yn} was generated as in (1), any L out of N + 1 vectors <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M93">View MathML</a> are linearly independent almost surely, which means that Proposition 3 can be invoked to ensure uniqueness of the group Lasso TV-AR model estimate. Notice from Figure 4 and Figure 5 that the <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M123">View MathML</a> 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.

thumbnailFigure 4. Synthetic data. From top to bottom: True signal, estimated of the TV-AR coefficients by group Lasso, group SCAD, and DP, respectively.

thumbnailFigure 5. Synthetic data. From top to bottom: True TV-AR 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 TV-AR model with L = 1 is adopted, and λ is selected to be λ*/10. Figure 7 depicts the TV-AR coefficients estimated by the group Lasso and the group SCAD. The estimated <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M117">View MathML</a> over time is displayed in Figure 8. Observe that the group Lasso captures the correct changes at time instants <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M124">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M125">View MathML</a> along with a small false change at time instant nf = 2770. Instead, the group SCAD retrieves two changes at time instants <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M124">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M125">View MathML</a> Notice also that group SCAD exhibits a reduced bias in the estimated <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M117">View MathML</a> which results in better estimates of the TV-AR coefficients.

thumbnailFigure 6. Piano sound comprising three monochromatic notes.

thumbnailFigure 7. Piano sound. Estimated of the TV-AR coefficients by group Lasso (top), and group SCAD (bottom).

thumbnailFigure 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 TV-AR model with L = 8 is adopted, and λ is selected to be λ*/10. The change of vocoid in the diphthong occurs approximately at time instant n1 = 1500, while the /o/ occurs approximately at n2 = 3000. Figure 10 shows the TV-AR 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 ║dn2 over time. The group SCAD detects four segments with change points at <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M126">View MathML</a>, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M127">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M128">View MathML</a>. 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 post-processing via either DP as advocated in [17], or, by simply peak picking.

thumbnailFigure 9. Speech signal: /ai//o/.

thumbnailFigure 10. Speech signal. Estimated of the TV-AR coefficients by group Lasso (top), and group SCAD (bottom).

thumbnailFigure 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 TV-AR speech segmentation [3-5,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 [3-5] 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.

thumbnailFigure 12. French speech signal and changes detected by the group SCAD (dot) and GLR (dash) change TV-AR 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 [3-5]. 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 [3-5], 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 TV-AR 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M129">View MathML</a>, let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M130">View MathML</a> denotes the LS estimates of the AR model in the kth segment, i.e., <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M131">View MathML</a>. The SPE is defined as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M132">View MathML</a>, and represents the error in approximating the original signal {yn} with a TV-AR model exhibiting abrupt changes at instants <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M129">View MathML</a>. The GLR segmentation entails SPEglr = 0.2638, while SPEg-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 block-coordinate descent algorithm is not available, a simulated test assessing its converge speed is performed. Figure 13 depicts the normalized step size amplitude variations, namely <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M133">View MathML</a>, across the iteration index (i) for the first LLA of the group SCAD (which amounts to a group Lasso since <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M114">View MathML</a> 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 <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M133">View MathML</a> 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 block-coordinate descent (implemented with Matlab) for the problem at hand, whereas advanced Bayesian techniques might require several hours (see e.g., [3]).

thumbnailFigure 13. Speed of the block-coordinate descent for the first LLA of the group SCAD.

8. Concluding remarks

Novel estimators were developed in this article for identification of piecewise-constant TV-AR 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 block-coordinate descent iterations. The latter incurs computational burden that scales linearly with the number of data samples, thus being particularly attractive for large-size problems. Regularization tuning issues are discussed along with conditions for uniqueness of the estimated piecewise-constant AR model. An alternative group smoothly-clipped 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 piecewise-constant 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 ARX-models using sum-of-norms 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 first-order optimality condition for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> to be the (unconstrained) minimum of (12), is that the subgradient of J(d) in (13) evaluated at <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> contains the zero vector [[36], p. 92]; i.e.,

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M134">View MathML</a>

(39)

Defining <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M181','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M181">View MathML</a> as w := XT(Xd - y) and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M135">View MathML</a>, with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M182','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M182">View MathML</a> for n = 0,1,..., N, the subgradient of J(d) evaluated at <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> is given by

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M136','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M136">View MathML</a>

(40)

where <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M137">View MathML</a> and

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M138','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M138">View MathML</a>

(41)

with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M189','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M189">View MathML</a> such that ║sn2 ≤ 1. Using (40) and (41), Equation (39) translates to the following conditions:

(c1) w0 = 0l; and,

(c2)<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M139','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M139">View MathML</a>

The change-free solution corresponds to having <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M87">View MathML</a> for n = 1,..., N. Thus, (c1) implies that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M140','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M140">View MathML</a>, which is uniquely satisfied by <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M86">View MathML</a>, since X0 has full column rank. Hence, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M86">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M87">View MathML</a> for n = 1,..., N hold if and only if (c2) is satisfied, which corresponds to ║wn2 ≤ λ for n = 1,..., N. Since <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M141','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M141">View MathML</a>, condition (c2) is satisfied if and only if <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M142','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M142">View MathML</a>.

Appendix 2: proof of lemma 1

Observe that

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M143','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M143">View MathML</a>

(42)

Notice that the first sub-sum in (42) comprises <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M183','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M183">View MathML</a> rank-1 matrices, the last sub-sum comprises s1 rank-1 matrices, while the gth sub-sum comprises <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M184','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M184">View MathML</a> rank-1 matrices for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M185','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M185">View MathML</a>. Since is such that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M186','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M186">View MathML</a>, and |sj - sk| ≥ L for each j k, and any L out of N + 1 vectors <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M144','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M144">View MathML</a> are linearly independent, each of the summands in (42) has rank L. Thus, it is possible to find L linearly independent vectors <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M145','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M145">View MathML</a> such that the first sub-sum in (42) equals to <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M146','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M146">View MathML</a> with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M147','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M147">View MathML</a>. Analogously, it is possible to find L linearly independent vectors <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M148','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M148">View MathML</a> such that the gth sub-sum in (42) can be written as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M149','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M149">View MathML</a> with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M150','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M150">View MathML</a> for <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M185','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M185">View MathML</a>. Finally, it is possible to find L linearly independent vectors <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M151','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M151">View MathML</a> such that the last sub-sum in (42) can be written as <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M152','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M152">View MathML</a> with <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M153','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M153">View MathML</a>. Thus, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M154','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M154">View MathML</a>, and since <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M155','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M155">View MathML</a> are (|| + 1) linearly independent vectors, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178">View MathML</a> has full-column rank.

Appendix 3: proof of proposition 3

Suppose that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M156','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M156">View MathML</a> are two solutions of (12) with the same g-support . Let <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178">View MathML</a> denotes the matrix obtained by selecting the columns of X relative to the nonzero indexes dictated by the set <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M187','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M187">View MathML</a>. The minimization of (13) over the vector having g-support in amounts to

<a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M157','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M157">View MathML</a>

(43)

From Lemma 1, <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M178">View MathML</a> has full column rank which implies that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M158','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M158">View MathML</a> is strictly convex, and so is <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M159','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M159">View MathML</a>. Thus, (43) admits a unique solution.

Since both <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M91">View MathML</a> and <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M156','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M156">View MathML</a> are g-supported in by hypothesis, and their restrictions to <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M187','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M187">View MathML</a> are equal to û, it follows readily that <a onClick="popup('http://asp.eurasipjournals.com/content/2012/1/70/mathml/M160','MathML',630,470);return false;" target="_blank" href="http://asp.eurasipjournals.com/content/2012/1/70/mathml/M160">View MathML</a>.

Acknowledgements

This work was supported by MURI (AFOSR FA9550-10-1-0567) grant.

References

  1. P Stoica, RL Moses, Introduction to Spectral Analysis (Prentice-Hall, New Jersey, 1997)

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

  3. N Dobigeon, J-Y 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)

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

  5. 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)

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

  7. M Basseville, IV Nikiforov, Detection of Abrupt Changes: Theory and Application (Prentice-Hall, Englewood Cliffs, NJ, USA, 1993)

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

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

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

  11. PJ Brockwell, RA Davis, Time Series: Theory and Methods (Springer-Verlag, New York, NY, USA, 1990)

  12. L Boysen, A Kempe, V Liebscher, A Munk, O Wittich, Consistencies and rates of convergence of jump-penalized least-squares estimators. Annals Stat 37(1), 157–183 (2009)

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

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

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

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

  17. Z Harchaoui, C Levy-Leduc, Catching change-points with Lasso. Proceedings of the Advanced Neural Information Processes Systems (Vancouver, Canada, 2008) 20, pp. 161–168

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

  19. 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)

  20. 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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  35. 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)

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