# Unidimensional semi-discrete partial optimal transport

Adrien Cances, Hugo Leclerc

► To cite this version:

Adrien Cances, Hugo Leclerc. Unidimensional semi-discrete partial optimal transport. 2025. hal-05252228

HAL Id: hal-05252228

https://hal.science/hal-05252228v1

Preprint submitted on 12 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# UNIDIMENSIONAL SEMI-DISCRETE PARTIAL OPTIMAL TRANSPORT

ADRIEN CANCES AND HUGO LECLERC

**ABSTRACT.** We study the semi-discrete formulation of one-dimensional partial optimal transport with quadratic cost, where a probability density is partially transported to a finite sum of Dirac masses of smaller total mass. This problem arises naturally in applications such as risk management, the modeling of crowd motion, and sliced partial transport algorithms for point cloud registration. Unlike higher-dimensional settings, the dual functional in the unidimensional case exhibits reduced regularity. To overcome this difficulty, we introduce a regularization procedure based on thickening the density along an auxiliary dimension. We prove that the maximizers of the regularized dual problem converge to those of the original dual problem, with quadratic rate in the introduced thickness. We further provide a numerical scheme that leverages the regularized functional, and we validate our analysis with simulations that confirm the quadratic convergence rate. Finally, we compare the semi-discrete and fully discrete settings, demonstrating that our approach offers both improved stability and computational efficiency for unidimensional partial transport problems.

## 1. INTRODUCTION

The standard definition of optimal transport requires both measures to have the same (finite) mass. In [5], Caffarelli and McCann introduced the notion of *partial optimal transport*<fnref n="1" />, which lifts this assumption and asks how to move a prescribed fraction of mass from the first measure to the second one, in the cheapest possible way. As a result, partial optimal transport allows to compare two measures of *different masses*. There is an extensive literature on the case where both the target and source measures are absolutely continuous. For the quadratic cost, existence and uniqueness of solutions, as well as regularity of the free boundaries — the boundaries of the *active regions*, on which sits the respective amounts of mass of the source and target measures that are actually transported — have been studied in the seminal work [5], and in [11, 16]. A few years later, Chen and Indrei study the free boundary regularity for a more general cost function in [6], and, in [8], Dávila and Kim inspect the evolution of the free boundaries and the monotonicity of the optimal potential as the fraction of mass that must be transported varies. At last, in 2018, the papers of Igbida and Nguyen [14, 15] address the partial optimal transport problem for a general cost and with generic source and target measures. The authors introduce a PDE of Monge-Kantorovich type which enables them to extend the uniqueness

<footnote>Date: September 11, 2025.</footnote>

<footnote n="1">¹Several authors (including Caffarelli and McCann) use the term *optimal partial transport*, but both expressions are found in the literature.</footnote>
and monotonicity results to this more general framework, and is also convenient for
numerical analysis and computations. It is interesting to note that partial optimal
transport is a special case of *unbalanced* optimal transport, obtained by using total
variation for the penalization terms. We refer to [7, Section 2.3] for a detailed
description of this result.

1.1. **Sliced partial optimal transport.** As its balanced version, partial optimal transport is a computationally expensive problem, especially in high dimension. This has led several authors [20, 3] to work on *sliced* partial optimal transport: the two measures are projected on all lines passing through the origin, and the corresponding partial transport costs are integrated. In the spirit of the sliced-Wasserstein distance, which was defined by Peyré et al. [20] and Rabin et al. [3], this definition leverages the fact that the problem is easier to solve in dimension one. However, while one-dimensional transport is trivial in the sense that for empirical measures, it amounts to a simple sort, one-dimensional partial transport is much more complex.

1.2. **Motivation for one-dimensional partial optimal transport.** In 2019, Bonneel et al. [2] have proposed an algorithm to solve the fully-discrete case. More precisely, their algorithm tackles the case where the two measures are (unweighted) sums of Dirac masses of $\mathbb{R}$, of the form $\sum_{x \in X} \delta_x$, with sets of points $X, \tilde{X}$ of different cardinalities $N \neq \tilde{N}$. The goal is to injectively assign every point of the smaller point cloud to a point of the larger one. The method can then be used to compute an approximation of the sliced partial optimal transport cost between two (unweighted) sets of points of different cardinalities. More recently, Schmitzer et al. [1] used a dual formulation to create an algorithm for the same problem, although slightly more general, since a tunable Lagrange parameter allows for leaving certain points of the source out of the mapping. Both methods have a quadratic worst-case time complexity, but in practice they achieve quasi-linear complexity. One of the applications highlighted in [2, 1] is an improved version of the Iterative Closest Point (ICP) algorithm for point cloud registration. The latter problem consists in estimating, among a certain class of transformations, the one that best aligns a given point cloud to another one (of different cardinality). While the ICP algorithm is quite efficient when considering rigid transformations (rotation and translation), it behaves in practice much less well for similarities (rotation, translation and scaling) because the scaling factor tends to go to zero as the iterations progress. Bonneel et al. showed that a sliced-OPT-based ICP helps tackle this issue by providing an injective matching at each iteration, while avoiding the intractability of full-fledged optimal transport on large problems in dimension $d \ge 2$. Schmitzer et al., on the other hand, developed an ICP-type algorithm which harnesses sliced partial optimal transport to be more robust to outliers, whose proportion is assumed to be known.

1.3. **Semi-discrete, one-dimensional, partial optimal transport.** In this paper, we focus on the semi-discrete setting. Specifically, the problem consists in transporting part of a given probability density to a (positively) weighted sum of Dirac measures whose total mass is less than one, so as to minimize the quadratic transport cost. This problem can be applied to risk management, see [9], where the goal is to solve a partial multimarginal optimal transport problem in which each one of the
D marginals is a univariate probability density $\rho_j$. If one searches an approximate solution as a discrete measure $\mu = \sum_{i=1}^N \alpha_i \delta_{y_i}$, where $y_1, \dots, y_N \in \mathbb{R}^D$ and where $\alpha_i, \dots, \alpha_N$ are non-negative weights summing to less than one, then one-dimensional semi-discrete partial optimal transport is a natural way to penalize the constraints. Indeed, the smaller the optimal partial transport cost from $\rho_j$ to $\pi_j^\mu = \sum_{i=1}^N \alpha_i \delta_{y_i}^j$, the closer the $j$th marginal of $\mu$ is to being dominated by $\rho_j$. This type of penalization already appears in [18] for solving a crowd motion problem via Lagrangian discretization, although in dimension two. The crowd is modeled as a probability density, constrained to be bounded by one to account for the incompressibility of individuals when densely packed, see for instance [19]. When passing to a discrete measure for numerical issues, this congestion constraint is dealt with using a penalty term proportional to the partial optimal transport cost between the discrete probability measure and the (unnormalized) Lebesgue measure on the domain. This ensures that the discrete measure in question is close to a density bounded by one.

1.4. **Precise setting.** Let us now introduce a few notations to better define the problem addressed in this work. We denote $\rho$ the probability density on $\mathbb{R}$ and $\sum_{i=1}^N \alpha_i \delta_{y_i}$ the discrete measure, where the $y_i$ are pairwise distinct points of $\mathbb{R}$ and where the $\alpha_i > 0$ sum to $\|\alpha\|_1 \in (0,1)$. The goal is to minimize the (quadratic) transport cost $\int |x-y|^2 d\gamma(x, y)$ over all measures $\gamma$ on $\mathbb{R} \times \mathbb{R}$ whose first marginal is bounded above by $\rho$ on all Borel sets and whose second marginal is $\mu$. By standard duality arguments, the problem is equal to its dual

$$ \max_{\psi \in \mathbb{R}^N} \left\{ \int_{\mathbb{R}} \min(0, \min_{1 \le i \le N} [|x - y_i|^2 - \psi_i]) d\rho(x) + \sum_{i=1}^N \alpha_i \psi_i \right\}, \quad (1.1) $$

which is of particular interest since the optimized dual variable $\psi$ lives in finite dimension. This kind of formulation was already used in [18] by Leclerc and co-authors, for semi-discrete optimal transport in dimension two or three. However, in dimension one, the problem is curiously much harder to solve numerically. We interpret this numerical difficulty as a consequence of the lack of regularity of the dual functional in dimension one, compared to dimension two or more (replacing absolute values with Euclidean norms), see Theorem 8 and Example 12. Indeed, whereas in the latter case the functional is of class $C^2$ on a large open containing the maximizers, in the former case the functional is only $C^1$, and its second derivative has singularities. In order to circumvent this difficulty, we minimize a regularized functional, which in fact corresponds to solving a two-dimensional transport problem, where the one-dimensional density $\rho$ is thickened by a width $2\epsilon$ along a second dimension. Thanks to the $C^2$ regularity of the new functional, the Hessian matrix is well-defined. The new problem reads

$$ \max_{\psi \in \mathbb{R}^N} -\epsilon^2 \int_{\mathbb{R}} f^*(\epsilon^{-2} \max_{1 \le i \le N} (\psi_i - |x - y_i|^2)) d\rho(x) + \sum_{i=1}^N \alpha_i \psi_i \quad (1.2) $$
where

$$f^*(t) = \begin{cases} 0 & \text{if } t < 0, \\ \frac{2}{3}t^{\frac{3}{2}} & \text{if } 0 \le t \le 1, \\ t - \frac{1}{3} & \text{if } t > 1. \end{cases} \qquad (1.3)$$

Our main theoretical result is the following theorem, which guarantees that the regularized problem converges at a fast rate to the initial one as the parameter $\epsilon$ goes to zero. Moreover, Example 18 shows that our result is sharp.

**Theorem 1.** Let $\psi^\epsilon \in \mathbb{R}^N$ be the maximizer of Problem (1.2). Under some regularity assumptions on the density, we have

$$\|\psi^\epsilon - \psi^*\| \lesssim \epsilon^2 \qquad (1.4)$$

where $\psi^* \in \mathbb{R}^N$ is the maximizer of Problem (1.1) and where $\lesssim$ hides a constant depending only on $\rho$, on the $y_i$ and on the $\alpha_i$.

The proof harnesses standard tools of semi-discrete optimal transport, and spectral graph theory in order to show that, in dimension $d \ge 2$ and assuming the probability density has convex, compact support, the dual functional of semi-discrete partial optimal transport is strictly concave on some subset containing the solution. In particular, this implies that the respective solutions $\psi^\epsilon$ of the regularized problems are unique. We emphasize that an important feature of partial optimal transport is its equivalence to a balanced version of the problem (where the two measures have equal mass), up to introducing a fictitious point and extending the cost function to zero for this particular point. This corresponds to imagining that the inactive part of the source measure — the part that is not sent to the target measure — is actually getting transported “for free” to the auxiliary point. This point of view was introduced by Caffarelli and McCann in [5, Section 2], in the case of absolutely continuous source and target measures, but the equivalence still holds in the more general framework of generic probability measures. As a consequence for semi-discrete partial transport, strict concavity of the dual functional amounts to connectedness of the (undirected) graph with vertices $\{1, \dots, N, \infty\}$ and weights given by the Hessian of the functional (the diagonal coefficients weigh the edges between $\infty$ and the other vertices). The reason why the dual variable $\psi$ is unique in semi-discrete partial optimal transport (assuming the density has connected support) is that the “hidden” component of $\psi$ associated to the fictitious point is implicitly set to zero. This eliminates the underdeterminacy of balanced semi-discrete optimal transport, for which adding the same constant to each component of the dual variable yields the same value.

To show quadratic convergence of $\psi^\epsilon$ towards $\psi^*$, we first differentiate in $\epsilon$ the optimality conditions satisfied by $\psi^\epsilon$ to get an ODE of the form $D^2K^\epsilon(\psi^\epsilon)\dot{\psi}^\epsilon + \partial_\epsilon[\nabla K^\epsilon](\psi^\epsilon)$. The Hessian at the maximizer can be bounded below by a positive-definite matrix independent of $\epsilon$, once again thanks to spectral graph theory, and more specifically to Laplacian matrices and a careful analysis of the second derivatives of $K^\epsilon$. As for the mixed derivative $\partial_\epsilon[\nabla K^\epsilon]$, it can be bounded above quite straightforwardly by $\epsilon$ up to a multiplicative constant, and the quadratic rate follows by integrating the resulting bound on $\dot{\psi}^\epsilon$.
**1.5. Outline.** After a brief overview of optimal transport, we introduce its semi-discrete partial counterpart, in its primal and dual formulations. We establish that the dual functional is of class $C^2$ for dimensions greater than one, whereas its Hessian has singularities in the case of dimension one. Next, we present the regularized problem and prove Theorem 1 by introducing an ODE that is satisfied by the family of solutions $\epsilon \mapsto \psi^\epsilon$. The final section is dedicated to numerics. After a brief description of the algorithm, we provide simulation results that illustrate the quadratic convergence rate predicted by Theorem 1, and apply our algorithm to a semi-discrete version of the sliced partial transport matching. At last, we compare the discrete-discrete versus semi-discrete settings of unidimensional optimal partial transport, both in terms of accuracy and execution time, via a benchmark of Boonnel's algorithm and ours.

## 2. SEMI-DISCRETE OPTIMAL TRANSPORT AND THE CASE OF PARTIAL TRANSPORT

We first introduce the Wasserstein distance (of order 2) for probability measures on $\mathbb{R}^d$, as well as its dual formulation and the specific framework of the semi-discrete case.

**Definition 2.** For two probability measures with finite second moment $\rho, \mu \in \mathcal{P}_2(\mathbb{R}^d)$, the 2-Wasserstein distance is defined as

$$W_2(\rho, \mu) = \left( \inf_{\gamma \in \Gamma(\rho, \mu)} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x - y\|^2 d\gamma(x, y) \right)^{\frac{1}{2}}, \quad (2.1)$$

where $\Gamma(\rho, \mu)$ denotes the set of couplings (or transport plans) of $\rho$ and $\mu$, that is, probability measures on $\mathbb{R}^d \times \mathbb{R}^d$ with first marginal $\rho$ and second marginal $\mu$.

