A General Statistical Physics Framework for Assignment Problems

Abstract

Linear assignment problems hold a pivotal role in combinatorial optimization, offering a broad spectrum of applications within the field of data sciences. They consist of assigning “agents” to “tasks” in a way that leads to a minimum total cost associated with the assignment. The assignment is balanced when the number of agents equals the number of tasks, with a one-to-one correspondence between agents and tasks, and it is and unbalanced otherwise. Additional options and constraints may be imposed, such as allowing agents to perform multiple tasks or allowing tasks to be performed by multiple agents. In this paper, we propose a novel framework that can solve all these assignment problems employing methodologies derived from the field of statistical physics. We describe this formalism in detail and validate all its assertions. A major part of this framework is the definition of a concave effective free energy function that encapsulates the constraints of the assignment problem within a finite temperature context. We demonstrate that this free energy monotonically decreases as a function of a parameter $\beta$ representing the inverse of temperature. As $\beta$ increases, the free energy converges to the optimal assignment cost. Furthermore, we demonstrate that when $\beta$ values are sufficiently large, the exact solution to the assignment problem can be derived by rounding off the elements of the computed assignment matrix to the nearest integer. We describe a computer implementation of our framework and illustrate its application to multi-task assignment problems for which the Hungarian algorithm is not applicable.

1. Introduction

An assignment problem can be interpreted as a problem of resource allocation in which a certain number of “tasks” are performed by some “agents” so that the total cost of pairing tasks with agents is minimized. There are many ways to constrain this problem, such as imposing one-to-one mapping (the balanced problem), searching for a specific number of assignments (the k-cardinality assignment problem), imposing constraints on the resources (the generalized assignment problem), or allowing agents to perform multiple tasks or reversely having some tasks performed by more than one agent. For reviews on the variety of assignment problems and the algorithms that have been proposed to solve them, we refer the reader to Refs. [1,2,3] and references therein. Assignment problems are fundamental combinatorial optimization problems. As such, solving such problems is a significant concern in operational research, economics, and data sciences, among other disciplines. Assignment problems have been and continue to be a focal point of investigation for mathematicians, statisticians, computer scientists, and even physicists. In this paper, we propose a general framework that regroups all the assignment problems mentioned above. We develop an approximate solution to this problem based on statistical physics and prove that it converges efficiently to an exact solution for non-degenerate as well as degenerate assignment problems.

Let us first define the assignment problem. Let T be the set of tasks and A the set of agents, with possibly different cardinalities $\left|T\right|=N_2$ and $\left|A\right|=N_1$. If we define $C(i,j)$ as the cost of assigning task i to the agent j, the balanced assignment problem can be stated as finding a transportation matrix G between T and A whose elements are non-negative integer numbers and that minimizes

(1)$U=\sum _{i\in T}\sum _{j\in A}C(i,j)G(i,j).$

In the special case of $N_1=N_2$, with a one-to-one assignment of tasks to agents (the balanced assignment problem), G corresponds to a permutation of $\{1,\dots ,N\}$. In the more general case, $N_1$ and $N_2$ can be different (unbalanced assignment), the number of pairings k may be imposed with $0<k\le min(N_1,N_2)$ (k-cardinality assignment problem), and/or tasks and agents may be assigned to more than one agents and tasks, respectively. In those general cases, G remains binary but may not be a permutation matrix (in the specific case of the k-cardinality problem, for example, G is a partial permutation matrix of rank k).

As mentioned above, assignment problems are basically combinatorial optimization problems that can be solved numerically using a linear programming model or, more specifically, a linear integer programming method. Those methods, however, are theoretically $NP$-complete. There are fast heuristics that can solve such problems, such as Dantzig’s simplex method [4]. The assignment problem, however, is a specialized case of integer linear programming that can be recast and solved by an exact polynomial algorithm. For example, the most common algorithm for solving the balanced assignment problem has its origin in the work of Jacobi [5]. It was “rediscovered” 60 years later by [6] and is now dubbed the Hungarian algorithm. It is a global algorithm that iteratively identifies assignments between agents and tasks. It is polynomial in running time ($O\left(N^4\right)$ or $O\left(N^3\right)$ depending on the implementation). While it is serial by nature, there are many recent efforts to parallelize this algorithm (see, for example [7,8,9], and references therein).

The Hungarian algorithm can be adapted to many unbalanced assignment problems. If there are $N_2$ tasks and $N_1$ agents, with $N_2>N_1$, for example, it is possible to add $N_2-N_1$ “dummy” agents, solve the balanced assignment problem of size $N_2$, and then only retain the assignments to the actual $N_1$ agents. If each agent can perform up to M tasks, the textbook solution is then to create M copies of each agent prior to using the Hungarian algorithm. This approach, however, may not lead to the expected solution. Consider a problem with fiv tasks and three agents, with each agent allowed to perform up to three tasks. The option of creating three copies of each agent and then solving the corresponding unbalanced assignment problem may lead to one agent performing two tasks, a second agent performing three tasks, and the last agent not performing any task, possibly violating a constraint that each agent needs to perform at least one task. Several methods have been proposed to circumvent this problem, either by dividing the problem into multiple subproblems [10,11] or by using random algorithms such as the ant colony algorithms [12,13]. It is unclear whether these algorithms can identify the optimal solution (they have been shown to fail in some cases, see for example [14]) or how they scale for large problems.

In this paper, we propose a radically different approach to solving a general class of assignment problems using continuous systems. Our approach is motivated by statistical physics. It is a full generalization of the method developed in a previous paper to solve balanced assignment problems [15]. In the following, we will use the generic term “multi-assignment” to represent an assignment problem with agents possibly performing multiple tasks, or with tasks assigned to multiple agents. In this paper, we aim to

• Introduce and substantiate a continuous framework for addressing potentially multifaceted assignment problems involving multiple tasks and/or multiple agents, leveraging principles from statistical physics;

• Demonstrate that, under generic circumstances where the multi-assignment problem possesses a unique solution, the aforementioned framework ensures convergence to that solution with arbitrarily high accuracy, in terms of both cost (referred to as energy) and assignment matrix,

• Present a modified approach capable of identifying at least one solution for degenerate multi-assignment problems that feature multiple solutions;

• Demonstrate that the implementation of this framework can be efficiently parallelized.

We emphasize that this formalism is not a mere adaptation but a full generalization of the framework we develop for solving assignment problems [16]. In particular, new features of this paper are as follows:

• A method to account for the fact that each task can be assigned to a variable number of agents and, on the other hand, that each agent can be assigned a variable number of tasks. To that end, we introduce indicator functions over the sets of tasks and agents that are optimized along with the transportation plan.

• The establishment of proofs of validity and convergence of our algorithm for those general multi-assignment problems. The main results are provided in the text, while the proofs themselves are relegated to the Appendices for more clarity. For the balanced assignment problem, these proofs rely heavily on the fact that the corresponding transportation matrices are permutation matrices that are the extreme points of the well-characterized convex set of doubly stochastic matrices (according to the Birkhoff–von Neumann theorem [17,18]). For more general multi-assignment problems, the transportation matrices also belong to a convex set, with properties akin to the Birkhoff–von Neumann theorem (these properties are discussed in Appendix A). The use of these properties to derive the convergence of our algorithm for solving multi-assignment problems is new.

The paper is organized as follows. In Section 2 Section 3 and Section 4, we describe in detail the general framework we propose for solving multi-assignment problems. Proofs of all important properties of this framework are provided in the appendices. Section 5 briefly describes the implementation of the method in a C++ program, Matching V1.0. In Section 6, we present some applications, along with a comparison to current algorithms.

2. A General Formulation of the Assignment Problem

We consider two sets of points, $S_1$ and $S_2$, with cardinalities $N_1$ and $N_2$, respectively. We represent the cost of transportation between $S_1$ and $S_2$ as a matrix C whose elements $C(i,j)$ are positive, for all $(i,j)\in [1,N_1]\times [1,N_2]$. We set the number of assignments between points in $S_1$ and $S_2$ to be a strictly positive integer number k, a constant given as input to the problem. A point i is assigned to $n_1\left(i\right)$ points in $S_2$, with $n_1\left(i\right)$ belonging to $[p_{min}\left(i\right),\dots ,p_{max}\left(i\right)]$, where $p_{min}\left(i\right)$ and $p_{max}\left(i\right)$ are non-negative integers satisfying $0\le p_{min}\left(i\right)\le p_{max}\left(i\right)\le N_2$. Similarly, a point j is assigned to $n_2\left(j\right)$ points in $S_1$, with $n_2\left(j\right)$ belonging to $[q_{min}\left(j\right),\dots ,q_{max}\left(j\right)]$, where $q_{min}\left(j\right)$ and $q_{max}\left(j\right)$ are non-negative integers integers satisfying $0\le q_{min}\left(j\right)\le q_{max}\left(j\right)\le N_1$. Using the traditional “task-agent” for the formulation of an assignment problem, $S_1$ represents the tasks, and $S_2$ the agents. $n_1\left(i\right)$ is the number of agents needed to perform task i, while $n_2\left(j\right)$ is the number of tasks that can be performed by agent j. $p_{min}\left(i\right)$ and $p_{max}\left(i\right)$ represent the minimum and maximum number of agents that are assigned to a task i, respectively, while $q_{min}\left(j\right)$ and $q_{max}\left(j\right)$ represent the minimum and maximum number of tasks that can be assigned to an agent j, respectively. We use a general definition of these boundaries, in that we allow each task and each agent to have their own limits. The multi-assignment problem is then framed as the search for an assignment matrix G that defines the correspondence between points in $S_1$ and points in $S_2$. This matrix is found by minimizing the assignment cost U, defined as

(2)$U(G,C)=\sum _{i,j}G(i,j)C(i,j).$

The summations encompass all i in $S_1$ and j in $S_2$. The objective is to locate the minimum of U given the values of $G(i,j)$ that adhere to the following constraints:

(3)$\begin{array}{c}\forall i,{\displaystyle \sum _j}G(i,j)=n_1\left(i\right)\hfill \\ \forall j,{\displaystyle \sum _i}G(i,j)=n_2\left(j\right)\hfill \\ {\displaystyle \sum _i\sum _j}G(i,j)=k\hfill \\ \forall (i,j),G(i,j)\in \{0,1\}\hfill \\ \forall \left(i\right),n_1\left(i\right)\in [p_{min}\left(i\right),\dots ,p_{max}\left(i\right)]\hfill \\ \forall \left(j\right),n_2\left(j\right)\in [q_{min}\left(j\right),\dots ,q_{max}\left(j\right)],\hfill \end{array}$

where k is the actual number of task–agent pairs. The numbers $p_{min}\left(i\right)$, $p_{max}\left(i\right)$, $q_{min}\left(j\right)$, $q_{max}\left(j\right)$, and k are problem-specific and given. They do satisfy some constraints, such as $0\le p_{min}\left(i\right)\le p_{max}\left(i\right)\le N_2$ and $0\le q_{min}\left(j\right)\le q_{max}\left(j\right)\le N_1$, as described above, as well as $0<k\le {\sum }_{i=1}^{N_1}{p}_{min}\left(i\right)$ and $k\le {\sum }_{j=1}^{N_2}{q}_{min}\left(j\right)$. We store those numbers in four vectors, $P_{min}=(p_{min}\left(1\right),\dots ,p_{min}\left(N_1\right))$, $P_{max}=(p_{max}\left(1\right),\dots ,p_{max}\left(N_1\right))$, $Q_{min}=(q_{min}\left(1\right),\dots ,q_{min}\left(N_2\right))$, $Q_{max}=(q_{max}\left(1\right),\dots ,q_{max}\left(N_2\right))$.

Equation (3) recap most of the standard linear assignment problems:

• (i). If $k=N_1=N_2$ and $p_{min}\left(i\right)=p_{max}\left(i\right)=q_{min}\left(j\right)=q_{max}\left(j\right)=1$ for all i and j, we recover the balanced assignment problem,

• (ii). If $p_{min}\left(i\right)=q_{min}\left(j\right)=0$ and $p_{max}\left(i\right)=q_{max}\left(j\right)=1$ for all i and j, the equations correspond to the k-cardinality assignment problem [19,20,21].

In Equation (3), $G(i,j)$, $n_1\left(i\right)$ and $n_2\left(j\right)$ are unknown, defining a total of $N_1{N}_2+N_1+N_2$ variables. The solution to the general assignment problem is given by the functions $n_1^{*}$ and $n_2^{*}$ on $S_1$ and $S_2$, respectively, which identify which tasks and which agents are involved in the optimal assignment, the transportation matrix $G^{*}$ that defines these correspondences, and the minimum assignment cost $U^{*}=U(G^{*},C)$.

Optimizing Equation (2) subject to the constraints Equation (3) constitutes a discrete optimization challenge, specifically an integer linear programming problem. To address this challenge, we employ a statistical physics methodology, transforming it into a temperature-dependent problem involving real variables. The integer optimal solution is then obtained as the temperature approaches zero.

2.1. Effective Free Energy for the General Assignment Problem

Solving the multi-assignment problem entails determining the minimum of the function defined in Equation (2) across the potential mappings between the two discrete sets of points under consideration. Identifying this minimum is tantamount to identifying the most probable state of the system it characterizes. This “system” encompasses the various binary transportation plans, also referred to as assignment matrices, between $S_1$ and $S_2$ that adhere to the constraints Equation (3). Each state within this system corresponds to assignment matrix G with its associated energy $U(G,C)$ as defined in Equation (2). The probability $P\left(G\right)$ linked with an assignment matrix G is given by

(4)$P\left(G\right)=\frac{1}{Z\left(\beta \right)}{e}^{-\beta U(G,C)}.$

In this equation, $\beta$ is the inverse temperature, namely $\beta =1/k_{B}T$, where $k_B$ is the Boltzmann constant and T is the temperature, and $Z\left(\beta \right)$ is the partition function computed over all states of the system. As this system is defined by the assignment matrices G and by the marginals of G, $n_1$, and $n_2$, the partition function is equal to

