Differentiable Equilibrium Computation with Decision Diagrams for Stackelberg Models of Combinatorial Congestion Games
Year: 2021

Abstract:
We address Stackelberg models of combinatorial congestion games (CCGs); we aim to optimize the parameters of CCGs so that the selfish behavior of non-atomic players attains desirable equilibria. This model is essential for designing such social infrastructures as traffic and communication networks. Nevertheless, computational approaches to the model have not been thoroughly studied due to two difficulties: (I) bilevel-programming structures and (II) the combinatorial nature of CCGs. We tackle them by carefully combining (I) the idea of differentiable optimization and (II) data structures called zero-suppressed binary decision diagrams (ZDDs), which can compactly represent sets of combinatorial strategies. Our algorithm numerically approximates the equilibria of CCGs, which we can differentiate with respect to parameters of CCGs by automatic differentiation. With the resulting derivatives, we can apply gradient-based methods to Stackelberg models of CCGs. Our method is tailored to induce Nesterov’s acceleration and can fully utilize the empirical compactness of ZDDs. These technical advantages enable us to deal with CCGs with a vast number of combinatorial strategies. Experiments on real-world network design instances demonstrate the practicality of our method.

text: We address Stackelberg models of combinatorial congestion games (CCGs); we aim to optimize the parameters of CCGs so that the selfish behavior of non-atomic players attains desirable equilibria. This model is essential for designing such social infrastructures as traffic and communication networks. Nevertheless, computational approaches to the model have not been thoroughly studied due to two difficulties: (I) bilevel-programming structures and (II) the combinatorial nature of CCGs. We tackle them by carefully combining (I) the idea of differentiable optimization and (II) data structures called zero-suppressed binary decision diagrams (ZDDs), which can compactly represent sets of combinatorial strategies. Our algorithm numerically approximates the equilibria of CCGs, which we can differentiate with respect to parameters of CCGs by automatic differentiation. With the resulting derivatives, we can apply gradient-based methods to Stackelberg models of CCGs. Our method is tailored to induce Nesterov’s acceleration and can fully utilize the empirical compactness of ZDDs. These technical advantages enable us to deal with CCGs with a vast number of combinatorial strategies. Experiments on real-world network design instances demonstrate the practicality of our method.

heading: 1 Introduction
text: Congestion games (CGs) [49] form an important class of non-cooperative games and appear in various resource allocation scenarios. Combinatorial CGs (CCGs) can model more complex situations where each strategy is a combination of resources. A well-known example of a CCG is selfish routing [50], where each player on a traffic network chooses an origin-destination path, which is a strategy given by a combination of some roads with limited width (resources). Computing the outcomes of players’ selfish behaviors (or equilibria) is essential when designing social infrastructures such as traffic networks. Therefore, how to compute equilibria of CCGs has been widely studied [5, 12, 57, 42].
In this paper, we are interested in the perspective of the leader who designs non-atomic CCGs. For example, the leader aims to optimize some traffic-network parameters (e.g., road width values) so that players can spend less traveling time at equilibrium. An equilibrium of non-atomic CCGs is characterized by an optimum of potential function minimization [41, 53]. Thus, designing CCGs can be seen as a Stackelberg game; the leader optimizes the parameters of CCGs to minimize an objective function (typically, the social-cost function) while the follower, who represents the population of selfish non-atomic players, minimizes the potential function. This mathematical formulation is called the Stackelberg model in the context of traffic management [46]. Therefore, we call our model with general combinatorial strategies a Stackelberg model of CCGs.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Stackelberg models of CCGs have been studied for cases where the potential function minimization has desirable properties. For example, Patriksson and Rockafellar [46] proposed a descent algorithm for traffic management using the fact that projections onto flow polyhedra can be done efficiently. However, many practical CCGs have more complicated structures. For example, in communication network design, each strategy is given by a Steiner tree (see [24, 42] and Section 1.1), and thus the projection (and even optimizing linear functions) is NP-hard. How to address such computationally challenging Stackelberg models of CCGs has not been well studied, despite its practical importance.
Inspired by a recent equilibrium computation method [42], we tackle the combinatorial nature of CCGs by representing their strategy sets with zero-suppressed binary decision diagrams (ZDDs) [39, 32], which are well-established data structures that provide empirically compact representations of combinatorial objects (e.g., Steiner trees and Hamiltonian paths). Although the previous method [42] can efficiently approximate equilibria with a Frank–Wolfe-style algorithm [18], its computation procedures break the differentiability of the outputs in the CCG parameters (the leader’s variables), preventing us from obtaining gradient information required for optimizing leader’s objective functions.
Our contribution is to develop a differentiable pipeline from leader’s variables to equilibria of CCGs, thus enabling application of gradient-based methods to the Stackelberg models of CCGs. We smooth the Frank–Wolfe iterations using softmin, thereby making computed equilibria differentiable with respect to the leader’s variables by automatic differentiation (or backpropagation). Although the idea of smoothing with softmin is prevalent [29, 38], our method has the following technical novelty:
• Our algorithm is tailored to induce Nesterov’s acceleration, making both equilibrium computation and backpropagation more efficient. To the best of our knowledge, the idea of simultaneously making iterative optimization methods both differentiable and faster is new.
• Our method consists of simple arithmetic operations performed with ZDDs as in Algorithm 2. This is essential for making our equilibrium computation accept automatic differentiation. The per-iteration complexity of our method is linear in the ZDD size.
Armed with these advantages, our method can work with CCGs that have an enormous number of combinatorial strategies. We experimentally demonstrate its practical usefulness in real-world network design instances. Our method brings benefits by improving the designs of social infrastructures.
Notation. Let [n] := {1, . . . , n}. For any S ⊆ [n], 1S ∈ {0, 1}n denotes a binary vector whose i-th entry is 1 if and only if i ∈ S. Let ∥ · ∥ be the ℓ2-norm.

