Dima Abi-Abdallah Rodriguez

Emmanuel Durand

Ludovic De Rochefort

Younes Boudjemline

Elie Mousseaux

Simultaneous pressure-volume measurements using optical sensors and MRI for left ventricle function assessment during animal experiment

Simultaneous pressure and volume measurements enable the extraction of valuable parameters for left ventricle function assessment. Cardiac MR has proven to be the most accurate method for volume estimation. Nonetheless, measuring pressure simultaneously during MRI acquisitions remains a challenge given the magnetic nature of the widely used pressure transducers. In this study we show the feasibility of simultaneous in vivo pressure-volume acquisitions with MRI using optical pressure sensors. Pressure-volume loops were calculated while inducing three inotropic states in a sheep and functional indices were extracted, using single beat loops, to characterize systolic and diastolic performance. Functional indices evolved as expected in response to positive inotropic stimuli. The end-systolic elastance, representing the contractility index, the diastolic myocardium compliance, and the cardiac work efficiency all increased when inducing inotropic state enhancement. The association of MRI and optical pressure sensors within the left ventricle successfully enabled pressure-volume loop analysis after having respective data simultaneously recorded during the experimentation without the need to move the animal between each inotropic state.

Introduction

The evaluation of the left ventricle (LV) performance is of high importance in physiologic investigation and clinical practice. The diastolic LV function can be assessed by evaluating relaxation properties and myocardium stiffness during the filling process. As for the systolic LV performance, it is governed by three principal factors, preload, afterload and the contractile or inotropic state of the myocardium. The first two determinants can be relatively easily measured or evaluated, but assessing contractile state is far more difficult. Any index intended to reflect contractility must be load-independent.

Common studies evaluate contractility using left ventricular stroke volume (SV), ejection fraction (EF) and cardiac output (CO). Though intuitive and simple, these parameters are load-dependent and consequently represent poor contractility indices. For instance, they are not good indicators of heart failure and cardiac dysfunction Maeder and Kaye [START_REF] Maeder | Heart failure with normal left ventricular ejection fraction[END_REF], Gaddam and Oparil [START_REF] Gaddam | Diastolic dysfunction and heart failure with preserved ejection fraction: rationale for raas antagonist/ccb combination therapy[END_REF]. Other studies rely on pressure measurements alone to assess the LV performance. They evaluate, for example, the maximum rate of pressure change (dP/dt max ) which is known to be sensitive to inotropic state and thus correlates with cardiac contractility. Nonetheless, its load-dependence Grossman et al. [START_REF] Grossman | Alterations in preload and myocardial mechanics in the dog and in man[END_REF], Kass et al. [START_REF] Kass | Comparative influence of load versus inotropic states on indexes of ventricular contractility: experimental and theoretical analysis based on pressure-volume[END_REF], Mason et al. [START_REF] Mason | Assessment of cardiac contractility: The relation between the rate of pressure rise and ventricular pressure during isovolumic systole[END_REF] makes it a weak contractility indicator. The peak decline of pressure (dP/dt min ), on the other hand, can quantify isovolumic relaxation, but it cannot be qualified as an intrinsic relaxation index since it is preload-dependent Cohn et al. [START_REF] Cohn | Maximal rate of pressure fall (peak negative dp/dt) during ventricular relaxation[END_REF] as well as afterload-dependent Weiss et al. [START_REF] Weiss | Hemodynamic determinants of the time-course of fall in canine left ventricular pressure[END_REF]. A widely used parameter for the relaxation quantification is the time constant τ which characterizes the pressure decay during isovolumic relaxation.

Despite being preload-independent Varma et al. [START_REF] Varma | Is tau a preload-independent measure of isovolumetric relaxation?[END_REF], τ cannot be considered an ideal index because of its afterload dependence.

No perfect assessment index can be obtained from volume or pressure alone. Indeed, volume or pressure measurements on their own are not sufficient to characterize the systolic performance, they cannot solely define contractility and cardiac response to inotropic agents. Simultaneous pressure and volume measurements are necessary to provide valuable functional parameters. Pressure indices can thus be coupled with volume information to negate load dependence. For example, rather than considering dP/dt max , the relationship of this value to end diastolic volume (EDV) can be computed by varying preload conditions. This relationship is linear and its slope provides a preload-independent contractility index Little [START_REF] Little | The left ventricular dp/dt max -end-diastolic volume relation in closed-chest dogs[END_REF]. The relaxation constant τ can also be evaluated over a range of afterloads and plotted against EDV Zile and Brutsaert [START_REF] Zile | New concepts in diastolic dysfunction and diastolic heart failure: Part I:diagnosis, prognosis, and measurements of diastolic function[END_REF].

In addition, a wide variety of indices that can be quantified by analyzing pressure-volume (PV) loops have been proposed to characterize the left ventricle systolic and diastolic performance. For instance a straightforward characterization of the myocardium stiffness can be obtained from the slope of the pressurevolume relationship during diastole. Furthermore, a measurement of the ventricular elastance when the contractile forces in the ventricle are at their peak, constitutes a good indicator of the ventricular contractility and systolic function. Known as end-systolic elastance E es , it is independent of preload, afterload and heart rate Brinke et al. [START_REF] Brinke | Single-beat estimation of the left ventricular end-systolic pressure-volume relationship in patients with heart failure[END_REF], Suga [START_REF] Suga | Global cardiac function: mechano-energetico-informatics[END_REF]. The arterial system, as afterload, can also be assessed from the PV loop, and, like the ventricle chamber, it can be characterized by its elastance E a . Studies have shown the importance of E a as a descriptor of the vascular load and its impact on cardiac performance Hayashida et al. [START_REF] Hayashida | Mechanical matching of the left ventricle with the arterial system in exercising dogs[END_REF], Kelly et al. [START_REF] Kelly | Effective arterial elastance as index of arterial vascular load in humans[END_REF], indicating that the ratio E a /E es quantifies the coupling between the ventricle and arterial system and governs ventriculoarterial matching. Additionally, PV loop area analysis can provide an evaluation of mechanical energies of a ventricular beat and an assessment of the LV efficiency.

To measure the ventricle volume, the conductance catheter technique has been used extensively, it is based on the electrical conductance of the blood contained in the cavity. This technique, nonetheless, is based on geometric assumptions, needs volume-dependent calibration and is limited by the non-linearity of the conductance-volume relation when large volume changes are involved Danton et al. [START_REF] Danton | Right ventricular volume measurement by conductance catheter[END_REF]. Furthermore, it relies on Ohm's law which might not be fully appropriate due to the non-uniformity of the composition of ventricle and blood, and uses correction parameters that are error prone due to ventricle geometry and wall thickness changes during contraction Jacoby et al. [START_REF] Jacoby | Direct comparison of magnetic resonance imaging and conductance microcatheter in the evaluation of left ventricular function in mice[END_REF]. Volume assessment using cine cardiac MR are proven to yield a more accurate and reliable estimate Jacoby et al. [START_REF] Jacoby | Direct comparison of magnetic resonance imaging and conductance microcatheter in the evaluation of left ventricular function in mice[END_REF], Winter et al. [START_REF] Winter | Left ventricular function in the post-infarct failing mouse heart by magnetic resonance imaging and conductance catheter: a comparative analysis[END_REF]. For PV loop construction, simultaneous volume and pressure measurements are necessary, but given the magnetic nature of standard pressure sensors, this is problematic during MRI. The pressure signal would be contaminated by the MR environment and the presence of the sensor would produce artifacts on the image.

In this paper we present a feasibility study for simultaneous pressure measurements using optical sensors during MRI to produce valid indices for the assessment of LV function. In vivo measurements on one sheep were performed and pressure-volume loops were derived by combining MRI-estimated ventricular volumes with simultaneous on-site pressure measurements for various inotropic states. PV loops were analyzed and functional parameters were computed to evaluate the response to inotropic agents.

Methods

Experimental Protocol

Animal preparation

Two female sheep were used in this study. An intervention on the first sheep was performed to compare the optical pressure sensor measurements against the Millar standard probes. A subsequent intervention on another sheep allowed simultaneous pressure and volume measurements. The ewes were placed in a sitting position to induce reflex akinesia then anesthetized with a slow injection of 1 g of diluted thiopental. Tracheal tube was inserted and forced ventilation was started with a concentration of 1.5-2% isoflurane to maintain anesthesia. To prevent clotting, acetylsalicylic acid (0.5 g) and heparin (3000 IU) were injected intravenously.

The carotid artery was catheterized and a vascular dilator was inserted under X-ray monitoring. The animals received humane care in compliance with the standards of the European Convention on Animal Care. The study was approved by a local institutional ethics committee (INRA, Paris, France). Qualified personnel supervised the procedures and adequate anesthesia was used to minimize unnecessary pain.

In vivo pressure measurements

An optical pressure sensor probe (model OPP-M, Opsens, Quebec, Canada) based on the white-light polarization interferometry technique was used. The probe tip, 0.4 mm in diameter and 0.5 mm in length, attached to a 10 m long optical fiber directly transduces pressure into an optical signal which is then sampled at 1kHz. This device is fully MRI compatible, and immune to radiofrequency effects. The sensor is linear, with a total error of 2 mmHg and a temperature sensitivity of 0.2 mmHg/°C. Furthermore, each probe has two specific constant calibration factors.The nude fiber was sheathed into a non-magnetic 4F catheter (C4F100D, Balt extrusion, Montmorency, France) 1 m in length. The sensor was positioned 5 mm before the catheter tip, facing a slot that had been manually cut in the catheter wall to ensure pressure transmission. The fiber was then glued into the catheter through another slot 5 mm above. A null point was set prior to inserting the sensor inside the sheep by zeroing the measurement at atmospheric pressure.

The validation of the optical probe against a reference pressure transducer (Millar, Houston, USA) was performed during an aortic intervention that served wider purposes than those stated in the present paper.

During the intervention, the mounted optical probe and the Millar sensor were inserted simultaneously through the dilator up to the aorta of a 56kg ewe. Pressure was recorded in vivo outside the magnet for about 5 minutes including 20 seconds during which ventilation was stopped to acquire a stable signal with no respiration modulation.

Simultaneous in vivo pressure and volume measurements