(5)$\begin{array}{ccc}\hfill Z_{\beta }& =& e^{-\beta {\mathcal{F}}_{\beta }}\hfill \\ & =& \sum _{G(i,j)=0}^1\sum _{n_1\left(i\right)=p_{min}\left(i\right)}^{p_{max}\left(i\right)}\sum _{n_2\left(j\right)=q_{min}\left(j\right)}^{q_{max}\left(j\right)}{e}^{-\beta U(G,C)},\hfill \end{array}$

where $\mathcal{F}\left(\beta \right)$ is the free energy of the system. The practicality of this free energy is limited due to the inability to compute it explicitly. We develop a method for approximating it through the saddle point approximation technique.

Taking into account the constraints in Equation (3), the partition function can be written as

(6)$\begin{array}{ccc}& & Z_{\beta }=\sum _{G(i,j)=0}^1\sum _{n_1\left(i\right)=p_{min}\left(i\right)}^{p_{max}\left(i\right)}\sum _{n_2\left(j\right)=q_{min}\left(j\right)}^{q_{max}\left(j\right)}{e}^{-\beta \sum _{i,j}C(i,j)G(i,j)}\hfill \\ & & \prod _i\delta \left(\sum _{j}G(i,j)-n_1\left(i\right)\right)\prod _j\delta \left(\sum _{i}G(i,j)-n_2\left(j\right)\right)\hfill \\ & & \delta \left(\sum _{i,j}G(i,j)-k\right).\hfill \end{array}$

The three sums impose that $G(i,j)$, $n_1\left(i\right)$, and $n_2\left(j\right)$ take integer values within some ranges defined by the constraints Equation (3). The additional constraints are imposed through the delta functions. We employ the Fourier representation of those delta functions, thereby introducing additional auxiliary variables x, $\lambda \left(i\right)$, and $\mu \left(j\right)$, with $i\in [1,N_1]$ and $j\in [1,N_2]$. The partition function is then given, after reorganization, by (up to a multiplicative constant) by

(7)$\begin{array}{ccc}\hfill Z_{\beta }=& & {\int }_{-\infty }^{+\infty }\prod _{i}d\lambda \left(i\right){\int }_{-\infty }^{+\infty }\prod _{j}d\mu \left(j\right){\int }_{-\infty }^{+\infty }dx{e}^{ß\beta kx}\hfill \\ & & \sum _{n_1\left(i\right)=p_{min}\left(i\right)}^{p_{max}\left(i\right)}\sum _{n_2\left(j\right)=q_{min}\left(j\right)}^{q_{max}\left(j\right)}{e}^{\beta \left(\sum _iß\lambda \left(i\right)n_1\left(i\right)+\sum _jß\mu \left(j\right)n_2\left(j\right)\right)}\hfill \\ & & \sum _{G(i,j)=0}^1{e}^{-\beta \sum _{i,j}G(i,j)(C(i,j)+ß\lambda \left(i\right)+ß\mu \left(j\right)+ßx)},\hfill \end{array}$

where $ß$ is the imaginary unit (${ß}^2=-1$). Note that we have rescaled the variables x, $\mathit{\lambda }$, and $\mathit{\mu }$ by a factor $\beta$ for consistency with the energy term. Conducting the summations over the variables $G(i,j)$, $n_1\left(i\right)$, and $n_2\left(j\right)$, we get

$Z_{\beta }={\int }_{-\infty }^{+\infty }\prod _{i}d\lambda \left(i\right){\int }_{-\infty }^{+\infty }\prod _{j}d{\mu }_j{\int }_{-\infty }^{+\infty }dx{e}^{-\beta F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)},$

where $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ is a functional, or effective free energy, that depends on the variables $\mathit{\lambda }$, $\mathit{\mu }$, and x and is defined by

(8)$\begin{array}{ccc}& & F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=\hfill \\ & & -\frac{1}{\beta }\sum _{i}ln\left(\frac{e^{ß\beta p_{min}\left(i\right)\lambda \left(i\right)}-e^{ß\beta (p_{max}\left(i\right)+1)\lambda \left(i\right)}}{1-e^{ß\beta \lambda \left(i\right)}}\right)\hfill \\ & & -\frac{1}{\beta }\sum _{j}ln\left(\frac{e^{ß\beta q_{min}\left(j\right)\mu \left(j\right)}-e^{ß\beta (q_{max}\left(j\right)+1)\mu \left(j\right)}}{1-e^{ß\beta \mu \left(j\right)}}\right)\hfill \\ & & -\frac{1}{\beta }\sum _{i,j}ln\left[1+e^{-\beta \left(C(i,j)+ß\lambda \left(i\right)+ß\mu \left(l\right)+ßx\right)}\right]-ßkx.\hfill \end{array}$

In contrast to the expression of the internal energy U defined in Equation (2) that depends on the $N_1{N}_2+N_1+N_2$ constrained binary variables $G(i,j)$, $n_1\left(i\right)$, and $n_2\left(j\right)$, the expression for the effective free energy $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ depends only on $N_1+N_2+1$ unconstrained variables, namely $\lambda \left(i\right)$, $\mu \left(j\right)$, and x. Below, we demonstrate how identifying the extremum of this function enables us to address the multi-assignment problem.

2.2. Optimizing the Effective Free Energy

Let $\overline{G}(i,j)$, $\overline{n_1}\left(i\right)$, and $\overline{n_2}\left(j\right)$ represent the expected values of $G(i,j)$, $n_1\left(i\right)$, and $n_2\left(j\right)$, respectively, in accordance with the Gibbs distribution specified in Equation (4). Computing these expected values directly is unfortunately not feasible as the partition function defined in Equation (8) lacks an analytical expression. Instead, we derive a saddle point approximation (SPA) by seeking extrema of the effective free energy with respect to the variables $\lambda$, $\mu$, and x:

(9)$\begin{array}{ccc}& & \frac{\partial F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial {\lambda }_i}=\frac{\partial F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial {\mu }_j}=\frac{\partial F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial x}=0.\hfill \end{array}$

After some rearrangements, these two equations can be written as

(10a)$\begin{array}{cc}\hfill & \forall i,\sum _{j}X(i,j)=d_1\left(i\right),\hfill \end{array}$

(10b)$\begin{array}{cc}\hfill & \forall j,\sum _{i}X(i,j)=d_2\left(j\right),\hfill \end{array}$

(10c)$\begin{array}{cc}\hfill & \sum _{i,j}X(i,j)=k,\hfill \end{array}$

where

(11)$\begin{array}{ccc}\hfill X(i,j)& =& h\left[\beta \left(C(k,l)+ß\lambda \left(i\right)+ß\mu \left(j\right)+ßx\right),\right]\hfill \\ \hfill d_1\left(i\right)& =& g\left(ß\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right)\right),\hfill \\ \hfill d_2\left(j\right)& =& g\left(ß\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right)\right),\hfill \end{array}$

and

$\begin{array}{ccc}\hfill h\left(x\right)=& & \frac{1}{e^x+1},\hfill \\ \hfill g(x,a,b)=& & \frac{(b-a+1)}{e^{(b-a+1)x}-1}+\frac{1}{1-e^x}+b.\hfill \end{array}$

Note that

• (i). $g(x,a,a)=a$. If $a=1$ and $N_1=N_2=1$, we recover Equation (11) from [16] corresponding to the balanced assignment problem.

• (ii). $g(x,0,1)=h(-x)$. In this case, Equation (11) refers to the k-cardinality problem.

As frequently encountered, the saddle point may be purely imaginary. In this instance, it is evident from Equations (10) and (11) that the variables $ß\lambda \left(i\right)$, $ß\mu \left(j\right)$, and $ßx$ must be real. Consequently, we will replace $\left\{ß\lambda \right(i),ß\mu (j),ßx\}$ by $\left\{\lambda \right(i),\mu (j),x\}$ below.

In order to analyze the saddle point approximation, SPA, it is necessary to verify the existence and evaluate the uniqueness of the extrema of the effective free energy. The following theorem proves that $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ is strictly concave.

Theorem 1.

The Hessian of the effective free energy $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ is negative definite if $p_{min}\left(i\right)<p_{max}\left(i\right)$ for all i and $q_{min}\left(j\right)<q_{max}\left(j\right)$, for all j, and negative semidefinite otherwise. The free energy function $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ is then strictly concave in the first case and concave in the second case.

Proof. 

See Appendix B.    □

The following property establishes a relationship between the solutions of the SPA system of equations and the expected values for the assignment matrix and for the indicator functions:

Proposition 1.

Let ${\overline{S}}_{\beta }$ be the expected state of the system at the temperature β with respect to the Gibbs distribution given in Equation (4). ${\overline{S}}_{\beta }$ is associated with an expected transportation plan ${\overline{G}}_{\beta }$ and expected indicator functions ${\overline{n_1}}_{\beta }$ and ${\overline{n_2}}_{\beta }$. Let ${\lambda }^{MF}\left(i\right)$, ${\mu }^{MF}\left(j\right)$, and $x^{MF}$ be the solutions of the system of Equation (10). Then the following identities hold:

(12) $\begin{array}{ccc}\hfill {\overline{G}}_{\beta }(i,j)& =& h\left(\beta (C(i,j)+{\lambda }^{MF}\left(i\right)+{\mu }^{MF}\left(j\right)+x^{MF})\right)\hfill \\ & =& X^{MF}(k,l)\hfill \\ \hfill {\overline{n}}_{1\beta }\left(i\right)& =& g(\beta {\lambda }^{MF}\left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))=d_1^{MF}\left(i\right)\hfill \\ \hfill {\overline{n}}_{2\beta }\left(j\right)& =& g(\beta {\mu }^{MF}\left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))=d_2^{MF}\left(j\right).\hfill \end{array}$

To indicate that the solutions are mean field solutions, we use the superscript $MF$.

Proof. 

The proof of this proposition is virtually identical to the proof of the same results for the assignment problem, found in Appendix B of [15].    □

For a given value of the parameter $\beta$, ${\overline{n}}_{1\beta }\left(i\right)$ and ${\overline{n}}_{2\beta }\left(j\right)$ are the numbers of elements of $S_2$ and the number of elements in $S_1$ that are in correspondence with i and j, respectively, and ${\overline{G}}_{\beta }$ forms a transportation plan between $S_1$ and $S_2$ that is optimal with respect to the free energy defined in Equation (8). Note that these values are mean values and fractional. The matrix ${\overline{G}}_{\beta }$ belongs to a polytope, which we define as ${\mathcal{U}}^k(P_{min}^{max},Q_{min}^{max})$ (see Appendix A for a full characterization of this polytope). We can link the mean field assignment matrix ${\overline{G}}_{\beta }$ to an optimal free energy $F_{\beta }^{MF}$ and to an optimal internal energy $U_{\beta }^{MF}={\sum }_{i,j}{\overline{G}}_{\beta }(i,j)C(i,j)$. These values serve as mean field approximations of the exact free energy and internal energy of the system, respectively. Here are the key properties of $U_{\beta }^{MF}$ and $F_{\beta }^{MF}$:

Proposition 2.

$F_{\beta }^{MF}$ and $U_{\beta }^{MF}$ are monotonic increasing and monotonic decreasing functions, respectively, of the parameter β.

Proof. 

See Appendix C for $F_{\beta }^{MF}$ and Appendix D for $U_{\beta }^{MF}$.    □

Theorem 1 and Proposition 2 outlined above underscore several advantages of the proposed framework, which reframes the multi-assignment problem as a temperature-dependent process. Firstly, at each temperature, the multi-assignment problem with a generic cost matrix is transformed into a concave problem with a unique solution. This problem exhibits linear complexity in the number of variables, contrasting with the quadratic complexity of the original problem. The concavity facilitates the utilization of straightforward algorithms for identifying a minimum of the effective free energy function (Equation (8)). Additionally, Equation (12) offers robust numerical stability for computing the transportation plan and the functions $n_1$ and $n_2$, owing to the characteristics of the functions $h\left(x\right)$ and $g\left(x\right)$. Lastly, the convergence with respect to temperature is monotonic.

3. Solving Generic Assignment Problems

In the previous section, we proposed an effective free energy, $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$, that depends on $N_1+N_2+1$ unconstrained variables and on the inverse of the temperature, $\beta$, whose extremum identifies the solution of a possibly multi-job, multi-agent assignment problem. We have shown that the trajectory of this extremum is monotonic, increasing with respect to $\beta$. We now relate these values to the optimal assignment energy $U^{*}$:

Theorem 2.

Let $F_{\beta }^{MF}$ and $U_{\beta }^{MF}$ be the mean field approximations of the exact free energy and internal energy of the system at the inverse temperature β. The optimal assignment energy $U^{*}$ is then given by,

(13) $\begin{array}{ccc}\hfill U^{*}& =& \underset{\beta \to +\infty }{lim}{F}_{\beta }^{MF},\hfill \\ \hfill U^{*}& =& \underset{\beta \to +\infty }{lim}{U}_{\beta }^{MF},\hfill \end{array}$

namely, the optimal assignment energy, $U^{*}$, is equal to the infinite limit with respect to the inverse temperature of both the mean field free energy and the internal energy.

Proof. 

See Appendix E.    □

This theorem illustrates that as the inverse temperature approaches infinity (or equivalently, as the temperature tends towards zero), the internal energy and free energy of the system converge to the optimal assignment energy. This confirms the validity of our statistical physics approach, particularly emphasizing the effectiveness of the saddle-point approximation. Note, however, that it does not define the behavior of the coupling matrix ${\overline{G}}_{\beta }=X^{MF}$. As ${\overline{G}}_{\beta }(i,j)=h\left(\beta (C(i,j)+{\lambda }^{MF}\left(i\right)+{\mu }^{MF}\left(j\right)+x^{MF})\right)$ and $0<h\left(x\right)<1$, the coupling matrix at any finite temperature is fractional. We need to show that as $\beta \to +\infty$, the corresponding matrix ${\overline{G}}_{\infty }$ does converge to the solution matrix $G^{*}$ and not to a fractional matrix that would lead to the same low energy $U^{*}$.

