Liang Yu
email: liang.yu@insa-lyon.fr

Jerome Antoni
email: jerome.antoni@insa-lyon.fr

Quentin Leclere
email: quentin.leclere@insa-lyon.fr

Combined Regularization Optimization For Separating Transient Signals From Strong Noise In Rolling Element Bearing Diagnosis

The problem of rolling element bearing diagnosis can be viewed as that of extracting transient signals from strong additive noise. In order to separate these two kinds of signals, the problem is formulated as the minimization of an objective function which consists of two terms, a data fitting term, and a regularization term. The data fitting term reflects the conservation of energy and the regularization term reflects the prior information (initial supposition) on the objective signals to be separated. Since the objective is to separate two kinds of signals, a combined regularization scheme (distinct priori information for each signal with a target for separation) is naturally investigated. In this paper, two kinds of strategies are considered : the first one is a combination of two sparse regularization terms, and the second one is a combination of a low rank regularization term and a sparse regularization term.

Introduction

Rolling element bearings are widely used in rotating machines and their failure is one of most frequent reason for machine breakdown [START_REF] Robert | Rolling element bearing diagnostics-a tutorial[END_REF]. Spectral analysis has for long been used as a basic tool for the task of detection and diagnosis ; nevertheless, being a tool dedicated to stationary signals, it is unable to capture changes in time. Time-frequency analysis has been widely used for rolling element bearing diagnosis [START_REF] Feng | Recent advances in time-frequency analysis methods for machinery fault diagnosis : A review with application examples[END_REF][3] [START_REF] Rubini | Application of the envelope and wavelet transform analyses for the diagnosis of incipient faults in ball bearings[END_REF] due to its good 'localization' capability (Short Time Fourier Transform (STFT) [START_REF] Cocconcelli | Stft based approach for ball bearing fault detection in a varying speed motor[END_REF]), 'Multi-Resolution' (wavelet [START_REF] Rubini | Application of the envelope and wavelet transform analyses for the diagnosis of incipient faults in ball bearings[END_REF]), 'Adaptivity (wavelet pack [START_REF] Zhao | Rolling bearing fault diagnosis using neural networks based on wavelet packet-characteristic entropy[END_REF])', which is useful for extracting machinery fault information from measurements. More recently, dictionary (dictionary provides a more flexible representation as compared with orthogonal transforms of signal [START_REF] Rubinstein | Dictionaries for sparse representation modeling[END_REF]) based sparse representation was brought in rolling element bearing diagnosis [START_REF] Liu | Adaptive feature extraction using sparse coding for machinery fault diagnosis[END_REF] [START_REF] Yang | Fault diagnosis of rolling element bearings using basis pursuit[END_REF] ; the principle is to select a small subset of functions (bases) from a general over-complete dictionary which provides a more flexible decomposition of the signal. In this paper, combined regularization structures (in contrast to single regularization in classical sparse representation model) are investigated for rolling element bearing diagnosis which resort to obtaining a better separation result. An overview of this paper is given as follows.

In order to separate two kinds of signals x = Ac 1 and y = Bc 2 from their mixture signal z = x + y, where N , a general constrained optimization problem (a similar and generalized structure can be found in mathematical optimization literature [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF] [11] [START_REF] He | An alternating direction-based contraction method for linearly constrained separable convex programming problems[END_REF]) is investigated in this paper :

A ∈ R N×N1 , B ∈ R N×N2 and c 1 ∈ R N1 , c 2 ∈ R N2 , x ∈ R N , y ∈ R N , z ∈ R
Minimize Reg 1 (c 1 ) + Reg 2 (c 2 )

Subject to Ac 1 + Bc 2 = z (1)
where Ac 1 + Bc 2 = z is the data fitting term, which indicates the conservation of energy, and Reg 1 (c 1 ) + Reg 2 (c 2 ) is the regularization term, which reflects prior information (initial supposition) assigned to each signal to be separated. The specific forms of the prior Reg 1 (c 1 ) and Reg 2 (c 2 ) are elaborated in the following parts. More specifically, in order to separate a transient signal x ∈ R N from strong additive stationary noise y ∈ R N in measurement z ∈ R N , two kinds of combined regularization forms are investigated :

• Sparse + Sparse Combination : it is assumed that x can be represented in an over-complete dictionary Ψ as x N×1 = Ψ N×M α M×1 , where M ≫ N and α M×1 is a sparse vector, while α is not sparse when y N×1 is decomposed as y N×1 = Ψ N×M αM×1 . Similarly, y N×1 can be represented by one over-complete dictionary Φ as y N×1 = Φ N×M β M×1 , where M ≫ N and β M×1 is a sparse vector, while, β is not sparse when x N×1 is decomposed as x N×1 = Φ N×M β M×1 . The signal separation problem is then formulated as :

Minimize α,β ∥α∥ L 0 + λ ∥β ∥ L 0 Subject to z N×1 = Ψ N×M α M×1 + Φ N×M β M×1 (2)
where L 0 is defined as total number of non-zero elements in a vector [START_REF] Elad | Sparse and redundant representations : from theory to applications in signal and image processing[END_REF], λ is a weight which serves the purpose of adding a trade-off between ∥α∥ L 0 and ∥β ∥ L 0 ; model ( 2) is named Morphological Component Analysis (MCA) in [14][15]. Unfortunately, the noise y N×1 can not be represented parsimoniously under one dictionary (term λ ∥β ∥ L 0 disappears), thus, in this paper, one modified model is proposed for separating the transient signal from strong additive noise :

Minimize α ∥α∥ L 0 Subject to z N×1 = Ψ N×M α M×1 + y N×1 ∥Ψα∥ L 2 ∥y∥ L 2 = c ( 3
)
where L 2 denotes the Euclidean norm, and the small constant c denotes the ratio between the transient signal and noise energies. It is noted that λ is a implicit ratio parameter in model ( 2) and c is a explicit parameter in model (3) ; -more details are argued in section 3.2 . • Low rank + Sparse Combination : Spectrogram is defined as the operation of calculating signals spectrogram matrix. For signal z ∈ R N , Spectrogram{z N×1 } = Spectrogram{x N×1 } + Spectrogram{y N×1 } can be represented as M P×Q = L P×Q + S P×Q . Where L P×Q corresponds to the stationary noise signal and is low rank, and S P×Q corresponds to the transient signal and is sparse. The signal separation problem is then formulated as : Minimize Rank(L)

+ λ ∥S∥ L0 Subject to L + S = M. (4)
It is evident that model ( 4) is similar to form (1), where Ac 1 = L , Ac 2 = S and Reg 1 (c 1 ) = Reg 1 (L), Reg 2 (c 2 ) = Reg 2 (S). The paper is organized as follows. The problem of separating a transient signal from strong additive noise is formally stated in section 2. The Sparse + Sparse combination and Low Rank + Sparse combination models are separately discussed in section 3 and section 4. Sections 5 and 6 presents simulation results and draws the conclusion respectively.

Problem Statement

The problem of separating a transient signal

x ∈ R N from strong additive noise y ∈ R N in measurement z ∈ R N is described as z = x + y (5)
In this work, the transient signal x is modelled as a generalised shot noise process

x(k) = ∑ k (X k * h(k -τ k )), k = 1...N
, where h(k) is the impulse response resulting from a single impact and {X k } and {τ k } are sequences of random variables which account for possibly random amplitudes and random occurrences of the impacts [START_REF] Antoni | The spectral kurtosis : application to the vibratory surveillance and diagnostics of rotating machines[END_REF]. Noise y is assumed to be of zero mean and wide-sense stationary. A simple illustrative example is shown in Fig. 1, where the transient signal (in green color) is merged into strong additive noise (in blue color) : the challenge is to separate the latter from the former. For the current problem (5), firstly, it is assumed that signal x can be represented in a over-complete dictionary Ψ as x = Ψα, and α is a sparse vector, which is obtained by the following equation :

Minimize α ∥α∥ L 0 Subject to x = Ψα (6)
where ∥α∥ L 0 is a prior hypothesis for the coefficients distribution of the transient signal x, and the dictionary Ψ can be built from existing time-frequency atoms or learned from a set of training samples (signals) [START_REF] Elad | Sparse and redundant representations : from theory to applications in signal and image processing[END_REF]. In this paper, dictionary Ψ is built from the discrete Gabor atoms [START_REF] Mallat | A wavelet tour of signal processing : the sparse way[END_REF] :

Ψ = {ψ m,l [n] = g(n -m)exp( i2πln N )} (m,l)∈R 2 (7)
where g(n) is the Hanning time window of period N with unit norm ∥g∥ = 1, centered at n = 0. The atom ψ m,l is translated in time sample m and frequency l, and each atom ψ m,l corresponds to a Heisenberg rectangle with a size (σ t , σ ω ) which satisfies the uncertainty principle σ t σ ω ≥ 1/2. Secondly, it is also known that the transient signal x is buried in strong noise y so that the measured signal reads z = x + y ; this condition is added into Eq. ( 6) to yield :

