Lei Yu

Pierre Maurel

Christian Barillot

Remi Gribonval

Quantifying CBF from Aterial Spin Labeling via Diverse-TI: sampling diversity or repetitions ?

Keywords: Arterial Spin Labeling, Cerebral Blood Flow, Compressive Matched Filter, Diverse-TI, Diversity, Repetition

Arterial Spin Labeling (ASL) is a noninvasive perfusion technique which allows the absolute quantification of Cerebral Blood Flow (CBF). The perfusion is obtained from the difference between images with and without magnetic spin labeling of the arterial blood and the captured signal is around 0.5-2% of the magnitude of the labeling images, so the noise is one of the main problems for further data analysis. Classical method, Mono-TI, for CBF quantification is averaging repetitions with only one Inversion Time (TI) -the time delay between labeling and acquisition to allow the labeled blood to arrive the imaging slice. It improves the robustness to noise, however, cannot compensate the variety of Arterial Arrival Time (AAT). In this paper, Diverse-TI is proposed to exploit different TI sampling instants (sampling diversity) to improve the robustness to variety of AAT and simultaneously average repetitions with each TI (sampling repetitions) to improve the robustness to noise. Generally, the sampling diversity is relatively small and can be considered as compressed measurements, thus the Compressive Matched Filter (CMF) enlightened from sparsity is exploited to directly reconstruct CBF and AAT directly from compressed measurements. Meanwhile, regarding the CBF quantification performance, the compromise between the sampling repetition and sampling diversity is discussed and the empirical protocol to determine the sampling diversity is proposed. Simulations are carried out to highlight our discussions.

Quantification du flux sanguin cérébral en ASL par Diverse-TI : diversité d'échantillonnage ou répétitions ?

Résumé : L'imagerie par marquage de spins artériels est une modalité d'imagerie par perfusion non-invasive qui permet la quantification absolue du flux sanguin cérébral (CBF). La perfusion est obtenue par soustraction de deux images : avant et sans marquage des spins artériels. Le signal capturé est alors d'environ 0.5-2% de l'amplitude du signal de chaque image. Le bruit est donc un des principaux problèmes. La méthode classique de quantification, Mono-TI, effectue la moyenne de plusieurs répétitions à un certain temps d'inversion (TI) correspondant au délai entre le marquage et l'acquisition. Cette moyenne améliore le rapport signal-sur-bruit, mais choisir un seul TI ne permet pas de prendre en compte les éventuelles variations du temps d'arrivée du bolus (AAT). Dans cet article, nous proposons une nouvelle technique, Diverse-TI , qui échantillonne à différents TI. Il s'agit alors de trouver un compromis entre la diversité d'échantillonnage et le nombre de répétitions par TI. En général la diversité d'échantillonnage reste relativement faible et on peut donc se placer dans le cadre d'acquisitions compressés (en anglais compressed sensing). La méthode Compressive Matched Filter (CMF) est donc utilisée pour estimer le CBF et l'AAT directement à partir des mesures. Concernant la quantification du CBF, le compromis entre la diversité d'échantillonnage et le nombre de répétitions par TI est discuté et un protocole expérimental est proposé.

Mots-clés : Marquage de spins artériels, flux sanguin cérébral, Compressive Matched Filter, diversité, répétition

Introduction

The Arterial Spin Labeling (ASL) is an MRI (Magnetic Resonance Imaging)-based perfusion technique which uses the magnetically tagged water as a freely diffusible tracer to measure perfusion non-invasive. This blood water is first labeled with a radio-frequency pulse in the neck of the patient. After a delay, called Inversion Time (TI), which allows the labeled blood to arrive in the brain, a labeled image of the brain is acquired. A control image is also acquired without labeling and the CBF estimation is done on the difference between the control and labeled image. For a fixed TI t, the perfusion signal is generally described by a kinetic model introduced by Buxton et. al. in [START_REF] Buxton | A general kinetic model for quantitative perfusion imaging with arterial spin labeling[END_REF]:

M ∆ (t) =      0, t ∈ [0, ∆) 2αf M 0b (t -∆)q p (t) exp(-t/T 1b ), t ∈ [∆, ∆ + τ ) 2αf M 0b τ q p (t) exp(-t/T 1b ), t ∈ [∆ + τ, ∞) (1)
where f is the CBF, α measures the labeling efficiency, M 0b is the equilibrium magnetization of arterial blood, ∆ is the Arterial Arrival Time (AAT) to the interesting slice, τ is the temporal width of the bolus, T 1b is the relaxation time in blood and q p (t) is considered to be approximately equal to 1. Specifically, in QUIPSS II [START_REF] W M Luh | QUIPSS II with thin-slice TI1 periodic saturation: a method for improving accuracy of quantitative perfusion imaging using pulsed arterial spin labeling[END_REF], a saturation pulse is given to the tagged bolus after a fixed period τ of the tagging time, thus the temporal width of the tagged bolus is known a priori.

The principal task of ASL analysis is quantifying the absolute value of CBF f . The traditional technique, which is called Mono-TI, only uses a single TI t after which the acquisition is gathered, where t is assumed to be bigger than ∆ + τ to guarantee a delay long enough to let the magnetic tagged blood arrive the interesting imaging slice. Then the quantification of CBF is a direct ratio between the magnetization difference M ∆ and the known terms, when t > ∆ + τ .

However, two major problems prevent the efficiency of Mono-TI. At first, as the amplitude of the difference ASL signal is usually around 0.5-2% of the control image magnitude, its SNR is thus not sufficient for further analysis. Commonly, numbers of sampling repetitions are acquired for single TI (typically more than 30) and then averaged to improve the SNR. Secondly, the assumption t > ∆ + τ is not easy to be guaranteed in Mono-TI technique. Since ∆ is varying between different patients, ages and physical situations. Even excluding these factors, the AAT is still an uncertain value due to the presence of laminar and turbulence flow, complicated vessel networks and cardiac pulsations [START_REF] Wu | The effects of flow dispersion and cardiac pulsation in arterial spin labeling[END_REF]. This fact leads that the magnetic labeled blood may not have reached the imaging slice if the sampling time is too small and thus leads to underestimation of CBF. One possible way is giving the TI large enough, however, the magnetization difference M ∆ might be too small if TI is too big and thus leads to very low SNR.

In this paper, another measuring procedure with different sampling times is investigated for the CBF quantification, where the collected ASL data at different sampling times are captured during separate ASL RF pulse periods. It is different from the classical Multi-TI [START_REF] Reinoud | Arterial spin labeling perfusion MRI at multiple delay times: a correlative study with H(2)(15)O positron emission tomography in patients with symptomatic carotid artery occlusion[END_REF][START_REF] Macintosh | Multiple inflow pulsed arterial spin-labeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF][START_REF] Wang | Pediatric perfusion imaging using pulsed arterial spin labeling[END_REF] and thus called Diverse-TI to avoid ambiguities. Multi-TI is a more recent and confidential sequence than Mono-TI and is not currently available on most MR-Scanners, while Diverse-TI proposed in this paper can easily be produced since it only uses regular ASL sequences. Before further analysis on real ASL data, the first and also the essential step is giving a protocol to concretely design the measuring procedures. It is clear that some issues should be preferentially considered:

(1) sampling repetition (for the same sampling time), which is intended to improve SNR; (2) sampling diversity (sampling at different times) which is intended to compensate for the inexact knowledge of parameters such as ∆, etc. For practical considerations, the most crucial criteria in real clinical studies is the total measuring time which is generally limited in a reasonable period.

Thus the values of sampling repetition and sampling diversity can not be designed as large as possible to improve CBF quantification performance. A method to explicitly guide the design of these parameters is applaudable. However, to the best knowledge of the authors, the existing papers rarely focus on this point.

Consequently, as a preliminary study, the contribution of this paper is twofold. First, it formulates the considered problem as an instance of Compressive Matched Filter (CMF) [START_REF] Mark A Davenport | Signal Processing With Compressive Measurements[END_REF]. Then, it numerically investigates the trade-offs between repetition and sampling diversity (for the same total measurement time) on CBF estimation with CMF. In a word, the main question investigated in this paper is: given a total measurement time, shall we favor repetitions or sampling diversity?

2 Diverse-TI Technique Since parameters including τ, T 1b , α, M 0b and q p (t) are (or assumed to be) known a priori, the signal model (1) can be simplified as M ∆ (t) = f × g(∆; t) where the time related term g(∆; t) will be called "wave form" and can be defined by (1) with f = 1.

-Sampling Repetition

In practice, an important level of noise is affecting ASL measurements, hence the signal captured at time t can be expressed as follows:

y(t) = M ∆ (t) + (t) = f × g(∆; t) + (t)
where (t) is the noise term and is assumed to be Gaussian in this paper. In our setting, the measure at time t is repeated R times (sampling repetition), and averaging all measures divides the noise variance by R and thus it will improve the SNR according to the following equation:

Q R = 10 log 10 R + Q 1 (2)
where Q 1 is the SNR of one repetition and Q R is the SNR after R repetitions.

-Sampling Diversity Define g(∆) = {g(∆, iT s )} N i=1 where T s is a reference regular sampling interval and t 1 ...t M (possibly irregular) time instants located on this regular grid, where M N . With different TI (sampling diversity) t 1 ...t M , the ASL data are respectively collected y {y(t 1 ), ..., y(t M )} T ∈ R M . Then we can write the sampling model in the form of linear operation:

y = Φ(f × g(∆) + )
where Φ ∈ R M ×N is a sensing matrix, which verifies ΦΦ T = I M .

At last, the task is to find the CBF parameter f and the optimal ∆ that best match the observations y. To fulfill it, Compressive Matched Filter (CMF) [START_REF] Mark A Davenport | Signal Processing With Compressive Measurements[END_REF] is exploited in the following sections.

Compressive Matched Filter

In the case considered in this paper, the CBF quantification problem can be formed as a multiple detection problem to distinguish the following hypothesis:

H i : y = Φ(f g(∆ i ) + ), with i ∈ {1, ...

, d}

Inria

where f is the CBF value, ∆ 1 , ..., ∆ d are d possible AAT values that should be assumed before and ∼ N (0, σ 2 I N ) is the white noise term. In this case, the evidence under different hypothesis, for i ∈ {1, ...d}, can be written

p(y|H i ) = exp -1 2 (y -Φf g(∆ i )) T (σ 2 ΦΦ T ) -1 (y -Φf g(∆ i )) (2π) N/2 |σ 2 ΦΦ T | 1/2
Then the final detection is carried out by finding the hypothesis with the biggest conditional probability:

(f s , H s ) = arg max f,Hi,i∈{1,...,d} p(y|H i )
Thus by differentiating w.r.t. f , we can obtain its estimation

fi = y T (ΦΦ T ) -1 Φg(∆ i ) g(∆ i ) T Φ T (ΦΦ T ) -1 Φg(∆ i ) (3)
and by taking logarithm we obtain an equivalent test that simplifies to

H s = arg max Hi,i∈{1,...,d} y - 1 2 Φ fi g(∆ i ) T (ΦΦ T ) -1 Φ fi g(∆ i ) = arg max Hi,i∈{1,...,d} 1 2 y
T (ΦΦ T ) -1 Φg(∆ i ) 2 g(∆ i ) T Φ T (ΦΦ T ) -1 Φg(∆ i ) ΦΦ T =I

====== arg max

Hi,i∈{1,...,d} y,

