Mariella Dimiccoli
email: mdimiccoli@iri.upc.edu

Lluís Garrido
email: lluis.garrido@ub.edu

Guillem Rodriguez-Corominas

Herwig Wendt
email: herwig.wendt@irit.fr

Graph Constrained Data Representation Learning for Human Motion Segmentation

Recently, transfer subspace learning based approaches have shown to be a valid alternative to unsupervised subspace clustering and temporal data clustering for human motion segmentation (HMS). These approaches leverage prior knowledge from a source domain to improve clustering performance on a target domain, and currently they represent the state of the art in HMS. Bucking this trend, in this paper, we propose a novel unsupervised model that learns a representation of the data and digs clustering information from the data itself. Our model is reminiscent of temporal subspace clustering, but presents two critical differences. First, we learn an auxiliary data matrix that can deviate from the initial data, hence confers more degrees of freedom to the coding matrix. Second, we introduce a regularization term for this auxiliary data matrix that preserves the local geometrical structure present in the high-dimensional space. The proposed model is efficiently optimized by using an original Alternating Direction Method of Multipliers (ADMM) formulation allowing to learn jointly the auxiliary data representation, a nonnegative dictionary and a coding matrix. Experimental results on four benchmark datasets for HMS demonstrate that our approach achieves significantly better clustering performance then state-of-the-art methods, including both unsupervised and more recent semi-supervised transfer learning approaches 1 .

Introduction

Human-motion segmentation (HMS) aims at breaking a continuous sequence of data depicting human actions and activities into a set of internally coherent temporal segments, and has emerged as a suitable first step in early processing of untrimmed videos for human action recognition [START_REF] Jhuang | Towards understanding action recognition[END_REF]. Despite the fact that HMS has been intensively investigated so far [START_REF] Tierney | Subspace clustering for sequential data[END_REF]15,[START_REF] Xia | Human motion segmentation via robust kernel sparse subspace clustering[END_REF], performance have not yet reached the level of accuracy required in real world applications.

Because labelling large amounts of videos for creating a 1 Project page at URL https://github.com/mdimiccoli/GCRL-for-HMS/.

training set for HMS is very expensive and time-consuming, the problem has been traditionally addressed through unsupervised learning techniques [22,2,[START_REF] Lv | Recognition and segmentation of 3-d human action using hmm and multi-class adaboost[END_REF]. Several approaches cast HMS as a clustering problem by relying on the framework of subspace clustering [START_REF] Vidal | Subspace clustering[END_REF]. The key idea of subspace clustering is to learn an effective representation in form of a coding matrix from which it is possible to construct an affinity matrix that allows to separate data points according to their underlying low-dimensional subspaces. These subspaces are assumed to correspond with the different motions. However, these clustering approaches typically do not take into account temporal continuity and, moreover, they are sensitive to noise. Some recent subspace clustering based approaches have focused on dealing with noise and outliers naturally present in the data, as well as with its sequential nature [6, [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF][START_REF] Tierney | Subspace clustering for sequential data[END_REF]15,[START_REF] Li | Robust subspace clustering with independent and piecewise identically distributed noise modeling[END_REF]. Outside the subspace clustering framework, other temporal clustering approaches work directly on the temporal data instead of on an underlying representation to estimate cluster labels [START_REF] Smyth | Probabilistic model-based clustering of multivariate and sequential data[END_REF]36,[START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Nevertheless, the main challenge of all these unsupervised approaches remains how to cope with the lack of prior knowledge, which could easily cause unpredictable segmentation output. Recently, this problem has been partially addressed by transfer subspace techniques, that leverage prior information from source data to improve clustering accuracy on the target data [31,[START_REF] Wang | Low-rank transfer human motion segmentation[END_REF]40]. These methods have indeed reported improved performance and currently represent the state of the art in HMS. Instead of leveraging information from source annotated data, here we propose to cope with the lack of prior knowledge in an unsupervised fashion by further exploiting the local geometric structure present in the original highdimensional space. Our a priori on the data samples is that they are drawn from an unknown underlying similarity graph, where nodes in the graph correspond to individual frames, edges to connections between nodes, and communities (i.e., groups of nodes that are interconnected by edges with large weights) roughly correspond to different motions. Our model aims at jointly learning the subspaces and a compatible representation for the similarity graph. Since we are interested in pushing the limits of state-of-the-art HMS, we build on the framework of Temporal Subspace Clustering (TSC) [15], that has proven to be the most effective subspace clustering technique for temporal data without transfer. The rational underlying TSC is to find a coding matrix whose affinity graph leads to good clustering with temporal coherence. This is ensured by a least squares regression formulation where the original data are approximated through a nonnegative dictionary and code under block diagonal and temporal Laplacian regularization terms that lend the model its global subspace structure and temporal smoothness, respectively. However, when the original data is noisy or corrupted, nuisances propagate also to the coding matrix. To avoid this, we allow the auxiliary data representation to vary, hence conferring to the coding matrix more degrees of freedom and robustness to subspace assumption violations and noise. Our prior for the auxiliary data matrix is that its columns act as representation vectors for the affinity matrix (graph) of the original data. Therefore, we introduce a graph regularization term, that controls the distance between the affinity graphs of the original and auxiliary data. This effectively allows the auxiliary data to cope with data nuisances while preserving the local geometrical structure of the original data. We propose an efficient optimization formulation for our model, based on the Alternating Direction Method of Multipliers (ADMM).

Our main contributions are summarized as follows: 1) we propose a novel approach for HMS, which jointly learns an auxiliary data matrix, a non-negative dictionary and a coding matrix under temporal, block-diagonal and graph constraints via an original ADMM formulation, 2) we introduce an original graph regularization term that preserves the geometrical structure of the original data in highdimensional space, 3) we present a comprehensive analysis of our model, including ablation study and sensitiveness analysis, 4) we achieve significant performance improvements (up to ⇡ 20% accuracy and ⇡ 5% NMI) over the state-of-the-art methods on four public benchmarks for HMS.

Related work

Subspace clustering. Subspace clustering assumes sampled data to be drawn from a union of multiple subspaces, hence breaking the assumption that all of the clusters in a dataset are found in the same (sub)set of dimensions. Under this assumption, the sampled data obey the so called self-expressiveness property, i.e., each data point in a union of subspaces can be represented by a linear combination of other points in the dataset. This can be formulated as

