Raphael F Garcia

Eelco Doornbos

Sean Bruinsma

Hélène Hebert

• Direct measurement of tsunami generated gravity waves in the thermosphere • Gravity wave arrival azimuth determined from GOCE measurements • Definition of a marker of thermospheric gravity waves inside GOCE data

Oceanic tsunami waves couple with atmospheric gravity waves, as previously observed through ionospheric and airglow perturbations. Aerodynamic velocities and density variations are computed from Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) accelerometer and thruster data during Tohoku-Oki tsunami propagation. High-frequency perturbations of these parameters are observed during three expected crossings of the tsunami-generated gravity waves by the GOCE satellite. From theoretical relations between air density and vertical and horizontal velocities inside the gravity wave, we demonstrate that the measured perturbations are consistent with a gravity wave generated by the tsunami and provide a way to estimate the propagation azimuth of the gravity wave. Moreover, because GOCE measurements can constrain the wave polarization, a marker (noted C 3 ) of any gravity wave crossing by the GOCE satellite is constructed from correlation coefficients between the observed atmospheric state parameters. These observations validate a new observation tool of thermospheric gravity waves generated by tsunamis above the open ocean.

Introduction

Long period tsunami waves couple at the ocean surface with atmospheric gravity waves propagating upward in the atmosphere [START_REF] Lognonné | Normal mode summation of seismograms and barograms in a spherical earth with realistic atmosphere[END_REF][START_REF] Occhipinti | Three-dimensional waveform modeling of ionospheric signature induced by the 2004 Sumatra tsunami[END_REF][START_REF] Hickey | Propagation of tsunami-driven gravity waves into the thermosphere and ionosphere[END_REF]. Due to kinetic energy conservation and to the exponential density decrease in the atmosphere, these waves are amplified and generate significant perturbations in the thermosphere/ionosphere system [START_REF] Peltier | On the possible detection of tsunamis by a monitoring of the ionosphere[END_REF]. Electron density perturbations due to these gravity waves were observed in GPS data [START_REF] Artru | Ionospheric detection of gravity waves induced by tsunamis[END_REF][START_REF] Rolland | Ionospheric gravity waves detected offshore Hawaii after tsunamis[END_REF][START_REF] Galvan | The 2009 Samoa and 2010 Chile tsunamis as observed in the ionosphere using GPS total electron content[END_REF][START_REF] Galvan | Ionospheric signatures of Tohoku-Oki tsunami of March 11, 2011: Model comparisons near the epicenter[END_REF][START_REF] Occhipinti | From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes[END_REF], ionospheric doppler sounders [START_REF] Liu | Giant ionospheric disturbances excited by the M9.3 Sumatra earthquake of 26 December 2004[END_REF][START_REF] Liu | Ionospheric disturbances triggered by the 11[END_REF], and ionospheric correction of altimetric satellites [START_REF] Occhipinti | Three-dimensional waveform modeling of ionospheric signature induced by the 2004 Sumatra tsunami[END_REF]. Recently, the effect of these gravity waves on the nightside airglow was modeled [START_REF] Hickey | Atmospheric airglow fluctuations due to a tsunami-driven gravity wave disturbance[END_REF] and observed [START_REF] Makela | Imaging and modeling the ionospheric airglow response over Hawaii to the tsunami generated by the Tohoku earthquake of 11 March 2011[END_REF], adding a new way to detect tsunami propagation in the thermosphere.

In this study, we demonstrate that accelerometer data from the GOCE (Gravity Field and Steady-State Ocean Circulation Explorer) satellite observed the gravity waves generated by the Tohoku-Oki tsunami (on 11 March 2011). The very low Earth orbit of the GOCE mission, at 270 km altitude, its drag compensated platform, and the unprecedented quality of its accelerometers [START_REF] Floberghagen | Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission[END_REF] allow us to recover very precisely the nongravitational forces exerted on the satellite and deduce air density and winds in the thermosphere along the satellite track [START_REF] Doornbos | Neutral density and crosswind determination from arbitrarily oriented multi-axis accelerometers on satellites[END_REF]. The nongravitational accelerations of the satellite are known to within 2 × 10 -12 m/s 2 [START_REF] Floberghagen | Mission design, operation and exploitation of the gravity field and steady-state ocean circulation explorer mission[END_REF]. The density and wind retrieval along satellite track have an uncertainty smaller than 1% and are provided with ten seconds sampling. Recently, [START_REF] Garcia | GOCE: The first seismometer in orbit around the Earth[END_REF] demonstrated that GOCE detected the infrasound created by the seismic waves of the same event. The analysis of GOCE data is pursued here in order to detect atmospheric perturbations due to the open ocean surface displacements generated by the tsunami. We first identify the expected windows in which GOCE crosses the tsunami-generated gravity wave and present the high-frequency perturbations in these windows. Then, a monochromatic plane wave analysis demonstrates that the aerodynamic velocities and air density perturbations are consistent with gravity waves created by the Tohoku-Oki tsunami. Finally, a marker of any gravity wave signal in GOCE data is designed by using expected correlations between aerodynamic velocities and time derivative of air density. The time of the gravity wave excitation by the tsunami is t1. The time of gravity wave arrival at the GOCE altitude (after about 1 h) is t2. Tsunami position at this time is provided on the right of the figure. Vg stands for the gravity wave group velocity and Dt is the time delay between the arrival of the tsunami wave at the ocean surface and the arrival of the gravity wave above that point at the altitude of GOCE. Adapted from [START_REF] Occhipinti | From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes[END_REF].

Observations of Atmospheric Gravity Waves Generated by Tohoku-Oki Tsunami

Identification of Expected Signals

The propagation time of the tsunami at the surface of the ocean is computed with Tsunami Travel Time (TTT) software [START_REF] Wessel | Analysis of observed and predicted tsunami travel times for the Pacific and Indian oceans[END_REF]. According to [START_REF] Occhipinti | From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes[END_REF], the horizontal and vertical gravity wave propagation from the ocean surface to the altitude of GOCE should be computed by taking into account the gravity wave dominant period and the altitude of the satellite. Figure 1 presents a scheme of the tsunami and associated gravity wave propagations. Gravity wave propagates in the atmosphere at group velocity Vg and horizontal phase velocity c x , whereas the tsunami wave propagates at velocity V tsunami on the ocean surface. The propagation time between the ocean surface and the GOCE altitude (t2 -t1) is about 1 h. In this case, the time delay (Dt) between the arrival of the tsunami wave at some point of the ocean surface and the arrival of the gravity wave at the altitude of GOCE above that point is fixed to 240 s, corresponding to the propagation of a gravity wave with a 10 min period up to the altitude of GOCE, according to [START_REF] Occhipinti | From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes[END_REF]. However, the uncertainty on this number is large because a proper estimate would require a full modeling of the gravity wave propagation, including wind effects. Then, the crossing points between the satellite track and the gravity wave are selected assuming that the gravity wave signals last during 1600 s following the expected first arrival time. This time range is in between 1 and 3 dominant periods of the signal. Figure 2 presents the result of this analysis for GOCE orbits following the Tohoku event. Three gravity wave crossing events are obtained on three consecutive ascending orbits following the earthquake.

GOCE Acceleration and Thruster Data

Figure 3 presents the high-frequency perturbations of thruster and aerodynamic accelerations measured during the three gravity wave crossing events. The high-frequency variations of thruster and Z and Y components of aerodynamic acceleration are obtained in the same way by subtracting a low-pass filtered version of the original signal with 650 s corner frequency or ≈1/8 of the orbit period. However, because of large amplitude high-frequency variations of the Y component of the aerodynamic acceleration at the pole, high-frequency variations of this signal are obtained by substracting a polynomial fit of order 7 along ±50 • latitude ascending branches. The crossing windows identified from the previous analysis present high-frequency variations of thruster and aerodynamic accelerations above the noise level in these midlatitude regions. However, these perturbations arrive at the satellite earlier than expected. We ascribe this arrival time discrepancy to the strong thermospheric winds blowing eastward and southward. A simple estimation of the wind effect is computed by taking into account the observation geometry (see Figure 4) and the wave parameters listed in Table 1 and assuming average wind and atmospheric parameters between 140 km altitude and GOCE's position. This minimum value of the wavefront shift due to wind is estimated to 241 km along the GOCE orbit, corresponding to an arrival on the GOCE data 30 s earlier than expected at the first gravity wave crossing.

Monochromatic Plane Wave Analysis

In the following we assume that the atmospheric gravity waves generated by tsunamis far from the seismic source can be fairly well approximated by a monochromatic plane wave excited by the movement of a tsunami wave of dominant period T w propagating along x direction. Under such an assumption, the gravity wave propagation from the ocean surface to the thermosphere can be described in two dimensions (coordinates x and z). The vertical winds are neglected, but horizontal winds along x direction (u 0x ) are included in our analysis. The velocity perturbations (U x and U z ), the pressure perturbations (P), and density perturbations (𝜌) are related to the background wind (u 0x ) and background density (𝜌 0 ) by the following equations, assuming incompressibility [START_REF] Nappo | An Introduction to Atmospheric Gravity Waves[END_REF]: (left to right) Close up on the geometry of the crossing events 1, 2, and 3, indicated in Figure 1. Satellite velocity (Vs ≈ 7.8 km/s, in black) is plotted with the gravity wave horizontal phase velocity (c x , in grey) and the background wind at the altitude of GOCE satellite (u 0 , in blue). See Table 1 for details. Contours of tsunami arrival time at ocean surface are also indicated.

