Samuel Auclair
email: s.auclair@brgm.fr

Xavier Goula

Jose-Antonio Antonio Jara

Yolanda Colom

Feasibility and Interest in Earthquake Early Warning Systems for Areas of Moderate Seismicity: Case Study for the Pyrenees

Keywords: earthquake early warning, real-time seismology, moderate seismicity, seismic risk, magnitude assessment

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

SUMMARY

Seismic risk in areas of moderate seismicity (AMS) is a reality and destructive earthquakes may occur there. This is why stakeholders and risk managers in these areas show an increasing interest in Earthquake Early Warning Systems (EEWS) and are wondering about the possibility to equip their territories. This question of feasibility has to be broached globally: distinguishing technical feasibility from notions of opportunity and utility, and thus has to involve both scientific experts and stakeholders. Moreover it is necessary to consider each situation specifically without any preconceptions, specifically concerning potential uses that can be made of early warnings: in particular, it is clear that when applied to AMS, the EEWS principles would translate into shorter lead-times than the ones expected for more seismic areas and would consequently require considering new ways of using early warnings. Through the experience conducted for the Pyrenees (France-Spain border) as part of the SISPyr project, this study aims at identifying responses to the potential value of EEWS as real-time seismic-risk mitigation tools for AMS as well as to offer a methodological framework to guide stake-holders in their assessment of opportunities to set up EEWS.

INTRODUCTION

As real-time seismology progresses, the principle of earthquake early warning becomes conceivable for AMS where destructive earthquakes may occur and that are continuously becoming more vulnerable to seismic risk as populations and critical infrastructures are concentrating there. Like in Switzerland [START_REF] Behr | Evaluation of Real-Time and Off-Line Performance of the Virtual Seismologist Earthquake Early Warning Algorithm in Switzerland[END_REF] and Spain [START_REF] Carranza | Earthquake early warning for southern Iberia: A P wave threshold-based approach[END_REF][START_REF] Romeu | An Earthworm based prototype of an EEWS for SW Iberia[END_REF], scientists are frequently called to study the feasibility of establishing these areas with EEWS based on existing seismological networks. This firstly requires investigating the technical feasibility of such systems in areas where recurrence intervals associated with destructive earthquakes are long, and then to assess their related theoretical performances. On the basis of these elements it is important to identify the potential end-users of EEWS in order to initiate a dialog and to broach with them the use that they could make of early warnings.

