Data Structures for Deviation Payoffs

TL;DR

This paper introduces optimized data structures for symmetric normal-form games that significantly accelerate the computation of symmetric mixed-strategy Nash equilibria, enabling practical analysis of large multi-player games.

Contribution

The paper presents novel data structures tailored for symmetric games, improving computational efficiency and extending applicability to role-symmetric and action-graph games.

Findings

Dramatic speedup in equilibrium computation

Efficient representation for large multi-player games

Extension to role-symmetric and action-graph games

Abstract

We present new data structures for representing symmetric normal-form games. These data structures are optimized for efficiently computing the expected utility of each unilateral pure-strategy deviation from a symmetric mixed-strategy profile. The cumulative effect of numerous incremental innovations is a dramatic speedup in the computation of symmetric mixed-strategy Nash equilibria, making it practical to represent and solve games with dozens to hundreds of players. These data structures naturally extend to role-symmetric and action-graph games with similar benefits.

Tables

Table 1. Table 1. For a given number of actions A 𝐴 A , the largest number of players P 𝑃 P before repetitions overflows a 64-bit integer.

$A$	2	3	4	5	6	7	8	9	10
$\max P$	67	44	36	32	29	27	26	25	25
$A$	11	12	13	14	15	16	17	18	$\ge \text{19}$
$\max P$	24	23	23	23	23	22	22	22	21

Equations

$$\text{Reps}\left(\text{@vec}c\right)=\left(\genfrac{}{}{0ex}{}{P-1}{\text{@vec}{c}_1,\text{@vec}{c}_2,\dots ,\text{@vec}{c}_A}\right)=\frac{(P-1)!}{\text{@vec}{c}_1!\text{@vec}{c}_2!\dots \text{@vec}{c}_A!}$$

$$\text{@vec}{u}_a(\text{@vec}\sigma )=\sum _{\text{@vec}s:\text{@vec}{s}_a>0}\underset{\text{@vec}\sigma }{\mathrm{Pr}}\left(\text{@vec}s\mid a\right)v_a\left(\text{@vec}s\mid a\right)$$

$$\text{@vec}{u}_a(\text{@vec}\sigma )$$

$$=\sum _{\text{@vec}c}{v}_a\left(\text{@vec}c\right)\text{Reps}\left(\text{@vec}c\right)\prod _{a^{\prime }}{\left(\text{@vec}{\sigma }_{a^{\prime }}\right)}^{\text{@vec}{c}_{a^{\prime }}}$$

$$\frac{\partial \text{@vec}{u}_a(\text{@vec}\sigma )}{\partial \text{@vec}{\sigma }_s}=\sum _{\text{@vec}c}{v}_a(\text{@vec}c)\text{Reps}\left(\text{@vec}c\right)\left(\text{@vec}{c}_s\right){\left(\text{@vec}{\sigma }_s\right)}^{\text{@vec}{c}_s-1}\prod _{a^{\prime }\ne s}{\left(\text{@vec}{\sigma }_{a^{\prime }}\right)}^{\text{@vec}{c}_{a^{\prime }}}$$

$$\text{@vec}{w}_a$$

$$\text{@vec}{\sigma }_a^{t+1}$$

$$g(\text{@vec}\sigma )=\sum _a\max \left(0,\text{@vec}{u}_a(\text{@vec}\sigma )-\text{@vec}u(\text{@vec}\sigma )\cdot \text{@vec}\sigma \right)$$

$${\nabla }_s(g)=\sum _a\left(\frac{\partial \text{@vec}{u}_a(\text{@vec}\sigma )}{\partial s}-\text{@vec}{u}_s(\text{@vec}\sigma )-\sum _{a^{\prime }}\text{@vec}{\sigma }_{a^{\prime }}\frac{\partial \text{@vec}{u}_{a^{\prime }}(\text{@vec}\sigma )}{\partial s}\right)1_{g_s}$$

$$\underset{\text{@vec}\sigma }{\mathrm{Pr}}\left(\text{@vec}s\mid a\right)=\text{Reps}\left(\text{@vec}s\mid a\right){\left(\text{@vec}{\sigma }_a\right)}^{\text{@vec}{s}_a-1}\prod _{a^{\prime }\ne a}{\left(\text{@vec}{\sigma }_{a^{\prime }}\right)}^{\text{@vec}{s}_{a^{\prime }}}$$

$$m^i\leftarrow \text{@vec}{s}_a^i\ne 0$$

$$\text{@vec}{c}^i\leftarrow \text{@vec}{s}^i-\text{@vec}{e}_a$$

$$p^i\leftarrow \text{Reps}\left(\text{@vec}{c}^i\right)\prod _{a^{\prime }}{\left(\text{@vec}{\sigma }_{a^{\prime }}\right)}^{\text{@vec}{c}_{a^{\prime }}^i}$$

$$\text{@vec}{u}_a(\text{@vec}\sigma )\leftarrow \sum _i{m}^i{p}^i{v}_a\left(\text{@vec}{c}^i\right)$$

$$p^j\leftarrow \text{Reps}\left(\text{@vec}{c}^j\right)\prod _{a^{\prime }}{\left(\text{@vec}{\sigma }_{a^{\prime }}\right)}^{\text{@vec}{c}_{a^{\prime }}^j}$$

$$\text{@vec}u(\text{@vec}\sigma )\leftarrow \sum _j{p}^j\text{@vec}{v}^j$$

$$\@vec{u}_{a}(\@vec{\sigma})\leftarrow\sum_{j}\Bigl{(}v_{a}\left(\@vec{c}^{j}\right)\textrm{Reps}\left(\@vec{c}^{j}\right)\Bigr{)}\prod_{a^{\prime}}\left(\@vec{\sigma}_{a^{\prime}}\right)^{\@vec{c}^{j}_{a^{\prime}}}$$

$$\text{Reps}\left(\text{@vec}{c}^j\right)v_a\left(\text{@vec}{c}^j\right)$$

$${\lambda }_a^j\leftarrow \log \left(\text{Reps}\left(\text{@vec}{c}^j\right)\right)+\log \left(v_a(\text{@vec}{c}^j)\right)$$

$${\gamma }^j\leftarrow \exp \left({\lambda }^j+\sum _{a^{\prime }}\text{@vec}{c}_{a^{\prime }}^j\log \text{@vec}{\sigma }_{a^{\prime }}\right)$$

$$\text{@vec}{u}_a(\text{@vec}\sigma )\leftarrow \sum _j{\gamma }^j$$

$$\frac{\partial \text{@vec}{u}_a(\text{@vec}\sigma )}{\partial \text{@vec}{\sigma }_s}\leftarrow \sum _j\frac{\text{@vec}{c}_s^j}{\text{@vec}{\sigma }_s}{\gamma }^j$$

$$\left[\begin{array}{c}3\\ 0\\ 0\end{array}\right]$$

$$\left[\begin{array}{c}1\\ 0\\ 2\end{array}\right]$$

$$\text{@vec}{u}_1(\text{@vec}\sigma )=$$

$$=$$

$$\mathrm{Pr}(⟨1,2,0⟩)\cdot v_1(⟨1,2,0⟩)+\mathrm{Pr}(⟨1,1,1⟩)\cdot v_1(⟨1,1,1⟩)+\mathrm{Pr}(⟨1,0,2⟩)\cdot v_1(⟨1,0,2⟩)$$

$$=$$

$$(1)(0.5^2)(1)+(2)(0.5\cdot 0.4)(1)+(1)(0.4^2)(1)$$

$$=$$

