K Samyn
email: k.samyn@brgm.fr

F Mathieu

A Bitri

A Nachbaur

L Closset

ÔØ Å ÒÙ× Ö ÔØ Integrated geophysical approach in assessing karst presence and sinkhole susceptibility along flood-protection dykes of

Keywords: River dyke, seismic surface waves, electrical resistivity tomography, sinkhole, karst

Non-invasive geophysical methods are often used for detecting near-surface defects and monitoring seepage in river dykes or dams. Between 2008 and 2009, a series of geophysical experiments were conducted to detect karst features below the dykes of the Loire River, Orléans, France. Multi-channel analysis of seismic surface waves (MASW) was used to obtain the shear wave velocity (Vs) profile of the subsurface below the dykes. As an effective approach for investigating the structure of the dykes, multi-channel resistivity surveys were also used to evaluate the electrical properties of material inside and under the dykes. This study discusses the exploration strategy and results for several sections of the dyke system.

Based on the experimental data, Vs contours are used to geometrically and qualitatively describe the subsurface under the dykes and identify areas of mechanical weakness corresponding to karst features, while resistivity contours allow distinguishing interactions between karst and the dykes through the identification of areas of downward material movement. Known and unknown anomalies are identified. A practical approach that combines seismic and electric resistivity results is proposed to assess karst presence and

Please cite this article as: Samyn, K., Mathieu, F., Bitri, A., Nachbaur, A., Closset, L., Integrated geophysical approach in assessing karst presence and sinkhole susceptibility along flood-protection dykes of the Loire River, Orléans, France, Engineering Geology (2014), doi: 10.1016/j.enggeo.2014.10.013 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. Introduction Shallow karst constitute a serious hazard for the public and existing construction due to the risk of sinkhole development and subsidence, which can cause damage and increase the costs of urban development. Detecting karst features, voids, surficial dissolution, fissuring, alteration and unconsolidated material by geophysical techniques is thus a major challenge for many infrastructure managers. As a matter of fact, reliable characterization of the dyke material are required for the verification of its stability. The importance of an observation of dyke stability became obvious in the light of major flood events, such as the Elbe flood of 2002 in central Europe. Survey campaigns on dykes often focus on geotechnical parameters obtained in situ by standard penetration and cone penetration tests, or in the laboratory based on borehole samples. A disadvantage of these tests is their selective character. An extensive number of tests are required to cover a dyke sufficiently. On the other hand, most nondestructive geophysical methods are able to survey large areas effectively. The dykes that were the subject of this study are located approximately 100 km south of Paris, on the banks of the Loire River, and are built on the highly karstified calcareous Beauce Formation of Aquitainian age substratum. As evidence of this, the Loiret springs are amongst the most important karst springs in France, and sinkholes regularly occur in their vicinity. These cavities were also the cause of unforeseen expenses during the construction of a tramway line and a bridge over the Loire River. Geophysical measurements conducted on the Loire River dykes are integrated into a general approach that aims to assess karst hazards. Karst systems are characterised by complex geomorphology and hydrology [START_REF] Ford | Karst Hydrogeology and Geomorphology[END_REF].

Voids can take the form of large cavities, super-imposed conduits and networks of open fractures or micrometric fissures. Karstification is commonly revealed by surficial features; in the uppermost layer of karstified zones, dissolution by meteoric water can enlarge fissures, thereby creating a zone of rapid infiltration filled with weathered, fractured and

A C C E P T E D

M A N U S C R I P T unconsolidated material, or epikarst. Epikarst development can cause sinkhole phenomena draining surface materials [START_REF] Palmer | Origin and morphology of limestone caves[END_REF][START_REF] Gunn | Encyclopedia of Caves and Karst Science[END_REF]. Some dangerous segments of the dykes need to be rebuilt, and others need to be strengthened. Because the nature and structure of the materials under the dykes vary from place to place, decisions about their treatment should first be supported by an understanding of the in situ subsurface conditions. The application of near-surface geophysics to river dyke characterisation has been studied in the past and has shown great promise [START_REF] Bitri | River dykes investigation and effects of the topography on seismic surface waves propagation[END_REF][START_REF] Chen | Integrated geophysical techniques in detecting hidden dangers in river embankments[END_REF][START_REF] Karl | Geotechnical characterization of a river dyke by surface waves[END_REF]. A

program to evaluate the condition of river dykes, supported by the local Regional Environment Department (DREAL), was executed from 2006 to 2011. The French Geological Survey (BRGM) and the DREAL funded a singular project to evaluate the probability of karst presence and the sinkhole susceptibility using a methodology based on a geophysical survey, under the dykes of the Loire River. The purpose of the project was therefore quite innovative as geophysical methodologies are usually used to characterize known hazardous events rather than try to predict them. Geophysical surveys covered a total length of more than 90 km of dykes with expected high probability of karst presence [START_REF] Donsimoni | Détermination de la présence de karst sous les levées domaniales du bassin de la Loire et réduction du risque d'effondrement de la digue lié à la présence de ces conduits souterrains naturels. Tranche 1: analyse géologique et hydrogéologique[END_REF]. The complexity of karst features produces multiple and timevariable geophysical signatures. Shear-wave velocity and resistivity are known to be very good indicators of the degree of weathering and fissuring of a draining or saturated epikarst, and electrical behaviour may be enhanced by the presence of shallow groundwater and strong water circulation. Providing that anthropogenic interference allows their use, electrical imagery is a proven method for indirectly detecting karst systems [START_REF] Guerin | Shallow karst exploration using MT-VLF and DC resistivity methods[END_REF][START_REF] Ioannis | Accurate subsurface characterization for highway applications using resistivity inversion methods[END_REF][START_REF] Kaufmann | Les effondrements karstiques du Tournaisis : genèse, évolution, localisation, prévention[END_REF][START_REF] Valois | Karst morphologies identified with geophysics around Saulges caves (Mayenne, France)[END_REF][START_REF] Valois | Assessment of doline geometry using geophysics on the Quercy plateau karst (South France)[END_REF]. However, the use of a combination of geophysical tools is always preferable for obtaining a better-constrained model of near-surface features [START_REF] Chen | Integrated geophysical techniques in detecting hidden dangers in river embankments[END_REF]. For example, the exploration of a karst network in the Swiss Jura was conducted by combining GPR, which provides a precise image

A C C E P T E D M A N U S C R I P T
of most surficial voids and structures, with gravimetry, which is sensitive to the presence of deeper and more extensive heterogeneities [START_REF] Beres | Integration of ground-penetrating radar and microgravimetric methods to map shallow caves[END_REF]. To fit the objective of the project, a survey coupling electrical resistivity tomography (ERT) and multichannel analysis of surface waves (MASW), suitable for providing evidence of areas of mechanical weakness [START_REF] Debeglia | Karst investigation using microgravity and MASW: application to Orléans, France[END_REF][START_REF] Ganji | Detection of underground obstacles by SASW method-numerical aspects[END_REF][START_REF] Matthews | The use of surface waves in the determination of ground stiffness profiles[END_REF][START_REF] Leparoux | Underground cavity detec-tions: a new method based on seismic Rayleigh waves[END_REF][START_REF] Park | Ground roll as a tool to image near-surface anomaly[END_REF][START_REF] Shtivelman | Surface wave section as a tool for imaging subsur-face inhomogeneities[END_REF], was therefore undertaken. Several anomalous zones indicating karst activities were detected with both geophysical methods used. Indices derived from a coupled empirical analyse of the geophysical parameters then allowed the assessment of a karst presence probability and a sinkhole susceptibility along the investigated dykes.

Supplementary geotechnical and observational data acquired from databases were used for calibration and as a test of the methodology.

Experimental site