One of the main scientific issues about seismic early warning deals with the ability to quickly assess the destructive potential of imminent seismic shaking. When considering regional EEWS, this problem amounts to quickly assessing the magnitude of earthquakes from the first seconds of the waveforms. While many authors have demonstrated the ability to rapidly assess magnitudes greater than 4 and up to 6 or 7 thanks to the first seconds of the P waves and they have established either regional (e.g. [START_REF] Lancieri | Magnitude scaling of early-warning parameters for the Mw 7.8 Tocopilla, Chile, earthquake and its aftershocks[END_REF] or global scaling relations (e.g. [START_REF] Zollo | A threshold-based earthquake early warning using dense accelerometer networks[END_REF], very few studies have considered lower magnitude events (e.g. [START_REF] Heidari | Magnitude-scaling relations using period parameters tau (sub c) and tau (sub p) (super max), for tehran region, Iran[END_REF]Kuyuk and Allen, 2013b). However, while particular attention must be paid to largest earthquakes (probably of the order of 6.5 in Pyrenees), which are liable to produce the greatest damage, smaller earthquakes should also be considered when working in AMS so as to be certain to be able to distinguish them during the real-time analysis.

On the other hand, early-warning should not be considered as an end in itself and its utility has to be determined with regard to its effective use for risk mitigation. Indeed EEWS may be divided into two essential components, which are the scientific analysis and the diffusion/use of the early warning [START_REF] Nigg | Risk communication and warning systems[END_REF]. Therefore, even for an early and reliable alert, if it is not relayed in a satisfactory way or if its recipients do not take appropriate actions, an EEWS could be justifiably considered as ineffective.

Nevertheless, end-user aspects are very often neglected during conception and setting up of EEWS, all the efforts being concentrated on technical questions concerning the emission of alerts. As an example, the Mexican "SAS" system, which is operating since 1991 and is able to provide alerts of nearly one minute in Mexico City, is often criticized for the absence of strategy in identification of end-users and their needs: this absence of strategy is judged as being the main limiting factor of the efficiency of the system [START_REF] Suárez | Performance evaluation of the seismic alert system (SAS) in Mexico City: A seismological and a social perspective[END_REF]. Therefore, it is indispensable not only to identify potential end-users of EEWS, but also to evaluate their perception of the system, their needs as well as their considered actions in response to early warnings.

Even though it is disconnected from the more technical aspects, this reflection is necessary to guide the conception of the system, from the choice of methods used for the data analysis to the way to use early warnings.

According to this principle and through the case study of the Pyrenees (see below), this study aims at identifying ways of responding to the potential value of EEWS as real-time seismic-risk mitigation tools for AMS, as well as to offer a methodological framework to guide stake-holders in their assessment of opportunities to set up EEWS. In the following section, we focused on "technical" feasibility aspects. The Pyrenean seismic network is first examined in order to assess its adaptability to early warning purposes. In particular, redundancy issues, network coverage, data processing and time latency of the existing real-time system are analyzed. Then, an assessment of the theoretical performance of the system is made, so as to estimate expected lead-times in the Pyrenees, thus providing a basis for reflections on how appropriate such a system may be in the area. In addition, performance criteria are proposed in order to facilitate this work for other AMS. Therefore, in the subsequent section, different rapid magnitude determination approaches are tested in order to check their adaptability to the Pyrenean context and to establish empirical relationships useful for this region. To that end, a waveform catalogue was constituted, gathering more than 2 400 records from 193 Pyrenean seismic events. Based on a survey of potential end-users from this region that aimed to identify their wishes in terms of early warning, the last section is dedicated to question of the usefulness of EEWS for the Pyrenees. The article ends with a discussion and some brief conclusions about applicability of seismic early warning for AMS in general.

GENERAL OVERVIEW OF THE CASE STUDY OF PYRENEES

The Pyrenees are a 400-km-long mountain range located in southwest Europe along the French-Spanish border, and constitute one of the most earthquake-prone regions of mainland France and Spain, with more than 400 M L ≥2.0 events per year, of which about ten are locally felt. In spite of their relatively moderate seismic activity compared to other European countries such as Romania and Italy, the Pyrenees have historically experienced numerous destructive earthquakes, including events in 1428, 1660, 1750 and 1967 which reached intensity VIII (MSK scale : [START_REF] Medvedev | Neue seismische Skala Intensity scale of earthquakes, 7. Tagung der Europäischen Seismologischen Kommission vom 24[END_REF] or more (Figure 1 -see SISFRANCE database -www.sisfrance.net, last accessed June 30 th 2014). The western part of the massif is characterized by a greater seismic activity concentrated along the North-Pyrenean Fault while the eastern part area shows a more diffuse seismicity [START_REF] Souriau | Analysis of the seismicity in the central part of the Pyrennes (France) and tectonic implications[END_REF][START_REF] Olivera | Els terratrèmols dels seglesXIV i XV a Catalunya[END_REF]. The existence of many events having caused intensity higher than VII (MSK) underlines the necessity to give attention to this area in terms of seismic risk [START_REF] Secanell | Probabilistic seismic hazard assessment of the Pyrenean region[END_REF].

The Pyrenean region is monitored by several seismological networks on both sides of the French -Spanish border, counting in total around 120 seismic stations (broad-band and strong-motion sensors taken together). Thanks to a progressive decrease of constraints associated to real-time seismology (generalization of low-cost robust data-transfer technologies, continuous increasing of data storage capacities, etc.), a growing proportion of these stations are progressively called to evolve toward real-time data transmission.

Moreover, a recent project called "SISPyr" (www.sispyr.eu, last accessed June 30 th 2014), involving (among others partners) the main owners of Pyrenean seismic stations, has notably allowed the establishment of a real-time pooling process of Pyrenean seismological data resulting in an improvement of the coverage of the massif.

A "rapid response system" (RRS) called "ISARD" [START_REF] Goula | Seismic rapid response system in the Eastern Pyrenees[END_REF] covers the eastern part of the Pyrenees, and is currently used by civil protection authorities of the Spanish autonomous region of Catalonia for crisis management purposes. It allows automatic rapid calculation and dissemination of loss assessments (both material and human losses). Moreover, the ground-motion and shaking intensity calculation code ShakeMap [START_REF] Wald | TriNet "ShakeMaps": Rapid generation of peak ground-motion and intensity maps for earthquakes in southern California: Earthquake Spectra[END_REF] has been recently implemented and adapted to the whole Pyrenean massif [START_REF] Bertil | ShakeMap implementation for Pyrenees in France-Spain border: regional adaptation and earthquake rapid response process[END_REF]www.sispyr.eu), and arouses the interest of the French and Spain civil protection authorities as a potential help for seismic crisis management. Therefore, it is quite natural that the question of the feasibility of using the current real-time monitoring network in the form of an EEWS has arisen. Of particular interest are the double issues of the technical feasibility and the potential use of such a tool within an AMS.

TECHNICAL FEASIBILITY OF AN EARLY WARNING SYSTEM

Adaptability and potential of the current real-time monitoring system for early-warning purposes

Effective network coverage for early warning

There is clear evidence that AMS such as the Pyrenees do not justify the setting-up and maintenance of seismological networks solely dedicated to earthquake early warning, and the notion of technical feasibility is, therefore, strongly interlocked with the one of the enhancement of pre-existing instrumentation. Indeed, conditional on efficient and reliable real-time data-transmission and processing technologies, both modern strong-motion and broad-band sensors currently used by monitoring seismological networks could be employed for early-warning purposes. In order to assess the adaptability of a given existing seismic monitoring network to early-warning purposes, we first define two distinct categories of technical requirements that have to be considered when implementing an EEWS (Table 1): those requirements that are essential to an EEWS ("basic requirements"), and others ("recommended requirements") that are not mandatory but that would significantly improve the reliability and efficiency of an EEWS and which may be critical for its operational use.

In the Pyrenees, most of the real-time stations fulfil the basic requirements listed in Table 1 because these stations are also used by regional/national seismological data-centers for "real-time" event detection and seismic alerts. In total, 52 stations fully accomplish the basic requirements and constitute that we call the current early-warning scenario called "EEWS1" (Table 4 in the appendix). Nevertheless, the recommended requirements are not entirely accomplished by all these stations. In particular, except for FNEB that consists of two sensors (one strong-motion plus one broad-band), there are no stations with spare data sensors or data acquisition, communication or power supply systems. To solve this problem, it is recommended to add redundancy to stations' equipment, especially to communication systems because their vulnerability is higher than the others. This redundancy also can be reached by increasing station density to allow some fault tolerance thereby reducing the effect of changes in the station distribution.

A key consideration when measuring the efficiency of any EEWS is the notion of "latency": in other words, "coverage" of the seismological network has to be no longer defined by the density of seismic stations in a given area, but rather by the one of the density of available seismic records in a given area at a given instant.

Therefore, studying the EEWS1 situation from this viewpoint, we have typical values of latency time (defined as the delay between the seismic signal being recorded at any seismic station and the time when the digitized signal is received at the data-center located in Barcelona) for 22 of the 52 considered stations.

Calculated over a time-length of 48h, most of these mean values are lower than 5 s, with an average latency time of 3.7 s. Then, the complete system latency time represented by the so-called "warning time" Tw (being the time at which the system is able to produce an alert) will be for a specific hypocenter the result of adding P-waves travel times to the closest stations and the recording time lengths, stations latency times and data processing times. Thus, under some hypotheses, this complete system latency may be regionally assessed in order to give an idea about the possible benefits of a virtual Pyrenean EEWS. Among these hypotheses, the most important one is the minimum number of stations required to perform the automatic real-time analysis (including detection, location, magnitude assessment and so forth). We decided to vary this parameter from one that corresponds to the "on-site" configuration up to four stations. In addition, a typical latency time of 3.7s (see previous paragraph) has been considered for those stations of Table 4 without any estimates for this parameter. This kind of analysis can lead to maps showing Tw as a function of the hypocentral location. As this time is not very useful in itself, we make it more explicit using the "blind zone" concept (see paragraph below).

Simulated performance analysis

When considering the application of EEWS, a very important issue is that of the "blind zone", which designates the area where a warning arrives after the destructive seismic waves (S waves followed by surface waves). Thus, the blind zone represents the area in which an EEWS is inefficient: its extension depends on many factors, such as epicentral location and focal depth of the earthquake, topology of the seismological network, latency time of close-field stations and speed of the calculation process. In a very similar way to the one used by Kuyuk and Allen (2013a) considering regional average values of seismic-wave velocity [V P =6km/s and V S =V P /1.75, [START_REF] Souriau | A new synthesis of Pyrenean seismicity and its tectonic implications[END_REF]] and a constant focal depth of 10 km, it is then possible to convert Tw-maps described previously into maps showing the spatial variation of the parameter R bz n,t . This parameter represents -for an EEWS using the first t seconds of the P-wave from at least n triggered stations to perform its analysis -the average radius of blind-zone as a function of the hypocentral location. Considering firstly that an analysis time of 2s seems long enough for the real-time assessment of magnitude in the Pyrenees (see next section), and secondly a minimum of three triggered stations, Figure 2 shows the distribution of parameter R bz 3,2 obtained for scenario EEWS1. This figure presents relatively short blind-zone radii of 30-40 km in areas with more stations, and at the same time underlines their strong spatial inhomogeneity as a function of the effective coverage of seismic monitoring.

Another thing that can be seen on Figure 2 is an increase of R bz n,t (here R bz 3,2

) very close to stations: this apparent paradox results from the observation that below a certain hypocentral distance, real-time stations are useless for early warning using fixed P-waves time-length for the analysis because of an insufficient time difference between P and S waves' arrivals for magnitude assessment. However, this effect can be easily bypassed considering shorter P-waves time-length in the epicentral area so that local magnitude can be conventionally computed for those of the stations located very close to the epicenter. R bz n,t maps may be very useful when considering the benefit of developing such a system in AMS such as the Pyrenees. Coupling these regional maps with a look at the regional seismicity is also a very valuable way to broach the question of the potential contribution over time of EEWS to real-time seismic risk management.

To that end, we generated a synthetic earthquake catalog using the "unified seismogenic zonation" PSHA and corresponding activity parameters proposed for the Pyrenees by [START_REF] Secanell | Probabilistic seismic hazard assessment of the Pyrenean region[END_REF]. The main steps taken to generate this synthetic earthquake catalog are as follows: 1) we first consider that earthquakes occurrence within a given seismotectonic zone follows a Poisson process and then we make use of the zonal annual activity rate  to calculate the number of events for each zone, 2) secondly we assume a Gutenberg-Richter relation [START_REF] Gutenberg | Seismicity of the Earth and Associated Phenomena[END_REF] for magnitude calculation using the b values, and 3) we finally assume that hypocenters occur randomly and uniformly within each zone. So as to be able to calculate robust statistics for large events, we considered an "observation" duration of 10 000 years for this catalog, resulting in 13 672 earthquakes with M L greater than 4.0 (with 907 5≤M L ≤6 events, and 55 events larger than M L 6) (Figure 2). This synthetic catalog has been created in order to easily identify trends about the usefulness of a potential Pyrenean EEWS. However, it is important to keep in mind that the results presented later in this article are strongly dependent on this catalog and might be slightly different if another catalog was used, for example using of a different seismotectonic zoning.

Considering earthquakes from the synthetic catalog (Figure 2) located inside the SISPyr project area as representative of the Pyrenean seismicity, we then calculate corresponding theoretical isoseismal curves by using the Pyrenean intensity prediction equation (IPE) of [START_REF] Bakun | Regional intensity attenuation models for france and the estimation of magnitude and location of historical earthquakes[END_REF]. These curves can be compared with warning times in order to deduce lead times (time interval between the arrival of the warning and the arrival time of the S waves) associated to given intensity levels and noted LT I n,t (Figure 3). While it is unlikely that the setting up of an EEWS justifies, in itself, the addition of many new seismic stations in areas such as the Pyrenees, it is reasonable to assume that most of the seismic networks covering AMS are likely to improve with time. In this way, we consider for the Pyrenees the on-going process of modernization of existing stations toward real-time data transmission, as well as the installation of new stations. Thus, in addition to the current real-time stations constituting the EEWS1 network, Pyrenees could be equipped in the coming years with an additional 30 more real-time stations usable for earthquake early warning:

-26 existing stations whose data-transmission in real-time is likely/possible in the medium-term ;

-4 stations whose installation is planned in the French part of Pyrenees through RESIF's permanent observation network (the French seismological component of EPOS, the European Plate Observing System -www.epos-eu.org).

We then consider all these stations as constituting an improved early warning scenario denoted "EEWS2".

Lastly, so as to assess the upper bounds of potentialities of EEWS for the Pyrenees, and even though such a configuration is unlikely in this region, we consider an "ideal" virtual network of 206 stations characterized by a triangular regular topology with mesh-size of 20 km as represented on Figure 2 and noted "EEWS3". On an inclusive basis, we assign to each station of network EEWS2 and EEWS3 a mean latency time of 2s, which is consistent with performances of current seismological data real-time transfer technologies [START_REF] Kuyuk | Designing a networkbased earthquake early warning algorithm for california; ElarmS-2[END_REF].

