Solving Sparse Assignment Problems on FPGAs

The assignment problem can be solved with well-studied algorithms from linear programming like the simplex method [33]. The simplex method starts by observing that the solution space of a linear programming problem is a polytope defined by the constraints of the linear program. Since the cost function is linear it can be shown that the optimum must lie on a vertex of this polytope. For the assignment problem, each possible mapping of agents to objects make up a vertex in a \((n+m)\)-dimensional polytope, where \(n\) is the number of agents and \(m\) is the number of objects. The Simplex method starts by picking an initial vertex. It then traverses the edges, from vertex to vertex, until the optimum is found. In the worst case, the algorithm must visit all vertices, giving it a computational complexity equivalent of a brute force method with \(O(2^N)\).

The assignment problem belongs to a special class of linear programming problems known as transportation problems. Transportation problems can be represented as graphs where the decision variables are the weights of the edges. Each decision variable only occurs twice in the constraints, once for the source node and once for the destination node. This additional structure is taken advantage in, e.g., the Streamlined Simplex [22].

In recent years, enterprise optimization software such as Gurobi [15] and CPLEX [6] have gained popularity. Here, the assignment problem is solved as a linear integer programming problem using Simplex-like algorithms.

The Hungarian method was introduced by Khun as a more efficient solution to the assignment problem than the traditional linear programming methods. In the Hungarian method, the assignment problem is modeled as a graph \(G=(V,E)\) and a reward matrix \(R\), where \(V = A \cup O\) is the set of vertices, \(E=\lbrace (a_i, o_j)\rbrace _{a_i \in A, o_j \in O}\) is the set of edges, and \(r_{ij}\) is the weight of edge \((a_i, o_j)\). A matching \(M\) is a subset of \(E\) such that no edges in \(M\) share a common vertex. We call \(M\) a perfect matching if it matches all the vertices of \(V\). Given a graph \(G\) and a matching \(M\), an augmenting path is a path in \(G\) that starts and ends at an unmatched vertex and contains alternating edges from \(M\) and \(E-M\). Augmenting paths can be used to iteratively find perfect matchings.

The assignment problem can thus be reformulated as finding the perfect matching with the maximum weight. Compared to the general class of transportation problems, the assignment problem has two additional features. First, the weights are binary, i.e., either 1 or 0. Second, if a particular edge between agent \(a_i\) and object \(o_j\) has the weight 1, then all other edges connected to both \(a_i\) and \(o_j\) has the weight 0. This information is not utilized by the linear programming methods and makes them comparatively slow [5]. The Hungarian method solves the assignment problem in polynomial with a computational complexity of \(O(N^3)\). The Hungarian method is a primal-dual algorithm where the starting point is a feasible solution. It works by maintaining and iteratively updating two data structures. The first is a labeling \(l\) of \(G\). A labeling is a function mapping the set of vertices to the real numbers. A valid labeling is defined as one where the the sum of the labels for two vertices is less than or equal to the weight of the edge connecting them. The second is a matching \(M\). It can be shown that given a labeling \(l\) and the equality subgraph \(G_l\), then a perfect matching \(M\) of \(G_l\) is also the maximum weight matching of \(G\). The Hungarian method starts with a matching \(M\) and a valid labeling \(l\). It then augments the matching until either it finds a perfect matching or no augmenting path exists. In the latter case, it improves the labeling function \(l\) before going back to looking for an augmenting path.

The Hungarian method is not suitable for parallelization as both the labeling improvement and finding an augmenting path are inherently sequential operations.

The LAPJV algorithm was introduced by Jonker and Volgenant [21] as an improvement of the Hungarian method. While the Hungarian method works by finding any feasible augmenting path, LAPJV looks for the shortest augmenting path. This is found using Dijkstra’s shortest path algorithm.

The Auction algorithm was proposed by Bertsekas [4] as a general solution to the assignment problem. The assignment problem is modeled as an auction and the reward \(r_{ij}\) is interpreted as how much agent \(a_i\) is willing to pay for the object \(o_j\).