The geophysical surveys were carried out at the top of the Loire River flood-protection dykes, near the city of Orléans, France. [START_REF] Donsimoni | Détermination de la présence de karst sous les levées domaniales du bassin de la Loire et réduction du risque d'effondrement de la digue lié à la présence de ces conduits souterrains naturels. Tranche 1: analyse géologique et hydrogéologique[END_REF] estimated the probability of karst presence over the entire Loire River and especially in this area (Fig. 1). In the light of this study, 90 km of dykes divided into 10 profiles with high probability of karst presence were then selected and investigated with geophysical surveys (Fig. 2). The geophysical data acquisition was realized in 2008 and spread over a period of one year. 50 destructive geotechnical boreholes were also available on the investigated segments at the top of the dykes, useful for the calibration of the geophysical responses.

These boreholes are described in further sections and were, by courtesy, provided by the local French roads and bridges laboratory (LRPC, Blois). The synthesis of boreholes data allowed to highlight the local lithology classification consisting in the upper sediments made of siliceous and argillaceous alluvia, 5 to 10 m thick, overlying the calcareous Beauce Formation (Aquitainian age), in which significant karst dissolution features occur. Karstification has resulted in many surface depressions (pits, sinkholes) and irregularities in both the upper surface of the Beauce Formation and at the ground

A C C E P T E D M A N U S C R I P T
surface (Fig. 3). These conduits, generally around 2 to 3 m wide by 1 to 3 m thick, plunge quite rapidly to depths of 15-20 m or more.

Methodology

Exploration strategy

