# Solving Structured Hierarchical Games Using Differential Backward Induction

Zun Li  University of Michigan, Ann Arbor
<EMAIL_ADDRESS>Feiran Jia  Pennsylvania State University
<EMAIL_ADDRESS>Aditya Mate  Harvard University
<EMAIL_ADDRESS>Shahin Jabbari  Drexel University
<EMAIL_ADDRESS>Mithun Chakraborty  University of Michigan, Ann Arbor
<EMAIL_ADDRESS>Milind Tambe  Harvard University
<EMAIL_ADDRESS>Yevgeniy Vorobeychik  Washington
University in St. Louis
<EMAIL_ADDRESS>
###### Abstract

From large-scale organizations to decentralized political systems,
hierarchical strategic decision making is commonplace. We introduce a novel
class of _structured hierarchical games (SHGs)_ that formally capture such
hierarchical strategic interactions. In an SHG, each player is a node in a
tree, and strategic choices of players are sequenced from root to leaves, with
root moving first, followed by its children, then followed by their children,
and so on until the leaves. A player’s utility in an SHG depends on its own
decision, and on the choices of its parent and _all_ the tree leaves. SHGs
thus generalize simultaneous-move games, as well as Stackelberg games with
many followers. We leverage the structure of both the sequence of player moves
as well as payoff dependence to develop a gradient-based back propagation-
style algorithm, which we call _Differential Backward Induction (DBI)_ , for
approximating equilibria of SHGs. We provide a sufficient condition for
convergence of DBI and demonstrate its efficacy in finding approximate
equilibrium solutions to several SHG models of hierarchical policy-making
problems.

##  1 Introduction

The COVID-19 pandemic has revealed considerable strategic tension among the
many parties involved in decentralized hierarchical policy-making. For
example, recommendations by the World Health Organization are sometimes
heeded, and other times discarded by nations, while subnational units, such as
provinces and urban areas, may in turn take a policy stance (such as on
lockdowns, mask mandates, or vaccination priorities) that is not congruent
with national policies. Similarly, in the US, policy recommendations at the
federal level can be implemented in a variety of ways by the states, while
counties and cities, in turn, may comply with state-level policies, or not,
potentially triggering litigation [15]. Central to all these cases is that,
besides this strategic drama, what ultimately determines infection spread is
how policies are implemented _at the lowest level_ , such as by cities and
towns, or even individuals. Similar strategic encounters routinely play out in
large-scale organizations, where actions throughout the management hierarchy
are ultimately reflected in the decisions made at the lowest level (e.g., by
the employees who are ultimately involved in production), and these lowest-
level decisions play a decisive role in the organizational welfare.

We propose a novel model of hierarchical decision making which is a natural
stylized representation of strategic interactions of this kind. Our model,
which we term _structured hierarchical games (SHGs)_ , represents each player
by a node in a tree hierarchy. The tree plays two roles in SHGs. First, it
captures the sequence of moves by the players: the root (the lone member of
level 1 of the hierarchy) makes the first strategic choice, its children
(i.e., all nodes in level 2) observe the root’s choice and follow, their
children then follow in turn, and so on, until we reach the leaf node players
who move upon observing their predecessors’ choices. Second, the tree
partially captures strategic dependence: a player’s utility depends on its own
strategy, that of its parent, and the strategies of _all of the leaf nodes_.
The sequence of moves in our model naturally captures the typical sequence of
decisions in hierarchical policy-making settings, as well as in large
organizations, while the utility structure captures the decisive role of leaf
nodes (e.g., individual compliance with vaccination policies), as well as
hierarchical dependence (e.g., employee dependence on a manager’s approval of
their performance, or state dependence on federal funding). Significantly, the
_SHG_ model generalizes a number of well-established models of strategic
encounters, including (a) simultaneous-move games (captured by a 2-level SHG
with the root having a single dummy action), (b) Stackelberg (leader-follower)
games (a 2-level game with a single leaf node) [37, 11], and (c) single-leader
multi-follower Stackelberg games (e.g., a Stackelberg security game with a
single defender and many attackers) [5, 8].

Our second contribution is a gradient-based algorithm for approximately
computing subgame-perfect equilibria of _SHGs_. Specifically, we propose
_Differential Backward Induction (DBI)_ , which is a backpropagation-style
gradient ascent algorithm that leverages both the sequential structure of the
game, as well as the utility structure of the players. As _DBI_ involves
simultaneous gradient updates of players in the same level (particularly at
the leaves), convergence is not guaranteed in general (as is also the case for
best-response dynamics [12]). Viewing _DBI_ as a dynamical system, we provide
a sufficient condition for its convergence to a stable point. Our results also
imply that in the special case of two-player zero-sum Stackelberg games, _DBI_
converges to a local Stackelberg equilibrium [11, 39].

Finally, we demonstrate the efficacy of DBI in finding approximate equilibrium
solutions to several classes of SHGs. First, we use a highly stylized class of
SHGs with polynomial utility functions to compare DBI with five baseline
gradient-based approaches from prior literature. Second, we use DBI to solve a
recently proposed game-theoretic model of 3-level hierarchical epidemic policy
making. Third, we apply DBI to solve a hierarchical variant of a public goods
game, which naturally captures the decentralization of decision making in
public good investment decisions, such as investments in sustainable energy.
Fourth, we evaluate DBI in the context of a hierarchical security investment
game, where hierarchical decentralization (e.g., involving federal government,
industry sectors, and particular organizations) can also play a crucial role.
In all of these, we show that DBI significantly outperforms the state of the
art approaches that can be applied to solve games with hierarchical structure.

Related Work  SHGs generalize both simultaneous-move games and Stackelberg
games with multiple followers [23, 5]. They are also related to graphical
games [20] in capturing utility dependence structure, although SHGs also
capture sequential structure of decisions. Several prior approaches use
gradient-based methods for solving games with particular structure. A
prominent example is generative adversarial networks (GANs), though these are
zero-sum games [14, 19, 31, 9, 29, 30]. Ideas from learning GANs have been
adopted in gradient-based approaches to solve multi-player general-sum games
[28, 4, 7, 17, 22, 25, 29]. However, all of these approaches assume a
simultaneous-move game. A closely-related thread to our work considers
gradient-based methods for bi-level optimization [24, 35]. Several related
efforts consider gradient-based learning in Stackelberg games, and also use
the implicit function theorem to derive gradient updates [2, 11, 32, 39, 38].
We significantly generalize these ideas by considering an arbitrary
hierarchical game structure.

Jia et al. [18] recently considered a stylized 3-level SHG for pandemic policy
making, and proposed several non-gradient-based algorithms for this problem.
We compare with their approach in Section 4.

##  2 Structured Hierarchical Games

#### Notation

We use bold lower-case letters to denote vectors. Let $f$ be a function of the
form
$f(\bm{x},\bm{y}):\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}\rightarrow\mathbb{R}^{d^{\prime\prime}}$.
We use $\nabla_{\bm{x}}f$ to denote the partial derivative of $f$ with respect
to $\bm{x}$. When there is functional dependency between $\bm{x}$ and
$\bm{y}$, we use $D_{\bm{x}}f$ to denote the total derivative of
$f(\bm{x},\bm{y}(\bm{x}))$ with respect to $\bm{x}$. We use
$\nabla^{2}_{\bm{x},\bm{x}}f$ and $\nabla^{2}_{\bm{x},\bm{y}}f$ to denote the
second-order partial derivatives and $D^{2}_{\bm{x},\bm{x}}f$ to denote the
second-order total derivative of $f$. For a mapping
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, we use $f^{t}(\bm{x})$ to denote
$t$ iterative applications of $f$ on $\bm{x}$. For mappings
$f_{1}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and
$f_{2}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, we define $(f_{1}\circ
f_{2})(\bm{x})\triangleq f_{1}(f_{2}(\bm{x}))$ and
$(f_{1}+f_{2})(\bm{x})\triangleq f_{1}(\bm{x})+f_{2}(\bm{x})$. Moreover, for a
given $\epsilon\in\mathbb{R}^{>0}$ and $\bm{x}\in\mathbb{R}^{d}$, we define
the $\epsilon$-ball around $\bm{x}$ as
$\mathbb{B}_{\epsilon}(\bm{x})=\\{\bm{x}^{\prime}\in\mathbb{R}^{d}\mid\|\bm{x}-\bm{x}^{\prime}\|_{2}<\epsilon\\}$.
Finally, $\bm{I}$ denotes an identity matrix.

#### Formal Model

A structured hierarchical game (SHG) $\mathcal{G}$ consists of the set
$\mathscr{N}$ of $n$ players. Each player $i$ is associated with a set of
actions $\mathcal{X}_{i}\subseteq\mathbb{R}^{d_{i}}$. The players are
partitioned across $L$ levels, where $\mathscr{N}_{l}$ is the set of $n_{l}$
players occupying level $l$. Let $l_{i}$ denote the level occupied by player
$i$. This _hierarchical_ structure of the game is illustrated in Figure 1
where players correspond to nodes and levels are marked by dashed boundaries.
The hierarchy plays two crucial roles: 1) it determines the order of moves,
and 2) it partly determines utility dependence among players. Specifically,
the temporal pattern of actions is as follows: level 1 has a single player,
the root, who chooses an action first, followed by all players in level 2
making simultaneous choices, followed in turn by players in level 3, and so on
until the _leaves_ in the final level $L$. Players of level $l$ only observe
the actions chosen by all players of levels $1,2,...,l-1$, but not their peers
in the same level.

Figure 1: Schematic representation of an SHG. The utility of player $i$ can
have direct functional dependence _only_ on the joint action of _all_ shaded
players.

So, for example, pandemic social distancing and vaccination policies in the US
are initiated by the federal government (including the Centers for Disease
Control and Prevention who acts as the root in our game model), with states
(second level nodes) subsequently instituting their own policies, counties
(third level nodes) reacting to these by determining their own, and behavior
of people (leaf nodes) ultimately influenced, but not determined, by the
guidelines and enforcement policies by the local county/city.

Next, we describe the utility structure of the game as entailed by the SHG
hierarchy. Each player $i$ in level $l_{i}>1$ (i.e., any node other than the
root) has a _unique parent_ in level $l_{i}-1$; we denote the parent of node
$i$ by $\textsc{Pa}(i)$. A player’s utility function is determined by 1) its
own action, 2) the action of its parent, and 3) the actions of _all_ players
in level $L$ (i.e., all _leaf_ players). To formalize, let $\bm{x}_{l}$ denote
the joint action profile of all players in level $l$. Player $i$’s utility
function then has the form $u_{i}(x_{i},\bm{x}_{L})$ if $l_{i}=1$,
$u_{i}(x_{i},x_{\textsc{Pa}(i)},\bm{x}_{L})$ if $1<l_{i}<L$, and
$u_{i}(x_{i},x_{\textsc{Pa}(i)},\bm{x}_{L,-i})$ if $l_{i}=L$, where
$\bm{x}_{L,-i}$ is the action profile of all players in level $L$ other than
$i$. For example, in our running pandemic policy example, the utility of a
county depends on both the policy and enforcement strategy of its state (its
_parent_) and on the ultimate pandemic spread and economic impact within it,
both determined largely by the behavior of the county residents (_leaf
nodes_). Note the considerable generality of the SHG model. For example, an
arbitrary simultaneous-move game is a SHG with 2 levels and a “dummy” root
node (utilities of all leaves depend on one another’s actions), and an
arbitrary Stackelberg game (e.g., Stackelberg security game), even with many
followers, can be modeled as a 2-level SHG with the leader as root and
followers as leaves. Furthermore, while we have defined SHGs with respect to
real-vector player action sets, it is straightforward to represent mixed
strategies of finite-action games in this way by simply using a softmax
function to map an arbitrary real vector into a valid mixed strategy.

Solution Concept  Since an SHG has important sequential structure, it is
natural to consider the _subgame perfect equilibrium (SPE)_ as the solution
concept [33]. Here, we focus on pure-strategy equilibria. To begin, we note
that in SHGs, the strategies of players in any level $l>1$ are, in general,
functions of the complete history of play in levels $1,\dots,l-1$, which we
denote by $h_{<l}=(\bm{x}_{1},\bm{x}_{2},\dotsc,\bm{x}_{l-1})$. Formally, a
(pure) strategy of a player $i$ is denoted by $s_{i}(h_{<l})$, which
deterministically maps an arbitrary history $h_{<l}$ into an action
$x_{i}\in\mathcal{X}_{i}$. A _Nash equilibrium_ of an SHG is then a strategy
profile $\bm{s}=(s_{1},\dotsc,s_{i},\ldots,s_{n})$ such that for all
$i\in\mathscr{N}$, $u_{i}(s_{i},\bm{s_{-i}})\geq
u_{i}(s_{i}^{\prime},\bm{s_{-i}})$ for all possible alternative strategies for
$i$, $s_{i}^{\prime}$. Here, we denote the realized payoff of $i$ from profile
$\bm{s}$ by $u_{i}(s_{i},\bm{s_{-i}})$. Next, we define a _level- $l$-subgame_
given $h_{<l}$ as an SHG that includes only players at levels $\geq l$, with
actions chosen in levels $<l$ fixed to $h_{<l}$. A strategy profile $\bm{s}$
is a _subgame perfect equilibrium_ of SHG if it is a Nash equilibrium of every
level-$l$-subgame of SHG for every $l$ and history $h_{<l}$. We prove in
appendix A that our definition of SPE is equivalent to the standard SPE in an
extensive-form representation of SHG.

While in principle we can compute an SPE of an SHG using backward induction,
this cannot be done directly (i.e., by complete enumeration of actions of all
players) as actions are real vectors. Moreover, even discretizing actions is
of little help, as the hierarchical nature of the game leads to exponential
explosion of the search space. We now present a gradient-based approach for
approximating SPE along the equilibrium path in an SHG that leverages the game
structure to derive backpropagation-style gradient updates.