The Auction method is not a primal-dual algorithm like the Hungarian method, and each iteration is not guaranteed to reduce the total cost. Optimality can be guaranteed if \(\epsilon\), the minimum raise, is less than a threshold. The algorithm is guaranteed to terminate if there exists a feasible solution. The relaxed constraint on monotonically improving solutions enables the Auction method to reach the optimum in fewer steps. Its computational complexity is \(O(N^2log(N))\) [29].

The Auction algorithm works by updating a price vector \(P\), which contains the current price for each object, and an assignment vector that stores information about which agents are assigned to which objects. The price vector is initialized to 0 and is increased as agents bid for objects. Each iteration consists of two phases: the bidding phase and the assignment phase. In the bidding phase, a subset of the agents \(A_{bid} \subseteq A\) calculate bids for a subset of the objects \(O_{bid} \subseteq O\). In the assignment phase, the highest bid is picked and the bidding agent is assigned to the object. There are two commonly used approaches to the bidding phase. The first alternative is referred to as the Jacobi variant. Here, all unassigned agents bid for a single unassigned object; i.e., \(A_{bid} = A_{unassigned}\) and \(O_{bid} = o_k\), where \(A_{unassigned}\) is the set of unassigned agents and \(o_k\) is a single unassigned object. Referring to the example in Table 1, the first iteration of the Jacobi method would consist of Alice, Bob, and Charlie bidding of object \(o_1\). The second alternative is referred to as the Gauss-Seidel variant. Here, a single unassigned agent bids for all the unassigned objects. This means that \(A_{bid} = a_k\) and \(O_{bid} = O\). Here, \(a_k\) is a single unassigned agent. The first iteration of the Gauss-Seidel variant would consist of Alice bidding for objects \(o_1,o_2,o_3\). The Jacobi variant can be implemented by assigning a processing element to each agent to compute the bids in parallel. The disadvantage of this approach is that as more agents are assigned, fewer of the processing elements are used. The Gauss-Seidel variant can be implemented by assigning a processing element to each object. It does not suffer from the same disadvantage as the Jacobi method, as the number of objects being bid on remains constant as the algorithm proceeds. For the rest of this work, the Gauss-Seidel variant is used.

The winning bid each round is the bid for the object that has the highest benefit to the agent. The benefit of object \(j\) for agent \(i\) is the difference between the agents’ reward \(R(i,j)\) for the object \(o_j\) and the current price \(P(j)\), computed as:

When the object with the highest benefit is found, the agent decides its bid for that object. The bid is made such that the object still has the highest benefit to the agent. This is achieved by raising the bid for the object with an amount equal to the difference between the highest and the second-highest benefit, expressed as:

The pseudo-code for the Auction algorithm is presented in Figure 1. Lines 1–3 initialize the arrays containing object prices and assignments and the queue containing the unassigned agents. Inside the loop, an unassigned agent is dequeued and agents’ benefits for each object are computed. Line 7 extracts the index of the object with the highest benefit, while Line 8 gets the highest and second-highest benefit. Line 9 calculates the raise according to Equation (5). Lines 14–17 unassign any previously assigned agent. This proceeds until agents are either assigned or disregarded due to not being able to create a positive raise (lines 10–13). An optimized C++ implementation is published in Reference [19].

Figure 3 shows the execution of the four iterations needed to solve the assignment problem in Table 1 using the Auction algorithm.

There are multiple ways of representing sparse matrices [30]. A simple tuple representation [28] was chosen, where the element and the column index are stored as a tuple. The original reward matrix \(R\) resides in the off-chip DRAM and each element of the reward matrix is represented by a 64-bit data word. The transformed vector \(C\) is the sparse representation of \(R\) and resides in the on-chip BRAM. For architecture with \(n\) PEs, each BRAM data word in \(C\) contains \(n\) tuples of object reward and column index as given:

The number of bits used to represent a single object reward is configurable. The words of the transformed BRAM vector \(C\) are padded such that the first non-zero element of each row of the original DRAM matrix \(R\) always occupies the first element \(e_1\) of a BRAM word \(c\). This means that each BRAM word only contains rewards from the same row of the original matrix. A row of the original matrix may span several words of the transformed vector \(C\). A separate array is used to map the row indices of the original matrix to indices of the transformed matrix.