Pressure measurements with the optical probes were conducted simultaneously with MRI acquisitions in another female sheep (60 kg). After four acquisitions at a baseline condition (heart rate ranging between 78 and 83 bpm), the sheep was infused dobutamine at a rate of 5 µg/kg/min. Heart rate increased to 145-150 bpm and three acquisitions were performed. Dobutamine rate was then decreased to 2.5 µg/kg/min (heart rate decreased to 125-131 bpm) and a last acquisition was performed before waking the animal up. The whole MR procedure lasted 2 hours and 20 minutes.

MRI acquisitions

Cardiac imaging was performed with a 1.5 T MR scanner (Philips Medical Systems, Achieva) with a 5-element SENSE cardiac array coil and ECG triggering. After calibration and scouting sequences, true FISP gated cine images ("balanced TFE") were acquired in 12 sequential 8-mm short-axis slices (2 mm interslice gap) from the apex to the atrial-ventricular ring, for 30 cardiac phases (TE/TR= 1.7ms/3.4ms, flip angle of 60°, field of view=320 x 240 mm 2 , matrix=256 x 192, voxel size 1.25mm, echo train length of 10, readout bandwidth 1042 Hz). Mechanical ventilation was continued during the MRI acquisition procedures without respiratory synchronization. Other sequences were also acquired for other purposes, not mentioned here.

Offline Pressure-Volume loop calculation

For each slice, the LV was delineated semi-automatically and its area was computed then multiplied by the slice thickness to derive an elementary volume. These slice volumes were then added up to define the total LV volume. This was performed for all cardiac phases, and an LV volume evolution in time was obtained. During each of the imaging sessions the ventricular pressure was recorded continuously. Since we only have one volume curve reflecting a somewhat average volume evolution during the image acquisition, an average pressure cycle was derived for the loop computation. Figure 1 summarizes the cycle calculation.

To delimit a heart cycle, pressure peaks were considered to be the end points of a pressure period. By Local maxima were first identified at times t peak,i , then individual periods T i were computed as being the times separating two consecutive peaks, finally pressure cycles were extracted between t peak,i -Ti 2 and t peak,i + Ti 2 . A mean pressure curve was thus obtained by averaging all cycles.

applying a double threshold on the pressure signal, we delimited individual windows each containing a pressure peak. Afterward the local maximum was identified inside each window, thus defining a cycle end point. Once all peaks were determined, individual periods T i were computed, and cycles were extracted between t peak,i -Ti 2 and t peak,i + Ti 2 . These cycles were then averaged to yield a mean pressure curve.

The volume and pressure curves were then associated by assuming that maximum volume is attained with maximum dP dt Mason et al. [START_REF] Mason | Assessment of cardiac contractility: The relation between the rate of pressure rise and ventricular pressure during isovolumic systole[END_REF].

Pressure-Volume Analysis

Selected indices were estimated using the PV loop in order to evaluate systolic and diastolic functions as well as mechanical efficiency during the three inotropic states. The end-systolic elastance was evaluated as a contractility index and the myocardium compliance was estimated during diastole to characterize the filling process. The effective arterial elastance, representing the afterload, was also evaluated to assess ventriculoarterial coupling.

End-systolic elastance as contractility index

The end-systolic elastance E es index is given by the slope of the end-systolic pressure-volume relationship (ESPVR). The ESPVR is defined as the line connecting the end-systolic points of a series of PV loops with different ventricular filling volumes. It is approximated as linear in the physiological range and is characterized by its slope, and a volume intercept V 0 , a theoretical volume for which no pressure is developed Grossman et al. [START_REF] Grossman | Contractile state of the left ventricle in man as evaluated from end-systolic pressure-volume relations[END_REF], Mehmel et al. [START_REF] Mehmel | The linearity of the end-systolic pressure-volume relationship in man and its sensitivity for assessment of left ventricular function[END_REF], Suga et al. [START_REF] Suga | Load independence of the instantaneous pressure-volume ratio of the canine left ventricle and effects of epinephrine and heart rate on the ratio[END_REF].

To obtain the ESPVR several PV loops are usually measured while gradually decreasing preload, using a balloon occlusion of the vena cava, for example. (fifth order polynomial instead of a cosine), showed that a good estimation of intercept volumes can be achieved.

In this study we estimated the ESPVR following the method proposed by Brinke et al. Brinke et al. [START_REF] Brinke | Single-beat estimation of the left ventricular end-systolic pressure-volume relationship in patients with heart failure[END_REF].

We computed a maximum isovolumic pressure P iso by using a fifth order polynomial fitting scheme to fit the left ventricle pressure curve. We started after end diastole and excluded from the fitting the pressure data points that lie after dP/dt max and before dP/dt min and those after the point where dP/dt increased above 10% of dP/dt min .

LV compliance during the filling phase

The pressure-volume relationship during diastole, known as the end-diastolic pressure-volume relationship (EDPVR), is generally exponential and depends on the elastic properties of the myocardium, the wall thickness and the ventricle's equilibrium volume V 0 Glantz and Kernoff [START_REF] Glantz | Muscle stiffness determined from canine left ventricular pressure-volume curves[END_REF] which is the intercept of the EDPVR with the volume axis. At any given instant of the EDPVR, the ventricle's stiffness modulus E = ∂P ∂V is given by the slope of the tangent to the curve at that point. We can also think in terms of compliance and compute C = ∂V ∂P as the slope reciprocal. The EDPVR is independent of the contractility state. When the inotropic state is enhanced (or decreased), the diastolic filling portion of the PV loop remains on this curve while sliding towards the left (or the right) Dickstein et al. [START_REF] Dickstein | Assessment of right ventricular contractile state with the conductance catheter technique in the pig[END_REF], Thomas and Weyman [START_REF] Thomas | Echocardiographic doppler evaluation of left ventricular diastolic function[END_REF]. This of course leads to compliance alterations in response to inotropic changes. In this study the EDPVR was fitted by P = α e β(V -V0) -1 , and the compliance was computed using the slope of the EDPVR during the slow filling phase, where viscous forces might be ignored, thus reflecting purely elastic properties of the ventricle. This phase is defined as the portion of diastole lying between the point where pressure starts increasing slowly while the volume continues to increase, and the point where pressure starts increasing rapidly while volume remains somewhat constant.

Effective Arterial Elastance

The arterial elastance is determined by the pressure in the ventricle at end systole, and by the amount of blood that the ventricle ejected into the arterial system such as E a = ESP/SV , where ESP denotes the end-systolic pressure. It can be represented on the PV loop diagram as the line connecting the end systolic coordinates and the end-diastolic volume point on the volume axis.

Cardiac work efficiency

The efficiency of energy transfer from the ventricle to the arterial system, called cardiac work efficiency (CWE), is given as the ratio of the stroke work to the total mechanical energy, CW E = SW/P V A. The total mechanical energy of a ventricular beat known as the pressure volume area (PVA) is equal to the sum of potential energy and stroke work, P V A = P E + SW ,Suga [START_REF] Suga | Global cardiac function: mechano-energetico-informatics[END_REF]. The stroke work (SW) is given by the area inside the PV loop, and the potential energy (PE) is represented by the area of the triangle between the point of end systole on the PV loop, V 0 and the end-systolic volume (ESV) intercept on the EDPVR.

CWE is reciprocally related to Ea Ees Kameyama et al. [START_REF] Kameyama | Energy conversion efficiency in human left ventricle[END_REF], it increases with contractility enhancement, and decreases with afterload increase Nozawa et al. [START_REF] Nozawa | Efficiency of energy transfer from pressure-volume area to external mechanical work increases with contractile state and decreases with afterload in the left ventricle of the anesthetized closed-chest dog[END_REF]. offset was corrected such that the mean error would be zero. Figure 3 shows the extracted average cycles of both Opsens and Millar with no ventilation. Pressure measured with the optical transducer agrees well with the standard Millar sensor.

Results

In vivo pressure measurements validation

MRI acquisitions and PV loops calculations

As shown in figure 4 the images were not affected by the presence of the pressure sensors, however some artifacts due to breathing motion can be noted. Pressure-volume loops were computed for all measurements covering three contractility states, figure 5 shows an example of the computed pressure and volume curves for a baseline state. Overall we obtained 8 PV loops : 4 for baseline, 1 with dobutamine 2.5, and 3 with dobutamine 5.

Baseline measurements

The computed mean baseline parameters were compared to values obtained in other studies dealing To produce the PV loop, the curves were matched so that maximum volume would coincide with the maximum first derivative of pressure. et al. [START_REF] Segers | Predicting systolic and diastolic aortic blood pressure and stroke volume in the intact sheep[END_REF] in order to situate our calculations with respect to usual parameter values. These values are given in table 1 , shown as (mean value ± standard deviation)

Contractility states comparison

Figure 6 shows the experimental time-line and heart rate evolution throughout the consecutive states.

The instantaneous frequencies were computed, over the pressure recordings, as the inverse of the time interval separating two pressure peaks, as shown in figure 1a. Heart rate (HR) increases gradually after the dobutamine administration starts, it stabilizes then decreases when the dobutamine rate is diminished. An HR increase can be observed as we go from one rate to another, as the remaining dobutamine might have been injected rapidly when flushing the tubes for the next infusion. Of course, functional changes in response to the inotropic agents, go beyond the obvious HR variation. Not only did the pressure cycle period vary, but its shape was drastically altered. Figure 7

Volume measurements

Unlike the conductance catheter technique, PV loop calculation using MRI imposes the use of an average pressure cycle vs average volume instead of instantaneous pressure vs instantaneous volume. Thus, the PV loop here reflects a global performance cycle and does not take into account an exact punctual correspondence between individual volume and a pressure measurements.Moreover the time synchronization between the volume and pressure cycle in this study was done based on bibliographic data, a more accurate match could be achieved by using the pressure signal to trigger the image acquisition and thus synchronize both waveforms.

Due to the lack of real-time volume recordings, it was impossible to obtain a series of load varying PV loops using MRI because a steady state over several cycles is needed in order to compute a volume. For the same reason, the assessment of contractility transit states is not possible using MRI. This remains, for now, an advantage in favor of the conductance technique, however, future advances in MRI may enable the overcoming of this limitation.

Pressure measurements