heading: 1.1 Problem setting
text: We introduce the problem setting and some assumptions. For simplicity, we describe the symmetric setting, although our method can be extended to an asymmetric setting, as explained in Appendix A.
Combinatorial congestion games (CCGs). Suppose that there is an infinite amount of players with an infinitesimal mass (i.e., non-atomic). We assume the total mass is 1 without loss of generality. Let [n] be a set of resources and let S ⊆ 2[n] be a set of all feasible strategies. We define d := |S|, which is generally exponential in n. Each player selects strategy S ∈ S. Let y ∈ [0, 1]n be a vector whose i-th entry indicates the total mass of players using resource i ∈ [n]. In other words, if we let z ∈ △d be a vector whose entry zS (S ∈ S) indicates the total mass of players choosing S, we have y := ∑ S∈S zS1S ∈ Rn. Therefore, y is included in convex hull C := {∑ S∈S zS1S
∣∣ z ∈ △d}, where △d := { z ∈ Rd
∣∣ z ≥ 0,∑S∈S zS = 1} is the (d− 1)-dimensional probability simplex. A player choosing S incurs cost cS(y) := ∑ i∈S ci(yi), where each ci : R → R is assumed to be strictly increasing; this corresponds to a natural situation where cost ci increases as i ∈ [n] becomes more congested. Each player selfishly selects a strategy to minimize his/her own cost.
Equilibrium and potential functions. We say z ∈ △d attains a (Wardrop) equilibrium if for every S ∈ S such that zS > 0, it holds that cS(y) ≤ minS′∈S cS′(y), where y := ∑ S∈S zS1S . That is, no one has an incentive to deviate unilaterally. Let f : Rn → R be a potential function defined as f(y) := ∑ i∈[n] ∫ yi 0 ci(u)du. From the first-order optimality condition, it holds that z attains an
equilibrium iff y = ∑
S∈S zS1S satisfies y = argminu∈C f(u). Note that minimizer y is unique since ci is strictly increasing, which means f is strictly convex (not necessarily strongly convex).
Stackelberg model of CCGs. We turn to the problem of designing CCGs. For i ∈ [n], let ci(yi;θ) be a cost function with parameters θ ∈ Θ. We assume ci(yi;θ) to be strictly increasing in yi for any θ ∈ Θ and differentiable with respect to θ for any y ∈ C. Let f(y;θ) =∑i∈[n] ∫ yi0 ci(u;θ)du be a parameterized potential function, which is strictly convex in y for any θ ∈ Θ. A leader who designs CCGs aims to optimize θ values so that an objective function, F : Θ× C → R, is minimized at an equilibrium of CCGs. Typically, F is a social-cost function defined as F (θ,y) = ∑ i∈[n] ci(yi;θ)yi, which represents the total cost incurred by all players. Since an equilibrium is characterized as a minimizer of potential function f , the leader’s problem can be written as follows:
minimize θ∈Θ F (θ,y) subject to y = argmin u∈C f(u;θ). (1)
Since minimizer y(θ) := y is unique, we can regard F (θ,y(θ)) as a function of θ. We study how to approximate derivatives of F (θ,y(θ)) with respect to θ for applying gradient-based methods to (1).
Example 1: traffic management. We are given a network with an origin-destination (OD) pair. Let [n] be the edge set and let S ⊆ 2[n] be the set of all OD paths. Each edge in the network has cost function ci(yi;θ), where θ controls the width of the roads (edges). A natural example of the cost functions is c(yi;θ) = yi/θi for θi > 0 (see, e.g., [46]), which satisfies the above assumptions, i.e., strictly increasing in yi and differentiable in θi. Once θ is fixed, players selfishly choose OD paths and consequently reach an equilibrium. The leader wants to find θ that minimizes social cost F at equilibrium, which can be formulated as a Stackelberg model of form (1). Note that although the Stackelberg model of standard selfish routing is well studied [46], there are various variants (e.g., routing with budget constraints [28, 42]) for which existing methods do not work efficiently.
Example 2: communication network design. We consider a situation where multi-site meetings are held on a communication network (see, e.g., [24, 42]). Given an undirected network with edge set [n] and some vertices called terminals, groups of people at terminals hold multi-site meetings, including people at all the terminals. Since each group wants to minimize the communication delays caused by congestion, each selfishly chooses a way to connect all the terminals, namely, a Steiner tree covering all the terminals. If we let ci(yi;θ) indicate the delay of the i-th edge, a group choosing Steiner tree S ∈ S incurs cost cS(y;θ). As with the above traffic-management example, the problem of optimizing θ to minimize the total delay at equilibrium can be written as (1).

