# Efficient model predictive control for nonlinear systems modelled by deep
neural networks

Jianglin Lan<EMAIL_ADDRESS>James Watt School of Engineering,
University of Glasgow, Glasgow G12 8QQ, United Kingdom

###### Abstract

This paper presents a model predictive control (MPC) for dynamic systems whose
nonlinearity and uncertainty are modelled by deep neural networks (NNs), under
input and state constraints. Since the NN output contains a high-order complex
nonlinearity of the system state and control input, the MPC problem is
nonlinear and challenging to solve for real-time control. This paper proposes
two types of methods for solving the MPC problem: the mixed integer
programming (MIP) method which produces an exact solution to the nonlinear
MPC, and linear relaxation (LR) methods which generally give suboptimal
solutions but are much computationally cheaper. Extensive numerical simulation
for an inverted pendulum system modelled by ReLU NNs of various sizes is used
to demonstrate and compare performance of the MIP and LR methods.

###### keywords:

Deep neural network, linear relaxation, mixed integer program, model
predictive control, nonlinear system.

††thanks: This work was supported by the Leverhulme Trust Early Career
Fellowship under Award ECF-2021-517.

##  1 Introduction

The advancement of deep learning techniques has stimulated much interest in
adopting neural networks (NNs) to power autonomous systems such as robots and
self-driving cars (Spielberg et al., 2019; Tang et al., 2022). In particular,
NNs have been demonstrated to be powerful in the modelling and control of
dynamic systems Moe et al. (2018). This paper focuses on NN-modelled control
systems, where deep NNs are used to model the nonlinear system dynamics and/or
uncertainties. Currently, the modelling is achieved by two main types of NNs:
static networks (e.g., multi-layer perception NNs) and dynamic networks (e.g.,
recurrent NNs, neuro-fuzzy NNs). The highly nonlinear nature of NNs imposes
big challenges on the control design of NN-modelled dynamic systems. The
design becomes even more challenging when using non-differentiable NNs such as
the ReLU (Rectified Linear Unit) NNs. This calls for control designs tailored
to NN-modelled dynamic systems.

There are a few existing control designs for NN-modelled dynamic systems. For
dynamic systems modelled with a single hidden layer NN, an internal model
control is designed in Nahas et al. (1992) and state feedback controllers are
developed in Nikolakopoulou et al. (2022); D’Amico et al. (2023) using the
linear matrix inequality technique. An iterative linear quadratic regulator is
proposed in Nagariya & Saripalli (2020) for dynamic systems modelled by a deep
NN. However, all the above designs consider only unconstrained controllers.

To provide optimal control actions under system constraints, model predictive
control (MPC) for NN-modelled dynamic systems has attracted much attention
(Ren et al., 2022). This line of research is of particular interest because
the combination of MPC and NN models opens the door for optimally control
complex dynamic systems, which is otherwise difficult or even impossible by
using first-principle models.

The MPC design for NN-modelled dynamic systems is non-trivial. In general, MPC
requires solving constrained optimisation problems, but this can be burdensome
for NN-modelled dynamics because the NN output normally includes high-order
nonlinearity of its input, the system state and control input. This is
particularly true for large-scale (deep) NNs with many hidden layers and
neurons. Several methods have been adopted to solve the nonlinear MPC for NN-
modelled dynamic systems, e.g., the quasi-Newton algorithm (Sørensen et al.,
1999), sequential quadratic programming method (Saint-Donat et al., 1991),
real-time-iteration scheme (Ławryńczuk, 2010; Patan, 2018; Wysocki &
Ławryńczuk, 2015), sampling method (Spielberg et al., 2021; Salzmann et al.,
2023), interior point line search filter method (Rokonuzzaman et al., 2021),
recurrent NN method (Yan & Wang, 2012), and the state observer method (Bonassi
et al., 2024). However, these methods are designed for small-scale NNs having
only one or two hidden layers. To design MPC for dynamic systems modelled by
deep NNs, sampling-based methods are developed in Williams et al. (2016);
Askari et al. (2022) to obtain approximate solutions through an interative
solving procedure.

This paper presents novel efficient solving methods for MPC tracking control
of deep ReLU NN-modelled nonlinear dynamic systems. The main contributions are
as follows:

  1. 1)

A dual-mode MPC is proposed to ensure output of the deep NN-modelled dynamic
systems track the given reference. Ways of determining the steady-state target
and terminal constraint set are detailed, which are lacking in the literature
(Williams et al., 2016; Askari et al., 2022; Wei & Liu, 2022).

  2. 2)

Three different methods are developed to solve the nonlinear MPC problems,
including the 1) MIP method by representing equivalently the NN activation
function as a set of mixed integer linear constraints, 2) linear relaxation
(LR) method by replacing the activation function by a triangle over-
approximation of each ReLU, and 3) enhanced LR (eLR) method which refines the
MPC cost function with a penalty of the deviation between the steady-state NN
output value and the over-approximated NN output.

  3. 3)

Extensive numerical simulation and ablation study are used to demonstrate and
compare performance of the three solving methods, in terms of output tracking
accuracy and computational efficiency.

The rest of the paper is structured as follows: Section 2 describes the
problem, Section 3 presents the MPC design with the solving methods in Section
4, Section 5 reports the simulation results and Section 6 draws the
conclusions.

##  2 Problem description and preliminary

Consider a discrete-time system represented by

$\displaystyle\begin{split}x(t+1)&=Ax(t)+Bu(t)+Df_{\text{nn}}(x(t),u(t)),\\\
y(t)&=Cx(t),\end{split}$ (1)

where $x\in\mathbb{R}^{n}$, $u\in\mathbb{R}^{m}$, and $y\in\mathbb{R}^{p}$ are
the vectors of system state, control inputs, and measured outputs,
respectively. $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$,
$D\in\mathbb{R}^{n\times s}$ and $C\in\mathbb{R}^{p\times n}$ are known
constant matrices. The pair $(A,B)$ is assumed to be controllable.
$f_{\text{nn}}(x(t),u(t))\in\mathbb{R}^{s}$ is a $(n+m)$-input $s$-output
feedforward NN capturing the nonlinear dynamics.

This paper considers an $(L+1)$-layer NN (with an input layer, $L-1$ hidden
layers and an output layer) defined as

$\displaystyle\begin{split}z_{0}&=[x(t)^{\top},~{}u(t)^{\top}]^{\top},\\\
z_{i}&=\phi(\hat{z}_{i}),~{}\hat{z}_{i}=W_{i}z_{i-1}+b_{i},~{}i\in[1,L-1],\\\
f_{\text{nn}}(z_{0})&=W_{L}z_{L-1}+b_{L},\end{split}$ (2)

where $\hat{z}_{i},z_{i}\in\mathbb{R}^{n_{i}\times 1}$ are the pre-activation
and post-activation vectors at the $i$-th hidden layer with $n_{i}$ neurons,
respectively. $W_{i}\in\mathbb{R}^{n_{i}\times n_{i-1}}$ and
$b_{i}\in\mathbb{R}^{n_{i}\times 1}$ are the $i$-th layer weight matrix and
bias vector, respectively. $\phi(\cdot)$ is the ReLU function, which is among
the most commonly used activation functions for NN models. The ReLU function
is defined as $\phi(\hat{z}_{i})=\max\\{\hat{z}_{i},0\\}$, where max is
applied to each element $\hat{z}_{i,j}$, $j\in[1,n_{i}]$, of $\hat{z}_{i}$. A
graphic illustration of a ReLU neuron is provided in Fig. 1(a).

This paper aims to design a controller $u(t)$, in the form of an MPC, such
that the output $y(t)$ track a given constant reference $y_{r}$ whilst
satisfying the constraints:

$\displaystyle u\in\mathcal{U}$
$\displaystyle:=\\{u\in\mathbb{R}^{m}\mid\underline{u}\leq u\leq\bar{u}\\},$
(3a) $\displaystyle x\in\mathcal{X}$
$\displaystyle:=\\{x\in\mathbb{R}^{n}\mid\underline{x}\leq x\leq\bar{x}\\},$
(3b)

with the known constant bounds $\underline{u}$, $\bar{u}$, $\underline{x}$ and
$\bar{x}$. Since ReLU is nonlinear (piecewise and non-differentiable), the NN
output $f_{\text{nn}}$ is a nonlinear function of its input, the system state
and control input, and the order of nonlinearity increases dramatically with
the hidden layers. This will impose computational challenges for solving the
MPC problem.

The proposed control design will need the layer bounds of the NN, including
the lower and upper bounds of the pre-activation vectors ($\hat{l}_{i}$ and
$\hat{u}_{i}$, $i\in[1,L-1]$), and the output ($\underline{f}_{\text{nn}}$ and
$\overline{f}_{\text{nn}}$). Given the constraints in (3), these bounds can be
computed using the interval arithmetic method as follows:

$\displaystyle l_{0}$
$\displaystyle\\!=\\![\underline{x}^{\top},\underline{u}^{\top}]^{\top},~{}u_{0}=[\overline{x}^{\top},\overline{u}^{\top}]^{\top},$
$\displaystyle\hat{l}_{i}$
$\displaystyle\\!=\\!W_{i}^{+}l_{i-1}+W_{i}^{-}u_{i-1},~{}\hat{u}_{i}=W_{i}^{-}l_{i-1}+W_{i}^{+}u_{i-1},$
$\displaystyle l_{i}$
$\displaystyle\\!=\\!\max(\hat{l}_{i},0),~{}u_{i}=\max(\hat{u}_{i},0),~{}i\in[1,L-1],$
$\displaystyle\\!\\!\\!\\!\underline{f}_{\text{nn}}$
$\displaystyle\\!=\\!W_{L}^{+}l_{L-1}\\!+\\!W_{L}^{-}u_{L-1},\overline{f}_{\text{nn}}=W_{L}^{-}l_{L-1}\\!+\\!W_{L}^{+}u_{L-1},$
(4)