Even though the dynamics of the pressure measurements performed here were shown to be accurate, these optical pressure sensors present an offset that drifts throughout the experiment, thus introducing an unknown pressure shift between measurements. This variable drift can exceed 2 mmHg/hour. To be able to compare the obtained loops and the extracted parameters, a correction was necessary to bring the pressure cycles to a common reference. To eliminate these relative offsets we relied on theoretical fact that the filling parts of all measured PV loops should lie on the same EDPVR. Hence, we computed that relation for one of the baseline loops (we arbitrarily chose the fourth) and then shifted the pressure cycles of the rest of the loops so that they would be supported by the estimated EDPVR. Of course this does not eliminate the "absolute" pressure offset that is added to the reference pressure measurement, but at least all loops had that same offset and could be compared. An attempt to estimate and correct the absolute offset was made by considering the EDPVR and ESPVR intersection point. Theoretically these relationships should have the same zero axis intercept. Thus, we could assume that the absolute offset would be eliminated by shifting the reference baseline loop so that EDPVR and ESPVR intersect at zero pressure. Of course the same shift would afterward be applied to the rest of the loops. This procedure would, however, not guarantee the elimination of the real offset. In fact it greatly depends on the ESPVR estimation, which itself is not very precise given that it is computed from a single beat. Besides given the instability of ESPVR and the equilibrium volume V 0 and their variability between repeated measurements Spratt et al. [START_REF] Spratt | The end-systolic pressure-volume relationship in conscious dogs[END_REF], the chosen loop probably does not reflect the real V 0 .

To sum up, the presented pressures here might not actually reflect the real ones, but present a residual offset that introduces errors in the calculations of some parameters, such as V 0 , P E, and E a . Other parameters such as E es , C , and SW , are independent of the pressure offset and thus can be assumed accurate even if there's an unknown shift of the PV loop. Regardless of whether the offset is substantial or not, comparison between the considered contractility states remains valid, since the relative offsets are corrected and the absolute shift applies equally to all states. Further studies will focus on modeling the sensors behavior to enable measurement correction.

Baseline PV loops indices

When comparing baseline measurements with other authors works we found that some parameters such as E es agree well with the literature range. Others, like the compliance and V 0 , differ significantly from what is obtained by other authors. These differences can be explained by the fact that the sheep studied by those authors are much smaller (weighting around 40 kg) than the one used in our experiments (60 kg).

Since the diastolic compliance decreases with the chamber size Diamond et al. [START_REF] Diamond | Diastolic pressure-volume relationship in the canine left ventricle[END_REF], Forrester et al. [START_REF] Forrester | Early increase in left ventricular complinace after myocardial infarction[END_REF],

smaller animals with smaller hearts have less compliant ventricles. Moreover, smaller volumes mean that the PV loops are shifted to the left, thus the ESPVR axis intercept is also shifted to the left, thus explaining the smaller V 0 .

Response to inotropic stimuli

Both end systolic elastance and diastolic compliance increased significantly in response to positive inotropic stimuli. Increasing the dobutamine concentration did not enhance the stroke work (the SV decreased with an insignificant increase in maximum pressure), it did however somewhat enhance the cardiac work efficiency. The SW varies very slightly with afterload but is very sensitive to preload Kass et al. [START_REF] Kass | Comparative influence of load versus inotropic states on indexes of ventricular contractility: experimental and theoretical analysis based on pressure-volume[END_REF] , and thus the study of the relationship between the EDV and the SW would be more suitable.

This relationship is called the preload recruitable stroke work (PRSW). It is a linear relationship Glower et al. [START_REF] Glower | Linearity of the frank-starling relationship in the intact heart : the concept of preload recruitable stroke work[END_REF] which slope reflects the inotropic state. It was also shown to be insensitive to afterload over the physiological range, and thus qualifies as a contractility index, although it might depend on cardiac geometry and heart rate Glower et al. [START_REF] Glower | Linearity of the frank-starling relationship in the intact heart : the concept of preload recruitable stroke work[END_REF]. Nonetheless, the PRSW is less sensitive to inotropic changes than E es and dP/dt max vs EDV Little et al. [START_REF] Little | Comparison of measures of left ventricular contractile performance derived from presure-volume loops in conscious dogs[END_REF]. Since preload conditions variations were not performed during this feasibility study, it was not possible to derive such relationships which rely on multiple PV loops.

A method to estimate the slope of PRSW from a single beat was proposed by Mohanraj et al. Mohanraj

and Feneley [START_REF] Mohanraj | Single-beat determination of preload recruitable stroke work relationship:derivation and evaluation in conscious dogs[END_REF], however their method did not exclude prior preload variation experiments. For their estimations they computed an empirical constant from a set of VCO experiments, which they used for the slope estimation of another group of the same species.

The ESPVR axis intercept V 0 varies slightly among states but not distinctly enough to be attributed to contractility changes. This volume should be independent of the inotropic state Suga et al. [38], Dickstein et al. [START_REF] Dickstein | Assessment of right ventricular contractile state with the conductance catheter technique in the pig[END_REF]. The changes observed here account for the instability of V 0 Spratt et al. [START_REF] Spratt | The end-systolic pressure-volume relationship in conscious dogs[END_REF] as well as measurement and estimation errors. Note that the V 0 range here is affected by the pressure offset and its correction.

The effective arterial elastance E a increases with inotropic enhancement as end systolic pressure (ESP) is increased. The ventriculoarterial coupling index E a /E es decreases when contractility increases. Studies have shown that Ea Ees = 0.5 ensures maximal CWE Starling [START_REF] Starling | Left ventricular-arterial coupling relations in the normal human heart[END_REF], Asanoi et al. [START_REF] Asanoi | Ventriculoarterial coupling in normal and failing heart in humans[END_REF]. On the other hand, given that the energy transfer from one elastic chamber to another is maximal when they both have the same elastance Sunagawa et al. [START_REF] Sunagawa | Optimal arterial resistance for the maximal stroke work studied in isolated canine left ventricle[END_REF], maximal SW would occur when E a = E es . It has been shown that the normal heart at rest operates at neither maximal CWE nor maximal SW, but in fact finds an optimal working point between the two, which is actually closer to maximal efficiency Starling [START_REF] Starling | Left ventricular-arterial coupling relations in the normal human heart[END_REF], Asanoi et al. [START_REF] Asanoi | Ventriculoarterial coupling in normal and failing heart in humans[END_REF].

In our baseline measurements, E a /E es is close to 1, indicating that the ventricle is working near maximal SW. With dobutamine, this index is approximately 0.6, thus the heart is functioning around the optimal working point between maximal SW and maximal CWE, while remaining closer to maximal efficiency.

Valid contractility index

The results obtained here showed that parameters computed using volume alone might not correctly reflect LV contractility. In fact, the stroke volume measurements confirmed that it is not a good contractility index. Despite the contractility enhancement, because of its load dependence, the stroke volume of the dobutamine loop is smaller than that of baseline. Since the preload (EDV) decreased, so has the SV. The same can be argued for the ejection fraction which did not reflect the inotropic increase induced by the dobutamine augmentation.

On the other hand, dP/dt max often used as a contractility index based on pressure measurements alone, behaved as expected in this study and did increase with inotropic stimuli. Even though this index is more sensitive to inotropic changes than E es , it depends on heart rate, and is still preload-dependent. Even when considering dP/dt max vs EDV to negate the preload effect, it was shown that E es remains a better index.

Kass et al. found that in the presence of both preload and afterload alterations, E es was relatively less affected than dP/dt max vs EDV . So, despite the E es lesser sensitivity to inotropic change, its minimal dependence on both types of load alterations, coupled with its adequate characterization of contractility, make it somewhat more advantageous than dP/dt max vs EDV Kass et al. [START_REF] Kass | Comparative influence of load versus inotropic states on indexes of ventricular contractility: experimental and theoretical analysis based on pressure-volume[END_REF].

The accuracy and load independence of E es is however questionable when obtained using single-beat estimation methods which rely on load-dependent elements such as dP/dt and end-diastolic volume. Sincepreload has not been changed during this study and only single beat methods were used, it would be more valid to view the estimated E es as a load-dependent approximation of the load-independent elastance Kjorstad et al. [START_REF] Kjorstad | Pressure-volume-based single-beat estimations cannot predict left ventricular contractility in vivo[END_REF]. Even though these methods might underestimate E es and do not agree closely with the vena cava occlusion (VCO) technique estimation, they are still reasonable approaches given their experimental advantages Brinke et al. [START_REF] Brinke | Single-beat estimation of the left ventricular end-systolic pressure-volume relationship in patients with heart failure[END_REF].

Conclusion

By analyzing the ventricular PV loop, several functional parameters can be extracted and used to assess the LV performance. In this study we showed the feasibility of PV loop computation using simultaneous MRI volume estimation and direct pressure measurements with optical sensors on a sheep in three inotropic states. The optical nature of the pressure sensors ensured the absence of pressure signal contamination as well as image artifacts. Functional parameters were extracted and compared to the values reported in the literature showing good agreement. Contractility indices such as end diastolic elastance, were possible to estimate despite the pressure sensors offset problem. Following up on this feasibility investigation, a more comprehensive study including several animals will be carried out.

List of abbreviations

CO

  Averaging over extracted periods

Figure 1 :

 1 Figure 1: Mean pressure cycle computationLocal maxima were first identified at times t peak,i , then individual periods T i were computed as being the times separating two consecutive peaks, finally pressure cycles were extracted between t peak,i -Ti 2 and t peak,i + Ti 2 . A mean pressure curve was thus obtained by averaging all cycles.

Figure 2

 2 Figure 2 shows two 20 second extracts of the recorded pressure signals in the aorta outside the MRI, one while the ventilation was off and another, almost 2 minutes later during normal ventilation. The Opsens

Figure 2 :

 2 Figure 2: Opsens and Millar pressure signals Extracts of pressure signals recorded simultaneously using OpSens and Millar are represented during a time period where the ventilation had been stopped, and later when the ventilation had been restarted. The OpSens offset was corrected so that the mean error would be null.

  with sheep Charan et al.[START_REF] Charan | Effect of increased coronary venous pressure on left ventricular function in sheep[END_REF], Lee et al.[START_REF] Lee | Ventricular restraint therapy for heart failure: The right ventricle is different from the left ventricle[END_REF], Bauer et al.[START_REF] Bauer | Left ventricular outflow tract mean systolic acceleration as a surrogate for the slope of the left ventricular end-systolic pressure-volume relationship[END_REF], Ratcliffe et al.[START_REF] Ratcliffe | Ventricular volume, chamber stiffness, and function after anteroapical aneurysm plication in the sheep[END_REF], Pilla et al.[START_REF] Pilla | Passive ventricular constraint to improve left ventricular function and mechanics in an ovine model of heart failure secondary to acute myocardial infarction[END_REF]

