PaperlayerView exact revision on arXiv

Source-linked AI summary

A convex numerical scheme for the Evans-Gangbo transport-density system, in one, two, and three dimensions

Babacar Mbaye Ndiaye, Diaraf Seck

arXiv:2609.29987v1math.NAcs.NA

TL;DR

The paper addresses how to recover particle travel distances and construct Monge maps beyond the earlier unfinished 2D treatment. It completes the 2D construction with direct and ray-based methods, then extends the framework to 3D, where the tested examples show recurring boundary-concentrated transport.

  • Problem

    Kantorovich duality identifies a potential but does not directly specify each particle’s travel distance, and an earlier report left the 2D optimal-map construction unfinished.

    Source: §1 · p. 1, §1.1 · p. 2

  • Method

    The paper completes the 2D construction using direct ODE integration and a ray-based method, cross-validating with strong duality, mass conservation, and regularization limits.

    Source: §8 · p. 24

  • Results

    Both 2D examples show boundary-concentrated mass; in 3D, the solver matches the 2D mass-weighted mean transport distance to all six reported digits and reproduces the same phenomena on a radial ball.

    Source: §8 · p. 24

  • Takeaways & Limitations

    The 3D extension preserves the underlying theory, with increased computational cost from discretization.

    Source: §8 · p. 24

    The idealized closed form is not a tight quantitative prediction for the large-pile, square-domain radial example; its quantitative validation instead uses symmetry, strong duality, and mass conservation.

    Source: §5.2 · p. 17

Abstract

from arXiv · show

We give a constructive numerical treatment of Monge's optimal transport problem in 1D, 2D, 3D, following Evans-Gangbo: the map is built from a Kantorovich potential $u$ and density $a$ solving $-\mathrm{div}(a\nabla u) = f^+-f^-$, via $\dot T = -a\nabla u(T)/[(1-t)f^+(T)+tf^-(T)]$, $T(0,x)=x$. In 1D this gives a closed form. In 2D--left unfinished in an earlier report--we obtain $(u,a)$ via a convex Beckmann flow problem (second-order-cone programming), cross-checked against a dual LP, then build the map by direct integration or a ray-based 1D "clock" ODE with explicit density per ray. Validation uses strong duality, mass conservation, and vanishing-regularization limits, on a translation-invariant case and two 2D examples (a radial cone, a checkerboard) checked by symmetry; both show boundary-concentrated mass and non-unique flux where supports touch. We extend to 3D (theory unchanged), matching the 2D answer to six digits and reproducing the same phenomena on a radial ball, exposing a limit of the 2D visualization style. Movies illustrating the transport maps are available as ancillary files on the arXiv page of this submission.

1 Introduction

The paper completes an earlier numerical program by constructing Monge transport maps from the Evans–Gangbo potential and density framework. It addresses the missing particle travel distances in Kantorovich’s approach and develops numerical methods for the 2D map.

1 Introduction

  • The paper completes an earlier report’s unfinished 2D optimal-map construction, using structural results from Evans–Gangbo to build a more accurate, validated solver.The structural results include straight transport rays, an ODE for density along a ray, and a regularized transfer-plan construction.Source: §1 · p. 1

1.1 The Monge problem

  • Kantorovich duality identifies a potential for the optimal plan but does not directly specify how far each particle should move.The Evans–Gangbo PDE/ODE method supplies this missing distance information.Source: §1.1 · p. 2

1.2 The Evans–Gangbo construction

  • The Evans–Gangbo construction obtains the potential and transport density through a p-Laplacian limit, then uses them in an ODE to define the transfer map.The theory requires transport rays and regularized approximations because the potential, density, and gradient need not be continuous.Source: §1.2 · p. 2

1.3 Related numerical work

  • The paper’s stated numerical-literature gap is coupling the Beckmann flow to Evans–Gangbo map reconstruction and extending that framework to three dimensions.The authors position their contribution as complementary to existing numerical approaches.Source: §1.3 · p. 3

