Diane Rivet

Michel Campillo

Nikolai M Shapiro

Victor Cruz-Atienza

Mathilde Radiguet

Nathalie Cotte

Vladimir Kostoglodov

Seismic evidence of nonlinear crustal deformation during a large slow slip event in Mexico

 ) was much larger than the volumetric deformation (~10 -6 ) at the depths probed (~5-20 km). Moreover, the time dependence of the velocity perturbation indicated that it was related to the strain rate rather than the strain itself. This suggests that during strong slow-slip events, the deformation of the overlying crust shows significant nonlinear elastic behavior.

Slow slip events (SSEs) recently discovered in many subduction zones are now considered as a very significant component of the strain release process and the seismic cycle because they modify the loading and releasing on the plate interface and can affect the recurrence of large thrust earthquakes. Here we present a study of the deformation that the SSEs induce in the overriding continental crust in the Guerrero segment of the Middle-America subduction zone (MASZ). In contrast to the remaining part of the MASZ, this region has not experienced large subduction thrust earthquakes for more than 100 years and, therefore, is referred as the Guerrero seismic gap [START_REF] Suarez | Geometry of subduction and depth of the seismogenic zone in the Guerrero gap, Mexico[END_REF]. At the same time, several large SSEs have been detected in Guerrero after installation of continuous GPS receivers [START_REF] Kostoglodov | A large silent earthquake in the Guerrero seismic gap, Mexico[END_REF][START_REF] Lowry | Transient fault slip in Guerrero, southern Mexico[END_REF][START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF].