Minimize α ∥α∥ L 0 Subject to z = Ψα + y (8)
where z = Ψα + y is the data fitting constrain. Equation ( 8) is not sufficient to accomplish the separation task, in order to approach the form of Eq. ( 2), thus, a prior hypothesis for the coefficient distribution of the noise y decomposition is searched. It is noted that noise y can not be represented parsimoniously under one determined Gabor dictionary Ψ, suggesting that the solution α = Ψ † y ( † denotes pseudo inverse of a matrix) of equation y = Ψ α is non-sparse (note that [START_REF] Shental | Sparse representation of white gaussian noise with application to l0-norm decoding in noisy compressed sensing[END_REF] claims that one realization of the noise can be represented parsimoniously under one realization of Gaussian dictionary, which is different with the Gabor dictionary Ψ that was used here), which is congruous with the noise being non-compressible. Because noise y can not be decomposed with sparse coefficients, it seems that Eq. ( 8) is enough for separating the transient signal x from the strong additive noise y. Nevertheless, the trivial solution of Eq. ( 8) is a zero vector, therefore, if the measurement signal is all noise, the coefficient for representing the transient signal will be the most sparse (zero vector). Thus, one more constrained condition is needed to guarantee a non-trivial solution to exist. Thirdly, it is observed that the noise energy σ y 2 of noise y is much larger than the signal energy σ x 2 of signal x, and it is assumed that the signal-to-noise ratio

σ x σ y = ∥x∥ L 2 ∥y∥ L 2 = ∥Ψα∥ L 2 ∥y∥ L 2
= c between the transient signal and noise energies is a small constant. It is assumed that the signal-to-noise ratio is a known parameter, which is considered as a constraint added to Eq. ( 8) to avoid a null sparse solution . Finally, all these constrains are collected together and problem (5) can be solved by the following equation :

Minimize α ∥α∥ L 0 Subject to z = Ψα + y ∥Ψα∥ L 2 ∥y∥ L 2 = c (9)

Matrix-Free Matching Pursuit Algorithm

There exists one difficulty to solve Eq. ( 9) which is explained here by an example. It is assumed that measurement z ∈ R 10000 , the dimension of the dictionary Ψ is 10000 × K and K is much large than 10000, which makes it is impossible to tackle this huge dimension matrix with a desktop computer. Thus, a matrixfree algorithm based on Gabor transform is used to solve this problem. With this algorithm, the dictionary is constructed implicitly by the Gabor transform. As a preparation, the discrete Gabor transform of a signal f of period Nw is represented as [START_REF] Mallat | A wavelet tour of signal processing : the sparse way[END_REF] f

[m, l] = GT ( f [n]) = Nw-1 ∑ n=0 f [n]g Nw [n -mR]exp( -i2πln Nw ) f [n] = IGT ( f [m, l]) = 1 Nw Nw-1 ∑ m=0 g Nw [n -mR] Nw-1 ∑ l=0 f [m, l]exp( i2πln Nw ) ( 10
)
where Nw is the window length (Hanning window is used here), R is the overlapping length between two continuous window blocks, and R is set to be a constant parameter in this paper for simplicity. An matrix free algorithm is proposed : Firstly, the measurement signal z is assigned to ŷ (ŷ is set to be 0 initially). Then execute steps 2 -6 until ∥ŷ∥ L2 ∥z∥ L2 < (1c). It is noteworthy that the trick is to calculate the Gabor transform C Nw×P = GT (ŷ, Nw) instead of calculating the matrix multiplication cN×1 = Ψ † ŷ, which avoids storing a huge data matrix in memory. In step

Set ŷ = z ; 2: While ∥ŷ∥ L2 ∥ẑ∥ L2 ≥ (1 -c) do 3: Compute C Nw×P = GT (ŷ, Nw) and find c * m * ,l * = max|C|, where 1 < m < P, 1 < l < Nw ; 4: xi = IGT (c * m * ,l * , Nw, m * , l * ) ; 5: ŷ = ŷ -xi ; 6: End While 7: Compute x = z -ŷ ;

Connection with Morphological Component Analysis

In this section, the connection between Morphological Component Analysis(MCA) [START_REF] Starck | Redundant multiscale transforms and their application for morphological component separation[END_REF][15] and the proposed model ( 9) is explained. For the MCA model in Eq. ( 2), there are three essential considerations :

• Signal x can be represented by dictionary Ψ and the coefficients α are sparse. According to the uncertainty principle [START_REF] Elad | Sparse and redundant representations : from theory to applications in signal and image processing[END_REF][18], any two different representations of the same signal using an arbitrary dictionary cannot be jointly sparse. Therefore, the coefficients under the dictionary Ψ for signal y is definitely dense.

For the current problem, the transient signal x is represented by the Gabor dictionary Ψ where the coefficients α are sparse, while noise y can not find a sparse representation under any dictionary. Thus, the coefficients under the dictionary Ψ for noise y is also dense. • Signal y can be represented by dictionary Φ where coefficients β are sparse. Similarly, the coefficients under the dictionary Φ for signal x are definitely dense according to the uncertainty principle. For the current problem, this part (λ ∥β ∥ L 0 ) disappears because the noise y can not find a sparse representation under any dictionary.

• A constant parameter λ is used to balance the ratio between the energy of signal x and the energy of signal y implicitly. The physical meaning is that the weight is assigned to the prior hypothesis and the these two combined prior hypothesis (regularization terms) terms are not equally important. When λ is increased, ∥β ∥ L0 is given more weight, which means that the signal-to-noise ratio is increased. In fact, it is important to note that if λ was assigned to ∥α∥ L0 (formally written as λ ∥α∥ L 0 + ∥β ∥ L 0 ), the opposite conclusion would have been obtained : when λ is increased, ∥α∥ L0 is given more weight, which means that the signal-to-noise ratio is decreased. For the current problem, the additional condition

∥Ψα∥ L 2 ∥y∥ L 2
= c is used to balance the ratio between the energy of the transient signal x and the energy of noise y explicitly. When c is increased, it suggests that the signal-to-noise ratio is increased.

Low Rank + Sparse Combination Model

For the rolling element bearing diagnosis problem [START_REF] Cocconcelli | Stft based approach for ball bearing fault detection in a varying speed motor[END_REF], it is evident that the transient signal is sparse compared to the noise y. Based on this assumption, Eq. ( 2) can be reformulated as follows [START_REF] Wright | Dense error correction via l1-minimization[END_REF] [START_REF] Ganesh | Dense error correction for low-rank matrices via principal component pursuit[END_REF] :

Minimize ∥β ∥ L0 + λ ∥x∥ L0 Subject to z = x + Φβ . ( 11
)
It is noted that the noise y can not be represented parsimoniously under one dictionary Φ, thus, Eq. ( 11) can be further reformulated as follows. Define matrix Y ([Y] m,l denotes the elements of matrix Y) the Gabor transform of noise signal y, that is

[Y] m,l = ŷ[m, l] = GT (y[n]
). The noise y is assumed to be stationary, thus L = E{|Y| 2 ml } can be recognized as a low rank matrix, as a matter of fact Rank(L) = 1 (see proof for more detail). Matrix X is defined as the Gabor transform of transient signal x, which denotes as

[X] m,l = x[m, l] = GT (x[n]
) and it is obvious that S = E{|X| 2 ml } is a sparse matrix. Matrix Z is defined as the Gabor transform of measured signal z, denoted as

[Z] m,l = ẑ[m, l] = GT (z[n]
) and the Spectrogram M = E{|Z| 2 ml } of signal z can be represented as M = L + S . L corresponds to stationary noise signal which is low rank, and S corresponds to the transient signal which is sparse. Thus, the signal separation problem can be solved by :

Minimize Rank(L) + λ ∥S∥ L0 Subject to L + S = M. ( 12
)
The minimization of the Rank and the L 0 norm function are both NP-hard problem, therefore Eq. ( 12) is relaxed into Eq. ( 13) entitled Robust Principal Component Analysis (RPCA) in the literature [START_REF] Emmanuel | Robust principal component analysis ?[END_REF] :

Minimize ∥L∥ * + λ ∥S∥ L1 Subject to L + S = M (13)
where ∥L∥ * := ∑ n k=1 σ k (L) is the nuclear norm of matrix L computed on the singular values σ 1 (L) ≥ σ 2 (L) ≥ ...σ n (L), and ∥S∥ L1 = ∑ i=P, j=Q i=1, j=1 |S i j |. Nevertheless, matrix S can not be directly used to reconstruct the transient signal x because the phase information is lost. However, the latter can be recovered by taking over the phase information of matrix Z, denoted by

[X] m,l = √ [S] m,l e ∠[Z] m,l and [Y] m,l = √ [L] m,l e ∠[Z] m,l .
It is noteworthy to point out that the sparse matrix X is often contaminated by noise when λ is chosen to be small (it suggests that some amount of noise is decomposed as sparse components). Therefore, a filter can be further used to refine the final results which is written as follows :