The celebrated Brenier theorem [4] states that under mild assumptions — for example if the source measure $\rho$ is absolutely continuous — the optimal transportation does not split mass. More precisely, any optimal plan is induced by an (optimal transport) map $T \in L^2(\rho)$, which is uniquely defined $\rho$-almost everywhere and (naturally) sends $\rho$ to $\mu$, in the sense that $T_{\#}\rho = \mu$. The notation $T_{\#}\rho$ stands for the pushforward of $\rho$ by $T$, defined by $[T_{\#}\rho](B) = \rho(T^{-1}(B))$ for every measurable set $B$, and we say that $T$ pushes $\rho$ forward to $T_{\#}\rho$. Any optimal transport plan thus writes $\gamma = (\mathrm{Id}, T)_{\#}\rho$, and the mass at $\rho$-almost every point $x$ is sent to a unique point $y = T(x)$.

**2.1. Kantorovich duality.** The optimal transport problem (2.1) is an infinite-dimensional linear program, and the Fenchel-Rockafellar theorem allows to establish strong duality: if $\rho$ and $\mu$ both have finite second moment, the problem is equal to its dual

$$W_2(\rho, \mu)^2 = \sup \left\{ \int_{\mathbb{R}^d} \varphi d\rho + \int_{\mathbb{R}^d} d\mu : \varphi, \psi \in \mathcal{C}_b(\mathbb{R}^d), \varphi \oplus \psi \le c \right\} \quad (2.2)$$

$$= \sup_{\psi \in \mathcal{C}_b(\mathbb{R}^d)} \left\{ \int_{\mathbb{R}^d} \psi^c d\rho + \int_{\mathbb{R}^d} \psi d\mu \right\} \quad (2.3)$$
where $\varphi \oplus \psi$ maps $(x, y) \in \mathbb{R}^d \times \mathbb{R}^d$ to $\varphi(x) + \psi(y)$, $c$ is the quadratic cost $c(x, y) = \|x - y\|^2$, and $\psi^c$ denotes the $c$-transform of $\psi$, defined by

$$ \psi^c(x) = \inf_{y \in \mathbb{R}^d} c(x, y) - \psi(y). \qquad (2.4) $$

The interested reader may refer to [22, Chapter 5] for a detailed study of the Kantorovich duality. In this work, we will often refer to the dual variable $\psi$ as the potential, whether or not it is optimal in (2.3). Note that most authors use this terminology for optimal $\psi$ only.

## 2.2. Semi-discrete transport.
In the semi-discrete case, we consider an absolutely continuous probability measure $\rho \in \mathcal{P}_2^{ac}(\mathbb{R}^d)$ with finite second moment, and a discrete probability measure $\mu = \sum_{i=1}^N \alpha_i \delta_{y_i}$, where the $y_i \in \mathbb{R}^d$ are pairwise distinct (without loss of generality) and where the $\alpha_i > 0$ sum to one. With a slight abuse of notation, throughout this work we also denote $\rho$ for the density of $\rho$. Using the fact that $\mu$ has finite support, we identify the potential with a vector $\psi \in \mathbb{R}^N$ and the dual problem reads $\sup_{\psi \in \mathbb{R}^N} \mathcal{K}(\psi)$, where the Kantorovich functional $\mathcal{K}: \mathbb{R}^d \to \mathbb{R}$ and the Laguerre cells are respectively defined by

$$ \mathcal{K}(\psi) = \sum_{i=1}^{N} \int_{\text{Lag}_i(\psi)} (\|x - y_i\|^2 - \psi_i) d\rho(x) + \sum_{i=1}^{N} \alpha_i \psi_i \qquad (2.5) $$

and

$$ \text{Lag}_i(\psi) = \{x \in \mathbb{R}^d : \|x - y_i\|^2 - \psi_i \le \|x - y_j\|^2 - \psi_j, \forall j \in [1, N]\}. \qquad (2.6) $$

The Laguerre cells form a tessellation of the space — that is, a cover with Lebesgue negligible pairwise intersection — into $N$ convex polytopes, each cell $\text{Lag}_i(\psi)$ being sent to the position $y_i$ of its corresponding Dirac mass. Indeed, the intersection of any two Laguerre cells belongs to a hyperplane orthogonal to the segment joining the two related Dirac masses, and is thus of Lebesgue measure zero. Since $\rho$ is absolutely continuous, the set of points belonging to more than one cell is $\rho$-negligible, and $\psi$ induces a transport map defined $\rho$-almost everywhere by $T_{\psi}|_{\text{Lag}_i(\psi)} \equiv y_i$. The corresponding transport plan is $\gamma = \sum_{i=1}^N (\rho \llcorner \text{Lag}_i(\psi)) \otimes \delta_{y_i}$. As a result, in the dual formulation mass cannot be split, and by strong duality the optimal plan is induced by a potential, hence by a map in $L^2(\rho)$.

**Remark 3.** The convexity of the Laguerre cells is very specific to the quadratic cost $c(x, y) = \|x - y\|^2$, since for this cost the sublevel sets of $c(\cdot, y_i) - c(\cdot, y_j) = \langle \cdot, y_j - y_i \rangle - \|y_j\|^2 + \|y_i\|^2$ are precisely half-spaces. For non-quadratic costs, the functions $c(\cdot, y_i) - c(\cdot, y_j)$ are not linear in $x$ anymore and so we cannot have convexity of both the sublevel and superlevel sets.

**Remark 4.** In semi-discrete optimal transport, by a slight abuse of notation, we also denote $\psi^c$ the $c$-transform of a vector $\psi \in \mathbb{R}^N$, which is defined by $\psi^c(x) = \min_{1 \le i \le N} \|x - y_i\|^2 - \psi_i$. Of course, it naturally corresponds to the $c$-transform of the function $\tilde{\psi}: \mathbb{R}^d \to \mathbb{R} \cup \{-\infty\}$ defined by $\tilde{\psi}(y_i) = \psi_i$ and $\tilde{\psi}(y) = -\infty$ for $y \notin \{y_1, \dots, y_N\}$.
2.3. **Partial transport.** To compare two finite measures with potentially different masses, we use a partial Wasserstein cost, which naturally extends the standard Wasserstein distance.

**Definition 5.** For two (positive) measures $\rho, \mu \in \mathcal{M}_{+,2}(\mathbb{R}^d)$ with finite second moment and such that $\rho(\mathbb{R}^d) \ge \mu(\mathbb{R}^d)$, the partial 2-Wasserstein cost is defined as

$$ \mathcal{T}_{\ge}(\rho, \mu) = \left( \inf_{\gamma \in \Gamma_{\max}(\rho, \mu)} \int_{\mathbb{R} \times \mathbb{R}} \|x - y\|^2 d\gamma(x, y) \right)^{\frac{1}{2}} \quad (2.7) $$

where $\Gamma_{\max}(\rho, \mu)$ denotes the set of measures $\gamma \in \mathcal{M}(\mathbb{R}^d \times \mathbb{R}^d)$ of mass equal to that of $\mu$, whose second marginal is $\mu$, and whose first marginal is dominated by $\rho$ in the sense that $\gamma(A \times \mathbb{R}^d) \le \rho(A)$ for all measurable $A$. The index max corresponds to the fact that $\Gamma_{\max}(\rho, \mu)$ is the set of subcouplings which have maximum mass, in the sense that as much mass as possible is transported from $\rho$ to $\mu$.

Since what used to be an equality constraint on the first marginal of $\gamma$ is now an inequality constraint, the associated Lagrange multiplier $\varphi$ must be non-positive, and the dual problem reads

$$ \mathcal{T}_{\ge}(\rho, \mu)^2 = \sup \left\{ \int_{\mathbb{R}^d} \varphi d\rho + \int_{\mathbb{R}^d} \psi d\mu : \varphi, \psi \in \mathcal{C}(\mathbb{R}^d), \varphi \oplus \psi \le c, \varphi \ge 0 \right\} \quad (2.8) $$

$$ = \sup_{\psi \in \mathcal{C}(\mathbb{R}^d)} \left\{ \int_{\mathbb{R}^d} \psi^{\tilde{c}} d\rho + \int_{\mathbb{R}^d} \psi d\mu \right\} \quad (2.9) $$

where $\psi^{\tilde{c}}$ denotes the $\tilde{c}$-transform of $\psi$, defined by $\psi^{\tilde{c}}(x) = \min\{\psi^c(x), 0\}$

2.4. **Semi-discrete partial transport.** We now make the assumption that the source density $\rho \in \mathcal{P}^{ac}(\mathbb{R}^d)$ is compactly supported and of mass normalized to one, and we take the target measure to be a discrete measure of mass strictly less than one. Namely, $\mu = \sum_{i=1}^N \alpha_i \delta_{y_i}$ where the $y_i \in \mathbb{R}^d$ are pairwise distinct and where the $\alpha_i > 0$ sum to some $\|\alpha\|_1 \in (0, 1)$. The potential is once again identified with a vector $\psi \in \mathbb{R}^N$, and for a fixed $\psi$, the optimal $\varphi$ is now the $\psi^{\tilde{c}}$-transform of $\psi$, defined by

$$ \psi^{\tilde{c}}(x) = \min \left\{ 0, \min_{1 \le i \le N} \|x - y_i\|^2 - \psi_i \right\}. $$

The Kantorovich dual thus writes $\sup_{\psi \in (0, \infty)^N} \mathcal{K}(\psi)$, where the dual functional $K: (0, +\infty) \to \mathbb{R}$ and the restricted Laguerre cells are respectively defined by

$$ \mathcal{K}(\psi) = \sum_{i=1}^{N} \int_{\mathrm{RLag}_i(\psi)} (\|x - y_i\|^2 - \psi_i) d\rho(x) + \sum_{i=1}^{N} \alpha_i \psi_i \quad (2.10) $$

and

$$ \mathrm{RLag}_i(\psi) = \mathrm{Lag}_i(\psi) \cap \overline{B}_{y_i}(\sqrt{\psi_i}). \quad (2.11) $$

Note that strict positivity of each component of $\psi$ is a necessary condition for the restricted Laguerre cells to all be non-empty. Also, depending on the potential $\psi$, the restricted Laguerre cells do not necessarily (and will most often not) cover spt$\rho$. When there is ambiguity, we will refer to the $\mathrm{Lag}_i(\psi)$ as the *unrestricted* Laguerre cells. We now prove the existence of a maximizer via a standard concavity argument.
**Proposition 6.** The functional $\mathcal{K}: (0, \infty)^N \to \mathbb{R}$ is concave and of class $C^1$, with gradient

$$ \nabla \mathcal{K}(\psi) = \alpha - G(\psi), \qquad (2.12) $$

where $G_i = \rho(\mathrm{RLag}_i)$. Therefore, $\psi$ is a maximizer of $\mathcal{K}$ if and only if every re-stricted Laguerre cell has the same mass as its corresponding Dirac mass in the target measure:

$$ \psi \in \operatorname*{argmax}_{\mathcal{K}} \Leftrightarrow G_i(\psi) = \alpha_i, \forall i \in [1, N]. \qquad (2.13) $$

Moreover, there exists a maximizer.

*Proof.* Thanks to equations (2.10) and (2.11), we can write

$$ \mathcal{K}(\psi) = \sum_{i=1}^{N} \int_{\mathrm{Lag}_i(\psi)} \min\{\|x - y_i\|^2 - \psi_i, 0\} d\rho(x) + \sum_{i=1}^{N} \alpha_i \psi_i, \qquad (2.14) $$

and the second sum is linear in $\psi$ (with constant gradient $\alpha$) so we focus on the first sum. Let $\psi, \psi' \in (0, \infty)^N$. On $\mathrm{Lag}_i(\psi')$ we have $\min\{\|x - y_i\|^2 - \psi'_i, 0\} \le \min\{\|x - y_j\|^2 - \psi'_j, 0\}$ for all $j$, so since for each potential the (unrestricted) Laguerre cells are pairwise disjoint and cover spt$\rho$, we have

