Student Member, IEEE Jordan Frecon

Member, IEEE Nelly Pustelnik

Fellow, IEEE Patrice Abry

Member, IEEE Laurent Condat
email: laurent.condat@gipsa-lab.grenoble-inp.fr

On-the-fly Approximation of Multivariate Total Variation Minimization

Keywords: Constrained optimization, Multivariate, Total Variation, On-the-fly algorithm

In the context of change-point detection, addressed by Total Variation minimization strategies, an efficient on-thefly algorithm has been designed leading to exact solutions for univariate data. In this contribution, an extension of such an onthe-fly strategy to multivariate data is investigated. The proposed algorithm relies on the local validation of the Karush-Kuhn-Tucker conditions on the dual problem. Showing that the nonlocal nature of the multivariate setting precludes to obtain an exact on-the-fly solution, we devise an on-the-fly algorithm delivering an approximate solution, whose quality is controlled by a practitioner-tunable parameter, acting as a trade-off between quality and computational cost. Performance assessment shows that high quality solutions are obtained on-the-fly while benefiting of computational costs several orders of magnitude lower than standard iterative procedures. The proposed algorithm thus provides practitioners with an efficient multivariate change-point detection on-the-fly procedure.

I. INTRODUCTION

T otal Variation (TV) has been involved in a variety of signal processing problems, such as nonparametric function estimation [START_REF] Mammen | Locally adaptive regression splines[END_REF], [START_REF] Davies | Local extremes, runs, strings and multiresolution[END_REF] or signal denoising [START_REF] Harchaoui | Catching change-points with Lasso[END_REF]- [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF]. The first contributions on this subject were formulated within the framework of taut string theory [START_REF] Mammen | Locally adaptive regression splines[END_REF], [START_REF] Davies | Local extremes, runs, strings and multiresolution[END_REF] while the term TV had first been introduced in image restoration [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF]. The equivalence between both formalisms has been clarified in [START_REF] Grasmair | The equivalence of the taut string algorithm and BVregularization[END_REF].

Formally, the univariate TV framework aims at finding a piece-wise constant estimate x ∈ R N of a noisy univariate signal y ∈ R N by solving the following non-smooth convex optimization problem,

x = arg min x=(x k ) 1≤k≤N 1 2 x -y 2 + λ N -1 k=1 |x k+1 -x k |, (1)
algorithms. On the one hand, solving this problem in the primal space requires to deal with the non-differentiability of the ℓ 1 -norm that is either handled by adding a small additional smoothing parameter [START_REF] Vogel | Computational Methods for Inverse Problems[END_REF] or by considering proximal algorithms [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF], [START_REF] Steidl | Splines in higher order TV regularization[END_REF]- [START_REF] Alais | Group fused lasso[END_REF]. On the other hand, one can make use of the Fenchel-Rockafellar dual formulation [START_REF] Steidl | Splines in higher order TV regularization[END_REF], [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] or Lagrangian duality [START_REF] Kim | An interiorpoint method for large-scale ℓ 1 -regularized least squares[END_REF], [START_REF] Wahlberg | An ADMM algorithm for a class of total variation regularized estimation problems[END_REF] that can be solved with quadratic programming techniques [START_REF] Steidl | Splines in higher order TV regularization[END_REF], [START_REF] Barbero | Fast Newton-type methods for total variation regularization[END_REF]. Both primal and dual solutions suffer from high computational loads, stemming from their iterative nature. To address the computational load issue, alternative procedures were investigated, such as the taut string algorithm of common use in the statistics literature [START_REF] Mammen | Locally adaptive regression splines[END_REF]. Very recently, elaborating on the dual formulation and thoroughly analysing the related Karush-Kuhn-Tucker (KKT) conditions, a fast algorithm has been proposed by one of the authors in [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF] to solve the univariate optimization problem [START_REF] Mammen | Locally adaptive regression splines[END_REF]. Compared to the taut string strategy, it permits to avoid running sum potentially leading to overflow values and thus numerical errors. Another specificity concerns its on-the-fly behavior that does not require the observation of the whole time sequence before a solution can be obtained. On-the-fly algorithms might be of critical interest for real-time monitoring such as in medical applications [START_REF] Godinez | On-line fetal heart rate monitor by phonocardiography[END_REF], [START_REF] Cannon | An algorithm for online detection of temporal changes in operator cognitive state using real-time psychophysiological data[END_REF].

Along another line, extension of the univariate optimization problem [START_REF] Mammen | Locally adaptive regression splines[END_REF] to multivariate purposes has been recently investigated in [START_REF] Vert | Fast detection of multiple change-points shared by many signals using group LARS[END_REF], [START_REF] Bleakley | The group fused Lasso for multiple change-point detection[END_REF], [START_REF] Hinterberger | Tube methods for BV regularization[END_REF]. The multivariate extension arises very naturally in numerous contexts, such as biomedical applications, for which the purpose is to extract simultaneous change points from multivariate data, e.g., EEG data [START_REF] Gramfort | Mixed-norm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. It also encompasses denoising of complex-valued data, which can naturally be interpreted as bivariate data. Multivariate optimization is known as the group fused Lasso in the statistics literature [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF], [START_REF] Liu | An efficient algorithm for a class of fused Lasso problems[END_REF]. From a Bayesian point of view, elegant solutions have been proposed in [START_REF] Dobigeon | Joint segmentation of multivariate astronomical time series: Bayesian sampling with a hierarchical model[END_REF], [START_REF] Picard | A statistical approach for array CGH data analysis[END_REF] and efficient iterative strategies have recently been proposed in [START_REF] Bleakley | The group fused Lasso for multiple change-point detection[END_REF], [START_REF] Hoefling | A path algorithm for the fused Lasso signal approximator[END_REF].

Mutivariate on-the-fly TV. In this context, the present contribution elaborates on [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF] to propose an on-the-fly algorithm solving the multivariate extension of [START_REF] Mammen | Locally adaptive regression splines[END_REF]. In Section II, the group fused Lasso problem is first detailed. It is then illustrated that the multivariate procedure has a non-local behavior as opposed to the local nature of the univariate problem [START_REF] Mammen | Locally adaptive regression splines[END_REF]. Consequently, any on-the-fly algorithm solving the multivariate minimization problem will only lead to an approximate solution. The KKT conditions resulting from the dual formulation of the multivariate problem are specified in Section III. From such conditions, a fast and on-the-fly, yet approximate, algorithm is derived in Section IV. The performance in terms of achieved solution and computational gain are presented in Section V. A video demonstrating the on-the-fly behavior of the algorithm is available at http://perso.ens-lyon.fr/jordan.frecon1 .

Notations. Let u = (u m,k ) 1≤m≤M,1≤k≤N ∈ R M×N denote a multivariate signal, where for every m ∈ {1, . . . , M }, u m = (u m,k ) 1≤k≤N ∈ R N stands for the m-th component while the k-th values will be shortened as u k = (u m,k ) 1≤m≤M ∈ R M . For every k ∈ {1, . . . , N }, use will also be made of the following functions:

abs(u k ) = (|u m,k |) 1≤m≤M ∈ R M and sgn(u k ) = sgn(u m,k ) 1≤m≤M ∈ R M .

II. LOCAL VS NON-LOCAL NATURE

We denote y the multivariate signal of interest. A multivariate extension of (1) reads:

x = arg min x=(x 1 ,...,x M ) 1 2 M m=1 xm -y m 2 + λ N-1 k=1 M m=1 |(Lxm) k | 2 , (2)
where λ > 0 denotes the regularization parameter and L ∈ R (N -1)×N denotes the first order difference operator, that is, for m ∈ {1, . . . , M } and k ∈ {1, . . . , N -1},