𝜕U x 𝜕t + u 0x 𝜕U x 𝜕x = - 1 𝜌 0 𝜕P 𝜕x (1) 𝜕U z 𝜕t + u 0x 𝜕U z 𝜕x = - 1 𝜌 0 𝜕P 𝜕z - 𝜌 𝜌 0 g (2) 𝜕U x 𝜕x + 𝜕U z 𝜕z = 0 (3)
𝜕𝜌 𝜕t + u 0x 𝜕𝜌 𝜕x + U z 𝜕𝜌 0 𝜕z = 0 ( 4
)
where g is the gravity acceleration. We ensure that the incompressibility assumption is valid by comparing, later in the text, the vertical wave numbers to the inverse of atmospheric scale height [START_REF] Lighthill | Waves in Fluids[END_REF]. The ratio of these two quantities is on the order of 2 which satisfies the "larger than one" condition but suggests that compressibility should be taken into account in future work. Assuming that the atmosphere can be considered locally as isothermal and with constant background wind along x direction, the background density 𝜌 0 decreases exponentially with scale height H and u 0x is constant. Then, a solution of these equations is derived by assuming that the perturbation variables (U x , U z , P, and 𝜌) are the product of a complex amplitude term by a monochromatic plane wave complex exponential of the form exp[i(k

x x + k z z -𝜔t)].
The two components of the wave vector are k x and k z , and 𝜔 = 2𝜋

T w is the wave angular frequency. After some obvious manipulations, we obtain the gravity wave dispersion relation in a constant background wind and isothermal atmosphere [START_REF] Nappo | An Introduction to Atmospheric Gravity Waves[END_REF]]

N 2 Ω 2 = 1 + k 2 z k 2 x (5)
where Ω = 𝜔k x u 0x is the intrinsic frequency of gravity wave, N = √ 𝛾-1 𝛾 g H is the Brunt-Väisälä frequency, and 𝛾 is the usual ratio of specific heats. In our case, the angular frequency 𝜔 is fixed by the dominant period (T w ) of the tsunami forcing the ocean surface through the following relation 𝜔 = 2𝜋 T w . It remains constant along the vertical propagation of the wave.

In addition, we obtain the following simple relation between the amplitudes of horizontal and vertical velocities inside the gravity wave:

U z = - k x k z U x (6)
From equation ( 4), we can also relate the vertical velocities to the horizontal derivative of density perturbations by the following expression:

U z = - ΩH k x 1 𝜌 0 𝜕𝜌 𝜕x (7)
and so Equations ( 6), (7), and (8) should be verified by the air density and aerodynamic velocities obtained from GOCE measurements if the measured perturbations are due to the gravity wave generated by the tsunami. However, possible discrepancies with measurements could be ascribed to the assumptions used to obtain these equations: plane wave propagation, no attenuation, and free propagation far from the gravity wave source. Future numerical studies should include 3-D spherical geometry, attenuation and wind effects, and bottom forcing by the propagating tsunami wave in order to reproduce the waveforms of observed air density and velocity variations.

U x = ΩH k z 1 𝜌 0 𝜕𝜌 𝜕x (8)

Application to GOCE Data

Figure 4 shows the observation geometry of the three crossings, indicating the satellite horizontal velocity (V S ), the horizontal phase velocity of the gravity wave (c x ), and the background wind direction (u 0 ) at the altitude of GOCE. Figure 5 defines all the parameters in the horizontal plane for a geometry corresponding to the first crossing window. V S (≈7.8 km/s) is more than 1 order of magnitude larger than both c x and u 0 . In order to demonstrate that the measured perturbations of air density and aerodynamic velocities verify equations ( 6), (7), and (8), the monochromatic plane wave parameters should be estimated. To do so, the Brunt-Väisälä frequency and atmospheric scale height are estimated, respectively, to N 2 ≈ 1.10 -4 rad 2 ∕s 2 and H ≈ 43.5 km from the NRLMSISE-00 model [START_REF] Picone | NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues[END_REF] at the time and location of the first crossing with GOCE satellite, and with the following input parameters: 81 day average of F 10.7 = 106.7, daily F 10.7 = 131 and AP = 37. The background wind (u 0 ) is obtained from HWM93 [START_REF] Hedin | Empirical wind model for the upper, middle and lower atmosphere[END_REF]. The dominant period T w of the tsunami and its propagation direction are extracted from numerical simulations of the tsunami at the ocean surface [START_REF] Hébert | Far-field simulation of tsunami propagation in the Pacific Ocean: Impact on the Marquesas Islands (French Polynesia)[END_REF]. Then the horizontal wave vector k x should be estimated. To do so, the time variable of GOCE data (t d ) is linked to the x coordinate along the gravity wave propagation direction by the formula x = V r Sx t d where V r Sx is the projection along x of the velocity of GOCE relative to the wavefront. With such relation, the apparent period of observed perturbations (T A ) can be used to determine the horizontal wavelength of the wave: 𝜆 x = V r Sx T A , and the horizontal wave vector k x = 2𝜋 𝜆 x . The satellite velocity relative to the wavefront is

V r Sx = V Sx -u 0x -c x (9)
where V Sx , u 0x , and c x are, respectively, the satellite velocity and the background wind in the direction of wave propagation (x) and the gravity wave phase velocity. The horizontal phase velocity depends on the intrinsic period (Ω), which is fixed by the tsunami wave period (T w ), and horizontal wave number (k x ). Because V S >> c x , unless the satellite orbit is almost exactly perpendicular to the wave propagation direction, we can consider V Sx >> c x and estimate the horizontal wavelength neglecting c x in a first step. Then c x is computed, and the process is iterated. So horizontal phase velocity and horizontal wave number are determined iteratively, starting with zero phase velocity, estimating V r Sx , then k x from T A , and finally c x from the dispersion relation, and repeating the process. The iteration process quickly converges to values presented in Table (1). Using equation ( 5), the vertical wave numbers are also estimated. Amplitudes are about twice the inverse of atmospheric scale height ( 1 H = 2.29e -5 m -1 ), ensuring that the incompressibility assumption is approximately satisfied.

Once the monochromatic plane wave parameters have been estimated, relations ( 6) and ( 7) are used to produce synthetic vertical velocities expected for a gravity wave generated by the tsunami from horizontal velocity along x and from density variations measured by GOCE according to

U x z d ) = - k x k z U d x (t d ) (10)
and

U 𝜌 z (t d ) = - ΩH k x V r sx 𝜕Δ𝜌 d ∕𝜌 𝜕t d (t d ) (11)
where U d x (t d ) and Δ𝜌 d 𝜌 (t d ) are, respectively, the projection along x of the horizontal velocity perturbation and the density variations measured by GOCE.

Validation of the Tsunami Gravity Wave Detection

The aerodynamic velocity is defined as the velocity of air particles with respect to the satellite. Its component in the cross-track direction (U d C ) and air density (𝜌 d ) are obtained simultaneously by an analysis of GOCE accelerations and thruster data in the horizontal plane following [START_REF] Doornbos | Neutral density and crosswind determination from arbitrarily oriented multi-axis accelerometers on satellites[END_REF]. The same algorithm is repeated along a vertical plane to estimate simultaneously the upward aerodynamic velocity (U d z ) and air density (𝜌 d ). Then, assuming that U d C is the projection in cross-track direction of the horizontal velocity along wave propagation direction U d x , these two are related through

U d x = U d C cos(𝜃)
, where 𝜃 is angle between the cross-track direction and the tsunami gravity wave propagation direction. The variations of these parameters are obtained according to the process described Figure 6 presents high-frequency perturbations of aerodynamic velocities and relative density at the three gravity wave crossing times. The three input measurements (U d C , U d z , and Δ𝜌 d 𝜌 ) present different amplitudes and phases. In particular, the relative density perturbations present frequencies lower than aerodynamic velocities and are not in phase with velocities. However, the synthetic vertical velocities (U x and U 𝜌 z ), computed assuming that waves were produced by the tsunami, are consistent both in phase and in amplitude (within a factor 2) between them and with the vertical velocities measured by GOCE. This consistency demonstrates that these waves were created by the tsunami wave at the ocean surface.

Determination of the Gravity Wave Propagation Direction

From equations ( 10), ( 11), and relations between the various parameters in the plane wave analysis, if we assume that the tsunami period and atmospheric parameters are known, the amplitude ratio

U d z U d C
can be used to constrain the azimuth of the gravity wave propagation direction (Φ w ). To do so, we repeated the above analysis by varying the arrival azimuth of the wave between 0 • and 360 • and plotting the logarithm of data to theoretical

U d z U d C
amplitude ratios defined by the following:

J(Φ w ) = | | | | | | log 10 ( | | | | | U d z U d C ∕ -k x k z cos(𝜃) | | | | | ) | | | | | | (12)
when this parameter is zero, the measured and theoretical amplitude ratios are equal. Due to the projection of U x on the cross-track direction, the horizontal velocities inside the wave cannot be estimated when the tsunami is propagating exactly along the GOCE track, which corresponds to 𝜃 values of 90 amplitude ratios have the same sign, or

U d z 𝜕Δ𝜌 d ∕𝜌 𝜕t d ∕ -ΩH k x V r sx > 0.
After excluding the azimuth range values for which the above conditions are not fulfilled, the results for the three crossing events are presented in Figure 7. The first crossing event presents a single minimum close to the expected azimuth. The other two crossings present two minimum values, and the expected azimuth is close by to one of them in both cases. This exercise demonstrates that knowing the tsunami forcing period and atmospheric state parameters, the propagation azimuth of the atmospheric gravity wave can be constrained, thanks to the vertical to cross-track velocities amplitude ratio.

Definition of a Simple Marker of Gravity Waves Inside GOCE Data: C 3

From the plane wave analysis presented before, we deduce that when GOCE crosses a gravity wave, cross-track aerodynamic velocity, vertical aerodynamic velocity, and time derivative of air density perturbations should be correlated. Consequently, we define a simple marker of these correlations inside GOCE data by computing the absolute value of the product between the three correlation coefficients