Those SSEs extended ≈150km along dip, rupturing a significant part of both the dipping portion and the 35 to 40km depth flat portion of the plate interface [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF][START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. While the SSEs have significant seismic moment (e.g. Mw ~7.5 for the 2006 event, according to [START_REF] Larson | The 2006 aseismic slow slip event in Guerrero, Mexico: New results from GPS[END_REF][START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF]), they do not radiate seismic waves that affect superficial soft soil layers by strong shaking [START_REF] Peng | Temporal changes of shallow seismic velocity around the Karadere-Duzce branch of the North anatolian fault and strong ground motion[END_REF][START_REF] Rubinstein | Evidence for widespread nonlinear strong ground motion in the M-W 6.9 Loma Prieta Earthquake[END_REF][START_REF] Sawazaki | Temporal change in site response caused by earthquake strong motion as revealed from coda spectral ratio measurement[END_REF]. Instead, they affect the Earth materials at depths only through slow stress changes. Therefore, analyzing the response of the Earth to SSEs might provide us with useful insights for understanding the deformation mechanisms within the deeper crustal layers.

We used continuous ambient seismic noise recorded by stations of the Meso-America Seismic Experiment (MASE), aligned along a North to South profile perpendicular to the coast (Figure 1), to recover the Green functions [START_REF] Shapiro | Emergence of broadband Rayleigh waves from correlations of the ambient seismic noise[END_REF] between pairs of stations at consecutive time intervals. The data recorded spans a 26-month period that included the 2006 Guerrero SSE (Figure 2). We computed cross-correlation functions (CCFs) between vertical components over 60-day windows that overlapped every 10 days, from January 2005 to July 2007.

We repeated this computation for different period bands. The time delay between a reference correlation function (RCF) stacked over the whole period and the 60-day CCFs can be related to the velocity changes within the propagating medium [Brenguier et al., 2008a;Brenguier et al., 2008b;[START_REF] Hadziioannou | Stability of monitoring weak changes in multiply scattering media with ambient noise correlation: Laboratory experiments[END_REF][START_REF] Wegler | Fault zone monitoring with passive image interferometry[END_REF]. Under the first-order assumption of a homogeneous perturbation in the crust, the relative time shift is related to the relative seismic velocity change by dv/v = -dt/t.

To ensure measurements that are independent of noise source variations, we took into account the travel-time changes only within the coda of the CCFs. The coda part of CCFs is made up of diffuse waves that scattered on heterogeneities of the crust and thus tend to lose the source signature. As a consequence travel time delay measured within the coda is less sensitive to source variations. For periods between 4 s and 10 s we also measured the travel-time changes of a CCF at a given time relative to the CCF at the previous time. We thus reduced the error related to the definition of a global reference cross-correlation function, which was generally the average cross-correlation in previous studies. For periods longer than 10 s, we also defined a seasonal reference by stacking the daily CCFs of the 2005 summer season (May to September, 2005). Applying these procedures allowed the seasonal bias to be limited, and thus revealed the robust features of the temporal changes in the seismic speed (see Supplementary Material for more details).

Global positioning system (GPS) records in Guerrero show the onset of the SSE in May 2006 [START_REF] Vergnolle | Slow slip events in Mexico revised from the processing of 11 year GPS observations[END_REF]. The MEZC GPS station, which is located approximately in the middle of the MASE profile, gave a good indication of when a perturbation in the crustal elastic properties was expected (Figure 2A). Figure 2B shows the average velocity changes along the entire profile, including all of the seismic stations, for periods ranging between 11 s and 15 s, which corresponded to the Rayleigh waves essentially sensitive to the middle-crust structure. We observed a rapid decrease in the seismic velocity starting in May 2006, with a maximum value of -0.2% of the initial velocity in June 2006; this was followed by a return to the initial velocity within 2 months.

We then investigated the geographical extension of the velocity perturbation by considering successive subsets of 11 neighboring stations from North to South. The result presented in Figure 2C shows that while the velocity perturbation related to the SSE affected all of the station subsets and had a North-to-South extension of at least 250 km, it was most pronounced in the northern part of the profile, rather far from the coast. We also observed the velocity change signature of the Atoyac earthquake (Mw 5.9) that occurred on April 14, 2007, a few kilometers off the Pacific coast [START_REF] Singh | A report on the Atoyac, Mexico, earthquake of 13[END_REF], which was seen as a decrease in the seismic velocity solely in the southern part of the MASE array (Figure 2C,D). These data demonstrate that noise-based measurements can discriminate localized speed variations and thus provide reliable information about the geographical extension of the velocity perturbation associated with the 2006 SSE. Our observations show that the SSE produced detectable changes in the elastic properties in the middle crust over a widespread area, which extended much farther inland than the slipping interface segment of the SSE (Figures 1,2C) [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF].

To constrain the extension of the velocity perturbation at depth, we quantified the seismic velocity changes in different frequency bands, from 5.0-6.9 s to 20-27 s (Figure 3B). It has been shown both theoretically and observationally that the seismic coda is mainly composed of surface waves [START_REF] Hennino | Observation of equipartition of seismic waves[END_REF][START_REF] Margerin | Energy partition of seismic coda waves in layered media: theory and application to Pinyon Flats Observatory[END_REF]. We therefore expect that the sensitivity of coda waves to velocity changes at depth depends on the period, i.e., shorter periods are sensitive to shallower structures, while longer periods sample the deeper crust. Velocity variations were measured between 1 s and 10 s during the Parkfield earthquake [Brenguier et al., 2008a] and the Wenchuan earthquake [START_REF] Chen | Distribution of seismic wave speed changes associated with the 12 May 2008 Mw 7.9 Wenchuan earthquake[END_REF], which indicated that they could be caused by perturbations in the shallow crustal layers due to strong co-seismic shaking. Unlike the velocity changes observed following regular earthquakes, we detected no measurable changes in velocity for periods shorter than 5 s during the 2006 SSE, which suggests that the SSE did not affect the superficial crustal layers. At the same time, we observed a strong perturbation of the velocity at periods between 7 s and 20 s, which indicates that this perturbation was caused by changes in the mechanical properties at depths ranging from 5 km to 20 km. The velocity change weakens at longer periods (i.e. for larger penetration of the waves). This suggests that the increasing pressure with depth reduces the sensitivity of the elastic waves to perturbation provoked by the SSE.

Our results show that recently developed noise-based passive seismic monitoring can detect perturbations of the elastic properties caused by relatively slow crustal deformation at depth. The period band in which we observed velocity changes associated with the SSE suggested a perturbation in the middle crust, rather than a localized change at the plate interface. Moreover, the perturbation of the velocity extended farther North than the SSE slipping zone at the interface (Figure 2). To better understand the nature of this observation, we computed the static strain field associated with the 2006 SSE. We used an elastic three-dimensional (3D) finite-difference code [START_REF] Olsen | ShakeOut-D: Ground motion estimates using an ensemble of large earthquakes on the southern San Andreas fault with spontaneous rupture propagation[END_REF] with the following model settings: a 2D velocity structure below the Guerrero province [START_REF] Iglesias | S wave velocity structure below central Mexico using high-resolution surface wave tomography[END_REF], the geometry of the subduction interface determined from receiverfunction analysis [START_REF] Perez-Campos | Horizontal subduction and truncation of the Cocos Plate beneath central Mexico[END_REF], and the final slip distribution of the 2006 SSE [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF]. Our numerical simulation shows that the SSE produced an extended increase in dilation in the middle crust and North of the slipping interface, with a maximum between the MEZC and IGUA GPS stations (Figure 2D). Both dilation and velocity changes affected the middle crust with similar geographical extensions, which suggested that the velocity change was related to a change in the dilation. The maximum relative change of dilation was, however, several orders of magnitude smaller than the relative velocity change observed, i.e. ~10 -6 versus ~10 -3 , respectively. The strain sensitivity of the seismic velocity change (dv/(v*d)) is thus estimated as ~10 3 , which corresponds to a stress sensitivity of the velocity change of 7 ×10 -3 MPa - 1 , assuming a bulk modulus of 5 ×10 4 MPa. This ratio remains within the range estimated from active seismic experiments in regions that have been affected by transient stress changes due to tidal loading, i.e. between 5 ×10 -3 and 2 MPa -1 [START_REF] Furumoto | Seismic wave velocity changes and stress build-up in the crust of the Kanto-Tokai region[END_REF][START_REF] Nishimura | Temporal changes in seismic velocity of the crust around Iwate volcano, Japan, as inferred from analyses of repeated active seismic experiment data from 1998 to 2003[END_REF][START_REF] Reasenberg | A precise, continuous measurement of seismic velocity for monitoring in situ stress[END_REF][START_REF] Yamamura | Long-term observation of in situ seismic velocity and attenuation[END_REF]. Another important characteristic of the observed crustal perturbation due to the SSE is that the seismic velocity change disappeared relatively quickly after reaching its maximum (i.e. ~2 months later), as compared to the crustal strain induced by the SSE. To better assess the temporal correlation between the strain and the seismic velocity changes, we simulated the 3D quasi-static time evolution of the elastic strain field in the crust using a SSE time-dependent slip model [START_REF] Radiguet | Spatial and temporal evolution of a long term slow slip event, the 2006 Guerrero Slow Slip Event[END_REF]. The results presented in Figure 3 show that the maximum velocity perturbation reached in June 2006 occurred during the phase of positive slip acceleration and was correlated with the maximum of the dilation rate.

The changes in the elastic modulus (i.e. seismic velocity changes) are too large to be explained by purely linear elastic behavior, as suggested by the large discrepancy between the values of the relative velocity changes and the crustal dilation. By examining Figure 3, it can be seen that the transient velocity perturbation did not occur at the maximum of the cumulative dilation expected from the elastic model, but instead at the time of the maximum strain change. These observations suggest that a non-elastic, or at least nonlinear elastic behavior, is reached at relatively low deformation level. The onset of the velocity change, or of the non-linear regime, corresponds to dilation of the order of 10 -6 . This also implies that the linear elastic modeling of the crustal strain associated with the SSE gives us only a rough approximation and that accounting for non-linear effects could be important for accurate quantitative description of these phenomena.

While the rheology of rocks at depth is subject to debate, our data can be understood in the light of recent theoretical and experimental studies. In particular, the dynamics of the rearrangement of cracks and pores at scales smaller than wavelengths used in the measurements might have an important role. This mechanism explains the strong decrease of elastic modulus observed in laboratory experiments on crustal samples under dynamic loading. As with our in situ measurements, the nonlinear elastic effects appear at strains as small as 10 -6 [START_REF] Johnson | Slow dynamics and anomalous nonlinear fast dynamics in diverse solids[END_REF][START_REF] Johnson | Nonlinear dynamics, granular media and dynamic earthquake triggering[END_REF][START_REF] Tencate | Slow dynamics in the nonlinear elastic response of Berea sandstone[END_REF] with a strain sensitivity of velocity change of the order of 10 2 . Laboratory experiments also showed that the initial and fast shear modulus reduction (i.e. seismic velocity reduction) was followed by a slower modulus recovery to its initial value, the so-called 'slow dynamics'. At the same time this recovery in the experiment is relatively fast (thousands of seconds) compared with the time resolution of our measurements, i.e. 60 days. The velocity reduction related to the SSE evolved through time very closely to the strain rate, and there was no velocity relaxation (recovery) much longer than 60 days. The confining pressure used in laboratory experiments was several orders below that expected in the middle crust.

However, the nonlinear elastic effects described above may be strongly enhanced by low effective pressure. Evidence of fluid content in the crust of the Guerrero province due to metamorphic dehydration of the subducting oceanic crust indicates a possible low effective pressure due to highly pressurized fluids within the crust [START_REF] Jodicke | Fluid release from the subducted Cocos plate and partial melting of the crust deduced from magnetotelluric studies in 9 southern Mexico: Implications for the generation of volcanism and subduction dynamics[END_REF][START_REF] Song | Subducting Slab Ultra-Slow Velocity Layer Coincident with Silent Earthquakes in Southern Mexico[END_REF]. Thus, the strain sensitivity of the velocity changes and the duration of the velocity recovery associated with the SSE suggest to interpret our observations in the framework of nonlinear mesoscopic elasticity, and that seismological records provide valuable insights into the deformation process at depth.

Figures Captions

Figure 1. Map of the Guerrero region and the 2006 SSE final slip distribution on the plate interface, as obtained from the inversion of the GPS cumulative displacements (8). Red triangles, positions of the Meso-America Seismic Experiment (MASE) broad-band stations. Blue triangles, positions of the permanent GPS stations. Gray patches, rupture zones of major earthquakes that have occurred in the last century. Black dashed line, the Middle America Trench. Thick gray 13 arrow, convergence rate between the Cocos and North America plates [START_REF] Demets | Effect of recent revisions to the geomagnetic reversal time-scale on estimates of current plate motions[END_REF].

Yellow dashed line, position of the velocity change and dilation profiles in Figures 2B and3.

Figure 2 .

 2 Figure 2. Correlation in time and space between the observed velocity changes and the SSE. (A) MEZC GPS station daily displacement time series for the West-East (WE) and South-North (SN) components. Yellow box, period of the 2006 slow-slip event. (B) Seismic velocity changes measured from the cross-correlation functions of the vertical-component continuous seismic records of the 26 MASE seismic stations at a period range of 11-15 s. (C) Seismic velocity changes at a period range of 11-15 s, as measured from the records of a North-to-South moving subset of 11 stations along the MASE array. Seismic velocity changes are represented as a function of the coastperpendicular distance of the center of the subset. (D) Map of the MASE array and dilation at 15km depth calculated from the 2006 SSE slip model. Black box, top north-moving subset of 11 MASE stations that corresponds to the first raw in panel (C). Gray dashed line represents the coastperpendicular distance in km.

Figure 3 .

 3 Figure 3. Comparison between the observed velocity changes computed from the entire seismic array and the modeled dilation. (A) Seismic velocity changes (red curve) measured in the period range of 10-13.6 s. Dilation (black curve) and dilation rate (increment over 20 days) (blue curve) computed at 15-km depth where the surface waves at periods of 10-13.6 s are mostly sensitive. (B) Seismic velocity changes measured at different period ranges; from 5-6.9 s to 20-27 s. (C) Increment of dilation over 20 days, and (D) total dilation, computed below the MASE array and averaged between 85 km and 220 km from the coast. (C) and (D) were computed from the SSE slip-propagation model, taking into account the secular inter-seismic slip model that tended to reduce the dilation with time.

Acknowledgements

Seismological data used in this study were collected during the Meso-American Seismic Experiment (MASE, http://www.tectonics.caltech.edu/mase/). We thank Paul Johnson, Shri Krishna Singh, Robert Clayton, Denis Legrand, Bérénice Froment, Céline Hadziioannou, Jean-Pierre Vilotte and Dimitri Zigone for discussions. This work was supported by the Agence Nationale de la Recherche (France) under the contract RA0000CO69 (G-GAP) and by the European Research Council (Advanced grant Whisper).


Richard L Weaver

Céline Hadziioannou

Éric Larose

Michel Campillo

Eric Larose

On the precision of noise correlation interferometry

Keywords:

Long duration noisy-looking waveforms such as those obtained in randomly multiply scattering and reverberant media are complex; they resist direct interpretation. Nevertheless, such waveforms are sensitive to small changes in the source of the waves or in the medium in which they propagate. Monitoring such waveforms, whether obtained directly or obtained indirectly by noise correlation, is emerging as a technique for detecting changes in media. Interpretation of changes is in principle problematic; it is not always clear whether a change is due to sources or to the medium. Of particular interest is the detection of small changes in propagation speeds. An expression is derived here for the apparent, but illusory, waveform dilation due to a change of source. The expression permits changes in waveforms due to changes in wavespeed to be distinguished with high precision from changes due to other reasons. The theory is successfully compared with analysis of a laboratory ultrasonic data set and a seismic data set from Parkfield California.

INTRODUCTION

) and later called "Coda wave interferometry" (Snieder et al.2002) compares coda waveforms from multiply scattered waves obtained under different circumstances or on different dates and detects changes in a medium. A multiply scattered

In many such cases the change is due to a uniform change of temperature, and thus a uniform change in wave velocity. To detect such changes, [START_REF] Weaver | Temperature Dependence of Diffuse Field Phase[END_REF] constructed a dilation correlation coefficient between waveforms φ 1 and φ 2 .

X( ) = φ 1 (t)φ 2 (t(1 + ))dt φ 2 1 (t)dt φ 2 2 (t(1 + ))dt (1)
X takes on a value of unity at = 0 if the two waveforms are identical. It will reach a value of unity at some characteristic value of if the two waveforms differ only by some temporal dilation.

The estimated degree of dilation between two waveforms is taken to be the value of at which X is maximum. X reaches a maximum of less than unity if the waveforms differ by more than dilation alone. Therefore, the value of X at its maximum, if it is less than unity, may be interpreted as a measure of the distortion between the waveforms.

An alternative formulation is Poupinet's doublet method [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras fault, California[END_REF], which breaks φ 1 and φ 2 into a series of short time windows at several distinct times t, and determines the apparent shift δt between them by examining conventional cross spectrum. δt as a function of t, and in particular its slope δt/t reveals a change in the medium. Poupinet developed the doublet method in which seismic signals from repeated seismic events could be compared to infer changes in the earth [START_REF] Poupinet | Monitoring velocity variations in the crust using earthquake doublets: An application to the Calaveras fault, California[END_REF]. [START_REF] Song | Seismological evidence for differential rotation of the Earth's inner core[END_REF] and Zhang et al. (1985) used this to show that certain earth crossing rays were shifted and distorted compared to versions some years earlier, indicating a relative rotation between the earth and its core.

The extensive literature in recent years on correlations of diffuse acoustic noise has reported theory and measurements in support of the notion that such correlations are essentially equal to the acoustic response that one would have at one receiver were there a source at the other [START_REF] Lobkis | On the emergence of the Greens function in the correlations of a diffuse field[END_REF][START_REF] Weaver | Diffuse waves in open systems and the emergence of the Green's function[END_REF][START_REF] Derode | Recovering the Green's function from field-field correlations in an open scattering medium[END_REF][START_REF] Snieder | Extracting the Green's function from the correlation of coda waves: A derivation based on stationary phase[END_REF][START_REF] Roux | Ambient noise cross correlation in free space: Theoretical approach[END_REF][START_REF] Gouédard | Cross-correlation of random fields: mathematical approach and applications[END_REF][START_REF] Tsai | On establishing the accuracy of noise tomography travel-time measurements in a realistic medium[END_REF]). More technically, what is recovered is the Green's function as filtered into the frequency band of the noise, whitened and symmetrized in time. Two different kinds of records can be correlated.

Sometimes it is coda that is correlated [START_REF] Campillo | Long-Range Correlations in the Diffuse Seismic Coda[END_REF]. Precision of noise correlation interferometry 3 use of controlled sources, and without waiting for a seismic event. Tomographic maps of seismic velocity with unprecedented resolution have been obtained [START_REF] Shapiro | High resolution surface wave tomography from ambient seismic noise[END_REF][START_REF] Sabra | Surface wave tomography from microseisms in southern California[END_REF].

The technique has even been applied on the moon [START_REF] Larose | Lunar subsurface investigated from correlation of seismic noise[END_REF]. Very commonly, the noise which is correlated is incompletely equipartitioned, such that the resulting correlation waveforms do not precisely correspond to the Green's function. A difference between noise correlation and Green's function can also be observed when one has not averaged enough raw data; the correlation may not have yet converged. Theoretical and applied work is ongoing in attempts to understand and correct for systematic errors due to these effects [START_REF] Weaver | On the correlation of non-isotropically distributed ballistic scalar diffuse waves[END_REF][START_REF] Froment | Estimation of the effect of non-isotropically distributed energy on the apparent arrival time in correlations[END_REF] ). Nevertheless, Hadziioannou et al. (2009) demonstrated that it is not necessary to reconstruct the Green's function to use correlations for monitoring purposes.

These two approaches have been combined into what may be termed noise-correlation interferometry [START_REF] Sabra | Extracting coda arrivals from cross-correlations of long period scattered waves during the Mount St. Helens 2004 eruption[END_REF][START_REF] Sens-Schönfelder | Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia[END_REF]Brenguier, Shapiro et al. 2008;Brenguier, Campillo et al. 2008) in which correlations of seismic noise taken in different circumstances are compared. The correlations may have been obtained from different samples of ambient noise, perhaps on different dates, or from the codas of different events. The correlations are of course never identical; they are often very different. One reason for a difference is that the source of the noise may be different (yet if the correlation has converged to the local Green's function, a change of noise source ought to have little effect). Continuous seismic sources can move and strengthen or weaken as weather changes at sea. It may also be that the correlation has not fully converged (ie, insufficient averaging has been done). A third possibility is that the local mechanical or acoustic environment may have evolved, in particular, the local wave speed(s) may have changed. It is this possibility that is of particular interest, as changes in seismic velocities are associated with relaxations after major seismic events (Brenguier, Campillo et al. 2008). In some cases changes in seismic velocity can be used to predict volcanic eruptions (Brenguier, Shapiro et al. 2008). Therefore, it is of great interest to be able to discern whether a change is due to a change in local environment or to a change in the character of the noise. The latter possibility is of some interest; the former is of great interest.

Our purpose here is to evaluate the precision with which wave speed changes can be evaluated. To do this we consider the case in which the two waveforms φ 1 (t) and φ 2 (t) differ only by noise so that the actual relative dilation without noise, is zero. We then ask for the apparent (non-zero in general) value of at which the corresponding X in equation ( 1) achieves its maximum. The next section calculates the root mean square ε 2 of this apparent, and erroneous, relative dilation. The subsequent sections compare this prediction with experiment.

DILATION CORRELATION COEFFICIENT

Here we examine the apparent waveform-dilation between two nominally identical signals. Theoretically, one ought to infer a relative dilation of zero, however, noise can corrupt the inference. Key to the following analysis is an understanding that the signals being discussed are like coda, in that they are statistically stationary with durations long compared to an inverse bandwidth. We take the two waveforms to have an identical part ψ(t), and to differ by noise 2µχ(t). In the limit µ → 0, the waveforms become identical and have no relative dilation. If µ = 0, there will be an apparent, but actually meaningless, temporal dilation between them. We wish to estimate this erroneous apparent relative dilation, and to identify any signatures that could be used to alert to the possibility of error.

Note that the common part ψ of the signals need not be the local Green's functions.

We split the difference between these two waveforms φ 1 and φ 2 , and define two signals ψ and χ;

φ 1,2 = ψ(t) ± µχ(t) (2)
The waveform dilation-correlation coefficient (1) between them is

X( , µ) = φ 1 (t(1 + /2)) φ 2 (t(1 -/2)) dt φ 2 1 (t(1 + /2)) dt φ 2 2 (t(1 -/2)) dt = 1 -2 /4 [ψ(t(1 + /2)) + µχ(t(1 + /2))] [ψ(t(1 -/2)) -µχ(t(1 -/2))] dt [ (ψ 2 + µ 2 χ 2 ) dt] 2 -4µ 2 [ χψ dt] 2 (3) = 1 -2 /4 N ( , µ) D(µ)
with N and D defined as, respectively, the numerator and the denominator of X. The integrations are typically taken over a finite time-window with tapered edges. We make the approximation that the change of variable t(1 + ε/2) → t and t(1 -ε/2) → t in the denominator only leaves a prefactor

1 -ε 2 /4.
The value of at which X achieves its maximum is the practitioner's estimate of the dilation between the waveforms φ 1 and φ 2 . It occurs at such that ∂X/∂ = 0, or,

0 = 1 -2 /4 D(µ) ∂X( , µ) ∂ = - N ( , µ) 4 + (1 -2 /4) ∂N ( , µ) ∂ . (4)
If the phase shift due to time dilation is much less than one oscillation, which implies tω 1 for all times t and frequencies ω of interest, it suffices to expand N ( , µ) through only the second power of :

N ( , µ) = ψ(t) + t 2 ψ(t) + t 2 2 8 ψ(t) + µχ(t) + µ t 2 χ(t) + µt 2 2 8 χ(t) F o r P e e r R e v i e

w

Precision of noise correlation interferometry 5

× ψ(t) - t 2 ψ(t) + t 2 2 8 ψ(t) -µχ(t) + µ t 2 χ(t) - µt 2 2 8 χ(t) dt. (5)
On collecting terms in N ( , µ) that are linear and quadratic in obtains:

N ( , µ) ∼ ψ(t) 2 -µ 2 χ(t) 2 dt + t 2 ψ(t) + µ t 2 χ(t) [ψ(t) -µχ(t)] dt + - t 2 ψ(t) + µ t 2 χ(t) [ψ(t) + µχ(t)] dt + t 2 2 8 ψ(t) + t 2 µ 2 8 χ(t) [ψ(t) -µχ(t)] dt + t 2 2 8 ψ(t) - t 2 µ 2 8 χ(t) [ψ(t) + µχ(t)] dt + t 2 ψ(t) + µ t 2 χ(t) - t 2 ψ(t) + µ t 2 χ(t) dt (6) = ψ(t) 2 -µ 2 χ(t) 2 dt + tµ ( χ(t)ψ(t) -χ(t) ψ(t) dt + 2 4 t 2 ψ(t)ψ(t) -µ 2 χ(t)χ(t) dt (7)
-

2 4 t 2 ψ(t) 2 -µ 2 χ(t) 2 dt.
The first term in 2 may be integrated by parts.

N ( , µ) ∼ ψ(t) 2 -µ 2 χ(t) 2 dt + tµ χ(t)ψ(t) -χ(t) ψ(t) dt - 1 2 2 t 2 ψ(t) 2 -µ 2 χ(t) 2 dt + 1 4 2 r (8)
where the quantity r is:

r = ψ(t) 2 -µ 2 χ(t) 2 dt -t ψ(t) 2 -µ 2 χ(t) 2 -t 2 ψ(t) ψ(t) -µ 2 χ(t) χ(t) , (9)
whose expectation is zero and whose typical value is much less (by a factor of t 2 ω 2 ) than the other coefficient of 2 in Eq. ( 8). For this reason we henceforth neglect it.

So finally, neglecting high order of ε,

∂N ( , µ)/∂ ∼ µt χ(t)ψ(t) -χ(t) ψ(t) dt - t 2 ψ(t) 2 -µ 2 χ(t) 2 dt. (10)
Equation ( 4) is satisfied for:

= n/d, (11) with
n = µ t χ(t)ψ(t) -χ(t) ψ(t) dt, d = t 2 ( ψ(t) 2 -µ 2 χ(t) 2 ) dt + 1 4 ψ(t) 2 -µ 2 χ(t) 2 dt.
Equation ( 11) is an expression for the apparent dilation induced by the difference 2µχ between the original waveforms φ 1 and φ 2 . Given specific ψ and χ, one could evaluate it. It will be more useful, however, to obtain statistical estimates for the apparent dilation given assumptions about the envelopes and spectra of ψ and χ. The numerator n has expectation zero, as χ and ψ are statistically unrelated.

Thus within the stated limit ωt << 1, differences φ 2 -φ 1 do not manifest as an apparent dilation and the expected dilation is zero.

VARIANCE ESTIMATION AND STATISTICAL ERROR

Given n = 0, one then seeks estimates for the root-mean-square of equation ( 11) in order to judge typical fluctuations around the expected zero. These will be made based on assumptions that ψ and χ are stationary, noise-like and Gaussian, with similar spectra, having central frequency ω c . ψ and χ have the same duration, long compared to the inverse of ω c . Without loss of generality it is also assumed that they have the same amplitudes ψ 2 = χ 2 = 1. They are taken to extend from a start time t 1 to an end time t 2 . Under these assumptions the denominator of ( 11) is estimated as :

d ≈ (1 -µ 2 ) 1 3 ω 2 c (t 3 2 -t 3 1 ) + 1 4 (t 2 -t 1 ) ≈ (1 -µ 2 ) 1 3 ω 2 c (t 3 2 -t 3 1 ). ( 12
)
The square of the numerator of ( 11) is

n 2 ≈ µ 2 tt ψ(t) χ(t) -ψ(t)χ(t) ψ(t ) χ(t ) -ψ(t )χ(t ) dt dt . (13)
On changing variables: t + t = 2τ , t -t = ξ and dropping the cross terms as having expectation zero, (13) becomes:

n 2 ≈ µ 2 τ 2 - ξ 2 4 ψ τ + ξ 2 χ τ + ξ 2 ψ τ - ξ 2 χ τ - ξ 2 + ψ τ + ξ 2 χ τ + ξ 2 ψ τ - ξ 2 χ τ - ξ 2 dτ dξ. (14)
Auto-correlation functions may be defined as

ψ τ + ξ 2 ψ τ - ξ 2 = ψ 2 (τ ) R ψ (ξ) = R(ξ), (15)
such that with similar expressions for χ. Then the expectation of the square of the numerator of ( 11) is:

ψ τ + ξ 2 ψ τ - ξ 2 ≈ ω 2 c ψ 2 (τ ) R ψ (ξ) = ω 2 c R(ξ), (16)
n 2 ≈ 2µ 2 τ 2 - ξ 2 4 ω 2 c R 2 (ξ) dτ dξ ≈ 2 µ 2 ω 2 c τ 2 dτ R 2 (ξ)dξ . ( 17
)
The first integral is merely (t 3 2 -t 3 1 )/3. The second requires knowing something of the spectra of ψ and χ, so we take these to be Gaussian and identical:

∼ exp(-(ω -ω c ) 2 T 2 ) + exp(-(ω + ω c ) 2 T 2 ).
T may be identified by noting that the -10 dB points are at ω c plus or minus ln 10/T . In this case R is related to the inverse Fourier transform of the power spectrum: R(ξ) = cos(ω c ξ) exp(ξ 2 /4T 2 ). Then the second integral in ( 17) is identified as T π/2.

Application of equations ( 11), ( 12), ( 17) requires that we also estimate the µ. The quantity µ is related to the maximum of the waveform dilation-correlation coefficient

X(0, µ) = N (0, µ) D(µ) = ψ(t) 2 -µ 2 χ(t) 2 dt [ ψ 2 + µ 2 χ 2 dt] 2 -4µ 2 [ χ ψ dt] 2 (18)
As χ and ψ are statistically independent, one estimates the following relation between the maximum of the dilation correlation coefficient and the parameter µ:

X = 1 -µ 2 1 + µ 2 (19)
Finally, the root mean square of the practitioner's (erroneous) estimate for the relative dilation between φ 1 and φ 2 is

rms = n 2 1/2 d = √ 1 -X 2 2X 6 π 2 T ω 2 c (t 3 2 -t 3 1 ) (20)
We recall that T is the inverse of the frequency bandwidth, t 1 and t 2 are begin and end time of the processed time-window in the coda, respectively, and ω c is the central pulsation. This expression scales inversely with the duration of the correlation waveform in units of the period, and inversely with the square root of the duration in units of the inverse bandwidth. In practice Eq. ( 20) can be very small.

The quantity ω(t 2 -t 1 ) represents the available time where coda waves are significantly larger than the noise; the duration of the waveform is in units of the period. The quantity T is the amount of time for one bit of information to be delivered, and corresponds roughly to the time of the initial source [START_REF] Derode | Ultrasonic pulse compression with one-bit time reversal through multiple scattering[END_REF]). Thus Eq. ( 20) can be recognized as scaling inversely the available time in the coda (this time is related to coda-Q), and inversely with the square root of the amount of information.

It also may be recognized that small X corresponding to waveforms φ 1 and φ 2 that are very different, permits the practitioners erroneous estimate of dilation to be large. It may be that lengthening the considered time interval t 2 -t 1 would increase the precision, however it could also diminish X: in principle there are trade-offs.

Application of Eq. ( 20) is straightforward. A practitioner's estimate of the relative dilation between two waveforms φ 1 and φ 2 may be compared to Eq. ( 20). Values in excess of Eq. ( 20) are consistent with the inference that the observed dilation is real. Changes in waveform source or other character should not generate apparent dilations in excess of Eq. ( 20). Furthermore, in absence of any actual dilation, estimates of of the order Eq. ( 20) will nevertheless be generated in practice. Such should be regarded as un-meaningful.

COMPARISON WITH EXPERIMENT

The prediction Eq. ( 20) has been compared to waveform dilation measurements in a laboratory ultrasonic experiment (Hadziioannou et al. 2009). For practical reasons, we mimicked here ambient noise correlation with diffuse waves correlation. Several piezoelectric sensors and sources were applied to a multiply scattering air-bubble filled gel (see Fig. 1). Sources and receivers were placed on opposite sides, 64 mm apart. Multiple scattering was strong: received waveforms f sr (t) from sources s to receivers r were coda-like, with envelopes that resembled the solution of a diffusion equation (fig. 2).

The auto-correlation of each f sr (t) was windowed between lapse times of 12.5 to 50 µsec, to yield the waveforms which we call g sr (τ ) (see fig. 3). Details of the experimental set-up are described in Hadziioannou et al (2009). The details are, however, unimportant here, as the present theory applies to any pair of coda-like waveforms φ 1 and φ 2 . The typical g sr is stationary over this interval and has a power spectrum centered on 2.35 MHz with -10 dB points at 1.7 and 3.0 MHz.

The tables below are formed by maximizing the dilation correlation coefficient X between sums φ = s g sr over different sets of sources {s}. Note that the φ are not Green's functions G rr , as the fields f sr (t) used to compose them were not fully equipartitioned. The excellent impedance match between the gel and the receivers prevented the field to be reflected back to the medium. The noise field thus lacked any components traveling from receiver to receiver. All tests were conducted at fixed temperature, the actual relative temporal dilation is therefore zero. Also, to mimic signals acquired at different date, we averaged correlations over different set of sources to eventually compare them. The addition of different sources results in an additional noise χ(t) in the correlations. The goal of the test is to measure the dilation induced by the difference in waveform due to a different source distribution.

Autocorrelation waveforms, like that illustrated in figure 3, in the interval from 12.5 to 50 µsec appear stationary. Thus we take t 1 =12.5 µsec, t 2 = 50 µsec, ω c = 15 rad/µsec; and T = 0.56 µsec and conclude from (20),

rms = 4 × 10 -4 √ 1 -X 2 2X (21)
Tables 1 and2 show two case studies. In the first case, autocorrelations calculated from the signals at a receiver r, as produced by eleven distinct sources s, were summed over to generate the reference waveform φ 1 = s g sr . For each of three comparison waveforms φ 2 , the same sum was done, keeping 21). A maximum value of X and the at which that X is maximum, are constructed for each of seven receivers (the seven columns)