(Lx m ) k = x m,k+1 -x m,k . (3)
Despite formal similarity, there is however a fundamental difference in nature between the univariate (M = 1) and multivariate (M > 1) cases: The former is intrinsically local [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF], [START_REF] Louchet | Total variation as a local filter[END_REF] while the latter is non-local2 . To make explicit such a notion, we have designed the following experiment, whose results are illustrated in Fig. 1. The results associated to the univariate (resp. bivariate) case are presented on the right plots (resp. left plots). A univariate signal y ∈ R N with N = 180, consisting of the additive sum of a piece-wise constant signal and white Gaussian noise (in gray, in Fig. 1, right top plot), is considered first. The solution of the minimization problem (1) is displayed in solid red lines in Fig. 1. Also, we search for the solution of the minimization problem (1) applied to two partitions of y, obtained by splitting it in half, i.e., y -= (y k ) 1≤k≤N/2 and y + = (y k ) N/2+1≤k≤N . The solutions x -and x + of (1) respectively associated to y -and y + are concatenated and displayed with dashed blue lines in Fig. 1. There is strictly no difference between x and the concatenation of x -and x + , as reported in Fig. 1 (bottom right plot), except for the segment that contains the concatenation point. The difference around the concatenation point is expected as x makes use of an information (the continuity between y -and y + ) that is not available to compute x -and x + . The fact that there is no difference elsewhere shows the local nature of the univariate solution to Problem (1).

This experiment is now repeated for M = 2 (as the simplest representative of M > 1). A bivariate signal y = (y 1 , y 2 ) ∈ R 2×N with N = 180, consisting of the additive sum of piecewise constant signals and white Gaussian noises (in gray, in Fig. 1, left plots, 1st and 3rd lines), is considered. Two partitions, obtained by splitting in half, y -= (y 1,k , y 2,k ) 1≤k≤N/2 and y + = (y 1,k , y 2,k ) N/2+1≤k≤N are also considered. The corresponding solutions of (2), applied to y, y -, y + , labeled x, x -and x + are obtained by means of the primal-dual algorithm proposed in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] with λ = 20. Solutions x -and x + of (2) respectively associated to y -and y + are concatenated and displayed with dashed blue lines in Fig. 1, while x is shown in red. Contrary to the case M = 1, differences between x and concatenated x -and x + , shown in black in bottom plots, differ unambiguously from 0 over the entire support of y, clearly showing the non-local nature of x when M > 1.

In the univariate case (Eq. ( 1)), the local nature of the solution permits to design an efficient taut string algorithm, that consists in finding the string of minimal length (i.e., taut string) that holds in the tube of radius λ around the antiderivative of y. The solution x of (1) is then obtained by computing the derivative of the taut string. An efficient strategy has been proposed in [START_REF] Davies | Local extremes, runs, strings and multiresolution[END_REF] in order to straightforwardly compute x by determining the points of contact between the taut string and the tube. Even though this approach can be generalized to multivariate signals, the detection of points of contact additionally requires the angle of contact between the taut string and the tube. However, this information is nonlocal and thus the on-the-fly minimization problem results in a challenging contact problem which can not be solved locally. This interpretation will be further discussed in Section III.

The non-local nature of the multivariate (M > 1) Problem [START_REF] Davies | Local extremes, runs, strings and multiresolution[END_REF] implies that one cannot expect to find an exact multivariate on-the-fly algorithm. Therefore, in the present work, we will derive an approximate on-the-fly algorithm that provides us a good-quality approximation of the exact solution to Problem [START_REF] Davies | Local extremes, runs, strings and multiresolution[END_REF]. A control parameter |Q|, defined in Section 3, will control the trade-off between the quality of the approximation and the computational cost.

III. MULTIVARIATE TOTAL VARIATION MINIMIZATION

A. Dual formulation

Fenchel-Rockafellar dual formulation3 of (2) reads :

minimize u∈R M ×(N -1) 1 2 M m=1 y m + L * u m 2 subject to (∀k = {1, . . . , N -1}) u k ≤ λ, (4)
where, for every m ∈ {1, . . . , M } and k = {2, . . . , N -2},

(L * u m ) k = u m,k-1 -u m,k (5)
and

(L * u m ) 1 = -u m,1 , (L * u m ) N = u m,N -1 . (6)
The optimal solutions u ∈ R M×(N -1) and x ∈ R M×N of the dual problem and of the primal problem respectively are related by

(∀m ∈ {1, . . . , M }) x m = y m + L * u m , (∀k ∈ {1, . . . , N -1}) u k ∈ -λ∂ • ( x k+1 -x k ). (7)
From [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF], we directly obtain the following necessary and sufficient conditions.

Proposition III.1 The solutions of the primal problem (2) and the dual problem (4) satisfy, for every m ∈ {1, . . . , M },

x m = y m + L * u m , (8)
and, for every k = {1, . . . , N -1},

if x k = x k+1 then u k ≤ λ, if x k = x k+1 then u k = -λ x k+1 -x k x k+1 -x k . (9)
The first condition corresponds to the configuration where every component keeps the same value from location k to k + 1. This configuration is illustrated in the bivariate case (M = 2) in Fig. 2 (left plot). The second condition models situations where some components of x admit change points between locations k and k + 1. An interesting configuration is that of non-simultaneous change points as illustrated in Fig. 2 (right plot). In the presence of noise, this second situation is rarely encountered. Thus, in the sequel, we will only consider simultaneous change points.

Remark III.2 Proposition III.1 for M = 1 leads to the usual KKT conditions associated to the minimization problem (1):

     if x k > x k+1 then u k = +λ, if x k < x k+1 then u k = -λ, if x k = x k+1 then u k ∈ [-λ, +λ]. (10)
The on-the-fly univariate TV algorithm proposed in [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF] is derived from Conditions [START_REF] Steidl | Splines in higher order TV regularization[END_REF].

u 1,k > 0 u k λ u 2,k > 0 u k = (0, λ) λ
Fig. 2: Comparing joint vs disjoint changes in the dual space.

Left: location k is suitable for a joint negative amplitude change on both components. Right: configuration suitable for introducing a negative amplitude change at k on the second component only.

B. Rewriting the KKT conditions

Contrary to Conditions [START_REF] Steidl | Splines in higher order TV regularization[END_REF], the multivariate conditions derived in Proposition III.1 are not directly usable in practice to devise an on-the-fly algorithm because x k+1x k is a priori unknown at instant k. Therefore, we propose to rewrite the second condition in ( 9) by means of auxiliary variables

( z k ) 1≤k≤N -1 such that if x k = x k+1 then u k ≤ λ, if x k = x k+1 then u k = -sign( x k+1 -x k ) • z k , ( 11
) with z k = λ abs( x k+1 -x k ) x k+1 -x k
and where • denotes the Hadamard product. Then, Proposition III.1, can be reformulated component-wise as follows.

Proposition III. [START_REF] Harchaoui | Catching change-points with Lasso[END_REF] The solutions of the primal problem (2) and of the dual problem (4) satisfy the following necessary and sufficient conditions. There exist nonnegative auxiliary variables ( z k ) 1≤k≤N -1 such that, for every m = {1, . . . , M } and k = {1, . . . , N -1},

     if x m,k > x m,k+1 then u m,k = + z m,k , if x m,k < x m,k+1 then u m,k = -z m,k , if x m,k = x m,k+1 then u m,k ∈ [-z m,k , + z m,k ], ( 12
) with z k = λ and x m = y m + L * u m .
Comparing Eqs. ( 10) and ( 12) highlights the similarity between the necessary conditions of the univariate and multivariate minimization problems: Conditions involving λ in the univariate case involve the auxiliary vector z in the multivariate one. The fact that z differs for each pair (m, k) can be interpreted in taut string procedures as the fact that the point of contact with the taut string may vary on the tube of radius λ. This significantly increases the difficulty of deriving an on-the-fly algorithm.