##  3 Differential Backward Induction

In this section, we describe our gradient-based algorithm, Differential
Backward Induction (DBI), for approximating an SPE (which we mean hereinafter
finding a joint-action profile $\bm{x}$ that constitutes a subgame-perfect
equilibrium path), and then analyze its convergence. Just as gradient ascent
does not, in general, identify a globally optimal solution to a non-convex
optimization problem, DBI in general yields a solution which only satisfies
first-order conditions (see Section 3.2 for further details). Moreover, we
leverage the structure of the utility functions to focus computation on an SPE
in which strategies of players are only a function of their immediate
parents.111Note that while we cannot guarantee that an SPE exists in SHGs in
general, let alone those possessing the assumed structure, we find
experimentally that our approach often yields good SPE approximations.

In this spirit, we define _local_ best response functions
$\phi_{i}:\mathbb{R}^{d_{\textsc{Pa}(i)}}\rightarrow\mathbb{R}^{d_{i}}$
mapping a player $i$’s parent’s action $x_{\textsc{Pa}(i)}$ to $i$’s action
$x_{i}$; note that the notation $\phi_{i}$ is distinct from $s_{i}$ above for
$i$’s strategy to emphasize the fact that $\phi_{i}$ is only locally optimal.
Now, suppose that a player $i$ is in the last level $L$. Local optimality of
$\phi_{i}$ implies that if $x_{i}=\phi_{i}(x_{\textsc{Pa}(i)})$, then
$\nabla_{x_{i}}u_{i}\left(\bm{x}_{i},\bm{x}_{\textsc{Pa}(i)},\bm{x}_{L,-i}\right)=0$
and
$\nabla^{2}_{x_{i},x_{i}}u_{i}\left(\bm{x}_{i},\bm{x}_{\textsc{Pa}(i)},\bm{x}_{L,-i}\right)\prec
0.$222For simplicity, we omit degenerate cases where
$\nabla^{2}_{\bm{x}_{i},\bm{x}_{i}}u_{i}=0$ and assume all local maxima are
strict.

Let $\phi_{l}$ denote the local best response for all the players in level $l$
given the actions of all players in level $l-1$. We can compose these local
best response functions to define the function
$\Phi_{l}:=\phi_{L}\circ\phi_{L-1}\circ\dotsc\circ\phi_{l+1}:\mathbb{R}^{d_{n_{l}}}\rightarrow\mathbb{R}^{d_{n_{L}}}$
i.e., the local best response of players in the last level $L$ given the
actions of the players in level $l$.333Note that in particular
$\Phi_{L}=\phi_{L}$. Then for any player $(i)$ in level $l_{i}<L$,
$D_{x_{i}}u_{i}\left(x_{i},x_{\textsc{Pa}(i)},\Phi_{l}\left(\langle
x_{i},\bm{x}_{l,-i}\rangle\right)\right)=0$ and
$D^{2}_{x_{i},x_{i}}u_{i}\left(x_{i},x_{\textsc{Pa}(i)},\Phi_{l}\left(\langle
x_{i},\bm{x}_{l,-i}\rangle\right)\right)\prec 0$, where $D_{x_{i}}$ is the
total derivative with respect to $x_{i}$ (as $\Phi_{l}(\langle
x_{i},\bm{x}_{l,-i}\rangle)$ is also a function of $x_{i}$). Note that the
functions $\phi$ and $\Phi$ are _implicit_ , capturing the functional
dependencies between actions of players in different levels at the local
equilibrium.

Throughout, we make the following standard assumption on the utility functions
[10, 39].

######  Assumption 1.

For any $x_{i}\in\mathcal{X}_{i}$, the second-order partial derivatives of the
form $\nabla^{2}_{x_{i},x_{i}}u_{i}$ are non-singular.

###  3.1 Algorithm

The DBI algorithm works in a bottom-up manner, akin to back-propagation: for
each level $l$, we compute the total derivatives (gradients) of the utility
functions and local best response maps ($\phi$, $\Phi$) based on analytical
expressions that we derive below. We then propagate this information up to
level $l-1$, as it is used to compute gradients for that level, and so on
until level 1. Algorithm 1 gives the full DBI algorithm. In this algorithm,
$\textsc{Chd}(i)$ denotes the set of children of player $i$ (i.e., nodes in
level $l_{i}+1$ for whom $i$ is the parent).

An SHG instance $\mathcal{G}$

Learning rate $\alpha$, maximum number of iterations $T$ for gradient update

A strategy profile

Randomly initialize
$\bm{x}^{0}=\langle\bm{x}_{1}^{0},\ldots,\bm{x}_{L}^{0}\rangle$

for $t=1,2,\dotsc,T$ do

for $l=L,L-1,\dotsc,1$ do

for $i=1,2,\dotsc,n_{l}$ do

if $l=L$ then

Back-propagate $D_{x_{i}}\Phi_{i}=\bm{I}$ to $\textsc{Pa}(i)$

Set $x^{t}_{i}\leftarrow x^{t-1}_{i}+\alpha\nabla_{x_{i}}u_{i}$

else

Compute $\nabla_{x_{i}}u_{i},\nabla_{\bm{x}_{L}}u_{i}$ at $\bm{x}^{t-1}$

Compute $D_{x_{i}}\phi_{j},\forall j\in\textsc{Chd}(i)$ (Eqn. (5))

Compute $D_{x_{i}}\Phi_{l}$ (Eqn. (4))

Back-propagate $D_{x_{i}}\Phi_{l}$ to $\textsc{Pa}(i)$

Compute
$D_{x_{i}}u_{i}=\nabla_{x_{i}}u_{i}+\nabla_{\bm{x}_{L}}u_{i}D_{x_{i}}\Phi_{l}$

Set $x^{t}_{i}\leftarrow x^{t-1}_{i}+\alpha D_{x_{i}}u_{i}$

Return $\bm{x}^{T}$

Algorithm 1 Differential Backward Induction (DBI)

DBI works in a backward message-passing manner, comparable to back-
propagation: after each player has computed its total derivative, it passes
(back-propagates) $D_{x_{i}}\Phi_{l}$ to its direct parent; this information
is, in turn, used by the parent to compute its own total derivative, which is
passed to its own parent, and so on.

Algorithm 1 takes the total derivates as given. We now derive closed-form
expressions for these. We start from the last level $L$. Given the actions of
players in level $L-1$, the total derivative of a player $i\in\mathscr{N}_{L}$
with respect to $x_{i}$ is

$D_{x_{i}}u_{i}\left(x_{i},x_{\textsc{Pa}(i)},\bm{x}_{L,-i}\right)=\nabla_{x_{i}}u_{i}.$
(1)

For a player $i$ in level $L-1$, the total derivative (at a local best
response) is

$\displaystyle D_{x_{i}}u_{i}(x_{i},x_{\textsc{Pa}(i)},$
$\displaystyle\phi_{L}(\langle x_{i},\bm{x}_{L-1,-i}\rangle))$
$\displaystyle=\nabla_{x_{i}}u_{i}+\left(\nabla_{\bm{x}_{L}}u_{i}\right)\left(D_{x_{i}}\phi_{L}\right),$
(2)

where $\nabla_{\bm{x}_{L}}u_{i}$ is a $1\times d_{n_{L}}$ vector and
$D_{x_{i}}\phi_{L}$ is a $d_{n_{L}}\times d$ matrix. The technical challenge
here is to derive the term $D_{x_{i}}\phi_{L}$ for $i\in\mathscr{N}_{L-1}$.
Recall that $\phi_{L}$ is the vectorized concatenation of the $\phi_{j}$
functions for $j\in\mathscr{N}_{L}$. Since the local best response strategy of
a player in level $L$ only depends on its parent in level $L-1$, the only
terms in $\phi_{L}$ that depend on $x_{i}$ are the actions of
$\textsc{Chd}(i)$ in level $L$. Consequently, it suffices to derive
$D_{x_{i}}\phi_{j}$ for $j\in\textsc{Chd}(i)$. Note that for these players
$j$, $\nabla_{x_{j}}u_{j}=0$ (by local optimality of $\phi_{L}$). We will use
this first-order condition to derive the expression for the total derivative
using the _implicit function theorem_.

######  Theorem 1 (Implicit Function Theorem (IFT) [10, Theorem 1B.1]).

. Let
$f(\bm{x}_{1},\bm{x}_{2}):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$
be a continuously differentiable function in a neighborhood of
$(\bm{x}_{1}^{*},\bm{x}_{2}^{*})$ such that
$f(\bm{x}_{1}^{*},\bm{x}_{2}^{*})=0$. Also suppose $\nabla_{\bm{x}_{2}}f$, the
Jacobian of $f$ with respect to $\bm{x}_{2}$, is non-singular at
$(\bm{x}_{1}^{*},\bm{x}_{2}^{*})$. Then around a neighborhood of
$\bm{x}_{1}^{*}$, we have a local diffeomorphism
$\bm{x}_{2}^{*}(\bm{x}_{1}):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ such that
$D_{\bm{x}_{1}}\bm{x}_{2}=-\left(\nabla_{\bm{x}_{2}}f\right)^{-1}\nabla_{\bm{x}_{1}}f$.

To use Theorem 1, we set $f=\nabla_{x_{j}}u_{j}$ (which satisfies the
conditions of Theorem 1 by Assumption 1), $\bm{x}_{1}=x_{i}$ and
$x_{2}=\bm{x}_{j}$ (recall that $j\in\textsc{Chd}(i)$). By IFT, there exists
$\phi_{j}(x_{i})$ such that
$D_{{x_{i}}}\phi_{j}=-(\nabla_{x_{j},x_{j}}^{2}u_{j})^{-1}\nabla_{x_{j},x_{i}}^{2}u_{j}.$
Define $\nabla^{2}_{j}:=\nabla_{x_{j},x_{i}}^{2}u_{j}$. Then

$\displaystyle\left(\nabla_{\bm{x}_{L}}u_{i}\right)\left(D_{x_{i}}\phi_{L}\right)$
$\displaystyle=-\sum_{j\in\textsc{Chd}(i)}\left(\nabla_{x_{j}}u_{i}\right)D_{x_{i}}\phi_{j}$
$\displaystyle=-\sum_{j\in\textsc{Chd}(i)}\left(\nabla_{x_{j}}u_{i}\right)(\nabla_{x_{j},x_{j}}^{2}u_{j})^{-1}\nabla^{2}_{j}.$

Plugging this into Equation (2), we obtain

$\displaystyle D_{x_{i}}u_{i}$
$\displaystyle\left(x_{i},x_{\textsc{Pa}(i)},\phi_{L}\left(\bm{x}_{L-1}\right)\right)$
$\displaystyle=\nabla_{x_{i}}u_{i}-\sum_{j\in\textsc{Chd}(i)}\left(\nabla_{x_{j}}u_{i}\right)(\nabla_{x_{j},x_{j}}^{2}u_{j})^{-1}\nabla^{2}_{j}.$
(3)

For a level $l<L-1$, the total derivative of player $i\in\mathscr{N}_{l}$ in a
local best response is
$D_{x_{i}}u_{i}=\nabla_{x_{i}}u_{i}+\left(\nabla_{\bm{x}_{L}}u_{i}\right)\left(D_{x_{i}}\Phi_{l}\right),$
where

$\displaystyle D_{x_{i}}\Phi_{l}$
$\displaystyle=\left(D_{\bm{x}_{l+1}}\Phi_{l+1}\right)\left(D_{x_{i}}\bm{x}_{l+1}\right)$
$\displaystyle=\sum_{j\in\textsc{Chd}(i)}\left(D_{x_{j}}\Phi_{l+1}\right)\left(D_{x_{i}}\phi_{j}\right).$
(4)

Applying IFT, we get

$D_{x_{i}}\phi_{j}=-(\nabla_{x_{j},x_{j}}^{2}u_{j})^{-1}\nabla_{x_{j},x_{i}}^{2}u_{j},$
(5)

for $j\in\textsc{Chd}(i)$. We can apply the above procedure recursively for
$D_{\bm{x}_{l+1}}\Phi_{l+1}$ to derive the total derivative for players
$i\in\mathscr{N}_{l}$ for $l<L-1$:

$\displaystyle D_{x_{i}}u_{i}=\nabla_{x_{i}}u_{i}$
$\displaystyle+\left(\sum\limits_{j\in\textsc{Leaf}(i)}(-1)^{L-l}\nabla_{x_{j}}u_{i}\right.$
$\displaystyle\left.\prod\limits_{\begin{subarray}{c}\eta\in\textsc{Path}(j\rightarrow
i)\end{subarray}}\left(\nabla^{2}_{x_{\eta},x_{\eta}}u_{\eta}\right)^{-1}\nabla^{2}_{x_{\eta},\bm{x}_{\textsc{Pa}(\eta})}u_{\eta}\right),$
(6)

where $\textsc{Path}(j\rightarrow i)$ is an ordered list of nodes (players)
lying on the unique path from $j$ to $i$, excluding $i$. Note that Equation
(3.1) is a generalization of Equation (3) where the Path only consists of the
leaf player.

While the above derivation assumes the $\phi$ and $\Phi$ functions are local
best responses, in our algorithm in each iteration we evaluate these
functional expressions for the total derivatives _at the current joint action
profile_. This significantly reduces computational complexity and ensures that
Algorithm 1 satisfies the first-order conditions upon convergence.

###  3.2 Convergence Analysis

As we remarked earlier, stable points of DBI are not guaranteed to be SPE just
as stable points of gradient ascent are not guaranteed to be globally optimal
with general non-convex objective functions. Furthermore, DBI algorithm
entails what are effectively iterative better-response updates by players, and
it is well-known that best response dynamic processes in games will in general
lead to cycles [28].

