Member, IEEE Denis Kouamé

Mathieu Biard

Member, IEEE Jean-Marc Girault

Aurore Bleuzen

Adaptive AR and Neuro-fuzzy approaches : access to cerebral particle signatures

Keywords: AR, emboli, Doppler, Detection, False alarm, Model, Neuro-Fuzzy, Non detection, score

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

S TROKEis the third leading cause of death and morbidity in the western world. In some countries the burden of stroke alone is estimated to be not less than 3 billion dollars per annum. The major factor in the total cost of stroke is the degree of related disability. Thus successful strategies to reduce the impact of stroke through early diagnosis will not only improve the quality of life of patients afflicted with stroke, but will also be very cost effective in many countries.

TransCranial Doppler (TCD) monitoring of signals from microemboli, a sophisticated technique for differentiation of stroke etiology and for early detection of impending stroke in carotid endarterectomy, plays an important role in primary and secondary stroke prevention. An embolus or foreign particle freely moving in the blood stream can be at the origin of sudden obstruction of an artery. This is referred to as an embolism. The consequences of cerebral embolism may be particularly severe, including cerebral infarction. Depending on its origin, an embolus can be a red blood cells aggregate or a piece of fat, a gas bubble, or any other foreign body, carried by the blood stream.

Detection of micro-emboli (small size emboli) and other emboli is important for several reasons such as preventing cerebrovascular accidents, finding the cause of embolism and validating the effectiveness of treatment. The main technique used to detect emboli is the recording of a transcranial ultrasound Doppler signal from the cerebral artery [START_REF] Spencer | Detection of middle cerebral artery emboli during carotid endarterectomy using transcranial doppler ultrasonography[END_REF]. The embolism signature in the blood stream is then assumed to be a non-predicted high intensity transient signal (HITS) superimposed on the Doppler signal backscattered by the blood. Most existing detection systems use an intensity measurement via the classical Fourier spectrogram, or any other time frequency distribution [START_REF] Cohen | Time-frequency distributions, a review[END_REF].

Commonly used embolus detection methods consist of comparing the ratio between blood energy distribution and that of the embolus with an empirical threshold [START_REF] Moehring | Pulse doppler ultrasound detection characterization and size estimation of emboli in flowing blood[END_REF], [START_REF] Moehring | Sizing emboli using pulse doppler ultrasound, part i: Verification of ebr model[END_REF]. As an embolus crossing the ultrasound beam can be heard, the first emboli detection methods were based on listening to ultrasound Doppler signals. These methods work very well and even outperform the automatic detection technique based on signal power distribution evaluation [START_REF] Van Zuilen | Automatic embolus detection compared with human experts : A doppler ultrasound study[END_REF] when dealing with large size (or quite audible) emboli. However, the risk of mistakes due to subjective interpretation in cases of micro-emboli led to investigation of more reliable methods which work when correctly, provided artifacts (which look like an embolus signature) are effectively rejected. Artifact rejection systems consist of either device-related systems (such as multiple gates or a multi-frequency system) or software-related systems, [START_REF] Brucher | Automatic embolus detection and artifact rejection with the first multifrequnecy transcranial doppler[END_REF], [START_REF] Devuyst | The matching pursuit : A new method of characterizing microembolic signal[END_REF]. Standard time-frequency distributions (TFD) were therefore studied for detection; the first being classical short time Fourier analysis which, although simple, is unable to distinguish small emboli. Other TFD techniques such as Wigner Ville distribution and its variants [START_REF] Smith | Processing doppler ultrasound signals from blood-borne emboli[END_REF], [START_REF] Roy | Comparaison of time-frequency estimators for peripheral embolus detection[END_REF], and wavelet transform [START_REF] Guetbi | New emboli detection methods[END_REF], [START_REF] Lui | Wavelet analysis of embolic heart sound detected by precordial doppler ultrasound during continuous venous air embolism in dogs[END_REF], [START_REF] Aydin | The use of the wavelet transform to describe embolic signals[END_REF], [START_REF] Krongold | Time-scale detection of microemboli in flowing blood with doppler ultrasound[END_REF] were also investigated, with limited success due to their intrinsic properties.

However, it has been shown that parametric modeling such as AutoRegressive (AR) modeling, combined with a sudden change detection technique [START_REF] Girault | Micro-emboli detection:an ultrasound doppler signal processing view point[END_REF], is one of the most reliable approaches to detect micro-emboli automatically. At the same time, it has been shown that integrating human expertise into a detection procedure provides good detection of emboli [START_REF] Fan | Automated embolus identification using a rule-based expert system[END_REF].

The common drawbacks of the last two techniques are twofold. First they give rise to a binary decision (existence or absence of embolus in the measurement area), although the embolism phenomenon is still not fully understood by physicians and no other means (e.g ultrasound brain scanner) is simultaneously available to check the reality of what is detected. Secondly they are mutually exclusive, whereas their advantages could be combined in a single approach. Human expertise has been taken into account, for this purpose, by using fuzzy logic. Thus we consider an adaptive fuzzy approach in order to fit the variability of embolus signatures. Some attempts at embolus detection have been made using neural networks [START_REF] Kemeny | Automatic embolus detection compared with human experts : A doppler ultrasound study[END_REF]. This study only used a classical neural network system such as "black-box", with its classical drawbacks (namely high number of iterations). Moreover, one important question, which was neither specific to the study reported by Kemeny et al [START_REF] Kemeny | Automatic embolus detection compared with human experts : A doppler ultrasound study[END_REF] nor solved in any automatic embolus detection, is the intrinsic evaluation of the reliability of the detection method. Here, we propose instead the use of both fuzzy logic and the neural network for optimization of membership function via the so-called neuro-fuzzy approach. A particular feature of our investigations is that, in order to take into account the possible uncertainty regarding the event detected and to provide an evaluation of the reliability of the event detected, we introduced a parameter score ranging from 0 to 1 as output of our system. After briefly revisiting the AR technique in the context of micro-embolus detection and introducing adaptive fuzzy modeling, the experimental results are used for comparison between parametric modeling and the fuzzy approach.

II. BACKGROUND

A. Doppler signal and parametric autoregressive modeling

Ultrasound (US) Doppler is the main method used for embolus detection. This technique consists in transmitting by means of a transducer, an ultrasound wave with frequency f 0 , to a selected area of a cerebral artery. After reception and demodulation, the resulting (analytic) complex number Doppler signal, see e.g. [START_REF] Jensen | Estimation of blood velocities using ultrasound. A signal processing approach[END_REF], denoted here x(t) may be defined by :

x(t) = I(t) + jQ(t) = K(t) exp [j(2πf d (t) + φ(t))] (1)
where j = √ -1, K(t) is the random magnitude of the Doppler signal depending on the characteristics of the transducers and the measurement area, and ϕ(t) is the random phase depending on the positions of the scatterers of the flowing fluid. The frequency shift f d , called Doppler frequency, is proportional to f 0 .

This Doppler signal contains a large amount of information concerning the moving target.

Unlike commonly used methods for embolus detection, parametric autoregressive (AR) modeling consists of working on a model of the signal and not directly on the signal. We briefly recall the classical results below. Considering a discrete time complex Doppler signal x(n) = x(t = nT e ) sampled at frequency 1/T e and assuming that it is the output of a p-order Autoregressive (AR) model, it can be expressed by :

x(n) = -a 1 (n)x(n-1)-a 2 (n)x(n-2)-...-a p (n)x(n-p)+η(n)
(2) where a i (n) are complex coefficients defining the AR model, p is order of the model ( number of coefficients) and η(n) is a complex white noise. For convenience, the previous equation is commonly expressed in matrix form as :

x(n) = ϕ T (n)θ(n) + η(n) (3)
where

ϕ T (n) = [-x(n -1), ..., -x(n -p)],
and

θ T (n) = [a 1 (n), ..., a p (n)].
It is well known that modeling the signal x as an AR process thus corresponds to obtaining the vector θ(n) from x which is an estimate of the vector θ(n). This estimation can be performed using, for example, the Recursive Least Squares (RLS) algorithm. Let

e(n) = x(n) -x(n)
be the prediction error. When the model ( 3) holds, the prediction error tends in probability, to white noise, when the number of observation data extends to infinity. Since the autocorrelation function (AF) of a white noise equals zero at any lag, except for the initial lag (n=0), the AF of prediction error therefore provides interesting Decision Information (DI), i.e information containing the embolus signature for this parametric method. Indeed, when an embolus crosses the sample volume, the prediction error will no longer be a white noise, and its AF at lag 1 (for example) will no longer be zero. The AF at lag 1 can be expressed by:

C n def = C(1) = 1 n n k=1 e(k)e(k -1) (4)
This can be estimated recursively each time n using a forgetting factor α (0 < α ≤ 1) by :

C n = αC n-1 + (1 -α)e(n)e(n -1) (5)
Here α = 0.9. Due to the above, C n will be almost zero for a normal Doppler signal, and the presence of an embolus will be characterized by a abrupt change. Therefore, to detect an embolus, we have to construct DI. Thus here :

DI =| C n |
The probability density function p(y) of this DI is [START_REF] Kouamé | Reliability evaluation of emboli detection using a statistical approach[END_REF] :

p(y) = 1 π ∞ 0 exp[-| y | cosh(u)]du , y = 0 (6)
An example of an embolus recorded on a patient is shown in fig. 1 to illustrate this.

In order to evaluate the reliability of embolus detection, such detection is performed in the framework of hypothesis testing. , where y is a particular value of the random variable Y . Denoting P (D i |H j ), the probability of deciding D i when H j is true, it follows that,

P ij def = P (D i |H j ) = Oi f Y |Hj (y|H j )dy. (7)
With these definitions we have

P F A = P 10 = 1 -P 00 (8) P N D = P 01 = 1 -P 11 ( 9
)
where PFA is the probability of false alarm and PND is the probability of non-detection evaluated using eq.( 8) and [START_REF] Roy | Comparaison of time-frequency estimators for peripheral embolus detection[END_REF]. In practice each hypothesis is characterized by decision information (DI), and belonging to O 0 or O 1 is represented by a threshold (λ). PND and PFA can thus be obtained by inverting eqs.( 8) and [START_REF] Roy | Comparaison of time-frequency estimators for peripheral embolus detection[END_REF]. PFA and PND are measures of reliability.

Finally the well known power spectrum density expression, from eq.( 3), denoted P x (n, f ) to account for its timedependent behavior, can be obtained for each time n as :

P x (n, f ) = σ 2 | 1 + p k=1 a k (n) exp(-2πjf k) | 2 ; (10)
where σ 2 is the power of the noise η and f is normalized frequency -0.5 ≤ f ≤ 0.5 and j = √ -1.