and the three choices for the set of sources described in the text (the three rows). The root mean square of those is compared with the predictions of theory. That X is of order 90% is consistent with one source in ten having changed.

the first ten sources unchanged. In order to deliberately change the waveform without dilation, the eleventh source is replaced with sources number twelve, thirteen and fourteen respectively. This was repeated for each of seven receivers. In each case we compare three waveforms φ 2 with the reference φ 1 and evaluate X( ). The table shows the maximum value of X( ), and the value of that did this, for each of the 21 cases. For each of the seven receivers we calculate the rms of these three . If the only changes were to the source of the noise field, and not the medium, one would expect no dilation, or = 0. Nevertheless, the differences in sources do generate apparent (feeble but noticeable) dilations ∼ . Theoretical and experimental rms(ε) are of the same order of magnitude. Theory, especially in light of the approximate modeling of the spectrum, may be said to have done a good job predicting the fluctuations.

In the second study (Tab. 2), four sources were held constant, and two were varied. Here the reference waveform was constructed from a sum over six sources s g sr ; each of the other three The dilation (×10 -3 ) as obtained by maximizing X for each of these cases Sources 1 to 4 & 5-6 -0.0800 0.0400 0.4600 -0.0200 0.1400 0.0400 0.1600 Sources 1 to 4 & 7-8 -0.0200 0.0400 0.0400 -0.3400 -0.0800 -0.1400 0.0800 Sources 1 to 4 & 9-10 -0.4600 0.1600 0.5600 -0.0400 0.0800 -0.0400 0.0400

