Random Assignment Problems on 2d Manifolds

The Euclidean assignment problem is a transportation problem between a set of N “red” points and a set of N “blue” points. Both sets are supposed to be on a given n-dimensional Riemannian manifold \(\Omega \). A transportation map is a bijective map , that is, a pairing among the red and blue points. A transportation cost \(w_{xy}\) is given for each pair , and the cost of the map T is the expression

(1)

We will also define as the minimum cost among the possible transportation maps. Given a probability distribution for and a probability distribution for , we use \(E_{\Omega }(N)\) as a shortcut for . For definiteness, we will assume in this paper that the points of and are uniformly distributed.

Within this geometric framework, it is natural to choose for \(w_{xy}\) a function of the distance d(x, y) between and , that is \(w_{xy} = f(d(x,y))\). We expect that, if the function f has some natural (monotonicity, smoothness,...) properties, the large-N behaviour of \(E_{\Omega }(N)\) (with the volume of \(\Omega \) kept constant) is dominated by the behaviour of f near zero. In turn, this suggests to concentrate only on power-law functions, \(f(d)=d^p\) for some \(p>0\), as any other detail of the function is either trivially rescaled, or washed out in the limit. The two cases most studied in the literature are \(p=1\) and \(p=2\), where a number of (different) useful extra features emerge. This paper makes no exception, and in fact we will only consider here the case \(p=2\), that is, we will set once and for all \(w_{xy} = d^2(x,y)\).

It is a longstanding question to understand the asymptotic behaviour, for large N, of \(E_{\Omega }(N)\), and, when \(n\ge 2\), the results are very partial for any manifold \(\Omega \), including the conceptually simplest ones (like the unit hypercube, or the unit hypertorus), and any value of p, including the special cases \(p=1\) and \(p=2\). In particular, in the two-dimensional case for \(p=2\) (see [1,2,3,4]), it has been proved that, as long as \(\Omega \) has unit volume, the leading term is \(\Omega \)-independent:

The main goal of the present paper is to show that, at least in the case \(n=2\), \(p=2\), the leading behaviour in N of \(E_{\Omega }(N)\) which is \(\Omega \)-dependent is a constant, which can be calculated exactly. More precisely, we are not able to establish a full perturbative expansion for \(E_{\Omega }(N)\), up to corrections o(1), for any \(\Omega \). Nonetheless, for all pairs \((\Omega , \Omega ')\), we predict that

with

where \(\{\lambda _i\}_{i \ge 1}\) is the set of eigenvalues of the Laplace–Beltrami operator on \(\Omega \) that are different from zero (if \(\Omega \) is a manifold with a boundary \(\partial \Omega \) of perimeter of order 1, it is the set of eigenvalues of the Laplace–Beltrami operator, with Neumann boundary conditions). That is, all terms in an asymptotic expansion of \(E_{\Omega }(N)\) which do not decrease with N must be ‘universal’.

One can notice that \(K_{\Omega }\) is a regularization of the trace of the Laplace–Beltrami operator. Another equivalent regularization is the so-called Robin mass \(R_\Omega ,\) see for instance [5, 6]. In particular Eq. (2) can be equivalently written as

The definition of \(R_{\Omega }\) is postponed to Sect. 3, while its relation with \(K_\Omega \) is described in Sect. 4.

Let us put this result in context, by summarising (part of) the state of the art for this problem. In contrast with the transportation problem for continuous measures, in this case the candidate optimal transportation maps are just the N! permutations, that is, the possible bijections between two sets of cardinality N, and in particular they are a finite set. For one given choice of the \(N^2\) weights \(w_{xy}\), the computational problem of finding one optimal map T,^{Footnote 1} and the associated cost is a well-studied problem, which turns out to be in the polynomial class [7,8,9]. Thus, the associated computational problem can be quickly solved.

This fact is in striking contrast with the problem in Probability Theory, of understanding the asymptotic of the average cost, on various domains \(\Omega \) and statistical ensembles for the red and blue point processes. This random version of the problem has attracted much attention both in Mathematics and Physics. In the Physics community, the interest has come from the analogy with ‘spin glasses’ in Statistical Mechanics, and a seminal contribution was given in the eighties by Orland [10], Mézard and Parisi [11], that considered the problem “in infinite dimension”, by introducing the so-called “random-link” approximation. This version of the problem was addressed using (non-rigorous) statistical physics techniques such as the replica theory and the cavity method [12]. Their original results were later put on rigorous ground [13,14,15].

The extension to finite-dimension of the random-link results is, however, quite challenging. A first attempt was carried on by Mézard and Parisi [16, 17] that showed how, for \(n>2\), the random-link result can be used as a zero-order approximation for the finite-dimensional solution, adding perturbatively a series of corrections. In the same years, a remarkable result was obtained by Ajtai and coworkers [18] for \(n=2\): they proved that, if the problem is considered on the unit square \(\Omega \equiv \mathcal {R}:=[0,1]^2\), then \(E_{\Omega }(N)=\mathcal O(\ln N)\).^{Footnote 2}

Recently, the forementioned result has been refined. In particular, by means of non-rigorous arguments, in Refs. [1, 2] it was claimed that, on the unit square \(\mathcal {R}\),

where \(c^\text {PP}_{\mathcal {R}}(N)=o(\ln N)\) (the factor 2 is for later convenience). This result has been later rigorously proved by Ambrosio and coworkers [3] and extended to any 2-dimensional compact closed manifold \(\Omega \) [4], where it is shown that

The latter paper also proves rigorous bounds for \(c^\text {PP}_\Omega \), namely that \(c^\text {PP}_\Omega (N)=\mathcal {O}(\sqrt{\ln N\ln \ln N})\). It has been recently conjectured that Eq. (6) holds also in the case of points generated from non-uniform densities [19].

In this paper we further investigate the problem of the estimation of \(c_\Omega (N)\). Extending the arguments given in [2], we argue that, on a generic two-dimensional manifold of unit area, the correction \(c_\Omega \) in Eq. (6) can be written as

where \(c_*^\text {PP}(N)\) does not depend on \(\Omega \). The index ‘PP’ is to denote that both the red and blue points are sampled with the ‘Poisson random process on \(\Omega \)’ (that is, are i.i.d. and uniform). Numerical investigations are compatible with the possibility that \(c_*^\text {PP}(N)\) is indeed a constant, and, under this hypothesis, we can give the constant \(c_*^\text {PP}\) the approximate value

Analogous claims and results hold for other variants of the problem, most notably when one set of points is still sampled with the Poisson random process, and the other set is either a deterministic regular grid (we investigate here the cases of square (S), triangular (T) or “Fibonacci” (F) [20, 21] grids), or, in the variant of the problem where T is the transportation between a discrete and a continuous measure, the uniform measure (U) on \(\Omega \). In these three new cases, the factor 2 in Eq. (6) disappears, and we have the similar structure

Let us stress again that the functions \(c_*\) are ‘universal’, in the sense that they do not depend on the choice of manifold \(\Omega \) (but they do depend on the choice of local randomness, e.g. among P, S, T, F, U), while the geometric correction \(K_{\Omega }\) depends on the choice of manifold, but is ‘universal’ in a different sense, as it is independent of the choice of local randomness (provided that the extra factor 2 in the P case is taken into account). Just as well as Eq. (2), such a decomposition is not at all granted a priori, and is somewhat surprising.

We also give numerical estimates of the associated values of \(c_*\),^{Footnote 3} under the hypothesis that they are indeed constant, namely

The paper is organized as follows. In Sect. 2 we define the random matching problems we are interested in. In Sect. 3 we present our functional approach for the derivation of the scaling of the optimal cost, including the finite-size corrections given in Eq. (7). For simplicity, we concentrate only on the Poisson–Poisson case. In Sect. 5 we apply our theory to different domains, giving an explicit computation of \(K_\Omega \) for all of them. In Sect. 6 we compare our predictions with numerical results obtained solving a large number of instances of the problem on the domains under investigation. Finally, in Sect. 7 we give our conclusions.

Let us consider a connected, two-dimensional smooth Riemannian manifold \(\Omega \) having finite volume and, if with a non-empty boundary, finite perimeter, with metric g. Given a system of local coordinates \((x^1,x^2)\) around a point \(p\in \Omega \), \(g=\sum _{ij}g_{ij}(p)\text {d}x^i\otimes \text {d}x^j\), and given two elements v, w in the tangent bundle in p, we will denote by \(\langle v,w\rangle _p :=\sum _{ij}g_{ij}(p) v_iw_j\). For the sake of generality, we will perform our analysis in the slightly more subtle case of \(\partial \Omega \ne \emptyset \) (the arguments below can be easily adapted to the case \(\partial \Omega =\emptyset \)). We will denote \(\text {d}\sigma \) the Riemannian measure for \(\Omega \), and we will assume the measure of \(\Omega \) to be equal to 1.^{Footnote 4} Moreover, we will denote by \(\delta _p(x)\text {d}\sigma (x)\) the unit measure concentrated in \(p\in \Omega \), that is, given a test function \(\varphi (x)\),

Suppose now that two sets of points are given on \(\Omega \), namely a set of N points , that we call the red points, and a set of N points , that we call the blue points. The assignment problem consists in assigning each red point to one blue point, in such a way that the resulting map T is a bijection, and a certain total cost function is minimized. As motivated in the introduction, the cost function is the sum of the costs for each pair \((X_i,Y_j)\) such that \(T(X_i)=Y_j\), and the cost of a pair (x, y) is the square of the Riemann distance between x and y of the selected pairs. In formulas, we have to find the optimal bijection such that

(10)

where

(11)

and d(x, y) is the Riemann distance between the points x and y, i.e., the infimum of the lengths of the curves that join the two points.

Note that each feasible T corresponds to a permutation \(\pi \) of N elements, so that \(T(X_i)=Y_{\pi (i)}\), and searching for the optimal map is equivalent to searching for the optimal permutation. If we introduce the two atomic measures

(12a)

(12b)

the optimal cost coincides with the 2-Wasserstein distance (squared) between the two empirical measures in Eq. (12), of which we shall recall here the definition (see e.g. [22]): given two probability measures \(\nu _1\) and \(\nu _2\) on \(\Omega \), their 2-Wasserstein distance is

where the infimum is taken over the set \(\Gamma (\nu _1,\nu _2)\) of all the joint probability distributions J with first and second marginal given by \(\nu _1\) and \(\nu _2\), respectively. It is well known (see for example [3, 23,24,25]) that, in our setting, the set of the optimal joint probability distributions J is a convex polytope, called Birkhoff polytope, whose extreme points are all and only the permutations \(\pi \) which are optimal within the probability distributions of the form \(J(x,y)=\sum _{i} \delta _{X_i}(x) \delta _{Y_{\pi (i)}}(y)\). Accordingly, the set of optimal maps \(T:\Omega \rightarrow \Omega \) pushing \(\nu _1\) to \(\nu _2\), i.e., those realising the infimum in the expression

coincides with the set of maps of the form \(T(X_i)=Y_{\pi (i)}\), with \(\pi \) optimal in the sense above. (The situation is much simpler when is absolutely continuous with respect to the Riemannian measure, as in this case the optimal transport map T would be unique. For a more complete discussion see also [24, ch. 9]).

The distance in Eq. (14) corresponds, up to a multiplicative constant, to the cost in Eq. (11) when and . Therefore

(15)

In the following, we will consider various statistical ensembles of pairs . At this point, many choices are possible. To be definite, we will always choose and to be independent of each other. We will also choose to be always what we shall call, with abuse of notation, the “Poisson random process on \(\Omega \) of size N”, that is, the \(Y_j\)’s are i.i.d., uniformly chosen on \(\Omega \) (w.r.t. the measure \(\text {d}\sigma \)).^{Footnote 5}

Poisson (P):

Also is given by a Poisson random process on \(\Omega \) of size N.

Uniform (U):

Together with the Poisson–Poisson, the most general and interesting case is the Uniform–Poisson case, in which the cost is the distance beetwen the Poisson random process and the uniform measure \(\text {d}\sigma \):

As a discrete approximation of this case, we can introduce various grid–Poisson assignment problem (GP), interesting by themselves:

Square grid (S):

when \(\Omega \) is a flat \(a \times b\) rectangular domain (possibly up to identification of the boundaries, e.g. as in a torus), and there exists a value k such that \(ka, kb \in \mathbb {N}\) and \(k^2ab=N\), then a natural choice is to fix to be the square grid of spacing 1/k. In the case of a torus, we can imagine identifying the horizontal sides of the fundamental rectangular region with a shift s. In this case the grid has no local defects when also \(ks \in \mathbb {N}\), and the modular parameter of the resulting surface is \(\tau =(s+i b)/a\), so that the set of points in the moduli space which can be realised by a grid with cardinality between N and \(N + N^{\frac{1}{2}+\epsilon }\) becomes dense everywhere in the limit of large N.

Triangular grid (T):

analogous to the square grid, in the case in which \(\Omega \) is a flat hexagon (possibly up to identification of the boundaries, e.g. as in a torus), with sizes (a, b, c, a, b, c) in cyclic order. Of course, this includes as special cases the regular triangle and hexagon, and the rhombus of angle \(\pi /3\). Now we require that there exists a value k such that \(ka, kb, kc \in \mathbb {N}\), and \(\frac{\sqrt{3}}{2} k^2(ab+bc+ca)=N\), and the natural choice is to fix to be the triangular grid of spacing 1/k. In the case of the torus, calling \(\omega =\hbox {e}^{2 \pi i/3}\), the associated modular parameter is \(\tau =\frac{\omega b + \omega ^2 c}{a + \omega ^2 c}\), so that also in this case the whole moduli space can be accessed by increasing N.

Fibonacci grid (F):

This is a subtle construction, adapted to the case of a sphere, see Fig. 1b, and based on stereographic projection from a Fibonacci spiral on the plane (from which the name), described in [20, 21]. The local aspect of this grid around one given point is somewhat intermediate between the one of a square and of a triangular lattice, with variations depending on the spherical coordinates of the point, and on the precise value of N. We will not enter in the detail of this construction, and the reader is referred to the forementioned papers.

a Fibonacci lattice on the unit disc. b Fibonacci lattice on the unit sphere. c Transportation between a set of \(N=144\) random blue points and a square grid of 9N red points on the flat torus (Color figure online)

Full size image

In Appendix A we give more details about the relation between grid-Poisson assignment problems and the UP problem, with an estimation of the convergence rate of the optimal cost in the former to the optimal cost in the latter.

As anticipated, we are interested in the study of the asymptotic behaviour in N of the average optimal transportation cost, for which we will adopt the general notation

(16)

where the average \({\mathbb E\left[ {\cdot }\right] }\) is taken over the pertinent statistical ensemble for the point processes.

As we said above, the study of \(E_{\Omega }(N)\) in any dimension \(n\ne 1\) or \(\infty \) (i.e., in the random-link model) seems rather difficult. A possible approach to the study of the asymptotic behaviour for large N is based on the fact that, in this limit, T(x) is expected to be ‘close’ to the identity map (more precisely, from [18] we expect that \(d(X_i,Y_j) =\mathcal O(\sqrt{(\ln N)/ N})\) for pairs of points which are paired by an optimal matching), and an expansion in this small parameter, at the first non-trivial order, might still capture the relevant features of the solution of the problem. In [1, 2, 26] this approach has been applied to the study of the problem when \(\Omega \) is the square, or the torus with modular parameter \(\tau =i\). The analysis leads to a result whose interpretation requires a regularization that takes into account the finite-N effects and avoid divergences, as we will see below.

In the approach in [1, 2, 26], a close analogy naturally emerges between the evaluation of the average optimal cost in the assignment problem and the evaluation of the electrostatic energy of 2N particles, N of each charge sign, pinned in random positions on \(\Omega \). This is a result a posteriori of the theory, as the obvious analogy just doesn’t hold as is (in the electrostatic problem, the energy is the sum of \(N(2N-1)\) pair contributions, not just N, which scale logarithmically with the distance of the pair, instead that quadratically). In a sense, the proposed linearization follows the opposite track of the suggestion by Born and Infeld [27] of a non-linear version of electrodynamics in order to solve the problem of divergencies. Similar ideas have been proposed recently by Brenier for fluid motion [28, 29].

3.1 Linearization

Let us review the arguments of [2], in their natural generalisation to a Riemannian manifold. We start by introducing, for each map \(T:\Omega \rightarrow \Omega \), the cost

(17)

Note that, at this stage, T is not a transportation map. First, because we haven’t still imposed the fact that the push-forward of is , and second, as specific to our transportation problem dealing with atomic measures, the ‘true’ optimal transportation map is only defined on the support of , which is not the whole \(\Omega \). We start by solving the first issue. An equivalent formulation of the constraint is that, for any function \(\phi :\Omega \rightarrow \mathbb {R}\), we must have

(18)

Again, it would be enough to consider functions \(\phi \) with the same support of , a fact of a certain relevance as it implies that, by expanding \(\phi \) over the appropriate basis of functions, we have only \(N-1\) independent constraints, instead that infinitely many, as it would be the case if were absolutely continuous with respect to the Lebesgue measure.

The idea is now to write down a Lagrangian that combines the cost expression in Eq. (11) with the condition in Eq. (18) as

(19)

where \(\phi \) plays the role of a Lagrange multiplier. The optimal map \(T^*\) satisfies the Euler–Lagrange equations obtained from the Lagrangian above (which turn out to be nonlinear).

We shall now use the idea that, for \(N\rightarrow +\infty \), we expect \(T(x)\rightarrow x\) for any \(x\in \Omega \), due to the fact that the matched pairs become infinitesimally close under the scaling in which \(|\Omega |\) is kept fixed. Then, there exists a vector field \(\mu (x)\) on \(\Omega \) such that, at the leading order

and

The direction of the field is the one of the geodesic curve realising the distance of T(x) from x. Pictorially

We shall now introduce

which is another perturbative parameter (when averaging over our statistical ensemble, monomials \(\mathbb {E}[\delta \nu (x_1) \delta \nu (x_2) \cdots \delta \nu (x_k)]\) have a definite scaling with N, and high powers are suppressed). The Lagrangian is approximated, in this limit, by its quadratic version,

where, by definition of gradient, \(\text {d}\phi (\mu )(x)=\langle \mu (x),\nabla \phi (x)\rangle _x\). We have also used the fact that, if the integrand is smooth enough, we can neglect the discrepancy between and \(\text {d}\sigma \) (while still treating more carefully \(\delta \nu (x)\)). Extremizing the new Lagrangian, and using that if \(\partial \Omega \ne \emptyset \) the field \(\mu \) is tangent to \(\partial \Omega \), we obtain the non-homogeneous linear equations

which are the linearization of the original Euler–Lagrange equations in the fields \(\mu \) and \(\phi \). In local coordinates:

where \(\partial _i\equiv \partial _{x^i}\), the tensor \(g^{ij}\) is the inverse of g and \(|g|=\det g\). The two equations imply the Poisson equation

to be solved with Neumann boundary conditions, since the flux of \(\mu (x)\) at the boundary is zero. Here \(-\Delta \) is the Laplace–Beltrami operator on \(\Omega \), i.e., in local coordinates

3.2 The Divergence of the Cost and the Problem of Regularization

The functional approach above tells us that \(\mu =-\nabla \phi \), where \(-\Delta \phi =\delta \nu \). We can use the fact that^{Footnote 6}

to write down an expression, valid for \(N\gg 1\), for a quantity \(\epsilon (x)\) that we shall call the cost density

so that \(E_{\Omega }(N)=\int \limits _\Omega \epsilon (x)\text {d}\sigma (x)\). Here we have introduced the Green function G(x, y) of \(-\Delta \) on the orthogonal complement of the locally constant functions. The Green function is a symmetric function that satisfies the equations

where \(\left. \partial _{n} G(x,y)\right| _{y\in \partial \Omega }\) is the normal derivative in x with respect to the boundary \(\partial \Omega \) of the domain. The equations above identify a unique Green function up to an additive constant: we will fix this constant adopting the convention

The obtained results have, however, a fundamental problem. The quantity in Eq. (29) is divergent for any \(x\in \Omega \). The responsibility for this fact comes from several sources, one of which is having treated the field \(\mu (x)\) as a continuous field, instead that a collection of N vectors, one per each point . This gives locally, in coordinates on the tangent space in \(X_i\), a field

(31)

where \(\hat{\mu }(x)\) is such that \(\hat{\mu }(X_i)\) is a finite quantity (here we have used the diagonal expression for the Green function G, see below Eq. (33)). Such an approximation could still be used through a delicate Cesàro limit, if we had to perform integrals in which \(\mu \) appears linearly. However our cost density is quadratic in these fields, and locally, at a formal level, for \(\delta \ll 1\),

(32)

that is, the appropriate result, which depends on the positions of the points and is finite at finite N, is shifted by a fixed but divergent quantity. Yet again, we observe a perfect analogy with 2-dimensional electrostatics, namely with the classical problem of the field self-energy for a distribution of point charges.

Analytically, we see this feature emerging from our result by observing that for \(d(x,y)\rightarrow 0\), the Green function behaves as [5, 6]

with a logarithmic divergence. We perform therefore a regularization of the logarithmic divergence, along the same lines of the classical treatment of electrostatics. Let us introduce \(\Omega _\delta (x)=\Omega \setminus B_\delta (x)\), where \(B_\delta (x)=\{y\in \Omega :d(x,y)<\delta \}\) is the ball of radius \(0<\delta \ll 1\) centered in x. We can introduce a regularized expression

and a corresponding “regularized cost”

(35)

where the second integral runs over the border \(\partial B_{\delta }(y):=\{x\in \Omega :d(x,y)=\delta \}\) of \(B_{\delta }(y)\), \(\text {d}\lambda (x)\) is the line element of \(\partial B_\delta (y)\) in x, and n is the outward normal to \(B_\delta \). By Eq. (30) and Eq. (30c) the first integral is infinitesimally small for \(\delta \rightarrow 0\). Therefore

(36)

For \(0<\delta \ll 1\), the inner integral can be estimated using the expression in Eq. (33), so that

(37)

We finally get

The integral of m(x)

is sometimes called Robin mass [5, 6].

Now, we suppose that the regularization by the parameter \(\delta \) acts in the same way on all geometries. Under this assumption we can compare two different geometries, \(\Omega \) and \(\Omega '\), obtaining the following conjecture.

Conjecture

Let \(\Omega \), \(\Omega '\) be two regular two-dimensional manifolds, then

In other words, the differences of the average cost among different manifolds, in the large-N limit, are expected to be regularization-independent, and, in addition to this, can be expressed in terms of the Robin’s masses of the Laplace–Beltrami Green function on \(\Omega \) and \(\Omega '\). The analytic evaluation of these differences will be the main object of our investigation, starting from Sect. 5. The remaining of this section is instead devoted to a further justification of the assumption at the basis of Eq. (40).

One problem at this point is that our regularization parameter \(\delta \) does not have a clear relation with the perturbative parameter \(N^{-1}\). In order to better understand what is the microscopic mechanism beyond the regularization, we observe that Eq. (38) can be formally written for \(\delta \rightarrow 0\) as

where the operator \(-\Delta ^{-1}\) is the inverse Laplace–Beltrami operator on \(\Omega \) (with Neumann boundary conditions, if the boundary exists)

As said above, a logarithmic divergence appears for \(\delta \rightarrow 0\) and both sides of Eq. (41) are infinite. By the Weyl law on the asymptotics of the eigenvalue counting function for the Laplace–Beltrami operator [30] we know that, for a 2-dimensional manifold with unit volume, and under Neumann boundary conditions, the leading behaviour of for large \(\lambda \) is^{Footnote 7}

(43)

Furthermore, the eigenfunction \(f_{\lambda }\) associated to a given value of \(\lambda \) ‘looks locally’ like a plane wave with wavelength \(1/\sqrt{\lambda }\). This has two consequences at the level of our approximations when passing from the complete Lagrangian, Eq. (19), to its quadratic approximation, Eq. (23). First, the Taylor expansion of \(\phi (T(x)) = \phi (x+\mu (x))\) around x, in the basis of the eigenfunctions \(\{f_{\lambda }\}\), is perturbative in the parameter \(\mathbb {E}(\mu (x)) \sqrt{\lambda }\), which we expect, from [18], to be of order \(\sqrt{\lambda \ln N /N}\). Second, if our basis is orthonormal for the measure \(\text {d}\sigma \), that is, \((f_{\lambda },f_{\rho })_{\text {d}\sigma }= \delta _{\lambda \rho }\) (assuming for simplicity of notation that the spectrum is non-degenerate), under the measure we get instead

(44)

The result of this analysis is that, if we decompose our fields in the basis of eigenfunctions of \(-\Delta \), we can neglect the corrections coming from the further terms in the Taylor expansion, and the discreteness of the measure, only for those eigenfunctions with \(\lambda \lesssim N\) (up to possible factors \((\ln N)^{\gamma }\) in the scaling). Conversely, in the regime \(\lambda \gtrsim N\) some unknown mechanism comes into play, and we expect that its effect is to dump the sum \(\hbox {tr}\Delta ^{-1}\) appearing in (41), possibly at a scale \(\lambda \lesssim N\). In [1], this unknown dumping mechanism is supposed to be encoded in a cut-off function \(F(\lambda /N)\). Now, as this function is related to the local expansion of the fields \(\mu \) and \(\phi \) at high frequencies, and as the relation between eigenvalue \(\lambda \) and local wavelength is universal, the function \(F(\lambda /N)\) must have one of the two flavours of universality: it shall not depend on the manifold \(\Omega \), while in general it must depend on the type of problem (among Poisson, various grids and uniform), that is, in a natural generalisation of the treatment of [1] to our setting, we should have some unknown functions \(F^{\bullet \text {P}} (\lambda /N)\), with \(\bullet \) being one among P, S, T, F or U, and no dependence on \(\Omega \). Going on with our analysis of the PP case (the reasoning can be repeated for all other cases similarly) and using \(F(\lambda /N)\) for \(F^\text {PP}(\lambda /N)\) for brevity, within the assumptions of [1] we should interpret the correspondence (41) above as

(45)

where \(\Lambda (\Omega )\) is the set of nonzero eigenvalues of the Laplace–Beltrami operator on \(\Omega \). Following the analysis already performed in [1], this gives

for any domain of unit measure, for some constant \(c_\Omega \) depending on the cut-off, which cannot be determined if F is not known. We recall that, as anticipated in the introduction, the leading term in Eq. (46) is the correct asymptotic cost, as rigorously proved in Refs. [3, 4, 33], the presence of a logarithm being known since the eighties [18].

Note that there is no guarantee that the cut-off function scales exactly as \(F(\lambda /N)\), as the mechanism beyond the dumping of the high-wavelength contributions, and the amount of this dumping, are not under control. It may well be, for example, that the function has the form \(F\big (\frac{\lambda }{N (\ln N)^\gamma }\big )\), which would give a variant of (46) in which instead of the constant term it will appear an universal term \(O(\gamma \ln \ln N).\)

All these arguments lead us to reformulate our conjecture as follows.

Conjecture

(Alternative formulation) Let \(\Omega \) be a regular two-dimensional manifold, then

where \(c_*(N){=o(\ln N)}\) is an universal function not depending on \(\Omega \).

Moreover, for \(\Omega \), \(\Omega '\) different regular manifolds,

Remark 1

It is important to remark that, from the simulations made in Sect. 4, we have evidence that the term \(c_*(N)\) can be chosen as a constant, i.e. that formula  (46) is compatible with our numerical results. Obviously it is not possibile to deduce that \(c_*(N)=c_*\) from numerical simulation only. Anyway, in case, we would get  (47) and \(c_{\Omega } = R_{\Omega }+c_*.\)

Remark 2

We notice that starting from Eq. (45) and comparing two manifolds, we get the interesting fact

(49)

The combination of the two integrals above can be rewritten as

(50)

The universality of Weyl law implies that the factor grows no faster than \(\sqrt{\lambda } \ln \lambda \), so that, even in absence of the function F (that is, in the limit of N large), the integral is convergent at infinity (and near zero is regularised by the spectral gap). This allows us to predict

(51)

The expression (45) is, annoyingly, a diverging expression depending on \(\Omega .\) A way of studying this expression is by introducing a regularization parameter \(\epsilon \) for these contributions, and then deducing an evaluation of (51) from a singular expansion in \(\epsilon \) around zero.

One standard way to perform this programme is the so-called zeta regularization [34]. Let us introduce the generating function

which is known to be absolutely convergent for \(\mathfrak {R}(s) > 1\), and in this case we recognise our scheme above under the identification \(s=1+\epsilon \). Then \(-\hbox {tr}\Delta ^{-1}\) can be regularized by looking at \(Z_{\Omega }(s)\) near \(s=1\) [35]

and by removing the pole at \(s=1\). That is, in equation (51),

For reasons that will appear clearer below, we will call Kronecker’s mass the constant \(K_\Omega \). Despite the fact that there seems to be no reason a priori to believe that \(K_{\Omega }\) and \(R_{\Omega }\) are related, it has been proved by Morpurgo that \(R_{\Omega }-K_{\Omega }\) is a universal constant (that is, it does not depend on \(\Omega \)), given by [5, 36,37,38]

where \(\gamma _{\text {E}}=0.57721\dots \) is the Euler–Mascheroni constant. In particular, this universality result is crucial in checking a posteriori that our two predictions (40) and (54), obtained by two different analyses, are consistent, and also implies that our Conjecture is equivalent to the statement of E. (51). The computation of the Kronecker’s mass is often easier than the Robin’s mass, as we will show below. For a few manifolds \(\Omega \), both computations, of \(R_{\Omega }\) and \(K_{\Omega }\), can be performed with relatively small effort, and we will do this, for pedagogical reasons, in order to illustrate with an example the forementioned general result.

Another way of performing our programme is to consider the regularized sum

for \(p>0\). This corresponds to a specific choice of function \(F(\lambda /N)\), provided that \(\epsilon \) is identified with \(\gamma \big /N^p\), for \(\gamma \) some constant. Also in this case we have, universally,

and this leads to the prediction

The analogue of the Morpurgo theorem reads in this case

as can be evinced by comparing the two regularizations for the diverging integral \(\frac{1}{4\pi }\int \nolimits _1^{\infty } \frac{\text {d}x}{x}\) (which, besides the fact that it is an integral rather than a sum, it has all the appropriate asymptotics properties implied by the Weyl law). Namely, for this choice we have

that is, \(K=0\), and

that is, \(W^{(p)}=-{\gamma _E}/{4p\pi }\).

Another regularization in the same spirit is through the regularized sums

with \(\theta (x)\) is the Heaviside step function, which, yet again, corresponds to a specific choice of function \(F(\lambda /N)\), provided that \(\epsilon \) is identified with \(\gamma /N\), for \(\gamma \) some constant. Also in this case we have a universal asymptotics

and this leads to the prediction

The analogue of the Morpurgo theorem reads in this case

as we have

that is, \(W^\text {sharp}=0\).

To verify our ansatz, we will compute the Kronecker’s mass and the Robin’s mass for different \(\Omega \). We will compare our analytic results with numerical simulations in Sect. 6. We will start considering flat manifolds having \(g(x)=\mathbb I\), and consider manifolds with uniform curvature starting from Sect. 5.5.

5.1 The Unit Rectangle

Let us start by considering the problem on the rectangle. We call \(\mathcal R(\rho )\) the rectangle \([0,\sqrt{\rho }] \times [0,1/\sqrt{\rho }]\), and we consider the Laplace–Beltrami operator with Neumann boundary conditions. The eigenfunctions of \(-\Delta \) on \(\mathcal R(\rho )\) are given by

The corresponding eigenvalues are

We proceed computing the Kronecker mass using the regularized function

Here we have adopted \(\tau =i\rho \), in compliance with standard notation for modular forms, and have introduced the lattice zeta function \(\zeta _\tau (s)\) defined in Appendix D. This calculation is readily performed thanks to a remarkable result due to Kronecker (and known as first limit formula of Kronecker), reported in Appendix D, Eq. (D.7), and that we repeat here:

where \(\eta (\tau )\) is the Dedekind \(\eta \) function. Kronecker’s formula allows us to immediately obtain^{Footnote 8}

that for \(\rho =1\) (unit square) simplifies to

We will see in the following that the first limit formula of Kronecker will allow us to extract the Kronecker’s mass for many types of flat domains: this explains our choice of ‘Kronecker mass’ for denoting \(K_\Omega \).

We can give a slightly more compact form to the function in (71), shifted by its minimal value above:

We shall make a remark on this expression. In the limit in which the rectangle is very elongated, we get^{Footnote 9}

that is the average cost for the Poisson–Poisson one-dimensional assignment problem on the segment of unit length [39]. This is not by accident. Indeed, in any rectangular domain we can evaluate the average energy of the permutation in which the k-th red point counting from the left is matched to the k-th blue point counting from the left. This configuration is optimal w.h.p. in the limit \(\rho \rightarrow +\infty \), and would be optimal, at any \(\rho \), if the vertical coordinates of all the points were equal. On the other side, a worst case is when all the vertical coordinates of red points are zero, and all the vertical coordinates of blue points are \(1/\sqrt{\rho }\), so that, calling \(E_{[0,1]}(N)\) the average energy for the 1-dimensional problem on the [0, 1] segment, we get

which, by substituting our scaling ansatz, gives

When we take a limit \(N \rightarrow \infty \), \(\rho \rightarrow \infty \) on a direction \(\rho \gg \sqrt{N}\), we thus get

which is consistent with (74). In the following we will encounter various other domains which allow a consistency check with a 1-dimensional limit. We will reach similar conclusions, without entering in the details of the estimates, as this is done by minor modifications of the reasonings presented here.

5.2 The Flat Torus

We shall now consider the problem on the flat torus \(\mathcal T(\tau )\). To describe the corresponding manifold, let us first consider the lattice of points on \(\mathbb {R}^2\)

generated by the matrix

(79)

corresponding to the base vectors

(80)

In such lattice it is possible to define fundamental parallelograms , containing no further lattice points in its interior or boundary. A fundamental parallelogram is given for example by

(81)

so that , see Fig. 2a. We will also use a shortcut adapted to rectangles,

a Pictorial representation of an assignment on a torus generated by quotient of \(\mathbb {R}^2\) with a periodic lattice, with fundamental parallelogram and the corresponding base vectors. b Contour plot of \( \mathfrak {I}(\tau ) |\eta (\tau )|^4\) in the complex plane \(\tau \). The triangoloid shape is the canonical fundamental region of the moduli space, given by \(|\tau | \le 1\), \(|\tau \pm 1| \ge 1\) and \(\mathfrak {I}(\tau )>0\)

Full size image

(82)

A torus \(\mathcal T\) is defined as a quotient between the complex plane and a lattice \(\Lambda \), \(\mathcal T:=\mathbb {R}^2/\Lambda \). In other words, each point is identified with the set of points \(\{x+\underline{\omega }\ n,\ n\in \mathbb {Z}^2\}\), the distance between two points in being the minimum distance between the elements of their equivalence classes. It is well known that two matrices \(\underline{\omega }\) and \(\underline{\omega }'\) identify the same lattice \(\Lambda \) and the same torus \(\mathcal T\) (although not the same fundamental domain ) if and only if \((\underline{\omega })^{-1} \underline{\omega }' \in \text {SL}(2,\mathbb {Z})\). For each \(\underline{\omega }\), we introduce the half-period ratio

Given a lattice \(\Lambda \) generated by \(\underline{\omega }\), it is possible to associate to it a dual lattice \(\Lambda ^*\) generated by \(\underline{\omega }^*\), such that \(\underline{\omega }^*\ \underline{\omega }^T=\mathbb I\), identity matrix, i.e.,

(84)

Each torus \(\mathcal T=\mathbb {R}^2/\Lambda \) is then naturally associated to a dual torus given by \(\mathcal T^*:=\mathbb {R}^2/\Lambda ^*\).

In the following, we will restrict, without loss of generality, to the case in which the fundamental parallelograms have unit area, choosing

(85)

such that \(\rho \in \mathbb {R}^+\) and \(\sigma \in \mathbb {R}\), and we will denote the corresponding torus by \(\mathcal T(\tau )\), where \(\tau :=\sigma +i\rho \) is the half-period ratio.

The Kronecker’s mass. Due to the periodicity conditions, the eigenfunctions of \(-\Delta \) on \(\mathcal T(\tau )\) have the form

for all \(k^*=\underline{\omega }^*\ k\in \Lambda ^*\), \(k=\left( {\begin{array}{c}n\\ -m\end{array}}\right) \in \mathbb {Z}^2\). The corresponding eigenvalue is

We can compute now the Kronecker mass using the regularized function

and removing the pole in \(s\rightarrow 1\), as discussed in Sect. 3. This calculation is readily performed, again thanks to the first limit formula of Kronecker, Eq. (D.7), which allows us to immediately obtain

In Fig. 2b we present a contour-plot of the related expression \(\mathfrak {I}(\tau ) |\eta (\tau )|^4\) in the complex plane \(\tau \), confined to the canonical fundamental region of the moduli space. In particular, this function diverges for \(\tau \rightarrow 0\) and has minimum at \(\tau =\exp (i \pi (\frac{1}{2} \pm \frac{1}{6}))\). This implies that, among all unit tori equipped with the flat metric, the “hexagonal” one, that is the one for which \(\tau \) is a sixth root of unity, is the one in which the average cost of the Euclidean Random Assignment Problem is minimised. More strikingly, as deduced from results in [5] which are in turn based on the results in [35], the hexagonal torus is minimal also among unit surfaces with non-uniform metric.

Example: the rectangular and rhomboidal tori. We shall call “rectangular torus” a torus in which the fundamental parallelogram is a rectangle. This case corresponds to \(\tau =i\rho \), with \(\rho >0\) real. Our formula specialises to

which is invariant under the map \(\rho \mapsto \rho ^{-1}\), as it should. In the region \(\rho \in (0,1]\) the lowest value is achieved at \(\rho =1\) (see also Fig. 2b), where

We can give a slightly more compact form to this function, shifted by its minimal value:

Similarly, we shall call “rhomboidal torus” a torus in which the fundamental parallelogram is a rhombus. This case corresponds to \(\tau =\hbox {e}^{i \theta }\), with \(0< \theta \le \pi \big /2\), and our formula specialises to

that is, again shifting by the value for the standard torus,

As was the case for the rectangle, the expression in Eq. (92), in the limit in which the torus is very “thin and long”, becomes

This happens to be the average cost for the Poisson–Poisson one-dimensional assignment problem on the circle of unit length [39], as was to be expected, by a reasoning analogous to the one presented for the case of the rectangle.

The Robin mass. Let us now evaluate, for the generic flat torus \(\mathcal T(\tau )\), the Robin mass \(R_{\mathcal T}(\tau )\). Calling \(z=z(x,y)= (x_1 -y_1) + i(x_2 - y_2)\), the Green’s function on the torus is given in this case by [40]

where \(\theta _1(z; \tau )\) is an elliptic \(\theta \) function. The Robin mass is obtained from

It is immediately seen that, in agreement with the Morpurgo theorem, Eq. (55) is satisfied.

5.3 Other Boundary Conditions on the Unit Rectangle

The unit rectangle and the rectangular torus are obtained starting from the fundamental domain in Eq. (82), and assuming respectively open (i.e., Neumann for the field \(\phi \)) and periodic boundary conditions. Other choices of boundary conditions are possible, which correspond to other classical surfaces, with or without boundary. Each choice leads to a different spectrum of the Laplacian, which in turn implies a different Kronecker mass (and, according to our theory, a different finite-size correction to the optimal cost of the assignment problem).

5.3.1 The Cylinder

Let us consider the domain and let us take periodic boundary conditions in the horizontal direction (i.e., the side of length \(\sqrt{\rho }\)) and Neumann boundary conditions in the vertical direction (i.e., the side of length \(1\big /\sqrt{\rho }\)), see Fig. 4a. The resulting surface is a cylinder, that we shall denote by \(\mathcal C(\rho )\). The eigenfunctions of \(-\Delta \) are the set of functions

The corresponding eigenvalues are therefore

Repeating the same type of calculations performed for the rectangle (that is, expressing the regularised sum as a combination of \(\zeta _{\tau }(s)\) (for some \(\tau \)’s) and \(\zeta (2s)\)), we obtain

so that

We also remark that

which is the cost density for the one-dimensional assignment problem with open boundary conditions (i.e., on the unit segment), while

which is the density of cost for the one-dimensional assignment problem with periodic boundary conditions (i.e., on the unit circle), again, as was to be expected. The nontrivial solution of Eq. \(K_{\mathcal C}(\rho )=K_{\mathcal C}\) is \(\rho =0.625352\dots \), while the minimum value of the mass occurs for \(\rho =0.793439\dots \)

Remark that the constants above do not appear in the study of [35], because in our context, in presence of a boundary, we should impose Neumann boundary conditions (while the authors of [35] only analyse the case of Dirichlet boundary conditions).

5.3.2 The Möbius Strip

Starting again from the rectangle , we can identify each point with all its images in \(\mathbb {R} \times [0,1\big /\sqrt{\rho }]\) generated by the map \((x,y)\rightarrow (x+\sqrt{\rho },1\big /\sqrt{\rho }-y)\). That is, if we see the surface as the fundamental rectangular domain , we impose open boundary conditions along the horizontal direction, and identify the two vertical sides after a ‘twist’ (that is, the top part on the left is glued to the bottom part on the right). The obtained domain \(\mathcal M(\rho )\) is called the Möbius strip, see Fig. 4b. The eigenfunctions of \(-\Delta \) are

The parity constraint implements the twisted identification of the strip. The corresponding eigenvalues are

Repeating now the usual arguments we get

so that

We also remark that

which is the average cost for the problem on the segment of length \(\frac{1}{2}\) (and the fact that the length is not 1 is related to the fact that the twisted boundary conditions are effectively ‘folding in two’ the segment), while

which is the cost for the problem on the unit circle. The non-trivial solution of Eq. \(K_{\mathcal M}(\rho ) = K_{\mathcal M}\) is found for \(\rho =4.1861\dots \), whereas the minimum is achieved at \( \rho =2.30422\dots \)

5.3.3 The Klein Bottle

If we identify both pairs of opposite sides of the rectangle, one pair (say, the horizontal sides) in the ordinary way, and the other pair in the twisted way as in the Möbius strip, we obtain the Klein bottle \(\mathcal K(\rho )\), see Fig. 4c. The eigenfunctions of \(-\Delta \) are in this case

and

Proceeding as above, one can obtain

so that, in particular,

We also remark that both one-dimensional limits coincide with the corresponding constructions for the Möbius strip, and indeed the limits of the analytical expressions are the same, as

while

Here \(K_{\mathcal K}(\rho ) = K_{\mathcal K}\) for \(\rho = 1.09673\dots \), whereas the minimum is obtained at \(\rho = 1.04689\dots \).

5.3.4 The Boy Surface

As a final example, let us take twisted boundary conditions for both pairs of opposite sides of . In this way we obtain the so-called Boy surface \(\mathcal B(\rho )\), see Fig. 4d. The eigenfunctions of \(-\Delta \) are

The calculation proceeds as in the other cases, giving

which is symmetric for \(\rho \leftrightarrow \frac{1}{\rho }\), as it must be. In particular

Now, both one-dimensional limits produce a domain corresponding to the segment of length \(\frac{1}{2}\), and indeed

5.4 The Disc and the Cone

Up to now, we have mostly solved the problem using the zeta regularization of the Laplacian, and relying on Kronecker’s first limit formula. Only for the case of the torus, we have also performed the calculation of the Robin mass, and verified the prediction of the Morpurgo theorem. In this section we will give the results for a geometry \(\Omega \) in which the calculation of the Robin mass is done with relatively small effort, as the Green function can be guessed through the method of images, while the calculation of the Kronecker mass would require a sum over maxima and minima of Bessel functions.^{Footnote 10} Let us introduce the notation \(\mathcal D_p(r)\) for the circular sector of radius r and angle \(\frac{2\pi }{p}\), see Fig. 6a,

The unit area condition implies \(2\pi r^2=p\). We considered the case \(p\in \mathbb {N}\), and we choose periodic boundary conditions in the angular direction: this is equivalent to say that we identify the two radii of the sector, obtaining in this way a cone of height \(r\sqrt{1-p^{-2}}\), see Fig. 6b. This surface is interesting, as it is the first example in our list of a surface with singular curvature, the conical singularity being at the vertex of the cone. We have argued in the introduction that, because of the scaling \(\sim \sqrt{N^{-1}\ln N}\) of the field \(\mu \), we expect that the same theory applies to the flat Euclidean space and to curved manifolds, as long as the curvature is non-singular. The case of surfaces with a finite number of conical singularities would require a different (although feasible) argument, and the verification of our theory on this family of surfaces (as well as the surfaces treated in Sect. 5.5) is an important validation of our predictions.

The Robin’s mass for this case is obtained in Appendix B, and it is equal to

where \(\phi (z)\) is the digamma function. In particular, for \(\alpha \rightarrow 2\pi \) we recover the case of the unit disc \(\mathcal D\equiv \mathcal D_{1}\):

The Kronecker’s mass is readily obtained using Eq. (55).

5.5 The Unit Sphere and the Real Projective Sphere

An example of transportation problem on the surface of the sphere \(\mathcal S^2\) has been considered in [42], where the problem of transporting a uniform mass distribution into a set of random points on \(\mathcal S^2\) is analyzed. As an example of applications of our approach to non-flat manifolds, here we consider the problem in our usual setting, i.e., a transportation between two atomic measures of random points. As in the previous cases, the information on the finite-size corrections is related to the spectrum of the Laplace–Beltrami operator on the manifold. It is well known that the eigenfunctions of \(-\Delta \) on the surface of a sphere of radius r are the spherical harmonics \(Y_{l,m}(\theta ,\phi )\) with \(l\in \mathbb {N}\) and \(m\in \mathbb {Z}\) with \(-l\le m\le l\). The corresponding eigenvalues are

that is, the eigenvalue \({l(l+1)}/{r^2}\) has multiplicity \(2l+1\), for the range of integers \(-l \le m \le l\).

We fix the surface area of the sphere to 1 (taking \(r=(4\pi )^{-1\big /2}\)) and we proceed using the zeta regularization, computing

In this case, after some algebra, we are led to use ‘just’ the version of the zeta regularization for the Riemann zeta function (which is much simpler than the Kronecker formula)

and we obtain

The Kronecker mass for the unit sphere is therefore

Alternatively, we can use one of the regularizations illustrated in Sect. 4. A convenient one is the evaluation of \(W^{(1\big /2)}\), which gives

Recalling that in our case \(r^2=(4\pi )^{-1}\), the final result is

that, in light of (59), allows to rederive equation (127).

The spherical lune. The calculation above can be extended to the spherical lune \(\mathcal S_k^2\), a surface on a sphere of radius r, \(4\pi r^2=k\), contained by two half great circles which meet at antipodal points with dihedral angle \(\frac{2\pi }{k}\), see Fig. 7. In Appendix C it is shown that the Kronecker’s mass corresponding to this manifold is

Choosing periodic boundary conditions on the two half great circles, the Kronecker’s mass takes an additional \(-\frac{1}{2\pi }\ln 2\) contribution, and we obtain

that reduces to (127) for \(k=1\), as it should.

The projective sphere. A variation of the problem on the (unit) sphere is the problem on the (unit) real projective sphere \(\mathcal {PS}^2\), that is, the sphere in which antipodal points are identified. The eigenfunctions of the Laplace–Beltrami operator are still the spherical harmonics \(Y_{l,m}(\theta ,\phi )\) with \(l\in \mathbb {N}\) and \(m\in \mathbb {Z}\), \(-l\le m\le l\), but we have to restrict ourselves to eigenfunctions that are invariant under the transformation \((\theta ,\pi ) \mapsto (\pi - \theta ,\phi + \pi )\), i.e., to even values of l. Working on the unit-area manifold accounts to have \(2 \pi r^2=1\). We get

so that the Kronecker’s mass is

Using a sharp cut-off instead

which, again by using (63) and (65), allows to rederive equation (133).

We have numerically investigated all the cases described above to check our predictions. To solve the assignment problem we have used an implementation of the Jonker-Volgenant algorithm [8]. For each domain \(\Omega \), we have computed the expected optimal cost averaging over at least \(10^4\) independent instances and different sizes N of the system, \(32\le N\le 1024\). In each case, we have estimated \(c^{\bullet \text {P}}_\Omega (N)\) assuming that they are indeed constant at the leading order, i.e., \(c^{\bullet \text {P}}_\Omega (N)\equiv c^{\bullet \text {P}}_\Omega \), via a least square regression. Here \(\bullet =\text {\{}P,S,T,F,U\}\). Our results are given in Table 1.

For each domain we have also computed \(c_*^{\bullet \text {P}}=c_\Omega ^{\bullet \text {P}}-K_\Omega \) that we expect to be domain-independent. In the PP case, our best estimation of \(c_*^\text {PP}\), obtained for the PP problem on the torus (see Fig. 3a) is

Numerically, however, we cannot rule out a weak N-dependence in \(c_*^\text {PP}\). In a similar way we have obtained the results given in (9), that we repeat here,

To verify our results in Eqs. (8) we have also proceeded in this way. Given two unit-area domains \(\Omega \), \(\bar{\Omega }\), and a given type of the problem (PP, SP, etc), we have computed, for each N, \(\Delta E^{\bullet \text {P}}_{\Omega ,\bar{\Omega }}(N)=c^{\bullet \text {P}}_{\Omega }(N)-c^{\bullet \text {P}}_{\bar{\Omega }}(N)\) and then extrapolated to \(N\rightarrow +\infty \). If the arguments above are correct, then

In Fig. 3b we plot in particular (Fig. 4)

a Numerical estimation for \(c_{\mathcal T}(N)\) for different values of N. Each data point is obtained averaging the optimal cost over at least \(10^6\) instances and then removing the leading \(\frac{1}{2\pi }\ln N\) term. The fit is obtained using a quadratic function in \(1\big /\sqrt{N}\). b Difference of average optimal costs for the assignment on the rectangle \(\mathcal R(\rho )\), on the torus \(\mathcal T(i\rho )\) and on the Boy surface \(\mathcal B(\rho )\) with the corresponding costs for \(\rho =1\). The numerical results, represented by the dots, are compared with the analytical prediction obtained from Kronecker’s masses

Full size image

Pictorial representation of the different boundary conditions considered in Sect. 5.3, with an example of assignment at \(N=3\) for each case. To obtain the corresponding surface, the red edges have to be considered joined in such a way that the directions of the arrows match (Color figure online)

Full size image

as functions of \(\rho \) for \(\bullet =\text {\{}P,S\}\). Similarly, in Fig. 5 we present our results for (the absolute value of)

Absolute difference of average optimal costs for the assignment on the cylinder \(\mathcal C(\rho )\), on the Möbius strip \(\mathcal M(\rho )\) and on the Klein bottle \(\mathcal K(\rho )\) with the corresponding costs for \(\rho =1\). The numerical results, represented by the dots, are compared with the analytical prediction obtained from Kronecker’s masses

Full size image

for \(\bullet =\text {\{}P,S\}\).^{Footnote 11} In Fig. 6c and in Fig. 7b we have also considered the differences of average optimal costs in the case of the circular sector and of the spherical lune. In all investigated cases we found a perfect agreement with our predictions.

Random points on a circular sector (a) and optimal assignment on the corresponding cone (b) (in red, the part of the domain boundary where periodic boundary conditions are imposed). c Comparison between numerical results and our theoretical prediction for \(E^\text {PP}_{\mathcal D_p}-E^\text {PP}_{\mathcal D}\) (Color figure online)

Full size image

a Optimal assignment on a spherical lune. b Comparison between numerical results and our theoretical prediction for \(E^\text {PP}_{\mathcal S^2_k}-E^\text {PP}_{\mathcal S^2}\) with periodic (smooth line) and Neumann (dashed line) boundary conditions.

Full size image

In this paper we have considered the assignment problem between two sets of random points on a generic two-dimensional smooth manifold of unit area. We have showed, by means of analytical arguments and numerical simulations, that the average optimal cost can be written as

if both sets of points are random (Poisson–Poisson case), and as

if one of the two sets is fixed on a grid (square, triangular, Fibonacci) or replaced with the uniform measure. In the equations above, \(K_\Omega \) is a precise quantity that can be obtained by a zeta-regularization of the trace of the inverse Laplace–Beltrami operator on \(\Omega \). The contributions \(c^{\bullet \text {P}}_*\) are instead independent on \(\Omega \) and related to the ‘local details’ of the problem (i.e., if the assignment is between random points, or with a grid, or with the uniform measure). We have given an exact computation of \(K_\Omega \) for different domains, and using different procedures.

The quantity \(c_*^{\bullet \text {P}}(N)\) shows a weak dependence on N (if no dependence at all): it has been proven indeed that \(c_*^\text {UP}=\mathcal {O}(\sqrt{\ln N\ln \ln N})\) [4], a bound that, because of triangular inequality, holds for all the cases that we have considered. Assuming that \(c_*^{\bullet \text {P}}(N)\) are constants, we have verified, within our numerical precision, their independence on \(\Omega \) in all considered transportation cases (Poisson–Poisson, grid–Poisson, uniform–Poisson).

Our results reduce the computation of the (leading) finite-size correction to the optimal cost to the computation of the \(\Omega \)-independent contributions \(c_*^{\bullet \text {P}}(N)\). These contributions are intrinsically dependent on the local nature of the problem (and therefore change if, for example, we fix on a grid one of the two sets of points) and are inherited by the regularization of the highest part of the spectrum of \(-\Delta \), as discussed in Sect. 3. What are the properties (and possibly the exact value, if they are constant) of \(c_*^{\bullet \text {P}}(N)\) remains an open question.

Finally, analogous results are expected to hold in higher dimension at the leading order. In particular, the analysis in [1] suggests that, for \(d > 2\), the local properties of the problem affect the coefficient of the leading term, with a finite-size correction depending on the spectrum of the Laplacian only. This case may also be investigable with our tools, as indeed versions of the Weyl law, which we crucially use in Sect. 3.2, exist in generic dimension. We leave the investigation of the higher dimensional problem for future works.