Figure 3 :

 3 Figure 3: Opsens and Millar average pressure cycles The mean pressure cycles of Opsens and Millar extracted during the period when the ventilation was stopped.

Figure 4 :

 4 Figure 4: Left ventricle short axis slices acquired with MRI at end systole MRI images used for LV volume estimation were not affected by the optical sensors but show some breathing motion artifacts.

Figure 5 :

 5 Figure 5: Associated LV pressure and volume curves Measured volume and average pressure cycle during baseline state.To produce the PV loop, the curves were matched so that maximum volume would coincide with the maximum first derivative of pressure.

Figure 6 :

 6 Figure 6: Heart rate evolution throughout the experiment Heart rate evolution vs the experimental time-line throughout the three consecutive inotropic states.

Figure 7 :

 7 Figure 7: Mean pressure cycles for the three inotropic states Extracted average cycles for baseline, dobutamine 2.5µg/kg/min and dobutamine 5µg/kg/min. The shape and duration of the pressure cycle is altered in response to positive inotropic stimuli.

Figure 8 :

 8 Figure8: End-systolic elastance for the three contractility states PV loops obtained during baseline, dobutamine 2.5µg/kg/min and dobutamine 5µg/kg/min. For each loop ESPVR is represented and the value of its slope E es is given in mmHg/ml .

  [START_REF] Suga | Load independence of the instantaneous pressure-volume ratio of the canine left ventricle and effects of epinephrine and heart rate on the ratio[END_REF], all ESPVR and EDPVR should meet at the same axis intercept. Positive inotropic interventions would increase E es and shift the ESPVR to the left, while depressors of the cardiac function would decrease E es shifting ESPVR to the right, nevertheless the volume axis intercept V 0 would not significantly change Spratt et al.

  The invasive nature of this procedure, limiting clinical applications, has motivated several authors to propose methods for computing the ESPVR from a single

PV loop. Takeuchi et al.

Takeuchi et al. [44]

simulated end-systolic isovolumic pressure by cosine fitting of the pressure curve, as done by Sunagawa et al.

Sunagawa et al. [43]

, to locate a theoretical maximum pressure P iso that would be attained if ejection did not take place. The E es is then computed as the slope of the line connecting the point defined by (EDV, P iso ) to the end-systolic point on the PV loop. Other authors relied on a time-varying elastance model to derive E es Senzaki et al.

[START_REF] Senzaki | Single-beat estimation of end-systolic pressure-volume relation in humans[END_REF]

, Shishido et al.

[START_REF] Shishido | Single-beat estimation of end-systolic elastance using bilinearly approximated time-varying elastance curve[END_REF]

, Chen et al.

[START_REF] Chen | Nonunvaasive single-beat determination of left ventricular end-systolic elastance in humans[END_REF]

. Recently Brinke et al. Brinke et al.

[START_REF] Brinke | Single-beat estimation of the left ventricular end-systolic pressure-volume relationship in patients with heart failure[END_REF]

using Takeuchi's technique with a modified fitting scheme

table 2 .

 2 Compared to baseline, E es increased significantly with positive inotropic stimuli. Figure

	shows mean pressure cycles, each representing

Table 1 :

 1 Comparison of PV loop extracted parameters with respect to the literature 8 shows an example of PV loops representing the three states (baseline measurement n°4, dobutamine Schmitt et al.[START_REF] Schmitt | Integrated assessment of diastolic and systolic ventricular function using diagnostic cardiac magnetic resonance catheterization: Validation in pigs and application in a clinical pilot study[END_REF] and Kuehne et al.Kuehne et al.[START_REF] Kuehne | Magnetic resonance imaging analysis of right ventricular pressure-volume loops: In vivo validation and clinical application in patients with pulmonary hypertension[END_REF] used liquid filled catheters to measure ventricular pressure down a 1m pressure line. The impedance of the transmission line reduces pressure amplitudes, and constitutes a low pass filter that attenuates high frequency components. A pressure sensor placed inside the ventricle would provide a more accurate representation of the actual LV pressure.

	The use of optical sensors in this study enabled pressure measurements directly in the ventricle without any
	interaction with the MRI acquisition.

2.5, dobutamine 5 measurement n°2) , as well their corresponding end-systolic elastances. Moreover, the cardiac work efficiency was enhanced by dobutamine administration. As expected, the diastolic compliance increased with inotropic enhancement. As contractility is increased, the PV loops shift to the left on the EDPVR, thus going towards lower slopes (i.e. stiffness).

4. Discussion

MRI is considered to be the gold standard for absolute cardiac volume measurements. Obtaining PV loops requires simultaneously recording pressure during the MRI sessions. Pattynama et al.Pattynama et al.

[START_REF] Pattynama | Magnetic-Resonance-Imaging Analysis Of Left-ventricular Pressure-Volume Relations -Validation With The Conductance Method At Rest And During Dobutamine Stress[END_REF]

used a Millar micromanometer to achieve such measurements, and argued that because it is made from brass the obtained pressure signals were not significantly altered, however, some image artifacts were observed.

Schmitt et al.

Table 2 :

 2 Functional indices extracted from the PV loops of the three inotropic states

Acknowledgment

The authors would like to thank Alain De Cesares with the Laboratoire d'Imagerie Fonctionnelle, U678, Paris, France, for his contribution to a part of data processing.

SV Stroke Volume

SW

Stroke Work

VCO Vena Cava Occlusion

Competing interests

The authors declare that they have no competing interest.

Declarations on funding and ethical approval

There was no funding for this research. This study did not involve any human subjects.

Authors contributions

DAR: participated in data acquisition, performed data processing, analysis and interpretation, and drafted the manuscript. ED: participated in the conception of the study, performed critical revision of the manuscript, and participated in data acquisition. LdR: participated in data acquisition. YB: participated in animal preparation and data acquisition. EM: participated in the conception of the study, and animal preparation. All authors read and approved the final manuscript


Charles-Alban Deledalle
email: charles-alban.deledalle@math.u-bordeaux.fr

Nicolas Papadakis
email: nicolas.papadakis@math.u-bordeaux.fr

Joseph Salmon
email: joseph.salmon@telecom-paristech.fr

On debiasing restoration algorithms: applications to total-variation and nonlocal-means

Bias in image restoration algorithms can hamper further analysis, typically when the intensities have a physical meaning of interest, e.g., in medical imaging. We propose to suppress a part of the bias the method bias -while leaving unchanged the other unavoidable part the model bias. Our debiasing technique can be used for any locally affine estimator including ℓ1 regularization, anisotropic total-variation and some nonlocal filters.

Introduction

Restoration of an image of interest from its single noisy degraded observation necessarily requires imposing some regularity or prior on the solution. Being often only crude approximations of the true underlying signal of interest, such techniques always introduce a bias towards the prior. However, in general, this is not the only source of bias. In many cases, even though the model was perfectly accurate, the method would remain biased. This part of the bias often emerges from technical reasons, e.g., when approaching an NP-hard problem by an easier one (typically, using the ℓ 1 convex relaxation of an ℓ 0 pseudo-norm).

It is well known that reducing bias is not always favorable in terms of mean square error because of the so-called bias-variance trade-off. It is important to highlight that a debiasing procedure is expected to re-inject part of the variance, therefore increasing the residual noise. Hence, the mean square error is not always expected to be improved by such techniques. Debiasing is nevertheless essential in applications where the image intensities have a physical sense and critical decisions are taken from their values. For instance, the authors of [START_REF] Denis De Senneville | Extended Kalman filtering for continuous volumetric MR-temperature imaging[END_REF] suggest using image restoration techniques to estimate a temperature map within a tumor tissue for real time automatic surgical intervention. In such applications, it is so crucial that the estimated temperature is not biased. A remaining residual noise is indeed favorable compared to an uncontrolled bias.

We introduce a debiasing technique that suppresses the extra bias -the method bias -emerging from the choice of the method and leave unchanged the bias that is due to the unavoidable choice of the model -the model bias. To that end, we rely on the notion of model subspace essential to carefully define different notions of bias. This leads to a mathematical definition of debiasing for any locally affine estimators that respect some mild assumptions.

Interestingly, our debiasing definition for the ℓ 1 synthesis (also known as LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] or Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]) recovers a well known debiasing scheme called refitting that goes back to the "Hybrid LASSO" [START_REF] Efron | Least angle regression[END_REF] (see [START_REF] Lederer | Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions[END_REF] for more details).

For the ℓ 1 analysis [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF], including the ℓ 1 synthesis but also the anisotropic total-variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], we show that debiasing can be performed with the same complexity as the primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] producing the biased estimate.

In other cases, e.g., for an affine version of the popular nonlocal-means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], we introduce an iterative scheme that requires only a few run of an algorithm of the same complexity as the original one producing the biased estimate.

Background

We consider observing f = f 0 + w ∈ R P a corrupted linear observation of an unknown signal u 0 ∈ R N such that f 0 = Φu 0 where Φ ∈ R N ×P is a linear operator and w is a random vector modeling the noise fluctuations. We assume that E[w] = 0 where E is the expectation operator. The linear operator Φ is a degrading operator typically with P N and with a non-empty kernel encoding some information loss such that the problem becomes ill-posed.

We focus on estimating the unknown signal u 0 . Due to the ill-posedness of the observation model, we consider variational approaches that attempt to recover u 0 from the single observation f as a solution of the optimization problem

u ⋆ f ∈ argmin u∈R N E(u, f ) . (1)
where E : R N ×R P → R is assumed to have at least one minimum. The objective E is typically chosen to promote some structure, e.g., smoothness, piece-wise constantness, sparsity, etc., that is captured by the so-called model subspace M * f . Providing u ⋆ f is uniquely defined and differentiable at f , we define M ⋆ f ⊆ R N as the tangent affine subspace at f of the mapping f → u ⋆ f , i.e.,