X = XZ,
where Z 2 R r⇥N is the representation coefficient or coding matrix. The goal of subspace clustering is to find a cod-ing matrix Z that minimizes the reconstruction error under some regularity constraints. The final clusters are obtained by applying spectral clustering to the affinity matrix of Z. Different variants of subspace clustering have been proposed that impose different constraints on Z for this purpose. Least Square Regression (LSR) [START_REF] Lu | Robust and efficient subspace segmentation via least squares regression[END_REF] uses an `2 norm regularizer for Z. Low-rank representations (LRR) [START_REF] Liu | Robust recovery of subspace structures by low-rank representation[END_REF] impose on Z to have small rank and are therefore more robust to corrupted data. Sparse subspace clustering (SSC) algorithms [6] enforce a sparsity constraint on Z and are also applicable to noisy data that do not cluster perfectly into subspaces. Probabilistic methods usually model the data points using a mixture of probabilistic PCAs and are often more robust to noise and outliers than LSR and LRR [START_REF] S Derin Babacan | Probabilistic low-rank subspace clustering[END_REF]. Information theoretic learning based framework such as [START_REF] Li | Robust subspace clustering with independent and piecewise identically distributed noise modeling[END_REF] have proved to be robust to piecewise identically distributed noise. One of the shortcomings of all of the above methods is that they neglect the temporal information that is implicitly encoded in sequential data.

Temporal data clustering. Temporal data clustering methods explore dynamic regularities underlying temporal data in an unsupervised learning fashion. Following the taxonomy introduced in [38], we can distinguish three general approaches to temporal data clustering, depending on the data dependency treatment: model-based [START_REF] Smyth | Probabilistic model-based clustering of multivariate and sequential data[END_REF]36], temporal proximity based [START_REF] Keogh | On the need for time series data mining benchmarks: a survey and empirical demonstration[END_REF] and representation based algorithms [START_REF] Dimitrova | Motion recovery for video content classification[END_REF][START_REF] Chen | Motion trajectory matching of video objects[END_REF][START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF]. Model-based and temporal proximity based approaches work directly on temporal data and deal with the temporal correlation during clustering analysis. A representative example in the context of HMS is Aligned Clustering analysis (ACA) [39], that extends kernel k-means and spectral clustering for time series clustering, formulating the problem within an optimization framework.

Instead, representation based approaches capture temporal data dependencies via a sparse representation. Within this class fall all methods that tried to adapt the subspace clustering framework [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF] to sequential data. Tierney et al. [START_REF] Tierney | Subspace clustering for sequential data[END_REF] proposed an Ordered Subspace Clustering (OSC) method that incorporates a neighbour penalty term to enforce consecutive frames to have similar representation. In Sequential Subspace Clustering [START_REF] Wu | Ordered subspace clustering with block-diagonal priors[END_REF] a quadratic normalizer is imposed on the sparse coefficients to model the temporal correlation among the data points. Clustering accuracy is enhanced by incorporating into the model a block-diagonal prior for the spectral clustering affinity matrix. Temporal Subspace clustering (TSC) [15] learns jointly a dictionary and a representation code with a Laplacian temporal regularization. The use of a dictionary allows for more expressive coding. These temporal clustering methods cannot explicitly and efficiently deal with corrupted or partially missing data. A probabilistic approach to TSC was proposed in [START_REF] Gholami | Probabilistic temporal subspace clustering[END_REF], that models temporal dependencies by using Gaussian Process (GP) priors [25] and is more robust to outliers.

Approach

Fig. 1 provides an overview of our approach for HMS. Our method take as input frame-level features extracted from the video sequence, and firstly computes an affinity graph G 0 in the original high-dimensional feature space. Then, it optimizes jointly a non-negative dictionary D and coding matrix Z, and the auxiliary data representation X under temporal, block-diagonal and graph based regularization constraints. The final segmentation is obtained by applying normalized cut on the affinity graph of Z.

Problem formulation

Preliminaries. Let us denote X = [x 1 , x 2 , ..., x N ] 2 R n⇥N +
time-series data corresponding to motion features with N time-steps in an n-dimensional Euclidean space. Without loss of generality, we assume that X is positive valued and normalized to the interval [0, 1]. Furthermore, we assume that X is drawn from a union of M subspaces where r is the number of atoms in the dictionary. The use of the dictionary in the above formulation instead of the selfrepresentation model (i.e., X ⇡ XZ) allows to learn a more expressive coding matrix Z from the data [15]. Z can then be used to construct an affinity graph for subspace clustering, and the temporal clustering results can be obtained by applying an efficient clustering algorithm, e.g. normalized cuts. Yet, noisy or corrupted data samples contained in X as well as portions of temporal segments that deviate from the subspace assumption may nonetheless propagate to Z and D and lead the model astray. Graph constrained data representation. Our key insight is therefore to let DZ approximate auxiliary data X 2 R n⇥N + instead of X, that we optimize jointly with D and Z. This lends the model extra flexibility to find more expressive codes Z. Our prior for the auxiliary data X is that they act as representation vectors for the affinity matrix (graph) of the original data X. We additionally assume that communities in the similarity graph roughly correspond to subspaces, where nodes correspond to individual frames and edges correspond to similarity values between frames. Thus, the auxiliary data points X approximately reproduce the local geometry of the original data points X but can otherwise move in R n⇥N + to find a representation that is more compatible with the subspace clustering model, and to unburden Z and D from dealing with corrupted data points.

{ m } M m=1 of unknown dimensions dim( m ) = M m (0 < M m < n),
The assumptions on X are expressed through a regularization term that is reminiscent of graph embedding [START_REF] Goyal | Graph embedding techniques, applications, and performance: A survey[END_REF] and penalizes deviations of the graph of the new data points X from the graph of X. Specifically, the graphs of X and X are represented by their respective R N ⇥N affinity matrices G 0 = S(X) and G = S( X). These are defined elementwise via a pairwise similarity, upon normalization to unit total weight,

(S(X)) kj = w kj /W, W = X (k,j)2E w kj w kj = exp (1 d(X(j), X(k)))/h , ( 1
)
where E is the set of pairwise relations, d(•, •) is the cosine distance and h the filtering parameter of the exponential function. The normalized pairwise similarities can be interpreted as empirical joint probability distributions between pairs of vertices p(k, j). Therefore, Eq. ( 1) defines a distribution p(•, •) over the space N ⇥ N . Then, to preserve first-order neighbor relations, we minimize the KL-divergence of the two joint probability distributions pG0 (•, •), p G(•, •) that correspond with the original and learnt auxiliary data graphs, respectively. Upon removal of constant terms, this yields the cross entropy (CE

) divergence CE(p G0 (•, •), p G(•, •)) = P (k,j) pG0 (k, j) log(p G(k, j
)); we denote this additional loss term as L G and name it graph constrained loss. Our CE term is easy to minimize and strongly convex, guaranteeing convergence of our algorithm, see Sec. 4.4. Model. With this graph constrained loss and the new data fidelity term on X, we obtain our final model:

min Z,D, X L G (S( X), G 0 ) + 0 || X DZ|| 2 F + 1 ||Z|| 2 F + 2 f (Z) s.t. Z 0, D 0, ||d i || 2 2  1, i = 1, . . . , r.
(2) The Frobenius norm on Z is used to enforce block diagonal structure in Z, see e.g. [START_REF] Lu | Subspace clustering by block diagonal representation[END_REF]. Moreover, f (Z) = tr(Z T L T Z) is a weighted Laplacian `2 regularization for temporal consistency defined as in [15]. The unit norm constraints on d i resolve the scale ambiguity between D and Z.

Optimization

To solve the objective function in (2), we devise an optimization algorithm based on ADMM. Upon introduction of a set of auxiliary variables Y , U and V , (2) is formulated equivalently as

min X,U,V,Y,D,Z L G (S( X), G 0 ) + 0 ||Y UV || 2 F + 1 ||V || 2 F + 2 f (V ) s.t. Y = X, U = D, V = Z, Z 0, D 0, ||d i || 2 2  1, i = 1, .
. . , r. The augmented Lagrangian for this problem is