In spite of these challenges, we provide sufficient conditions for the DBI
algorithm to converge to a stable point. In particular, in the rest of this
section, we first show that the gradient updates of DBI can be written as a
dynamical system and characterize the conditions in which this system will
converge to an stable point (Proposition 1). We then show how DBI can be tuned
(in terms of learning rate in Proposition 2, number of iterations in
Proposition 4 and initializations in Proposition 4) to converge to such stable
points when they exists. While the set of stable points and approximate SPEs
are not necessarily the same, we empirically show that DBI is effective in
converging to SPEs.

To begin, we observe that the gradient updates in DBI can be interpreted as a
discrete dynamical system, $\bm{x}^{t+1}=F(\bm{x}^{t})$, with
$F(\bm{x}^{t})=(\bm{I}+\alpha G)(\bm{x}^{t})$ where $G$ is an update gradient
vector. This discrete system can be viewed as an approximation of a continuous
limit dynamical system $\dot{\bm{x}}=G(\bm{x})$ (i.e., letting
$\alpha\rightarrow 0$). A standard solution concept for such dynamical systems
is a _locally asymptotic stable point (LASP)_.

######  Definition 1 ([13]).

A continuous (or discrete) dynamical system $\dot{\bm{x}}=G(\bm{x})$ (or
$\bm{x}^{t+1}=F(\bm{x}^{t})$) has a locally asymptotic stable point (LASP)
$\bm{x}^{*}$ if
$\exists\epsilon>0,\lim_{t\rightarrow\infty}\bm{x}^{t}=\bm{x}^{*},\forall\bm{x}^{0}\in\mathbb{B}_{\epsilon}(\bm{x}^{*})$.

There are well-known necessary and sufficient conditions for the existence of
an LASP.

######  Proposition 1 (Characterization of LASP [40, Theorem 1.2.5, Theorem
3.2.1]).

A point $\bm{x}^{*}$ is an LASP for the continuous dynamical system
$\dot{\bm{x}}=G(\bm{x})$ if $G(\bm{x}^{*})=0$ and all eigenvalues of Jacobian
matrix $\nabla_{\bm{x}}G$ at $\bm{x}^{*}$ have negative real parts.
Furthermore, for any $\bm{x}^{*}$ such that $G(\bm{x}^{*})=0$, if
$\nabla_{\bm{x}}G$ has eigenvalues with positive real parts at $\bm{x}^{*}$,
then $\bm{x}^{*}$ cannot be an LASP.

Note that an LASP of DBI is an action profile of all players that satisfies
the first-order conditions, i.e., it has the property that no player can
improve their utility through a local gradient update. While the existence of
an LASP depends on game structure, we show that under Assumption 1, and as
long as the sufficient conditions for LASP existence in Proposition 1 are
satisfied, DBI converges to LASP. We defer all the omitted proofs to the long
version.Appendix C.

######  Proposition 2.

Let $\lambda_{1},\dotsc,\lambda_{d}$ denote the eigenvalues of the updating
Jacobian $\nabla_{\bm{x}}G$ at an LASP $\bm{x}^{*}$ and define
$\lambda^{*}=\arg\max_{i\in[d]}Re(\lambda_{i})/\left\lvert\lambda_{i}\right\rvert^{2}$,
where $Re$ is the real part operator. Then with a learning rate
$\alpha<-2Re(\lambda^{*})/\left\lvert\lambda^{*}\right\rvert^{2}$, and an
initial point $\bm{x}^{0}\in\mathbb{B}_{\epsilon}(\bm{x}^{*})$ for some
$\epsilon>0$ around $\bm{x}^{*}$, DBI converges to an LASP. Specifically, if
the choice of learning rate equals $\alpha^{*}$ and the modulus of matrix
$\rho(\bm{I}+\alpha^{*}\nabla_{\bm{x}}G)=1-\kappa<1$, then the dynamics
converge to $\bm{x}^{*}$ with the rate of $O((1-\kappa/2)^{t})$.

Proposition 2 states that there exists a region such that, if the initial
point is in that region, then DBI will converge to an LASP. We next show that
if we assume first-order Lipschitzness for the update rule, then we can also
characterize the region of initial points which converge to an LASP.

######  Proposition 3.

Suppose $G$ is $L$-Lipschitz.444Formally, this means that $\exists L>0$ such
that
$\forall\bm{x},\bm{x}^{\prime}\in\mathcal{X},\|G(\bm{x})-G(\bm{x}^{\prime})\|_{2}\leq
L\|\bm{x}-\bm{x}^{\prime}\|_{2}$. Then for all
$\bm{x}^{0}\in\mathbb{B}_{\kappa/2L}(\bm{x}^{*})$, $\epsilon>0$ and after $T$
rounds of gradient update, DBI will output a point
$\bm{x}^{T}\in\mathbb{B}_{\epsilon}(\bm{x}^{*})$ as long as
$T\geq\lceil\frac{2}{\kappa}\log\left\lVert\bm{x}^{0}-\bm{x}^{*}\right\rVert/\epsilon\rceil$
where $\kappa$ is as defined in Proposition 2.

Figure 2: Convergence behaviors on (a) a $(1,1,1)$ game with 1-d actions (b) a
$(1,1,2)$ game with 1-d actions (c) a $(1,1,1)$ game with 3-d actions.

We further show that through random initialization, the probability of
reaching a saddle point is 0, which means that with probability 1, DBI
converges to an LASP in which players are playing _local_ best responses.

######  Proposition 4.

Suppose $G$ is $L$-Lipschitz. Let $\alpha<1/L$ and define the saddle points of
the dynamics $G$ as
$\mathcal{X}^{*}_{sad}=\\{\bm{x}^{*}\in\mathcal{X}\mid\bm{x}^{*}=(\bm{I}+\alpha
G)(\bm{x}^{*}),\rho((\bm{I}+\alpha\nabla_{\bm{x}}G)(\bm{x}^{*}))>1\\}$. Also
let
$\mathcal{X}^{0}_{sad}=\\{\bm{x}^{0}\in\mathcal{X}\mid\lim_{t\rightarrow\infty}(\bm{I}+\alpha
G)^{t}(\bm{x}^{0})\in\mathcal{X}^{*}_{sad}\\}$ denote the set of initial
points that converge to a saddle point. Then $\mu(\mathcal{X}^{0}_{sad})=0$,
where $\mu$ is Lebesgue measure.

While our convergence analysis does not guarantee convergence to an
approximate SPE, our experiments show that DBI is in fact quite effective in
doing so in practice.

Figure 3: Solution qualities on (a) a $(1,1,1)$ game with 1-d actions (b) a
$(1,1,2)$ game with 1-d actions (c) a $(1,1,1)$ game with 3-d actions.

##  4 Experiments

In this section, we empirically investigate the following questions: (1) the
convergence rate of DBI, (2) the solution quality of DBI, (3) the behavior of
DBI in games where we can verify global stability. All our code is written in
python. We ran our experiments on an Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
to obtain the results in Sections C.2, and on an Intel(R) Core(TM) i9-9820X
CPU @ 3.30GHz for the rest of the experiments. 555Code available at
https://github.com/jtongxin/SHG_DBI.

We evaluate the performance in terms of quality of equilibrium approximation
as a function of the number of iterations of a given algorithm, or its running
time. Ideally, given a collection of actions $\bm{x}$ played by players along
the (approximate) equilibrium path computed, we wish to find the largest
utility gain any player can have by deviating from this path, which we denote
by $\epsilon(\bm{x})$. However, this computation is impossible in our setting,
as it would need to consider all possible histories as well, whereas our
approach and alternatives only return $\bm{x}$ along the path of play
(moreover, considering all possible histories is itself intractable).

Therefore, we consider two heuristic alternatives. The first, which we call
_local SPE regret_ , runs DBI for every player $i$ starting with $\bm{x}$, and
returns the greatest benefit that any player can thereby obtain; we use this
in Section C.2. In the rest of this section, we use the second alternative,
which we call _global SPE regret_. It considers for each player $i$ in level
$l$ a discrete grid of alternative actions, and uses best response dynamics to
compute an approximate SPE of the level-$(l+1)$ subgame to evaluate player
$i$’s utility for each such deviation. This approach then returns the highest
regret among all players computed in this way.

Our evaluation considers three SHG scenarios. We begin by comparing DBI to a
number of baselines on simple, stylized SHG models, then move on to three
complex hierarchical game models motivated by concrete applications.

###  4.1 Polynomial Games

We begin by considering instances of SHGs to which we can readily apply
several state-of-the-art baselines, allowing us a direct comparison to
previous work. Specifically, we consider 3 SHG instances with different game
properties: (a) a three-level chain structure (or the $(1,1,1)$ game) with 1-d
actions (b) a “$\Yup$" shape tree (or the $(1,1,2)$ game) with 1-d action
spaces, and (c) and $(1,1,1)$ game with 3-d actions. In all the games, the
payoffs are polynomial functions of $\bm{x}$ with randomly generated
coefficients (we can think of these as proxies for a Taylor series
approximation of actual utility functions). The exact coefficient of these
polynomial functions as well as an analysis of the running time of each method
can be found in Appendix C.

We compare DBI with the following five baselines: 1) simultaneous partial
gradient ascent (SIM) [7, 28], 2) symplectic gradient dynamics with or 3)
without alignment (SYM_ALN and SYM, respectively) [4], 4) consensus
optimization (CO) [30], and 5) Hamilton gradient (HAM) [26, 1]. SIM, SYM_ALN,
SYM, CO and HAM are all designed to compute a local Nash equilibrium [4, 7].

We start by comparing convergence behavior of DBI to the baselines. We run all
algorithm with the same initial point and learning rate. The results are in
Figure 2 where we plot the $L_{2}$ norm of total gradient for each of the
algorithms (Y axis) against the number of iterations (X axis).

In all cases, DBI converges to a critical point that meets the first-order
conditions while the baseline algorithms fail to do so in most cases. In
Figures 2(a) and (c), all baselines have converged to a point with finite norm
for the total gradients. In (b), however, only CO and HAM converge to a
stationary point while SIM, SYM, SYM_ALN all diverge. For scenario (b), DBI
appears to be on an inward spiral to a critical point. We further check the
second-order condition (see Appendix C) and verify that DBI has actually
converged to local maxima of individual payoffs in all three games.

Next, we investigate solution quality in terms of _local regret_ of DBI
compared to baselines. As shown in Figure 3, across all three game instances,
DBI outputs a profile of actions (along the path of play) with near-zero local
regret while other algorithm fail to do so.

###  4.2 Decentralized Epidemic Policy Game

Next, we consider DBI for solving a class of games inspired by hierarchical
decentralized policy-making in the context of epidemics such as COVID-19 [18].
The hierarchy has levels corresponding to the (single) federal government,
multiple states, and county administrations under each state. Each player’s
action (policy) is a scalar in $[0,1]$ that represents, for example, the
extent of social distancing recommended or mandated by a player (e.g., a
state) for its administrative subordinates (e.g., counties). Crucially, these
subordinates have considerable autonomy about setting their own policies, but
incur a non-compliance cost for significantly deviating from recommendations
made by the level immediately above (of course, non-compliance costs are not
relevant for the root player). The full cost function of each player
additionally includes an infection prevalence within the geographic territory
of interest to the associated entity (e.g., within the state), as well as the
socio-economic cost of the policy itself. To summarize, the total cost for
each player is a combination of the infection cost, socio-economic cost as
well as the non-compliance cost (when applicable). However, different players
can have different combinations of these cost (through player-specific weights
for each of the costs) that can lead to strategic tensions between the players
(see Appendix Cfor details).

---
Figure 4: Global regret for the decentralized epidemic policy game. Top and
bottom columns correspond to result for games with 2 and 3 levels,
respectively.

Since the actions are in a one-dimensional compact space and the depth of the
hierarchy is at most 3, our baseline is the best response dynamics (BRD)
algorithm proposed by Jia et al. [18] (detailed in Appendix C), and we use
_global regret_ as a measure of efficacy in comparing the proposed DBI
algorithm with BRD. The results of this comparison are shown in Figures 4 and
5 for two-level (government and states) and three-level (government, states,
counties) variants of this game. We consider two-level games with 20 and 50
leaves (states), and three-level games with 2 players in level 2 (states) and
4 and 10 leaves (counties).

---
Figure 5: Running time for the decentralized epidemic policy game. Top and
bottom columns correspond to result for games with 2 and 3 levels,
respectively.

As we can see in Figure 4, BRD can have poor convergence behavior in terms of
global regret, whereas DBI appears to converge quite reliably to a path of
play with a considerably lower global regret. Notably, the improvement in
solution quality becomes more substantial as we increase the game complexity
either in terms of scale (number of leaves) or in terms of the level of
hierarchy (moving from 2- to 3-level games).

Running time (in seconds) demonstrates the relative efficacy of DBI even
further (see Figure 5). In particular, observe the significant increase in the
running time of BRD as we increase the number of leaves. In contrast, DBI is
far more scalable: indeed, even more than doubling the number of players
appears to have little impact on its running time. Moreover, BRD is several
orders of magnitude slower than DBI for the more complex games.

###  4.3 Hierarchical Public Goods Games

Figure 6: Performance ($\epsilon$) in the Public Goods Game; the scatter
points show the results of BRD with discretization factors $0.5,0.2,0.1,0.05$,
and best response rounds $2,3$.

