Primal Dual Methods for Wasserstein Gradient Flows

Combining the classical theory of optimal transport with modern operator splitting techniques, we develop a new numerical method for nonlinear, nonlocal partial differential equations, arising in models of porous media, materials science, and biological swarming. Our method proceeds as follows: first, we discretize in time, either via the classical JKO scheme or via a novel Crank–Nicolson-type method we introduce. Next, we use the Benamou–Brenier dynamical characterization of the Wasserstein distance to reduce computing the solution of the discrete time equations to solving fully discrete minimization problems, with strictly convex objective functions and linear constraints. Third, we compute the minimizers by applying a recently introduced, provably convergent primal dual splitting scheme for three operators (Yan in J Sci Comput 1–20, 2018). By leveraging the PDEs’ underlying variational structure, our method overcomes stability issues present in previous numerical work built on explicit time discretizations, which suffer due to the equations’ strong nonlinearities and degeneracies. Our method is also naturally positivity and mass preserving and, in the case of the JKO scheme, energy decreasing. We prove that minimizers of the fully discrete problem converge to minimizers of the spatially continuous, discrete time problem as the spatial discretization is refined. We conclude with simulations of nonlinear PDEs and Wasserstein geodesics in one and two dimensions that illustrate the key properties of our approach, including higher-order convergence our novel Crank–Nicolson-type method, when compared to the classical JKO method.

     

Clustering via finite nonparametric ICA mixture models

We propose a novel extension of nonparametric multivariate finite mixture models by dropping the standard conditional independence assumption and incorporating the independent component analysis (ICA) structure instead. This innovation extends nonparametric mixture model estimation methods to situations in which conditional independence, a necessary assumption for the unique identifiability of the parameters in such models, is clearly violated. We formulate an objective function in terms of penalized smoothed Kullback–Leibler distance and introduce the nonlinear smoothed majorization-minimization independent component analysis algorithm for optimizing this function and estimating the model parameters. Our algorithm does not require any labeled observations a priori; it may be used for fully unsupervised clustering problems in a multivariate setting. We have implemented a practical version of this algorithm, which utilizes the FastICA algorithm, in the R package icamix. We illustrate this new methodology using several applications in unsupervised learning and image processing.

     

An Interpolating Distance Between Optimal Transport and Fisher–Rao Metrics

This paper defines a new transport metric over the space of nonnegative measures. This metric interpolates between the quadratic Wasserstein and the Fisher–Rao metrics and generalizes optimal transport to measures with different masses. It is defined as a generalization of the dynamical formulation of optimal transport of Benamou and Brenier, by introducing a source term in the continuity equation. The influence of this source term is measured using the Fisher–Rao metric and is averaged with the transportation term. This gives rise to a convex variational problem defining the new metric. Our first contribution is a proof of the existence of geodesics (i.e., solutions to this variational problem). We then show that (generalized) optimal transport and Hellinger metrics are obtained as limiting cases of our metric. Our last theoretical contribution is a proof that geodesics between mixtures of sufficiently close Dirac measures are made of translating mixtures of Dirac masses. Lastly, we propose a numerical scheme making use of first-order proximal splitting methods and we show an application of this new distance to image interpolation.     

Iterated hard-thresholding for linear inverse problems with sparsity constraints