Experimental root mean square dilation (×10 -3 ) all sets 0.27 0.098 0.419 0.198 0.104 0.087 0.106

Theoretical root mean square (×10 -3 ) from Eq. ( 21) all sets 0.19 0.14 0.17 0.12 0.14 0.12 0.11

Table 2. As in table 1, but for sources that differ by more: fewer sources are kept fixed (four) and more sources are changing (two). This results in as maller values of X and a larger value of error (uncertainty).

waveforms was constructed by replacing sources number five and six in that sum with two others.

Again theory may be said to have done a good job: the rms theoretical predictions accurately fit the actual experimental errors within 40%. This means that Eq. 20 properly predicts the order of magnitude of the error.

COMPARISON WITH SEISMIC DATA FROM PARKFIELD

We also analyze data from seismic measurements near Parkfield, California. Brenguier et al. (2008) showed that correlation waveforms obtained from ambient seismic noise over a period of five years from 2002 to 2007 changed in a manner consistent with a decrease of the seismic velocity after the earthquake of 2004 (Fig. 4). This decreased velocity then relaxed like log(t) after the earthquake.

While they used the doublet technique, we have re-analyzed their data using the dilation coefficient (see Eq. 1). For each of 78 receiver pairs, we compared the 1550-day average correlation waveform with the correlation waveforms constructed from each of 1546 overlapping 5-day segments. The whitening operation before correlation ensures that the spectrum of the correlations is constant. Note that direct arrivals are never processed. A representative correlation waveform is shown in figure 5. Each such waveform was windowed between -50 and -20 seconds, and again from 20 to 50 seconds (thus excluding direct Rayleigh arrivals and emphasizing the multiply scattered diffuse part of the signal for which the theory was developed). As in the previous section, the details of the measurements are available elsewhere (Brenguier, Campillo et al. 2008) but are unimportant for the present purposes. An X and an were deduced for each day. Power spectra were centered on 0.5 Hz, with -10dB shoulders at 0.1 and 0.9 Hz. These numbers permit the evaluation of (20):