1.4 Structure of this paper

  • The paper proceeds from closed-form 1D examples to 2D construction and validation, further 2D tests, and then 3D examples.The 2D construction includes a convex reformulation, an independent linear program, direct ODE integration, and a ray-based method.Source: §1.4 · p. 3

2 Dimension one: closed-form solutions

In one dimension, the transport set is an interval and the transport ODE can be solved explicitly. The examples illustrate the potential, density, and resulting map for linear and cubic data.

2 Dimension one: closed-form solutions

  • On an interval transport set, the transport ODE (8) can be solved explicitly.The section recalls two examples on [−10, 10].Source: §2 · p. 3

2.1 Example: f(x) = x

  • For f(x) = x, the solution has u(x) = x and Du ≡ 1; its transport density is nonnegative and vanishes at x = ±10.The endpoint behavior agrees with the stated general property of transport density along a ray.Source: §2.1 · p. 3
  • The f(x) = x example’s optimal map is obtained by explicitly solving equation (8).The expression for the map is not included in the supplied passage.Source: §2.1 · p. 4
  • Symbolic differentiation verifies the f(x) = x map against equation (8) and resolves a sign ambiguity in the ODE.The density and map are also shown in Figures 1 and 2.Source: §2.1 · p. 4

2.2 Example: f(x) = x3

  • Figure 1 shows f and a(x) for f(x) = x and f(x) = x^3; in both cases, a vanishes at the transport-ray endpoints.The ray is [−10, 10].Source: Figure 1 · p. 4
  • Figure 2 shows reflection-like endpoint maps that move the positive part of f onto its negative part.Source: Figure 2 · p. 5

3 Dimension two: numerical construction

For f=x-y, the 2D construction recasts density recovery as Beckmann minimization, computes the potential independently, and builds maps by direct or ray-based ODE integration.

3.1 A convex reformulation of the transport-density equation

  • The transport-density system and Beckmann minimization are equivalent: either solution yields the other with the same optimal value, and a minimizing flux recovers a=|v|.The Beckmann problem is discretized as an SOCP on a staggered grid, where the discrete divergence theorem conserves mass to machine precision.Source: §3.1 · pp. 5–7

3.2 The Kantorovich potential via an independent linear program

  • 0.11779 is the objective of both the independent potential LP and Beckmann solve, matching to a relative gap of about 10^-8 and numerically confirming strong duality.The computed potential is nearly linear, with level sets parallel to y=x and transport rays perpendicular to that diagonal.Source: §3.1 · p. 5, §3.2 · p. 8, Figure 3 · p. 8

3.3 Direct integration of the transport ODE

  • Direct integration uses bilinear interpolation and vectorized explicit RK4; the resulting particles move along straight lines perpendicular to the diagonal.The reported implementation reduced the runtime for 700 particles over 300 steps from more than two minutes to under one second.Source: §3.3 · p. 10, Figure 6 · p. 10

3.4 A sign correction in equation (8)

  • Symbolic differentiation selects the minus sign in equation (8), resolving a transcription ambiguity and agreeing with the original Evans–Gangbo equation [2, eq. 1.23].Source: §3.4 · p. 10

3.5 Ray-based construction

  • Flux-derived directions avoid inconsistent discrete gradients; rays receive endpoint-constrained densities and a scalar clock ODE, a better-conditioned alternative to full vector integration.For this linear f and nearly linear u, the density along each ray is computed in closed form with a shooting correction that makes it vanish at both endpoints.Source: §3.5 · p. 11

4 Validation

Validation compares ray-based density with the Beckmann estimate, examines regularization and cross-checks two map constructions, then tests mass conservation. The results show close agreement in most checks, alongside a grid-related limitation in the regularization study.