where $W_{i}^{+}=\max(W_{i},0)$ and $W_{i}^{-}=\min(W_{i},0)$.

##  3 Model predictive control design

(a) ReLU neuron

(b) MIP encoding

(c) Linear relaxation

Figure 1: Illustration of (a) the ReLU neuron $z_{i,j}=\phi(\hat{z}_{i,j})$
with $\hat{l}_{i,j}\leq\hat{z}_{i,j}\leq\hat{u}_{i,j}$ and its (b) MIP
encoding (the neuron is inactive if $\delta_{i,j}=0$ and active if
$\delta_{i,j}=1$) and (c) linear relaxation (the area between the blue and red
lines).

The proposed MPC takes a dual-mode form Rossiter (2003) with a state feedback
controller to pre-stablise the linear steady-state dyanmics of the system (1),
which has not been considered in the relevant litrature. The overall
controller is designed as

$u(t)=u_{\text{s}}(t)+u_{\text{c}}(t),$ (5)

where $u_{\text{s}}(t)$ and $u_{\text{c}}(t)$ are the steady-state feedback
controller and the constraint-enforcing controller, respectively.

Steady-state feedback controller design. When the system output $y(t)$ track
the given output reference $y_{r}$ accurately, the system (1) reaches the
steady state $(x^{*},u^{*})$ that satisfies the following relation:

$\displaystyle\begin{split}x^{*}&=Ax^{*}+Bu^{*}+Df_{\text{nn}}(x^{*},u^{*}),\\\
y^{*}&=Cx^{*},~{}x^{*}\in\mathcal{X},~{}u^{*}\in\mathcal{U}.\end{split}$ (6)

In view of (6), the target steady state $(x^{*},u^{*})$ can be solved from the
following optimisation problem:

$\displaystyle(x^{*},u^{*})$
$\displaystyle:=\arg\underset{x_{r},u_{r}}{\min}\|u_{r}\|^{2}_{R_{s}}$ s.t.
$\displaystyle\begin{bmatrix}I_{n}-A&-B\\\
C&0\end{bmatrix}\begin{bmatrix}x_{r}\\\
u_{r}\end{bmatrix}=\begin{bmatrix}Df_{\text{nn}}(x_{r},u_{r})\\\
y_{r}\end{bmatrix},$ (7a) $\displaystyle
x_{r}\in\mathcal{X},~{}u_{r}\in\mathcal{U},$ (7b)

where $R_{s}\succ 0$ is a given weight matrix. This problem is nonlinear
because (7a) depends on the ReLU function. It can be approximately solved by
using the fminsearchbnd method (D’Errico, 2023), which extends fminsearch in
Matlab to incorporate constraints. This problem can be solved exactly by using
the MIP method to be described in Section 4, but with high computational
burden. Designing more efficient methods for solving this problem exactly are
left for future study.

Subtracting (6) from (1) gives the tracking error system

$\delta x(t+1)=A\delta x(t)+B(u(t)-u^{*})+D\delta f_{\text{nn}}(t),$ (8)

where $\delta x=x-x^{*}$ and $\delta
f_{\text{nn}}=f_{\text{nn}}(x,u)-f_{\text{nn}}(x^{*},u^{*})$.

The steady-state feedback controller $u_{\text{s}}(t)$ is designed as

$u_{\text{s}}(t)=K\delta x(t)+u^{*},$ (9)

where $\delta x(t)=x(t)-x^{*}$ and $K\in\mathbb{R}^{m\times n}$ is a constant
gain. In this paper $K$ is designed such that $A_{s}=A+BK$ is Schur stable by
using the linear quadratic regulator (LQR) method (Anderson & Moore, 2007)
with the given state and input weights $Q\in\mathbb{R}^{n\times n}\succeq 0$
and $R\in\mathbb{R}^{m\times m}\succ 0$.

Substituting the controller (5) into (8) yields

$\delta x(t+1)=A_{s}\delta x(t)+Bu_{\text{c}}(t)+D\delta f_{\text{nn}}(t),$
(10)

where $u_{\text{c}}(t)$ is designed below through the MPC framework.

MPC controller design. The MPC design needs the terminal state set
$\mathcal{X}_{f}$ to ensure recursive feasibility. For the system (1), the set
$\mathcal{X}_{f}$ can be constructed based on (10) by setting
$u_{\text{c}}(t)=0$, as constraints are inactive at steady state. Since the
system matrix $A_{s}$ is Schur stable, the terminal set is defined as a robust
positively invariant (RPI) set for (10).

Let $\mathcal{X}_{\delta}$ be a RPI set for the error system (10) with the
disturbance term $D\delta f_{\text{nn}}(t)$. Given the constraints in (3) and
the steady-state target $(x^{*},u^{*})$, the intervals of $\delta x(t)$ and
$\delta f_{\text{nn}}(t)$ can be derived as
$[\underline{x}-x^{*},\bar{x}-x^{*}]$ and
$[\underline{f}_{\text{nn}}-f_{\text{nn}}(x^{*},u^{*}),\overline{f}_{\text{nn}}-f_{\text{nn}}(x^{*},u^{*})]$,
where $\underline{f}_{\text{nn}}$ and $\overline{f}_{\text{nn}}$ are computed
in (2). Based on the bounds of $\delta x(t)$ and $\delta f_{\text{nn}}(t)$,
the RPI set $\mathcal{X}_{\delta}$ is constructed using (Kolmanovsky &
Gilbert, 1998, Algorithm 6.1) and given as $\mathcal{X}_{\delta}=\\{\delta
x(t)\in\mathbb{R}^{n}\mid H\delta x(t)\leq h\\}$. Hence, the terminal set
$\mathcal{X}_{f}$ for the system (1) with the proposed controller (5) is
constructed as

$\mathcal{X}_{f}=\\{x(t)\in\mathbb{R}^{n}\mid Hx(t)\leq h+Hx^{*}\\}.$ (11)

The MPC policy at time $t$ is defined as $u_{\mathrm{mpc}}(t)=c^{*}(1|t)$
which is the first element of the optimal policy sequence
$\\{c(k|t)\\}_{k=1}^{N}$ solved from

$\displaystyle\hskip
20.00003pt\underset{\\{c(k|t)\\}_{k=1}^{N},\\{x(k|t)\\}_{k=1}^{N+1}}{\min}V_{t}$
$\displaystyle\mathrm{s.t.}~{}x(k+1|t)=Ax(k|t)+Bu(k|t)+Df_{\text{nn}}(k|t),\\!$
(12a) $\displaystyle\hskip 10.00002ptu(k|t)=K(x(k|t)-x^{*})+u^{*}+c(k|t),$
(12b) $\displaystyle\hskip
10.00002ptu(k|t)\in\mathcal{U},~{}x(1|t)=x(t),~{}x(k|t)\in\mathcal{X},$ (12c)
$\displaystyle\hskip 10.00002ptx(N+1|t)\in\mathcal{X}_{f},$ (12d)
$\displaystyle\hskip
10.00002ptz_{0}(k)=[x(k|t)^{\top},~{}u(k|t)^{\top}]^{\top},$ (12e)
$\displaystyle\hskip 10.00002ptz_{i}(k)=\phi(\hat{z}_{i}(k)),~{}i\in[1,L-1],$
(12f) $\displaystyle\hskip 10.00002ptf_{\text{nn}}(k|t)=\hat{z}_{L-1}(k),$
(12g) $\displaystyle\hskip 10.00002pt\underline{f}_{\text{nn}}\leq
f_{\text{nn}}(k|t)\leq\overline{f}_{\text{nn}},~{}k\in[1,N],$ (12h)

where the cost function is designed as
$V_{t}=\sum_{k=1}^{N}(\|x(k|t)-x^{*}\|^{2}_{Q}+\|u(k|t)-u^{*}\|^{2}_{R})+\|x(N+1|t)-x^{*}\|^{2}_{P}.$
$N$ is the prediction horizon. $Q$ and $R$ are those weight matrices used in
the LQR design for $K$. $P$ is the solution to the Riccati equation
$A^{\top}PA-P-(A^{\top}PB)(R+B^{\top}PB)^{-1}(B^{\top}PA)+Q=0$.
$\mathcal{X}_{f}$ is the terminal constraint set to ensure recursive
feasibility. $\Phi$ is a user-given weight matrix.

By implementing the proposed controller (9), stability of the tracking error
system under the constraints can be proved using the standard Lyapunov
function methods as in Rossiter (2003); Rawlings et al. (2017) and the details
are omitted here.

Due to the high-order nonlinear constraints in (12f), the optimisation problem
(12) is generally hard to solve. Several methods will be proposed in the next
section to efficiently solve it. Due to the use of concrete bound propagation,
the interval arithmetic method in (2) generates a conservative outer-
approximation of $f_{\text{nn}}(x(t),u(t))$, resulting in a small terminal
constraint set. Computing a tighter interval for $f_{\text{nn}}(x(t),u(t))$
would be useful to enlarge $\mathcal{X}_{f}$ and improve recursive feasibility
of the MPC problem. A method to achieve this will be presented later in
Section 4.

##  4 Methods for solving the MPC problem

This section describes methods for solving the formulated nonlinear MPC
problem (12), including an exact MIP method and two linear relaxation (LR)
methods.

###  4.1 MIP method

