Kevin Bleakley
email: kevbleakley@gmail.com

Jean-Philippe Vert
email: jean-philippe.vert@mines.org

The group fused Lasso for multiple change-point detection

published or not. The documents may come

The group fused Lasso for multiple change-point detection

Introduction

Finding the place (or time) where most or all of a set of one-dimensional signals (or profiles) jointly change in some specific way is an important question in several fields. A common situation is when we want to find change-points in a multidimensional signal, e.g., in audio and image processing [START_REF] Desobry | An online kernel change detection algorithm[END_REF][START_REF] Harchaoui | A regularized kernel-based approach to unsupervised audio segmentation[END_REF], to detect intrusion in computer networks [START_REF] Tartakovsky | A novel approach to detection of intrusions in computer networks via adaptive sequential and batch-sequential changepoint detection methods[END_REF][START_REF] Lévy-Leduc | Detection and localization of change-points in high-dimensional network traffic data[END_REF], or in financial and economics time series analysis [START_REF] Talih | Structural learning with time-varying components: tracking the crosssection of financial time series[END_REF]. Another important situation is when we are confronted with several 1-dimensional signals which we believe share common change-points, e.g., genomic profiles of a set of patients. The latter application is increasingly important in biology and medicine, in particular for the detection of copy-number variation along the genome [START_REF] Picard | A statistical approach for array CGH data analysis[END_REF], or the analysis of microarray and genetic linkage studies [START_REF] Zhang | Detecting simultaneous change-points in multiple sequences[END_REF]. The common thread in biological applications is the search for data patterns shared by a set of individuals, such as cancer patients, at precise places on the genome; in particular, sudden changes in measured values. As opposed to the segmentation of multidimensional signals such as speech, where the dimension is fixed and collecting more data means having longer profiles, the length of signals in genomic studies (i.e., the number of probes measured along the genome) is fixed for a given technology while the number of signals (i.e., the number of individuals) can increase when we collect data about more patients. From a statistical point of view, it is therefore of interest to develop methods that identify multiple change-points shared by several signals that can benefit from increasing the number of signals.

There exists a vast literature on the change-point detection problem [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF][START_REF] Brodsky | Nonparametric Methods in Change-Point Problems[END_REF]. Here we focus on computationally efficient methods to segment a multidimensional signal by approximating it with a piecewiseconstant one, using quadratic error criteria. It is well-known that, in this case, the optimal segmentation of a p-dimensional signal of length n into k segments can be obtained in O(n 2 pk) by dynamic programming [START_REF] Yao | Estimating the number of change-points via schwarz criterion[END_REF][START_REF] Birgé | Gaussian model selection[END_REF][START_REF] Lavielle | Detection of multiple change-points in multivariate time series[END_REF]. However, the quadratic complexity in n is prohibitive in applications such as genomics, where n can be in the order of 10 5 to 10 7 with current technology. An alternative to such global procedures, which estimate change-points as solutions of a global optimization problem, are fast local procedures such as binary segmentation [START_REF] Vostrikova | Detection of disorder in multidimensional stochastic processes[END_REF], which detect breakpoints by iteratively applying a method for single change-point detection to the segments obtained after the previous change-point is detected. While such recursive methods can be extremely fast, in the order of O(np log(k)) when the single change-point detector is O(np), quality of segmentation is questionable when compared with global procedures [START_REF] Lavielle | Adaptive detection of multiple change-points in asset price volatility[END_REF].

For p = 1 (a single signal), an interesting alternative to these global and local procedures is to express the optimal segmentation as the solution of a convex optimization problem, using the (convex) total variation instead of the (non-convex) number of jumps to penalize a piecewise-constant function in order to approximate the original signal [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF][START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF]. The resulting piecewise-constant approximation of the signal, defined as the global minimum of the objective function, benefits from theoretical guaranties in terms of correctly detecting change-points [START_REF] Harchaoui | Catching change-points with lasso[END_REF][START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF], and can be implemented efficiently in O(nk) or O(n log(n)) [START_REF] Friedman | Pathwise coordinate optimization[END_REF][START_REF] Harchaoui | Catching change-points with lasso[END_REF][START_REF] Hoefling | A path algorithm for the Fused Lasso Signal Approximator[END_REF].

In this paper we propose an extension of total-variation based methods for single signals to the multidimensional setting, in order to approximate a multidimensional signal with a piecewise-constant signal with multiple change-points. We define the approximation as the solution of a convex optimization problem which involves a quadratic approximation error penalized by the sum of the Euclidean norms of the multidimensional increments of the function. The problem can be reformulated as a group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], which we show how to solve exactly and efficiently. Alternatively, we provide an approximate yet often computationally faster solution to the problem using a group LARS procedure [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. In the latter case, using the particular structure of the design matrix, we can find the first k change-points in O(npk), thus extending the method of [START_REF] Harchaoui | Catching change-points with lasso[END_REF] to the multidimensional setting.

Unlike most previous theoretical investigations of change-point methods (e.g., [START_REF] Harchaoui | Catching change-points with lasso[END_REF][START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF]), we are not interested in the case where the dimension p is fixed and the length of the profiles n increases, but in the opposite situation where n is fixed and p increases. Indeed, this corresponds to the case in genomics where, for example, n would be the fixed number of probes used to measure a signal along the genome, and p the number of samples or patients analyzed. We want to design a method that benefits from increasing p in order to identify shared change-points, even though the signal-to-noise ratio may be very low within each signal. As a first step towards this question, we give conditions under which our method is able to consistently identify a single change-point as p increases. We also show by simulation that the method is able to correctly identify multiple change-points as p → +∞, validating its relevance in practical settings.

The paper is organized as follows. After fixing notation in Section 2, we present the group fused Lasso method in Section 3. We propose two efficient algorithms to solve it in Section 4, and discuss its theoretical properties in Section 5. Lastly, we provide an empirical evaluation of the method and a comparison with other methods in the study of copy number variations in cancer in Section 6. A preliminary version of this paper was published in [START_REF] Vert | Fast detection of multiple change-points shared by many signals using group LARS[END_REF].

Notation

For any two integers u ≤ v, we denote by [u, v] the interval {u, u + 1, . . . , v}. For any u × v matrix M we note M i,j its (i, j)-th entry, and M = u i=1 v j=1 M 2 i,j its Frobenius norm (or Euclidean norm in the case of vectors). For any subsets of indices

A = a 1 , . . . , a |A| ∈ [1, u] |A| and B = b 1 , . . . , b |B| ∈ [1, v] |B| , we denote by M A,B the |A| × |B| matrix with entries M a i ,b j for (i, j) ∈ [1, |A|] × [1, |B|].
For simplicity we will use • instead of [1, u] or [1, v], i.e., A i,• is the i-th row of A and A •,j is the j-th column of A. We note 1 u,v the u × v matrix of ones, and I p the p × p identity matrix.

Formulation

We consider p real-valued profiles of length n, stored in an n × p matrix Y . The i-th profile Y •,i = (Y 1,i , . . . , Y n,i ) is the i-th column of Y . We model each profile as a piecewise-constant signal corrupted by noise, and assume that change-point locations tend to be shared across profiles. Our goal is to detect these shared change-points, and benefit from the possibly large number p of profiles to increase the statistical power of change-point detection.

Segmentation with a total variation penalty

When p = 1 (a single profile), a popular method to find change-points in a signal is to approximate it by a piecewise-constant function using a quadratic error criterion, i.e., to solve

min U ∈R n Y -U 2 subject to n-1 i=1 δ(U i+1 -U i ) ≤ k , ( 1
)
where δ is the Dirac function, equal to 0 if its argument is null, 1 otherwise. In other words, (1) expresses the best approximation of Y by a piecewise-constant profile U with at most k jumps. It is well-known that (1) can be solved in O(n 2 k) by dynamic programming [START_REF] Yao | Estimating the number of change-points via schwarz criterion[END_REF][START_REF] Birgé | Gaussian model selection[END_REF][START_REF] Lavielle | Detection of multiple change-points in multivariate time series[END_REF]. Although very fast when n is of moderate size, the quadratic dependency in n renders it impractical in current computers when n reaches millions or more, which is often the case in many application such as segmentation of genomic profiles.

An alternative to the combinatorial optimization problem ( 1) is to relax it to a convex optimization problem, by replacing the number of jumps by the convex total variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], i.e., to consider:

min U ∈R n 1 2 Y -U 2 + λ n-1 i=1 | U i+1 -U i | . (2)
For a given λ > 0, the solution U ∈ R n of ( 2) is again piecewise-constant. Recent work has shown that (2) can be solved much more efficiently than (1): [START_REF] Friedman | Pathwise coordinate optimization[END_REF] proposed a fast coordinate descent-like method, [START_REF] Harchaoui | Catching change-points with lasso[END_REF] showed how to find the first k change-points iteratively in O(nk), and [START_REF] Hoefling | A path algorithm for the Fused Lasso Signal Approximator[END_REF] proposed a O(n ln(n)) method to find all change-points. Adding penalties proportional to the ℓ 1 or ℓ 2 norm of U to (2) does not change the position of the change-points detected [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF], and the capacity of TV denoising to correctly identify change-points when n increases has been investigated in [START_REF] Harchaoui | Catching change-points with lasso[END_REF][START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF].

Here, we propose to generalize TV denoising to multiple profiles by considering the following convex optimization problem, for Y ∈ R n×p :

min U ∈R n×p 1 2 Y -U 2 + λ n-1 i=1 U i+1,• -U i,• . (3)
The second term in (3) can be considered a multidimensional TV: it penalizes the sum of Euclidean norms of the increments of U , seen as a time-dependent multidimensional vector, and reduces to the classical 1-dimensional TV when p = 1. Intuitively, when λ increases, this penalty will enforce many increment vectors U i+1,• -U i,• to collapse to 0, just like the total variation in (2) in the case of 1dimensional signals. This implies that the positions of non-zero increments will be the same for all profiles. As a result, the solution to (3) provides an approximation of the profiles Y by an n × p matrix of piecewise-constant profiles U which share change-points. While (3) is a natural multidimensional generalization of the classical TV denoising method (2), we more generally investigate the following variant:

min U ∈R n×p 1 2 Y -U 2 + λ n-1 i=1 U i+1,• -U i,• d i , (4)
where (d i ) i=1,...,n-1 are position-dependant weights which affect the penalization of the jump differently at different positions. While (4) boils down to (3) for uniform weights d i = 1, i = 1, . . . , n -1, we will see that the unweighted version suffers from boundary effects and that position-dependent schemes such as:

∀i ∈ [1, n -1], d i = n i(n -i) , (5)
are both theoretically and empirically better choices.

To illustrate the grouping effect of the penalty in (4), Figure 1 compares the segmentation of three simulated profiles obtained with and without enforced sharing of change-points across profiles. We simulated three piecewise-constant signals corrupted by independent additive Gaussian noise. All profiles have length 500 and share the same 5 change-points, though with different amplitudes, at positions 38, 139, 268, 320 and 397. On the left-hand side, we show the first 5 change-points captured by TV denoising with weights (5) applied to each signal independently. On the right, we show the first 5 change-points captured by formulation [START_REF] Lévy-Leduc | Detection and localization of change-points in high-dimensional network traffic data[END_REF]. We see that the latter formulation finds the correct change-points, whereas treating each profile independently leads to errors. For example, the first two change-points have a small amplitude in the second profile and are therefore very difficult to detect from the profile only, while they are very apparent in the first and third profiles.

Reformulation as a group Lasso problem