M ⋆ f = u ⋆ f +Im[J ⋆ f ] = u ∈ R N ; ∃z ∈ R P , u = u ⋆ f + J ⋆ f z with J ⋆ f = ∂u ⋆ f ∂f f (2)
where J ⋆ f is the Jacobian operator at f of the mapping f → u ⋆ f (see [START_REF] Vaiter | Model selection with low complexity priors[END_REF] for an alternative but related definition of model subspace). When u ⋆ f ∈ Im[J ⋆ f ], the model subspace restricts to the linear vector subspace

M ⋆ f = Im[J ⋆ f ].
In the rest of the paper, u ⋆ f is assumed to be differentiable at f 0 and for almost all f . Example 1 The least square estimator constrained to the affine subspace

C = b + Im[A], b ∈ R N and A ∈ R N ×Q , is a particular instance of (1) where E(u, f ) = ||Φu -f || 2 + ι C (u) (3)
and for any set C, ι C is its indicator function:

ι C (u) = 0 if u ∈ C, +∞ otherwise.
The solution of minimum Euclidean norm is unique and given by

u ⋆ f = b + A(ΦA) + (f -Φb) (4)
where for a matrix M , M + is its Moore-Penrose pseudo-inverse. The affine constrained least square restricts the solution u ⋆ f to the affine model subspace

M ⋆ f = b + Im[A(ΦA) t ] (as Im[M + ] = Im[M t ]). Taking C = R N with for instance Q = N , A = Id and b = 0, leads to an unconstrained solution u ⋆ f = Φ + f whose model subspace is M ⋆ f = Im[Φ t
] reducing to R N when Φ has full column rank. Example 2 The Tikhonov regularization (or Ridge regression) [START_REF] Tikhonov | On the stability of inverse problems[END_REF][START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF] is another instance of (1) where, for some parameter λ > 0 and matrix

Γ ∈ R L×N , E(u, f ) = 1 2 ||Φu -f || 2 + λ 2 ||Γ u|| 2 . ( 5
)
Provided Ker Φ∩Ker Γ = {0}, u ⋆ f is uniquely defined as u ⋆ f = (Φ t Φ+λΓ t Γ ) -1 Φ t f which has a linear model subspace given by M ⋆ f = Im[Φ t ]. Example 3
The hard thresholding [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], used when Φ = Id and f 0 is supposed to be sparse, is a solution of (1) where, for some parameter λ > 0,

E(u, f ) = 1 2 ||u -f || 2 + λ 2 2 ||u|| 0 , (6)
where ||u|| 0 = # {i ∈ [P ] ; u i = 0} counts the number of non-zero entries of u and [P ] = {1, . . . , P }. The hard thresholding operation writes

(u ⋆ f ) I f = f I f and (u ⋆ f ) I c f = 0 (7)
where

I f = {i ∈ [P ] ; |f i | > λ} is the support of u ⋆ f , I c f is the complement of I f on [P ],
and for any vector v, v I f is the sub-vector whose elements are indexed by I f . As u ⋆ f is piece-wise differentiable, its model subspace is only defined for almost all f as

M ⋆ f = {u ∈ R N ; u I c f = 0} = Im[Id I f ],
where for any matrix M , M I f is the sub-matrix whose columns are indexed by I f . Note that

Id I f ∈ R N ×#I f .
Example 4 The soft thresholding [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], used when Φ = Id and f 0 is supposed to be sparse, is another particular solution of (1) where

E(u, f ) = 1 2 ||u -f || 2 + λ||u|| 1 , (8)
with ||u|| 1 = i |u i | the ℓ 1 norm of u. The soft thresholding operation writes

(u ⋆ f ) I f = f I f -λ sign(f I f ) and (u ⋆ f ) I c f = 0 , (9)
where I f is defined as above, and, as for the hard thresholding:

M ⋆ f = Im[Id I f ].

Bias of reconstruction algorithms

Due to the ill-posedness of our observation model and without any assumptions on u 0 , one cannot ensure the noise variance to be reduced while keeping the solution u ⋆ f unbiased. Recall that the statistical bias is defined as the difference

Statistical bias = E[u ⋆ f ] -u 0 . ( 10
)
An estimator is said unbiased when its statistical bias vanishes. Unfortunately the statistical bias is difficult to manipulate when f → u ⋆ f is non linear. We therefore restrict to a definition of bias at f 0 = Φu 0 as the error u ⋆ f0 -u 0 . Note that when f → u ⋆ f is affine, both definitions match (the expectation being linear). Most methods are biased since, without assumptions, u 0 cannot be guaranteed to be in complete accordance with the model subspace, i.e., u 0 / ∈ M ⋆ f0 . It is then important to distinguish techniques that are only biased due to a problem of modeling to the ones that are biased due to the method. We then define the model bias and the method bias as the quantities

u ⋆ f0 -u 0 = u ⋆ f0 -Π M ⋆ f 0 (u 0 ) Method bias -Π (M ⋆ f 0 ) ⊥ (u 0 ) Model bias , (11)
where for any set S, Π S denotes the orthogonal projection on S and S ⊥ denotes its orthogonal set. We now define a methodically unbiased estimator as follows.

Definition 1 An estimator u ⋆ f is methodically unbiased if

∀u 0 ∈ R N , u ⋆ f0 = Π M ⋆ f 0 (u 0 )
We also define the weaker concept of weakly unbiased estimator as follows.

Definition 2 An estimator u ⋆ f is weakly unbiased if ∀u 0 ∈ M ⋆ f0 , u ⋆ f0 = u 0 . The quantity u ⋆ f0 -u 0 for u 0 ∈ M ⋆ f0 is called the weak bias of u ⋆ f at u 0 .
Remark that a methodically unbiased estimator is also weakly unbiased.

Examples. The unconstrained least-square estimator is methodically unbiased since

u ⋆ f0 = Φ + f 0 = Φ + Φu 0 = Π Im[Φ t ] (u 0 ) = Π M ⋆ f 0
(u 0 ). Moreover, being linear, it becomes statistically unbiased whenever Φ has full column rank since Φ + Φ = Id. However the constrained least-square estimator is only weakly unbiased: its methodical bias only vanishes when u 0 ∈ M ⋆ f0 , i.e., when there exists t 0 ∈ R Q such that u 0 = b+A(ΦA) t t 0 . The hard thresholding is also methodically unbiased remarking that u ⋆ f0 is the orthogonal projection on M ⋆ f0 = Im[Id I f 0 ]. Unlike the unconstrained least-square estimator, Tikhonov regularization has a non zero weak bias. The soft thresholding is also known to be biased [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF] and its weak bias is given by -λId I f 0 sign(f 0 ) I f 0 . Often, estimators are said to be unbiased when they are actually only weakly unbiased.

Definitions of debiasing

Given an estimate u ⋆ f of u 0 , we define a debiasing of u ⋆ f as follows. Definition 3 An estimator ũ⋆ f of u 0 is a weak debiasing of u ⋆ f if it is weakly unbiased and M⋆ f = M ⋆ f for almost all f , with M⋆ f the model subspace of ũ⋆ f at f . Moreover, it is a methodical debiasing if it is also methodically unbiased.

Examples. The unconstrained least square estimator is a methodical debiasing of the Tikhonov regularization, since it is a methodically unbiased estimator of u 0 and they share the same model subspace. The hard thresholding is a methodical debiasing of the soft thresholding, for the same reasons.

A good candidate for debiasing u ⋆ f is the constraint least squares on M ⋆ f :

ũ⋆ f = u ⋆ f + U ⋆ f (ΦU ⋆ f ) + (f -Φu ⋆ f ) ∈ argmin u∈M ⋆ f ||Φu -f || 2 (12)
where

U ⋆ f ∈ R N ×n with n = rank[J ⋆ f ] is a matrix whose columns form a basis of Im[J ⋆ f ]. Let V ⋆ f ∈ R n×P be a matrix such that J ⋆ f = U ⋆ f V ⋆ f .
The following theorem shows that under mild assumptions this choice corresponds to a debiasing of u ⋆ f .

Theorem 1 Assume that f → u ⋆ f is locally affine for almost all f and that Φ is invertible on M ⋆ f . Then ũ⋆ f defined in Eq. ( 12) is a weak debiasing of u ⋆ f .

Proof. Since f → u ⋆ f is locally affine, f → U ⋆ f can be chosen locally constant. Deriving [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF] for almost all f leads to the Jacobian J⋆ f of ũ⋆ f given by

J⋆ f = ∂ ũ⋆ f ∂f = ∂u ⋆ f ∂f + U ⋆ f (ΦU ⋆ f ) + ∂f ∂f -Φ ∂u ⋆ f ∂f = J ⋆ f + U ⋆ f (ΦU ⋆ f ) + (Id -ΦJ ⋆ f ) = U ⋆ f V ⋆ f + U ⋆ f (ΦU ⋆ f ) + (Id -ΦU ⋆ f V ⋆ f ) = U ⋆ f (ΦU ⋆ f ) + , (13)
since ΦU ⋆ f has full column rank due to the assumption that Φ is invertible on

M ⋆ f . It follows that M⋆ f = ũ⋆ f + Im[ J⋆ f ] = u ⋆ f + U ⋆ f (ΦU ⋆ f ) + (f -Φu ⋆ f ) + Im[U ⋆ f (ΦU ⋆ f ) + ] (14) = u ⋆ f + Im[U ⋆ f (ΦU ⋆ f ) + ] = u ⋆ f + Im[U ⋆ f ] = M ⋆ f , (15)
since ΦU ⋆ f has full column rank. Moreover, for any u 0 ∈ M ⋆ f0 , the equation Φu = f 0 has a unique solution u = u 0 in M ⋆ f0 since Φ is invertible on M ⋆ f0 . Hence, ũ⋆ f0 = u 0 is the unique solution of [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF], which concludes the proof.

The next proposition shows that the condition "Φ invertible on M ⋆ f " can be dropped when looking at u ⋆ f and ũ⋆ f through Φ. The debiasing becomes furthermore methodical.

Proposition 1 Assume f → u ⋆
f is locally affine for almost all f . Taking ũ⋆ f defined in Eq. ( 12), then the predictor