L ⇢ =L G (S( X), G 0 ) + 0 ||Y UV || 2 F + 1 ||V || 2 F + 2 tr(V L T V ) + h⇤ X , Y Xi + h⇤ U , U Di + h⇤ V , V Zi + ⇢ 2 ||Y X|| 2 F + ⇢ 2 ||U D|| 2 F + ⇢ 2 ||V Z|| 2 F s.t. Z 0, D 0, ||d i || 2 2  1, i = 1, . . . , r. (3)
The ADMM algorithm that solves this optimization problem (3) is given by alternately minimizing

L ⇢ w.r.t. each variable Y, V, U, X, Z, D individually [3]. Update for V, U .
The updates for V, U are similar to those of the TSC algorithm of Li et al. [15], but need to be slightly modified to account for the new term L G and the additional regularization parameter 0 . Specifically, the update for V is obtained as the solution to the Sylvester equation

(2 0 U T U + 1 I + ⇢I)V + 2 V L T = 2 0 U T Y ⇤ U + ⇢Z (4)
and the update for U has the closed-form solution

U = (2 0 Y V T ⇤ U + ⇢D)(2 0 V V T + ⇢I) 1 .
Update for Z, D. The updates for V, U are identical to those of the TSC algorithm of Li et al. [15] and are given by the expressions

Z = F + (V + ⇤ V /⇢), D = F + (U + ⇤ U /⇢),
where

(F + (A)) ij = max(A ij , 0). Update for Y . We cancel the gradient of (3) w.r.t. Y r Y L ⇢ = 2 0 Y 2 0 UV + ⇤ X + ⇢(Y X)
and find the update for

Y Y = (2 0 + ⇢) 1 (2 0 UV ⇤ X + ⇢ X).
Update for X. The update for X does not have a closed form solution since one has to solve

X = arg min X L G (S( X), G 0 )+h⇤ X , Y Xi+ ⇢ 2 ||Y X|| 2 F .
This expression is easy to solve using gradient descent. The expression for the gradient is given by

r X L ⇢ = r X L G (S( X), G 0 ) ⇤ X ⇢(Y X),
and thus essentially requires the calculation of the gradient of the cross-entropy loss L G composed with the similarity function S. Both Y and X are normalized to [0, 1], as X is assumed to be.

Update for the Lagrange multipliers ⇤.

Finally, the Lagrange multipliers are updated as usual as

⇤ U ⇤ U + ⇢(U D), ⇤ V ⇤ V + ⇢(V Z), ⇤ X ⇤ X + ⇢(Y X).
The details for solving Eq. (3) via the above ADMM formulation are summarized in Algorithm 1. Note that the normalization of X and X, in addition to the strong convexity property of the CE function, allows to ensure mild convergence conditions for the graph regularization term, see Sec. 4.4.

Clustering

To produce the temporal segmentation results, we applied the Normalized Cut algorithm [START_REF] Elhamifar | Sparse subspace clustering[END_REF] to the affinity matrix A of the learnt coding matrix Z. We proceed as in the TSC algorithm [15] and define each element of A as the cosine distance between the respective pair of the learnt code vectors

(A(Z)) kj = z T j z k ||z T j || 2 ||z k || 2 . ( 5
U = (2 0Y V T U + D)(2 0V V T + I) 1 Z = F+(V + V / ) D = F+(U + U / ) -Update for Y and X graph embedding step Y = (2 0 + ) 1 (2 0U V X + X) Y =normalize(Y ) normalize to [0, 1] X = arg min X LG(S( X), G0) + X , Y X + 2 ||Y X|| 2 F X =normalize( X) normalize to [0, 1] -Update Lagrange Multiplier U U + (U D) V V + (V Z) X X + (Y X) end Output: Z R r N + coding matrix
Algorithm 1: Solving problem (3) via ADMM

Clustering

To produce the temporal segmentation results, we applied the Normalized Cut algorithm [START_REF] Elhamifar | Sparse subspace clustering[END_REF] to the affinity matrix A of the learnt coding matrix Z. We proceed as in the TSC algorithm [15] and define each element of A as the cosine distance between the respective pair of the learnt code vectors

(A(Z)) kj = z T j z k ||z T j || 2 ||z k || 2 . ( 5
)

Experimental results

Datasets

We used four human motion datasets to evaluate our proposed model, which we detail next. Weizmann action dataset (Weiz). The Weiz dataset [9]

Experimental setup

Dataset. We used the n = 324 dimensional HOG features, extracted frame by frame and concatenated together, pub-

Experimental setup

Dataset. We used the n = 324 dimensional HOG features, extracted frame by frame and concatenated together, publicly available for these benchmark datasets2 , and used the same setting as in [START_REF] Wang | Low-rank transfer human motion segmentation[END_REF]. Evaluation metrics and parameters settings. In all the experiments, we use the clustering accuracy (ACC) and normalized mutual information (NMI) as the evaluation metrics [34]. The clustering accuracy is defined as Each color correponds to a cluster. GT is the ground truth.

cluster label c to a ground truth label g. The NMI is defined as

NMI(⌦, C) = I(⌦, C) (H(⌦) + H(C)/2 ,
where I denotes the mutual information and H the entropy, ⌦ is the set of clusters and C is the set of classes. Higher values indicate better performance for the two metrics. We set the similarity bandwidth parameter h in (2) empirically to h = 0.0015, which leads to visually satisfactory results for the graph affinity matrices G. The hyperparameters ( 0 , 1 , 2 ) of our method are tuned independently for each dataset by grid search. Sensitivity to hyperparameter tuning is reported in Section 4.4 below.

Comparative results

We compare to several state-of-the-art subspace clustering based methods such as LRR [START_REF] Liu | Robust recovery of subspace structures by low-rank representation[END_REF], SSC [START_REF] Elhamifar | Sparse subspace clustering: Algorithm, theory, and applications[END_REF], LSR [START_REF] Lu | Robust and efficient subspace segmentation via least squares regression[END_REF], temporal data clustering approaches such as TSC [15], OSC [START_REF] Tierney | Subspace clustering for sequential data[END_REF], and more recent Transfer Subspace Clustering based methods including Transfer Subspace Segmentation (TSS) [31], Low-rank Transfer Subspace (LTS) [START_REF] Wang | Low-rank transfer human motion segmentation[END_REF], Multimutual transfer subspace learning (MTSL) [40]. For the lat- ter class of methods, we show in Tab. 1 the best results obtained independently of the source data as reported in [START_REF] Wang | Low-rank transfer human motion segmentation[END_REF] and [40]. We computed the best performance on the four benchmark datasets using the method [START_REF] Li | Robust subspace clustering with independent and piecewise identically distributed noise modeling[END_REF] that we named Robust Subspace Clustering (RSC), and report average values over five random seeds obtained by using the code made available by the authors3 . Bold indicates the best performance. Our method consistently obtains considerably better performance, also when compared to the transfer subspace clustering-based methods (TSS, LTS, MTS). Qualitative results are shown in Fig. 2, where we show the estimated clusters for the first sequences of the four datasets, respectively. It can be appreciated that the segmentation obtained with the proposed approach is highly consistent and well aligned with the ground truth. Results obtained with TSC and LTS 4 are less satisfactory and prone to produce spurious erroneous short segments, in particular close to segment transitions (e.g., at the center of the Weizmann dataset sequence), as well as merge segments into one cluster (e.g., at the center of the MAD dataset sequence).

Model study

Graph constrained representation. We illustrate how our method leverages on the local geometric structure of the original data X to produce an improved auxiliary data representation X. Specifically, we compare, for the first subject of the Weizmann dataset, the data X and the graph constrained data representation X for different relative weight of the graph penalty L G in (3) (realized by multiplying/dividing i , i = 0, 1, 2, by a common factor). In Fig. 3 we plot a subset of 4 arbitrarily chosen features of X and X (plotted as time series, left plots), the 50 first features of X (center plots) and a visualization of X by t-SNE [29] (right plots); the vertical bars indicate true motion boundaries. It can be observed that the features X resemble de-noised versions of the original data X, that become smoother and smoother as the relative weight of L G in (3) is decreased and thus X is more free to move from X, while at the same time retaining the global motion feature evolution and, in particular, sharp transitions at motion boundaries (cf., e.g., motion segment 4). The t-SNE visualization corroborates such observations and reveals that the clusters get more and more clearly separated in the high dimensional space. Temporal versus Graph regularization. Fig. 4 provides a more quantitative illustration of the effect of the graph penalty L G (results obtained for the first subject of the Weizmann dataset). Specifically, we plot the norm of the temporal second derivative of X as a measure of smoothness (left plot) and of the norm of X X as a measure of data fidelity (right plot) as functions of temporal regularization parameter 2 , and for different relative weight of L G in (3) (blue, red and yellow color). The results for this experiment indicate that reducing the relative weight of L G leads to auxiliary representations X that are smoother and differ more from X, as expected. The same holds true for increasing 2 . However, it can be appreciated that reducing the relative weight of the graph penalty term (i.e., letting X move more freely) is more effective in controlling the temporal smoothness of X than tuning 2 , while at the same time retaining a significantly tighter fit with the original data X. For instance, decreasing the weight of L G by a factor 10 increases smoothness by a factor 10 0.7 ⇡ 5 8 (left plot) at only 10 0.05 ⇡ 10 15% increased deviation from the original data, while increasing 2 by two orders of magnitude leads to modest smoothness increase of 10 0.3 ⇡ 2 at more than 10 0.08 ⇡ 20% larger deviation of X from X. Noise modeling. In Fig. 5 we illustrate on the first subject of the Weizmann dataset the effectiveness of our model in dealing with noise. Specifically, we add centred Gaussian noise of varying global noise level (standard deviation) to the frames of the sequence before HOG feature extraction. Moreover, the noise variance within each frame is modeled piece-wise constant according to: i) a fixed geometry (staircase from left to right image boundary, results reported in Fig. 5 (left)) and ii) quadrants whose size changes randomly from one frame to the other, emulating temporal burstiness of feature noise (Fig. 5 (right)). Results (NMI) are given for three relative weights for the graph penalty L G and compared with TSC and LTS and indicate that the graph constrained data representation X effectively absorbs parts of the data corruptions and leads to significantly better clustering performance than the TSC and LTS approach that operate directly on the corrupted data X. Parameter sensitivity analysis. In our approach, three key regularization parameters, i.e., 0 , 1 and 2 , need to be manually tuned. To study sensitivity to their variation on the model output, we fix the value of one parameter and vary the other two parameters over 2 to 3 orders of magnitude. Results for the Weizmann dataset are shown in Fig. 6 (left). Clearly, very small values for the temporal regularization weight 2 worsen the performance. This is to be expected because data temporal coherence is not modeled any longer when 2 ! 0. Otherwise, the results are consistent in a large range of parameter values. In particular, our proposed method obtains better NMI performance when 0 2 [0.1, 1.2] and when 2 5. Moreover, similarly to the observations in [15], the performance of our model is found to vary little with the parameter 1 that regularizes the block-diagonal constraint. Most importantly, these results demonstrate that every term in our model is useful for improving performance.