Similarly to scenario EEWS1, the theoretical performances of networks EEWS2 and EEWS3 have been calculated to assess the impact of the strengthening of real-time seismological monitoring of the Pyrenees in terms of seismic early warning. In order to easily compare results obtained from each of these EEWS configurations, Figure 4 represents the portion of synthetic earthquakes leading to positive lead-times for areas affected by different levels of macroseismic intensity (i.e. blind zones smaller than the considered isoseismal) as a function of magnitude, for EEWS using the first 2s of the P-waves on the records of the first three triggered stations. Evidently, it can be noticed that a densification of the EEWS seismic network improves performance. For an exposure to a given level of intensity, passing from configuration EEWS1 to EEWS2 or EEWS3 leads, in most cases, to a significant decrease in the minimum magnitude for which the first cases of positive lead times appear: thus, the minimum threshold of positive lead times corresponding to intensity VI goes from M L 5.6 for configuration EEWS1 to M L 5.2 for configurations EEWS2/EEWS3. Owing to the fact that network EEWS2 is locally as dense as EEWS3, this value of minimal magnitude changes only slightly between configurations EEWS2 and EEWS3. On the other hand, the distribution of EEWS3 stations being regular, performances associated to that configuration are more homogeneous at the scale of the whole Pyrenees, which in turn results in thinner curves than the ones corresponding to EEWS1 and EEWS2, whose dispersion indicates the spatial heterogeneity of the networks (Figure 4).

Discussion

