Vincent Brault
email: vincentbrault@agroparistech.fr

Julien Chiquet

Céline Lévy-Leduc

Fast Detection of Block Boundaries in Block-Wise Constant Matrices

Keywords: Change-points, High-dimensional sparse linear model, HiC experiments

We propose a novel approach for estimating the location of block boundaries (change-points) in a random matrix consisting of a block wise constant matrix observed in white noise. Our method consists in rephrasing this task as a variable selection issue. We use a penalized least-squares criterion with an 1-type penalty for dealing with this problem. We first provide some theoretical results ensuring the consistency of our change-point estimators. Then, we explain how to implement our method in a very efficient way. Finally, we provide some empirical evidence to support our claims and apply our approach to data coming from molecular biology which can be used for better understanding the structure of the chromatin.

Introduction

Detecting automatically the block boundaries in a block wise constant matrix corrupted with noise is a very important issue which may have several applications. One of the main situations in which this problem occurs is the detection of chromosomal regions having close spatial location in the nucleus. Detecting such regions will improve our understanding of the influence of the chromosomal conformation on the cells functioning. The data provided by the most recent technology called HiC consist of a list of pairs of locations along the chromosome which are often summarized as a square matrix such that each entry corresponds to the number of interactions between two positions along the chromosome, see [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]. Since this matrix can be modeled as a block wise matrix corrupted by some additional noise, it is of particular interest to design an efficient and fully automated method to find the block boundaries of large matrices, which may typically have several thousands of rows and columns, in order to identify the interacting chromosomal regions.

A large literature is dedicated to the change-point detection issue for onedimensional data. This problem can be addressed from a sequential (online) [START_REF] Tartakovsky | Sequential Analysis: Hypothesis Testing and Changepoint Detection[END_REF] or from a retrospective (off-line) [START_REF] Brodsky | Non-parametric statistical diagnosis: problems and methods[END_REF] point of view. Many off-line approaches are based on the dynamic programming algorithm which retrieves K change-points within n observations of a one-dimensional signal with a complexity of O(Kn 2 ) in time [START_REF] Kay | Fundamentals of statistical signal processing: detection theory[END_REF]. Such a complexity is however prohibitive for dealing with very large data sets. In this situation, [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF] proposed to rephrase the change-point estimation issue as a variable selection problem. This approach has also been extended by [START_REF] Vert | Fast detection of multiple change-points shared by many signals using group lars[END_REF] to find shared change-points between several signals. To the best of our knowledge no method has been proposed for addressing the case of two-dimensional data where the number of rows and columns may be very large (n × n ≈ 5000 × 5000, namely 2 × 10 7 observations). The only statistical approach proposed for retrieving the change-point positions in the two-dimensional framework is the one devised by [START_REF] Lévy-Leduc | Two-dimensional segmentation for analyzing hi-c data[END_REF] but it is limited to the case where the blockwise matrix is assumed to be blockwise constant on the diagonal and constant outside the diagonal blocks.

It has first to be noticed that the classical dynamic programming algorithm cannot be applied in such a framework since the Markov property does not hold anymore. Moreover, the group-lars approach of [START_REF] Vert | Fast detection of multiple change-points shared by many signals using group lars[END_REF] cannot be used in this framework since it would only provide change-points in columns and not in rows. As for the generalized Lasso recently devised by [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] or the two dimensional fused Lasso of [START_REF] Hoefling | A path algorithm for the fused lasso signal approximator[END_REF], they are very helpful for image denoising but do not give access to the change-point positions since they are not derived to provide a partitioning of a matrix in rectangular blocks.

The paper is organized as follows. In Section 2, we first describe how to rephrase the problem of two-dimensional change-point estimation as a high dimensional sparse linear model and give some theoretical results which prove the consistency of our change-point estimators. In Section 3, we describe how to efficiently implement our method. In Section 4, we provide experimental evidence of the relevance of our approach on synthetic and real data coming from molecular biology.

Statistical Framework

Statistical Modeling

In this section, we explain how the two-dimensional retrospective change-point estimation issue can be seen as a variable selection problem. Our goal is to estimate

t 1 = (t 1,1 , . . . , t 1,K 1 ) and t 2 = (t 2,1 , . . . , t 2,K 2 ) from the random matrix Y = (Y i,j ) 1≤i,j≤n defined by Y = U + E, (1)
where U = (U i,j ) is a blockwise constant matrix such that