C 3 = | | | | CC(U C ; Uz) × CC(U C ; d(Δ𝜌∕𝜌) dt ) × CC(Uz; d(Δ𝜌∕𝜌) dt ) | | | | (13)
where CC(f , g)(t d ) is the correlation coefficient between functions f and g over a given window length at time t d centered on this window. In this computation, the high-frequency relative density perturbations are computed as described above, and the high-frequency variations of aerodynamic velocities are obtained by simply removing a trend inside each correlation window before computing the correlation coefficient. When GOCE is crossing a gravity wave, C 3 should be close to 1 and 0 from random perturbations. In particular, C 3 decreases drastically as soon as one of the three correlation coefficients is far from ±1. Figure 8 presents the input data and C 3 computations during three orbits after the Tohoku event. C 3 presents values above 0.2 at the three crossing times with the tsunami-generated gravity wave, meaning that the logarithmic average of the three correlation coefficients is larger than 0.6. Another gravity wavefront may have been crossed 260 min after the quake on the other side of the planet during the descending branch. These waves are not related to the quake and tsunami, and this observation reveals the potential of C 3 marker for gravity wave detection.

Conclusions

High-frequency perturbations of GOCE accelerometer and thruster data were shown to be due to atmospheric gravity waves created by the Tohoku-Oki tsunami propagating in the Pacific Ocean. Our demonstration relies on the aerodynamic velocities and density perturbations computed from GOCE data and uses the theoretical relations between these parameters expected for gravity waves generated by the tsunami. Despite limitations due to the plane wave hypothesis and to assumptions on gravity wave excitation and propagation, the GOCE data are in agreement with these theoretical relations. In addition, if tsunami forcing frequency and atmospheric state parameters are known, aerodynamic velocities allow estimating the propagation azimuth of the tsunami-generated gravity wave. Moreover, a new marker of gravity wave signals inside GOCE data (C 3 ) is designed by the product of the three correlation coefficients computed between vertical and cross-track aerodynamic velocities and the time derivative of relative air density perturbations. This marker has the potential to detect any gravity wave signal inside GOCE data.

The detection of tsunami-generated gravity waves by GOCE provides a new observation tool for these waves. Compared to GPS ionospheric sounding and nightside airglow emissions, it is not limited by the integration along the line of sight, the complex response of the wave marker (either electron density or airglow emission) to the neutral wave, or the observation conditions related to the background ionosphere state. Future modeling studies of tsunami-generated gravity waves should include lateral atmospheric variations and background winds in order to properly reproduce the arrival times at the satellite altitude and the observed waveforms.

Figure 1 .

 1 Figure1. Scheme of tsunami and gravity wave propagation in the vertical plane along tsunami propagation direction (hereafter noted as x). The time of the gravity wave excitation by the tsunami is t1. The time of gravity wave arrival at the GOCE altitude (after about 1 h) is t2. Tsunami position at this time is provided on the right of the figure. Vg stands for the gravity wave group velocity and Dt is the time delay between the arrival of the tsunami wave at the ocean surface and the arrival of the gravity wave above that point at the altitude of GOCE. Adapted from[START_REF] Occhipinti | From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes[END_REF].

Figure 2 .

 2 Figure 2. Expected crossing points (circles) between GOCE tracks and the gravity waves generated by the tsunami. Contour lines and color code represent the propagation time of the tsunami at the ocean surface computed by TTT (thick lines every hour). Circles filled with the tsunami propagation time color code are superimposed on the satellite tracks (black lines) in order to indicate the gravity wave crossing areas. The three crossing events are numbered in chronological order.

Figure 3 .

 3 Figure3. to right) Signals at the three crossing windows (dashed lines) identified by the previous analysis. (top to bottom) High-frequency variations of GOCE thruster (in nanonewton) and spacecraft body-fixed Z and Y components of aerodynamic accelerations (in meter per second squared). Black lines are signals along ascending branches corresponding to the three crossing events. Grey lines are signals for three ascending branches before and two after the crossing events.

Figure 4 .

 4 Figure 4.(left to right) Close up on the geometry of the crossing events 1, 2, and 3, indicated in Figure1. Satellite velocity (Vs ≈ 7.8 km/s, in black) is plotted with the gravity wave horizontal phase velocity (c x , in grey) and the background wind at the altitude of GOCE satellite (u 0 , in blue). See Table1for details. Contours of tsunami arrival time at ocean surface are also indicated.

Figure 6 .

 6 Figure 6. Crossing events 1, 2, and 3. (top) For each event, velocities and relative density variations deduced from GOCE accelerations and thruster data. (bottom) Comparison between vertical velocity deduced from GOCE accelerations (grey line) and synthetic vertical velocities deduced from horizontal velocity (thin black line) and from density variations (thick black line).

Figure 7 .

 7 Figure 7. Cost function (J(Φ w ), open circles) for the three gravity wave crossing events as a function of gravity wave propagation azimuth (Φ w , in degrees). Vertical dashed lines indicate expected values listed in Table Vertical thick lines indicate singularities due to gravity wave propagation along GOCE orbit.

Figure 8 .)

 8 Figure 8. (top) Vertical (Uz) and cross-track (U C ) aerodynamic velocities (in meter per second) and high-frequency variations of the relative density (×1000) and its time derivative (×10 000). (middle) Absolute values of CC(U C ; Uz), CC(U C ; d(Δ𝜌∕𝜌) dt ), and CC(Uz; d(Δ𝜌∕𝜌) dt ) correlation coefficients. (bottom) C 3 marker. Correlation window is 4 min length. Vertical dashed lines indicate the tsunami gravity wave crossing windows.

Table 1 .

 1 Wave and Geometry Parameters of Each Gravity Wave Crossing by GOCE Figure5. Observation geometry in the horizontal plane, aligned with the direction of propagation of the gravity wave generated by the tsunami. The satellite velocity (V S ) is not to scale. The horizontal phase velocity of the gravity wave, the horizontal wind, and the cross-track wind velocity measured by GOCE are c x , u 0 , and U c , respectively. The angles between gravity wave propagation direction and cross-track direction and northward direction are 𝜃 and Φ w .

	Crossing	T w	k x Azimuth T A	u 0	Meridional	𝜆 x	c x	k z
	Number	(min)	Φ w (deg)	(s)	Zonal	(m/s)	(km)	(m/s)	(rad/m)
	1	10	110	40	92	-56	177	295	-3.94e-5
	2	25	190	55	100	-59	423	282	-3.87e-5
	3	35	240	90	104	-47	244	116	-4.83e-5

  • and 270 • , or Φ w values around 170 • and 350 • . The following additional constraints allow restricting the parameter space of Φ w consistent with a gravity wave propagation: data and theoretical

	U d z U d C k z cos(𝜃) > 0; Ω is positive (necessary for gravity wave propagation); k z is real (necessary for gravity amplitude ratios have the same -k x U d z wave propagation); data and theoretical sign, or U d z U d ∕ C 𝜕Δ𝜌 d ∕𝜌
	𝜕t d

Acknowledgments

This study was funded by CNES/TOSCA through space research scientific projects and by ESA. We acknowledge Frédéric Nahan who started this analysis during his master's degree research training period at IRAP. The processing of GOCE data for retrieval of thermosphere density and wind data was funded by ESA's Support To Science Element (STSE) program.


Pierre Latouche

Pierre-Alexandre Mattei
email: pierrealex.mattei@gmail.com.

Charles Bouveyron

Julien Chiquet

Combining a Relaxed EM Algorithm with Occam's Razor for Bayesian Variable Selection in High-Dimensional Regression

Keywords: EM algorithm, high-dimensional data, linear regression, Occam's razor, spikeand-slab, variable selection

We address the problem of Bayesian variable selection for high-dimensional linear regression. We consider a generative model that uses a spike-and-slab-like prior distribution obtained by multiplying a deterministic binary vector, which traduces the sparsity of the problem, with a random Gaussian parameter vector. The originality of the work is to consider inference through relaxing the model and using a type-II log-likelihood maximization based on an EM algorithm. Model selection is performed afterwards relying on Occam's razor and on a path of models found by the EM algorithm. Numerical comparisons between our method, called spinyReg, and state-of-the-art high-dimensional variable selection algorithms (such as lasso, adaptive lasso, stability selection or spike-and-slab procedures) are reported. Competitive variable selection results and predictive performances are achieved on both simulated and real benchmark data sets. An original regression data set involving the prediction of the number of visitors of the Orsay museum in Paris using bike-sharing system data is also introduced, illustrating the efficiency of the proposed approach. An R package implementing the spinyReg method is currently under development and is available at https://r-forge.r-project.org/projects/spinyreg.

Introduction

Over the past decades, parsimony has emerged as a very natural way to deal with highdimensional data spaces [START_REF] Candès | Mathematics of sparsity (and a few other things)[END_REF]. In the context of linear regression, finding a parsimonious parameter vector can both prevent overfitting, make an ill-posed problem (such as a "large p, small n" situation) tractable, and allow to interpret easily the data by finding which predictors are relevant. The problem of finding such predictors is referred to as sparse regression or variable selection and has mainly been considered either by likelihood penalization of the data, or by using Bayesian models.

Penalized likelihood