In Figure 6, an example of a transformation between a problem matrix \(R\) and its sparse representation \(C\) is given. An architecture with 4 PEs is assumed for simplicity. This means that each BRAM word stores up to four valid, non-zero object rewards and their column indices. The reward matrix \(R\) is sparse with the majority of the rewards being zero. The bottom right corner has indices (0,0). Each element \(c_i\) in the transformed vector \(C\) consists of four tuples \(e_1 \ldots e_4\). The individual tuples consist of the element value and its column index in the original matrix \(R\). The transformed vector \(C\) is not fully dense, because the first element of each row in \(R\) is aligned to the beginning of a data word in \(C\). This is done to efficiently read out an entire row of the original input matrix by reading a single BRAM data word. In the example, elements 1 and 2 of \(C\) (zero-indexed starting from the bottom) both belong to row 1 of \(R\). When the number of valid elements in a row in \(R\) exceeds the number of PEs in our architecture the elements are spread out over multiple words. The row index of each element in \(C\) is stored in an array that maps each row of \(R\) to an element in \(C\).

There are three major costs associated with using a sparse representation. First, the translation from dense to sparse representation must be performed. In the SaS Auction architecture, the translation is done in hardware, and this incurs an area cost as opposed to the time cost associated with a software implementation. Second, a sparse representation removes the static relationship between PE and object price. In SoA Auction, PE0 always calculates the benefit of object 0 for the various agents. This means that each PE can store the prices locally and thus can all the prices be accessed in parallel at low cost. To achieve parallel lookup of all PEs, SaS Auction duplicates the BRAM with object prices for each PE. Last, sparse representation adds overhead, because two values are needed for each element: the column index and the element value. The memory needed to store a sparse matrix is given by:where \(w_{reward}\) is the bits needed to represent a reward, \(n_{max}\) is the maximum problem size, and \(r_{sparse}\) is the rate of non-zero elements. The memory requirements for a dense representation is as follows:

For an FPGA implementation, both \(Mem_{sparse}\), \(n_{max}\) and \(w_{reward}\) are fixed at compile-time. For a given combination of these parameters, there is an upper bound on the size and density of a problem matrix that can be stored.

SoA Auction contains a simple memory interface that either can read the rewards from an off-chip main memory or an on chip BRAM. The latter necessitates that the SW has transferred the data to BRAM first. SaS Auction implements a DMA in Dram2Bram, which moves the data from off-chip DRAM to on-chip BRAM and compresses the data on the fly.

With sparsity-aware processing, the static mapping of object price to PE no longer exists. A object price might be needed by any PE. To enable single-cycle price-lookup for all PEs, SaS Auction stores a copy of the prices for each PE.

SaS Auction implements speculative processing. While SoA Auction stalls execution in the PE stage, SaS Auction continues with possibly stale prices. The AssignmentEngine will detect a strong data dependency and reschedule the agent that used stale prices.

First, a cost model for the state-of-the-art implementation is derived. The cost of stalling execution at the PE stage depends on the number of clock cycles needed by the preceding agent to calculate and commit its bid. It can be expressed as:where \(d_\mathrm{PE}\), \(d_\mathrm{SearchTask}\) and \(d_\mathrm{AssignmentEngine}\) refer to the latency (or depth) of the respective modules, given in Table 2. Note that these values are based on our implementation inspired by Zhu et al. [36], where no data regarding the latencies of the original implementation is available.

Clearly, the stalling cost is independent of the input problem dimension and only depends on the number of PEs. However, the total processing latency of one agent also depends on the size of input problem dimension relative to the number of PEs. Thus, the relative cost of stalling decreases the more underdimensioned the accelerator is.

The speculation cost is defined as the additional number of clock cycles needed to process an agent if it is not stalled at the PE stage, but rather being speculatively processed. This cost is a stochastic variable, as it depends on whether there is a data dependency between the agent and any unresolved agents. The goal is thus to find the expected value of the speculative cost. The misspeculation cost is first derived. It is the penalty of dispatching an agent from the PE stage if there is an unresolved strong dependency and is defined as:

The constant cost of \(C_{stall}\) is the worst-case expected cost of redoing the loop consisting of PE, SearchTask, and AssignmentEngine stages once the agent is ready in the PE-stage again. The variable \(n_\mathrm{bubbles}\) is the number of pipeline bubbles occurring between the misspeculation and the next instant the agent is ready in the PE-stage. It is estimated as:

If \(n_\mathrm{unassigned}\) is larger than the number of clock cycles from the AssignmentEngine to the PE stage, then there are no bubbles. The expected costs of stalling and speculating are expressed as:where \(p_\mathrm{miss}\) is the probability of a misspeculation, i.e., the probability that the agent at the PE stage has a strong dependency on an unresolved agent.

To estimate the probability of a strong data dependencies, \(p_{miss}\), some simplifications must be made. It is assumed that the most beneficial object \(o_{pick}\) for an agent \(a_k\) is an independent random variable drawn from a uniform distribution. This can be expressed as follows:where \(m\) is the number of objects.

The probability of a strong dependency can thus be expressed as the probability of an agent picking the same object as any of the unresolved agents.where \(n\) is the number of unresolved agents and \(m\) is the number of objects. The actual number of unresolved agents \(n\) depends several factors.

First, it has an upper bound defined by the number of PEs in the accelerator. The number of PEs decides the number of pipeline stages between reading the prices in the PE stage and committing a bid in the AssignmentEngine. In an architecture with 8 PEs, there are 7 pipeline stages between PE and the commit operation.

Second, it also depends on how many pipeline stages each agent spans. In a scenario with 8 PEs and agents with 64 valid object rewards, the agent will span 8 pipeline stages. This means that when a new agent arrives at the PE stage, there will be at most one unresolved agent that would span all the 7 pipeline stages between the PE stage and the commit stage. Consider another scenario, also with 8 PEs, but this time with agents with at most 8 valid object rewards. This implies that each agent only spans a single pipeline stage. In this scenario, when a new agent arrives at the PE stage, there might be as many as 7 unresolved agents.

Last, it also depends on how many unassigned agents are left. When the agent of the previous scenario arrives at the PE stage. The number of unresolved agents could be anywhere from 0–7. The actual number depends on how full the pipeline is, which depends on how many unassigned agents are left.

Finally, using the derived probability of misspeculation \(p_{miss}\), Equation (11) is plotted in Figure 10 for an accelerator with 8 PEs with different input problem sizes. The grey horizontal line represents the cost of stalling, which is constant and independent of the number of unresolved agents and the number of potential pipeline bubbles. The expected cost of speculating depends on both \(n_\mathrm{unresolved}\) and the potential number of bubbles \(n_\mathrm{bubbles}\). Thus, Figure 10 shows the range of expected cost of speculating. The maximum expected cost occurs in the case of an empty pipeline behind the agent in the PE stage. The minimum cost arises when the pipeline is full and there are enough unassigned agents such that \(n_\mathrm{bubbles}=0\).

Clearly, for bigger problems the probability of misspeculation is very small and thus also the expected cost of speculation. Intuitively, increasing number of unresolved agents leads to increased chance of misspeculation and thus increased expected cost. Given the assumption of uniform probability distribution, a static speculation scheme is sufficient. In static speculation, as opposed to dynamic speculation, the pipeline always proceeds in the PE stage, regardless of the number of unresolved agents or the number of expected pipeline bubbles. Such a static speculation scheme is implemented by SaS Auction.

Figure 12 shows the results of evaluating SaS and SoA Auction on the random sparse dataset. Each plot is divided into six subplots. Each subplot shows problems with the same input dimensions, which are given in the subplot title. The y-axes denote the processing time in nanoseconds, while the x-axes show the number of valid object-rewards per agent, i.e., the sparsity level. There are five sparsity levels in each subplot; e.g., in the 600 \(\times\) 600 subplot there are results for problems with 6, 12, 60, 120, and 180 valid object rewards per agent. For each sparsity level, the reported value is the average processing time for that particular sparsity level.