heading: 1.2 Related work
text: Problems of form (1) arise in many fields, e.g., Stackelberg games [55], mathematical programming with equilibrium constraints [36], and bilevel programming [10, 13], which have been gaining attention in machine learning [17, 15]. Optimization problems with bilevel structures are NP-hard in most cases [22] (tractable cases include, e.g., when follower’s problems are unconstrained and strongly convex [19], which does not hold in our case). Thus, how to apply gradient-based methods occupies central interest [15, 20]. In our Stackelberg model of CCGs, in addition to the bilevel structure, the follower’s problem is defined on combinatorial strategy sets S , further complicating it. Therefore, unlike the above studies, we focus on how to address such difficult problems by leveraging computational tools, including ZDDs [39] and automatic differentiation [34, 21].
Stackelberg models often arise in traffic management. Although many existing studies [46, 35, 6, 9] analyze theoretical aspects utilizing instance-specific structures (e.g., compact representations of flow polyhedra), applications to other types of realistic CCGs remain unexplored. By contrast, as with the previous method [42], our method is built on versatile ZDD representations of strategy sets, and the derivatives with respect to CCG parameters can be automatically computed with backpropagation. Thus, compared to the methods studied in traffic management, ours can be easily applied to and works efficiently with a broad class of realistic CCGs with complicated combinatorial strategies.
Our method is inspired by the emerging line of work on differentiable optimization [4, 60, 2, 51]. For differentiating outputs with respect to the parameters of optimization problems, two major approaches have been studied [20]: implicit and iterative differentiation. The first approach applies the implicit function theorem to equation systems derived from the Karush–Kuhn–Tucker (KKT) condition (akin to the single-level reformulation approach to bilevel programming). In our case, this approach is too expensive since the combinatorial nature of CCGs generally makes the KKT equation system exponentially large [16]. Our method is categorized into the second approach, which computes
a numerical approximation of an optimum with iterations of differentiable steps. This idea has yielded success in many fields [14, 37, 45, 7, 17]. Concerning combinatorial optimization, although differentiable methods for linear objectives are well studied [38, 60, 47, 8], no differentiable method has been developed for convex minimization on polytopes of, e.g., Steiner trees or Hamiltonian paths; this is what we need for dealing with the potential function minimization of CCGs. To this end, we use a Frank–Wolfe-style algorithm and ZDD representations of combinatorial objects.

heading: 2 Differentiable iterative equilibrium computation
text: We consider applying gradient-based methods (e.g., projected gradient descent) to problem (1). To this end, we need to compute the following gradient with respect to θ ∈ Θ ⊆ Rk in each iteration:
∇F (θ,y(θ)) = ∇θF (θ,y(θ)) +∇y(θ)⊤∇yF (θ,y(θ)),
where ∇θF (θ,y(θ)) and ∇yF (θ,y(θ)) denote the gradients with respect to the first and second arguments, respectively, and ∇y(θ) is the n× k Jacobian matrix.1 The computation of ∇y(θ) is the most challenging part and requires differentiating y(θ) = argminu∈C f(u;θ) with respect to θ. We employ the iterative differentiation approach for efficiently approximating ∇y(θ).

heading: 2.1 Technical overview
text: For computing equilibrium y(θ), Nakamura et al. [42] solved potential function minimization with a variant of the Frank–Wolfe algorithm [33], whose iterations can be performed efficiently by using compact ZDD representations of combinatorial strategies. To the best of our knowledge, no other equilibrium computation methods can deal with various CCGs that have complicated combinatorial strategies, e.g., Steiner trees. Hence we build on [42] and extend their method to Stackelberg models.
First, we review the standard Frank–Wolfe algorithm. Starting from x0 ∈ C, it alternately computes st = argmins∈C⟨∇f(xt;θ), s⟩ and xt+1 = (1 − γt)xt + γtst, where γt is conventionally set to 2 t+2 . As shown in [18, 27], xT has an objective error of O(1/T ). Thus, we can obtain numerical approximation yT (θ) = xT of equilibrium y(θ) such that f(yT (θ);θ) − f(y(θ);θ) ≤ O(1/T ). For obtaining gradient ∇yT (θ), however, the above Frank–Wolfe algorithm does not work (neither does its faster variant used in [42]). This is because st = argmins∈C⟨∇f(xt;θ), s⟩ is piecewise constant in θ, which makes ∇yT (θ) zero almost everywhere and undefined at some θ. To resolve this issue, we develop a differentiable Frank–Wolfe algorithm by using softmin. We denote the softmin operation by µS(c) (detailed below). Since softmin can be seen as a differentiable proxy for argmin, one may simply replace st = argmins∈C⟨∇f(xt;θ), s⟩ with st = µS(ηt∇f(xt;θ)), where ηt > 0 is a scaling factor. Actually, the modified algorithm yields an O(1/T ) convergence by setting ηt = Ω(t) (see [27, Theorem 1]). This modification, however, often degrades the empirical convergence of the Frank–Wolfe algorithm, as demonstrated in Section 4.1. We, therefore, consider leveraging softmin for acceleration while keeping the iterations differentiable. Based on an accelerated Frank–Wolfe algorithm [59], we compute yT (θ) as in Algorithm 1, and obtain ∇yT (θ) by applying backpropagation. Furthermore, in Section 3, we explain how to efficiently compute µS(c) by using a ZDD-based technique [52]; importantly, its computation procedure also accepts the backpropagation.
While our work is built on the existing methods [59, 52, 42], none of them are intended to develop differentiable methods. A conceptual novelty of our work is its careful combination of those methods for developing a differentiable and accelerated optimization method, with which we can compute ∇yT (θ). This enables the application of gradient-based methods to the Stackelberg models of CCGs.