$$\sum\Biggl{(}$$

$$\sum_{rows}\Biggl{(}$$

$$\Biggr{)}\;\;=\;\;\begin{bmatrix}0.63\\ -0.5\\ -0.36\end{bmatrix}$$

$$\sum _{rows}\left(\cdot \right)=\left[\begin{array}{c}0.63\\ -0.5\\ -0.36\end{array}\right]$$

Full text

1. Introduction

Symmetric games, where all players face identical incentives, play a central role in game-theoretic analysis. Many of the classic examples used to teach game-theoretic concepts Leyton-Brown and Shoham (2008); Shoham and Leyton-Brown (2008) such as prisoner’s dilemma and rock-paper-scissors are symmetric, and the original paper proposing the concept of a Nash equilibrium Nash (1951) addressed symmetry and proved that symmetric games must have symmetric equilibria. Symmetric games are quite common in the recent multi-agent systems literature Wang et al. (2021); Boix-Adserà et al. (2020). The role of symmetry becomes especially prominent in games with a large number of players, because it is often only with the help of player symmetry that incentives can even be tractably described, and multiple distinct sub-fields of computational game theory Wellman (2020); Caines (2021) rely heavily on foundational assumptions about player symmetry to scale up game-theoretic analysis.

Despite this importance, data structures for representing and solving large symmetric games have received a paucity of attention in the research literature. Most libraries for representing and solving games, particularly Gambit McKelvey et al. (2006), but also GameTracer Blum et al. (2006) and QuantEcon quantecon.org (2022), include few tools for efficiently solving large symmetric games. One notable exception is the work of Jiang et al. (2011) on action-graph games, which incorporate player symmetry, and which have been partially implemented in Gambit. However, action-graph games focus primarily on theoretical compactness in the action space, and only consider expected utility calculations enough to ensure they take polynomial time. The data structures that have been implemented are not specifically optimized for equilibrium computation under player symmetry, meaning that solving large symmetric games remains inefficient.

In the sub-fields that scale up symmetric game analysis, approximations play a central role. Simulation-based games Wellman (2020) frequently employ player reduction methods Wiedenbeck and Wellman (2012) that use a reduced game with a very small number of players to replace a game with dozens or hundreds of players in the hope that analysis of the reduced-game will approximately hold in the large game. Mean-field games Caines (2021) give up on representing discrete players and replace them with an average effect that aggregates over a continuum. Other approaches aim to identify underlying representational compactness Ficici et al. (2008); Li and Wellman (2020) or to avoid explicitly constructing game models in favor of sampling-based techniques Gemp et al. (2022); Sokota et al. (2019). These approximations are sometimes unavoidable, when the representation is too large or the payoff computations are too slow for a game to even be enumerated, but analysts often end up resorting to approximation when faced with any non-trivial number of players. The result is a substantial gap between the very small games that can be represented and solved exactly and the very large games where approximation is necessary. We bridge this gap by designing efficient data structures that make it practical to exactly solve much larger instances.

In this work, we present a detailed exploration of data structure improvements for symmetric normal-form games, with specific focus on the task of computing symmetric mixed-strategy Nash equilibria. We argue, following Cheng et al. (2004), that this is by far the most compelling solution concept for symmetric games, because symmetric equilibria are guaranteed to exist, and it provides greater intuitive explanatory power for symmetries of the game to be reflected in its equilibria. While computing symmetric Nash equilibria can be hard in the worst case Daskalakis et al. (2009); Conitzer and Sandholm (2002), incomplete local-search algorithms such as replicator dynamics and gradient descent are often highly successful in practice.

To facilitate algorithms for identifying symmetric equilibria, we aim for data structures that optimize the calculation of deviation payoffs (and their derivatives). For a symmetric mixed strategy employed by all players, the deviation payoff vector gives, for each action, the expected utility a single agent would receive if they deviated unilaterally to that action. For a large fraction of the algorithms used to compute Nash equilibria, the ability to evaluate deviation payoffs for symmetric mixed strategies is both necessary and sufficient; for most other algorithms, deviation payoffs plus their derivatives suffice. We show that focusing the design of symmetric game data structures toward deviation payoffs leads to a number of optimizations that jointly yield a dramatic improvement in the practical efficiency of solving large symmetric games.

1.1. Contributions

We describe seven distinct upgrades to classic data structures for representing symmetric normal-form games. Two of these (sections 3.3 and 3.4) constitute asymptotic improvements: changing from $P$-player profiles to $(P-1)$-opponent configurations reduces the number of stored payoff entries (with $A$ actions) from $A\binom{P+A-1}{P}$ to $A\binom{P+A-2}{P-1}$, and pre-computing probability weights to avoid repeated multinomial calculations accelerates equilibrium searches by a factor of $A$. Four of the upgrades focus on vectorization, resulting in better constants and enabling SIMD acceleration. And we find—as has been well-established in the settings of neural networks and scientific computation—that such improvements can qualitatively change the scope of the computational tools. The remaining upgrade (section 3.6) follows because the overall expansion-of-scope enables us to analyze games so large that the probability calculations can overflow 64-bit integers, necessitating a switch to a log-space representation of payoffs and probabilities.

Our main result is a roughly ten-thousand-fold speedup in the running time of equilibrium computation algorithms for symmetric games with many players. This makes it possible to run practical-but-incomplete search methods like replicator dynamics on any symmetric game that can fit in memory, and can also facilitate other slower equilibrium-search techniques. Our results effectively close the gap between small and large symmetric games, relegating approximation techniques to only those games too large to be represented explicitly.

Two open-source libraries implement most of the data structures and algorithms we discuss. The gameanalysis.jl111https://github.com/Davidson-Game-Theory-Research/gameanalysis.jl repository provides simple Julia implementations of the best CPU and GPU versions of our symmetric game data structure, and also includes all of the code for experiments in this paper. The gameanalysis.py222https://github.com/egtaonline/gameanalysis module implements in Python the data structure variants for role-symmetric games that we discuss in Section 5.2 and provides a richer command-line interface with numerous tools for solving role-symmetric games.

2. Background

2.1. Terminology

In a symmetric game, all players have the same action set and identical incentives. We therefore refer to the number of players $P$, but rarely distinguish individual players. We call the number of actions $A$, and will often index the actions by $a\in\{1,\ldots,A\}$.

A profile specifies an action for each player, and in a symmetric game, we can represent a profile by an integer-vector $\@vec{s}$ specifying a non-negative number of players selecting each action. We denote the entries in a vector with a subscript, so $\@vec{s}_{a}$ is the player-count for action $a$. We use a superscript $i$ or $j$ to index vectors in a collection: for example $\@vec{s}^{i}$; to avoid ambiguity, exponents are generally applied to a quantity in parentheses. We will also distinguish a profile from an opponent configuration $\@vec{c}$, which differs only in that $\sum_{a}\@vec{s}_{a}=P$, while $\sum_{a}\@vec{c}_{a}=P-1$. We refer to the configuration resulting from removing action $a$ from profile $\@vec{s}$ as $(\@vec{s}\mid a)$, since it will often appear in probability calculations and other contexts where it is given that one player selects action $a$. In terms of the integer-vector representation, $(\@vec{s}\mid a)$ subtracts 1 from dimension $a$ of $\@vec{s}$. A symmetric game’s payoffs can be expressed as the value achieved by a player choosing action $a$ when opponents play configuration $\@vec{c}$, which we denote by $v_{a}\left(\@vec{c}\right)$ normally, or by $v_{a}\left(\@vec{s}\mid a\right)$, when working in terms of profile $\@vec{s}$.

A mixed strategy specifies a probability distribution over actions, and a mixed strategy used by all players is called a symmetric mixed-strategy profile. We denote this with the variable $\@vec{\sigma}$, and will often abbreviate the term by referring to a mixture. When computing probabilities, we will frequently refer to the number of asymmetric arrangements corresponding to a symmetric configuration, which we call repetitions. We denote this quantity by $\textrm{Reps}\left(\@vec{c}\right)$ or $\textrm{Reps}\left(\@vec{s}\mid a\right)$, and calculate it with the following multinomial:

[TABLE]

2.2. Deviation Payoffs

When analyzing a symmetric game, we are most often interested in computing symmetric mixed-strategy Nash equilibria. For many algorithms that compute such equilibria a necessary and sufficient condition is the ability to compute deviation payoffs, and for most other algorithms, deviation payoffs plus deviation derivatives suffice. We begin by formally defining these terms, and then describe their application to our preferred equilibrium computation algorithms in section 2.3.

Given a symmetric mixed-strategy profile $\@vec{\sigma}$, we define the deviation payoff $\@vec{u}_{a}(\@vec{\sigma})$ for action $a$ as the expected utility one player would receive if they played $a$ while all opponents randomized according to $\@vec{\sigma}$. This expectation is often expressed as a sum over all profiles in which $a$ is played of the probability that profile occurs times the payoff to $a$ in that profile:

[TABLE]

but can be stated much more cleanly using configurations:

[TABLE]

The deviation payoff vector $\@vec{u}(\@vec{\sigma})$ collects these values for all actions $a\in\{1,\ldots,A\}$. We call the partial derivatives of deviation payoffs with respect to mixture probabilities deviation derivatives. Specifically, $\frac{\partial\@vec{u}_{a}(\@vec{\sigma})}{\partial\@vec{\sigma}_{s}}$ gives the change in the deviation payoff for action $a$ as the probability of action $s$ is varied. Again this can be expressed in terms of profiles, but is more straightforward in terms of configurations: 333Note that efficient computation of $\left(\@vec{\sigma}_{s}\right)^{\@vec{c}_{s}-1}$ can result in numerical errors for mixtures where $\@vec{\sigma}_{s}=0$. This sort of error can be avoided here and elsewhere with no real loss of precision by lower-bounding mixture probabilities at machine-epsilon.

[TABLE]

We can define most other standard game-theoretic quantities for symmetric games in terms of deviation payoffs. The expected utility experienced by all players when following a symmetric mixed strategy is given by the dot product $u(\@vec{\sigma})=\@vec{u}(\@vec{\sigma})\cdot\@vec{\sigma}$. The regret of a mixture is $\textrm{reg}(\@vec{\sigma})=\max_{a}\left(\@vec{u}_{a}(\@vec{\sigma})-u(\@vec{\sigma})\right)$. A symmetric Nash equilibrium is a mixture with $\textrm{reg}(\@vec{\sigma})=0$, while an approximate Nash equilibrium has $\textrm{reg}(\@vec{\sigma})\leq\varepsilon$, for suitable $\varepsilon$.

2.3. Equilibrium Computation

Computing deviation payoffs and/or their derivatives is the key step for a number of algorithms that identify Nash equilibria in symmetric games. We describe two of the most practical algorithms here: replicator dynamics, which depends on deviation payoffs, and gradient descent on sum-of-gains, which depends on deviation derivatives. Our data structures can also efficiently support a number of other Nash algorithms, including fictitious play Brown (1951), Scarf’s simplicial subdivision Scarf (1967), and the global Newton method of Govindan and Wilson (2003), as well as some algorithms for correlated equilibria Jiang and Leyton-Brown (2011). The Nash Algorithms Appendix444See Appendix C: Nash Algorithms. https://arxiv.org/abs/2302.13232 presents further details on how some of these algorithms depend on deviation payoffs.

Replicator dynamics

Taylor and Jonker (1978) is often presented as a rule governing evolutionary dynamics, but can also be viewed as an algorithm for computing symmetric Nash equilibria. Expressed in terms of deviation payoffs, replicator dynamics starts from some initial mixture $\@vec{\sigma}^{0}$ at $t=0$, and performs iterative updates of the form:

[TABLE]

This update assumes that all payoffs are non-negative; a positive affine transformation can be applied to any game to ensure this assumption holds (and to adjust the effective step-size).

Gradient descent

is a classic local-search algorithm for minimizing differentiable functions. We can easily define a function whose minima correspond to Nash equilibria based on the following sum of deviation gains:

[TABLE]

Then we can iteratively take steps in the direction of $-\nabla g(\@vec{\sigma})$. The elements of this gradient vector are given by:

[TABLE]

Where $1_{g_{s}}$ is an indicator variable for ${\@vec{u}_{s}(\@vec{\sigma})>\@vec{u}(\@vec{\sigma})\cdot\@vec{\sigma}}$. So deviation payoffs and deviation derivatives suffice to compute the gain gradient. When performing gradient descent, the mixture resulting from $\@vec{\sigma}-\nabla g(\@vec{\sigma})$ may not lie on the probability simplex, so it is necessary to project each step back onto the simplex, which we do using the method from Wang and Carreira-Perpinán (2013).

Neither replicator dynamics nor gradient descent is guaranteed to identify a Nash equilibrium. However, since these algorithms are many orders of magnitude faster than complete algorithms like simplicial subdivision or global Newton, it is practical to re-run them from many initial mixtures. We find that both algorithms tend to identify multiple equilibria in large symmetric games and that the sets of equilibria they find are often only partially overlapping. On the other hand, we find that fictitious play and other best-response-based updates are ineffective on large symmetric games. Therefore in practice, we recommend repeatedly running both replicator dynamics and gradient descent, filtering by regret to isolate $\varepsilon$-Nash mixtures, and merging the resulting equilibrium sets. An example trace from running both algorithms on a 100-player, 3-action Gaussian mixture game (represented using our data structures) appears in Figure 1.

3. Data Structure Improvements

The classic payoff-matrix representation of a normal-form game has a dimension for each player, and a size along each dimension equal to the number of actions available to that player. In a symmetric game, it suffices to store just one player’s payoffs, so if there are $P$ players and $A$ actions, a symmetric game can be represented by a symmetric tensor of size $A^{P}$.

Because symmetric tensors also arise in other settings, it is worth considering whether generic techniques for symmetric tensors would suffice for efficiently representing symmetric games. In particular, Schatz et al. (2014) propose a block-representation of symmetric tensors and a cache-conscious algorithm for multiplying with them. In principle, their representation has the same level of asymptotic compression as the data structures presented here. However, in Figure 2, we compare the memory footprint of the Schatz et al. symmetric-tensor representation (as implemented by the SymmetricTensors.jl library) against the largest of our symmetric-game data structures (pre-computed repetitions), showing that while the block-symmetric representation is better than storing the full payoff matrix, it is significantly larger than any of our data structure variants.

3.1. Payoff Dictionary

We therefore turn to purpose-built data structures for storing symmetric games. The classic approach is to store a mapping from profiles to payoffs. In Gambit McKelvey et al. (2006) the implementation of symmetric action-graph games stores this mapping in a trie, while others have used a database Cassell and Wellman (2012). Such a mapping is also implicit in many theoretical descriptions, for example when relating the size of a compact representation to the number of profiles Jiang et al. (2011, 2017). In our experiments, our baseline data structure stores this mapping using a hash table that maps a vector of integers representing the profile $\@vec{s}$ to a vector of floats that stores the payoff $v_{a}(\@vec{s})$ for each played action with $\@vec{s}_{a}>0$. With $P$ players and $A$ actions, this table will store $\binom{P+A-1}{P}$ profiles.

Calculating a deviation payoff $u_{a}(\@vec{\sigma})$ using this mapping requires iterating through all profiles to compute the sum from Equation (2), where the probabilities are given by this expression:

[TABLE]

$\textrm{Reps}\left(\@vec{s}\mid a\right)$, calculated according to the multinomial from Equation (1), gives the number of asymmetric orderings of the symmetric configuration $\left(\@vec{s}\mid a\right)$, while the remaining terms gives the probability of the opponents jointly playing one such asymmetric ordering. The Worked Examples Appendix shows the full mapping and a detailed walk-through of the deviation payoff calculation on an example 3-player, 3-strategy symmetric game for this version of the data structure as well as for all subsequent variants. We strongly recommend stepping through these examples to help resolve any confusion that arises from the data structure and algorithm descriptions that follow.555See Appendix A: Worked Examples. https://arxiv.org/abs/2302.13232

3.2. Array Vectorization

The first idea for improving this deviation payoff calculation is to store the profiles $\@vec{s}$ and the payoffs $v$ in a pair of two-dimensional arrays with parallel structure, denoted by a shared index $i$ over all profiles. In both arrays, each row corresponds to an action and each column corresponds to a profile, resulting in arrays of size $A\times\binom{P+A-1}{P}$. Extracting column $i$ of the profiles-array gives a profile vector $\@vec{s}^{i}$ with the count for each action. The corresponding column of the payoffs-array stores the payoff $v_{a}\left(\@vec{s}^{i}\mid a\right)$ for each action $a$ where that profile has a non-zero count.

Of note, the array representation here (and in all subsequent data structure variants) hinders direct look-up of payoffs for a pure-strategy profile using the mapping. However, this trade-off is clearly worthwhile for several reasons:

1. (1)

Calculating mixed-strategy deviation payoffs is a much tighter bottleneck than looking up pure-strategy payoffs. 2. (2)

Symmetric mixed-strategy equilibria in symmetric games are more often relevant than asymmetric pure-strategy equilibria (and are guaranteed to exist). 3. (3)

The profile-to-index mapping can be computed by a ranking algorithm for combinations-with-replacement Loehr (2017), or can be stored alongside the arrays for fast lookup with only linear additional memory.

Using this array representation, we can vectorize each of the steps in the calculation of the deviation payoff $\@vec{u}_{a}(\@vec{\sigma})$. We describe each of the following operations in terms of profile $i$, and broadcast those operations across array columns. First, we can compute a mask $m$ identifying the profiles where action $a$ is played:

[TABLE]

Then for each such profile, we can remove a player choosing $a$ to get the opponent-configuration $\@vec{c}^{i}=(\@vec{s}^{i}|a)$ by subtracting an indicator vector $\@vec{e}_{a}$ for action $a$:

[TABLE]

For each profile, the probability $p^{i}=\Pr\left(\@vec{c}^{i}\right)$ of the opponent configuration is calculated by a multinomial to compute repetitions along with a broadcast-and-reduce to apply exponents and take the product, resulting in an array of configuration probabilities:

[TABLE]

Finally the deviation payoff for action $a$ can be computed by extracting row $a$ of the payoffs array, multiplying element-wise by the configuration probabilities, masking out the invalid configurations, and summing over the remaining profiles:55footnotemark: 5

[TABLE]

Since these steps can be performed using only arithmetic primitives and indexed reduction, the computation is vectorizable in any numeric library. Further vectorization over actions is in principle possible, but we defer this to later sections. The degree of improvement from array vectorization will vary with the language used. More importantly, it sets the stage for our subsequent innovations.

3.3. Pre-computed Repetitions

The next improvement comes from recognizing that each time we compute a deviation-payoff vector $\@vec{u}(\@vec{\sigma})$, we make use of the $\textrm{Reps}\left(\@vec{s}\mid a\right)$ value for every profile-action pair. Since equilibrium search involves many, many deviation-payoff calculations, we can save significant work by pre-computing these repetitions and storing them in a third parallel array. Then the probability calculation in Equation (9) can extract row $a$ of this repetitions array and element-wise multiply by the array of products. Our gameanalysis.jl library11footnotemark: 1 and the Worked Examples Appendix55footnotemark: 5 each include a version of this data structure that computes $\@vec{u}(\@vec{\sigma})$ in a single pass instead of computing $\@vec{u}_{a}(\@vec{\sigma})$ separately for each action, but we defer the detailed explanation of this vectorization to the following variant where it is dramatically simplified.

3.4. Opponent Configurations

A closer inspection of these parallel array structures reveals some redundancies: first the payoffs array contains a number of meaningless entries, since whenever $\@vec{s}_{a}=0$ the corresponding $v_{a}$ is masked out of the calculations.55footnotemark: 5 Second, the repetitions array contains a number of redundant entries: whenever two profiles differ by a single player’s action, that is $(\@vec{s}\mid a)=(\@vec{s}\,^{\prime}\mid a^{\prime})$ we will end up storing identical entries for $\textrm{Reps}\left(\@vec{s}\mid a\right)$ and $\textrm{Reps}\left(\@vec{s}\,^{\prime}\mid a^{\prime}\right)$.

Both of these issues can be avoided if we switch to storing opponent-configurations, with a shared index $j$ over all configurations. This gives us two parallel arrays of size $A\times\binom{P+A-2}{P-1}$ for configurations and payoffs, and a $1\times\binom{P+A-2}{P-1}$ repeats-array. The configurations-array now stores each configuration $\@vec{c}^{j}$ over $P-1$ opponents. Each corresponding column of the repeats-array stores $\textrm{Reps}\left(\@vec{c}^{j}\right)$, and the same column of the payoffs-array stores, for each action $a$, the payoff $v_{a}(\@vec{c}^{j})$ for a player choosing $a$ when their opponents play according to $\@vec{c}^{j}$.

This also simplifies the deviation-payoffs calculation. We compute probability $p^{j}=\Pr(\@vec{c}^{j})$ of a configuration in same way as before, except we can skip the steps from Equations 7 and 8, instead simply accessing configuration $\@vec{c}^{j}$ from storage:

[TABLE]

The key difference in this equation is that we compute the probability of a configuration $j$, instead of the probability of a configuration derived from a masked profile $i$. This change removes the need to mask payoffs and means that the configuration probabilities are identical for all deviations, meaning we do not have to perform a separate probability calculation for each action. This single array of configuration probabilities in turn simplifies vectorization over deviation actions, easily giving us the whole deviation-payoff vector at once. To get deviation payoffs, we multiply the configuration probabilities $p^{j}$ by the payoffs-array (broadcasting over actions), and sum over the configurations dimension.

[TABLE]

Here $\@vec{v}^{j}$ refers to a column of the payoffs-array, and summing these configuration-payoff vectors gives us a vector of deviation payoffs. The code for this version is the easiest to understand, so despite our recommended further improvements, we include it in the Julia Code Appendix.666Implementation shown in Appendix B: Julia Code. https://arxiv.org/abs/2302.13232

3.5. Pre-Weighting by Repetitions

The next improvement comes from the realization that by associativity we can re-group Equation (10) as follows.

[TABLE]

Which means that even though repetitions are logically part of the probability calculation, we can simplify our computations by storing them with the payoff values. Specifically, we can combine the repeats-array into the payoffs-array by multiplying each payoff value by that configuration’s repetitions, so that the entry in row $a$, column $j$ of the payoffs array stores the value:

[TABLE]

This operation can be performed once in advance, saving space and speeding up all subsequent deviation payoff calculations.55footnotemark: 5

3.6. Log Transformation

The combined effect of all the improvements so far is sufficient to allow reasonably efficient deviation payoff calculations for games with over 100 players, as long as the number of actions is kept small, but this poses a new problem: $\binom{32}{6,6,6,7,7}>2^{63}$, meaning that the repetitions overflow a 64-bit integer for some profiles with $P=33$ and $A=5$. We can solve this problem by working with log-probabilities, and this incidentally produces a slight computational speed-up by making use of arithmetic operations that consume fewer processor cycles when we transform exponents into multiplication and multiplication into addition. Specifically, we can store the natural (or other base) log of the repetition-weighted payoffs $\lambda$, and to calculate repetitions, we can use a log-gamma function to avoid overflows.

[TABLE]

This will not work if payoffs are negative, but since any positive affine transformation of utilities has no effect on incentives, we can transform any game into one with non-negative payoffs. In fact, we find it useful under all of our game representations to transform the payoffs into a standardized range, to simplify hyperparameter tuning for various equilibrium-computation algorithms. And of course these calculations can be vectorized as before with the exponential applied element-wise to an array with contribution-values for each configuration.

We can now compute deviation payoffs using this representation by first computing each configuration’s contribution $\gamma^{j}$ to the deviation payoff in log-space.

[TABLE]

and then summing over all configurations.55footnotemark: 5

[TABLE]

3.7. GPU Acceleration

Now that our deviation-payoff calculations are done as a sequence of simple mathematical operations broadcast across large arrays of floating-point data, an obvious way to accelerate them is to move the computation to a graphics processor. Most modern programming languages have libraries that make GPU translation manageable for operations as simple as ours, and specifically in Julia the translation is trivial, requiring only that we choose an array data type that moves the data to the GPU; the code for CPU and GPU versions of our final data structure is otherwise identical.66footnotemark: 6

3.8. Batch Processing

The final improvement we implemented requires no change to the preceding data structure, but instead uses it more efficiently. If memory allows, we can take even greater advantage of SIMD operations on the GPU by computing deviation payoffs for a batch containing multiple mixtures.66footnotemark: 6 This is useful because many of the algorithms we use to compute Nash equilibria are forms of local search where it can help to initialize them from many starting points. For example, in both replicator dynamics and gain-gradient descent, we generally restart the algorithm as many as 100 times from different initial mixtures, so by parallelizing deviation-payoff calculations, we can parallelize multiple restarts of the algorithm. This adds an extra dimension to the intermediate arrays constructed by the deviation payoff computation, giving them dimension actions $\times$ configurations $\times$ mixtures. This is mainly helpful with small games, since large games can occupy all GPU execution units with the deviation-payoff calculation for a single mixture. If batch processing is used on large games, the number of mixtures must be carefully chosen to not exceed available GPU memory.

3.9. Deviation Derivatives

For any of the data structure variants we describe, it is possible to calculate deviation derivatives, but here we describe the approach only for the final log-transformed variants.

Computing the Jacobian with respect to the mixture probability $\@vec{\sigma}_{s}$ only requires a small addition to the deviation payoff contributions computed in equation (12); we multiply by $\frac{\@vec{c}_{s}}{\@vec{\sigma}_{s}}$ over a new dimension $s$, before summing across configurations.66footnotemark: 6

[TABLE]

This calculation can be fully vectorized using similar operations but adds a dimension since we have a partial derivative for each pair of actions. This extra dimension is the main reason gradient descent runs slower than replicator dynamics in our experiments.

4. Validation and Experiments

4.1. Size Limits

The first constraint on the size of the data structures comes from integer overflows when calculating repetitions. Table 1 shows the upper limit in terms of $P$ and $A$ beyond which at least one profile will overflow a signed 64-bit integer. For example, if we want to represent a 5-action game with more than 32 players, we must use the log-transform to ensure that we calculate probabilities correctly. Note that using unsigned integers usually makes no difference and never allows more than one extra player.

This raises the question of whether we sacrifice precision when using the log-transform. To test this, we generated 10 random bipartite action-graph games for every combination of $2\leq P\leq 512$ and $2\leq A\leq 20$ for which the size of our data structures would be below 1GiB. We then represented these games using the CPU (64-bit) and GPU (32-bit) versions of our data structure, and calculated deviation payoffs for 1000 mixtures per game and recorded the largest error in the calculated deviation payoffs. These errors were tiny, except in games with hundreds of players in the 32-bit representation, and even then the scale of the errors was at least 7 orders of magnitude smaller than the payoff range of the game. Full results are shown in the Supplementary Figures Appendix.777See Appendix D: supplementary figures. https://arxiv.org/abs/2302.13232

With these results, the question of what games we can represent comes down to what can fit in (main or GPU) memory. Figure 3 shows the size of the GPU data structure as $P$ increases for various values of $A$. A similar figure for 64-bit arrays appears in Appendix D.77footnotemark: 7 Note that calculating deviation payoffs requires storing intermediate arrays of similar size to the data structure, so game size should be restricted to below half of available memory (or smaller if we want to calculate deviation payoffs in batches). Note that these sizes mean that it is entirely reasonable to store a 100-player, 6-action game.

4.2. Timing Comparisons

Figures 4 and 5 show our main results: our data structures produce significant speedup in computing deviation payoffs and therefore running Nash-finding algorithms. Timing comparisons were run using an AMD Threadripper 2990WX CPU and NVIDIA RTX 2080-Ti GPUs. Figure  4 shows the time to compute deviation payoffs using each of the data structures described in the previous section. The first several lines end at $P=32$, since this is near the int-overflow limit for $A=4$. The batch-processing line ends at $P=128$, because operating on a batch of 64 mixtures requires much more memory. Note that our proposed data structure can compute deviation payoffs in a game with 384 players faster than the baseline data structure can handle $P=12$. Many more deviation-payoff timing comparisons including 6- and 8-action games appear in Appendix D.77footnotemark: 7

Figure 5 shows a similar result for Nash-finding, but focusing only on our best data structure, and using batch sizes adapted to available memory. First, note that replicator dynamics using our data structure outperforms the previous implementation by four orders of magnitude. Gradient descent is consistently an order of magnitude slower than replicator dynamics, making it impractical and not generally implemented with older data structures. But with our data structure that slowdown can often be acceptable because gradient descent and replicator dynamics frequently identify some distinct equilibria (as illustrated in Figure 1), and so in cases where replicator dynamics fails or more equilibria are desired it is entirely reasonable to run both algorithms. For example, with $A=4$ and $P=512$, we can perform 100 runs of replicator dynamics or gradient descent in a reasonable time-frame. Timing results shown in Appendix D for $A=6$ and $A=8$ are broadly similar.77footnotemark: 7

5. Extensions

5.1. Many Strategies

The size of our data structures and therefore the time to compute deviation payoffs scales reasonably similarly with respect to the number of players $P$ and the number of actions $A$. The configurations and payoffs arrays both have size $A\times\binom{P+A-2}{P-1}$, which is dominated by the binomial that grows exponentially in the smaller of $P$ and $A$. As a result, it is generally possible to represent games with a large number of players but a small number of actions or a large number of actions but a small number of players, with the caveat that actions contribute an extra linear factor.

Unfortunately, equilibrium computation does not scale equally well in players and actions. Identifying symmetric mixed-strategy Nash equilibria means searching for specific points in the mixed-strategy simplex, whose dimension grows with the number of actions. The various complete algorithms are generally hopelessly slow with large numbers of actions, and the local search methods we prefer are much less likely to identify an equilibrium on any given run. This means that solving symmetric games with many players is far more feasible than solving ones with many actions.

To handle games with many actions, we can take inspiration from Porter et al. (2008), and search for equilibria with small support sets. This requires solving a subgame defined by restricting the game to a subset of the actions, and then checking whether those subgame solutions have beneficial deviations in the full action set. To facilitate this, we can use a variant of our data structure where we include additional rows in the payoffs-array representing actions outside the support set, but make no change to the configurations-array. Then for a mixture in the support set, we can compute deviation payoffs for all actions against that mixture by expanding the step where we broadcast the multiplication of configuration probabilities and payoffs. If we determine that the current subgame set does not contain a Nash equilibrium, we can update the support set replacing weaker strategies with stronger ones, as measured by deviation gain against the subgame’s candidate solutions (or other metrics). This will iteratively search the space of supports and eventually identify a small-support equilibrium if any exist. This approach of constructing complete subgames with extra payoff information for unilateral deviations is already widely used in empirical game settings Wellman (2020) where determining payoffs for a single configuration is expensive.

5.2. Role-Symmetric Games

Outside of fully-symmetric games, many large normal-form games exhibit role-symmetry, where players can be grouped into $R$ roles—like buyers and sellers or attackers and defenders—and within each role, all players are indistinguishable. To solve role symmetric games, analysts typically search for role-symmetric Nash equilibria, where all players in a given role play the same mixed strategy. Role symmetric mixed-strategy profiles can be represented by a vector that concatenates the mixed strategy played by each role, and the deviation payoff vector of the same dimension gives the expected utility of a unilateral deviator selecting each action. These deviation payoffs and their derivatives remain central to computing role-symmetric equilibria, and almost all of the techniques we propose for improving representations of symmetric games also apply under role symmetry.

Unfortunately, representing opponent configurations gets trickier with more than one role, because the set of opponents is different for players belonging to different roles. For example, in a game with 10 buyers and 10 sellers, when computing deviation payoffs for buyer-actions, opponent configurations include 9 buyers and 10 sellers, whereas for seller actions, there are 10 buyer-opponents and 9 seller-opponents. We could resolve this by building separate opponent-configuration tables for each role $r$, where the arrays for role $r$ are based on configurations with $P_{r}-1$ opponents in the same role and $P_{r^{\prime}}$ opponents in other roles. This gives us $R$ pairs of configuration- and payoff-arrays where the arrays for role $r$ have $\binom{P_{r}+A_{r}-2}{P_{r}-1}\left(\prod_{r^{\prime}\neq r}\binom{P_{r^{\prime}}+A_{r^{\prime}}-1}{P_{r^{\prime}}}\right)$ columns, by $\left(\sum_{r=1}^{R}A_{r}\right)$ rows for configurations and $A_{r}$ rows for payoffs.

We could also dispense with the opponent-configuration approach and instead store full $P$-player profiles. This results in one profile-array and one payoff-array that each have size $\left(\sum_{r=1}^{R}A_{r}\right)\times\left(\prod_{r=1}^{R}\binom{P_{r}+A_{r}-1}{P_{r}}\right)$, but requires us to return to masking-based deviation-payoff computations. Thus for multi-role games, we have a choice between slightly smaller storage requirements for a profile-based data structure or slightly simpler computations for a configuration-based representation; the gameanalysis.py library22footnotemark: 2 employs the former option. Under either of these approaches, all of our other optimizations can still apply, and as long as the number of roles is small both options provide for reasonably efficient deviation-payoff calculations.

5.3. Action-Graph Games

Action-graph games Jiang et al. (2011) represent certain games compactly by storing, for each node, a mapping from neighborhood configurations to payoffs. Each action (for any player) corresponds to a node in a graph, and an action’s payoff depends only on the number of players choosing each of the adjacent actions. This means that action-graph games are role-symmetric, with any group of players who share an action set belonging to the same role; when all players have the same action set, action-graph games are symmetric.

In the symmetric case, we can extend our data structures to efficiently compute deviation payoffs in action-graph games with a pair of arrays for each action storing that action’s neighborhood configurations and repetition-weighted payoffs. The key augmentation is to include, as part of the representation of an opponent-configuration, an extra “action” capturing the number of players choosing any action outside the neighborhood. Then when computing the deviation payoff for an action, we can sum the probabilities for each non-adjacent action to get the out-of-neighborhood probability. This means that all of our data structure improvements can be applied to symmetric AGG-$\emptyset$ games.

For role-symmetric AGG-$\emptyset$ games, our role-symmetric data structure variants could be applied after splitting up action nodes shared by more than one role. For action-graph games with contribution-independent function nodes, many of our data-structure improvements can be applied, but the representation of these games tends to be sufficiently compact that the potential gains may be small. The gameanalysis.jl library11footnotemark: 1 implements symmetric action-graph games with a bipartite graph between actions and contribution-independent functions as one of the tools for generating interesting random symmetric game instances. This implementation makes use of vectorization, opponent configurations, pre-computed repetitions, and the log transformation.

5.4. Multi-Threading

The next extension we would like to explore would focus on improving the CPU version of our data structure to capture some of the parallelism achievable on the GPU. At present, our deviation-payoff calculations are single-threaded, and while we might not expect to outperform SIMD array operations, there should still be room for significant speedup on multi-core systems from parallelizing a batch of deviation payoff computations across threads, especially for games in the multi-gigabyte size-range that can stress GPU memory but still fit in system RAM.

6. Conclusion

Our validation experiments show that the data structures we propose are capable of storing (Section 4.1) and efficiently computing (Section 4.2) deviation payoffs in symmetric games with very large numbers of players. By using incomplete-but-effective search algorithms we are consistently able to identify symmetric Nash equilibria. In combination with the iterative exploration approach described in Section 5.1 we can also find small-support equilibria in games with many actions. These results dramatically out-class existing solvers, enabling analysis of a much wider range of symmetric games, and closing the gap between micro- and macro-scale symmetric games.

Appendix A Appendix A: Worked Examples

We will illustrate the various iterations of our data structure and the corresponding deviation payoff calculations using the following example of a 3-player, 3-action symmetric game. In this game, each player selects an action $a\in\{1,2,3\}$, and if a player’s action $a$ is unique, their payoff is $a$. If two players select action $a$, they both receive payoff $-a$; if all three players select $a$, they each receive 0. The standard normal-form representation of this game is shown in Figure 6.

For each data structure, we will show the calculation of $\@vec{u}_{1}(\@vec{\sigma})$ for $\@vec{\sigma}=\left<0.1,0.5,0.4\right>$, that is the deviation payoff of action 1 when both opponents play a mixture that places a 10% probability on 1, a 50% probability on 2 and a 40% probability on 3. In cases where the computation naturally vectorizes over actions, we will show how to compute $\@vec{u}(\@vec{\sigma})$.

A.1. Payoff Dictionary

If we represent our example game using a mapping from profiles to payoffs, the mapping will store the following $\binom{P+A-1}{P}=10$ pairs. We denote extraneous entries in the payoff vector by a dot.

[TABLE]

To calculate $\@vec{u}_{1}(\@vec{\sigma})$, we begin from Equation (2), giving

[TABLE]

So this computation loops over all six of the profiles where action 1 is played. Calculating $\@vec{u}_{2}(\@vec{\sigma})$ and $\@vec{u}_{3}(\@vec{\sigma})$ requires similar loops over different subsets of the profiles.

A.2. Array Vectorization

If we use parallel arrays to store the profiles and payoffs, our example game is stored in the following $3\times 10$ arrays. Note that each column of the tables corresponds to a profile/payoff pair from the dictionary version.

Profiles:

3

2

2

1

1

1

0

0

0

0

0

1

0

2

1

0

3

2

1

0

0

0

1

0

1

2

0

1

2

3

Payoffs:

0

-1

-1

1

1

1

$\cdot$

$\cdot$

$\cdot$

$\cdot$

$\cdot$

2

$\cdot$

-2

2

$\cdot$

0

-2

2

$\cdot$

$\cdot$

$\cdot$

3

$\cdot$

3

-3

$\cdot$

3

-3

0

The deviation-payoff calculation $\@vec{u}_{1}(\@vec{\sigma})$ is fundamentally the same as in the previous section, but performed with element-wise array operations instead of explicit loops. For clarity, we will show some of the intermediate steps here. First, we mask out the entries where action 1 is not played, decrement the count for action 1, and use the resulting profiles to exponentiate the probabilities from $\@vec{\sigma}$.

[TABLE]

$\mbox{\LARGE.}^{\wedge}$

2

1

1

0

0

0

0

1

0

2

1

0

0

0

1

0

1

2

$=$

​​.01

.1

.1

1

1

1

1

.5

1

​​.25

.5

1

1

1

.4

1

.4

​​.16

Where the $\mbox{\LARGE.}^{\wedge}$ operator indicates element-wise exponentiation broadcast across the rows: that is, element $r,c$ of the output matrix comes from raising element $r$ of the input vector to the power of element $r,c$ in the input matrix. Then we can get the probability of each profile by taking a product down each column and multiplying element-wise by an array of repeats, which we compute by applying a multinomial function to each column of the same masked-and-decremented profiles array.

[TABLE]

$\;=\;$

​​.01

.1

​​.08

​​.25

.4

​​.16

We can then apply the same mask to the payoffs array, slice out the row for action 1, element-wise multiply with the profile-probabilities array, and sum to get the deviation payoff.

[TABLE]

[TABLE]

\Biggr{)}\;=\;0.63