rms = 2, 4 × 10 -3 √ 1 -X 2 2X ( 22
)
Figures 6 (top and bottom) show the mean (over the 78 receiver pairs) values of X and between each of the 1546 overlapping 5-day correlation waveforms, φ 1 , and the correlation waveform φ 2 as obtained by averaging over the entire 5 year period. Except for the two events on days 152 and 437, and the slow relaxation after the latter, the dilation appears constant, with daily random fluctuations of order 10 -4 . A correlation coefficient X of 0.8 predicts a rms fluctuation of 10 -3 (Eq. ( 22)). On averaging over 78 pairs, this prediction is reduced by a factor √ 78, to 1.1 × 10 -4 , consistent with the observed fluctuations in . In light of the approximations, in particular that of modeling the spectrum as Gaussian and the waveform as stationary, we count this as excellent agreement.

The discontinuities in at December 22, 2003 and September 28, 2004 are of particular interest.

The latter is coincident with the Parkfield earthquake. Jumps in dilation on those dates by ∼ 0.8×10 -3

were interpreted (Brenguier, Campillo et al. 2008) as decrease of local seismic wavespeed. But one might wish to entertain the hypothesis that these jumps are due to a change in the source of the noise.

To examine the question, we evaluated X and using correlation waveforms φ 1 as averaged over a 70-day period before each event as a reference and correlation waveforms φ 2 obtained over a series of 5 day spans after the events. The relative dilation across the events are the same as seen in figure 6 (top), of order 5 × 10 -4 . The values of X for these pairs of waveforms varied between 60% and 70%. According to equation ( 22) divided by √ 78, the value of X would have had to be below 33% if this large and apparent dilation were to be due to a random function with no actual dilation. The relative dilation between correlation waveforms before and after the event is therefore due to changes in seismic Green's function, and not to changes in the source of the waves.

SUMMARY

Waveforms constructed by noise correlation can be extraordinarily sensitive to changes in material properties. Such waveforms are in principle affected by both changes in noise sources and changes in the acoustic properties of the medium in which the waves propagate. It has been shown here that long-duration diffuse waveforms permit changes in the source of the noise to be distinguished with high precision from changes due to a temporal dilation.

An expression was derived for the rms of the apparent dilation measured on two waveforms, when there is no actual dilation between the two. This apparent dilation can be an effect of e.g. a change in noise sources. The rms value thus allows us to distinguish between an erroneous dilation measurement due to waveform change, and a physical wavespeed change in the medium.

We have tested the validity of the rms value using data from laboratory experiments, and we find that the theory predicts errors well.

  Coda waves consist in a long duration random looking signals that follows the main arrivals from a strong seismic source; coda waves are due to single and/or multiple scattering. More commonly the diffuse noise is due to ambient seismic waves from continuously acting sources such as human activity or ocean storms. Much recent literature reports constructions of the earth's seismic response between two seismograph stations, without the

  L.Weaver et al.

  Precision of noise correlation interferometry7

Figure 1 .

 1 Figure 1. Experimental setup with ultrasound. An air-gel mixturemimics a multiple scattering medium. Coda waves sensed by the receivers (see Fig. 2) are processed like ambient seismic noise : they are autocorrelated and compared from one date to another.

Figure 2 .

 2 Figure 2. A typical signal f sr (t) in the gel resembles noise, under an envelope which is a solution of a diffusion equation.

Figure 3 .

 3 Figure3. A typical autocorrelation g sr , of the signal from one of the sources to one of the receivers. The interval from 12.5 to 50 µsec was selected for dilation coefficient evaluation.

Figure 4 .

 4 Figure 4. Map of the seismic stations used in our study. They are part of the High Resolution Seismic Network operated by the Berkeley Seismological Laboratory.

Figure 5 .

 5 Figure 5. A typical daily correlation waveform from the Parkfield data set. Dilations were constructed by comparing waveforms like this as windowed from 50 to +50 seconds, with the direct signal from 20 to +20 omitted.

Figure 6 .

 6 Figure 6. Top: Dilation averaged over 78 receiver pairs, using a 5-day sliding window. The grey squares indicate the 70-day reference windows. The best fit varied weakly and stochastically over this period, with two notable jumps, after December 22, 2003 and after the Parkfield earthquake on September 28, 2004. The latter jump was followed by a slow recovery. Fluctuations have an rms strength of about 10 -4 . Bottom: Dilation coefficient X. the maximum value of the dilation coefficient X, averaged over 78 receiver pairs.