By definition, the $(i,j)$-neuron (i.e., the $j$-th ReLU neuron at the $i$-th
layer) can be either active (when $z_{i,j}=\hat{z}_{i,j}$) or inactive (when
$z_{i,j}=0$). An extra binary variable $\delta_{i,j}$ can be introduced to
indicate the neuron state: When $\delta_{i,j}=1$, the neuron is active; When
$\delta_{i,j}=0$, the neuron is inactive. Hence, each ReLU function
$z_{i,j}=\phi(\hat{z}_{i,j})$ is split into two cases (as illustrated in Fig.
1(b)):

$\displaystyle\delta_{i,j}=1\implies
z_{i,j}=\hat{z}_{i,j},~{}\delta_{i,j}=0\implies z_{i,j}=0.$ (13)

By using (13), the nonlinear ReLU constraints in (12f) can be represented as a
set of linear constraints as follows:

$\displaystyle z_{i}(k)\geq
0,~{}z_{i}(k)\geq\hat{z}_{i}(k),~{}z_{i}(k)\leq\hat{z}_{i}(k)-\hat{l}_{i}\odot(1-\delta_{k,i}),$
$\displaystyle
z_{i}(k)\leq\hat{u}_{i}\odot\delta_{k,i},~{}\delta_{k,i}\in\\{0,1\\}^{n_{i}},~{}i\in[1,L-1],$
(14)

where the bounds $\hat{l}_{i}$ and $\hat{u}_{i}$ are computed in (2) and
$\odot$ is the element-wise product. It can be easily shown that (4.1) is
always equivalent to (13) providing that
$\hat{l}_{i}\leq\hat{z}_{i}(k)\leq\hat{u}_{i}$, while conservativeness of the
bounds $\hat{l}_{i}$ and $\hat{u}_{i}$ is not a matter.

By applying (4.1), the MPC problem (12) is equivalently reformulated as the
following MIP problem:

$\displaystyle\underset{\\{c(k|t)\\}_{k=1}^{N},\\{x(k|t)\\}_{k=1}^{N+1},\\{\delta_{k,i}\\}_{k=1}^{N}}{\min}V_{t}$
$\displaystyle\\!\\!\mathrm{s.t.}~{}$ $\displaystyle\eqref{eq:origin mpc
const1}-\eqref{eq:origin mpc const4},\eqref{eq:origin mpc
const6},\eqref{eq:origin mpc const7},\eqref{MIP relaxation},k\in[1,N].$ (15a)

In the literature, a similar MIP encoding of ReLU neurons is used to design a
safe controller for NN-based dynamic systems (Wei & Liu, 2022) and analyse the
stability and feasibility of NN control systems (Karg & Lucia, 2020). The MIP
reformulation does not introduce any conservatism. But its computational
complexity grows rapidly with the number of binary variables and the
prediction horizon, which limits its scalability to large NNs and its
applicability in real-time control for deep NN-modelled dynamic systems. This
motivates the proposal of more computationally efficient solving methods
below.

###  4.2 Linear relaxation (LR) method

The LR method consists in using a triangle relaxation to over-approximate the
$(i,j)$-neuron as follows (see an illustrative example in Fig. 1(c)):

$z_{i,j}\geq 0,~{}z_{i,j}\geq\hat{z}_{i,j},~{}z_{i,j}\leq
a_{i,j}(\hat{z}_{i,j}-\hat{l}_{i,j})+\phi(\hat{l}_{i,j}),$ (16)

where
$a_{i,j}=(\phi(\hat{u}_{i,j})-\phi(\hat{l}_{i,j}))/(\hat{u}_{i,j}-\hat{l}_{i,j})$
when $\hat{l}_{i,j}\neq\hat{u}_{i,j}$ and $a_{i,j}=0$ otherwise. Note that an
over-approximation is necessary to ensure that if the computed MPC policy
satisfies the state and input constraints, then it is still the case when
applying the policy to the original system model. This is, however, not
achievable by using the sampling-based methods (Williams et al., 2016; Askari
et al., 2022) or using approximation of the ReLU such as the polynomial
approximation (Kochdumper et al., 2022).

Replacing the ReLU constraints (12f) with (16) yields the LR method:

$\displaystyle\underset{\\{c(k|t)\\}_{k=1}^{N},\\{x(k|t)\\}_{k=1}^{N+1}}{\min}V_{t}$
$\displaystyle\\!\\!\mathrm{s.t.}~{}$ $\displaystyle\eqref{eq:origin mpc
const1}-\eqref{eq:origin mpc const4},\eqref{eq:origin mpc
const6},\eqref{eq:origin mpc const7},$ (17a) $\displaystyle\eqref{LR
relaxation},j\in[1,n_{i}],i\in[1,L-1],k\in[1,N-1].$ (17b)

The problem (17) is relatively easier to solve than the MIP problem (15), but
there is usually a large relaxation gap, especially for large-scale NNs, which
leads to undesirable output tracking performance. A method is further
developed below to improve the solution quality of the LR method.

###  4.3 Enhanced linear relaxation (eLR) method

Due to the use of over-approximation of ReLU in the LR method, the computed
MPC cannot ensure the system reach the target steady state, leading to
inaccurate output tracking. To improve the tracking accuracy, the cost
function is redesigned as

$\displaystyle
J_{t}=\sum_{k=1}^{N}(\|x(k|t)-x^{*}\|^{2}_{Q}+\|u(k|t)-u^{*}\|^{2}_{R}$
$\displaystyle+\|f_{\text{nn}}(k|t)-f_{\text{nn}}(x^{*},u^{*})\|^{2}_{\Phi})+\|x(N+1|t)-x^{*}\|^{2}_{P},$
(18)

where the term $\|f_{\text{nn}}(k|t)-f_{\text{nn}}(x^{*},u^{*})\|^{2}_{\Phi}$
with a user-given weight matrix $\Phi$ is introduced to minimise the deviation
of the over-approximated NN output $f_{\text{nn}}(k|t)$ from its steady-state
value $f_{\text{nn}}(x^{*},u^{*})$. By minimising this new cost function, the
obtained MPC policy can minimise the total deviation of the system from its
steady-state dynamics in (6). Due to introduction of the new term, it holds
that $J_{t}\geq V_{t}$.

Applying the new cost function in (4.3) to (17) results in the enhanced LR
(eLR) method:

$\displaystyle\underset{\\{c(k|t)\\}_{k=1}^{N},\\{x(k|t)\\}_{k=1}^{N+1}}{\min}J_{t}$
$\displaystyle\\!\\!\mathrm{s.t.}~{}$ $\displaystyle\eqref{eq:origin mpc
const1}-\eqref{eq:origin mpc const4},~{}\eqref{eq:origin mpc
const6},~{}\eqref{eq:origin mpc const7},$ (19a) $\displaystyle\eqref{LR
relaxation},j\in[1,n_{i}],i\in[1,L-1],k\in[1,N-1].$ (19b)

Note that when the optimal solution to (19) satisfies $J_{t}^{*}=0$, the
obtained MPC ensures that the system reaches the target steady-state and
accurate output tracking.

###  4.4 Discussions

Proposition 4.1 provides a comparison of the solution quality (optimality or
tracking accuracy) $S_{\text{method}}$ of the MIP, LR, and eLR methods against
the original MPC (12).

######  Proposition 4.1.

The solution quality of the proposed methods satisfy: $S_{\text{LR}}\leq
S_{\text{eLR}}\leq S_{\text{MIP}}=S_{\text{origin}}.$

The LR and enhanced LR methods are of the same size and will have roughly the
same computational complexity, while the MIP method is generally much more
computationally expensive. A further demonstration of this is to be through
simulation in Section 5.

Computational complexity of the proposed solving methods can be reduced via
directly encoding the stable (either inactive or active) neurons using
equality constraints. For each $(i,j)$-neuron, the pre-activation bounds
$\hat{l}_{i,j}$ and $\hat{u}_{i,j}$ can be used to identity its status as
follows: 1) If $\hat{u}_{i,j}\leq 0$, the neuron is known to be strictly
inactive and thus $z_{i,j}=0$, and 2) If $\hat{l}_{i,j}\geq 0$, the neuron is
known to be strictly active and thus $z_{i,j}=\hat{z}_{i,j}$. Applying these
to the proposed methods can reduce the required numbers of binary variables or
triangle relaxation. This can then improve the computational efficiency
without compromising the solution quality. It is also worth noting that more
stable neurons can be identified if less conservative pre-activation bounds
are used. The symbolic bound propagation method (Henriksen & Lomuscio, 2020)
can be an efficient way to achieve this and it will be explored in the future
work.

This paper focuses on developing efficient algorithm for solving the MPC
problem, so the NN $f_{\text{nn}}$ is assumed to capture the system nonlinear
dynamics and uncertainty. If a separate term is used to lump the total effects
of uncertainty and diturbance, the proposed MPC can be extended as robust MPC
based on the concept of tube-MPC (Langson et al., 2004), in which the proposed
solving algorithms will still be applicable.

##  5 Numerical example

Consider an inverted pendulum system

$\displaystyle x_{1}(t+1)$ $\displaystyle=x_{1}(t)+t_{s}x_{2}(t),$
$\displaystyle x_{2}(t+1)$
$\displaystyle=x_{2}(t)+\frac{t_{s}g\sin(x_{1}(t))}{l}-\frac{t_{s}cx_{2}(t)}{ml^{2}}+\frac{t_{s}u(t)}{ml^{2}},$
$\displaystyle y(t)$ $\displaystyle=x_{1}(t),$