We first establish bounds on the internal energy and free energy at the SPA. Let us define

(14)$\begin{array}{ccc}\hfill A(N_1,N_2)=& & N_1{N}_{2}ln\left(2\right)+\sum _{i=1}^{N_1}ln(p_{max}\left(i\right)-p_{min}\left(i\right)+1)\hfill \\ & & +\sum _{j=1}^{N_2}ln(q_{max}\left(j\right)-q_{min}\left(j\right)+1).\hfill \end{array}$

Theorem 3.

$F_{\beta }^{MF}$ and $U_{\beta }^{MF}$, the mean field approximations of the exact free energy and internal energy of the system at the inverse temperature β, satisfy the following inequalities:

(15) $\begin{array}{ccc}\hfill U^{*}-\frac{A(N_1,N_2)}{\beta }& \le F_{\beta }^{MF}\le & U^{*},\hfill \end{array}$

(16) $\begin{array}{ccc}\hfill U^{*}& \le U_{\beta }^{MF}\le & U^{*}+\frac{A(N_1,N_2)}{\beta },\hfill \end{array}$

where $U^{*}$ is the optimal assignment energy and $A(N_1,N_2)$ is defined in Equation (14).

Proof. 

See Appendix F.    □

Now, let us assume that the multi-assignment problem linked with the $N_1\times N_2$ cost matrix C possesses a unique optimal assignment matrix, denoted as $G^{*}$. We have the following theorem:

Theorem 4.

Let ${\overline{G}}_{\beta }$ be the coupling matrix solution of the SPA equations at the inverse temperature β. We assume that the multi-assignment problem has a unique solution, with $G^{*}$ being the associated coupling matrix. If Δ is the difference in total cost between the optimal solution and the second-best solution, then

(17) $\begin{array}{c}\hfill \underset{i,j}{max}\left|{\overline{G}}_{\beta }(i,j)-G^{*}(i,j)\right|\le \frac{A(N_1,N_2)}{\beta \Delta }.\end{array}$

Proof. 

See Appendix G.    □

This theorem validates that, in the generic case in which the solution to the assignment problem is unique, the converged solution matrix ${\overline{G}}_{\infty }$ is this unique solution.

4. Solving Degenerate Assignment Problems

The method we propose is basically a relaxation approach to the general assignment problem. Indeed, we build a collection of transportation matrices ${\overline{G}}_{\beta }$ that belong to ${\mathcal{U}}^k(P_{min}^{max},Q_{min}^{max})\setminus {\mathcal{P}}^k(P_{min}^{max},Q_{min}^{max})$ (see Appendix A). Entries of these matrices are non-integer, strictly in the interval $(0,1)$. If the general assignment problem is known to have a unique integer solution, we have shown that these matrices converge to an element of ${\mathcal{P}}^k(P_{min}^{max},Q_{min}^{max})$ when $\beta \to +\infty$. The question remains as to what happens when the problem is degenerate, i.e., when it may have multiple integer solutions.

The general assignment problem is a linear programming problem. Checking if such a problem is degenerate is unfortunately often $NP$-complete ([22,23]). The degeneracies occur due to the presence of cycles in the linear constraints, i.e., in the cost matrix for the assignment problem. If this is the case, we propose randomly perturbing that matrix to bring it back to the generic problem. Megiddo and colleagues [24] have shown that an $\epsilon$-perturbation of a degenerate linear programming problem reduces this problem to a non-degenerate one.

5. Matching: A Program for Solving Assignment Problems

We have implemented the multi-assignment framework described here in a C++ program, Matching, that is succinctly described in Algorithm 1.

Matching relies on an iterative process in which the parameter $\beta$ is gradually increased. At each value of $\beta$, the nonlinear system of equations defined by Equation (10) is solved. We write this system as

$\begin{array}{c}\hfill \left\{\begin{array}{c}{\mathbf{A}}_{\lambda }=0,\hfill \\ {\mathbf{A}}_{\mu }=0,\hfill \\ A_x=0,\hfill \end{array}\right.\end{array}$

where $\mathbf{A}=({\mathbf{A}}_{\lambda },{\mathbf{A}}_{\mu },A_x)$ is a vector of predicates defined as

$\begin{array}{ccc}\hfill A_{\lambda }\left(i\right)=& & \sum _{j}h\left(\beta (C(i,j)+\lambda \left(i\right)+\mu \left(j\right)+x)\right)\hfill \\ & & -g(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right)),\hfill \\ \hfill A_{\mu }\left(j\right)=& & \sum _{i}h\left(\beta (C(i,j)+\lambda \left(i\right)+\mu \left(j\right)+x)\right)\hfill \\ & & -g(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right)),\hfill \\ \hfill A_x=& & \sum _{i,j}h\left(\beta (C(i,j)+\lambda \left(i\right)+\mu \left(j\right)+x)\right)-k.\hfill \end{array}$

This system has $N_1+N_2+1$ equations, with the same number of variables. It is solved using an iterative Newton–Raphson method (for details, see, for example [16,25]). Once the SPA system of equations is solved, the assignment matrix ${\overline{G}}_{\beta }$, the functions ${\overline{n}}_{1\beta }$, ${\overline{n}}_{2\beta }$, and the corresponding transportation energy $U^{MF}\left(\beta \right)$ are computed. The iterations over $\beta$ are stopped if the mean field energy $U_{\beta }^{MF}$ does not change anymore (within a tolerance TOL usually set to ${10}^{-6}$). Otherwise, $\beta$ is increased, and the current values of $\mathit{\lambda }$, $\mathit{\mu }$ and x are used as input for the following iteration. At convergence, the values of the assignment matrix are rounded to the nearest integer (indicated as in the output of Algorithm 1). The minimal energy is then computed using the corresponding integer matrix.

The primary computational expense of this algorithm arises from solving the nonlinear set of equations corresponding to the SPA at every $\beta$ value. We use for this purpose an iterative Newton–Raphson method (see [15]).

Algorithm 1 Matching: a temperature-dependent framework for solving the multi-assignment problem

Input:$N_1$ and $N_2$, the number of agents and tasks; $p_{min}\left(i\right)$ and $p_{max}\left(i\right)$, the minimum and maximum number of agents needed to perform task i, and $q_{min}\left(j\right)$ and $q_{max}\left(j\right)$, the minimum and maximum number of tasks that can be performed by an agent j; k, the expected number of assignments; the cost matrix C. Initial value ${\beta }_0$ for $\beta$, the inverse of the temperature Initialize: Initialize arrays $\mathit{\lambda }$ and $\mathit{\mu }$ to 0 and initialize $x=0$. Set $STEP=\sqrt{10}$. for$i=1,\dots$ until convergence do     (1). Initialize ${\beta }^i=STEP\ast {\beta }^{i-1}$.     (2). Solve non linear Equation (10) for $\mathit{\lambda }$, $\mathit{\mu }$ and x at saddle point     (3). Compute corresponding ${\overline{G}}_{\beta }^i$, ${\overline{n}}_{1\beta }^i$, ${\overline{n}}_{2\beta }^i$, and $U_{\beta }^{i,MF}$     (4). Check for convergence: if $|U_{\beta }^{i,MF}-U_{\beta }^{i-1,MF}|<TOL$, stop end for Output: The converged assignment matrix $⌊{\overline{G}}_{\beta }⌉$, the functions $⌊{\overline{n}}_{1\beta }⌉$, $⌊{\overline{n}}_{2\beta }⌉$ over the agents and tasks, and the minimal associated cost $U_{\beta }$.

As for any annealing scheme, the initial temperature, or, in our case, the initial value ${\beta }_0$, is a parameter that significantly impacts the efficiency of the algorithm. Setting ${\beta }_0$ to be too small (i.e., a high initial temperature) will lead to inefficiency as the system will spend a significant amount of time at high temperatures, while setting ${\beta }_0$ too high will require many steps to converge at the corresponding low temperature, thereby decreasing the efficiency brought about by annealing. The value of $\beta$ scales the cost matrix C and as such is related to the range of this matrix, more specifically to its largest value, $C_{max}$. We found that setting ${\beta }_0{C}_{max}=1$ provides satisfactory annealing efficiency. The value for $STEP$ was chosen empirically.

6. Examples

The framework we propose is general: it allows us to solve balanced and unbalanced assignment problems, including those that allow for multi-agent or multi-task assignments. We illustrate our method on the latter types of problems, as, to our knowledge, there are currently no solutions to those that are guaranteed to reach the optimal assignment.

As a first illustration of our framework, we ran Matching on the cost matrix $C_1$ given in Table 1. This matrix has been used in previous studies of multi-agent assignment problems [10,11,13,14,26,27]. Matching was run with the following parameters: $p_{min}=0$ for all agents (i.e., some agents may stay idle), $p_{max}=4$ for all agents (i.e., an agent can perform up to 4 tasks), $q_{min}=q_{max}=1$ (i.e., a task is only performed by one agent), and $k=8$ (all tasks are performed). In Figure 1, we show the corresponding trajectories of the internal energy $U_{\beta }^{MF}$ and free energy $F_{\beta }^{MF}$.

As expected, the internal energy is monotonically decreasing while the free energy is monotonically increasing, and both converge to the same value, 1440. Theorem 3 provides bounds on those energy values. Note that it can be rewritten as

$\begin{array}{c}\hfill U_{\beta }^{MF}-\frac{A(N_1,N_2)}{\beta }\le U^{*}\le U_{\beta }^{MF},\\ \hfill F_{\beta }^{MF}\le U^{*}\le U_{\beta }^{MF}+\frac{A(N_1,N_2)}{\beta },\end{array}$

i.e., at each value of $\beta$, we have bounds based on internal energy and free energy for the actual optimal cost of the assignment problem. These bounds are shown as shaded areas in Figure 1. Note that the widths of the intervals defined by the bounds are inversely proportional to $\beta$ and therefore decrease monotonically as $\beta$ increases.

In Table 2, we compare the assignment found by Matching on $C_1$ with $p_{max}$ set to 4 (i.e., up to four tasks per agent), and $p_{min}=0$ (agents can be idle) or $p_{min}=1$ (each agent must execute at least one task). The removal of the latter constraint leads to a better optimal assignment in which agent 3 remains idle. When we maintain this constraint, our results are the same as those found by Wang et al. [13].

As a second more general test of our method, we compared the results of the framework we propose with the results of our own version of the Hungarian algorithm, as well as with those found in previous studies. We used both the matrix C1 and the matrix C2 given in Table 3 in those computing experiments.

We ran Matching as follows. For all tasks j, we set $q_{min}\left(j\right)=q_{max}\left(j\right)=1$, i.e., a task is only performed by one agent, and all tasks are expected to be executed. The latter is further enforced by setting k to be the total number of tasks (8 for $C_1$ and 10 for $C_2$). For an agent i, we either set all $p_{min}\left(i\right)$ to be 1 (in which case all agents must perform at least one task) or 0 (in which case some agents may not be assigned). We considered all $p_{max}\left(i\right)$ to be equal to a given value p and ran Matching with p varying from 1 to k. Note that both matrices $C_1$ and $C_2$ contain cycles. As such, optimal assignments based on those matrices are not unique. To remove the degeneracy, we added random uniform noise between $[0,ϵ]$ with $ϵ=0.001$ to all values of the cost matrices (see Section 4 for details).

The Hungarian algorithm remains a standard for solving balanced as well as unbalanced assignment problems. It needs to be adapted, however, when solving multi-task assignment problems. The corresponding textbook solution that was used, for example, in Ref. [14] is to make copies or clones of the agents, solve the augmented unbalanced assignment problem using the Hungarian algorithm, and then regroup all clones of an agent to define the tasks that it is assigned to. We consequently created multiple versions of the cost matrices $C_1$ and $C_2$, with all agents copied p times, where p defines the maximum number of tasks that an agent can perform. We ran our own version of the Hungarian algorithm on those matrices. This version is adapted from the serial code written by C. Ma and is publicly available at https://github.com/mcximing/hungarian-algorithm-cpp, (accessed on 1 May 2021).

Comparisons of the results of Matching, of the Hungarian algorithm, and of previous studies are provided in Table 4 for the cost matrix C1 and in Table 5 for the cost matrix C2.

The results of Matching have been proven to be optimal (see Section 3). As expected, those results are strongly dependent on the constraints imposed on the problem: we can find assignments with a lower total cost if we allow agents to remain idle by setting $p_{min}$ to 0. For example, for the cost matrix $C_1$, in the extreme case with $p_{min}=0$ and $p_{max}=8$ for all agents, we find that all eight tasks have been assigned to agent 5, while the other agents have no task to perform, with a total cost of 1400. If instead we set $p_{min}=1$ and keep $p_{max}=8$, we find that agent 5 is assigned tasks 1, 2, 5, and 6; agent 1 is assigned task 3; agent 2 is assigned task 8; agent 3 is assigned task 4; and agent 4 is assigned task 7, for a total cost of 1450. If we instead run the Hungarian algorithm with agents that have been cloned to allow for multi-task assignments, we find optimal assignments that match those found by Matching with $p_{min}$ set to 0 for both cost matrices. For matrix $C_1$, for example, if each agent is represented with 8 clones, we find that agent 5 is assigned tasks 1, 2, 4, 5, 6, 7, and 8 while agent 1 is assigned task 3, for a total cost of 1400. This illustrates two points. First, as expected, the optimal assignment is not unique, as we find two different assignments with the same cost, 1400. Second, and more importantly, this shows the difficulty of applying the Hungarian algorithm for such problems. It works fine if we want to find the optimal assignment without consideration of constraints, such as the constraint each agent needs to perform at least one task, but does not fit if such a constraint is necessary.