[F] m,l = { 1 if |[X] m,l | > G|[Y] m,l | 0 if |[X] m,l | ≤ G|[Y] m,l | (14)
where G is a parameter to control the amplitude gain ratio between |[X] m,l | and |[Y] m,l |, similar filter was used in singing-voice separation scenario [START_REF] Huang | Singing-voice separation from monaural recordings using robust principal component analysis[END_REF]. The purpose of this filter is to reduce noise in the sparse matrix X.

An illustrative example is seen in Fig. 2 : Fig. 2 (e) shows a matrix L that is rank-1, however, sparse matrix S is corrupted by noise in Fig. 2 (f). X is also corrupted by noise since the difference between S and X is only that the phase information is added to X by [X] m,l = √ [S] m,l e ∠[Z] m,l as compared to S (the noisy X figure is not shown in this paper). Figure . 3 shows the results after the filter operation F ⊙ X, where ⊙ denotes Hadamard product. The parameter λ in Eq. ( 13) is a trade-off between the low rank matrix and sparse matrix. When λ is decreased , the decomposed matrix L has lower rank and S is less sparse, and vice versa. In Fig. 2 (a)(b), λ is chosen to be 0.4, it is observed that no sparse components are decomposed (S is a zero matrix). In Fig. 2 (c)(d), λ is decreased to be 0.1, it is found that more sparse components are decomposed for matrix S. When λ is further decreased to be 0.01, the rank-1 matrix L is fully decomposed, but the matrix S is no longer sparse due to the noise corruption resulting in the utilization of filter. Finally, the algorithm procedure for separating transient signal from noise is presented as follows :

Algorithm 2 Algorithm Procedure For Low Rank + Sparse Combination Model 1: Calculate spectrogram matrix M of measured signal z by Eq. ( 10).

2: Decompose spectrogram matrix M into low rank real matrix L and sparse real matrix S by Eq. ( 13). 4: Filter the noisy X by operation X = F ⊙ X, where X is the filtered result. 5: Calculate Inverse Gabor Transform (IGT) of sparse complex matrix X to reconstruct transient signal x by Eq. [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF].

It is noteworthy to mention that without the expected value, matrix L in Eq. ( 12) is not a low rank matrix in practice. However, it is found that Eq. ( 13) is still close to the reality since the eigenvalues of matrix L decrease exponentially with the numerical index (more details about this are omitted here for the sake of brevity).

Simulation Experiment

• Simulation setup : a simulation setup to produce one transient signal x that is buried in the strong noise y is built by using the following equation :

s[n] = Nm ∑ m=0 δ (n -mD T ) for n = 1, ...N (15)
H[z] = 1 1 -2cos(2π • 0.37)rz -1 + r 2 z -1 (16)

Conclusion

In this paper, two kinds of combined regularization structures have been investigated to separate transient signals from strong additive noise in rolling element bearings diagnosis. First, the Sparse + Sparse combined It is found that this model can not be directly used in current problem, thus, a modified Sparse + Sparse combined regularization is proposed instead. Second, Low Rank + Sparse combined regularization is investigated and the connection with the Sparse + Sparse combined regularization is elucidated (two models are investigated in a uniform framework). The following topics can be further investigated :

• Is the Gabor dictionary good enough to represent the transient signal ? More specifically, how to choose a better dictionary ? • How to choose the regularization parameter wisely since the separation performance highly depends on it ?

A Proof of Low Rank Property For Matrix L

Let us define the autocovariance R 2y

[n -n ′ ] = E{y[n]y * [n ′ ]} and k = n -n ′ L = [Y 2 ] m,l = E{|Y| 2 ml } = ∑ n ∑ n ′ ψ * ml [n ′ ]ψ ml [n]E{y[n]y * [n ′ ]} = ∑ n ∑ n ′ g N [n -m]g N [n ′ -m]exp{ i2πl(n -n ′ ) N }R 2y [n -n ′ ] = ∑ n ∑ k g N [n -m]g N [n -k -m]exp{ i2πlk N }R 2v [k] = ∑ k ( ∑ n g N [n -m]g N [n -k -m] R NN [k] )exp{ i2πlk N }R 2y [k] = ∑ k R NN [k]exp{ i2πlk N }R 2y [k] = ∑ k R NN [k]R 2y [k]exp{ i2πlk N } (17)
The final results depends only on l. Thus, the matrix L is low rank.

B Accelerated Proximal Gradient Algorithm

In this paper, Accelerated Proximal Gradient (APG) algorithm [START_REF] Lin | Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix[END_REF] is used to solved Eq. ( 13), which is rewritten as a penalized formulation as follows :

Minimize µ(∥L∥ * + λ ∥S∥ L1 ) + 1 2 ∥M -L -S∥ 2 F ( 18
)
where µ is a regularization parameter which is used as a trade-off between the regularization term ∥L∥ * + λ ∥S∥ L1 and the data fitting term 1 2 ∥M -L -S∥ 2 F . It is also observed that the regularization term is a nondifferentiable function and the data fitting term is a smooth function with Lipschitz gradient (Lipschitz constant L f = 2). A soft thresholding operator S λ µ L f for any matrix Q is defined as :

S λ µ L f (Q) = max(|q| - λ µ L f )sgn(q) (19)
where q denotes each entry of the matrix Q and sgn() is the sign function. A soft thresholding D µ L f for any matrix Q is defined as :

D µ L f (Q) = US µ L f (Σ)V T ( 20
)
where matrix Q = UΣV T is decomposed by Singular Values Decomposition (SVD). Therefore, the accelerated proximal gradient algorithm can be iteratively written as algorithm 3. The iteration stops when SC ≤ ε, where the stopping criterion SC and C L , C S are defined as :

SC = ∥[C L ,C S ]∥ F L f max{1, ∥[C L ,C S ]∥ F } C L = L f (G k L -L k+1 ) + (L k+1 + S k+1 -G k L -G k S ) C S = L f (G k S -S k+1 ) + (L k+1 + S k+1 -G k L -G k S ) (21)
where [C L ,C S ] denotes the matrix constructor operator that vertically concatenates C L and C S , and ∥∥ F denotes the Frobenius Norm. For the simulation in this paper, the parameters are chosen fixed as : µ = 10 -3 , η = 0.9 and ε = 10 -7 .

Algorithm 3 Accelerated Proximal Gradient Algorithm Input:

Spectrogram matrix M and parameter λ and µ Output:

Low rank real matrix L and sparse real matrix S 1: Set L 0 = L 1 = 0, S 0 = S 1 = 0, t 0 = t 1 = 1, µ 0 = ∥M∥ F 2: While not converged do 3:

G k L = L k + t k-1 -1 t k (L k -L k-1 ) 4: G k S = S k + t k-1 -1 t k (S k -S k-1 ) 5: L k+1 = D µ k L f [G k L -L f -1 (M -G k L -G k S )] 6: S k+1 = S λ µ k L f [G k S -L f -1 (M -G k L -G k S )] 7: t k+1 = 1 2 (1 + √ 1 + 4t 2 k ) 8: µ k+1 = max(η µ k , µ)

FIGURE 1 -

 1 FIGURE 1 -Transient Signal x + Additive Stationary Noise y = Measured Signal z

3 :

 3 Recover sparse complex matrix X by [X] m,l = √ [S] m,l e ∠[Z] m,l .

Firstly, an impulse

  signal s[n] (time duration is 1s, sampling frequency is 10 4 Hz) with length N = 10 4 is produced by Eq. (15) where D T = 5 × 10 2 and Nm = 19. Secondly, the impulse signal s[n] is passed through a filter H[z] in Eq. (16) with r = 0.98 to synthesize the transient signal x[n]. Finally, the simulated measured signal z[n] is synthesized by adding the transient signal x[n] to Gaussian noise y[n] ∼ N (0, 0.5). • Simulation results for Sparse + Sparse model : in this simulation experiment, the parametric analysis focuses on the window length Nw and signal-to-noise ratio parameter c. Firstly, the Hanning window length is fixed to Nw = 256, and the ratio c is increased from 0.05 and 0.10 to 0.15. It is observed that more and more transient signal components are separated , see Fig. (4)(a)(c)(e). Secondly, the Hanning window length is set to Nw = 1024, and similar phenomenon is also observed, see Fig. (4)(b)(d)(f). It is also observed that the separation results are noisier with the increase of window length from Nw = 256 to Nw = 1024. • Simulation results for Low Rank + Sparse model : Fig. 5 shows the separation results with different parameters by Low Rank + Sparse model. In Fig. 5 (a)(b)(c)(d), window length Nw = 64 and gain G = 6 are set to be constant, and it is observed that more and more transient signal is captured with the decrease of λ from 0.02 to 0.001. It is noted that Fig. 5 (a) fails to capture the transient signal since too many sparse components in X are recovered (λ is set too small) , which can be verified by Fig. 3 (a).

Figure 5 (

 5 e)(f) show the results when the widow length is increased to Nw = 256.

FIGURE 2 -FIGURE 3 -FIGURE 4 -FIGURE 5 -

 2345 FIGURE 2 -Spectrogram M is decomposed into a low rank matrix L and a sparse matrix S with different parameters λ : (a) Low rank decomposition L with λ = 0.4, (b) Sparse decomposition S with λ = 0.4, (c) Low rank decomposition L with λ = 0.1, (d) Sparse decomposition S with λ = 0.1, (e) Low rank decomposition L with λ = 0.01, (f) Sparse decomposition S with λ = 0.01,

  3, c * m * ,l * is the element with the highest value in matrix C, m * and l * denote the position index in the matrix C, which is formally written as c * m

* ,l * = max|C| (|C| denotes absolute value for each entry in matrix C). The calculation xi = d i c * m * ,l * (d i denote the i-th atom in dictionary Ψ) is implemented by inverse Gabor transform xi = IGT (c * m * ,l * , Nw, m * , l * ) of the corresponding time sample m * and frequency l * . Algorithm 1 Matrix Free Matching Pursuit Algorithm Input: Measured signal z ; Signal-to-noise ratio c ; Window length Nw ; Output: Estimated transient signal x ; 1:


André L F De Almeida
email: andre.almeida@ieee.org

Alain Y Kibangou
email: alain.kibangou@ujf-grenoble.fr

DISTRIBUTED LARGE-SCALE TENSOR DECOMPOSITION

Keywords: Tensor decompositions, large-scale data, distributed computation

Canonical Polyadic Decomposition (CPD), also known as PARAFAC, is a useful tool for tensor factorization. It has found application in several domains including signal processing and data mining. With the deluge of data faced in our societies, large-scale matrix and tensor factorizations become a crucial issue. Few works have been devoted to large-scale tensor factorizations. In this paper, we introduce a fully distributed method to compute the CPD of a large-scale data tensor across a network of machines with limited computation resources. The proposed approach is based on collaboration between the machines in the network across the three modes of the data tensor. Such a multi-modal collaboration allows an essentially unique reconstruction of the factor matrices in an efficient way. We provide an analysis of the computation and communication cost of the proposed scheme and address the problem of minimizing communication costs while maximizing the use of available computation resources.

INTRODUCTION

From Internet to large research infrastructures, the volume of data generated by our societies is continuously increasing. A deluge faced by the producers of these data as well as their users. The big data issue is a significant scientific challenge that requires deep investigations in both engineering and fundamental science. Everyone is concerned and it is urgent to get answers to questions such as how to store these huge amount of data? How to process and analyze them?

Low-rank matrix factorization has received a particular attention in recent years, since it is fundamental to a variety of mining tasks that are increasingly being applied to massive datasets. In large applications, matrix factorizations can involve matrices with billions of entries. At this massive scale, distributed algorithms for matrix factorization are essential to achieve reasonable performance [START_REF] Gemulla | Large-Scale Matrix Factorization with Distributed Stochastic Gradient Descent[END_REF]. However, in many disciplines, data inherently has more than two axes of variation and can be arranged as tensors (i.e. multi-way arrays). Computing tensor decompositions of multi-way datasets is particularly useful to extract hidden patterns and structure in data analytics problems. Specifically, CPD (Canonical Polyadic Decomposition) also known as PARAFAC (Parallel factor analysis) is an extension of a low rank matrix decomposition to tensors.

In order to compute CPD, several algorithms have been proposed in the literature, which can be classified into three main categories: alternating algorithms, derivative based algorithms, and non-iterative algorithms (see e.g. [START_REF] Tomasi | A comparison of algorithms for fitting the PARAFAC model[END_REF][START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF][START_REF] Tichavsky | Simultaneous Search for all Modes in Multilinear Models[END_REF]). Recently, considering that the data tensor can be serially acquired or the underlying process can be time-varying, adaptive algorithms have been proposed in [START_REF] Nion | Adaptive algorithms to track the PARAFAC decomposition of a third-order tensor[END_REF]. And more recently, assuming that the data tensor is spatially spread out, distributed algorithms have been introduced in [START_REF] Kibangou | Distributed PARAFAC based DS-CDMA blind receiver for wireless sensor networks[END_REF] and [START_REF] Kibangou | Distributed Khatri-Rao space-time coding and decoding for cooperative networks[END_REF]. These algorithms are based on an alternating least square (ALS) algorithm that proceeds iteratively and minimizes a criterion (that is usually quadratic) with respect to individual factors one by one. Then, distributed algorithms are obtained using the concept of average consensus that has been extensively studied in computer science and is a central topic for load balancing (with divisible tasks) in parallel computers [START_REF] Xiao | Distributed average consensus with least-mean-square deviation[END_REF].

The need for large-scale tensor computations is increasing and there is a huge gap to be filled. In contrast to large-scale matrix factorization, very few works are devoted to large-scale tensors. Two ways have been recently considered in the literature. The first one consists in exploiting sparseness of tensors. In [START_REF] Kang | GigaTensor: scaling tensor analysis up by 100 times-Algorithms and discoveries[END_REF], the GigaTensor algorithm is designed in order to minimize the number of floating point operations and to handle the size of intermediate data in order to overcome the intermediate data explosion problem. For large-scale tensors the intermediate data explosion problem arises when Khatri-Rao matrix products are implemented in a naive way. In [START_REF] Kang | GigaTensor: scaling tensor analysis up by 100 times-Algorithms and discoveries[END_REF] a smart ordering of computations is proposed and explicit computation of Khatri-Rao matrix product is avoided. In addition, a way to implement the proposed scheme in the distributed computing frameworks MapReduce and Hadoop is proposed by the authors and then, from extensive simulations, linear scalability on the number of machines is claimed. Note that the data explosion problem can also be handled by exploiting sparseness of both data and latent factors as in [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF][START_REF] Acar | Scalable tensor factorizations for incomplete data[END_REF][START_REF] Sidiropoulos | Multi-way compressed sensing for sparse low-rank tensors[END_REF].

The second class of methods consists in a subdivision of the large-scale tensor into smaller ones. Then factor matrices are reconstituted from the estimated sub-factors.

The work in [START_REF] Papalexakis | ParCube: Sparse parallelizable tensor decompositions[END_REF] is motivated by the success of random sampling-based matrix algorithms. The large-scale tensor is under-sampled several times, then the different sub-tensors are processed in parallel and eventually the results are combined in a clever way. In order to overcome permutation indeterminacies, all the different sub-tensors are enforced to overlap in a common set of indices in all the three modes of the tensor. The merging operation is guaranteed to be successful only if the random samples (sub-tensors) fulfill CPD identifiability conditions. In [START_REF] Phan | PARAFAC algorithms for large-scale problems[END_REF], sub-division of the original tensor is achieved in a deterministic way; a grid decomposition is proposed. CPD is first performed for each sub-tensor in order to get a an alternative representation that allows avoiding Khatri-Rao matrix products. Then, grid-CPD performs as a decentralized computation method where at each step parallel machines have to communicate with a central server. There is no direct collaboration between the involved machines. In [START_REF] De Almeida | Distributed Computation of Tensor Decompositions in Collaborative Networks[END_REF], the authors have introduced the notion of collaboration in the grid. In this set-up, the computation is really distributed. The central server only initiates the process by sending data to the machines in the network and merges the results at the end of the computation process. During the computation process, the central server can allocate its own resources to other tasks. Moreover, thanks to this full collaboration and by introducing overlapping between sub-tensors, uniqueness of the decomposition with indeterminacies as in the standard CPD can be guaranteed. However, collaboration was restricted to two modes and overlapping induced loss of efficiency (increase of the dimensions and the number of sub-tensors).

This paper introduces a new way for computing a grid-CPD of a large-scale data tensor that ensures essential uniqueness and perfect reconstruction of the factor matrices. The main idea of the proposed approach is to allow collaboration between the machines in the network across the three modes of the tensor. To cope with the multi-way nature of the data, multi-graphs are used to model the network. Due to direct collaboration between machines, fulfilling CPD uniqueness conditions is not required for each sub-tensor. We provide an analysis of the computation and communication cost of the proposed scheme and address the problem of minimizing communication costs while maximizing the use of available computation resources. Notations: Vectors are written as boldface lower-case letters (a,b,• • • ), matrices as boldface upper-case letters (A,B,• • • ), and tensor as calligraphic letters (X , Y, • • • ). A T stands for the transpose of A whereas A H stands for its complex conjugate. The operator • denotes the outer product between vectors, while ⊙ stands for the Khatri-Rao (columnwise Kronecker) product. The Hadamard (element-wise) product is denoted by the * symbol. For a three-way tensor X ∈ C I×J ×K , I, J and K are referred to as the mode-1, mode-2 and mode-3 dimensions of X , respectively. The operator ⊔1 denotes the mode-1 concatenation of any two tensors X and Y having the same mode-2 and mode-3 dimensions. This operation is similarly defined to denote concatenation along the other modes.

PROBLEM SETTING

Consider a large-scale tensor X ∈ C I×J ×K , i.e. a tensor with a total number of entries IJK that can easily be of order of billions. The CPD of X , denoted by X = [[A, B, C]], is given by [START_REF] Harshman | Foundation of the PARAFAC procedure: models and conditions for an "explanatory" multimodal factor analysis[END_REF][START_REF] Caroll | Analysis of individual differences in multidimensional scaling via an N-way generalization of "Eckart-Young" decomposition[END_REF]:

X = R r=1 ar • br • cr, (1)
where

A = [a1, . . . , aR] ∈ C I×R , B = [b1, . . . , bR] ∈ C J ×R and C = [c1, . . . , cR] ∈ C K×R
are the factor matrices associated with the CPD of X , and R denotes the tensor rank. A sufficient condition for the essential uniqueness of the CPD was first established by Kruskal in [START_REF]Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics[END_REF], and states that A, B and C can be uniquely estimated (up to column permutation and scaling) if kA+kB+kC ≥ 2R + 2, where k (•) denotes the Kruskal-rank, or "k-rank", of its matrix argument. Shortly, k (A) is equal to r if every set of r columns of A is linearly independent. In what follows, we work under the assumption that CPD uniqueness for X holds. Usually, the matrices A, B, and C can be estimated by computing an approximation to the rank-R CPD of X , i.e.:

min A,B,C X - R r=1 ar • br • cr 2 F . (2)
The traditional way for carrying out such a computation is to resort to an alternating least squares algorithm (ALS) where estimations of A, B, and C are obtained by alternatively solving three different linear least squares (LS) problems. Despite its conceptual simplicity, ALS involved Khatri-Rao matrix products and multiplications of matrices that cannot be handled when the tensor exhibits large dimensions since computation and storage resource can easily become insufficient. Apart from complexity issues, ALS-based algorithms present a very slow convergence due to the huge number of unknowns involved in the estimation of the large factor matrices A, B and C. Consider a central server having at its disposal a large-scale data tensor X to be factorized. This server is connected to a dense network of L machines (cores or independent computers) with limited processing powers. Our goal is to perform the CPD of X across the network by using all the in-network computation resources in an optimal way while making the central server available to other tasks. The central server and the L machines, where computations are achieved, are nodes of a network connected through reliable high speed data links, so that data exchanges may be considered noise-free. The central server first subdivides the tensor and assigns each sub-tensor to a given machine. Once computation is achieved, the L machines send their results back to the central server where global factor matrices are built. In contrast to the grid-CPD in [START_REF] Phan | PARAFAC algorithms for large-scale problems[END_REF], the central server is not involved in the computation process and available to perform other tasks at the same time. The proposed set-up is therefore fully distributed while that of [START_REF] Phan | PARAFAC algorithms for large-scale problems[END_REF] is simply decentralized.

PROPOSED SOLUTION

The central server generates different (possibly overlapped) data sub-tensors

X (ℓ 1 ,ℓ 2 ,ℓ 3 ) ∈ C I ℓ 1 ×J ℓ 2 ×K ℓ 3 , ℓ1 = 1, . . . , L1, ℓ2 = 1, . . . , L2, ℓ3 = 1, . . . , L3
, of much smaller dimensions, by partitioning X in a dense 3-D tensor grid. The dimensions of the sub-tensors are related as follows: I1 + . . . + IL 1 ≥ I, J1 +. . .+JL 2 ≥ J, K1 +. . .+KL 3 ≥ K. Equality in a given mode occurs if and only if the sub-tensors do not overlap in that mode. The sub-tensors X (ℓ 1 ,ℓ 2 ,ℓ 3 ) can be concatenated along different pair of modes to form mode-1, mode-2, and mode-3 sub-tensors

X (ℓ 1 ) 1 ∈ C I ℓ 1 ×J ×K , X (ℓ 2 ) 2 ∈ C I×J ℓ 2 ×K , X (ℓ 3 ) 3
∈ C I×J ×K ℓ 3 , which are defined, respectively, as

X (ℓ 1 ) 1 = X (ℓ 1 ,1) 1 ⊔3 • • • ⊔3 X (ℓ 1 ,L 3 ) 1 (3)
where

X (ℓ 1 ,ℓ 3 ) 1 = X (ℓ 1 ,1,ℓ 3 ) ⊔2 • • • ⊔2 X (ℓ 1 ,L 2 ,ℓ 3 ) , X (ℓ 2 ) 2 = X (ℓ 2 ,1) 2 ⊔1 • • • ⊔1 X (ℓ 2 ,L 1 ) 2 (4)
where

X (ℓ 2 ,ℓ 1 ) 2 = X (ℓ 1 ,ℓ 2 ,1) ⊔3 • • • ⊔3 X (ℓ 1 ,ℓ 2 ,L 3 ) , X (ℓ 3 ) 3 = X (ℓ 3 ,1) 3 ⊔2 • • • ⊔2 X (ℓ 3 ,L 2 ) 3 (5)
where

X (ℓ 3 ,ℓ 2 ) 3 = X (1,ℓ 2 ,ℓ 3 ) ⊔1 • • • ⊔1 X (L 1 ,ℓ 2 ,ℓ 3 ) .
An interesting fact resulting from the CPD trilinearity is as follows: if the CPD of X is essentially unique and provided that k A (ℓ 1 ) = kA, k B (ℓ 2 ) = kB, and k C (ℓ 3 ) = kC, each mode-i sub-tensor, i = 1, 2, 3, admits an exact CPD given by:

X (ℓ 1 ) 1 = [[A (ℓ 1 ) , B, C]], X (ℓ 2 ) 2 = [[A, B (ℓ 2 ) , C]] X (ℓ 3 ) 3 = [[A, B, C (ℓ 3 ) ]]. (6)
Now, the central server assigns a sub-tensor X (ℓ 1 ,ℓ 2 ,ℓ 3 ) to each of the L available machines. To guarantee a one-to-one mapping between the sub-tensors and the machines, each machine is uniquely labeled by the triad (ℓ1, ℓ2, ℓ3) and is associated with the corresponding sub-tensor X (ℓ 1 ,ℓ 2 ,ℓ 3 ) . Therefore the total number L of machines is equal to L1L2L3. According to [START_REF] Nion | Adaptive algorithms to track the PARAFAC decomposition of a third-order tensor[END_REF], from their respective sub-tensors, the machines having the same i-th coordinate, i = 1, 2, 3, can collaborate to estimate their common sub-factor matrix in mode-i. Therefore the adjacency relations between the machines evolve according to the mode of interest. To capture this behavior, we adopt multi-layer graphs, or simply multi-graphs, for representing the network [START_REF] Tang | Community detection via heterogeneous interaction analysis[END_REF][START_REF] Dong | Clustering with Multi-Layer Graphs: A Spectral Perspective[END_REF]. For this purpose, let G be a multi-layer graph which contains 3 individual graph layers G (i) , i = 1, 2, 3, where each layer G (i) is a undirected graph consisting of a common vertex set V , with cardinality L, and a specific edge set E (i) . Each layer is associated with a mode of the tensor to be decomposed. In addition, each layer is constituted with Li connected components of L/Li vertices, each connected component being a star graph 1 . More precisely, the nodes labeled (l1, l2, l3) and (l ′ 1 , l ′ 2 , l ′ 3 ) belongs to the same component of G (i) if and only if li = l ′ i . They are adjacent only if one of them is an internal node and the other is a leaf of a tree. Therefore, each layer can be viewed as a set of parallel star sub-networks, each star sub-network constituting a connected component of the graph. As it is common for tensor decomposition, data processing is to be carried out in an alternating way. The layers of the multi-graph are to be considered one after the other. In time, the network topology switches periodically between parallel star subnetworks (see Fig. 1).

X (1,1,1) X (1,1,L3) X (L1,1,L3) X (L1,L2,L3) X (1,L2,L3) X (L1,L2 ,1)
X (1,L2,1)

X (L1,1,1) (l1, 1, 1) (l1, l2, l3) (l1, L2 -1, l3 -1) (l1, L2, L3 -1)
(l1, L2, L3)

Fig. 1.

Left: Subdivision of a large-scale tensor in small sub-tensors X (l 1 ,l 2 ,l 3 each one being associated to a machine labeled (ℓ1, ℓ2, ℓ3). Right: ℓ1-th star subnetwork of layer 1 associated to the tensor mode 1.

By exploiting the degrees of freedom provided by the three layers G (i) , i = 1, 2, 3, of the multi-graph G, the problem consists in finding A (ℓ 1 ) , ℓ1 = 1, . . . , L1, B (ℓ 2 ) , ℓ2 = 1, . . . , L2, and C (ℓ 3 ) , ℓ3 = 1, . . . , L3 by solving the three following sets of LS problems: min

A (ℓ 1 ) X (ℓ 1 ) 1 -[[A (ℓ 1 ) , B, C]] 2 F , ℓ1 = 1, . . . , L1, (7) min B (ℓ 2 ) X (ℓ 2 ) 2 -[[A, B (ℓ 2 ) , C]] 2 F , ℓ2 = 1, . . . , L2, (8) min C (ℓ 3 ) X (ℓ 3 ) 3 -[[A, B, C (ℓ 3 ) ]] 2 F , ℓ3 = 1, . . . , L3. (9)
Each set of LS problems is associated with disjoint subsets of machines and therefore can be solved in parallel. For instance, the L1 LS problems in mode-1 are solved independently by L1 subsets of machines working in parallel. The same applies to the two other 1 The choice of a star topology is not restrictive. According to the available resource any other kind of topology can be adopted. LS problems. Moreover, fast computations can be performed on small sub-tensors at each machine. Let us define X

(ℓ 1 ,ℓ 2 ,ℓ 3 ) (1) ∈ C J ℓ 2 K ℓ 3 ×I ℓ 1 , X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (2) ∈ C K ℓ 3 I ℓ 1 ×J ℓ 2 , and X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (3)
∈ C I ℓ 1 J ℓ 2 ×K ℓ 3 as the matrices obtained by unfolding ("matricization" of) the sub-tensors X

(ℓ 1 ) 1 , X (ℓ 2 ) 2
and X (ℓ 3 ) 3 along mode-1, mode-2 and mode-3, respectively. Note that these matrices collect different rearrangements of the same data contained in the sub-tensor X (ℓ 1 ,ℓ 2 ,ℓ 3 ) . With these definitions, we can recast each one of the LS problems ( 7)-( 9) as summations over smaller LS subproblems, as follows:

J A (ℓ 1 ) = min A (ℓ 1 ) L 2 ℓ 2 =1 L 3 ℓ 3 =1 X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (1) -(B (ℓ 2 ) ⊙ C (ℓ 3 ) )A (ℓ 1 )T 2 F J (ℓ 2 ,ℓ 3 ) A (ℓ 1 )
,

J B (ℓ 2 ) = min B (ℓ 2 ) L 1 ℓ 1 =1 L 3 ℓ 3 =1 X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (2) -(C (ℓ 3 ) ⊙ A (ℓ 1 ) )B (ℓ 2 )T 2 F J (ℓ 1 ,ℓ 3 ) B (ℓ 2 )
,

J C (ℓ 3 ) = min C (ℓ 3 ) L 1 ℓ 1 =1 L 2 ℓ 2 =1 X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (3) -(A (ℓ 1 ) ⊙ B (ℓ 2 ) )C (ℓ 3 )T 2 F J (ℓ 1 ,ℓ 2 ) C (ℓ 3 )
.

The solution that minimizes J A (ℓ 1 ) (resp. J B (ℓ 2 ) and J A (ℓ 3 ) ) can then be written as:

A (ℓ 1 )T =   1 L2L3 ℓ 2 ,ℓ 3 Γ (ℓ 2 ,ℓ 3 ) A (ℓ 1 )   -1   1 L2L3 ℓ 2 ,ℓ 3 Ψ (ℓ 2 ,ℓ 3 ) A (ℓ 1 )   = Γ A (ℓ 1 ) -1 Ψ A (ℓ 1 ) (10)
B (ℓ 2 )T =   1 L1L3 ℓ 1 ,ℓ 3 Γ (ℓ 1 ,ℓ 3 ) B (ℓ 2 )   -1   1 L1L3 ℓ 1 ,ℓ 3 Ψ (ℓ 1 ,ℓ 3 ) B (ℓ 2 )   = Γ B (ℓ 2 ) -1 Ψ B (ℓ 2 ) (11)
C (ℓ 3 )T =   1 L1L2 ℓ 1 ,ℓ 2 Γ (ℓ 1 ,ℓ 2 ) C (ℓ 3 )   -1   1 L1L2 ℓ 1 ,ℓ 2 Ψ (ℓ 1 ,ℓ 2 ) C (ℓ 3 )   = Γ C (ℓ 3 ) -1 Ψ C (ℓ 3 ) (12)
where Γ (ℓ 2 ,ℓ 3 )

A (ℓ 1 ) = B (ℓ 2 )H B (ℓ 2 ) * C (ℓ 3 )H C (ℓ 3 ) ∈ C R×R , Γ (ℓ 1 ,ℓ 3 ) B (ℓ 2 ) = C (ℓ 3 )H C (ℓ 3 )H * A (ℓ 1 )H A (ℓ 1 ) ∈ C R×R , Γ (ℓ 1 ,ℓ 2 ) C (ℓ 3 ) = A (ℓ 1 )H A (ℓ 1 )H * B (ℓ 2 )H B (ℓ 2 ) ∈ C R×R , Ψ (ℓ 2 ,ℓ 3 ) A (ℓ 1 ) = (B (ℓ 2 ) ⊙ C (ℓ 3 ) ) H X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (1) ∈ C R×I ℓ 1 , Ψ (ℓ 1 ,ℓ 3 ) B (ℓ 2 ) = (C (ℓ 3 ) ⊙ A (ℓ 1 ) ) H X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (2) ∈ C R×J ℓ 2 , and Ψ (ℓ 1 ,ℓ 2 ) C (ℓ 3 ) = (A (ℓ 1 ) ⊙ B (ℓ 2 ) ) H X (ℓ 1 ,ℓ 2 ,ℓ 3 ) (3) ∈ C R×K ℓ 3 .
The estimation of {A (1) , . . . , A (L 1 ) }, {B (1) , . . . , B (L 2 ) }, and {C (1) , . . . , C (L 3 ) } is carried out in three steps. First, the L1 connected components of G (1) operate in parallel to estimate A (1) , . . . , A (L 1 ) . Second, the L2 connected components of G (2) operate in parallel to estimate B (1) , . . . , B (L 2 ) . Finally, the L3 connected components of G (3) operate in parallel to estimate C (1) , . . . , C (L 3 ) . Thanks to the star topology adopted herein, the averaging operations involved in [START_REF] Kang | GigaTensor: scaling tensor analysis up by 100 times-Algorithms and discoveries[END_REF], [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF], and ( 12) can be carried out in two steps 2 . First, the leaves of the star sub-network send their data to the internal node of the corresponding sub-network and then the averaged quantities are sent back to the leaves.

Identifiability issues:

Strictly local computations of a CPD of the (ℓ1, ℓ2, ℓ3)-th node requires satisfying the following necessary condition for the identifiability of the triplet {A (ℓ 1 ) , B (ℓ 2 ) , C (ℓ 3 ) } at each node [START_REF] Liu | Cramér-Rao lower bounds for low-rank decomposition of multidimensional arrays[END_REF]:

min(J ℓ 2 K ℓ 3 , I ℓ 1 K ℓ 3 , I ℓ 1 J ℓ 2 ) ≥ R. (13)
This condition may not always hold, especially for small sub-tensors where the dimensions are too small compared to the number R of hidden factors to be identified from the global data tensor. By allowing collaboration along the three modes, identifiability conditions are improved by relaxing the necessary constraints on the sub-tensor dimensions at each node. In this case, the necessary condition turns out to be min

ℓ 2 ,ℓ 3 J ℓ 2 K ℓ 3 , ℓ 1 ,ℓ 3 I ℓ 1 K ℓ 3 , ℓ 1 ,ℓ 2 I ℓ 1 J ℓ 2 ≥ R, (14)
which is clearly less restrictive than the previous one. Moreover, thanks to the multi-mode cooperation between the machines in the proposed scheme, permutation and scaling indeterminacies are similar to those occurring in the standard CPD. Thus, the reconstruction of the global factor matrices at the central server can be done without resorting to additional structural constraints such as anchor [START_REF] Papalexakis | ParCube: Sparse parallelizable tensor decompositions[END_REF] or overlapping [START_REF] De Almeida | Distributed Computation of Tensor Decompositions in Collaborative Networks[END_REF] sub-tensors.

Communication and computation costs:

The overall computation cost of the proposed scheme is similar to that of the centralized scheme (standard ALS applied to the large-scale tensor). However, the cost per machine can be significantly lowered. For instance, for a large cubic tensor (I = J = K, I > 10 2 ) and considering a uniform sampling of the tensor so that I l = J l = K l = mI, with m < 1, the communication and computation costs are given in Table 1: Now let us define by q the ratio between the power required for achieving one floating operation and that for transmitting one real value. Assuming that q << I 2 , the ratio between the per machine consumptions E d and Ec in the distributed and the centralized schemes, respectively, is given by: E d Ec ≈ (1 + q 2I 2 )m ≈ m < 1. This means that the reduction of power consumption per machine follows the reduction of the tensor dimensions. From Table 1, we can note that the additional power consumption is due to communications. Indeed at each ALS iteration, each machine (ℓ1, ℓ2, ℓ3) has to transmit three R × R matrices (Γ A (ℓ 1 ) , Γ B (ℓ 2 ) , Γ C (ℓ 3 ) ) and three rectangular matrices (Ψ

A (ℓ 1 ) ∈ C R×I ℓ 1 ,Ψ B (ℓ 2 ) ∈ C R×J ℓ 2 Ψ C (ℓ 3 ) ∈ C R×K ℓ 3
). The overall communication cost per iteration is then given by:

C = 3R 2 L + R ℓ 1 ,ℓ 2 ,ℓ 3 (I ℓ 1 + J ℓ 2 + K ℓ 3 ).
2 For arbitrary topologies, exact averaged quantities can be obtained using the finite-time average consensus protocol designed in [START_REF] Kibangou | Graph Laplacian based Matrix Design for Finite-time Distributed Average Consensus[END_REF]. In particular, it has been shown that consensus can be achieved in two steps when the topology is restricted to be strongly regular [START_REF] Kibangou | Step-Size Sequence Design for finite-time Average Consensus in Secure Wireless Sensor Networks[END_REF].

We can conclude that the communication cost grows linearly with the number of machines, meaning that the proposed scheme is scalable in terms of communication cost.

Moreover, by appropriately selecting the dimensions of the sub-tensors, one can significantly reduce the communication cost per node.

Optimizing the sub-tensors' dimensions: Given a number L of available machines, one have to minimize the communication cost per machine so that the necessary identifiability conditions are fulfilled. The dimensions of the sub-tensors and the parameters Li are obtained by solving the following integer programming problem:

minimize 1 L 1 L 2 L 3 ℓ 1 ,ℓ 2 ,ℓ 3 (I ℓ 1 + J ℓ 2 + K ℓ 3 ) s.t. ℓ 1 I ℓ 1 ≥ I, ℓ 2 J ℓ 2 ≥ J, ℓ 3 K ℓ 3 ≥ K, L1L2L3 ≤ L ℓ 2 ,ℓ 3 J ℓ 2 K ℓ 3 ≥ R, ℓ 1 ,ℓ 3 I ℓ 1 K ℓ 3 ≥ R, ℓ 1 ,ℓ 2 I ℓ 1 J ℓ 2 ≥ R.
This problem is NP-hard. A more tractable solution can be obtained if we adopt a uniform sampling for each mode, i.e. I ℓ 1 = m1I, J ℓ 2 = m2J, and K ℓ 3 = m3K. Table 2 gives the optimal configuration for a cubic tensor obtained by solving the above problem relaxed to the case of uniform sampling for each mode. Note that as the number of available machines increases a more efficient use of the resources is obtained and the communication cost per node is reduced, corroborating the relevance of optimizing the sampling of the data tensor across the three modes.

CONCLUSION

Performing CPD for large-scale tensors is particularly challenging. In this paper, we have proposed a fully distributed method that allows computing CPD across a network of machines with limited computation resources. In the proposed scheme, the central server does not interfere in the computation process and can achieve other tasks at the same time. The machines in the network collaborate with each other according to the three modes of the tensor. Such a multi-modal collaboration is captured by a multi-graph whose layers have several connected components. Due to this full collaboration, improved necessary identifiability conditions are achieved and the global factor matrices can be reconstructed in an efficient way. However, a communication overhead is introduced. Hopefully, the communication cost grows linearly with the number of machines and the dimensions of the sub-tensors. Additionally, optimizing the sampling of the tensor across the three modes yields a significant reduction of this overhead while maximizing the use of the available computation resources. Extension to tensors of higher orders is straightforward although the proposed optimization problem becomes more challenging. Future works include resources optimization for heterogeneous networks of machines having different computation capabilities, extension to other tensor factorization, and exploitation of sparsity and structural properties of the factor matrices.

Table 1 .

 1 Complexity evaluation per machine and per iteration for an I × I × I tensor of rank R.

	Method	Computation cost Communication cost
	Centralized	6RI 3	-
	Distributed	6mRI 3	3R 2 + 3mRI

Table 2 .

 2 Optimal number and dimensions of sub-tensors for different numbers of available machines for a 10 3 ×10 3 ×10 3 tensor ( Cm= communication cost per machine; R= rate of machines used).

	♯ machines	Sub-tensor size	(L1, L2, L3)	Cm	R (%)
	20	500 × 334 × 334	(2, 3, 3)	1468	90
	40	334 × 334 × 250	(3, 3, 4)	1218	90
	60	334 × 250 × 200	(3, 4, 5)	1084	100
	80	250 × 250 × 200	(4, 4, 5)	1000	100
	100	250 × 200 × 200	(4, 5, 5)	950	100

This work is partially supported by the CNPq and the FUNCAP-INRIA project TeMP.


Alexandre Bilger

Christian Duriez

Stéphane Cotin

Computation and Visualization of Risk Assessment in Deep Brain Stimulation Planning

Keywords: Deep Brain Stimulation, risk planning, brain shift, simulation

Deep Brain Stimulation is a neurosurgical approach for the treatment of pathologies such as Parkinson's disease. The basic principle consists in placing a thin electrode in a deep part of the brain. To safely reach the target of interest, careful planning must be performed to ensure that no vital structure (e.g. blood vessel) will be damaged during the insertion of the electrode. Currently this planning phase is done without considering the brain shift, which occurs during the surgery once the skull is open, leading to increased risks of complications. In this paper, we propose a method to compute the motion of anatomical structures induced by the brain shift. This computation is based on a biomechanical model of the brain and the cerebro-spinal fluid. We then visualize in a intuitive way the risk of damaging vital structures with the electrode.

Introduction

Deep Brain Stimulation is a neurosurgical technique to treat medication-resistant motion disorders symptoms, such as Parkinson's disease and essential tremor, or affective symptoms such as major depression. It consists in electrically stimulating a specific structure in the deep brain tissue. The surgical procedure involves the implantation of one (unilateral) or two (bilateral) stimulating electrodes. The electrodes are then left in the brain, attached to the skull and linked to a neurostimulator, which sends the electrical impulses.

Before the operation, the surgeon needs to determine the target coordinates in the patient frame, as well as a trajectory to reach the target. The path is linear and has to follow some surgical constraints, such as avoiding vital structures (e.g. blood vessels and ventricles), since the dissection of these structures could cause severe surgical complications. Currently the planning phase relies entirely on pre-operative data (which are a combination of the patient pre-operative image and an atlas containing detailed anatomical information). However, intra-operative brain deformation, called brain shift, can alter the pre-operative planning, since vital structures may have moved in the path of the electrode. This risk is limited by selecting a trajectory not too close from vital structures, in case of brain shift. In practice, this is done by defining a circle (typically with a radius Figure 1. Considering a 5mm safety margin (in white) around an electrode path in classic pre-operative planning, the close blood vessel risks to be damaged if it shifts toward the trajectory of 5mm), centered along the planned trajectory. This safety margin also includes errors due to the fusion and registration methods used during the planning (see Fig. 1). This geometric approach does not correctly describe the complexity of tissue motion during the brain shift. In particular, the brain shift induces an anisotropic tissue motion, while describing the safety margin as a circle assumes an constant isotropic motion.

We propose an advanced, more accurate method which better accounts for the complexity of the brain shift when selecting a safe trajectory. It consists in visualizing a volume around a candidate trajectory, which includes information about the brain tissue motion, based on a biomechanical simulation of brain shift. Our intention is to improve the current planning procedure without changing the surgeon's habit. In our method, the surgeon will check if a vessel is inside a complex closed shape, rather than a simple circle. In the physics-based approach, it means the vessel will shift towards the trajectory, and can potentially be damaged during the electrode insertion.

To achieve this goal, we propose a framework that reproduces the intra-operative brain tissue deformation. It includes a physics-based model of the brain, mechanicals interactions between the brain and the skull, and the influence of the cerebro-spinal fluid (CSF). The framework is based on our previous work [START_REF] Bilger | Biomechanical simulation of electrode migration for deep brain stimulation[END_REF]. Similar works modeled the brain shift phenomena, but a few included the interaction with the CSF [START_REF] Clatz | Patient-specific biomechanical model of the brain: application to Parkinson's disease procedure[END_REF]. Other works such as [START_REF] Wittek | Brain shift computation using a fully nonlinear biomechanical model[END_REF] use external virtual forces to register intra-operative images. A precise prediction of brain shift is difficult as some parameters (e.g. CSF loss volume) are unknown at the time of the planning. That is why most of the papers present methods to model the brain shift with intra-operative data. Our method does not predict the exact brain shift, but prevents risky electrode trajectories with a risk volume notion. The following section shows how we model the brain shift, then we present an intuitive way to visualize a risk to intersect vital structure with the electrode.

Methods & Materials

Numerical Simulation of Brain Shift

Our simulation relies on a physics-based model of the brain tissue deformation, the contact response with the skull and the falx cerebri, and the interaction with the CSF. The brain deformation is computed using a non-linear geometric finite element method model, with a linear constitutive law [START_REF] Felippa | A unified formulation of small-strain corotational finite elements: I. Theory[END_REF]. The model is adapted for large displacements and small deformations, and is fast to compute for a clinical use. Although this model is linear, compared to more complex laws [START_REF] Wittek | Algorithms for computational biomechanics of the brain[END_REF], it is a good trade-off between precision and computation time. Both brain hemispheres are meshed with approximatively 2500 linear tetrahedral elements. Because Young's modulus ranges from 0.6 kPa [START_REF] Clatz | Robust nonrigid registration to capture brain shift from intraoperative MRI[END_REF] to 18 kPa [START_REF] S K Kyriacou | Nonlinear elastic registration of brain images with tumor pathology using a biomechanical model[END_REF] in the literature, we use a mean value of 6000 Pa. Brain is a nearly incompressible medium, therefore Poisson's ratio is set to 0.45. The simulation computes the contacts between the brain, the inner part of the skull and the falx cerebri. They are detected using Layered Depth Images [START_REF] Heidelberger | Detection of collisions and self-collisions using image-space techniques[END_REF] and controlled using constraints solved with Lagrange Multipliers [START_REF] Allard | Volume contact constraints at arbitrary resolution[END_REF]. A null displacement constraint is set on both hemispheres near the area of the brainsterm. Figure 2 illustrates the different constraints and interactions in our simulation.

The main cause of brain shift is a CSF loss, happening after the skull opening [START_REF] Slotty | The impact of brain shift in deep brain stimulation surgery: observation and obviation[END_REF]. In rest position, the fluid forces and the gravity force are balanced. But after the skull opening, fluid leaks and fluid forces decrease. As the brain is more influenced by its weight than by fluid forces, this results in a brain shift. In order to account for this important phenomenon, we include a model of the CSF in our framework. Despite the development of computational fluid dynamics adapted for the CSF (e.g. [START_REF] Kurtcuoglu | Computational fluid dynamics for the assessment of cerebrospinal fluid flow and its coupling with cerebral blood flow[END_REF]), a force model is sufficient in a brain shift simulation context, that is why we use the following model of fluid forces [START_REF] Clatz | Patient-specific biomechanical model of the brain: application to Parkinson's disease procedure[END_REF] that apply on the brain surface:

f CSF = S ρ g h(P) dS
where ρ is the density of CSF (1000 kg/m 3 ), g is the norm of the gravity and h is the distance between a point P on the mesh and the fluid surface. To compute the impact of the CSF onto the finite element model of the brain, the force f CSF is applied onto every immersed triangle S of the brain surface mesh.

To compute the brain deformation in the presence of CSF and under the boundary conditions described previously, we solve the following static differential system of nonlinear equations:

K(u) u = f CSF + g + H T λ
where K is the stiffness matrix that depends on the displacement u between initial configuration and deformed configuration, g is the acceleration due to gravity and H T λ gathers constraints response resulting from unilateral contacts and Dirichlet boundary conditions.

Risk Visualization

To assess the risk of perforating with vital structures during the insertion of the electrode, we propose to visualize the result of the brain deformation in a natural way. To this end, we introduce the notion of risk volume associated to an electrode path. This risk volume is inspired from current clinical practice. In classic pre-operative planning, a constant safety margin is defined around the planned trajectory (see Fig. 1). In practice, this safety margin defines a cylinder, centered around the electrode axis, with a radius of a few millimeters. This safety margin encompasses various uncertainties in the location of the target and other anatomical structures, due to image registration errors, as well as a potential brain shift. In this paper, we propose to compute more precisely the safety margin as a risk volume that depends on the brain tissue motion, rather than a simple geometry. This risk volume is a combination of an omnidirectional geometric safety margin representing the unpredictable image processing errors, and a physics-based directional safety margin accounting for brain deformation. The computation of the physics-based safety margin is based on our biomechanical simulation, therefore the risk volume handles: the trajectory angle in the patient frame, the depth in the brain tissue and the orientation of the patient's head compared to gravity direction.

As mentioned previously, our approach consists in visualizing the risk that a structure might shift on the trajectory at the time of the surgery. We use this assumption in our algorithm to compute the risk volume associated to a given trajectory: let V be the risk volume. Now, let T be the trajectory, i.e. the segment defined between the target and a point on the surface of the skull, and P ∈ T the points defining the trajectory segment. We also define the displacement field Φ joining the undeformed configuration and the deformed configuration due to the brain shift. We now compute the points P 0 such as ∀P ∈ T, Φ(P 0 ) = P, hence P 0 = Φ -1 (P). V , T , and an example of P and P 0 are depicted in figures 3 and 4. To compute P 0 , we use the interpolation functions N e i of an element e, with i ∈ [START_REF] Allard | Volume contact constraints at arbitrary resolution[END_REF][START_REF] Clatz | Robust nonrigid registration to capture brain shift from intraoperative MRI[END_REF] in case of linear tetrahedrons. Let us find the element e such as P ∈ e, then using the shape function of the tetrahedral finite element model, we have P 0 = ∑ 4 i=1 N e i (P) x e,i 0 .

Here {x e,i 0 |i ∈ [START_REF] Allard | Volume contact constraints at arbitrary resolution[END_REF][START_REF] Clatz | Robust nonrigid registration to capture brain shift from intraoperative MRI[END_REF]} is the coordinates of the nodes in the element e of the underformed brain. Finally, we obtain a volume risk

V = {x ∈ R 3 | x -x p ≤ r, ∀x p ∈ [P, Φ -1 (P)]
, ∀P ∈ T }, with r an error parameter handling the errors due to the model and the image processing uncertainties. The following section presents our results.

Results

In this section, we present the results of our method, applied on a high-fidelity anatomical 3D model of the brain. It includes brain tissues, skull, falx cerebri, ventricles, vessels and the target (subthalamic nucleus). In our tests, we assume the target coordinates are already defined. First results are shown in figures 5 and 6. In these tests, we compare the physics-based risk volume and the geometry-based risk volume. We notice the physicsbased volume is larger, making the trajectory selection more restrictive (see Fig. 8 compared to Fig. 1). Moreover, we observe the restriction acts mainly in the direction of the brain shift. Finally, we notice that the dimension of the volume in the direction of the brain shift varies depending on the depth in the brain tissues (see Fig. 5 and6). This is due to the fact that the tissues on the surface have a larger motion than the deep tissues.

Our method applies in the pre-operative planning step, at which point we do not know the CSF volume that will be lost during the surgery. Besides, there is no precise model to anticipate it. The CSF loss depends on many parameters: orientation of the patient's head, burr hole location, anatomy, pathology etc. For this reason, we ask the surgeon to select an amount of CSF loss, expressed in percentage from 0% (no CSF loss) to 100% (all the CSF of the intracranial cavity has leaked out). Different CSF losses have been tested, and are depicted in figure 7. It appears that the more the brain shift, the larger the risk volume, as expected. Although, the amount of CSF loss is difficult to evaluate, we have determined from the literature [START_REF] Elias | Cortical and subcortical brain shift during stereotactic procedures[END_REF] that a realistic spectrum of percentage is between 0%-40% (knowing that 40% is a huge value that happens in some rare cases). Depending on the pathology and its known relationship with the amount of brain shift (e.g. atrophy of the brain), the surgeon can adjust the value of CSF loss to better estimate the intra-operative scenario. It is also known that the surgical technique has an impact on the CSF loss (e.g. dura mater opening). There is obviously a duality between the brain shift and the CSF loss volume, therefore we can also define a brain shift amplitude based on a displacement, rather than a CSF loss.

Conclusion

In conclusion, we propose a method to intuitively visualize a risk assessment in DBS planning. The computation is based on a physics-based brain shift simulation, handling the trajectory angle, the depth in the brain tissue, and the patient's head orientation. Morevoer, with a more advanced biomechanical model of the brain, we will be able to include more variations in the risk volume: heterogeneity and anisotropy of the brain tissue, interactions with the vascular network.

As there is a duality between brain shift and CSF loss, our next step is to use intraoperative images to retrieve the brain state with the fitting parameters in our brain shift simulation.

Figure 2 .

 2 Figure 2. Schematic representation of the brain model and the set of constraints describing its interactions with the environment. In this pre-operative configuration, the patient is in the supine position (see the gravity vector direction -→ g ). The cerebro-spinal fluid surrounds the brain tissue and acts on it with pressure forces. The resultant force balances the weight (see the --→ f CSF compared to the weight force m -→ g ). The brain is under a null displacement constraint near the brainstem area. The illustration also shows the potential contacts between the brain and the skull.

Figure 3 .

 3 Figure 3. Schematic representation of the brain shift phenomenon. Due to the CSF loss (compared to figure 2), the fluid forces do not counterbalance the gravity force (see the two vectors). This leads to a brain shift, and the structures located in the vicinity of the electrode trajectory may move towards the electrode, increasing the risk of damaging vital structures. The risk volume V (section 1.2) is represented with the associated displacement field Φ.

Figure 4 .

 4 Figure 4. 2d illustration showing an example of a point P on the trajectory, and its associated point P 0 = Φ -1 (P). V is represented depending on r.

Figure 5 .

 5 Figure 5. 3D representation of the physics-based risk volume in the cerebral environment. The risk volume is computed for a trajectory of an electrode implanted in the right hemisphere, represented in transparent.

Figure 6 .

 6 Figure 6. Close-up in the 3D representation of the physics-based risk volume (darker shape), compared to the geometry-based safety margin (cylinder).

Figure 7 .

 7 Figure 7. Variation of the risk volume depending on the CSF loss volume: from 10% lost (left), to 40% (right). The same trajectory for the four simulations is drawn in white.

Figure 8 .

 8 Figure 8. With our visualization of brain shift risk, vessels which can shift toward a trajectory are detected.

Corresponding Author: E-mail: alexandre.bilger@inria.fr.

Acknowledgments

The authors would like to thank the French Research Agency (ANR) for funding this study through the Acoustic project (ANR 2010 BLAN 020901). The authors would also like to thank the Acoustic members for their participation in the useful discussions about this work: Éric Bardinet, Tiziano D'Albis, Caroline Essert, Sara Fernandez Vidal, Noura Hamze, Pierre Jannin, Gregory Laheurte and Jérôme Yelnik.


Alonso Sánchez
email: sanchezsec@lirmm.fr

Chao Liu
email: liu@lirmm.fr

Nabil Zemiti
email: zemiti@lirmm.fr

Etienne Dombre
email: dombre@lirmm.fr

Philippe Poignet
email: poignet@lirmm.fr

Control of meso-robots for endoluminal surgery

Keywords: safety, design framework, adaptive control, dualquaternion, bilateral teleoperation, surgical robotics, endoluminal surgery I

This work deals with the design, the implementation and the experimental validation of a forcereflecting teleoperation architecture for robotic assisted endoluminal surgery. On that account, a meso-scale surgical robot was used as target device to be teleoperated. This mesorobot was developed during the ARAKNES (Array of Robots Augmenting the Kinematics of Endoluminal Surgery) FP7 project, in which the key idea consisted on transferring the robotic degrees of freedom that are required for surgery inside the peritoneal cavity.

INTRODUCTION

In order to build the desired teleoperation architecture, several works were carried out:

1. First, the specifics of designing surgical robots were investigated in order to propose a systematic design framework for surgical robotic components that could facilitate the technology transfer from research centers to the operating room.

2. Second, by applying the proposed framework, a generic real-time software architecture for surgical robots was proposed and implemented. This software architecture serves as support for the global teleoperation architecture.

3. Third, the latter teleoperation architecture was developed. This architecture uses adaptive control techniques in which:

A) An active state observer (AOB) is employed. Moreover, the original Kalman equations (proper of the AOB) were revised and reformulated in order to approximate a Markov chain Monte-Carlo method that allows to better consider, in real-time, nonlinear uncertainties and disturbances that affect the mesorobotic system. B) A realistic viscoelastic Kelvin-Boltzmann softtissue environment model is also used. C) All fundamental (robot) control equations are formulated and performed in the dual-quaternion domain. To that end, one of the first dual quaternion based control toolboxes was implemented.

Such combination of tools and techniques, allowed to demonstrate improved performance of the system with respect to previous works in the field. In this context, two teleoperation schemes of types ``position--position'' and ``position--force'' were proposed in order to fulfill the specific hardware requirements of ARAKNES. Both schemes were analysed in terms of their transparency and stability, i.e. the two parameters that allow quantifying the performance of teleoperated systems.

In addition, having implemented a stable real-time teleoperation control system (and the associated control toolboxes) allowed in turn to study the effects of scaling the positions and the forces, which are respectively commanded and felt by the surgeon. These effects were examined from the point of view of the surgeon's performance and also of the teleopoeration system's performance. In both cases, improvements were observed and it was additionally shown that the availability of online estimated environment's parameters could be exploited to offer surgeons new assistance capabilities that facilitate tasks such as detection of contact and/or discrimination of tissues.

Finally, the possibility of performing wireless control of surgical robots was also explored and experimentally demonstrated. The purpose of this latter study was to increase the mobility of the meso-robot and, possibly, to improve its mechanical performance and/or degree of miniaturization.