The most natural sparsity-inducing penalty, the 0 -pseudonorm, is linked to the Akaike information criterion [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF] and to optimal subset selection. As proven by [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF], it unfortunately leads to an NP-hard optimization problem that is intractable as soon as the number of predictors exceeds a few dozens. To overcome this restriction, convex relaxation of the 0 -pseudonorm, that is, 1 -regularization, have become a basic tool in modern statistics. The most spread formulation of the 1 -penalized linear regression was introduced by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] as the "least absolute shrinkage and selection operator" (lasso) and by [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] as "basis pursuit" in a signal processing framework. Several algorithms allow fast computations of the lasso, even when the number of predictors largely exceeds the number of observations. Among them is the popular least angle angle regression algorithm (LARS) [START_REF] Efron | Least angle regression[END_REF]. The Dantzig selector, introduced by [START_REF] Candès | The Dantzig selector: Statistical estimation when p is much larger than n[END_REF] as a refined 1 -regularization problem, gives good variable selection performances while simply involving the resolution of a linear program. However, as proved by [START_REF] Zhao | On model selection consistency of lasso[END_REF], the crude lasso is not model-consistent unless some cumbersome conditions on the design matrix. Moreover, [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] showed that it can be sensitive to highly correlated predictors and [START_REF] Pötscher | On the distribution of penalized maximum likelihood estimators: The lasso, scad, and thresholding[END_REF] warned that its distributional properties can be surprisingly complex. A large number of proposals have been made to enhance the lasso as a selection operator. The adaptive lasso of [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] is a weighted version enjoying nice oracle properties that works extremely well in practice. "Bolasso", introduced by [START_REF] Bach | Bolasso: model consistent lasso estimation through the bootstrap[END_REF], achieves model consistency by combining the lasso with a bootstrap step. In a similar fashion, the stability selection of [START_REF] Meinshausen | Stability selection[END_REF] applies many lasso procedures with randomized weights on subsamples of the original data. This technique leads to an effective model selection, even in the presence of correlated predictors.

Bayesian modelling

Bayesian models have also been widely studied in a variable selection context (see [START_REF] O'hara | A review of bayesian variable selection methods: what, how and which[END_REF] for a recent review). However, most Bayesian techniques have difficulties in treating the case where the number of observations is smaller than the number of predictors (the so called "large p, small n" situation), mostly because of the exponential growth of the number of possible models (p predictors lead to 2 p models). Another drawback is the fact that the most classical linear regression prior, Zellner's g-prior (for example reviewed and improved by [START_REF] Liang | Mixtures of g-priors for bayesian variable selection[END_REF]), involves to invert the Fisher information matrix which is impossible in a "large p, small n" situation. Even though some regularization attempts of the g-prior have been made by [START_REF] Baragatti | A study of variable selection using g-prior distribution with ridge parameter[END_REF], the most efficient high-dimensional Bayesian techniques essentially rest on spike-and-slab procedures.

Spike-and-slab models, first introduced by [START_REF] Mitchell | Bayesian variable selection in linear regression (with discussion)[END_REF], use mixtures of two distributions as priors for the regression coefficients: a thin one, corresponding to irrelevant predictors (the spike, typically a Dirac law or a Gaussian distribution with small variance) and a thick one, corresponding to the relevant variables (the slab, typically a uniform or Gaussian distribution of large variance). Notably, the refined spike-and-slab model of Ishwaran and Rao (2005a) or the PAC-Bayesian approach of [START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF] have been particularly efficient even in very high-dimensional settings. Markov chain Monte Carlo (MCMC) methods have been usually chosen to select models with the highest posterior distributions. MCMC techniques, reviewed for example by [START_REF] Robert | Monte Carlo statistical methods[END_REF], have an important computational cost and may suffer, as underlined by [START_REF] O'hara | A review of bayesian variable selection methods: what, how and which[END_REF], from poor mixing properties in the case of spike-and-slablike priors. A few deterministic methods have also recently been proposed to tackle this issue. The expectation propagation (EP) algorithm was applied to perform approximate inference for group feature selection with a spike-and-slab model by [START_REF] Hernández-Lobato | Generalized spike-and-slab priors for bayesian group feature selection using expectation propagation[END_REF]. The expectation maximization (EM) algorithm was used by [START_REF] Ročková | Emvs: The em approach to bayesian variable selection[END_REF] in the case of a hierarchical Bayesian model or by Yengo et al. (2014a) in the case of a multi-slab empirical Bayes framework.

Our approach

As an alternative, our approach uses spike-and-slab-like priors induced by a binary vector which segregates the relevant from the irrelevant predictors. Such vectors, introduced by [START_REF] George | Variable selection via gibbs sampling[END_REF] have been widely used in the Bayesian literature, but have always been considered as random parameters. In most Bayesian contexts like the (hierarchical) ones of [START_REF] George | Variable selection via gibbs sampling[END_REF] and Ishwaran and Rao (2005b) or the (empirical Bayes) one of [START_REF] George | Calibration and empirical bayes variable selection[END_REF], such a binary vector would be classically endowed with a product of Bernoulli prior distributions. In a PAC-Bayesian perspective, more complex prior distributions used for example by [START_REF] Alquier | Pac-bayesian bounds for sparse regression estimation with exponential weights[END_REF] or [START_REF] Rigollet | Exponential screening and optimal rates of sparse estimation[END_REF] led to precise oracle inequalities and competitive predictive performances. In our work, the originality is to consider a deterministic binary vector, and to relax it in order to rely on an EM algorithm. This relaxed procedure allows us to find a family of p models, ordered by sparsity. Model selection is performed afterwards by maximizing the marginal likelihood over this family of models. This way to treat some parameters in a Bayesian way, and others in a frequentist one, is particularly motivated by the unifying multi-level inference approach advocated by [START_REF] Guyon | Model selection: Beyond the bayesian/frequentist divide[END_REF] and by recent advances in Bayesian theory on the merging partly frequentist empirical Bayes methods and classical hierarchical Bayesian approaches [START_REF] Scott | Bayes and empirical-bayes multiplicity adjustment in the variable-selection problem[END_REF][START_REF] Petrone | Bayes and empirical bayes: do they merge?[END_REF].

The remainder of this document is organized as follows. In Section 2, a sparse generative model is defined and the general properties of its posterior distribution are exhibited.

Section 3 shows how a relaxation of this model is considered in order to perform inference through an EM algorithm. Section 4 explains the model selection procedure of our approach and gives details about Occam's razor automatic selection as well as a link with classical frequentist penalized estimators. In Section 5, a new algorithm, called "spinyReg", for variable selection in high-dimensional regression is introduced. Section 6 presents a benchmark comparison between spinyReg and classical frequentist and Bayesian variable selection procedures, real and simulated data sets are considered. In Section 7, an original high-dimensional regression database, called "OrsayVelib", is introduced and is used to demonstrate the efficiency of our approach.

Notation

Vectors and matrices are denoted by bold cases. Given a vector x ∈ R p , we define its Euclidean norm as

||x|| 2 = ( p i=1 |x i | 2 ) 1/2
, his support as Supp(x) = {i ∈ {1, ..., p}, x i = 0}, and its 0 -pseudonorm as ||x|| 0 = #Supp(x), where x i denotes the i-th coordinate of x.

We write M n,p the set of real matrices of dimension n × p. Given a vector x ∈ R n , we denote diag(x) the matrix of M n,n with diagonal x. For two matrices A and B of M n,p , we define their Hadamard product as A B = (a ij b ij ) i≤n,j≤p where a ij and b ij respectively denote the (i, j)-th coordinate of A and B. The identity matrix of dimension n is denoted by I n . Given a binary vector z ∈ {0, 1} p , we denote z the binary vector of {0, 1} p whose support is exactly the complement of Supp(z). Given a binary vector z ∈ {0, 1} p and a matrix A ∈ M n,p , we denote A z the extracted matrix of A where only the columns corresponding to the nonzero indexes of z have been kept. Given a mean vector µ ∈ R n and a positive definite covariance matrix S ∈ M n , the density of the normal distribution is denoted N ( • ; µ, S). Given a real number y, δ y denotes the Dirac function with mass at y.

A sparse generative model

This section introduces a sparse generative model based on a spike-and-slab-like prior, and describes the general properties of its posterior distribution. Links with related models are also discussed.

The model

     Y = Xβ + ε β = z w, (1)
where Y ∈ R n is the set of n observed responses, X ∈ M n,p (R) is the design matrix with p input variables. The vector ε is a noise term with p(ε|γ) = N (ε; 0, I n /γ). A prior distribution p(w|α) = N (w; 0, I p /α) with an isotropic covariance matrix is further assumed. Moreover, we denote by z ∈ {0, 1} p a binary deterministic parameter vector, whose nonzero entries correspond to the active variables of the regression model. It is worth noticing that such modeling induces a spike-and-slab-like prior distribution for β:

p(β|z, α) = p j=1 p(β j |z j , α) = p j=1 δ 0 (β j ) 1-z j N (β j ; 0, 1/α) z j .
(2)

However, we emphasize that, contrary to standard spike-and-slab models [START_REF] Mitchell | Bayesian variable selection in linear regression (with discussion)[END_REF] which assume a Bernoulli prior distribution over z, we see z here as a deterministic parameter to be inferred from the data. As we shall see in Section 3, this allows us to work with a marginal log-likelihood which involves an Occam's razor term, allowing model selection afterwards. In the same spirit, we do not put any prior distribution on γ nor α. Finally, the graphical model is presented in Figure 1 and we denote by q = p j=1 z j the number of relevant variables and Z = diag(z).

Posterior distribution

From now on, to simplify notations, the dependency on X in conditional distributions will be omitted.

Proposition 1. The posterior distribution of w given the data is given by p(w|Y, Z, α, γ) = N (w; m, S),

where S = (γZX T XZ + αI p ) -1 and m = γSZX T Y.

Proof. Using Bayes' rule, we have log p(w|Y, Z, α, γ) = log p(Y|w, Z, γ) + log p(w|α)

+ K 1 = - γ 2 Y -XZw 2 2 - α 2 w 2 2 + K 2 = - γ 2 w T ZX T XZw + γw T ZX T Y - α 2 w 2 2 + K 3 = - 1 2 w T S -1 w + w T S -1 m + K 3 .
where K 1 , K 2 and K 3 are quantities that do not depend on w. Therefore p(w|Y, Z, α, γ) = N (w, m, S).

The vector m is the maximum a posteriori (MAP) estimate of β. Next proposition assures that it recovers the support of the parameter vector. Moreover, its nonzero coefficients correspond to ridge estimates with regularization parameter α/γ of the model where only the q predictors corresponding to the support of z have been kept.

Proposition 2. We have Supp(m) = Supp(z) almost surely and

m z = (X T z X z + α γ I p ) -1 X T z Y. (4)
Proof. Using (3), one can write

S -1 m = γZX T XZm + αm = γZX T Y,
which leads, by separating the lines corresponding to zero and nonzero coefficients of z, to m z = 0 and to (4). Notice that m z = 0 implies Supp(m) ⊂ Supp(z).

The vector m z therefore corresponds to the ridge estimator of the model where only the q predictors corresponding to the support of z have been kept. As a particular instance of a strictly convex bridge estimator, the coefficients of m z are almost surely nonzero [START_REF] Fu | Penalized regressions: the bridge versus the lasso[END_REF], Theorem 1), therefore Supp(m) ⊂ Supp(z) implies that m and z have almost surely same support.