A.3. Pre-computed Repetitions

Our next improvement came from pre-computing $\textrm{Reps}\left(\@vec{s}\mid a\right)$ for each profile-action pair, and storing them in a third parallel array. Here is the resulting array for our example game.

Repeats:

1

2

2

1

2

1

$\cdot$

$\cdot$

$\cdot$

$\cdot$

$\cdot$

1

$\cdot$

2

2

$\cdot$

1

2

1

$\cdot$

$\cdot$

$\cdot$

1

$\cdot$

2

2

$\cdot$

1

2

1

The deviation-payoff calculation for a single action using this version would identical to the one just described, except that this array of configuration-repetitions:

[TABLE]

would now be simply a masked slice out of the repeats-array, avoiding the need to re-calculate multinomials for every deviation-payoff computation. But we also discussed how to compute the full vector of deviation payoffs $\@vec{u}(\@vec{\sigma})$ at once. The computation begins by raising the mixture probabilities to each of the profile counts.

[TABLE]

$\mbox{\LARGE.}^{\wedge}$

3

2

2

1

1

1

0

0

0

0

0

1

0

2

1

0

3

2

1

0

0

0

1

0

1

2

0

1

2

3

$\;\;=\;\;$

​​​.001

​​.01

​​.01

.1

.1

.1