In Figure 12(a) the core processing times of SaS and SoA Auction with different number of PEs is shown. The core processing time is the time spent doing actual computation, neglecting the time spent transferring the problem matrices from DRAM to the on-chip memory. The core processing time is measured by the accelerator itself, which counts the cycles spent computing the solution. There are several things to note.

As expected, the potential speedup when taking advantage of speculation and sparsity-aware processing is big, as high as \(50\times\) for big sparse problems. It is also clear that the optimal number of PEs is related to the sparsity of the problem. The 50 \(\times\) 50 subplot shows that SaS 8PE performs best while there are less than 8 valid objects per agent. Then SaS 16PE performs better until we have more than 16 valid objects per agent when SaS 32PE is best.

It is also clear that the speedup of SaS relative SoA is generally higher for sparse than for dense problems. For large and sparse input problems the potentially speedup is enormous. This is mostly due to the sparsity-aware processing; e.g., for a 600 \(\times\) 600 problem each agent in SoA 8PE will span 75 pipeline stages, regardless of the sparsity. When the problem is fully dense, the only remaining speedup is due to the dependency speculation.

Figure 12(b) depicts the total processing times of SaS and SoA Auction. The total processing time includes the time spent transferring the rewards from DRAM to the on-chip-memory and the time spent writing the final results back. Clearly, the relative speedup of the SaS accelerators is limited by the fact that a major portion of the runtime is spent moving data from the DRAM to the on-chip BRAM. This effect is most visible for highly sparse problems. Consider the 600 \(\times\) 600 problems with only 6 valid object rewards per agent. The core processing time of SoA 8PE is 50\(\times\) higher than that of SaS 8PE. However, the total processing time is only 1.8\(\times\) higher. The reason is that these highly sparse problems require fewer iterations to solve than dense problems. This results in a situation where the core processing time only accounts for 2% of the total processing time.

Another interesting observation is that neither the core, nor the total, processing time is monotonically increasing as the problems get denser. Rather, the processing time peaks at a relatively low density. Figure 11 shows the number of iterations needed by the C++ (SW) implementation to solve a dataset of random sparse 50 \(\times\) 50 problems. The x-axis shows the level of density and the y-axis shows the number of iterations. Surprisingly, the processing time peaks at 22% density with an average core processing time of 7,173 ns before monotonically decreasing as density increases. At 100% density the core processing time has decreased to an average of 4,964 ns. The mechanism behind this is not clear to the authors.

Figure 12(c) compares the total processing time of SaS and SoA Auction with SW implementations. The y-axis now uses a logarithmic scale. As described in prior work, there is enormous potential for speedup: from 100\(\times\) speedup for the small 10 \(\times\) 10 problems to 1,000\(\times\) speedup for the sparse 600 \(\times\) 600 problems. It is also clear that a vanilla C++ auction algorithm implementation outperforms industrial optimization software based on linear programming. However, for larger problems the speedup is less significant.

Figure 13 shows the results of evaluating SaS and SoA Auction on the MOT dataset. There are a few subtleties that make this plot different than the previous ones. Each subplot now shows problems with the same number of rows (i.e., agents). The subplots have titles like 16xY, which means that the problems have 16 agents, however, the number of objects vary from 1 to over 300. The x-axes show the number of valid object-rewards per agent; this is no longer directly tied to the sparsity level, since the number of non-valid object-rewards varies. The values on the x-axes now indicate limits; e.g., for the 16xY subplot, the first bar-group is the average processing time for problems with 1–8 valid object-rewards, the second is for 8–20 valid object-rewards, and so on.

As expected, the processing time of SaS Auction increases as the number of valid object-rewards increases. This is in line with the results from the randomly sparse dataset. Interestingly, this is not the case for SoA Auction, where the processing time does not seem to strongly correlate with the number of valid object-rewards. The reason is that the processing time of SoA Auction mostly only depend on the input dimensions of the problem. The input distribution of the MOT problems are such that the agents seldom bid for the same objects. Thus, the sparsity level is not an indicator of how many iterations are needed and thus to overall processing time. The SaS Auction, however, can take advantage of the sparsity and drastically reduce the core processing time.