Links with spike-and-slab models

Let us briefly link the proposed model to typical spike-and-slab models. The corresponding frameworks [START_REF] Mitchell | Bayesian variable selection in linear regression (with discussion)[END_REF][START_REF] Hernández-Lobato | Generalized spike-and-slab priors for bayesian group feature selection using expectation propagation[END_REF] would add a hierarchical layer above the model of Figure 1 by using a multivariate Bernoulli prior of the form

p(z) = p j=1 τ z j j (1 -τ j ) 1-z j ,
where τ = (τ 1 , ..., τ p ) ∈ [0, 1] p . However, as emphasized by [START_REF] Scott | Bayes and empirical-bayes multiplicity adjustment in the variable-selection problem[END_REF], the estimation of τ using empirical Bayes techniques can be extremely delicate and is likely to lead to poor variable selection performances. For instance, Hernández-Lobato et al.

(2013) underline the fact that, in the case of their spike-and-slab model, the maximization of the evidence led to a sub-optimal choice of the hyper-parameter τ , and therefore to poor variable selection. To avoid such drawbacks, the use of Bernoulli priors are not considered in this paper.

Inference

This section now focuses on inferring the model proposed above. To this end, w is seen as a latent variable while Z = diag(z), α, γ are parameters to be estimated from the data (X, Y) using an empirical Bayes framework.

Inference strategy and relaxation

The estimators of z, α and γ will be the ones that maximize the evidence (or type-II likelihood ) of the data:

p(Y|X, z, α, γ) = R p p(Y|X, w, z, α, γ)p(w|α)dw. (5)
Seing w as a latent variable, a natural optimization procedure is the expectation-maximization (EM) algorithm introduced by [START_REF] Dempster | Maximum likelihood for incomplete data via the em algorithm[END_REF]. However, the maximization of ( 5) would be problematic for two reasons -both linked to the discreteness of the model parameter. First, because the optimization problem in z is combinatorial and 2 p values of z are possible. Then, because in this case, the parameter space is partly discrete and all theoretical convergence properties of the EM algorithm require a continuous parameter space [START_REF] Wu | On the convergence properties of the em algorithm[END_REF][START_REF] Mclachlan | The EM Algorithm and Extensions[END_REF].

To overcome these issues, we propose to use a simple relaxation by replacing the model parameter by a vector z relaxed in [0, 1] p . This relaxation allows us to efficiently maximize the new, relaxed version of ( 5) using an EM approach.

From now on, and until the end of this section, we will only consider the relaxed model with z relaxed ∈ [0, 1] p . In order to simplify notations, we denote Z = diag(z relaxed ).

E-step

At the E-step of the relaxed EM algorithm, one has to compute the expectation of the complete data log-likelihood E w (log p(Y, w, |Z, α, γ)) with respect to the posterior distribution p(w|Y, Z, α, γ). Consequently, the parameters S and m of the Gaussian posterior (3) have to be computed at each step. Notice that these two parameters also allow us to compute a convenient expression of the evidence.

Proposition 3. The type-II log-likelihood is given by

log p(Y|Z, α, γ) = - n 2 log(2π) + n 2 log(γ) + p 2 log(α) - γ 2 Y 2 2 + 1 2 log det S + 1 2 m T S -1 m. (6)
Proof. By directly computing the integrand of (5), we find

log p(Y|Z, α, γ) = - n 2 log(2π) + n 2 log(γ) + p 2 log(α) + log w∈R p 1 (2π) p exp - γ 2 Y T Y + γY T XZw - γ 2 w T ZX T XZw - α 2 w T w dw, which leads to log p(Y|Z, α, γ) = - n 2 log(2π) + n 2 log(γ) + p 2 log(α) - γ 2 Y 2 2 + log w∈R p 1 (2π) p exp - 1 2 w T S -1 w + w T S -1 m dw,
which allows us to conclude.

Notice that, by replacing Z by diag(z), the expression (6) remains valid in the nonrelaxed binary case.

M-step

At the M-step, the expectation of the complete data log-likelihood E w (log p(Y, w|Z, α, γ))

with respect to p(w|Y, Z, α, γ), is maximized over Z, α, γ.

Proposition 4. Denoting Σ = S + mm T , the expected complete data log-likelihood is given by

E w (log p(Y, w|Z, α, γ)) = n 2 log(γ) - γ 2 Y T Y - α 2 Tr(Σ) + p 2 log(α) - p + n 2 log(2π) + γz relaxed T (m (X T Y)) - γ 2 z relaxed T (X T X Σ)z relaxed . (7)
Proof. We have log p(Y, w|Z, α, γ) = log p(Y|w, Z, α, γ) + log p(w|α). Thus, since both the prior on w and the noise are Gaussian, we can write

log p(Y, w|Z, α, γ) = n 2 log γ+ p 2 log(α)- p + n 2 log(2π)- γ 2 (Y-XZw) T (Y-XZw)- α 2 w T w.
Therefore, by expanding and computing the expectation of the expression, we find :

E w (log p(Y, w|Z, α, γ)) = n 2 log(γ) + p 2 log(α) - p + n 2 log(2π) - γ 2 Y T Y - γ 2 E w (w T ZX T XZw) + γY T XZE w (w) - α 2 E w (w T w).
From (4), we have E w (w) = m and, by using the properties of the trace operator,

E w (w T w) = E w (Tr(ww T )) = Tr(E w (ww T )) = Tr(S + mm T ) = Tr(Σ).
Thus, we will also have

E w (w T ZX T XZw) = E w (Tr(ZX T XZww T )) = Tr(ZX T XZΣ).
Moreover, since Z = diag(z relaxed ), we can compute

Y T XZm = z relaxed T (m (X T Y))
and

Tr(ZX T XZΣ) = z relaxed T (X T X Σ)z relaxed .
By replacing the values of the terms we have just computed, we eventually find the appropriate value of the evidence.

Maximizing the expectation of the complete data log-likelihood ( 7) with respect to the parameter γ, α, z relaxed leads to the following M-step updates.

Proposition 5. The values of γ, α, z relaxed maximizing (7) are

γ-1 = 1 n Y T Y + z T (X T X Σ)z -2z T (m (X T Y)) (8) α = p Tr(Σ) (9) ẑrelaxed = argmax u∈[0,1] p - 1 2 u T (X T X Σ)u + u T (m (X T Y)) (10)
Notice that the z relaxed update ( 10) is a quadratic program (QP) which is strictly convex if, and only if Σ X T X is positive definite. In fact, the next proposition assures that it is the case if and only if X has no null column. Therefore, in all practical cases, the objective function of this program is strictly convex and fast convex optimization procedures such as the L-BFGS-B method of [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] can be used.

Proposition 6. The matrix X T X Σ is positive definite if and only if X has no null column.

Proof. According to the Schur product theorem (Bapat and Raghavan, 1997, chap. 3), since X T X and Σ are positive semidefinite, X T X Σ is also positive semidefinite. Therefore, X T X Σ is positive definite if and only if its determinant is different from zero.

If one of the columns of X is null, then the same column of Σ X T X is also null and det(Σ X T X) = 0. The proposed condition is therefore necessary.

If none of the columns x 1 , ...x p of X are null, then Oppenheim's inequality [START_REF] Oppenheim | Inequalities connected with definite hermitian forms[END_REF][START_REF] Markham | Oppenheim's inequality for positive definite matrices[END_REF]

leads to det(Σ X T X) ≥ ||x 1 || 2 2 ...||x p || 2 2 det(Σ). ( 11
)
Since Σ = S + m T m, the determinant matrix lemma assures that

det(Σ) = (1 + m T S -1 m) det(S),
and, since S and S -1 are positive definite, det(S) > 0 and m T S -1 m ≥ 0. Therefore, we find

det(Σ) = (1 + m T S -1 m) det(S) ≥ det(S) > 0,
which, combined to (11), leads to det(Σ X T X) > 0. The condition is therefore also sufficient.

Pseudo-code

Algorithm 1 presents a pseudo-code for the EM algorithm of the relaxed model.

Links with automatic relevance determination

Interestingly, this relaxed model is somehow related to the automatic relevance determination (ARD) which uses a prior of the form p(β|a) = N (0; diag(a)) and for which the most classical way of inference is also an EM algorithm [START_REF] Mackay | Comparison of approximate methods for handling hyperparameters[END_REF][START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF].

However, our method avoids several drawbacks of this technique. First, we do not assume any hyperprior on z relaxed while Tipping (2001) uses a product of flat Gamma priors. More importantly, as pointed out by [START_REF] Wipf | A new view of automatic relevance determination[END_REF], the convergence of the EM algorithm is extremely slow and not theoretically guaranteed in the case of the ARD model. However, with our approach, since we only need the order of the coefficients

S = γ(ZX T XZ + αI p ) -1 ; m = γSZX T Y ; Σ = S + mm T ; // M-step
Compute α and γ using ( 8) and ( 9);

Compute ẑrelaxed using ( 10) and the L-BFGS-B method;

until convergence of the evidence;

of z relaxed (see Section 4), we do not have to wait for the full convergence of this parameter.