1

1

1

1

1

.5

1

​​.25

.5

1

​​​.125

​​.25

.5

1

1

1

.4

1

.4

​​.16

1

.4

​​.16

​​​.064

We can then take a product down each column, giving an array of the $\prod_{a}$ terms from Equation (6). We then want to multiply this by the repeats-array (with zeros filled in for the missing elements) element-wise divided by the mixture probabilities:

[TABLE]

Note that the first six entries in the top row of this array correspond to the $\Pr_{\@vec{\sigma}}(\@vec{s}\mid 1)$ array from before. Likewise, the non-zero entries in the second and third rows give the various $\Pr_{\@vec{\sigma}}(\@vec{s}\mid 2)$ and $\Pr_{\@vec{\sigma}}(\@vec{s}\mid 3)$ probabilities. We can now element-wise multiply by the payoffs array (again filling in zeros for missing entries), then sum across the rows to get the vector of deviation payoffs.

[TABLE]

[TABLE]

[TABLE]

A.4. Opponent Configurations

If we move to storing arrays based on configurations over $P-1$ opponents instead of profiles over $P$ players, we shrink our arrays to $\binom{P+A-2}{P-1}=6$ columns, and can eliminate the actions dimension from the repeats-array. The resulting tables for our example game are as follows.

Configurations:

2

1

1

0

0

0

0

1

0

2

1

0

0

0

1

0

1

2

Payoffs:

0

-1

-1

1

1

1

2

-2

2

0

-2

2

3

3

-3

3

-3

0

Repeats:

1

2

2

1

2

1

With this representation, computing a vector of deviation payoffs is simplified considerably. We begin once again by raising mixture probabilities to the configuration counts.

[TABLE]

$\mbox{\LARGE.}^{\wedge}$

2

1

1

0

0

0

0

1

0

2

1

0

0

0

1

0

1

2

$=$

​​.01

.1

.1

1

1

1

1

.5

1

​​.25

.5

1

1

1

.4

1

.4

​​.16

Then we take the product down each column and multiply column-wise by the repeats-array.

[TABLE]

$\;=\;$

​​.01

.1

​​.08

​​.25

.4

​​.16

Note that while the steps so far appear identical to those given for payoff-arrays above, they are simpler because no masking is required in the configurations array and the repeats-array has been pre-computed. From here, we column-wise multiply by the payoffs array and then sum across rows.

[TABLE]

[TABLE]

A.5. Weighted Payoffs

The next trick was to avoid storing the repeats-array by re-arranging the multiplication. By associativity, we can multiply the repeats by the payoffs instead of by the probabilities. This can be done in advance, with the result that we store the following weighted payoffs.

