Finite-size corrections in the random assignment problem

TL;DR

This paper analytically derives finite-size corrections to the average optimal cost in the random assignment problem using the replica formalism, revealing how these corrections vary with different cost distributions.

Contribution

It provides the first analytical derivation of finite-size corrections for a broad class of cost distributions in the assignment problem, including power-law, Gamma, and delta functions.

Findings

Corrections change sign and scaling when moving from power-law to Gamma distributions.

Numerical solutions confirm the analytical predictions.

Behavior near delta-function distribution is analyzed.

Abstract

We analytically derive, in the context of the replica formalism, the first finite size corrections to the average optimal cost in the random assignment problem for a quite generic distribution law for the costs. We show that, when moving from a power-law distribution to a $\Gamma$ distribution, the leading correction changes both in sign and in its scaling properties. We also examine the behavior of the corrections when approaching a $\delta$-function distribution. By using a numerical solution of the saddle-point equations, we provide predictions that are confirmed by numerical simulations.

Figures

Tables

Table 1. Table 1: Numerical values of the rescaled corrections appearing in Eqs. ( 13 ) for different values of r 𝑟 r .

$r$	${\eta }_0^{\frac{1}{r+1}}{\hat{E}}_r$	${\left(N{\eta }_0^2\right)}^{\frac{1}{r+1}}\frac{{\eta }_0\Delta {\hat{F}}_r^{\left(1\right)}}{{\eta }_1}$	$N{\eta }_0^{\frac{1}{r+1}}\Delta {\hat{F}}_r^T$	$N{\eta }_0^{\frac{1}{r+1}}\Delta {\hat{F}}_r^F$
-0.5	1.125775489	-2.777285153	-3.917446075	-1.192663973
-0.4	1.334614017	-2.952484269	-3.665262242	-1.250475151
-0.3	1.471169704	-2.921791666	-3.324960744	-1.222990786
-0.2	1.558280634	-2.784084499	-2.984917100	-1.157857158
-0.1	1.612502443	-2.600804197	-2.675513663	-1.079610016
0	1.644934067	-2.404113806	-2.404113806	-1
0.1	1.662818967	-2.215821874	-2.168528577	-0.924257491
0.2	1.671039856	-2.038915744	-1.966713438	-0.854501434
0.3	1.672729262	-1.877696614	-1.792481703	-0.791231720
0.4	1.670005231	-1.732453452	-1.641566768	-0.734262435
0.5	1.664311154	-1.602337915	-1.510248399	-0.683113178
0.6	1.656639222	-1.486024319	-1.395391897	-0.637204338
0.7	1.647677145	-1.382051819	-1.294397704	-0.595951390
0.8	1.637905005	-1.288993419	-1.205124002	-0.558807473
0.9	1.627659755	-1.205532353	-1.125808312	-0.525279810
1	1.617178636	-1.130489992	-1.054997763	-0.494933215
2	1.519733739	-0.670341811	-0.626403698	-0.303146650
3	1.446919560	-0.461144035	-0.431759755	-0.211631545
4	1.393163419	-0.346056113	-0.324185048	-0.159938240
5	1.352087648	-0.274505368	-0.257174804	-0.127356338
6	1.319651066	-0.226200326	-0.211931870	-0.105200594
7	1.293333076	-0.191617643	-0.179566694	-0.089276830
8	1.271505390	-0.165752490	-0.155385461	-0.077340947
9	1.253073980	-0.145742887	-0.136697943	-0.068095120
10	1.237277174	-0.129842072	-0.121861122	-0.060741591

Table 2. Table 2: Numerical results for the average optimal cost for different values of r 𝑟 r and theoretical predictions. The value E ^ r P superscript subscript ^ 𝐸 𝑟 P \hat{E}_{r}^{\text{P}} from Ref.  Houdayer  et al. ( 1998 ) , due to a different convention adopted in that paper, is obtained as E ^ r P = ( 2 ​ π r + 1 2 ​ Γ ​ ( r + 1 ) Γ ​ ( r + 1 2 ) ​ Γ ​ ( r + 2 ) ) 1 r + 1 ​ β num ​ ( r + 1 ) superscript subscript ^ 𝐸 𝑟 P superscript 2 superscript 𝜋 𝑟 1 2 Γ 𝑟 1 Γ 𝑟 1 2 Γ 𝑟 2 1 𝑟 1 subscript 𝛽 num 𝑟 1 \hat{E}_{r}^{\text{P}}=\left(\frac{2\pi^{\frac{r+1}{2}}\Gamma(r+1)}{\Gamma\left(\frac{r+1}{2}\right)\Gamma(r+2)}\right)^{\frac{1}{r+1}}\beta_{\text{num}}(r+1) from Table II therein. The data for E ^ r P superscript subscript ^ 𝐸 𝑟 P \hat{E}_{r}^{\text{P}} from Ref.  Lee and Orlin ( 1993 ) have been obtained via a linear fit, using a fitting function in the form of Eq. ( 105 ).

$r$	${\hat{E}}_r^{\text{P}}$	${\lambda }_r{\hat{E}}_r^{\Gamma }$	${\hat{E}}_r^{\text{P}}$ Houdayer et al. (1998)	${\hat{E}}_r^{\text{P}}$ Lee and Orlin (1993)	Th. prediction
-0.5	4.5011(3)	4.504(1)	–	–	4.503101957
-0.4	3.12611(5)	3.1268(2)	–	–	3.126825159
-0.3	2.4484(1)	2.4488(3)	–	–	2.448788557
-0.2	2.0593(5)	2.0593(3)	–	–	2.059601452
-0.1	1.8127(3)	1.8126(2)	–	–	1.812767212
0	1.64500(5)	1.6449(2)	1.645(1)	1.6450(1)	1.644934067
0.1	1.5245(2)	1.5253(9)	–	–	1.524808331
0.2	1.4356(2)	1.4357(5)	–	–	1.435497487
0.3	1.3670(1)	1.3670(4)	–	–	1.367026464
0.4	1.31323(6)	1.3132(3)	–	–	1.3132296
0.5	1.27007(8)	1.2697(4)	–	–	1.270107121
0.6	1.2350(1)	1.2348(3)	–	–	1.234960167
0.7	1.20585(6)	1.2062(6)	–	–	1.205907312
0.8	1.18143(3)	1.1812(7)	–	–	1.181600461
0.9	1.16099(8)	1.1605(6)	–	–	1.161050751
1	1.14344(7)	1.1433(4)	1.143(2)	–	1.14351798
2	1.05371(1)	1.054(1)	1.054(1)	1.054(1)	1.053724521
3	1.02311(1)	1.0288(9)	1.0232(1)	1.0236(2)	1.023126632
4	1.009690(4)	1.010(3)	1.0098(1)	–	1.009736514
5	1.00303(2)	1.005(3)	1.00306(8)	1.0026(8)	1.003027802