In practice, in all the experiments that we carried out, we only had to perform less than a few hundreds of iterations of the algorithm to obtain convergence of the evidence in order to perform variable selection. Notice that the fact that the evidence converges faster than the parameters of the model is a quite general property of EM algorithms [START_REF] Xu | On convergence properties of the em algorithm for gaussian mixtures[END_REF]. Moreover, conversely to ARD-like models, our model additionally includes a "ridge parameter" α which, according to Occam's razor (see Section 4), also controls the sparsity.

This also leads to an objective function different from the classical ARD one.

Model selection

In practice, the vector z relaxed has to be binarized in order to select the relevant input variables. A common choice would consist in relying on a threshold τ such that z j is set to 1 if z j ≥ τ , and to 0 otherwise. However, numerical experiments showed that such a procedure would lead to poor estimates of z. In order to perform an efficient variable selection, we will use the outputs of the relaxed EM algorithm to create a path of models and, relying on Occam's razor, we will afterward maximize the type-II likelihood over this path to finally select the relevant variables.

Occam's Razor

One of the key advantages of the approach proposed is that it maximizes a marginal loglikelihood, which automatically penalizes the model complexity by adding a term to the sum of squared errors.

Proposition 7. Up to unnecessary additive constants, the negative type-II log-likelihood can be written as

-log p(Y|z, α, γ) = -log p(Y|m, z, γ) + pen(z, α, γ) = γ 2 ||Y -X z m z || 2 2 + pen(z, α, γ) (12)
where

pen(z, α, γ) = -log p(m|α) - 1 2 log det S (13) = α 2 m 2 2 - log α 2 m 0 - 1 2 log det(γX T z X z + αI q ) a.s (14)
is the Occam factor.

Proof. First, replacing w by m in the log-likelihood leads to

log p(Y|m, Z, α, γ) = - n 2 log(2π) + n 2 log(γ) - γ 2 Y 2 2 - γ 2 m T ZX T XZm + γY T XZm therefore, since m T S -1 m = γm T ZX T Y = γY T XZm, we have log p(Y|m, Z, α, γ) = - n 2 log(2π) + n 2 log(γ) - γ 2 Y 2 2 - γ 2 m T ZX T XZm + m T S -1 m. Furthermore, log p(m|α) = -p 2 log(2π) + p 2 log(α) -α 2 m T m.
By summing the terms of the right-hand side of (12), we find the same expression of the type-II log-likelihood as in (6), which proves (12). To prove ( 14), let us note that

- 1 2 log det S = 1 2 log det(γZX T XZ + αI p ) = log α 2 (p -z 0 ) - 1 2 log det(γX T z X z + αI q ).
Then, since z 0 = m 0 almost surely (see Proposition 2), we find

- 1 2 log det S = log α 2 (p -m 0 ) - 1 2 log det(γX T z X z + αI q ) a.s.
which leads to ( 14).

The sparse generative model therefore automatically adds a 0 -2 penalty to the likelihood of the model at the MAP value of w. This is somehow similar to the "elastic net" penalty of [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], combined with a penalty linked to the volume of the gaussian posterior N (w; m, S). Notice that, when α is small, the Occam factor will be extremely sparsity-inducing but the coefficients will have a large variance. When α is close to 1, this penalty will lead to moderately sparse but notably shrinked solution. Moreover, if we write λ = (α -log α)/2 and κ = α/(α -log α), we obtain almost surely the expression

pen(z, α, γ) = λ (1 -κ) m 0 + κ m 2 2 - 1 2 log det(γX T z X z + αI q ),
involving a convex combination of the 0 and 2 penalties in an elastic net fashion. The elastic net can therefore be seen as some kind of strictly convex approximation of Occam's automatic penalty.

Interestingly, the term pen(z, α, γ) exactly corresponds to Occam's razor described by [START_REF] Mackay | Bayesian interpolation[END_REF] and detailed by Bishop (2006, chap. 4). Such a term has been widely used for model selection purposes. Let us emphasize that pen(z, α, γ) is related to the penalization term of the Bayesian information criterion (BIC). Indeed, if a broad Gaussian prior distribution for the vector w is considered and if the corresponding matrix S is assumed to have full rank, then Occam's razor is approximately (-1/2)q log n. Contrary to BIC which relies on an asymptotic Laplace approximation, we obtained here an analytical expression of the evidence.

In our model, the minimization of Equation ( 12) assures that the selected model realizes a tradeoff between the log-likelihood and an automatic penalty term. Note that a PAC-Bayesian study of the performance of Occam's penalty -similarly to the study of BIC-like penalties by [START_REF] Bunea | Aggregation for gaussian regression[END_REF] for instance -would be particularly interesting. However, to the best of our knowledge, such a work has yet to be done.

SpinyReg: an algorithm for sparse regression

We called our algorithm, which successively runs the Algorithm 1 and performs model selection over the path of models using Algorithm 2, spinyReg.

Prediction

The spinyReg algorithm is essentially a variable selection algorithm. In order to perform prediction, the natural estimator of the model is ẑ where m = γ(γdiag(ẑ)X T Xdiag(ẑ) + αI p ) -1 diag(ẑ)X T Y.

However, as it was stated at the end of Subsection 3.1, this estimator is exactly the ridge estimator performed on a small model where only the predictors corresponding to nonzero coefficients of ẑ are kept. Since we do not wait for the full convergence of the parameters in the EM algorithm, we would rather recommend to perform an ordinary least squares (OLS) estimation or a ridge regression with only a small amount of regularization on the same small model. This is the choice we made in the numerical simulations hereafter.

Initialization

The choice of initialization z relaxed = (1, ..., 1) appears particularly natural because it helps to avoid the unwanted apparition of true zero coefficients in z relaxed . Indeed, if a coefficient of z relaxed by the M-step update (10), then it can not go back to a positive value. This behavior is typical of ARD-like iterative procedures [START_REF] Mackay | Comparison of approximate methods for handling hyperparameters[END_REF][START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF].

Contrary to ARD models, we do not need true zeros in the vector z relaxed . Therefore, another solution to avoid their apparition would be to perform the quadratic program (10) over [η n , 1 -η n ] were (η n ) n≤1 is a vanishing real sequence. The resulting algorithm would be a generalized EM (GEM) algorithm satisfying Wu's convergence conditions [START_REF] Wu | On the convergence properties of the em algorithm[END_REF], contrary to the classical EM algorithm for ARD [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF][START_REF] Wipf | A new view of automatic relevance determination[END_REF]. However, because we do not wait for the convergence of z relaxed , setting the initial coefficients at 1 is sufficient in practice to avoid true zeros. Regarding the parameter α, the form of the Occam factor suggests that using a small value such as α = 10 -3 will lead to sparse solutions. This is the choice we made in the numerical simulations hereafter.

Computational cost

At each iteration, the most expensive step is the inversion of the p × p matrix S during the E-step. It would imply a O(p 3 ) complexity, not allowing us to deal with high-dimentional data. However, using the Woodbury identity, one can write when p > n,

S = 1 α I p + 1 α 2 ZX T 1 γ I n + 1 α XZ 2 X T -1 (XZ) .
Thus, the final computational cost has therefore a O(p 2 min(n, p)) complexity, which is more suitable for high-dimensional problems.

Overall, MCMC-based Bayesian variable selection methods for regression have a very large computational cost. To the best of our knowledge, the fastest efficient spike-andslab algorithm for linear regression is the EP procedure of [START_REF] Hernández-Lobato | Generalized spike-and-slab priors for bayesian group feature selection using expectation propagation[END_REF].

Each iteration of the EP algorithm costs O(n 2 p) operations, and in practice it needs more iterations than our relaxed EM algorithm to converge. complexity of the LARS algorithm is O(pqn + pq 2 + q 3 ) [START_REF] Bach | Optimization with sparsity-inducing penalties[END_REF]. SpinyReg therefore realizes a complexity tradeoff between slow MCMC Bayesian techniques and fast 1 -based methods.

Let us also emphasize that, whereas frequentist methods use cross-validation to optimize the prediction performance, spinyReg automatically estimates its hyper-parameters. In particular, its inference procedure includes the estimation of the penalty term α which is linked to the sparsity level. Therefore, the computational cost of spinyReg has to be compared to the one of 1 -based methods with the cross-validation included.

Path of Models

We rely on ẑrelaxed to find a path of models which are likely to have a high evidence. We build a path by assuming that the larger the coefficients of ẑrelaxed are, the more likely they are to correspond to relevant variables.

We define the set of vectors (ẑ (k) ) k≤p as the binary vectors such that, for each k, the k top coefficients of ẑrelaxed are set to 1 and the others to 0. For example, ẑ(1) contains only zeros and a single 1 at the position of the highest coefficient of ẑrelaxed . The set of vectors (ẑ (k) ) k≤p defines a path of models to look at for model selection. Note that this path allows us to deal with a family of p models (ordered by sparsity) instead of 2 p , allowing our approach to deal with a large number of input variables. Thus, the evidence is evaluated for all ẑ(k) and the number q of relevant variables is chosen such that the evidence is maximized: q = argmax 1≤k≤p p(Y|ẑ (k) , α, γ) and ẑ = ẑ(q) .

An introductory example

We consider here an introductory example which aims at highlighting the main features of the proposed approach. For this experiment, the Toeplitz simulation setup is used with p = 30, q = 5, ρ = 0.25, α = 1 and γ = 1. From this setup, two data sets were simulated with respectively n = 100 and n = 30 observations. The second setting corresponds to a difficult scenario where n = p whereas the first one should be easier to fit. Notice that the dimensionality is kept relatively low mainly for visualization purpose. Figure 3 shows the evolution of the evidence of the relaxed model along the iterations of the EM algorithm.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 0 10 20 30

Variables Zrelaxed q q q q q q q q q q q q q q q q q q q q q q -200 -180

-160 0 5 10 15 20