B. Neuro-Fuzzy approach

In this section, we also provide a brief introduction to the concepts used to study a problem with the neuro-fuzzy approach. Introduced in 1965 by Zadeh [START_REF] Zadeh | Fuzzy set[END_REF], the fuzzy approach is based on fuzzy reasoning or approximate reasoning, which is an inference procedure used to derive conclusions from a set of fuzzy if-then rules, as follows : "IF conditions THEN conclusion". Conditions and conclusion are of the type :" x is A", "y is B", respectively and x and y are variables representing for example input and output of the system under consideration. A and B are referred to as linguistic terms such as for instance "LARGE", "LOW",... and are characterized by membership functions µ A and µ B . For a particular value x 0 of x, it can be said that "x is A" with a degree of truth µ A (x 0 ). An important characteristic of such a system is that, due to its structure, it is immediately possible to insert human expertise through the rules. Direct use of this fuzzy reasoning may need a high number of cumbersome manual settings in complex problems. Neuro-fuzzy modeling was introduced [START_REF] Lin | Neural-network-based fuzzy logic control and decision system[END_REF], [START_REF] Jang | Anfis : Adaptive-network-based fuzzy inference system[END_REF] to account for this and fit the possible change in the system. Using our application, we consider here the case of the Sugeno model. That means that conclusions of rules are crisp linear combinations of variables. Many choices are possible [START_REF] Jang | Anfis : Adaptive-network-based fuzzy inference system[END_REF]. Our choice is explained below.

For the purpose of illustration, consider the example of two rules (R1 and R2), two inputs (x 1 and x 2 ) and one output (z) below. R1 :

IF x 1 is A 11 and x 2 is A 12 THEN f 1 = p 1 x 1 + q 1 x 2 + r 1 R2 : IF x 1 is A 21 and x 2 is A 22 THEN f 2 = p 2 x 1 + q 2 x 2 + r 2
The resulting output is z = µ1f1+µ2f2 µ1+µ2 ; where µ i = µ Ai1 (x 1 )× µ Ai2 (x 2 ) with i=1,2, as shown in (fig. 2-a).

The symbol × may be either a product symbol or any other T-norm symbol [START_REF] Zadeh | Fuzzy set[END_REF], [START_REF] Lin | Neural-network-based fuzzy logic control and decision system[END_REF], [START_REF] Jang | Anfis : Adaptive-network-based fuzzy inference system[END_REF]. T-norm is the operation, performed in order to provide the value of the output membership function value, when different fuzzy values, say µ i , i = 1, 2, .. result from different input membership functions (see example below). The T-norm used here is a Min-Max operator. This system can be modeled as a four layer neural network (fig. 2-b).

This can be generalized for a system with N rules as in (fig. 3). Thus a general neuro-fuzzy system is an equivalent four layer network (fig. 3) for a fuzzy system of which the ith rule is : N inputs four layers neural network system.

Ri : IF x 1 is A 1i and ... x n is A ni THEN f i = p ni x n + ... + p 1i x 1 + p 0i

Sugeno Fuzzy inference system with two rules (a) and its equivalent neural net model (b).

where "x i is A ij " is evaluated by µ ij (x i ). µ ij is a the membership function which is typically gaussian (of course other types can be used) with mean a ij and variance b ij .

Given a set of rules, the neuro-fuzzy technique adjusts the parameters of the system under consideration through the four layers defined as : Layer 1. This is the input layer. Inputs are x i , i = 1...n Layer 2. Each node corresponds to evaluation of degrees of truth. o , where f j is the consequent part of the j th rule : f j = p nj x n + ... + p 1j x 1 + p 0j a ij and b ij are referred to as premise parameters and p ij as consequent parameters. Basically, estimating and adjusting parameters is based on a hybrid algorithm in which the consequent parameters are identified by the least squares method. Many alternative algorithms have been developed to speed up the convergence, [START_REF] Onder Efe | A novel optimization procedure for training of fuzzy inference systems by combining variable structures systems technique and levenberg marquardt algorithm[END_REF], [START_REF] Hagan | Training feedforward networks with markardt algorithm[END_REF], [START_REF] Lera | Neighborhood based levenberg marquardt algorithm for neural network training[END_REF], [START_REF] Palit | Efficient training algorithm for takagisugeno type neuro-fuzzy network[END_REF]. In order to speed up the convergence rate and mostly to overcome the problem of selection of number of rules not solved by these techniques, we used the δ-operator Recursive Least Squares method [START_REF] Kouamé | Parametric estimation: Improvement of the rls approach using a differential approach[END_REF], [START_REF] Kouamé | Delta high order cumulant-based recursive instrumental variable algorithm[END_REF]. The problem of selection of number of rules, which is opened even in the fast algorithm, is solved here by using the matrix decomposition U DV H technique [START_REF] Niu | A factored form of instrumental variable algorithm[END_REF].

(2) ij = µ ij (x i ) = exp(-( xi-aij bij ) 2 ); i = 1, ...n, j = 1, ...
Due to the complexity of our application, and even if it is certain that the high intensity transient signal (HITS) detected is an embolus signature, in practice, there is no absolute guarantee that it is in fact an embolus signature. This is true both for expert physician and for all automatic detection systems.

To take this into account, we decided to provide a PFA for the parametric method for each HITS detected. For the neuro-fuzzy method, instead of giving a binary decision (absence/presence of embolus) as output of our system we gave a measure of detection, i.e. a score between 0 and 1.

This simple idea is, to our knowledge, new. With conventional approaches used for embolus detection, it is not possible to provide such an output (which is binary in these conventional approaches). The score output we propose here -a kind of "re-fuzzification" of the output -is both a decision and self evaluation of the decision.

This thus means that a score above 0.7 reveals the presence of an embolus, with a high level of reliability (70%), and a score below 0.3 reveals absence of embolus. We thus have an uncertainty range between 0.3 and 0.7 which has to be accounted for the physician.

III. APPLICATION

In order to validate the above techniques, we used experimental in vitro data obtained using "blood mimicking fluid" (BMF) embedded in a "Tygon" tube and circulated by a pump. Emboli were simulated by acrylic particles of different sizes (diameters from 200µm to 400µm). The BMF had acoustic properties similar to those of blood. Doppler signals were recorded using a Transcranial Doppler system, (WAKI 2 from ATYS MEDICAl).

Signals were recorded with the following experimental protocol:

* probe frequency: 2M Hz * Pulse Repetition Frequency (PRF) = 6KHz * The Tygon tube was located so that explored depth was 54mm. * burst length: 8.3 mm * sample length: 4.59 mm The inputs of the neuro-fuzzy approach were the estimated neuro-fuzzy parameters and the characteristics of the signal. For example, for a single gate system these characteristics were:

-DI (Decision information, see section II-A), with a range from 0 to 50 dB above the detection threshold. -(HITS) duration, (range 0 to 300ms). -Amaxn/Amaxp, i.e. the ratio between maximum power spectrum density in the domain of negative and positive frequencies, respectively, (range 0 to 50 dB).

f maxn, i.e. the normalized frequency of the maximum power spectrum density in the domain of negative frequencies, (range -0.5 to 0). -f maxp, i.e. the normalized frequency of the maximum power spectrum density in the domain of positive frequencies, (range 0 to 0.5). It is important to note here that all these parameters were computed from the parametric model, section II-A. A gaussian membership function was used for each of these parameters .

The output of the system, in this approach was the score. The inputs of the parametric method were only the above characteristics of the Doppler signal and the output was the decision.

Here are some additional practical details about our neuro-fuzzy system. The "expert knowledge" is available knowledge on emboli. This knowledge helps to choose the input parameters and relevant rules. For example most of the artifacts are bi-directional that is, they have frequency components in both positive and negative frequency domains. So inputs "fmaxn" and "fmaxp" were useful for dealing with direction of detected events. Here are the details of our system. 1) We recorded a set of artifact signals (29 in the example shown in fig. 4) and another set of acrylic signals (12 in example shown in fig. 4).

Score evaluated using neuro-fuzzy methods. (a) training data (40 signals). (b) Whole data (130 signals).

For each of these 41 signals we had five input parameters and one linear output (score). The artifact signals output scores were set to 0 and those of acrylic were set to 1. So we had an input matrix X of 41x5 values and an output vector Y of 41x1 values consisting of relevant zeros and ones. The number of rules per input was estimated. Here we had for DI, HITS duration, Amaxn/Amaxp, fmaxn, fmaxn respectively 3,3,3,2,2.

Then we chose gaussian membership functions. And from the above inputs and outputs, the membership parameters and output parameter were estimated using the hybrid algorithm. For illustration, the membership functions of Input "Amxn/Amaxp" are shown in fig. 5 In fig. 5 (a) and (b) are shown respectively initial membership functions and membership functions after optimization. When these parameters were estimated and the system was asked to find the scores for the above signals, the answer was of course correct, that is, 0 for artifact and 1 for acrylic signals.

2) Then from the whole data (including the training data), five input parameters per signal are computed, and the relevant scores were estimated.

IV. RESULTS

Several signals were recorded consisting of different types of artifact, together with acrylic particles. A typical in vitro circulating acrylic Doppler signal is shown in fig. 6, with its relevant DI.

A typical in vivo artifact signal is shown in fig. 7.

Two sizes (diameters) of acrylic particle (240µm and 300µm) were used first.

A set of 130 signals consisting of a mixture of artifacts (slight taps on the transducers) and acrylic particles were Our system was trained only once throughout. Then to test the generalization of the model learnt, 1,000 independent signals consisting of acrylic particles of different sizes (diameters 200µm, 230µm, 230µm, 260µm, 280µm, 300µm) and artifacts of different types were recorded and presented to the system. The results shown in fig. 4 are representative of the different tests made with acrylic particles of different sizes. Artifacts were detected at scores close to zero. The scores for acrylic particles were always higher than 0.7, except for fewer than 5% of cases.

We cannot in the strict sense talk about PFA and PND, however for purpose of comparison we defined for neuro-fuzzy approach: P F A = N F A N R where N F A was number of artifacts detected with score > 0.2 and N R was the number of tests. P N D = N N D N R where N N D was number of emboli (acrylic particles) detected with score < 0.7 and NR was the number of tests.

These led to P F A ≃ 0 and P N D ≤ 5%. This is summarized in Table.I. This means that this algorithm missed less than 5% of the acrylic signal data. In all the cases, when no acrylic was present, the decision provided by the algorithm was correct.