Electrical resistivity tomography (ERT) has become widely used in environmental and engineering applications. The method is popular in the investigation of groundwater (e.g., [START_REF] Benson | Mapping groundwater contamination using DC resistivity and VLF geophysical methods-A case study[END_REF][START_REF] Khair | A model study of the effect of salination on groundwater resistivity[END_REF][START_REF] Saksa | Case study: mapping of groundwater conditions at Olkiluoto site with electrical and electromagnetic boreholes[END_REF][START_REF] Sandberg | Examples of resolution improvement in geoelectrical boreholes applied to groundwater Investigations[END_REF] as well as monitoring leaks in dyke dams (e.g., [START_REF] Johansson | Seepage monitoring in an earth embankment dam by repeated resistivity measurements[END_REF]Sjödahl et al., 2003). A multi-channel data acquisition system is suitable for obtaining the variation of resistivity with depth, using different electrode configurations along a line. The ERT approach was used to delineate the resistivity profile of the dyke bodies, with the aim of characterising their structure and detecting near-surface features.

Dispersive phase velocities of surface waves can be inverted to obtain shear-wave velocities (Vs) (e.g., [START_REF] Aki | Quantitative Seismology[END_REF]. The Vs profile of near-surface materials is one of the key parameters required in engineering. Multi-channel analysis of surface waves (MASW) [START_REF] Park | Multichannel analysis of surface-waves[END_REF][START_REF] Xia | Estimation of nearsurface shear-wave velocity by inversion of Rayleigh wave[END_REF] has shown great promise in detecting shallow voids and tunnels, mapping bedrock surfaces, and delineating fracture systems in rock [START_REF] Miller | Multichannel analysis of surface waves to map bedrock[END_REF][START_REF] Xia | Construction of vertical section of near-surface shearwave velocity from ground roll[END_REF][START_REF] Xia | Delineation of a collapse feature in a noisy environment using a multichannel surface wave technique[END_REF]. The Vs of materials is an indicator of the seismic stiffness of the material.

Current improvements in data acquisition, processing, and inversion (e.g., Xia et al., 2002aXia et al., , 2002bXia et al., , 2002c[START_REF] Xia | Inversion of high frequency surface-waves with fundamental and higher modes[END_REF][START_REF] Xia | Delineation of a collapse feature in a noisy environment using a multichannel surface wave technique[END_REF][START_REF] Xia | Resolution of high-frequency Rayleigh-wave data[END_REF] provide confidence in the assessment of Vs. In this study, the MASW approach was used to delineate the Vs profile of the dykes, with the aim of detecting mechanical weakness in the deep limestone caused by karstification.

A C C E P T E D M A N U S C R I P T

ERT acquisition and processing

2D resistivity surveys were carried out using a CC resistivity meter, the OhmMapper (Geometrics Inc., CA, USA), which uses line antennas as its capacitive sensors to acquire apparent resistivity data. The only array implemented by the OhmMapper is a line electrode dipole-dipole array which can properly image disturbances associated to karst features due to a high spatial resolution [START_REF] Dahlin | A numerical comparison of 2D resisitivity imaging with 10 electrode arrays[END_REF]. However other arrays could also provide reliable results for this type of application but no towable acquisition device with capacitive sensors currently exists with other electrode arrays to provide an acceptable cost effectiveness ratio for such a large investigation area. The system is a capacitively coupled resistivity meter that measures the electrical properties of rock and soil without the cumbersome ground stakes used in traditional resistivity surveys. A simple coaxial cable array with transmitter and receiver sections was pulled along the ground by a four-wheel drive. Data collection is many times faster with this system than with systems using conventional DC resistivity, and the system allowed an acquisition rate of approximately 8 km/day. The emitter on the back is connected to the receiver by an electrically insulating rope. The value of the intensity of the current injected by the transmitter is transmitted to the receivers via Bluetooth. Each receiver measures an electric potential whose value is sent to the recorder, together with the amperage.

The recorder stores data at a frequency set by the user. For each position of the transmitter, the electrical signal was simultaneously recorded by 4 receiver locations. Dipole receivers are connected to each other and transmit the measured potential to the recorder. The distance between transmitter and receiver, as well as the length of the transmitter and receiver dipoles, determines the investigation depth. In this experiment, the line antennas (electrodes) with a nominated length of 5 m were used for both the transmitter and receivers. Therefore, the dipole length of the line electrode array was 10 m. Four receivers were used in the measurement with a rope length of 5 m. The corresponding available dipole distances were

A C C E P T E D M A N U S C R I P T
then 15 m, 20 m, 25 m and 30 m. The apparent resistivity data were first filtered, removing data with standard deviations higher than 3%. A visual inspection was then made to remove outlier data before being smoothed using X2IPI [START_REF] Bobatchev | X2IPI software[END_REF] to reduce geological noise (Fig. 4a andb). The apparent resistivities were then inverted to a default depth of 12 m using Res2Dinv [START_REF] Loke | A comparison of Gauss-Newton and quasi-Newton methods in resistivity imaging inversion[END_REF][START_REF] Loke | A comparison of smooth and blocky inversion methods in 2D electrical imaging surveys[END_REF] to obtain inverted resistivity geoelectric cross-sections. Each model consisted of 9 horizontal layers. Inversions with five iterations, robust constraint and a vertical to horizontal flatness filter ratio of 0.3, yielded satisfactory convergence, as is indicated by an RMS error ranging between 0.86% and 4.8% between measured and modeled apparent resistivity values for all surveys. To avoid an over-interpretation of the results, inverted resistivity sections were systematically truncated to a depth properly constrained by data of 8 m which can therefore be regarded as a reliable depth for interpretation.

MASW acquisition and processing

Seismic data were acquired following a roll-along procedure. To increase the speed and efficiency of MASW data recording and to thereby keep acquisition costs down, a new type of multichannel seismic cable was designed and manufactured to reach an acquisition rate of approximately 3 km/day. It consists of 24 geophones, at fixed 2-m intervals. Each geophone is attached to a single self-orienting, gimbal-mounted, vertical geophone with a 10 Hz resonance frequency. To help ensure proper coupling, each gimbal geophone is housed in a heavy casing (~1 kg). To damp the motion of the sensor around its rotational axis, the inside of the casing is filled with viscous oil. The seismic cable is towed behind a vehicle. A 24channel Geometrics Stratavizor seismograph was used to record the impact from a weightdrop electronic seismic source. The source-to-nearest-receiver offset was 4 m. Therefore, the total length of the seismic line was 50 m. The source stations were separated by 10 m

A C C E P T E D M A N U S C R I P T
increments along the segments of the dykes. Various tests (van der [START_REF] Van Der Veen | Land streamer for shallow seismic data acquisition: evaluation of gimbal mounted geophones[END_REF] applying these techniques have demonstrated that high-quality data, comparable to planted spike geophone data, can be acquired with gimbal geophones (Fig. 4c). Each 24-trace shot gather was analysed with Surfseis v2.0 [START_REF] Park | Multichannel analysis of surface-waves[END_REF], facilitating the use of MASW with continuous profiling techniques. The first step consists of determining the dispersion curve (i.e., the phase velocity variation versus frequency) for each shot. The curve is calculated by using the slant-stack method in a common-shot gather, followed by a 1D Fourier transform over the intercept time (McMechan and Yeldin, 1981;[START_REF] Moktar | Seismic velocity and Q model for the shallow structure of the Arabian shield from short period Rayleigh waves[END_REF], Park et al., 1998).

The dispersion curve is directly obtained by picking the maximum values of the fundamental propagation mode of the obtained spectrum (Fig. 4d). The frequency range of 8 to 60 Hz was globally obtained for the dispersion curves of the whole seismic survey. The phase velocity and the thickness of layers for each dispersion curve were inverted using a 10 layers model to obtain a 1D Vs (shear wave velocity) profile with depth. Different approaches have been proposed for the inversion procedure. The more widely used approach is the linearised iterative least-squares technique [START_REF] Gabriels | In-situ measurements of shear wave velocity in sediments with higher mode Rayleigh waves[END_REF][START_REF] Moktar | Seismic velocity and Q model for the shallow structure of the Arabian shield from short period Rayleigh waves[END_REF][START_REF] Socco | Surface-wave method for near-surface characterization: a tutorial[END_REF][START_REF] Xia | Estimation of nearsurface shear-wave velocity by inversion of Rayleigh wave[END_REF]. The linearised iterative least-squares technique used here is adapted from [START_REF] Hermann | Computer Programs in Seismology[END_REF]. The inversion was stopped when a satisfactory convergence was reached, as is indicated by a maximum accepted RMS error in phase velocity between calculated and measured dispersion curves of 5 m/s. By placing each 1-D Vs profile at a surface location corresponding to the middle of the receiver line, a 2-D surface versus depth Vs map is constructed through a kriging interpolation scheme.

4.

Results of ERT

Thirty-five destructive geotechnical boreholes were performed on the investigated segments at the top of the dykes, allowing calibration of the electrical responses. Information on these

A C C E P T E D M A N U S C R I P T
surveys was, by courtesy, provided by the local French roads and bridges laboratory (LRPC, Blois). The synthesis of these calibrations allowed development of a lithology classification based on the inverted resistivity. Resistivities below 75 Ohm.m were found to be in good overall agreement with a clayey silt lithology (A) encountered by boreholes (Fig. 5a).

Resistivities ranging from 75 to 200 Ohm.m were found to be in good overall agreement with a clayey sand lithology (SA) with a slightly variable clay fraction (SAa+ or SAs+) (Fig. 5b).

Finally, resistivities over 200 Ohm.m were found to be in good overall agreement with clean sand and gravel lithology, corresponding to the Loire River deposits (SG) (Fig. 5c).

Disorders within the dyke bodies have also been detected by sudden and large changes in resistivity. These resistivity anomalies are generally limited to a few hundred metres in size, the dimensions of which are consistent with historical breaches, whose average length was evaluated at 190 m [START_REF] Lino | Méthodologie de diagnostic des digues appliquée aux levées de la Loire moyenne[END_REF]. These resistivity anomalies are highlighted either by a conductor body (A or SA) in a resistant bank (SG) or, conversely, by a resistant body (SG) in a conductor bank (A or SA), and they have been compared to historical disorders available in the French National database of underground cavities (BD cavités).

According to this comparison, the resistivity anomalies in the dyke body can have different origins, such as former breaches, sinkhole filling or former reinforcement works on the dykes. Figure 6 shows examples of such detected disorders in the dyke body.

Disorders affecting the base of the dykes (e.g., the alluvial cover) are also highlighted by sudden and large changes in resistivity. The extension of these resistivity anomalies is generally less than in the dyke bodies, roughly a dozen metres to a hundred metres long.

These disorders are highlighted not only by a sudden change in resistivity but also by a turning-down pattern of the interface characterising the base of the dyke body. These anomalies are likely to correspond to areas of potential sediment flow toward a karst network.

Figure 7 shows examples of this type of anomaly.

A C C E P T E D M A N U S C R I P T

Results of MASW

Calibration of the MASW was made more difficult by the small number available destructive geotechnical boreholes (15) that reached the underlying karstified limestone. However, synthesis of these calibrations with available boreholes allowed development of a lithology classification based on the inverted Vs (Fig. 8). Vs ranging from 100 to 300 m/s was found to be in good overall agreement with the Loire River overburden deposits (Se) encountered by boreholes. Vs ranging from 300 to 450 m/s were found to be in good overall agreement with a strongly karstified marly limestone (CA). Vs ranging from 450 to 600 m/s was found to be in good overall agreement with a more competent limestone-marl (C). Finally, Vs greater than 600 m/s was found to be in good overall agreement with the competent limestone-marl bedrock (S).

This classification allowed karstified areas at the top of the marly limestone formations to be detected by sudden and significant decreases in Vs. These Vs anomalies are generally limited in dimension to a few hundred metres, revealing the presence of karst networks instead of isolated karst channels, and are in good overall agreement with historical karst-related surface disorders of the BD cavités. The associated low-velocity anomalies may be due to weathering of the medium (CA) because of karstification. Figure 9 shows examples of this type of anomaly.

Assessing probability of karst presence

As Values of dVs above -10 m/s were not taken into account.

Assessing sinkhole susceptibility: a combined interpretation

All the dykes were built on the alluvial floodplain of the Loire River. The thickness of the alluvial deposits is therefore an important parameter to take into account, as it can have a significant impact on the stability of these structures. Indeed, the thicker the alluvial cover, the more it will be able to fill any karst voids and thus attenuate the effect at the base of the dyke. Conversely, if the cover is thin, the filling of karst voids will only be partially provided by the alluvial cover, and the body of the dykes will be exposed to disorders.

A C C E P T E D M A N U S C R I P T

The complementary use of MASW and ERT investigations allows highlighting multiple correlated Vs and resistivity anomalies on the segments, which may correspond to different stages of sinkhole propagation from the deep karst network to the surface through the alluvial cover (Figs. 12 to 14).

In the western part of the segment shown in figure 12a, a low Vs anomaly, most likely due to the presence of a karst network in depth, is identified. In addition, neither electrical flow anomalies nor historical karst-related surface disorders are observed in line with this Vs anomaly. Hence, this case highlights the development of a karstification process at the top of the marly limestone formation, without any impact on the overlying materials, and can be represented via the conceptual diagram of figure 12b.

Between abscissa -550 and -150 m of the segment shown in figure 13a, a vast low Vs anomaly, most likely due to the presence of a karst network in depth, is identified. This anomaly has good consistency with electrical flow anomalies observed in the alluvial deposits at the same abscissa. On the other hand, no karst-related surface disorders were identified in this area of the BD. This case could correspond to a flow of alluvial materials into an existing deep karst system (Fig. 13b). This type of scenario highlights a potential future collapse area.

At the centre of the segment shown in Figure 14a, 3 low Vs anomalies, most likely due to the presence of karst networks in depth, are identified. These anomalies have good consistency with electrical flow anomalies observed in the alluvial deposits, and many historical karst-related surface disorders with karst origins (BD) are in line with these Vs anomalies. This type of scenario highlights the final stage of a karst sinkhole (Fig. 14b).

Although electric investigations performed with the OhmMapper system do not reach the marly limestone bedrock where karstification occurs, some resistivity anomalies affecting the underbody of the dyke are likely to mark not only the possible presence of karst but also

A C C E P T E D M A N U S C R I P T

their consequences on near-surface formations. Potential collapse of these formations is demonstrated by patterns of sediment flow into the karst system. Thus, flowing sediment electrical anomalies affecting the underbody of the dike, linked to a strong karst presence probability extracted from Vs values may be reported as future collapse areas and constitute a good indicator for the assessment of sinkhole susceptibility. However, the electrical resistivity is known to be quite sensitive to the nature of the geological formation (clay content, water content, sand etc.). As mentionned in the section 4 through the calibration of the electrical cross-sections, resistivities are relevant for different types of alluvial deposits (from conductive clayey sediments to resistive sand and gravels). Resistivity values linked to the flowing sediment electrical anomalies can thus not be uniquely related to a weathering index or to the intensity of a sinkhole phenomenon but rather to the nature of the sediment which flows to the depth.

A statistical index extracted from the resistivity data of the electrical anomalies would doubtless more highlight geological changes in the flowing sediments than changes in the intensity of the sinkhole phenomenon. A sinkhole occurrence susceptibility level was then derived considering only a conditional index 0 or 1 according to the presence or not of an electrical anomaly (Fig. 15). An index of 1 was then added to the karst presence susceptibility level where flowing sediment electrical anomalies are identified and 0 was added elsewhere.

Assessing the sinkhole susceptibility level over the entire set of investigated segments then became possible (Fig. 12, 13, 14 and 16a).

Discussion and Conclusion

To assess the reliability of the methodology, the obtained sinkhole susceptibility level was interpolated using a kriging technique over the entire set of investigated segments to obtain a

A C C E P T E D M A N U S C R I P T

sinkhole susceptibility map. The map shows sinkhole susceptibility on an empirical scale based on geophysical indicators. It is first compared to the projected karst-related surface disorders on the dykes from the BD cavités (Fig. 16a). Karstic caves in plan view resemble ink blots or Rorschach patterns. Ramiform irregular cavities and galleries wander threedimensionally with branches extending outward, up to possibly hundred meters, from the main areas of development [START_REF] Palmer | Origin and morphology of limestone caves[END_REF]. Thereby, assuming that a disorder related to the same karstic network will not necessarily appear directly on the measured line, a search area of 100 m was used for the interpolation and as a buffer area on either side of dykes for the selection of database events. There is an acceptable consistency between areas of high sinkhole susceptibility and the projected locations of karst-related surface disorders, 60 % of them being located in high susceptibility areas (  2). This provides a consistency check of the sensibility of the methodology to karst phenomena, although the susceptibility level should be more related to the prediction of future disorders than to the detection of former ones.

Secondly, the sinkhole susceptibility map is compared to the projected sinkholes extracted from the French National database of land movements (BD mvt) which occurred after 2008, i. e. after the acquisition of geophysical data (Fig 16a). Here again, a buffer area of 100 m on either side of dykes was used for the selection of database events. One can see that there is a better correlation between these two elements, 86% of the sinkholes which occured after the geophysical survey being located within high susceptibility areas (  2) (Fig. 16c). Hence, it appears too that some areas which were not disturbed at the time of the geophysical acquisition, especially at the north of St-Pryvé-St-Mesmin city and at the north-east of Sandillon city, present sinkholes few years later. This observation highlights the effective ability of the methodology to predict areas of potential future collapse and provides a reliability check of the presented sinkhole susceptibility map. However, this statistical analysis should be carefully interpreted according to the few number of sinkholes which

A C C E P T E D M A N U S C R I P T
occurred after the geophysical survey. Let's note that the low and very low susceptibility categories do not indicate that sinkholes have not or cannot occur in these areas as shown by the sinkhole that occured after the geophysical survey at the east of the Sandillon city in an area of low sinkhole susceptibility level. This could be explained by the variability of the collapse phenomenon kinetics, the sudden occurrence of a sinkhole due to a surficial layer breaking being undetectable by the methodology based on the presence of precursor resistivity anomalies near the surface.

As a consequence, a limitation of the susceptibility map is the possibility that other causative factors not considered in our analysis may result in sinkhole occurence. Local hydrology and precipitations are two of these factors. This map is not intended for use at scales other than the published scale and such use may result in erroneous conclusions regarding sinkhole hazard. The map is designed for regional planning use to determine areas where local sinkhole hazards may exist on the dykes and where more detailed planning maps and geotechnical studies are needed, particularly in areas of high and moderate susceptibility.

ERT and MASW coupling is a promising approach to detecting and characterising karst structures in urban environments, with the two methods being complementary in depth of investigation and origin of detected anomalies. MASW surveying was effective at mapping the interface between overburden alluvial deposits and the marly limestone bedrock. MASW surveys successfully highlighted zones of low velocity and therefore mechanical weaknesses corresponding to karstified areas in limestones. Such detection allowed the assessment of a karst presence probability level with a better resolution (10 m) than former studies in the same area. For comparison, [START_REF] Donsimoni | Détermination de la présence de karst sous les levées domaniales du bassin de la Loire et réduction du risque d'effondrement de la digue lié à la présence de ces conduits souterrains naturels. Tranche 1: analyse géologique et hydrogéologique[END_REF] estimated karst presence probability over the entire Loire River and especially in the study area with a resolution of 500 m using 3 indirect geological and hydrogeological indicators: the sensibility of the substratum to the dissolution phenomenon according to its nature, structure and fracturation degree criteria, the

A C C E P T E D M A N U S C R I P T
thickness of the alluvial sedimentary cover and the water circulation mode between the alluvial and karst aquifer. In addition to the resolution, the exposed methodology also has the advantage of providing a probability level based on direct observations of the subsurface state, making the evaluation more reliable.

Resistivity surveying at the top of the dykes was also useful to detect potential flow of nearsurface materials into karst at depth. By integrating results obtained using these two geophysical techniques, it was possible to highlight areas of interaction between karst and the near-surface formations that were previously poorly understood. Providing this characterization of near surface materials using the resistivity survey, the integrated results allowed evaluating a sinkhole susceptibility empirical level at the field scale which confers to the methodology its predictive purpose in opposition to the usual characterization aim of geophysical applications. This non-intrusive, rapid and inexpensive approach seems promising for mapping long segments of river dykes and aid engineers in predicting any potential dangers or hidden hazards. The associated sinkhole susceptibility level obtained by incrementing the karst presence probability level using the electrical resistivity anomaly conditional index.

ACCEPTED MANUSCRIPT 44

Highlights -A series of geophysical experiments were conducted to detect karst features below a dyke system -The mapping of weak areas in limestone allowed the assessment of a karst presence probability level -Resistivity surveying was useful to detect potential flow of near-surface materials into karsts -By integrating results, interaction between karsts and the near-surface formations was highlighted -It then allowed mapping of these sinkhole-susceptible areas with a better accuracy than former studies

  mentioned in the previous paragraph, potentially karstified areas are highlighted by a sudden drop in Vs, revealing mechanical weakness associated with karst weathering. The methodology for assessing the probability of karst presence is then based on a step-by-step analysis of the lateral variation of the Vs parameter and is illustrated in figure 10. The first A C C E P T E D M A N U S C R I P T step of the methodology consists of extracting the Vs parameter at a given depth in marly limestone. To do this for each segment, a depth of -10 m below the calibrated or assumed top of marly limestone was chosen to ensure that the extraction is reliable and representative of the marly limestone lithology Vs values. Next, an absolute difference (dVs) relative to the average Vs value of each segment at this depth was computed. This step eliminates largescale Vs variations associated with changes in the lithology and allows retention of only local changes in the parameter. Smooth dVs curves were then obtained by applying a 10 m moving average filter to eliminate variations induced by 2D interpolation of Vs profiles. Finally, an extraction of negative dVs values was made. In the end, this step allows the identification of areas that highlight lower mechanical parameters than the mean of the segment. The statistical evaluation of dVs values over the entire studied area allows merging these values in a histogram (Fig. 11) representing their occurrence frequency according to the dVs classes. Such an analysis permits identifying 3 groups of different dVs, each most likely linked to different weathering steps of limestone during the karstification process. Based on these three groupings, karst presence probability was assigned a level of 1 for dVs ranging from -10 to -45 m/s, 3 for dVs ranging from -45 to -95 m/s and 5 for dVs over -95 m/s.

Figure 1 :

 1 Figure captions

Figure 2 :

 2 Figure 2: a) Location map of the study area showing geophysical profiles in black solid lines.

Figure 3 :

 3 Figure 3: Photography of a sinkhole at the top of a Loire River dyke.

Figure 4 :

 4 Figure 4: ERT typical a) raw and b) filtered apparent resistivity pseudo-section obtained on a

Figure 5 :

 5 Figure 5: Typical examples of electrical resistivity calibration for the clayey silt (a), clayey

Figure 6 :

 6 Figure 6: Typical examples of interpreted electrical resistivity sections, highlighting

Figure 7 :

 7 Figure 7: Typical examples of interpreted electrical resistivity sections, highlighting sudden

Figure 8 :

 8 Figure 8: Typical examples of Vs calibration. The lithological description of the destructive

Figure 9 :

 9 Figure 9: Typical examples of interpreted MASW sections. Vs anomalies highlight potential

Figure 10 :

 10 Figure 10: successive steps of the karst presence probability assessment methodology. a)