Weighted Payoffs:

0

-2

-2

1

2

1

2

-4

4

0

-4

2

3

6

-6

3

-6

0

The deviation payoff calculation proceeds identically to above, except that we can skip one step and the final operation before the sum becomes:

[TABLE]

A.6. Log Probabilities

The last major change to the data structure allows us to work with log-probabilities. Unfortunately, we cannot simply log-transform the weighted payoffs array, because it often has negative entries. However, since a positive affine transformation of a game’s payoffs has no effect on incentives, we can re-scale the payoffs to permit this transformation. Somewhat arbitrarily, we chose to re-scale the payoffs so that the minimum value becomes $10^{-5}$ and the maximum value becomes $10^{3}$. Applying this transformation to the opponent-configuration version of the payoff array gives:

Scaled Payoffs:

500.00

333.33

333.33

666.67

666.67

666.67

833.33

166.67

833.33

500.00

166.67

833.33

1000.0

1000.0

1.0e-5

1000.0

1.0e-5

500.00

We can then multiply by the repeats and take the log-transform, resulting in the following array, which is stored alongside the opponent-configurations array.

Weighted Log-Payoffs:

6.215

6.502

6.502

6.502

7.195

6.502

6.725

5.809

7.419

6.215

5.809

6.725

6.908

7.601