U i,j = μ k, if t 1,k-1 ≤ i ≤ t 1,k -1 and t 2, -1 ≤ j ≤ t 2, -1,
with the convention t 1,0 = t 2,0 = 1 and t 1,K 1 +1 = t 2,K 2 +1 = n + 1. An example of such a matrix U is displayed in Figure 1 (left). The entries E i,j of the matrix E = (E i,j ) 1≤i,j≤n are iid zero-mean random variables. With such a definition the Y i,j are assumed to be independent random variables with a blockwise constant mean.

Let T be a n × n lower triangular matrix with nonzero elements equal to one and B a sparse matrix containing null entries except for the B i,j such that (i, j) ∈ {t 1,0 , . . . , t 1,K 1 } × {t 2,0 , . . . , t 2,K 2 }. Then, (1) can be rewritten as follows:

Y = TBT + E, (2)
where T denotes the transpose of the matrix T. For an example of a matrix B, see Figure 1 (right). Let Vec(X) denotes the vectorization of the matrix X formed by stacking the columns of X into a single column vector then Vec(Y) = Vec(TBT ) + Vec(E). Hence, by using that Vec(AXC) = (C ⊗ A)Vec(X), where ⊗ denotes the Kronecker product, (2) can be rewritten as:

Y = X B + E, (3)
where

Y = Vec(Y), X = T ⊗ T, B = Vec(B)
B(λ n ) = Argmin B∈R n 2 Y -X B 2 2 + λ n B 1 , ( 4
)
where u 2 2 and u 1 are defined for a vector 4) is related to the popular Least Absolute Shrinkage and Selection Operator (LASSO) in least-square regression. Thanks to the sparsity enforcing property of the 1 -norm, the estimator B of B is expected to be sparse and to have non-zero elements matching with those of B. Hence, retrieving the positions of the non zero elements of B thus provides estimators of (t 1,k ) 1≤k≤K 1 and of (t 2,k ) 1≤k≤K 2 . More precisely, let us define by A(λ n ) the set of active variables:

u in R N by u 2 2 = N i=1 u 2 i and u 1 = N i=1 |u i |. Criterion (
A(λ n ) = j ∈ {1, . . . , n 2 } : B j (λ n ) = 0 .
For each j in A(λ n ), consider the Euclidean division of (j -1) by n, namely (j -1) = nq j + r j then

t 1 = ( t 1,k ) 1≤k≤| A1(λn)| ∈ {r j + 1 : j ∈ A(λ n )}, t 2 = ( t 2, ) 1≤ ≤| A2(λn)| ∈ {q j + 1 : j ∈ A(λ n )} where t 1,1 < t 1,2 < • • • < t 1,| A1(λn)| , t 2,1 < t 2,2 < • • • < t 2,| A2(λn)| . (5) In (5), | A 1 (λ n )| and | A 2 (λ n )| correspond
to the number of distinct elements in {r j : j ∈ A(λ n )} and {q j : j ∈ A(λ n )}, respectively.

As far as we know, neither thorough practical implementation nor theoretical grounding have been given so far to support such an approach for change-point estimation in the two-dimensional case. In the following section, we give theoretical results supporting the use of such an approach.

Theoretical Results

In order to establish the consistency of the estimators t 1 and t 2 defined in (5), we shall use assumptions (A1-A4). These assumptions involve the two following quantities

I min = min 0≤k≤K 1 |t 1,k+1 -t 1,k | ∧ min 0≤k≤K 2 |t 2,k+1 -t 2,k |, J min = min 1≤k≤K 1 ,1≤ ≤K 2 +1 |μ k+1, -μ k, | ∧ min 1≤k≤K 1 +1,1≤ ≤K 2 |μ k, +1 -μ k, |,
which corresponds to the smallest length between two consecutive change-points and to the smallest jump size between two consecutive blocks, respectively.

(A1) The random variables (E i,j ) 1≤i,j≤n are iid zero mean random variables such that there exists a positive constant β such that for all ν in R, E[exp(νE 1,1 )] ≤ exp(βν 2 ).

(A2) The sequence (λ n ) appearing in ( 4) is such that (nδ n J min ) -1 λ n → 0, as n tends to infinity.

(A3) The sequence (δ n ) is a non increasing and positive sequence tending to zero such that nδ n J min 2 / log(n) → ∞, as n tends to infinity.

(A4) I min ≥ nδ n .

Proposition 1. Let (Y i,j ) 1≤i,j≤n be defined by ( 1) and t 1,k , t 2,k be defined by [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF]. Assume that

| A 1 (λ n )| = K 1 and that | A 2 (λ n )| = K 2 ,
with probabilty tending to one, then, P max

1≤k≤K 1 t 1,k -t 1,k ≤ nδ n ∩ max 1≤k≤K 2 t 2,k -t 2,k ≤ nδ n -→ n→∞ 1. (6)
The proof of Proposition 1 is based on the two following lemmas. The first one comes from the Karush-Kuhn-Tucker conditions of the optimization problem stated in ( 4). The second one allows us to control the supremum of the empirical mean of the noise. Lemma 1. Let (Y i,j ) 1≤i,j≤n be defined by [START_REF] Breiman | Classification and Regression Trees[END_REF]. Then, U = X B, where X and B are defined in ( 3) and ( 4) respectively, is such that

n k=rj +1 n =qj +1 Y k, - n k=rj +1 n =qj +1 U k, = λ n 2 sign( B j ), if B j = 0, (7)
n k=rj +1 n =qj +1 Y k, - n k=rj +1 n =qj +1 U k, ≤ λ n 2 , if B j = 0, (8)
where q j and r j are the quotient and the remainder of the euclidean division of (j -1) by n, respectively, that is (j -1) = nq j + r j . In [START_REF] Kay | Fundamentals of statistical signal processing: detection theory[END_REF], sign denotes the function which is defined by sign(x) = 1, if x > 0, -1, if x < 0 and 0 if x = 0. Moreover, the matrix U, which is such that U = Vec( U), is blockwise constant and satisfies [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF].

U i,j = μ k, , if t 1,k-1 ≤ i ≤ t 1,k -1 and t 2, -1 ≤ j ≤ t 2, -1, k ∈ {1, . . . , | A 1 (λ n )|}, ∈ {1, . . . , | A 2 (λ n )|}, where the t 1,k , t 2,k , A 1 (λ n ) and A 2 (λ n ) are defined in
Lemma 2. Let (E i,j ) 1≤i,j≤n be random variables satisfying (A1). Let also

(v n ) and (x n ) be two positive sequences such that v n x 2 n / log(n) → ∞, then P ⎛ ⎝ max 1≤rn<sn ≤n |rn-sn|≥vn (s n -r n ) -1 sn-1 j=rn E n,j ≥ x n ⎞ ⎠ -→ n→∞ 0, the result remaining valid if E n,j is replaced by E j,n .
The proofs of Proposition 1, Lemmas 1 and 2 can be seen as a natural extension of the results of [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF].

Implementation

In order to identify a series of change-points we look for the whole path of solutions in (4), i.e., { B(λ), λ min < λ < λ max } such that | Â(λ max )| = 0 and | Â(λ min )| = s with s a predefined maximal number of activated variables. To this end it is natural to adopt the famous homotopy/LARS strategy [START_REF] Efron | Least angle regression[END_REF][START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF]. Such an algorithm identifies in Problem (4) the successive values of λ that correspond to the activation of a new variable, or the deletion of one that became irrelevant. However, the existing implementations do not apply here since the size of the design matrix X -even for reasonable n -is challenging both in terms of memory requirement and computational burden. To overcome these limitations, we need to take advantage of the particular structure of the problem. In the following lemmas, we show that the most involving computations in the LARS can be made extremely efficiently thanks to the particular structure of X . Lemma 3. For any vector v ∈ R n 2 , computing X v and X v requires at worse 2n 2 operations. Lemma 4. Let A = {a 1 , . . . , a K } and for each j in A let us consider the Euclidean division of j -1 by n given by j -1 = nq j + r j , then

X X A,A 1≤k, ≤K = ((n -(q a k ∨ q a )) × (n -(r a k ∨ r a ))) 1≤k, ≤K . ( 9
)
Moreover, for any non empty subset A of distinct indices in 1, . . . , n 2 , the matrix X A X A is invertible. Remark 1. We were able to obtain a closed-form expression of the inverse (X A X A ) -1 for some special cases of the subset A, namely, when the quotients/ratios associated with the Euclidean divisions of the elements of A are endowed with a particular ordering. For addressing any general problem though, we rather solve system involving X A X A by means of a Cholesky factorization which is updated along the homotopy algorithm. These updates correspond to adding or removing an element at a time in A and are performed efficiently as stated in Lemma 5.

These lemmas are the building blocks for our LARS implementation given in Algorithm 1, where we detail the leading complexity associated with each part. The global complexity is in O(mn 2 + ms 2 ) where m is the final number of steps in the while loop. These steps include all the successive additions and deletions needed to reach s, the final targeted number of active variables. At the end of day, we have m block-wise prediction Ŷ associated with the series of m estimations of B(λ). Concerning the memory requirements, we only need to store the n × n data matrix Y once. Indeed, since we have at our disposal the analytic form of any sub matrix extracted from X X , we never need to compute neither store this large n 2 × n 2 matrix. This paves the way for quickly processing data with thousands of rows and columns.

Numerical Experiments

Synthetic Data

The goal of this section is to assess the statistical and numerical performances of our methodology. We generated observations according to Model [START_REF] Breiman | Classification and Regression Trees[END_REF] where U is a symmetric blockwise constant matrix defined by

μ k, k∈{1,...K 1 +1}, ∈{1,...K 2 +1} = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 0 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ , ( 10
)
and the E i,j are zero mean i.i.d. Gaussian random variables of variance σ 2 where σ is in {1, 2, 5}. Some examples of data generated from this model with n = 500, K 1 = K 2 = 4 can be found in Figure 2 for the different samples (in gray) and the median of the mean square errors in thick as a function of the number of active variables s defined in Algorithm 1. We can see from this figure that even in the high noise level case, the mean square error is small. Moreover, the ROC curves displayed in the bottom part of Figure 2 ensure that the changepoints in rows: (t 1,k ) 1≤k≤K 1 are properly retrieved with a very small error rate even in high noise level frameworks. The same results hold for the change-points in columns: (t 2,k ) 1≤k≤K 2 but are not displayed in order to save space.

Since, to the best of our knowledge, no two-dimensional method are available, we propose to compare our approach to an adaptation of the CART procedure of [START_REF] Breiman | Classification and Regression Trees[END_REF] and to an adaptation of [START_REF] Harchaoui | Multiple change-point estimation with a total variation penalty[END_REF] (HL) dedicated to univariate observations. We adapt the CART methodology by using the successive boundaries provided by CART as change-points for the two-dimensional data. The associated ROC curve is displayed with '•' in Figure 3. For adapting the HL methodology, we apply it to each row of the data matrix and for each λ, we obtain the change-points of each row. The change-points appearing in the different rows are claimed to be change-points for the two-dimensional data either if they appear at least in one row (the associated ROC curve for this approach called HL1 is displayed with '+' in Figure 3) or if they appear in ([n/2] + 1) rows (the associated ROC curve for this approach called HL2 is displayed with ' ' in Figure 3). Since the procedures HL1 and HL2 are much slower than ours, the ROC curves are displayed for matrices of size 250 × 250. We can see from Figure 3 that our method outperforms the other ones.

Numerical Performances. We implemented Algorithm 1 in C++ using the library armadillo for linear algebra [START_REF] Sanderson | Armadillo: An open source C++ linear algebra library for fast prototyping and computationally intensive experiments[END_REF] and also provide an interface to the R platform [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] through the R package blockseg which is available from the Comprehensive R Archive Network (CRAN). All experiments were conducted on Linux workstation with Intel Xeon 2.4 GHz processor and 8 GB of memory.

We generated data as in Model [START_REF] Osborne | A new approach to variable selection in least squares problems[END_REF] for different values of n: n ∈ {100, 250, 500, 1000, 2500, 5000} and different values of the maximal number of activated variables: s ∈ {50, 100, 250, 500, 750}. The median runtimes obtained from 4 replications (+ 2 for warm-up) are reported in Figures 4. Left of Figure 4 (resp. right) gives the runtimes in seconds as a function of s (resp. of n). These results give experimental evidence for the theoretical complexity O(mn 2 + ms 2 ) that we established in Section 3 and thus for the computational efficiency of our approach: applying blockseg to matrices containing 10 7 entries takes less than 2 minutes.

Model Selection. In practice, we take s = K 2 max where K max is an upper bound for K 1 and K 2 . For choosing the final change-points we shall adapt the wellknown stability selection approach devised by [START_REF] Meinshausen | Stability selection[END_REF]. More precisely, we randomly choose M times n/2 columns and n/2 rows of the matrix Y and for each set of observations thus generated we select s = K 2 max active variables. Finally, after the M data resamplings, we keep the change-points which appear a number of times larger than a given threshold. By the definition of the change-points given in (5), a change-point t 1,k or t 2, may appear several times in a given set of resampled observations. Hence, the score associated with each change-point corresponds to the sum of the number of times it appears in each of the M resamplings. To evaluate the performances of this methodology, we generated observations according to the model defined in Section 4.1 with s = 225 and M = 100. The results are given in Figure 5 which displays the score associated to each changepoint for a given matrix Y (top). We can see from the top part of Figure 5 some spurious change-points appearing near from the true change-point positions. In order to identify the most representative change-point in a given neighborhood, we keep the one with the largest score among a set of contiguous candidates. The result of such a post-processing is displayed in the second and third rows of Figure 5. More precisely the boxplots associated to the estimation of K 1 (resp. the histograms of the estimated change-points in rows) are displayed in the middle (resp. bottom) part of Figure 5 for different threshold (resp. when the threshold is equal to T = 30% of the largest score). We can see from these figures that when thresh is in the interval [20,40] the number and the location of the change-points are very well estimated even in the high noise level case.

Application to HiC Data

In this section, we applied our methodology to publicly available data (http:// chromosome.sdsc.edu/mouse/hi-c/download.html) already studied by [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]. More precisely, we studied the interaction matrices of Chromosomes 1 and 19 of the mouse cortex at a resolution 40 kb and we compared the number and the location of the estimated change-points found with our approach with those obtained by [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF] on the same data since no ground truth is available. The matrices of these interaction matrices are displayed in Figure 6. We can see from this figure that modeling these matrices as block wise constant matrices corrupted with white seems to be relevant. We display in Figure 7 the number of change-points in rows found by our approach as a function of the threshold thresh used in our adaptation of the stability selection approach presented in the previous section. We also display in this figure a red line corresponding to the number of change-points found by [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF]. We also compute the two parts of the Hausdorff distance for the change-points in rows which is defined by

d t B , t = max d 1 t B , t , d 2 t B , t , (11)
where t and t B are the change-points in rows found by our approach and [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF], respectively. In [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF],

d 1 (a, b) = sup b∈b inf a∈a |a -b| , ( 12
)
d 2 (a, b) = d 1 (b, a) . ( 13
)
More precisely, Figure 8 displays the boxplots of the d 1 and d 2 parts of the Hausdorff distance without taking the supremum in orange and blue, respectively. We can observe from Figure 8 that some differences indeed exist between the segmentations produced by the two approaches but that the boundaries of the blocks are quite close when the number of estimated change-points are the same, which is the case when thresh = 1.8% (left) and 10% (right).

In the case where the number of estimated change-points are on a par with those of [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF], we can see from Figure 9 that the change-points found with our strategy present a lot of similarities with those found by the HMM based approach of [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF].

Our method also gives access to the estimated change-point positions for different values of the thresholds. However, contrary to our method, the approach of [START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF] can only deal with binned data at the resolution of several kilobases of nucleotides. The very low computational burden of our strategy paves the way for processing data collected at a very high resolution, namely at the nucleotide resolution, which is one of the main current challenges of molecular biology.

Conclusion

In this paper, we proposed a novel approach for retrieving the boundaries of a block wise constant matrix corrupted with noise by rephrasing this problem as a variable selection issue. Our approach is implemented in the R package blockseg which will be available from the Comprehensive R Archive Network (CRAN). In the course of this study, we have shown that our method has two main features which make it very attractive. Firstly, it is very efficient both from the theoretical and practical point of view. Secondly, its very low computational burden makes its use possible on very large data sets coming from molecular biology.

Fig. 1 .

 1 Fig. 1. Left: An example of a matrix U with n = 16, K 1 = 3 and K 2 = 4. Right: The matrix B associated to this matrix U.

Lemma 5 .

 5 Assume that we have at our disposal the Cholesky factorization of X A X A . The updated factorization on the extended set A ∪ {j} only requires solving an |A|-size triangular system, with complexity O(|A| 2 ). Moreover, the downdated factorization on the restricted set A\ {j} requires a rotation with negligible cost to preserve the triangular form of the Cholesky factorization after a column deletion.

Algorithm 1 . 2 );

 12 Fast LARS for two-dimensional change-point detection Input: data matrix Y, maximal number of active variables s. // Initialization Start with no change-point A ← ∅, B = 0; Compute current correlations ĉ = X Y with Lemma 3 ; // O(n 2 ) while λ > 0 or |A| < s do // Update the set of active variables Determine next change-point(s) by setting λ ← ĉ ∞ and A ← {j : ĉj = λ}; Update the Cholesky factorization of X A XA with Lemma 4; // O(|A| 2 ) // Compute the direction of descent Get the unormalized direction wA ← X •A X•A Normalize wA ← α wA with α ← 1/ w A sign(ĉA); Compute the equiangular vector uA = XAwA and a = X uA with Lemma 3; // O(n 2 ) // Compute the direction step Find the maximal step preserving equicorrelation γin ← min + Find the maximal step preserving the signs γout ← min + j∈A -BA/wA ; The direction step that preserves both is γ ← min(γin, γout); Update the correlations ĉ ← ĉγa and BA ← BA + γwA accordingly ; // O(n) // Drop variable crossing the zero line if γout < γin then Remove existing change-point(s) A ← A\ j ∈ A : Bj = 0 ; Downdate the Cholesky factorization of X A XA; // O(|A|) Output: Sequence of triplet (A, λ, B) recorded at each iteration.

Fig. 2 .

 2 Fig. 2. Top: Examples of 500 × 500 matrices Y generated from the model described in Section 4.1. Middle: Mean square errors n -2 B -B 2 2 for the different realizations in gray and the median of the mean square errors in thick line as a function of the number of nonzero elements in B for each scenario. Bottom: ROC curves for the estimated change-points in rows.

Fig. 3 .Fig. 4 .

 34 Fig.3. ROC curves for the estimated change-points in rows for our method (' '), HL1 ('+'), HL2 (' ') and CART ('•'). timings (seconds, log-scale)

Fig. 5 .

 5 Fig. 5. Top: Scores associated to each estimated change-points for different values of σ; the true change-point positions in rows and columns are located at 101, 201, 301 and 401. Middle: Boxplots of the estimation of K 1 for different values of σ and thresh after the post-processing step. The horizontal line corresponds to the true value of K 1 . Bottom: Histograms of the estimated change-points in rows for different values of σ after the post-processing step with thresh=30%.

Fig. 6 .

 6 Fig. 6. Raw interaction matrices of Chromosome 1 (left) and Chromosome 19 (right) for the mouse cortex. The darkest entries correspond to the lowest values.

Fig. 7 .

 7 Fig.7. Number of change-points in rows found by our approach as a function of the threshold ('•') in % for the interaction matrices of Chromosome 1 (left) and Chromosome 19 (right) of the mouse cortex. The horizontal line corresponds to the number of change-points found by[START_REF] Dixon | Topological domains in mammalian genomes identified by analysis of chromatin interactions[END_REF].

Fig. 8 .

 8 Fig. 8. Boxplots for the infimum parts of the Hausdorff distances d1 and d2 between the change-points found by [3] and our approach for the Chromosome 1 (left) and the Chromosome 19 (right) of the mouse cortex for the different thresholds in %. In each plot, the boxplot on the left corresponds to d1 and the boxplot on the right corresponds to d2.

Fig. 9 .

 9 Fig. 9. Topological domains detected by [3] (upper triangular part of the matrix) and by our method (lower triangular part of the matrix) from the interaction matrix of Chromosome 1 (left) and Chromosome 19 (right) of the mouse cortex with a threshold giving 232 (resp. 85) estimated change-points in rows and columns.

  Figure 10 displays the different change-point locations that can be obtained for these different values of the threshold.

Fig. 10 .

 10 Fig. 10. Plots of the estimated change-points locations (x-axis) for different thresholds from 0.5% to 50% by 0.5%. The estimated change-point locations associated to threshold which are multiples of 5% are displayed with black points.

  Statistical Performances. For each value of σ in {1, 2, 5}, we generated 1000 matrices following the model described in Section 4.1 with (t1,k ) 1≤k≤K 1 = ([nk/(K 1 +1)]+1) 1≤k≤K 1 and (t 2,k ) 1≤k≤K 2 = ([nk/(K 2 +1)]+1) 1≤k≤K 2 ,where [x] denotes the integer part of x and K 1 = K 2 = 4 and n = 500. Figure 2 (middle) displays the mean square error n -2 B -B 2 2

(Top)

.