Figure 11 :

 11 Figure 11: Histogram of dVs values in the entire studied area. The statistical Gaussian

Figure 12 :

 12 Figure 12: a) An example of correlations between the Vs and electrical resistivity section on

Figure 13 :

 13 Figure 13: a) An example of correlations between the Vs and electrical resistivity section on

Figure 14 :

 14 Figure 14: a) An example of correlations between the Vs and electrical resistivity section on

Figure 15 :

 15 Figure 15: a) An example of the karst presence probability level obtained along the segment

Figure 16

 16 Figure 16: a) Sinkhole susceptibility map interpolated over the entire set of investigated

Figure 4 A

 4 Figure 4

Figure 6 A

 6 Figure 6

Figure 9 A

 9 Figure 9

Acknowledgments

We thank the LRPC Blois (Laboratoire Régional des Ponts et Chaussées) for having made available their geotechnical surveys of the Loire River dykes. We also thank the DREAL Centre (Direction Régionale de l'Environnement, de l'Aménagement et du Logement) for providing useful information to the geophysicists. We are very grateful to the BRGM technicians for their field work. This work was supported by BRGM-DREAL Centre public research program 2006-CEN-18. We thank the anonymous reviewers for their constructive comments which helped to improve the manuscript.

A C C E P T E D

M A N U S C R I P T