where $x_{1}(t)$ is the angular displacement, $x_{2}(t)$ is the angular
velocity, $t_{s}$ is the sampling time, $m$ is the mass of pendulum, $l$ is
the distance from pivot to centre of mass of the pendulum, $c$ is the
rotational friction coefficient, and $g$ is the gravitational constant. The
parameters used in the simulation are $t_{s}=0.1~{}\mathrm{s}$, $m=1$, $l=1$,
$g=9.8$ and $c=0.01$. The input and state constraints are: $|u(t)|\leq
3,~{}|x(t)|\leq[\pi/2,~{}5]^{\top}$.

Suppose the value of the rotational friction coefficient $c$ is unknown for
the control design. The pendulum system is reformulated as

$\begin{split}x(t+1)&=Ax(t)+Bu(t)+Df(x(t)),\\\ y(t)&=Cx(t),\end{split}$ (20)

where $A=\begin{bmatrix}1&t_{s}\\\ 0&1\end{bmatrix}$, $B=\begin{bmatrix}0\\\
\frac{t_{s}}{ml^{2}}\end{bmatrix}$, $D=\begin{bmatrix}0\\\
\frac{t_{s}}{l}\end{bmatrix}$, $C=[1~{}0]$ and
$f(x(t))=g\sin(x_{1}(t))-(cx_{2}(t))/(ml).$

A three-input one-output NN $f_{\text{nn}}(x(t),u(t))$ is trained to model the
unknown nonlinear dynamics $f(x(t))$. Then the pendulum system can be
represented in the form of (1). The simulations are conducted in Matlab on a
Linux machine with an AMD Ryzen Threadripper PRO 5955WX 16-cores CPU and 128
GB RAM. The optimisation problems are modelled using YALMIP (Löfberg, 2004)
and solved by Gurobi (Gurobi Optimization, LLC, 2023). The weight matrices
used in the control design are: $R_{s}=1$,
$Q=[1.0\mathrm{e}5,0;0,1.0\mathrm{e}2]$, $R=1$, and $\Phi=100$. The initial
state is $x(0)=[0;0]$ and the output reference is $y_{r}=\pi/5$.

Figure 2: Tracking performance and runtime for 3-layer NNs with different
number of neurons at the hidden layer. Figure 3: Example output tracking and
control input for the 3-layer NN with 50 neurons in Case 1. Figure 4: Tracking
performance and runtime for NNs with different number of hidden layers each
having 10 neurons. Figure 5: Tracking performance and runtime for 3-layer NNs
with 50 neurons under different prediction horizons.

Comparison is made for the MIP method in (15), the LR method in (17), and the
eLR method in (19). The performance metrics to compare are the steady-state
relative tracking error $|y-y_{r}|/y_{r}\times 100$ and the maximum runtime
for solving the MPC problems at each sampling step. Three simulation cases are
used to make a comprehensive comparison:

  * •

Case 1: Network width test. This case simulates 3-layer NNs where the number
of neurons at the only hidden layer is 10, 20, 50, 100, or 200, repectively.
The prediction horizon is $N=1$. The results are depicted in Fig. 2 and an
example of the outputs and control inputs for the 50 neurons NN is provided in
Fig. 3.

  * •

Case 2: Network depth test. This case uses NNs having 1, 2, 4 or 5 hidden
layers with 10 neurons at each hidden layer. The prediction horizon is $N=1$.
The results are reported in Fig. 4.

  * •

Case 3: Prediction horizon test. This case considers a 3-layer NN with 50
neurons. The prediction horizon is $N=2,3,4,5$, respectively. The results are
reported in Fig. 5.

In all the three cases, the tracking accuracy from the best to the worst are:
MIP, eLR, and LR. Moreover, the relative tracking errors of the LR method grow
rapidly with increase in the network width (i.e., the number of neurons) and
depth (i.e., the number of hidden layers). But the tracking errors of the MIP
and eLR methods are almost zero in all simulations.

As expected, the relaxation methods, LR and eLR methods, need around the same
runtime in all the simulations, while the MIP method is much more
computationally expensive. Moreover, the runtime of MIP grows much faster than
the other methods when the NNs become larger (with more neurons or hidden
layers) or the prediction horizon becomes longer. The runtime of MIP method
exceeds the sampling time 0.1 s when the 3-layer NNs have more than 50 neurons
per layer (as shown in Fig. 2), when the NNs have more than 2 hidden layers
(as shown in Fig. 4), or when the prediction horizon is longer than 2 (as
shown in Fig. 5). However, the runtime for the relaxation methods remains
lower than the sampling time across all the NNs and simulations.

##  6 Conclusion

A dual-mode MPC is designed to ensure output set point tracking of NN-modelled
dynamic systems in the presence of state and input constraints. Both exact MIP
method and linear relaxation (LR) methods are proposed to solve the formulated
nonlinear MPC problem. Theoretic analysis and simulation evaluation show that
the enhanced LR can enjoy both the benefits of low computational burden from
the LR method and accurate output tracking from the MIP method. This suggests
the enhanced LR method is more suitable for real-time control of deep NN-
modelled dynamic systems. Future work will be extending the proposed control
design for NNs with other types of activation functions.

## References

  * Anderson & Moore (2007) Anderson, B. D., & Moore, J. B. (2007).  Optimal control: linear quadratic methods.  Courier Corporation.
  * Askari et al. (2022) Askari, I. et al. (2022).  Sampling-based nonlinear MPC of neural network dynamics with application to autonomous vehicle motion planning.  In Proc. ACC (pp. 2084–2090).  IEEE.
  * Bonassi et al. (2024) Bonassi, F., La Bella, A., Farina, M., & Scattolini, R. (2024).  Nonlinear mpc design for incrementally iss systems with application to gru networks.  Automatica, 159, 111381.
  * D’Amico et al. (2023) D’Amico, W., La Bella, A., & Farina, M. (2023).  An incremental input-to-state stability condition for a class of recurrent neural networks.  IEEE Trans. Autom. Control, 69, 2221–2236.
  * D’Errico (2023) D’Errico, J. (2023).  fminsearchbnd, fminsearchcon.  URL: https://www.mathworks.com/matlab central/fileexchange/8277-fminsearchbnd-fmin searchcon accessed: 2023-02-28.
  * Gurobi Optimization, LLC (2023) Gurobi Optimization, LLC (2023).  Gurobi Optimizer Reference Manual.  URL: https://www.gurobi.com.
  * Henriksen & Lomuscio (2020) Henriksen, P., & Lomuscio, A. (2020).  Efficient neural network verification via adaptive refinement and adversarial search.  In Proc. ECAI20 (pp. 2513–2520).  IOS Press.
  * Karg & Lucia (2020) Karg, B., & Lucia, S. (2020).  Stability and feasibility of neural network-based controllers via output range analysis.  In Proc. IEEE Conf. Decis. Control (pp. 4947–4954).  IEEE.
  * Kochdumper et al. (2022) Kochdumper, N., Schilling, C., Althoff, M., & Bak, S. (2022).  Open-and closed-loop neural network verification using polynomial zonotopes.  arXiv preprint arXiv:2207.02715, .
  * Kolmanovsky & Gilbert (1998) Kolmanovsky, I., & Gilbert, E. G. (1998).  Theory and computation of disturbance invariant sets for discrete-time linear systems.  Mathematical Problems in Engineering, 4, 317–367.
  * Langson et al. (2004) Langson, W., Chryssochoos, I., Raković, S., & Mayne, D. Q. (2004).  Robust model predictive control using tubes.  Automatica, 40, 125–133.
  * Ławryńczuk (2010) Ławryńczuk, M. (2010).  Suboptimal nonlinear predictive control based on multivariable neural hammerstein models.  Appl. Intell., 32, 173–192.
  * Löfberg (2004) Löfberg, J. (2004).  YALMIP: A toolbox for modeling and optimization in MATLAB.  In Proc. CACSD.  volume 3.
  * Moe et al. (2018) Moe, S., Rustad, A. M., & Hanssen, K. G. (2018).  Machine learning in control systems: An overview of the state of the art.  In Proc. Innov. Appl. Artif. Intell. Conf. (pp. 250–265).  Springer.
  * Nagariya & Saripalli (2020) Nagariya, A., & Saripalli, S. (2020).  An iterative LQR controller for off-road and on-road vehicles using a neural network dynamics model.  In Proc. IEEE IV (pp. 1740–1745).  IEEE.
  * Nahas et al. (1992) Nahas, E., Henson, M., & Seborg, D. (1992).  Nonlinear internal model control strategy for neural network models.  Comput. Chem. Eng., 16, 1039–1057.
  * Nikolakopoulou et al. (2022) Nikolakopoulou, A., Hong, M. S., & Braatz, R. D. (2022).  Dynamic state feedback controller and observer design for dynamic artificial neural network models.  Automatica, 146, 110622.
  * Patan (2018) Patan, K. (2018).  Two stage neural network modelling for robust model predictive control.  ISA Trans., 72, 56–65.
  * Rawlings et al. (2017) Rawlings, J. B., Mayne, D. Q., & Diehl, M. (2017).  Model predictive control: theory, computation, and design volume 2.  Nob Hill Publishing Madison, WI.
  * Ren et al. (2022) Ren, Y. M. et al. (2022).  A tutorial review of neural network modeling approaches for model predictive control.  Comput. Chem. Eng., (p. 107956).
  * Rokonuzzaman et al. (2021) Rokonuzzaman, M., Mohajer, N., Nahavandi, S., & Mohamed, S. (2021).  Model predictive control with learned vehicle dynamics for autonomous vehicle path tracking.  IEEE Access, 9, 128233–128249.
  * Rossiter (2003) Rossiter, J. A. (2003).  Model-based predictive control: a practical approach.  CRC press.
  * Saint-Donat et al. (1991) Saint-Donat, J., Bhat, N., & McAvoy, T. J. (1991).  Neural net based model predictive control.  Int. J. Control, 54, 1453–1468.
  * Salzmann et al. (2023) Salzmann, T. et al. (2023).  Real-time neural MPC: Deep learning model predictive control for quadrotors and agile robotic platforms.  IEEE Robot. Autom. Lett., 8, 2397–2404.
  * Sørensen et al. (1999) Sørensen, P. H., Nørgaard, M., Ravn, O., & Poulsen, N. K. (1999).  Implementation of neural network based non-linear predictive control.  Neurocomputing, 28, 37–51.
  * Spielberg et al. (2021) Spielberg, N. A., Brown, M., & Gerdes, J. C. (2021).  Neural network model predictive motion control applied to automated driving with unknown friction.  IEEE Trans. Control Syst. Technol., 30, 1934–1945.
  * Spielberg et al. (2019) Spielberg, N. A. et al. (2019).  Neural network vehicle models for high-performance automated driving.  Sci. Robot., 4, eaaw1975.
  * Tang et al. (2022) Tang, Y. et al. (2022).  Perception and navigation in autonomous systems in the era of learning: A survey.  IEEE Trans. Neural Netw. Learn. Syst., .
  * Wei & Liu (2022) Wei, T., & Liu, C. (2022).  Safe control with neural network dynamic models.  In Learning for Dynamics and Control Conference (pp. 739–750).  PMLR.
  * Williams et al. (2016) Williams, G. et al. (2016).  Information theoretic MPC using neural network dynamics.  In Proc. NeurIPS.
  * Wysocki & Ławryńczuk (2015) Wysocki, A., & Ławryńczuk, M. (2015).  Jordan neural network for modelling and predictive control of dynamic systems.  In Proc. MMAR (pp. 145–150).  IEEE.
  * Yan & Wang (2012) Yan, Z., & Wang, J. (2012).  Model predictive control of nonlinear systems with unmodeled dynamics based on feedforward and recurrent neural networks.  IEEE Trans. Industr. Inform., 8, 746–756.


