A General Statistical Physics Framework for Assignment Problems

: 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 β representing the inverse of temperature. As β increases, the free energy converges to the optimal assignment cost. Furthermore, we demonstrate that when β 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.


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 |T| = N 2 and |A| = 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 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, . . ., 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 ≤ 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(N 4 ) or O(N 3 ) 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 "multiassignment" 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 sections 2 to 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.

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 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 (i) points in S 2 , with n 1 (i) belonging to [p min (i), . . ., p max (i)], where p min (i) and p max (i) are non-negative integers satisfying 0 ≤ p min (i) ≤ p max (i) ≤ N 2 .Similarly, a point j is assigned to n 2 (j) points in S 1 , with n 2 (j) belonging to [q min (j), . . ., q max (j)], where q min (j) and q max (j) are non-negative integers integers satisfying 0 ≤ q min (j) ≤ q max (j) ≤ 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 (i) is the number of agents needed to perform task i, while n 2 (j) is the number of tasks that can be performed by agent j. p min (i) and p max (i) represent the minimum and maximum number of agents that are assigned to a task i, respectively, while q min (j) and q max (j) 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 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: where k is the actual number of task-agent pairs.The numbers p min (i), p max (i), q min (j), q max (j), and k are problem-specific and given.They do satisfy some constraints, such as 0 ≤ p min (i) ≤ p max (i) ≤ N 2 and 0 ≤ q min (j) ≤ q max (j) ≤ N 1 , as described above, as well as j=1 q min (j).We store those numbers in four vectors, P min = (p min (1), . . ., p min (N 1 )), P max = (p max (1), . . ., p max (N 1 )), Q min = (q min (1), . . ., q min (N 2 )), Q max = (q max (1), . . ., q max (N 2 )).
Equation (3) recap most of the standard linear assignment problems: (i) If k = N 1 = N 2 and p min (i) = p max (i) = q min (j) = q max (j) = 1 for all i and j, we recover the balanced assignment problem, (ii) If p min (i) = q min (j) = 0 and p max (i) = q max (j) = 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 (i) and n 2 (j) 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 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 temperaturedependent problem involving real variables.The integer optimal solution is then obtained as the temperature approaches zero.

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(G) linked with an assignment matrix G is given by In this equation, β is the inverse temperature, namely β = 1/k B T, where k B is the Boltzmann constant and T is the temperature, and Z(β) 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 where F (β) 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 The three sums impose that G(i, j), n 1 (i), and n 2 (j) 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, λ(i), and µ(j), with i ∈ [1, N 1 ] and j ∈ [1, N 2 ].The partition function is then given, after reorganization, by (up to a multiplicative constant) by where ß is the imaginary unit (ß 2 = −1).Note that we have rescaled the variables x, λ, and µ by a factor β for consistency with the energy term.Conducting the summations over the variables G(i, j), n 1 (i), and n 2 (j), we get where F β (λ, µ, x) is a functional, or effective free energy, that depends on the variables λ, µ, and x and is defined by 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 (i), and n 2 (j), the expression for the effective free energy F β (λ, µ, x) depends only on N 1 + N 2 + 1 uncon-strained variables, namely λ(i), µ(j), and x.Below, we demonstrate how identifying the extremum of this function enables us to address the multi-assignment problem.

Optimizing the Effective Free Energy
Let G(i, j), n 1 (i), and n 2 (j) represent the expected values of G(i, j), n 1 (i), and n 2 (j), 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 λ, µ, and x: After some rearrangements, these two equations can be written as where and 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.
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 β (λ, µ, x) is strictly concave.Theorem 1.The Hessian of the effective free energy F β (λ, µ, x) is negative definite if p min (i) < p max (i) for all i and q min (j) < q max (j), for all j, and negative semidefinite otherwise.The free energy function F β (λ, µ, x) 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 S β be the expected state of the system at the temperature β with respect to the Gibbs distribution given in Equation ( 4).S β is associated with an expected transportation plan G β and expected indicator functions n 1 β and n 2 β .Let λ MF (i), µ MF (j), and x MF be the solutions of the system of Equation (10).Then the following identities hold: n 2β (j) = g(βµ MF (j), q min (j), q max (j)) = d MF 2 (j).
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 β, n 1β (i) and n 2β (j) 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 G β 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 G β belongs to a polytope, which we define as U k (P max min , Q max min ) (see Appendix A for a full characterization of this polytope).We can link the mean field assignment matrix G β to an optimal free energy F MF β and to an optimal internal energy U MF 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 MF Theorem 1 and Proposition 2 outlined above underscore several advantages of the proposed framework, which reframes the multi-assignment problem as a temperaturedependent 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(x) and g(x).Lastly, the convergence with respect to temperature is monotonic.

