Kai Dang

Maureen Clerc

Clair Vandersteen

Nicolas Guevara

Dan Gnansia

In situ validation of a parametric model of electrical field distribution in an implanted cochlea

Cochlear implants have been proved to be an effective treatment for patients with sensorineural hearing loss. Among all the approaches that have been developed to design better cochlear implants, 3D model-based simulation stands out due to its detailed description of the electric field which helps reveal the electrophysiological phenomena inside the cochlea. With the advances in the cochlear implant manufacturing technology, the requirement on simulation accuracy increases. Improving the simulation accuracy relies on two aspects: 1) a better geometrical description of the cochlea that is able to distinguish the subtle differences across patients; 2) a comprehensive and reliable validation of the created 3D model. In this paper, targeting at high precision simulation, we propose a parametric cochlea model which uses micro-CT images to adapt to different cochlea geometries, then demonstrate its validation process with multichannel stimulation data measured from a implanted cochlea. Comparisons between the simulation and validation data show a good match under a variety of stimulation configurations. The results suggest that the electric field distribution is affected by the geometric characteristics of each individual cochlea. These differences can be correctly reflected by simulations based on a 3D model tuned with personalized data.

I. INTRODUCTION

As a highly successful neural prostheses, cochlear implant have been used for many years to help patients who suffer from sensorineural hearing loss recover the sense of hearing. The cochlear implant works by directly stimulating the residual auditory nerve fibers inside the cochlea according to their naturally developed tonotopic organization [START_REF] Wilson | Cochlear implants: current designs and future possibilities[END_REF], [START_REF] Zeng | Cochlear implants: system design, integration, and evaluation[END_REF]. Despite the significant improvements that have been made in the past decades, the performance of modern cochlear implants is still far behind the healthy cochlea. A suspected reason is the precision of stimulation. In a healthy cochlea, an auditory nerve fiber is activated by the hair cell to which it is connected. While in an implanted cochlea, the auditory nerve fibers are activated by electrical stimulations originating from a certain distance. The tissue and liquid between the current source and the nerve fiber diffuse the stimulating current, reducing its spatial selectivity [START_REF] Tang | Spatial channel interactions in cochlear implants[END_REF]. As a result, the excited nerve fibers spreads and the number of uncorrelated information channels decreases [START_REF] Mcdonnell | A channel model for inferring the optimal number of electrodes for future cochlear implants[END_REF].

Efforts have been made to improve spatial selectivity with optimized the stimulation strategies. Since the implants only have a limited capability of measuring electric field related data in the patients' cochleas, 3D model based simulations are often used to fill the blanks. Until now, based on the computation of the electric field [START_REF] Briaire | 3d mesh generation to solve the electrical volume conduction problem in the implanted inner ear[END_REF], [START_REF] Sibella | Modeling and computation of electric potential field distribution generated in cochlear tissues by cochlear implant stimulations[END_REF], [START_REF] Choi | A novel approach to compute the impedance matrix of a cochlear implant system incorporating an electrode-tissue interface based on finite element method[END_REF], the simulation approach has been used to optimize the current focusing strategy [START_REF] Bonham | Current focusing and steering: Modeling, physiology, and psychophysics[END_REF], to study the consequences of neural degeneration [START_REF] Briaire | The consequences of neural degeneration regarding optimal cochlear implant position in scala tympani: a model approach[END_REF], and to model the electric compound action potential (ECAP) of the auditory nerve [START_REF] Choi | Modeling ecap in cochlear implants using the fem and equivalent circuits[END_REF], [START_REF] Briaire | Unraveling the electrically evoked compound action potential[END_REF].

Besides current diffusion, the performance of the implant on each patient is also largely affected by individual factors such as the cochlea shape, position of the electrode array and the neuron survival [START_REF] Briaire | The consequences of neural degeneration regarding optimal cochlear implant position in scala tympani: a model approach[END_REF], [START_REF] Finley | Role of electrode placement as a contributor to variability in cochlear implant outcomes[END_REF]. For the generalized 3D cochlear models developed in the previous studies, the effects of these factors are hard to estimate. To address this problem and propose patient-specific 3D electric field simulation in the cochlea, we propose an adaptive cochlear model based on parametric descriptions. Using micro-CT images, this model is quickly adapted to the scanned cochlea in terms of the cochlea shape and the position of the implanted electrode array.

In previous researches, cochlear model are often validated through indirect approaches. For instance, comparing the simulated input impedance of each electrode with the impedance measurement of a real cochlea implant is a commonly used method [START_REF] Choi | A novel approach to compute the impedance matrix of a cochlear implant system incorporating an electrode-tissue interface based on finite element method[END_REF], [START_REF] Kalkman | Place pitch versus electrode location in a realistic computational model of the implanted human cochlea[END_REF]. However, since the input impedance is measured between pairs of electrodes, this validation method does not cover multi-channel stimulation, which is commonly performed by the actual implants [START_REF] Bierer | Identifying cochlear implant channels with poor electrode-neuron interface: partial tripolar, single-channel thresholds and psychophysical tuning curves[END_REF]. In this paper, we also demonstrate a novel validation approach which directly measures the electric potential during a variety of stimulation configurations. The potential data is recorded in the cochlea of an implanted human specimen with high sampling rate and amplitude resolution. The acquired data is compared to the simulation results in terms of stimulation patterns, stimulating positions and measuring positions.

II. METHOD

A. Electric field simulation

Patient-specific geometric details, such as the placement of the electrode array inside the cochlea, are closely correlated to the performance of a cochlear implant [START_REF] Finley | Role of electrode placement as a contributor to variability in cochlear implant outcomes[END_REF]. Yet the exact causes of performance variations are still not clear. Motivated by the lack of a simulation tool on this aspect, we developed a parametrical modeling approach that is able to adapt to a specified cochlea through micro-CT images.

Simulation of the electric field is equivalent to solving a Poisson equation (1) in 3D:

∇ • (σ∇ϕ) = 0 (1)
in which ϕ is the 3D electric potential distribution and σ represents the conductivity distribution of the tissues. In our simulation, equation 1 is solved via a symmetric boundary element method (BEM) approach [START_REF] Gramfort | OpenMEEG: opensource software for quasistatic bioelectromagnetics[END_REF]. The BEM is economical compared to FEM because it only requires surface meshes as geometrical input. In our parametric model, these meshes are created by rotating the cross-section of the cochlea in a cylindrical coordinate system. The output of the simulation is the value of the electric potential at a set of fixed positions within the cochlea. These virtual measurement points are organized along 200 virtual auditory nerve fibers, stretching from the base to the apex of the cochlea. Each fiber contains 50 measurement points, located from the peripheral process to the central axon, and depicted as red curves in Fig 1.

B. Geometrical adaptation

Our experimental setup involved an EVO electrode array (Neurelec-Oticon Medical) implanted in the right cochlea of a human cadaver. The EVO array has a total length of 25mm, with 20 equally spaced cylindrical platinum contacts. The implanted temporal bone was placed into a Micro-CT which provided a high resolution 3D scan of the cochlea (resolution 24.8µm in each direction). As the first step of adapting the 3D model, the same cylindrical coordinate system used by the parametric model is defined in the scan data using the methods from [START_REF] Verbist | Consensus panel on a cochlear coordinate system applicable in histological, physiological and radiological studies of the human cochlea[END_REF], [START_REF] Clark | A scalable model for human scala-tympani phantoms[END_REF]. The polar axis inside the scan is shown in Fig 2(a) as the yellow arrow, while the longitudinal axis is pointing from the centre of the modiolus to the apex of the cochlea, parallel to the view direction. 4) that create the cross-section of ST through interpolation.

C. Intra-cochlear stimulation and acquisition platform

We devised an original experiment setup for in-situ validation, by implanting by cochleostomy the left cochlea of a human cadaver. A home-designed stimulation and acquisition platform allowed to measure the intra-cochlear potential resulting from stimulation with an EVO electrode array. On a multichannel electrode array, different channels can be activated at the same time to create a desired electric field pattern inside the cochlea [START_REF] Bierer | Identifying cochlear implant channels with poor electrode-neuron interface: partial tripolar, single-channel thresholds and psychophysical tuning curves[END_REF], [START_REF] Landsberger | Reducing current spread using current focusing in cochlear implant users[END_REF]. To record this field pattern, our system combines the stimulation research platform with a multichannel recording circuit. The system diagram is shown in Fig 5 . The final ST cross-section is created from the interpolation In the stimulation part of the system, a stimulation chip is connected to a 20 channel electrode array and a reference as in the actual implant. The chip belongs to a research platform, which generates different stimulation sequences according to user's configuration. The reference electrode is placed between the skull and the scalp at the same side of the head as the implanted cochlea. It is a platinum spherical contact with a diameter of 1.5mm.

