# Augmented Lagrangian Method for Optimal Partial Transportation

Augmented Lagrangian Method for Optimal Partial Transportation Abstract The use of augmented Lagrangian algorithm for optimal transport problems goes back to Benamou & Brenier (2000, Acomputational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math., 84, 375–393), in the case where the cost corresponds to the square of the Euclidean distance. It was recently extended in Benamou & Carlier (2015, Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations. J. Optim. Theory Appl., 167, 1–26), to the optimal transport with the Euclidean distance and Mean-Field Games theory and in Benamou et al. (2017, A numerical solution to Monge’s problem with a Finsler distance cost ESAIM: M2AN), to the optimal transportation with Finsler distances. Our aim here is to show how one can use this method to study the optimal partial transport problem with Finsler distance costs. To this aim, we introduce a suitable dual formulation of the optimal partial transport, which contains all the information on the active regions and the associated flow. Then, we use a finite element discretization with the FreeFem++ software to provide numerical simulations for the optimal partial transportation. A convergence study for the potential together with the flux and the active regions is given to validate the approach 1. Introduction The theory of optimal transportation deals with the problem to find the optimal way to move materials from a given source to a desired target in such a way to minimize the work. The problem was first proposed and studied by G. Monge in 1781, and then L. Kantorovich made fundamental contributions to the problem in the 1940s by relaxing the problem into a linear one. Since the late 1980s, this subject has been investigated under various points of view with many applications in image processing, geometry, probability theory, economics, evolution partial differential equations (PDEs) and other areas. For more information on the optimal mass transport problem, we refer the reader to the pedagogical books (Villani, 2003; Ambrosio et al., 2005; Villani, 2009; Santambrogio, 2015). The standard optimal transport problem requires that the total mass of the source is equal to the total mass of the target (balance condition of mass) and that all the materials of the source must be transported. Here, we are interested in the optimal partial transportation. That is the case where the balance condition of mass is excluded, and the aim is to transport effectively a prescribed amount of mass from the source to the target. In other words, the optimal partial transport problem aims to study the practical situation, where only a part of the commodity (respectively, consumer demand) of a prescribed total mass $$\mathbf{m}$$ needs to be transported (respectively, fulfilled). This generalized problem brings out additional variables called active regions. The problem was first studied theoretically in Caffarelli & McCann (2010) (see also Figalli, 2010) in the case where the work is proportional to the square of the Euclidean distance. Recently in Igbida & Nguyen (2017), we give a complete theoretical study of the problem in the case where the work is proportional to a Finsler distance $$d_{F}$$ (covering by the way the case of the Euclidean distance), where $$d_{F}$$ is given as follows (see Section 2)   dF(x,y):=infξ∈Lip([0,1];Ω¯){∫01F(ξ(t),ξ˙(t))dt:ξ(0)=x,ξ(1)=y}. Concerning numerical approximations for the optimal partial transport, Barrett & Prigozhin (2009) studied the case of the Euclidean distance by using an approximation based on nonlinear approximated PDEs and Raviart–Thomas finite elements. Benamou et al. (2015) and Chizat et al. (2016) introduced general numerical frameworks to approximate solutions to linear programs related to the optimal transport (including the optimal partial transport). Their idea is based on an entropic regularization of the initial linear programs. This is a static approach to optimal transport-type problems and needs to use (approximated) values of $$d_{F}(x, y)$$. In this article, we use a different approach (based mainly on Benamou & Brenier, 2000; Benamou & Carlier, 2015; Igbida & Nguyen, 2017) to compute the solution of the optimal partial transport problem. We first show how one can directly reformulate the unknown quantities (variables) of the optimal partial transport into an infinite-dimensional minimization problem of the form:   minϕ∈VF(ϕ)+G(Λϕ), where $$\mathcal{F}, \mathcal{G}$$ are l.s.c., convex functionals and $${\it{\Lambda}}\in \mathcal{L}(V, Z)$$ is a continuous linear operator between two Banach spaces. Thanks to peculiar properties of $$\mathcal{F}$$ and $$\mathcal{G}$$ in our situation, an augmented Lagrangian method is applied effectively in the same spirit of Benamou & Carlier (2015) and Benamou et al. (2017). We show that, for the computation, we just need to solve linear equations (with a symmetric positive definite coefficient matrix) or to update explicit formulations. It is worth noting that this method uses only elementary operations without evaluating $$d_{F}$$. This article is organized as follows: In the next section, we introduce the optimal partial transport problem and its equivalent formulations with a particular attention to the Kantorovich dual formulation. In Section 3, we give a finite-dimensional approximation of the problem, and show that primal-dual solutions of the discretized problems converge to the ones of original continuous problems. The details of the ALG2 algorithm is given in Section 4. Some numerical examples are presented in Section 5. We terminate the article by an appendix, where we give proofs of some facts we need in this article. 2. Partial transport and its equivalent formulations Let $${\it{\Omega}}$$ be a connected bounded Lipschitz domain and $$F$$ be a continuous Finsler metric on $$\overline{{\it{\Omega}}}$$, i.e., $$F: \overline{{\it{\Omega}}}\times \mathbb{R}^N \longrightarrow [0, +\infty)$$ is continuous and $$F(x, .)$$ is convex, positively homogeneous of degree $$1$$ in the sense   F(x,tv)=tF(x,v)∀t>0,v∈RN. We assume moreover that $$F$$ is nondegenerate in the sense that there exist positive constants $$M_{1}, M_{2}$$, such that   M1|v|≤F(x,v)≤M2|v|∀x∈Ω¯,v∈RN. Let $$\mu, \nu\in \mathcal M^+_{b}(\overline{{\it{\Omega}}})$$ be two Radon measures on $$\overline{{\it{\Omega}}}$$ and $$\mathbf{m}_{\max}:=\min\{\mu(\overline{{\it{\Omega}}}), \nu(\overline{{\it{\Omega}}})\}.$$ Given a total mass $$\mathbf{m}\in [0, \mathbf{m}_{\max}]$$, the optimal partial transport problem (or partial Monge–Kantorovich problem, PMK for short) aims to transport effectively the total mass $$\mathbf{m}$$ from a supply subregion of the source $$\mu$$ into a subregion of the target $$\nu$$. The set of subregions of mass $$\mathbf{m}$$ is given by   Subm(μ,ν):={(ρ0,ρ1)∈Mb+(Ω¯)×Mb+(Ω¯):ρ0≤μ,ρ1≤ν,ρ0(Ω¯)=ρ1(Ω¯)=m}. An element $$(\rho_{0}, \rho_{1})\in Sub_{\mathbf{m}}(\mu, \nu)$$ is called a couple of active regions. As for the optimal transport, one can work with different kinds of cost functions for the optimal partial transport, i.e., in the formulation (2.1) below, $$d_{F}(x, y)$$ can be replaced by a general measurable cost function $$c(x, y)$$. However, in this article, we focus on the case where the cost $$c=d_{F}$$. So let us state the problem directly for $$d_{F}$$. The PMK problem (Barrett & Prigozhin, 2009; Caffarelli & McCann, 2010; Figalli, 2010; Igbida & Nguyen, 2017) aims to minimize the following problem   min{K(γ):=∫Ω¯×Ω¯dF(x,y)dγ:γ∈πm(μ,ν)}, (2.1) where $$d_{F}$$ is the Finsler distance on $$\overline{{\it{\Omega}}}$$ associated with $$F$$, i.e.,   dF(x,y):=inf{∫01F(ξ(t),ξ˙(t))dt:ξ(0)=x,ξ(1)=y,ξ∈Lip([0,1];Ω¯)}. $${\large{\mathbf{\pi}}}_{\mathbf{m}}(\mu, \nu)$$ is the set of transport plans of mass $$\mathbf{m}$$, i.e.,   πm(μ,ν):={γ∈Mb+(Ω¯×Ω¯):(πx#γ,πy#γ)∈Subm(μ,ν)}. Here, $$\pi_x \#\gamma$$ and $$\: \pi_y \# \gamma$$ are the first and second marginals of $$\gamma$$. An optimal $$\gamma^*$$ is called an optimal plan and $$(\pi_x \#\gamma^*, \pi_y \# \gamma^*)$$ is called a couple of optimal active regions. Following Igbida & Nguyen (2017), to study the PMK problem we use its dual problem that we call the dual partial Monge–Kantorovich problem (DPMK). To this aim, we consider $$Lip_{d_{F}}$$ the set of $$1-$$Lipschitz continuous functions w.r.t. $$d_{F}$$ given by   LipdF:={u:Ω¯⟶R|u(y)−u(x)≤dF(x,y)∀x,y∈Ω¯}. Then, the connection between the PMK and DPMK problems is summarized in the following theorem. Theorem 2.1 Let $$\mu, \nu\,{\in}\,\mathcal M^+_b(\overline{{\it{\Omega}}})$$ be Radon measures and $$\mathbf{m}\,{\in}\, [0,\mathbf{m}_{\max}]$$. The partial Monge–Kantorovich problem has a solution $$\sigma^*\,{\in}\, {\large{\mathbf{\pi}}}_{\mathbf{m}}(\mu, \nu)$$ and   K(σ∗)=max{D(λ,u):=∫Ω¯ud(ν−μ)+λ(m−ν(Ω¯)):λ≥0 and u∈LdFλ}, (2.2) where   LdFλ:={u∈LipdF:0≤u(x)≤λ for any x∈Ω¯}. Moreover, $$\sigma\in {\large{\mathbf{\pi}}}_{\mathbf{m}}(\mu, \nu)$$ and $$(\lambda, u)\in \mathbb{R}^+\times L^\lambda_{d_{F}}$$ are solutions, respectively, if and only if   u(x)=0 for (μ−πx#σ)-a.e. x∈Ω¯,u(x)=λ for (ν−πy#σ)-a.e. x∈Ω¯and u(y)−u(x)=dF(x,y) for σ-a.e. (x,y)∈Ω¯×Ω¯. Proof. The proof follows in the same way of Theorem 2.4 in Igbida & Nguyen (2017), where the authors study the case $${\it{\Omega}}=\mathbb{R}^N.$$ □ The DPMK problem (2.2) contains all the information concerning the optimal partial mass transportation. However, for the numerical approximation of the optimal partial transportation and to use the augmented Lagrangian method, we need to rewrite the problem into the form   infϕ∈VF(ϕ)+G(Λϕ). To do that, we consider the polar function $$F^*$$ of $$F,$$ which is defined by   F∗(x,p):=sup{⟨v,p⟩:F(x,v)≤1} for x∈Ω¯,p∈RN. Note that $$F^*(x, .)$$ is not the Legendre–Fenchel transform. It is easy to see that $$F^*$$ is also a continuous, nondegenerate Finsler metric on $$\overline{{\it{\Omega}}}$$ and   ⟨v,p⟩≤F∗(x,p)F(x,v)∀x∈Ω¯,v,p∈RN. Remark 2.2 Using the polar function $$F^*$$, we can characterize the set $$Lip_{d_{F}}$$ as (see the appendix if necessary)   LipdF={u:Ω¯⟶R|u is Lipschitz continuous and F∗(x,∇u(x))≤1 a.e. x∈Ω}. Thanks to this remark, the DPMK problem (2.2) can be written as   max{D(λ,u):0≤u(x)≤λ,u is Lipschitz continuous, F∗(x,∇u(x))≤1 a.e. x∈Ω}. Moreover, we have Theorem 2.3 Under the assumptions of Theorem 2.1, setting $$V:=\mathbb{R}\times C^1(\overline{{\it{\Omega}}})$$ and $$Z:=C(\overline{{\it{\Omega}}})^{N}\times C(\overline{{\it{\Omega}}})\times C(\overline{{\it{\Omega}}}),$$ we have   K(σ∗)=−inf{F(λ,u)+G(Λ(λ,u)):(λ,u)∈V}, (2.3) where $${\it{\Lambda}}\in \mathcal{L}(V, Z)$$ is given by   Λ(λ,u):=(∇u,−u,u−λ)∀(λ,u)∈V, and $$\mathcal{F}: V\longrightarrow (-\infty, +\infty]$$, $$\mathcal{G}: Z\longrightarrow (-\infty, +\infty]$$ are the l.s.c. convex functions given by   F(λ,u) :=−∫Ω¯ud(ν−μ)−λ(m−ν(Ω¯))∀(λ,u)∈V;G(q,z,w) :={0if z(x)≤0,w(x)≤0,F∗(x,q(x))≤1∀x∈Ω¯+∞otherwise  for (q,z,w)∈Z. To prove this theorem, we need the following lemma. Lemma 2.4 Let $$\lambda\,{\ge}\, 0$$ be fixed. For any $$u\,{\in}\, L^\lambda_{d_{F}}$$, there exists a sequence of smooth functions $$u_{\varepsilon}\in C^\infty_{c}(\mathbb{R}^N) \bigcap L^\lambda_{d_{F}}$$ such that $$u_{\varepsilon} \rightrightarrows u$$ uniformly on $$\overline{{\it{\Omega}}}$$. The result of the lemma is more or less known in some cases (see Igbida & Ta Thi, 2017 for the case where the function $$u$$ is null on the boundary). The proof in the general case is quite technical and will be given in the appendix. Proof of Theorem 2.3. Thanks to Remark 2.2 and Lemma 2.4, we have   −inf(λ,u)∈VF(λ,u)+G(Λ(λ,u)) =sup{∫Ω¯ud(ν−μ)+λ(m−ν(Ω¯)):λ≥0,u∈C1(Ω¯)∩LdFλ} =max{D(λ,u):λ≥0 and u∈LdFλ}. Using the duality (2.2), the proof is completed. □ To end this section, we prove the following result that will be useful for the proof of the convergence of our discretization. Theorem 2.5 Under the assumptions of Theorem 2.1, we have   −inf(λ,u)∈VF(λ,u)+G(Λ(λ,u))=min{∫Ω¯F(x,Φ|Φ|(x))d|Φ|:(Φ,θ0,θ1)∈Ψm(μ,ν)}, (2.4) where   Ψm(μ,ν):={(Φ,θ0,θ1)∈Z∗=Mb(Ω¯)N×Mb(Ω¯)×Mb(Ω¯):θ0≥0,θ1≥0,θ1(Ω¯)=ν(Ω¯)−m  and −∇⋅Φ=ν−θ1−(μ−θ0) with Φ.n=0 on ∂Ω}. Actually, the minimal flow-type formulation   min{∫Ω¯F(x,Φ|Φ|(x))d|Φ|:(Φ,θ0,θ1)∈Ψm(μ,ν)} (2.5) introduces the Beckmann problem (see Beckmann, 1952) for the optimal partial transport with Finsler distance costs. See here that in the balanced case, i.e., $$\mathbf{m}=\mu(\overline{{\it{\Omega}}})=\nu(\overline{{\it{\Omega}}})$$, the formulation (2.5) becomes   min{∫Ω¯F(x,Φ|Φ|(x))d|Φ|:Φ∈Mb(Ω¯)N,−∇⋅Φ=ν−μ with Φ.n=0 on ∂Ω}. (2.6) An optimal solution $${\it{\Phi}}$$ of the problem (2.6) is called an optimal flow of transporting $$\mu$$ onto $$\nu$$. As known from the optimal transport theory, the optimal flow gives a way to visualize the transportation. To prove Theorem 2.5, we will use the well-known duality arguments. For convenience, let us recall here the Fenchel–Rockafellar duality. Let us consider the problem   infϕ∈VF(ϕ)+G(Λϕ), (2.7) where $$\mathcal{F}: V\,{\longrightarrow}\, (-\infty, +\infty]$$ and $$\mathcal{G}: Z\,{\longrightarrow}\, (-\infty, +\infty]$$ are convex, l.s.c. and $${\it{\Lambda}}\,{\in}\, \mathcal{L}(V, Z)$$ the space of linear continuous functions from $$V$$ to $$Z.$$ Using $$\mathcal{F}^{*}$$ and $$\mathcal{G}^{*}$$ the conjugate functions (given by the Legendre–Fenchel transformation) of $$\mathcal{F}$$ and $$\mathcal{G}$$, respectively, and $${\it{\Lambda}}^{*}$$ is the adjoint operator of $${\it{\Lambda}},$$ it is not difficult to see that   supσ∈Z∗(−F∗(−Λ∗σ)−G∗(σ))≤infϕ∈VF(ϕ)+G(Λϕ), where $$Z^*$$ is the topological dual space associated with $$Z$$. This is the so called weak duality. For the strong duality, which corresponds to equality we have the following well-known result. Proposition 2.6 (cf. Ekeland & Teman, 1976) In addition, assume that there exists $$\phi_{0}$$ such that $$\mathcal{F}(\phi_0)\,{<}+\infty$$, $$\mathcal{G}({\it{\Lambda}} \phi_{0}) <+\infty,$$$$\mathcal{G}$$ being continuous at $${\it{\Lambda}} \phi_0$$. Then the Fenchel–Rockafellar dual problem   supσ∈Z∗(−F∗(−Λ∗σ)−G∗(σ)) (2.8) has at least a solution $$\sigma\in Z^*$$ and $$\inf$$ (2.7) = $$\max$$ (2.8). Moreover, in this case, $$\phi$$ is a solution to the primal problem (2.7) if and only if   −Λ∗σ∈∂F(ϕ) and σ∈∂G(Λϕ). (2.9) Proof of Theorem 2.5. We work with the uniform convergence for the spaces $$C(\overline{{\it{\Omega}}})^N$$, $$C(\overline{{\it{\Omega}}})$$ and the norm $$\|u\|_{C^1}:=\max\{\|u\|_{\infty}, \|\nabla u\|_{\infty}\}$$ for $$C^1(\overline{{\it{\Omega}}})$$. It is not difficult to see that the hypotheses of Proposition 2.6 are satisfied. Now, let us compute the Fenchel–Rockafellar dual problem of (2.3). Since $$\mathcal{F}$$ is linear, $${\mathcal{F}}^*(-{\it{\Lambda}}^*({\it{\Phi}}, \theta^0, \theta^1))$$ is finite (and always equals to $$0$$) if and only if   −Λ∗(Φ,θ0,θ1)=−(m−ν(Ω¯),ν−μ) in V∗ i.e.,   ⟨Φ,∇u⟩−⟨θ0,u⟩+⟨θ1,u−λ⟩=λ(m−ν(Ω¯))+⟨ν−μ,u⟩∀(λ,u)∈V. This implies that   ∫Ω¯∇udΦ=∫Ω¯ud(ν−θ1)−∫Ω¯ud(μ−θ0) for all u∈C1(Ω¯) and   −λ∫Ω¯dθ1=λ(m−ν(Ω¯))∀λ∈R. These mean that   −∇⋅Φ=ν−θ1−(μ−θ0) with Φ.n=0 on ∂Ω and   θ1(Ω¯)=ν(Ω¯)−m. We also have   G∗(Φ,θ0,θ1)={∫Ω¯F(x,Φ|Φ|(x))d|Φ| if θ0≥0,θ1≥0+∞ otherwise  for any (Φ,θ0,θ1)∈Z∗. Then the proof follows by Proposition 2.6. □ Remark 2.7 The optimality relations (2.9) reads   {−∇⋅Φ=ν−θ1−(μ−θ0) and Φ⋅n=0 on ∂Ωθ1(Ω¯)=ν(Ω¯)−m⟨Φ,∇u⟩≥⟨Φ,q⟩∀q∈C(Ω¯),F∗(x,q(x))≤1∀x∈Ω¯λ∈R+,u∈C1(Ω¯)⋂LdFλu=0,θ0-a.e. in Ω¯u=λ,θ1-a.e. in Ω¯.  In fact, the optimality condition $$-{\it{\Lambda}}^* \sigma \in \partial \mathcal{F}(\phi)$$ gives the first two equations and $$\sigma\in \partial \mathcal{G}({\it{\Lambda}} \phi)$$ gives the last four equations. Moreover, if $${\it{\Phi}}\in L^1({\it{\Omega}})^N$$ then the condition   ⟨Φ,∇u⟩≥⟨Φ,q⟩∀q∈C(Ω¯),F∗(x,q(x))≤1∀x∈Ω¯ can be replaced by   F(x,Φ(x))=⟨∇u(x),Φ(x)⟩ a.e. x∈Ω. (2.10) However, it is not clear in general that $${\it{\Phi}}$$ belongs to $$L^1({\it{\Omega}})^N$$. In the case where $${\it{\Omega}}$$ is convex and $$F(x, v):=|v|$$ the Euclidean norm (or some other uniformly convex and smooth norms), the $$L^p$$ regularity results are known under suitable assumptions on $$\mu$$ and $$\nu$$ (see, e.g., Feldman & McCann, 2002; De Pascale et al., 2004; De Pascale & Pratelli, 2004; Santambrogio, 2009). To our knowledge, the case of general Finsler metrics is still an open question. In the case where $${\it{\Phi}}$$ is a vector-valued measure, the condition (2.10) should be adapted to the tangential gradient. Rigorous formulations using the tangential gradient with respect to a measure, as well as rigorous proofs in the general case, can be found in the article by Igbida & Nguyen (2017) with $${\it{\Omega}}=\mathbb{R}^N$$. It is expected that $$\theta^0\le \mu$$ and $$\theta^1\le \nu$$ for optimal solutions $$({\it{\Phi}}, \theta^0, \theta^1)$$ of the minimal flow formulation (2.5). This is the case whenever $$\mathbf{m}\in [(\mu\wedge \nu)(\overline{{\it{\Omega}}}), {\mathbf{m}}_{\max}]$$, where $$\mu\wedge \nu$$ is the common mass measure of $$\mu$$ and $$\nu$$, i.e., if $$\mu, \nu\in L^1({\it{\Omega}})$$, then $$\mu \wedge \nu\in L^1({\it{\Omega}})$$ and   (μ∧ν)(x)=min{μ(x),ν(x)} for a.e. x∈Ω. In general, the measure $$\mu \wedge \nu$$ is defined by (see Ambrosio et al., 2000)   μ∧ν(A)=inf{μ(A1)+ν(A2):disjoint Borel setsA1,A2,such that A1∪A2=A}. Proposition 2.8 Let $$\mathbf{m}\in [(\mu\wedge \nu)(\overline{{\it{\Omega}}}), {\mathbf{m}}_{\max}]$$ and $$({\it{\Phi}}, \theta^0, \theta^1)\in Z^*$$ be an optimal solution of (2.5). Then $$\theta^0\le \mu$$ and $$\theta^1\le \nu$$. Moreover, $$(\mu-\theta^0, \nu -\theta^1)$$ is a couple of optimal active regions and $${\it{\Phi}}$$ is an optimal flow of transporting $$\mu-\theta^0$$ onto $$\nu -\theta^1$$. Proof. The proof follows in the same way as Theorem 5.21 and Corollary 5.20 in Igbida & Nguyen (2017). □ Our next work is to compute an approximation of $${\it{\Phi}}$$ (in fact, approximations of $${\it{\Phi}}, u, \lambda, \theta^0, \theta^1$$). To do that, we will apply an augmented Lagrangian method to the DPMK problem (2.2). 3. Discretization and convergence Coming back to the DPMK problem (2.2), our aim now is to give, by using a finite element approximation, the discretized problem associated with (2.2). To begin with, let us consider regular triangulations $$\mathcal{T}_{h}$$ of $$\overline{{\it{\Omega}}}$$. For a fixed integer $$k\ge 1$$, $$P_{k}$$ is the set of polynomials of degree less or equal $$k$$. Let $$E_{h}\subset H^{1}({\it{\Omega}})$$ be the space of continuous functions on $$\overline{{\it{\Omega}}}$$ and belonging to $$P_{k}$$ on each triangle of $$\mathcal{T}_{h}$$. We denote by $$Y_{h}$$ the space of vectorial functions such that their restrictions belong to $$(P_{k-1})^N$$ on each triangle of $$\mathcal{T}_{h}$$. Let $$f=\nu-\mu$$ and $$f_{h}\in E_{h}$$ such that $$\{f_{h}\}$$ converges weakly* to $$f$$ in $$\mathcal{M}_{b}(\overline{{\it{\Omega}}})$$. Considering the finite-dimensional spaces   Vh=R×Eh,Zh=Yh×Eh×Eh, we set   Λh(λ,u) :=(∇u,−u,u−λ)∈Zh for (λ,u)∈Vh,Fh(λ,u) :=−⟨u,fh⟩−λ(m−ν(Ω¯))∀(λ,u)∈Vh and   Gh(q,z,w):={0if z≤0,w≤0,F∗(x,q(x))≤1 for a.e. x∈Ω+∞otherwise  for (q,z,w)∈Zh. Then the finite-dimensional approximation of (2.2) reads   inf(λ,u)∈VhFh(λ,u)+Gh(Λh(λ,u)). (3.1) The following result shows that this is a suitable approximation of (2.2). Theorem 3.1 Assume that $$\mathbf{m}\,{<}\, \nu(\overline{{\it{\Omega}}})$$. Let $$(\lambda_{h}, u_{h})\,{\in}\, V_{h}$$ be an optimal solution to the approximated problem (3.1) and $$({\it{\Phi}}_{h}, \theta^0_h, \theta^1_h)$$ be an optimal dual solution to (3.1). Then, up to a subsequence, $$(\lambda_{h}, u_h)$$ converges in $$\mathbb{R}\times C(\overline {\it{\Omega}})$$ to $$(\lambda, u)$$ an optimal solution of the DPMK problem (2.2) and $$({\it{\Phi}}_h, \theta^0_h, \theta^1_h)$$ converges weakly* in $$\mathcal M_b(\overline{\it{\Omega}})^N\times \mathcal M_b(\overline{\it{\Omega}})\times \mathcal M_b(\overline{\it{\Omega}})$$ to $$({\it{\Phi}}, \theta^0, \theta^1)$$ an optimal solution of (2.5). Proof. Since $$\mathbf{m}\,{<}\,\nu(\overline{{\it{\Omega}}})$$, $$\{\lambda_h\}$$ is bounded in $$\mathbb{R}$$ and $$\{u_h\}$$ is bounded in $$(C(\overline{{\it{\Omega}}}), \|.\|_{\infty})$$. From the nondegeneracy of $$F$$ and the definitions of $$\mathcal{F}_{h}, \mathcal{G}_{h}, {\it{\Lambda}}_{h}$$, we have that $$\{u_h\}$$ is equi-Lipschitz and   uh(y)−uh(x)≤dF(x,y)∀x,y∈Ω¯. Using the Ascoli–Arzela theorem, up to a subsequence, $$u_h \rightrightarrows u$$ uniformly on $$\overline{{\it{\Omega}}}$$ and $$\lambda_h \to \lambda$$. Obviously, $$\lambda\ge 0$$ and $$u\in L^\lambda_{d_{F}}$$. Now, by the optimality of $$(\lambda_h, u_h)$$ and of $$({\it{\Phi}}_h, \theta^0_h, \theta^1_h)$$, we have   −Λh∗(Φh,θh0,θh1)=−(m−ν(Ω¯),fh) in Vh∗ and   Fh(λh,uh)+Gh(Λh(λh,uh))=−Fh∗(−Λh∗(Φh,θh0,θh1))−Gh∗(Φh,θh0,θh1). More concretely,   ⟨Φh,∇v⟩−⟨θh0,v⟩+⟨θh1,v−s⟩=s(m−ν(Ω¯))+⟨fh,v⟩∀(s,v)∈Vh, (3.2)  θh0≥0,θh1≥0,θh1(Ω¯)=ν(Ω¯)−m (3.3) and   ⟨uh,fh⟩+λh(m−ν(Ω¯))=sup{⟨q,Φh⟩:q∈Yh,F∗(x,q(x))≤1 a.e. x∈Ω}. (3.4) In (3.2), taking $$v=0$$ and $$s=1$$ (respectively, $$v=s=1$$), we see that $$\{\theta^1_h\}$$ (respectively, $$\{\theta^0_h\}$$) is bounded in $$\mathcal M_{b}(\overline{{\it{\Omega}}})$$. Moreover, using (3.4) and the boundedness of $$(\lambda_h, u_h)$$ we deduce that $$\{{\it{\Phi}}_h\}$$ is bounded in $$\mathcal{M}_b(\overline{{\it{\Omega}}})^N.$$ So, up to a subsequence,   (Φh,θh0,θh1)⇀(Φ,θ0,θ1) in Mb(Ω¯)N×Mb(Ω¯)×Mb(Ω¯)−w∗. Using (3.2) and (3.3), it is clear that $$({\it{\Phi}}, \theta^0, \theta^1)$$ satisfies   ⟨Φ,∇v⟩−⟨θ0,v⟩+⟨θ1,v−s⟩=s(m−ν(Ω¯))+⟨f,v⟩∀(s,v)∈V and   θ0≥0,θ1≥0,θ1(Ω¯)=ν(Ω¯)−m, i.e., $$({\it{\Phi}}, \theta^0, \theta^1)$$ is feasible for the minimal flow problem (2.5). Next, let us show the optimality of $$(\lambda, u)$$ and of $$({\it{\Phi}}, \theta^0, \theta^1)$$, i.e.,   ∫Ω¯F(x,Φ|Φ|(x))d|Φ|=⟨u,ν−μ⟩+λ(m−ν(Ω¯)). (3.5) We fix $$q\,{\in}\, C(\overline{{\it{\Omega}}})^N$$ such that $$F^*(x, q(x))\,{\le}\, 1 \quad\forall x\,{\in}\, \overline{{\it{\Omega}}}$$ and we consider $$q_h\,{\in}\, Y_{h}$$ such that $$\|q_{h}-q\|_{L^{\infty}({\it{\Omega}})} \to 0$$ as $$h\to 0$$. We see that   F∗(x,qh(x))=F∗(x,q(x))+F∗(x,qh(x))−F∗(x,q(x))≤1+O(h) a.e. x∈Ω. By taking $$\frac{q_{h}}{1+ O(h)}$$, we can assume that $$q_{h}\,{\in}\, Y_h, F^*(x, q_h(x))\,{\le}\, 1 \text{ a.e. } x\,{\in}\, {\it{\Omega}}$$ and $$\|q_{h}-q\|_{L^{\infty}({\it{\Omega}})} \,{\to}\, 0$$ as $$h\to 0$$. Using (3.4), we have   ⟨q,Φ⟩ =⟨qh,Φh⟩+⟨q,Φ−Φh⟩+⟨q−qh,Φh⟩ ≤sup{⟨qh,Φh⟩:qh∈Yh,F∗(x,qh(x))≤1 a.e. x∈Ω}+O(h) =⟨uh,fh⟩+λh(m−ν(Ω¯))+O(h). Letting $$h\to 0$$, we get   ⟨q,Φ⟩≤⟨u,ν−μ⟩+λ(m−ν(Ω¯)) for any q∈C(Ω¯)N,F∗(x,q(x))≤1∀x∈Ω¯. Taking supremum in $$q$$, we obtain   ∫Ω¯F(x,Φ|Φ|(x))d|Φ|≤⟨u,ν−μ⟩+λ(m−ν(Ω¯)). At last, thanks to the duality equality (2.4), this implies (3.5), the optimality of $$(\lambda, u)$$ and of $$({\it{\Phi}}, \theta^0, \theta^1)$$. □ Remark 3.2 In the case $$\mathbf{m}=\mathbf{m}_{\max}$$ (called the unbalanced transport), the DPMK problem has a simpler formulation. So for the purpose of implementation, we distinguish the two cases: the partial transport and the unbalanced transport. In the unbalanced case, let us assume that $$\mathbf{m}=\mathbf{m}_{\max}=\nu(\overline{{\it{\Omega}}})$$$$(\text{i.e., }\, \mu(\overline{{\it{\Omega}}})\ge \nu(\overline{{\it{\Omega}}}))$$, the DPMK problem (2.2) can be written as   maxu∈LipdF,u≥0∫Ω¯ud(ν−μ). (3.6) By using $$V_h=E_h,\: Z_h=Y_h \times E_h, {\it{\Lambda}}_h u =(\nabla u, -u)$$ and   Gh(q,z)={0 if z≤0,F∗(x,q(x))≤1 a.e. x∈Ω+∞ otherwise  a finite-dimensional approximation can be given by   infu∈Vh−⟨u,fh⟩+Gh(Λhu). (3.7) As in Theorem 3.1, we can prove the convergence of this finite-dimensional approximation to the original one (3.6). More precisely, we have Proposition 3.3 Assume that $$\mathbf{m}=\nu(\overline{{\it{\Omega}}})$$. Let $$u_{h}\in V_{h}$$ be an optimal solution to the approximated problem (3.7) and $$({\it{\Phi}}_{h}, \theta^0_h)$$ be an optimal dual solution to (3.7). Then, up to a subsequence and translation by constant, $$u_{h}$$ converges to $$u$$ an optimal solution of the DPMK problem (3.6) and $$({\it{\Phi}}_h, \theta^0_h)$$ converges to $$({\it{\Phi}}, \theta^0)$$ an optimal solution of (2.5) with $$\theta^1=0$$. The proof of this proposition is similar to the proof of Theorem 3.1. 4. Solving the discretized problems Our task now is to solve the finite-dimensional problems (3.1) and (3.7). First, let us recall the augmented Lagrangian method we are dealing with. 4.1 ALG2 method Assume that $$V$$ and $$Z$$ are two Hilbert spaces. Let us consider the problem   infϕ∈VF(ϕ)+G(Λϕ), (4.1) where $$\mathcal{F}: V\longrightarrow (-\infty, +\infty]$$ and $$\mathcal{G}: Z\longrightarrow (-\infty, +\infty]$$ are convex, l.s.c. and $${\it{\Lambda}}\in \mathcal{L}(V, Z)$$. We introduce a new variable $$q\in Z$$ to the primal problem (4.1) and we rewrite it in the form   inf(ϕ,q)∈V×Z:Λϕ=qF(ϕ)+G(q). The augmented Lagrangian is given by   L(ϕ,q;σ):=F(ϕ)+G(q)+⟨σ,Λϕ−q⟩+r2|Λϕ−q|2r>0. The so called ALG2 algorithm is given as follows: For given $$q_{0}, \sigma_{0}\in Z$$, we construct the sequences $$\{\phi_i\}, \{q_i\}$$ and $$\{\sigma_i\}, i=1, 2, ...,$$ by Step 1: Minimizing $$\inf\limits_{\phi} L(\phi, q_{i}; \sigma_{i})$$, i.e.,   ϕi+1∈argminϕ∈V{F(ϕ)+⟨σi,Λϕ⟩+r2|Λϕ−qi|2}. Step 2: Minimizing $$\inf\limits_{q\in Z} L(\phi_{i+1}, q; \sigma_{i})$$, i.e.,   qi+1∈argminq∈Z{G(q)−⟨σi,q⟩+r2|Λϕi+1−q|2}. Step 3: Update the multiplier $$\sigma$$,   σi+1=σi+r(Λϕi+1−qi+1). For the theory of this method and its interpretation, we refer the reader to Gabay & Mercier (1976), Glowinski et al. (1981), Fortin & Glowinski (1983), Glowinski & Le Tallec (1989) and Eckstein & Bertsekas (1992). Here, we recall the convergence result of this method which is enough for our discretized problems. Theorem 4.1 (cf. Eckstein & Bertsekas, 1992, Theorem 8) Fixed $$r>0$$, assuming that $$V=\mathbb{R}^n, Z=\mathbb{R}^m$$ and that $${\it{\Lambda}}$$ has full column rank. If there exists a solution to the optimality relations (2.9), then $$\{\phi_{i}\}$$ converges to a solution of the primal problem (2.7) and $$\{\sigma_{i}\}$$ converges to a solution of the dual problem (2.8). Moreover, $$\{q_{i}\}$$ converges to $${\it{\Lambda}} \phi^{*}$$, where $$\phi^{*}$$ is the limit of $$\{\phi_{i}\}$$. The proof of this result in the case of finite-dimensional spaces $$V$$ and $$Z$$ can be found in Eckstein & Bertsekas (1992). The result holds true in infinite-dimensional Hilbert spaces under additional assumptions. One can see Fortin & Glowinski (1983) and Glowinski & Le Tallec (1989) for more details in this direction. Next, we use the ALG2 method for the discretized problems. To simplify the notations, let us drop out the subscript $$h$$ in $$(\lambda_h, u_h)$$ and $$({\it{\Phi}}_h, \theta^0_h, \theta^1_h).$$ Thanks to Remark 3.2, we treat separately the case where $$\mathbf{m}=\nu(\overline{{\it{\Omega}}})$$ and the case where $$\mathbf{m}<\nu(\overline{{\it{\Omega}}}).$$ 4.2 Partial transport ($$\mathbf{m}<\nu(\overline{{\it{\Omega}}})$$) Given $$(q_{i}, z_{i}, w_{i}), ({\it{\Phi}}_{i}, \theta^0_i, \theta^1_i)$$ at the iteration $$i,$$ we compute Step 1:   (λi+1,ui+1) ∈argmin(λ,u)∈VhFh(λ,u)+⟨(Φi,θi0,θi1),Λh(λ,u)⟩+r2|Λh(λ,u)−(qi,zi,wi)|2 =argmin(λ,u)∈Vh−⟨u,fh⟩−λ(m−ν(Ω¯))+⟨Φi,∇u⟩+⟨θi0,−u⟩+⟨θi1,u−λ⟩  +r2|∇u−qi|2+r2|u+zi|2+r2|u−λ−wi|2. Step 2:   (qi+1,zi+1,wi+1) ∈argmin(q,z,w)∈ZhGh(q,z,w)−⟨(Φi,θi0,θi1),(q,z,w)⟩+r2|Λh(λi+1,ui+1)−(q,z,w)|2 =argmin(q,z,w)∈ZhI[F∗(.,q(.))≤1](q)+I[z≤0](z)+I[w≤0](w)−⟨Φi,q⟩−⟨θi0,z⟩−⟨θi1,w⟩  +r2|∇ui+1−q|2+r2|ui+1+z|2+r2|ui+1−λi+1−w|2. Step 3: Update the multiplier   (Φi+1,θi+10,θi+11)=(Φi,θi0,θi1)+r(∇ui+1−qi+1,−ui+1−zi+1,ui+1−λi+1−wi+1). Before giving numerical results, let us take a while to comment the above iteration. Overall, Step 1 is a quadratic programming. Step 2 can be computed easily in many cases and Step 3 updates obviously. We denote by $$\text{Proj}_{C}(.)$$ the projection onto a closed convex subset $$C$$. In Step 1, we split the computation of the couple $$(\lambda_{i+1}, u_{i+1})$$ into two steps: We first minimize w.r.t. $$u$$ to compute $$u_{i+1}$$ and then we use $$u_{i+1}$$ to compute $$\lambda_{i+1}$$. More precisely, we proceed for Step 1 as follows: (1) For $$u_{i+1}$$,   ui+1 ∈argminu∈Eh−⟨u,fh⟩+⟨Φi,∇u⟩+⟨θi0,−u⟩+⟨θi1,u⟩  +r2|∇u−qi|2+r2|u+zi|2+r2|u−λi−wi|2. This is equivalent to   r⟨∇ui+1,∇v⟩+2r⟨ui+1,v⟩ =⟨fh,v⟩−⟨Φi,∇v⟩+⟨θi0,v⟩−⟨θi1,v⟩  +r⟨qi,∇v⟩−r⟨zi,v⟩+r⟨λi+wi,v⟩∀v∈Eh. Remark here that the equation is linear with a symmetric positive definite coefficient matrix. (2) For $$\lambda_{i+1}$$, it is computed explicitly   λi+1 ∈argmins∈R−s(m−ν(Ω¯))+⟨θi1,ui+1−s⟩+r2⟨ui+1−s−wi,ui+1−s−wi⟩  =−ν(Ω¯)−m−∫Ω¯θi1+r∫Ω(wi−ui+1)r∫Ω1. In Step 2, the variables $$q, z, w$$ are independent. So, we solve them separately: (1) For $$z_{i+1}$$ and $$w_{i+1}$$, if we choose $$P_2$$ finite element for $$z_{i+1}$$ and $$w_{i+1}$$, at vertex $$x_k$$,   zi+1(xk) =Proj[r∈R:r≤0](−ui+1(xk)+θi0(xk)r) =min(−ui+1(xk)+θi0(xk)r,0) (4.2) and   wi+1(xk) =Proj[r∈R:r≤0](ui+1(xk)−λi+1+θi1(xk)r) =min(ui+1(xk)−λi+1+θi1(xk)r,0). (4.3) (2) For $$q_{i+1}$$, if we choose $$P_{1}$$ finite element for $$q_{i+1}$$ then at each vertex $$x_{l}$$  qi+1(xl)=ProjBF∗(xl,.)(∇ui+1(xl)+Φi(xl)r), (4.4) where $$B_{F^*(x, .)}:=\left\{q\in \mathbb{R}^N: F^*(x, q)\le 1\right\}$$ the unit ball for $$F^*(x, .)$$. It remains to explain how we compute the projection onto $$B_{F^*(x_l, .)}$$. This issue is recently discussed in Benamou et al. (2017) for Riemann-type Finsler distances and for crystalline norms. For the convenience of the reader, we retake here the case where the unit ball of $$F(x, .)$$ is (not necessarily symmetric) convex polytope. For short, we ignore the dependence of $$x$$ in $$F$$ and $$F^*$$. Given $$d_1, ..., d_k\ne 0$$ such that, for any $$0\ne v\in \mathbb{R}^N$$, $$\max\limits_{1\le i\le k}\left\{\langle v, d_{i}\rangle\right\}> 0$$. We consider the nonsymmetric Finsler metric given by   F(v):=max1≤i≤k{⟨v,di⟩} for any v∈RN. (4.5) It is not difficult to see that the unit ball $$B^*$$ corresponding to $$F^*$$ is exactly the convex hull of $$\{d_i\}$$,   B∗=conv(di,i=1,...,k). Thus, we need to compute the projection onto the convex hull of finite points. In dimension 2, the projection onto $$B^*$$ can be performed as follows: Compute the successive vertices $$S_1, ..., S_n$$. If $$q\,{\notin}\,B^*$$ then compute the projections of $$q$$ onto the segments $$[S_{i}, S_{i+1}]$$ and compare among these projections to chose the right one. Another way is as the one in Benamou et al. (2017): Compute outward orthogonal vectors $$v_{1}, ..., v_n$$ (Fig. 1). If $$q$$ belongs to $$[S_i, S_{i+1}]+\mathbb{R}_+v_i$$ then the projection coincides with the one on the line through $$S_i, S_{i+1}$$. If $$q$$ belongs to the sector $$S_{i}+\mathbb{R}_+v_{i-1}+\mathbb{R}_+v_{i}$$, the projection is $$S_i$$. Fig. 1. View largeDownload slide Illustration of the projection. Fig. 1. View largeDownload slide Illustration of the projection. 4.3 Unbalanced transport ($$\mathbf{m}=\nu(\overline{{\it{\Omega}}})$$) Thanks to Remark 3.2, we can reduce the algorithm in this particular case by ignoring the variable $$\lambda$$. With similar considerations for $${\it{\Lambda}}_h u=(\nabla u, -u)$$, we get the following iteration: Step 1:   ui+1 ∈argminu∈Eh−⟨u,fh⟩+⟨Φi,∇u⟩+⟨θi0,−u⟩+r2|∇u−qi|2+r2|u+zi|2. Equivalently,   r⟨∇ui+1,∇v⟩+r⟨ui+1,v⟩=⟨fh,v⟩−⟨Φi,∇v⟩+⟨θi0,v⟩+r⟨qi,∇v⟩−r⟨zi,v⟩∀v∈Eh. (4.6) Step 2: (1) For $$z_{i+1}$$, choosing $$P_{2}$$ finite element for $$z_{i+1}$$ then at each vertex $$x_k$$,   zi+1(xk)=Proj[r∈R:r≤0](−ui+1(xk)+θi0(xk)r)=min(−ui+1(xk)+θi0(xk)r,0). (4.7) (2) For $$q_{i+1}$$, choosing $$P_{1}$$ finite element, at vertex $$x_l$$,   qi+1(xl)=ProjBF∗(xl,.)(∇ui+1(xl)+Φi(xl)r). Step 3: $$({\it{\Phi}}_{i+1}, \theta^0_{i+1})=({\it{\Phi}}_{i}, \theta^0_{i})+r(\nabla u_{i+1 }- q_{i+1}, -u_{i+1}-z_{i+1}).$$ 5. Numerical experiments For the numerical implementation, we use the FreeFem++ software (Hecht, 2012) and base on Benamou & Brenier (2000) and Benamou & Carlier (2015). We use $$P_{2}$$ finite element for $$u_{i}, z_{i}, w_{i}, \theta^0_{i}, \theta^1_{i}$$ and $$P_{1}$$ finite element for $${\it{\Phi}}_{i}, q_{i}$$. 5.1 Stopping criterion In computational version, the measures $$\mu$$ and $$\nu$$ are approximated by non-negative regular functions that we denote again by $$\mu$$ and $$\nu$$. We use the following stopping criteria: For the partial transport: (1) $$\text{MIN-MAX}:=\min\left\{\min\limits_{\overline{{\it{\Omega}}}} u(x), \lambda-\max\limits_{\overline{{\it{\Omega}}}} u(x), \min\limits_{\overline{{\it{\Omega}}}} \theta^0(x), \min\limits_{\overline{{\it{\Omega}}}} \theta^1(x)\right\}\!\!.$$ (2) $$\text{Max-Lip}:=\sup\limits_{\overline{{\it{\Omega}}}} F^*(x, \nabla u(x)).$$ (3) $$\text{DIV}:=\|\nabla \cdot{\it{\Phi}} + \nu -\theta^1-\mu +\theta^0\|_{L^{2}}.$$ (4) $$\text{DUAL}:=\|F(x, {\it{\Phi}}(x)) -{\it{\Phi}}(x)\cdot \nabla u\|_{L^{2}}$$. (5) $$\text{MASS}:= \left |\int (\nu -\theta^1)\,{\rm d}x - \mathbf{m}\right |$$. For the unbalanced transport: We change (1) $$\text{MIN-MAX}:=\min\left\{\min\limits_{\overline{{\it{\Omega}}}} u(x), \min\limits_{\overline{{\it{\Omega}}}} \theta^0(x)\right\}$$. (2) $$\text{DIV}:=\|\nabla \cdot {\it{\Phi}} + \nu-\mu +\theta^0\|_{L^{2}}.$$ We expect $$\text{MIN-MAX}\ge 0, \text{Max-Lip}\le 1$$; DIV, DUAL and MASS are small. 5.2 Some examples In all the examples below, we take $${\it{\Omega}}=[0, 1]\times [0, 1]$$. We test for the Riemannian case and the crystalline case. For the latter, we consider the Finsler metric of the form $$F(x, v)=\max\limits_{1\le i\le k}\left\{\langle v, d_i\rangle\right\}$$ with given directions $$d_1, ..., d_k$$ such that for any $$0\ne v\in \mathbb{R}^2$$,   max1≤i≤k{⟨v,di⟩}>0. 5.2.1 For the unbalanced transport Example 5.1 Taking $$\mu=3{\mathcal{L}^{2}}_{{\it{\Omega}}}$$ and $$\nu=\delta_{(0.5, 0.5)}$$ the Dirac mass at $$(0.5, 0.5)$$. The Finsler metric is the Euclidean one. The optimal flow is given in Fig. 2. The stopping criterion at each iteration is given in Fig. 3. Fig. 2. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, v)=|v|$$. Fig. 2. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, v)=|v|$$. Fig. 3. View largeDownload slide Stopping criterion at each iteration. Fig. 3. View largeDownload slide Stopping criterion at each iteration. Example 5.2 We take $$\mu$$ and $$\nu$$ as in the previous example and the Finsler metric given by $$F(x, v)\,{:=}|v_1|\,{+}\,|v_2|$$ for $$v\,{=}\,(v_1, v_2)\,{\in}\, \mathbb{R}^2$$. This corresponds to the crystalline norm with $$d_1\,{=}\,(1, 1), d_2\,{=}(-1, 1), d_3\,{=}\,(-1, -1) \ \text{and} \ d_4\,{=}\,(1, -1)$$. The optimal flow is given in Fig. 4 and the stopping criterion at each iteration is given in Fig. 5. Fig. 4. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, (v_1, v_2))=|v_1|+|v_2|$$. Fig. 4. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, (v_1, v_2))=|v_1|+|v_2|$$. Fig. 5. View largeDownload slide Stopping criterion at each iteration. Fig. 5. View largeDownload slide Stopping criterion at each iteration. 5.2.2 For the partial transport Example 5.3 Taking $$\mu\,{=}\,4\chi_{[(x-0.3)^2+(y-0.2)^2 <0.03]}$$, and $$\nu\,{=}\,4\chi_{[(x-0.7)^2+(y-0.8)^2<0.03]}.$$ The mass of the transport is $$\mathbf{m}\,{:=}\,\frac{\nu(\overline{{\it{\Omega}}})}{2}.$$ We test for different Finsler metrics. On each figure below, the subfigure at left illustrates the unit ball of $$F$$ and the subfigure at right gives the numerical result (see Figs 6–9). The stopping criteria are summarized in Table 1. Table 1 Stopping criteria for $$800$$ iterations Case  DIV  DUAL  MASS  MIN–MAX  Max–Lip  Time execution (s)  1  2.48182e-05  9.5294e-06  0.000161361  $$-$$0.0149942  1.00068  357  2  3.38395e-05  5.58717e-05  0.000195881  $$-$$0.00120123  1.00248  867  3  7.44768e-05  5.5997e-05  6.66404e-06  $$-$$0.00272389  1.00351  1269  4  6.33726e-05  3.20691e-05  0.000120909  $$-$$0.0104915  1.02572  1123  Case  DIV  DUAL  MASS  MIN–MAX  Max–Lip  Time execution (s)  1  2.48182e-05  9.5294e-06  0.000161361  $$-$$0.0149942  1.00068  357  2  3.38395e-05  5.58717e-05  0.000195881  $$-$$0.00120123  1.00248  867  3  7.44768e-05  5.5997e-05  6.66404e-06  $$-$$0.00272389  1.00351  1269  4  6.33726e-05  3.20691e-05  0.000120909  $$-$$0.0104915  1.02572  1123  Fig. 6. View largeDownload slide Case 1: $$F(x, v)=|v|$$. Fig. 6. View largeDownload slide Case 1: $$F(x, v)=|v|$$. Fig. 7. View largeDownload slide Case 2: The crystalline case with $$d_1=(1, 1), d_2=(-1, 1), d_3=(-1, -1)$$ and $$d_4=(1, -1)$$. Fig. 7. View largeDownload slide Case 2: The crystalline case with $$d_1=(1, 1), d_2=(-1, 1), d_3=(-1, -1)$$ and $$d_4=(1, -1)$$. Fig. 8. View largeDownload slide Case 3: The crystalline case with $$d_1\,{=}\,(1, 0), d_2\,{=}\,(\frac{1}{5}, \frac{1}{5}), d_3\,{=}\,(-\frac{1}{5}, \frac{1}{5}), d_4\,{=}\,(-\frac{1}{5}, -\frac{1}{5})$$ and $$d_5\,{=}\,(\frac{1}{5}, -\frac{1}{5})$$ makes the transport more expensive in the direction of the vector $$(1, 0)$$. Fig. 8. View largeDownload slide Case 3: The crystalline case with $$d_1\,{=}\,(1, 0), d_2\,{=}\,(\frac{1}{5}, \frac{1}{5}), d_3\,{=}\,(-\frac{1}{5}, \frac{1}{5}), d_4\,{=}\,(-\frac{1}{5}, -\frac{1}{5})$$ and $$d_5\,{=}\,(\frac{1}{5}, -\frac{1}{5})$$ makes the transport more expensive in the direction of the vector $$(1, 0)$$. Fig. 9. View largeDownload slide Case 4: The crystalline case with $$d_1\,{=}\,(1, -1), d_2\,{=}\,(1, -\frac{4}{5}), d_3\,{=}\,(-\frac{4}{5}, 1), d_4\,{=}\,(-1, 1)$$ and $$d_5\,{=}\,(-1, -1)$$ makes the transport cheaper in the direction of the vector $$(1, 1)$$. Fig. 9. View largeDownload slide Case 4: The crystalline case with $$d_1\,{=}\,(1, -1), d_2\,{=}\,(1, -\frac{4}{5}), d_3\,{=}\,(-\frac{4}{5}, 1), d_4\,{=}\,(-1, 1)$$ and $$d_5\,{=}\,(-1, -1)$$ makes the transport cheaper in the direction of the vector $$(1, 1)$$. Example 5.4 Let $$\mu\,{=}\,2\chi_{[(x-0.2)^2\,{+}\,(y-0.2)^2 <0.03]}+2\chi_{[(x-0.6)^2+(y-0.1)^2 <0.01]}$$ and $$\nu\,{=}\,2\chi_{[(x-0.6)^2+(y-0.8)^2 <0.03]}$$. In this example, we take the Euclidean norm and we let $$\mathbf{m}$$ vary by taking the values $$\mathbf{m}_{i}\,{=}\,\frac{i}{6}$$$$\min\{\mu({\it{\Omega}}), \nu({\it{\Omega}})\},$$$$i=1, ..., 6.$$ The results are given in Fig. 10. Fig. 10. View largeDownload slide Optimal flows. Fig. 10. View largeDownload slide Optimal flows. Acknowledgements The authors are grateful to J. D. Benamou and G. Carlier who provide some codes of ALG2 on the link https://team.inria.fr/mokaplan/software/. Some parts of our codes are inspired from their work. Appendix Our aim here is to show Lemma A.1 that gives a smooth approximation of 1-$$d_{F}$$ Lipschitz continuous function for continuous nondegenerate Finsler metrics $$F$$. This result is more or less known in some particular cases. However, we could not find any rigorous proofs for the general case in the literature. Lemma A.1 Let $${\it{\Omega}}$$ be a connected bounded Lipschitz domain and $$F$$ be a continuous nondegenerate Finsler metric on $$\overline{{\it{\Omega}}}$$. For any Lipschitz continuous function $$u$$ on $$\overline{{\it{\Omega}}}$$ satisfying   F∗(x,∇u(x))≤1 a.e. x∈Ω, (A.1) there exists a sequence of functions $$u_{\varepsilon}\in C^{\infty}_{c}(\mathbb{R}^N)$$ such that   F∗(x,∇uε(x))≤1∀x∈Ω¯, and   uε⇉u uniformly on Ω¯. Note that $$F$$ and $$F^*$$ are defined only in $$\overline{{\it{\Omega}}}$$ and that the gradient of $$u$$ is controlled only inside of $${\it{\Omega}}$$ by (A.1). If we use the standard convolution to define $$u_{\varepsilon}$$, the value of $$u_{\varepsilon}(x)$$ is affected by the value of $$u(y)$$ outside of $$\overline{{\it{\Omega}}}$$ which remains uncontrolled. To overcome this difficulty, if $$x$$ is near the boundary, we move it a little into inside of $${\it{\Omega}}$$ before taking the convolution. To this aim, we use the smooth partition of unity tool to deal with approximation of $$u$$ near the boundary. Proof. Set   ∀x∈RN,u~(x):={u(x) if x∈Ω¯0 otherwise.  Step 1: Fix $$z\in \partial {\it{\Omega}}$$. Since $${\it{\Omega}}$$ is a Lipschitz domain, there exist $$r_{z}>0$$ and a Lipschitz continuous function $$\gamma_{z}:\mathbb{R}^{N-1} \longrightarrow \mathbb{R}$$ such that (up to rotating and relabeling if necessary)   Ω∩B(z,rz)={x|xN>γz(x1,...,xN−1)}∩B(z,rz). Set $$U_{z}:={\it{\Omega}} \cap B(z, \frac{r_{z}}{2})$$. For any $$x\in \mathbb{R}^N$$, taking   xzε:=x+ελzen, (A.2) where we choose a sufficiently large fixed $$\lambda_{z}$$ and all small $$\varepsilon$$, say fixed $$\lambda_{z} \ge \text{Lip} (\gamma_{z})+1$$, $$0<\varepsilon< \frac{r_{z}}{2(\lambda_{z}+1)}.$$ By this choice and the Lipschitz property of $$\gamma_{z}$$, we see that   B(xzε,ε)⊂Ω∩B(z,rz) for all x∈Uz. (A.3) Defining   u~ε(x):=∫RNρε(y)u~(xzε−y)dy=∫B(xzε,ε)ρε(xzε−y)u~(y)dy for all x∈RN, (A.4) where $$\rho_{\varepsilon}$$ is the standard mollifier on $$\mathbb{R}^N$$. Obviously, $$\tilde{u}_{\varepsilon} \in C^{\infty}_{c}(\mathbb{R}^N)$$. Using (A.3), (A.4) and the continuity of $$u$$ on $$\overline{{\it{\Omega}}}$$, we get   u~ε⇉u on U¯z. Step 2: Now, using the compactness of $$\partial {\it{\Omega}}$$ and $$\partial {\it{\Omega}} \,{\subset}\, \bigcup\limits_{z\,{\in}\, \partial {\it{\Omega}}} B(z, \frac{r_{z}}{2})$$, there exist $$z_{1}, ..., z_{n}\in \partial {\it{\Omega}}$$ such that   ∂Ω⊂⋃i=1nB(zi,rzi2). For short, we write $$r_{i}, U_{i}, x_{i}$$ instead of $$r_{z_{i}}, U_{z_{i}}, x_{z_{i}}$$. Taking an open set $$U_{0}\Subset {\it{\Omega}}$$ such that   Ω¯⊂⋃i=1nB(zi,ri2)⋃U0. Let $$\{\phi\}^{n}_{i=0}$$ be a smooth partition of unity on $$\overline{{\it{\Omega}}}$$, subordinate to $$\left\{U_{0}, B(z_{1}, \frac{r_{1}}{2}), ..., B(z_{n}, \frac{r_{n}}{2})\right\}$$, that is,   {ϕi∈Cc∞(RN)0≤ϕi≤1,∀i=0,...,nsupp(ϕi)⋐B(zi,ri2)∀i=1,...,n, supp(ϕ0)⋐U0∑i=0nϕi(x)=1 for all x∈Ω¯.  Because of Step 1, there exist $$\tilde{u}^1_\varepsilon, ..., \tilde{u}^n_\varepsilon \in C^{\infty}_{c}(\mathbb{R}^N)$$ such that   u~εi⇉u on U¯ii=1,...,n. For $$i=0,$$ since $$U_{0}\Subset {\it{\Omega}}$$, we can take $$\tilde{u}^0_\varepsilon:= \rho_{\varepsilon}\star \tilde{u} \in C^{\infty}_{c}(\mathbb{R}^N)$$ and $$\tilde{u}^0_{\varepsilon} \rightrightarrows u$$ on $$\overline{U}_{0}$$. Set   uε:=11+Cε+w(ε)∑i=0nϕiu~εi, where $$C$$ is chosen later and   w(ε):=sup{|F∗(x,p)−F∗(y,p)|:x,y∈Ω¯,|x−y|≤Mε,|p|≤‖∇u‖L∞} with constant $$M:=\max\limits_{1\le i\le n}\{\lambda_{z_{i}}+1\}$$, $$\lambda_{z_{i}}$$ is given in Step 1. We show that $$u_{\varepsilon}$$ satisfies all the desired properties. By the construction, $$u_{\varepsilon} \in C^{\infty}_{c}(\mathbb{R}^N)$$ and   uε⇉∑i=0nϕiu=u on Ω¯. At last, we show that $$F^*(x, \nabla u_{\varepsilon} (x))\le 1 \quad\forall x\in \overline{{\it{\Omega}}}$$. Indeed, for any $$x\in {\it{\Omega}}$$, if $$x\in U_{i}, i=1, ..., n$$ (near the boundary of $${\it{\Omega}}$$), we move $$x$$ a bit into inside of $${\it{\Omega}}$$ to $$x^\varepsilon_i:=x^\varepsilon_{z_i}$$ (see (A.2) and (A.3)), if $$x\in U_{0}$$, set $$x_0^{\varepsilon}=x$$. We have   ∇uε(x) =11+Cε+w(ε)(∑i=0n∇ϕi(x)u~εi(x)+∑i=0nϕi(x)∇u~εi(x)) =11+Cε+w(ε)(∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy  +∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)∇u(y)dy). The first sum on the right-hand side has a small norm. Indeed, using the fact that   ∑i=0n∇ϕi(x)u(x)=0 for all x∈Ω, we have   ∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy=∑i=0n∇ϕi(x)(∫B(xiε,ε)ρε(xiε−y)u(y)dy−u(x)). (A.5) Moreover,   |∫B(xiε,ε)ρε(xiε−y)u(u)dy−u(x)| ≤|∫B(xiε,ε)ρε(xiε−y)(u(y)−u(xiε))dy|+|u(xiε)−u(x)| ≤C1ε∀i=0,...,n, where the constant $$C_{1}$$ depends only on Lip($$\gamma_{z_{i}}$$) and the Lipschitz constant of $$u$$ on $$\overline{{\it{\Omega}}}$$. Thus, by combining this with (A.5),   |∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy|≤C2ε∀x∈Ω, where $$C_{2}$$ depends only on $$C_{1}$$ and $$\|\nabla \phi_{i}\|_{L^{\infty}}$$. Using the nondegeneracy of $$F$$, we have   F∗(x,∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy)≤C3ε for all x∈Ω. Fixed any $$x\in {\it{\Omega}}$$, if $$y\in B(x_i^{\varepsilon}, \varepsilon)$$ then $$|x-y|\le |x-x_i^{\varepsilon}|+|x_i^{\varepsilon}-y|\le M\varepsilon$$. So we obtain   F∗(x,∇uε(x)) ≤11+Cε+w(ε)[F∗(x,∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy)  +F∗(x,∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)∇u(y)dy)] ≤11+Cε+w(ε)(C3ε+∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)F∗(x,∇u(y))dy) ≤11+Cε+w(ε)[C3ε+∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)F∗(y,∇u(y))dy  +∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)(F∗(x,∇u(y))−F∗(y,∇u(y)))dy] ≤C3ε+1+w(ε)1+Cε+w(ε) ≤1(choose a constant C≥C3). By the continuity of $$\nabla u_{\varepsilon}$$ and of $$F^*$$, we also have $$F^*(x, \nabla u_{\varepsilon}(x))\le 1 \, \quad\forall x\in \overline{{\it{\Omega}}}.$$ □ Proposition A.2 Let $$F$$ be a continuous nondegenerate Finsler metric on a connected bounded Lipshitz domain $${\it{\Omega}}$$. We have   LipdF={u:Ω¯⟶R|u is Lipschitz continuous and F∗(x,∇u(x))≤1 a.e. x∈Ω}:=BF∗. (A.6) As a consequence, for any 1-$$d_{F}$$ Lipschitz continuous function $$u$$, there exists a sequence of 1-$$d_{F}$$ Lipschitz continuous functions $$u_{\varepsilon}\in C^{\infty}_{c}(\mathbb{R}^N)$$ and $$u_{\varepsilon} \rightrightarrows u$$ uniformly on $$\overline{{\it{\Omega}}}$$. Lemma A.3 We have $$Lip_{d_{F}}\subset \mathcal{B}_{F^*}.$$ Proof. Let $$u\in Lip_{d_{F}}$$. Then $$u$$ is Lipschitz and $$u$$ is differentiable a.e. in $${\it{\Omega}}$$. Let $$x\in {\it{\Omega}}$$ be any point where $$u$$ is differentiable. We have, for any $$v\in \mathbb{R}^N$$,   ⟨∇u(x),v⟩F(x,v) =limh→0u(x+hv)−u(x)F(x,hv) ≤lim suph→0dF(x,x+hv)F(x,hv) ≤lim suph→0∫01F(x+thv,hv)dtF(x,hv)=1. Hence, $$F^*(x, \nabla u(x))\le 1$$. So $$u\in \mathcal{B}_{F^*}$$. □ Lemma A.4 We have $$\mathcal{B}_{F^*}\subset Lip_{d_{F}}.$$ Proof. Fix any $$u\in \mathcal{B}_{F^*}$$. Case 1: If $$u$$ is smooth then $$F^*(x, \nabla u(x))\le 1 \quad\forall x\in \overline{{\it{\Omega}}}$$. For any $$x, y \in \overline{{\it{\Omega}}}$$ and any Lipschitz curve $$\xi$$ in $$\overline{{\it{\Omega}}}$$ joining $$x$$ and $$y$$, we have   u(y)−u(x) =∫01∇u(ξ(t))ξ˙(t)dt ≤∫01F∗(ξ(t),∇u(ξ(t)))F(ξ(t),ξ˙(t))dt ≤∫01F(ξ(t),ξ˙(t))dt. Hence $$u\in Lip_{d_{F}}$$. Case 2: For general Lipschitz continuous function $$u$$ satisfying $$F^*(x, \nabla u(x))\,{\le}\, 1 \, \text{ a.e. }\, x\in {\it{\Omega}}$$, thanks to Lemma A.1, there exist $$u_{\varepsilon} \in \mathcal{B}_{F^*}\bigcap C^{\infty}_{c}(\mathbb{R}^N)$$ such that $$u_{\varepsilon} \,{\rightrightarrows}\, u$$ on $$\overline{{\it{\Omega}}}$$. According to Case 1 above, $$u_{\varepsilon}\in Lip_{d_{F}}$$. Since $$u_{\varepsilon} \rightrightarrows u$$ on $$\overline{{\it{\Omega}}}$$, we obtain $$u\in Lip_{d_{F}}$$. □ Proof of Proposition A.2. The proof follows by Lemma A.3 and Lemma A.4. □ Proof of Lemma 2.4. Since $$0\le u\le \lambda$$, the sequence $$u_{\varepsilon}$$ in the proof of Lemma A.1 satisfies $$0\le u_{\varepsilon}\le \lambda$$. So $$u_{\varepsilon}\in C^{\infty}_{c}(\mathbb{R}^N) \cap L^\lambda_{d_{F}}$$ and $$u_{\varepsilon} \rightrightarrows u$$ on $$\overline{{\it{\Omega}}}$$. □ Remark A.5 The results still hold true if $${\it{\Omega}}$$ is connected, bounded and has the segment property. References Ambrosio, L., Fusco, N. & Pallara, D. ( 2000) Functions of Bounded Variation and Free Discontinuity Problems.  Oxford Mathematical Monographs. Oxford: Oxford University Press. Ambrosio, L., Gigli, N. & Savaré, G. ( 2005) Gradient Flows in Metric Spaces and in the Space of Probability Measures. Lectures in Mathematics.  Basel: Birkháuser. Barrett, J. W. & Prigozhin, L. ( 2009) Partial L1 Monge-Kantorovich problem: Variational formulation and numerical approximation. Interface. Free Bound.,  11, 201– 238. Google Scholar CrossRef Search ADS   Beckmann, M. ( 1952) A continuous model of transportation. Econometrica,  20, 643– 660. Google Scholar CrossRef Search ADS   Benamou, J. D. & Brenier, Y. ( 2000) A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math.,  84, 375– 393. Google Scholar CrossRef Search ADS   Benamou, J. D. & Carlier, G. ( 2015) Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations. J. Optim. Theory Appl. , 167, 1– 26. Google Scholar CrossRef Search ADS   Benamou, J. D., Carlier, G., Cuturi, M., Nenna, L. & Peyré, G. ( 2015) Iterative Bregman projections for regularized transportation problems. SIAM J. Sci. Comput.,  37, A1111– A1138. Google Scholar CrossRef Search ADS   Benamou, J. D., Carlier, G. & Hatchi, R. ( 2017) A numerical solution to Monge’s problem with a Finsler distance cost. ESAIM: M2AN,  DOI: http://dx.doi.org/10.1051/m2an/2016077. Caffarelli, L. & McCann, R. J. ( 2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Ann. Math.,  171, 673– 730. Google Scholar CrossRef Search ADS   Chizat, L., Peyré, G., Schmitzer, B. & Vialard, F. X. ( 2016) Scaling algorithms for unbalanced transport problems. arXiv preprint arXiv:1607.05816. De Pascale, L., Evans, L. C. & Pratelli, A. ( 2004) Integral estimates for transport densities. Bull. London Math. Soc.,  36, 383– 395. Google Scholar CrossRef Search ADS   De Pascale, L. & Pratelli, A. ( 2004) Sharp summability for Monge transport density via interpolation. ESAIM Contr. Optim. Ca. Va. , 10, 549– 552. Google Scholar CrossRef Search ADS   Eckstein, J. & Bertsekas, D. P. ( 1992) On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators. Math. Program.,  55, 293– 318. Google Scholar CrossRef Search ADS   Ekeland, I. & Teman, R. ( 1976) Convex Analysis and Variational Problems.  Amsterdam-New York: North-Holland American Elsevier. Feldman, M. & McCann, R. J. ( 2002) Uniqueness and transport density in Monge’s mass transportation problem. Calc. Var. Partial Differ. Equ. , 15, 81– 113. Google Scholar CrossRef Search ADS   Figalli, A. ( 2010) The Optimal Partial Transport Problem. Arch. Rational Mech. Anal. , 195, 533– 560. Google Scholar CrossRef Search ADS   Fortin, M. & Glowinski, R. ( 1983) Augmented Lagrangian methods: Applications to the Numerical Solution of Boundary-Value Problems , vol 15. Amsterdam: North-Holland Publishing. Gabay, D. & Mercier, B. ( 1976) A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Comput. Math. Appl. , 2, 17– 40. Google Scholar CrossRef Search ADS   Glowinski, R., Lions, J. L. & Tre1molie1res, R. ( 1981) Numerical Analysis of Variational Inequalities.  Amsterdam: North-Holland Publishing. Glowinski, R. & Le Tallec, P. ( 1989) Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics,  vol. 9. Philadelphia: SIAM. Google Scholar CrossRef Search ADS   Hecht, F. ( 2012) New development in freefem++. J. Numer. Math.,  20, 251– 266. Google Scholar CrossRef Search ADS   Igbida, N. & Nguyen, V. T. ( 2017) Optimal partial mass transportation and obstacle Monge-Kantorovich equation (submitted for publication). Igbida, N. & Ta Thi, N. N. ( 2017) Sub-gradient Diffusion Operator. J. Differ. Equ . 262, 3837– 3863. Google Scholar CrossRef Search ADS   Santambrogio, F. ( 2009) Absolute continuity and summability of transport densities: simpler proofs and new estimates. Calc. Var. Partial Differ. Equ.,  36, 343– 354. Google Scholar CrossRef Search ADS   Santambrogio, F. ( 2015) Optimal Transport for Applied Mathematicians . Basel: Birkhäuser. Google Scholar CrossRef Search ADS   Villani, C. ( 2003) Topics in Optimal Transportation . Graduate Studies in Mathematics, vol. 58. Providence: American Mathematical Society. Villani, C. ( 2009) Optimal Transport, Old and New . Grundlehren des Mathematischen Wissenschaften (Fundamental Principles of Mathematical Sciences), vol. 338. Berlin, Heidelberg: Springer. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# Augmented Lagrangian Method for Optimal Partial Transportation

, Volume 38 (1) – Jan 1, 2018
28 pages

/lp/ou_press/augmented-lagrangian-method-for-optimal-partial-transportation-VgteHYjKL6
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drw077
Publisher site
See Article on Publisher Site

### Abstract

Abstract The use of augmented Lagrangian algorithm for optimal transport problems goes back to Benamou & Brenier (2000, Acomputational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math., 84, 375–393), in the case where the cost corresponds to the square of the Euclidean distance. It was recently extended in Benamou & Carlier (2015, Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations. J. Optim. Theory Appl., 167, 1–26), to the optimal transport with the Euclidean distance and Mean-Field Games theory and in Benamou et al. (2017, A numerical solution to Monge’s problem with a Finsler distance cost ESAIM: M2AN), to the optimal transportation with Finsler distances. Our aim here is to show how one can use this method to study the optimal partial transport problem with Finsler distance costs. To this aim, we introduce a suitable dual formulation of the optimal partial transport, which contains all the information on the active regions and the associated flow. Then, we use a finite element discretization with the FreeFem++ software to provide numerical simulations for the optimal partial transportation. A convergence study for the potential together with the flux and the active regions is given to validate the approach 1. Introduction The theory of optimal transportation deals with the problem to find the optimal way to move materials from a given source to a desired target in such a way to minimize the work. The problem was first proposed and studied by G. Monge in 1781, and then L. Kantorovich made fundamental contributions to the problem in the 1940s by relaxing the problem into a linear one. Since the late 1980s, this subject has been investigated under various points of view with many applications in image processing, geometry, probability theory, economics, evolution partial differential equations (PDEs) and other areas. For more information on the optimal mass transport problem, we refer the reader to the pedagogical books (Villani, 2003; Ambrosio et al., 2005; Villani, 2009; Santambrogio, 2015). The standard optimal transport problem requires that the total mass of the source is equal to the total mass of the target (balance condition of mass) and that all the materials of the source must be transported. Here, we are interested in the optimal partial transportation. That is the case where the balance condition of mass is excluded, and the aim is to transport effectively a prescribed amount of mass from the source to the target. In other words, the optimal partial transport problem aims to study the practical situation, where only a part of the commodity (respectively, consumer demand) of a prescribed total mass $$\mathbf{m}$$ needs to be transported (respectively, fulfilled). This generalized problem brings out additional variables called active regions. The problem was first studied theoretically in Caffarelli & McCann (2010) (see also Figalli, 2010) in the case where the work is proportional to the square of the Euclidean distance. Recently in Igbida & Nguyen (2017), we give a complete theoretical study of the problem in the case where the work is proportional to a Finsler distance $$d_{F}$$ (covering by the way the case of the Euclidean distance), where $$d_{F}$$ is given as follows (see Section 2)   dF(x,y):=infξ∈Lip([0,1];Ω¯){∫01F(ξ(t),ξ˙(t))dt:ξ(0)=x,ξ(1)=y}. Concerning numerical approximations for the optimal partial transport, Barrett & Prigozhin (2009) studied the case of the Euclidean distance by using an approximation based on nonlinear approximated PDEs and Raviart–Thomas finite elements. Benamou et al. (2015) and Chizat et al. (2016) introduced general numerical frameworks to approximate solutions to linear programs related to the optimal transport (including the optimal partial transport). Their idea is based on an entropic regularization of the initial linear programs. This is a static approach to optimal transport-type problems and needs to use (approximated) values of $$d_{F}(x, y)$$. In this article, we use a different approach (based mainly on Benamou & Brenier, 2000; Benamou & Carlier, 2015; Igbida & Nguyen, 2017) to compute the solution of the optimal partial transport problem. We first show how one can directly reformulate the unknown quantities (variables) of the optimal partial transport into an infinite-dimensional minimization problem of the form:   minϕ∈VF(ϕ)+G(Λϕ), where $$\mathcal{F}, \mathcal{G}$$ are l.s.c., convex functionals and $${\it{\Lambda}}\in \mathcal{L}(V, Z)$$ is a continuous linear operator between two Banach spaces. Thanks to peculiar properties of $$\mathcal{F}$$ and $$\mathcal{G}$$ in our situation, an augmented Lagrangian method is applied effectively in the same spirit of Benamou & Carlier (2015) and Benamou et al. (2017). We show that, for the computation, we just need to solve linear equations (with a symmetric positive definite coefficient matrix) or to update explicit formulations. It is worth noting that this method uses only elementary operations without evaluating $$d_{F}$$. This article is organized as follows: In the next section, we introduce the optimal partial transport problem and its equivalent formulations with a particular attention to the Kantorovich dual formulation. In Section 3, we give a finite-dimensional approximation of the problem, and show that primal-dual solutions of the discretized problems converge to the ones of original continuous problems. The details of the ALG2 algorithm is given in Section 4. Some numerical examples are presented in Section 5. We terminate the article by an appendix, where we give proofs of some facts we need in this article. 2. Partial transport and its equivalent formulations Let $${\it{\Omega}}$$ be a connected bounded Lipschitz domain and $$F$$ be a continuous Finsler metric on $$\overline{{\it{\Omega}}}$$, i.e., $$F: \overline{{\it{\Omega}}}\times \mathbb{R}^N \longrightarrow [0, +\infty)$$ is continuous and $$F(x, .)$$ is convex, positively homogeneous of degree $$1$$ in the sense   F(x,tv)=tF(x,v)∀t>0,v∈RN. We assume moreover that $$F$$ is nondegenerate in the sense that there exist positive constants $$M_{1}, M_{2}$$, such that   M1|v|≤F(x,v)≤M2|v|∀x∈Ω¯,v∈RN. Let $$\mu, \nu\in \mathcal M^+_{b}(\overline{{\it{\Omega}}})$$ be two Radon measures on $$\overline{{\it{\Omega}}}$$ and $$\mathbf{m}_{\max}:=\min\{\mu(\overline{{\it{\Omega}}}), \nu(\overline{{\it{\Omega}}})\}.$$ Given a total mass $$\mathbf{m}\in [0, \mathbf{m}_{\max}]$$, the optimal partial transport problem (or partial Monge–Kantorovich problem, PMK for short) aims to transport effectively the total mass $$\mathbf{m}$$ from a supply subregion of the source $$\mu$$ into a subregion of the target $$\nu$$. The set of subregions of mass $$\mathbf{m}$$ is given by   Subm(μ,ν):={(ρ0,ρ1)∈Mb+(Ω¯)×Mb+(Ω¯):ρ0≤μ,ρ1≤ν,ρ0(Ω¯)=ρ1(Ω¯)=m}. An element $$(\rho_{0}, \rho_{1})\in Sub_{\mathbf{m}}(\mu, \nu)$$ is called a couple of active regions. As for the optimal transport, one can work with different kinds of cost functions for the optimal partial transport, i.e., in the formulation (2.1) below, $$d_{F}(x, y)$$ can be replaced by a general measurable cost function $$c(x, y)$$. However, in this article, we focus on the case where the cost $$c=d_{F}$$. So let us state the problem directly for $$d_{F}$$. The PMK problem (Barrett & Prigozhin, 2009; Caffarelli & McCann, 2010; Figalli, 2010; Igbida & Nguyen, 2017) aims to minimize the following problem   min{K(γ):=∫Ω¯×Ω¯dF(x,y)dγ:γ∈πm(μ,ν)}, (2.1) where $$d_{F}$$ is the Finsler distance on $$\overline{{\it{\Omega}}}$$ associated with $$F$$, i.e.,   dF(x,y):=inf{∫01F(ξ(t),ξ˙(t))dt:ξ(0)=x,ξ(1)=y,ξ∈Lip([0,1];Ω¯)}. $${\large{\mathbf{\pi}}}_{\mathbf{m}}(\mu, \nu)$$ is the set of transport plans of mass $$\mathbf{m}$$, i.e.,   πm(μ,ν):={γ∈Mb+(Ω¯×Ω¯):(πx#γ,πy#γ)∈Subm(μ,ν)}. Here, $$\pi_x \#\gamma$$ and $$\: \pi_y \# \gamma$$ are the first and second marginals of $$\gamma$$. An optimal $$\gamma^*$$ is called an optimal plan and $$(\pi_x \#\gamma^*, \pi_y \# \gamma^*)$$ is called a couple of optimal active regions. Following Igbida & Nguyen (2017), to study the PMK problem we use its dual problem that we call the dual partial Monge–Kantorovich problem (DPMK). To this aim, we consider $$Lip_{d_{F}}$$ the set of $$1-$$Lipschitz continuous functions w.r.t. $$d_{F}$$ given by   LipdF:={u:Ω¯⟶R|u(y)−u(x)≤dF(x,y)∀x,y∈Ω¯}. Then, the connection between the PMK and DPMK problems is summarized in the following theorem. Theorem 2.1 Let $$\mu, \nu\,{\in}\,\mathcal M^+_b(\overline{{\it{\Omega}}})$$ be Radon measures and $$\mathbf{m}\,{\in}\, [0,\mathbf{m}_{\max}]$$. The partial Monge–Kantorovich problem has a solution $$\sigma^*\,{\in}\, {\large{\mathbf{\pi}}}_{\mathbf{m}}(\mu, \nu)$$ and   K(σ∗)=max{D(λ,u):=∫Ω¯ud(ν−μ)+λ(m−ν(Ω¯)):λ≥0 and u∈LdFλ}, (2.2) where   LdFλ:={u∈LipdF:0≤u(x)≤λ for any x∈Ω¯}. Moreover, $$\sigma\in {\large{\mathbf{\pi}}}_{\mathbf{m}}(\mu, \nu)$$ and $$(\lambda, u)\in \mathbb{R}^+\times L^\lambda_{d_{F}}$$ are solutions, respectively, if and only if   u(x)=0 for (μ−πx#σ)-a.e. x∈Ω¯,u(x)=λ for (ν−πy#σ)-a.e. x∈Ω¯and u(y)−u(x)=dF(x,y) for σ-a.e. (x,y)∈Ω¯×Ω¯. Proof. The proof follows in the same way of Theorem 2.4 in Igbida & Nguyen (2017), where the authors study the case $${\it{\Omega}}=\mathbb{R}^N.$$ □ The DPMK problem (2.2) contains all the information concerning the optimal partial mass transportation. However, for the numerical approximation of the optimal partial transportation and to use the augmented Lagrangian method, we need to rewrite the problem into the form   infϕ∈VF(ϕ)+G(Λϕ). To do that, we consider the polar function $$F^*$$ of $$F,$$ which is defined by   F∗(x,p):=sup{⟨v,p⟩:F(x,v)≤1} for x∈Ω¯,p∈RN. Note that $$F^*(x, .)$$ is not the Legendre–Fenchel transform. It is easy to see that $$F^*$$ is also a continuous, nondegenerate Finsler metric on $$\overline{{\it{\Omega}}}$$ and   ⟨v,p⟩≤F∗(x,p)F(x,v)∀x∈Ω¯,v,p∈RN. Remark 2.2 Using the polar function $$F^*$$, we can characterize the set $$Lip_{d_{F}}$$ as (see the appendix if necessary)   LipdF={u:Ω¯⟶R|u is Lipschitz continuous and F∗(x,∇u(x))≤1 a.e. x∈Ω}. Thanks to this remark, the DPMK problem (2.2) can be written as   max{D(λ,u):0≤u(x)≤λ,u is Lipschitz continuous, F∗(x,∇u(x))≤1 a.e. x∈Ω}. Moreover, we have Theorem 2.3 Under the assumptions of Theorem 2.1, setting $$V:=\mathbb{R}\times C^1(\overline{{\it{\Omega}}})$$ and $$Z:=C(\overline{{\it{\Omega}}})^{N}\times C(\overline{{\it{\Omega}}})\times C(\overline{{\it{\Omega}}}),$$ we have   K(σ∗)=−inf{F(λ,u)+G(Λ(λ,u)):(λ,u)∈V}, (2.3) where $${\it{\Lambda}}\in \mathcal{L}(V, Z)$$ is given by   Λ(λ,u):=(∇u,−u,u−λ)∀(λ,u)∈V, and $$\mathcal{F}: V\longrightarrow (-\infty, +\infty]$$, $$\mathcal{G}: Z\longrightarrow (-\infty, +\infty]$$ are the l.s.c. convex functions given by   F(λ,u) :=−∫Ω¯ud(ν−μ)−λ(m−ν(Ω¯))∀(λ,u)∈V;G(q,z,w) :={0if z(x)≤0,w(x)≤0,F∗(x,q(x))≤1∀x∈Ω¯+∞otherwise  for (q,z,w)∈Z. To prove this theorem, we need the following lemma. Lemma 2.4 Let $$\lambda\,{\ge}\, 0$$ be fixed. For any $$u\,{\in}\, L^\lambda_{d_{F}}$$, there exists a sequence of smooth functions $$u_{\varepsilon}\in C^\infty_{c}(\mathbb{R}^N) \bigcap L^\lambda_{d_{F}}$$ such that $$u_{\varepsilon} \rightrightarrows u$$ uniformly on $$\overline{{\it{\Omega}}}$$. The result of the lemma is more or less known in some cases (see Igbida & Ta Thi, 2017 for the case where the function $$u$$ is null on the boundary). The proof in the general case is quite technical and will be given in the appendix. Proof of Theorem 2.3. Thanks to Remark 2.2 and Lemma 2.4, we have   −inf(λ,u)∈VF(λ,u)+G(Λ(λ,u)) =sup{∫Ω¯ud(ν−μ)+λ(m−ν(Ω¯)):λ≥0,u∈C1(Ω¯)∩LdFλ} =max{D(λ,u):λ≥0 and u∈LdFλ}. Using the duality (2.2), the proof is completed. □ To end this section, we prove the following result that will be useful for the proof of the convergence of our discretization. Theorem 2.5 Under the assumptions of Theorem 2.1, we have   −inf(λ,u)∈VF(λ,u)+G(Λ(λ,u))=min{∫Ω¯F(x,Φ|Φ|(x))d|Φ|:(Φ,θ0,θ1)∈Ψm(μ,ν)}, (2.4) where   Ψm(μ,ν):={(Φ,θ0,θ1)∈Z∗=Mb(Ω¯)N×Mb(Ω¯)×Mb(Ω¯):θ0≥0,θ1≥0,θ1(Ω¯)=ν(Ω¯)−m  and −∇⋅Φ=ν−θ1−(μ−θ0) with Φ.n=0 on ∂Ω}. Actually, the minimal flow-type formulation   min{∫Ω¯F(x,Φ|Φ|(x))d|Φ|:(Φ,θ0,θ1)∈Ψm(μ,ν)} (2.5) introduces the Beckmann problem (see Beckmann, 1952) for the optimal partial transport with Finsler distance costs. See here that in the balanced case, i.e., $$\mathbf{m}=\mu(\overline{{\it{\Omega}}})=\nu(\overline{{\it{\Omega}}})$$, the formulation (2.5) becomes   min{∫Ω¯F(x,Φ|Φ|(x))d|Φ|:Φ∈Mb(Ω¯)N,−∇⋅Φ=ν−μ with Φ.n=0 on ∂Ω}. (2.6) An optimal solution $${\it{\Phi}}$$ of the problem (2.6) is called an optimal flow of transporting $$\mu$$ onto $$\nu$$. As known from the optimal transport theory, the optimal flow gives a way to visualize the transportation. To prove Theorem 2.5, we will use the well-known duality arguments. For convenience, let us recall here the Fenchel–Rockafellar duality. Let us consider the problem   infϕ∈VF(ϕ)+G(Λϕ), (2.7) where $$\mathcal{F}: V\,{\longrightarrow}\, (-\infty, +\infty]$$ and $$\mathcal{G}: Z\,{\longrightarrow}\, (-\infty, +\infty]$$ are convex, l.s.c. and $${\it{\Lambda}}\,{\in}\, \mathcal{L}(V, Z)$$ the space of linear continuous functions from $$V$$ to $$Z.$$ Using $$\mathcal{F}^{*}$$ and $$\mathcal{G}^{*}$$ the conjugate functions (given by the Legendre–Fenchel transformation) of $$\mathcal{F}$$ and $$\mathcal{G}$$, respectively, and $${\it{\Lambda}}^{*}$$ is the adjoint operator of $${\it{\Lambda}},$$ it is not difficult to see that   supσ∈Z∗(−F∗(−Λ∗σ)−G∗(σ))≤infϕ∈VF(ϕ)+G(Λϕ), where $$Z^*$$ is the topological dual space associated with $$Z$$. This is the so called weak duality. For the strong duality, which corresponds to equality we have the following well-known result. Proposition 2.6 (cf. Ekeland & Teman, 1976) In addition, assume that there exists $$\phi_{0}$$ such that $$\mathcal{F}(\phi_0)\,{<}+\infty$$, $$\mathcal{G}({\it{\Lambda}} \phi_{0}) <+\infty,$$$$\mathcal{G}$$ being continuous at $${\it{\Lambda}} \phi_0$$. Then the Fenchel–Rockafellar dual problem   supσ∈Z∗(−F∗(−Λ∗σ)−G∗(σ)) (2.8) has at least a solution $$\sigma\in Z^*$$ and $$\inf$$ (2.7) = $$\max$$ (2.8). Moreover, in this case, $$\phi$$ is a solution to the primal problem (2.7) if and only if   −Λ∗σ∈∂F(ϕ) and σ∈∂G(Λϕ). (2.9) Proof of Theorem 2.5. We work with the uniform convergence for the spaces $$C(\overline{{\it{\Omega}}})^N$$, $$C(\overline{{\it{\Omega}}})$$ and the norm $$\|u\|_{C^1}:=\max\{\|u\|_{\infty}, \|\nabla u\|_{\infty}\}$$ for $$C^1(\overline{{\it{\Omega}}})$$. It is not difficult to see that the hypotheses of Proposition 2.6 are satisfied. Now, let us compute the Fenchel–Rockafellar dual problem of (2.3). Since $$\mathcal{F}$$ is linear, $${\mathcal{F}}^*(-{\it{\Lambda}}^*({\it{\Phi}}, \theta^0, \theta^1))$$ is finite (and always equals to $$0$$) if and only if   −Λ∗(Φ,θ0,θ1)=−(m−ν(Ω¯),ν−μ) in V∗ i.e.,   ⟨Φ,∇u⟩−⟨θ0,u⟩+⟨θ1,u−λ⟩=λ(m−ν(Ω¯))+⟨ν−μ,u⟩∀(λ,u)∈V. This implies that   ∫Ω¯∇udΦ=∫Ω¯ud(ν−θ1)−∫Ω¯ud(μ−θ0) for all u∈C1(Ω¯) and   −λ∫Ω¯dθ1=λ(m−ν(Ω¯))∀λ∈R. These mean that   −∇⋅Φ=ν−θ1−(μ−θ0) with Φ.n=0 on ∂Ω and   θ1(Ω¯)=ν(Ω¯)−m. We also have   G∗(Φ,θ0,θ1)={∫Ω¯F(x,Φ|Φ|(x))d|Φ| if θ0≥0,θ1≥0+∞ otherwise  for any (Φ,θ0,θ1)∈Z∗. Then the proof follows by Proposition 2.6. □ Remark 2.7 The optimality relations (2.9) reads   {−∇⋅Φ=ν−θ1−(μ−θ0) and Φ⋅n=0 on ∂Ωθ1(Ω¯)=ν(Ω¯)−m⟨Φ,∇u⟩≥⟨Φ,q⟩∀q∈C(Ω¯),F∗(x,q(x))≤1∀x∈Ω¯λ∈R+,u∈C1(Ω¯)⋂LdFλu=0,θ0-a.e. in Ω¯u=λ,θ1-a.e. in Ω¯.  In fact, the optimality condition $$-{\it{\Lambda}}^* \sigma \in \partial \mathcal{F}(\phi)$$ gives the first two equations and $$\sigma\in \partial \mathcal{G}({\it{\Lambda}} \phi)$$ gives the last four equations. Moreover, if $${\it{\Phi}}\in L^1({\it{\Omega}})^N$$ then the condition   ⟨Φ,∇u⟩≥⟨Φ,q⟩∀q∈C(Ω¯),F∗(x,q(x))≤1∀x∈Ω¯ can be replaced by   F(x,Φ(x))=⟨∇u(x),Φ(x)⟩ a.e. x∈Ω. (2.10) However, it is not clear in general that $${\it{\Phi}}$$ belongs to $$L^1({\it{\Omega}})^N$$. In the case where $${\it{\Omega}}$$ is convex and $$F(x, v):=|v|$$ the Euclidean norm (or some other uniformly convex and smooth norms), the $$L^p$$ regularity results are known under suitable assumptions on $$\mu$$ and $$\nu$$ (see, e.g., Feldman & McCann, 2002; De Pascale et al., 2004; De Pascale & Pratelli, 2004; Santambrogio, 2009). To our knowledge, the case of general Finsler metrics is still an open question. In the case where $${\it{\Phi}}$$ is a vector-valued measure, the condition (2.10) should be adapted to the tangential gradient. Rigorous formulations using the tangential gradient with respect to a measure, as well as rigorous proofs in the general case, can be found in the article by Igbida & Nguyen (2017) with $${\it{\Omega}}=\mathbb{R}^N$$. It is expected that $$\theta^0\le \mu$$ and $$\theta^1\le \nu$$ for optimal solutions $$({\it{\Phi}}, \theta^0, \theta^1)$$ of the minimal flow formulation (2.5). This is the case whenever $$\mathbf{m}\in [(\mu\wedge \nu)(\overline{{\it{\Omega}}}), {\mathbf{m}}_{\max}]$$, where $$\mu\wedge \nu$$ is the common mass measure of $$\mu$$ and $$\nu$$, i.e., if $$\mu, \nu\in L^1({\it{\Omega}})$$, then $$\mu \wedge \nu\in L^1({\it{\Omega}})$$ and   (μ∧ν)(x)=min{μ(x),ν(x)} for a.e. x∈Ω. In general, the measure $$\mu \wedge \nu$$ is defined by (see Ambrosio et al., 2000)   μ∧ν(A)=inf{μ(A1)+ν(A2):disjoint Borel setsA1,A2,such that A1∪A2=A}. Proposition 2.8 Let $$\mathbf{m}\in [(\mu\wedge \nu)(\overline{{\it{\Omega}}}), {\mathbf{m}}_{\max}]$$ and $$({\it{\Phi}}, \theta^0, \theta^1)\in Z^*$$ be an optimal solution of (2.5). Then $$\theta^0\le \mu$$ and $$\theta^1\le \nu$$. Moreover, $$(\mu-\theta^0, \nu -\theta^1)$$ is a couple of optimal active regions and $${\it{\Phi}}$$ is an optimal flow of transporting $$\mu-\theta^0$$ onto $$\nu -\theta^1$$. Proof. The proof follows in the same way as Theorem 5.21 and Corollary 5.20 in Igbida & Nguyen (2017). □ Our next work is to compute an approximation of $${\it{\Phi}}$$ (in fact, approximations of $${\it{\Phi}}, u, \lambda, \theta^0, \theta^1$$). To do that, we will apply an augmented Lagrangian method to the DPMK problem (2.2). 3. Discretization and convergence Coming back to the DPMK problem (2.2), our aim now is to give, by using a finite element approximation, the discretized problem associated with (2.2). To begin with, let us consider regular triangulations $$\mathcal{T}_{h}$$ of $$\overline{{\it{\Omega}}}$$. For a fixed integer $$k\ge 1$$, $$P_{k}$$ is the set of polynomials of degree less or equal $$k$$. Let $$E_{h}\subset H^{1}({\it{\Omega}})$$ be the space of continuous functions on $$\overline{{\it{\Omega}}}$$ and belonging to $$P_{k}$$ on each triangle of $$\mathcal{T}_{h}$$. We denote by $$Y_{h}$$ the space of vectorial functions such that their restrictions belong to $$(P_{k-1})^N$$ on each triangle of $$\mathcal{T}_{h}$$. Let $$f=\nu-\mu$$ and $$f_{h}\in E_{h}$$ such that $$\{f_{h}\}$$ converges weakly* to $$f$$ in $$\mathcal{M}_{b}(\overline{{\it{\Omega}}})$$. Considering the finite-dimensional spaces   Vh=R×Eh,Zh=Yh×Eh×Eh, we set   Λh(λ,u) :=(∇u,−u,u−λ)∈Zh for (λ,u)∈Vh,Fh(λ,u) :=−⟨u,fh⟩−λ(m−ν(Ω¯))∀(λ,u)∈Vh and   Gh(q,z,w):={0if z≤0,w≤0,F∗(x,q(x))≤1 for a.e. x∈Ω+∞otherwise  for (q,z,w)∈Zh. Then the finite-dimensional approximation of (2.2) reads   inf(λ,u)∈VhFh(λ,u)+Gh(Λh(λ,u)). (3.1) The following result shows that this is a suitable approximation of (2.2). Theorem 3.1 Assume that $$\mathbf{m}\,{<}\, \nu(\overline{{\it{\Omega}}})$$. Let $$(\lambda_{h}, u_{h})\,{\in}\, V_{h}$$ be an optimal solution to the approximated problem (3.1) and $$({\it{\Phi}}_{h}, \theta^0_h, \theta^1_h)$$ be an optimal dual solution to (3.1). Then, up to a subsequence, $$(\lambda_{h}, u_h)$$ converges in $$\mathbb{R}\times C(\overline {\it{\Omega}})$$ to $$(\lambda, u)$$ an optimal solution of the DPMK problem (2.2) and $$({\it{\Phi}}_h, \theta^0_h, \theta^1_h)$$ converges weakly* in $$\mathcal M_b(\overline{\it{\Omega}})^N\times \mathcal M_b(\overline{\it{\Omega}})\times \mathcal M_b(\overline{\it{\Omega}})$$ to $$({\it{\Phi}}, \theta^0, \theta^1)$$ an optimal solution of (2.5). Proof. Since $$\mathbf{m}\,{<}\,\nu(\overline{{\it{\Omega}}})$$, $$\{\lambda_h\}$$ is bounded in $$\mathbb{R}$$ and $$\{u_h\}$$ is bounded in $$(C(\overline{{\it{\Omega}}}), \|.\|_{\infty})$$. From the nondegeneracy of $$F$$ and the definitions of $$\mathcal{F}_{h}, \mathcal{G}_{h}, {\it{\Lambda}}_{h}$$, we have that $$\{u_h\}$$ is equi-Lipschitz and   uh(y)−uh(x)≤dF(x,y)∀x,y∈Ω¯. Using the Ascoli–Arzela theorem, up to a subsequence, $$u_h \rightrightarrows u$$ uniformly on $$\overline{{\it{\Omega}}}$$ and $$\lambda_h \to \lambda$$. Obviously, $$\lambda\ge 0$$ and $$u\in L^\lambda_{d_{F}}$$. Now, by the optimality of $$(\lambda_h, u_h)$$ and of $$({\it{\Phi}}_h, \theta^0_h, \theta^1_h)$$, we have   −Λh∗(Φh,θh0,θh1)=−(m−ν(Ω¯),fh) in Vh∗ and   Fh(λh,uh)+Gh(Λh(λh,uh))=−Fh∗(−Λh∗(Φh,θh0,θh1))−Gh∗(Φh,θh0,θh1). More concretely,   ⟨Φh,∇v⟩−⟨θh0,v⟩+⟨θh1,v−s⟩=s(m−ν(Ω¯))+⟨fh,v⟩∀(s,v)∈Vh, (3.2)  θh0≥0,θh1≥0,θh1(Ω¯)=ν(Ω¯)−m (3.3) and   ⟨uh,fh⟩+λh(m−ν(Ω¯))=sup{⟨q,Φh⟩:q∈Yh,F∗(x,q(x))≤1 a.e. x∈Ω}. (3.4) In (3.2), taking $$v=0$$ and $$s=1$$ (respectively, $$v=s=1$$), we see that $$\{\theta^1_h\}$$ (respectively, $$\{\theta^0_h\}$$) is bounded in $$\mathcal M_{b}(\overline{{\it{\Omega}}})$$. Moreover, using (3.4) and the boundedness of $$(\lambda_h, u_h)$$ we deduce that $$\{{\it{\Phi}}_h\}$$ is bounded in $$\mathcal{M}_b(\overline{{\it{\Omega}}})^N.$$ So, up to a subsequence,   (Φh,θh0,θh1)⇀(Φ,θ0,θ1) in Mb(Ω¯)N×Mb(Ω¯)×Mb(Ω¯)−w∗. Using (3.2) and (3.3), it is clear that $$({\it{\Phi}}, \theta^0, \theta^1)$$ satisfies   ⟨Φ,∇v⟩−⟨θ0,v⟩+⟨θ1,v−s⟩=s(m−ν(Ω¯))+⟨f,v⟩∀(s,v)∈V and   θ0≥0,θ1≥0,θ1(Ω¯)=ν(Ω¯)−m, i.e., $$({\it{\Phi}}, \theta^0, \theta^1)$$ is feasible for the minimal flow problem (2.5). Next, let us show the optimality of $$(\lambda, u)$$ and of $$({\it{\Phi}}, \theta^0, \theta^1)$$, i.e.,   ∫Ω¯F(x,Φ|Φ|(x))d|Φ|=⟨u,ν−μ⟩+λ(m−ν(Ω¯)). (3.5) We fix $$q\,{\in}\, C(\overline{{\it{\Omega}}})^N$$ such that $$F^*(x, q(x))\,{\le}\, 1 \quad\forall x\,{\in}\, \overline{{\it{\Omega}}}$$ and we consider $$q_h\,{\in}\, Y_{h}$$ such that $$\|q_{h}-q\|_{L^{\infty}({\it{\Omega}})} \to 0$$ as $$h\to 0$$. We see that   F∗(x,qh(x))=F∗(x,q(x))+F∗(x,qh(x))−F∗(x,q(x))≤1+O(h) a.e. x∈Ω. By taking $$\frac{q_{h}}{1+ O(h)}$$, we can assume that $$q_{h}\,{\in}\, Y_h, F^*(x, q_h(x))\,{\le}\, 1 \text{ a.e. } x\,{\in}\, {\it{\Omega}}$$ and $$\|q_{h}-q\|_{L^{\infty}({\it{\Omega}})} \,{\to}\, 0$$ as $$h\to 0$$. Using (3.4), we have   ⟨q,Φ⟩ =⟨qh,Φh⟩+⟨q,Φ−Φh⟩+⟨q−qh,Φh⟩ ≤sup{⟨qh,Φh⟩:qh∈Yh,F∗(x,qh(x))≤1 a.e. x∈Ω}+O(h) =⟨uh,fh⟩+λh(m−ν(Ω¯))+O(h). Letting $$h\to 0$$, we get   ⟨q,Φ⟩≤⟨u,ν−μ⟩+λ(m−ν(Ω¯)) for any q∈C(Ω¯)N,F∗(x,q(x))≤1∀x∈Ω¯. Taking supremum in $$q$$, we obtain   ∫Ω¯F(x,Φ|Φ|(x))d|Φ|≤⟨u,ν−μ⟩+λ(m−ν(Ω¯)). At last, thanks to the duality equality (2.4), this implies (3.5), the optimality of $$(\lambda, u)$$ and of $$({\it{\Phi}}, \theta^0, \theta^1)$$. □ Remark 3.2 In the case $$\mathbf{m}=\mathbf{m}_{\max}$$ (called the unbalanced transport), the DPMK problem has a simpler formulation. So for the purpose of implementation, we distinguish the two cases: the partial transport and the unbalanced transport. In the unbalanced case, let us assume that $$\mathbf{m}=\mathbf{m}_{\max}=\nu(\overline{{\it{\Omega}}})$$$$(\text{i.e., }\, \mu(\overline{{\it{\Omega}}})\ge \nu(\overline{{\it{\Omega}}}))$$, the DPMK problem (2.2) can be written as   maxu∈LipdF,u≥0∫Ω¯ud(ν−μ). (3.6) By using $$V_h=E_h,\: Z_h=Y_h \times E_h, {\it{\Lambda}}_h u =(\nabla u, -u)$$ and   Gh(q,z)={0 if z≤0,F∗(x,q(x))≤1 a.e. x∈Ω+∞ otherwise  a finite-dimensional approximation can be given by   infu∈Vh−⟨u,fh⟩+Gh(Λhu). (3.7) As in Theorem 3.1, we can prove the convergence of this finite-dimensional approximation to the original one (3.6). More precisely, we have Proposition 3.3 Assume that $$\mathbf{m}=\nu(\overline{{\it{\Omega}}})$$. Let $$u_{h}\in V_{h}$$ be an optimal solution to the approximated problem (3.7) and $$({\it{\Phi}}_{h}, \theta^0_h)$$ be an optimal dual solution to (3.7). Then, up to a subsequence and translation by constant, $$u_{h}$$ converges to $$u$$ an optimal solution of the DPMK problem (3.6) and $$({\it{\Phi}}_h, \theta^0_h)$$ converges to $$({\it{\Phi}}, \theta^0)$$ an optimal solution of (2.5) with $$\theta^1=0$$. The proof of this proposition is similar to the proof of Theorem 3.1. 4. Solving the discretized problems Our task now is to solve the finite-dimensional problems (3.1) and (3.7). First, let us recall the augmented Lagrangian method we are dealing with. 4.1 ALG2 method Assume that $$V$$ and $$Z$$ are two Hilbert spaces. Let us consider the problem   infϕ∈VF(ϕ)+G(Λϕ), (4.1) where $$\mathcal{F}: V\longrightarrow (-\infty, +\infty]$$ and $$\mathcal{G}: Z\longrightarrow (-\infty, +\infty]$$ are convex, l.s.c. and $${\it{\Lambda}}\in \mathcal{L}(V, Z)$$. We introduce a new variable $$q\in Z$$ to the primal problem (4.1) and we rewrite it in the form   inf(ϕ,q)∈V×Z:Λϕ=qF(ϕ)+G(q). The augmented Lagrangian is given by   L(ϕ,q;σ):=F(ϕ)+G(q)+⟨σ,Λϕ−q⟩+r2|Λϕ−q|2r>0. The so called ALG2 algorithm is given as follows: For given $$q_{0}, \sigma_{0}\in Z$$, we construct the sequences $$\{\phi_i\}, \{q_i\}$$ and $$\{\sigma_i\}, i=1, 2, ...,$$ by Step 1: Minimizing $$\inf\limits_{\phi} L(\phi, q_{i}; \sigma_{i})$$, i.e.,   ϕi+1∈argminϕ∈V{F(ϕ)+⟨σi,Λϕ⟩+r2|Λϕ−qi|2}. Step 2: Minimizing $$\inf\limits_{q\in Z} L(\phi_{i+1}, q; \sigma_{i})$$, i.e.,   qi+1∈argminq∈Z{G(q)−⟨σi,q⟩+r2|Λϕi+1−q|2}. Step 3: Update the multiplier $$\sigma$$,   σi+1=σi+r(Λϕi+1−qi+1). For the theory of this method and its interpretation, we refer the reader to Gabay & Mercier (1976), Glowinski et al. (1981), Fortin & Glowinski (1983), Glowinski & Le Tallec (1989) and Eckstein & Bertsekas (1992). Here, we recall the convergence result of this method which is enough for our discretized problems. Theorem 4.1 (cf. Eckstein & Bertsekas, 1992, Theorem 8) Fixed $$r>0$$, assuming that $$V=\mathbb{R}^n, Z=\mathbb{R}^m$$ and that $${\it{\Lambda}}$$ has full column rank. If there exists a solution to the optimality relations (2.9), then $$\{\phi_{i}\}$$ converges to a solution of the primal problem (2.7) and $$\{\sigma_{i}\}$$ converges to a solution of the dual problem (2.8). Moreover, $$\{q_{i}\}$$ converges to $${\it{\Lambda}} \phi^{*}$$, where $$\phi^{*}$$ is the limit of $$\{\phi_{i}\}$$. The proof of this result in the case of finite-dimensional spaces $$V$$ and $$Z$$ can be found in Eckstein & Bertsekas (1992). The result holds true in infinite-dimensional Hilbert spaces under additional assumptions. One can see Fortin & Glowinski (1983) and Glowinski & Le Tallec (1989) for more details in this direction. Next, we use the ALG2 method for the discretized problems. To simplify the notations, let us drop out the subscript $$h$$ in $$(\lambda_h, u_h)$$ and $$({\it{\Phi}}_h, \theta^0_h, \theta^1_h).$$ Thanks to Remark 3.2, we treat separately the case where $$\mathbf{m}=\nu(\overline{{\it{\Omega}}})$$ and the case where $$\mathbf{m}<\nu(\overline{{\it{\Omega}}}).$$ 4.2 Partial transport ($$\mathbf{m}<\nu(\overline{{\it{\Omega}}})$$) Given $$(q_{i}, z_{i}, w_{i}), ({\it{\Phi}}_{i}, \theta^0_i, \theta^1_i)$$ at the iteration $$i,$$ we compute Step 1:   (λi+1,ui+1) ∈argmin(λ,u)∈VhFh(λ,u)+⟨(Φi,θi0,θi1),Λh(λ,u)⟩+r2|Λh(λ,u)−(qi,zi,wi)|2 =argmin(λ,u)∈Vh−⟨u,fh⟩−λ(m−ν(Ω¯))+⟨Φi,∇u⟩+⟨θi0,−u⟩+⟨θi1,u−λ⟩  +r2|∇u−qi|2+r2|u+zi|2+r2|u−λ−wi|2. Step 2:   (qi+1,zi+1,wi+1) ∈argmin(q,z,w)∈ZhGh(q,z,w)−⟨(Φi,θi0,θi1),(q,z,w)⟩+r2|Λh(λi+1,ui+1)−(q,z,w)|2 =argmin(q,z,w)∈ZhI[F∗(.,q(.))≤1](q)+I[z≤0](z)+I[w≤0](w)−⟨Φi,q⟩−⟨θi0,z⟩−⟨θi1,w⟩  +r2|∇ui+1−q|2+r2|ui+1+z|2+r2|ui+1−λi+1−w|2. Step 3: Update the multiplier   (Φi+1,θi+10,θi+11)=(Φi,θi0,θi1)+r(∇ui+1−qi+1,−ui+1−zi+1,ui+1−λi+1−wi+1). Before giving numerical results, let us take a while to comment the above iteration. Overall, Step 1 is a quadratic programming. Step 2 can be computed easily in many cases and Step 3 updates obviously. We denote by $$\text{Proj}_{C}(.)$$ the projection onto a closed convex subset $$C$$. In Step 1, we split the computation of the couple $$(\lambda_{i+1}, u_{i+1})$$ into two steps: We first minimize w.r.t. $$u$$ to compute $$u_{i+1}$$ and then we use $$u_{i+1}$$ to compute $$\lambda_{i+1}$$. More precisely, we proceed for Step 1 as follows: (1) For $$u_{i+1}$$,   ui+1 ∈argminu∈Eh−⟨u,fh⟩+⟨Φi,∇u⟩+⟨θi0,−u⟩+⟨θi1,u⟩  +r2|∇u−qi|2+r2|u+zi|2+r2|u−λi−wi|2. This is equivalent to   r⟨∇ui+1,∇v⟩+2r⟨ui+1,v⟩ =⟨fh,v⟩−⟨Φi,∇v⟩+⟨θi0,v⟩−⟨θi1,v⟩  +r⟨qi,∇v⟩−r⟨zi,v⟩+r⟨λi+wi,v⟩∀v∈Eh. Remark here that the equation is linear with a symmetric positive definite coefficient matrix. (2) For $$\lambda_{i+1}$$, it is computed explicitly   λi+1 ∈argmins∈R−s(m−ν(Ω¯))+⟨θi1,ui+1−s⟩+r2⟨ui+1−s−wi,ui+1−s−wi⟩  =−ν(Ω¯)−m−∫Ω¯θi1+r∫Ω(wi−ui+1)r∫Ω1. In Step 2, the variables $$q, z, w$$ are independent. So, we solve them separately: (1) For $$z_{i+1}$$ and $$w_{i+1}$$, if we choose $$P_2$$ finite element for $$z_{i+1}$$ and $$w_{i+1}$$, at vertex $$x_k$$,   zi+1(xk) =Proj[r∈R:r≤0](−ui+1(xk)+θi0(xk)r) =min(−ui+1(xk)+θi0(xk)r,0) (4.2) and   wi+1(xk) =Proj[r∈R:r≤0](ui+1(xk)−λi+1+θi1(xk)r) =min(ui+1(xk)−λi+1+θi1(xk)r,0). (4.3) (2) For $$q_{i+1}$$, if we choose $$P_{1}$$ finite element for $$q_{i+1}$$ then at each vertex $$x_{l}$$  qi+1(xl)=ProjBF∗(xl,.)(∇ui+1(xl)+Φi(xl)r), (4.4) where $$B_{F^*(x, .)}:=\left\{q\in \mathbb{R}^N: F^*(x, q)\le 1\right\}$$ the unit ball for $$F^*(x, .)$$. It remains to explain how we compute the projection onto $$B_{F^*(x_l, .)}$$. This issue is recently discussed in Benamou et al. (2017) for Riemann-type Finsler distances and for crystalline norms. For the convenience of the reader, we retake here the case where the unit ball of $$F(x, .)$$ is (not necessarily symmetric) convex polytope. For short, we ignore the dependence of $$x$$ in $$F$$ and $$F^*$$. Given $$d_1, ..., d_k\ne 0$$ such that, for any $$0\ne v\in \mathbb{R}^N$$, $$\max\limits_{1\le i\le k}\left\{\langle v, d_{i}\rangle\right\}> 0$$. We consider the nonsymmetric Finsler metric given by   F(v):=max1≤i≤k{⟨v,di⟩} for any v∈RN. (4.5) It is not difficult to see that the unit ball $$B^*$$ corresponding to $$F^*$$ is exactly the convex hull of $$\{d_i\}$$,   B∗=conv(di,i=1,...,k). Thus, we need to compute the projection onto the convex hull of finite points. In dimension 2, the projection onto $$B^*$$ can be performed as follows: Compute the successive vertices $$S_1, ..., S_n$$. If $$q\,{\notin}\,B^*$$ then compute the projections of $$q$$ onto the segments $$[S_{i}, S_{i+1}]$$ and compare among these projections to chose the right one. Another way is as the one in Benamou et al. (2017): Compute outward orthogonal vectors $$v_{1}, ..., v_n$$ (Fig. 1). If $$q$$ belongs to $$[S_i, S_{i+1}]+\mathbb{R}_+v_i$$ then the projection coincides with the one on the line through $$S_i, S_{i+1}$$. If $$q$$ belongs to the sector $$S_{i}+\mathbb{R}_+v_{i-1}+\mathbb{R}_+v_{i}$$, the projection is $$S_i$$. Fig. 1. View largeDownload slide Illustration of the projection. Fig. 1. View largeDownload slide Illustration of the projection. 4.3 Unbalanced transport ($$\mathbf{m}=\nu(\overline{{\it{\Omega}}})$$) Thanks to Remark 3.2, we can reduce the algorithm in this particular case by ignoring the variable $$\lambda$$. With similar considerations for $${\it{\Lambda}}_h u=(\nabla u, -u)$$, we get the following iteration: Step 1:   ui+1 ∈argminu∈Eh−⟨u,fh⟩+⟨Φi,∇u⟩+⟨θi0,−u⟩+r2|∇u−qi|2+r2|u+zi|2. Equivalently,   r⟨∇ui+1,∇v⟩+r⟨ui+1,v⟩=⟨fh,v⟩−⟨Φi,∇v⟩+⟨θi0,v⟩+r⟨qi,∇v⟩−r⟨zi,v⟩∀v∈Eh. (4.6) Step 2: (1) For $$z_{i+1}$$, choosing $$P_{2}$$ finite element for $$z_{i+1}$$ then at each vertex $$x_k$$,   zi+1(xk)=Proj[r∈R:r≤0](−ui+1(xk)+θi0(xk)r)=min(−ui+1(xk)+θi0(xk)r,0). (4.7) (2) For $$q_{i+1}$$, choosing $$P_{1}$$ finite element, at vertex $$x_l$$,   qi+1(xl)=ProjBF∗(xl,.)(∇ui+1(xl)+Φi(xl)r). Step 3: $$({\it{\Phi}}_{i+1}, \theta^0_{i+1})=({\it{\Phi}}_{i}, \theta^0_{i})+r(\nabla u_{i+1 }- q_{i+1}, -u_{i+1}-z_{i+1}).$$ 5. Numerical experiments For the numerical implementation, we use the FreeFem++ software (Hecht, 2012) and base on Benamou & Brenier (2000) and Benamou & Carlier (2015). We use $$P_{2}$$ finite element for $$u_{i}, z_{i}, w_{i}, \theta^0_{i}, \theta^1_{i}$$ and $$P_{1}$$ finite element for $${\it{\Phi}}_{i}, q_{i}$$. 5.1 Stopping criterion In computational version, the measures $$\mu$$ and $$\nu$$ are approximated by non-negative regular functions that we denote again by $$\mu$$ and $$\nu$$. We use the following stopping criteria: For the partial transport: (1) $$\text{MIN-MAX}:=\min\left\{\min\limits_{\overline{{\it{\Omega}}}} u(x), \lambda-\max\limits_{\overline{{\it{\Omega}}}} u(x), \min\limits_{\overline{{\it{\Omega}}}} \theta^0(x), \min\limits_{\overline{{\it{\Omega}}}} \theta^1(x)\right\}\!\!.$$ (2) $$\text{Max-Lip}:=\sup\limits_{\overline{{\it{\Omega}}}} F^*(x, \nabla u(x)).$$ (3) $$\text{DIV}:=\|\nabla \cdot{\it{\Phi}} + \nu -\theta^1-\mu +\theta^0\|_{L^{2}}.$$ (4) $$\text{DUAL}:=\|F(x, {\it{\Phi}}(x)) -{\it{\Phi}}(x)\cdot \nabla u\|_{L^{2}}$$. (5) $$\text{MASS}:= \left |\int (\nu -\theta^1)\,{\rm d}x - \mathbf{m}\right |$$. For the unbalanced transport: We change (1) $$\text{MIN-MAX}:=\min\left\{\min\limits_{\overline{{\it{\Omega}}}} u(x), \min\limits_{\overline{{\it{\Omega}}}} \theta^0(x)\right\}$$. (2) $$\text{DIV}:=\|\nabla \cdot {\it{\Phi}} + \nu-\mu +\theta^0\|_{L^{2}}.$$ We expect $$\text{MIN-MAX}\ge 0, \text{Max-Lip}\le 1$$; DIV, DUAL and MASS are small. 5.2 Some examples In all the examples below, we take $${\it{\Omega}}=[0, 1]\times [0, 1]$$. We test for the Riemannian case and the crystalline case. For the latter, we consider the Finsler metric of the form $$F(x, v)=\max\limits_{1\le i\le k}\left\{\langle v, d_i\rangle\right\}$$ with given directions $$d_1, ..., d_k$$ such that for any $$0\ne v\in \mathbb{R}^2$$,   max1≤i≤k{⟨v,di⟩}>0. 5.2.1 For the unbalanced transport Example 5.1 Taking $$\mu=3{\mathcal{L}^{2}}_{{\it{\Omega}}}$$ and $$\nu=\delta_{(0.5, 0.5)}$$ the Dirac mass at $$(0.5, 0.5)$$. The Finsler metric is the Euclidean one. The optimal flow is given in Fig. 2. The stopping criterion at each iteration is given in Fig. 3. Fig. 2. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, v)=|v|$$. Fig. 2. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, v)=|v|$$. Fig. 3. View largeDownload slide Stopping criterion at each iteration. Fig. 3. View largeDownload slide Stopping criterion at each iteration. Example 5.2 We take $$\mu$$ and $$\nu$$ as in the previous example and the Finsler metric given by $$F(x, v)\,{:=}|v_1|\,{+}\,|v_2|$$ for $$v\,{=}\,(v_1, v_2)\,{\in}\, \mathbb{R}^2$$. This corresponds to the crystalline norm with $$d_1\,{=}\,(1, 1), d_2\,{=}(-1, 1), d_3\,{=}\,(-1, -1) \ \text{and} \ d_4\,{=}\,(1, -1)$$. The optimal flow is given in Fig. 4 and the stopping criterion at each iteration is given in Fig. 5. Fig. 4. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, (v_1, v_2))=|v_1|+|v_2|$$. Fig. 4. View largeDownload slide Optimal flow for $$\mu=3, \; \nu=\delta_{(0.5, 0.5)}, \: F(x, (v_1, v_2))=|v_1|+|v_2|$$. Fig. 5. View largeDownload slide Stopping criterion at each iteration. Fig. 5. View largeDownload slide Stopping criterion at each iteration. 5.2.2 For the partial transport Example 5.3 Taking $$\mu\,{=}\,4\chi_{[(x-0.3)^2+(y-0.2)^2 <0.03]}$$, and $$\nu\,{=}\,4\chi_{[(x-0.7)^2+(y-0.8)^2<0.03]}.$$ The mass of the transport is $$\mathbf{m}\,{:=}\,\frac{\nu(\overline{{\it{\Omega}}})}{2}.$$ We test for different Finsler metrics. On each figure below, the subfigure at left illustrates the unit ball of $$F$$ and the subfigure at right gives the numerical result (see Figs 6–9). The stopping criteria are summarized in Table 1. Table 1 Stopping criteria for $$800$$ iterations Case  DIV  DUAL  MASS  MIN–MAX  Max–Lip  Time execution (s)  1  2.48182e-05  9.5294e-06  0.000161361  $$-$$0.0149942  1.00068  357  2  3.38395e-05  5.58717e-05  0.000195881  $$-$$0.00120123  1.00248  867  3  7.44768e-05  5.5997e-05  6.66404e-06  $$-$$0.00272389  1.00351  1269  4  6.33726e-05  3.20691e-05  0.000120909  $$-$$0.0104915  1.02572  1123  Case  DIV  DUAL  MASS  MIN–MAX  Max–Lip  Time execution (s)  1  2.48182e-05  9.5294e-06  0.000161361  $$-$$0.0149942  1.00068  357  2  3.38395e-05  5.58717e-05  0.000195881  $$-$$0.00120123  1.00248  867  3  7.44768e-05  5.5997e-05  6.66404e-06  $$-$$0.00272389  1.00351  1269  4  6.33726e-05  3.20691e-05  0.000120909  $$-$$0.0104915  1.02572  1123  Fig. 6. View largeDownload slide Case 1: $$F(x, v)=|v|$$. Fig. 6. View largeDownload slide Case 1: $$F(x, v)=|v|$$. Fig. 7. View largeDownload slide Case 2: The crystalline case with $$d_1=(1, 1), d_2=(-1, 1), d_3=(-1, -1)$$ and $$d_4=(1, -1)$$. Fig. 7. View largeDownload slide Case 2: The crystalline case with $$d_1=(1, 1), d_2=(-1, 1), d_3=(-1, -1)$$ and $$d_4=(1, -1)$$. Fig. 8. View largeDownload slide Case 3: The crystalline case with $$d_1\,{=}\,(1, 0), d_2\,{=}\,(\frac{1}{5}, \frac{1}{5}), d_3\,{=}\,(-\frac{1}{5}, \frac{1}{5}), d_4\,{=}\,(-\frac{1}{5}, -\frac{1}{5})$$ and $$d_5\,{=}\,(\frac{1}{5}, -\frac{1}{5})$$ makes the transport more expensive in the direction of the vector $$(1, 0)$$. Fig. 8. View largeDownload slide Case 3: The crystalline case with $$d_1\,{=}\,(1, 0), d_2\,{=}\,(\frac{1}{5}, \frac{1}{5}), d_3\,{=}\,(-\frac{1}{5}, \frac{1}{5}), d_4\,{=}\,(-\frac{1}{5}, -\frac{1}{5})$$ and $$d_5\,{=}\,(\frac{1}{5}, -\frac{1}{5})$$ makes the transport more expensive in the direction of the vector $$(1, 0)$$. Fig. 9. View largeDownload slide Case 4: The crystalline case with $$d_1\,{=}\,(1, -1), d_2\,{=}\,(1, -\frac{4}{5}), d_3\,{=}\,(-\frac{4}{5}, 1), d_4\,{=}\,(-1, 1)$$ and $$d_5\,{=}\,(-1, -1)$$ makes the transport cheaper in the direction of the vector $$(1, 1)$$. Fig. 9. View largeDownload slide Case 4: The crystalline case with $$d_1\,{=}\,(1, -1), d_2\,{=}\,(1, -\frac{4}{5}), d_3\,{=}\,(-\frac{4}{5}, 1), d_4\,{=}\,(-1, 1)$$ and $$d_5\,{=}\,(-1, -1)$$ makes the transport cheaper in the direction of the vector $$(1, 1)$$. Example 5.4 Let $$\mu\,{=}\,2\chi_{[(x-0.2)^2\,{+}\,(y-0.2)^2 <0.03]}+2\chi_{[(x-0.6)^2+(y-0.1)^2 <0.01]}$$ and $$\nu\,{=}\,2\chi_{[(x-0.6)^2+(y-0.8)^2 <0.03]}$$. In this example, we take the Euclidean norm and we let $$\mathbf{m}$$ vary by taking the values $$\mathbf{m}_{i}\,{=}\,\frac{i}{6}$$$$\min\{\mu({\it{\Omega}}), \nu({\it{\Omega}})\},$$$$i=1, ..., 6.$$ The results are given in Fig. 10. Fig. 10. View largeDownload slide Optimal flows. Fig. 10. View largeDownload slide Optimal flows. Acknowledgements The authors are grateful to J. D. Benamou and G. Carlier who provide some codes of ALG2 on the link https://team.inria.fr/mokaplan/software/. Some parts of our codes are inspired from their work. Appendix Our aim here is to show Lemma A.1 that gives a smooth approximation of 1-$$d_{F}$$ Lipschitz continuous function for continuous nondegenerate Finsler metrics $$F$$. This result is more or less known in some particular cases. However, we could not find any rigorous proofs for the general case in the literature. Lemma A.1 Let $${\it{\Omega}}$$ be a connected bounded Lipschitz domain and $$F$$ be a continuous nondegenerate Finsler metric on $$\overline{{\it{\Omega}}}$$. For any Lipschitz continuous function $$u$$ on $$\overline{{\it{\Omega}}}$$ satisfying   F∗(x,∇u(x))≤1 a.e. x∈Ω, (A.1) there exists a sequence of functions $$u_{\varepsilon}\in C^{\infty}_{c}(\mathbb{R}^N)$$ such that   F∗(x,∇uε(x))≤1∀x∈Ω¯, and   uε⇉u uniformly on Ω¯. Note that $$F$$ and $$F^*$$ are defined only in $$\overline{{\it{\Omega}}}$$ and that the gradient of $$u$$ is controlled only inside of $${\it{\Omega}}$$ by (A.1). If we use the standard convolution to define $$u_{\varepsilon}$$, the value of $$u_{\varepsilon}(x)$$ is affected by the value of $$u(y)$$ outside of $$\overline{{\it{\Omega}}}$$ which remains uncontrolled. To overcome this difficulty, if $$x$$ is near the boundary, we move it a little into inside of $${\it{\Omega}}$$ before taking the convolution. To this aim, we use the smooth partition of unity tool to deal with approximation of $$u$$ near the boundary. Proof. Set   ∀x∈RN,u~(x):={u(x) if x∈Ω¯0 otherwise.  Step 1: Fix $$z\in \partial {\it{\Omega}}$$. Since $${\it{\Omega}}$$ is a Lipschitz domain, there exist $$r_{z}>0$$ and a Lipschitz continuous function $$\gamma_{z}:\mathbb{R}^{N-1} \longrightarrow \mathbb{R}$$ such that (up to rotating and relabeling if necessary)   Ω∩B(z,rz)={x|xN>γz(x1,...,xN−1)}∩B(z,rz). Set $$U_{z}:={\it{\Omega}} \cap B(z, \frac{r_{z}}{2})$$. For any $$x\in \mathbb{R}^N$$, taking   xzε:=x+ελzen, (A.2) where we choose a sufficiently large fixed $$\lambda_{z}$$ and all small $$\varepsilon$$, say fixed $$\lambda_{z} \ge \text{Lip} (\gamma_{z})+1$$, $$0<\varepsilon< \frac{r_{z}}{2(\lambda_{z}+1)}.$$ By this choice and the Lipschitz property of $$\gamma_{z}$$, we see that   B(xzε,ε)⊂Ω∩B(z,rz) for all x∈Uz. (A.3) Defining   u~ε(x):=∫RNρε(y)u~(xzε−y)dy=∫B(xzε,ε)ρε(xzε−y)u~(y)dy for all x∈RN, (A.4) where $$\rho_{\varepsilon}$$ is the standard mollifier on $$\mathbb{R}^N$$. Obviously, $$\tilde{u}_{\varepsilon} \in C^{\infty}_{c}(\mathbb{R}^N)$$. Using (A.3), (A.4) and the continuity of $$u$$ on $$\overline{{\it{\Omega}}}$$, we get   u~ε⇉u on U¯z. Step 2: Now, using the compactness of $$\partial {\it{\Omega}}$$ and $$\partial {\it{\Omega}} \,{\subset}\, \bigcup\limits_{z\,{\in}\, \partial {\it{\Omega}}} B(z, \frac{r_{z}}{2})$$, there exist $$z_{1}, ..., z_{n}\in \partial {\it{\Omega}}$$ such that   ∂Ω⊂⋃i=1nB(zi,rzi2). For short, we write $$r_{i}, U_{i}, x_{i}$$ instead of $$r_{z_{i}}, U_{z_{i}}, x_{z_{i}}$$. Taking an open set $$U_{0}\Subset {\it{\Omega}}$$ such that   Ω¯⊂⋃i=1nB(zi,ri2)⋃U0. Let $$\{\phi\}^{n}_{i=0}$$ be a smooth partition of unity on $$\overline{{\it{\Omega}}}$$, subordinate to $$\left\{U_{0}, B(z_{1}, \frac{r_{1}}{2}), ..., B(z_{n}, \frac{r_{n}}{2})\right\}$$, that is,   {ϕi∈Cc∞(RN)0≤ϕi≤1,∀i=0,...,nsupp(ϕi)⋐B(zi,ri2)∀i=1,...,n, supp(ϕ0)⋐U0∑i=0nϕi(x)=1 for all x∈Ω¯.  Because of Step 1, there exist $$\tilde{u}^1_\varepsilon, ..., \tilde{u}^n_\varepsilon \in C^{\infty}_{c}(\mathbb{R}^N)$$ such that   u~εi⇉u on U¯ii=1,...,n. For $$i=0,$$ since $$U_{0}\Subset {\it{\Omega}}$$, we can take $$\tilde{u}^0_\varepsilon:= \rho_{\varepsilon}\star \tilde{u} \in C^{\infty}_{c}(\mathbb{R}^N)$$ and $$\tilde{u}^0_{\varepsilon} \rightrightarrows u$$ on $$\overline{U}_{0}$$. Set   uε:=11+Cε+w(ε)∑i=0nϕiu~εi, where $$C$$ is chosen later and   w(ε):=sup{|F∗(x,p)−F∗(y,p)|:x,y∈Ω¯,|x−y|≤Mε,|p|≤‖∇u‖L∞} with constant $$M:=\max\limits_{1\le i\le n}\{\lambda_{z_{i}}+1\}$$, $$\lambda_{z_{i}}$$ is given in Step 1. We show that $$u_{\varepsilon}$$ satisfies all the desired properties. By the construction, $$u_{\varepsilon} \in C^{\infty}_{c}(\mathbb{R}^N)$$ and   uε⇉∑i=0nϕiu=u on Ω¯. At last, we show that $$F^*(x, \nabla u_{\varepsilon} (x))\le 1 \quad\forall x\in \overline{{\it{\Omega}}}$$. Indeed, for any $$x\in {\it{\Omega}}$$, if $$x\in U_{i}, i=1, ..., n$$ (near the boundary of $${\it{\Omega}}$$), we move $$x$$ a bit into inside of $${\it{\Omega}}$$ to $$x^\varepsilon_i:=x^\varepsilon_{z_i}$$ (see (A.2) and (A.3)), if $$x\in U_{0}$$, set $$x_0^{\varepsilon}=x$$. We have   ∇uε(x) =11+Cε+w(ε)(∑i=0n∇ϕi(x)u~εi(x)+∑i=0nϕi(x)∇u~εi(x)) =11+Cε+w(ε)(∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy  +∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)∇u(y)dy). The first sum on the right-hand side has a small norm. Indeed, using the fact that   ∑i=0n∇ϕi(x)u(x)=0 for all x∈Ω, we have   ∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy=∑i=0n∇ϕi(x)(∫B(xiε,ε)ρε(xiε−y)u(y)dy−u(x)). (A.5) Moreover,   |∫B(xiε,ε)ρε(xiε−y)u(u)dy−u(x)| ≤|∫B(xiε,ε)ρε(xiε−y)(u(y)−u(xiε))dy|+|u(xiε)−u(x)| ≤C1ε∀i=0,...,n, where the constant $$C_{1}$$ depends only on Lip($$\gamma_{z_{i}}$$) and the Lipschitz constant of $$u$$ on $$\overline{{\it{\Omega}}}$$. Thus, by combining this with (A.5),   |∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy|≤C2ε∀x∈Ω, where $$C_{2}$$ depends only on $$C_{1}$$ and $$\|\nabla \phi_{i}\|_{L^{\infty}}$$. Using the nondegeneracy of $$F$$, we have   F∗(x,∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy)≤C3ε for all x∈Ω. Fixed any $$x\in {\it{\Omega}}$$, if $$y\in B(x_i^{\varepsilon}, \varepsilon)$$ then $$|x-y|\le |x-x_i^{\varepsilon}|+|x_i^{\varepsilon}-y|\le M\varepsilon$$. So we obtain   F∗(x,∇uε(x)) ≤11+Cε+w(ε)[F∗(x,∑i=0n∇ϕi(x)∫B(xiε,ε)ρε(xiε−y)u(y)dy)  +F∗(x,∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)∇u(y)dy)] ≤11+Cε+w(ε)(C3ε+∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)F∗(x,∇u(y))dy) ≤11+Cε+w(ε)[C3ε+∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)F∗(y,∇u(y))dy  +∑i=0nϕi(x)∫B(xiε,ε)ρε(xiε−y)(F∗(x,∇u(y))−F∗(y,∇u(y)))dy] ≤C3ε+1+w(ε)1+Cε+w(ε) ≤1(choose a constant C≥C3). By the continuity of $$\nabla u_{\varepsilon}$$ and of $$F^*$$, we also have $$F^*(x, \nabla u_{\varepsilon}(x))\le 1 \, \quad\forall x\in \overline{{\it{\Omega}}}.$$ □ Proposition A.2 Let $$F$$ be a continuous nondegenerate Finsler metric on a connected bounded Lipshitz domain $${\it{\Omega}}$$. We have   LipdF={u:Ω¯⟶R|u is Lipschitz continuous and F∗(x,∇u(x))≤1 a.e. x∈Ω}:=BF∗. (A.6) As a consequence, for any 1-$$d_{F}$$ Lipschitz continuous function $$u$$, there exists a sequence of 1-$$d_{F}$$ Lipschitz continuous functions $$u_{\varepsilon}\in C^{\infty}_{c}(\mathbb{R}^N)$$ and $$u_{\varepsilon} \rightrightarrows u$$ uniformly on $$\overline{{\it{\Omega}}}$$. Lemma A.3 We have $$Lip_{d_{F}}\subset \mathcal{B}_{F^*}.$$ Proof. Let $$u\in Lip_{d_{F}}$$. Then $$u$$ is Lipschitz and $$u$$ is differentiable a.e. in $${\it{\Omega}}$$. Let $$x\in {\it{\Omega}}$$ be any point where $$u$$ is differentiable. We have, for any $$v\in \mathbb{R}^N$$,   ⟨∇u(x),v⟩F(x,v) =limh→0u(x+hv)−u(x)F(x,hv) ≤lim suph→0dF(x,x+hv)F(x,hv) ≤lim suph→0∫01F(x+thv,hv)dtF(x,hv)=1. Hence, $$F^*(x, \nabla u(x))\le 1$$. So $$u\in \mathcal{B}_{F^*}$$. □ Lemma A.4 We have $$\mathcal{B}_{F^*}\subset Lip_{d_{F}}.$$ Proof. Fix any $$u\in \mathcal{B}_{F^*}$$. Case 1: If $$u$$ is smooth then $$F^*(x, \nabla u(x))\le 1 \quad\forall x\in \overline{{\it{\Omega}}}$$. For any $$x, y \in \overline{{\it{\Omega}}}$$ and any Lipschitz curve $$\xi$$ in $$\overline{{\it{\Omega}}}$$ joining $$x$$ and $$y$$, we have   u(y)−u(x) =∫01∇u(ξ(t))ξ˙(t)dt ≤∫01F∗(ξ(t),∇u(ξ(t)))F(ξ(t),ξ˙(t))dt ≤∫01F(ξ(t),ξ˙(t))dt. Hence $$u\in Lip_{d_{F}}$$. Case 2: For general Lipschitz continuous function $$u$$ satisfying $$F^*(x, \nabla u(x))\,{\le}\, 1 \, \text{ a.e. }\, x\in {\it{\Omega}}$$, thanks to Lemma A.1, there exist $$u_{\varepsilon} \in \mathcal{B}_{F^*}\bigcap C^{\infty}_{c}(\mathbb{R}^N)$$ such that $$u_{\varepsilon} \,{\rightrightarrows}\, u$$ on $$\overline{{\it{\Omega}}}$$. According to Case 1 above, $$u_{\varepsilon}\in Lip_{d_{F}}$$. Since $$u_{\varepsilon} \rightrightarrows u$$ on $$\overline{{\it{\Omega}}}$$, we obtain $$u\in Lip_{d_{F}}$$. □ Proof of Proposition A.2. The proof follows by Lemma A.3 and Lemma A.4. □ Proof of Lemma 2.4. Since $$0\le u\le \lambda$$, the sequence $$u_{\varepsilon}$$ in the proof of Lemma A.1 satisfies $$0\le u_{\varepsilon}\le \lambda$$. So $$u_{\varepsilon}\in C^{\infty}_{c}(\mathbb{R}^N) \cap L^\lambda_{d_{F}}$$ and $$u_{\varepsilon} \rightrightarrows u$$ on $$\overline{{\it{\Omega}}}$$. □ Remark A.5 The results still hold true if $${\it{\Omega}}$$ is connected, bounded and has the segment property. References Ambrosio, L., Fusco, N. & Pallara, D. ( 2000) Functions of Bounded Variation and Free Discontinuity Problems.  Oxford Mathematical Monographs. Oxford: Oxford University Press. Ambrosio, L., Gigli, N. & Savaré, G. ( 2005) Gradient Flows in Metric Spaces and in the Space of Probability Measures. Lectures in Mathematics.  Basel: Birkháuser. Barrett, J. W. & Prigozhin, L. ( 2009) Partial L1 Monge-Kantorovich problem: Variational formulation and numerical approximation. Interface. Free Bound.,  11, 201– 238. Google Scholar CrossRef Search ADS   Beckmann, M. ( 1952) A continuous model of transportation. Econometrica,  20, 643– 660. Google Scholar CrossRef Search ADS   Benamou, J. D. & Brenier, Y. ( 2000) A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem. Numer. Math.,  84, 375– 393. Google Scholar CrossRef Search ADS   Benamou, J. D. & Carlier, G. ( 2015) Augmented Lagrangian methods for transport optimization, mean field games and degenerate elliptic equations. J. Optim. Theory Appl. , 167, 1– 26. Google Scholar CrossRef Search ADS   Benamou, J. D., Carlier, G., Cuturi, M., Nenna, L. & Peyré, G. ( 2015) Iterative Bregman projections for regularized transportation problems. SIAM J. Sci. Comput.,  37, A1111– A1138. Google Scholar CrossRef Search ADS   Benamou, J. D., Carlier, G. & Hatchi, R. ( 2017) A numerical solution to Monge’s problem with a Finsler distance cost. ESAIM: M2AN,  DOI: http://dx.doi.org/10.1051/m2an/2016077. Caffarelli, L. & McCann, R. J. ( 2010) Free boundaries in optimal transport and Monge-Ampere obstacle problems. Ann. Math.,  171, 673– 730. Google Scholar CrossRef Search ADS   Chizat, L., Peyré, G., Schmitzer, B. & Vialard, F. X. ( 2016) Scaling algorithms for unbalanced transport problems. arXiv preprint arXiv:1607.05816. De Pascale, L., Evans, L. C. & Pratelli, A. ( 2004) Integral estimates for transport densities. Bull. London Math. Soc.,  36, 383– 395. Google Scholar CrossRef Search ADS   De Pascale, L. & Pratelli, A. ( 2004) Sharp summability for Monge transport density via interpolation. ESAIM Contr. Optim. Ca. Va. , 10, 549– 552. Google Scholar CrossRef Search ADS   Eckstein, J. & Bertsekas, D. P. ( 1992) On the Douglas–Rachford splitting method and the proximal point algorithm for maximal monotone operators. Math. Program.,  55, 293– 318. Google Scholar CrossRef Search ADS   Ekeland, I. & Teman, R. ( 1976) Convex Analysis and Variational Problems.  Amsterdam-New York: North-Holland American Elsevier. Feldman, M. & McCann, R. J. ( 2002) Uniqueness and transport density in Monge’s mass transportation problem. Calc. Var. Partial Differ. Equ. , 15, 81– 113. Google Scholar CrossRef Search ADS   Figalli, A. ( 2010) The Optimal Partial Transport Problem. Arch. Rational Mech. Anal. , 195, 533– 560. Google Scholar CrossRef Search ADS   Fortin, M. & Glowinski, R. ( 1983) Augmented Lagrangian methods: Applications to the Numerical Solution of Boundary-Value Problems , vol 15. Amsterdam: North-Holland Publishing. Gabay, D. & Mercier, B. ( 1976) A dual algorithm for the solution of nonlinear variational problems via finite element approximation. Comput. Math. Appl. , 2, 17– 40. Google Scholar CrossRef Search ADS   Glowinski, R., Lions, J. L. & Tre1molie1res, R. ( 1981) Numerical Analysis of Variational Inequalities.  Amsterdam: North-Holland Publishing. Glowinski, R. & Le Tallec, P. ( 1989) Augmented Lagrangian and Operator-Splitting Methods in Nonlinear Mechanics,  vol. 9. Philadelphia: SIAM. Google Scholar CrossRef Search ADS   Hecht, F. ( 2012) New development in freefem++. J. Numer. Math.,  20, 251– 266. Google Scholar CrossRef Search ADS   Igbida, N. & Nguyen, V. T. ( 2017) Optimal partial mass transportation and obstacle Monge-Kantorovich equation (submitted for publication). Igbida, N. & Ta Thi, N. N. ( 2017) Sub-gradient Diffusion Operator. J. Differ. Equ . 262, 3837– 3863. Google Scholar CrossRef Search ADS   Santambrogio, F. ( 2009) Absolute continuity and summability of transport densities: simpler proofs and new estimates. Calc. Var. Partial Differ. Equ.,  36, 343– 354. Google Scholar CrossRef Search ADS   Santambrogio, F. ( 2015) Optimal Transport for Applied Mathematicians . Basel: Birkhäuser. Google Scholar CrossRef Search ADS   Villani, C. ( 2003) Topics in Optimal Transportation . Graduate Studies in Mathematics, vol. 58. Providence: American Mathematical Society. Villani, C. ( 2009) Optimal Transport, Old and New . Grundlehren des Mathematischen Wissenschaften (Fundamental Principles of Mathematical Sciences), vol. 338. Berlin, Heidelberg: Springer. © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations

Abstract access only

18 million full-text articles

Print

20 pages / month

PDF Discount

20% off