Next, we consider hierarchical public goods games. A conventional networked
public goods game endows each player $i$ with a utility function
$u_{i}(x_{i},x_{-i})=a_{i}+b_{i}x_{i}+\sum_{j}g_{ji}x_{i}x_{j}-c_{i}(x_{i})$,
where $g_{ji}$ is the impact of player $j$ on player $i$ (often represented as
a weighted edge on a network), and $x_{i}\in[0,1]$ the level of investment in
the public good by player $i$ [6]. We construct a 3-level hierarchical variant
of such games by starting with the karate club network [43] which represents
friendships among 34 individuals. Level-2 nodes are obtained by partitioning
the network into two (sub)clubs, with leaves (level-3 nodes) representing all
the individuals. The utility of level-2 nodes is the sum of utilities of
individual members of associated clubs, with the utility of the root being the
sum of the utilities of all individuals. Furthermore, we introduce non-
compliance costs with investment policies in the level immediately above, as
we did in the decentralized epidemic policy game (Section 4.2). Further
details on the exact form of the utility functions and parameters of the games
are provided in Appendix C.4.

Figure 6 presents the global regret as a function of running time for DBI
(black line) and BRD with different levels of discretization (dots). We
observe that DBI yields considerably lower regret in these games than BRD even
as we discretize the latter finely. Moreover, DBI reaches smaller regret
orders of magnitude faster than BRD.

(a) Security game $\kappa=0.1$.

(b) Security game $\kappa=0.5$.

Figure 7: Results on $(1,3,6)$ hierarchical security games. (a) $\kappa=0.1$
and (b) $\kappa=0.5$; legend is shared.

###  4.4 Hierarchical Security Games

In the final set of experiments, we evaluate DBI on a hierarchical extension
of interdependent security games [3]. In these games, $n$ defenders can each
invest $x_{i}\geq 0$ in security. If defender $i$ is attacked, the probability
that the attack succeeds is $1/(1+x_{i})$. Furthermore, defenders are
interdependent, so that a successful attack on defender $i$ cascades to
defender $j$ with probability $q_{ji}$. In the variant we adopt, the attacker
strategy is a uniform distribution over defenders (e.g., the “attacker” is
just nature, with attacks representing stochastic exogenous failures). The
utility of the defender is the probability of surviving the attack less the
cost of security investment.

We extend this simultaneous-move game to a hierarchical structure consisting
of one root player (e.g., government), three level-2 players (e.g., sectors),
and six leaf players (e.g., organizations). The policy-makers in the first two
levels of the game recommend an investment policy to the level below, and aim
to maximize total welfare (sum of utilities) among the leaf players in their
subtrees. Just as in both hierarchical epidemic and public goods games,
whenever a player in level $l$ does not act according to the recommendation of
their parent in level $l-1$, they incur a non-compliance cost. Complete model
details are deferred to Appendix C.5. We conduct experiments with two weights
$\kappa$ that determine the relative importance of non-compliance costs in the
decisions of non-root players in the game: $\kappa\in\\{0.1,0.5\\}$.

Figures 7(a) and 7(b) present the results of comparing DBI with BRD on this
class of games, where BRD is again evaluated with different levels of action
space discretization (note, moreover, that in this setting discretizing
actions is not enough, since these are unbounded, and we also had to impose an
upper bound). We can observe that for either value of $\kappa$, DBI yields
high-quality SPE approximation (in terms of global SPE regret) far more
quickly than BRD. In particular, when we use relatively coarse discretization,
BRD is approximately an order of magnitude slower, and yields significantly
higher regret. In contrast, if we use finer discretization for BRD, global
regret for BRD and DBI becomes comparable, but now BRD is several orders of
magnitude slower. For example, DBI converges within several seconds, whereas
if we discretize $x_{i}$ into multiples of 0.02, BRD takes nearly 2 hours,
while discretization at the level of 0.01 results in BRD taking nearly 7
hours.

##  5 Conclusion

We introduced a novel class of hierarchical games, proposed a new game-
theoretic solution concept and designed an algorithm to compute it. We assume
a specific form of utility dependency between players and our solution concept
only guarantees local stability. Improvement on each of these two fronts is an
interesting direction for future work.

Given the generality of our framework, our approach can be used for many
applications characterized by a hierarchy of strategic agents e.g., pandemic
policy making. However, our modeling requires the full knowledge of the true
utility functions of all players and our analysis assumes full rationality for
all the players. Although the model we have addressed here is already
challenging, these assumptions are unlikely to hold in many real-world
applications. Therefore, further analysis is necessary to fully gauge the
robustness of our approach before deployment.

## Acknowledgments

This work was supported in part by the US Army Research Office under MURI
grant # W911NF-18-1-0208.

## References

  * Abernethy et al. [2021] Jacob Abernethy, Kevin Lai, and Andre Wibisono.  Last-iterate convergence rates for min-max optimization: convergence of Hamiltonian gradient descent and consensus optimization.  In _Thirty-Second International Conference on Algorithmic Learning Theory_ , pages 3–47, 2021.
  * Amin et al. [2016] Kareem Amin, Satinder Singh, and Michael Wellman.  Gradient methods for Stackelberg security games.  In _Thirty-Second Conference on Uncertainty in Artificial Intelligence_ , pages 2–11, 2016.
  * Bachrach et al. [2013] Yoram Bachrach, Moez Draief, and Sanjeev Goyal.  Contagion and observability in security domains.  In _2013 51st Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 1364–1371. IEEE, 2013.
  * Balduzzi et al. [2018] David Balduzzi, Sebastien Racaniere, James Martens, Jakob Foerster, Karl Tuyls, and Thore Graepel.  The mechanics of $n$-player differentiable games.  In _Thirty-Fifth International Conference on Machine Learning_ , pages 354–363, 2018.
  * Basilico et al. [2016] Nicola Basilico, Stefano Coniglio, and Nicola Gatti.  Methods for finding leader-follower equilibria with multiple followers.  In _Fifteenth International Conference on Autonomous Agents and Multiagent Systems_ , pages 1363–1364, 2016.
  * Bramoullé and Kranton [2007] Yann Bramoullé and Rachel Kranton.  Public goods in networks.  _Journal of Economic theory_ , 135(1):478–494, 2007.
  * Chasnov et al. [2020] Benjamin Chasnov, Lillian Ratliff, Eric Mazumdar, and Samuel Burden.  Convergence analysis of gradient-based learning in continuous games.  In _Thirty-Sixth Conference on Uncertainty in Artificial Intelligence_ , pages 935–944, 2020.
  * Coniglio et al. [2020] Stefano Coniglio, Nicola Gatti, and Alberto Marchesi.  Computing a pessimistic stackelberg equilibrium with multiple followers: The mixed-pure case.  _Algorithmica_ , 82(5):1189–1238, 2020.
  * Daskalakis and Panageas [2018] Constantinos Daskalakis and Ioannis Panageas.  The limit points of (optimistic) gradient descent in min-max optimization.  In _Thirty-Second International Conference on Neural Information Processing Systems_ , pages 9256–9266, 2018.
  * Dontchev and Rockafellar [2009] Asen Dontchev and Tyrrell Rockafellar.  _Implicit functions and solution mappings_.  Springer, 2009.
  * Fiez et al. [2020] Tanner Fiez, Benjamin Chasnov, and Lillian Ratliff.  Implicit learning dynamics in Stackelberg games: Equilibria characterization, convergence analysis, and empirical study.  In _Thirty-Seventh International Conference on Machine Learning_ , pages 3133–3144, 2020.
  * Fudenberg et al. [1998] Drew Fudenberg, Fudenberg Drew, David K Levine, and David K Levine.  _The theory of learning in games_ , volume 2.  MIT press, 1998.
  * Galor [2007] Oded Galor.  _Discrete dynamical systems_.  Springer Science & Business Media, 2007.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  In _Twenty-Seventh Conference on Neural Information Processing Systems_ , pages 2672–2680, 2014.
  * Hill and Varone [2021] Michael Hill and Frederic Varone.  _The public policy process_.  Routledge, 2021.
  * Horn and Johnson [2012] Roger Horn and Charles Johnson.  _Matrix analysis_.  Cambridge University Press, 2012.
  * Ibrahim et al. [2020] Adam Ibrahim, Waıss Azizian, Gauthier Gidel, and Ioannis Mitliagkas.  Linear lower bounds and conditioning of differentiable games.  In _Thirty-Seventh International Conference on Machine Learning_ , pages 4583–4593, 2020.
  * Jia et al. [2021] Feiran Jia, Aditya Mate, Zun Li, Shahin Jabbari, Mithun Chakraborty, Milind Tambe, Michael Wellman, and Yevgeniy Vorobeychik.  A game-theoretic approach for hierarchical policy-making.  _CoRR_ , abs/2102.10646, 2021.
  * Jin et al. [2020] Chi Jin, Praneeth Netrapalli, and Michael Jordan.  What is local optimality in nonconvex-nonconcave minimax optimization?  In _Thirty-Seventh International Conference on Machine Learning_ , pages 4880–4889, 2020.
  * Kearns et al. [2001] Michael Kearns, Michael Littman, and Satinder Singh.  Graphical models for game theory.  In _Seventeenth Conference on Uncertainty in Artificial Intelligence_ , pages 253–260, 2001.
  * Kelley [1955] John Kelley.  _General topology_.  Springer, 1955.
  * Letcher [2021] Alistair Letcher.  On the impossibility of global convergence in multi-loss optimization.  In _Ninth International Conference on Learning Representations_ , 2021\.
  * Leyffer and Munson [2005] Sven Leyffer and Todd Munson.  Solving multi-leader-follower games.  _Preprint ANL/MCS-P1243-0405_ , 2005.
  * Li and Başar [1987] Shu Li and Tamer Başar.  Distributed algorithms for the computation of noncooperative equilibria.  _Automatica_ , 23(4):523–533, 1987.
  * Lin et al. [2020] Tianyi Lin, Zhengyuan Zhou, Panayotis Mertikopoulos, and Michael Jordan.  Finite-time last-iterate convergence for multi-agent learning in games.  In _Thirty-Seventh International Conference on Machine Learning_ , pages 6161–6171, 2020.
  * Loizou et al. [2020] Nicolas Loizou, Hugo Berard, Alexia Jolicoeur-Martineau, Pascal Vincent, Simon Lacoste-Julien, and Ioannis Mitliagkas.  Stochastic Hamiltonian gradient methods for smooth games.  In _Thirty-Seventh International Conference on Machine Learning_ , pages 6370–6381, 2020.
  * Mas-Colell et al. [1995] Andreu Mas-Colell, Michael Dennis Whinston, Jerry R Green, et al.  _Microeconomic theory_ , volume 1.  Oxford university press New York, 1995.
  * Mazumdar et al. [2020] Eric Mazumdar, Lillian Ratliff, and Shankar Sastry.  On gradient-based learning in continuous games.  _SIAM Journal on Mathematics of Data Science_ , 2(1):103–131, 2020.
  * Mertikopoulos and Zhou [2019] Panayotis Mertikopoulos and Zhengyuan Zhou.  Learning in games with continuous action sets and unknown payoff functions.  _Mathematical Programming_ , 173(1):465–507, 2019\.
  * Mescheder et al. [2017] Lars Mescheder, Sebastian Nowozin, and Andreas Geiger.  The numerics of GANs.  In _Thirty-First International Conference on Neural Information Processing Systems_ , pages 1823–1833, 2017.
  * Nagarajan and Kolter [2017] Vaishnavh Nagarajan and Zico Kolter.  Gradient descent GAN optimization is locally stable.  In _Thirty-First International Conference on Neural Information Processing Systems_ , pages 5591–5600, 2017.
  * Nguyen et al. [2020] Thanh Nguyen, Arunesh Sinha, and He He.  Partial adversarial behavior deception in security games.  In _Twenty-Nineth International Joint Conference on Artificial Intelligence_ , pages 283–289, 2020.
  * Osborne and Rubinstein [2004] Martin J Osborne and Ariel Rubinstein.  _An introduction to game theory_.  Oxford university press New York, 2004.
  * Selten [1975] R Selten.  Reexamination of the perfectness concept for equilibrium points in extensive games.  _International Journal of Game Theory_ , 4(1):25–55, 1975.
  * Shaban et al. [2019] Amirreza Shaban, Ching-An Cheng, Nathan Hatch, and Byron Boots.  Truncated back-propagation for bilevel optimization.  In _Twenty-Second International Conference on Artificial Intelligence and Statistics_ , pages 1723–1732, 2019.
  * Shub [1987] Michael Shub.  _Global stability of dynamical systems_.  Springer, 1987.
  * Von Stackelberg [1952] Heinrich Von Stackelberg.  _The theory of the market economy_.  Oxford University Press, 1952.
  * Wang et al. [2021] Kai Wang, Lily Xu, Andrew Perrault, Michael Reiter, and Milind Tambe.  Coordinating followers to reach better equilibria: End-to-end gradient descent for stackelberg games.  _CoRR_ , abs/2106.03278, 2021.
  * Wang et al. [2020] Yuanhao Wang, Guodong Zhang, and Jimmy Ba.  On solving minimax optimization locally: A follow-the-ridge approach.  In _Eighth International Conference on Learning Representations_ , 2020.
  * Wiggins [2003] Stephen Wiggins.  _Introduction to applied nonlinear dynamical systems and chaos_.  Springer Science & Business Media, 2003.
  * Wilder et al. [2020] Bryan Wilder, Marie Charpignon, Jackson Killian, Han-Ching Ou, Aditya Mate, Shahin Jabbari, Andrew Perrault, Angel Desai, Milind Tambe, and Maimuna Majumder.  Modeling between-population variation in COVID-19 dynamics in Hubei, Lombardy, and New York City.  _Proceedings of the National Academy of Sciences_ , 117(41):25904–25910, 2020.
  * Wolfram [1999] Stephen Wolfram.  _The Mathematica_.  Cambridge University Press, 1999.
  * Zachary [1977] Wayne W Zachary.  An information flow model for conflict and fission in small groups.  _Journal of anthropological research_ , 33(4):452–473, 1977.