The optimal solutions found by Matching when $p_{min}$ is set to 1 match those found by Wang et al. [13] for both matrices and are better than those found by Majumdar et al. [26] for matrix $C_2$. Note that the latter are obtained using either an ant colony algorithm or a genetic algorithm. Both are probabilistic algorithms that do not guarantee that the true minimum is found.

7. Conclusions

We have developed a general framework for solving balanced and unbalanced and constrained and unconstrained assignment problems. Given two sets of points $S_1$ and $S_2$ with the possibly different cardinalities $N_1$ and $N_2$, constraints on the number of possible assignments for each element of $S_1$ and of $S_2$, and a cost matrix defining the penalties of on assignment, we have constructed a concave free energy parameterized by temperature that captures those constraints. The definition of this free energy is general enough that it includes balanced and unbalanced cases. Its extremum establishes an optimal assignment between the two sets of points. We have demonstrated that this free energy consistently decreases as a function of $\beta$ (the inverse of temperature) towards the optimal assignment cost. Moreover, we have established that for sufficiently large $\beta$ values, the precise solution to the generic multi-assignment problem can be directly derived through straightforward rounding to the nearest integer of the elements of the assignment matrix linked to this extremum. Additionally, we have developed a method that guarantees convergence for addressing degenerate assignment problems.

The formalism introduced in this paper was designed to generalize the Hungarian algorithm for a large range of assignment problems. We expect it to be useful for a much larger set of problems, especially those associated with graphs. Graphs capture data relationships and, as such, are essential to numerous applications. They are particularly useful in domains such as Web search [28], neural [29] and social network analysis [30], gene networks [31], etc. More generally, they are designed to represent complex knowledge. The scale of modern graphs leads to a need to develop more efficient methods to process very large graphs. The formalism we propose is well adapted to tackle this problem for applications such as bipartite or simple graph matching. We will also consider extensions to higher-order matching problems, such as k-matching problems [32] (for example the three-assignment problem [33]) that are known to be $NP$-complete [34]. These problems have their own applications for graph analyses.

Author Contributions

Conceptualization, P.K. and H.O.; methodology, P.K. and H.O.; software, P.K.; validation, P.K.; formal analysis, P.K. and H.O.; investigation, P.K. and H.O.; writing—original draft preparation, P.K. and H.O.; writing—review and editing, P.K. and H.O. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

No new data were created for this study.

Acknowledgments

The work discussed here originated from visits by P.K. at the Institut de Physique Théorique, CEA Saclay, France, during the falls of 2022 and 2023. He thanks them for their hospitality and financial support.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figure and Tables

Figure 1. Convergence of the internal energy [Forumla omitted. See PDF.] (A) and free energy [Forumla omitted. See PDF.] (B) as a function of [Forumla omitted. See PDF.] when solving the multi-task assignment problem with the cost matrix [Forumla omitted. See PDF.], such that each agent can perform up to 4 tasks, may be idle, and all 8 tasks are performed. On both panels, we show the bounds on the expected value for the optimal cost [Forumla omitted. See PDF.] as shaded areas (see text for details).

Table 1

The cost matrix $C_1$.

	Tasks
Agents	${\mathit{T}}_{\mathbf{1}}$	${\mathit{T}}_{\mathbf{2}}$	${\mathit{T}}_{\mathbf{3}}$	${\mathit{T}}_{\mathbf{4}}$	${\mathit{T}}_{\mathbf{5}}$	${\mathit{T}}_{\mathbf{6}}$	${\mathit{T}}_{\mathbf{7}}$	${\mathit{T}}_{\mathbf{8}}$
$A_1$	300	250	180	320	270	190	220	260
$A_2$	290	310	190	180	210	200	300	190
$A_3$	280	290	300	190	190	220	230	260
$A_4$	290	300	190	240	250	190	180	210
$A_5$	210	200	180	170	160	140	160	180

Table 2

Assignments for $C_1$ with $p_{max}=4$.

	${\mathit{p}}_{\mathbf{min}}=0$	${\mathit{p}}_{\mathbf{min}}=1$	Wang a
	$A_1\to T_3$	$A_1\to T_3$	$A_1\to T_3$
	$A_2\to T_4,T_8$	$A_2\to T_8$	$A_2\to T_8$
	$A_3$	$A_3\to T_4$	$A_3\to T_4$
	$A_4\to T_7$	$A_4\to T_7$	$A_4\to T_7$
	$A_5\to T_1,T_2,T_5,T_6$	$A_5\to T_1,T_2,T_5,T_6$	$A_5\to T_1,T_2,T_5,T_6$
Total cost	1440	1450	1450

^a Assignment based on an ant colony algorithm [13].

Table 3

The cost matrix $C_2$.

	Tasks
Agents	${\mathit{T}}_{\mathbf{1}}$	${\mathit{T}}_{\mathbf{2}}$	${\mathit{T}}_{\mathbf{3}}$	${\mathit{T}}_{\mathbf{4}}$	${\mathit{T}}_{\mathbf{5}}$	${\mathit{T}}_{\mathbf{6}}$	${\mathit{T}}_{\mathbf{7}}$	${\mathit{T}}_{\mathbf{8}}$	${\mathit{T}}_{\mathbf{9}}$	${\mathit{T}}_{\mathbf{10}}$
$A_1$	10	2	14	9	6	7	21	32	18	11
$A_2$	7	12	9	3	5	6	9	16	54	12
$A_3$	4	8	6	12	21	9	21	14	45	13
$A_4$	21	9	12	9	32	10	19	25	16	10
$A_5$	10	12	30	15	12	17	30	12	12	9
$A_6$	15	7	34	17	7	16	14	17	9	5

Table 4

Solving the multi-task assignment problems defined by the cost matrix $C_1$.

Method $\setminus \mathit{p}$ a	2	3	4	5	6	7	8
Matching, $p_{min}=1$	1520	1470	1450	1450	1450	1450	1450
Matching, $p_{min}=0$	1520	1470	1440	1420	1410	1400	1400
Hungarian	1520	1470	1440	1420	1410	1400	1400
Majumdar GA b	1520	1470	1450	NA c	NA	NA	NA
Wang AC d	1520	1470	1450	NA	NA	NA	NA

^a Each agent is allowed to execute up to p tasks. Note that as there are 8 tasks and 5 agents only, the minimum value of p is 2 if all the tasks must be executed. ^b Genetic algorithm [26]. ^c Not available in either [26] or [13]. ^d Ant colony algorithm [13].

Table 5

Solving the multi-task assignment problems defined by the cost matrix $C_2$.

Method $\setminus \mathit{p}$ a	2	3	4	5	6	7	8
Matching, $p_{min}=1$	66	65	65	65	65	65	65
Matching, $p_{min}=0$	66	62	61	61	61	61	61
Hungarian	66	62	61	61	61	61	61
Majumdar GA b	66	67	66	74	NA c	NA	NA
Wang AC d	66	65	65	NA	NA	NA	NA

^a Each agent is allowed to execute up to p tasks. Note that as there are 10 tasks and 6 agents only, and the minimum value of p is 2 if all the tasks must be executed. ^b Genetic algorithm [26]. ^c Not available in either [26] or [13]. ^d Ant colony algorithm [13].

Appendix A. A Generalized Birkhoff–von Neumann Theorem

An $N\times N$ matrix $A=\left(a\right(i,j\left)\right)$ is said to be doubly stochastic if and only if it satisfies the following conditions, $$\begin{array}{ccc}\hfill a(i,j)& \ge & 0,\hfill \\ \hfill \sum _{i=1}^{N}a(i,j)& =& 1,\sum _{j=1}^{N}a(i,j)=1,\hfill \end{array}$$ for all $(i,j)\in {[1,N]}^2$. The set ${\Omega }_N$ of doubly stochastic matrices is a convex polytope whose vertices are the permutation matrices with the same size. This is expressed by the following theorem, established by [17,18]:

Theorem A1.

An $N\times N$ matrix A is doubly stochastic if and only if it can be written as a weighted sum of permutation matrices, i.e., $$\begin{array}{c}\hfill A=\sum _{\pi \in {\Pi }_N}{a}_{\pi }\pi ,\end{array}$$

where ${\Pi }_N$ is the set of permutation matrices of size N, $a_{\pi }$ is a positive real number and ${\displaystyle \sum _{\pi \in {\Pi }_N}}{a}_{\pi }=1$.

Doubly stochastic matrices and their properties associated with the Birkhoff–von Neumann theorem above proved useful to establish convergence properties of statistical physics frameworks for solving the balanced assignment problem [15,35].

For the more general assignment problem considered here, however, we need to consider a different but related set of matrices, defined as follows. Let N and M be two natural numbers. Let $R_{min}=(r\left(1\right),\dots ,r\left(N\right))$, $R_{max}=(R\left(1\right),\dots ,R\left(N\right))$, $C_{min}=(c\left(1\right),\dots ,c\left(M\right))$, and $C_{max}=(C\left(1\right),\dots ,C\left(M\right))$ be four non-negative integral vectors satisfying $$\begin{array}{cc}& \forall i,0\le r\left(i\right)\le R\left(i\right),\\ & \forall j,0\le c\left(j\right)\le C\left(j\right).\end{array}$$

Let A be a non-negative matrix of size $N\times M$, and let us denote the row sum vector of A as $RA$ and the column sum vector of A as $CA$. Let us also define $\sigma \left(A\right)$ as the sum of all elements of A, i.e., $\sigma \left(A\right)={\sum }_{i=1}^N{\sum }_{j=1}^{M}A(i,j)$. The transportation polytope $\mathcal{U}(R_{min}^{max},C_{min}^{max})$ is the set of $N\times M$ matrices A that satisfy $$\begin{array}{cc}& \forall (i,j),0\le A(i,j)\le 1,\\ & \forall i,r\left(i\right)\le RA\left(i\right)\le R\left(i\right),\\ & \forall j,c\left(j\right)\le CA\left(j\right)\le C\left(j\right).\end{array}$$

The transportation polytope ${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$ is the set of $N\times M$ matrices A that satisfy $$\begin{array}{c}\hfill {\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})=\left\{A\in \mathcal{U}(R_{min}^{max},C_{min}^{max})|\sigma \left(A\right)=k\right\}\end{array}$$