Variables Evidence

Toeplitz setup with ρ = 0.25, p = 30 and n = 100 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 0 10 20 30

Variables Zrelaxed q q q q q q q q q q q q q q q q q q q q q q -70 show the values of evidence computed on the path of models.

Regarding the first example, one can see that the five largest values of ẑrelaxed actually correspond to the five active variables. This confirms that spinyReg succeeds here in finding the relevant variables in the regression model. The second panel confirms that spinyReg would select five variables among the 30 original ones. On this quite simple example, spinyReg yields a true positive rate (TPR) equals to 1 and a false positive rate (FPR) equals to 0.

For the second and much more difficult situation (bottom row of Figure 4), the estimated values for z relaxed are less discriminative. Indeed, the values of ẑrelaxed are smaller than in the simpler case. However, even though the ranking of variables induced by ẑrelaxed respects the partition between active and inactive variables, Occam's razor leads to a too conservative choice and misses one active variable. On this more difficult data set, spinyReg yields a true positive rate (TPR) equals to 0.8 and a false positive rate (FPR) equals to 0.

Benchmark study on simulated data

We now compare the performance of spinyReg with three of the most recent and popular variable selection methods based on 1 regularization: the lasso of [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], the adaptive lasso of [START_REF] Zou | The adaptive lasso and its oracle properties[END_REF] and the stability selection of [START_REF] Meinshausen | Stability selection[END_REF]. We also added two very recent spike-and-slab approaches: the multi-slab framework of CLERE (Yengo et al., 2014a) and the EP procedure of [START_REF] Hernández-Lobato | Generalized spike-and-slab priors for bayesian group feature selection using expectation propagation[END_REF].

To this end, we simulated 100 data sets for each of the three simulations schemes (uniform, Toeplitz and blockwise), for three data set sizes (n = p/2, n = p, n = 2p) and two values for the correlation parameter (ρ = 0.25 and ρ = 0.75). The other simulation parameters were p = 100, q = 40, α = 1 and γ = 1. The measures used to evaluate the method performances are the prediction mean square error on test data (MSE, hereafter), the F-score (the harmonic mean of precision and recall, which provides a good summary of variable selection performances) and the estimated value of q (number of relevant predictors).

Lasso and Stability selection were trained using the R package quadrupen [START_REF] Grandvalet | Sparsity by worst-case quadratic penalties[END_REF]. We used the package parcor [START_REF] Kraemer | Regularized estimation of large-scale gene regulatory networks using gaussian graphical models[END_REF] to train the adaptive lasso and the package clere (Yengo et al., 2014b) to train CLERE. The spike-and-slab approach of Hernández-Lobato et al. ( 2013), which uses expectation propagation, will be subsequently denoted SSEP and was trained using the code available on the authors' web

n = p/2 n = p n = 2p
F-score q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 lasso adalasso stabsel clere ssep spiny q q q q q q q 0.00 0.25 0.50 0.75 1.00 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 lasso adalasso stabsel clere ssep spiny

MSE

q q q q 0 5 10 15 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q q q q q q q q q q q q q 0 5 10 15 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q q q q q q q q q q 0 5 10 15 lasso adalasso stabsel clere ssep spiny q q q q q q 0 20 40 60 80 lasso adalasso stabsel clere ssep spiny q 0 20 40 60 80 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q 0 20 40 60 80 lasso adalasso stabsel clere ssep spiny We present here only the results for two simulation setups: the "blockwise" one with ρ 0.75 and the "Toeplitz" one with ρ = 0.25. All the other results are available as supplementary material. Note that similar conclusions can be drawn on these other scenarios.

Figure 5 presents the F-score, MSE and q of the 6 studied methods for the blockwise simulation setup with ρ = 0.75 and for the three data set sizes, while Figure 6 presents these measures for the Toeplitz simulation setup with ρ = 0.25 and for the three data set sizes.

n = p/2 n = p n = 2p
F-score q q q q 0.00 0.25 0.50 0.75 1.00 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q q 0.00 0.25 0.50 0.75 1.00 lasso adalasso stabsel clere ssep spiny q q q q q q q 0.00 0.25 0.50 0.75 1.00 lasso adalasso stabsel clere ssep spiny MSE q q q q q q 0 5 10 15 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q q q q q q q q 0 5 10 15 lasso adalasso stabsel clere ssep spiny q q q q q q q q q q q q q 0 5 10 15 lasso adalasso stabsel clere ssep spiny q q q q 0 20 40 60 80 lasso adalasso stabsel clere ssep spiny q q q q q q 0 20 40 60 80 lasso adalasso stabsel clere ssep spiny q q q q q 0 20 40 60 80 lasso adalasso stabsel clere ssep spiny The first row of Figure 5 and Figure 6 gives the F-score. This measure allows us to figure out how the methods behave in terms of detection of the relevant variables. We can see that spinyReg and SSEP outperform other methods and have close variable selection performances. SpinyReg appears to be at his best in the "n = p/2" case on these runs.

The second row of Figure 5 and Figure 6 Let us finally highlight that the medium prediction rank of spinyReg is the second best, behind the adaptive lasso. Let us also emphasize that all frequentist methods were trained using cross-validation which optimizes prediction performance. Conversely, SSEP, CLERE and spinyReg automatically estimate their hyper-parameters. In particular, the inference procedure of spinyReg includes the estimations of the penalty term α which is linked to the sparsity level.

7 Prediction of the frequentation of the Orsay museum using bike-sharing data

In this section, we introduce a new regression problem, which aims at predicting the number of visitors of the Orsay museum (Paris) using the activity of the Paris bike-sharing system (Vélib' ).

Predicting a touristic index using open data

The emergence of open data systems has brought about a surge of complex data illustrating various social behaviors. In this challenging context, the analysis of bike-sharing systems (BSSs) provides a new insight into the touristic patterns of a city. We therefore wanted to see how well, in a city like Paris, bike-sharing data could predict a touristic index, such as the number of visitors of an important museum.

With nearly three million annual visitors, the Orsay museum is one of the ten most visited museums in the world [START_REF] Skeggs | Special report, visitor figures 2013[END_REF] 7).

This data set, thereafter called the "OrsayVelib" database, has several interesting aspects:

• While most "large p, small n" regression problems inherit their dimensionality from genomics or signal processing, this data set is purely related to social sciences. This illustrates the fact that modern social data can also lead to high-dimensional challenging statistical problems.

• Since the variables are the Vélib' stations, a sparse solution can be easily interpretable and visualizable. We would expect the relevant predictors to correspond -at least to some extend -to stations used by the visitors of the Orsay museum. In particular, the behavior of the stations closest to the museum are expected to be of important interest. For visualization purposes, one can plot on a map the location of the selected variables, being able to efficiently interpret the selection.

• The learning/test segregation of the data harshly punishes overfitting. Indeed, while September 2014 (the learning month) corresponded to exceptionally good weather conditions in Paris, whereas October had some rainy days. Since BSS data are naturally heavily linked to the weather, this means that overfitting algorithms will struggle with predicting the number of predictors on rainy days (such as October 8th). This interesting behavior is exhibited in the next subsection.

To illustrate the behavior of the data, Figure 7 provides the curve of the number of visitors during the learning and test phases and Figure 8 shows the loadings of four Vélib' stations during the first week of September. Two of these stations correspond to touristic areas with different behaviors: one is the closest one to the Orsay museum and one is one of the closest ones to the Eiffel tower. The other two correspond to large railway stations (which also happen to be large subway stations). We will show in the next subsection that

Figure 1 :

 1 Figure 1: Graphical representation of the sparse generative model.

Algorithm 1 :

 1 EM algorithm for the relaxed model Input: X,Y Output: z relaxed Initialize γ = 1, α = 1, z relaxed = (1, ..., 1); repeat // E-step

Figure 3 :

 3 Figure 2: Covariance structures for the simulation setup.

  ρ = 0.25, p = 30 and n = 30

Figure 4 :

 4 Figure 4: Variable selection with spinyReg on the two introductory examples (p = 30 and n = 150 or n = 30). The left panels present the values of ẑrelaxed (dark blue) and the actual binary values of z (pale blue). The right panels show the values of evidence computed on the path of models.

Figure 4

 4 Figure 4 presents the results of the application of spinyReg on those two data sets. The left panels present in dark blue the values of ẑrelaxed (sorted in decreasing order) and the corresponding true values of z (pale blue points) used in the simulations. The right panels

Figure 5 :

 5 Figure 5: Scenario "blockwise" with ρ = 0.75.

Figure 6 :

 6 Figure 6: Scenario "Toeplitz" with ρ = 0.25.