##  Appendix A Omitted Details from Section 2

######  Proposition 5.

The SPE notion we defined for an SHG model corresponds exactly to the SPE of
its extensive-form game (EFG) representation.

###### Proof.

First we show how to construct the EFG representation of an SHG. Since within
level $l$ players act simultaneously, we can designate a canonical player
order for player in level $l$, say we label them as
$(l,1),(l,2),\dotsc,(l,n_{l})$. Then in the EFG representation, we
sequentialize the playing order as always be
$(1,1),(2,1),(2,2),\dotsc,(3,1),\dotsc,(L,1),\dotsc,(L,n_{L})$. For every
$(l,1),l=1,\dotsc,L$, its information sets are all in the form of some
singleton $\\{h_{<l}\\}$ for some history $h_{<l}$. While for $(l,n_{l}),l>1$,
its information sets are sets that consist of some $h_{<l}$ concatenated with
all possible moves of $(l,1),\dotsc,(l,n_{l-1})$. Then according the the
classical definition of the subgame notion [34, 27], the only well-defined
subgames are those rooted at $(l,1),\forall l$ since these are the only ones
with singleton information set. So every level-$l$ subgames corresponds to a
subgame in the EFG representation, and every subgame in the EFG representation
maps to a level-$l$ subgame of the EFG by investigating the player index. Then
the SPE for both representations are also equivalent. ∎

##  Appendix B Omitted Details from Section 3.2

###### Proof of Proposition 2.

The learning dynamics of DBI can be written as $\bm{x}^{t}=(\bm{I}+\alpha
G)(\bm{x}^{t-1})$. Since $G$ has eigenvalues $\lambda_{1},\dotsc,\lambda_{d}$
at $\bm{x}^{*}$, the matrix $\bm{I}+\alpha\nabla_{\bm{x}}G$ at a stationary
point $\bm{x}^{*}$ has eigenvalues
$1+\alpha\lambda_{1},\dotsc,1+\alpha\lambda_{d}$. Since the set of LASPs is
non-empty, $Re(\lambda_{i})<0$ for all $\lambda_{i}$. Then for the choice of
$\alpha^{*}$ in the Proposition, the modulus of the Jacobian

$\rho(\bm{I}+\alpha^{*}\nabla_{\bm{x}}G)=\max_{i\in[d]}\left\lvert
1+\alpha^{*}\lambda_{i}\right\rvert=\max_{i\in[d]}\sqrt{1+2\alpha^{*}Re(\lambda_{i})+(\alpha^{*})^{2}\left\lvert\lambda_{i}\right\rvert^{2}}<1.$

Proofs for the convergence rate of $O((1-\kappa/2)^{t})$ can be found in
previous work (e.g., Fiez et al. [11, Proposition F.1] or Wang et al. [39,
Proposition 4])). For completeness, we provide a proof. Since
$\rho((\bm{I}+\alpha\nabla_{\bm{x}}G)(\bm{x}^{*}))=1-\kappa$, according to
Horn and Johnson [16, Lemma 5.6.10], there exists a matrix norm
$\left\lVert\cdot\right\rVert$ such that
$\left\lVert\bm{I}+\alpha\nabla_{\bm{x}}G\right\rVert<1-\kappa+\epsilon$, for
$\forall\epsilon>0$. We choose $\epsilon=\frac{\kappa}{4}$ The Taylor
expansion of $(\bm{I}+\alpha G)$ at $\bm{x}^{*}$ is

$\displaystyle(\bm{I}+\alpha G)(\bm{x})=(\bm{I}+\alpha
G)(\bm{x}^{*})+(\bm{I}+\alpha\nabla_{\bm{x}}G)(\bm{x}^{*})(\bm{x}-\bm{x}^{*})+R(\bm{x}-\bm{x}^{*}),$

where $R(\bm{x}-\bm{x}^{*})=o(\left\lVert\bm{x}-\bm{x}^{*}\right\rVert)$. Let
$R_{1}(\bm{x}-\bm{x}^{*})=\frac{1}{\alpha}R(\bm{x}-\bm{x}^{*})$, Then we have
$\lim_{\bm{x}\rightarrow\bm{x}^{*}}\frac{R_{1}(\bm{x}-\bm{x}^{*})}{\left\lVert\bm{x}-\bm{x}^{*}\right\rVert}=0$.
Then we can choose $\delta>0$ such that $\left\lVert
R_{1}(\bm{x}-\bm{x}^{*})\right\rVert\leq\tfrac{\kappa}{4}\left\lVert\bm{x}-\bm{x}^{*}\right\rVert$
when $\left\lVert\bm{x}-\bm{x}^{*}\right\rVert<\delta$.

$\displaystyle\left\lVert G(\bm{x})-G(\bm{x}^{*})\right\rVert$
$\displaystyle\leq\left\lVert\nabla_{\bm{x}}G(\bm{x}^{*})(\bm{x}-\bm{x}^{*})\right\rVert+\left\lVert
R_{1}(\bm{x}-\bm{x}^{*})\right\rVert$
$\displaystyle\leq\left\lVert\nabla_{\bm{x}}G(\bm{x}^{*})\right\rVert\left\lVert(\bm{x}-\bm{x}^{*})\right\rVert+\frac{\kappa}{4}\left\lVert\bm{x}-\bm{x}^{*}\right\rVert$
$\displaystyle\leq(1-\frac{\kappa}{2})\left\lVert\bm{x}-\bm{x}^{*}\right\rVert.$

This shows the operator $(\bm{I}+\alpha G)$ is a contract mapping with
contraction constant $(1-\tfrac{\kappa}{2})$. Therefore the convergence rate
is $O((1-\kappa/2)^{t})$ ∎

Before we present the proof of Proposition 4, we state the following lemma.

######  Lemma 1.

The update gradient vector $G$ is L-Lipchitz if and only if
$\left\lVert\nabla_{\bm{x}}G\right\rVert\leq L$ at all $\bm{x}\in\mathcal{X}$.

###### Proof.

First we prove for "if" direction. Consider
$\bm{x}_{1},\bm{x}_{2}\in\mathcal{X}$,

$\displaystyle G(\bm{x}_{2})-G(\bm{x}_{1})$
$\displaystyle=\int_{0}^{1}\nabla_{\bm{x}}G(x_{1}+\tau(\bm{x}_{2}-\bm{x}_{1}))(\bm{x}_{2}-\bm{x}_{1})d\tau$
$\displaystyle\Rightarrow\left\lVert G(\bm{x}_{2})-G(\bm{x}_{1})\right\rVert$
$\displaystyle=\left\lVert\int_{0}^{1}\nabla_{\bm{x}}G(x_{1}+\tau(\bm{x}_{2}-\bm{x}_{1}))(\bm{x}_{2}-\bm{x}_{1})d\tau\right\rVert$
$\displaystyle\leq\left\lVert\int_{0}^{1}\nabla_{\bm{x}}G(x_{1}+\tau(\bm{x}_{2}-\bm{x}_{1}))d\tau\right\rVert\left\lVert\bm{x}_{2}-\bm{x}_{1}\right\rVert$
$\displaystyle\leq\int_{0}^{1}\left\lVert\nabla_{\bm{x}}G(x_{1}+\tau(\bm{x}_{2}-\bm{x}_{1}))\right\rVert
d\tau\left\lVert\bm{x}_{2}-\bm{x}_{1}\right\rVert$ $\displaystyle\leq
L\left\lVert\bm{x}_{2}-\bm{x}_{1}\right\rVert.$

Then we prove for the "only if" direction. Take $\epsilon>0$, then for
$\forall\bm{x}_{1},\bm{x}_{2}\in\mathcal{X}$,

$\displaystyle\left\lVert\int_{0}^{\epsilon}\left(\nabla_{\bm{x}}G(\bm{x}_{1}+\tau\bm{x}_{2})\right)\cdot\bm{x}_{2}d\tau\right\rVert$
$\displaystyle=\left\lVert
G(\bm{x}_{1}+\epsilon\bm{x}_{2})-G(\bm{x}_{1})\right\rVert\leq\epsilon
L\left\lVert\bm{x}_{2}\right\rVert$
$\displaystyle\Rightarrow\lim_{\epsilon\rightarrow
0}\frac{\left\lVert\int_{0}^{\epsilon}\left(\nabla_{\bm{x}}G(\bm{x}_{1}+\tau\bm{x}_{2})\right)\cdot\bm{x}_{2}d\tau\right\rVert}{\epsilon\left\lVert\bm{x}_{2}\right\rVert}$
$\displaystyle=\frac{\left\lVert\nabla_{\bm{x}}G(\bm{x}_{1})\cdot\bm{x}_{2}\right\rVert}{\left\lVert\bm{x}_{2}\right\rVert}\leq
L.$

Since it holds for any $\bm{x}_{2}$, it must be
$\left\lVert\nabla_{\bm{x}}G(\bm{x}_{1})\right\rVert\leq L$. And since it
applies for any $\bm{x}_{1}\in\mathcal{X}$, this completes the proof. ∎

###### Proof of Proposition 4.

From the proof from Proposition 2 notice that in order to make
$\left\lVert\bm{x}^{t}-\bm{x}^{*}\right\rVert\leq\epsilon$ we only have to let
$t\geq\lceil\frac{2}{\kappa}\log\left\lVert\bm{x}^{0}-\bm{x}^{*}\right\rVert/\epsilon\rceil$
since now
$\left\lVert\bm{x}^{t}-\bm{x}^{*}\right\rVert\leq(1-\frac{\kappa}{2})^{t}\left\lVert\bm{x}^{0}-\bm{x}^{*}\right\rVert\leq\exp\left({-\kappa/2t}\right)\leq\epsilon$.
Now we need to characterize the region of initial points that can converge to
$\bm{x}^{*}$ by characterizing the maximum possible radius of such initial
point to $\bm{x}^{*}$. Recall in the proof of Proposition 2, this is captured
by the. parameter $\delta$. On one hand by using the Lipschitzness, we can
bound the residual function by

$\displaystyle\left\lVert
R_{1}(\bm{x}-\bm{x}^{*})\right\rVert\leq\int^{1}_{0}\left\lVert\bm{I}+\alpha\nabla_{\bm{x}}G(\bm{x}^{*}+\tau(\bm{x}-\bm{x}^{*}))-(\bm{I}+\alpha\nabla_{\bm{x}}G(\bm{x}^{*}))\right\rVert\left\lVert\bm{x}-\bm{x}^{*}\right\rVert
d\tau\leq\frac{L}{2}\left\lVert\bm{x}-\bm{x}^{*}\right\rVert^{2}.$

One the other hand to maintain this convergence rate we should let
$\left\lVert
R_{1}(\bm{x}-\bm{x}^{*})\right\rVert\leq\tfrac{\kappa}{4}\left\lVert\bm{x}-\bm{x}^{*}\right\rVert$.
Then we simply let
$\leq\tfrac{L}{2}\left\lVert\bm{x}-\bm{x}^{*}\right\rVert^{2}\leq\tfrac{\kappa}{4}\left\lVert\bm{x}-\bm{x}^{*}\right\rVert$
we get an initial point should satisfy
$\left\lVert\bm{x}-\bm{x}^{*}\right\rVert\leq\tfrac{\kappa}{2L}$. ∎

To prove Proposition 4, we need a few more machinery.

######  Proposition 6.

With $\alpha<1/L$, $\bm{I}+\alpha G$ is a diffeomorphism.

###### Proof.

First we show $(\bm{I}+\alpha G)$ is invertible. Suppose
$\bm{x}_{1},\bm{x}_{2}\in\mathcal{X}$ such that $\bm{x}_{1}\not=\bm{x}_{2}$
and $(\bm{I}+\alpha G)(\bm{x}_{1})=(\bm{I}+\alpha G)(\bm{x}_{2})$. Then
$\bm{x}_{1}-\bm{x}_{2}=\alpha(G(\bm{x}_{2})-G(\bm{x}_{1}))$. And by
$\alpha<1/L$ we have $\left\lVert\bm{x}_{1}-\bm{x}_{2}\right\rVert\leq\alpha
L\left\lVert\bm{x}_{1}-\bm{x}_{2}\right\rVert<\left\lVert\bm{x}_{1}-\bm{x}_{2}\right\rVert$,
which is a contraction.

Next we show we show its invert function is well-defined on any point of
$\mathcal{X}$. Notice that
$\rho(\alpha\nabla_{\bm{x}}G)\leq\left\lVert\alpha\nabla_{\bm{x}}G\right\rVert\leq\alpha
L<1$. And notice that since the eigenvalues of $\bm{I}+\alpha\nabla_{\bm{x}}G$
is the eigenvalues of $\alpha\nabla_{\bm{x}}G$ plus 1, the only way that makes
$\det(\bm{I}+\alpha\nabla_{\bm{x}}G)=0$ is to have one of the eigenvalues of
$\alpha\nabla_{\bm{x}}$ to be $-1$. This contradicts to
$\rho(\alpha\nabla_{\bm{x}}G)<1$. Therefore by the implicit function theorem,
$(\bm{I}+\alpha G)$ is a local diffeomorphism on any point of $\mathcal{X}$,
and therefore $(\bm{I}+\alpha G)^{-1}$ is well defined on $\mathcal{X}$. ∎

######  Theorem 2 (Center and Stable Manifold [36, Theorem III.7, Chapter 5]).