Φũ ⋆ f of f 0 = Φu 0 is equal to Π ΦM ⋆ f (f ) and is a methodical debiasing of Φu ⋆ f . Proof. Since ΦU ⋆ f (ΦU ⋆ f ) + = Π Im[ΦU ⋆ f ] , we have Φũ ⋆ f = Π ΦM ⋆ f (f ).
As the orthogonal projector on its own model space, it is methodically unbiased. Moreover

Im[Π Im[ΦU ⋆ f ] ] = Im[ΦU ⋆ f ],
hence Φũ ⋆ f and Φu ⋆ f share the same model subspace.

Remark 1 As an immediate consequence, the debiasing of any locally affine denoising algorithm is a methodical debiasing, since Φũ ⋆ f = ũ⋆ f .

We focus in the next sections on the debiasing of estimators without explicit expression for M ⋆ f , meaning that Eq. ( 12) cannot be used directly. We first introduce an algorithm for the case of ℓ 1 analysis relying on the computation of the directional derivative J ⋆ f f . We propose next a general approach, applied to an affine nonlocal estimator, that requires J ⋆ f δ for randomized directions δ.

5 Debiasing the ℓ 1 analysis minimization

From now on, the dependency of all quantities with respect to the observation f will be dropped for the sake of simplicity. Given a linear operator Γ ∈ R L×N , the ℓ 1 analysis minimization reads, for λ > 0, as

E(u, f ) = 1 2 ||Φu -f || 2 + λ||Γ u|| 1 . (16)
Provided Ker Φ ∩ Ker Γ={0}, there exists a solution given implicitly, see [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF], as

u ⋆ = U (ΦU ) + f -λU (U t Φ t ΦU ) -1 U t (Γ t ) I s I (17)
for almost all f and where

I = {i ; (Γ u ⋆ ) i = 0} ⊆ [L] = {1, .
. . , L} is called the co-support of the solution, s = sign(Γ u ⋆ ), U = U ⋆ f is a matrix whose columns form a basis of Ker[Id t

I cΓ ] and ΦU has full column rank. Note that s I and U are locally constant almost everywhere since the co-support is stable with respect to small perturbations [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF]. It then follows that the model subspace is implicitly defined as

M ⋆ = Im[U ] = Ker[Id t I cΓ ],
and so, the ℓ 1 analysis minimization suffers from a weak bias equal to -λU

(U t Φ t ΦU ) -1 U t (Γ t ) I s I . Given that u ⋆ ∈ Im[U ]
and it is locally affine, its weak debiased solution is defined for almost all f as ũ⋆ = U (ΦU ) + f .

(

) 18
The ℓ 1 synthesis [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF] consists in taking Γ = Id, hence U = Id I , so [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] becomes

u ⋆ I = (Φ I ) + f -λ((Φ I ) t Φ I ) -1 s I and u ⋆ (I) c = 0 . ( 19
)
Its model subspace is implicitly defined as M ⋆ = Im[Id I ], its weak bias is -λId I ((Φ I ) t Φ I ) -1 s I and its weak debiasing is ũ⋆ = Id I (Φ I ) + f . Subsequently, taking Φ = Id leads to the soft-thresholding presented earlier.

The anisotropic Total-Variation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is a particular instance of ( 16) where u 0 ∈ R N can be identified to a d-dimensional discrete signal, for which Γ ∈ R L×N , with L = dN , is the concatenation of the discrete gradient operators in each canonical directions. In this case I is the set of indexes where the solution has discontinuities (non-null gradients) and M ⋆ is the space of piece-wise constant signals sharing the same discontinuities as the solution. Its weak bias reveals a loss of contrast: a shift of intensity on each piece depending on its surrounding and the ratio between its perimeter and its area, as shown, e.g., in [START_REF] Strong | Edge-preserving and scale-dependent properties of total variation regularization[END_REF]. Note that the so-called staircasing effect of TV regularization is encoded in our framework as a model bias, and is therefore not reduced by our debiasing technique. Strategies devoted to the reduction of this effect have been studied in, e.g., [START_REF] Louchet | Total variation as a local filter[END_REF].

Since in general u ⋆ has no explicit solutions, it is usually estimated thanks to an iterative algorithm that can be expressed as a sequence u k converging to u ⋆ . The question we address is how to compute ũ⋆ in practice, i.e., to evaluate Eq. ( 18), or more precisely, how to jointly build a sequence ũk converging to ũ⋆ .

We propose a technique that relies on the observation that, given [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF], for almost all f , the Jacobian J ⋆ of u ⋆ at f applied to f , leads to Eq. ( 18), i.e.,

J ⋆ [f ] = U (ΦU ) + f = ũ⋆ (20)
since U and s I are locally constant [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF]. We so define a sequence ũk which is, up to a slight modification, the closed-form derivation of the primal-dual sequence u k of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Most importantly, we provide a proof of its convergence towards ũ⋆ . Note that other debiasing techniques could be employed for the ℓ 1 analysis, e.g., using iterative hard-thresholding [START_REF] Herrity | Sparse approximation via iterative thresholding[END_REF][START_REF] Blumensath | Iterative thresholding for sparse approximations[END_REF], refitting techniques [START_REF] Efron | Least angle regression[END_REF][START_REF] Lederer | Trust, but verify: benefits and pitfalls of least-squares refitting in high dimensions[END_REF], postrefinement techniques based an Bregman divergences and nonlinear inverse scale spaces [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF][START_REF] Burger | Nonlinear inverse scale space methods[END_REF][START_REF] Xu | Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising[END_REF] or with ideal spectral filtering in the analysis sense [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF].

Primal-dual algorithm

Before stating our main result, let us recall some of the properties of primaldual techniques. Dualizing the ℓ 1 analysis norm u → λ||Γ u|| 1 , the primal problem can be reformulated as the following saddle-point problem

z ⋆ = argmax z∈R L min u∈R N 1 2 Φu -f 2 + Γ u, z -ι B λ (z) (21)
where z ⋆ ∈ R L is the dual variable, and

B λ = {z ; ||z|| ∞ λ} is the ℓ ∞ ball.
First order primal-dual optimization. Taking στ < 1

Γ 2 2
, θ ∈ [0, 1] and initializing (for instance,) u 0 = v 0 = 0 ∈ R N , z 0 = 0 ∈ R L , the primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] applied to problem [START_REF] Tikhonov | On the stability of inverse problems[END_REF] reads

   z k+1 = Π B λ (z k + σΓ v k ), u k+1 = (Id + τ Φ t Φ) -1 u k + τ (Φ t f -Γ t (z k+1 )) , v k+1 = u k+1 + θ(u k+1 -u k ), (22)
where the projection of z over B λ is done component-wise as

Π B λ (z) i = z i if |z i | λ, λ sign(z i ) otherwise. ( 23
)
The primal-dual sequence u k converges to a solution u ⋆ of ( 16) [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. We assumed here that u ⋆ verifies ( 17) with ΦU full-column rank. This could be enforced as shown in [START_REF] Vaiter | Local behavior of sparse analysis regularization: Applications to risk estimation[END_REF], but it did not seem to be necessary in our experiments.

Debiasing algorithm

As pointed out earlier, the debiasing of u ⋆ consists in applying the Jacobian matrix J ⋆ at f to f itself. This idea leads to the proposed debiasing algorithm that constructs a sequence of debiased iterates from the original biased primaldual sequence with initialization ũ0 = ṽ0 = 0 ∈ R N , z0 = 0 ∈ R L as follows

   zk+1 = Π z k +σΓ v k (z k + σΓ ṽk ), ũk+1 = (Id + τ Φ t Φ) -1 ũk + τ (Φ t f -Γ t zk+1 ) , ṽk+1 = ũk+1 + θ(ũ k+1 -ũk ), ( 24
)
where

Π z k +σΓ v k (z i ) = zi if |z k + σΓ v k | i λ + β, 0 otherwise.
with β 0. Note that when β = 0, deriving z k , u k and v k for almost all f at f in the direction f using the chain rule leads to the sequences zk , ũk and ṽk respectively (see also [START_REF] Deledalle | Stein unbiased gradient estimator of the risk (SUGAR) for multiple parameter selection[END_REF]). However, as shown in Theorem 2, it is important to choose β > 0 to guarantee the convergence of the sequence1 .

Theorem 2 Let α > 0 be the minimum non zero value

2 of |Γ u ⋆ | i for all i ∈ [L].
Choose β such that ασ > β > 0. The sequence ũk defined in (24) converges to the debiasing ũ⋆ of u ⋆ .

Before turning to the proof of this theorem, let us introduce a first lemma.

Lemma 1

The debiasing ũ⋆ of u ⋆ is the solution of the saddle-point problem

min ũ∈R N max z∈R L Φũ -f 2 + Γ ũ, z -ι FI (z), (25)
where ι FI is the indicator function of the convex set

F I = p ∈ R L ; p I = 0 .
Proof. As ΦU has full column rank, the debiased solution is the unique solution of the constrained least square estimation problem ũ⋆ = U (ΦU ) [START_REF] Xu | Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising[END_REF] by the projection onto F I , i.e.,

+ f = argmin ũ∈U ⋆ Φũ -f 2 . ( 26
) Remark that ũ ∈ U ⋆ = Ker[Id t I cΓ ] ⇔ (Γ ũ) I c = 0 ⇔ ι F I c (Γ ũ) = 0, where F I c = p ∈ R L ; p I c = 0 . Using Fenchel transform, ι F I c (Γ ũ) = max z Γ ũ, z -ι * F I c (z), where ι * F I c is the convex conjugate of ι F I c . Observing that ι FI = ι * F I c concludes the proof. Given Lemma 1, replacing Π z k +σΓ v k in
Π FI (z) I c = zI c and Π FI (z) I = 0 , (27)
leads to the primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] applied to problem (25) which converges to the debiased estimator ũ⋆ . It remains to prove that the projection Π z k +σΓ v k defined in [START_REF] Xu | Iterative regularization and nonlinear inverse scale space applied to wavelet-based denoising[END_REF] converges to Π FI in finite time.