Felipe Yanez

Francis Bach

Primal-Dual Algorithms for Non-negative Matrix Factorization with the Kullback-Leibler Divergence

Non-negative matrix factorization (NMF) approximates a given matrix as a product of two non-negative matrices. Multiplicative algorithms deliver reliable results, but they show slow convergence for high-dimensional data and may be stuck away from local minima. Gradient descent methods have better behavior, but only apply to smooth losses such as the leastsquares loss. In this article, we propose a first-order primal-dual algorithm for non-negative decomposition problems (where one factor is fixed) with the KL divergence, based on the Chambolle-Pock algorithm. All required computations may be obtained in closed form and we provide an efficient heuristic way to select step-sizes. By using alternating optimization, our algorithm readily extends to NMF and, on synthetic examples, face recognition or music source separation datasets, it is either faster than existing algorithms, or leads to improved local optima, or both.

Introduction

The current development of techniques for big data applications has been extremely useful in many fields including data analysis, bioinformatics and scientific computing. These techniques need to handle large amounts of data and often rely on dimensionality reduction; this is often cast as approximating a matrix with a low-rank element.

Non-negative matrix factorization (NMF) is a method that aims at finding part-based, linear representations of non-negative data by factorizing it as the product of two low-rank nonnegative matrices [START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error[END_REF][START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF]. In 2000, two multiplicative algorithms for NMF were introduced by Lee and Seung, one that minimizes the conventional least-squares error, and other one that minimizes the generalized Kullback-Leibler (KL) divergence [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]. These algorithms extend to other losses and have been reported in different applications, e.g., face recognition [START_REF] Wang | Non-negative matrix factorization framework for face recognition[END_REF], music analysis (Févotte et al., 2009), and text mining [START_REF] Guduru | Text mining with support vector machines and non-negative matrix factorization algorithms[END_REF]. An important weakness of multiplicative algorithms is their slow convergence rate in high-dimensional data and their susceptibility to become trapped in poor local optima [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF]. Gradient descent methods for NMF provide additional flexibility and fast convergence [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF][START_REF] Kim | Fast projection-based methods for the least squares nonnegative matrix approximation problem[END_REF][START_REF] Gillis | Nonnegative Matrix Factorization: Complexity, Algorithms and Applications[END_REF]. These methods have been extensively studied for the minimization of the least-squares error [START_REF] Lin | Projected gradient methods for nonnegative matrix factorization[END_REF][START_REF] Kim | Fast projection-based methods for the least squares nonnegative matrix approximation problem[END_REF].

The goal of this paper is to provide similar first-order methods for the KL divergence, with updates as cheap as multiplicative updates. Our method builds on the recent work of [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF] which consider the alternating direction method of multipliers (ADMM) adapted to this problem. We instead rely on the Chambolle-Pock algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], which may be seen as a linearized version of ADMM, and thus we may reuse some of the tools developed by [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF] while having an empirically faster algorithm.

Contributions

The main contributions of this article are as follows:

-We propose a new primal-dual formulation for the convex KL decomposition problem in Section 3.1, and an extension to the non-convex problem of NMF by alternating minimization in Section 3.5.

-We provide a purely data-driven way to select all step-sizes of our algorithm in Section 3.3.

-In our simulations in Section 4 on synthetic examples, face recognition or music source separation datasets, our algorithm is either faster than existing algorithms, or leads to improved local optima, or both.

-We derive a cheap and efficient implementation (Algorithm 2). Matlab code is available online at: anonymized website 2 Problem Formulation Two widely used cost functions for NMF are the conventional least-squares error (not detailed herein), and the generalized KL divergence

Let V ∈ R n×m
D(V WH) = - m i=1 n j=1 V ij log (WH) ij V ij + 1 + m i=1 n j=1 (WH) ij . (1)
In this work, only the KL divergence is considered. Therefore, the optimization problem is as follows: minimize

W,H ≥ 0 D(V WH).
(2)

We recall that the previous problem is non-convex in both factors simultaneously, whereas convex in each factor separately, i.e., the non-negative decomposition (ND) problems, minimize

W ≥ 0 D(V WH) (3)
and minimize

H ≥ 0 D(V WH), (4)
are convex.

We now present two algorithms for NMF, multiplicative updates [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF], and the ADMM-based approach [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF]. [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] introduced two multiplicative updates algorithms for NMF. One minimizes the conventional least-squares error, and the other one minimizes the KL divergence.

Multiplicative updates

The NMF problem, for both losses, is a non-convex problem in W and H simultaneously, but convex with respect to each variable taken separately; this make alternating optimization techniques, i.e., solving at each iteration two separate convex problems, very adapted: first fixing H to estimate W, and then fixing W to estimate H [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF]Févotte et al., 2009). The multiplicative updates algorithms (like ours) follow this approach.

For the KL divergence loss, the multiplicative update rule [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] for W and H is as follows and may be derived from expectation-maximization (EM) for a certain probabilistic model [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF][START_REF] Févotte | Nonnegative matrix factorizations as probabilistic inference in composite models[END_REF]:

W ia ← W ia n µ=1 H aµ V iµ /(WH) iµ n ν=1 H aν , and
H aµ ← H aµ m i=1 W ia V iµ /(WH) iµ m k=1 W ka .
The complexity per iteration is O(rmn).

2.2 Alternating direction method of multipliers (ADMM) [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF] propose an ADMM technique to solve Problem (2) by reformulating it as minimize D(V X)

subject to X = YZ Y = W, Z = H W ≥ 0, H ≥ 0.
The updates for the primal variables W, H, X, Y and Z are as follows and involve certain proximal operators for the KL loss which are the same as ours in Section 3.2:

Y ⊤ ← ZZ ⊤ + I -1 ZX ⊤ + W ⊤ + 1 ρ Zα ⊤ X -α ⊤ Y Z ← Y ⊤ Y + I -1 Y ⊤ X + H + 1 ρ Y ⊤ α X -α Z X ← (ρYZ -α X -1) + (ρYZ -α X -1) 2 + 4ρV 2ρ W ← Y + 1 ρ α Y + H ← Z + 1 ρ α Z + .
Note that the primal updates require solving linear systems of sizes r × r, but that the overal complexity remains O(rmn) per iteration (the same as multiplicative updates).

The updates for the dual variables α X , α Y and α Z are then:

α X ← α X + ρ (X -YZ) α Y ← α Y + ρ (Y -W) α Z ← α Z + ρ (Z -H) .
This formulation introduces a regularization parameter, ρ ∈ R + , that needs to be tuned (in our experiments we choose the best performing one from several candidates).

Our approach has the following differences: (1) we aim at solving alternatively convex problems with a few steps of primal-dual algorithms for convex problems, as opposed to aiming at solving directly the non-convex problem with an iterative approach, (2) for the convex decomposition problem, we have certificates of guarantees, which can be of used for online methods for which decomposition problems are repeatedly solved [START_REF] Lefèvre | Online algorithms for nonnegative matrix factorization with the Itakura-Saito divergence[END_REF] and (3) we use a different splitting method, namely the one of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], which does not require matrix inversions, and which allows us to compute all step-sizes in a data-driven way.

Proposed Method

In this section we present a formulation of the convex KL decomposition problem as a first-order primal-dual algorithm (FPA), followed by the proposed NMF technique.

Primal and dual computation

We consider a vector a ∈ R p + and a matrix K ∈ R p×q + as known parameters, and x ∈ R q + as an unknown vector to be estimated, where the following expression holds, a ≈ Kx, and we aim at minimizing the KL divergence between a and Kx. This is equivalent to a ND problem as defined in Problems (3) and (4), considering a as a column of the given data, K as the fixed factor, and x as a column of the estimated factor, i.e., in Problem (3) a and x are column vectors of V ⊤ and W ⊤ with the same index and K is H ⊤ , and in Problem (4) a and x are columns of V and H with the same index and K is W.

The convex ND problem with KL divergence is thus minimize

x∈R q + - p i=1 a i (log(K i x/a i ) + 1) + p i=1 K i x, (5)
which may be written as

minimize x∈X F (Kx) + G(x), (6)
with

F (z) = - p i=1 a i (log(z i /a i ) + 1) G(x) = ½ x 0 + p i=1 K i x.
Following [START_REF] Pock | An algorithm for minimizing the mumfordshah functional[END_REF]; [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], we obtain the dual problem

maximize y∈Y -F * (y) -G * (-K * y), with F * (y) = sup z y ⊤ z -F (z) = - p i=1 a i log (-y i ) G * (y) = sup x y ⊤ x -G(x) = ½ y K ⊤ 1 .
We then get the dual problem maximize

K ⊤ (-y) K ⊤ 1 a ⊤ log (-y) . (7)
In order to provide a certificate of optimality, we have to make sure that the constraint K ⊤ (-y) K ⊤ 1 is satisfied. Therefore, when it is not satisfied, we project as follows:

y ← y/ max{K ⊤ (-y) ⊘ K ⊤ 1},
where ⊘ represents the entry-wise division operator.

Primal-dual algorithm

The general FPA framework of the approach proposed by Chambolle and Pock for Problem ( 6) is presented in Algorithm 1.

Algorithm 1: First-order primal-dual algorithm.

Select K ∈ R p×q + , x ∈ R q + , σ > 0, and τ > 0; Set x = x old = x, and y = Kx;

for N iterations do y ← prox σF * (y -σK x); x ← prox τ G (x -τ K * y); x ← 2x -x old ; x old ← x; end return x ⋆ = x.
Algorithm 1 requires the computation of proximal operators prox σF * (y) and prox τ G (x). These are defined as follows:

prox σF * (y) = arg min v v -y 2 2σ + F * (v) , and
prox τ G (x) = arg min u u -x 2 2τ + G(u) .
For further details, see [START_REF] Boyd | Convex Optimization[END_REF][START_REF] Rockafellar | Convex Analysis[END_REF].

Using the convex functions F * and G, we can easily solve the problems for the proximal operators and derive the following closed-form solution operators

prox σF * (y) = 1 2 y -y • y + 4σa , and
prox τ G (x) = x -τ K ⊤ 1 + .
The detailed derivation of these operators may be found in the Appendix, the first one was already computed by [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF].

Automatic heuristic selection of σ and τ

In this section, we provide a heuristic way to select σ and τ without user intervention, based on the convergence result of Chambolle and Pock (2011, Theorem 1), which states that (a) the step-sizes have to satisfy τ σ K 2 ≤ 1, where K = max{ Kx : x ∈ X with x ≤ 1} is the largest singular value of K; and (b) the convergence rate is controlled by the quantity

C = y 0 -y ⋆ 2 2σ + x 0 -x ⋆ 2 2τ ,
where (x ⋆ , y ⋆ ) is an optimal primal/dual pair. If (x ⋆ , y ⋆ ) was known, we could thus consider the following minimization problem with the constraint τ σ K 2 ≤ 1:

min σ,τ y 0 -y ⋆ 2 2σ + x 0 -x ⋆ 2 2τ ⇐⇒ min σ y 0 -y ⋆ 2 2σ + x 0 -x ⋆ 2 2 σ K 2 .
Applying first order conditions, we find that

σ = y 0 -y ⋆ x 0 -x ⋆ 1 K and τ = x 0 -x ⋆ y 0 -y ⋆ 1 K .
However, we do not know the optimal pair (x ⋆ , y ⋆ ) and we use heuristic replacements. That is, we consider the unknown constants α and β, and assume that x ⋆ = α1 and y ⋆ = β1 solve Problems ( 5) and ( 7). Letting (x 0 , y 0 ) = (0, 0) we have

x 0 -x ⋆ = |α| √ q and y 0 -y ⋆ = |β| √ p.
Plugging x ⋆ to Problem (5), we are able to find that α = 1 ⊤ a 1 ⊤ K1 > 0. Now, using optimality conditions:

y ⋆ • (Kx ⋆ ) = -a, we obtain β = -1.
The updated version of the parameters is:

σ = p q 1 α K and τ = q p α K .
Finally, an automatic heuristic selection of step sizes σ and τ is as follows:

σ = √ p p i=1 K i 1 √ q K p i=1 a i and τ = √ q p i=1 a i √ p K p i=1 K i 1 .
Note the invariance by rescaling of a and K.

Implementation

The proposed method is based on Algorithm 1. The required parameters to solve each ND problem are

• Problem (3): a ← V ⊤ i K ← H ⊤ x ← W ⊤ i σ ← m r 1 ⊤ H1 1 ⊤ (V ⊤ ) i H τ ← r m 1 ⊤ (V ⊤ ) i 1 ⊤ H1 H • Problem (4): a ← V i K ← W x ← H i σ ← n r 1 ⊤ W1 1 ⊤ V i W τ ← r n 1 ⊤ V i 1 ⊤ W1 W .
The previous summary treats each ND problem by columns. For algorithmic efficiency, we work directly with the matrices, e.g., a ∈ R n×m + instead of R n + . We also include normalization steps such that the columns of W have sums equal to 1. The stopping criteria is enabled for maximum number of iterations (access to data) and for duality gap tolerance.

Extension to NMF

A pseudo-code of the first-order primal-dual algorithm for non-negative matrix factorization can be found in Algorithm 2. It corresponds to alternating between minimizing with respect to H and minimizing with respect to W. A key algorithmic choice is the number of inner iterations iter N D of the convex method, which we consider in Section 4.

The running-time complexity is still O(rnm) for each inner iterations. Note moreover, that computing the largest singular value of H or W (required for the heuristic selection of step-sizes everytime we switch from one convex problem to the other) is of order O(r max{m, n}) and is thus negligible compared to the iteration cost.

Extension to topic models

