Hervé Michallet
email: herve.michallet@legi.cnrs.fr

B Gerben Ruessink
email: b.g.ruessink@uu.nl

Mariana Vieira

Lima Matias Da Rocha

Anouk De Bakker

Dominic A Van Der A

Andrea Ruju
email: andrea.ruju@unican.es

Paulo A Silva
email: psilva@ua.pt

Nadia Sénéchal
email: n.senechal@epoc.u-bordeaux1.fr

Vincent Marieu

Marion Tissier

Lima Matias

Gerben Ruessink

Mariana V L Rocha

Dominic Van Der

Rafael Almar
email: rafael.almar@ird.fr

Tiago Abreu

Florent Birrien

Laure Vignal

Eric Barthélemy

Dominique Mouazé
email: dominique.mouaze@unicaen.fr

Rodrigo Cienfuegos

Peter Wellens
email: peter.wellens@deltares.nl

GLOBEX: Wave dynamics on a shallow sloping beach

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

INTRODUCTION

Morphological change in the nearshore zone arises from the complex interplay of water motion, sediment transport, and the morphology itself. The governing water motion, predominantly short (2 -15 s) waves and wave-induced processes, contains strongly non-linear aspects. This non-linearity is obvious from, for example, the increasingly non-sinusoidal shape of short waves as they approach the shore through the shoaling and breaking zone (e.g., [START_REF] Elgar | Observations of bispectra of shoaling surface gravity waves[END_REF], and the generation of infragravity waves (e.g., Herbers et al., 1994Herbers et al., , 1995)), oscillatory motions in the sea surface with periods of 20 to 200 s. The non-linearity in the water motion is crucial to sediment transport. Nonsinusoidal short waves are the dominant mechanism responsible for onshore sediment transport under mild weather conditions (e.g., [START_REF] O'donoghue | Flow tunnel measurements of velocities and sand flux in oscillatory sheet flow for well-sorted and graded sands[END_REF][START_REF] Van Der | Measurements of sheet flow transport in acceleration-skewed oscillatory flow and comparison with practical formulations[END_REF][START_REF] Ruessink | Observations of velocities, sand concentrations, and fluxes under velocity-asymmetric oscillatory flows[END_REF] and may also contribute to off-shore transport [START_REF] Ribberink | Sediment transport in oscillatory boundary layers in cases of rippled beds and sheet flow[END_REF][START_REF] Grasso | Sediment transport associated with morphological beach changes forced by irregular asymmetric, skewed waves[END_REF]. Infragravity waves can dominate the water motion close to the shore during more adverse weather (e.g., [START_REF] Guza | Swash oscillations on a natural beach[END_REF][START_REF] Ruessink | Observations of swash under highly dissipative conditions[END_REF][START_REF] Ruggiero | Wave run-up on a high-energy dissipative beach[END_REF][START_REF] Sénéchal | Wave runup during extreme storm conditions[END_REF][START_REF] De Bakker | Shoreline dissipation of infragravity waves[END_REF] and are important to beach and dune erosion (e.g., [START_REF] Russell | Mechanisms for beach erosion during storms[END_REF][START_REF] Thiel De Vries | Analysis of dune erosion processes in large-scale flume experiments[END_REF].

Wave non-linearity is of paramount importance on both steep and low-sloping beaches, but has different characteristics in each case. While on steep beaches very strong non-linear processes happen locally and on a short time scale (as demonstrated by narrow surf zones with plunging or collapsing breakers), the non-linearity on low-sloping beaches is significant because it builds up over a long period of time along a cross-shore extensive area. This second case of strong non-linearity is not well documented; most of our understanding and predictive ability of nearshore processes have resulted from theoretical work, laboratory experiments, and field studies on beaches with slopes steeper than about 1:40. Within this context, a team including researchers from 9 European institutes in France, Spain, Portugal, United Kingdom and the Netherlands, as well as an institute from Chile, obtained funding in the framework of the Hydralab IV programme to perform the GLOBEX project -Gently sLOping Beach Experimentin the 100-m long, small-scale Scheldegoot (Scheldt Wave Flume) of Deltares, Delft, The Netherlands.

The overall aim of the project was to collect a high-resolution (in space and time) data set of the crossshore evolution of short and infragravity waves on a 1:80 concrete beach for a range of wave conditions, focusing on infragravity-wave dynamics, wave propagation, and boundary-layer dynamics. A detailed account of the GLOBEX project, including the description of the experimental set-up, the experimental programme, and the instrumentation used, has been presented at Coastal Dynamics 2013 (see [START_REF] Ruessink | GLOBEX: Wave dynamics on a gently sloping laboratory beach[END_REF]. A general overview of the cross-shore transformation of various wave properties was also provided. Additional preliminary results have been presented at Coastal Dynamics 2013 covering the various WPs in GLOBEX [START_REF] Almar | On the use of the Radon transform in studying wave dynamics[END_REF][START_REF] De Bakker | Infragravity-wave dissipation on a low-sloping beach[END_REF][START_REF] Rocha | Nonlinearities of short and long waves across the shoaling, surf and swash zones: large-scale physical model results[END_REF][START_REF] Ruju | Transient swash motions on a gently sloping beach[END_REF][START_REF] Tissier | Role of short and long-wave interaction on wave celerity in the surf zone of a low-sloping beach[END_REF][START_REF] Van Der | Wave boundary layer dynamics on a low sloping laboratory beach[END_REF]. In this paper, after recalling briefly the experimental procedure, we present various aspects of data processing applied on a bichromatic wave condition.

METHODS

Experimental set-up

The laboratory experiments were performed in the Scheldegoot (Scheldt Wave Flume) of Deltares, Delft, The Netherlands in April 2012. The flume is 110 m long, 1 m wide and 1.2 m high; it has glass windows along most of its length. The waves were generated with a piston-type wave maker equipped with Active Reflection Compensation (ARC) to minimize reflections from the wave paddle. An impermeable concrete beach with a 1:80 slope was constructed, with its toe at 16.57 m from the wave maker (in rest position, x = 0 m), see Figure 1. All experiments were run with a still water depth of 0.85 m over the horizontal section; this implies that the still-water shoreline was at x = 84.57 m. The material that was laying loose on the concrete bed before the flume was filled with water had a median grain size D50 of 0.75 mm; D10 and D90 amounted to 0.49 and 1.15 mm, respectively.

Figure 1. Bed elevation Z versus cross-shore distance x in the Scheldegoot during the GLOBEX project. Here, x = 0 is the wave-maker location, while Z = 0 corresponds to the still water level. The symbols indicate the location of the measurements used in this paper: the 190 dots are the positions of the wave gauges; the green (red) circles are the positions and heights above the bed of the EMCMs (ADV); the triangle is the location of the LDA. The bed level was determined with a terrestrial laser scanner before the flume was filled with water.

Experimental programme

The programme comprised 8 wave conditions grouped in 3 test series with varying complexity (Table 1). All wave-paddle steering signals in Series A and B were made with second-order wave generation.

In Series A, the signals were based on a JONSWAP spectrum with a peak enhancement factor γ of 3.3 for A1 and A2, and of 20 for A3 to provide a narrow-banded (in frequency) spectrum for the swell waves. The bichromatic frequencies in Series B were chosen in order to get the same waves in each wave-packet (6 waves per packet for B1, 10 waves per packet for B2 and B3). The packets of condition B3 had a larger modulation compared to that of B1 and B2. Since the wave height and steepness were small, the use of first-order wave generation was considered adequate for C1. The signals of A and B had a total duration of 75 minutes; the signal of C1 was 30 minutes long.

Instruments and measurements

A wide suite of instruments was deployed during the experiments:  Waves -Water-surface elevation and flow velocities were measured using 21 wave gauges (WGs) and 5 electromagnetic current meters (EMCMs) mounted on movable trolleys. All data from the WGs and EMCMs were collected at 128 Hz. The height of the EMCMs above the bed varied between 0.01 m and 0.3 m, depending on the location within the flume (see Figure 1). In addition, the water-surface elevation was also imaged at 20 Hz using 2 video cameras looking sideways through the glass windows and one video camera looking obliquely from above.  Swash -Instruments focusing specifically on swash motions comprised a single downwardlooking video camera (sampling at 20 Hz), a swash-wire mounted at about 1cm above the bed (sampling at 128 Hz), a RIEGL VZ-400 terrestrial laser scanner, and an acoustic Doppler velocimeter (ADV).  Boundary-layer dynamics -Equipment to measure horizontal and vertical velocities inside the wave bottom boundary layer comprised a Deltares Laser Doppler Velocimeter (LDV) and a Dantec Laser Dopper Anemometry (LDA) system at x = 60.6 m and x = 69.3 m, respectively. Only results from the LDA system are presented here. The vertical translation precision of the LDA was 0.05 mm. The output of a WG not included in the aforementioned 21 WGs was used to measure sea-surface elevation at the location of the LDA. Silver-coated hollow-glass spheres with a diameter of 10 μm were used to seed the flow, except during the first few days of the experiment when seeding material was unavailable.

PRELIMINARY RESULTS

Various

Short and long wave propagation

Figure 2a shows a 40-second long time-space diagram of the free-surface elevation of condition B3, together with the approximate location of the uprush and backwash motion as determined with the swash wire. The low-pass filtered (cut-off at (f 1 +f 2 )/4) signal shown in Figure 2b exhibits a lowfrequency standing-wave pattern. Following [START_REF] Guza | Swash oscillations on a natural beach[END_REF], incoming and reflected waves are deduced from the analysis of co-located WGs and EMCMs (Figure 2c,d). The separation method is likely only qualitative in both breaking (x ~ 60 m) and surf zones, but it nevertheless gives a good indication of the long-wave propagation. The reflected long-wave energy at the toe of the slope is of the same order of that of the incoming bound long wave, for this high-groupiness bichromatic case. This is in marked contrast with the A series for which most of the incoming infragravity-wave energy is dissipated at the beach (see De Bakker et al., 2013).

Time-series of the free-surface elevation at the first gauge (x = 6.96 m), at the breaking point (x = 63.58 m), in the innermost part of the surf zone and in the swash zone are plotted in Figure 3. The lowfrequency modulation of 23.8 s period is seen to evolve along the beach, developing a distinct sawtooth shape in the inner surf zone. The short waves riding the long-wave crests are travelling faster than the short-wave crests in the long-wave troughs, due to the water-depth modulation. The bore formation mechanism is thus illustrated as short waves piling up at the long-wave front (see also [START_REF] Tissier | Role of short and long-wave interaction on wave celerity in the surf zone of a low-sloping beach[END_REF]. Subsequent long-wave dissipation is likely induced by breaking at the steep front [START_REF] Van Dongeren | Shoaling and shoreline dissipation of low-frequency waves[END_REF].

The wave spectrum evolution along the beach is shown in Figure 4. Before reaching the sloping part (at x = 16.57 m), the wave energy is mainly concentrated in the two first-order components of the bichromatic, f 1 and f 2 , and in three second-order terms: the bound long wave at (f 2 -f 1 ) and the higher harmonics (2f 1 ) and (f 1 +f 2 ). During shoaling, the energy of higher-order harmonics increases:

(2f 1 -f 2 )
and (2f 2 -f 1 ) first, then (3f 1 -f 2 ), (3f 1 ) and (2f 1 +f 2 ). These are produced by interactions between first-and second-order waves [START_REF] Elgar | Observations of bispectra of shoaling surface gravity waves[END_REF]. Higher-order harmonics become increasingly important before the break point (x ~ 60 m). In the breaking zone, the main components are the first to lose their energy. Notably, in the inner surf zone (x ~ 75 m), high-order harmonics tend to dominate first-order ones.

Figure 3. Free-surface elevation of the wave group at various cross-shore positions: black lines are ensemble average over 173 groups and the gray shading indicates the standard deviation. The red lines show the propagation at the linear wave speed.

3.2

Comparison with numerical modeling One of the aims of the project is to provide a data set that can be used for comparison with numerical models. [START_REF] Ruju | Numerical analysis of run-up oscillations under dissipative conditions[END_REF] have presented SWASH numerical results focusing on run-up oscillations for the A series. The Boussinesq-type Serr1D model [START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part I: model development and analysis[END_REF][START_REF] Cienfuegos | A fourth-order compact finite volume scheme for fully nonlinear and weakly dispersive Boussinesq-type equations. Part II: boundary conditions and validation[END_REF] has been shown to reproduce correctly wave non-linearities for random waves over complex beach topographies [START_REF] Michallet | Kinematics of waves propagating and breaking on a barred beach[END_REF]. Here it is used for comparison with case B3. Time-series in the breaking, inner surf and swash zones are shown in Figure 5. The overall bore formation is correctly reproduced by the model. The short wave phases are well captured in the surf zone while discrepancies in height may appear for individual waves, leading to larger differences in the swash zone. Noteworthy, the computed velocities in the swash zone compare well with the ADV measurements when the water depth is larger than 1 cm. Figure 4. Logarithm of energy density E d as a function of frequency f and cross-shore position x for condition B3. The energy is concentrated at the multiples of (f 2 -f 1 ) for which the lines have been thickened into stripes for readability. Serr1D relies on several parameters for computing breaking inception and further dissipation (see [START_REF] Cienfuegos | A wave-breaking model for Boussinesq-type equations including roller effects in the mass conservation equation[END_REF]. The cross-shore evolution of statistical quantities is plotted for two different sets of parameters in Figure 6. The model is easily reproducing some general trends, including setdown and set-up, breaking location, overall energy decay, long wave build-up, dissipation and reflection. In particular, the standard deviation of the low-pass filtered signals exhibits a standingwave pattern with nodes close to zero, thus confirming that the reflected free long wave has about the same amplitude as the incoming bound long wave in this case.

The wave skewness and asymmetry in the surf zone are more sensitive to the choice of the parameters, such as the limiting front slope for breaking initiation and termination. This emphasizes that a good description of individual wave breaking is crucial for reproducing a correct phasing of the short-wave crests (see bottom of Figure 5) that may thus influence wave non-linearity variations (see the four bottom plots in Figure 6). This also highlights the need for investigating further the individual wave propagation characteristics, including a correct description of the individual wave celerity [START_REF] Tissier | Field measurements and nonlinear prediction of wave celerity in the surf zone[END_REF]2013).

As shown on the right side of Figure 6, the trends for depth-integrated velocities computed by the model are similar to those of the velocities measured at one depth (close to the bed in general, see Figure 1). Noteworthy, the experimental data of the mean velocity are very scattered (top-right in Figure 6). (Not shown here, the mean cross-shore velocities are much less scattered for the random wave cases.) This scatter might be attributed to the presence of secondary motions generated by the transverse waves observed for series B (see also the discussion in [START_REF] Van Der | Wave boundary layer dynamics on a low sloping laboratory beach[END_REF]. These spurious waves are likely generated at the break point, at a depth (h b ≈ 0.35 m) where the transverse mode-1seiching can be excited at frequency (f 1 +f 2 ). Figure 6. Cross-shore variation of statistical quantities for surface displacement (left) and velocities (right): measurements (dots) vs Serr1D results for two sets of breaking parameters (red and green lines). From top to bottom: mean, standard deviation (low-pass filtered are darker dots and thicker lines), skewness, asymmetry.

Bispectral analysis

In order to better understand the cross-shore evolution of the wave non-linear properties, a bispectrum analysis was conducted. The method is described in [START_REF] Elgar | Observations of bispectra of shoaling surface gravity waves[END_REF] and is not recalled here. It is aimed at finding the three-wave couplings defined by frequency triads (f m , f n , f m +f n ) with n ≤ m. The sum over all triads of the real (imaginary) part of the bispectrum is equal to the wave skewness (asymmetry). The components that constitute wave skewness at four cross-shore locations are shown in Figure 7 along with the wave spectrum and the rate of dissipation of the Fourier components. Largely positive (hot colors) or negative (cold) values of Sk m,n indicate strong triad interactions, while hot (cold) colors in the top panels indicate gain (loss) of energy.

In the shoaling zone (Figure 7a), the triad at (m=10, n=10) corresponding to (f 1 , f 1 , 2f 1 ) is the main constituent of the wave skewness, with noticeable secondary positive contributions at (f 2 , f 1 , f 1 + f 2 ) and (2f 1 , f 1 , 3f 1 ) and a negative contribution at (f 1 , f 2 -f 1 , f 2 ). In terms of spectral evolution, the main feature seen at this stage is a decay of the main components f 1 and f 2 and an increase of the higher harmonics 2f 1 and 3f 1 . After the break point, the number of triads contributing to Sk considerably increases as shown in Figure 7b. The main bispectrum component (circled) is now at (f 2 , f 1 , f 1 +f 2 ), although it still represents less than 10% of the total Sk at this stage (Sk = 1.8 at x = 67 m in Figure 6). The outer surf zone is also characterized by an energy decay of the main frequency components (e.g. f 1 , f 2 , 2f 1 ). A notable energy increase that is not linked to any triad interaction is seen at frequencies corresponding to m = 25 and 35. The energy transfer to these frequencies is thus likely due to four-wave interactions. In the inner surf zone (Figure 7c) most of the wave components are efficiently interacting, with the prevalence of interactions involving low-frequency waves at n < 10, i.e. with at least one of the frequencies of each triad being smaller than f 1 . The peaks in the spectrum are more evenly distributed in the highfrequency domain (m > 10) and a notable energy increase in the low-frequency range is observed at m = 1 and m = 8. Finally, interactions involving only low-frequency waves (m+n < 10) are prevailing in the swash zone (Figure 7d), and the main peaks at m = 1, 2, 3 are seen to lose their energy.

As an attempt to summarize these observations, the frequency domain is separated into several subdomains. The contributions of the different triads are integrated over each of these sub-domains.

f m ≥ f 1 , f n < f 1 , f m +f n ), (green) only high frequencies (f m ≥ f 1 , f n ≥ f 1 , f m +f n ), (red) only low frequencies (f m < f 1 , f n < f 1 , f m +f n < f 1 ), ( light

Bottom boundary layer measurements

The procedure used for the study of the boundary layer has been extensively presented in Van der A et al. (2013), along with the estimation of the bed roughness and results on the vertical variation of wave non-linear parameters. As an example, Figure 9f shows for the B3 wave group that the asymmetry observed outside the boundary layer is transferred into skewness close to the bed. This is a key process for explaining how steep-front waves may promote sediment transport [START_REF] Henderson | Nearshore sandbar migration predicted by an eddy-diffusivity boundary layer model[END_REF][START_REF] Abreu | Bed shear stress under skewed and asymmetric oscillatory flows[END_REF][START_REF] Berni | Surf zone cross-shore boundary layer velocity asymmetry and skewness: an experimental study on a mobile bed[END_REF].

As the reproduction of identical bichromatic wave groups was achieved, the turbulent velocity fluctuations can be extracted through ensemble averaging, as shown by the density spectra in Figure 9a

CONCLUSION

The GLOBEX experiment took place over an approximately 1 month period in April 2012. Almost 25 scientists from 10 institutes were involved, and 13 staff from Deltares performed technical support.

The experiments resulted in a highly detailed data set of the cross-shore evolution of short and infragravity waves from well seaward of the surf zone up to and including the leading edge of the swash motion on a fixed 1:80 sloping beach. In addition, high-resolution observations of the horizontal and vertical flow in the wave bottom boundary layer were collected. Examples of data analysis were briefly presented here for a single bichromatic wave condition with the aim of illustrating the wealth of processes that can be studied with the GLOBEX data set.

  results were discussed in De Bakker et al. (2013) on infragravity-wave dissipation; Rocha et al. (2013) on non-linearity of short and long waves; Almar et al. (2013) and Tissier et al. (2013) on short-wave celerity; Van der A et al. (2013) on wave boundary layer dynamics; and Ruju et al. (2013) on swash motions. Here, measurements obtained with the WGs and EMCMs are used to provide some features of the observed cross-shore evolution of the short and infragravity waves in series B3. In addition, boundary layer measurements for the same wave condition will be shown in section 3.4.

Figure 2 .

 2 Figure 2. Forty-second time-space diagram of the free-surface elevation measured by the wave gauges during condition B3: total signal (a), low-pass filtered (b), incoming (c), and reflected (d) lowfrequency waves. Warm (red) colours correspond to values above 0 (the still water level), cold (blue) colours to values below 0.The gray line fluctuating around x ~ 87 m is the shore-line determined with the swash wire.

Figure 5 .

 5 Figure 5. Time series of surface elevation (left) and velocities (right) in (from top to bottom) the breaking (x = 63.2 m), inner surf (x = 79.5 m) and swash zones (x = 85.4 m): model results (green lines) vs measurements (grey dots). Model velocities are depth integrated while velocities were measured 1 cm above the bed.

Figure 7 .

 7 Figure 7. Spectral and bispectral analysis at four cross-shore locations: shoaling (a), outer surf (b), inner surf (c), and swash (d) zones. In each panel: (top) amplitude of the Fourier components (solid line) and their variations with x (in color) against m = f m /(f 2 -f 1 ) and (bottom) real part of the bispectrum against m and n, where n = f n /(f 2 -f 1 ) with (f m , f n , f m +f n ) the frequency triad considered (only triads with large both bicoherence and bispectrum amplitude are shown). Note that, for the present case B3, f 1 /(f 2f 1 )=10 and f 2 /(f 2 -f 1 )=11.

  Figure 8. (a) Standard deviation of free surface elevation (solid line), rms of Fourier (magenta) main components f 1 and f 2 , (yellow) components larger than f 2 , (orange) components smaller than f 1 . (b) Wave skewness (solid line), sum of the real part of the bispectrum of triads involving (blue) only one low frequency (f m ≥ f 1 , f n < f 1 , f m +f n ), (green) only high frequencies (f m ≥ f 1 , f n ≥ f 1 , f m +f n ), (red) only low frequencies (f m < f 1 , f n < f 1 , f m +f n < f 1 ), (light blue) two low frequencies and one high frequency (f m < f 1 , f n < f 1 , f m +f n ≥ f 1 ). (c) Same as in (b) for the wave asymmetry. (d) Corresponding averaged biphase.

  Figure 8. (a) Standard deviation of free surface elevation (solid line), rms of Fourier (magenta) main components f 1 and f 2 , (yellow) components larger than f 2 , (orange) components smaller than f 1 . (b) Wave skewness (solid line), sum of the real part of the bispectrum of triads involving (blue) only one low frequency (f m ≥ f 1 , f n < f 1 , f m +f n ), (green) only high frequencies (f m ≥ f 1 , f n ≥ f 1 , f m +f n ), (red) only low frequencies (f m < f 1 , f n < f 1 , f m +f n < f 1 ), (light blue) two low frequencies and one high frequency (f m < f 1 , f n < f 1 , f m +f n ≥ f 1 ). (c) Same as in (b) for the wave asymmetry. (d) Corresponding averaged biphase.

  ,b. Since most of the peaks at the main frequencies (around 0.4 Hz) disappear through the procedure, low-frequency variations are notable, particularly outside the boundary layer at z = 10 mm. The vertical profiles of the fluctuations are shown in Figure9d. While the magnitude of the vertical component is relatively uniform over the 1-cm range investigated, the horizontal component tends to decrease towards the bed. Noteworthy, the standard deviation of the horizontal velocity follows a loglaw trend for elevations in the range 0.3 -1 mm (Figure9c). The measurements in the first 0.3 mm above the bed are likely affected by local roughness. Finally, Figure9eillustrates the on-shore mean current known as boundary-layer streaming[START_REF] Longuet-Higgins | Mass transport in water waves[END_REF][START_REF] Yu | Sediment transport under wave groups: relative importance between nonlinear wave shape and nonlinear boundary layer streaming[END_REF][START_REF] Kranenburg | Sand transport beneath waves: The role of progressive wave streaming and other free surface effects[END_REF]

Figure 9 .

 9 Figure 9. Example of boundary-layer measurements: (a) spectral density of horizontal velocities at three elevations above the bed and (b) with the ensemble average signal removed (the dashed line is the -5/3 slope); vertical profiles of: (c) standard deviation of horizontal and vertical velocities, (d) their fluctuations, (e) mean values; (f) skewness and asymmetry of u. The blank symbols stand for measurements made without seeding material.

Table 1 .

 1 Overview of wave conditions

	Test series A
	H s (m) T p (s)			Remark
	A1 0.10	1.58			JONSWAP, γ=3.3; prototype: H Remark
	B1 0.09	0.01	6/15	7/15	1/(f 2 -f 1 ) = 15 s;	(f 1 +f 2 )/2 = 0.433 Hz
	B2 0.09	0.01	0.42	0.462	1/(f 2 -f 1 ) = 23.81 s; (f 1 +f 2 )/2 = 0.441 Hz
	B3 0.07	0.03	0.42	0.462	1/(f 2 -f 1 ) = 23.81 s; (f 1 +f 2 )/2 = 0.441 Hz
	Test Series C
	H (m)	T (s)			Remark
	C1 0.02	23.8
	C2 0.20	2.25

s = 2 m; T p = 7 s A2 0.20 2.25 JONSWAP, γ=3.3; prototype: H s = 4 m; T p = 10 s A3 0.10 2.25 JONSWAP, γ=20; prototype: H s = 2 m; T p = 10 s Test Series B a 1 (m) a 2 (m) f 1 (Hz) f 2 (Hz) Ran for the wave bottom boundary layer study only

ACKNOWLEDGEMENT

We are greatly indebted to the Deltares technicians and staff members (Frank van Eeden, John Coolegem, Pieter Pasterkamp, Frans de Vreede, Richard Boele, Rob Hoffman, Jaap de Schipper, Paul Meys, Job Waaijerink, Hans Deutekom, Jelle Molenaar, Ivo Wenneker and Irv Elshoff) for their excellent support during the experiments. It was a pleasure to work with you! We also wish to thank all other GLOBEX team members (Céline Berni, Marcel van Maarseveen, Philippe Bonneton, Javier L. Lara, Giovanni Coco, and Iñigo Losada Rodriguez) for their help in preparing the original proposal or collecting the GLOBEX data set, and for making the stay in Delft such a pleasant and memorable occasion. This work has been supported by European Community's Seventh Framework Programme through the grant to the budget of the Integrating Activity HYDRALAB IV within the Transnational Access Activities, Contract no. 261520.


Samuel Vaiter
email: vaiter@cmap.polytechnique.fr

Gabriel Peyré
email: peyre@ceremade.dauphine.fr

Jalal Fadili
email: jalal.fadili@greyc.ensicaen.fr

Model Consistency of Partly Smooth Regularizers

Keywords: Regularization, regression, inverse problems, model consistency, partial smoothness, sensitivity analysis, sparsity, low-rank. 1

This paper studies least-square regression penalized with partly smooth convex regularizers. This class of penalty functions is very large and versatile, and allows to promote solutions conforming to some notion of low-complexity. Indeed, such penalties/regularizers force the corresponding solutions to belong to a low-dimensional manifold (the so-called model) which remains stable when the penalty function undergoes small perturbations. Such a good sensitivity property is crucial to make the underlying low-complexity (manifold) model robust to small noise. In a deterministic setting, we show that a generalized "irrepresentable condition" implies stable model selection under small noise perturbations in the observations and the design matrix, when the regularization parameter is tuned proportionally to the noise level. We also prove that this condition is almost necessary for stable model recovery. We then turn to the random setting where the design matrix and the noise are random, and the number of observations grows large. We show that under our generalized "irrepresentable condition", and a proper scaling of the regularization parameter, the regularized estimator is model consistent. In plain words, with a probability tending to one as the number of measurements tends to infinity, the regularized estimator belongs to the correct low-dimensional model manifold. This work unifies and generalizes a large body of literature, where model consistency was known to hold, for instance for the Lasso, group Lasso, total variation (fused Lasso) and nuclear/trace norm regularizers. We show that under the deterministic model selection conditions, the forward-backward proximal splitting algorithm used to solve the penalized least-square regression problem, is guaranteed to identifiy the model manifold after a finite number of iterations. Lastly, we detail how our results extend from the quadratic loss to an arbitrary smooth and strictly convex loss function. We illustrate the usefulness of our results on the problem of low-rank matrix recovery from random measurements using nuclear norm minimization.

Introduction

Problem Statement

We consider the following observation model y = Xβ 0 + w, where X ∈ R n×p is the design matrix (in statistics or machine learning) or the forward operator (in signal and imaging sciences), β 0 ∈ R p is the vector to recover and w ∈ R n is the noise. The design can be either deterministic or random, and similarly for the noise w.

Regularization is now a central theme in many fields including statistics, machine learning and inverse problems. It allows one to impose on the set of candidate solutions some prior structure on the object x 0 to be estimated. We therefore consider a proper, lowersemicontinuous (lsc) and convex function J : R p → R ∪ {+∞} to promote such a prior. Without loss of generality, we also assume that J is non-negative. This then leads to solving the following convex optimization problem min

β∈R p J(β) + 1 2λ ||Xβ -y|| 2 , (1)
where λ > 0 is the so-called regularization parameter used to balance the amount of regularization and loss.

To simplify the notations, we introduce the following "canonical" parameters

θ = (µ, u, Γ) = λ n , X * y n , X * X n ∈ Θ = R + × R p × R p×p
and we denote ε = X * w n = u -Γβ 0 .

In the following, we assume that y ∈ Im(X) and thus u ∈ Im(Γ). Obviously, this does not entail any loss of generality, as the loss term can always be written as 1 2λ ||Xβ -P Im(X) y|| 2 + 1 2λ || P Im(X) ⊥ y|| 2 , where P T the orthogonal projection on T . With these new parameters, the original problem (1) now reads min

β∈R p E(β, θ) = J(β) + 1 2µ Γβ, β - 1 µ β, u + 1 2µ Γ + u, u . (P θ )
where A + stands for the Moore-Penrose pseudo-inverse of a matrix A. With these notations, E is a function on R p × Θ. We also consider the constrained problem min β∈R p {E(β, θ 0 ) = J(β) + ι Hu (β)} where H u = {β ∈ R p ; Γβ = u} (P θ 0 )

where θ 0 = (0, u, Γ), and ι C is the indicator function of the non-empty closed convex set C, i.e. ι C (β) = 0 if β ∈ C and ι C (β) = +∞ otherwise. Problem (P θ 0 ) can be viewed as a limit of (P θ ) as µ → 0 + .

At this stage, it is worth mentioning that though we focus here, for simplicity of exposition, on the squared loss x → 1 2 ||y -Xx|| 2 , our results generalize to more general smooth losses, see Section 3.4 for further details.

The goal of this paper is to assess the recovery performance of (P θ ), i.e. to understand how close are the properties of the recovered solution of (P θ ) to those of β 0 . More precisely, we focus here on the low-noise regime, i.e. when ε is small enough, and we investigate stability in 2 sense, but also, and more importantly, the identifiability of the correct lowdimensional manifold associated to β 0 . This unifies and extend a large body of literature, including sparsity and low-rank regularization, which turn to be a very special case of the powerful theory of partly smooth regularization.

Notations

We recall some basic ingredients from differential geometry that are essential to our exposition. A good source on smooth manifold theory is [START_REF] Lee | Smooth manifolds[END_REF]. A set M ⊂ R p is a C 2 -smooth manifold around a point β ∈ R p , if β ∈ M and M consists locally around β of the solutions of some C 2 -smooth equations with linearly independent gradients. In this case, the tangent space of M at β is denoted T β (M). We define the tangent model subspace as

T β = VectHull(∂J(β)) ⊥ ,
where VectHull(C) = ρ(β -β ) ; (β, β ) ∈ C 2 , ρ ∈ R is the linear hull of the convex set C ⊂ R p . We denote ri(C) (resp. rbd(C))) the relative interior of C (resp. relative boundary), i.e. its interior (boundary) for the topology of its affine hull (the smallest affine space containing C).

For a function J, dom(J) = {β ∈ R p ; J(β) < +∞} is its domain. We denote ∂J(β) the subdifferential at β of the proper, lsc and convex function J. Geometrically, when β ∈ dom(J), ∂J(β) (if non-empty) is the set of gradients of the affine minorants of J supporting it at β. The subdifferential ∂J(β) is a closed convex set.

For a linear space T , we denote P T the orthogonal projection on T and for a matrix Γ ∈ R p×p , Γ T = P T Γ P T .

Partly-smooth Functions

Toward the goal of studying the recovery guarantees of problem (P θ ), our central assumption will be that J is a partly smooth function. Partial smoothness of functions was originally defined by Lewis (2003a). Our definition hereafter specializes it to the case of proper, lsc and convex functions.

Definition 1 Let J be a proper, lsc convex function, and

β ∈ R p such that ∂J(β) = ∅. J is partly smooth at β relative to a set M containing β if (i) (Smoothness) M is a C 2 -manifold around β and J restricted to M is C 2 around β. (ii) (Sharpness) The tangent space T β (M) is T β . (iii) (Continuity) The set-valued mapping ∂J is continuous at β relative to M.
J is said to be partly smooth relative to a set M if M is a manifold and J is partly smooth at each point β ∈ M relative to M. J is said to be locally partly smooth at β relative to a set M if M is a manifold and there exists a neighbourhood U of β such that J is partly smooth at each point β ∈ M ∩ U relative to M.

Note that in the previous definition, M needs only to be defined locally around β, and it can be shown to be locally unique thanks to prox-regularity of proper, lsc and convex functions, see (Hare and Lewis, 2004, Corollary 4.2).

Loosely speaking, a partly smooth function behaves smoothly as we move on the identifiable manifold, and sharply if we move normal to the manifold.

Remark 2 (Discussion of the properties) Since J is proper, lsc and convex, it is subdifferentially regular at any point in its domain, and in particular at β. Therefore, the Clarke regularity property (Lewis, 2003a, Definition 2.7(ii)) is automatically verified. In view of (Lewis, 2003a, Proposition 2.4(i)-(iii)), the sharpness property ((ii)) is equivalent to (Lewis, 2003a, Definition 2.7(iii)). The continuity property ((iii)) is equivalent to the fact that ∂J is inner semicontinuous at β relative to M, that is: for any sequence β n in M converging to β and any η ∈ ∂J(β), there exists a sequence of subgradients η n ∈ ∂J(β n ) converging to η. This equivalent characterization will be very useful in the proof of our main result.

Examples in Imaging and Machine Learning

We describe below some popular examples of partly smooth regularizers that are routinely used in machine learning, statistics, signal and image processing. We first expose basic building blocks (sparsity, group sparsity) and then show how the machinery of partial smoothness enables a powerful calculus to create new priors (using post-composition with a linear operator, spectral lifting, positive linear combinations and separable priors).

1 sparsity. One of the most popular non-quadratic convex regularization is the 1 norm J(β) = p i=1 |β i |, which promotes sparsity. Indeed, it is easy to check that J is partly smooth at β relative to the subspace

M = T β = {u ∈ R p ; supp(u) ⊆ supp(β)} .
The use of sparse regularizations has been popularized in the signal processing literature under the name basis pursuit method [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] and in the statistics literature under the name Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF].

1 -2 group sparsity. To better capture the sparsity pattern of natural signals and images, it is useful to structure the sparsity into non-overlapping blocks/groups B such that b∈B b = {1, . . . , p}. This group structure is enforced by using typically the mixed 1 -2 norm J(β) = b∈B ||β b ||, where β b = (β i ) i∈b ∈ R |b| . We refer to [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]Bach, 2008b) and references therein for more details. Unlike the 1 norm, and except the case |b| = 1, the 1 -2 norm is not polyhedral, but is still partly smooth at β relative to the linear manifold defined as

M = T β = β ; supp B (β ) ⊆ supp B (β)
where supp

B (β) = {b ; β b = 0} .
Spectral functions. The natural spectral extension of sparsity to matrix-valued data β ∈ R p 0 ×p 0 (where p = p 2 0 ) is to impose a low-rank prior, which should be understood as sparsity of the singular values. Denote β = U β diag(Λ β )V * β an SVD decomposition of β, where Λ β ∈ R p 0 + . The nuclear norm is defined as

J(β) = ||β|| * = ||Λ β || 1 . (2)
It has been used for instance in machine learning applications (Bach, 2008b), matrix completion [START_REF] Recht | Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization[END_REF][START_REF] Candès | Exact matrix completion via convex optimization[END_REF] and phase retrieval (Candès et al., 2013).

The nuclear norm can be shown to be partly smooth at x relative to the manifold [START_REF] Lewis | Alternating projections on manifolds[END_REF], Example 2)

M = β ; rank(β ) = rank(β) . (3)
More generally, if j : R p 0 → R is a permutation-invariant closed convex function, then one can consider the function J(β) = j(Λ β ) which can be shown to be a convex function as well (Lewis, 2003b). When restricted to the linear space of symmetric matrices, j is partly smooth at Λ β for a manifold m Λ β , if and only if J is partly smooth at β relative to the manifold

M = U diag(Λ)U * ; Λ ∈ m Λ β , U ∈ O p 0 ,
where O p 0 ⊂ R p 0 ×p 0 is the group of orthogonal matrices. This result is proved in (Daniilidis et al., 2014a, Theorem 3.19), building upon the work of (Daniilidis et al., 2014b) on manifold smoothness transfer under spectral lifting. This result can be extended to non-symmetric (possibly rectangular) matrices by requiring that j is an absolutely permutation-invariant closed convex function, see (Daniilidis et al., 2014a, Theorem 5.3). The nuclear norm || • || * is a special case where j(Λ) = ||Λ|| 1 .

Analysis regularizers. If J 0 : R q → R ∪ {+∞} is a proper lsc convex function and D ∈ R p×q is a linear operator, an analysis regularizer (following the terminology introduced in [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]) is of the form

J(β) = J 0 (D * β).
Such a prior controls the low-complexity (as measured by J 0 ) of the correlations between the columns of D and β. A popular example is when taking J 0 = || • || 1 and D * a finite-difference approximation of the gradient of an image. This defines the (anisotropic) total variation, which promotes piecewise constant images, and is popular in image processing [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. The fused Lasso [START_REF] Tibshirani | Sparsity and smoothness via the fused Lasso[END_REF] corresponds to J 0 being the 1 -norm and D * is the concatenation of the identity and a finite-difference operator. To cope with correlated covariates in linear regression, it was devised in [START_REF] Grave | Trace lasso: a trace norm regularization for correlated designs[END_REF][START_REF] Richard | Intersecting singularities for multi-structured estimation[END_REF] to use a family of analysis-type priors where J 0 = || • || * is the nuclear norm. If J 0 is partly smooth at α = D * β for the manifold M 0 α , then it is shown in (Lewis, 2003a, Theorem 4.2) that J is partly smooth at β relative to the manifold

M = β ∈ R p ; D * β ∈ M 0
α . provided that the following transversality condition holds (Lee, 2003, Theorem 6.30(a))

Ker(D) ∩ T α (M 0 α ) ⊥ = {0} ⇐⇒ Im(D * ) + T α (M 0 α ) = R N . Moreover, the co-dimension of M in R p equals the co-dimension of M 0 α in R q .
Mixed regularization. Starting from a family of proper, lsc and convex functions {J } ∈L , L = {1, . . . , L}, it is possible to design a convex function as J (β) = ∈L ρ J (β), where ρ > 0 are weights. A popular example is to impose both sparsity and low rank of a matrix, by using J 1 = || • || 1 and J 2 = || • || * , see for instance [START_REF] Golbabaee | Hyperspectral Image Compressed Sensing Via Low-Rank And Joint-Sparse Matrix Recovery[END_REF][START_REF] Oymak | Simultaneously structured models with application to sparse and low-rank matrices[END_REF].

Suppose that ∈L ri(dom(J )) = ∅. Let S ⊆ R p be a C 2 -manifold. If each J is partly smooth at β relative to a submanifold M ⊆ S, then it can be shown that J is also partly smooth at

β for M = ∈L M ,
with the proviso that the submanifolds M intersect transversally, i.e.

∈L z = 0 and z ∈ T β (M ) ⊥ for each ∈ L ⇒ z ∈ T β (S) ⊥ for each ∈ L .
Moreover, the co-dimension of M (in S) equals the sum of the co-dimensions of M . This assertion is a weaker version of (Lewis, 2003a, Corollary 4.8), since we use convexity and closedness of the functions J . For the case where L = 2, the above transversality condition reads (Lee, 2003, Theorem 6.30(b))

T β (M 1 ) ⊥ ∩ T β (M 2 ) ⊥ = T β (S) ⊥ ⇐⇒ T β (M 1 ) + T β (M 2 ) = T β (S) . (4)
Separable Regularization. Let {J } ∈L , L = {1, . . . , L}, be a family of proper lsc convex functions. If J is partly smooth at β relative to a manifold M β , then the separable function (Lewis, 2003a, Proposition 4.5).

J ({β } ∈L ) = ∈L J (β ) is partly smooth at (β 1 , . . . , β L ) relative to M 1 β 1 × • • • × M L β L
One fundamental problem that has attracted a lot of interest in the recent years in data processing involves decomposing an observed object into a linear combination of components/constituents β , ∈ L. One instance of such a problem is image decomposition into texture and piece-wise-smooth (cartoon) parts, see e.g. [START_REF] Starck | Image decomposition via the combination of sparse representatntions and variational approach[END_REF][START_REF] Aujol | Image decomposition into a bounded variation component and an oscillating component[END_REF][START_REF] Peyré | Learning the morphological diversity[END_REF] and references therein. Another example of decomposition is principal component pursuit, proposed in [START_REF] Candès | Robust principal component analysis?[END_REF], to decompose a matrix which is the superposition of a low-rank component and a sparse component. In this case

J 1 = || • || 1 and J 2 = || • || * .

Main results

In the following, we denote T = T β 0 , e = P T (∂J(β 0 )) ∈ R p . Before stating our main contributions, we first introduce a central object of this paper, which controls the stability of M when the signal to noise ratio is large enough.

Definition 3 (Linearized pre-certificate) For some matrix Γ ∈ R p×p , assuming ker(Γ)∩ T = {0}, we define η Γ = ΓΓ + T e.

In Section 5.2, we will investigate the connection between η Γ and the so-called minimal norm certificate.

3.1 Deterministic model consistency.

We first consider the case where X and w (or equivalently Γ and u) are fixed and deterministic. Our main contribution is the following theorem, which shows the robustness of the manifold M associated to β 0 to small perturbations on both the observations and the design matrix, provided that µ (or equivalently λ) is well chosen.

Theorem 4 Assume that J is locally partly smooth at β 0 relative to M and that there exists Γ ∈ R p×p such that ker( Γ) ∩ T = {0}, and η Γ ∈ ri(∂J(β 0 )).

(5)

Then, there exists a constant C > 0 such that if

max ||Γ -Γ||, ||ε||µ -1 , µ C, (6)
the solution β θ of (P θ ) is unique and satisfies

β θ ∈ M and ||β θ -β 0 || = O(||ε||). ( 7
)
This theorem is proved in Section 5.3.

Remark 5 (Inverse problems) A typical case of application of this result is in inverse problems that are encountered in various disciplines in science and engineering, such as in signal and image processing. In such a setting, the forward operator X is generally fixed and known, and one then takes Γ = Γ = X * X/n.

Remark 6 (Uncertain design/forward operator) If only a noisy version of the forward operator (in inverse problems) or the design (in regression) is available then this can also be handled by Theorem 4. This scenario has been considered for sparse recovery (i.e. J the 1 -norm) by several authors for sparse linear regression and compressed sensing, see e.g. [START_REF] Herman | General deviants: An analysis of perturbations in compressed sensing[END_REF][START_REF] Rosenbaum | Sparse recovery under matrix uncertainty[END_REF][START_REF] Loh | High-dimensional regression with noisy and missing data: Provable guarantees with nonconvexity[END_REF].

Remark 7 (Random setting) In statistics or machine learning, one considers a regression problem where the design X and the noise w are random, under the asymptotic regime where the number of observations n, i.e. number of rows X, grows large, so that Γ only reach Γ in the limit n → +∞. See Theorem 11 for details.

Remark 8 (Identification of the manifold) Theorem 4 guarantees that, under some hypotheses on β 0 and θ, β θ belongs to M. For all the regularizations considered in Section 2.1, it turns out that actually M β θ = M. This is because, for any

(β, β ) with β ∈ M β close enough to β, one has M β = M β .
The following proposition, proved in Section 5.6, shows that Theorem 4 is in some sense sharp, since the hypothesis η Γ ∈ ri(∂J(β 0 )) (almost) characterizes the stability of M.

Proposition 9 Suppose that β 0 is the unique solution of (P (0, Γβ 0 , Γ) ) and that ker( Γ) ∩ T = {0}, and η Γ / ∈ ∂J(β 0 ). ( 8)

Then there exists C > 0 such that if (6) holds, then any solution β θ of (P θ ) for µ > 0 satisfies β θ / ∈ M.

In the particular case where ε = 0 (no noise) and Γ = Γ, this result shows that the manifold M cannot be correctly identified by solving P (µ,Γβ 0 ,Γ) for any µ > 0 small enough.

Remark 10 (Critical case) The only case not covered by either Theorem 4 or Proposition 9 is when η Γ ∈ rbd(∂J(β 0 )). In this case, one cannot conclude in general, since depending on the noise w, one can have either stability or non-stability of M. We refer to (Vaiter et al., 2013b) where an example illustrates this situation for the 1-D total variation J = ||D * DIF • || 1 (here D * DIF is a discretization of the 1-D derivative operator).

3.2 Probabilistic model consistency.

We now turn to study consistency of our estimator. In this section, we work under the classical setting where p and β 0 are fixed as the number of observations n → ∞. We consider that the design matrix and the noise are random. More precisely, the data

(ξ i , w i ) are random vectors in R p × R, i = 1, • • • , n
, where ξ i is the i-th row of X, are assumed independent and identically distributed (i.i.d.) samples from a joint probability distribution such that E (w i |ξ i ) = 0, finite fourth-order moments, i.e. E w 4 i < +∞ and E ||ξ i || 4 < +∞. Note that in general, w i and ξ i are not necessarily independent. It is possible to extend our result to other distribution models by weakening some of the assumptions and strenghthening others, see e.g. [START_REF] Knight | Asymptotics for Lasso-Type Estimators[END_REF][START_REF] Zhao | On model selection consistency of Lasso[END_REF]Bach, 2008b). Let's denote Γ = E(ξ * ξ) ∈ R p×p , where ξ is any row of X. We do not make any assumption on invertibility of Γ.

To make the discussion clearer, the canonical parameters θ will be indexed by n. The estimator β θn obtained by solving (P θn ) for a sequence θ n is said to be consistent for β 0 if, lim n→+∞ Pr (β θn is unique) → 1 and β θn converges to β 0 in probability. The estimator is said to be model consistent if lim n→+∞ Pr (β θn ∈ M) → 1, where M is the manifold associated to β 0 .

The following result ensures model consistency for certain scaling of µ n . It is proved in Section 5.5

Theorem 11 If conditions (5) hold and

µ n = o(1) and µ -1 n = o(n 1/2 ). ( 9
)
Then the estimator β θn of β 0 obtained by solving (P θn ) is model consistent.

Remark 12 (Sharpness of the criterion) One can also state a probabilistic equivalent to Proposition 8. That is, if β 0 is the unique solution of (P 0, Γβ 0 , Γ), and conditions (8) and (9) hold, then the estimator β θn of β 0 defined by solving (P θn ) cannot be model consistent.

Algorithmic Implications

A popular iterative scheme to compute a solution of (P θ ) is the Forward-Backward splitting algorithm. A comprehensive treatment of the convergence properties of this algorithm, and other proximal splitting schemes, can be found in the monograph [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. Starting from some β 0 ∈ R p , the algorithm implements the following iteration

β k+1 = Prox τ k µJ (β k + τ k (u -Γβ k )) ,
where the step size satisfies 0 < τ τ k τ < 2/||Γ||, and the proximity operator is defined, for γ > 0, as

Prox γJ (β) = argmin β ∈R p 1 2 ||β -β || 2 + γJ(β ).
The following theorem shows that the Forward-Backward algorithm correctly identifies the manifold M after a finite number of iterations.

Theorem 13 Suppose that the assumptions of Theorem 4 hold. Then, there exists k 0 large enough, such that for all k k 0 , the Forward-Backward iterates satisfy β k ∈ M.

Proof Inspection of the proof of Theorem 4 shows that the solution β θ of (P θ ), which is unique, is such that the vector η θ = u-Γβ θ λ satisfies η θ ∈ ri(∂J(β θ )) when ( 5) and ( 6) hold. Moreover, as β θ ∈ M, β θ is near β 0 , and J is locally partly smooth at β 0 relative to M, it is also partly smooth at β θ relative to the same manifold M. Altogether, this implies that the assumptions of (Liang et al., 2014, Theorem 3.1) are fulfilled and the manifold identification claim follows.

This result sheds light on the convergence behaviour of this algorithm in the favourable case where condition (5) holds and (||Γ -Γ||, ||ε||/µ, µ) are sufficiently small.

General Loss Functions

For the sake of simplicity, we have described our contributions with the squared loss function u ∈ R n → 1 2 ||y -u|| 2 . Our results, however, extend readily to the case of more general loss functions of the form F (u, y). In the following ∇ 2 1 F (u, y) denotes the Hessian of F with respect to the first variable evaluated at (u, y).

We thus consider the variational problem

min β F (Xβ, y) + λJ(β) ,
where the loss function F : R n × R n → R fulfills the following assumptions:

(A.1) For any y ∈ R n , F (•, y) ∈ C 2 (R n ) is σ m -strongly convex and ∇ 1 F (•, y) is σ M -Lipschitz continuous, σ m > 0 and σ M > 0. (A.
2) The gradient of F with respect to the first variable,

∇ 1 F (u, y), is such that ∇ 1 F (u, u) = 0.
Any loss function of the form F (u, y) = G(u) -u, y , where G is a C 2 strongly convex function and and its gradient is Lipschitz-continuous, satisfies assumptions (A.1). Assumption (A.2) is quite natural for a data fidelity term, and is fulfilled for instance for some losses in the exponential family. In this setting, Theorem 4 (and in a similar way our other contributions) remains valid, and one simply needs to replace condition (6) by

max || Γ -Γ||, ||ε||µ -1 , µ C, (10)
where now

Γ = 1 n X * ∇ 2 1 F (y, y)X and ε = 1 n X * ∇ 2 1 F (y, y)w ,
where ∇ 2 1 F (y, y) is the Hessian with respect to the first variable (assumed to be positive definite by assumption (A.1)) taken at (y, y). A detailed treatment on the way to adapt the proofs to handle such a generic loss is provided in Section 5.4.

Relation to Previous Works

Works on linear convergence rates. Following the pioneer work [START_REF] Burger | Convergence rates of convex variational regularization[END_REF]) (who study convergence in terms of Bregman divergence), there is a large amount of works on the study conditions under which ||β θ -β 0 || = O(||ε||) (so-called linear convergence rate) where β θ is any solution of (P θ ), see for instance the book [START_REF] Scherzer | Variational Methods in Imaging[END_REF] for an overview of these results. The initial work of [START_REF] Grasmair | Necessary and sufficient conditions for linear convergence of 1 -regularization[END_REF] proves a sharp criteria to ensure linear convergence rate for the 1 norm, and this approach is further extended to arbitrary convex functions by [START_REF] Grasmair | Linear convergence rates for Tikhonov regularization with positively homogeneous functionals[END_REF] and [START_REF] Fadili | Stable recovery with analysis decomposable priors[END_REF], who respectively proved linear convergence rates in terms of the penalty J and 2 -norm.

These works show that if ker(Γ) ∩ T = {0} and ∃η ∈ Im(Γ) ∩ ri(∂J(β 0 ))

(which is often called the source condition), then linear convergence rate holds. Note that condition (5) implies ( 11), but it is stronger. Indeed, condition (11) does not ensure model consistency (7), which is a stronger requirement. Model consistency requires, as we show in our work, the use of a special certificate, the minimal norm certificate ηΓ , which is equal to

η Γ if η Γ ∈ ri(∂J(β 0 )) (see Proposition 17).
Works on model consistency. Theorem 4 is a generalization of a large body of results in the literature. For the Lasso, i.e. J = ||•|| 1 , and when Γ = Γ, to the best of our knowledge, this result was initially stated in [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF]. In this setting, the result (7) corresponds to the correct identification of the support, i.e. supp(β θ ) = supp(β 0 ). Condition (5) for J = || • || 1 is known in the statistics literature under the name "irrepresentable condition", see e.g. [START_REF] Zhao | On model selection consistency of Lasso[END_REF]. [START_REF] Knight | Asymptotics for Lasso-Type Estimators[END_REF] have shown estimation consistency for Lasso for fixed p and β 0 and asymptotic normality of the estimates. The authors in [START_REF] Zhao | On model selection consistency of Lasso[END_REF] proved Theorem 11 for J = || • || 1 , though under slightly different assumptions on the covariance and noise distribution. A similar result was established in [START_REF] Jia | On model selection consistency of the elastic net when p n[END_REF] for the elastic net, i.e. J = || • || 1 + ρ|| • || 2 2 for ρ > 0. In (Bach, 2008a) and (Bach, 2008b), the author has shown Theorem 11 for two special cases, namely the group Lasso nuclear/trace norm minimization. Note that these previous works assume that the asymptotic covariance Γ is invertible. We do not impose such an assumption, and only require the weaker restricted injectivity condition ker( Γ) ∩ T = {0}. In a previous work (Vaiter et al., 2013b), we have proved an instance of Theorem 4 when Γ = Γ and J(β) = ||D * β|| 1 , where D ∈ R p×q is an arbitrary linear operator. This covers as special cases the discrete anisotropic total variation or the fused Lasso. This result was further generalized in (Vaiter et al., 2013a) when Γ = Γ, and J belongs to the class of partly smooth functions relative to affine manifolds M, i.e. M = β + T β . Typical instances encompassed in this class are the 1 -2 norm, or its analysis version, as well as nonnegative polyhedral functions including the ∞ norm. Note that the nuclear norm (and composition of it with linear operators as studied for instance in [START_REF] Grave | Trace lasso: a trace norm regularization for correlated designs[END_REF][START_REF] Richard | Intersecting singularities for multi-structured estimation[END_REF]), whose manifold is not affine, does not fit into the framework of (Vaiter et al., 2013a), while it is covered by Theorem 4. J. D. [START_REF] Lee | On model selection consistency of m-estimators with geometrically decomposable penalties[END_REF] investigated a class of geometrically decomposable penalties, for which they extended the irrepresentable condition and used it to establish 2 -consistency and model consistency. This class of penalties turns out to be a very special case of ours. Lastly, a similar result was proved in [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF] for an infinite dimensional sparse recovery problem over space of Radon measures, when J is the total variation of a measure (not to be confused with the total variation semi-norm mentioned above). In this setting, an interesting finding is that, when ηX * X ∈ ri(∂J(β 0 )), ηX * X is not equal to η X * X but to a different certificate (called "vanishing derivative" certificate by [START_REF] Duval | Exact support recovery for sparse spikes deconvolution[END_REF]) that can also be computed by solving a linear system.

Compressed sensing Condition (5) is often used when X is drawn from the Gaussian matrix ensemble to asses the performance of compressed sensing recovery with 1 norm, see [START_REF] Wainwright | Sharp thresholds for high-dimensional and noisy sparsity recovery using 1 -constrained quadratic programming (lasso)[END_REF][START_REF] Dossal | Sharp support recovery from noisy random measurements by l1 minimization[END_REF]. This is extended to a more general family of decomposable norms (including in particular 1 -2 norms and the nuclear norm) in [START_REF] Candès | Simple bounds for recovering low-complexity models[END_REF], but only in the noiseless setting. Our result shows that this analysis extends to the noisy setting as well, and ensures model consistency at high signal to low noise levels. The same condition is used to asses the performance of matrix completion (i.e. the operator X is a random masking operator) in a noiseless setting [START_REF] Candès | Exact matrix completion via convex optimization[END_REF][START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF]. It was also used to ensure 2 robustness of matrix completion in a noisy setting [START_REF] Candès | Matrix completion with noise[END_REF], and our findings shows that these results also ensure rank consistency for matrix completion at high signal to low noise levels.

Sensitivity analysis. Sensitivity analysis is a central theme in variational analysis. Theorems 4 can be understood as a sensitivity analysis of the minimizers of E at the point (β 0 , θ 0 ). Classical sensitivity analysis of non-smooth optimization problems seeks conditions to ensure smoothness of the mapping θ → β θ where β θ is a minimizer of f (•, θ), see for instance [START_REF] Mordukhovich | Sensitivity analysis in nonsmooth optimization. Theoretical Aspects of Industrial Design[END_REF][START_REF] Rockafellar | Variational analysis[END_REF][START_REF] Bonnans | Perturbation analysis of optimization problems[END_REF]. This is usually guaranteed by the non-degenerate source condition and restricted injectivity condition (11), which, as already reviewed above, ensures linear convergence rate, and hence Lipschitz behaviour of this mapping. The result captured by Theorem 4 goes one step further, by assessing that M β 0 is a stable manifold (in the sense of [START_REF] Wright | Identifiable surfaces in constrained optimization[END_REF]), since the minimizer β θ is unique and remains in M β 0 for θ close to θ 0 . Our starting point for establishing Theorem 4 is the inspiring work of Lewis (2003a) who first introduced the notion of partial smoothness and showed that this broad class of functions enjoys a powerful calculus and sensitivity theory. For convex functions (which is the setting considered in our work), partial smoothness is closely related to U -V-decompositions developed in [START_REF] Lemaréchal | The u-lagrangian of a convex function[END_REF]. In fact, the behaviour of a partly smooth function and of its minimizers (or critical points) depend essentially on its restriction to this manifold, hence offering a powerful framework for sensitivity analysis theory. In particular, critical points of partly smooth functions move stably on the manifold as the function undergoes small perturbations (Lewis, 2003a;[START_REF] Lewis | Partial smoothness, tilt stability, and generalized hessians[END_REF]. A important and distinctive feature of Theorem 4 is that, partial smoothness of J at β 0 relative to M transfers to E(•, θ) for λ > 0, but not when λ = 0 in general. In particular, (Lewis, 2003a, Theorem 5.7) does not apply to prove our claim.

Case Study: Nuclear Norm Regularization

In this section, we illustrate the usefulness of our model consistency results to derive a sharp manifold stability analysis for the nuclear norm (a.k.a trace norm) regularization. As detailed in Section 3.5, previous consistency results due to Bach (2008b) only apply to the overdetermined setting, while our result tackles arbitrary design X by only requiring the weaker injectivity condition (5). For simplicity of exposition, we consider recovery of square matrices of size p = p 0 × p 0 , but the same holds for arbitrary rectangular matrices.

Irrepresentability Criterion IC

The nuclear norm, defined in (2), turns out to be the tightest convex relation of the rank function on the spectral ball. It is then the best convex candidate to enforce a low-rank prior [START_REF] Fazel | Matrix rank minimization with applications[END_REF]. It is moreover partly smooth at any β 0 ∈ R p 0 ×p 0 relative to the manifold M of fixed rank r = rank(β 0 ) defined in (3).

Let β = U diag(Λ)V * be a reduced rank-r SVD decomposition of β 0 , with V, U ∈ R p 0 ×r with orthonormal columns and Λ ∈ (R * + ) r . The subdifferential of the nuclear norm at β 0 reads (see for instance [START_REF] Candès | Simple bounds for recovering low-complexity models[END_REF])

∂|| • || * (β 0 ) = η ∈ R p 0 ×p 0 ; η T = e and ||η S || 1 , ( 12
)
where ||η|| is the operator norm, T = T β (M), S = T ⊥ and e = P T (∂J(β 0 )), with

T = U A * + BV * ; A, B ∈ R p 0 ×r and e = U V *
and S is the subspace of matrices spanned by the family (wz * ), where w (resp. z) is any vector orthogonal to U (resp. V ).

The relative interior of ∂|| • || * (β 0 ) is formed by subgradients η satisfying the inequality in (12) strictly. Thus, condition (5) in the case Γ = Γ takes the analytical form

η Γ ∈ ri(∂J(β 0 )) ⇐⇒ IC(β 0 ) < 1 where IC(β 0 ) = || P S ΓΓ + T e|| . ( 13
)
The value of IC(β 0 ) can then be easily computed. Loosely speaking, the smaller the quantity 1-IC(β 0 ) is, the further η Γ is from the relative boundary of ∂J(β 0 ), and in turn the smaller the stability constant controlling ||β θ -β 0 ||/||ε|| in ( 7) is.

Recovery from Gaussian Measurements

Bounding IC for an arbitrary operator X and matrix β 0 is in general difficult. It is however possible to leverage tools from random matrix theory to obtain sharp upperbounds when X is drawn from certain matrix ensembles. This strategy has been deployed to study matrix completion problems, see for instance [START_REF] Candès | Exact matrix completion via convex optimization[END_REF]; [START_REF] Candès | The power of convex relaxation: Near-optimal matrix completion[END_REF]. Another problem on which we now focus is when X is drawn from the standard Gaussian ensemble, i.e. its entries are independent identically distributed from N (0, 1). The following result, proved by [START_REF] Candès | Simple bounds for recovering low-complexity models[END_REF], shows that IC(β 0 ) < 1 with high probability as soon as n is larger than 6rp 0 (up to negligible terms).

Proposition 14 ( [START_REF] Candès | Simple bounds for recovering low-complexity models[END_REF]

, Theorem 1.2) Let β 0 ∈ R p 0 ×p 0 such that rank(β 0 ) = r. If n δr(6p 0 -5r) ( 14
)
for some δ > 1, then IC(β 0 ) < 1 with probability at least 1 -2e (1-δ)p 0 /8 .

Combining this result with Theorem 4, this shows that under the scaling ( 14) of (n, p 0 , r), one obtains with high probability on the design matrix a rank-consistent estimation of the unknown matrix β 0 , which is (to the best of our knowledge) a novel result.

Figure 1 illustrates this result by computing the average (over 25 Monte Carlo replications) values of IC(β 0 ) for either a varying n or rank r. The shaded area corresponds to ±3× standard deviation across the 25 replications, and the dashed vertical line indicates the transition predicted by ( 14). This suggests numerically that the upper-bound ( 14) is indeed sharp.

Proofs

Uniqueness sufficient condition

Proposition 15 Let J be a proper lsc convex function. For a point β, assume that ker(Γ) ∩ T β = {0}, and η Γ ∈ ri(∂J(β)).

Then β is the unique minimizer of (P θ ) (resp. (P 0,Γβ,Γ )).

Proof This is a consequence of (Vaiter et al., 2013a, Corollary 1). Though their result was stated for J finite-valued convex, it remains valid when it is proper lsc and convex. Indeed, in this case, J is subdifferentially regular at β (Rockafellar and Wets, 1998, Example 7.27). Moreover, ∂J(β) = ∅ by assumption, and thus the directional derivative at β is proper, sublinear and closed, and it is the support of ∂J(β) (Rockafellar and Wets, 1998, Theorem 8.30). Continuing the proof as in (Vaiter et al., 2013a, Corollary 1) shows the claim.

Minimal norm certificate and linearized pre-certificate

In this section, we establish the connections between the minimal norm certificate and the linearized pre-certificate η Γ. Recall that J is proper lsc convex function, and we suppose that Im( Γ) ∩ ∂J(β 0 ) = ∅ (so-called range or source condition in the inverse problem 14). (b) IC(β 0 ) as a function of r for a fixed n = 0.6p 2 0 . The vertical dashed line shows the threshold r = (3p 0 -9p 2 0 -5n)/5 indicated by ( 14). community). The latter is equivalent to the fact that β 0 is a minimizer of (P θ 0 ). This is straightforward to see by writing the first-order optimality condition of this convex program.

Definition 16 (Minimal norm certificate) The minimal norm certificate is the vector ηΓ = Γz Γ, where zΓ = argmin

Γz∈∂J(β 0 ) ||z||. ( 15
)
This certificate is uniquely defined as the constraint set is non-empty closed and convex, and the solution of the minimization problem, which is the projection of the origin on it, is obviously unique.

Proposition 17 Assume that ker( Γ) ∩ T = {0}. Then,

η Γ ∈ ri(∂J(β 0 )) =⇒ ηΓ = η Γ, ( 16
) ηΓ ∈ ri(∂J(β 0 )) =⇒ ηΓ = η Γ. ( 17
)
Under either of these conditions, β 0 is the unique minimizer to (P 0, Γβ 0 , Γ).

Proof Proof of (16) Under condition ker(X) ∩ T = {0}, we have from the definition of Γ+ T , that

z Γ = Γ+ T e = argmin z ||z|| subject to ΓT z = e (18)
and thus η Γ = Γz Γ .

Clearly, the constraint set of problem (18) includes that of (15), which entails

||z Γ|| ||z Γ|| .
If η Γ ∈ ri(∂J(β 0 )), then z Γ is also a feasible point of problem ( 15) and thus

||z Γ|| ||z Γ|| .
Altogether, we get that ||z Γ|| = ||z Γ|| and, since zΓ is the unique minimizer of (15), we get that zΓ = z Γ, which implies that ηΓ = η Γ.

Proof of (17) Let S = T ⊥ . Problem ( 15) can be conviniently rewritten as zΓ = argmin z ||z|| subject to ΓT z = e and ΓS z ∈ P S (∂J(β 0 )) .

The fact that ηΓ = Γz Γ ∈ ri(∂J(β 0 )) implies P S ηΓ = P S Γz Γ ∈ ri(P S ∂J(β 0 )), and thus, the second constraint in the last problem is inactive. We then recover problem (18), which in turn implies that ηΓ = η Γ.

Proof of uniqueness. See Proposition 15.

Proof of Theorem 4

In order to prove Theorem 4, we consider any sequence

θ k = (µ k , u k = Γ k x 0 + ε k , Γ k ) k where X k ∈ R n k ×p . Assume that Γ k , ε k µ -1 k , µ k -→ ( Γ, 0, 0) . (19)
Then proving Theorem 4 boils down to showing that for k large enough, the solution β k of (P θ k ) is unique and satisfies β k ∈ M.

Constrained problem. We consider the following non-smooth, in general non-convex, constrained minimization problem

β k ∈ Argmin β∈M∩K E(β, θ k ) ( 20
)
where K is an arbitrary fixed convex compact neighbourhood of β 0 .

The following key lemma establishes the convergence of β k to β 0 .

Lemma 18 Under conditions (5) and (19), β k → β 0 .

Proof We denote ||u|| 2 Γ = Γu, u for any positive semidefinite matrix Γ. Under condition (5), Proposition 15 implies that β 0 is the unique solution of (P 0, Γβ 0 , Γ). By optimality of β k one has E(β k , θ k ) E(β 0 , θ k ) and hence

1 2 ||β k || 2 Γ k -β k , Γ k β 0 + ε k + µ k J(β k ) 1 2 ||β 0 || 2 Γ k -β 0 , Γ k β 0 + ε k + µ k J(β 0 )
which is equivalently stated as

1 2 ||β k -β 0 || 2 Γ k -β k -β 0 , ε k + µ k J(β k ) µ k J(β 0 ). ( 21
)
Since β k ∈ K, the sequence (β k ) k is bounded, and we let β be any accumulation point. Using ( 19), that J is non-negative and J(β k ) are bounded, we have lim sup

k→∞ (µ k J(β k )) lim k→∞ µ k lim sup k→∞ J(β k ) = 0 and lim inf k→∞ (µ k J(β k )) lim k→∞ µ k lim inf k→∞ J(β k ) = J(β ) lim k→∞ µ k = 0 ,
and thus lim k→∞ (µ k J(β k )) = 0. Consequently, passing to the limit in (21), using ( 19), and continuity of the inner product and the norm, shows that ||β -β 0 || 2 Γ 0, or equivalently Γβ = Γβ 0 , i.e. β is a feasible point of (P 0, Γx 0 , Γ). Furthermore, since

1 2 ||β k -β 0 || 2 Γ k 0, (21) yields -β k -β 0 , ε k µ k + J(β k ) J(β 0 ).
Passing again to the limit, using lower semicontinuity of J, (19) and continuity of the inner product, we then get

J(β ) lim inf k→∞ J(β k ) = lim inf k→∞ -β k -β 0 , ε k µ k + J(β k ) lim sup k→∞ -β k -β 0 , ε k µ k + J(β k ) = lim sup k→∞ J(β k ) J(β 0 ) .
Combining this with the previous claim on feasibility of β for (P 0, Γx 0 , Γ) allows to conclude that β is a solution of (P 0, Γx 0 , Γ). Since β 0 is unique, this leads to β = β 0 .

We now aim at showing that for k large enough, β k is the unique solution of (P θ k ).

Convergence of the tangent model subspace. By definition of the constrained problem (20), β k ∈ M. Moreover, since E(•, θ k ) is partly smooth at β 0 relative to M, the sharpness property Definition 1((ii)) holds at all nearby points in the manifold M (Lewis, 2003a, Proposition 2.10). Thus as soon as k is large enough, we have T k = T β k (M). Using the fact that M is of class C 2 , we get

T k = T β k (M) -→ T β 0 (M) = T (22)
when ( 19) holds, where the convergence should be understood over the Grassmannian of linear subspaces with the same dimension (or equivalently, as the convergence of the projection operators P T k → P T ). Since ker( Γ) ∩ T = {0}, (22) implies that for k large enough, when (19) holds, ker

(Γ k ) ∩ T k = {0}, (23)
which we assume from now on.

First order condition. Let B be the Euclidean unit ball in R n . Take K = β 0 + rB for r > 0 sufficiently large. For any δ > 0, ∃k δ such that ∀k k δ , β k ∈ β 0 + δB according to Lemma 18. Thus, for k large enough, i.e. δ sufficiently small, we indeed have β k ∈ int(K). Furthermore, it is easy to see that ι K is locally partly smooth at β 0 relative to K, and thus is partly smooth at β k relative to K for k large enough. Moreover, local partial smoothness of J at β 0 relative to M entails that J is also partly smooth at β k relative to M. Therefore, the sum rule (Lewis, 2003a, Corollary 4.6) (the transversality condition is satisfied as K is full-dimensional and β k ∈ int(K), see ( 4)) shows that, for all sufficiently large k, J + ι K is locally partly smooth at β k relative to M ∩ K, and then so is E(•, θ k ) + ι K by the smooth perturbation rule (Lewis, 2003a, Corollary 4.7). Therefore, (Lewis, 2003a, Proposition 2.4(a)-(b)) applies, and it follows that β k is a critical point of ( 20) if, and only if,

0 ∈ Aff(∂E(β k , θ k ) + N K (β k )) = Γ k β k -u k µ k + Aff(∂J(β k )) = Γ k β k -u k µ k + e β k + T k ⊥ .
The first equality comes from the fact that E(•, θ) is a closed convex function, and that the normal cone of K at β k vanishes on the interior points of K, and the second one from the decomposability of the subdifferential. Projecting this relation onto T k , we get, since

e β k ∈ T k , P T k (Γ k β k -u k ) + µ k e β k = 0. ( 24
)
Convergence of the primal variables. Since both β k and β 0 belong to M, and partial smoothness implies that M is a manifold of class C 2 around each of them, we deduce that each point in their respective neighbourhoods has a unique projection on M Poliquin et al. (2000). In particular, β k = P M (β k ) and β 0 = P M (β 0 ). Moreover, P M is of class C 1 near β k (Lewis and Malick, 2008, Lemma 4). Thus, C 2 differentiability shows that

β k -β 0 = P M (β k ) -P M (β 0 ) = D P M (β k )(β k -β 0 ) + R(β k ) where R(β k ) = O(||β k -β 0 || 2
) and where D P M (β k ) is the derivative of P M at β k . Using (Lewis and Malick, 2008, Lemma 4), and recalling that T k = T β k (M) by the sharpness property, the derivative D P M (β k ) is given by D P M (β k ) = P T k . Inserting this in (24), we get

P T k Γ k (P T k (β k -β 0 ) + R(β k )) -P T k ε k + µ k e β k = 0. (25)
Using ( 23), Γ k,T k has full rank, and thus

β k -β 0 = Γ + k,T k (ε k -µ k e β k -Γ k R(β k )) , (26)
where we also used that

T k ⊥ ⊂ ker(Γ + k,T k ). One has Γ + k,T k → Γ+ T so that Γ + k,T k Γ k = O(1) and Γ + k,T k = O(1)
. Altogether, we thus obtain the bound

||β k -β 0 || = O (||ε k ||, µ k ) . ( 27
)
Convergence of the dual variables. We define

η k = u k -Γ k β k µ k
. Arguing as above, and using (26) we have

µ k η k = ε k + Γ k (β 0 -β k ) = ε k -Γ k Γ + k,T k (ε k -µ k e β k -Γ k R(β k )) = ε k -Γ k P T k Γ + k,T k (ε k -µ k e β k -Γ k R(β k )) = P V ⊥ T k ε k + P V T k Γ k R(β k ) + µ k Γ k Γ + k,T k e β k ,
where we denoted V T k = Im(Γ k P T k ), and used that Im(Γ + k,T k ) ⊂ T k . We thus arrive at

||η k -η Γ|| = O ||ε k ||µ -1 k , ||Γ k Γ + k,T k e β k -η Γ||, ||Γ k ||||β k -β 0 || 2 µ -1 k .
Since M is a C 2 manifold, and by partial smoothness (J is

C 2 on M), we have β → e β is C 1 on M, one has ||e β k -e|| = O(||β k -β 0 ||). (28)
Using the triangle inequality, we get

||Γ k Γ + k,T k -ΓΓ + T || ||Γ + k,T k ||||Γ k -Γ|| + || Γ||||Γ + k,T k -Γ+ T ||. Again, since Γ + k,T k → Γ+ T , we have ||Γ + k,T k || = O(1)
. Moreover, A → A + is smooth at A = Γ T along the manifold of matrices of constant rank, and M is a C 2 manifold near β 0 . Thus

||Γ + k,T k -Γ+ T || = O(||Γ k,T k -ΓT ||) = O(||Γ k -Γ||, || P T k -P T ||) = O(||Γ k -Γ||, ||β k -β 0 ||).
This shows that

||Γ k Γ + k,T k -ΓΓ + T || = O(||Γ k -Γ||, ||β k -β 0 ||). (29)
Putting ( 28) and ( 29) together implies

||Γ k Γ + k,T k e β k -η Γ|| = O(||Γ k -Γ||, ||β k -β 0 ||).
Altogether, we get the bound

||η k -η Γ|| = O ||ε k ||µ -1 k , ||β k -β 0 ||, ||Γ k -Γ||, ||Γ k ||||β k -β 0 || 2 µ -1 k .
Since ||β k -β 0 || is bounded according to (27), we arrive at

||η k -η Γ|| = O ||Γ k -Γ||, ||ε k ||µ -1 k , µ k . (30)
Convergence inside the relative interior. Using the hypothesis that η Γ ∈ ri(∂J(β 0 )), we will show that for k large enough,

η k ∈ ri(∂J(β k )). (31)
Let us suppose this does not hold. Then there exists a sub-sequence of η k , that we do not relabel for the sake of readability of the proof, such that

η k ∈ rbd(∂J(β k )) . (32)
According to (30) and Lemma 18, under (19), (β k , η k ) → (β 0 , η Γ). Condition (32) is equivalently stated as, for each k

∃z k ∈ T ⊥ β k , ∀ η ∈ ∂J(β k ), z k , η -η k 0, (33)
where one can impose the normalization ||z k || = 1 by positive-homogeneity. Up to a subsequence (that for simplicity we still denote z k with a slight abuse of notation), since z k is in a compact set, we can assume z k approaches a non-zero cluster point z .

Since

T ⊥ β k → T ⊥ because M is a C 2 manifold, one has that z ∈ T ⊥ . We now show that ∀ v ∈ ∂J(β 0 ), z , η -η Γ 0. (34)
Indeed, let us consider any v ∈ ∂J(β 0 ). In view of the continuity property in Definition 1((iii)) ∂J is continuous at β 0 along M, so that since β k → β 0 there exists v k ∈ ∂J(β k ) with v k → v. Applying (33) with η = v k gives z k , v k -η k 0. Taking the limit k → +∞ in this inequality leads to (34), which contradicts the fact that η Γ ∈ ri(∂J(β 0 )). In view of ( 31) and ( 23), using Proposition 15 shows that β k is the unique solution of (P θ k ).

General Loss Function

Using , all our results remain valid by appropriately adapting the proofs. We now detail the necessary arguments to adapt the proof of Theorem 4 to a generic loss function satisfying assumptions (A.1)-(A.2).

Proof of Proposition 15 It follows from assumption (A.1) that F (•, y) is strictly convex, and the uniqueness follows from (Liang et al., 2014, Theorem A.1).

Proof of Lemma 18 Problem (20) now reads

β k ∈ Argmin β∈M∩K F (X k β, y k ) + λ k J(β) . Optimality of β k entails F (X k β k , y k ) + λ k J(β k ) F (X k β 0 , y k ) + λ k J(β 0 ) .
By assumptions (A.1)-(A.2), we have the following useful inequalities for any u ∈ R n , see e.g. (Nesterov, 2004, p. 57 and 64) σ m 2 ||y -u|| 2 F (u, y) -F (y, y) = F (u, y) -F (y, y) -∇F (y, y), u -y σ M 2 ||y -u|| 2 .

It then follows that

F (X k β k , y k ) -F (X k β 0 , y k ) σ m 2 ||y k -X k β k || 2 - σ M 2 ||w k || 2
and therefore

λ k J(β 0 ) σ m 2 ||y k -X k β k || 2 - σ M 2 ||w k || 2 + λ k J(β k ) σ m 2σ M ||y k -X k β k || 2 ∇ 2 1 F (y k ,y k ) - σ M 2 ||w k || 2 + λ k J(β k ) σ m 2σ M ||β k -β 0 || 2 X * k ∇ 2 1 F (y k ,y k )X k - σ m σ M β k -β 0 , X * k ∇ 2 1 F (y k , y k )w k - σ M 2 1 - σ 2 m σ 2 M ||w k || 2 + λ k J(β k ) ,
where we used strong convexity of assumption (A.1) in the second and third inequalities. Dividing both sides by 1/P we obtain

σ m 2σ M ||β k -β 0 || 2 Γk - σ m σ M β k -β 0 , εk - σ M 2 1 - σ 2 m σ 2 M ||n -1/2 w k || 2 + µ k J(β k ) µ k J(β 0 ) ,
Let τ = dist(η Γ, ∂J(β 0 )) = inf η∈∂J(β 0 ) ||η -η Γ||. Since ∂J(β 0 ) is a non-empty, closed and convex set, the infimum is attained and one has τ > 0 since η Γ ∈ ∂J(β 0 ).

We now prove the claim by contradiction. Let β j be a solution of (P θ j ) such that (6) holds at θ j for j sufficiently large (taking C smaller if necessary so that C < τ ), and suppose that β j ∈ M. Thus, β j is also a solution of (20) for θ j , whence it follows that η j ∈ η Γ + CB. Using the triangle inequality, we then get dist(η j , ∂J(β 0 )) > τ -C > 0 .

(35)

Now, in view of the continuity property in Definition 1((iii)), we have ∂J(β k ) → ∂J(β 0 ) along M. This is equivalent, since ∂J(β 0 ) is closed and using (Rockafellar and Wets, 1998, Corollary 4.7), to dist(η, ∂J(β k )) → dist(η, ∂J(β 0 )) for every η ∈ R n , i.e.

∀δ > 0, ∃k 0 , ∀k k 0 , | dist(η, ∂J(β k )) -dist(η, ∂J(β 0 ))| < δ, ∀η ∈ R n .

In particular, as β j is a minimizer of (P θ j ) for j large enough, we have η j ∈ ∂J(β j ), and thus dist(η j , ∂J(β 0 )) < δ, leading to a contradiction with (35). Hence, β j ∈ M.

Conclusion

In this paper, we provided a very general and principled analysis of the recovery performance when partly smooth functions are used to regularize linear inverse/regression problems. This class of functions encompass all popular regularizers used in the literature. The generality of our results is unprecedented since for the first time, a unified analysis is provided together with a generalized "irrepresentable condition" to guarantee consistent identification of the low-complexity manifold underlying the original object. Our work also shows that model consistency is not only of theoretical interest, but also has algorithmic and practical consequences. Indeed, after a finite number of iterations, the iterates of the proximal splitting algorithm used to solve the original optimization problem (here the Forward-Backward), are guaranteed to lie on the original manifold. This opens the door to acceleration by switching to a higher-order smooth optimization method, exploiting the smoothness of the partly smooth objective function along the identified smooth model manifold.

Figure 1 :

 1 Figure 1: Curves of IC(β 0 ) (central solid line: average, blue shaded band: ±3× standard deviation) for 25 realizations of X ∈ R n×p 2 0 from the standard Gaussian ensemble, where p 0 = 10 3 , and a random β 0 = AB * of rank r where A, B ∈ R p 0 ×r are Gaussian matrices. (a) IC(β 0 ) as a function of n for a fixed r = 4. The vertical dashed line shows the threshold n = r(6p 0 -5r) indicated by (14). (b) IC(β 0 ) as a function of r for a fixed n = 0.6p 2 0 . The vertical dashed line shows the threshold r = (3p 0 -9p 2 0 -5n)/5 indicated by (14).

Acknowledgements

The authors would like to thank Vincent Duval and Jérôme Malick for fruitful discussions. This work has been supported by the European Research Council (ERC project SIGMA-Vision).

where now

Changing (19) to Γk , εk µ -1 k , µ k -→ ( Γ, 0, 0), which entails implicitly that n -1/2 w k → 0, and arguing as in the rest of the proof of the lemma allows to conclude that β k → β 0 .

Proof of Theorem 4 C 2 -continuity of F allows to use the smooth perturbation rule to conclude that partial smoothness of J is preserved upon adding F . Condition (24) now becomes

Dividing by n, plugging this expansion back into the above first-order (criticality) condition, and grouping the O(.) terms, condition (25) becomes

where

Then with the new notations ( Γk , εk ) in place of (Γ k , ε k ), one sees that the proof continues unchanged.

Proof of Theorem 11

It is sufficient to check that ( 6) is in force with probability 1 as n → +∞. Owing to classical results on convergence of sample covariances, which apply thanks to the assumption that the fourth order moments are finite, we get Γ n -Γ = O P n -1/2 and 1 n ξ i , w = O P n -1/2 , where we used the assumption that E ( ξ i , w ) = 0. As p is fixed, it follows that ||Γ n -Γ|| = O P n -1/2 and ||ε n || = O P n -1/2 . Thus under the scaling (9), we get

which indeed converges to 0 in probability. This concludes the proof.

Proof of Proposition 9

Let (β k ) k be a sequence of solutions to the constrained problem (20). Since β 0 is the unique minimizer to (P (0, Γβ 0 , Γ) ) and ( 5) is satisfied, η Γ is well-defined. Moreover, arguing as in the proof of Lemma 18 and Theorem 4, under condition (6), we have (β k , η k ) → (β 0 , η Γ), and η k ∈ η Γ + CB.