In order to measure the electric potential generated by the stimulation, the 20 wires from the electrode array are also connected to an analog multiplexer that leads to the recording part of the system. Since an open circuit is easily corrupted by interference from the environment, the output of the multiplexer is connected to the reference via a 10kΩ resistor to create a loop in the circuit. The two leads of the resistor are connected to the input of a differential probe, which keeps the stimulation part isolated from outside while extracting the differential signal to the digital oscilloscope. Finally, the digital oscilloscope samples the signal at 5MHz sampling rate and 16bits amplitude resolution. The sampled data is sent to PC for storage and analysis. During the recording process, each stimulation waveform is programmed to repeat 20 times. In each repetition, an analog multiplexer connects one of the 20 electrodes to the input of the recording part, so waveforms conducted by the cochlear tissue can be recorded at the electrode positions. In a modern cochlear implant system, the different strategies to stimulate auditory nerve fibers involve changing the waveform, the duration and intensity of stimulation and the spatial layout of the stimulation electrodes. Those variable parameters are configured as follows in the stimulation experiment:

1) Stimulation waveform: During the stimulation, each electrode can be assigned to one of the following modes: output / passive discharge or high impedance. The stimulation waveform is created by switching between the modes given above. work at different modes at the same time creates a spatial stimulation pattern. In this experiment, we measure the electric potential with the following spatial patterns:

• Monopolar. Only one electrode is assigned to output mode while the rest are set to high impedance mode.

• Bipolar. A pair of neighbouring electrodes are assigned to the output mode but one sends the opposite stimulation waveform from the other. All the other electrodes are set to high impedance.

• BP + n. Same as BP, but the two stimulating electrodes have n electrodes in between. Here we set n to 3.

• Tripolar. 3 neighbouring electrodes are assigned to output mode. The middle one sends the stimulation waveform with a +1 gain factor while its two neighbours send the same waveform at -1/2 gain.

Besides the configurations given above, due to the circuit topology given in Fig 5, the reference electrode is always assigned to passive discharge mode during simulations.

For each spatial stimulation pattern, all waveforms with the amplitude and duration configurations above will be used. In the monopolar mode, each electrode is used once as the stimulating electrode. In other modes, the stimulation is repeated at least on 2 different positions on the electrode array.

The combination of parameters given above makes a comprehensive measurement scenario which covers most of the realistic stimulation settings.

III. RESULTS AND DISCUSSION

The in-situ validation discussed here concerns the spatial current distribution. The output of the model is the electric potential simulated along 200 auditory nerve fibers. Fig 7 .  gives an example of the output when performing bipolar stimulation between the 12th and 13th electrodes.

To compare the simulation with the potential measured on the implanted electrodes, we chose the measurement point closest to the electrode array on each nerve fiber, and create a 200 point potential plot along the scala tympani. Since there is gain error between the simulation output and the measurements, each simulation output zoomed with a factor to match the amplitude of the measurements.

The global error between the simulation and the measurements can be described using the Relative Difference Measure (RDM) computed in the following way:

RDM (g s , g m ) = g s ∥g s ∥ - g m ∥g m ∥ (2)
In equation 2, g s is the simulation output while g m is the potential measured during the stimulation. The results are given in table I.

Plots of the potential distribution over different spatial stimulation patterns is shown in Fig 8 . Four different patterns are compared: monopolar, bipolar, tripolar and BP+3 . The results show a generally good match between the simulation output and the measured potential. However, the details of the error indicate that the simulation accuracy can be further improved with more accurate geometry or even conductivity adaption. For example, all the measured potentials have a small step between the 10th and the 11th electrode, suggesting the affect of an unknown geometry factor in that position. The validation process shows the importance of considering geometrical details in order to acquire higher simulation accuracy with 3D cochlear models, which was done here by adapting the model to micro-CT scan data. The validated model will provide more insight into the interaction between the stimulation pattern and the cochlear geometry, and help the design of personalized cochlear implants in the future.

Fig. 1 .

 1 Fig. 1. The parametric cochlear model with different parts marked out by color: (scala tympani(ST) and scala vestibuli(SV), light green), (nerve tissue, yellow); (electrode array, light blue); (metal electrodes, dark blue) and (nerve fibers, red). Only one out of 5 nerve fibers is plotted for image clarity.

Fig. 2 .

 2 Fig. 2. Position of the electrodes extracted from micro-CT scan (a). Yellow arrow is the polar axis of the cylindrical coordinate system, green lines are the locations of the cross-sections images shown in Fig 3. Compared with 3D electrode array model in (b). The blue and the red array each represents the position of the electrode array before and after adaption.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Vertical cross-sections of the cochlea at four angular coordinates. Green curves are the manually segmented ST and SV, separated by the estimated position of the basilar membrane (yellow line). The red curve is the apex of the cochlea where the ST and SV cannot be separated.

Fig. 5 .

 5 Fig. 5. Diagram of the stimulation -recording system.

  Fig 6 plots a stimulation waveform used by the experiment. A relatively long period of passive discharge follows the output mode in order to absorb the residual unbalanced charge caused by the error of the stimulation chip.

2 )Fig. 6 .

 26 Fig. 6. Plateau stimulation waveform. Colors indicate the work modes of the electrode: (black, high impedance); (blue, output); (green, passive discharge).

Fig. 7 .

 7 Fig. 7. 3D electric field distribution of a bipolar stimulation.

Fig. 8 .

 8 Fig. 8. 3D electric field distribution of a bipolar stimulation.

TABLE I

 I

	.	RDM OVER DIFFERENT STIMULATION PATTERNS
	Spatial pattern	Monopolar	Bipolar	Tripolar	BP+3	Average
	RDM	0.155	0.148	0.179	0.134	0.154


Guillaume Carlier
email: carlier@ceremade.dauphine.fr

Xavier Dupuis
email: xdupuis@luiss.it

An iterated projection approach to variational problems under generalized convexity constraints

Keywords: Principal-agent problem, b-convexity constraint, convexity constraint, convex envelopes, iterated projections, Dykstra's algorithm Mathematics subject classification: 49M25, 65K15, 90C25

. Our method also turns out to be simple for convex envelope computations.

Introduction

Variational problems subject to a convexity constraint arise in several different contexts such as mathematical economics [START_REF] Rochet | Ironing, sweeping, and multidimensional screening[END_REF], Newton's least resistance problem [START_REF] Brock | A symmetry problem in the calculus of variations[END_REF], [START_REF] Buttazzo | Minimum problems over sets of concave functions and related questions[END_REF], [START_REF] Lachand-Robert | An example of non-convex minimization and an application to Newton's problem of the body of least resistance[END_REF], optimal transport for the quadratic cost [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] or shape optimization [START_REF] Buttazzo | Shape optimization problems over classes of convex domains[END_REF], [START_REF] Lamboley | Regularity and singularities of optimal convex shapes in the plane[END_REF]. Existence of minimizers is generally not an issue since the set of convex functions have good local compactness properties in most reasonable functional spaces. However, on the one hand, the convexity constraint makes it difficult to write optimality conditions in the form of a tractable Euler-Lagrange PDE (see [START_REF] Lions | Identification du cône dual des fonctions convexes et applications[END_REF]) which could be used for instance to derive regularity results (see [START_REF] Carlier | Regularity of solutions for some variational problems subject to a convexity constraint[END_REF]). On the other hand, handling the convexity constraint numerically in a consistent and efficient way is also a challenging problem which has received a lot of attention in the last fifteen years. Choné and Le Meur [START_REF] Choné | Non-convergence result for conformal approximation of variational problems subject to a convexity constraint[END_REF] have first identified specific difficulties, one of which being that one cannot use conformal convex finite-elements since they form, as the meshsize vanishes, a set of functions which is not (by far!) dense in the set of convex functions. From this negative result, a lot of proposals have been made: use of interpolates of convex functions [START_REF] Carlier | A numerical approach to variational problems subject to convexity constraint[END_REF], supporting hyperplanes [START_REF] Ekeland | An algorithm for computing solutions of variational problems with global convexity constraints[END_REF], finite-differences and semi-definite programming for nonnegativity of the Hessian [START_REF] Aguilera | Approximating optimization problems over convex functions[END_REF], polyhedral approximation and directional convexity in the spirit of wide stencils for nonlinear PDEs [START_REF] Oberman | A numerical method for variational problems with convexity constraints[END_REF], splitting of the directional convexity constraints and proximal algorithms [START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF]. . . Whatever method is used, some subtle tradeoff has to be made between provable convergence, accuracy and the computational cost resulting from the number of convexity constraints enforced at the discretized level (typically O(N 2 ) with an N points grid). A major breakhtrough (for two dimensions) has been made recently by Mirebeau [START_REF] Mirebeau | Adaptive, anisotropic and hierarchical cones of discrete convex functions[END_REF] who introduced a hierarchy of subcones of the cone of interpolates of convex functions and an adaptative refinement strategy leading typically on a grid with N points to essentially only O(N ln 2 (N)) convexity constraints.

