Large-scale semi-discrete optimal transport
with distributed Voronoi diagrams

Bruno Lévy¹, Nicolas Ray², Quentin Mérigot¹, Hugo Leclerc¹

1. Centre Inria Saclay, Université Paris Saclay, CNRS, Labo. de Maths. d'Orsay

2. Centre Inria de l'Université de Lorraine, CNRS, LORIA

bruno.levy@inria.fr

nicolas.ray@inria.fr

quentin.merigot@universite-paris-saclay.fr

hugo.leclerc@universite-paris-saclay.fr

January, 2025

Abstract

In this article, we propose a numerical method to solve semi-discrete
optimal transport problems for gigantic pointsets ($10^8$ points and more).
By pushing the limits by several orders of magnitude, it opens the path
to new applications in cosmology, fluid simulation and data science to
name but a few. The method is based on a new algorithm that com-
putes (generalized) Voronoi diagrams in parallel and in a distributed way.
First we make the simple observation that the cells defined by a subgraph
of the Delaunay graph contain the Voronoi cells, and that one can de-
duce the missing edges from the intersections between those cells. Based
on this observation, we introduce the Distributed Voronoi Diagram algo-
rithm (DVD) that can be used on a cluster and that exchanges vertices
between the nodes as need be. We also report early experimental re-
sults, demonstrating that the DVD algorithm has the potential to solve
some giga-scale semi-discrete optimal transport problems encountered in
computational cosmology.

# 1 Introduction and previous work

## 1.1 Optimal transport in computational physics

Optimal transport is the study of displacements that minimize a total integrated
transport cost — such as the squared Euclidean distance multiplies by the trans-
ported mass, while satisfying a mass conservation constraint. This problem was
initially studied by Gaspard Monge in 1784 [Mon84], but the first important
step towards the mathematical resolution of the transport problem is due to
Leonid Kantorovich in the 1940s. A quantum leap in understanding the structure of the problem was made in the 90s by Yann Brenier, with his celebrated polar factorization theorem [Bre91], that revealed rich connections with other fields of mathematics and physics. We refer the reader to [Vil03, Vil09, San15] and references herein. Another important outcome of optimal transport theory is the definition of transport distances, often referred to as "Wasserstein distances." These metrics enable comparisons between objects of different na-tures and dimensions, such as point clouds, surfaces, and probability densities. This versatility has been used in numerous recent developments in data science, artificial intelligence and inverse problems. Motivated by the wide spectrum of possible application, significant efforts have been dedicated over the past 20 years to the efficient numerical resolution of optimal transport, resulting in substantial progress [PC19, MT21, LS18].

In physics, optimal transport has interesting connections with the least ac-
tion principle [BB00] and gradient flows [JKO98]. These connections were ex-
plored to design new modelling, simulation and inverse methods in various fields
such as fluid simulation [GM18, dGWH<sup>+</sup>15, Lév22, QLdGJ22], cosmological
reconstruction [BFH<sup>+</sup>03, FMMS02, LMv21, vLM22, NSLM22, NPL<sup>+</sup>23], alter-
native theories of gravity [Bre15, LBM24], crowd dynamics [BC13, BCMO16,
LMSS20] and climatology [CP84, EBC<sup>+</sup>22, BCM24] to name but a few. How-
ever, most numerical implementation of these methods depend on an efficient
way of solving the equation that underlies the optimal transport problem, namely
the Monge-Ampère equation in the quadratic case. This equation can be viewed
as a generalization of the Poisson equation, but solving large-scale instances
poses significantly greater challenges. This difficulty is explained by the non-
linearity of the Monge-Ampère equation, but also by its degenerate ellipticity,
which translates into a convexity constraint that must be satisfied by the solu-
tion [BFO14]. The Fast Fourier transform had a tremendous impact by making
it computationally easy to solve the Poisson equation. Is it possible to find an
efficient and scalable computational method that would play the same role for
the Monge-Ampère equation and optimal transport ?

A possible candidate is the “entropic regularization” idea [Cut13], that uses
a softened iterative projection method in dual space, leading to spectacularly
efficient solvers for a fixed value of the regularization parameter, which opened
the way to many applications of optimal transport in machine learning and data
science [PC19]. To enable scaling up to problems of very large size, recent works
also explore ways of distributing computations, such as multiple domain decom-
position [MS24, MFF<sup>+</sup>17], and combination of clever multiresolution techniques
(on the point sets) and scaling techniques (on the regularization parameter) to
decrease the number of iterations [FSV<sup>+</sup>19]. However, despite these improve-
ments, the complexity of each iteration of the Sinkhorn algorithm that is used for
solving entropy-regularized optimal transport between unstructured point cloud
is quadratic in the size of the problem, prohibiting very large-scale applications.

In the application domain targeted by the present work, namely computa-
tional cosmology, one needs to solve optimal transport problems that are very
singular and where the source and target measures are of different nature. The
source measure represents the initial condition of the Universe, considered to be uniform<fnref n="1" /> and modeled as the Lebesgue measure over some domain. The target measure, on the other hand, is a finite set, and may contain up to $10^8$ or $10^{10}$ weighted Dirac masses. Each weighted Dirac mass corresponds to a galaxy cluster, modeled as a punctual mass, which makes sense at the cosmological scale. In these datasets, there are huge variations of density, typically 5 orders of magnitudes or more: there are zones where matter clustered a lot, creating highly dense clusters of galaxies, and other areas that are nearly completely empty of matter. This is illustrated in Figure 4, p. 21). This makes entropy-regularized methods difficult to use, because their efficient implementation depend on regular grids. The cell size is conditioned by the shortest distance between two Dirac masses, which would require in our case extremely fine grid resolution, making space requirement and computation time prohibitive. It may be possible to develop a version working with adaptively subdivided grids, but it is out of reach of current solvers.

To keep the accuracy required to represent these highly heterogeneous distributions of matter, it is possible to use so-called *semi-discrete* optimal transport [Mér11, KMT19, Lév15]. Semi-discrete optimal transport exploits the generality of optimal transport theory, that is written in a general mathematical language (measure theory), that deals not only functions (probability densities), but also much more irregular objects, such as weighted point sets, modeled as measures supported on finite sets. Remarkably, the very same mathematical language can be used from the mathematical modeling to computer implementation, without having mathematical details “lost in translation”, as summarized in the next subsection.

## 1.2 From Monge's problem to semi-discrete OT

**Monge's problem** We consider our domain $X$ to be a compact subset of the Euclidean space $\mathbb{R}^d$, even though the results that we present below hold on more general spaces, such as tori, which allow to account for periodic boundary conditions. On this set $X$, we consider two non-negative measures $\mu$ and $\nu$ with the same total mass, i.e. $\int_X d\mu = \int_X d\nu$. These measure describe a distribution of some mass in $X$, which might be diffuse when $\mu$ has a density or highly concentrated when, for instance, $\mu$ consists in one or several Dirac masses. We call *transport map between $\mu$ and $\nu$* any map $T: X \to X$ which deforms the mass distribution distributed by $\mu$ onto $\nu$. More formally, this means that for any measureable subset $B$ of $X$, we want amount of mass of $\nu$ contained in $B$ to be equal to the amount of mass of $\mu$ contained in the pre-image $T^{-1}(B)$, i.e. $\mu(T^{-1}(B)) = \nu(B)$. If $T$ is a transport map between $\mu$ and $\nu$, we denote $T_{\#}\mu = \nu$. The quadratic version of Monge's problem [Mon84] consists in finding

<footnote n="1">¹it is supposed that the Universe had a uniform density at its early ages, which seems to be confirmed by the cosmic microwave background that is very homogeneous</footnote>
a transport map $T: X \to X$ that minimizes the transport cost:

$$ \inf_{T:X \to X|T_{\#}\mu=\nu} \int_X \|T(\mathbf{x}) - \mathbf{x}\|^2 d\mu(\mathbf{x}) \quad (1) $$

We note at this point that this problem can be very ill-posed, and in particular that the set of admissible transport maps might be empty.

**Kantorovich's relaxation** Most numerical methods for solving optimal transport rely on (the dual of) a relaxed version of Monge's problem, which was proposed by Kantorovich in the 1940s. The relaxation is given by

$$ \inf_{\gamma \in \Gamma(\mu, \nu)} \int_{X \times X} \|\mathbf{x} - \mathbf{y}\|^2 d\gamma(\mathbf{x}, \mathbf{y}) \quad (2) $$

where the set $\Gamma(\mu, \nu)$ is the set of non-negative measures over $X \times X$ with marginals $\mu$ and $\nu$, that is, satisfying the constraints

$$ \Pi_1\#\gamma = \mu, \quad \Pi_2\#\gamma = \nu $$

where $\Pi_i : (x_1, x_2) \in X \times X \to x_i \in X$ are the projection maps. A measure $\gamma$ that satisfies these two marginal constraints is called a *transport plan* between $\mu$ and $\nu$. It corresponds to the notion of coupling between the measures $\mu$ and $\nu$ in probability theory.

The measure $\gamma$ is supported on the Cartesian product $X \times X$ and may be thought of as a possibly diffuse approximation of the graph of a transport map. The quantity $\gamma(A \times B)$ is the amount of mass that is transported from the set $A \subset X$ to the set $B \subset X$ by the transport plan $\gamma$. The marginal constraint $\Pi_1\#\gamma = \mu$ implies that $\gamma(A \times X) = \mu(A)$ and encodes the fact that all the mass emanating from the set $A$ has to be sent by $\gamma$ to some point in $X$. The constraint $\Pi_2\#\gamma = \nu$ is interpreted similarly.

In the relaxed problem (2) the objective function and the constraints are both linear in the unknown transport plan $\gamma$, i.e. Kantorovich's relaxation is an infinite-dimensional linear programming problem.

**Kantorovich dual and c-transform** The linear and linearly constrained nature of Kantorovich's problem calls for taking a look at the dual formulation:

$$ \sup_{\phi, \psi} \left[ \int_X \phi(\mathbf{x}) d\mu + \int_X \psi(\mathbf{y}) d\nu \right] \quad (3) \\ \text{subject to: } \phi(\mathbf{x}) + \psi(\mathbf{y}) \le \|\mathbf{x} - \mathbf{y}\|^2 \quad \forall \mathbf{x}, \mathbf{y} \in X \times X $$

where $\phi$ and $\psi$, the Lagrange multipliers associated with the marginal constraints, are both functions from $X$ to $\mathbb{R}$. It is interesting to notice that as compared to the primal problem, that depends on an unknown transport plan $\gamma$ on $X \times X$, the dual problem depends on two unknown functions $\phi, \psi: X \to \mathbb{R}$. Assuming that the measures $\mu$ and $\nu$ are probability densities, it is known since the work of Brenier [Bre91] that Monge's problem (1) admits a solution $T$. The
optimal transport map $T$ from $\mu$ to $\nu$ and its inverse $T^{-1}$ are connected to the maximizers of the dual Kantorovich problem through the following formulas:

$$ 
\begin{aligned}
T(\mathbf{x}) &= \mathbf{x} - \frac{1}{2}\nabla\phi(\mathbf{x}) \\
T^{-1}(\mathbf{y}) &= \mathbf{y} - \frac{1}{2}\nabla\psi(\mathbf{y})
\end{aligned}
\qquad (4) $$

With computer implementation in mind, considering the dual problem seems to offer significant advantages: it appears much easier to work with two functions, $\phi, \psi: X \to \mathbb{R}$, rather than a measure $\gamma$ supported on $X \times X$. If the space $X$ is discretized using $N$ points, the dual problem approach involves $2N$ variables, each function being represented as a vector with $N$ entries, compared to the $N^2$ variables potentially required to discretize the transport plan $\gamma$. However, in this formulation of the dual problem, the number of constraints becomes quadratic. This challenge can be addressed using the concept of the $c$-transform, a variant of the concept of convex conjugation or Legendre-Fenchel transform. Specifically,<fnref n="2" /> if $\psi: X \to \mathbb{R}$ is given, the function $\phi$ that maximizes the dual cost $\int \psi d\mu$ while satisfying the constraints $\phi(\mathbf{x}) + \psi(\mathbf{y}) \le \|\mathbf{x} - \mathbf{y}\|^2$ can be constructed easily: one should take $\phi = \psi^c$, where

$$ \psi^c(\mathbf{x}) = \inf_{\mathbf{y} \in X} \|\mathbf{x} - \mathbf{y}\|^2 - \psi(\mathbf{y}). \qquad (5) $$

Using this construction of the $c$-transform, the dual problem becomes the following problem, that depends on a single function $\psi: X \to \mathbb{R}$:

$$ \sup_{\psi} \mathcal{K}(\psi) = \int_X \psi^c(\mathbf{x}) d\mu + \int_X \psi(\mathbf{y}) d\nu \qquad (6) $$

From an implementation viewpoint, the main advantage of this formulation is that we have turned the (constrained) linear programming problem (3) into an unconstrained and concave maximization problem. We have however, moved the difficulty induced by constraints in (3) into the computation of the $c$-transform $\psi^c$ in the unconstrained problem (6).

**Semi-discrete optimal transport** Up to now we did not make any assumption regarding the measures $\mu$ and $\nu$ besides the assumptions that they have the same total mass. Semi-discrete optimal transport is a special case, where $\mu$ is a density over the set $X$ and $\nu = \sum_{i=1}^N \nu_i \delta_{x_i}$ is a finitely supported measure, i.e. a weighted sum of Dirac masses, satisfying the mass balance $\sum_i \nu_i = \int_X d\mu$. We also assume that the points $x_i$ are distinct. In this particular setting, the Kantorovich dual (6) can be further simplified. Indeed, a function $\psi$ on the finite support of the measure $\mu$ can be fully described by its values at the points

<footnote n="2">²The problem is symmetric, one could swap the roles of $\phi$ and $\psi$. Here we consider that the second Lagrange multiplier $\psi$ is given and that the first one $\phi$ is deduced as the $c$-transform of $\phi$, in that order. As we shall see further on, in the semi-discrete setting, $\psi$ has a finite support and can be represented as a vector of $N$ values.</footnote>
$\mathbf{x}_i$, denoted $\psi_i$, thus by a vector of $N$ scalar values $(\psi_i)_{1\le i\le N}$. We therefore consider the function $K: \mathbb{R}^N \to \mathbb{R}$ defined by

$$ 
\begin{aligned}
K(\psi) &= \int_X \psi^c(\mathbf{x})d\mu(\mathbf{x}) + \sum_{i=1}^N \psi_i \nu_i \\
&= \int_X \min_i [\|\mathbf{x} - \mathbf{x}_i\|^2 - \psi_i] d\mu(\mathbf{x}) + \sum_i \psi_i \nu_i
\end{aligned}
$$

Since the minimum inside the integral is taken over a finite set, it is very natural to introduce the minimization diagram of the functions $\|\mathbf{x} - \mathbf{x}_i\|^2 - \psi_i$, which is often called a Laguerre diagram in computational geometry. More precisely, we define the ith Laguerre cell to be

$$ \text{Lag}_i^\psi = \{ \mathbf{x} \in \mathbb{R}^d \mid \|\mathbf{x} - \mathbf{x}_i\|^2 - \psi_i \le \|\mathbf{x} - \mathbf{x}_j\|^2 - \psi_j \quad \forall j \ne i \}. $$

With this definition at hand, one might replace the integral over $X$ in (7) into a sum of integral over Laguerre cells:

$$ K(\psi) = \sum_i \int_{\text{Lag}_i^\psi} (\|\mathbf{x} - \mathbf{x}_i\|^2 - \psi_i) d\mathbf{x} + \sum_i \psi_i \nu_i \quad (7) $$

In the setting that we just described, the Kantorovich dual $K$ depends on the vector $(\psi_i)_{1\le i\le N}$. The Kantorovich potential $\phi$, on the other hand, is continuous and defined on the whole domain $X$. However, it is defined as the c-transform $\psi^c$ of $\psi$, and is therefore also naturally parameterized by the same vector. We summarize some of the known results about this functional and its relation to semi-discrete optimal transport into a theorem, which follows from [AHA98, §5], and whose proof can be found in [MT21, Theorem 40].

**Theorem 1.** Assume that all the points $\mathbf{x}_1, \dots, \mathbf{x}_N$ are in distinct. The function $K: \mathbb{R}^N \to \mathbb{R}$ is concave, continuously differentiable, and its gradient is given by

$$ \frac{\partial K}{\partial \psi_i} = \nu_i - \int_{\text{Lag}_i^\psi} d\mu. \quad (8) $$

Moreover, $K$ admits at least a global maximizer. If $\psi$ is this maximizer and $\phi(\mathbf{x}) = \psi^c(\mathbf{x}) = \min_i \|\mathbf{x} - \mathbf{x}_i\|^2 - \psi_i$ its c-transform, the optimal transport map between $\mu$ and $\nu$ is given by $T: \mathbf{x} \mapsto \mathbf{x} - \frac{1}{2}\nabla\phi(\mathbf{x})$. This map $T$ is piecewise constant, equal to $\mathbf{x}_i$ on the ith Laguerre cell $\text{Lag}_i^\psi$.

Thanks to this theorem, we could think of solving the dual of the Kantorovich problem using gradient ascent of the functional $K$, as suggested in [AHA98, §5]. However, it turns out that because of the special structure of the functional $K$, evaluating $K(\psi)$ or $\nabla K(\psi)$ is as costly as evaluating the second derivative $D^2K(\psi)$. We will therefore rely on a second order Newton method for the maximization of $K$ rather than a first-order method. To do so, we will restrict
the set of potentials $\psi$ that are admissible in the iterations of the algorithm.
We will denote

$$ \mathrm{Adm}_{\varepsilon} = \left\{ \psi \in \mathbb{R}^N \mid \forall i, \int_{\mathrm{Lag}_i^{\psi}} d\mu \geq \varepsilon \right\}, \quad \mathrm{Adm} = \bigcup_{\varepsilon>0} \mathrm{Adm}_{\varepsilon}. $$

Thus, $\mathrm{Adm}$ is the set of potentials $\psi$ such that all Laguerre cells contains a positive fraction of the mass of $\mu$. Under some regularity assumptions on $\mu$ and $X$ one can show that $K$ is twice differentiable and give an explicit expression for its second derivatives. The following theorem is a special case of [KMT19], see also [Lév15]. It shows that the functional $K$ is twice continuously differentiable on the admissibility set $\mathrm{Adm}$:

**Theorem 2.** Assume that $\mu$ has a continuous density on a convex compact domain $X$, which is bounded from below by a positive constant. Then, $K$ is twice continuously differentiable on $\mathrm{Adm}$ and its second partial derivatives are

$$ \begin{align} \frac{\partial^2 K}{\partial \psi_i \partial \psi_j}(\psi) &= \frac{1}{2} \frac{1}{\|\mathbf{x}_j - \mathbf{x}_j\|} \int_{\mathrm{Lag}_{i,j}^{\psi}} \mu(x) d\mathrm{vol}^{d-1}(x) \quad \text{if } j \neq i \tag{9} \\ \frac{\partial^2 K}{\partial \psi_i^2}(\psi) &= -\sum_{j \neq i} \frac{\partial^2 K}{\partial \psi_i \partial \psi_j}(\psi) \nonumber \end{align} $$

where $\mathrm{Lag}_{i,j}^{\psi}$ denotes the $(d-1)$-dimensional intersection between the $i$th and $j$th Laguerre cells, and where the integral in the first formula should be interpreted as an integral with respect to the surface measure.

We note that the formula (9) is reminiscent of a finite-volume discretization of the Poisson equation on the mesh defined by the Laguerre cells [EGH00, §3], which may not be surprising as the equation $\nabla K(\psi) = 0$ is a discretized Monge-Ampère equation, as explained in [MT21, Remark 18]. From formula (9), we can directly see that the matrix $D^2K(\psi)$ of second-order partial derivatives is weakly diagonally dominant with negative diagonal entries, and therefore has non-positive eigenvalues. We could also have derived this property from the concavity of $K$. Exploiting in addition the connectedness of the graph over the set $\{1, \dots, N\}$, where two vertices $i \neq j$ are connected if $D^2K(\psi)_{i,j} > 0$, we can verify that the kernel of the matrix $D^2K(\psi)$ is spanned by the constant vector $(1, \dots, 1)$. This shows that the restriction of $K$ on $\mathrm{Adm} \cap \{\psi \in \mathbb{R}^N \mid \psi_0 = 0\}$ is strongly concave$^3$.

At this point, we have done nothing else than rewriting the generic optimal transport problem in the specific case where $\mu$ is the Lebesgue measure and $\nu$ is an empirical measure. Let us now explain how to compute the optimal transport in practice, that is, how to determine the vector $(\psi_i)_{1 \le i \le N}$. The

<footnote>$^3$where the constraint $\psi_0 = 0$ (for instance) determines the translational degree of freedom along the $(1, \dots, 1)$ vector. Intuitively, the vector of Lagrange multipliers $\psi$ may be thought of as a *potential*, defined up to a translational constant $(C, \dots, C)$.</footnote>
Kantorovich dual that depends on these parameters is concave, twice differentiable and strongly convex on Adm. Moreover, the solution $\psi$ is characterized by $\int_{\text{Lag}_i^\psi} d\mu = \nu_i > 0$, and therefore also belongs to this set Adm. The Newton method that we will use will rely on a specific backtracking which ensures that its iterates remain in the set Adm, as explained in Algorithm 2. More specifically, the algorithm iteratively halves the descent parameter $\alpha$ until two criteria are met: the volume of the smallest cell needs to be larger than a threshold $a_0$ (first condition in line 3), and the norm of the gradient needs to decrease sufficiently (second condition in line 3). The threshold $a_0$ for the minimum cell volume corresponds to (half) the minimum cell volume for $\psi = 0$ (also called Voronoi diagram) and minimum prescribed area $\nu_i$.

**Algorithm 1. Semi-discrete optimal transport**

**input:** a pointset $(\mathbf{x}_i)_{i=1}^N$ and masses $(\nu_i)_{i=1}^N$
**output:** the Laguerre diagram $\{\text{Lag}_i^\psi\}_{i=1}^N$ such that $|\text{Lag}_i^\psi| = \nu_i \quad \forall i$

(0) $\psi \leftarrow 0$

(1) **while** $\|\nabla K\|_\infty < \epsilon$

(2) **solve for p** in $[\nabla^2 K(\psi)]p = -\nabla K(\psi)$

(3) find descent parameter $\alpha$ (see Algorithm 2)

(4) $\psi \leftarrow \psi + \alpha p$

(5) **end while**

**Algorithm 2. Kitagawa-Mérigot-Thibert descent (KMT)**

**input:** current values of $(\psi_i)_{i=1}^N$ and Newton direction $p$

**output:** descent parameter $\alpha$ determining the next iterate $\psi \leftarrow \psi + \alpha p$

(1) $\alpha \leftarrow 1$

(2) **loop**

(3) **if** $\inf_i |\text{Lag}_i^{\psi+\alpha p}| > a_0$ **and** $\|\nabla K(\psi + \alpha p)\| \le (1-\alpha/2)\|\nabla K(\psi)\|$

(4) **then exit loop**

(5) $\alpha \leftarrow \alpha/2$

(6) Compute the Laguerre diagram $(\text{Lag}_i^{\psi+\alpha p})_{i=1}^N$

(7) **end loop**

where $a_0 = \frac{1}{2} \min (\inf_i |\text{Lag}_i^{\psi=0}|, \inf_i \nu_i)$.

In the numerical algorithm for semi-discrete optimal transport above, the two dominant costly operations are computing Laguerre diagrams and solving linear systems. By examining the coefficients of the linear system given in Equation 9, we have already noticed a strong similarity with a Poisson system. Like in a Poisson system, the Hessian $\nabla^2 K$ is semi-definite positive. It has the
vector $[1, 1... 1]^t$ in its kernel. In such a configuration, the conjugate gradient algorithm converges to the unique solution of the system orthogonal to the kernel, which is also the minimum-norm one (see e.g. [Hay20]). To accelerate the convergence of the conjugate gradient algorithm, an algebraic multigrid preconditioner [Dem19] is highly effective, as we noticed in [LBM24] and confirmed experimentally. In this setting, the lion's share of the computational cost is the computation of Laguerre diagrams: for a typical reconstruction problem in cosmology, it takes 60% of the computation time (whereas matrix assembly and linear solve take 25% and 15% respectively). Besides computational cost, another important obstacle to scaling up is the memory cost of the data structure used to store the Laguerre diagrams. Both aspects call for a distributed algorithm, that could run on a cluster. We shall now review the standard algorithm to compute (generalized) Voronoi diagrams, then we will see how to derive a distributed version, that can run on a cluster.

## 1.3 Voronoi diagrams and algorithms to compute them

Voronoi diagrams are ubiquitously used in various domains of computational sciences, including geometry processing, data analysis, machine learning and computational physics to name but a few. The properties of Voronoi diagrams and their various generalizations were studied [Aur91, EM90] and used to derive efficient computational algorithm, see surveys in [BY98] and [OBSC00]. To efficiently represent a Voronoi diagram in a computer, one can use the dual combinatorial structure, called the Delaunay triangulation. Historically, one of the first methods, proposed in the 80s by Sloan [Slo87], transforms an arbitrary triangulation into the Delaunay triangulation by flipping the edges until all the triangles are valid, that is, such that their circumcircles are empty of any other point. The most widely used algorithm to compute Voronoi diagram was proposed in the 80s, simultaneously by Bowyer and Watson [Bow81, Wat81]. It also operates on the Delaunay triangulation. It inserts the points one by one in the triangulation, and iteratively discards the triangles that violate the empty circumcircle condition and creates new ones, connecting the newly inserted points to the boundary of the zone where triangles were discarded, which is more efficient than flipping edges.

Numerical precision is an important aspect to be considered with uttermost care when implementing these algorithms. The most important and subtle part of a geometric algorithm is a set of functions, called geometric predicates, that take as argument geometric objects (for instance two points $x_1$ and $x_2$) and that returns a discrete value (for instance, +1 if $x_1$ is above $x_2$, -1 if $x_1$ is below $x_2$ and 0 if they have the same altitude). In general, these geometric predicates are more complicated, and correspond to the sign of a polynomial in the points coordinates. The main difficulty is that the floating point numbers manipulated in computers have limited precision. If no care is taken, the computed signs are not exact, and depending on the order of the operations, the predicate may indicate that $x_1$ is above $x_2$ in one part of the algorithm, and that
$x_1$ is below $x_2$ in another part, leading to inconsistencies in the constructed combinatorial structure. For this reason, different techniques were proposed to compute the exact sign of the polynom, and several software packages were proposed, by Shewchuk [She97] (who uses arithmetic expansions), in CGAL [The23] (that uses a combination of arithmetic filters, interval arithmetic and exact numbers), PCK [Lév16] and GEOGRAM [Lmc23], that use a combination of arithmetic filters and arithmetic expansions.

To speed-up the computation of Voronoi diagrams, a parallel algorithm based on the “security radius” theorem was proposed [LB12], as well as a version that runs on GPUs for Voronoi and Laguerre diagrams [RSLL18, BAR$^+21$. To compute Delaunay triangulations and Voronoi diagram with gigantic pointsets, several techniques were proposed, based on decomposing the input domain into regions. The “streaming” algorithm proposed in [ILSS06] inputs a stream of points and outputs a stream of triangles. It minimizes the globally stored information based on “finalization” tags inserted in the point stream, indicating that a point is the last one in its region, which allows one to discard the information that is no-longer needed. Another approach consists in first computing the Delaunay triangulation in each region independently, then determine the points that should be exchanged with the neighboring regions in order to get the correct triangulation / Voronoi cells [PMP14]. We mention also the AREPO code [Spr10], used for large scale cosmological simulations, that is based on a variant of Sloan’s edge flipping algorithm, distributed over a cluster, using an algorithm similar to the one in [PMP14]. As done in [CMYB19], a similar result can be obtained by using the *star splaying* algorithm [She05], that can restore the Delaunay condition in a nearly-Delaunay triangulation. It represents a mesh as a collection of (possibly incoherent) vertex stars, and restores coherency between them.

As explained later, the method proposed in this article shares some similarities with the *stars splaying* idea. The main difference is the path of reasoning that we followed, that is nearly equivalent, but that is structured around the Voronoi cells (instead of the 1-ring neighborhoods). This path of reasoning naturally leads to several algorithms for parallel and distributed Voronoi diagrams. As far as the distributed algorithm is concerned, the approach results in a smaller number of points exchanged between the regions.

With the goal of making it possible for semi-discrete optimal transport (see Section 3) to scale up beyond $10^8$ particles, what follows introduces new algorithms to compute (generalized) Voronoi diagrams in parallel and in a distributed way. The idea is based on the simple observation that the cells defined by a subgraph of the Delaunay graph contain the Voronoi cells, and that their intersections define the edges of the Delaunay graph (Section 2.1). From this observation, we deduce two algorithms, a parallel one (PVD for Parallel Voronoi Diagram), described in Section 2.2, that can be used in multicore shared memory machines, and a distributed one (DVD for Distributed Voronoi Diagram),
<figure><img src="image_1.png" /><figcaption>Figure 1: Left: the intersection of four half-spaces $\Pi^+(i, j)$ (their boundaries are highlighted in red). The area shaded in blue corresponds to the excluded half-spaces. The white region corresponds to the intersection between four half-spaces. It contains the Voronoi cell. Center: The Voronoi cell $Vor(\mathbf{x}_i)$ corresponds to the intersection of all the half-spaces $\Pi^+(i, j)$. Among them, some touch the boundary of $Vor(\mathbf{x}_i)$ and are contributing (in red), and some are non-contributing (in blue). Right: the Delaunay graph (in red) corresponds to the set of edges $(i - j)$ such that $\Pi^+(i, j)$ is contributing to $Vor(\mathbf{x}_i)$ (or equivalently, $\Pi^+(j, i)$ is contributing to $Vor(\mathbf{x}_j)$).</figcaption></figure>

described in Section 2.3, that computes a (potentially gigantic) Voronoi diagram on a cluster, and exchanges vertices between the nodes as need be. The distributed Voronoi Diagram algorithm (DVD) can be used to design a giga-scale semi-discrete optimal transport algorithm (see Section 3 for early results).

# 2 Distributed Voronoi Diagram

## 2.1 The Voronoi diagram

**Definition 1** (Voronoi diagram). The generalized Voronoi diagram $Vor(F)$, also called lower envelope or minimization diagram, associated to a set of $N$ continuous functions $F_1, \dots, F_N$ on $\mathbb{R}^d$ is the set of generalized Voronoi cells defined by

$$Vor_i(F) = \{\mathbf{x} \in \mathbb{R}^d | \forall j \neq i, F_i(\mathbf{x}) \leq F_j(\mathbf{x})\}.$$

The $i$th generalized Voronoi cells can be recovered by intersecting the $N - 1$ dominance regions $\Pi_{i,j}^+(F)$ for $j \neq i$ where

$$\Pi_{i,j}^+(F) = \{\mathbf{x} \in \mathbb{R}^d | F_i(\mathbf{x}) \leq F_j(\mathbf{x})\}. \qquad (10)$$

Finally we call bisector the set of points $\Pi(i, j)(F) = \{\mathbf{x} | F_i(\mathbf{x}) = F_j(\mathbf{x})\}$.

If $\mathbf{x}_1, \dots, \mathbf{x}_N$ are $N$ distinct points in $\mathbb{R}^d$ (called generating points or sites), and if $F_i(\mathbf{x}) = \|\mathbf{x} - \mathbf{x}_i\|^2$, then the generalized Voronoi cell $Vor_i(F)$ coincides with the usual Voronoi cell associated with the point $\mathbf{x}_i$. If in addition we are given $N$ scalars $\psi_1, \dots, \psi_N$ and if $F_i(\mathbf{x}) = \|\mathbf{x} - \mathbf{x}_i\|^2 - \psi_i$, then the minimization
diagram is called a *Laguerre* or *power* diagram. In both cases, the dominance regions are halfspaces and the general Voronoi cells are convex polyhedra. These notions are illustrated in Figure 1.

**Definition 2** (Voronoi diagram associated to a graph). We call a graph over {$1, \dots, N$} any subset $\mathcal{E}$ of the product space {$1, \dots, N$}^2. An edge $i \to j$ belongs to the graph if and only if the pair $(i, j)$ belongs to $\mathcal{E}$. The Voronoi cell associated to the graph $\mathcal{E}$ is defined by:

$$ \mathrm{Vor}_i^{\mathcal{E}}(F) = \bigcap_{j \text{ s.t. } (i,j) \in \mathcal{E}} \Pi_{i,j}^{+}(F), \qquad (11) $$

and the Voronoi diagram $\mathrm{Vor}^{\mathcal{E}}(F)$ associated to the graph $\mathcal{E}$ is the set of cells $\mathrm{Vor}_i^{\mathcal{E}}(F)$.

One can make the following (trivial yet useful) observation:

**Observation 1.** Given an arbitrary set $\mathcal{E}$ of oriented edges, the Voronoi regions with respect to the graph $\mathcal{E}$ contain the corresponding Voronoi regions:

$$ \mathrm{Vor}_i(F) \subseteq \mathrm{Vor}_i^{\mathcal{E}}(F) $$

In general, the inclusion is strict, and the regions $\mathrm{Vor}_i^{\mathcal{E}}(F)$ may overlap (not only over their boundaries).

*Proof.* This property is directly deduced from the definition of $\mathrm{Vor}_i(F)$ and $\mathrm{Vor}_i^{\mathcal{E}}(F)$:

$$ \mathrm{Vor}_i(F) = \bigcap_{j \neq i} \Pi_{i,j}^{+}(F) \subseteq \bigcap_{j | (i,j) \in \mathcal{E}} \Pi_{i,j}^{+}(F). \qquad \square $$

In terms of these definitions, starting from an arbitrary oriented graph $\mathcal{E}$, our goal is to find a new graph $\mathcal{E}'$ that contains all the edges that are missing in $\mathcal{E}$ for $\mathrm{Vor}^{\mathcal{E}'}(F)$ to be the Voronoi diagram $\mathrm{Vor}(F)$. To do so, we make another simple observation:

**Observation 2.** Let $\mathcal{E}$ be a graph over {$1, \dots, N$}, and define

$$ \mathcal{E}' = \{(i, j) \in \{1, \dots, N\}^2 \mid \mathrm{Vor}_i^{\mathcal{E}}(F) \cap \mathrm{Vor}_j^{\mathcal{E}}(F) \neq \emptyset\}, $$

Then,

$$ \mathrm{Vor}^{\mathcal{E}'}(F) = \mathrm{Vor}(F). $$

*Proof.* We already know that the cell $\mathrm{Vor}_i^{\mathcal{E}'}(F)$ relative to the graph $\mathcal{E}'$ contains the Voronoi cell $\mathrm{Vor}_i(F)$. We now need to prove inclusion in reverse order. Consider a point $\mathbf{x}$ in $\mathrm{Vor}_i^{\mathcal{E}'}(F)$. We shall show that $F_i(\mathbf{x})$ is smaller than $F_j(\mathbf{x})$ for all $j$'s:

* For an edge $(i, j)$ of the graph $\mathcal{E}'$ we have, by definition of $\mathrm{Vor}_i^{\mathcal{E}'}(F)$:

$$ \forall j \in \{1, \dots, N\}, (i, j) \in \mathcal{E}' \Rightarrow F_i(\mathbf{x}) \le F_j(\mathbf{x}). $$
* Consider now an edge $(i, j)$ that is not in $\mathcal{E}'$. By definition of $\mathcal{E}'$, we have $\mathcal{V}or_i^{\mathcal{E}'}(F) \cap \mathcal{V}or_j^{\mathcal{E}'}(F) = \emptyset$, thus $\mathcal{V}or_i^{\mathcal{E}'}(F) \cap \mathcal{V}or_j(F) = \emptyset$ (because $\mathcal{V}or_j(F) \subseteq \mathcal{V}or_j^{\mathcal{E}'}(F)$), hence there is no Voronoi cell $\mathcal{V}or_j(F)$ that contains $\mathbf{x}$:

$$ \forall j \in \{1, \dots, N\} \setminus \mathcal{N}_i^{\mathcal{E}'}, \mathbf{x} \notin \mathcal{V}or_j(F) $$

where: $\mathcal{N}_i^{\mathcal{E}'} = \{i\} \cup \{j \mid (i, j) \in \mathcal{E}'\}$.

Remembering that the Voronoi diagram is a covering of $\mathbb{R}^d$, there exists a Voronoi cell $\mathcal{V}or_k(F)$ that contains $\mathbf{x}$, with $k$ in the complement of $\{1, \dots, N\} \setminus \mathcal{N}_i^{\mathcal{E}'}$, that is, $k \in \mathcal{N}_i^{\mathcal{E}'}$. Hence, $F_j(\mathbf{x}) \ge \inf_{k \in \mathcal{N}_i^{\mathcal{E}'} \setminus \{k\}} F_k(\mathbf{x})$.

Remembering that $\mathbf{x}$ is in $\mathcal{V}or_i^{\mathcal{E}'}(F)$ and by definition of $\mathcal{V}or_i^{\mathcal{E}'}(F)$, we have $\inf_{k \in \mathcal{N}_i^{\mathcal{E}'} \setminus \{k\}} F_k(\mathbf{x}) = F_i(\mathbf{x})$ hence $F_j(\mathbf{x}) \ge F_i(\mathbf{x})$.

□

## 2.2 PVD: Parallel Voronoi Diagram

From Observation 2, one can deduce a new algorithm to compute a Voronoi diagram:

**Algorithm 3. Parallel Voronoi Diagram (PVD)**

input: a pointset $\mathbf{X}$
output: the Voronoi diagram $\{\mathcal{V}or_i\}_{i=1}^N$

(1) $\mathcal{E} \leftarrow kNN(\mathbf{X})$

(2) compute the cells $(\mathcal{V}_i^{\mathcal{E}})_{i=1}^N$

(3) $\mathcal{E} \leftarrow \mathcal{E} \cup \{(i-j) \mid \mathcal{V}_i^{\mathcal{E}} \cap \mathcal{V}_j^{\mathcal{E}} \neq \emptyset\}$

(4) re-compute the cells $(\mathcal{V}_i^{\mathcal{E}})_{i=1}^N$

where:

* in step (1), the k-nearest neighbor graph (kNN) is extracted from a kD-tree (Figure 2-A). Both the kD-tree construction and extraction of the k-nearest neighbors is computed in parallel. The reader is referred to [Sam90] for a survey on spatial data structures (including kNNs, KD-trees, and axis-aligned bounding box trees used later in step 3). The number $k$ of neighbors in the kNN is chosen based on the average number of neighbors in a Delaunay triangulation (6 in 2D, 16 in 3D);
<figure><img src="image_2.png" /><figcaption>Figure 2: A: k-NN graph; B: the cells deduced from the k-NN graph. They mostly correspond to the Voronoi diagram, except some intersections (in red); C: detecting the intersections with an Axis-Aligned Bounding Box tree; D: adding an edge in the graph for each pair of intersecting cells; E: the cells now correspond to Voronoi cells; F: the Delaunay graph contains the edges that correspond to the contributing bisectors. One gets the Voronoi cells as soon as the graph contains all these edges (as is the case at step D).</figcaption></figure>
* in step (2), the cells $V_i^\mathcal{E} = \bigcap_{(i \to j) \in \mathcal{E}} \Pi^+(i, j)$ are computed in parallel. In the case of a Voronoi diagram or Laguerre diagram, the dominance regions $\Pi^+(i, j)$ are half-spaces. Their intersection are efficiently computed by representing the cells in dual form [Aur91, L'22]. The cells computed at step (2) are shown in Figure 2-B. As one can see, the set of cells $V_i^\mathcal{E}$ deduced from the k-NN graph $\mathcal{E}$ is very near the Voronoi diagram, except for a few overlaps between the cells (in red);

* in step (3), we use an axis-aligned bounding box tree to determine a set of candidate intersection pairs (Figure 2-C). Then all the candidate intersections are processed in parallel. Each time the bounding boxes of a pair of cells $V_i^\mathcal{E}$ and $V_j^\mathcal{E}$ have an intersection, the edges $(i-j)$ and $(j-i)$ are appended to $\mathcal{E}$ (in red in Figure 2-D). This will insert too many edges, but any supergraph will give the correct result;

* in step (4), each re-computed cell $V_i^\mathcal{E}$ corresponds to the Voronoi cell $V_{or_i}$ (since $Del(\mathbf{X}) \subseteq \mathcal{E}$). The so-computed Voronoi diagram is shown in Figure 2-E. The Delaunay graph is shown in Figure 2-F. All its edges are contained either in the k-NN graph (Figure 2-A) or in the edges coming from cell intersections (Figure 2-D).

The Delaunay triangulation, if needed, is deduced from the combinatorics of the Voronoi diagram. In 3D, a tetrahedron $(i, j, k, l)$ is generated each time the intersection between the dominance regions $\Pi^+(i, j)$, $\Pi^+(i, k)$, $\Pi^+(i, l)$ define a vertex of $V_{or_i}$. To avoid obtaining the tetrahedra four times, one can keep only the one for which $i$ is the smallest index. To make sure the same tetrahedra are generated from different Voronoi cell, one needs to take care of configurations with cospherical points, that generate Voronoi vertices of degree potentially larger than 4 (that is, incident to more than 4 Voronoi regions), hence several possible tetrahedralizations. A simple solution is to use the so-called classical symbolic perturbation approach [EM90]. Computing a Voronoi diagram requires one to deduce some combinatorics from the relative positions of geometric objects, classifying them as 'above' or 'below' for instance. This is often expressed as the sign of polynomial function of the point's coordinates called a predicate. The configurations where the polynomial is zero - for instance co-spherical points in Voronoi diagrams - require a special treatment. One possibility is to symbolically apply a permutation to each point's coordinate, and use the sign of the first non-zero term of the Taylor expansion. The only predicate used by the computation of the cells $V_i^\mathcal{E}$ requires to classify a Voronoi vertex, that is, the intersection of three bisectors, relative to a fourth bisector, which is equivalent to the in_sphere predicate (for Voronoi diagrams), or its weighted generalization (for Laguerre diagrams). Symbolically-perturbed implementations of these predicates are readily available [She97, Lév16, Lmc23, The23].

This algorithm is very similar to star splaying [She05], that also starts from a subset of the Delaunay graph, and finds the missing edges. In other words, one may describe this algorithm as “star splaying seen from the point of view of the Voronoi cells”. It adds an alternative way of considering the neighborhoods,
in addition to the three points of view usually associated with star splaying ("star of Delaunay vertex", "convex hull of ray" and "cross-section of cone"). The four points of view are formally equivalent, but focusing on Voronoi cells has three practical consequences:

* it may be a question of taste, but one may find the explanation in terms of Voronoi cells shown in Figure 2 easier to follow: simply put, the dual of "star splaying" is "Voronoi clipping" or "Voronoi trimming";

* instead of storing all the vertex's neighborhoods, one can store only the graph, and recompute the Voronoi cells as needed. While this has an additional computational cost, it saves substantial amounts of memory;

* by considering the Voronoi cells explicitly, one can exploit their relations with the boundary of the domain to optimize point exchanges in the distributed algorithm, as explained in the next section.

## 2.3 DVD: Distributed Voronoi Diagram

We now suppose that the domain is partitioned into $M$ regions $(R_k)_{k=1}^M$, and that the pointset $\mathbf{X}$ is partitioned accordingly, into $M$ subsets $\mathbf{X}_k = \mathbf{X} \cap R_k$. In addition, we suppose that each region $R_k$ is represented and processed by a different node of a cluster. Each node $k$ is going to exchange some points with other nodes in order to ensure that it obtains all the necessary information to compute all the Voronoi cells of the points $\mathbf{X}_k$.

**Algorithm 4.** By-region parallel Voronoi Diagram

**Input:** the regions $\{R_k\}_{k=1}^M$ and the pointsets $\{\mathbf{X}_k\}_{k=1}^M$

**Output:** $M$ graphs $\mathcal{E}_k$, such that $\forall i = V_{or_i} = V_{\mathcal{E}_k} \ \forall i$ such that $\mathbf{x}_i \in R_k$

(1) **for** $k = 1...M$, $\mathcal{E}_k \leftarrow \text{Del}(\mathbf{X}_k)$

(2) **for** $k = 1...M$, $\mathbf{Y}_k \leftarrow \{\mathbf{x}_i \in R_k \mid l \neq k \text{ and } V_{\mathcal{E}_k} \cap R_k \neq \emptyset\}$

(3) **for** $k = 1...M$, $\mathcal{E}_k \leftarrow \text{Del}(\mathbf{X}_k \cup \mathbf{Y}_k)$

(4) **for** $k = 1...M$, $\mathbf{Z}_k \leftarrow \{\mathbf{x}_j \mid \exists l \neq k, \exists (i \to j) \in \mathcal{E}_l, \mathbf{x}_i \in \mathbf{X}_k, \mathbf{x}_j \in \mathbf{X}_l\}$

(5) **for** $k = 1...M$, $\mathcal{E}_k \leftarrow \text{Del}(\mathbf{X}_k \cup \mathbf{Y}_k \cup \mathbf{Z}_k)$

As in the previous section, in step (1), the algorithm starts by computing a graph $\mathcal{E}$, but this time the graph is distributed over a set of $M$ nodes, $\mathcal{E} = \mathcal{E}_1 \cup ... \cup \mathcal{E}_M$, and each component $\mathcal{E}_k$ stored in a node corresponds to the Delaunay graph $\text{Del}(\mathbf{X}_k)$ of the vertices of $\mathbf{X}$ that fall in the region $R_k$ associated with
<figure><img src="image_3.png" /><figcaption>Figure 3: A simple example showing why two steps of point exchange are needed in the distributed Voronoi diagram algorithm.</figcaption></figure>

that node. To compute the Delaunay graph, one may use the standard Bowyer-Watson algorithm [Bow81, Wat81], or Algorithm 3 in the previous section. Once the Delaunay graph $\mathcal{E}_k$ is computed, most cells $V_i^{\mathcal{E}_k}$ correspond to the Voronoi cells $V_{or_i}$ (since $\mathcal{E}_k$ is the Delaunay graph of $\mathbf{X}_k$), except for some cells on the boundary, that are larger than the Voronoi cells (Observation 1). The edges $(i-j)$ that are missing in $\mathcal{E}_k$ correspond to pairs of cells with non-empty intersections (Observation 2), stored in different nodes. To have in each $\mathcal{E}_k$ all the edges of the Delaunay graph, we need to determine both the relation $i$ is clipped by $j$ (that corresponds to the oriented edge $(i \to j)$ and the relation $i$ is clipping $j$ (that corresponds to the oriented edge $(j \to i)$), hence we have an algorithm that operates in two phases (steps (2,3) and steps (4,5)): in step (2), each node $k$ gathers the points $\mathbf{x}_i \in R_l$ with a cell $V_i^{E_k}$ that has a non-empty intersection with $R_k$ (to detect the missed is clipped by relations). In step (3), the graphs are updated. In step (4), the so-discovered new edges are symmetrized, by sending back the neighbors (is clipping relations). Each graph $\mathcal{E}_k$ obtained in the node $k$ at step (5) has all the edges to properly compute the Voronoi cells of the points in $R_k$.

**Why do we need the two steps (2) and (4) ?** Our goal is to represent, in a distributed manner, a graph $\mathcal{E} = \bigcup_k \mathcal{E}_k$ that has all the edges of the Delaunay
graph $\mathcal{Del}(\mathbf{X})$. For each (unoriented) edge $(i \to j) \in \mathcal{Del}(\mathbf{X})$, we need to find
the two oriented edges $(i \to j)$ (that corresponds to the relation *i is clipped by j*) and $(j \to i)$ (that corresponds to the relation *i is clipping j*) in all the nodes
that contain $\mathbf{x}_i$ or $\mathbf{x}_j$, hence two steps are needed. Let us examine them more
closely:

* Step (2) ensures that all for all $\mathbf{x}_i$ in a region $R_k$, all the edges $(j \to i)$ such that $(i-j) \in \mathcal{Del}(\mathbf{X})$ are in the $\mathcal{E}_{l,l \neq k}$ (that is, the *i is clipping j* relation, for the $x_j$'s in the other regions);

* Step (4) ensures that for all $\mathbf{x}_i$ in a region $R_k$, all the edges $(i \to j)$ such that $(i-j) \in \mathcal{Del}(\mathbf{X})$ are in $\mathcal{E}_k$ (that is, the *i is clipped by j* relation, for the $x_j$'s in the other regions).

Hence step (4) is simply a combinatorial operation that *symmetrizes* the ori-
ented edges discovered at step (2).

Note that observations 1, 2 hold in arbitrary dimension, as well as algorithm 3. To help grasping an intuition of what the algorithm does, we shall now examine the effect of the two steps on a particular 2D example, illustrated in Figure 3. One wants to compute the Voronoi diagram shown on the left, using two regions $R_1$ and $R_2$. The boundary between $R_1$ and $R_2$ is shown in red. The sets of points $\mathbf{X}_1$ and $\mathbf{X}_2$ are initialized with the points contained by regions $R_1$ and $R_2$ respectively, shown in Fig. 3-A. In the first step (Fig. 3-B),

* Region $R_1$ contains $\mathbf{x}_1$. The cell $V_1^\epsilon$ of $\mathbf{x}_1$ has an intersection with the boundary, hence $\mathbf{x}_1$ is sent to $R_2$.

* Region $R_2$ contains $\mathbf{x}_2, \mathbf{x}_3$ and $\mathbf{x}_4$. The cells of $\mathbf{x}_2$ and $\mathbf{x}_3$ have an intersection with the boundary, and are sent to $R_1$.

In the second step (Fig. 3-C), region $R_1$ knows $\mathbf{x}_1$, $\mathbf{x}_2$ and $\mathbf{x}_3$. The cell $V_1^\epsilon$ is not correct, because it does not know (yet) its neighbor $\mathbf{x}_4$, but the neighboring relation $(1-4)$ is known by $R_2$, that sends $\mathbf{x}_4$ to $R_1$.

With this by-region granularity, Algorithm 4 can be implemented in a clus-
ter. All the variables with index *k* are stored in a specific node of the cluster.
The five loops in steps (1)-(5) are executed in parallel, by all the nodes in the
cluster. The only inter-node communications are at steps (2) and (4) (plus the
initial broadcast and partition of **X** into the *R*<sub>*k*</sub>'s). In algorithm 4, at steps (2)
and (4) a node *k* gathers points from potentially all the other nodes. From an
implementation point of view, it is more natural to have the nodes intercon-
nected with only the neighboring nodes through communication channels. In
this setting, the node that corresponds to region *R*<sub>*k*</sub> is connected to the nodes
of the adjacent regions. During the execution of the algorithm, each node *sends*
the points where they are needed. This alternative point of view leads to the
following algorithm, executed by each node in parallel (the loops on the nodes
*k* are omitted):
<caption>**Algorithm 5.** *Distributed Voronoi Diagram (algorithm for a node k)*</caption>

**Input:** the regions $\{R_k\}_{k=1}^M$ and the pointsets $\{\mathbf{X}_k\}_{k=1}^M$
**Output:** $M$ graphs $\mathcal{E}_k$, such that $\text{Vor}_i = \mathcal{V}_{i}^{\mathcal{E}_k} \forall \mathbf{x}_i \in R_k$

(1) $\mathcal{E}_k \leftarrow \text{Del}(\mathbf{X}_k)$

(2) **for each** $\mathbf{x}_i \in \mathbf{X}_k$

(3) **for each** $R_l$ neighbor of $R_k$

(4) **if** $\mathcal{V}_{i}^{\mathcal{E}_k} \cap R_l \neq \emptyset$ **then** send $\mathbf{x}_i$ to neighbor node $l$

(5) $\mathbf{Y}_k \leftarrow$ receive points from neighbor nodes

(6) **sync**

(7) $\mathcal{E}_k \leftarrow \text{Del}(\mathbf{X}_k \cup \mathbf{Y}_k)$

(8) **for each** $(i \to j) \in \mathcal{E}$

(9) $l \leftarrow \text{Region}(\mathbf{x}_i)$ ; $m \leftarrow \text{Region}(\mathbf{x}_j)$

(10) **if** $l = k$ **and** $m \neq k$ **then** send $\mathbf{x}_j$ to neighbor node $m$

(12) $\mathbf{Z}_k \leftarrow$ receive points from neighbor nodes

(13) **sync**

(14) $\mathcal{E}_k \leftarrow \text{Del}(\mathbf{X}_k \cup \mathbf{Y}_k \cup \mathbf{Z}_k)$

where **sync** is a synchronization point for all the threads (one waits there until everybody has reached this point), and where $\text{Region}(l)$, $\text{Region}(m)$ denote the region that sent $l$ (resp. $m$) at step (5) (each node keeps track of who sent each point).

The algorithm detailed above, running on each cluster node $k$, sends a vertex to another node $l$ each time the Voronoi cell of the vertex has an intersection with the region $R_l$ associated with the node $l$, which is straightforward to determine if the set of regions form a regular grid. In our specific case (cosmological simulation), each region contains a large number of points filling the entire region (even if density can vary a lot locally), see Figure 4. While the situation where a Voronoi cell traverses several regions may theoretically happen in general, in our specific application, assuming that it does not happen is a viable assumption. Hence it is possible to make an optimization, simplifying the connection graph between the nodes of the cluster: Each cluster node $k$, that represents region $R_k$ is typically connected to nodes corresponding to the regions adjacent to $R_k$. If regions are bounded by polyhedra, adjacent regions can have a face, an edge or a vertex in common with $R_k$. Typically, one can use a regular grid, with up to four regions incident to each edge and up to eight regions incident to each vertex. If a candidate Voronoi cell encroaches a facet, then the corresponding point is sent to the region connected to the other side of the facet. If the candidate Voronoi cell encroaches an edge or a vertex, then the point is sent to all the regions incident to the edge or the vertex. Note that modified with this optimization, Algorithm 5 does not work if a cell traverses a region completely. If one wants to support this condition, besides not doing the optimization / using the unmodified algorithm, one can also execute the two points exchange
steps of the algorithm on the reduced cluster node graph until each point was
sent to all the regions it has a non-empty intersection with.

The algorithm also works for Laguerre diagrams, but one needs to take care:
in a Laguerre diagram, the cell of $\mathbf{x}_i$ does not necessarily contain $\mathbf{x}_i$, and can
even have no intersection with $R_k$ (it can be totally on the other side, inside
$R_l$).

This algorithm is similar to the one in [PMP14], that shares the same global
structure (and the same restriction that a cell should not traverse a region), with
the two steps to recover the missing edges in the Delaunay graph. The main
difference is in Step (4): instead of considering the cells, their method is based
on the vertices, that should satisfy the empty ball property. The vertex common
to cells $V_i^\epsilon, V_j^\epsilon, V_k^\epsilon$ is a Voronoi vertex if the circumscribed ball to the triangle
$(i, j, k)$ does not contain any point of $\mathbf{x}_i$. Hence, the points $\mathbf{x}_i, \mathbf{x}_j, \mathbf{x}_k$ are sent to
the regions that have a non-empty intersection with the ball (in 3D, there are
4 cells meeting at a vertex, and it is the circumscribed ball to a tetrahedron).
Infinite cells require a special treatment, and consider the cone defined by the
edges (triangles in 3D) on the border of the Delaunay graph. On the one hand,
computations are simpler than in our algorithm that needs computing Voronoi
cells explicitly. On the other hand, the objects considered for sending a point
$\mathbf{x}_i$ to the neighboring cells are larger than in our case. In our case, the cost
of sending more point than necessary is larger than the cost of computing the
cells. Moreover, in our cosmological application shown in the next section, we
can have a huge variation of pointset density, potentially creating pencil-shaped
Voronoi cells on the border. Hence it is important to be able to determine which
elements of the region boundary they encroach in order to avoid sending the
point to too many nodes.

# 3 Large-scale optimal transport for cosmology

## 3.1 An example in early Universe reconstruction

We shall now see how the algorithm behaves in a practical scenario, in Early Universe Reconstruction. Figure 4-A shows the result of a cosmological simulation of structure formation, in a cube of 300 Mpc/h edge length in cosmological units<fnref n="4" />. The figure shows a thin slice (5 Mpc/h thick) of the cube. Each point corresponds to a galaxy cluster. The goal of cosmological reconstruction [FMMS02, BFH+03, LMv21, vLM22, NSLM22] is to trace the trajectories of these galaxy clusters back in time, supposing that they started from a uniform initial condition. Under some simplifying assumptions (convexity of the integrated potential), it can be restated as a semi-discrete optimal transport

<footnote n="4">⁴300 Mpc/h stands for 300 mega parsecs (Mpc) in co-moving coordinates (/h). Co-moving coordinates factor-out the expansion of the Universe. From this 300 Mpc/h threshold, the large-scale structure of the Universe starts to be homogeneous. For reference, the nearest star is 1.399 parsecs away from the sun, and the radius of the observable Universe is around 14.25 giga parsecs.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 4: A: Thin-slice in a cosmological simulation of a 300 Mpc/h cube. Ultimately, the goal will be to compute optimal transport with 10 billion points. The cube will be decomposed into 5×5×5 subvolumes (one of them highlighted). For now, the algorithm is tested and analyzed with periodic boundary conditions (B) in one of the subvolumes (C) with 100 million points. The reconstructed initial condition (D).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Newton iterations of semi-discrete optimal transport with periodic boundary conditions</figcaption></figure>
problem. The goal, in the future, is to apply such optimal-transport-based re-
construction to massive pointsets, with billions of particles. When reconstruct-
ing from observational data, such gigantic-scale problems appear when modeling
the unobserved dark matter as additional particles [NSLM22]. They are also
needed in cosmological simulation of modified gravity [LBM24]. Typically, for
our cube of 300 Mpc/h, it would be interesting to be able to solve optimal-
transport problems with 2560³ points, which would give sufficient resolution to
track fine-scaled dynamics. The idea will be to run the Distributed Voronoi
Diagram on a PC cluster with 125 machines, each machine will be responsible
for a cube of 60 Mpc/h with 100 million particles (512³), depicted as a small
square in the corner of Figure 4-A. Such a large-scale experiment will be the
subject of later work. For now, we conducted a "proof-of-concept" experiment:
as often done in cosmology and material science, to simulate a large volume of a
globally homogeneous medium, one can use periodic boundary conditions. For
instance, one could imagine considering the small cube in Figure 4 and connect
its opposite faces before reconstructing the motion of the galaxy clusters. One
then obtains an initial condition that paves the space, as the 2D example shown
in Figure 4-B.

Figure 4-C shows the final timestep of a cosmological simulation, with 100 million points, in a cube of 60 Mpc/h edge length (same size of the small cube in Figure 4-A) with a periodic boundary condition. We made this simulation with a slightly modified version of [CTP94] (modifications are to make it work with 512³ points, original version was limited to 256³). Starting from this pointset, one can reconstruct the initial condition using optimal transport. In the end, one obtains a reconstructed initial condition, the union of the Laguerre cells shown as the yellow shape. This shape paves the space, its curvy boundaries correspond to matter that crossed the periodic boundaries. It represented by a Laguerre diagram (Figure 4-D), and each Laguerre cells correspond to the matter that lumped to form one of the galaxy clusters. The closeup in Figure 4-D shows what these Laguerre cells look like in the central region of the cube.

Computing Voronoi and Laguerre diagrams with periodic boundary condi-
tions is just a particular case of the algorithm presented in this article. In
a 2D version, one considers 9 regions, indexed by “translation vectors” with
coefficients in {−1, 0, 1}:

<table>
  <tr>
    <td>R<sub>-1,1</sub></td>
    <td>R<sub>0,1</sub></td>
    <td>R<sub>1,1</sub></td>
  </tr>
  <tr>
    <td>R<sub>-1,0</sub></td>
    <td>R<sub>0,0</sub></td>
    <td>R<sub>1,0</sub></td>
  </tr>
  <tr>
    <td>R<sub>-1,-1</sub></td>
    <td>R<sub>0,-1</sub></td>
    <td>R<sub>1,-1</sub></td>
  </tr>
</table>

In practice, only $R_{0,0}$ is stored. In steps (4) and (10) of Algorithm 5, it “sends points to itself”, by creating new copies of the points, translated according to which facet/edge/corner of the unit cube was traversed.

The example shown in Figure 4-C,D is a good “stress test” for semi-discrete
optimal transport. Besides the number of points that starts to be large, and the
additional difficulty of the periodic boundary conditions [BPR23], this optimal
transport problem is especially challenging, because at this small scale (relative
<table><caption>Table 1: Statistics of the Laguerre diagram computations during the Newton iterations of a semi-discrete optimal transport problem with 134M points and periodic boundary conditions. All times are in s.</caption><thead><tr><th rowspan="2">Iter.</th><th colspan="4">phase I</th><th colspan="3">phase II</th><th rowspan="2">total<br>t</th></tr><tr><th>npts cross</th><th>npts out</th><th>t<sub>class</sub></th><th>t<sub>insert</sub></th><th>npts</th><th>t<sub>class</sub></th><th>t<sub>insert</sub></th></tr></thead><tbody><tr><td>0</td><td>646800</td><td>0</td><td>2.2</td><td>0.5</td><td>233050</td><td>3.4</td><td>0.3</td><td>22.2</td></tr><tr><td>1</td><td>646894</td><td>1</td><td>2.1</td><td>0.5</td><td>233042</td><td>3.2</td><td>0.3</td><td>18.2</td></tr><tr><td>2</td><td>648153</td><td>455</td><td>2.0</td><td>0.5</td><td>232706</td><td>3.2</td><td>0.3</td><td>18.1</td></tr><tr><td>3</td><td>650365</td><td>2577</td><td>2.2</td><td>0.6</td><td>232284</td><td>3.2</td><td>0.4</td><td>19.2</td></tr><tr><td>4</td><td>656861</td><td>13925</td><td>2.0</td><td>0.6</td><td>232124</td><td>3.2</td><td>0.3</td><td>17.7</td></tr><tr><td>5</td><td>665621</td><td>33767</td><td>2.1</td><td>0.6</td><td>232375</td><td>3.2</td><td>0.3</td><td>18.3</td></tr><tr><td>6</td><td>679112</td><td>74305</td><td>2.0</td><td>0.6</td><td>233205</td><td>3.4</td><td>0.3</td><td>19.6</td></tr><tr><td>7</td><td>698330</td><td>154415</td><td>2.2</td><td>0.6</td><td>234937</td><td>3.1</td><td>0.3</td><td>18.8</td></tr><tr><td>8</td><td>727751</td><td>293985</td><td>1.6</td><td>0.7</td><td>237624</td><td>3.2</td><td>0.3</td><td>17.8</td></tr><tr><td>9</td><td>769940</td><td>518739</td><td>2.0</td><td>0.9</td><td>243229</td><td>3.2</td><td>0.4</td><td>22.0</td></tr><tr><td>10</td><td>773413</td><td>540357</td><td>1.9</td><td>0.9</td><td>243777</td><td>3.4</td><td>0.3</td><td>18.9</td></tr><tr><td>11</td><td>780948</td><td>583207</td><td>2.2</td><td>0.8</td><td>245062</td><td>3.1</td><td>0.4</td><td>19.0</td></tr><tr><td>12</td><td>787834</td><td>624921</td><td>2.1</td><td>0.8</td><td>246301</td><td>3.4</td><td>0.3</td><td>19.7</td></tr><tr><td>13</td><td>838950</td><td>984704</td><td>2.1</td><td>1.0</td><td>257811</td><td>3.4</td><td>0.3</td><td>22.2</td></tr><tr><td>14</td><td>874334</td><td>1306647</td><td>2.3</td><td>1.2</td><td>269045</td><td>3.5</td><td>0.4</td><td>21.9</td></tr><tr><td>15</td><td>921254</td><td>1898587</td><td>2.0</td><td>1.3</td><td>289628</td><td>3.5</td><td>0.4</td><td>18.7</td></tr><tr><td>16</td><td>952246</td><td>2407000</td><td>1.9</td><td>1.5</td><td>306489</td><td>3.9</td><td>0.5</td><td>21.5</td></tr><tr><td>17</td><td>978500</td><td>2873474</td><td>2.0</td><td>1.6</td><td>320861</td><td>4.1</td><td>0.4</td><td>19.8</td></tr><tr><td>18</td><td>1024563</td><td>3768530</td><td>2.2</td><td>1.5</td><td>344896</td><td>4.5</td><td>0.4</td><td>20.2</td></tr><tr><td>19</td><td>1057950</td><td>4580610</td><td>2.0</td><td>1.8</td><td>363521</td><td>4.7</td><td>0.4</td><td>21.1</td></tr><tr><td>20</td><td>1129350</td><td>6117650</td><td>2.3</td><td>2.6</td><td>394322</td><td>4.7</td><td>0.5</td><td>22.8</td></tr><tr><td>21</td><td>1188534</td><td>7459069</td><td>2.1</td><td>2.7</td><td>417273</td><td>5.2</td><td>0.5</td><td>23.9</td></tr><tr><td>22</td><td>1234116</td><td>8635647</td><td>2.5</td><td>3.1</td><td>435938</td><td>5.4</td><td>0.5</td><td>24.4</td></tr><tr><td>23</td><td>1266707</td><td>9666107</td><td>2.3</td><td>3.5</td><td>451221</td><td>5.5</td><td>0.6</td><td>24.5</td></tr><tr><td>24</td><td>1322035</td><td>11457655</td><td>2.3</td><td>3.8</td><td>478621</td><td>5.8</td><td>0.6</td><td>29.5</td></tr><tr><td>25</td><td>1359821</td><td>12802087</td><td>2.2</td><td>4.3</td><td>499038</td><td>5.8</td><td>0.6</td><td>25.5</td></tr><tr><td>26</td><td>1411051</td><td>14836994</td><td>2.3</td><td>5.8</td><td>532127</td><td>6.3</td><td>0.3</td><td>27.3</td></tr><tr><td>27</td><td>1411345</td><td>14853396</td><td>2.5</td><td>4.9</td><td>532418</td><td>6.3</td><td>0.3</td><td>27.1</td></tr><tr><td>28</td><td>1413260</td><td>14984207</td><td>2.7</td><td>4.9</td><td>534741</td><td>6.7</td><td>0.3</td><td>31.1</td></tr><tr><td>29</td><td>1421313</td><td>15472725</td><td>2.6</td><td>5.3</td><td>542958</td><td>6.6</td><td>0.3</td><td>32.1</td></tr><tr><td>30</td><td>1432308</td><td>16207255</td><td>2.4</td><td>6.3</td><td>557217</td><td>6.6</td><td>0.4</td><td>28.8</td></tr><tr><td>31</td><td>1443959</td><td>16951863</td><td>2.7</td><td>5.7</td><td>574707</td><td>6.8</td><td>0.3</td><td>33.2</td></tr><tr><td>32</td><td>1444741</td><td>16962135</td><td>2.5</td><td>5.8</td><td>575423</td><td>6.9</td><td>0.4</td><td>32.9</td></tr><tr><td>33</td><td>1444765</td><td>16962171</td><td>2.4</td><td>5.3</td><td>575465</td><td>7.0</td><td>0.4</td><td>32.1</td></tr><tr><td>34</td><td>1444765</td><td>16962171</td><td>2.4</td><td>5.4</td><td>575466</td><td>6.8</td><td>0.4</td><td>32.6</td></tr></tbody></table>
<table><caption>Table 2: Optimal transport for 134M points with periodic boundary conditions. Timing breakdown</caption><tr><td>total time</td><td>100%</td><td>1:12:37</td></tr><tr><td>Laguerre</td><td>29.8%</td><td>0:20:59</td></tr><tr><td>Linear solve</td><td>37.89%</td><td>0:27:31</td></tr><tr><td>Eval gradient</td><td>1.3%</td><td>58s</td></tr><tr><td>Assemble Hessian</td><td>28.1%</td><td>0:20:29</td></tr><tr><td>Misc</td><td>3.65%</td><td>0:02:39</td></tr></table>

<figure><img src="image_6.png" /><figcaption>Figure 6: Top panel: Average cell difference and worst cell difference in function of Newton iterations (in percent of target cell measure, logarithmic scale). Bottom panel: number of substeps in each Newton iteration.</figcaption></figure>
to the size of the visible Universe), the variations of density are considerable: distances between the points vary by 5 orders of magnitude. Some regions are very dense, and some regions are nearly empty. This has several consequences for the iterative algorithm that computes the Laguerre diagram. As shown in the 2D example in Figure 5, it starts with a Voronoi diagram, then evolves the weights of a Laguerre diagram to find the (unique) Laguerre diagram such that the areas of the Laguerre cells match the masses of the points. As can be seen, the large variations of points density makes the first iterations quite “explosive”. Then, as one can see, the Laguerre cells “travel” in space, to reach the position where the matter that corresponds to each cluster started from. As one can see, a Laguerre cell does not necessarily contain its generating point (unlike Voronoi cells). Another consequence of the high variations of density is that the matter that lumped to form the structures traveled a lot (up to 15% of cube edge length). In figure 4-D, they correspond to the parts of the yellow shape that crossed the boundary of the cube. In terms of our algorithm, it means that a significant number of points will be sent through the boundaries of the domain. The same behavior is expected when we will be running the algorithm on a cluster grid, with different pointsets in each cube.

In Table 1, we report different timings and statistics on Algorithm 5 executed in the context of semi-discrete optimal transport for early Universe reconstruction, with the pointset shown in Figure 4-B. The algorithm converged after 34 Newton iterations, which took 1h 15min in total, on a dual-CPU AMD EPYC 9754 (256 cores per CPU). For each iteration, the table reports the number of points with a Laguerre cell that crosses the boundary of the domain (npts cross), the number of points completely on the other side of the boundary (npts out), the time to classify the points ($t_{class}$), to insert the points into the triangulation ($t_{insert}$), for both phases of the algorithm (steps (1) to (6) and (7) to (14) in Algorithm 5 p. 19), as well as the total time for computing the Laguerre diagram (this includes the time to compute the initial triangulation in $R_{0,0,0}$). As can be seen, each iteration takes longer and longer, which is because the Laguerre cells travel a longer and longer distance (like in Figure 5), and a larger and larger number of points needs to be sent. Phase II takes a shorter time as compared to Phase I, because it is a purely combinatorial operation, whereas Phase I needs to compute the Laguerre cells and their relations with the boundary. However, execution time remains reasonable, even in the last iterations, that need to send more than 10% of the points. Figure 6 shows the convergence and the number of KMT substeps used by each Newton iteration. The average difference between a cell’s measure $|\text{Lag}_i^\psi|$ and the target measure $\nu_i$ (expressed in % of target cell measure), as well as the maximum difference in function of the iteration are both displayed. As can be seen, the first Newton iterations are challenged by the huge difference of density, and use a tiny $\alpha = 2^{-18}$ descent parameter. Starting from the 30$^th$ iteration, the algorithm enters a regime with quadratic speed of convergence, and the measure error of the worst cell quickly decreases. Table 2 shows the timing breakdown of the optimal transport computation. The most time consuming phases are the linear solves, the computation of the Laguerre diagrams and matrix assemblies.
# 4 Conclusions and future works

The numerical experiments presented here are a “proof of concept”: this “mock-up” numerical experiment demonstrates the feasibility of the approach for the large-scale experiment illustrated in Figure 4-A: the distributed Voronoi diagram algorithm exchanges a larger and larger number of points as Laguerre cells travel through the boundaries of the periodic domain. In the large-scale scenario, a similar number of points will be exchanged (but this time, with neighboring cluster nodes instead of the same cube seen from periodic boundaries). The experiment also demonstrates that the algorithm can handle the huge variations of density, that challenge both its geometric part (Laguerre diagram) and numeric part (linear systems).

The software prototype is useful as it is for optimal-transport reconstruction with periodic boundary conditions, and was successfully applied it to problems with sizes ranging from 16M to 300M points [LMv21, vLM22, NSLM22, NPL+23]. Now the goal is to scale-up to problems of 1 billion points and above. This will make it possible to work with a volume of 300 Mpc/h (Figure 4-A) with the same resolution, and better measure different phenomena at a fine scale. This means running a cube of $5 \times 5 \times 5$ interconnected instances of the 60 Mpc/h reconstruction shown here on a cluster, which will be done in future works. Besides the implementation of the DVD algorithm on a cluster, this will also require a distributed linear solver (see [PSY+15] and references herein). Efforts are also underway to continue the development of a version where cells are computed in a fully parallel fashion, bypassing the Bowyer-Watson algorithm, and leading naturally to other ways of managing separate diagrams on multiple machines.

**Acknowledgements** The authors wish to thank Roya Mohayaee, Sebastian von Hausegger and Farnik Nikakhtar for testing, using and reporting bugs with the early version of the algorithm. We wish also to thank Rémi Flamary for discussions, as well as the anonymous reviewers for the suggestions they made that helped improving this article.

B. Lévy’s work is supported by the COSMOGRAM-launchpad Inria exploratory action (AeX). Q. Mérigot and H. Leclerc acknowledge the support of the Agence nationale de la recherche, through the PEPR PDE-AI project, contract number ANR-23-PEIA-0004. Experiments presented in this article were carried out using the Grid’5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr)

# References

[AHA98] Franz Aurenhammer, Friedrich Hoffmann, and Boris Aronov. Minkowski-type theorems and least-squares clustering. *Algorithmica*, 20:61–76, 1998.
[Aur91] Franz Aurenhammer. Voronoi diagrams—a survey of a fundamental geometric data structure. *ACM Comput. Surv.*, 23(3):345–405, sep 1991.

[BAR+21] Justine Basselin, Laurent Alonso, Nicolas Ray, Dmitry Sokolov, Sylvain Lefebvre, and Bruno Lévy. Restricted power diagrams on the GPU. *Comput. Graph. Forum*, 40(2):1–12, 2021.

[BB00] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the monge-kantorovich mass transfer problem. *Numerische Mathematik*, 84(3):375–393, 2000.

[BC13] Lorenzo Brasco and Guillaume Carlier. Congested traffic equilibria and degenerate anisotropic pdes. *Dyn. Games Appl.*, 3(4):508–522, 2013.

[BCM24] J.-D. Benamou, C.J. Cotter, and H. Malamut. Entropic optimal transport solutions of the semigeostrophic equations. *Journal of Computational Physics*, 500:112745, 2024.

[BCMO16] Jean-David Benamou, Guillaume Carlier, Quentin Mérigot, and Édouard Oudet. Discretization of functionals involving the monge-ampère operator. *Numerische Mathematik*, 134(3):611–636, 2016.

[BFH+03] Y. Brenier, U. Frisch, M. Hénon, G. Loeper, S. Matarrese, R. Mohayaee, and A. Sobolevskii. Reconstruction of the early universe as a convex optimization problem. *Mon. Not. R. Astron. Soc.*, 346(501), 2003.

[BFO14] Jean-David Benamou, Brittany D Froese, and Adam M Oberman. Numerical solution of the optimal transportation problem using the monge-ampère equation. *Journal of Computational Physics*, 260:107–126, 2014.

[Bow81] Adrian Bowyer. Computing dirichlet tessellations. *Comput. J.*, 24(2):162–166, 1981.

[BPR23] D.P. Bourne, M. Pearce, and S.M. Roper. Geometric modelling of polycrystalline materials: Laguerre tessellations and periodic semi-discrete optimal transport. *Mechanics Research Communications*, 127:104023, 2023.

[Bre91] Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. *Communications on Pure and Applied Mathematics*, 44:375–417, 1991.

[Bre15] Yann Brenier. A double large deviation principle for monge-ampère gravitation. working paper or preprint, March 2015.
[BY98] Jean-Daniel Boissonnat and Mariette Yvinec. Algorithmic geometry. Cambridge University Press, 1998.

[CMYB19] Laurent Caraffa, Pooran Memari, Murat Yirci, and Mathieu Brédif. Tile & merge: Distributed delaunay triangulations for cloud computing. In Chaitanya K. Baru, Jun Huan, Latifur Khan, Xiaohua Hu, Ronay Ak, Yuanyuan Tian, Roger S. Barga, Carlo Zaniolo, Kisung Lee, and Yanfang (Fanny) Ye, editors, 2019 IEEE International Conference on Big Data (IEEE BigData), Los Angeles, CA, USA, December 9-12, 2019, pages 1613–1618. IEEE, 2019.

[CP84] M. Cullen and R. Purser. An extended lagrangian theory of semi-geostrophic frontogenesis. J. of the Atmospheric Sciences, pages 1477–1497, 1984.

[CTP94] Hugh M. P. Couchman, Peter A Thomas, and Frazer R. Pearce. Hydra: An adaptive-mesh implementation of pppm-sph. arXiv: Astrophysics, 1994.

[Cut13] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In C.J. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K.Q. Weinberger, editors, Advances in Neural Information Processing Systems, volume 26. Curran Associates, Inc., 2013.

[Dem19] D. Demidov. Amgcl: An efficient, flexible, and extensible algebraic multigrid implementation. Lobachevskii J. of Math., 40(5):535–546, May 2019.

[dGWH+15] Fernando de Goes, Corentin Wallez, Jin Huang, Dmitry Pavlov, and Mathieu Desbrun. Power particles: an incompressible fluid solver based on power diagrams. ACM Trans. Graph., 34(4):50:1–50:11, 2015.

[EBC+22] C.P. Egan, D.P. Bourne, C.J. Cotter, M.J.P. Cullen, B. Pelloni, S.M. Roper, and M. Wilkinson. A new implementation of the geometric method for solving the easy slice equations. Journal of Computational Physics, 469:111542, 2022.

[EGH00] Robert Eymard, Thierry Gallouët, and Raphaèle Herbin. Finite volume methods. Handbook of numerical analysis, 7:713–1018, 2000.

[EM90] Herbert Edelsbrunner and Ernst Peter Mücke. Simulation of simplicity: A technique to cope with degenerate cases in geometric algorithms. ACM Trans. Graph., 9(1):66–104, jan 1990.

[FMMS02] U. Frisch, S. Matarrese, R. Mohayaee, and A. Sobolevskii. A reconstruction of the initial conditions of the universe by optimal mass transportation. Nature, 417(260), 2002.
[FSV+19] Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shunichi Amari, Alain Trouve, and Gabriel Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In *The 22nd International Conference on Artificial Intelligence and Statistics*, pages 2681–2690, 2019.

[GM18] Thomas O. Gallouët and Quentin Mérigot. A lagrangian scheme à la brenier for the incompressible euler equations. *Found. Comput. Math.*, 18(4):835–865, 2018.

[Hay20] Ken Hayami. Convergence of the conjugate gradient method on singular systems, 2020.

[ILSS06] Martin Isenburg, Yuanxin Liu, Jonathan Richard Shewchuk, and Jack Snoeyink. Streaming computation of delaunay triangulations. *ACM Trans. Graph.*, 25(3):1049–1056, 2006.

[JKO98] Richard Jordan, David Kinderlehrer, and Felix Otto. The variational formulation of the fokker-planck equation. *SIAM Journal on Mathematical Analysis*, 29(1):1–17, 1998.

[KMT19] Jun Kitagawa, Quentin Mérigot, and Boris Thibert. Convergence of a newton algorithm for semi-discrete optimal transport. *J. Eur. Math. Soc.*, 21(9), 2019.

[L'22] Bruno Lévy. Partial optimal transport for a constant-volume lagrangian mesh with free boundaries. *J. Comput. Phys.*, 451(C), feb 2022.

[LB12] Bruno Lévy and Nicolas Bonneel. Variational anisotropic surface meshing with voronoi parallel linear enumeration. In Xiangmin Jiao and Jean-Christophe Weill, editors, *Proceedings of the 21st International Meshing Roundtable, IMR 2012, October 7-10, 2012, San Jose, CA, USA*, pages 349–366. Springer, 2012.

[LBM24] Bruno Lévy, Yann Brenier, and Roya Mohayaee. Monge ampère gravity: from the large deviation principle to cosmological simulations through optimal transport, 2024.

[Lév15] Bruno Lévy. A numerical algorithm for $L_2$ semi-discrete optimal transport in 3d. *ESAIM M2AN (Mathematical Modeling and Analysis)*, 2015.

[Lév16] Bruno Lévy. Robustness and efficiency of geometric programs: The predicate construction kit (PCK). *Comput. Aided Des.*, 72:3–12, 2016.

[Lév22] Bruno Lévy. Partial optimal transport for a constant-volume lagrangian mesh with free boundaries. *J. Comput. Phys.*, 451:110838, 2022.
[Lmc23] Inria (Bruno Lévy and multiple contributors). Geogram: a programming library of geometric algorithms. http://github.com/BrunoLevy/geogram, 2014-2023.

[LMSS20] Hugo Leclerc, Quentin Mérigot, Filippo Santambrogio, and Federico Stra. Lagrangian discretization of crowd motion and linear diffusion. *SIAM J. Numer. Anal.*, 58(4):2093–2118, 2020.

[LMv21] Bruno Levy, Roya Mohayaee, and Sebastian von Hausegger. A fast semidiscrete optimal transport algorithm for a unique reconstruction of the early Universe. *Monthly Notices of the Royal Astronomical Society*, 506(1):1165–1185, 06 2021.

[LS18] Bruno Lévy and Erica L. Schwindt. Notions of optimal transport theory and how to implement them on a computer. *Comput. Graph.*, 72:135–148, 2018.

[Mér11] Quentin Mérigot. A multiscale approach to optimal transport. *Comput. Graph. Forum*, 30(5):1583–1592, 2011.

[MFF+17] Rahul Mourya, André Ferrari, Rémi Flamary, Pascal Bianchi, and Cédric Richard. Distributed deblurring of large images of wide field-of-view, 2017.

[Mon84] Gaspard Monge. Mémoire sur la théorie des déblais et des remblais. *Histoire de l'Académie Royale des Sciences (1781)*, pages 666–704, 1784.

[MS24] Ismael Medina and Bernhard Schmitzer. Flow updates for domain decomposition of entropic optimal transport. *arXiv preprint arXiv:2405.09400*, 2024.

[MT21] Quentin Merigot and Boris Thibert. Optimal transport: discretization and algorithms. In *Handbook of numerical analysis*, volume 22, pages 133–212. Elsevier, 2021.

[NPL+23] Farnik Nikakhtar, Nikhil Padmanabhan, Bruno Lévy, Ravi K. Sheth, and Roya Mohayaee. Optimal transport reconstruction of biased tracers in redshift space. *Phys. Rev. D*, 108:083534, Oct 2023.

[NSLM22] Farnik Nikakhtar, Ravi K. Sheth, Bruno Lévy, and Roya Mohayaee. Optimal Transport Reconstruction of Baryon Acoustic Oscillations. *Phys. Rev. Lett.*, 129(25):251101, December 2022.

[OBSC00] Atsuyuki Okabe, Barry Boots, Kokichi Sugihara, and Sung Nok Chiu. *Spatial Tessellations: Concepts and Applications of Voronoi Diagrams*. Series in Probability and Statistics. John Wiley and Sons, Inc., 2nd ed. edition, 2000.
[PC19] Gabriel Peyré and Marco Cuturi. Computational optimal transport. *Found. Trends Mach. Learn.*, 11(5-6):355–607, 2019.

[PMP14] Tom Peterka, Dmitriy Morozov, and Carolyn Phillips. High-performance computation of distributed-memory parallel 3d voronoi and delaunay tessellation. In *SC '14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis*, pages 997–1007, 2014.

[PSY$^{+}$15] Jongsoo Park, Mikhail Smelyanskiy, Ulrike Meier Yang, Dheevatsa Mudigere, and Pradeep Dubey. High-performance algebraic multi-grid solver optimized for multi-core based distributed parallel systems. SC '15, New York, NY, USA, 2015. Association for Computing Machinery.

[QLdGJ22] Ziyin Qu, Minchen Li, Fernando de Goes, and Chenfanfu Jiang. The power particle-in-cell method. *ACM Trans. Graph.*, 41(4):118:1–118:13, 2022.

[RSLL18] Nicolas Ray, Dmitry Sokolov, Sylvain Lefebvre, and Bruno Lévy. Meshless voronoi on the GPU. *ACM Trans. Graph.*, 37(6):265, 2018.

[Sam90] Hanan Samet. *The design and analysis of spatial data structures*. Addison-Wesley Longman Publishing Co., Inc., USA, 1990.

[San15] Filippo Santambrogio. Optimal transport for applied mathematicians. *Birkhäuser, NY*, 55(58-63):94, 2015.

[She97] Jonathan Richard Shewchuk. Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates. *Discrete & Computational Geometry*, 18(3):305–363, October 1997.

[She05] Richard Shewchuk. Star splaying: an algorithm for repairing delaunay triangulations and convex hulls. In *Proceedings of the Twenty-First Annual Symposium on Computational Geometry*, SCG '05, page 237–246, New York, NY, USA, 2005. Association for Computing Machinery.

[Slo87] S.W. Sloan. A fast algorithm for constructing delaunay triangulations in the plane. *Advances in Engineering Software (1978)*, 9(1):34–55, 1987.

[Spr10] Volker Springel. E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh. *Monthly Notices of the Royal Astronomical Society*, 401(2):791–851, 01 2010.

[The23] The CGAL Project (multiple authors). *CGAL User and Reference Manual*. CGAL Editorial Board, 5.6 edition, 2023.
[Vil03] Cédric Villani. *Topics in optimal transportation*. Graduate studies in mathematics. American Mathematical Society, Providence (R.I.), 2003.

[Vil09] Cédric Villani. *Optimal transport : old and new*. Grundlehren der mathematischen Wissenschaften. Springer, Berlin, 2009.

[vLM22] Sebastian von Hausegger, Bruno Lévy, and Roya Mohayaee. Accurate Baryon Acoustic Oscillations Reconstruction via Semidiscrete Optimal Transport. *Phys. Rev. Lett.*, 128(20):201302, May 2022.

[Wat81] David Watson. Computing the n-dimensional delaunay tessellation with application to voronoi polytopes. *Comput. J.*, 24(2):167-172, 1981.