######  Corollary 4.14 (Multi-resource without Regret).

Let $\mathcal{I}=\\{(A_{i},R_{i})_{i=1}^{m}\\}$ be an $m$-player multi-
resource allocation game where $A_{i}=P_{k_{i}}(n_{i,1})\times\cdots\times
P_{k_{i}}(n_{i,B})$ for $n_{1},\cdots n_{B},k\in\mathbb{Z}^{+}$. Suppose the
reward of the first player is bounded by $L_{\max}$,
$\max_{i\in[m],b\in[B]}n_{i,b}=n$, $\max_{i\in[m]}k_{i}=k$. Then there is a
no-regret learning algorithm for $\mathcal{I}$ with regret:

$\operatorname{Reg}_{T}\leq
O\left(L_{\max}\sqrt{T}\cdot\left(\sqrt{\sum_{b=1}^{B}\log(n_{B})}+\sqrt{\log(1/\eta)}\right)\right)$

with probability at least $1-\eta$ that runs in time

$O\left(TkB(n+1)^{2B}+mTk^{2}B\cdot(n+1)^{B}\right)$

assuming query access to
$r_{i,h}:([0,n_{i,1}])\times\cdots\times([0,n_{i,B}])\times
A_{-i}\mapsto\mathbb{R}$.

######  Corollary 4.15 (Equilibrium Computation for Multi-Resource Allocation
Games).

Let $\mathcal{I}=\\{(A_{i},R_{i})_{i=1}^{m}\\}$ be an $m$-player multi-
resource allocation game where $A_{i}=P_{k_{i}}(n_{i,1})\times\cdots\times
P_{k_{i}}(n_{i,B})$ for $n_{1},\cdots n_{B},k\in\mathbb{Z}^{+}$. Suppose the
reward of the first player is bounded by $L_{\max}$,
$\max_{i\in[m],b\in[B]}n_{i,b}=n$, $\max_{i\in[m]}k_{i}=k$. There exists an
algorithm to compute an $\varepsilon$-approximate CCE (Nash if two-player
zero-sum) with probability at least $1-\eta$ in time

$O\left(m(n+1)^{2B}k^{2}B^{2}L_{\max}^{2}\varepsilon^{-2}\log(mn/\eta)+m^{2}(n+1)^{B}k^{3}B^{2}L_{\max}^{2}\varepsilon^{-2}\cdot\log(mn/\eta)\right)$

assuming query access to
$r_{i,h}:([0,n_{i,1}])\times\cdots\times([0,n_{i,B}])\times
A_{-i}\mapsto\mathbb{R}$.

## Acknowledgements

We thank Thuy Duong Vuong, Nima Anari, and Kuikui Liu for helpful discussions
on MCMC-based approaches for Dueling games and Colonel Blotto, and Saeed
Seddighin for discussions on algorithms for Colonel Blotto and related
problems. We also thank Maxwell Fishelson for enlightening discussion of
optimistic hedge and its variants.