Even though regional EEWS should ideally rely on seismic networks conceived to optimize robustness and fastness of early-warnings such as the one of the Irpinia area in Italy (ISNet, cf. [START_REF] Weber | An advanced seismic network in the southern Apennines (Italy) for seismicity investigations and experimentation with earthquake early warning[END_REF], a look at the SISPyr real-time network shows that the existing stations may be used for early warning purposes.

Indeed, the current SISPyr real-time network still satisfies an important part of EEWS requirements such as the real-time centralization of data. Moreover, even though it still has a limited coverage, the real-time network is composed of numerous seismic stations. In addition to its densification, one of the main ways of improving the network in order to implement a Pyrenean-EEWS would resort to low-latency data transmission modes, preferring to "push" data by "sample by sample" sending rather than waveform packets as it is the case today.

Then, comparing intensity predictions related to regional seismicity with theoretical blind-zone extensions shows that in some cases the current network could be used efficiently to provide early warnings of a few seconds up to around 10s in areas were intensities greater or equal to VI (and thus damages) are expected, depending of the minimum number of stations used for the analysis as well as the location of the earthquakes. Nevertheless, Figure 2 also indicates that performances of an EEWS based on the SISPyr network (configuration EEWS1) would not be homogeneous in the whole of the Pyrenees due to differences in seismic monitoring coverage. As a consequence, such a system would be much more efficient for earthquakes occurring east of the Pyrenean massif.

In the best-case scenario, if the Pyrenees were equipped with an EEWS based on a dense and homogeneous network with a low data-transmission latency (i.e. configuration EEWS3), areas potentially suffering damage (i.e. with intensities greater than or equal to VI) will likely have positive lead times for more than 70% of earthquakes with M L ≥5.6, while blind zones associated to the more frequent events with magnitudes around 5 are always larger than the damaged areas. In that case, lead times of about 10s are expected far from the epicenter, and it is highly probable that the associated intensities would be fairly weak (IV to VI). It should be borne in mind, however, that the considered IPE does not take into account neither potential site effects nor finite fault effects, which could cause damage at greater epicentral distances, in areas that could benefit from longer lead times.

Testing procedures for rapid assessment of earthquake magnitude

Waveforms catalogue

To be able to test current early-warning procedures for estimating magnitude in real time, a catalogue of representative seismic signals of the Pyrenean seismicity was first compiled. For the Pyrenees, an attempt was made to build a catalogue of waveforms, not statistically representative of seismicity all along the Pyrenean chain, but containing as many records as possible corresponding to all the magnitude ranges to be considered in order to better constrain the regression. In practice, the catalogue is composed of vertical strong-motion and broad-band seismograms, provided by the various organizations that operate stations in the region. Finally, the waveform catalogue used in this study is made up of more than 2 400 vertical records corresponding to 193 events with local magnitudes ranging between 2.0 and 5.0 (Figure 5).

Data processing and analysis

Four scaling parameters have been selected to be applied in a test stage in the Pyrenees.  c [START_REF] Kanamori | Real-time seismology and earthquake damage mitigation[END_REF] and  p max [START_REF] Allen | The potential for earthquake early warning in southern california[END_REF], that can be considered as representing respectively the predominant and the effective period of the P-waves over a fixed time window, are assumed to be empirically log-linearly correlated with magnitude (equation 1), while logarithms of peak displacement (Pd, [START_REF] Wu | Magnitude estimation using the first three seconds p-wave amplitude in earthquake early warning[END_REF]) and peak velocity (Pv, [START_REF] Wurman | Toward earthquake early warning in northern California[END_REF]) of the first seconds of the P-waves depend both on magnitude and hypocentral distance (R) (equation 2).

log10 y = a + b.Ml (equation 1)
where y represents either  c or  p max (in seconds), M L is the local magnitude, and a and b are constants to be determined through a best-fit regression analysis.

log10 y = a + b.Ml + c. log10 R (equation 2)
where y represents either Pd (in cm) or Pv (in cm/s), and a, b and c are constants to be determined through a best-fit regression analysis.

So as to calculate these scaling parameters over the waveform catalog, records have been first corrected for the instrumental response and then the P-wave arrival manually picked on the unfiltered vertical records.

Afterwards, after a simple or double integration process (depending of the type of instrument considered), records have been bandpass-filtered between 1 and 50 Hz using a Butterworth filter. Then, Pd, Pv,  c and  p max have been measured from the bandpass-filtered displacement and velocity records over a time window varying between 1 and 3 s after the first P-wave arrival. In order to avoid the "contamination" of the analyzed time-window with S-wave arrivals due to a time interval separating the onsets of the P and S waves shorter than the analysis duration, we reject all records where P and S arrivals are not far enough apart so as to compute the various selected parameters. Rather than systematically picking the S-wave onset, we consider the procedure proposed by [START_REF] Wurman | Toward earthquake early warning in northern California[END_REF] consisting in merely computing the simulated arrival time for the S waves and retaining only those records in which the interval between P-and S-wave arrivals are greater than or equal to the duration of the analysis.

In order to guarantee the quality of our analysis, we also reject all the noisy records considering a minimum signal/noise ratio (SNR) equal to 60 for  c , 30 for  p max , and 10 for peak parameters P d and P v . Finally, for the sake of homogeneity, we focus our analysis on data recorded at a maximum epicentral distance of 100 km, which covers the maximum distance considered for early warning applications in the Pyrenees.

Results & discussion

Once the scaling-parameters have been computed using the Pyrenean earthquake waveform catalogue according to the procedure described above, they have been compared with reference magnitudes given by the Spanish IGN institute, in order to establish empirical relations enabling magnitude to be estimated in real time from an analysis of the first few seconds of the P wave. The peak parameters P d and P v being a function not only of magnitude but also of hypocentral distance, we have normalized them to a reference distance (fixed to 10 km) as suggested by [START_REF] Zollo | Earthquake magnitude estimation from peak amplitudes of very early seismic signals on strong motion records[END_REF] in order to remove the dependency on distance. Thus, we can then establish correlations between magnitude and normalized peak values P d 10 and P v 10

. To reduce the scatter as much as possible, it is better to study the values of the parameters under consideration averaged for each event rather than the results obtained station by station [START_REF] Wu | Rapid assessment of damage potential of earthquakes in Taiwan from the beginning of P waves[END_REF]. To do so, we have not considered a mean of each of the scaling parameters, but rather a mean of their decimal logarithms, which are assumed to be linearly correlated with magnitude. Thus, and to avoid assigning too much weight to certain seismic traces, only those earthquakes for which we have at least two traces satisfying the selection criteria previously presented have been retained for analysis, thereby providing one mean value per event. An example of results obtained with a 3s time-window analysis are shown on Figure 6 while Table 2 summarizes all empirical relations established for the Pyrenees.

From a qualitative standpoint, examining Table 2 allows us to suggest that initial magnitude estimates are able to be derived from a very short analysis interval, which can subsequently be refined in the framework of an evolving approach. For example, shortening the analysis interval from three to two seconds for the case of a single station situated at the epicenter leads to a reduction of the blind zone by about 5 km. Moreover, it is interesting to note that the empirical correlations found for the Pyrenees are globally coherent with those found in other more seismically-active regions of the world [START_REF] Allen | The potential for earthquake early warning in southern california[END_REF][START_REF] Lancieri | A Bayesian approach to the realtime estimation of magnitude from the early P and S wave displacement peaks[END_REF][START_REF] Zollo | A threshold-based earthquake early warning using dense accelerometer networks[END_REF].

Furthermore, it is shown that although each of the parameters  c ,  p max , P d 10 and P v 10 would seem individually to be correlated with earthquake magnitude, some appear to be better proxies than others. Thus, from a statistical standpoint, a look on the coefficient of determination R 2 from Table 2 directly shows that the most satisfactory proxies for the Pyrenees are, in order, P d 10 and P v 10 (which are very similar by nature and cannot be considered as independent parameters), the frequency parameter  c , and lastly the frequency parameter

 p max
. Moreover, it should be noted that, as underlined by Kuyuk and Allen (2013b), the correlation observed on Figure 6 between  p max and M<3 is possible because the high-frequencies intrinsic to smaller earthquakes have been preserved in our analysis, while the 3Hz low-pass filter initially proposed by [START_REF] Allen | The potential for earthquake early warning in southern california[END_REF] for the calculation of  p max removes them.

Due to the limited range of magnitudes considered in our analysis, it is not surprising to observe good correlation coefficients for analysis intervals shortened to 1 s (Table 2), as the events considered correspond to relatively short rupture times. In other words, the portion of the signal being analyzed, however short it may be, bears the signature of most, if not all the rupture, and accordingly of the magnitude. A study by [START_REF] Murphy | Estimating earthquake magnitude with early arrivals: A test using dynamic and kinematic models[END_REF] showed that a 1-s analysis was long enough to assess moment magnitudes smaller than 6.0 (2s for M<6.5; 3s for M<7.0). For larger magnitudes, the authors observe a saturation of the proxy P d . Theoretically, and taking into consideration a probable maximum magnitude of 6.5 for earthquakes in the Pyrenees, a 2-s analysis of the P wave would appear to suffice to determine the magnitude of Pyrenean events in real-time.

USEFULNESS OF A PYRENEAN EEWS

Another very important issue on the evaluation of the feasibility of an EEWS deals with the question of the end-users in order to assess if such a system could answer an existing need or not. In particular, the question of "How useful is an early warning for earthquakes associated to long recurrence intervals?" is preponderant and strongly linked to the potential end-users' risk aversion. By its conception, it is consequently imperative that an EEWS should answer requirements of end-users in such a way to be adapted to their needs, and not the contrary as it is often the case. Henceforth, the approach of consultation of potential end-users is necessary so as to better understand their requirements, their expectations (or conversely their absence of expectations), as well as their needs in terms of early warning. To that end, a survey of potential French endusers had been carried out in order to evaluate their wishes in terms of earthquake early warning.

Owing to the simulated performances of a virtual Pyrenean EEWS presented in chapter 4, the Pyrenean context only allows obtaining early warnings with short lead times. Rather than conducting an inquiry in an "open" way addressing all potential end-users, we then favor focusing on actors still well accustomed to crisis management and to the taking of preventive measures, represented by the managers of critical infrastructures. Indeed, these actors are likely to be interested in short seismic early warnings and are the most likely to act as a consequence. In addition, it is important to notice that the operational establishment of an EEWS is generally mainly conditioned on the criticality of exposed elements. As a result, our survey has focused on a limited number of industries and managers of critical networks and dams. Selection of targeted infrastructures was made considering, in priority, 'facilities at risk' according to the French seismic zonation (decree No. 2010-1255 of 22 th October 2010). In addition, this list of addressees was completed with administrators of electricity, gas and high-speed train networks. Based on telephone interviews with persons in charge of security of each site/facility, this survey finally allowed gathering the opinion of eleven participants, covering the main regional managers of energy networks (electricity and gas), high-speed railways, two managers of dams representing together nearly 40 dams, as well as six industrial sites.

Results and discussion

Briefly summarized in Table 3, results of the survey highlight a favorable welcome of the idea to provide the Pyrenees with an EEWS in order to be warned of incoming strong motions, even though their destructive potential remains moderate. Indeed, despite the fact that the majority of conceivable early mitigation actions considered by participants are likely to need more time than the potential available lead times to be fully efficient, their early application seems to be considered as a way to reduce the risk of accidents as well as minimize economic losses, which are very important issues when considering strong interdependencies between critical infrastructures. Indeed, like Wu et al. (2013) who define utility curves for each « early action » as a function of available lead times, it seems to be too restrictive to consider that only fullyimplemented actions before the arrival of strong shaking are useful: that is to say, for example, that the expected benefit of an early process aiming at stopping a high-speed train is maximal if the train is still stopped when strong motions arrive, but remains positive since the speed has been significantly reduced.

Otherwise, it can be noticed that while reduction of consequences thanks to the use of an EEWS is expected by most of participants, it does not constitute a necessary condition to use it since the early warning seems to be perceived as precious information in itself, showing that the utility of EEWS does not entirely rely on the ability of users to take early actions.

It should be noted that, due to their own experience, participants tend to favor manual or semi-automatic actions over fully-automatic actions in order to react on a case-by-case basis. Nevertheless, discussing individually with each of participants let us conclude that additional automatic actions compatible with very short lead-times could be conceivable by means of a shared stage of reflection requiring their accompaniment by specialists of EEWS. This points out the importance to educate potential end-users so they can understand the potential and limitations of EEW alerts.

Another aspect that can be pointed out from this survey is that -not surprisingly -potential Pyrenean EEWS end-users express a relatively strong demand for reliability of the system. Nevertheless a lesson can be learned that a slight priority is generally given to security since missed alerts are judged to be slightly more important than false alerts: that observation comes essentially from companies with a line of business identified as "industry", for which this trend is particularly well marked (for these participants, the impact of missed alerts are on average one level higher than those of false alerts). Indeed, uncertainty being inevitably associated to each early warning, EEWS end-users have necessarily to define their own level of acceptability:

does an addressee of an early warning prefer to be warned, in a secure sense, at the risk of taking useless (costly) prevention actions in case of false alert, or does he prefer to be warned, in a certain sense, at the risk of missing relevant alerts associated with high uncertainty. It is usually qualitatively considered that the greater the effect of damage (socio-economic, environmental, etc.), the lower the acceptability of the situation, which in terms of early warning is translated by favoring false alerts over missed alarms [START_REF] Marzocchi | Probabilistic eruption forecasting and the call for an evacuation[END_REF]. It is then not surprising that so-called "industrial" users be favorable to "security" as they are submitted to the EU SEVESO Directives on industrial accidents, which in France are associated with a very low acceptance of societal risk. This requirement of reliability of early warnings implies the necessity for new EEWS to qualify the reliability of each alert by means of probabilistic assessments indicating probability of exceeding threshold values of given Intensity Measures (IMs): then these probabilities could be used within automatic Decision Support Systems (DSS) so as to perform real-time Cost-Benefit or Multi-Attribute analysis deciding the suitable actions [START_REF] Le Guenan | Participatory decision making for operational earthquake forecasting and Earthquake early warning[END_REF].

CONCLUSION AND PERSPECTIVES

Aiming at dealing with the question of adaptability of the principle of EEWS to AMS without any preconceptions and bringing together technical aspects with operational needs of end-users, this explorative study does not allow concluding definitively on the interest to endow AMS, such as the Pyrenees, with EEWS. However, it highlighted several important points:

1. The methodological framework of early assessment of magnitude seems to be adapted for AMS, making possible the discrimination of numerous smaller earthquakes and less frequent damaging events.

2. Thanks to their continuous modernization, current seismological networks in AMS, such as the ones in the Pyrenees, may be considered as a basis for the development of EEWS, when their densification suggests the possibility of greatly improving lead times.

3. There is clear evidence that the Pyrenees -and probably most AMS -contain stakeholders interested in using seismic early warnings in alternative ways.

In spite of encouraging results, it is, therefore, important to point out that unless these alternative ways are considered, many barriers still remain for the practically establishment of EEWS in AMS: strong attenuation of destruction effects with distance, important installation and maintenance costs of efficient monitoring/communication/transmission solutions, difficulty in calibrating the system, difficulty in raising awareness among end-users in the long-term, necessity to identify and to maintain fast automatic mitigation-actions and so forth.

This being said, several ways offer new alternatives allowing enlarging the range of applications of earthquake early warnings. The first one comes from the field of new technologies and, more specifically, from micro-electronic strong-motion sensors (so-called MEMS). As suggested by [START_REF] Evans | Performance of several low-cost accelerometers[END_REF], some of these MEMS sensors are sensitive and reliable enough to be used for a wide range of applications. Like the Quake Catcher Network [START_REF] Cochran | The Quake-Catcher Network: Citizen science expanding seismic horizons[END_REF] that makes use of MEMS data to produce ShakeMaps, the possibility of acquiring and networking many inexpensive sensors is giving birth to new approaches for seismic early warning [START_REF] Fleming | The selforganizing seismic early warning information network (SOSEWIN)[END_REF]Wu et al., 2013), which could possibly bypass some of the current barriers to the setting-up of EEWS in AMS.