Probabilistic latent semantic analysis [START_REF] Hofmann | Probabilistic latent semantic analysis[END_REF] or latent Dirichlet allocation [START_REF] Blei | Latent dirichlet allocation[END_REF], generative probabilistic models for collections of discrete data, have been extensively used in text analysis. Their formulations are related to ours in Problem ( 5), where we just need to include an additional constraint: 1 ⊤ x = 1. In this sense, if we modify G, i.e., G(x) = ½{1 ⊤ x = 1; x 0} + 1 ⊤ Kx, we can use Algorithm 1 to find the latent topics. It is important to mention that herein prox τ G (x) does not have a closed solution, but can be efficiently solved with dedicated methods for orthogonal projections on the simplex [START_REF] Maculan | A linear-time median-finding algorithm for projecting a vector on the simplex of R n[END_REF]).

Experimental Results

The proposed technique was tested on synthetic data, the CBCL face images database and a music excerpt from a recognized jazz song by Louis Armstrong & His Hot Five. The performance of the proposed first-order primal-dual algorithm (FPA) was compared against the traditional multiplicative updates algorithm (MUA) by [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] and the contemporary alternating direction method of multipliers (ADMM) by [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF]. The three techniques were implemented in Matlab.

Synthetic data

A given matrix V of size n = 200 and m = 1000 is randomly generated from the uniform distribution U (0, 750). The low-rank element was set to r = 15. Initializations W 0 and H 0 are defined by the standard normal distribution's magnitude plus a small offset.

ND problem

To examine the accuracy of our method, we first apply Algorithm 2 to convex ND problems for fixed values of n, m and r, solving separately Problems (3) and (4). For both problems, we set Algorithm 2: Proposed technique.

Select V ∈ R n×m + , W 0 ∈ R n×r + ,
1 ⊤ H1 1 ⊤ V ⊤ H 1, τ = r m 1 ⊤ V ⊤ 1 ⊤ H1 H 1, and H(-χ ⊤ ) ≤ H1; for iter N D iterations do χ ⊤ ← χ ⊤ -σ • WH ⊤ ; χ ⊤ ← 1 2 χ ⊤ -χ ⊤ • χ ⊤ + 4σ • V ⊤ ; W ⊤ ← W ⊤ -τ • H χ ⊤ + 1 + ; W⊤ ← 2W ⊤ -W ⊤ old ; W ⊤ old ← W ⊤ ; end Normalize H and set σ = n r 1 ⊤ W1 1 ⊤ V W 1, τ = r n 1 ⊤ V 1 ⊤ W1 W 1, and W ⊤ (-χ) ≤ W ⊤ 1; for iter N D iterations do χ ← χ -σ • W H ; χ ← 1 2 χ - √ χ • χ + 4σ • V ; H ← H -τ • W ⊤ (χ + 1) + ; H ← 2H -H old ; H old ← H; end end return W ⋆ = W, and H ⋆ = H.
the number of iterations of the traditional MUA and contemporary ADMM to 1200, as well as for the proposed FPA. Optimal factors W ⋆ and H ⋆ are obtained by running 5000 iterations of the MUA, starting from W 0 and H 0 . For the first ND problem, we fix H to H ⋆ and estimate W starting from W 0 ; for the second one, we fix W to W ⋆ and estimate H from H 0 . The optimal regularization parameter of ADMM, the tuning parameter that controls the convergence rate, is ρ = 0.15 (small values imply larger step sizes, which may result in faster convergence but also instability). Figure 1 (a-b) present us the distance to optimum of MUA and ADMM, as well as for the primal and dual of our technique that reveals strong duality. The FPA and ADMM algorithm converge to the same point, whereas the MUA exhibits slow convergence. Note moreover the significant advantage towards our algorithm FPA, together with the fact that we set automatically all step-sizes. In Figure 1 (c-d), we illustrate the distance to optimum versus time of the three methods.

NMF problem

The setting is slightly different as in the ND experiment, we increased the problem dimension to n = 250, m = 2000 and r = 50, and repeat both previously presented experiments. For all methods, we set the number of iterations to 3000. The parameter iter N D indicates the number of iterations to solve each ND problem. We set iter N D to 5 iterations. To have a fair comparison between algorithms, for FPA, the number of iterations means access to data, i.e., we use 5 iterations to solve Problem (3) (as well as for Problem (4)), and repeat this 600 times. The optimal regularization parameter of the ADMM is here ρ = 1.

In Figure 2 we present the objective function of the three algorithm for the non-convex Problem (2). The MUA initially reports high decrement in the objective, but as time increases it exhibits evident slow tail convergence. The FPA primal objective decreases dramatically in only seconds (few iterations), and furthermore, it presents fast tail convergence achieving the lowest objective value. The ADMM has poor initial performance, but then achieves an optimal value similar to the one obtained by FPA. In order to show that FPA converges faster and with lower computational cost, we store the cost function values and computation times at each iteration. The total time required by the FPA was 190s, whereas 205s by the ADMM and 473s by the MUA. Then we analyze the ADMM and MUA for the same time 190s (the vertical dotted line in Figure 2 (b)), i.e., 2786 and 1211 iterations, respectively: the competitive algorithms have a significantly larger cost function for the same running time. The result of this experiment is illustrated in Figure 2 (b). The results considering the objective function versus iteration number may be found in the Appendix.

NMF with warm restarts

The problem dimension is n = 150, m = 2000 and r 1 = 50. We run 3000 iterations of each method using initializations W 0 and H 0 ; then we increase ten times the low-rank element, r 2 = 100; and finally run 2000 more iterations, producing W 2 and H 2 . The idea is to use as initializations the estimations obtained after the first 3000 iterations, W 1 and H 1 , considering that the low-rank element changed. A trivial solution could be to include random entries so that W 1 and H 1 have the proper dimensions, but that increases the objective value, diminishing the estimations. On the other hand, if we include zero entries so that W 1 and H 1 have the proper dimensions, we would be in a saddle-point where none of the algorithms could perform. However, if we set only one factor with zero entries, [W 1 , c1] ∈ R n×r 2 with c = 0, and the other one with non-zero values, [H 1 ; ν] ∈ R r 2 ×m , we still maintain the same last objective value and perform FPA. In this situation, MUA cannot perform either (because of the absorbing of zeros), therefore we try some values of c to run the algorithm. Figure 3 illustrates the proposed experiment. Notice that as c → 0, the MUA starts to get stuck in a poor local optima, i.e., the one obtained with W 1 and H 1 . ADMM has a similar behavior as FPA, therefore, it is not displayed.

MIT-CBCL Face Database #1

We use the CBCL face images database [START_REF] Sung | Learning and Example Selection for Object and Pattern Recognition[END_REF] composed of m = 2429 images of size n = 361 pixels. The low-rank element was set to r = 49. Figure 4 shows samples from the database.

Our next experiment is to determine the optimal the number of iterations for the current database. Therefore, we run 3000 iterations of FPA, using 3, 5, 10 and 15 iterations for the ND problem. The MUA and ADMM (ρ=50) algorithms are performing as well. Figure 5 illustrates the decay of the objective function of the FPA, MUA and ADMM algorithms. We appreciate that setting the number of iterations to 3 yield to over-alternation, whereas using 15 or even more iterations result in an under-alternating method. Using 10 iterations reveal good performance, but the best trade-off is obtained with 5 iterations. Therefore, we set iter N D = 5, i.e., the number of iterations to solve Problem (3) and Problem (4). All following results in the MIT-CBCL Face Database #1 are with the same setting.

Finally, in Figure 6 (a) we present the objective function of the three algorithm for the nonconvex Problem (2), where all algorithms perform similarly. However, in the zoomed Figure 6 (b) we can appreciate that the MUA presents the slowest convergence, whereas the proposed method the fastest one. The results considering the objective function versus iteration number may be found in the Appendix.

Music excerpt from the song "My Heart (Will Always Lead Me Back to You)"