## References

  * AASV [21] Yeganeh Alimohammadi, Nima Anari, Kirankumar Shiragur, and Thuy-Duong Vuong, _Fractionally log-concave and sector-stable polynomials: counting planar matchings and more_ , Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, 2021, pp. 433–446.
  * AB [21] Alexandr Andoni and Daniel Beaglehole, _Learning to hash robustly, guaranteed_ , arXiv preprint arXiv:2108.05433 (2021).
  * ABL [14] Jean-Yves Audibert, Sébastien Bubeck, and Gábor Lugosi, _Regret in online combinatorial optimization_ , Mathematics of Operations Research 39 (2014), no. 1, 31–45.
  * ADF+ [22] Ioannis Anagnostides, Constantinos Daskalakis, Gabriele Farina, Maxwell Fishelson, Noah Golowich, and Tuomas Sandholm, _Near-optimal no-regret learning for correlated equilibria in multi-player general-sum games_ , Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, 2022, pp. 736–749.
  * ADH+ [19] Amirmahdi Ahmadinejad, Sina Dehghani, Mohammadtaghi Hajiaghayi, Brendan Lucier, Hamid Mahini, and Saeed Seddighin, _From duels to battlefields: Computing equilibria of blotto and other games_ , Mathematics of Operations Research 44 (2019), no. 4, 1304–1325.
  * AJK+ [21] Nima Anari, Vishesh Jain, Frederic Koehler, Huy Tuan Pham, and Thuy-Duong Vuong, _Entropic independence i: modified log-sobolev inequalities for fractionally log-concave distributions and high-temperature ising models_ , arXiv preprint arXiv:2106.04105 (2021).
  * AL [20] Vedat Levi Alev and Lap Chi Lau, _Improved analysis of higher order random walks and applications_ , arXiv preprint arXiv:2001.02827 (2020).
  * ALO [20] Nima Anari, Kuikui Liu, and Shayan Oveis Gharan, _Spectral independence in high-dimensional expanders and applications to the hardcore model_ , arXiv preprint arXiv:2001.00303 (2020).
  * ALOV [19] Nima Anari, Kuikui Liu, Shayan Oveis Gharan, and Cynthia Vinzant, _Log-concave polynomials ii: high-dimensional walks and an fpras for counting bases of a matroid_ , Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, 2019, pp. 1–12.
  * Art [94] W Brian Arthur, _Inductive reasoning and bounded rationality_ , The American economic review 84 (1994), no. 2, 406–411.
  * ARV [08] Heiner Ackermann, Heiko Röglin, and Berthold Vöcking, _On the impact of combinatorial structure on congestion games_ , Journal of the ACM (JACM) 55 (2008), no. 6, 1–22.
  * Aum [74] Robert J. Aumann, _Subjectivity and correlation in randomized strategies_ , Journal of Mathematical Economics 1 (1974), no. 1, 67–96.
  * BAEJ [20] Enric Boix-Adserà, Benjamin L. Edelman, and Siddhartha Jayanti, _The multiplayer colonel blotto game_ , Proceedings of the 21st ACM Conference on Economics and Computation (2020).
  * BB [19] Mourad Baïou and Francisco Barahona, _Faster algorithms for security games on matroids_ , Algorithmica 81 (2019), no. 3, 1232–1246.
  * BBD+ [19] Soheil Behnezhad, Avrim Blum, Mahsa Derakhshan, MohammadTaghi Hajiaghayi, Christos H. Papadimitriou, and Saeed Seddighin, _Optimal strategies of blotto games: Beyond convexity_ , 2019.
  * BDD+ [17] Soheil Behnezhad, Sina Dehghani, Mahsa Derakhshan, MohammadTaghi HajiAghayi, and Saeed Seddighin, _Faster and simpler algorithm for optimal strategies of blotto game_ , Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence (2017), 369–375.
  * BHVW [11] Petter Brändén, James Haglund, Mirkó Visontai, and David G Wagner, _Proof of the monotone column permanent conjecture_ , Notions of Positivity and the Geometry of Polynomials, Springer, 2011, pp. 63–78.
  * Bor [53] Emile Borel, _The theory of play and integral equations with skew symmetric kernels_ , Econometrica 21 (1953), no. 1, 97.
  * BR [21] Yakov Babichenko and Aviad Rubinstein, _Settling the complexity of nash equilibrium in congestion games_ , Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, 2021, pp. 1426–1437.
  * Bri [88] E Oran Brigham, _The fast fourier transform and its applications_ , Prentice-Hall, Inc., 1988.
  * BŠVV [04] Ivona Bezáková, Daniel Štefankovic, Vijay V Vazirani, and Eric Vigoda, _Approximating the permanent in o (n7) time_.
  * CBL [06] Nicolo Cesa-Bianchi and Gábor Lugosi, _Prediction, learning, and games_ , Cambridge university press, 2006.
  * CBL [12] , _Combinatorial bandits_ , Journal of Computer and System Sciences 78 (2012), no. 5, 1404–1422.
  * CGM [19] Mary Cryan, Heng Guo, and Giorgos Mousa, _Modified log-sobolev inequalities for strongly log-concave distributions_ , 2019 IEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), IEEE, 2019, pp. 1358–1370.
  * CHK [17] Alon Cohen, Tamir Hazan, and Tomer Koren, _Tight bounds for bandit combinatorial optimization_ , Conference on Learning Theory, PMLR, 2017, pp. 629–642.
  * CLV [21] Zongchen Chen, Kuikui Liu, and Eric Vigoda, _Optimal mixing of glauber dynamics: Entropy factorization via high-dimensional expansion_ , Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, 2021, pp. 1537–1550.
  * CPV [95] Charles J Colbourn, J Scott Provan, and Dirk Vertigan, _The complexity of computing the tutte polynomial on transversal matroids_ , Combinatorica 15 (1995), no. 1, 1–10.
  * DFG [21] Constantinos Daskalakis, Maxwell Fishelson, and Noah Golowich, _Near-optimal no-regret learning in general games_ , Advances in Neural Information Processing Systems 34 (2021), 27604–27616.
  * DGP [06] Constantinos Daskalakis, Paul W. Goldberg, and Christos H. Papadimitriou, _The complexity of computing a nash equilibrium_ , Proceedings of the thirty-eighth annual ACM symposium on Theory of computing - STOC 06 (2006).
  * DGP [09] Constantinos Daskalakis, Paul W Goldberg, and Christos H Papadimitriou, _The complexity of computing a nash equilibrium_ , SIAM Journal on Computing 39 (2009), no. 1, 195–259.
  * DJM [17] Martin Dyer, Mark Jerrum, and Haiko Müller, _On the switch markov chain for perfect matchings_ , Journal of the ACM (JACM) 64 (2017), no. 2, 1–33.
  * DSDMDB [06] Bart De Schuymer, Hans De Meyer, and Bernard De Baets, _Optimal strategies for equal-sum dice games_ , Discrete applied mathematics 154 (2006), no. 18, 2565–2576.
  * Eck [87] Roger Eckhardt, _Stan ulam, john von neumann, and the monte carlo method_ , Los Alamos Science 15 (1987), no. 131-136, 30.
  * FS [96] Yoav Freund and Robert E. Schapire, _Game theory, on-line prediction and boosting_ , Proceedings of the ninth annual conference on Computational learning theory - COLT ’96 (1996).
  * FS [99] , _Adaptive game playing using multiplicative weights_ , Games and Economic Behavior 29 (1999), no. 1-2, 79–103.
  * GP [09] Russell Golman and Scott E. Page, _General blotto: Games of allocative strategic mismatch_ , Public Choice 138 (2009), no. 3/4, 279–299.
  * Gro [50] Oliver Gross, _The symmetric blotto game [1]_ , Rand Corporation, 1950.
  * GW [50] Oliver Gross and Robert Wagner, _A continuous colonel blotto game_ , Tech. report, Rand Project Air Force Santa Monica Ca, 1950.
  * Har [08] Sergiu Hart, _Discrete colonel blotto and general lotto games_ , International Journal of Game Theory 36 (2008), no. 3, 441–460.
  * HKW [09] David P Helmbold and Manfred K Warmuth, _Learning permutations with exponential weights._ , Journal of Machine Learning Research 10 (2009), no. 7.
  * IKL+ [11] Nicole Immorlica, Adam Tauman Kalai, Brendan Lucier, Ankur Moitra, Andrew Postlewaite, and Moshe Tennenholtz, _Dueling algorithms_ , Proceedings of the forty-third annual ACM symposium on Theory of computing, 2011, pp. 215–224.
  * Ins [16] Insider, _scientists say they can predict two-party outcomes after solving a 95-year-old game theory problem._ , 2016.
  * JSV [04] Mark Jerrum, Alistair Sinclair, and Eric Vigoda, _A polynomial-time approximation algorithm for the permanent of a matrix with nonnegative entries_ , Journal of the ACM (JACM) 51 (2004), no. 4, 671–697.
  * JVV [86] Mark R Jerrum, Leslie G Valiant, and Vijay V Vazirani, _Random generation of combinatorial structures from a uniform distribution_ , Theoretical computer science 43 (1986), 169–188.
  * Kle [05] Robert David Kleinberg, _Online decision problems with large strategy sets_ , Ph.D. thesis, Massachusetts Institute of Technology, 2005.
  * KO [20] Tali Kaufman and Izhar Oppenheim, _High order random walks: Beyond spectral gap_ , Combinatorica (2020), 1–37.
  * KPT [09] Robert Kleinberg, Georgios Piliouras, and Éva Tardos, _Multiplicative updates outperform generic no-regret learning in congestion games_ , Proceedings of the forty-first annual ACM symposium on Theory of computing, 2009, pp. 533–542.
  * KR [12] Dan Kovenock and Brian Roberson, _Coalitional colonel blotto games with application to the economics of alliances_ , Journal of Public Economic Theory 14 (2012), no. 4, 653–676.
  * KR [21] , _Generalizations of the general lotto and colonel blotto games_ , Economic Theory 71 (2021), no. 3, 997–1032.
  * Kva [07] Dmitriy Kvasov, _Contests with limited resources_ , Journal of Economic Theory 136 (2007), no. 1, 738–748.
  * KWK+ [10] Wouter M Koolen, Manfred K Warmuth, Jyrki Kivinen, et al., _Hedging structured concepts._ , COLT, Citeseer, 2010, pp. 93–105.
  * Las [02] Jean-François Laslier, _How two-party competition treats minorities_ , Review of Economic Design 7 (2002), no. 3, 297–307.
  * LP [02] Jean-François Laslier and Nathalie Picard, _Distributive politics and electoral competition_ , Journal of Economic Theory 103 (2002), no. 1, 106–130.
  * LW [94] Nick Littlestone and Manfred K Warmuth, _The weighted majority algorithm_ , Information and computation 108 (1994), no. 2, 212–261.
  * LY [94] Richard J Lipton and Neal E Young, _Simple strategies for large zero-sum games with applications to complexity theory_ , Proceedings of the twenty-sixth annual ACM symposium on Theory of computing, 1994, pp. 734–740.
  * Mag [05] Avner Magen, _Lecture notes from csc2411 spring 2005_ , 2005.
  * Mew [16] Mewright, _Umd-led team first to solve well-known game theory scenario_ , Feb 2016.
  * MM [15] Scott T Macdonell and Nick Mastronardi, _Waging simple wars: a complete characterization of two-battlefield blotto equilibria_ , Economic Theory 58 (2015), no. 1, 183–216.
  * Nas [51] John Nash, _Non-cooperative games_ , Annals of mathematics (1951), 286–295.
  * OvS [13] Georg Ostrovski and Sebastian van Strien, _Payoff performance of fictitious play_ , arXiv preprint arXiv:1308.4049 (2013).
  * PRG [22] Vianney Perchet, Philippe Rigollet, and Thibaut Le Gouic, _An algorithmic solution to the blotto game using multi-marginal couplings_ , 2022.
  * Rob [06] Brian Roberson, _The colonel blotto game_ , Economic Theory 29 (2006), no. 1, 1–24.
  * Ros [73] Robert W Rosenthal, _A class of games possessing pure-strategy nash equilibria_ , International Journal of Game Theory 2 (1973), no. 1, 65–67.
  * RS [13] Sasha Rakhlin and Karthik Sridharan, _Optimization, learning, and games with predictable sequences_ , Advances in Neural Information Processing Systems 26 (2013).
  * Rub [16] Aviad Rubinstein, _Settling the complexity of computing approximate two-player nash equilibria_ , 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), IEEE, 2016, pp. 258–265.
  * Sch [96] James C Schatzman, _Accuracy of the discrete fourier transform and the fast fourier transform_ , SIAM Journal on Scientific Computing 17 (1996), no. 5, 1150–1166.
  * SLS [14] Galina Schwartz, Patrick Loiseau, and Shankar S Sastry, _The heterogeneous colonel blotto game_ , 2014 7th international conference on NETwork Games, COntrol and OPtimization (NetGCoop), IEEE, 2014, pp. 232–238.
  * Sze [17] Dávid Szeszlér, _Security games on matroids_ , Mathematical Programming 161 (2017), no. 1, 347–364.
  * Tam [11] Milind Tambe, _Security and game theory: algorithms, deployed systems, lessons learned_ , Cambridge university press, 2011.
  * Tho [17] Caroline Thomas, _N-dimensional blotto game with heterogeneous battlefield values_ , Economic Theory 65 (2017), no. 3, 509–544.
  * Tho [18] , _N-dimensional blotto game with heterogeneous battlefield values_ , Economic Theory 65 (2018), no. 3, 509–544.
  * VA [06] Berthold Vöcking and R Aachen, _Congestion games: Optimization in competition._ , ACiD, 2006, pp. 9–20.
  * VL [21] Dong Quan Vu and Patrick Loiseau, _Colonel blotto games with favoritism: Competitions with pre-allocations and asymmetric effectiveness_ , 2021.
  * VLS [18] Dong Quan Vu, Patrick Loiseau, and Alonso Silva, _Efficient computation of approximate equilibria in discrete colonel blotto games_ , Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence (2018).
  * VLS [19] , _Combinatorial bandits for sequential learning in colonel blotto games_ , 2019 IEEE 58th Conference on Decision and Control (CDC), IEEE, 2019, pp. 867–872.
  * WK [08] Manfred K Warmuth and Dima Kuzmin, _Randomized online pca algorithms with regret bounds that are logarithmic in the dimension_ , Journal of Machine Learning Research 9 (2008), no. Oct, 2287–2320.

