Section: Research Program

Numerical Tools

The above continuous models require a careful discretization, so that the fundamental properties of the models are transferred to the discrete setting. Our team aims at developing innovative discretization schemes as well as associated fast numerical solvers, that can deal with the geometric complexity of the variational problems studied in the applications. This will ensure that the discrete solution is correct and converges to the solution of the continuous model within a guaranteed precision. We give below examples for which a careful mathematical analysis of the continuous to discrete model is essential, and where dedicated non-smooth optimization solvers are required.

Geometric Discretization Schemes

Discretizing the cone of convex constraints.

(Participants: J-D. Benamou, G. Carlier, J-M. Mirebeau, Q. Mérigot) Optimal transportation models as well as continuous models in economics can be formulated as infinite dimensional convex variational problems with the constraint that the solution belongs to the cone of convex functions. Discretizing this constraint is however a tricky problem, and usual finite element discretizations fail to converge.

Our expertise: Our team is currently investigating new discretizations, see in particular the recent proposal  [68] for the Monge-Ampère equation and  [154] for general non-linear variational problems. Both offer convergence guarantees and are amenable to fast numerical resolution techniques such as Newton solvers. Since  [68] explaining how to treat efficiently and in full generality Transport Boundary Conditions for Monge-Ampère, this is a promising fast and new approach to compute Optimal Transportation viscosity solutions. A monotone scheme is needed. One is based on Froese Oberman work  [128], a new different and more accurate approach has been proposed by Mirebeau, Benamou and Collino  [66]. As shown in  [113], discretizing the constraint for a continuous function to be convex is not trivial. Our group has largely contributed to solve this problem with G. Carlier   [104], Quentin Mérigot  [157] and J-M. Mirebeau   [154]. This problem is connected to the construction of monotone schemes for the Monge-Ampère equation.

Goals: The current available methods are 2-D. They need to be optimized and parallelized. A non-trivial extension to 3-D is necessary for many applications. The notion of $c$-convexity appears in optimal transport for generalized displacement costs. How to construct an adapted discretization with “good” numerical properties is however an open problem.

Numerical JKO gradient flows.

(Participants: J-D. Benamou, G. Carlier, J-M. Mirebeau, G. Peyré, Q. Mérigot) As detailed in Section 2.3, gradient Flows for the Wasserstein metric (aka JKO gradient flows  [137]) provides a variational formulation of many non-linear diffusion equations. They also open the way to novel discretization schemes. From a computational point, although the JKO scheme is constructive (it is based on the implicit Euler scheme), it has not been very much used in practice numerically because the Wasserstein term is difficult to handle (except in dimension one).

Our expertise:

Solving one step of a JKO gradient flow is similar to solving an Optimal transport problem. A geometrical a discretization of the Monge-Ampère operator approach has been proposed by Mérigot, Carlier, Oudet and Benamou in [65] see Figure 4. The Gamma convergence of the discretisation (in space) has been proved.

Goals: We are also investigating the application of other numerical approaches to Optimal Transport to JKO gradient flows either based on the CFD formulation or on the entropic regularization of the Monge-Kantorovich problem (see section 3.2.3). An in-depth study and comparison of all these methods will be necessary.

Sparse Discretization and Optimization

From discrete to continuous sparse regularization and transport.

(Participants: V. Duval, G. Peyré, G. Carlier, Jalal Fadili (ENSICaen), Jérôme Malick (CNRS, Univ. Grenoble)) While pervasive in the numerical analysis community, the problem of discretization and $\Gamma$-convergence from discrete to continuous is surprisingly over-looked in imaging sciences. To the best of our knowledge, our recent work  [120], [121] is the first to give a rigorous answer to the transition from discrete to continuous in the case of the spike deconvolution problem. Similar problems of $\Gamma$-convergence are progressively being investigated in the optimal transport community, see in particular  [105].

Our expertise: We have provided the first results on the discrete-to-continous convergence in both sparse regularization variational problems  [120], [121] and the static formulation of OT and Wasserstein barycenters  [105]

Goals: In a collaboration with Jérôme Malick (Inria Grenoble), our first goal is to generalized the result of  [120] to generic partly-smooth convex regularizers routinely used in imaging science and machine learning, a prototypal example being the nuclear norm (see  [179] for a review of this class of functionals). Our second goal is to extend the results of  [105] to the novel class of entropic discretization schemes we have proposed  [64], to lay out the theoretical foundation of these ground-breaking numerical schemes.