heading: 2.2 Details of Algorithm 1
text: We compute yT (θ) with Algorithm 1. Note that st, ct, and xt depend on θ, which is not explicitly indicated for simplicity. The most crucial part is Step 5, where we use softmin rather than argmin to
1Although the derivatives of y(θ) may not be unique, we abuse the notation and write ∇y(θ) for simplicity. As we will see shortly, we numerically approximate y(θ) with yT (θ) whose derivative, ∇yT (θ), exists uniquely. Therefore, when discussing our iterative differentiation method, we can ignore the abuse of notation.
Algorithm 1 Differentiable Frank–Wolfe-based equilibrium computation 1: s0 = c0 = 0, x−1 = x0 = µS(c0), and αt = t (t = 0, . . . , T ) 2: for t = 1, . . . , T : 3: st = st−1 − αt−1xt−2 + (αt−1 + αt)xt−1 4: ct = ct−1 + ηαt∇f ( 2 t(t+1)st;θ )
▷ ∇f(y;θ)i = ci(yi;θ) is differentiable in θ 5: Compute xt = µS(ct) with Algorithm 2 ▷ Differentiable softmin computation
return yT (θ) = 2T (T+1) ∑T t=1 αtxt
make the output differentiable in θ. Specifically, given any c ∈ Rn, we compute µS(c) as follows:
µS(c) := ∑ S∈S 1S exp
( −c⊤1S )∑ S′∈S exp (−c⊤1S′) .
Intuitively, each entry in c represents the cost of each i ∈ [n], and c⊤1S represents the cost of S ∈ S . We consider a probability distribution over S defined by softmin with respect to costs {c⊤1S}S∈S , and then marginalize it. The resulting vector is a convex combination of {1S}S∈S and thus always included in C. In the context of graphical modeling, this operation is called marginal inference [58]. Here, µS is defined by a summation over S, and explicitly computing it is prohibitively expensive. Section 3 details how to efficiently compute µS by leveraging the ZDD representations of S.

heading: 2.3 Convergence guarantee
text: Algorithm 1 is designed to induce Nesterov’s acceleration [43, 59] and achieves an O(1/T 2) convergence, which is faster than the O(1/T ) convergence of the original Frank–Wolfe algorithm. Theorem 1. Fix θ ∈ Θ and assume f(·;θ) to be L-smooth on Rd, i.e., Φ(z) := f(∑S∈S zS1S ;θ) (∀z ∈ Rd) satisfies Φ(z′) ≤ Φ(z) + ⟨∇Φ(z), z′ − z⟩+ L2 ∥z′ − z∥2 for all z, z′ ∈ Rd.2 If we let η ∈ [ 1CL , 14L ] for some C > 4, Algorithm 1 returns yT (θ) such that
f(yT (θ);θ)− f(y(θ);θ) ≤ O ( CL ln d
T 2
) .
We present the proof in Appendix B. In essence, softmin can be seen as a dual mirror descent step with the Kullback–Leibler divergence, and combining it with a primal gradient descent step yields the acceleration [3]. Although the acceleration technique itself is well studied, it has not been explored in the context of differentiable optimization. To the best of our knowledge, simultaneously making iterative optimization methods both differentiable and faster is a novel idea. This observation can be beneficial for developing other fast differentiable iterative algorithms. Experiments in Section 4.1 demonstrate that the acceleration indeed enhances the convergence speed in practice.
Note that the faster convergence enables us to more efficiently compute both yT (θ) and ∇yT (θ). The latter is because Algorithm 1 with a smaller T generates a smaller computation graph, which determines the computation complexity of the backpropagation for obtaining ∇yT (θ). Therefore, Algorithm 1 is suitable as an efficient differentiable pipeline between θ and yT (θ).