Table 1 .

 1 Comparison of best-fit waveform dilations with the predictions of equation (

						Precision of noise correlation interferometry	9
	X for seven receivers and three different choices for the set of sources
	Sources 1 to 11 & 12 0.9312 0.9169 0.8893 0.8458 0.8226 0.8852 0.8683
	Sources 1 to 11 & 13 0.9394 0.8833 0.9083 0.8464 0.8872 0.8631 0.8928
	Sources 1 to 11 & 14 0.9458 0.9009 0.8730 0.7942 0.8322 0.8396 0.7979
	The dilation (×10 -3 ) as obtained by maximizing X for each of these cases
	Sources 1 to 11 & 12	0.06	0.04	-0.16	0.08	-0.04	-0.10	0.18
	Sources 1 to 11 & 13	-0.04	-0.04	-0.08	0.06	-0.14	-0.08	0.10
	Sources 1 to 11 & 14	-0.16	0.08	-0.12	0.00	-0.14	-0.24	-0.12
	all sets	F Experimental root mean square dilation (×10 -3 ) o r 0.1013 0.0566 0.1244 0.0577 0.1166 0.1571 0.1376
		P e Theoretical root mean square (×10 -3 ) from Eq. (21)
	all sets		0.07	0.09	r e 0.09 0.12	0.11	0.10	0.11
						R
						e
						v i e
								w

Geophysical Journal International

ACKNOWLEDGMENTS

Most of this work was done while the first author was a visitor at the Laboratoire de Géophysique Interne et Tectonophysique (Université J. Fourier). This work was partially funded by the ANR JCJC08 SISDIF grant, and by the ERC Whisper Advanced grant. All the seismic data used in this study came from the Parkfield HRSN and were collected by the Berkeley Seismological Laboratory (BSL).


Gaël Varoquaux
email: gael.varoquaux@inria.fr

Alexandre Gramfort
email: alexandre.gramfort@inria.fr

Bertrand Thirion
email: bertrand.thirion@inria.fr

Small-sample brain mapping: sparse recovery on spatially correlated designs with randomization and clustering

Functional neuroimaging can measure the brain's response to an external stimulus. It is used to perform brain mapping: identifying from these observations the brain regions involved. This problem can be cast into a linear supervised learning task where the neuroimaging data are used as predictors for the stimulus. Brain mapping is then seen as a support recovery problem. On functional MRI (fMRI) data, this problem is particularly challenging as i) the number of samples is small due to limited acquisition time and ii) the variables are strongly correlated. We propose to overcome these difficulties using sparse regression models over new variables obtained by clustering of the original variables. The use of randomization techniques, e.g. bootstrap samples, and clustering of the variables improves the recovery properties of sparse methods. We demonstrate the benefit of our approach on an extensive simulation study as well as two fMRI datasets.

Introduction

Functional brain imaging, for instance using functional MRI, is nowadays central to human neuroscience research and the corresponding medical applications. Learning statistical links between the observed brain images and the corresponding subject's behavior can be formulated as a machine learning problem [START_REF] Mitchell | Learning to decode cognitive states from brain images[END_REF]. Indeed, prediction from neuroimaging data has lead to impressive results, such as brain reading, e.g. guessing which image a subject is However, the main goal of a functional neuroimaging study is most often not prediction per se, but brain mapping, that is identifying the brain regions involved in the cognitive processing of an external stimuli. With fMRI, the data at hand is the brain activity amplitude measured on a grid on voxels. It is classically modeled as a linear effect driven by the cognitive task performed by the subjects. Detecting regions of the brain active in a specific task can be formalized as identifying the non-zero coefficients of a linear model predicting the external stimuli from the neuroimaging data.

From a statistical standpoint, this estimation is very challenging because of the dimensionality of the problem: the behavior must be linked to the brain activity measured on more than 20 000 voxels, often 50 000, with a few hundred observations. For this reason, the problem is most-often tackled as a mass-univariate approach: a model is fitted separately on each voxel. Detection of active voxels then suffers of a multiple comparison problem; the detection power decreases linearly with the number of variables tested. Yet, only a few brain regions are active for a given task and the linear model is sparse. There is great potential interest to use sparse recovery techniques [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF], that can recover active voxels suffering only a loss in detection power sub-linear in the number of voxels. However, with brain mapping as with many other experimental fields, the design matrix is imposed by the problem, and due to the strong correlations across regressors, univariate approaches are often more effective than multivariate approaches [START_REF] Haury | The influence of feature selection methods on accuracy, stability and interpretability of molecular signatures[END_REF].

The main contribution of this paper is to propose an efficient sparse recovery procedure well-suited to the specificities of brain mapping situations: spatially-clustered weights in very high-dimensional correlated designs. Unlike previous work in fMRI [START_REF] Carroll | Prediction and interpretation of distributed neural activity with sparse models[END_REF], we focus on recovery and not predictive power. We also provide a detailed empirical study of sparse estimation in the case of large spatial configurations for the weight and correlated designs. The organization of the paper is the following: in section 2 we review results related to small sample sparse estimation, in section 3, we expose our contributed method, and in section 4 we report experimental results on spatially-correlated synthetic data and brain images.

Preliminaries: small-sample estimation of sparse linear models

The problem at hand is the recovery of sparse models in the presence of noise. Formally, given n observed brain images composed of p voxels, X ∈ R n×p and the corresponding behavioral variable y ∈ R n related to the stimulus, the data and the target variable are related by a linear model y = X β + e with e ∈ R n the measurement noise and β ∈ R p the coefficients of the model to recover. More specifically, we want to find the sparsity pattern of β,

S = supp(β) = {i ∈ [1, . . . p] s.t. β i = 0}.
For the application at hand, n will typically be a few hundreds of images while p can be thousands of voxels.

In general, the sample complexity -that is the number of samples requires to bound the estimation error-of a linear model with n observations and p variables is O p . However, under certain conditions, for sparse signals, i.e. linear models for which k ≪ p coefficients are non-zeros, only O k log(p -k) observations are required to identify these coefficients [START_REF] Wainwright | Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1-constrained quadratic programming[END_REF][START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF], and the corresponding estimation problem can be solved using tractable algorithms such as ℓ 1 -penalized square loss regression -the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]-or greedy approaches [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF].

Two types of conditions govern the success of the recovery of a k-sparse vector from noisy measurements, i.e. model consistency of the estimator: i) subsets of the columns of design matrix X larger than k should be well conditioned, as for instance implied by the Restricted Isometry Property (RIP) [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF]. In particular, on the signal subspace, the design matrix X S should be sufficiently well conditioned. ii) the regressors on the signal subspace X S should not be too correlated to regressors on the noise subspace X S , as formalized by Tropp's Exact Recovery Condition (ERC) [START_REF] Tropp | Greed is good: Algorithmic results for sparse approximation[END_REF] or for the ℓ 1 -penalized regression the irrepresentable condition [START_REF] Zhao | On model selection consistency of lasso[END_REF], or mutual incoherence [START_REF] Wainwright | Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1-constrained quadratic programming[END_REF]. In addition, the smallest non-zero coefficient of β must scale as the standard deviation of noise and the inverse of the minimum singular value of X S . The number of observations necessary and sufficient to almost surely recover the sparsity grows as n min = 2 θ k log(p -k) where θ depends on the various factors listed previously [START_REF] Wainwright | Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ1-constrained quadratic programming[END_REF]. In the specific case of X i.i.d Gaussian distributed with identity covariance, θ ≈ 1. Importantly, the Lasso estimator yields at most n nonzero coefficients, and for n below n min , sparse methods will often fail dramatically to recover a sparsity pattern close to the ground truth.

A common cause for failure of sparse recovery is multicolinearity: the presence of strong correlations in the design X. In practice, if several columns of X are strongly correlated, sparse estimators will often select arbitrarily one of them and not the others, leading for instance to a high rate of false negatives in the support estimation if they all correspond to the signal subspace. In the extreme case of exactly collinear regressors, the Lasso, which is a non-strictly convex optimization, does not admit a unique solution. For this reason, the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] adds a strongly convex term to the Lasso in the form of a ℓ 2 penalty. This term has the effect of grouping together correlated regressors. It opens the door to sparse recovery in relaxed conditions, in particular with ill-conditioned signal-subspace design matrices X S [START_REF] Jia | On model consistency of the elastic net when p ≫ n[END_REF]. In addition, while the sample complexity of the elastic net scales similarly to that of the Lasso, the former can select more than n coefficients. Another approach to ensuring strong convexity in the presence of correlated variables is based on mixed norms, imposing an ℓ 2 penalty on sub-groups of features known a priori to covary [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. When the group structure is unknown, [START_REF] Grave | Trace lasso: a trace norm regularization for correlated designs[END_REF] proposed the trace-lasso, adapting the penalization to the design matrix. The limitation of this approach is that in the small-sample limit, the design matrix can only define a small number of strictly convex directions. For image data, correlation across neighboring voxels can be modeled using overlapping groups [START_REF] Jacob | Group lasso with overlap and graph lasso[END_REF].

Another challenge of high-dimensional sparsity recovery from noisy observations is the difficulty to control false positives, i.e. inclusion in the estimated support Ŝ of variables not present in the true support S. For the Lasso, theoretical bounds on false detections are conditional on choosing the right value for the parameter λ, controlling the amount of ℓ 1 regularization. This optimum depends on the noise level and structure, and is generally considered challenging to set in the small-sample settings. For these reasons, [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF] and [START_REF] Meinshausen | Stability selection[END_REF] introduce resampled estimators based on the Lasso. As it is well known in mass-univariate analysis, such sam-pling of the posterior yields control of the probability of false selection with weaker constraints on the model, here λ. In addition, different variables in a strongly correlated group can be selected across the resampling estimates, and as a result these schemes can in theory recover correlated variables and more variables than observations. Another Lasso variant, the randomized Lasso [START_REF] Meinshausen | Stability selection[END_REF] introduces a random perturbation of the design X by rescaling the regressors. It achieves sparse recovery under a RIPlike condition, bounds on sparse eigenvalues of the design, without the need of an irrepresentable-like condition imposing a strong separation between signal and noise regressors. Provided λ is chosen large enough, the randomized lasso yields good control of the inclusion of noise variables even with a small number of observations.

While most of the theoretical results on sparse recovery have been established for the square loss, these techniques can carry over to other losses, for instance the logistic loss, which is strongly convex as the square loss, but is well suited to classification problems, rather than regression. [START_REF] Bach | Self-concordant analysis for logistic regression[END_REF] extends non-asymptotic leastsquare results to logistic regression by approximating it as a weighted least square and show that similar conditions on the design apply for sparse recovery with ℓ 1 penalization. Small sample learning rates for ℓ 1penalized logistic regression were established earlier by [START_REF] Ng | Feature selection, l1 vs. l2 regularization, and rotational invariance[END_REF] based on the non rotational invariance of the ℓ 1 ball. This general argument gives a simple necessary condition for sub-linear sample complexity of empirical risk minimizers: the corresponding optimization problem should not be rotational invariant. In addition the rotational asymmetries of the learning problem should correspond to a good representation of the data. Indeed, adapting the representation of a signal such that its decomposition is more sparse is key to the performance of sparse methods.

Randomization and clustering for sparse recovery in brain imaging

Our approach to sparse recovery on spatially correlated images can be summarized by i) clustering highly-correlated variables before sparse regression ii) randomizing the clustering and sparse estimation. Our work is based on the randomized lasso [START_REF] Meinshausen | Stability selection[END_REF], which we first briefly present.