Φg(∆ i ) Φg(∆ i ) 2 (4)
Above all, using (4) and (3), we can respectively obtain the estimation of arrival time ∆ = ∆ s and the estimation of CBF f = f s .

Mono-TI technique v.s. Diverse-TI technique

In fact, M = 1 is exactly the Mono-TI technique. In order to well illustrate the performance of Mono-TI technique, a concrete simulation is carried on. In the following, the parameters of ASL model (1) are arbitrarily set as: M 0b = 100, α = 0.9, T 1b = 1.2s and τ = 1s, and Fig. 1(left) shows one example when AAT ∆ = 1 and CBF f = 10 mL/(100g)/min. Obviously, given constant noise variance σ, the highest SNR of single repetition Q 1 is reached when t i = ∆ + τ while it will decrease as t i increasing.

The assumption on TI t i > ∆ + τ is supposed to be verified. At first, t i can not be much smaller than ∆ + τ otherwise the tagged blood can not reached the interesting slice on time, as shown by the red diamond point in Fig. 1(left) when t i = 1s; meanwhile, it can not be too big otherwise the captured signal will be too small which makes the SNR Q 1 very low, as shown by the green circle point in Fig. 1(left) when t i = 4.5s. Then the CBF f can be calculated directly by f = y(t i )/g(∆; t i ), however, the AAT ∆ can not be estimated in Mono-TI technique.

Then Mono-TI technique is carried out in the following procedure. Fix the CBF f = 10 mL/(100g)/min, give one AAT value ∆, sample at the time t i with noise corruption, and then one single repetition of the ASL signal is simulated. Without loss of generality, the time used for one single repetition is assumed T S = 5s (this value is set identically for all simulations). Then as an example, we can set the total time T = 30min, so according to (5) the Mono-TI technique allows R = 360 repetitions. According to (2), the SNR can reach to Q R = 30.56dB with Q 1 = 5dB. However, this value is very sensible to the sampling location, or equivalently, the AAT value. This fact can be well illustrated by Fig. 1, where the result of Mono-TI is obtained by fixing t i = 1s, 3s or 4.5s, then ranging AAT ∆ from 0.1s to 3s (the performance shown in Fig. 1(right) is measured through the quantification SNR averaged over 1000 restarts). In Fig. 1, the result when t i = 1s shows the situation when TI is too short, t i = 4.5s shows the situation when TI is too long, even with proper TI t i = 3s, the performance of quantification might largely decrease when AAT and the TI are mismatch. In the framework of Diverse-TI , following the CMF algorithm, we first assess a large interval where the value of AAT should locate. In this simulation, this interval is set to ∆ ∈ [0.1, 3]. We choose the number1 of sampling diversity to M = 10, and the sampling repetitions are set identically according to (5) with the same total measuring time as Mono-TI, then the CBF quantification performances through CMF algorithm with different AAT value are obtained, as shown in Fig. 1. Diverse-TI can give constant performance which is also better than Mono-TI with a majority of AAT.

As a conclusion, in the case of M = 1, the CBF quantification quality might be better than M > 1, however, its performance is much sensitive to the sampling time, which is not easy to be designed practically. Moreover, as a supplement, Diverse-TI can also estimate the AAT.

Experiments and simulations

Without loss of generality, the time used for one single sampling repetition can be assumed to be T S , including the process of one time control imaging and tagging imaging. Then the total measuring time can be expressed as:

T ∝ M • R • T S ( 5
)
where M indicates the level of sampling diversity i.e. the number of TI and R the number of sampling repetitions. For practical considerations, the most crucial criteria in real clinical studies is the total measuring time which is generally limited in a reasonable period. In one aspect, the sampling repetitions improve the signal SNR; while in the other aspect, the sampling diversity promote the robustness to the variety of AAT. Consequently, the compromise between the sampling diversity and sampling repetitions needs to be determined. In other words, given the maximum total time T , how to optimally choose R and M to reach better CBF quantification?

Acquisition

A healthy volunteer was scanned on QUIPSS II with image size 64 × 64 and the volume size 3 × 3 × 4mm. Fourteen slices were acquired with TIs for the lowest slice T I 1 /∆T I/T I N = 1200ms/100ms/2200ms and an extra delay of 45ms between slices, i.e. with diversity M = 11. Meanwhile, 40 repetitions (control-labeled pairs) for each TI were captured, i.e. with repetitions R = 40. All 40 repetitions were co-registered to compensate for the patient motion with the first one using 3D rigid-body transformation implemented in SPM-Toolbox.

To estimate the CBF and AAT, the parameters in the perfusion model to the sequence of Diverse-TI difference images are set according to the most used in the literatures: tissue relaxation time T 1 = 900ms, blood relaxation time T b = 1500ms, blood-tissue ratio coefficient λ = 0.9 and the magnetization efficiency α = 0.9. The time width of the labeled blood in QUIPSS II is fixed as τ = 700ms.

Performance Assessment

The performance assessment to simulations can be easily calculated by referring the ground truth of CBF or AAT, and then the quantification SNR for CBF and AAT can be defined:

Q f (M, R) = 20 log 10 f true f (M, R) -f true and Q ∆ (M, R) = 20 log 10 ∆ true ∆(M, R) -∆ true
with f (M, R) the estimation of CBF when diversity is M and repetition is R, and f true , ∆ true are respectively the groundtruth of CBF and AAT. It is worth mentioning that, instead of voxelwised SNR, the quantification SNRs (both for CBF and for AAT) are computed over the whole perfusion images.

However, the ground truth is not available for real data, and thus an approximation to the ground truth is computed from all available data (M = 11 and R = 40) using CMF. Then the performance are assessed regarding to the approximated ground truth, noted by f all = f (11, 40) and ∆ all = ∆(11, 40). At the end, the performance assessment for CBF quantification with the real data is carried out by calculating the following reconstruction SNR:

Q f (M, R) = 20 log 10 f all f (M, R) -f all
Accordingly, for AAT:

Q ∆ (M, R) = 20 log 10 ∆ all ∆(M, R) -∆ all
Yu & others Remark 5.1 Although the approximated ground truth can not perfectly represent the real ground truth, we can still assume that the expected quantification with more data will be closer to the real ground truth than that estimated from less data.

Sampling diversity and repetitions with random sampling locations

In section 4, it has been shown that M = 1 is worse than M > 1, while in this subsection, we will propose a method to answer the question that how many sampling time locations are enough? Similarly, the total time T is the crucial parameter in clinical studies. Consequently, given a fixed total time T , the relationship between the CBF quantification performance and the number of sampling locations is essential to determine the required M .

In this simulation, the CBF f is fixed equal to 10mL/(100g)/min and the AAT ∆ is randomly choosing from the interval [0.1, 2], then the ASL signal are generated and corrupted by noises. Then to simulate the sampling procedure, the total time T is fixed respectively to 5, 10, 30, 60min and the number of sampling time locations M is varying from 2 to 30, thus the number of repetitions R is determined via (5), then samples are collected. Since ∆ ∈ [0.1, 2] is assumed, the CMF can be used to carry out the CBF quantification and AAT estimation, and the SNR of CBF quantification and AAT estimation are respectively shown in Fig. 2 (b) and(c). Meanwhile, the relationship between sampling repetition and sampling diversity can also be drawn according to [START_REF] Macintosh | Multiple inflow pulsed arterial spin-labeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF], as shown in Fig. 2 (a). From the results, as expected, we can find that as M increases, the performance of Diverse-TI technique is improving before reaching an asymptotic value. From this simulation, we can find that the convergent point is around M = 8, as shown by the performance of CBF (or AAT) quantification in bottom figures of Fig. 2. Then according to [START_REF] Macintosh | Multiple inflow pulsed arterial spin-labeling reveals delays in the arterial arrival time in minor stroke and transient ischemic attack[END_REF], the repetitions at one sample time R can be calculated. It is worth mentioning that this simulation is only an example of choosing the number of multiple sampling time locations, for concrete applications, the curve in Fig. 2 should be re-drawn with appropriate parameters, and then the best value should be re-selected.

Best Sampling Locations

From the perfusion model, the perfusion signal magnitude has a peak at some time position. Assuming that the measuring model is Gaussian, one can find that the measuring SNR will be Inria the highest at the peak position and thus give the best reconstruction performance for CBF. On the other hand, the estimation of arrival time ∆ is very sensitive to the first sampling point, for instance y(1) = 0 and y(1) = 0 will give very different results for ∆. Consequently, to eliminate the influence of the sampling locations, we first determine the best sampling locations, where the samples will give the best performance at each diversity and repetition level.

In this experiment, the estimations are obtained from all possible combinations of sampling positions at every sampling diversity, for instance, when M = 2, it will have 11 2 possible sampling patterns. Then we can choose the best sampling locations as the sampling patterns with the best estimation performance, either for CBF or for AAT. The results are listed in Tab. 1. To make a compromise between the best sampling locations for CBF and for AAT, we consider the following value as the criterion:

Q f +∆ = Q f + Q ∆
Then choose the biggest Q f +∆ to determine the best sampling locations both for CBF and AAT, as shown in the last column in Tab. 1.

Remark 5.2 Considering the results listed in Tab. 1, we can get another table by counting the number of time locations selected as the best sampling locations, as shown in Tab. 3. From Tab. 3, we can conclude that in order to obtain better performance, the optimal way to choose the sampling time positions are the following: first capture a sample at the starting point, then some samples in the middle and one sample at the end.

In the following experiments, we will fix the sampling locations following the best sampling locations shown in the last column of Tab. 1 for each diversity level, then investigate the effects of different repetitions.

RR n°8258

Performance Evaluation Map with Best Sampling Locations

In section 5.3, we have discussed that the level of diversity can be determined according to the performance evolution curves, as shown in Fig. 5. And from the result, we can easily assess the optimal diversity and optimal repetitions for a given total measuring time. In this section, we will turn to use the performance evolution map to determine an optimal compromise, where we can find the optimal diversity or repetitions much easier.

Experiments on simulated data

In order to make the simulated data much closer to the real case, we first generate the AAT values from the estimation of AAT from real perfusion data, as shown in Fig. 3. Then the simulated perfusion signals are generated with the generated AAT and a fixed CBF value (it is set to 80 here). Thus for each level of diversity, the best TIs are chosen from Tab. 1, which means part of the data has been used to estimate the CBF and AAT. Meanwhile, different levels of noises are added with respect to different level of repetitions, and the values of noise SNR can be referred by [START_REF] W M Luh | QUIPSS II with thin-slice TI1 periodic saturation: a method for improving accuracy of quantitative perfusion imaging using pulsed arterial spin labeling[END_REF]. Then CMF is exploited to do the CBF and AAT quantification, for every level of diversity and every level of repetition. And the quantification SNRs of each level of diversity and repetitions for CBF Q f (M, R) and for AAT Q ∆ (M, R) can then be calculated. Finally, we re-run the simulations 100 times with different draw of noises, and correspondingly, 100 quantification SNRs for CBF and for AAT are obtained, respectively denoted by Q

(1)

f (M, R), ..., Q (100) f (M, R)
and

Q (1) ∆ (M, R), ..., Q (100) ∆
(M, R). Then the mean value of the quantification SNR for CBF and for AAT can be calculated as following:

Qf (M, R) = 1 100 100 i=1 Q (i) f (M, R) Q∆ (M, R) = 1 100 100 i=1 Q (i) ∆ (M, R) (6)
And the map of Qf (M, R) and Q∆ (M, R) are respectively shown in Fig. 4 (a) and (c). Meanwhile, the standard variation of the quantification SNR for CBF and for AAT can also be calculated as following:

Qf (M, R) = 1 100 100 i=1 (Q (i) f (M, R) -Qf (M, R)) 2 Q∆ (M, R) = 1 100 100 i=1 (Q (i) ∆ (M, R) -Q∆ (M, R)) 2 (7)
And the map of Qf (M, R) and Q∆ (M, R) are respectively shown in Fig. 4 (b) and (d). After that, we can then plot the isoheight of the total measuring time for different diversity and repetitions. As shown in Fig. 4, three isoheight curves are given respectively to T ∝ 100, 200, 300.

First, we can focus on the mean value of the quantification SNR for CBF and AAT, and we can find that the results are similar to the results presented in section 5.3. The conclusion of both results is that a diversity M around 5 can give similar results as when the diversity M much bigger. One more conclusion from Fig. 4 is that the performance will drop down when the sampling repetitions are decreasing.

Secondly, let us focus on the standard variation of the quantification SNR for CBF and AAT. For CBF, as shown in Fig. 4 (b), the variation is much lower when the diversity M ≤ 5. While for AAT, as shown in Fig. 4 (d), the sharp transition of variation appears when the diversity M is from 5 to 4. It is shown that the lower diversity gives the more robust quantification both for CBF and AAT.

Above all, we can conclude from the simulated results that the best diversity level is M = 5, and then the corresponding repetition level can be chosen according to the total measuring time T .

Remark 5.3 One abnormal phenomenon is that lower diversity gives more robust quantification

RR n°8258

for fixed repetitions. It is due to the reason that even though higher diversity indeed improves partial of the voxels, there still exist lots of voxels of which the performance cannot be improved by only increasing the diversity. On the other hand, the performances with low diversity for all voxels are low. Consequently, increasing the diversity might increase the variation of the quantification SNR.

Another abnormal phenomenon is that the variation of the SNR for AAT is increasing as increasing the repetitions for fixed diversity. It is due to the fact that for fixed diversity, by increasing repetitions can improve the performance for some of the voxels, while others do not change to much. Thus the average SNR is increasing while the standard variation is also increasing. In this subsection, we will use the similar procedure to obtain the performance evaluation map for real data. It is worth claiming that the references of CBF and AAT of real data are their quantifications by using all data, i.e. the quantification when M = 11 and R = 40.

Experiments on real data

Then the one-running process i is described as following: Firstly, the sampling locations for each diversity level are determined by the best TIs for both CBF and AAT, as shown in Tab. 1. Denote Ω the set of index of sampling repetitions, thus Ω = {1, ..., 40}. Then for repetition level R, there will be P (R) = |Ω| R different selections and the selected samples are randomly chosen from the set Ω. After that, the quantification of CBF and AAT can be calculated by CMF for different diversity levels and repetition levels, i.e.

Q (i) f (M, R) and Q (i)
∆ (M, R). Finally, we re-do the one-running process 100 times with randomly selection of repetitions, i.e. i is ranging from 1 to 100. The performance evaluation map can be obtained by averaging the quantification SNR according to [START_REF] Wang | Pediatric perfusion imaging using pulsed arterial spin labeling[END_REF] and calculating its standard variation according to [START_REF] Mark A Davenport | Signal Processing With Compressive Measurements[END_REF], as shown in Fig. 5. In addition, the CBF and AAT quantifications from real data with different diversity and repetition levels are given in Fig. 6.

The first impression of the performance evolution maps for real data (shown in Fig. 5) and the simulated data (shown in Fig. 4) is that there exist lots of differences:

• Fixing the repetition level, for real data, increasing the diversity will always improve the average SNR for both CBF and AAT; while for simulated data, the average SNR will increase until the diversity reaches M = 5 then stays almost constant as increasing the diversity.

• Considering the standard variation of performance SNR for AAT, for real data, increasing repetition level will decrease the variation; while for simulated data, the variation is increasing (in general) as increasing the repetition level.

For the first difference, it is due to the fact that the Buxton's model used here is not proper for the real data, as shown in Fig. 7. We first generate the perfusion signal with Harbe's model with Gaussian distribution, then the Buxton's model is exploited for CBF and AAT quantification by CMF. The resulted performance evaluation map (Fig. 7) is very similar to those of the real data (Fig. 5). In other words, we can also conclude from the results that Buxton's model is not precise for the real perfusion data. The reason of the second difference is because of the following. The value |Ω| R is varying as increasing the repetitions and if the number of re-running times are bigger than P (R), there will exist some identical selections from Ω and thus it will influent the standard variation map. Particularly, the variation when R = 40 will equal to 0.

Unfortunately, the results from the real data can not clearly give any clues on determining the compromise between diversity level and repetition level under fixed total measuring time.

Discussion and Conclusion

In this report, we first introduced the Diverse-TI, which exploits samples at multiple TIs to quantify the CBF and AAT. Many repetitions are required for each TIs sampling point to improve the data fidelity. Then the question is how many repetitions and how many TIs should we determine to get as highest performance as possible? In this report, we proposed to use the performance evaluation map to determine the compromise between the diversity level and the repetition level.

The simulated results can clearly give us the best way for design the diversity level and the repetition level. The best diversity level is M = 5, and then the corresponding repetition level can be chosen according to the total measuring time T .

The performance evaluation map from the real data can show us that increasing either diversity or repetitions can improve the performance. However, it cannot clearly give us any clues on the compromise of the diversity level and the repetition level.

Even though the performance evaluation map can give us some clues on designing the sampling procedure for Diverse-TI, the relationships between the ultimate results and the following aspects are still unclear and need further investigations: quantification algorithm, approximated references and perfusion model.

The relation with the exploited quantification algorithm

In this report, we did not analyze the performance of quantification algorithm, i.e. CMF, which is assumed to be good enough here to fulfill the quantification. In fact, many algorithms can be exploited here. Consequently, corresponding to different quantification algorithms, the resulted performance evaluation maps might be different, but it will not completely change the evaluation map. Thus the variation with respect to different algorithms requires further investigations and we believe that using the similar procedure, one can get it's performance evaluation map, with which the compromise between diversity and repetitions can be optimally determined.

The relation with the approximated references

In fact, the reference (the ground-truth) of the CBF and AAT in real world is unknown. Consequently, instead of the ground-truth, the approximated reference estimated from the full data set is exploited in the experiments. However, the fidelity of the approximated reference is not investigated in this report. In fact, the approximated reference exploited here can be assumed to be close to the ground-truth. Even though the error introduced by the approximated reference inherently exists, comparing to the quantification error arisen by reducing the sampling diversity or repetitions, this error is relatively small and thus the corresponding performance evaluation map will not change too much. Fig. 8 presents the performance evaluation map using the approximated reference for simulated perfusion data. Comparing to Fig. 4, we can find that the results obtained with the approximated reference are similar to the results with the correct ground-truth and the similar conclusion on the comprise between diversity and repetition can be made. Consequently, in some sense, it shows that the approximated reference does not influent too much on the final result.

The relation with the perfusion model

The last issue is about the perfusion signal model. As shown by the mismatch of the results between simulation and reality, we found that the misuse of the perfusion signal model will lead very different performance evaluation maps. Consequently, it is better to use more precise perfusion signal model for the real data.

As a conclusion, we proposed a simple multi-TIs sampling procedure for ASL, Diverse-TI, which is superior to the Mono-TI method. Meanwhile, the performance evaluation map is proposed to determine the diversity level and repetition level in the Diverse-TI framework. Although we can not easily find the optimal compromise between diversity and repetition by using the performance evaluation map from the real data, the simulated results give us some clues on designing the two essential parameters.

Figure 1 :

 1 Figure 1: (Top) Kinetic Model of ASL with ∆ = 1 (τ = 1s) and f = 10 mL/(100g)/min (dark solid curve) and its noisy version with average SNR 5dB (grey dashed curve). (Bottom) Comparison between Mono-TI and Diverse-TI (M = 10) methods for different AAT with total time T = 30.

Figure 2 :

 2 Figure 2: Performance Evolution with respect to different sampling diversities M via CMF using different total time.

Figure 3 :Figure 4 :

 34 Figure 3: The distribution of estimated AAT from real perfusion data.

Figure 5 :

 5 Figure 5: The performance evaluation map with best sampling locations on real perfusion data.

Inria

  AAT quantification for real data

Figure 6 :

 6 Figure 6: The estimated CBF and AAT maps from real data corresponding to different diversity and repetition levels (M, R).

  Mean of SNR for AAT

Figure 7 :

 7 Figure 7: The performance evaluation map when the perfusion models used for generating data and used for quantification are mismatch.

Inria

  Figure8: The performance evaluation map obtained with best sampling locations and with approximated references on simulated perfusion data.

Table 1 :

 1 The best sampling locations. The relation between location index and sampling time is shown in Tab. 2.

	M	Best for CBF	Best for AAT	CBF & AAT
	2	5,10	2,9	2,9
	3	2,5,10	2,7,11	2,6,11
	4	2,6,8,11	2,5,7,11	2,6,7,11
	5	3,5,6,8,11	1,5,7,9,11	1,5,6,8,11
	6	3,6,8,9,10,11	1,5,6,7,9,11	1,6,7,9,10,11
	7	3,6,7,8,9,10,11	1,3,5,7,8,9,11	1,3,6,7,9,10,11
	8	3,5,6,7,8,9,10,11	1,3,4,6,7,8,9,11	1,3,6,7,8,9,10,11
	9	1,3,5,6,7,8,9,10,11	1,3,5,6,7,8,9,10,11	1,3,5,6,7,8,9,10,11
	10 1,2,3,5,6,7,8,9,10,11 1,3,4,5,6,7,8,9,10,11 1,2,3,4,6,7,8,9,10,11

Table 2 :

 2 Location index and its corresponding sampling time.

	Location Index Sampling Time (seconds)
	1	1.2
	2	1.3
	. . .	. . .
	11	2.2

Table 3 :

 3 Number of time locations selected as the best sampling locations.

	Time Locations 1 2 3 4 5 6 7 8 9 10 11
	For CBF	2 3 6 0 6 7 4 7 5 7	7
	For AAT	6 3 4 2 6 4 8 4 7 2	8
	CBF & AAT	6 4 4 1 2 8 6 4 6 5	8

RR n°8258

This value is chosen according to empirical result, referring Fig.

Inria

Inria

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399


Frédéric Nataf

Absorbing boundary conditions and perfectly matched layers in wave propagation problems

Keywords: artificial boundary condition, perfectly matched layer AMS classification. 35L05, 65Mxx

In this article we discuss different techniques to solve numerically wave propagation phenomena in unbounded domains. We present in a unified and simple way the two ways to restrict the computation to a finite domain: absorbing (or artificial) boundary conditions (ABC) and perfectly matched layers (PML). The intent is to give the possibility to the reader to grasp easily similarities and differences between these two truncation techniques. It should also allow the reader to adapt a truncation technique to the peculiarities of his physical modeling.

Introduction

