Source-linked AI summary
Convergence of a Newton algorithm for semi-discrete optimal transport
Jun Kitagawa, Quentin Mérigot, Boris Thibert
TL;DR
Semi-discrete optimal transport has methods with convergence analyses that are slow and faster Newton-type methods without guarantees. This paper introduces a damped Newton algorithm under cost-regularity and source-connectivity assumptions. It proves global linear and local rate-1+α convergence, while noting that its convergence estimate deteriorates as the target size grows.
Problem
Semi-discrete optimal transport algorithms have traded convergence-speed analysis for practical speed, leaving fast Newton-type methods without convergence guarantees.
Source: §1 · p. 1Method
The paper analyzes a damped Newton algorithm under cost regularity and a weighted Poincaré-Wirtinger condition on the source density.
Source: §1 · p. 2Results
The damped Newton algorithm converges globally with a linear rate and locally with rate 1+α under the stated assumptions.
Source: §1.7 · p. 7Takeaways & Limitations
The analysis establishes global convergence with optimal convergence rates for the damped Newton algorithm under its assumptions.
Source: §1.7 · p. 7The convergence estimate degrades as the number of target points grows, and a space-discretization of the proposed continuous approach remains open.
Source: §1.7 · p. 7
Abstract
from arXiv · showhide
Many problems in geometric optics or convex geometry can be recast as optimal transport problems: this includes the far-field reflector problem, Alexandrov's curvature prescription problem, etc. A popular way to solve these problems numerically is to assume that the source probability measure is absolutely continuous while the target measure is finitely supported. We refer to this setting as semi-discrete optimal transport. Among the several algorithms proposed to solve semi-discrete optimal transport problems, one currently needs to choose between algorithms that are slow but come with a convergence speed analysis (e.g. Oliker-Prussner) or algorithms that are much faster in practice but which come with no convergence guarantees Algorithms of the first kind rely on coordinate-wise increments and the number of iterations required to reach the solution up to an error of $\epsilon$ is of order $N^3/\epsilon$, where $N$ is the number of Dirac masses in the target measure. On the other hand, algorithms of the second kind typically rely on the formulation of the semi-discrete optimal transport problem as an unconstrained convex optimization problem which is solved using a Newton or quasi-Newton method. The purpose of this article is to bridge this gap between theory and practice by introducing a damped Newton's algorithm which is experimentally efficient and by proving the global convergence of this algorithm with optimal rates. The main assumptions is that the cost function satisfies a condition that appears in the regularity theory for optimal transport (the Ma-Trudinger-Wang condition) and that the support of the source density is connected in a quantitative way (it must satisfy a weighted Poincar\'e-Wirtinger inequality).
1. Introduction
Semi-discrete optimal transport models an absolutely continuous source and a finitely supported target, but existing algorithms trade convergence guarantees for practical speed. The paper introduces a damped Newton algorithm intended to bridge that gap, using Laguerre cells to represent how prices alter assignments.
1. Introduction
- Existing semi-discrete OT methods trade speed with guarantees against faster Newton or quasi-Newton methods without convergence guarantees.Coordinate-wise methods require O(N^3/ε) iterations to reach error ε, where N is the number of target Dirac masses.Source: §1 · p. 1
1.1. Semi-discrete optimal transport.
- Semi-discrete optimal transport assumes an absolutely continuous source probability measure and a finitely supported target measure.The setting is used to model problems in geometric optics and convex geometry, including far-field reflectors and Alexandrov’s curvature prescription.Source: §1 · p. 1
1.2. Laguerre tessellation and economic interpretation.
- Laguerre cells assign locations by minimizing travel cost plus bakery price, so raising a bakery’s price shrinks its cell.Under the twist condition, the induced transport map is uniquely defined almost everywhere, and each cell’s source mass determines the mass sent to its bakery.Source: §1.2 · p. 3, Figure 1 · p. 3
1.3. Kantorovich’s functional.
The paper uses Kantorovich’s functional to express semi-discrete transport through price vectors and Laguerre-cell masses. It proposes a damped Newton iteration and establishes that suitable prices induce an optimal map matching the target measure.
1.3. Kantorovich’s functional.
- For any target probability measure on the finite set, a price function exists whose induced transport map pushes the source density to that target.The map induced by Laguerre cells is optimal for its push-forward measure.Source: §1.3 · p. 3
- Kantorovich’s functional is concave and C1-smooth, and its maximizers correspond to optimal transport maps matching the target measure.The functional and its gradient are invariant to adding a constant to all prices.Source: §1.3 · p. 4
1.4. Damped Newton algorithm.
- The proposed algorithm uses a pseudo-inverse Newton direction and halves its step until the specified line-search condition is met.It updates the price vector by the selected damped step.Source: §1.4 · p. 4, Step 1 · p. 5, Step 3 · p. 5
- Global convergence and convergence-speed estimates rely on regularity and strong monotonicity of the Laguerre-cell mass map.The paper links regularity to cost geometry and density regularity, and strong monotonicity to a weighted Poincaré-Wirtinger inequality for the source support.Source: §1.4 · p. 4
1.5. Regularity of Kantorovich’s functional and MTW condition.
The paper connects cost geometry to the smoothness of Kantorovich’s functional and uses a weighted Poincaré-Wirtinger condition to obtain strict concavity. These properties provide the regularity and monotonicity ingredients used in the convergence analysis.
1.5. Regularity of Kantorovich’s functional and MTW condition.
- Higher-order regularity of Kantorovich’s functional depends on cost geometry; the paper uses a discretized Loeper condition as a sufficient condition.The condition is related to the Ma-Trudinger-Wang condition in smooth optimal-transport regularity theory.Source: §1.5 · pp. 5–6
- Under (Reg), (Twist), and (QC), with a compact c-convex domain and a Hölder-continuous density, Kantorovich’s functional is C2,α.The stated regularity estimate is uniform when all Laguerre cells have mass bounded below by a positive constant.Source: §1.5 · p. 6
1.6. Strong concavity of Kantorovich’s functional.
- A disconnected density support can disconnect the Hessian’s induced graph and obstruct strict concavity; the paper assumes a weighted Poincaré-Wirtinger inequality to address this.The inequality is imposed on the density and source set, not simply as a requirement that the support be simply connected.Source: §1.6 · p. 6
- Under (Reg), (Twist), and (QC), a continuous density satisfying (PW) makes Kantorovich’s functional strictly concave on the zero-mean price vectors with positive-mass cells.A quantitative strong-concavity statement assumes a positive lower bound on the mass of every Laguerre cell.Source: §1.6 · p. 7
1.7. Convergence result.
Theorem 1.5 establishes global convergence of the damped Newton algorithm with linear rate and local rate 1 + α under regularity, cost, and weighted Poincaré–Wirtinger assumptions. The result extends prior convergence strategies to semi-discrete optimal transport, while the damping lower bound weakens as the target support grows.
1.7. Convergence result.
- The damped Newton algorithm converges globally at a linear rate and locally at rate 1 + α under the theorem’s stated assumptions.These include (Reg), (Twist), (QC), a compact c-convex set X, a C0,α density, and a positive Poincaré–Wirtinger constant.Source: §1.7 · p. 7
- The theorem requires no convexity or c-convexity of the source density’s support, unlike other numerical methods for Monge–Ampère equations [6, 24].The paper gives an annular radial-density example whose Poincaré–Wirtinger constant is positive.Source: §1.7 · p. 7
- The theorem’s positive lower bound on the damping parameter degrades as the number of target masses N grows.Controlling it for large N is described as plausible but indirect; the space-discretization of the proposed continuous approach remains open.Source: §1.7 · p. 7
Comparison to previous work.
- Reported Newton and quasi-Newton methods for discretized Monge–Ampère or transport problems [27, 11, 6] lack convergence proofs.The paper distinguishes its convergence theorem from those reports of good practical behavior.Source: Comparison to previous work · pp. 7–8
- Loeper and Rapetti’s damped Newton method for quadratic transport on the torus requires densities bounded below, whereas Theorem 1.5 does not.Their damping parameter cannot be determined explicitly from the data, and their proof concerns continuous densities.Source: Comparison to previous work · p. 8
- The convergence strategy extends Mirebeau’s damped-Newton approach by replacing Dirichlet conditions with optimal-transport boundary conditions, allowing general probability densities, and accommodating more general costs.Mirebeau’s method established global convergence without rates.Source: Comparison to previous work · p. 8
- The paper’s uniform concavity result answers a Gangbo–McCann conjecture quantitatively for costs, measures, and domains satisfying Theorem 5.1’s assumptions.The conjecture concerned convergence of the negative gradient flow of the Kantorovich functional.Source: Comparison to previous work · p. 8
2. Kantorovich’s functional
The section formulates semi-discrete optimal transport through the Kantorovich functional and establishes its concavity, smoothness, and gradient. Its proof uses optimality of the induced transport map and continuity of the cell-mass map to identify a continuous supergradient selection.
2. Kantorovich’s functional
- For bounded source densities and a finitely supported target measure, the Kantorovich functional is concave, C1-smooth, and has the stated gradient.The section adapts the variational formulation from to costs satisfying (Reg’) and (Twist’).Source: §2 · p. 9
- For any potential ψ, the map Tψ is optimal for the cost c between ρ and its pushforward Tψ#ρ.This optimality supports the variational argument for concavity.Source: §2 · p. 9
- Under (Twist’) and (Reg’), the map G from potentials to cell masses is continuous for a density on compact X.The proof controls changes in Laguerre cells as the potential varies.Source: §2 · p. 10
2.1. Proof of Theorem 1.1.
- The proof obtains concavity by integrating the pointwise transport inequality, with equality for Tψ, to construct a supergradient of Φ.Continuity of G then supplies a continuous selection from the superdifferential.Source: §2.1 · pp. 10–11
3. Local regularity in a c-exponential chart
The section gives sufficient geometric conditions for regularity of a cell-mass map in a c-exponential chart. Non-degeneracy and transversality yield C1,α regularity, with transversality controlling singular boundary intersections through a Hausdorff-measure bound.
3. Local regularity in a c-exponential chart
- The chart’s regularity problem is to find sufficient conditions for ˆG near the origin when the functions defining the cells are quasi-convex.The section treats compact convex ˆX and a continuous probability density.Source: §3 · p. 11
Non-degeneracy.
- The non-degeneracy condition bounds the gradients of the defining functions away from zero and is necessary for continuity of ˆG, even when N = 1.The section states that the conditions are met for functions arising from transport costs satisfying Loeper’s condition.Source: Non-degeneracy · p. 11, §3.1 · p. 11
Transversality.
- Transversality ensures a lower bound on angles between adjacent cell-boundary facets as λ varies near λ0.The condition applies where two defining-function boundaries meet or where one meets ∂ˆX.Source: Transversality · pp. 11–12
- Non-degeneracy and transversality near λ0 imply that the map ˆG is C1,α on the specified cube, with partial derivatives given by the theorem.The stated hypotheses include a C0,α probability density on ˆX.Source: Transversality · p. 12
3.2. Sketch of proof.
- The C1,α estimate depends on the transversality constant εtr but not on Ttr.The sketch attributes the derivative formula to the co-area formula and the Hölder estimate to control of singular intersections.Source: Transversality · p. 12, §3.2 · pp. 12–13
3.3. A control on the (d −2)–Hausdorff measure of singular points.
- Transversality yields a uniform upper bound on the (d − 2)-Hausdorff measure of the singular set Σ(λ).The proof reduces the control to a bound for τ-singular points of a compact convex body.Source: §3.3 · p. 13
Step 2.
- The singular set can be covered by Nε balls of radius ε with Nε ≤ C(d, diam(K))/(τε^(d−2)).The argument uses packing and covering numbers to bound the Hausdorff measure.Source: Step 2 · p. 14
- The singular set Σ(λ) lies within the τ-singular points of K(λ) for τ = εtr, linking the convex-body estimate to the regularity proof.The conclusion follows by applying the bound from Proposition 3.3.Source: §3.3 · p. 13, Step 2 · p. 14
4. C2,α regularity of Kantorovich’s functional
This section establishes uniform C²,α regularity of Kantorovich’s functional on potentials whose Laguerre cells have mass above ε, under the stated assumptions, and develops an alternative upper-transversality bound for interior cell-boundary points, its proof geometry, and cost conditions implying the estimates used.
4. C2,α regularity of Kantorovich’s functional
- Φ is uniformly C2,α on Kε, with Hessian given by (1.8) and a norm bound depending only on ∥ρ∥∞, ε, diam(X), and constants specified in Remark 4.1.The theorem assumes compact c-convex X, ρ ∈ Pac(X) ∩ C0,α(X), α ∈ (0, 1], and (Reg), (Twist), and (QC). Kε consists of potentials for which every Laguerre cell has mass greater than ε.Source: §4 · pp. 14–15
- Even with (Reg), (Twist), and (QC), Φ need not be differentiable when a Laguerre cell has zero volume, so a positive cell-volume lower bound is necessary for second differentiability.The Euclidean example with Y = {−1, 0, 1} and ρ = 1 yields nondifferentiability at τ = 1/2.Source: §4 · p. 15
4.4. Alternative upper transversality estimates.
- An alternative upper transversality estimate is proposed when Y is sampled from a bounded target domain Λ and the stated exponential-map and convexity conditions hold.The estimate is introduced as Proposition 4.8; the supplied passages state its assumptions but not its bound.Source: Outline · p. 16, §4.4 · p. 23
4.1. Lower transversality estimates.
- For ψ ∈ Kε/2, Proposition 4.2 bounds how close two unit outward normals to a Laguerre cell can be to opposite directions, using convexity and the cell’s positive volume.The proof bounds the cell volume above using its diameter and the angle between the normals, then compares this with the lower bound implied by ε.Source: §4.1 · p. 16
4.2. Upper transversality estimates.
- The first quantitative upper transversality estimate uses the condition number of the differential of exponential maps and does not require Y to be sampled from a continuous domain.The section also announces a second estimate in §4.4.Source: §4.2 · p. 17
Notation.
- Using the positive-mass cell bound, Proposition 4.5 derives transversality estimates for the functions built from the cost, transferring normal-vector bounds across exponential charts.The argument uses the lower transversality result and lemmas on how diffeomorphisms and linear maps affect normal cones and angles.Source: Notation · pp. 18–19
4.3. Proof of Theorem 4.1.
- For any y0 ∈ Y and ψ0 ∈ Kε, Φ is C2,α-regular on Kε; the claimed dependency of ∥Φ∥C2,α(Kε) follows from equations (4.31)–(4.32) and Theorem 3.1.Source: §4.3 · p. 22
4.4. Alternative upper transversality estimates.
- Proposition 4.8 gives an alternative bound in (T) when p0 is interior to X, including for Laguerre cells not meeting ∂X, without requiring Ccond.The bound’s stated advantage is that it does not require knowledge of the condition number.Source: §4.4 · p. 23
- Both upper-transversality proofs use convex sublevel sets to rule out a zero-mass Laguerre cell, but the continuous proof rotates gradients instead of changing coordinates.The discrete proof quantifies the estimate using a cone argument or the lower-transversality bound after a coordinate change, which can introduce the condition number; the continuous proof estimates rotation error.Source: §4.4 · p. 24
- Under standard cost and domain assumptions, (QC) and (4.34) hold; the cited results also establish their equivalence in the Euclidean case.The assumptions include bounded smooth domains, a C4 cost, twist conditions, convexity of the exponential-map domains, nonzero det D2_xy c, and (A3w); the cited results are [23, Theorem 3.2] and [20, Theorem 4.10], with equivalence noted in.Source: §4.4 · p. 26
5. Strong concavity of Kantorovich’s functional
Under the stated regularity and weighted Poincaré–Wirtinger assumptions, the section establishes strong concavity of the Kantorovich functional by relating its Hessian to a weighted graph Laplacian. The proof connects continuous and discrete Cheeger bounds to control the Hessian’s first nonzero eigenvalue.
5. Strong concavity of Kantorovich’s functional
- Theorem 5.1 gives a positive lower bound for the first nonzero eigenvalue of the Kantorovich functional’s Hessian under (Reg), (Twist), (QC), and (PW).The result assumes X is compact and c-convex and ρ is a continuous probability density.Source: §5 · p. 27
- The support of ρ need not be c-convex; an annular radial measure can have a positive Poincaré–Wirtinger constant.The upper bound on the largest nonzero Hessian eigenvalue degrades as N grows because ε is of order 1/N; the stated reverse isoperimetric estimate is also acknowledged as crude.Source: §5 · p. 27
5.1. Poincar´e inequality and continuous Cheeger constant.
- A finite weighted Poincaré–Wirtinger constant implies that the weighted Cheeger constant of (X, ρ) is positive.The argument passes through Lipschitz domains and functions of bounded variation.Source: §5.1 · p. 28
5.2. Cheeger constant of a graph.
- The weighted graph Cheeger inequality bounds the first nonzero Laplacian eigenvalue using its Cheeger constant and minimum weighted degree.The graph Laplacian is the matrix used to connect the discrete Cheeger estimate to the functional’s Hessian.Source: §5.2 · p. 29
5.3. Proof of Theorem 5.1.
- The Hessian −D2Φ(ψ) is the Laplacian of a graph whose edge weights are determined by shared Laguerre-cell facets.The proof bounds the graph Cheeger constant and minimum degree from below using continuous weighted Cheeger control and cell-mass assumptions.Source: §5.3 · pp. 29–30
5.1. Combining (5.42), (5.43) and |A|ρ = P
- The resulting graph-eigenvalue lower bound yields the claimed Hessian eigenvalue bound on the subspace orthogonal to constants.The Hessian’s kernel is the constant functions when the induced graph is connected.Source: §5.1 · p. 31
6. Convergence of the damped Newton algorithm
The section establishes global convergence guarantees for a damped Newton algorithm, addressing the lack of guarantees for faster practical approaches. Under regularity and strong monotonicity conditions, it gives global linear and local superlinear convergence, then illustrates the algorithm on a nonconvex, non-simply-connected source.
6. Convergence of the damped Newton algorithm
- The section derives the semi-discrete result by applying the general damped Newton theorem to the gradient of the Kantorovich functional.The preceding sections supply the required regularity and strict concavity properties.Source: §6 · p. 31, §6.2 · p. 33
6.1. General damped Newton algorithm.
- The damped Newton algorithm converges globally at a linear rate and locally at a superlinear rate, quadratic when α = 1.The general result applies to maps satisfying regularity and uniform monotonicity assumptions.Source: §6.1 · p. 32
6.2. Proof of Theorem 1.5.
- For the Kantorovich functional, regularity and the weighted Poincaré–Wirtinger assumption establish the conditions needed to apply the general convergence result.The paper invokes Theorem 4.1 for uniform C1,α regularity and Theorem 5.1 for uniform monotonicity.Source: §6.2 · p. 33
6.3. Numerical results.
- 25 iterations reached machine-precision error on a quadratic-cost example with a piecewise-linear source density and a uniform-grid target.The source support is neither simply connected nor convex.Source: §6.3 · pp. 33–34
Appendix A. A weighted Poincar´e-Wirtinger inequality
Appendix A constructs a radial probability density supported on an annulus and shows that it satisfies a weighted Poincaré-Wirtinger inequality. The proof combines inequalities on the sphere and on a radial segment.
Appendix A. A weighted Poincar´e-Wirtinger inequality
- A radial probability density with profile zero on [0,r] and concave on [r,R] satisfies the weighted Poincaré-Wirtinger inequality for some positive constant.The example has annular, and therefore not simply connected, support.Source: Appendix A · p. 34
- The proof uses Poincaré-Wirtinger inequalities on the sphere and on [0,R] weighted by ρ; the segment inequality follows from Theorem 2.1 in and ρ's concavity.Source: Appendix A · p. 34
- Polar coordinates split the estimate into sphere and radial contributions, which the proof bounds using the two inequalities.Source: Appendix A · pp. 34–36
B.1. Existence of partial derivatives.
This section establishes the existence of partial derivatives for ˆG by showing continuity of their integral expressions. The proof uses coarea and controls changes in the facets defining those expressions.
B.1. Existence of partial derivatives.
- The finite difference of ˆG in direction e_i is rewritten with the coarea formula, reducing derivative existence to continuity of an integral over level-set facets.Source: §B.1 · p. 37
- Continuity of ˆg near λ_i yields the partial derivative formula (3.14) for ˆG.Taking σ = ˆρ/∥∇f_i∥ makes ˆg continuous under the non-degeneracy condition and f_i ∈ C1,1.Source: §B.1 · p. 37
- Proposition B.1 establishes continuity by comparing facet integrals at t and s and bounding contributions from boundary-crossing subsets.The estimates use transversality and constants depending on the regularity and geometric bounds listed in the proof.Source: §B.1 · pp. 38, 39, 42
- A flow along the level-set direction maps corresponding facets between levels and supports the comparison of their integrals.Source: §B.1 · pp. 38–39
B.2. Continuity of partial derivatives.
Proposition B.3 establishes uniform continuity of H on Q. Its proof combines continuity under changes in the ith variable with a Lipschitz bound for changes in another variable.
B.2. Continuity of partial derivatives.
- H is uniformly continuous on Q with a modulus controlled by ∥f_i∥C1,1, diam(ˆX), ε_nd, ε_tr, and ∥σ∥∞.Source: §B.2 · p. 42
- For changes in a variable other than i, the proof bounds the corresponding function h as Lipschitz, then combines coordinate-wise estimates to obtain continuity of H.Source: §B.2 · pp. 42–43
B.3. Proof of Theorem 3.1.
The proof obtains C0,α regularity for each partial derivative of ˆG on Q from the C0,α regularity of ˆρ.
B.3. Proof of Theorem 3.1.
- Each partial derivative of ˆG is C0,α on Q, with its C0,α constant controlled by the stated data.Source: §B.3 · p. 43