Proof (Theorem 2). First consider i ∈ I, i.e., |Γ u ⋆ | i > 0. By assumption on α, |Γ u ⋆ | i ≥ α > 0. Necessary z ⋆ i = λ sign(Γ u ⋆ ) i in order to maximize (21). Hence, |z ⋆ + σΓ u ⋆ | i ≥ λ + σα. Using the triangle inequality shows that λ + σα |z ⋆ + σΓ u ⋆ | i |z ⋆ -z k | i + σ|Γ u ⋆ -Γ v k | i + |z k + σΓ v k | i . (28)
Choose ε > 0 sufficiently small such that σα -ε(1 + σ) > β. From the convergence of the primal-dual algorithm of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], the sequence (z

k , u k , v k ) converges to (z ⋆ , u ⋆ , u ⋆ ). Therefore, for k large enough, |z ⋆ -z k | i < ε, |Γ u ⋆ -Γ v k | i < ε, and
|z k + σΓ v k | i λ + σα -ε(1 + σ) > λ + β . (29)
Next consider i ∈ I c , i.e., |Γ u ⋆ | i = 0, where by definition |z ⋆ | i λ. Using again the triangle inequality shows that

|z k + σΓ v k | i |z k -z ⋆ | i + σ|Γ v k -Γ u ⋆ | i + |z ⋆ | i . (30)
Choose ε > 0 sufficiently small such that ε(1

+ σ) < β. As (z k , u k , v k ) → (z ⋆ , u ⋆ , u ⋆ ), for k large enough, |z k -z ⋆ | i < ε, |Γ v k -Γ u ⋆ | i < ε, and
|z k + σΓ v k | i < λ + ε(1 + σ) λ + β . (31)
It follows that for k sufficiently large |z k +σΓ v k | i λ+β if and only if i ∈ I c , and hence Π z k +σKv k (z) = Π FI (z). As a result, all subsequent iterations of (24) will solve (25), and hence from Lemma 1 this concludes the proof of the theorem.

Debiasing other affine estimators

In most cases, U ⋆ cannot be computed in reasonable memory load and/or time, such that Eq. ( 12) cannot be used directly. However, the directional derivative, i.e., the application of J ⋆ to a direction δ, can in general be obtained with an algorithm of the same complexity as the one providing u ⋆ . If one can compute the directional derivatives for any direction, a general iterative algorithm for the computation of ũ⋆ can be derived as given in Algorithm 1.

The proposed technique relies on the fact that given n = dim(M ⋆ ) uniformly random directions δ 1 , . . . , δ n on the unit sphere of R P , J ⋆ δ 1 , . . . , J ⋆ δ n forms a basis of Im[J ⋆ ] almost surely. Given this basis, the debiased solution can so be retrieved from [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF]. Unfortunately, computing the image of the usually large number n of random directions can be computationally prohibitive.

The idea is to approach the debiased solution by retrieving only a low dimensional subspace of M ⋆ leading to a small approximation error. Our greedy heuristic is to chose random perturbations around the current residual (the strength of the perturbation being controlled by a parameter ε). As soon as ε > 0, the algorithm converges in n iterations as explained above. But, by focusing in directions guided by the current residual, the algorithm refines in priority the directions for which the current debiasing gets significantly away from the data f , i.e., directions that encodes potential remaining bias. Hence, the debiasing can be very effective even though a small number of such directions has been explored. We notice in our experiments that with a small value of ε, this strategy leads indeed to a satisfying debiasing, close to convergence, reached in a few iterations. The nonlocal-means example. The block-wise nonlocal-means proposed in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] can be rewritten as an instance of the minimization problem (1) with

E(u, f ) = 1 2 i,j w i,j ||P i u -P j f || 2 with w i,j = ϕ ||P i f -P j f || 2 2σ 2 (32)
Algorithm 1 General debiasing pseudo-algorithm for the computation of ũ⋆ .

Inputs: f ∈ R P , u ⋆ ∈ R N , δ ∈ R P → J ⋆ δ ∈ R N , ε > 0. Outputs: ũ⋆ ∈ R N and U ∈ R N×n ′ with n ′ n an orthonormal family of Im[J ⋆ ] Initialize U ← [ ] repeat until ũ⋆ reaches convergence Generate δ ← η/||η||, η ∼ NP (0, Id) (perturbation ensuring convergence) Compute u ′ ← J ⋆ (f -Φũ ⋆ + εδ) (perturbed image of the current residual) Compute e ← u ′ -U (U t u ′ ) (projection on the orthogonal of the current M ⋆ ) Update U ← [U e/||e||] Update ũ⋆ ← u ⋆ + U ((ΦU ) + (f -Φu ⋆ )) until where i ∈ [n 1 ] × [n 2 ] spans the whole image domain, j -i ∈ [-s, s] × [-s,
s] spans a limited search window domain and σ 2 is the noise variance. We denote by P i the linear operator extracting the patch at pixel i of size (2p + 1) × (2p + 1). Note that we assume periodical conditions such that all quantities remain inside the image domain. The kernel ϕ : R + → [0, 1] is a decreasing function which is typically a decay exponential function. Taking ϕ piece-wise constant3 , leads to computing u ⋆ and its Jacobian at f applied to δ for almost all f as follow u ⋆ i = j wi,j f j j wi,j and (J ⋆ δ) i = wi,j δ j j wi,j with wi,j =

k w i-k,j-k (33)
where k ∈ [-p, p] × [-p, p] spans the patch domain. Note that the values of w and w can be obtained by discrete convolutions leading to an algorithm with complexity in O(N s 2 ), independent of the half patch size p.

With such a choice of ϕ, the block-wise nonlocal filter becomes a piece-wise affine mapping of f and hence Algorithm 1 applies.

Numerical experiments and results

Figure 1 gives an illustration of TV used for denoising a 1D piece-wise constant signal in [0, 192] and damaged by additive white Gaussian noise (AWGN) with a standard deviation σ = 10. Even though TV has perfectly retrieved the support of ∇u 0 with one more extra jump, the intensities of some regions are biased. Our debiasing is as expected unbiased for every region.

Figure 2 gives an illustration of our debiasing of 2D anisotropic TV used for the restoration of an 8bits approximately piece-wise constant image damaged by AWGN with σ = 20. The observation operator Φ is a Gaussian convolution kernel of bandwidth 2px. TV introduced a significant loss of contrast, typically for the thin contours of the drawing, which are re-enhanced by our debiased result.

Figure 3 gives an illustration of our iterative debiasing for the block-wise nonlocal-means algorithm used in a denoising problem for an 8bits image enjoying many repetitive patterns and damaged by AWGN with σ = 20. Convergence has been considered as reached after 4 iterations only. Our debiasing provides favorable results with many enhanced details compared to the biased result.

Conclusion

We have introduced in this paper a mathematical definition of debiasing which has led to an effective debiasing technique that can remove the method bias that does not arise from the unavoidable choice of the model. This debiasing technique simply consists in applying a least-square estimation constrained to the model subspace chosen implicitly by the original biased algorithm. Numerical experiments have demonstrated the efficiency of our technique in retrieving the correct intensities while respecting the structure of the original model subspace. Our technique is nevertheless limited to locally affine estimators. Isotropic total variation, structured sparsity or nonlocal-means with smooth kernels are not yet handled by our debiasing technique, and left for future work.

Fig. 1 .

 1 Fig. 1. Solutions of 1D-TV and our debiasing on a piece-wise constant signal.

Fig. 2 .

 2 Fig. 2. (left) Blurry image f = Φu0 +w, (center) TV u ⋆ , (right) debiased ũ⋆ .

Fig. 3 .

 3 Fig. 3. (left) Noisy image f = u0 +w, (center) nonlocal-means u ⋆ , (right) debiased ũ⋆ .

In practice, β can be chosen as the smallest positive floating number.

If |Γ u ⋆ |i = 0 for all i ∈ [L], the result remains true for any α > 0.

For instance, by quantification on a subset of predefined values in [0, 1].


Patrick Amestoy

Alfredo Buttari

Guillaume Joslin

Excellent Jean-Yves L'

Mohamed Sid-Lakhdar

Clément Weisbecker

Michele Forzan

Cristian Pozza

Rémy Perrin

Valène Pellissier

Jean-Yves L'excellent

Remy Perrin

Shared memory parallelism and low-rank approximation techniques applied to direct solvers in FEM simulation

Keywords: Eddy currents, FEMs, linear systems, parallel algorithms, sparse matrices

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

I N 3-D finite element simulation of induction heating processes, the solution time is a limiting factor in the design and optimization of new devices. Time-harmonic electromagnetic (EM) problems coupled with thermal (TH) problems are solved in sequence, and the linear system solution in the EM problem is often the bottleneck. For this reason, solver performances have been verified only for the solution of the EM problem with constant material properties in a quasistationary time-harmonic regime. Although iterative methods are widely applied for the solution of this kind of problem, in this paper the performances of direct solvers are investigated. Therefore, an efficient direct method to solve large sparse complex matrices should exploit parallelization and reduce memory consumption. In this paper, a reduction of memory requirements through low-rank (LR) techniques and factorization times through shared-memory parallelism in multifrontal massively parallel sparse direct solver (MUMPS) will be discussed [START_REF]MUMPS[END_REF]. Memory usage and factorization times are compared with another popular sparse direct solver, parallel direct solver (PARDISO) [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF]. Matrices arise from the modelization of induction heating industrial devices. Heating of a susceptor by pancake coils and gear induction hardening is taken as test benchmarks [START_REF] Dughiero | Multi-crystalline silicon ingots growth with an innovative induction heating directional solidification furnace[END_REF]- [START_REF] Lupi | Modelling single-and double-frequency induction hardening of gear-wheels[END_REF]. Geometric model design, meshing, and matrix building are performed by a commercial software [START_REF]CEDRAT[END_REF]. Starting from the same geometry (Pancake or Gear), meshes are gradually refined to solve problems of different sizes leading to matrices ranging from 320 k to 2.5 M degrees of freedom.

II. DESCRIPTION OF THE FIRST BENCHMARK