heading: 2.4 Implementation consideration: how to choose η and T
text: While Theorem 1 suggests setting η to 14L or less, this choice is often too conservative in practice. Thus, we should search for η values that bring high empirical performances. When using Algorithm 1 as a subroutine of gradient-based methods, it is repeatedly called to solve similar equilibrium computation instances. Therefore, a simple and effective way for locating good η values is to apply Algorithm 1 with various η values to example instances, as we will do in Section 4.1. We expect the empirical performance to improve with a line search of η, which we leave for future work.
2Smoothness parameter L defined on Rd can be, in general, exponentially large in n, albeit constant in T . How to alleviate the dependence on L remains an open problem.
s t 1 2 3 4 5
1 Root r 2
3 3
45
⊥ ⊤
To check whether the Frank–Wolfe algorithm has converged or not, we usually use the Frank–Wolfe gap [27], an upper-bound on an objective error. To obtain high-quality solutions, we terminate the algorithm when the gap becomes sufficiently small. In our case, however, our purpose is to obtain gradient information ∇yT (θ), and yT (θ) with a small T sometimes suffices to serve the purpose. By using a small T , we can reduce the computation complexity. Experiments in Section 4.2 show that the performance of a projected gradient method that uses ∇yT (θ) is not so sensitive to T values.

heading: 3 Efficient softmin computation with decision diagrams
text: This section describes how to efficiently compute µS(c) by leveraging ZDD representations of S. The main idea is to apply a technique called weight pushing [40] (or path kernel [56]) to ZDDs. A similar idea was used for obtaining efficient combinatorial bandit algorithms [52], but this research does not use it to develop differentiable algorithms. Our use of weight pushing comes from another important observation: it consists of simple arithmetic operations that accept reverse-mode automatic differentiation with respect to c, as shown in Algorithm 2. In other words, Algorithm 2 does not use, e.g., | · | or argmin. Therefore, ZDD-based weight pushing can be incorporated into the pipeline from θ to yT (θ) without breaking the differentiability.

heading: 3.1 Zero-suppressed binary decision diagrams
text: Given set family S ⊆ 2[n], we explain how to represent it with ZDD ZS = (V,A), a DAG-shaped data structure (see, e.g., Figure 1). The node set, V, has two terminal nodes ⊤ and ⊥ (they represent true and false, respectively) and non-terminal nodes. There is a single root, r ∈ V \ {⊥,⊤}. Each v ∈ V \ {⊥,⊤} has label lv ∈ [n] and two outgoing arcs, 1- and 0-arcs, which indicate whether lv is chosen or not, respectively. Let c0v , c 1 v ∈ V denote two nodes pointed by 0- and 1-arcs, respectively,
outgoing from v. For any v ∈ V \ {⊤,⊥}, let Rv ⊆ 2A be the set of all directed paths from v to ⊤. We define R := ⋃v∈V\{⊤,⊥} Rv. For any R ∈ R, let X(R) := { lv ∈ [n] ∣∣ (v, c1v) ∈ R}, i.e., labels of tails of 1-arcs belonging to R. ZDD ZS represents S as a set of r–⊤ paths: S = {X(R) | R ∈ Rr }. There is a one-to-one correspondence between S ∈ S and R ∈ Rr, i.e., S = X(R). Figure 1 presents an example of ZDD ZS , where S is the family of simple s–t paths. For example, S = {2, 5} ∈ S is represented in ZS by the red path, R ∈ Rr, with labels {1, 2, 3, 5, ⊤}. The labels of the tails of the 1-arcs form X(R) = {2, 5}, which equals S. We define the size of ZS = (V,A) by |ZS | := |V|. Note that |A| ≤ 2× |ZS | always holds. In general, the ZDD sizes and the complexity of constructing ZDDs can be exponential in n. Fortunately, many existing studies provide efficient methods for constructing compact ZDDs. One such method is the frontier-based search [30], which is based on Knuth’s Simpath algorithm [32]. Their method is particularly effective when S is a family of network substructures such as Hamiltonian paths, Steiner trees, matchings, and cliques. Furthermore, the family algebra [39, 32] of ZDDs enables us to deal with various logical constraints. Using those methods, we can flexibly construct ZDDs for various complicated combinatorial structures, e.g., Steiner trees whose size is at most a certain value.
Moreover, we can sometimes theoretically bound the ZDD sizes and the construction complexity. For example, if S consists of the aforementioned substructures on network G = (V,E) with a constant pathwidth, the ZDD sizes and the construction complexity are polynomial in |E| [30, 26].