C. Approximate solution

If we first assume that z is known and such that, for every k ∈ {1, . . . , N -1}, z k = λ, the primal problem associated to Conditions [START_REF]Proximal splitting methods in signal processing[END_REF] reads

min x M m=1 1 2 y m -x m 2 + N -1 k=1 z m,k |(Lx m ) k | (13)
and can be interpreted as M univariate minimization problems having time-varying regularization parameters ( z m ) 1≤m≤M .

The proposed approximation consists in restricting the estimation of z to a predefined set Q = {ζ (1) , . . . , ζ (|Q|) } chosen such as for every q ∈ {1, . . . , |Q|},

ζ q = (ζ (q) m ) 1≤m≤M ∈ R M satisfies ζ (q) = λ.
The most naive strategy would consist in solving M univariate minimization problems for every |Q| candidate values of z, i.e., find for every m = {1, . . . , M } and q = {1, . . . , |Q|},

x (q) m = arg min xm 1 2 y m -x m 2 + ζ (q) m Lx m 1 (14)
and to devise a method to pick the solution amongst the |Q| candidates. For instance, the one that maximizes some quality criterion f , i.e.,

x = x (q * ) with q * = arg max 1≤q≤|Q| f ( x (q) ). (15)
Although it benefits from parallel on-the-fly implementations, this situation would correspond to a constant estimate z = ζ (q * ) 1 N . Therefore, changes in the mean would be processed independently on all components and group-sparsity would not be enforced.

In order to benefit from an on-the-fly implementation and to enforce group-sparsity, we propose an algorithmic solution based on a piece-wise constant estimator of z detailed in the next section.

IV. ALGORITHMIC SOLUTION

In the following, we first extend the on-the-fly algorithm proposed in [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF] to the multivariate case, with z assumed to be known a priori. This strong assumption, unrealistic in pratice, permits to describe clearly the behaviour of the multivariate on-the-fly algorithm. Then, we will focus on the question of the automated and on-the-fly estimation of z taking its values in Q, which consequently introduce a parameter |Q| controlling the quality of the approximation. The main steps of the on-the-fly algorithm are summarized in Algorithm 1. It is based on the range control of both unknown primal and dual solutions x and u by lower and upper bounds updated with the incoming data stream.

The design of Algorithm 1 results in specifying Rule 1 and Rule 2 allowing respectively to detect a change point and to find suitable change-point locations according to Proposition III.3.

A. Ideal case with z known 1) Lower and upper bounds:

According to Proposition III.3, the solution of the primal problem, the solution of the dual problem and the auxiliary variable have to satisfy, for every k ∈ {0, . . . , N -1},

     u k+1 = y k+1 + u k -x k+1 , abs( u k+1 ) ≤ z k+1 , z k+1 = λ. ( 16
)
Algorithm 1: On-the-fly scheme for multivariate TV with u 0 = u N = 0. Considering the two first conditions, the prolongation condition x k+1 = x k leads to

y k+1 ≥ x k -z k+1 -u k , y k+1 ≤ x k + z k+1 -u k . ( 17
)
Following the solution proposed for the univariate case derived in [START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF], one can check that ( 17) is satisfied by reasoning on lower and upper bounds of u k and x k . For every k ∈ {1, . . . , N -1}, we define the lower and upper bounds of x k , labeled x k and x k respectively, as:

x k ≤ x k ≤ x k , (18)
and we set u k and u k as follows

(∀m ∈ {1, . . . , M }) u m,k = u m,k if x m,k = x m,k , u m,k = u m,k if x m,k = x m,k , (19)
where u k and u k appear to be the upper and lower bounds of x k respectively, i.e.

u k ≤ u k ≤ u k , (20)
as detailed in Appendix VII-A.

2) Updating rules & Rule 1:

The prolongation condition x k+1 = x k , which has led to [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], becomes

y k+1 ≥ x k -z k+1 -u k , y k+1 ≤ x k + z k+1 -u k . (21)
If the latter condition, labeled as Rule 1, holds, then according to the primal-dual relation, we perform the update of the lower and upper bounds at location k + 1 as follows:

u k+1 = y k+1 + u k -x k , u k+1 = y k+1 + u k -x k , (22)
and

x k+1 = x k , x k+1 = x k . ( 23
)
Remark IV.1 Equivalently, one can systematically update primal (resp. dual) bounds according to [START_REF] Wahlberg | An ADMM algorithm for a class of total variation regularized estimation problems[END_REF] (resp. ( 23)) and verify that the following rewriting of the prolongation condition [START_REF] Kim | An interiorpoint method for large-scale ℓ 1 -regularized least squares[END_REF] holds:

u k+1 ≥ -z k+1 , u k+1 ≤ + z k+1 . (24)
3) Signal prolongation & Rule 2:

If Rule 1 (i.e. Condition [START_REF] Kim | An interiorpoint method for large-scale ℓ 1 -regularized least squares[END_REF] or equivalently [START_REF] Godinez | On-line fetal heart rate monitor by phonocardiography[END_REF]) holds, then the assumption x k+1 = x k is valid. However, the upper and lower bounds may have to be updated in order to be consistent with

u k+1 ∈ [-z k+1 , + z k+1 ].
According to [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], this condition requires to verify that the following Rule 2 holds:

u k+1 ≤ + z k+1 , u k+1 ≥ -z k+1 . ( 25
)
For every m ∈ {1, . . . , M }, three configurations can be encountered:

• When both Conditions (25) are satisfied, the bounds are left unchanged.

• When u m,k+1 = u m,k + y m,k+1 -x m,k > + z m,
ν m = x m,k + u m,k+1 -z m,k+1 k -k 0 + 1 , (26)
we propose the following updates

(∀j ∈ {k 0 , . . . , k + 1}) x m,j = ν m , u m,k+1 = + z m,k+1 . (27)
• When u m,k+1 <z m,k+1 , then it results that the upper bound ν m ≡ x m,j (∀j ∈ {k 0 , . . . , k + 1}) has been overevaluated. Similarly, since u m,k+1 is lower bounded by z m,k+1 , we can show that the upper bound

ν m = x m,k + u m,k+1 + z m,k+1 k -k 0 + 1 , (28)
permits to ensure the consistency of the following updates

(∀j ∈ {k 0 , . . . , k + 1}) x m,j = ν m , u m,k+1 = -z m,k+1 . (29)
4) Estimate of the change point k rupt : When Rule 1 does not hold, a change point has to be created. For every m ∈ {1, . . . , M }, we can distinguish three cases:

• When u m,k+1 = u m,k + y m,k+1 -x m,k < -z m,k+1 ,
then, since u m,k is bounded, it means that x m,k is overevaluated and therefore a negative amplitude change has to be introduced on the m-th component in the time index set {k 0 , . . . , k} in order to decrease its value. Following Proposition III.3 and Eq. ( 20), the set of locations κ m suitable for a change-point on the m-th component reads:

κ m = {j ∈ {k 0 , . . . , k} | u m,j = + z m,j }. ( 30
)
Such locations correspond to the indexes where the value of the bound u m,j has been updated in order to be consistent with the condition u m,j ∈ [-z m,j , z m,j ] (see the previous paragraph) • When u m,k+1 > + z m,k+1 , then a positive amplitude change has to be introduced in the m-th component within the time index {k 0 , . . . , k}. The set of locations suitable for a change-point on the m-th component reads:

κ m = {j ∈ {k 0 , . . . , k} | u m,j = -z m,j }. ( 31
)
This set of locations corresponds to indexes where the value of the bound u m,j was updated in order to be consistent with u m,j ∈ [-z m,j , z m,j ]. • Else, when component m does satisfy [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF], then we set κ m = {k 0 , . . . , k}. The change-point location k rupt corresponds to the last location suitable for introducing the adequate amplitude change on each component, i.e.,

k rupt = max j∈∩ M m=1 κm j. ( 32
)
Once the change point location has been specified, we are able to assign a value to ( x j ) k0≤j≤krupt . When a negative amplitude change is detected on the m-th component, we set

(∀j ∈ {k 0 , . . . , k rupt }) x m,j = x m,k+1 , (33)
in consistence with [START_REF] Alais | Group fused lasso[END_REF]. Similarly, when a positive amplitude change is detected, we set

(∀j ∈ {k 0 , . . . , k rupt }) x m,j = x m,k+1 . (34)

5) Starting a new segment:

When a segment has been created, we start the detection of a new segment considering k 0 = k rupt + 1 as long as k 0 < N .

According to [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF] and by definition of the bounds, for every k ∈ {1, . . . , N }

x k = y k -u k + u k-1 , x k = y k -u k + u k-1 . (35)
In particular, for k = k 0 , combining ( 12), ( 18), ( 19) and ( 22) allows us to find the following initialization procedure

u k0 = + z k0 , u k0 = -z k0 , x k0 = y k0 -z k0 + u k0-1 , x k0 = y k0 + z k0 + u k0-1 , (36)
where the value of u k0-1 is given according to Proposition III.3. In addition, according to the writing of ( 16), u 0 = 0.

B. Estimation of the auxiliary multivariate vector z

In order to describe the generic behavior of the multivariate on-the-fly algorithm, we have so far assumed z to be known a priori. We now focus on the simultaneous estimation of the multivariate vector z and of the multivariate signal x.

To provide an on-the-fly approximate solution, we propose:

• to build a piece-wise constant estimator z of z,

• to only consider amplitude changes jointly on all components m ∈ {1, . . . , M }.

1) Piece-wise constant estimator of z:

The proposed estimate is assumed to be constant between each changepoint with values belonging to the predefined set Q defined in Section III-C. For each candidate value ζ (q) with q ∈ {1, . . . , |Q|}, we create upper and lower bounds labeled u

(q) k , u (q) k , x (q)
k , and x (q) k . They are initialized at each new segment location k 0 and are updated independently according to [START_REF] Wahlberg | An ADMM algorithm for a class of total variation regularized estimation problems[END_REF] and ( 23) until the prolongation condition

u (q) k+1 ≥ -ζ (q) , u (q) k+1 ≤ +ζ (q) , (37)
based on [START_REF] Godinez | On-line fetal heart rate monitor by phonocardiography[END_REF], does not hold anymore. In the following, we investigate how to modify the algorithm described in Section IV-A, to account for the automated selection of z in Q. The resulting algorithm is reported in Algorithm 2.

2) Estimate of the change point k (q)

rupt :

For every q ∈ {1, . . . , |Q|}, we create change points as described in Section IV-A4. The main difference consists in the restriction to simultaneous change points. As detailed after Proposition III.1, non-simultaneous changes have a zero probability to occur. The restriction to simultaneous change-points will thus not impact the solution. It results that if there exists at least one component m ∈ {1, . . . , M } such that u

(q) m,k+1 < -ζ (q) m (resp. u (q) m,k+1 > ζ (q) m ), then κ (q) m = {j ∈ {k 0 , . . . , k} | u (q)
m,j = +ζ (q) m } (38)

(resp. κ (q) m = {j ∈ {k 0 , . . . , k} | u (q) m,j = -ζ (q) m }), (39)
and, ∀m -= m such that u

(q) m-,k+1 + u (q)
m-,k+1 < 0, then

κ (q) m-= {j ∈ {k 0 , . . . , k} | u (q) m-,j = +ζ (q) m } (40)
or, ∀m + = m such that u

(q) m+,k+1 + u (q)
m+,k+1 ≥ 0, then

κ (q) m+ = {j ∈ {k 0 , . . . , k} | u (q) m+,j = +ζ (q) m }. (41)
A bivariate example of these configurations where the second component breaks Condition (37) is provided in Fig. 3. The change-point location k

(q)
rupt and the assignment of x (q) on the current segment follow (32), [START_REF] Hoefling | A path algorithm for the fused Lasso signal approximator[END_REF] and [START_REF] Louchet | Total variation as a local filter[END_REF].

3) Estimate of the change point k rupt : According to the previous paragraph, the piece-wise estimation procedure leads to several possible change-point locations (at most |Q|). Here we select the solution indexed by q * with tightest bounds x (q * ) and x (q * ) , i.e.,

q * ∈ Argmin 1≤q≤|Q| x (q) k (q) rupt -x (q) k (q) rupt σ -1 2 , ( 42
) with σ = diag(σ 1 , . . . , σ M ), (43)
where, for every m ∈ {1, . . . , M }, σ m stands for the standard deviation of y m . The factor σ -1 permits to ensure that every component contributes equally to the criterion (42). When the

u (q) 2,k+1 ▽ △ u (q) 2,k+1 ⊲ u (q) 1,k+1 u (q) 1,k+1 ⊳ u (q) 2,k+1 ▽ △ ⊳ ⊲ u (q) 2,k+1 u (q) 1,k+1 u (q) 1,k+1
q) 1 = ζ (q) 2 = λ/ √ 2. Since u (q) 2,k+1 > ζ (q) (
2 , condition (37) is violated. The left (resp. right) plot displays the configuration u

(q) 1,k+1 + u (q) 1,k+1 < 0 (resp. u (q) 1,k+1 + u (q) 1,k+1 ≥ 0) described in Section IV-B2.
minimizer of ( 42) is not unique, we select the index q * yielding the largest k (q * ) rupt . In other words, we choose the set of auxiliary variables which permits to hold the prolongation condition (37) as long as possible.

Therefore, it finally leads to an index q * which permits to estimate k rupt = k (q * ) rupt and,

(∀j ∈ {k 0 , . . . , k rupt }) zj = ζ (q * ) , x j = x (q * ) j . ( 44
)
The starting location for the next segment is then, k 0 = k rupt + 1, and the algorithm iterates as long as k 0 < N .

4) Starting a new segment:

Let us consider the location k 0 of a new segment. For every q ∈ {1, . . . , |Q|}, the initialization step can be recast into u (q) k0 = +ζ (q) , u

(q) k0 = -ζ (q) , x (q) k0 = y k0 -ζ (q) + u k0-1 , x (q) k0 = y k0 + ζ (q) + u k0-1 ,
(45) with u 0 = 0.

Remark IV.2 The initialization step (45) implicitly depends on the estimation of z made on the last segment through the term u k0-1 . Simulations have shown that (45) may lead to an inconsistent solution x as soon as z has been poorly estimated on a segment. Empirically, a better approximation of the iterative solution is observed if each segment is treated independently, i.e., u (q) k0 = +ζ (q) , u

(q) k0 = -ζ (q) , x (q) k0 = y k0 -ζ (q) , x (q) k0 = y k0 + ζ (q) .
(46)

V. PERFORMANCE ASSESSMENT

A. Experimental setting

Unless specified otherwise, we consider that data consist of a M -multivariate piece-wise constant signal x ∈ R N (solid black), to which a centered Gaussian noise ǫ is additively superimposed:

y = x + ǫ ∈ R M×N .
Signal x is generated as follows. First the length of each segment is drawn according to a folded Gaussian distribution N (12.5, 16.25). Then, for each m ∈ {1, . . . , M }, the amplitudes of the corresponding changes are drawn independently from a Gaussian distribution N (2, 0.4). Regularization parameter λ > 0.

Predefined set Q = {ζ (1) , . . . , ζ (|Q|) }.

Starting location k0 = 1. while k0 < N do for q ← 1 to |Q| do Set k ← k0 Initialize primal/dual bounds according to (46) while (37) is satisfied do

Set k ← k + 1 for m ← 1 to M do Update primal/dual bounds if u (q) m,k+1 > +ζ (q) m or u (q) m,k+1 < -ζ (q)
m then Revise the update of primal/dual bounds Estimate k (q) rupt and ( x

(q) j ) k 0 ≤j≤k (q) rupt according to Section (IV-B2) Estimate krupt ∈ (k (q) rupt ) 1≤q≤|Q| according to Section IV-B3 Estimate ( xj) k 0 ≤j≤krupt and (zj) k 0 ≤j≤krupt according to (44) Set k0 ← krupt + 1 Result: Solution x
The exact minimizer of (2), labeled x, is computed by means of the ADMM algorithm proposed in [START_REF] Wahlberg | An ADMM algorithm for a class of total variation regularized estimation problems[END_REF]. Iterations are stopped when the relative criterion error is lower than 10 -10 . The proposed solution computed with the predefined set Q is denoted x approx,Q .

In a second set of experiments (see V-D), the proposed onthe-fly algorithmic solution will be compared to an on-the-fly ADMM solution.

B. Design of Q

We propose to compare solutions x approx,Q obtained with two different sets Q = {ζ (1) , . . . , ζ (|Q|) } in the bivariate case (i.e., M = 2) for N = 10 4 . For both configurations, we choose (∀q ∈ {1, . . . , |Q|}) ζ (q) = (λ cos(θ q ), λ sin(θ q )) (47) with θ q ∈ [0, π/2]. The first solution consists to homogeneously cover the ℓ 2 ball such that, for some some parameter R ∈ N * , θ q = qπ/2 R+1 and |Q| = R-1 q ′ =0 2 q ′ . The second solution draws a set of the same size whose values (θ q ) 1≤q≤|Q| follow a uniform distribution on [0, π/2].

Two experimental settings are investigated. In the first one, y 1 is one order of magnitude larger than y 2 (Fig. 4, left plots) whereas in the second one, both are of the same order of magnitude (Fig. 4, right plots).

Estimation performances in terms of mean squared error

MSE( x approx,Q , x) = E[ 1 N x approx,Q -x 2 ]
(where E stands for the sample mean estimator computed over 100 realizations) are reported on the first line. It shows that a random covering of the ℓ 2 -ball provides solutions as good as the homogeneous covering up to the limit of |Q| small.

On the 2nd and 3rd lines, the distributions of θ q * , where q * has been selected by criterion (42), are reported for |Q| = 127. These histograms show the impact of the relative amplitude of the components on the distribution θ q * : components with same order of magnitude yield a symmetric distribution while unbalanced components yield an asymmetric distribution. For instance, in Fig. 4 (right plots), it appears more meaningful to draw θ q according to a Gaussian distribution than to a uniform distribution. Therefore, if one has some knowledge of components amplitudes, this can be incorporated to better design the set Q. This will also decrease the computational cost discussed in section V-D.

In the following, we restrict ourselves to a random covering of the ℓ 2 ball.

C. Offline performance

In this section, we focus on the comparison of offline performance, extended for M = 10, for two different signal-to-noise ratios (SNRs), namely 4dB and 10dB. impact of |Q| on x approx,Q . For a single realization of noise, x approx,Q and x are plotted Fig. 5 for λ = 29, adjusted to provide the best visual (qualitative) performance. Solution x approx,Q for |Q| = 5 × 10 4 (light orange) provides a visually better approximation of x (dashed blue) than for |Q| = 10 3 (mixed red).

Estimation performance x approx,Q vs. x. The quality of the approximation is further quantified Fig. 6 in terms of MSE( x approx,Q , x) as a function of λ for different |Q|.

It shows that the MSE systematically decreases when |Q| increases. Further, on the examples considered here and depending on λ, using |Q| ≥ 10 4 no longer yields significantly improved solutions, thus showing that the selection of |Q| does not require a complicated tuning procedure.

Estimation performance x vs. x and x approx,Q vs. x. Let us now compare the absolute quality of the solutions against x. MSE( x, x) and MSE( x approx,Q , x) for different |Q|, are reported in Fig. 7. MSEs are consistent with the previous paragraph: it shows that increasing |Q| up to a certain value permits to significantly lower the MSE. However, x has a lower estimation error than x approx,Q .

D. Online performance

In this section we focus on the comparison between two online solutions. The first one is derived from the proposed on-the-fly algorithm whereas the second one is based on an iterative algorithm.

Comparison is made for different values of λ on a signal x ∈ R M×N (N = 400) to which a Gaussian noise is superimposed such that SNR = 3dB. Performance are provided for M = 2 and M = 5 components.

Proposed online solution x online,Q . As the time step k increases, x approx,Q is only computed up to the last k 0 and the algorithm has not yet output a solution on {k 0 + 1, . . . , k}. In that sense, the solution is said to be "on-the-fly". However, a solution x online,Q , providing an online approximation of x, can be output up to k by imposing limit conditions at k. Windowed iterative solution x win,K . We consider a naive online ADMM version, where at each time step k a solution x win,K is computed by optimizing over the previous K points. The choice of K is of critical importance. On the one hand, if this value is too small, the observer may miss amplitude changes in the multivariate data stream. On the other hand, if the window size is too large, the computational cost may be too high to handle any online observation. Three window sizes have been investigated, respectively K = 20, 50 and 80.

Computational cost. Comparisons of median computational costs per incoming sample (in seconds), over 10 realizations of noise, are reported Fig. 8 (left plots) as functions of λ.

As expected, we observe that the computational cost does increase along with the size of Q. Therefore, |Q| acts as a trade-off between the computational cost and the MSE. However, the computational cost of x online,Q is still several orders of magnitude lower than the one associated to the online ADMM. Interestingly, computational costs are comparable for M = 2 (top left plot) and M = 5 (bottom left plot). Note that a warm-up starting strategy for online ADMM only reduces by a factor two the computational cost with respect to the implementation displayed in Fig. 8.

The computational cost of x online,Q could still be reduced in two ways. First, one could design the set Q according to a priori knowledge of components amplitudes (see V-B). Second, one could also benefit from the separable form of the algorithm and compute solutions x (q) in parallel for every q ∈ {1, . . . , |Q|}.

Change-point detection accuracy. The Jaccard index J(α, β) ∈ [0, 1] between any α and β ∈ [0, 1] N is defined as [START_REF] Jaccard | Distribution de la flore alpine dans le bassin des Dranses et dans quelques régions voisines[END_REF], [START_REF] Hamon | Discovering the structure of complex networks by minimizing cyclic bandwidth sum[END_REF]

J(α, β) = N i=1 min(α i , β i ) N 1≤i≤N αi>0,βi>0 αi+βi 2 + 1≤i≤N βi=0 α i + 1≤i≤N αi=0 β i .
(48) It varies from 0, when α ∩ β = ∅, up to 1 when α = β. The Jaccard index is a demanding measure: As an example, if β ∈ {0, 1} N is the truth and if α ∈ {0, 1} N has correctly identified half non-zero values of β but has misidentified the other half, then J(α, β) = 1/3.