The last experiment is to decompose a 108-second-long music excerpt from "My Heart (Will Always Lead Me Back to You)" by Louis Armstrong & His Hot Five in the 1920s (Févotte et al., 2009). The time-domain recorded signal is illustrated in Figure 7. The recording consists of a trumpet, a double bass, a clarinet, a trombone, and a piano. The recorded signal is original unprocessed mono material contaminated with noise. The signal was downsampled to 11025 kHz, yielding 1,19×10 6 samples. The Fourier Transform of the recorded signal was computed using a sinebell analysis window of length 23 ms with 50% overlap between two frames, leading to m = 9312 frames and n = 129 frequency bins. Additionally, we set r = 10. Figure 8 illustrates the previously described spectrogram. The decomposition of the song is produced by the three algorithms. We initialize them with the same random values W 0 and H 0 . For a fair competition, the number of iterations is set to 5000 for MUA and ADMM, and for our algorithm FPA we consider it as access to data, i.e., we use 5 iterations for the ND, repeating it 1000 times. For comparison, we measure the computation time of the three techniques. FPA has a run time of 13 min, whereas the ADMM (ρ = 10) one of 15 min and the MUA one of 80 min. In this experiment, Figure 9 illustrates the evolution of the objective of the three techniques along time. Initially the MUA obtained the lowest objective value, but as previously discussed, as the number of iterations increases the MUA starts exhibiting evident slow tail convergence and since approximately 100s it is reached by the FPA and shows no further substantial decrement, i.e., it gets stuck in a worse local optima. FPA converges to a slight lower cost value, overpassing MUA. Finally, ADMM reveals a slow initial performance on this dataset, but later converges to a similar point as the previous algorithms. The results considering the objective function versus iteration number may be found in the Appendix.

Conclusion

We have presented an alternating projected gradient descent technique for NMF that minimizes the KL divergence loss; this approach solves convex ND problems with the FPA. Our approach demonstrated faster convergence than the traditional MUA by [START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF] and contemporary ADMM by [START_REF] Sun | Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence[END_REF]. The FPA introduces a new parameter, the number of iterations for each convex ND problem. Experiments reveal that the number of iterations is mostly bounded between 3 and 10 iterations, which leads to a trivial tuning by the user. Therefore, our algorithm affords reasonable simplicity, where further user manipulation is not required. Finally, an extension to latent Dirichlet allocation and probabilistic latent semantic indexing can be easily implemented using our proposed method, thus allowing to go beyond the potential slowness of the expectation-maximization (EM) algorithm.

prox τ G (x) = arg min u 0 u -x 2 2τ + G(u) = arg min u 0 u -x 2 2τ + n i=1 K i u = x -τ K ⊤ 1 + .
Synthetic data: additional results

NMF problem

A given matrix V of size n = 250 and m = 2000 is randomly generated from the uniform distribution U (0, 750). The low-rank element was set to r = 50. For the three methods, we set the number of iterations to 3000. We set iter N D to 5 iterations. The optimal tuning parameter of the ADMM is here ρ = 1. In Figure 10 we present the objective function versus iteration number of the three algorithms for the non-convex NMF problem. MIT-CBCL Face Database #1: additional results

ND problem

We solve convex ND problems for fixed values of n, m and r, setting the number of iterations of all algorithms to 1500. Optimal factors W ⋆ and H ⋆ are obtained by running 5000 iterations of the MUA. The optimal tuning parameter of the ADMM is here ρ = 0.1. Figure 11 (a-b) presents us the distance to optimum of the MUA and ADMM, as well as for the primal and dual of our technique that reveals strong duality in all experiments. In Figure 11 (c-d), we illustrate the distance to optimum versus time of the three methods.

NMF problem

For all methods, we set the number of iterations to 3000. We set iter N D to 5 iterations. The optimal tuning parameter of the ADMM is here ρ = 50. In Figure 12 we present the objective function versus iteration number of the three algorithms.

Features learned from the CBCL face image database

The features learned from the CBCL face image database obtained with the three algorithms is presented in Figure 13. The figure reveals the parts-based learned by the algorithm, i.e. W. "My Heart (Will Always Lead Me Back to You)": additional results

ND problem

We solve convex ND problems for fixed values of n, m and r, setting the number of iterations of all algorithms to 2000. Optimal factors W ⋆ and H ⋆ are obtained by running 5000 iterations of the MUA. The optimal tuning parameter of the ADMM is here ρ = 0.5. Figure 14 (a-b) presents us the distance to optimum of the MUA and ADMM, as well as for the primal and dual of our technique that reveals strong duality in all experiments. In Figure 14 (c-d), we illustrate the distance to optimum versus time of the three methods.

NMF problem

For all methods, we set the number of iterations to 5000. We set iter N D to 5 iterations. The optimal tuning parameter of the ADMM is here ρ = 10. In Figure 15 we present the objective

Features learned from the song

The decomposition of the song by Louis Armstrong and band obtained with the proposed FPA is presented in Figure 16, revealing the parts-based learned by the algorithm, i.e., W. The timedomain signal is recovered from Wiener filtering (Févotte et al., 2009).

+

  denote the n × m given matrix formed by m non-negative column vectors of dimensionality n. Considering r ≤ min(n, m), let W ∈ R m×r + and H ∈ R r×n + be the matrix factors such that V ≈ WH.

  and H 0 ∈ R r×m + ; Set W = W = W old = W 0 , H = H = H old = H 0 , and χ = WH; while stopping criteria not reached do Normalize W and set σ = m r

  Estimate H given W ⋆ .

Figure 1 :

 1 Figure 1: ND on synthetic data. (a-b) Distance to optimum versus iteration number. Distance to optimum reveals the difference between the values of the objective function and optimal point, p ⋆ . In the case of the dual function values, the distance to optimum is the difference between p ⋆ and the dual points. (c-d) Distance to optimum versus time.

  Objective function versus time (zoomed).

Figure 2 :

 2 Figure 2: NMF on synthetic data. Recall that the dual function is not presented due to the non-convexity of the NMF problem.

Figure 3 :

 3 Figure 3: NMF with warm restarts on synthetic data. Value of the objective function at each iteration.

Figure 4 :

 4 Figure 4: MIT-CBCL Face Database #1 samples.

Figure 5 :

 5 Figure 5: NMF on the CBCL database. Value of the objective function at each iteration solving Problem (2) varying the number of iterations to solve each ND problem.

  Objective function versus time (zoomed).

Figure 6 :

 6 Figure 6: NMF on the CBCL face image database.

Figure 7 :

 7 Figure 7: Time-domain recorded signal.

Figure 8 :

 8 Figure 8: Log-power spectrogram.

  Objective function versus time (zoomed).

Figure 9 :

 9 Figure 9: NMF on an excerpt of Armstrong's song.

Figure 10 :

 10 Figure 10: NMF on synthetic data. It is important to recall that the dual function is not presented due to the non-convexity of the NMF problem.

  Estimate H given W ⋆ .

Figure 11 :

 11 Figure 11: Distance to optimum versus (a-b) iteration number, and (c-d) time.

Figure 12 :

 12 Figure 12: NMF on the CBCL face image database.

Figure 13 :

 13 Figure 13: Features learned from the CBCL face image database.

Figure 14 :

 14 Figure 14: Distance to optimum versus (a-b) iteration number, and (c-d) time.

Figure 15 :

 15 Figure 15: NMF on an excerpt of Armstrongs song.

Acknowledgements

This work was partially supported by a grant from the European Research Council (ERC SIERRA 239993).

Appendix Derivation of proximal operators

The definition of the proximal operator of F * and G, i.e., (I + σ∂F * ) -1 (y) and (I + τ ∂G) -1 (x), respectively, is as follows [START_REF] Pock | An algorithm for minimizing the mumfordshah functional[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]:

where ∂F * and ∂G are the subgradients of the convex functions F * and G.

To facilitate the computation of (I + σ∂F * ) -1 (y), we consider Moreau's identity

Let us consider the variable v ∈ Y, and using Moreau's identity, we can compute

Applying first order conditions to obtain the minimum:

Finally, the proximal operator is as follows:

For the second proximal operator, we consider u ∈ X and compute prox τ G (x) as