##  Appendix A No-regret Learning and Equilibrium Computation

######  Lemma A.1 (Lemma 4.1 from [22] Rephrased).

Let $\mathcal{I}=\\{(A_{i},R_{i})_{i=1}^{m}\\}$ be an $m$-player game played
repeatedly for $T$ rounds. Denote $\mathbf{s}^{(t)}_{i}$ as the mixed strategy
chosen by the $i$-th player, and $a^{(t)}_{i}\sim\mathbf{s}^{(t)}_{i}$ as the
action sampled. Assume the $i$-th player follows an algorithm which computes
$\mathbf{s}^{(t)}_{i}$ solely based on $a^{(t^{\prime})}_{j}$ for
$t^{\prime}<t$ and $j\in[m]\backslash\\{i\\}$. Furthermore, suppose the
following is true

$\sup_{b^{(t)}_{j}\in A_{j}\text{ for }t\in[T],j\neq
i}\mathbb{E}\left[\max_{e\in
A_{i}}\sum_{t=1}^{T}R_{i}(e,b^{(t)}_{-i})-\sum_{t=1}^{T}R_{i}(a^{(t)}_{i},b^{(t)}_{-i})\right]\leq
B.$

Then, for all $\delta\in(0,1)$, with probability at least $1-\delta$, it holds

$\max_{e\in
A_{i}}\sum_{t=1}^{T}R_{i}(e,a^{(t)}_{-i})-\sum_{t=1}^{T}R_{i}(a^{(t)}_{i},a^{(t)}_{-i})\leq
B+L_{\max}\sqrt{T/2\log(1/\delta)}.$

at the end of the repeated play for the $i$-th player.

###### Proof of Lemma 2.10.

If one follows exactly from the Randomized Weighted Majority algorithm, one
has the guarantee that

$\sup_{b^{(t)}_{j}\in A_{j}\text{ for }t\in[T],j\neq
i}\mathbb{E}_{a^{(t)}_{i}\sim RM^{(t)}_{i}\left(\beta\right)}\left[\max_{e\in
A_{i}}\sum_{t=1}^{T}R_{i}(e,b^{(t)}_{-i})-\sum_{t=1}^{T}R_{i}(a^{(t)}_{i},b^{(t)}_{-i})\right]\leq
L_{\max}\sqrt{T\log N}.$

In reality, since we are performing $\kappa$-approximate sampling, we have
$\operatorname{TV}\left(RM^{(t)}_{i}\left(\beta\right),\mathbf{s}^{(t)}_{i}\right)\leq\kappa$.
Since the reward of the game is bounded by $L_{\max}$, we have

$\sup_{b^{(t)}_{j}\in A_{j}\text{ for }t\in[T],j\neq
i}\mathbb{E}_{a^{(t)}_{i}\sim\mathbf{s}^{(t)}_{i}}\left[\max_{e\in
A_{i}}\sum_{t=1}^{T}R_{i}(e,b^{(t)}_{-i})-\sum_{t=1}^{T}R_{i}(a^{(t)}_{i},b^{(t)}_{-i})\right]\leq
L_{\max}\sqrt{T\log N}+\kappa\cdot L_{\max}\cdot T.$

Then, by Lemma A.1, it then holds

$\frac{1}{T}\max_{e\in
A_{i}}\sum_{t=1}^{T}R_{i}(e,a^{(t)}_{-i})-\frac{1}{T}\sum_{t=1}^{T}R_{i}(a^{(t)}_{i},a^{(t)}_{-i})\leq
L_{\max}\sqrt{\log N/T}+\kappa\cdot
L_{\max}+L_{\max}\sqrt{\log(1/\delta)/(2T)}.$

with probability at least $1-\delta$. Setting
$T=L_{\max}^{2}\varepsilon^{-2}\log(N/\delta)$ and $\kappa=\min(1/2,)$ then
gives the average regret is bounded by $O(\varepsilon)$. ∎

###### Proof of Corollary 2.11.

By Lemma 2.10, if we set
$T=C\cdot\left(L_{\max}^{2}\varepsilon^{-2}\log(Nm/\eta)\right)$,
$\delta=\varepsilon/(CL_{\max})$ for a sufficient large constant and simulate
the repeated game playing for $T$ rounds where each player makes her decision
based on $\delta$-RWM, the regret of the $i$-th player is bounded by
$\varepsilon$ with probability at least $1-\delta/m$. By union bound, this
holds for all players simultaneously. Ths results then follows from Theorem
2.8. ∎

##  Appendix B Bit-complexity and Stability of Numeric Operations

To avoid un-necessary technical details of the bit-complexities of numbers and
time complexities of algebraic operations, the algorithmic results in the main
body of this work are stated in the _Algebraic Computation Model_. In
particular, we assume additions, subtractions, multiplications, divisions,
exponentiation and comparisons can be carried between real numbers in constant
time, and that the computing device has query access to the digits of the real
numbers. We remark that, for all the games we study, if the rewards are
rational numbers with bounded bit complexities, our algorithms can all be
implemented exactly in the RAM computation model with their runtime increased
by at most polynomial factors. Unsurprisingly, if one is more careful with the
numeric precision needed and maintains only multiplicative approximations of
each algebraic operations, our algorithms can be implemented in the RAM model
losing only poly-logarithmic factors. In this section, we discuss some
standard techniques to this end for the reader interested in any
implementation of our algorithms.

We first discuss how we represent and perform algebraic operations on numbers
whose absolute values are exponentially large or small. While writing these
numbers down exactly is costly, for the purpose of $\delta$-approximate
sampling, it is actually sufficient t keep $\operatorname{poly}(\delta)$
multiplicative approximations of these numbers. Fortunately, these
approximations can indeed be represented much more succinctly using the
scientific notations. For convenience, for a number $a\in\mathbb{R}^{+}$, we
will call $\tilde{a}$ a $\delta$-approximation of $a$ if we have
$(1-\delta)a\leq\tilde{a}\leq(1+\delta)a$ and a one-sided
$\delta$-approximation if we have $(1-\delta)a\leq\tilde{a}\leq a$.

######  Fact B.1.

Given $a\in\mathbb{R}^{+}$ satisfying $\exp(-q)\leq a\leq\exp(q)$ for
$q\in\mathbb{Z}^{+}$, let $\tilde{a}$ be $a$ written in scientific form
keeping $\Theta(\log(1/\delta))$ many significant figures. Then, $\tilde{a}$
is a one-sided $\delta$-approximation of $a$ and can be represented using
$\Theta(\log(\delta^{-1})+\log q)$ many bits.

Instead of performing exact arithmetic computations, we can perform
‘approximate’ arithmetic operations on real numbers in all our algorithms.

######  Claim B.2.

Let $a,b\in\mathbb{R}^{+}$ be two numbers in scientific notations with $s$
significant figures.

  * •

one-sided $\delta$-approximation of additions and multiplications can be
performed in time $O(s)$ and $O(s\log s)$ respectively.

  * •

$\delta$-approximation of division can be performed in time
$O(s+\log(1/\delta))$.

  * •

Given $1\leq\alpha\leq 2$ that has $s$ significant figures and
$i\in\mathbb{Z}^{+}\cup\\{0\\}$, a one-side $\delta$-approximation of
$\alpha^{i}$ can be computed in time
$O(s\cdot\log^{2}i\cdot\log(\delta^{-1}))$.

###### Proof.

The first two claims follow from the definition of (one-sided)
$\delta$-approximation. We proceed to show that one can perform approximate
exponentiation efficiently. In particular, we argue that $\alpha^{i}$ can be
computed fairly accurately via fast exponentiation while keeping
$C\cdot\log(i/\delta)$ significant figures throughout the computation for some
large enough constant $C$. By doing so, we can make sure the approximation to
$\alpha$ is a one-sided $\xi$-approximation where $\xi=\delta/i^{c}$ for some
large enough constant $c$. Consequently, the approximation of $\beta^{j}$ for
any $j\in[i]$ that is a power of $2$ is within $(1\pm
4^{\log_{2}(i)}\cdot\xi)$. It then follows $\alpha^{i}$ can be approximated
within $(1\pm\xi\cdot O(\log i)\cdot 4^{\log_{2}(i)})=(1\pm\delta)$ when $c$
is sufficiently large. ∎

Unsurprisingly, the output of applying a series of arithmetic operations will
be within multiplicative factors of the result obtained by replacing each
operation with its approximate counterpart.

######  Fact B.3.

Given a variable $y$ that is the result of $V$ arithmetic operations including
Addition, Multiplication and Division on the inputs
$x_{1},\cdots,x_{n}\in\mathbb{R}^{+}$ in scientific notations with $s$
significant figures, let $\hat{y}$ be the variable obtained by replacing all
the arithmetic operations with their $\delta/(10V)$-approximate counterparts
for small enough $\delta$. Then, $\hat{y}$ will be a $\delta$-approximation of
$y$. Moreover, if only additions and multiplications are used, the
approximation is one-sided i.e. $\hat{y}\leq y$.

Careful readers may find that subtraction is excluded when we discuss
approximate algebraic operations in Claim B.2. For two numbers
$a,b\in\mathbb{R}^{+}$ and $\hat{a},\hat{b}$ be their $\delta$ approximations
counterparts, $\hat{a}-\hat{b}$ may be wildly different from $a-b$ when $a$ is
substantially larger than $b$. Yet, subtraction between real numbers is indeed
used in two different places. Firstly, subtractions occur in Discrete Fourier
Transform, which is used in Proposition 4.5 to compute the convolution between
functions. The numeric stability of DFT varies among different implementations
and depends on a number of subtle factors (See [66]). If numeric stability
indeed becomes an issue in the actual implementation, one can fallback to
evaluate the convolution in the brute-force manner, which increases the
complexity from $O(nk\log n)$ to $O(n^{2}k)$.