Suppose $\bm{x}^{*}=h(\bm{x}^{*})$ is a critical point for the $C^{r}$ local
diffeomorphism $h:\mathcal{X}\rightarrow\mathcal{X}$. Let
$\mathcal{X}=\mathcal{X}_{s}\oplus\mathcal{X}_{u}$, where $\mathcal{X}_{s}$ is
the stable center eigenspace belonging to those eigenvalues of
$\nabla_{\bm{x}}h(\bm{x}^{*})$ whose modulus is no greater than 1, and
$\mathcal{X}_{s}$ is the unstable eigenspace belonging to those whose modulus
is greater than 1. Then there exists a $C^{r}$ embeded disk
$W^{cs}_{loc}(\bm{x}^{*})$ that is tangent to $\mathcal{X}_{s}$ at
$\bm{x}^{*}$ called the local stable center manifold. Moreover there
$\exists\epsilon>0$,
$h(W^{cs}_{loc}(\bm{x}^{*}))\bigcap\mathbb{B}_{\epsilon}(\bm{x}^{*})\subset
W^{cs}_{loc}$ and
$\bigcap_{t=0}^{\infty}h^{-t}(\mathbb{B}_{\epsilon}(\bm{x}^{*}))\subset
W^{cs}_{loc}(\bm{x}^{*}).$

###### Proof of Proposition 4.

For $\forall\bm{x}^{*}\in\mathcal{X}^{*}_{sad}$, let $\epsilon(\bm{x}^{*})>0$
be the radius of neighborhood provided by Theorem 2 for diffeomorphism
$h=\bm{I}+\alpha G$ and point $\bm{x}^{*}$. Then define
$\mathscr{B}=\bigcup_{\bm{x}^{*}\in\mathcal{X}^{*}_{sad}}\mathbb{B}_{\epsilon(\bm{x}^{*})}(\bm{x}^{*})$.
And since $\mathcal{X}$ is a subset of Euclidean space so it is second-
countable. By Lindelöf’s Lemma [21], which stated that every open cover there
is a countable subcover, we can actually write
$\mathscr{B}=\bigcup_{i=1}^{\infty}\mathbb{B}_{\epsilon(\bm{x}^{*}_{i})}(\bm{x}^{*}_{i})$
for a countable family of saddle points
$\\{\bm{x}^{*}_{i}\\}_{i=1}^{\infty}\subseteq\mathcal{X}^{*}_{sad}$. Therefore
for $\forall\bm{x}^{0}\in\mathcal{X}^{0}_{sad}$ that converges to a saddle
point, it must converge to $\bm{x}^{*}_{i}$ for some $i$. And $\exists
T(\bm{x}^{0})>0$, such that $\forall
t>T(\bm{x}^{0}),h^{t}(\bm{x}^{0})\in\mathbb{B}_{\epsilon(\bm{x}^{*}_{i})}(\bm{x}^{*}_{i})$.
From Theorem 2 we have $h^{t}(\bm{x}^{0})\in W^{sc}_{loc}(\bm{x}^{*}_{i})$.
Since $h$ is a diffeomorphism on $\mathcal{X}$ we have $\bm{x}^{0}\in
h^{-t}(W^{sc}_{loc}(\bm{x}_{i}^{*}))$. We furthermore union over all finite
time step
$\bm{x}^{0}\in\bigcup_{t=0}^{\infty}h^{-t}(W^{sc}_{loc}(\bm{x}_{i}^{*}))$.
Then we have
$\mathcal{X}^{0}_{sad}\subseteq\bigcup_{i=1}^{\infty}\bigcup_{t=0}^{\infty}h^{-t}(W^{sc}_{loc}(\bm{x}_{i}^{*}))$.
For each $i$ since $\bm{x}^{*}_{i}$ is a saddle point, it has an eigenvalue
greater than 1, so the dimension of unstable eigenspace
$dim(\mathcal{X}_{u}(\bm{x}^{*}_{i}))\geq 1$. Therefore the dimension of
$W^{sc}_{loc}(\bm{x}_{i}^{*})$ is less than full dimension. This leads to
$\mu(W^{sc}_{loc}(\bm{x}_{i}^{*}))=0$ for any $i$. Again since $h$ is a
diffeomorphism, $h^{-1}$ is locally Lipschitz which are null set preserving.
Then $\mu(h^{-t}(W^{sc}_{loc}(\bm{x}_{i}^{*})))=0$ for $\forall i,t$. And
since the countable union of measure zero sets is measure zero,
$\mu(\bigcup_{i=1}^{\infty}\bigcup_{t=0}^{\infty}h^{-t}(W^{sc}_{loc}(\bm{x}_{i}^{*})))=0$.
So $\mu(\mathcal{X}^{0}_{sad})=0$. ∎

##  Appendix C Omitted Details from Section 4

###  C.1 Derivation of Total Second-Order Derivatives

In our experiments we may need to check the second-order total derivative at a
point. This section provides our approach to calculate such quantity.

For player of level $L$,
$D^{2}_{\bm{x}_{l,i},\bm{x}_{l,i}}u_{l,i}=\nabla^{2}_{\bm{x}_{l,i},\bm{x}_{l,i}}u_{l,i}$.

For level $l<L$, notice that $D_{\bm{x}_{l,i}}u_{l,i}$ is an explicit function
of the action variables of players that are descendants of $(l,i)$. Denote
this set of players as $\mathit{DES}(l,i)$ (excluding $(l,i)$), then we have
$D_{\bm{x}_{l,i},\bm{x}_{l,i}}^{2}u_{l,i}=\nabla_{\bm{x}_{l,i}}(D_{\bm{x}_{l,i}}u_{l,i})+\sum_{(\lambda,\eta)\in\mathit{DES}(l,i)}\nabla_{\bm{x}_{\lambda,\eta}}(D_{\bm{x}_{l,i}}u_{l,i})D_{\bm{x}_{l,i}}\bm{x}_{\lambda,\eta}$.

Here
$D_{\bm{x}_{l,i}}\bm{x}_{\lambda,\eta}=(-1)^{\lambda-l}\prod\limits_{\begin{subarray}{c}(\lambda^{\prime},\eta^{\prime})\in\\\
\mathit{PATH}((\lambda,\eta)\\\
\rightarrow(l,i))\end{subarray}}\left(\nabla^{2}_{\bm{x}_{\lambda^{\prime},\eta^{\prime}},\bm{x}_{\lambda^{\prime},\eta^{\prime}}}u_{\lambda^{\prime},\eta^{\prime}}\right)^{-1}\nabla^{2}_{\bm{x}_{\lambda^{\prime},\eta^{\prime}},\bm{x}_{\mathit{PA}(\lambda^{\prime},\eta^{\prime}})}u_{\lambda^{\prime},\eta^{\prime}}$.

The terms $\nabla_{\bm{x}_{l,i}}(D_{\bm{x}_{l,i}}u_{l,i})$ and
$\nabla_{\bm{x}_{\lambda,\eta}}(D_{\bm{x}_{l,i}}u_{l,i})$ involve thrice
derivatives. In practice we use either numerical finite-difference methods to
approximate this value or use symbolic algebra feature in software tools such
as Mathematica [42] to calculate the analytical form of this derivative.

###  C.2 Details of Polynomial Games Experiments

#### Setup

We compare the following algorithms, such that the gradient update vectors are
as follows

DBI: $G^{\mathit{DBI}}=(D_{\bm{x}_{1}}u_{1},\dotsc,D_{\bm{x}_{n}}u_{n})$. SIM:
$G^{\mathit{SIM}}=(\nabla_{\bm{x}_{1}}u_{1},\dotsc,\nabla_{\bm{x}_{n}}u_{n})$.
SYM: $G^{\mathit{SYM}}=G^{\mathit{SIM}}+G^{\mathit{SIM}}A^{\mathit{SIM}}$,
where $A^{\mathit{SIM}}=(J^{\mathit{SIM}}-(J^{\mathit{SIM}})^{T})/2$, and
$J^{\mathit{SIM}}$ is the Jacobian matrix of $G^{\mathit{SIM}}$. HAM:
$G^{\mathit{HAM}}=-\nabla_{\bm{x}}\left\lVert
G^{\mathit{SIM}}\right\rVert^{2}$. CO:
$G^{\mathit{CO}}=G^{\mathit{SIM}}+\gamma G^{\mathit{HAM}}$, where we let
$\gamma=0.1$ for all experiments. SYM_ALN:
$G^{\mathit{SYM}}=G^{\mathit{SIM}}+\zeta G^{\mathit{SIM}}A^{\mathit{SIM}}$,
where $\zeta$ is the sign of
$\frac{1}{2d}(G^{\mathit{HAM}}(G^{\mathit{SIM}})^{T})(G^{\mathit{HAM}}(G^{\mathit{SIM}}A^{\mathit{SIM}})^{T})+0.1$.

#### Experimental Details for (1, 1, 1), 1-d Action Games

We use $u_{1}(x,z)=-7x^{2}+9xz+x-z$,
$u_{2}(x,y,z)=-2y^{2}-4yz-10x^{2}+2xz-3z^{2}+4y+7x-8z-8xyz$, and,
$u_{3}(y,z)=-10z^{2}-9yz+9y^{2}-5z-2y$. Here $(x,y,z)$ are action variables
for players 1, 2, 3. We use the learning rate of $\alpha=1e-5$ for all
algorithms. We verify that DBI converges to an LSPE
$(x,y,z)=(-0.34,1.85,-1.08)$.

#### Experimental Details for (1, 1, 2), 1-d Action Games

We use $u_{1}(x,y,z)=-2x^{2}-3xy+y^{2}+5x+7y+3xz-10yz+5xyz-6z$,
$u_{2}(w,x,y,z)=2w^{2}-wx-3wy-5x^{2}+9xy+2y^{2}+3w+5x-4y+5z^{2}+8wz+7xz-9yz-10z$,
$u_{3}(w,y,z)=-5y^{2}-8yz+z^{2}+8y-9z-2wy-4wz-w^{2}-8wyz-2w$, and,
$u_{4}(w,y,z)=-10z^{2}-2yz+5y^{2}-7z-6y-3wz-8wy-10wyz+5w$ Here $(x,w,y,z)$ are
action variables for players 1, 2, 3, 4. We use $\alpha=4e-6$ for all
algorithms. We verify that DBI converges to an LSPE
$(x,w,y,z)=(4.70,-2.13,10.27,9.93)$.

#### Experimental Details for (1, 1, 1), 3-d Action Games

We use
$u_{1}(\bm{x},\bm{z})=-7\sum_{i=1}^{3}x_{i}^{2}+9(\sum_{i=1}^{3}x_{i})(\sum_{i=1}^{3}z_{i})+\sum_{i=1}^{3}x_{i}-\sum_{i=1}^{3}z_{i}$,
$u_{2}(\bm{x},\bm{y},\bm{z})=-2\sum_{i=1}^{3}y_{i}^{2}-4(\sum_{i=1}^{3}y_{i})(\sum_{i=1}^{3}z_{i})-10(\sum_{i=1}^{3}x_{i})(\sum_{i=1}^{3}z_{i})+2\sum_{i=1}^{3}x_{i}^{2}-3\sum_{i=1}^{3}z_{i}^{2}+4(\sum_{i=1}^{3}x_{i})(\sum_{i=1}^{3}y_{i})(\sum_{i=1}^{3}z_{i})+7(\sum_{i=1}^{3}x_{i})-8(\sum_{i=1}^{3}y_{i})-8(\sum_{i=1}^{3}z_{i})$,
and,
$u_{3}(\bm{y},\bm{z})=-10\sum_{i=1}^{3}z_{i}^{2}-9(\sum_{i=1}^{3}y_{i})(\sum_{i=1}^{3}z_{i})+9\sum_{i=1}^{3}y_{i}^{2}-5\sum_{i=1}^{3}y_{i}-2\sum_{i=1}^{3}z_{i}$,
where $\bm{x},\bm{y},\bm{z}$ are action variables of players 1, 2, 3. We use
the learning rate $\alpha=1e-5$ for all algorithms. We verify that DBI
converges to an LSPE
$(\bm{x},\bm{y},\bm{z})=((-0.39,-0.39,-0.39),(0.29,0.29,0.29),(-0.58,-0.58,-0.58))$.

The time cost of each algorithm is measured in Table (1).

| DBI | SIM | SIM_ALN | SYM | CO | HAM
---|---|---|---|---|---|---
(1,1,1)-1d | 1.65 | 1.26 | 5.92 | 1.61 | 4.91 | 4.58
(1,1,2)-1d | 2.93 | 2.07 | 11.9 | 3.57 | 14.9 | 13.3
(1,1,1)-3d | 12.4 | 10.9 | 60.1 | 9.63 | 49.4 | 53.9
Table 1: Comparison of the running times of DBI and baselines from experiment
1.

#### A note on the physical significance of the above utility functions

Although our methodology is applicable in principle to any set of agents’
utility function satisfying the definition of SHGs, it is reasonable to ask
whether there exists a game model, grounded in some real-world domain, that
exhibits the SHG hierarchy as well as the polynomial multi-variable utility
structure above. To that end, consider the COVID-19 Game delineated below. In
the experiments described below, each agent (a policy-maker in a hierarchy
such as the federal government, state governments, and county governments) has
a non-dimensional action restricted to closed interval $[0,1]$, interpreted as
a social distancing factor due to policy intervention (either recommended or
implemented). However, in general, a policy intervention (hence, and agent’s
action) may take the form of a multi-dimensional real vector, e.g. factor for
social distancing, factor for vaccination roll-out, etc.; hence, in general,
the overall cost (negative payoff) of each agent is a function of $2+n_{L}$
vectors rather than $2+n_{L}$ scalars (correspoding to the agent, its parent,
and all leaves) as in Equation (7). Although the actual functional form is not
necessarily polynomial, we can always take a polynomial approximation in our
model via curve-fitting or a Taylor series expansion. The above polynomial
functions would then belong to the resulting class of functions; in our
experiments for Convergence Analysis experiments are randomly chosen since the
focus of this paper is not on modeling any real-world multi-agent interaction
scenario to a high degree of accuracy but on assessing the performance of our
algorithm on reasonable models consistent with the SHG structure. This line of
thinking can be also applied to the utility functions used in the experiments
on the relationship between LASP and LSPE described at the end of the
appendix.