Denote $\mathcal{P}(R_{min}^{max},C_{min}^{max})$ as the set of all matrices in $\mathcal{U}(R_{min}^{max},C_{min}^{max})$ whose entries are either 0 or 1 with a similar definition for ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$ with respect to ${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$.

Relaxed solutions to the general assignment problem considered in this paper belong to ${\mathcal{U}}^k(P_{min}^{max},Q_{min}^{max})$, while integer solutions belong to ${\mathcal{P}}^k(P_{min}^{max},Q_{min}^{max})$, where $P_{min}$, represents the vectors containing the constraints on the number of tasks that can be assigned to an agent or the number of agents that can be assigned to a task. In Koehl [36], we show the following theorem:

Theorem A2.

${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$ is the convex hull of all matrices in ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$.

Theorem A2 can be stated as any matrix in ${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$ that can be written as a linear combination of matrices in ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$. This result will be useful for some of the proofs below.

Remark A1.

• (i) 

  If $R_{min}=R_{max}=1_N$, $C_{min}=C_{max}=1_M$ (where $1_N$ and $1_N$ are vectors of one of size N and M, respectively, $N=M=k$, $\mathcal{U}(R_{min}^{max},C_{min}^{max})$ is the set of doubly stochastic matrices, ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$ is the set of permutation matrices, and Theorem A2 is then equivalent to the Birkhoff–von Neumann theorem for doubly stochastic matrices [17,18].

• (ii) 

  If $R_{min}=0_n$, $C_{min}=0_m$, $R_{max}=1_n$, $C_{max}=1_m$, ${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$ is the set of doubly sub-stochastic matrices with sum k and ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$ is the set of sub-permutation matrices of rank k; a specific version of Theorem A2 was established by Mendelssohn and Dulmage for square matrices ([37]), and later by Brualdi and Lee for rectangular matrices ([38]).

Appendix B. Proof of Theorem 1: Concavity of the Effective Free Energy

The free energy associated with the multi-assignment problem can be written as $$\begin{array}{ccc}\hfill F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=& & -\frac{1}{\beta }\sum _{i}G(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))\hfill \\ & & -\frac{1}{\beta }\sum _{j}G(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))\hfill \\ & & +\frac{1}{\beta }\sum _{i,j}H\left(\beta \left(C\right(i,j)+\lambda (i)+\mu (l)+x)\right)-kx,\hfill \end{array}$$ with:$$\begin{array}{ccc}\hfill H\left(x\right)& =& -ln(1+e^{-x}),\hfill \\ \hfill G(x,a,b)& =& ln\left(\frac{e^{ax}-e^{(b+1)x}}{1-e^x}\right).\hfill \end{array}$$ where a and b are two positive integers with $0\le a\le b$. The functions $H\left(x\right)$, its derivatives, and some associated functions have been fully characterized in [15]. In parallel, the functions $G(x,a,b)$, its derivatives, and associated functions are characterized in Appendix H.

We prove that this free energy is weakly concave by showing that its Hessian $HE$ is negative-definite. $HE$ is a symmetric matrix of size $(N_1+N_2+1)\times (N_1+N_2+1)$, such that its rows and columns correspond to all $N_1$$\lambda$ values first, followed by all $N_2$$\mu$ values, and finally to the value x. Let $h\left(x\right)=H^{\prime }\left(x\right)$, and let $h^{\prime }$ be its derivative, i.e., (A1)$$h^{\prime }\left(x\right)=-\frac{e^x}{{(1+e^x)}^2}.$$

In [15], we have shown that $h^{\prime }\left(x\right)\in [\frac{-1}{4},0)\forall x\in \mathbb{R}$, i.e., that $h^{\prime }\left(x\right)$ is always strictly negative. Similarly, let $g\left(x\right)=G^{\prime }\left(x\right)$, and let $g^{\prime }$ be its derivative, i.e., (A2)$$g^{\prime }(x,a,b)=-\frac{{(b-a+1)}^2{e}^{(b-a+1)x}}{{(e^{(b-a+1)x}-1)}^2}+\frac{e^x}{{(e^x-1)}^2}$$

In Appendix H, we show that when $b=a$, $g^{\prime }(x,a,a)=0\forall x\in \mathbb{R}$, and, when $b>a$, $g^{\prime }(x,a,b)$ is strictly positive $\forall x\in \mathbb{R}$.

We define the matrix $X^{\prime }$ and the vector ${\mathbf{d}}_{\mathbf{1}}^{\prime }$ and ${\mathbf{d}}_{\mathbf{2}}^{\prime }$ such that $$\begin{array}{ccc}\hfill X^{\prime }(i,j)& =& h^{\prime }(\beta (C(i,j)+\lambda \left(i\right)+\mu \left(j\right)+x),\hfill \\ \hfill d_1^{\prime }\left(i\right)& =& g^{\prime }(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right)),\hfill \\ \hfill d_2^{\prime }\left(j\right)& =& g^{\prime }(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right)).\hfill \end{array}$$

From Equation (10), we obtain $$\begin{array}{ccc}\hfill HE(i,i^{\prime })& =& \frac{{\partial }^2{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial \lambda \left(i\right)\partial \lambda \left(i^{\prime }\right)}=\beta {\delta }_{i{i}^{\prime }}\left(\sum _j{X}^{\prime }(i,j)-d_1^{\prime }\left(i\right)\right),\hfill \\ \hfill HE(i,j)& =& \frac{{\partial }^2{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial \lambda \left(i\right)\partial \mu \left(j\right)}=\beta X^{\prime }(i,j),\hfill \\ \hfill HE(i,N)& =& \frac{{\partial }^2{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial \lambda \left(i\right)\partial x}=\beta \sum _j{X}^{\prime }(i,j),\hfill \\ \hfill HE(j,j^{\prime })& =& \frac{{\partial }^2{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu }\right)}{\partial \mu \left(j\right)\partial \mu \left(j^{\prime }\right)}=\beta {\delta }_{j{j}^{\prime }}\left(\sum _i{X}^{\prime }(i,j)-d_2^{\prime }\left(j\right)\right),\hfill \\ \hfill HE(j,N)& =& \frac{{\partial }^2{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial \mu \left(j\right)\partial x}=\beta \sum _i{X}^{\prime }(i,j),\hfill \\ \hfill HE(N,N)& =& \frac{{\partial }^2{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial x^2}=\beta \sum _i\sum _j{X}^{\prime }(i,j),\hfill \end{array}$$ where $\delta$ are Kronecker functions, the indices i and $i^{\prime }$ belong to $[1,N_1]$ and the indices j and $j^{\prime }$ belong to $[1,N_2]$, and we have defined $N=N_1+N_2+1$.

Let $\mathbf{x}=({\mathbf{x}}_{\mathbf{1}},{\mathbf{x}}_{\mathbf{2}},x_3)$ be an arbitrary vector of size N. The quadratic form $Q\left(\mathbf{x}\right)={\mathbf{x}}^{T}H\mathbf{x}$ is equal to $$\begin{array}{ccc}\hfill Q\left(\mathbf{x}\right)=& & \sum _{i,i^{\prime }}{x}_1\left(i\right)HE(i,i^{\prime })x_1\left(i^{\prime }\right)+\sum _{j,j^{\prime }}{x}_2\left(j\right)HE(j,j^{\prime })x_2\left(j^{\prime }\right)+\hfill \\ & & x_{3}HE(N,N)x_3+2\sum _{i,j}{x}_1\left(i\right)HE(i,j)x_2\left(j\right)+\hfill \\ & & 2\sum _i{x}_1\left(i\right)HE(i,N)x_3+2\sum _j{x}_2\left(j\right)HE(j,N)x_3\hfill \\ \hfill =& & \beta \sum _{i,j}{(x_1\left(i\right)+x_2\left(j\right)+x_3)}^2{X}^{\prime }(i,j)\hfill \\ & & -\beta \left(\sum _i{x}_1{\left(i\right)}^2{d}_1^{\prime }\left(i\right)+\sum _j{x}_2{\left(j\right)}^2{d}_2^{\prime }\left(j\right)\right).\hfill \end{array}$$

As $X^{\prime }(i,j)$ is based on the function $h^{\prime }$ that is negative and $d_1^{\prime }\left(i\right)$, and $d_2^{\prime }\left(j\right)$ are based on the function $g^{\prime }$ that is positive, the summands in the equation above are all negative for all $i\in [1,N_1]$ and $j\in [1,N_2]$, and therefore, $Q\left(\mathbf{x}\right)$ is negative for all vector $\mathbf{x}$. The Hessian H is negative, semidefinite.

To check for definitiveness, let us note first that as $Q\left(\mathbf{x}\right)$ is a sum of negative terms, it is 0 if and only if all the terms are equal to 0. As the function $h^{\prime }\left(x\right)$ is strictly negative, this means that $\forall (i,j)$, (A3)$$\begin{array}{ccc}\hfill {(x_1\left(i\right)+x_2\left(j\right)+x_3)}^2& =& 0.\hfill \end{array}$$

For the two other terms, we consider two cases:

• (i)

  $p_{min}\left(i\right)<p_{max}\left(i\right)\forall i$ and $q_{min}\left(j\right)<q_{max}\left(j\right)\forall j$. Then $g^{\prime }(x,p_{min}\left(i\right),p_{max}\left(i\right))<0$ and $g^{\prime }(x,q_{min}\left(j\right),q_{max}\left(j\right))<0$. This leads to $\forall (i,j)$: (A4)$$\begin{array}{ccc}\hfill x_1{\left(i\right)}^2& =& 0,\hfill \\ \hfill x_2{\left(j\right)}^2& =& 0,\hfill \end{array}$$

  Equations (A3) and (A4) are satisfied when all $x_1\left(i\right)$, $x_2\left(j\right)$, and when $x_3$ are zero, namely that $\mathbf{x}=\mathbf{0}$. Therefore in this case, H is negative, definite, and the free energy $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ is strictly concave.

• (ii)

  For all i$p_{min}\left(i\right)=p_{max}\left(i\right)$ or for all j$q_{min}\left(j\right)=q_{max}\left(j\right)$. Then, either $d_1^{\prime }\left(x\right)=0\forall x\in {\mathbb{R}}^{N_1}$ or $d_2^{\prime }\left(x\right)=0\forall x\in {\mathbb{R}}^{N_2}$, or both. There are then non-zero solutions to the equation $Q\left(\mathbf{x}\right)=0$. The Hessian is then only semidefinite.

Appendix C. Monotonicity of $F_{\beta }^{MF}$

The effective free energy $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ defined in Equation (8) is a function of the cost matrix C and of real unconstrained variables $\lambda \left(i\right)$, $\mu \left(j\right)$, and x. For the sake of simplicity, for any $(i,j)\in [1,N_1]\times [1,N_2]$, we define $$\begin{array}{ccc}\hfill y(i,j)& =& C(i,j)+\lambda \left(i\right)+\mu \left(j\right)+x,\hfill \\ \hfill y^{MF}(i,j)& =& C(i,j)+{\lambda }^{MF}\left(i\right)+{\mu }^{MF}\left(j\right)+x^{MF}.\hfill \end{array}$$

The effective free energy is then $$\begin{array}{ccc}\hfill F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=& & -\frac{1}{\beta }\sum _{i}G(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))\hfill \\ & & -\frac{1}{\beta }\sum _{j}G(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))\hfill \\ & & +\frac{1}{\beta }\sum _{i,j}H\left(\beta y(i,j)\right)-kx.\hfill \end{array}$$

As written above, $F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)$ is a function of the independent variables $\beta$, $\lambda \left(i\right)$, $\mu \left(j\right)$, and x. However, under the saddle point approximation, the variables $\lambda \left(i\right)$, $\mu \left(j\right)$, and x are constrained by the conditions $$\begin{array}{c}\hfill \frac{\partial F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial \lambda \left(i\right)}=\frac{\partial F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial \mu \left(j\right)}=\frac{\partial F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)}{\partial x}=0,\end{array}$$ and the free energy under those constraints is written as $F_{\beta }^{MF}$. In the following, we use the notations $\frac{d{F}_{\beta }^{MF}}{d\beta }$ and $\frac{\partial F_{\beta }^{MF}}{\partial \beta }$ to differentiate between the total derivative and partial derivative of $F_{\beta }^{MF}$ with respect to $\beta$, respectively. It can be shown easily that (see Appendix C of [15]) $$\begin{array}{c}\hfill \frac{d{F}_{\beta }^{MF}}{d\beta }=\frac{\partial F_{\beta }\left({\mathit{\lambda }}^{MF},{\mathit{\mu }}^{MF},x^{MF}\right)}{\partial \beta },\end{array}$$ namely that the total derivative with respect to $\beta$ is, in this specific case, equal to the corresponding partial derivative, which is easily computed to be $$\begin{array}{ccc}& & \frac{\partial F_{\beta }^{MF}}{\partial \beta }=\hfill \\ & & \frac{1}{{\beta }^2}\sum _{i}G(\beta {\lambda }^{MF}\left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))-\hfill \\ & & \frac{1}{\beta }{\lambda }^{MF}\left(i\right)g(\beta {\lambda }^{MF}\left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))+\hfill \\ & & \frac{1}{{\beta }^2}\sum _{j}G(\beta {\mu }^{MF}\left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))-\hfill \\ & & \frac{1}{\beta }{\mu }^{MF}\left(j\right)g(\beta {\mu }^{MF}\left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))\hfill \\ & & -\frac{1}{{\beta }^2}\sum _{i,j}\left[H\left(\beta y^{MF}\left(i\right)\right)-\beta y^{MF}(i,j)h\left(\beta y^{MF}(i,j)\right)\right].\hfill \end{array}$$

Let $t\left(x\right)=H\left(x\right)-xh\left(x\right)$. The function $t\left(x\right)$ is continuous and defined over all real values x and is bounded above by 0, i.e., $t\left(x\right)\le 0\forall x\in \mathbb{R}$. Similarly, let us define $l(x,a,b)=G(x,a,b)-xg(x,a,b)$. In Appendix H, we show that $l(x,a,b)$ is positive over $\mathbb{R}$. As (A5)$$\begin{array}{ccc}\hfill \frac{d{F}_{\beta }^{MF}}{d\beta }=& & \frac{1}{{\beta }^2}\sum _{i}l(\beta {\lambda }^{MF}\left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))+\hfill \\ & & \frac{1}{{\beta }^2}\sum _{j}l(\beta {\mu }^{MF}\left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))\hfill \\ & & -\frac{1}{{\beta }^2}\sum _{i,j}t\left(\beta y^{MF}(i,j)\right),\hfill \end{array}$$ we conclude that $$\begin{array}{c}\hfill \frac{d{F}_{\beta }^{MF}}{d\beta }\ge 0,\end{array}$$ namely that $F_{\beta }^{MF}$ is a monotonically increasing function of $\beta$.

Appendix D. Monotonicity of $U_{\beta }^{MF}$

Let $$\begin{array}{c}\hfill U_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=\sum _{i,j}C(i,j)\overline{G}(i,j),\end{array}$$ and let the corresponding mean field approximation of the internal energy at the saddle point, $$\begin{array}{c}\hfill U_{\beta }^{MF}=U_{\beta }\left({\mathit{\lambda }}^{MF},{\mathit{\mu }}^{MF},x^{MF}\right).\end{array}$$

Before computing $\frac{d{U}_{\beta }^{MF}}{d\beta }$, we prove the following property

Proposition A1.

(A6) $$\begin{array}{ccc}\hfill U_{\beta }^{MF}& =& F_{\beta }^{MF}+\beta \frac{d{F}_{\beta }^{MF}}{d\beta },\hfill \end{array}$$ i.e., it extends the well-known relationship between the free energy and the average energy to their mean field counterparts.

Proof. 

The proof of this proposition is virtually identical to the proof of the same results for the assignment problem, found in Appendix C of [15]. □

Based on the chain rule, $$\begin{array}{ccc}\hfill \frac{d{U}_{\beta }^{MF}}{d\beta }=& & \frac{\partial U_{\beta }^{MF}}{\partial \beta }+\sum _i\frac{\partial U_{\beta }^{MF}}{\partial \lambda \left(i\right)}\frac{\partial \lambda \left(i\right)}{\partial \beta }+\sum _j\frac{\partial U_{\beta }^{MF}}{\partial \mu \left(j\right)}\frac{\partial \mu \left(j\right)}{\partial \beta }+\hfill \\ & & \frac{\partial U_{\beta }^{MF}}{\partial x}\frac{\partial x}{\partial \beta }.\hfill \end{array}$$

Using Proposition A1, we find that the partial derivatives of $U_{\beta }^{MF}$ with respect to $\lambda$, $\mu$, and x are all zeros, and therefore, $$\begin{array}{ccc}\hfill \frac{d{U}_{\beta }^{MF}}{d\beta }=& & \frac{\partial U_{\beta }^{MF}}{\partial \beta }.\hfill \end{array}$$

Using again Proposition A1, we obtain $$\begin{array}{ccc}\hfill \frac{d{U}_{\beta }^{MF}}{d\beta }& =& \frac{d{F}_{\beta }^{MF}}{d\beta }+\frac{d}{d\beta }\left(\beta \frac{d{F}_{\beta }^{MF}}{d\beta }\right)\hfill \\ & =& 2\frac{\partial F_{\beta }^{MF}}{\partial \beta }+\beta \frac{\partial }{\partial \beta }\left(\frac{\partial F_{\beta }^{MF}}{\partial \beta }\right).\hfill \end{array}$$