To model numerically a wave propagation phenomena in a domain which extends to infinity, it is common practice to truncate the computational domain to a finite domain. The newly formed external boundary is somewhat artificial. To reduce undesired numerical reflections, there is a need for a special treatment at these boundaries. The solution must be as close as possible to the restriction of the solution that would have been computed in the original physical infinite domain. The most commonly used grid truncation techniques are so-called ABC, Artificial (or Absorbing) Boundary Conditions, and perfectly matched layer, PML, formulations. ABC technique was first formalized by B. Engquist and E. Majda in a seminal 1977 paper in [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF]. The PML concept was introduced by J.-P. Berenger in a 1994 paper, see [START_REF] Berenger | A Perfectly Matched Layer for the absorption of electromagnetic waves[END_REF]. ABC techniques are more general than PML. However, PML (which is an absorbing region rather than a boundary condition per se) can provide orders-of-magnitude lower reflections. Both methods are widely used in acoustic, electromagnetic or elastodynamic simulation software codes. Their precise formulations and efficiency depend strongly on the underlying partial differential equations that are discretized. For example, a stable and quite efficient PML for the non dissipative wave equation may surprisingly lead to unstable computations if used in conjunction with a damped wave equation. We present in this article, a unified and simple approach to the design and study of both ABC and PML formulations. The intent is to give the possibility to the reader to grasp easily similarities and differences between these two techniques. It should also allow the reader to adapt a truncation technique to the peculiarities of his physical modeling.

ABC

Homogeneous Dirichlet or Neumann boundary conditions for the wave equation lead to a total reflection of an impinging wave on the artificial boundary of a computational domain. We present here the method introduced by B. Engquist and A. Majda in 1977 (see [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF]) to design absorbing conditions that are easy to implement and yield a small reflection coefficient on the truncation boundary. The wave equation is solved in a domain exterior to an open bounded subset Ω of R 2 denoted by Ω c :

∂ 2 u ∂t 2 -c 2 ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 = f (t, x, y), (x, y) ∈ Ω c , t > 0 u(0, x, y) = u 0 (x, y) in Ω c ∂u ∂t (0, x, y) = u 1 (x, y) in Ω c u = 0 on ∂Ω u satisfies a radiation boundary condition at infinity (2.1)
The origin of the coordinate system is such that Ω is a subset of the half-space of negative x. We assume that the initial conditions u 0 and u 1 and the right hand side f have a compact support in this half space. The computational domain is truncated on its right by the artificial boundary x = 0, denoted Γ. It is not possible to use as an ABC a Dirichlet boundary condition on Γ since the solution is not known on Γ. The idea to design the ABC by finding a relation satisfied by the solution on Γ that does not depend on the values of the right hand sides. This integrodifferential relation defines an exact ABC which will be approximated by a partial differential equation.

Exact ABC

We first study the solution u in the right half space x ≥ 0 independently of the boundary condition satisfied at x = 0. We have:

∂ 2 u ∂t 2 -c 2 ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 = 0, x ≥ 0, y ∈ R, t > 0 u(0, x, y) = 0, x ≥ 0, y ∈ R ∂u ∂t (0, x, y) = 0, x ≥ 0, y ∈ R
In order to use the Fourier transform in time, we first prolongate u by zero for negative times. Since the initial conditions are zero for x ≥ 0, the prolongated function ũ still satisfies the homogeneous wave equation We take the Fourier transform in t and y (dual variables are ω and k respectively) of this equation:

∂ 2 ũ ∂t 2 -c 2 ∂ 2 ũ ∂x 2 + ∂ 2 ũ ∂y 2 = 0, x > 0, t, y ∈ R.
-ω 2 û -c 2 ∂ 2 û ∂x 2 + c 2 k 2 û = 0.
For a given (ω, k), it is an ordinary differential equation whose solution may be sought in the form û(ω, x, k) = j α j (ω, k)exp(λ j (ω, k)x). We have

λ 2 = k 2 - ω 2 c 2 .
We have two possibilities:

• evanescent waves, |k| > ω c : λ ± = ± k 2 -ω 2 c 2 • propagative waves, |k| < ω c : λ ± = ±i ω c 1 - c 2 k 2 ω 2 (i 2 = -1
) By using inverse Fourier transform, the general form for ũ is:

ũ(t, x, y) = |k|> ω c α -(ω, k)e -k 2 -ω 2 c 2 x e i(ωt+ky) dk dω (2.2) + |k|< ω c α -(ω, k)e i(ωt-ω c 1-c 2 k 2 ω 2 x+ky) dk dω (2.
3)

+ |k|> ω c α + (ω, k)e k 2 -ω 2 c 2
x e i(ωt+ky) dk dω (2.4)

+ |k|< ω c α + (ω, k)e i(ωt+ ω c 1-c 2 k 2 ω 2 x+ky) dk dω (2.5)
Next, we show that α + (ω, k) ≡ 0 for all (ω, k). In (2.4) indeed, the term e

k 2 -ω 2 c 2 x
tends to infinity as x tends to infinity. Solution ũ is bounded, we must have α + (ω, k) = 0 for |k| > ω c . Moreover, integral (2.5) corresponds to a linear combination of waves F. Nataf propagating in the direction of negative x. It is impossible that waves originate from infinity, only the boundary Γ creates waves. As a result, we have α + (ω, k) = 0 for |k| < ω c . Thus we have:

ũ(t, x, y) = |k|> ω c α -(ω, k)e -k 2 -ω 2 c 2
x e i(ωt+ky) dk dω

+ |k|< ω c α -(ω, k)e i(ωt-ω c 1-c 2 k 2 ω 2 x+ky) dk dω (2.6)
By taking x = 0 in this expression, we see that α -(ω, k) = û(ω, 0, k). By differentiating (2.6) with respect to x, we have thus at x = 0

∂ ũ ∂x (t, 0, y) - û(ω, 0, k)λ -(ω, k)e i(ωt+ky) dk dω = 0. (2.7)
where

λ -(ω, k) =        -k 2 - ω 2 c 2 for |k| > ω c -i ω c 1 - c 2 k 2 ω 2 for |k| < ω c .
We assume that relation (2.7) can be used as a boundary condition. It is then an exact ABC since solution u to problem (2.1) has a unique solution and satisfies (2.7) as well.

Approximation of the exact ABC

The exact ABC (2.7) involves inverse and direct Fourier transforms in both time t and y through the multiplication by a non polynomial function in the Fourier variables ω and k. Due to the square root in the formula for λ -, the inverse Fourier transform of (2.7) is not a differential equation. Its implementation would be thus costly. Moreover, for a varying velocity field c(x, y) or a curved boundary, the previous computation is not valid and the analytic form for the exact absorbing boundary condition is not known. For these reasons, integrodifferential equation (2.7) is approximated by a partial differential equation. They are easy to implement and have a straightforward extension to variable coefficients and/or curved boundaries.

In the Fourier space, (2.7) reads

∂ û ∂x (ω, 0, k) -λ -(ω, k) û(ω, 0, k) = 0. (2.8)
For reasons that will be given later (see § 3), small values of c k ω correspond to waves whose propagation direction is close to the normal to the artificial boundary, see figure 2. We approximate (2.8) by a Taylor expansion in c k ω of order zero and centered at zero:

∂ û ∂x (ω, 0, k) + i ω c 1 - c 2 k 2 ω 2 û(ω, 0, k) = 0 ≈ ∂ û ∂x (ω, 0, k) + i ω c û(ω, 0, k) = 0.
The inverse Fourier transform yields a zero-order approximation to the exact ABC (2.7): ∂u ∂t (t, 0, y) + c ∂u ∂n (t, 0, y) = 0 (zero-order ABC) (2.9)

where we replaced x by the outward normal to the boundary n in order to have a more intrinsic form . This is a partial differential equation which enables the dropping of the superscript˜. It is easy to check that (2.9) is an exact ABC for the one dimensional constant coefficient wave equation.

A higher order ABC is obtained with a Taylor expansion of order two (

√ 1 -2 ≈ 1 - 2 2 ): ∂ û ∂x (ω, 0, k) + i ω c 1 - c 2 k 2 ω 2 û(ω, 0, k) = 0 (2.10) ≈ ∂ û ∂x (ω, 0, k) + i ω c û(ω, 0, k) + c k 2 2iω û(ω, 0, k) = 0. (2.11)
The term c k 2 2iω is not polynomial in ω and thus does not correspond (under inverse Fourier Transform) to a partial differential equation on the artificial boundary. In order to fix this difficulty, we have two possibilities. The first one consists simply in multiplying relation (2.11) by iω before taking an inverse Fourier transform:

1 c ∂ 2 u ∂t∂n + 1 c 2 ∂ 2 u ∂t 2 - 1 2 ∂ 2 u ∂τ 2 = 0 (Order 2) (2.12)
where τ is the tangential direction on the artificial boundary. The second possibility is to introduce an auxiliary function φ = c k 2 2iω û(ω, 0, k) or equivalently, iω φ = c k 2 2 û(ω, 0, k). We take the inverse Fourier transform and get:

     1 c ∂u ∂t + ∂u ∂n + φ = 0 - 1 2
∂ 2 u ∂τ 2 = 1 c ∂φ ∂t (Order 2) (2.13)
For φ = 0, this is equivalent to condition (2.9). Condition (2.13) is thus an incremental modification of (2.9) and is easier to implement than (2.12), see [START_REF] Robert | Absorbing boundary conditions for difference approximations to the multidimensional wave equation[END_REF] for more details.

In order to get higher order ABCs a natural idea consists in using Taylor expansions of higher and higher orders. For instance, one could use a fourth order expansion

√ 1 -2 ≈ 1 - 2 2 - 4 8
. The corresponding ABC would read in the Fourier space:

∂ û ∂x (ω, 0, k) + i ω c û(ω, 0, k) + c k 2 2iω û(ω, 0, k) + c 3 k 4 8iω 3 û(ω, 0, k) = 0. (2.14)
As a matter of fact, the boundary value problem for the wave equation with the inverse Fourier transform of this boundary condition is not well-posed. It does not satisfy

F. Nataf
Kreiss condition [START_REF] Kreiss | Initial boundary value problems for hyperbolic systems[END_REF] or Shapiro-Lopatinski condition. In order to write stable high order ABCs it was proposed in [START_REF] Engquist | Absorbing Boundary Conditions for the Numerical Simulation of Waves[END_REF] to use Padé approximations:

1 -2 ≈ 1 - 3 2 4 1 - 2 4
, see [START_REF] Robert | Absorbing boundary conditions for difference approximations to the multidimensional wave equation[END_REF] as well for more details.

Plane waves analysis of an ABC

The analysis consists in computing for the various interface conditions the reflection of a plane wave as a function of its angle of incidence with respect to the normal to the interface. This gives an intuitive insight into the behavior of the ABCs we have introduced in the previous section. It enables to quantify the superiority of the ABC of order 2 (2.12) over the zeroth order one (2.9). It is also justified by the fact both the incident wave and the scattered one can be decomposed into a linear combination of plane waves as in formula (2.2)-(2.5). The same kind of analysis will be performed in section 5 for the PML method. The incident wave is written as follows:

θ Artifical boundary

B(u)=0

Onde incidente

Incident wave

Figure 2. Incident wave with an angle θ with respect to the normal to the artificial boundary

u I = e -i(ωt - ω c 1 - c 2 k 2 ω 2 x + ky) .
The angle of incidence θ with respect to the normal of the artificial boundary is such that:

tan θ = ck ω 1 -c 2 k 2 ω 2 and cos θ = 1 - c 2 k 2 ω 2 .
The reflected wave reads