#### Additional Experiments on the Relationship Between LASP and LSPE

In this additional experiment, We will study the chance of having a convergent
algorithm as well as the probability of finding an equilibrium upon
convergence on several classes of SHGs. We design game classes of SHGs in a
way that every game class $\mathcal{F}$ has the same parameter space on both
its topology and payoff structure. For each class, we generate $N$ instances,
and for each instance, we first find the set of critical points by numerically
solving the first order condition $D_{\bm{x}_{i}}u_{i}=0,\forall(i)$. For each
of these critical points, we determine whether it is an LASP of DBI by
checking whether the Jacobian of updating gradient $G$ have eigenvalues all of
which have negative real part. If this is the case, we classify the critical
point as an LASP according to Proposition 1. Finally, for each of the LASPs,
we further check whether it is an LSPE by checking whether
$D_{\bm{x}_{i},\bm{x}_{i}}^{2}u_{i}\prec 0,\forall(i)$. Finally, across the
$N$ instances we compute the fraction of games for which an LASP was found (%
LASP), and for those with an LASP, the fractions of the instances where an
LASP is also an LSPE (% LSPE). We call these two numbers the measure
properties of a game class.

We denote by our game classes by $\mathcal{F}_{\text{sub}}^{\text{super}}$
where the subscript is a vector listing the number of players in each SHG
level from top to bottom, and the superscript indicates the parameters of the
game. For example $\mathcal{F}^{C}_{1,1}$ is an SHG class with two levels and
one player in each level (a Stackelberg game) with $C\in\mathbb{R}^{+}$
determining payoffs as follows: $u_{i}=\sum_{\alpha+\beta\leq
4,\alpha,\beta\in\mathbb{N}}c_{i,\alpha,\beta}x^{\alpha}y^{\beta}$, where
$x,y\in\mathbb{R}$ are 1-d action variables for the player ($i=1$) and the
second level-$2$ player ($i=2$), respectively; the player-specific
coefficients $c_{i,\alpha,\beta}$ are integers generated uniformly in $[-C,C]$
for different values of $C$. When $C=\infty$, we generate the coefficients
from a continuous uniform distribution in $[-1,1]$. See Appendix C for more
information (e.g., the utility structure) of the other game classes.

$\mathcal{F}$ |  $\mathcal{F}^{1}_{1,1}$ |  $\mathcal{F}^{10}_{1,1}$ |  $\mathcal{F}^{\infty}_{1,1}$ |  $\mathcal{F}^{1}_{1,2}$ |  $\mathcal{F}^{10}_{1,2}$ |  $\mathcal{F}^{\infty}_{1,2}$ |  $\mathcal{F}^{1}_{1,1,1}$ |  $\mathcal{F}^{10}_{1,1,1}$ |  $\mathcal{F}^{\infty}_{1,1,1}$ |  $\mathcal{F}^{1}_{1,1,2}$ |  $\mathcal{F}^{10}_{1,1,2}$ |  $\mathcal{F}^{\infty}_{1,1,2}$
---|---|---|---|---|---|---|---|---|---|---|---|---
% LASP |  52.6 |  57.3 |  59.7 |  32.0 |  35.0 |  36.8 |  47.0 |  51.0 |  50.3 |  13.6 |  21.3 |  21.2
% LSPE |  88.8 |  87.9 |  87.1 |  65.6 |  64.3 |  60.0 |  66.5 |  64.3 |  67.0 |  62.1 |  60.7 |  60.2
Table 2: Results for $\mathcal{F}^{C}_{1,1}$ $\mathcal{F}^{C}_{1,2}$,
$\mathcal{F}^{C}_{1,1,1}$ and $\mathcal{F}^{C}_{1,1,2}$ averaged over
$N=10^{5}$ instances.

The results are shown in Table 2. First, we notice for the same game topology,
$C$ does not appear to substantially affect the measure property. Parameter
$C$ essentially controls the granularity of a uniform discrete distribution
and approaches a uniform continuous distribution as it becomes large. The
measure properties will also be quite similar as $C$ grows. The second
observation is that as the topology becomes much more complex, it is less
probable for DBI to find a stable point. The relevant probabilities degrade
from $52\%\sim 59\%$ for Stackelberg games to $13\%\sim 21\%$ for 4-player
3-level games. This suggests a limitation of our algorithm in facing complex
game topologies with more intricate back-propagation. However, we observe that
the probability of an LASP being an LSPE does not decay speedily, as it
achieves $87\%\sim 88\%$ for the structure of $(1,1)$, $60\%\sim 67\%$ for
$(1,2)$, $(1,1,2)$ and $(1,1,1)$.

###  C.3 Details of COVID-19 Game Experiments

#### Subgame Perfect Equilibrium in an SHG

We here elaborate how we evaluate the regret of a profile. First let’s imagine
when solving a two-level games with $n_{2}$ players at the second-level. Given
$(1,1)$’s action choice, the problem of computing SPE in the second-level is
equivalent to computing a Nash equilibrium in the second-level. However, an
exact pure Nash equilibrium may not exist. So which profile should we choose
to propagate back to $(1,1)$? To resolve this issue, we define
$\varepsilon$-Nash equilibrium.

######  Definition 2.

For a simultaneous-move game, a profile $\bm{x}^{*}$ is an $\varepsilon$-Nash
if for any player $n$,
$\forall\bm{x}_{n}^{\prime}\in\mathcal{X}_{n},u_{n}(\bm{x}_{n}^{\prime},\bm{x}^{*}_{-n})\leq
u_{n}(\bm{x}^{*})+\varepsilon$.

In another word, an $\varepsilon$-Nash is a profile where for every player a
unilateral deviation cannot offer benefit more than $\varepsilon$ while fixing
other’s profile. In the context of our example for the two level game, given
$(1,1)$’s action, we select the profile with the minimum $\varepsilon$ of the
simultaneous-move game defined on level 2 as an SPE back to $(1,1)$.

We now generalize this example to formally define the notion of
$\varepsilon$-SPE in an SHG. First let us define
$\Phi_{l,i}(\bm{x})\in\mathbb{R}^{d_{L}}$ that returns the equilibrated
profile at level $L$. In this profile, leaves that are not descendants of
$(l,i)$ are fixed in $\bm{x}$, while $\mathit{LEAF}(l,i)$ moved to a profile
that corresponding to an SPE of $\mathcal{G}_{l,i}(\bm{x})$ with an minimum
$\varepsilon$. Then we define the $\varepsilon_{l,i}(\bm{x})$ as
$\max_{\bm{x}^{\prime}_{l,i}\in\mathcal{X}_{l,i}}u_{l,i}(\bm{x}_{l,i}^{\prime},\bm{x}_{\textsc{Pa}(l,i)},\Phi(\bm{x}_{l,i}^{\prime}))-u_{l,i}(\bm{x}_{l,i},\bm{x}_{\textsc{Pa}(l,i)},\Phi(\bm{x}_{l,i}))$
and define $\varepsilon(\bm{x})=\max_{l,i}\varepsilon(\bm{x})$ as the
$\varepsilon$ of profile $\bm{x}$ in an SHG $\mathcal{G}$.

An SHG instance $\mathcal{G}$

Best response iterations $T$

procedure Search($\bm{x},(l,i)$)

for $\bm{x}_{l,i}^{\prime}\in\mathcal{X}_{l,i}$ do

$u_{l,i}^{re-
eq}(\bm{x}_{l,i}^{\prime}),\varepsilon_{\textsc{Des}(l,i)}(\bm{x}_{l,i}^{\prime}),\bm{x}^{re-
eq(l,i)}(\bm{x}_{l,i}^{\prime})\leftarrow\textsc{Re-
Eq}((l,i),\bm{x}:\bm{x}_{l,i}\rightarrow\bm{x}_{l,i}^{\prime})$

Find $\bm{x}_{l,i}^{*}\leftarrow\arg\max_{\bm{x}_{l,i}^{\prime}}u_{l,i}^{re-
eq}(\bm{x}_{l,i}^{\prime})$

return $\bm{x}_{l,i}^{*},u_{l,i}^{re-eq}(\bm{x}_{l,i}^{*}),\bm{x}^{re-
eq(l,i)}(\bm{x}_{l,i}^{*})$

end procedure

procedure Re-Eq($\mathcal{G},\bm{x},(l,i)$)

$\varepsilon\leftarrow 0$

if $l<L$ then

$\bm{x},\varepsilon\leftarrow\textsc{SHG\\_Solve}((l,i),\bm{x})$ return
$u_{l,i}(\bm{x}),\varepsilon,\bm{x}$

end procedure

procedure SHG_Solve($(l,i),\bm{x}$)

$\bm{x}^{0}\leftarrow\bm{x}$

$\forall(l+1,j)\in\textsc{Chd}(l,i)$, replace its action in $\bm{x}^{0}$ with
other random initialization

for $t=1,2,\dotsc,T$ do $\triangleright$ level-wise best response

for $(l+1,j)\in\textsc{Chd}(l,i)$ do

$u_{l+1,j}^{re-eq},\varepsilon_{\textsc{Des}(l+1,j)},\bm{x}^{re-
eq(l+1,j)}\leftarrow\textsc{Re-Eq}(\mathcal{G},\bm{x}^{t-1},(l+1,j))$

$\bm{x}_{l+1,j}^{\prime},u_{l+1,j}^{re-eq-max},\bm{x}^{re-eq-
max(l+1,j)}\leftarrow\textsc{Search}(\bm{x}^{t-1},(l+1,j))$

$\varepsilon_{l+1,j}^{t-1}=\max\\{\varepsilon_{\textsc{Des}(l+1,j)},u_{l+1,j}^{re-
eq-max}-u_{l+1,j}^{re-eq}\\}$

$\bm{x}^{t}\leftarrow\bm{x}^{t-1}$

for $\forall(l+1,j)\in\textsc{Chd}(l,i)$ do

Replace dimensions of $\bm{x}^{t}$ belonging to $\textsc{Des}(l+1,j)$ with the
ones in $\bm{x}^{re-eq-max(l+1,j)}$

$\varepsilon_{l+1}^{t-1}\leftarrow\max_{i}\varepsilon_{l+1,j}^{t-1}$

$t^{*}\leftarrow\arg\min_{t}\epsilon_{l+1}^{t}$

return $\bm{x}^{t^{*}},\varepsilon^{t^{*}}_{l+1}$

end procedure

procedure Compute_$\varepsilon$($\bm{x}$)

for $(l,i)$ do

$\bm{x}^{*}_{l,i},u_{l,i}^{*},\bm{x}^{re-
eq(l,i)}\leftarrow\textsc{Search}(\bm{x},(l,i))$

$\varepsilon_{l,i}\leftarrow u_{l,i}^{*}-u_{l,i}(\bm{x})$ return
$\max_{l,i}\varepsilon_{l,i}$

end procedure

Algorithm 2 Procedures for computing $\varepsilon$ of profile $\bm{x}$ and an
algorithm for computing an approximate SPE

We next provide an approximate algorithm (see Algorithm 2) to compute
$\epsilon$ in an SHG. Pay attention here that
$\bm{x}:\bm{x}_{i}\rightarrow\bm{x}_{i}^{\prime}$ means replacing $\bm{x}_{i}$
of $\bm{x}$ with $\bm{x}_{i}^{\prime}$. The functions in Algorithm 2 operate
as follows. The procedure Search takes a given joint profile $\bm{x}$, player
index $(l,i)$ and returns a best response profile for a given player $(l,i)$.
Our approach is that for each action in
$\bm{x}_{l,i}^{\prime}\in\mathcal{X}_{l,i}$, we re-equilibriate subgames
$\mathcal{G}_{l,i}(\bm{x}:\bm{x}_{l,i}\rightarrow\bm{x}_{l,i}^{\prime})$, and
then compute the corresponding payoff. In our actual implementation, we
discretize $\mathcal{X}_{l,i}$ in a bucket of grid points, and search within
such bucket. The procedure Re-Eq returns the re-equilibrated profile of
$\mathcal{G}_{l,i}(\bm{x})$ given $(l,i)$ and $\bm{x}$. The procedure
SHG_Solve solve a simultaneous-move game at level $l+1$. It applies an
iterative best response approach for $\textsc{Chd}(l,i)$ to generate diverse
profiles, and select the one with minimum possible $\varepsilon$. In each
iteration, for each $(l+1,j)\in\textsc{Chd}(l,i)$, it computes its best
response action against the previous joint profile $\bm{x}^{t-1}$. And then in
the next iteration it replace those descendant-profiles of $(l+1,j)$ in
$\bm{x}^{t-1}$ by computed re-equilibrated profiles when $(l+1,j)$ selected
its best response. Then it will return the joint profile with the minimum
$\varepsilon$ found so-far.

To solve the whole game, we just call $\textsc{SHG\\_Solve}((0,0),\bm{x})$,
where $\bm{x}$ is some other action profile. To compute the $\varepsilon$ of a
given profile, we just call Compute_ $\varepsilon$, where it just compute the
maximum unilateral deviation for every player using Search to compute the best
response action and payoff.

#### Structured Game Model Inspired by COVID-19 Policy-Making

We will now describe in detail the particular subclass of SHGs that we studied
in our experiments reported in Section 4.2. This class is based on the SHG
proposed in Jia et al. [18] which we describe in detail here. The exposition
is in terms of a cost function
$\mathtt{C}_{i}(\bm{x}_{i},\bm{x}_{\textsc{Pa}(i)},\bm{x}_{L})$ for each
player, which is more natural in this context, rather than the payoff function
$u_{l,i}$ introduced in Section 2, with the understanding that
$u_{i}\equiv-\mathtt{C}_{i}$.