The Jaccard index is used to measure the similarity between change-point locations of x and those obtained during the computation of x win,K and x online,Q . To this end, we consider the change-point indicator vector r = (r i ) 1≤i≤N of x (as well as r win,K and r online,Q respectively associated to x win,K and x online,Q ), defined as

r i = 1, if x has a change-point at location i, 0, otherwise. (49)
In order to incorporate a tolerance level on change-point locations, r, r win,K and r online,Q are first convolved with a Gaussian kernel of size 10 with a standard deviation of 3. J( r win,K , r) and J( r online,Q , r) are averaged over 10 realizations of noise and reported in Fig. 8 (right plots) as functions of λ for different set size |Q| and window size K.

Performance show that J( r online,Q , r) ≥ J( r win,K , r) for almost all λ and |Q|. Therefore, x online,Q provides a better online detection of change-points of x. It also show that J( r online,Q , r) does not vary significantly with |Q| but slightly decreases with M . Indeed, as M increases, the prolongation condition (37) is more likely to be violated, thus leading to more change-points.

VI. CONCLUSION

In this contribution, we have developed an algorithm which provides an on-the-fly approximate solution to the multivariate total variation minimization problem. Besides a thorough examination of the KKT conditions, the key-step of the algorithm lies in updating and controlling the range of the upper and lower bounds of the dual solution within a tube of radius λ. An on-the-fly derivation is achieved by means of an auxiliary vector z, which needs to be estimated, providing information on the angle of contact with the tube. The latter estimation strongly affects the quality of the solution and the proposed on-the-fly estimation of z is currently achieved by assigning a value chosen within a predefined set Q. It has been shown that the size of Q permits to achieve a desired trade-off between the targeted quality of the solution and the application-dependent affordable computational cost. In addition, the proposed method could also be extended to other ℓ 1,p penalization norms in the right-hand side of (2), for p > 1.

However one would still face the issue of estimating z which would have to lie within a ℓ p ball of radius λ. Under current interest is the investigation of how to estimate z in the case where the assumption of piece-wise constant behavior is a priori relaxed.

VII. APPENDIX

A. Proof of Equation [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF] According to the primal-dual relation [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF], for every m ∈ {1, . . . , M } and k ∈ {1, . . . , N -1},

u m,k = y m,k + u m,k-1 -x m,k , (50)
and, according to [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF], u m,ku m,k ≤ 0. The arguments are similar for proving that u m,k > u m,k .

then updating rules of x m,k , specified in [START_REF] Barbero | Fast Newton-type methods for total variation regularization[END_REF], have underevaluated its value ν m . To modify the lower bounds (x m,j ) k0≤j≤k+1 , on the one hand, we consider the cumulative sum of the observations which, according to [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF], leads to

Fig. 1 :

 1 Fig. 1: Non-local vs. local nature. Left: bivariate TV (upper plots: first component, lower plots: second component). Right: univariate TV. Observations y (gray), solution x (red), concatenation of solutions x-and x+ (dashed blue).

  k+1 , then the updating rules specified in[START_REF] Barbero | Fast Newton-type methods for total variation regularization[END_REF] have under-evaluated the bound ν m ≡ x m,j (∀j ∈ {k 0 , . . . , k + 1}) where k 0 denotes the last starting location of a new segment. Since u m,k+1 is upper-bounded by + z m,k+1 and, that for such a value it can be shown (see Appendix VII-B) that

Fig. 3 :

 3 Fig. 3: Example of configurations leading to the detection of a change-point. In this example M = 2, ζ

Algorithm 2 :

 2 On-the-fly Multivariate TV Data: Multivariate signal y = (y 1 , . . . , y M ) ∈ R M ×N .

Fig. 4 :

 4 Fig. 4: Influence of the design of Q.Comparison over 100 realizations of an homogeneous covering of the ℓ2-ball (blue) against a random covering (red). Two experimental settings are considered depending on whether if y 1 is one order of magnitude larger than y 2 (left) or not (right). Top: MSE( xapprox,Q, x) for different set sizes |Q|. Bottom (2nd and 3rd lines): distributions of θq * where q * has been selected by criterion (42) for |Q| = 127.

Fig. 5 :

 5 Fig. 5: Qualitative impact of |Q| on xapprox,Q. For visibility, only 3 components out of M = 10 are displayed. SNR = 4dB. xapprox,Q for |Q| = 5 × 10 4 is more satisfying than for |Q| = 10 3 since it has more discontinuities in common with x.

Fig. 6 :Fig. 7 :

 67 Fig. 6: Estimation performance xapprox,Q vs x. MSE( xapprox,Q, x) for different |Q|. SNR is set to 4dB (resp. 10dB) on left plot (resp. right plot).

  and by definition of the lower and upper bounds of x m,k and u m,k , we haveu m,k = y m,k + u m,k-1x m,k , (51) u m,k = y m,k + u m,k-1x m,k .(52)By subtracting (51) from (50) we obtainu m,ku m,k = x m,kx m,k

Fig. 8 :

 8 Fig. 8: Online Performance. The proposed solution x online,Q is displayed in solid line while the online ADMM solution xwin,K is displayed in dashed line. Performance for M = 2 (resp. M = 5) are illustrated top (resp. bottom). Left: median computational cost per incoming sample (in seconds). Right: J( rwin,K, r) and J( r online,Q , r) for different values of |Q| and K.

  j = u m,k+1u m,k0 + (kk 0 + 1)x m,k+1 , (55) and thus, if u m,k+1 = + z m,k+1 , would lead tok+1 j=k0+1 y m,j = z m,k+1u m,k0 + (kk 0 + 1)ν m , (56)by definition of x m,k+1 = ν m . On the other hand, the updating rules (22) and (23) have led tou m,k+1 = u m,k0 + k+1 j=k0+1 y m,j -(kk 0 + 1)x m,k . (57)The combinaison of (56) and (57) leads toν m = x m,k + -u m,k0 + u m,k+1z m,k+1 + u m,k0 kk 0 + 1 .(58)Because x m,k have been under-evaluated and by definition u m,k0 ≤ u m,k0 , we can propose the following valueν m = x m,k + u m,k+1z m,k+1 kk 0 + 1 ,(59)in order to adjust the lower bounds, i.e.,(∀j ∈ {k 0 , . . . , k + 1}) x m,j = ν m .(60)In addition, as a result of u m,k+1 ∈ [-z m,k+1 , + z m,k+1 ] and according to the inequality (20), we set u m,k+1 = + z m,k+1 .(61)

  Multivariate signal y = (y 1 , . . . , y M ) ∈ R M ×N .

	Regularization parameter λ > 0.
	Starting location k0 = 1.
	while k0 < N do
	Set k ← k0 Initialize primal/dual bounds
	while Rule 1 is satisfied do
	Set k ← k + 1 for m ← 1 to M do Update primal/dual bounds
	if Rule 2 is not satisfied then
	Revise the update of primal/dual bounds
	Estimate the change point krupt
	Estimate ( xj) k 0 ≤j≤krupt Set k0 ← krupt + 1 Result: Solution xapprox

Data:

A toolbox will be provided at the time of publication.

In this article, we denote a problem as local if the solution at a given location does not depend on the signal located earlier (later) than the previous (next) change-point.

Note that, the usual dual formulation and the resulting stationarity conditions would involve u rather than -u. The choice made in this article enables us to be consistent with the results obtained in[START_REF] Condat | A direct algorithm for 1D total variation denoising[END_REF] for the univariate case.