Otherwise, indications given by potential Pyrenean end-users (Table 3) concerning on the one hand the value of the early warning in itself (independently of the way this warning could be used) and, on the other hand, the interest in engaging early application of actions in order to reduce risk of accidents and minimize economic losses, leads to a reconsideration of preconceived ideas that users may have about the utility of EEWS. Thus, it seems to be suitable to bring down the barrier between on one hand the « early stage » joining artificially and very closely early warnings with early actions, and in the other hand the « rapid stage » joining rapid-warnings and rapid-response, allowing the taking advantage of the early assessment of earthquake's focal parameters in order to improve rapidity of rapid response systems and, in doing so, acting in the first seconds after the earthquake. Should we move towards the concept of "Earthquake Early

Response Systems" for AMS? This could result in early impact assessments, e.g. ShakeMaps or PAGER [START_REF] Wald | Development of the U.S. Geological Survey's PAGER system (Prompt Assessment of Global Earthquakes for Response)[END_REF], that are of primary importance for players in emergency management, enabling them to achieve very quickly an overall view of the damage situation and to implement appropriate response strategies.

Lastly, since one of the roles of engineering seismologists working on real-time seismology is to provide decision-makers and stakeholders with a methodological framework to guide their assessment of EEWS, this study constitutes a contribution of general interest. Indeed, it generalizes the blind-zone radius parameter proposed in 2013 by Kuyuk and Allen through R bz n,t

, and introduces the lead-time parameter LT I n,t

. These parameters can be used in order to assess potential performances of an EEWS, respectively, in terms of blind-zone radius and lead time, for different configurations of the system (number "n" of triggered stations and duration "t" of analysis of the P-wave) and protection objectives (intensity "I"). Therefore, this type of "a priori" feasibility analysis may be particularly useful in guiding the decision to move or not to the establishment of a regional EEWS. Then, the aforesaid parameters could, for example, be useful in the attempt to achieve a prototype EEWS that is currently under development in the SW of Spain, where large Table 3. Synthetic view of answers to the main questions of the survey Figure 1. Localization of the studied area and of the main regional seismotectonic features [START_REF] Terrier | Document d'information sur le risque sismique à destination des médias[END_REF]. The historical Pyrenean seismicity [START_REF] Secanell | Probabilistic seismic hazard assessment of the Pyrenean region[END_REF] is represented as well as the synthetic earthquake catalog for an "observation" duration of 10 000 years and M≥5.0.

APPENDIX

Figure 3 .

 3 Figure 3.a clearly underlines that, unless considering a major earthquake characterized by a broad epicentral

  adequate for the kind of seismic signal to record in terms of bandwidth, sample rate, dynamic range, adequate signal/noise ratios, robust enough for the continuous processing of data coming from real-time stations -Centralized analysis by a data-center or notcentralized thanks to a pre-analysis at each

Figure 2 .

 2 Figure 2. Simulated performances of a virtual Pyrenean EEWS based on the EEWS1 SISPyr network: extension of the blind zone as a function of the hypocentral location with an EEWS using a minimum of three records with 2 seconds of the P-wave (Rbz 3,2 ). EEWS1 and EEWS2 seismic networks are also represented .The SISPyr project area is represented by the purple polygon. Topology of the ideal network "EEWS3" is schematically shown in the bottom left-hand corner of the figure.

Figure 3

 3 Figure 3. a. Mean lead-time per magnitude-range corresponding to isoseismal curves I=VI (which corresponds more or less to the threshold of apparition of damages) induced by earthquakes of the synthetic catalog (see Figure2), considering a virtual Pyrenean EEWS based on the EEWS1 network (analysis of the first two seconds of the P-waves on the records of the first triggered stationsranging from one to four: respectively LTI=VI 1,2 , LTI=VI 2,2 , LTI=VI3,2 and LTI=VI 4,2 ). b. Repartition of lead-times when considering at least three triggered stations(LTI=VI 3,2 ).

Figure 4 .

 4 Figure 4. Percentage of earthquakes from the synthetic catalog leading to positive lead-times for areas with given values of intensity, and for different configurations of EEWS (scenarios EEWS1 in red, EEWS2 in purple and EEWS3 in green) -(analysis of the first two seconds of the P-waves on the records of the first three triggered stations).

Figure 5 .

 5 Figure 5. Representation of the study data-set with the distribution of the seismic signals as a function of magnitude and epicentral distance. Top right-hand corner: map of earthquakes included in the catalogue and of corresponding recording seismological stations.

Figure 6 .

 6 Figure 6. A, B, C and D. Empirical correlations linking proxy parameters c, p max , Pd 10 and Pv 10 (means per event) with the local magnitude of Pyrenean earthquakes. Calculations carried out on the signals in the catalogue that satisfy the selection criteria (analysis interval set at 3 s). The straight line obtained by linear regression is shown in black, while the grey dashed lines indicate the confidence interval at 95 % for a new observation.

Table 1 .

 1 Technical requirements for an EEWS seismological network

	y	Analysis length	a		b		R 2	
		1 s	-1.6014	± 0.0689	0.2566	± 0.0198	0.88	0.07
	 c	2 s	-1.5267	± 0.0896	0.2326	± 0.0254	0.80	0.09
		3 s	-1.4870	± 0.1154	0.2198	± 0.0319	0.71	0.10
		1 s	-1.1360	± 0.0611	0.1354	± 0.0178	0.72	0.07
	 p max	2 s	-1.0750	± 0.0577	0.1246	± 0.0168	0.70	0.07
		3 s	-1.1286	± 0.0564	0.1413	± 0.0164	0.76	0.07
		1 s	-8.6609	± 0.1726	0.9279	± 0.0509	0.93	0.21
	P d 10	2 s	-8.7822	± 0.2049	1.0007	± 0.0604	0.92	0.00
		3 s	-8.4182	± 0.1831	0.9049	± 0.0527	0.93	0.19
		1 s	-6.5038	± 0.1803	0.7849	± 0.0531	0.90	0.22
	P v 10	2 s	-6.6589	± 0.1766	0.8536	± 0.0521	0.92	0.00
		3 s	-6.3537	± 0.1582	0.7821	± 0.0455	0.93	0.17

Table 2 .

 2 Summary of the empirical scaling relations established for Pyrenees, with log10 y = a + b.ML, R 2 the coefficient of determination representing the percentage of the y data that can be explained by the log-linear relationship with ML, and  the standard deviation of this relationship.

	Questions	Answers
	Line of business	Industry	Energy	Dam	Transportati on
		55%	18%	18%	9%
	Do you think that your installations are	Yes	No
	exposed to seismic risk?	91%	9%
	In

your opinion, what could be maximal consequences of a major regional earthquake on your installations?

	Without	Minor	Notable	Generalized
	effects	damages	damages	damages
	0%	35%	45%	18%

Does an early warning announcing coming destructive strong motions seem to you useful for your company?

		Yes	No	Without opinion
		73%	9%	18%
	-If yes, from what level of expected effects do you need to	Slight damages	Moderate damages	Major damages
	be alerted?	36%	45%	9%
	Could you precise minimum lead time that you think to be enough to take	<5s	5-10s	10-20s	>20s
	automatic actions to put in security your installations?	0%	9%	18%	64%
	What

kinds of actions seem to you conceivable for your installations with different lead-times?

		-stopping circulation of chemicals fluids (complete stop after few seconds)
		≥10s: Depressurization of (surface) installations and closing of wells
		≥20s: Closure of valves
	Energy	≥20s: Semi-automatic targeted local closing of networks
	Dam	≥10s: Isolating of upstream water by mean of automatic closing of front gates, and
		stopping of production sets (semi-automatic actions)
	Transportation ≥20s: Warning of the traffic control center in order to get under way semi-
		automatic stop of high-speed trains
	In your opinion, could consequences of	Yes	No	Without opinion
	strong earthquake be reduced thanks to the use of a EWS?	82%	9%	9%
	-If yes, for what kind of consequences?	Material damages	Risk of accident over	Financial loss (due loss of data, etc.) to activity stop,	Other
				60%	80%	40%	20%
	If
		≥5s:	-warning of the staff in order to move away from critical
	Industry	installations/processes
		-stopping critical chemical reactions by introduction of security bars (complete
		stop after 1h)