It is well-known that the 1-dimensional TV denoising problem (2) can be reformulated as a Lasso regression problem by an appropriate change of variable [START_REF] Harchaoui | Catching change-points with lasso[END_REF]. We now show that our generalization (4) can be reformulated as a group Lasso regression problem, which will be convenient for theoretical analysis and implementation [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. To this end, we make the change of variables (β, γ) ∈ R (n-1)×p × R 1×p given by:

γ = U 1,• , β i,• = U i+1,• -U i,• d i for i = 1, . . . , n -1 .
In other words d i β i,j is the jump between the i-th and the (i + 1)-th positions of the j-th profile. We immediately get an expression for U as a function of β and γ:

U 1,• = γ , U i,• = γ + i-1 j=1 d j β j,• for i = 2, . . . , n .
This can be rewritten in matrix form as

U = 1 n,1 γ + Xβ , ( 6
)
where X is the n × (n -1) matrix with entries X i,j = d j for i > j, and 0 otherwise. Making this change of variable, we can re-express (4) as follows:

min β∈R (n-1)×p ,γ∈R 1×p 1 2 Y -Xβ -1 n,1 γ 2 + λ n-1 i=1 β i,• . (7)
For any β ∈ R (n-1)×p , the minimum in γ is attained with γ = 1 1,n (Y -Xβ)/n. Plugging this into (7), we get that the matrix of jumps β is solution of

min β∈R (n-1)×p 1 2 Ȳ -Xβ 2 + λ n-1 i=1 β i,• , (8)
where Ȳ and X are obtained from Y and X by centering each column. Equation ( 8) is now a classical group Lasso regression problem [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], with a specific design matrix X and groups of features corresponding to the rows of the matrix β. The solution β of ( 8) is related to the solution U of our initial problem (4) by equation ( 6).

Implementation

Although [START_REF] Lévy-Leduc | Detection and localization of change-points in high-dimensional network traffic data[END_REF] and [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF] are convex optimization problems that can in principle be solved by generalpurpose solvers [START_REF] Boyd | Convex Optimization[END_REF], we want to be able to work in dimensions that reach millions or more, making this computationally difficult. In particular, the design matrix X in ( 8) is a non-sparse matrix of size n × (n -1), and cannot even fit in a computer's memory when n is large. Moreover, we would ideally like to obtain solutions for various values of λ, corresponding to various numbers of change-points, in order to be able to select the optimal number of change-points using statistical criteria. In the single profile case (p = 1), fast implementations in O(nk) or O(n ln n) have been proposed [START_REF] Friedman | Pathwise coordinate optimization[END_REF][START_REF] Harchaoui | Catching change-points with lasso[END_REF][START_REF] Hoefling | A path algorithm for the Fused Lasso Signal Approximator[END_REF]. However, none of these methods is applicable directly to the p > 1 setting since they all rely on specific properties of the p = 1 case, such as the fact that the solution is piecewise-affine in λ and that the set of change-points is monotically decreasing with λ.

In this section we propose two algorithms to respectively exactly or approximately solve (4) efficiently. We adopt the algorithms suggested by [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] to solve the group Lasso problem (8) and show how they can be implemented very efficiently in our case due to the particular structure of the regression problem. We have placed in Annex A several technical lemmas which show how to efficiently perform several operations with the given design matrix X that will be used repeatedly in the implementations proposed below.

Exact solution by block coordinate descent

A first possibility to solve the group Lasso problem ( 8) is to follow a block coordinate descent approach, where each group is optimized in turn with all other groups fixed. It can be shown that this strategy converges to the global optimum, and is reported to be stable and efficient [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF][START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF]. As shown by [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], it amounts to iteratively applying the following equation to each block i = 1, . . . , n -1 in turn, until convergence:

β i,• ← 1 γ i 1 - λ S i + S i , (9)
where

γ i = X•,i 2 = i(n-i)d 2 i /n and S i = X⊤ •,i Ȳ -Xβ -i
, and where β -i denotes the (n-1)×p matrix equal to β except for the i-th row β -i i,• = 0. The convergence of the procedure can be monitored by the Karush-Kuhn-Tucker (KKT) conditions:

-X⊤ •,i Ȳ -Xβ + λβ i,• β i,• = 0 ∀β i,• = 0 , -X⊤ •,i Ȳ -Xβ ≤ λ ∀β i,• = 0 . ( 10
)
Since the number of blocks n can be very large and we expect only a fraction of non-zero blocks at the optimum (corresponding to the change-points), we implemented this block coordinate descent with an active set strategy. In brief, a set of active groups A corresponding to non-zero groups is maintained, and the algorithm alternates between optimizing β over the active groups in A and updating A by adding or removing groups based on violation of the KKT conditions. The resulting pseudo-code is shown in Algorithm 1. The inner loop (lines 3-7) corresponds to the optimization of β on the current active groups, using iteratively block coordinate descent [START_REF] Brodsky | Nonparametric Methods in Change-Point Problems[END_REF]. After convergence, groups that have been shrunk to 0 are removed from the active set (line 8), and the KKT conditions are checked outside of the active set (lines 9-10). If they are not fulfilled, the group that most violates the conditions is added to the active set (line 11), otherwise the current solution satisfies all KKT conditions and is therefore the global optimum (line 13).

Although it is difficult to estimate the number of iterations needed to reach convergence for a certain level of precision, we note that by Lemma 5 (Annex A), computation of X⊤ Ȳ in line 1 can be done in O(np), and each group optimization iteration (lines 3-7) requires computing X⊤

•,i X •,A (line 5), done in O(|A|) (see Lemma 6 in Annex A), then computing S i (line 5) in O(|A|p) and soft-thresholding (line 6) in O(p). The overall complexity of each group optimization iteration is therefore O(|A|p). Since each group in A must typically be optimized several times, we expect complexity that is at least quadratic in |A| and linear in p for each optimization over an active set A (lines 3-7). To check optimality of a solution after optimization over an active set A, we need to compute X⊤ Xβ (line 9) which takes O(np) (see Lemma 7, Annex A). Although it is difficult to upper bound the number of iterations needed to optimize over A, this shows that a best-case complexity to find k change-points, if we correctly add groups one by one to the active set, would be O(npk) to check k times the KKT conditions and find the next group to add, and O(pk 3 ) in total if each optimization over an active set A is in O(p|A| 2 ). In Section 6, we provide some empirical results on the behavior of this block coordinate descent strategy.

Group fused LARS implementation

Since exactly solving the group Lasso with the method described in Section 4.1 can be computationally intensive, it may be of interest to find fast, approximate solutions to [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. We propose to implement a strategy based on the group LARS, proposed in [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] as a good way to approximately find the regularization path of the group Lasso. More precisely, the group LARS approximates the solution path of ( 8) with a piecewise-affine set of solutions and iteratively finds change-points. The resulting algorithm is presented here as Algorithm 2, and is intended to approximately solve [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. Change-points are added one by one (lines 4 and 8), and for a given set of change-points the solution moves straight along a descent direction (line 6) with a given step (line 7) until a new change-point is added (line 8). We refer to [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF] for more details and justification for this algorithm.

While the original group LARS method requires storage and manipulation of the design matrix [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], implausible for large n here, we can again benefit from the computational tricks provided in Annex A Pick i ∈ A.

5:

Compute S i ← C i,• -X⊤ •,i Xβ -i . 6:
Update β i,• according to (9). Remove inactive groups: A ← A\ {i ∈ A : β i,• = 0}.

9:

Check KKT: S ← C -X⊤ Xβ.

10:

û ← argmax i / ∈A S i,• 2 , M = S û,• 2 .
11:

if M > λ 2 then 12:
Add a new group: A ← A ∪ {û}. First change-point : û ← argmin j∈[1,n-1] ĉj,• , A ← {û}.

5:

end if

6: Descent direction: compute w ← X⊤ •,A X•,A -1
ĉA,• , then a = X⊤ XA w.

7:

Descent step: for each u ∈ [1, n -1] \A, find if it exists the smallest positive solution α u of the second-order polynomial in α:

ĉu,• -αa u,• 2 = ĉv,• -αa v,• 2 ,
where v is any element of A.

8:

Next change-point:

û ← argmin j∈[1,n-1] ĉj,• , A ← A ∪ {û}. 9:
Update ĉ ← ĉα ûa. 10: end for

Theoretical analysis

In this section, we study theoretically to what extent the estimator (4) recovers correct change-points. The vast majority of existing theoretical results for offline segmentation and change-point detection con-sider the setting where p is fixed (usually p = 1), and n increases (e.g., [START_REF] Harchaoui | A regularized kernel-based approach to unsupervised audio segmentation[END_REF]). This typically corresponds to cases where we can sample a continuous signal with increasing density, and wish to locate more precisely the underlying change-points as the density increases.

We propose a radically different analysis, motivated notably by applications in genomics. Here, the length of profiles n is fixed for a given technology, but the number of profiles p can increase when more samples or patients are collected. The property we would like to study is then, for a given change-point detection method, to what extent increasing p for fixed n allows us to locate more precisely the changepoints. While this simply translates our intuition that increasing the number of profiles should increase the statistical power of change-point detection, and while this property was empirically observed in [START_REF] Zhang | Detecting simultaneous change-points in multiple sequences[END_REF], we are not aware of previous theoretical results in this setting. In particular we are interested in the consistency of our method, in the sense that it should correctly detect the true change-points if enough samples are available.

Consistent estimation of a single change-point

As a first step towards the analysis of this "fixed n increasing p" setting, let us assume that the observed centered profiles Ȳ are obtained by adding noise to a set of profiles with a single shared change-point between positions u and u + 1, for some u ∈ [1, n -1]. In other words, we assume that

Ȳ = Xβ * + W ,
where β * is an (n -1) × p matrix of zeros except for the u-th row β * u,• , and W is a noise matrix whose entries are assumed to be independent and identically distributed with respect to a centered Gaussian distribution with variance σ 2 . In this section we study the probability that the first change-point found by our procedure is the correct one, when p increases. We therefore consider an infinite sequence of jumps β * u,i i≥1

, and letting β2 p = 1 p p i=1 (β * u,i ) 2 , we assume that β2 = lim p→∞ β2 p exists and is finite. We first characterize the first selected change-point as p increases.

G i = d 2 i i(n -i) n σ 2 + β2 d 2 i d 2 u n 2 × i 2 (n -u) 2 if i ≤ u , u 2 (n -i) 2 otherwise. ( 11
)
When p → +∞, the first change-point selected by the group fused Lasso ( 4) is in argmax i∈[1,n-1] G i with probability tending to 1.

Proof of this result is given in Annex B. From it we easily deduce conditions under which the first change-point is correctly found with increasing probability as p increases. Let us first focus on the unweighted group fused Lasso (3), corresponding to the setting d i = 1 for i = 1, . . . , n -1.

Theorem 2. Let α = u/n be the position of the change-point scaled in the interval [0, 1], and

σ2 α = n β2 (1 -α) 2 (α -1 2n ) α -1 2 -1 2n . ( 12
)
If σ 2 < σ2 α , the probability that the first change-point selected by the unweighted group fused Lasso ( 3) is the correct one tends to 1 as p → +∞. When σ 2 > σ2 α , it is not the correct one with probability tending to 1.

This theorem, the proof of which can be found in Annex C, deserves several comments.

• To detect a change-point at position u = αn, the noise level σ 2 must not be larger than the critical value σ2 α given by ( 12), hence the method is not consistent for all positions. σ2 α decreases monotonically from α = 1/2 to 1, meaning that change-points near the boundary are more difficult to detect correctly than change-points near the center. The most difficult change-point is the last one (u = n -1) which can only be detected consistently if σ 2 is smaller than

σ2 1-1/n = 2 β2 n + o(n -1
).

• For a given level of noise σ 2 , change-point detection is asymptotically correct for any α ∈ [ǫ, 1ǫ], where ǫ satisfies σ 2 = σ2 1-ǫ , i.e.,

ǫ = σ 2 2n β2 + o(n -1/2 ) .
This shows in particular that increasing the profile length n increases the relative interval (as a fraction of n) where change-points are correctly identified, and that we can get as close as we want to the boundary for n large enough.

• When σ 2 < σ2 α , the correct change-point is found consistently when p increases, showing the benefit of the accumulation of many profiles.

Theorem 2 shows that the unweighted group fused Lasso (3) suffers from boundary effects, since it may not correctly identify a single change-points near the boundary is the noise is too large. In fact, Lemma 1 tells us that if we miss the correct change-point position, it is because we estimate it more towards the middle of the interval (see proof of Theorem 2 for details). The larger the noise, the more biased the procedure is. We now show that this issue can be fixed when we consider the weighted group fused Lasso (4) with well-chosen weights.

Theorem 3. The weighted group fused Lasso [START_REF] Lévy-Leduc | Detection and localization of change-points in high-dimensional network traffic data[END_REF] with weights given by ( 5) correctly finds the first change-point at any position with probability tending to 1 as p → +∞.

The proof of Theorem 3 is postponed to Annex D. It shows that the weighting scheme (5) cancels the effect of the noise and allows us to consistently estimate any change-point, independently of its position in the signal, as the number of signals increases.

Consistent estimation of a single change-point with fluctuating position

An interesting variant of the problem of detecting a change-point common to many profiles is that of detecting a change-point with similar location in many profiles, allowing fluctuations in the precise location of the change-point. This can be modeled by assuming that the profiles are random, and that the i-th profile has a single change-point of value β i at position U i , where (β i , U i ) i=1,...,p are independent and identically distributed according to a distribution P = P β ⊗ P U (i.e., we assume β i independent from U i ). We denote β2 = E P β β 2 and

p i = P U (U = i) for i ∈ [1, n -1]. Assuming that the support of P U is [a, b] with 1 ≤ a ≤ b ≤ n -1,
the following result extends Theorem 2 by showing that the first change-point discovered by the unweighted group fused Lasso is in the support of P U under some condition on the noise level, while the weighted group fused Lasso correctly identifies a change-point in the support of P U asymptotically without conditions on the noise.

→ +∞. If 1/2 < α m or α M < 1/2, let σ2 P U = n β2 (1 -Eα) 2 + var(α) 2 ×    αm-1 2n αm-1 2 -1 2n if α m > 1 2 , 1-1 2n -α M 1 2 -α M -1 2n if α M < 1 2 . ( 13
)
The probability that the first selected change-point is in the support of P U tends to 1 when σ 2 < σ2 P U . When σ 2 > σ2 P U , it is outside of the support of P U with probability tending to 1.

2. The weighted group fused Lasso (4) with weights given by ( 5) finds the first change-point in the support of P U with probability tending to 1 as p → +∞, independently of σ 2 and of the support of P U .

This theorem, the proof of which is postponed to Annex E, illustrates the robustness of the method to fluctuations in the precise position of the change-point shared between profiles. Although this situation rarely occurs when we are considering classical multidimensional signals such as financial time series or video signals, it is likely to be the rule when we consider profiles coming from different biological samples, where for example we can expect frequent genomic alterations at the vicinity of important oncogenes or tumor suppressor genes. Although the theorem only gives a condition on the noise level to ensure that the selected change-point lies in the support of the distribution of change-point locations, a precise estimate of the location of the selected change-point as a function of P U , which generalizes Lemma 1, is given in the proof.

The case of multiple change-points

While the theoretical results presented above focus on the detection of a single change-point, the real interest of the method is to estimate multiple change-points. The extension of Theorem 2 to this setting is, however, not straightforward and we postpone it for future efforts. We conjecture that the group fused Lasso estimator can, under certain conditions, consistently estimate multiple change-points. More precisely, in order to generalize the proof of Theorem 2, we must analyze the path of the vectors (ĉ i,• ), and check that, for some λ in (3) or ( 4), they reach their maximum norm precisely at the true changepoints. The situation is more complicated than in the single change-point case since, in order to fulfill the KKT optimality conditions, the vectors (ĉ i,• ) must hit a hypersphere at each correct change-point, and must remain strictly within the hypersphere between consecutive change-points. This can probably be ensured if the noise level is not too high (like in the single change-point case), and if the positions corresponding to successive change-points on the hypersphere are far enough from each other, which could be ensured if two successive change-points are not too close to each other, and are in sufficiently different directions. Although the weighting scheme (5) ensures consistent estimation of the first changepoint independently of the noise level, it may however not be sufficient to ensure consistent estimation of subsequent change-points.

Although we propose no theoretical results besides these conjectures for the case of multiple changepoints, we provide experimental results below that confirm that, when the noise is not too large, we can indeed correctly identify several change-points, with probability of success increasing to 1 as p increases.

Estimating the number of change-points

The number of change-points detected by the group fused Lasso in the multidimensional signal depends on the choice of λ in (3) and (4). In practice, we propose the following scheme in order to estimate a segmentation and the number of change-points. We try to select a λ that over-segments the multidimensional signal, that is, finds more change-points that we would normally expect for the given type of signal or application. Then, on the set of k change-points found, we perform post-processing using a simple least-squares criteria. Briefly, for each given subset of k ′ ≤ k change-points, we approximate each signal between successive change-points with the mean value of the points in that interval; then, we calculate the total sum of squared errors (SSE) between the set of real signals and these piecewiseconstant approximations to them. Though it may appear computationally intensive or even impossible to do this for all subsets of k ′ ≤ k change-points, a dynamic programming strategy (e.g., [START_REF] Picard | A statistical approach for array CGH data analysis[END_REF]) means that the best subset of k ′ ≤ k change-points can be calculated for all k ′ ∈ {1, . . . , k} in O(k 3 ).

It then remains to choose the "best" k ′ ∈ {1, . . . , k} using, for example, a model-selection strategy. The optimal SSE for k ′ + 1 (which we may call SSE(k ′ + 1) to ease notation), will be smaller than SSE(k ′ ) but at a certain point, adding a further change-point will have no physical reality, it only improves the SSE due to random noise. Here, we implemented a method proposed in [START_REF] Lavielle | Using penalized contrasts for the change-point problem[END_REF][START_REF] Picard | A statistical approach for array CGH data analysis[END_REF] where we first normalize the SSE for k ′ = 1, . . . , k into a score J(k ′ ) such that J(1) = k and J(k) = 1, in such a way the J(k ′ ) has an average slope of -1 between 1 and k; we then try to detect a kink in the curve by calculating the discrete second derivative of J(k ′ ), and selecting the k ′ after which this second derivative no longer rises above a fixed threshold (typically 0.5).

Experiments

In this section we test the group fused Lasso on several simulated and real data sets. All experiments were run under Linux on a machine with two 4-core Intel Xeon 3.16GHz processors and a total of 16Gb of RAM. We have implemented the group fused Lasso in MATLAB; the package GFLseg is available for download1 .

Speed trials

In a first series of experiments, we tested the behavior of the group fused Lasso in terms of computational efficiency. We simulated multidimensional profiles with various lengths n between 2 4 and 2 23 , various dimensions p between 1 and 2 15 , and various number of shared change-points k between 1 and 2 7 . In each case, we first ran the iterative weighted group fused LARS (Section 4.2) to detect successive change-points, and recorded the corresponding λ values. We then ran the exact group fused Lasso implementation by block coordinate descent (Section 4.1) on the same λ values. Figure 2 shows speed with respect to increasing one of p, n and k while keeping the other two variables fixed, for both implementations. The axes are log-log, so the slope gives the exponent of the complexity (resp. n, p and k). For the weighted group fused LARS, linearity is clearest for k, whereas for n and p, the curves are initially sub-linear, then slightly super-linear for extremely large values of n and p. As these time trials reach out to the practical limits of current technology, we see that this is not critical -on average, even the longest trials here took less than 200 seconds. The weighted fused group Lasso results are perhaps more interesting, as it is harder to predict in advance the practical time performance of the algorithm. Surprisingly, when increasing n (p and k fixed) or increasing p (n and k fixed), the group fused Lasso eventually becomes as fast the iterative, deterministic group fused LARS. This suggests that at the limits of current technology, if k is small (say, less than 10), the potentially superior performance of the Lasso version (see later) may not even be punished by a slower run-time with respect to the LARS version. We suggest that this may be due to the Lasso optimization problem becoming relatively "easier" to solve when n or p increases, as we observed that the Lasso algorithm converged quickly to its final set of change-points. The main difference between the Lasso and LARS performance appears when the number of change-points increases: with respective empirical complexities cubic and linear in k, as predicted by the theoretical analysis, Lasso is already 1,000 times slower than LARS when we seek 100 change-points.

Accuracy for detection of a single change-point

Next, we tested empirically the accuracy the group fused Lasso for detecting a single change-point. We first generated multidimensional profiles of dimension p, with a single jump of height 1 at a position u, for different values of p and u. We added to the signals an i.i.d. Gaussian noise with variance σ2 α = 10.78, the critical value corresponding to α = 0.8 in Theorem 2. We ran 1000 trials for each value of u and p, and recorded how often the group fused Lasso with or without weights correctly identified the change-point. According to Theorem 2, we expect that, for the unweighted group fused Lasso, for 50 ≤ u < 80 there is convergence in accuracy to 1 when p increases, and for u > 80, convergence in accuracy to zero. This is indeed what is seen in Figure 3 (left panel), with u = 80 the limit case between the two different modes of convergence. The center panel of Figure 3 shows that when the default weights (5) are added, convergence in accuracy to 1 occurs across all u, as predicted by Theorem 3. In addition, the right-hand-side panel of Figure 3 shows results for the same trials except that change-point locations can vary uniformly in the interval u ± 2. We see that, as predicted by Theorem 4, the accuracy of the weighted group fused Lasso remains robust against fluctuations in the exact change-point location.

Accuracy for detecting multiple change-points

To investigate the potential for extending the results to the case of many shared change-points, we further simulated profiles of length n = 100 with a change-point at all of positions 10, 20, . . . , 90. We consider dimensions p between 1 and 500. Jumps at each change-point of each profile were drawn from a Gaussian with mean 0 and variance 1; we then added centered Gaussian noise with σ 2 ∈ {0.05, 0.2, 1} to each position in each profile. For each value of p and σ 2 , we ran one hundred trials of both implementations, with or without weights, and recorded the accuracy of each method, defined as the percentage of trials where the first 9 change-points detected by the method are exactly the 9 true change-points. Results are presented in Figure 4 (from left to right, resp. σ 2 = 0.05, 0.2, 1). Clearly, the group fused Lasso outperforms the group fused LARS, and the weighted version of each algorithm outperforms the unweighted version. Although the group LARS is usually considered a reliable alternative to the exact group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], this experiment shows that the exact optimization by block coordinate descent may be worth the computational burden if one is interested in accurate group selection. It also demonstrates that, as we conjectured in Section 5.3, the group fused Lasso can consistently estimate multiple change-points as the number of profiles increases.

Application to gain and loss detection

We now consider a possible application of our method for the detection of regions with frequent gains (positive values) and losses (negative values) among a set of DNA copy number profiles, measured by array comparative genomic hybridization (aCGH) technology [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF]. We propose a two-step strategy for this purpose: first, find an adequate joint segmentation of the signals; then, check the presence of gain or loss on each interval of the segmentation by summarizing each profile by its average value on the interval. Note that we do not assume that all profiles share exactly the same change-points, but merely see the joint segmentation as an adaptive way to reduce the dimension and remove noise from data.

In practice, we used group fused LARS on each chromosome to identify a set of 100 candidate change-points, and selected a subset of them by post-processing as described in Section 5.4. Then, in each piecewise-constant interval between successive shared change-points, we calculate the mean of the positive segments (shown in green in Figures 5(a) and 6(c)) and the mean of the negative segments (shown in red). The larger the mean of the positive segments, the more likely we are to believe that a region harbors an important common gain; the reasoning is analogous for important common losses and the mean of the negative segments. Obviously, many other statistical tests could be carried out to detect frequent gains and losses on each segment, once the joint segmentation is performed.

We compare this method for detecting regions of gain and loss with the state-of-the-art H-HMM method [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF], which has been shown to outperform several other methods in this setting. As [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF] have provided their algorithm online with several of their data sets tested in their article, we implemented our method and theirs (H-HMM) on their benchmark data sets.

In the first data set in [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF], the goal is to recover two regions -one amplified, one deleted, that are shared in 8 short profiles, though only 6 of the profiles exhibit each of the amplified or deleted regions. Performance is measured by area under ROC curve (AUC), following [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF]. Running H-HMM with the default parameters, we obtained an AUC (averaged over 10 trials) of 0.96 ± .01, taking on average 60.20 seconds. The weighted group fused LARS, asked to select 100 breakpoints and followed by dynamic programming, took 0.06 seconds and had an AUC of 0.97. Thus, the performance of both methods was similar, though weighted group fused LARS was around 1000 times faster.

The second data set was a cohort of lung cancer cell lines originally published in [START_REF] Coe | Differential disruption of cell cycle pathways in small cell and non-small cell lung cancer[END_REF][START_REF] Garnis | High resolution analysis of non-small cell lung cancer cell lines by whole genome tiling path array CGH[END_REF]. As in [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF], we concentrated on the 18 NSCLC adenocarcinoma (NA) cell lines. Figure 5 shows the score statistics obtained on Chromosome 8 when using either weighted group fused LARS or H-HMM. Weighted group fused LARS first selected 100 candidate change-points per chromosome, then followed optimization of the number of change-points by dynamic programming, took in total 4.7 seconds and finally selected 260 change-points. In contrast, H-HMM took 38 minutes (100 iterations, as given in the code provided by the authors). The H-HMM scores should look like those shown in Figure 4 (top panel) of [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF]; the difference is either due to the stochastic nature of the algorithm or using a different number of iterations than given in the sample code by the authors. In any case, at the MYC locus (near 13 × 10 7 bp), both methods strongly suggest a common gained region. However, the supposed advantage of H-HMM to very sparsely predict common gains and losses is not clear here; for example, it gives high common gain confidence to several fairly large genomic regions between 9 and 14 ×10 7 bp.

Application to bladder tumor aCGH profiles

We further considered a publicly available aCGH data set of 57 bladder tumor samples [START_REF] Stransky | Regional copy numberindependent deregulation of transcription in cancer[END_REF]. Each aCGH profile gave the relative quantity of DNA for 2215 probes. We removed the probes corresponding to sex chromosomes because the sex mismatch between some patients and the reference made the computation of copy number less reliable, giving us a final list of 2143 probes. Results are shown in Figure 6. 97 change-points were selected by the weighted group fused LARS; this took 1.1 seconds (Figure 6(c)). The H-HMM method (Figure 6(d)) took 13 minutes for 200 iterations (after 100 iterations convergence had not occured). We used the comprehensive catalogue of common -1 -0.5 0 0.5 shows the result of superimposing the smoothed versions of the 57 bladder tumor aCGH profiles obtained using weighted group fused LARS followed by dimension-selection. 6(c) shows the result of transforming the set of smoothed outputs into "scores" for amplification/deletion (see Section 6.4) and 6(d) the corresponding output for the H-HMM method [START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF]. Vertical black lines indicate chromosome boundaries. genomic alterations in bladder cancer provided in Table 2 in [START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF] to validate the method and compare with H-HMM. Our method (Figure 6(c)) concurred with the known frequently-amplified chromosome arms 20q, 8q, 19q, 1q, 20p, 17q, 19p, 5p, 2p, 10p, 3q and 7p, and frequently-lost 9p, 9q, 11p, 10q, 13q, 8p, 17p, 18q, 2q, 5q, 18p, 14q and 16q. The only known commonly-lost region which showed unconvincing common loss here was 6q. As for the H-HMM method (Figure 6(d)), it selects a small number of very small regions of gain and loss, which are difficult to verify with respect to the wellknown frequently amplified arms in [START_REF] Blaveri | Bladder cancer stage and outcome by array-based comparative genomic hybridization[END_REF]. As is suggested, the method may therefore be useful for selecting the precise location of important genes. However, as can be seen in Figure 6(a)-(b), many, but not all, alterations are much larger than those found with H-HMM, and where for example there are clearly several localized gains and losses in chromosome 8, H-HMM finds nothing at all. Perhaps the complexity of rearrangements in chromosome 8 is not easily taken into account by the H-HMM algorithm. Note finally that the weighted group fused LARS was more than 700 times faster than H-HMM.

Conclusion

We have proposed a framework that extends total-variation based approximation to the multidimensional setting, and developed two algorithms to solve the resulting convex optimization problem either exactly or approximately. We have shown theoretically and empirically that the group fused Lasso can consistently estimate a single change-points, and observed experimentally that this property is likely to hold also when several change-points are present. In particular, we observed both theoretically and empirically that increasing the number of profiles is highly beneficial to detect approximatively shared change-points, an encouraging property for biological applications where the accumulation of data measured on cohorts of patients promises to help in the detection of common genomic alterations.

Although we do not assume that all profiles have the same change-points, we estimate only shared change-points. In other words, we try to estimate the union of the change-points present in the profiles. This can be useful by itself, eg, for dimension reduction. If we wanted to detect change-points of individual profiles, we may either post-process the results of the group fused Lasso, or modify the formulation by, e.g., adding a TV penalty to each profile in addition to the group lasso penalty. Similarly, for some applications, we may want to add a ℓ 1 /ℓ 2 norm to the group fused Lasso objective function in order to constrain some or all signals to be frequently null. Finally, from a computational point of view, we have proposed efficient algorithms to solve an optimization problem (4) which is the proximal operator of more general optimization problems where a smooth convex functional of U is minimized with a constraint on the multidimensional TV penalty; this paves the way to the efficient minimization of such functionals using, e.g., accelerated gradient methods [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

Annex A: Computational lemmas

In this Annex we collect a few results useful to carry out the fast implementations claimed in Section 4. Remember that the n × (n -1) matrix X defined in (6) is defined by X i,j = d j for i > j, 0 otherwise. Since the design matrix X of the group Lasso problem (8) is obtained by centering each column of X to zero mean, its columns are given by: ∀i

= 1, . . . , n -1, X•,i =     i n -1 d i , . . . , i n -1 d i i , i n d i , . . . , i n d i n-i     ⊤ . (14)
We first show how to compute efficiently X⊤ R for any matrix R: Lemma 5. For any R ∈ R n×p , we can compute C = X⊤ R in O(np) operations and memory as follows:

1. Compute the n × p matrix r of cumulative sums r i,• = i j=1 R j• by the induction:

• r 1,• = R 1,• . • For i = 2, . . . , n, r i,• = r i-1,• + R i,• . 2. For i = 1, . . . , n -1, compute C i,• = d i (ir n,• /n -r i,• ) .
Proof. Let us denote D the (n -1) × (n -1) diagonal matrix with entries D i,i = d i . By Lemma 6, we know that X⊤ X = DV D, with V i,j = min(i, j) [nmax(i, j)] /n, for 1 ≤ i, j ≤ n -1. Since step 1 computes R = DR and step 5 computes C = DU , we just need to show that the U computed in step 4 satisfies U = V R to conclude that C = DV R = DV DR = X⊤ XR. By step 4, U is defined by the relation U i,• -U i-1,• = S -T i,• for i = 1, . . . , n -1 (with the convention U 0,• = 0), therefore we just need to show that (V i,• -V i-1,• ) R = S -T i,• for i = 1, . . . , n -1 to conclude. For 0 ≤ j < i ≤ n -1, we note that V i,j = j(n-i)/n (with the convention V 0,• = 0) and V i-1,j = j(n-i+1)/n, and therefore V i,j -V i-1,j = -j/n. For 1 ≤ i ≤ j ≤ n -1, we have V i,j = i(nj)/n and V i-1,j = (i -1)(nj)/n and therefore V i,j -V i-1,j = 1j/n. Combining these expressions we get, for i = 1, . . . , n -1:

(V i,• -V i-1,• ) R = - n-1 j=1 j Rj,• n + n-1 j=i Rj,• = S -T i,• ,
where S and T are defined in steps 2 and 3. This concludes the proof that C = X⊤ XR.

Next we show that X⊤ X -1 has a tridiagonal structure, resulting in fast matrix multiplication.

∆ i = d -1 a i+1 R i+1,• -d -1 a i R i,• a i+1 -a i .
2. Compute the successive rows of C according to:

C 1,• = d -1 a 1 R 1,• a 1 -∆ 1 , C i,• = d -1 a i (∆ i-1 -∆ i ) for i = 2, . . . , |A| -1 , C |A|,• = d -1 a |A| ∆ |A|-1 + R |A|,• n -a |A| . (16)
Proof. Let us denote V = X⊤ •,A X•,A . By Lemma 6 we know that, for 1 ≤ i ≤ j ≤ |A|, V i,j = d a i d a j a i (na j ) n .

V being symmetric semi-separable, one can easily check that V is invertible and admits as inverse a tridiagonal matrix with the following entries [START_REF] Baranger | Matrices tridiagonales symétriques et matrices factorisables[END_REF]:

V -1 i,i = d -2 a i 1 a i -a i-1 + 1 a i+1 -a i for i = 1, . . . , |A|, V -1 i,i+1 = V -1 i+1,i = - d -1 a i d -1 a i+1 a i+1 -a i for i = 1, . . . , |A| -1, (17)
where by convention we define a 0 = 0 and a |A|+1 = n. This tri-diagonal structure allows successive rows of C to be expressed as a sum of just a few terms. More precisely, for 1 < i < |A|, we obtain:

C i,• = - d -1 a i-1 d -1 a i R i-1,• a i -a i-1 + d -2 a i R i,• 1 a i -a i-1 + 1 a i+1 -a i - d -1 a i d -1 a i+1 R i+1,• a i+1 -a i = d -1 a i d -1 a i R i,• -d -1 a i-1 R i-1,• a i -a i-1 + d -1 a i R i,• -d -1 a i+1 R i+1,• a i+1 -a i = d -1 a i (∆ i-1 -∆ i ) .
Similarly, for i = 1 and i = |A| we easily recover [START_REF] Tibshirani | Sparsity and smoothness via the fused lasso[END_REF].

Annex B: Proof of Lemma 1

The solution of ( 4) is constant, i.e., corresponds to β = 0 (no change-point), as long as the KKT conditions [START_REF] Yao | Estimating the number of change-points via schwarz criterion[END_REF] are satisfied for β = 0. This translates to X⊤

On the other hand, by [START_REF] Lavielle | Adaptive detection of multiple change-points in asset price volatility[END_REF] we have for any i ∈ [1, n -1],

X⊤ W i,• = d i   i j=1 i n -1 W j,• + n j=i+1 i n W j,•   . Since E W ⊤ i,• W j,• = δ i,j σ 2 I p ,
where δ i,j is the Dirac function, we have for 1 ≤ i ≤ j ≤ n -1:

E X⊤ W ⊤ i,• X⊤ W j,• = d i d j i i n -1 j n -1 + (j -i) i n j n -1 + (n -j) i n j n σ 2 I p = d i d j i (n -j) n σ 2 I p . (19)
In summary, we have shown that ĉ is jointly Gaussian with E (ĉ i,• ) given by [START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF] and covariance between ĉi,• and ĉj,• given by [START_REF] Friedman | Pathwise coordinate optimization[END_REF].

In particular, if we denote F i = ĉi,• 2 , then, for i ≤ u, F i n/ d 2 i i(ni)σ 2 follows a non-central χ 2 distribution with p degrees of freedom and non-centrality parameter p β2 p d 2 u i(nu) 2 / n(ni)σ 2 . In particular,

EF i = p βp 2 d 2 i d 2 u i 2 (n -u) 2 n 2 + pd 2 i i(n -i) n σ 2 ,
where

σ2 m = n β2 E(1 -α) 2 (α m -1 2n ) α m -1 2 -1 2n
.

This shows that G a > G a-1 if and only if σ 2 < σ2 m . Since b > n/2, we also know that û ≤ b, i.e., û ∈ [a, b] in that case. The case 1 ≤ a ≤ b < n/2 can be treated similarly. To conclude the proof it suffices to observe that E(1α) 2 = (1 -Eα) 2 + var(α) .

• In the weighted case d i = n i(n-i) for i = 1, . . . , n -1, we obtain from ( 26) and ( 27) that for i ≤ a:

G i = β2 i n -i E n -U U + σ 2 . ( 29
)
This is always an increasing function of i on [1, a], showing that the maximum of G i can not be strictly smaller than a. By symmetry, it can also never be larger than b, from which we conclude that it is always between a and b, i.e., in the support of P U .

Figure 1 :

 1 Figure 1: First 5 change-points detected on three simulated profiles by TV denoising of each profile (left) and by joint TV denoising (right).

Algorithm 1

 1 Block coordinate descent algorithmRequire: centered data Ȳ , regularization parameter λ.1: Initialize A ← ∅, β = 0, C ← X⊤ Ȳ .

Algorithm 2

 2 end loop to efficiently run the fast group LARS method. Computing X⊤ Ȳ in line 1 can be done in O(np) using Lemma 5. To compute the descent direction (line 6), we first compute w in O(|A|p) using Lemma 8, then a in O(np) using Lemma 7. To find the descent step (line 7), we need to solve n polynomial equations of degree 2, the coefficients of which are computed in O(p), resulting in a O(np) complexity. Overall the main loop for each new change-point (lines 2-10) takes O(np) in computation and memory, resulting in O(npk) complexity in time and O(np) in memory to find the first k change-points. We provide in Section 6 empirical results that confirm this theoretical complexity. Group fused LARS algorithm Require: centered data Ȳ , number of breakpoints k. 1: Initialize A ← ∅, ĉ ← X⊤ Ȳ . 2: for i = 1 to k do 3: if i=1 then 4:

Lemma 1 .

 1 Assume, without loss of generality, that u ≥ n/2, and let, for i ∈ [1, n -1],

Theorem 4 . 1 .

 41 Let α = U/n be the random position of the change-point on [0, 1] and α m = a/n and α M = b/n the position of the left and right boundaries of the support of P U scaled to [0, 1]. If 1/2 ∈ (α m , α M ), then for any noise level σ 2 , the probability that the first change-point selected by the unweighted group fused Lasso (3) is in the support of P U tends to 1 as p

Figure 2 :

 2 Figure 2: Speed trials for group fused LARS (top row) and Lasso (bottom row). Left column: varying n, with fixed p = 10 and k = 10; center column: varying p, with fixed n = 1000 and k = 10; right column: varying k, with fixed n = 1000 and p = 10. Figure axes are log-log. Results are averaged over 100 trials.

Figure 3 :

 3 Figure 3: Single change-point accuracy for the group fused Lasso. Accuracy as a function of the number of profiles p when the change-point is placed in a variety of positions u = 50 to u = 90 (left and centre plots, resp. unweighted and weighted group fused Lasso), or: u = 50 ± 2 to u = 90 ± 2 (right plot, weighted with varying change-point location), for a signal of length 100.

Figure 4 :

 4 Figure 4: Multiple change-point accuracy. Accuracy as a function of the number of profiles p when change-points are placed at the nine positions {10, 20, . . . , 90} and the variance σ 2 of the centered Gaussian noise is either 0.05 (left), 0.2 (center) and 1 (right). The profile length is 100.

Figure 5 :

 5 Figure5: Joint scores for a set of 18 NSCLC adenocarcinoma cell lines. 5(a) using weighted group fused LARS; 5(b) using H-HMM with the actual code provided by[START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF].

Figure 6 :

 6 Figure6: Bladder cancer profiles.6(a) shows one of the original 57 profiles and its associated smoothed version. 6(b) shows the result of superimposing the smoothed versions of the 57 bladder tumor aCGH profiles obtained using weighted group fused LARS followed by dimension-selection. 6(c) shows the result of transforming the set of smoothed outputs into "scores" for amplification/deletion (see Section 6.4) and 6(d) the corresponding output for the H-HMM method[START_REF] Shah | Modeling recurrent DNA copy number alterations in array CGH data[END_REF]. Vertical black lines indicate chromosome boundaries.

Lemma 8 . 1 R

 81 For any set A = a 1 , . . . , a |A| of distinct indices with 1 ≤ a 1 < . . . < a |A| ≤ n -1, the matrix X⊤ •,A X•,A is invertible, and for any |A| × p matrix R, the matrix C = X⊤ •,A X•,A can be computed in O(|A|p) in time and memory by 1. For i = 1, . . . , |A| -1, compute

Ȳ

  ≤ λ for all i. The first change-point occurs when λ = max i X⊤•,iȲ , and the change-point is precisely located in the position i that reaches the maximum. Therefore the first change-point is the row with the largest Euclidean norm of the matrix:ĉ = X⊤ Ȳ = X⊤ Xβ * + X⊤ W .The entries of the matrix ĉ are therefore jointly Gaussian. Since only the u-th row β u,• of β is non-zero, we getE(ĉ) = X⊤ Xβ * = X⊤ X•,u β * u,• . Using Lemma 6 we compute:E(ĉ i,• ) = X⊤ Xβ * i,• = d i d u i(n-u) n β * u,• for 1 ≤ i ≤ u , d i d u u(n-i) n β * u,• for u ≤ i ≤ n -1 .

Available at http://cbio.ensmp.fr/GFLseg

Proof. Using [START_REF] Lavielle | Adaptive detection of multiple change-points in asset price volatility[END_REF] we obtain the i-th row of C = X⊤ R, for i = 1, . . . , n -1, as follows:

Next, we show how to compute efficiently submatrices of the (n -1) × (n -1) matrix X⊤ X.

, denoting u = min(a i , b j ) and v = max(a i , b j ), we easily get from ( 14) an explicit formula for V i,j , namely,

The next lemma provides another useful computational trick to compute efficiently X⊤ XR for any matrix R:

j=i Rj,• by the induction:

) by the induction:

Each step in Lemma 7 has complexity O(np) in memory and time, leading to an overall complexity in O(np) to compute X⊤ XR. We note that if R is row-sparse, i.e., is several rows of R are null, then the first two steps have complexity O(sp), where s is the number of non-zero rows in R. Although this does not change the overall complexity to compute X⊤ XR, this leads to a significant speed-up in practice when s ≪ n. and since lim p→+∞ βp 2 = β2 , we get that F i /p converges in probability to

A similar computation shows that for i ≥ u, F i /p converges in probability to

Note that ( 20) and ( 21) are equivalently defined in [START_REF] Birgé | Gaussian model selection[END_REF]. Now, let V = argmax i∈[1,n-1] G i . For any v ∈ V and j / ∈ V , the probability of the event F v > F j tends to 1, because G v > G j . By the union bound the probability of the event max j / ∈V F i < max v∈V F v also converges to 1, showing that the probability to select a change-point in V converges to 1 as p → +∞.

Annex C: Proof of Theorem 2

By Lemma 1, we know that the first change-point selected by ( 4) is in argmax i∈ [1,n] G i with probability tending to 1 as p increases, where G i is defined in [START_REF] Birgé | Gaussian model selection[END_REF]. We will therefore asymptotically select the correct change-point u if and only if G u = max i∈[1,n-1] G i . Remember we assume, without lack of generality, that u ≥ n/2. For u ≤ i ≤ n -1, we observe that G i given by ( 21) is a decreasing function of i as a sum of two decreasing functions. Therefore, it always holds that G u = max i∈[u,n-1] G i , and we just need to check whether or not G u = max i∈[1,u] G i holds.

For i ∈ [1, u], G i given by ( 20) is a second-order polynomial of i, which is equal to 0 at i = 0 and strictly positive for i

where α = u/n and

.

This shows that, when α > 1/2 + 1/(2n), G u > G u-1 if and only if σ < σ. On the other hand, when α = 1/2 or 1/2 + 1/(2n), we have always that G u > G u-1 .

Annex D: Proof of Theorem 3

As for the proof of Theorem 2, we need to check whether or not G u = max i∈[1,n-1] G i , where G i is defined in [START_REF] Birgé | Gaussian model selection[END_REF], to deduce whether the method selects the correct change-point u or a different position with probability tending to 1 when p increases. Substituting weights d i defined in [START_REF] Talih | Structural learning with time-varying components: tracking the crosssection of financial time series[END_REF] into G i , we obtain:

It is then easy to see that ( 23) is increasing on [1, u], and decreasing on [u, n -1], showing that we always have argmax i∈[1,n-1] G i = u. The result then follows from Lemma 1.

Annex E: Proof of Theorem 4

Following the proof of Lemma 1, let us estimate

we first observe by [START_REF] Rinaldo | Properties and refinements of the fused lasso[END_REF] that

Therefore,

and by independence of β i and U i :

Since (β i , U i ) i=1,...,p are independent of the noise, we obtain that F i /p converges in probability to

As in Lemma 1 we can conclude that the method will select the position

G i with probability tending to 1 as p increases.

Let us now assume that the support of P U is an interval [a, b] (corresponding to a possible range of fluctuation of a change-point). Then, we observe that for i ≤ a, G i in (26) reduces to

Let us now consider the two possible weighting schemes.

• In the unweighted case d i = 1 for i = 1, . . . , n -1, we obtain from ( 27) that for i ≤ a:

While the first term in ( 28) is strictly increasing on [0, a], the second term moves the maximum of G i towards n/2. This shows that the maximum of G i is always at least a when a ≤ n/2. By symmetry, it is also always smaller or equal to b when b ≥ n/2. When n/2 ∈ [a, b], we deduce that for any σ 2 > 0, û ∈ [a, b]. Otherwise, let us suppose without lack of generality that n/2 < a ≤ b. Then, G i being quadratic on [0, a] and equal to 0 at 0, the maximum of G i will not occur before a if and only if G a-1 < G a . A computation similar to the one in the proof of Theorem 2 shows that G a -G a-1 = 2 σ2 mσ 2 α m -