We describe an iterative algorithm for the minimization of Tikhonov type functionals which involve sparsity constraints in form of ℓ^p -penalties which have been proposed recently for the regularization of ill-posed problems. In contrast to the well-known algorithm considered by Daubechies, Defrise and De Mol, it uses hard instead of soft thresholding. This hard thresholding algorithm is based on the generalized conditional gradient method. General results on the convergence of the generalized conditional gradient method enable us to prove strong convergence of the iterates. Furthermore we are able to establish convergence rates of O (n^–1/2) and O (λⁿ) for p = 1 and 1 < p ≤ 2 respectively. (© 2008 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

Alternate step gradient method*

The Barzilai and Borwein (BB) gradient method does not guarantee a descent in the objective function at each iteration, but performs better than the classical steepest descent (SD) method in practice. So far, the BB method has found many successful applications and generalizations in linear systems, unconstrained optimization, convex-constrained optimization, stochastic optimization, etc. In this article, we propose a new gradient method that uses the SD and the BB steps alternately. Hence the name “alternate step (AS) gradient method.” Our theoretical and numerical analyses show that the AS method is a promising alternative to the BB method for linear systems. Unconstrained optimization algorithms related to the AS method are also discussed. Particularly, a more efficient gradient algorithm is provided by exploring the idea of the AS method in the GBB algorithm by Raydan (1997).

To establish a general R-linear convergence result for gradient methods, an important property of the stepsize is drawn in this article. Consequently, R-linear convergence result is established for a large collection of gradient methods, including the AS method. Some interesting insights into gradient methods and discussion about monotonicity and nonmonotonicity are also given.     

Sufficient Descent Riemannian Conjugate Gradient Methods

This paper considers sufficient descent Riemannian conjugate gradient methods with line search algorithms. We propose two kinds of sufficient descent nonlinear conjugate gradient method and prove that these methods satisfy the sufficient descent condition on Riemannian manifolds. One is a hybrid method combining a Fletcher–Reeves-type method with a Polak–Ribière–Polyak-type method, and the other is a Hager–Zhang-type method, both of which are generalizations of those used in Euclidean space. Moreover, we prove that the hybrid method has a global convergence property under the strong Wolfe conditions and the Hager–Zhang-type method has the sufficient descent property regardless of whether a line search is used or not. Further, we review two kinds of line search algorithm on Riemannian manifolds and numerically compare our generalized methods by solving several Riemannian optimization problems. The results show that the performance of the proposed hybrid methods greatly depends on the type of line search used. Meanwhile, the Hager–Zhang-type method has the fast convergence property regardless of the type of line search used.

     

A projected Newton method in a Cartesian product of balls

We formulate a locally superlinearly convergent projected Newton method for constrained minimization in a Cartesian product of balls. For discrete-time,N-stage, input-constrained optimal control problems with Bolza objective functions, we then show how the required scaled tangential component of the objective function gradient can be approximated efficiently with a differential dynamic programming scheme; the computational cost and the storage requirements for the resulting modified projected Newton algorithm increase linearly with the number of stages. In calculations performed for a specific control problem with 10 stages, the modified projected Newton algorithm is shown to be one to two orders of magnitude more efficient than a standard unscaled projected gradient method.This work was supported by the National Science Foundation, Grant No. DMS-85-03746.     

Simplified versions of the conditional gradient method

ABSTRACT

We suggest simple modifications of the conditional gradient method for smooth optimization problems, which maintain the basic convergence properties, but reduce the implementation cost of each iteration essentially. Namely, we propose an adaptive step-size procedure without any line-search and inexact solution of the direction finding subproblem. Preliminary results of computational tests confirm efficiency of the proposed modifications.     

Evolution Models for Mass Transportation Problems

We present a survey on several mass transportation problems, in which a given mass dynamically moves from an initial configuration to a final one. The approach we consider is the one introduced by Benamou and Brenier in [5], where a suitable cost functional F(??, v), depending on the density ?? and on the velocity v (which fulfill the continuity equation), has to be minimized. Acting on the functional F various forms of mass transportation problems can be modeled, as for instance those presenting congestion effects, occurring in traffic simulations and in crowd motions, or concentration effects, which give rise to branched structures.     

A Survey on dynamical transport distances

In this paper, we review some transport models based on the continuity equation, starting with the so-called Benamou − Brenier formula, which is nothing but a fluid mechanics reformulation of the Monge − Kantorovich problem with cost c(x, y) = |x − y|². We discuss some of its applications (gradient flows, sharp functional inequalities, etc.), as well as some variants and generalizations to dynamical transport problems, where interaction effects among mass particles are considered. Bibliography: 43 titles.     

An optimal control problem in image processing

As a starting point, we present a control problem in mammographic image processing which leads to non-standard penalty terms and involves a degenerate parabolic PDE which has to be controlled in the coefficients. We then discuss the classical conditional gradient method from constrained optimization and propose a generalization for non-convex functionals which covers the conditional gradient method as well as the recently proposed iterative shrinkage method of Daubechies, Defrise and De Mol for the solution of linear inverse problems with sparsity promoting penalty terms. We prove that this new algorithm converges. This also gives a deeper understanding of the iterative shrinkage method. Further, we show an application to the above-mentioned control problem in image processing. (© 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

On the use of the energy norm in trust-region and adaptive cubic regularization subproblems

We propose a multi-time scale quasi-Newton based smoothed functional (QN-SF) algorithm for stochastic optimization both with and without inequality constraints. The algorithm combines the smoothed functional (SF) scheme for estimating the gradient with the quasi-Newton method to solve the optimization problem. Newton algorithms typically update the Hessian at each instant and subsequently (a) project them to the space of positive definite and symmetric matrices, and (b) invert the projected Hessian. The latter operation is computationally expensive. In order to save computational effort, we propose in this paper a quasi-Newton SF (QN-SF) algorithm based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update rule. In Bhatnagar (ACM TModel Comput S. 18(1): 27–62, 2007), a Jacobi variant of Newton SF (JN-SF) was proposed and implemented to save computational effort. We compare our QN-SF algorithm with gradient SF (G-SF) and JN-SF algorithms on two different problems – first on a simple stochastic function minimization problem and the other on a problem of optimal routing in a queueing network. We observe from the experiments that the QN-SF algorithm performs significantly better than both G-SF and JN-SF algorithms on both the problem settings. Next we extend the QN-SF algorithm to the case of constrained optimization. In this case too, the QN-SF algorithm performs much better than the JN-SF algorithm. Finally we present the proof of convergence for the QN-SF algorithm in both unconstrained and constrained settings.     

Sur quelques limites de la physique des particules chargées vers la (magnéto)hydrodynamiqueOn some limits in charged particle physics towards (magneto)hydrodynamic equations

We discuss the connection between different scalings limits of the quantum-relativistic Dirac–Maxwell system. In particular we give rigorous results for the quasi-neutral/non-relativistic limit of the Vlasov–Maxwell system: we obtain a magneto-hydro-dynamic system when we consider the magnetic field as a non-relativistic effect and we obtain the Euler equation when we see it as a relativistic effect. A mathematical key is the modulated energy method. To cite this article: Y. Brenier et al., C. R. Acad. Sci. Paris, Ser. I 334 (2002) 239–244.     

Implementation of an optimal first-order method for strongly convex total variation regularization

We present a practical implementation of an optimal first-order method, due to Nesterov, for large-scale total variation regularization in tomographic reconstruction, image deblurring, etc. The algorithm applies to μ-strongly convex objective functions with L-Lipschitz continuous gradient. In the framework of Nesterov both μ and L are assumed known—an assumption that is seldom satisfied in practice. We propose to incorporate mechanisms to estimate locally sufficient μ and L during the iterations. The mechanisms also allow for the application to non-strongly convex functions. We discuss the convergence rate and iteration complexity of several first-order methods, including the proposed algorithm, and we use a 3D tomography problem to compare the performance of these methods. In numerical simulations we demonstrate the advantage in terms of faster convergence when estimating the strong convexity parameter μ for solving ill-conditioned problems to high accuracy, in comparison with an optimal method for non-strongly convex problems and a first-order method with Barzilai-Borwein step size selection.     

Représentation du cône polaire des fonctions convexes et applicationsRepresentation of the polar cone of convex functions and applications

Following Y. Brenier, we give a representation of the polar cone of the set K of the gradient of convex functions, implying the set of measure-preserving maps. This can also be formulated in terms of doubly stochastic measures, and has a geometrical characterization.We deduce an Euler–Lagrange equation and regularity results for some minimization problems in the set K. To cite this article: G. Carlier, T. Lachand-Robert, C. R. Acad. Sci. Paris, Ser. I 335 (2002) 571–576.     

Linear Convergence Rates for Variants of the Alternating Direction Method of Multipliers in Smooth Cases

In the present paper, we propose a novel convergence analysis of the alternating direction method of multipliers, based on its equivalence with the overrelaxed primal–dual hybrid gradient algorithm. We consider the smooth case, where the objective function can be decomposed into one differentiable with Lipschitz continuous gradient part and one strongly convex part. Under these hypotheses, a convergence proof with an optimal parameter choice is given for the primal–dual method, which leads to convergence results for the alternating direction method of multipliers. An accelerated variant of the latter, based on a parameter relaxation, is also proposed, which is shown to converge linearly with same asymptotic rate as the primal–dual algorithm.     

Sample size selection in optimization methods for machine learning

This paper presents a methodology for using varying sample sizes in batch-type optimization methods for large-scale machine learning problems. The first part of the paper deals with the delicate issue of dynamic sample selection in the evaluation of the function and gradient. We propose a criterion for increasing the sample size based on variance estimates obtained during the computation of a batch gradient. We establish an complexity bound on the total cost of a gradient method. The second part of the paper describes a practical Newton method that uses a smaller sample to compute Hessian vector-products than to evaluate the function and the gradient, and that also employs a dynamic sampling technique. The focus of the paper shifts in the third part of the paper to L ₁-regularized problems designed to produce sparse solutions. We propose a Newton-like method that consists of two phases: a (minimalistic) gradient projection phase that identifies zero variables, and subspace phase that applies a subsampled Hessian Newton iteration in the free variables. Numerical tests on speech recognition problems illustrate the performance of the algorithms.     

Riemannian proximal gradient methods

In the Euclidean setting the proximal gradient method and its accelerated variants are a class of efficient algorithms for optimization problems with decomposable objective. In this paper, we develop a Riemannian proximal gradient method (RPG) and its accelerated variant (ARPG) for similar problems but constrained on a manifold. The global convergence of RPG is established under mild assumptions, and the O(1/k) is also derived for RPG based on the notion of retraction convexity. If assuming the objective function obeys the Rimannian Kurdyka–?ojasiewicz (KL) property, it is further shown that the sequence generated by RPG converges to a single stationary point. As in the Euclidean setting, local convergence rate can be established if the objective function satisfies the Riemannian KL property with an exponent. Moreover, we show that the restriction of a semialgebraic function onto the Stiefel manifold satisfies the Riemannian KL property, which covers for example the well-known sparse PCA problem. Numerical experiments on random and synthetic data are conducted to test the performance of the proposed RPG and ARPG.

     

Shape inverse problem for Stokes–Brinkmann equations

In this paper, we propose an imaging technique for the detection of porous inclusions in a stationary flow governed by Stokes–Brinkmann equations. We introduce the velocity method to perform the shape deformation, and derive the structure of shape gradient for the cost functional based on the continuous adjoint method and the function space parametrization technique. Moreover, we present a gradient-type algorithm to the shape inverse problem. The numerical results demonstrate the proposed algorithm is feasible and effective for the quite high Reynolds numbers problems.     

Block coordinate descent for smooth nonconvex constrained minimization

We present a new algorithm for large-scale unconstrained minimization that, at each iteration, minimizes, approximately, a quadratic model of the objective function plus a regularization term, not necessarily based on a norm. We prove convergence assuming only gradient continuity and complexity results assuming Lipschitz conditions. For solving the subproblems in the case of regularizations based on the 3-norm, we introduce a new method that quickly obtains the approximate solutions required by the theory. We present numerical experiments.