Solving Generic Assignment Problems
In the previous section, we proposed an effective free energy, F β (λ, µ, x), that depends on N 1 + N 2 + 1 unconstrained variables and on the inverse of the temperature, β, 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 β.We now relate these values to the optimal assignment energy U * : Theorem 2. Let F MF β and U 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, 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 )) and 0 < h(x) < 1, the coupling matrix at any finite temperature is fractional.We need to show that as β → +∞, the corresponding matrix G ∞ 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 ln(q max (j) − q min (j) + 1).( 14) Theorem 3. F MF β and U MF β , the mean field approximations of the exact free energy and internal energy of the system at the inverse temperature β, satisfy the following inequalities: 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 × N 2 cost matrix C possesses a unique optimal assignment matrix, denoted as G * .We have the following theorem: Theorem 4. Let G β 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 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 G ∞ is this unique solution.

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 Ḡβ that belong to U k (P max min , Q max min ) ∖ P k (P max min , Q max min ) (see Appendix A).Entries of these matrices are noninteger, 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 P k (P max min , Q max min ) when β → +∞.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 ε-perturbation of a degenerate linear programming problem reduces this problem to a non-degenerate one.

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 β is gradually increased.At each value of β, the nonlinear system of equations defined by Equation ( 10) is solved.We write this system as where A = (A λ , A µ , A x ) is a vector of predicates defined as −g(βµ(j), q min (j), q max (j)), 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 G β , the functions n 1β , n 2β , and the corresponding transportation energy U MF (β) are computed.The iterations over β are stopped if the mean field energy U MF β does not change anymore (within a tolerance TOL usually set to 10 −6 ).Otherwise, β is increased, and the current values of λ, µ 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 β value.We use for this purpose an iterative Newton-Raphson method (see [15]).
Algorithm 1 Matching: a temperature-dependent framework for solving the multiassignment problem Input: N 1 and N 2 , the number of agents and tasks; p min (i) and p max (i), the minimum and maximum number of agents needed to perform task i, and q min (j) and q max (j), 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 β 0 for β, the inverse of the temperature Initialize: Initialize arrays λ and µ to 0 and initialize x = 0. Set STEP = √ 10. for i = 1, . . .until convergence do (1) Initialize (2) Solve non linear Equation (10) for λ, µ and x at saddle point As for any annealing scheme, the initial temperature, or, in our case, the initial value β 0 , is a parameter that significantly impacts the efficiency of the algorithm.Setting β 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 β 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 β 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 β 0 C max = 1 provides satisfactory annealing efficiency.The value for STEP was chosen empirically.

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 MF β and free energy F MF β .(B) as a function of β when solving the multi-task assignment problem with the cost matrix C 1 , 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 U * as shaded areas (see text for details).
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 i.e., at each value of β, 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 β and therefore decrease monotonically as β 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.
Total cost 1440 1450 1450 a Assignment based on an ant colony algorithm [13].We ran Matching as follows.For all tasks j, we set q min (j) = q max (j) = 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 (i) 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 (i) 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. 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].
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.

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 β (the inverse of temperature) towards the optimal assignment cost.Moreover, we have established that for sufficiently large β 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.
Let A be a non-negative matrix of size N × 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 σ(A) as the sum of all elements of A, i.e., σ(A) = ∑ N i=1 ∑ M j=1 A(i, j).The transportation polytope U (R max min , C max min ) is the set of N × M matrices A that satisfy The transportation polytope U k (R max min , C max min ) is the set of N × M matrices A that satisfy as the set of all matrices in U (R max min , C max min ) whose entries are either 0 or 1 with a similar definition for P k (R max min , C max min ) with respect to U k (R max min , C max min ).Relaxed solutions to the general assignment problem considered in this paper belong to U k (P max min , Q max min ), while integer solutions belong to P k (P max min , Q max min ), 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 can be stated as any matrix in U k (R max min , C max min ) that can be written as a linear combination of matrices in P k (R max min , C max min ).This result will be useful for some of the proofs below.
where 1 N and 1 N are vectors of one of size N and M, respectively, N = M = k, U (R max min , C max min ) is the set of doubly stochastic matrices, P k (R max min , C max min ) is the set of permutation matrices, and Theorem A2 is then equivalent to the Birkhoff-von Neumann theorem for doubly stochastic matrices [17,18].
is the set of doubly sub-stochastic matrices with sum k and P k (R max min , C max min ) 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 with: where a and b are two positive integers with 0 ≤ a ≤ b.The functions H(x), 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) × (N 1 + N 2 + 1), such that its rows and columns correspond to all N 1 λ values first, followed by all N 2 µ values, and finally to the value x.Let h(x) = H ′ (x), and let h ′ be its derivative, i.e., (A1) In [15], we have shown that , that h ′ (x) is always strictly negative.Similarly, let g(x) = G ′ (x), and let g ′ be its derivative, i.e., In Appendix H, we show that when b = a, g ′ (x, a, a) = 0 ∀x ∈ R, and, when b > a, g ′ (x, a, b) is strictly positive ∀x ∈ R.
We define the matrix X ′ and the vector d ′ 1 and d ′ 2 such that From Equation ( 10), we obtain where δ are Kronecker functions, the indices i and i ′ belong to [1, N 1 ] and the indices j and j ′ belong to [1, N 2 ], and we have defined x 3 ) be an arbitrary vector of size N.The quadratic form Q(x) = x T Hx is equal to As X ′ (i, j) is based on the function h ′ that is negative and d ′ 1 (i), and d ′ 2 (j) are based on the function g ′ that is positive, the summands in the equation above are all negative for all i ∈ [1, N 1 ] and j ∈ [1, N 2 ], and therefore, Q(x) is negative for all vector x.The Hessian H is negative, semidefinite.
To check for definitiveness, let us note first that as Q(x) is a sum of negative terms, it is 0 if and only if all the terms are equal to 0. As the function h ′ (x) is strictly negative, this means that ∀(i, j), For the two other terms, we consider two cases: (i) p min (i) < p max (i) ∀i and q min (j) < q max (j) ∀j.Then g ′ (x, p min (i), p max (i)) < 0 and g ′ (x, q min (j), q max (j)) < 0. This leads to ∀(i, j): x 1 (i) 2 = 0, x 2 (j) 2 = 0, (A4) Equations (A3) and (A4) are satisfied when all x 1 (i), x 2 (j), and when x 3 are zero, namely that x = 0. Therefore in this case, H is negative, definite, and the free energy F β (λ, µ, x) is strictly concave.(ii) For all i p min (i) = p max (i) or for all j q min (j) = q max (j).Then, either d or both.There are then non-zero solutions to the equation Q(x) = 0.The Hessian is then only semidefinite.