u R = R(θ)e -i(ωt + ω c 1 - c 2 k 2 ω 2 x + ky) .
Function u I + R(θ)u R solves the homogeneous wave equation. The reflection coefficient R(θ) is computed by inserting the expression of u I + R(θ)u R in conditions (2.9) or (2.12).

For the zero-th order ABC (2.9), we have

R 0 (θ) = cos θ -1 cos θ + 1 ; R < 10 -2 for |θ| ≤ 11.45 •
and for the second order ABC (2.12):

R 2 (θ) = - cos θ -1 cos θ + 1 2 ; R < 10 -2 for |θ| ≤ 35 • .
The reflection coefficient for the first order ABC is smaller than one percent for a cone of angle 11.45 • whereas the reflection of the second order reaches this magnitude for a larger cone of angle 35 • . The latter ABC is sometimes referred to as a wide angle ABC.

PML

Historically, the first idea was to border the computational by a dissipative zone where the waves are damped. The dissipative equations were imposed a priori. The problem with this approach lies in the interface between the physical zone and the "sponge" layer. It led to reflection. This matching problem was first solved by Berenger in 1994 [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. He introduced perfectly matched layers (PML) for use with Maxwell's equations, and since that time there have been several related reformulations of PML for both Maxwell's equations and for other wave equations. Berenger's original formulation is based on a splitting of the electromagnetic fields into two unphysical fields in the PML region. In [START_REF] Chew | A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates[END_REF] PMLs were shown to correspond to a coordinate transformation in which the coordinate normal to the artificial boundary is mapped to complex numbers. This is actually an analytic continuation of the wave equation into complex coordinates, replacing propagating (oscillating) waves by exponentially decaying waves. The idea also consists in imposing the behavior of the solution in the PML and then retrieve the corresponding equations and interface conditions. In this way, it is possible to design a reflectionless PML. Thus key property of a PML that distinguishes it from an ordinary absorbing material is that it is designed so that waves incident upon the PML from a non-PML medium do not reflect at the interface. This property allows the PML to strongly absorb outgoing waves from the interior of a computational region without reflecting them back into the interior.

F. Nataf

Helmholtz equation

Let f : R 2 → C be a function with a compact support in the half-space of negative x, c > 0 the speed of sound and ω a non zero real number. We start with a Helmholtz problem defined in the whole space with f as the right hand-side:

(-ω 2 -c 2 ∆)(u) = f in R 2
with Sommerfeld radiation condition at infinity. We consider a PML that will be located in the positive x half-space. Let u 0 be the solution at x = 0, u 0 (y) := u(0, y).

We know from formula (2.6) that the solution to the Helmholtz equation in the right half-space x > 0 has the form:

u(x, y) = R û0 (k) exp(-i ω c 1 - c 2 k 2 ω 2 x) exp(iky) dk (4.1)
The PML must damp the propagative modes of the solution which correspond to the zone c|k| ≤ |ω|. For this, we first consider an analytic continuation of the solution with respect to the x coordinate. The complex variable will be denoted by z ∈ C. The analytic continuation is formally defined as:

ũ(z, y) = R û0 (k) exp(-i ω c 1 - c 2 k 2 ω 2 z) exp(iky) dk (4.2)
We now make a change of coordinate w.r.t to the first coordinate. Let σ(x) be a function such that σ(x) ≥ 0 for x ≥ 0. We define a path in the complex plane s : R → C:

s(x) := x + c iω x 0 σ(a) da (4.3)
and the function U : R 2 → C:

U (x, y) :=      u(x, y) if x < 0 ũ(s(x), y) if x > 0
We first investigate the behavior of U (x, y) for positive x.

Let û(x, k) := û0 (k) exp(-i ω c 1 - c 2 k 2 ω 2 x) .
We have with this notation:

U (x, y) = R û(x, k) exp -1 - c 2 k 2 ω 2 x 0 σ(a) da exp(iky) dk
We see that U is exponentially decreasing for c|k| < |ω| as x goes to infinity. The propagative modes of u have been turned into vanishing modes. As for the vanishing modes of u (c|k| > |ω|), they are still exponentially decreasing as x goes to infinity although now with a sinusoidal modulation. In summary, U contains only vanishing modes as x → ∞.

We now have to find out the equations satisfied by U . First of all, since U (x, y) = u(x, y) for x < 0, we have:

(-ω 2 -c 2 ∆)(U ) = f if x < 0 . (4.4)
We now determine the equation satisfied by u in the right half-space x > 0. We have:

∂U ∂x (x, y) = s (x) ∂ ũ ∂z (s(x), y) .
This yields:

1 s (x) ∂U ∂x (x, y) = ∂ ũ ∂z (s(x), y) .
Derivating again with respect to x, we get:

1 s (x) ∂ ∂x ( 1 s (x) ∂U ∂x )(x, y) = ∂ 2 ũ ∂z 2 (s(x), y) .
As for the y derivatives, we have for any integer m:

∂ m U ∂y m (x, y) = ∂ m ũ ∂y m (s(x), y) .
From formula (4.2), ũ satisfies the equation:

(-ω 2 -c 2 ∂ 2 ∂z 2 -c 2 ∂ 2 ∂y 2 )( ũ) = 0 .
Therefore, the last three equations yield the equation on U on x > 0:

(-ω 2 -c 2 1 s (x) ∂ ∂x ( 1 s (x) ∂ ∂x ) -c 2 ∂ 2 ∂y 2 )(U ) = 0 if x > 0 . (4.5)
We introduce the x-pml derivative:

∂ pml x := 1 s (x) ∂ ∂x = iω iω + cσ(x) ∂ ∂x .
We now turn to the interface conditions at x = 0. Since U satisfies a second order equation both half-spaces, we need two matching conditions on the interface x = 0.

F. Nataf

By definition of U and since s(0) = 0, we have that U is continuous along the interface x = 0: U (0 -, y) = U (0 + , y) . (

We consider now the derivative with respect to x. We know that

∂ pml x U (0, y) = ∂ ũ ∂z (s(0), y) .
By formula (4.2), we have ∂ ũ ∂z (s(0), y) = ∂u ∂x (0, y). Finally, since s(0) = 0, we get:

∂u ∂x (0, y) = ∂ pml x U (0, y) ,
which yields the interface condition at x = 0 for U

∂ x U (0 -, y) = ∂ pml x U (0 + , y) . (4.7)
In practice, it is necessary to consider a bounded PML zone, for instance a vertical slab (0, δ) × R. On the new boundary x = δ, we have to impose a boundary condition, for instance a homogeneous boundary condition. Finally, we get the following boundary value problem:

(-ω 2 -c 2 ∆)(v) = f in R -× R (-ω 2 -c 2 iω iω + cσ(x) ∂ ∂x ( iω iω + cσ(x) ∂ ∂x .) -c 2 ∂ 2 ∂y 2 )(v) = 0, 0 < x < δ v(0 -, y) = v(0 + , y) ∀y ∈ R ∂ x (v)(0 -, y) = iω iω + cσ(x) ∂ ∂x (v)(0 + , y) ∀y ∈ R v(δ, y) = 0 ∀y ∈ R (4.8)

The wave equation

In order to design a PML for the wave equation, we should take the inverse Fourier transform of the PML of the Helmholtz equation (4.8). The terms in (iω + cσ(x)) -1 correspond to non differential operators. The way to bypass this difficulty is, similarly to what is done in (2.13), to introduce two auxiliary unknowns in the PML region: ŵ(ω, x, y) and φ(ω, x, y) defined as:

ŵ(ω, x, y) = iω iω + cσ(x) ∂ û ∂x (ω, x, y) φ(ω, x, y) = iω iω + cσ(x) ∂ ŵ ∂x (ω, x, y) (4.9)
or equivalently:

(iω + cσ(x)) ŵ(ω, x, y) = iω ∂ û ∂x (ω, x, y) (iω + cσ) φ(ω, x, y) = iω ∂ ŵ ∂x (ω, x, y) (4.10)
At this point, an inverse Fourier transform in time yields partial differential equations. The time dependent counterpart of equation (4.8) reads:

∂ 2 v ∂t 2 -c 2 (x, y)∆v = f (t, x, y), t > 0, x < 0, y ∈ R ∂w ∂t + cσ(x)w = ∂ 2 v ∂t∂x , t > 0, 0 < x < δ, y ∈ R ∂φ ∂t + cσ(x)φ = ∂ 2 w ∂t∂x , t > 0, 0 < x < δ, y ∈ R - 1 c 2 (x, y) ∂ 2 v ∂t 2 + φ - ∂ 2 v ∂y 2 = 0, t > 0, 0 < x < δ, y ∈ R v(t, 0 -, y) = v(t, 0 -, y) t > 0, y ∈ R ∂v ∂x (t, 0 -, y) = w(t, 0 + , y) t > 0, y ∈ R (4.11)
together with initial conditions at t = 0: v(0, x, y) = u 0 (x, y), ∂v/∂t(0, x, y) = u 1 (x, y) and since u 0 has compact support in R -× R, we can take zero initial conditions on the auxiliary unknowns in the PML region w and φ.

Computation of the reflection coefficient of a PML

The analysis consists in computing the reflection of a plane wave as a function of its angle of incidence with respect to the normal to the interface. This gives an intuitive insight into the behavior of the PML and the impact of the damping parameter σ and the thickness of the PML. It is also justified by the fact that both the incident wave and the scattered one can be decomposed into a linear combination of plane waves as in formula (2.2)-(2.5). The same kind of analysis was performed in section 3 for the ABCs. It enables thus to compare quantitatively ABCs with the PML technique. In this section, the function σ takes a constant value still denoted by σ. As in (4.8), the Helmholtz domain is the left half-space x < 0. The PML is a band of width δ, (0, δ) × R. We compute the reflection coefficient for a PML with a homogeneous Dirichlet at the boundary x = δ of the PML. We make a planar wave analysis by considering homogeneous equations but instead of having a radiation condition at minus infinity, we consider a planar wave "coming" from minus infinity. Let

λ H = iω c 1 -c 2 k 2 ω 2 λ pml H = (1 + cσ iω )λ H
It can be checked that the functions e ±λ H x are solutions to the homogeneous Helmholtz equations and e ±λ pml H x are solutions to the homogeneous PML equations. We shall make use of the following formulas:

   ∂ x (e λ H x ) = λ H e λ H x ∂ pml x (e λ pml H x ) = iω iω + cσ λ pml H e λ pml H x = λ H e λ pml H x (5.1)
The solution u H to the homogeneous Helmholtz equation in the left domain has the form ûH = α H e λ H x + β H e -λ H x and in the PML zone:

ûpml = α pml e λ pml H x + β pml e -λ pml H x
Using (5.1), the coupling conditions (4.6)-(4.7) at the interface between the Helmholtz and the PML yield:

α H + β H = α pml + β pml λ H (α H -β H ) = λ H (α pml -β pml ) (5.2)
The homogeneous Dirichlet boundary condition in the PML at x = δ yields:

α pml e λ pml H δ + β pml e -λ pml H δ = 0 (5.3)
We study the damping of the PML by considering β H to be given. This corresponds to an ingoing wave from the physical (in space) medium and moving towards the interface between the physical medium and the PML medium. The three other quantities (α H , α pml , β pml ) are determined by the interface conditions. The reflection caused by a finite width PML is given by the ratio R := |α H /β H |. A direct computation shows that:

R = e -2 1-c 2 k 2 ω 2 σδ (5.4)
Using the physical meaning of 1 -c 2 k 2 ω 2 as the cosine of the angle θ between the ingoing wave and the normal to the boundary, R reads: R = (e -2σδ ) cos(θ)

(5.5)

We see here the main advantage of the PML approach over the ABC approach. In order to decrease the reflection coefficient, it is sufficient to increase the width of the PML and/or to take a large damping coefficient σ whereas for the ABC approach, it is necessary to implement more and more complex absorbing boundary conditions. At this point, one could think that it is enough to take a very large damping parameter σ. Actually, it creates a numerical reflection that is proportional to σ∆x. A remedy consists in taking a variable damping σ(x) = Cx 2 /δ and tune the parameters C and δ. In practice, δ is roughly a wavelength 2πc/ω and C is 20 -40. Of course, taking a smaller ∆x (finer mesh) will improve the results as well since we come closer to the continuous limit.

Conclusion

We have designed absorbing boundary conditions and perfectly matched layers for some scalar second order partial differential operators: the Helmholtz equation and the wave equation. The construction of the ABC is based on Fourier symbol computations of an exact ABC which is then approximated by a rational fraction. By introducing auxiliary unknowns, this approximation yields, in the physical space, a system of partial differential equations. This technique is quite general and can be applied to systems of partial differential equations, like the isotropic or anisotropic elasticity system with or without anisotropy, Oseen equations (linearized Navier-Stokes equations), Maxwell system, . . .. As for PML, their designs rely on a complex coordinate coordinate change that will be stable in the sense that propagative modes are turned into evanescent ones and that evanescent modes remain evanescent. For complex systems of PDEs with anisotropy, it is not always easy to find such a change of coordinates. Thus, although PMLs are usually better fitted to numerical computations than ABCs, their designs are not always feasible and are still an active area of research.

We have considered the design of absorbing boundary conditions and perfectly matched layers in a common Fourier framework and thus restricted ourselves to rectilinear boundaries for scalar equations. Actually, these concepts are also valid for curved boundaries and/or systems of partial differential equations, see for example [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF][START_REF] Givoli | Numerical methods for problems in infinite domains[END_REF][START_REF] Barucq | Micro-differential boundary conditions modelling the absorption of acoustic waves by 2D arbitrarily-shaped convex surfaces[END_REF][START_REF] Tsynkov | Artificial boundary conditions for the numerical simulation of unsteady acoustic waves[END_REF][START_REF]Exact nonreflecting boundary condition for elastic waves[END_REF][START_REF] Grote | Nonreflecting boundary conditions for Maxwell's equations[END_REF] and references therein.

An appealing feature of ABCs and PMLs is that these notions are building blocks for apparently unrelated applications such as domain decomposition methods [START_REF] Benamou | A domain decomposition method for the Helmholtz equation and related optimal control problems[END_REF], approximate factorizations preconditioners for Helmholtz type problems [START_REF] Engquist | Sweeping preconditioner for the Helmholtz equation: moving perfectly matched layers[END_REF] or inverse problems via time reversal techniques [START_REF] Assous | Time Reversed Absorbing Condition: Application to inverse problem[END_REF]. From a general point of view, it is due to the fact that ABCs or PMLs are cheap approximations to non local operators.

Figure 1 .

 1 Figure 1. Truncation of the computational domain to a half-plane

Author information

Frédéric Nataf, Laboratory J.L. Lions, UPMC and CNRS, Paris, France. E-mail: nataf@ann.jussieu.fr


Anthony Bourrier

Rémi Gribonval

Patrick Pérez

COMPRESSIVE GAUSSIAN MIXTURE ESTIMATION

Keywords: Gaussian mixture estimation, compressive sensing, database sketch, compressive learning

When fitting a probability model to voluminous data, memory and computational time can become prohibitive. In this paper, we propose a framework aimed at fitting a mixture of isotropic Gaussians to data vectors by computing a low-dimensional sketch of the data. The sketch represents empirical moments of the underlying probability distribution. Deriving a reconstruction algorithm by analogy with compressive sensing, we experimentally show that it is possible to precisely estimate the mixture parameters provided that the sketch is large enough. Our algorithm provides good reconstruction and scales to higher dimensions than previous probability mixture estimation algorithms, while consuming less memory in the case of numerous data. It also provides a privacy-preserving data analysis tool, since the sketch doesn't disclose information about individual datum it is based on.

INTRODUCTION

When processing a data collection for learning or for other tasks, it is often useful to fit a probability model to the data to get a more compact representation of the data. The fitted model can be used to represent data in a concise way, to feed learning algorithms that work on densities, to extract features or, simply, to uncover underlying structures.

However, fitting a mixture model often requires extensive access to the data to estimate parameters through the iterations of an optimization algorithm. Depending on the size of the data, this may not be possible: to compute the estimation while meeting the memory and computational limits of the hardware, one must typically randomly subsample the data.

In this paper, we consider the idea of computing a sketch of the data, which we define as a concatenation of empirical moments computed in one pass over the data. Given this compressed representation, we wish to estimate a model that fits the data. Figure 1 illustrates this idea. Related works that define sketches of data collections are discussed in section 2.

In order to solve our estimation problem, we search for a model which has a sketch close to the sketch of the data collection. Section 3 is dedicated to the definition of our problem with a special focus on Gaussian Mixture Models (GMM). To solve it, we exploit the analogy of the problem to compressive sensing and derive in section 4 an algorithm analog to Iterative Hard Thresholding (IHT) [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF].

We show experimentally in section 5 that it is possible to estimate precisely the GMM density in some cases which will be described later, with a lower memory consumption than EM. We give conclusions and discuss outlooks in section 6.

n x1 . . . xN Learning set (size N n) A =⇒ m ẑ Database sketch (size m) L =⇒ K θ Learned parameters (size K)
Fig. 1. Illustration of the proposed framework. The learning set is compressed to a sketch ẑ by a sketching operator A and the mixture parameters θ are estimated by a learning algorithm L using only the sketch. Intuitively, we should have K ≤ m for the algorithm to work and m N n for the method to allow computational and memory gains.

RELATED WORK

Algorithms using sketches can be found in the database literature [START_REF] Gilbert | How to summarize the universe: Dynamic maintenance of quantiles[END_REF][START_REF] Cormode | An improved data stream summary: The count-min sketch and its applications[END_REF]. In this case, a sketch is a compressed representation of the data and can be updated whenever an element is added or removed from the database without reprocessing all the data. A popular application is the search for frequent items in a data stream, also called heavy hitters [START_REF] Cormode | Methods for finding frequent items in data streams[END_REF].

Histogram fitting using a sketch has been explored [START_REF] Thaper | Dynamic multidimensional histograms[END_REF] in the case of n-dimensional discrete vectors. In this case, the sketch is an accumulated random projection of the vectors. This method does not scale to high dimensions, the construction of the histogram from the sketch having complexity which is exponential in n.

Inverse problems on density mixtures have also been studied [START_REF] Bunea | Spades and mixture models[END_REF][START_REF] Bertin | Adaptive Dantzig density estimation[END_REF]. Given data drawn from a mixture of candidate functions, both papers propose to cast the reconstruction as the optimization of a sparsity-inducing cost function on the vector of mixture coefficients. Both methods require the considered set of candidate densities to be finite and the elements of this set to be incoherent, i.e., different from each other. These assumptions do not hold for GMM: the centroids of the Gaussians can vary continuously and be arbitrarily close to one another.

Compressed representations of data vectors based on random projections for linear classification have been studied in [START_REF] Calderbank | Compressed learning : Universal sparse dimensionality reduction and learning in the measurement domain[END_REF][START_REF] Maillard | Compressed least-squares regression[END_REF]. These compressed representations aim at replacing a vector x by Mx where M is a dimensionality-reducing matrix, but do not compress the whole data set to a size that is independent of the number N of vectors in the set.

In [START_REF] Pitkow | Compressive neural representation of sparse, highdimensional probabilities[END_REF], a compressed representation of sparse probability distributions over multidimensional binary vectors is studied. The com-pressed representation is obtained by measuring the underlying distribution through random perceptrons and shown to carry enough information for one to be able to reconstruct the initial distribution with a 1 minimization algorithm.

PROBLEM STATEMENT

Let F ⊂ L 1 (R n ) be a set of probability densities, i.e., ∀f ∈ F, f is positive and ||f ||1 = R n f (x)dx = 1. Let X = {x1, . . . , xN } ⊂ R n be vectors drawn i.i.d. from a mixture p = k s=1 αsfs, with ∀s, fs ∈ F , αs ≥ 0 and k s=1 αs = 1. We want to find αs and fs from X .

Given σ > 0, we consider the case where F is a family of isotropic Gaussian densities of covariance matrix σI, that is:

F = fµ : x → 1 (2π) n 2 σ n exp x -µ 2 2 2σ 2 , µ ∈ R n . ( 1
)
Given m frequencies ω1, . . . , ωm, the sketching operator A is defined as:

F → C m (2) q → (h1(q), . . . , hm(q)) , (3)
where hj(q) = R n q(x)e -i x,ω j dx and F is the linear span of F . The sketch of a density mixture q is thus obtained by sampling the characteristic function of q at m frequencies. Given the empirical data X drawn from p, an empirical estimate of Ap, noted ẑ, is given by ẑ = ĥ1(X ), . . . , ĥm(X ) with ĥj(X ) = 1 N N i=1 e -i x i ,ω j . Reconstructing p as a k-sparse mixture of functions in F can be cast as the following minimization problem:

p = argmin q∈Σ k ẑ -Aq 2 2 , (4)
where Σ k := {f = k s=1 αsfs|∀s, fs ∈ F ∧ αs ≥ 0}. Note that we do not enforce the condition k s=1 αs = 1. In practice, we will see that the parameters can still be retrieved only with the positivity constraint.

This minimization problem is reminiscent of compressive sensing, where one wants to reconstruct a sparse finite-dimensional vector from a number of measures that is lower than its dimension. Typically, these measures are linear combinations of the entries of the vector. The difference is that here, the vector is infinite-dimensional and is sparse in a basis (F ) which is extremely coherent, i.e., composed of vectors that can be very close to each other (if µ ≈ ν, then fµ ≈ fν ).

ALGORITHM

To address the estimation problem (4), we propose an algorithm analogous to Iterative Hard Thresholding (IHT) [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF].

Iterative Hard Thresholding

IHT is designed to reconstruct a k-sparse signal x of dimension n from m measurements given by y = Mx, where M is a m × n matrix (with m < n). At each iteration, IHT updates an estimate x of x, decreasing the objective function φ : x → 1 2 y -Mx 2 2 while ensuring the k-sparsity of x. This is performed in two steps:

1. The n-dimensional gradient of φ, noted g, is computed.

2. The update is given by x ← H k (x-λg), where λ is a descent step and H k is a hard thresholding operator which keeps only the k entries of the vector with largest module and sets the others to 0.

Continuous case

In our case, the signal to reconstruct is p and our estimate p is parametrized by a vector α ∈ R k of positive weights and by the support Γ = { μ1 , . . . , μk } ⊂ R n corresponding to the means of the current estimated Gaussians. The current residual is defined by r = ẑ -Ap. Some differences with IHT require a modification of the procedure. In IHT, the signal one wants to reconstruct is supposed sparse in a finite base of vectors, and the gradient of the objective function thus have a finite number of entries which measure the infinitesimal shift of the objective function if a certain entry of the vector is shifted. In our case, the density belongs to the infinite-dimensional space F , and there are as many directions in which the current estimation of the density can be shifted as the cardinality of R n . The equivalent of the gradient is the function:

g p(µ) = lim t→0 1 2t ẑ -A(p + tfµ) 2 2 -ẑ -Ap 2 = -Afµ, r . (5)
Instead of seeking k global extrema as in IHT, we seek a certain number M of local minima of g p by initializing randomly a gradient descent algorithm. The random initialization will be discussed in section 5. These local minima are added to the current support Γ and the hard thresholding is performed by selecting the highest k components of the projection of ẑ on the cone generated by Γ, i.e., the set of linear combinations of elements of Γ with positive coefficients.

Since the family F is extremely coherent, the local minima we find may be shifted from the true mean vectors because of the imprecision induced by the other components of the mixture. This is why we add a gradient descent step after the search for local minima: this way we reduce further the objective function after having found the local minima.

Algorithm 1 describes the overall procedure. It mainly consists of three steps by iteration:

1. M local minima of µ → g p(µ) are sought with a gradient descent algorithm with random initialization and are added to the current support.

2. The sketch ẑ is projected on this support with a positivity constraint on the coefficients. Only the k highest coefficients and the corresponding vectors of the support are kept.

3. A gradient descent algorithm is applied to further decrease the objective function with respect to the weights and support vectors.

Algorithms 2, 3, 4 and 5 describe the subfunctions.

Memory usage

Let's consider that n, k and m are much larger than For i = 1 to M Find a local minimum νi of the function: ν ∈ R n → -Afν , r with a gradient descent algorithm, initialized randomly. End For Return ν1, . . . , νM . memory requirement of EM is O(N (k+n)) to store both the vectors and their probabilities to belong to each current component of the mixture. The compressed algorithm allows to memory savings as soon as m + kn k+n N . Since kn m, this condition is nearly equivalent to m N .

EXPERIMENTS

Experimental setup

To evaluate the behavior of our algorithm, we conducted experiments on vectors drawn from a mixture of k isotropic Gaussians with identity covariance matrices (σ = 1). In each case, we drew weights uniformly on the simplex 1 , and we chose the Gaussian means by drawing random vectors, each entry being drawn from a probability law of density N (0, 1).

The experiments were performed in the following way: after the choice of the probability distribution p, we drew N random vectors from this probability distribution and computed the empirical sketch of the distribution in one pass of the data, which was then discarded from hard memory (see section 5.2). We then applied the reconstruction algorithm to the sketch to get an approximated mixture p. The random initialization for the reconstruction algorithm is detailed in section 5.3.

To measure the quality of the reconstruction, we considered a symmetrized version of the Kullback-Leibler (KL) divergence between the true mixture and the estimated one, defined by D(p||p) + D(p||p), where D denotes the regular KL divergence. To approximate this divergence, we drew N = 10 5 points (yi) N i=1 i.i.d. from p and computed 1 N N i=1 ln p(y i ) p(y i ) + p(y i ) p(y i ) ln p(y i ) p(y i )

. We also considered an approximated Hellinger distance, measured by the quantity 1 -1 N N i=1 p(y i ) p(y i ) . The KL divergence ranges from 0 to +∞ while the Hellinger distance ranges from 0 to 1. In both 1 We also performed experiments where all the weights were equal to Find a local minimum (α , µ 1 , . . . , µ k ) of the function:

R k × (R n ) k → R+ (β, ν1, . . . , ν k ) → ẑ -[Afν 1 , . . . , Afν k ]β 2, using a gradient descent algorithm initialized at α, µ 1 , . . . , µ k . Return α , {µ 1 , . . . , µ k }.
cases, lower values mean closer distributions.

Choice of the frequencies

Note that if p is a mixture of probability densities taken in F , then its Fourier transform decays as exp -σ 2 2 ω 2 2 , giving an upper bound for the entries of the sketch.

Since we do not want the elements of the sketch to be too small (they should be informative about p), we chose the frequencies ωj randomly with a probability law N (0, 1 σ 2 Id) (thus proportional to the upper bound we found for the values of the sketch). Here, since σ = 1, this law was N (0, Id).

Heuristic for random initialization

The search for local minima in algorithm 2 was initialized randomly by exploiting a measure performed during the construction of the sketch: during the single pass on the data, the norms of the vectors are computed and the maximum of the norms R = max x∈X ||x||2 is computed. This measure has a negligible effect on the computation time of the sketch and delimits a ball in which the centers of the Gaussians are very probably contained.

We performed the random initialization by drawing a direction uniformly on the unit sphere and multiplying this unit vector by a scalar uniformly drawn in [0; R]. Table 1. Comparison between our compressed estimation algorithm and an EM algorithm in terms of precision of the estimation and memory usage (in megabytes). Experiments were performed with n = 20, k = 10, m = 1000. In each cell, the value is a median on 10 experiments with the standard deviation for the precision measures.

Results

is the Hellinger distance such that 80% of the experiments lead to a smaller Hellinger distance. We can essentially observe a gradually increasing measure of the Hellinger distance as the number of mixture components rises. For the considered parameters range, choosing m = 10kn, i.e., choosing m so that it contains 10 times more values than the number of parameters to estimate, leads to a Hellinger distance smaller than 0.03 for 80% of the cases.

Table 1 compares our algorithm with a standard EM algorithm in the case where n = 20, k = 10, m = 1000 for values of dataset size N ranging from 10 3 to 10 5 . For each case, we can see that the precision of the estimation increases with the number of samples. In the compressed case, this can be explained by the fact that the components of the sketch are better estimated with more points. We notice that the memory used for EM is proportional to the number N of samples in the dataset, while the memory required by the compressed algorithm does not depend on this parameter, which leads to a substantial improvement in memory usage for N ≥ 10 4 . Even with this reduced memory cost, the compressed algorithm is able to provide a precision comparable to the precision of the EM algorithm.

CONCLUSION AND OUTLOOKS

We proposed a framework for density estimation exploiting a sketch of the data instead of the data itself. This sketch is computed in one pass of the data, which can be discarded on the fly, leading to memory savings in the case of numerous data and to privacy-keeping estimation, since the individual data points cannot be recovered from the sketch. We detailed how this framework can be instantiated for the estimation of a mixture of isotropic Gaussians by deriving an algorithm by analogy with IHT. We experimentally showed that it provides precise reconstruction compared to standard EM while having memory requirements that do not depend on the number of data vectors. Future work includes generalization to non-isotropic Gaussians, richer choice of sketching functions and other choices of initialization methods for the search of local minima, as well as a theoretical understanding of the well-posedness and stability of the considered problems.

On the experimental side, it would be interesting to see how other methods for large-scale learning such as online EM [START_REF] Liang | Online EM for unsupervised models[END_REF] or stochastic gradient [START_REF] Bottou | The tradeoffs of large scale learning[END_REF] can also be compared with our framework.

1

 1

k

  and this didn't alter the conclusions drawn from the experiments.

Algorithm 3 Algorithm 4 Algorithm 5

 345 Proj cone(v, Γ = {u1 . . . uK }) Solve the following convex optimization problem: a = argmin β∈R K + ||v -Uβ||2, with U = [u1, . . . , uK ]. Return a. Hard threshold(a, Γ = {u1, . . . , uK }, k) Let ai 1 , . . . , ai k be the k highest entries of a. Return (ai 1 , . . . , ai k ), {ui 1 , . . . , ui k }. Shift support(ẑ, α, Γ = {µ 1 , . . . , µ k }, k)

Figure 2

 2 Figure 2 visually illustrates the behavior of the algorithm on a simple mixture of 4 Gaussians in dimension 2. N = 10 3 points were drawn from this mixture and used to compute a m = 30-dimensional sketch. As shown in the figure, the parameters are precisely estimated without referring to the initial data. The symmetric KL divergence and Hellinger distance are respectively 0.026 and 0.003.Figure 3 illustrates the reconstruction quality of our algorithm in dimension 10 for different values of mixture components k and sketch sizes m in terms of Hellinger distance. For each sketch size m ranging from 200 to 2000 with range 200, k was chosen to range from m/200 to m/10 with step m/200. For each choice of parameters, 10 experiments were performed and the depicted value

Figure 3

 3 Figure 2 visually illustrates the behavior of the algorithm on a simple mixture of 4 Gaussians in dimension 2. N = 10 3 points were drawn from this mixture and used to compute a m = 30-dimensional sketch. As shown in the figure, the parameters are precisely estimated without referring to the initial data. The symmetric KL divergence and Hellinger distance are respectively 0.026 and 0.003.Figure 3 illustrates the reconstruction quality of our algorithm in dimension 10 for different values of mixture components k and sketch sizes m in terms of Hellinger distance. For each sketch size m ranging from 200 to 2000 with range 200, k was chosen to range from m/200 to m/10 with step m/200. For each choice of parameters, 10 experiments were performed and the depicted value

Fig. 2 .

 2 Fig.2. Real and reconstructed centroids (respectively represented as circles and squares) of 4 Gaussians in dimension 2 from 10 3 points drawn from the mixture. To estimate the 12 real parameters of the mixture, the data was compressed to a complex-valued sketch of dimension 30, represented to the right as a 60-dimensional real signal.

Fig. 3 .

 3 Fig. 3. Quality of reconstruction in dimension n = 10, with N = 10 4 points, measured as a Hellinger distance (darker is better). Each square corresponds to 10 experiments, and the depicted values are the values of the Hellinger distance under which 80% of performed experiments are placed.

  Compressive isotropic Gaussian mixture parameter estimation Input: Sketch ẑ, operator A, target sparsity k, integer M . Initialize Γ = ∅, r = ẑ. repeat Set ν1, . . . , νM ← Find min(A, r, M ).

	Set Γ	← Γ ∪ {ν1, . . . , νM }.
	Set α	← Proj cone(ẑ, Γ ).
	Set α, Γ	← Hard threshold( α , Γ , k).
	Set α, Γ	← Shift support(ẑ, α, Γ, k).
	Set r	← ẑ -k j=1 αjAfμ j .
	until Stopping criterion is satisfied
	Return α, Γ.
	Algorithm 2 Find min(A, r, M )

1. If we suppose that optimization algorithms only use first-order quantities, their memory costs are dominated by O(kn). The computation of the cost function of Algorithm 5 requires O(km). The storage of the operator A (via the frequencies ωj) requires O(mn).

Therefore, the total memory usage is O((k + n)m + kn) and does not depend on the number N of vectors. In comparison, the Algorithm 1

Acknowledgments

This work was supported in part by the European Research Council, PLEASE project (ERC-StG-2011-277906).


Bassem Mouawad
email: bassem.mouawad@insa-lyon.fr

Maher Soueidan
email: bmaher.soueidan@insa-lyon.fr

Damien Fabrègue
email: cdamien.fabregue@insa-lyon.fr

Cyril Buttay
email: dcyril.buttay@insa-lyon.fr

Vincent Bley
email: evincent.bley@laplace.univ-tlse.fr

Bruno Allard
email: fbruno.allard@insa-lyon.fr

Mechanical Study of Copper Bonded at Low Temperature using Spark Plasma Sintering Process

Keywords: Copper bonding, Spark Plasma Sintering, Tensile test

Bonding of high purity polished copper was investigated using the Spark Plasma Sintering technique (SPS) showing the effect of SPS parameters (surface roughness, time, temperature and pressure) on the bonding strength behaviour. Mechanical characterization of the bonded samples was performed at room temperature using tensile test. Two surfaces roughnesses were studied (un-polished and polished samples). It was found that the bonding strength varied from 50 MPa to 233 MPa for un-polished and polished surfaces respectively The tensile strength of the used bulk copper-rod was found to be 365 MPa, while most results are over 122 MPa (a third of the bulk value).

Introduction

The work presented here is aimed at developing a new interconnect technology for power electronic devices. Such new technology must provide a very low electrical and thermal resistance. Furthermore, the manufacturing conditions must be compatible with the semiconductor dies (temperature, pressure) and with the industrial constraints (time). Therefore, we present here a technique in which pieces of copper are joined together at temperature lower than 300 °C, with a short turnaround cycle (30 min max) and a relatively low pressure (16-77 MPa).

Several experiments on Cu/Cu direct bonding have been reported [START_REF] Fan | Copper Wafer Bonding[END_REF][START_REF] Chen | Copper Bonded Layers Analysis and Effects of Copper Surface Conditions on Bonding Quality for Three-Dimensional Integration[END_REF][START_REF] Li | Low Temperature Copper to Copper Direct Bonding[END_REF]. Although the methods are simple, they require either high external pressure and annealing or toxic chemical cleaning processes to achieve a reasonable bonding strength. For example, Cu wafers exhibit good bonding properties when Cu/Cu contacts are carried out at 400 °C and 4 bar for 30 min, followed by an annealing cure at 400 °C for 30 min in N 2 atmosphere [START_REF] Chen | Microstruture Examination of Copper Wafer Bonding[END_REF]. Another solution, the so-called Surface Activated Bonding (SAB) process was used to bond Cu-coated wafers at room temperature without any wet chemical process. SAB is a process that joins similar or dissimilar materials by means of the adhesion force between atoms of two atomically clean surfaces in an ultrahigh vacuum UHV at room temperature [START_REF] Takagi | Surface activated bonding of silicon wafers at room temperature[END_REF][START_REF] Takagi | Room-Temperature Bonding of Si Wafers to Pt Films on SiO 2 or LiNbO 3 Substrates Using Ar-Beam Surface Activation[END_REF]. The same authors also reported that the bonding strength of Cu/Cu interface is approximatively 6.47 MPa [START_REF] Kim | Room temperature Cu-Cu direct bonding using surface activated bonding method[END_REF]. J. W. Elmer & al [START_REF] Elmer | Diffusion bonding and brazing of high purity copper for linear collider accelerator structures[END_REF] have presented a diffusion bonding of high purity copper using a conventional furnace. A series of diffusion bonds was done to determine the relationship between bond strength and bonding parameters in order to determine a range of bonding conditions. Bonding strengths of 65 MPa and 190 MPa were found using uniaxial tensile testing for samples bonded at 400 and 800°C respectively (bonding time 60 min). In general, high temperature annealing (>400°C) is required to increase the bonding strength. However, semiconductor devices have been designed to meet the requirement of the existing soldering technology (a few seconds above 300°C). Therefore, an increase in the process temperature could require heavy modification in the dies manufacturing parameters to avoid possible degradation. In this paper, we report the Cu/Cu direct bonding using Spark Plasma Sintering (SPS) where we focus on the bonding issue by evaluating the bonding strength as a function of temperature, surface roughness, pressure and bonding time. The samples were then cleaned as follows (i) ultrasonic cleaning in acetone to remove the lubricant, (ii) ultrasonic cleaning in trichloroethylene, (iii) ultrasonic cleaning in microelectronic grade ethanol, (iv) the samples were immersed in 5% HF for 1 min at room temperature for oxide elimination (v) followed by rinsing in DI water and blow drying with nitrogen gas.

Bonding Process. The Cu/Cu bonding was carried out with a Spark Plasma Sintering (SPS) system (FCT HPD 25). As with a classical sintering press, a SPS system applies both pressure and heat to the samples. The difference is that the heat is provided by a pulsed DC-current circulating through the graphite die (or through the sample if it is conductive), whereas the conventional sintering press uses external heating elements. As a result in the case of the SPS, the heat is generated much closer to the sample, offering very fast temperature dynamics: more than 100 K/min for the system used in this study, using current pulses of several kA. In each experiment, two copper pieces were introduced in a graphite die of 20 mm inner diameter (Fig. 1). The bonding was performed in vacuum (10 2 torr). Temperature, pressure and dwell time were set according to Fig. 2. The temperature was measured using a thermocouple through a hole in the graphite die, as shown in Fig. 1. The heating rate was fixed to 100 K/min and the holding time was varied from 5 to 20 minutes (Fig. 2).

Characterization. After bonding the samples were machined for the mechanical test (tensile test), as described in Fig. 3. Mechanical testing of the bonded samples was performed at room temperature using a MTS Adamel Lhomargy-DY35 with a 20kN cell force. The displacement rate was set up to 1mm/min. Copper being ductile, the sample quickly enters plastic deformation during the tensile test (Fig. 4). This results in a reduction of the actual section of the sample as it elongates. In this document, a nominal stress value is calculated using the initial surface of the tensile sample. It will thus constitute a lower bound of the true stress experimented on the sample. It should be mentioned that all the samples fractured along the bond line (due to the design of the sample). Surface microstructure observations were studied by a Leitz optical microscope and a scanning electron microscope (SEM Philips XL20).

Results and discussion

Surface roughness. A series of Cu/Cu bonding was made to determine the relationship between bond strength and bonding parameters such as bonding temperature, holding pressure and dwell time, in order to determine a range of bonding conditions that could be used to produce satisfactory joints. 1) For this purpose, two polished rods were bonded at 300°C for 20 minutes and a pressure of 16 MPa. The same conditions were used to bond together two un-polished rods. Both experiments resulted in a proper bond. On both assemblies, the bonding interface was not visible, obviously leading to the conclusion that SPS process has enabled to create a Cu/Cu joint with this set of parameters. Tensile tests, however, revealed two different behaviours for the polished and unpolished samples. In the case of un-polished sample, fracture occured without significant plasticity, and at a relatively low maximal force (50 MPa). For the polished sample, fracture occured after a lot of plastic deformation, exhibits a large ductility, and a high maximum force (233 MPa).

Temperature, time, pressure. Table 1 summarizes a list of the experiments, with the corresponding tensile strengths. All these experiments were performed with polished samples for different temperature, holding pressure and dwell time. The first eight experiments show a variation of time (5 and 20 min), holding pressure (16 and 77 MPa, which correspond to the pressure applied on a 20mm-diameter sample) and the bonding temperature (200 to 300 °C). Fig. 5 shows the influence of different parameters on the bonding strength.

The couples (t,T), (t,P), (T,P) and (t,T,P) express the conjugated effect of time, temperature and pressure between each other and they are given a normalized values (-1,1). It can be seen that time and pressure has the expected effect: an increasing in either of these parameters results in a stronger bond. The negligible influence of temperature and the negative effect of the simultaneous increase of time and temperature can firstly be surprising, but can actually be explained by the following considerations: in these tests, a good bond could be achieved for any temperature or time value. However, the application of temperature to the samples also constitutes a heat treatment of the copper. The effect of such treatment on work-hardened copper is a change in its microstructure. For example, increasing the temperature from 200 °C to 300 °C would involve further annealing recovery and recrystallization of the sample, decreasing its overall mechanical properties. Fig. 6 SEM image of the fracture surface (after tensile test), for a sample bonded at 300 °C for 20 minutes and with a pressure of 16 MPa) Even with the limited variation range for all three factors (5-20 minutes, 200-300 °C, and 16-77 MPa), substantial differences were found in tensile strength. In most cases, however, a very good bond was formed: the tensile strength of the used bulk copper-rod was found to be 365 MPa, while most results in table 1 are over 122 MPa (a third of the bulk value).