Moreover, Fig. 6 (right) shows the performance of our method for different values of the filtering parameter h of the similarity kernel S that tunes the data graph representations and hence the coupling of the original and the auxiliary data graphs G 0 and G. Segmentation performance are found to remain nearly constant over more than two orders of magnitude for h, which demonstrates that our approach is robust also to the precise choice of value for this parameter. Convergence analysis. As the problem is non-convex, the global convergence of ADMM cannot be guaranteed theoretically. However, following [START_REF] Xu | An alternating direction algorithm for matrix completion with nonnegative factors[END_REF] the convergence property of ADMM can be shown if the first-order necessary conditions are satisfied. As we iterate over alternating the minimization of U, V, D, Z with X fixed, and the minimization of Y, X, the convergence of the algorithm is guaranteed if each of these minimization steps converges. The convergence of U, V, D, Z with X fixed has been proved in [15]. The convergence of the minimization of the graph regularization term is guaranteed by our design choice: the crossentropy is a strongly convex function with respect the variables Y and X, X is always positive since we normalize the features before performing gradient descent and therefore the mild conditions for the convergence of this term are satisfied. We show the convergence curve that demonstrates the convergence of our algorithm (i.e., the update of the Lagrange multiplier associated to the graph regularization term over time) in Fig. 7. Complexity analysis. The most time-consuming parts in our model during optimization is the update of V in TSC that has complexity O(r 2 N ). Denoting by K the number of iterations at convergence, the overall computational com-plexity of our algorithm is O(Kr 2 N ), thus equal to TSC. Ablation study. To validate the effectiveness of the proposed model, we compare the performance of our method with that obtained when the graph based regularization term is removed and X = X is fixed in the fidelity term. This corresponds to the TSC model of [15], but with different values for the parameters i , i = 0, 1, 2, individually optimized here for each dataset. Quantitative results are reported in Tab. 2 for all four datasets5 and confirm that our approach is highly effective and leads to consistent and significant improvements (NMI increase of up to +9%).

Conclusions

We proposed a novel subspace clustering based approach for HMS which jointly learns an auxiliary data representation, a non-negative dictionary and a coding matrix under temporal, block-diagonal and graph constraints via an original ADMM formulation. The rationale underlying our approach is to operate on auxiliary data X instead of the original time-series data X, lending the model extra flexibility for finding more expressive codes Z, while preserving the local geometrical structure of the original data in highdimensional space. We present a comprehensive analysis of our model, including several illustrations of the model behaviour, model design justification, ablation study, parameter sensitivity and convergence analysis. Furthermore, we achieve significant performance improvements of up to ⇡ 20% (accuracy) and ⇡ 5% (NMI) over the state-of-theart methods on four public benchmarks for HMS.

Figure 1 :

 1 Figure 1: Overview of the proposed framework for HMS.

  and aim at grouping the N vectors into p(p M ) sequential motion segments, and cluster these p seg-ments into their respective subspaces. Given a nonnegative dictionary (i.e., set of atoms) D = [d 1 , d 2 , . . . , d r ] 2 R n⇥r + and a coding matrix Z 2 R r⇥N + , the time series data can be approximately represented as X ⇡ DZ,

  (d) Results on MAD dataset N number of frames n feature dimension Input : X R n N + time-series data / * Compute data similarity graph * / X =normalize(X) normalize data to [0, 1] G0 = S(X) R N N + / * Main loop * / while not converged do -Update for U,V,D,Z TSC step V solution to Sylvester equation (4)