Appendix C. Monotonicity of F MF β
The effective free energy F β (λ, µ, x) defined in Equation ( 8) is a function of the cost matrix C and of real unconstrained variables λ(i), µ(j), and x.For the sake of simplicity, for any (i, The effective free energy is then G(βµ(j), q min (j), q max (j)) H(βy(i, j)) − kx.
As written above, F β (λ, µ, x) is a function of the independent variables β, λ(i), µ(j), and x.However, under the saddle point approximation, the variables λ(i), µ(j), and x are constrained by the conditions and the free energy under those constraints is written as F MF β .In the following, we use the notations with respect to β, respectively.It can be shown easily that (see Appendix C of [15]) namely that the total derivative with respect to β is, in this specific case, equal to the corresponding partial derivative, which is easily computed to be G(βµ MF (j), q min (j), q max (j)) − 1 β µ MF (j)g(βµ MF (j), q min (j), q max (j)) H(βy MF (i)) − βy MF (i, j)h(βy MF (i, j)) .
Let t(x) = H(x) − xh(x).The function t(x) is continuous and defined over all real values x and is bounded above by 0, i.e., t(x) ≤ 0 ∀x ∈ 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 R. As l(βµ MF (j), q min (j), q max (j)) we conclude that and let the corresponding mean field approximation of the internal energy at the saddle point, Before computing dU MF β dβ , we prove the following property Proposition A1.
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, Using Proposition A1, we find that the partial derivatives of U MF β with respect to λ, µ, and x are all zeros, and therefore, Using again Proposition A1, we obtain Using Equation (A5), where t(x) and l(x, a, b) are the functions defined above, and t ′ (x) and l ′ (x, a, b) are their derivatives with respect to x.Let us define r(x) = xt ′ (x) and s(x, a, b) = xl ′ (x, a, b).Then, Note that r(x) = x 2 e x (1 + e x ) 2 .Therefore, r(x) is positive, bounded below by 0. In Appendix H, we show that s(x, a, b) is negative, bounded above 0. Therefore, The free energy is given by βµ(j), q min (j), q max (j)) H(βy(i, j)) − kx, while the internal energy is C(i, j)X(i, j).
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) − λ(i) − µ(j) − x, we obtain (after some reorganization), where t(x) and l(x, a, b) are defined above.Let us define Then, 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 where X MF β is the solution the the SPA system of equations.At a finite inverse temperature β, X MF β is strictly non-integral.In addition, X MF β satisfies constraints on row sums and column sums that make it an element of U k (R max min , C max min ); see appendix A for details.Based on Theorem A2, X MF β can be written as a linear combination of matrices in P k (R max min , C max min ), As U * is the minimum matching cost over all possible matrices in P k (R max min , C max min ), we have for all π ∈ P k (R max min , C max min ).Combining Equations (A13) and (A14), we obtain from which we conclude that at each β, Therefore U * ≤ U MF (∞), and consequently U * ≤ F MF (∞), based on Equation (A12).
Appendix E.4.U * ≥ F MF (∞) Proof.Let us first recall the definition of the free energy, Note these two properties of limits: when a and b are non-negative integers.Therefore, Let us consider now a matrix π in P k (R max min , C max min ) (see above).We can write where RS and CS are the row sums and column sums of π, 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: ∑ i,j y(i, j)π(i, j) ≥ ∑ i|y(i,j)≤0 y(i, j).(A17) As π belongs to P k (R max min , C max min ), its row sums satisfy the constraints: p min (i) ≤ RS(i) ≤ p max (i). (A18) We multiply this equation separately for positive and negative λ(i).We find ∑ i|λ(i)>0 p min (i)λ(i) ≤ ∑ i|λ(i)>0 λ(i)RS(i) ≤ ∑ i|λ(i)>0 p max (i)λ(i) p min (i)λ(i). Therefore, p max (i)λ(i) − ∑ i|λ(i)<0 p min (i)λ(i).q max (j)µ(j) − ∑ j|µ(j)<0 q min (j)µ(j). (A20) Combining Equations (A15)-(A17), (A19) and (A20), we obtain ∑ i,j C(i, j)π(i, j) ≥ lim β→+∞ F β (λ, µ, x).(A21) Equation (A21) is valid for all matrices π in P k (R max min , C max min ).It is therefore valid for the optimal π * that solves the general assignment problem.Since U * = ∑ i,j C(i, j)π * (i, j), we have As this equation is true for all λ, µ, and x, it is true in particular for λ = λ MF , µ = µ MF , and x = x MF , leading to We have shown that U * ≤ F MF (∞) and F MF (∞) ≤ U * , and therefore U * = F MF (∞).The corresponding result for the internal energy, U * = U MF (∞) follows directly from Equation (A12).a π = 1.
We assume that the general assignment problem considered has a unique solution.We want to prove that max i,j G β (i, j) − G * (i, j) ≤ A(N 1 ,N 2 ) β∆ , where G * is the optimal solution of the assignment problem, ∆ = 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 A(N 1 , N 2 ) β∆ < Ḡβ (i, j) − G * (i, j) .
In the first case, a π π(i, j).In conclusion, in both cases, we have In Theorem 3, we have shown that Therefore, as ∆ is strictly positive.We have shown that 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 ′ (x, a, a) = 0, l(x, a, a) = 0, and s(x, a, a) = 0.In the following, we only consider a < b.
and monotonic decreasing functions, respectively, of the parameter β.Proof.See Appendix C for F MF β and Appendix D for U MF β .

( 4 )
and U i,MFβ Check for convergence: if |U i,MF β − U i−1,MF β | < TOL, stop end for Output:The converged assignment matrix G β , the functions n 1β , n 2β over the agents and tasks, and the minimal associated cost U β .

Figure 1 .
Figure 1.Convergence of the internal energy U MF β (A) and free energy F MF β to differentiate between the total derivative and partial derivative of F MF β with all a π ∈ [0, 1] and∑ π∈P k (R max min ,C max min )

Table 1 .
The cost matrix C 1 .

Table 2 .
Assignments for C 1 with p max = 4.

Table 3 .
The cost matrix C 2 .

Table 4 .
Solving the multi-task assignment problems defined by the cost matrix C 1 .

Table 5 .
Solving the multi-task assignment problems defined by the cost matrix C 2 .