heading: 3.2 Details of Algorithm 2 and computation complexity
text: Algorithm 2 computes x = µS(c) for any c = (c1, . . . , cn)⊤ ∈ Rn. First, it computes {Bv}v∈V in a bottom-up topological order of ZS . Note that Bv = ∑ S∈{X(R) | R∈Rv } exp(− ∑ i∈S ci) holds. Then it computes {Pv}v∈V. Each Pv indicates the probability that a top-down random walk starting from root node r reaches v ∈ V, where we choose 0-arc (1-arc) with probability p0 (p1). From B⊥ = 0 and the construction of {Bv}v∈V, the random walk never reaches ⊥, and its trajectory R ∈ Rr recovers X(R) ∈ S with a probability proportional to exp(−∑i∈X(R) ci) = exp(−c⊤1X(R)). Therefore, by summing the probabilities of reaching v and choosing a 1-arc outgoing from v for each i ∈ [n] as in Step 9, we obtain x = µS(c). In practice, we recommend implementing Algorithm 2 with the log-sum-exp technique and double-precision computations for numerical stability.
Algorithm 2 runs in O(|ZS |) time, and thus Algorithm 1 takes O((n+C∇+ |ZS |)T ) time, where C∇ is the cost of computing ∇f . From the cheap gradient principle [21], the complexity of computing ∇F (θ,yT (θ)) with backpropagation is almost the same as that of computing F (θ,yT (θ)). That is, smaller ZDDs make the computation of both yT (θ) and ∇yT (θ) faster. Therefore, our method significantly benefits from the empirical compactness of ZDDs. Note that we can construct ZDD ZS in a preprocessing step; once we obtain ZS , we can reuse it every time µS(·) is called.

heading: 4 Experiments
text: Section 4.1 confirms the benefit of acceleration to empirical convergence speed. Section 4.2 demonstrates the usefulness of our method via experiments on communication network design instances. Section 4.3 presents experiments with small instances to see whether our method can empirically find globally optimal θ. Due to space limitations, we present full experimental results in Appendix C.
All the experiments were performed using a single thread on a 64-bit macOS machine with 2.5 GHz Intel Core i7 CPUs and 16 GB RAM. We used C++11 language, and the programs were compiled by Apple clang 12.0.0 with -O3 -DNDEBUG option. We used Adept 2.0.5 [23] as an automatic differentiation package and Graphillion 1.4 [25] for constructing ZDDs, where we used a beamsearch-based path-width optimization method [26] to specify the traversal order of edges. The source code is available at https://github.com/nttcslab/diff-eq-comput-zdd.
Problem setting. We address Stackelberg models for optimizing network parameters θ ∈ Rn, where [n] represents an edge set. We focus on two situations where combinatorial strategies S ⊆ 2[n] are Hamiltonian cycles and Steiner trees. The former is a variant of the selfish-routing setting, and the latter arises when designing communication networks as in Section 1.1. Note that in both settings, common operations on C, e.g., projection and linear optimization, are NP-hard. We use two types of cost functions: fractional cost ci(yi;θ) = di(1 + C × yi/(θi + 1)) and exponential cost ci(yi;θ) = di(1 + C × yi exp(−θi)), where di ∈ (0, 1] is the length of the i-th edge (normalized so that maxi∈[n] di = 1 holds) and C > 0 controls how heavily the growth in yi (congestion) affects cost ci. We set C = 10. Note that edge i with a larger θi is more tolerant to congestion. The leader aims to minimize social cost F (θ,y(θ)) = ∑ i∈[n] ci(yi(θ);θ)yi(θ). In realistic situations, the leader cannot let all edges have sufficient capacity due to budget constraints. To model this situation, we impose a constraint on θ by defining Θ = { θ ∈ Rn≥0
∣∣ θ⊤1 = n}, where 1 is the all-one vector. Datasets. Table 1 summarizes the information about datasets and ZDDs used in the experiments. For the Hamiltonian-cycle setting (HAMILTON), we used att48 (ATT) and dantzig42 (DANTZIG) datasets in TSPLIB [48]. Following [11, 44], we obtained networks in Figure 4 using Delaunay triangulation [54]. For the Steiner-tree setting (STEINER), we used Uninett 2011 (UNINETT) and TW Telecom (TW) networks of Internet Topology Zoo [31]. We selected terminal vertices as shown in Figure 4. We can see in Table 1 that ZDDs are much smaller than the strategy sets. As mentioned in Section 3.2, we can construct ZDDs in a preprocessing step, and the construction times were so short as to be negligible compared with the times taken for minimizing the social cost (see Section 4.2). Therefore, we do not take the construction times into account in what follows.

heading: 4.1 Empirical convergence of equilibrium computation
text: We studied the empirical convergence of Algorithm 1 with acceleration (w/ A), where we let η = 0.05, 0.1, 0.2, and 0.5. We applied it to the minimization problems of form miny∈C f(y;θ), where f is a potential function defined by cost function ci(yi;θ) (fractional or exponential). We let θ = 1.
For comparison, we used two kinds of baselines. One is a differentiable Frank–Wolfe algorithm without acceleration (w/o A), which just replaces argmin with softmin as explained in Section 2.1. To guarantee the convergence of the modified algorithm, we let ηt = η0 × t (η0 = 0.1, 1.0, and 10.0). The other is the standard non-differentiable Frank–Wolfe algorithm (FW) implemented as in [27].
Figure 2 shows how quickly the Frank–Wolfe gap [27], which is an upper bound of the objective error, decreased as the number of iterations and the computation time increased. w/ A and w/o A tend to be faster and slower than FW, respectively. That is, Algorithm 1 (w/ A) becomes both differentiable and faster than the original FW, while the naive modified one (w/o A) becomes differentiable but slower. As in the TW results, however, w/ A with a too large η (η = 0.5) sometimes failed to be accelerated; this is reasonable since Theorem 1 requires η to be a moderate value. Thus, if we can locate appropriate η, Algorithm 1 achieves faster convergence in practice. As discussed in Section 2.4, we can search for η by examining the empirical convergence for various η values, as we did above.