Another evidence of the good bond quality is visible in Fig. 6: the fracture surface is characteristic of ductile fracture with the presence of dimples. This means that the sample has been heavily deformed, which is also exhibited by the large strain value attained during the tensile tests (as can be seen for example in figure 4). Similar patterns were found over the entire fracture surface and for all samples.

Fig. 1 Fig. 2 .

 12 Fig.1The graphite die used to bond the copper rods, with the two 20 mmdiameter pressing cylinders

Fig. 3

 3 Fig. 3 Shape (a) and photograph (b) of the samples used for the tensile test Fig. 4 Tensile test result for one of the sample (sample 2 in table1) For this purpose, two polished rods were bonded at 300°C for 20 minutes and a pressure of 16 MPa. The same conditions were used to bond together two un-polished rods. Both experiments resulted in a proper bond. On both assemblies, the bonding interface was not visible, obviously leading to the conclusion that SPS process has enabled to create a Cu/Cu joint with this set of parameters. Tensile tests, however, revealed two different behaviours for the polished and unpolished samples. In the case of un-polished sample, fracture occured without significant plasticity, and at a relatively low maximal force (50 MPa). For the polished sample, fracture occured after a lot of plastic deformation, exhibits a large ductility, and a high maximum force (233 MPa).Temperature, time, pressure. Table1summarizes a list of the experiments, with the corresponding tensile strengths. All these experiments were performed with polished samples for different temperature, holding pressure and dwell time. The first eight experiments show a variation of time (5 and 20 min), holding pressure (16 and 77 MPa, which correspond to the pressure applied on a 20mm-diameter sample) and the bonding temperature (200 to 300 °C). Fig.5shows the influence of different parameters on the bonding strength.The couples (t,T), (t,P), (T,P) and (t,T,P) express the conjugated effect of time, temperature and pressure between each other and they are given a normalized values (-1,1). It can be seen that time and pressure has the expected effect: an increasing in either of these parameters results in a stronger bond. The negligible influence of temperature and the negative effect of the simultaneous increase of time and temperature can firstly be surprising, but can actually be explained by the following considerations: in these tests, a good bond could be achieved for any temperature or time value. However, the application of temperature to the samples also constitutes a heat treatment of the copper. The effect of such treatment on work-hardened copper is a change in its microstructure. For example, increasing the temperature from 200 °C to 300 °C would involve further annealing recovery and recrystallization of the sample, decreasing its overall mechanical properties.

Fig. 5

 5 Fig. 5 Influence of different parameters on engineering stress

Table 1

 1 Strength bonding and hardness as function of time and temperature for SPS copper bonding with polished surface

	Sample	Time [min]	Temperature [°C]	Pressure [MPa]	Engineering stress [MPa]
	1	5	300	77	181
	2	20	300	77	239
	3	20	300	16	170
	4	5	300	16	132
	5	5	200	77	106
	6	20	200	77	270
	7	20	200	16	261
	8	5	200	16	No bond (0)
	9	20	250	16	150
	10	20	250	16	141
	11	20	250	16	180
	12	20	250	16	169
	13	20	250	16	183

Summary

The high purity copper-to-copper bonding technique using SPS is presented. The un-polished sample showed considerably less bonding strength than the polished sample. A bonding strength higher than previously reported in the literature was found (between a third and more than two thirds of the bulk copper tensile strength). This bond is compatible with the requirements of the application: a semiconductor die can sustain temperatures up to 300 °C for a few minutes and pressures in the order of a few tens of MPa.