4.1 Exact density along a ray vs. the grid estimate

  • About 10^-3 RMS difference separates the closed-form ray density from the grid-based Beckmann estimate across six representative rays.The compared density range is [0, 0.35] on a 44 × 44 discretization.Source: §4.1 · pp. 11–12, Figure 7 · p. 12

4.2 Regularization limit δ →0

  • At the smallest tested δ values, mean transport distance plateaus near 0.722, about 2% above target, rather than approaching it monotonically.The authors attribute the persistent bias to the fixed 44×44 grid; refining the spatial grid alongside δ is needed for a genuine convergence study.Source: §4.2 · p. 12, Figure 8 · p. 13

4.3 Cross-validation of the two constructions

  • The direct and ray/clock maps agree within about 10^-4 for most particles, while about 1% near the boundary differ by up to 0.4.The authors trace the boundary-localized disagreement to an artifact in the discretized Beckmann flux; Figure 9 shows the maps side by side.Source: §4.3 · p. 13, Figure 9 · p. 13

4.4 Mass conservation

  • Mass-conservation checks for four test functions have relative errors below 2.2% at the working resolution.The check uses a 44 × 44 grid and 500 Monte Carlo source points; Table 1 reports the errors.Source: §4.4 · p. 14, Table 1 · p. 14

5 Two additional two-dimensional test cases

The two additional 2D examples test curved transport directions under rotational and discrete symmetry. They expose limitations of the potential discretization and show boundary-focused routing where source and target supports touch.

5.1 A more accurate potential solve for curved solutions

  • The eight-neighbor Lipschitz stencil yields a 1.7% relative duality gap for the radial example, unchanged from 44 to 80 grid points per side.A rotationally invariant second-order-cone constraint reduces the gap from 1.97 × 10^-3 at 48 points to 4.9 × 10^-4 at 96.Source: §5.1 · pp. 14–15

5.2 A radial example: the cone f = c −r2

  • The radial example’s rotational symmetry supplies structural checks unavailable in the translation-invariant case, whose transport direction is constant.The radial potential is expected to be radially symmetric and its optimal flux radial; Figures 10 and 11 show the data and computed fields.Source: §5 · p. 14, §5.2 · p. 15, Figure 10 · p. 15
  • At the radial interface r = ρ0, touching source and target supports allow non-unique minimal-flux routing, while the potential remains comparatively robust.The flux’s angular deviation from the radial direction is about 6–11° and does not shrink under grid refinement.Source: §5.2 · p. 16
  • The idealized disk’s radial map is not a tight prediction for the square: relative errors range from 15% at small r to under 1% near r = ρ0.The square’s excavation mass is concentrated mainly in its four corners; the example’s Beckmann and potential-SOCP objectives are 0.014869 and 0.014886, respectively.Source: §5.2 · pp. 16–17, Figure 12 · p. 17