$$ \begin{aligned} & \sum_{i=1}^{N} \int_{\mathrm{Lag}_i(\psi')} \min\{\|x - y_i\|^2 - \psi'_i, 0\} d\rho(x) \\ & \le \sum_{j=1}^{N} \int_{\mathrm{Lag}_j(\psi)} \min\{\|x - y_j\|^2 - \psi'_j, 0\} d\rho(x) = \sum_{j=1}^{N} \int_{\mathrm{RLag}_j(\psi)} (\|x - y_j\|^2 - \psi'_j) d\rho(x) \\ & = \sum_{j=1}^{N} \int_{\mathrm{RLag}_j(\psi)} (\|x - y_j\|^2 - \psi_j) d\rho(x) + \sum_{j=1}^{N} \rho(\mathrm{RLag}_j(\psi))(\psi_j - \psi'_j). \end{aligned} $$

Adding the linear terms of $\mathcal{K}(\psi')$ yields $\mathcal{K}(\psi') \le \mathcal{K}(\psi) + \langle \alpha - G(\psi), \psi' - \psi \rangle$, where $G_i(\psi) = \rho(\mathrm{RLag}_i(\psi))$, and it follows that $\mathcal{K}: (0, \infty)^N \to \mathbb{R}$ is concave and that its superdifferential at $\psi$ contains $\alpha - G(\psi)$. Furthermore, by dominated convergence the $\rho$-measure of each restricted Laguerre cell is continuous in $\psi$, and we deduce that $\mathcal{K}$ is indeed differentiable with continuous gradient $\nabla \mathcal{K} = \alpha - G$. To prove the existence of a maximizer, take a maximizing sequence $(\psi^n)$ and suppose by contradiction that it is not bounded. Then, up to extraction of a subsequence, there exists an index $i$ such that $(\psi_i^n)$ monotonically converges to $-\infty$ or $+\infty$. In the first case, the corresponding restricted Laguerre cell is eventually empty, which is a contradiction; in the second case, the complement $\mathbb{R}^d \setminus \cup_{i=1}^N \mathrm{RLag}_i(\psi^n) = \mathbb{R}^d \setminus \cup_{i=1}^N \bar{B}_{y_i}(\sqrt{\psi_i^n})$ of the union of the cells has $\rho$-measure converging to zero, which is also a contradiction, since summing the optimality conditions (2.13) yields $1 - \sum_{i=1}^N G_i(\psi) = 1 - \sum_{i=1}^N \alpha_i > 0$. $\square$

**Remark 7.** Unlike in the balanced case, the functional of semi-discrete partial optimal transport is not invariant up to addition of an identical constant to each component of the potential. This is due to the target and source measures having different mass. However, one can retrieve this property by introducing a fictitious
point $y_\infty$ to which the inactive mass of $\rho$ will be sent, and by defining a new cost function

$$c(x, y_i) = \begin{cases} \|x - y_i\|^2 & \text{if } i \in [1, N], \\ 0 & \text{if } i = \infty. \end{cases}$$

In this setting, the potential has an additional component $\psi_\infty$, and the $i^{th}$ restricted Laguerre cell of the partial transport setting is now the $i^{th}$ unrestricted Laguerre cell for these new potential $\underline{\psi} = (\psi_1, \dots, \psi_N, \psi_\infty)$ and cost function $c$. Indeed, the constraint indexed by $j = \infty$ states that $\|x - y_i\|^2 - \psi_i \le 0 - \psi_\infty$, i.e. that $x$ is in the closed ball of center $y_i$ and of radius $\sqrt{\psi_i - \psi_\infty}$. The additional Laguerre cell is sent to the auxiliary point $y_\infty$, and is the closure of the complement of the $N$ open balls $B_{y_i}(\sqrt{\psi_i - \psi_\infty})$. The reason for which the dual functional $K$ in (2.10) — or equivalently in (2.14) — is not invariant up to addition of the same constant to each component of the potential is that we implicitly set $\psi_\infty = 0$ in the balanced transport formulation that we have just described. Note that in this formulation, the additional Laguerre cell is not convex, and in fact not even connected.

In this article, we will mainly deal with the one-dimensional case. Let us emphasize that semi-discrete partial transport in 1D is somewhat degenerate compared to higher dimensions. This is due to the fact that the boundary of each restricted Laguerre cell contains only finitely many points — two points exactly, the left and right extremities — and that, in some cases, their velocity is not linear in the (arbitrarily small) variation of the potential. As will be apparent in section 2.5, this is due to the fact that these extremities are respectively defined as a maximum and a minimum, and that the maximum or minimum of two smooth quantities is in general not smooth when the latter coincide. We will thus be lead to consider the unilateral partial derivatives of each cell's boundaries.

**Theorem 8.** Fix $d \ge 2$ and suppose that $\rho$ is compactly supported and absolutely continuous, with continuous density on the support. In addition, suppose that the intersection of $\text{spt}\rho$ with any sphere centered at one of the $y_i$ is $\rho H^{d-1}$-negligible, as well as the intersection of $\text{spt}\rho$ with any hyperplane intersecting its interior and orthogonal to some segment $[y_i, y_j]$. Then the functional $K$ is $C^2$ on the convex open set $\mathcal{D}$ of potentials such that each restricted Laguerre cell has strictly positive $\rho$-measure, as well as the complement of their union. The second derivatives are given by

$$\partial_{\psi_i, \psi_j}^2 K(\psi) = \frac{1}{2\|y_i - y_j\|} \int_{\text{RLag}_i(\psi) \cap \text{RLag}_j(\psi)} \rho(x) d\mathcal{H}^{d-1}(x) \quad (2.15)$$

for $i \ne j$, and

$$\partial_{\psi_i, \psi_i}^2 K(\psi) = -\frac{1}{2\sqrt{\psi_i}} \int_{\text{RLag}_i(\psi) \cap S_{y_i}(\sqrt{\psi_i})} \rho(x) d\mathcal{H}^{d-1}(x) - \sum_{j \ne i} \partial_{\psi_i, \psi_j}^2 K(\psi). \quad (2.16)$$

In particular, $K$ is strictly concave on $\mathcal{D}$ and the optimal potential is unique.

**Remark 9.** The assumptions on the support of $\rho$ guarantee for instance that $\partial_{\psi_i, \psi_j}^2 K$ does not jump when the spherical part of the boundary of $\text{RLag}_i$ crosses the support.
A simple case for which the assumptions on sptρ are met is when the support of ρ is a bounded convex polygon.

In dimension $d \ge 2$, the Kantorovich functional (2.10) is of class $C^2$ whenever the density $\rho$ has a bounded polytope as its support and is continuous on this set, while it is only of class $C^1$ in the unidimensional setting. Before dealing with the case of dimension one, let us state the precise result we derived for the multidimensional case, the proof of which is in Appendix A.

*Edges of a cell.* The key notion to understand (at least conceptually) the $C^2$ regularity of the dual functional in dimension greater than one is that of *edges* of a (restricted Laguerre) cell. We call *edge* of a cell any intersection of the cell with two other cells, or with one other cell and the inactive part (the auxiliary cell introduced in Remark 7). The edges of a given cell are precisely the points of its boundary whose respective velocities are ill-defined. Indeed, if a point belongs to an edge, then the two (or more) facets it belongs to have different velocity vectors as we vary the potential $\psi$. Roughly speaking, in dimension $d \ge 2$, the differentiability of the mass of a cell with respect to the potential is due to the fact that its edges are of dimension $d-2$, and are thus not “seen” by the $(d-1)$-dimensional Hausdorff measure.

2.5. **The unidimensional case.** As mentioned above, in the one-dimensional case the gradient of $\mathcal{K}$ is not differentiable. However, it does have a weaker kind of regularity, as stated in the following lemma, the proof of which is reported to Appendix B. For convenience, we assume (without loss of generality) that the positions of the Dirac masses are labeled in increasing order: $y_1 < y_2 < \cdots < y_N$. In this case, provided every (unrestricted) Laguerre cell has non-empty interior, we have $\text{Lag}_i = [z_{i-1}, z_i]$, where the positions

$$z_i(\psi) = \frac{y_i + y_{i+1}}{2} - \frac{\psi_{i+1} - \psi_i}{2(y_{i+1} - y_i)} \quad (2.17)$$

are the boundaries between consecutive Laguerre cells. Naturally, we set $z_0 = -\infty$ and $z_N = +\infty$. The restricted Laguerre cells can then be written $\text{RLag}_i = [a_i, b_i]$, where

$$\begin{cases} a_i(\psi) = \max\{z_{i-1}(\psi), y_i - \sqrt{\psi_i}\}, \\ b_i(\psi) = \min\{z_i(\psi), y_i + \sqrt{\psi_i}\}. \end{cases} \quad (2.18)$$

We will make the following assumption on the source measure.

**Assumption 1.** The probability measure $\rho \in \mathcal{P}(\mathbb{R})$ is absolutely continuous, its support is a bounded interval, and on this interval its density $\rho$ is both continuous and bounded away from zero and infinity, i.e. there exists finite strictly positive constants $\rho_{\min}, \rho_{\max}$ such that $\rho_{\min} \le \rho(x) \le \rho_{\max}$ for every $x \in \text{spt}\rho$.

**Lemma 10.** Suppose that $\rho$ satisfies Assumption 1. Denote $\mathcal{D}$ the open set of potentials for which every restricted Laguerre cell has nonzero $\rho$-measure, as well as the complement of their union. At each $\psi \in \mathcal{D}$, the function $G$ has well-defined unilateral directional derivatives $\partial_v^+ G(\psi) := \lim_{t \to 0^+} (G(\psi+tv) - G(\psi))/t$ with respect to any vector $v \in \mathbb{R}^N$. Moreover, any such derivative can be written $\partial_v^+ G(\psi) = H(\psi, v)v$ for some tridiagonal, weakly diagonally dominant matrix $H(\psi, v)$.
Each $H(\psi, v)$ can in addition be chosen to be symmetric, so that it is block diagonal, with all blocks being symmetric, tridiagonal, irreducible, weakly diagonally dominant matrices; see Lemma 29 in Appendix B. At last, we show in Lemma 31 of the same section that each of these blocks has at least one strictly diagonally dominant coefficient, and the matrix $H(\psi, v)$ is thus non-singular by Taussky's theorem, see for instance [13, Corollary 6.2.27]. This allows for the following uniqueness result.

**Proposition 11.** If $\rho$ satisfies Assumption 1, then the dual functional $K: (0, \infty)^N \to \mathbb{R}$ admits a unique maximizer $\psi^*$, and this maximizer naturally belongs to the set $D$ defined in Lemma 10.

*Proof.* It is straightforward to see that any maximizer must be in $D$. Suppose that $\psi^{(1)}$ and $\psi^{(2)}$ are two distinct maximizers. Then, by concavity, every convex combination of these is also a maximizer, so setting $v = \psi^{(2)} - \psi^{(1)}$ we have $\partial_v^+ G(\psi^*) = 0$. But the latter quantity writes $H(\psi^{(1)}, v)v$, where $H(\psi^{(1)}, v)$ is a non-singular matrix, so $v$ must be zero and we have $\psi^{(1)} = \psi^{(2)}$. $\square$

**Example 12.** A concrete example of non-differentiability can be constructed as follows. Take two Dirac masses, respectively at $y_1 = 0$ and $y_2 = 1$, and $\rho$ the Lebesgue measure on $[0, 1]$. The two unrestricted Laguerre cells intersect at $z_1(\psi) = \frac{1}{2}(1 - \psi_2 + \psi_1)$. Consider the potential $\psi^0 = (\frac{1}{4}, \frac{1}{4})$, for which $z_1(\psi^0)$ coincides with $\sqrt{\psi_1^0} = \frac{1}{2}$, and let $\psi^t = \psi^0 + t(1, 0)$, $t \in \mathbb{R}$. The derivatives of $z_1(\psi^t)$ and $\sqrt{\psi_1^t}$ at $t = 0$ are respectively $\frac{1}{2}$ and 1. Thanks to equation (2.18), at $t = 0$ the quantity $b_1(\psi^t)$ has left-hand derivative $-1$ and right-hand derivative $\frac{1}{2}$. Since $a_1 \equiv 0$ and $\rho \equiv 1$ on $[0, 1]$, we deduce that $G_1(\psi^t)$ has the same left-hand and right-hand derivatives at $t = 0$ as $b_1(\psi^t)$, which shows that $G$ is not differentiable at $\psi^0$.

In order to numerically solve (balanced) semi-discrete optimal transport, Mérigot et al. [17] introduced a damped Newton algorithm; in particular they use a variable step size approach (backtracking line search) to make sure that, throughout the descent, the Laguerre cells all have a mass greater than a given threshold. To adapt this algorithm to the partial version of the problem, we have added a natural extra constraint: the complement of the union of the Laguerre cells must also have mass greater than the chosen threshold. Unfortunately, the resulting algorithm is not as efficient as in the balanced case: the convergence crucially depends on the initial point and the gradient descent often involves long plateau phases, with the algorithm occasionally getting stuck. These difficulties are likely linked to the non-$C^2$ nature of the dual functional. In order to avoid this kind of issues, we introduce a regularization of the partial optimal transport problem, which we detail in the section below.

### 3. REGULARIZATION OF SEMI-DISCRETE PARTIAL TRANSPORT IN DIMENSION ONE

The partial optimal transport problem (2.7) can be written

$$ \inf_{\gamma \in \Gamma_{\max}(\rho, \mu)} \int_{\mathbb{R} \times \mathbb{R}} |x-y|^2 d\gamma(x,y) = \inf_{\sigma \in \mathcal{M}_+(\mathbb{R})} W_2^2(\sigma, \mu) + F_0(\sigma), \quad (3.1) $$
where the functional $F_0: \mathcal{M}_+(\mathbb{R}) \to \mathbb{R} \cup \{+\infty\}$ defined by

$$F_0(\sigma) = \begin{cases} \int_{\mathbb{R}} \chi_{[0,1]}\left(\frac{d\sigma}{d\rho}\right) d\rho & \text{if } \sigma \ll \rho, \\ +\infty & \text{otherwise,} \end{cases}$$

ensures that $\sigma$ is dominated by $\rho$. Here, $\chi_{[0,1]}$ denotes the characteristic function of the interval $[0, 1]$ in the sense of convex analysis. In this formulation, the optimal $\sigma$ corresponds to the active part of $\rho$, that is, the part of $\rho$ that actually gets sent to $\mu$ in the partial transport problem. In the semi-discrete framework, the optimal $\sigma$ is simply the restriction of $\rho$ to the union of the restricted Laguerre cells.

**3.1. The regularized functional.** In order to regularize Problem (3.1), we may replace $F_0$ by a functional $F$ of the form

$$F(\sigma) = \begin{cases} \int_{\mathbb{R}} f\left(\frac{d\sigma}{d\rho}\right) d\rho & \text{if } \sigma \ll \rho, \\ +\infty & \text{otherwise,} \end{cases}$$

where $f: \mathbb{R} \to [0, +\infty]$ is a convex function approximating $\chi_{[0,1]}$ (in a sense which will be clarified below), and try to solve

$$\inf_{\sigma \in \mathcal{M}_+(\mathbb{R})} W_2^2(\sigma, \mu) + F(\sigma). \quad (3.2)$$

**Assumption 2.** Throughout the rest of this subsection, we make the following assumptions on $f$, which guarantee that $F$ be both convex and lower semi-continuous for the weak convergence of measures, and that the domain of $F$ is included in the set of positive measures. The fact that $f$ is strictly convex, with domain $[0, 1]$ and subderivatives bounded below by 1 guarantees that $f^*$ has full domain and is convex $C^1$, non-decreasing, superlinear, with non-decreasing derivative upper bounded by 1.

(1) $f: \mathbb{R} \to [0, +\infty]$ is proper, lower semi-continuous, and strictly convex on its domain.

(2) $\text{dom} f = [0, 1]$.

(3) The subderivatives of $f$ are all at least 1.

By classical arguments, the dual problem of (3.2) reads $\sup_{\psi \in \mathbb{R}^N} \mathcal{K}^f(\psi)$, where the regularized functional $\mathcal{K}^f: \mathbb{R}^d \to \mathbb{R}$ is defined by

$$\mathcal{K}^f(\psi) = - \sum_{i=1}^{N} \int_{\text{Lag}_i(\psi)} f^*(\psi_i - |x - y_i|^2) d\rho(x) + \sum_{i=1}^{N} \alpha_i \psi_i. \quad (3.3)$$

Moreover, the following proposition — which holds in particular when $\rho$ and $f$ respectively satisfy Assumption 1 and Assumption 2 — states that the regularized functional $\mathcal{K}^f$ is concave and indeed twice continuously differentiable on a large domain containing its maxima.

**Proposition 13.** Suppose that $\rho$ is an absolutely continuous probability measure and that $f$ satisfies Assumption 2. Then, the functional $\mathcal{K}^f$ is both concave and of class $C^1$ on $\mathbb{R}^N$, with gradient

$$\nabla \mathcal{K}^f(\psi) = \alpha - G^f(\psi), \quad (3.4)$$
where the quantity

$$G_i^f(\psi) = \int_{\text{Lag}_i(\psi)} (f^*)'(\psi_i - |x - y_i|^2) d\rho(x) \quad (3.5)$$

will be referred to as the regularized mass of the $i^{th}$ Laguerre cell.

If in addition $\rho$ has compact, convex support, with continuous density on this set, and $(f^*)'$ is piecewise $C^1$, then $K^f$ is of class $C^2$ on the open set $\tilde{D}$ of potentials such that all Laguerre cells have strictly positive $\rho$-measure.. On this set, the Hessian matrix $D^2K^f(\psi)$ is (symmetric) tridiagonal with coefficients

$$[D^2K^f(\psi)]_{i,i+1} = \frac{1}{2(y_{i+1} - y_i)} (f^*)'(\psi_i - |z_i(\psi) - y_i|^2) \rho(z_i(\psi))$$

and

$$[D^2K^f(\psi)]_{i,i} = - \int_{\text{Lag}_i(\psi)} (f^*)''(\psi_i - |x - y_i|^2) d\rho(x) - [D^2K^f(\psi)]_{i,i-1} - [D^2K^f(\psi)]_{i,i+1}.$$

**Remark 14.** Since the domain of $f$ is the set of non-negative numbers and $f$ is increasing on this set, we have $f^*(t) \ge -f(0)$ for every $t \in \mathbb{R}$, with equality if $t \le 0$. The regularizing functions we will define in section 3.2 will be zero at $x = 0$, so their respective convex conjugates will be non-negative, with support in $[0, +\infty)$.

*Proof.* The second sum in the right-hand side of (3.3) is linear in $\psi$, so we only consider the first sum, which writes $\int_{\mathbb{R}} f^*(\max_{1 \le i \le N} \psi_i - |x - y_i|^2) d\rho(x)$. Since $f^*$ is continuous and $\rho$ has compact support, this quantity is differentiable thanks to dominated convergence, and its first derivatives are given by (3.5), Continuity of these derivatives is a direct consequence of the dominated convergence theorem. To check that $K^f$ is convex, we conveniently write the first integral as $\int_{\mathbb{R}} f^*(-\psi^c) d\rho$. The $c$-transform is convex in $\psi$ as a minimum of affine linear functions, so since $f^*$ is convex and non-decreasing, the integral in question is a convex function of $\psi$, and concavity of $K^f$ follows.

Let us now suppose that $(f^*)'$ is piecewise $C^1$. Thanks to Leibniz integral rule, we get the desired expression for the second derivatives of $K^f$, and these are continuous thanks to the additional assumption that the density $\rho$ is continuous on its support. □

## 3.2. A family of regularizing functions.
Let us now introduce a convenient family of functions $(f_\varepsilon)_{\varepsilon>0}$ approximating $\chi_{[0,1]}$ and satisfying all the assumptions in Proposition 13. Define $f_\varepsilon(t) = \varepsilon^2 f(t)$, where

$$f(t) = \begin{cases} \frac{1}{3}t^3 & \text{if } 0 \le t \le 1, \\ +\infty & \text{otherwise,} \end{cases} \quad (3.6)$$

so that, taking the derivatives to be zero on the negative half-line,

$$f^*(t) = \begin{cases} \frac{3}{2}t^{\frac{2}{3}} & \text{if } 0 < t < 1, \\ t - \frac{1}{3} & \text{if } t > 1, \\ 0 & \text{if otherwise,} \end{cases} \quad (3.7)$$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Graphs of $(f_{\varepsilon}^*)'$ for different values of $\varepsilon$, as well as the limit graph when $\varepsilon \to 0$, which corresponds to the unregularized problem.</figcaption></figure>

$$ (f^*)'(t) = \begin{cases} \min\{\sqrt{t}, 1\} & \text{if } t > 0, \\ 0 & \text{otherwise,} \end{cases} \qquad (f^*)''(t) = \begin{cases} \frac{1}{2\sqrt{t}} & \text{if } 0 < t < 1, \\ 0 & \text{otherwise.} \end{cases} \qquad (3.8) $$

The derivatives of $f_{\varepsilon}$ then read

$$ f_{\varepsilon}^{*}(t) = \varepsilon^{2} f^{*}(\varepsilon^{-2}t), \quad (f_{\varepsilon}^{*})'(t) = (f^{*})'(\varepsilon^{-2}t), \quad (f_{\varepsilon}^{*})''(t) = \varepsilon^{-2}(f^{*})''(\varepsilon^{-2}t). $$

Figure 1 shows the graph of $(f_{\varepsilon}^*)'$ for different values of $\varepsilon$, as well as their pointwise limit $1_{(0,+\infty)}$ as $\varepsilon \to 0$. Note that replacing $(f^*)'$ by this pointwise limit in the right-hand side of (3.5) yields the mass of the restricted Laguerre cell (with respect to $\rho$). Similarly, taking $f^* = \lim_{\varepsilon \to 0} f_{\varepsilon}^* = 1_{\ge 0}$Id in (3.3), we recover the unregularized functional for the multi-dimensional partial transport problem.

*Notations.* For convenience, we write write $K^{\varepsilon}$, $G^{\varepsilon}$ for $K^{f^{\varepsilon}}$, $G^{f^{\varepsilon}}$ respectively.

**3.3. Two-dimensional interpretation.** The choice of $f^{\varepsilon}$ was motivated by the fact that the corresponding regularized problem is a two-dimensional version of the original one. Indeed, define $\rho^{\varepsilon}$ the probability measure on $\mathbb{R}^2$ with density $\rho^{\varepsilon}(x) = (2\varepsilon)^{-1}\rho(x^1)1_{|x^2|<\varepsilon}$ for $x = (x^1, x^2)$, whose support is the rectangle $\text{spt}\rho \times [-\varepsilon, \varepsilon]$. The discrete measure $\mu$ is naturally extended to $\mathbb{R}^2$ by identifying $y_i$ to $z_i = (y_i, 0)$. Since the $z_i$ are all aligned on the first axis, a straightforward computation shows that the Laguerre cells of the two-dimensional problem are the vertical strips $\text{Lag}_i^{(2)}(\psi) = \text{Lag}_i(\psi) \times \mathbb{R}$ and that the corresponding restricted Laguerre cells are $\text{RLag}_i^{(2)}(\psi) = (\text{Lag}_i(\psi) \times \mathbb{R}) \cap \bar{B}_{z_i}^{(2)}(\sqrt{\psi_i})$. As a result, we can write $\text{RLag}_i^{(2)}(\psi)$ as the following disjoint union of vertical segments

$$ \text{RLag}_i^{(2)}(\psi) = \bigsqcup_{x^1 \in \text{Lag}_i(\psi)} \{x^1\} \times \bar{B}_0(\sqrt{\psi_i - |x^1 - y_i|^2}). $$

Since $\rho^{\varepsilon}$ is concentrated on the horizontal strip $\mathbb{R} \times (-\varepsilon, \varepsilon)$, Fubini's theorem yields

$$ \begin{align*} \rho^{\varepsilon}(\text{RLag}_i^{(2)}(\psi)) &= \frac{1}{2\varepsilon} \int_{\text{Lag}_i(\psi)} 2 \min\{\sqrt{\psi_i - |x^1 - y_i|^2}, \varepsilon\} \rho(x^1) dx^1 \\ &= \int_{\text{Lag}_i(\psi)} \min\{\varepsilon^{-2}(\psi_i - |x^1 - y_i|^2), 1\} d\rho(x^1) \\ &= \int_{\text{Lag}_i(\psi)} (f_{\varepsilon}^*)'(\psi_i - |x^1 - y_i|^2) d\rho(x^1), \end{align*} $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Visualization of the two-dimensional problem corresponding to the regularization. The long gray rectangle is the support of $\rho^\epsilon$ and the colored areas are the two-dimensional restricted Laguerre cells.</figcaption></figure>

and this is precisely the regularized mass of the one-dimensional Laguerre cell. Figure 2 illustrates this two-dimensional interpretation of the regularized problem.

**Remark 15.** Since the Laguerre of the two-dimensional problem are of the form $\text{Lag}_i^{(2)}(\psi) = \text{Lag}_i(\psi) \times \mathbb{R}$, they form a tessellation the plan into $N$ vertical strips, whose corresponding indices are increasing from left to right.

We can, then, establish the following convergence result.

**Proposition 16.** The functional $-K^\epsilon \Gamma$-converges to $-K$ as $\epsilon \to 0$.

*Proof.* Suppose that $\epsilon_n \xrightarrow{n \to \infty} 0$ and $\psi^n \xrightarrow{n \to \infty} \psi$, and denote $f_n := f_{\epsilon_n}$. By dominated convergence, we have

$$ \int_{\mathbb{R}} \mathbb{1}_{\text{Lag}_i(\psi^n)}(x) (f_n^*)'(\psi_i^n - |x-y_i|^2) \rho(x) dx \xrightarrow{n \to \infty} \int_{\mathbb{R}} \mathbb{1}_{\text{Lag}_i(\psi)}(x) \mathbb{1}_{\{|x-y_i|^2 \le \psi_i\}}(x) \rho(x) dx. $$

Indeed, the integrand is dominated by the integrable density $\rho$, the set $\{x \in \mathbb{R} : |x - y_i|^2 = \psi_i\}$ is negligible, and for $x$ such that $|x - y_i|^2 \ne \psi_i$ we have

$$ (f_n^*)'(\psi_i^n - |x - y_i|^2) = (f^*)'(\epsilon_n^{-2}(\psi_i^n - |x - y_i|^2)) \\ \rightarrow \begin{cases} \lim_{+\infty} f' = 1 & \text{if } |x - y_i|^2 < \psi_i, \\ \lim_{-\infty} f' = 0 & \text{if } |x - y_i|^2 > \psi_i. \end{cases} \quad \square $$

## 4. AN ODE CHARACTERIZATION OF THE REGULARIZED PROBLEM

The maximizer $\psi^\epsilon$ of $K^\epsilon$ is characterized by $\nabla K^\epsilon(\psi^\epsilon) = 0$, that is, $G^\epsilon(\psi^\epsilon) = \alpha$, where we recall that $G^\epsilon(\psi)$ is the vector of regularized masses of the Laguerre cells corresponding to $\psi$. Thanks to these optimal conditions, the implicit function theorem yields an ODE satisfied by $\epsilon \mapsto \psi^\epsilon$ which, combined with some estimates on the derivatives of $G^\epsilon$, allows to quantify the convergence of $\psi^\epsilon$ to the solution $\psi^*$ of the unregularized problem. For every $\epsilon > 0$, the uniqueness of $\psi^\epsilon$ is guaranteed by Theorem 8 and Assumption 1.

**Theorem 17.** Suppose that $\rho$ satisfies Assumption 1. Then there exists $\epsilon_0 > 0$ such that for any $\epsilon \in (0, \epsilon_0)$, the regularized functional $K^\epsilon$ has a unique maximizer $\psi^\epsilon \in \mathbb{R}^N$. The function $\epsilon \mapsto \psi^\epsilon$ is of class $C^1$ on $(0, \epsilon_0)$, and on this interval we have

$$ [D_\psi G^\epsilon(\psi^\epsilon)] \dot{\psi}^\epsilon + (\partial_\epsilon G^\epsilon)(\psi^\epsilon) = 0, \quad (4.1) $$
where $DG^\epsilon = -D^2K^\epsilon$. Yet again on this interval, it holds

$$ \|\psi^\epsilon - \psi^*\| \lesssim \epsilon^2, \qquad (4.2) $$

where $\psi^* \in (0, +\infty)^N$ is the maximizer of the unregularized functional $K$ and where $\lesssim$ hides a constant depending only on $\rho$, on the $y_i$ and on the $\alpha_i$.

Moreover, this behavior is widely observed in numerical tests, as illustrated in Figure 8. The following example shows that the quadratic convergence rate obtained in this section is tight.

**Example 18.** To see that the quadratic rate given in Theorem 17 is tight, consider the simple case where $\rho$ is the uniform probability measure on $[0, 1]$ and where the discrete measure is $\mu = \frac{1}{2}\delta_0$. The optimal $\psi \in \mathbb{R}^1$ for the unregularized problem is $\psi^* = \frac{1}{4}$, and the optimal regularized mass of the Laguerre cell for $\epsilon > 0$ writes

$$ \frac{1}{2} = G^\epsilon(\psi^\epsilon) = \frac{1}{2}\sqrt{\psi^\epsilon - \epsilon^2} + \frac{\psi^\epsilon}{2\epsilon} \arcsin\left(\frac{\epsilon}{\sqrt{\psi^\epsilon}}\right) = \sqrt{\psi^\epsilon} - \frac{\epsilon^2}{6\sqrt{\psi^\epsilon}} + o(\epsilon^2) $$

asymptotically as $\epsilon \to 0^+$, since $\sqrt{\psi^\epsilon} \to \frac{1}{2}$. We thus have $\sqrt{\psi^\epsilon} - \frac{1}{2} \sim \frac{\epsilon^2}{3}$ and so

$$ \psi^\epsilon - \psi^* = \left(\sqrt{\psi^\epsilon} + \frac{1}{2}\right) \left(\sqrt{\psi^\epsilon} - \frac{1}{2}\right) \sim \frac{\epsilon^2}{3}. $$

**4.1. Notations.** In this subsection, we exploit the 2D interpretation of the regularization introduced above, so the notations $\text{Lag}_i$ and $\text{RLag}_i$ refer respectively to the two-dimensional unrestricted and restricted Laguerre cells. In consideration of Remark 7, the complement of $\cup_{i=1}^N \text{RLag}_i$ will be denoted $\text{RLag}_\infty$. To further lighten the notation, we also introduce the following definitions for the quantities related to the optimal potential $\psi^\epsilon$. For $i \in [1, N] \cup \{\infty\}$, $\text{RLag}_i^\epsilon$ will denote the intersection of $\text{RLag}_i(\psi^\epsilon)$ with the rectangle $\text{spt}\rho_\epsilon = \text{spt}\rho \times [-\epsilon, \epsilon]$, and for $i \neq \infty$, $B_i^\epsilon$ will be the open disk of center $(y_i, 0)$ and of radius $\sqrt{\psi_i^\epsilon}$. The Hessian matrix $-DG^\epsilon(\psi^\epsilon) = -D^2K^\epsilon(\psi^\epsilon)$ and the mixed derivative $-(\partial_\epsilon G^\epsilon)(\psi^\epsilon) = -\partial_\epsilon(\nabla K^\epsilon)(\psi^\epsilon)$ will respectively be denoted $H^\epsilon$ and $q^\epsilon$, so that equation (4.1) writes $H^\epsilon \dot{\psi}^\epsilon = q^\epsilon$. Note that there is a minus sign in front of both quantities, so that in particular $H^\epsilon$ is positive semi-definite. At last, we denote $\alpha_{\min} = \min\{\alpha_1, \dots, \alpha_N\}$ and $\alpha_\infty = 1 - \sum_{i=1}^N \alpha_i$.

**4.2. Proof of the theorem.** In order to prove Theorem 17, we first need to control the minimum eigenvalue of $H^\epsilon$ — which is a positive semi-definite matrix since $K^\epsilon$ is concave — independently of $\epsilon$, which is the main difficulty. The argument is based on results of spectral graph theory, and we thus introduce the concept of Laplacian matrices. Only simple (i.e. without loops), *undirected* graphs will be considered throughout the paper.

**Definition 19.** A Laplacian matrix is any symmetric matrix $M = (M_{i,j})_{0 \le i,j \le N}$ with non-negative off-diagonal entries such that each line sums to zero. It corresponds to the unique weighted graph with vertices $\{1, \dots, N, \infty\}$ and weights $w_{i,j} = |M_{i,j}|$, where for convenience we identify the indices 0 and $\infty$. Two vertices i and j are called neighbors if $|M_{i,j}| > 0$, and we say the matrix $M$ connected if the corresponding graph is connected.
In this work, the indices of a square matrix with $N + 1$ rows will always range from 0 to $N$, whereas they will range from 1 to $N$ for a square matrix with $N$ rows. As will be apparent in the definition of equation (4.4), in our case the index zero corresponds to the fictitious point that was the subject of Remark 7. We therefore sometimes use the symbol $\infty$ for this index. The well known properties given in the following lemma can be found in Spielman's book [21, Chapter 3].

**Lemma 20.** Let $M \in \mathbb{R}^{(N+1) \times (N+1)}$ be a Laplacian matrix. Then $M$ is symmetric positive semi-definite and its smallest eigenvalue is 0, with $1$ (the vector of ones) an associated eigenvector. We thus denote $0 = \lambda_0(M) \le \lambda_1(M) \le \cdots \le \lambda_N(M)$ its eigenvalues with multiplicities. The second value is given by

$$ \lambda_1(M) = \min_{\substack{\|v\|=1 \\ v \perp 1}} v^T M v, \qquad (4.3) $$

and is non-zero if and only if the weighted graph associated to $M$ is connected.

Because of the correspondence mentioned above between Laplacian matrices and weighted graphs on $\{1, \dots, N, \infty\}$, for any such graph $G$ we denote $0 = \lambda_0(G) \le \lambda_1(G) \le \cdots \le \lambda_N(G)$ the eigenvalues with multiplicities of its Laplacian matrix. The following lemma states that the second smallest eigenvalue is non-decreasing in the weights of the graph.

**Lemma 21.** A Laplacian matrix $M$ is said to dominate another Laplacian matrix $M'$ if $|M_{i,j}| \ge |M'_{i,j}|$ for all distinct indices $i \ne j$. If this is the case, then $\lambda_1(M) \ge \lambda_1(M')$.

*Proof.* Suppose $M$ dominates $M'$. Then $M - M'$ is semi-definite positive, so writing $M = M' + (M - M')$ and distributing the minimum in (4.3) yields the desired inequality. $\square$

Now, let $M^\epsilon \in \mathbb{R}^{(N+1) \times (n+1)}$ be the tridiagonal, symmetric matrix defined by

$$ \left\{ \begin{array}{ll} M^\epsilon_{i,i+1} = \frac{-1}{2(y_{i+1} - y_i)} \int_{RLag^\epsilon_i \cap RLag^\epsilon_{i+1}} \rho^\epsilon(x) d\mathcal{H}^1(x) & \text{for } i \in [1, N-1], \\ M^\epsilon_{i,0} = \frac{-1}{2\sqrt{\psi_i^\epsilon}} \int_{RLag^\epsilon_i \cap \partial B^\epsilon_i} \rho^\epsilon(x) d\mathcal{H}^1(x) & \text{for } i \in [1, N], \\ M^\epsilon_{i,i} = -\sum_{j \ne i} M^\epsilon_{i,j} & \text{for } i \in [0, N], \end{array} \right. \qquad (4.4) $$

with all other coefficients null, and call $G^\epsilon$ the corresponding weighted graph on $\{1, \dots, N, \infty\}$. Indeed, it is straightforward to check that $M^\epsilon$ is a Laplacian matrix. Note that by Theorem 8 and Proposition 13, the submatrix obtained by removing the line and the column of index 0 — which correspond the auxiliary point — is precisely the tridiagonal matrix $H^\epsilon$, the minimal eigenvalue of which we seek to bound from below. Thanks to the following lemma, it suffices to find a lower bound for $\lambda_1(M^\epsilon)$.

**Lemma 22.** Let $M = (M_{i,j})_{0 \le i,j \le N}$ be a Laplacian matrix, and let $\tilde{M}$ be the submatrix obtained by removing the line and the column of index 0. Then $\lambda_{\min}(\tilde{M})$, the
smallest eigenvalue of $\tilde{M}$, satisfies

$$ \lambda_{\min}(\tilde{M}) \geq \frac{\lambda_1(M)}{N+1}. $$

*Proof.* Denote $P$ the orthogonal projection from $\mathbb{R}^{N+1}$ onto $1^{\perp} \subseteq \mathbb{R}^{N+1}$, let $u \in \mathbb{R}^N$ and define $\underline{u} = (0, u_1, \dots, u_N) \in \mathbb{R}^{N+1}$. Since $M1 = 0$, we have

$$ u^T \tilde{M} u = \underline{u}^T M \underline{u} = (P \underline{u})^T M (P \underline{u}) \geq \lambda_1(M) \|P \underline{u}\|^2. \quad (4.5) $$

Moreover, using the additional notation $e = (N+1)^{-1/2}1 \in \mathbb{S}^{N+1}$, we find

$$ \|P\underline{u}\|^2 = \|\underline{u} - (e^T \underline{u})e\|^2 = \|\underline{u}\|^2 - (e^T \underline{u})^2 $$

and $(e^T \underline{u})^2 = (N+1)^{-1}(\sum_{i=1}^N u_i)^2 \leq N(N+1)^{-1}\|u\|^2$ by Cauchy-Schwarz inequality, so $\|P\underline{u}\|^2 \geq (N+1)^{-1}\|u\|^2$. Inequality (4.5) and the characterization of the smallest eigenvalue of $\tilde{M}$ as the minimum of the Rayleigh quotient $\frac{u^T \tilde{M} u}{\|u\|^2}$ over all non-zero vectors $u$ yields the desired inequality. $\square$

Our strategy for getting a lower bound on $\lambda_1(M^\epsilon)$ which is independent of $\epsilon$ is to compare the Laplacian matrix $M^\epsilon$ to that of a uniformly weighted graph. In order to do so, we need to control from below the absolute values of the non-zero off-diagonal entries $|M^\epsilon_{i,i+1}|$ and $|M^\epsilon_{i,0}|$. This in turn requires us to control the components of $\psi^\epsilon$. The object of the following lemma is to show that for $\epsilon$ small enough, these components are uniformly bounded away from zero and infinity. Roughly speaking, the radii $\sqrt{\psi^\epsilon_i}$ of the balls $B^\epsilon_i$ cannot be two small because the cells have mass bounded below and because the density is upper bounded; they cannot be too large either because the auxiliary cell is not empty.

**Lemma 23.** If $0 < \epsilon < 1$, then for all $i \in [1, N]$ we have

$$ r \leq \sqrt{\psi_i^\epsilon} \leq R, $$

where $r = \frac{\alpha_{\min}}{2\rho_{\max}}$ and $R = (1 + \text{diam}(\text{spt}\rho \cup \{y_1, \dots, y_N\})^2)^{\frac{1}{2}}$. Moreover, for any $0 < \epsilon < r$, the cells $\text{RLag}_i^\epsilon$ all intersect the horizontal lines of ordinate $\pm\epsilon$. In consideration of these two properties, we define $\epsilon_0 = \min(1, r)$.

*Proof.* We first prove the lower bound on $\sqrt{\psi_i^\epsilon}$. To do so, we note that $\text{RLag}_i^\epsilon$ is included in a rectangle of width $2\epsilon$ and of height $2\sqrt{\psi_i^\epsilon}$. Since $G_i^\epsilon(\psi^\epsilon) = \alpha_i$ and since $\rho^\epsilon$ is bounded above by $\frac{\rho_{\max}}{2\epsilon}$, we get $\alpha_i \leq \frac{\rho_{\max}}{2\epsilon} \times 2\epsilon \times 2\sqrt{\psi_i^\epsilon}$, which yields the lower bound. We now derive the upper bound. The cell $\text{RLag}_\infty^\epsilon$ associated to the auxiliary point has non-zero $\rho^\epsilon$-measure $\alpha_\infty > 0$, so it must intersect the support of $\rho^\epsilon$. Fix $(x, w)$ a point in this intersection. It satisfies $\psi_i^\epsilon \leq |x - y_i|^2 + |w - 0|^2$ for all $i \in [1, N]$. The first square in the left-hand side is bounded by the squared diameter of $\text{spt}\rho \cup \{y_1, \dots, y_N\}$, while the second one is bounded by $\epsilon^2$, and the upper bound follows.

Let us now show that for $\epsilon$ small enough, the optimal cells all intersect the horizontal lines $\mathbb{R} \times \{\pm\epsilon\}$. Fix $\epsilon \in (0, r)$ and suppose by contradiction that, for some $i \in [1, N]$, the cell $\text{RLag}_i^\epsilon$ does not intersect these horizontal lines. Then the $\rho^\epsilon$-measure of the cell is bounded above by $\frac{\rho_{\max}}{2\epsilon}$ times the area of the two disk segments shown in
<figure><img src="image_4.png" /><figcaption>FIGURE 3. The area of each of the two disk segments in blue is $\omega^2(\arcsin(\varepsilon/\omega) - \varepsilon/\omega\sqrt{1-(\varepsilon/\omega)^2})$. A simple computation shows that for $u \in [0, 1]$ we have $\arcsin(u) - u\sqrt{1-u^2} \le \pi u^3/2$, so that the latter area is at most $\pi\varepsilon^3/(2\omega)$. A rougher estimate consists in applying $1-u^2 \le \sqrt{1-u^2}$ to $u = \varepsilon/\omega$, which yields the upper bound $2\varepsilon^3/\omega$ for the area.</figcaption></figure>

Figure 3 (taking $\omega = \sqrt{\psi_i^\varepsilon}$). Recalling that $\sqrt{\psi_i^\varepsilon} \ge r$, each of these disk segments has area bounded above by $\frac{\pi\varepsilon^3}{2r}$, so we have $\alpha_{\min} \le \rho^\varepsilon(\text{RLag}_i^\varepsilon) \le 2 \times \frac{\rho_{\max}}{2\varepsilon} \times \frac{\pi\varepsilon^3}{2r} = \frac{\pi\rho_{\max}}{2r}\varepsilon^2$. As a result, $\varepsilon \ge \sqrt{\frac{2r\alpha_{\min}}{\pi\rho_{\max}}} = \sqrt{\frac{4}{\pi}}r > r$, which is a contradiction. Note that since we only assumed $\varepsilon < r$, the rougher bound mentioned in Figure 3 actually suffices to get the wanted contradiction. $\square$

**Lemma 24.** If $0 < \varepsilon < \varepsilon_0$, then for all $i \in [1, N-1]$ we have

$$|M_{i,i+1}^{\varepsilon}| \ge \beta \quad \text{or} \quad \begin{cases} |M_{i,0}^{\varepsilon}| \ge \beta, \\ |M_{i+1,0}^{\varepsilon}| \ge \beta, \end{cases}$$

where $\beta = \frac{\rho_{\min}}{4R}$.

From this lemma, we construct as follows a graph $\mathcal{G}_u^\varepsilon$ on $\{1, \dots, N, \infty\}$ with weights in $\{0, 1\}$. For any two distinct vertices $i \ne j$, we put a weight equal to one between $i$ and $j$ if and only if $|M_{i,j}^\varepsilon| \ge \beta$. All the other weights are zero. Figure 4 provides an example of a possible graph constructed in this way. Lemma 24, implies that $\mathcal{G}_u^\varepsilon$ is connected. The proof amounts to showing that if two consecutive vertices $i$ and $i+1$ are not connected in the graph, then both are connected to the auxiliary vertex $\infty$, cf Figure 5.

*Proof of Lemma 24.* For $i \in [1, N]$, define

$$\begin{cases} R_i^\varepsilon = \text{RLag}_i^\varepsilon \cap \partial B_i^{\varepsilon,+}, \\ L_i^\varepsilon = \text{RLag}_i^\varepsilon \cap \partial B_i^{\varepsilon,-}, \end{cases}$$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Example of a graph $\mathcal{G}_u^\epsilon$, where the connected components of $\mathcal{G}_u^\epsilon \setminus \{\infty\}$ are $\{1, 2, 3, 4\}$, $\{5, 6, 7\}$ and $\{8, 9\}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Illustration for the sets $R_i^\epsilon$, $E_i^\epsilon$ and $L_{i+1}^\epsilon$, respectively in blue, in red, and in green. The dotted vertical line corresponds to the boundary between the two unrestricted Laguerre cells.</figcaption></figure>

where $\partial B_i^{\epsilon,+} := \partial B_i^{\epsilon} \cap ([y_i, +\infty) \times \mathbb{R})$ and $\partial B_i^{\epsilon,-} := \partial B_i^{\epsilon} \cap ((-\infty, y_i]) \times \mathbb{R})$ are the right-hand and left-hand sides of the ball $B_i^\epsilon$. In addition, let us denote

$$ E_i^\epsilon = \begin{cases} \mathrm{RLag}_1^\epsilon \cap (\{\min \mathrm{spt} \rho\} \times \mathbb{R}) & \text{if } i=0, \\ \mathrm{RLag}_i^\epsilon \cap \mathrm{RLag}_{i+1}^\epsilon & \text{if } i \in [1, N-1], \\ \mathrm{RLag}_N^\epsilon \cap (\{\max \mathrm{spt} \rho\} \times \mathbb{R}) & \text{if } i=N. \end{cases} $$

Since $\epsilon < r$, Lemma 23 ensures that the restricted Laguerre cells all intersects the horizontal lines $\mathbb{R} \times \{\pm\epsilon\}$. As a result, for every $i \in [1, N-1]$ we have

$$ \left\{ \begin{array}{ll} \mathcal{H}^1(E_i^\epsilon) + \mathcal{H}^1(R_i^\epsilon) & \geq 2\epsilon, \\ \mathcal{H}^1(E_i^\epsilon) + \mathcal{H}^1(L_{i+1}^\epsilon) & \geq 2\epsilon, \end{array} \right. \qquad (4.6) $$

as illustrated in Figure 5. Recall that $\rho^\epsilon$ is bounded below by $\frac{\rho_{\min}}{2\epsilon}$ on its support, that the components of $\psi^\epsilon$ are all bounded above by $R$, and that the distance between the positions of two consecutive Dirac masses is at most $R$. From (4.4) and from the definition of $\beta$ we deduce that

$$ \left\{ \begin{array}{ll} |M_{i,i+1}^\epsilon| & \geq \beta \mathcal{H}^1(E_i^\epsilon)/\epsilon, \\ |M_{i,0}^\epsilon| & \geq \beta \mathcal{H}^1(R_i^\epsilon)/\epsilon, \\ |M_{i+1,0}^\epsilon| & \geq \beta \mathcal{H}^1(L_{i+1}^\epsilon)/\epsilon. \end{array} \right. $$

These inequalities combined with (4.6) allow to conclude the proof. $\square$

Lemma 24 implies that for $\epsilon$ small enough, $M^\epsilon$ is uniformly connected. More precisely, the weights of its associated graph $\mathcal{G}^\epsilon$ are bounded below by the constant $\beta$ times that of the uniformly weighted graph $\mathcal{G}_u^\epsilon$ introduced after Lemma 24. We
<figure><img src="image_7.png" /><figcaption>FIGURE 6. The restricted Laguerre cells represented on the support of $\rho^\epsilon$. The vertical segments in red correspond to the sets $E_i^\epsilon$, while the blue area is the relative complement of the union of the Laguerre cells in spt$\rho^\epsilon$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. The blue area in the left drawing is at most equal to the blue area in the right drawing.</figcaption></figure>

can thus harness a connectivity inequality by Fiedler on simple unweighted graphs in order to get a lower bound for the smallest non-zero eigenvalue of $M^\epsilon$. This, in turn, implies that the matrices $H^\epsilon$ are uniformly positive definite, thanks to Lemma 22.

**Remark 25.** Since the potential $\psi^\epsilon$ will be shown to converge to $\psi^*$ (Theorem 17), for $\epsilon$ small enough the graph $G_u^\epsilon$ is in fact independent of $\epsilon$. Indeed, two consecutive vertices $i$ and $i+1$ will be connected if and only if the corresponding optimal cells for the unregularized multi-dimensional problem intersect, and otherwise both will be connected to $\infty$. Note however that the proof we provide for Theorem 17 does not rely on this property.

A graph theory inequality by Fiedler. The following result can be found in [10, Paragraph 4.3]. Let $G$ be a graph with $n$ vertices, and let $e(G)$ denote its edge connectivity, i.e. the minimum number of edges whose deletion disconnects the graph. Then the second smallest eigenvalue $\lambda_1(G)$ of its Laplacian matrix satisfies

$$ \lambda_1(G) \geq 2e(G) (1 - \cos(\pi n^{-1})). \quad (4.7) $$
**Proposition 26.** For $0 < \varepsilon < \varepsilon_0$ we have

$$ \lambda_{\min}(H^{\varepsilon}) \geq \frac{4\beta}{N+1} \sin^2\left(\frac{\pi}{2N+2}\right) \qquad (4.8) $$

where we recall that

$$ 4\beta = \rho_{\min} \left( 1 + \operatorname{diam} \left( \operatorname{spt} \rho \cup \{y_1, \dots, y_N\} \right)^2 \right)^{-\frac{1}{2}}. $$

*Proof.* By Lemma 24 and Lemma 21 we have $\lambda_1(M^\varepsilon) \geq \beta\lambda_1(\mathcal{G}_u^\varepsilon)$. Since $\mathcal{G}_u^\varepsilon$ is connected, its edge connectivity is at least 1, so Fiedler's connectivity inequality combined with the identity $1 - \cos(2\theta) = 2\sin^2(\theta)$ yields $\lambda_1(\mathcal{G}_u^\varepsilon) \geq 4\sin^2\left(\frac{\pi}{2N+2}\right)$. Lemma 22 allows to conclude the proof. $\square$

We now express the derivative of $G_\varepsilon$ with respect to $\varepsilon$, and give an estimate for its value at the optimal potential $\psi^\varepsilon$. The main theorem will then follow from the implicit function theorem, combined with the uniform bound 4.8 and the integration of $\varepsilon \mapsto \dot{\psi}^\varepsilon$ from 0 to $\varepsilon$.

**Proposition 27.** Suppose that $0 < \varepsilon < \varepsilon_0$. Then, for every $i \in [1, N]$, the partial derivative of $G_i^\varepsilon$ with respect to $\varepsilon$ writes

$$ (\partial_{\varepsilon} G_i^{\varepsilon})(\psi) = -\frac{1}{\varepsilon} \int_{\mathrm{RLag}_i(\psi)} \rho^{\varepsilon}(x) dx + \int_{\mathrm{RLag}_i(\psi) \cap (\mathbb{R} \times \{-\varepsilon, \varepsilon\})} \rho^{\varepsilon}(x) d\mathcal{H}^1(x). \quad (4.9) $$

In particular, the function $(\varepsilon, \psi) \mapsto (\partial_\varepsilon G^\varepsilon)(\psi)$ is jointly continuous on the open set $\Omega$ of pairs $(\varepsilon, \psi) \in (0, \varepsilon_0) \times \mathbb{R}^N$ such that all components of $G^\varepsilon(\psi)$ are non-zero, as well as $1 - \sum_{i=1}^N G_i^\varepsilon(\psi)$. Moreover, for every $0 < \varepsilon < \varepsilon_0$ we have

$$ \|q^{\varepsilon}\| \leq \frac{4\rho_{\max}^2}{\alpha_{\min}} N^{\frac{1}{2}} \varepsilon. \qquad (4.10) $$

*Proof.* The formula for $(\partial_\varepsilon G_i^\varepsilon)(\psi)$ follows from Fubini's theorem, which allows us to write

$$ \begin{aligned} G_i^\varepsilon(\psi) &= \int_{\mathrm{RLag}_i(\psi) \cap (\mathbb{R} \times [-\varepsilon, \varepsilon])} d\rho^\varepsilon \\ &= \frac{1}{2\varepsilon} \int_{-\varepsilon}^{\varepsilon} \left( \int_{\mathrm{RLag}_i(\psi) \cap (\mathbb{R} \times \{w\})} \rho(x) d\mathcal{H}^1(x) \right) dw. \end{aligned} $$

The joint continuity of $\partial_\varepsilon G^\varepsilon$ in $(\varepsilon, \psi)$ then immediately follows from the dominated convergence theorem. Let us now bound $q^\varepsilon = -(\partial_\varepsilon G_i^\varepsilon)(\psi^\varepsilon)$. We take $\varepsilon < \varepsilon_0$, so that every cell contains a point of ordinate $\varepsilon$. Leveraging the definition of $\rho^\varepsilon$, we can write the second term in the right-hand side of (4.9) as $1/\varepsilon$ times the $\rho^\varepsilon$-measure of the rectangle which is the convex hull of $\mathrm{RLag}_i^\varepsilon \cap \partial(\mathbb{R} \times [-\varepsilon, \varepsilon])$. More precisely,

$$ \begin{aligned} q_i^\varepsilon &= \frac{1}{\varepsilon} [\rho^\varepsilon(\mathrm{RLag}_i^\varepsilon) - \rho^\varepsilon(\mathrm{RLag}_i^\varepsilon \cap H_i^\varepsilon)] \\ &= \frac{1}{\varepsilon} \rho^\varepsilon(\mathrm{RLag}_i^\varepsilon \setminus H_i^\varepsilon), \end{aligned} $$

where $H_i^\varepsilon$ is the vertical strip of width $2\sqrt{\psi_i^\varepsilon - \varepsilon^2}$ centered at abscissa $y_i$. The set $\mathrm{RLag}_i^\varepsilon \setminus H_i^\varepsilon$ is contained in the union of two disjoint disk segments of base $2\varepsilon$ and
height $\sqrt{\psi_i^\varepsilon} - \sqrt{\psi_i^\varepsilon - \varepsilon^2}$, as illustrated in Figure 3 when choosing $\omega = \sqrt{\psi_i^\varepsilon}$. Since the height of the disk segments is bounded above by $\varepsilon^2/\sqrt{\psi_i^\varepsilon}$, their respective area is at most $2\varepsilon^3/\sqrt{\psi_i^\varepsilon}$. The density $\rho^\varepsilon$ is moreover bounded by $\rho_{\max}/2\varepsilon$, and by Lemma 23 we have $\sqrt{\psi_i^\varepsilon} \ge r = \alpha_{\min}/(2\rho_{\max})$, so in the end

$$|q_i^\varepsilon| \le 2 \times \frac{\rho_{\max}}{2\varepsilon} \times \frac{2\varepsilon^3}{\sqrt{\psi_i^\varepsilon}} = \frac{4\rho_{\max}^2}{\alpha_{\min}}\varepsilon,$$

and the bound on $\|q^\varepsilon\|$ follows. $\square$

*Proof of Theorem 17.* We apply the implicit function theorem to the function $G: \Omega \to \mathbb{R}^N$ defined by $G(\varepsilon, \psi) = G^\varepsilon(\psi)$, where $\Omega$ is the open set defined in Proposition 27. Thanks to this proposition and Theorem 8, $G$ has partial first derivatives in both $\varepsilon$ and $\psi$, and these first derivatives are jointly continuous on $\Omega$, so $G$ is $C^1$ on this set. Note that continuity of the first derivative in $\psi$ is a direct consequence of the dominated convergence and of the fact that the derivative in question may be written

$$\partial_{\psi_{i+1}} G_i(\varepsilon, \psi) = (4(y_{i+1} - y_i)\varepsilon)^{-1} \int_{\mathbb{R}} \rho(z_i(\psi)) \mathbb{1}_{[-w_i(\psi), w_i(\psi)]}(x^2) \mathbb{1}_{[-\varepsilon, \varepsilon]}(x^2) dx^2,$$

where $w_i(\psi) = \sqrt{\psi_i - (z_i(\psi) - y_i)^2}$ if the radicand is non-negative, $w_i(\psi) = 0$ otherwise. Since the Jacobian matrix $DG^\varepsilon(\psi^\varepsilon)$ is invertible for every $\varepsilon \in (0, \varepsilon_0)$, the implicit function theorem implies that $\varepsilon \mapsto \psi^\varepsilon$ is of class $C^1$ on this interval, with derivative given by (4.1).

We now prove the convergence of $\psi^\varepsilon$ to $\psi^*$ as $\varepsilon \to 0$. Let $(\varepsilon_n)$ be a sequence of strictly positive reals converging to zero. By Lemma 23, the $\psi^{\varepsilon_n}$ are uniformly bounded so there exists a converging subsequence. Since the functionals $-K^\varepsilon \Gamma$ converge to $-K$, the limit of the subsequence is necessarily $\psi^*$, the unique maximizer of $K$. This holds for any convergent subsequence of $(\psi^{\varepsilon_n})$, so the whole sequence converges to $\psi^*$.

Now that the convergence of $\psi^\varepsilon$ to $\psi^*$ is established, let us bound $\|\psi^\varepsilon - \psi^*\|$. For $\varepsilon \in (0, \varepsilon_0)$, equation (4.1) and the inequalities (4.8) and (4.10) imply that

$$\|\dot{\psi}^\varepsilon\| \le \frac{1}{\lambda_{\min}(H^\varepsilon)} \|q^\varepsilon\| \\ \le \frac{N+1}{4\beta \sin^2\left(\frac{\pi}{2N+2}\right)} \times \frac{2\rho_{\max}}{r} N^{\frac{1}{2}}\varepsilon.$$

As a result,

$$\|\psi^\varepsilon - \psi^*\| \le \int_0^\varepsilon \|\dot{\psi}^\tau\| d\tau \le 2C \int_0^\varepsilon \tau d\tau = C\varepsilon^2$$

where, using the definitions of $\beta$ and $r$,

$$C = \frac{(N+1)\rho_{\max}N^{\frac{1}{2}}}{8\beta \sin^2\left(\frac{\pi}{2N+2}\right) r} = \frac{\rho_{\max}R(N+1)N^{\frac{1}{2}}}{2\sin^2\left(\frac{\pi}{2N+2}\right) \rho_{\min}r} = \frac{\rho_{\max}^2 R(N+1)N^{\frac{1}{2}}}{\alpha_{\min}\rho_{\min}\sin^2\left(\frac{\pi}{2N+2}\right)},$$

and we recall that $R = (1 + \operatorname{spt}\rho \cup \{y_1, \dots, y_N\})^2^{\frac{1}{2}}$. $\square$
## 5. NUMERICS

The aim of this section is to illustrate the algorithm's behavior, particularly in terms of accuracy and execution speed. The solver is based on Newton's algorithm applied to the regularized dual problem. To avoid leaving the domain on which the functional is $C^2$, backtracking has been implemented: if the potentials proposed at a given step lead to empty cells, we go back to the previous step and we apply a relaxation coefficient twice smaller. Each step of the algorithm writes

$$ \psi^{(k+1)} = \psi^{(k)} + \eta^{(k)} d^{(k)}, \qquad (5.1) $$

where $d^{(k)} = -DG^{\epsilon}(\psi^{(k)})[G^{\epsilon}(\psi^{(k)}) - \alpha]$ is the direction of descent and $\eta^{(k)}$ is the step size given by the backtracking line search.

### 5.1. Implementation details.
The implementation naturally exploits the properties of the one-dimensional case. First, the restricted Laguerre cells in 1D are straightforward to compute once the positions of the Dirac masses are sorted, see equation (2.18) in Appendix B. Second, the matrix from Newton's system is both symmetrical and tridiagonal. It is therefore possible to store the complete system using three vectors (two for the matrix, and one for the second member). However, in practice, to increase efficiency the solution of the system is calculated at the same time as the cells are constructed. This allows to store only one intermediate vector. As in classical semi-discrete optimal transport, the pure Newton's algorithm (i.e. with constant step size $\eta^{(k)} \equiv 1$) may not converge. We therefore use a backtracking strategy, similar to that of [17], to ensure that the iterates do not reach an invalid state (empty cells, negative potentials, etc.). The library developed for this publication is available at https://github.com/sdot-team/usdot.

### 5.2. Regularization error.
The impact of regularization on the computed $\psi$ is illustrated in Figure 8. Although this example corresponds to a specific numerical setting, similar behavior was consistently observed across all our experiments: the error scales with the square of $\epsilon$ and can easily become negligibly small. In particular, these simulations illustrate the tight convergence rate provided by Theorem 17.

### 5.3. Shape registration example.
Bonneel et al. [3] introduced a shape registration algorithm based on a 1D discrete-discrete optimal transport algorithm, known as the Fast Iterative Sliced Transport (FIST) algorithm. This method leverages 1D point correspondences to generate unidimensional displacement proposals, which are then used as combinations to propose movements in the target space. Although these proposed displacements may not perfectly match those obtained through direct optimal transport in the target space, the approximation can still effectively capture global motions, such as rigid body transformations.

This algorithm has been adapted for the semi-discrete setting, the main difference being the representation of the target shape. In our case, the projections can be defined as almost any kind of functions, for instance sum of polynomials, Gaussians, etc. It is still possible to start from a set of points, for instance by binning the projections in histograms, but the semi-discrete setting offers the possibility of starting from more generic and more precise representations. For instance, if the target shape is defined by a 3D triangular mesh, one simply has to sum the projection of
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Impact of the regularization parameter $\varepsilon$ on the potentials $\psi$. For this example, the density function was a Gaussian centered around 0, with $\sigma = 1$, the discrete measure is supported on 15 Dirac masses which were positioned randomly in [-1, +1] with a uniform distribution, and the ratio of the total mass of the Dirac masses to that of the Gaussian was set at 1/2. As can be seen, the error on the potentials is globally proportional to the square of the $\varepsilon$ parameter, in accordance with the theorem 17.</figcaption></figure>

the triangles into a piecewise affine function. One may consider that the triangles have some thickness to avoid Dirac masses in borderline cases.

For the semi-discrete setting, the unidimensional displacement proposals are defined by the displacements between the Dirac masses and the barycenter of the corresponding cells. Figure 9 illustrates the two numerical experiments used to benchmark the semi-discrete approach for the FIST algorithm.

In these examples, computation time is strongly influenced by the distance to the solution. Nonetheless, in general, we observe that semi-discrete problems are approximately 4 to 5 times slower than fully discrete problems solved using Bonneel's algorithm ([3]). This additional cost can be seen as the trade-off for increased generality and potentially higher accuracy. However, there is clearly room for further optimization, particularly to accelerate convergence. For instance, in the rigid body registration examples, the Newton method requires an average of 63 iterations — a very high number. While strategies to reduce this iteration count have been explored, they have not yet been implemented within the regularization framework proposed in this work.

**5.4. Synthetic example.** In this example, we explore the gains in accuracy that come from being able to use generic densities, without having to discretize them into sums of Dirac masses. The density $\rho$ for this example is equal to $(1 - |x|)_+$. We start with 100 randomly placed points in the interval [-1, 1]. They are associated with uniform weights, such that the sum is equal to 3/4 (i.e. 3/4 of the mass of
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Illustration of an iteration of the FIST algorithm for two semi-discrete settings. For the left image, the unidimensional densities are computed using the projected mass of the triangles, which gives a piecewise affine function. For the right image, unidimensional densities are histogram of the projected points. The black items (lines and dots) represent the target shape. The red dots represent the position of the Dirac masses (the source points). The blue lines represent the displacement from the Dirac positions to the estimated 3D barycenters, that come from the combinations of the 1D barycenters with their respective directions.</figcaption></figure>

$\rho$). The barycenters are then computed, with the semi-discrete method, and with the discrete-discrete method with a growing number of points for the discretization of $\rho$. For the discrete-discrete setting, the spot library has been used (https://github.com/nbonneel/spot). It is highly optimized in terms of execution time. Nevertheless, the library requires that there be a correspondence for each Dirac. We therefore repeated the list of initial Dirac masses in order to obtain the right mass ratio for the partial transport problem.

Figure 10 shows the accuracy losses due to the discretization of $\rho$ used to reduce to a discrete-discrete problem, while Figure 11 shows the relative timings to compute the barycenters using a the discrete-discrete and semi-discrete settings, respectively. As expected, in this case, the discrete-discrete leads to discretization errors. Of course, increasing the accuracy leads to increased execution times, and can quickly exceed to time needed to solve the exact problem with a semi-discrete setting.

## ACKNOWLEDGMENTS

The authors would like to warmly thank Quentin Mérigot and Luca Nenna for their guidance and support in the elaboration of this work. This work has been partially supported by the UDOPIA doctoral program, as well as the Agence nationale de la recherche, through the ANR project GOTA (ANR-23-CE46-0001) and the PEPR PDE-AI project (ANR-23-PEIA-0004).
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Accuracy losses for a discrete-discrete setting vs a semi-discrete setting for $\rho = (1 - |x|)_+$ with 100 Dirac masses. The accuracy is measured as the Euclidean distance between the exact barycenters (computed using semi-discrete optimal transport), and the ones computed using the points taken from the discretization of $\rho$ with a discrete-discrete setting.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Execution time to solve the partial optimal transport problem with a discrete-discrete (DD) setting and with a semi-discrete (SD) setting. The Dirac masses for the count in the abscissa are the one that are used to discretize $\rho$, to be able to use the discrete-discrete setting</figcaption></figure>
REFERENCES

[1] Yikun Bai, Bernhard Schmitzer, Matthew Thorpe, and Soheil Kolouri. Sliced optimal partial transport. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, pages 13681–13690, 2023.

[2] Nicolas Bonneel and David Coeurjolly. Spot: sliced partial optimal transport. *ACM Transactions on Graphics (TOG)*, 38(4):1–13, 2019.

[3] Nicolas Bonneel, Julien Rabin, Gabriel Peyré, and Hanspeter Pfister. Sliced and radon wasserstein barycenters of measures. *Journal of Mathematical Imaging and Vision*, 51:22–45, 2015.

[4] Yann Brenier. Polar factorization and monotone rearrangement of vector-valued functions. *Communications on pure and applied mathematics*, 44(4):375–417, 1991.

[5] Luis A Caffarelli and Robert J McCann. Free boundaries in optimal transport and monge-ampere obstacle problems. *Annals of mathematics*, pages 673–730, 2010.

[6] Shibing Chen and Emanuel Indrei. On the regularity of the free boundary in the optimal partial transport problem for general cost functions. *Journal of Differential Equations*, 258(7):2618–2632, 2015.

[7] Lenaic Chizat. *Unbalanced optimal transport: Models, numerical methods, applications*. PhD thesis, Université Paris sciences et lettres, 2017.

[8] Gonzalo Dávila and Young-Heon Kim. Dynamics of optimal partial transport. *Calculus of Variations and Partial Differential Equations*, 55(5):116, 2016.

[9] Hamza Ennaji, Quentin Mérigot, Luca Nenna, and Brendan Pass. Robust risk management via multi-marginal optimal transport. *Journal of Optimization Theory and Applications*, pages 1–28, 2024.

[10] Miroslav Fiedler. Algebraic connectivity of graphs. *Czechoslovak mathematical journal*, 23(2):298–305, 1973.

[11] Alessio Figalli. The optimal partial transport problem. *Archive for rational mechanics and analysis*, 195(2):533–560, 2010.

[12] Gerald B Folland. *Real analysis: modern techniques and their applications*. John Wiley & Sons, 1999.

[13] Roger A Horn and Charles R Johnson. *Matrix analysis*. Cambridge university press, 2012.

[14] Noureddine Igbida et al. Optimal partial mass transportation and obstacle monge-kantorovich equation. *Journal of Differential Equations*, 264(10):6380–6417, 2018.

[15] Noureddine Igbida and Van Thanh Nguyen. Augmented lagrangian method for optimal partial transportation. *IMA Journal of Numerical Analysis*, 38(1):156–183, 2018.

[16] Emanuel Indrei. Free boundary regularity in the optimal partial transport problem. *Journal of Functional Analysis*, 264(11):2497–2528, 2013.

[17] Jun Kitagawa, Quentin Mérigot, and Boris Thibert. Convergence of a newton algorithm for semi-discrete optimal transport. *Journal of the European Mathematical Society*, 21(9):2603–2651, 2019.

[18] Hugo Leclerc, Quentin Mérigot, Filippo Santambrogio, and Federico Stra. Lagrangian discretization of crowd motion and linear diffusion. *SIAM Journal on Numerical Analysis*, 58(4):2093–2118, 2020.

[19] Bertrand Maury, Aude Roudneff-Chupin, Filippo Santambrogio, and Juliette Venel. Handling congestion in crowd motion modeling. *arXiv preprint arXiv:1101.4102*, 2011.

[20] Julien Rabin, Gabriel Peyré, Julie Delon, and Marc Bernot. Wasserstein barycenter and its application to texture mixing. In *Scale Space and Variational Methods in Computer Vision: Third International Conference, SSVM 2011, Ein-Gedi, Israel, May 29–June 2, 2011, Revised Selected Papers 3*, pages 435–446. Springer, 2012.

[21] Daniel Spielman. Spectral and algebraic graph theory. *Yale lecture notes, draft of December*, 4:47, 2019.

[22] Cédric Villani et al. *Optimal transport: old and new*, volume 338. Springer, 2008.
## APPENDIX A. PROOF OF THEOREM 8

Thanks to Proposition 6, it is enough to prove that the functions $G_i : \psi \mapsto \int_{RLag_i(\psi)} \rho(x)dx$ are $C^1$ on $\mathcal{D}$. Notice that $\mathcal{D}$ is both open and convex, since the inequalities $\|x - y_i\|^2 - \psi_i \le \|x - y_j\|^2 - \psi_j$ and $\|x - y_i\|^2 \le \psi_i$ are respectively linear and convex in $(x, \psi)$. Fix $\psi^0 \in \mathcal{D}$ and let $\psi^t = \psi^0 + tv$, where $v \in \mathbb{R}^N$ is a non-zero vector. For the sake of readability, we denote $V_j^t = RLag_j(\psi^t)$, $V_\infty^t = RLag_\infty(\psi^t)$, and $L_j^t = Lag_j(\psi^t)$.

**Polar coordinates.** Let us express the integral in polar coordinates centered at the Dirac position $y_i$. Using [12, Theorem 2.49], we have

$$ 
\begin{aligned}
G_i(\psi^t) &= \int_{\mathbb{R}^d} \mathbb{1}_{V_i^t}(x) \rho(x) dx = \int_{\mathbb{S}^{d-1}} \left( \int_0^\infty \mathbb{1}_{V_i^t}(s\theta) \rho(y_i + s\theta) s^{d-1} ds \right) d\mathcal{H}^{d-1}(\theta) \\
&= \int_{\mathbb{S}^{d-1}} \left( \int_0^{\omega_i^t(\theta)} \rho(y_i + s\theta) |s|^{d-1} ds \right) d\mathcal{H}^{d-1}(\theta), \quad (A.1)
\end{aligned}
$$

where $\omega_i^t(\theta)$ is the supremum of the set of real numbers $s \in \mathbb{R}$ for which $y_i + s\theta \in V_i^t$, or 0 if no such $s$ exists. The last equality is obtained by making the change of variable $\theta \leftrightarrow -\theta$ in half of the integral, and summing back the two halves. Note that equation (A.1) holds even when $y_i$ is not in the interior of $V_i^t$, because the $\omega_i^t(\theta)$ can be negative.

**Cover of the sphere.** We now define a cover of the sphere $\mathbb{S}^{d-1}$ by $N + 1$ sets with $\mathcal{H}^{d-1}$-negligible pairwise intersections. These sets correspond to the different facets of the boundary of the cell (some of which may be empty), except for the last set, which corresponds to the lines passing through $y_i$ without intersecting the cell. Namely, for each $j \in \{1, \dots, N, \infty\}$ distinct from $i$, we let $\Theta_{i,j}^t$ be the set of unit vectors $\theta \in \mathbb{S}^{d-1}$ such that $y_i + \omega_i^t(\theta) \in V_j^t$. As for $\Theta_{i,i}^t$, we define it as the set of unit vectors which do not belong to any other $\Theta_{i,j}^t$. It is immediate to check that

$$ \forall \theta \in \Theta_{i,j}^t, \quad \omega_i^t(\theta) = \begin{cases} \frac{\|y_j - y_i\|^2 - (\psi_j^t - \psi_i^t)}{2\langle y_j - y_i, \theta \rangle} & \text{if } j \notin \{i, \infty\} \\ \sqrt{\psi_i^t} & \text{if } j = \infty, \\ 0 & \text{if } j = i. \end{cases} \quad (A.2) $$

The first expression is derived by developing the squares in the equality $\|(y_i + \omega\theta) - y_i\|^2 - \psi_i^t = \|(y_j + \omega\theta) - y_j\|^2 - \psi_j^t$, which defines the hyperplane of points $x = y_i + \omega\theta$ for which traveling to $y_i$ costs exactly as much as traveling to $y_j$. The same goes for the second expression, except that we replace the second member of the equality by zero.

**Differentiating under the integral sign.** We now make the temporary additional assumption that $y_i$ is not contained in any of the $N$ iso-hypersurfaces defined by $c(x, y_i) - \psi^0 = c(x, y_j) - \psi_j^t$, $j \in \{1, \dots, N, \infty\} \setminus \{i\}$. Thanks to this assumption, for $t$ in a small enough neighborhood of zero, each line passing through $y_i$ intersects the boundary of the cell in at most two points. In particular, since the facets of the cell vary continuously in $t$, every $\theta$ in the relative interior of some $\Theta_{i,j}^0$ is also in $\Theta_{i,j}^t$.
for each $t$ in a small enough neighborhood of 0 (which depends on $\theta$). As a result, $t \mapsto \omega_i^t(\theta)$ is differentiable around 0 for every such $\theta$, and hence for $\mathcal{H}^{d-1}$-almost every $\theta \in \mathbb{S}^{d-1}$, and equation (A.2) can be used to find the derivative. We can therefore differentiate under the integral sign in (A.1) and find

$$ 
\begin{aligned}
\partial_{t=0}[G_i(\psi^t)] &= \int_{\Theta_i^0} [\partial_t|_{t=0}\omega_i^0(\theta)] \rho(y_i + \omega_i^0(\theta)\theta) |\omega_i^0(\theta)|^{d-1} d\mathcal{H}^{d-1}(\theta) \\
&= \sum_{j \in \{1, \dots, N\} \setminus \{i\}} \frac{-(v_j - v_i)}{2\langle y_j - y_i, \theta \rangle} \int_{\Theta_{i,j}^0} \rho(y_i + \omega_{i,j}^0(\theta)\theta) |\omega_{i,j}^0(\theta)|^{d-1} d\mathcal{H}^{d-1}(\theta) \\
&\quad + \frac{v_i}{2\sqrt{\psi_i^0}} \int_{\Theta_{i,\infty}^0} \rho(y_i + \omega_{i,\infty}^0(\theta)\theta) |\omega_{i,\infty}^0(\theta)|^{d-1} d\mathcal{H}^{d-1}(\theta).
\end{aligned}
$$

Since all terms are linear in $v$ and since $v$ was an arbitrary vector, we conclude that $G_i$ is differentiable at $\psi^0$. The changes of variables $\theta \mapsto y_i + \omega_{i,j}^0(\theta)\theta$ allow to recover the expressions given in (2.15) and (2.16). Indeed, thanks to the additional assumption, these maps are diffeomorphisms from their respective facet $V_i^t \cap V_j^t$, $j \neq i$ of the cell onto their corresponding set $\Theta_{i,j}^t$ of unit vectors. The continuity in $\psi \in \mathcal{D}$ of the facet integrals appearing in (2.15) and (2.16) can be established as in the proof of [17, Proposition B.1] by Mérigot et al.

Since nonzero mass of every Laguerre cell implies strict positivity of all the components of the potential, the $\psi \in \mathcal{D}$ which do not satisfy the additional assumption are contained in the $N-1$ hyperplanes of $\mathbb{R}^N$ of equations $\psi_j - \psi_i = \|y_j - y_i\|^2$. By classical arguments, $G_i$ is thus $\mathcal{C}^1$ on the whole domain $\mathcal{D}$, with first partial derivatives given by (2.15) and (2.16), up to a minus sign. In fact, the singular $\psi$ are an artifact due to the polar coordinates being centered at $y_i$. Centering these coordinates at some arbitrary point in the interior of the cell (e.g. the $\rho$-barycenter of the cell), one can directly prove differentiability of $G_i$ on all of $\mathcal{D}$. At the cost, of course, of a less concise expression for $\omega_i^t(\theta)$ when $\theta \in \Theta_{i,\infty}^t$, because projecting a sphere on a ball is less simple when they are not concentric.

**Strict concavity.** Since $G$ is $\mathcal{C}^1$ on $\mathcal{D}$, the functional $\mathcal{K}$ is $\mathcal{C}^2$ on the same domain. Thanks to (2.15) and (2.16), for $\psi \in \mathcal{D}$ the Hessian matrix $D^2\mathcal{K}(\psi)$ is irreducible, weakly diagonally dominant, and has at least one strictly dominant diagonal coefficient — take any cell $RLag_i$ which has non-degenerate intersection with $spt\rho \setminus \cup_{j \neq i} RLag_j$. The invertibility of $D^2\mathcal{K}$ then follows from Taussky's theorem [13, Theorem II], and in particular, the optimal potential is unique.

## APPENDIX B. PROOF OF LEMMA 10

As explained in section 2.5, for $\psi \in \mathcal{D}$, the (unrestricted) Laguerre cells write $Lag_i = [z_{i-1}, z_i]$, where

$$ z_i(\psi) = \begin{cases} -\infty & \text{if } i=0, \\ \frac{(y_{i+1}^2 - \psi_{i+1}) - (y_i^2 - \psi_i)}{2(y_{i+1} - y_i)} & \text{if } 0 < i < N, \\ +\infty & \text{if } i = N, \end{cases} $$
and we can therefore express the restricted Laguerre cells as $RLag_i = [a_i, b_i]$, where

$$ \begin{cases} a_i(\psi) = \max\{z_{i-1}(\psi), y_i - \sqrt{\psi_i}\}, \\ b_i(\psi) = \min\{z_i(\psi), y_i + \sqrt{\psi_i}\}. \end{cases} $$

**Lemma 28.** Fix $\psi \in \mathcal{D}$ and let $v \in \mathbb{R}^N$ be non-zero. Then

$$ \partial_v^+ a_i(\psi) = \begin{cases} \frac{-(v_i - v_{i-1})}{2(y_i - y_{i-1})} & \text{if } z_{i-1}(\psi) > y_i - \sqrt{\psi_i}, \\ \frac{-v_i}{2\sqrt{\psi_i}} & \text{if } z_{i-1}(\psi) < y_i - \sqrt{\psi_i}, \\ \max\left\{\frac{-(v_i - v_{i-1})}{2(y_i - y_{i-1})}, \frac{-v_i}{2\sqrt{\psi_i}}\right\} & \text{if } z_{i-1}(\psi) = y_i - \sqrt{\psi_i}, \end{cases} \quad (B.1) $$

and

$$ \partial_v^+ b_i(\psi) = \begin{cases} \frac{v_i}{2\sqrt{\psi_i}} & \text{if } z_i(\psi) > y_i + \sqrt{\psi_i}, \\ \frac{-(v_{i+1} - v_i)}{2(y_{i+1} - y_i)} & \text{if } z_i(\psi) < y_i + \sqrt{\psi_i}, \\ \min\left\{\frac{-(v_{i+1} - v_i)}{2(y_{i+1} - y_i)}, \frac{v_i}{2\sqrt{\psi_i}}\right\} & \text{if } z_i(\psi) = y_i + \sqrt{\psi_i}. \end{cases} \quad (B.2) $$

Moreover, the unilateral directional derivative of the mass of the $i^{th}$ restricted Laguerre cell is

$$ \partial_v^+ G_i = \rho(b_i) \partial_v^+ b_i - \rho(a_i) \partial_v^+ a_i, \quad (B.3) $$

so we can write $\partial_v^+ G(\psi) = H(\psi, v)v$ where $H(\psi, v) \in \mathbb{R}^{N \times N}$ is a tridiagonal matrix whose rows are all weakly diagonally dominant.

There may be several possibilities for the coefficients of $H(\psi, v)$, so we fix coefficients to be those which are directly readable from (B.1) and (B.2). A case where there might be an ambiguity is for instance the case where $z_{i-1}(\psi) = y_i - \sqrt{\psi_i}$ and where the two quantities in the maximum in (B.1) are equal, in which case we choose to read the coefficients from the quantity on the right in the maximum. We use the same convention for the other case of ambiguity, namely when $z_i(\psi) = y_i + \sqrt{\psi_i}$ and both quantities in the minimum in (B.2) are equal, in which case we read the coefficients from the second argument of the minimum.

**Proof.** The expressions of $\partial_v^+ a_i(\psi)$ and $\partial_v^+ b_i(\psi)$ come from the fact that

$$ \partial_v^+ z_i(\psi) = \frac{-(v_{i+1} - v_i)}{2(y_{i+1} - y_i)} \quad \text{and} \quad \partial_v^+ \sqrt{\psi_i} = \frac{v_i}{2\sqrt{\psi_i}}. $$

As for the expression of $\partial_v^+ G_i$, it follows directly from the fact that $G_i = F_\rho(b_i) - F_\rho(a_i)$, where $F_\rho$ is the cumulative probability distribution of $\rho$. At last, from (B.1), (B.2), and (B.3), we deduce that $H(\psi, v)$ is tridiagonal, with weakly diagonally dominant rows. $\square$

**Lemma 29.** Fix $\psi \in \mathcal{D}$ and let $v \in \mathbb{R}^N$ be non-zero. Then the matrix $H(\psi, v)$ defined in Lemma 28 is symmetric. As a result, it writes

$$ H(\psi, v) = \begin{pmatrix} H_1 & & \\ & \ddots & \\ & & H_r \end{pmatrix} \quad (B.4) $$
where each $H_j$ is a symmetric, tridiagonal, irreducible, diagonally dominant matrix. Of course, the right-hand side also depends on $\psi$ and $v$, but we do not mention it in the notation for the sake of readability.

*Proof.* Since $H(\psi, v)$ is tridiagonal, it suffices to prove that its subdiagonal and superdiagonal are equal. We thus fix an index $1 \le i \le N-1$. Since the restricted Laguerre cells of index $i$ and $i+1$ are non-empty, we must have $y_i - \sqrt{\psi_i} < z_i(\psi)$ and $z_i(\psi) < y_{i+1} + \sqrt{\psi_{i+1}}$. Denoting $\delta = y_{i+1} - y_i$, we have

$$ \begin{cases} z_i(\psi) - (y_i - \sqrt{\psi_i}) = \frac{1}{2\delta} [\delta^2 + 2\sqrt{\psi_i}\delta - (\psi_{i+1} - \psi_i)], \\ (y_{i+1} + \sqrt{\psi_{i+1}}) - z_i(\psi) = \frac{1}{2\delta} [\delta^2 + 2\sqrt{\psi_{i+1}}\delta + (\psi_{i+1} - \psi_i)], \end{cases} \quad (B.5) $$

and the roots of the two left-hand side polynomials are respectively $-\sqrt{\psi_i}$ and $\pm\sqrt{\psi_i} - \sqrt{\psi_{i+1}}$, so since the right-hand sides are positive we deduce $|\sqrt{\psi_{i+1}} - \sqrt{\psi_i}| < \delta$.

Similarly,

$$ \begin{cases} z_i(\psi) - (y_i + \sqrt{\psi_i}) = \frac{1}{2\delta} [\delta^2 - 2\sqrt{\psi_i}\delta - (\psi_{i+1} - \psi_i)], \\ (y_{i+1} - \sqrt{\psi_{i+1}}) - z_i(\psi) = \frac{1}{2\delta} [\delta^2 - 2\sqrt{\psi_{i+1}}\delta + (\psi_{i+1} - \psi_i)], \end{cases} \quad (B.6) $$

and the roots of the left hand-sides are respectively $\sqrt{\psi_i} \pm \sqrt{\psi_{i+1}}$ and $\pm\sqrt{\psi_i} + \sqrt{\psi_{i+1}}$, so unless $\delta = \sqrt{\psi_i} + \sqrt{\psi_{i+1}}$ the two polynomials are non-zero and have the same sign. In other words, we have either

$$ y_i + \sqrt{\psi_i} < z_i(\psi) < y_{i+1} - \sqrt{\psi_{i+1}} \quad (B.7) $$

or

$$ y_{i+1} - \sqrt{\psi_{i+1}} < z_i(\psi) < y_i + \sqrt{\psi_i}. \quad (B.8) $$

**Case 1** If (B.7) holds, then for any $\psi'$ close enough to $\psi$, $b_i$ does not depend on $\psi'_{i+1}$ and $a_{i+1}$ does not depend on $\psi'_i$, so $H_{i+1,i}(\psi, v) = H_{i,i+1}(\psi, v) = 0$.

**Case 2** If (B.8) holds, then on a neighborhood of $\psi$ we have $b_i = a_{i+1} = z_i$, and so a straightforward computation yields $H_{i+1,i}(\psi, v) = H_{i,i+1}(\psi, v) = -\rho(z_i(\psi)) \frac{1}{2(y_{i+1}-y_i)}$, thanks to the expression of $z_i$.

**Case 3** Let us now handle the degenerate case where $\delta = \sqrt{\psi_i} + \sqrt{\psi_{i+1}}$. Then $z_i(\psi) = y_i + \sqrt{\psi_i} = y_{i+1} - \sqrt{\psi_{i+1}}$, which implies that $b_i(\psi) = a_{i+1}(\psi)$, and thanks to (2.18) we find

$$ \partial_v^+ b_i(\psi) = \min \left\{ \frac{v_i - v_{i+1}}{2(\sqrt{\psi_i} + \sqrt{\psi_{i+1}})}, \frac{v_i}{2\sqrt{\psi_i}} \right\} $$

and

$$ \partial_v^+ a_{i+1}(\psi) = \max \left\{ \frac{v_i - v_{i+1}}{2(\sqrt{\psi_i} + \sqrt{\psi_{i+1}})}, \frac{-v_{i+1}}{2\sqrt{\psi_{i+1}}} \right\}. $$
To conclude, we notice that

$$ \frac{v_i - v_{i+1}}{2(\sqrt{\psi_i} + \sqrt{\psi_{i+1}})} < \frac{v_i}{2\sqrt{\psi_i}} \iff \begin{cases} 0 < \frac{v_i}{\sqrt{\psi_i}} + \frac{v_{i+1}}{\sqrt{\psi_{i+1}}} \\ \frac{v_i - v_{i+1}}{2(\sqrt{\psi_i} + \sqrt{\psi_{i+1}})} > \frac{-v_{i+1}}{2\sqrt{\psi_{i+1}}} \end{cases} $$

and this series of equivalence still holds if all strict inequality signs are reversed. As a consequence, the coefficient of $v_{i+1}$ in $\partial_v^+ b_i(\psi)$ is the opposite of that of $v_i$ in $\partial_v^+ a_{i+1}(\psi)$, and so once again $H_{i+1,i}(\psi, v) = H_{i,i+1}(\psi, v)$. □

**Remark 30.** In the proof of Lemma 29, case **Case 1** corresponds to the situation where $RLag_i$ and $RLag_{i+1}$ touch. On the contrary, case **Case 2** is the situation where they touch and the inequality $z_i < y_i + \sqrt{\psi_i}$ is strict. At last, case **Case 3** is the degenerate situation where the two cells touch but $z_i = y_i + \sqrt{\psi_i}$, so that in any neighborhood of $\psi$ there is a $\psi'$ such that they do not touch and another such that they do.

**Lemma 31.** Each $H_j$ in the decomposition (B.4) has a strictly dominant diagonal coefficient.

*Proof.* We recall that

$$ \partial_v^+ G_i = \rho(b_i) \left[ \frac{-(v_{i+1} - v_i)}{2(y_{i+1} - y_i)} \middle| \frac{v_i}{2\sqrt{\psi_i}} \right] + \rho(a_i) \left[ \frac{v_i - v_{i-1}}{2(y_i - y_{i-1})} \middle| \frac{v_i}{2\sqrt{\psi_i}} \right] $$

where the notation $[c|d]$ stands for either expression $c$ or expression $d$. As a result, line $i$ has a strictly dominant diagonal coefficient if and only if one of $\rho(a_i)$ or $\rho(b_i)$ is non-zero and the square bracket coefficient next to it takes the value on the right.

Let us first consider a line of $H(\psi, v)$ with index $1 < i < N$. Then $\rho(a_i) > 0$ and $\rho(b_i) > 0$, since the density is strictly positive in the interior of $spt\rho$, which is convex. As a result, line $i$ is strictly diagonally dominant if and only if its $(i-1)^{th}$ or $(i+1)^{th}$ coefficient (or both) is zero. This ensures that if there are two blocks $H_j$ or more, each of these blocks has a strictly diagonally dominant line — its first or its last.

Let us now tackle the case where the decomposition consists of a single block, that is, $H(\psi, v)$ is already irreducible. Then the restricted Laguerre cells are all packed together: $b_i = a_{i+1}$ for every $1 \le i < N$. Since the complement of their union — the additional cell, associated to the auxiliary point introduced in Remark 7 — has non-zero $\rho$-measure, we must have $\min(spt\rho) < a_1$ or $b_N < \max(spt\rho)$. As a result, the density is non-zero at one of $a_1$ or $b_N$, so at least one of the first and last lines is strictly diagonally dominant. □

ADRIEN CANCES (adrien.cances@universite-paris-saclay.fr), UNIVERSITÉ PARIS-SACLAY, LABORATOIRE DE MATHÉMATIQUES D'ORSAY, PARMA, 91405, ORSAY, FRANCE

HUGO LECLERC (hugo.leclerc@universite-paris-saclay.fr), UNIVERSITÉ PARIS-SACLAY, LABORATOIRE DE MATHÉMATIQUES D'ORSAY, PARMA, 91405, ORSAY, FRANCE