It turns out that in optimal transport (for non quadratic transportation costs) and in principal-agent problems in economics (see section 2 for a brief presentation), the case of the convexity constraint is somehow special and corresponds to a very particular choice of the transport cost (quadratic) or of the valuation function (bilinear). The general form of the principal-agent problem typically involves a (valuation) function b of two variables x and y (in R d , say) and the relevant constraint on the admissible functions is that they are b-convex (a b-convex function is a function u of the variable x which can be written as a supremum of the form u(x) = sup y {b(x, y)v(y)} so that when b is the scalar product, we recover the convexity constraint). One then has to minimize a certain integral functional among b-convex functions; general existence results can be proved but not much more in general and this has been a serious limitation for understanding principal-agent models in several dimensions. It is only recently with the pathbreaking work of Figalli, Kim and McCann [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF] that some conditions (intimately related to the conditions of Ma, Trudinger and Wang [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF] for the regularity of optimal transport maps) were identified to make the principal-agent problem a convex program and, in the first place, the set of b-convex functions a convex set (see section 4 below for an elementary presentation in a special case). This raised the hope to use convex optimization algorithms to solve numerically some convex minimization problems posed over the set of b-convex functions, provided it is convex. To the best of our knowledge, the present paper is the first one which addresses b-convexity numerically (even for a quite resticted class of b). What makes the analysis of Figalli, Kim and McCann amenable to a computational approach is that, under some conditions detailed in section 3, the fact that u is b-convex and that q = ∇u can be rewritten as a single condition:

u(x ′ ) -u(x) ≥ Γ b (x ′ , x, q(x)), ∀x ′ , x,
for a certain function Γ b depending (in a somehow indirect way) on the primitive datum b and which, provided it is convex in its last argument, immediately implies convexity of the set of b-convex functions. Now, discretizing this constraint on an N-points grid gives N 2 convex constraints. We will restrict ourselves to projection problems, i.e. quadratic minimization problems, so that, once discretized on a grid with N points, we face a projection on the intersection of convex sets, each of whom being given by one of the constraints above corresponding to two grid points x ′ and x. This can be solved with Dykstra's iterative projection algorithm (see the original papers by Dyskstra [START_REF] Dykstra | An algorithm for restricted least squares regression[END_REF] for the case of convex cones, Boyle and Dykstra [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF] for general closed convex sets; see also the extension to general Bregman projections by Bauschke and Lewis [START_REF] Bauschke | Dykstra's algorithm with Bregman projections: a convergence proof[END_REF] or to monotone operators by Bauschke and Combettes [START_REF] Bauschke | A Dykstra-like algorithm for two monotone operators[END_REF], Combettes [START_REF] Combettes | Iterative construction of the resolvent of a sum of maximal monotone operators[END_REF]). Each step of the algorithm consists in a projection onto only one of the convex sets; if these elementary projections are easy (it depends on the complexity and geometry of Γ b , note that they are explicit in the case of the convexity constraint), Dykstra's algorithm can effectively be used even if it is computationally costly.

In the case of the convexity constraint, the use of Dykstra's iterative projection algorithm is not really competitive with most of the approaches recalled above. Nevertheless, as far as we know, it is the only one which can be quite directly adapted to the case of the b-convexity constraint. It is also worth noting that our numerical approach presents some similarities with the one of Ekeland and Moreno-Bromberg [START_REF] Ekeland | An algorithm for computing solutions of variational problems with global convexity constraints[END_REF] in the sense that it replaces finite families of supporting hyperplanes by finite families of graphs of the nonlinear functions Γ b (., x j , q j ). It also shares at least in spirit some similarities with the method of Mérigot and Oudet [START_REF] Mérigot | Handling convexity-like constraints in variational problems[END_REF]: both methods somehow split the huge number of constraints into smaller subsets on which projections can be performed quickly.

The paper is organized as follows. In section 2, we motivate the study of variational problems with a b-convexity constraint by the principal-agent problem in economics. Section 3 recalls the structural conditions of Figalli, Kim and McCann under which b-convex functions form a convex set and thanks to which existence and uniqueness of a minimizer can be established. In section 4, we specify a particular class of b satisfying the conditions of Figalli, Kim and McCann, which can be seen as perturbations of the scalar product and for which we give an elementary and self-contained presentation. In section 5, we restrict ourselves even further to the case of quadratic objectives, i.e. to projection problems, and give some examples, in particular related to envelope computations. In section 6 we prove a convergence result for discretization of these projection problems. Section 7 explains how to use Dykstra's algorithm in this context. Finally, numerical results are presented in section 8.

Motivations

The principal-agent problem

To motivate variational problems with a b-convexity constraint, let us start by an informal description of the so-called principal-agent problem in economics. Let us consider a population of heterogeneous agents, each of whom has some vector of characteristics x in the sense that an agent with type x derives a utility b(x, y)p from consuming a good of type y for the price p. The set of possible agents types is denoted X, the set of possible goods types is denoted Y and b is a given function from X × Y → R. We now consider a monopolist (the principal ), proposing to the population of agents a contract menu that is a pair of functions x ∈ X → (y(x), p(x)) ∈ Y ×R specifying for each type of agent x the good and the price intended for her. The monopolist cannot observe directly the characteristics of agents, so the contract menu should be consistent with the fact that (y(x), p(x)) is really prefered by type x; this is the so-called incentive-compatibility condition:

b(x, y(x)) -p(x) ≥ b(x, y(x ′ )) -p(x ′ ), ∀x, x ′ . (2.1)
Moreover the agents have access to an outside option given by a fixed pair (y 0 , p 0 ) and thus accept a contract from the monopolist only if b(x, y(x))p(x) ≥ u 0 (x) := b(x, y 0 )p 0 , ∀x.

(2.

2)

The aim of the monopolist is then to minimize his total cost, which is given by an expression of the form X φ(x, p(x), y(x))dx, over the (y(.), p(.)) satisfying conditions (2.1) and (2.2). For instance, φ can be given by φ(x, p, y) = (c(y)p)ρ(x) where c(y) is the production cost of y and ρ represents the density of agents types, but other models are possible. Defining u(x) := b(x, y(x))-p(x) for every x, we see that the incentive-compatibility condition can be rewritten as: [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF], by assuming that the relation q = ∂ x b(x, y) can be inverted in the sense that q = ∂ x b(x, y) ⇐⇒ y = y b (x, q) for some map y b , then one can actually deduce y(x) from the knowledge of q(x) := ∇u(x) by the relation y(x) = y b (x, q(x)). Replacing y(x) by y b (x, q(x)) in (2.3) and defining Γ b (x ′ , x, q) := b(x ′ , y b (x, q))b(x, y b (x, q)), we obtain the following reformulation of the incentive-compatibility constraint:

u(x ′ ) -u(x) ≥ b(x ′ , y(x)) -b(x, y(x)), ∀x, x ′ . ( 2
u(x ′ ) -u(x) ≥ Γ b (x ′ , x, q(x)), ∀x, x ′ . (2.4)
Note that (2.4) implies that, if u is differentiable at x, then ∇u(x) = ∂ x b(x, y b (x, q(x)), i.e. q(x) = ∇u(x) so that constraint (2.4) encodes the relation q = ∇u. Using the variables (q, u) instead of (y, p) (so that p(x) = b(x, y(x))-u(x) and q(x) = ∂ x b(x, y(x)), y(x) = y b (x, q(x))), the monopolist's program rewrites inf (u,q) X L(x, u(x), q(x))dx subject to u ≥ u 0 and (2.4), with L(x, u, q) := φ(x, b(x, y b (x, q))u, y b (x, q)). Of course, a crucial role is played here by the map y b but also by the function Γ b (x ′ , x, q) := b(x ′ , y b (x, q))b(x, y b (x, q)) which appears in the right-hand side of the incentivecompatibility condition written in the form (2.4).

The case of the convexity constraint

Variational problems subject to a convexity constraint arise in different applied settings: economics with Rochet-Choné model [START_REF] Rochet | Ironing, sweeping, and multidimensional screening[END_REF] (which corresponds to the principal-agent problem with a bilinear b), but also Newton's least resistance problem (see [START_REF] Brock | A symmetry problem in the calculus of variations[END_REF], [START_REF] Buttazzo | Minimum problems over sets of concave functions and related questions[END_REF], [START_REF] Lachand-Robert | An example of non-convex minimization and an application to Newton's problem of the body of least resistance[END_REF]). We shall also see in section 5 how to relate such variational problems with the computations of convex envelopes, a problem with its own interest. Such problems consist in minimizing, given a convex domain X of R d , an integral functional of the form X L(x, u(x), ∇u(x))dx among convex functions. This is in fact equivalent to minimize X L(x, u(x), q(x))dx subject to the constraint

u(x ′ ) -u(x) ≥ (x ′ -x) • q(x), ∀x ′ , x
which automatically implies convexity of u as well as q ∈ ∂u hence q = ∇u a.e.. In other words, the constraints that u is convex and q = ∇u are a particular case of the non-local condition (2.4) with b(x, y) = x • y so that y b (x, q) = q and Γ b (x ′ , x, q) = (x ′x) • q. The idea of writing the convexity constraints in this way for numerical purposes was first used by Ekeland and Moreno-Bromberg [START_REF] Ekeland | An algorithm for computing solutions of variational problems with global convexity constraints[END_REF].

Existence and uniqueness

Let us now consider the following assumptions:

(A1) b ∈ C 1 (X × Y )
where X and Y are open subsets of R d with X convex and bounded, and b is uniformly semiconvex with respect to x, i.e. there exists λ ≥ 0 such that for every y ∈ Y , x → b(x, y) + λ 2 |x| 2 is convex on X, (A2) b is twisted (or satisfies the generalized Spence-Mirrlees condition), i.e. for every

x ∈ X, y → ∂ x b(x, y) is a diffeomorphism from Y to Q x := ∂ x b(x, Y
), the set Q x is closed and convex, and its graph Graph(Q (A3) For all (x ′ , x) ∈ X 2 the map

) := {(x, q) ∈ R d × R d , x ∈ X, q ∈ Q x } is
q ∈ Q x → Γ b (x ′ , x, q) := b(x ′ , y b (x, q)) -b(x, y b (x, q)), (3.1)
is convex and continuous on Q x .

Assumption (A1) ensures that every (finite) b-convex function u is semi-convex and therefore:

• u is subdifferentiable at every point, i.e. for every x ∈ X,

∂u(x) := {q ∈ R d : u(x + h) ≥ u(x) + q • h + o(h)} = ∅.
• u is differentiable at every point except possibly a set of Hausdorff dimension at most d -1 (hence of zero measure).

The general form of the variational problems we are interested in is:

inf (u,q) J(u, q) := X L(x, u(x), q(x))dx (3.2) subject to q(x) ∈ Q x , for a.e. x ∈ X, (3.3)
and u(x ′ )u(x) ≥ Γ b (x ′ , x, q(x)), for a.e. x ∈ X and all x ′ ∈ X,

where Γ b is defined by (3.1). Note that inequality (3.4) is an equality when x ′ = x so, under assumption (A2), whenever u is differentiable at x, we have q(x) = ∇u(x). Thanks to assumption (A3), first outlined by Figalli, Kim and McCann [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF], the constraint

(3.4) is convex in (u, q) (equivalently in (u, ∇u)); if, in addition L(x, •, •) is convex, then the minimization problem (3.2)-(3.4
) is a convex program. Therefore we have the following well-posedness result (existence essentially follows from [START_REF] Carlier | A general existence result for the principal-agent problem with adverse selection[END_REF] and uniqueness from [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF]; for the sake of completeness we give a short proof).

Theorem 3.1. Assume that

• b satisfies (A1), (A2) and (A3),
• L is a lsc integrand such that, for a.e. x, (u, q) ∈ R × Q x → L(x, u, q) is strictly convex, and it satisfies for some C > 0 the coercivity condition

L(x, u, q) ≥ C(|u| + |q| -1), for a.e. x ∈ X and all (u, q) ∈ R × Q x , (3.5)
• there exists (u 0 , q 0 ) satisfying (3.4)-(3.3) such that J(u 0 , q 0 ) < +∞.

Then the problem (3.2)-(3.4) admits a unique solution.

Proof. Let (u n , q n ) be a minimizing sequence for (3.2)- (3.4). By (3.5), (u n ) n is bounded in W 1,1 (X). Since the functions (u n ) are all λ-convex (with the same λ), taking a subsequence if necessary, we may assume that u n converges locally uniformly to some λ-convex u and q n = ∇u n converges to q := ∇u a.e. (see for instance [START_REF] Carlier | Calculus of variations with convexity constraint[END_REF]). Then by Fatou's Lemma, we deduce that J(u, q) is the infimum of the problem (3.2)- (3.4). Note that for a.e. x, q(x) ∈ Q x since this set is closed. It remains to show that (u, q) satisfies the constraint (3.4), but this is obvious by passing to the limit in the inequality

u n (x ′ ) -u n (x) ≥ Γ b (x ′ , x, q n (x))
which holds for all x ′ , all n and all x in a set of full measure on which we may further assume that q n (x) converges to q(x). This shows existence; uniqueness directly follows from the convexity of the constraints and the strict convexity of J.

To show consistency of discrete approximations in section 6, we shall need the following elementary result:

q ∈ Q x , one has y b (x, q) = y b (x, ∂ x b(x, y b (x, q))), (3.6)
and

Γ b (x ′ , x, q) -Γ b (x, x, q) = Γ b (x ′ , x, ∂ x b(x, y b (x, q))). (3.7) Proof. Let q := ∂ x b(x, y b (x, q)). By definition, ∂ x b(x, y b (x, q)) = q, which implies y b (x, q) = y b (x, q) (i.e. proves (3.6)) since b is twisted. Then (3.7) immediately fol- lows since Γ b (x ′ , x, q) -Γ b (x, x, q) = b(x ′ , y b (x, q)) -b(x, y b (x, q)) = b(x ′ , y b (x, q)) - b(x, y b (x, q)) = Γ b (x ′ , x, q).

A tractable specification for b

Checking assumption (A3) directly on b is not easy in practice since it involves y b via Γ b .

There are however (a few) known examples, which can be found in [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF], in connection with the seminal work of Ma, Trudinger and Wang [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF] on the regularity of optimal transport maps. For the sake of tractability, we shall restrict ourselves to a class of particular b for which computations can be performed explicitly up to a certain point. Namely, from now on, we will only consider perturbations of the scalar product. More precisely, we consider the following specification: Note that (4.2) is trivially satisfied in the following two cases:

(B1) b is of the form b(x, y) := x • y + f (x)g(y) (4.1) with (f, g) ∈ C 1 (X) × C 1 (R d )
• if g is Lipschitz on R d and ∇f L ∞ (X) ∇g L ∞ (R d ) < 1; • if ∇f (x) • ∇g(y) = 0 for all (x, y) ∈ X × R d (which is in particular the case if for instance f (x) = f (x 1 ) and g(y) = g(y 2 , . . . , y d )).
It is known that suitable perturbations of the scalar product of the form (4.1) satisfy the convexity assumption (A3) (see [START_REF] Figalli | When is multidimensional screening a convex program?[END_REF], [START_REF] Loeper | On the regularity of solutions of optimal transportation problems[END_REF], [START_REF] Ma | Regularity of potential functions of the optimal transportation problem[END_REF]). We shall give in the next paragraph an elementary and self-contained proof as well as some properties of y b and Γ b under assumption (B1) (which is slightly weaker than the one considered in the previous references).

Properties of y b and Γ b

Note first that (A1) is satisfied with λ = 0 since f is convex and g ≥ 0, therefore bconvex functions are convex as suprema of convex functions. Next, given (x, q) ∈ X ×R d , y = y b (x, q) is obtained by solving

y = q -g(y)∇f (x),
or equivalently by solving the scalar equation

λ -g(q -λ∇f (x)) = 0 (4.3)
and setting y = q -λ∇f (x); (4.3) has a unique nonnegative root since its left-hand side, as a function of λ, has derivative not less than 1 + κ > 0 by (4.2) and value -g(q) ≤ 0 for λ = 0. This shows that (A2) holds with

Q x = R d . Moreover, y b (x, q) = q -g b (x, q)∇f (x) (4.4)
where

g b (x, q) := g(y b (x, q)) (4.5)
is the unique root of (4.3). Thanks to (4.2), one readily checks that

0 ≤ g b (x, q) ≤ g(q) 1 + κ . (4.6)
By (4.2), (4.6) and the fact that g is locally Lipschitz, one also deduces that for any compact subset K of R d there exists a positive constant C K such that, for all x, x ′ ∈ X and q, q ′ ∈ K,

|g b (x, q) -g b (x ′ , q ′ )| ≤ C K (|q -q ′ | + |∇f (x) -∇f (x ′ )|). (4.7)
which implies the continuity of g b , as well as that of y b and Γ b by (4.4). Let us now check that the convexity assumption (A3) is satisfied. By direct computation,

Γ b (x ′ , x, q) = (x ′ -x) • q + D f (x ′ , x)g b (x, q)
where D f is the Bregman divergence associated to f :

D f (x ′ , x) := f (x ′ ) -f (x) -∇f (x) • (x ′ -x),
which is nonnegative since f is convex. The convexity of Γ b with respect to q thus amounts to showing that g b (x, q) is convex with respect to q, which can be done as follows. Let (x, q 0 , q 1 , t)

∈ X × R d × R d × [0, 1], and set q t := (1 -t)q 0 + tq 1 , λ 0 := g b (x, q 0 ), λ 1 := g b (x ′ , q ′ ), and λ t = (1 -t)λ 0 + tλ 1 . Then g b (x, q t ) is the root of λ → λ -g(q t -λ∇f (x))
, which is increasing with respect to λ by (4.2) and nonnegative at λ = λ t by convexity of g. We deduce that g b (x, q t ) ≤ λ t , which shows that (A3) is satisfied.

Recall that if u is b-convex, then it is convex hence differentiable a.e. on X. Moreover if we denote by A ⊂ X the set of points at which u is differentiable, then ∇u is continuous on A (see [START_REF] Rockafellar | Convex analysis[END_REF]). The constraints that u is b-convex and q = ∇u can be expressed as

u(x ′ ) -u(x) ≥ Γ b (x ′ , x, q(x)), ∀(x ′ , x) ∈ X × A. Now if x ∈ X \ A,
one can choose a sequence x n in A, remaining away from ∂X and converging to x. Since u is locally Lipschitz, we may further assume that the bounded sequence q n := ∇u(x n ) converges to some q. Up to redefining q at x by setting q(x) = q, and using the fact that Γ b is continuous, we see that there is no loss of generality in imposing that the inequality u(x ′ )u(x) ≥ Γ b (x ′ , x, q(x)) actually holds for every

(x ′ , x) ∈ X 2 .

Further specification for the numerical examples

We will consider in our simulations (section 8) a special case where g b (and thus y b and Γ b ) can be explictly computed. This special case is given by g(y) := 1 + |y| 2 and f any C 1 convex function with Lipschitz constant on X strictly less than 1 (for instance f (x) = 1 + |x| 2 in our simulations). In this case, the solution g b (x, q) to (4.3) is the unique positive root of the quadratic equation

λ 2 = 1 + |q -λ∇f (x)| 2 ,
and thus has the following explicit expression

g b (x, q) = [(q, ∇f (x)) 2 + (1 -|∇f (x)| 2 )(1 + |q| 2 )] 1/2 -q • ∇f (x) 1 -|∇f (x)| 2 .
Note that explicit expressions are also obtained for any f and g satisfying (4.2) with ∇f (x) • ∇g(y) = 0 for all (x, y) ∈ X × R d . Indeed, in this case, λ → g(q -λ∇f (x)) is constant hence g b (x, q) = g(q). This is in particular the case when f (x) = f (x 1 ) and g(y) = g(y 2 , . . . , y d ).

A tractable specification for L

We have addressed in section 4 the choice of a tractable class of b so that conditions (A1), (A2) and (A3) are satisfied. To make the problem (3.2) tractable numerically, we shall restrict ourselves to projection problems, i.e. to quadratic Lagrangians L. From now on, we assume that b is a perturbation of the scalar product of the form (4.1) which satisfies assumption (B1), and we consider the quadratic problem

inf (u,q) J(u, q) := X α(x) 2 |q(x) -q 0 (x)| 2 + β(x) 2 |u(x) -u 0 (x)| 2 dx (5.1) subject to u(x ′ ) -u(x) ≥ Γ b (x ′ , x, q(x)), ∀(x ′ , x) ∈ X 2 , (5.2)
where

(u, q) ∈ C(X) ∩ H 1 (X) × L 2 (X, R d ) and (u 0 , q 0 , α, β) ∈ C(X, R) × C(X, R d ) × C(X, (0, +∞)) 2 .
We denote

K b := {(u, q) ∈ C(X) ∩ H 1 (X) × L 2 (X, R d ) : (u, q) satisfies (5.2)}. ( 5.3)
We have already seen that if (u, q) ∈ K b , then necessarily u is b-convex and q = ∇u, and explained at the end of section 4.1 why the b-convexity constraint can be imposed on the whole of X × X. We know from theorem 3.1 that (5.1)-(5.2) admits a unique solution which we denote (u, q) = (u, ∇u). Of course, there is no particular extra difficulty in considering additional pointwise convex constraints on u or q such as u

≥ v, u ≤ v, u = v on ∂X or q ∈ Q with Q a closed convex subset of R d .
The reason why we consider only quadratic Lagrangians L, i.e. projection problems, is that we wish to keep numerical computations as simple as possible. It is in fact possible (but computationally costly) to use Dykstra's algorithm with general Bregman distances to address convex but non quadratic Lagrangians, but we leave this for future research. We may also view the projection problem as a first step to address more general Lagrangians for instance by a projected gradient method, each step of which amounts to solving a problem of the form (5.1)-(5.2). We present in the next paragraphs some examples which, we believe, motivate directly the choice of a quadratic L. Let us also remark that we require α and β to be everywhere strictly positive, that is we ask L to be strongly convex in both variables u and q; this is due to the fact that in Dykstra's iterative projection algorithm, the condition q = ∇u is not enforced during the iterations, so to guarantee convergence, this nondegeneracy condition is necessary.

Rochet and Choné principal-agent model

In [START_REF] Rochet | Ironing, sweeping, and multidimensional screening[END_REF], Rochet and Choné considered the principal-agent problem in the case of a bilinear utility, a quadratic production cost c(y) := 1 2 |y| 2 and 0 reservation utility. The Rochet and Choné problem then consists in the quadratic minimization problem with a convexity and a nonnegativity constraints:

inf X [ 1 2 |∇u(x)| 2 -x • ∇u(x) + u(x)]ρ(x)dx : u convex, u ≥ 0
which can equivalently be rewritten as

inf X [ 1 2 |q(x)| 2 -x • q(x) + u(x)]ρ(x)dx : (u, q) ∈ K b , u ≥ 0 (5.4)
with b(x, y) = x • y and K b defined by (5.3). This is not exactly a projection problem of the form (5.1)-(5.2) because there is no explicit term in u 2 in (5.4), so in our simulations we will regularize it adding an extra quadratic term with a (small) parameter ε > 0.

Namely, the regularized version of (5.4) then reads

inf X [ 1 2 |q(x)| 2 -x • q(x) + ε 2 u 2 (x) + u(x)]ρ(x)dx : (u, q) ∈ K b , u ≥ 0 , (5.5)
which corresponds to a projection problem with q 0 = x, u 0 = -ε -1 , α = ρ, β = ερ, and an additional nonnegativity constraint on u.

Convex envelopes

Let X be a bounded open convex subset of R d and u 0 ∈ C(X) ∩ H 1 (X), the convex envelope, u * * 0 of u 0 is by definition the largest convex function which is pointwise below u 0 in X. Obviously, u * * 0 minimizes X (uu 0 ) 2 among all convex functions u such that u ≤ u 0 in X, but this is not a projection problem of the form (5.1)-(5.2) because here we have no term in |∇u| 2 . One could of course regularize this minimization problem by adding ε 2 X |∇u| 2 (we will actually do it to approximate more general envelopes in the next paragraph). It turns out however that the convex envelope also minimizes a Dirichlet integral. Indeed, let u be a H 1 convex function such that u ≤ u 0 and u = u * * 0 on ∂X; then by convexity of the square and an integration by parts,

1 2 X (|∇u| 2 -|∇u * * 0 | 2 ) ≥ X ∆u * * 0 (u * * 0 -u),
and the latter is nonnegative since u ≤ u * * 0 , and u * * 0 being convex, it has a nonnegative Laplacian. Hence u * * 0 is also the minimizer of the Dirichlet integral 1 2 X |∇u| 2 among all convex functions which agree with u * * 0 on ∂X and are below u 0 on X. Therefore provided that u * * 0 = u 0 on ∂X, the convex envelope u * * 0 can be obtained by solving

inf X α 2 |∇u| 2 + β 2 |u -u 0 | 2 : u convex, u ≤ u 0 , (u -u 0 )| ∂X = 0
for every positive constants α and β. Of course, to make this observation relevant in practice to compute u * * 0 , we should already know u * * 0 on ∂X. It is however well-known that u 0 = u * * 0 on ∂X (equivalently, u 0 agrees with a convex function on ∂X) in the following cases:

• if ∂X contains no segment, by the following formula for the convex envelope:

u * * 0 (x) = inf d+1 i=1 λ i u 0 (x i ) : λ i ≥ 0, x i ∈ X, d+1 i=1 λ i = 1, d+1 i=1 λ i x i = x ,
which implies that u * * 0 = u 0 on ∂X if every boundary point of X is extreme;

• if the dimension d = 1, as a particular case of the previous one;

• if u 0 is the sum of a convex function and a function with compact support in X.

It is worth mentioning here some purely PDE approaches to the computation of convex envelopes. Vese [START_REF] Vese | A method to convexify functions via curve evolution[END_REF] showed that the convex envelope can be obtained as the limit for large time of a solution of a nonlinear parabolic equation (also see [START_REF] Carlier | Exponential convergence for a convexifying equation[END_REF] for exponential convergence results) and Oberman [START_REF] Oberman | The convex envelope is the solution of a nonlinear obstacle problem[END_REF], [START_REF] Oberman | Computing the convex envelope using a nonlinear partial differential equation[END_REF] developed a direct approach based on the observation that convex envelopes solve a nonlinear elliptic equation.

More general envelopes

Since by definition, b-convex functions are stable by suprema, one can also define the b-convex envelope of a function u 0 as the largest b-convex function below u 0 . If b is as in section 4 (or more generally if b(., y) is convex for every y), then b-convex functions are convex hence subharmonic; the same argument as for convex envelopes thus enables to conclude that if u 0 coincides with its b-convex envelope on ∂X, then the b-convex envelope of u 0 can also be obtained by the quadratic minimization problem

inf 1 2 X α 2 |q| 2 + β 2 |u -u 0 | 2 : (u, q) ∈ K b , u ≤ u 0 , (u -u 0 )| ∂X = 0 . (5.6)
Let us mention two cases where the convex envelope of u 0 coincides with u 0 on ∂X:

• if u 0 is the sum of a b-convex function and a nonnegative function which vanishes on ∂X;

• if the dimension d = 1 and u 0 has compact support in X (see below).

Lemma 5.1. Assume that b satisfies (B1) with X = (0, 1) and let u 0 ∈ C c ((0, 1)).

Then u 0 coincides with its b-convex envelope at 0 and 1.

Proof. Let ε > 0 be such that

u 0 = 0 on [0, ε] ∪ [1 -ε, 1].
Let us denote by v the b-convex envelope of u 0 ; by definition v ≤ u 0 up to the boundary of X. To prove equality at 0 and 1, it is enough to construct a b-convex function w such that w ≤ u 0 and w(0) = w(1) = 0. Define w by w(x) := max{w 1 (x), w 2 (x)} where

w 1 (x) := b(x, y 1 ) -b(0, y 1 ), w 2 (x) := b(x, y 2 ) -b(1, y 2 ),
with y 1 < 0 < y 2 and |y i | ≥ M, i = 1, 2, with M to be chosen properly. By (4.2), we have w

′ 1 (x) ≤ y 1 (1 + κ) + f ′ ∞ |g(0)| ≤ -M(1 + κ) + f ′ ∞ |g(0)| where κ > -1. Then choosing M large enough, w 1 (x) ≤ -M 2 (1 + κ)x for every x ∈ [0, 1].
If we further restrict M so that -M 2 (1 + κ) < ε -1 min u 0 , we get w 1 ≤ u 0 . In a similar way for such an M, we also have w 2 ≤ u 0 . Hence v ≥ w, so that v(0) ≥ w(0) = 0 and v(1) ≥ w(1) = 0.

In the general case where one cannot take for granted that u 0 coincides with its b-convex envelope on ∂X, one can still approximate the b-convex envelope by solving for a small ε > 0 the problem

inf 1 2 X ε 2 |q| 2 + 1 2 |u -u 0 | 2 : (u, q) ∈ K b , u ≤ u 0 .
(5.7)

Discretization and convergence

We now address the discretization of (5.1)-(5.2) and our aim in this section is to prove a Γ-convergence result. To prove Γ-convergence, it will be useful to have that Lipschitz functions are dense in energy for the problem (5.1)-(5.2), which requires the extra assumption that g is Lipschitz on R d (this would not be needed for variants of (5.1)-(5.2) with additional constraints that bound uniformly q). Note first that since here b-convex functions are convex, they are locally Lipschitz hence

K b ⊂ W 1,∞ loc (X) × L ∞ loc (X). Let (u, q) ∈ K b , M > 0 and (u M , q M ) be defined by u M (x ′ ) := sup x∈X, |q(x)|≤M {u(x) + Γ b (x ′ , x, q(x))}, q M := ∇u M . (6.1) By construction u M is b-convex, u M ≤ u, u M = u on A M := {x ∈ X : |q(x)| ≤ M};
since q is locally bounded, u M converges locally uniformly to u and q M converges to q almost everywhere. If x / ∈ A M , let x n ∈ A M be such that

u M (x) = lim n u(x n ) + Γ b (x, x n , q(x n ));
since |q(x n )| ≤ M, passing to a subsequence, we may assume that (x n , q n , y n ) = (x n , q(x n ), y b (x n , q(x n ))) converges to some (x, q, y) with x ∈ X, |q| ≤ M and y = y b (x, q) since y b is continuous by (4.7). Let now x ′ ∈ X; since u M (x ′ ) ≥ u(x n ) + Γ b (x ′ , x n , q n ), using the continuity of Γ b , we get u M (x ′ )u M (x) ≥ Γ b (x ′ , x, q) -Γ b (x, x, q), so that if u M is differentiable at x (which is the case almost everywhere), we necessarily have

∇u M (x) = q M (x) = ∂ x b(x, y b (x, q)). Lemma 3.2 thus gives u M (x ′ ) -u M (x) ≥ Γ b (x ′ , x, q) -Γ b (x, x, q) = Γ b (x ′ , x, q M (x)),
so that (u M , q M ) ∈ K b . Moreover, recalling (4.4) and (4.5), we have

q M (x) = ∂ x b(x, y b (x, q)) = y b (x, q) + g b (x, q))∇f (x) = q + g b (x, q)(∇f (x) -∇f (x)),
so that by (4.6),