Using Equation (A5), $$\begin{array}{ccc}\hfill {\beta }^2\frac{d{U}_{\beta }^{MF}}{d\beta }& & =\sum _i\beta \lambda \left(i\right)l^{\prime }(\beta {\lambda }^{MF}\left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))\hfill \\ & & +\sum _j\beta {\mu }^{MF}\left(j\right)l^{\prime }(\beta {\mu }^{MF}\left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))\hfill \\ & & -\sum _{i,j}\beta y^{MF}(i,j)t^{\prime }\left(\beta y^{MF}(i,j)\right),\hfill \end{array}$$ where $t\left(x\right)$ and $l(x,a,b)$ are the functions defined above, and $t^{\prime }\left(x\right)$ and $l^{\prime }(x,a,b)$ are their derivatives with respect to x. Let us define $r\left(x\right)=x{t}^{\prime }\left(x\right)$ and $s(x,a,b)=x{l}^{\prime }(x,a,b)$. Then, $$\begin{array}{ccc}\hfill {\beta }^2\frac{\partial U_{\beta }^{MF}}{\partial \beta }=& & \sum _{i}s(\beta {\lambda }^{MF}\left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))+\hfill \\ & & \sum _{j}s(\beta {\mu }^{MF}\left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))-\hfill \\ & & \sum _{i,j}r\left(\beta y^{MF}(i,j)\right).\hfill \end{array}$$

Note that ${\displaystyle r\left(x\right)=\frac{x^2{e}^x}{{(1+e^x)}^2}}$. Therefore, $r\left(x\right)$ is positive, bounded below by 0. In Appendix H, we show that $s(x,a,b)$ is negative, bounded above 0. Therefore, $$\begin{array}{c}\hfill \frac{d{U}_{\beta }^{MF}}{d\beta }=\frac{\partial U_{\beta }^{MF}}{\partial \beta }\le 0,\end{array}$$ and the function $U_{\beta }^{MF}$ is a monotonically decreasing function of $\beta$.

Appendix E. Proof of Theorem 2: Convergence of the Mean Field Free Energy and the Internal Energy to the Optimal Assignment Cost

For simplicity in notation, we define $F^{MF}(\infty )={lim}_{\beta \to +\infty }{F}_{\beta }^{MF}$ and $U^{MF}(\infty )={lim}_{\beta \to +\infty }{U}_{\beta }^{MF}$.

Appendix E.1. Defining Entropy

Recall that we have defined $$\begin{array}{ccc}\hfill y(i,j)& =& C(i,j)+\lambda \left(i\right)+\mu \left(j\right)+x,\hfill \\ \hfill X(i,j)& =& h\left(\beta y\right(i,j\left)\right),\hfill \\ \hfill d_1\left(i\right)& =& g(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right)),\hfill \\ \hfill d_2\left(j\right)& =& g(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right)).\hfill \end{array}$$

The free energy is given by $$\begin{array}{ccc}\hfill F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)& & =-\frac{1}{\beta }\sum _{i}G(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))\hfill \\ & & -\frac{1}{\beta }\sum _{j}G(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))\hfill \\ & & +\frac{1}{\beta }\sum _{i,j}H\left(\beta y(i,j)\right)-kx,\hfill \end{array}$$ while the internal energy is $$\begin{array}{c}\hfill U_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=\sum _{i,j}C(i,j)X(i,j).\end{array}$$

By adding and subsequently subtracting the internal energy in the equation defining the free energy, and then replacing $C(i,j)$ with $y(i,j)-\lambda \left(i\right)-\mu \left(j\right)-x$, we obtain (after some reorganization), $$\begin{array}{ccc}& & F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=U_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)\hfill \\ & & -\frac{1}{\beta }\sum _{i}l(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))\hfill \\ & & -\frac{1}{\beta }\sum _{j}l(\beta \mu \left(j\right),q_{min}\left(j\right),q_{max}\left(j\right))+\frac{1}{\beta }\sum _{i,j}t\left(\beta y(i,j)\right)\hfill \\ & & +x\left(\sum _{i,j}X(i,j)-k\right)+\sum _i\lambda \left(i\right)\left(\sum _{j}X(i,j)-d_1\left(i\right)\right)\hfill \\ & & +\sum _j\mu \left(j\right)\left(\sum _{i}X(i,j)-d_2\left(j\right)\right),\hfill \end{array}$$ where $t\left(x\right)$ and $l(x,a,b)$ are defined above. Let us define (A7)$$\begin{array}{c}\hfill S_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=\sum _{i}l(\beta \lambda \left(i\right),p_{min}\left(i\right),p_{max}\left(i\right))+\\ \hfill \sum _{j}l(\beta \mu \left(j\right),q_{min}\left(i\right),q_{max}\left(i\right))-\sum _{i,j}t\left(\beta y(i,j)\right).\end{array}$$

Then, (A8)$$\begin{array}{ccc}& & F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=U_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)-\frac{1}{\beta }{S}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)\hfill \\ & & +x\left(\sum _{i,j}X(i,j)-k\right)+\sum _i\lambda \left(i\right)\left(\sum _{j}X(i,j)-d_1\left(i\right)\right)+\hfill \\ & & \sum _j\mu \left(j\right)\left(\sum _{i}X(i,j)-d_2\left(j\right)\right).\hfill \end{array}$$

The form of the free energy given in Equation (A8) has an intuitive physical interpretation. The first term is the original unbalanced assignment energy, the second is -T times an entropy term (defined in Equation (A7)), and the third, fourth, and fifth terms impose the constraints via Lagrange multipliers. At the saddle point, these constraints are satisfied, and the free energy has the form (A9)$$\begin{array}{c}\hfill F_{\beta }^{MF}=U_{\beta }^{MF}-\frac{1}{\beta }{S}_{\beta }^{MF}.\end{array}$$

Let us find the bound for the entropy term. $t\left(x\right)$ is negative, bounded below by $-ln\left(2\right)$. Similarly, in Appendix H, we show that $l(x,a,b)$ is positive, bounded above by $ln(b-a+1)$. Let us define (A10)$$\begin{array}{ccc}\hfill A(N_1,N_2)=& & N_1{N}_{2}ln\left(2\right)+\hfill \\ & & \sum _{i=1}^{N_1}ln(p_{max}\left(i\right)-p_{min}\left(i\right)+1)+\hfill \\ & & \sum _{j=1}^{N_2}ln(q_{max}\left(j\right)-q_{min}\left(j\right)+1).\hfill \end{array}$$

Then, the entropy satisfies the following constraints:(A11)$$\begin{array}{c}\hfill 0\le S_{\beta }^{MF}\le A(N_1,N_2).\end{array}$$

Appendix E.2. F^MF (∞) = U^MF (∞)

Proof. 

Using Equation (A9), after rearrangements, we obtain $$\begin{array}{c}\hfill U_{\beta }^{MF}-\frac{1}{\beta }A(N_1,N_2)\le F_{\beta }^{MF}\le U_{\beta }^{MF}.\end{array}$$

Taking the limits when $\beta \to +\infty$, we obtain (A12)$$\begin{array}{c}\hfill F^{MF}(\infty )=U^{MF}(\infty ).\end{array}$$

□

Appendix E.3. U * ≤ F^MF (∞)

Proof. 

Let $U^{MF}\left(\beta \right)$ be the mean field internal energy at the inverse temperature $\beta$: $$\begin{array}{c}\hfill U_{\beta }^{MF}=\sum _{i,j}C(i,j)X_{\beta }^{MF}(i,j),\end{array}$$ where $X_{\beta }^{MF}$ is the solution the the SPA system of equations. At a finite inverse temperature $\beta$, $X_{\beta }^{MF}$ is strictly non-integral. In addition, $X_{\beta }^{MF}$ satisfies constraints on row sums and column sums that make it an element of ${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$; see appendix A for details. Based on Theorem A2, $X_{\beta }^{MF}$ can be written as a linear combination of matrices in ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$, $$\begin{array}{c}\hfill X_{\beta }^{MF}=\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\pi ,\end{array}$$ with all $a_{\pi }\in [0,1]$ and ${\sum }_{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }=1$. Therefore, (A13)$$\begin{array}{ccc}\hfill U_{\beta }^{MF}& =& \sum _{i,j}C(i,j)X_{\beta }^{MF}(i,j)\hfill \\ & =& \sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\sum _{i,j}C(i,j)\pi (i,j).\hfill \end{array}$$

As $U^{*}$ is the minimum matching cost over all possible matrices in ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$, we have (A14)$$\begin{array}{c}\hfill \sum _{i,j}C(i,j)\pi (i,j)\ge U^{*},\end{array}$$ for all $\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$. Combining Equations (A13) and (A14), we obtain $$\begin{array}{ccc}\hfill U_{\beta }^{MF}& \ge & \sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }{U}^{*}\hfill \\ & \ge & \left(\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\right)U^{*},\hfill \end{array}$$ from which we conclude that at each $\beta$, $$\begin{array}{c}\hfill U^{*}\le U_{\beta }^{MF}.\end{array}$$

Therefore $U^{*}\le U^{MF}(\infty )$, and consequently $U^{*}\le F^{MF}(\infty )$, based on Equation (A12). □

Appendix E.4. U * ≥ F^MF (∞)

Proof. 

Let us first recall the definition of the free energy, $$\begin{array}{ccc}\hfill F_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)& & =-\frac{1}{\beta }\sum _{i}ln\left(\frac{e^{\beta p_{min}\left(i\right)\lambda \left(i\right)}-e^{\beta (p_{max}\left(i\right)+1)\lambda \left(i\right)}}{1-e^{\beta \lambda \left(i\right)}}\right)\hfill \\ & & -\frac{1}{\beta }\sum _{j}ln\left(\frac{e^{\beta q_{min}\left(j\right)\mu \left(j\right)}-e^{\beta (q_{max}\left(j\right)+1)\mu \left(j\right)}}{1-e^{\beta \mu \left(j\right)}}\right)\hfill \\ & & -\frac{1}{\beta }\sum _{i,j}ln\left[1+e^{-\beta y(i,j)}\right]-kx.\hfill \end{array}$$