The randomized lasso The randomized lasso consists in randomly perturbing the design matrix by taking only a fraction of the training samples and ran-domly scaling each variable, in our case each voxel. By repeating the later procedure and then counting how often each variable is selected across the repetitions, each variable can be assigned a score. Higher scores denote variables likely to belong to the true support.

Let i ∈ [1 . . . l] denote the repetition and β i be the corresponding estimated coefficients. The design matrix X i is formed by a random fraction π of the training data. Each column of X i is then randomly scaled using a Bernoulli distribution to 1 or to 1 -α with equal probability. The procedure is a subsampling of the data and a random perturbation of each column. The stability score of each voxel v is then the percentage of the repetitions for which the voxel has a non-zero weight as estimated by sparse regression, i.e. is used for the prediction:

a v = #{i s.t. v ∈ supp(β i )}/l ∈ [0, 1]
. The estimated support is then defined as {v s.t. a v ≥ τ }. In the following, we set α = 0.5 and π = 75% as suggested by [START_REF] Meinshausen | Stability selection[END_REF], and use l = 200.

The randomized lasso does not fully address the problem of large groups of correlated variables. Indeed, the selection score will undoubtedly tend to zero as the size of group increase, and be confounded by selection scores of uninformative variables due to chance.

Adding randomized clustering To improve the stability of the estimation we propose to work with clusters of correlated variables. The motivation is that clustering variables and replacing groups of variables by their mean value will reduce the number of variables and improve the properties of the resulting design X red : less variables and less correlation between them. As the variables are voxels defined over a 3 dimensional grid and as correlation between voxels is strongly local, we can further constrain the clustering procedure to cluster only neighboring voxels. Following [START_REF] Michel | A supervised clustering approach for fMRI-based inference of brain states[END_REF], we use spatially-constrained Ward hierarchical clustering to cluster the voxels in q spatially connected regions in which the signal is averaged. When merging two singleton clusters, the Ward criteria chooses the samples with largest cross-covariance. The benefit of this algorithm is thus that it is a bottom-up approach that captures well local correlations into spatial clusters. In addition, it gives a multi-scale representation: a tree of nested regions.

We run the clustering algorithm in the randomization loop of the randomized lasso: after sub-sampling the observations and rescaling the variables. As the hierarchical tree is estimated each time on a random fraction of the data, the tree is different for every randomization. Note that a similar randomization on trees by resampling is performed in the Random Forests algo-Algorithm 1 Randomized ward lasso Require: Input: datasets X ∈ R n×p , target variable y ∈ R n , number of clusters q, penalization λ, number of resampling l, scaling α ∈ [0, 1], sub-sampling fraction π ∈ [0, 1]

1: for i = 1 to l do 2: Sub-sample: X ← X J , ỹ ← y J where J ⊂ {1 . . . n}, card(J) = π n, X ∈ R π n×p , ỹ ∈ R π n 3: Randomize feature scaling: X ← X•w, where w ∈ R p , w i = 1 -αµ i , µ i ∼ Bern(0.5) 4:
Cluster features, and use the clusters mean: Xred ← ward( X), Xred ∈ R π n×q 5:

Estimate β i red ∈ R q from Xred and ỹ with sparse linear regression 6:

Label initial features with estimated coefficient on corresponding cluster:

β i ← ward -1 (β i red ), β i ∈ R p 7: end for 8: return stability scores a v = #{k s.t. v ∈ supp(β i )}/l, a v ∈ [0, 1]
rithm [START_REF] Breiman | Random forests[END_REF]. It is motivated by the highvariance of estimators based on greedy tree construction. Similarly, in consensus clustering [START_REF] Monti | Consensus clustering: a resampling-based method for class discovery and visualization of gene expression microarray data[END_REF] resampling is applied to the construction of clusters to account for the sensibility of the algorithms, that are most often based on greedy approaches to nonconvex problems. After variable clustering, the sparse estimator can be fitted on a q-dimensional dataset. A variable is marked as selected in repetition k if it belongs to a cluster with a non-zero weight. As the estimated supp(β i ) is in R q , we will write v ∈ supp(ward -1 (β i )), where ward -1 amounts to assigning to the initial features the coefficient estimated for the center of the cluster that they belong to. The full estimation procedure, that we denote by randomized ward lasso, is summarized in Algorithm 1.

The main benefit of the additional clustering step is to reduce the correlations in the design matrix, therefore improving the behavior of sparse methods. Indeed, as correlated features are merged, the conditioning of submatrices of X red improves, and thus the conditions for control of false detections with the randomized lasso are fulfilled for q small enough. The key to the success of our method is this tight control: with high probability, clusters selected during each iteration all contain an original non zero feature. When back-projecting the selection frequency in the non-clustered space, because of the randomness of the cluster shapes in the different iterations, the selection frequency score can outline structures of a smaller size than the cluster size. There is an inherent trade-off in the choice of q: large clusters will most likely not match the geometry of the support of β but are more easily selected by sparse regression and are therefore more stable. Note that the separation between the signal and noise subspace may not be improved by the clustering, but it is not central to the success of the randomized lasso.

Importantly, thanks to the clustering, our method can select more variables than the number of observations, which would be impossible with a standard Lasso-type estimator and difficult with only randomization. Formally, the randomized ward lasso can be seen as applying a sparse estimator on random projections, learned by clustering. Importantly, those random projections are made of a small number of voxels, and thus are consistent with the assumed sparse representation of our problem. In addition, as they are not performed in rotationally-invariant directions, the resulting estimator is not rotationally invariant, and we suggest that it can achieve sub-linear sample complexity.

Finally, the randomized ward lasso is computationally tractable in the large p settings. Indeed, agglomerative clustering can be implemented in linear time with regards to the number of possible agglomerations tested [START_REF] Müllner | Modern hierarchical, agglomerative clustering algorithms[END_REF]. As we constrain clustering on the spatial grid, this number of agglomeration grows only linearly with p. The cost of a Lasso estimate on an n × p problem is O n p k if it selects k non-zero coefficients. Thus, for l resampling iterations q = κ p clusters, κ ∈ [0, 1], the computational complexity of the randomized ward lasso scales as O l κ 2 n p k . Importantly for our applications, the scaling in p is linear. This is to be contrasted with sparse methods imposing a spatial regularization, such as group Lasso with overlapping spatial groups.

Experimental results

Synthetic data

Data generation We compare different sparse estimators on simulated data for which the ground truth is known. We simulate data according to a simple model reflecting the key elements of brain imaging data: local correlations of the design matrix X, i.e. the observed brain images, and spatially clustered non-zero weights. Specifically, we work on a (32 × 64) grid of 2048 features. We generate an i.i.d Gaussian design X that we spatially smooth with a 2D Gaussian filter of standard deviation σ. The weight vectors β to recover have k = 64 non-zero coefficients. We split them into spatial clusters of size c, evenly-separated on the spatial grid. Non-zero coefficients are uniformly sampled between β min = .2 and 1 + β min . Finally, the target variable y is generated according to the linear model with additive Gaussian i.i.d noise. The amplitude of the noise is set to ensure that X explains 80% of the variance of y under the true model.

The two parameters of our procedure, the smoothing and the cluster size, control to which point the synthetic data violate the conditions for recovery. Sufficient smoothing can render arbitrarily ill conditioned groups of regressors located in a spatial neighborhood, thus violating sparse eigenvalue properties or the RIP. In particular, if the clusters are large, the design matrix restricted to the signal subspace X S is ill conditioned, whereas if they are small, the non-zero coefficients are well separated and X S is well conditioned. Finally, as the smoothing increases, so does the coupling between the signal and noise subspaces, although this coupling is less important for large clusters.

We investigate sample sizes n = 128 and 256, which are standard numbers for fMRI datasets. Note that with no smoothing, elements of X are Gaussian i.i.d distributed, with identity covariance, and n min ≈ 1000.

Success metrics

We compare the ability to recover the true support of different estimators. Each estimator yields a score for each variable. The higher the score, the more likely is the variable to be in the support. We use precision-recall curves to quantify the ability of each estimator to recover the non-zero coefficients as the threshold on the scores is varied. We summarize the precision-recall by the area under the curve (AUC). An AUC of .5 is chance, while 1 is perfect recovery: there exists a threshold such that active and non-active features are perfectly separated. In practice, we consider that an AUC above 0.95 is a near-perfect recovery, and an AUC above .75 as usable. To match the context of real data where ground truth is unknown, note that we report performance for automatically selected parameters, and not best possible performance, unlike e.g. [START_REF] Grave | Trace lasso: a trace norm regularization for correlated designs[END_REF].

Results

We investigated two non-convex approaches: a greedy algorithm, orthogonal matching pursuit [START_REF] Mallat | Matching pursuits with timefrequency dictionaries[END_REF], and an iterative approach, the adaptive Lasso [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF]. We do not report the results as they did not achieve useful recovery. This failure can be understood because of the low number of samples, below n min in non smoothed design, and the fragility of non-convex estimators when working with correlated designs. In addition, we studied screening based on F-tests, the elastic net, the randomized lasso, and our contribution, the randomized ward lasso. Each time, we set the regularization parameters using 6-fold cross-validation to maximize the explained variance of left out data.