|q M (x)| ≤ |q| + 2 ∇f L ∞ (X) 1 + κ g(q). (6.2)
Lemma 6.1. Assume that b satisfies (B1) with the extra assumption that g is Lipschitz on R d . Let (u, q) ∈ K b and (u M , q M ) be defined by (6.1). Then there exists a constant C (independent of M and (u, q)) such that |q M | ≤ C(1 + |q|), which implies that u M converges to u in H 1 as M → +∞ and that the minimum of J over K b coincides with its infimum over

K b ∩ (W 1,∞ (X) × L ∞ (X)).
Proof. Only the inequality |q M | ≤ C(1 + |q|) has to be shown: the H 1 convergence will directly follow by Lebesgue's dominated convergence Theorem. If |q(x)| ≤ M, there is nothing to prove since q M = ∇u M = ∇u = q a.e. in A M . If |q(x)| > M and u M is differentiable at x, then by (6.2), the fact that |q| ≤ M and the assumption that g is Lipschitz, we have

|q M (x)| ≤ C(1 + M) ≤ C(1 + |q(x)|)
for some constant C.

Discretization

We now discretize (5.1)-(5.2) as follows. Denoting by λ the Lebesgue measure on X, we consider for each discretization parameter h > 0, a finite subset X h = {x h i , i = 1, . . . , N h } of X in such a way that the sequence of discrete measures