The results for the parametric method are summarized in Table .II. PFA was less than 5% and PND was less than 7% for a given detection threshold. PFA and PND were defined as follows : P F A = N F A N R where N F A was number of false alarm obtained for N R tests. P N D = N N D N R where N N D was number of non detections obtained for N R tests. These estimations of PFA and PND provide results similar to direct computations [START_REF] Kouamé | Reliability evaluation of emboli detection using a statistical approach[END_REF].

This means that the algorithm detected something (acrylic particles) whereas nothing was present in the signal in less than 5% of the (1,000) data used. This algorithm missed less than 7% of the acrylic signal data.

It is important to note that, as for most conventional systems, artifact rejection was performed by another procedure. Two main procedures are usually used. The first uses two or more gates (measurement areas), one of which can be located at an area where it is unlikely to detect an embolus. If the event detected is present simultaneously at all the gates, it is an artifact. The second is based on estimating the direction of the detected event. If the detected event is bi-directional then it is an artifact.

The performance of the two embolus detection techniques presented was very similar. The advantage of our neuro-fuzzy technique was that processing of artifacts was immediate, and only one gate was needed to perform embolus detection and artifact rejection. Thus this detection system does not require additional gates as in conventional detection systems. In vivo validation of this system is now being performed.

V. CONCLUSION

A specific neuro-fuzzy approach is proposed here in the framework of embolus detection. This technique was compared to automated embolus detection based on a parametric AR method using in vitro data. Although the performance of the two techniques was similar, the neuro-fuzzy technique presents the advantage of being able to perform detection using only one gate. This technique is thus a promising method to detect emboli effectively with a low cost system. In vivo validation of this system is currently being performed.

Jean-Marc Girault Jean-Marc Girault (S'98 M'00) received the Master degree of signal processing, biological and medical imaging from the University of Angers, France, in 1996 and the Ph.D. degree in "sciences de l'ingénieur" from the University of Tours. He is lecturer at the University of Tours in signal processing in the I.U.T. of Blois. His research area concerns ultrasound signal processing in biomedical applications and particularly in emboli detection.

Aurore Bleuzen

Fig. 1 .

 1 Fig. 1. Typical in vivo circulating embolus signatures in symptomatic patients. (a) real part Doppler signal of a large size embolus, (b) and (c) are spectrogram and DI, respectively, (d) real part Doppler signal of a smaller size embolus, (e) and (f) are spectrogram and DI, respectively.

  Fig. 2.

  Fig. 3.

, m Layer 3 .ij ; i = 1 ,

 31 Each node performs involvements through the T-norm operation : o ...n, j = 1, ..., m Layer 4. This output node performs defuzzification : z =

  Fig. 4.

Fig. 5 .

 5 Fig. 5. Illustration of membership function behavior using Input "Amxn/Amaxp": (a) and (b) are respectively initial membership functions and membership function after optimization. m.function1, m.function2 and m.function3 are membership function # 1, 2 and 3

Fig. 6 .

 6 Fig. 6. Typical in vitro circulating acrylic (of 240µm) signature. (a) real part Doppler signal, (b) spectrogram, (c) DI.

Fig. 7 .

 7 Fig. 7. Typical in vivo artifact signature. (a) real part Doppler signal, (b) spectrogram, (c) DI.

  Aurore Bleuzen received the Ph. D. degree in biological and medicine images and signals in 2003 from the University of Tours, France. She works as radiological practitioner in the Technological Innovation Center (CIT) of ultrasound of the hospital center of Tours, France. Her research interest includes the applications of contrast agents in ultrasound, the emboli detection and the development of high resolution ultrasound system.

  .

				Amaxn/Amxp initial membership functions
	Degree of membership	0.2 0.4 0.6 0.8 1	m.fontion1		m.fontion2		m.fontion3
		0
			0.04	0.06	0.08	0.1	0.12	0.14	0.16
			Amaxn/Amxp membership functions after optimization
		1	m.fontion1
	Degree of membership	0.2 0.4 0.6 0.8		m.fontion2				m.fontion3
		0
			0.04	0.06	0.08	0.1	0.12	0.14	0.16
					input3

Manuscript received February , 2005; revised October, 2005. This work was supported by European Community (UMEDS Project) and French Government.The authors are with LUSSI FRE CNRS 2448 University of Tours, 10 boulevard Tonnellé F37032 Tours Cedex France. E-mail: kouame@univtours.fr

ACKNOWLEDGEMENT

The authors wish to thank anonymous reviewers for their careful readings.


Emilie Chouzenoux

Jean-Christophe Pesquet

Anisia Florescu

A STOCHASTIC 3MG ALGORITHM WITH APPLICATION TO 2D FILTER IDENTIFICATION

Keywords: stochastic approximation, optimization, subspace algorithms, memory gradient methods, descent methods, recursive algorithms, majorization-minimization, filter identification, Newton method, sparsity, machine learning, adaptive filtering

Stochastic optimization plays an important role in solving many problems encountered in machine learning or adaptive processing. In this context, the second-order statistics of the data are often unknown a priori or their direct computation is too intensive, and they have to be estimated on-line from the related signals. In the context of batch optimization of an objective function being the sum of a data fidelity term and a penalization (e.g. a sparsity promoting function), Majorize-Minimize (MM) subspace methods have recently attracted much interest since they are fast, highly flexible and effective in ensuring convergence. The goal of this paper is to show how these methods can be successfully extended to the case when the cost function is replaced by a sequence of stochastic approximations of it. Simulation results illustrate the good practical performance of the proposed MM Memory Gradient (3MG) algorithm when applied to 2D filter identification.

INTRODUCTION

We consider a sequence of random variables (X n , y n ) n 1 taking their values in R N ×Q × R Q , defined on a probability space (Ω, F, P). Our objective is to solve the following minimization problem:

minimize h∈R N F (h) (1)
where

(∀h ∈ R N ) F (h) = 1 2 E y n -X ⊤ n h 2 + Ψ(h). (2)
Throughout this paper, E(•) denotes the mathematical expectation, • is the Euclidean norm, and Ψ is a function from R N to R, which plays the role of a regularization function. In particular, this function may be useful to incorporate some prior knowledge about h, e.g. some sparsity requirement, possibly in some transformed domain. We assume here that the following wide-sense stationarity properties hold:

(∀n ∈ N * ) E( y n 2 ) = ̺ (3) E(X n y n ) = r (4) E(X n X ⊤ n ) = R (5)
where ̺ ∈]0, +∞[, r ∈ R N , and R ∈ R N ×N is a symmetric positive semi-definite matrix.

Many optimization algorithms can be devised to solve Problem [START_REF] Nocedal | Numerical Optimization[END_REF] depending on the assumptions made on Ψ [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Afonso | An augmented Lagrangian aproach to the constrained optimization formulation of imaging inverse problems[END_REF]. In this work, we will be interested in Majorize-Minimize (MM) subspace algorithms [START_REF] Zibulevsky | ℓ 2 -ℓ1 optimization in signal and image processing[END_REF]. These approaches proceed by building at each iteration a simple majorant (e.g. a quadratic majorant) of the cost-function, which is minimized in a subspace of low dimension. This subspace is often restricted to the gradient computed at the current iterate and to a memory part (e.g. the difference between the current iterate and a previous one). In a number of recent works [START_REF] Chouzenoux | A majorizeminimize subspace strategy for subspace optimization applied to image restoration[END_REF][START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ2-ℓ0 image regularization[END_REF][START_REF] Florescu | A majorize-minimize memory gradient method for complex-valued inverse problem[END_REF], these algorithms are shown to provide fast numerical solutions to optimization problems involving smooth functions, in particular in the case of large-scale problems. Note that, although our approach will assume that Ψ is a differentiable function, it has been shown that tight approximations of nonsmooth penalizations such as ℓ 1 (resp. ℓ 0 ) functions, namely ℓ 2 -ℓ 1 (resp. ℓ 2 -ℓ 0 ) functions, can be employed and are often quite effective in practice [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ2-ℓ0 image regularization[END_REF][START_REF] Florescu | A majorize-minimize memory gradient method for complex-valued inverse problem[END_REF]. Another advantage of the class of optimization methods under investigation is that their convergence can be established under some technical assumptions, even in the case when Ψ is a nonconvex function (see [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ2-ℓ0 image regularization[END_REF] for more details).

One of the difficulties encountered in machine learning or adaptive processing is that Problem (1) cannot be directly solved since the second-order statistical moments ̺, r and R are often unknown a priori or their direct computation is too intensive, and they have thus to be estimated on-line from the related time series. In the simple case when Ψ = 0, the classical Recursive Least Squares (RLS) algorithm can be used for this purpose [START_REF] Haykin | Adaptive Filter Theory[END_REF]. When Ψ is nonzero, stochastic approximation algorithms have been developed such as the celebrated stochastic gradient descent (SGD) algorithm [START_REF] Bottou | Stochastic learning[END_REF]. This algorithm has been at the origin of a tremendous amount of works. It is known to be robust and easy to implement, but its con-vergence speed may be relatively slow. Various extensions of this algorithm have been developed to alleviate this problem [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF][START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF], to make it adaptive, or to improve its performance when estimating sparse vectors [START_REF] Paleologu | Sparse adaptive filters for echo cancellation[END_REF][START_REF] Murakami | A sparse adaptive filtering using time-varying soft-thresholding techniques[END_REF]. When Ψ ∝ • 1 , an on-line variant of the RLS algorithm was designed in [START_REF] Angelosante | Online adaptive estimation of sparse signals: where RLS meets the ℓ1 -norm[END_REF] which relies on a coordinate descent approach.

Designing Majorize-Minimize optimization algorithms in a stochastic context constitutes a challenging task since most of the existing works have been focused on batch optimization procedures, and the related convergence proofs usually rely on deterministic tools. We can however mention a few recent works [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF][START_REF]Stochastic Majorization-Minimization algorithms for large-scale optimization[END_REF] where stochastic MM algorithms are investigated for general loss functions under an independence assumption on the involved random variables, but without introducing any search subspace. Works which are more closely related to ours are those based on Newton or quasi-Newton stochastic algorithms [START_REF] Birge | A stochastic Newton method for stochastic quadratic programs with recourse[END_REF][START_REF] Bordes | SGD-QN: Careful quasi-Newton stochastic gradient descent[END_REF][START_REF] Yu | A stochastic quasi-Newton method for online convex optimization[END_REF][START_REF] Byrd | A stochastic quasi-Newton method for large-scale optimization[END_REF], in particular the approaches in [START_REF] Yu | A stochastic quasi-Newton method for online convex optimization[END_REF][START_REF] Byrd | A stochastic quasi-Newton method for large-scale optimization[END_REF] provide extensions of BFGS algorithm, but proving the convergence of these algorithms requires some specific assumptions. Like BFGS algorithm, MM subspace methods use a memory of previous estimates so as to accelerate the convergence.

In Section 2, we show how Problem (1) can be reformulated in a learning context. The MM strategy which is proposed in this work is described in Section 3.1. In Section 3.2, we give the form of the resulting recursive algorithm and, in Section 3.3, we evaluate its computational complexity. In Section 4, we show the good performance of the proposed stochastic Majorize-Minimize Memory Gradient (3MG) algorithm for solving a two-dimensional filter identification problem. Some conclusions are drawn in Section 5.

PROBLEM FORMULATION

In a learning context, function F can be replaced by a sequence (F n ) n 1 of stochastic approximations of it, which are defined as: for every n ∈ N * ,

(∀h ∈ R N ) F n (h) = 1 2n n k=1 y k -X ⊤ k h 2 + Ψ(h) = 1 2 ρ n -r ⊤ n h + 1 2 h ⊤ R n h + Ψ(h) (6)
where ρ n , r n , and R n are the following classical sample estimates of ̺, r, and R:

ρ n = 1 n n k=1 y k 2 (7) r n = 1 n n k=1 X k y k (8) R n = 1 n n k=1 X k X ⊤ k . (9)
Our objective in the next section will be to propose an efficient method for minimizing F n , for every n ∈ N * .

PROPOSED METHOD

Majorization property

At each iteration n ∈ N * , we propose to replace F n by a surrogate function Θ n (•, h n ) based on the current estimate h n (computed at the previous iteration). More precisely, a tangent majorant function is chosen such that

(∀h ∈ R N ) F n (h) Θ n (h, h n ) (10) F n (h n ) = Θ n (h n , h n ). ( 11
)
For the so-defined MM strategy to be worthwhile, the surrogate function has to be built in such a way that its minimization is simple. For this purpose, we will assume that the regularization function Ψ has the following form:

(∀h ∈ R N ) Ψ(h) = 1 2 h ⊤ V 0 h-v ⊤ 0 h+ S s=1 ψ s ( V s h-v s ) (12) where v 0 ∈ R N , V 0 ∈ R N ×N
is a symmetric positive semidefinite matrix, and, for every s ∈ {1, . . . , S}, v s ∈ R Ps , V s ∈ R Ps×N , and ψ s : R → R. In addition, the following assumptions will be made:

Assumption 1. (i) R + V 0 is a positive definite matrix.
(ii) For every s ∈ {1, . . . , S}, ψ s is a lower-bounded differentiable function and lim t→0 t =0 ψs (t)/t ∈ R, where ψs denotes the derivative of ψ s .

(iii) For every s ∈ {1, . . . , S}, ψ s ( √ .) is concave on

[0, +∞[. (iv) There exists υ ∈ [0, +∞[ such that (∀s ∈ {1, . . . , S}) (∀t ∈]0, +∞[) 0 υ s (t) υ, where (∀t ∈ [0, +∞[) υ s (t) = ψs (t)/t. 1
These assumptions are satisfied by a wide class of functions Ψ, in particular quadratic regularization functions, ℓ 2 -ℓ 1 functions, and various forms of smooth ℓ 2 -ℓ 0 functions [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ2-ℓ0 image regularization[END_REF].

Note that, for every n ∈ N * , the gradient of F n is given by

(∀h ∈ R N ) ∇F n (h) = A n (h)h -c n (h) (13)
1 The function is extended by continuity when t = 0.

where

A n (h) = R n + V 0 + V ⊤ Diag b(h) V ∈ R N ×N (14) c n (h) = r n + v 0 + V ⊤ Diag b(h) v ∈ R N (15) V = [V ⊤ 1 . . . V ⊤ S ] ⊤ ∈ R P ×N (16) v = [v ⊤ 1 . . . v ⊤ S ] ⊤ ∈ R P (17)
with

P = P 1 + • • • + P S , and b(h) = b i (h) 1 i P ∈ R P is such that (∀s ∈ {1, . . . , S}) (∀p ∈ {1, . . . , P s }) b P1+•••+Ps-1+p (h) = υ s ( V s h -v s ). (18)
We have then the following result:

Proposition 1. Under Assumptions 1(ii)-1(iv), for every n ∈ N * and h ∈ R N , a tangent majorant of F n at h is (∀h ′ ∈ R N ) Θ n (h ′ , h) = F n (h) + ∇F n (h) ⊤ (h ′ -h) + 1 2 (h ′ -h) ⊤ A n (h)(h ′ -h) ( 19
)
where A n (h) is given by [START_REF] Angelosante | Online adaptive estimation of sparse signals: where RLS meets the ℓ1 -norm[END_REF].

The proposed MM subspace algorithm consists of defining the following sequence of random vectors (h n ) n 1 :

(∀n ∈ N * ) h n+1 ∈ arg min h∈span Dn Θ n (h, h n ) (20)
where span D n is the vector subspace delineated by the columns of matrix D n ∈ R N ×Mn , and h 1 has to be set to an initial value. For example, we can choose, for every n ∈ N * ,

D n = [-∇F n (h n ), h n , h n -h n-1 ] if n > 1 [-∇F n (h 1 ), h 1 ] if n = 1 (21)
which yields the 3MG algorithm. Note that a similar choice of subspace can be found in optimization algorithms such as TWIST [START_REF] Bioucas-Dias | A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF]. A common assumption for subspace algorithms which will be adopted subsequently is that ∇F n (h n ) belongs to span D n .

Recursive MM strategy

By setting, for every n ∈ N, h n+1 = D n u n where u n is an R Mn -valued random vector, we deduce from ( 13), ( 19) and ( 20) that

u n = B † n D ⊤ n A n (h n )h n -∇F n (h n ) = B † n D ⊤ n c n (h n ) (22)
where

B n = D ⊤ n A n (h n )D n (23)
and (•) † is the pseudo-inverse operation. It is important to note that, as B n is of dimension M n ×M n where M n is small (typically M n = 3), this pseudo-inversion is not costly. This constitutes the main advantage of the proposed approach.

Let us now introduce the intermediate variables:

(∀n ∈ N * ) D R n = R n D n ∈ R N ×Mn (24) D V0 n = V 0 D n ∈ R N ×Mn (25) D V n = V D n ∈ R P ×Mn (26) D A n = A n+1 (h n+1 ) D n ∈ R N ×Mn . ( 27
)
By using ( 8), ( 9), (13) ( 14), ( 15), ( 22), ( 23), and by performing recursive updates of (r n ) n 1 and (R n ) n 1 , Algorithm 1 is obtained.

Algorithm 1 Stochastic MM subspace method

r 0 = 0, R 0 = 0 Initialize D 0 , u 0 h 1 = D 0 u 0 , D R 0 = 0, D V0 0 = V 0 D 0 , D V 0 = V D 0 For all n = 1, . . .                     r n = r n-1 + 1 n (X n y n -r n-1 ) c n (h n ) = r n + v 0 + V ⊤ Diag b(h n ) v D A n-1 = (1 -1 n )D R n-1 + 1 n X n (X ⊤ n D n-1 ) +D V0 n-1 + V ⊤ Diag b(h n ) D V n-1 ∇F n (h n ) = D A n-1 u n-1 -c n (h n ) R n = R n-1 + 1 n (X n X ⊤ n -R n-1 ) Set D n using ∇F n (h n ) D R n = R n D n , D V0 n = V 0 D n , D V n = V D n B n = D ⊤ n D R n + D V0 n + V ⊤ Diag b(h n ) D V n u n = B † n D ⊤ n c n (h n ) h n+1 = D n u n

Complexity

Since M n is small, the complexity of a direct implementation of this algorithm, evaluated in terms of multiplications at iteration n, is of the order of

N P (3M n + 1) + N (4M n + Q)/2
when N is large. However, this complexity can be reduced if matrices V 0 or V have a specific structure. In particular, if they are null matrices, the algorithm has the same order of complexity as the classical recursive least squares algorithm. Since the criterion then reduces to a quadratic function, Sherman-Morrison-Woodbury formula can be used in order to calculate iteratively the minimizer on the whole space in an efficient manner. The computational complexity can also be reduced by taking advantage of the specific form of matrices (D n ) n 1 . For example, if the subspace is chosen according to [START_REF] Bioucas-Dias | A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF], for every n > 1,

D V n = [-V ∇F n (h n ), V h n , V h n -V h n-1 ]. (28)
On the other hand,

V h n = V D n-1 u n-1 = D V n-1 u n-1 , (29)
which shows that a recursive formula holds to compute the last two components of D V n in (28). The initial complexity of 3P N multiplications is thus reduced to (P + 3)N . Similar recursive procedures can be employed to compute (D V0 n ) n>1 and (D Rn n ) n>1 .

APPLICATION TO 2D FILTER IDENTIFICATION

Problem statement

We now demonstrate the efficiency of the proposed stochastic algorithm in a filter identification problem. Consider the following observation model:

y = S(h)x + w, (30)
where x ∈ R L and y ∈ R L represent the original and degraded version of a given image, h ∈ R N is the vectorized version of an unknown two-dimensional blur kernel, S is the linear operator which maps the kernel to its associated Hankel-block Hankel matrix form, and w ∈ R L represents a realization of an additive noise. When the images x and y are of very large scale, finding an estimate h ∈ R N of the blur kernel can be very memory consuming, and one can expect good estimation performance by learning the blur kernel through a sweep of blocks of the dataset. Let us denote by X ∈ R L×N the matrix such that S(h)x = Xh. Then, we propose to define h as a solution to [START_REF] Nocedal | Numerical Optimization[END_REF], where, for all n ∈ N * , y n ∈ R Q and X ⊤ n ∈ R Q×N , are subparts of y and X, respectively, corresponding to Q ∈ {1, . . . , L} lines of this vector/matrix. For the regularization term Ψ, we consider, for every s ∈ {1, . . . , N } (S = N ), an isotropic penalization on the gradient between neighboring coefficients of the blur kernel, i.e., P s = 2 and

V s = ∆ h s ∆ v s ⊤
, where ∆ h s ∈ R N (resp. ∆ v s ∈ R N ) is the horizontal (resp. vertical) gradient operator applied at pixel s. The smoothness of h is then enforced by choosing, for every s ∈ {1, . . . , S} and u ∈ R, ψ s (u) = λ 1 + u 2 /δ 2 with (λ, δ) ∈]0, +∞[ 2 . Finally, in order to guarantee the existence of a unique minimizer, the strong convexity of F is imposed by taking v 0 = 0 and V 0 = τ I N , where τ is a small positive value (typically τ = 10 -10 ).

Simulation results

The original image, presented in Figure 1(a), is the San Diego image, of size 1024 × 1024 pixels, available at http://sipi.usc.edu/database/. The original blur kernel h with size 21 × 21, and the resulting blurred image, which has been corrupted with a zero-mean Gaussan noise with standard deviation σ = 0.03 (blurred signal-to-noise ratio equal to 24.8 dB), are displayed in Figures 1(b)(c). Figure 1(d) presents the estimated kernel, using the proposed stochastic algorithm with the subspace given by [START_REF] Bioucas-Dias | A new TwIST: two-step iterative shrinkage/thresholding algorithms for image restoration[END_REF]. Parameters (λ, δ) were adjusted so as to minimize the normalized root mean square estimation error, here equal to 0.087. Figure 2 illustrates the variations of the estimation error with respect to the computation time for the proposed algorithm, the SGD algorithm with a decreasing stepsize proportional to n -1/2 , and the regularized dual averaging (RDA) method with a constant stepsize from [START_REF] Xiao | Dual averaging methods for regularized stochastic learning and online optimization[END_REF], when running tests on an Intel(R) Core(TM) i7-3520M @ 2.9GHz using a Matlab 7 implementation. Note that for the latter two algorithms, the stepsize parameter was optimized manually so as to obtain the best performance in terms of convergence speed. Finally, note that stochastic 3MG and RDA algorithms were observed to provide asymptotically the same estimation quality, whatever the size of the blocks. In this example, the best trade-off in terms of convergence speed is obtained for Q = 64 × 64.

CONCLUSION

In this work, we have proposed a stochastic MM Memory Gradient algorithm for on-line penalized least squares estimation problems. The method makes it possible to use large-size datasets the second-order moments of which are not known a (dashed-dotted red line) and RDA algorithm with constant stepsize (dashed blue line).

priori. We have shown that the proposed algorithm is of the same order of complexity as the classical RLS algorithm and that its computational cost can be even reduced by taking advantage of specific forms of the search subspace. The good numerical performance of the proposed algorithm has been demonstrated in the context of 2D filter identification for large size images. In our future work, a theoretical analysis of the convergence properties of the proposed method will be conducted. In addition, we plan to apply this technique to system identification or inverse modeling using adaptive filters.

Fig. 1 .

 1 Fig. 1. (a) Original image. (b) Blurred and noisy image. (c) Original blur kernel. (d) Estimated blur kernel, with relative error 0.087.

Fig. 2 .

 2 Fig. 2.Comparison of stochastic 3MG algorithm (solid black line), SGD algorithm with decreasing stepsize ∝ n 1/2 (dashed-dotted red line) and RDA algorithm with constant stepsize (dashed blue line).


F Abboud

E Chouzenoux

J.-C Pesquet

J.-H Chenot

L Laborelli

A HYBRID ALTERNATING PROXIMAL METHOD FOR BLIND VIDEO RESTORATION

Keywords: Blind deconvolution, video processing, optical flow, proximal algorithms, convex optimization, regularization

Old analog television sequences suffer from a number of degradations. Some of them can be modeled through convolution with a kernel and an additive noise term. In this work, we propose a new blind deconvolution algorithm for the restoration of such sequences based on a variational formulation of the problem. Our method accounts for motion between frames, while enforcing some level of temporal continuity through the use of a novel penalty function involving optical flow operators, in addition to an edge-preserving regularization. The optimization process is performed by a proximal alternating minimization scheme benefiting from theoretical convergence guarantees. Simulation results on synthetic and real video sequences confirm the effectiveness of our method.

INTRODUCTION

The visual contents of videos from the previous decades are often in a degraded state so that they no longer correspond to the actual public demand. The observed degradations result from several physical phenomena which happened during the sensing, transmission, and recording processes [START_REF] Chenot | Restoration of archived television programmes for digital broadcasting[END_REF][START_REF] Naranjo | Flicker reduction in old films[END_REF][START_REF] Pizurica | Noise reduction in video sequences using wavelet-domain and temporal filtering[END_REF][START_REF] Kokaram | MCMC for joint noise reduction and missing data treatment in degraded video[END_REF][START_REF] Kokaram | Motion Picture Restoration: Digital Algorithms for Artefact Suppression in Degraded Motion Picture Film and Video[END_REF][START_REF] Tendero | Blind uniform motion blur deconvolution for image bursts and video sequences based on sensor characteristics[END_REF]. Most of them can be summarized into two main categories. The first type of degradation is of random nature and appears in images as noise, mainly caused by electronic devices. The second one is deterministic and results in blur and oscillations, whose common causes are camera defocus, lens aberrations, relative camera-scene motion, diffusion in sensors, physical or electronic transmission problems. In this paper, we address the problem of restoring a video sequence (y t ) 1 t T , composed of T successive frames having N pixels, related to an original unknown sequence (x t ) 1 t T through the degradation model

(∀t ∈ {1, . . . , T }) y t = S(h)x t + w t , (1)
where h ∈ R P is the vectorized version of an unknown convolution kernel, S is the linear operator which maps the kernel to its associated Hankel-block Hankel matrix form, and w t ∈ R N represents an unknown additive noise. The objective of video restoration is to recover

(x t ) 1 t T from (y t ) 1 t T ,
The research leading to these results has received funding from the EC Seventh Framework Programme (FP7/2007-2013) under Grant Agreement n o 600827 while h is unknown. When T = 1, one recovers the well known still image blind deconvolution problem, that has received much attention [START_REF] Almeida | Multi-resolution approach for blind deblurring of natural images[END_REF][START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF][START_REF] Bolte | Alternating proximal algorithm for blind image recovery[END_REF][START_REF] Komodakis | MRF-based blind image deconvolution[END_REF][START_REF] Levin | Understanding and evaluating blind deconvolution algorithms[END_REF][START_REF] Krishnan | Blind deconvolution using a normalized sparsity measure[END_REF]. However, the specific multiframe case when T > 1 remains scarcely studied in the literature. Among existing works, variational approaches based on alternating minimization strategies have been proposed in [START_REF] Soulez | Blind deconvolution of video sequences[END_REF][START_REF] Sroubek | Superresolution and blind deconvolution of video[END_REF]. Both methods deal with the case of a time-varying positive blur operator, possibly combined with a decimation [START_REF] Sroubek | Superresolution and blind deconvolution of video[END_REF]. The video restoration process then consists of minimizing a least-squares criterion augmented with regularization terms on the sought images and kernels. In [START_REF] Soulez | Blind deconvolution of video sequences[END_REF], the time variations are assumed to be small, which allows the use of penalty functions on the difference of consecutive frames, in addition to a total variation-like prior on each image. Moreover, a parametric model is used on the convolution operator, so as to enforce a specific form of the kernel, defined a priori (e.g. a Lorentzian function). In [START_REF] Sroubek | Superresolution and blind deconvolution of video[END_REF], motion effects are taken into account through a modification of the observation model, while quadratic penalty functions are chosen for the images and the kernels with the aim to simplify the minimization process. It is also worth mentioning the work in [START_REF] Sroubek | Robust multichannel blind deconvolution via fast alternating minimization[END_REF] which deals with a variant of Problem (1) where the convolution kernel is time-varying while the same scene is observed (up to some motion effects) in a given group of frames.

The contributions of this paper are (i) a versatile formulation of the blind video sequence restoration problem that includes an original strategy to account for motion effects, and (ii) the proposition of a novel efficient optimization algorithm with proven convergence to solve this problem. Our developments are grounded on the use of recent advances in nonsmooth optimization, and more specifically on a new hybrid proximal alternating method. The paper is organized as follows: Section 2 presents the penalized criterion we minimize. In Section 3 we introduce the proposed algorithm. Then, Section 4 illustrates the applicability of our method through a set of experiments for video archive contents restoration.

PROBLEM FORMULATION

Minimization problem

Let us define x = (x t ) 1 t T ∈ R T N . Estimates of h and x are obtained by solving the minimization problem:

minimize x∈R T N ,h∈R P (G(x, h) = Φ(x, h) + Θ(h) + Ψ(x)) . (2)
Hereabove, Φ is a data fidelity term, which is chosen in accordance with the noise statistics. In this paper, we consider the case when w = (w t ) 1 t T is an additive Gaussian noise, so that Φ is equal to a least squares data criterion, defined, for every x ∈ R T N and h ∈ R P as

Φ(x, h) = 1 2 T t=1 S(h)x t -y t 2 . (3)
Functions Θ and Ψ are regularizing terms acting, respectively, on the kernel and the image sequence. The design of suitable penalization functions plays a major role on the quality of the video restoration process. This will be discussed in the next two sections.

Regularization strategy for the kernel

Prior information on the kernel h is incorporated by choosing:

(∀h ∈ R P ) Θ(h) = ι H (h) + ρ h , (4)
where ρ > 0, and ι H is the indicator function of a closed convex set H, i.e.,

ι H (h) = 0 if h ∈ H +∞ otherwise. (5)
Here, the constrained set is the polyhedral set defined as

H = h = (h p ) 1 p P ∈ R P : P p=1 h p = 1, (∀p ∈ {1, . . . , P }) h min,p h p h max,p . (6)

Regularization strategy for the video sequence

We consider a composite penalization function on the sequence x, whose first term introduces an a priori knowledge on each image independently, while the second term penalizes the whole sequence, i.e.

Ψ(x) = T t=1 Ψ t (x t ) + M(x). (7)
For all t ∈ {1, . . . , T }, Ψ t is defined as

(∀x t ∈ R N ) Ψ t (x t ) = η tv(x t ) + ι [xmin,xmax] N (x t ), ( 8
)
where η is a positive regularization constant, and tv denotes the total variation function from [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. Furthermore, in order to favor the similarity between the successive video frames while taking into account existing motions, we define

M(x) = 1 2 T t=1 ∈Vt β ,t x t -M →t x 2 , (9)
where (β ,t ) ,t are positive weights, the index set V t defines the neighborhood of t (i.e. a set of indices ∈ V t for which | -t| is small), and M →t ∈ R N ×N is the linear operator modeling the motion compensation process between the reference frame t and a neighboring frame [START_REF] Liu | Human-assisted motion annotation[END_REF].

PROPOSED ALGORITHM

Optimization tools

Let us first recall the notion of proximity operator relative to a metric.

Definition 1 Let f : R N → (-∞, +∞] be a convex, proper, lower semicontinuous function, let U ∈ R N ×N be a symmetric positive definite matrix. For every x ∈ R N , the variational problem

minimize z∈R N f (z) + 1 2 x -z 2 U , (10)
admits a unique solution, which is denoted by prox U ,f (x).

The so-defined operator prox U ,f : R N → R N is the proximity operator of f relative to the metric induced by U .

Hereabove, • U denotes the weighted norm defined by . U = •, U • , where •, • is the usual Euclidean scalar product. Note that, if U is the identity matrix, one recovers the usual proximity operator prox f : R N → R N , which is at the core of numerous convex optimization algorithms (see [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Parikh | Proximal algorithms[END_REF] for a tutorial).

We also recall the following characterization of smooth functions:

Definition 2 Let f : R N → R be a differentiable function. Its gradient ∇f is said L-Lipschitzian continuous if, for every (x, z) ∈ (R N ) 2 , ∇f (x) -∇f (z) L x -z . (11)

Hybrid proximal alternating strategy

Although the objective function G is nonconvex because of the coupling in the data fidelity term (3), it remains yet convex with respect to each variable x 1 , . . . , x T and h, as the regularization terms are convex. A standard approach for solving (2) is thus to adopt an alternating minimization strategy, where, at each iteration, G is minimized with respect to one variable while the others remain fixed. This approach, sometimes referred to as Block Coordinate Descent [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF] or nonlinear Gauss-Seidel method [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF], has been widely used in the context of blind deconvolution [START_REF] Almeida | Multi-resolution approach for blind deblurring of natural images[END_REF][START_REF] Komodakis | MRF-based blind image deconvolution[END_REF][START_REF] Soulez | Blind deconvolution of video sequences[END_REF][START_REF] Sroubek | Superresolution and blind deconvolution of video[END_REF][START_REF] Sroubek | Robust multichannel blind deconvolution via fast alternating minimization[END_REF]. However, its convergence is only guaranteed under restrictive assumptions [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]. Therefore, alternative strategies based on proximal tools have been proposed in [START_REF] Bolte | Alternating proximal algorithm for blind image recovery[END_REF][START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF][START_REF] Xu | A block coordinate descent method for regularized multi-convex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Luo | On the convergence of the coordinate descent method for convex differentiable minimization[END_REF][START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems. An approach based on the Kurdyka-Łojasiewicz inequality[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]. They consist of replacing, at each iteration, the minimization step by either a (single) proximal step [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF][START_REF] Attouch | Proximal alternating minimization and projection methods for nonconvex problems. An approach based on the Kurdyka-Łojasiewicz inequality[END_REF] or a Forward-Backward step [START_REF] Xu | A block coordinate descent method for regularized multi-convex optimization with applications to nonnegative tensor factorization and completion[END_REF][START_REF] Luo | On the convergence of the coordinate descent method for convex differentiable minimization[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF][START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF], giving rise, respectively, to the so-called proximal (resp. forward-backward) alternating algorithms. In this work, we propose to adopt the following hybrid alternating algorithm:

Algorithm 1 Video blind deconvolution algorithm

For every k ∈ N, (λ 1 t , . . . , λ k T , µ k ) ∈ (0, +∞) T +1 . Initialize with x 0 1 , . . . , x 0 T ∈ R N , and h 0 ∈ R P . Iterations:

For k = 0, 1, . . .         For t = 1, . . . , T     xt,k = x k+1 1 , . . . , x k+1 t-1 , x k t , . . . , x k T , x k t = x k t -λ k t ∇ xt Φ(x t,k , h k ) + ∇ xt M(x t,k ) , x k+1 t = prox λ k t Ψt x k t , h k+1 = prox µ k (Θ+Φ(x k+1 ,•)) (h k ).
In the above algorithm, for every k ∈ N, for every t ∈ {1, . . . , T }, ∇ xt Φ(x t,k , h k ) is the partial gradient of Φ with respect to x t computed at (x t,k , h k ). The vector ∇ xt M(x t,k ) is defined in a similar way.

Convergence analysis

The convergence of Algorithm 1 requires the design of a proper strategy to determine the stepsize parameters λ k t 1 t T and µ k at each iteration k ∈ N. First, let us state two properties, which are related to the quadratic form of functions Φ and M.

Proposition 1 Let µ ∈ (0, +∞) and x ∈ R T N . Let us define the symmetric definite positive matrix

A(x) = T t=1 X t X t + µ -1 I P , ( 12
)
where I P is the identity matrix of R P , and, for every t ∈ {1, . . . , T },

X t ∈ R N ×P is such that S(h)x t = X t h. Then, for every h ∈ R P , prox µ(Θ+Φ(x,•)) (h) = prox A(x),Θ (h -A(x) -1 ∇ h Φ(x, h)), (13)
and, for every x ∈ R T N and h ∈ R P ,

Φ(x, h ) + ∇ h Φ(x , h ) (h -h ) + 1 2 h -h 2 A(x) Φ(x, h). ( 14
)
Proposition 2 Let t ∈ {1, . . . , T } and h ∈ R P . The partial gradient function x t → ∇ xt Φ(x, h) is Lipschitzian continuous, with Lipschitz modulus

L t (h) = S(h) 2 + ∈Vt β ,t + β t, M t→ 2 . ( 15
)
Combining Propositions 1, 2 and the convergence properties of the BC-VMFB (Block Coordinate Variable Metric Forward Backward) algorithm derived in [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF], the convergence of Algorithm 1 can be proved as expressed by the following result:

Theorem 1 Let (x k , h k ) k∈N be a sequence generated by Algorithm 1, where, for all k ∈ N, λ k t = θL t (h k ) -1 with θ ∈ (0, 2) and L t (•) defined in [START_REF] Sroubek | Robust multichannel blind deconvolution via fast alternating minimization[END_REF]. Then the sequence

(x k , h k ) k∈N converges to a critical point ( x, h) of G. More- over, (G(x k , h k )) k∈N is a nonincreasing sequence converg- ing to G( x, h).

Practical implementation

The proximal steps in Algorithm 1 are not explicit, and subiterations are needed. In practice, the proximity operator of Θ+ Φ(x, •) is computed with Dykstra algorithm [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF], while the Dual Forward-Backward algorithm [START_REF] Combettes | Proximity for sums of composite functions[END_REF] is used to compute the proximity operators of functions Ψ t , for t ∈ {1, . . . , T }. The proposed algorithm is thus used under an inexact form for which the conclusions of Theorem 1 still hold.

EXPERIMENTAL RESULTS

Synthetic data

We first demonstrate the practical performance of our method on a synthetic video restoration example. We extract T = 15 frames of size N = 256×256 from the video sequence Claire, available at http://www.cipr.rpi.edu/resource/ sequences/cif.html. To generate the observed video y = (y t ) 1 t T , we degraded the original one with the horizontal convolution kernel displayed in Figure 2, with size P = 41, which corresponds to a realistic model of those affecting old analog television sequences. The video is further corrupted with a white additive zero-mean Gaussian noise. Parameters ρ, η and (β ,t ) ,t were adjusted to maximize the SNR between the original and reconstructed sequences. Neighboring frames such that | -t| = 1 have been taken into account in the regularization term M. The motion matrices (M →t ) ,t have been estimated from the degraded sequence y, using the optical flow estimation algorithm from [START_REF] Liu | Human-assisted motion annotation[END_REF]. Subpixel motions have been ignored, so as to obtain simple expressions for M →t ,t and M →t 2 ,t . Algorithm 1 provides estimates ( x, h) displayed in Figures 12, with SNR = 32.2 dB, and relative quadratic error on the kernel equal to 0.006. To emphasize the efficiency of the motion-based regularization strategy, note that, when taking β ,t ≡ 0, we observe a decrease of the estimation quality, as the SNR of the restored sequence is equal to 31.2 dB and the kernel estimation error to 0.016.

Real data

We now apply our method to a real video sequence (for which no ground truth is available) containing T = 5 frames of N = 720 × 576 pixels, extracted from a French broadcast archive programme "Au théâtre ce soir". This sequence, provided by INA, results from a recording taken from a radio frequency analog terrestrial link, affected by multiple paths. Algorithm 1 is employed to restore the luminance component of the YCrCb representation of the sequence, the two chrominance components remaining untouched. To account for the interlacing, the odd and even lines of the video are restored separately, and then gathered at the end of the restoration process. Examples of degraded and restored images, as well as the estimated kernel, are displayed in Figures 345. One can observe on Figure 4 the ability of the proposed method to highly improve the sharpness of the image. Moreover, the knowledge of an estimation of the convolution kernel allows us to better analyze the main source of the video degradation process, here mainly due to oscillation effects probably caused by the analog transmission.

CONCLUSION

In this work, we have presented a new variational method for blind deconvolution of video sequences. Our approach relies on the minimization, through a proximal alternating algorithm, of a penalized criterion that accounts for the temporal continuity between the video frames. Our experimental results on both synthetic and real data showed that our method leads to well restored images, as well as satisfactory estimations of the blur kernel. It should be emphasized that the versatility of the proposed algorithm makes it possible to handle various forms of blur kernels, and to minimize convex and nonconvex data fidelity and penalty functions. Moreover, using variants of this algorithm could open the gate to faster implementation strategies based on parallel computing architectures.

Fig. 1 .

 1 Fig. 1. 3rd and 12th frames of the Claire sequence: Noisy convolved images, SNR = 25.6 dB (top) and restored images (bottom) with the proposed algorithm, SNR = 32.2 dB.

Fig. 2 .

 2 Fig. 2. Original kernel (black) and estimated kernel (red), with relative error 0.006.

Fig. 3 .

 3 Fig. 3. 2nd frame of the real INA sequence: Degraded image (top) and restored image (bottom) with the proposed algorithm.

Fig. 4 .

 4 Fig. 4. Zoom on a part of the above images.

Fig. 5 .

 5 Fig. 5. Estimated kernel


Y Marnissi

A Benazza-Benyahia

E Chouzenoux

J.-C Pesquet

MAJORIZE-MINIMIZE ADAPTED METROPOLIS-HASTINGS ALGORITHM. APPLICATION TO MULTICHANNEL IMAGE RECOVERY

Keywords: MCMC methods, Langevin diffusion, Majorize-Minimize, MMSE, multichannel image restoration

One challenging task in MCMC methods is the choice of the proposal density. It should ideally provide an accurate approximation of the target density with a low computational cost. In this paper, we are interested in Langevin diffusion where the proposal accounts for a directional component. We propose a novel method for tuning the related drift term. This term is preconditioned by an adaptive matrix based on a Majorize-Minimize strategy. This new procedure is shown to exhibit a good performance in a multispectral image restoration example.

INTRODUCTION

Recovering the signal of interest from degraded observations embedded in an additive noise is a key issue for many applications such as remote sensing imaging. In this respect, a Bayesian framework can be adopted to compute the Minimum Mean Squared Estimator (MMSE). However, it is not always possible to derive a closed expression of the posterior distribution involved in the MMSE. To alleviate this problem, Markov Chain Monte Carlo (MCMC) approaches have been developed. They consist of constructing an irreducible Markov chain whose stationary distribution is the unknown posterior distribution. Building the Markov chain corresponds to a specific way of exploring the state space. For this purpose, Metropolis-Hastings (MH) algorithms have been intensively used [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF]. The key issue however is the choice of a proposal density. Recent methods such as those based on Langevin diffusion have incorporated a directional component for the proposal [START_REF] Roberts | Exponential convergence of Langevin distributions and their discrete approximations[END_REF]. More precisely, two parameters (a stepsize and a scale matrix) are introduced to guide the directional component. The problem of setting the scale matrix must be carefully addressed, especially for high dimensional problems. Several solutions have been considered [START_REF] Roberts | Exponential convergence of Langevin distributions and their discrete approximations[END_REF][START_REF] Pillai | Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions[END_REF][START_REF] Stuart | Conditional path sampling of SDEs and the Langevin MCMC method[END_REF]. In this work, we propose a novel approach for choosing the scale matrix based on a Majorize-Minimize (MM) strategy. The paper is organized as follows: In Section 2, we formulate the problem and we give a brief overview of the Langevin diffusion process. In Section 3, we describe the new MM adapted MH algorithm. Section 4 is devoted to experimental results for multicomponent image restoration. Finally, some concluding remarks are drawn in Section 5

RELATED WORKS

Problem statement

In this paper, we address a wide array of problems where the vector x in R Q of samples of an unknown signal X is estimated from the observation vector z in R N given the following observation model:

z = Hx + w (1)
where the matrix H ∈ R N ×Q corresponds to a linear degradation operator, and w in R N is the vector of a Gaussian additive noise samples. In this work, we adopt a Bayesian framework and we aim at computing the MMSE x MMSE = E π X [x] where the posterior distribution π X is related to the prior distribution p X of the unknown vector X:

π X (x) ∝ p X (x)p(z|x). (2)
Generally, the computation of x MMSE involves integrals that are both analytically and numerically intractable. The Monte Carlo approach is a classical alternative solution which consists of simulating a sufficient number of i.i.d. random variables from the posterior distribution π X and approximating the MMSE estimator by the empirical average over all these samples. However, the target posterior is often complex and does not present a closed form, so that direct sampling is not always possible. To alleviate this difficulty, MCMC methods have been developed. They consist of building an irreducible Markov chain whose stationary distribution is π X . The asymptotic state of the chain is then considered as a sample of the target distribution. From an initial state x 0 , the problem reduces to exploring the state space according to the transition probabilities that characterize the Markov chain. Different ways of moving from a state to another have been reported. Among them, much attention was paid to MH algorithms that generate a random walk according to a pro-posal density and implement a method for rejecting proposed moves [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF].

Langevin diffusion

The choice of the proposal distribution is crucial as it impacts the statistical properties of the resulting Markov chain especially for complex and high-dimensional target distributions. Advanced MH methods introduce a directional component for the proposal. In this respect, Langevin diffusion strategies adjust the state transition by accounting for the gradient direction of the target density. We have thus, for every t ∈ N,

x t+1 = x t + ε 2 b(x t ) + ε σ(x t )n t+1 , (3)
where ε > 0 is the stepsize resulting from Euler's discretization of the diffusion, (n t ) t∈N are realizations of a zero-mean white noise, σ(x t ) is a positive definite matrix and, b(x) = (b i (x)) Q i=1 is a drift term. The latter is defined as follows:

b i (x) = 1 2 N j=1 A ij (x) ∂ log π X (x) ∂x j +|A(x)| 1 2 N j=1 ∂ ∂x j A ij (x)|A(x)| -1 2 , (4)
where A(x) = σ(x)σ ⊤ (x) and |A(x)| denotes the determinant of this matrix. It can be proved that the Langevin process has π X as its stationary distribution and is more likely to accept proposed values than a standard random walk. Indeed, the gradient information of the target distribution allows the chain to be guided toward regions of higher probability, where most of the samples should lie and hence, enables to achieve high acceptance rates. To this end, it is worth noting that the two scale parameters play an important role: ε determines the length of proposed jumps whereas A controls the direction. Three classes of algorithms have been developed from this diffusion depending on the choice of A.

Choice of the scale matrix

The standard Metropolis adjusted Langevin algorithm (MALA) is the simplest form of this diffusion when A equals I Q , the identity matrix of R Q [2]:

(∀t ∈ N) x t+1 = x t + ε 2 2 ∇ log π X (x t ) + ε n t+1 . (5)
However, it should be noted that a bad adjustment of ε can significantly affect the convergence rate especially for large scale problems [START_REF] Pillai | Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions[END_REF]. For this reason, many methods focus on how to choose a suitable stepsize such that the asymptotic average acceptance rate is bounded away from zero for high dimensions [START_REF] Pillai | Optimal scaling and diffusion limits for the Langevin algorithm in high dimensions[END_REF][START_REF] Roberts | Optimal scaling of discrete approximations to Langevin diffusions[END_REF]. Another approach consists of accelerating the algorithm by preconditioning the proposal density using a given constant scale matrix [START_REF] Stuart | Conditional path sampling of SDEs and the Langevin MCMC method[END_REF]. However, there is no clear guiding strategies for the selection of such a constant matrix.

Recent algorithms [START_REF] Martin | A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion[END_REF][START_REF] Zhang | Quasi-Newton methods for Markov chain Monte Carlo[END_REF][START_REF] Bui-Thanh | A scaled stochastic Newton algorithm for Markov chain Monte Carlo simulations[END_REF][START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF] propose adaptive procedures where A is tuned according to the past behavior of the Markov chain resorting to some deterministic optimization tools. For example, when setting A to the inverse of the Hessian matrix of log π X and, assuming a locally constant curvature, the term involving the derivatives of the scale matrix in (4) reduces to zero. Consequently, the computation of the drift term b becomes a scaled Newton step for minimizinglog π X . Thus, a new sample of the Newton-based MCMC is more likely drawn from a highly probable region and then more likely accepted, which can speed up the convergence of the sampling process [START_REF] Martin | A stochastic Newton MCMC method for large-scale statistical inverse problems with application to seismic inversion[END_REF][START_REF] Zhang | Quasi-Newton methods for Markov chain Monte Carlo[END_REF][START_REF] Bui-Thanh | A scaled stochastic Newton algorithm for Markov chain Monte Carlo simulations[END_REF]. However, in practice, this method has a high computational load since it requires the computation of the Hessian matrix and its inverse at each iteration. This is particularly critical for large scale problems and/or when the Hessian matrix is not definite positive. One appealing solution is to replace the Hessian by a scale matrix that can efficiently accelerate the algorithm with a lower computational cost. In particular, many methods have proposed the Fisher information matrix as a preconditioning matrix in the Langevin diffusion [START_REF] Girolami | Riemann manifold Langevin and Hamiltonian Monte Carlo methods[END_REF][START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF] which can be interpreted as the discretization of the MALA algorithm directly on a natural Riemannian manifold where the parameters live. In this work, we propose a new approach where the scale matrix of the Langevin diffusion is chosen according to a Majorize-Minimize strategy.

MAJORIZE-MINIMIZE ADAPTED METROPOLIS-HASTINGS

We focus on the case when the minus-log of the target density function Θ =log π X can be expressed as:

(∀x ∈ R Q ) Θ(x) = Φ(Hx -z) + Ψ(x), (6)
where z ∈ R N , H ∈ R N ×Q , Φ is a continuous coercive differentiable function with an L-Lipschitzian gradient and

Ψ(x) = S s=1 ψ s ( V s x -c s ), (7)
where (∀s ∈ {1, ..., S}) V s ∈ R Ps×Q , c s ∈ R Ps and (ψ s ) 1 s S is a set of positive continuous functions satisfying the following properties:

• (∀s ∈ {1, ..., S}) ψ s is a differentiable function,

• (∀s ∈ {1, ..., S}) ψ s ( √ •) is concave over R + ,
• (∀s ∈ {1, ..., S}) (∃ ωs > 0) such that (∀u > 0) 0 ψs (u) ωs u and lim t→0 ψs (u)/u < ∞. The minimization of (6) using the MM approach consists of performing successive minimizations of its tangent majorant functions [START_REF] Hunter | A tutorial on MM algorithms[END_REF]. Let

x ′ ∈ R Q . A function f is said to be a tangent majorant function of Θ at x ′ provided that f (x ′ , x ′ ) = Θ(x ′ ), f (x, x ′ ) Θ(x) (∀x ∈ R Q ). (8)
Let us assume the existence, for every x ′ ∈ R Q , of a positive definite matrix Q(x ′ ) ∈ R Q×Q such that the following quadratic function, defined for every

x ∈ R Q , f (x, x ′ ) = Θ(x ′ )+(x-x ′ ) ⊤ ∇Θ(x ′ )+ 1 2 (x-x ′ ) ⊤ Q(x ′ )(x-x ′ ) (9
) is a tangent majorant of (6) at x ′ . Then, the MM optimization algorithm reduces to building a sequence (x t ) t∈N through the following scheme:

(∀t ∈ N) x t+1 = x t + ε 2 2 Q -1 (x t )∇ log π X (x t ), (10)
with ε ∈ (0, √ 2]. According to the majorization properties (8), the MM update rule [START_REF] Vacar | Langevin and Hessian with Fisher approximation stochastic sampling for parameter estimation of structured covariance[END_REF] will produce a monotically decreasing sequence (Θ(x t )) t∈N that converges to a local minimum of Θ. We take up this idea to speed up the Langevin diffusion by using the inverse of the curvature matrix Q(x t ) as a scale matrix in (4). Similarly to Newton-based MCMC methods, the drift term, assuming zero curvature changes, proposes, from a current state x t , a state with a higher value of log π X , resulting from an iteration of MM algorithm on log π X . Then, the obtained proposal reduces to a noisy version of a MM iteration for minimizinglog π X . Since the deterministic MM optimization approach can suffer from convergence to a local minimum in the nonconvex case, the addition of the noise component can solve this issue. The resulting 3MH sampler is described by Algorithm 1.

Algorithm 1: Majorize-Minimize adapted Metropolis-Hastings algorithm

0. Initialize x 0 , t = 0, ε ∈ (0, √ 2] 1. Compute A(x t ) = Q -1 (x t ) and g(x t ) = ∇ log π X (x t ) 2. Generate x * ∼ q(x t , •), where q(x t , •) = N x t + ε 2 2 A(x t )g(x t
), ε 2 A(x t ) 3. Accept with probability α(x t , x * ) = min 1, π X (x * )q(x * , x t ) π X (x t )q(x t , x * ) 3. Set t ← t + 1 and go to 1 until convergence.

There remains to define a set of suitable preconditioning matrices {Q(x)} x∈R Q . According to [START_REF] Chouzenoux | A majorize-minimize subspace approach for ℓ 2 -ℓ 0 image regularization[END_REF], convex quadratic tangent majorants of ( 6) can be obtained by using

(∀x ∈ R Q ) Q(x) = µH ⊤ H + V ⊤ diag{ω(x)}V + ζ I Q , ( 11
) where µ ∈ [L, +∞[, V = V ⊤ 1 , . . . , V ⊤ S ⊤ and ω(x) = (ω i (x))
P i=1 is such that, for all s ∈ {1, . . . , S}, p ∈ {1, . . . , P s },

ω P1+P2+...+Ps-1+p (x) = ψs ( V s x -c s ) V s x -c s . (12)
Moreover, ζ 0 is a constant that ensures the invertibility of Q(x) for every x ∈ R Q . In the context of large scale problems, the inversion of the curvature matrix [START_REF] Hunter | A tutorial on MM algorithms[END_REF] at each iteration may become intractable. We thus also propose to resort to the following alternative choice described in [START_REF] Chouzenoux | Variable metric forward-backward algorithm for minimizing the sum of a differentiable function and a convex function[END_REF], which can be understood as a diagonal approximation of ( 11):

(∀x ∈ R Q ) Q(x) = (µ H 2 + ζ)I Q + Diag P ⊤ ω(x) , (13
) where 1 Q is the unit vector of R Q and P ∈ R P ×Q , with P = s P s , is the matrix whose elements are given by

(∀i ∈ {1, . . . , P })(∀j ∈ {1, . . . , Q})P i,j = |V i,j | Q k=1 |V i,k |.

EXPERIMENTAL RESULTS

Objective The experiments we carried out deal with the recovery of a multicomponent image with B components degraded by a blur modelled by a linear operator D and an additive Gaussian noise w with covariance matrix Λ. The restoration is performed in the wavelet transform domain: our objective is to compute the MMSE of x ∈ R Q through the linear model defined in [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] with H = DF * where F * denotes a linear synthesis wavelet operator. Note that the wavelets coefficients are grouped into M subbands of size Q m , m ∈ {1, . . . , M } and, for each subband m, we can extract the set of vectors (x m,q ) Qm q=1 containing the wavelet coefficients located at the same spatial position q through all the B channels using a B × Q permutation matrix P m,q such as x m,q = P m,q x. Prior distribution Similarly to [START_REF] Marnissi | Generalized multivariate exponential power prior for wavelet-based multichannel image restoration[END_REF], we assume that the vectors (x m,q ) Qm q=1 are realizations of a random vector with a Generalized Multivariate Exponential Power (GMEP) distribution whose multivariate probability density function p GMEP is defined, for every u in R B , by

p GMEP (u; θ m ) = C m |Σ m | -1/2 g u ⊤ Σ -1 m u; β m , δ m , (14)
where

θ m = {β m > 0, δ m > 0, Σ m }, for every t ∈ R + , g(t; β m , δ m ) = exp -1 2 (t + δ m ) βm , and C m is the associ- ated normalization constant. Σ m is related to the covariance matrix Γ m through Σ m = K 2 βm,δm Γ m with K 2 βm,δm = B ∞ 0 t B 2 -1 e -1 2 (t+δm) βm dt ∞ 0 t B 2 e -1 2 (t+δm) βm dt . ( 15
)
The GMEP is an elliptical distribution that reflects the sparsity of the coefficients and, following [START_REF] Sanchez Manzano | Multivariate exponential power distributions as mixtures of normal distributions with Bayesian applications[END_REF] it can be proved that it is a scale mixture of Gaussian distributions when β m < 1. Its mixing density is expressed as follows:

h βm,δm (ν) = 2 B 2 Γ( B 2 )ν B-3 S βm,δm ( 1 2 ν -2 , 2 -1 βm ) ∞ 0 t B 2 -1 e -1 2 (t+δm) βm dt (16)
where, for every α 1 ∈ (0, 1), α 2 > 0, σ > 0, S α1,α2 (•, σ) = e -1 2 α2 S α1 (•, σ), and S α1 (•, σ) is the alpha-stable density whose characteristic function is exp(-σ α1 |•| α1 e -i π 2 α1 sign(•) ) [START_REF] Sanchez Manzano | Sequences of elliptical distributions and mixtures of normal distributions[END_REF].

Proposed priors for the hyperparameters In the following, we suppose that (δ m ) 1 m M is known. We also assume that, for every m ∈ {1, . . . , M }, parameters β m and Γ m are independent and we denote by p βm and p Γm their respective prior density functions. Since wavelet coefficient images have leptokurtic histograms [START_REF] Simoncelli | Bayesian denoising of visual images in the wavelet domain[END_REF], we set p βm = U (0, 1). We use an inverse Wishart prior for Γ m , with parameters fixed according to a prior knowledge about Γ m .

Posterior distributions

The posterior distributions of the GMEP hyperparameters have complicated form and there is no practical way for designing algorithms to simulate samples from them especially for the scale matrix. Since for every m ∈ {1, . . . , M }, β m ∈ (0, 1), we propose to exploit the fact that GMEP is a scale mixture of normal distributions. Then, there exists a vector v m = (v m,q ) Qm q=1 of random variables v m,q such that K βm,δm v m,q ∼ h βm,δm and, for all q ∈ {1, . . . , Q m }, x m,q is drawn independently from a zero mean Gaussian distribution with covariance matrix v 2 m,q Γ m . Hence, the posterior distributions of Γ m reduces to an inverse Wishart distribution and the posterior distribution of β m is related to a product of densities of stable distributions [START_REF] Sanchez Manzano | Multivariate exponential power distributions as mixtures of normal distributions with Bayesian applications[END_REF].

Let Ω = H ⊤ Λ -1 H + M m=1 Qm q=1 v -2
m,q P ⊤ m,q Γ -1 m P m,q then the posterior distribution of x reduces to a normal distribution with mean µ = Ω -1 H ⊤ Λ -1 z and covariance matrix Ω -1 . Note that sampling from high-dimensional Gaussian distributions is often very difficult since matrix factorization (Cholesky, QR, square root) is not always possible because of its high computation cost and/or memory requirements. Some solutions have been proposed for some special structures of the covariance matrix (circular, sparse,. . . ) [START_REF] Chellappa | Classification of textures using Gaussian Markov random fields[END_REF][START_REF] Orieux | Sampling highdimensional Gaussian distributions for general linear inverse problems[END_REF]. In this work, we propose to use a step of the 3MH algorithm with the diagonal curvature matrix. The interest of this alternative solution is that it does not require any assumption on the structure of the covariance matrix. Note that we follow the method proposed in [START_REF] Sanchez Manzano | Multivariate exponential power distributions as mixtures of normal distributions with Bayesian applications[END_REF] for the sampling of v m,q . Moreover, MH steps are used to generate samples from the posterior distribution β m .

Results

The test image is a remote sensing multispectral SPOT image of size 128 × 128 with three components (B = 3) and, corresponding to a scene depicting the city of Tunis. This image is artificially blurred with a cosine blur FTM and corrupted with a zero-mean white Gaussian noise with a variance adjusted so as to correspond to an initial averaged Blurred Signal-to-Noise Ratio (BSNR) equal to 21.64 dB. We apply a 2-resolution wavelet orthonormal decomposition us-ing a Symmlet wavelet transform of order 8. We run the hybrid Gibbs sampler for 8,000 iterations, reject the 6,000 first ones as a burn-in, and take the last 2,000 results as samples for the target data. Fig. 1 provides the evolution of improved SNR (ISNR) with respect to the computational time using different algorithms to generate π X , when performing tests on an Intel Core i7 CPU, @ 3.00 GHz and using a Matlab 7.12 implementation. It can be observed that the 3MH algorithm reaches stability faster than MALA. In fact, MALA algorithm requires less time per iteration but our algorithm converges in a significantly smaller number of iterations. The obtained samples of the wavelet coefficients are then used to compute the empirical MMSE estimator for the original image. We have used the hyperparameters estimated by the Gibbs sampler to run the method described in [START_REF] Marnissi | Generalized multivariate exponential power prior for wavelet-based multichannel image restoration[END_REF] which computes the MAP estimate with a GMEP prior using the MM Memory Gradient algorithm. Table 1 reports the results obtained for the different components in terms of SNR, BSNR and Structural Similarity Index (SSIM). It can be observed that the MMSE estimator shows better performance than the MAP. This can also be observed on Fig. 2 showing the visual improvement for the first component of the image.

CONCLUSION

In this paper, we have proposed a new MCMC algorithm that can be considered as a scaled MALA where the scale matrix is adapted at each iteration with a MM strategy. We have applied this algorithm to compute the MMSE estimator of a multicomponent image from its blurred version. Experimental results indicate the good performance of this new MCMC method. Note that the proposed approach can be applied to general models where the posterior distribution is non-Gaussian.

Fig. 1 .

 1 Fig. 1. Convergence speed of 3MH, MALA and MH.

Fig. 2 .

 2 Fig. 2. (a) Original image, (b) degraded version of the first component (SNR = 21.72 dB, SSIM= 0.729), (c) restored version using the MAP estimator (SNR =24.95 dB, SSIM = 0.860), (d) restored version using the MMSE estimator (SNR =25.20 dB, SSIM = 0.872).

Table 1 .

 1 Restoration results b = 1 b = 2 b = 3 Average

	Initial	BSNR 22.32 20.29 22.30 21.64 SNR 21.72 19.56 21.96 21.08 SSIM 0.729 0.761 0.720 0.737
	MAP	BSNR 26.76 24.16 26.48 25.80 SNR 24.95 22.33 25.19 24.16 SSIM 0.860 0.863 0.843 0.855
	MMSE	BSNR 27.34 24.75 27.06 26.38 SNR 25.20 22.59 25.51 24.43 SSIM 0.872 0.874 0.855 0.867
	1.5
	1
	0.5
	0 0	100	200	300	400	500