Note these two properties of limits: $$\begin{array}{c}\hfill \underset{\beta \to +\infty }{lim}\frac{ln(1+e^{-a\beta })}{\beta }=\left\{\begin{array}{cc}0\hfill & \mathrm{if}a\ge 0\hfill \\ -a\hfill & \mathrm{if}a\le 0,\hfill \end{array}\right.\end{array}$$ and $$\begin{array}{c}\hfill \underset{\beta \to +\infty }{lim}\frac{1}{\beta }ln\left(\frac{e^{xa\beta }-e^{x(b+1)\beta }}{1-e^{x\beta }}\right)=\left\{\begin{array}{cc}xb\hfill & \mathrm{if}x\ge 0\hfill \\ xa\hfill & \mathrm{if}x\le 0,\hfill \end{array}\right.\end{array}$$ when a and b are non-negative integers. Therefore, (A15)$$\begin{array}{ccc}& & \underset{\beta \to +\infty }{lim}{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right)=\hfill \\ & & -\sum _{i\left|\lambda \right(i)\ge 0}{p}_{max}\left(i\right)\lambda \left(i\right)-\sum _{i\left|\lambda \right(i)<0}{p}_{min}\left(i\right)\lambda \left(i\right)\hfill \\ & & -\sum _{j\left|\mu \right(j)\ge 0}{q}_{max}\left(j\right)\mu \left(j\right)-\sum _{j\left|\mu \right(j)<0}{q}_{min}\left(j\right)\mu \left(j\right)\hfill \\ & & +\sum _{(i,j)\left|y\right(i,j)\le 0}y(i,j)-kx.\hfill \end{array}$$

Let us consider now a matrix $\pi$ in ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$ (see above). We can write (A16)$$\begin{array}{ccc}\hfill \sum _{i,j}C(i,j)\pi (i,j)=& & \sum _{i,j}y(i,j)\pi (i,j)-\sum _i\lambda \left(i\right)RS\left(i\right)\hfill \\ & & -\sum _j\mu \left(j\right)CS\left(j\right)-kx,\hfill \end{array}$$ where $RS$ and $CS$ are the row sums and column sums of $\pi$, respectively. For each index i, the summand included in the first term on the right is always larger than or equal to the sum of all the corresponding terms that are negative: (A17)$$\begin{array}{c}\hfill \sum _{i,j}y(i,j)\pi (i,j)\ge \sum _{i\left|y\right(i,j)\le 0}y(i,j).\end{array}$$

As $\pi$ belongs to ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$, its row sums satisfy the constraints: (A18)$$\begin{array}{c}\hfill p_{min}\left(i\right)\le RS\left(i\right)\le p_{max}\left(i\right).\end{array}$$

We multiply this equation separately for positive and negative $\lambda \left(i\right)$. We find $$\begin{array}{c}\hfill \sum _{i\left|\lambda \right(i)>0}{p}_{min}\left(i\right)\lambda \left(i\right)\le \sum _{i\left|\lambda \right(i)>0}\lambda \left(i\right)RS\left(i\right)\le \sum _{i\left|\lambda \right(i)>0}{p}_{max}\left(i\right)\lambda \left(i\right)\\ \hfill \sum _{i\left|\lambda \right(i)<0}{p}_{max}\left(i\right)\lambda \left(i\right)\le \sum _{i\left|\lambda \right(i)<0}\lambda \left(i\right)RS\left(i\right)\le \sum _{i\left|\lambda \right(i)<0}{p}_{min}\left(i\right)\lambda \left(i\right).\end{array}$$

Therefore, (A19)$$\begin{array}{c}\hfill -\sum _i\lambda \left(i\right)RS\left(i\right)\ge -\sum _{i\left|\lambda \right(i)\ge 0}{p}_{max}\left(i\right)\lambda \left(i\right)-\sum _{i\left|\lambda \right(i)<0}{p}_{min}\left(i\right)\lambda \left(i\right).\\ \hfill \end{array}$$

Similarly, we can show that (A20)$$\begin{array}{c}\hfill -\sum _j\mu \left(j\right)CS\left(j\right)\ge -\sum _{j\left|\mu \right(j)\ge 0}{q}_{max}\left(j\right)\mu \left(j\right)-\sum _{j\left|\mu \right(j)<0}{q}_{min}\left(j\right)\mu \left(j\right).\\ \hfill \end{array}$$

Combining Equations (A15)–(A17), (A19) and (A20), we obtain (A21)$$\begin{array}{c}\hfill \sum _{i,j}C(i,j)\pi (i,j)\ge \underset{\beta \to +\infty }{lim}{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right).\end{array}$$

Equation (A21) is valid for all matrices $\pi$ in ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$. It is therefore valid for the optimal ${\pi }^{*}$ that solves the general assignment problem. Since $U^{*}={\sum }_{i,j}C(i,j){\pi }^{*}(i,j)$, we have $$\begin{array}{c}\hfill U^{*}\ge \underset{\beta \to +\infty }{lim}{F}_{\beta }\left(\mathit{\lambda },\mathit{\mu },x\right).\end{array}$$

As this equation is true for all $\mathit{\lambda }$, $\mathit{\mu }$, and x, it is true in particular for $\mathit{\lambda }={\mathit{\lambda }}^{MF}$, $\mathit{\mu }={\mathit{\mu }}^{MF}$, and $x=x^{MF}$, leading to $$\begin{array}{c}\hfill U^{*}\ge \underset{\beta \to +\infty }{lim}{F}_{\beta }^{MF}=F^{MF}(\infty ).\end{array}$$

□

We have shown that $U^{*}\le F^{MF}(\infty )$ and $F^{MF}(\infty )\le U^{*}$, and therefore $U^{*}=F^{MF}(\infty )$. The corresponding result for the internal energy, $U^{*}=U^{MF}(\infty )$ follows directly from Equation (A12).

Appendix F. Proof of Theorem 3: Bounds on the Internal Energy and Free Energy

Appendix F.1. Bounds on the Free Energy

Proposition A1 from Appendix D states that $$\begin{array}{ccc}\hfill \beta \frac{d{F}_{\beta }^{MF}}{d\beta }& =& -F_{\beta }^{MF}+U_{\beta }^{MF}.\hfill \end{array}$$

Using this equation and the relationship between free energy, energy, and entropy at SPA (see Equation (A22)), we obtain $$\begin{array}{c}\hfill \frac{d{F}_{\beta }^{MF}}{d\beta }=\frac{1}{{\beta }^2}{S}_{\beta }^{MF}.\end{array}$$

From the bounds on the entropy (Equation (A11)), $$\begin{array}{c}\hfill 0\le \frac{d{F}_{\beta }^{MF}}{d\beta }\le \frac{A(N_1,N_2)}{{\beta }^2}.\end{array}$$

By integrating over $\beta$ between $+\infty$ and $\beta$, $$\begin{array}{c}\hfill 0\le F^{MF}(\infty )-F_{\beta }^{MF}\le \frac{A(N_1,N_2)}{\beta }.\end{array}$$

Finally, as $F^{MF}(\infty )=U^{*}$, (A22)$$\begin{array}{c}\hfill U^{*}-\frac{A(N_1,N_2)}{\beta }\le F_{\beta }^{MF}\le U^{*}.\end{array}$$

Appendix F.2. Bounds on the Energy

As $U_{\beta }^{MF}=F_{\beta }^{MF}+\frac{1}{\beta }{S}_{\beta }^{MF}$, using the bounds on the free energy and on the entropy, we obtain (A23)$$\begin{array}{c}\hfill U_{\beta }^{MF}\le U^{*}+\frac{A(N_1,N_2)}{\beta }.\end{array}$$

In addition, $U_{\beta }^{MF}$ is monotonic, decreasing, with the limit $U^{*}$ as $\beta \to +\infty$, $U^{*}\le U_{\beta }^{MF}$. Therefore, $$\begin{array}{c}\hfill U^{*}\le U_{\beta }^{MF}\le U^{*}+\frac{A(N_1,N_2)}{\beta }.\end{array}$$

Appendix G. Proof of Theorem 4: Bounds on the Transportation Matrix ${\overline{G}}_{\beta }$

This proof is inspired by the proof of Theorem 6 in Appendix 2 of [35] and by Appendix F of [15].

We first recall that ${\overline{G}}_{\beta }$ belongs to ${\mathcal{U}}^k(R_{min}^{max},C_{min}^{max})$ (see above), with $\sigma \left({\overline{G}}_{\beta }\right)=k$. As such, it can be written as a linear combination of the matrices $\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$, $$\begin{array}{c}\hfill {\overline{G}}_{\beta }=\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\pi ,\end{array}$$ with all $a_{\pi }\in [0,1]$ and $\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }=1$.

We assume that the general assignment problem considered has a unique solution. We want to prove that $\underset{i,j}{max}\left|{\overline{G}}_{\beta }(i,j)-G^{*}(i,j)\right|\le \frac{A(N_1,N_2)}{\beta \Delta }$, where $G^{*}$ is the optimal solution of the assignment problem, $\Delta =U^{2*}-U^{*}$ is the difference in total cost between the second best solution and the optimal solution, and $A(N_1,N_2)$ is defined in Equation (A10). We use for that a proof by contradiction. We assume that there exists a pair $(i,j)$ such that $$\begin{array}{c}\hfill \frac{A(N_1,N_2)}{\beta \Delta }<\left|{\overline{G}}_{\beta }(i,j)-G^{*}(i,j)\right|.\end{array}$$

Let us denote $B(i,j)=\left|{\overline{G}}_{\beta }(i,j)-G^{*}(i,j)\right|$. As $G^{*}$ is a binary matrix, $G^{*}(i,j)=0$ or $G^{*}(i,j)=1$.

In the first case, $$\begin{array}{ccc}\hfill B(i,j)& =& {\overline{G}}_{\beta }(i,j)\hfill \\ & =& \sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\pi (i,j).\hfill \end{array}$$

Since $G^{*}$ belongs to ${\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})$, it is included in the decomposition of ${\overline{G}}_{\beta }$, and therefore, $$\begin{array}{ccc}\hfill B(i,j)& =& a_{G^{*}}{G}^{*}(i,j)+\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\pi (i,j)\hfill \\ & =& \sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\pi (i,j)\hfill \\ & <& \sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max}))-\left\{G^{*}\right\}}{a}_{\pi }=1-a_{G^{*}},\hfill \end{array}$$ where the final equality follows from the fact that the sum of all coefficients a is equal to 1.

In the second case, $G^{*}(i,j)=1$, $$\begin{array}{ccc}\hfill B(i,j)& =& 1-{\overline{G}}_{\beta }(i,j)\hfill \\ & =& 1-\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\pi (i,j).\hfill \end{array}$$

Again, as $G^{*}$ is included in the decomposition of ${\overline{G}}_{\beta }$, $$\begin{array}{ccc}\hfill B(i,j)& =& 1-a_{G^{*}}{G}^{*}(i,j)-\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\pi (i,j)\hfill \\ & =& 1-a_{G^{*}}-\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\pi (i,j)\hfill \\ & <& 1-a_{G^{*}},\hfill \end{array}$$ where the final inequality follows from the fact that ${\sum }_{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\pi (i,j)$ is positive.

In conclusion, in both cases, we have (A24)$$\begin{array}{c}\hfill \frac{A(N_1,A_2)}{\beta \Delta }<1-a_{G^{*}}.\end{array}$$

Now, let us look at the energy associated with ${\overline{G}}_{\beta }$:$$\begin{array}{ccc}\hfill U_{\beta }^{MF}& =& \sum _{i,j}C(i,j){\overline{G}}_{\beta }(i,j)\hfill \\ & =& \sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})}{a}_{\pi }\sum _{i,j}C(i,j)\pi (i,j)\hfill \\ & =& a_{G^{*}}{U}^{*}+\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\sum _{i,j}C(i,j)\pi (i,j)\hfill \\ & \ge & a_{G^{*}}{U}^{*}+\left(\sum _{\pi \in {\mathcal{P}}^k(R_{min}^{max},C_{min}^{max})-\left\{G^{*}\right\}}{a}_{\pi }\right)U^{2*}\hfill \\ & \ge & a_{G^{*}}{U}^{*}+(1-a_{G^{*}})U^{2*}\hfill \\ & \ge & U^{*}+(1-a_{G^{*}})\Delta .\hfill \end{array}$$

In Theorem 3, we have shown that $$\begin{array}{c}\hfill U^{*}\le U_{\beta }^{MF}\le U^{*}+\frac{A(N_1,N_2)}{\beta }.\end{array}$$

Therefore, $$\begin{array}{c}\hfill U^{*}+(1-a_{G^{*}})\Delta \le U^{*}+\frac{A(N_1,N_2)}{\beta },\end{array}$$ i.e., (A25)$$\begin{array}{c}\hfill (1-a_{G^{*}})\le \frac{A(N_1,N_2)}{\beta \Delta },\end{array}$$ as $\Delta$ is strictly positive.

We have shown that $\frac{A(N_1,N_2)}{\beta \Delta }<1-a_{G^{*}}$ (Equation (A24)) and $(1-a_{G^{*}})\le \frac{A(N_1,N_2)}{\beta \Delta }$ (Equation (A25)); i.e., we have reached a contradiction. Our hypothesis is wrong, and therefore $\underset{i,j}{max}\left|{\overline{G}}_{\beta }(i,j)-G^{*}(i,j)\right|\le \frac{A(N_1,N_2)}{\beta \Delta }$.

Appendix H. Properties of the General Function G(x,a,b) and Its Derivatives

Theorem A3.

Let $G(x,a,b)$ be the function defined by, $$\begin{array}{ccc}\hfill G(x,a,b)& =& ln\left(\frac{e^{ax}-e^{(b+1)x}}{1-e^x}\right).\hfill \end{array}$$

where a and b are two positive integers with $0\le a<b$. Let $g\left(x\right)=G^{\prime }\left(x\right)$, $l(x,a,b)=G(x,a,b)-xg(x,a,b)$, and $s(x,a,b)=x{l}^{\prime }(x,a,b)$. We have the following properties:

• (1) 

  $g^{\prime }(x,a,b)>0\forall x\in \mathbb{R}$.

• (2) 

  $0\le l(x,a,b)\le ln(b-a+1)\forall x\in \mathbb{R}$.

• (3) 

  $s(x,a,b)\le 0\forall x\in \mathbb{R}$.

Note that in Theorem A3, we consider a, b non-negative, with $a<b$. The result can be trivially expanded to the case $a=b$, with $G(x,a,a)=ax$ and therefore $g(x,a,a)=a$, $g^{\prime }(x,a,a)=0$, $l(x,a,a)=0$, and $s(x,a,a)=0$. In the following, we only consider $a<b$.

Appendix H.1. The Function G(x,a,b)

Note first that $$\begin{array}{ccc}\hfill G(x,a,b)& =& ln\left(\frac{e^{ax}-e^{(b+1)x}}{1-e^x}\right)\hfill \\ & =& ln\left(\sum _{k=a}^b{e}^{kx}\right).\hfill \end{array}$$ $G(x,a,b)$ is therefore defined, continuous, and differentiable over $\mathbb{R}$. Some special values:$$\begin{array}{ccc}\hfill G(0,a,b)& =& ln(b+a-1),\hfill \\ \hfill \underset{x\to -\infty }{lim}G(x,a,b)& =& -\infty ,\hfill \\ \hfill \underset{x\to +\infty }{lim}G(x,a,b)& =& +\infty .\hfill \end{array}$$

Appendix H.2. g(x,a,b) and Its Derivative g ′ (x,a,b)

Note that $$\begin{array}{ccc}\hfill g(x,a,b)=& & \frac{(b-a+1)}{e^{(b-a+1)x}-1}+\frac{1}{1-e^x}+b,\hfill \\ \hfill g^{\prime }(x,a,b)=& & -\frac{{(b-a+1)}^2{e}^{(b-a+1)x}}{{(e^{(b-a+1)x}-1)}^2}+\frac{e^x}{{(e^x-1)}^2},\hfill \end{array}$$ g and its derivative $g^{\prime }$ are defined when $e^{(b-a+1)x}-1\ne 0$ and $e^x-1\ne 0$, i.e., when $x\ne 0$. They can be extended by continuity at 0:$$\begin{array}{ccc}\hfill g(0,a,b)=& & \frac{a+b}{2},\hfill \\ \hfill g^{\prime }(0,a,b)=& & \frac{{(b-a+1)}^2-1}{12}.\hfill \end{array}$$

We have $g^{\prime }(0,a,b)>0$ as $b-a+1>1$. In addition, $$\begin{array}{c}\hfill \underset{x\to -\infty }{lim}g(x,a,b)=a,\\ \hfill \underset{x\to +\infty }{lim}g(x,a,b)=b,\\ \hfill \underset{x\to -\infty }{lim}{g}^{\prime }(x,a,b)=0,\\ \hfill \underset{x\to +\infty }{lim}{g}^{\prime }(x,a,b)=0.\end{array}$$

Figure A1 shows a few examples of g and $g^{\prime }$ for different values of a and b.

Figure A1. Examples of [Forumla omitted. See PDF.] (A) and its derivative [Forumla omitted. See PDF.] (B) for [Forumla omitted. See PDF.] (black), [Forumla omitted. See PDF.], (red), and [Forumla omitted. See PDF.] (blue).

We now prove property 1 of Theorem A3.

Proof. 