Figure 2 :

 2 Figure 2: Clustering results on all datasets (sequences 1).Each color correponds to a cluster. GT is the ground truth.

  We used four human motion datasets to evaluate our proposed model, which we detail next. Weizmann action dataset (Weiz). The Weiz dataset [9] consists of 90 video sequences, each with a different person performing 10 different types of actions in an outdoor environment. Keck gesture dataset (Keck). The Keck dataset[START_REF] Lin | Recognizing actions by shape-motion prototype trees[END_REF] consists of 168 video sequences with 4 people performing 14 military gestures with perturbations in the background. The same gesture is repeated three times by each person. Each person is treated as a motion segmentation task.UT interaction dataset (UT).The UT dataset [24] consists of 20 videos, each containing 6 classes of human-human interactions. Each video sequence is around 60 seconds long. Multi-model Action Detection dataset (MAD). The MAD dataset [11] consists of 40 sequences captured from 20 subjects. Each subject performs 35 activities (e.g., jumping, kicking) continuously in two different indoor environments, each corresponding to a different sequence. Each person is treated as a motion segmentation task.

  Figure 2: Clustering results on all datasets (sequences 1).Each color correponds to a cluster. GT is the ground truth.

Figure 3 :

 3 Figure 3: Illustration of X for different weight for the graph penalty (from top to bottom). Left: arbitrary set of 4 features, Center: first 50 features. Right: t-SNE visualization.

Figure 4 :

 4 Figure 4: Frobenius norm of the temporal second derivative (left) and difference with data X (right) of auxiliary representation X for different values of the regularization term L G and temporal regularization 2 .

Figure 5 :

 5 Figure 5: Clustering performance (Weiz, subject 1) upon addition of centred piece-wise iid Gaussian noise to each frame, as a function of noise level: fixed piece-wise geometry (left) and randomly varying piece-wise geometry (right).

Figure 6 :

 6 Figure 6: Sensitivity analysis for 0 , 1 , 2 in (3), and for the graph filtering parameter h.

Figure 7 :

 7 Figure 7: Convergence analysis: ||Y X|| F (left) and ||Y X|| 1 (right) as a function of iteration number (average over UT dataset for 5 seeds).

Table 1 :

 1 ) Clustering comparison results in terms of ACC and NMI on four HMS datasets. For the transfer subpaces learning approaches (TSS, LTS and MTS), we reported the best results and the source dataset used to obtain them between parenthesis (M,K, U, and W stands for MAD, Keck, UT, and Weiz respectively). The best clustering results are denoted in bold.

	Method	ACC	NMI	Method	ACC	NMI	Method	ACC	NMI	Method	ACC	NMI
	LRR	0.4382 0.3638	LRR	0.4862 0.4297	LRR	0.4051 0.4162	LRR	0.2249 0.2397
	RSC	0.4112 0.4894	RSC	0.3485 0.3252	RSC	0.3664 0.1881	RSC	0.3730 0.3418
	OSC	0.7047 0.5216	OSC	0.5931 0.4393	OSC	0.6877 0.5846	OSC	0.5589 0.4327
	SSC	0.6009 0.4576	SSC	0.3858 0.3137	SSC	0.4998 0.4389	SSC	0.4758 0.3817
	LSR	0.5093 0.5091	LSR	0.4548 0.4894	LSR	0.4322 0.5183	LSR	0.3667 0.3979
	TSC	0.6111 0.8199	TSC	0.4781 0.7129	TSC	0.5340 0.7593	TSC	0.5556 0.7721
	TSS(M) 0.6208 0.8509	TSS(M) 0.5395 0.8049	TSS(W) 0.5944 0.7878	TSS(K) 0.5792 0.8286
	LTS(K) 0.6391 0.8599	LTS(M) 0.5509 0.8226	LTS(M) 0.6299 0.8128	LTS(U) 0.5980 0.8211
	MTS(K) 0.6436 0.8371	MTS(M) 0.6010 0.8270	MTS(M) 0.6433 0.8239	MTS(U) 0.6163 0.8314
	Ours 0.8501 0.9053	Ours 0.7864 0.8325	Ours	0.8700 0.8256	Ours 0.8297 0.8471
	± 0.0048 ± 0.0017 (a) Results on Weiz dataset	±0.0069 ±0.0041 (b) Results on Keck dataset	± 0.0022 ± 0.0018 (c) Results on UT dataset		± 0.0042 ± 0.0028

Table 2 :

 2 Ablation study results.

https://github.com/wanglichenxj/Low-Rank-Transfer-Human-Motion-Segmentation

https://github.com/YuanmanLi/github-MWEE

https://github.com/wanglichenxj/Low-Rank-Transfer-Human-Motion-Segmentation

Results obtained using 5 random seeds and grid search to achieve best performance in terms of NMI, which explains quantitative differences with results from[START_REF] Wang | Low-rank transfer human motion segmentation[END_REF] 40] shown in Tab.1.

Acknowledgements

Work partially funded by projects MINECO/ERDF RyC, PID2019-110977GA-I00, JAEINT19-EX-0014, RED2018-102511-T.


Charles-Gérard Lucas

Patrice Abry

Herwig Wendt
email: herwig.wendt@irit.fr

Gustavo Didier
email: gdidier@tulane.edu

Bootstrap for testing the equality of selfsimilarity exponents across multivariate time series

Because of the ever-increasing collections of multivariate data, multivariate selfsimilarity has become a widely used model for scale-free dynamics, with successful applications in numerous different fields. Multivariate selfsimilarity exponent estimation has therefore received considerable attention, with notably an original procedure recently proposed and based on the eigenvalues of the covariance random matrices of the wavelet coefficients at fixed scales. Expanding on preliminary work aiming to test for the equality of the selfsimilarity exponents in bivariate time series, we propose and study here a truly multivariate procedure that permits, from a single observation of multivariate time series, to test for the equality of several, possibly many, selfsimilarity exponents. It is based on an original bootstrap procedure, applied in a multivariate time-scale domain and designed to effectively capture the scale-dependent joint covariance structure of multivariate wavelet coefficients as well as the associated wavelet eigenvalue structures. Extensive simulations conducted on synthetic data, modeled by operator fractional Brownian motions, the reference multivariate selfsimilarity model, permit to show that the proposed multivariate timescale domain bootstrap based test yields the targeted significance level under the null hypothesis (all selfsimilarity exponents are equal) and to assess the power of the test for several alternative hypotheses. This analysis leads us to conclude that the proposed test for the equality of multivariate selfsimilarity exponents is effective and ready for use on (a single time series of) real data.

I. INTRODUCTION