There are $L=3$ levels in the hierarchy such that player $(1,1)$ represents
the federal government (or, simply, government), the players $(2,i)$,
$i\in\\{1,2,\dots,n_{2}\\}$ are state governments (or, simply, states), and
the players $(3,i)$, $i\in\\{1,2,\dots,n_{3}\\}$ are county governments (or,
simply, counties) partitioned into groups such that each group shares a single
state as a parent.

Each player $(i)$ takes a bounded, scalar action $x_{i}\in[0,1]$ which is a
social-distancing factor that (multiplicatively) reduces the proportion of
post-intervention contacts among individuals — a lower number implies a
stronger policy intervention, hence a lower number of infections but a higher
cost of implementation (see below). The actions taken by counties represent
policies that get actually implemented (hence directly impact the realized
cost of every player — one of the defining attributes of SHGs) while those
taken by the government and states are recommendations. Similar to Jia et al.
[18], we also study a restricted variant where each county is non-strategic
and constrained to comply with the action (recommendation) of its parent-
state, effectively reducing the model to a $2$-level hierarchy. We call this
special case a two-level game (Figure 4(a)) and the more general model a
three-level game (Figure 4(b)).

The cost function of each player $(i)$ has, in general, three components: a
_policy impact cost_ $\mathtt{C}^{\mathrm{inc}}_{i}(\bm{x}_{L})$ which we will
elaborate on below; a _policy implementation cost_
$\mathtt{C}^{\mathrm{dec}}_{i}(\bm{x}_{L})$, e.g. economic and psychological
costs of a lockdown; and, for each player in levels $l>1$, a _non-compliance
cost_ $\mathtt{C}^{\mathrm{NC}}_{i}(\bm{x}_{i},\bm{x}_{\textsc{Pa}(i)})$, a
penalty incurred by a policy-maker for deviating from the recommendation of
its parent in the hierarchy (e.g., a fine, litigation costs, or reputation
harm).

Let $N_{(3,i)}>0$ denote the fixed population under the jurisdiction of county
$(3,i)$ for every $i\in\\{1,2,\dots,n_{3}\\}$. By construction, the population
under each state is given
$N_{(2,i)}=\sum_{(3,j)\in\textsc{Chd}(2,i)}N_{(3,j)}$ and the population under
the government is $N_{(1,1)}=\sum_{i=1}^{n_{2}}N_{(2,i)}$. We next define the
expressions for each component of the cost function.

#### Policy Impact Cost:

This cost component is a quadratic closed-form approximation to the agent-
based model introduced by Wilder et al. [41]. This is inspired by the
infection cost computation approach in Jia et al. [18] but they used a
different closed-form approximation. For each each county $(3,i)$, let
$N^{init}_{(3,i)}$ denote the number of infected individuals within the
population of the county prior to policy intervention; thus, the number of
post-intervention susceptible individuals is
$(N_{(3,i)}-N^{init}_{(3,i)})x_{3,i}$. Another parameter in the game is the
the transport matrix
$R=\\{r_{aa^{\prime}}\\}_{a,a^{\prime}\in{(3,1),(3,2),\dots,(3,n_{3})}}$,
where $r_{aa^{\prime}}\geq 0$ is the proportion of the population of county
$a^{\prime}$ that is active in county $a$ in the absence of an intervention.
Thus, in the number of post-intervention infected individuals of county
$a^{\prime}$ that is active in county $a$ is
$r_{aa^{\prime}}N^{init}_{a^{\prime}}x_{a^{\prime}}$. The last parameter in
the model is $M$, the average number of contacts with active individuals that
a susceptible individual makes, and finally $\mu$ is the probability that a
susceptible individual gets infected upon contact with an active infected
individual is $\mu\in(0,1)$.

Putting these together, the policy impact cost is defined by the fraction of
post-intervention infected individuals in county $a=(3,i)$,
$i\in\\{1,2,\dots,n_{3}\\}$:

$\mathtt{C}^{\mathrm{inc}}_{a}(\bm{x}_{L})=\mu
Mx_{a}\frac{N_{a}-N_{a}^{init}}{N_{a}^{2}}(\sum_{a^{\prime}}r_{aa^{\prime}}N_{a^{\prime}}^{init}x_{a^{\prime}}).$

For a higher-level player $(i)$,

$\mathtt{C}^{\mathrm{inc}}_{i}(\bm{x}_{L})=\frac{1}{N_{(i)}}\sum_{a\in\textsc{Chd}(i)}N_{a}\mathtt{C}^{\mathrm{inc}}_{a}(\bm{x}_{L}).$

#### Policy Implementation Cost:

For each county $(3,i)$, the policy implementation cost is given by

$\mathtt{C}^{\mathrm{dec}}_{3,i}(\bm{x}_{L})=1-x_{3,i}.$

For a higher-level player $(i)$,

$\mathtt{C}^{\mathrm{dec}}_{i}(\bm{x}_{L})=\frac{1}{N_{(i)}}\sum_{a\in\textsc{Chd}(i)}N_{a}\mathtt{C}^{\mathrm{dec}}_{a}(\bm{x}_{L}).$

#### Non-Compliance Cost:

The non-compliance cost of player $(i)$ for $l\in\\{2,3\\}$ is given by
Euclidean distance between its action and that of its parent:

$\displaystyle\mathtt{C}^{\mathrm{NC}}_{i}(x_{i},x_{\textsc{Pa}(i)})=(x_{i}-x_{\textsc{Pa}(i)})^{2}.$

Finally, each player $(i)$ for $l>1$ has an idiosyncratic set of weights
$\kappa_{i}\geq 0$ and $\eta_{i}\geq 0$ that trade its three cost components
off against each other via a convex combination, and account for differences
in ideology; the overall cost of such a player is given by

$\displaystyle\mathtt{C}_{i}(x_{i},x_{\textsc{Pa}(i)},\bm{x}_{L})=\kappa_{i}\mathtt{C}^{\mathrm{inc}}_{i}(\bm{x}_{L})+\eta_{i}\mathtt{C}^{\mathrm{dec}}_{i}(\bm{x}_{L})+(1-\kappa_{i}-\eta_{i})\mathtt{C}^{\mathrm{NC}}_{i}(x_{i},x_{\textsc{Pa}(i)}).$
(7)

The player $(1,1)$ obviously has no non-compliance issues, hence it has only
one weight $\kappa_{1,1}>0$, its overall cost being

$\displaystyle\mathtt{C}_{1,1}(x_{i},\bm{x}_{L})=\kappa_{1,1}\mathtt{C}^{\mathrm{inc}}_{1,1}(\bm{x}_{L})+(1-\kappa_{1,1})\mathtt{C}^{\mathrm{dec}}_{1,1}(\bm{x}_{L}).$

In our experiments, we set $r_{aa^{\prime}}=1/n_{3}$ for every pair of
counties $(a,a^{\prime})$, $M=20$ and $\mu=0.3$.

#### Experimental Setup and Further Results

We formally define the algorithm BRD as the procedure
$\textsc{SHG\\_Solve}((0,0),\bm{x})$ for some randomly initilized $\bm{x}$. We
start by defining the concept of one full algorithm iteration for each of DBI
and BRD. For DBI(1, 20), DBI(1, 50), DBI(1, 2, 4), DBI(1, 2, 10), one algrithm
iteration consists of 50 steps of gradient ascent, with a learning rate of
0.01. For BRD(1, 20) and BRD(1, 50), one algorithm iteration consists of 20
iterations of level-wise best response during the recursive procedure in
SHG_Solve; for BRD(1, 2, 4) and BRD(1, 2, 10), it corresponds to 500 and 200
iterations of level-wise best response, respectively. For DBI we adopt a
projector operator that project the resulted action into the nearest point in
$[0,1]$.

We discretize each action space uniformly into 101 grid points for two-level
games, and 11 grid points for three-level games. We let $T=100$ for BRD(1, 20)
and BRD(1, 50) and $T=20$ for BRD(1, 2, 4) and BRD(1, 2, 10). In the three-
level experiments, the $\kappa$ is set to be $0.5$ for counties and the states
and $0.8$ for the government. The $\eta$ is set to be $0.2$ for the states,
$0.3$ for the counties in (1,2,4) setting , and $0.2$ for (1,2,10) experiment.
In the two-level experiments, the $\kappa$ is set to be $0.2$ for the
government, $0.5$ for counties and the states. The $\eta$ is set to be $0.2$
for counties and states.

(a) Two-level games. (b) Three-level games.

Figure 8: The run-time (in secs) on COVID-19 Game Models. T = 20 for BRD
experiments.

To further study the scalability of the BRD and DBI algorithms, we compare the
run-time that (1) the DBI algorithm converges and (2) the BRD algorithm
terminates. The convergence of DBI is defined as the convergence of the action
profile. When DBI converges, the action profile remains unchanged. In the
projected gradient descent method, if all the gradients go to zero, DBI
converges. Besides, it is also likely to hit the boundary of the constraints.
The DBI is possible to converge with a non-zero gradient norm. The BRD
algorithm terminates either T achieves, or $\epsilon$ goes to zero for all
players. Under these conditions, the two algorithms find their optimal
solutions. Figure 8 demonstrates the run-time results. We conduct each
experiment four times with different random seeds. In two-level problems, the
DBI algorithm is more than two times faster than BRD algorithms. Although the
action spaces are discretized to 11 grid points for three-level games rather
than 101 grids, DBI algorithms still perform better. In practice, discretizing
the action spaces in 101 grids for three-level games is computational
intensively. The performance and run-time of the BRD algorithm are more
dependent on randomization and the initial points for the best responses at
each level. When we face many players or multiple levels, the DBI algorithm is
a natural choice. The DBI algorithm is significantly more efficient and more
stable than the BRD algorithm.

###  C.4 Details of Hierarchical Public Goods Game Experiments

We extend the described public good game on Zachary’s Karate club network to a
(1-2-34) hierarchical game by introducing the non-compliance costs similar to
the COVID-19 game model. The overall utility of the club individuals should be
the combination of public good utility $u_{i}$ and
$\text{C}_{i}^{\text{NC}}(x_{i},x_{\text{PA(i)}})=(x_{i}-x_{\text{PA(i)}})^{2}$.

$\displaystyle
U_{3,i}(x_{i},x_{\text{PA(i)}},\mathbf{x}_{L})=(1-\kappa_{3,i})u_{i}(x_{i},\mathbf{x}_{L})-\kappa_{3,i}\text{C}_{i}^{\text{NC}}(x_{i},x_{\text{PA(i)}}).$
(8)

In the second level, the administrator and instructor’s utility is composed of
the total public good utility of their group members, and their own non-
compliance costs.

$\displaystyle
U_{2,i}(x_{i},x_{\text{PA(i)}},\mathbf{x}_{L})=(1-\kappa_{2,i})\sum_{j\in\text{CHD}(i)}u_{j}(x_{j},\mathbf{x}_{L})-\kappa_{2,i}\text{C}_{i}^{\text{NC}}(x_{i},x_{\text{PA(i)}}).$
(9)

The root player’s utility is considered as the social welfare of the whole
club.

$\displaystyle
U_{1,1}(x_{i},x_{\text{PA(i)}},\mathbf{x}_{L})=\sum_{j}u_{j}(x_{j},\mathbf{x}_{L})$
(10)

In our experiments, we set $\kappa_{2,i}=\kappa_{3,j}=0.5,\ \forall i,j$.
Other parameters of the public good utility are set to be
$a_{i}=0,b_{i}=1,c_{i}=6,\ \forall i$. We use learning rate $0.1$ in DBI. We
project the results to $[0,1]$. To compare with the DBI, we use BRD with
discretized factors $0.5,0.2,0.1,0.05$, and best response rounds $2,3$.

###  C.5 Details of Hierarchical Security Game Experiments

Let $a=(a_{1},\ldots,a_{n})$ be the attacker strategy, the probability
distribution over which defenders are attacked. In our experiments, $a$ is set
to follow a logit distribution $\mathrm{softmax}(\lambda(1-\mathbf{x}_{L}))$,
with defenders having lower security investment more likely to be attacked.
Let $c_{i}(x_{i})=c_{i}x_{i}$ be the cost of security investment for defender
$i$, which is assumed to be linear. The utility of defender $i$ is then as
follows:

$u_{i}(x_{i},\mathbf{x}_{L})=a_{i}\frac{x_{i}}{1+x_{i}}+\sum_{j\neq
i}a_{j}(1-q_{ji}\frac{1}{1+x_{j}}\frac{1}{1+x_{i}})-c_{i}(x_{i})$

.

We then build the (1,3,6) structured hierarchical game by introducing the non-
compliance costs
$\text{C}_{i}^{\text{NC}}(x_{i},x_{\text{PA(i)}})=(x_{i}-x_{\text{PA(i)}})^{2}$.
The utility of each agent in different levels are the same as that of public
good game (Eq. 8, 9, 10).

In the experiments, we consider the agents who almost obey the parent’s
suggestion ($\kappa=0.1$) and the somewhat selfish agents ($\kappa=0.5$). The
influence probability $q_{ji},\ \forall i,j$ are set to be equal among the
agents, and other parameters are set to $c_{i}=0.2$, $\lambda=5$. We use a
learning rate of $0.1$ in DBI. To compare with the DBI, we use BRD with
discretized factors $0.2,0.1,0.05$ and best response rounds $2,3$. Since there
is no upper bound of the strategy profile, we set $x_{i}\leq 1.0$ according to
the performance. Then we search the strategy space $x\in[0,1]$ in the BRD
execution.