5.3 A checkerboard example: f = (x −1

  • For the checkerboard, the Beckmann and potential-SOCP objectives are 0.016657 and 0.016667, with a relative gap of 5.8 × 10^-4.The map also satisfies swap and central symmetry, and mass-conservation errors on four test functions range from 1.2% to 4.5%.Source: §5.3 · p. 18
  • The checkerboard flux hugs the square boundary, sending the northeast pile to the northwest excavation and the southwest pile to the southeast.The plotted routing avoids the center rather than mixing the two flows.Source: §5.3 · p. 18, Figure 14 · p. 18, Figure 15 · p. 19

6 Dimension three: a first extension

The paper extends the Evans–Gangbo construction to 3D without conceptual changes to the theory. A z-independent cube problem tests the implementation against the known 2D solution.

6 Dimension three: a first extension

  • The z-independent unit-cube test reduces to repeated copies of the known 2D problem, providing a check before genuinely 3D-varying data.The theory and Beckmann reformulation already apply in R^n.Source: §6 · p. 19

6.1 Discretization

  • The 3D discretization adds a z-face flux array and a three-component norm to the Beckmann and potential problems.On a 24 × 24 × 24 grid, the two SOCP solves take about 30 and 25 seconds, respectively, on one core.Source: §6.1 · p. 20

6.2 Validation

  • 0.707107 is the 3D solver’s mean transport distance, matching the known 2D answer to all six reported digits.The mean trajectory z-displacement is 1.9 × 10^-11; at N = 24, the duality gap is 6.0 × 10^-3 and shrinks under refinement.Source: §6.2 · pp. 20–21

7 A genuinely three-dimensional example: the radial ball

The 3D radial-ball example exhibits a larger, refinement-shrinking duality gap and nonradial flux. A ball-domain shell model captures the transport qualitatively, but not tightly, because the cube’s corner and edge region contains most of the excavation mass.

Example setup

  • The pile ball is tangent to all six cube faces, creating a more boundary-touching configuration than the 2D radial example and intensifying the touching-support phenomenon.The example is the genuinely 3D-varying analogue of the 2D radial cone.Source: §7 · p. 21

Strong duality and structural checks

  • The duality gap is 1.9 × 10−2 at N = 24, down from 2.7 × 10−2 at N = 20; flux directions deviate from radial by 17–28° on average.The gap exceeds those in the 2D radial example (1.1 × 10−3) and z-invariant 3D case (6.0 × 10−3); potential angular standard deviation grows from 0.003 at r = 0.06 to 0.046 at r = 0.46.Source: §7 · p. 22

An idealized semi-analytic comparison

  • The idealized ball-domain shell comparison has relative errors from 19% to under 4%, because the cube’s corner and edge region holds most of the excavation mass.The ball of radius ρ0 = 0.5 occupies π/6 ≈ 52% of the cube’s volume; matching shell masses gives a quintic in s, solved by 1D root-finding, with idealized radius L = 5c/3 ≈ 0.645.Source: §7 · p. 22

Transport summary

  • The total pile mass is 0.052131, and the mean transport distance is 0.353572.Mean distance is reported as objective divided by mass.Source: §7 · p. 23

Visualization

  • The central z-slice figure shows potential and in-slice flux with the pile boundary marked, while the final scatter plot shows target concentration toward cube corners and edges.The static source/target scatter plot is substantially clearer than the accompanying movie for this volumetric transport.Source: Figure 18 · p. 22, Figure 20 · p. 23, §7 · p. 24

8 Conclusion and future work

The paper completes its 2D construction and extends the same machinery to 3D, validating it on known and symmetric cases. The examples show boundary-concentrated transport and non-unique flux when supports touch, while motivating further numerical and visualization work.

Conclusion

  • The 2D construction completes the optimal Monge transfer plan left as future work in, using direct ODE integration and a ray-based method, with validation by duality, mass conservation, and regularization limits.For f(x, y) = x − y, the map moves perpendicular to the diagonal and decouples into one-dimensional problems.Source: §8 · p. 24
  • In the radial-cone and checkerboard examples, mass concentrates near the domain boundary; touching supports coincide with non-unique minimal flux despite unique potential and optimal cost.The radial case concentrates near corners, while the checkerboard concentrates along edges; touching supports violate hypothesis (2.4) of.Source: §8 · p. 24
  • The 3D method reproduces the known 2D mass-weighted mean transport distance to all six reported digits and confirms z-invariance, then shows boundary concentration and flux non-uniqueness in a radial-ball case.In 3D, the pile is tangent to the domain boundary; a partially transparent source/target scatter plot communicates the transport more clearly than the movie.Source: §8 · p. 24

Future work

  • Future work includes mesh-refinement studies, handling touching-support and boundary-layer discrepancies, extending the ray/clock method, developing scalable solvers, and improving volumetric rendering.The current ray/clock construction is worked out only for (13), and the generic SOCP method’s N^3 scaling makes finer 3D examples a motivation for more scalable solvers.Source: §8 · p. 25
Loading 2609.29987v1…