Polynomial optimization for grid-free regularization.

(Participants: G. Peyré, V. Duval, C. Poon) There has been a recent spark of attention of the imaging community on so-called “grid free” methods, where one tries to directly tackle the infinite dimensional recovery problem over the space of measures, see for instance  [96], [120]. The general idea is that if the range of the imaging operator is finite dimensional, the associated dual optimization problem is also finite dimensional (for deconvolution, it corresponds to optimization over the set of trigonometric polynomials).

Our expertise: We have provided in  [120] a sharp analysis of the support recovery property of this class of methods for the case of sparse spikes deconvolution.

Goals: A key bottleneck of these approaches is that, while being finite dimensional, the dual problem necessitates to handle a constraint of polynomial positivity, which is notoriously difficult to manipulate (except in the very particular case of 1-D problems, which is the one exposed in  [96]). A possible, but very costly, methodology is to ressort to Lasserre's SDP representation hierarchy  [142]. We will make use of these approaches and study how restricting the level of the hierarchy (to obtain fast algorithms) impacts the recovery performances (since this corresponds to only computing approximate solutions). We will pay a particular attention to the recovery of 2-D piecewise constant functions (the so-called total variation of functions regularization  [171]), see Figure 3 for some illustrative applications of this method.

First Order Proximal Schemes

$L^2$ proximal methods.

(Participants: G. Peyré, J-D. Benamou, G. Carlier, Jalal Fadili (ENSICaen)) Both sparse regularization problems in imaging (see Section 2.4) and dynamical optimal transport (see Section 2.3) are instances of large scale, highly structured, non-smooth convex optimization problems. First order proximal splitting optimization algorithms have recently gained lots of interest for these applications because they are the only ones capable of scaling to giga-pixel discretizations of images and volumes and at the same time handling non-smooth objective functions. They have been successfully applied to optimal transport  [60], [158], congested optimal transport  [89] and to sparse regularizations (see for instance  [168] and the references therein).

Our expertise: The pioneering work of our team has shown how these proximal solvers can be used to tackle the dynamical optimal transport problem  [60], see also  [158]. We have also recently developed new proximal schemes that can cope with non-smooth composite objectives functions  [168].

Goals: We aim at extending these solvers to a wider class of variational problems, most notably optimization under divergence constraints  [62]. Another subject we are investigating is the extension of these solvers to both non-smooth and non-convex objective functionals, which are mandatory to handle more general transportation problems and novel imaging regularization penalties.

Bregman proximal methods.

(Participants: G. Peyré G. Carlier, L. Nenna, J-D. Benamou, L. Nenna, Marco Cuturi (Kyoto Univ.)) The entropic regularization of the Kantorovich linear program for OT has been shown to be surprisingly simple and efficient, in particular for applications in machine learning  [118]. As shown in  [64], this is a special instance of the general method of Bregman iterations, which is also a particular instance of first order proximal schemes according to the Kullback-Leibler divergence.

Our expertise: We have recently  [64] shown how Bregman projections  [80] and Dykstra algorithm  [56] offer a generic optimization framework to solve a variety of generalized OT problems. Carlier and Dupuis  [101] have designed a new method based on alternate Dykstra projections and applied it to the principal-agent problem in microeconomics. We have applied this method in computer graphics in a paper accepted in SIGGRAPH 2015  [176]. Figure 9 shows the potential of our approach to handle giga-voxel datasets: the input volumetric densities are discretized on a ${100}^3$ computational grid.

Goals: Following some recent works (see in particular  [110]) we first aim at studying primal-dual optimization schemes according to Bregman divergences (that would go much beyond gradient descent and iterative projections), in order to offer a versatile and very effective framework to solve variational problems involving OT terms. We then also aim at extending the scope of usage of this method to applications in quantum mechanics (Density Functional Theory, see  [114]) and fluid dynamics (Brenier's weak solutions of the incompressible Euler equation, see  [81]). The computational challenge is that realistic physical examples are of a huge size not only because of the space discretization of one marginal but also because of the large number of marginals involved (for incompressible Euler the number of marginals equals the number of time steps).