Let us define $m=b-a+1$. As $a<b$, we have $m>1$. Let x be a real number; then, $$\begin{array}{ccc}\hfill g^{\prime }(x,a,b)=& & -\frac{m^2{e}^{mx}}{{(e^{mx}-1)}^2}+\frac{e^x}{{(e^x-1)}^2},\hfill \end{array}$$ which we can rewrite as $$\begin{array}{ccc}& & g^{\prime }(x,a,b)\hfill \\ & =& \frac{1}{4}\left(-\frac{m^2}{{sinh}^2\left(\frac{mx}{2}\right)}+\frac{1}{{sinh}^2\left(\frac{x}{2}\right)}\right)\hfill \\ & =& \frac{1}{4}\left(\frac{{sinh}^2\left(\frac{mx}{2}\right)-m^2{sinh}^2\left(\frac{x}{2}\right)}{{sinh}^2\left(\frac{mx}{2}\right){sinh}^2\left(\frac{x}{2}\right)}\right)\hfill \\ & =& \frac{1}{4}\left(\frac{(sinh\left(\frac{mx}{2}\right)-msinh\left(\frac{x}{2}\right))(sinh\left(\frac{mx}{2}\right)+msinh\left(\frac{x}{2}\right))}{{sinh}^2\left(\frac{mx}{2}\right){sinh}^2\left(\frac{x}{2}\right)}\right).\hfill \end{array}$$

From this formulation, it is clear that $g^{\prime }(x,a,b)$ is even. As $g^{\prime }(0,a,b)$ is strictly positive (see above), we only need to consider $x>0$. Let us define $y=x/2$. For $y>0$, the denominator of $g^{\prime }(x,a,b)$ and the term $sinh\left(my\right)+msinh\left(y\right)$ are strictly positive. The only remaining term is $sinh\left(my\right)-msinh\left(y\right)$, which we call $f\left(y\right)$. We have $$\begin{array}{ccc}\hfill f^{\prime }\left(y\right)& =& m\left(cosh\left(my\right)-cosh\left(y\right)\right)\hfill \\ & =& 2msinh\left(\right(m+1\left)y\right)sinh\left(\right(m-1\left)y\right).\hfill \end{array}$$

As $m>1$, $sinh\left(\right(m+1\left)y\right)$ and $sinh\left(\right(m-1\left)y\right)$ are strictly positive when $y>0$: $f\left(y\right)$ is monotonically increasing for $y\ge 0$. As $f\left(0\right)=0$, we obtain $f\left(y\right)>0$ for $y>0$. Therefore $g^{\prime }(x,a,b)>0$ for all $x>0$, and since it is an even function and $g^{\prime }(0,a,b)>0$, we find that $\forall x\in \mathbf{R}{g}^{\prime }(x,a,b)>0$. □

Appendix H.3. The Function l(x,a,b) = G(x,a,b)-xg(x,a,b)

Recall that $$\begin{array}{ccc}& & l(x,a,b)=G(x,a,b)-xg(x,a,b).\hfill \end{array}$$

We nowprove property 2 of Theorem A3.

Proof. 

We prove first that $l(x,a,b)$ is even. Note that $$\begin{array}{ccc}\hfill G(-x,a,b)& =& ln\left(\frac{e^{-ax}-e^{-(b+1)x}}{1-e^{-x}}\right)\hfill \\ & =& ln\left(\frac{e^{(b+1)x}-e^{ax}}{e^x-1}\right)+ln\left(\frac{e^x}{e^{(b+1)x}{e}^{ax}}\right)\hfill \\ & =& -bx-ax+G(x,a,b),\hfill \end{array}$$ and $$\begin{array}{ccc}\hfill g(-x,a,b)& =& \frac{(b-a+1)}{e^{-(b-a+1)x}-1}+\frac{1}{1-e^{-x}}+b\hfill \\ & =& \frac{(b-a+1)e^{(b-a+1)x}}{1-e^{(b-a+1)x}}+\frac{e^x}{e^x-1}+b\hfill \\ & =& a+b-g(x,a,b).\hfill \end{array}$$

Therefore, $$\begin{array}{ccc}\hfill l(-x,a,b)& =& G(-x,a,b)+xg(-x,a,b)\hfill \\ & =& G(x,a,b)-xg(x,a,b)\hfill \\ & =& l(x,a,b).\hfill \end{array}$$

Let now x be a positive real number. Then, $$\begin{array}{ccc}\hfill l^{\prime }(x,a,b)& =& G^{\prime }(x,a,b)-g(x,a,b)-x{g}^{\prime }(x,a,b).\hfill \\ & =& -x{g}^{\prime }(x,a,b)\hfill \end{array}$$

We know that $g^{\prime }(x,a,b)$ is positive for all real numbers x (property 1 of Theorem A3); therefore, $l^{\prime }(x,a,b)$ is negative on ${\mathbb{R}}^{+}$ and $l(x,a,b)$ is monotonic decreasing on this domain. Since: $$\begin{array}{ccc}\hfill l(0,a,b)& =& ln(b-a+1),\hfill \\ \hfill \underset{x\to +\infty }{lim}l(x,a,b)& =& 0,\hfill \end{array}$$ we have $0\le l(x,a,b)\le ln(b-a+1)$ for $x\in {\mathbb{R}}^{+}$, and since $l(x,a,b)$ is even, $0\le l(x,a,b)\le ln(b-a+1)$ for $x\in \mathbb{R}$, which concludes the proof of property 2 of Theorem A3. □

Appendix H.4. The Function s(x,a,b) Is Negative

We now prove the property 3 of A3.

Proof. 

Recall that $$\begin{array}{ccc}\hfill l^{\prime }(x,a,b)& =& -x{g}^{\prime }(x,a,b).\hfill \end{array}$$

Therefore, $$\begin{array}{ccc}\hfill s(x,a,b)& =& -x^2{g}^{\prime }(x,a,b).\hfill \end{array}$$

We know that $g^{\prime }(x,a,b)$ is positive for all real numbers x (property 1 of Theorem A3); therefore, $s(x,a,b)$ is negative on $\mathbb{R}$. □

References

1. Dell’Amico, M.; Toth, P. Algorithms and codes for dense assignment problems: The state of the art. Discret. Appl. Math.; 2000; 100, pp. 17-48. [DOI: https://dx.doi.org/10.1016/S0166-218X(99)00172-9]

2. Pentico, D.W. Assignment problems: A golden anniversary survey. Eur. J. Oper. Res.; 2007; 176, pp. 774-793. [DOI: https://dx.doi.org/10.1016/j.ejor.2005.09.014]

3. Burkard, R.; Dell’Amico, M.; Martello, S. Assignment Problems; Society for Industrial and Applied Mathematics (SIAM): Philadelphia, PA, USA, 2009.

4. Dantzig, G. Origins of the simplex method. A History of Scientific Computing; Association for Computing Machinery: New York, NY, USA, 1990; pp. 141-151.

5. Jacobi, C. De investigando ordine systematis aequationum differentialum vulgarium cujuscunque. J. Reine Angew. Math.; 1890; 94, pp. 292-320.

6. Kuhn, H. The Hungarian method for the assignment problem. Nav. Res. Logist.; 1955; 2, pp. 83-97. [DOI: https://dx.doi.org/10.1002/nav.3800020109]

7. Date, K.; Nagi, R. GPU-accelerated Hungarian algorithms for the Linear Assignment Problem. Parallel Comput.; 2016; 57, pp. 52-72. [DOI: https://dx.doi.org/10.1016/j.parco.2016.05.012]

8. Lopes, P.A.; Yadav, S.S.; Ilic, A.; Patra, S.K. Fast block distributed CUDA implementation of the Hungarian algorithm. J. Parallel Distrib. Comput.; 2019; 130, pp. 50-62. [DOI: https://dx.doi.org/10.1016/j.jpdc.2019.03.014]

9. Yadav, S.S.; Lopes, P.A.C.; Ilic, A.; Patra, S.K. Hungarian algorithm for subcarrier assignment problem using GPU and CUDA. Int. J. Commun. Syst.; 2019; 32, e3884. [DOI: https://dx.doi.org/10.1002/dac.3884]

10. Kumar, A. A modified method for solving the unbalanced assignment problems. Appl. Math. Comput.; 2006; 176, pp. 76-82. [DOI: https://dx.doi.org/10.1016/j.amc.2005.09.056]

11. Yadaiah, V.; Haragopal, V. A New Approach of Solving Single Objective Unbalanced Assignment Problem. Am. J. Oper. Res.; 2016; 6, pp. 81-89. [DOI: https://dx.doi.org/10.4236/ajor.2016.61011]

12. Costa, D.; Hertz, A. Ants can color graphs. J. Oper. Res. Soc.; 1997; 48, pp. 295-305. [DOI: https://dx.doi.org/10.1057/palgrave.jors.2600357]

13. Wang, L.; He, Z.; Liu, C.; Chen, Q. Graph based twin cost matrices for unbalanced assignment problem with improved ant colony algorithm. Results Appl. Math.; 2021; 12, 100207. [DOI: https://dx.doi.org/10.1016/j.rinam.2021.100207]

14. Betts, N.; Vasko, F.J. Solving the unbalanced assignment problem: Simpler is better. Am. J. Oper. Res.; 2016; 6–9, 296. [DOI: https://dx.doi.org/10.4236/ajor.2016.64028]

15. Koehl, P.; Orland, H. Fast computation of exact solutions of generic and degenerate assignment problems. Phys. Rev. E; 2021; 103, 042101. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.042101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34005932]

16. Koehl, P.; Delarue, M.; Orland, H. Physics approach to the variable-mass optimal-transport problem. Phys. Rev. E; 2021; 103, 012113. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.012113]

17. Birkhoff, G. Tres observaciones sobre el algebra lineal. Univ. Nac. Tucuman. Ser. A; 1946; 5, pp. 147-154.

18. Von Neumann, J. A certain zero-sum two-person game equivalent to the optimal assignment problem. Contrib. Theory Games; 1953; 2, pp. 5-12.

19. Dell’Amico, M.; Martello, S. The k-cardinality assignment problem. Discret. Appl. Math.; 1997; 76, pp. 103-121. [DOI: https://dx.doi.org/10.1016/S0166-218X(97)00120-0]

20. Dell’Amico, M.; Lodi, A.; Martello, S. Efficient algorithms and codes for k-cardinality assignment problems. Discret. Appl. Math.; 2001; 110, pp. 25-40. [DOI: https://dx.doi.org/10.1016/S0166-218X(00)00301-2]

21. Volgenant, A. Solving the k-cardinality assignment problem by transformation. Eur. J. Oper. Res.; 2004; 157, pp. 322-331. [DOI: https://dx.doi.org/10.1016/S0377-2217(03)00205-4]

22. Chandrasekaran, R.; Kabadi, S.; Murty, K. Some NP-complete problems in linear programming. Oper. Res. Lett.; 1982; 1, pp. 101-104. [DOI: https://dx.doi.org/10.1016/0167-6377(82)90006-2]

23. Greenberg, H.J. An analysis of degeneracy. Nav. Res. Logist. Q.; 1986; 33, pp. 635-655. [DOI: https://dx.doi.org/10.1002/nav.3800330409]

24. Megiddo, N.; Chandrasekaran, R. On the ε-perturbation method for avoiding degeneracy. Oper. Res. Lett.; 1989; 8, pp. 305-308. [DOI: https://dx.doi.org/10.1016/0167-6377(89)90014-X]

25. Koehl, P.; Delarue, M.; Orland, H. Finite temperature optimal transport. Phys. Rev. E; 2019; 100, 013310. [DOI: https://dx.doi.org/10.1103/PhysRevE.100.013310] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31499816]

26. Majumdar, J.; Bhunia, A.K. An alternative approach for unbalanced assignment problem via genetic algorithm. Appl. Math. Comput.; 2012; 218, pp. 6934-6941. [DOI: https://dx.doi.org/10.1016/j.amc.2011.12.070]

27. Rabbani, Q.; Khan, A.; Quddoos, A. Modified Hungarian method for unbalanced assignment problem with multiple jobs. Appl. Math. Comput.; 2019; 361, pp. 493-498. [DOI: https://dx.doi.org/10.1016/j.amc.2019.05.041]

28. Heist, N.; Hertling, S.; Ringler, D.; Paulheim, H. Knowledge Graphs on the Web—An Overview. arXiv; 2020; arXiv: 2003.00719

29. Veličković, P. Everything is connected: Graph neural networks. Curr. Opin. Struct. Biol.; 2023; 79, 102538. [DOI: https://dx.doi.org/10.1016/j.sbi.2023.102538] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36764042]

30. Camacho, D.; Panizo-LLedot, A.; Bello-Orgaz, G.; Gonzalez-Pardo, A.; Cambria, E. The four dimensions of social network analysis: An overview of research methods, applications, and software tools. Inf. Fusion; 2020; 63, pp. 88-120. [DOI: https://dx.doi.org/10.1016/j.inffus.2020.05.009]

31. Ideker, T. Network genomics. Ernst Scher. Res. Found. Workshop; 2007; 61, pp. 89-115.

32. Pierskalla, W.P. The multidimensional assignment problem. Oper. Res.; 1968; 16, pp. 422-431. [DOI: https://dx.doi.org/10.1287/opre.16.2.422]

33. Balas, E.; Saltzman, M.J. Facets of the three-index assignment polytope. Discret. Appl. Math.; 1989; 23, pp. 201-229. [DOI: https://dx.doi.org/10.1016/0166-218X(89)90014-0]

34. Spieksma, F.C. Multi index assignment problems: Complexity, approximation, applications. Nonlinear Assignment Problems: Algorithms and Applications; Springer: Berlin/Heidelberg, Germany, 2000; pp. 1-12.

35. Kosowsky, J.; Yuille, A. The invisible hand algorithm: Solving the assignment problem with statistical physics. Neural Netw.; 1994; 7, pp. 477-490. [DOI: https://dx.doi.org/10.1016/0893-6080(94)90081-7]

36. Koehl, P. Extreme points of general transportation polytopes. arXiv; 2024; arXiv: 2404.16791

37. Mendelsohn, N.S.; Dulmage, A.L. The convex hull of sub-permutation matrices. Proc. Am. Math. Soc.; 1958; 9, pp. 253-254. [DOI: https://dx.doi.org/10.1090/S0002-9939-1958-0095128-8]

38. Brualdi, R.A.; Lee, G.M. On the truncated assignment polytope. Linear Algebra Its Appl.; 1978; 19, pp. 33-62. [DOI: https://dx.doi.org/10.1016/0024-3795(78)90004-6]