your company adhered to a EWS, at which point would it be important to notice a reduction of consequences in case of earthquake?

	Not important at all		Extremely important
	1	2	3	4	5
	0%	27%	27%	36%	0%
	In

your opinion, what would be the impact of a false alert for your company?

		No impact				Major impact
		1	2	3	4	5
		9%	0%	27%	45%	9%
	In your opinion, what would be the	No impact				Major impact
	impact of a missed alert for your company?	1	2	3	4	5
		0%	9%	18%	45% 18%

Table 4 .

 4 List of stations considered for the early warning scenarios. IDs with an asterisk indicate planed stations whose name and coordinates are likely to change. * Whilst the typical latency time corresponds to the mean value coming from the 48h analysis, the value indicated for station PYLO had been increased by one standard-deviation because of its very low mean value.

	Station			Sensor	Typical latency
		Owner	Lat.	Long.	RT
	ID			type	time (s)
	ARBS	IEA	42.43 1.53 BB	Yes 3.50
	ATE	OMP	43.09 -0.70 BB	Yes 6.59
	AVIN	ICGC	41.85 1.97 ACC	Yes -
	CARA	ICGC	42.71 0.82 BB	Yes -
	CARF	OMP	42.72 2.11 BB	No -
	CAVN	ICGC	41.88 0.75 BB	Yes 3.51
	CBEU	ICGC	42.26 2.68 BB	Yes 3.15
	CBRU	ICGC	42.28 2.18 BB	Yes 3.98
	CCAS	ICGC	41.88 2.90 BB	Yes 2.71
	CELS	ICGC	41.69 2.50 ACC	Yes -
	CEST	ICGC	42.60 1.25 BB	Yes 3.44
	CFON	ICGC	41.76 2.43 BB	Yes 3.84
	CLLI	ICGC	42.48 1.97 BB	Yes 3.17
	CORG	ICGC	42.23 1.32 BB	Yes 3.53
	CORI	ICGC	41.97 2.05 BB	Yes 3.23
	CPAL	ICGC	42.31 3.16 BB	Yes 3.11
	CSOR	ICGC	42.37 1.13 BB	Yes 3.22
	CTRE	ICGC	42.32 0.77 BB	Yes 2.92
	EALK	IGN	43.22 -1.51 BB	Yes -

magnitude earthquakes have occurred in the past (e.g. 1755 Lisbon earthquake, or M7.8, February 1969 earthquake) offshore San Vicente (Romeu et al., 2014). However, this cannot replace the operational testing phase of an EEWS since many aspects such as source/propagation/site/etc. effects, instruments noise level, software robustness, online communication latency, etc., cannot be accounted for by the framework described in this paper.

ACKNOWLEDGEMENTS

This work was supported by the POCTEFA European program Interreg IVA 2007-2013 France-Spain-Andorra, as well as by the ALERTES-RIM project funded by the Spanish Ministry of Economy and Competitiveness. We thank the people that helped in the database creation (M. Calvet from OMP and E. Nus from ICGC), as well as J. Baron for her upstream work on data processing for magnitude scaling-parameters and J. Douglas for his comments that greatly helped us during this study. Special thanks to E. Buforn for her invitation to contribute to this special issue and to J. Douglas for proof reading the final version. Finally, we are very grateful to the three anonymous reviewers whose comments greatly improved the quality of the manuscript.


Édouard Grave

Laurent El Ghaoui

Fast imbalanced binary classification: a moment-based approach

In this paper, we consider the problem of imbalanced binary classification in which the number of negative examples is much larger than the number of positive examples. The two mainstream methods to deal with such problems are to assign different weights to negative and positive points or to subsample points from the negative class. In this paper, we propose a different approach: we represent the negative class by the two first moments of its probability distribution (the mean and the covariance), while still modeling the positive class by individual examples. Therefore, our formulation does not depend on the number of negative examples, making it suitable to highly imbalanced problems and scalable to large datasets. We demonstrate empirically, on a protein classification task and a text classification task, that our approach achieves similar statistical performance than the two mainstream approaches to imbalanced classification problems, while being more computationally efficient.

Introduction

Many real world classification problems are highly imbalanced, meaning that examples drawn from the positive class are significantly outnumbered by negative examples. Consider the problem of building a text classifier to automatically detect news article about renewable energies. Only a small portion of published news articles are about this topic, and the corresponding classification problem is thus highly imbalanced. Detection problems in computer vision, such as finding cars in images, also lead to extremely imbalanced classification problems. Indeed, the negative, or background, class is much more likely than the class of interest.

In the case of highly imbalanced datasets, a naive strategy that classifies all the examples as negative will achieves a very low classification error, since the vast 1 majority of examples are indeed negative. Most classifiers that minimize the classification error thus lead to decision boundaries that are skewed toward the positive class, and thus to a large false negative rate. For most imbalanced problems, this is not acceptable, since the interesting class is the positive class, and having a low false negative rate is important. Another challenge pertaining to imbalanced datasets is the fact that the number of examples from the negative class is often gigantic. This large amount of negative data points are thus the bottleneck of the optimization algorithm.

Related work

Many different approaches have been proposed to deal with imbalanced datasets, and the corresponding literature is too large to be summarized here. We invite the interested reader to look at the extensive review of the subject by [START_REF] He | Learning from imbalanced data[END_REF].