Another place where subtractions are used is in Algorithm 2 to compute range
sum of piece-wise function efficiently. As such, we need a numerically more
stable technique for performing range sum query in place of the prefix sum
technique. In particular, given the succinct description of a $q$-piecewise
constant function $f:\\{0\\}\cup[n]\mapsto\mathbb{R}^{+}$ and $\delta\in(0,1)$
beforehand, we want to perform some preprocessing in time
$q\cdot\operatorname{polylog}(n,q,\delta^{-1})$ and then answer a series of
queries of the form $\sum_{i=a}^{b}f(i)$ within $(1\pm\delta)$ multiplicative
factors in time $\operatorname{polylog}(n,q,\delta^{-1})$.

######  Claim B.4.

Given the succinct description $D_{f}$ of a $q$-piecewise constant function
$f:\\{0\\}\cup[n]\mapsto\mathbb{R}$ where the function values contain at most
$s$ significant figures, there exists an algorithm Range-Sum-Query which
performs some preprocess in time
$O(qs)+\operatorname{polylog}(q,\delta^{-1},n)$, and can compute one-sided
$\delta$-approximation to query of the form $\sum_{i=a}^{b}f(i)$ in time
$\operatorname{polylog}(n,q,\delta^{-1})$.

###### Proof.

Let $D_{f}=\\{(a_{1},b_{1},y_{1}),\cdots(a_{q},b_{q},y_{q})\\}$. In the
preprocessing step, we first compute the range sum of all intervals
$[a_{i},b_{i}]$ (approximately). Denote the results as an array
$[s_{1},\cdots,s_{q}]$. Then, we build a segment tree with the array, where
nodes store the approximate range sum of intervals of lengths that are powers
of $2$. This takes $O(q\log(q))$ arithmetic operations. The data structure
then allows us to answer $\sum_{i=a}^{b}s_{i}$ with $O(\log(q))$ arithmetic
operations with enough accuracy. Then, when we receives a query
$\sum_{i=a}^{b}f(i)$. We first binary search for the intervals of $f$ that
$a,b$ fall into respectively. Next, one uses the pre-built segment tree to
answer the range sum of any intervals that are strictly contained in $[a,b]$
and then adds the sum of remaining elements. It is then not hard to see that
the pre-processing step takes time at most
$O(qs)+\operatorname{polylog}(q,\delta^{-1},n)$, and answering each query
takes time $\operatorname{polylog}(q,\delta^{-1},n)$. ∎

Lastly, we discuss the building block of sampling: sampling from multinomial
distributions. Typically, our algorithm computes a vector $w_{1},\cdots,w_{n}$
and then samples from the multinomial distribution $X$ where $\Pr[X=i]\propto
w_{i}$. This is simple in the Algebraic computation model as one can easily
reduce this to sampling from uniform distributions over real intervals (which
can be done at assumed unit cost). In particular, one first computes the
prefix sum $W_{1},\cdots W_{n}$. Then, one samples $z$ from the uniform
distribution over the interval $[0,W_{n}]$ and returns $j$ for
$W_{j-1}<z<W_{j}$. This clearly takes at most time $O(n)$. In the bit-
complexity model, we can nonetheless achieve approximate sampling from
arbitrary multinomial distributions with similar runtime.

######  Claim B.5.

Given a weight vector $(w_{1},\cdots,w_{n})$ in scientific notations with $s$
significant figures, $\delta$-approximate sampling from the multinomial
distribution $X$ such that $\Pr[X=i]\propto w_{i}$ can be done in time
$\widetilde{O}(n\cdot(\log(\delta^{-1})+s))$.

###### Proof.

After reading the input, one first truncates to make sure each $w_{i}$ has at
most $O(\log(n/\delta))$ many significant figures as that is already enough
for the specified ‘sampling accuracy’. After that, all arithmetic operations
will be carried out with their $c\cdot\delta/n$ approximations for some
sufficiently small constant $c$. In the next step, one normalizes the weight
vector and rounds each $w_{i}$ to their nearest multiple of $c\cdot\delta/n$.
Doing so changes the distribution by at most $c\cdot\delta$ in total variation
distance. One can then multiply all $w_{i}$ by a factor of $n/(c\cdot\delta)$
to make everything an integer. Finally, one can do the same thing as sampling
in the Algebraic computation model: computing the prefix sums and reducing the
problem to sampling from uniform distributions, now over integer intervals.
The integers in the interval can be at most $n^{2}/(c\cdot\delta)$ so the
runtime is dominated by the preliminary computations performed. ∎

We note that in many cases we actually require a slightly more complicated
sampling procedure where we wish to sample from a $q$-piecewise support-$n$
multinomial with $q\ll n$. This can be done similarly in time
$\widetilde{O}((q+\log(n))(\log(\delta^{-1})+s))$ by first sampling one of the
$q$ piecewise intervals by the above technique, then sampling uniformly within
the interval.

##  Appendix C Implementing Glauber Dynamics

This section is devoted to proving Theorem 3.8, which we repeat here for
convenience.

######  Theorem C.1 (RWM on Matroids).

Let $\mathcal{I}=\\{(A_{i},R_{i})_{i=1}^{m}\\}$ be an $m$-player game on a
size-$n$ ground set $\Omega$. If $A_{i}$ consists of the bases of a rank-$k$
matroid, is linear, and is collision-sensitive with support $q$, then it is
possible to implement $\delta$-RWM${}^{T}_{\beta}$ in time

$O\left(k_{i}T(CO+q\log(n)+mk_{-i}T\log(n))\log\left(\frac{k_{i}\log(n)+L_{\max}T\log(\beta^{-1})}{\delta}\right)\right),$

assuming access to a $q$-piecewise succinct description of
$\operatorname{NC}_{i}$ encoded under an ordering of $\Omega$ and a
contraction oracle matching the same ordering.

###### Proof.

Since $A_{i}$ is given by the bases of a rank-$k$ matroid, the single-step
Glauber Dynamics on $\Omega^{w}$ mix in time

$T(GD(\Omega^{w}),\delta)\leq
O\left(k\log\left(\frac{\log(|A_{i}|/w_{*})}{\delta}\right)\right)$

for any external field $w\in\mathbb{R}_{+}^{\Omega}$. In our setting, we have
that $|A_{i}|\leq n^{k}$, and $w_{*}\geq\beta^{L_{\max}T}$, and the process
needs to be repeated once per round bringing the complexity to:

$O\left(kT\log\left(\frac{k\log(n)+L_{\max}T\log(\beta^{-1})}{\delta}\right)\right)$

times the implementation time of a single step of Glauber dynamics.

It is left to bound this cost. To implement a step of GD in the $t$-th round
of optimistic hedge, we first remove a uniformly random element from our
current basis, then re-sample from the conditional distribution. The first
step can easily be implemented in $O(\log(k))$ time. The latter step requires
more care. Let $\hat{e}$ denote the $(k-1)$-size set resulting from the down-
step of the walk. Query the contraction oracle on $\hat{e}$ and call the
resulting set $S_{e}\subset E$. Notice that by definition, the conditional
measure of any $x\in S_{e}$ is proportional to $\beta^{-r_{t}(x)}$ where

$r_{t}(x)=\sum\limits_{j=1}^{t}R_{i}^{\Omega}(x,s^{(j)}).$

Thus to perform the conditional sampling efficiently, it is sufficient to
compute the external field for each element in $S_{e}$ and sample from the
corresponding multinomial distribution.

While implementing this naively would require time at least $|S_{e}|$ to check
the weight of each element in the conditional distribution, this can be
circumvented via our assumption that our game is collision-sensitive with
bounded support. In particular, assume for the moment we have access to a
succinct description for the vertex-wise total rewards $r_{t}(v)$ that is
$(q+tmk_{-i})$-piecewise, and that the output of the contraction oracle
respects the order of the description (we will argue this can be constructed
efficiently shortly). As a result, the total rewards in $S_{e}$ are
$(q+tmk_{-i})$-piecewise as well. This means that using query access to
CO,131313Formally we are also assuming here one has query access to the size
of the output of the contraction oracle. Note this can be easily implemented
in $\text{polylog}(n)$ time even if one does not assume such access. we can
build a succinct description for total rewards on the elements in $|S_{e}|$
(labeled by their index in CO). Sampling from the corresponding multinomial
distribution in the algebraic computation model then takes $O((q+tmk_{-i}))$
time, and one can then feed the sampled index into CO to receive the correct
vertex. Altogether, a single step of GD can therefore be implemented in
$O(CO+(q+tmk_{-i})\log(n))$ time assuming access to the appropriate
description of total rewards.

Finally, we argue we can construct and maintain the succinct descriptions of
the vertex-wise reward functions over $T$ rounds efficiently. Recall we start
with an $q$-piecewise succinct description for the no-collision vertex-wise
reward values. In each round, at most $mk_{-i}$ new elements of $\Omega$ are
introduced into the history, and since the game is collision-sensitive the
resulting succinct description of rewards is at most $(q+tmk_{-i})$-piecewise
in the $t$-th round as desired. The computational cost stems from noting that
it is actually sufficient just to update the rewards for vertices which have
appeared in the opponent history (and the number of rounds in which it has
appeared). During look-up, computing the total reward for any vertex $v$ that
has appeared $t$ times can be computed in $O(1)$ time by simply adding the
stored value $(T-t)R_{i}(v,s)$ for any $s\notni v$. The cost of building the
succinct description is therefore asymptotically dominated by the sampling
procedure above, which gives the final complexity. ∎