Univariate selfsimilarity. The selfsimilarity paradigm [START_REF] Samorodnitsky | Stable non-Gaussian random processes[END_REF] provides a relevant framework for the analysis and modeling of temporal dynamics in a wide range of applications very different in nature. Fractional Brownian motion (fBm), defined as the only Gaussian, self-similar stochastic process with stationary increments [START_REF] Mandelbrot | Fractional Brownian motion, fractional noises and applications[END_REF], [START_REF] Pipiras | Long-Range Dependence and Self-Similarity[END_REF], is probably the most well-known selfsimilarity model. It quantifies temporal dynamics in data by means of a single selfsimilarity parameter H ∈ (0, 1). The estimation of H is thus central in applications, where it is involved in various classical signal processing tasks (e.g., detection, classification, diagnosis,. . . ), and a large body of literature has been dedicated to this topic, see, e.g., [START_REF] Pipiras | Long-Range Dependence and Self-Similarity[END_REF]- [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF] and references therein. It is well documented that wavelet analysis plays a central role in this context, and has been shown to lead to theoretically well-grounded, accurate and robust estimations for H [START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF], [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]. Most of the work and success in applications remained limited to the modeling and analysis of univariate, scalar-valued time series so far. Yet, in an increasing number of areas, the available data naturally consist of several, jointly collected time series that call for adequate multivariate selfsimilarity models. Selfsimilarity in multivariate time series. Operator fractional Brownian motion (ofBm) was recently proposed as a model for multivariate selfsimilarity [START_REF] Maejima | Operator-self-similar stable processes[END_REF]- [START_REF]Exponents, symmetry groups and classification of operator fractional Brownian motions[END_REF]. Practically speaking, this multivariate extension of fBm can be interpreted as a mixture of M possibly correlated fBm, each with possibly distinct selfsimilarity exponents H m , m = 1, . . . , M . For this model, classical univariate-like statistical methods such as entry-wise sample covariance matrix analysis, Fourier or wavelet spectra do not yield relevant results. Recently, a statistical procedure was devised that instead relies on the scaling properties of the eigenvalues of the scale dependent sample wavelet covariance matrices [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]- [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. It has been shown that this approach is effective for the joint estimation of the vector of selfsimilarity exponents H = (H 1 , . . . , H M ). The use of this novel tool in real-world applications immediately leads to a fundamental question: are observed differences in the selfsimilarity exponents (H 1 , . . . , H M ) statistically significant, or are they within the natural statistical fluctuations of the estimation? A rigorous answer to this question could be achieved by formulating it as a statistical test of the hypothesis H 0 :

H 1 = H 2 = . . . = H M against the alternative H 1 : not H 0 . Yet, so far, this has only been considered in restricted settings, non mixed ofBm in [START_REF] Wendt | Multivariate Hadamard self-similarity: testing fractal connectivity[END_REF], mixed but non correlated ofBm in [START_REF] Abry | Two-step wavelet-based estimation for Gaussian mixed fractional processes[END_REF] or for pairs of time series only in [START_REF] Wendt | Wavelet domain bootstrap for testing the equality of bivariate self-similarity exponents[END_REF], but never for ofBm with multiple, possibly many, correlated and mixed components. Goals, contributions and outline. The present work constructs and assesses a bootstrap driven statistical test, from a single finite size observation of multivariate time series, for the equality of the selfsimilarity exponents estimated on the M ≥ 2 components. As a first significant contribution, we expand on the wavelet estimation procedure originally proposed in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] (cf. Section II-B) to cope with the estimation bias induced by the so-named eigenvalue repulsion effect [START_REF] Anderson | An Introduction to Random Matrices[END_REF]- [START_REF] Wendt | Bootstrap-based bias reduction for the estimation the self-similarity exponents multivariate time series[END_REF] (cf. Section III-A).

Building up on the wavelet domain bootstrap proposed in [START_REF] Wendt | Wavelet domain bootstrap for testing the equality of bivariate self-similarity exponents[END_REF], the second main contribution consists of the bootstrap scheme, designed to preserve the time-scale joint (multivariate) dependence structure of the wavelet coefficients and thus the joint statistics of the (H 1 , . . . , H M ), rather than focusing on the univariate component-wise dependence structure only as in [START_REF] Wendt | Wavelet domain bootstrap for testing the equality of bivariate self-similarity exponents[END_REF] (cf. Section III-B).

Monte Carlo simulations, conducted using OfBm synthetic data (cf. Section II-A for definitions) with M = 4 and 10 components, and for several sample sizes, are used to show that the proposed multivariate wavelet domain bootstrap procedure reproduces well the targeted test significance level under the null hypothesis (cf. Section IV-B), and to study the power of the proposed test under several alternative hypothesis (cf. Section IV-C). Our results demonstrate that the proposed test effectively permits to test for the equality of selfsimilarity across components from a single observation of multivariate data. OfBm synthesis, selfsimilarity exponent estimation and test procedures are implemented by the authors and will be made available at the time of publication.

II. MULTIVARIATE SELFSIMILARITY

A. Model: operator fractional Brownian motion

For simplicity, we make use of a subclass of time reversible ofBm, see [START_REF] Maejima | Operator-self-similar stable processes[END_REF]- [START_REF]Exponents, symmetry groups and classification of operator fractional Brownian motions[END_REF] for an introduction of ofBm in all generality. It is a natural multivariate extension of fBm consisting of a multivariate Gaussian self-similar process with stationary increments.

Let X {X H1 (t), . . . ,X H M (t)} t∈R be a collection of M fBm components, each characterized by their selfsimilarity exponents H = (H 1 , . . . , H M ), 0 < H 1 ≤ . . . ≤ H M < 1. These components are possibly correlated according to a pointwise M × M covariance matrix Σ X . If we denote by σ 2 m and ρ m,m the variances of each component and their respective correlation coefficients, the entries of Σ X can be written as (Σ X ) m,m = σ m σ m ρ m,m . It is in general not possible to choose H and Σ X independently, see [START_REF] Didier | Integral representations and properties of operator fractional Brownian motions[END_REF].

Further, let denote W a real-valued and invertible M ×M matrix (the mixing matrix). Then, the M -variate stochastic process, referred to as ofBm, is defining as a linear mixture of these M correlated fBm:

Y {Y H,Σ X ,W 1 (t), . . . , Y H,Σ X ,W M (t)} t∈R W {X H1 (t), . . . , X H M (t)} t∈R = W X. (1)
It can be shown that [START_REF] Samorodnitsky | Stable non-Gaussian random processes[END_REF] implies that Y has the following key properties: i) its increments are stationary and ii) it satisfies the multivariate selfsimilarity relation:

∀a > 0 : { Y H,Σ X ,W 1 (t), . . . , Y H,Σ X ,W M (t) } t∈R fdd { a H (Y H,Σ X ,W 1 (t/a), . . . , Y H,Σ X ,W M (t/a)) } t∈R , (2)
where fdd stands for the convergence of finite dimensional distributions, a H +∞ k=0 log k (a)H k /k!, and H W diag(H)W -1 denotes the so-named selfsimilarity exponent matrix. When W is diagonal, the multivariate selfsimilarity relation simplifies to component-wise selfsimilarity relations as in [START_REF] Amblard | of multivariate fractional motion[END_REF].