Jordan Frecon received the Master degree in Physics from the Ecole Normale Supérieure de Lyon (ENS de Lyon), France, in 2013. He is currently working toward a Ph.D. degree in physics and signal processing from the ENS de Lyon. His area of interest includes inverse problems in signal processing, data-driven estimation techniques based on variational approaches and optimization tools for scaling analysis.

Nelly


Jordan Frecon

Nelly Pustelnik

Herwig Wendt
email: herwig.wendt@irit.fr

Laurent Condat
email: laurent.condat@gipsa-lab.grenoble-inp.fr

Patrice Abry

Multifractal-based texture segmentation using variational procedure

Keywords: Local regularity, multifractal spectrum, segmentation, convex optimization, wavelet leaders

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Recent contributions in image processing highlighted the need of segmentation techniques for scale-free textures analysis [START_REF] Nafornita | Semilocal Hurst estimation via generalised lasso and dualtree complex wavelets[END_REF][START_REF] Frecon | Multivariate optimization for multifractal-based texture segmentation[END_REF][START_REF] Pustelnik | Local regularity, wavelet leaders and total variation based procedures for texture segmentation[END_REF]. The possible applications go from texture medical images such as the bone study [START_REF] Benhamou | Fractal analysis of radiographic trabecular bone texture and bone mineral density: two complementary parameters related to osteoporotic fractures[END_REF] to art investigations [START_REF] Abry | When Van Gogh meets Mandelbrot: Multifractal classification of paint-ing's texture[END_REF][START_REF] Johnson | Pursuing automated classification of historic photographic papers from raking light photomicrographs[END_REF].

Scale-free behavior is captured with local regularity, technically measured via the concept of Hölder exponent [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF]. First studies based on the estimation of this quantity can be traced back to [START_REF] Roux | A waveletbased method for multifractal image analysis. III. Applications to high-resolution satellite images of cloud structure[END_REF] where local regularity was assumed homogeneous throughout the image (i.e., monofractal). Recent contributions have considered a more realistic model, where the local regularity may be heterogeneous throughout the image (i.e, piecewise monofractal). This further increases the complexity of the estimation procedure since it additionnaly amounts in segmenting the image into a priori unknown regions where the local regularity can be considered homogeneous.. Among numerous techniques for image segmentation, some efficient variational approaches have recently been designed relying on the use of the total variation [START_REF] Chambolle | A convex approach to minimal partitions[END_REF][START_REF] Nafornita | Semilocal Hurst estimation via generalised lasso and dualtree complex wavelets[END_REF][START_REF] Pustelnik | Local regularity, wavelet leaders and total variation based procedures for texture segmentation[END_REF][START_REF] Nelson | Semilocal scaling exponent estimation with box-penalty constraints and total-variation regularization[END_REF]. The good results of Work supported by GdR 720 ISIS under the junior research project GALILEO, ANR AMATIS grant #112432, 2010-2014, and CNRS Imag'in project under grant 2015OPTIMISME. these estimation and segmentation techniques (on simulated and real data) thus pave the way for considering more complex piecewise scale-free models. This is the subject of the present contribution.

We aim to go further by proposing a segmentation method for piecewise multifractal processes analysis, thus allowing a richer modeling of real-world textures. The multifractal formalism detailled Section 2 relies on the local estimation of multiscale quantities in place of the local regularity. It is combined with a joint vectorial segmentation procedure whose algorithmic solution is briefly summarized in Section 3. Estimation performance on synthetic results is reported in Section 4. Comparisons with state-of-the-art methods are also provided.

MULTIFRACTAL ANALYSIS

Piecewise mono/multi-fractal. We denote X = (X ) 1≤ ≤N the scale-free image to analyze, which has N pixels. Its local regularity can be quantified by the Hölder exponent [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF] denoted by h = (h ) 1≤ ≤N . While large h points to a locally smooth portion of the field, low h indicates local high irregularity. Two examples of X are provided in Figure 1(b)-(c).

On the one hand, in Figure 1(b), X models a piecewise monofractal process having two different values of h , i.e.

(∀ ∈ Ω

A ) h = h A , (1)
(∀ ∈ Ω B ) h = h B , (2)
with h A < h B and Ω = Ω A ∪ Ω B being separable in two distinct areas, i.e., Ω A ∩ Ω B = ∅, according to the mask presented in Figure 1 (a). A smoother behavior is thus observed for ∈ Ω A than for ∈ Ω B . On the other hand, in Figure 1 (c), X models a piecewise multifractal process where h may locally vary both within Ω A and Ω B . Therefore, the multifractal spectrum D(h), which describes local regularity fluctuations, is not reduced to a Dirac (cf. e.g., [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF] for details), and The aim of multifractal analysis is to estimate D(h). For practical purposes, the multifractal spectrum can often be approximated as a parabola:

(∀ ∈ Ω A ) h can be described by D A (h) (3) (∀ ∈ Ω B ) h can be described by D B (h). (4)
D(h) = 2 + (h -c 1 ) 2 /(2c 2 ).
Here, we follow the efficient procedure proposed in [START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF] based on wavelet-leader coefficients [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF].

Wavelet-leader coefficients. We denote d

(m) j,k = X, ψ (m) j,k
the (L 1 -normalized) 2D discrete wavelet coefficients of X at location k = 2 -j , at scale 2 j with j ∈ {1, . . . , J}, and where m stands for the horizontal/vertical/diagonal subband. For a detailed definition of the 2D-DWT, readers are referred to e.g., [START_REF] Mallat | A wavelet tour of signal processing[END_REF]. The wavelet leader coefficient L j,k , located around position = 2 j k, is defined as the local supremum of all wavelet coefficients taken within a spatial neighborhood across all finer scales 2 j ≤ 2 j , that is,

L j,k = sup m=1,2,3, λ j ,k ⊂Λ j,k |d (m) j ,k |, (5)
with λ j,k = [k2 j , (k+1)2 j ), and Λ j,k = p∈{-1,0,1} 2 λ j,k+p . An illustration is provided in Figure 2 where L j,k is indicated with a black cross and the neighborhood Λ j,k is displayed in red.

Hölder exponent. The Hölder exponent can be obtained by a linear regression of wavelet-leader coefficients as follows:

h k = j∈{1,...,J} w j,k log 2 L j,k , (6)
where the w j,k model regression weights [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF].

Multifractal spectrum. The multifractal spectrum can be obtained by multiscale quantities C 1,j ∈ R N and C 2,j ∈ R N defined as the sample estimates of the first and second cumulant of log 2 L j at each given scale 2 j . It has been shown that C 1,j and C 2,j are related to the multifractal spectrum D(h) via the coefficients c 1 and c 2 as follows [START_REF] Wendt | Wavelet leaders and bootstrap for multifractal analysis of images[END_REF]: To evaluate changes in multifractal spectrum, one could naturally consider estimating C 1,j and C 2,j locally in a neighborhood of each pixel , estimate the corresponding local parameters c 1 and c 2 , and then perform a vectorial segmentation of (c 1 , c 2 ). However, this relies on the strong assumption that real-world textures follow precisely the scaling behaviors prescribed in Eqs. ( 7) and ( 8) above. In the present contribution, it has been chosen to relax this requirement and to directly perform a vectorial segmentation over multiscale quantities C 1,j and C 2,j , possibly smoothed.

EC 1,j = c 0 1 + c 1 ln 2 j , (7)
EC 2,j = c 0 2 + c 2 ln 2 j . (8)

VECTORIAL SEGMENTATION

We propose to follow segmentation procedure ideas that perform well for piecewise mono/multi-fractal estimation [START_REF] Frecon | Multivariate optimization for multifractal-based texture segmentation[END_REF].