heading: 4.2 Stackelberg models for designing communication networks
text: We consider minimizing social cost F (θ,y(θ)). To this end, roughly speaking, we should assign large θi values to edges with large yi values. We applied the projected gradient method with a step size of 5.0 to problem (1). We approximated ∇F (θ,y(θ)) by applying automatic differentiation to F (θ,yT (θ)), where yT (θ) was computed by Algorithm 1 with T = 100, 200, and 300.
To the best of our knowledge, no existing methods can efficiently deal with the problems considered here due to the complicated structures of S . Therefore, as a baseline method, we used the following iterative heuristic. Given current y(θ), we replace θ with θ + δ(y(θ)− ȳ(θ)1), where δ > 0 and ȳ(θ) = 1n ∑ i∈[n] yi(θ). That is, we increase/decrease θi if edge i is used more/less than average. We then project θ onto Θ and compute y(θ) with Algorithm 1 (T = 300). If F (θ,y(θ)) value does not decrease after the above update, we restart from a random point in Θ.
Figure 3 compares our method (η = 0.1) and the baseline on ATT and TW instances, where both started from θ = 1 and continued to update θ for two minutes. Our method found better θ values than the baseline. The results only imply the empirical tendency, and our method is not guaranteed to find globally optimal θ. Nevertheless, experiments in Section 4.3 show that it tends to find a global optimum at least for small instances. Figure 4 shows the yT (θ) and θ values obtained by our method with η = 0.1 and T = 300 for fractional-cost instance. We confirmed that large θi values were successfully assigned to edges with large yi values. Those results demonstrate that our method is useful for addressing Stackelberg models of CCGs with complicated combinatorial strategy sets S.

heading: 4.3 Experiments on empirical convergence to global optimum
text: We performed additional experiments on small instances to see whether the projected gradient method used in Section 4.2 can empirically find θ that is close to being optimal. We used small selfish-routing instances, where a graph is given by the left one in Figure 1 and the edges are numbered from 1 to 5 in that order. We let strategy set S be the set of all simple s-t paths. The cost functions and feasible region Θ were set as with those in the above sections. Our goal is to minimize social cost F (θ,y(θ)).
As in Section 4.2, we computed y(θ) using Algorithm 1 with η = 0.1 and T = 300, and performed the projected gradient descent to minimize F (θ,y(θ)), where gradient ∇F (θ,y(θ)) was computed with backpropagation. On the other hand, to obtain (approximations of) globally optimal θ, we performed an exhaustive search over the feasible region, where the step size was set to 0.05. Regarding computation times, the projected gradient method converged in less than 30 iterations, which took less than 20 ms, while the exhaustive search took about 1000 seconds.
Results on fractional costs. A global optimum found by the exhaustive search was θ = (0, 2.5, 0, 0, 2.5), whose social cost F (θ,y(θ)) was 6.444. Our method started from θ = 1, whose social cost was 7.000, and returned θ = (1.25, 1.25, 0, 1.25, 1.25) with social cost 6.444. Although the solution is different from that of the exhaustive search, both attain the identical social cost. Thus, the solution returned by the projected gradient method is also globally optimal.
Results on exponential costs. A global optimum found by the exhaustive search was θ = (0, 2.5, 0, 0, 2.5) with social cost 3.517. Our method started from θ = 1 with social cost 5.678 and reached θ = (0, 2.5, 0, 0, 2.5) with social cost 3.517. Along the way, the projected gradient method was about to be trapped in θ = (1.25, 1.25, 0, 1.25, 1.25) with social cost 4.865, which seems to be a saddle point. However, it successfully got out of there and reached the global optimum.

heading: 5 Conclusion and discussion
text: We proposed a differentiable pipeline that connects CCG parameters to their equilibria, enabling us to apply gradient-based methods to the Stackelberg models of CCGs. Our Algorithm 1 leverages softmin to make the Frank–Wolfe algorithm both differentiable and faster. ZDD-based softmin computation (Algorithm 2) enables us to efficiently deal with complicated CCGs. It also naturally works with automatic differentiation, offering an easy way to compute desired derivatives. Experiments confirmed the accelerated empirical convergence and practicality of our method.
An important future direction is further studying theoretical aspects. From our experimental results, ∇yT (θ) is expected to converge to ∇y(θ), although its theoretical analysis is very difficult. Recently, some relevant results have been obtained for simple cases where iterative optimization methods are written by a contraction map defined on an unconstrained domain [1, 20]. In our CCG cases, however, we need to study iterative algorithms that numerically solve the constrained potential minimization, which requires a more profound understanding of iterative differentiation approaches. Another interesting future work is to make linearly convergent Frank–Wolfe variants [33] differentiable.
Finally, we discuss limitations and possible negative impacts. Our work does not cover cases where minimizer y(θ) of potential functions is not unique. Since the complexity of our method mainly depends on the ZDD sizes, it does not work if ZDDs are prohibitively large, which can happen when strategy sets consist of the substructures of dense networks. Nevertheless, many real-world networks are sparse, and thus our ZDD-based method is often effective, as demonstrated in experiments. At a meta-level, optimizing social infrastructures in terms of a single objective function (e.g., the social cost) may lead to an extreme choice that is detrimental to some individuals. We hope our method can provide a basis for designing social infrastructures that are beneficial for all.