B. Wavelet based joint estimation of H m , m = 1, . . . , M
The central goal of selfsimilarity analysis is to estimate the vector of exponents H = (H 1 , . . . , H M ) from the Mvariate time series Y . It is only when W and thus the selfsimilarity matrix H are diagonal that each component Y m is a fBm and H can be inferred with conventional, univariatelike approaches [START_REF] Wendt | Multivariate Hadamard self-similarity: testing fractal connectivity[END_REF], [START_REF] Amblard | of multivariate fractional motion[END_REF]. Otherwise, each component is a mixture of fBms, and such estimations are strongly biased [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. To obtain consistent estimates for H in all generality, a statistical procedure based on the multivariate discrete wavelet transform was proposed [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]. Multivariate wavelet coefficients are defined as

(D(2 j , k)) D Y (2 j , k) = (D Y1 (2 j , k), . . . , D Y M (2 j , k)), ∀k ∈ Z, ∀j ∈ {j 1 , . . . , j 2 }, with ∀m ∈ {1, . . . , M }: D Ym (2 j , k) = 2 -j/2 ψ 0 (2 -j t - k)|Y m (t)
and ψ 0 is a mother wavelet [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF], [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The procedure relies on first estimating wavelet coefficient M × M sample covariance matrices at each scale 2 j independently:

S(2 j ) 1 n j nj k=1 D(2 j , k)D(2 j , k) * , w = 1, . . . , 2 j , (3)
with N the sample size, j the octave and n j = 2 -j N . Second, eigen-decompositions of these S(2 j ) are performed at each scale independently, yielding eigenvalues {λ 1 (2 j ), . . . , λ M (2 j )}. Then, it was shown in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF], [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] that the wavelet estimator ( Ĥ0 1 , . . . , Ĥ0 M ) for (H 1 , . . . , H M ), defined as weighted linear regressions of log 2 λ m (2 j ) against analysis scales 2 j1 , . . . , 2 j2 ,

Ĥ0 m =   j2 j=j1 w j log 2 λ m (2 j )   2 - 1 2 , ∀m = 1, . . . , M, (4)
with w j classical weights satisfying j jw j = 1 and j w j = 0 (cf. [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]), converges asymptotically in probability, N → ∞, under mild assumptions, ( Ĥ0 1 , . . . , Ĥ0 M )

P → (H 1 , . . . , H M ).

III. MULTIVARIATE WAVELET DOMAIN DE-BIASED ESTIMATION AND BOOTSTRAP HYPOTHESIS TEST

A. De-biased estimation

The estimation procedure described above shows very satisfactory practical performance but can suffer from a bias in estimation in certain situations that impairs testing the equality in H. It is indeed well known that eigenvalue decomposition applied to estimated covariance matrices is impaired by the socalled repulsion effect: the gap between estimated eigenvalues is on average larger than the gap between actual values, and more importantly, two equal eigenvalues are not estimated equal. This bias is all the more pronounced as the sample size available for covariance estimation is small. In wavelet transforms, the number of wavelet coefficients actually available for covariance estimation (cf. Eq. 3) decreases as scale increases. This yields a scale dependent bias in the estimation of λ m (2 j ) that results in a bias in Ĥ0 m .

To reduce such bias, one first original contribution of the present work is to use for wavelet coefficient covariance matrix estimation the same number of wavelet coefficients at each scale. Specifically, a reference scale 2 J = 2 j2 is chosen here as the upper limit of the range where the linear fit is performed. For scale 2 j , a collection of 2 j-j2 covariance matrices are estimated from non-overlapping wavelet coefficients snapshots of equal size n j2 = 2 -J N , for w = 1, . . . , 2 j-j2 :

S (w) (2 j ) 1 n j2 wnj 2 k=1+(w-1)nj 2 D(2 j , k)D(2 j , k) * . (5)
The eigenvalues of each of these matrices S (w) (2 j ), denoted as {λ

(w) 1 (2 j ), . . . , λ (w)
M (2 j )}, have similar repulsion at all scales j ∈ {j 1 , . . . , j 2 }. Their logarithms are averaged

ϑ m (2 j ) 2 j2-j 2 j-j 2 w=1 log 2 (λ (w) m (2 j )) (6)
and a new wavelet estimator ( Ĥ1 , . . . , ĤM ) for (H 1 , . . . , H M ) is defined by means of weighted linear regressions of ϑ m (2 j ) against analysis scales 2 j1 , . . . , 2 j2 as

Ĥm =   j2 j=j1 w j ϑ m (2 j )   2 - 1 2 , ∀m = 1, . . . , M. (7)

B. Multivariate wavelet domain bootstrap test

Asymptotic normality and χ 2 statistics. It was shown in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]- [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF], that under mild hypothesis, the wavelet estimator ( Ĥ0 1 , . . . , Ĥ0 M ) was characterized by jointly Gaussian multivariate statistics. Numerical simulations suggest that the new wavelet estimator proposed here, ( Ĥ1 , . . . , ĤM ) for (H 1 , . . . , H M ) also displays jointly Gaussian multivariate statistics. Such empirical evidence will be theoretical studied elsewhere and jointly Gaussian multivariate statistics are assumed to hold for the rest of the present work.

Given joint normality for the selfsimilarity exponents H = (H 1 , . . . , H M ), a natural statistic to test the null hypothesis

H 0 (H 1 = H 2 = . . . = H M ) reads: T = ( Ĥ -Ĥ 1 M ) T Σ -1 Ĥ ( Ĥ -Ĥ 1 M ) (8)
where Σ Ĥ denotes the M × M covariance matrix of Ĥ, while Ĥ stands for the average across components:

Ĥ = 1/M M m=1
Ĥm and 1 M = (1, . . . , 1). Asymptotic joint normality of Ĥ = ( Ĥ1 , . . . , ĤM ) implies that the test statistic T follows a χ 2 distribution with M -1 degrees of freedom under the null hypothesis H 0 .