Figure 9 :

 9 Figure 9: Ground truth (dashed line) and predicted values for the number of visitors at each hour.

  provides the MSE values for the studied methods. Most of the methods perform well except stability selection and CLERE when n ≤ p. In particular, spinyReg has the best prediction performance for n = p/2 with the

		Prostate (n = 77, p = 8)	Eyedata (n = 96, p = 200)
		MSE×100 Selected variables MSE×100 Selected variables
	Lasso	63.6 ± 21.8	3.33 ± 0.877	1.26 ± 0.964	16.7 ± 5.56
	Adalasso	58.4 ± 15.9	4.42 ± 1.57	1.50 ± 1.248	2.4 ± 0.700
	Stability Selection 61.6 ± 14.4	1.94 ± 0.239	1.58 ± 0.850	1.7 ± 0.823
	Clere	59.8 ± 19.7	2.87 ± 0.825	-	-
	SSEP	56.6 ± 15.0	2.76 ± 0.474	-	-
	SpinyReg	58.3 ± 15.4	3.34 ± 0.607	1.25 ± 0.920	143 ± 9
		OzoneI (n = 162, p = 134)	DiabetesI (n = 353, p = 64)
		MSE	Selected variables MSE/1000 Selected variables
	Lasso	18.9 ± 4.96	10.3 ± 2.27	3.22 ± 0.407	7.43 ± 2.41
	Adalasso	16.84 ± 4.48	8.32 ± 3.16	3.02 ± 0.395	9.31 ± 2.25
	Stability Selection 17.9 ± 5.25	9.68 ± 1.10	2.97 ± 0.387	7.77 ± 0.423
	Clere	19.6 ± 5.48	5.43 ± 2.55	3.15 ± 0.384	2.33 ± 0.587
	SSEP	29.6 ± 10.2	74.8 ± 5.45	3.70 ± 0.647	62.0 ± 1.36
	SpinyReg	18.9 ± 5.46	10.79 ± 2.69	3.13 ± 0.376	8.5 ± 1.45
		Table 1: Results on real-world data sets
	and reminds the interesting analogy between the Occam factor (12) used in spinyReg and
	the elastic net penalty.

  Njato[START_REF] Njato Randriamanamihaga | Clustering the vélib'dynamic origin/destination flows using a family of poisson mixture models[END_REF]. The predictive variables that will interest us for our regression problem are the percentages of parked bikes (or loadings) for all the Vélib' stations of Paris. These percentages are available through the open data API provided by At each hour, the number of visitors present in the museum constitutes the response variable of our regression problem. The predictors are the loadings at each hour of the p = 1158 Vélib' stations in Paris. Only the hours corresponding to opening days (from 8am to 6pm, except Mondays) of the museum are kept. The month of September 2014 constitutes the learning set (with n = 316 observations), and the first two weeks of October 2014 the test set (see Figure

	JCDecaux 1 .
	7.2 The "OrsayVelib" database

. Known for having the vastest collection of impressionist paintings in the world, it holds for example Manet's Le Déjeuner sur l'herbe or Van Gogh's Nuit étoilée sur le Rhône. The frequentation of the museum at each hour was given as a courtesy by the museum services.

The Paris bike-sharing system, called Vélib', was launched by JCDecaux and the city of Paris in 2007 and is nowadays certainly the most active BSS in Europe. Statistical studies of the Vélib' system have been for example conducted by

[START_REF] Bouveyron | The discriminative functional mixture model for the analysis of bike sharing systems[END_REF]

  Figure 7: Number of visitors during the learning and test phases. Only opening hours of the museum (8am to 7pm, from Tuesday to Sunday) are shown.Figure 8: Loadings of four Vélib' stations during the first week of September. Only opening hours of the museum (8am to 7pm, from Tuesday to Sunday) are shown.. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07 oct. 08 oct. 09 oct. 10 oct. 11 oct. 12 oct. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07 oct. 08 oct. 09 oct. 10 oct. 11 oct. 12 . 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07 oct. 08 oct. 09 oct. 10 oct. 11 oct. 12

												Number of visitors at the Orsay museum
									Learning Ridge (MSE=143.61, nvar=1158)		Test
		6000
				4000
	visitors	4000	0 2000 sept. 30 octvisitors
		2000
											Lasso (MSE=132.08, nvar=167)
		0
				Aug 31 4000	Sep 03	Sep 05	Sep 07	Sep 10	Sep 12	Sep 14	Sep 17	Sep 19	Sep 21	Sep 24	Sep 26	Sep 28	Oct 01	Oct 03	Oct 05	Oct 08	Oct 10	Oct 12
			0 2000 sept. 30 visitors
								Orsay museum	SSEP (MSE=146.4, nvar=1146)	Eiffel tower
				1.00											1.00
			0.00 0.25 0.50 0.75 4000 Bike station loading 0 2000 sept. 30 octvisitors								Bike station loading	0.00 0.25 0.50 0.75
					sept. 30 oct. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07			sept. 30 oct. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07
				1.00 6000		Paris Nord railway station SpinyReg (MSE=127.36, nvar=45) Paris Montparnasse railway station
			visitors	0.25 0.50 0.75 4000 Bike station loading 2000										Bike station loading	0.1 0.2 0.3
				0
				0.00	0.0 sept. 30 oct. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07 oct. 08 oct. 09 oct. 10 oct. 11 oct. 12
					sept. 30 oct. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07			sept. 30 oct. 01 oct. 02 oct. 03 oct. 04 oct. 05 oct. 07

The real time data are available at https://developer.jcdecaux.com/ with an API key.

Algorithm 2: Model selection algorithm

Input: X,Y,α,γ,ẑ relaxed

Output: z for k = 1..p do Compute ẑ(k) ; q = argmax 1≤k≤p p(Y|ẑ (k) , α, γ) ; ẑ = ẑ(q) ;

Numerical comparisons

In this section, we illustrate the behavior of spinyReg on simulated and real data sets, and compare it to the most efficient state-of-the-art methods.

Simulation setup

In order to consider a wide range of scenarios, we use three different simulation scenarios:

"uniform", "Toeplitz" and "blockwise". The simulation of the parameter w and of the noise ε is common for the three schemes: w ∼ N (0, I p /α) and ε ∼ N (0, I n /γ). The design matrix X is simulated according to a Gaussian distribution with zero mean and a covariance matrix R depending on the chosen scheme. The correlation structure of R = (r ij ) i,j=1,...,p is as follows:

• "uniform": r ii = 1 for all i = 1, ...p and r ij = ρ for i, j = 1, . . . , p and i = j,

• "Toeplitz": r ii = 1 for all i = 1, ...p and r ij = ρ |i-j| for i, j = 1, . . . , p and i = j,

• "blockwise": R = diag(R 1 , ..., R 4 ) is a 4-blocks diagonal matrix where R is such that r ii = 1 and r ij = ρ for i, j = 1, . . . , p/4 and i = j.

These three correlation structures are represented on Figure 2.

Then, Z is simulated by randomly picking q active variables among p. The predictive vector Y is finally computed according to Equation (1).

highly correlated blockwise case.

The last row of Figure 5 and Figure 6 gives the number q of active variables estimated by the 6 methods. We remind that the actual number of active variables is q = 40 for these simulations (represented by the dashed lines on Figure 5). It is worth noticing that lasso has a clear tendency to overestimate the number of active variables, particularly when n becomes large. Conversely, stability selection has the opposite behavior and underestimates q. Its very conservative behavior has the advantage that it avoids false-positives. It turns out that spinyReg provides consistently a good estimate of the actual value of q.

Study on classical regression data sets

We now consider four real-world data sets: the classical prostate data set used for example by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], the eyedata data set of [START_REF] Scheetz | Regulation of gene expression in the mammalian eye and its relevance to eye disease[END_REF], which contains gene expression data of mammalian eye tissue samples, the OzoneI data set included in the spikeslab package [START_REF] Ishwaran | spikeslab: Prediction and variable selection using spike and slab regression[END_REF] and which uses the ozone data set of [START_REF] Breiman | Estimating optimal transformations for multiple regression and correlation[END_REF] with some additional interactions and the DiabetesI data set which is also available in the spikeslab package and uses the diabetes data set of [START_REF] Efron | Least angle regression[END_REF] with some additional interactions. Applying the same methods as before, we trained our data randomly using 80% of the observations and computed the test error on the remaining data. Repeating this procedure 100 times, we computed the mean and the standard deviation of the test error and of the number of variables selected. Results are reported in Table 1. We did not compute the test error for methods which did not succeed in selecting variables.

We can see that spinyReg obtains competitive predictive results on all data sets. Moreover, we can note that it is less conservative than most other algorithms. On the challenging eyedata data set for example, while the two other Bayesian methods fail to select at least one variable, spinyReg selects three quarters of the predictors and has the lowest MSE.

Results

We applied the algorithms of Section 6 to the OrsayVelib database. Since the sparsity of this regression problem is not absolutely certain, we also added a non-sparse method to the benchmark: ridge regression with a cross-validated regularization parameter. The test errors and sparsity patterns obtained are detailed in Table 2 (for the sake of clarity, only the five best methods are displayed). One can notice that spinyReg has the lowest generalization error and that it selects fewer variables than its competitors.

Figure 9 allows to compare the true number of visitors during the test phase with the predicted values of the four methods. We can notice that, as expected, all algorithms struggle with October 8th, which was a rainy day. On this specific day, spinyReg is (especially in the afternoon) the closest one to the truth. In a similar fashion, spinyReg is the only method that accurately predicts the small augmentation of the first three days of October.

Eventually, one can plot the location of the selected variables on the map of Paris. For the sake of clarity, we only did it for the two best methods: lasso and spinyReg. Figure 10 presents the maps of selected stations by both methods. Green dots correspond to positive coefficients and red dots to negative coefficients. The size indicates the magnitude of the coefficient (the larger the dot, the larger the absolute value of the coefficient). The black dot corresponds to the location of the Orsay museum.

The lasso selection appears to be very broad and difficult to interpret. In particular, the lasso does not select the closest station to the museum. Conversely, the spinyReg selection is more interpretable: one can see that it does select the closest stations to the museum, and that their regression coefficients are positive (which means that these stations are q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 48.81 lon lat q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 48.81 likely to be full when the museum is crowded). Around the neighborhood of the museum, there is a ring of stations with almost exclusively negative coefficients (Eiffel tower, Paris

Nord and Montparnasse railway stations, place de la Bastille) which can be interpreted as stations from where the visitors of the museum rent their bikes. Beyond this ring, the selected stations essentially correspond to popular public parks (bois de Vincennes, parc

Montsouris, parc André Citroën, bois de Boulogne). This is not surprising since their frequentation is also linked to the touristic activity of the city.

As a summary, spinyReg both succeeds in providing an interpretable selection of Vélib' stations while having the most effective prediction performance.

Conclusion

We considered the problem of Bayesian variable selection for high-dimensional linear regression through a sparse generative model. The sparsity is induced by a deterministic binary

SUPPLEMENTARY MATERIALS

Additional benchmark results: Boxplots corresponding to all scenarios described in Section 6 (.pdf file)