heading: Acknowledgements
text: The authors thank the anonymous reviewers for their valuable feedback, corrections, and suggestions. This work was partially supported by JST ERATO Grant Number JPMJER1903 and JSPS KAKENHI Grant Number JP20H05963.

Review summaries:
This paper studies the problem of designing parameters of combinatorial congestion games, where selfish non-atomic agents choose an optimal route/network to optimize their own objectives. The combinatorial congestion games are a special type of potential games. The equilibrium of a congestion game with given parameters can be solved by optimizing the potential function. The problem of optimizing the game parameters is therefore a bilevel optimization problem with an equilibrium (potential function optimization problem) involved in the inner problem. Benefited from the recent development of differentiable optimization, the authors propose to solve the bilevel optimization problem by backpropagating the gradient from the argmax of the inner problem to the outer objective to run gradient descent. This end-to-end gradient descent approach can avoid solving the bilevel optimization problem with complex constraints by directly running end-to-end gradient descent.  The challenge involved in this end-to-end differentiable optimization approach is how to efficiently compute the gradient of a complex optimization problem (equilibrium computation/potential function optimization) with high-dimensional constraints. The authors leverage the Frank-Wolfe algorithm to solve the optimization problem and to get the corresponding gradient of the optimal solution. The use of the Frank-Wolfe algorithm decomposes the optimization problem into an iterative algorithm composed of multiple minimization problems, where each of them is relaxed to a softmin in order to make them differentiable. To solve the softmin problem with complex constraints, they use zero-suppressed binary decision diagrams (ZDD) to compactly represent the decision variables, which can speed up the optimization step while keeping the differentiability. The authors also utilize Nesterov's acceleration to improve the convergence of the Frank-Wolfe algorithm to $O(1/T^2)$, with the differentiability maintained.

The authors of this paper aim to find a combinatorial congestion game that optimizes a given objective over some parameterized class of games, assuming that the players in the game attain a Wardrop equilibrium. The problem is formulated using a framework of Stackelberg bi-level optimization models, where the designer (leader) chooses the game's parameters, and the followers then play the induced game. The authors propose to apply gradient-based methods to approximate the solution. To this end, it remains critical to be able to compute the followers' equilibrium and estimate its gradient, i.e., the changes in the equilibrium when the parameters of the game shift. For this purpose, the authors introduce a differentiable variant of the Frank–Wolfe algorithm. Instead of the strict minimization included in the original version, the authors suggest using softmin. Computing the softmin may be computationally demanding because it sums over all possible strategies, which is (in general) a set exponential in the problem size. More efficient computation can be achieved using the so-called zero-suppressed binary decision diagrams (ZDDs), a data structure enabling calculation of the softmin in time linear in its size. Because the size of the ZDDs is often many orders of magnitude smaller than the cardinality of the set of strategies, ZDDs offer a significant speedup in practice. In the last part of the manuscript, the authors present empirical results achieved with their algorithm and its variants on two classes of games per each strategy representation (Hamiltonian cycles and Steiner trees). The results show that the algorithm can converge to Wardrop equilibria faster than baselines, and the gradient descent running on top of it consequently designs better congestion games than the heuristic baselines. 

This paper focus on the Stackelberg model of combinatorial congestion games (CCGs), which is important for the case involved the leader and multiple self-interested followers. The main contribution of this paper is proposing a fully differentiable framework which can take the advantage of automatic differentiation.  My main concern is the contribution of this paper. There are two key techiniques in the proposed algorithms, softmin and Zero-suppressed binary decision diagrams (ZDD), where softmin is a widely used technique to make argmax operation to be differetiable and the ZDD technique is proposed in the 2020 AAAI paper, Practical Frank–Wolfe Method with Decision Diagramsfor Computing Wardrop Equilibrium of Combinatorial Congestion Games. I would like to ask authors to clearly state the difference of the techniques in this paper and the one in the 2020 AAAI paper.  Minor issues and typos:  "practical utility" in line 60, "utility" -> usage or application? "Our method introduces benefits by improving the designs of social infrastructures" is strange.

The problem of computing equilibria in CGs is consider an easy task whenever the strategy space has a compact representation, since it turns out that the users minimize a convex function. When having to decide on some parameters of the CG in order to optimize the performance at equilibrium we get a bi-level optimization problem, where in the upper level the designer sets the parameters and in the lower level the users optimize their function reaching an equilibrium. Solving this bi-level problem is regarded to be far more difficult. This work deals with the latter problem and combines optimization techniques to give an algorithm that performs well in practice.