​​-10.82

6.908

​​-10.82

6.215

Using the log-transform, the multiplication operations in our deviation payoff calculation become additions and the exponents become multiplications. As a first step, we need to take the natural log of our mixture probabilities.

[TABLE]

Then we column-wise multiply the log probabilities by the configurations.

[TABLE]

$\mbox{\LARGE.}\times$

2

1

1

0

0

0

0

1

0

2

1

0

0

0

1

0

1

2

$\;\;=\;\;$

​​-4.605

​​-2.303

​​-2.303

0

0

0

0

​​-0.693

0

​​-1.386

​​-0.693

0

0

0

​​-0.916

0

​​-0.916

​​-1.833

Next we sum down each column and add the weighted payoffs.

[TABLE]

Now we can exponentiate to undo the log transform and then sum over the rows.

[TABLE]

$\;=$

​5.0000

​33.333

​26.667

​166.67

​266.67

​106.67

​8.3333

​16.667

​66.667

​125.00

​66.667

​133.33

​10.000

​100.00

​8.0e-7

​250.00

​4.0e-6

​80.000

[TABLE]

This leaves us with deviation payoffs for the re-scaled game from $10^{-5}$ to $10^{3}$, which are fine (in fact preferable) for computing equilibria, but if we want payoffs from the original game, we can undo the affine transformation, resulting in the usual answer of $\left[0.63\;-\!0.50\;-\!0.36\right]^{\mathrm{T}}$.

A.7. Further Improvements

Subsequently, we proposed GPU acceleration, which performs a mathematically identical calculation to that shown for log-probabilities above, but takes advantage of SIMD instructions. We also proposed operating on a batch of mixed strategies stored as a 2-dimensional array. This is harder to show in a paper-and-pencil example, because it adds a third dimension to several of the arrays, and really amounts only to doing several instances of the same computation using instructions broadcast along that extra dimension. The Julia code appendix shows how the implementations differ.

Appendix B Appendix B: Julia Code

This appendix shows a small subset of the code implementing our data structures and Nash computation algorithms, focusing on those parts that best illustrate the key concepts discussed in the paper. For more complete (and runnable) implementations, see the code supplement.

B.1. Data Structure Implementations

First, we show the deviation_payoffs function for the data structure variant from the Opponent Configurations section. This is not the best version of the data structure, but it most-clearly exhibits the core idea of the vectorized deviation-payoff calculation. The @reduce macro comes from the TensorCast.jl package.

function deviation_payoffs(game::SymGame_OppConfig, mixture::AbstractVector) @reduce probs[c] := prod(a) mixture[a] ^ game.config_table[a,c] @reduce dev_pays[a] := sum(c) game.payoff_table[a,c] * probs[c] * game.repeat_table[c] return dev_pays end

Next we show the struct definition for the best version of the data structure, called SymGame_GPU. The offset and scale attributes are saved so that if necessary we can undo the affine transformation used to normalize the payoffs and present payoffs on their original scale. Of note, all of the code for the SymGame_CPU variant is identical, except that it uses Arrays in place of CUDA.CuArrays.

struct SymGame_GPU <: SymmetricGame num_players::Int32 num_actions::Int32 config_table::CUDA.CuArray{Float32,2} payoff_table::CUDA.CuArray{Float32,2} offset::Float32 scale::Float32 end

The constructor for the SymGame_GPU struct takes a payoff_generator function as input. This is a function that, when passed a configuration $\@vec{c}$, returns a vector of payoffs $\@vec{v}$ containing $v_{a}(\@vec{c})$ for each action $a$. CwR is an alias for Combinatorics.with_replacement_combinations. The call to logmultinomial uses SpecialFunctions.loggamma to compute the log of the multinomial from Equation (1). The call to set_scale determines the positive affine transformation that will normalize the payoffs to the range const MIN_PAYOFF = 1e-5 through const MAX_PAYOFF = 1e3.

function SymGame_GPU(num_players, num_actions, payoff_generator) num_configs = multinomial(num_players-1, num_actions-1) config_table = Array{Float32}(undef, num_actions, num_configs) payoff_table = Array{Float32}(undef, num_actions, num_configs) log_repeat_table = Array{Float32}(undef, 1, num_configs) for (c,config) in enumerate(CwR(1:num_actions, num_players-1)) prof = counts(config, 1:num_actions) config_table[:,c] = prof log_repeat_table[c] = logmultinomial(prof...) payoff_table[:,c] = payoff_generator(prof) end (offset, scale) = set_scale(minimum(payoff_table), maximum(payoff_table)) payoff_table = log.(normalize(payoff_table, offset, scale)) .+ log_repeat_table SymGame_GPU(num_players, num_actions, config_table, payoff_table, offset, scale) end