The key issue is that the matrix Σ Ĥ is a priori unknown and therefore needs to be replaced with an estimate (see [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]- [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF] for theoretical elements). Because it is intended here that the proposed test can be put at work from a single observation of the multivariate data, Σ Ĥ cannot be estimated by averaging across realizations, thus calling for a bootstrap procedure [START_REF]Bootstrap Techniques for Signal Processing[END_REF] in the spirit of those developed in [START_REF] Wendt | Wavelet domain bootstrap for testing the equality of bivariate self-similarity exponents[END_REF], [START_REF] Wendt | Bootstrap-based bias reduction for the estimation the self-similarity exponents multivariate time series[END_REF], [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. Wavelet domain bootstrap. To take advantage of the short term temporal and across-scale dependencies of the wavelet coefficients, the bootstrap procedure is implemented in the wavelet domain rather than in the time domain. The vector wavelet coefficients D(2 j , k), k = 1, . . . , n j are resampled jointly (i.e., in a multivariate way) by a block-bootstrap procedure which preserves their scale-time dependence structure, as opposed to a resampling for each component independently, which would not preserve the cross-component dependencies [START_REF]Resampling Methods for Dependent Data[END_REF]. Technically, for each scale 2 j , R block bootstrap resamples D * (r) j = (D * (r) (2 j , 1), . . . , D * (r) (2 j , n j )), r = 1, . . . , R are drawn with replacement, from card(Y )/L B overlapping blocks of size L B , (D(2 j , k), . . . , D(2 j , k + L B -1)), k = 1, . . . , n j .

From each bootstrap sample of wavelet coefficients D * (r) j , bootstrap estimates S * (r) (2 j ), ϑ * (r) m (2 j ) and Ĥ * (r) m are produced using Eqs. [START_REF] Beran | Long-Memory Processes: Probabilistic Properties and Statistical Methods[END_REF][START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Flandrin | Wavelet analysis and synthesis of fractional Brownian motion[END_REF]. Finally, from the R vector bootstrap samples Ĥ * (r) = ( Ĥ * (r) 1 , . . . , Ĥ * (r) M ), r = 1, . . . , R, a covariance matrix Σ Ĥ * is estimated that is expected to approximate the true covariance matrix of H = (H 1 , . . . , H M ). Bootstrap test formulation. From the estimates Ĥ and the bootstrap covariance matrix Σ Ĥ * , the following test statistics is constructed:

T * = ( Ĥ -Ĥ 1 M ) T Σ-1 Ĥ * ( Ĥ -Ĥ 1 M ). ( 9
)
If the covariance structure is well approximated by the blockbootstrap procedure, the bootstrap statistic T * should asymptotically be well approximated by a χ 2 distribution (as verified by numerical results in Section IV-B), hence yielding the following test decision:

d α = 1 (H 0 rejected) : T * ≥ F -1 χ 2 (M -1) (1 -α), (10
) with F χ 2 (M -1) the chi-square cumulative distribution function with M -1 degrees of freedom and α the significance level.

IV. TEST PERFORMANCE ASSESSMENT

A. Monte Carlo Simulation

To assess the accuracy of the proposed bootstrap test procedure and actual reproduction of a targeted significance level by the test under the null hypothesis, and to quantify the power of the test under several representative alternative hypotheses, we conduct Monte Carlo simulations with synthetic data.

Specifically, the numerical simulations make use of N M C = 1000 independent realizations of synthetic M -variate ofBm, with M = 4 or M = 10 components, and with three different sample sizes N = 2 14 , 2 16 , and 2 18 . The covariance matrix Σ X is set to a Toeplitz matrix, with off-diagonal entries all set to ρ = 0.5. The mixing matrix W is chosen randomly in the space of invertible M × M matrices. Both matrices are kept fixed for all the experiments.

Under the null hypothesis H 0 , the selfsimilarity exponents are set to H 1 = . . . = H M = 0.7. Several significance levels α ∈ {0.005, 0.01, 0.05, 0.1, 0.15} are studied.

For the alternative hypothesis, two scenarios are tested. In Scenario1, there are only two different values of H, denoted H 1 and H 2 , amongst the M components. In Scenario1a,

(M 1 , M 2 , M 3 , M 4 ) = (2, 3, 2, 3) (M = 10) and (M 1 , M 2 , M 3 , M 4 ) = (1, 1, 1, 1) (M = 4) with selfsimilarity exponents (H 1 = 0.5, H 2 = H 1 + ∆H, H 3 = H 1 + 2∆H, H 4 = H 1 + 3∆H), ∆H ∈ [0, 0.13].
The significance level is set to α = 0.05.

For estimation, the multivariate wavelet transform is computed with the least asymmetric Daubechies3 wavelets, estimation scales are set to j 1 = 6 and j 2 = log 2 (N ) -5. For the test, R = 500 bootstrap resamples are used with L B = 6 corresponding to the size of the Daubechies3 mother-wavelet time support.

B. Test behavior under the null hypothesis

Fig. 1 reports quantile-quantile plots of the empirical distribution of T * obtained from the N M C independent realizations against the theoretical chi-square distribution with M -1 degrees of freedom. It shows, for three different samples sizes N and two different number of components M , an excellent agreement, hence validating the heuristic assumption in Section III-B that the test statistic T * follows χ 2 distribution. Furthermore, Fig. 2 shows that the empirical distributions of the corresponding p-values are approximately uniform, as should theoretically be the case under the null hypothesis H 0 .

To finish with, Fig. 3 shows, for the three different samples sizes N , and the two different number of components M , that the achieved significance levels, obtained as averages of the binary decisions of the bootstrap test across realizations, reproduce very satisfactorily the targeted test significance.

Overall, these results strongly suggest that, for large ranges of sample sizes N and numbers of components M , the bootstrapped statistics T * reproduce well the distribution of the true statistics T and follow well a χ 2 distribution with M -1 degrees of freedom, and that the proposed bootstrap procedure test is effective in controlling the significance level (type-I error) of the test.

C. Power of the test under different alternative hypothesis

Fig. 4 report the empirical assessment (by average across independent realizations of the test decisions) of the power of the proposed bootstrap-based χ 2 tests as functions of ∆H = |H 2 -H 1 |, for Scenario1 and Scenario2, respectively. In both scenarios, and as expected, the power of the test increases with sample size and with ∆H. From an elementary study of the χ 2 statistics T , and given that the variance in the estimation of H essentially decreases as 1/N (see [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]- [START_REF] Abry | Wavelet eigenvalue regression for n-variate operator fractional Brownian motion[END_REF]), it can be deduced that the power of the test increases as ∆H • √ N . In addition, Fig. 4 (top plots) shows that, for a scenario with two different H and a fixed ∆H, the test power decreases when the number of components M increases.

Further, and interestingly, comparing the results of Sce-nario1a and 1b shows that, for a fixed ∆H, test powers are larger when the groups of components with same H have equal sizes, compared to when groups of components have unbalanced sizes, which can be explained as follows. Under alternative hypothesis, the test statistics T * follows a noncentral χ 2 distribution with M -1 degrees of freedom, with a noncentrality parameter that explicitly depends on the values of the true values H 1 , . . . , H M , as well as on the covariance and mixing matrices Σ X and W . The test power, directly stemming from the distance between the central and noncentral χ 2 distributions, both with M -1 degrees of freedom, is thus fully determined by this noncentrality parameter. For instance, for Scenario1, calculations not detailed here show that the noncentrality parameter is essentially proportional to (M 1 × M 2 )/(M 1 + M 2 ) and thus decreases when |M 2 -M 1 | increases, and so does the test power. Fig. 4 (bottom plots) finally shows that the test also has significant power when there are more than two different values of H.

V. CONCLUSIONS AND PERSPECTIVES

We have devised and assessed a multivariate wavelet domain bootstrap procedure that permits to test, from a finite size single observation of a multivariate time series, whether or not the selfsimilarity exponents are equal. Monte Carlo simulations, based on operator fractional Brownian motion, the reference process to model multivariate selfsimilarity, enabled us to show that the procedure is effective in controlling an a priori prescribed significance level and that the proposed test benefits from significant power under several alternative hypothesis. The proposed test is ready for application to realworld multivariate time series. OfBm synthesis, selfsimilarity exponent estimation and testing procedures are implemented by the authors and will be made available at the time of publication. The work opens interesting perspectives related to i) the assessment of how many different values for H actually exist amongst the M components and ii) the study of large dimensional cases, where the number of components M and sample size N become comparable in magnitude.

Fig. 1 :

 1 Fig. 1: Test statistics under the null hypothesis. Quantilequantile plot of the bootstrap statistic T * against a theoretical χ 2 distribution with M -1 degrees of freedom.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Distributions of p-values under the null hypothesis. Empirical distributions of the p-values of the proposed bootstrap test compared to the Uniform Distribution.

Fig. 4 :

 4 Fig. 4: Power of the test. Empirical assessment of the test powers, obtained as average of the decision of the bootstrap test, as a function of ∆H, for M = 4 (left) and M = 10 (right), for different sample sizes, and with equal or different sizes in the groups of equal components, for Scenario1 (top) and Scenario2 (bottom). Significance level is set to α = 0.05 (magenta dashed line).

ACKNOWLEDGMENT

Supported by PhD Grant DGA/AID (no 01D20019023), ANR-16-CE33-0020 MultiFracs, ANR-18-CE45-0007 MUTATION.