A first class of methods for imbalanced learning is based on sampling: the idea is to sample a balanced training set from the original unbalanced set of examples. Such methods are based on undersampling the negative class (Kubat et al., 1997;[START_REF] Batista | A study of the behavior of several methods for balancing machine learning training data[END_REF], or on (synthetic) oversampling of positive examples [START_REF] Chawla | Smote: Synthetic minority over-sampling technique[END_REF][START_REF] Batista | A study of the behavior of several methods for balancing machine learning training data[END_REF].

A second class of methods, refered to as cost sensitive learning, is based on assigning different misclassification costs to the negative and the positive examples [START_REF] Fan | Adacost: misclassification cost-sensitive boosting[END_REF][START_REF] Drummond | Exploiting the cost (in) sensitivity of decision tree splitting criteria[END_REF][START_REF] Zadrozny | Cost-sensitive learning by costproportionate example weighting[END_REF].

Finally, closely related to our approach, one-class support vector machines were also considered by [START_REF] Raskutti | Extreme re-balancing for svms: a case study[END_REF] in the case of extremely imbalanced datasets and for task of text classification by [START_REF] Manevitz | One-class svms for document classification[END_REF].

Contributions

In this article, we propose a new formulation for solving the problem of imbalanced binary classification. Our method is inspired by the minimax approach to classification proposed by [START_REF] Lanckriet | A robust minimax approach to classification[END_REF]: we model the negative class by moments of its probability distribution instead of examples. Thanks to this approach our formulation does not depend on the number of negative examples, making our method suitable to highly imbalanced problems and making it scalable to large datasets. On the other hand, we propose to model the positive class by individual examples, avoiding the pitfalls mentioned above about false negative rate. More precisely, we make the following contributions:

• we propose a new formulation for the problem of imbalanced binary classification, inspired by the model of [START_REF] Lanckriet | A robust minimax approach to classification[END_REF], where the negative class is represented by its mean and its variance, while the positive class is represented by individual examples (section 2);

• we give a geometric interpretation of our formulation, obtained through Lagrangian duality (section 3);

• we show that our approach is strongly related to one-class support vector machines, as introduced by Schölkopf et al. ( 2001) (section 4);

• we demonstrate that our approach is competitive with the two mainstream approaches to imbalanced classification problems (undersambling and asymmetric cost function) on two classification problems (sections 6 and 7).

Problem formulation

In this section, we propose a new formulation for solving the problem of imbalanced binary classification. We assume that the number of negatively labeled points is much larger than the number of positively labeled points. In order to cope with this large number of negative examples, we propose to model the negative class by its distribution instead of a set of examples. More precisely, we will represent the probability distribution of the negative class by its two first moments, the mean and the covariance. On the other hand, we still represent the positive class by examples.

Let (x i ) i∈{1,...,n} be a set of n positive training examples and let x and Σ be the mean and the covariance of the probability distribution of the negative class. In the following, we will always assume that the covariance matrix Σ is positive definite. This is not a strong assumption, since we can always add a small regularization term λI d to the covariance matrix.

Our goal is to find the hyperplane (w, b) such that all the positive examples are correctly classified:

w ⊤ x i -b ≥ 0, ∀i ∈ {1, ..., n},
while maximizing the probability of correctly classifying examples drawn from the negative class, with respect to all distributions with mean x and covariance Σ:

max w,b inf x∼(x,Σ) P r(w ⊤ x -b ≤ 0), (1)
where x ∼ (x, Σ), refers to the class of probability distributions with mean x and covariance Σ. In other word, our goal is to maximizie the specificity of the separating hyperplane, while correctly classifying all the positive examples.

According to the following lemma from [START_REF] Lanckriet | A robust minimax approach to classification[END_REF], which is a consequence of a theorem by [START_REF] Marshall | Multivariate chebyshev inequalities[END_REF], the condition 1 has a geometric characterization:

Lemma 1. Let x ∈ R d , Σ ∈ R d×d a positive definite matrix, w ∈ R d such that w = 0, b ∈ R and α ∈ [0, 1[. Then, the condition inf x∼(x,Σ) P r x ⊤ w -b ≤ 0 ≥ α, holds if and only if b -x⊤ w ≥ κ(α) √ w ⊤ Σw, where κ(α) = α 1-α .
Using this result, the optimal hyperplane we are looking for is thus the optimal solution of the following optimization problem:

max α,w,b α s.t. b -x⊤ w ≥ κ(α) √ w ⊤ Σw and x ⊤ i w -b ≥ 0. Since the function κ : α → α 1-α is increasing on [0, 1[, this problem is equivalent to the program: max κ,w,b κ s.t. b -x⊤ w ≥ κ √ w ⊤ Σw and x ⊤ i w -b ≥ 0.
We can further eliminate b from this program, finally obtaining:

max κ,w κ s.t. x ⊤ i w -x⊤ w ≥ κ √ w ⊤ Σw. ( 2
)
This problem is a convex program [START_REF] Boyd | Convex optimization[END_REF], which has an interesting geometric interpretation that will discuss in the next section. We will also show that this formulation is strongly related to one-class support vector machines [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF].

Dual problem and geometric interpretation

In this section, we derive the dual of the problem defined in equation 2 through Lagragian duality. This will allow us to give a nice geometric interpretation of our approach. Assuming that Σ is invertible, we can perform the change of variable u = Σ 1/2 w. Then, our problem becomes:

max κ,u κ s.t. z ⊤ i u -z⊤ u ≥ κ √ u ⊤ u,
where z i = Σ -1/2 x i and z = Σ -1/2 x. Exploiting the homogeneity in variable u, we can impose the constraint u 2 = 1, leading to the program

max κ,u κ s.t. z ⊤ i u -z⊤ u ≥ κ and u 2 = 1.
This problem is equivalent to

max κ,u κ s.t. z ⊤ i u -z⊤ u ≥ κ and u 2 ≤ 1.
Then, introducing dual variables p i , it can be expressed as

max κ,u min p i ≥0 κ + i p i (z i -z) ⊤ u -κ s.t. u 2 ≤ 1.
Using Sion's minimax theorem [START_REF] Sion | General Minimax Theorems[END_REF], we can invert the min and the max operators, obtaining:

min p i ≥0 max κ,u κ + i p i (z i -z) ⊤ u -κ s.t. u 2 ≤ 1.
Then, if i p i = 1, the max over κ is equal to +∞. The previous program is thus equivalent to

min p i ≥0 max u i p i (z i -z) ⊤ u s.t. u 2 ≤ 1 and i p i = 1.
Finally, using the definition of the dual norm, and the fact that the dual norm of the ℓ 2 -norm is itself, we get

min p i ≥0 i p i z i -z 2 s.t. i p i = 1.
We can now replace z i and z by their expressions, to obtain the dual problem

min p i ≥0 i p i x i -x Σ -1 s.t. i p i = 1,
where • Σ -1 is the norm defined by:

x Σ -1 = √ x ⊤ Σ -1 x.
Thus, the dual problem correspond to finding the orthogonal projection of the mean vector x of the negative class onto the convex hull of the positive examples, for the scalar product defined by the inverse covariance matrix Σ -1 . Given an optimal solution p * of the dual problem, the correponding primal solution is equal to

w * = Σ -1 (Xp * -x).
Thus, the corresponding separating hyperplane is orthogonal to the difference between x and its projection on the convex hull of the positive points x i , for the scalar product defined by the inverse covariance matrix Σ -1 . This geometric interpretation is illustrated in Figure 1.

Relation to one-class support vector machines

In this section, we show that our approach is strongly related to the one-class support vector machine formulation introduced by [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF]. First, let us restate our formulation for solving the problem of imbalanced binary classification:

max κ,w κ s.t. x ⊤ i w -x⊤ w ≥ κ √ w ⊤ Σw.
Exploiting the homogeneity in w again, we now impose that κ √ w ⊤ Σw = 1, leading to the formulation

max w,κ κ s.t. (x i -x) ⊤ w ≥ 1 and κ = 1 √ w ⊤ Σw .
We can now eliminate the variable κ, and since the function x → 1/ √ x is decreasing on R * + , we obtain the equivalent program:

min w w ⊤ Σw s.t. (x i -x) ⊤ w ≥ 1.
This formulation is very close to the formulation of hard margin, one-class support vector machines, with two important differences: first, we do not minimize the ℓ 2 -norm of w, but the Mahalanobis norm corresponding to the covariance matrix of the negative class distribution. Second, we do not try to separate the positively labeled points from the origin, but from the mean of the negative class distribution.

Similarly to support vector machines, the constraint that all positive examples should be correctly classified might be unrealistic in practice. We thus propose to relax these constraints by introducing the penalized slack variables ξ i ≥ 0. We then obtain the following convex program:

min w 1 2 w ⊤ Σw + C n i=1 ξ i (3) s.t (x i -x) ⊤ w ≥ 1 -ξ i , ξ i ≥ 0,
where C ∈ R + is a tradeoff parameter. This formulation is the one that we will use in practice, and will be refered to as moment-based imbalanced binary classifier, or MIBC.

Again, this formulation is strongly related to the formulation of soft margin, one-class support vector machines. Indeed, by making the change of variables u = Σ 1/2 w, and by introducing vectors z i such that z i = Σ -1/2 (x ix), we obtain the formulation

min w 1 2 u ⊤ u + C n i=1 ξ i s.t. z ⊤ i u ≥ 1 -ξ i , ξ i ≥ 0
which is exactly the formulation of one-class support vector machines.

Solving the optimization problem

In this section, we discuss how to optimize our proposed formulation for imbalanced binary classification. As said earlier, one major difference between our approach and (one-class) support vector machines is the fact that we minimize the Mahalanobis norm corresponding to the covariance matrix Σ of the negative class distribution. Thus, the dimension of the problem will be an important factor in the choice of the optimization algorithm for our method. We will thus discuss optimization algorithms for small dimensional problems first (problems where the dimension is smaller than the number of positive data points), before moving to high dimensional problems. But first, let us derive the dual of the problem introduced in equation 3, as it will be usefull in our discussion.

Dual problem

One of the classical ways to compute the solution of a support vector machine is to solve the dual problem [START_REF] Platt | Sequential minimal optimization: A fast algorithm for training support vector machines[END_REF]. We will follow a similar approach. Thus, we now compute the dual problem of our our moment-based imbalanced binary classifier, as defined in equation 3. Let us define the vectors xi = x ix.

Introducing Lagrange multipliers α ∈ R n and ν ∈ R n , such that α i ≥ 0 and ν i ≥ 0, the Lagrangian corresponding to the problem 3 is equal to

L(w, ξ, α, ν) = 1 2 w ⊤ Σw + C n i=1 ξ i - i α i x⊤ i w -1 + ξ i -ν ⊤ ξ.
To find the dual function, we minimize the Lagrangian over w and ξ. Minimizing over ξ, we find that the dual function is equal to -∞, unless C -α i -ν i = 0, in which case, we are left with

L(w, ξ, α, ν) = 1 2 w ⊤ Σw - i α i x⊤ i w -1 .
Minimizing over w, we then obtain that

w = Σ -1 Xα, where X = [x 1 , ..., xn ] = [x 1 -x, ..., x n -x].
Replacing w by its optimal value and combining the constraints ν i ≥ 0 with C -α i -ν i = 0, we finally obtain the dual problem

max 0≤α≤C - 1 2 α ⊤ X⊤ Σ -1 Xα + 1 ⊤ α. (4)
Unsurprisingly, this dual is very similar to the dual of support vector machines. Thus, most algorithms that were proposed to solve the dual of SVMs could be used to solve this dual.

Small dimensional problems

Similarly to support vector machines, the dual problem obtained in equation 4 is a quadratic program of dimension n, where n is the number of positive points. We first need to compute the inverse of the covariance matrix Σ. This operation has a complexity of O(d 3 ), where d is the dimension of the data. Then, finding the optimal solution of the dual problem has a complexity of O(n 3 ). Thus, the overall complexity is O(d 3 + n 3 ), and this approach is only applicable in the case of small dimensional problem. In the following, we will assume that the covariance matrix Σ has some structure that can be exploited to solve high dimensional problems.

High dimensional problems: factor models

In this section, we describe how to speed up the computation in the case of a factor model, that is when Σ = D + FF ⊤ , where D is a positive definite diagonal matrix and F ∈ R d×k with k ≪ d. We will assume, without loss of generality, that D = λI d . Replacing Σ by its new expression, the dual problem obtained in equation 4 then become:

max 0≤α≤C - 1 2 α ⊤ X⊤ (λI d + FF ⊤ ) -1 Xα + α ⊤ 1.
Using the Woodbury matrix identity, we then obtain the equivalent problem

max 0≤α≤C - 1 2λ α ⊤ X⊤ I d -F(λI k + F ⊤ F) -1 F ⊤ Xα + α ⊤ 1.
We thus have replaced the inverse of a d × d matrix by the inverse of a k × k matrix. The most expensive step to compute the matrix appearing in the dual cost function is to compute the matrix product F ⊤ X. The complexity of this operation is O(kdn). Then, solving the dual has a complexity of O(n 3 ).

Another potential approach to speed up the optimization algorithm for high dimensional settings is to directly estimate a sparse inverse covariance matrix [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF]d'Aspremont et al., 2008). Recently, [START_REF] Rolfs | Iterative thresholding algorithm for sparse inverse covariance estimation[END_REF] and [START_REF] Hsieh | Big & quic: Sparse inverse covariance estimation for a million variables[END_REF] have proposed efficient algorithms to solve this problem, making it possible to estimate sparse inverse covariance matrices in high dimensional spaces.

Finally, it should be noted that it is often the case that several problems with the same covariance matrix have to be solved, for example when doing a grid search for parameter selection. Thus, the computation or the estimation of the inverse of the covariance matrix (or the computation of a factor model) is often amortized over the resolution of several problems.

Small scale experiments

In this section, we present experiments performed on small dimensional datasets. We compare our moment-based imbalanced binary classifier (MIBC) with two standard strategies for imbalanced classification problems, using SVMs: undersambling the negative class and using different costs for the negative and positive examples. We will use the liblinear [START_REF] Fan | Liblinear: A library for large linear classification[END_REF]) implementation of linear support vector machines, and an implementation of our approach using Mosek1 . It is thus important to note that our implementation is based on a general purpose quadratic programming solver. A custom implementation, for example based on the SMO algorithm, should greatly improve the efficiency of our moment-based imbalanced binary classifier.

Dataset 1: Basic statistics about the different datasets.

Datasets

We evaluated the different methods on four datasets introduced by [START_REF] Radivojac | Classification and knowledge discovery in protein databases[END_REF], which are publicly available2 . These datasets correspond to protein classification problems, such as predicting protein phosphorylation sites (PHOST, PHOSS, PHOSY) or predicting binding regions (CAM). Following the approach proposed by [START_REF] Radivojac | Classification and knowledge discovery in protein databases[END_REF]

Discussion

We report areas under the ROC curve for the four datasets in Table 2, computational times in Table 3 and ROC curves for two datasets in Figure 2 (PHOST and PHOSY). We performed a paired samples t-test to determine if our results are statistically significant.

First, we observe that our moment-based imbalanced binary classifier always outperforms the undersampling approach, while performing at least as well as the cost sensitive method. Second, the two datasets on which our method outperforms the asymmetric cost function SVM (PHOST and PHOSY) correspond to the highest ratio of number of negative to positive examples (64 and 37 respectively). This seems to indicate that our method is particularly adapted to highly imbalanced datasets. Finally, our method is computationally more efficient, leading to speedup between 1.4 and 4.8 over cost-sensitive SVM, while obtaining as good or even better statistical performances. We remind our reader that we implemented our method using Python and Mosek, and it is thus certainly possible to get much better performances.

Application to text classification

In this section, we report experiments performed on the task of text classification. We will follow the same methodology as described in the section 6.2. Since bagof-words representations of textual documents live in high dimensional spaces, we propose to replace the full covariance matrix of the negative class by its diagonal.

Dataset

We use the REUTERS RCV1 dataset, introduced by [START_REF] Lewis | Rcv1: A new benchmark collection for text categorization research[END_REF], which is a classical test bed for text classification methods. Each document of the corpus is tagged with respect to three different category sets: topics, industries and regions.

We consider classification problems that consist in classifying documents that are labeled with a given topic label v.s. the rest of the documents. There are 104 different topics, and we will thus consider only a subset of the 104 possible classification tasks. Since we want to focus on highly imbalanced classification problems, we set the ratio of negative examples to positive examples to 1, 000.

Discussion

We report areas under the ROC curve in Table 4, computational times in Table 5 and ROC curves in Figure 3. We performed a paired samples t-test to determine if our results are statistically significant.

We observe that our moment-based imbalanced binary classifier achieves similar statistical performances than the cost-sensitive method, while generally outperforming the undersampling approach. In particular, we observe that on the classification problem TOPIC 25 V.S. REST our method achieves a much lower false

Figure 1 :

 1 Figure 1: Geometric interpretation of our proposed formulation for imbalanced binary classification.

Figure 2 :

 2 Figure 2: ROC curves, averaged over twenty experiments, on the PHOST and PHOSY datasets.

Figure 3 :

 3 Figure 3: ROC curves, averaged over ten experiments, on the REUTERS RCV1 dataset.

Table

		# positive # negative ratio
	PHOSS	613	10, 798	17
	PHOST	140	9, 051	64
	PHOSY	136	5, 103	37
	CAM	942	17, 974	19

Table 2 :

 2 , we keep the 150 features which are the most correlated to the class labels. The ratio of negative examples to positive ones varies from 17.6 to 64.6 on the different datasets. Basic statistics about those are given in Table1.For each dataset, we use 50% of the examples as training set, 20% as validation set and 30% as test set. For all methods, we chose C in the set {10 5 , ..., 10 -4 }. When undersampling the negative class, we keep as many negative examples as positive examples. For the asymmetric cost function method, we consider the following ratios between the weights for the positive examples and negative examples: {1.0, r/4, r/2, r, 2r}, where r is the ratio of number of negative examples to the number positive ones. We replicated the experiments over twenty random splits of the data. Areas under the ROC curve (with confidence intervals), averaged over twenty experiments, on the protein classification tasks. † indicates that our method is significantly better than the two others (with p-value p < 0.01).

	6.2. Methodology

Table 3 :

 3 Computational times, in milliseconds, required to solve one protein classification problem, averaged over twenty experiments.

Table 4 :

 4 Areas under the ROC curve (with confidence intervals), averaged over ten experiments, on the text classification tasks. Differences between our momentbased imbalanced binary classifier and the subsampling method are statistically significant (with p-value p < 0.01).

	Topic This work Cost-sensitive Sampling
	2 9 25 33 59 84	89.7 ± 1.0 96.1 ± 0.7 95.1 ± 0.8 96.0 ± 0.4 96.1 ± 0.4 96.9 ± 0.8	89.9 ± 1.4 96.3 ± 0.8 94.3 ± 1.6 95.7 ± 0.6 95.9 ± 1.4 96.4 ± 1.5	87.7 ± 1.2 94.1 ± 1.3 93.7 ± 1.2 93.9 ± 0.7 95.0 ± 0.6 96.3 ± 0.9

www.mosek.com

www.informatics.indiana.edu/predrag/publications.htm

Acknowledgments

Édouard Grave is supported by a grant from Inria (Associated-team STATWEB).

positive rate, for a true positive rate equal to 1.0 (See Figure 3). Finally, our approach to imbalanced classification is much more computationally efficient than a SVM with asymmetric costs, leading to speed-up between 21 and 36.

Conclusion

In this paper, we introduced a new approach to imbalanced classification problems, refered to as moment-based imbalanced binary classifier (MIBC). We proposed to cope with the large number of negative examples by modeling the negative class by the mean and the variance of the corresponding probability distribution. On the other hand, we still model the positive class by individual examples. We show that our formulation is strongly related to one-class support vector machines. We empirically demonstrated that our approach is competitive with other methods for imbalanced classification, based on subsampling and cost sensitive learning. In particular, we demonstrated that our approach leads to much lower computational times than cost-sensitive learning, while obtaining as good, or even better, classification performances. On the other hand, our method achieves better statistical performances than subsampling.

As future work, we would like to apply our technique to detection problems in computer vision, such as finding cars or pedestrians in images. We would also like to implement efficient optimization algorithms for our approach, in order to fully exploit the much reduced computational complexity of our approach, compared to cost-sensitive support vector machines.