Equations

$$E(\pi ):=\sum _{e\in \pi }{w}_e$$

$$ϵ(\pi ):=\frac{1}{\nu (\mathcal{G})}\sum _{e\in \pi }{w}_e.$$

$$E({\pi }^{\ast })=\underset{\pi \in \mathcal{M}}{\min }E(\pi ).$$

$${\pi }_{ij}=\{\begin{array}{ll}1& \text{for }e=(i,j)\in \pi \\ 0& \text{otherwise},\end{array}$$

$$\sum _{i=1}^N{\pi }_{ij}=\sum _{i=1}^N{\pi }_{ji}=1\forall j\in [N].$$

$$E(\pi )=\sum _{i=1}^N\sum _{j=1}^N{\pi }_{ij}{w}_{ij}.$$

$$\overline{E}:=\overline{E({\pi }^{\ast })}=\overline{\underset{\pi \in \mathcal{M}}{\min }\sum _{i=1}^N\sum _{j=1}^N{\pi }_{ij}{w}_{ij}},$$

$$Z[w]=\sum _{\pi }\left[\prod _{j=1}^N\delta \left(1-\sum _{i=1}^N{\pi }_{ij}\right)\delta \left(1-\sum _{i=1}^N{\pi }_{ji}\right)\right]e^{-\beta E(\pi )},$$

$$\overline{F}$$

$$\overline{E}$$

$${\rho }_r(w)=w^r\sum _{k=0}^{\infty }{\eta }_k(r)w^k,r>-1,{\eta }_0(r)\ne 0.$$

$${\hat{E}}_r=\underset{N\to \infty }{\lim }\frac{1}{N^{\frac{r}{r+1}}}{\overline{E}}_r=\frac{r+1}{{\left[{\eta }_0\Gamma (r+1)\right]}^{\frac{1}{r+1}}}{J}_r^{(r+1)}$$

$$J_r^{(\alpha )}:={\int }_{-\infty }^{+\infty }{\hat{G}}_r(-u)D_u^{\alpha }{\hat{G}}_r(u)du$$

$${\hat{G}}_r(l)={\int }_{-l}^{+\infty }\frac{(l+y{)}^r}{\Gamma (r+1)}{e}^{-{\hat{G}}_r(y)}dy$$

$${\hat{E}}_r(N)={\hat{E}}_r+\sum _{k=1}^{[r]+1}\Delta {\hat{F}}_r^{(k)}+\Delta {\hat{F}}_r^T+\Delta {\hat{F}}_r^F+o\left(\frac{1}{N}\right),$$

$$\Delta {\hat{F}}_r^{(k)}=$$

$$\Delta {\hat{F}}_r^T=$$

$$\Delta {\hat{F}}_r^F=$$

$$\Delta {\hat{F}}_r^{(1)}=-\frac{{\eta }_1}{N^{\frac{1}{r+1}}}\frac{r+1}{{\eta }_0[{\eta }_0\Gamma (r+1){]}^{\frac{2}{r+1}}}{J}_r^{(r)}.$$

$${\rho }_r^{\Gamma }(w):=\frac{w^r{e}^{-w}\theta (w)}{\Gamma (r+1)},$$

$${\rho }_r^{\text{P}}(w):=(r+1)w^r\theta (w)\theta (1-w),$$

$${\eta }_k^{\Gamma }(r)=\frac{1}{\Gamma (r+1)}\frac{(-1{)}^k}{k!},k\ge 0,$$

$${\eta }_k^{\text{P}}(r)=(r+1){\delta }_{k,0},k\ge 0.$$

$$\overline{\ln Z}=\underset{n\to 0}{\lim }\frac{\overline{Z^n}-1}{n}.$$

$${\int }_0^{2\pi }{e}^{ik\lambda }d\lambda =2\pi {\delta }_{k0}.$$

$$\begin{gathered} Z^n[w]=\left[\prod _{a=1}^n\prod _{i=1}^N{\int }_0^{2\pi }\frac{d{\lambda }_i^a}{2\pi }{\int }_0^{2\pi }\frac{d{\mu }_i^a}{2\pi }{e}^{i({\lambda }_i^a+{\mu }_i^a)}\right]\times \\ \times \prod _{i=1}^N\prod _{j=1}^N\prod _{a=1}^n\left[1+e^{-i({\lambda }_j^a+{\mu }_i^a)-\beta w_{ij}}\right]. \end{gathered}$$

$$\begin{gathered} \prod _{a=1}^n\left[1+e^{-i({\lambda }_j^a+{\mu }_i^a)-\beta w_{ij}}\right] \\ =\sum _{\alpha \in \mathcal{P}([n])}{e}^{-\beta |\alpha |w_{ij}-i{\sum }_{a\in \alpha }({\lambda }_j^a+{\mu }_i^a)} \\ =1+\sum _{p=1}^n{e}^{-\beta p{w}_{ij}}\sum _{\begin{array}{c}\alpha \in \mathcal{P}([n])\\ |\alpha |=p\end{array}}{e}^{-i{\sum }_{a\in \alpha }({\lambda }_j^a+{\mu }_i^a)}, \end{gathered}$$

$$g_{\alpha }\equiv g_p:={\int }_0^{+\infty }{\rho }_r(w)e^{-\beta pw}dw.$$

$$\underset{N\to \infty }{\lim }\frac{1}{N^{\frac{r}{r+1}}}\overline{F}=$$

$$\underset{N\to \infty }{\lim }\frac{1}{N^{\frac{r}{r+1}}}\overline{E}=$$

$$\beta =\hat{\beta }{N}^{\frac{1}{r+1}},$$

$$\begin{gathered} {\hat{g}}_p:=N{g}_p=N{\int }_0^{+\infty }{\rho }_r(w)e^{-p\hat{\beta }{N}^{\frac{1}{r+1}}w}dw \\ =\sum _{k=0}^{+\infty }\frac{1}{N^{\frac{k}{r+1}}}\frac{{\eta }_k\Gamma (k+r+1)}{(\hat{\beta }p{)}^{k+r+1}}. \end{gathered}$$

$$\begin{gathered} \overline{Z^n}=\left[\prod _{a=1}^n\prod _{i=1}^N{\int }_0^{2\pi }\frac{d{\lambda }_i^a}{2\pi }{\int }_0^{2\pi }\frac{d{\mu }_i^a}{2\pi }{e}^{i({\lambda }_i^a+{\mu }_i^a)}\right]\prod _{i,j}^N\left(1+\frac{T_{ij}}{N}\right) \\ =\left[\prod _{a=1}^n\prod _{i=1}^N{\int }_0^{2\pi }\frac{d{\lambda }_i^a}{2\pi }{\int }_0^{2\pi }\frac{d{\mu }_i^a}{2\pi }{e}^{i({\lambda }_i^a+{\mu }_i^a)}\right]\times \\ \times \exp \left[\frac{1}{N}\sum _{i=1}^N\sum _{j=1}^N\left(T_{ij}-\frac{T_{ij}^2}{2N}\right)+o\left(\frac{1}{N^2}\right)\right], \end{gathered}$$

$$T_{ij}:=\text{sideset}{}^{\prime }{\sum }_{\alpha \in \mathcal{P}([n])}{\hat{g}}_{\alpha }{e}^{-i{\sum }_{a\in \alpha }({\lambda }_j^a+{\mu }_i^a)}$$

Full text

I Introduction

Matching is a classical problem in combinatorial optimization Papadimitriou and Steiglitz (1998); Lovász and Plummer (2009). It can be defined on any graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is its set of vertices and $\mathcal{E}$ its set of edges. A matching on $\mathcal{G}$ is a set of non adjacent edges of $\mathcal{G}$, that is, a set of edges that do not share a common vertex. A matching is maximal when the addition of any new edge to the set makes it no longer a matching. A matching is said to be maximum when it has the maximal cardinality among the maximal matchings. All the maximum matchings have the same cardinality $\nu(\mathcal{G})$, which is called the matching number of the graph $\mathcal{G}$. A perfect matching (or 1-factor) is a matching that matches all vertices of the graph. That is, every vertex of the graph is incident to exactly one edge of the matching. Every perfect matching is maximum and hence maximal. A perfect matching is also a minimum-size edge cover. We will denote by $\mathcal{M}$ the set of perfect matchings.

Suppose now that we can assign a cost $w_{e}\geq 0$ to each edge $e\in\mathcal{E}$. For each perfect matching $\pi\in\mathcal{M}$ we define a total cost (energy)

[TABLE]

and a mean cost for edge

[TABLE]

The matching problem consists in the choice of the optimal perfect matching $\pi^{*}\in\mathcal{M}$ which minimizes the total cost

[TABLE]

Let us now consider the matching problem on the complete bipartite graph $\mathcal{K}_{N,M}$, in which the vertex set is the union of two disjoint sets $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ of cardinality $N$ and $M$, respectively. Defining $[n]\coloneqq\{1,\dots,n\}$, we identify $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ with $[N]$ and $[M]$, respectively, and the edge set is therefore the set of all couples $e=(i,j)$ with $i\in[N]$ and $j\in[M]$. In this case the matching problem is usually called the (general) assignment problem and the matching number is $\nu(\mathcal{K}_{N,M})=\min\{N,M\}$. In the following we will concentrate on the $N=M$ case, in which a perfect matching $\pi$ is a permutation in the symmetric group $\mathcal{S}_{N}$ and can be represented by a square matrix with entries $\pi_{ij}\in\{0,1\}$ for all $i\in[N]$ and $j\in[N]$ such that

[TABLE]

with the constraints

[TABLE]

The matching cost associated with $\pi$ can be written as

[TABLE]

From the point of view of computational complexity, matching problems are simple problems, being in the P complexity class, as Kuhn (1955) proved with his celebrated Hungarian algorithm for the assignment problem. Very fast algorithms are nowadays available both to find perfect matchings and to solve the matching problem on a generic graph Edmonds (1965); Edmonds and Karp (1972); Micali and Vazirani (1980); Lovász and Plummer (2009).

The properties of the solution of a matching problem on a given ensemble of realizations are often interesting as well Steele (1997). In the random assignment problem, for example, we consider the matching problem on $\mathcal{K}_{N,N}$, whereas the costs for all edges are independent random variables, identically distributed with a common law $\rho$. Each particular choice $\mathcal{W}=\{w_{e}\}_{e\in\mathcal{E}}$ for the set of edge costs is called an instance of the problem. In this random version of the problem, we are interested in the typical properties of the optimal matching. In particular we will concentrate on the asymptotic behavior for large $N$ of the average optimal cost

[TABLE]

where we have denoted by an over bar the average over all possible instances (i.e., the average over the disorder). In the same way, we can consider the matching problem with random weights on the complete graph $\mathcal{K}_{2N}$, having $2N$ vertices such that each one of them is connected to all the others. We simply call this variation of the problem random matching problem. Both the random matching problem and the random assignment problem have been solved by Mézard and Parisi (1985) by means of the replica trick. The random assignment problem and the random matching problem have also been generalized to the Euclidean case, in which the weights in $\mathcal{W}$ are functions of the distances between points associated with the vertices of the graph and the points are assumed to be randomly generated on a certain Euclidean domain Shor (1985); Mézard and Parisi (1988); Yukich (1998); Lucibello et al. (2017). Due to the underlying Euclidean structure, dimensionality plays an important role in the scaling of the optimal cost of random Euclidean matching problems Mézard and Parisi (1988); Caracciolo et al. (2014), and correlation functions can be introduced and calculated Boniolo et al. (2014); Caracciolo and Sicuro (2015a). Euclidean matching problems proved to be deeply connected with Gaussian stochastic processes Boniolo et al. (2014); Caracciolo and Sicuro (2014) and with the theory of optimal transport Caracciolo and Sicuro (2015b). In the latter context, Ambrosio et al. (2016) rigorously derived the asymptotic behavior of the average optimal cost for the two-dimensional random Euclidean assignment problem, previously obtained in Ref. Caracciolo et al. (2014) using a proper scaling ansatz. For a recent review on random Euclidean matching problems, see Ref. Sicuro (2017).

Remarkably enough, after the seminal works of Kirkpatrick et al. (1983), Orland (1985), and Mézard and Parisi, the application of statistical physics techniques to random optimization problems proved to be extremely successful in the study of the typical properties of the solutions, but also in the development of algorithms to solve a given instance of the problem Mézard and Montanari (2009); Bapst et al. (2013). In formulating a combinatorial problem as a model in statistical mechanics, an artificial inverse temperature $\beta$ is introduced to define a Boltzmann weight $\exp\left(-\beta E\right)$ for each configuration. Of course, configurations of minimal energy are the only ones to contribute in the limit of infinite $\beta$. For example, in the assignment problem, the corresponding partition function for each instance is

[TABLE]

where the “energy” $E(\pi)$ is given by (6). Thermodynamic information is obtained from the average total free energy

[TABLE]

In this paper we apply the formalism above to the random assignment problem, where the costs of all the edges are taken to be independent and identically distributed random variables with probability distribution density $\rho_{r}(w)$ such that, in the neighborhood of the origin, $\rho_{r}$ can be written as

[TABLE]

In the previous expression, $\eta_{k}(r)$ are coefficients (possibly dependent on $r$) of the Maclaurin series expansion of the function $\rho_{r}(w)w^{-r}$, which is supposed to be analytic in the neighborhood of the origin. The constraint $r>-1$ is required to guarantee the integrability of the distribution near the origin. By the general analysis performed in Refs. Orland (1985); Mézard and Parisi (1985), which we will resume in Section II, the average cost, in the asymptotic regime of an infinite number $N$ of couples of matched points, will depend on the power $r$ that appears in Eq. (11) only, aside from a trivial overall rescaling related to $\eta_{0}$. More precisely, if $\overline{E}_{r}$ is the average optimal cost obtained using the law $\rho_{r}$, then

[TABLE]

and it is analytically known for $r=0$ and, after a proper rescaling of its variable, in the $r\to\infty$ limit.

Our main results concern the finite-size corrections to the average optimal costs, and they will be presented in Section III, extending the classical achievements in Refs. Mézard and Parisi (1987); Parisi and Ratiéville (2002). In particular, we obtain the expansion

[TABLE]

In our discussion, we will consider in particular two probability distribution densities, namely, the Gamma distribution

[TABLE]

defined on $\mathds{R}^{+}$, and the power-law distribution

[TABLE]

defined on the compact interval $[0,1]$. In the previous expressions we have denoted by $\theta(w)$ the Heaviside theta function on the real line. Observe that, for the distribution $\rho_{r}^{\Gamma}$, we have

[TABLE]

whereas in the case of $\rho_{r}^{\text{P}}$,

[TABLE]

The case $r=0$ has already been considered by Mézard and Parisi (1987) and subsequently revised and corrected by Parisi and Ratiéville (2002). In the case analyzed in their works, the contributions $\Delta\hat{F}_{0}^{(1)}$, $\Delta\hat{F}_{0}^{T}$, and $\Delta\hat{F}_{0}^{F}$ are of the same order. This is not true anymore for a generic distribution with $r\neq 0$. As anticipated, a relevant consequence of our evaluation is that, if $\eta_{1}\neq 0$, for $r>0$ the most important correction comes from $\Delta\hat{F}_{r}^{(1)}$ and scales as $N^{-\frac{1}{r+1}}$. It follows that, in order to extrapolate to the limit of an infinite number of points, the best choice for the law for random links (in the sense of the one that provides results closer to the asymptotic regime) is the pure power law $\rho_{r}^{\text{P}}$, where only analytic corrections in inverse power of $N$ are present. Such a remark is even more pertinent in the limit when $r\to\infty$ at a fixed number of points, where the corrections $\Delta\hat{F}_{r}^{(k)}$ become of the same order of the leading term. Indeed, the two limits $r\to\infty$ and $N\to\infty$ commute only if the law $\rho_{r}^{\text{P}}$ is considered.

The paper is organized as follows. In Section II we review, in full generality, the calculation of the replicated partition function of the random assignment problem. In Section III we evaluate the finite-size corrections, discussing the different contributions and proving Eqs. (13). In Section IV we evaluate the relevant $r\to\infty$ case, pointing out the noncommutativity of this limit with the thermodynamic limit. In Section V we provide the numerical values of the necessary integrals and we compare our prediction with a Monte Carlo simulation for different values of $r$. In Section VI we summarize and give our conclusions.

II The replicated action

In the present section we perform a survey of the classical replica computation for the random assignment problem, following the seminal works of Mézard and Parisi (1985, 1987) (for a slightly different approach see also Ref. Orland (1985)), but we do not adopt their choice to replace $\beta$ with $\beta/2$. As anticipated in Section I, the computation of the average of $\ln Z$ goes through the replica trick Edwards and Anderson (1975)

[TABLE]

In other words, in order to compute $\overline{\ln Z}$ we introduce $n$ noninteracting replicas of the initial system, denoted by the index $a\in[n]$. For each $i\in[N]$, $2n$ replicated fields $\{\lambda_{i}^{a},\mu_{i}^{a}\}_{a=1,\dots,n}$ appear to impose the constraints in Eq. (5), using the relation

[TABLE]

We obtain

[TABLE]

Let ${\cal P}([n])$ be the set of subsets of the set $[n]$ and for each subset $\alpha\in{\cal P}([n])$ let $|\alpha|$ be its cardinality. Then

[TABLE]

where we have extracted the contribution from the empty set in the sum, which is $1$, and we have partitioned the contribution from each subset of replicas in terms of their cardinality. This expression is suitable for the average on the costs. From the law $\rho_{r}$ we want to extract the leading term for large $\beta$ of the contribution of each subset $\alpha\in{\cal P}([n])$ with $|\alpha|=p$. In particular, we define

[TABLE]

Due to the fact that short links only participate in the optimal configuration, approximating $\rho_{r}(w)\sim\eta_{0}w^{r}$, we obtain that the minimal cost for each matched vertex is of the order $N^{-\frac{1}{r+1}}$ so the total energy $E$ and the free energy should scale as $N^{\frac{r}{r+1}}$, that is, the limits

[TABLE]

are finite. This regime can be obtained by considering in the thermodynamic limit

[TABLE]

where $\hat{\beta}$ is kept fixed. As a consequence we set

[TABLE]

The replicated partition function can be written therefore as

[TABLE]

with

[TABLE]

where in the sum $\sum^{\prime}$ on subsets the empty set is excluded. If we introduce, for each subset $\alpha\in{\cal P}([n])$, the quantities

[TABLE]

we can write

[TABLE]

We perform now a Hubbard–Stratonovich transformation, neglecting $o(N^{-2})$ terms in the exponent in Eq. (27), obtaining

[TABLE]

up to higher order terms in the exponent. Now, let us observe that

[TABLE]

Introducing the function of $v_{\alpha}$,

[TABLE]

and the order parameters

[TABLE]

we can write

[TABLE]

a form that is suitable to be evaluated, in the asymptotic limit for large $N$, by means of the saddle-point method. It is immediately clear that $\Delta S^{T}$ contains a contribution to the action that is $O(N^{-1})$ and therefore it can be neglected in the evaluation of the leading contribution. It follows that the stationarity equations are of the form

[TABLE]

The application of the saddle point method gives

[TABLE]

where $\bm{\Omega}$ is the Hessian matrix of $S[Q]$ and $Q^{\text{sp}}$ is the saddle-point solution. As we will show below, the contribution $\ln\det\bm{\Omega}[Q^{\text{sp}}]$ provides finite-size corrections to the leading contribution of the same order of the corrections in $N\Delta S^{T}[Q^{\text{sp}}]$.

II.1 Replica symmetric ansatz and limit of vanishing number of replicas

To proceed with our calculation, we adopt, as usual in the literature, a replica symmetric ansatz for the solution of the saddle-point equations. A replica symmetric solution is of the form

[TABLE]

In particular this implies that $Y_{\alpha}=0$. In order to analytically continue to $n\to 0$ the value at the saddle-point of $S$ in Eq. (35b), let us first remark that under the assumption in Eq. (38)

[TABLE]

Moreover, as shown in Appendix A,

[TABLE]

where

[TABLE]

In conclusion, under the replica symmetric ansatz in Eq. (38), the functional to be minimized is

[TABLE]

A variation with respect to $q_{k}$ gives the saddle-point equation

[TABLE]

which is to say

[TABLE]

This implies that

[TABLE]

These formulas are for a general law $\rho_{r}$. Observe also that the expression of $\hat{g}_{p}$ is not specified. For finite $r$ and $N\to\infty$, Eq. (26) simplifies as

[TABLE]

We will restrict the analysis to the case in which Eq. (46) holds. Then Eq. (44) becomes

[TABLE]

with

[TABLE]

In Eq. (47), and in the following, we introduce the subindex $r$ to stress the dependence of $G$ and of the thermodynamical functionals on $r$. The average cost is therefore

[TABLE]

Using the fact that (see Appendix B)

[TABLE]

if we introduce

[TABLE]

in the limit $\hat{\beta}\to+\infty$, the function $\hat{G}_{r}$ satisfies Eq. (12c) and the value of $\hat{E}_{r}$ is the one reported in Eq. (12a). In particular, at fixed $r$, if we consider the two laws $\rho_{r}^{\text{P}}$ and $\rho_{r}^{\Gamma}$, the ratio between the corresponding average optimal costs is given by

[TABLE]

In the case $r=0$, we have the classical result by Mézard and Parisi (1985)

[TABLE]

a result that was later obtained with a cavity approach by Aldous (2001). For the evaluation of the integral, see Appendix D.

III Finite size corrections

The evaluation of the first-order corrections for a finite number of points has been considered in Refs. Mézard and Parisi (1987); Parisi and Ratiéville (2002) in the $r=0$ case. For this particular choice and assuming a distribution law $\rho^{\Gamma}_{0}$, a much stronger conjecture was proposed by Parisi (1998) and later proved by Linusson and Wästlund (2004); Nair et al. (2005), that is, for every $N$,

[TABLE]

For large $N$, Parisi’s formula implies

[TABLE]

Using instead the law $\rho_{0}^{\text{P}}$ (uniform distribution on the interval) we have Mézard and Parisi (1987); Parisi and Ratiéville (2002)

[TABLE]

from which we see that corrections for both laws are analytic, with the same inverse power of $N$, but different coefficients.

In their study of the finite-size corrections, the authors of Ref. Parisi and Ratiéville (2002) show that, in their particular case, there are two kind of finite-size corrections. The first one comes from the application of the saddle-point method, giving a series of corrections in the inverse powers of $N$. This contribution is the sum of two terms. The first term in this expansion corresponds to the contribution of the $\Delta S^{T}$ term given in Eq. (35c) appearing in the exponent in Eq. (37). The second term is related to the fluctuations, also appearing in Eq. (37), involving the Hessian of $S$. The second kind of corrections, instead, is due to the particular form of the law $\rho_{r}(w)$ for the random links and in particular to the series expansion in Eq. (26). This contribution can be seen at the level of the action $S$ in Eq. (35b), being

[TABLE]

In full generality, the expansion of $1/\hat{g}_{\alpha}$ generates a sum over terms each one of order $N^{-\frac{k}{r+1}}$ with $k\geq 1$. All these corrections are $o(N^{-1})$ for $r\in(-1,0)$, whereas the corrections obtained from the contributions with $1\leq k\leq r+1$ are of the same order as the analytic term, or greater, for $r\geq 0$. In particular, if $\eta_{1}\neq 0$, for $r>0$ the $k=1$ term provides the leading correction, scaling as $N^{-\frac{1}{r+1}}$. It is also evident that all these corrections are absent if $\eta_{k}=0$ for $k\geq 1$, as it happens in the case of the $\rho_{r}^{\text{P}}$ law.

III.1 Correction due to $\eta_{1}$

Let us consider the $r\geq 0$ case and let us restrict ourselves to the $k=1$ term, of order $N^{-\frac{1}{r+1}}$ in Eq. (13). Its contribution to the total free energy is given by

[TABLE]

where we already made a replica symmetric assumption and considered the $n\to 0$ limit. Imposing the saddle-point relation in Eq. (43) and using the limit in Eq. (46), we obtain

[TABLE]

To put the expression above in the form presented in Eq. (13e), observe that

[TABLE]

a structure that can be more useful for numerical evaluation, at least for $r$ integer. In this equation we have used Eq. (119) and we have introduced the Riemann–Liouville fractional derivative

[TABLE]

where $\Omega\coloneqq(-\infty,t)$ is the domain of integration (see the monographs in Refs. Samko et al. (1993); Podlubny (1999) for further details).

III.2 Correction due to the saddle point approximation

Let us now consider the corrections due to the saddle point approximation. The first contribution is expressed by $\Delta S^{T}$, given in Eq. (35c). In the replica symmetric hypothesis, we have that

[TABLE]

We can write the corresponding correction to the free energy as

[TABLE]

In Appendix E we show that the previous quantity can be written as

[TABLE]

Another type of finite-size correction comes from the fluctuations around the saddle point (Parisi and Ratiéville, 2002, Section B.3), related to the Hessian matrix $\bm{\Omega}$ appearing in Eq. (37). The evaluation of the contribution of the Hessian matrix is not trivial and it has been discussed by Mézard and Parisi (1987) and later by Parisi and Ratiéville (2002). They proved that the whole contribution comes from a volume factor due to a non trivial metric $\hat{\bm{\Omega}}$ obtained from $\bm{\Omega}$ imposing the replica symmetric assumption and such that

[TABLE]

The $n\times n$ matrix $\hat{\bm{\Omega}}$ can be written as

[TABLE]

where $\mathbb{I}_{n}$ is the $n\times n$ identity matrix and we have introduced the quantities

[TABLE]

and $\bm{\Pi}$ is a projection matrix on the constant vector defined as

[TABLE]

where $\mathbb{J}_{n}$ is the $n\times n$ matrix with all entries equal to $1$. The matrix $\bm{\Pi}$ has one eigenvalue equal to $1$ and $n-1$ eigenvalues equal to [math]. It follows that, because the two matrices $\bm{\Pi}$ and $\mathbb{I}_{n}$ obviously commute, $\hat{\bm{\Omega}}$ has one eigenvalue equal to $a_{0}+(n-1)a_{1}$ and $n-1$ eigenvalues equal to $a_{0}-a_{1}$. Its determinant is therefore simply given by

[TABLE]

In the limit of $n\to 0$ we easily get

[TABLE]

for all values of $r$. Similarly,

[TABLE]

so

[TABLE]

Therefore,

[TABLE]

In conclusion, integrating by parts and using Eq. (121), we obtain

[TABLE]

III.3 Application: the $r=0$ case

The results obtained in the $r=0$ case, analyzed by Parisi and Ratiéville (2002), can be easily recovered. From the general expression in Eq. (13), by setting $r=0$, we get

[TABLE]

where we have used the results discussed in the Appendix D for the two integrals involved in the expression above. Eqs. (56) and (57) are obtained using Eqs. (16) and (17), respectively.

IV The limiting case $r\to+\infty$

In this section we concentrate on the limiting case in which $r\to+\infty$. We can easily verify that, in the weak sense,

[TABLE]

so all the weights become equal to unity. We expect therefore that

[TABLE]

independently of $N$. The average cost obtained using $\rho^{\Gamma}_{r}$ instead diverges and it is therefore more interesting to consider the modified law

[TABLE]

According to our general discussion, we have that

[TABLE]

implying that, independently of $N$,

[TABLE]

and therefore

[TABLE]

in the limit of infinite $N$. In particular,

[TABLE]

It follows that, even though the two laws $\rho^{\text{P}}$ and $\rho^{\gamma}$ both converge to the same limiting distribution, according to our formulas, the corresponding average costs are not the same. This is due to the fact that the two limits $N\to+\infty$ and $r\to+\infty$ do not commute for the law $\rho^{\gamma}$, because of the presence of $O(N^{-\frac{k}{r+1}})$ corrections that give a leading contribution if the $r\to\infty$ limit is taken first.

To look into more details in the $r\to+\infty$ limit, we find it convenient, when looking at the saddle-point solution, to perform a change of variables, following the approach in Refs. Houdayer et al. (1998); Parisi and Ratiéville (2001), that is, writing

[TABLE]

then Eq. (12c) becomes

[TABLE]

so, in the $r\to+\infty$ limit

[TABLE]

If we set $\mathcal{G}_{\infty}(x)=a\operatorname{e}^{x}$, with

[TABLE]

we recover

[TABLE]

From Eq. (12a) with the change of variable in Eq. (84), we get

[TABLE]

so

[TABLE]

in agreement with the previous results.

Let us now evaluate the integrals appearing in the finite-size corrections in Eq. (13). Let us first start with the $\Delta F^{T}_{r}$ and the $\Delta F_{r}^{F}$ corrections. From the definition, for large $r$,

[TABLE]

where

[TABLE]

and therefore for large $r$,

[TABLE]

so we get

[TABLE]

a contribution that vanishes as $r^{-1}$ both for the law $\rho^{\text{P}}_{r}$ and for the law $\rho^{\gamma}_{r}$ (indeed we know that, in this case, all corrections must vanish when $r\to\infty$ at fixed $N$). In particular, if we consider the law $\rho^{\Gamma}_{r}$, we have, in the $r\to\infty$ limit,

[TABLE]

Similarly, for large $r$, we have that

[TABLE]

and therefore

[TABLE]

Instead, if we consider $\Delta\hat{F}^{(1)}_{r}$, we have that

[TABLE]

finally obtaining

[TABLE]

so that, considering the law $\rho^{\gamma}_{r}$, if we send $r\to\infty$ before taking the limit $N\to\infty$, $\Delta\hat{F}_{+\infty}^{(1)}\sim O(1)$ and we get a new contribution to the average optimal cost

[TABLE]

a series where we miss the contributions of order $N^{\frac{k}{r+1}}$ for $k\geq 2$, and that we know will sum to $1$.

V Numerical results

In this section we discuss some numerical results. First, we present a numerical study of our theoretical predictions obtained in the previous sections. Second, we compare with numerical simulations, in which the random assignment problem is solved using an exact algorithm.

The evaluation of all quantities in Eq. (13) depends on the solution of Eq. (12c). We solved numerically this equation for general $r$ by a simple iterative procedure. In particular, for $r>0$ we generated a grid of $2K-1$ equispaced points in an interval $\left[-y_{\text{max}},y_{\text{max}}\right]$ and we used a discretized version of the saddle-point equation in Eq. (12c) in the form

[TABLE]

We imposed as the initial function $\hat{G}_{r}^{[0]}$ of the iterative procedure

[TABLE]

We observed that the quantity

[TABLE]

decays exponentially with $s$ and therefore convergence is very fast. For our computation, we used typically $30$ iterations.

For $r<0$ the term $(l+y)^{r}$ in the saddle-point equation is divergent in $y=-l$ and Eq. (101) cannot be adopted. We have therefore rewritten the saddle-point equation using an integration by parts, obtaining

[TABLE]

After discretizing the previous equation, we used the same algorithm described for the $r\geq 0$ case (for a discussion on the uniqueness of the solution of Eq. (12c), see Ref. Salez (2015)).

In Table 1 we present our numerical results for the quantities

[TABLE]

for different values of $r$. Observe that the quantities appearing in the expansion in Eqs. (13) can be calculated using these values for any $\rho_{r}$ at given $r$, in addition to simple prefactors depending on the chosen distribution $\rho_{r}$.

In order to test our analysis for correction terms, we performed direct Monte Carlo sampling on a set of instances. Previous simulations have been reported, for example, in Refs. Mézard and Parisi (1985); Brunetti et al. (1991); Houdayer et al. (1998); Lee and Orlin (1993). In our setting, each realization of the matching problem has been solved by a C++ implementation of the Jonker-Volgenant algorithm Jonker and Volgenant (1987).

We first evaluated the asymptotic average optimal costs $\hat{E}_{r}^{\text{P}}$ and $\hat{E}_{r}^{\Gamma}$, obtained, for different values of $r$, using the laws $\rho_{r}^{\text{P}}$ and $\rho_{r}^{\Gamma}$, respectively. In the case of the law $\rho_{r}^{\text{P}}$, the asymptotic estimate for $\hat{E}_{r}^{\text{P}}$ has been obtained using the fitting function

[TABLE]

with $\alpha_{r}^{\text{P}}$ and $\beta_{r}^{\text{P}}$ fitting parameters to be determined, $\alpha_{r}^{\text{P}}$ corresponding to the value of the average optimal cost in the $N\to\infty$ limit. For a given value of $r$, we averaged over $I_{N}$ instances for each value of $N$ accordingly with the table below.

[TABLE]

Similarly, the asymptotic average optimal cost $\hat{E}_{r}^{\Gamma}$ has been obtained using a fitting function in the form

[TABLE]

We adopted therefore a three-parameter fitting function, constructed according to Eq. (13) including the finite-size correction up to $o(N^{-1})$ for $r\geq 0$ and up to $O(N^{-2})$ for $\frac{1}{2}\leq r<2$. As in the case before, the asymptotic estimation for $\hat{E}_{r}^{\Gamma}$ is given by $\alpha^{\Gamma}_{r}$. Our data were obtained extrapolating the $N\to\infty$ limit from the average optimal cost for different values of $N$. The investigated sizes and the number of iterations were the same adopted for the evaluation of $\hat{E}_{r}^{\text{P}}$. To better exemplify the main differences in the finite-size scaling between the $\rho_{r}^{\text{P}}$ case and the $\rho_{r}^{\Gamma}$ case, we have presented the numerical results for $r=1$ in Fig. 1. In the picture, it is clear that the asymptotic value $\hat{E}_{1}^{\text{P}}=\frac{1}{\sqrt{2}}\hat{E}_{1}^{\Gamma}$ is the same in the two cases, as expected from Eq. (52), but the finite-size corrections are different both in sign and in their scaling properties.

In Table 2 we compare the results of our numerical simulations with the ones in the literature (when available) for both $\hat{E}_{r}^{\text{P}}(N)$ and $\lambda_{r}\hat{E}_{r}^{\Gamma}(N)$, $\lambda_{r}$ being defined in Eq. (52). In Fig. 2 we plot our theoretical predictions and the numerical results that are presented in Table 2.

Let us now consider the finite-size corrections. In the case of the $\rho_{r}^{\text{P}}$ law, the $O(N^{-1})$ corrections are given by $\Delta\hat{F}_{r}^{T}+\Delta\hat{F}_{r}^{F}$ and no nonanalytic corrections to the leading term appear. We obtain the finite-size corrections from the data used for Table 2, using Eq. (105) but fixing $\alpha_{r}^{\text{P}}$ to the average optimal cost $\hat{E}_{r}^{\text{P}}$ given by the theoretical prediction in Table 1 and therefore with one free parameter only, namely, $\beta_{r}^{\text{P}}$. In Fig. 3 we compare our predictions for $\Delta\hat{F}_{r}^{T}+\Delta\hat{F}_{r}^{F}$, deduced by the values in Table 1, with the results of our numerical simulations for different values of $r$.

In the case of the $\rho_{r}^{\Gamma}$ law with $r>0$, the first correction to the average optimal cost is given by $\Delta\hat{F}_{r}^{(1)}$, whereas $\Delta\hat{F}_{r}^{(1)}$ is $o(N^{-1})$ for $r<0$. Again, this correction can be obtained by a fit of the same data used to extrapolate the average optimal cost, fixing the fitting parameter $\alpha_{r}^{\Gamma}$ in Eq. (106) to the theoretical prediction, and performing a two parameters fit in which the quantity $\gamma_{r}^{\Gamma}$ appearing in Eq. (106) corresponds to $\Delta\hat{F}_{r}^{(1)}$. In Fig. 4 we compare our prediction for $\Delta\hat{F}_{r}^{(1)}$, given in Table 1, with the results of our fit procedure for $\gamma_{r}^{\Gamma}$ for $-\frac{1}{2}<r\leq 5$. Observe that the numerical evaluation of the single contribution $\Delta\hat{F}_{r}^{(1)}$ is not possible for $r=0$. In this case, the result of our fit for the $O(N^{-1})$ correction was $\beta^{\Gamma}_{r}+\gamma_{r}^{\Gamma}=-0.97(4)$, to be compared with the theoretical prediction $N(\Delta\hat{F}^{(1)}_{0}\!+\!\Delta\hat{F}^{F}_{0}\!+\!\Delta\hat{F}^{T}_{0})=-0.998354732\dots$.

VI Conclusions

In the present paper we have discussed the finite-size corrections in the random assignment problem for a generic distribution law $\rho_{r}(w)$ for the weights in the form of Eq. (11). We have shown that, for $r>0$ and $\eta_{1}\neq 0$, the first finite-size correction scales as $O(N^{-\frac{1}{r+1}})$ and it is proportional to $\eta_{1}$. In particular, the integrals $J_{r}^{(0)}$, $J_{r}^{(r)}$ and $J_{r}^{(r+3)}$ are positive quantities (see Appendix C). Therefore, independently of $r$, the corrections $\Delta\hat{F}_{r}^{T}$ and $\Delta\hat{F}_{r}^{F}$ are negative, while $\Delta\hat{F}_{r}^{(1)}$ has opposite sign with respect to $\eta_{1}$, so that, for example, it is positive for the law $\rho_{r}^{\Gamma}$, while it vanishes for the law $\rho_{r}^{\text{P}}$. We also provided a general expression for the coefficients of the $O(N^{-\frac{1}{r+1}})$ and $O(N^{-1})$ corrections. Moreover, we have shown that, if $\lim_{r\to+\infty}\rho_{r}=\rho$, then in general

[TABLE]

We have finally numerically verified our results, by a numerical integration of our formulas and a comparison with simulations.

The $O(N^{-\frac{k}{r+1}})$ corrections appearing in Eq. (13), for $2\leq k\leq[r]+1$, remain to be computed. As discussed above, in the $r\to\infty$ limit, it is expected that all these finite size corrections contribute to the leading term, justifying the noncommutativity of the limits in Eq. (107).

VII Acknowledgments

The authors thank Giorgio Parisi for many discussions and suggestions. The authors also thank Carlo Lucibello and Marco Gherardi for discussions on the numerical solution of the integral equations involving $\hat{G}_{r}$. M.P.D. thanks the LCM for the extensive use of their computational resources. The work of G.S. was supported by the Simons Foundation (Grant No. 454949, Giorgio Parisi).

Appendix A Evaluation of $z[Q]$ on the saddle point and analytic continuation for $n\to 0$

Let us evaluate now the quantity $z[Q]$ on the saddle point. Using the fact that, for any analytic function $f$,

[TABLE]

we can write

[TABLE]

where $\bm{\alpha}=\{\alpha_{i}\}_{i}$ and $\alpha_{i}\in{\cal P}([n])$ are disjoint subsets whose union is $[n]$; however,

[TABLE]

To perform the analytic prolongation, we prove now that, if $f(0)=1$, then

[TABLE]

This fact can be seen observing that, for $n\to 0$,

[TABLE]

By the change of variable $s=\operatorname{e}^{l}$, Eq. (111) follows.

Appendix B Asymptotic behaviour of the function $B_{r}$

In this appendix we study the asymptotic behavior for large $\lambda$ of the function $B_{r}(\lambda x)$. By definition in Eq. (48)

[TABLE]

so that

[TABLE]

where we have used the fact that

[TABLE]

is the Bessel function of zeroth order of the first kind.

Appendix C Some properties of the function $\hat{G}_{r}$

In this appendix we give some properties of the function $\hat{G}_{r}$, defined by the integral equation (12c). From the definition, we have that, for $0\leq\alpha<\beta+1$ and $r>-1$,

[TABLE]

Observe that

[TABLE]

In this equation we have used the fact that, for $0\leq\alpha<\beta+1$, we have Podlubny (1999)

[TABLE]

In particular, for $\alpha=r$ we have the simple relation

[TABLE]

Moreover, for $0\leq\alpha<r+1$,

[TABLE]

From Eq. (119)

[TABLE]

The relations above imply that

[TABLE]

Similarly, for $0<k<r+1$ an integer,

[TABLE]

For large $l$ we have

[TABLE]

As anticipated, an exact solution is available in the $r=0$ case. In particular, for $r=0$, the second derivative

[TABLE]

is an even function of $l$,

[TABLE]

with the constant $c=1$ by evaluating the left-hand side in the limit of infinite $l$ and

[TABLE]

Then we have that

[TABLE]

which means that

[TABLE]

where we have used the initial condition at $l=0$, that is, because of Eq. (128),

[TABLE]

Appendix D Evaluation of the integrals in the $r=0$ case

To explicitly evaluate some of the integrals above, let us introduce the polygamma function

[TABLE]

which satisfies the recursion relation

[TABLE]

which, for a positive integer argument and assuming $m\geq 1$, leads to

[TABLE]

For $m=0$ this implies

[TABLE]

with $\gamma_{E}$ is Euler’s gamma constant and

[TABLE]

are the harmonic numbers.

With these considerations in mind and using Eq. (131), we have that

[TABLE]

Then we compute

[TABLE]

We remark now that

[TABLE]

Applying now the identity

[TABLE]

discovered by Euler, we recover the result obtained by Parisi and Ratiéville (2002)

[TABLE]

To finally evaluate $J^{(3)}_{0}$, we remark now that

[TABLE]

Then, as

[TABLE]

we have

[TABLE]

Appendix E Calculation of $\Delta F_{r}^{T}$

To evaluate explicitly $\Delta F_{r}^{T}$, let us start from Eq. (64),

[TABLE]

and, in order to perform the sum over $s$, we introduce integral representations

[TABLE]

Observing now that the value $k=1$ can be included in the sum over $k$ and defining

[TABLE]

we can write

[TABLE]

This implies that, for $h\to\infty$,

[TABLE]

which is exactly Eq. (65).