For elastic net, we also set by cross-validation the amount of ℓ 2 regularization, and for randomized ward lasso, the number of clusters. Setting penalties by minimizing test error gives no guaranties on estimation error in terms of recovery, but alternative methods like information-based criterion (IC) such as BIC or AIC do not hold in theory for small-sample settings. In our simulations cross-validation was superior over IC.

Results can be seen on figure 1. With no smoothing, all methods fail to give satisfactory recovery, although very large clusters can be partially recovered by randomized ward lasso with n = 256. For large cluster sizes, adding a small amount of correlation to the design improves the recovery. Indeed, the added correlation is consistent with the sparsity pattern to recover. Elastic net and randomized lasso can then recover even small clusters. Due to the small sample size, no method performs well in case of small clusters without any smoothing, nor in the case of small clusters and large smoothing, which corresponds to an ill-posed deconvolution problem. In the large cluster and large smoothing case, the univariate approach outperforms all others. In this situation, sub-matrices of the design matrix are too ill-conditioned for sparse methods to succeed, but the effective number of degrees of freedom to be estimated is strongly reduced by the inter-feature correlations: the univariate approach does not perform many independent hypothesis tests and thus the number of false positive is reduced. In a large region of the parameter space explored, for clusters composed of 8 or more features, and smoothing standard deviation of the order of a pixel, our proposed method, the randomized ward lasso, outperforms other methods. We note that in the small sample case (n = 128) parameter selection is particularly difficult for this method. Interestingly, we find that, in the settings that we explore, the randomized lasso is most often outperformed by the elastic net, suggesting that the elastic net's relaxation of the conditions on sparse eigenvalues is better suited to the problem than the randomized lasso's relaxation of the irrepresentable condition. Finally, we confirmed empirically that, in the randomized ward lasso, randomizing the clusters at each iteration was important: fully-random clusters outperform using a non-random clusters, and our randomized approach performs best (Fig. 2).

FMRI data

To assess the performance of our approach on real fMRI data, we investigated two datasets. The first one from [START_REF] Tom | The neural basis of loss aversion in decision-making under risk[END_REF] is a gambling task where the subject is asked to accept or reject gambles that offered a 50/50 chance of gaining or losing money. Each gamble has an amount that can be used as target in a regression setting. The second dataset from [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF] is a visual object recognition task. Each object category can then be used as target in a classification setting. In this setting, we use a sparse logistic regression model instead of a Lasso. We refer to [START_REF] Tom | The neural basis of loss aversion in decision-making under risk[END_REF] and [START_REF] Haxby | Distributed and overlapping representations of faces and objects in ventral temporal cortex[END_REF] for a detailed description of the experimental protocol. Data are publicly available on http://openfmri.org.

Regression After standard preprocessing (slice timing, motion correction, first level analysis with a general linear model, inter-subject spatial normalization), the regression dataset consists of 17 subjects with 48 fMRI observations per subject. Only the gain condition was used (see [START_REF] Jimura | Analyses of regional-average activation and multivoxel pattern information tell complementary stories[END_REF]). The 48 observations contain 6 repetitions of the experiment with 8 levels of gain. FMRI volumes were downsampled to 4×4×4 mm voxels. When learning from the full set of 17 subjects, the training data therefore consist of 816 samples with approximately 26 000 voxels. The prediction here is inter-subject: the estimator learns on some subjects and predicts on left out subjects.

Classification In the object recognition experiment, 5 subjects were asked to recognize 8 different types of objects (face, house, cat, shoe, scissors etc.). We focus here on the binary classification task that consists in predicting whether the subject was viewing one object or the other (e.g. a house or a face). The data consist of 12 sessions, each with 9 volumes per category. When considering N sessions for training, data consist of 18N volumes each containing about 30 000 voxels.

The problem is here intra-subject: some sessions are used for learning and prediction is performed on left out sessions of the same subject.

As for such real data the ground truth is not known, we test these qualitative observations by training a predictive model on the voxels with the highest scores as in [START_REF] Haury | The influence of feature selection methods on accuracy, stability and interpretability of molecular signatures[END_REF]. The rationale is that if voxels with the highest F-scores contain more false positives than voxels with the highest stability scores then a predictive model trained with the first set should give significantly worse performance than a model trained with the second one. To validate this hypothesis, we train a ℓ 2 logistic regression model on the best voxels as scored by the methods, for different object contrasts and numbers of training sessions. The number of selected voxels and the regularization parameter is optimized by crossvalidation. For T training sessions (4 ≤ T ≤ 10) we obtain a prediction score on 12-T left out sessions, unseen by the feature-selection method. For completeness, we also benchmark a linear SVM, as well as ℓ 1 and ℓ 2 logistic regressions trained on the whole brain. Experiments were performed with the scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], using LIBLINEAR [START_REF] Fan | LIBLINEAR: A library for large linear classification[END_REF] for logistic regression and LIBSVM [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] for SVM.

Results Figure 3 shows maps of F values, the scores obtained with standard stability selection, and the scores obtained using our approach. The scores are not thresholded for visualization purposes, meaning that the zeros obtained are actually zeros. On both datasets, one can observe that despite being fairly noisy F-tests highlight well localized regions of the brain. Very similar regions are clearly outlined by our approach, whereas the standard stability selection procedure yields more spatially scattered selection scores. We also ran elastic-net in the regression settings (not displayed). The cross-validation lead setting a compromise between ℓ 2 and ℓ 1 penalty in favor of ℓ 1 : λ (0.1ℓ 2 + 0.9ℓ 1 ). As a result, the maps were overly sparse, displaying single scattered voxels. Finally, in the regres- sion and classification settings, we find ≈ 1000 and 2700 non-zero scores, which corresponds to a ratio n/n min respectively ≈ 25 and 780 for n min = 2p log(p -k). While visual inspection of the figures suggests that thresholding the scores at 0 probably leads to false detections and thus over estimating k, the estimation problems tackled here can be safely considered as well beyond the reach of conventional sparse estimators in terms of number of samples.

These prediction results are reported in Fig. 4. First we observe that a linear SVM and an ℓ 2 -logistic regression yield very similar results which is expected due to the similarity between the logistic loss and the hinge loss. Second, as this learning problem displays many noninformative features, using an ℓ 1 penalization actually helps the prediction. Finally, we observe that using the best voxels as scored by the randomized ward logistic leads to significantly better performance than when using the voxels with the highest F-score. This suggests that it actually achieves a better identification of the voxels involved in the cognitive task.

Conclusion

Motivated by the problem of brain mapping using functional MRI data we address the problem of support identification when working with strongly spatially cor-related designs as well as low number of samples. We propose to use a clustering of variables as well as randomization techniques, in order to position the problem in settings where a good recovery is achievable, i.e., low correlation between variables in the support and limited number of variables compared to the number of samples. Further informed by the strong spatial correlations between voxels in fMRI data we constrain spatially the clusters. This constraint injects additional prior in the estimation to facilitate the recovery of spatially contiguous supports. Our formulation enables the use of computationally-efficient sparse linear regression models and yields an overall algorithmic complexity linear in the number of features.

Results on simulations highlight the settings in which our approach outperforms techniques like elastic net or univariate feature screening. Results on two publicly available fMRI datasets demonstrate how our approach is able to zero out some non-informative regions while outlining relevant clusters of voxels. Although the true support is unknown for real data, prediction scores obtained on the best voxels outperforms alternative supervised methods suggesting that our approach is not only better in terms of support identification but can also significantly improve prediction on such data.

The success of our approach with very small sample sizes lies in exploiting the characteristics of the prob- lem tackled: recovery of contiguous clusters from image measurements with local correlations. Indeed by clustering variables, we exploit the fact that their observed correlation shares common structure with the sparsity pattern to recover. Further work include investigating other clustering algorithms and comparing to more computationally costly but convex estimators such as using overlapping group penalties.

Figure 1 .

 1 Figure 1. Recovery power on simulated data for p = 2048, n = 128 (left) or n = 256 (right) as a function of the cluster size and the smoothing. The recovery power is measured by the area under the precision-recall curve (AUC). Top row: contours at 0.95 and 0.75 of the AUC for different methods. Bottom row: best performing method.

Figure

  Figure 2. Recovery on simulated data for n = 128 as a function of the clustering approach: learning clusters on the original X, a random Gaussian i.i.d signal, or our randomization approach.

Figure 3 .

 3 Figure 3. Score maps (unthresholded) as estimated by different methods on fMRI datasets. Left column: regression settings. Right column: classification settings (contrast face vs house) using 4 sessions of training data.

Figure 4 .

 4 Figure4. Average classification scores obtained using full brain linear SVM, ℓ1 and ℓ2-logistic regressions and ℓ2logistic regression defined over the most significant voxels given by F-test (F & L2-Logistic) or randomized ward logistic. The average is computed over 6 binary classification tasks (face-house, face-cat, bottle-scissors, cathouse, shoe-house, shoe-face). P-values, computed with a Wilcoxon 1-sided paired, test if the randomized ward logistic prediction scores are larger than the ones obtained with the other methods. Randomized ward logistic significantly outperforms all the other methods.