To assess the performance of different sparse solvers for complex linear systems, the model of a pancake inductor designed to heat a graphite susceptor has been taken as the first benchmark problem. This system is used in a furnace for the production of solar grade silicon by the directional solidification system I-DSS [START_REF] Dughiero | Multi-crystalline silicon ingots growth with an innovative induction heating directional solidification furnace[END_REF]. In Fig. 1, the model geometry is presented and Fig. 2 shows its mesh. The numerical solution of the EM problem has been carried out by applying the well known A, AV formulation and using the complex representation of the sinusoidal quantities [START_REF] Dughiero | A translational coupled electromagnetic and thermal innovative model for induction welding of tubes[END_REF], [START_REF] Bíró | Edge element formulations of eddy current problems[END_REF]. The magnetic vector potential coupled with the scalar electric potential formulation has been chosen because it guarantees a precise solution of the eddy current distribution, mostly in models where there is no highly permeable conducting region [START_REF] Bíró | Edge element formulations of eddy current problems[END_REF], [START_REF] Mifune | A fast solver for FEM analyses using the parallelized algebraic multigrid method[END_REF]. Usually this formulation originates systems with more degrees of freedom than the so-called scalar one (used in the second benchmark). This model has been solved using four different types of discretizations, with a number of elements ranging from 600 k up to 2.9 M volume elements, leading to complex sparse systems from 320 k to 1.5 M unknowns. Typical eddy currents distribution is shown in Fig. 3.

III. DESCRIPTION OF THE SECOND BENCHMARK

Simulations of induction hardening process is one of the most challenging tasks for multi-physics simulation. A complete process description requires not only to compute the EM-TH coupled problem, but also to couple these results with other different physics, able to calculate the metallurgical transformation and residual stresses [START_REF] Lupi | Modelling single-and double-frequency induction hardening of gear-wheels[END_REF]. To obtain reliable results of the entire process, the solution of EM-TH coupled problem must be as accurate as possible. To properly describe the EM field distribution, a fine and possibly mapped mesh has to be built mostly in the heat affected region. In general, the model should consider the non-linear material properties of the steel, and consequently the solution requires to apply an iterative method. Furthermore, to properly describe the dependence upon temperature of material properties, the timedependent TH problem, part of coupled solution, has to be solved using several time steps.

Sometimes, all these requirements lead to a very long computation time. For the purpose of this paper, only the EM steady-state problem is solved with linear magnetic properties. The model used as benchmark describes a slice of the whole system and it is shown in Fig. 4. The 3-D eddy current problem is solved through T -T 0-φ formulation where the magnetic field within conducting regions is

H = T + T 0 -∇ (1)
where T (A/m) is the unknown electric vector potential and φ (A) is the total magnetic scalar potential. T 0 (A/m) represents the given imposed electric vector potential, which describes a prescribed current distribution, fed by an external source. Then, the induced current density can be derived from the Ampère's law [14]

J = ∇ × (T + T 0 ). (2)
The model mesh, shown in Fig. 5, contains 5 × 10 5 volume elements. With this discretization, the dimension of linear system is 370 k or 2.5 M when, respectively, first-and secondorder elements are used. Typical eddy currents distribution is shown in Fig. 6.

IV. TESTING ENVIRONMENT

The hardware platform for the evaluation is a dual Intel Xeon (X5670, six cores) with 96 Gb RAM. MUMPS version 4.10.0 [START_REF]MUMPS[END_REF] and another popular tool Intel MKL parallel direct solver (PARDISO) version 10.3.9 are compared, both in the in-core (IC) and out-of-core (OOC) mode. The effect of swapping temporary data on virtual memory is also evaluated. Configuration details are reported in Table I.

V. PRELIMINARY RESULTS

Although MUMPS is based on a message passing model for parallel computation, a single MPI process with multithreaded BLAS libraries is used for tests, as a better performance with such a configuration has been observed on multicore architectures. Both MUMPS and PARDISO are executed using 12 threads and the corresponding factorization times are shown in Fig. 7 and Table II. The IC factorization times are slightly different with MUMPS and PARDISO, leading to very close behavior on 12 cores (see Fig. 7). Memory usage for IC executions is also very similar (see Table II). In the OOC mode, MUMPS is slightly slower than PARDISO on small matrices, but uses less memory since PARDISO does not The missing values for PARDISO in OOC (in Table II and Fig. 8) are due to an error during the factorization step in OOC and with multithreads. Timings on largest problems (Fig. 8) show that MUMPS is less penalized by memory swapping than PARDISO.

VI. IMPROVEMENTS BY LR APPROXIMATION TECHNIQUES

A LR matrix can be represented in a form which decreases its memory requirements and the complexity of involved basic linear algebra operations, such as matrix-matrix products. This is formalized in [START_REF] Bebendorf | Hierarchical Matrices-A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF].

Let A be a matrix of size m × n. Let k ε be the approximated numerical rank of A at accuracy ε. A is said to be a LR matrix if there exist three matrices W of size m × k ε , Z of size n × k ε , and E of size m × n such that

A = W • Z T + E (3)
where E 2 ≤ ε and k ε (m + n) < mn. k ε is commonly referred to as the numerical rank at precision ε and can be computed, together with W and Z with a rank-revealing QR factorization. LR approximation techniques are based upon the idea to ignore E and simply represent A as the product of W and Z T , at accuracy ε.

In practice, matrices coming from applicative problems are not LR, which means that they cannot be directly approximated. However, LR approximations can be performed on sub-blocks defined by an appropriately chosen partitioning of matrix indexes [START_REF] Bebendorf | Hierarchical Matrices-A Means to Efficiently Solve Elliptic Boundary Value Problems[END_REF], [START_REF] Börm | Efficient Numerical Methods for Non-Local Operators[END_REF]. Theoretical studies based on mathematical properties of the underlying operators have shown that variable sets that are far away in the domain are likely to have weak interactions which translates into the fact that the corresponding matrix block has a LR. To benefit from this property, a format called block LR (BLR) is used [START_REF] Amestoy | Improving multifrontal methods by means of low-rank approximations techniques[END_REF]. This format can be exploited within internal data structures of the multifrontal method used in the MUMPS software to decrease the memory consumption and the operation count of the solver [START_REF] Amestoy | Grouping variables in frontal matrices to improve low-rank approximations in a multifrontal solver[END_REF].

In Table III, results show that the method is more efficient on large problems, which is a good property in the context of large scale computing (note that there is no compression for Pancake 2 with ε = 10 -10 ). This property has also been observed in many other applications. For Pancake 2, significant gains are obtained with a LR threshold of 10 -8 ; the memory footprint is reduced by almost a factor of three and the complexity is almost halved. In this case, a few steps of iterative refinement are performed to recover full precision from the original approximated precision of 10 -10 . For Pancake 4 with LR threshold precision of 10 -14 , a good compression is obtained naturally without any loss of accuracy. Note that by choosing a more aggressive threshold ε, the BLR format can be used to efficiently produce loosely approximated factors that can be used as effective preconditioners for iterative solvers.

VII. IMPROVEMENTS OF SHARED-MEMORY PARALLELISM IN MUMPS

To solve a sparse linear system, the multifrontal method transforms the initial sparse matrix into a (elimination) tree of much smaller dense matrices. This tree is traversed in a topological order and each node is computed following a partial LU decomposition. Then, a two-pass solve operation is applied on each node of the tree to find the solution of the original system.

The structure of the tree offers an inner parallelism, called tree parallelism, in which different processes work on data subsets, stored on different nodes of a network. This kind of parallelism is already exploited by MUMPS in distributed memory environments [START_REF] Amestoy | Grouping variables in frontal matrices to improve low-rank approximations in a multifrontal solver[END_REF]. However, in shared-memory environments only node parallelism is applied; many processes collaborate on the decomposition, working on the same node. In this approach, threaded BLAS libraries are preferred to parallelize dense matrix operations. To make a step further, the fork-join model of parallelism has been implemented in a new code, based on OpenMP programming. Our first goal was then to exploit tree parallelism in shared-memory environments, by adopting algorithms commonly used in distributed memory environments and by rearranging them to a multithreading model. Therefore, the so-called AlgL0 algorithm consists in finding a separating layer in the tree, called L0, such that node parallelism will still be applied above it, but tree parallelism will be applied under it through OpenMP.

The use of adequate memory mapping policies allows to improve the performance of MUMPS on symmetric multiprocessor and non-uniform memory access architectures. The local alloc policy, which consists in mapping the memory pages on the local memory of the processor that first touches them, is applied on data structures used under L0, in order to achieve a better data locality and cache exploitation. The interleave policy, which consists in allocating the memory pages on all memory banks in a round-robin fashion such that the allocated memory is spread over all the physical memory, has been used on data structures above L0 in order to improve the bandwidth.

Both MUMPS 4.10.0 and MUMPS 4.10.0 with new AlgL0 algorithm have been tested on the set of test matrices, with the stated memory allocation techniques. As reported in Table IV, this approach brings a remarkable reduction of computational time on all tested matrices. The benefit tends to decrease on large matrices because the fraction of workload in the top of the tree (above L0) gradually increases in comparison with the workload in the bottom of the tree. Furthermore, this gain should rise through the use of the interleave policy.

VIII. CONCLUSION

Significant improvements of parallel sparse direct solver MUMPS have been successfully tested in 3-D numerical simulation of induction heating industrial applications. Results show a remarkable reduction of both memory usage and number of performed operations. Furthermore, both tree and node parallelism are exploited to reduce the solution time.

Fig. 1 .

 1 Fig. 1. Geometry of the pancake coil.

Fig. 2 .

 2 Fig. 2. Mesh of the pancake coil.

Fig. 3 .

 3 Fig. 3. Eddy current distribution in the graphite susceptor. Values range from 1.2 to 120 kA/m 2 .

Fig. 4 .

 4 Fig. 4. Figure represents only a slice of the whole inductor workpiece geometry because of symmetries.

Fig. 5 .

 5 Fig. 5. Mesh of inductor, workpiece, and clampers.

Fig. 6 .

 6 Fig. 6. Eddy current distribution in the gear. Values range from 200 kA/m 2 to 300 MA/m 2 .

Fig. 8 .

 8 Fig. 8. IC (with swap) and OOC factorization times for large cases.

TABLE I TESTED

 I MATRICES Fig.7. IC factorization time (12 cores, small cases).

TABLE II MEMORY

 II

NEEDED FOR FACTORIZATION (GB, LARGE CASES)

TABLE IV TIME

 IV SAVE BY USING MUMPS 4.10.0 WITH AlgL0 ALGORITHM