The vectorial segmentation procedure is based on a convex relaxation of the piecewise constant Potts segmentation problem [START_REF] Chambolle | A convex approach to minimal partitions[END_REF][START_REF] Storath | Fast partitioning of vector-valued images[END_REF]. Here we consider the extension to joint vectorial segmentation proposed in [START_REF] Chambolle | A convex approach to minimal partitions[END_REF] Problem formulation. The labeling procedure of Y into Q level sets can be formalized as a minimization problem where

Q + 1 binary images Θ = (θ q ) 1≤q≤(Q+1) ∈ R N ×Q+1 are estimated such that min Θ Q q=1 (θ q -θ q+1 ) M m=1 (Y m -v q,m ) 2 +λ Q q=1 TV(θ q -θ q+1 ) subject to      θ 1 = 1, θ Q+1 = 0, 1 ≥ θ 2 ≥ . . . ≥ θ Q ≥ 0, (9)
with λ > 0 and where TV denotes the usual total-variation penalization as defined in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], i.e., for every θ ∈ R N ,

TV(θ) = N =1 (Dθ) 2 , (10)
where D ∈ R 2N ×N stands for the discrete horizontal/vertical difference operator and thus (Dθ) ∈ R 2 . The choice of v q,m ∈ R will be discussed further.

The Q resulting labelling areas (Ω 1 , . . . , Ω Q ) are obtained from the binary images (θ q -θ q+1 ) 1≤q≤Q as follows:

(∀ ∈ Ω) θ q, -θ q+1, = 1, if ∈ Ω q , 0, otherwise. (11)
The first term in ( 9) is a data fidelity term allowing to impose similar properties within each area Ω q . The second term imposes the regularity for each labelling area Ω q . The smaller is λ the higher is the granularity of Ω q .

Algorithmic solution. In order to efficiently estimate Θ we use the algorithmic strategy proposed in [START_REF] Condat | Segmentation d'image par optimisation proximale[END_REF] that consists in using a proximal splitting method coupled with an efficient strategy to compute the involved proximal operators. The reader could refer to [START_REF] Chambolle | A convex approach to minimal partitions[END_REF][START_REF] Condat | Segmentation d'image par optimisation proximale[END_REF] for details regarding the algorithmic strategy for solving (9).

EXPERIMENTS

Experimental setting. Performance of the proposed segmentation procedure are assessed on synthetic data, numerically produced by inclusion of a 2D MRW patch Ω A [START_REF] Robert | Gaussian multiplicative chaos revisited[END_REF] into a 2D-MRW background Ω B with different multifractal parameters. Patch and background (Figure 3 (b)) have been normalized to ensure that the local variance does not depend on the image location. An illustration is provided in Figure 3 (a) where the background parameters (c 1 , c 2 ) = (0.7, 0.1) and the patch parameters (c 1 , c 2 ) = (0.4, 0.005).

The wavelet leader coefficients are estimated using a standard 2D-DWT with orthonormal tensor product Daubechies mother wavelets with 2 vanishing moments.

In our simulation, we have set Q = 2, λ = 1 and J = 3. For every component m ∈ {1, . . . M }, (v q,m ) 1≤q≤Q are initially chosen to be equally distributed between the minimum and maximum values of Y m . Then, an alternate minimization of ( 9) and re-estimation of the (v q,m ) 1≤q≤Q is performed 5 times.

Scalar vs. vectorial segmentation. The scalar segmentation of Y = h, originally envisaged for the analysis of piecewise monofractal processes [START_REF] Pustelnik | Local regularity, wavelet leaders and total variation based procedures for texture segmentation[END_REF], is illustrated in Figure 3-(c)). We observe that it yields poor results. Indeed, since the process is multifractal, every h of the support of D(h) is present in any open interval of the trajectory for finite resolution data. The local estimation of h is hence not very meaningful. In addition, since D A and D B overlap, the sole feature h can not permit to discriminate between Ω A and Ω B . This limitation thus shows the need of investigating multiscale quantities related to the multifractal spectrum, namely C 1 and C 2 . Such a strategy has been previously investigated in [START_REF] Frecon | Multivariate optimization for multifractal-based texture segmentation[END_REF] with a different algorithmic solution involving a different segmentation for each component m ∈ {1, . . . , M }.

Segmentation performance. For the same patch displayed in Figure 3 (b), we have considered eight configurations reported in Figure 4 and defined as follow. The top (resp. bottom) plots correspond to a large (resp. small) difference of c 1 between Ω A and Ω B . The left (resp. right) plots model large (resp. small) difference of c 2 between Ω A and Ω B . In addition, we have further investigated the impact of c 2 : results in orange correspond to two tight multifractal spectrum D A and D B (i.e., small c 2 on both Ω A and Ω B ) whereas the blue one models one tight D A (i.e., small c 2 on Ω A ) and a widespread D B (i.e., large c 2 on Ω B ). Estimation performance are quantified in terms of misclassified pixels percentage over 20 realizations for all the different inputs Y . Unsurprisingly, using Y = h always leads to poor segmentation performance. In addition, all these experiments reproduce the expected behavior that the larger is the difference of c 1 (resp. c 2 ) between Ω A and Ω B , the better are the segmentation performance associated to Y = C 1 (resp. Y = C 2 ). However, a closer inspection shows that a larger difference in c 1 does not necessarily lead to better segmentation results for Y = C 2 (see the orange line in both left plots). Overall, we observe that there is always an interest in combining the information of both C 1 and C 2 .

Finally, it is worth noticing that the scalar segmentation of h is only 2 times faster than the vectorial ones either based on C 1 , C 2 or (C 1 , C 2 ). Experimentally, it takes less than 1 minute per image of size N = 2 9 × 2 9 .

CONCLUSIONS AND PERSPECTIVES

In this work, we have designed an analysis procedure for dealing with piecewise multifractal processes. We have shown the need of considering multiscale quantities C 1 and C 2 rather than the sole Hölder exponent usually considered for piecewise monofractal processes. An efficient joint vectorial segmentation procedure is proposed and yields satisfactory performance when applied to (C 1 , C 2 ).

However, the proposed segmentation procedure rely on the strong assumption that θ q,1 -θ q+1,1 = . . . = θ q,Mθ q+1,M . Therefore, the performance are very sensitive to the arbitrary choice and order of the level sets v q,m . In order to alleviate this limitation, different regularization terms are under current investigation to provide a more flexible joint vectorial segmentation strategy.

Fig. 1 .

 1 Fig. 1. Examples of scale-free textures.

Fig. 2 .

 2 Fig. 2. waveleter-leader coefficients. The wavelet leader coefficient (black cross) is defined as the local supremum of all wavelet coefficients taken within a spatial neighborhood across all finer scales (red)

  for image labelling. In what follows, Y = (Y m ) 1≤m≤M ∈ R N ×M can either model h (i.e., M = 1), C 1 or C 2 (i.e., M = J), or their concatenation (C 1 , C 2 ) (i.e., M = 2J).

Fig. 3 .

 3 Fig. 3. Illustration of the segmentation results. The piecewise multifractal image X is presented in (a) and has been generated from the mask (b). The scalar segmentation result obtained from h is displayed in (c). The vectorial segmentation results based on C 1 , C 2 and (C 1 , C 2 ) are respectively presented in (d), (e) and (f).

Fig. 4 .

 4 Fig. 4. Segmentation performance. Comparison of the segmentation performance is quantified in terms of mean misclassification rate depending on the choice of the input Y represented along the x-axis (error bar indicate the standard deviation). Each plot represent two configurations (blue and orange) of c 1 and c 2 on Ω A and Ω B .