λ h := λ(X) N h N h i=1 δ x h i approaches λ, i.e. λ h * ⇀ λ as h → 0. (6.3)
We now consider the following discretized version of (5.1)-(5.2):

inf (u h ,q h )∈K h b J h (u h , q h ) := λ(X) N h N h i=1 α(x h i ) 2 |q h i -q 0 (x h i )| 2 + β(x h i ) 2 |u h i -u 0 (x h i )| 2 (6.4)
where

K h b := {(u h i , q h i ) ∈ (R × R d ) N h : u h i -u h j ≥ Γ b (x h i , x h j , q h j ), ∀(i, j) ∈ {1, . . . , N h } 2 }. (6.5)
Note that (6.4)-(6.5) is a finite-dimensional projection problem, and under our assumptions K h b is a closed and convex set so that existence and uniqueness of a minimizer which we denote by (u h , q h ) is straightforward.

Elements of K h b can easily be extended to elements of K b admissible for the continuous problem as follows: Lemma 6.2. Given (u h , q h ) ∈ K h b , let us extend u h to X by setting

u h (x) := max i=1,••• ,N h {u h i + Γ b (x, x h i , q h i )}, ∀x ∈ X.
For each x ∈ X, select an index i among those for which the maximum above is achieved which minimizes the distance between x and x h i and set

q h (x) := ∂ x b(x, y b (x h i , q h i )), so that q h (x h i ) = ∂ x b(x h i , y b (x h i , q h i )) = q h i .
Then the obtained extension (which by abusing notations we continue to denote by (u h , q h )) belongs to K b .

Proof. First observe that (u h , q h ) ∈ W 1,∞ (X) × L ∞ (X). Let x ∈ X and let i be an index for which

u h (x) = u h i + Γ b (x, x h i , q h i ) and q h (x) = ∂ x b(x, y b (x h i , q h i )). For x ′ ∈ X, we have u h (x ′ ) ≥ u h i + Γ b (x ′ , x h i , q h i ) hence u h (x ′ ) -u h (x) ≥ Γ b (x ′ , x h i , q h i ) -Γ b (x, x h i , q h i ) = Γ b (x ′ , x, q h (x))
where the last equality follows from Lemma 3.2.

Extending (u h , q h ) ∈ K h b as in Lemma 6.2, one can write

J h (u h , q h ) = X L(x, u h (x), q h (x))dλ h (x)
whereas the functional J for the continuous problem is, for all (u, q) ∈ K b ,

J(u, q) = X L(x, u(x), q(x))dx
where

L(x, u, q) := α(x) 2 |q -q 0 (x)| 2 + β(x) 2 |u -u 0 (x)| 2 , ∀(x, u, q) ∈ X × R × R d .

Convergence

This section is devoted to a detailed Γ-convergence proof. For the Γ-liminf inequality, the following result will be useful: Lemma 6.3. Let (u h , q h ) ∈ K h b extended as in Lemma 6.2 be such that

sup h q h L 2 (λ h ) < +∞. (6.6)
Assume that u h converges uniformly on compact subsets of X to some function u and that q h converges a.e. to q := ∇u. Then

1. (u, q) ∈ K b , 2. for every F ∈ C b (X × R d ), one has lim h→0 λ(X) N h N h i=1 F (x h i , q h i ) = lim h→0 X F (x, q h (x))dλ h (x) = X F (x, ∇u(x))dx, (6.7)
3. the following Γ-liminf inequality holds:

lim inf h J h (u h , q h ) ≥ J(u, q). (6.8) Proof. Let us define the discrete measure γ h on X × X × R d by γ h := λ h ⊗ ((id, q h ) # λ h ), i.e. for every F ∈ C b (X × X × R d ), X×X×R d F (x ′ , x, q)dγ h (x ′ , x, q) := X×X F (x ′ , x, q h (x))dλ h (x ′ )dλ h (x).
The family of measures γ h is tight because its first and second marginals are λ h which satisfies (6.3) and its third marginals have bounded second moments by (6.6). Up to a subsequence, we may therefore assume that γ h narrowly converges to some measure γ. Note that γ is necessarily of the form λ ⊗ θ where θ has first marginal λ hence can be disintegrated as θ(dx, dq) = λ(dx)θ x (dq), so that

lim h→0 X×X F (x ′ , x, q h (x))dλ h (x ′ )dλ h (x) = X×X× R d F (x ′ , x, q)dθ x (q) dx ′ dx for every F ∈ C b (X ×X ×R d ). Since (u h , q h ) ∈ K h b , we have u h (x ′ ) -u h (x) ≥ Γ b (x ′ , x, q) γ h -a.e.. Let then ψ ∈ C c (X × X × R d ) with ψ ≥ 0; since (x ′ , x, q) → ψ(x ′ , x, q)(u h (x ′ ) - u h (x) -Γ b (x ′ , x, q)) converges uniformly to ψ(x ′ , x, q)(u(x ′ ) -u(x) -Γ b (x ′ , x, q)
) and γ h converges narrowly to γ, we have

X×X×R d ψ(x ′ , x, q)(u(x ′ ) -u(x) -Γ b (x ′ , x, q))dγ(x ′ , x, q) ≥ 0 so that the continuous function u(x ′ )-u(x)-Γ b (x ′ , x, q
) is nonnegative for every x ′ ∈ X, a.e. x ∈ X and θ x -a.e. q. If x is a point of differentiability of u, this implies in particular that q = ∇u(x), so θ x necessarily coincides with the Dirac mass at ∇u(x). We then have γ = λ ⊗ ((id, ∇u) # λ) and the whole family γ h converges narrowly to γ as h → 0, hence

lim h→0 λ(X) 2 N 2 h N h j=1 N h i=1 F (x h j , x h i , q h i ) = X×X F (x ′ , x, ∇u(x))dx ′ dx for every F ∈ C b (X × X × R d )
, which in particular proves that (6.7) holds. Fixing M > 0, we have by (6.7) and (6.6)

X α(x) 2 min(M, |∇u(x) -q 0 (x)| 2 )dx ≤ lim inf h X α(x) 2 |q h (x) -q 0 (x)| 2 dλ h (x) ≤ C
(6.9) for some C independent of h. Letting M go to ∞, this shows that q = ∇u ∈ L 2 and thus (u, q) ∈ K b (the b-convexity inequality has already been established). In a similar way, taking a family of cutoff functions

ψ k ∈ C c (X), 0 ≤ ψ k ≤ 1 with ψ k = 1 on {x ∈ X : dist(x, ∂X) ≥ 1/
k} and ψ k ≤ ψ k+1 , using the uniform convergence of u h to u on compact subsets of X, we get

X ψ k (x) β(x) 2 (u(x) -u 0 (x)) 2 dx ≤ lim inf h X β(x) 2 (u h (x) -u 0 (x)) 2 dλ h (x). (6.10)
Passing to the supremum in M in (6.9) and in k in (6.10) finally establishes the Γ-liminf inequality (6.8):

lim inf h J h (u h , q h ) ≥ lim inf h X α(x) 2 |q h (x) -q 0 (x)| 2 dλ h (x) + lim inf h X β(x) 2 (u h (x) -u 0 (x)) 2 dλ h (x) ≥ J(u, q).
Remark 6.4. If in Lemma 6.3 one further assume that u h is uniformly Lipschitz, i.e.

sup h q h L ∞ < +∞,
then the measures γ h in the proof are supported by a fixed compact set and thus one can use a quadratic-test function in (6.7) which actually gives

X α(x) 2 |∇u(x) -q 0 (x)| 2 dx = lim h X α(x) 2 |q h (x) -q 0 (x)| 2 dλ h (x).
In a similar way, in this case u h converges uniformly to u in X and we also have

X β(x) 2 (u(x) -u 0 (x)) 2 dx = lim h X β(x) 2 (u h (x) -u 0 (x)) 2 dλ h (x).
Hence with a uniform bound on q h , the last statement in Lemma 6.3 can be streghtened to: lim h J h (u h , q h ) = J(u, q). Theorem 6.5. Assume that b satisfies (B1) with the extra assumption that g is Lipschitz on R d . Let (u h , q h ) be the solution of (6.4)-(6.5), extended as in Lemma 6.2, and let (u, q) be the solution of (5.1)-(5.2). Then u h converges locally uniformly to u and q h converges in L 2 loc and a.e. to q as h → 0.

Proof. Since (u h , q h ) solves (6.4)-(6.5) it is easy to see that there is a constant C such that, for all h,

q h L 2 (λ h ) + u h L 2 (λ h ) ≤ C. (6.11)
We first claim that for every ω convex such that ω ⊂⊂ X,

sup h ∇u h L ∞ (ω) = C ω < +∞.
Indeed, if it was not the case, using the fact that the all functions u h are convex and arguing as in the proof of Theorem 1 in [START_REF] Carlier | A numerical approach to variational problems subject to convexity constraint[END_REF], this would imply that, up to a subsequence, the norm of all elements of ∂u h converge to +∞ uniformly on some open subset of X, which would also imply a uniform explosion for q h since q h (x) ∈ ∂u h (x), and thus this would contradict (6.11) and (6.3). Passing to a subsequence if necessary we may therefore assume that for some convex function u, one has u h → u uniformly on compact subsets of X, q h → q = ∇u a.e. and in L 2 loc . By Lemma 6.3, (u, q) ∈ K b and J(u, q) ≤ lim inf h J h (u h , q h ). To prove that (u, q) = (u, q) (and then that the whole family converges, not only a subsequence), it is enough by Lemma 6.1 to show that J(u, q) ≤ J(v, p) for every

(v, p) ∈ K b ∩(W 1,∞ (X)×L ∞ (X)). Let (v, p) ∈ K b ∩ W 1,∞ (X) × L ∞ (X)
, then extend as in Lemma 6.2 the discrete data (v(x h i ), p(x h i )) i=1,...,N h and denote by (v h , p h ) this extension. It is easy to check that p h is uniformly bounded and v h converges uniformly to v, then using Remark 6.4 and the optimality of (u h , q h ) we get

J(u, q) ≤ lim inf h J h (u h , q h ) ≤ lim inf h J h (v h , p h ) = J(v, p),
which gives the desired result.

7 Numerical method

Dykstra's iterative projection algorithm

We focus on the discretized problem (6.4)- (6.5), where for notational simplicity we drop the discretization index h and set Γ ij (q) := Γ b (x i , x j , q):

inf (u,q) N k=1 α k 2 q k -q 0 k 2 + β k 2 u k -u 0 k 2 subject to u i -u j ≥ Γ ij (q j ), ∀(i, j){1, . . . , N} 2 .
Defining the convex subsets of R N × R dN C i,j := (u, q) : u iu j ≥ Γ ij (q j ) , ∀(i, j){1, . . . , N} 2 , and the weighted squared Euclidean distance in R N × R dN D α,β ((u, q), (u ′ , q ′ )) :=

N k=1 α k 2 |q k -q ′ k | 2 + β k 2 |u k -u ′ k | 2 ,
the discrete problem is to find the projection P α,β C (u 0 , q 0 ), solution in R N × R dN to inf (u,q) D α,β (u, q), (u 0 , q 0 ) subject to (u, q) ∈ C := (i,j)

C i,j .

This projection problem can be solved iteratively by Dykstra's algorithm [START_REF] Bauschke | Dykstra's algorithm with Bregman projections: a convergence proof[END_REF][START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF]. Let us now recall it for a projection problem onto a closed convex subset K of R m which can be written as an intersection of elementary closed convex subsets K 1 , . . . , K L of R m :

K := L l=1 K l .
Given z 0 , we wish to compute P α,β K (z 0 ) = argmin z∈K D α,β (z, z 0 ). First extend by Lperiodicity the familly of convex sets (K l ) l=1,••• ,L by setting

K l+kL := K l , ∀(k, l) ∈ N × {1, . . . , L}.
Then initialize the algorithm by setting

θ -L+1 = • • • = θ -1 = θ 0 = 0
and starting from the point z 0 we want to project, update z n and θ n for n ≥ 1 by:

z n = P α,β Kn z n-1 + θ n-L , θ n = z n-1 -z n + θ n-L .
The fact that the sequence z n converges to P α,β K (z 0 ) has been established by Boyle and Dyskstra [START_REF] Boyle | A method for finding projections onto the intersection of convex sets in Hilbert spaces[END_REF]; Bauschke and Lewis [START_REF] Bauschke | Dykstra's algorithm with Bregman projections: a convergence proof[END_REF] have extended this result to projection problems with general Bregman distances.

Remark 7.1. For the sake of simplicity, we have omitted above the case of pointwise convex constraints on u or q such as u i ≥ 0, u i ≤ v i , |q i | ≤ r i , or q i ∈ R d + . Such constraints can be handled easily with Dykstra's algorithm by adding elementary convex sets on which the projections are totally explicit.

Remark 7.2. The main drawback of Dystra's algorithm in our context is the a priori very large number N 2 of elementary convex sets on which we have to project. One can however significantly reduce the computational cost if C = ∩ (i,j)∈I C i,j where the cardinality of I ⊂ {1, . . . , N} 2 is smaller than N 2 . This reduction of the number of constraints can be done when the b-convexity constraint propagates. This is in particular the case in dimension 1 under the so-called Spence-Mirrlees condition ∂ 2 xy b > 0 and for ordered grid points (x i < x j for i < j), where it is enough to impose the bconvexity on neighbour points (j = i ± 1), resulting in a number of elementary convex sets O(N) instead of O(N 2 ). It is indeed easy to see that the constraint propagates in this case: set y i := y b (x i , q i ) and first remark that the two "local" constraints u i+1u i ≥ b(x i+1 , y i )b(x i , y i ) and u iu i+1 ≥ b(x i , y i+1 )b(x i+1 , y i+1 ) together with Spence-Mirrlees condition imply that y i+1 ≥ y i ; then summing the two consecutive constraints u i+1 -u i ≥ b(x i+1 , y i )-b(x i , y i ) and u i+2 -u i+1 ≥ b(x i+2 , y i+1 )-b(x i+1 , y i+1 ) and using the fact that b(x i+2 , y i+1 )b(x i+1 , y i+1 ) ≥ b(x i+2 , y i )b(x i+1 , y i ) again by Spence-Mirrlees condition and the fact that y i+1 ≥ y i , we obtain that u i+2u i ≥ b(x i+2 , y i )b(x i , y i ), which is exactly what we mean by propagation of the constraint.

Remark 7.3. In a similar way, for the convexity constraint b(x, y) = x • y, Γ b (x ′ , x, q) = (x ′x) • q and for regular grid points, the number of convexity constraints u iu j ≥ (x ix j ) • q j can be significantly reduced (for instance for aligned points, it is enough to consider consecutive points x i and x j ). We refer to [START_REF] Carlier | A numerical approach to variational problems subject to convexity constraint[END_REF] and the recent efficient approach of Mirebeau [START_REF] Mirebeau | Adaptive, anisotropic and hierarchical cones of discrete convex functions[END_REF] for details.

Elementary projections

Lagrangian relaxation

In Dyskstra's algorihm, we have to compute at each step an elementary projection onto a single convex C i,j , P α,β C i ,j (q, û) := P i,j (q, û) for various (q, û). The elementary projection (q, u) = P i,j (q, û) is obviously given by q k = qk for k = j and u k = ûk for k = i, j with (q j , u i , u j ) obtained by solving the low-dimensional projection problem inf

(u i ,u j ,q j ) α j 2 |q j -qj | 2 + β i 2 |u i -ûi | 2 + β j 2 |u j -ûj | 2 subject to u i -u j ≥ Γ ij (q j ).
If ûiûj ≥ Γ ij ( qj ), then this projection problem is trivial and (u i , u j , q j ) = (û i , ûj , qj ), so only the case ûiûj < Γ ij ( qj ) where the constraint is binding requires some specific attention and can conveniently be addressed by Lagrangian relaxation as follows. Since the problem above is convex and qualified (the constraint is linear in u i and u j ), there is no duality gap with its dual problem, which is

sup λ≥0 l ij (λ)
where l ij (λ) is defined, for λ ≥ 0, by inf

(u i ,u j ,q j ) α j 2 |q j -qj | 2 + β i 2 |u i -ûi | 2 + β j 2 |u j -ûj | 2 + λ (-u i + u j + Γ ij (q j )) . (7.1)
The dual function l ij is concave -as an infimum of affine functions -and differentiablesince the Lagrange problem (7.1) has a unique solution for any λ ≥ 0, that we denote by

(q j (λ), u i (λ), u j (λ)) -with l ′ ij (λ) = -u i (λ) + u j (λ) + Γ ij (q j (λ)). Note that u i (λ), u j (λ) are given explicitly by u i (λ) = ûi + λ β i , u j (λ) = ûj - λ β j , (7.2)
whereas q j (λ) is the unique solution (not explicit in general) to α j (q jqj ) + λ∇Γ ij (q j ) = 0, i.e. q j (λ) = id + λ α j ∇Γ ij -1 (q j ).

The solution of the projection problem (u i , u j , q j ) in the case where ûiûj < Γ ij ( qj ) is then (u i (λ), u j (λ), q j (λ)) where λ > 0 solves

l ′ ij (λ) = -u i (λ) + u j (λ) + Γ ij (q j (λ) = 0. (7.4)
Note that (7.4) is monotone with respect to λ since l ij is concave.

In practice, we have experimented two numerical methods to solve the elementary projection problem in the case λ > 0. The first one is a duality method: it consists in solving the dual problem via (7.4):

1. Given λ > 0, compute q j (λ) from (7.3) (by Newton's method if necessary) and u i (λ), u j (λ) from (7.2) (explicitly).

Evaluate l

′ ij (λ) = -u i (λ) + u j (λ) + Γ ij (q j (λ))
and stop if λ is satisfying (7.4) up to a fixed tolerance.

3. Update λ > 0 by one step of a dichotomy method to solve l ′ ij (λ) = 0. The second method that we have experimented (and actually used for the numerical results in section 8) is a primal-dual method: it consists in finding simultaneously the primal-dual solution (q j (λ), u i (λ), u j (λ), λ) to the elementary projection problem via the optimality conditions for λ > 0. This amounts to computing simultaneously (q j (λ), λ) as the unique solution to

   α j (q j -qj ) + λ∇Γ ij (q j ) = 0 -(û i + λ β i ) + (û j - λ β j ) + Γ ij (q j ) = 0, (7.5)
by Newton's method. See e.g. [START_REF] Bonnans | Numerical optimization[END_REF] for variants.

Special cases

In the convexity constraint case, Γ ij (q) = (x ix j ) • q, the solution q j (λ) to equation (7.3) is explicitly given by

q j (λ) = qj - λ α j (x i -x j ).
Plugging it, together with u i (λ), u j (λ) given by (7.2), into the optimality equation (7.4), we get that the solution to the dual problem when λ > 0 (i.e. when (q j , ûi , ûj ) is not feasible for the primal problem) is

λ = -û i + ûj + Γ ij (q j ) 1 β i + 1 β j + |x i -x j | 2 α j .
Therefore the elementary projections are explicit in this case.

In the b-convexity constraint case with b(x, y) := x • y + f (x)g(y), recall that

Γ ij (q) = (x i -x j ) • q + D f (x i , x j )g b (x j , q).
When g b is known explicitly (as it is in the special cases described in 4.2) and twice differentiable, the elementary projection problems can be solved as explained above.

In practice, we solve (7.5) by Newton's method initialized with the explicit solution (q j , λ) to the convexity constraint case. In the more general case where g b is not known explicitly, one could imagine to solve simultaneously (4.3) to evaluate it numerically, but we have not explored this way.

Numerical results

We present here various examples of projection problems of the form (5.1)-(5.2) (sometimes on convex subsets of C corresponding to additional pointwise convex constraints) and their numerical resolution by the method described in sections 6 and 7. Since the main limitation to this method (up to the numerical resolution of (4.3) to get g b ) is the fact that Γ b has to be explicit (very likely y b too), we have restricted ourselves to two cases for b:

• the bilinear case b(x, y) := x • y, which corresponds to the convexity constraint;

• a perturbation of it given by b(x, y) := x • y + f (x)g(y) with f (x) := 1 + |x| 2 , g(y) := 1 + |y| 2 .

In 3d figures, the z-axis is directed downward. The first example is precisely the projection problem (5.1)-(5.2) with α = β = 1 and q 0 = ∇u 0 . We have solved numerically such problems in dimension d = 2, see Figure 1 and 2 for the H 1 ((0, 1) 2 ) projection of two different u 0 .

b-convex envelopes in dimension one

We next have computed convex and b-convex envelopes in dimension d = 1 using the formulation (5.6) (both for convex and b-convex envelopes) with α = β = 1, q 0 = 0 and a nonconvex u 0 depicted in Figure 3.

Approximated b-convex envelopes in dimension two

We now present two-dimensional approximated convex/b-convex envelopes, using the regularized formulation (5.7) with different values of the regularization parameter ε > 0, respresented in Figure 4.

  closed; this enables one to define the map y b on Graph(Q) by y = y b (x, q) ⇐⇒ ∂ x b(x, y) = q.

  where X is an open bounded convex subset of R d (here Y = R d ), f and g are convex, g ≥ 0, and inf (x,y)∈X×R d ∇f (x) • ∇g(y) =: κ > -1. (4.2)

8. 1 H

 1 1 projection on b-convex functions b-convex functions Projection on b-convex functions

Figure 1 -

 1 Figure 1 -H 1 ((0, 1) 2 ) projections of -x(1x)(2x -1) 2 y(1y)(2y -1) 2

  b-convex functions Projection on b-convex functions

Figure 2 -

 2 Figure 2 -H 1 ((0, 1) 2 ) projections of -x(1x)(y -12 )

Figure 3 -

 3 Figure 3 -Envelopes of a function in dimension 1

Figure 4 -

 4 Figure 4 -Approximated envelopes of a function in dimension 2

  .

  Distribution of the types of products q = ∇u

Figure 5 -

 5 Figure 5 -Rochet and Choné problem regularized with ε = 0.01

(ANR-12-MONU-0013), OPTIFORM (ANR-12-BS01-0007) and from INRIA through the "action exploratoire" MOKAPLAN and wish to thank J.-D. Benamou for stimulating discussions. G.C. gratefully acknowledges the hospitality of the Mathematics and Statistics Department at UVIC (Victoria, Canada) and support from the CNRS.