Next, we show the function that computes deviation payoffs on a batch of mixtures. This differs from the OppConfig version above in that most operations are performed under a log transformation, and that the computation is being broadcast across an additional mixture-batch dimension. Note that the returned deviation payoffs are in the normalized scale, because this helps with hyperparameter tuning for various Nash algorithms. A call to denormalize(game, payoffs) can transform them back to the game’s original payoff scale.

function deviation_payoffs(game::SymGame_GPU, mixture_batch::CUDA.CuArray{Float32,2}) log_mixtures = log.(mixture_batch .+ eps32) @reduce log_probs[m,c] := sum(a) log_mixtures[a,m] * game.config_table[a,c] @reduce dev_pays[a,m] := sum(c) exp(game.payoff_table[a,c] + log_probs[m,c]) return Array(dev_pays) end

Finally, we show the function to compute deviation derivatives for a single mixture. Note that the mixture .+ eps32 operation ends up only changing array entries equal to zero, and is slightly more efficient than explicitly filtering for those entries.

function deviation_derivatives(game::SymGame_GPU, mixture::CUDA.CuArray{Float32,1}) mixture = mixture .+ eps32 log_mixture = log.(mixture) @reduce log_probs[c] := sum(a) log_mixture[a] * game.config_table[a,c] deriv_configs = game.config_table ./ mixture @reduce J[a,s] := sum(c) exp(game.payoff_table[a,c] + log_probs[c]) * deriv_configs[s,c] return J end

The code supplement implements many additional features for this and other data structure variants, including batched deviation derivatives, single-mixture deviation payoffs, and a number of other convenience functions.

B.2. Nash Algorithm Implementations

Here, we present code for the two algorithms we have found to be the most practical for identifying Nash equilibria in large symmetric games: replicator dynamics and gradient descent on sum-of-gains. Both of these functions can operate on either a vector representing a single mixture or an array representing a batch of mixtures, thanks to multiple dispatch to vector- or matrix-versions of the deviation_payoffs and gain_gradients functions, respectively.

For replicator dynamics, the offset input allows for the possibility that the game’s payoffs are not normalized, in which case it should be less than or equal to the minimum payoff in the game to prevent negative probabilities. This parameter can also be used as a quasi-step-size, since a larger negative offset will diminish the relative magnitude of payoff differences. The simplex_normalize function simply divides by the sum to ensure probabilities add to 1.

function replicator_dynamics(game::SymmetricGame, mix::AbstractVecOrMat; iters::Integer=1000, offset::Real=0) for i in 1:iters mix = simplex_normalize(mix .* (deviation_payoffs(game, mix) .- offset)) end return mix end

For gain descent, we allow either a fixed step size or schedule of step sizes. The simplex_normalize function implements the projection algorithm given by Wang and Carreira-Perpinán (2013), and also operates on either a single mixture or a batch.

function gain_descent(game::SymmetricGame, mix::AbstractVecOrMat; iters::Integer=1000, step_size::Union{Real,AbstractVector}=1e-6) if step_size isa Real step_size = [step_size for i in 1:iters] end for i in 1:iters mix = simplex_project(mix .- step_size[i]*gain_gradients(game, mix)) end return mix end

Here is the implementation of the gain_gradients function for a single mixture. The version operating on a batch of mixtures also appears in the code supplement.

function gain_gradients(game::SymmetricGame, mixture::AbstractVector) dev_pays = deviation_payoffs(game, mixture) mixture_EV = mixture’ * dev_pays dev_jac = deviation_derivatives(game, mixture) util_grads = (mixture’ * dev_jac)’ .+ dev_pays gain_jac = dev_jac .- util_grads’ gain_jac[dev_pays .< mixture_EV,:] .= 0 return dropdims(sum(gain_jac, dims=1), dims=1) end

The code supplement also contains implementations of fictitious play and iterated better-response, along with a number other functions including regret-filtering and deduplication for processing candidate equilibria returned by these methods.

Appendix C Appendix C: Nash Algorithms

In this appendix we provide a description of a few additional Nash-finding algorithms, formulated in terms of deviation payoffs. Note that this is by no means a complete survey of Nash-finding algorithms.

Fictitious play

Brown (1951) operates by iteratively best-responding to the opponents’ empirical distribution of past actions as though that represented a mixed strategy. In the context of a symmetric game, we can formulate FP as keeping a count of the number of times each action has been “played”, and normalizing those counts to produce a mixed strategy. Then on each iteration, the count for the best response to the current mixture is “played”, incrementing that action’s count. Because we need to specify some initial distribution for the first step, we begin with a mixture $\@vec{\sigma}^{0}$, and generate the initial counts based on some constant multiple of those probabilities: $\@vec{w}\leftarrow c\cdot\@vec{\sigma}^{0}$.

[TABLE]

Scarf’s simplicial subdivision

algorithm Scarf (1967) is a generic method for finding the fixed point of a continuous function on the simplex. Since the existence of Nash equilibria is proven by constructing a continuous function whose fixed-points correspond to equilibria, we can apply Scarf’s algorithm to that function. Similar to the sum-of-gains $g$ we defined for gradient descent, we can define the gain for a particular action, and use it to construct the better-response function $B$.

[TABLE]

The function $B$ outputs a mixture that increases the relative probability of any actions that have deviation payoffs better than the expected utility of $\@vec{\sigma}$. Nash equilibria of the game correspond exactly to fixed-points of $B$, so we can use Scarf’s algorithm to identify a fixed point and therefore a Nash equilibrium. Scarf’s algorithm traces a path through sub-spaces of the simplex to identify one that must contain a fixed point. Then within the identified sub-simplex, further subdivisions can be searched recursively to narrow in on a fixed point.

An even simpler algorithm that makes use of $B$ is to simply apply it repeatedly. We refer to this algorithm as iterated better-response. This algorithm tends to behave like a smoothed version of fictitious play, because it up-weights several better responses instead of just a single best response. Figure 13 in the Appendix D shows traces of FP and IBR on the same game. Note that these algorithms follow very similar trajectories, except that FP always moves in a straight line toward one of the corners of the simplex, since it only up-weights one action, while IBR can curve, by up-weighting multiple actions at once, and by slightly different amounts on consecutive iterations.

For all of the incomplete algorithms described (replicator dynamics, gradient descent, fictitious play, and iterated better-response), we typically run them for a fixed number of iterations, from a grid of equally-spaced starting mixtures. The resulting mixtures may not be Nash equilibria, due to non-convergence or local-optimality. We therefore filter out any endpoints with regret above a predetermined $\varepsilon$. It is also common for several runs to converge to essentially the same $\varepsilon$-equilibrium, so we filter the results to eliminate “duplicate” points within some threshold for Euclidean or Chebyshev distance.

Appendix D Appendix D: Supplementary Figures

In this appendix we provide, for completeness, a number of slight variations on the figures appearing in the main text.

Most of the simplex plots shown here are based on random additive-sine games. We generate these games as bipartite action-graph games with 200 function nodes. Each action will be connected to a random subset of the function nodes by an Erdős-Renyi random-graph process. Then each function node is connected to every action with a random weight. Function nodes are contribution-independent and simply count the number of agents in their neighborhood. They then output a value based on a random sin-function added to a random low-degree polynomial. Each action then computes its value as a weighted sum over the function outputs.

However, our data structures make it easy to specify even more interesting classes of random games. As an example, we generated random Gaussian mixture games for some of our experiments. To generate a GMG, we specify a number of functions per action, and for each, we choose a random centroid from a Dirichlet distribution, scaled to the size-$P$ integer simplex, a random covariance matrix from a Lewandowski-Kurowicka-Joe distribution, and a random multiplicative scaling factor (since these Gaussians need not represent valid probability distributions). Then to determine the payoffs for a given configuration/action pair, we evaluate each Gaussian function for that action at the given configuration, multiplied by the scaling factor, and sum the results. An example regret heatmap for a 3-action Gaussian mixture game is shown in Figure 14.